# Adaptive optimal estimation of irregular mean and covariance functions

Steven Golovkine111Renault Technocentre & Univ. Rennes, Ensai, CNRS, CREST -
UMR 9194, F-35000 Rennes, France<EMAIL_ADDRESS>Nicolas
Klutchnikoff222Univ. Rennes, CNRS, IRMAR - UMR 6625, F-35000 Rennes, France;
<EMAIL_ADDRESS>Valentin Patilea333Univ. Rennes, Ensai,
CNRS, CREST - UMR 9194, F-35000 Rennes, France<EMAIL_ADDRESS>
###### Abstract

We propose straightforward nonparametric estimators for the mean and the
covariance functions of functional data. Our setup covers a wide range of
practical situations. The random trajectories are, not necessarily
differentiable, have unknown regularity, and are measured with error at
discrete design points. The measurement error could be heteroscedastic. The
design points could be either randomly drawn or common for all curves. The
definition of our nonparametric estimators depends on the local regularity of
the stochastic process generating the functional data. We first propose a
simple estimator of this local regularity which takes strength from the
replication and regularization features of functional data. Next, we use the
“smoothing first, then estimate” approach for the mean and the covariance
functions. The new nonparametric estimators achieve optimal rates of
convergence. They can be applied with both sparsely or densely sampled curves,
are easy to calculate and to update, and perform well in simulations.
Simulations built upon a real data example on household power consumption
illustrate the effectiveness of the new approach.

Key words: Functional data analysis; Hölder exponent; Local polynomials;
Minimax optimality

MSC2020:  62R10; 62G05; 62M09

##  1 Introduction

Motivated by a large number of applications, there is a great interest in
models for observation entities in the form of a sequence of measurements
recorded intermittently at several discrete points in time. Functional data
analysis (FDA) considers such data as being values on the trajectories of a
stochastic process, recorded with some error, at discrete random times. The
mean and the covariances functions play a critical role in FDA.

To formalize the framework, let $\mathcal{T}\subset\mathbb{R}$ be a non-
degenerate bounded, open interval, typically $(0,1)$. Data consist of random
realizations of sample paths from a second-order stochastic process
$X=(X_{t}:t\in\mathcal{T})$ with continuous trajectories. The mean and
covariance functions are $\mu(t)=\mathbb{E}(X_{t})$ and

$\Gamma(s,t)=\mathbb{E}\left\\{[X_{s}-\mu(s)][X_{t}-\mu(t)]\right\\}=\mathbb{E}\left(X_{s}X_{t}\right)-\mu(s)\mu(t),\qquad
s,t\in\mathcal{T},$

respectively. If the independent realizations
$X^{(1)},\dotsc,X^{(i)},\dotsc,{X^{(N)}}$ of $X$ were observed, the ideal
estimators would be

$\widetilde{\mu}_{N}(t)=\frac{1}{N}\sum_{i=1}^{N}X^{(i)}_{t},\qquad
t\in\mathcal{T},$

and

$\widetilde{\Gamma}_{N}(s,t)=\frac{1}{N-1}\sum_{i=1}^{N}\\{X^{(i)}_{s}-\widetilde{\mu}_{N}(s)\\}\\{X^{(i)}_{t}-\widetilde{\mu}_{N}(t)\\},\qquad
s,t\in\mathcal{T}.$ (1)

In real applications, the curves are rarely observed without error and never
at each value $t\in\mathcal{T}$. This is why we consider the following common
and more realistic setup. For each $1\leq i\leq N$, and given a positive
integer $M_{i}$, let ${T^{(i)}_{m}}\in\mathcal{T}$, $1\leq m\leq M_{i}$, be
the observation times for the curve $X^{(i)}$. The observations associated
with a curve, or trajectory, $X^{(i)}$ consist of the pairs
$({Y^{(i)}_{m}},{T^{(i)}_{m}})\in\mathbb{R}\times\mathcal{T}$ where
${Y^{(i)}_{m}}$ is defined as

${Y^{(i)}_{m}}={X^{(i)}}({T^{(i)}_{m}})+{\varepsilon^{(i)}_{m}},\qquad 1\leq
m\leq M_{i},\;\;1\leq i\leq N,$ (2)

and ${\varepsilon^{(i)}_{m}}$ is an independent (centered) error variable.
Here and in the following we use the notation $X^{(i)}_{t}$ for the value at a
generic point $t\in\mathcal{T}$ of the realization $X^{(i)}$ of $X$, while
$X^{(i)}({T^{(i)}_{m}})$ denotes the measurement at ${T^{(i)}_{m}}$ of this
realization.

A commonly used idea is to build feasible versions of
$\widetilde{\mu}_{N}(\cdot)$ and $\widetilde{\Gamma}_{N}(\cdot,\cdot)$ using
nonparametric estimates of $X^{(i)}_{t}$ and $X^{(i)}_{s}X^{(i)}_{t}$, such as
obtained by smoothing splines or local polynomials. This approach, usually
called ‘smoothing first, then estimate’ or ‘two-stage procedure’, have been
considered, amongst others, by Hall et al.,  (2006) and Zhang and Chen,
(2007). In general, the sample trajectories are required to admit at least
second-order derivatives over the whole domain $\mathcal{T}$. Li and Hsing,
(2010), Zhang and Wang,  (2016) and Zhang and Wang,  (2018) propose an
alternative local linear smoothing approach where the estimators are
determined by suitably weighting schemes which involve the whole sample of
curves. This idea exploits the so-called replication and regularization
features of functional data (see Ramsay and Silverman, , 2005, ch. 22). In
this alternative approach, the regularity assumptions are imposed on the mean
and covariance functions, which are required to admit second, or higher, order
derivatives over the whole domain. Since, in general, the mean and covariance
functions are more regular than the sample trajectories, the approach based on
weighting schemes using all the sample curves might be preferable. However, in
some cases, for instance in energy, chemistry and physics, astronomy and
medical applications, the mean and covariance functions could be quite
irregular, of unknown irregularity.

Cai and Yuan,  (2011) and Cai and Yuan,  (2010) derive the optimal rates of
convergence, in the minimax sense, for the mean and covariance functions,
respectively, and propose optimal estimators. The optimal estimator of the
mean function proposed by Cai and Yuan,  (2011) is a smoothing spline
estimator which could be built only if the regularity of the sample paths is
known. For the covariance function, Cai and Yuan,  (2010) used the
representation of the covariance function in a tensor product reproducing
kernel Hilbert space (RKHS) space. Next, they derived estimators for
$\Gamma(s,t)$ using a low dimension version of this representation obtained by
a regularization procedure, provided the values $M_{i}$ are not very
different. This procedure does not require a given regularity, but involve
some numerical optimization. See also Wong and Zhang,  (2019). The minimax
optimal rates for the mean and covariance functions are defined by the sum of
two types of terms. One corresponds to the rate of convergence of the
$\widetilde{\mu}_{N}(\cdot)$ and $\widetilde{\Gamma}_{N}(\cdot,\cdot)$, which
is the standard rate of convergence for empirical means and covariances. The
other contribution to the optimal rates is given by the differences between
$\widetilde{\mu}_{N}(\cdot)$ and $\widetilde{\Gamma}_{N}(\cdot,\cdot)$ and
their feasible versions. The optimal rates of the differences depend on the
sample trajectories regularity. The reason is that the minimax lower bounds
should also take into account, the case where the trajectories have the same
regularity as the function to be estimated.

The estimation of the mean and covariance functions presents another specific
feature. The minimax optimal rates of convergence depend on the nature of the
measurement times ${T^{(i)}_{m}}$. For now, two situations were investigated
in the literature. On the one hand, the so-called _independent design_ case
where, given the $M_{i}$’s, the ${T^{(i)}_{m}}$ are obtained as a random
sample of size $M_{1}+\cdots+M_{N}$ from the same continuous distribution. On
the other hand, the so-called _common design_ case where the $M_{i}$ are all
equal to some integer value $\mathfrak{m}$, and the ${T^{(i)}_{m}}$, $1\leq
m\leq\mathfrak{m}$, are the same across the curves $X^{(i)}$. In both cases,
the best rates for the nonparametric estimators depend on the regularity of
the sample trajectories. These rates also depend on the number of different
observation times ${T^{(i)}_{m}}$, that is equal to $M_{1}+\cdots+M_{N}$ with
independent design, and equal to $\mathfrak{m}$ with common design. In other
words, the replication feature of functional data is less impactful with
common design. See Cai and Yuan,  (2011) for the case of the mean function,
and Cai and Yuan,  (2010) and Cai and Yuan,  (2016) for the covariance
function case.

In this paper, we propose simple “smoothing first, then estimate” type
methods, based on $1-$dimensional smoothing, which achieve minimax optimal
rates of convergence. The process is allowed to have a piecewise constant,
unknown regularity. Our method does not require complex numerical
optimization. It applies in the same way with common and independent design
situations, and allows for general heteroscedastic measurement errors
${\varepsilon^{(i)}_{m}}$. Moreover, our approach is suitable with both
sparsely or densely sampled curves. The definition of sparse and dense regimes
is recalled in Section 2.

Let ${\widehat{X}^{(i)}}$ be a suitable nonparametric estimator of the curve
$X^{(i)}$ applied to the $M_{i}$ pairs $({Y^{(i)}_{m}},{T^{(i)}_{m}})$, for
instance a local polynomial estimator. What will make this estimator suitable
is that it takes into account, the regularity of the process $X$ and the final
estimation purpose, that is the mean or the covariance function. These
features can be achieved in an easy, data-driven way, as will be explained
below. For now, with at hand the ${\widehat{X}^{(i)}}$’s tuned for the mean
function estimation, let us define

$\widehat{\mu}_{N}(t)=\frac{1}{N}\sum_{i=1}^{N}{\widehat{X}^{(i)}}_{t},\qquad
t\in\mathcal{T}.$ (3)

For the covariance function, we distinguish the diagonal from the non-diagonal
points. For now, with at hand the ${\widehat{X}^{(i)}}$’s tuned for the
covariance function estimation, and for some $\mathfrak{d}\geq 0$ determined
using the data, let us define

$\widehat{\Gamma}_{N}(s,t)=\frac{1}{N}\sum_{i=1}^{N}{\widehat{X}^{(i)}}_{s}{\widehat{X}^{(i)}}_{t}-\widehat{\mu}_{N}(s)\widehat{\mu}_{N}(t),\qquad
s,t\in\mathcal{T},\;\;|s-t|>\mathfrak{d}.$ (4)

Moreover, for $0\leq t-s\leq\mathfrak{d}$,

$\widehat{\Gamma}_{N}(s,t)=\widehat{\Gamma}_{N}(u-\mathfrak{d}/2,u+\mathfrak{d}/2),\quad\text{with
}\;u=(s+t)/2,$ (5)

and $\widehat{\Gamma}_{N}(t,s)=\widehat{\Gamma}_{N}(s,t)$. It is well known
that the variance function $\Gamma(s,s)$ induces a singularity when estimating
the covariance function $\Gamma(\cdot,\cdot)$. See, for instance, Zhang and
Wang,  (2016), Remark 4. This singularity causes little problem when studying
pointwise rates of convergence and confidence intervals. In that cases, one
could set $\mathfrak{d}=0$ and use a modified covariance function estimator
only for the diagonal $\widehat{\Gamma}_{N}(s,s)$. However, the diagonal
singularity could deteriorate the rate of convergence of
$\widehat{\Gamma}_{N}$ in the integrated squared norm. In general, this rate
deterioration is propagated and affects the rates of convergence for the
estimators of the eigenvalues and eigenfunctions of the covariance operator.
The idea underlying (5) is natural and the diagonal tuning parameter
$\mathfrak{d}$ decreases to zero according to a data-driven rule that we
provide in the following.

Although the methodology we propose is general and can be used with different
types of smoothers, we focus on the case where the ${\widehat{X}^{(i)}}_{t}$
are obtained by local polynomials with a compactly supported kernel. In this
case tuning the ${\widehat{X}^{(i)}}$’s means to suitably determine the rate
of decrease and the constant defining the bandwidth. In our case, this is done
completely data-driven by a one variable minimization of a new, suitable risk
function.

To the best of our knowledge, there is no contribution following the
‘smoothing first, then estimate’ idea, which considers estimators of the
curves $X^{(i)}$ adapted to their regularity and to the purpose of estimating
mean or covariance functions. It is clear that trajectory-by-trajectory
adaptive optimal smoothing, for instance using the Goldenshluger and Lepski,
(2011) method, in general yields sub-optimal rates of convergence for
$\widehat{\mu}_{N}(t)$ and $\widehat{\Gamma}_{N}(s,t)$. The reason is that
trajectory-by-trajectory smoothing ignores the information contained in the
other $N-1$ curves in the sample generated according to the same stochastic
process $X$. See Cai and Yuan,  (2011) for a discussion on the differences
with the usual nonparametric rates. One can also use cross-validation for
choosing the bandwidth with the suitably weighting schemes, such as proposed
by Li and Hsing,  (2010) or Zhang and Wang,  (2016). However, this would
require significant computational effort, and, to the best of our knowledge,
the idea has not yet received a theoretical justification. Using the
replication and regularization features of functional data, we propose a new
effective estimator for the local regularity of the process $X$, a
probabilistic concept which determines the analytic regularity of the
trajectories of $X$. The replication feature of the functional data makes the
concept of local regularity of the process a more meaningful parameter than
the usual curve regularity, which is an analytic concept designed for a single
function. Our local regularity estimator, inspired by, but different from, the
estimator introduced by Golovkine et al.,  (2020), combines information both
across and within curves. Moreover, it allows for general heteroscedastic
measurement errors, does not involve any optimization and is obtained after a
fast, possibly parallel, computation. With at hand the local regularity
estimator, we derive the suitable estimators ${\widehat{X}^{(i)}}_{t}$, and
finally our optimal mean and covariance functions estimators. The smoothing
parameter used to build the ${\widehat{X}^{(i)}}_{t}$ depends on $M_{i}$ and
$N$, but can be easily computed given the estimate of the local regularity of
$X$. It is worth noting that our estimators have good theoretical and finite
sample properties without using a bivariate nonparametric estimator of
$X^{(i)}_{s}X^{(i)}_{t}$. We explain this somehow unexpected result by the
particular separable structure of the bivariate functions $(s,t)\mapsto
X^{(i)}_{s}X^{(i)}_{t}$ which can be estimated at the faster rate for
univariate functions.

In Section 2 we provide insight on why the local regularity of the process $X$
is a natural feature to be considered. Moreover, we explain why the ‘smoothing
first, then estimate’ approach could achieve optimal rates when the regularity
of $X$ is known. In Section 3, we formally define the local regularity of the
process $X$ and discuss the relationship between this regularity and the
decrease rate for the eigenvalues of the covariance operator. Moreover, we
introduce a new estimator for this regularity and present exponential bounds
for the concentration under mild conditions. In particular, both independent
and common designs are allowed, and the process regularity is allowed to be
piecewise constant. Section 3 ends with a discussion on the relationship
between the process regularity and the trajectories’ analytical regularity. In
Section 4, we use the regularity estimate to build sharp bounds of the
pointwise quadratic risk function between our estimators and the unfeasible
estimators $\widetilde{\mu}_{N}$ and $\widetilde{\Gamma}_{N}$, respectively.
The bounds depend on quantities which could be estimated by sample averages.
Minimizing the risk bounds with respect to the bandwidth, we derive the
optimal bandwidth for the local polynomial estimates of the trajectories.
These estimates are further used to optimally estimate the mean and covariance
functions. The finite sample performance of the new estimators is illustrated
in Section 5 using simulated samples generated according to the setup of a
real data set on the power consumption of households. Some conclusions and
discussions are gathered in Section 6. Few proofs are relegated to the
Appendix. A Supplementary Material contains more proofs, technical arguments
and simulation results.

