# Inference and Computation for Sparsely Sampled Random Surfaces

Tomas Masak  Tomas Rubin  Victor M. Panaretos  Ecole Polytechnique Fédérale de
Lausanne  Institut de Mathématiques
Ecole Polytechnique Fédérale de Lausanne
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Non-parametric inference for functional data over two-dimensional domains
entails additional computational and statistical challenges, compared to the
one-dimensional case. Separability of the covariance is commonly assumed to
address these issues in the densely observed regime. Instead, we consider the
sparse regime, where the latent surfaces are observed only at few irregular
locations with additive measurement error, and propose an estimator of
covariance based on local linear smoothers. Consequently, the assumption of
separability reduces the intrinsically four-dimensional smoothing problem into
several two-dimensional smoothers and allows the proposed estimator to retain
the classical minimax-optimal convergence rate for two-dimensional smoothers.
Even when separability fails to hold, imposing it can be still advantageous as
a form of regularization. A simulation study reveals a favorable bias-variance
trade-off and massive speed-ups achieved by our approach. Finally, the
proposed methodology is used for qualitative analysis of implied volatility
surfaces corresponding to call options, and for prediction of the latent
surfaces based on information from the entire data set, allowing for
uncertainty quantification. Our cross-validated out-of-sample quantitative
results show that the proposed methodology outperforms the common approach of
pre-smoothing every implied volatility surface separately.

62G05,

65P05,

Separability,

implied volatility surfaces,

functional data,

###### keywords:

[class=AMS]

###### keywords:

, and t1Research supported by a Swiss National Science Foundation grant.

###### Contents

  1. 1 Introduction
  2. 2 Methodology
    1. 2.1 Model and Observation Scheme
    2. 2.2 Motivation
    3. 2.3 Estimation of the Model Components
    4. 2.4 Prediction
  3. 3 Asymptotic Properties
  4. 4 Simulation Study
  5. 5 Data Analysis: Implied Volatility Surfaces
    1. 5.1 Quantitative Comparison
  6. 6 Conclusions and Future Directions
  7. A Computational Details
    1. A.1 Weighting Scheme
    2. A.2 Marginalization on a Grid
    3. A.3 Implementation Details
    4. A.4 Setup of the Simulation Study
    5. A.5 Additional Simulation Results
  8. B Proofs of Formal Statements
    1. B.1 Explicit Formula for Local Polynomial Regression
    2. B.2 Kernel Averages of $m$-dependent Data
    3. B.3 Proof of Proposition 1
    4. B.4 Proof of Theorem 1 and Corollary 1
    5. B.5 Proof of Proposition 2
    6. B.6 Proof of Theorem 2

##  1 Introduction

The term _random surfaces_ refers to continuous data on a two-dimensional
domain. Such data sets consist of multiple independent replications of some
underlying process $X=(X(t,s),t\in\mathcal{T},s\in\mathcal{S})$, forming a
random sample $X_{1},\ldots,X_{N}$, and can be tackled in the context of
functional data analysis (FDA), see [29, 30]. In practice, surface-valued data
cannot be observed as infinite-dimensional objects, one instead observes only
a finite number of measurements per surface. The classical split in literature
(see [41] for a comprehensive overview) consists of two measurement regimes:
the dense sampling and the sparse sampling. The former considers the setting
where the surface data are recorder densely enough (usually on a grid) so they
can be worked with as if they were truly infinite dimensional, possibly after
a pre-smoothing step, while retaining the same asymptotic properties as truly
infinite dimensional data [16]. The sparse regime, however, requires a
different approach. Here, the surfaces are observed only on a small (varying)
number of irregular locations, which also vary across the surfaces and are
corrupted with measurement errors. If such observations are gridded, which is
common for computational reasons, the resulting matrix-valued samples contain
many missing entries. Densely observed random surfaces arise frequently in
medical imaging, linguistics [28], or climate studies [15]. Examples of
sparsely observed random surfaces include longitudinal studies (where only a
part of a functional profile is measured at each visit) [24], geolocalized
data [38, 40, 36], or financial data [12, 21].

Regardless of the sampling regime, any functional data analysis will begin by
estimating two fundamental objects: the mean function
$\mu(t,s)=\mathbb{E}X(t,s)$ and the covariance kernel
$c(t,s,t^{\prime},s^{\prime})=\operatorname{Cov}(X(t,s),X(t^{\prime},s^{\prime}))$.
These ought to be estimated non-parametrically, since availability of
replicated observations should allow so. While the mean has the same dimension
as the underlying process itself, the covariance is an object of a higher
dimension, and its estimation thus suffers from the curse of dimensionality.
The size of a general covariance often poses both computational and
statistical issues, e.g. restricting the resolution of the grid one can handle
[1, 26].

Separability of the covariance is a popular non-parametric assumption to
reduce the statistical and computational burden, which arises when working on
a multi-dimensional domain. Well-known and often coupled together with
parametric models in the field of spatial statistics (see e.g. [8], and
references therein), separability has lately attracted a lot of attention also
from the non-parametric viewpoint of FDA [2, 5, 25, 10]. Intuitively, it
allows one to decouple the dimensions, regarding the random surfaces as an
enlarged ensemble of random curves instead. Specifically, separability assumes
that the full covariance can be written as a product of covariance kernels
corresponding to individual dimensions, i.e.:

$c(t,s,t^{\prime},s^{\prime})=a(t,t^{\prime})b(s,s^{\prime})$

for some bivariate kernels $a=a(t,t^{\prime})$ and $b=b(s,s^{\prime})$. Thus
separability decomposes the four-dimensional covariance kernel into two
kernels, each of the same dimensionality as the mean function.

Though separability cannot always be justified in the context of applied
problems [32], it is nevertheless frequently imposed because of the
computational advantages it entails [15, 13, 28]. And, while wrongfully
assuming separability introduces a bias, it may still lead to improved
estimation and prediction due to an implicit bias-variance trade-off. What is
certain is that separability greatly reduces computational costs of both the
estimation and prediction tasks.

To the best of our knowledge, the aforementioned properties of separability
have been mostly explored in the case of densely observed data. In the sparse
regime, there is no procedure for the non-parametric estimation of a separable
covariance. Methodology designed for the densely observed surfaces does not
apply, and methods designed for sparsely observed curves cannot feasibly be
extended to surfaces, as we will later argue. Since sparse data are generally
associated with higher computational complexity (extra costs associated with
smoothing) as well as higher statistical complexity (variance is inflated due
to sparse measurements burdened by noise), it appears that the gains stemming
from separability could be even larger in the sparse regime.

The goal of this paper is to leverage the separability assumption to reduce
complexity of covariance estimation down to that of mean estimation when
working under the sparse regime. The method of choice for sparsely observed
functional data on one-dimensional domains is the PACE approach [37], which is
based on kernel regression smoothers. A naive generalization of PACE to a two-
dimensional domain would entail a computationally infeasible local linear
smoothing step in four-dimensional space. Instead, we demonstrate how to make
careful use of separability to collapse the four-dimensional smoother into
several two-dimensional surface smoothers. Consequently, the estimation of the
mean, covariance, and noise level all become of similar computational
complexity. Our asymptotic theory shows that the estimators match the minimax
optimal convergence rates for two dimensional smoothing problems. Moreover,
our simulation study demonstrates the enormous computational gains offered by
separability, the statistical gains when the ground-truth is separable, and a
favorable bias-variance trade-off when the truth is in fact not separable.

Finally, we illustrate the benefits of our approach on a qualitative analysis
of implied volatility surfaces corresponding to call options. Here, one
surface corresponds to a fixed asset (e.g. a stock), for which the right to
buy it for an agreed-upon _strike_ price (first dimension) at a future _time
to expiration_ (second dimension) is traded. The value of implied volatility
at given strike and time to expiration is derived directly from observed
market data (the option prices), by the well-known and commonly used Black-
Scholes formula [3, 27]. The implied volatilities are preferred over the
option prices, since they are dimensionless quantities, allow for a direct
comparison of different assets, and are well familiar to practitioners. Since
the options are traded only for a finite number of strikes and times to
expiration, which vary across different surfaces, the observed data consist of
a sparse ensemble. The interpolation of such an ensemble is a typical
objective in financial mathematics as the latent implied volatility surfaces
are of interest for other tasks such as prediction (forecasting) of option
prices [20]. The common practice is to interpolate or smooth the measurements
for every surface independently. For example, Cont and Da Fonseca [6] utilized
pre-smoothing by local polynomial regression, evaluating an individual
smoother for every single surface. This approach may, however, pose issues
when the available sparse measurements for a given surface are concentrated
only on a subset of the domain, which is often the case. Once the predicted
surfaces are fed into a subsequent predictive models, the naively extrapolated
parts of the surface are given the same weight as the more reliable
interpolated parts. Instead, we advocate for the idea of “borrowing strength”
for the purpose of predicting the latent surfaces via best linear unbiased
prediction using the information from the entire data set, which also allows
for uncertainty quantification [37]. Under separability, we only need to use
two-dimensional surface smoothing, which is the case of the pre-smoothing
approach as well. At the same time, the proposed methodology outperforms the
pre-smoothing approach in terms of prediction error.

##  2 Methodology

###  2.1 Model and Observation Scheme

We assume the existence of i.i.d. latent surfaces
$X_{n}\in\mathcal{L}^{2}([0,1]^{2})$, $n=1,\ldots,N$, which are mean-square
continuous with continuous sample paths. We denote the mean function as
$\mu=\mu(t,s)$, where

$\mu(t,s)=\mathbb{E}X_{1}(t,s),\qquad t,s\in[0,1],$

and the covariance kernel as $c=c(t,s,t^{\prime},s^{\prime})$, where

$c(t,s,t^{\prime},s^{\prime})=\mathbb{E}\left[\left(X_{1}(t,s)-\mu(t,s)\right)\left(X_{1}(t^{\prime},s^{\prime})-\mu(t^{\prime},s^{\prime})\right)\right],\qquad
t,s,t^{\prime},s^{\prime}\in[0,1].$

We think of the first dimension as being temporal, denoted by variable $t$,
and the second dimension as being spatial, denoted by variable $s$, though
this convention is only made for the purposes of presentation.

The crucial assumption in this paper is that of separability of the
covariance:

  1. (A1)

The covariance kernel of the random surfaces $X_{1},\dots,X_{N}$ satisfies

$c(t,s,t^{\prime},s^{\prime})=a(t,t^{\prime})b(s,s^{\prime}),\qquad
t,s,t^{\prime},s^{\prime}\in[0,1],$ (2.1)

for some purely temporal covariance $a=a(t,t^{\prime})$ and some purely
spatial covariance $b=b(s,s^{\prime})$.

The process $X\in\mathcal{L}^{2}([0,1]^{2})$ is separable if it is, for
example, an outer product of two independent univariate processes (a purely
temporal one and a purely spatial one). In that case, apart from the
covariance, the mean function also separates into a product of a purely
temporal and a purely spatial functions. However, a process can have a
separable covariance even when it is not itself separable [32], for example
the mean function may not be separable. We do not assume separability of the
latent process itself, we only assume separability of its covariance in the
sense of (2.1).

We work under the sparse sampling regime, where every surface is observed only
at a finite number of irregularly distributed locations, and those
measurements are corrupted by independent additive errors. For the $n$-th
latent surface $X_{n}$, the number of measurements $M_{n}$ as well as the
locations of the measurements
$\\{(t_{nm},s_{nm})\;|\;m=1,\ldots,M_{n}\\}\subset[0,1]^{2}$ are considered
random, and the observations are given by the errors-in-measurements model
[37, 23, 41]:

$Y_{nm}=X_{n}(t_{nm},s_{nm})+\varepsilon_{nm},\quad
m=1,\ldots,M_{n},\;n=1,\ldots,N,$ (2.2)

where $\varepsilon_{nm}$ are i.i.d. with $\mathbb{E}\,\varepsilon_{nm}=0$ and
$\operatorname{Var}(\varepsilon_{nm})=\sigma^{2}>0$ being the noise level.

###  2.2 Motivation

In this section, we assume for simplicity that the mean is zero, and provide a
heuristic description of how one might estimate the separable covariance
(2.1). Note that we have

$\operatorname{Cov}(Y_{nm},Y_{nm^{\prime}})=a(t_{nm},t_{nm^{\prime}})b(s_{nm},s_{nm^{\prime}})+\sigma^{2}\mathds{1}_{[m=m^{\prime}]},\quad
n=1,\ldots,N,\;m,m^{\prime}=1,\ldots,M_{n}.$

Consider the “raw covariances” $G_{nmm^{\prime}}:=Y_{nm}Y_{nm^{\prime}}$.
Ignoring the assumption of separability, one could attempt to “lift” the PACE
approach [37] up to higher dimensions. This amounts to plotting the raw
covariances as a scatter plot in a four-dimensional domain, discarding the
diagonal covariances burdened by noise, and using a surface smoother to obtain
an estimator of the covariance. We refer to this procedure as 4D smoothing,
see Section A.3 for details. However, there are two issues with 4D smoothing.
The curse of dimensionality results in an estimator of poor quality, unless
surfaces are observed relatively densely and many replications are available.
And, especially when the latter is true, the computational costs of smoothing
in a higher dimension can be excessive. We make the assumption of separability
mainly to cope with these two issues, which is often the case in the
literature already when working with fully observed data [15, 28]. We do not
see separability as a critical modeling assumption, but rather as a
regularization, which possibly introduces a bias. Separability reduces both
the statistical and the computational complexity of the covariance estimation
task. This is always important when working with random surfaces, whatever
their mode of observation, but becomes particularly crucial when working with
sparsely observed surfaces.

In the following, we provide a heuristic on how separability can be used to
our advantage in the sparse observation regime. Assuming zero mean for now, we
have

$\mathbb{E}G_{nmm^{\prime}}=\mathbb{E}Y_{nm}Y_{nm^{\prime}}=a(t_{nm},t_{nm^{\prime}})b(s_{nm},s_{nm^{\prime}})+\sigma^{2}\mathds{1}_{[m=m^{\prime}]}.$

Imagine for a moment that the spatial kernel $b=b(s,s^{\prime})$ is known, and
consider the set of values

$\left\\{\frac{Y_{nm}Y_{nm^{\prime}}}{b(s_{nm},s_{nm^{\prime}})}\;\Bigg{|}\;m,m^{\prime}=1,\ldots,M_{n},m\neq
m^{\prime},n=1,\ldots,N\right\\}.$ (2.3)

The expectation of every point in this set is $a(t_{nm},t_{nm^{\prime}})$, so
we can chart these points in a scatter plot as

$\left(t_{nm},t_{nm^{\prime}},\frac{Y_{nm}Y_{nm^{\prime}}}{b(s_{nm},s_{nm^{\prime}})}\right),$

and use a two-dimensional surface smoother to obtain an estimator of
$a=a(t,t^{\prime})$.

Correspondingly, if one knew the temporal kernel $a=a(t,t^{\prime})$ instead,
the set of values (2.3) (with $a$ in the denominator instead of $b$) could be
arranged against $s_{nm}$ and $s_{nm^{\prime}}$ to obtain an estimator of
$b=b(s,s^{\prime})$. When neither the temporal kernel $a$ nor the spatial $b$
is known, one can start with a fixed $b$ and iterate between updates of $a$
and $b$, smoothing a scatterplot once per every single update.

However, there are two issues with such an approach. Firstly, for small
denominators, the corresponding points on the scatterplot are not reliable,
and using them as they are can have a severe negative impact on estimation
quality. Secondly, unless very few observations per surface are available, the
procedure above is computationally very demanding, see Section A.5.

To cope with these issues, we use weights for the surface smoother and utilize
gridding, i.e. split the domain into disjoint intervals and work on a grid.
While, gridding can significantly reduce computations already on a univariate
domain [37], the gains are much bigger in higher dimensions. In the following
section, we introduce our methodology in full from the theoretical
perspective, while computational aspects are deferred to Section A.3.

###  2.3 Estimation of the Model Components

We use local linear regression surface smoothers [11] to formalise the
heuristic described in the previous section and estimate the components of the
model from Section 2.1, i.e. the mean $\mu=\mu(t,s)$, the temporal kernel
$a=a(t,t^{\prime})$, the spatial kernel $b=b(s,s^{\prime})$, and the noise
level $\sigma^{2}$.

By applying a surface smoother to the set
$\\{(x_{k},y_{k},z_{k})\;|\;k=1,\ldots,M\\}\subset\mathbb{R}^{3}$ with given
weights $\\{w_{k}\;|\;k=1,\ldots,M\\}$, we understand calculating
$\widehat{\gamma}_{0}=\widehat{\gamma}_{0}(x,y)$ as the minimizer of the
weighted sum of squares

$(\widehat{\gamma_{0}},\widehat{\gamma_{1}},\widehat{\gamma_{2}})=\operatornamewithlimits{arg\,min}_{\gamma_{0},\gamma_{1},\gamma_{2}}\sum_{k=1}^{M}K\left(\frac{x-x_{k}}{h_{1}}\right)K\left(\frac{y-y_{k}}{h_{2}}\right)w_{k}\Big{[}z_{k}-\gamma_{0}-\gamma_{1}(x-x_{k})-\gamma_{2}(y-y_{k})\Big{]}^{2}$
(2.4)

for every fixed $(x,y)\in[0,1]^{2}$, where $K(\cdot)$ is a smoothing kernel
function and $h_{1},h_{2}>0$ are bandwidths. Throughout the paper, we use the
Epanechnikov kernel, utilize cross-validation to select the bandwidths, and
mention weights only when they are not all equal.

First, we estimate the mean by applying the surface smoother to the set

$\\{(t_{nm},s_{nm},Y_{nm})\;|\;m=1,\ldots,M_{n},\;n=1,\ldots,N\\}.$ (2.5)

Denote the resulting estimator by $\widehat{\mu}=\widehat{\mu}(t,s)$.

Next, consider the “raw” covariances
$G_{nmm^{\prime}}=\big{[}Y_{nm}-\widehat{\mu}(t_{nm},s_{nm})\big{]}\big{[}Y_{nm^{\prime}}-\widehat{\mu}(t_{nm^{\prime}},s_{nm^{\prime}})\big{]}.$
We begin by applying the surface smoother to the set

$\\{(t_{nm},t_{nm},G_{nmm^{\prime}})\;|\;m,m^{\prime}=1,\ldots,M_{n},\,m\neq
m^{\prime},\;n=1\ldots,N\\}$

