Functional time series prediction under partial observation of the future curve

TL;DR

This paper introduces a novel partial functional prediction method that leverages both complete and partial trajectory data to improve future function predictions in functional time series analysis.

Contribution

The paper develops the PFP method, which incorporates partial observations and includes an automatic tuning criterion, with proven convergence rates.

Findings

PFP outperforms existing methods in simulation studies.

Incorporating partial data improves prediction accuracy.

Method is effective in environmental and traffic flow data analysis.

Abstract

This paper tackles one of the most fundamental goals in functional time series analysis which is to provide reliable predictions for future functions. Existing methods for predicting a complete future functional observation use only completely observed trajectories. We develop a new method, called partial functional prediction (PFP), which uses both completely observed trajectories and partial information (available partial data) on the trajectory to be predicted. The PFP method includes an automatic selection criterion for tuning parameters based on minimizing the prediction error, and the convergence rate of the PFP prediction is established. Simulation studies demonstrate that incorporating partially observed trajectory in the prediction outperforms existing methods with respect to mean squared prediction error. The PFP method is illustrated to be superior in the analysis of environmental data and traffic flow data.

Tables

Table 1. Table 1: fFPE and PMSE values of the three methods under different simulation setups.

${\sigma }_1$
${\kappa }_1$	${\kappa }_2$	fFPEPFP	PMSEPFP	PMSEts	fFPEr	PMSEr
1.8	0.0	0.2024	0.2097	0.8442	0.3269	0.3431
0.8	0.0	0.2003	0.2112	0.8396	0.2664	0.2763
0.2	0.0	0.1928	0.2018	0.8286	0.1938	0.2025
0.4	0.4	0.2038	0.2123	0.8388	0.2309	0.2392
0.0	0.8	0.2058	0.2115	0.8419	0.2647	0.2685
${\sigma }_2$
${\kappa }_1$	${\kappa }_2$	fFPEPFP	PMSEPFP	PMSEts	fFPEr	PMSEr
1.8	0.0	0.5554	0.5801	1.2269	1.1012	1.1668
0.8	0.0	0.5455	0.5640	1.2112	0.7011	0.7431
0.2	0.0	0.5302	0.5561	1.1813	0.5287	0.5536
0.4	0.4	0.5711	0.5985	1.2593	0.6128	0.6391
0.0	0.8	0.5740	0.5907	1.2631	0.6995	0.7127

Table 2. Table 2: The minimum values of fFPE for different pairs of p , d 𝑝 𝑑 p,d , when τ = 8 : 00 : 𝜏 8 00 \tau=8\colon 00 , d x = 6 subscript 𝑑 𝑥 6 d_{x}=6 , d y = 9 subscript 𝑑 𝑦 9 d_{y}=9 , p = 0 𝑝 0 p=0 ; when τ = 12 : 00 : 𝜏 12 00 \tau=12\colon 00 , d x = 7 subscript 𝑑 𝑥 7 d_{x}=7 , d y = 8 subscript 𝑑 𝑦 8 d_{y}=8 , p = 0 𝑝 0 p=0 ; when τ = 16 : 00 : 𝜏 16 00 \tau=16\colon 00 , d x = 8 subscript 𝑑 𝑥 8 d_{x}=8 , d y = 8 subscript 𝑑 𝑦 8 d_{y}=8 , p = 1 𝑝 1 p=1 , d = 3 𝑑 3 d=3 .

$\tau =8:00$
	$d=1$	$d=2$	$d=3$	$d=4$	$d=5$	$d=6$	$d=7$	$d=8$
$p=0$	0.5993	0.5993	0.5993	0.5993	0.5993	0.5993	0.5993	0.5993
$p=1$	0.6278	0.6380	0.6330	0.6494	0.6459	0.6452	0.6462	0.6635
$p=2$	0.6349	0.6591	0.6568	0.6695	0.6742	0.6965	0.7659	0.7933
$p=3$	0.6357	0.6739	0.6412	0.6520	0.6542	0.7130	0.7966	0.8346
$\tau =12:00$
	$d=1$	$d=2$	$d=3$	$d=4$	$d=5$	$d=6$	$d=7$	$d=8$
$p=0$	0.4184	0.4184	0.4184	0.4184	0.4184	0.4184	0.4184	0.4184
$p=1$	0.4274	0.4344	0.4260	0.4498	0.4655	0.4605	0.4485	0.4417
$p=2$	0.4292	0.4460	0.4515	0.4691	0.4933	0.4962	0.5263	0.5441
$p=3$	0.4268	0.4610	0.4385	0.4636	0.4830	0.5045	0.5657	0.5774
$\tau =16:00$
	$d=1$	$d=2$	$d=3$	$d=4$	$d=5$	$d=6$	$d=7$	$d=8$
$p=0$	0.1446	0.1446	0.1446	0.1446	0.1446	0.1446	0.1446	0.1446
$p=1$	0.1517	0.1494	0.1431	0.1472	0.1444	0.1447	0.1525	0.1514
$p=2$	0.1519	0.1490	0.1436	0.1494	0.1453	0.1474	0.1717	0.1744
$p=3$	0.1514	0.1510	0.1535	0.1625	0.1580	0.1675	0.2004	0.1925

Table 3. Table 3: Prediction MSE of the two methods.

Method	$\tau =8:00$	$\tau =12:00$	$\tau =16:00$
moving block	0.56194	0.34591	0.20138
PFP method	0.34789	0.26852	0.10722

Table 4. Table 4: Prediction MSE of the three methods (See Table 2 for the explanation of notations)

$h$	${\text{PMSE}}_n$	${\text{PMSE}}_s$	${\text{PMSE}}_a$	${\text{PMSE}}_i$
1	0.1508980	0.3540901	0.2307175	0.4642256
2	0.2680703	0.4757538	0.6123183	0.6500848
3	0.2309391	0.2576980	0.5273938	0.6132668
4	0.4849306	0.4972889	0.8758383	0.9419479
5	0.3512944	0.4108830	0.9005394	0.9464275
6	0.2363455	0.3411949	0.9953703	0.9732108
7	0.2317724	0.2619626	0.9681887	0.9455279
8	0.2184283	0.2406333	0.9389497	0.8657568
9	0.2376993	0.2154614	0.9931003	0.8339264
10	0.1853883	0.2210090	1.1429001	0.9116115

Equations

$$Y_k(t)=S_k(t)+{ϵ}_k(t),k=1,2,\dots ,$$

$$Y_k(t)=S_k(t)+{ϵ}_k(t)+e_k(t),k=1,\dots ,n.$$

$${\hat{\mu }}_n(t)=\frac{1}{n}\sum _{k=1}^n{Y}_k(t),{\hat{C}}_n(x)=\frac{1}{n}\sum _{k=1}^n⟨Y_k-{\hat{\mu }}_n,x⟩(Y_k-{\hat{\mu }}_n),x\in H.$$

$$Y_k-\mu =\sum _{h=1}^p{\Phi }_h(Y_{k-h}-\mu )+e_k,$$

$$X(s)={\mu }_X(s)+\sum _{i=1}^{\infty }{\xi }_i{\varphi }_i(s)\text{mbox}andY(t)={\mu }_Y(t)+\sum _{j=1}^{\infty }{\zeta }_j{\psi }_j(t),$$

$$Y(t)={\mu }_Y(t)+\int \beta (t,s)(X(s)-{\mu }_X(s))ds+\tilde{e}(t),$$

$${\zeta }_i=\sum _{j=1}^{d_x}{\beta }_{ij}{\xi }_j+{\eta }_i,$$

$$\mathrm{fFPE}(d_x,d_y)=\frac{n+d_x}{n}\mathrm{tr}({\hat{\Sigma }}_{\eta })+\sum _{\ell >d_y}{\hat{\lambda }}_{\ell }^Y,$$

$$Z(s){|}_{[0,\tau ]}=\sum _{j=1}^{\infty }{\xi }_j^{(\tau )}{\varphi }_j^{(\tau )}(s)\text{mbox}andZ(t){|}_{(\tau ,1]}=\sum _{i=1}^{\infty }{\zeta }_i^{(\tau )}{\psi }_i^{(\tau )}(t),$$

$$Z(t){|}_{(\tau ,1]}={\int }_0^{\tau }{\beta }_{\tau }(t,s)Z(s){|}_{[0,\tau ]}ds+\tilde{e}(t).$$

$$Z(t){|}_{(\tau ,1]}=\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{\beta }_{\tau ,ij}{\xi }_j^{(\tau )}{\psi }_i^{(\tau )}(t)+\tilde{e}(t),$$