##  2 From unfeasible to feasible optimal estimators

Let us first provide insight into the reason why the local regularity of the
process generating the curves is a meaningful concept, and why our approach
can achieve good performance. For this purpose, we analyze the difference
$\widehat{\mu}_{N}(t)-\widetilde{\mu}_{N}(t)$, $s,t\in\mathcal{T},$ but
similar apply to the covariance function estimation.

The data $({Y^{(i)}_{m}},{T^{(i)}_{m}})\in\mathbb{R}\times\mathcal{T}$ are
generated according to model (2) with

${\varepsilon^{(i)}_{m}}=\sigma({T^{(i)}_{m}},{X^{(i)}}({T^{(i)}_{m}})){e^{(i)}_{m}},\quad
1\leq m\leq M_{i},\;1\leq i\leq N,$ (6)

where the ${X^{(i)}}$ are independent trajectories of $X$, ${e^{(i)}_{m}}$ are
independent copies of a centered variable $e$ with unit variance, and
$\sigma(t,x)$ is some unknown bounded function which account for possibly
heteroscedastic measurement errors. The integers $M_{1},\dotsc,M_{N}$
represent an independent sample of an integer-valued random variable $M$ with
expectation $\mathfrak{m}$ which increases with $N$. Thus,
$M_{1},\dotsc,M_{N}$ is the $N$th line of a triangular array of integer
numbers. In the independent design case, for each $1\leq i\leq N$, the
observation times ${T^{(i)}_{m}}$ are random realizations of a variable
$T\in\mathcal{T}$. We assume that the realizations of $X$, $M$ and $T$ are
mutually independent. Let $\mathcal{T}_{obs}^{(i)}$ denote the set of
observation times ${T^{(i)}_{m}}$, $1\leq i\leq M_{i}$, over the trajectory
$X^{(i)}$. In the common design case, $M\equiv\mathfrak{m}$, and the
$\mathcal{T}_{obs}^{(i)}$ are the same for all $i$. Thus, if not stated
differently, the issues discussed in this section apply to both independent
design and common design cases.

Let

$\mathbb{E}_{i}(\cdot)=\mathbb{E}(\cdot\mid
M_{i},\mathcal{T}_{obs}^{(i)},X^{(i)})\quad\text{ and
}\quad\mathbb{E}_{M,T}(\cdot)=\mathbb{E}\left(\cdot\mid
M_{i},\mathcal{T}_{obs}^{(i)},1\leq i\leq N\right).$

For any $t\in\mathcal{T}$, we consider a generic linear nonparametric
estimator

${\widehat{X}^{(i)}}_{t}=\sum_{m=1}^{M_{i}}{Y^{(i)}_{m}}W_{m}^{(i)}(t),\quad
1\leq i\leq N.$ (7)

The weights $W_{m}^{(i)}(t)$ are defined as functions of the elements in
$\mathcal{T}_{obs}^{(i)}$. The example we keep in mind is that of local
polynomials which we investigate in detail in Section 4. Let

${\widehat{X}^{(i)}}_{t}-X^{(i)}_{t}=B_{t}^{(i)}+V_{t}^{(i)},\quad
t\in\mathcal{T},$ (8)

where

$B_{t}^{(i)}:=\mathbb{E}_{i}\left[{\widehat{X}^{(i)}}_{t}\right]-X^{(i)}_{t}\quad\text{
and }\quad
V_{t}^{(i)}:={\widehat{X}^{(i)}}_{t}-\mathbb{E}_{i}\left[{\widehat{X}^{(i)}}_{t}\right]=\sum_{m=1}^{M_{i}}{\varepsilon^{(i)}_{m}}W_{m}^{(i)}(t).$

The pairs of random variables $(B_{t}^{(i)},V_{t}^{(i)})$, $1\leq i\leq N$,
are independent and we could reasonably assume that they are squared
integrable for all $t$.

For illustration purposes, let us suppose that ${\widehat{X}^{(i)}}_{t}$ in
(7) are well defined for all $1\leq i\leq N$. The general case is presented in
Section 4. Then, for the mean, we can write

$\widehat{\mu}_{N}(t)-\widetilde{\mu}_{N}(t)=\frac{1}{N}\sum_{i=1}^{N}B_{t}^{(i)}+\frac{1}{N}\sum_{i=1}^{N}V_{t}^{(i)}.$
(9)

All the variables ${\varepsilon^{(i)}_{m}}$ are centered and conditionally
independent, with bounded conditional variance, given all $M_{i}$,
$\mathcal{T}_{obs}^{(i)}$ and $X^{(i)}$. Thus, given $M_{1},\ldots,M_{N}$ and
$\mathcal{T}_{obs}^{(1)},\ldots,\mathcal{T}_{obs}^{(N)}$,

$\mathbb{E}_{M,T}\left[\left\\{N^{-1}\sum_{i=1}^{N}V_{t}^{(i)}\right\\}^{2}\right]=N^{-1}\mathbb{E}_{M,T}\left[N^{-1}\sum_{i=1}^{N}\left\\{V_{t}^{(i)}\right\\}^{2}\right]\\\
\leq N^{-1}\sup_{x}\sigma^{2}(t,x)\times
N^{-1}\sum_{i=1}^{N}\left\\{\max_{m}\left|W_{m}^{(i)}(t)\right|\times\sum_{m=1}^{M_{i}}\left|W_{m}^{(i)}(t)\right|\right\\}.$
(10)

For local polynomials with bandwidth $h$, under some mild conditions, the rate
of decrease of the right-hand side in the last display, given the design, is
$O_{\mathbb{P}}\left((N\mathfrak{m}h)^{-1}\right)$.

On the bias part, let us suppose for the moment that the trajectories are not
differentiable. By Cauchy-Schwarz inequality, we have

$\mathbb{E}_{M,T}\left[\left\\{N^{-1}\sum_{i=1}^{N}B_{t}^{(i)}\right\\}^{2}\right]\leq
N^{-1}\sum_{i=1}^{N}\mathbb{E}_{M,T}\left[\left\\{B_{t}^{(i)}\right\\}^{2}\right]\\\
\leq
N^{-1}\sum_{i=1}^{N}\left\\{\sum_{m=1}^{M_{i}}\left|W_{m}^{(i)}(t)\right|\times\sum_{m=1}^{M_{i}}\mathbb{E}_{M,T}\left(\left\\{{X^{(i)}}({T^{(i)}_{m}})-X^{(i)}_{t}\right\\}^{2}\mid\mathcal{T}_{obs}^{(i)}\right)|W_{m}^{(i)}(t)|\right\\}.$
(11)

It now becomes clear that the rate of the square of the bias term in the
difference $\widehat{\mu}_{N}(t)-\widetilde{\mu}_{N}(t)$ is determined by the
second-order moment of the increments $X^{(i)}_{s}-X^{(i)}_{t}$. If, for $s$
and $t$ which are close,

$\mathbb{E}\left(\left\\{X^{(i)}_{s}-X^{(i)}_{t}\right\\}^{2}\right)\approx
L_{0}^{2}|t-s|^{2H_{0}},$ (12)

with some constants $0<H_{0}\leq 1$ and $L_{0}$, then the rate of the right-
hand side quantity in (11) can be bounded by

$N^{-1}\sum_{i=1}^{N}\left\\{\sum_{m=1}^{M_{i}}\left|W_{m}^{(i)}(t)\right|\times\sum_{m=1}^{M_{i}}L_{0}^{2}\left|{T^{(i)}_{m}}-t\right|^{2H_{0}}|W_{m}^{(i)}(t)|\right\\}.$
(13)

For local polynomials with bandwidth $h$, this leads to
$O_{\mathbb{P}}\left(h^{2H_{0}}\right).$ We call $H_{0}$ the Hölder exponent
of $X^{(i)}$, whenever the trajectories are not differentiable.

With $\delta-$times differentiable trajectories, $\delta\geq 1$ integer, the
condition (12) should be considered with the difference of the values of the
$\delta-$th derivative of $X^{(i)}$ instead of the $X^{(i)}$ itself. Then,
following the lines of Proposition 1.13 of Tsybakov,  (2009), we use the
Taylor expansion and the property
$\sum_{m=1}^{M_{i}}({T^{(i)}_{m}}-t)^{d}W_{m}^{(i)}(t)=0$ which is satisfied
by any integers $0\leq d\leq\delta$ and $1\leq i\leq N$, in the case of local
polynomials of order $\delta$. The rate in (11) then becomes
$O_{\mathbb{P}}\left(h^{2(\delta+H_{\delta})}\right)$, where
$H_{\delta}\in(0,1]$ denotes the Hölder exponent of the $\delta-$th
derivatives of $X^{(i)}$.

Gathering facts, we deduce that, in the case of $\delta-$times differentiable
trajectories, $\delta\geq 0$, for local polynomial of order $\delta$ and $h$
of order $(N\mathfrak{m})^{-1/(1+2\\{\delta+H_{\delta}\\})}$, as long as this
bandwidth order does not yield degenerate estimates, we obtain

$\mathbb{E}_{M,T}\left[\\{\widehat{\mu}_{N}(t)-\widetilde{\mu}_{N}(t)\\}^{2}\right]=O_{\mathbb{P}}\left((N\mathfrak{m})^{-\;\frac{2(\delta+H_{\delta})}{1+2(\delta+H_{\delta})}}\right).$
(14)

Thus, given the local regularity $\delta+H_{\delta}$, the estimator
$\widehat{\mu}_{N}(t)$ can achieve the minimax optimal rate for the estimation
of the mean function $\mu(t)$. See Cai and Yuan,  (2011).

Let us note that in some cases, in particular with local polynomials, the
estimator defined in (7) could be degenerate, _i.e._ , the weights
$W_{m}^{(i)}(t)$ are not well defined because $h$ is too small. The
trajectories for which this happens could change with $t$. Then,
$\widehat{\mu}_{N}(t)$ is defined as an average over the trajectories for
which the estimator (7) is not degenerate. This can more likely happen in the
so-called _sparse_ regime, where $\mathfrak{m}^{2(\delta+H_{\delta})}\ll N$. A
similar phenomenon occurs with estimators determined by suitably weighting
schemes, see for instance equation (2.1) in Li and Hsing,  (2010), or equation
(2.3) in Zhang and Wang,  (2016). However, in the independent case, one could
benefit from the replication feature of functional data, because only a
fraction of trajectories will yield non degenerate estimators
${\widehat{X}^{(i)}}_{t}$. The size of this fraction plays a central role in
the sparse regime. This crucial aspect is taken into account in Sections 4.1
and 4.2, where we choose the bandwidths while penalizing the number of
trajectories which yield degenerate estimators.

The case of common design requires some special attention. For simplicity, let
us assume the common design points are equidistant and consider the local
polynomials are built with the kernel supported on $[-1,1]$. In this case, the
bandwidth cannot have a rate smaller than $\mathfrak{m}^{-1}$, otherwise all
the weights $W_{m}^{(i)}(t)$ could all be equal to zero. This means that with
a common design, the optimal bandwidth is given by the minimization of
$h^{2(\delta+H_{\delta})}+(N\mathfrak{m}h)^{-1}$ under the constraint that
$\mathfrak{m}h$ stays away from zero. Without loss of generality, we could set
$h=k/\mathfrak{m}$ with $k$ a positive integer and search $k$ which minimizes
$h^{2(\delta+H_{\delta})}+(N\mathfrak{m}h)^{-1}$. Balancing the two terms, one
expects the optimal $k/\mathfrak{m}$ to have the rate
$(N\mathfrak{m})^{-1/\\{1+2(\delta+H_{\delta})\\}}$. If
$\mathfrak{m}^{2(\delta+H_{\delta})}$ is larger than $N$, _i.e._ in the so-
called _dense_ regime, the optimal $k$ is well defined and
$k\sim(\mathfrak{m}^{2(\delta+H_{\delta})}/N)^{1/\\{1+2(\delta+H_{\delta})\\}}$
and, with this optimal choice,
$\mathbb{E}_{M,T}\left[\\{\widehat{\mu}_{N}(t)-\widetilde{\mu}_{N}(t)\\}^{2}\right]=o_{\mathbb{P}}\left(N^{-1}\right).$
If $\mathfrak{m}^{2(\delta+H_{\delta})}\ll N$, then the constraint that $k\geq
1$ becomes binding, and it is no longer possible to balance the squared bias
term and the variance term. The rate of $h^{2(\delta+H_{\delta})}$ dominates
the rate $(N\mathfrak{m}h)^{-1}$. The minimal rate for
$\mathbb{E}_{M,T}\left[\\{\widehat{\mu}_{N}(t)-\widetilde{\mu}_{N}(t)\\}^{2}\right]$
then corresponds to $k=1$, and is
$O_{\mathbb{P}}\left(\mathfrak{m}^{-2(\delta+H_{\delta})}\right).$ We obtain a
similar conclusion with differentiable trajectories. Gathering facts, we
recover the optimal rate for mean estimation with common design, that is
$O_{\mathbb{P}}\left(\mathfrak{m}^{-2(\delta+H_{\delta})}+N^{-1}\right),$ as
found by Cai and Yuan,  (2011). Finally, let us recall the somehow surprising
message from Cai and Yuan,  (2011), p. 2332: the interpolation is rate optimal
when $\mathfrak{m}^{2(\delta+H_{\delta})}\gg N$ in the case of common design;
smoothing does not improve converges rates. Our contribution on this aspect is
a data-driven rule for the practitioner which completes this theoretical fact
about the interpolation. The adaptive bandwidth rule proposed in Section 4
automatically chooses between smoothing and interpolation with finite sample
sizes.

We learn from the above that the ‘smoothing first, then estimate’ approach can
lead to optimal rates of convergence for estimating the mean function with
independent and common design, as derived by Cai and Yuan,  (2011), provided
the regularity of the process $X$ is known. We will show that achieving
optimal rates using the local regularity is also possible for the covariance
function. This probabilistic concept of regularity, summarized by (12), will
be formally introduced in the next section. Next, we introduce a simple
estimator of the local regularity of $X$. It will be shown that, under some
mild conditions, this estimator concentrates around the true local regularity
faster than a suitable negative power of $\log(\mathfrak{m})$. This suffices
to guarantee that our mean and covariance functions estimators achieve the
same rates as when the local regularity is known.

Let us end this section with a discussion of the differences with the
weighting schemes approach, as for instance considered by Li and Hsing,
(2010) and Zhang and Wang,  (2016). If the regularity of $\mu(\cdot)$ is
known, then one could define $B_{t}^{(i)}$ and $V_{t}^{(i)}$ in (8) centering
by the mean function instead of the trajectory $X^{(i)}_{t}$. Then, one could
derive the rate of
$\mathbb{E}\left[\\{\widehat{\mu}_{N}(t)-\mu(t)\\}^{2}\right]$ and find the
bandwidth which minimizes this rate, exactly as done in Li and Hsing,  (2010)
and Zhang and Wang,  (2016) where $\mu(\cdot)$ is assumed to be twice
differentiable. However, the estimation of the regularity of $\mu(\cdot)$
remains an open problem.