to obtain a preliminary estimator of $a=a(t,t^{\prime})$, denoted by
$\widehat{a}_{0}=\widehat{a}_{0}(t,t^{\prime})$. Then we use this preliminary
estimator to calculate a proxy of $b=b(s,s^{\prime})$, like we described in
the previous section. Namely, we apply the surface smoother to the set

$\left\\{\left(s_{nm},s_{nm^{\prime}},\frac{G_{nmm^{\prime}}}{\widehat{a}_{0}(t_{nm},t_{nm^{\prime}})}\right)\;\Bigg{|}\;m,m^{\prime}=1,\ldots,M_{n},\,m\neq
m^{\prime},\;n=1\ldots,N\right\\}$ (2.6)

using weights $\\{\widehat{a}_{0}^{2}(t_{nm},t_{nm^{\prime}})\\}$ to obtain
$\widehat{b}_{0}=\widehat{b}_{0}(s,s^{\prime})$. If the denominator in the set
above is ever zero, we remove the corresponding point from the set. Note that
since the weights are equal exactly to the denominators squared, it makes
sense even formally that such a point is never considered for the surface
smoother. As the next step, we refine our estimator of $a$ by applying the
surface smoother to the set

$\left\\{\left(t_{nm},t_{nm^{\prime}},\frac{G_{nmm^{\prime}}}{\widehat{b}_{0}(s_{nm},s_{nm^{\prime}})}\right)\;\Bigg{|}\;m,m^{\prime}=1,\ldots,M_{n},\,m\neq
m^{\prime},\;n=1\ldots,N\right\\}$ (2.7)

using weights $\\{\widehat{b}_{0}^{2}(s_{nm},s_{nm^{\prime}})\\}$ from which
we obtain $\widehat{a}=\widehat{a}(t,t^{\prime})$. Finally, we refine the
estimator of $b$ by applying the surface smoother to set (2.6) with
$\widehat{a}_{0}$ replaced by $\widehat{a}$, resulting in the estimator
$\widehat{b}=\widehat{b}(s,s^{\prime})$.

Finally, once both the mean and the separable covariance have been estimated,
it remains to estimate the noise level $\sigma^{2}$, which is of interest e.g.
for the purposes of prediction. We begin by applying the surface smoother to
the set

$\left\\{\left(t_{nm},s_{nm},G_{nmm}\right)\;|\;m=1,\ldots,M_{n},\;n=1\ldots,N\right\\}$

to obtain $\widehat{V}=\widehat{V}(t,s)$. Note that since
$\mathbb{E}G_{nmm}\approx a(t_{nm},t_{nm})b(s_{nm},s_{nm})+\sigma^{2}$, we can
estimate $\sigma^{2}$ by

$\widehat{\sigma}^{2}=4\int_{1/4}^{3/4}\int_{1/4}^{3/4}\big{[}\widehat{V}(t,s)-\widehat{a}(t,t)\widehat{b}(s,s)\big{]}dtds,$

where (similarly to [37]) we integrate only along the middle part of the
domain to mitigate boundary issues.

Figure 1: Workflow of the proposed estimation procedure. We estimate firstly
the mean from the data $\\{Y_{nm}\\}$, then the separable covariance (in
several steps) from the raw covariances $\\{G_{nmm^{\prime}}\\}$, and finally
the noise level. A surface smoother over a 2D domain is utilized in every step
(once per a single thin arrow).

The workflow of the estimation scheme described above is visualised in Figure
1. The main novelty of our approach lies in the part where the separable
covariance is estimated. Separability allows us to reduce dimensionality of
the problem. Hence only two-dimensional smoothing is required, while a
straightforward multi-dimensional generalization of e.g. the PACE approach
[37] not utilizing separability would require four-dimensional smoothing to
estimate the covariance.

The estimation of the separable terms can be viewed as an iterative procedure,
where either $a$ or $b$ is kept fixed while the other term is being updated.
The initializing estimate $\widehat{a}_{0}$ is also obtained in this way,
starting from $\widehat{b}_{0}\equiv 1$. A natural question is whether one
should iterate this process until convergence, or simply stop after a single
step, and use e.g. $\widehat{a}_{0}$ as the estimator of $a$. The approach we
advocate for uses exactly two steps (for both $a$ and $b$). The reason is the
following. As will be shown in Section 3, the asymptotic distribution of
$\widehat{b}$ does not depend on $\widehat{a}_{0}$, and similarly for
$\widehat{a}$. This fact follows from separability. However, given the
motivation in the previous section, one can anticipate the finite sample
performance of $\widehat{a}$ to be better than that of $\widehat{a}_{0}$,
which is verified in our simulation study. One can think of the first step as
estimating the optimal weights consistently, and the second step as using
those consistently estimated weights to produce the estimators.

One of the distinctive features of our methodology is the use of this explicit
weighing scheme, where the weights for each of the two covariance kernels
depend on the other covariance kernel. We explain why this is done in the
following example. A more precise justification for the specific (quadratic)
form of the weights, which requires some additional background, is deferred to
Section A.1.

######  Example 1.

Assume we observe zero-mean surfaces sparsely, and one of these surfaces, say
$X_{1}$, is observed at four locations only, as depicted in Figure 2. First,
let us describe the 4D smoothing estimator at a fixed location
$(t,s,t^{\prime},s^{\prime})$, when the bandwidths are also fixed. $X_{1}$
contributes to $\widehat{c}(t,s,t^{\prime},s^{\prime})$ only if there is a
pair of two locations, where $X_{1}$ is observed, such that one of the
locations is close to $(t,s)$, and the other is close to
$(t^{\prime},s^{\prime})$. In this case, closeness in time, resp. space, is
measured by $h_{t}$, resp. $h_{s}$, which control the bandwidth of the
smoothing kernel. In Figure 2 (left), only a single pair of locations
contributes to estimation at $(t,s,t^{\prime},s^{\prime})$.

Now, let us contrast this to a single step in the proposed estimating
procedure. Assume that $b=b(s,s^{\prime})$ is fixed in the current step, and
we are estimating $a=a(t,t^{\prime})$. In this step, the spatial dimensions
are not explicitly considered, we are performing smoothing only in the
temporal dimensions. Hence the product of any two locations, where $X_{1}$ is
observed, contributes to $\widehat{a}$, as long as the locations are close to
$(t,s)$ and $(t^{\prime},s^{\prime})$ in the temporal domain. In the situation
displayed in Figure 2 (right), this leads to four contributing raw covariance
pairs. In other words, when estimating the temporal part of the covariance at
$(t,s,t^{\prime},s^{\prime})$, we can consider even raw covariances, _which
are spatially far_ from $(t,s,t^{\prime},s^{\prime})$. This is meaningful due
to separability. The adopted weighting scheme then ensures that raw
covariances arising from points which are spatially distant are appropriately
weighted.

To sum up, 4D smoothing can be understood as averaging over information about
$c(t,s,t^{\prime},s^{\prime})$ captured in raw covariances, whose locations
are close to $(t,s,t^{\prime},s^{\prime})$. Under separability, however, the
proposed methodology borrows information in a different manner, always
allowing for more freedom in one dimension or the other, depending on which
dimension is currently held fixed.

|  |  |
---|---|---|---
Figure 2: A single random surface is observed at four locations (depicted by
“$\mathsf{x}$”), and these observations contribute differently to estimation
of the covariance at a fixed location $(t,s,t^{\prime},s^{\prime})$ in the
case of 4D smoothing _( left)_ and one step of the proposed approach leading
to the estimator of $a(t,t^{\prime})$ _( right)_. The gray areas depict active
neighborhoods and the dashed lines depict the contributing raw covariances
(products of the values at the connected locations).

###  2.4 Prediction

Another objective of our methodology is the recovery of the latent surfaces
based on the sparse and noisy observations thereon. The prediction method we
are going to present in this section follows the principle of borrowing
strength across the entire data set, an expression framed by Yao, Müller and
Wang [37]. Specifically, consider the training data set composed of the sparse
observations made on random surfaces $X_{1},\dots,X_{N}$, i.e.
$\\{Y_{nm}:m=1,\dots,M_{n},\,n=1,\dots,N\\}$, and a new random surface
$X^{new}$ observed under the same sampling protocol as (2.2):

$Y^{new}_{m}=X^{new}(t^{new}_{m},s^{new}_{m})+\varepsilon^{new}_{m},\qquad
m=1,\dots,M^{new}.$ (2.8)

We assume that $X^{new}$ comes from the same population as
$X_{1},\dots,X_{N}$. In fact, it may be set (together with its sparse
measurements) to one of the training surfaces $X_{1},\dots,X_{N}$ if the task
is to predict one of those.

Our prediction method is calibrated on all the observations, and this
information is used for the prediction of $X^{new}$. This contrasts to the
pre-smoothing step often used in the functional data literature [29, 30],
typically in the dense regime, where the prediction (usually some kind of a
smoother) of $X^{new}$ is based only on $\\{Y^{new}_{m}:m=1,\dots,M^{new}\\}$.

Let $X^{new}$ be a random surface with the mean $\mu(t,s),\,t,s\in[0,1],$ and
the separable covariance kernel $a(t,t^{\prime})b(s,s^{\prime})$,
$t,t^{\prime},s,s^{\prime}\in[0,1]$, observed through sparse measurements
(2.8). Then the best linear unbiased predictor of the latent surface $X^{new}$
given the sparsely observed data
$\mathbb{Y}^{new}=(Y_{1}^{new},\ldots,Y_{M^{new}}^{new})$, denoted as
$\Pi(X^{new}|\mathbb{Y}^{new})$, is given by (cf. [19]):

$\Pi(X^{new}(t,s)|\mathbb{Y}^{new})=\mu(t,s)+\operatorname{Cov}(X^{new}(t,s),\mathbb{Y}^{new})\left[\operatorname{Var}(\mathbb{Y}^{new})\right]^{-1}(\mathbb{Y}^{new}-\mathbb{E}\mathbb{Y}^{new}),\qquad
t,s\in[0,1],$ (2.9)

where

$\displaystyle\operatorname{Cov}(X^{new}(t,s),\mathbb{Y}^{new})$
$\displaystyle=\Big{(}a(t,t^{new}_{m})b(s,s^{new}_{m})\Big{)}_{m=1}^{M^{new}}\in\mathbb{R}^{M^{new}},\qquad
t,s\in[0,1],$ (2.10) $\displaystyle\operatorname{Var}(\mathbb{Y}^{new})$
$\displaystyle=\Big{(}a(t^{new}_{m},t^{new}_{m^{\prime}})b(s^{new}_{m},s^{new}_{m^{\prime}})+\sigma^{2}\mathds{1}_{[m=m^{\prime}]}\Big{)}_{m,m^{\prime}=1}^{M^{new}}\in\mathbb{R}^{M^{new}\times
M^{new}},$ (2.11) $\displaystyle\mathbb{E}\mathbb{Y}^{new}$
$\displaystyle=\Big{(}\mu(t^{new}_{m},s^{new}_{m})\Big{)}_{m=1}^{M^{new}}\in\mathbb{R}^{M^{new}}.$

Formula (2.9) contains the unknown mean surface $\mu(\cdot,\cdot)$, covariance
kernels $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ as well as the measurement error
variance $\sigma^{2}$. In reality, we estimate these quantities from the
training data set, say sparsely observed measurements on $X_{1},\dots,X_{N}$,
and plug-in the estimates $\widehat{\mu},\widehat{a},\widehat{b}$, and
$\widehat{\sigma}^{2}$ into (2.9). We shall denote such predictor as
$\widehat{\Pi}(X^{new}|\mathbb{Y}^{new})$. We show in the following section
(Theorem 2) that the predictor $\widehat{\Pi}(X^{new}|\mathbb{Y}^{new})$
converges to its theoretical counterpart $\Pi(X^{new}|\mathbb{Y}^{new})$ as
the number of training samples $N$ grows to infinity.

In the rest of this section, we turn our attention to the construction of
confidence bands under a Gaussian assumption.

  1. (A2)

The random surface $X^{new}$ is a Gaussian random element in
$\mathcal{L}^{2}([0,1]^{2})$ and the measurement error ensemble
$\\{\varepsilon_{m}^{new}\\}_{m=1}^{M^{new}}$ is a Gaussian random vector.

First note that under the assumption (A2), the best linear unbiased predictor
(2.9) actually corresponds to the conditional expectation
$\mathbb{E}[X^{new}(t,s)|\mathbb{Y}^{new}]$. Furthermore, we may calculate the
conditional covariance structure as

$\operatorname{Cov}\left(X^{new}(t,s),X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)\\\
=a(t,t^{\prime})b(s,s^{\prime})-\operatorname{Cov}(X^{new}(t,s),\mathbb{Y}^{new})\left[\operatorname{Var}(\mathbb{Y}^{new})\right]^{-1}\left[\operatorname{Cov}(X^{new}(t^{\prime},s^{\prime}),\mathbb{Y}^{new})\right]^{\top},\qquad
t,t^{\prime},s,s^{\prime}\in[0,1].$ (2.12)

Moreover, we denote
$\widehat{\operatorname{Cov}}\left(X^{new}(t,s),X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)$
the empirical counterpart to (2.12), where the unknown quantities $a$, $b$ and
$\sigma^{2}$ are replaced by their estimators.

Fixing $(t,s)\in[0,1]^{2}$ and $\alpha\in(0,1)$, the $(1-\alpha)$-confidence
interval for $X^{new}(t,s)$ is given by

$\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})\pm
u_{1-\alpha/2}\sqrt{\widehat{\operatorname{Cov}}\left(X^{new}(t,s),X^{new}(t,s)|\mathbb{Y}^{new}\right)}$
(2.13)

where $u_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard Gaussian
law. The point-wise confidence band is then constructed by connecting the
intervals (2.13) for all $(t,s)\in[0,1]^{2}$.

The construction of the simultaneous confidence band is more involved, and we
shall use the technique proposed by Degras [9]. Define the conditional
correlation kernel

$\widehat{\operatorname{Corr}}\left(X^{new}(t,s),X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)=\frac{\widehat{\operatorname{Cov}}\left(X^{new}(t,s),X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)}{\sqrt{\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)\widehat{\operatorname{Var}}\left(X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)}}$
(2.14)

if the division on the right-hand side makes sense and zero otherwise, and
where
$\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)=\widehat{\operatorname{Cov}}\left(X^{new}(t,s),X^{new}(t,s)|\mathbb{Y}^{new}\right)$.

Then, construct the simultaneous confidence band by connecting the intervals

$\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})\pm\widehat{z}_{1-\alpha}\sqrt{\widehat{\operatorname{Cov}}\left(X^{new}(t,s),X^{new}(t,s)|\mathbb{Y}^{new}\right)},$
(2.15)

where $\widehat{z}_{1-\alpha}$ is the $(1-\alpha)$-quantile of the law of

$\widehat{W}=\sup_{t,s\in[0,1]^{2}}\left|\widehat{Z}(t,s)\right|$ (2.16)

with $\widehat{Z}$ being a Gaussian random surface with
$\operatorname{Cov}(\widehat{Z}(t,s),\widehat{Z}(t^{\prime},s^{\prime}))=\widehat{\operatorname{Corr}}\left(X^{new}(t,s),X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)$
for $t,t^{\prime},s,s^{\prime}\in[0,1]$, where $\widehat{\operatorname{Corr}}$
is the from (2.14). Therefore, by the definition,
$\mathbb{P}(\sup_{(t,s)\in[0,1]^{2}}|\widehat{Z}(t,s)|\leq\widehat{z}_{1-\alpha})=1-\alpha$.
Numerical calculation of this quantile is explained by Degras [9], who also
concludes that $\widehat{z}_{1-\alpha}<u_{1-\alpha/2}$. Therefore the point-
wise confidence band is always enveloped by the simultaneous confidence band,
as expected.

The asymptotic coverage of the point-wise band (2.13) and the simultaneous
band (2.15) is verified in Theorem 2 in the following section.

##  3 Asymptotic Properties

In this section, we establish consistency and convergence rates of the
estimators $\widehat{\mu}=\widehat{\mu}(t,s)$,
$\widehat{a}=\widehat{a}(t,t^{\prime})$, and
$\widehat{b}=\widehat{b}(s,s^{\prime})$, as well as the measurement error
variance $\sigma^{2}$.

The following assumptions refine the sparse observation scheme introduced in
Section 2.1.

  1. (B1)

The counts of measurements per surface $M_{n}$ are independent identically
distributed random variables with the law $M_{n}\sim\mathcal{M}>0$ such that
$\mathbb{P}(\mathcal{M}>1)>0$ and $\mathcal{M}\leq M^{max}$ where
$M^{max}\in\mathbb{N}$ is a constant.

  2. (B2)

The measurement locations $(t_{nm},s_{nm}),\,m=1,\dots,M_{n},\,n=1,\dots,N$,
are independent identically distributed random variables generated from the
density $f_{(t,s)}(\cdot,\cdot)$ on $[0,1]^{2}$. The density
$f_{(t,s)}(\cdot,\cdot)$ is assumed to be twice continuously differentiable
and positive on $[0,1]^{2}$.

  3. (B3)

The counts $(M_{n})$, the locations $(t_{nm},s_{nm})$, and the latent surfaces
$(X_{n})$ are independent.

The following two assumptions are required for consistent estimation of the
mean surface.

  1. (B4)

The mean surface $\mu(\cdot,\cdot)$ is twice continuously differentiable on
$[0,1]^{2}$.

  2. (B5)

There exists $\rho>2$ such that the random surface $X_{1}$ and the measurement
error $\varepsilon_{11}$ satisfy

$\sup_{(t,s)\in[0,1]^{2}}\mathbb{E}\left[|X_{1}(t,s)|^{\rho}\right]<\infty,\qquad\mathbb{E}\left[|\varepsilon_{11}|^{\rho}\right]<\infty.$

  3. (B6)

The bandwidths $h_{\mu,1},h_{\mu,2}$ for the mean estimator satisfy $(\log
N)/(Nh_{\mu,1}h_{\mu,2})=o(1)$ and, furthermore, we assume that they decay
with the same rate: $h_{\mu,1}\asymp h$ and $h_{\mu,2}\asymp h$ as
$N\to\infty$. The statement $x_{n}\asymp x^{\prime}_{n}$ as $n\to\infty$ for
two sequences $\\{x_{n}\\}$ and $\\{x^{\prime}_{n}\\}$ is understood as
$\lim_{n\to\infty}x_{n}/x^{\prime}_{n}\in(0,\infty)$, i.e. $x_{n}$ and
$x^{\prime}_{n}$ differ asymptotically up to a multiplicative constant

