Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
69,280 characters · 16 sections · 63 citation commands
Conformal Prediction Bands for Two-Dimensional Functional Time Series
Data observed on a two-dimensional domain arise naturally across several disciplines, motivating an increasing demand for dedicated analysis techniques. Functional data analysis (FDA) (ramsay2005functional) is naturally apt to represent and model this kind of data, as it allows preserving their continuous nature, and provides a rigorous mathematical framework. Among the others, Zhou analyzed temperature surfaces, presenting two approaches for Functional Principal Component Analysis (FPCA) of functions defined on a non-rectangular domain, Munoz focuses on image processing using FDA, while a novel regularization technique for Gaussian random fields on a rectangular domain has been proposed by Raket and applied to 2D electrophoresis images. Another bivariate smoothing approach in a penalized regression framework has been introduced by Ivanescu, allowing for the estimation of functional parameters of two-dimensional functional data. As shown by Gervini, even mortality rates can be interpreted as two-dimensional functional data.
Whereas in all the reviewed works functions are assumed to be realization of iid or at least exchangeable random objects, to the best of our knowledge there is no literature focusing on forecasting time-dependent two-dimensional functional data. In this work, we focus on time series of surfaces, representing them as two-dimensional Functional Time Series (FTS).
A two-dimensional Functional Time Series is an ordered sequence $Y_1, \dots, Y_T$ of random variables with values in a functional Hilbert space $\mathbb{H}$, characterized by some sort of temporal dependency. More formally, we consider a probability space $(\Omega, \mathcal{F},\mathbb{P})$, and define a random function at time $t$ as $Y_t: \Omega \rightarrow \mathbb{H}$, measurable with respect to the Borel $\sigma$-algebra $\mathcal{B}(\mathbb{H})$. In the rest of the article we consider functions belonging to $\mathbb{H} = \mathcal{L}^2([c,d] \times [e,f])$, with $c,d,e,f \in \mathbb{R}$, $c<d$, $e<f$ We stress the fact that, from a theoretical point of view, our methodology can be applied to functions defined on a generic subset of $\mathbb{R}^2$, however, for simplicity and without loss of generality, we will only consider rectangular domains.
Given a realization of the stochastic process $\{Y_t\}_{t= \, \dots, N}$, we aim to forecast the next surface and quantify the uncertainty around the predicted function. Whereas uncertainty quantification in the context of univariate FTS forecasting has received great attention in the statistical community in recent decades, no attempts have been made to extend them to functions defined on a bidimensional domain.
Most of the research tackling univariate FTS forecasting has focused on adaption of the Bootstrap to the functional setting (see e.g. Hyndman_Shang, Rossini and hernandez2021simultaneous). However, Bootstrap is a very computationally intensive procedure, especially in the infinite-dimensional context of functional data. In this work, we instead focus on Conformal Prediction (CP), a versatile nonparametric approach to prediction. The first appearance of such technique dates back to gammerman, and it has been later presented in greater detail in Vovk2005AlgorithmicLI and Balasubramanian. An extensive and unified review of the theory of Conformal inference can be found in Zeni2020ConformalPA. The attractiveness of Conformal Prediction relies on its great versatility, which allows to couple it with any predictive algorithm, in order to obtain distribution free prediction sets. Throughout this work, we resort to Inductive Conformal Prediction, also known as Split Conformal Prediction (Papadopoulos2002InductiveCM). Such modification of the original Transductive Conformal method is not only computationally efficient, but also necessary in high-dimensional frameworks like the functional data one. It should be noted that the main drawback of the Full Conformal approach is the need of retraining the prediction algorithm for every possible candidate realization $y$. In practice, in multivariate problems, where $y$ lies in $\mathbb{R}^p$, one runs the above routine for several candidates $y$ over a $p$-dimensional regular grid. While such approach is prohibitive for high-dimensional spaces, since computational times grow exponentially with $p$, it becomes unfeasible in a functional setting, in which $y$ lies in an infinite-dimensional space. Employing Split Conformal inference along with nonconformity scores tailored to functional data (diquigiovanni2021importance) allows to obtain prediction sets in closed form.
Whereas CP theory has been originally developed under the assumption of exchangeability, chernozhukov2018exact reframed the CP framework in the context of randomization inference, proving approximate validity of the resulting prediction sets under weak assumptions on the conformity scores and on the ergodicity of the time series. Later, diquigiovanni2021_FTS adapted such methodology to allow for Functional Time Series in a Split Conformal setting. We extend such method to two-dimensional functional data.
Since we want to quantify prediction uncertainty, we also need to provide a forecasting technique for 2D Functional Time Series. We start by taking into account the literature on unidimensional FTS forecasting (bosq, antoniadis, horvath2012inference, Aue, jiao), extending the theory of Functional Autoregressive Processes (FAR) to the bivariate setting, proposing estimation techniques for the FAR(1), and comparing them in an order to assess how forecasting performances influence the amplitude of prediction bands.
The rest of this paper is as follows: we illustrate Conformal Prediction for two-dimensional functional data in Section (ref), providing theoretical guarantees of the resulting prediction bands. In Section (ref), we introduce forecasting algorithms for two-dimensional Functional Time Series, proposing an extension of the FAR(1) for two-dimensional functional data. Such forecasting algorithms are then compared by means of the resulting prediction bands in a simulation study in Section (ref). Finally, in Section (ref) we employ the developed techniques to obtain forecasts and prediction bands of real data, predicting day by day the Black Sea level. Section (ref) concludes.
Consider a time series $Z_1, \dots, Z_T$ of regression pairs $Z_t = (X_t, Y_t)$, with $t=1, \dots, T$. Let $Y_t$ be a random variable with values in $\mathbb{H}$, while $X_t$ is a set of random covariates at time $t$ belonging to a measurable space. Notice that $X_t$ is a generic set of regressors, which may contain both exogenous and endogenous variables. Later in the manuscript, we will consider $X_t$ to contain only the lagged version of the function $Y_t$, namely $Y_{t-1}$. Given a significance level $\alpha \in [0,1]$, we aim to design a procedure that outputs a prediction set $\mathcal{C}_{T,1-\alpha}(X_{T+1})$ for $Y_{T+1}$ based on $Z_1, \dots, Z_{T}$ and $X_{T+1}$, with unconditional coverage probability close to $1-\alpha$. More formally, we define $\mathcal{C}_{T,1-\alpha}(X_{T+1})$ to be a valid prediction set if:
We would like to construct a specific type of prediction sets, commonly known as prediction bands, formally defined as:
where $B_n(u,v) \subseteq \mathbb{R}$ is an interval for each $(u,v) \in [c,d] \times [e,f]$. The convenience of such type of prediction sets in applications is extensively motivated in the literature (see e.g. Pintado, Lei1 and diquigiovanni2021importance), since a prediction set of this kind can be visualized easily, a property that is instead not guaranteed if the prediction region is a generic subset of $\mathbb{H}$.
Let $z_1, \dots, z_T$ be realizations of $Z_1, \dots, Z_T$. As the name suggests, Split Conformal inference is based on a random split of data into two disjoint sets: let $\mathcal{I}_1$, $\mathcal{I}_2$ be a random partition of $\{1, \dots, T\}$, such that $|\mathcal{I}_1|=m$, $|\mathcal{I}_2|=l$, $m,l \in \mathbb{N}$ $m,l > 0$, $m+l=T$. Historical observations $z_1, \dots, z_T$ are divided into a training set $\{z_h,\, h \in \mathcal{I}_1\}$, used for model estimation, and a calibration set $\{z_h,\, h \in \mathcal{I}_2\}$, used in an out-of-sample context to measure the nonconformity of a new candidate function. The choice of the split ratio and the type of split is non-trivial and has motivated discussion in the statistical community. We fix the training-calibration ratio equal to 1 and perform a random split, and refer to (ref) for a more extensive discussion on the topic.
We then introduce a nonconformity measure $\mathcal{A}$, which is a measurable function with values in $\mathbb{R} \cup \{+\infty\}$. The role of $\mathcal{A}(\{z_h, \, h \in \mathcal{I}_1\}, z)$ is to quantify the nonconformity of a new datum $z$ with respect to the training set $\{z_h, \, h \in \mathcal{I}_1\}$. The choice of the nonconformity measure is crucial to find prediction bands (ref) in closed form. We employ the following nonconformity score, introduced by diquigiovanni2021importance, extended here to two-dimensional functional data:
where $z=(x_{T+1}, y)$, $g_{\mathcal{I}_1}$ is a point predictor built from the training set $\mathcal{I}_1$, and $s_{\mathcal{I}_1}$ is a modulation function, which is a positive function depending on $\mathcal{I}_1$ that allows for prediction bands with non-constant width along the domain. Section (ref) and (ref) discuss the estimation of $g_{\mathcal{I}_1}$. The functional standard deviation is employed as modulation function $s_{\mathcal{I}_1}$, allowing for wider bands in the parts of the domain where data show high variability and narrower and more informative prediction bands in those parts characterized by low variability. For an extensive discussion on the optimal choice of modulation function, we refer to diquigiovanni2021importance.
Consider now a candidate function $y \in \mathbb{H}$ and define the augmented dataset as $Z_{(y)} = \{Z_t\}_{t=1}^{T+1}$, where:
The key idea of the methodology proposed by chernozhukov2018exact and extended by diquigiovanni2021_FTS is to generate several randomized versions of $Z_{(y)}$ through a specifically tailored permutation scheme, and compute nonconformity scores on each of them. We then decide whether to include $y$ in the prediction region, by comparing the nonconformity score of $Z_{(y)}$ with that of its permuted replicas.
In order to obtain such replicas, we aim to define a family $\Pi$ of index permutations $\pi: \{1,\dots,T+1\} \rightarrow \{1,\dots,T+1\}$, that keeps unchanged the training set indices $\mathcal{I}_1$, and modifies only $\{\mathcal{I}_2, T+1\}$, namely the indices of the calibration set and the next time step. We first introduce a function $\lambda: \{\mathcal{I}_2, T+1\} \rightarrow \{1,\dots,l+1\}$ such that $\lambda(t)$ returns the $t$-th element of the ordered set $\{\mathcal{I}_2,T+1\}$. Fix now a positive integer $b \in \{1, \dots l+1\}$ such that $\frac{l+1}{b} \in \mathbb{N}$ and define a family $\tilde{\Pi}$ of index permutations that acts on the set $\{1,\dots,l+1\}$. Each $\tilde{\pi}_i \in \tilde{\Pi}$ is required to be a bijection $\tilde{\pi}_i : \{1,\dots,l+1\} \rightarrow \{1,\dots,l+1\}$, for $i = 1, \dots, \frac{l+1}{b}$. We consider the non-overlapping blocking permutation scheme proposed by chernozhukov2018exact, dividing data in blocks of size $b$, in such a way that each permutation is unique in $\tilde{\Pi}$:
By definition, we have that $|\tilde{\Pi}|= \frac{l+1}{b}$, and $\tilde{\Pi}$ forms an algebraic group, containing the identity transformation $\tilde{\pi}_1$. It is then straightforward to introduce the family $\Pi$ of index permutations acting on $\{1,\dots,T+1\}$. Each $\pi_i \in \Pi$, with $i=1,\dots,\frac{l+1}{b}$ is defined as:
Figure (ref) reports an example of the families of permutation $\Pi$ and $\tilde{\Pi}$. We refer to $Z^{\pi}_{(y)} = \{Z_{\pi(t)}\}_{t=1}^{T+1}$ as the randomized version of $Z_{(y)}=\{Z_t\}_{t=1}^{T+1}$, and define the randomization p-value as:
where nonconformity scores $S(Z_{(y)})$ and $S(Z^{\pi}_{(y)})$ are defined as:
The idea is to apply permutations, modifying the order of observations in the calibration set, while at the same time preserving the dependence between them, thanks to the block structure of $\Pi$. For each $\pi$, we compute the nonconformity score of $Z^{\pi}_{(y)}$. The p-value (ref) of a test candidate value $y$ is then determined as the proportion of randomized versions $Z^{\pi}_{(y)}$ with a higher or equal nonconformity score than the one of the original augmented dataset $Z_{(y)}$. Notice that $p(y)$ is a measure of the conformity of the candidate function $y$ with respect to the permutation family $\Pi$. It is then natural to include in the prediction set only functions $y$ with an “high" conformity level. Given a significance level $\alpha \in [b/(l+1),1]$, the prediction region is hence obtained by test inversion:
As a remark, it should be noted that if $\alpha \in (0,b/(l+1))$ the resulting prediction set coincides with the entire space $\mathbb{H}$ (diquigiovanni2021importance).
The advantage of using the Split Conformal method along with the conformity measure (ref) relies on the possibility to find the prediction set in closed form. By defining $k^s$ as the $\lceil (|\Pi|+1)(1-\alpha) \rceil$th smallest value of $\{S(Z_{(y)}^{\pi}), \pi \in \Pi \setminus \pi_1\}$, we derive:
The prediction band is therefore:
If regression pairs are exchangeable, the proposed method retains exact, model-free validity (chernozhukov2018exact, Theorem 1). When such assumption is not met, one can instead guarantee approximate validity of the proposed approach under weak assumptions on the nonconformity score and the ergodicity of the time series. This result is illustrated in great detail by Theorem 2 of chernozhukov2018exact and Theorem 1 of diquigiovanni2021_FTS. We report here the latter, with a slightly modified notation.
Let $Z = Z_{(Y_{T+1})}$, where the candidate function $y$ is now substituted by the random function $Y_{T+1}$. Let $\mathcal{A}^*$ be an oracle nonconformity measure, inducing oracle nonconformity score $S^*$. Define $F$ to be the cumulative (unconditional) distribution function of the oracle nonconformity scores, namely $F(x)=\mathbb{P}(S^*(Z_{(y)}^{\pi})<x)$ and $\hat{F}$ the empirical counterpart, obtained by applying permutations $\pi \in \Pi$: $\hat{F}(x) := \frac{1}{|\Pi|} \sum_{\pi \in \Pi} \mathbbm{1}\{S^*(Z_{(y)}^{\pi}) < x\}$. Let $\{\delta_{1\bar{l}},\delta_{2\bar{m}},\gamma_{1\bar{l}},\gamma_{2\bar{m}}\}$ be sequences of numbers converging to zero.
The first condition concerns the approximate ergodicity of $\hat{F}$, a condition which holds for strongly mixing time series using blocking permutation $\Pi$ defined in (ref) (chernozhukov2018exact). The other conditions are requirements for the quality of approximation of $S^*(Z^{\pi})$ with $S(Z^{\pi})$. Intuitively, $\delta_{2\bar{m}}^2$ bounds the discrepancy between the nonconformity scores and their oracle counterparts. Such condition is related to the quality of the point prediction and to the choice of the employed nonconformity measure.
In order to obtain CP band with empirical coverage close to the nominal one, the choice of an accurate point predictor is important. As mentioned before, whereas in the typical i.i.d. case finite-sample unconditional coverage still holds when the model is heavily misspecified (diquigiovanni2021importance), in the time series context a strong model misspecification may compromise the coverage guarantees and not only the efficiency of the resulting prediction bands (chernozhukov2018exact, diquigiovanni2021_FTS). For this reason, it is important to consider models that are consistent with the functional nature of the observations and that can adequately deal with their infinite dimensionality. We build on top of the literature on functional autoregressive processes in Hilbert spaces, extending them for the first time to temporarily evolving surfaces. We narrow the forecasting methodology to the FAR(p), with $p=1$ because of its wide success in the literature (hernandez2021simultaneous, Papadopoulos2002InductiveCM and Aue). Whereas in the scalar context it is often beneficial to consider lags $p$ greater than one, given the intrinsic high dimensionality of functional data, we would rather fit a biased but simpler model than an unbiased but more complicated model. This issue is enhanced in the two-dimensional context, because of the extra dimension in the domain of the function, and for this reason, we consider only the case $p=1$. We introduce the Functional Autoregressive model of order 1 in Section (ref) and propose estimation techniques in Section (ref).
The most popular statistical model used to capture temporal dependence between functional observations is the functional autoregressive process. The theory of functional autoregressive processes in Hilbert spaces is developed in the pioneering monograph of bosq and a comprehensive collection of statistical advancements for the FAR model can be found in horvath2012inference.
A sequence of mean zero random functions $\{Y_t\}_{t =1}^{T} \subset \mathbb{H}$ follows a non-concurrent Functional Autoregressive Process of order 1 if:
where $\{ \varepsilon_{t} \}_{t \in \mathbb{N}}$ is a sequence of iid mean-zero innovation errors with values in $\mathbb{H}$ satisfying $\mathbb{E}[||\varepsilon_{t}||^2] < +\infty$ and $\Psi$ is a linear bounded operator from $\mathbb{H}$ to $\mathbb{H}$. We consider $\Psi$ to be a Hilbert-Schmidt operator with kernel $\psi$, in such a way that:
In order to ensure existence of a stationary solution of (ref), one has to require that $\exists j_0 \in \mathbb{N}$ such that $||\Psi||^{j_0} < 1$ (bosq, Lemma 3.1).
Proceeding similarly to horvath2012inference, and adopting an approach akin to the Yule-Walker estimation in the scalar setting, we propose the following estimator of $\Psi$:
where $\xi_1, \dots ,\xi_M$ are the first M normalized functional principal components (FPC's), $\lambda_1, \dots, \lambda_M$ are the corresponding eigenvalues, and $\langle x, \xi_1 \rangle, \dots, \langle x, \xi_M \rangle$ are the scores of $x$ along the FPC's. (ref) illustrates two different estimation techniques for $\xi_i$ and $\lambda_i$, one based on a discretization of functions on a fine grid and the other designed starting from an expansion of data on a finite basis system. We further refer to (ref) for more details on the derivation of estimator (ref) and for an extensive discussion on how to adapt it to the Conformal Prediction setting, where $\Psi$ is estimated from the training set $\mathcal{I}_1$ only.
Another forecasting procedure based on FPC's has been proposed by Aue for one-dimensional functional data and is here extended to the two-dimensional setting. Calling once again $\xi_1, \dots, \xi_M$ the first $M$ functional principal components, we decompose the Functional Time Series as follows:
where $\bm{Y}_t=[\langle Y_t, \xi_1 \rangle, \dots, \langle Y_t, \xi_M \rangle]^T$ contains the projection scores, $\bm{\xi}(u,v) = [\xi_1(u,v),\dots,\xi_M(u,v)]^T$ collects the evaluated principal components, and $e_t(u,v)$ is the approximation error due to the expansion's truncation on the first $M$ principal components. Neglecting the approximation error $e_t$, one can prove that the vector $\bm{Y}_t$ follows a multivariate autoregressive process of order 1 (VAR(1)). Plugging in the estimated FPCs $\hat{\xi}_1, \dots, \hat{\xi}_M$, we can estimate the parameters of the resulting VAR(1) model using standard techniques of multivariate statistics and forecast $\hat{\bm{Y}}_{T+1}$ based on historical data $\bm{Y}_1,\dots,\bm{Y}_T$. The predicted function $\hat{Y}_{T+1}$ is then reconstructed as:
We finally introduce a model that may appear simplistic, since it does not exploit the possible time dependence between functions' values in different points of the domain, but that in practical applications provides satisfying results. The prediction method assumes an autoregressive structure in each location $(u,v)$ of the domain, ignoring the dependencies between different points. We call this model a concurrent FAR(1):
where $\psi_{u,v} \in \mathbb{R}$ and $t=2,\dots,T$. Supposing to have observed functional data $y_1, \dots, y_T$ on a common two-dimensional grid $\{(u_i,v_j), i=1,\dots,N_1, j=1,\dots,N_2 \}$, we can estimate $\psi_{u_i,v_j}$ for each location $(u_i,v_j)$.
The goal of this section is twofold: we aim to assess the quality of the proposed CP bands and evaluate different point predictors in terms of the resulting prediction regions. Since this work is focused on uncertainty quantification, we compare forecasting performances by means of the resulting Conformal Prediction bands. Firstly and foremost, we estimate the unconditional coverage by computing the empirical unconditional coverage in order to compare it with the nominal confidence level $1-\alpha$. In the second place, we consider the size of the prediction bands, since a small prediction region is preferable as it includes subregions of the sample space where the probability mass is highly concentrated (Lei1) and it is typically more informative in practical applications.
We employ as a data generating process a FAR(1) model in order to evaluate the estimation routines presented in Section (ref). In order to benchmark forecasting performances, we examine the forecasting methods against a naive one: $\hat{Y}_{T+1}=Y_T$. By including a forecasting algorithm that is not coherent with the data generating process, we can illustrate how the presented CP procedure performs when a good point predictor $g_{\mathcal{I}_1}$ is not available. Although as reported in Section (ref) a sufficiently accurate forecasting algorithm is necessary to guarantee asymptotic validity, we notice that in the simulations CP bands remain valid even when such assumption is not met.
To obtain further insights, we include the performances obtained by assuming perfect knowledge of the operator $\Psi$. For ease of reference, we list here the forecasting algorithms, introducing some convenient notation.
When it is required (namely in FAR(1)-EK, FAR(1)-EK+, FAR(1)-VAR), FPCA is performed using the discretization approach, as motivated in (ref), truncating the representation to the first 4 harmonics.
In Section (ref), we fix the size $b$ of the blocking scheme (ref) equal to 1 and let the sample size $T$ take values $19,49,99,499$. Secondly, in Section (ref), we instead fix the sample size equal to $119$ and repeat the simulations with $b=1,3,6$. As usually done in the time series setting, the first observation is taken into account as a covariate only and does not enter neither the training set nor the calibration set. The proportion of data in the training and in the calibration set are hence equal to one half of the remaining observations: $m=l=(T-1)/2$. For each value of $T$, we repeat the procedure by considering $N=1000$ simulations. Simulations are implemented in the R Programming Language (R).
In order to simulate a sequence of functions $\{Y_t\}_{t=1,\dots,T}$ from a FAR(1), we assume that observations lie in a finite dimensional subspace of the function space $\mathbb{H}$. Without loss of generality, throughout this section we consider functions in $\mathbb{H}=\mathcal{L}^2([0,1]\times[0,1])$. $\mathbb{H}$ is spanned by orthonormal basis functions $\phi_1, \dots, \phi_M$, with $M\in \mathbb{N}$ representing the dimension of such subspace. Therefore, we have:
where $\bm{\phi}(u,v)=[\phi_1(u,v), \dots, \phi_M(u,v)]^T \in \mathbb{R}^M, \, \forall (u,v) \in [0,1]\times[0,1]$, $\bm{Y}_t, \,\bm{\varepsilon}_t \subset \mathbb{R}^M, \, \forall t = 1, \dots, M$ and $\bm{\Psi} \in \mathbb{R}^{M \times M}$ and $\bm{W} \in \mathbb{R}^{K \times K}$, defined as $\bm{W} := \int_c^d \int_e^f \bm{\phi}(u,v) \bm{\phi}(u,v)^T du dv$. It follows that:
The basis system $\phi_1, \dots, \phi_M$ is constructed as the tensor product basis of two cubic B-spline systems $\{g_i\}_{i = 1,\dots, M_1}$, $\{h_j\}_{j = 1,\dots, M_2}$, both defined on $[0,1]$. We set $M_1=M_2=5$, so that $M=25$. For a discussion on the tensor product basis system, we refer to (ref) The matrix $\bm{\Psi}$ is defined as $\bm{\Psi} :=0.7 \frac{\tilde{\bm{\Psi}}}{||\tilde{\bm{\Psi}}||_F}$, with $\tilde{\bm{\Psi}}$ having diagonal values equal to $0.8$ and out-diagonal elements equal to 0.3 Innovation errors $\bm{\varepsilon}_t$ are independently sampled from a multivariate normal distribution, with mean zero and covariance matrix $\bm{\Sigma}$ having diagonal elements equal to 0.5 and out-diagonal entries equal to 0.3. (ref) depicts an example of a simulated Functional Autoregressive Process of order one. A GIF of the time-evolving FAR(1) process can be found on \href{https://github.com/Niccolo-Ajroldi/ARMA-Surfaces/blob/main/Pics/FAR.gif}{GitHub}.
Notice that simulations have been designed in such a way to generate a stationary process. This condition is important to guarantee the existence of a solution to the FAR(1) equation (ref), and is here guaranteed by setting $||\Psi|| < 1$, which satisfies the sufficient condition for stationary presented in Lemma 3.1 of bosq. One can indeed prove that, if relation (ref) holds, then $||\Psi||=||\bm{\Psi}||_F$, where $||.||$ is the usual operatorial norm and $||.||_F$ denotes the Frobenius norm. In this way, FAR(1) estimation techniques are well-defined. For what concerns the theoretical assumptions of the CP scheme, proving the hypothesis of Theorem (ref) is difficult, in particular in the context of functional data. To the best of our knowledge, no test for strongly mixing two-dimensional time series has been proposed, and testing the bounds on the oracle nonconformity scores is even more challenging. However, we aim to show here how the CP procedure can still be applied to obtain valid and efficient prediction bands.
We first fix the size $b$ of the blocking scheme equal to 1 and let the sample size $T$ take values $19,49,99,499$. We replicate the experiments with different significance levels: $\alpha=0.1$, $\alpha=0.2$ and $\alpha=0.3$, in order to assess how the confidence level influences the coverage and the width of the prediction bands.
(ref) shows the empirical coverage, together with the related 99% confidence interval. Empirical coverage is computed as the fraction of the $N=1000$ replications in which $y_{T+1}$ belongs to $\mathcal{C}_{T,1-\alpha}(x_{T+1})$, and the confidence interval is reported in order to provide insights into the variability of the results, rather than to draw inferential conclusions about the unconditional coverage. Notice that different point predictors might intrinsically have dissimilar coverages, consequently this analysis aims to compare forecasting algorithm in terms of their predictive performances. We can appreciate that the 99% confidence interval for the empirical coverage almost always includes the nominal confidence level, regardless of the sample size at disposal. The only exception is obtained with $\alpha=0.3$ and $T=19$. In this case, the method produces very narrow prediction bands ((ref)), that result in an empirical coverage smaller than the nominal one. This behavior however disappears as soon as the sample size $T$ increases. It is also interesting to notice that, even when an accurate forecasting algorithm $g_{\mathcal{I}_1}$ is not available (namely with the Naive predictor), the proposed CP procedure still outputs prediction regions with a high unconditional coverage.
Similarly to diquigiovanni2021importance, we define the size of a two-dimensional prediction band as the volume between the upper and the lower surfaces that define the prediction band:
Measuring the size of the correspondent prediction bands, we can compare the efficiency of different forecasting routines. We stress the fact that distinct point predictors may guarantee potentially different coverage levels. For this reason, it is crucial to first evaluate the empirical coverage of the resulting prediction bands and only afterward investigate their size. (ref) reports boxplots with prediction bands' size for the $N=1000$ simulations and for different values of $\alpha$ and $T$. Bands' size tends to decrease as long as the number of observations $T$ increases, hence improving the efficiency of prediction sets. Moreover, the size tends to decrease when the confidence level $1-\alpha$ increases. As expected, Naive predictor provides larger prediction bands, particularly in the large sample settings. On the other hand, FAR(1)-EK and FAR(1)-EK+, both based on the estimation of the autoregressive operator $\Psi$, provide the tightest prediction bands, not only when numerous observations are available, but also in small sample sizes scenario. Notice also that eigenvalue correction slightly improves the performances of FAR(1)-EK+ wrt FAR(1)-EK, especially when few samples are available. We acknowledge that, when $T=19$, VAR-efpc performs remarkably worse than the other methods. We argue that this performance gap might be caused by the simultaneous OLS estimation of the underling VAR(1) equations, which might provide biased estimates if the sample size is small. However, when the sample size increases, such forecasting algorithm performs comparably with the already mentioned FAR(1)-EK and FAR(1)-EK+. Finally, although the Conformal Prediction bands produced by the oracle predictor are obviously the most performing one, both FAR(1)-EK and FAR(1)-EK+ provide CP bands with coverage and size comparable to the theoretically perfect oracle forecasting method.
This time, we fix the sample size $T$ and let the blocking scheme $b$ vary, in order to determine how the validity and efficiency of the resulting prediction bands are influenced by such parameter. Once again, we repeat the experiments for $\alpha=0.1$, $\alpha=0.2$, $\alpha=0.3$, whereas the value of $T$ is fixed and equal to 119, which provides a good balance between scenarios with small and large sample sizes. Analogous results have been found by letting the value of $T$ vary. Results are reported in (ref) and (ref). Once again, in all circumstances the empirical coverage is close to the nominal one, confirming validity of CP bands even for higher values of $b$. Moreover, one can notice that, as already pointed out by diquigiovanni2021_FTS, the band size tends to decreases when $b$ decreases, thus providing more efficient prediction regions. We argue that this behaviour is related to the inverse proportionality between the blocking scheme size $b$ and the dimension of permutation family $|\Pi|$. Finally, notice that larger prediction bands attain as expected a larger empirical coverage, and that larger values of $\alpha$ (smaller confidence level $1-\alpha$) result in smaller prediction bands. A comparison of forecasting algorithms performances validates the considerations in Section (ref).
In this section, we aim to illustrate the application potential of the proposed methodology on a proper case study. We analyze a data set from Copernicus Climate Change Service (\href{https://climate.copernicus.eu/}{C3S}), a project operated by the European Center for Medium-Range Weather Forecasts (\href{https://www.ecmwf.int/}{ECMWF}), collecting daily sea level anomalies of the Black Sea in the last twenty years (BS_data). Sea level anomalies are measured as the height of water over the mean sea surface in a given time and region. Specifically, altimetry instruments give access to the sea surface height (SSH) above the reference ellipsoid, which is calculated as the difference between the orbital altitude of the satellite and the measured altimetric distance of the satellite from the sea (see (ref)). Starting from this information, Sea Level Anomaly (SLA) is defined as the anomaly of the signal around the Mean Sea Surface component (MSS), which is computed with respect to a 20-year reference period (1993-2012). Observations are collected on a spatial raster, with a resolution of $0.125^{\circ}$ both on the longitude and on the latitude axis. Since observations are collected on a geoid, the domain lies on a two-dimensional manifold, however, because both longitude and latitude ranges are very small ($14^{\circ}$ and $7^{\circ}$ respectively), we assume data to be observed on a bidimensional grid. The resulting lattice can hence be considered as the Cartesian product of a grid on the longitude axis made by $N_1 = 120$ points and a latitude grid of $N_2 = 56$ points. We refer to $(u_i,v_j)$, with $i=1,\dots,N_1$ and $j=1,\dots,N_2$ as the $(i,j)$-th point of such two-dimensional mesh. Since the Black Sea does not have a rectangular shape, we model data as realization of random surfaces defined on the rectangle circumscribed to the perimeter of the sea, but identically equal to zero outside of it. As a consequence, being $\mathcal{B}$ the set of points internal to the perimeter of the Black Sea, we slightly redefine the non conformity measure (ref) to become:
where $\mathcal{R}(u,v)$ is defined as:
Hereafter, we will consider the time series $\{SLA_t\}$, without making explicit the dependence on the bivariate domain.
If possible, one would preferably forecast directly the time series of Sea Level Anomalies ($SLA_t$). However, in order to estimate the FAR(1) process, we would also like to guarantee stationarity of the time series at our disposal, and given the particular nature of the dataset, we expect observations to exhibit a strong seasonality as well as a linear trend. Indeed, both tide gauge and altimetry observations show that sea level trends in the Black Sea vary over time (black_sea_changes, black_sea_linear_trend). Tsimplis estimated a rise in the mean sea level of $2.2$ mm/year from 1960 to the early 1990s, while long-track altimetry data indicate that sea level rose at a rate of $13.4 \pm 0.11$ mm/year over 1993–2008 (black_sea_Ginzburg).
In order to further investigate this issue, we should proceed by testing the functional time series $\{SLA_t\}_t$ for stationarity. However, while for one-dimensional Functional Time Series one could resort to the tests proposed by FTS_stationarity or Aue_fts_test, to the best of our knowledge no stationarity test for two-dimensional functional time series has been developed. For such reason, and aware of the limits of this approach, we resort to analyzing stationarity of univariate time series $SLA_t(u_i,v_j)$, where each $(u_i,v_j)$ represents a grid point in the lattice. We stress the fact that stationarity of individual time series does not guarantee stationarity of the underlying functional process, and this constitutes only a necessary and not sufficient condition. Therefore, the goal is not to derive inferential results on the stationarity of the process, but rather to describe the evolution of the process by means of its individual component, and to potentially obtain a better time series to work with. We test each univariate time series for stationarity, using the Augmented Dickey Fuller (ADF) test. We report in (ref) a grid map of the p-values, for $\{SLA_t\}_t$, $\{\Delta SLA_t\}_t$ and $\{\Delta^2 SLA_t\}_t$ , where $\Delta$ and $\Delta^2$ denote respectively one and two differentiations. Despite the fact that we can't make inferential conclusions on the stationarity of the process based on the individual tests, we can see that the original time series exhibit a very non-stationary behavior, and after one differentiation there are many non-stationary locations. After two differentiations, all the individual time series can be confidently considered stationary. We argue that the Functional Time Series might exhibit a behavior similar to one described in terms of stationarity, and hence proceed by differentiating twice the process, defining
We forecast the differentiated time series $Y_t$, and obtain prediction bands for it, using the methodology presented in Section (ref). However, in order to provide a better insight into the prediction problem, we need to retrieve results pertaining to the original time series $SLA_t$. Specifically, we apply the conformal machinery to the differentiated time series, calling $\hat{Y}_{T+1}$ the forecasted function, and obtaining the prediction band for $Y_{T+1}$:
Exploiting the fact that:
we define the prediction band for $SLA_{T+1}$ as:
Finally, notice once again that testing the hypothesis underlying the CP procedure is not possible in this setting. We are nevertheless interested in applying the CP scheme together with the forecasting techniques in order to verify their empirical performances.
The case study employs a rolling estimation framework which recalculates the model parameters on a daily basis and consequently shifts and recomputes the entire training, calibration and test windows by 24 hours, as shown in (ref). As before, we use a random split of data in the training and calibration sets, with split proportion equal to 50%. The significance level $\alpha$ is once again fixed equal to 0.1. For each of the 1000 days we aim to predict, we build the corresponding prediction band based on the information provided by the last 99 days, thereby fixing $T=99$. Choosing this sample size provides accurate forecasts and thus small prediction bands, while maintaining reasonable computational times. The size of the blocking scheme is fixed to 1, since, as motivated in Section (ref), this choice produces the narrowest prediction bands, preserving at the same time satisfactory performance in terms of empirical coverage. The rolling window is shifted 1000 times, thus iterating for almost three years the forecasting scheme. More specifically, and to allow for reproducibility of subsequent results, we consider a rolling window ranging from 01/01/2017 to 04/01/2020.
The point predictors used throughout this application are those described in Section (ref). The number of Functional Principal Components is set equal to 8. For each shift of the rolling window and for each forecasting algorithm, we check if $SLA_{T+1}$ belongs to $\tilde{\mathcal{C}}_{T,1-\alpha}(x_{T+1})$, and save the size of the corresponding prediction band. After collection of results, we calculate the average coverage, and use it to compare performances of the different point predictors in this scenario.
(ref) outlines observed and forecasted surfaces obtained for one of the day in the rolling window, as long as the lower and upper bounds defining the prediction band. For the sake of simplicity, we display only results obtained using the FAR(1)-EK estimator (ref), since, as discussed below, it provides on average the narrowest prediction bands. In order to allow for a more insightful analysis, and to further investigate the evolution of the surfaces, we implemented a dedicated \href{https://niccolo-ajroldi.shinyapps.io/Black-Sea-Forecasting/}{Shiny App} available online where results can be explored. We report in (ref) the average coverage of CP bands obtained across 1000 predictions. As in Section (ref), we pair such quantity with a 99% confidence interval. Notice that in this case the confidence interval may be biased, due to the inevitable correlation between data used to construct it, however, we include it in order to assess the dispersion of the average coverage around the mean. Coherently with the results of the simulation study, we can appreciate that in all cases the prediction bands capture the observed surface $y_{T+1}$ approximately $(1-\alpha)\%$ of the times, regardless of the forecasting algorithm used.
For what concerns the size of the prediction bands, the Naive predictor produces by far the widest ones (see (ref)), and, this fact does not reflect in a greater coverage compared to the other methods. On the other hand, prediction bands obtained with autoregressive forecasting algorithms provide narrower prediction regions. Among these, we can see that the non-concurrent FAR(1) is the most performing one, regardless of the way in which it is estimated (namely with FAR(1)-EK, FAR(1)-EK+ or FAR(1)-VAR). Nevertheless, also the concurrent FAR(1) model provides very tight prediction bands, almost comparable with the ones produced by the non-concurrent prediction algorithm.
We are also interested in analyzing the pointwise properties of CP bands in this scenario. Therefore, we display in (ref) a map of the pointwise coverage of the prediction bands, obtained using FAR(1)-EK. We can appreciate, as expected, a high empirical coverage across the entire domain, emphasizing once again the peculiarity of our approach, which guarantees global coverage of the prediction surfaces, reflected by an obvious pointwise coverage higher than the nominal one. We report in (ref) the average width of CP bands, which denotes a peculiar pattern, likely caused by data collection routines. Indeed, we can see from (ref), that a similar behaviour observed in the map of pointwise width can be found by plotting the standard deviation of original data. This is coherent with the employed CP framework, since the size of prediction bands depends on the amplitude of the functional standard deviation.
In conclusion, this case study confirms the validity of our procedure and proves how a FAR(1) model significantly improves the predictive efficiency even in this more complex scenario.
In this work, we introduce a mathematical framework for probabilistic forecasting of two-dimensional Functional Time Series. Leveraging the CP scheme developed by chernozhukov2018exact and diquigiovanni2021_FTS and adapting it to the 2D-FTS setting, we propose technique for quantifying uncertainty when predicting time evolving surfaces. In order to provide point prediction of surfaces, we model functions through a Functional autoregressive process of order one, extending the mathematical theory of autoregressive processes to allow for bivariate functions. Estimations techniques for the FAR(1) are presented and compared We test the benefits and limits of the proposed approach, first on synthetic data and then on a real novel time series dataset, collecting daily observations of sea level anomaly over the Black Sea. Empirical results proved the validity of the methodology on non-synthetic data. We acknowledge that in applying the proposed procedure to the case study, we had to introduce some simplifications due to the novelty of the subject and the limited amount of work on two-dimensional Functional Time Series. We hope that this work will encourage the development of novel analysis techniques for 2D functional data, such as stationarity tests and other forecasting tools. Finally, throughout the work, we limited the analysis to the FAR(p), with $p=1$, as motivated in (ref). One may consider a FAR(p) model, with $p>1$. Whereas it is not straightforward to extend the estimation of $\Psi_M$ (ref) to a FAR(p) model with $p>1$, both the concurrent estimator (ref) and the estimator based on an expansion of FPC's, may be easily adapted to the case $p>1$. The problem can in fact be rendered as a standard functional linear regression and solved by using the many off-the-shelf methods apt to the task and present in the literature chiou2016multivariate. Finally, notice that thanks to the flexibility of our approach, a modification of such kind can be easily achieved by changing the point predictor $g_{\mathcal{I}_1}$ only, while preserving the desirable properties of prediction bands.
The present research has been partially supported by MUR, grant Dipartimento di Eccellenza 2023-2027 and by Accordo Attuativo ASI-POLIMI “Attività di Ricerca e Innovazione" n. 2018-5-HH.0, collaboration agreement between the Italian Space Agency and Politecnico di Milano. M.F. acknowledges the support of the JRC Centre of Advanced Studies CSS4P - “Computational Social Science for Policy".