# Forecasting Environmental Data: An example to ground-level ozone
concentration surfaces

Alexander<EMAIL_ADDRESS>[ Nazarii Salish label=e2
[<EMAIL_ADDRESS>Cognite AS  Universidad Carlos III de Madrid,

###### Abstract

Environmental problems are receiving increasing attention in socio-economic
and health studies. This in turn fosters advances in recording and data
collection of many related real-life processes. Available tools for data
processing are often found too restrictive as they do not account for the rich
nature of such data sets. In this paper, we propose a new statistical
perspective on forecasting spatial environmental data collected sequentially
over time. We treat this data set as a surface (functional) time series with a
possibly complicated geographical domain. By employing novel techniques from
functional data analysis we develop a new forecasting methodology. Our
approach consists of two steps. In the first step, time series of surfaces are
reconstructed from measurements sampled over some spatial domain using a
finite element spline smoother. In the second step, we adapt the dynamic
functional factor model to forecast a surface time series. The advantage of
this approach is that we can account for and explore simultaneously spatial as
well as temporal dependencies in the data. A forecasting study of ground-level
ozone concentration over the geographical domain of Germany demonstrates the
practical value of this new perspective, where we compare our approach with
standard functional benchmark models.

Functional Time Series,

Factor Models,

Environmental Data,

###### keywords:

,

##  1 Introduction

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., Aue, Norinho
and Hörmann, 2015, Mu, Wang and Wang, 2018), water quality control (see, e.g.,
Cañedo-Argüelles et al., 2016 and Estévez et al., 2019), temperature (see,
e.g., Kuenzer, Hörmann and Kokoszka, 2021) 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, Ramsay and Silverman, 2005 and Kokoszka and Reimherr, 2017 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, Ramsay
(2002), Sangalli, Ramsay and Ramsay (2013), Ettinger, Perotto and Sangalli
(2016), Mu, Wang and Wang (2018) and Ferraccioli et al. (2021). In the last
decade, we have also seen a rise of literature addressing modelling and
forecasting of functional time series (see e.g., Bosq, 2000, Hyndman and
Shang, 2009, Hörmann and Kokoszka, 2010, Hörmann and Kokoszka, 2012, and
Salish and Gleim, 2019). 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., Delicado et al., 2010). 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, Giraldo, Delicado and Mateu (2010), Giraldo,
Delicado and Mateu (2011) and Nerini, Monestiez and Manté (2010) extend the
methodology of kriging from classical geostatistics to the functional setting
while Yamanishi and Tanaka (2003) 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

$\displaystyle\left\\{X_{t}(\boldsymbol{s}):\boldsymbol{s}\in{\cal D},\ t\geq
1\right\\}.$ (1)

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 1).

Figure 1: Spatial locations of measurement stations and triangular mesh. The
left panel shows the geographical locations of ozone measurement stations
together with the border of Germany. The middle panel shows the Delaunay
triangulation where the nodes are placed at ozone sample stations. The right
panel shows the corrected Delaunay triangular mesh where triangles that cover
the exterior of Germany have been removed.

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 (1)
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 Ramsay (2002) and Sangalli, Ramsay and Ramsay (2013). 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., Hays, Shen and Huang, 2012, Liebl, 2013, Aue, Norinho and Hörmann, 2015
and Otto and Salish, 2022). 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 3 for a more detailed
description). As discussed and shown in Aue, Norinho and Hörmann (2015) and
Otto and Salish (2022) 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 Ramsay and
Silverman (2005) 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 Otto and Salish (2022) 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
Agency111https://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.222https://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 2 describes how
smooth surfaces can be obtained from noisy discrete spatial measurements by
using a finite element spline smoother. Section 3 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 4. Section 5 concludes.

##  2 Surface Reconstruction

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 Ramsay (2002), Sangalli, Ramsay and Ramsay
(2013) and Ferraccioli et al. (2021).

The spline estimate $\widetilde{X}_{t}$ of $X_{t}$ is defined as a minimizer
of the following functional

$\displaystyle J_{\lambda}({\cal
X})=\sum_{i=1}^{N}\left(X^{*}_{t}({s}_{i})-{\cal
X}({s}_{i})\right)^{2}+\lambda\int_{{\cal D}}\left(\Delta{\cal
X}\right)^{2}\text{ for }t=1,...,T,$ (2)

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 (2) 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., Folland, 1995 or Ramsay, 2002), where for
any ${\cal X}\in H^{2}({\cal D})$

$\displaystyle\Delta{\cal X}(\boldsymbol{s})=\Delta{\cal
X}(x,y):=\frac{\partial^{2}{\cal X}}{\partial
x^{2}}(x,y)+\frac{\partial^{2}{\cal X}}{\partial y^{2}}(x,y).$

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 Sangalli, Ramsay and Ramsay (2013) a unique solution to the
minimization problem in (2) 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:

  * Step 1.

Partition the spatial domain ${\cal D}$ into disjoint sets;

  * Step 2.

Construct separate polynomial functions (of second order) on each piece of
partition, which forms finite elements.

  * Step 3.

Use the union of functions constructed in Step 2 to closely approximate the
solution to (2).

_Implementation 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 1 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}}$.

Figure 2: Reference triangular element. This figure shows a reference triangle
with corresponding nodal points of a triangular finite reference element.

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 2. 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 3.

Figure 3: Quadratic finite element shape functions. This figure shows the six
shape functions of a triangular finite reference element.

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 4 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.

Figure 4: Quadratic finite element basis function. This figure shows a plot of
a quadratic finite element basis function associated with the nodal point
$(0,0)$.

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

$\displaystyle\boldsymbol{A}_{K}$ $\displaystyle:=\int_{\triangle_{{\cal
D}}}\left(\boldsymbol{\phi}_{K}\right)\left(\boldsymbol{\phi}_{K}\right)^{\mathsf{T}},$
$\displaystyle\boldsymbol{B}_{K}$ $\displaystyle:=\int_{\triangle_{{\cal
D}}}\left(\boldsymbol{\phi}_{K}^{(x)}\right)\left(\boldsymbol{\phi}_{K}^{(x)}\right)^{\mathsf{T}}+\left(\boldsymbol{\phi}_{K}^{(y)}\right)\left(\boldsymbol{\phi}_{K}^{(y)}\right)^{\mathsf{T}},$

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 (2) as

$\displaystyle\widetilde{X}_{t}(\boldsymbol{s})=\sum_{k=1}^{K}\widetilde{X}(\boldsymbol{\xi}_{k})\phi_{k}(\boldsymbol{s})=\left(\boldsymbol{\widetilde{X}}_{K}\right)^{\mathsf{T}}\boldsymbol{\phi}_{K}(\boldsymbol{s}).$

Hence, the minimization problem is equivalent to searching for a $K\times 1$
vector $\boldsymbol{\widetilde{X}}_{K}$. As shown in Sangalli, Ramsay and
Ramsay (2013) solving for $\boldsymbol{\widetilde{X}}_{K}$ is then equivalent
to solving the system of linear equations given by

$\displaystyle\left(\begin{array}[]{cc}-\boldsymbol{D}_{K}&\lambda\boldsymbol{B}_{K}\\\
\lambda\boldsymbol{B}_{K}&\lambda\boldsymbol{A}_{K}\end{array}\right)\left(\begin{array}[]{c}\boldsymbol{\widetilde{X}}_{K}\\\
\boldsymbol{Z}_{K}\end{array}\right)=\left(\begin{array}[]{c}-\boldsymbol{X}_{K}^{*}\\\
\boldsymbol{0}_{K}\end{array}\right),$ (9)

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).

######  Remark 1 (Choosing the smoothing parameter).

One method commonly used in the literature for choosing the smoothing
parameter $\lambda$ is to consider minimization of the generalized cross-
validation criterion given by

$\displaystyle
GCV(\lambda)=\frac{1}{N\left(1-\mathrm{tr}(\boldsymbol{S}_{N})/N\right)^{2}}\left(\boldsymbol{X}^{*}_{N}-\boldsymbol{S}_{N}\boldsymbol{X}^{*}_{N}\right)^{\mathsf{T}}\left(\boldsymbol{X}^{*}_{N}-\boldsymbol{S}_{N}\boldsymbol{X}^{*}_{N}\right),$
(10)

where $\boldsymbol{S}_{N}$ is $N\times N$ matrix corresponding to the first
$N$ rows and $N$ columns of $\boldsymbol{S}_{2K}$ and $\boldsymbol{S}_{2K}$ is
the inverse of the matrix on the left-hand side of (9), i.e.,

$\boldsymbol{S}_{2K}=\left(\begin{array}[]{cc}-\boldsymbol{D}_{K}&\lambda\boldsymbol{B}_{K}\\\
\lambda\boldsymbol{B}_{K}&\lambda\boldsymbol{A}_{K}\end{array}\right)^{-1}.$

(See e.g., Ramsay and Silverman, 2005 for further details on this topic).

##  3 Forecasting Surfaces Time Series

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