The common decay rate assumption in (B6) is not really required for our
asymptotic theory; it is imposed to simplify the statements on the convergence
rates. The same argument applies to the other bandwidths below. Because all
smoothing in our methodology is restricted to two dimensions, the bandwidths
are indeed expected to decay with the same rate, and we may assume that they
differ asymptotically up to a multiplicative constant.

In order to estimate the covariance kernels $a(\cdot,\cdot)$ and
$b(\cdot,\cdot)$ we need the following assumptions:

  1. (B7)

The covariance kernels $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ are twice
continuously differentiable on $[0,1]^{2}$.

  2. (B8)

There exists $\rho^{\prime}>2$ such that the random surface $X_{1}$ and the
measurement error $\varepsilon_{11}$ satisfy

$\sup_{(t,t^{\prime},s,s^{\prime})\in[0,1]^{4}}\mathbb{E}\left[|X_{1}(t,s)X_{1}(t^{\prime},s^{\prime})|^{\rho^{\prime}}\right]<\infty,\qquad\mathbb{E}\left[|\varepsilon_{11}|^{2\rho^{\prime}}\right]<\infty.$

  3. (B9)

The bandwidths $h_{a}$ and $h_{b}$ used for smoothing the covariance kernel
$a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ satisfy $(\log N)/(Nh_{a}^{2})=o(1)$ and
$(\log N)/(Nh_{b}^{2})=o(1)$ as $N\to\infty$ and, for simplicity of the
convergence rates statements, we assume $h_{a}\asymp h$ and $h_{b}\asymp h$ as
$N\to\infty$ where $h$ is from assumption (B6).

  4. (B10)

The true value of the covariance kernel $b(\cdot,\cdot)$, of the separable
model (2.1) satisfies

$\Theta\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\int_{0}^{1}\int_{0}^{1}b(s,s^{\prime})f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}\neq
0$ (3.1)

where $f_{s}(s)=\int_{0}^{1}f_{(t,s)}(t,s)dt$ is the marginal density of the
random location $s_{11}$.

While the other assumptions are standard in the smoothing literature,
assumption (B10) might be surprising, especially considering it is not
symmetric between $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$. The reason behind
this asymmetry is that our estimation methodology starts by smoothing the raw
covariances $G_{nmm^{\prime}}$ against $(t_{nm},t_{nm^{\prime}})$ in order to
produce the preliminary estimator $\widehat{a}_{0}(\cdot,\cdot)$. The
condition (3.1) ensures that the estimator $\widehat{a}_{0}(\cdot,\cdot)$
converges to a nonzero quantity, see the constant $\Theta$ in Theorem 1. By
contrast, this issue is not present in the follow-up steps. Due to positive
semi-definitness of $b(\cdot,\cdot)$, the constant $\Theta$ can only be zero
if all eigenfunctions of $b(\cdot,\cdot)$ are orthogonal to the marginal
sampling density $f_{s}$. This cannot happen e.g. unless $B$ is exactly low-
rank, and with all the eigenfunctions changing signs. From the practical
perspective, the condition is merely a technicality.

The estimation of the noise level $\sigma^{2}$ furthermore requires the
following assumption.

  1. (B11)

The bandwidths $h_{V,1},h_{V,2}$ for the smoother $\widehat{V}(\cdot,\cdot)$
satisfy $(\log N)/(Nh_{V,1}h_{V,2})=o(1)$ and, for simplicity of the
convergence rates statements, we assume $h_{V,1}\asymp h$ and $h_{V,2}\asymp
h$ as $N\to\infty$ where $h$ is from assumption (B6).

The mean surface asymptotic theory is presented as the following proposition.
Note that the separability assumption (A1) is not required here.

######  Proposition 1.

Under the assumptions (B1) – (B6):

$\sup_{(t,s)\in[0,1]^{2}}\left|\widehat{\mu}(t,s)-\mu(t,s)\right|=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)\qquad\text{as}\quad N\to\infty.$

Our main asymptotic result, the consistency and the convergence rates for the
separable model components (2.1), is presented in the following theorem.

######  Theorem 1.

Under the assumptions (A1), (B1) – (B10):

$\displaystyle\sup_{(t,t^{\prime})\in[0,1]^{2}}\left|\widehat{a}(t,t^{\prime})-\Theta
a(t,t^{\prime})\right|$
$\displaystyle=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right),$ (3.2)
$\displaystyle\sup_{(s,s^{\prime})\in[0,1]^{2}}\left|\widehat{b}(s,s^{\prime})-\frac{1}{\Theta}b(s,s^{\prime})\right|$
$\displaystyle=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right),$ (3.3)

as $N\to\infty$, where $\Theta$ is defined in (3.1).

The separable decomposition (2.1) is not identifiable, because a constant can
multiply one component while dividing the other, i.e.
$a(t,t^{\prime})b(s,s^{\prime})=\left[\lambda
a(t,t^{\prime})\right]\left[(1/\lambda)b(s,s^{\prime})\right]$,
$t,t^{\prime},s,s^{\prime}\in[0,1]$, for any $\lambda\in(0,\infty)$. Therefore
we can only aim to recover the covariance kernels $a(\cdot,\cdot)$ and
$b(\cdot,\cdot)$ up to a multiplicative constant and its reciprocal,
respectively. The number $\Theta$ in statements (3.2) and (3.3) plays the role
of such a constant and depends on the initialization of the algorithm, in our
case on the fact that the first estimator $\widehat{a}_{0}$ smooths the raw
covariances $G_{nmm^{\prime}}$ without any weighting. Still, the product
$\widehat{a}(t,t^{\prime})\widehat{b}(s,s^{\prime})$,
$t,t^{\prime},s,s^{\prime}\in[0,1]$,, estimates consistently the covariance
structure $c(t,s,t^{\prime},s^{\prime})=a(t,t^{\prime})b(s,s^{\prime})$,
$t,t^{\prime},s,s^{\prime}\in[0,1]$,, which is summarised in the following
corollary.

######  Corollary 1.

Under the assumptions (A1), (B1) – (B10):

$\sup_{(t,s,t^{\prime},s^{\prime})\in[0,1]^{4}}\left|\widehat{a}(t,t^{\prime})\widehat{b}(s,s^{\prime})-a(t,t^{\prime})b(s,s^{\prime})\right|=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)$

as $N\to\infty$.

Finally, the asymptotic behaviour of the noise level $\sigma^{2}$ is given as
the following proposition.

######  Proposition 2.

Under the assumptions (A1), (B1) – (B11):

$\widehat{\sigma}^{2}=\sigma^{2}+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)\qquad\text{as}\quad N\to\infty.$

This completes the asymptotic theory for our estimsators, and we now turn to
prediction. The following theorem shows that the predictor
$\widehat{\Pi}(X^{new}|\mathbb{Y}^{new})$ defined in Section 2.4 converges –
as the sample size grows to infinity – to its oracle counterpart (2.9), which
assumes the knowledge of the true distribution of the data. Moreover, the
theorem also proves the asymptotic coverage of the point-wise and simultaneous
confidence bands (2.13) and (2.15).

######  Theorem 2.

Under the assumptions (A1), (B1) – (B11):

$\sup_{(t,s)\in[0,1]^{2}}\left|\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})-\Pi(X^{new}(t,s)|\mathbb{Y}^{new})\right|=o_{\mathbb{P}}\left(1\right),\qquad\text{as}\quad
N\to\infty,$ (3.4)

conditionally on $\mathbb{Y}^{new}$.

Assuming further (A2) and fixing $\alpha\in(0,1)$:

$\displaystyle\forall(t,s)\in[0,1]^{2}\lim_{N\to\infty}\mathbb{P}\left(\left|\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)\right|\leq
u_{1-\alpha}\sqrt{\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}\,\middle|\,\mathbb{Y}^{new}\right)$
$\displaystyle=1-\alpha,$
$\displaystyle\lim_{N\to\infty}\mathbb{P}\left(\sup_{(t,s)\in[0,1]^{2}}\left[\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)\right]^{-1/2}\left|\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)\right|\leq\widehat{z}_{1-\alpha}\,\middle|\,\mathbb{Y}^{new}\right)$
$\displaystyle=1-\alpha.$

The rates established in this section manifest the statistical consequences of
separability. Corollary 1 shows that the complete covariance structure is
estimated with the rate $O_{\mathbb{P}}(\sqrt{(\log N)/(Nh^{2})}+h^{2})$,
which is the known optimal minimax convergence rate [11] for two dimensional
non-parametric regression. By steps similar to our proofs presented in Section
B, it could be shown that the empirical covariance smoother yields the
convergence rate $O_{\mathbb{P}}(\sqrt{(\log N)/(Nh^{4})}+h^{2})$. The
empirical covariance smoother’s convergence rate is thus slower than the one
found in Corollary 1, achieved via the separable model.

##  4 Simulation Study

We explore the finite sample performance of the proposed methodology by means
of a moderate simulation study (total runtime of about one thousand CPU
hours). Computational efficiency (relatively small runtimes) is achieved by
working on a $20\times 20$ grid, like described in Section A.2. Every surface
$X_{1},\ldots,X_{100}$ is first sampled fully on this grid as a zero-mean
matrix-variate Gaussian with covariance $C$ (to be specified), superposed with
noise (zero-mean i.i.d. Gaussian entries with variance $\sigma^{2}$), and then
sub-sampled in a way that only a fraction of the entries, selected at random,
is retained. The covariance $C$ is always standardized to have trace one, and
$\sigma^{2}$ is chosen such that the gridded white noise process is also trace
one.

Methods compared. We compare the proposed separable estimator
$\widehat{c}=\widehat{a}\cdot\widehat{b}$ against the non-separable empirical
estimator obtained by local linear smoothing in four dimensions (4D
smoothing), and also against the best separable approximation (BSA) [13, 26]
obtained from the fully observed and noise-free surfaces. We also compare the
proposed estimator against its one-step version
($\widehat{c}=\widehat{a}_{0}\cdot\widehat{b}_{0}$, cf. Figure 1), while it is
shown in Section A.5 that more than two steps lead to a similar performance as
the proposed two-step estimator.

Covariance choices. We consider four specific choices for the covariance: (a)
a separable covariance with Fourier basis eigenfunctions and power decay of
the eigenvalues; (b) Brownian sheet covariance, which is also separable; (c)
the parametric covariance introduced by Gneiting [14], which is non-separable;
and (d) a superposition of the Fourier covariance from (a) with another
separable covariance, having shifted Legendre polynomial eigenfunctions and
power decay of the eigenvalues. We discuss these choices in detail in Section
A.4. However, for the purpose of presenting the simulation results, it is only
important to point out the following. Firstly, while the Fourier setup (a) and
the Brownian setup (b) are separable, the Gneiting setup (c) and the Fourier-
Legendre setup (d) are non-separable. Secondly, while the Brownian setup (b)
and Gneiting setup (c) lead to rather flat covariances, the Fourier setup (a)
and the Fourier-Legendre setup (d) lead to quite wiggly (though infinitely
smooth) covariances.

Sparsity. In all simulation setups, we consider different percentages of the
entries observed $p\in\\{2,5,10,20,40,70\\}$. Since the grid size is $20\times
20$, this means for example that for $p=2$ we have $2/100\cdot 20^{2}=8$
observations per surface. For all the setups and percentages, we report
relative estimation errors $\|\widehat{C}-C\|_{2}/\|C\|_{2}$, where
$\widehat{C}$ is an estimator obtained by one of the four methods above (4D
smoothing, one-step, proposed, or BSA). The results are shown in Figure 3.
Every reported error was calculated as an average of 100 Monte Carlo runs.

(a) Fourier – separable, wiggly |  | (b) Brownian – separable, flat
---|---|---
|  |
(c) Gneiting – non-separable, flat |  | (d) Fourier-Legendre – non-separable, wiggly
|  |
Figure 3: Relative estimation errors depending on percentages of the surfaces
observed $p$ for 4 ground truth covariance choices (a)-(d) and 4 methods
compared. BSA provides a baseline, having access to full surfaces and hence
not depending on $p$. For 4D smoothing, only results for small $p$ are
reported.

Error components. The reported estimation errors can be thought of having four
components: (i) asymptotic bias, which is zero if the true covariance is
separable, i.e. in cases (a) and (b); (ii) error due to finite number of
samples ($N=100$); (iii) error due to sparse observations (i.e. not observing
the full surfaces); and (iv) noise contribution. BSA errors are always free of
the latter two, providing a baseline. The effect of not observing the full
surfaces is displayed for different values of $p$. Finally, the noise
contamination prevents smoothing approaches to reach the performance of BSA
even for $p$ large. Although our methodology explicitly handles noise, the
finite sample performance is better with noise-free data, which only BSA has
access to.

Results. There is a number of comments to be made about the results in Figure
3:

  1. 1.

In the setups where the covariance is flat, i.e. (b) and (c), the one-step and
the proposed approaches work the same, and 4D smoothing also works relatively
well. These two setups are simple in a sense, because information can be
borrowed quite efficiently via smoothing, regardless of whether the truth is
separable or not. Still, the proposed approach utilizing separability does not
perform worse than 4D smoothing even in the non-separable case (c), having the
advantage of being much faster to obtain. For $p=10$ the proposed estimator
takes only a couple of seconds while 4D smoothing takes about 40 minutes even
at this relatively small size of data.

  2. 2.

When the true covariance is wiggly, the proposed methodology clearly
outperforms 4D smoothing, both for the separable truth (a) and the non-
separable truth (d). The reason is that smoothing procedures are not very
efficient in this case, and borrowing strength via separability is imperative.

  3. 3.

The reason why error curves for 4D smoothing are only calculated up to $p=10$
is the computational cost of smoothing in higher dimensions. In fact,
performing 4D smoothing when $p=10$ took more time than calculationg of all
the remaining results combined (cross-validations included). The runtimes are
reported in Section A.5.

######  Remark 1.

We used a standard cross-validation strategy to choose bandwidths for the
separable model. However, for 4D smoothing, cross validation is not feasible –
using it would increase the total runtime of our simulation study to over one
thousand CPU _days_. Therefore, we chose bandwidths for 4D smoothing based on
those cross-validated for a separable model. This approach is taken throughout
the paper. See Section A.3 for details and Section A.5 for an empirical
demonstration of the effectiveness of this choice.

##  5 Data Analysis: Implied Volatility Surfaces

A European call option is a contract granting its holder the right, but not
the obligation, to buy an underlying asset, for example a stock, for an
agreed-upon strike price at a defined expiration time. Finding a model and
deriving a pricing formula for the fair price of a European call option was a
milestone problem in quantitative finance and stochastic calculus. [3, 27]
solved this problem and under the so-called the Black-Scholes-Merton model
they shown that the fair price of the European call option on a non-divident
paying asset is given by the Black-Scholes formula [20]:

$C_{t}^{BS}(m,\tau,\sigma_{S})=S_{t}F_{N(0,1)}(d_{1})-\kappa
e^{-rt}F_{N(0,1)}(d_{2}),$ (5.1) $d_{1}=\frac{-\log
m+\tau(r+\sigma_{S}^{2}/2)}{\sigma_{S}\sqrt{t}},\qquad d_{2}=\frac{-\log
m+\tau(r-\sigma_{S}^{2}/2)}{\sigma_{S}\sqrt{t}},$

where $m=\kappa/S_{t}$ is the moneyness defined as the ratio of the strike
$\kappa$ and the current underlying asset price $S_{t}$ at the current time
$t$, $\tau=T-t$ denotes the time to expiration, $\sigma_{S}$ is the volatility
parameter in the Black-Scholes-Merton model, $r$ is the risk-free interest
rate, and $F_{N(0,1)}(\cdot)$ denotes the cumulative distribution function of
the standard normal distribution. The only unknown quantity among the inputs
in (5.1) is the volatility $\sigma_{S}$.

Besides calculating the fair price of an option given the (estimated/realised)
volatility, the Black-Scholes formula (5.1) can be used in reverse: having
observed the market price of the option, denoted as $C_{t}^{*}(m,\tau)$, find
the value of $\sigma_{t}^{BS}(m,\tau)$ that solves the equation

$C_{t}^{BS}(m,\tau,\sigma_{t}^{IV}(m,\tau))=C_{t}^{*}(m,\tau).$

It can be shown that such value $\sigma_{t}^{IV}(m,\tau)>0$, called the
implied volatility, exists uniquely for each triplet of $m>0$, $\tau>0$, and
$C_{t}^{*}(m,\tau)>0$. Now, if the market indeed followed the Black-Scholes-
Merton model and the investors were rational, the implied volatility
$\sigma_{t}^{IV}(m,\tau)$ for various $m$ and $\tau$ would be constant.
However, this is not true for real market data pointing to the shortcomings of
the Black-Scholes-Merton model. Despite these shortcomings, the Black-Scholes
formula (5.1) is widely used for transforming the observed option prices into
an ensemble of implied volatilities in a bijective manner. The advantage of
considering the implied volatilities as opposed to the market prices of the
options is that the implied volatility surfaces tend to be smoother and
comparable across assets. Thus we can take advantage of the functional data
analysis framework.

In contrast to the European options, an American call option grants the right
to buy the underlying asset anytime until the expiration time $T$. The pricing
of American options on possibly dividend paying stocks is more complicated
because the pricing involves the optimal stopping problem. In general, no
closed form solution exists and numerical algorithms are required [7].
Likewise, the observed market option prices can be transformed into implied
volatilities.

In this section we consider the options data offered by DeltaNeutral [39].
This free data set contains the end-of-day prices as well as the calculated
implied volatilities for options on U.S. Equities markets. The data covers the
period from January 2003 until April 2019 but limits each month to contain the
daily options data on only one symbol (a stock or an index). The currently
included symbol changes every month and the options on some of the symbols are
American while some are European. For each month we pick randomly only one
trading day with the data on the currently available symbol and discard the
other trading days. Therefore the sample we analyse contains 196 snapshots
with option prices and implied volatilities. We discard the non-liquid options
and consider the contracts with the log-moneyness $\log
m=\log(K/S_{t})\in[-0.5,0.5]$ and the time to expiration $\tau=T-t\in[14,365]$
(in days). Moreover, we take the logarithm of the implied volatilities to
transform them from the domain $(0,\infty)$ onto the real line. To reduce
computational costs we round the log-moneyness and the time to expiration to
fall on a common $50\times 50$ grid, c.f. Section A.