$$Y_k(t)=S_k(t)+{ϵ}_k(t),t\in [0,1],$$

$$\begin{gathered} Y_k(t_j)=S_k(t_j)+{ϵ}_k(t_j)+e_k(t_j),k=1,\dots ,n, \\ j\ge 1. \end{gathered}$$

$${ϵ}_k(t){|}_{(\tau ,1]}={\int }_0^{\tau }{\beta }_{\tau }(t,s){ϵ}_k(s){|}_{[0,\tau ]}ds+{\tilde{e}}_k(t).$$

$${ϵ}_k(s){|}_{[0,\tau ]}=\sum _{j=1}^{\infty }{\xi }_{kj}^{(\tau )}{\varphi }_j^{(\tau )}(s)\text{mbox}and{ϵ}_k(t){|}_{(\tau ,1]}=\sum _{i=1}^{\infty }{\zeta }_{ki}^{(\tau )}{\psi }_i^{(\tau )}(t).$$

$$Y_k(t_j)=S_k(t_j)+{ϵ}_k(t_j)+e_k(t_j),k\in \mathbb{N},j\ge 1,$$

$${\mathcal{G}}_{\tau }(t,s)=\text{cov}\{{\widetilde{ϵ}}_k(t){|}_{(\tau ,1]},{\widetilde{ϵ}}_k(s){|}_{[0,\tau ]}\},$$

$$C_{ϵ,1}(t,s)=\text{cov}\left\{{\widetilde{ϵ}}_k(t){|}_{[0,\tau ]},{\widetilde{ϵ}}_k(s){|}_{[0,\tau ]}\right\},C_{ϵ,2}(t,s)=\text{cov}\left\{{\widetilde{ϵ}}_k(t){|}_{(\tau ,1]},{\widetilde{ϵ}}_k(s){|}_{(\tau ,1]}\right\},$$

$${\hat{\mathcal{G}}}_{\tau }(t,s)=\frac{1}{n}\sum _{k=1}^n{\widetilde{ϵ}}_k(t){|}_{(\tau ,1]}{\widetilde{ϵ}}_k(s){|}_{[0,\tau ]},$$

$${\hat{C}}_{ϵ,1}(t,s)=\frac{1}{n}\sum _{k=1}^n{\widetilde{ϵ}}_k(t){|}_{[0,\tau ]}{\widetilde{ϵ}}_k(s){|}_{[0,\tau ]},{\hat{C}}_{ϵ,2}(t,s)=\frac{1}{n}\sum _{k=1}^n{\widetilde{ϵ}}_k(t){|}_{(\tau ,1]}{\widetilde{ϵ}}_k(s){|}_{(\tau ,1]}.$$

$$C_{ϵ,1}(t,s)=\sum _{j=1}^{\infty }{\theta }_j{\varphi }_j(t){\varphi }_j(s),C_{ϵ,2}(t,s)=\sum _{i=1}^{\infty }{\omega }_i{\psi }_i(t){\psi }_i(s).$$

$$G={\int }_{\tau }^1{\int }_0^{\tau }{\mathbf{V}}_2^T(t){\mathcal{G}}_{\tau }(t,s){\mathbf{V}}_1(s)dsdt,\hat{G}={\int }_{\tau }^1{\int }_0^{\tau }{\hat{\mathbf{V}}}_2^T(t){\hat{\mathcal{G}}}_{\tau }(t,s){\hat{\mathbf{V}}}_1(s)dsdt,$$

$$\displaystyle B_{11}=\left(\begin{tabular}[]{ccc}$\beta_{\tau,11}$&$\cdots$&$\beta_{\tau,1d_{x}}$\\ $\vdots$&$\ddots$&$\vdots$\\ $\beta_{\tau,d_{y}1}$&$\cdots$&$\beta_{\tau,d_{y}d_{x}}$\end{tabular}\right),$$

$${\hat{\mathbf{V}}}_1(s)=({\hat{\varphi }}_1(s),\dots ,{\hat{\varphi }}_{d_x}(s)),{\hat{\mathbf{V}}}_2(t)=({\hat{\psi }}_1(t),\dots ,{\hat{\psi }}_{d_y}(t)),$$

$$C_y^{-1}{i}^{-\alpha }\le {\lambda }_i\le C_y{i}^{-\alpha },C_1^{-1}{i}^{-{\alpha }_1}\le {\theta }_i\le C_1{i}^{-{\alpha }_1},C_2^{-1}{i}^{-{\alpha }_2}\le {\omega }_i\le C_2{i}^{-{\alpha }_2},$$

$${\delta }_{1i}\equiv {\theta }_i-{\theta }_{i+1}\ge C_1{i}^{-{\alpha }_1-1},{\delta }_{2i}\equiv {\omega }_i-{\omega }_{i+1}\ge C_2{i}^{-{\alpha }_2-1}.$$

$$\sum _{\begin{array}{c}1\le |r_1|<n\end{array}}\sum _{\begin{array}{c}1\le i\le n-|r_1|\end{array}}$$

$$\sum _{\begin{array}{c}1\le |r_1|<n\\ 1\le |r_2|<n\end{array}}\sum _{\begin{array}{c}1\le i\le n-|r_1|\\ 1\le j\le n-|r_2|\end{array}}$$

$$\displaystyle\hskip 48.36958pt\mathcal{Y}_{j}(t_{2},s_{2})\mathcal{Y}_{j+r_{2}}(t_{2},s_{2})\bigg{\}}\bigg{|}dt_{1}ds_{1}dt_{2}ds_{2}=O(n^{2}),$$

$$\sum _{\begin{array}{c}1\le |r_1|<n\\ 1\le |r_2|<n\\ 1\le |r_3|<n\end{array}}\sum _{\begin{array}{c}1\le i\le n-|r_1|\\ 1\le j\le n-|r_2|\\ 1\le k\le n-|r_3|\end{array}}$$

$$\displaystyle\mathcal{Y}_{j+r_{2}}(t_{2},s_{2})\mathcal{Y}_{k}(t_{3},s_{3})\mathcal{Y}_{k+r_{3}}(t_{3},s_{3})\bigg{\}}\bigg{|}dt_{1}ds_{1}dt_{2}ds_{2}dt_{3}ds_{3}=O(n^{3}).$$

$$E{\Vert \int ({\hat{\beta }}_{\tau }-{\beta }_{\tau }){\widetilde{ϵ}}_{n+1}{|}_{[0,\tau ]}\Vert }^2\le \text{mbox}const.n^{\rho }$$

