Recent Developments in Complex and Spatially Correlated Functional Data

Israel Martínez-Hernández∗111 Statistics Program, King Abdullah University of
Science and Technology, Thuwal 23955-6900, Saudi Arabia.
∗ Corresponding author. E-mail<EMAIL_ADDRESS>
This publication is based on research supported by the King Abdullah
University of Science and Technology (KAUST) Office of Sponsored Research
(OSR) under Award No: OSR-2018-CRG7-3742.  and Marc G.
Genton${}^{\ref{note1}}$

Abstract

As high-dimensional and high-frequency data are being collected on a large
scale, the development of new statistical models is being pushed forward.
Functional data analysis provides the required statistical methods to deal
with large-scale and complex data by assuming that data are continuous
functions, e.g., realizations of a continuous process (curves) or continuous
random field (surfaces), and that each curve or surface is considered as a
single observation. Here, we provide an overview of functional data analysis
when data are complex and spatially correlated. We provide definitions and
estimators of the first and second moments of the corresponding functional
random variable. We present two main approaches: The first assumes that data
are realizations of a functional random field, i.e., each observation is a
curve with a spatial component. We call them spatial functional data. The
second approach assumes that data are continuous deterministic fields observed
over time. In this case, one observation is a surface or manifold, and we call
them surface time series. For these two approaches, we describe software
available for the statistical analysis. We also present a data illustration,
using a high-resolution wind speed simulated dataset, as an example of the two
approaches. The functional data approach offers a new paradigm of data
analysis, where the continuous processes or random fields are considered as a
single entity. We consider this approach to be very valuable in the context of
big data.

Keywords: Functional data; Functional random field; Manifold data; Spatial
functional data; Spatial statistics; Spatio-temporal statistics; Surface data.

Short title: Complex and Spatially Correlated Functional Data

##  1 Introduction

The statistical analysis of large, complex, and high-dimensional data has
become a significant challenging problem. Due to the rapid development of
complex, performant technologies, data can now be collected on a large scale,
resulting in high-dimensional and high-frequency data, sometimes necessitating
high-performance computing, which is often a limitation for practitioners; see
Galeano and Peña (2019) for a general view of data science and big data. Among
various approaches, functional data analysis (FDA) provides statistical
methods to handle large-scale and complex data (Chen et al., 2017; Giraldo et
al., 2018). For a general introduction to FDA, the reader is referred to
Ramsay and Silverman (2005), Ferraty and Vieu (2006), Horváth and Kokoszka
(2012), and Kokoszka and Reimherr (2017). FDA assumes that observations
(called functional data) have characteristics that vary along a continuum,
e.g., curves or surfaces. Thus, FDA deals with data that are defined on a
space that is intrinsically infinite-dimensional.

The approach of FDA has several advantages as a realization of a continuous
process or a continuous random field can be considered as functional data. In
this case, the stationarity of the process (random field) is not needed since
FDA treats the whole curve (surface) as a single entity. Thus, FDA is part of
object data analysis (Menafoglio and Secchi, 2017). FDA is useful when the
number of variables, $p$, is bigger than the sample size, $n$: $p\gg n$. In
particular, FDA can analyze longitudinal data. Smoothness is an important
property of functional data, in contrast with multivariate data analysis,
where smoothness has no meaning. Thus, FDA extracts additional information
contained in a continuous function or in its derivative. Although in practice,
each continuous functions, say $y_{i}(v)$, is observed on a finite set of
points, the continuity is obtained with smoothing techniques. In the process
of smoothing, FDA does not require observed $y_{i}(v)$ on a regular grid, that
is, a sample of $y_{i}(v)$ and $y_{j}(v)$ can be collected on a different set
of points; $\mathbf{v}_{i}=\\{v_{1_{i}},\ldots,v_{n_{i}}\\}$ and
$\mathbf{v}_{j}=\\{v_{1_{j}},\ldots,v_{n_{j}}\\}$, respectively. In general,
the methods of FDA are essentially nonparametric and can model complex and
spatially correlated data.

Functional data can also have a spatial component, because data are collected
somewhere at some time (Haining, 2003, p.15). If the functional data are
curves with a spatial component, we call them spatial functional data (SFD).
Thus, a dataset of SFD has the form
$y(\mathbf{s}_{1};v),\ldots,y(\mathbf{s}_{n};v)$, where $\mathbf{s}_{i}\in D$
represents the locations in a given region $D$, and $v$ represents the
continuous parameter of the functional data. For example,
$y(\mathbf{s}_{i};v)$ can be the daily wind speed observed at location
$\mathbf{s}_{i}$, with $v$ indicating the time within a day, or
$y(\mathbf{s}_{i};v)$ can be the spectrum of brain activity at location
$\mathbf{s}_{i}$, with $v$ representing the frequency. If the continuous
parameter $v$ represents time, as for the example of wind speed, then SFD can
be related to spatio-temporal data where the temporal dependence is captured
through the continuity of the curve.

The combination of FDA and spatial statistics provides a powerful tool to deal
with complex and large spatial data. This combination is attracting interest,
and much research is focused on this topic. Nerini et al. (2010) proposed a
spatial functional linear model, and they analyzed data in Oceanography. Zhou
et al. (2010) proposed mixed effects models for spatially correlated
hierarchical functional data. Ruiz-Medina (2011) extended the spatial
autoregressive processes and the spatial moving average processes to the
Hilbert space. Giraldo et al. (2012) proposed a methodology for clustering
spatially correlated functional data; see also, Jiang and Serban (2012) and
Romano et al. (2017). Staicu et al. (2010) proposed a methodology for
functional models with a hierarchical structure where the functions at the
lowest hierarchy level are spatially correlated. Delicado et al. (2010), Ruiz-
Medina (2012), and Mateu and Romano (2017) provided surveys of SFD. Menafoglio
and Secchi (2017) presented a review of complex and spatially dependent data,
such as curves and surfaces. Some references to SFD with a Bayesian
perspective are Baladandayuthapani et al. (2008), Zhang et al. (2016), Song
and Mallick (2019), and Rekabdarkolaee et al. (2019).

Functional data can also have a complex domain, e.g., a two-dimensional
Euclidean domain. Spatial data can be considered as functional data with the
same domain as the corresponding random field, for instance if observations
are dense over the region or if the domain is not a subset of the Euclidean
space; see, e.g., Alfeld et al. (1996), Wahba (1981), and Gneiting (2013) for
spatial data over complex domains. With the FDA approach, a realization of a
random field is considered as a single point observation of the functional
data. Similarly to the case when data are curves, the continuous surface (or
the manifold) needs to be estimated. For that estimation, one can use a tensor
product of univariate B-splines (Eilers and Marx, 1996; Wood, 2006; Qingguo
and Longsheng, 2010; Xiao et al., 2013). Another way is to approximate the
continuous data using finite elements analysis; see Ramsay (2002), Duchamp and
Stuetzle (2003), and Sangalli et al. (2013). These types of functional data
can also be observed over time; they are then called a functional time series.
Functional time series can be related to spatio-temporal data when
observations are dense over the domain. Thus, the approach of functional time
series can be used in complex and large spatio-temporal data.

In this paper, our goal is to provide a review of complex and spatially
correlated functional data, with two approaches using either spatial
functional data or surface (manifold) data. In both approaches, we focus on
covariance functions and modeling. The remainder of our paper is organized as
follows: In Section 2, we present basic concepts of functional random fields
that include the mean, the covariance, and other important concepts related to
the covariance. In this section, we also present different estimators of the
various objects defined. In Section 3, we describe how to model small-scale
and large-scale variations of spatial functional data, as well as
corresponding methods to estimate the parameters involved. Section 4 (first
approach) presents a brief description of the concept of functional kriging.
In Section 5, we describe the second approach based on surface data, which can
be considered as an alternative to the analysis of spatio-temporal data, and
describe how to model and estimate the continuous surface. In Section 6, we
present some software available for the analysis of spatial functional data
and surface (manifold) data. In Section 7, we provide an example of the two
approaches using a high-resolution wind speed simulated dataset in Dumat Al
Jandal, Saudi Arabia. Section 8 ends the paper with some discussions.

##  2 Functional Random Fields

###  2.1 Basic Concepts

In this section, we introduce the basic concepts of SFD. Let $(\Omega,F,P)$ be
a probability space. Without loss of generality, we assume that the domain
$\mathcal{T}$ of the curves is $\mathcal{T}=[0,1]$, and let
$\mathcal{H}=L^{2}([0,1])$ be the Hilbert space of square integrable functions
defined on $[0,1]$, equipped with the inner product $\langle
f,g\rangle=\int_{0}^{1}\\!f(v)g(v)\mathrm{d}v$. We denote by
$\|\cdot\|_{\mathcal{H}}$ the norm in $\mathcal{H}$ induced by the inner
product. A random variable $X:\Omega\to\mathcal{H}$ taking values in the
Hilbert space $\mathcal{H}$ is called a functional random variable (Ramsay and
Silverman, 2005; Ferraty and Vieu, 2006; Horváth and Kokoszka, 2012). Let
$D\subset\mathbb{R}^{2}$ be a fixed study area (e.g., a country). A random
field $\\{X(\mathbf{s}):\mathbf{s}\in D\\}$ taking values in $\mathcal{H}$ is
called a functional random field, that is, for each location $\mathbf{s}\in
D$, $X(\mathbf{s}):\Omega\to\mathcal{H}$ is a functional random variable.

We denote by $X(\mathbf{s}_{0};v)$ the functional random variable at a fixed
location $\mathbf{s}_{0}$ and $v\in[0,1]$, and we denote by
$X(\mathbf{s}_{0};v_{0})$ the scalar random variable obtained by evaluating
$X(\mathbf{s};v)$ at $\mathbf{s}=\mathbf{s}_{0}$ and $v=v_{0}$. Lastly, we use
$f$ to denote a function in $\mathcal{H}$.

Let $X(\mathbf{s};v)$ be a functional random field with
$\mathbb{E}(\|X(\mathbf{s};v)\|_{\mathcal{H}})<\infty$, for all $\mathbf{s}\in
D$. The mean $\mu(\mathbf{s};v):=\mathbb{E}\\{X(\mathbf{s};v)\\}$ of
$X(\mathbf{s};v)$ is defined as an element of $\mathcal{H}$ such that

$\langle\mu(\mathbf{s};\cdot),f\rangle=\mathbb{E}\langle
X(\mathbf{s};\cdot),f\rangle,\quad\forall f\in\mathcal{H},$

where the dot in $(\mathbf{s};\cdot)$ indicates the integrated variable in
$[0,1]$. This implies that
$\mathbb{E}\\{X(\mathbf{s};v_{0})\\}=\mu(\mathbf{s};v_{0})$ for almost all
$v_{0}\in[0,1]$.

The covariance is one of the most studied objects in spatial statistics, due
to its relevance for prediction. If
$\mathbb{E}\\{\|X(\mathbf{s};v)\|^{2}_{\mathcal{H}}\\}<\infty$, then the
covariance operator at locations $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ is
defined as an operator
$C(\mathbf{s}_{1},\mathbf{s}_{2};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}}):\mathcal{H}\to\mathcal{H}$
such that

$\displaystyle C(\mathbf{s}_{1},\mathbf{s}_{2};f)(\cdot)$
$\displaystyle=\mathbb{E}[\langle
X(\mathbf{s}_{1};\cdot)-\mu(\mathbf{s}_{1};\cdot),f\rangle\\{X(\mathbf{s}_{2};\cdot)-\mu(\mathbf{s}_{2};\cdot)\\}]$
$\displaystyle=\int_{0}^{1}\sigma(\mathbf{s}_{1},\mathbf{s}_{2};u,\cdot)f(u)\mathrm{d}u,\quad
f\in\mathcal{H},$ (2.1)

where
$\sigma(\mathbf{s}_{1},\mathbf{s}_{2};u_{0},v_{0}):=\mathbb{E}[\\{X(\mathbf{s}_{1};u_{0})-\mu(\mathbf{s}_{1};u_{0})\\}\\{X(\mathbf{s}_{2};v_{0})-\mu(\mathbf{s}_{2};v_{0})\\}]$
is the point-wise covariance, called the kernel of
$C(\mathbf{s}_{1},\mathbf{s}_{2};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$.
This definition can be written in terms of a tensor operation as follows:
$C(\mathbf{s}_{1},\mathbf{s}_{2};f)=\mathbb{E}[\\{X(\mathbf{s}_{1};\cdot)-\mu(\mathbf{s}_{1};\cdot)\\}\otimes\\{X(\mathbf{s}_{2};\cdot)-\mu(\mathbf{s}_{2};\cdot)\\}(f)]$.

A common assumption, in practice, is the stationarity condition of a process,
which is defined as follows:

######  Definition 1 (Weak stationarity)

A functional random field $X(\mathbf{s};v)$ is said to be (weakly) stationary
if

  1. 1.

$\mathbb{E}(\|X(\mathbf{s};v)\|_{\mathcal{H}}^{2})<\infty$,

  2. 2.

$\mu(\mathbf{s};v)=\mu(v)$, that is, the mean does not depend on the location
$\mathbf{s}$, and

  3. 3.

$C(\mathbf{s}_{1}+\mathbf{h},\mathbf{s}_{2}+\mathbf{h};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})=C(\mathbf{s}_{1},\mathbf{s}_{2};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$,
for all $\mathbf{s}_{1},\mathbf{s}_{2},\mathbf{h}\in D$.

The last condition is equivalent to the property that the covariance operator
depends only on the increments $\mathbf{s}_{1}-\mathbf{s}_{2}$. This means,
there exists a covariance operator
$\tilde{C}(\mathbf{s};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}}):\mathcal{H}\to\mathcal{H}$
such that

$C(\mathbf{s}_{1},\mathbf{s}_{2};f)=\tilde{C}(\mathbf{s}_{1}-\mathbf{s}_{2};f),\quad
f\in\mathcal{H},$

and so the variance operator can be written as
$\mathrm{Var}\\{X(\mathbf{s};v)\\}=C(\mathbf{s},\mathbf{s};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})=\tilde{C}(\mathbf{0};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$.
Thus, for convenience, we write
$C(\mathbf{0};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$
to denote the variance operator $C(\mathbf{s},\mathbf{s};f)$ of the stationary
functional random field $X(\mathbf{s};v)$.

Now, we define the concept of isotropy for SFD.

######  Definition 2 (Isotropy)

A stationary functional random field $X(\mathbf{s};v)$ is said to be isotropic
if there exists a covariance operator
$\tilde{C}_{0}(h;\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}}):\mathcal{H}\to\mathcal{H}$
such that

$C(\mathbf{s}_{1},\mathbf{s}_{2};f)=\tilde{C}_{0}(h;f),\quad f\in\mathcal{H},$