Figure 4 shows two observations in our samples. The snapshot of Qualcomm Inc
(QCOM) features a cummulation of observation at the short expiration. Here, it
happens five times that two raw observations fall in the same pixel on the
common $50\times 50$ grid. In these few cases we calculate the average of the
two observations in each pair. Due to smoothness, the option prices (and hence
the implied volatilities) attain very similar values and thus this rounding
and averaging does not change the conclusions of our analysis. We have
observed this fact by fitting the model on finer grids while the estimates
remained similar.

(a) Dell on 01/19/2006 | (b) Qualcomm Inc on 02/07/2018 | (c) mean surface
---|---|---
|  |
Figure 4: Two sample snapshots of the considered log implied volatility
surfaces corresponding to the call options on the stocks of Dell Technologies
Inc on 01/19/2006 _( a)_ and Qualcomm Inc on 02/07/2018 _( b)_, and the mean
surface of the implied volatility gained from pooling all the data together _(
c)_.

The estimated mean surface is displayed on the right-hand side of Figure 4.
The mean surface captures the typical feature of the implied volatility
surfaces: the volatility smile [20]. The implied volatility is typically
greater for the options with moneyness away from 1, while this aspect is more
significant for shorter times to expiration.

Figure 5 displays the estimates of the separable covariance components by our
methodology presented in Section 2.3. The moneyness component demonstrates the
highest marginal variability at the center of the covariance surface, meaning
that the log implied volatility oscillates the most for the options with the
log-moneyness around 0 (i.e. moneyness 1). The marginal variance is lower as
the log-moneyness departs from 0. The eigendecomposition plot of the moneyness
covariance kernel shows that the most of variability is explained by a nearly
constant function with a small bump at the log-moneyness 0. The second leading
eigenfunction adjusts the peak at the log-moneyness around 0 to a greater
extent than the first eigenfunction. The covariance kernel corresponding to
the time to expiration variable is smoother and demonstrates slightly higher
marginal variability at shorter expiration. This phenomenon is well known for
implied volatility [20]. The eigendecomposition of this covariance kernel
indicates that the log implied volatility variation is mostly driven by the
constant function while the second leading eigenfunction adjusts the slope of
the surface for varying time to expiration.

|
---|---
|
Figure 5: Top-left: The estimated covariance kernel
$\widehat{a}=\widehat{a}(t,t^{\prime})$ corresponding to the time to
expiration variable. Top-right: The three leading eigenfunctions of the
spectral decomposition of the covariance kernel
$\widehat{a}=\widehat{a}(t,t^{\prime})$. Bottom-left and botton-right: The
same as above but for the estimated covariance kernel
$\widehat{b}=\widehat{b}(s,s^{\prime})$ corresponding to the log-moneyness
variable.

Figure 6 demonstrates our prediction techniques presented in Section 2.4
together with the 95 % simultaneous confidence band. We recall that the
confidence band aims to capture the latent smooth random surface itself, while
our raw observations are modelled by adding an error term. Therefore, the raw
data are not guaranteed to be covered in the confidence band.

|
---|---
Figure 6: Two views on prediction based on the call options written on the
stock of Dell Technologies Inc on 01/19/2006. The circles depict the available
sparse observations, the ribbons depict the predicted latent surface by the
method of Section 2.4, where the covariance structure was assumed separable,
and finally the transparent gray surfaces depict the 95 % simultaneous
confidence band for the latent log implied volatility surface.

###  5.1 Quantitative Comparison

The prediction method outline in Section 2.4 requires as an input the pairwise
covariances regardless whether they have been estimated by the separable
estimator $\widehat{a}(t,t^{\prime})\widehat{b}(s,s^{\prime})$ or the 4D
smoother $\widehat{c}(t,s,t^{\prime},s^{\prime})$. As the benchmark for our
comparison, we choose the locally linear kernel smoother [11] applied
individually for each surface as such smoothers constitute a usual pre-
processing step [6]. We will refer to this predictor as pre-smoothing. In this
section, we demonstrate that the predictive performance is comparable for both
covariance estimator strategies (the separable and the 4D smoother), and that
both of these approaches are superior to pre-smoothing. Moreover, the
separable smoother is substantially faster than the 4D smoother.

We compare the prediction error by performing a 10-fold cross-validation,
where the covariance structure is fitted always on varying 90 % of the
surfaces, with the remaining 10 % used for out-of-sample prediction. In the
set that is held out for prediction, we select some of the sparse observations
and predict them based on the remaining observations on that surface. We use
the following hold-out patterns:

  1. (a)

Leave one chain out. Since the options are quotes always for a range of
strikes, they constitute features known as option chains (c.f. Figure 4) where
multiple option prices (or equivalently implied volatilities) are available
for a fixed time to expiration. For those surfaces that include at least two
such chains, we remove gradually each chain and predict it based on the other
chains. Therefore the number of prediction tasks performed on a single surface
is equal to the number of chains observed per that surface.

  2. (b)

Predict in-the-money. Predict implied volatilities for below-average moneyness
(i.e. moneyness $m\leq 1$) based on the out-of-the-money observations
(moneyness $m>1$).

  3. (c)

Predict out-of-the-money. Predict implied volatilities for above-average
moneyness (i.e. moneyness $m\geq 1$) based on the in-the-money observations
(moneyness $m<1$).

  4. (d)

Predict short maturities. Predict the implied volatility for options with the
time to maturity $\tau<183$ [days] based on the implied volatility of the
options with the time to maturity $\tau\geq 183$ [days].

  5. (e)

Predict long maturities. Predict the implied volatility for options with the
time to maturity $\tau>183$ [days] based on the implied volatility of the
options with the time to maturity $\tau\leq 183$ [days].

All the prediction strategies are performed only for those surfaces where both
the discarded part and the kept part are non-empty. We measure the prediction
error on surface with the index $n$ (in the test partition within the $K$-fold
cross-validation) by the following root mean square error criterion, relative
to the pre-smoothing benchmark:

$RMSE^{\text{method}}(n)=\sqrt{\frac{\sum_{m\in
M_{n}^{\text{discarded}}}\left((\widehat{\Pi}^{\text{method}}(X(t_{nm},s_{nm})|\mathbb{Y}_{n}^{\text{kept}}))-Y_{nm}\right)^{2}}{\sum_{m\in
M_{n}^{\text{discarded}}}\left((\widehat{\Pi}^{\text{pre-
smooth}}(X(t_{nm},s_{nm})|\mathbb{Y}_{n}^{\text{kept}}))-Y_{nm}\right)^{2}}}$
(5.2)

where $Y_{nm},m=1,\dots,M_{n}$ are the implied volatility observations on the
$n$-th surface, $M_{n}^{\text{discarded}}\subset\\{1,\dots,M_{n}\\}$ denotes
the set of observations’ indexes discarded for the $n$-th surface,
$\mathbb{Y}_{n}^{\text{kept}}$ are the vectorized implied volatility
observations that were kept to be conditioned on. The predictor
$\Pi^{\text{pre-smooth}}$ denotes the pre-smoothing based on the observations
$\mathbb{Y}_{n}^{\text{kept}})$. The predictors $\Pi^{\text{method}}$ for
method being either the separable smoother or the 4D smoother constitute the
proposed predictors in this article with the covariance structure estimated by
either of the two smoothers. These predictors are always trained only on the
training partition (90 % of the surfaces) within the $10$-fold cross-
validation scheme. Note that this out-of-sample comparison is adversarial for
the proposed approach, because when predicting a fixed surface, the
measurements on that surface (and another 10% of measurements total) are not
used for the mean and covariance estimation. In practice, we naturally utilize
all available information. However, for the hold-out comparison study here,
that would require frequent re-fitting of the covariance, which would not be
computationally feasible, in particular for 4D smoothing.

Figure 7: Boxplots of root mean square errors relative to the pre-smoothing
benchmark (5.2) for prediction method of Section 2.4 with the covariance
estimated by 4D smoothing or the proposed separable approach, and different
hold-out patterns: (a) leave one chain out; (b) predict in-the-money; (c)
predict out-of-the-money; (d) predict short maturities; and (e) predict long-
maturities. Numbers inside the boxes provide numerical values of the median.
For a given method, RMSE value 3 means that the given method is 3-times worse
than the benchmark, while RMSE value 1/3 corresponds to 3-fold improvement.

Figure 7 presents the results under the five hold-out patterns in form of
boxplots created from the relative errors (5.2). We see that the prediction
errors based on estimated covariances, be it the separable smoother or the 4D
smoother, are typically smaller than the pre-smothing benchmark with the only
exception of the 4D smoother in the hold-out pattern (b). The predictive
performances of the separable and the 4D smoother are comparable, but they
differ a lot in terms of runtime. It typically takes 30 seconds to calculate
the separable smoother (including a cross-validation based selection of the
smoothing bandwidths), while the 4D smoother takes around 3 hours. The latter
runtime is moreover without considering any automatic selection of the
bandwidths, because such would be computationally infeasible. Hence we use the
bandwidths selected by the separable model, as described in Remark 1. The
calculations are performed on a quite coarse grid of size $20\times 20$. The
calculations on a dense grid, such as $50\times 50$ used in the qualitative
analysis in Section 5 are not feasible for the 4D smoother.

Therefore we conclude that the separable smoother approach enjoys a better
predictive performance than the pre-smoothing benchmark and – while having
having similar predictive performance as the predictor based on 4D smoothing –
is computationally much faster than the said competitor. In fact, it requires
two-dimensional smoothers only, just as the pre-smoothing benchmark.

##  6 Conclusions and Future Directions

In practice, covariances are often non-separable [15, 1, 2, 10, 28, 32], and
assuming separability induces a bias. The variance stemming from sparse
measurements and noise contamination is, however, often of a of larger
magnitude, thus sanctioning separability as a means to achieve a better bias-
variance trade-off. Moreover, separability entails faster computation and
lower storage requirements. As demonstrated above, these advantages are
pronounced in the sparse regime.

When data are observed fully and no smoothing is used (i.e. simple averages
are calculated at the grid points instead), calculation of the estimator takes
about 0.03 seconds in case of separability and 0.12 seconds in general, with
our grid size $20$ and sample size $N=100$. With larger grid size, the
acceleration gained by assuming separability naturally magnifies. But still,
separability is assumed mostly to save memory rather than time, in the fully
observed case [1]. On the other hand, with sparsely observed data and kernel
smoothing deployed, the speed-up factor stemming from separability can be of
the order of hundreds already with rather small grid sizes. The point we
advocate is that – with sparse observations, kernel smoothing deployed, and
compared to observations on a grid – the computational savings offered by
separability are the same in terms of memory, but are much more profound in
terms of runtime. This is due to the extra costs associated with smoothing.

Our analysis of implied volatility surfaces provides some insights into the
statistical dependencies of such sparsely observed data, and our quantitative
comparison demonstrated the prediction performance of our prediction method
over the pre-smoothing benchmark. We have formed our data set collecting
options across various symbols and timestamps with the reasoning that the
volatility across these come also from one population, which is debatable.
This study should be seen rather as a proof-of-concept on how to “borrow
strength” across the data set offered by DeltaNeutral [39], and how this
approach can be used to decrease the prediction error.

Having said that, we believe that our methodology can provide even better
results when considering a more homogeneous population such as the time series
of the implied volatility surfaces related to a single fixed symbol/asset. In
this case, pre-smoothing is usually required [6, 21] for forecasting of such
time series. Our methodology could avoid the pre-smoothing step in this case
by predicting the surfaces while borrowing the information across the entire
data set. Furthermore, our methodology could be easily tailored to predicting
principal components scores by conditional expectation (similarly to [37]),
which could be used for forecasting by a vector autoregression.

Combining separable covariances and time series of sparsely observed surfaces
hints at other directions of future work. For example, Rubín and Panaretos
[33] showed how to estimate the spectral density operator non-parametrically
from sparsely observed functional time series. Estimating separable spectral
density operator for sparsely observed surface-valued time series seems to be
within reach and is likely to provide predictions that benefit from the
information across time, and thus reducing the prediction error even more.

##  Appendix A Computational Details

###  A.1 Weighting Scheme

In this section, we provide a heuristic justification for the quadratic choice
of weights in the smoothers for the estimation of the covariance kernels
$a(\cdot,\cdot)$ and $b(\cdot,\cdot)$, such as (2.6) or (2.7). Then, we
provide a more precise justification, showing that the quadratic choice of
weights corresponds to the optimal choice, when data are observed densely.

We begin with the heuristic justification. The quadratic weights can be
motivated by the connection to weighted least squares. We recall that weighted
least squares are used for linear regression models where the model errors are
not necessarily i.i.d. Their covariance matrix is assumed to be a diagonal
matrix known up to a multiplicative constant:

$\mathbf{y}=\mathbf{X}\beta+\varepsilon,\qquad\mathbb{E}[\varepsilon|\mathbf{X}]=0,\qquad\operatorname{Var}(\varepsilon|\mathbf{X})=\sigma_{\varepsilon}^{2}\operatorname{diag}(\mathbf{v}),\qquad\sigma_{\varepsilon}^{2}>0,$
(A.1)

where $\mathbf{y}=(y_{1},\dots,y_{I})^{\top}$ is the response,
$\mathbf{X}=(\mathbf{x}_{1}^{\top},\dots,\mathbf{x}_{I}^{\top})^{\top}$ is the
model matrix, and $\operatorname{diag}(\mathbf{v})$ denotes the diagonal
matrix with the known vector $\mathbf{v}=(v_{1},\dots,v_{I})^{\top}$ on its
diagonal. The regression coefficients $\beta$ in the model (A.1) are estimated
by the weighted least squares:

$\widehat{\beta}=\operatornamewithlimits{arg\,min}_{\beta}\sum_{i=1}^{I}w_{i}\left(y_{i}-\mathbf{x}_{i}\beta\right)^{2},\qquad\text{where}\quad
w_{i}=\frac{1}{v_{i}},\,\,i=1,\dots,I.$ (A.2)

Consider the surface smoother of

$\left\\{\left(s_{nm},s_{nm^{\prime}},\frac{G_{nmm^{\prime}}}{\alpha(t_{nm},t_{nm^{\prime}})}\right)\;\Bigg{|}\;m,m^{\prime}=1,\ldots,M_{n},\,m\neq
m^{\prime},\;n=1\ldots,N\right\\}$

where $\alpha(t,t^{\prime}),\,t,t^{\prime}\in[0,1]$, is a fixed deterministic
kernel. The kernel smoothing technique we deploy is based on fitting a linear
regression locally. In view of model (A.1) we want to assess the variance of
the response $G_{nmm^{\prime}}/\alpha(t_{nm},t_{nm^{\prime}})$ to improve the
estimation procedure:

$\operatorname{Var}\left(\frac{G_{nmm^{\prime}}}{\alpha(t_{nm},t_{nm^{\prime}})}\right)=\frac{1}{\alpha^{2}(t_{nm},t_{nm^{\prime}})}\operatorname{Var}\left(G_{nmm^{\prime}}\right).$
(A.3)

The variance of $G_{nmm^{\prime}}$ is unknown and therefore cannot be used to
improve the estimation. Still, we observe in equation (A.3) that the variance
is multiplied by the reciprocal of $\alpha^{2}(t_{nm},t_{nm^{\prime}})$.
Therefore, we would define the weights for the weighted least squares (A.2) as
$w_{i}=\alpha^{2}(t_{nm},t_{nm^{\prime}})$, to utilize the knowledge we
actually have.

Let us now describe the connection between the quadratic weighting scheme and
the case of fully observed surfaces. With fully observed surfaces, the
separable model can be estimated via the generalized power iteration method,
where a single step is given by the partial inner product between the
empirical covariance and the previous step [26]. For example, when
$b=b(s,s^{\prime})$ is fixed, one step of the power iteration method is given
by [26, Proposition 1]

$\widehat{a}(t,t^{\prime})=\int_{0}^{1}\int_{0}^{1}b(s,s^{\prime})\widehat{c}_{N}(t,s,t^{\prime},s^{\prime})dsds^{\prime}\bigg{/}\int_{0}^{1}\int_{0}^{1}b^{2}(s,s^{\prime})dsds^{\prime},$
(A.4)

where $\widehat{c}_{N}$ is the empirical covariance estimator. In this
section, we demonstrate that, with fully observed data and with no smoothing
conducted, the estimation methodology of Section 2.3 corresponds to the power
iteration step (A.4).

Firstly, assume that $b$ is fixed, and we are using the surface smoother on
the set of points (2.7) to obtain $\widehat{a}$. Assume that the $n$-th
surface is observed twice at the temporal location $t$, i.e. at two locations
$(t,s_{1})$ and $(t,s_{2})$, and once more in a general location
$(t^{\prime},s^{\prime})$. Let us denote the raw covariance corresponding to
the $n$-th surface and locations $(t,s)$ and $(t^{\prime},s^{\prime})$
explicitly by $G_{n}(t,s,t^{\prime},s^{\prime})$. Then, two values are
available for the location $(t,t^{\prime})$ in set (2.7):

$\frac{G_{n}(t,s_{1},t^{\prime},s^{\prime})}{b(s_{1},s^{\prime})}\quad\&\quad\frac{G_{n}(t,s_{2},t^{\prime},s^{\prime})}{b(s_{2},s^{\prime})}.$

The corresponding weights are $b^{2}(s_{1},s^{\prime})$ and
$b^{2}(s_{2},s^{\prime})$, respectively. If the bandwidth is small enough, and
no other observations are available for this location,
$\widehat{a}(t,t^{\prime})$ is calculated as a weighted average:

$\left[b^{2}(s_{1},s^{\prime})\frac{G_{n}(t,s_{1},t^{\prime},s^{\prime})}{b(s_{1},s^{\prime})}+b^{2}(s_{2},s^{\prime})\frac{G_{n}(t,s_{2},t^{\prime},s^{\prime})}{b(s_{2},s^{\prime})}\right]\Big{/}\big{[}b^{2}(s_{1},s^{\prime})+b^{2}(s_{2},s^{\prime})\big{]}$