$X_{t}(s)=\mu(s)+\sum_{l=1}^{L}{x}_{l,t}\psi_{k}(s)+\epsilon_{t}(s)\quad
s\in{\cal D},$ (11)

where a $L\times 1$ vector of factors
$\mathbf{x}_{t}=[x_{1,t},...,x_{L,t}]^{\prime}$ 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)]^{\prime}$, are unobserved
deterministic terms. The product $\Psi^{\prime}\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 (11) 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
$\|x\|=\langle x,x\rangle^{1/2}$.

######  Assumption 1.

x

  * (a)

The loading surfaces $\\{\psi_{l}\\}_{l=1}^{L}$ form an orthonormal system,
that is, $\langle\psi_{m},\psi_{l}\rangle=0$ and $\|\psi_{l}\|=1$, for all
$m,l=1,\ldots,L$ with $m\neq l$;

  * (b)

The factors satisfy $E[\mathbf{x}_{t}]=0$,
$E[\mathbf{x}_{t}\mathbf{x}_{t}^{\prime}]=\text{diag}(\lambda_{1},\ldots,\lambda_{L})$,
where $\lambda_{1}>\ldots>\lambda_{L}>0$, and
$E\left[\|\mathbf{x}_{t}\|_{2}^{4}\right]<\infty$;

  * (c)

$E\left[x_{L,t}x_{l,t-i}\right]\neq 0$ for some $l=1,\ldots,L$ and
$i\in\mathbb{N}$;

######  Assumption 2.

For any $r,s\in{\cal D}$, $l=1,\ldots,L$, and $t$, the following conditions
hold true:

  * (a)

$\\{\epsilon_{t}\\}$ is a $H$-martingale difference sequence, that is,
$\\{\epsilon_{t}\\}$ is adapted to the natural filtration
$\mathcal{A}_{\tau}=\sigma\left(\epsilon_{\tau},\ \tau\leq t\right)$ with
$E[\epsilon_{t}(r)|\ \mathcal{A}_{t-1}]=0$, $E[\epsilon_{t}(r)\epsilon_{t}(s)\
|\ \mathcal{A}_{t-1}]=\delta(r,s)$, and $\sup_{r,s\in{\cal
D}}|\delta(r,s)|<\infty$. Furthermore, the eigenvalues $\\{\zeta_{l}\\}$ of
the integral operator with kernel $\delta(r,s)$ satisfy
$\zeta_{l}>\zeta_{l+1}$ for all $l\in\mathbb{N}$, and
$E\left\|\epsilon_{t}\right\|^{\kappa}<\infty$ for some $\kappa>4$;

  * (b)

$E\big{[}\epsilon_{t}(r)\langle\epsilon_{t},\psi_{l}\rangle\big{]}=0$;

  * (c)

$E\big{[}\langle\epsilon_{t},x\rangle^{2}\big{]}<\lambda_{K}\|x\|$ for all
$x\in H$;

  * (d)

The common and the idiosyncratic component are weakly dependent, where

$\lim_{T\to\infty}\sup_{s\in{\cal
D}}\sup_{s\in\mathbb{N}}E\left\|\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathbf{x}_{t}\epsilon_{t-s}(s)\right\|_{2}^{2}<\infty.$

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 Otto and Salish, 2022 for the
theoretical discussion of this result). This, in turn, indicates that all
latent components of model (11) can be easily estimated. Section 3.1 guides
the practical implementation of this step. Further, Assumptions 1 and 2 do not
impose a particular modelling structure on the factors allowing for some
flexibility when choosing a forecasting approach. In Section 3.2 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}$.

###  3.1 Estimation of latent components

The estimation of all latent components in model (11), 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

$\widehat{\mu}(s)=\frac{1}{T}\sum_{t=1}^{T}X_{t}(s),\quad s\in{\cal D},$

and the sample covariance function as

$\widehat{c}(s,r)=\frac{1}{T}\sum_{t=1}^{T}\big{(}X_{t}(s)-\widehat{\mu}(s)\big{)}\big{(}X_{t}(r)-\widehat{\mu}(r)\big{)},\quad
s,r\in{\cal D}.$

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 Otto and Salish (2022), Theorem 2 these
estimators provide consistent estimates. For practical implementation of this
step, we refer a reader to “fda” package of Ramsay and Silverman (2005).

###  3.2 Forecasting Methodology

The best $h$-step ahead predictor (in the mean square error sense) of the
process of interest $X_{t}$ is given as,

$X_{T+h|T}(s)=E\big{[}X_{T+h}(s)|X_{T},X_{T-1},...\big{]}=\mu(s)+\Psi(s)^{\prime}\mathbf{x}_{T+h|T},$

where $\mathbf{x}_{T+h|T}$ denotes an $h$-step ahead predictor of the factors,
which in turn is given as

$\mathbf{x}_{T+h|T}=E\big{[}\mathbf{x}_{T+h}|\mathbf{x}_{T},\mathbf{x}_{T-1},\ldots\big{]},$
(12)

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 3.1 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, Lütkepohl, 2005 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}$

$\mathbf{x}_{T+1|T}=\mathbf{A}_{1}\mathbf{x}_{T}+,...,+\mathbf{A}_{p}\mathbf{x}_{T-p+1},$
(13)

where $\mathbf{A}_{1},...,\mathbf{A}_{p}$ denote the $L\times L$ coefficient
lag matrices and $p$ is the number of relevant lags. Predictor (13) 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 3.1, 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 (13).

Estimators of $L$ and $p$ are obtained by employing the information criterion
from Aue, Norinho and Hörmann (2015) (ANH) or one of the variants of the
information criterion proposed in Otto and Salish (2022) (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., Lütkepohl, 2005,
Section 4.3.3). 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:

  * Step 1.

Use reconstructed surfaces $\\{\widetilde{X}_{t}\\}_{t=1}^{T}$ to compute
estimates of latent factors $\widehat{x}_{l,t}$ (as described in Section 3.1);

  * Step 2.

Use estimated latent factors $\\{\widehat{x}_{l,t}\\}$ and one of the
available information criteria (e.g., ANH or OS) to estimate the number of
factors, $L$, and the number of lags, $p$.

  * Step 3.

Estimate VAR model of order $\widehat{p}$ for $\widehat{L}\times 1$
dimensional time series of vectors
$\widehat{\mathbf{x}}_{t}{(\widehat{L})}=\left[\widehat{x}_{1,t},...,\widehat{x}_{\widehat{L},t}\right]^{\prime}$
(with $\widehat{p}$ and $\widehat{L}$ obtained in Step 2). Use estimated
matrices $\widehat{\mathbf{A}}_{1},...,\widehat{\mathbf{A}}_{p}$ to estimate
(13) as:

$\widehat{\mathbf{x}}_{T+1|T}=\widehat{\mathbf{A}}_{1}\widehat{\mathbf{x}}_{T}(\widehat{L})+,...,+\widehat{\mathbf{A}}_{p}\widehat{\mathbf{x}}_{T-\widehat{p}+1}(\widehat{L}).$

_Non-Linear case._ When there is little a priori knowledge about the shape of
(12) 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 (13) 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., Györfi, 2002, Chapter 6 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 (13) with

$\mathbf{x}_{T+1|T}=m(\mathbf{x}_{T},...,\mathbf{x}_{T-p+1}).$ (14)

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 Biau, Cerou and Guyader, 2010, Kudraszow and Vieu,
2013 and Kara et al., 2017) we can adopt techniques already available for
classical time series (see e.g., Yakowitz, 1987). The estimator of $m(\cdot)$
and therefore the estimator of the forecast $\mathbf{x}_{T+1|T}$ is obtained
following the next basic steps.

  * Step 1.

Use reconstructed surfaces $\\{\widetilde{X}_{t}\\}_{t=1}^{T}$ to compute
estimates of latent factors $\widehat{x}_{l,t}$ (as described in Section 3.1);

  * Step 2.

Determine by cross-validation the relevant values for the number of
neighbours, $K$, $p$ and $L$. To increase computational speed $p$ and $L$ can
be approximated using one of the available information criteria ANH or OS.
Following theoretical results in Yakowitz (1987) we set the maximum number of
$K_{max}=T^{4/5}$ in a search for the optimal (in MSE sense) number of
neighbours, denoted as $\widehat{K}$.

  * Step 3.

Compute distances between most recent data point vector
$\mathbf{z}_{T}\equiv[\widehat{\mathbf{x}}_{T}(\widehat{L}),...,\widehat{\mathbf{x}}_{T-\widehat{p}+1}(\widehat{L})]^{\prime}$
and each vector in the rest of sample
$\\{\mathbf{z}_{t}\\}_{t=\widehat{p}}^{T-1}$. Here,
$\mathbf{z}_{t}=[\widehat{\mathbf{x}}_{t}(\widehat{L}),...,\widehat{\mathbf{x}}_{t-\widehat{p}+1}(\widehat{L})]^{\prime}$
and
$\widehat{\mathbf{x}}_{t}{(\widehat{L})}=\left[\widehat{x}_{1,t},...,\widehat{x}_{\widehat{L},t}\right]^{\prime}$,
were $\widehat{p}$ and $\widehat{L}$ are determined in Step 2. The Minkowsky
distance is a typical choice to measure a distance between the feature vector
$\mathbf{z}_{T}$ and the rest of vectors in the sample. Denote the index set
of $\widehat{K}$ closest neighbors to the feature vector $\mathbf{z}_{T}$ by
$\mathcal{I}_{T}$.

  * Step 4.