##  3 Local regularity estimator

We will allow the sample paths of the real-valued random process
$X=(X_{t}:t\in\mathcal{T})$ to have piecewise constant regularity. That means
that there exists a finite, but unknown, partition of $\mathcal{T}$ in non-
degenerate intervals, such that, in the interior of each interval of the
partition, the regularity of the process is the same. To allow this
generality, we first need to formally introduce the notion of _local
regularity_ of the process $X$. Given this type of regularity, the Kolmogorov
Extension Theorem allows us to determine the analytic regularity of the
trajectories of $X$. The details are provided in Section 3.4.

The most convenient way to obtain the theoretical properties of our ‘smoothing
first, then estimate’ approach for the mean and covariance functions, will be
to consider that the estimator of the local regularity of $X$ is built from a
separate, independent sample of ${N_{0}}$ trajectories. We will call the data
obtained from this separate random sample of trajectories, generated by the
same stochastic process $X$, the _learning sample_. We will show that, in
theory, the size ${N_{0}}$ of the learning sample can be as small as any
positive power of $\mathfrak{m}$, regardless of the type of design,
independent or common. Our simulation experiences reveal that in practice one
can confidently use the same sample for learning the local regularity of $X$
and for estimating the mean and covariance functions.

###  3.1 Local regularity in quadratic mean

Let $\mathcal{O}_{*}$ be an open subinterval of $\mathcal{T}$, with length
$\Delta_{*}=\operatorname{diam}(\mathcal{O}_{*})$. Whenever it exists, let
$\nabla^{d}X_{t}$ denote the $d-$th derivative, $d\geq 1$, of the generic
curve $X_{t}$ at $t\in\mathcal{T}$. By definition $\nabla^{0}X_{t}\equiv
X_{t}$. The following assumptions ensure that, locally on $\mathcal{O}_{*}$,
the derivatives exist up to the integer order $\delta\geq 0$ and that the
$\delta$-th derivative is regular in quadratic mean.

Assumptions. Let $\delta\in\mathbb{N}$ and $0<H_{\delta}\leq 1$.

  1. (H1)

With probability 1, for any $d\in\\{0,\dotsc,\delta\\}$ the $d$-th order
derivative $\nabla^{d}X_{t}$ of $X_{t}$ exists for all $t\in\mathcal{O}_{*}$,
and satisfies:

$0<\underline{a}_{d}=\inf_{u\in\mathcal{O}_{*}}\mathbb{E}\left[(\nabla^{d}X_{u})^{2}\right]\leq\sup_{u\in\mathcal{O}_{*}}\mathbb{E}\left[(\nabla^{d}X_{u})^{2}\right]=\overline{a}_{d}<\infty.$

  2. (H2)

Two positive constants $S_{\delta}$ and $\beta_{\delta}$ exist such that:

$\left|\mathbb{E}\left[(\nabla^{\delta}X_{t}-\nabla^{\delta}X_{s})^{2}\right]-L_{\delta}^{2}|t-s|^{2H_{\delta}}\right|\leq
S_{\delta}^{2}|t-s|^{2H_{\delta}}\Delta_{*}^{2\beta_{\delta}},\qquad
s,t\in\mathcal{O}_{*}.$

  3. (H3)

$\mathfrak{a}>0$ and $\mathfrak{A}>0$ exist such that, for any
$d\in\\{0,\dotsc,\delta\\}$ and any $p\geq 1$:

$\mathbb{E}\left[|\nabla^{d}X_{t}-\nabla^{d}X_{s}|^{2p}\right]\leq\frac{p!}{2}\mathfrak{a}\mathfrak{A}^{p-2},\qquad
s,t\in\mathcal{O}_{*}.$

######  Definition 1.

For any $\delta\in\mathbb{N}$, $0<H_{\delta}\leq 1$ and $L_{\delta}>0$, the
class $\mathcal{X}(\delta+H_{\delta},L_{\delta};\mathcal{O}_{*})$ is the set
of stochastic processes indexed by $t\in\mathcal{O}_{*}$ for which the
conditions (H1) to (H3) hold true. The quantity $\alpha=\delta+H_{\delta}$ is
the local regularity of the process on $\mathcal{O}_{*}$, while $L_{\delta}$
is the Hölder constant of the $\delta-$th derivative of the trajectories.

See, _e.g._ , Blanke and Vial,  (2014) for some examples and references on
processes satisfying the mild condition in (H2), which we impose on the
$\delta-$th derivative of the trajectories. Examples include stationary or
stationary increment processes. For some common processes with the ordered
eigenvalues of the covariance operator such that, for some $\nu>1$,
$\lambda_{j}\sim j^{-\nu}$, $j\geq 1$, one has $\delta+H_{\delta}=(\nu-1)/2$.
Details are provided in the Supplementary Material. The condition in (H3)
serves to derive the exponential bound for the concentration of the local
regularity estimator. By definition, when the local regularity $\alpha$ is an
integer, $\delta=\alpha-1$ and $H_{\delta}=1$. In Lemma LABEL:lem:regularity-1
in the Supplement we describe the embedding structure of the spaces
$\mathcal{X}(\cdot,\cdot;\mathcal{O}_{*})$. More precisely, we show that
whenever $\delta\geq 1$, for any $1\leq d\leq\delta-1$, the process $X$
restricted to $\mathcal{O}_{*}$ belongs to
$\mathcal{X}(d+1,L_{d};\mathcal{O}_{*})$ for some $L_{d}>0$.

###  3.2 The local regularity estimation method

Let us assume that $X$ restricted to $\mathcal{O}_{*}$ belongs to
$\mathcal{X}(\delta+H_{\delta},L_{\delta};\mathcal{O}_{*})$ for some
$\delta\in\mathbb{N}$, $0<H_{\delta}\leq 1$ and $L_{\delta}>0$. Our first goal
is to construct an estimator of $\alpha=\delta+H_{\delta}$. To do so, we first
have to estimate $H_{d}$ for any $d=0,\dotsc,\delta$. For simplicity, let us
denote

$\theta_{d}(s,t)=\mathbb{E}\left[(\nabla^{d}X_{t}-\nabla^{d}X_{s})^{2}\right]\quad
s,t\in\mathcal{O}_{*}.$

Since the restriction of the process $X$ belongs to
$\mathcal{X}(d+H_{d},L_{d};\mathcal{O}_{*})$, with
$H_{d}=\mathbf{1}_{\\{d\neq\delta\\}}+H_{\delta}\mathbf{1}_{\\{d=\delta\\}}$,
we have

$\theta_{d}(s,t)\approx L_{d}^{2}|t-s|^{2H_{d}}\quad\text{if }\Delta_{*}\text{
is small}.$ (15)

Now, let $t_{1}$ and $t_{3}$ be such that
$[t_{1},t_{3}]\subset\mathcal{O}_{*}$ and $t_{3}-t_{1}=\Delta_{*}/2$. Denote
by $t_{2}$ the middle point of $[t_{1},t_{3}]$. It is easily seen that

$H_{d}\approx\tilde{H}_{d}=\frac{\log(\theta_{d}(t_{1},t_{3}))-\log(\theta_{d}(t_{1},t_{2}))}{2\log(2)}\quad\text{if
}\Delta_{*}\text{ is small}.$ (16)

Consider a learning sample of ${N_{0}}$ curves, generated according to (2),
which yield

$({Y^{(n)}_{m}},{T^{(n)}_{m}})\in\mathbb{R}\times\mathcal{T},\quad 1\leq m\leq
M_{n},\;1\leq n\leq{N_{0}}.$

Then, given a nonparametric estimator $\widetilde{\nabla^{d}X}_{t}$ of
${\nabla^{d}X}_{t}$, for $t\in\mathcal{O}_{*}$, we define a natural estimator
of $\tilde{H}_{d}$, and thus of $H_{d}$, as

$\hat{H}_{d}=\frac{\log\big{(}\hat{\theta}_{d}(t_{1},t_{3})\big{)}-\log\big{(}\hat{\theta}_{d}(t_{1},t_{2})\big{)}}{2\log(2)},$
(17)

where

$\hat{\theta}_{d}(s,t)=\frac{1}{{N_{0}}}\sum_{n=1}^{{N_{0}}}\left(\widetilde{\nabla^{d}X}^{(n)}_{t}-\widetilde{\nabla^{d}X}^{(n)}_{s}\right)^{2}.$
(18)

By the definition of the local regularity,
$\delta=\min\\{d\in\mathbb{N}:H_{d}<1\\}$, which suggests defining

$\hat{\delta}=\min\\{d\in\mathbb{N}:\hat{H}_{d}<1-\varphi(\mathfrak{m})\\},$
(19)

for some decreasing function $\varphi(\cdot)$ which is defined later.
Typically it could decrease as fast as a negative power of
$\log(\mathfrak{m})$. The estimator of the local regularity is then

$\hat{\alpha}=\hat{\delta}+\hat{H}_{\hat{\delta}}.$ (20)

Since (17) can be easily modified when $N_{0}$ increases, $\hat{\alpha}$ can
be updated with little effort. Our local regularity estimator is related to
the estimator introduced by Golovkine et al.,  (2020). However, here we
propose smoothing the trajectories to calculate $\hat{H}_{d}$ for $d=0$, while
Golovkine et al.,  (2020) directly uses the noisy measurements of the
trajectories. To derive the properties of the local regularity estimator, for
simplicity, we suppose that $\mathfrak{m}$ is given. In applications, it could
be estimated by the average of the realizations of $M$.

###  3.3 Concentration properties of the local regularity estimator

The quality of the estimators $\hat{\delta}$ and $\hat{H}_{d}$ depends on the
quality of the generic nonparametric estimators $\widetilde{\nabla^{d}X}$ of
$\nabla^{d}X$. To quantify their behavior, we consider the local
$\mathbb{L}^{p}$-risk

$R_{p}(d)=R_{p}(d;\mathcal{O}_{*})=\sup_{t\in\mathcal{O}_{*}}\mathbb{E}\left(|\xi_{d}(t)|^{p}\right),\quad\text{where
}\quad\xi_{d}(t)=\widetilde{\nabla^{d}X}_{t}-{\nabla^{d}X}_{t}.$

Our method applies with any type of nonparametric estimator
$\widetilde{\nabla^{d}X}$ (local polynomials, splines,…) as soon as, for any
$p\in\mathbb{N}$, its $\mathbb{L}^{p}$-risk is suitably bounded. The following
mild condition is satisfied by common estimators, see for instance Theorem 1
in Gaïffas,  (2007) for the case of local polynomials.

Assumptions.

  1. (LP1)

There exist two positive constants $\mathfrak{c}$ and $\mathfrak{C}$ such that

$R_{2p}(d)\leq\frac{p!}{2}\mathfrak{c}\mathfrak{C}^{p-2},\qquad\forall p\geq
1,\;d\in\\{0,\dotsc,\delta\\}.$

We can now derive an exponential bound for the concentration of all the
estimators $\hat{H}_{d}$, $d\in\\{0,\dotsc,\delta\\}$. To make this
exponential bound useful for deriving optimal rates for our estimators of the
mean and covariance functions, we will require the largest risk among
$R_{2}(0),\ldots,R_{2}(\delta)$ to tend to zero as $\mathfrak{m}$ increases to
infinity.

######  Theorem 1.

Assume that $X$ restricted to $\mathcal{O}_{*}$ belongs to
$\mathcal{X}(\delta+H_{\delta},L_{\delta};\mathcal{O}_{*})$, for some integer
$\delta\geq 0$ and $0<H_{\delta}\leq 1$, and that (LP1) holds. Assume also
that there exists $\tau>0$ and $B>0$ such that:

$\rho_{*}=\max_{d\in\\{0,\dotsc,\delta\\}}R_{2}(d)\leq B\mathfrak{m}^{-\tau}.$

Let $0<\gamma<1$ and $\Gamma>0$, and consider

$\Delta_{*}=2\exp\left(-\log^{\gamma}(\mathfrak{m})\right)\quad\text{and}\quad\varphi(\mathfrak{m})=\log^{-\Gamma}(\mathfrak{m}).$

Then, for any $\mathfrak{m}$ larger than some constant $\mathfrak{m}_{0}$
depending on $B$, $\tau$, $\gamma$, $\Gamma$, $H_{\delta}$, $\beta_{\delta}$
and for some constant $\mathfrak{f}$, we have

$\mathbb{P}\left(|\hat{\alpha}-\alpha|>\varphi(\mathfrak{m})\right)\leq
8(1+\delta)\exp\left(-\mathfrak{f}{N_{0}}\varphi^{2}(\mathfrak{m})\Delta_{*}^{4H_{\delta}}\right).$

The proof of Theorem 1 follows the lines of the results of Golovkine et al.,
(2020) and is thus omitted. The three quantities $\rho_{*}$, $\Delta_{*}$ and
$\varphi(\mathfrak{m})$ are required to decrease to zero, as $\mathfrak{m}$
tends to infinity, in such a way that
$\rho_{*}/\Delta_{*}+\Delta_{*}/\varphi(\mathfrak{m})\rightarrow 0$. We
propose $\Gamma=2$ and $\gamma=1/2$. The choices of the rates for $\rho_{*}$,
$\Delta_{*}$ and $\varphi(\mathfrak{m})$ satisfy some additional requirements.
First, it will be shown below that, in order to achieve optimal rates of
convergence for the mean and covariance estimators, the local regularity has
to be estimated with a concentration rate $\varphi(\mathfrak{m})$ faster than
$\log^{-1}(\mathfrak{m})$. This is a consequence of the identity
$\mathfrak{m}^{1/\log(\mathfrak{m})}=e$ for any $\mathfrak{m}>1$ and of a mild
condition on the rates of increase of $N$ and $\mathfrak{m}$; see (32) below.
Second, we want to allow for reasonable rates of increase for $N_{0}$, the
size of the learning set. In Theorem 1, ${N_{0}}$ can increase as fast as an
arbitrary positive power of $\mathfrak{m}$. Third, since $\tau>0$ could be
arbitrarily small, the rate imposed on the nonparametric estimators
$\widetilde{\nabla^{d}X}$ of $\nabla^{d}X$ is a very mild consistency
requirement which is achieved by the common estimators, with random or fixed
design, under mild conditions, in particular on the distribution of the
$M_{i}$. See, for instance, Tsybakov,  (2009) and Belloni et al.,  (2015). In
particular, the required rate for the $\widetilde{\nabla^{d}X}$ can be
obtained under general forms of heteroscedasticity.

###  3.4 From the regularity of the process to the regularity of the
trajectories

Let us now connect the probabilistic concept of local regularity with the
regularity of the sample paths considered as functions. Let $\mathcal{O}_{*}$
be some open subinterval of $\mathcal{T}$. We say that a function is $\beta$
times differentiable on $\mathcal{O}_{*}$ if the function has an up to
$\lfloor\beta\rfloor-$order derivative and the $\lfloor\beta\rfloor-$th
derivative is Hölder continuous with exponent $\beta-\lfloor\beta\rfloor$.
(For a real number $a$, let $\lfloor a\rfloor$ denotes the largest integer not
exceeding $a$.) Let us call a Hölder space of local regularity $\beta$, the
set of functions which are $\beta$ times differentiable on $\mathcal{O}_{*}$.