$$\displaystyle\vee\left\{\begin{array}[]{rcl}n^{\{(\frac{\alpha_{1}-2\gamma+2}{\tau_{1}}-1)\vee(\frac{4\alpha_{1}-2\beta_{1}-2\gamma+2}{\tau_{1}}-1)\}},&&\mbox{if}\ {\alpha_{1}-2\gamma>-1,2\alpha_{1}-\beta_{1}-\gamma>-1}\\ n^{\frac{1}{\tau_{1}}-1}\log(n)\vee n^{\frac{4\alpha_{1}-2\beta_{1}-2\gamma+2}{\tau_{1}}-1},&&\mbox{if}\ {\alpha_{1}-2\gamma=-1,2\alpha_{1}-\beta_{1}-\gamma>-1}\\ n^{\{(\frac{1}{\tau_{1}}-1)\vee(\frac{4\alpha_{1}-2\beta_{1}-2\gamma+2}{\tau_{1}}-1)\}},&&\mbox{if}\ {\alpha_{1}-2\gamma<-1,2\alpha_{1}-\beta_{1}-\gamma>-1}\\ n^{\frac{\alpha_{1}-2\gamma+2}{\tau_{1}}-1},&&\mbox{if}\ {\alpha_{1}-2\gamma>-1,2\alpha_{1}-\beta_{1}-\gamma\leq-1}\\ n^{\frac{1}{\tau_{1}}-1}\log(n),&&\mbox{if}\ {\alpha_{1}-2\gamma=-1,2\alpha_{1}-\beta_{1}-\gamma\leq-1}\\ n^{\frac{1}{\tau_{1}}-1},&&\mbox{if}\ {\alpha_{1}-2\gamma<-1,2\alpha_{1}-\beta_{1}-\gamma\leq-1}\end{array}\right..$$

Full text

Keywords: Dimension reduction, Functional principal component, Final prediction error, Functional time series, Intra-day fully functional linear regression model, Long-term and short-term dynamics, Updating prediction.

1 Introduction

Functional data is collected in sociological, environmental, transportation, biological and clinical research. In this paper, we analyze the daily trajectories of the pollutant PM10 in Graz, Austria displayed in Figure 1. These are fine particulate matter with diameter less than 10 micrometers, measured in $\mu$g$/m^{3}$. The task in this paper is to develop a new method that predicts the unobserved part of future trajectories using all the past daily trajectories and the partially observed trajectory of the curve to be predicted.

As noted, functional data are often collected over many natural consecutive time intervals. The PM10 dataset discussed above consists of many daily trajectories. One of the interesting aspects of functional time series is that the trajectories may share similar behavior and that these trajectories may be temporally correlated. Functional time series data is a time series of functions which may be described by $(Y_{k}(t)\ \colon\ k\in\mathbb{N})$, $\mathbb{N}$ denoting the integers, where the observations in the sequence are realizations of random functions $Y_{k}(t)$, and $t$ is in some domain $\mathcal{U}$ which here taken to be the unit interval $[0,1]$. The object $(Y_{k}(t)\ \colon\ k\in\mathbb{N})$ will be referred to as a functional time series.

Prediction methods using only completely observed curves have been discussed in the literature. However, one limitation of these methods is that they are not tailored for the case where there is partially observed data that is available for predicting the remaining part of the curve. Such existing methods often focus on the functional auto-regressive model (FAR). In Bosq (2000), a one-step ahead predictors based on a functional form of the Yule–Walker equations for FAR($1$) processes is derived. Non-parametric kernel predictors were proposed in Besse, Cardot and Stephenson (2000). FAR($1$) curve prediction based on linear wavelet methods were developed in Antoniadis, Paparoditis and Sapatinas (2006). Another line of interesting research in Kargin and Onatski (2008) is the predictive factor method which seeks to replace functional principal components with directions that are most relevant for predictions. Didericksen, Kokoszka and Zhang (2012) evaluated several competing prediction models in a comparative simulation study and report that the Bosq (2000) method has the best overall prediction performance based on mean squared error (MSE) and averaged distance. Aue, Dubart Norinho and Hörmann (2015) proposed a prediction method that first projects the random functions onto the space spanned by the first few functional principal components. A VAR model is then used to predict the truncated eigenscores of the future curve and then transform the predicted eigenscores to functional data by truncated Karhunen-Loéve transformation. In Baranowski and Fryzlewicz (2020), the authors proposed the AMAR method for univariate time series where the AR coefficients are assumed to be constant over different segments. Existing full-curve prediction methods for functional time series only incorporate the dynamics across functions, but not the intra-function information. Thus partial observation is not utilized to improve the prediction. To overcome this limitation, we will propose a method that will incorporate information both across and within trajectories.

In contrast to complete curve prediction, our proposed partial functional prediction (PFP) method yields predictions that utilize both completely observed and partially observed trajectory. Fully functional regression methods for updating time series prediction have been proposed in Chiou (2012), and Yao and Shen (2017). Motivated by the need to predict traffic flow, functional mixture methods were developed which combine fully functional regression, functional clustering, and discrimination. Shang (2017) also considered the fully functional regression method and implemented a moving block approach to update functional time series prediction. However, these approaches have some limitations. The moving block method is developed for prediction over the full interval $\mathcal{U}$, but not for the subinterval of the unobserved partial trajectory. Fully functional regression methods do not take the cross-curve information into account, which would be problematic when there existed strong cross-curve correlation. The PFP prediction method uses all available data – both complete trajectories and partial trajectories. In contrast to methods that use only complete curves, the PFP method has the added flexibility because it can update the prediction according to different time of the interval $\mathcal{U}$. Consequently, by incorporating partially observed trajectory, the prediction error over the forecasting time interval is anticipated to be smaller.

The algorithm for the proposed prediction method is a stepwise procedure and can be summarized as follows. For smooth trajectories, the curves are decomposesd into two parts:

[TABLE]

where $S_{k}(t)$ is the random signal function, and $\epsilon_{k}(t)$ is the zero-mean $i.i.d$ innovation function across $k$. An example of this decomposition is FAR($p$) process, where $S_{k}(t)=\sum\limits_{h=1}^{p}\Phi_{h}(Y_{k-h})$, $\Phi_{h}$’s are the functional autoregressive operators and $Y_{k}$’s are zero-mean random functions without loss of generality. The innovation functions $(\epsilon_{k}(t)\colon k\in\mathbb{N})$ are assumed to be $i.i.d$ for the identifiability of signals and innovations. For $\tau\in(0,1]$, assume that the sub-function over the interval $[0,\tau]$, denoted by $Y_{n+1}|_{[0,\tau]}$, has already been observed and that the goal is to predict $Y_{n+1}|_{(\tau,1]}$. The first step in the proposed method is to use existing functional time series models (e.g. FAR models) to calculate the fitted functions $\widehat{Y}_{k}(t)$ for all $t\in[0,1]$, and obtain the residual functions $\tilde{\epsilon}_{k}(t)=Y_{k}(t)-\widehat{Y}_{k}(t)$. The residual functions are then separated into two segments $\tilde{\epsilon}_{k}|_{[0,\tau]}$ and $\tilde{\epsilon}_{k}|_{(\tau,1]}$ at the current time $\tau$ and then the fully functional regression of $\tilde{\epsilon}_{k}|_{(\tau,1]}$ on $\tilde{\epsilon}_{k}|_{[0,\tau]}$ is employed to predict $\epsilon_{n+1}|_{(\tau,1]}$. The predicted function ${\hat{\epsilon}}_{n+1}|_{(\tau,1]}$ is then used to update the prediction of the $(\tau,1]$ block of the $(n+1)$-th function. The final prediction is derived to be $\widehat{Y}^{u}_{n+1}|_{(\tau,1]}=\widehat{Y}_{n+1}|_{(\tau,1]}+{\hat{\epsilon}}_{n+1}|_{(\tau,1]},$ which is the summation of predictions at each step, where $\widehat{Y}_{n+1}$ is the full-curve prediction. In the case where the observed trajectories are noisy, they are decomposed into three parts:

[TABLE]

In addition to the two steps, a third step is included to extract the dependence information in the random error $e_{k}(t)$ which represents short-term dynamics (i.e., correlation across $t$ within a trajectory). Note that by properly accounting for the dependence between and across trajectories, the proposed PFP method achieves a non-trivial reduction on prediction error.

The remainder of this paper is organized as follows. Section 2 reviews the functional autoregressive model which is widely applied for “next-interval” prediction, and also the fully functional linear regression model and its application in partial functional prediction. In this section, limitations of the these methods are also discussed. In Section 3, we develop the prediction algorithm, and derive the convergence rate of the partial prediction. Simulation studies are discussed in Section 4. In Section 5, the PFP method is validated in the analysis of PM10 concentration and also traffic flow trajectories. Section 6 concludes the article. Additional simulation results and technical details are included in the supplementary material.

2 Functional Models for Prediction

The two popular classes of models for analyzing correlated functional data are functional autoregressive model (FAR) and fully functional regression model. FAR models have been employed for analyzing dependent functional time series data. Fully functional regression models have been used to uncover linear associations between two sets of functions. In this paper, we develop the partial functional prediction (PFP) method which is inspired by the principles of these two approaches.

2.1 Preliminaries

Let $(Y_{k}(t)\colon k\in\mathbb{N})$ be an arbitrary stationary functional time series satisfying the following assumptions:

• (A.1)

All random functions are defined on a common probability space. We assume the observations of $(Y_{k}(t)\colon k\in\mathbb{N})$ are elements of the Hilbert space $H=L^{2}[0,1]$ equipped with the inner product $\langle x,y\rangle=\int_{0}^{1}x(t)y(t)dt$, and $E\{\int_{0}^{1}Y_{k}^{2}(t)dt\}<\infty$.

• (A.2)

The mean function $\mu(t)$ and the covariance operator $C(x)$ of $(Y_{k}(t)\colon k\in\mathbb{N})$ are well defined, and are defined as $\mu(t)=E\{Y(t)\}\mbox{ and }C(x)(t)=\int_{0}^{1}c(t,s)x(s)ds$, where $c(t,s)=E\{Y(t)Y(s)\}$. By the Mercer’s theorem, assume $C(x)$ admits the decomposition $C(x)=\sum\limits_{j=1}^{\infty}\lambda_{j}\langle v_{j},x\rangle v_{j},$ where $(\lambda_{j}\colon j\in\mathbb{N})$ are the eigenvalues in strictly descending order and $(v_{j}\colon j\in\mathbb{N})$ are the corresponding normalized eigenfunctions, so that $C(v_{j})=\lambda_{j}v_{j}$ and $\|v_{j}\|=1$.

The set $(v_{j}(t)\colon j\in\mathbb{N})$ forms orthonormal bases of $L^{2}[0,1]$. Then, by the Karhunen-Loève theorem, $Y_{k}(t)$ allows for the representation $Y_{k}(t)=\mu(t)+\sum\limits_{j=1}^{\infty}\langle Y_{k}-\mu,v_{j}\rangle v_{j}(t).$ The coefficients $\langle Y_{k}-\mu,v_{j}\rangle$ in this expansion are called the fPC scores of $Y_{k}(t)$. Consider now the setting where $Y_{1}(t),\ldots,Y_{n}(t)$ are fully-observed functional time series with unknown mean function $\mu(t)$ and covariance operator $C(x)(t)$. Here, the estimators of $\mu(t)$ and $C(x)(t)$ are

[TABLE]

Remark 1**.**

In (A.1), all functions are defined in the same probability space. Thus, there is common information from all the functions in the random sample that has to be extracted. The eigenfunctions $(v_{j}(t)\colon j\in\mathbb{N})$ in Assumption (A.2) form a series of orthonormal bases for the space $H$. Compared with other orthonormal basis, the eigenfunctions give the best approximation of functions for a given fixed number of basis functions. Under the PFP method, the prediction of functional time series is derived in the subspace spanned by the first few eigenfunctions associated with the largest eigenvalues.

2.2 Multivariate technique of predicting FAR($p$) process

There are multiple existing methods for the prediction of functional time series. In Aue et al. (2015), a dimension-reduction method for prediction of stationary functional time series is proposed which is demonstrated to provide competitive prediction results. To determine the order of the FAR model and the dimension of the auxiliary projected eigenspace, the functional final prediction error criterion was proposed. The FAR($p$) process is defined by the stochastic recursion

[TABLE]

where $(e_{k}\colon k\in\mathbb{N})$ are centered, independent and identically distributed innovations in $L^{2}_{H}$ and $\Phi_{h}\colon H\to H$ are bounded linear operators so that the above recursive equation has a unique causal solution. The multivariate prediction algorithm proceeds in three steps.

Step 1.

Fix the dimension $d$, obtain the empirical fPC scores $y_{k,j}^{e}=\langle Y_{k}-\mu,\hat{v}_{j}\rangle$ for each observation $Y_{k}$, $k=1,\ldots,n$, $j=1,\ldots,d$, and the emperical $d$-variate fPC score vectors ${\bf{Y}}_{k}^{e}=(y_{k,1}^{e},\ldots,y_{k,d}^{e})^{\prime},$ $k=1,\ldots,n$.

Step 2.

Fix the order $p$, and construct a VAR model ${\bf{Y}}_{k}=\sum\limits_{j=1}^{p}\bm{\Phi}_{j}{\bf{Y}}_{k-j}+{\bf E}_{k}$ for the eigenscore vectors to produce the prediction $\widehat{\bf{Y}}_{n+1}^{e}=(\hat{y}_{n+1,1}^{e},\ldots,\hat{y}_{n+1,d}^{e})^{\prime}.$ Durbin–Levinson and innovations algorithm can be readily applied here, given the vectors ${\bf{Y}}_{1}^{e},\ldots,{\bf{Y}}_{n}^{e}$.

Step 3.

The predicted eigenscores are transformed back to a functional trajectory. This transformation is achieved by the truncated Karhunen-Loève representation $\widehat{Y}_{n+1}=\hat{\mu}+\hat{y}_{n+1,1}^{e}\hat{v}_{1}+\cdots+\hat{y}_{n+1,d}^{e}\hat{v}_{d}.$

Based on the predicted fPC scores and the estimated eigenfunctions, the resulting $\widehat{Y}_{n+1}(t)$ is then used as the $1$-step ahead prediction of $Y_{n+1}(t)$.

Remark 2**.**

This method is easy-to-implement for functional time series prediction. However, it only incorporates the correlation across trajectories and is, therefore, inferior to the PFP method when only partial trajectories are to be predicted. The optimal values of $p$ and $d$ are obtained by minimizing the final functional prediction error as developed in Aue et al. (2015).

2.3 Fully Functional Regression Model

2.3.1 General setting

In a fully functional regression model, both the explanatory and response components are functions. Consider random explanatory functions $X(s)$ and response functions $Y(t)$. Denote the mean functions by $\mu_{X}(s)=\mathrm{E}\{X(s)\}$ and $\mu_{Y}(t)=\mathrm{E}\{Y(t)\}$, and their covariance functions by $C_{X}(s_{1},s_{2})=\mathrm{cov}(X(s_{1}),X(s_{2}))$  and  $C_{Y}(t_{1},t_{2})=\mathrm{cov}(Y(t_{1}),Y(t_{2}))$. Suppose that the Karhunen-Loève expansions of the trajectories $X$ and $Y$ are

[TABLE]

where $\xi_{i}$’s and $\phi_{i}$’s ($\zeta_{j}$’s and $\psi_{j}$’s) are the fPC scores and eigenfunctions of $C_{X}$ ($C_{Y}$). The fully functional linear regression model is written as

[TABLE]

where $\tilde{e}(t)$’s are $i.i.d.$ innovation functions with mean zero, and the bivariate regression kernel $\beta(t,s)$ is assumed to be continuous and square-integrable. By the basis expansion, $\beta(t,s)$ admits the representation $\beta(t,s)=\sum\limits_{i=1}^{\infty}\sum\limits_{j=1}^{\infty}\beta_{ij}\phi_{j}(s)\psi_{i}(t).$

Without any loss of generality, the mean function of $X$ and $Y$ are now assumed to be both zero. Replacing $Y(t)$ and $X(s)$ with their Karhunen-Loéve representation leads to $\sum\limits_{i=1}^{\infty}\zeta_{i}\psi_{i}(t)=\sum\limits_{i=1}^{\infty}\sum\limits_{j=1}^{\infty}\beta_{ij}\xi_{j}\psi_{i}(t)+\tilde{e}(t).$ For any $i\in\mathbb{N}_{+}$, taking the inner product with $\psi_{i}(t)$ on both sides yields $\zeta_{i}=\sum\limits_{j=1}^{\infty}\beta_{ij}\xi_{j}+u_{i},$ where $u_{i}=\langle\tilde{e},\psi_{i}\rangle$. In practice, only the first few ($d_{x}$) fPCs are adopted as predictors, so we consider the following equation

[TABLE]

where $\eta_{i}=u_{i}+\sum\limits_{j>d_{x}}\beta_{ij}\xi_{j}$. Then, incorporating the first $d_{y}$ fPC scores of response functions resembles $d_{y}$ multivariate linear regression models. The estimators of $\beta_{ij}$ are obtained by fitting a regression model to the $d_{y}$-dimensional eigenscore vectors of the responses against the $d_{x}$-dimensional eigenscore vectors of the explanatory functions as presented in Equation (2-1). Thus, the eigenscores $\xi$’s and $\zeta$’s are first estimated and then the $\beta_{ij}$’s are estimated by fitting multiple multivariate linear regression models. As for prediction, the first step is to predict the eigenscores of $Y$ with the fitted multiple multivariate linear regression models. Then construct the predicted curve $\widehat{Y}$ by the truncated Karhunen-Loève expansion $\widehat{Y}(t)=\hat{\mu}_{Y}(t)+\sum\limits_{j=1}^{d_{y}}\hat{\zeta_{j}}\hat{\psi}_{j}(t)$.

2.3.2 Dimension selection for fully functional regression model

A common approach is to select the minimum number of principal components so that proportion of variance explained by the functional principal components exceeds a pre-specified threshold. However, the main goal here is prediction and thus it might not be always appropriate to select principal components that explain a large portion of variance. We consider a new criterion for selecting the best dimensions of eigen-spaces of predictors and responses. Here, we propose to choose the dimensions by minimizing the mean squared error (MSE) of prediction, which is asymptotically equivalent to the functional final prediction error (fFPE)

[TABLE]

where ${\lambda}^{Y}_{\ell}$ is the $\ell$-th eigenvalue of $C_{Y}$, ${\Sigma}_{\eta}$ is the covariance matrix of the random vector $(\eta_{1},\ldots,\eta_{d_{y}})$, and $\widehat{\Sigma}_{\eta}$ is the unbiased estimator of ${\Sigma}_{\eta}$. The best $d_{x}$ and $d_{y}$ is the minimizer the fFPE function. The consistency of the fFPE criterion is demonstrated in Section 3.4, and the procedure is derived and explained in detail in the supplementary material.

2.3.3 Intra-curve prediction with functional regression

Define $Z$ being a zero-mean random function in $L^{2}[0,1]$. In a regression setting for intra-curve prediction, for any $0<\tau<1$, the sub-curve over $[0,\tau]$, denoted as $Z(s)|_{[0,\tau]}=(Z(s)\colon s\in[0,\tau])$, serves as the explanatory function, and the sub-curve over $(\tau,1]$, denoted as $Z(t)|_{(\tau,1]}=(Z(t)\colon t\in(\tau,1])$, serves as the response function. Suppose that the Karhunen-Loève expansion of the two functional variables are

[TABLE]

where the notations $\xi^{(\tau)}_{j}$, $\phi^{(\tau)}_{j}(s)$, $\zeta^{(\tau)}_{i}$ and $\psi^{(\tau)}_{i}(t)$ are defined analogously to those on the entire domain $[0,1]$, but they correspond to the sub-domains $[0,\tau]$ or $(\tau,1]$. Then we consider a fully functional linear regression model

[TABLE]

Here, given any fixed value of $\tau$, assume the bivariate regression function $\beta_{\tau}(t,s)$ to be continuous and square integrable. Suppose now that $\beta_{\tau}(t,s)=\sum\limits_{i=1}^{\infty}\sum\limits_{j=1}^{\infty}\beta_{\tau,ij}\phi^{(\tau)}_{j}(s)\psi^{(\tau)}_{i}(t)$, then the functional regression model is expressed as

[TABLE]

where $\beta_{\tau,ij}$ are the unknown regression parameters that will be estimated. Under the continuity assumption on $\beta_{\tau}(t,s)$ along with $\tau$, it follows that $\beta_{\tau,ij}$ is also continuous over $\tau\in(0,1]$ for any $i$ and $j$.

Remark 3**.**

Chiou (2012) employed this method to update prediction for traffic flow trajectories. However, as trajectories in a functional times series can be correlated, it is necessary to include appropriate steps that account for the dependence across trajectories.

3 Partial Functional Prediction (PFP)

3.1 Decomposition of functional time series

As has been discussed in the introduction, an well-established model for smooth functional time series is referred to the following decomposition framework,

[TABLE]

where $S_{k}(t)$ is the signal correlated to the previous trajectories and independent with the future innovations by causality, and $\epsilon_{k}(t)$ is the innovation function that is independent with the previous trajectories. If the observed trajectories are contaminated by random noise, we decompose the functional time series into three parts:

[TABLE]

Note that $\epsilon_{k}(t_{j})$ and $e_{k}(t_{j})$ have different roles in the sense that, $e_{k}(t_{j})$ represents random error accounting for the short-term pointwise dynamics and the roughness of functions, while $\epsilon_{k}(t_{j})$ represents functional innovation accounting for the long-term functional dynamics. In practice, the observations are available only at pre-specified discrete grids, so here we use $t_{j}$ instead of $t$. We propose a step-wise procedure to predict all components. The predicted components are combined to form the final prediction.

3.2 Smooth functions

Suppose that we have observed $Y_{1},\ldots,Y_{n}$, and $Y_{n+1}|_{[0,\tau]}$. The updated prediction of the curve over $(\tau,1]$ is given by $\widehat{Y}^{u}_{n+1}|_{(\tau,1]}=\widehat{Y}_{n+1}|_{(\tau,1]}+\hat{\epsilon}_{n+1}|_{(\tau,1]},$ where $\widehat{Y}_{n+1}$ is the “next-interval” prediction of $Y_{n+1}$ and $\hat{\epsilon}_{n+1}|_{(\tau,1]}$ is the intra-curve prediction of the $(n+1)$-th innovation function over $(\tau,1]$.

To predict ${\epsilon}_{n+1}|_{(\tau,1]}$, we consider a fully functional regression model, where $(\epsilon_{k}(s)|_{[0,\tau]}\colon k\in\mathbb{N})$ serve as the predictors and $(\epsilon_{k}(t)|_{(\tau,1]}\colon k\in\mathbb{N})$ serve as the responses,

[TABLE]

By the Karhunen-Loève expansion,

[TABLE]

The innovation function is unobserved, so the functional regression model is applied to the prediction residuals $\tilde{\epsilon}_{k}=Y_{k}-\widehat{Y}_{k}$ ($k\leq n$) of the previous functions. Replacing the unknown terms with the estimated (predicted) values and adopting the first $d_{x}$ and $d_{y}$ fPCs for predictors and responses respectively lead to the prediction of ${\epsilon}_{n+1}(t)|_{(\tau,1]}$ as follows ${\hat{\epsilon}}_{n+1}(t)|_{(\tau,1]}=\sum\limits_{i=1}^{d_{y}}\sum\limits_{j=1}^{d_{x}}\hat{\beta}_{\tau,ij}\hat{\xi}^{(\tau)}_{n+1,j}\hat{\psi}^{(\tau)}_{i}(t),$ and the final prediction of $Y_{n+1}|_{(\tau,1]}$ is $\widehat{Y}^{u}_{n+1}|_{(\tau,1]}=\widehat{Y}_{n+1}|_{(\tau,1]}+{\hat{\epsilon}}_{n+1}|_{(\tau,1]}.$

The updated prediction $\widehat{Y}^{u}_{n+1}|_{(\tau,1]}$ can be regarded as the complete curve prediction $\widehat{Y}_{n+1}|_{(\tau,1]}$ adjusted by the intra-curve prediction of the $(\tau,1]$ block of the residual function ${\hat{\epsilon}}_{n+1}|_{(\tau,1]}$. The prediction steps are now summarized by the following algorithm.

Step 1.

Fix $d$ and $p$, apply “next-interval” prediction to obtain the prediction of the entire curve $\widehat{Y}_{n+1}$ and the prediction residuals $(\tilde{\epsilon}_{k}(t)\colon k=p+1,\ldots,n)$, $\tilde{\epsilon}_{n+1}(t)|_{[0,\tau]}$.

Step 2.

Segment the prediction residual functions in Step 1 at “current time” $\tau$. Treat the first parts $(\tilde{\epsilon}_{k}|_{[0,\tau]})_{k=p+1}^{n}$ as the predictors, and the second parts $(\tilde{\epsilon}_{k}|_{(\tau,1]})_{k=p+1}^{n}$ as the responses. Fix $d_{x}$ and $d_{y}$, and apply intra-curve functional regression of $(\tilde{\epsilon}_{k}|_{(\tau,1]})_{k=p+1}^{n}$ on $(\tilde{\epsilon}_{k}|_{[0,\tau]})_{k=p+1}^{n}$, and use the fitted model to obtain the prediction of the $(\tau,1]$ block of the $(n+1)$-th innovation function ${\hat{\epsilon}}_{n+1}|_{(\tau,1]}$.

Step 3.

Add the $(\tau,1]$ segment of the complete predicted curve $\widehat{Y}_{n+1}$ and ${\hat{\epsilon}}_{n+1}|_{(\tau,1]}$ to get the final prediction $\widehat{Y}^{u}_{n+1}|_{(\tau,1]}=\widehat{Y}_{n+1}|_{(\tau,1]}+{\hat{\epsilon}}_{n+1}|_{(\tau,1]}$.

3.3 Noisy functions

In this section, we consider functional data as noisy sampled points from a collection of consecutive trajectories. In practice, the observed functional time series is observed at discrete time grids, thus the observed trajectories can be rough. This may be due to measurement errors or sparsely-spaced observation time grids. As discussed in Yao et al. (2005), the random error leads to biased fPC scores. To prevent this problem, the raw trajectories are preprocessed by smoothing. However, in the random error $(e_{k}(t_{j}),k\in\mathbb{N},~{}j\geq 1)$, which is not smooth, there can still exist short-term temporal dependence, thus additional steps are necessary to utilize the dependence information in $(e_{k}(t_{j}),k\in\mathbb{N},~{}j\geq 1)$ to improve the prediction.

As previously discussed, any noisy functional time series $(Y_{k}(t),k\in\mathbb{N})$ can be decomposed into three parts,

[TABLE]

where $S_{k}(t_{j})$ is the smooth signal depending on the smooth part of the past functions, $\epsilon_{k}(t_{j})$ is the independent smooth innovation function, and $e_{k}(t_{j})$ is the random error. Let $f_{k}(t_{j})=S_{k}(t_{j})+\epsilon_{k}(t_{j})$ represent the smooth part of the functional time series, which can be estimated by smoothing techniques (e.g. kernel smoothing, basis smoothing, local polynomial regression). Now fit the ARMA model to the smoothing residuals $(r_{k}(t_{j}))_{k,j\geq 1}$, defined as $r_{k}(t_{j})=Y_{k}(t_{j})-\tilde{f}_{k}(t_{j})$, where $(\tilde{f}_{k}(t)\colon k\in\mathbb{N})$ are the smoothed trajectories. For noisy trajectories, there are additional steps in the algorithm.

Step 4.

Produce the smoothed curve $(\tilde{f}_{k}(t)\colon k\in\mathbb{N})$ with any smoothing technique and denote the smoothing residuals to be $(r_{k}(t_{j}))_{k,j\geq 1}$. Apply ARMA model to the smoothing residuals to predict the future random errors $\hat{r}_{n+1}(t_{j})$, $t_{j}>\tau$.

Step 5.

Combine the prediction of the smooth part in Step 3, $\hat{f}_{n+1}(t)$, and the prediction of the smoothing residuals to obtain the final prediction $\widehat{Y}^{u}_{n+1}(t_{j})=\hat{f}_{n+1}(t_{j})+\hat{r}_{n+1}(t_{j}).$

This adjustment is necessary if the observed trajectories are significantly rough and the dependence of $(e_{k}(t_{j}))_{k,j\geq 1}$ across $t_{j}$ is pronounced. The auto-correlation of the smoothing residuals typically decays much faster than that of the original time series, which indicates that the long-term dynamics (e.g. seasonal trend) is removed.

3.4 Convergence rate and criterion consistency

In this section, we develop the upper bound of the convergence rate of the partial prediction and the consistency of the fFPE criterion. The convergence rate is related to the fPC scores of the current partial residual $\tilde{\epsilon}_{n+1}(t)|_{[0,\tau]}$, the correlation between the two segments of prediction residuals, the decay rates of the eigenvalues of the covariance function of the first-part residual $\tilde{\epsilon}_{k}(t)|_{[0,\tau]}$ and the second-part residual $\tilde{\epsilon}_{k}(t)|_{(\tau,1]}$, where $\tilde{\epsilon}_{k}(t)=Y_{k}-\mu-\sum\limits_{h=1}^{p}\Phi_{h,d}(Y_{k-h,d}-\mu),$ and $\Phi_{h,d}$, $Y_{k,d}$ are the $d$-dimensional projection of $\Phi_{h}$, $Y_{k}$ on $\{v_{j}\colon j=1,\ldots,d\}$. Define

[TABLE]

[TABLE]

and the estimators are

[TABLE]

[TABLE]

By the Mercer’s theorem, suppose that $C_{\epsilon,1}(t,s)$ and $C_{\epsilon,2}(t,s)$ admit the following spectral decomposition

[TABLE]

Consider the intra-curve regression model $\tilde{\epsilon}_{k}(t)|_{(\tau,1]}=\int\beta_{\tau}(t,s)\tilde{\epsilon}_{k}(s)|_{[0,\tau]}ds+\delta_{k}(t)$ and represent $\beta_{\tau}(t,s)$ with $\phi_{j}(s)$’s and $\psi_{i}(t)$’s as $\beta_{\tau}(t,s)=\sum\limits_{i,j=1}^{\infty}\beta_{\tau,ij}\psi_{i}(t)\phi_{j}(s).$ Notationally let $\bm{V}_{1}(s)=(\phi_{1}(s),\ldots,\phi_{d_{x}}(s)),\bm{V}_{2}(t)=(\psi_{1}(t),\ldots,\psi_{d_{y}}(t))$ and

[TABLE]

[TABLE]

Note that $G$ and $\widehat{G}$ are $d_{y}\times d_{x}$ matrices. A truncated estimator of the unknown coefficient $\beta_{\tau}(t,s)$ can be represented as $\hat{\beta}_{\tau}(t,s)=\widehat{\bm{V}}_{2}(t)\widehat{B}_{11}\widehat{\bm{V}}^{T}_{1}(s),$ where

[TABLE]

and $\widehat{B}_{11}=\widehat{G}\widehat{\Theta}^{-1}.$ We represent $\tilde{\epsilon}_{n+1}(t)|_{[0,\tau]}$ as $\tilde{\epsilon}_{n+1}(t)|_{[0,\tau]}=\sum\limits_{j\geq 1}x_{j}\phi_{j}(t)$, where $x_{j}=\langle\tilde{\epsilon}_{n+1}|_{[0,\tau]},\phi_{j}\rangle$. Here $\tilde{\epsilon}_{n+1}(t)|_{[0,\tau]}$ is treated as a fixed function. The following assumptions are made.

• (A.3)

With some $\alpha,\alpha_{1},\alpha_{2}>1$, $C_{y},C_{1},C_{2}>0$, the eigenvalue sequences $(\lambda_{i}\colon i\in\mathbb{N})$, $(\theta_{i}\colon i\in\mathbb{N})$ and $(\omega_{i}\colon i\in\mathbb{N})$ are strictly decreasing, and satisfy the conditions below

[TABLE]

• (A.4)

There exists some positive constants $\beta_{1}$, $\beta_{2}$, and $\gamma$, so that $\lvert G_{ij}\rvert\equiv\lvert\langle\int\mathcal{G}_{\tau}\phi_{j},\psi_{i}\rangle\rvert\leq\mbox{const.}i^{-\beta_{2}}j^{-\beta_{1}}$ and $\lvert x_{j}\rvert\leq\mbox{const.}j^{-\gamma},$ where “const.” denotes a constant.

• (A.5)

$\gamma>1/2$, $\alpha_{1}-\beta_{1}<-\frac{1}{2}$, $\alpha_{2}-\beta_{2}<-\frac{1}{2}$, $\alpha_{1}>\max\{1,(2\gamma-3)/4\}$, $\alpha_{2}>1$.

• (A.6)

Define $\mathcal{Y}_{k}(t,s)=\tilde{\epsilon}_{k}(t)\tilde{\epsilon}_{k}(s)-E\{\tilde{\epsilon}_{k}(t)\tilde{\epsilon}_{k}(s)\}$, and we assume that $p$ and $d$ are selected so that

[TABLE]

• (A.7)

Each $Y_{k}$ admits the representation $Y_{k}=f(u_{k},u_{k-1},\ldots),$ where $(u_{k}\colon k\in\mathbb{N})$ are $i.i.d$ elements taking values in a measurable space $S$, and $f$ is a measurable function $f\colon S^{\infty}\to H$. Moreover, we assume that if $(u^{\prime}_{k}\colon k\in\mathbb{N})$ is an independent copy of $(u_{k}\colon k\in\mathbb{N})$ defined on the same probability space, then letting $Y_{k}^{(m)}=f(u_{k},u_{k-1},\ldots,u_{k-m+1},u^{\prime}_{k-m},u^{\prime}_{k-m-1},\ldots),$ we have $\sum\limits_{k=1}^{\infty}(E\|Y_{k}-Y_{k}^{(k)}\|_{4}^{4})^{1/4}<\infty$, where $\|\cdot\|_{4}$ signifies the $\ell_{4}$-norm.

From Assumption (A.3) and (A.4), it follows that $\lvert\beta_{\tau,ij}\rvert\leq\mbox{const.}i^{-\beta_{2}}j^{\alpha_{1}-\beta_{1}}.$ Assumption (A.3) assures the identifiability of principal components and consistency of estimation (see e.g. Gohberg, I. and Krupnik, N., 1992). Assumption (A.4) quantifies the correlation between the two parts of residual functions. Assumption (A.5) assures the square integrability of $\tilde{\epsilon}_{n+1}|_{[0,\tau]}$ and $\beta_{\tau}(t,s)$, and $\alpha_{2}-\beta_{2}<-1/2$ is an analogical assumption of $\alpha_{1}-\beta_{1}<-1/2$, which indicates that if the predictor and response are switched, the intra-curve regression model is still well defined. $\alpha_{1}>(2\gamma-3)/4$ restricts that the fPC scores of the current partial residual does not decay much faster than average, and the assumptions $\alpha_{1},\alpha_{2}>1,\gamma>1/2$ come from the assumption $Y_{k}\in L^{2}[0,1]$. Assumption (A.6) and (A.7) assure that the weak dependence across the prediction residual functions does not influence the convergence rate, which are naturally true for $m$-dependent processes.

The upper bound of the convergence rate of the partial prediction is presented in the following theorem.

Theorem 1**.**

Set $d_{x}=O(n^{1/\tau_{1}})$ and $d_{y}=O(n^{1/\tau_{2}})$, where $\tau_{1}>2(\alpha_{1}+1)$ and $\tau_{2}>\max\{2(\alpha_{2}+1),(2\alpha_{2}+3)\tau_{1}/(2\tau_{1}-1)\}$. Under Assumptions (A.3)—(A.7), we can find $\alpha_{1}$ and $\alpha_{2}$ (where $\alpha_{1},\alpha_{2}$ are related to $\tau$), such that $\alpha_{1}>\alpha$ and $\alpha_{2}>\alpha$ for any $\tau\in(0,1)$. Let $\rho=\max\{\frac{2(\alpha_{1}-\beta_{1}-\gamma+1)}{\tau_{1}}\vee(\frac{-2\beta_{2}+1}{\tau_{2}})\vee(\frac{1}{\tau_{1}}+\frac{2\alpha_{2}+3}{\tau_{2}}-2)\vee(\frac{4\alpha_{1}-2\gamma+4}{\tau_{1}}-2)\}$, then

[TABLE]

Remark 4**.**

In Theorem 1, “$a\vee b$” signifies $\max\{a,b\}$. A practical implication of $\alpha_{1}>\alpha$ and $\alpha_{2}>\alpha$ is that, compared to the entire function, partial trajectories can be approximated with less principal components, while equal proportion of variance is explained. The assumptions on $\tau_{1}$ and $\tau_{2}$ guarantee the convergence of prediction and uniform consistency. Otherwise, if $d_{x}$ and $d_{y}$ were selected too large, the prediction would diverge.

The following theorem demonstrates the consistency of the fFPE criterion in Section 2.3.2. Define $\Delta=n^{-1}\sum\limits_{k=1}^{n}\tilde{\epsilon}_{k}(t)\tilde{\epsilon}_{k}(s)-E\{\tilde{\epsilon}_{k}(t)\tilde{\epsilon}_{k}(s)\}$, then if the following assumptions hold, we can develop the consistency property of the fFPE criterion.

• (A.8)

$E\|\Delta\|^{16}=O(n^{-8})$, $n^{-1}d^{3\alpha_{2}+2}_{y}\to 0$, $n^{-1}d_{y}\left(\sum\limits_{i=1}^{d_{y}}\delta^{-4}_{2i}\right)^{1/2}\to 0$, and

[TABLE]

Theorem 2**.**

Let $\hat{\epsilon}^{u}_{n+1}|_{(\tau,1]}$ be the prediction of $\tilde{\epsilon}_{n+1}|_{(\tau,1]}$ based on $C_{\epsilon,1}$ and $C_{\epsilon,2}$, and $\tilde{\epsilon}^{u}_{n+1}|_{(\tau,1]}$ be the prediction of $\tilde{\epsilon}_{n+1}|_{(\tau,1]}$ based on $\widehat{C}_{\epsilon,1}$ and $\widehat{C}_{\epsilon,2}$. Then, under Assumptions (A.3), (A.6) and (A.8), we have

[TABLE]

uniformly for $d_{x}$ and $d_{y}$ as $n\to\infty.$

Remark 5**.**

In practice, $\Phi_{h}(\cdot)$’s are unknown and need to be estimated. However, if Assumptions (A.3)—(A.8) still hold for the estimated residuals, the result still follows.

4 Simulation

4.1 General setting

To analyze the finite sample properties of PFP, a comparative simulation study was conducted. The PFP method was tested on simulated FAR models. In each simulation test, 400 trajectories were generated by a FAR($p$) model ($p=1$ or $2$). The $(\tau,1]$ blocks of the last 20 trajectories were predicted. The corresponding mean squared error (MSE) of prediction was computed, as well as the fFPE value for comparison. This procedure was repeated for 100 times for each simulation setup.

The simulation was conducted under the context of a $D$-dimensional functional space which is spanned by the Fourier basis functions $\bm{\nu}=(\nu_{1},\nu_{2},\ldots,\nu_{D})$ on the unit interval $[0,1]$ that correspond to the first $(D-1)/2$ fundamental Fourier frequencies. Thus, any simulated function admits the representation $x(t)=\sum\limits_{j=1}^{D}c_{j}\nu_{j}(t)$ with coefficients $\bm{c}=(c_{1},\ldots,c_{D})^{\prime}$. Then for any linear operator $\Psi\colon H\to H$,

[TABLE]

where ${\bf\Psi}$ is a $D\times D$ matrix with elements $\{\langle\Psi(\nu_{j}),\nu_{j^{\prime}}\rangle\}_{j,j^{\prime}=1}^{D}$. The innovation function admits the basis expansion $\epsilon_{k}(t)=\sum\limits_{j=1}^{D}a_{k,j}\nu_{j}(t)$, where $a_{k,j}$’s are $i.i.d.$ normal random variables with mean zero and standard deviation $\sigma_{j}$. Two sets of standard deviations used here are ${\bf\sigma}_{1}=(j^{-1}\colon j=1,\ldots,D)$ and ${\bf\sigma}_{2}=(1.2^{-j}\colon j=1,\ldots,D)$.

4.2 Prediction comparison for smooth trajectories

In this section, we show the comparison of partial functional prediction with Aue et al. (2015)’s method and the intra-curve functional regression method on FAR(2) processes $Y_{k}=\Psi_{1}(Y_{k-1})+\Psi_{2}(Y_{k-2})+\epsilon_{k}$. The operators were generated such that $\Psi_{1}=\kappa_{1}\Psi$ and $\Psi_{2}=\kappa_{2}\Psi$. (Here, note that $\kappa_{2}=0$ yields a FAR(1) process). The operator matrix ${\bf\Psi}$ was generated at random, with each element following a normal distribution with mean zero and variance $\sigma_{jj^{\prime}}$, and then scaled by its $l_{2}$-norm. We set $\sigma_{jj^{\prime}}$ to be $({\bf\sigma}_{i}{\bf\sigma}_{i}^{\prime})_{jj^{\prime}}$ to ensure the simulated functions satisfying Riemann-Lebesgue Lemma, and set $D=15$.

In each simulation setup, the MSEs of prediction $\int_{\tau}^{1}\{Y_{n+1}(t)-\widehat{Y}^{u}_{n+1}(t)|_{(\tau,1]}\}^{2}dt$ of PFP and the two competitor methods were computed. The fFPE values were also calculated for PFP and the intra-curve regression, which were recorded to be close to the corresponding MSE of prediction. Results for five pairs of values $(\kappa_{1},\kappa_{2})$ are provided in Table 1. For clarity, the following abbreviations are introduced: PMSE=“prediction mean squared error”, fFPE=“final functional prediction error”. The subscript: $ts$=“FAR model prediction (Aue et al. (2015))”, $r$=“intra-curve regression prediction”.

Clearly, when the signal-noise ratio is high, PFP outperforms the other methods. Otherwise, the performances of PFP and the intra-curve regression method are similar, and are both better than full-curve prediction. The fFPE value and the prediction MSE are always very close for different situations, which numerically justifies the practical applicability of the fFPE criterion. The results of bootstrap prediction interval are included in the supplementary material.

4.3 Prediction comparison for noisy trajectories

Noisy functional time series were simulated by adding errors following AR($1$) process to the smooth functional time series. Specifically, the simulated functions admit the decomposition $Y_{k}(t_{j})=S_{k}(t_{j})+e_{k}(t_{j}),\ j=1,\ldots,48,$ where $S_{k}(t_{j})$ is the smooth curve obtained from the simulated FAR($1$) process ($\kappa_{1}=1.8$, $\kappa_{2}=0$), and $e_{k}(t_{j})$ is the random error. The “current time” $\tau=0.5$.

To incorporate the random errors into prediction, two other approaches were also considered: 1). ARIMA model, and 2). applying only Step 1–3 to the original time series without smoothing (or equivalently, smoothing the functions by linear interpolation). The average prediction MSE of the following 5 grids ($1\leq h\leq 5$) of the last 20 trajectories are shown in Table 2. The simulation experiments indicate that the ARIMA model should be the “last-resort” method for long-term prediction. The ARIMA model may provide reasonable short term prediction, and one can use this approach to predict the rough errors. However, if the error term is incorporated into PFP in the smooth case by linear interpolation, the prediction will deteriorate. This is because that the estimation of the actual fPC scores is biased, and this error propagates to the estimated FAR model.

