§ INTRODUCTION

With the increase of multidimensional data availability and modern computing power, statistical models for spatial
and spatio-temporal data are developing at a rapid pace. Hence, there is a need for stable and reliable, yet updated and
efficient, software packages. In this section, we briefly discuss multidimensional data in climate and environmental studies as well as statistical software for space-time data.

§.§ Multidimensional data

Large multidimensional data sets often arise when climate and environmental phenomena are observed at the global
scale over extended periods. In climate studies, relevant physical variables are observed on a three-dimensional (3D)
spherical shell (the atmosphere) while time is the fourth dimension.
For instance, measurements are obtained by radiosondes flying from ground level up to the stratosphere [Fassò et al., 2014], by interferometric sensors aboard
satellites [Finazzi et al., 2018] or by laser-based methods, such as Light Detection and Ranging (LIDAR) [Negri et al., 2018].
In this context, statistical modelling of multidimensional data requires describing and exploiting the spatio-temporal correlation of the underlying phenomenon or data-generating process.
This is done using explanatory variables and multidimensional latent variables with covariance functions defined over a convenient spatio-temporal support.
When considering 3D$\times$T data (4D for brevity), covariance functions defined over the 4D support may be adopted.
However, these covariance functions often have a complex form [Porcu et al., 2018].
Moreover, when estimating the model parameters or making inferences,
very large covariance matrices (though they may be sparse) are implied.

In large climate and environmental applications, 4D data are rarely collected at high frequency in all spatial and temporal dimensions.
Often, only one dimension is sampled at high frequency while the remaining dimensions are sampled sparsely.
Radiosonde data, for instance, are sparse over the Earth's sphere, but they are dense along the vertical dimension, providing atmospheric profiles.
This suggests that handling all spatial dimensions equally (e.g. using a 3D covariance function) may not be the best option from a modelling or computational perspective,
and a data reduction technique may be useful instead. In this paper, the functional data analysis (FDA) approach
[Ramsay and Silverman, 2007] is adopted to model the relationship between measurements along the profile, while the remaining dimensions
are handled following the classic spatio-temporal data modelling approach using only 2D spatial covariance functions.

§.§ Statistical software

Various software programmes are available for considering data on a plane or in a two-dimensional (2D) Euclidean space.
The choice is more restricted when considering multidimensional or non-Euclidean spaces arising from atmospheric or remote sensing spatio-temporal data observed
on the surface of a sphere and over time.

For example, Figure <ref> depicts
the spatial locations of measurements collected globally in a single day through radio sounding, as discussed in Section <ref>.
Space is three-dimensional, and measurements are repeated over time at the same spatial locations over the Earth's surface but at different pressure values.

Radio sounding data example. Each dot represents the spatial location of a measurement taken by a radiosonde. Dots of the same colour belong to the same radiosonde. (Pressure axis not in scale).

The spBayes package [Finley et al., 2015] handles large spatio-temporal data sets, but space is only 2D. The documentation of the spacetime [Pebesma, 2012] and gstat [Pebesma and Heuvelink, 2016] packages does not explicitly address the multidimensional case, but, according to [Gasch et al., 2015], both packages have some capabilities to handle the 3D$\times$T.
However, we want to avoid working with 3D spatial covariance functions or sample spatio-temporal variograms.
Fixed rank kriging ([Cressie and Johannesson, 2008] [Cressie and Johannesson, 2008]) implemented in the R
package FRK [Zammit-Mangion, 2018] handles spatial and  spatio-temporal data both on the Euclidean plane and on the
surface of the sphere. FRK implements a set of tools for data gridding and basis function computation, resulting in
efficient dimension reduction, allowing it to handle large satellite data sets [Cressie, 2018].
It is based on a spatio-temporal random effects (SRE) model estimated by the
expectation-maximisation (EM) algorithm.
Recent extensions to FRK include the use of multi-resolution basis functions [Tzeng and Huang, 2018].

A second package based on SRE and the EM algorithm is D-STEM v1 [Finazzi and Fassò, 2014].
This package implements an efficient state-space approach for handling the temporal dimension and a heterotopic multivariate response approach that is useful when correlating heterogeneous networks [Fassò and Finazzi, 2011, Calculli et al., 2015].

D-STEM v1 has been successfully used in various medium-to-large applications, proving that the EM algorithm implementation, being mainly based on closed-form iterations, is quite stable.
These applications include air quality assessment in the metropolitan areas of Milan, Teheran and Beijing [Fassò, 2013, Taghavi-Shahri et al., 2019, Wan et al., 2020];
multivariate spatio-temporal modelling at the country and continental levels in Europe [Finazzi et al., 2013, Fassò et al., 2016];
time series clustering [Finazzi et al., 2015];
emulators of atmospheric dispersion modelling systems [Finazzi et al., 2019]; and
near real-time prediction of earthquake parameters [Finazzi, 2020].

A brief, non-exhaustive list of other models and/or software packages for advanced spatial data modelling is presented below, according to the principal technique, allowing the handling of large data sets.
In general, these techniques aim at avoiding the Cholesky decomposition of large and dense covariance matrices.

Some approaches, including FRK and D-STEM v1, leverage sparse variance–covariance matrices. Others exploit the sparsity of the precision matrix, thanks to a spatial Markovian assumption. This class includes the R packages LatticeKrig [Nychka et al., 2015, Nychka et al., 2016], INLA
[Blangiardo et al., 2013, Lindgren and Rue, 2015, Bivand et al., 2015, Rue et al., 2014] and the multi-resolution approximation approach of [Katzfuss, 2017],
which uses the predictive process and the state space representation [Jurek and Katzfuss, 2018] to model spatio-temporal data. Low-rank models are another
popular approach used by spBayes. Finally, the R package laGP [Gramacy, 2016], based on a machine learning approach, implements an efficient nearest neighbour prediction-oriented method.
[Heaton et al., 2018] develop an interesting spatial prediction competition considering a large data set and involving the above-mentioned approaches.

We observe that, although some of the software packages mentioned above consider both space and time, to the best of our knowledge, none of them handles a spatio-temporal FDA approach for data sets of the kind discussed in <ref>.

In this paper, we present D-STEM v2, which is a MATLAB package, extending D-STEM v1. The new version introduces modelling of functional data indexed over space and time.
Moreover, new complexity reduction techniques have been added for both model estimation and dynamic mapping, which are especially useful for large data sets.

The rest of the paper is organised as follows.
Section <ref> introduces the methodology adopted in this paper and, in particular, the data modelling approach and the complexity-reduction techniques.
Section <ref> describes the D-STEM v2 software in terms of the MATLAB classes used to define the data structure, model fitting and diagnostics and kriging.
This is followed by an illustration of the software use through two case studies.
The first one, discussed in Section <ref>, considers high-frequency spatio-temporal ozone data in Beijing.
The second one, in Section <ref>, considers modelling of global atmospheric temperature profiles and exploits the complexity-reduction capabilities of the new package. Finally, concluding remarks are provided in Section <ref>.

§ METHODOLOGY

This section discusses the methodology behind the modelling and the complexity-reduction techniques implemented in D-STEM v2
when dealing with functional space-time data sets. Moreover, model estimation, validation and dynamic kriging are briefly discussed.

§.§ Model equations

Let $\bm{s}=(s_{lat},s_{lon})^{ \top }$ be a generic spatial location on the
Earth's sphere, $\mathbb{S}^{2}$, and $t\in \mathbb{N} $ a discrete time index.
It is assumed that the function of interest, $f\left(  \bm{s},h,t\right )$, with domain $\mathcal{H}=\left[  h_{1},h_{2}\right]  \subset \mathbb{R}$, can be observed at any $\left(  \bm{s},t\right)  $ and $h\in\mathcal{H}$ through noisy measurements $y(\bm{s},h,t)$ according to the following model:

\begin{align}
y(\bm{s},h,t)  &  =f\left(  \bm{s},h,t\right)  +\varepsilon
f\left(  \bm{s},h,t\right)   &  =\bm{x}(\bm{s},h,t)^{ \top }\bm{\beta}\left(  h\right)  +\bm{\phi}(h)^{\top}\bm{z}(\bm{s}%
\bm{z}(\bm{s},t)  &  =\bm{G}\bm{z}(\bm{s}%
,t-1)+\bm{\eta}(\bm{s},t). \label{eq:model_line3}%
\end{align}

This model is referred to as the functional hidden dynamic geostatistical model (f-HDGM). In Equation (<ref>), $\varepsilon$ is a zero-mean Gaussian
measurement error independent in space and time with functional variance
$\sigma_{\varepsilon}^{2}\left(  h\right)  $, implying
that $\varepsilon$ is heteroskedastic across the domain $\mathcal{H}$. The variance is modelled as

\[
\log(\sigma_{\varepsilon}^{2}\left(  h\right)  )=\bm{\phi}(h)^{\top
\]

where $\bm{\phi}(h)$ is a $p\times1$ vector of basis functions evaluated at
$h$, while $\bm{c}_{\varepsilon}$ is a vector of coefficients to be estimated. In Equation (<ref>), $\bm{x}(\bm{s},h,t)$ is a $b \times 1$ vector of covariates while $\bm{\beta}\left(  h\right)  =(\beta_{1}(h),...,\beta_{b}(h))^{\top}$ is the vector of functional parameters modelled as

\[
\beta_{j}(h)=\bm{\phi}(h)^{\top}\bm{c}_{\beta,j}, j=1,...,b,
\]

and $\bm{c}_{\beta}=\left( \bm{c}_{\beta,1}^{\top},...,\bm{c}_{\beta,b}^{\top}\right) ^{\top}$ is a $pb\times1$  vector of coefficients that needs to be estimated. Additionally, $\bm{z}(\bm{s},t)$ is a $p\times1$ latent space-time
variable with Markovian dynamics given in Equation (<ref>).
The matrix $\bm{G}$ is a diagonal transition matrix with diagonal elements in the $p \times 1$ vector
$\bm{g}$. The innovation vector $\bm{\eta}$ is obtained from a
multivariate Gaussian process that is independent in time but correlated across space
with matrix spatial covariance function given by

\[
\bm{\Gamma}(\bm{s},\bm{s}^{\prime};\bm{\theta})=diag\left(
\rho(\bm{s},\bm{s}^{\prime};\bm{\theta}_{p})\right),
\]

where $\bm{v}=\left(  v_{1},...,v_{p}\right)  ^{\top}$ is a vector of variances and $\rho(\bm{s},\bm{s}^{\prime};\bm{\theta}_{j})$ is a valid spatial correlation function for locations $\bm{s},\bm{s}^{\prime}\in\mathbb{S}^{2}$, parametrised by $\bm{\theta}_{j}$, and $\bm{\theta}=(\bm{\theta}_{1},...,\bm{\theta}_{p})^{\top}$. The unknown model parameter vector is given by $\bm{\psi}=\left(
\bm{c}_{\varepsilon}^{\top},\bm{c}_{\beta}^{\top},\bm{g}^{\top},\bm{v}^{\top},\bm{\theta}^{\top}\right) ^{\top} $.

Note that, in order to ease the notation, the same $p$-dimensional basis functions $\bm{\phi}(h)$ are used to model
$\sigma_{\varepsilon}^{2}$, $\beta_{j}$ and $\bm{\phi}(h)^{\top}\bm{z}(\bm{s},t)$ in Equations
(<ref>)-(<ref>). In practice, D-STEM v2 allows one to specify a different number
of basis functions for each model component. Also note that $\varepsilon$ is not a pure measurement error since it
also accounts for model misspecification.
Finally, the covariates $\bm{x}(\bm{s},h,t)$ are assumed to be known without error for any $\bm{s}$, $h$ and $t$,
and thus they do not need a basis function representation.

§.§ Basis function choice

Choosing basis functions essentially means choosing the basis type and the number of basis functions. D-STEM v2 currently supports Fourier bases and B-spline bases. The former guarantee that the function is periodic in the domain $\mathcal{H}$, while the latter are not (in general) periodic
but have higher flexibility in describing functions with a complex shape.
Whichever basis function type is adopted, the number $p$ of basis functions must be fixed before model estimation. Usually, a high $p$ implies
a better model $R^2$, but over-fitting may be an issue. Moreover, special care must be taken when choosing the number of basis functions for
$\bm{\phi}(h)^{\top}\bm{z}(\bm{s},t)$. The classic FDA approach suggests fixing a high number of basis functions and adopting penalisation to avoid over-fitting.
In our context, this is not viable since the covariance matrices involved in model estimation have dimension $n^3p^3 \times n^3p^3$. Since $n$ is usually large, a large $p$ would make model estimation unfeasible, especially if the number of time points $T$ is also high.
When using B-spline basis, a small $p$ implies that the location of knots along the domain $\mathcal{H}$ also matters and may affect the model fitting performance. Ideally, $p$ and knot locations are chosen using a model validation technique (see <ref>) by trying different combinations of $p$ and knot locations.
If, due to time constraints, this is not possible, equally spaced knots are a convenient option.

§.§ Model estimation

The estimation of $\bm{\psi}$ and the latent space-time variable $\bm{z}(\bm{s},t)$ is based on the maximum likelihood approach considering profile data observed at spatial locations $\mathcal{S}=\{\bm{s}_{i},i=1,...,n\}$ and time points $t=1,...,T$.

At a specific location $\bm{s}_{i}$ and time $t$, $q_{i,t}$ measurements are taken at points $\bm{h}_{\bm{s}_{i},t}=\left(  h_{i,1,t},...,h_{i,q_{i,t},t}\right)  ^{\top}$ and collected in the vector

\[
\bm{y}_{\bm{s}_{i},t}=(y(\bm{s}_{i},h_{i,1,t},t),...,y(\bm{s}_{i},h_{i,q_{i,t},t},t))^{\top},
\]

here called the observed profile.

Although D-STEM v2 allows for varying $q_{i,t}$, for ease of notation, it is assumed here that
all profiles include exactly $q$ measurements, although $\bm{h}_{\bm{s}_{i},t}$ may be different across profiles.
Profiles observed at time $t$ across spatial locations $\mathcal{S}$ are then stored in the $nq\times1$
vector $\bm{y}_{t}=(\bm{y}_{s_{1},t}^{\top},...,\bm{y}_{s_{n},t}^{\top})^{\top}$.
Applying model (<ref>)-(<ref>) to the defined data above, we have the following matrix representation:

\begin{align*}
\bm{y}_{t}  &  =\tilde{\bm{X}}_{t}\bm{c}_{\bm{\beta}}+\bm{\Phi}_{\bm{z},t}\bm{z}%
\bm{z}_{t}  &  =\tilde{\bm{G}}\bm{z}_{t-1}+\bm{\eta}_{t},
\end{align*}

where $\tilde{\bm{X}}_{t}=\bm{X}_{t}\bm{\Phi}_{\bm{\beta},t}$ is a $nq\times bp$
matrix, with $\bm{X}_{t}$ the matrix of covariates and $\bm{\Phi}_{\bm{\beta},t}$
the basis matrix for $\bm{\beta}$. $\bm{\Phi}_{\bm{z},t}$ is the $nq\times np$ basis
matrix for the latent $np\times1$ vector
}(\bm{s}_{1},t)^{\top},...,\bm{\eta}(\bm{s}_{n},t)^{\top})^{\top}$ is the $np\times1$ innovation vector, while $\bm{\varepsilon}_{t}\ $ is the $nq\times1$ vector of measurement errors.
Additionally, $\tilde{\bm{G}}= \bm{I}_{n} \otimes \bm{G}$ is the $np\times np$
diagonal transition matrix.

The complete-data likelihood function $L(\bm{\psi};\bm{Y},\bm{Z})$ can
be written as

\[
\]

where $\bm{Y}=\left(  \bm{y}_{1},...,\bm{y}_{T}\right)  $,
$\bm{Z}=\left(  \bm{z}_{0},\bm{z}_{1},...,\bm{z}_{T}\right)
$, $\bm{\psi}_{\bm{z}}=\left(  \bm{g}^{\top},\bm{v}^{\top},\bm{\theta}^{\top}\right) ^{\top} $, $\bm{\psi}_{\bm{y}}=\left(
\bm{c}_{\varepsilon}^{\top},\bm{c}_{\beta}^{\top}\right) ^{\top} $, and $\bm{z}_{0}$ is the Gaussian initial vector with parameter $\bm{\psi}_{\bm{z}_{0}}$. Maximum likelihood estimation is based on an extension of the EM algorithm detailed in
[Calculli et al., 2015]. The model parameter set $\bm{\psi}$ is initialised
with starting values $\bm{\psi}^{\left\langle 0\right\rangle }$ and then
updated at each iteration $\iota$ of the EM algorithm.

The algorithm terminates if any of the following conditions is satisfied:

\[
\max_{l}\left\vert \psi _{l}^{\left\langle \iota \right\rangle }-\psi
_{l}^{\left\langle \iota -1\right\rangle }\right\vert /\left\vert \psi
_{l}^{\left\langle \iota \right\rangle }\right\vert <\epsilon _{1}
\]