Once the $\widehat{K}$ closest neighbours of $\mathbf{z}_{T}$ are identified
their subsequent values are used to obtain the final estimator

$\widehat{\mathbf{x}}_{T+1|T}=\widehat{m}(\widehat{\mathbf{x}}_{T},...,\widehat{\mathbf{x}}_{T-p+1})=\frac{1}{\widehat{K}}\sum_{t\in\mathcal{I}_{T}}w_{t}\
\widehat{\mathbf{x}}_{t+1}(\widehat{L}),$ (15)

where $\\{w_{t}\\}_{t\in\mathcal{I}_{T}}$ denote weights. For instance, the
KNN estimator (15) can be used with equal weights $w_{i}=1/\widehat{K}$ or
weights can be set inversely proportional to the distance between the feature
vector $z_{T}$ and neighbours.

###  3.3 Benchmark Forecasting Methods in FDA

The key method in functional data analysis used for forecasting is the
functional autoregressive model (FAR). Since the seminal monograph by Bosq
(2000) this model has seen a lot of development and have become the main
toolbox for the analysis of functional time series. (See, e.g., Besse, Cardot
and Stephenson, 2000, Mas, 2007, Kargin and Onatski, 2008, Kokoszka and
Reimherr, 2013 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 Horváth and Kokoszka, 2012 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

$X_{t}=\rho\left(X_{t-1}\right)+\varepsilon_{t}.$ (16)

Model (16) satisfies standard regularity conditions on the autoregressive
operator $\rho(\cdot)$ and errors $\\{\varepsilon_{t}\\}$ such that
$\\{X_{t}\\}$ is strictly stationary333In particular, $\\{\varepsilon_{t}\\}$
is a strong $H$-white noise and the autoregressive operator $\rho$ satisfies
$\|\rho^{k}\|_{\mathcal{L}}<1\text{ for some }k\geq 1$. For more details cf.
Bosq, 2000, Theorem 3.1.. Then the one-step ahead forecast for process
$\\{X_{t}\\}$ is given as

$X_{T+1|T}=E\left[X_{T+1}|X_{T},X_{T-1},...\right]=\rho\left(X_{T}\right).$

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 (16) 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

$\displaystyle C(x)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{\infty}\xi_{l}\langle\phi_{l},x\rangle\phi_{\ell},$
(17)

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 $\|C^{-1}(\phi_{l})\|=\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 (17). That is, we consider a truncated inverse covariance

$C^{-1}_{L}=\sum_{l=1}^{L}\xi_{l}^{-1}\langle\phi_{l},x\rangle\phi_{l}.$ (18)

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 (18). 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

$\widehat{\rho}(x)=\frac{1}{T-1}\sum_{t=1}^{T-1}\sum_{k,l=1}^{L}\widehat{\xi}_{l}^{-1}\langle
x,\widehat{\phi}_{l}\rangle\langle X_{t},\widehat{\phi}_{l}\rangle\langle
X_{t+1},\widehat{\phi}_{k}\rangle\widehat{\phi}_{k}.$

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 3.1.
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 2.

Two other standard benchmarks are commonly employed in FDA for comparative
analysis (see e.g., Didericksen, Kokoszka and Zhang, 2012). The first one is
the _mean predictor_ (MP), where the predictor is calculated as the mean of
the sample,

$\widehat{X}_{T+1|T}=\sum_{t=1}^{T}\widetilde{X}_{t},$ (19)

and the second one in the _naive predictor_ (NP) give as

$\widehat{X}_{T+1|T}=\widetilde{X}_{T}.$ (20)

##  4 An Application to Ozone Concentration Surfaces

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.

###  4.1 Data

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 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.

###  4.2 Surface Reconstruction

Figure 5: Ozone concentration surface. This figure shows ozone surfaces over
Germany containing the highest value of ozone concentration in the
reconstructed time series of surfaces, observed on May 31 in the east of
Germany.

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 2. 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 1 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
(2). The value of the smoothing parameter $\lambda$ has been determined by
minimizing the generalized cross-validation criterion presented in (10) which
resulted in a parameter value of $\lambda=0.01$. Results in this section were
obtained using the fda package in
Matlab.444https://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 5 provides a
graphical illustration of the corresponding surface.

###  4.3 Modeling and Forecasting

In the second part of our analysis, we proceed with modelling and forecasting
ozone surfaces. We first estimate the factors loading surfaces in (11) 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 Otto and Salish (2022) and Aue, Norinho and
Hörmann (2015). Results are reported in Table 1 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 Otto and Salish (2022)
information criterion, $\widehat{K}_{\text{HQ}}$ and $\widehat{p}_{\text{HQ}}$
are estimated quantities from the HQ variant Otto and Salish (2022)
information criterion, and $\widehat{K}_{\text{ANH}}$ and
$\widehat{p}_{\text{ANH}}$ obtained from Aue, Norinho and Hörmann (2015).

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 6,
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.

Table 1: Estimated numbers of factors and lags | $K_{\text{BIC}}$ | $K_{\text{HQ}}$ | $K_{\text{ANH}}$ | $p_{\text{BIC}}$ | $p_{\text{HQ}}$ | $p_{\text{ANH}}$
---|---|---|---|---|---|---
(full data) | 2 | 5 | 6 | 3 | 3 | 3
(first 200 days) | 2 | 2 | 5 | 3 | 3 | 3

Figure 6: Estimated Loading surfaces. The right panel plots the first
estimated loading surface and the left plots the second one.

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:

$MSE={\frac{1}{N}\sum_{i=1}^{N}\left(\widetilde{X}_{\tau}(\boldsymbol{s}_{i}-\widehat{X}_{\tau}(\boldsymbol{s}_{i}))\right)^{2}},$

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 3.3.

Figure 7: Boxplots of the prediction MSE. DFFM(BIC) - denotes the DFFM model
where parameters where selected with the BIC variant of the OS information
criterion; DFFM(HQ) - denotes the DFFM model with the HQ variant of the OS
information criterion; DFFM(ANH) - denotes the DFFM model with the BIC variant
of the ANH information criterion; DFFM(KNN) - denotes the DFFM model with K-NN
method; FAR -stands for the Functional Autoregression and NP for the Naive
Predictor.

The results are reported in Figure 7 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 9 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 8 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 10 and 11 in the appendix for a complete picture.

Figure 8: Events with exceeded WHO standards of ozone concentration and their
forecasts Top row illustrates the corresponding figures on August 24 and the
bottom row on August 26.

##  5 Conclusion and discussion

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
Tavakoli, Nisol and Hallin (2021). 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.

The second author gratefully acknowledges the financial support from Juan de
la Cierva Incorporación, IJC2019-041742-I.

## References

  * Aue, Norinho and Hörmann (2015) [author] Aue, AlexanderA., Norinho, Diogo DubartD. D. and Hörmann, SiegfriedS. (2015). On the prediction of stationary functional time series. Journal of the American Statistical Association 110 378–392.
  * Besse, Cardot and Stephenson (2000) [author] Besse, Philippe CP. C., Cardot, HervéH. and Stephenson, David BD. B. (2000). Autoregressive forecasting of some functional climatic variations. Scandinavian Journal of Statistics 27 673–687.
  * Biau, Cerou and Guyader (2010) [author] Biau, G.G., Cerou, F.F. and Guyader, A.A. (2010). Rates of convergence of the functional k-nearest neighbor estimate. IEEE Transactions on Information Theory 56 2034-2040.
  * Bosq (2000) [author] Bosq, D.D. (2000). Linear Processes in Function Spaces. Springer, New York.
  * Cañedo-Argüelles et al. (2016) [author] Cañedo-Argüelles, M.M., Hawkins, C. P.C. P., Kefford, B. J.B. J., Schäfer, R. B.R. B., Dyack, B. J.B. J., Brucet, S.S., Buchwalter, D.D., Dunlop, J.J., Frör, O.O., Lazorchak, J.J., Coring, E.E., Fernandez, H. R.H. R., Goodfellow, W.W., Achem, A. L. GonzálezA. L. G., Hatfield-Dodds, S.S., Karimov, B. K.B. K., Mensah, P.P., Olson, J. R.J. R., Piscart, C.C., Prat, N.N., Ponsá, S.S., Schulz, C. J.C. J. and Timpano, A. J.A. J. (2016). Saving freshwater from salts. Science 351 914-916. 10.1126/science.aad3488
  * Delicado et al. (2010) [author] Delicado, P.P., Giraldo, R.R., Comas, C.C. and Mateu, J.J. (2010). Statistics for spatial functional data: some recent contributions. Environmetrics 21 224–239.
  * Didericksen, Kokoszka and Zhang (2012) [author] Didericksen, DevinD., Kokoszka, PiotrP. and Zhang, XiX. (2012). Empirical Properties of Forecasts with the Functional Autoregressive Model. Computational Statistics 27 285–298.
  * Estévez et al. (2019) [author] Estévez, EdurneE., Rodríguez-Castillo, TamaraT., González-Ferreras, Alexia MaríaA. M., Cañedo-Argüelles, MiguelM. and Barquín, JoséJ. (2019). Drivers of spatio-temporal patterns of salinity in Spanish rivers: a nationwide assessment. Philosophical Transactions of the Royal Society B: Biological Sciences 374 20180022\. 10.1098/rstb.2018.0022
  * Ettinger, Perotto and Sangalli (2016) [author] Ettinger, B.B., Perotto, S.S. and Sangalli, L. M.L. M. (2016). Spatial regression models over two-dimensional manifolds. Biometrika 103 71-88. 10.1093/biomet/asv069
  * Ferraccioli et al. (2021) [author] Ferraccioli, FedericoF., Arnone, EleonoraE., Finos, LivioL., Ramsay, James O.J. O. and Sangalli, Laura M.L. M. (2021). Nonparametric density estimation over complicated domains. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 83 346-368. https://doi.org/10.1111/rssb.12415
  * Folland (1995) [author] Folland, G. B.G. B. (1995). Introduction to Partial Differential Equations. Princeton: Princeton University Press.
  * Giraldo, Delicado and Mateu (2010) [author] Giraldo, R.R., Delicado, P.P. and Mateu, J.J. (2010). Continuous time-varying kriging for spatial prediction of functional data: An environmental application. Journal of Agricultural, Biological and Environmental Statistics 15 66–82.
  * Giraldo, Delicado and Mateu (2011) [author] Giraldo, R.R., Delicado, P.P. and Mateu, J.J. (2011). Ordinary kriging for functional-valued spatial data. Environmental and Ecological Statistics 18 411–426.
  * Györfi (2002) [author] Györfi, Kohler M. Krzyzak A. Walk H.K. M. K. A. W. H. L. (2002). A Distribution-Free Theory of Nonparametric Regression. Springer Series in Statistics, Springer-Verlag, New York.
  * Hays, Shen and Huang (2012) [author] Hays, SpencerS., Shen, HaipengH. and Huang, Jianhua ZJ. Z. (2012). Functional dynamic factor models with application to yield curve forecasting. The Annals of Applied Statistics 6 870–894.
  * Hörmann and Kokoszka (2010) [author] Hörmann, S.S. and Kokoszka, P.P. (2010). Weakly dependend functional data. The Annals of Statistics 38 1845-1884.
  * Horváth and Kokoszka (2012) [author] Horváth, L.L. and Kokoszka, P.P. (2012). Inference for Functional Data with Applications. Springer, New York.
  * Hörmann and Kokoszka (2012) [author] Hörmann, SiegfriedS. and Kokoszka, PiotrP. (2012). 7 - Functional Time Series. In Time Series Analysis: Methods and Applications, (TataT. Subba Rao, SuhasiniS. Subba Rao and C. R.C. R. Rao, eds.). Handbook of Statistics 30 157-186. Elsevier. https://doi.org/10.1016/B978-0-444-53858-1.00007-7
  * Hyndman and Shang (2009) [author] Hyndman, Rob J.R. J. and Shang, Han LinH. L. (2009). Forecasting functional time series. Journal of the Korean Statistical Society 38 199-211. https://doi.org/10.1016/j.jkss.2009.06.002
  * Kara et al. (2017) [author] Kara, Lydia-ZaitriL.-Z., Laksaci, AliA., Rachdi, MustaphaM. and Vieu, PhilippeP. (2017). Data-driven kNN estimation in nonparametric functional data analysis. Journal of Multivariate Analysis 153 176-188. https://doi.org/10.1016/j.jmva.2016.09.016
  * Kargin and Onatski (2008) [author] Kargin, VladislavV. and Onatski, AlexeiA. (2008). Curve forecasting by functional autoregression. Journal of Multivariate Analysis 99 2508–2526.
  * Kokoszka and Reimherr (2013) [author] Kokoszka, PiotrP. and Reimherr, MatthewM. (2013). Determining the order of the functional autoregressive model. Journal of Time Series Analysis 34 116–129. 10.1111/j.1467-9892.2012.00816.x
  * Kokoszka and Reimherr (2017) [author] Kokoszka, PiotrP. and Reimherr, MatthewM. (2017). Introduction to Functional Data Analysis. Chapman and Hall/CRC, Boca Raton.
  * Kudraszow and Vieu (2013) [author] Kudraszow, Nadia L.N. L. and Vieu, PhilippeP. (2013). Uniform consistency of kNN regressors for functional variables. Statistics & Probability Letters 83 1863-1870.
  * Kuenzer, Hörmann and Kokoszka (2021) [author] Kuenzer, ThomasT., Hörmann, SiegfriedS. and Kokoszka, PiotrP. (2021). Principal Component Analysis of Spatially Indexed Functions. Journal of the American Statistical Association 116 1444-1456. 10.1080/01621459.2020.1732395
  * Liebl (2013) [author] Liebl, DominikD. (2013). Modeling and forecasting electricity spot prices: A functional data perspective. The Annals of Applied Statistics 7 1562 – 1592. 10.1214/13-AOAS652
  * Lütkepohl (2005) [author] Lütkepohl, HelmutH. (2005). New introduction to multiple time series analysis. Springer.
  * Mas (2007) [author] Mas, AndréA. (2007). Weak convergence in the functional autoregressive model. Journal of Multivariate Analysis 98 1231 - 1261. http://dx.doi.org/10.1016/j.jmva.2006.05.010
  * Mu, Wang and Wang (2018) [author] Mu, JingruJ., Wang, GuannanG. and Wang, LiL. (2018). Estimation and inference in spatially varying coefficient models. Environmetrics 29 e2485. e2485 env.2485. https://doi.org/10.1002/env.2485
  * Nerini, Monestiez and Manté (2010) [author] Nerini, D.D., Monestiez, P.P. and Manté, C.C. (2010). Cokriging for spatial functional data. Journal of Multivariate Analysis 101 409–418.
  * Otto and Salish (2022) [author] Otto, SvenS. and Salish, NazariiN. (2022). Dynamic Factor Model for Functional Time Series: Identification, Estimation, and Prediction. arXiv:2201.02532.
  * Ramsay (2002) [author] Ramsay, T.T. (2002). Spline smoothing over difficult regions. Journal of the Royal Statistical Society 64 307–319.
  * Ramsay and Silverman (2005) [author] Ramsay, J. O.J. O. and Silverman, B. W.B. W. (2005). Functional Data Analysis, 2nd ed. Springer, New York.
  * Salish and Gleim (2019) [author] Salish, NazariiN. and Gleim, AlexanderA. (2019). A moment-based notion of time dependence for functional time series. Journal of Econometrics 212(2) 377–392.
  * Sangalli, Ramsay and Ramsay (2013) [author] Sangalli, L. M.L. M., Ramsay, J. O.J. O. and Ramsay, T. O.T. O. (2013). Spatial spline regression models. Journal of the Royal Statistical Society, Series B 75 681–703.
  * Tavakoli, Nisol and Hallin (2021) [author] Tavakoli, ShahinS., Nisol, GillesG. and Hallin, MarcM. (2021). Factor Models for High-Dimensional Functional Time Series.
  * Yakowitz (1987) [author] Yakowitz, S.S. (1987). Nearest-neighbour methods for time series analysis. Journal of Time Series Analysis 8 235-247.
  * Yamanishi and Tanaka (2003) [author] Yamanishi, Y.Y. and Tanaka, Y.Y. (2003). Geographically weighted functional multiple regression analysis: a numerical investigation. Journal of the Japanese Society of Computational Statistics 15 307–317.

### Supplementary Material

In this section we provide additional results from the empirical application

Figure 9: Boxplots of the prediction MSE. Additionally to the results reported
in Figure 7 includes mean predictor. Figure 10: Events with exceeded WHO
standards of ozone concentration and their forecasts on August 24 Figure 11:
Events with exceeded WHO standards of ozone concentration and their forecasts
on August 26


# Towards Quantum Computing Timelike Hadronic Vacuum Polarization and Light-
by-Light Scattering: Schwinger Model Tests

João Barata<EMAIL_ADDRESS>Physics Department, Brookhaven National
Laboratory, Upton, NY 11973, USA  Kazuki Ikeda<EMAIL_ADDRESS>Co-design Center for Quantum Advantage  Center for Nuclear Theory, Department
of Physics and Astronomy, Stony Brook University, Stony Brook, New York
11794-3800, USA  Swagato Mukherjee<EMAIL_ADDRESS>Physics Department,
Brookhaven National Laboratory, Upton, NY 11973, USA  Jonathan Raghoonanan
<EMAIL_ADDRESS>Physics Department, Brookhaven National Laboratory, Upton,
NY 11973, USA

###### Abstract

Hadronic vacuum polarization (HVP) and light-by-light scattering (HLBL) are
crucial for evaluating the Standard Model predictions concerning the muon’s
anomalous magnetic moment. However, direct first-principle lattice gauge
theory-based calculations of these observables in the timelike region remain
challenging. Discrepancies persist between lattice quantum chromodynamics
(QCD) calculations in the spacelike region and dispersive approaches relying
on experimental data parametrization from the timelike region. Here, we
introduce a methodology employing 1+1-dimensional quantum electrodynamics
(QED), i.e. the Schwinger Model, to investigate the HVP and HLBL. To that end,
we use both tensor network techniques, specifically matrix product states, and
classical emulators of digital quantum computers. Demonstrating feasibility in
a simplified model, our approach sets the stage for future endeavors
leveraging digital quantum computers.

##  I Introduction

The determination of the muon’s anomalous magnetic moment is one of the most
ambitious and important high-precision programs related to the Standard Model
(SM) of particle physics [1, 2]. The quantum corrections to the classical
(Dirac) term to the muon’s magnetic moment couple to all three sectors, QED,
Weak, and QCD, of the SM and, therefore, any discrepancies between experiment
and theory give a clear signal of beyond the SM physics. The comparisons
between the experimental measurements and the theoretical predictions are
still being actively studied.

The dominant SM contribution to the anomalous moment comes from QED, and it
has been computed to $\mathcal{O}(\alpha_{em}^{5})$, with low theoretical
uncertainty, where $\alpha_{em}$ is the QED coupling constant. Weak
corrections have been evaluated up to two loops but give a small numerical
contribution since they are suppressed by powers of the muon to the $W$ boson
mass ratio. By far the most complex and least controlled contributions
originate from the QCD sector, where intrinsically non-perturbative effects
come into play, see e.g. [3, 4]. At the leading order in $\alpha_{em}$, these
contributions come from the hadronic vacuum polarization (HVP), associated to
photons fluctuating into low energy hadronic states
($\mathcal{O}(\alpha_{em}^{2})$, see Fig. 1), and hadronic light-by-light
(HLBL) scattering, which involves photons combining into QCD loops
annihilating to light ($\mathcal{O}(\alpha_{em}^{3})$, see Fig. 1).

Figure 1: Left: HVP contribution to the muon’s magnetic moment. Right: The
HLBL contributions to the muon’s magnetic moment (bottom), and the generic
form of a matrix element contributing to the photon state (top).
Electromagnetic currents are explicitly indicated, while the QCD elements are
denoted by gray circles.

State of the art theoretical calculations of the QCD contributions involve
either the use of dispersive methods or euclidean lattice approaches [4].
Dispersive methods make use of the fact that a virtual photon decaying to
hadrons can be related to $e^{+}e^{-}$ annihilation to QCD states, and that
the hadronic spectral functions can be extracted from QCD decays. As a result,
this is an intrinsically data driven approach. In contrast, in lattice based
calculations one can directly extract correlators of electromagnetic currents,
which are related to the HVP and HLBL. Although lattice QCD is the only true
first-principle approach available, it is limited to computations of
correlators of space-like separated currents.

In this work we propose using quantum computers (QCs) to circumvent some of
the issues found in Euclidean lattice calculations. The main motivation to
using QCs is that they allow for direct computations in Minkowski time without
the presence of a sign problem for any spacetime separation. In the current
stage of QC hardware, however, the available quantum devices are still too
small and error-prone to effectively perform any meaningful calculations in
QCD, and any precision calculation, if possible, will only be achieved in some
far future [5, 6, 7]. As a result, here we illustrate the approach in 1+1-d
QED, i.e. the Schwinger model [8]. This theory shares several phenomenological
properties with 3+1-d QCD (see e.g. [9, 10, 11]), but it is technically
simpler to implement on the lattice [12, 13]. In particular, as we detail
below, it can be efficiently analyzed using (classical) tensor network
methods, see e.g. [14], as a precursor to implementations in large scale
digital QCs [15, 16]. We recall that, despite the classical nature of tensor
networks, their implementation mimics several aspects faced when using QCs,
and thus both approaches are intimately connected, as we illustrate below.

Our calculations are performed in the massive Schwinger model; for the HVP the
main object of our computation is the current-current correlator (see Fig. 1)

$\displaystyle\Pi^{\mu\nu}(\vec{x},\vec{x}_{0})\equiv\langle
TJ^{\mu}(t,x)J^{\nu}(t_{0},x_{0})\rangle\,,$ (1)

where $t$ denotes time and $x$ the spatial direction, while $\vec{x}=(t,x)$.
In QCD $J^{\mu}$ would correspond to the electromagnetic current for the
appropriate combination of quark flavors; in the single flavor Schwinger model
it is simply given by

$\displaystyle J^{\mu}(t,x)=\bar{\psi}(t,x)\gamma^{\mu}\psi(t,x)\,,$ (2)

where $\psi$ is a continuum two-component fermionic field, and $\mu=0,1$. The
HVP tensor is then defined in momentum space as

$\displaystyle\Pi^{\mu\nu}(Q)$ $\displaystyle=\int d^{2}\vec{x}\int
d^{2}\vec{x}_{0}e^{i\vec{Q}\cdot(\vec{x}-\vec{x}_{0})}\langle
J^{\mu}(t,x)J^{\nu}(t_{0},x_{0})\rangle$
$\displaystyle=(Q^{\mu}Q^{\nu}-\delta^{\mu\nu}Q^{2})\Pi(Q^{2})\,,$ (3)

where we used that the HVP only depends on $\vec{x}-\vec{x}_{0}$.111For
further discussion on the hadronic tensor and related quantities in the
quantum computing context see e.g. [17, 18, 19, 20]. Similarly, the HLBL
contribution is related to four point functions of the electromagnetic
currents:

$\displaystyle\mathcal{C}^{\mu\alpha\beta\nu}=\langle
TJ^{\mu}(\vec{x})J^{\alpha}(\vec{w})J^{\beta}(\vec{y})J^{\nu}(\vec{x}_{0})\rangle\,,$
(4)

corresponding to the simplest realization of Fig. 1, when
$\hat{Q}^{\alpha\beta}=J^{\alpha}J^{\beta}$.

In this paper, we evaluate Eqs. (I) and (4) from real-time lattice
simulations. As alluded before, the computations are performed using
(classical) tensor network methods, which do not suffer from sign problems
related to real-time evolution. More, they offer a first step towards the
implementation in quantum computers, where larger scale and higher dimensional
quantum simulations should eventually become feasible in the near future. To
that end, we first review some basic properties of the Schwinger model in
section II, and then present numerical results from tensor network simulations
for the extraction of the HVP and HLBL tensors in sections III, IV. We further
illustrate our approach by providing a calculation using a QC emulator in a
small computational lattice. We conclude in section V by summarizing our
results and discussing future extensions of this work.

##  II The Schwinger model

The Schwinger model corresponds to the theory of quantum electrodynamics in
1+1-d [8]. In the temporal gauge, i.e. $A^{0}=0$, the Hamiltonian of this
model is given by

$\displaystyle H$ $\displaystyle=\int dx\,\frac{1}{2}E^{2}(x)$
$\displaystyle+\bar{\psi}(x)(-i\gamma^{1}\partial_{1}+g\gamma^{1}A_{1}(x)+m)\psi(x)\,,$
(5)

where $\psi$ denotes the single flavor two component spinor field with mass
$m$, while $A^{\mu}$ is the gauge field, of which only $A^{1}$ appears
explicitly in $H$. The electric field, $E=-\partial_{t}A^{1}$, carries all the
energy associated to the gauge degrees of freedom. Note that in two dimensions
there is no magnetic field, and, in fact, the gauge field is not dynamical and
acts as a potential. To make explicit the static character of the gauge field,
one can integrate out the remaining field component using Gauss’s law:
$\partial_{x}E=g\psi^{\dagger}\psi$.

The theory can be discretized on a lattice, following the Kogut-Susskind
prescription for fermions [21, 22]. In practice, this is done by introducing
the one component spinor $\chi_{n}$, which is related to the continuum one by

$\displaystyle\psi(\tilde{n})\to\frac{1}{\sqrt{a}}\begin{pmatrix}\chi_{2\tilde{n}}\\\
\chi_{2\tilde{n}-1}\end{pmatrix}\,.$ (6)

with $a$ the staggered lattice spacing.222We denote the naive lattice spacing
$\tilde{a}=2a$. Here we use $\tilde{n}$ as the index on the physical lattice;
we shall use $n$ as the site index on the staggered lattice. For each
$\tilde{n}$ there are two values of $n$. The gauge field is introduced in
terms of the reduced electric field $L\equiv E/g$ and can be integrated out by
using the discrete Gauss’ law generator

$\displaystyle G_{n}$ $\displaystyle=\delta
L_{n}-\Bigg{(}\chi_{n}^{\dagger}\chi_{n}-\frac{1}{2}(1-(-1)^{n})\Bigg{)}\,,$
(7)

where $\delta L_{n}\equiv L_{n}-L_{n-1}$. In the charge zero sector, one has
that any physical state satisfies $G|\psi\rangle=0$.

The single component spinors can be mapped to spin operators via a Jordan-
Wigner transform (JWt)

$\displaystyle\chi_{n}=S_{n}\sigma^{-}_{n}\,,\quad\chi^{\dagger}_{n}=S^{\dagger}_{n}\sigma^{+}_{n}\,;$
(8)

where the string operator reads $S_{n}\equiv\prod_{k<n}[(-i)\sigma^{z}_{k}]$,
and $\sigma^{x,y,z,+,-}_{n}$ denote the standard Pauli matrices acting on the
site $n$. We define $\sigma^{\pm}=(\sigma^{x}\pm i\sigma^{y})/2$. Using the
lattice discretization described and employing the JWt, one finds that the
continuum Hamiltonian in Eq. (II) maps to the spin Hamiltonian

$\displaystyle H$
$\displaystyle=\frac{g^{2}a}{2}\sum_{n=1}^{N-1}\left[\frac{1}{2}\sum_{k=1}^{n}(\sigma^{z}_{k}+(-1)^{k})\right]^{2}+\sum_{n=1}^{N}m(-1)^{n}\frac{\sigma^{z}_{n}}{2}$
$\displaystyle+\frac{1}{2a}\sum_{n=1}^{N-1}\sigma^{+}_{n}\sigma^{-}_{n+1}+\sigma^{+}_{n+1}\sigma^{-}_{n}\,,$
(9)

where we imposed open boundary conditions, such that for the $N$th lattice
point $L_{N}=\chi_{N+1}=0$. For the $\gamma$ matrices we used the basis
$\gamma^{0}=\sigma^{z}$, $\gamma^{1}=-i\sigma^{y}$, and
$\gamma^{5}=-\sigma^{x}$.

Finally, having provided the lattice version of the Schwinger model, we turn
our attention to the form of the current operators on the lattice. In what
follows, we shall only need the temporal component of the electric charge
density, $J^{0}$, which reads

$\displaystyle\psi^{\dagger}(\tilde{n})\psi(\tilde{n})$
$\displaystyle=\frac{1}{a}\left(\chi_{2\tilde{n}}^{\dagger}\chi_{2\tilde{n}}+\chi_{2\tilde{n}-1}^{\dagger}\chi_{2\tilde{n}-1}\right)$
$\displaystyle=\frac{1}{a}(\sigma^{+}_{2\tilde{n}}\sigma^{-}_{2\tilde{n}}+\sigma^{+}_{2\tilde{n}-1}\sigma^{-}_{2\tilde{n}-1})\,.$
(10)

##  III Hadronic vacuum polarization

Having introduced the Schwinger model, we now discuss how to extract the HVP
in this theory. In general, one can rewrite the position space form of the
current-current correlator as

$\displaystyle\Pi^{\mu\nu}=\langle
J^{\mu}(t,x)J^{\nu}(0,x_{0})\rangle\equiv\langle\psi^{\mu}_{l}|\psi^{\nu}_{r}\rangle\,,$
(11)

with

$\displaystyle|\psi_{l}^{\mu}\rangle$
$\displaystyle=j^{\mu}(\tilde{n})e^{-iHt}|\Omega\rangle\,,$
$\displaystyle|\psi_{r}^{\nu}\rangle$
$\displaystyle=e^{-iHt}j^{\nu}(\tilde{n}_{0})|\Omega\rangle\,,$ (12)

where $j^{\mu}$ denotes the discretized version of the current operator and
$|\Omega\rangle$ is the ground state of the Schwinger model at finite $g$ and
$m$. To perform the Fourier transform necessary to obtain Eq. (I), we use

$\displaystyle\Pi^{\mu\nu}(Q^{0},Q^{1})=\int
dx\,dt\,e^{iQ^{0}t}e^{-iQ^{1}x}\Pi^{\mu\nu}(t,x)$
$\displaystyle=\sum_{j=0}^{N_{x}-1}\sum_{k=0}^{N_{t}-1}\,\Delta t\,\Delta
x\,e^{iQ^{0}k\Delta t}e^{-iQ^{1}j\Delta x}\Pi^{\mu\nu}(k\Delta t,j\Delta
x)\,,$ (13)

where $N_{t}$ and $N_{x}$ denote the number of points taken in the time and
spatial directions, respectively. Here $\Delta t$ and $\Delta x$ correspond to
the spacing in the time and spatial directions; the respective momentum modes
read $Q^{0}=0,\frac{2\pi}{N_{t}\Delta t},\cdots\frac{2\pi}{N_{t}\Delta
t}(N_{t}-1)$ and $Q^{1}=0,\frac{2\pi}{N_{x}\Delta
x},\cdots\frac{2\pi}{N_{x}\Delta x}(N_{x}-1)$.

We start by discussing the form of Eq. (III) in the exactly solvable limits of
the Schwinger model. The first limit corresponds to the case where $m=0$ and
the theory can be shown to become that of free massive bosons, with a mass
$m_{b}=g/\sqrt{\pi}$ [10]. In this case one has that [23]

$\displaystyle\Pi(Q^{2})\Big{|}_{m=0}=\frac{1}{\pi\left(Q^{2}+m_{b}^{2}\right)}\,.$
(14)

The HVP is regular for spacelike $Q$, but has a pole at timelike separations
at the value of the boson mass. Another interesting limit corresponds to the
case where one decouples the fermion sector from the gauge field, which can be
achieved in the $g\to 0$ limit, or by treating $A^{\mu}$ as an external field.
In this case, one has that [23]

$\displaystyle\Pi(Q^{2})\Big{|}_{g=0}$ $\displaystyle=\frac{1}{\pi
Q^{2}}+\frac{4m^{2}\arctan\left(\frac{\sqrt{Q^{2}}}{\sqrt{-4m^{2}-Q^{2}}}\right)}{\pi\sqrt{Q^{6}}\sqrt{-4m^{2}-Q^{2}}}$
$\displaystyle\stackrel{{\scriptstyle m\to 0}}{{\to}}\frac{1}{\pi Q^{2}}\,.$
(15)

In this case there is a pole associated to the massless limit of the theory,
while when $m\neq 0$, there is a branch cut along the timelike axis related to
the particle production threshold.

Finally, to go beyond these two cases, when $m$ and $g$ are taken finite, one
has to extract $\Pi^{\mu\nu}$ from the lattice, which we follow to do using
tensor network methods. In the what follows we use $a\cdot g=0.50$, $a\cdot
m=0.25$, and we consider the lattice improved action by shifting the mass term
as $m_{\rm lat}=m_{\rm bare}-g^{2}a/8$ [24]. We will focus on $\Pi^{00}$; the
extraction of the two other components poses no extra difficulties. In
particular, note that $\Pi(Q^{2})=\Pi^{00}(Q^{2})/Q_{1}^{2}$ for $Q$ not
timelike.

Figure 2: HVP extracted from the tensor network simulations using MPS. The
quantities are plotted in terms of the lattice spacing $a$.

In Fig. 2 we show the HVP extracted from a tensor network simulation. We used
the tensor network package ITensor [25] in Julia, adopting a matrix product
state (MPS) form to represent the quantum states. The time evolution of the
states in Eq. (III) is achieved using the time-dependent variational principle
(TDVP) algorithm [26, 27], while the preparation of the initial state of the
system is performed using the density matrix renormalization group (DMRG)
algorithm [28, 29]. We utilized a lattice with $N=120$ staggered sites, and we
extract the observables in the center of the lattice, minimizing edge effects,
which are suppressed by the finite fermion mass. We use six values of spatial
separation for $x-x_{0}=0,1,\cdots 5$ in units of the lattice spacing, and
evolve the states up sufficiently large times, $0<t-t_{0}<5$ in lattice
spacing units, reducing the time steps to ensure that the trotterization of
the evolution operator does not significantly affect the final distribution.
The result shown exhibits the expected feature for the parameters used, i.e. a
decay in both the spacelike and timelike directions, with a peak around the
origin. Nonetheless, we note that the obtained curve’s quality depends greatly
on the numerical Fourier transform one has to perform. Nonetheless, this
operation is purely classical, and even when using a QC, it would be realized
the same way as illustrated here. Thus, the only requirement to obtain a
better numerical result for $\Pi^{00}$, lies on increasing (decreasing) the
lattice size (spacing) and collecting more data points.

In Fig. 3 we show the real-time quantum simulation of $\Pi^{00}$ (in position
space) for $\Delta x=x-x_{0}=2$ in lattice units, in a smaller system with
$N=20$ qubits.

Figure 3: Time evolution of $\text{Re}(\Pi^{00})$, extracted from classical
emulator of a digital quantum computing, for $\Delta x=x-x_{0}=2$ in lattice
units; the remaining parameters are the same as the ones used in Fig. 2.

We again recall that the transformation to momentum space can be done the same
way as shown for the above MPS example. However, one important difference
compared to tensor network methods relates to the extraction of the
expectation value of hermitian operators. In the example shown in Fig. 3 this
is done by using the simple interference measurement protocol illustrated in
Fig. 4. Note that by measuring the ancilla qubit of
$\frac{Z_{m}U}{\sqrt{2}}\ket{\Omega}\ket{0}_{c}+\frac{UZ_{n}}{\sqrt{2}}\ket{\Omega}\ket{1}_{c}$
in the $X$-basis, one can compute
$\text{Re}\bra{\Omega}Z_{m}(t)Z_{n}\ket{\Omega}$, where $\ket{\Omega}$ is the
ground state of the Schwinger model, $Z_{m}(t)=U^{\dagger}Z_{m}U$, $Z_{m}$
would correspond to the discretized version of the current operator, and
$U=e^{-itH}$. A Qiskit code to compute this function is provided in [30].

Figure 4: Measurement protocol used to extract $\Pi^{00}$. Here, $X$ and $H$
operators denote the $\sigma^{x}$ and Hadamard gates, respectively.

##  IV HLBL and hadronic structure of the photon

Having computed the HVP, we now consider the HLBL contribution. This involves
the (Fourier transform of) expectation values of the form shown in Eq. (4),
see also Fig. 1, which can be extracted in a similar fashion to the HVP,
despite the more complex numerical calculation. We note that the HLBL
contribution is closely related to the correlators determining the photon’s
hadronic structure [31, 32, 33]. In this case, and in QCD, one is interested
in computing the expectation value of a quark/gluon local operator, in between
virtual photon states, see Fig 1. For example, at leading order in
$\alpha_{em}$ and in 3+1-d, this results in

$\displaystyle\mathcal{C}^{\mu\mu_{1}\cdots\mu_{n}\nu}(\vec{x}-\vec{x}_{0})\varepsilon^{\mu}(\lambda)\varepsilon^{*,\mu}(\lambda)\,,$
(16)

where $\mathcal{C}^{\mu\mu_{1}\cdots\mu_{n}\nu}=\langle
J^{\mu}(\vec{x})\hat{O}^{\mu_{1}\cdots\mu_{n}}(0)J^{\nu}(\vec{x}_{0})\rangle$
is illustrated in Fig. 1, with $\hat{O}$ denoting the local operator being
measured, and $\varepsilon^{\mu}$ is the photon polarization vector, which has
$d-1$ physical polarizations. For $d=1$ there are no physical polarizations
(the electric field is not propagating), and the contraction in Eq. (16)
vanishes exactly. Nevertheless, it is still possible to compute the correlator
$\mathcal{C}$; to that end we consider the simplest case where the local
operator is given by the product of two electromagnetic currents:

$\displaystyle\mathcal{C}^{\mu\mu_{1}\mu_{2}\nu}_{4}(\vec{x}-\vec{x}_{0})=\langle
TJ^{\mu}(\vec{x})J^{\mu_{1}}(0)J^{\mu_{2}}(0)J^{\nu}(\vec{x}_{0})\rangle\,,$
(17)

corresponding to Eq. (4), with $\vec{w}=\vec{y}=0$.

In Fig. 5 we show the evaluation of Eq. (17) using a MPS tensor network, for
$\Delta x=x-x_{0}=2,4$ in units of $\tilde{a}$ as a function of time
($t_{0}=0$). We use a lattice with $N=60$ sites in a set up analog to the HVP
case, and the intermediate current operators are evaluated in the middle of
the lattice. Numerically, we observe that the extraction of this expectation
value is considerably more complex than the HVP case, with the bond dimensions
of the intermediate links in the MPS growing much faster with time. This
limits the amount of data that can be extracted efficiently in our simple
evaluation. We, therefore, refrain from performing the final Fourier transform
to momentum space.

Figure 5: Four point current correlator extracted from a tensor network
simulation, following Eq. (17) with $\mu=\nu=\mu_{1}=\mu_{2}=0$. Here $\Delta
x$ denotes the spatial separation, while the horizontal axis denotes the time
evolution of the expectation value.

##  V Conclusion and discussion

In this work we have explored the potential use of quantum computing methods
and technologies to extract the HVP from timelike momenta transfers. These
calculations can be extended to access the HBLB as well as the hadronic
structure of the photon. While these computations are related to the structure
of hadrons [17, 34, 35], we note an important difference between studying the
photon structure and the partonic distributions inside hadrons: the latter
will require preparing a hadronic target state, which is usually hard to
obtain through quantum computing methods; the former has as external states
virtual photons, which can be easily created from the vacuum by the insertion
of electromagnetic currents.

In order to illustrate the main points of our discussion, we have provided a
small system size tensor network calculation of the HVP and HLBL in the
Schwinger model. The major difficulty in the extraction of these quantities
lies in the need to Fourier transform to momentum space, requiring a finer
spatial and time grid than the ones employed here. We hope to address these
numerical issues in the future, with more precise numerical calculation. We
further note that when translating these computations to a real quantum
computer, one may face further difficulties, such as the extraction of the
expectation value of (hermitian) observables [17], which are not present using
tensor networks; this is one of the important differences between these two
approaches. We further illustrated this point, by extracting the position
space HVP using a quantum simulator.

The methods employed here surpass some of the issues faced using either
dispersive relations or euclidean lattice methods, but they require a
significant progress in a formulation of QCD which could be implemented in
real quantum computers. We hope coming developments in the field will allow
for the questions and methods posed here to be tested in a context closer to
QCD or in other field theories, see e.g. [36, 37, 38, 39, 40, 41, 42, 43, 44,
45].

###### Acknowledgements.

We thank A. Florio, R. Pisarski, A. V. Sadofyev, and R. Szafron for helpful
discussions related to this work. This material is based upon work supported
by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics
and National Quantum Information Science Research Centers, Co-design Center
for Quantum Advantage (C2QA) under contract number DE-SC0012704. The work of
J.R. was supported by the U.S. Department of Energy, Office of Science, Office
of Workforce Development for Teachers and Scientists (WDTS) under the Science
Undergraduate Laboratory Internships Program (SULI) and the Brookhaven
National Laboratory (BNL), Physics Dept. under the BNL Supplemental
Undergraduate Research Program (SURP).

## References

  * [1] T. Aoyama, et al., The anomalous magnetic moment of the muon in the Standard Model, Phys. Rept. 887 (2020) 1–166.  arXiv:2006.04822, doi:10.1016/j.physrep.2020.07.006.
  * [2] F. Jegerlehner, A. Nyffeler, The Muon g-2, Phys. Rept. 477 (2009) 1–110.  arXiv:0902.3360, doi:10.1016/j.physrep.2009.04.003.
  * [3] N. Asmussen, A. Gérardin, A. Nyffeler, H. B. Meyer, Hadronic light-by-light scattering in the anomalous magnetic moment of the muon, SciPost Phys. Proc. 1 (2019) 031.  arXiv:1811.08320, doi:10.21468/SciPostPhysProc.1.031.
  * [4] F. Jegerlehner, Muon g – 2 theory: The hadronic part, EPJ Web Conf. 166 (2018) 00022.  arXiv:1705.00263, doi:10.1051/epjconf/201816600022.
  * [5] C. W. Bauer, et al., Quantum Simulation for High-Energy Physics, PRX Quantum 4 (2) (2023) 027001.  arXiv:2204.03381, doi:10.1103/PRXQuantum.4.027001.
  * [6] A. Di Meglio, et al., Quantum Computing for High-Energy Physics: State of the Art and Challenges. Summary of the QC4HEP Working Group (7 2023).  arXiv:2307.03236.
  * [7] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2 (2018) 79.  arXiv:1801.00862, doi:10.22331/q-2018-08-06-79.
  * [8] J. S. Schwinger, On gauge invariance and vacuum polarization, Phys. Rev. 82 (1951) 664–679.  doi:10.1103/PhysRev.82.664.
  * [9] A. Casher, J. B. Kogut, L. Susskind, Vacuum polarization and the absence of free quarks, Phys. Rev. D 10 (1974) 732–745.  doi:10.1103/PhysRevD.10.732.
  * [10] S. R. Coleman, More About the Massive Schwinger Model, Annals Phys. 101 (1976) 239.  doi:10.1016/0003-4916(76)90280-3.
  * [11] T. Banks, L. Susskind, J. B. Kogut, Strong Coupling Calculations of Lattice Gauge Theories: (1+1)-Dimensional Exercises, Phys. Rev. D 13 (1976) 1043.  doi:10.1103/PhysRevD.13.1043.
  * [12] C. J. Hamer, W.-h. Zheng, J. Oitmaa, Series expansions for the massive Schwinger model in Hamiltonian lattice theory, Phys. Rev. D 56 (1997) 55–67.  arXiv:hep-lat/9701015, doi:10.1103/PhysRevD.56.55.
  * [13] H. J. Rothe, Lattice Gauge Theories : An Introduction (Fourth Edition), Vol. 43, World Scientific Publishing Company, 2012.  doi:10.1142/8229.
  * [14] M. C. Bañuls, et al., Simulating Lattice Gauge Theories within Quantum Technologies, Eur. Phys. J. D 74 (8) (2020) 165.  arXiv:1911.00003, doi:10.1140/epjd/e2020-100571-8.
  * [15] R. C. Farrell, M. Illa, A. N. Ciavarella, M. J. Savage, Quantum Simulations of Hadron Dynamics in the Schwinger Model using 112 Qubits (1 2024).  arXiv:2401.08044.
  * [16] R. C. Farrell, M. Illa, A. N. Ciavarella, M. J. Savage, Scalable Circuits for Preparing Ground States on Digital Quantum Computers: The Schwinger Model Vacuum on 100 Qubits (8 2023).  arXiv:2308.04481.
  * [17] H. Lamm, S. Lawrence, Y. Yamauchi, Parton physics on a quantum computer, Phys. Rev. Res. 2 (1) (2020) 013272.  arXiv:1908.10439, doi:10.1103/PhysRevResearch.2.013272.
  * [18] N. Mueller, A. Tarasov, R. Venugopalan, Deeply inelastic scattering structure functions on a hybrid quantum computer, Phys. Rev. D 102 (1) (2020) 016007.  arXiv:1908.07051, doi:10.1103/PhysRevD.102.016007.
  * [19] T. Li, X. Guo, W. K. Lai, X. Liu, E. Wang, H. Xing, D.-B. Zhang, S.-L. Zhu, Partonic collinear structure by quantum computing, Phys. Rev. D 105 (11) (2022) L111502.  arXiv:2106.03865, doi:10.1103/PhysRevD.105.L111502.
  * [20] W. Qian, R. Basili, S. Pal, G. Luecke, J. P. Vary, Solving hadron structures using the basis light-front quantization approach on quantum computers, Phys. Rev. Res. 4 (4) (2022) 043193.  arXiv:2112.01927, doi:10.1103/PhysRevResearch.4.043193.
  * [21] L. Susskind, Lattice fermions, Phys. Rev. D 16 (1977) 3031–3039.  doi:10.1103/PhysRevD.16.3031.
URL https://link.aps.org/doi/10.1103/PhysRevD.16.3031

  * [22] J. B. Kogut, L. Susskind, Hamiltonian Formulation of Wilson’s Lattice Gauge Theories, Phys. Rev. D 11 (1975) 395–408.  doi:10.1103/PhysRevD.11.395.
  * [23] D. N. Blaschke, R. Carballo-Rubio, E. Mottola, Fermion Pairing and the Scalar Boson of the 2D Conformal Anomaly, JHEP 12 (2014) 153.  arXiv:1407.8523, doi:10.1007/JHEP12(2014)153.
  * [24] R. Dempsey, I. R. Klebanov, S. S. Pufu, B. Zan, Discrete chiral symmetry and mass shift in the lattice Hamiltonian approach to the Schwinger model, Phys. Rev. Res. 4 (4) (2022) 043133.  arXiv:2206.05308, doi:10.1103/PhysRevResearch.4.043133.
  * [25] M. Fishman, S. R. White, E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases (2022) 4doi:10.21468/SciPostPhysCodeb.4.
URL https://scipost.org/10.21468/SciPostPhysCodeb.4

  * [26] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pizorn, H. Verschelde, F. Verstraete, Time-Dependent Variational Principle for Quantum Lattices, Phys. Rev. Lett. 107 (2011) 070601.  arXiv:1103.0936, doi:10.1103/PhysRevLett.107.070601.
  * [27] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94 (2016) 165116.  doi:10.1103/PhysRevB.94.165116.
URL https://link.aps.org/doi/10.1103/PhysRevB.94.165116

  * [28] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69 (1992) 2863–2866.  doi:10.1103/PhysRevLett.69.2863.
URL https://link.aps.org/doi/10.1103/PhysRevLett.69.2863

  * [29] S. R. White, Density-matrix algorithms for quantum renormalization groups, Phys. Rev. B 48 (1993) 10345–10356.  doi:10.1103/PhysRevB.48.10345.
URL https://link.aps.org/doi/10.1103/PhysRevB.48.10345

  * [30] K. Ikeda, Qiskit implementation of two-point correlation function.
URL https://github.com/IKEDAKAZUKI/Qiskit-
Tutorial/blob/main/Quantum_computation_of_n_point_correpation_functions.ipynb

  * [31] B. Jager, M. Stratmann, W. Vogelsang, Longitudinally polarized photoproduction of inclusive hadrons at fixed-target experiments, Eur. Phys. J. C 44 (2005) 533–543.  arXiv:hep-ph/0505157, doi:10.1140/epjc/s2005-02380-0.
  * [32] X.-d. Ji, C.-w. Jung, Photon structure functions from quenched lattice QCD, Phys. Rev. D 64 (2001) 034506.  arXiv:hep-lat/0103007, doi:10.1103/PhysRevD.64.034506.
  * [33] X.-d. Ji, C.-w. Jung, Studying hadronic structure of the photon in lattice QCD, Phys. Rev. Lett. 86 (2001) 208.  arXiv:hep-lat/0101014, doi:10.1103/PhysRevLett.86.208.
  * [34] S. Grieninger, K. Ikeda, I. Zahed, Quasi-parton distributions in massive QED2: Towards quantum computation (4 2024).  arXiv:2404.05112.
  * [35] S. Grieninger, I. Zahed, Quasi-fragmentation functions in the massive Schwinger model (6 2024).  arXiv:2406.01891.
  * [36] R. Dempsey, I. R. Klebanov, S. S. Pufu, B. T. Søgaard, Lattice Hamiltonian for Adjoint QCD2 (11 2023).  arXiv:2311.09334.
  * [37] A. M. Czajka, Z.-B. Kang, Y. Tee, F. Zhao, Studying chirality imbalance with quantum algorithms (10 2022).  arXiv:2210.03062.
  * [38] S. Grieninger, K. Ikeda, D. E. Kharzeev, I. Zahed, Entanglement in massive schwinger model at finite temperature and density, Phys. Rev. D 109 (2024) 016023.  doi:10.1103/PhysRevD.109.016023.
URL https://link.aps.org/doi/10.1103/PhysRevD.109.016023

  * [39] J. Barata, W. Gong, R. Venugopalan, Realtime dynamics of hyperon spin correlations from string fragmentation in a deformed four-flavor Schwinger model (8 2023).  arXiv:2308.13596.
  * [40] N. Klco, E. F. Dumitrescu, A. J. McCaskey, T. D. Morris, R. C. Pooser, M. Sanz, E. Solano, P. Lougovski, M. J. Savage, Quantum-classical computation of Schwinger model dynamics using quantum computers, Phys. Rev. A 98 (3) (2018) 032331.  arXiv:1803.03326, doi:10.1103/PhysRevA.98.032331.
  * [41] N. Butt, S. Catterall, Y. Meurice, R. Sakai, J. Unmuth-Yockey, Tensor network formulation of the massless Schwinger model with staggered fermions, Phys. Rev. D 101 (9) (2020) 094509.  arXiv:1911.01285, doi:10.1103/PhysRevD.101.094509.
  * [42] G. Magnifico, M. Dalmonte, P. Facchi, S. Pascazio, F. V. Pepe, E. Ercolessi, Real Time Dynamics and Confinement in the $\mathbb{Z}_{n}$ Schwinger-Weyl lattice model for 1+1 QED, Quantum 4 (2020) 281.  arXiv:1909.04821, doi:10.22331/q-2020-06-15-281.
  * [43] D. E. Kharzeev, Y. Kikuchi, Real-time chiral dynamics from a digital quantum simulation, Phys. Rev. Res. 2 (2) (2020) 023342.  arXiv:2001.00698, doi:10.1103/PhysRevResearch.2.023342.
  * [44] R. Belyansky, S. Whitsitt, N. Mueller, A. Fahimniya, E. R. Bennewitz, Z. Davoudi, A. V. Gorshkov, High-energy collision of quarks and mesons in the schwinger model: From tensor networks to circuit qed, Phys. Rev. Lett. 132 (2024) 091903.  doi:10.1103/PhysRevLett.132.091903.
URL https://link.aps.org/doi/10.1103/PhysRevLett.132.091903

  * [45] J. W. Pedersen, E. Itou, R.-Y. Sun, S. Yunoki, Quantum Simulation of Finite Temperature Schwinger Model via Quantum Imaginary Time Evolution, PoS LATTICE2023 (2024) 220.  arXiv:2311.11616, doi:10.22323/1.453.0220.