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.
64,300 characters · 12 sections · 71 citation commands
Covariance Function Estimation for High-Dimensional Functional Time Series with Dual Factor Structures
\centerline{\bf Abstract}
We propose a flexible dual functional factor model for modelling high-dimensional functional time series. In this model, a high-dimensional fully functional factor parametrisation is imposed on the observed functional processes, whereas a low-dimensional version (via series approximation) is assumed for the latent functional factors. We extend the classic principal component analysis technique for the estimation of a low-rank structure to the estimation of a large covariance matrix of random functions that satisfies a notion of (approximate) functional “low-rank plus sparse" structure; and generalise the matrix shrinkage method to functional shrinkage in order to estimate the sparse structure of functional idiosyncratic components. Under appropriate regularity conditions, we derive the large sample theory of the developed estimators, including the consistency of the estimated factors and functional factor loadings and the convergence rates of the estimated matrices of covariance functions measured by various (functional) matrix norms. Consistent selection of the number of factors and a data-driven rule to choose the shrinkage parameter are discussed. Simulation and empirical studies are provided to demonstrate the finite-sample performance of the developed model and estimation methodology.
{\em Keywords}: Covariance operator, functional factor model, functional time series, generalised shrinkage, high dimensionality, PCA, sparsity.
\setcounter{equation}{0}
A fundamental problem of increasing interest in modelling time series of random functions is to estimate their second-order characteristics, such as the covariance, auto-covariance and spectral density operators B00, HK10, HK12. Understanding these characteristics are not only important for understanding the randomness of the corresponding processes themselves, but also crucial to subsequent down-stream applications such as functional principal component analysis HK10, HK12, PT13, HKH15. While most of the existing literature focuses on a fixed number of functional time series, it is increasingly common to collect a large number of them in a diverse range of fields. For example, in climatology, temperature curves are routinely recorded in hundreds of weather stations, while in finance, stock return curves are typically available for thousands of stocks.
The main aim of this paper is to estimate the covariance structure of functional time series, when the dimensionality of the data (the number of functional time series) is comparable to or greater than the sample size (the length of these time series). To overcome the curse of dimensionality arising from this setup, a common practice is to make an approximate sparsity assumption on this structure. Consequently, regularisation methods such as thresholding or adaptive thresholding can be applied FGQ23. However, the sparsity restriction only works when these functional time series are at most weakly correlated, ruling out many interesting cases where they can be highly correlated. In practice, it is widely recognised that multivariate functional time series are often influenced by common functions over the temporal dimension, leading to strong cross-sectional dependence. For example, the rainfall curves collected in many stations may be affected by the common weather pattern in the region DHP19, and the intraday return curves are often driven by latent market and industry factors.
To accommodate cross-sectional dependence, an idea is to use factor models in the context of functional time series. Indeed, the approximate factor model has been extensively studied in the literature for panel time series CR83, BN02, SW02, FLM13. In recent years, there have been some attempts to extend this model to the functional time series setting, broadly classified into two directions. The first category can be loosely referred to as functional factor models for fixed-dimensional problems. Specifically, SH12 and KMZ15 consider a factor model with real-valued factor loadings and functional factors, whereas KMRT18 and MGG22 propose a factor model with real-valued factors and functional loadings. These functional factor models are essentially {\em low-dimensional}, because the developed methodology is not applicable to large-scale functional time series when the dimension is comparable to the sample size. The second direction is the {\em high-dimensional} functional factor model. GSY19 first apply the classic functional principal component analysis to each functional time series before modelling the component scores via a factor model. However, this method may lead to information loss in the dimension reduction stage, which subsequently may result in inaccurate factor number estimation. To address these, two types of functional factor models with different construction of common components have been proposed in the recent literature. TNH23a, TNH23b introduce a high-dimensional functional factor model with functional loadings and real-valued factors, whereas GQW21 propose a different version with real-valued loadings and functional factors.
In this paper, we propose a novel dual functional factor model with functional loadings and functional latent factors, in which the latent factors are modelled through series approximation. Towards this, we relax some restrictive assumptions imposed on existing functional factor models studied in the aforementioned literature. The proposed model includes those in TNH23a, TNH23b as special cases. We allow the covariance functions for the idiosyncratic components to be approximately sparse in a functional sense. Thus, the covariance function (or operator) of the high-dimensional functional process under our setup admits an approximate “low-rank plus sparse" functional covariance structure FLM13, extending the functional sparsity condition assumed in FGQ23 to allow latent factors.
The main estimation methodology can be seen as a functional analog of the POET method introduced by FLM13. For the factor loadings and common factors, we propose the use of a functional version of principal component analysis (PCA). For the covariance function of the idiosyncratic components, we apply a functional generalised shrinkage to a preliminary estimate of it. Under regularity conditions, we derive the mean squared convergence of the estimated factors, uniform convergence of the estimated functional factor loadings, and convergence of the estimated functional covariance matrices (measured by different matrix norms). These convergence rates depend on the dimension, the length of functional time series, the number of factors, and the error of sieve functional approximation. The Monte-Carlo simulation study shows that the developed methodology has reliable finite-sample performance. An empirical application to the cumulative intraday returns (CIDR) curves of the S&P 500 index data confirms the usefulness of the model and the accuracy of the suggested estimation approach.
After completing a preliminary draft of the paper, we found a recent working paper by LQW23 who consider a similar problem. It is worth comparing the two papers before concluding the introductory section. First, LQW23 estimate the large functional covariance matrix based on the high-dimensional functional factor models of TNH23a, TNH23b and GQW21, whereas our estimation is built on the dual functional factor model structure. More discussion and comparison between the models are provided in Section (ref) below. Second, we allow the number of factors (via the functional sieve approximation) to be divergent and the functional observations from different subjects to be defined on different domains, relaxing some restrictions implicitly imposed in LQW23. Third, we derive general convergence theory under the flexible model framework and the obtained convergence rates are comparable to those in LQW23 where the factor number is fixed and the sieve approximation error is zero. Finally, we use a different criterion to consistently estimate the factor number and choose the tuning parameter in functional shrinkage. In particular, we propose a modified cross-validation to select the shrinkage parameter, taking into account the temporal dependence of the functional data over a long time span.
The rest of the paper is organised as follows. Section (ref) introduces the dual functional factor model framework. Section (ref) describes the main estimation methodology and Section (ref) presents the large sample theory for the developed estimators. Section (ref) discusses practical issues for implementation and reports the numerical studies. Section (ref) concludes the paper. Proofs of the main theorems are given in Appendix A and proofs of some technical lemmas are available in Appendix B. Appendix C briefly reviews the $\tau$-mixing dependence and the concentration inequality. Throughout the paper, for a separable Hilbert space ${\mathscr H}$ defined as a set of real measurable functions $x(\cdot)$ on a compact set ${\mathbb C}$ such that $\int_{\mathbb C} x^2(u)du<\infty$, the inner product of $x_1,x_2\in{\mathscr H}$ is $\langle x_1,x_2\rangle=\int_{\mathbb C}x_1(u)x_2(u)du$, and the norm of $x\in{\mathscr H}$ is $\Vert x\Vert_2=\langle x,x\rangle^{1/2}$. Denote by ${\mathscr H}_1\times {\mathscr H}_2$ the Cartesian product of ${\mathscr H}_1$ and ${\mathscr H}_2$. Let $\Vert \cdot\Vert_O$ and $\Vert \cdot\Vert_{\rm S}$ respectively denote the operator and Hilbert-Schmidt norms for continuous linear operators. Denote the Euclidean norm of a vector by $\vert \cdot\vert_2$, and the operator, Frobenius and maximum norms of a matrix by $\Vert \cdot\Vert$, $\Vert\cdot\Vert_F$ and $\Vert\cdot\Vert_{\max}$, respectively. For $A=(a_1,\cdots, a_p)\in{\mathscr H}^p$ (the $p$-fold Cartesian product of ${\mathscr H}$), we define $|||A|||_2^2=\sum_{i=1}^p\int_{\mathbb C}a_i^2(u)du$. Let $a_{n}\sim b_{n}$, $a_{n}\propto b_{n}$ and $a_{n}\gg b_{n}$ denote that $a_{n}/b_{n}\rightarrow1$, $0<\underline{c}\leq a_{n}/b_{n}\leq\overline{c}<\infty$ and $b_{n}/a_{n}\rightarrow0$, respectively. Write “with probability approaching one" as “{\em w.p.a.1}" for brevity.
\setcounter{equation}{0}
Let $X_t=(X_{1t},\cdots,X_{Nt})^{^\intercal}$, $t=1,\cdots,T$, with $X_{it}=\left(X_{it}(u):\ u\in{\mathbb C}_i\right)\in{\mathscr H}_i$, where ${\mathscr H}_i$ is a Hilbert space defined as a set of measurable and square-integrable functions on a bounded set ${\mathbb C}_i$. Throughout the paper, we use $i$ to denote the subject index and $t$ to denote the time index. We propose the following functional factor model:
where $F_{tk}=\left(F_{tk}(u):\ u\in{\mathbb C}_k^\ast\right)\in{\mathscr H}_k^\ast$, ${\mathbb C}_k^\ast$ is a bounded set which may be different from ${\mathbb C}_i$, ${\mathscr H}_k^\ast$ is defined similarly to ${\mathscr H}_i$ but replacing ${\mathbb C}_i$ by ${\mathbb C}_k^\ast$, ${\mathbf B}_{ik}$ is a continuous linear operator from ${\mathscr H}_k^\ast$ to ${\mathscr H}_i$, and $\varepsilon_{it}=\left(\varepsilon_{it}(u):\ u\in{\mathbb C}_i\right)\in{\mathscr H}_i$. Neither the factor loading operator ${\mathbf B}_{ik}$ nor functional factor $F_{kt}$ is observable. The number $k_\ast$ is unknown but assumed to be finite. Model ((ref)) is a natural extension of the conventional factor model to a panel of functional time series. It follows from ((ref)) that each functional observation $X_{it}$ is decomposed into the functional common and idiosyncratic components: $\chi_{it}$ and $\varepsilon_{it}$. Write $F_t=(F_{t1},\cdots,F_{tk_\ast })^{^\intercal}$, $\varepsilon_t=(\varepsilon_{1t},\cdots,\varepsilon_{Nt})^{^\intercal}$, $\chi_t=(\chi_{1t},\cdots,\chi_{Nt})^{^\intercal}$ and \[ {\mathbf B}=\left(
\right), \] an $N\times k_\ast$ matrix of continuous linear operators from ${\mathscr H}_F:={\mathscr H}_1^\ast\times\cdots\times {\mathscr H}_{k_\ast}^\ast$ to ${\mathscr H}_X:={\mathscr H}_1\times\cdots\times {\mathscr H}_N$. Then model ((ref)) can be written in the following vector/matrix form:
To facilitate comparison with the functional factor models proposed in the recent literature, we next consider ${\mathbf B}_{ik}$ as a linear integral operator with kernel $B_{ik}=\left(B_{ik}(u,v):\ u\in{\mathbb C}_i, v\in{\mathbb C}_k^\ast\right)$, i,e., \[ {\mathbf B}_{ik}(x)(u)=\int_{{\mathbb C}_k^\ast} B_{ik}(u,v)x(v)dv,\ \ x\in{\mathscr H}_k^\ast. \] Model ((ref)) can be re-written as
TNH23a introduce the following model:
where $B_{ik}^\dag(\cdot)$ is the functional factor loading and $F_t^\dag=\left(F_{t1}^\dag,\cdots,F_{tk_\dag}^\dag\right)^{^\intercal}$ is a vector of real-valued factors. It is easy to see that model ((ref)) is a special case of ((ref)). Letting the matrix ${\mathbf B}$ of operators in ((ref)) be replaced by a matrix of real-valued factor loadings, we obtain the model in GQW21:
where $B_{ik}^\ddag$ is the real-valued factor loading and $F_{tk}^\ddag(\cdot)$ is the functional factor. The factor numbers $k_\ast$ in (ref) and $k_\ddag$ in (ref) are assumed to be fixed positive integers. As in HG18 and TNH23a, TNH23b, we allow the functional observations $X_{it}(\cdot)$, $i=1,\cdots,N$, and the latent factors $F_{tk}(\cdot)$, $k=1,\cdots,k_\ast$, to be defined on different domains, i.e., ${\mathbb C}_i$ and ${\mathbb C}_k^\dag$ may vary over $i$ and $k$. For instance, the collected high-dimensional functional time series may be a combination of time series with values in spaces with different dimensions (for different subjects). In contrast, GQW21 implicitly assume that $X_{it}(\cdot)$ and $F_{tk}^\ddag(\cdot)$ are defined on the same domain.
For the model in ((ref))--((ref)), we further impose a low-dimensional functional factor condition on the latent factor $F_t$ by assuming the following series approximation:
where $\Phi_k(\cdot)$ is a $q$-dimensional vector of basis functions, $G_t$ is a $q$-dimensional vector of stationary random variables, $\eta_{tk}\in{\mathscr H}_k^\ast$ is the sieve approximation error satisfying Assumption (ref)(iii), and $q$ is a positive integer which may slowly diverge to infinity. Note that the set of basis functions in ((ref)) is allowed to vary over $k$. Model ((ref)) extends the models studied in KMRT18 and MGG22 for univariate functional time series to multivariate ones.
Letting \[\Lambda_i(u)=\sum_{k=1}^{k_\ast} {\mathbf B}_{ik}\Phi_k(u)=\sum_{k=1}^{k_\ast}\int_{{\mathbb C}_k^\ast}B_{ik}(u,v)\Phi_k(v)dv\] and \[ \chi_{it}^\eta(u)=\sum_{k=1}^{k_\ast}{\mathbf B}_{ik}\eta_{tk}(u)=\sum_{k=1}^{k_\ast}\int_{{\mathbb C}_k^\ast}B_{ik}(u,v)\eta_{tk}(v)dv,\] the dual functional factor model structure combining ((ref)) and ((ref)) leads to
where $\chi_{it}^\ast(u)=\Lambda_{i}(u)^{^\intercal}G_{t}$ and $\varepsilon_{it}^\ast(u)=\chi_{it}^\eta(u)+\varepsilon_{it}(u)$. By Assumption (ref)(ii)--(iii) in Section (ref) below and the Cauchy-Schwarz inequality, we may show that $\chi_{it}^\eta$, the common component driven by the series approximation error, converges to zero {\em w.p.a.1}. If $\eta_{tk}\equiv0$, model ((ref)) reduces to ((ref)) as in TNH23a, TNH23b but with a possibly diverging number of factors. Our main interest is to estimate the large contemporaneous functional covariance structure under the general functional factor model framework ((ref)).
\setcounter{equation}{0}
Without loss of generality, assume that $X_{it}, \varepsilon_{it}, F_{kt}$ and $G_t$ all have zero mean. Define ${\mathbf C}_X={\sf E}\left(X_tX_t^{^\intercal}\right)=\left({\mathbf C}_{X,ij}\right)_{N\times N}$, where ${\mathbf C}_{X,ij}={\sf E}\left(X_{it}X_{jt}\right)$ is the covariance operator between $X_{it}$ and $X_{jt}$ with the corresponding kernel $C_{X,ij}(u,v)={\sf E}\left[X_{it}(u)X_{jt}(v)\right]$. Assuming $\chi_t$ and $\varepsilon_t$ are uncorrelated, we may decompose ${\mathbf C}_X$ as
where ${\mathbf C}_{\chi,ij}$ is the covariance operator between the functional common components $\chi_{it}$ and $\chi_{jt}$, and ${\mathbf C}_{\varepsilon,ij}$ is the covariance operator between the functional idiosyncratic components $\varepsilon_{it}$ and $\varepsilon_{jt}$\footnote{Throughout the paper, we use ${\mathbf C}_{X,ij}, {\mathbf C}_{\chi,ij}$ and ${\mathbf C}_{\varepsilon,ij}$ to denote the operators, and $C_{X,ij}(u,v), C_{\chi,ij}(u,v)$ and $C_{\varepsilon,ij}(u,v)$ to denote their respective kernels.}. Let $C_X=\left(C_{X,ij}\right)_{N\times N}, C_\chi=\left(C_{\chi,ij}\right)_{N\times N}$ and $C_\varepsilon=\left(C_{\varepsilon,ij}\right)_{N\times N}$ be matrices of covariance functions (or kernels) corresponding to ${\mathbf C}_X$, ${\mathbf C}_\chi$ and ${\mathbf C}_\varepsilon$, respectively.
It follows from ((ref)) that $C_{\chi,ij}(u,v)={\sf E}[\chi_{it}(u)\chi_{jt}(v)]$ can be approximated by
where $\Sigma_G$ is a $q\times q$ covariance matrix of latent factors $G_t$ and $\Lambda_i(u)$ is defined as in ((ref)). Hence $C_\chi$ can be approximated by $C_\chi^\ast=\left(C_{\chi,ij}^\ast\right)_{N\times N}$, a low-rank functional covariance matrix. On the other hand, since the functional idiosyncratic components are weakly cross-sectionally correlated, it is sensible to impose a functional version of the approximate sparsity restriction on $C_\varepsilon$, i.e.,
where $\varpi_N$ denotes a positive number which depends on $N$, and “$\Sigma\succeq0$" denotes that $\Sigma$ is a positive semi-definite matrix of covariance functions, i.e., $$\sum_{i=1}^N\sum_{j=1}^N \int_{u\in{\mathbb C}_i}\int_{v\in{\mathbb C}_j}x_i(u)\Sigma_{ij}(u,v)x_j(v)dudv\geq0,\ \ \forall\ x=(x_1,\cdots,x_N)^{^\intercal}\in{\mathscr H}_X.$$ Combining the above arguments, we obtain the functional low-rank plus sparse structure for
which may serve as a proxy of $C_X$. Under some mild conditions such as Assumption (ref)(iii), we may show $C_X^\ast$ converges to $C_X$ (under the matrix maximum norm), using Lemma (ref). With ((ref)), we may extend the POET method in FLM13 to estimate $C_X^\ast$ in the high-dimensional functional data setting.
Due to the general functional structure for the common components in ((ref))--((ref)), it is practically infeasible to directly estimate the factor loading operators and functional factors in the latent structure. Instead, motivated by ((ref))--((ref)), we may make use of the low-dimensional functional factor model ((ref)) and estimate $\Lambda_i(\cdot)$ and $G_t$ (subject to rotation) in the latent $\chi_{it}^\ast$ which approximates $\chi_{it}$. The latter can be done by extending the PCA technique BN02, SW02 to a large panel of functional observations. Assume the number of real-value factors in ((ref)) is known for the time being. We will discuss on how to determine this number in Section (ref). We estimate the factor loading functions and real-valued factors by minimising the following least squares objective function: \[ \sum_{t=1}^T \left\Vert X_{t}-\lambda g_t\right\Vert_{N,2}^2=\sum_{t=1}^T\sum_{i=1}^N \left\Vert X_{it}-\lambda_i g_t\right\Vert_2^2, \] where $\Vert\cdot\Vert_{N,2}$ denotes the norm of elements in ${\mathscr H}_X$, $\lambda_i=\left(\lambda_{i}(u):\ u\in{\mathbb C}_i\right)$ is a $q$-dimensional vector of functions and $g_t$ is a $q$-dimensional vector of numbers. Minimisation of the above least squares objective function can be achieved via the eigenanalysis of
Consider the following identification condition in the PCA algorithm:
which are similar to those in BN02 and FLM13. With ((ref)), we let $\widetilde G=\big(\widetilde G_1,\cdots,\widetilde G_T\big)^{^\intercal}$ as a matrix consisting of the eigenvectors (multiplied by root-$T$) corresponding to the $q$ largest eigenvalues of $\Delta$ defined in ((ref)). The factor loading functions are estimated as \[\widetilde\Lambda_i=\left(\widetilde\Lambda_i(u): u\in{\mathbb C}_i\right)=\frac{1}{T}\sum_{t=1}^TX_{it}\widetilde G_t,\ \ i=1,\cdots,N,\] via the least squares, using the normalisation restriction $\frac{1}{T}\sum_{t=1}^T\widetilde G_t\widetilde G_t^{^\intercal}=I_q$ by ((ref)). Consequently, the low-rank matrix of covariance functions $C_\chi^\ast$ is estimated by
We next turn to the estimation of $C_\varepsilon$. Letting $ \widetilde\varepsilon_{it}=X_{it}-\widetilde\Lambda_{i}^{^\intercal}\widetilde G_t$ be the approximation of $\varepsilon_{it}$, it is natural to estimate $C_\varepsilon$ by \[ \widehat{C}_\varepsilon=\left(\widehat{C}_{\varepsilon,ij}\right)_{N\times N}\ \ {\rm with}\ \ \widehat{C}_{\varepsilon,ij}=\left(\widehat C_{\varepsilon,ij}(u,v):\ u\in{\mathbb C}_i,\ v\in{\mathbb C}_j\right)=\frac{1}{T}\sum_{t=1}^T \widetilde\varepsilon_{it}\widetilde\varepsilon_{jt}. \] However, this matrix of conventional sample covariance functions often performs poorly when the number $N$ is comparable to or larger than $T$. To address this problem, we adopt the generalised shrinkage and estimate $C_\varepsilon$ by
where $s_\rho$ is a functional thresholding operator satisfying (i) $\Vert s_\rho(C)\Vert_{\rm S} \leq \Vert C\Vert_{\rm S}$ for any covariance function (or operator) $C$; (ii) $\Vert s_\rho(C)\Vert_{\rm S}=0$ if $\Vert C\Vert_{\rm S} \leq \rho$; and (iii) $\Vert s_\rho(C)-C\Vert_{\rm S}\leq \rho$, where $\rho$ is a user-specified tuning parameter controlling the level of shrinkage. A modified cross-validation method will be given in Section (ref) to determine $\rho$, taking into account the temporal dependence of high-dimensional functional data. We may further replace the universal thresholding by an adaptive functional thresholding as suggested by FGQ23 and LQW23, in which case the theory to be developed in Section (ref) remains valid.
Combining $\widetilde{C}_\chi$ and $\widetilde{C}_\varepsilon$, we finally obtain the estimate of $C_X^\ast$ or $C_X$:
\setcounter{equation}{0}
In this section, we first give some regularity conditions and then present the convergence properties for the estimates developed in Section (ref).
We derive the convergence properties in the so-called large panel setting, i.e., when $N$ and $T$ diverge to infinity jointly. We start with the convergence property for the PCA estimators of factors and functional factor loadings. Let $V$ be a $q\times q$ diagonal matrix with the diagonal elements being the first $q$ largest eigenvalues of $\frac{1}{T}\Delta$ (arranged in the decreasing order), and define the following $q\times q$ rotation matrix: \[R=V^{-1}\left(\frac{1}{T}\widetilde G^{^\intercal}G\right)\left[\frac{1}{N}\sum_{i=1}^N \int_{u\in{\mathbb C}_i}\Lambda_i(u)\Lambda_i(u)^{^\intercal}du\right],\] where $G=\left(G_1,\cdots,G_T\right)^{^\intercal}$ and $\widetilde G=\left(\widetilde G_1,\cdots,\widetilde G_T\right)^{^\intercal}$. We may show that this random rotation matrix $R$ is asymptotically invertible, see Lemma (ref).
\setcounter{prop}{0}
The following theorem gives the convergence rates for the estimated covariance functions for the functional idiosyncratic components.
\setcounter{theorem}{1}
The following theorem states the uniform convergence property for $\widetilde{C}_{X}$ defined in ((ref)).
Due to the low-rank plus sparse functional matrix structure for $C_X^\ast$ and $C_X$, we cannot derive the uniform convergence property in the functional version of matrix $\ell_1$ norm as in ((ref)). To address this problem, FLM13 recommends measuring the relative error of large low-rank plus sparse covariance matrix estimation for a high-dimensional random vector. However, it seems difficult to directly extend their relative error measurement to the setting of high-dimensional functional data due to the following reasons. First, the inverse of the large covariance matrix is often involved in defining the relative error measurement, but it seems difficult to compute the inverse of the large matrix of covariance operators ${\mathbf C}_X$. Second, although the uniform approximation properties ((ref)) and ((ref)) hold in the functional version of the matrix maximum norm, it is non-trivial to derive similar approximation properties under the relative error measurement. Although it is difficult to provide a theoretical justification via the relative error measurement, in the simulation study we will report some results on two types of its discrete approximation, illustrating the finite-sample performance of $\widetilde{C}_{X}$.
\setcounter{equation}{0}
In this section, we first introduce an easy-to-implement criterion to consistently estimate the factor number $q$ and discuss selection of the tuning parameter $\rho$ in the functional shrinkage. Then, we present Monte-Carlo simulation and empirical studies.
The functional PCA algorithm proposed in Section (ref) requires accurate estimation of $q$, the dimension of latent $G_t$. There have been extensive studies on the factor number selection in the context of high-dimensional real-valued time series. For example, BN02 introduces some information criteria to consistently estimate the factor number, and LY12 and AH13 propose a simple ratio criterion by comparing ratios of the estimated eigenvalues. The factor number is often assumed to be fixed and does not change with the size of data in the aforementioned literature. In this paper, we introduce a modified information criterion to consistently estimate the factor number $q$, which can diverge to infinity slowly\footnote{LLS17 consider a similar scenario in the factor number selection for real-valued panel time series.}. Let $\nu_k(\Delta/T)$ be the $k$\textsuperscript{th} largest eigenvalue of $\Delta/T$ with $\Delta$ defined in ((ref)), and define
where $\phi_{N,T}$ in the penalty parameter and $q_{\max}$ is a user-specified positive integer. A similar criterion is also adopted by AX17 to determine the number of high-frequency latent factors. The selection criterion in ((ref)) amends BN02's information criterion which replaces $\nu_k(\Delta/T)$ by summation of $\nu_j(\Delta/T)$ over $j>k$ and does not require “-1" adjustment. It is worth pointing out that the amended information criterion (ref) is easier to implement and its theoretical justification is straightforward. Proposition (ref) below derives its consistency property.
\setcounter{prop}{0}
The shrinkage estimation of $C_\varepsilon$ is often sensitive to the choice of the tuning parameter $\rho$. BL08 recommend the cross-validation method to select $\rho$ when the observations are independent and real-valued, see also FGQ23 for the extension to functional-valued observations. However, since the functional time series observations are serially correlated over time, satisfying the functional factor model structure, the cross-validation method may no longer work well in our setting. We modify the traditional cross-validation as follows.
The above selection criterion is introduced by CLL19 for large covariance matrix estimation of weakly dependent real-valued time series. The reason for leaving out $K$ observations between the two subsamples in each rolling window is to make these two subsamples have negligible correlation. In practice, for weakly dependent functional time series, we may set $K=10$.
We start with the description of the data generating process. For simplicity, let $k_\ast=1$ in the fully functional factor model ((ref)) or ((ref)). The functional factor process $\{F_t: t=1,2,\cdots\}$ is defined by ((ref)), i.e., \[ F_{t}(u)=\Phi(u)^{^\intercal}G_t+\frac{1}{q}\eta_t^\ast(u),\ \ u\in{\mathbb C}=[0,1], \] where $\eta_{t}^\ast(\cdot)$ is generated from a Brownian bridge, $\Phi(\cdot)$ is a $q$-dimensional vector of Fourier basis functions, $G_t$ follows a VAR$(1)$ model: \[ G_t = {\mathbf A}G_{t-1}+\zeta_t \] with ${\mathbf A}$ being a $q\times q$ coefficient matrix with its element $A_{jk} = 0.25^{|j-k|+1}$ for $j,k = 1,\dots,q$ and $\zeta_t$ being independently generated by a $q$-dimensional standard normal distribution. Each function is generated on a common set of 21 equally-spaced grid points on $[0,1]$. The number of real-valued factors, $q$, is set as $5,10$ or $15$. We simulate the factor loadings $B_{i}(u,v)$ via \[ B_{i}(u,v) = \sum^{50}_{j=1}c_{i,j}\phi_j(u)\phi_j(v), \] where $c_{i,j}$ is independently generated from a standard normal distribution, and $\{\phi_1,\cdots,\phi_{50}\}$ is a set of Fourier basis functions.
As in FGQ23, we simulate the functional idiosyncratic term as \[ \varepsilon_{it}(u)=\sum_{j=1}^{50}\frac{1}{j}\theta_{it,j}\phi_j(u),\quad i=1,\cdots,N,\quad t=1,\cdots,T, \] where $\phi_j$ is the Fourier basis function and $\theta_t=(\theta_{1t}^{^\intercal},\cdots,\theta_{Nt}^{^\intercal})^{^\intercal}$ with $\theta_{it}=(\theta_{it,1},\cdots,\theta_{it,50})^{^\intercal}$ are independently generated from a multivariate Gaussian distribution with mean zero and covariance matrix $\Sigma\in {\mathscr R}^{50N\times 50N}$. Write $\Sigma=(\Sigma_{jk})_{N\times N}$ with $\Sigma_{jk}\in {\mathscr R}^{50\times 50}$ being the $(j,k)$ block, $1\leq j,k\leq N$. The functional sparsity pattern on $C_\varepsilon$ may be characterised by a sparsity structure in $\Sigma$. In the simulation study, we define $\Sigma_{jk} = \omega_{jk}\Omega$ with $\Omega=\text{diag}(1^{-2},\cdots,50^{-2})$ and \[ \omega_{jk} = \left\{
\right. \] where $I(\cdot)$ denotes the binary indicator function.
We apply the developed estimation method to the simulated data and compute the difference between the true and estimated functional covariance matrices over 200 replications. As in FGQ23, we compute the functional version of $\ell_1$ and $\ell_2$ matrix norms for the functional idiosyncratic covariance matrix estimation (see Table (ref)), and the matrix maximum norms (see Table (ref)), providing the finite-sample justification for Theorem (ref). These functional $\ell_1$, $\ell_2$ and maximum matrix norms are defined as
We consider four shrinkage functions in the functional shrinkage technique: hard thresholding (Hard), soft thresholding (Soft), SCAD and adaptive lasso (Alasso), and include the sample covariance function estimation (sample) as a benchmark. The tuning parameter involved in the functional shrinkage is selected via the modified cross-validation introduced in Section (ref). Table (ref) reports the functional $\ell_1$- and $\ell_2$-norm estimation errors for functional idiosyncratic covariance matrices. Both the $\ell_1$- and $\ell_2$-norm estimation errors decrease significantly when $T$ increases from 100 to 200; the $\ell_1$-norm shrinkage estimation errors slightly increase when $N$ increases from $50$ to 200 but $T$ is fixed; whereas the increase of $\ell_2$-norm estimation errors is more substantial. The estimation performance is generally stable as $q$ increases from 5 to 15. The use of functional shrinkage significantly outperforms the naive sample covariance function without shrinkage. In particular, the adaptive lasso performs best among the four shrinkage functions. Table (ref) further compares the functional max-norm estimation errors between the adaptive lasso\footnote{The results are almost the same for the other three shrinkage methods.} and sample covariance function. Unlike Table (ref), their estimation performance is very close.
To measure the covariance matrix estimation accuracy of the observed functional observations simulated by the dual functional factor model, due to the existence of spiked eigenvalues, we cannot adopt the functional $\ell_1$- and $\ell_2$-norm estimation errors both of which are divergent as $N$ increases. As discussed in Remark (ref), we define the following two relative-error measurements:
whose results are reported in Table (ref). The relative-error measurements increase when $N$ increases, but decrease when $T$ increases from $100$ to $200$. The adoptive lasso again outperforms the other three shrinkage methods, which is consistent with the finding in Table (ref).
We assess the performance of the amended information criterion in Table (ref) by reporting percentages of accurately estimating the true number of factors. When $N=100$ and $200$, we achieve perfect accuracy in selecting the factor number; when $N=50$, the proposed information criterion may occasionally over-estimate the factor number, which has negligible impact on the subsequent covariance function estimation.
We apply the developed method to analyse the functional covariance structure of the S&P 500 data containing $473$ common stocks traded in the years 2019 and 2020, which are available at the Refinitiv Datascope\footnote{\url{https://select.datascope.refinitiv.com/DataScope/}}. We consider their CIDR curves from 2 January 2019 to 31 December 2020. The S&P 500 index comprises some of the largest companies traded in New York Stock Exchange. After removing the public holidays and half-day trading days (like Christmas Eve), there are $T_1=242$ trading days in 2019 and $T_2=228$ in 2020. For each trading day, we consider $5$-minute resolution data between 9:30 and 16:00 Eastern Standard Time, and obtain $78$ data points. For asset $i$, let $P_{it}(u_j)$ be the intraday 5-minute close price at time $u_j$ on trading day $t$, and construct a sequence of CIDRs RWZ20: \[ X_{it}(u_j) = 100\times [\ln P_{it}(u_j) - \ln P_{it}(u_1)], \] where $\ln(\cdot)$ denotes the natural logarithm, $j=2,3,\cdots,78$ and $i=1,2,\cdots,473$. The way CIDR is constructed removes the effect of the starting price. The linear interpolation algorithm in HAB+23 is adopted to convert discrete data points into a continuous function. In Figure (ref), we display the CIDRs for Apple Inc., one of the most liquid stocks, in 2019 and 2020.
We fit the dual functional factor model to the S&P 500 time series and apply the estimation methodology developed in Section (ref). We consider not only the entire time period 2019--2020, but also the two calendar years 2019 and 2020 separately. With the modified information criterion, we select the number of real-valued factors (i.e., $G_t$) to be one for 2019, two for 2020, and two for 2019--2020. As recommended in the simulation, we use the functional shrinkage with adaptive lasso to estimate the functional idiosyncratic covariance matrix. With the amended cross-validation in Section (ref), we select the tuning parameter as 0.160 for 2019, 0.230 for 2020, and 0.114 for 2019--2020.
With the estimated covariance functions $\widetilde{C}_{X,ij}$ between stocks, we may further compute the correlation functions whose Hilbert-Schmidt norms are comparable over stocks. Define \[ \widetilde{\sf cor}(X_i, X_j) = \frac{\widetilde{C}_{X,ij}}{\Vert \widetilde{C}_{X,ii}\Vert_S \Vert\widetilde{C}_{X,jj}\Vert_S}, \] where $\Vert\cdot\Vert_S$ denotes the Hilbert-Schmidt norm of a covariance function (or operator). The correlation functions between the functional idiosyncratic components are computed in a similar way. The Hilbert-Schmidt norms of the estimated correlation functions are used to plot the heat maps displayed in Figure (ref). The heat maps show that the functional correlation (or covariance) structure among the stocks is dense whereas that among the functional idiosyncratic components (after removing the common factors) is very sparse, justifying the proposed (approximate) functional low-rank plus sparse structure. Furthermore, we also note that the correlations among the stocks are stronger in 2020 than in 2019, which may be due to a stronger co-movement after the declaration of the COVID-19 pandemic in early 2020.
\setcounter{equation}{0}
In this paper, we introduce a general dual functional factor model framework to tackle high-dimensional functional time series, extending recent proposals in the literature GQW21, TNH23a, TNH23b. The main model combines a high-dimensional fully functional factor model for the functional observations and a low-dimensional one for the latent functional factors. Using a sieve approximation for the functional factors, we approximate the large matrix of covariance functions (for the observed functional time series) by the functional low-rank plus sparse structure, making it feasible to extend FLM13's POET method to the large-scale functional data setting. A functional version of principal component analysis is proposed to estimate the functional factor loadings and common factors, and a functional shrinkage technique is adopted to estimate the covariance structure of the functional idiosyncratic components. Under some mild assumptions, we establish the convergence properties of the developed estimators. An amended information criterion is proposed to consistently estimate the factor number whereas a modified cross-validation is used to select the shrinkage parameter. Simulation study is provided to demonstrate reliable finite-sample performance of the developed methodology and the empirical application confirms the rationality of adopting the functional low-rank plus sparse covariance (or correlation) structure for CIDR curves of S&P 500 index.
Leng’s research is partly supported by EPSRC (EP/X009505/1). Li's research is partly supported by the Leverhulme Research Fellowship (RF-2023-396), Australian Research Council Discovery Project (DP230102250) and Heilbronn Institute for Mathematical Research. Shang's research is partly supported by the Australian Research Council Discovery Project (DP230102250). Xia's research is partly supported by the National Natural Science Foundation of China (72033002).