5 Analysis of PM10 and Traffic flow trajectories

5.1 Analysis of PM10 concentration

One goal of this section is to analyze the trajectories of PM10 concentration which broadly refers to particulate matter with an aerodynamic diameter of less than 10$\mu m$ in ambient air. The trajectories were measured every 30 minutes in Graz, Austria. Prior to applying the proposed PFP prediction method, the data was segmented according to the day of the week. The 48 observations for each day were combined into a vector. Visual inspection of the data revealed several extreme outliers around New Year’s Eve known to be caused by firework activities. The corresponding week was removed from the sample. To stabilize the variance, the square root transformation was applied to the recordings. The remaining 175 transformed trajectories are displayed in Figure 2.

Then the discrete vectors were transformed into functional objects using 10 cubic B-spline basis functions and the least squares fitting. To remove weekly seasonality, the mean functions were computed for each day of the week and were used to centralize the trajectories before prediction.

5.1.1 Prediction of the smoothed PM10 concentration

Suppose that the current time in a day is $\tau$, where $\tau\in(0,24]$, and there is partially-observed curve whose latter trajectory is to be predicted (i.e., the observation is only available over $[0,\tau]$ and the remainder of the curve on $(\tau,24]$ needs to be predicted). The one-step ahead prediction was conducted and the corresponding fFPE was computed. The minimum values of fFPE (over $d_{x}$ and $d_{y}$) are shown in Table 3 corresponding to different values of $p$ and $d$. Figure 3 shows the partial prediction of two randomly selected curves as $\tau=8\colon 00,\ 12\colon 00,\ 16\colon 00$ respectively. Note that the prediction residual functions are not necessarily centered at zero, and thus the mean function needs to be adjusted when computing the intra-curve prediction. The final prediction is

