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.
114,256 characters · 14 sections · 36 citation commands
On the Modelling and Prediction of High-Dimensional Functional Time Series
\if11 { \spacingset{1.25}
\affil[1]{\it Joint Laboratory of Data Science and Business Intelligence, Southwestern University of Finance and Economics, Chengdu, Sichuan, China} \affil[2]{\it Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China} \affil[3]{\it Business School, University of Sydney, Sydney, Australia} \affil[4]{\it Faculty of Business and Economics, The University of Hong Kong, Hong Kong } \affil[5]{\it Department of Statistics, London School of Economics, London, U.K. }
\setcounter{Maxaffil}{0}
} \fi \if01 {
} \fi \spacingset{1.5}
\noindentKeywords: Dimension reduction; Eigenanalysis; Functional thresholding; Hilbert--Schmidt norm; Permutation; Segmentation transformation.
\spacingset{1.69} {0.2\baselineskip} {0.2\baselineskip} {0.2\baselineskip} {0.2\baselineskip}
Functional time series typically refers to continuous-time records that are naturally divided into consecutive time intervals, such as days, months or years. With recent advances in data collection technology, multivariate or even high-dimensional functional time series arise ubiquitously in many applications, including daily pollution concentration curves over different locations, annual temperature curves at different stations, annual age-specific mortality rates for different countries, and intraday energy consumption trajectories from different households. Those data can be represented as a $p$-dimensional functional time series ${\mathbf Y}_t(u)=\{Y_{t1}(u), \ldots, Y_{tp}(u)\}^{{ \mathrm{\scriptscriptstyle \top} }}$ defined on a compact set $u\in{\mathcal U}$, and we observe ${\mathbf Y}_t(\cdot)$ for $t=1, \ldots, n$. In this paper we tackle the high-dimensional settings when the dimension $p$ is comparable to, or even greater than, the sample size $n$, which poses new challenges in modelling and forecasting ${\mathbf Y}_{t}(\cdot)$.
By assuming ${\mathbf Y}_t(\cdot)$ is stationary, a conventional approach is first to extract features by performing dimension reduction for each component series $Y_{tj}(\cdot)$ separately via, e.g. functional principal component analysis (FPCA) or dynamic FPCA Bathia2010,hormann2015, and then to model $p$ vector time series by, e.g., regularized vector autoregressions guo2021 or factor model Gao2019. However, more effective dimension-reduction can be achieved by pulling together the information from different component series in the first place. This is in the same spirit of multivariate FPCA Chiou2014,Happ2018 (for fixed $p$) and sparse FPCA hu2021, though those approaches make no use of the information on the serial dependence which is the most relevant for future prediction.
To achieve more effective dimension reduction and better predictive performance, we propose in this paper a two-step approach. Our first step is a segmentation transformation step in which we seek for a linear transformation ${\mathbf Y}_t(\cdot)={\mathbf A} {\mathbf Z}_t(\cdot)$, where ${\mathbf A}$ is a $p\times p$ invertible constant matrix, such that the transformed series ${\mathbf Z}_t(\cdot) = \{{\mathbf Z}_{t}^{(1)}(\cdot)^{ \mathrm{\scriptscriptstyle \top} }, \dots,{\mathbf Z}_{t}^{(q)}(\cdot)^{ \mathrm{\scriptscriptstyle \top} }\}^{ \mathrm{\scriptscriptstyle \top} }$ can be segmented into $q$ groups ${\mathbf Z}_{t}^{(1)}(\cdot), \dots,{\mathbf Z}_{t}^{(q)}(\cdot)$, and curve subseries ${\mathbf Z}_t^{(i)}(\cdot)$ and ${\mathbf Z}_t^{(j)}(\cdot)$ are uncorrelated at all time lags for any $i \neq j $, i.e., \[ {\rm Cov}\{ {\mathbf Z}_t^{(i)}(u), {\mathbf Z}_{t+k}^{(j)}(v) \} ={\mathbf 0}\,, \quad (u,v) \in {\mathcal U}^2 \;\; {\rm and}\;\; k=0, \pm 1, \pm 2, \ldots\,. \] Hence each ${\mathbf Z}_t^{(i)}$ can be modelled and forecasted separately as far as the linear dynamics is concerned. Under the stationarity assumption, the estimation of the transformation matrix ${\mathbf A}$ boils down to the eigenanalysis of a positive definite matrix defined by the double integral of quadratic forms in the autocovariance functions of ${\mathbf Y}_t(\cdot)$. An additional permutation on the components of ${\mathbf Z}_t(\cdot)$ will be specified in order to identify the latent group structure.
Our second step is to identify a finite-dimensional dynamic structure for each transformed subseries ${\mathbf Z}_{t}^{(l)}(\cdot)$ separately, which is based on a latent decomposition
where ${\mathbf X}^{(l)}_t(\cdot)$ represents the dynamics of ${\mathbf Z}^{(l)}_t(\cdot)$, $\boldsymbol{\varepsilon}_t^{(l)}(\cdot) $ is white noise with $\mathbb{E} \{ \boldsymbol{\varepsilon}^{(l)}_t(u) \} ={\mathbf 0}$ and $\mathbb{E}\{\boldsymbol{\varepsilon}^{(l)}_t(u)\boldsymbol{\varepsilon}^{(l)}_s(v)^{ \mathrm{\scriptscriptstyle \top} }\}={\bf0}$ for any $(u,v)\in {\mathcal U}^2$ and $t\ne s$, and $\{{\mathbf X}_t^{(l)}(\cdot)\}_{t=1}^n$ are uncorrelated with $\{\boldsymbol{\varepsilon}_t^{(l)}(\cdot)\}_{t=1}^n$. Furthermore we assume that the dynamic structure of ${\mathbf X}_t^{(l)}(\cdot)$ admits a vector time series presentation via a variational multivariate FPCA. For given $\{{\mathbf Z}_t^{(l)}(\cdot)\}_{t=1}^n$, the standard multivariate FPCA performs dimension reduction based on the eigenanlysis of the sample covariance function of ${\mathbf Z}_{t}^{(l)}(\cdot),$ which cannot be used to identify the finite-dimensional dynamic structure of ${\mathbf X}_{t}^{(l)}(\cdot)$ due to the contamination of $\boldsymbol{\varepsilon}^{(l)}_{t}(\cdot).$ Inspired by the fact that the lag-$k$ ($k \ne 0$) autocovariance function of ${\mathbf Z}_{t}^{(l)}(\cdot)$ automatically filters out the white noise, our variational multivariate FPCA is based on the eigenanalysis of a positive-definite matrix defined in terms of its nonzero lagged autocovariance functions; leading to a low-dimensional vector time series which bears all the dynamic structure of ${\mathbf X}_t^{(l)}(\cdot)$, and consequently, also that of ${\mathbf Z}_{t}^{(l)}(\cdot)$. This is possible as the number of components in each ${\mathbf Z}_{t}^{(l)}(\cdot)$ is usually small in practice. Finally, owing to the one-to-one linear transformation in the segmentation step, the good predictive performance of ${\mathbf Z}_t(\cdot)$ can be easily carried back to ${\mathbf Y}_t(\cdot).$
Our paper makes useful contributions on multiple fronts. Firstly, the segmentation transformation in the first step transforms the serial correlations across different series into the autocorrelations within each of the identified $q$ subseries. This not only avoids the direct modelling of the $p$ functional time series together, but also makes each of those transformed subseries more serially correlated and, hence, more predictable. As the serial correlations across different series are valuable for future prediction, the segmentation provides an effective way to use the information. Note that the prediction directly based on a multivariate ARMA-type model with even a moderately large dimension is not recommendable, as the gain from using the autocorrelations across different component series is often cancelled off by the errors in estimating too many parameters. Furthermore, even in the special case with $q = 1$, our decorrelation transformation can effectively push the cross-autocorrelations that are previously spread over $p$ components into a block-diagonally dominate structure, where the cross-autocorrelations along the block diagonal are significantly stronger than those off the diagonal. This still leads to reasonably good segmentation by retaining the strong within-group cross-autocorrelations while ignoring the weak between-group cross-autocorrelations and, as evidenced by simulations in Section (ref), results in more accurate future predictions than those based on models without transformation. Therefore, the proposed transformation can always be used as an initial step in modelling high-dimensional functional time series.
Secondly, though the segmentation transformation is motivated from the decorrelation idea of chang2018 for vector time series, its adaption to the functional setting introduces additional methodological and theoretical complexities and requires innovative advancements in both methodology and theory due to the intrinsic infinite-dimensionality of functional data. A simple extension of chang2018 would be to apply their method to the $p$-dimensional vector ${\mathbf Y}_t(u)$ on each evaluation grid value $u$ followed by aggregation, which fails to account for the smoothness and continuity of the functional nature of observed data. In contrast, our proposal on ${\mathbf Y}_t(\cdot)$ implements novel integral-based normalization and utilizes double integral over $(u,v) \in {\mathcal U}^2$ to fully leverage the autocovariance information, thus leading to more efficient estimation. Moreover, when performing permutation on the components of the transformed series, our method relies on Hilbert--Schmidt norm to measure the magnitude of bivariate functions, which introduces extra theoretical complexities compared to the absolute value measure used in chang2018. Finally, we develop a novel functional thresholding procedure, which guarantees the consistency of our estimation under high-dimensional scaling. Its theoretical analysis involves establishing novel inequalities between functional versions of matrix norms.
Thirdly, the nonzero lagged autocovariance-based dimension reduction approach in the second step makes the good use of the serial dependence information in our estimation, which is most relevant in time series prediction. On the method side, our proposed variational multivariate FPCA extends the univariate method of Bathia2010 by incorporating the cross-autocovariance. This extension addresses a crucial gap in dimension-reduction techniques, enabling us to accommodate multivariate functional time series. Importantly, when $p$ is fixed or moderately large, such method can be directly applied to the observed curve series ${\mathbf Y}_t(\cdot)$ for dimension reduction and forecasting purposes. On the theory side, we demonstrate that our proposal exhibits appealing convergence properties despite the additional transformation and estimation errors arisen from the first step, which are not involved in Bathia2010. By comparison, standard (multivariate) FPCA methods under ((ref)) suffer from inconsistent estimation and less efficient dimension reduction.
Existing research on functional time series has mainly focused on adapting the univariate or low-dimensional multivariate time series methods to the functional domain. An incomplete list of the relevant references includes Bathia2010, cho2013, aue2015, hormann2015, aue2018 and li2019. Following the recent emergence of high-dimensional functional time series data, there has been a wave of significant advancements aimed at addressing its complexities. Notable developments include functional factor models Gao2019,Tavakoli2023, functional dependence analysis guo2021, functional clustering tang2022, statistical inference for mean functions zhou2023, sparse vector functional autoregressions chang2024, and graphical PCA tan2024.
The rest of the paper is organized as follows. In Section (ref), we develop the methods employed in the first step, i.e. the segmentation transformation, the permutation and the functional thresholding. Section (ref) specifies the variational multivariate FPCA method used in the second dimension reduction step. We investigate the associated theoretical properties of the proposed methods in Section (ref). The finite-sample performance of our methods is examined through extensive simulations in Section (ref). Section (ref) applies our proposal to two real datasets, revealing its superior predictive performance over most frequently used competitors.
{\it Notation.} Denote by $I(\cdot)$ the indicator function. For a positive integer $m$, write $[m] = \{1,\dots, m\}$ and denote by ${\mathbf I}_m$ the identity matrix of size $m \times m$. For $x, y \in \mathbb{R},$ we use $x \vee y = \max(x,y).$ For two positive sequences $\{a_n\}$ and $\{b_n\}$, we write $a_n\ll b_n$ or $b_n\gg a_n$ if $\lim \sup_{n \to \infty} a_n/b_n = 0$. For a $p \times q$ real matrix ${\mathbf E}$, denote by ${\mathbf E}^{{ \mathrm{\scriptscriptstyle \top} }}$ its transpose, and write ${\mathbf E}^{\otimes2} = {\mathbf E} {\mathbf E}^{{ \mathrm{\scriptscriptstyle \top} }}$ and $\|{\mathbf E}\|_{2}= \lambda_{\max}^{1/2}({\mathbf E}^{{ \mathrm{\scriptscriptstyle \top} }}{\mathbf E}),$ where $\lambda_{\max}({\mathbf M})$ denotes the largest eigenvalue of the matrix ${\mathbf M}.$ Let $L_2({\mathcal U})$ be the Hilbert space of square integrable functions defined on ${\mathcal U}$ and equipped with the inner product $\langle f,g \rangle=\int_{{\mathcal U}} f(u)g(u)\,{\rm d}u$ for $f,g \in L_2({\mathcal U})$ and the induced norm $\|\cdot\|=\langle \cdot,\cdot \rangle^{1/2}.$ For any $B$ in ${\mathbb S}\equiv L_2({\mathcal U} \times {\mathcal U}), $ we denote the Hilbert--Schmidt norm by $\|B\|_{{\mathcal S}} = \{\int_{{\mathcal U}} \int_{{\mathcal U}} B^2(u,v)\, {\rm d}u {\rm d}v\}^{1/2}.$
We consider the following linear decomposition of ${\mathbf Y}_t(u)$:
where $q\in[p]$ is an unknown positive integer, ${\mathbf A}=({\mathbf A}_1,\ldots,{\mathbf A}_q)$ is a $p\times p$ unknown loading matrix, and ${\mathbf Z}_t(u)=\{{\mathbf Z}_{t}^{(1)}(u)^{ \mathrm{\scriptscriptstyle \top} },\ldots{\mathbf Z}_t^{(q)}(u)^{ \mathrm{\scriptscriptstyle \top} }\}^{ \mathrm{\scriptscriptstyle \top} }$ is a latent $p$-dimensional functional time series such that ${\rm Cov}\{{\mathbf Z}_{t}^{(l)}(u), {\mathbf Z}_{s}^{(l')}(v)\} = \boldsymbol{0}$ for all $t,s\in[n]$, $l \neq l'$ and $(u,v)\in {\mathcal U}^2$. Such linear decomposition possesses three key properties:
Our goal is then to find a linear decomposition (ref) for ${\mathbf Y}_t(\cdot)$, where each group ${\mathbf Z}_t^{(l)}(\cdot)$ for $l \in [q]$ cannot be further divided into smaller uncorrelated subgroups. This allows us to model each ${\mathbf Z}_t^{(l)}(\cdot)$ separately, as there are no cross-correlations among them at all time lags. We formalize the inseparability for each ${\mathbf Z}_t^{(l)}(\cdot)$ as Condition (ref) in Section (ref), which in turn defines the number of groups $q$ and the segmentation of ${\mathbf Z}_t(\cdot)$ in (ref). In Section (ref), we will present the estimation of the number of groups $q$, the linear spaces ${\mathcal C}({\mathbf A}_l)$ and the associated transformed subseries ${\mathbf Z}_t^{(l)}(\cdot)$ of group size $p_l$ for $l \in [q]$. Before that, let us firstly illustrate the validity and benefit of the linear decomposition (ref), i.e. segmentation transformation, in predicting multivariate functional time series with a real-life example. As we will demonstrate, such a decomposition (ref) is commonly achieved with a relatively large $q$ in practice. This effectively reduces the modelling burden while retaining the full linear dynamics of the original curve series $\mathbf{Y}_t(\cdot),$ thus leading to more accurate predictions.
We now consider how to find the segmentation transformation under (ref). Assume that $\max_{i \in [p]}\int_{{\mathcal U}} \mathbb{E}\{Z_{ti}^2(u)\}\,{\rm d}u = O(1)$. Define $\boldsymbol \Sigma_{y,k}(u, v) = {\rm Cov}\{{\mathbf Y}_{t}(u), {\mathbf Y}_{t+k}(v)\}$ and $\boldsymbol \Sigma_{z,k}(u, v) = {\rm Cov}\{{\mathbf Z}_{t}(u), {\mathbf Z}_{t+k}(v)\}.$ Without loss of generality, we focus on the orthogonal transformations only, i.e., ${\mathbf A}^{{ \mathrm{\scriptscriptstyle \top} }}{\mathbf A} = {\mathbf A}{\mathbf A}^{{ \mathrm{\scriptscriptstyle \top} }}={\mathbf I}_p$, as we can replace $({\mathbf Y}_t, {\mathbf Z}_t)$ in ((ref)) by $({\mathbf V}_y^{-1/2} {\mathbf Y}_t, {\mathbf V}_z^{-1/2} {\mathbf Z}_t)$ with ${\mathbf V}_y = \int_{{\mathcal U}}\boldsymbol \Sigma_{y,0} (u,u)\,{\rm d}u$ and ${\mathbf V}_z = \int_{{\mathcal U}}\boldsymbol \Sigma_{z,0} (u,u)\,{\rm d}u$. Then ${\mathbf A}$ is replaced by ${\mathbf V}_y^{-1/2} {\mathbf A} {\mathbf V}_z^{1/2}$ which is an orthogonal matrix as
Due to the unobservability of ${\mathbf Z}_t$, we can take ${\mathbf V}_z^{-1/2} {\mathbf Z}_t$ as ${\mathbf Z}_t$ since they share the same block structure. In practice, we can replace observations ${\mathbf Y}_t$ by $\widehat {\mathbf V}_y^{-1/2} {\mathbf Y}_t$, where $\widehat {\mathbf V}_y $ is a consistent estimator of ${\mathbf V}_y$.
For a given integer $k_0\ge 1,$ let
Then both ${\mathbf W}_y$ and ${\mathbf W}_x$ are non-negative definite. According to (ref), it holds that $ \boldsymbol \Sigma_{y,k}(u, v) ={\mathbf A} \boldsymbol \Sigma_{z,k}(u, v){\mathbf A}^{ \mathrm{\scriptscriptstyle \top} }\,,$ where $\boldsymbol \Sigma_{z,k}(u,v)$ is block-diagonal with blocks on the main diagonal of sizes $p_1\times p_1, \ldots, p_q\times p_q$. Due to ${\mathbf A} {\mathbf A}^{ \mathrm{\scriptscriptstyle \top} }= {\mathbf I}_p$, by ((ref)),
As all $\boldsymbol \Sigma_{z,k}(u, v)$ for $k\ge 0$ and $(u, v)\in {\mathcal U}^2$ are block-diagonal matrices of the same sizes, so is ${\mathbf W}_z$. Perform the eigenanalysis for each of $q$ blocks on the main diagonal of ${\mathbf W}_z$ separately, leading to $q$ orthogonal matrices of sizes $p_l \times p_l$ for $l\in [q]$. The columns of each of those orthogonal matrices are the $p_l$ orthonormal eigenvectors from the corresponding eigenanalysis. We form a $p\times p$ block diagonal orthogonal matrix $\boldsymbol{\Gamma}_z$ with those $q$ orthogonal matrices along the main block diagonal. Then the columns of $\boldsymbol{\Gamma}_z$ are the orthonormal eigenvectors of ${\mathbf W}_z$, i.e.,
where ${\mathbf D}$ is a diagonal matrix consisting of the $p$ eigenvalues. Then by ((ref)) and (ref), $ {\mathbf W}_y {\mathbf A} \boldsymbol{\Gamma}_z = {\mathbf A} {\mathbf W}_z \boldsymbol{\Gamma}_z = {\mathbf A} \boldsymbol{\Gamma}_z {\mathbf D}$. Thus the columns of $ \boldsymbol{\Gamma}_y\equiv{\mathbf A} \boldsymbol{\Gamma}_z$ are the orthonormal eigenvectors of ${\mathbf W}_y$. Combining this with ((ref)) yields that $ \boldsymbol{\Gamma}_y^{ \mathrm{\scriptscriptstyle \top} }{\mathbf Y}_t(\cdot) = \boldsymbol{\Gamma}_z^{ \mathrm{\scriptscriptstyle \top} }{\mathbf A}^{ \mathrm{\scriptscriptstyle \top} }{\mathbf Y}_t(\cdot) = \boldsymbol{\Gamma}_z^{ \mathrm{\scriptscriptstyle \top} }{\mathbf Z}_t(\cdot)\,. $ Since $\boldsymbol{\Gamma}_z$ is a block-diagonal orthogonal matrix with $q$ blocks, $\boldsymbol{\Gamma}_z^{ \mathrm{\scriptscriptstyle \top} }{\mathbf Z}_t(\cdot)$ effectively applies orthogonal transformation within each of the $q$ groups of ${\mathbf Z}_t(\cdot)$. Thus $\boldsymbol{\Gamma}_z^{ \mathrm{\scriptscriptstyle \top} }{\mathbf Z}_t(\cdot)$ is of the same segmentation structure of ${\mathbf Z}_t(\cdot)$, i.e. knowing $\boldsymbol{\Gamma}_z^{ \mathrm{\scriptscriptstyle \top} }{\mathbf Z}_t(\cdot)$ is as good as knowing the latent segmentation of ${\mathbf Z}_t(\cdot).$ By (ref), we have ${\mathbf Z}_t(\cdot) = {\mathbf A}^{{ \mathrm{\scriptscriptstyle \top} }} {\mathbf Y}_t(\cdot)$. Hence $\boldsymbol{\Gamma}_y$ can be taken as the required transformation matrix ${\mathbf A}$.
Let $\widehat \boldsymbol \Sigma_{y,k}(u,v)$ be some consistent estimator of $\boldsymbol \Sigma_{y,k}(u,v)$ for $k \in \{0\}\cup [k_0]$, to be specified in Section (ref) below. We define an estimator of ${\mathbf W}_y$ as
and calculate its orthonormal eigenvectors $\widehat \boldsymbol{\eta}_1,\dots,\widehat \boldsymbol{\eta}_p$. Let $\widehat \boldsymbol{\Gamma}_y=(\widehat{\boldsymbol{\eta}}_1,\ldots,\widehat{\boldsymbol{\eta}}_p)$. Then the required transformation matrix ${\mathbf A}$ can be estimated by a (latent) column-permutation of $\widehat \boldsymbol{\Gamma}_y$. More specifically, put
We propose below a data-driven procedure to divide the $p$ components of $\widehat{{\mathbf Z}}_t(\cdot)$ into $\hat q$ uncorrelated groups.
Recall $ {\mathbf Z}_t(\cdot)=\{Z_{t1}(\cdot),\ldots,Z_{tp}(\cdot)\}^{ \mathrm{\scriptscriptstyle \top} }$ with $\boldsymbol \Sigma_{z,k}(\cdot,\cdot)=\{\Sigma_{z,k,ij}(\cdot,\cdot)\}_{i,j\in[p]}.$ For two curve series $Z_{ti}(\cdot)$ and $ Z_{tj}(\cdot)$ within the same group, one would expect that their lag-$k$ cross-autocovariance function $\Sigma_{z,k,ij}(u,v)$ to be significantly different from zero for some integer $k$ and $(u,v)\in {\mathcal U}^2$, thus leading to at least one large $\|\Sigma_{z,k,ij}\|_{{\mathcal S}}$ for some integer $k$. Based on $\widehat{{\mathbf Z}}_t(\cdot)$ defined as (ref), we let $ \widehat{\boldsymbol \Sigma}_{z,k}(u,v)\equiv\{\widehat{\Sigma}_{z,k,ij}(u,v)\}_{i,j\in[p]}=\widehat \boldsymbol{\Gamma}_y^{ \mathrm{\scriptscriptstyle \top} } \widehat{\boldsymbol \Sigma}_{y,k}(u,v)\widehat \boldsymbol{\Gamma}_y$ for any $(u,v)\in{\mathcal U}^2$. Given a fixed integer $m\geq0$, we define the maximum cross-autocovariance over the lags between prespecified $-m$ and $m$ as
for any pair $(i,j)\in[p]^2$ such that $i < j$, and regard $\widehat Z_{ti}(\cdot)$ and $\widehat Z_{tj}(\cdot)$ from the same group if $\widehat T_{ij}$ takes some large value. To be specific, we rearrange $\aleph=p(p-1)/2$ values of $\widehat T_{ij}$ ($1 \leq i < j \leq p$) in the descending order $\widehat T_{(1)} \geq \cdots \geq \widehat T_{(\aleph)}$ and compute
for some $\delta_{n}>0$. Corresponding to $\widehat T_{(1)},\dots, \widehat T_{(\hat \varrho)},$ we identify $\hat \varrho$ pairs of cross-correlated curves. To divide the $p$ components of $\widehat{{\mathbf Z}}_t(\cdot)$ into several uncorrelated groups, we can first start with $p$ groups with each $\widehat Z_{tj}(\cdot)$ in one group and then repeatedly merge two groups if two cross-correlated curves are split over the two groups. The iteration is terminated until all the cross-correlated pairs are within one group. Hence we obtain the estimated group structure of $\widehat {\mathbf Z}_t(\cdot)$ with the number of the final groups $\hat q$ being the estimated value for $q$. Denote by $\widehat {\mathbf Z}_t^{(l)}(\cdot)$ the estimated $l$-th group for $l \in [\hat q]$. The estimated transformation matrix $ \widehat {\mathbf A} = (\widehat {\mathbf A}_1, \ldots, \widehat {\mathbf A}_{\hat q})$ can then be found by reorganizing the order of $(\widehat{\boldsymbol{\eta}}_1,\ldots,\widehat{\boldsymbol{\eta}}_p)$ such that
The estimate $\widehat{\boldsymbol \Sigma}_{y,k}(u,v)$ plays a key role in Section (ref). Let $\bar{{\mathbf Y}}(u)=n^{-1}\sum_{t=1}^n{\mathbf Y}_t(u).$ A natural candidate for $\widehat{\boldsymbol \Sigma}_{y,k}(u,v)$ is the sample version of $\boldsymbol \Sigma_{y,k}(u,v)$ defined as
When $p^2/n\rightarrow0$, $\widehat \boldsymbol \Sigma{}^{{\scriptscriptstyle\textup{S}}}_{y,k}(u,v)$ is a valid estimator for $\boldsymbol \Sigma_{y,k}(u,v)$. However when $p$ grows faster than $n^{1/2}$, it does not always hold that $\|\widehat \boldsymbol \Sigma{}^{{\scriptscriptstyle\textup{S}}}_{y,k}(u,v)-\boldsymbol \Sigma_{y,k}(u,v)\|_2 \to 0$ in probability. Under the high-dimensional scenario, the orthogonality of ${\mathbf A}$ naturally results in the magnitude of many of its entries being small, leading to certain sparsity on ${\mathbf A}$ which will then pass onto the autocovariance functions $\boldsymbol \Sigma_{y,k}(\cdot,\cdot)$, as $\boldsymbol \Sigma_{y,k}(\cdot, \cdot) = {\mathbf A} \boldsymbol \Sigma_{z,k}(\cdot, \cdot){\mathbf A}^{{ \mathrm{\scriptscriptstyle \top} }}.$
Inspired by the spirit of threshold estimator for large covariance matrix bickel2008, we apply the functional thresholding rule, which combines the functional generalizations of hard thresholding and shrinkage with the aid of the Hilbert--Schmidt norm of functions, on the entries of the sample autocovariance function $\widehat \boldsymbol \Sigma{}^{{\scriptscriptstyle\textup{S}}}_{y,k}(u,v)=\{\widehat\Sigma_{y,k,ij}^{{\scriptscriptstyle\textup{S}}}(u,v)\}_{i,j\in[p]} $ in ((ref)). This leads to the estimator
where $\omega_k \geq 0$ is the thresholding parameter at lag $k.$ Taking $\widehat \boldsymbol \Sigma_{y,k}$ in ((ref)) as ${\mathcal T}_{\omega_k}(\widehat \boldsymbol \Sigma{}^{{\scriptscriptstyle\textup{S}}}_{y,k})$ yields
Our second step is to represent (linear) dynamic structure of each ${\mathbf Z}_t^{(l)}(\cdot)$ in terms of a vector time series via representation ((ref)). The key idea is to identify the finite decomposition for ${\mathbf X}_t^{(l)}(\cdot)$. For $(u, v) \in {\mathcal U}^2$ and $k\ge 0$, let $\boldsymbol{\mu}^{(l)}(u)= \mathbb{E}\{{\mathbf X}_t^{(l)}(u)\}$ and
Then the multivariate Karhunen-Lo\`eve decomposition for ${\mathbf X}^{(l)}_t(\cdot)$ serving as the foundation of multivariate FPCA Chiou2014, Happ2018 admits the form
where $\lambda_1^{(l)} \ge \lambda_2 ^{(l)}\ge \dots \ge 0$ are the ordered eigenvalues of ${\mathbf M}_0^{(l)}(\cdot,\cdot)$, $\boldsymbol{\varphi}_1^{(l)}(\cdot), \boldsymbol{\varphi}_2^{(l)}(\cdot), \dots $ are the corresponding orthonormal eigenfunctions satisfying $\int_{{\mathcal U}} \boldsymbol{\varphi}_{j}^{(l)}(u)^{ \mathrm{\scriptscriptstyle \top} } \boldsymbol{\varphi}_{k}^{(l)}(u) \,{\rm d}u = I(j=k)$, and $ \xi_{tj}^{(l)} = \int_{\mathcal U} \boldsymbol{\varphi}_j^{(l)}(u) ^{ \mathrm{\scriptscriptstyle \top} }\{ {\mathbf X}_t^{(l)}(u) - \boldsymbol{\mu}^{(l)}(u) \}\,{\rm d}u $ with $\mathbb{E}\{\xi_{tj}^{(l)}\}=0$ and ${\rm Cov}\{\xi_{tj}^{(l)},\xi_{tk}^{(l)}\} = \lambda_j^{(l)}I(j=k).$
When ${\mathbf X}_t^{(l)}(\cdot)$ is $r_l$-dimensional in the sense that $\lambda_{r_l}^{(l)}>0$ and $\lambda_{r_l+1}^{(l)}=0,$ the dynamics of ${\mathbf X}^{(l)}_t(\cdot)$ is entirely determined by that of $r_l$-vector time series $\boldsymbol{\xi}_t^{(l)} = \{\xi_{t1}^{(l)}, \ldots, \xi_{tr_l}^{(l)}\}^{ \mathrm{\scriptscriptstyle \top} }$. Unfortunately, under the latent decomposition (ref), i.e.,
the standard multivariate FPCA based on ((ref)) is inappropriate as ${\mathbf X}_t^{(l)}(\cdot)$ is unobservable and we cannot provide a consistent estimator for ${\mathbf M}_0^{(l)}(u, v)$ based on ${\mathbf Z}_t^{(l)}(\cdot)$ due to the fact ${\rm Cov}\{{\mathbf Z}_{t}^{(l)}(u),{\mathbf Z}_{t}^{(l)}(v)\}={\mathbf M}_0^{(l)}(u, v) + {\rm Cov}\{\mbox{\boldmath$\varepsilon$}_t^{(l)}(u),\mbox{\boldmath$\varepsilon$}_t^{(l)}(v)\}.$
Now we introduce the variational multivariate FPCA based on a variational multivariate Karhunen-Lo\`eve decomposition for ${\mathbf X}_t^{(l)}(\cdot)$. Motivated from the fact ${\rm Cov}\{{\mathbf Z}_{t}^{(l)}(u),{\mathbf Z}_{t+k}^{(l)}(v)\}={\mathbf M}_k^{(l)}(u, v)$ for any $k \geq 1,$ for a prespecified small integer $k_0 \ge 1$, we define
Similar to ${\mathbf M}_0^{(l)}$, ${\mathbf K}^{(l)}$ is also non-negative definite and admits a spectral decomposition
where $\theta_1^{(l)} \ge \theta_2^{(l)} \ge \cdots \ge 0$ are the eigenvalues of ${\mathbf K}^{(l)}$, and $\boldsymbol{\psi}_1^{(l)}(\cdot), \boldsymbol{\psi}_2^{(l)}(\cdot), \ldots $ are the corresponding orthonormal eigenfunctions.
Proposition (ref) shows that, under the expansion ((ref)), ${\mathbf K}^{(l)}$ has exactly $r_l$ nonzero eigenvalues , and the dynamic space spanned by $\{\boldsymbol{\psi}_1^{(l)}(\cdot), \dots, \boldsymbol{\psi}_{r_l}^{(l)}(\cdot)\}$ remains the same as that spanned by $\{\boldsymbol{\varphi}_1^{(l)}(\cdot), \dots, \boldsymbol{\varphi}_{r_l}^{(l)}(\cdot)\}.$ Therefore, ${\mathbf X}_t^{(l)}(\cdot)$ can be expanded using $r_l$ basis functions $\boldsymbol{\psi}_1^{(l)}(\cdot), \dots, \boldsymbol{\psi}_{r_l}^{(l)}(\cdot),$ i.e.,
where the basis coefficients $\zeta^{({l})}_{tj} = \int_{\mathcal U} \boldsymbol{\psi}_j^{(l)}(u)^{ \mathrm{\scriptscriptstyle \top} } \{ {\mathbf X}_t^{(l)}(u) - \boldsymbol{\mu}^{(l)}(u) \}\, {\rm d}u$. Note that we take the sum in defining ${\mathbf K}^{(l)}(u,v)$ in ((ref)) to accumulate the information from different lags, and there is no information cancellation as each term in the sum is non-negative definite. An additional advantage for using the nonzero lagged autocovariance-based decomposition is that the identified directions $\boldsymbol{\psi}_1^{(l)}(\cdot),\ldots,\boldsymbol{\psi}_{r_l}^{(l)}(\cdot)$ catch the most significant serial dependence, which leads to the most efficient dimension reduction and is thus advantageous for prediction.
Noting that ${\mathbf Z}_t^{(l)}(\cdot)$ is not directly observable, we can only estimate ${\mathbf M}_k^{(l)}$ and ${\mathbf K}^{(l)}$ based on $\hat p_l$-vector of estimated transformed curve subseries $\widehat {\mathbf Z}_t^{(l)}(\cdot)= \{\widehat Z_{t1}^{(l)}(\cdot), \dots, \widehat Z_{t\hat p_l}^{(l)}(\cdot)\}^{{ \mathrm{\scriptscriptstyle \top} }}$ obtained in the segmentation transformation step. With the aid of (ref), for $k \in \{0\}\cup[k_0]$, put
It is easy to see from ((ref)) that $\widehat {\mathbf M}_k^{(l)}(u, v)$ is a reasonable estimator for ${\mathbf M}_k^{(l)}(u, v)$ when $k \geq1$, as it filters out white noise $\mbox{\boldmath$\varepsilon$}_t^{(l)}(\cdot)$ automatically. It is noteworthy that ((ref)) requires the consistent estimators for $\boldsymbol \Sigma_{y,k}(u,v)$. Its implementation under the high-dimensional setting can thus be done by setting $\widehat\boldsymbol \Sigma_{y,k}(u,v)={\mathcal T}_{\omega_k}(\widehat \boldsymbol \Sigma{}^{{\scriptscriptstyle\textup{S}}}_{y,k})(u,v)$ defined in (ref).
To estimate $\boldsymbol{\psi}_j^{(l)}(\cdot)$ and $\zeta^{(l)}_{tj}$ in ((ref)), we perform eigenanalysis of the estimator for ${\mathbf K}^{(l)},$
leading to the eigenvalues $\hat \theta_1^{(l)} \ge \hat \theta_2^{(l)} \ge \cdots \ge 0$, and the corresponding orthonormal eigenfunctions $\widehat \boldsymbol{\psi}_1^{(l)}(\cdot), \widehat\boldsymbol{\psi}_2^{(l)}(\cdot), \dots$. To estimate $r_l$ (i.e. the number of nonzero eigenvalues), we take the commonly-adopted ratio-based estimator for $r_l$ as:
for some $\tilde\delta_{n}>0$. Under some regularity conditions, such defined $\hat r_l$ is a consistent estimator for $r_l$; see Theorem (ref) in Section (ref). In practice, since $\tilde \delta_n$ is usually unknown, we instead adopt $ \hat r_l = \arg \max_{j \in [c_r(n-k_0)]} \hat \theta_j^{(l)} / \hat \theta_{j+1}^{(l)}$, where $c_r \in (0, 1)$ is a prescribed constant aiming to avoid fluctuations due to the ratios of extreme small values.
Let $\hat \zeta_{tj}^{(l)}=\int_{{\mathcal U}} \widehat \boldsymbol{\psi}_j^{(l)}(u)^{ \mathrm{\scriptscriptstyle \top} } \{\widehat {\mathbf Z}_t^{(l)}(u)- \widebar{\mathbf Z}^{(l)}(u)\}\, {\rm d}u$ for $t\in[n]$, $j\in[\hat r_l]$ and $l \in [\hat q]$. We can fit a model for the $\hat r_l$-dimensional vector time series $\widehat\boldsymbol{\zeta}_{t}^{(l)}=\{\hat \zeta_{t1}^{(l)}, \ldots, \hat \zeta_{t \hat r_l}^{(l)}\}^{ \mathrm{\scriptscriptstyle \top} }$ with $t\in[n]$ to obtain its $h$-step ahead prediction $\mathring \boldsymbol{\zeta}_{n+h}^{(l)}$ and then recover the $h$-step ahead functional prediction as
We finally obtain the $h$-step ahead prediction $ \widehat {\mathbf A} \mathring {\mathbf Z}_{n+h}(\cdot)$ for original functional time series, where $\widehat {\mathbf A} = (\widehat {\mathbf A}_1, \ldots, \widehat {\mathbf A}_{\hat q})$ and $\mathring {\mathbf Z}_{n+h} (\cdot)=\{\mathring{\mathbf Z}_{n+h}^{(1)} (\cdot)^{{ \mathrm{\scriptscriptstyle \top} }}, \dots, \mathring {\mathbf Z}_{n+h}^{(\hat q)}(\cdot)^{{ \mathrm{\scriptscriptstyle \top} }}\}^{{ \mathrm{\scriptscriptstyle \top} }}.$
This section presents theoretical analysis of our two-step estimation procedure. To ease presentation, we focus on the high-dimensional scenario and develop the theoretical results based on the estimator $ {\mathcal T}_{\omega_k}(\widehat \boldsymbol \Sigma{}^{{\scriptscriptstyle\textup{S}}}_{y,k})(u,v)$ in ((ref)). To simplify notation, we use $B $ to denote the linear operator induced from the kernel function $B \in {\mathbb S}, $ i.e., for any $f \in L_2({\mathcal U}),$ $B(f)(\cdot)=\int_{{\mathcal U}}B(\cdot\,,v)f(v)\,{\rm d}v \in L_2({\mathcal U}).$ Denote the $p$-fold Cartesian product $\mathbb{H}=L_2({\mathcal U}) \times \dots \times L_2({\mathcal U}).$ For any ${\boldsymbol f}, {\boldsymbol g}\in \mathbb{H},$ we denote the inner product by $ \mbox{$ \langle {\boldsymbol f}, {\boldsymbol g} \rangle $} =\int_{{\mathcal U}}{\boldsymbol f}(u)^{{ \mathrm{\scriptscriptstyle \top} }}{\boldsymbol g}(u)\,{\rm d}u $ with the induced norm $\|\cdot\| = \mbox{$ \langle \cdot, \cdot \rangle $} ^{1/2},$ and use ${\mathbf B} $ to denote the linear operator induced from the kernel matrix function ${\mathbf B} = (B_{ij})_{m_1 \times m_2}$ with each $B_{ij} \in {\mathbb S},$ i.e., for any ${\boldsymbol f} \in \mathbb{H},$ ${\mathbf B}({\boldsymbol f})(\cdot) = \int_{\mathcal U} {\mathbf B}(\cdot\,, v) {\boldsymbol f}(v)\,{\rm d} v \in \mathbb{H}$. We write $ \|{\mathbf B}\|_{{\mathcal S},\infty} = \max_{i\in [m_1]} \sum_{j = 1}^{m_2}\|B_{ij}\|_{{\mathcal S}}.$ Before imposing the regularity conditions, we firstly define the functional version of sub-Gaussianity that facilities the development of non-asymptotic results for Hilbert space-valued random elements.
Write $ \boldsymbol \Sigma{}_{y,k}(u,v)=\{\Sigma_{y,k,ij}(u,v)\}_{i,j\in[p]}$. Conditions (ref)(ii) and (ref)(iii) guarantee the covariance-stationarity of $\{{\mathbf Y}_t(\cdot)\}$ and imply that $\max_{j\in[p]}\rmintop\nolimits_{{\mathcal U}} \Sigma_{y,0,jj}(u,u)\,{\rm d}u = O(1).$ Condition (ref) places a finite upper bound on the functional stability measure, which characterizes the effect of small decaying eigenvalues of $\boldsymbol \Sigma_{y,0}$ on the numerator of ((ref)), thus being able to handle infinite-dimensional functional objects $Y_{tj}(\cdot).$ See its detailed discussion in guo2021. Conditions (ref) and (ref) are essential to derive $\max_{i,j\in[p]}\| {\widehat\Sigma{}^{{\scriptscriptstyle\textup{S}}}_{y,k,ij} - \Sigma_{y,k,ij}}\|_{\mathcal S} = O_{\rm p}\{\mathcal{M}_y (n^{-1}{\log p})^{1/2}\}$ for $\widehat\Sigma{}^{{\scriptscriptstyle\textup{S}}}_{y,k,ij}$ involved in (ref), which plays a crucial rule in our theoretical analysis.
The parameters $s_1$ and $s_2$ determine the row and column sparsity levels of ${\mathbf A},$ respectively. The row sparsity with small $s_1$ entails that each component of ${\mathbf Y}_t(\cdot)$ is a linear combination of a small number of components in ${\mathbf Z}_t(\cdot),$ while the column sparsity with small $s_2$ corresponds to the case that each $Z_{tj}(\cdot)$ has impact on only a few components of ${\mathbf Y}_t(\cdot).$ The parameter $\alpha$ also controls the sparsity level of ${\mathbf A}$ with a smaller value yielding a sparser ${\mathbf A}$. Write
Lemma (ref) in the supplementary material reveals that the functional sparsity structures in columns/rows of $\boldsymbol \Sigma_{y,k}(\cdot,\cdot)$ are determined by $s_1s_2p_{\dagger}$ with smaller values of $s_1$, $s_2$ and $p_{\dagger}$ yielding functional sparser $\boldsymbol \Sigma_{y,k}(\cdot,\cdot)$.
Recall that ${\mathbf W}_z=\textrm{diag}({\mathbf W}_{z,1},\ldots,{\mathbf W}_{z,q})$ in ((ref)) is a block-diagonal matrix, where ${\mathbf W}_{z,l}$ is a $p_l\times p_l$ matrix. We further define
where $\Lambda(\cdot)$ denotes the set of eigenvalues of the matrix, and assume $\rho>0$.
We first establish the group recovery consistency of the segmentation step. To do this, we reformulate the permutation in Section (ref) in an equivalent graph representation way. Recall $\boldsymbol{\Gamma}_y = {\mathbf A} \boldsymbol{\Gamma}_z$ and $\boldsymbol{\Gamma}_z$ is a block-diagonal orthogonal matrix with the main block sizes $p_1,\dots, p_q$. Write $ \boldsymbol{\Gamma}_z = \text{diag}(\boldsymbol{\Gamma}_{z,1}, \dots, \boldsymbol{\Gamma}_{z,q})$. Since ${\mathbf A} = ({\mathbf A}_1,\dots, {\mathbf A}_q),$ we have $\boldsymbol{\Gamma}_y\equiv(\boldsymbol{\eta}_1,\ldots,\boldsymbol{\eta}_p) = ({\mathbf A}_1\boldsymbol{\Gamma}_{z,1}, \ldots, {\mathbf A}_q\boldsymbol{\Gamma}_{z,q})$. The columns of $\boldsymbol{\Gamma}_y$ are naturally partitioned in to $q$ groups $G_1, \dots, G_q$, where $ G_l = \{\boldsymbol{\eta}_{\sum_{l'=0}^{l-1} p_{l'}+1}, \dots, \boldsymbol{\eta}_{\sum_{l'=0}^{l} p_{l'}}\}$ with $p_0=0$. To simplify the notation, we just write
Recall that the columns of such defined $\boldsymbol{\Gamma}_y$ are the eigenvectors of ${\mathbf W}_y$. For $\rho$ defined in (ref), if $\|\widehat{{\mathbf W}}_y-{\mathbf W}_y\|_2 \leq \rho/5$, by Lemma (ref) in the supplementary material, there exists an orthogonal matrix ${\mathbf H} = \text{diag} ({\mathbf H}_1, \dots, {\mathbf H}_q)$ with ${\mathbf H}_l\in\mathbb{R}^{p_l\times p_l}$ for each $l \in [q]$ and a column permutation matrix ${\mathbf R}$ for $\widehat \boldsymbol{\Gamma}_y$, such that $ \widehat \boldsymbol{\Gamma}_y {\mathbf R} \equiv (\widehat \boldsymbol{\Pi}_1, \dots, \widehat \boldsymbol{\Pi}_q)$ with $\widehat \boldsymbol{\Pi}_l \in \mathbb{R}^{p\times {p_l}}$, and
If the $p$ eigenvalues of ${\mathbf W}_y$ are distinct, ${\mathbf H}$ is a diagonal matrix with elements in the diagonal being $1$ or $-1$. Write $\boldsymbol{\Gamma}_y {\mathbf H} = ({\mathbf A}_1 \boldsymbol{\Gamma}_{z,1} {\mathbf H}_1, \dots, {\mathbf A}_q \boldsymbol{\Gamma}_{z,q} {\mathbf H}_q)\equiv (\boldsymbol{\gamma}_1,\dots, \boldsymbol{\gamma}_p)$. For each $l\in[q]$, we can define a graph $(G_l,E_l)$ such that $(i,j)\in E_l$ if and only if $\max_{|k| \leq m}\|\boldsymbol{\gamma}_{i}^{{ \mathrm{\scriptscriptstyle \top} }}\boldsymbol \Sigma_{y,k}\boldsymbol{\gamma}_{j}\|_{\mathcal S}\neq 0$.
Condition (ref) ensures that the group $G_l$ is inseparable at the minimal signal level $\varsigma$ given the transformation ${\mathbf A}_l \boldsymbol{\Gamma}_{z,l} {\mathbf H}_l$ for each $l\in[q]$, and facilitates the specifications of the true number of groups $q$ and the associated segmentation structure under (ref). Define $ T_{ij} =\max_{|k| \leq m}\|{\boldsymbol{\gamma}}_{i}^{ \mathrm{\scriptscriptstyle \top} } \boldsymbol \Sigma_{y,k}{\boldsymbol{\gamma}}_{j}\|_{\cal S}$ and $\varrho=\sum_{l=1}^q|E_l|$. Rearrange $\aleph=p(p-1)/2$ values of $ T_{ij}$ ($1 \leq i < j \leq p$) in the descending order, $ T_{(1)} \geq \cdots \geq T_{(\aleph)}. $ We then have $T_{(i)} \geq \varsigma$ for $i \in [\varrho]$ and $T_{(i)} = 0$ for $i \geq \varrho+1$. Denote by $E = \{(i,j): T_{ij} \geq T_{(\varrho)}, 1 \leq i < j \leq p\}$ the edge set of $G = [p]$ under the transformation $\boldsymbol{\Gamma}_y {\mathbf H}.$ The true segmentation $\{G_1, \dots, G_q\}$ in ((ref)) can then be identified by splitting $(G,E)$ into $q$ isolated subgraphs $(G_1, E_1),\ldots,(G_q,E_q)$, where $q$ represents the true number of uncorrelated groups.
Recall that with the aid of $\widehat \boldsymbol{\Gamma}_y$, the estimated segmentation is obtained via the ratio-based estimator $\hat \varrho$ as defined in ((ref)). To be specific, we build an estimated graph $(G, \widetilde E)$ with vertex set $G=[p]$ and edge set $ \widetilde E = \{(i,j): \widehat T_{ij} \geq \widehat T_{( \hat \varrho)}, 1 \leq i < j \leq p\} $, and split it into $\hat{q}$ isolated subgraphs $(\widetilde G_{1}, \widetilde E_{1}),\ldots,(\widetilde{G}_{\hat{q}},\widetilde{E}_{\hat{q}})$. Note that $p$ columns of $\widehat \boldsymbol{\Gamma}_y=(\widehat{\boldsymbol{\eta}}_1,\ldots,\widehat{\boldsymbol{\eta}}_p) $ correspond to the ordered eigenvalues $\lambda_1(\widehat {\mathbf W}_y)\geq\cdots\geq\lambda_p(\widehat {\mathbf W}_y)$. Write $ \widehat \boldsymbol{\Gamma}_y {\mathbf R} \equiv (\widehat \boldsymbol{\gamma}_1,\dots, \widehat \boldsymbol{\gamma}_p)$ and let $\pi:[p] \to [p]$ denote the permutation associated with ${\mathbf R}$, i.e., $\widehat \boldsymbol{\gamma}_i =\widehat \boldsymbol{\eta}_{\pi(i)}$. Based on the permutation mapping $\pi$, we let $ \widehat G_{l} = \{\pi^{-1}(i): i \in \widetilde G_{l}\}$ for $l\in [\hat q]$.
Theorem (ref) gives the group recovery consistency of our segmentation step. We next evaluate the errors in estimating ${\cal C}({\mathbf A}_l)$ for $l\in[q]$. Based on the estimated groups $\{\widehat G_1, \dots, \widehat G_{\hat q}\}$, we reorganize the order of $(\widehat{\boldsymbol{\gamma}}_1,\ldots,\widehat{\boldsymbol{\gamma}}_p) = (\widehat \boldsymbol{\eta}_{\pi(1)}, \dots, \widehat \boldsymbol{\eta}_{\pi(p)})$ and define $\widehat {\mathbf A}_l$ in (ref) as $\widehat {\mathbf A}_l = (\widehat \boldsymbol{\gamma}_i)_{i \in \widehat G_l}$ for $l \in [\hat q]$. We consider a general discrepancy measure chang2015,chang2018 between two linear spaces ${\mathcal C}({\mathbf E}_1)$ and ${\mathcal C}({\mathbf E}_2)$ spanned by the columns of ${\mathbf E}_1 \in \mathbb{R}^{p\times \tilde p_1}$ and ${\mathbf E}_2 \in \mathbb{R}^{p\times \tilde p_2}$, respectively, with ${\mathbf E}_i^{{ \mathrm{\scriptscriptstyle \top} }} {\mathbf E}_i ={\mathbf I}_{\tilde p_i}$ for $i \in[2]$ as
Then $ D\{{\mathcal C}( {\mathbf E}_1),{\mathcal C}( {\mathbf E}_2)\}$ is equal to 0 if and only if ${\mathcal C}({\mathbf E}_1) \subset {\mathcal C}({\mathbf E}_2)$ or ${\mathcal C}({\mathbf E}_2) \subset {\mathcal C}({\mathbf E}_1),$ and to 1 if and only if the two spaces are orthogonal.
Theorem (ref) presents the uniform convergence rate for $\min_{j\in[\hat q]} D\{{\mathcal C}( {\mathbf A}_l),{\mathcal C}( \widehat {\mathbf A}_j)\}$ over $l\in [q]$. The rate is faster for smaller values of $\{s_1, s_2, p_\dagger, {\cal M}_y,\alpha\},$ while enlarging the minimum eigen-gap between different blocks (i.e., larger $\rho$) reduces the difficulty of estimating each ${\mathcal C}( {\mathbf A}_l).$
Supported by Theorems (ref) and (ref), our subsequent theoretical results are developed by assuming that the group structure of ${\mathbf Z}_t(\cdot)$ is correctly identified or known, i.e., $\hat{q}=q$ and $\widehat{G}_l=G_l$ for each $l$. We now turn to investigate the theoretical properties of the dimension reduction step. Inherited from the segmentation step, ${\mathbf Z}_{t}^{(1)}(\cdot),\ldots,{\mathbf Z}_t^{(q)}(\cdot)$ rely on the specific form of ${\mathbf A} = ({\mathbf A}_1, \dots, {\mathbf A}_q),$ and thus is not uniquely defined. Yet intuitively, we only require a certain transformation matrix to make our subsequent analysis related to $\widehat \boldsymbol{\eta}_1,\ldots,\widehat \boldsymbol{\eta}_p$ mathematically tractable. Based on ((ref)), we define $\boldsymbol{\Pi}_{l} = {\mathbf A}_l \boldsymbol{\Gamma}_{z,l} {\mathbf H}_l$ and it holds that ${\mathcal C}(\boldsymbol{\Pi}_{l}) = {\mathcal C}({\mathbf A}_l)$ for each $l\in[q]$. Let ${\mathbf Z}_{t}^{(l)}(\cdot) =\boldsymbol{\Pi}_{l}^{{ \mathrm{\scriptscriptstyle \top} }}{\mathbf Y}_t(\cdot)$. Recall (ref) and (ref). The primary goal of the second dimension reduction step is to identify each $r_l$ and to estimate the associated dynamic space ${\mathcal C}_l =\text{span}\{\boldsymbol{\psi}_{1}^{(l)}(\cdot), \dots, \boldsymbol{\psi}_{r_l}^{(l)}(\cdot)\}$. Recall that $\{\hat \theta_{j}^{(l)},\widehat\boldsymbol{\psi}_{j}^{(l)}(\cdot)\}_{j \geq 1}$ are the eigenvalue/eigenfunction pairs of $\widehat {\mathbf K}^{(l)}(\cdot,\cdot)$ defined in ((ref)) with $\widehat {\mathbf A}_l = (\widehat \boldsymbol{\gamma}_i)_{i \in G_l}$ and the dimension $r_l$ is fixed for all $l \in [q]$. Our asymptotic results are based on the following regularity condition:
Theorem (ref) shows that $r_l$ can be correctly identified with probability tending to one uniformly over $l \in [q].$ Let $\widehat{{\mathcal C}}_l = \text{span}\{\widehat \boldsymbol{\psi}_{1}^{(l)}(\cdot),\dots, \widehat \boldsymbol{\psi}_{\hat r_l}^{(l)}(\cdot)\}$ be the dynamic space spanned by $\hat r_l$ estimated eigenfunctions. To measure the discrepancy between ${\mathcal C}_l$ and $\widehat {\mathcal C}_l,$ we introduce the following metric. For two subspaces ${\mathcal C}({\mathbf b}_1) = \text{span} \{{\mathbf b}_{11}(\cdot), \dots,{\mathbf b}_{1\tilde r_1}(\cdot) \}$ and ${\mathcal C}({\mathbf b}_2)= \text{span} \{{\mathbf b}_{21}(\cdot), \dots,{\mathbf b}_{2\tilde r_2}(\cdot) \}$ satisfying $\langle {\mathbf b}_{ij}, {\mathbf b}_{ik}\rangle=I(j=k)$ for each $i \in [2],$ the discrepancy measure between ${\mathcal C}({\mathbf b}_1)$ and ${\mathcal C}({\mathbf b}_2)$ is defined as $$ \widetilde D\{{\mathcal C}({\mathbf b}_1),{\mathcal C}({\mathbf b}_2)\} = \sqrt{1-\frac{1}{\max(\tilde r_1,\tilde r_2)}\sum_{j=1}^{\tilde r_1} \sum_{k=1}^{\tilde r_2}\langle {\mathbf b}_{1j}, {\mathbf b}_{2k}\rangle^2}\, \in [0,1]\,, $$ which equals 0 if and only if ${\mathcal C}({\mathbf b}_1)\subset{\mathcal C}({\mathbf b}_2)$ or ${\mathcal C}({\mathbf b}_2)\subset{\mathcal C}({\mathbf b}_1)$ and 1 if and only if two spaces are orthogonal.
We conduct a series of simulations to illustrate the finite sample performance of the proposed methods. To simplify the data-generating process, we consider a relaxed form of ((ref)) as
with no orthonormality restriction on the transformation matrix $\widecheck {\mathbf A} = (\widecheck {\mathbf A}_1,\dots,\widecheck {\mathbf A}_q)$. The $p$-dimensional transformed functional time series $\widecheck {\mathbf Z}_t(\cdot)$ is formed by $q$ uncorrelated groups $\{\widecheck {\mathbf Z}_t^{(l)}(\cdot): l \in [q]\}$, where each $\widecheck {\mathbf Z}_t^{(l)}(\cdot)$ arises as the sum of dynamics $ \widecheck{\mathbf X}_{t}^{(l)}(\cdot)$ and white noise $ \check\boldsymbol \varepsilon_{t}^{(l)}(\cdot).$ Based on (ref) in Section (ref), ((ref)) can then be easily reformulated as ((ref)) by setting
where ${\mathbf V}_{\check y} = \int_{{\mathcal U}}{\rm Cov}\{\widecheck{\mathbf Y}_{t}(u), \widecheck{\mathbf Y}_{t} (u)\}\,{\rm d}u$ and ${\mathbf V}_{\check z} = \int_{{\mathcal U}}{\rm Cov}\{\widecheck{\mathbf Z}_{t}(u), \widecheck{\mathbf Z}_{t} (u)\}\,{\rm d}u$. Then the orthonormality of ${\mathbf A}$ is satisfied.
Write $\check\boldsymbol \varepsilon_{t}(\cdot)= \{\check\boldsymbol \varepsilon_{t}^{(1)}(\cdot)^{ \mathrm{\scriptscriptstyle \top} }, \dots,\check\boldsymbol \varepsilon_{t}^{(q)}(\cdot)^{ \mathrm{\scriptscriptstyle \top} } \}^{ \mathrm{\scriptscriptstyle \top} } \equiv \{\check \varepsilon_{t1}(\cdot), \dots, \check \varepsilon_{tp}(\cdot)\}^{{ \mathrm{\scriptscriptstyle \top} }}.$ We generate each curve component of $\check\boldsymbol{\varepsilon}_t(\cdot)$ independently by $\check\varepsilon_{tj}(\cdot) = \sum_{l=1}^{10} {2^{-(l-1)}}e_{tjl} \psi_{l}(\cdot)$ for $j \in [p],$ where $e_{tjl}$'s are sampled independently from $\mathcal{N}(0,1)$ and $\{\psi_l(\cdot)\}_{l=1}^{10}$ is a 10-dimensional Fourier basis function. The finite-dimensional dynamics $\widecheck{\mathbf X}_t(\cdot)= \{\widecheck{\mathbf X}_{t}^{(1)}(\cdot)^{ \mathrm{\scriptscriptstyle \top} }, \dots,\widecheck{\mathbf X}_{t}^{(q)}(\cdot)^{ \mathrm{\scriptscriptstyle \top} } \}^{ \mathrm{\scriptscriptstyle \top} }$ with prescribed group structure is generated based on some $5$-dimensional curve dynamics $\vartheta_{tg}(\cdot) = \sum_{l = 1}^5\kappa_{tgl} \psi_{l}(\cdot)$ for $g\in[20]$. The basis coefficients $\boldsymbol{\kappa}_{tg}=(\kappa_{tg1}, \dots, \kappa_{tg5})^{{ \mathrm{\scriptscriptstyle \top} }}$ are generated from a stationary VAR model $\boldsymbol{\kappa}_{tg}={\mathbf U}_g\boldsymbol{\kappa}_{(t-1)g}+ {\mathbf e}_t$ for each $g$. To guarantee the stationarity of $\boldsymbol{\kappa}_{tg},$ we generate ${\mathbf U}_g = \iota\widecheck {\mathbf U}_g/\rho(\widecheck {\mathbf U}_g)$ with $\iota \sim$ $\text{Uniform}[0.5,1]$ and $\rho(\widecheck {\mathbf U}_g)$ being the spectral radius of $\widecheck {\mathbf U}_g \in {\mathbb R}^{5 \times 5},$ the entries of which are sampled independently from $\text{Uniform}[-3,3]$. The components of the innovation ${\mathbf e}_t$ are sampled independently from ${\cal N}(0,1).$ We will specify the exact forms of $\widecheck{\mathbf X}_t(\cdot)$ under the fixed and large $p$ scenarios in Sections (ref) and (ref), respectively. The white noise sequence $\check \boldsymbol{\varepsilon}_t(\cdot)$ ensures that $\widecheck {\mathbf Z}_t(\cdot)$ as well as ${\mathbf Z}_t(\cdot)$ share the same group structure as $\widecheck {\mathbf X}_t(\cdot).$ Unless otherwise stated, we set $k_0 = m = 5$ and $c_r = c_\varrho =0.75$ in our procedure, as our simulation results suggest that our procedure is robust to the choices of these parameters.
We consider the following three examples of $\widecheck{\mathbf X}_t(\cdot)=\{\check X_{t1}(\cdot), \dots, \check X_{tp}(\cdot)\}^{{ \mathrm{\scriptscriptstyle \top} }}$ with different group structures for $p\in\{6, 10, 15\}$ based on independent $\vartheta_{t1}(\cdot), \dots, \vartheta_{t5}(\cdot).$
\newcounter{bean} \setcounter{bean}{0}
Therefore, $\widecheck {\mathbf X}_t(\cdot)$ consists of $q=3, 4$ and $5$ uncorrelated groups of curve subseries in Examples (ref), (ref) and (ref), respectively, where the number of component curves per group is $p_l=l$ for $l\in [q].$ The $p$-dimensional observed functional time series $\widecheck {\mathbf Y}_t(\cdot) = \{\check Y_{t1}(\cdot), \dots, \check Y_{tp}(\cdot)\}^{{ \mathrm{\scriptscriptstyle \top} }}$ for $t\in [n]$ is then generated by ((ref)) with the entries of $\widecheck {\mathbf A}$ sampled independently from $\text{Uniform}[-3,3]$. To obtain $h$-step ahead prediction of $\widecheck {\mathbf Y}_t(\cdot),$ we integrate the segmentation and dimension reduction steps respectively in Sections (ref) and (ref) into the VAR estimation as outlined in Algorithm (ref). For each of the three examples introduced above, we select
with $ \bar {\mathbf Y}_*(\cdot) = (n-h-k)^{-1}\sum_{t=1}^{n-h-k} \widetilde {\mathbf Y}_{t}(\cdot)$, for the quantities involved in Step (ref) of Algorithm (ref). We refer to the segmentation-(Variational-multivariate-FPCA)-and-VAR-based Algorithm (ref) with selections of $\widehat {\mathbf V}_{\check y}^{(h)}$ in ((ref)) and $\widehat \boldsymbol \Sigma_{\tilde{y},k}^{(h)}(u,v) $ in ((ref)) as SegV hereafter.
The performance of our two-step proposal is examined in terms of linear space estimation, group identification and post-sample prediction. For ${\mathbf A}=({\mathbf A}_1,\dots, {\mathbf A}_q)$ specified in ((ref)), with the aid of (ref), define $ f(l) = \arg\min_{j\in[\hat q]} D^2\{{\mathcal C}( {\mathbf A}_l),{\mathcal C}( \widehat {\mathbf A}_j)\}$ for each $l\in[q]$. We then call $\widehat {\mathbf A} = (\widehat {\mathbf A}_1,\dots, \widehat {\mathbf A}_{\hat q})$ an effective segmentation of ${\mathbf A}$ if (i) $1<\hat q \leq q$, and (ii) $\text{rank}(\widehat {\mathbf A}_{l'})=\sum_{l\in[q]:\,f(l) = l' }\text{rank}({\mathbf A}_l) $ for each $l' \in [\hat{q}]$. The intuition is as follows. The effective segmentation implies that each identified group in $\widehat {\mathbf Z}_t(\cdot)$ contains at least one, but not all, groups in $ {\mathbf Z}_t(\cdot).$ Since our main target is to forecast $\widecheck {\mathbf Y}_t(\cdot)$ based on the cross-serial dependence in $\{{\mathbf Z}_t^{(l)}(\cdot):l \in [q]\},$ this segmentation result is effective in the sense that the linear dynamics in ${\mathbf Z}_t(\cdot)$ is well kept in $\{\widehat {\mathbf Z}_t^{(l)}(\cdot): l \in [\hat q]\}$ without any contamination or damage and a mild dimension reduction is achieved with $\hat q>1.$ For the special case of complete segmentation ($\hat q = q$), we use the maximum and averaged estimation errors for $(\widehat {\mathbf A}_1,\dots, \widehat {\mathbf A}_{\hat q})$, respectively, defined as $ {\rm MaxE} = \max_{l \in [q]} D^2\{{\mathcal C}( {\mathbf A}_l),{\mathcal C}( \widehat {\mathbf A}_{f(l)})\}$ and ${\rm AvgE} = q^{-1}\sum_{l = 1}^q D^2\{{\mathcal C}( {\mathbf A}_l),{\mathcal C}( \widehat {\mathbf A}_{f(l)})\} $ to assess the ability of our method in fully recovering the spanned spaces ${\mathcal C}({\mathbf A}_1), \dots, {\mathcal C}({\mathbf A}_q).$ Note that ${\mathbf A}$ in ((ref)) can not be easily computed, as the true ${\mathbf V}_{\check y}$ and ${\mathbf V}_{\check z}$ are hard to find even for simulated examples. For $\widecheck {\mathbf A} $ specified in ((ref)), let $\widetilde {\mathbf A} ={\mathbf V}_{\check y}^{-1/2} \widecheck{\mathbf A} \equiv (\widetilde{\mathbf A}_1,\dots, \widetilde{\mathbf A}_q)$ with $\widetilde {\mathbf A}_l = {\mathbf V}_{\check y}^{-1/2} \widecheck{\mathbf A}_l$. Since ${\mathbf V}_{\check z}$ is a block-diagonal matrix, then ${\mathcal C}(\widetilde {\mathbf A}_l) = {\mathcal C}({\mathbf A}_l)$ for $l \in [q]$. Hence, we can replace ${\mathcal C}({\mathbf A}_l)$ by ${\mathcal C}( \{\widehat{\mathbf V}_{\check y}^{(h)}\}^{-1/2}\widecheck {\mathbf A}_l)$ to obtain the approximations of MaxE and AvgE in our simulations.
To evaluate the post-sample predictive accuracy, we define the mean squared prediction error (MSPE) as
with $v_1, \dots, v_N$ being equally-spaced points in $[0,1],$ and compute the relative prediction error as the ratio of MSPE in ((ref)) to that under the `oracle' case. In the oracle case, we apply the procedure in Section (ref) directly on each true $\{\widecheck {\mathbf Z}_{t}^{(l)}(\cdot)\}_{t=1}^{n-h}$ to achieve the $h$-step ahead prediction for $\{\widecheck {\mathbf Z}_{n}^{(l)}(\cdot) :l \in [ q]\},$ denoted by $\{\breve {\mathbf Z}_{n}^{(l)}(\cdot) :l \in [ q]\},$ and further obtain the $h$-step ahead prediction $\widecheck {\mathbf A} \{\breve {\mathbf Z}_{n}^{(1)} (\cdot)^{{ \mathrm{\scriptscriptstyle \top} }}, \dots, \breve {\mathbf Z}_{n}^{( q)}(\cdot)^{{ \mathrm{\scriptscriptstyle \top} }}\}^{{ \mathrm{\scriptscriptstyle \top} }}$ for the original curves $\widecheck {\mathbf Y}_t(\cdot).$ By comparison, we also implement an {\bf uni}variate functional prediction method on each $\check Y_{tj}(\cdot)$ separately by performing {uni}variate dimension reduction Bathia2010, then predicting vector time series based on the best fitted {\bf V}AR model and finally recovering functional prediction (denoted as UniV).
We generate $n\in \{200, 400, 800, 1600\}$ observations with $N=30$ for each example and replicate each simulation 500 times. Table (ref) provides numerical summaries, including the relative frequencies of the effective segmentation with $\hat q=q$ and $\hat q \geq q-1,$ and the estimation errors for $\widehat {\mathbf A}=(\widehat {\mathbf A}_1,\dots, \widehat {\mathbf A}_{\hat q})$ under the complete segmentation case. As one would expect, the proposed method provides higher proportions of effective segmentation and lower estimation errors as $n$ increases, and performs fairly well for reasonably large $n$ as $p$ increases. For $(p,n)=(6,200),$ we observe $62.6 \%$ complete segmentation with ${\rm AvgE}$ as low as $0.079.$ Furthermore, the proportions of effective segmentation with $\hat q \geq q-1$ are above $93\%$ for $n \geq 200.$ Similar results can be found for cases of $(p,n) = (10,800+)$ and $(15, 1600),$ whose proportions of effective segmentation with $\hat q \geq q-1$ remain higher than $87.4\%$ and $83.2\%,$ respectively. Table (ref) also reports the relative one-step ahead prediction errors. It is evident that SegV significantly outperforms UniV in all settings, demonstrating the effectiveness of our proposed segmentation transformation and dimension reduction in predicting future values. Although the proportions of complete segmentation are not high when $p=15$, the corresponding proportions of $\hat q \geq q-1$ become satisfactorily higher, and SegV performs similarly to the oracle case with its relative prediction errors being closer to 1 as $n$ increases.
Under a large $p$ scenario, a natural question to ask is whether the segmentation method based on the classical estimation for autocovariance functions of $\widetilde {\mathbf Y}_t(\cdot)$ (denoted as NonT) as (ref) in Section (ref) still performs well, and if not, whether a satisfactory improvement is attainable via the functional-thresholding estimation (denoted as FunT) developed in Section (ref). To this end, we generate $\widecheck {\mathbf Y}_t(\cdot)$ from ((ref)) with $p\in\{30, 60\}$ and $n\in\{200, 400\}$. Specifically, we let $\check X_{t(3l-2)}(\cdot)=\vartheta_{tl}(\cdot),$ $\check X_{t(3l-1)}(\cdot)=\vartheta_{(t+1)l}(\cdot),$ $\check X_{t(3l)}(\cdot)=\vartheta_{(t+2)l}(\cdot)$ for $l\in [q].$ This setting ensures $q$ uncorrelated groups of curve subseries in $\widecheck {\mathbf X}_t(\cdot)$ with $p_l=3$ component curves per group and hence $q=10$ and $20$ correspond to $p=30$ and $60,$ respectively. Let the $p \times p$ transformation matrix $\widecheck {\mathbf A} = \boldsymbol{\Delta}_1 + \delta \boldsymbol{\Delta}_2.$ Here $\boldsymbol{\Delta}_1 = \text{diag} \{\boldsymbol{\Delta}_{11},\dots, \boldsymbol{\Delta}_{1 (p/6)}\}$ with elements of each $\boldsymbol{\Delta}_{1i} \in \mathbb{R}^{6 \times 6}$ being sampled independently from $\text{Uniform}[-3,3]$ for $i \in [p/6],$ and $\boldsymbol{\Delta}_2$ is a matrix with two randomly selected nonzero elements from $\text{Uniform}[-1,1]$ each row. We set $\delta\in\{0.1, 0.5\}$. It is notable that our setting results in a very high-dimensional learning task in the sense that the intrinsic dimension $30\times 5=150$ or $60 \times 5=300$ is large relative to the sample size $n=200$ or $400.$
We assess the performance of NonT and FunT in discovering the group structure. The optimal thresholding parameters $\hat\omega_k$ in FunT are selected by the five-fold cross-validation (see Remark (ref)), and $ {\mathbf V}_{\check y}$ in the normalization step is estimated by $\widehat {\mathbf V}_{\check y}^{(0)}$ given in ((ref)), as the threshold version of $\widehat {\mathbf V}_{\check y}^{(0)}$ might not be positive definite. In practice, when $p$ is large, FunT may lead to segmentation with a small $\hat q,$ indicating that some groups of $\{\widehat {\mathbf Z}_{t}^{(l)}(\cdot): l \in [\hat q]\}$ contain multiple groups in $\{ {\mathbf Z}_{t}^{(l)}(\cdot): l \in [ q]\}.$ To ease the modelling burden of complex VAR process, we may consider performing further segmentation transformation on the estimated groups by repeating FunT $R$ times. To be precise, the $i$-th round of segmentation transformation via FunT is performed within each group discovered in the $(i-1)$-th round with $c_\varrho = 1$ for $i \in [R],$ and hence $(\widehat {\mathbf A}_1, \dots, \widehat {\mathbf A}_{\hat q})$ is updated after each iteration. Table (ref) reports the relative frequencies of the effective segmentation for NonT and FunT with $R\in\{1,5,10\}$. Finally, we apply FunT-based SegV (denoted as FTSegV) combined with the $R$-round segmentation transformation for $R\in\{1,5,10\}$ in Step (ref) of Algorithm (ref), and compare their one-step ahead predictive performance with UniV and SegV. Table (ref) summarizes the relative prediction errors for all five comparison methods.
Several conclusions can be drawn from Tables (ref) and (ref). Firstly, the performance of SegV severely deteriorates under the high-dimensional setting, as this procedure fails to detect any effective segmentation, resulting in elevated prediction errors. By comparison, FTSegV exhibits superior predictive ability over SegV and UniV. In particular, for large $n,$ e.g., $n = 400$, FTSegV does a reasonably good job in recovering the group structure of ${\mathbf Z}_t(\cdot)$ and performs comparably well to the oracle method with the relative prediction errors lower than $1.149$ in all scenarios. Secondly, comparing the results for $n = 200$ among different $R,$ we observe an interesting phenomenon that even though the relative frequencies of effective segmentation for FunT drop as $R$ increases, implying that some groups in $\{\widehat {\mathbf Z}_{t}^{(l)}(\cdot): l \in [\hat q]\}$ are split incorrectly before forecasting, the prediction errors stay low and slightly decrease as shown in Table (ref). This is not surprising, since further segmentation based on FunT yields fewer parameters to be estimated in VAR models and thus benefits the forecasting accuracy even if a few small but significant cross-covariances of ${\mathbf Z}_t(\cdot)$ are ignored. Such finding highlights the success of FTSegV and its $R$-round segmentation in the sense that although FTSegV may not be able to accurately recover the group structure in ${\mathbf Z}_t(\cdot)$ for a small $n$, it achieves an appropriate dimension reduction to provide significant improvement in high-dimensional functional prediction.
To further illustrate the advantage of our proposed segmentation transformation in predicting high-dimensional functional time series, we simulate data from a more generalized functional time series framework instead of strictly adhering to (ref). Specifically, we consider the vector functional autoregressive (VFAR) model of order 1,
where $\boldsymbol \epsilon_t(\cdot) = \{\epsilon_{t1}(\cdot), \dots, \epsilon_{tp}(\cdot) \}^{{ \mathrm{\scriptscriptstyle \top} }}$ are independently sampled from a $p$-dimensional vector of mean zero Gaussian processes, independent of ${\mathbf Y}_{t-1}(\cdot)$, and ${\mathbf Q} = (Q_{ij})_{i,j \in [p]}$ is the functional transition matrix with each $Q_{ij} \in \mathbb{S}$. See Section (ref) of the supplementary material for the detailed data-generating process.
We compare the predictive performance of three competing methods. The first VFAR method is developed by knowing the true data-generating process through VFAR model. We relegate the detailed prediction procedure to Section (ref) of the supplementary material. We next consider two segmentation-based prediction methods:
Table (ref) reports one-step ahead MSPEs for three methods with different values of $p.$ As anticipated, the performance of VFAR deteriorates severely as $p$ increases, demonstrating that the joint model suffers from the high-dimensionality, even when the true model is known. Meanwhile, both segmentation-based prediction methods exhibit improved predictive performance, with Seg+Z notably outperforming Seg+Y, particularly in scenarios with large $p.$ It is crucial to emphasize that the improvement of Seg+Z over Seg+Y is attributed to the decorrelation transformation. Table (ref) in the supplementary material provides further insights into the impact of transformation, where $\mathring q_y$ and $\mathring q_z$ denote the numbers of the identified groups using Seg+Y and Seg+Z, respectively. Interestingly, Seg+Z yields more groups than Seg+Y while retaining the same amount of strongly connected pairs. This observation indicates that the decorrelation transformation effectively pushes the cross-autocorrelations that were previously spread over $p$ components into a block-diagonally dominate structure, where the cross-autocorrelations along the block diagonal are significantly stronger than those off the diagonal. Such enhancement of within-group autocorrelations, along with the reduction of cross-autocorrelations between the groups, leads to reasonably good segmentation by only retaining the strong within-group cross-autocorrelations while ignoring the weak between-group cross-autocorrelations, and thus yields more accurate future predictions.
In this section, we apply our proposed SegV and FTSegV to two real data examples arising from different fields. Our main goal is to evaluate the post-sample predictive accuracy of both methods. By comparison, we also implement componentwise univariate prediction method (UniV) and the multivariate prediction method of Gao2019 (denoted as GSY) to jointly predict $p$ component series by fitting a factor model to estimated scores obtained via eigenanalysis of the long-run covariance function hormann2015. It is worth mentioning that the joint prediction model VmV (see Example (ref)) completely fail due to high dimensionality, so we do not report their results here. To evaluate the effectiveness of the segmentation transformation and its impact on prediction, we forge two other segmentation cases, namely under-segmentation and uni-segmentation, for both SegV and FTSegV (denoted as Under.SegV, Uni.SegV, Under.FTSegV and Uni.FTSegV, respectively). Denote by $\{\widehat G_l:l \in [\hat q]\}$ the segmented groups of $\{\widehat {\mathbf Z}_{t}^{(l)}(\cdot): l \in [\hat q]\}$ discovered in Step (ref) of Algorithm (ref) (seen also as correct-segmentation). The under-segmentation updates $\{\widehat G_l:l \in [\hat q]\}$ by merging two groups $\widehat{G}_{l_1}$ and $\widehat{G}_{l_1'}$ together before subsequent analysis, where $\arg\max_{(i,j):\,i\in\widehat{G}_l,\,j\in\widehat{G}_{l'},\, 1 \leq l\neq l' \leq \hat q}\widehat{T}_{ij} \in \widehat{G}_{l_1}\times \widehat{G}_{l_1'}$ with $\widehat{T}_{ij}$ defined in ((ref)). The uni-segmentation, on the other hand, regards each curve component of $\{\widehat {\mathbf Z}_{t}^{(l)}(\cdot): l \in [\hat q]\}$ as an individual group and then applies UniV componentwisely. For a fair comparison, the orders of VAR models adopted in all SegV/FTSegV-related methods and UniV are determined by the AIC criterion, while GSY is implemented using the R package \verb"ftsa".
To examine the predictive performance, we apply an expanding window approach to the observed data $\check Y_{tj}(v_i)$ for $t\in[n], j\in[p], i\in[N]$. We first split the dataset into a training set and a test set respectively consisting of the first $n_1$ and the remaining $n_2$ observations. For any positive integer $h,$ we implement each comparison method on the training set $\{\check Y_{tj}(v_i): t\in[n_1], j\in[p], i\in[N]\}$ and obtain its $h$-step ahead prediction, denoted as $\hat Y_{(n_1+h)j}^{(h)}(v_i)$, based on the fitted model. We then increase the training size by one, i.e. $\{\check Y_{tj}(v_i):t\in[n_1+1], j\in[p], i\in[N]\},$ refit the model and compute the next $h$-step ahead prediction $\hat Y_{(n_1+1+h)j}^{(h)}(v_i)$ for $ j\in[p], i\in[N].$ Repeat the above procedure until the last $h$-step ahead prediction $\hat Y_{nj}^{(h)}(v_i)$ is produced. Finally, we compute the $h$-step ahead MAPE and MSPE as
The first dataset, analyzed in tang2022, contains age-specific and gender-specific mortality rates for developed countries during 1965 to 2013 ($n=49$). See Table (ref) in the supplementary material for the list of $p = 29$ countries after removing certain countries with missing data. Following the proposal of tang2022, we model the log transformation of the mortality rate of people aged $v_i=i-1$ living in the $j$-th country during year $1964+t$ as a random curve $\check Y_{tj}(v_i)$ ($t\in[49],$ $j\in [29],$ $i \in [101]$) and perform smoothing for observed mortality curves via smoothing splines. We divide the smoothed dataset into the training set of size $n_1=34$ and the test set of size $n_2=15.$ Since the smoothed curve series exhibit weak autocorrelations when lags are beyond 3 and the training size is relatively small, we use $k_0=m=3$ in our procedure for this example.
Table (ref) reports the MAPEs and MSPEs for females and males. Several obvious patterns are observable. Firstly, our proposed methods, SegV and FTSegV, provide the best predictive performance uniformly for both females and males, and all $h$. This demonstrates the effectiveness of reducing the number of parameters via the segmentation transformation in predicting high-dimensional functional time series. Secondly, although the cases of under- and uni-segmentation are inferior to the correct-segmentation case, they significantly outperform UniV and GSY. Note that the improvement of Uni.SegV over UniV reveals the capability of the transformation matrix $\widehat {\mathbf A}$ to effectively decorrelate the original curves, thereby leading to more accurate predictions. One may also notice that, Uni.SegV does not perform as well as SegV and Under.SegV. In most cases, the transformed curve series exhibits $\hat{q}=26$ groups, with $25$ groups of size 1 and one large group of size 4; see Figures (ref)--(ref) in the supplementary material. The limitation of Uni.SegV thus becomes apparent as it fails to account for the cross-serial dependence within the large group, resulting in less accurate predictions. This finding again confirms the effectiveness of our procedure, in particular, the within-group cross-autocorrelations is also valuable in forecasting future values.
Our second dataset contains energy consumption readings (in kWh) taken at half hourly intervals for thousands of London households, and is available at \url{https://data.london.gov.uk/dataset/smartmeter-energy-use-data-in-london-households}. In our study, we select households with flat energy prices during the period between December 2012 and May 2013 ($n=182$) after removing samples with too many missing records, and hence construct $4000$ samples of daily energy consumption curves observed at $N=48$ equally spaced time points following the proposal of cho2013. To alleviate the impact of randomness from individual curves, we randomly split the data into $p$ groups of equal size, then take the sample average of curves within each group and finally smooth the averaged curves based on a 15-dimensional Fourier basis. We target to evaluate the $h$-day ahead predictive accuracy for the $p$-dimensional intraday energy consumption averaged curves in May 2013 based on the training data from December 2012 to the previous day. The eight comparison methods are built in the same manner as Section (ref) with $k_0= m = 5.$
Table (ref) presents the mean prediction errors for $h\in\{1,2,3\}$ and $p\in\{40, 80\}.$ A few trends are apparent. Firstly, the prediction errors for $p=80$ are higher than those for $p=40$ as higher dimensionality poses more challenges in prediction. Secondly, likewise in previous examples, SegV and FTSegV attain the lowest prediction errors in comparison to five competing methods under all scenarios. All segmentation-based methods consistently outperform UniV and GSY by a large margin. Thirdly, despite being developed for high-dimensional functional time series prediction, GSY provides the worst result in this example.