Also, for the purposes of the surface smoother, using the two points
separately with their separate quadratic weights is equivalent to using the
weighted average with the weight
$b^{2}(s_{1},s^{\prime})+b^{2}(s_{2},s^{\prime})$.

When the temporal slice $t$ of the $n$-th surface is observed fully, the
weighted averaging can be done continuously:

$\widehat{a}(t,t^{\prime})=\int_{0}^{1}b(s,s^{\prime})G_{n}(t,s,t^{\prime},s^{\prime})ds\bigg{/}\int_{0}^{1}b^{2}(s,s^{\prime})ds.$

When the temporal slice $t^{\prime}$ of the $n$-th surface is also observed
fully, the weighted average becomes

$\widehat{a}(t,t^{\prime})=\int_{0}^{1}\int_{0}^{1}b(s,s^{\prime})G_{n}(t,s,t^{\prime},s^{\prime})dsds^{\prime}\bigg{/}\int_{0}^{1}\int_{0}^{1}b^{2}(s,s^{\prime})dsds^{\prime}.$

When this is true for all $N$ surfaces, the result is averaged over all the
independent realizations, and we arrive directly to (A.4), since
$\widehat{c}_{N}(t,s,t^{\prime},s^{\prime})=\frac{1}{N}\sum_{n=1}^{N}G_{n}(t,s,t^{\prime},s^{\prime})$.

Altogether, our estimation procedure can be thought of (due to the specific
weighting scheme used) as a sparse version of the generalized power iteration
method of [26]. This link has important computational implications, which are
discussed in the following section.

###  A.2 Marginalization on a Grid

By marginalization, we mean preparation of the raw covariances for the 2D
smoothing step, i.e. charting the raw covariances either in time or in space
and weighting them as in formulas (2.7) and (2.6). In this section, we will
show that scatter points can be pooled together during the marginalization
process to save computations during the subsequent smoothing step, when the
data are observed (or rounded to) a common grid. The associated computational
advantages are discussed in the following section.

Assume for the remainder of this section that data from the measurement scheme
(2.2) arrive as matrices
$\mathbf{Y}_{1},\ldots,\mathbf{Y}_{N}\in\mathbb{R}^{d_{1}\times d_{2}}$ with
only some of their entries known, i.e. most of the entries are missing. The
marginal covariance kernels $a=a(t,t^{\prime})$ and $b=b(s,s^{\prime})$ are
replaced by matrices $\mathbf{A}\in\mathbb{R}^{d_{1}\times d_{1}}$ and
$\mathbf{B}\in\mathbb{R}^{d_{2}\times d_{2}}$, respectively. We assume again
for simplicity that the mean $\mu=\mu(t,s)$ is zero. The raw covariances then
form a tensor
$\mathbf{G}_{n}=\mathbf{Y}_{n}\otimes\mathbf{Y}_{n}\in\mathbb{R}^{d_{1}\times
d_{2}\times d_{1}\times d_{2}}$ with entries
$\mathbf{G}_{n}[i,j,i^{\prime},j^{\prime}]=\mathbf{Y}_{n}[i,j]\mathbf{Y}_{n}[i^{\prime},j^{\prime}]$.

Again, like in the previous section, assume that $\mathbf{B}$ is fixed, and we
are using the surface smoother on the discrete equivalent to set (2.7), i.e.

$\bigg{\\{}\left(i,i^{\prime},\frac{\mathbf{G}_{n}[i,j,i^{\prime},j^{\prime}]}{\mathbf{B}[j,j^{\prime}]}\right)\bigg{|}\;\mathbf{Y}_{n}\text{
observed at }(i,j)\text{ and
}(i^{\prime},j^{\prime}),\;(i,j)\neq(i^{\prime},j^{\prime}),\;n=1,\ldots,N\bigg{\\}}$
(A.5)

to obtain $\widehat{A}$. Like in the previous section, assume $\mathbf{Y}_{n}$
was observed at locations at $[i,j_{1}]$, $[i,j_{2}]$ and
$[i^{\prime},j^{\prime}]$ where no two locations are the same. As explained in
the previous section, it is equivalent for the surface smoother to replace the
corresponding two values from (A.5), i.e.

$\frac{\mathbf{G}_{n}[i,j_{1},i^{\prime},j^{\prime}]}{\mathbf{B}[j_{1},j^{\prime}]}\quad\&\quad\frac{\mathbf{G}_{n}[i,j_{2},i^{\prime},j^{\prime}]}{\mathbf{B}[j_{2},j^{\prime}]}$

with weights $\mathbf{B}[j_{1},j^{\prime}]$ and
$\mathbf{B}[j_{2},j^{\prime}]$, by a single value

$\big{(}\mathbf{B}[j_{1},j^{\prime}]\mathbf{G}_{n}[i,j_{1},i^{\prime},j^{\prime}]+\mathbf{B}[j_{2},j^{\prime}]\mathbf{G}_{n}[i,j_{2},i^{\prime},j^{\prime}]\big{)}\big{/}\big{(}\mathbf{B}^{2}[j_{1},j^{\prime}]+\mathbf{B}^{2}[j_{2},j^{\prime}]\big{)}$
(A.6)

with the aggregated weight
$\mathbf{B}^{2}[j_{1},j^{\prime}]+\mathbf{B}^{2}[j_{2},j^{\prime}]$.

Let $\mathbf{y}_{n,i}$ (resp. $\mathbf{y}_{n,i^{\prime}}$) denote the $i$-th
(resp. $i^{\prime}$-th) column of $\mathbf{Y}_{n}$. Let $\mathbf{q}_{i}$
denote the identifier of whether the entries of $\mathbf{y}_{n,i}$ were
observed (and similarly $\mathbf{q}_{i^{\prime}}$), i.e.

$\mathbf{q}_{i}[l]=\begin{cases}1,\;l\in\\{j_{1},j_{2}\\}\\\
0,\;\text{otherwise}\end{cases}\quad\&\quad\mathbf{q}_{i^{\prime}}[l]=\begin{cases}1,\;l=j^{\prime}\\\
0,\;\text{otherwise}.\end{cases}$

Then, value (A.6) can be calculated as
$\mathbf{y}_{n,i}^{\top}\mathbf{B}\mathbf{y}_{n,i}/\mathbf{q}_{i}^{\top}\mathbf{B}_{2}\mathbf{q}_{i^{\prime}}$
with the aggregated weight given by
$\mathbf{q}_{i}^{\top}\mathbf{B}_{2}\mathbf{q}_{i^{\prime}}$, where
$\mathbf{B}_{2}$ is the entry-wise square of $\mathbf{B}$. Naturally, this can
be generalized to the case when arbitrary number of entries in the $i$-th and
$i^{\prime}$-th columns of $\mathbf{Y}_{n}$ are observed. But more
importantly, it can be also generalized to account for different pairs of
columns of $\mathbf{Y}_{n}$ at the same time.

Let $\mathbf{Q}_{n}$. Then the contribution of the $n$-th surface
$\mathbf{Y}_{n}$ into set (A.5) can be calculated at once as

$\mathbf{Y}_{n}^{\top}\widetilde{\mathbf{B}}\mathbf{Y}_{n}/\mathbf{Q}^{\top}\widetilde{\mathbf{B}}_{2}\mathbf{Q}$
(A.7)

where $\widetilde{\mathbf{B}}$ is $\mathbf{B}$ with the diagonal values
replaces by zeros and $\widetilde{\mathbf{B}}_{2}$ is the entry-wise square of
$\widetilde{\mathbf{B}}$. The diagonal values of $\mathbf{B}$ are replaced by
zeros as described in order to discard products of the type
$\mathbf{Y}^{2}[i,j]$, which are burdened with noise.

The situation is analogous in the other step, when $\mathbf{A}$ is fixed and
$\mathbf{B}$ is calculated. The whole procedure of estimating the separable
covariance based on gridded sparse measurements is outlined in Algorithm 1.

Algorithm 1 Estimation of separable model from sparsely observed zero-mean
surfaces.

Input

$\mathbf{Y}_{1},\ldots,\mathbf{Y}_{N}\in(\mathbb{R}\cup\\{\diamond\\})^{d_{1}\times
d_{2}}$, where $\diamond$ represents a missing value

$\mathbf{Q}_{n}:=\mathds{1}_{[\mathbf{Y}_{n}\neq\diamond]}\in\\{0,1\\}^{d_{1}\times
d_{2}}$, for $n=1,\ldots,N$

replace all $\diamond$ entries in $\mathbf{Y}_{1},\ldots,\mathbf{Y}_{N}$ by
zeros

$\mathbf{A}=\big{(}1\big{)}_{i,j=1}^{d_{1}\times d_{1}}$

repeat

for

$n=1,\ldots,N$

$\widetilde{\mathbf{B}}:=\mathbf{B}$ with diagonal entries replaced by zeros

$\widetilde{\mathbf{B}}_{2}:=$ entry-wise square of $\widetilde{\mathbf{B}}$

$\mathbf{W}_{n}:=\mathbf{Q}_{n}\widetilde{\mathbf{B}}_{2}\mathbf{Q}_{n}^{\top}$

$\mathbf{Z}_{n}:=\mathbf{Y}_{n}\widetilde{\mathbf{B}}\mathbf{Y}_{n}^{\top}$
entry-wise divided by $\mathbf{W}_{n}$

end for

$\mathbf{A}:=$ surface smoother of
$\\{\mathbf{Z}_{1},\ldots,\mathbf{Z}_{N}\\}$ with
$\\{\mathbf{W}_{1},\ldots,\mathbf{W}_{N}\\}$ as the smoothing weights

for

$n=1,\ldots,N$

$\widetilde{\mathbf{A}}:=\mathbf{A}$ with diagonal entries replaced by zeros

$\widetilde{\mathbf{A}}_{2}:=$ entry-wise square of $\widetilde{\mathbf{A}}$

$\mathbf{W}_{n}:=\mathbf{Q}_{n}^{\top}\widetilde{\mathbf{A}}_{2}\mathbf{Q}_{n}$

$\mathbf{Z}_{n}:=\mathbf{Y}_{n}^{\top}\widetilde{\mathbf{A}}\mathbf{Y}_{n}$
entry-wise divided by $\mathbf{W}_{n}$

end for

$\mathbf{B}:=$ surface smoother of
$\\{\mathbf{Z}_{1},\ldots,\mathbf{Z}_{N}\\}$ with
$\\{\mathbf{W}_{1},\ldots,\mathbf{W}_{N}\\}$ as the smoothing weights

until convergence (or only twice)

Output

$\mathbf{A},\mathbf{B}$

###  A.3 Implementation Details

Separability offers reductions in both time and memory complexities already
when data are observed fully [1, 26]. In this section, we argue that
computational gains of separability are even greater, when data are observed
sparsely and kernel smoothing is used.

Kernel smoothers are known to be computationally demanding. To directly
evaluate a kernel smoother in $d_{1}$ locations using $d_{2}$ observations
takes $\mathcal{O}(d_{1}d_{2})$ operations. Table 1 shows these quadratic
complexities in our situation, explained below. Assume that $N$ surfaces were
observed on a grid of size $d\times d$ relatively densely (i.e. a fixed
percentage of the grid was observed – this is not unrealistic since one often
chooses the grid size in such a way), and an unbounded kernel was used. The
quadratic complexity of kernel smoother translates into estimating a general
covariance by a surface smoother in $\mathcal{O}(Nd^{8})$ operations, because
all $\mathcal{O}(Nd^{4})$ raw covariances have to be accessed at every single
one of $\mathcal{O}(d^{4})$ grid points. When we consider $N$ fixed, the
resulting complexity in $d$, i.e. $\mathcal{O}(d^{8})$, is huge. Under
separability, not using the marginalization procedure, the complexity is
$\mathcal{O}(d^{6})$, because $\mathcal{O}(d^{4})$ raw covariances have to be
accessed at $\mathcal{O}(d^{2})$ grid points. With marginalization, i.e. using
formula (A.7), the time complexity drops down to $\mathcal{O}(d^{4})$, because
the number of raw covariances that has to be accessed at every grid point
decreases to $\mathcal{O}(d^{2})$.

In practice, the quadratic complexity of kernel smoothers becomes intractable
and there exist many computational approaches to reduce the burden. Most
notably, the fast Fourier transform can be used on equispaced domains to
reduce quadratic complexity to log-linear [34], effectively cutting down the
powers of $d$ in the first row of Table 1. Many other accelerating approaches
exist, see e.g. [31, 22], and references therein. However, software
availability utilizing these computationally efficient approaches is rather
limited, and this is particularly true for multi-dimensional problems.

We do not provide our own implementation of kernel smoothing. For the proposed
approach, we implement Algorithm 1, which uses a “surface smoother”. To this
end, we utilize local linear smoothers provided in the fdapace package [37,
4]. We also utilize internal functions from fdapace to perform cross-
validation for the choice of bandwidths.

Table 1: Complexities for covariance estimation of a random surface observed on a $d\times d$ grid. Complexity | Separability | Separability w/o | 4D smoothing
---|---|---|---
|  | Marginalization |
Time | $\mathcal{O}(d^{4})$ | $\mathcal{O}(d^{6})$ | $\mathcal{O}(d^{8})$
Memory | $\mathcal{O}(d^{2})$ | $\mathcal{O}(d^{4})$ | $\mathcal{O}(d^{4})$

For 4D smoothing, which we consider only for comparison, we use the np package
[18], which is to the best of our knowledge the only R [35] package able to
perform local linear polynomial regression surface smoothing in more than two
dimensions. The 4D smoothing estimator requires smoothing in four dimensions.
We provide the formula here for completeness. For the set of points
$\\{(x_{k},y_{k},x_{k}^{\prime},y_{k}^{\prime})|k=1,\ldots,M\\}$ and weights
$\\{w_{k}|k=1,\ldots,M\\}$, the smoothed surface
$\widehat{\gamma}_{0}=\widehat{\gamma}_{0}(x,y,x^{\prime},y^{\prime})$ is
calculated as

$\begin{split}(\widehat{\gamma_{0}},\widehat{\gamma_{1}},\widehat{\gamma_{2}},\widehat{\gamma_{3}},\widehat{\gamma_{4}})=\operatornamewithlimits{arg\,min}_{\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}}\sum_{k=1}^{M}&K\left(\frac{x-x_{k}}{h_{1}}\right)K\left(\frac{y-y_{k}}{h_{2}}\right)K\left(\frac{x^{\prime}-x_{k}^{\prime}}{h_{3}}\right)K\left(\frac{y^{\prime}-y_{k}^{\prime}}{h_{4}}\right)w_{k}\cdot\\\
&\cdot\Big{[}z_{k}-\gamma_{0}-\gamma_{1}(x-x_{k})-\gamma_{2}(y-y_{k})-\gamma_{3}(x^{\prime}-x_{k}^{\prime})-\gamma_{4}(y^{\prime}-y_{k}^{\prime})\Big{]}^{2}\end{split}$
(A.8)

for every fixed $(x,y,x^{\prime},y^{\prime})\in[0,1]^{4}$, where $K(\cdot)$ is
a smoothing kernel function and $h_{1},h_{2},h_{3},h_{4}>0$ are bandwidths.
Even though the np package implements cross-validation to choose the
bandwidths, we found the computational burden to be huge and the performance
rather poor in our setups, see Section A.5. Hence, whenever we use 4D
smoothing, we fix the unknown bandwidths as chosen by cross-validation for the
proposed separable model. While this intuitively leads to smaller than optimal
bandwidths, we found out in our simulation study that bandwidths are governed
mainly by smoothness of the underlying covariance rather than the number of
points per surface available. Since the smoothness of a four-dimensional
covariance and its separable proxy is similar, optimal bandwidths chosen for
the proposed (separable) approach seem to be reasonable for 4D smoothing as
well, and this was verified in our experiments. Regardless, we can hardly
afford other strategy for choosing the four bandwidths for 4D smoothing. Even
the sophisticated combination of cross-validation and optimization provided in
the np package leads huge runtimes in our setups, see Section A.5.

###  A.4 Setup of the Simulation Study

Here, we fully specify the four covariances chosen for our simulation study
(Section 4).

(a)

In the Fourier scenario, $a=a(t,t^{\prime})$ and $b=b(s,s^{\prime})$ are
chosen to be the same, such that they have the trigonometric basis as their
eigenfunctions and power decay of their eigenvalues, resulting in a rather
wiggly univariate covariance displayed in Figure 8 (left). The covariance is
then set as $c(t,s,t^{\prime},s^{\prime})=a(t,t^{\prime})\,b(s,s^{\prime})$,
resulting in a separable covariance.

(b)

In the Brownian scenario, $a=a(t,t^{\prime})$ and $b=b(s,s^{\prime})$ are both
chosen as the covariance of the Wiener process, i.e.
$a(t,t^{\prime})=\min(t,t^{\prime})$ and $b(s,s^{\prime})=\min(s,s^{\prime})$,
resulting in a rather flat covariance displayed in Figure 8 (center). The
covariance is then set as
$c(t,s,t^{\prime},s^{\prime})=a(t,t^{\prime})\,b(s,s^{\prime})$, i.e. it is
separable again.

(c)

In the Gneiting scenario, the covariance has the following parametric form:

$c(t,s,t^{\prime},s^{\prime})=\frac{\sigma^{2}}{(a^{2}|t-t^{\prime}|^{2\alpha}+1)^{\tau}}\exp\left(\frac{b^{2}|s-s^{\prime}|^{2\gamma}}{(a^{2}|t-t^{\prime}|^{2\alpha}+1)^{\beta\gamma}}\right),$
(A.9)

where $a=b=\tau=\alpha=\gamma=\sigma^{2}=1$ and $\beta=0.7$. This covariance
is non-separable [14], but it is rather flat.

(d)

In the Fourier-Legendre scenario, we choose $a_{1}=a_{1}(t,t^{\prime})$ and
$b_{1}=b_{1}(s,s^{\prime})$ as the Fourier univariate covariances specified
above. Furthermore, $a_{2}=a_{2}(t,t^{\prime})$ and
$b_{2}=b_{2}(s,s^{\prime})$ are both chosen as rank-4 covariances with shifted
Legendre basis as their eigenfunctions, resulting in rather wiggly univariate
covariances (see Figure 8, right). The covariance is then chosen as
$c(t,s,t^{\prime},s^{\prime})=a_{1}(t,t^{\prime})\,b_{1}(s,s^{\prime})+a_{2}(t,t^{\prime})\,b_{2}(s,s^{\prime})$,
resulting in a non-separable covariance.