[TABLE]

where $\hat{\mu}_{n+1}$ is the estimated weekday mean function, and $\hat{\mu}_{e}$ is the estimated mean function of the prediction residuals in the first step.

5.1.2 Comparison with moving block method

Shang (2017) proposed a functional time series prediction method, called the moving block method, to update the prediction with switching $\tau$. The time support is shifted forward by $\tau$. Specifically, the $(\tau,1]$ block of the $m$-th curve is combined with the $[0,\tau]$ block of the $(m+1)$-th curve to form a new function. The new functions are a recombination of the original functional time series with the loss of the $[0,\tau]$ part of the first curve, which typically lays trivial effect on the prediction. The full-curve prediction method is then applied to the new functional time series, and the $[0,\tau]$ block of the predicted function is the updated prediction.

Table 4 includes the prediction MSE of the last 20 trajectories by PFP and the moving block method. It is noted that PFP robustly outperforms the moving block method over a broad range of values of $\tau$. The result is not unexpected since the moving block method actually belongs to “next-interval” prediction method, which provides complete curve prediction, while PFP aims to produce prediction only for the unobserved block, so the prediction error of the unobserved block provided by PFP should be smaller than that of the moving block method. Indeed one of the advantages of PFP over the moving block is that it directly uses the intra-curve variation, that is, the partially observed trajectory is directly treated as a part of the trajectory of interest rather than treating it artificially as part of the previous curve. A severe limitation of the moving block method is that the partially observed curve for the current trajectory is artificially forced to be a part of the previous one. This principle might work for some data settings but will not be reasonable for many other biological settings where the start of the curve has a well defined meaning (such as the onset of a stimulus presentation or a shock in biological experiments).