where $h=\|\mathbf{s}_{1}-\mathbf{s}_{2}\|$, for all
$\mathbf{s}_{1},\mathbf{s}_{2}\in D$.

In spatial statistics, the variogram plays an important role to make
inference. The extension of this concept to functional random fields is as
follows: The variogram operator $\Gamma$ is defined as the variance operator
of the difference between the functional random field at two locations
$\mathbf{s}_{1},\mathbf{s}_{2}$, i.e.,
$\Gamma(\mathbf{s}_{1},\mathbf{s}_{2};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}}):=\frac{1}{2}\mathrm{Var}\\{X(\mathbf{s}_{1};v)-X(\mathbf{s}_{2};v)\\}$.
If the functional random field $X(\mathbf{s};v)$ has a finite second moment,
then we have
$2\Gamma(\mathbf{s}_{1},\mathbf{s}_{2};f)=C(\mathbf{s}_{1},\mathbf{s}_{1};f)+C(\mathbf{s}_{2},\mathbf{s}_{2};f)-C(\mathbf{s}_{1},\mathbf{s}_{2};f)-C(\mathbf{s}_{2},\mathbf{s}_{1};f),$
for $f\in\mathcal{H}$. Thus, if the functional random field is stationary,
then there exists a variogram operator
$\tilde{\Gamma}(\mathbf{s};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}}):\mathcal{H}\to\mathcal{H}$
such that it satisfies

$\Gamma(\mathbf{s}_{1},\mathbf{s}_{2};f)=\tilde{\Gamma}(\mathbf{s}_{1}-\mathbf{s}_{2};f)=C(\mathbf{0};f)-C(\mathbf{s}_{1}-\mathbf{s}_{2};f),\quad
f\in\mathcal{H},$ (2.2)

where the corresponding kernel is
$\gamma(\mathbf{s}_{1}-\mathbf{s}_{2};u,v)=\sigma(\mathbf{0};u,v)-\sigma(\mathbf{s}_{1}-\mathbf{s}_{2};u,v)$.
Furthermore, if $X(\mathbf{s};v)$ is isotropic, then there exists
$\tilde{\Gamma}_{0}(h;\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}}):\mathcal{H}\to\mathcal{H}$
such that

$\Gamma(\mathbf{s}_{1},\mathbf{s}_{2};f)=\tilde{\Gamma}_{0}(h;f),\quad
f\in\mathcal{H},$

where $h=\|\mathbf{s}_{1}-\mathbf{s}_{2}\|$, for all
$\mathbf{s}_{1},\mathbf{s}_{2}\in D$.

###### Remark 1

Unlike the finite dimensional case (multivariate random field) where the
covariance is composed of matrices, here, the covariance is composed of
operators, since the space of functional data is intrinsically infinite-
dimensional.

Another way to describe the second-order spatial dependence of the functional
random field is by using a “global” measure. This global measure is the trace-
covariogram $\sigma_{\mathrm{tr}}:D\times D\to\mathbb{R}$ (Giraldo et al.,
2011; Menafoglio et al., 2013) defined as

$\sigma_{\mathrm{tr}}(\mathbf{s}_{1},\mathbf{s}_{2})=\mathbb{E}\\{\langle
X(\mathbf{s}_{1};\cdot)-\mu(\mathbf{s}_{1};\cdot),X(\mathbf{s}_{2};\cdot)-\mu(\mathbf{s}_{2};\cdot)\rangle\\}=\int_{0}^{1}\\!\sigma(\mathbf{s}_{1},\mathbf{s}_{2};v,v)\mathrm{d}v.$
(2.3)

The trace-covariogram computes the covariance of the inner product of the
functional random field at two locations. Thus, it summarizes the covariance
on the diagonal, and so, in general, it depends only on the locations. If the
functional random field is stationary, then there exists
$\tilde{\sigma}_{\mathrm{tr}}:D\to\mathbb{R}$ such that it depends only on the
separation vector $\mathbf{s}_{1}-\mathbf{s}_{2}$, that is,

$\sigma_{\mathrm{tr}}(\mathbf{s}_{1},\mathbf{s}_{2})=\tilde{\sigma}_{\mathrm{tr}}(\mathbf{s}_{1}-\mathbf{s}_{2}).$

In addition, if the functional random field is isotropic, then there exists
$\tilde{\sigma}_{\mathrm{tr},0}:\mathbb{R}\to\mathbb{R}$ such that

$\sigma_{\mathrm{tr}}(\mathbf{s}_{1},\mathbf{s}_{2})=\tilde{\sigma}_{\mathrm{tr},0}(h),$

where $h=\|\mathbf{s}_{1}-\mathbf{s}_{2}\|$.

Similarly, the trace-variogram is defined in terms of the inner product of the
difference, i.e.,

$\gamma_{\mathrm{tr}}(\mathbf{s}_{1}-\mathbf{s}_{2})=\frac{1}{2}\mathbb{E}\left\\{\langle
X(\mathbf{s}_{1};\cdot)-X(\mathbf{s}_{2};\cdot),X(\mathbf{s}_{1};\cdot)-X(\mathbf{s}_{2};\cdot)\rangle\right\\}\\\
-\frac{1}{2}\|\mu(\mathbf{s}_{1};\cdot)-\mu(\mathbf{s}_{2};\cdot)\|_{\mathcal{H}}^{2}.$

We observe that, if $X(\mathbf{s};v)$ is stationary, then
$\sigma(\mathbf{s};u,v)=\sigma(\mathbf{0};u,v)-\gamma(\mathbf{s};u,v).$ Thus,
the trace-variogram satisfies

$\gamma_{\mathrm{tr}}(\mathbf{s}_{1}-\mathbf{s}_{2})=\sigma_{\mathrm{tr}}(\mathbf{0})-\sigma_{\mathrm{tr}}(\mathbf{s}_{1}-\mathbf{s}_{2}).$

Trace-covariogram and trace-variogram are also important for optimization
problems. Especially if we want to use the criterion of minimizing equations
of the form $\mathbb{E}(\langle X,Y\rangle)$, as in (2.7) below.

###  2.2 Estimation

Now, we describe estimators of the mean $\mu$, the covariance $C$, and the
variogram $\Gamma$. For this purpose, we assume that $X(\mathbf{s};v)$ is an
isotropic (stationary) functional random field with mean $\mu(v)$ and
covariance operator $C$. Let $x(\mathbf{s}_{1};v),\ldots,x(\mathbf{s}_{n};v)$
be observations of the functional random field $X(\mathbf{s};v)$. We assume
that the observations $x(\mathbf{s}_{i};v)$ are given in the functional form.
Although in real data, $x(\mathbf{s}_{i};v)$ are observed on a finite set of
points $v_{i1},\ldots,v_{im}$, the continuous curves should be estimated (see
Ramsay and Silverman, 2005).

The main feature of spatial data is that “nearby” data look similar, and an
estimator must take into account such spatial dependence. Otherwise, it will
not have desirable properties, such as consistency.

####  2.2.1 Mean estimation

We describe two different approaches to obtain an estimator of the mean
$\mu(v)$ (Gromenko et al., 2012). A model of the mean can be written as

$X(\mathbf{s};v)=\mu(v)+\varepsilon(\mathbf{s};v),$ (2.4)

where $\varepsilon(\mathbf{s};v)$ is an isotropic functional random field with
zero mean and covariance operator $C$.

The first approach is similar to the kriging method. Specifically, this is
defined as a weighting of the observed curves:

$\hat{\mu}(v)=\sum_{i=1}^{n}w_{i}x(\mathbf{s}_{i};v),$ (2.5)

where the weights $w_{i}$ are estimated by solving the optimization problem

$\min_{w_{1},\ldots,w_{n}}\mathbb{E}\\{\langle\hat{\mu}-\mu,\hat{\mu}-\mu\rangle\\}=\min_{w_{1},\ldots,w_{n}}\mathbb{E}\left\\{\left\|\sum_{i=1}^{n}w_{i}x(\mathbf{s}_{i};\cdot)-\mu\right\|^{2}_{\mathcal{H}}\right\\},$
(2.6)

subject to the condition $\sum_{i=1}^{n}w_{i}=1$. Using the Lagrange
multiplier method, this leads to solve

$\sum_{i=1}^{n}w_{i}=1,\quad\sum_{i=1}^{n}w_{i}\,\sigma_{\mathrm{tr},\varepsilon}(\mathbf{s}_{i},\mathbf{s}_{j})-\lambda=0,\quad
j=1,\ldots,n,$ (2.7)

where $\sigma_{\mathrm{tr},\varepsilon}(\mathbf{s}_{1},\mathbf{s}_{2})$ is the
trace-covariogram of $\varepsilon(\mathbf{s};v)$. Thus, the estimation problem
(2.6) becomes estimating the matrix
$\\{\sigma_{\mathrm{tr},\varepsilon}(\mathbf{s}_{i},\mathbf{s}_{j})\\}_{i,j=1}^{n}$.
Since the functional random field $\varepsilon(\mathbf{s};v)$ is unobserved, a
common approach is to use an iterative procedure. At the first iteration, an
initial estimator of $\mu(v)$ is obtained by assuming that
$\varepsilon(\mathbf{s};v)$ is spatially uncorrelated, i.e.,
$\hat{\mu}_{0}(v)=\frac{1}{n}\sum_{i=1}^{n}x(\mathbf{s}_{i};v)$. Next,
$\hat{\mu}_{0}(v)$ is subtracted from the data $x(\mathbf{s}_{i};v)$, then an
initial estimator of
$\sigma_{\mathrm{tr},\varepsilon}(\mathbf{s}_{i},\mathbf{s}_{j})$ is obtained
as described below in (2.11). At the second iteration, the mean is re-
estimated by solving (2.7), with the initial information of
$\sigma_{\mathrm{tr},\varepsilon}(\mathbf{s}_{i},\mathbf{s}_{j})$. This
process is repeated until convergence.

The second approach uses finite basis functions similarly to the cokriging
method on the coefficients of the basis functions, see Goulard and Voltz
(1993), Nerini et al. (2010), and Giraldo et al. (2011). Let
$\\{\eta_{1}(v),\ldots,\eta_{K}(v)\\}$ be basis functions, e.g., Fourier basis
functions or B-spline basis functions. Then, the observed curves are
approximated as

$x(\mathbf{s}_{i};v)\approx\sum_{k=1}^{K}z_{k}(\mathbf{s}_{i})\eta_{k}(v),\quad
i=1,\ldots,n,$ (2.8)

where $z_{k}(\mathbf{s}_{i})=\langle x(\mathbf{s}_{i};\cdot),\eta_{k}\rangle$,
which is a scalar for each $i$ and $k$. Let $Z_{k}(\mathbf{s})$ be the
corresponding scalar random field with realization
$\\{z_{k}(\mathbf{s}_{i})\\}_{i=1}^{n}$. Then, by using (2.4), the mean can be
approximated as

$\mu(v)\approx\sum_{k=1}^{K}\mathbb{E}\\{Z_{k}(\mathbf{s})\\}\eta_{k}(v),$

and $\mathbb{E}\\{Z_{k}(\mathbf{s})\\}$ should be estimated. For this, from
(2.4) we have that

$Z_{k}(\mathbf{s}_{i})=\mu^{z}_{k}+\varepsilon^{k}(\mathbf{s}_{i}),\quad
i=1,\ldots,n,$

where
$\mu_{k}^{z}=\langle\mu,\eta_{k}\rangle=\mathbb{E}\\{Z_{k}(\mathbf{s})\\}$,
and
$\varepsilon^{k}(\mathbf{s}_{i})=\langle\varepsilon(\mathbf{s}_{i};\cdot),\eta_{k}\rangle$,
for each $k=1,\ldots,K$ . Thus, $Z_{k}(\mathbf{s}_{i})$ is an isotropic
(stationary) scalar random field, for $k=1,\ldots,K$. Once $\mu^{z}_{k}$ is
estimated, we obtain the estimator of the mean

$\hat{\mu}(v)=\sum_{k=1}^{K}\hat{\mu}^{z}_{k}\eta_{k}(v).$

To estimate $\mu^{z}_{k}$, we can use the ordinary least squares (OLS)
estimator, by minimizing

$\sum_{i=1}^{n}\left\\{z_{k}(\mathbf{s}_{i})-w_{i}\mu^{z}_{k}\right\\}^{2},$

with $w_{i}=1$ for all $i=1,\ldots,n$, or a weighted least squares (WLS)
estimator, i.e, $w_{i}$ are estimated with the information of the covariance
matrix of $\varepsilon^{k}(\mathbf{s})$. Also, we can use a generalized least
squares (GLS) estimator, which can perform better than OLS or WLS; see
Schabenberger and Gotway (2005) and Cressie (2015) for the development of
these estimators.

Another way to estimate $\mu^{z}_{k}$ is by using Maximum Likelihood
Estimators (MLEs). For this, assume that, for each $i$, $X(\mathbf{s}_{i};v)$
is a Gaussian process in $[0,1]$. Then, for each $k$,
$\mathbf{z}_{k}=\\{z_{k}(\mathbf{s}_{1}),\ldots,z_{k}(\mathbf{s}_{n})\\}^{\mathsf{T}}$
is a realization of $\mathbf{Z}_{k}\sim{\cal
N}_{n}(\mu^{z}_{k}\mathbf{1},\bm{\Sigma}(\bm{\theta})),$ where
$\bm{\Sigma}(\bm{\theta})$ is the covariance of $\varepsilon^{k}(\mathbf{s})$
described by some valid parametric covariance model, e.g., Matérn,
exponential, or spherical. Then, the log likelihood is

$l(\bm{\theta},\mu^{z}_{k})=-\frac{n}{2}\log(2\pi)-\frac{n}{2}\log\,\mathrm{det}\\{\bm{\Sigma}(\bm{\theta})\\}-\frac{1}{2}(\mathbf{z}_{k}-\mu^{z}_{k}\mathbf{1})^{\mathsf{T}}\bm{\Sigma}(\bm{\theta})^{-1}(\mathbf{z}_{k}-\mu^{z}_{k}\mathbf{1}).$
(2.9)

To read more details of this estimator, see Stein (1999).

The coefficients in estimator (2.5) can be also considered as operators
$\omega_{i}:\mathcal{H}\to\mathcal{H}$. In this case the estimator of the mean
takes the form
$\hat{\mu}(v)=\sum_{i=1}^{n}\omega_{i}\\{x(\mathbf{s}_{i};\cdot)\\}(v)$, where
the coefficients, which are operators, need to be estimated under a certain
constraint on $\sum_{i=1}^{n}\omega_{i}$, see Nerini et al. (2010). This
approach takes into account the information of the whole curve to define a
specific weight for each point $v\in[0,1]$. Thus, it is expected to obtain
better results than estimator (2.5), but the estimation procedure can be
complicated.

