Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
57,300 characters · 11 sections · 67 citation commands
Forecasting Environmental Data: An example to ground-level ozone concentration surfaces
Studies on the prediction of environmental data sampled across time and over geographical areas/regions have received wide attention across many fields. Examples of such studies include air quality control (see, e.g., aue2015, Mu2018), water quality control (see, e.g., Canedo2016 and Edurne2018), temperature (see, e.g., Kuenzer2021) and many more. Although the analysis and prediction of such environmental data are available in the literature, new conceptual approaches are required because of the increasing complexity of the collected data. The first challenge that appears on the way of constructing an appropriate forecasting methodology lies in the fact that data collection is often scattered over some geographical region with complicated shapes and irregular boundaries. The second challenge is that such data usually interact in spatial and temporal aspects and both have to be taken into account to provide reliable forecasts. In this paper, we adopt recent advances in functional data analysis (FDA) to address these challenges and construct a new forecasting methodology.
Prediction of environmental data with FDA methods is not new in the literature (see e.g, RamsaySilverman2005 and KokoszkaReimherr2017 for an overview of FDA). However, most of the applications study only one of the key aspects of the data, either the spatial one or the temporal one. Current advances in FDA allow us to reconstruct two dimensional manifold (or surface) from the environmental data sets sampled across complicated domains for one time period. For a review of the available methods see, for instance, Ramsay2002, SangalliRamsayRamsay2013, Ettinger2016, Mu2018 and Ferraccioli2020. In the last decade, we have also seen a rise of literature addressing modelling and forecasting of functional time series (see e.g., Bosq2000, HYNDMAN2009, HoermannKokoszka2010, HORMANN2012, and salish2019). The majority of the existing contributions in the forecasting of functional time series have been tailored to one-dimensional manifolds rather than surfaces. The closest literature in FDA that takes into account both spatial and temporal aspects of data is related to the analysis of spatial functional processes. (see e.g., Delicadoetal2010). Within this literature, it is assumed that one observes a family of $N$ functions $( {\cal Y}_{\boldsymbol{s}_i})_{i=1}^N$ at locations $\boldsymbol{s}_1,\ldots,\boldsymbol{s}_N$, which are thought to exhibit spatial dependencies. The main interest is to exploit this dependence structure when predicting a function ${\cal Y}_{\boldsymbol{s}^*}(\cdot)$ at some unmonitored site $\boldsymbol{s}^*$. Several approaches have been suggested that are mainly based on the framework of functional linear regression models. For example, GiraldoDelicadoMateu2009, GiraldoDelicadoMateu2011 and NeriniMonestiezMante2010 extend the methodology of kriging from classical geostatistics to the functional setting while YamanishiTanaka2003 consider additional functional covariates.
In this paper, we propose to take a new approach to process and forecast spatial data collected sequentially over time. This approach is based on a synthesis of the tools available for spatial and temporal analysis of functional data. We treat an environmental data set, observed at a set of locations of a domain ${\cal D}$ over time period $T$, as a realization of a time series of surfaces defined as
For instance, the data that motivated this research consists of daily measurements of ground-level ozone concentration. For the case of Germany, daily measurements of ozone concentration are available at $N=171$ stations for 365 days (see the left panel of Figure (ref)).
It is expected that the ozone concentration is distributed continuously over the geographical territory of Germany, whereas information is available at the $N$ sampling locations at a regular time interval. Our main target is to obtain forecasts of a future surface $X_{T+1}(\boldsymbol{s})$, $\boldsymbol{s} \in {\cal D}$, given an available surface time series $\{X_t(\boldsymbol{s})\}_{t=1}^T$. In the case of ground-level ozone surfaces, such forecasts are especially useful for air quality control and identification of (short-term or long term) future areas that are at the high risk of elevated level of pollution. Note that the prediction problem of (ref) is conceptually different from the prediction of a spatial functional process. While we are interested in time forecast (reflecting the classical time series problem), the prediction of a spatial functional process is spatial in nature.
Given the typical structure of an environmental data set, we first have to reconstruct a corresponding surface time series from available discrete measurements. Further, as we are working with spatial data, we are interested in taking the geographical boundaries into account. Therefore, in the first step of our approach, we employ a finite element smoothing approach suggested in Ramsay2002 and SangalliRamsayRamsay2013. Once a surface time series is reconstructed from discrete observations we proceed with the forecasting step. To this end, we adapt the dynamic functional factor model (DFFM) to the settings of the surface time series. The DFFM is proven to be useful in many applications to one-dimensional manifold times series (See e.g., hays2012, Liebl2013, aue2015 and OttoSalish2021). The advantage of this forecasting approach over existing methods is the simplicity of its application and the flexibility of the modelling framework. The DFFM model allows decomposing rather complicated surface time series into a finite number of scalar factors that have predictable power and a reminder term (see Section (ref) for a more detailed description). As discussed and shown in aue2015 and OttoSalish2021 the simplicity of the application lies in the fact that estimates of the factors are obtained via standard functional principal component analysis. (For instance, R and Matlab packages for FDA by RamsaySilverman2005 can be used for this purpose). The flexibility of the model comes from the fact that, after the factors are estimated, any suitable multivariate forecasting technique can be used to obtain a forecast of the factors and therefore of the surface $X_{T+1}(\cdot)$. As shown by the identification results in OttoSalish2021 the behaviour of the factors is restricted only by weak stationarity (with some moment restrictions) indicating that nonlinear multivariate forecasting techniques can potentially be employed. Therefore, in this paper we propose two specifications of the adapted DFFM model to surface time series: one that accounts for a linear behaviour of the factors (using the vector autoregressive model) and the other for nonlinear (using K-nearest neighbours framework).
We demonstrate the practical value of this new modelling and forecasting perspective with an application to ground-level (or tropospheric) ozone, which is a harmful air pollutant, over the geographical domain of Germany. The World Health Organization (WHO) and the European Environment Agency\footnote{\url{https://www.eea.europa.eu/themes/air/health-impacts-of-air-pollution}} report that exposure to high concentrations of this pollutant can cause breathing problems and cardiovascular diseases. It is estimated that ground-level ozone in 2019 was linked to 16800 premature deaths in the EU.\footnote{\url{https://www.eea.europa.eu/publications/air-quality-in-europe-2021/health-impacts-of-air-pollution}} Hence, it is one of the pollutants which is closely monitored by environmental agencies. This makes the prediction of the concentration level in different geographical locations especially useful as a preventive tool. In particular, predictions (in space and time) can help to detect and anticipate local areas where air quality thresholds might be exceeded and additional measures could/should be implemented to meet the standards set by European Union. Our framework can predict well episodes where the ozone layer exceeding the threshold of WHO standards. In addition to the forecasts, it is capable of providing an estimation of the areas that contribute most to the variability of the ozone surface in Germany (via estimated loading functions). On the final note, we also supplemented our study with standard benchmark models commonly employed in functional time series analysis to provide a full comparative forecasting analysis. The DFFM model accounting for the possibly nonlinear behaviour of factors showed the best forecasting performance in this study.
The remainder of the paper is organized as follows. Section (ref) describes how smooth surfaces can be obtained from noisy discrete spatial measurements by using a finite element spline smoother. Section (ref) discusses the dynamic functional factor model and details the estimation and forecasting steps. A forecasting study of ground-level ozone concentration surfaces over the geographical domain of Germany is reported in section (ref). Section (ref) concludes.
Our main object of interest is a time series of surfaces, $\{X_t({s})\}_{t=1}^T$, over some spatial domain ${s}\in{\cal D}$. However, in practice we only observe surface's value, along with some noise, at discrete locations $\{{s}_i\}_{i=1}^N$ on the domain ${\cal D}$, where each spatial location $\boldsymbol{s}_i$ is represented by pairs $(x_i,y_i)$ of longitude and latitude coordinates. As a consequence, some statistical smoothing procedure is required to approximate $X_t$ for each $t$. In this paper we resort to one of the most established methods in the literature, referred to as finite element spline smoother and discuss its implementation. This method is in particular attractive as it allows to account for complicated geometry of the domain of interest. Further details on this approach can be found in Ramsay2002, SangalliRamsayRamsay2013 and Ferraccioli2020.
The spline estimate $\widetilde{X}_{t}$ of $X_t$ is defined as a minimizer of the following functional
where data points are given as $\{{s}_i,X^*_t({s}_i)\}_{i=1}^N$, with $X^*_t$ being a noisy observation of $X_t$. Functional $J_{\lambda}({\cal X})$ is a penalized sum of squared errors, where the first term measures the fit to the data and the second term governs the roughness (or smoothness) of the obtained approximation. It is sufficient for the purpose of this paper to measure roughness through a (suitably defined) notion of squared second derivatives. We thus restrict our attention to $L^2$-functions that have square-integrable derivatives up to second order and denote the corresponding functional space by $H^2({\cal D})$.
We require that smoothing problem (ref) is independent of the underlying spatial coordinate system so that the roughness penalty should be invariant under translation and rotation of the spatial coordinates. This is ensured if a differential operator $\Delta$ in the second term is comprised of polynomials of the Laplacian operator (see e.g., Folland1995 or Ramsay2002), where for any ${\cal X} \in H^2({\cal D})$
Furthermore, parameter $\lambda>0$ is a smoothing parameter, for which large values are used to provide smoother estimates and small values are used to estimate with a better fit.
As shown by SangalliRamsayRamsay2013 a unique solution to the minimization problem in (ref) exists if one imposes additional boundary conditions. That is, the minimization is solved over ${\cal X} \in H^2_{n0}({\cal D})$, where $H^2_{n0}({\cal D})$ denotes the space of $L^2$-functions that have square-integrable partial derivatives up to second order and assume zero normal derivatives at the boundary. The solution in practice is approximated using finite element spline smoother, which proceeds in three steps:
\noindentImplementation Details:
Step 1 is typically done by triangulation of the domain ${\cal D}$, where the spatial locations $\boldsymbol{s}_1, \boldsymbol{s}_2, \ldots, \boldsymbol{s}_N$ correspond to vertices of the resulting triangular mesh. How to choose the set of triangles in practice is difficult and many possibilities have been offered in the literature. For our purposes, the Delaunay triangular mesh seems the most appropriate as it avoids thin triangles and favours triangles that are as equiangular as possible. As long as no four or more nodes lie on a common circle, the Delaunay triangulation is uniquely defined and implementation routines are readily available for most statistical software. Figure (ref) shows an example of the Delaunay triangular mesh for ozone measurement stations in Germany. The triangulation of the domain ${\cal D}$ in the remainder of the paper is denoted as $\triangle_{{\cal D}}$.
We wish to approximate a surface $X_t$ over $\triangle_{{\cal D}}$ (for each $t$) through the polynomials of second order over any triangle while being continuous over edges and vertices. For this purpose in Step 2, we construct triangular finite elements, each consisting of a single triangle from $\triangle_{{\cal D}}$, a set of six nodes and associated nodal basis functions. To construct a quadratic polynomial over a triangle, the function value is specified at six nodal points, which are the vertices and the midpoints of each edge of a triangle as indicated in Figure (ref). With each of the nodal points, we associate a shape function which is a second order polynomial that takes the value one at one local nodal point and the value zero at all other local nodal points. The six shape functions that are constructed in such a way are plotted in Figure (ref).
Let us denote the union of nodal points of the triangulation by $\boldsymbol{\xi}_k$, $k = 1,\ldots,K$. For ease of notation, we number the nodal points in such a way to have the spatial locations $\boldsymbol{s}_i$, $i=1,\ldots,N$, correspond to the first $N$ nodal points. We associate with each node $\boldsymbol{\xi}_k$, $k = 1,\ldots,K$, a nodal basis function $\phi_k$. Each nodal basis function is constructed piecewise over each triangle in $\triangle_{{\cal D}}$. More precisely, $\phi_k$ takes the value of a shape function if this shape function has value 1 at the node $\boldsymbol{\xi}_k$ and $\phi_k$ is zero otherwise. In Figure (ref) we plot the resulting finite element nodal basis function associated with the nodal point $(0,0)$ which is shared by the four triangles. This set of $K$ basis functions spans a function space that we denote by $H^1(\triangle_{{\cal D}})$, i.e. the space of continuous functions on ${\cal D}$ that are piecewise quadratic polynomials when restricted to some triangular finite element.
In Step 3 we obtain a computable representation of the finite element spline approximation to the solution $\widetilde{X}$. Denote by $\boldsymbol{\phi}_K := \left( \phi_1, \phi_2, \ldots, \phi_K \right)^{\mathsf{T}}$ the $K$-vector of spatial basis functions $\phi_k$. Let us furthermore denote the $K$-vectors of partial derivatives of the nodal basis functions with respect to spatial $x$ and $y$ coordinates as $\boldsymbol{\phi}_K^{(x)} := \left( \partial \phi_1 / \partial x, \ldots, \partial \phi_K / \partial x \right)^{\mathsf{T}}$ and $\boldsymbol{\phi}_K^{(y)} := \left( \partial \phi_1 / \partial y, \ldots, \partial \phi_K / \partial y \right)^{\mathsf{T}}$ and define $K\times K$ matrices
and $\boldsymbol{D}_{K}$ is a $K\times K$ diagonal matrix that has $i$-th diagonal element 1 if the $i$-th node is a data point and 0 otherwise. By construction of the finite elements in Step 2 we can write an approximation of the solution $ \widetilde{X}_t$ to minimization problem (ref) as
Hence, the minimization problem is equivalent to searching for a $K\times 1$ vector $\boldsymbol{\widetilde{X}}_K$. As shown in SangalliRamsayRamsay2013 solving for $\boldsymbol{\widetilde{X}}_K$ is then equivalent to solving the system of linear equations given by
where $[\boldsymbol{ \widetilde{X}}_K, \boldsymbol{Z}_K ]^{\mathsf{T}}$ is a vector of unknowns and $\boldsymbol{X}_K^{*}:= \left( X^*(\boldsymbol{s}_1), X^*(\boldsymbol{s}_2), \ldots, \linebreak[0] X^*(\boldsymbol{s}_N),0, \ldots, 0 \right)^{\mathsf{T}}$ is the $K\times 1$ vector which has on the first $N$ entries the observations $\boldsymbol{X}_t^{*}((\boldsymbol{s}_i)$ for $i=$ and zeros otherwise (note that the nodal points $\boldsymbol{\xi}_k$ are numbered such that the first $N$ nodes correspond to the spatial locations of measurement stations).
As discussed in the introduction, once the surface time series is reconstructed, we are interested in predicting its future values $X_{T+h}$ for $h>0$, given observed values at $t=1,...,T$. For this purpose, we adopt the dynamic functional factor model (DFFM, hereafter) to the settings where the data points are given as surfaces over some domain. We also review at the end of this section other available techniques in FDA literature and standard benchmarks.
The main attractive feature of the DFFM is that it allows decomposing the original surface time series into a low-dimensional predictable component and an idiosyncratic component with no predictive power. The model is written as
where a $L\times 1$ vector of factors $\mathbf{x}_{t}=[x_{1,t},...,x_{L,t}]'$ describe the dynamic nature of the original process $X_{t}$, whereas an intercept function, $\mu(\cdot)$, and the vector of loading surfaces, $\Psi(\cdot)=[\psi_{1}(\cdot),...,\psi_{L}(\cdot)]'$, are unobserved deterministic terms. The product $\Psi'\mathbf{x}_{t}$ is referred to in the literature as the common component and $\{\epsilon_{t}\}$ denote a series of idiosyncratic error terms that do not have any predictive power. We consider model (ref) under the fairly general assumption specified below.
Recall $L^2({\cal D})$ denote the space of square-integrable functions defined over the spatial domain ${\cal D}$ and equipped with the inner product $\langle x,y\rangle = \int\limits_{{\cal D}} x(s) y(s) \mathrm{d}s$ and the norm $\Vert x\Vert=\langle x,x\rangle^{1/2}$.
This set of assumptions allows us to identify latent components of the model as functional principal components extracted from the global covariance function of $\{ X_{t}\}$ defined as $c(s,r):=\lim\limits_{T\to\infty}\frac{1}{T}Cov\left[X_t(s)X_t(r)\right]$ for $s,r\in{\cal D}$. (See Theorem 1 in OttoSalish2021 for the theoretical discussion of this result). This, in turn, indicates that all latent components of model (ref) can be easily estimated. Section (ref) guides the practical implementation of this step. Further, Assumptions (ref) and (ref) do not impose a particular modelling structure on the factors allowing for some flexibility when choosing a forecasting approach. In Section (ref) we consider two possibilities: (i) factors $\mathbf{x}_{t}$ are forecasted with linear vector autoregressive model and (ii) with k-nearest neighbours to account for potential nonlinearities in the dynamics of $\mathbf{x}_{t}$.
The estimation of all latent components in model (ref), such as $\mu$, $\{\psi\}_{l=1}^L$ and $\mathbf{x}_{t}$, is done via the first two sample moments of the original process $X_t$. More precisely, let the first sample moments of $X_t$ be denoted as
and the sample covariance function as
Covariance function $\widehat{c}(r,s)$ is a kernel of the sample covariance integral operator $\widehat C_Y$, which in turn has the eigenvalues $\widehat \lambda_1 \geq \widehat \lambda_2 \geq \ldots, \geq \widehat \lambda_T \geq 0$ and corresponding orthonormal eigenfunctions $\widehat \psi_1, \ldots, \widehat \psi_T$. These functions are used as estimators of loading functions $\psi_1, \ldots, \psi_L$. Finally, the projections $\widehat x_{l,t} = \langle X_t - \widehat \mu, \widehat \psi_l \rangle$ are used as an estimator of factor $x_{l,t}$. As shown in OttoSalish2021, Theorem 2 these estimators provide consistent estimates. For practical implementation of this step, we refer a reader to “{fda}” package of RamsaySilverman2005.
The best $h$-step ahead predictor (in the mean square error sense) of the process of interest $X_t$ is given as,
where $\mathbf{x}_{T+h|T}$ denotes an $h$-step ahead predictor of the factors, which in turn is given as
with $\mathbf{x}_{T+j|T} = \mathbf{x}_{T+j}$ if $j \leq 0$. Since estimates of $\mu(s)$ and $\Psi(s)$ are already discussed in Section (ref) to obtain the final predictor of the surface $X_{T+h|T}$ it remains to estimate $ \mathbf{x}_{T+h|T}$. For the simplicity of the discussion and without loss of generality we will proceed to discuss a case with one-step ahead forecasts.
Linear case. If dynamic elements of the common component are modelled linearly then we can resort to the classical multivariate time series literature, which gives us plenty of methodological approaches to obtain forecast $\mathbf{x}_{T+1|T}$. (see e.g, luetkepohl2005 for the review of multivariate time series literature). For instance, it is common to approximate the dynamics of $\mathbf{x}_{t}$ with the Vector Autoregressive Model (VAR). This will result in the following simplified expression of $\mathbf{x}_{T+1|T}$
where $\mathbf{A}_1,...,\mathbf{A}_p$ denote the $L\times L$ coefficient lag matrices and $p$ is the number of relevant lags. Predictor (ref) can be easily obtained by estimating a VAR model of order $p$ using any statistical software capable of processing time series data. The estimation of the latent factors $\{\mathbf{x}_{t}\}_{t=1}^T$, as discussed in Section (ref), is obtained via principal components of $\widehat{c}(s,r)$. Hence, it only remains to specify the relevant number of factors, $L$, and the number of lags, $p$, necessary to estimate the forecast in (ref).
Estimators of $L$ and $p$ are obtained by employing the information criterion from aue2015 (ANH) or one of the variants of the information criterion proposed in OttoSalish2021 (OS). The main difference between these criteria is similar to the difference between AIC and BIC (or HQ) in the classical time series analysis (See e.g., luetkepohl2005). That is, while the first one (ANH) is tailored to minimize the sample mean squared error of the forecasts, the second one selects $L$ and $p$ consistently. Hence, both criteria are of interest to our application.
To facilitate the practical implementation of this forecasting routine we provide below a summary of the main steps:
Non-Linear case. When there is little a priori knowledge about the shape of (ref) or there is strong evidence about a non-linear relationship between ${\mathbf{x}}_{t}$ and its past a different forecasting technique to the one in (ref) might be appropriate/required. Because of its flexibility, a nonparametric k-nearest neighbours (KNN) method is often used when linearity assumption is not imposed. The KNN has been studied extensively in the statistical literature (see e.g., Gyorfi2002 for a general discussion). Although the theoretical treatment of KNN estimates is found difficult the efficiency and simplicity of its application made it prominent in the applied literature. Hence, in this paper, we also consider an estimator of forecast ${\mathbf{x}}_{T+1|T}$ using the KNN principle. That is, we replace (ref) with
and our goal is to provide an estimation of $m(\cdot)$ locally.
The DFFM model allows us to model dynamics of the original process $\{X_t\}$ via $L$-dimensional multivariate process ${\mathbf{x}}_{t}$ . Hence, instead of developing a new KNN method for functional (time dependent) data (as for instance in Biau2010, KUDRASZOW2013 and KARA2017) we can adopt techniques already available for classical time series (see e.g., Yakowitz1987). The estimator of $m(\cdot)$ and therefore the estimator of the forecast $\mathbf{x}_{T+1|T}$ is obtained following the next basic steps.
The key method in functional data analysis used for forecasting is the functional autoregressive model (FAR). Since the seminal monograph by Bosq2000 this model has seen a lot of development and have become the main toolbox for the analysis of functional time series. (See, e.g., besse2000, Mas2007, kargin2008, KokoszkaReimherr2013 etc.). Hence, it is natural to use it as the main benchmark model for forecasting functional/surface time series. We follow closely the exposition in HorvathKokoszka2012 to describe the model, its estimation and forecasting procedure. For brevity, we assume that process $\{X_t(s)\}$ for $s\in{\cal D}$ is zero mean such that the FAR model is given as
Model (ref) satisfies standard regularity conditions on the autoregressive operator $\rho(\cdot)$ and errors $\{\varepsilon_{t}\}$ such that $\{X_t\}$ is strictly stationary\footnote{In particular, $\{\varepsilon_{t}\}$ is a strong $H$-white noise and the autoregressive operator $\rho$ satisfies $\Vert\rho^k\Vert_\mathcal{L}<1 \text{ for some }k\geq1$. For more details cf. Bosq2000.}. Then the one-step ahead forecast for process $\{X_t\}$ is given as
Hence, it is necessary to estimate operator $\rho\left(\cdot\right)$ in order to calculate the forecast $X_{T+1|T}$.
The estimation of the autoregressive operator $\rho(\cdot)$ faces the well known ill-posed inverse problem. Denote the covariance and the first autocovariance operators of process $\{X_t\}$ as $C(x)=E\left[\langle X_t,x\rangle X_t\right]$ and $\Gamma_1(x)=E\left[\langle X_t,x\rangle X_{t-1}\right]$, then from ((ref)) we have that operator equation $\Gamma_1=\rho C$ holds and formally gives the solution $\rho=\Gamma_1 C^{-1}$. However, the covariance operator $C$ does not have a bounded inverse on the entire space $L^2({\cal D})$. To see this point consider the spectral representation
where $\{\xi_{l}\}$ is the sequence of eigenvalues of $C(x)$ and $\{\phi_{\ell}\}$ is the sequence of the corresponding eigenfunctions. It follows that $C^{-1}(y)=\sum_{\ell=1}^\infty \xi_{l}^{-1}\langle \phi_{l},y\rangle \phi_{l}$, which is defined if all $\xi_{l}$ are positive. Further, $C^{-1}$ is unbounded since $\Vert C^{-1}(\phi_{l})\Vert=\xi_{l}^{-1}$, where $\xi_{l}^{-1}\to\infty$ as $l\to\infty$. As a consequence, estimating the bounded operator $\rho\left(\cdot\right)$ through the relationship $\rho = \Gamma_1 C^{-1}$ is difficult and some form of regularization has to be employed. One of the practical solutions in the literature to this problem is to use only the first $L$ summand in (ref). That is, we consider a truncated inverse covariance
To compute the final estimator of $\rho\left(\cdot\right)$ components of $\Gamma_1 C^{-1}$ are replaced by their sample counterparts accounting for the ill-posed problem as in (ref). To be more specific $\Gamma_1$ is estimated by $\widehat{\Gamma}_1=\frac{1}{T-1} \sum_{t=1}^{T-1} \langle {X}_t, x \rangle {X}_{t+1}$ and $C^{-1}$ is replaced by estimator of $C^{-1}_L$ given as $ \widehat{C}^{-1}_L=\sum_{l=1}^{L} \hat{\xi}_l^{-1} \langle x, \hat{\phi}_l \rangle \hat{\phi}_l$. This leads to the following expression of the estimator
On the final note to obtain the estimator $\widehat{\rho}(x)$ and the forecast $X_{T+1|T}$, the principle components $\{\widehat{\xi}_{l}\}$ and $\{\widehat{\phi}_{l}\}$ have to be estimated as discussed in Section (ref). Further, the series of surfaces $\{X_{t}\}$ is not observed directly and have to be reconstructed first. Hence, to calculate the final expression of $\widehat{\rho}(x)$ observations $\{X_{t}\}$ should be replaced with approximated ones $\{\widetilde{X}_{t}\}$, as discussed in Section (ref).
Two other standard benchmarks are commonly employed in FDA for comparative analysis (see e.g., didericksen2012). The first one is the mean predictor (MP), where the predictor is calculated as the mean of the sample,
and the second one in the naive predictor (NP) give as
In this section, we showcase the performance of the new framework in a forecasting study of ground-level ozone concentration. Ground level ozone is a harmful air pollutant as opposed to the stratospheric ozone located in the upper layers of the atmosphere which shields the earth from ultraviolet rays. The ground-level ozone is not directly emitted but formed as a byproduct of chemical reactions following the emission of primary air pollutants such as carbon monoxide, nitrogen oxide and methane. Its negative impact on human health and the biosphere in general is well documented by epidemiologists and toxicologists, which makes ozone one of the closely monitored and regulated air pollutants. For instance, the EU member countries are obliged to ensure and maintain compliance with targeted values of ozone set by the Ambient Air Quality Directives. Hence, reliable statistical tools are required by environmental agencies and regulators to monitor the air quality and to make sure the critical values are not exceeded. In what follows, we demonstrate that the predictive methodology proposed in this paper can serve this purpose as a preventive tool.
The main objective of this section is to provide forecasts of ozone concentration over the geographical region of Germany. This aspect is covered by evaluating the performance of our framework in comparison to the available alternatives in FDA literature (i.e., FAR model, mean predictor and naive predictor). We also show that the DFFM model is helpful not only in forecasting but also in estimating the regions that contribute most to the variation of ozone concentration levels.
The data that motivated this research consists of daily measurements of ozone concentration made available through AirBase - the European air quality database provided by the European Environment Agency (see this \href{http://www.eea.europa.eu/data-and-maps/data/airbase-the-european-air-quality-database-7#tab-data-by-country}{link} for more information and access to the raw data). This is a public database containing air quality monitoring information throughout Europe. We analyzed raw data of daily ozone concentration for Germany measured in $\mu g/m^3$, which consists of measurements at $1656$ stations dating back as far as $1984/01/01$. However, not all stations have been operated continuously and continuous measurements are available at $N=171$ stations for the year $2011$. As the raw data set was very large (ca. $5$GB), consisting of roughly $20,000$ text files that also contain recordings of other air pollutants, the data was first parsed with a Python script to generate a dataset that could be analyzed further.
In the first step, we reconstructed the sample of $T=365$ surfaces $( \widetilde{X}_t )_{t=1}^T$ from discrete observations at sample measurement stations using the FEM spline approach presented in section (ref). To deal with the complicated shape of the domain dictated by Germany's geographical borders, we obtain the Delaunay triangular mesh as the convex hull of the 171 measurement stations. As some of the triangles cover the exterior of Germany they were manually removed. See Figure (ref) where we plot locations of the stations (left panel), the Delaunay mesh (middle panel) and the final cleaned mesh (right panel).
Once the triangulation is obtained we proceed with the minimization problem (ref). The value of the smoothing parameter $\lambda$ has been determined by minimizing the generalized cross-validation criterion presented in (ref) which resulted in a parameter value of $\lambda=0.01$. Results in this section were obtained using the fda package in Matlab.\footnote{https://www.psych.mcgill.ca/misc/fda/downloads/FDAfuns/} The inspection of the obtained time series of surfaces reveal 28 episodes, in which ozone concentration exceeded the WHO standards of 100 $\mu g/m^3$ in some regions of Germany. The highest level is observed on May 31 in the region surrounding Dresden reaching 123.74 $\mu g/m^3$. Figure (ref) provides a graphical illustration of the corresponding surface.
In the second part of our analysis, we proceed with modelling and forecasting ozone surfaces. We first estimate the factors loading surfaces in (ref) as well as their number. This step allows us not only to understand and adequately describe the dynamic nature of our data set (necessary for forecasting) but also localize areas that contribute to the variability of the ozone layer over Germany. To be more precise, while factors and their dynamics are useful for forecasting, the loading surfaces supply us with insights on the spatial variability of ozone concentration. For this purpose, we implement the information criterion developed in OttoSalish2021 and aue2015. Results are reported in Table (ref) for the full sample and a subsample used for model training (see details below). $\widehat K_{\text{BIC}}$ and $\widehat p_{\text{BIC}}$ are the estimated numbers of factors and lags from the BIC variant OttoSalish2021 information criterion, $\widehat K_{\text{HQ}}$ and $\widehat p_{\text{HQ}}$ are estimated quantities from the HQ variant OttoSalish2021 information criterion, and $\widehat K_{\text{ANH}}$ and $\widehat p_{\text{ANH}}$ obtained from aue2015.
We conclude that in all cases at least the first two common components are required to capture the dynamics of the data. These two components combined explain 89% of the variability in the data (with the first one explaining around 82%). A closer look at the loadings surfaces plotted in Figure (ref), reveal that the areas that contribute the most to ozone variability are located in the belt connecting Frankfurt and Berlin according to the first loading surface. The second loading discriminates contribution to ozone concentration between south-west Germany vs north-east.
To evaluate the quality of forecasts obtained by the DFFM model the following scenario was considered. We split the available data set into two parts: the first 200 observations are used for model training and estimation and the remaining 165 observations for assessment of the forecasting performance. We proceed recursively. That is, first we produce a 1-step ahead forecast from the training sample of 200 observations. Then the training sample is updated to 201 observations to obtain the next 1-step ahead forecast and repeat this process until the end of the sample. In each step when the training sample is updated all parameters of the DFFM model are reestimated, which includes factors, loading functions, their number, number of lags and number of neighbours. This gives us 165 forecasts for ozone concentration surface, quality of which are assessed by the mean squared error:
where $\widetilde{X}_\tau(\boldsymbol{s}_i)$ denotes surface over Germany, $\widehat{X}_\tau(\boldsymbol{s}_i)$ denotes a forecast of $\widetilde{X}_\tau(\boldsymbol{s}_i)$ and $H$ denotes the number of obtained forecasts. To provide a comparative analysis we also include all standard benchmark models reviewed in Section (ref).
The results are reported in Figure (ref) in the form of box plots of RMSE. The results for the mean predictor are excluded from this figure since they were comparably large, making a graphical comparison of the remaining models impossible. (see Figure (ref) in the appendix for an illustration of the mean predictor performance). We report that the DFFM model with an K-NN forecasting framework has the smallest dispersion of the RMSE and the lowest average RMSE, making it the most attractive model for forecasting ozone concentration over Germany. To be more precise, the average MSE of the DFFM model with K-NN is 89.83, whereas the closest competitor, the DFFM model with ANH information criterion, has MSE equal to 103.53.
Furthermore, the part of the sample reserved for the forecast comparison contains two events when WHO standards for ozone concentration were exceeded in the area surrounding Dresden on August 24 and August 26. It is interesting to note that both the DFFM and the FAR models are able to predict well these events and their locations. Figure (ref) illustrates both occasions with the real surface and forecasted ones. We have reported only the DFFM model with the K-NN method as the other variants of the DFFM model provide similar outcomes. See Figures (ref) and (ref) in the appendix for a complete picture.
In this paper, we contribute to the literature on environmental data analysis. Given the mounting evidence of the manifold effects that climate change has almost on every aspect of life from health over biodiversity to the economy, new tools for processing environmental data collected over some geographical region and over time is becoming increasingly demanded.
To address this problem, we introduce a new methodological approach to environmental data based on new techniques available in Functional Data Analysis. This new approach models simultaneously the spatial and temporal aspects of data and consequently helps to prevent a potential loss of information concealed, for instance, by aggregation over time and/or space. We have illustrated the usefulness of this modelling perspective in the application to ground-level ozone and found several interesting findings. First, the new approach provides more accurate predictions than the available benchmark models. Second, it predicts well local areas, where air quality standards are exceeded. In addition, it provides an estimation of the areas that contribute most to the variability of ozone surfaces.
This paper focuses on a single air pollutant. A particular interesting direction for future research is the extension of our modelling framework to the multivariate settings, where several major air pollutants are modelled and monitored simultaneously. A potential difficulty in this regard is the theoretical background and practical implementation of a multivariate dynamic functional factor model. Although there are no results available yet in this direction a promising theoretical development was recently proposed in tavakoli2021. Another attractive direction of research is an application of the new methodology to different environmental data sets such as water pollution and insolation maps. It will be of interest to see if our modelling approach can provide additional information or improve forecasting performance when compared to the existing approaches.