5.1.3 Prediction of the original trajectories

Since the PM10 trajectories are not smooth and present seasonal dynamics, it is natural to implement the PFP method for the noisy case. The prediction results are compared with ARIMA model prediction, and PFP for smooth case was also implemented for comparison. We also applied linear interpolation when smoothing the original trajectories to incorporate the random errors, and then applied PFP for smooth case to finalize the prediction.

The current time $\tau$ was assumed to be $12\colon 00$, say the first 24 values were observed. The prediction methods were applied to predict the $h$-step ahead point values for the last 25 trajectories, where $1\leq h\leq 10$. Table 5 shows the prediction error of the three methods. Figure 4 shows part of the original centered time series and the corresponding smoothing residuals, and it is noted that after removing the smoothed functions, the residuals have no obvious seasonal trend compared with the original time series.

Table 5 indicates that there exists dependence across the smoothing residuals since Step 4 and 5 significantly improve the prediction. These results demonstrate that PFP captures both the short-term dynamics (across smoothing residuals) and long-term dynamics (across and within smoothed functions). The ARIMA model can only give good predictions for the short-term predicted values but cannot give accurate predictions if we are interested in the long-term future. Linear interpolation does not perform well since the random errors contaminate the smooth signals and lead to bias in the estimated functional principal components.