By Assumption (H1), for almost all realizations of the process $X$, the
derivatives of the sample path exist up to order $\delta$. Under a suitable
moment condition on the increments the derivative process $\nabla^{\delta}X$,
and using the refined version of Kolmogorov’s criterion stated in Revuz and
Yor,  (2013), it can be proven that, for any
$\delta<\beta<\delta+H_{\delta}=\alpha$, the sample path of the process $X$
restricted to $\mathcal{O}_{*}$ belong to the Hölder space of exponent $\beta$
over $\mathcal{O}_{*}$. Technical details are provided in the Supplementary
Material. As an example, the Brownian motion has a local regularity equal to
1/2 on any open interval. Almost surely, the sample paths of the Brownian
motion belong to any Hölder space of local regularity $\beta<1/2$, but cannot
have a Hölder continuity of order $\alpha$, $\alpha\geq 1/2$. Hence, the
probability theory indicates that imposing assumptions on the regularity of
the sample paths could be a delicate issue. Indeed, even for some widely used
examples, this regularity is not well defined in the sense required by the
nonparametric statistics theory. Since the sample paths have a regularity
which can be arbitrarily close to the (local) regularity of the process $X$ as
defined above, the probabilistic concept of local regularity seems more
appropriate for establishing the rates of convergence for the mean and
covariance estimators.

##  4 Optimal mean and covariance function estimators

We now explain how to suitably select the bandwidths for the local polynomial
smoothing of the trajectories, and next build mean and covariance function
estimates. Our data-driven adaptive bandwidth rules lead to optimal rate
estimates whenever the estimator of the local regularity concentrates to the
true value faster than $\log^{-1}(\mathfrak{m})$. Theorem 1 then guarantees
that the suitable rate of the bandwidth can be achieved with high probability.

Motivated by the applications, we allow the process $X$ to have a piecewise
constant regularity.

Assumptions.

  1. (E1)

There exist $K\geq 1$ and $\min\mathcal{T}=t_{0}<t_{1}<\ldots
t_{K-1}<t_{K}=\max\mathcal{T}$ such that, for each
${\mathcal{T}_{k}^{\circ}}=(t_{k-1},t_{k})$,

$X\text{ restricted to }\mathcal{T}_{k}^{\circ}\text{ belongs to
}\mathcal{X}(\delta_{k}+H_{\delta_{k}},L_{\delta_{k}};{\mathcal{T}_{k}^{\circ}}),$

for some integer $\delta_{k}\geq 0$, $0<H_{\delta_{k}}\leq 1$ and
$L_{\delta_{k}}>0$.

For any $t\in\mathcal{T}$, let $\alpha_{t}$ be the local regularity from
Definition 1 corresponding to a small neighborhood of $t$, and let
$L_{\delta}$, with $\delta=\lfloor\alpha_{t}\rfloor$, denote the corresponding
Hölder constant. By our assumptions, $\alpha_{t}$ is well defined everywhere
except a finite number of points in $\mathcal{T}$. For a small neighborhood of
points like $t_{k}$, the values of $\alpha_{t}$ are different for $t$ to the
left and to the right of $t_{k}$. However, for practical purposes, any
convention would work, such as for instance, setting $\alpha_{t_{k}}$ equal to
the average of the neighboring values. Hereafter therefore we assume that
$\alpha_{t}$ is well-defined everywhere on $\mathcal{T}$.

Hereafter, $\hat{\alpha}_{t}$ will be the estimator of $\alpha_{t}$ defined in
(20) and (17), built with an independent learning sample, and with
$\varphi(\mathfrak{m})$ replaced by $\log^{-2}(\widehat{\mathfrak{m}})$, where
$\widehat{\mathfrak{m}}=N^{-1}\sum_{i=1}^{N}M_{i}$. If $t_{1}<\ldots<t_{K-1}$
are known, then for each $t\in\mathcal{T}_{k}^{\circ}$, $\hat{\alpha}_{t}$ can
be defined as in (20) and (17) with $t_{2}$ the middle point in
${\mathcal{T}_{k}^{\circ}}$. In practice the $t_{1},\ldots,t_{K-1}$ are likely
not given. One could then estimate the local regularity on a grid on points in
$\mathcal{T}$, and let each point on the grid play the role of $t_{2}$ in
(17). Such implementation is used in Section 5.

For each $1\leq i\leq N$, we use a local polynomials (LP) approach to build
suitable nonparametric estimators ${\widehat{X}^{(i)}}$. For
${t}\in\mathcal{T}$, if the local regularity $\alpha_{t}$ is available, using
the measurements $({Y^{(i)}_{m}},{T^{(i)}_{m}})$, $1\leq m\leq M_{i}$, of a
generic trajectory $X^{(i)}$, we consider the $LP(\lfloor\alpha_{t}\rfloor)$
estimator defined by

${\widehat{X}^{(i)}}_{t}={\widehat{X}^{(i)}}_{t}(h_{t})=\sum_{m=1}^{M_{i}}{Y^{(i)}_{m}}W_{m}^{(i)}(t),\qquad
1\leq i\leq N,$ (21)

with a suitable bandwidth $h_{t}$ which depends on $\alpha_{t}$, and

$W_{m}^{(i)}(t)=W_{m}^{(i)}(t;h)=\frac{1}{M_{i}h}U^{\top}(0)A_{M_{i}}^{(i)}(t,h)^{-1}U\left(\frac{{T^{(i)}_{m}}-{t}}{h}\right)\\!K\\!\left(\frac{{T^{(i)}_{m}}-{t}}{h}\right),\quad
1\leq m\leq M_{i},$ (22)

with

$A_{M_{i}}^{(i)}(t,h)=\frac{1}{M_{i}h}\sum_{m=1}^{M_{i}}U\left(\frac{{T^{(i)}_{m}}-{t}}{h}\right)U^{\top}\left(\frac{{T^{(i)}_{m}}-{t}}{h}\right)K\left(\frac{{T^{(i)}_{m}}-{t}}{h}\right).$

Here, for any $z\in\mathbb{R}$,
$U(z)=U(z;\alpha_{t})=(1,z,\dotsc,z^{\lfloor\alpha_{t}\rfloor}/\lfloor\alpha_{t}\rfloor!)^{\top}$
and $K:\mathbb{R}\to\mathbb{R}_{+}$ is a continuous kernel with the support in
$[-1,1]$.

For a given ${t}\in\mathcal{T}$ and each $1\leq i\leq N$, the weights
$W_{m}^{(i)}(t)$, $1\leq m\leq M_{i}$ are well defined as soon as the matrix
$A_{M_{i}}^{(i)}(t,h)$ is invertible. When $A_{M_{i}}^{(i)}(t,h)$ is not
invertible, we consider the smoothing of the curve $i$ at point $t$ as
degenerate and the curve $i$ should not be considered in the construction of
the estimator of $\mu(t)$. A similar reasoning applies to the covariance
estimator. In the case of common design, for each $t$, the number of
degenerate estimates ${\widehat{X}^{(i)}}_{t}$ is either equal to $N$ or to
zero. In the independent design case, this number could be any integer between
0 and $N$. A suitable bandwidth rule should be penalizing for the number of
curves which are not considered for the estimation. In the following sections
we propose a natural way to penalize which adapts to the sparse and dense
regimes. Moreover, the two types of designs are automatically handled.

###  4.1 Adaptive mean estimation

Let $k_{0}$ be some integer, and let $\mathbf{1}\\{\cdot\\}$ denote the
indicator function. For any ${t}\in\mathcal{T}$, define

$w_{i}(t;h)=1\quad\text{if
}\;\sum_{m=1}^{M_{i}}\mathbf{1}\\{|{T^{(i)}_{m}}-t|\leq h\\}\geq
k_{0},\quad\text{ and }\;w_{i}(t;h)=0\text{ otherwise},$ (23)

and let

$\mathcal{W}_{N}(t;h)=\sum_{i=1}^{N}w_{i}(t;h).$

Our adaptive mean function estimator is

$\widehat{\mu}_{N}^{*}(t)=\widehat{\mu}_{N}(t;h_{\mu}^{*})\quad\text{ with
}\quad\widehat{\mu}_{N}(t;h)=\frac{1}{\mathcal{W}_{N}(t;h)}\sum_{i=1}^{N}w_{i}(t;h){\widehat{X}^{(i)}}_{t}.$
(24)

Here, ${\widehat{X}^{(i)}}_{t}$ is the local polynomial estimator
$LP(\lfloor\hat{\alpha}_{t}\rfloor)$ defined in (21) with some suitable
bandwidth $h_{\mu}^{*}$ which is defined below. The mean estimator
$\widehat{\mu}_{N}(t;h)$ is a practical version of that defined in (3) which
takes into account that some trajectories could have less than $k_{0}$
observation times between $t-h_{\mu}^{*}$ and $t+h_{\mu}^{*}$. The threshold
defined by $k_{0}$ avoids considering degenerate ${\widehat{X}^{(i)}}_{t}$.
For this purpose, it has to be greater than or equal to
$\lfloor\hat{\alpha}_{t}\rfloor+1$; more details are provided in the
Supplementary Material. The normalization of the mean estimator by
$\mathcal{W}_{N}(t;h)$ is also implicitly used in the definition of the
estimators proposed by Li and Hsing,  (2010) and Zhang and Wang,  (2016).

To introduce our bandwidth rule, for any $h>0$, $\alpha>0$, let

$c_{i}(t;h)=\sum_{m=1}^{M_{i}}\left|W_{m}^{(i)}(t;h)\right|\quad\text{and}\quad
c_{i}(t;h,\alpha)\\!=\sum_{m=1}^{M_{i}}\left|({T^{(i)}_{m}}-t)/h\right|^{\alpha}\left|W_{m}^{(i)}(t;h)\right|,$
(25)

and

$\overline{C}_{1}(t;h,\alpha)=\\!\frac{1}{\mathcal{W}_{N}(t;h)}\sum_{i=1}^{N}w_{i}(t;h)c_{i}(t;h)c_{i}(t;h,\alpha).$
(26)

With the Nadaraya-Watson (NW) estimator, all the $c_{i}(t;h)$ are equal to 1.
Moreover,

$\overline{C}_{1}(t;h,\alpha)\approx\int|u|^{\alpha}K(u)du.$ (27)

Using the equivalent kernels idea, see section 3.2.2 in Fan and Gijbels,
(1996), the same approximation could be used in the case of local linear
estimators. The accuracy of the approximation (27) could be high since it
involves the ${T^{(i)}_{m}}$ close to $t$ for all the curves with
$w_{i}(t;h)=1$. Next, using the rule $0/0=0$, let

$\mathcal{N}_{i}(t;h)=\frac{w_{i}(t;h)}{\max_{1\leq m\leq
M_{i}}|W_{m}^{(i)}(t;h)|}\quad\text{and}\quad\mathcal{N}_{\mu}(t;h)=\left[\frac{1}{\mathcal{W}_{N}^{2}(t;h)}\sum_{i=1}^{N}w_{i}(t;h)\frac{c_{i}(t;h)}{\mathcal{N}_{i}(t;h)}\right]^{-1}.$
(28)

With the NW estimator, $\mathcal{N}_{\mu}(t;h)$ is equal to
$\mathcal{W}_{N}(t;h)$ times the harmonic mean of $\mathcal{N}_{i}(t;h)$, over
the curves with $w_{i}(t;h)=1$. Moreover, under the mild condition (33) below,
$\mathcal{N}_{i}(t;h)=O_{\mathbb{P}}(\mathfrak{m}h)$.

Quantities like $c_{i}(t;h)$, $c_{i}(t;h,\alpha)$ and
$\max_{m}|W_{m}^{(i)}(t,h)|$ are commonly used to bound the risk of LP
estimators, see also Section 2 above. For deriving theoretical results,
$c_{i}(t;h)$ and $c_{i}(t;h,\alpha)$ are usually bounded by a constant.
Meanwhile, the maximum of the absolute values of the LP weights could be
bounded by a suitable constant divided by the number of observation times
between $t-h_{\mu}^{*}$ and $t+h_{\mu}^{*}$. See Guillou and Klutchnikoff,
(2011) or Tsybakov,  (2009). In the context of functional data, we directly
use the $c_{i}(t;h)$, $c_{i}(t;h,\alpha)$ and $\mathcal{N}_{i}(t;h)$ and thus
entirely exploit the information contained in all the curves. The
computational complexity remains low and the reward is a sharper risk bound, a
better bandwidth choice and an improved mean estimator in finite samples.

We define the bandwidth for computing $\widehat{\mu}_{N}^{*}(t)$ such that it
minimizes the mean squared difference between $\widehat{\mu}_{N}(t;h)$ and
$\widetilde{\mu}_{N}(t)$. This leads us to define the optimal bandwidth

$h_{\mu}^{*}=h_{\mu}^{*}(t)=\arg\min_{h>0}\mathcal{R}_{\mu}(t;h),$ (29)

with

$\mathcal{R}_{\mu}(t;h)=\\!q_{1}^{2}h^{2\hat{\alpha}_{t}}+\frac{q_{2}^{2}}{\mathcal{N}_{\mu}(t;h)}+q_{3}^{2}\left[\frac{1}{\mathcal{W}_{N}(t;h)}-\frac{1}{N}\right],$
(30)

and

$q_{1}^{2}=\frac{\overline{C}_{1}(t;h,2\hat{\alpha}_{t})}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\hat{L}_{\delta}^{2}\;,\quad
q_{2}^{2}=\sigma^{2}_{\max}\;,\quad q_{3}^{2}=Var(X_{t}),$ (31)

where $\sigma_{\max}$ is a bound for the function $\sigma(t,x)$ in (6) and
$\hat{L}_{\delta}$ is an estimate of the Hölder constant $L_{\delta}$ from
Assumption (E1). In Section 5, we propose a simple procedure to build
$\hat{L}_{\delta}$ based on the preliminary nonparametric estimates of the
sample paths used for $\hat{\alpha}_{t}$. We show in the Appendix that
$2\mathcal{R}_{\mu}(t;h)$ is a sharp bound for
$\mathbb{E}_{M,T}\left[\left\\{\widehat{\mu}_{N}(t;h)-\widetilde{\mu}_{N}(t)\right\\}^{2}\right]$.
The maximization of $\mathcal{R}_{\mu}(t;h)$ can be easily performed on a grid
of $h$ values. With some additional effort, $\mathcal{R}_{\mu}(t;h)$ can also
be minimized with respect to $k_{0}$ in (23) over a small set of integers
greater than $\lfloor\hat{\alpha}_{t}\rfloor$.

The bandwidth rule (29) could be used with both independent and common design.
With common design, the ${T^{(i)}_{m}}\equiv T_{m}$ and
$W_{m}^{(i)}(t;h)\equiv W_{m}(t;h)$ no longer depend on $i$ and the solution
$h_{\mu}^{*}$ will always be a value in the set of $h$ such that
$\mathcal{W}_{N}(t;h)=N$. Moreover, whenever $\mathcal{W}_{N}(t;h)=N$,