In the context of a nonstationary functional random field, the mean can depend
on the location $\mathbf{s}$, in this case, $\mu(\mathbf{s};v)$ can be
represented by a linear model
$\mu(\mathbf{s};v)=\sum_{l=0}^{L}a_{l}(\mathbf{s})f_{l}(v)$, where $f_{l}(v)$,
$l=0,\ldots,L$, are elements of $\mathcal{H}$ independent of the spatial
location $\mathbf{s}$, and $a_{0}(\mathbf{s})=1$ (Menafoglio et al., 2013;
Caballero et al., 2013). The latter approach is described with more details in
Section 3.1.

####  2.2.2 Covariance estimation

Now, we describe how to estimate the covariance. We assume that
$\mathbb{E}\\{X(\mathbf{s};v)\\}=0$. Let $N(h)$ be a set of pairs of indexes
defined as $N(h)=\\{(i,j):\|\mathbf{s}_{i}-\mathbf{s}_{j}\|=h\\}$, and let
$|N(h)|$ be the cardinality of $N(h)$. The empirical covariance operator of
the functional random field is defined as

$\hat{C}(h;f)=\frac{1}{|N(h)|}\sum_{(i,j)\in N(h)}x(\mathbf{s}_{i})\otimes
x(\mathbf{s}_{j})(f).$

In practice, the distance between $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ is not
considered to be exactly $h$, instead
$\|\mathbf{s}_{i}-\mathbf{s}_{j}\|\in(h-\delta,h+\delta)$, with $\delta>0$.
Also, since it is almost impossible to obtain $\hat{C}(h;f)$ for all $h$,
discretized values of $h$, $h_{1},\ldots,h_{m}$, are computed.