5.2 Analysis of traffic flow trajectories

We now analyze the traffic flow data that was collected by a dual loop vehicle detector near the Shea-San Tunnel on National Highway 5 in Taiwan in 2009 (shared by Chiou (2012)). It refers to the vehicle count per minute over 15-min time intervals (96 observations for each day). There are 92 days of observed trajectories in total, and the goal is to predict the unobserved block of the last 12 trajectories. Figure 5 shows the raw daily trajectories and smoothed daily trajectories (smoothed with 21 B-splines).

Chiou (2012) proposed a functional mixture prediction method for independent trajectories which were first classified into several clusters. Intra-curve regression model within each cluster was applied to predict the unknown block in each potential cluster. The predictions in each cluster were combined to form the final prediction. Here, we used the first 80 trajectories as the training set to determine the cluster membership by subspace projection cluster algorithm (see Chiou and Li (2007)), and the last 12 trajectories are classified based on only the $[0,\tau]$ block.

In the testing set (the last 12 trajectories), for a sample $Y_{k}(t)$ observed up to $\tau$, we used the mean integrated prediction error (abbreviated as MIPE) to measure the performance of different methods

[TABLE]

Figure 6 shows the MIPE of the two methods. The result shows that proposed PFP method yields lower MIPE compared to the functional mixture prediction method. Although the functional mixture prediction method can work well in some cases, it has some limitations. First, the method classifies the future curve only based on the observed part, however, when the observed part is not very representative of the whole curve, the curve to be predicted is likely to be classified into a wrong cluster, which will potentially increase the prediction error. Furthermore, applying functional linear regression in each cluster actually reduces the sample size, resulting in larger estimation error.