$\overline{C}_{1}(t;h,2\hat{\alpha}_{t})\\!=\sum_{m=1}^{\mathfrak{m}}\\!\left|W_{m}(t;h)\right|\sum_{m=1}^{\mathfrak{m}}\\!\left|\frac{T_{m}\\!-\\!t}{h}\right|^{2\hat{\alpha}_{t}}\\!\\!\left|W_{m}(t;h)\right|\quad\text{and}\quad\mathcal{N}_{\mu}(t;h)\\!=\frac{N\sum\limits_{m=1}^{\mathfrak{m}}\\!|W_{m}(t;h)|}{\max\limits_{1\leq
m\leq\mathfrak{m}}\\!|W_{m}(t;h)|}.$

In a data-driven way, $h_{\mu}^{*}$ automatically chooses between
interpolation and smoothing.

The following result states that our estimator $\widehat{\mu}_{N}^{*}(t)$
defined by (24) and (29) achieves the optimal rate. For simplicity, we assume
that

$\limsup_{N,\mathfrak{m}\rightarrow\infty}\\{\log(N)/\log(\mathfrak{m})\\}<\infty,$
(32)

a technical condition which matches general situations found in applications.
We also impose the following mild technical condition in the independent
design case:

$\exists c_{L},C_{U}>0\;\text{ such that }\;c_{L}\leq
M_{i}\mathfrak{m}^{-1}\leq C_{U},\;\;\text{ for all }N\text{ and }1\leq i\leq
N.$ (33)

Moreover, in the common design case, where $M_{i}\equiv\mathfrak{m}$ and the
$T_{1}^{(i)},\ldots,T_{\mathfrak{m}}^{(i)}$ are not changing with $i$, we
suppose that:

$\exists C_{U}\geq 1\;\text{ such that }\;\;\max_{1\leq
m\leq\mathfrak{m}-1}\\{T_{m+1}^{(i)}-T_{m}^{(i)}\\}\leq C_{U}\min_{1\leq
m\leq\mathfrak{m}-1}\\{T_{m+1}^{(i)}-T_{m}^{(i)}\\}.$ (34)

######  Theorem 2.

Assume that ${T^{(i)}_{m}}$ are either independently drawn, have a continuous
density which is bounded away from zero and (33) holds true, or
${T^{(i)}_{m}}$ represent the points of a common design satisfying (34).
Moreover, (32) and Assumption (E1) hold true. For $t\in\mathcal{T}$, let
$\hat{\alpha}_{t}$ be an estimator of $\alpha_{t}<1$ computed on a separate
sample such that
$\hat{\alpha}_{t}-\alpha_{t}=o_{\mathbb{P}}(\log^{-1}(\mathfrak{m}))$. Then,
the estimator $\widehat{\mu}_{N}^{*}(t)=\widehat{\mu}_{N}(t;h_{\mu}^{*})$
defined by (24) and (29) satisfies

$\widehat{\mu}_{N}^{*}(t)-\widetilde{\mu}_{N}(t)=O_{\mathbb{P}}\left((N\mathfrak{m})^{-\;\frac{\alpha_{t}}{1+2\alpha_{t}}}\right)\;\;\;\;\text{and}\;\;\;\;\widehat{\mu}_{N}^{*}(t)-\mu(t)=O_{\mathbb{P}}\left((N\mathfrak{m})^{-\;\frac{\alpha_{t}}{1+2\alpha_{t}}}+N^{-1/2}\right),$

in the independent design case. Meanwhile, with the common design,

$\widehat{\mu}_{N}^{*}(t)-\widetilde{\mu}_{N}(t)=O_{\mathbb{P}}\left(\max\left\\{(N\mathfrak{m})^{-\;\frac{\alpha_{t}}{1+2\alpha_{t}}},\mathfrak{m}^{-\;\alpha_{t}}\right\\}\right)=O_{\mathbb{P}}\left(\mathfrak{m}^{-\;\alpha_{t}}\right),$

and

$\widehat{\mu}_{N}^{*}(t)-\mu(t)=O_{\mathbb{P}}\left(\max\left\\{(N\mathfrak{m})^{-\;\frac{\alpha_{t}}{1+2\alpha_{t}}},\mathfrak{m}^{-\;\alpha_{t}}\right\\}+N^{-1/2}\right)=O_{\mathbb{P}}\left(\mathfrak{m}^{-\;\alpha_{t}}\right).$

The rates achieved by $\widehat{\mu}_{N}^{*}(t)$ are the best one could expect
in view of the results of Cai and Yuan,  (2011). To avoid additional technical
arguments, we only prove our result for a local regularity less than 1 and
Nadaraya-Watson estimators ${\widehat{X}^{(i)}}_{t}$. We conjecture that it
also true for $\alpha_{t}\geq 1$, but we leave formal justification for future
work. The difference between the common and independent designs comes from the
fact that, in order to avoid degenerate mean estimator, the bandwidth cannot
decrease faster than $\mathfrak{m}^{-1}$.

###  4.2 Adaptive covariance function estimates

For any $s,{t}\in\mathcal{T}$, $s\neq t$, define

$w_{i}(s,t;h)=w_{i}(s;h)w_{i}(t;h)\qquad\text{and}\qquad\mathcal{W}_{N}(s,t;h)=\sum_{i=1}^{N}w_{i}(s,t;h),$
(35)

with $w_{i}(s;h)$ and $w_{i}(t;h)$ as in (23). Our adaptive covariance
function estimator is

$\widehat{\Gamma}^{*}_{N}(s,t)=\widehat{\Gamma}_{N}(s,t;h_{\Gamma}^{*})\quad\text{with}\quad\widehat{\Gamma}_{N}(s,t;h)=\widehat{\gamma}_{N}(s,t;h)-\widehat{\mu}_{N}^{*}(s)\widehat{\mu}_{N}^{*}(t),$
(36)

where $\widehat{\mu}_{N}^{*}(s)$, $\widehat{\mu}_{N}^{*}(t)$ are defined
according to (24) with the corresponding bandwidths, and

$\widehat{\gamma}_{N}(s,t;h)=\frac{1}{\mathcal{W}_{N}(s,t;h)}\sum_{i=1}^{N}w_{i}(s,t;h){\widehat{X}^{(i)}}_{s}{\widehat{X}^{(i)}}_{t}.$
(37)

Here, ${\widehat{X}^{(i)}}_{s}$ and ${\widehat{X}^{(i)}}_{t}$ are the local
polynomial estimators $LP(\lfloor\hat{\alpha}_{s}\rfloor)$ and
$LP(\lfloor\hat{\alpha}_{t}\rfloor)$ built, respectively, with some suitable
bandwidth $h_{\Gamma}^{*}$ which is defined below. This covariance function
estimator is a practical version of that defined in (4). The normalization of
the mean estimator by $\mathcal{W}_{N}(s,t;h)$ is also implicitly used in the
definition of the estimators proposed by Li and Hsing,  (2010) and Zhang and
Wang,  (2016).

We define the bandwidth for computing $\widehat{\gamma}_{N}(s,t;h)$, and
eventually $\widehat{\Gamma}_{N}^{*}(s,t)$, such that it minimizes the mean
squared difference between $\widehat{\gamma}_{N}(s,t;h)$ and the unfeasible
estimator
$\widetilde{\gamma}_{N}(s,t)=N^{-1}\sum_{i=1}^{N}X^{(i)}_{s}X^{(i)}_{t}$ of
$\mathbb{E}(X^{(i)}_{s}X^{(i)}_{t})$. To this aim, we define modified versions
of $\mathcal{N}_{i}(t;h)$ and $\mathcal{N}_{\mu}(t;h)$, see (28), taking into
account only the curves with $w_{i}(s,t;h)=1$:

$\mathcal{N}_{i}(t|s;h)=\frac{w_{i}(s,t;h)}{\max_{1\leq m\leq
M_{i}}|W_{m}^{(i)}(t,h)|},$ (38)

and

$\mathcal{N}_{\Gamma}(t|s;h)=\left[\frac{1}{\mathcal{W}_{N}^{2}(s,t;h)}\sum_{i=1}^{N}w_{i}(s,t;h)\frac{c_{i}(t;h)}{\mathcal{N}_{i}(t|s;h)}\right]^{-1},$
(39)

where the $c_{i}(t;h)$ are defined as in (25). This idea leads us to define
the optimal bandwidth as

$h_{\Gamma}^{*}=h_{\Gamma}^{*}(s,t)=\arg\min_{h>0}\\{\mathcal{R}_{\Gamma}(s|t;h)+\mathcal{R}_{\Gamma}(t|s;h)\\},$
(40)

with

$\mathcal{R}_{\Gamma}(t|s;h)=\mathfrak{q}_{1}^{2}h^{2\hat{\alpha}_{t}}+\frac{\mathfrak{q}_{2}^{2}}{\mathcal{N}_{\Gamma}(t|s;h)}+\mathfrak{q}_{3}^{2}\left[\frac{1}{\mathcal{W}_{N}(s,t;h)}-\frac{1}{N}\right].$
(41)

In the last equation, $\mathfrak{q_{\ell}}$, $1\leq l\leq 3$, are defined by:

$\mathfrak{q}_{1}^{2}=2\mathbb{E}(X_{s}^{2})\frac{\overline{\mathfrak{C}}_{1}(t;h,2\hat{\alpha}_{t})}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\hat{L}_{\delta}^{2}\qquad\mathfrak{q}_{2}^{2}=\sigma^{2}_{\max}\;\mathbb{E}(X_{s}^{2}),\qquad\mathfrak{q}_{3}^{2}=\frac{Var(X_{s}X_{t})}{2},$
(42)

where

$\overline{\mathfrak{C}}_{1}(t|s;h,\alpha)=\frac{\sum_{i=1}^{N}w_{i}(s,t;h)c_{i}(t;h)c_{i}(t;h,\alpha)}{\mathcal{W}_{N}(s,t;h)}.$
(43)

With NW or local linear estimators,

$\overline{\mathfrak{C}}_{1}(t|s;h,\alpha)\approx\int|u|^{\alpha}K(u)du.$ (44)

The details for (44) are provided in the Supplementary material. The
definition of $\mathcal{R}_{\Gamma}(s|t;h)$ is the same as in (41) with the
occurrences of $s$ and $t$ switched. We show in the Supplementary Material
that the function of $h$ minimized in (40) is a sharp bound for
$\mathbb{E}_{M,T}\left[\left\\{\widehat{\gamma}_{N}(s,t;h)-\widetilde{\gamma}_{N}(s,t)\right\\}^{2}\right]/2$.
The sum of the first two terms in the expressions of
$\mathcal{R}_{\Gamma}(s|t;h)$ and $\mathcal{R}_{\Gamma}(t|s;h)$ represents the
quadratic risk of our estimator of $\mathbb{E}(X_{s}X_{t})$ compared to the
unfeasible one based on the true values $X^{(i)}_{s}X^{(i)}_{t}$ from the
curves yielding non-degenerate estimates
${\widehat{X}^{(i)}}_{s}{\widehat{X}^{(i)}}_{t}$. Like for the mean function,
the third term in (41) penalizes for the number of curves which are dropped
when calculating our estimator. The minimization of
$\mathcal{R}_{\Gamma}(s|t;h)+\mathcal{R}_{\Gamma}(t|s;h)$ can be performed on
a grid of values $h$. The minimization could be also done over a set of values
$k_{0}$ used in (23).

Like for the mean function, with obvious modifications, the definition (40)
could be used with both independent and common design. Indeed, with common
design, the solution $h_{\Gamma}^{*}$ will always be a value in the set of $h$
such that $\mathcal{W}_{N}(s,t;h)=N$. In a completely data-driven way,
$h_{\Gamma}^{*}$ will automatically choose between interpolation and
smoothing.

######  Theorem 3.

Assume that the conditions of Theorem 2 are met for $s,t\in\mathcal{T}$ such
that $s\neq t$. Moreover,
$\sup_{t\in\mathcal{T}}\mathbb{E}(X_{t}^{4})<\infty$. Let
$\alpha(s,t)=\min\\{\alpha_{s},\alpha_{t}\\}$ and assume $\alpha(s,t)<1$. Then
the estimator
$\widehat{\Gamma}_{N}^{*}(s,t)=\widehat{\Gamma}_{N}^{*}(s,t;h_{\Gamma}^{*})$
defined by (36) and (40) satisfies

$\widehat{\Gamma}_{N}^{*}(s,t)-\widetilde{\Gamma}_{N}(s,t)=O_{\mathbb{P}}\left((N\mathfrak{m})^{-\;\frac{\alpha(s,t)}{1+2\alpha(s,t)}}\right)$

and

$\widehat{\Gamma}_{N}^{*}(s,t)-\Gamma(s,t)=O_{\mathbb{P}}\left((N\mathfrak{m})^{-\;\frac{\alpha(s,t)}{1+2\alpha(s,t)}}+N^{-1/2}\right),$

in the independent design case. Meanwhile with the common design,

$\widehat{\Gamma}_{N}^{*}(s,t)\\!-\widetilde{\Gamma}_{N}(s,t)=O_{\mathbb{P}}\\!\left(\max\left\\{(N\mathfrak{m})^{-\;\frac{\alpha(s,t)}{1+2\alpha(s,t)}},\mathfrak{m}^{-\;\alpha(s,t)}\right\\}\right)=O_{\mathbb{P}}\left(\mathfrak{m}^{-\;\alpha(s,t)}\right),$

and

$\widehat{\Gamma}_{N}^{*}(s,t)-\Gamma(s,t)=O_{\mathbb{P}}\left(\max\left\\{(N\mathfrak{m})^{-\;\frac{\alpha(s,t)}{1+2\alpha(s,t)}},\mathfrak{m}^{-\;\alpha(s,t)}\right\\}+N^{-1/2}\right)=O_{\mathbb{P}}\left(\mathfrak{m}^{-\;\alpha(s,t)}\right).$

The rates achieved by $\widehat{\Gamma}_{N}^{*}(s,t)$ are the best one could
expect in view of the results of Cai and Yuan,  (2010). For now, we only prove
our result for the case $\min\\{\alpha_{s},\alpha_{t}\\}<1$ and Nadaraya-
Watson estimators ${\widehat{X}^{(i)}}_{t}$. We conjecture that it also true
in the general case.

###  4.3 The estimator on the diagonal band of the covariance function

As mentioned in (4) and (5), we propose to use the estimator (37) only outside
the diagonal set $\\{(s,t):|s-t|\leq\mathfrak{d}\\}$, for some suitable
$\mathfrak{d}>0$. It remains to give a data-driven rule for choosing
$\mathfrak{d}$ decreasing to zero, and propose an estimator for
$\mathbb{E}(X_{s}X_{t})$ when $s$ and $t$ are in the diagonal set. Let
${\mathcal{T}_{k}^{\circ}}$, $1\leq k\leq K$, as defined in Assumption (E1).
With a piecewise constant regularity, the value of $\mathfrak{d}$ depends on
$k$.

To understand how to build a covariance estimator which is optimal in
integrated squared norm, for $1\leq k\leq K$ and
$s,t\in{\mathcal{T}_{k}^{\circ}}$ such that $s\leq t$, let $u=(s+t)/2$ and