The usual approach in the scalar (multivariate) random field case is to fit a
valid parametric model to the empirical covariance by least squares methods
(Cressie, 2015). However, in spatial functional data, for each $h$,
$\hat{C}(h;\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$
is an operator. Thus, it requires new mathematical developments to define a
“covariance model” in this context. Finite basis functions have been used to
overcome the modeling of the covariance function in $\mathcal{H}$ (Nerini et
al., 2010). Recall that
$C(h;\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$
is represented by its kernel $\sigma(h;u,v)$ (see (2.1)). Now, we assume that
$x(\mathbf{s}_{i};v)$ is approximated by a finite set of basis functions, as
in (2.8). Then, we obtain

$\sigma(h;u,v)=\mathbb{E}\\{X(\mathbf{s+h})X(\mathbf{s})\\}=\bm{\eta}^{\mathsf{T}}(u)\mathbb{E}\\{\mathbf{Z}(\mathbf{s+h})\mathbf{Z}^{\mathsf{T}}(\mathbf{s})\\}\bm{\eta}(v),$

where $\bm{\eta}(v)=\\{\eta_{1}(v),\ldots,\eta_{K}(v)\\}^{\mathsf{T}}$,
$\mathbf{Z}(\mathbf{s})=\\{Z_{1}(\mathbf{s}),\ldots,Z_{K}(\mathbf{s})\\}^{\mathsf{T}}$,
and $Z_{k}(\mathbf{s})=\langle X(\mathbf{s};\cdot),\eta_{k}\rangle$ is the
scalar random field with mean zero, and so $\mathbf{Z}(\mathbf{s})$ is a
multivariate random field with covariance

$\bm{\Sigma}(h)=\mathrm{Cov}\\{\mathbf{Z}(\mathbf{s+h}),\mathbf{Z}(\mathbf{s})\\}=\\{\bm{\Sigma}_{kl}(h)\\}_{k,l=1}^{K},$

where
$\bm{\Sigma}_{kl}(h)=\mathrm{Cov}\\{Z_{k}(\mathbf{s+h}),Z_{l}(\mathbf{s})\\}$
for $k,l=1,\ldots,K$, $h=\|\mathbf{h}\|$. Thus, estimating the covariance
operator $C$ is equivalent to estimating the covariance of a multivariate
random field $\mathbf{Z}(\mathbf{s})$, which can be modeled with a valid
covariance function $\bm{\Sigma}(h;\bm{\theta})$ (Genton and Kleiber, 2015).

If $\bm{\eta}(v)$ is a set of orthogonal basis functions, then we can estimate
each marginal-covariance functions $\bm{\Sigma}_{kk}(h)$, for $k=1,\ldots,K$,
separately. In addition, if $X(\mathbf{s};v)$ is assumed to be a Gaussian
process, we can use the log likelihood (2.9) with $\mu^{z}_{k}=0$, i.e.,
maximizing $l(\bm{\theta},0)$ as function of $\bm{\theta}$. Also, the
restricted maximum likelihood (REML) estimation can be used in this case, see
Stein (1999).

Once the estimator $\hat{\bm{\Sigma}}(h;\bm{\theta})$ is obtained, then we
define an estimator of $\sigma(h;u,v)$ as follows,

$\hat{\sigma}(h;u,v)=\bm{\eta}^{\mathsf{T}}(u)\hat{\bm{\Sigma}}(h;\bm{\theta})\bm{\eta}(v)$
(2.10)

and then,
$\hat{C}(h;f)(\cdot)=\int_{0}^{1}\hat{\sigma}(h;u,\cdot)f(u)\mathrm{d}u$.

In the estimation of the covariance of $Z_{k}$, we can also use a
nonparametric approach (Hall et al., 1994; Hall and Patil, 1994). Bayesian
approaches can be found in, e.g., Banerjee et al. (2015) and Diggle and
Ribeiro (2007).

As mentioned before, the trace-covariogram is a measure that describes
dependence globally, in the sense that it integrates the kernel
$\sigma(\mathbf{s};u,v)$ on the diagonal. Also, the trace-covariogram appears
in several optimization problems (Delicado et al., 2010; Giraldo et al., 2011;
Menafoglio et al., 2013), that make its estimation important. The empirical
trace-covariogram is defined as

$\hat{\sigma}_{\mathrm{tr}}(h)=\frac{1}{|N(h)|}\sum_{(\mathbf{s}_{i},\mathbf{s}_{j})\in
N(h)}\int_{0}^{1}x(\mathbf{s}_{i};v)x(\mathbf{s}_{j};v)\mathrm{d}v.$ (2.11)

This empirical estimator is computed on discrete values $h_{1},\ldots,h_{m}$
of $h$. Then, we can fit any covariance model by least squares methods to the
empirical estimates
$\hat{\sigma}_{\mathrm{tr}}(h_{1}),\ldots,\hat{\sigma}_{\mathrm{tr}}(h_{m})$.

We observe that estimating the trace-covariogram is much easier than
estimating the covariance $C$, where the inner products $\langle
x(\mathbf{s}_{i};\cdot),x(\mathbf{s}_{j};\cdot)\rangle=\int_{0}^{1}x(\mathbf{s}_{i};v)x(\mathbf{s}_{j};v)\mathrm{d}v$
can be computed using the R (R Core Team, 2019) package fda (Ramsay et al.,
2018).

####  2.2.3 Variogram estimation

One of the advantages of using the variogram is the robustness under
misspecification of the mean. Also, the variogram can be used to estimate the
covariance. The empirical variogram operator of the functional random field is
defined as

$\hat{\Gamma}(h;f)=\frac{1}{2|N(h)|}\sum_{(i,j)\in
N(h)}\\{x(\mathbf{s}_{i})-x(\mathbf{s}_{j})\\}\otimes\\{x(\mathbf{s}_{i})-x(\mathbf{s}_{j})\\}(f),\quad
f\in\mathcal{H}.$

Similarly to the estimator of the covariance operator, if a basis function
$\bm{\eta}(v)$ is assumed, then the corresponding kernel $\gamma$ is obtained
as

$\displaystyle 2\gamma(h;u,v)$
$\displaystyle=\mathbb{E}[\\{X(\mathbf{s+h})-X(\mathbf{s})\\}^{2}]$
$\displaystyle=\bm{\eta}^{\mathsf{T}}(u)\mathbb{E}[\\{\mathbf{Z}(\mathbf{s}+\mathbf{h})-\mathbf{Z}(\mathbf{s})\\}\\{\mathbf{Z}(\mathbf{s}+\mathbf{h})-\mathbf{Z}(\mathbf{s})\\}^{\mathsf{T}}]\bm{\eta}(v),$

where $h=\|\mathbf{h}\|$. Thus, we need to estimate the variogram of the
multivariate random field
$\mathbf{Z}(\mathbf{s})=\\{Z_{1}(\mathbf{s}),\ldots,Z_{K}(\mathbf{s})\\}^{\mathsf{T}}$,
to obtain the estimator of the kernel $\gamma(h;u,v)$.

The corresponding empirical trace-variogram is defined as

$\hat{\gamma}_{\mathrm{tr}}(h)=\frac{1}{2|N(h)|}\sum_{(i,j)\in
N(h)}\int_{0}^{1}\\{x(\mathbf{s}_{i};v)-x(\mathbf{s}_{j};v)\\}^{2}\mathrm{d}v.$
(2.12)

Then, following the common method used in spatial statistics, a variogram
model is fitted to
$\hat{\gamma}_{\mathrm{tr}}(h_{1}),\ldots,\hat{\gamma}_{\mathrm{tr}}(h_{m})$.

The principal component analysis is generally important in statistics because
of its applicability to dimensional reduction techniques. In the context of
functional data, a function $\zeta\in\mathcal{H}$ is called eigenfunction of
the operator
$C(\mathbf{s},\mathbf{s};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$,
if $C(\mathbf{s},\mathbf{s};\zeta)=\lambda\zeta$ with $\lambda$ a positive
real number. The estimation of eigenfunctions can be obtained from the
estimated covariance operator
$\hat{C}(\mathbf{s},\mathbf{s};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$.
If we are only interested in the eigenfunctions, we can obtain the estimators
without estimating the covariance operator, which allows us to reduce
computational costs (Zhou et al., 2010; Gromenko et al., 2012). Also see Liu
et al. (2017) for functional principal component analysis of spatial
functional data.

##  3 Modeling Functional Random Fields

Spatial statistics studies the variations among the observed data at different
locations. The spatial variation is generally described through the mean and
the covariance (Haining, 2003). The mean represents the large-scale
variations, and the covariance represents the small-scale variations.

In this section, we describe the statistical models for spatial functional
data. We denote the observed functional data as
$y(\mathbf{s}_{1};v),\ldots,y(\mathbf{s}_{n};v)$, and we denote by
$Y(\mathbf{s};v)$ the corresponding functional random field with realization
$y(\mathbf{s}_{i};v)$, $i=1,\ldots,n$. The functional random field
$X(\mathbf{s};v)$ denotes a stationary functional random field with covariance
function $C_{X}$, and $\varepsilon(\mathbf{s};v)$ denotes a functional white
noise, i.e., $\mathbb{E}\\{\varepsilon(\mathbf{s},v)\\}=0$, with a covariance
such that
$C_{\varepsilon}(\mathbf{s}_{1},\mathbf{s}_{2};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})=0$,
if $\mathbf{s}_{1}\neq\mathbf{s}_{2}$.

###  3.1 Large-scale Variation

Regression models with covariates in terms of the locations can be used to
model large-scale spatial variations. Such regression models must account for
spatial dependence. In general, estimators obtained from regression models are
smooth functions defined in $D$. In the context of spatial functional data,
functional regression models in $\mathcal{H}$ extend models of finite-
dimensional data to model large-scale variations. In Caballero et al. (2013),
Menafoglio et al. (2013), and Reyes et al. (2015), covariates are assumed to
be separable in the spatial component and the continuity of the data. In this
case, the model for large-scale variations is

$Y(\mathbf{s};v)=\sum_{l=0}^{L}a_{l}(\mathbf{s})f_{l}(v)+\varepsilon(\mathbf{s};v),$
(3.1)

where $f_{l}\in\mathcal{H}$ are independent of $\mathbf{s}$,
$a_{0}(\mathbf{s}):=1$, $\\{a_{l}(\mathbf{s})\\}_{l=1}^{L}$ are known scalar
regressors, and $\varepsilon(\mathbf{s};v)$ is the functional white noise. For
example, with $L=5$, we could specify the scalar regressors
$a_{l}(\mathbf{s})$ as

$a_{1}(\mathbf{s})=s_{1},\,a_{2}(\mathbf{s})=s_{2},\,a_{3}(\mathbf{s})=s_{1}s_{2},\,a_{4}(\mathbf{s})=s_{1}^{2},\mbox{
and }a_{5}(\mathbf{s})=s_{2}^{2},$

where $\mathbf{s}=(s_{1},s_{2})$ denotes the coordinates of a spatial
location. In this example, $f_{0}(v),\ldots,f_{5}(v)$ need to be estimated.
Thus, model (3.1) is a functional regression model with scalar covariates
$a_{l}(\mathbf{s})$, $l=0,1,\ldots,L$.

The mean of the functional random field $Y(\mathbf{s};v)$ in (3.1) is

$\mathbb{E}\\{Y(\mathbf{s};v)\\}=\mu(\mathbf{s};v)=\sum_{l=0}^{L}a_{l}(\mathbf{s})f_{l}(v),$

and the covariance $C_{Y}$ of $Y(\mathbf{s};v)$ is $C_{Y}=C_{\varepsilon}$,
which is zero at $(\mathbf{s}_{i},\mathbf{s}_{j})$ if $i\neq j$. Then, the
spatial variation is described through the mean of the functional random field
$Y(\mathbf{s};v)$. The covariates $a_{l}(\mathbf{s})$ capture the spatial
dependence, and the set of functions $\\{f_{l}(v)\\}_{l=0}^{L}$ carries the
continuity of the functional data.

The estimators of $f_{l}(v)$ in model (3.1) can be obtained using the OLS
method. The matrix form of the model can then be written as

$\mathbf{y}(v)=\mathbf{Af}(v)+\bm{\varepsilon}(v),$

where
$\mathbf{y}(v)=\\{y(\mathbf{s}_{1};v),\ldots,y(\mathbf{s}_{n};v)\\}^{\mathsf{T}}$,
$\mathbf{A}=\\{a_{l}(\mathbf{s}_{i})\\}_{i,l}$, $i=1,\ldots,n$,
$l=0,1,\ldots,L$, is the design matrix,
$\mathbf{f}(v)=\\{f_{0}(v),\ldots,f_{L}(v)\\}^{\mathsf{T}}$, and
$\bm{\varepsilon}(v)=\\{\varepsilon(\mathbf{s}_{1};v),\dots,\varepsilon(\mathbf{s}_{n};v)\\}^{\mathsf{T}}$.
Then, the OLS estimator is obtained by solving the optimization problem

$\min_{f_{0},\ldots,f_{L}\in\mathcal{H}}\sum_{i=1}^{n}\left\|y(\mathbf{s}_{i};\cdot)-\sum_{l=0}^{L}a_{l}(\mathbf{s}_{i})f_{l}\right\|^{2}_{\mathcal{H}}.$
(3.2)

Under some conditions (Menafoglio et al., 2013), (3.2) admits a unique
solution

$\hat{\mathbf{f}}(v)=(\mathbf{A}^{\mathsf{T}}\mathbf{A})^{-1}\mathbf{A}^{\mathsf{T}}\mathbf{y}(v).$

Thus, the drift estimator is obtained as

$\hat{\bm{\mu}}(v)=\mathbf{A}(\mathbf{A}^{\mathsf{T}}\mathbf{A})^{-1}\mathbf{A}^{\mathsf{T}}\mathbf{y}(v).$
(3.3)

Since the estimator $\hat{\bm{\mu}}(v)$ in (3.3) is a linear combination of
the observed curves, $\hat{\bm{\mu}}(v)$ inherits the continuity property of
$\mathbf{y}(v)$.

The basis functions approach offers another alternative to obtain estimators
of $\\{f_{l}(v)\\}$ (Reyes et al., 2015). Each component of the model (3.1)
can be assumed to be in the space generated by finite basis functions, i.e.,

$y(\mathbf{s}_{i};v)=\sum_{k=1}^{K}z_{ik}\eta_{k}(v),\,f_{l}(v)=\sum_{k=1}^{K}b_{lk}\eta_{k}(v),\mbox{
and }\varepsilon(\mathbf{s}_{i};v)=\sum_{k=1}^{K}e_{ik}\eta_{k}(v),$

for $i=1,\ldots,n$, and $l=0,1,\ldots,L$, where
$\bm{\eta}(v)=\\{\eta_{1}(v),\ldots,\eta_{K}(v)\\}^{\mathsf{T}}$ is the basis
function. In this case, the matrix form of the model (3.1) is

$\mathbf{Z}\bm{\eta}(v)=\mathbf{AB}\bm{\eta}(v)+\mathbf{E}\bm{\eta}(v),$

where $\mathbf{Z}=\\{z_{ik}\\}$, $\mathbf{B}=\\{b_{lk}\\}$, and
$\mathbf{E}=\\{e_{ik}\\}$, $i=1,\ldots,n$, $l=0,1,\ldots,L$, and
$k=1,\ldots,K$. The corresponding normal equation is

$\mathbf{A}^{\mathsf{T}}\mathbf{Z}\mathbf{J}_{\bm{\eta}}=\mathbf{A}^{\mathsf{T}}\mathbf{AB}\mathbf{J}_{\bm{\eta}},$

where
$\mathbf{J}_{\bm{\eta}}=\int\bm{\eta}(v)\bm{\eta}^{\mathsf{T}}(v)\mathrm{d}v$.
The solution for $\mathbf{B}$ is found by vectorizing the normal equation,
that is

$\mathrm{vec}(\hat{\mathbf{B}})=(\mathbf{J}_{\bm{\eta}}^{\mathsf{T}}\otimes\mathbf{A}^{\mathsf{T}}\mathbf{A})^{-1}\mathrm{vec}(\mathbf{A}^{\mathsf{T}}\mathbf{Z}\mathbf{J}_{\bm{\eta}}).$
(3.4)

Consequently, $\hat{\bm{\mu}}(v)=\mathbf{A}\hat{\mathbf{B}}\bm{\eta}(v)$.

In the context of explicative modeling, we can use other spatial functional
data as covariates to describe the mean $\mu(\mathbf{s};v)$ (Ignaccolo et al.,
2014). For example, the mean can be modeled as
$\mu(\mathbf{s};v)=\beta_{0}(v)+\sum_{p=1}^{P}\beta_{p}(v)U_{p}(\mathbf{s};v)$,
where $\\{U_{p}(\mathbf{s};v)\\}_{p=1}^{P}$ are the functional covariates, and
$\beta_{0}(v),\beta_{1}(v),\ldots,\beta_{P}(v)$ are the functional parameters
to be estimated.

###  3.2 Small-scale Variation

Small-scale variations are usually represented through the covariance
structure. The modeling of covariance is one of the most studied subjects in
spatial statistics (Stein, 1999; Cressie, 2015; Genton and Kleiber, 2015). For
finite-dimensional data, the exponential, the Gaussian, and the Matérn are
examples of parametric covariance function models.

In this section, we assume that the mean is constant over locations, and
without loss of generality set to be zero. A model to describe the small-scale
variations is

$Y(\mathbf{s};v)=X(\mathbf{s};v)+\varepsilon(\mathbf{s};v),$ (3.5)

where $\varepsilon(\mathbf{s};v)$ represents the functional white noise, and
is assumed to be uncorrelated with $X(\mathbf{s};v)$.

The mean of $Y(\mathbf{s};v)$ in model (3.5) is zero, and the covariance is
such that

$C_{Y}(\mathbf{s}_{i},\mathbf{s}_{j};f)=C_{X}(\mathbf{s}_{i},\mathbf{s}_{j};f)+\mathds{1}(i=j)C_{\varepsilon}(\mathbf{s}_{i},\mathbf{s}_{j};f),\quad
f\in\mathcal{H},$ (3.6)

where $\mathds{1}(\cdot)$ is the indicator function. Thus, we need to consider
the additional term
$C_{\varepsilon}(\mathbf{s},\mathbf{s};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$
in the estimation. Similarly, as in Section 2.2, a basis function approach can
be used to obtain the estimators.

Let
$\mathbf{Z}(\mathbf{s})=\\{Z_{1}(\mathbf{s}),\ldots,Z_{K}(\mathbf{s})\\}^{\mathsf{T}}$
be the multivariate random field obtained from the projection of
$Y(\mathbf{s};v)$ onto the basis functions $\\{\eta_{k}\\}_{k=1}^{K}$, that
is, each component of $\mathbf{Z}(\mathbf{s})$ is defined as
$Z_{k}(\mathbf{s})=\langle Y(\mathbf{s};\cdot),\eta_{k}\rangle$, which are
scalar random fields with mean zero. From (3.5), for each $k=1,\ldots,K$, we
have that the process $Z_{k}(\mathbf{s})$ is such that

$Z_{k}(\mathbf{s})=\langle
X(\mathbf{s}),\eta_{k}\rangle+\langle\varepsilon(\mathbf{s}),\eta_{k}\rangle.$

The variance of $Z_{k}(\mathbf{s})$ is

$\displaystyle\mathbb{E}\\{\langle Y(\mathbf{s};\cdot),\eta_{k}\rangle\langle
Y(\mathbf{s};\cdot),\eta_{k}\rangle\\}$ $\displaystyle=\langle
C_{Y}(\mathbf{s},\mathbf{s};\eta_{k}),\eta_{k}\rangle$ $\displaystyle=\langle
C_{X}(\mathbf{s},\mathbf{s};\eta_{k}),\eta_{k}\rangle+\langle
C_{\varepsilon}(\mathbf{s},\mathbf{s};\eta_{k}),\eta_{k}\rangle,$

where the last equality is obtained using (3.6). This implies that the random
field $Z_{k}(\mathbf{s})$ has a nugget effect $\langle
C_{\varepsilon}(\mathbf{s},\mathbf{s};\eta_{k}),\eta_{k}\rangle$ that should
be considered when fitting a covariance (variogram) model.

The covariance estimator of $Y(\mathbf{s};v)$ is obtained from (2.10), after
estimating the covariance
$\bm{\Sigma}(h)=\mathrm{Cov}\\{\mathbf{Z}(\mathbf{s}_{1}),\mathbf{Z}(\mathbf{s}_{2})\\}$,
$h=\|\mathbf{s}_{1}-\mathbf{s}_{2}\|$, from the data $\\{\langle
y(\mathbf{s}_{i};\cdot),\eta_{k}\rangle\\}_{i=1}^{n}$, $k=1,\ldots,K$. These
ideas are extended to estimating the variogram of $Y(\mathbf{s};v)$ by
estimating the matrix-variogram of $\mathbf{Z}(\mathbf{s})$.

In the case of the trace-variogram of $Y(\mathbf{s};v)$, from (3.5) we have
that the mean of the inner product satisfies $\mathbb{E}\\{\langle
Y(\mathbf{s}_{i};\cdot),Y(\mathbf{s}_{j};\cdot)\rangle\\}=\mathbb{E}\\{\langle
X(\mathbf{s}_{i};\cdot),X(\mathbf{s}_{j};\cdot)\rangle\\}+\mathbb{E}\\{\langle\varepsilon(\mathbf{s}_{i};\cdot),\varepsilon(\mathbf{s}_{j};\cdot)\rangle\\}$.
That is, the trace-covariogram of $Y(\mathbf{s};v)$ is such that

$\sigma_{\mathrm{tr,Y}}(\mathbf{s}_{1},\mathbf{s}_{2})=\sigma_{\mathrm{tr,X}}(\mathbf{s}_{1},\mathbf{s}_{2})+\mathds{1}(\mathbf{s}_{1}=\mathbf{s}_{2})\sigma_{\mathrm{tr,\varepsilon}}(\mathbf{s}_{1},\mathbf{s}_{2}).$

Then, when estimating $\sigma_{\mathrm{tr,Y}}(\mathbf{s}_{1},\mathbf{s}_{2})$
using $\\{\langle
y(\mathbf{s}_{i};\cdot),y(\mathbf{s}_{j};\cdot)\rangle\\}_{i,j}$, as described
in Section 2.2.2, one should consider the nugget effect
$\sigma_{\mathrm{tr,\varepsilon}}(\mathbf{s},\mathbf{s})$. The same is true in
the case of the trace-variogram.

###  3.3 Large-scale and Small-scale Variations

Datasets often have both a trend component (large-scale variation) and a
spatial variability (small-scale variation), e.g., temperature data show an
increasing tendency and a small variability around this tendency. In the
context of functional data correlated only in time, Martínez-Hernández and
Genton (2020) proposed a method to estimate trend using tensor product
surfaces.

A model for spatial functional data can be written as

$Y(\mathbf{s};v)=\mu(\mathbf{s};v)+X(\mathbf{s};v)+\varepsilon(\mathbf{s};v).$

As before, the mean can be expressed as
$\mu(\mathbf{s};v)=\sum_{l=0}^{L}a_{l}(\mathbf{s})f_{l}(v)$. Similarly to
(3.1), the model to estimate the parameters of the mean can be written as

$Y(\mathbf{s};v)=\sum_{l=0}^{L}a_{l}(\mathbf{s})f_{l}(v)+\epsilon(\mathbf{s};v),$
(3.7)

where the residual
$\epsilon(\mathbf{s};v):=X(\mathbf{s};v)+\varepsilon(\mathbf{s};v)$ is now a
functional random field with mean zero and covariance
$C_{\epsilon}(\mathbf{s}_{1},\mathbf{s}_{2};\mathbin{\mathchoice{\vbox{\hbox{\set@color$\displaystyle\centerdot$}}}{\vbox{\hbox{\set@color$\textstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptstyle\centerdot$}}}{\vbox{\hbox{\set@color$\scriptscriptstyle\centerdot$}}}})$.
Thus, unlike model (3.1), in which the residuals are not correlated, model
(3.7) has spatially correlated residuals. Because of this correlation, we can
use the GLS method instead of the OLS method (Menafoglio et al., 2013). Let
$\bm{\Sigma}$ be the trace-variogram matrix of $\epsilon(\mathbf{s}_{i};v)$ at
different distances of the locations $\mathbf{s}_{i}$. Then, the GLS estimator
of $\mu(\mathbf{s}_{i};v)$ is

$\hat{\bm{\mu}}(v)=\mathbf{A}(\mathbf{A}^{\mathsf{T}}\bm{\Sigma}^{-1}\mathbf{A})^{-1}\mathbf{A}^{\mathsf{T}}\bm{\Sigma}^{-1}\mathbf{y}(v),$
(3.8)

where $\mathbf{A}$ is the design matrix, and $\mathbf{y}(v)$ the evaluation of
$y(\mathbf{s};v)$ at $n$ locations, both defined in Section 3.1.

The trace-variogram matrix $\bm{\Sigma}$ cannot be estimated directly from
$\epsilon(\mathbf{s}_{i};v)$, because we do not observe them. The approach
commonly used, to solve this problem, is to use the following iterative
procedure:

  1. 1.

Compute the initial estimator $\hat{\mu}_{0}(\mathbf{s}_{i};v)$ as in (3.3),
assuming that $\\{\epsilon(\mathbf{s}_{i};v)\\}_{i=1}^{n}$ are not spatially
correlated.

  2. 2.

Compute the residuals
$\hat{\epsilon}(\mathbf{s}_{i};v)=y(\mathbf{s}_{i};v)-\hat{\mu}_{0}(\mathbf{s}_{i};v)$,
for $i=1,\ldots,n$.

  3. 3.

Estimate the initial empirical trace-variogram
$\hat{\gamma}_{\mathrm{tr},0}(h)$ as in (2.12), using
$\\{\hat{\epsilon}(\mathbf{s}_{i};v)\\}_{i=1}^{n}$, and then, obtain an
estimator $\hat{\bm{\Sigma}}_{0}$ of $\bm{\Sigma}$ by fitting a parametric
model with a nugget effect.

  4. 4.

Re-estimate the mean to get $\hat{\mu}_{1}(\mathbf{s}_{i};v)$ by using (3.8)
with $\hat{\bm{\Sigma}}_{0}$.

  5. 5.

Repeat steps $2-4$ until convergence.

Once the mean $\mu(\mathbf{s};v)$ is estimated, it is removed from the data.
Then, the covariance is estimated as in Section 3.2, using the spatial
functional data
$\\{y(\mathbf{s}_{i};v)-\hat{\mu}(\mathbf{s}_{i};v)\\}_{i=1}^{n}$.

The iterative procedure can also be performed with the finite basis function
$\bm{\eta}(v)$. In step $1$, we can use (3.4), and in step $4$, use the
estimator
$\mathrm{vec}(\hat{\mathbf{B}})=(\mathbf{J}_{\bm{\eta}}^{\mathsf{T}}\otimes\mathbf{A}^{\mathsf{T}}\hat{\bm{\Sigma}}_{0}\mathbf{A})^{-1}\mathrm{vec}(\mathbf{A}^{\mathsf{T}}\hat{\bm{\Sigma}}_{0}\mathbf{Z}\mathbf{J}_{\bm{\eta}})$
to obtain $\hat{\mu}_{1}$.

##  4 Kriging for Functional Random Fields

In spatial statistics, the concept of kriging (co-kriging for the multivariate
setting) is a synonym of optimal interpolation. The main goal is to be able to
predict at locations where data are not observed. This predictor is a linear
combination of the observed data, such that it is the best linear unbiased
predictor under squared loss. Here, we briefly mention the concept of kriging
and we redirect readers to the references provided.

Let $\mathbf{s}_{0}\in D$ be the location at which the curve will be
predicted. Let $\\{\Psi_{1},\ldots,\Psi_{n}\\}$ be linear operators from
$\mathcal{H}$ to $\mathcal{H}$. In general, kriging can be defined as

$\hat{x}(\mathbf{s}_{0};v)=\sum_{i=1}^{n}\Psi_{i}\\{x(\mathbf{s}_{i};v)\\},$
(4.1)

where the coefficients $\Psi_{i}$ are obtained by minimizing the square norm
of the error prediction, $\hat{x}(\mathbf{s}_{0};v)-X(\mathbf{s}_{0};v)$, with
an additional constraint of unbiasedness. That is,

$\displaystyle\min_{\Psi_{1},\ldots,\Psi_{n}}$
$\displaystyle\mathbb{E}\\{\langle\hat{x}(\mathbf{s}_{0};\cdot)-X(\mathbf{s}_{0};\cdot),\hat{x}(\mathbf{s}_{0};\cdot)-X(\mathbf{s}_{0};\cdot)\rangle\\}$
s.t.
$\displaystyle\mathbb{E}\\{\hat{x}(\mathbf{s}_{0};v)\\}=\mathbb{E}\\{X(\mathbf{s}_{0};v)\\}.$

A particular case of the coefficients operators $\Psi_{i}$ in (4.1) is the so-
called kernel operators, which are defined as
$\Psi_{i}(f)(v)=\int_{0}^{1}\\!\lambda_{i}(v,u)f(u)\mathrm{d}u$,
$f\in\mathcal{H}$. In this case, the estimation is through functions
$\lambda_{i}(v,u)$. Other cases are $\Psi_{i}(f)(v)=\lambda_{i}(v)f(v)$ and
$\Psi_{i}(f)(v)=w_{i}f(v)$, with $w_{i}$ scalars. The latter corresponds to a
simple ponderation of the observed curves, i.e.,
$\hat{x}(\mathbf{s}_{0};v)=\sum_{i=1}^{n}w_{i}x(\mathbf{s}_{i};v)$. All these
cases can be fitted into ordinary or universal kriging. For ordinary kriging,
see Goulard and Voltz (1993), Nerini et al. (2010), Giraldo et al. (2010), and
Giraldo et al. (2011), as well as the review paper by Delicado et al. (2010).
For universal kriging, see Menafoglio et al. (2013), Caballero et al. (2013),
Reyes et al. (2015), and Menafoglio et al. (2016). For co-kriging
(multivariate functional random fields), see Bohorquez et al. (2017) and
Grujic et al. (2018). An alternative approach to kriging, which is based on a
tensor function space, can be found in Aguilera-Morillo et al. (2017).

##  5 Surface Time Series

In many phenomena, data can be collected in the form of a surface, called
surface data in our case, or manifold data for more complex structures. For
example, we can have data arising from neuroimaging (Lila et al., 2016), from
two-dimensional time-frequency domains (Aston et al., 2017), from satellite
images (Zhang et al., 2011), and functional data with two-dimensional domain
(Crainiceanu et al., 2011; Morris et al., 2011). Here, we consider surface
data as functional, that is, the atoms of the functional random variable are
continuous surfaces. Surface data provide an alternative approach to analyzing
spatial data, where the continuous realization of a random field is considered
as a unit. This approach can have computational advantages, especially if the
locations, where data are observed, are dense in space. Particularly, spatio-
temporal data (Cressie and Wikle, 2011) can be considered as surface data that
are observed over time (surface time series). Surface data capture the spatial
dependence through the continuity of the surface (see, e.g., Bernardi et al.,
2017). Moreover, the approach of surface data can be applied to a nonplanar
spatial domain, such as a sphere or a general two-dimentional manifold (Dassi
et al., 2015; Wilhelm et al., 2016; Ettinger et al., 2016; Menafoglio and
Secchi, 2017; Greco et al., 2018). Kriging can be applied to these complex
domains by using an appropriate distance, but the covariance models do not
necessarily guarantee a positive definite covariance, e.g., the Matérn
covariance (Gneiting, 2013). Here, we focus on spatio-temporal data, where at
each time point we observe a surface.

The functional time series approach to spatial statistics has been studied by
Ruiz-Medina et al. (2007). In Aston et al. (2017), a tensor product Hilbert
space was considered to propose a separability test for the covariance
operators of random surfaces.

###  5.1 Basic Concepts

Let $(\Omega,F,P)$ be a probability space, and let $\mathcal{H}$ be the
Hilbert space defined as the set of functions with domain
$D\subset\mathbb{R}^{2}$,
$\mathcal{H}=\\{f:D\to\mathbb{R}:\int_{D}|f(\mathbf{s})|^{2}\mathrm{d}\mathbf{s}<\infty\\}$,
and with the inner product $\langle
f,g\rangle=\int_{D}f(\mathbf{s})g(\mathbf{s})\mathrm{d}\mathbf{s}$. The norm
induced by the inner product is denoted by $\|\cdot\|_{\mathcal{H}}$. Thus, a
random variable $X:\Omega\to\mathcal{H}$ is a functional random variable with
a surface as atom. We denote by $X(\mathbf{s})$ this functional random
variable with $\mathbf{s}=(s_{1},s_{2})\in D$.

A functional time series is a sequence of functional random variables
$\\{X_{t}(\mathbf{s});t\in\mathbb{Z}\\}$ in $\mathcal{H}$. Bosq (2000) is a
monograph on linear processes in function spaces, including functional time
series in a Hilbert space. Hörmann and Kokoszka (2012) reviewed functional
time series. Also, see Ramsay and Silverman (2005), Ferraty and Vieu (2006),
and Horváth and Kokoszka (2012) for a general introduction to functional data
analysis. In this paper, we assume that
$\mathbb{E}\\{\|X_{t}(\mathbf{s})\|_{\mathcal{H}}^{2}\\}<\infty$ for all
$t\in\mathbb{Z}$.

The mean of the surface time series $\\{X_{t}(\mathbf{s})\\}$ is defined as
$\mu_{t}(\mathbf{s})=\mathbb{E}\\{X_{t}(\mathbf{s})\\}$, where
$\mu_{t}(\mathbf{s})$ is such that $\mathbb{E}(\langle
X_{t},f\rangle)=\langle\mu_{t},f\rangle$ for all $f\in\mathcal{H}$. The
covariance function at lag $h\in\mathbb{N}$ is defined as

$C_{X_{t-h},X_{t}}(f)=\mathbb{E}\\{\langle
X_{t-h}-\mu_{t-h},f\rangle(X_{t}-\mu_{t})\\},\quad f\in\mathcal{H}.$

This covariance function can be expressed as

$C_{X_{t-h},X_{t}}(f)(\cdot)=\int\sigma_{t-h,t}(\cdot,\mathbf{s})f(\mathbf{s})\mathrm{d}\mathbf{s},$

where
$\sigma_{t-h,t}(\mathbf{s}_{1},\mathbf{s}_{2})=\mathrm{Cov}\\{X_{t-h}(\mathbf{s}_{1}),X_{t}(\mathbf{s}_{2})\\}$.
The stationarity condition is important for statistical inference.

######  Definition 3 (Weak stationarity)

A surface time series $\\{X_{t}(\mathbf{s});t\in\mathbb{Z}\\}$ is said to be
(weakly) stationary if

  1. 1.

$\mathbb{E}\\{X_{t}(\mathbf{s})\\}=\mu(\mathbf{s})$ for all $t\in\mathbb{Z}$,
and

  2. 2.

$C_{X_{t_{1}+h},X_{t_{2}+h}}(f)=C_{X_{t_{1}},X_{t_{2}}}(f)$ for all
$t_{1},t_{2}\in\mathbb{Z}$, $h\in\mathbb{N}$, and $f\in\mathcal{H}$.

If the surface time series is stationary, we write $C_{h}$ for the covariance
functions instead of $C_{X_{t},X_{t+h}}$. The definition of stationarity does
not require stationarity over the space $D$, e.g., for each $t$, data can be a
realization of a nonstationary random field. In general, the covariance
function $C_{0}$ describes only the spatial dependence, whereas $C_{h}$, for
$h\neq 0$, describes the dependency over time of the surface time series.

A surface white noise is a stationary surface time series with a zero mean and
a covariance function $C_{h}=0$, if $h\neq 0$. Thus, the surface white noise
can have a spatial correlation at each time point, but not across time.

Similarly as before, the eigenfunctions are defined as functions
$\zeta\in\mathcal{H}$ such that

$C_{0}(\zeta)(\mathbf{s})=\lambda\zeta(\mathbf{s}),$

where $\lambda$ is positive and is the corresponding eigenvalue. Moreover, the
covariance operator $C_{0}$ can be decomposed in terms of the eigenfunctions,
that is

$C_{0}(f)(\mathbf{s})=\sum_{j=1}^{\infty}\lambda_{j}\langle\zeta_{j},f\rangle\zeta_{j}(\mathbf{s}),$

where $\zeta_{j}$, $j=1,2,\ldots$, are the eigenfunctions of $C_{0}$ with
eigenvalues $\lambda_{j}$. The eigenvalues are such that
$\sum_{j=1}^{\infty}\lambda_{j}=\mathbb{E}\\{\|X_{0}(\mathbf{s})\|^{2}_{\mathcal{H}}\\}<\infty$.
This operator $C_{0}$ is nuclear and therefore Hilbert-Schmidt. Similarly as
in the finite dimensional case, eigenfunctions are important to reduce
dimensionality.

###  5.2 Estimation

Now, we describe estimators of the mean $\mu(\mathbf{s})$ and the covariance
operator $C_{h}$. Let $\\{x_{t}(\mathbf{s})\\}_{t=1}^{T}$ be a realization of
a stationary surface time series $X_{t}(\mathbf{s})$ with mean
$\mu(\mathbf{s})$. The sample mean is defined as

$\hat{\mu}(\mathbf{s})=\frac{1}{T}\sum_{t=1}^{T}x_{t}(\mathbf{s}).$ (5.1)

The sample mean $\hat{\mu}(\mathbf{s})$ is an unbiased estimator of
$\mu(\mathbf{s})$.

The empirical covariance at lag $h$ of $X_{t}(\mathbf{s})$ is defined as

$\hat{C}_{h}(f)(\mathbf{s})=\frac{1}{T-h}\sum_{t=1}^{T-h}\langle
x_{t}-\hat{\mu},f\rangle\\{x_{t+h}(\mathbf{s})-\hat{\mu}(\mathbf{s})\\},\quad
f\in\mathcal{H},$ (5.2)

and the corresponding empirical kernel is defined as

$\hat{\sigma}_{h}(\mathbf{s}_{1},\mathbf{s}_{2})=\frac{1}{T-h}\sum_{t=1}^{T-h}\\{x_{t+h}(\mathbf{s}_{1})-\hat{\mu}(\mathbf{s}_{1})\\}\\{x_{t}(\mathbf{s}_{2})-\hat{\mu}(\mathbf{s}_{2})\\}.$

The empirical covariance operator $\hat{C}_{h}$ is an unbiased estimator of
$C_{h}$, see Bosq (2000). For papers related to the mean and the covariance
functions, refer to previous studies by authors Hörmann and Kokoszka (2010),
Horváth et al. (2013), and Horváth et al. (2014).

In time series analysis, there is no direct modeling of the covariance
function. Instead, a model for the process is proposed, and the covariance
function is derived from the model. We adopt this idea in the next section
(Section 5.3).

###  5.3 Modeling

This section discusses two topics: the continuous estimation of the surface
and the modeling of the continuous surface series.

####  5.3.1 Estimating the continuous surface

In practice, data are observed on a finite set of points, i.e., for each
$t=1,\ldots,T$, we observe $n_{t}$ points
$\\{y_{t,i}=y_{t}(\mathbf{s}_{i})\\}_{i=1}^{n_{t}}$ of the functional data
$Y_{t}(\mathbf{s})$ on a set of points
$\\{\mathbf{s}_{1},\ldots,\mathbf{s}_{n_{t}}\\}\subset D$, and possibly with
measurement errors. Thus, for each $t$, the unknown surface (deterministic
field) $y_{t}(\mathbf{s})$ needs to be estimated. This continuous surface
estimate can be associated with kriging in classical spatial data analysis.

Because the procedure to estimate $y_{t}(\mathbf{s})$ is independent of $t$,
we drop the subindex $t$ in the sequel, and we consider $n$ as the sample
size. Thus, a model of $y(\mathbf{s})$ can be written as

$y_{i}=y(\mathbf{s}_{i})+\varepsilon_{i},$

where $\\{\varepsilon_{i}\\}_{i=1}^{n}$ represent the measurement errors that
are spatially uncorrelated. The function $y(\mathbf{s})$ describes the spatial
structure of the phenomenon being studied.

To estimate $y(\mathbf{s})$, one can extend the smoothness techniques of the
curves described in Ramsay and Silverman (2005) to surfaces (manifolds). In
particular, one can extend the spline smoothing. The extension of spline
smoothing to surfaces is an important research area. Extensions have been done
over Euclidean domains and non-Euclidean domains, including the spherical
domain. One extension is to use the tensor product of univariate B-splines
(Eilers and Marx, 1996; Wood, 2006; Qingguo and Longsheng, 2010; Xiao et al.,
2013). In this case, an estimator of $y(\mathbf{s})$ has the form

$\hat{y}(\mathbf{s})=\sum_{k=1}^{K_{1}}\sum_{l=1}^{K_{2}}\theta_{kl}\eta_{k}(s_{1})\nu_{l}(s_{2}),$

where $\\{\eta_{k}\\}_{k=1}^{K_{1}}$ and $\\{\nu_{l}\\}_{l=1}^{K_{2}}$ are
B-splines basis functions for $s_{1}$ and $s_{2}$ coordinates, respectively,
with $\mathbf{s}=(s_{1},s_{2})$. In the estimation procedure, smoothness
properties are imposed through a penalization term. To read about some
approaches of bivariate smoothing, see Ruppert et al. (2003), Lai and
Schumaker (2007), and Wood (2017).

In general, the estimation of $y(\mathbf{s})$ can be formulated as the
minimization of the sum of squared errors with a penalization term. The
penalization term measures the roughness of the fitted surface and can carry
partial information of $y(\mathbf{s})$. That is, the optimization problem is
written as

$\sum_{i=1}^{n}\\{y_{i}-y(\mathbf{s}_{i})\\}^{2}+\lambda P(y),$ (5.3)

where $P(y)$ is the penalization term, and $\lambda$ the smoothness parameter
which controls the smoothness of the estimated surface. A popular penalization
is the thin-plate energy, which is defined as
$P(y)=\int\\{(\frac{\partial^{2}y}{\partial
s^{2}_{1}})^{2}+2(\frac{\partial^{2}y}{\partial s_{1}\partial
s_{2}})^{2}+(\frac{\partial^{2}y}{\partial
s^{2}_{2}})^{2}\\}\mathrm{d}\mathbf{s}$. The resulting estimator is called the
thin plate splines (Duchon, 1977). A Bayesian adaptive thin plate spline was
proposed in Yue and Speckman (2010). Another example of penalization involves
the Laplacian, that is $P(y)=\int(\frac{\partial^{2}y}{\partial
s^{2}_{1}}+\frac{\partial^{2}y}{\partial s^{2}_{2}})^{2}\mathrm{d}\mathbf{s}$
(Wood et al., 2008; Sangalli et al., 2013). The definition of the penalty term
depends on each specific problem. For example, in the Laplacian case, the
unique penalty parameter $\lambda$ that controls both directions $s_{1}$ and
$s_{2}$ implies an isotropic smoothing. In contrast, if
$P(y)=\lambda_{1}P_{1}(y)+\lambda_{2}P_{2}(y)$, where $P_{i}(y)$ is a penalty
term in the $i$-th coordinate, then it results in anisotropic smoothing.

In general, the penalty term can be defined in terms of a partial differential
equation (PDE). For example, $P(y)=\int(Ly-u)^{2}\mathrm{d}\mathbf{s}$, where
$L$ is a differential operator and $Ly=u$ is a PDE (see, e.g., Azzimonti et
al., 2015; Sangalli, 2019). The advantage of the penalty term with PDE is that
it can handle complex domains with boundary conditions or interior holes, and
varies depending on the phenomena being studied. The PDE is such that it
contains information about the phenomena, and it regularizes the estimation
with values of $\lambda$. The solution of (5.3) may not have a closed form,
but it can be approximated by using finite elements analysis, see Ramsay
(2002), Duchamp and Stuetzle (2003), and Sangalli et al. (2013).

We describe the solution of (5.3) using the finite elements analysis
technique. Let $\mathcal{M}$ be a mesh of $D$. Let
$\\{\phi_{1}(\mathbf{s}),\ldots,\phi_{K}(\mathbf{s})\\}$ be basis functions
that are piece-wise polynomials associated with the mesh $\mathcal{M}$. Then,
the estimator of $y(\mathbf{s})$ is assumed to have the form

$\hat{y}(\mathbf{s})=\sum_{k=1}^{K}\beta_{k}\phi_{k}(\mathbf{s}),$

where the coefficients $\bm{\beta}=(\beta_{1},\ldots,\beta_{K})^{\mathsf{T}}$
need to be estimated. Let $\mathbf{y}=(y_{1},\ldots,y_{n})^{\mathsf{T}}$ be
the observed values over $D$, and let $\mathbf{P}$ be the discretization of
the penalty. Then, the estimator of $\bm{\beta}$ has the form

$\hat{\bm{\beta}}=\left(\varPhi^{\mathsf{T}}\varPhi+\lambda\mathbf{P}\right)^{-1}\varPhi^{\mathsf{T}}\mathbf{y},$

where $\varPhi=\\{\phi_{k}(\mathbf{s}_{i})\\}_{i,k=1}^{n,K}$ is the $n\times
K$ matrix which represents the evaluation of each basis function at the
locations at which data are observed.

The time component can also be considered in the PDE. Arnone et al. (2019)
proposed general forms of time-dependent PDEs in the context of time dependent
surface data.

Additional studies of data over complex domains have been published. Wang and
Ranalli (2007) proposed a modified thin plate spline over a complex domain;
Lindgren et al. (2011) linked Gaussian fields via stochastic partial
differential equations, where the solution is found using a finite elements
analysis; Scott-Hayward et al. (2014) proposed a complex region spatial
smoother using the geodesic distance; and Menafoglio et al. (2018) proposed a
methodology for spatial fields of object data over complex domains.

####  5.3.2 Functional autoregressive models

Here, we assume that the functional time series consist of continuous surfaces
(deterministic fields) that can be estimated with the methods described in
Section 5.3.1. In the context of functional time series, the most popular
model is the functional autoregressive model of order $P$, FAR$(P)$. A surface
time series $\\{X_{t}(\mathbf{s})\\}$ follows the FAR$(P)$ model if
$X_{t}(\mathbf{s})=\sum_{p=1}^{P}\Psi_{p}(X_{t-p})(\mathbf{s})+W_{t}(\mathbf{s})$,
where each coefficient $\Psi_{p}:\mathcal{H}\to\mathcal{H}$ is an operator,
and $\\{W_{t}(\mathbf{s})\\}$ is a surface white noise. In practice, the order
$P$ needs to be estimated (Kokoszka and Reimherr, 2013). Here, we assume $P=1$
to illustrate the ideas and to simplify the notations.

Let $\\{y_{t}(\mathbf{s})\\}_{t=1}^{T}$ be the surface data observed over time
$t=1,\ldots,T$, and assume that it is a stationary surface time series. Then,
the dependency over time can be modeled by using the FAR$(1)$ process, that is

$\displaystyle y_{t}(\mathbf{s})$
$\displaystyle=\mu(\mathbf{s})+x_{t}(\mathbf{s}),$ (5.4) $\displaystyle
x_{t}(\mathbf{s})$
$\displaystyle=\Psi(x_{t-1})(\mathbf{s})+W_{t}(\mathbf{s}),$ (5.5)

where $\mu(\mathbf{s})$ represents the surface mean, i.e., the large-scale
variation for all $t$. The unobserved $x_{t}(\mathbf{s})$ follows a stationary
FAR$(1)$ process with mean zero, and $\\{W_{t}(\mathbf{s})\\}$ is a surface
white noise. The surface white noise $\\{W_{t}(\mathbf{s})\\}$ can be
interpreted as the surface data components that describe the small-scale
variation for each time $t$, which are not correlated over time. The
dependency over time is driven by the operator $\Psi$.

Since $\\{y_{t}(\mathbf{s})\\}$ is assumed to be stationary, the estimation of
$\mu(\mathbf{s})$ can be obtained as in (5.1). After removing the mean, the
rest of the analysis is performed on the process
$x_{t}(\mathbf{s})=y_{t}(\mathbf{s})-\hat{\mu}_{t}(\mathbf{s})$.

Now, we focus on the estimation of the coefficient operator $\Psi$. Let
$C_{h}$ be the covariance operator of the FAR$(1)$ process $X_{t}(\mathbf{s})$
with realization $\\{x_{t}(\mathbf{s})\\}$. Then, the covariance operator of
$X_{t}(\mathbf{s})$ satisfies

$C_{1}(f)=\Psi\\{C_{0}(f)\\},\quad f\in\mathcal{H}.$ (5.6)

Moreover, for $h\in\mathbb{N}$, we have that
$C_{h}(f)=\Psi^{h}\\{C_{0}(f)\\}$. So if the coefficient $\Psi$ and $C_{0}$
are known, then we can compute the covariance function at any lag
$h\in\mathbb{N}$. Thus, the estimation of $\Psi$ is crucial.

To obtain an estimator of $\Psi$, we can use the estimation of $C_{0}$ and
$C_{1}$ in (5.6), defining $\Psi(f)=C_{1}\\{C_{0}^{-1}(f)\\}$, if $C_{0}$ is
invertible. In principle, we can always use the estimator (5.2) for $h=0,1$,
and compute the inverse $\hat{C}^{-1}_{0}$. However, when the sample size $T$
tends to infinity, $\hat{C}^{-1}_{0}$ becomes unbounded (Cardot et al., 1999).
This is because $C_{0}$ is a compact operator (Bosq, 2000). Thus, it is
necessary to use some regularization methods to obtain $C_{0}^{-1}$. That is,
$(C_{0}+\alpha_{T})^{-1}$ is computed instead of $C_{0}^{-1}$ where
$\alpha_{T}>0$ and $\alpha_{T}\downarrow 0$. Alternatively, $C_{0}^{-1}$ can
be approximated by using only the first $k$ eigenfunctions corresponding to
the largest eigenvalues, that is,
$C_{0}^{-1}=\sum_{j=1}^{k}\lambda_{j}^{-1}\zeta_{j}\otimes\zeta_{j}$, see Bosq
(2000) and Kokoszka and Reimherr (2017). Let $\hat{C}_{0}^{-1}$ denote an
estimator of the inverse operator, either using some regularization method,
finite eigenfunctions or other methods (Martínez-Hernández et al., 2019).
Then, the estimator of $\Psi$ is defined as

$\hat{\Psi}(f)=\hat{C}_{1}\\{\hat{C}_{0}^{-1}(f)\\},\quad f\in\mathcal{H}.$

Once $\mu(\mathbf{s})$ and $\Psi$ are estimated, then the one-step-ahead
prediction is obtained as

$\hat{y}_{T+1}(\mathbf{s})=\hat{\mu}(\mathbf{s})+\hat{x}_{T+1}(\mathbf{s}),$

where $\hat{x}_{T+1}(\mathbf{s})=\hat{\Psi}(x_{T})(\mathbf{s}).$ The
estimators of $\mu$, $\Psi$, and $C_{h}$ can be explicitly expressed in terms
of the basis functions
$\\{\phi_{1}(\mathbf{s}),\ldots,\phi_{K}(\mathbf{s})\\}$ from the finite
elements technique. Alternative approaches can be used to predict data, for
example one can extend the ideas described in Hyndman and Ullah (2007) and Aue
et al. (2015) to surfaces.

The approach of the surface time series provides alternatives to analyze
spatio-temporal data. The advantage of this approach is that it can handle
data collected on a large scale for each time point, over a general domain,
Euclidean, or non-Euclidean. With this approach, the process is modeled
instead of the covariance, which is convenient when the classical covariance
models are not guaranteed to be positive definite.

##  6 Statistical Software

Here, we mention some packages available in R for SFD and surface data. The
package fda (Ramsay et al., 2018) provides several commands to analyze and
construct continuous functions. It contains several options of basis
functions, such as Fourier basis functions and spline basis functions. The
basis functions can be used to estimate the continuous functions as in (2.8).
Once the coefficients are obtained in (2.8), we can use packages for classical
spatial data, e.g., spatial (Venables and Ripley, 2002), gstat (Pebesma,
2004), RandomFields (Schlather et al., 2019), fields (Nychka et al., 2017),
geoR (Ribeiro Jr and Diggle, 2018), ExaGeoStatR (Abdulah et al., 2019) for
large datasets, and spBayes Finley et al. (2015) for Bayesian analysis of
hierarchical multivariate models. Thus, when data are expressed in terms of
basis functions as in (2.8), we can combine the fda package and the packages
for spatial data to obtain estimators of the mean and the covariance functions
described in Section 2.2.

The package geofd by Giraldo et al. (2015) implements kriging of functional
data described in Section 4. The curves observed are pre-processed by fitting
Fourier or B-splines basis functions. Also, this package provides a command to
compute the trace-variogram defined in (2.12). Another package related to
Section 4 is fdagstat (Grujic and Menafoglio, 2017). This package implements
kriging, cokriging, and universal kriging, and includes the large-scale
variation described in Sections 3.1 and 3.3.

For the surface data described in Section 5, one can use the package mgcv by
Wood (2017), which allows us to smooth surfaces. The package fdaPDE (Lila et
al., 2019) implements smoothing with PDE penalization described in Sangalli et
al. (2013) and Azzimonti et al. (2015). INLA Rue et al. (2009) can be used to
estimate continuous surfaces. The package Manifoldgstat by Sartori and
Torriani (2019) implements kriging for manifold-valued random fields.

Some visualization tools for functional data and functional time series are
the functional boxplots (Sun and Genton, 2011, 2012) implemented in the fda
package; for functional images and surfaces, surface boxplots are used (Genton
et al., 2014). These tools are based on an order induced by a depth notion for
functional data.

##  7 Data Analysis

In this section, we illustrate general ideas of modeling spatial functional
data and surface data, without a deep statistical analysis of the data. Our
goal is to provide a general example, using only available packages.

We use wind data simulated from the Weather Research and Forecasting (WRF)
model by Yip (2018). Each measurement corresponds to hourly wind speed from
$2009$ to $2014$, in a $115\times 115$ km region centered in Dumat Al Jandal,
Saudi Arabia. That is where the first wind farm of the country is being built.
These wind speeds are simulated on a regular grid of points in space, namely
at $5$-km resolution.

###  7.1 Spatial Functional Data Approach

Our data are the daily wind speed, where $y(\mathbf{s}_{i};v_{j})$ is the wind
speed at location $\mathbf{s}_{i}$ for hour $v_{j}$, for
$v_{j}=1,2,\ldots,24$. With the functional approach, we consider
$y(\mathbf{s}_{i};\cdot)$ as a single object, assuming continuity over time.

Figure 1: Example of two estimations of the continuous curve at two different
locations $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ separated by $10$ km with
latitude fixed. We use $23$ cubic B-spline basis functions.

To illustrate the continuous estimation of the curves, we focus on data
observed on June $22$, 2010. For that specific day, we fit $23$ cubic B-spline
basis functions for each location, i.e.,
$\hat{y}(\mathbf{s}_{i};v)=\sum_{k=1}^{23}z_{k}(\mathbf{s}_{i})\eta_{k}(v)$,
where $\\{\eta_{k}\\}_{k=1}^{23}$ are B-spline basis functions. The continuous
curves are required to be smooth; thus, we consider a penalization term on the
second derivative of $\hat{y}(\mathbf{s}_{i};v)$ with respect to $v$.
Specifically, the estimator is such that

$y(\mathbf{s}_{i};v_{j})=\hat{y}(\mathbf{s}_{i};v_{j})+e_{ij},$

where the coefficients $z_{k}(\mathbf{s}_{i})$ of $\hat{y}(\mathbf{s}_{i};v)$
are obtained as the solution that minimizes

$\sum_{j=1}^{24}\\{y(\mathbf{s}_{i};v_{j})-\hat{y}(\mathbf{s}_{i};v_{j})\\}^{2}+\lambda\int\hat{y}^{\prime\prime}(\mathbf{s}_{i};v)\mathrm{d}v,$

with the smoothing parameter $\lambda$ fixed. We select the optimal smoothing
parameter by generalized cross-validation through the $n=529$ locations.
Figure 1 shows an example of the curve estimation at two locations. These two
locations are separated by $10$ km at the same latitude.

Once continuous curves are estimated with basis functions, we can compute the
mean (large-scale variation). Since data are measured on a regular grid of
points, and if each corresponding coefficient random field
$\\{Z_{k}(\mathbf{s}_{i})\\}_{i=1}^{n}$, with realization
$\\{z_{k}(\mathbf{s}_{i})\\}_{i=1}^{n}$, is stationary, then we can use the
empirical mean to estimate each $\mathbb{E}\\{Z_{k}(\mathbf{s})\\}$. Under
this scenario, Figure 2 (left) shows the result of the estimated mean curve.

On the other hand, if the coefficient random fields
$\\{Z_{k}(\mathbf{s}_{i})\\}_{i=1}^{n}$ are nonstationary, then
$\mathbb{E}\\{Z_{k}(\mathbf{s})\\}$ can depend on the location $\mathbf{s}$.
Under this framework, we use the surf.gls command defined in the spatial
package. Figure 2 (right) shows the estimator of
$\mathbb{E}\\{Z_{1}(\mathbf{s})\\}$ for the first random field coefficient
$Z_{1}(\mathbf{s})$. The surf.gls command uses the GLS method to obtain the
estimator with an exponential model as the covariance function. From Figure 2
(right), we observe some evidence that the mean depends on the location
$\mathbf{s}$.

Figure 2: Left: mean curve if the coefficient random fields
$\\{Z_{k}(\mathbf{s})\\}$ are stationary. Right: trend estimation of the first
coefficient random field $Z_{1}(\mathbf{s})$ if the coefficient random fields
$\\{Z_{k}(\mathbf{s})\\}$ are nonstationary.

Figure 3: Trace-variogram of the daily wind speed. Left: empirical. Right:
fitted model (blue).

For the modeling of small-scale variations, we use the trace-variogram defined
in (2.12). With the assumption of isotropy, Figure 3 (left) shows boxplots of
the corresponding empirical trace-variogram for different distances. The
fit.tracevariog command in the geofd package allows us to fit four covariance
models: spherical, exponential, Gaussian, and Matérn. In our data, the
exponential model is the best model in terms of minimizing the sum of squares
errors. Figure 3 (right) shows the fitted model. Finally, ordinary kriging
predictors can be obtained with the command okfd in the geofd package.

We have used daily wind speed to illustrate general results. Similarly, we can
use monthly or yearly wind speed with $y(\mathbf{s}_{i};v)$ representing the
entire year, and repeat the procedure.

###  7.2 Surface Data Approach

Figure 4: Estimation of surface data for $t=1$. Left: daily average of wind
speed on a regular grid of points. Right: the continuous surface estimated
with finite elements analysis.

For illustration purposes, we consider daily average data of the hourly wind
speed, that is, $y_{t}(\mathbf{s}_{i})$ represents the daily average wind
speed on day $t$, for $t=1,\ldots,365$ in the year $2010$.

To estimate the continuous surface $y_{t}(\mathbf{s})$ described in Section
5.3.1, we use the smooth.FEM.basis command in the fdaPDE package. It uses the
finite elements analysis and penalizes with the Laplacian. Figure 4 shows one
surface data for a specific day, $t=1$. In the left panel, we plot the surface
data observed on the $115\times 115$ km region centered in Dumah Al Jandal,
Saudi Arabia. In the right panel, we plot the estimated continuous surface.
The estimated continuous functions
$\\{\hat{y}_{1}(\mathbf{s}),\ldots,\hat{y}_{365}(\mathbf{s})\\}$ is the
surface time series, and can be modeled as in (5.4) and (5.5). This approach
allows us to successfully forecast the next day of surface wind data,
$\hat{x}_{T+1}(\mathbf{s})$ and $\hat{y}_{T+1}(\mathbf{s})$, as well as
understanding its temporal dependence.

Here, we consider a rectangular region domain centered at Dumah Al Jandal,
Saudi Arabia, but could also, similarly, consider the domain as being the
entire kingdom of Saudi Arabia or the entire world, and study the statistical
properties of the surface time series. In general, the approach of functional
data can be used in more complex and large datasets.

##  8 Discussion

In this paper, we have provided an overview of functional data analysis in the
case of spatially correlated data. We presented two main approaches, one when
data are curves observed over space (spatial functional data), and the other
when spatio-temporal random fields are considered as a surface (or manifold)
time series. These two approaches present a new paradigm of data analysis, in
which the continuous processes or random fields are considered as a single
entity. Although software packages are still limited for statistical analysis,
we believe that this mode of thinking can be valuable in the context of big
data.

It is a welcomed fact that users can have different tools for data analysis,
either with the surface time series approach or with classical spatio-temporal
techniques. The choice will depend on the data, on the phenomena being
studied, and on the scientific questions of interest. We appreciate all the
effort and work done in both directions. We have attempted to collect all
significant references, and we apologize for any unmentioned works.

Similarly, as functional data with temporal dependence arise, we can also
consider spatio-temporal functional random fields (Bel et al., 2011; Lee et
al., 2015) meaning that a functional time series is observed at each location.
For example, $y_{t}(\mathbf{s}_{i};v)$ can represent the functional time
series of daily wind speed at location $\mathbf{s}_{i}$, where $t$ is the day
index, and $v$ is the time within a day. Multivariate spatial functional data
can also be considered, i.e., at each location, we can observe different
functional data, such as temperature, precipitation, and humidity.

We conclude that the functional approach opens new areas of research to
develop methodologies and theories for the analysis of complex and large
spatio-temporal datasets.

## Acknowledgments

We thank Carolina Euan (KAUST), Joydeep Chowdhury (ISI), Laura Sangalli
(PoliMi), Soumya Das (KAUST), and Zhuo Qu (KAUST) for comments on this paper.
We are grateful to Professor Georgiy Stenchikov’s group, the Atmospheric and
Climate Modeling group at KAUST, for producing and providing the high-
resolution WRF dataset.

## References

  * Abdulah et al. (2019) Abdulah, S., Y. Li, J. Cao, H. Ltaief, D. E. Keyes, M. G. Genton, and Y. Sun (2019).  ExaGeoStatR: A package for large-scale geostatistics in R.  arXiv:1908.06936.
  * Aguilera-Morillo et al. (2017) Aguilera-Morillo, M. C., M. Durbán, and A. M. Aguilera (2017).  Prediction of functional data with spatial dependence: a penalized approach.  Stochastic Environmental Research and Risk Assessment 31(1), 7–22.
  * Alfeld et al. (1996) Alfeld, P., M. Neamtu, and L. L. Schumaker (1996).  Fitting scattered data on sphere-like surfaces using spherical splines.  Journal of Computational and Applied Mathematics 73(1), 5–43.
  * Arnone et al. (2019) Arnone, E., L. Azzimonti, F. Nobile, and L. M. Sangalli (2019).  Modeling spatially dependent functional data via regression with differential regularization.  Journal of Multivariate Analysis 170, 275–295.  Special Issue on Functional Data Analysis and Related Topics.
  * Aston et al. (2017) Aston, J. A. D., D. Pigoli, and S. Tavakoli (2017).  Tests for separability in nonparametric covariance operators of random surfaces.  The Annals of Statistics 45(4), 1431–1461.
  * Aue et al. (2015) Aue, A., D. D. Norinho, and S. Hörmann (2015).  On the prediction of stationary functional time series.  Journal of the American Statistical Association 110(509), 378–392.
  * Azzimonti et al. (2015) Azzimonti, L., L. M. Sangalli, P. Secchi, M. Domanin, and F. Nobile (2015).  Blood flow velocity field estimation via spatial regression with PDE penalization.  Journal of the American Statistical Association 110(511), 1057–1071.
  * Baladandayuthapani et al. (2008) Baladandayuthapani, V., B. K. Mallick, M. Y. Hong, J. R. Lupton, N. D. Turner, and R. J. Carroll (2008).  Bayesian hierarchical spatially correlated functional data analysis with application to colon carcinogenesis.  Biometrics 64(1), 64–73.
  * Banerjee et al. (2015) Banerjee, S., B. P. Carlin, and A. E. Gelfand (2015).  Hierarchical Modeling and Analysis for Spatial Data (Second ed.), Volume 135 of Monographs on Statistics and Applied Probability.  CRC Press, Boca Raton, FL.
  * Bel et al. (2011) Bel, L., A. Bar-Hen, R. Petit, and R. Cheddadi (2011).  Spatio-temporal functional regression on paleoecological data.  Journal of Applied Statistics 38(4), 695–704.
  * Bernardi et al. (2017) Bernardi, M. S., L. M. Sangalli, G. Mazza, and J. O. Ramsay (2017).  A penalized regression model for spatial functional data with application to the analysis of the production of waste in Venice province.  Stochastic Environmental Research and Risk Assessment 31(1), 23–38.
  * Bohorquez et al. (2017) Bohorquez, M., R. Giraldo, and J. Mateu (2017).  Multivariate functional random fields: prediction and optimal sampling.  Stochastic Environmental Research and Risk Assessment 31(1), 53–70.
  * Bosq (2000) Bosq, D. (2000).  Linear Processes in Function Spaces, Volume 149 of Lecture Notes in Statistics.  Springer-Verlag, New York.  Theory and Applications.
  * Caballero et al. (2013) Caballero, W., R. Giraldo, and J. Mateu (2013).  A universal kriging approach for spatial functional data.  Stochastic Environmental Research and Risk Assessment 27(7), 1553–1563.
  * Cardot et al. (1999) Cardot, H., F. Ferraty, and P. Sarda (1999).  Functional linear model.  Statistics & Probability Letters 45(1), 11–22.
  * Chen et al. (2017) Chen, K., X. Zhang, A. Petersen, and H.-G. Müller (2017).  Quantifying infinite-dimensional data: Functional data analysis in action.  Statistics in Biosciences 9(2), 582–604.
  * Crainiceanu et al. (2011) Crainiceanu, C. M., B. S. Caffo, S. Luo, V. M. Zipunnikov, and N. M. Punjabi (2011).  Population value decomposition, a framework for the analysis of image populations.  Journal of the American Statistical Association 106(495), 775–790.
  * Cressie and Wikle (2011) Cressie, N. and C. K. Wikle (2011).  Statistics for Spatio-Temporal Data.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ.
  * Cressie (2015) Cressie, N. A. C. (2015).  Statistics for Spatial Data (Revised ed.).  Wiley Classics Library. John Wiley & Sons, Inc., New York.  Paperback edition of the 1993 edition.
  * Dassi et al. (2015) Dassi, F., B. Ettinger, S. Perotto, and L. M. Sangalli (2015).  A mesh simplification strategy for a spatial regression analysis over the cortical surface of the brain.  Applied Numerical Mathematics. An IMACS Journal 90, 111–131.
  * Delicado et al. (2010) Delicado, P., R. Giraldo, C. Comas, and J. Mateu (2010).  Statistics for spatial functional data: some recent contributions.  Environmetrics 21(3-4), 224–239.
  * Diggle and Ribeiro (2007) Diggle, P. J. and P. J. Ribeiro, Jr. (2007).  Model-based Geostatistics.  Springer Series in Statistics. Springer, New York.
  * Duchamp and Stuetzle (2003) Duchamp, T. and W. Stuetzle (2003).  Spline smoothing on surfaces.  Journal of Computational and Graphical Statistics 12(2), 354–381.
  * Duchon (1977) Duchon, J. (1977).  Splines minimizing rotation-invariant semi-norms in sobolev spaces.  In W. Schempp and K. Zeller (Eds.), Constructive Theory of Functions of Several Variables, Berlin, Heidelberg, pp. 85–100. Springer Berlin Heidelberg.
  * Eilers and Marx (1996) Eilers, P. H. C. and B. D. Marx (1996).  Flexible smoothing with $B$-splines and penalties.  Statistical Science 11(2), 89–121.  With comments and a rejoinder by the authors.
  * Ettinger et al. (2016) Ettinger, B., S. Perotto, and L. M. Sangalli (2016).  Spatial regression models over two-dimensional manifolds.  Biometrika 103(1), 71–88.
  * Ferraty and Vieu (2006) Ferraty, F. and P. Vieu (2006).  Nonparametric Functional Data Analysis.  Springer Series in Statistics. Springer, New York.  Theory and Practice.
  * Finley et al. (2015) Finley, A. O., S. Banerjee, and A. E. Gelfand (2015).  spBayes for large univariate and multivariate point-referenced spatio-temporal data models.  Journal of Statistical Software 63(13), 1–28.
  * Galeano and Peña (2019) Galeano, P. and D. Peña (2019).  Data science, big data and statistics.  TEST 28(2), 289–329.
  * Genton et al. (2014) Genton, M. G., C. Johnson, K. Potter, G. Stenchikov, and Y. Sun (2014).  Surface boxplots.  Stat 3(1), 1–11.
  * Genton and Kleiber (2015) Genton, M. G. and W. Kleiber (2015).  Cross-covariance functions for multivariate geostatistics.  Statistical Science 30(2), 147–163.  With comments and a rejoinder by the authors.
  * Giraldo et al. (2018) Giraldo, R., S. Dabo-Niang, and S. Martínez (2018).  Statistical modeling of spatial big data: An approach from a functional data analysis perspective.  Statistics & Probability Letters 136, 126–129.  The role of Statistics in the era of big data.
  * Giraldo et al. (2010) Giraldo, R., P. Delicado, and J. Mateu (2010).  Continuous time-varying kriging for spatial prediction of functional data: An environmental application.  Journal of Agricultural, Biological, and Environmental Statistics 15(1), 66–82.
  * Giraldo et al. (2011) Giraldo, R., P. Delicado, and J. Mateu (2011).  Ordinary kriging for function-valued spatial data.  Environmental and Ecological Statistics 18(3), 411–426.
  * Giraldo et al. (2012) Giraldo, R., P. Delicado, and J. Mateu (2012).  Hierarchical clustering of spatially correlated functional data.  Statistica Neerlandica 66(4), 403–421.
  * Giraldo et al. (2015) Giraldo, R., P. Delicado, and J. Mateu (2015).  geofd: Spatial Prediction for Function Value Data.  R package version 1.0.
  * Gneiting (2013) Gneiting, T. (2013).  Strictly and non-strictly positive definite functions on spheres.  Bernoulli 19(4), 1327–1349.
  * Goulard and Voltz (1993) Goulard, M. and M. Voltz (1993).  Geostatistical Interpolation of Curves: A Case Study in Soil Science, pp. 805–816.  Dordrecht: Springer Netherlands.
  * Greco et al. (2018) Greco, F., M. Ventrucci, and E. Castelli (2018).  P-spline smoothing for spatial data collected worldwide.  Spatial Statistics 27, 1–17.
  * Gromenko et al. (2012) Gromenko, O., P. Kokoszka, L. Zhu, and J. Sojka (2012).  Estimation and testing for spatially indexed curves with application to ionospheric and magnetic field trends.  The Annals of Applied Statistics 6(2), 669–696.
  * Grujic and Menafoglio (2017) Grujic, O. and A. Menafoglio (2017).  fdagstat, an R package.  R package version 1.0.
  * Grujic et al. (2018) Grujic, O., A. Menafoglio, G. Yang, and J. Caers (2018, Jul).  Cokriging for multivariate hilbert space valued random fields: application to multi-fidelity computer code emulation.  Stochastic Environmental Research and Risk Assessment 32(7), 1955–1971.
  * Haining (2003) Haining, R. (2003).  Spatial Data Analysis: Theory and Practice.  Cambridge University Press.
  * Hall et al. (1994) Hall, P., N. I. Fisher, and B. Hoffmann (1994).  On the nonparametric estimation of covariance functions.  The Annals of Statistics 22(4), 2115–2134.
  * Hall and Patil (1994) Hall, P. and P. Patil (1994).  Properties of nonparametric estimators of autocovariance for stationary random fields.  Probability Theory and Related Fields 99(3), 399–424.
  * Hörmann and Kokoszka (2010) Hörmann, S. and P. Kokoszka (2010).  Weakly dependent functional data.  The Annals of Statistics 38(3), 1845–1884.
  * Horváth and Kokoszka (2012) Horváth, L. and P. Kokoszka (2012).  Inference for Functional Data with Applications.  Springer Series in Statistics. Springer, New York.
  * Horváth et al. (2013) Horváth, L., P. Kokoszka, and R. Reeder (2013).  Estimation of the mean of functional time series and a two-sample problem.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75(1), 103–122.
  * Horváth et al. (2014) Horváth, L., P. Kokoszka, and G. Rice (2014).  Testing stationarity of functional time series.  Journal of Econometrics 179(1), 66–82.
  * Hyndman and Ullah (2007) Hyndman, R. J. and M. S. Ullah (2007).  Robust forecasting of mortality and fertility rates: a functional data approach.  Computational Statistics & Data Analysis 51(10), 4942–4956.
  * Hörmann and Kokoszka (2012) Hörmann, S. and P. Kokoszka (2012).  Chap 7 - functional time series.  In T. S. Rao, S. S. Rao, and C. Rao (Eds.), Time Series Analysis: Methods and Applications, Volume 30 of Handbook of Statistics, pp. 157–186. Elsevier.
  * Ignaccolo et al. (2014) Ignaccolo, R., J. Mateu, and R. Giraldo (2014).  Kriging with external drift for functional data for air quality monitoring.  Stochastic Environmental Research and Risk Assessment 28(5), 1171–1186.
  * Jiang and Serban (2012) Jiang, H. and N. Serban (2012).  Clustering random curves under spatial interdependence with application to service accessibility.  Technometrics 54(2), 108–119.
  * Kokoszka and Reimherr (2013) Kokoszka, P. and M. Reimherr (2013).  Determining the order of the functional autoregressive model.  Journal of Time Series Analysis 34(1), 116–129.
  * Kokoszka and Reimherr (2017) Kokoszka, P. and M. Reimherr (2017).  Introduction to Functional Data Analysis.  Texts in Statistical Science Series. CRC Press, Boca Raton, FL.
  * Lai and Schumaker (2007) Lai, M.-J. and L. L. Schumaker (2007).  Spline Functions on Triangulations.  Encyclopedia of Mathematics and its Applications. Cambridge University Press.
  * Lee et al. (2015) Lee, D.-J., Z. Zhu, and P. Toscas (2015).  Spatio-temporal functional data analysis for wireless sensor networks data.  Environmetrics 26(5), 354–362.
  * Lila et al. (2016) Lila, E., J. A. D. Aston, and L. M. Sangalli (2016).  Smooth principal component analysis over two-dimensional manifolds with an application to neuroimaging.  The Annals of Applied Statistics 10(4), 1854–1879.
  * Lila et al. (2019) Lila, E., L. M. Sangalli, J. Ramsay, and L. Formaggia (2019).  fdaPDE: Functional Data Analysis and Partial Differential Equations; Statistical Analysis of Functional and Spatial Data, Based on Regression with Partial Differential Regularizations.  R package version 0.1-6.
  * Lindgren et al. (2011) Lindgren, F., H. Rue, and J. Lindström (2011).  An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73(4), 423–498.
  * Liu et al. (2017) Liu, C., S. Ray, and G. Hooker (2017).  Functional principal component analysis of spatially correlated data.  Statistics and Computing 27(6), 1639–1654.
  * Martínez-Hernández and Genton (2020) Martínez-Hernández, I. and M. G. Genton (2020).  Nonparametric trend estimation in functional time series with application to annual mortality rates.  arXiv:2001.04660.
  * Martínez-Hernández et al. (2019) Martínez-Hernández, I., M. G. Genton, and G. González-Farías (2019).  Robust depth-based estimation of the functional autoregressive model.  Computational Statistics & Data Analysis 131, 66–79.  High-dimensional and functional data analysis.
  * Mateu and Romano (2017) Mateu, J. and E. Romano (2017).  Advances in spatial functional statistics.  Stochastic Environmental Research and Risk Assessment 31(1), 1–6.
  * Menafoglio et al. (2018) Menafoglio, A., G. Gaetani, and P. Secchi (2018).  Random domain decompositions for object-oriented kriging over complex domains.  Stochastic Environmental Research and Risk Assessment 32(12), 3421–3437.
  * Menafoglio et al. (2016) Menafoglio, A., O. Grujic, and J. Caers (2016).  Universal kriging of functional data: Trace-variography vs cross-variography? Application to gas forecasting in unconventional shales.  Spatial Statistics 15, 39–55.
  * Menafoglio and Secchi (2017) Menafoglio, A. and P. Secchi (2017).  Statistical analysis of complex and spatially dependent data: A review of object oriented spatial statistics.  European Journal of Operational Research 258(2), 401–410.
  * Menafoglio et al. (2013) Menafoglio, A., P. Secchi, and M. Dalla Rosa (2013).  A universal kriging predictor for spatially dependent functional data of a Hilbert space.  Electronic Journal of Statistics 7, 2209–2240.
  * Morris et al. (2011) Morris, J. S., V. Baladandayuthapani, R. C. Herrick, P. Sanna, and H. Gutstein (2011).  Automated analysis of quantitative image data using isomorphic functional mixed models, with application to proteomics data.  The Annals of Applied Statistics 5(2A), 894–923.
  * Nerini et al. (2010) Nerini, D., P. Monestiez, and C. Manté (2010).  Cokriging for spatial functional data.  Journal of Multivariate Analysis 101(2), 409–418.  Statistical Methods and Problems in Infinite-dimensional Spaces.
  * Nychka et al. (2017) Nychka, D., R. Furrer, J. Paige, and S. Sain (2017).  fields: Tools for spatial data.  R package version 9.9.
  * Pebesma (2004) Pebesma, E. J. (2004).  Multivariable geostatistics in S: the gstat package.  Computers & Geosciences 30, 683–691.
  * Qingguo and Longsheng (2010) Qingguo, T. and C. Longsheng (2010).  B-spline estimation for spatial data.  Journal of Nonparametric Statistics 22(1-2), 197–217.
  * R Core Team (2019) R Core Team (2019).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Ramsay and Silverman (2005) Ramsay, J. O. and B. W. Silverman (2005).  Functional Data Analysis (Second ed.).  Springer Series in Statistics. Springer, New York.
  * Ramsay et al. (2018) Ramsay, J. O., H. Wickham, S. Graves, and G. Hooker (2018).  fda: Functional Data Analysis.  R package version 2.4.8.
  * Ramsay (2002) Ramsay, T. (2002).  Spline smoothing over difficult regions.  Journal of the Royal Statistical Society. Series B 64(2), 307–319.
  * Rekabdarkolaee et al. (2019) Rekabdarkolaee, H. M., C. Krut, M. Fuentes, and B. J. Reich (2019).  A Bayesian multivariate functional model with spatially varying coefficient approach for modeling hurricane track data.  Spatial Statistics 29, 351–365.
  * Reyes et al. (2015) Reyes, A., R. Giraldo, and J. Mateu (2015).  Residual kriging for functional spatial prediction of salinity curves.  Communications in Statistics. Theory and Methods 44(4), 798–809.
  * Ribeiro Jr and Diggle (2018) Ribeiro Jr, P. J. and P. J. Diggle (2018).  geoR: Analysis of Geostatistical Data.  R package version 1.7-5.2.1.
  * Romano et al. (2017) Romano, E., A. Balzanella, and R. Verde (2017).  Spatial variability clustering for spatially dependent functional data.  Statistics and Computing 27(3), 645–658.
  * Rue et al. (2009) Rue, H., S. Martino, and N. Chopin (2009).  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  Journal of the Royal Statistical Society. Series B 71(2), 319–392.
  * Ruiz-Medina (2011) Ruiz-Medina, M. (2011).  Spatial autoregressive and moving average hilbertian processes.  Journal of Multivariate Analysis 102(2), 292–305.
  * Ruiz-Medina (2012) Ruiz-Medina, M. (2012).  New challenges in spatial and spatiotemporal functional statistics for high-dimensional data.  Spatial Statistics 1, 82–91.
  * Ruiz-Medina et al. (2007) Ruiz-Medina, M. D., R. Salmerón, and J. M. Angulo (2007).  Kalman filtering from POP-based diagonalization of ARH(1).  Computational Statistics & Data Analysis 51(10), 4994–5008.
  * Ruppert et al. (2003) Ruppert, D., M. P. Wand, and R. J. Carroll (2003).  Semiparametric Regression, Volume 12 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, Cambridge.
  * Sangalli (2019) Sangalli, L. M. (2019).  A novel approach to the analysis of spatial and functional data over complex domains.  Quality Engineering 0(0), 1–10.
  * Sangalli et al. (2013) Sangalli, L. M., J. O. Ramsay, and T. O. Ramsay (2013).  Spatial spline regression models.  Journal of the Royal Statistical Society. Series B 75(4), 681–703.
  * Sartori and Torriani (2019) Sartori, I. and L. Torriani (2019).  Manifoldgstat, an R package.  R package.
  * Schabenberger and Gotway (2005) Schabenberger, O. and C. A. Gotway (2005).  Statistical Methods for Spatial Data Analysis.  Texts in Statistical Science Series. Chapman & Hall/CRC.
  * Schlather et al. (2019) Schlather, M., A. Malinowski, M. Oesting, D. Boecker, K. Strokorb, S. Engelke, J. Martini, F. Ballani, O. Moreva, J. Auel, P. J. Menck, S. Gross, U. Ober, P. Ribeiro, B. D. Ripley, R. Singleton, B. Pfaff, and R Core Team (2019).  RandomFields: Simulation and Analysis of Random Fields.  R package version 3.3.6.
  * Scott-Hayward et al. (2014) Scott-Hayward, L. A. S., M. L. Mackenzie, C. R. Donovan, C. G. Walker, and E. Ashe (2014).  Complex region spatial smoother (CReSS).  Journal of Computational and Graphical Statistics 23(2), 340–360.
  * Song and Mallick (2019) Song, J. J. and B. Mallick (2019).  Hierarchical Bayesian models for predicting spatially correlated curves.  Statistics 53(1), 196–209.
  * Staicu et al. (2010) Staicu, A.-M., C. M. Crainiceanu, and R. J. Carroll (2010).  Fast methods for spatially correlated multilevel functional data.  Biostatistics 11, 177–194.
  * Stein (1999) Stein, M. L. (1999).  Interpolation of Spatial Data. Some Theory for Kriging.  Springer Series in Statistics. Springer-Verlag, New York.
  * Sun and Genton (2011) Sun, Y. and M. G. Genton (2011).  Functional boxplots.  Journal of Computational and Graphical Statistics 20(2), 316–334.
  * Sun and Genton (2012) Sun, Y. and M. G. Genton (2012).  Adjusted functional boxplots for spatio-temporal data visualization and outlier detection.  Environmetrics 23(1), 54–64.
  * Venables and Ripley (2002) Venables, W. N. and B. D. Ripley (2002).  Modern Applied Statistics with S (Fourth ed.).  New York: Springer.  ISBN 0-387-95457-0.
  * Wahba (1981) Wahba, G. (1981).  Spline interpolation and smoothing on the sphere.  SIAM Journal on Scientific and Statistical Computing 2(1), 5–16.
  * Wang and Ranalli (2007) Wang, H. and M. G. Ranalli (2007).  Low-rank smoothing splines on complicated domains.  Biometrics 63(1), 209–217.
  * Wilhelm et al. (2016) Wilhelm, M., L. Dedè, L. M. Sangalli, and P. Wilhelm (2016).  IGS: an IsoGeometric approach for smoothing on surfaces.  Computer Methods in Applied Mechanics and Engineering 302, 70–89.
  * Wood (2017) Wood, S. (2017).  Generalized Additive Models: An Introduction with R (2 ed.).  Chapman and Hall/CRC.
  * Wood (2006) Wood, S. N. (2006).  Low-rank scale-invariant tensor product smooths for generalized additive mixed models.  Biometrics 62(4), 1025–1036.
  * Wood et al. (2008) Wood, S. N., M. V. Bravington, and S. L. Hedley (2008).  Soap film smoothing.  Journal of the Royal Statistical Society. Series B 70(5), 931–955.
  * Xiao et al. (2013) Xiao, L., Y. Li, and D. Ruppert (2013).  Fast bivariate $P$-splines: the sandwich smoother.  Journal of the Royal Statistical Society. Series B 75(3), 577–599.
  * Yip (2018) Yip, C. M. A. (2018).  Statistical characteristics and mapping of near-surface and elevated wind resources in the Middle East.  Ph.D. Thesis, King Abdullah University of Science and Technology.
  * Yue and Speckman (2010) Yue, Y. and P. L. Speckman (2010).  Nonstationary spatial Gaussian Markov random fields.  Journal of Computational and Graphical Statistics 19(1), 96–116.  With supplementary data available online.
  * Zhang et al. (2011) Zhang, J., M. K. Clayton, and P. A. Townsend (2011).  Functional concurrent linear regression model for spatial images.  Journal of Agricultural, Biological, and Environmental Statistics 16(1), 105–130.
  * Zhang et al. (2016) Zhang, L., V. Baladandayuthapani, H. Zhu, K. A. Baggerly, T. Majewski, B. A. Czerniak, and J. S. Morris (2016).  Functional CAR models for large spatially correlated functional datasets.  Journal of the American Statistical Association 111(514), 772–786.
  * Zhou et al. (2010) Zhou, L., J. Z. Huang, J. G. Martinez, A. Maity, V. Baladandayuthapani, and R. J. Carroll (2010).  Reduced rank mixed effects models for spatially correlated hierarchical functional data.  Journal of the American Statistical Association 105(489), 390–400.