6 Conclusion

The contribution of this paper is a new functional prediction methodology that provides an update on the prediction given that the curve to be predicted is partially observed. The PFP method is motivated by the idea that the updated prediction should be a projection onto the $\sigma$-algebra expanded by the past observed trajectories and the partial observation. The prediction algorithm is a step-wise procedure, and can be applied to both smooth and non-smooth functions. In non-smooth case, the functional techniques can be applied for removing the seasonal trend, and univariate time series models can be applied to predict the smoothing residuals more effectively.

There are already several prediction methods for functions which we summarize here. In the functional time series prediction method (e.g. Aue et al. (2015)), the “next-interval” prediction only considers the big picture of the next function. In the setting where there is available partial observation, it is more natural to use all the available information (in particular, intra-curve information) to predict the unobserved part in order to improve prediction. The primary limitation of the existing full-curve prediction methods is that they do not incorporate this available information. Another method, the moving block method (Shang, 2017), is essentially “next-interval” prediction method, so it has the same limitations discussed above. Another limitation of this method is that it is unnatural to arbitrarily assign starting and ending points of a curve especially in studies where such points are explicitly determined (e.g., start of the day; start of a trial in an experiment). As for the fully functional regression method (see e.g. Ramsay and Silverman (2005)), while it is commended for incorporating intra-curve information, its limitation is that it does not take into account the correlation across trajectories. This is a serious issue when the dependence across functions is non-trivial and the past trajectories are highly informative for predicting future trajectories. The method of Chiou et al. (2012) is an extension of functional regression which smartly and intuitively combines functional regression and clustering. The limitation of this approach is that when the partial observation does not give a strong indication of cluster membership, the classification will not be reliable and this could lead to serious prediction errors. It can be challenging to classify time series with low signal-to-noise ratio or with short time series length. Moreover, the sample size is potentially reduced (per cluster) as we need to do estimation for each cluster separately. On the general approach of functional time series prediction after smoothing the curve by linear interpolation, although this is a way to jointly incorporate the information of both long-term and short-term dynamics, random errors will be included in the obtained curve and that can result in bias.

The PFP method has several advantages. Since functional data are usually obtained in consecutive time intervals, the time series structure (e.g., autocorrelation) ubiquitously exists in functional data, and thus the PFP method has a broad range of applicability. Additionally, it is flexible, in the sense that, users can decide whether or not the time series structure should be taken into account when predicting data. The PFP method is developed for predicting the unobserved block, so comparing with the full-curve prediction (including the moving block method), the prediction error of the unobserved block should be smaller. The simulation studies and real data analysis demonstrate that PFP consistently yields superior and competitive prediction.