$\widetilde{D}(\mathfrak{d}):=\iint_{t-\mathfrak{d}\leq s\leq
t}\left\\{\widetilde{\Gamma}_{N}(u-\mathfrak{d}/2,u+\mathfrak{d}/2)-\widetilde{\Gamma}_{N}(s,t)\right\\}^{2}dsdt.$

Let $\alpha$ denote the local regularity corresponding to $u=(s+t)/2$. Under
mild assumptions on the moments of $X_{t}$ and $X_{t}-X_{s}$, we show in the
Appendix that

$\widetilde{D}(\delta)=O_{\mathbb{P}}\left(\mathfrak{d}^{2\alpha+1}\right).$
(45)

Thus, for $(s,t)$ inside the diagonal band, it suffices to use
$\widehat{\Gamma}_{N}(u-\mathfrak{d}/2,u+\mathfrak{d}/2)$ defined according to
(36) when $s\leq t$ and use the symmetry of the covariance function when
$s>t$. That is, in order to estimate the covariance function at a point from
the diagonal set, we simply apply the covariance function estimator, designed
for outside the diagonal set, for the closest point on the boundary of the
diagonal band. The rate (45) indicates that in order to make this a suitable
estimator, $\mathfrak{d}$ should have a rate of decrease slower than the power
$(2\alpha+1)^{-1}$ of the optimal rate achievable by a nonparametric estimator
of the covariance. If $\hat{\alpha}$ is the estimate of $\alpha$, we can take

$\mathfrak{d}=\left\\{N^{-2}\sum_{i=1}^{N}(1/M_{i})\right\\}^{c}\quad\text{
for some
}\quad\frac{2\hat{\alpha}}{\\{2\hat{\alpha}+1\\}^{2}}<c<\frac{1}{2\hat{\alpha}+1},$
(46)

for instance $c=\\{2\hat{\alpha}+1/2\\}/\\{2\hat{\alpha}+1\\}^{2}$.

##  5 Empirical study

###  5.1 Implementation aspects

The risks $\mathcal{R}_{\mu}$ and $\mathcal{R}_{\Gamma}$ defined in (30) and
(41), respectively, depend on $L^{2}_{\delta}$ and the conditional variance
bound $\sigma^{2}_{\max}$. In view of (15), if $t_{2}$ is the midpoint of
$[t_{1},t_{3}]$,

$L_{\delta}^{2}\approx\frac{1}{2}\left(\frac{\theta_{\delta}(t_{2},t_{3})}{\lvert
t_{3}-t_{2}\rvert^{2(\alpha-\delta)}}+\frac{\theta_{\delta}(t_{1},t_{2})}{\lvert
t_{2}-t_{1}\rvert^{2(\alpha-\delta)}}\right),$ (47)

provided $t_{3}-t_{1}=\Delta_{*}/2$ is small. Given the estimate
$\hat{\alpha}_{t}$ and the estimators $\hat{\theta}_{\delta}(t_{2},t_{3})$ and
$\hat{\theta}_{\delta}(t_{1},t_{2})$ as in (17), with
$\delta=\lfloor\hat{\alpha}_{t}\rfloor$, we then define a natural estimator of
$L^{2}_{\delta}$ as

$\widehat{L}_{\delta}^{2}\approx\frac{1}{2}\left(\frac{\hat{\theta}_{\lfloor\hat{\alpha}_{t}\rfloor}(t_{2},t_{3})}{\lvert
t_{3}-t_{2}\rvert^{2(\hat{\alpha}_{t}-\lfloor\hat{\alpha}_{t}\rfloor)}}+\frac{\hat{\theta}_{\lfloor\hat{\alpha}_{t}\rfloor}(t_{1},t_{2})}{\lvert
t_{2}-t_{1}\rvert^{2(\hat{\alpha}_{t}-\lfloor\hat{\alpha}_{t}\rfloor)}}\right).$
(48)

To estimate the conditional variance bound, let us first consider the case
where $\sigma^{2}(t,x)$ does not depend on $x$. In this case, one can compute

$\widehat{\sigma}^{2}(t)=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{2|\mathcal{S}_{i}|}\sum_{m\in\mathcal{S}_{i}}\left[Y_{m}^{(i)}-Y_{m-1}^{(i)}\right]^{2},$
(49)

where $\mathcal{S}_{i}$ is a subset of indices $m$ for the $i-$th trajectory.
When the variance of the errors is considered constant, $\mathcal{S}_{i}$ can
be the set $\\{2,3,\ldots,M_{i}\\}$. When the variance depends on ${t}$, one
could define $\mathcal{S}_{i}$ as the set of indices corresponding to the
$K_{0}$ values $T_{m}^{(i)}$ closest to $t$. The theory allows a choice such
as
$K_{0}=\lfloor\hat{\mathfrak{m}}\exp\left(-\\{\log\log\hat{\mathfrak{m}}\\}^{2}\right)\rfloor$,
with $\hat{\mathfrak{m}}=N^{-1}\sum_{i=1}^{N}M_{i}$. Then $\sigma^{2}_{\max}$
could be $\max_{t\in\mathcal{T}}\widehat{\sigma}^{2}(t)$, and this choice was
used in our empirical investigation. When the variance of the errors also
depends on the realizations $X_{u}$, in general it is no longer possible to
consistently estimate $\sigma^{2}(u,X_{u})$. However, the simple inspection of
the squared differences between the observations $Y_{m}^{(i)}$ and the
presmoothed values allows to build a bound $\sigma^{2}_{\max}$. In our
simulation experiments, $\sigma^{2}_{\max}$ obtained from the estimate (49),
with $\mathcal{S}_{i}$ given by $K_{0}$, works well.

###  5.2 Simulation design

Our simulation study is based on the Household Active Power Consumption
dataset which was sourced from the UC Irvine Machine Learning Repository
(https://archive.ics.uci.edu/ml/datasets/Individual+household+electric+power+consumption).
This dataset contains diverse energy related features gathered in a house
located near Paris, every minute between December 2006 and November 2010. In
total, its represents around $2$ million data points. Here, we are only
interested in the daily voltage and we only consider the days without missing
values in the measurements. The extracted dataset contains $708$ voltage
curves with an uniform common design with $1440$ points. In Figure
LABEL:fig:sample_from_data in the Supplementary Material, we plot the
logarithms of the ordered eigenvalues  of a smoothed version of the empirical
covariance matrix. Moreover, we present the empirical mean curve and three
randomly selected curves from the sample. The time is normalized such that
$\mathcal{T}=[0,1]$. We also present in that figure the estimates of the local
regularity obtained with the $708$ voltage curves. Most of the estimates are
smaller than 1.

The 708 voltage curves are used to build a mean function $\mu(\cdot)$, and a
covariance function $\Gamma(\cdot,\cdot)$. We also derive a conditional
variance function $\sigma^{2}(t,x)$ for the noise which we plot in Figure
1(a). Next, we generate samples of independent trajectories from the Gaussian
process characterized by these $\mu(\cdot)$ and $\Gamma(\cdot,\cdot)$. Their
local regularity is approximately equal to 0.7. Finally, we add the
heteroscedastic noise. A random sample of curves generated according to our
simulation setup are plotted in Figure 1(b). The details on the construction
of our simulation setup are provided in the Supplementary Material.

We consider four experiments, each of them replicated 500 times. In
_Experiment 1_ , $N\in\\{40,100,200\\}$, $\mathfrak{m}\in\\{40,100,200\\}$,
and $(1-p)\mathfrak{m}\leq M_{i}\leq(1+p)\mathfrak{m}$ with $p=0.2$.
_Experiment 2_ is the same as _Experiment 1_ except that $p=0.5$. In
_Experiment 3_ we modify _Experiment 1_ by multiplying the conditional
variance of the error terms by $2$. Finally, in _Experiment 4_ we consider the
setup from _Experiment 1_ but with a higher regularity for the mean function.
In the Supplementary Material we present plots of the mean functions and
trajectories from _Experiment 4_.

For the presmoothing, we use the locpoly function in R. Given that most of the
local regularity estimates are between 0 and 1, our mean and covariance
functions estimators are built with the NW smoother. Moreover, the value
$k_{0}$ in (23) is set equal to 2, and we use the biweight kernel
$K(t)=(15/16)\left(1-t^{2}\right)^{2}\mathbf{1}_{[-1,1]}(t)$. The estimates
$\hat{\alpha}_{t}$ and the estimates of the mean and covariance functions are
obtained using the same data. That means we did not use a _learning sample_
for $\hat{\alpha}_{t}$. An implementation of the method used in the four
experiments is available as a R package on Github at the URL adress:
https://github.com/StevenGolovkine/funestim. More details on the
implementation are also available in the Supplementary Material.

###  5.3 Mean estimation

Concerning the estimation of the mean, the estimates $\hat{\alpha}_{t}$ are
computed using (17), on a uniform grid of 50 points $t_{2}$ between 0.05 and
0.95, with $t_{3}-t_{1}=\Delta_{*}/2=\exp(-\log^{1/2}(\hat{\mathfrak{m}}))$.
For each value of the 50 estimates $\hat{\alpha}_{t}$, we compute the optimal
bandwidths $h^{*}_{\mu}$ by minimization with respect to $h$ over a
logarithmic grid of 151 points.

Our mean estimator is compared to that of Cai and Yuan,  (2011), denoted
$\widehat{\mu}_{CY}$, and Zhang and Wang,  (2016), denoted
$\widehat{\mu}_{ZW}$. To compute $\widehat{\mu}_{CY}$, we use smooth.splines
in R, with the $M_{1}+\ldots+M_{N}$ data points
$({Y^{(i)}_{m}},{T^{(i)}_{m}})$. To obtain $\widehat{\mu}_{ZW}$, we use the R
package fdapace, see Carroll et al.,  (2021). To compare the accuracy of the
estimators, we use the ISE risk with respect to the target. For any $f$ and
$g$ real-valued functions defined on $[0,1]$, the ISE is defined as

$\text{ISE}(f,g)=\lVert f-g\rVert^{2}=\int_{[0,1]}\\{f(t)-g(t)\\}^{2}dt.$ (50)

We approximate the integral using the mean estimates on a uniform grid of 101
points and the trapezoidal rule. For each configuration $N$, $\mathfrak{m}$,
$p$, and each of the 500 samples, we compute the $\text{ISE}s$ with respect to
the infeasible $\widetilde{\mu}$, and the $\text{ISE}s$ with respect to the
mean function $\mu$ used for generating the samples. The 101 bandwidth values
used for our estimator are obtained from the 50 optimal bandwidths
$h^{*}_{\mu}$ by linear interpolation.

The results obtained in _Experiment 1_ are plotted in the Figure 2, on a
logarithmic scale. The results for the other three experiments are relegated
to the Supplementary Material. Our mean function estimator reveals good
performance. It provides a much more accurate estimate of the infeasible
empirical mean function $\widetilde{\mu}$. When compared to the true mean
function, except the cases $N\in\\{100,200\\}$ and $\mathfrak{m}=40$, our
estimator outperforms the competitors. In that cases, our estimator,
$\widehat{\mu}_{CY}$ and $\widehat{\mu}_{ZW}$ have similar performance. The
fact that the advantage of our estimator wanes in these cases is explained by
the poorer performance of the infeasible empirical mean function. Similar
conclusions could be drawn from the _Experiments 2_ to _4_. In the setup with
a more regular mean function, the advantage of our estimator diminishes for
all pairs $N,\mathfrak{m}$.

(a) The conditional variance function $\sigma^{2}(t,x)$

(b) The mean function $\mu(t)$ (bold line) and three randomly generated noisy
trajectories

Figure 1: The simulation setup in _Experiment 1_

(a) ISE with respect to the empirical mean $\tilde{\mu}$ in each simulated
sample

(b) ISE with respect to the true mean function $\mu$

Figure 2: Results from _Experiment 1_ on the log-scale

###  5.4 Covariance estimation

Concerning the estimation of the covariance, the estimates $\hat{\alpha}_{t}$
are computed using (17), on a uniform grid of 10 points $t_{2}$ between 0.05
and 0.95, with
$t_{3}-t_{1}=\Delta_{*}/2=\exp(-\log^{1/2}(\hat{\mathfrak{m}}))$. For each
values $s,t$ in the grid, we compute the optimal bandwidths
$h_{\Gamma}^{*}(s,t)$ by minimization over a logarithmic grid of 41 points
between 0.01 and 0.1. Our covariance estimator is compared to the ones from
Cai and Yuan,  (2010), denoted $\widehat{\Gamma}_{CY}$, and from Zhang and
Wang,  (2016), denoted $\widehat{\Gamma}_{ZW}$. We compute
$\widehat{\Gamma}_{CY}$ using the R package ssfcov, see Cai and Yuan,  (2010).
For $\widehat{\Gamma}_{ZW}$, we use the R package fdapace, see Carroll et al.,
(2021). The risk we consider is defined as

$\text{ISE}(f,g)=\lVert
f-g\rVert^{2}=\int_{[0,1]}\int_{[0,1]}\\{f(s,t)-g(s,t)\\}^{2}dsdt,$ (51)

approximated by the trapezoidal rule applied with a uniform grid of $101\times
101$ points. The $101\times 101$ bandwidth values used for our estimator are
obtained from the $10\times 10$ optimal bandwidths $h^{*}_{\Gamma}$ by linear
interpolation. For each configuration $N$, $\mathfrak{m}$, $p$, and each of
the 500 samples, we compute the $\text{ISE}s$ with respect to the infeasible
$\widetilde{\Gamma}$, and the $\text{ISE}s$ with respect to the covariance
function $\Gamma$ used for generating the samples.

The results obtained in _Experiment 1_ are plotted in the Figure 3, on a
logarithmic scale. The results for the other three experiments are relegated
to the Supplementary Material. We were not able to calculate
$\widehat{\Gamma}_{CY}$ when $M=200$ and $N=200$, and therefore some cases are
not reported. Our estimator provides the most accurate approximation of
$\widetilde{\Gamma}$. Our estimator and $\widehat{\Gamma}_{ZW}$ show better
accuracy for estimating $\Gamma$ in all cases considered. Meanwhile, our
estimator performs similarly or slightly better than $\widehat{\Gamma}_{ZW}$.
The advantage of our approach increases with $N$. Let us point out the clear
advantage for our estimator in terms of computation times. For instance, with
$N=100$ and $M=100$, our covariance estimator is computed about three times
faster that of Zhang and Wang,  (2016) (with fixed bandwidth) and about thirty
times faster than that of Cai and Yuan,  (2010). More details on the
computation time are provided in the Supplementary Material.

(a) ISE with respect to the empirical covariance $\tilde{\Gamma}$ in each
simulated sample

(b) ISE with respect to the true covariance function $\Gamma$

Figure 3: Results from _Experiment 1_ on the log-scale

##  6 Discussion and conclusions

We propose new nonparametric estimators for the mean and covariance function.
They are built using a novel ‘smoothing first, then estimate’ strategy based
on univariate local polynomials. The main novelty comes from the fact that the
optimal bandwidths for the local polynomials are selected by minimization of
suitable penalized quadratic risks. The penalized risks for the mean and the
covariance functions are quite similar, could be easily built from data and be
optimized on a grid of bandwidths. What distinguishes them from the usual sum
between the squared bias and the variance, is a penalty for the fact that not
all the curves have enough observation points to be included in the final
estimator. Removing curves from the nonparametric estimators of the mean and
covariance functions is an aspect which characterizes practically all
smoothing-based approaches. Indeed, to entirely benefit from the replication
feature of functional data, one has to determine the amount of smoothing for
the mean and covariance estimation using all the curves. In this case some
curves could present too few observation points and thus will be dropped. This
is more likely to happen in the so-called sparse regime. To our best
knowledge, our bandwidth choice procedure is the first attempt to explicitly
account for this aspect. We thus build estimators which achieve the optimal
rates of convergence in a completely adaptive, data-driven way. The
theoretical results are derived under very mild conditions. In particular, the
curves could be observed with heteroscedastic errors at discrete observations
points. These points could be common to all curves or they could change
randomly from one curve to another. In the case of the common observation
points, our procedure automatically chooses between smoothing and
interpolation, the latter being known to be rate optimal, but is not
necessarily the best solution with finite samples.

Our nonparametric estimation approach relies on a probabilistic concept of
local regularity for the sample paths of the process generating the curves. In
some common examples, this local regularity is related to the polynomial
decrease rate of the eigenvalues of the covariance operator, a characteristic
of the data generating process widely used in the literature and usually
supposed to be known. The local regularity also determines the regularity of
the trajectories, the usual concept used in nonparametric regression. It is
well-known that the optimal rates, in the minimax sense, for estimating the
mean and covariance functions depend on the regularity of the trajectories.
Moreover, the so-called sparse and dense regimes, commonly invoked in
functional data analysis, are defined using the regularity of the
trajectories, which usually is supposed to be known. We therefore propose a
novel simple estimator of the local regularity of the process and we use it to
build our penalized quadratic risk. Applied to real data, the local regularity
estimator reveals that the regularity of the trajectories could be quite far
from what is usually assumed in the existing theoretical contributions.

Our method performs quite well in simulations and outperforms the main
competitors when the mean and covariance functions have a regularity close to
that of the trajectories. The approach is still satisfactory when the mean or
the covariance functions are more regular than the trajectories. The reason is
that, in some sense, our nonparametric estimators are as close as possible to
the empirical mean and covariance, respectively, which are the ideal
estimators if the trajectories were observed at any point without error. In
the case where the mean and covariance function are smoother than the
trajectories, our penalized quadratic risk should be built using the mean or
covariance functions’ regularity instead of the trajectories’ regularity.
However, for now, the estimation of the regularity of the mean or covariance
function remains an open problem.

##  Appendix A Details on the definition (30)

To explain our empirical risk bound $\mathcal{R}_{\mu}(t;h)$ defined in (30),
let

$\widetilde{\mu}_{W}(t;h)=\frac{1}{\mathcal{W}_{N}(t;h)}\sum_{i=1}^{N}w_{i}(t;h)X^{(i)}_{t},$
(A.1)

be the unfeasible estimator of $\mu(\cdot)$ using only the curves for which
${\widehat{X}^{(i)}}_{t}$ is well-defined. Using the definitions, see (8), we
then have,

$\displaystyle\mathbb{E}_{M,T}$
$\displaystyle\left[\left\\{\widetilde{\mu}_{N}(t)-\widehat{\mu}_{N}(t;h)\right\\}^{2}\right]$
$\displaystyle=\mathbb{E}_{M,T}\left[\left\\{\widetilde{\mu}_{N}(t)-\widetilde{\mu}_{W}(t;h)-\frac{1}{\mathcal{W}_{N}(t;h)}\sum_{i=1}^{N}w_{i}(t;h)\left(B^{(i)}_{t}+V^{(i)}_{t}\right)\right\\}^{2}\right]$
$\displaystyle\leq
2\mathbb{E}_{M,T}\left[\left\\{\widetilde{\mu}_{N}(t)-\widetilde{\mu}_{W}(t;h)\right\\}^{2}\right]+2\mathbb{E}_{M,T}\left[\left\\{\frac{1}{\mathcal{W}_{N}(t;h)}\sum_{i=1}^{N}w_{i}(t;h)\left(B^{(i)}_{t}+V^{(i)}_{t}\right)\right\\}^{2}\right]$
$\displaystyle=:2E_{1}+2E_{2}.$

In the following, for simplicity, we write $w_{i}$ and $\mathcal{W}_{N}$
instead of $w_{i}(t;h)$ and $\mathcal{W}_{N}(t;h)$, respectively. Since

$\widetilde{\mu}_{W}(t;h)-\widetilde{\mu}_{N}(t)=\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}\left\\{X_{t}^{(i)}-\mu(t)\right\\}\left\\{w_{i}-\frac{\mathcal{W}_{N}}{N}\right\\},$