\[ \left\vert
L(\bm{\psi}^{\left\langle \iota\right\rangle };\bm{Y}%
)-L(\bm{\psi}^{\left\langle \iota-1\right\rangle };\bm{Y})\right\vert
/\left\vert L(\bm{\psi}^{\left\langle \iota\right\rangle };\bm{Y})\right\vert

\[ \iota>\iota^{\ast}, \]

where $\psi _{l}^{\left\langle \iota \right\rangle }$ is the generic element of $\bm{\psi}^{\left\langle \iota \right\rangle }$ at the
$\iota\text{-}th$ iteration, $L(\bm{\psi}^{\left\langle \iota\right\rangle };\bm{Y})$ is the
observed-data likelihood function evaluated at $\bm{\psi}^{\left\langle
\iota\right\rangle }$, $0<\epsilon_1\ll1$ and $0<\epsilon_2\ll1$ are small positive numbers

(e.g. $10^{-4}%
$), while $\iota^{\ast}$ is a user-defined positive integer number (e.g. $100$) to limit the
iterations in the case of convergence failure of the EM algorithm.

Note that $\mathcal{S}$ is not time-varying, which means that spatial locations are fixed. This could be a limit in applications where spatial locations change for each $t$. On the other hand, missing profiles are allowed; that is, $\bm{y}_{\bm{s}_{i},t}$ may be a vector of $q$ missing values at some $t$. In the extreme case, a given spatial location $\bm{s}_{i}$ has only one profile over the entire period (if all the profiles are missing, the spatial location can be dropped from the data set).
<cit.> explains how the likelihood function of a state-space model changes in the case of a missing observation vector and how the EM estimation formulas are derived. Missing data handling in D-STEM v2 is based on the same approach.

§.§ Partitioning

At each iteration of the EM algorithm, the computational complexity of the E-step is
$O\left(  Tn^{3}p^{3}\right)  $, which may be unfeasible if $n$ is large.
When necessary, D-STEM v2 allows one to use a partitioning approach [Stein, 2013] for model
The spatial locations $\mathcal{S}$ are divided into $k$ partitions,
and $\bm{z}_{t}$ is partitioned conformably, namely, $\bm{z}_{t}=\left(  \bm{z}_{t}^{(1)\top},...,\bm{z}%
_{t}^{(k)\top}\right)  ^{\top}$.
Hence, the likelihood function becomes

\[
L\left(  \bm{\psi}_{\bm{y}};\bm{y}_{t}\mid\bm{z}_{t}\right)  \cdot%
L\left(  \bm{\psi}_{\bm{z}_{0}};\bm{z}_{0}^{(j)}\right)  \cdot%
L\left(  \bm{\psi}_{\bm{z}};\bm{z}_{t}^{(j)}\mid\bm{z}_{t-1}%
\]

From the EM algorithm point of view, this implies
that the E-step is independently applied to each partition, possibly in parallel. When all partitions are equal in size,
the computational complexity reduces to $\mathcal{O}\left( Tkr^{3}p^{3}\right)$, with $r$ as the partition size.

Geographical partitioning, constructed aggregating proximal locations, is a natural choice for environmental applications.
Given the number of partitions $k$, the k-means algorithm applied to spatial coordinates
provides a geographical partitioning of $\mathcal{S}$.
However, the number of points in each partition is not controlled, and a heterogeneous partitioning may arise.
If some subsets are very large and others are small, the reduction in computational complexity given above is far from being achieved.
This can easily happen, for example, when $\mathcal{S}$ is a global network constrained by continent shapes.

For this reason, D-STEM v2 provides a heuristically modified k-means algorithm that encourages
partitions with similar numbers of elements.
The algorithm optimises the following objective function:

\begin{equation}
\sum_{j=1}^{k}\sum_{\bm{s}\in\mathcal{S}_{j}}d\left(  \bm{s},\bm{c}_{j}\right)  +\lambda\sum_{j=1}^{k}\left(
r_{j}-\frac{n}{k}\right)  ^{2}, \label{eq:k-means}%
\end{equation}

where $\lambda\ge0$, $\mathcal{S}_{j} \subset \mathcal{S}$ is the set of coordinates in the $j\text{-}th$
partition, $d$ is the geodesic distance on the sphere $\mathbb{S}^{2}$ and $\bm{c}_{j}$ and $r_{j}$
are the centroid and the number of elements in the $j\text{-}th$ partition, respectively.

The second term in (<ref>) accounts for the variability of the partition sizes and acts as a penalisation for heterogeneous partitionings.
Clearly, when $\lambda=0$, the above-mentioned objective function gives the classic k-means algorithm.
For high values of $\lambda$, solutions with similarly sized partitions are favoured.

Unfortunately, an optimality theory for this algorithm has not yet been developed, and the choice of $\lambda$ is left to the user. Nonetheless, it may be a useful tool to define a partitioning that is appropriate for the application at hand with regard to computing time and geographical properties.

§.§ Variance-covariance matrix estimation

The EM algorithm provides a point estimate of the parameter vector $\bm{\psi}$
but no uncertainty information. Building on <cit.>, D-STEM v2
estimates the variance–covariance matrix
by means of the observed Fisher information matrix, $\mathbf{I}_{T}$, namely

To understand its computational cost, note that the information matrix given above may be written as a sum:

For large data sets, each matrix $\mathbf{i}_t$ may be expensive to compute, and the total computational cost is linear in $T$, provided missing data are evenly distributed in time.
This results in a time-consuming task with a computational burden even higher than that for model estimation.
For this reason, D-STEM v2 makes it possible to approximate $\hat{\Sigma}_{\bm{\psi},T}$ using a truncated information matrix, namely:

\begin{equation}
\tilde{\Sigma}_{\bm{\psi},t^*}=(\frac{T}{t^*}\mathbf{I}_{t^*})^{-1},
\label{eq:Fisher_approximated}
\end{equation}

which reduces the computational burden by a factor of $1-t^*/T$.

$\tilde{\Sigma}_{\bm{\psi},t^*} \rightarrow \hat{\Sigma}_{\bm{\psi},T}$ for $t^* \rightarrow T$,
the truncation time $t^*$ is chosen to control the approximation error in $\hat{\Sigma}_{\bm{\psi}}$. In particular, $t^*$ is the first integer such that

\begin{equation}
\frac{\left\Vert \tilde{\Sigma}_{\bm{\psi},t}-\tilde{\Sigma}_{\bm{\psi},t-1}%
\right\Vert_{F} }{\left\Vert \tilde{\Sigma}_{\bm{\psi},t}\right\Vert_{F} }\leq
\delta,\label{eq:varcov_approximated}
\end{equation}

where $\left\Vert { \cdot }\right\Vert_{F}$ is the Frobenius norm, and $\delta$ may be defined by the user.

Generally speaking, the behaviour of $\hat{\Sigma}_{\bm{\psi},T}$ for large $T$ and, hence, the behaviour of $\tilde{\Sigma}_{\bm{\psi},t}$ relays on stationarity and ergodicity of the underlying stochastic process; see, for example,
<cit.> and references therein.

To have operative guidance for the user, let us assume first that no missing values are present, the information matrix is well-conditioned and the covariates have no isolated outliers or extreme trends.
In this case, away from the borders $t\cong1$ and $t\cong T$,
the observed conditional information $\mathbf{i}_t$ has a relatively smooth stochastic behaviour, and the approximation in (<ref>) is expected to be satisfactory at the level defined by $\delta$.
Conversely, if some data are missing at time $t$, the information $\mathbf{i}_t$ is reduced accordingly. If the missing pattern is random over time, this is not an issue.
But, in the unfavourable case with a high percentage of missing data mostly concentrated at the end the time series, $t\cong T$, the above approximation may over-estimate the information and under-estimate the variances of the parameter estimates.

§.§ Dynamic kriging

In this paper, dynamic kriging refers to evaluating the following

\begin{align}
\hat{f} \left(  \bm{s},h,t\right) &= \mathbb{E}_{\hat{\bm{\psi}}}\left(  f\left(  \bm{s},h,t\right)  \mid
\bm{Y}\right), \label{eq:krig_exp}\\
\VAR\left( \hat{f} \left(  \bm{s},h,t\right) \right) &= \mathbb{V}_{\hat{\bm{\psi}}}\left(  f\left(  \bm{s},h,t\right)  \mid\bm{Y}\right),  \label{eq:krig_var}%
\end{align}

for any $\bm{s}\in\mathbb{S}^{2}$, $h\in\mathcal{H}$ and $t=1,...,T$. A common approach is to map the kriging estimates on a regular pixelation $\mathcal{S}^{\ast}=\left\{  \bm{s}_{1}^{\ast},...,\bm{s}_{m}^{\ast}\right\}  $. This may be a time-consuming task when $m$ and/or $n$ and/or $T$ are large. To tackle this problem, D-STEM v2 allows one to exploit a nearest-neighbour approach, where the conditioning term in Equations (<ref>) and
(<ref>) is not $\bm{Y}$, but the data at the spatial locations $\mathcal{S}_{\sim j}$, where
$\mathcal{S}_{\sim j}\subset\mathcal{S}$ is the set of the $\tilde{n}\ll n$ nearest
spatial locations to $\bm{s}_{j}^{\ast}$. The use of the nearest-neighbour approach
is justified by the so-called screening effect. Even when the spatial correlation function exhibits
long-range dependence, it can subsequently be assumed that $y$ at spatial location $\bm{s}$ is nearly independent of spatially distant
observations when conditioned on nearby observations <cit.>.

For computational efficiency, D-STEM v2 performs kriging for blocks of pixels.
To do this, $\mathcal{S}^{\ast}$ is partitioned in $u$ blocks $\mathcal{S}^{\ast}=\left\{\mathcal{S}_{1}^{\ast},...,\mathcal{S}_{u}^{\ast}\right\}$, and kriging is
done on each block $\mathcal{S}_{l}^{\ast}$, $l=1,...,u$, with $u\ll m$ controlled by the user.
For each target block $\mathcal{S}_{l}^{\ast}$, the conditioning term in Equations
(<ref>) and (<ref>) is given by the data
observed at $\mathcal{\tilde{S}}_{l}=\bigcup\nolimits_{j\in \mathcal{J}_{l}}\mathcal{S}_{\sim j},%
\mathcal{J}_{l}=\left\{ j:s_{j}^{\ast }\in \mathcal{S}_{l}^{\ast}\right\}$.
Note that, if $\mathcal{S}_{l}^{\ast}$ is dense and $\mathcal{S}$ is sparse (namely $n\ll m$), then $\mathcal{\tilde{S}}_{l}$ is not
much larger than $\mathcal{S}_{\sim j}$ since most of the spatial locations in $\mathcal{S}_{l}^{\ast}$
tend to have the same neighbours $\mathcal{S}_{\sim j}$.

§.§ Validation

D-STEM v2 allows one to implement an out-of-sample validation by partitioning the original spatial locations
$\mathcal{S}$ into subsets $\mathcal{S}_{est}$ and $\mathcal{S}_{val}$. Data at $\mathcal{S}%
_{est}$ are used for model estimation while data at $\mathcal{S}%
_{val}$ are used for validation. Once the model is estimated, the kriging
formula in Equation (<ref>) is used to predict at $\mathcal{S}_{val}$ for all times $t$ and heights $\bm{h}$. The following validation mean squared errors are
then computed

\begin{align*}
MSE_{t}  &  =\frac{1}{P_{1}}\sum_{\bm{s}\in\mathcal{S}_{val}}%
\sum_{h\in\bm{h}_{\bm{s},t}}\left(  y\left(  \bm{s},h,t\right)  -\hat{y}\left(  \bm{s},h,t\right)  \right)  ^{2},\\
MSE_{\bm{s}}  &  =\frac{1}{P_{2}}\sum_{t=1}^{T}\sum_{h\in
\bm{h}_{\bm{s},t}}\left(  y\left(  \bm{s},h,t\right)
-\hat{y}\left(  \bm{s},h,t\right)  \right)  ^{2},\\
MSE_{h}  &  =\frac{1}{P_{3}}\sum_{t=1}^{T}\sum_{\bm{s}\in\mathcal{S}_{val}}\left(  y\left(  \bm{s},h,t\right)
-\hat{y}\left(  \bm{s},h,t\right)  \right)  ^{2},%
\end{align*}

where $\hat{y}\left(  \bm{s},h,t\right)  $ is obtained from Equation (<ref>), while $P_{1}$, $P_{2}$ and $P_{3}$ are the number of
terms in each sum.

When $\bm{h}_{\bm{s},t}$ varies across the profiles, D-STEM v2 provides a binned MSE by splitting the continuous domain $\mathcal{H}$ into $B$ equally spaced intervals. Let $H^*_r$ be the set of observation points in the $r\text{-}th$ interval, let $n_r$ be the corresponding observation number and let $\bar{h}_r = \frac{1}{n_r} \sum_{h\in H^*_r} h$ be the mean of points in $b\text{-}th$ interval. Then, the $MSE_{\bar{h}_r}$ is computed by

\begin{align*}
MSE_{\bar{h}_r}  &  =\frac{1}{P_{4}} \sum_{h\in H^*_r} \sum_{t=1}^{T}\sum_{\bm{s}\in\mathcal{S}_{val}}\left(  y\left(  \bm{s},h,t\right)
-\hat{y}\left(  \bm{s},h,t\right)  \right)  ^{2},%
\end{align*}

where $P_{4}$ is the total number of observations in the $b\text{-}th$ interval.

D-STEM v2 also provides the validation $R^2$ with respect to time

\begin{align*}
R^{2}_{t}  &  =1 - \frac{MSE_{t}}{\VAR\left( \{y\left(  \bm{s},h,t\right), \bm{s}\in\mathcal{S}_{val}, h\in\bm{h}_{\bm{s},t} \} \right)}.
\end{align*}

and the analogous validation $R^2$ with respect to location $\bm{s}$ and $h_r$.

§ SOFTWARE

This section starts by briefly describing the modelling capabilities of D-STEM v2 inherited by the previous version for dealing with spatio-temporal data sets. Then, it focuses on the D-STEM v2 classes and methods, which implement estimation, validation and dynamic mapping of the model presented in Section <ref>. Although some of the classes are already available in D-STEM v1, they are listed here for completeness.

§.§ Software description

D-STEM v1 implemented a substantial number of models. The dynamic coregionalisation model (DCM, [Finazzi and Fassò, 2014] [Finazzi and Fassò, 2014]) and the hidden dynamic geostatistical model (HDGM, [Calculli et al., 2015] [Calculli et al., 2015])
are suitable for modelling and mapping multivariate space-time data collected from
unbalanced monitoring networks. Model-based clustering (MBC, [Finazzi et al., 2015] [Finazzi et al., 2015])  has been introduced
for clustering time series, and it is suitable for large data sets with spatially registered time series.
Moreover, the emulator model [Finazzi et al., 2019] is based on a Gaussian emulator, and it is exploited for
modelling the multivariate output of a complex physical model.

In addition, D-STEM v2 (available at <github.com/graspa-group/d-stem>) provides the functional version of HDGM, denoted by f-HDGM, which handles modelling and mapping of functional space-time data, following the methodology of Section <ref>. For implementing f-HDGM,
D-STEM v2 relies on the MATLAB version of the fda package [Ramsay et al., 2018], which is automatically downloaded and installed by D-STEM v2.

§.§ Data format

Two data formats are available to define observations for the f-HDGM. One is the internal format used by the D-STEM v2 classes, and the other one is the user format based on the more user-friendly table data type implemented in recent versions of MATLAB.
The latter permits storing measurement profiles, covariate profiles, coordinates, timestamps and units of measure in a single object. The internal format is not discussed here.

Considering a table in the user format, each row includes the profiles collected at a given spatial location and time point.
The column labels are defined as follows:
columns Y and Y_name are used for the dependent variable $y$ and its name as a string field, respectively;
the column with prefix X_h_ is used for the values of the domain $h$; eventually, columns with prefix X_beta_ are used for covariates $\bm{x}$.
These tables have only one column for $y$ and only one column for $h$. Instead, we can have any number $b \geqslant 0$ of covariate columns. Additionally, the table has columns X_coordinate and Y_coordinate for  spatial location $\bm{s}$ and column Time for the timestamp. Units of measure are stored in the Properties.VariableUnits property of the table columns and used in outputs and plots. Units for X_coordinate and Y_coordinate can be deg for degrees, m for meters and km for kilometres. Geodetic distance is used when the unit is deg; otherwise, the Euclidean distance is used.

At the table row corresponding to location $\bm{s}_i$ and time $t$, the elements related to $y$ and $\bm{x}$ are vectors with $q_{i,t}$ elements.
Vectors related to $y$ may include missing data (NaN). If $y$ is entirely missing for a
given $\left(  \bm{s},t\right)  $, the row must be removed from the table.
Since spatial locations $\mathcal{S}$ are fixed in time, and as their number $n$ is
determined by the number of unique coordinates in the table,
profiles observed at different time points but the same spatial location $\bm{s}$ must have
the same coordinates.

§.§ Software structure

In D-STEM, a hierarchical structure of object classes and methods is used to handle data definition, model definition and estimation, validation, dynamic kriging and the related plotting capabilities.
The structure is schematically given below.
Further details on the use of each class are given within the two case studies in this paper, while class constructors, methods and property details can be obtained in MATLAB using the command

doc <class_name>.

§.§.§ Data handling

The stem_data class allows the user to define the data used in f-HDGM models, mainly through the following objects and methods.

* Objects of stem_data

* stem_modeltype: model type (DCM, HDGM, MBC, Emulator or f-HDGM); note that model type is needed here because the data structure varies among the different models;

* stem_fda: basis functions specification;

* stem_validation (optional): definition of the learning and testing datasets for model validation.

* Methods and Properties of stem_data

* kmeans_partitioning: data partitioning for parallel EM computations of Section <ref>;
this method is applied to a stem_data object, and its output is used by the EM_estimation method in the stem_model class below;

* shape (optional): structure with geographical borders used for mapping.

* Internal Objects of stem_data

* stem_varset: observed data and covariates;

* stem_gridlist: list of stem_grid objects

* stem_grid: spatial locations coordinates;

* stem_datestamp: temporal information.

Interestingly, stem_misc.data_formatter is a helper method, which is useful for building stem_varset objects starting from data tables. Its class,  stem_misc, is a miscellanea static class implementing other methods for various intermediate tasks not discussed here for brevity.

§.§.§ Model building

The stem_model class is used to define, estimate, validate and output a f-HDGM, mainly through the following objects and methods.

* Objects of stem_model

* stem_data: defined above;

* stem_par: model parameters;

* stem_EM_result: container of the estimation output, after EM_estimate;

* stem_validation_result (optional): container of validation output, available only if stem_data contains the stem_validation object;

* stem_EM_options (optional): model estimation options; it is an input of the  EM_estimate method below.

* Methods of stem_model

* EM_estimate: computation of parameter estimates;

* set_varcov: computation of the estimated variance-covariance matrix;

* plot_profile: plot of functional data;

* print: print estimated model summary;

* beta_Chi2_test: testing significance of covariates;

* plot_par: plot functional parameter;

* plot_validation: plot MSE validation.

§.§.§ Kriging

The kriging handling is implemented with two classes. The first is the stem_krig class, which implements the kriging spatial interpolation.

* Objects of stem_krig

* stem_krig_data: mesh data for kriging;

* stem_krig_options: kriging options;

* Methods of stem_krig

* kriging: computation of kriging, the output is a stem_krig_result object.

The second is the stem_krig_result class, which stores the kriging output and implements the methods for plotting the kriging output.

* Methods of stem_krig_result

* surface_plot: mapping of kriging estimate and their standard deviation for fixed $h$;

* profile_plot: method for plotting the kriging function and the variance-covariance matrix for a fixed space and time.

Although at first reading the user could prefer a single object for both input and output of the kriging, these objects may be quite large, making the current approach more flexible.

§ CASE STUDY ON OZONE DATA

This section illustrates how to make inferences on an f-HDGM for ground-level high-frequency air
quality data collected by a monitoring network. In particular, hourly ozone ($O_{3}$, in $\mu g/m^3$) measured
in Beijing, China, is considered.

§.§ Air quality data

Ground-level $O_{3}$ is an increasing public concern due to its
essential role in air pollution and climate change. In China, $O_{3}$ has
become one of the most severe air pollutants in recent years
[Wang et al., 2017].

In this case study, the aim is to model hourly $O_{3}$ concentrations from 2015 to 2017 with respect to temperature and ultraviolet radiation (UVB) across Beijing.
Concentration and temperature data are available at twelve monitoring stations (Figure <ref>). Hourly UVB data are obtained from the ERA-Interim product of the European Centre for Medium-Range Weather Forecasts (ECMWF) at a grid size of $0.25^{\circ} \times 0.25^{\circ}$ over the city.

locations of the twelve stations in Beijing [Kahle and Wickham, 2013].

To describe the diurnal cycle of $O_{3}$, which peaks in the afternoon and reaches a minimum at night-time, the 24 hours of the day are used as domain $\mathcal{H}$ of the basis functions, while the time index $t$ is on the daily scale. Moreover, due to the circularity of
time, Fourier basis functions are adopted, which implies that $\beta_{j}\left(
h\right)  $, $\sigma_{\varepsilon}^{2}\left(  h\right)  $ are periodic functions.

The measurement equation for $O_{3}$ is

\begin{equation}
_{temp}(h)+x_{uvb}\left(  t\right)  \beta_{uvb}(h)+\bm{\phi}(h)^{\top}\bm{z}(\bm{s},t)+\varepsilon(\bm{s},h,t), \label{eq:O3meaurement}%
\end{equation}

where $\bm{s}$ is the generic spatial location, $h\in\left[  0,24\right)
$ is the time within the day expressed in hours and $t=1,...,1096$ is the day
index over the period 2015–2017. Based on a preliminary analysis, the number of basis functions for $\beta_{j}\left(  h\right)  $, $\sigma_{\varepsilon}^{2}\left(  h\right)  $
and $\bm{\phi}(h)^{\top}\bm{z}(\bm{s},t)$ is chosen to be $5$, $5$
and $7$, respectively.

§.§ Software implementation

This paragraph details the implementation of the D-STEM v2 in three
aspects: model estimation, validation and kriging. Relevant
scripts are demo_section4_model_estimate.m,
demo_section4_validation.m and demo_section4_kriging.m,
respectively, which are available in the supplementary material. All the scripts can be executed by choosing the option number from 1 to 3 in the demo_menu_user.m script.

§.§.§ Model estimation

This paragraph describes the demo_section4_model_estimate.m script
devoted to the estimation of the model parameters and of their
variance–covariance matrix.

The data set needed to perform this case study is stored as a MATLAB table in the user format of Section <ref> and named Beijing_O3. It can be loaded from the corresponding file as follows:

load ../Data/Beijing_O3.mat;

In the Beijing_O3 table, each row refers to a fixed space-time point and gives a 24-element hourly ozone profile with the corresponding conformable covariates, which are: a constant, temperature and UVB.

The following lines of code specify the model type and the basis functions, which are stored in an object of class stem_fda:

o_modeltype = stem_modeltype('f-HDGM');
input_fda.spline_type = 'Fourier';
input_fda.spline_range = [0 24];
input_fda.spline_nbasis_z = 7;
input_fda.spline_nbasis_beta = 5;
input_fda.spline_nbasis_sigma = 5;
o_fda = stem_fda(input_fda);

When using a Fourier basis, spline_nbasis_z must be set to a positive odd
Meanwhile, spline_nbasis_beta and/or spline_nbasis_sigma must be left empty, if $\bm{\beta}(h)\equiv\bm{\beta}$ and/or $\sigma
_{\varepsilon}^{2}(h)\equiv\sigma_{\varepsilon}^{2}$ are constant functions.

The next step is to define an object of class stem_data, which specifies the model type and contains the basis function object and the data from the Beijing_O3 table, transformed in the internal data format. This is done using the intermediate input_data structure:

input_data.stem_modeltype = o_modeltype;

input_data.data_table = Beijing_O3;

input_data.stem_fda = o_fda;
o_data = stem_data(input_data);

Then, an object of class stem_model is created by using both information on data, stored in the o_data object, and on parametrisation, contained in the stem_par object named o_par:

o_par = stem_par(o_data,'exponential');
o_model = stem_model(o_data, o_par);

To facilitate visualisation, the method plot_profile of class stem_model shows the $O_3$ profile data at location (lat0, lon0), in the days between t_start and t_end (Figure <ref>):

lat0 = 40; lon0 = 116;
t_start = 880; t_end = 900;
o_model.plot_profile(lat0, lon0, t_start, t_end);

Before running the EM algorithm, the model parameters need to be initialised.
This is done using the method get_beta0 of class stem_model, which
provides the starting values for $\bm{\beta}$, and the method
get_coe_log_sigma_eps0 for the case of a functional $\sigma_{\varepsilon}^{2}(h)$.
Next, the method set_initial_values of the o_model object is called to complete the initialisation of model parameters:

n_basis = o_fda.get_basis_number;
o_par.beta = o_model.get_beta0();
o_par.sigma_eps = o_model.get_coe_log_sigma_eps0();
o_par.theta_z = ones(1, n_basis.z)*0.18;
o_par.G = eye(n_basis.z)*0.5;
o_par.v_z = eye(n_basis.z)*10;

Note that the theta_z parameter must be provided in the same unit of measure as the spatial coordinates.

$O_3$ concentrations at location 39.92 latitude and 116.19 longitude for 21 days beginning on 29 May 2017. Left: each dot is a concentration measurement. The colour of the dot depicts the concentration. Right: each graph is a daily concentration profile.

Before model estimation, EM exiting conditions $\epsilon_1$ (exit_toll_par), $\epsilon_2$ (exit_toll_loglike) and $\iota^{\ast}$ (max_iterations) introduced in Section <ref> can be optionally defined as follows:

o_EM_options = stem_EM_options();
o_EM_options.exit_toll_par = 0.0001;
o_EM_options.exit_toll_loglike = 0.0001;
o_EM_options.max_iterations = 200;

Model estimation is started by calling the method EM_estimate of the
o_model object, with the optional o_EM_options object passed as an input argument.
After model estimation, the variance–covariance matrix of the estimated parameters is evaluated by calling the method set_varcov, with the optional approximation level $\delta$ of Equation (<ref>) passed as an input parameter.
Finally, set_logL computes the observed data log-likelihood.

delta = 0.001;

All the relevant estimation results are found in the internal
stem_EM_result object, which can be accessed as a property of the
o_model object as follows:

Figure <ref> is produced by calling the plot_par method and
shows the estimated $\bm{\beta}_{0}(h)$, $\bm{\beta}_{temp}(h)$, $\bm{\beta}_{uvb}(h)$, and $\sigma_{\varepsilon}^{2}(h)$.
Thanks to the use of a Fourier basis, the functions are periodic with a period of one day. In the plot
of $\sigma_{\varepsilon}^{2}(h)$, the unexplained portion of $O_{3}$ variance,
$\sigma_{\varepsilon}^{2}(h)$, is small during daylight hours, which is
consistent with the results of [Dohan and Masschelein, 1987].

When the confidence bands of parplot contain zero, it
may be useful to test the significance of the covariates. By calling the method beta_Chi2_test,
the results of $\chi^{2}$ tests are obtained, and they are reported in Table <ref>. Although $\beta_{uvb}$ is
close to $0$ in the morning, all fixed effects are highly significant overall. The model output is shown in the
MATLAB command window by calling the print method.

$\beta_{0}(h),~\beta_{temp}(h),~\beta_{uvb}(h)$ and
$\sigma_{\epsilon}^{2}(h)$, with $90\%,~95\%,~99\%$ confidence bands, respectively, shown through the different shades.

$\chi^{2}$ statistic     $p$ value
Constant     1r136.33     1r0
Temperature     1r14266.07     1r0
UVB     1r2094.34     1r0

$\chi^{2}$ tests for significance of covariates.

§.§.§ Validation

This paragraph describes the script demo_section4_validation.m, which implements validation. Compared
to the code in demo_section4_model_estimate.m, it only differs in
providing an object of class stem_validation.

To create the object called o_validation, the name of the validation
variable is needed as well as the indices of the validation stations.  Moreover, if the size of the nearest neighbour
set for each kriging site (nn_size) is not provided as the third input
argument in the stem_validation class constructor, D-STEM v2 uses all the
remaining stations. For example, a validation data set with three stations is constructed as follows:

S_val = [1,7,10];
input_data.stem_validation = stem_validation('O3', S_val);

The validation statistics, computed by EM_estimate, are saved in the internal object
stem_validation_result, which can be accessed as a property of the o_model object. The stem_validation_result object contains the estimated $O_{3}$ residuals for the above-mentioned validation stations as well as the validation mean square errors and $R^2$, as defined in Section <ref>.

§.§.§ Kriging

This paragraph describes the demo_section4_kriging.m script, which applies the
approach of Section <ref> to
the estimated model to map the $O_{3}$ concentrations over Beijing city.

The first step is to create an object of class
stem_grid, which collects the information about the regular grid of
pixels $\mathcal{B}$ to be used for mapping.
Then, an object of class stem_krig_data is created, where the o_krig_grid object is passed as an input argument:

load ../Output/ozone_model;

step_deg = 0.05;
lat = 39.4:step_deg:41.1;
lon = 115.4:step_deg:117.5;
[lon_mat,lat_mat] = meshgrid(lon,lat);

krig_coordinates = [lat_mat(:) lon_mat(:)];
o_krig_grid = stem_grid(krig_coordinates,'deg','regular','pixel',...
o_krig_data = stem_krig_data(o_krig_grid);

Two comments on the above lines follow. First, since the grid in the o_krig_grid object is regular, the dimensions of the grid (size(lat_mat), $35 \times 43$), must be provided as well as the shape of the pixels and the spatial resolution of the grid, which is $0.05^{\circ} \times 0.05 ^{\circ}$. Second, the above step using the stem_krig_data constructor may appear redundant at first glance. Indeed, it is needed for compatibility with other model types for which, in addition to the stem_grid object, other information is also necessary for the stem_krig_data constructor.

Next, the stem_krig_options class provides some options for kriging.
By default, the output is back-transformed in the original unit of measure if the observations have been log-transformed and/or standardised. The back_transform property enables handling this.
Moreover, the no_varcov property must be set to 1 to avoid the time-consuming computation of the kriging variance.
Eventually, the block_size property is used to define the number of spatial locations in $\mathcal{S}%

o_krig_options = stem_krig_options();
o_krig_options.back_transform = 0;
o_krig_options.no_varcov = 0;
o_krig_options.block_size = 30;

After storing the map of Beijing boundaries into the o_model object, the latter is used with o_krig_data to create an object of class stem_krig.
This and o_krig_options together contain all information for kriging, which is obtained by the corresponding kriging method:

o_model.stem_data.shape = shaperead('../Maps/Beijing_adm1.shp');
o_krig = stem_krig(o_model, o_krig_data);
o_krig_result = o_krig.kriging(o_krig_options);

Note that this task may be time consuming for large grids. The kriging output saved in the o_krig_result object gives the latent process estimate $z_t$ and its variance.
The surface_plot and profile_plot methods may be used to obtain and plot $\hat{f}(\bm{s},h,t)$ of Equation (<ref>). In this case, the user has to provide the corresponding covariate (X_beta) for the scale/vector h, time t or location $\bm{s}$ (lon0, lat0) of interest.

Specifically, the surface_plot method is used to display the $O_{3}$ map using h, t, X_beta as input arguments. In the case of unavailable X_beta, the mapping concerns the component $\bm{\phi}(h)^{\top}\bm{z}(\bm{s},t)$. Loaded from the homonym file, the array X_beta_t_100 refers to time $t=100$ and hour $h=10.5$ and has the dimension $35 \times 43 \times 3$. Maps of $O_{3}$ concentrations and their standard deviation are shown in Figure <ref>.

load ../Data/kriging/X_beta_t_100;
t = 100;
h = 10.5;
[y_hat, diag_Var_y_hat] = o_krig_result.surface_plot(h, t, X_beta_t_100);

$O_3$ concentrations and their standard deviation at 10:30 am ($h = 10.5$), on 10 April 2015, where 12 stations are marked with black stars.

On the other hand, the profile_plot method is used to display the $O_{3}$ profile at a given spatial location $\bm{s}$ (lon0, lat0) and time t. Still, the profile plot concerns the component $\bm{\phi}(h)^{\top}\bm{z}(\bm{s},t)$ if X_beta is not provided. After loading the X_beta_h (dimension $25\times 3$) from the homonym file, this method represents the profile of $O_{3}$ concentrations and their variance–covariance matrix as in Figure <ref>:

load ../Data/kriging/X_beta_h;
h = 0:24;
lon0 = 116.25;
lat0 = 40.45;
t = 880;
[y_hat, diag_Var_y_hat] = o_krig_result.profile_plot(h, lon0, lat0, ...
t, X_beta_h);

Note that the prediction in Equation (<ref>) and the variance in Equation (<ref>) are stored in the output arguments y_hat, and diag_Var_y_hat, respectively.

$O3$ concentrations with $90\%, 95\%, 99\%$ confidence bands (different shadings), and their variance–covariance at latitude 40.45, longitude 116.25 on 29 May 2017.

§ CASE STUDY ON CLIMATE DATA

In order to show the complexity-reduction capabilities of D-STEM v2, a data set of temperature vertical profiles collected by the radiosondes of the Universal Radiosonde Observation Program (RAOB) is now considered.
The profiles are observed over the Earth's
sphere, and they are misaligned, that is, each profile differs in terms of the number
of observations and altitude above the ground of each observation.
Additionally, the computation burden is higher due to the higher number of
spatial locations at which profiles are observed.

§.§ RAOB data

Radiosondes are routinely launched from stations all over the world to
measure the state of the upper troposphere and lower stratosphere. Data collected
by radio sounding have applications in weather prediction and climate studies.

Temperature data from 200 globally distributed stations collected daily
during January 2015 at 00:00 and 12:00 UTC are considered here.
Each profile consists of a given number of measurements taken at different
pressure levels. Since the weather balloon carrying the radiosonde usually
explodes at an unpredictable altitude, the profile measurements
are misaligned across the profiles and have different pressure ranges. A
functional data approach is natural in this case since the underlying temperature
profile can be seen as a continuous function sampled at some pressure levels.

Figure <ref> depicts the spatial locations of temperature measurements
taken on 1 January 2015 at 00:00 UTC. This demo data set, which only covers one month,
includes around $10^5$ data points. When the full data set is used in climate studies,
the number of data points grows to around $10^8$. In this case, a recent server machine with multiple CPUs
with at least 256 GB of RAM is required for model estimation and kriging.

The focus of the case study is on the difference between the radiosonde
measurement and the output of the ERA-Interim global atmospheric reanalysis
model provided by ECMWF. In particular, the aim is to study the spatial
structure of the this difference in 4D space, where the
dimensions are latitude, longitude, altitude and time.

The model for temperature $y$ is as follows

\[
y \left(  \bm{s},h,t\right)  =x_{ERA}\left(  \bm{s},h,t\right)
\beta_{ERA}\left(  h\right)  +\bm{\phi}\left(  h\right)  ^{\top}\bm{z}\left(  \bm{s},t\right)  +\varepsilon\left(  \bm{s}%
\]

where $h\in\left[ 50, 925\right]  $ $hPa$ is the pressure level, while
$t=1,...,62$ is a discrete time index for January 2015. Figure <ref> shows the temperature measurements at a given station, where 50 and 925 $hPa$ correspond approximately to 25 and 1.3 km, respectively.

Temperature at location 5.25 latitude and -3.93 longitude in January 2015. Left: each dot is a temperature measurement. The colour of the dot depicts the temperature. Right: each graph is a temperature vertical profile collected through radio sounding.

§.§ Software implementation

This section details the software implementation of the case study described above as in script

demo_section5.m, which can be also executed in the demo_menu_user.m script. To avoid repetition, only the relevant parts of the script that differ from the case study of Section <ref> are reported and commented on here. In particular, data
loading and instantiation of the stem_model object are not described.

§.§.§ Model estimation

The problem of vertical misalignment of the measurements is completely transparent to the user and is handled by the internal stem_misc.data_formatter method when creating the stem_data object.
Note that the dimension of the matrices in o_varset depends on $q$, the maximum number of measurements in each profile.
To prevent out-of-memory problems, it is advisable to avoid data sets in which only a few profiles have a large number of
measurements, which could result in large matrices in
o_varset, with most of the elements set to NaN.

B-spline bases are used, since, in this application, vertical profiles are not periodic with respect to the pressure domain. The corresponding object of class stem_fda is created in the following way:

spline_order = 2;
rng_spline = [50,925];
knots_number = 5;
knots = linspace(rng_spline(1),rng_spline(2),knots_number);

input_fda.spline_type = 'Bspline';

input_fda.spline_order = spline_order;
input_fda.spline_knots = knots;

input_fda.spline_range = rng_spline;
o_fda = stem_fda(input_fda);

Note that the knots are equally spaced along the functional range. In general,
however, non-equally spaced knots can be provided, and each model component
(i.e. $\sigma_{\varepsilon}^{2}$, $\beta_{j}$ and $\bm{\phi}(h)^{\top}\bm{z}(\bm{s}%
,t)$) can have a different set of knots. This is obtained using spline_order and spline_knots with additional suffixes _sigma, _beta, _z.

Although this data set is not large, the demo shows how to
enable the partitioning discussed in Section <ref>. First, the spatial locations are partitioned using the modified k-means algorithm:

k = 5;
trials = 100;
lambda = 5000;
partitions = o_data.kmeans_partitioning(k, trials, lambda);

where k is the number of partitions, trials is the number of times when the k-means
algorithm is executed starting from randomised centroids and lambda is $\lambda$ in Equation (<ref>).

At the end of the k-means algorithm, data are internally re-ordered for parallel
computing. Model estimation is done after creating and setting an object of class
stem_EM_options. To do this, the output of the kmeans_globe method is passed to the
partitioning_block_size property of the o_EM_options object.
Additionally, for parallel computing, the number of workers must be set to a
value higher than 1. In general, this could be any number up to the number of
cores available on the machine.

o_EM_options = stem_EM_options();
o_EM_options.partitions = partitions;
o_EM_options.workers = 2;

The three validation MSEs defined in Section <ref> are shown in Figure <ref> and <ref>. To generate these figures, the method plot_validation is called with vertical = 1, which provides “atmospheric profile” plots with $h$ on the vertical axis:

vertical = 1;

(Left) Validation MSE with respect to the $\bar{h}$ coloured by the number of observations $n_b$, and (Right) the validation MSE with respect to time $t$.

Validation MSE for the thirty-three stations, where the stations used for estimation are marked with blue stars.

§.§.§ Kriging

Interpolation across space and over time is done as in Section
<ref>. However, complexity reduction is enabled by adopting the
nearest neighbour approach detailed in Section <ref>.

To do this, a class constructor is first called, where the block_size is used to define the number of spatial locations in $\mathcal{S}_{l}^{\ast}$, and then nn_size is used to define $\tilde{n}$. Additionally, setting
o_krig_options.workers makes it possible to do the kriging over the $u$ blocks in
parallel using up to the allocated number of workers:

o_krig_options = stem_krig_options();
o_krig_options.block_size = 150;
o_krig_options.nn_size = 10;
o_krig_options.workers = 2;

Finally, kriging predictions and standard errors are mapped for a given $h\in\mathcal{H}$ and time $t$:

h = 875.3;
t = 12;
o_krig_result.surface_plot(h, t);

Since covariates are not provided to the surface_plot method, the plots are on the component
$\bm{\phi}(h)^{\top}\bm{\hat{z}}(\bm{s},t)$, namely, the difference between RAOB and ERA-Interim and its standard deviation. The output
of the above code is depicted in Figures <ref> and <ref>.

$\bm{\phi}(h)^{\top}\bm{\hat{z}}(\bm{s},t)$ at pressure $875.3$ $hPa$, and 12:00 am on 6 January 2015, where $200$ stations are shown as black stars.

Standard deviation of $\bm{\phi}(h)^{\top}\bm{\hat{z}}(\bm{s},t)$ at pressure $875.3$ $hPa$, and 12:00 am on 06 January 2015, where $200$ stations are shown as black stars.

§ CONCLUDING REMARKS

This paper introduced the package D-STEM v2 through two case studies of spatio-temporal modelling of functional data. It is shown that, in addition to maximum likelihood estimation, Hessian approximation and kriging for large
data sets, D-STEM v2 also develops several data-handling capabilities, allows for automatic construction of relevant objects and provides graphical output. In particular, it provides high-quality global maps and two kinds of functional
plotting: the traditional x–y plot and the vertical profile plot, which is popular, for example, in atmospheric data analysis.
In this regard, model validation and kriging are straightforward.

D-STEM v2 fills a gap in functional geostatistics. In fact, although statistical methods for georeferenced functional data have been recently
developed (e.g. [Ignaccolo et al., 2014]), standard geostatistical packages do not consider functional data, especially in the spatio-temporal context.

The successful use of D-STEM v1 in a number of applications proved that the EM algorithm implementation is quite stable. Now, due to improvements in computational efficiency, the new D-STEM v2 has the capability to handle large data sets.
Moreover, thanks to the approximated variance–covariance matrix, it is possible to compute standard errors for all model parameters relatively fast and avoid the large number of iterations typically required by an MCMC approach for making inferences.

However, a limit of the EM algorithm is its limited flexibility to changes in the model equations. Indeed, changes in parametrisation or latent variable structure usually require deriving new closed-form estimation formulas and changing the software accordingly. Moreover, changes in covariance functions are not easy to handle.

Computationally, the main limit of D-STEM v2 is in the number $p$ of basis functions that can be handled.
Even if partitioning is exploited in $k$ blocks of size $r$, computational complexity is $\mathcal{O}\left(Tkr^{3}p^{3}\right)$, meaning that $p$ cannot be large.

Currently, the authors are working on a new version, which makes it possible to handle multivariate functional space-time data and user-defined spatial covariance functions, which will make D-STEM v2 a valid and comprehensive alternative to the Gaussian process regression models (fitrgp) implemented in the Statistics and Machine Learning Toolbox of MATLAB.

§ ACKNOWLEDGMENTS

We would like to thank China's National Key Research Special Program Grant (No. 2016YFC0207701) and the Center for Statistical Science at Peking University.

[Bivand et al., 2015]
Bivand RS, Gómez-Rubio V, Rue H (2015).
Spatial Data Analysis with R-INLA with Some
Journal of Statistical Software, 63(20), 1–31.

[Blangiardo et al., 2013]
Blangiardo M, Cameletti M, Baio G, Rue H (2013).
Spatial and Spatio-Temporal Models with
Spatial and Spatio-Temporal Epidemiology, 4, 33–49.

[Calculli et al., 2015]
Calculli C, Fassò A, Finazzi F, Pollice A, Turnone A (2015).
Maximum Likelihood Estimation of the Multivariate Hidden
Dynamic Geostatistical Model with Application to Air Quality in Apulia,
Environmetrics, 26(6), 406–417.

[Cressie, 2018]
Cressie N (2018).
Mission CO2ntrol: A Statistical Scientist's Role in Remote
Sensing of Atmospheric Carbon Dioxide.
Journal of the American Statistical Association,
113(521), 152–168.

[Cressie and Johannesson, 2008]
Cressie N, Johannesson G (2008).
Fixed Rank Kriging for Very Large Spatial Data Sets.
Journal of the Royal Statistical Society B, 70(1),

[Dohan and Masschelein, 1987]
Dohan JM, Masschelein WJ (1987).
The Photochemical Generation of Ozone: Present
Ozone: Science $\&$ Engineering, 9(4), 315–334.

[Fassò, 2013]
Fassò A (2013).
Statistical Assessment of Air Quality Interventions.
Stochastic Environmental Research and Risk Assessment,
27(7), 1651–1660.
<DOI: 10.1007/s00477-013-0702-5>.

[Fassò and Finazzi, 2011]
Fassò A, Finazzi F (2011).
Maximum Likelihood Estimation of the Dynamic
Coregionalization Model with Heterotopic Data.
Environmetrics, 22(6), 735–748.

[Fassò et al., 2016]
Fassò A, Finazzi F, Ndongo F (2016).
European Population Exposure to Airborne Pollutants Based on
a Multivariate Spatio-Temporal Model.
Journal of agricultural, biological, and environmental
statistics, 21(3), 492–511.

[Fassò et al., 2014]
Fassò A, Ignaccolo R, Madonna F, Demoz B, Franco-Villoria M (2014).
Statistical Modelling of Collocation Uncertainty in
Atmospheric Thermodynamic Profiles.
Atmospheric Measurement Techniques, 7(6), 1803–1816.

[Finazzi, 2020]
Finazzi F (2020).
Fulfilling the Information Need after an Earthquake:
Statistical Modelling of Citizen Science Seismic Reports for Predicting
Earthquake Parameters in Near Realtime.
Journal of the Royal Statistical Society: Series A,
183(3), 857–882.

[Finazzi and Fassò, 2014]
Finazzi F, Fassò A (2014).
D-STEM: A Software for the Analysis and Mapping of
Environmental Space-Time Variables.
Journal of Statistical Software, 62(6), 1–29.

[Finazzi et al., 2018]
Finazzi F, Fassò A, Madonna F, Negri I, Sun B, Rosoldi M (2018).
Statistical Harmonization and Uncertainty Assessment in the
Comparison of Satellite and Radiosonde Climate Variables.
arXiv preprint arXiv:1803.05835.

[Finazzi et al., 2015]
Finazzi F, Haggarty R, Miller C, Scott M, Fassò A (2015).
A Comparison of Clustering Approaches for the Study of the
Temporal Coherence of Multiple Time Series.
Stochastic Environmental Research and Risk Assessment,
29(2), 463–475.

[Finazzi et al., 2019]
Finazzi F, Napier Y, Scott M, Hills A, Cameletti M (2019).
A Statistical Emulator for Multivariate Model Outputs with
Missing Values.
Atmospheric Environment, 199, 415 – 422.

[Finazzi et al., 2013]
Finazzi F, Scott EM, Fassò A (2013).
A Model-Based Framework for Air Quality Indices and
Population Risk Evaluation, with an Application to the Analysis of Scottish
Air Quality Data.
Journal of the Royal Statistical Society: Series C,
62(2), 287–308.

[Finley et al., 2015]
Finley A, Banerjee S, Gelfand A (2015).
spBayes for Large Univariate and Multivariate
Point-Referenced Spatio-Temporal Data Models.
Journal of Statistical Software, Articles, 63(13),
ISSN 1548-7660.

[Furrer et al., 2006]
Furrer R, Genton MG, Nychka D (2006).
Covariance Tapering for Interpolation of Large Spatial
Journal of Computational and Graphical Statistics,
15(3), 502–523.

[Gasch et al., 2015]
Gasch CK, Hengl T, Gräler B, Meyer H, Magney TS, Brown DJ (2015).
Spatio-Temporal Interpolation of Soil Water, Temperature,
and Electrical Conductivity in 3D+T: The Cook Agronomy Farm Data Set.
Spatial Statistics, 14, 70–90.

[Gramacy, 2016]
Gramacy RB (2016).
laGP: Large-Scale Spatial Modeling via Local
Approximate Gaussian Processes in R.
Journal of Statistical Software, 72(1), 1–46.

[Heaton et al., 2018]
Heaton MJ, Datta A, Finley AO, Furrer R, Guinness J, Guhaniyogi R, Gerber F,
Gramacy RB, Hammerling D, Katzfuss M, Lindgren F, Nychka DW, Sun F,
Zammit-Mangion A (2018).
A Case Study Competition among Methods for Analyzing Large
Spatial Data.
Journal of Agricultural, Biological and Environmental
Statistics, pp. 1–28.

[Ignaccolo et al., 2014]
Ignaccolo R, Mateu J, Giraldo R (2014).
Kriging with External Drift for Functional Data for Air
Quality Monitoring.
Stochastic Environmental Research and Risk Assessment,
28(5), 1171–1186.

[Jurek and Katzfuss, 2018]
Jurek M, Katzfuss M (2018).
Multi-Resolution Filters for Massive Spatio-Temporal Data.
arXiv preprint arXiv:1810.04200.

[Kahle and Wickham, 2013]
Kahle D, Wickham H (2013).
ggmap: Spatial visualization with ggplot2.
The R Journal, 5(1), 144–161.

[Katzfuss, 2017]
Katzfuss M (2017).
A Multi-Resolution Approximation for Massive Spatial
Journal of the American Statistical Association,
112(517), 201–214.

[Lindgren and Rue, 2015]
Lindgren F, Rue H (2015).
Bayesian Spatial Modelling with R-INLA.
Journal of Statistical Software, 63(19), 1–25.

[Negri et al., 2018]
Negri I, Fassò A, Mona L, Papagiannopoulos N, Madonna F (2018).
Modeling Spatiotemporal Mismatch for Aerosol Profiles.
In Quantitative Methods in Environmental and Climate Research,
pp. 63–83. Springer-Verlag.

[Nychka et al., 2015]
Nychka D, Bandyopadhyay S, Hammerling D, Lindgren F, Sain S (2015).
A Multiresolution Gaussian Process Model for the Analysis of
Large Spatial Datasets.
Journal of Computational and Graphical Statistics,
24(2), 579–599.

[Nychka et al., 2016]
Nychka D, Hammerling D, Sain S, Lenssen N (2016).
LatticeKrig: Multiresolution Kriging Based on Markov
Random Fields.
University Corporation for Atmospheric Research, Boulder, CO, USA.
R package version 7.0.

[Pebesma, 2012]
Pebesma E (2012).
spacetime: Spatio-temporal Data in R.
Journal of statistical software, 51(7), 1–30.

[Pebesma and Heuvelink, 2016]
Pebesma E, Heuvelink G (2016).
Spatio-temporal Interpolation Using gstat.
RFID Journal, 8(1), 204–218.

[Porcu et al., 2018]
Porcu E, Alegria A, Furrer R (2018).
Modeling Temporally Evolving and Spatially Globally
Dependent Data.
International Statistical Review, 86(2), 344–377.

[Ramsay and Silverman, 2007]
Ramsay JO, Silverman BW (2007).
Applied Functional Data Analysis: Methods and Case Studies.

[Ramsay et al., 2018]
Ramsay JO, Wickham H, Graves S, Hooker G (2018).
fda: Functional Data Analysis.
R package version 2.4.8,

[Rue et al., 2014]
Rue H, Martino S, Blangiardo FL, Simpson D, Riebler A, Krainski ET (2014).
INLA: Functions which Allow to Perform Full Bayesian
Analysis of Latent Gaussian Models using Integrated Nested Laplace
R package version 0.0-1404466487,

[Shumway and Stoffer, 2017]
Shumway RH, Stoffer DS (2017).
Time Series Analysis and Its Applications: with R

[Stein, 2002]
Stein ML (2002).
The Screening Effect in Kriging.
The Annals of Statistics, 30(1), 298–323.

[Stein, 2013]
Stein ML (2013).
Statistical Properties of Covariance Tapers.
Journal of Computational and Graphical Statistics,
22(4), 866–885.

[Taghavi-Shahri et al., 2019]
Taghavi-Shahri S, Fassò A, Mahaki B, Amin H (2019).
Concurrent Spatiotemporal Daily Land Use Regression Modeling
and Missing Data Imputation of Fine Particulate Matter Using Distributed
Space-Time Expectation Maximization.
Atmospheric Environment, 224, 1–11.

[Tzeng and Huang, 2018]
Tzeng S, Huang HC (2018).
Resolution Adaptive Fixed Rank Kriging.
Technometrics, 60(2), 198–208.

[Wan et al., 2020]
Wan Y, Xu M, Huang H, Chen S (2020).
A Spatio-Temporal Model for the Analysis and Prediction of
Fine Particulate Matter Concentration in Beijing.
Environmetrics, accepted.

[Wang et al., 2017]
Wang T, Xue L, Brimblecombe P, Lam YF, Li L, Zhang L (2017).
Ozone Pollution in China: A Review of Concentrations,
Meteorological Influences, Chemical Precursors, and Effects.
Science of the Total Environment, 575, 1582–1596.

[Zammit-Mangion, 2018]
Zammit-Mangion A (2018).
FRK: Fixed Rank Kriging.
R package version 0.2.2,

get arXiv to do 4 passes: Label(s) may have changed. Rerun