Fourier | Brownian | Legendre
---|---|---
|  |
Figure 8: Univariate covariances used as building blocks in our simulation
study.

###  A.5 Additional Simulation Results

Firstly, we show runtimes for the Fourier scenario from our simulation study
in Figure 4. The runtimes look similarly for any of the remaining scenarios
(not reported). To demonstrate effectiveness of the marginalization procedure
described in detail in Section A.2, we also show runtimes for the non-pooled
procedure, considering all raw covariances in sets 2.7 or 2.6 as separate
points for the purposes of smoothing. It leads to the same results as the
proposed approach but, as more and more points per surface are observed, the
number of scatter points supplied to the smoothing procedure increases
rapidly, which increases the runtimes. But more importantly, at the edge of
computational feasibility for the 4D smoothing approach (i.e. with the
percentage $p=10$) the proposed separable estimator is calculated about 200
times faster than the 4D smoothing estimator.

Secondly, it was observed in Section 4 that the proposed approach outperforms
the one-step version of our estimator $\widehat{a}_{0}\cdot\widehat{b}_{0}$,
cf. Figure 1. While the proposed approach can be seen as a two-step version, a
natural question arises whether a multi-step version of the estimator could
not be much better. Figure 9 compares the estimation errors achieved by the
proposed approach and by a three-step approach in all four scenarios
considered in Section 4. The third step offers a significant improvement in
only one of the scenarios, and even then the improvement is relatively small
compared to the improvements achieved in Section 4.

|
---|---
Figure 9: Left: Runtimes for Fourier simulations. Right: The proposed approach
(two-step) compared the proposed approach with included third step for all the
four scenarios considered in Section 4.

Finally, Table 2 shows relative estimation errors of the 4D smoothing approach
in all four scenarios used in Section 4, but only with the smallest considered
percentage $p=2$. The cross-validated choice of bandwidths is compared against
the choice of bandwidth suggested by the separable model (which is used in
Section 4). It is clear that cross-validation fails here, because even in the
simplest Brownian scenario, cross-validated relative error is larger that one.
The reason for that is likely the following. Cross-validation for 4D
smoothing, as implemented in the np package, does not evaluate its objective
function on a grid. In order to reduce the computational burden, the cross-
validation objective is optimized is a step-wise manner, until a stopping
criterion is met. The stopping criterion is set as a _tolerance_ (defaults to
$10^{-8}$), and the iterative optimization is stopped once both the change in
objective value and change in the bandwidths is smaller than the tolerance for
two consecutive iterations. While this saves computation time compared to
creating a grid over potential bandwidth values and evaluating the objective
function in all the grid points, the optimization can get stuck in a local
minimum. This is the reason why the cross-validated errors in Table 2 are so
high. The sampling pattern is very sparse with $p=2$, and there likely are
many local minima. However, the cross-validation still requires fitting the
covariance for different values of the bandwidths. We tried to obtain results
for larger values of $p$ as well, however we ran out of time (with a single
task) at 70 hours with $p=5$ even with the tolerance decreased to $10^{-2}$.
Hence, we have no other choice but to use the bandwidths suggested by the
separable model also for 4D smoothing.

Table 2: Relative estimation errors of 4D smoothing with cross-validated bandwidths and bandwidths suggested by the separable model for the four scenarios considered in Section 4 with $p=2$ percentages of the surfaces observed. | (a) Fourier | (b) Brownian | (c) Gneiting | (d) Fourier-Legendre
---|---|---|---|---
cross-validation | 1.07 | 1.06 | 1.04 | 1.39
separable choice | 0.04 | 0.26 | 1.34 | 1.1

##  Appendix B Proofs of Formal Statements

###  B.1 Explicit Formula for Local Polynomial Regression

Our estimators introduced in Section 2 are based on local polynomial
regression techniques and are defined as minimizers of (weighted) least
squares problems (2.4). It turns out that the minimizers for these point-wise
optimization problems admit a unique solution given by an explicit formula. In
this section we recall this formula for a general local linear polynomial
regression with possibly exogenous weights which will be later used in the
proofs of the asymptotic behaviour of our estimators.

The local linear surface smoother of the generic set
$\\{(x_{k},y_{k},z_{k})\;|\;k=1,\ldots,M\\}\subset\mathbb{R}^{3}$ given
weights $\\{w_{k}\;|\;k=1,\ldots,M\\}$ is defined as the solution of the least
squares problem (2.4). It turns out that this minimizer to this optimization
problem admits a unique solution:

$\widehat{\gamma}_{0(x,y)}=\Psi_{1}(x,y)\left[\Psi_{2}(x,y)\right]^{-1},\qquad(x,y)\in[0,1]^{2},$
(B.1)

where for each $(x,y)\in[0,1]^{2}$ and $p,q\in\mathbb{N}_{0}$ we define

$\displaystyle\Phi_{1}(x,y)$
$\displaystyle=S_{20}(x,y)S_{02}(x,y)-\left[S_{11}(x,y)\right]^{2},$
$\displaystyle\Phi_{2}(x,y)$
$\displaystyle=S_{10}(x,y)S_{02}(x,y)-S_{01}(x,y)S_{11}(x,y),$
$\displaystyle\Phi_{3}(x,y)$
$\displaystyle=S_{01}(x,y)S_{20}(x,y)-S_{10}(x,y)S_{11}(x,y),$
$\displaystyle\Psi_{1}(x,y)$
$\displaystyle=\Phi_{1}(x,y)Q_{00}(x,y)-\Phi_{2}(x,y)Q_{10}(x,y)-\Phi_{3}(x,y)Q_{01}(x,y),$
$\displaystyle\Psi_{2}(x,y)$
$\displaystyle=\Phi_{1}(x,y)S_{00}(x,y)-\Phi_{2}(x,y)S_{10}(x,y)-\Phi_{3}(x,y)S_{01}(x,y),$
$\displaystyle S_{pq}(x,y)$
$\displaystyle=\frac{1}{M}\sum_{k=1}^{M}\left(\frac{x-x_{k}}{h_{1}}\right)^{p}\left(\frac{y-y_{k}}{h_{2}}\right)^{q}\frac{1}{h_{1}h_{2}}K\left(\frac{x-x_{k}}{h_{1}}\right)K\left(\frac{y-y_{k}}{h_{2}}\right)w_{m},\qquad
0\leq p+q\leq 2,$ (B.2) $\displaystyle Q_{pq}(x,y)$
$\displaystyle=\frac{1}{M}\sum_{k=1}^{M}\left(\frac{x-x_{k}}{h_{1}}\right)^{p}\left(\frac{y-y_{k}}{h_{2}}\right)^{q}\frac{1}{h_{1}h_{2}}K\left(\frac{x-x_{k}}{h_{1}}\right)K\left(\frac{y-y_{k}}{h_{2}}\right)w_{k}z_{k},\qquad
0\leq p+q\leq 1,$ (B.3)

where $h_{1}>0$ and $h_{2}>0$ are smoothing bandwidths in the first and the
second dimension respectively.

The formula (B.1) is derived by differentiating the weighted least squares
(2.4) and finding the solution to the normal equations. It is based on the
standard steps used in the local regression literature, e.g. [11][§3.1] or
[33][§B.2].

###  B.2 Kernel Averages of $m$-dependent Data

Thanks the explicit formula (B.1) we may reduce the asymptotic behaviour
assessment to the investigation of the terms (B.2) and (B.3). In this section
we review the general asymptotic framework for the asymptotics of these kernel
averages and hence the framework for the asymptotics of (B.1). We shall use
the general theory developed by [17] who derived a toolbox for strong mixing
time series data where the regressors attain values in possibly unbounded
sets. Here we recall this result and write down a simplified version
sufficient for our data.

######  Theorem 3 ([17]).

Let $\\{(U_{k},V_{k},Z_{k})\\}_{k\in\mathbb{Z}}\in\mathbb{R}^{3}$ be a
strictly stationary sequence of random vectors and consider the averages of
the form

$\Xi_{k}(u,v)=\frac{1}{kh_{1}h_{2}}\sum_{i=1}^{k}\left(\frac{u-U_{i}}{h_{1}}\right)^{p}\left(\frac{v-V_{i}}{h_{2}}\right)^{q}K\left(\frac{u-U_{i}}{h_{1}}\right)K\left(\frac{v-V_{i}}{h_{2}}\right)Z_{i}$
(B.4)

where $p,q\in\mathbb{N}_{0}$.

  1. (C1)

$K(\cdot)$ is the Epanechnikov kernel, i.e.
$K(u)=(3/4)(1-u^{2})\mathds{1}_{[|u|<1]}$.

  2. (C2)

$(U_{k},V_{k})$ attain values in the set $[0,1]^{2}$.

  3. (C3)

$\\{(U_{k},V_{k},Z_{k})\\}$ is an $m$-dependent sequence, i.e. for each
$k\in\mathbb{Z}$, the random vectors $(\dots,U_{k},V_{k},Z_{k})$ and
$(U_{k+m},V_{k+m},Z_{k+m},U_{k+m+1},V_{k+m+1},Z_{k+m+1},\dots)$ are
independent.

  4. (C4)

There exists $s>2$ such that
$(u,v)\mapsto\mathbb{E}[|Z_{1}|^{s}|U_{1}=u,V_{1}=v]$ is bounded.

  5. (C5)

The probability density function of $(U_{1},V_{1})$ is twice continuously
differentiable.

  6. (C6)

The smoothing bandwidth satisfies $(\log k)/(kh_{1}h_{2})=o(1)$ as
$k\to\infty$.

Then the kernel averages (B.4) satisfy

$\sup_{(u,v)\in[0,1]^{2}}\left|\Xi_{k}(u,v)-\mathbb{E}\Xi_{k}(u,v)\right|=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
k}{kh_{1}h_{2}}}\right),\qquad\text{as}\quad k\to\infty.$ (B.5)

Theorem 3 is essentially a special case of [17][Thm 2] when the consider
sequence is defined on a bounded domain and is $m$-dependent. We also assume
that the two dimensional smoothing kernel (B.4) is the product of two one
dimensional Epanechnikov kernels. Note that the function $(u,v)\mapsto
u^{p}v^{q}K(u)K(v)$ satisfies Hansen’s conditions on the smoothing kernel.

The only generalisation where Theorem 3 deviates from [17][Thm 2] is that we
allow the smoothing bandwidth $(h_{1},h_{2})$ to attain different values in
different directions. The proof of such generalisation, while having the
smoothing kernel as a product of one-dimensional kernels, follows the lines of
the proof [17][Thm 2] where $h^{d}$, with $d=2$, is replaced by $h_{1}h_{2}$.

The following corollary goes one step further and incorporates the convergence
of $\mathbb{E}\Xi_{k}(u,v)$ into the statement (B.5).

######  Corollary 2.

Under the assumption of Theorem 3 consider the function

$M(u,v)=c_{p}c_{q}\mathbb{E}\left[Z_{1}\middle|U_{1}=u,V_{1}=v\right]f_{U_{1},V_{1}}(u,v),\qquad
x,y\in[0,1],$ (B.6)

where $f_{U_{1},V_{1}}(\cdot,\cdot)$ denotes the probability density function
of $(U_{1},V_{1})$ and $c_{r}=\int x^{r}K(x)dx$ for $r\in\mathbb{N}_{0}$.
Moreover:

  1. (D1)

Assume that the function $M(\cdot,\cdot)$ is twice continuously differentiable
on $[0,1]^{2}$.

Then the kernel averages (B.4) satisfy

$\sup_{(u,v)\in[0,1]^{2}}\left|\Xi_{k}(u,v)-M(u,v)\right|=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
k}{kh_{1}h_{2}}}+h_{1}^{2}+h_{2}^{2}\right),\qquad\text{as}\quad k\to\infty.$
(B.7)

###### Proof.

We start by decomposing the supremum (B.7) into a stochastic and a
deterministic part

$\sup_{(u,v)\in[0,1]^{2}}\left|\Xi_{k}(u,v)-M(u,v)\right|\leq\sup_{(u,v)\in[0,1]^{2}}\left|\Xi_{k}(u,v)-\mathbb{E}\Xi_{k}(u,v)\right|+\sup_{(u,v)\in[0,1]^{2}}\left|\mathbb{E}\Xi_{k}(u,v)-M(u,v)\right|.$
(B.8)

The first term on the right-hand side of (B.8) is of order
$O_{\mathbb{P}}(\sqrt{(\log n)/(nh_{1}h_{2})})$ by Theorem 3. The expectation
in the second term on the right-hand side of (B.8) is developed as

$\displaystyle\mathbb{E}\Xi_{k}(u,v)$
$\displaystyle=\mathbb{E}\left[\left(\frac{u-U_{1}}{h_{1}}\right)^{p}\left(\frac{v-V_{1}}{h_{2}}\right)^{q}\frac{1}{h_{1}h_{2}}K\left(\frac{u-U_{1}}{h_{1}}\right)K\left(\frac{v-V_{1}}{h_{2}}\right)\mathbb{E}\left[Z_{11}\middle|U_{1},V_{1}\right]\right]$
$\displaystyle=\mathbb{E}\left[\left(\frac{u-U_{1}}{h_{1}}\right)^{p}\left(\frac{v-V_{1}}{h_{2}}\right)^{q}\frac{1}{h_{1}h_{2}}K\left(\frac{u-U_{1}}{h_{1}}\right)K\left(\frac{v-V_{1}}{h_{2}}\right)M(U_{1},V_{1})\right]$
$\displaystyle=\iint\left(\frac{u-x}{h_{1}}\right)^{p}\left(\frac{v-y}{h_{2}}\right)^{q}\frac{1}{h_{1}h_{2}}K\left(\frac{u-x}{h_{1}}\right)K\left(\frac{v-y}{h_{2}}\right)M(x,y)f_{U_{1},V_{1}}(x,y)dxdy$
$\displaystyle=\iint\left(\tilde{x}\right)^{p}\left(\tilde{y}\right)^{q}K\left(\tilde{x}\right)K\left(\tilde{y}\right)M(u+h_{1}\tilde{x},v+h_{2}\tilde{y})f_{U_{1},V_{1}}(u+h_{1}\tilde{x},v+h_{2}\tilde{y})d\tilde{x}d\tilde{y}.$
(B.9)

Applying the Taylor expansion of order 2 in the right-hand side of (B.9) and
using the assumptions (C5) and (D1), the second term on the right-hand side of
(B.8) is of order $O(h_{1}^{2}+h_{2}^{2})$. ∎

###  B.3 Proof of Proposition 1

Tailoring the generic smoother (2.4) to the mean surface estimator (2.5), we
arrive at the customized versions of (B.2) and (B.3):

$\displaystyle S_{pq}^{\mu}(t,s)$
$\displaystyle=\frac{1}{\sum_{n=1}^{N}M_{n}}\sum_{n=1}^{N}\sum_{m=1}^{M_{n}}\left(\frac{t-t_{nm}}{h_{\mu,1}}\right)^{p}\left(\frac{s-s_{nm}}{h_{\mu,2}}\right)^{q}\cdot$
(B.10)
$\displaystyle\qquad\cdot\frac{1}{h_{\mu,1}h_{\mu,2}}K\left(\frac{t-t_{nm}}{h_{\mu,1}}\right)K\left(\frac{s-s_{nm}}{h_{\mu,2}}\right),\quad
0\leq p+q\leq 2.$ (B.11) $\displaystyle Q_{pq}^{\mu}(t,s)$
$\displaystyle=\frac{1}{\sum_{n=1}^{N}M_{n}}\sum_{n=1}^{N}\sum_{m=1}^{M_{n}}\left(\frac{t-t_{nm}}{h_{\mu,1}}\right)^{p}\left(\frac{s-s_{nm}}{h_{\mu,2}}\right)^{q}\cdot$
(B.12)
$\displaystyle\qquad\cdot\frac{1}{h_{\mu,1}h_{\mu,2}}K\left(\frac{t-t_{nm}}{h_{\mu,1}}\right)K\left(\frac{s-s_{nm}}{h_{\mu,2}}\right)Y_{nm},\quad
0\leq p+q\leq 1.$ (B.13)

We assess the asymptotic behaviour of (B.10) and (B.12) in the following
lemmas.

######  Lemma 1.

Under the assumptions (B1) – (B6),

$\sup_{(t,s)\in[0,1]^{2}}\left|Q^{\mu}_{pq}(t,s)-M_{[Q^{\mu}_{pq}]}(t,s)\right|=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh_{1}h_{2}}}+h_{\mu,1}^{2}+h_{\mu,2}^{2}\right),\qquad\text{as}\quad
N\to\infty,$ (B.14)

for each $0\leq p+q\leq 1$ and where

$M_{[Q^{\mu}_{00}]}(t,s)=\mu(t,s)f_{(t,s)}(t,s),\qquad
M_{[Q^{\mu}_{10}]}(t,s)=M_{[Q^{\mu}_{01}]}(t,s)=0,\qquad t,s\in[0,1].$ (B.15)

###### Proof.

Define the sequence of random vectors
$\\{(U_{k},V_{k},Z_{k})\\}_{k=1}^{\infty}$ by putting $\\{t_{nm}\\}$,
$\\{s_{nm}\\}$ and $\\{Y_{nm}\\}$ in order such that

$\displaystyle\\{U_{1},U_{2},\dots\\}$
$\displaystyle=\\{t_{11},t_{12},\dots,t_{1m_{1}},t_{21},\dots,t_{2m_{2}},t_{31},\dots\\},$
$\displaystyle\\{V_{1},V_{2},\dots\\}$
$\displaystyle=\\{s_{11},s_{12},\dots,s_{1m_{1}},s_{21},\dots,s_{2m_{2}},s_{31},\dots\\},$
$\displaystyle\\{Z_{1},Z_{2},\dots\\}$
$\displaystyle=\\{Y_{11},Y_{12},\dots,Y_{1m_{1}},Y_{21},\dots,Y_{2m_{2}},Y_{31},\dots\\}.$
(B.16)

The sequence$\\{(U_{k},V_{k},Z_{k})\\}_{k=1}^{\infty}$ satisfies the
assumption of Theorem 3 and Corollary 2, namely strict stationarity is by
assumptions (B1) – (B3), is $M^{max}$-dependent by assumption (B1), and the
conditions (C4), (C5), (C6), (D1) are satisfied by the assumptions (B5), (B2),
(B6), (B5) respectively. Therefore the sequence of kernel averages