and the trajectories are drawn independently, we have

$E_{1}=\frac{Var(X_{t})}{\mathcal{W}_{N}^{2}}\sum_{i=1}^{N}\left\\{w_{i}-\frac{\mathcal{W}_{N}}{N}\right\\}^{2}\\\
=Var(X_{t})\left\\{\frac{1}{\mathcal{W}_{N}}-\frac{1}{N}\right\\}.$ (A.2)

For $E_{2}$, let us first look at the bias part. By the arguments used in the
proof of Proposition 1.13 in Tsybakov,  (2009) and the Cauchy-Schwarz
inequality, we have

$\displaystyle\mathbb{E}_{M,T}\left[\left\\{\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}B^{(i)}_{t}\right\\}^{2}\right]$
$\displaystyle\leq\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}\mathbb{E}_{M,T}\left[\left\\{B^{(i)}_{t}\right\\}^{2}\right]$
$\displaystyle\mkern-126.0mu\leq\frac{1}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\times\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}\left\\{\sum_{m=1}^{M_{i}}\left|W_{m}^{(i)}(t)\right|\right.$
$\displaystyle\qquad\left.\times\sum_{m=1}^{M_{i}}\mathbb{E}\left(\left\\{\nabla^{\lfloor\hat{\alpha}_{t}\rfloor}{X^{(i)}}({T^{(i)}_{m}})-\nabla^{\lfloor\hat{\alpha}_{t}\rfloor}X^{(i)}_{t}\right\\}^{2}\mid\mathcal{T}_{obs}^{(i)}\right)|W_{m}^{(i)}(t)|\right\\}$
$\displaystyle\mkern-126.0mu=\frac{L_{\delta}^{2}h^{2\hat{\alpha}_{t}}\\{1+o_{\mathbb{P}}(1)\\}}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\times\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}\left\\{\sum_{m=1}^{M_{i}}\left|W_{m}^{(i)}(t)\right|\times\sum_{m=1}^{M_{i}}\left|\frac{{T^{(i)}_{m}}-t}{h}\right|^{2\hat{\alpha}_{t}}|W_{m}^{(i)}(t)|\right\\}$
$\displaystyle\mkern-126.0mu\eqqcolon\frac{L_{\delta}^{2}h^{2\hat{\alpha}_{t}}}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\times\overline{C}_{1}(t;h,2\hat{\alpha}_{t})\times\\{1+o_{\mathbb{P}}(1)\\}.$

For the first equality we use the condition (H2) on the event
$\\{\lfloor\hat{\alpha}_{t}\rfloor=\delta\\}$. By Theorem 1, this event has an
exponentially small probability. When the kernel function is supported on
$[-1,1]$, we could use the bound

$\mathbb{E}_{M,T}\left[\left\\{\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}B^{(i)}_{t}\right\\}^{2}\right]\\\
\leq\frac{L_{\delta}^{2}h^{2\hat{\alpha}_{t}}\\{1+o_{\mathbb{P}}(1)\\}}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\times\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}\left\\{\sum_{m=1}^{M_{i}}\left|W_{m}^{(i)}(t)\right|\right\\}^{2}.$
(A.3)

In the case $\lfloor\hat{\alpha}_{t}\rfloor=0$, with the Nadaraya-Watson
estimator, we could use the more refined bound

$\displaystyle\mathbb{E}_{M,T}\left[\left\\{\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}B^{(i)}_{t}\right\\}^{2}\right]$
$\displaystyle\mkern-90.0mu\leq\frac{L_{0}^{2}h^{2\hat{\alpha}_{t}}\\{1+o_{\mathbb{P}}(1)\\}}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\times\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}\left\\{\sum_{m=1}^{M_{i}}\left|\frac{{T^{(i)}_{m}}-t}{h}\right|^{2\hat{\alpha}_{t}}W_{m}^{(i)}(t)\right\\}$
$\displaystyle\mkern-90.0mu\approx\frac{L_{0}^{2}h^{2\hat{\alpha}_{t}}\\{1+o_{\mathbb{P}}(1)\\}}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\times\int|u|^{2\hat{\alpha}_{t}}K(u)du.$

Using the equivalent kernels idea, see section 3.2.2 in Fan and Gijbels,
(1996), the bound on the last line of the last display could be used in the
case of local linear estimators. To complete the bound for $E_{2}$, note that
by construction, $\mathbb{E}_{M,T}\left\\{V_{t}^{(i)}B_{t}^{(i)}\right\\}=0$
and

$\mathbb{E}_{M,T}\left\\{V_{t}^{(i)}B_{t}^{(j)}\right\\}=\mathbb{E}_{M,T}\left\\{V_{t}^{(i)}V_{t}^{(j)}\right\\}=0,\quad\forall
1\leq i\neq j\leq N.$

Up to negligible terms, we can then write

$\displaystyle E_{2}$
$\displaystyle\leq\mathbb{E}_{M,T}\left[\left\\{\frac{1}{\mathcal{W}_{N}}\sum_{i=1}^{N}w_{i}B^{(i)}_{t}\right\\}^{2}\right]+\frac{1}{\mathcal{W}_{N}^{2}}\sum_{i=1}^{N}w_{i}\mathbb{E}_{M,T}\left[\left\\{V_{t}^{(i)}\right\\}^{2}\right]$
(A.4) $\displaystyle\leq
h^{2\hat{\alpha}_{t}}\frac{L_{\delta}^{2}}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\overline{C}_{1}(t;h,2\hat{\alpha}_{t})+\frac{\sigma^{2}_{\max}}{\mathcal{W}_{N}^{2}}\sum_{i=1}^{N}w_{i}\left\\{\max_{m}\left|W_{m}^{(i)}(t;h)\right|\times\sum_{m=1}^{M_{i}}\left|W_{m}^{(i)}(t;h)\right|\right\\}$
(A.5)
$\displaystyle=h^{2\hat{\alpha}_{t}}\frac{L_{\delta}^{2}}{\lfloor\hat{\alpha}_{t}\rfloor!^{2}}\overline{C}_{1}(t;h,2\hat{\alpha}_{t})+\frac{\sigma^{2}_{\max}}{\mathcal{W}_{N}^{2}}\sum_{i=1}^{N}w_{i}\frac{c_{i}(t;h)}{\mathcal{N}_{i}(t;h)},$
(A.6)

where $\sigma_{\max}$ is a bound for the function $\sigma(t,x)$ in (6).

##  Appendix B Proofs

###### Proof of Theorem 2.

For simplicity, $\widehat{\mu}_{N}^{*}$ is built with $k_{0}=1$ and the
uniform kernel. The lines below adapt to other choices of $k_{0}$ and
$K(\cdot)$ at the price of more involved technical arguments. First, let us
prove that

$\frac{1}{\mathcal{W}_{N}(t;h)}-\frac{1}{N}\leq\max\left\\{h^{2\alpha_{t}},\;\mathcal{N}^{-1}_{\mu}(t;h)\right\\}O_{\mathbb{P}}(1),$
(A.7)

provided $N\mathfrak{m}h\rightarrow\infty$ and $h\rightarrow 0$. Note that if
$\mathcal{W}_{N}(t;h)=0$, then necessarily $\mathcal{N}_{\mu}(t;h)=0$. The
property (A.7) is implied by the following ones: there exist two constants
$\mathfrak{c}_{1},\mathfrak{c}_{2}>0$ such that

$\frac{1}{\mathcal{N}_{\mu}(t;h)}\geq\frac{\mathfrak{c}_{1}\\{1+o_{\mathbb{P}}(1)\\}}{N\mathfrak{m}h},$
(A.8)

and

$\frac{1}{\mathcal{W}_{N}(t;h)+1}\leq\max\left\\{\frac{1}{N+1},\;\frac{\mathfrak{c}_{2}\\{1+o_{\mathbb{P}}(1)\\}}{N\mathfrak{m}h}\right\\}.$
(A.9)

Indeed, the latter two properties imply

$\displaystyle\frac{1}{\mathcal{W}_{N}(t;h)+1}-\frac{1}{N+1}$
$\displaystyle=\max\left\\{0,\;\frac{1}{N}\left(\frac{1}{\mathfrak{m}h}-1\right)\right\\}O_{\mathbb{P}}(1)$
$\displaystyle\leq\max\left\\{h^{2\alpha_{t}},\;\frac{1}{\mathcal{N}_{\mu}(t;h)}\right\\}O_{\mathbb{P}}(1),$

and from this we obtain (A.7) because

$\\{\mathcal{W}_{N}(t;h)+1\\}^{-1}-\\{N+1\\}^{-1}=\left\\{\mathcal{W}_{N}^{-1}(t;h)-N^{-1}\right\\}\frac{N\\{N+1\\}}{\mathcal{W}_{N}(t;h)\\{\mathcal{W}_{N}(t;h)+1\\}}.$

To justify (A.8) and (A.9), we omit the arguments $t$ and $h$. For instance,
we write $\mathcal{W}_{N}$ and $\mathcal{N}_{\mu}$ instead of
$\mathcal{W}_{N}(t;h)$ and $\mathcal{N}_{\mu}(t;h)$, respectively.

Using the fact that the harmonic mean is less than or equal to the mean we
obtain

$\frac{1}{\mathcal{N}_{\mu}(t;h)}\geq\frac{c_{i}}{\sum_{i=1}^{N}w_{i}\mathcal{N}_{i}},$
(A.10)

with $c_{i}=c_{i}(t;h)$ and $\mathcal{N}_{i}=\mathcal{N}_{i}(t;h)$ defined in
(25) and (28), respectively. In the case of $\widehat{\alpha}_{t}<1$, for all
$i$ we have $c_{i}\equiv 1$. To justify (A.8) it suffices to prove that there
exists a sequence of variables $C_{N}$ such that almost surely
$0<\liminf_{N}C_{N}<\limsup_{N}C_{N}<\infty$, and

$C_{N}\\{1+o_{\mathbb{P}}(1)\\}=\frac{\sum_{i=1}^{N}w_{i}\mathcal{N}_{i}}{N\mathfrak{m}h}.$
(A.11)

For this purpose, let us notice that in the case of an NW estimator with a
uniform kernel, when $k_{0}=1$,

$\sum_{i=1}^{N}w_{i}\mathcal{N}_{i}=\sum_{i=1}^{N}\sum_{m=1}^{M_{i}}\mathbf{1}\\{|{T^{(i)}_{m}}-t|\leq
h\\}.$

Thus $\sum_{i=1}^{N}w_{i}\mathcal{N}_{i}$ is a binomial random variable with
$M_{1}+\ldots+M_{N}$ trials and success probability

$p(h)=p(h;t)=\int_{t-h}^{t+h}f_{T}(s)ds\approx 2f_{T}(t)h,$

provided $f_{T}$, the density of the ${T^{(i)}_{m}}$, is continuous at $t$.
Since $N\mathfrak{m}h\rightarrow\infty$, by Bernstein’s inequality the
property (A.11) holds true with

$C_{N}=2f_{T}(t)\frac{M_{1}+\ldots+M_{N}}{N\mathfrak{m}}.$

Condition (33) guarantees that almost surely $C_{N}$ stays away from zero and
infinity. The property (A.8) follows.

For justifying (A.9), let us first recall that by definition,
$\mathcal{W}_{N}(t;h)\leq N$. Thus it remains to show that

$\frac{1}{\mathcal{W}_{N}(t;h)+1}\leq\frac{\mathfrak{c}_{2}\\{1+o_{\mathbb{P}}(1)\\}}{N\mathfrak{m}h}.$

For this purpose, it suffices to show that there exists a positive constant
$\mathfrak{c}_{\mathcal{W}}>0$ such that

$\min\\{N,\mathfrak{c}_{\mathcal{W}}N\mathfrak{m}h\\{1+o_{\mathbb{P}}(1)\\}\\}\leq\mathcal{W}_{N}(t;h).$
(A.12)

Let $\mathbb{P}_{M}(\cdot)=\mathbb{P}(\cdot\mid M_{i},1\leq i\leq N)$. We now
note that

$p_{i}=p_{i}(t):=\mathbb{P}_{M}(w_{i}=1)=1-[1-p(h)]^{M_{i}}.$

Next, note that the variable $\mathcal{W}_{N}(t;h)$ is a sum of $N$
independent Bernoulli variables with probabilities $p_{i}$. In the case where
$\mathfrak{m}h\geq\underline{c}>0$ for some constant $\underline{c}$, since
$\\{1-p(h)\\}^{-1/p(h)}>e$ for any $h<1/2$, we can write

$p_{i}=1-[1-p(h)]^{M_{i}}\geq
1-[1-p(h)]^{c_{L}\mathfrak{m}}>1-\exp(-c_{L}\mathfrak{m}p(h))\\\ \approx
1-\exp(-2f_{T}(t)c_{L}\mathfrak{m}h)\geq
1-\exp(-2f_{T}(t)c_{L}\underline{c}),$ (A.13)

and the approximations are valid as soon as $h\rightarrow 0$. Thus, in the
case where $\mathfrak{m}h\geq\underline{c}>0$, $\liminf p_{i}>0$. Then,
Bernstein’s inequality yields the property (A.12). More precisely,
$\mathcal{W}_{N}(t;h)$ will increase at a rate at least as fast as
$[1-\exp(-2f_{T}(t)c_{L}\underline{c})]N$. In the case
$\mathfrak{m}h\rightarrow 0$ we have

$p_{i}=\mathbb{P}_{M}(w_{i}=1)\approx 1-\exp(-M_{i}p(h))\approx M_{i}p(h)\geq
2c_{L}f_{T}(t)\mathfrak{m}h.$

Again, Bernstein’s inequality guarantees that

$\frac{\mathcal{W}_{N}(t;h)}{p_{1}+\ldots+p_{N}}=1+o_{\mathbb{P}}(1).$

Since $p_{1}+\ldots+p_{N}\geq 2c_{L}f_{T}(t)\times N\mathfrak{m}h$, condition
(A.12) follows with $\mathfrak{c}_{\mathcal{W}}=2c_{L}f_{T}(t)$.

Finally, to complete the proof in the independent design case, it suffices
first to notice that from above, we have

$\max\left\\{h^{2\alpha_{t}},\;\mathcal{N}^{-1}_{\mu}(t;h)\right\\}=O_{\mathbb{P}}\left(h^{2\alpha_{t}}+(N\mathfrak{m}h)^{-1}\right),$

which is minimized by $h$ with the rate
$(N\mathfrak{m})^{-2\alpha_{t}/\\{2\alpha_{t}+1\\}}$. Next, condition (32)
guarantees that $\log(N\mathfrak{m})/\log(\mathfrak{m})$ is bounded, and thus
$h^{2\hat{\alpha}_{t}}=h^{2\alpha_{t}}\\{1+o_{\mathbb{P}}(1)\\}$ whenever
$\hat{\alpha}_{t}-\alpha_{t}=o_{\mathbb{P}}(\log^{-1}(\mathfrak{m}))$. For the
rate of $\widehat{\mu}_{N}^{*}(t)-\mu(t)$ we simply add the rate of
$\widetilde{\mu}_{N}(t)-\mu(t)$.

With a common design $\mathcal{W}_{N}(t;h)$ can only take the values 0 or
$N^{-1}$. Thus the penalty introduced by $\mathcal{W}_{N}(t;h)^{-1}-N^{-1}$
plays a different role. It constrains the bandwidth to be greater than or
equal to the lengths of the intervals $[T_{m}^{(i)},T_{m+1}^{(i)}]$ including
$t$. By condition (34), this means that the rate of convergence of
$\widehat{\mu}_{N}^{*}(t)-\widetilde{\mu}_{N}(t)$ could not be faster than
$O_{\mathbb{P}}\left(\mathfrak{m}^{-\;2\alpha_{t}}\right)$. This aspect is
automatically included in the definition of $\mathcal{R}_{\mu}(t;h)$ because,
under the constraint $\mathfrak{m}h\geq c_{L}/C_{U}$,

$O_{\mathbb{P}}\left(\max\left\\{h^{2\alpha_{t}},\;\mathcal{N}^{-1}_{\mu}(t;h)\right\\}\right)=O_{\mathbb{P}}\left(\max\left\\{h^{2\alpha_{t}},\;(N\mathfrak{m}h)^{-1},\;N^{-1}\right\\}\right)=O_{\mathbb{P}}\left(\mathfrak{m}^{-\;2\alpha_{t}}\right).$

Finally, $\alpha_{t}$ can be replaced by $\hat{\alpha}_{t}$ using the
arguments from the independent design case. ∎

The proof of Theorem 3 follows the lines of that of the proof of Theorem 2 and
is left to the Supplementary Material.

###### Proof of equation (45)..

To prove the rate of $\widetilde{D}(\mathfrak{d})$, we use the assumptions:
$\sup_{t\in\mathcal{T}}\mathbb{E}(X_{t}^{4})<\infty$ and there exists a
constant $c$ such that

$\mathbb{E}(\\{X_{s}-X_{t}\\}^{4})\leq
c\mathbb{E}^{2}(\\{X_{s}-X_{t}\\}^{2}),\quad\forall s,t\in\mathcal{T}.$

Omitting the integration domain, we can now write

$\displaystyle\mathbb{E}[\widetilde{D}(\mathfrak{d})]$ $\displaystyle\leq
2\iint\mathbb{E}\left[\left(\frac{1}{N}\sum_{i=1}^{N}\left\\{X^{(i)}_{s}X^{(i)}_{t}-X^{(i)}_{\frac{s+t-\mathfrak{d}}{2}}X^{(i)}_{\frac{s+t+\mathfrak{d}}{2}}\right\\}\right)^{2}\right]dsdt$
$\displaystyle\quad+2\\!\iint\\!\mathbb{E}\\!\left[\left(\\!\left\\{\\!\frac{1}{N}\sum_{i=1}^{N}X^{(i)}_{s}\\!\right\\}\\!\left\\{\\!\frac{1}{N}\sum_{i=1}^{N}X^{(i)}_{t}\\!\right\\}-\\!\left\\{\\!\frac{1}{N}\sum_{i=1}^{N}X^{(i)}_{\frac{s+t-\mathfrak{d}}{2}}\\!\right\\}\left\\{\\!\frac{1}{N}\sum_{i=1}^{N}X^{(i)}_{\frac{s+t+\mathfrak{d}}{2}}\\!\right\\}\right)^{\\!2}\right]\\!dsdt$
$\displaystyle=2\widetilde{D}_{1}+2\widetilde{D}_{2}.$

By the Cauchy-Schwarz inequality,

$\mathbb{E}\left[\left\\{X^{(i)}_{s}X^{(i)}_{t}-X^{(i)}_{\frac{s+t-\mathfrak{d}}{2}}X^{(i)}_{\frac{s+t+\mathfrak{d}}{2}}\right\\}^{2}\right]=O\left(\mathfrak{d}^{2\alpha}\right),$

where $\alpha$ is the local regularity on ${\mathcal{T}_{k}^{\circ}}$ to which
$s$, $t$, $(s+t+\mathfrak{d})/2$ and $(s+t-\mathfrak{d})/2$ eventually belong.
We deduce that $\widetilde{D}_{1}=O\left(\mathfrak{d}^{2\alpha+1}\right).$ On
the other hand, by repeatedly application of the Cauchy-Schwarz inequality and
the moment conditions imposed above,

$\displaystyle\widetilde{D}_{2}$ $\displaystyle\leq
2\iint\mathbb{E}^{1/2}\\!\left[\\!\left(\frac{1}{N}\sum_{i=1}^{N}\\{X^{(i)}_{s}\\!-X^{(i)}_{\frac{s+t-\mathfrak{d}}{2}}\\}\right)^{\\!\\!4}\right]\\!\mathbb{E}^{1/2}\left[\left(\frac{1}{N}\sum_{i=1}^{N}X^{(i)}_{t}\right)^{\\!\\!4}\right]dsdt$
$\displaystyle\quad+2\iint\\!\mathbb{E}^{1/2}\left[\\!\left(\frac{1}{N}\sum_{i=1}^{N}X^{(i)}_{\frac{s+t-\mathfrak{d}}{2}}\right)^{\\!4}\right]\\!\mathbb{E}^{1/2}\\!\left[\\!\left(\frac{1}{N}\sum_{i=1}^{N}\\{X^{(i)}_{t}\\!-X^{(i)}_{\frac{s+t+\mathfrak{d}}{2}}\\}\right)^{\\!\\!4}\right]\\!dsdt$
$\displaystyle=O\left(\mathfrak{d}^{2\alpha+1}\right).$

Gathering facts, we deduce that
$\widetilde{D}(\mathfrak{d})=O_{\mathbb{P}}\\!\left(\mathfrak{d}^{2\alpha+1}\right).$
∎

## Acknowledgements

The authors wish to thank Groupe Renault and the ANRT (French National
Association for Research and Technology) for their financial support via the
CIFRE convention no. 2017/1116. Valentin Patilea gratefully acknowledges
support from the Joint Research Initiative “Models and mathematical processing
of very large data” under the aegis of Risk Foundation, in partnership with
MEDIAMETRIE and GENES, France.

## Supplementary Material

In the Supplementary Material, we provide some additional technical arguments,
proofs and simulation results. In Section LABEL:sec:reg_loc, we provide
details on the relation between local regularity and the decreasing of the
eigenvalues of the covariance operator, while in Section
LABEL:sec:reg_traj_sm, we detail the connection between local reglarity as a
probabilistic concept and the regularity of the sample paths. In Section
LABEL:sec:lemma_sm, we describe the embedding structure of the spaces
$\mathcal{X}(\cdot,\cdot;\mathcal{O}_{*})$. In Section LABEL:sec:details_sm,
we provide details on some quantities and equations from the main text. In
Section LABEL:sec:proof3_sm, we prove Theorem 3. Details on our simulation
experiments and additional results are gathered in Section
LABEL:sec:add_results_sm.

## References

  * Belloni et al., (2015) Belloni, A., Chernozhukov, V., Chetverikov, D., and Kato, K. (2015).  Some new asymptotic theory for least squares series: Pointwise and uniform results.  J. Econometrics, 186(2):345 – 366.
  * Blanke and Vial, (2014) Blanke, D. and Vial, C. (2014).  Global smoothness estimation of a Gaussian process from general sequence designs.  Electron. J. Stat., 8(1):1152–1187.
  * Cai and Yuan, (2010) Cai, T. and Yuan, M. (2010).  Nonparametric covariance function estimation for functional and longitudinal data.  University of Pennsylvania and Georgia Institute of Technology.
  * Cai and Yuan, (2011) Cai, T. T. and Yuan, M. (2011).  Optimal estimation of the mean function based on discretely sampled functional data: Phase transition.  Ann. Statist., 39(5):2330–2355.
  * Cai and Yuan, (2016) Cai, T. T. and Yuan, M. (2016).  Minimax and adaptive estimation of covariance operator for random variables observed on a lattice graph.  J. Amer. Statist. Assoc., 111(513):253–265.
  * Carroll et al., (2021) Carroll, C., Gajardo, A., Chen, Y., Dai, X., Fan, J., Hadjipantelis, P. Z., Han, K., Ji, H., Mueller, H.-G., and Wang, J.-L. (2021).  fdapace: Functional Data Analysis and Empirical Dynamics.  R package version 0.5.6.
  * Fan and Gijbels, (1996) Fan, J. and Gijbels, I. (1996).  Local polynomial modelling and its applications.  Monographs on statistics and applied probability. Chapman & Hall, London.
  * Gaïffas, (2007) Gaïffas, S. (2007).  Sharp estimation in sup norm with random design.  Statist. Probab. Lett., 77(8):782–794.
  * Goldenshluger and Lepski, (2011) Goldenshluger, A. and Lepski, O. (2011).  Bandwidth selection in kernel density estimation: Oracle inequalities and adaptive minimax optimality.  Ann. Statist., 39(3):1608–1632.
  * Golovkine et al., (2020) Golovkine, S., Klutchnikoff, N., and Patilea, V. (2020).  Learning the smoothness of noisy curves with application to online curve estimation. arxiv 2009.03652.
  * Guillou and Klutchnikoff, (2011) Guillou, A. and Klutchnikoff, N. (2011).  Minimax pointwise estimation of an anisotropic regression function with unknown density of the design.  Math. Methods Statist., 20(1):30–57.
  * Hall et al., (2006) Hall, P., Müller, H.-G., and Wang, J.-L. (2006).  Properties of principal component methods for functional and longitudinal data analysis.  Ann. Statist., 34(3):1493–1517.
  * Li and Hsing, (2010) Li, Y. and Hsing, T. (2010).  Uniform convergence rates for nonparametric regression and principal component analysis in functional/longitudinal data.  Ann. Statist., 38(6):3321–3351.
  * Ramsay and Silverman, (2005) Ramsay, J. and Silverman, B. W. (2005).  Functional Data Analysis.  Springer Series in Statistics. Springer-Verlag, New York, 2 edition.
  * Revuz and Yor, (2013) Revuz, D. and Yor, M. (2013).  Continuous Martingales and Brownian Motion.  Springer Science & Business Media.
  * Tsybakov, (2009) Tsybakov, A. B. (2009).  Introduction to Nonparametric Estimation.  Springer Series in Statistics. Springer New York.
  * Wong and Zhang, (2019) Wong, R. K. and Zhang, X. (2019).  Nonparametric operator-regularized covariance function estimation for functional data.  Comput. Statist. Data Anal., 131:131 – 144.
  * Zhang and Chen, (2007) Zhang, J.-T. and Chen, J. (2007).  Statistical inferences for functional data.  Ann. Statist., 35(3):1052–1079.
  * Zhang and Wang, (2016) Zhang, X. and Wang, J.-L. (2016).  From sparse to dense functional data and beyond.  Ann. Statist., 44(5):2281–2321.
  * Zhang and Wang, (2018) Zhang, X. and Wang, J.-L. (2018).  Optimal weighting schemes for longitudinal and functional data.  Statist. Probab. Lett., 138:165 – 170.