$\Xi^{\mu}_{pq,k}(t,s)=\frac{1}{kh_{\mu,1}h_{\mu,2}}\sum_{i=1}^{n}\left(\frac{t-U_{i}}{h_{1}}\right)^{p}\left(\frac{s-V_{i}}{h_{2}}\right)^{q}K\left(\frac{t-U_{i}}{h_{1}}\right)K\left(\frac{s-V_{i}}{h_{2}}\right)Z_{i}$

satisfies

$\sup_{(t,s)\in[0,1]^{2}}\left|\Xi^{\mu}_{pq,k}(t,s)-M_{[Q^{\mu}_{pq}]}(t,s)\right|=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
k}{kh_{1}h_{2}}}+h_{\mu,1}^{2}+h_{\mu,2}^{2}\right),\qquad\text{as}\quad
k\to\infty,$

and the formulae (B.15) follow from the definition of
$\\{(U_{k},V_{k},Z_{k})\\}_{k=1}^{\infty}$ and definition (B.6). Since the
sequence $\\{Q^{\mu}_{pq}(t,s)\\}_{N=1}^{\infty}$ is a subsequence of
$\\{\Xi^{\mu}_{pq,k}(t,s)\\}_{k=1}^{\infty}$ and $k=k(N)\asymp N$ as
$N\to\infty$, the convergence rate (B.14) holds as well. ∎

######  Lemma 2.

Under the assumptions (B1) – (B3),

$\sup_{(t,s)\in[0,1]^{2}}\left|S^{\mu}_{pq}(t,s)-M_{[S^{\mu}_{pq}]}(t,s)\right|=\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh_{1}h_{2}}}+h_{\mu,1}^{2}+h_{\mu,2}^{2}\right),\qquad\text{as}\quad
N\to\infty,$ (B.17)

for each $0\leq p+q\leq 1$ and where

$\displaystyle M_{[S^{\mu}_{00}]}(t,s)$
$\displaystyle=\mu(t,s)f_{(t,s)}(t,s),$ $\displaystyle
M_{[S^{\mu}_{11}]}(t,s)$
$\displaystyle=M_{[S^{\mu}_{10}]}(t,s)=M_{[S^{\mu}_{01}]}(t,s)=0,$ (B.18)
$\displaystyle M_{[S^{\mu}_{20}]}(t,s)$
$\displaystyle=M_{[S^{\mu}_{02}]}(t,s)=c_{2}f_{(t,s)}(t,s),$ $\displaystyle
c_{2}$ $\displaystyle=\int x^{2}K(x)dx.\qquad t,s\in[0,1].$ (B.19)

###### Proof.

The proof of this lemma follows essentially the same lines as the proof of
Lemma 1. In the definition of the sequence
$\\{(U_{k},V_{k},Z_{k})\\}_{k=1}^{\infty}$ we put $Z_{k}=1$, for all
$k\in\mathbb{N}$, on the line (B.16). The formulae (B.18) and (B.19) are
verified analogously by the definition (B.6). ∎

###### Proof of Proposition 1.

We are now ready to combine the above and proof Proposition 1. Following the
formulae presented in Section B.1, we have for

$\displaystyle\Phi^{\mu}_{1}(t,s)$
$\displaystyle=S^{\mu}_{20}(t,s)S^{\mu}_{02}(t,s)-\left[S^{\mu}_{11}(t,s)\right]^{2},$
$\displaystyle\Phi^{\mu}_{2}(t,s)$
$\displaystyle=S^{\mu}_{10}(t,s)S^{\mu}_{02}(t,s)-S^{\mu}_{01}(t,s)S^{\mu}_{11}(t,s),$
$\displaystyle\Phi^{\mu}_{3}(t,s)$
$\displaystyle=S^{\mu}_{01}(t,s)S^{\mu}_{20}(t,s)-S^{\mu}_{10}(t,s)S^{\mu}_{11}(t,s),$
$\displaystyle\Psi^{\mu}_{1}(t,s)$
$\displaystyle=\Phi_{1}(t,s)Q^{\mu}_{00}(t,s)-\Phi_{2}(t,s)Q^{\mu}_{10}(t,s)-\Phi_{3}(t,s)Q^{\mu}_{01}(t,s),$
$\displaystyle\Psi^{\mu}_{2}(t,s)$
$\displaystyle=\Phi_{1}(t,s)S^{\mu}_{00}(t,s)-\Phi_{2}(t,s)S^{\mu}_{10}(t,s)-\Phi_{3}(t,s)S^{\mu}_{01}(t,s),$

their asymptotic behaviour

$\displaystyle\Phi^{\mu}_{1}(t,s)$
$\displaystyle=\left(c_{2}f_{(t,s)}(t,s)\right)^{2}+\mathcal{O}_{\mathbb{P}}\left(r^{\mu}_{N}\right),$
$\displaystyle\Phi^{\mu}_{2}(t,s)$
$\displaystyle=\mathcal{O}_{\mathbb{P}}\left(r^{\mu}_{N}\right),\qquad\qquad\Phi^{\mu}_{3}(t,s)=\mathcal{O}_{\mathbb{P}}\left(r^{\mu}_{N}\right),$
$\displaystyle\Psi^{\mu}_{1}(t,s)$
$\displaystyle=\left(c_{2}\right)^{2}\left(f_{(t,s)}(t,s)\right)^{3}\mu(t,s)+\mathcal{O}_{\mathbb{P}}\left(r^{\mu}_{N}\right),$
$\displaystyle\Psi^{\mu}_{2}(t,s)$
$\displaystyle=\left(c_{2}\right)^{2}\left(f_{(t,s)}(t,s)\right)^{3}+\mathcal{O}_{\mathbb{P}}\left(r^{\mu}_{N}\right),$

uniformly in $(t,s)\in[0,1]^{2}$, as $N\to\infty$, where
$r^{\mu}_{N}=\sqrt{(\log
N)/(Nh_{\mu,1}h_{\mu,2})}+h_{\mu,1}^{2}+h_{\mu,2}^{2}$. Hence

$\widehat{\mu}(t,s)=\Psi^{\mu}_{1}(t,s)/\left[\Psi^{\mu}_{2}(t,s)\right]^{-1}=\mu(t,s)+\mathcal{O}_{\mathbb{P}}\left(r^{\mu}_{N}\right),\qquad\text{as}\quad
N\to\infty,$

by the uniform version of Slutsky’s theorem and by the fact that
$f_{(t,s)}(\cdot,\cdot)\neq 0$ on $[0,1]^{2}$ by (B2). Thanks to the
assumption (B6), the convergence rate simplifies to the common rate $h$ and

$\widehat{\mu}(t,s)=\mu(t,s)+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right),\qquad\text{as}\quad N\to\infty.$

∎

###  B.4 Proof of Theorem 1 and Corollary 1

######  Lemma 3.

Assume the conditions (A1), (B1) – (B9) and fix a deterministic twice
continuously differentiable kernel
$\beta(s,s^{\prime}),\,s,s^{\prime}\in[0,1]$ such that
$\iint[\beta(s,s^{\prime})]^{2}dsds^{\prime}>0$. Then the smoother
$\widehat{\alpha}(t,t^{\prime}),\,t,t^{\prime}\in[0,1]$, obtained by smoothing
the set

$\left\\{\left(t_{nm},t_{nm^{\prime}},\frac{G_{nmm^{\prime}}}{\beta(s_{nm},s_{nm^{\prime}})}\right)\;\Bigg{|}\;m,m^{\prime}=1,\ldots,M_{n},\,m\neq
m^{\prime},\;n=1\ldots,N\right\\}$ (B.20)

using weights $\\{\beta^{2}(s_{nm},s_{nm^{\prime}})\\}$ admits the following
asymptotics

$\widehat{\alpha}(t,t^{\prime})=a(t,t^{\prime})\frac{\iint\beta(s,s^{\prime})b(s,s^{\prime})f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}}{\iint\left[\beta(s,s^{\prime})\right]^{2}f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}}+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)$ (B.21)

uniformly in $(t,t^{\prime})\in[0,1]^{2}$ as $N\to\infty$. We recall that
$a(t,t^{\prime}),\,t,t^{\prime}\in[0,1],$ and
$b(s,s^{\prime}),\,s,s^{\prime}\in[0,1]$, on the right-hand side of (B.21) are
the true separable covariance structure components (2.1).

###### Proof.

We start the proof by the analysis of a simplified case. Suppose that we know
the mean surface $\mu(\cdot,\cdot)$ and define the raw covariances accordingly

$\tilde{G}_{nmm^{\prime}}=(Y_{nm}-\mu(t_{nm},s_{nm}))(Y_{nm^{\prime}}-\mu(t_{nm^{\prime}},s_{nm^{\prime}})).$
(B.22)

Construct the smoother of the set (B.20) where we replace $G_{nmm^{\prime}}$
by $\tilde{G}_{nmm^{\prime}}$. Such smoother, denoted as
$\tilde{\alpha}(\cdot,\cdot)$, is given by the formula

$\displaystyle\tilde{\alpha}(t,t^{\prime})$
$\displaystyle=\Psi^{\tilde{\alpha}}_{1}(t,t^{\prime})/\left[\Psi^{\tilde{\alpha}}_{2}(t,t^{\prime})\right]^{-1},$
(B.23) $\displaystyle\Psi^{\tilde{\alpha}}_{1}(t,t^{\prime})$
$\displaystyle=\Phi_{1}(t,t^{\prime})Q^{\tilde{\alpha}}_{00}(t,t^{\prime})-\Phi_{2}(t,t^{\prime})Q^{\tilde{\alpha}}_{10}(t,t^{\prime})-\Phi_{3}(t,t^{\prime})Q^{\tilde{\alpha}}_{01}(t,t^{\prime}),$
$\displaystyle\Psi^{\tilde{\alpha}}_{2}(t,t^{\prime})$
$\displaystyle=\Phi_{1}(t,t^{\prime})S^{\tilde{\alpha}}_{00}(t,t^{\prime})-\Phi_{2}(t,t^{\prime})S^{\tilde{\alpha}}_{10}(t,t^{\prime})-\Phi_{3}(t,t^{\prime})S^{\tilde{\alpha}}_{01}(t,t^{\prime}),$
$\displaystyle\Phi^{\tilde{\alpha}}_{1}(t,t^{\prime})$
$\displaystyle=S^{\tilde{\alpha}}_{20}(t,t^{\prime})S^{\tilde{\alpha}}_{02}(t,t^{\prime})-\left[S^{\tilde{\alpha}}_{11}(t,t^{\prime})\right]^{2},$
$\displaystyle\Phi^{\tilde{\alpha}}_{2}(t,t^{\prime})$
$\displaystyle=S^{\tilde{\alpha}}_{10}(t,t^{\prime})S^{\tilde{\alpha}}_{02}(t,t^{\prime})-S^{\tilde{\alpha}}_{01}(t,t^{\prime})S^{\tilde{\alpha}}_{11}(t,t^{\prime}),$
$\displaystyle\Phi^{\tilde{\alpha}}_{3}(t,t^{\prime})$
$\displaystyle=S^{\tilde{\alpha}}_{01}(t,t^{\prime})S^{\tilde{\alpha}}_{20}(t,t^{\prime})-S^{\tilde{\alpha}}_{10}(t,t^{\prime})S^{\tilde{\alpha}}_{11}(t,t^{\prime}),$
$\displaystyle S_{pq}^{\tilde{\alpha}}(t,t^{\prime})$
$\displaystyle=\frac{1}{\sum_{n=1}^{N}M_{n}(M_{n}-1)}\sum_{n=1}^{N}\sum_{\begin{subarray}{c}m,m^{\prime}=1\\\
m\neq
m^{\prime}\end{subarray}}^{M_{n}}\left(\frac{t-t_{nm}}{h_{a}}\right)^{p}\left(\frac{t^{\prime}-t_{nm^{\prime}}}{h_{a}}\right)^{q}\cdot$
$\displaystyle\qquad\cdot\frac{1}{h_{a}^{2}}K\left(\frac{t-t_{nm}}{h_{a}}\right)K\left(\frac{t^{\prime}-t_{nm^{\prime}}}{h_{a}}\right)\left[\beta(s_{nm},s_{nm^{\prime}})\right]^{2},\qquad
0\leq p+q\leq 2,$ (B.24) $\displaystyle Q_{pq}^{\tilde{\alpha}}(t,t^{\prime})$
$\displaystyle=\frac{1}{\sum_{n=1}^{N}M_{n}(M_{n}-1)}\sum_{n=1}^{N}\sum_{\begin{subarray}{c}m,m^{\prime}=1\\\
m\neq
m^{\prime}\end{subarray}}^{M_{n}}\left(\frac{t-t_{nm}}{h_{a}}\right)^{p}\left(\frac{t^{\prime}-t_{nm^{\prime}}}{h_{a}}\right)^{q}\cdot$
$\displaystyle\qquad\cdot\frac{1}{h_{a}^{2}}K\left(\frac{t-t_{nm}}{h_{a}}\right)K\left(\frac{t^{\prime}-t_{nm^{\prime}}}{h_{a}}\right)\beta(s_{nm},s_{nm^{\prime}})\tilde{G}_{nmm^{\prime}},\qquad
0\leq p+q\leq 1.$ (B.25)

The asymptotic behaviour of $Q_{pq}^{\tilde{\alpha}}$ and
$S_{pq}^{\tilde{\alpha}}(t,t^{\prime})$ is assessed similarly as the surface
smoother in Lemmas 1 and 2. We proceed again with defining the
$[M^{max}]^{2}$-dependent sequences $\\{(U_{k},V_{k},Z_{k})\\}_{k=1}^{\infty}$
by putting the pairs
$(t_{nm},t_{nm^{\prime}}),\,m,m^{\prime}=1,\dots,M_{n},\,m\neq
m^{\prime},\,n=1,\dots,N$ into the sequence
$\\{(U_{k},V_{k})\\}_{k=1}^{\infty}$ such that we set $U_{k}=t_{nm}$ and
$V_{k}=t_{nm^{\prime}}$ while starting from the data from the first surface
($n=1$), then proceeding with $n=2$ etc.

For the asymptotics of $Q_{pq}^{\tilde{\alpha}}$, define

$Z^{Q}_{k}=\beta(s_{nm},s_{nm^{\prime}})\tilde{G}_{nmm^{\prime}}$ (B.26)

where $s_{nm},s_{nm^{\prime}},\tilde{G}_{nmm^{\prime}}$ correspond to that
sparse observation which was assigned to $(U_{k},V_{k})$. We use Theorem 3 and
Corollary 2 thanks to the assumptions (B7) – (B9). Moreover, we verify

$\displaystyle\mathbb{E}\left[Z^{Q}_{k}\middle|U_{k}=t,V_{k}=t^{\prime}\right]$
$\displaystyle=\mathbb{E}\left[\beta(s_{nm},s_{nm^{\prime}})\left(X_{n}(t,s_{nm})+\varepsilon_{nm}\right)\left(X_{n}(t,s_{nm^{\prime}})+\varepsilon_{nm^{\prime}})\right)\middle|U_{k}=t,V_{k}=t^{\prime}\right]$
$\displaystyle=\mathbb{E}\left[\beta(s_{nm},s_{nm^{\prime}})a(t,t^{\prime})b(s_{nm},s_{nm^{\prime}})\right]$
$\displaystyle=a(t,t^{\prime})\iint\beta(s,s^{\prime})b(s,s^{\prime})f_{s}(s)f_{s}(s^{\prime})dsds^{\prime},$

where $f_{s}(s)=\int f_{(t,s)}(t,s)dt$ is the marginal density of the random
position $s_{11}$. Likewise, denote $f_{t}(t)=\int f_{(t,s)}(t,s)ds$ is the
marginal density of the random position $t_{11}$. Then

$Q_{00}^{\tilde{\alpha}}(t,t^{\prime})=a(t,t^{\prime})f_{t}(t)f_{t}(t^{\prime})\iint\beta(s,s^{\prime})b(s,s^{\prime})f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh_{a}^{2}}}+h_{a}^{2}\right),$

uniformly in $(t,t^{\prime})\in[0,1]^{2}$ as $N\to\infty$ and where
$c_{r}=\int x^{r}K(x)dx,\,r\in\mathbb{N}$.

Similarly to the analysis above we assess the asymptotics of
$S_{pq}^{\tilde{\alpha}}$. Instead of the definition in (B.26) we set here
$Z^{S}_{k}=\left[\beta(s_{nm},s_{nm^{\prime}})\right]^{2}$ and calculate

$\mathbb{E}\left[Z^{S}_{k}\middle|U_{k}=t,V_{k}=t^{\prime}\right]=\iint[\beta(s,s^{\prime})]^{2}f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}.$

Hence

$S_{pq}^{\tilde{\alpha}}(t,t^{\prime})=c_{p}c_{q}f_{t}(t)f_{t}(t^{\prime})\iint[\beta(s,s^{\prime})]^{2}f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh_{a}^{2}}}+h_{a}^{2}\right)$

uniformly in $(t,t^{\prime})\in[0,1]^{2}$ as $N\to\infty$. By the assumptions
on the kernel $\beta(\cdot,\cdot)$, the uniform Slutsky theorem, the formula
(B.23), and the fact that $h_{a}\asymp h$ (assumption (B9)):

$\tilde{\alpha}(t,t^{\prime})=a(t,t^{\prime})\frac{\iint\beta(s,s^{\prime})b(s,s^{\prime})f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}}{\iint\left[\beta(s,s^{\prime})\right]^{2}f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}}+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)$

uniformly in $(t,t^{\prime})\in[0,1]^{2}$ as $N\to\infty$.

It remains to comment on the difference $\tilde{\alpha}(t,t^{\prime})$ and
$\widehat{\alpha}(t,t^{\prime})$, i.e. when the empirical mean
$\widehat{\mu}(\cdot,\cdot)$ is supplied into the raw covariances
$G_{nmm^{\prime}}$. Since

$\displaystyle G_{nmm^{\prime}}$ $\displaystyle=\tilde{G}_{nmm^{\prime}}$
$\displaystyle+\left(\mu(t_{nm},s_{nm})-\widehat{\mu}(t_{nm},s_{nm})\right)\left(Y_{nm^{\prime}}-\widehat{\mu}(t_{nm^{\prime}},s_{nm^{\prime}})\right)$
$\displaystyle+\left(\mu(t_{nm^{\prime}},s_{nm^{\prime}})-\widehat{\mu}(t_{nm^{\prime}},s_{nm^{\prime}})\right)\left(Y_{nm}-\widehat{\mu}(t_{nm},s_{nm})\right)$
$\displaystyle+\left(\mu(t_{nm},s_{nm})-\widehat{\mu}(t_{nm},s_{nm})\right)\left(\mu(t_{nm^{\prime}},s_{nm^{\prime}})-\widehat{\mu}(t_{nm^{\prime}},s_{nm^{\prime}})\right)$

we conclude by Proposition 1 that

$G_{nmm^{\prime}}=\tilde{G}_{nmm^{\prime}}+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)$ (B.27)

uniformly across all $n,m,m^{\prime}$ as $N\to\infty$. Therefore the claim
(B.21) follows. ∎

######  Corollary 3.

Assume the conditions (A1), (B1) – (B9) and consider a random kernel
$\widehat{\beta}(\cdot,\cdot)$ such that

$\widehat{\beta}(s,s^{\prime})=\beta(s,s^{\prime})+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)$ (B.28)

uniformly in $(s,s^{\prime})\in[0,1]^{2}$ as $N\to\infty$, where
$\beta(s,s^{\prime}),\,s,s^{\prime}\in[0,1]$ is a deterministic twice
continuously differentiable kernel such that
$\iint[\beta(s,s^{\prime})]^{2}dsds^{\prime}>0$. Then the smoother
$\widehat{\alpha}(t,t^{\prime}),\,t,t^{\prime}\in[0,1]$, obtained by smoothing
the set

$\left\\{\left(t_{nm},t_{nm^{\prime}},\frac{G_{nmm^{\prime}}}{\widehat{\beta}(s_{nm},s_{nm^{\prime}})}\right)\;\Bigg{|}\;m,m^{\prime}=1,\ldots,M_{n},\,m\neq
m^{\prime},\;n=1\ldots,N\right\\}$

using weights $\\{\widehat{\beta}^{2}(s_{nm},s_{nm^{\prime}})\\}$ admits the
same asymptotics as in the previous lemma:

$\widehat{\alpha}(t,t^{\prime})=a(t,t^{\prime})\frac{\iint\beta(s,s^{\prime})b(s,s^{\prime})f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}}{\iint\left[\beta(s,s^{\prime})\right]^{2}f_{s}(s)f_{s}(s^{\prime})dsds^{\prime}}+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)$

uniformly in $(t,t^{\prime})\in[0,1]^{2}$ as $N\to\infty$.

###### Proof.

The proof of this corollary follows the same lines as the proof of Lemma 3. We
define

$\displaystyle S_{pq}^{\widehat{\alpha}}(t,t^{\prime})$
$\displaystyle=\frac{1}{\sum_{n=1}^{N}M_{n}(M_{n}-1)}\sum_{n=1}^{N}\sum_{\begin{subarray}{c}m,m^{\prime}=1\\\
m\neq
m^{\prime}\end{subarray}}^{M_{n}}\left(\frac{t-t_{nm}}{h_{a}}\right)^{p}\left(\frac{t^{\prime}-t_{nm^{\prime}}}{h_{a}}\right)^{q}\cdot$
$\displaystyle\qquad\cdot\frac{1}{h_{a}^{2}}K\left(\frac{t-t_{nm}}{h_{a}}\right)K\left(\frac{t^{\prime}-t_{nm^{\prime}}}{h_{a}}\right)\left[\widehat{\beta}(s_{nm},s_{nm^{\prime}})\right]^{2},\qquad
0\leq p+q\leq 2,$ $\displaystyle Q_{pq}^{\widehat{\alpha}}(t,t^{\prime})$
$\displaystyle=\frac{1}{\sum_{n=1}^{N}M_{n}(M_{n}-1)}\sum_{n=1}^{N}\sum_{\begin{subarray}{c}m,m^{\prime}=1\\\
m\neq
m^{\prime}\end{subarray}}^{M_{n}}\left(\frac{t-t_{nm}}{h_{a}}\right)^{p}\left(\frac{t^{\prime}-t_{nm^{\prime}}}{h_{a}}\right)^{q}\cdot$
$\displaystyle\qquad\cdot\frac{1}{h_{a}^{2}}K\left(\frac{t-t_{nm}}{h_{a}}\right)K\left(\frac{t^{\prime}-t_{nm^{\prime}}}{h_{a}}\right)\widehat{\beta}(s_{nm},s_{nm^{\prime}})\tilde{G}_{nmm^{\prime}},\qquad
0\leq p+q\leq 1,$

as analogues of (B.24) and (B.25). Thanks to the assumption (B.28), the
difference in asymptotically negligible

$\displaystyle S_{pq}^{\widehat{\alpha}}(t,t^{\prime})$
$\displaystyle=S_{pq}^{\tilde{\alpha}}(t,t^{\prime})+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right),$ $\displaystyle
Q_{pq}^{\widehat{\alpha}}(t,t^{\prime})$
$\displaystyle=Q_{pq}^{\tilde{\alpha}}(t,t^{\prime})+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right),$

uniformly in $(t,t^{\prime})\in[0,1]^{2}$ as $N\to\infty$. The rest of the
proof follows from the proof of Lemma 3. ∎

We are now ready to prove our main result.

###### Proof of Theorem 1.

The proof is now quite a simple application of Lemma 3 and Corollary 3. First
note that even though these results are formulated for the estimation of the
covariance kernel $a(\cdot,\cdot)$, they can be likewise applied for the
estimation of $b(\cdot,\cdot)$ due to their symmetry in the separable model
(2.1).

The estimator $\widehat{a}_{0}(\cdot,\cdot)$ is realised by smoothing the raw
covariances $G_{nmm^{\prime}}$ without any weights, thus corresponding to the
initial guess $\beta(s,s^{\prime})\equiv 1,\,s,s^{\prime}\in[0,1]$. Therefore
its asymptotic behaviour is by Lemma 3:

$\widehat{a}_{0}(t,t^{\prime})=\Theta
a(t,t^{\prime})+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)$

uniformly in $(t,t^{\prime})\in[0,1]^{2}$ as $N\to\infty$ where $\Theta$ is
defined in (3.1).

Now, applying Corollary 3 three times and by the assumption (B10) we obtain

$\displaystyle\widehat{b}_{0}(s,s^{\prime})$
$\displaystyle=\frac{1}{\Theta}b(s,s^{\prime})+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right),$ $\displaystyle\widehat{a}(t,t^{\prime})$
$\displaystyle=\Theta
a(t,t^{\prime})+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right),$ $\displaystyle\widehat{a}(s,s^{\prime})$
$\displaystyle=\frac{1}{\Theta}b(s,s^{\prime})+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right),$

uniformly in $(t,t^{\prime})\in[0,1]^{2}$ or $(s,s^{\prime})\in[0,1]^{2}$, as
$N\to\infty$. ∎

###### Proof of Corollary 1.

This corollary follows directly by applying Theorem 1 onto the right hand side
of:

$\left|\widehat{a}(t,t^{\prime})\widehat{b}(s,s^{\prime})-a(t,t^{\prime})b(s,s^{\prime})\right|\leq\left|\widehat{a}(t,t^{\prime})-\Theta\widehat{a}(t,t^{\prime})\right|\left|\widehat{b}(s,s^{\prime})\right|+\left|a(t,t^{\prime})\right|\Theta\left|\widehat{b}(s,s^{\prime})-\frac{1}{\Theta}b(s,s^{\prime})\right|.$

∎

###  B.5 Proof of Proposition 2

The noise level estimator asymptotic behaviour is treated analogously to
previous estimators of the mean surface $\mu(\cdot,\cdot)$ and the covariance
kernels $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$.

The estimator $\widehat{V}(t,s)$ is formed by smoothing the raw covariances
$G_{nmm}$ agains $(t_{nm},s_{nm})$. for $m=1,\dots,M_{n},\,n=1,\dots,N$.
Therefore we form the sequence of vectors
$\\{(U_{k},V_{k},Z_{k})\\}_{k=1}^{\infty}$ by putting $\\{t_{nm}\\}$,
$\\{s_{nm}\\}$ and $\\{\tilde{G}_{nmm}\\}$ (defined in (B.22)) in order such
that

$\displaystyle\\{U_{1},U_{2},\dots\\}$
$\displaystyle=\\{t_{11},t_{12},\dots,t_{1m_{1}},t_{21},\dots,t_{2m_{2}},t_{31},\dots\\},$
$\displaystyle\\{V_{1},V_{2},\dots\\}$
$\displaystyle=\\{s_{11},s_{12},\dots,s_{1m_{1}},s_{21},\dots,s_{2m_{2}},s_{31},\dots\\},$
$\displaystyle\\{Z_{1},Z_{2},\dots\\}$
$\displaystyle=\\{\tilde{G}_{111},\tilde{G}_{122},\dots,\tilde{G}_{1m_{1}m_{1}},\tilde{G}_{211},\dots,\tilde{G}_{2m_{2}m_{2}},\tilde{G}_{311},\dots\\}.$

By following the steps of the proof of Lemma 1 or Lemma 3. Verifying
$\mathbb{E}\left[Z_{1}\middle|U_{1}=t,V_{1}=s\right]=a(t,t)b(s,s)+\sigma^{2}$
for $t,s\in[0,1]$, the asymptotic equivalence (B.27), and the assumption (B11)
implies

$\widehat{V}(t,s)=a(t,t)b(s,s)+\sigma^{2}+\mathcal{O}_{\mathbb{P}}\left(\sqrt{\frac{\log
N}{Nh^{2}}}+h^{2}\right)$

uniformly in $(t,s)\in[0,1]^{2}$ as $N\to\infty$.

This fact, together with Corollary 1 reduced to
$(t,s,t^{\prime},s^{\prime})=(t,s,t,s)$ implies the statement of Proposition
2.

###  B.6 Proof of Theorem 2

By Proposition 1 and Theorem 1, the model components $\mu,a,b,\sigma^{2}$ are
estimated consistently. Moreover, consider all the following statements
conditionally on $\mathbb{Y}^{new}$. Consequently,

$\widehat{\operatorname{Var}}(\mathbb{Y}^{new})\stackrel{{\scriptstyle
def}}{{=}}\left(\widehat{a}(t^{new}_{m},t^{new}_{m^{\prime}})\widehat{b}(s^{new}_{m},s^{new}_{m^{\prime}})+\widehat{\sigma}^{2}\mathds{1}_{[m=m^{\prime}]}\right)_{m,m^{\prime}=1}^{M^{new}}\stackrel{{\scriptstyle\mathbb{P}}}{{\to}}\operatorname{Var}(\mathbb{Y}^{new}),\qquad\text{as}\quad
N\to\infty,$

in the matrix space $\mathbb{R}^{M^{new}\times M^{new}}$. Due to continuity of
the matrix inversion and the fact that $\operatorname{Var}(\mathbb{Y}^{new})$
is positive definite,

$\left[\widehat{\operatorname{Var}}(\mathbb{Y}^{new})\right]^{-1}\stackrel{{\scriptstyle\mathbb{P}}}{{\to}}\left[\operatorname{Var}(\mathbb{Y}^{new})\right]^{-1}\qquad\text{as}\quad
N\to\infty.$

Moreover

$\widehat{\operatorname{Cov}}(X^{new}(t,s),\mathbb{Y}^{new})\stackrel{{\scriptstyle
def}}{{=}}\left(\widehat{a}(t,t^{new}_{m})\widehat{b}(s,s^{new}_{m})\right)_{m=1}^{M^{new}}=\operatorname{Cov}(X^{new}(t,s),\mathbb{Y}^{new})+o_{\mathbb{P}}\left(1\right),\qquad\text{as}\quad
N\to\infty,$

in the supremum norm over $(t,s)\in[0,1]^{2}$. Therefore, together with the
consistency of $\widehat{\mu}$ in the supremum norm, we conclude the statement
(3.4).

Assuming (A2), we conclude by the similar steps as above that

$\sup_{(t,s,t^{\prime},s^{\prime})\in[0,1]^{4}}\left|\widehat{\operatorname{Cov}}\left(X^{new}(t,s),X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)-\operatorname{Cov}\left(X^{new}(t,s),X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)\right|=o_{\mathbb{P}}\left(1\right),$
(B.29)

as $N\to\infty$.

Fixing $(t,s)\in[0,1]^{2}$ we have the conditional distribution given
$\mathbb{Y}^{new}$

$\frac{\Pi(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)}{\operatorname{Var}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}\sim
N\left(0,1\right)$

where the denominator is positive for all $t,s\in[0,1]$. Therefore

$\mathbb{P}\left(\left|\Pi(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)\right|\leq
u_{1-\alpha}\sqrt{\operatorname{Var}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}\,\middle|\,\mathbb{Y}^{new}\right)=1-\alpha.$

Now, since

$\frac{\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)}{\sqrt{\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}}\stackrel{{\scriptstyle
d}}{{\to}}N(0,1),\qquad\text{as}\quad N\to\infty.$

where $d$ denotes the convergence in distribution and therefore

$\lim_{N\to\infty}\mathbb{P}\left(\left|\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)\right|\leq
u_{1-\alpha}\sqrt{\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}\,\middle|\,\mathbb{Y}^{new}\right)=1-\alpha.$

It remains to justify the asymptotic coverage of the simultaneous confidence
band. By the constriction of the simultaneous confidence bands à la [9],
reviewed in Section 2.4, we have

$\mathbb{P}\left(\sup_{(t,s)\in[0,1]^{2}}\left|\Pi(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)\right|\leq
z_{1-\alpha}\sqrt{\operatorname{Var}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}\,\middle|\,\mathbb{Y}^{new}\right)=1-\alpha$

where the quantile $z_{1-\alpha}$ is calculated from the law of
$W=\sup_{(t,s)\in[0,1]^{2}}|Z(t,s)|$ where the true (non-estimated)
correlations are used:
$\operatorname{Cov}(Z(t,s),Z(t^{\prime},s^{\prime}))=\operatorname{Corr}\left(X^{new}(t,s),X^{new}(t^{\prime},s^{\prime})|\mathbb{Y}^{new}\right)$
with $t,t^{\prime},s,s^{\prime}\in[0,1]$. Recall that we denote the empirical
analogue of this law as $\widehat{W}$ already defined in (2.16).

In other words

$\sup_{(t,s)\in[0,1]^{2}}\left|\frac{\Pi(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)}{\sqrt{\operatorname{Var}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}}\right|\sim
W,\qquad\text{conditionally on}\,\mathbb{Y}^{new},$

and therefore

$\sup_{(t,s)\in[0,1]^{2}}\left|\frac{\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)}{\sqrt{\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}}\right|\stackrel{{\scriptstyle
d}}{{\to}}W,\qquad\text{as}\quad N\to\infty,\quad\text{conditionally
on}\,\mathbb{Y}^{new}.$

Now, if $c_{n}(\cdot,\cdot,\cdot,\cdot)\to c(\cdot,\cdot,\cdot,\cdot)$
uniformly (cf. (B.29)), then $N(0,c_{n})\stackrel{{\scriptstyle
d}}{{\to}}N(0,c)$. Therefore $\widehat{W}\stackrel{{\scriptstyle d}}{{\to}}W$
and thus $\widehat{z}_{1-\alpha}\to z_{1-\alpha}$ where
$\widehat{z}_{1-\alpha}$ and $z_{1-\alpha}$ are the quantiles calculated from
the law of $\widehat{W}$ and $W$ respectively. We conclude the proof by
observing

$\mathbb{P}\left(\sup_{(t,s)\in[0,1]^{2}}\frac{\left|\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)\right|}{\sqrt{\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}}\leq\widehat{z}_{1-\alpha}\,\middle|\,\mathbb{Y}^{new}\right)\\\
=\mathbb{P}\left(\sup_{(t,s)\in[0,1]^{2}}\frac{\left|\widehat{\Pi}(X^{new}(t,s)|\mathbb{Y}^{new})-X^{new}(t,s)\right|}{\sqrt{\widehat{\operatorname{Var}}\left(X^{new}(t,s)|\mathbb{Y}^{new}\right)}}\frac{z_{1-\alpha}}{\widehat{z}_{1-\alpha}}\leq
z_{1-\alpha}\,\middle|\,\mathbb{Y}^{new}\right)\to
1-\alpha,\qquad\text{as}\quad N\to\infty.$

## References

  * Aston, Pigoli and Tavakoli [2017] [author] Aston, John ADJ. A., Pigoli, DavideD. and Tavakoli, ShahinS. (2017). Tests for separability in nonparametric covariance operators of random surfaces. The Annals of Statistics 45 1431–1461.
  * Bagchi and Dette [2017] [author] Bagchi, PramitaP. and Dette, HolgerH. (2017). A test for separability in covariance operators of random surfaces. Annals of Statistics (to appear), available at arXiv:1710.08388.
  * Black and Scholes [1973] [author] Black, FischerF. and Scholes, MyronM. (1973). The pricing of options and corporate liabilities. Journal of political economy 81 637–654.
  * Chen et al. [2020] [author] Chen, YY., Carroll, CC., Dai, XX., Fan, JJ., Hadjipantelis, PZP., Han, KK., Ji, HH., Müller, HGH. and Wang, JLJ. (2020). fdapace: Functional Data Analysis and Empirical Dynamics.
  * Constantinou, Kokoszka and Reimherr [2017] [author] Constantinou, PanayiotisP., Kokoszka, PiotrP. and Reimherr, MatthewM. (2017). Testing separability of space-time functional processes. Biometrika 104 425–437.
  * Cont and Da Fonseca [2002] [author] Cont, RamaR. and Da Fonseca, JoséJ. (2002). Dynamics of implied volatility surfaces. Quantitative finance 2 45–60.
  * Cox, Ross and Rubinstein [1979] [author] Cox, John CJ. C., Ross, Stephen AS. A. and Rubinstein, MarkM. (1979). Option pricing: A simplified approach. Journal of financial Economics 7 229–263.
  * De Iaco, Myers and Posa [2002] [author] De Iaco, SS., Myers, Donald ED. E. and Posa, DD. (2002). Nonseparable space-time covariance models: some parametric families. Mathematical Geology 34 23–42.
  * Degras [2011] [author] Degras, David AD. A. (2011). Simultaneous confidence bands for nonparametric regression with functional data. Statistica Sinica 1735–1765.