Large-Scale Curve Time Series with Common Stochastic Trends

Large-Scale Curve Time Series with Common Stochastic Trends

\centerline{\bf Abstract}

This paper studies high-dimensional curve time series with common stochastic trends. A dual functional factor model structure is adopted with a high-dimensional factor model for the observed curve time series and a low-dimensional factor model for the latent curves with common trends. A functional PCA technique is applied to estimate the common stochastic trends and functional factor loadings. Under some regularity conditions we derive the mean square convergence and limit distribution theory for the developed estimates, allowing the dimension and sample size to jointly diverge to infinity. We propose an easy-to-implement criterion to consistently select the number of common stochastic trends and further discuss model estimation when the nonstationary factors are cointegrated. Extensive Monte-Carlo simulations and two empirical applications to large-scale temperature curves in Australia and log-price curves of S&P $500$ stocks are conducted, showing finite-sample performance and providing practical implementations of the new methodology.

{\em Keywords}: Common trends, Curve time series, Factor models, Functional PCA, High dimensionality.

Introduction

\setcounter{equation}{0}

The past few decades have seen notable developments in modeling curve time series, a sequence of random curves or functions often defined within a bounded set. Many estimation, inference and forecasting techniques have been proposed to tackle curve time series B00,RS05,HK12,PJ25a which arise in a variety of areas such as climatology, transportation, finance, demography and health sciences. One may reduce the infinite dimension of curve time series to a finite dimension through a functional version of principal component analysis (PCA), and subsequently apply classic time series models such as VAR Lu06 to the finite-dimensional time series which retain much of the dynamic sample information. The existing literature often assumes the curve time series to be stationary, thereby facilitating theory development using standard asymptotics. Stationarity may be too restrictive in many applications and is often rejected when we test practical curve time series data. For example, CKP16 find evidence of a unit root structure for intra-month distribution curves of S&P 500 index returns; ARS18 reject the null hypothesis of stationarity for Australian temperature curves; LRS23 detect a nonstationary feature in US treasury yield curves; and PJ25b find unit root behavior in Engel curve data for leisure, health, and food expenditure among ageing seniors in Singapore.

There have been some attempts in recent years to relax the stationarity restriction for curve time series. HKR14 introduce a functional KPSS test KPSS92 for stationarity of curve time series; CKP16 study nonstationary time series of state density curves by decomposing an infinite-dimensional Hilbert space into the nonstationary I(1) and stationary subspaces; BSS17 establish the Granger-Johansen representation theorem for I(1) autoregressive curve processes, which has been further extended by BS20 and FP20 to I(2) and more general I($d$) autoregressive curve processes; LRS23 introduce a nonstationary fractionally integrated curve time series framework, covering the nonstationary I(1) curve as a special case; NSS23 propose a variance ratio-type test to determine the dimension of the nonstationary subspace of the cointegrated curve time series; and PJ25b develop ADF and semiparametric unit root tests for curve time series autoregression. The aforementioned literature limits attention to a single curve time series with nonstationarity. P25 considers rank selection in vector autoregression with multiple curve time series but in a parametric setting. In practice, we often have to jointly model a large number of curve time series driven by some common stochastic trends which are usually latent. For example, thousands of stock return curve time series in financial markets may be driven by latent market and industry factors; curve time series of temperature and rainfall recorded in hundreds of weather stations may be affected by common weather patterns in the region. Hence, it is imperative for adequate empirical modeling to develop a flexible curve time series framework that accommodates large dimensionality and unobserved factors as well as nonstationarity.

The approximate factor model has proven to be an effective tool for analyzing large-scale real-valued panel data CR83, BN02. BN04 proposed a so-called PANIC method under the factor model framework to test for unit roots in the idiosyncratic components and determine the number of common stochastic trends that are present among the cointegrated nonstationary factors. That work was extended in BC09 to accommodate structural breaks. B04 used classic PCA to estimate common stochastic trends and factor loadings, assuming the idiosyncratic components to be stationary over time; and BLL21 examined an approximate factor model for nonstationary panels with the primary concern of impulse-response function estimation with cointegrated factors within a vector error-correction model (VECM) specification.

The main focus of the present paper centers on the interaction of recent advances in nonstationary curve time series and large-dimensional approximate factor models. The goal is to build a fully functional factor model approach designed for curve time series with common stochastic trends. There has been increasing interest in extending the approximate factor model to curve time series under stationarity conditions. For a single or a small number of curve time series, HSH12, KMZ15 and KMRT18 consider low-dimensional functional factor models, where either factors or factor loadings take functional values. For large-scale curve time series with the cross-sectional size increasing with the temporal dimension, GQW21 consider a high-dimensional functional factor model with functional factors and real-valued loadings, whereas TNH23a,TNH23b introduce a different functional factor model with functional loadings and real-valued factors, proposing functional PCA to estimate the functional common and idiosyncratic components. In addition, LLSX24 recently introduced a dual functional factor model for high-dimensional stationary curve time series, providing estimates of the functional covariance structure. As far as we know, there is no literature on high-dimensional factor models for nonstationary curve time series. The present paper employs such a framework, adopting a dual functional factor model structure that admits common stochastic trends in high-dimensional curve time series, thereby allowing practical implementation with many financial market and climatic curve time series that manifest nonstationary behavior.

With a high-dimensional functional factor structure for large-scale curve time series, we decompose each functional observation into common and idiosyncratic components. In particular, we define the common component via an integral operator and allow both the factors and factor loadings to be functional, giving a more flexible structure than those in GQW21 and TNH23a, TNH23b. For the latent factor curves with common stochastic trends, we impose another functional factor model structure via multivariate series approximation, where the number of real-valued stochastic trends is allowed to diverge slowly to infinity. As in TNH23b, functional PCA methods are used to estimate the common stochastic trends and functional factor loadings. Under some technical but justifiable assumptions, we derive the mean square convergence of the estimated common trends, where the convergence rate relies on the dimension, time series length and number of common trends. To facilitate inference we establish limit theory for the estimated common trends and functional factor loadings, extending Theorems 2 and 3 in B04 to large-scale curve time series with a diverging number of common trends. In particular, super-fast convergence is established for the functional factor loading estimate and its limit distribution is derived as if the common stochastic trends were known.

Practical implementation of functional PCA, like other PCA applications, requires consistent estimation of the number of common stochastic trends. A suitable information-based selection criterion is employed for this purpose, modifying existing criteria that have been extensively studied in the literature BN02. The proposed criterion is easily implemented and consistency follows straightforwardly. A more general model setting is also considered in which the integrated factors are themselves cointegrated and the idiosyncratic components may be nonstationary. In this setting a functional version of PANIC is proposed for estimating factors (via first-order differences) and functional factor loadings.

Extensive simulations are conducted to assess numerical performance of the methods in finite samples. The findings reveal that when the factors are full-rank integrated and functional idiosyncratic components are stationary, functional PCA estimates of common stochastic trends and functional factor loadings are more efficient than those obtained via functional PANIC; but when the integrated factors are rank-reduced and functional idiosyncratic components are nonstationary, functional PANIC estimation continues to work well but functional PCA estimation is inconsistent. In the empirical applications, functional PCA is used to analyze temperature curve time series for Australia over the period 1943-2022, and functional PANIC is used to analyze log-price curves of S&P $500$ stocks from January 2023 to November 2023. The empirical results confirm the existence of common stochastic trends for both these datasets of large-scale curve time series.

The rest of the paper is organized as follows. Section (ref) introduces the dual functional factor model framework. Section (ref) describes functional PCA estimation and provides the relevant theory. Section (ref) proposes a modified information criterion to estimate the number of common stochastic trends. Section (ref) discusses model estimation with cointegrated factors. Sections (ref) and (ref) present the simulation study and the empirical applications. Section (ref) concludes. Proofs of the main asymptotic theorems are given in Appendix A. Some useful technical lemmas with proofs are in Appendix B. Throughout the paper, we define a separable Hilbert space ${\mathscr H}$ as a set of real measurable functions $f(\cdot)$ on a compact set ${\mathbb C}$ such that $\int_{\mathbb C} f^2(u)du<\infty$, with the inner product of $f_1$ and $f_2$ as $\langle f_1,f_2\rangle=\int_{\mathbb C}f_1(u)f_2(u)du$, and the norm as $\Vert f\Vert=\langle f,f\rangle^{1/2}$. We further generalize $\langle \cdot,\cdot\rangle$ to handle vectors or matrices of functions: for ${\boldsymbol F}_1=(f_{1,ki})$ and ${\boldsymbol F}_2=(f_{2,kj})$ which are $k_1\times k_2$ and $k_1\times k_3$ matrices of functions, define $\langle \boldsymbol F_1^{^\intercal},\boldsymbol F_2\rangle$ as a $k_2\times k_3$ matrix whose $(i,j)$-th entry is $\sum_{k=1}^{k_1}\int f_{1,ki}(u)f_{2,kj}(u)du$, and for a vector of functions $\boldsymbol F$, write $\Vert \boldsymbol F\Vert=\langle \boldsymbol F^{^\intercal},\boldsymbol F\rangle^{1/2}$. We also use the notation $\Vert\cdot\Vert$ as the Euclidean norm of a vector and the operator norm of a matrix or a continuous linear operator whenever no ambiguity arises and let $\Vert\cdot\Vert_F$ be the Frobenius norm of a matrix. 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. For brevity “with probability approaching one' is written “{\em w.p.a.1}".

Dual functional factor model structure

\setcounter{equation}{0}

$N$-vectors of curve time series $Z_t=(Z_{1t},\cdots,Z_{Nt})^{^\intercal}$, where $Z_{it}=\left(Z_{it}(u):\ u\in{\mathbb C}_{i}\right)\in{\mathscr H}_{i}$ for $t=1,\cdots,T$ are observed, with ${\mathscr H}_{i}$ being a separable Hilbert space defined as a set of measurable and square-integrable functions on a bounded set ${\mathbb C}_{i}$. We may decompose $Z_{it}$ into a common component and an idiosyncratic component as follows

where $\chi_{it}=\left(\chi_{it}(u):\ u\in{\mathbb C}_{i}\right)$ whose dynamic patterns are driven by some latent factors, and $\varepsilon_{it}=\left(\varepsilon_{it}(u):\ u\in{\mathbb C}_{i}\right)$ is allowed to be correlated over $i$ and $t$. The functional factor model ((ref)) is similar to the classic approximate factor model studied in CR83 and BN02 with the exception that all the components in ((ref)) take functional values. But the formulation of the common component $\chi_{it}$ is non-trivial. Broadly speaking, two different ways have been recommended in the recent literature to define $\chi_{it}$: GQW21 construct $\chi_{it}$ as a product of real-valued factor loadings and functional factors, whereas TNH23a define $\chi_{it}$ as a product of functional factor loadings and real-valued factors. The factor number is assumed fixed in GQW21 and TNH23a,TNH23b to achieve dimension reduction in large-dimensional curve time series modeling. Both approaches involve real-valued components, either as parametric factor loadings or as real-valued factors. As in LLSX24, we introduce a more flexible functional factor model, constructing $\chi_{it}$ via an integral operator and allowing both the factors and factor loadings to be functional.

Let $F_t=(F_{1t},\cdots,F_{kt})^{^\intercal}$, where $F_{jt}=\left(F_{jt}(u):\ u\in{\mathbb C}_{j}^\ast\right)\in{\mathscr H}_{j}^\ast$ and ${\mathscr H}_{j}^\ast$ is defined similarly to ${\mathscr H}_{i}$ but with ${\mathbb C}_{i}$ replaced by a possibly different bounded set ${\mathbb C}_{j}^\ast$. For $i=1,\cdots,N$ and $j=1,\cdots,k$, we let ${\cal B}_{ij}$ be a linear (kernel) integral operator defined by \[ {\cal B}_{ij}f(u)=\int_{{\mathbb C}_{j}^\ast} B_{ij}(u,v)f(v)dv,\ \ f\in{\mathscr H}_{j}^\ast,\ \ u\in{\mathbb C}_{i}, \] where $B_{ij}=\left(B_{ij}(u,v):\ u\in{\mathbb C}_{i}, v\in{\mathbb C}_{j}^\ast\right)$ denotes the kernel of the linear operator ${\cal B}_{ij}$, and write $\chi_{it}$ as

Combining ((ref)) and ((ref)), we obtain a fully functional factor model structure which is more general than those in GQW21 and TNH23a,TNH23b. Neither the factor loading operator ${\cal B}_{ij}$ nor functional factor $F_t$ is known a priori. As in HG18, we allow the curve time series observations $Z_{it}$, $i=1,\cdots,N$, and the latent functional factors $F_{jt}$, $j=1,\cdots,k$, to be defined on different domains, i.e., ${\mathbb C}_{i}$ and ${\mathbb C}_{j}^\ast$ may vary over $i$ and $j$. The number of functional factors is unknown but assumed to be a finite positive integer.

Although the linear integral operator provides a flexible structure for functional common components, it makes the estimation of functional factors and factor loadings challenging. To address the difficulty we impose a low-dimensional functional factor model representation for $F_t$ via the following series expansion

where $\Phi_j=\left(\phi_{j1},\cdots,\phi_{jq}\right)^{^\intercal}$ is a $q$-dimensional vector of deterministic basis functions, $G_t$ is a $q$-dimensional vector of nonstationary real-valued factors, $\eta_{jt}$ denotes the series approximation error which can be either stationary or nonstationary, and $q$ is a positive integer which may slowly diverge to infinity. Model ((ref)) extends the low-dimensional functional factor model studied in KMRT18 and MGG22 to multivariate curve time series. It is also similar to the multivariate Karhunen-Lo\`eve representation in HG18 if $\left(\phi_{1l},\cdots,\phi_{kl}\right)^{^\intercal}$ is an orthonormal basis vector of eigenfunctions and $q$ is set as the truncation parameter. In the present paper, the integrated factor $G_t$ is generated by

where $L$ is the lag operator and $\{\xi_t\}$ is a sequence of stationary $I(0)$ random vectors. Without loss of generality, we assume that the initial value $G_0=(G_{10},\cdots,G_{q0})^{^\intercal}$ satisfies $\max_{1\leq j\leq q}|G_{j0}|=O_P(1)$.

Writing

with the high-dimensional factor structure ((ref)) and ((ref)) for the observed curve time series and the low-dimensional factor structure ((ref)) for the latent factor curves, we obtain

where $\Lambda_i=\left(\Lambda_{i}(u):\ u\in{\mathbb C}_{i}\right)$ and $\chi_{it}^{\eta}=\left(\chi_{it}^{\eta}(u):\ u\in{\mathbb C}_{i}\right)$. For practical purposes our main interest lies in estimating $\Lambda_i$ and $G_t$, and determining $q$, the number of common stochastic trends. In the context of stationary curve time series, TNH23a, TNH23b's high-dimensional functional factor model can be seen as a special case of ((ref)) with $\chi_{it}^{\eta}\equiv0$ and $q$ being a finite positive integer. Model ((ref)) also extends the nonstationary factor model in B04, BN04 and BLL21 from real-valued time series to more general curve time series.

Functional PCA estimation methodology and theory

\setcounter{equation}{0}

This section introduces functional PCA methodology to estimate $\Lambda_i$ and $G_t$. PCA has been commonly used to estimate factors and factor loadings (subject to appropriate rotation) in the standard factor model for a large panel of real-valued time series BN02, SW02, B04, BN04, BLL21. We extend the technique to high-dimensional nonstationary curve time series. Here we assume the number of common stochastic trends is known and $G_t$ is a vector of full-rank integrated variables. Section (ref) introduces an easy-to-implement criterion to estimate $q$ and Section (ref) considers the more general setting of cointegrated factors.

Functional PCA

Since $\Lambda_i$ and $G_t$ are not identifiable in the functional factor model ((ref)), identification restrictions are imposed in the functional PCA algorithm using

where ${\boldsymbol I}_q$ is a $q\times q$ identity matrix. These identification conditions are comparable to those used by B04 for the traditional factor model. Since $G_t$ is integrated, the normalization rate in ((ref)) is $T^2$ instead of $T$ (for stationary time series). Eigenanalysis is conducted on the matrix

with $\widetilde {\boldsymbol G}=(\widetilde G_1,\cdots,\widetilde G_T)^{^\intercal}$ a $T\times q$ matrix consisting of the eigenvectors scaled by $T$, corresponding to the $q$ largest eigenvalues of $\widetilde{\boldsymbol\Omega}$. The functional factor loadings are subsequently estimated as

using least squares and the first restriction in ((ref)).

Mean square convergence of $\widetilde{G}_t$

The following assumptions are needed to develop the convergence theory of $\widetilde G_t$ and $\widetilde{\Lambda}_i$.

Let ${\boldsymbol V}_{NT}$ be a $q\times q$ diagonal matrix with its diagonal elements being the $q$ largest eigenvalues of $\frac{1}{T^2}\widetilde{\boldsymbol\Omega}$ (arranged in the decreasing order), and define the following $q\times q$ rotation matrix

where ${\boldsymbol G}=\left(G_1,\cdots,G_T\right)^{^\intercal}$. By Proposition (ref) and Lemma (ref), the limits as $\{N,T\} \to \infty$ of ${\boldsymbol V}_{NT}$ and ${\boldsymbol H}_{NT}$ are ${\boldsymbol V}_{0}$ and ${\boldsymbol H}_{0}$, respectively, where ${\boldsymbol V}_0$ is a $q\times q$ diagonal matrix with the diagonal elements being the eigenvalues of ${\boldsymbol\Sigma}_\Lambda^{1/2}(\int_0^1 B_\xi(u) B_\xi(u)^{^\intercal}du){\boldsymbol\Sigma}_\Lambda^{1/2}$ (arranged in decreasing order) and ${\boldsymbol H}_0={\boldsymbol V}_0^{-1/2}{\boldsymbol W}_0^{^\intercal}{\boldsymbol \Sigma}_\Lambda^{1/2}$ where ${\boldsymbol W}_0$ is a matrix consisting of the eigenvectors of ${\boldsymbol\Sigma}_\Lambda^{1/2}(\int_0^1 B_\xi(u) B_\xi(u)^{^\intercal}du){\boldsymbol\Sigma}_\Lambda^{1/2}$. The following proposition derives the mean square convergence property for $\widetilde{G}_t$.

\setcounter{prop}{0}

Limit distribution of $\widetilde G_t$ and $\widetilde\Lambda_i$

To conduct inference on the estimated common trends and functional loadings, limit theory is needed, for which the following additional conditions are employed.

\setcounter{theorem}{1}

We next turn to the limit distribution theory of the functional factor loading estimate $\widetilde\Lambda_i$, which requires the following conditions.

Estimation of the factor number

\setcounter{equation}{0}

In practice, the number of latent nonstationary factors is unknown. Implementation of the functional PCA proposed in Section (ref) requires a consistent estimation of $q$. There have been extensive studies on determining the number of factors in the conventional factor model for real-valued time series. BN02 and B04 propose some information criteria to consistently estimate the factor number for a large panel of stationary and nonstationary time series; LY12 and AH13 recommend an easy-to-implement ratio criterion where ratios of consecutive estimated eigenvalues are compared; T18 and BT22 estimate factor numbers by randomised sequential testing using estimated eigenvalues. For the present setting, we here employ an easy-to-implement information criterion to consistently estimate the number of common stochastic trends.

Let $\widetilde\nu_i=\nu_i(\widetilde{\boldsymbol\Omega}/T^2)$ denote the $i$-th eigenvalue of $\widetilde{\boldsymbol\Omega}/T^2$. We start with the following proposition on the asymptotic orders of $\widetilde\nu_i$, which motivate the selection criterion.

\setcounter{prop}{0}

In particular, define the criterion

where the penalty parameter $\rho_{NT}$ satisfies some mild restrictions (see Theorem (ref)) and $q_{\max}$ is a user-specified upper bound of the factor number. This criterion was used in AX17 in the context of factor models for high-dimensional and high-frequency financial data. It can be viewed as a modification of the information criterion proposed in BN02 and B04, which replaces $\widetilde\nu_j$ by summation of $\widetilde\nu_i$ over $i>j$ and which does not require a “-1" adjustment. The modified information criterion ((ref)) is easier to implement and proof of consistency is simpler.

\setcounter{theorem}{1}

Estimation with cointegrated factors

\setcounter{equation}{0}

This section considers the case where $G_t$ is cointegrated, i.e., ${\boldsymbol\Sigma}_\xi$ has reduced rank $q^\dagger$ with $1\leq q^\dagger\leq q-1$. Following PP88, PP89, there exists a $q\times q$ orthogonal matrix ${\boldsymbol P}=({\boldsymbol P}_1, {\boldsymbol P}_2)$ with ${\boldsymbol P}_1$ and ${\boldsymbol P}_2$ dimensioned $q\times q^\dagger$ and $q\times q^{\ddag}$, such that

where $q^{\ddag}=q-q^\dagger$ is called the cointegrating rank, $G_{t1}$ is a $q^\dagger$-dimensional vector of full-rank integrated variables and $G_{t2}$ is a $q^{\ddag}$-dimensional vector of stationary time series. The rotation ((ref)) successfully separates out stationary and nonstationary components, the latter of which drive the common stochastic trends for large-scale curve time series. We may use the functional PCA method as in Section (ref) to estimate both the integrated factors $G_{t1}$ and stationary factors $G_{t2}$ B04. However, it would require consistent estimation of $q$ and $q^\dag$ (or $q^{\ddag}$), and different normalization rates for $G_{t1}$ and $G_{t2}$. In this section, we propose a different approach which is a functional version of the PANIC method. PANIC was introduced by BN04 for high-dimensional real-valued time series, allowing some idiosyncratic components to be nonstationary BLL21.

As in BN04, taking differences on both sides of ((ref)) gives

where $z_{it}=\Delta Z_{it}$, $\xi_t=\Delta G_t$ and $\varepsilon_{it}^\dagger=\Delta(\chi_{it}^{\eta}+\varepsilon_{it})$. Model ((ref)) can be seen as a functional factor model for high-dimensional stationary curve time series LLSX24. However, weak cross-section dependence of $\varepsilon_{it}^\dagger$ over $i$ may not be satisfied due to the presence of series approximation errors in the functional common components. Throughout this section, we only require $\Delta\varepsilon_{it}$ to be stationary over $t$, which implies that $\varepsilon_{it}$ may be integrated.

We next estimate the functional factor loadings $\Lambda_i$ and stationary factor vector $\xi_t$. Since $\Lambda_i$ and $\xi_t$ are not identifiable in the functional factor model ((ref)), we employ the following identification restrictions in the functional PCA algorithm

Unlike ((ref)), the adjusted normalization rate $T-1$ is used for the stationary factors $\xi_t$. Define

Conducting the eigenanalysis of $\widehat{\boldsymbol\Omega}$, we obtain $\widehat {\boldsymbol \xi}=(\widehat \xi_2,\cdots,\widehat \xi_T)^{^\intercal}$ as a matrix consisting of the eigenvectors scaled by $\sqrt{T-1}$, corresponding to the $q$ largest eigenvalues of $\widehat{\boldsymbol\Omega}$. It follows from Proposition 4.1 in LLSX24 that, under some mild conditions,

where \[ {\boldsymbol H}_{NT}^\dagger=\left({\boldsymbol V}_{NT}^{\dagger}\right)^{-1}\left(\frac{1}{T-1}\widehat{\boldsymbol \xi}^{^\intercal}{\boldsymbol \xi}\right)\left[\frac{1}{N}\sum_{i=1}^N \int_{u\in{\mathbb C}_i}\Lambda_i(u)\Lambda_i(u)^{^\intercal}du\right], \] ${\boldsymbol V}_{NT}^\dagger={\sf diag}\{\widehat{\lambda}_1,\cdots,\widehat{\lambda}_q\}$ with $\widehat\lambda_j$ being the $j$-th largest eigenvalue of $\frac{1}{T-1}\widehat{\boldsymbol\Omega}$, and ${\boldsymbol\xi}=\left(\xi_2,\cdots,\xi_T\right)^{^\intercal}$. Evidently, the mean square convergence rate in ((ref)) is slower than that of ((ref)). Using the first restriction in ((ref)), the factor loading functions are estimated as

via least squares. Furthermore, we can estimate the (original) cointegrated factors by

We next discuss estimation of $q$ and $q^\ddag$. The information criterion ((ref)) needs modification to consistently estimate $q$. Specifically, we define

where $\rho_{NT}^\dagger$ satisfies \[ \rho_{NT}^\dagger\rightarrow0,\quad q\left(N^{-1/2}+\delta_q\right)+T^{-1/2}=o\left(\rho_{NT}^\dagger\right), \] which differ from ((ref)) in Theorem (ref). It follows from Proposition 5.1 in LLSX24 that ${\sf P}(\widehat q=q)\rightarrow1$. To estimate the cointegrating rank $q^\ddag$, we adopt the information criterion introduced by CP09, CP12 and modified for a curve time series context in P25, which is robust to weak dependence and time-varying variances in the errors. Assume the following VECM structure:

where $\alpha_0$ and $\beta_0$ are two $q\times q^{\ddag}$ matrices, and $v_t$ is stationary satisfying the conditions in CP09 or heterogeneously distributed as assumed in CP12. For each $j = 1, \cdots, \widehat q-1$, we estimate the $\widehat{q}\times j$ matrices $\alpha_0$ and $\beta_0$ via reduced-rank regression (RRR), giving $\widehat\alpha(j)$ and $\widehat\beta(j)$, and subsequently define \[ \widehat{\boldsymbol\Sigma}(j)=\frac{1}{T-1}\sum_{t=2}^T\left[\widehat\xi_t-\widehat\alpha(j)\widehat\beta(j)^{^\intercal} \widehat G_t\right]\left[\widehat\xi_t-\widehat\alpha(j)\widehat\beta(j)^{^\intercal} \widehat G_t\right]^{^\intercal}, \] as the residual covariance matrix, with $\widehat{\boldsymbol\Sigma}(0)=\frac{1}{T-1}\sum_{t=2}^T\widehat\xi_t\widehat\xi_t^{^\intercal}$. Cointegrating rank is selected as

with $\rho_{T}^{\ddagger}=\log T$ corresponding to the Bayesian information criterion (BIC) and $\rho_{T}^{\ddagger}=2\log \log T$ corresponding to the HQ criterion HQ79. Limit theory and consistency for these criteria, as well as the inconsistency of the related AIC criterion, are provided in P25.

Simulations

\setcounter{equation}{0}

This section reports the findings of two simulation studies designed to examine the finite-sample performance of our proposed methods. In each example, we first assess the estimation performance given that the number of common stochastic trends (or the cointegrating rank) is known and then examine the performance of various information criteria defined in Sections (ref) and (ref). To quantify the assessment of functional PCA, we compute the approximation errors for factors and factor loadings as follows

where $\widetilde\xi_t=\widetilde G_t-\widetilde G_{t-1}$ for $t=2,\cdots,T$. Similarly, for functional PANIC estimation the measurements are defined as

{\bf Example 6.1}.\ \ Consider $G_t=\sum_{s=1}^t\xi_s$ with $\xi_t$ following a VAR(1) model given by \[ \xi_t = \mathbf{A}\xi_{t-1} +\epsilon^{\xi}_t, \] where $\mathbf{A}$ is a $q\times q$ diagonal companion matrix and $\epsilon^{\xi}_t$'s denote the innovations. As in TNH23b, the diagonal entries of $\mathbf{A}$ were randomly drawn from a uniform distribution ${\cal U}[-1,1]$ and the matrix rescaled to have operator norm $0.8$. The innovations were independently drawn from a $q$-variate standard normal distribution. The initial 100 observations of the VAR(1) process $\{\xi_t\}$ were discarded to ensure data stability and independence of initial conditions. Letting $\phi_1,\cdots,\phi_{51}$ be $51$ orthonormal basis functions on $[0,1]$, we generated \[ \eta_{t1}(u)=\sum_{j=1}^{51}b^{\eta}_{tj}\phi_j(u)/j^2,\quad t=1,\cdots,T, \] where \[ b^{\eta}_{tj}=\frac{1}{\sqrt{T}}\left[\sum_{s=1}^t\epsilon_{tj}^{\eta}-\left(\frac{t}{T}\right)\sum_{s=1}^T\epsilon_{tj}^{\eta}\right] \] and the $\epsilon_{tj}^{\eta}$'s were independently drawn from the standard normal distribution. Latent factor curves were generated via (ref) with $k=1$ and ${\Phi}_1(u)=[\phi_1(u),\cdots,\phi_q(u)]^{^\intercal}$, i.e.,

where the factor $1/q$ in the series approximation errors serves the purpose of enhancing the signal-to-noise ratio. Factor loading functions were simulated as

where the $b_{i,j_1j_2}$'s were independently generated from the uniform distribution ${\cal U}[0,3]$ over $i$, $j_1$ and $j_2$.

The functional idiosyncratic components $\varepsilon_{it}(u)$ were generated by

where, for each $t$, $\mathbf{b}_t^{\varepsilon}=(b^{\varepsilon}_{1t,1},b^{\varepsilon}_{1t,2},\cdots,b^{\varepsilon}_{Nt,51})^{^\intercal}\in {\mathbb R}^{51N}$ were independently drawn from ${\cal N}({\bf 0},\ {\boldsymbol\Sigma}_{b})$ with $\ {\boldsymbol\Sigma}_{b}$ a block covariance matrix with ($i,j$)-block \[ \Sigma_{b,ij}=\max\{0,(1-|i-j|/10)\}\cdot {\sf diag}(1^{-2},\cdots,51^{-2}),\quad 1\leq i,j\leq N. \]

Combining (ref)--(ref), the nonstationary curve observations $Z_{it}$ were generated as $$

$$ where

Table (ref) reports the logarithms of the approximation errors for common stochastic trends obtained by functional PCA and PANIC, i.e., $\log({\sf AE}(\widetilde G))$ and $\log({\sf AE}(\widehat G))$, respectively. In Table (ref) for functional PCA estimation performance, a consistent reduction in the approximation errors is observed with an increase in $N$ across all values of $T$ and $q$. A slight increase in the approximation errors is noted with an increase in $T$ but the magnitude of increase is insignificant compared with standard deviations. Focusing on the three diagonal cases in each block of the table, i.e., ($N=100$, $T=200$), ($N=200$, $T=300$), and ($N=300$, $T=400$), we observe that the approximation errors diminish as both $T$ and $N$ approach infinity, confirming the joint convergence of functional PCA (see Proposition (ref)). In addition, as $q$ increases, the logarithms of the approximation errors decrease generally, which may be partly attributed to dominant loadings on the first few basis functions in the definitions of $\varepsilon_{i t}(u)$ and $\Lambda_i(u)$, see (ref) and (ref), and reduction of the sieve approximation errors (when $q$ increases). The results of functional PANIC follow a similar pattern, as evident in Table (ref). The approximation errors decrease as $N$ increases, but are insensitive to increasing $T$. Although the functional PANIC estimates converge when $N$ and $T$ jointly diverge to infinity, they generally exhibit higher approximation errors than functional PCA, indicating that the latter is more efficient when $G_t$ is of full rank.

Table (ref) reports logarithms of approximation errors for factor loading functions obtained through functional PCA and PANIC, i.e., $\log({\sf AE}(\widetilde \Lambda))$ and $\log({\sf AE}(\widehat \Lambda))$. In Table (ref) for functional PCA, a consistent decrease is evident in the approximation errors as $T$ increases across all combinations of $N$ and $q$. In contrast when $N$ varies the approximation errors remain relatively stable. This reversal of roles between $N$ and $T$ in comparison to Table (ref) reveals an interesting pattern, which was also observed in TNH23b for stationary curve time series. When $N$ and $T$ are fixed, the approximation errors increase as $q$ increases, which differs from the evolving pattern observed in Table (ref). In Table (ref) for functional PANIC, the approximation errors decrease when $T$ and $N$ increase. As in Table (ref), functional PANIC also exhibits higher approximation errors than functional PCA, which again shows that functional PCA converges faster than functional PANIC in this example.

Table (ref) reports the numbers for underestimation (in square brackets), correct-estimation, and over-estimation (in round brackets) for $q$ (the number of full-rank stochastic trends) over $1000$ replications. For functional PCA the number of stochastic trends is correctly estimated in most trials. The underestimation numbers are zero across all combinations of $(N,T,q)$, whereas overestimation numbers are generally small (i.e., $<3\%$). For functional PANIC, the correct estimation numbers are again close to $1000$. The underestimation numbers are zero except when $q$ is large but $N$ and $T$ are small ($q=15$, $N=100$, $T=200$). The overestimation numbers are small, decreasing rapidly when $N$ and $T$ increase. These outcomes suggest that the two information criteria ((ref)) and ((ref)) perform accurately in determining the number of common stochastic trends when either functional PCA or PANIC is adopted.

{\bf Example 6.2}.\ \ The next example has $G_t$ generated from the VECM (ref), where $\alpha_0$ and $\beta_0$ are $4\times q^{\ddag}$ matrices to be defined later, and $v_t$ follows a VARMA(1,1) process\footnote{The initial 100 observations of the VARMA(1,1) process are discarded to ensure data stability over time.}: \[ v_t=0.4v_{t-1}+\epsilon^{v}_t+0.4\epsilon^{v}_{t-1} \] with $\epsilon^{v}_{t}$ independently drawn from ${\cal N}({\bf0},{\sf diag}(1.25,0.75,1.4,0.6))$. Similar to CP09, we consider the following four scenarios for $(\alpha_0,\beta_0,q^\ddag)$:

where ${\boldsymbol O}$ is a $4\times 4$ null matrix, \[ {\boldsymbol R}_1=

,\quad {\boldsymbol R}_2=

,\quad {\rm and}\quad {\boldsymbol R}_3=

. \] The functional idiosyncratic components are generated by $\varepsilon_{it}=\sum_{s=1}^t\varepsilon_{is}^\dagger$ with $\varepsilon_{it}^\dagger$ simulated according to ((ref)) where $\varepsilon_{it}(u)$ is replaced by $\varepsilon_{it}^\dagger(u)$. Finally, we generate the nonstationary curve observations: \[ Z_{it}(u)=\Lambda_i(u)^{^\intercal}G_{t}+\varepsilon_{i t}(u), \] where $\Lambda_i(u)$ is generated in the same way as in (ref) with $q=4$.

Table (ref) reports logarithms of approximation errors for the stochastic trends in levels obtained by functional PCA and PANIC. In Table (ref) for functional PCA, we observe a significant increase in approximation errors with the expansion of the time series length ($T$) across all values of $N$ and $q^{\ddagger}$. When $N$ increases, the pattern for approximation errors is not the same in different settings. Focusing on the three diagonal cases in each block of the table, i.e., ($N=100$, $T=200$), ($N=200$, $T=300$) and ($N=300$, $T=400$), we observe that the approximation errors still increase as both $T$ and $N$ diverge, suggesting that functional PCA estimates are inconsistent due to violation of the full-rank condition in Assumption (ref). In Table (ref) for functional PANIC, we observe an increase in approximation errors when $T$ increases, but decreases in approximation errors when $N$ increases. Furthermore, the decreasing approximation errors in the three diagonal cases as $N$ and $T$ increase indicates that functional PANIC consistently estimates the stochastic trends in levels.

Table (ref) reports logarithms of approximation errors for the stochastic trends in differences obtained by functional PCA and PANIC. Functional PCA, as observed in Table (ref), is inconsistent when $q^{\ddagger}>0$. However, when $q^{\ddagger}=0$, the approximation errors decrease as $N$ grows but are stable with respect to $T$, and consequently decrease as both $N$ and $T$ tend to infinity. Functional PANIC, reported in Table (ref), has decreasing approximation errors with expansion of $N$ across all values of $T$ and $q^\ddag$. The approximation errors slightly increase as $T$ increases. Focusing on the three diagonal cases in each block of the table, we observe that the approximation errors generally diminish as both $T$ and $N$ increase, suggesting that functional PANIC can consistently estimate increments of the stochastic trends.

Table (ref) reports logarithms of approximation errors for the functional factor loadings, i.e., $\log({\sf AE}(\widetilde \Lambda))$ and $\log({\sf AE}(\widehat \Lambda))$. For the functional PCA results in Table (ref), the patterns of approximation errors evolving with $N$ and $T$ observed within each block of the table, indicate that the factor loading estimates via functional PCA are inconsistent. This may be due to inconsistency of the functional PCA in estimating the cointegrated factors. For the functional PANIC results in Table (ref), the approximation errors decrease as $T$ increases and remain stable when $N$ varies. Consequently, the approximation errors via the functional PANIC decrease as $N$ and $T$ jointly diverge.

Table (ref) reports numbers of underestimation, correct-estimation, and overestimation of cointegrating ranks under the BIC and HQ criteria. When the sample size is small ($N=100$ or $200$ and $T=200$ or $300$) and $q^{\ddag}=2$ or 3, BIC tends to underestimate cointegrating rank and thereby choose more parsimonious models. This observation aligns with the finding of CP09. It is worth pointing out that $N$ and $T$ play different roles in cointegrating rank estimation. An increase in $N$ results in reduced approximation errors for the stochastic trends and consequently increases the numbers of correct estimation (of $q^\ddag$) in most scenarios. In contrast, an increase in $T$ leads to larger approximation errors, as seen in Table (ref), but simultaneously contributes to more accurate cointegrating rank estimation (when the cointegrated factors were known), which is assured by theorems in CP09,CP12. When $T$ increases to $400$, the performance of BIC improves significantly. It follows from Table (ref) that the HQ criterion exhibits a notable tendency to over-estimate cointegrating rank, especially when $q^{\ddag}$ is small. The HQ criterion outperforms BIC only when $q^{\ddag}=3$.

Real data analysis

\setcounter{equation}{0}

This section presents two empirical applications of our methods. The first example studies a dataset of temperature curves and functional PCA is used. The second example studies a dataset of stock price curves and functional PANIC is employed. A notable distinction between the two applications lies in the properties of the nonstationary idiosyncratic components. For the temperature data collected from different weather stations, it is less likely for the idiosyncratic components to be nonstationary, given that temperature records at individual locations generally do not deviate significantly from the global or regional temperature patterns -- see OU03 and the references therein. In contrast, for the stock price data, the presence of nonstationary idiosyncratic components is highly probable due to the dynamic nature of financial markets, influenced by the occurrence of firm-specific information and announcements. It would be unrealistic to expect stock prices to exhibit a nonstationary factor structure with stationary idiosyncratic components. Such scenarios could create numerous hedging opportunities among randomly selected stocks, a phenomenon not observed in real-world financial markets.

Minimum temperatures in Australia

We applied functional PCA to yearly minimum temperature curves in Australia. The initial dataset collected from the Australian Bureau of Meteorology at \url{http://www.bom.gov.au} comprised daily minimum temperature observations. This dataset was considered by ARS18 and NSS23 in studying nonstationarity in temperature dynamics. For each calendar year, we employed a smoothing algorithm of Rpackagefda using 51 Fourier basis functions to create a curve representing minimum temperatures throughout the year\footnote{This smoothing process was applied only if the number of observations exceeded 200 days in a year. Otherwise, the data for the weather station was removed for that year.}. To deal with missing observations in an entire year, we adjust the calculation of $\widetilde{\Omega}_{ts}$ in (ref), only using weather stations that have observations in both years $t$ and $s$. Our analysis focuses on weather stations that initiated observations before 1943 and continued beyond 2022. Consequently, we have 36 stations with 80 years of observations and $Z_{it}(u)$ denotes the minimum temperature of weather station $i$ on year $t$ at day $u$.

Figure (ref) shows estimates of the stochastic trends and their loadings. For illustration, the estimates of stochastic trends were scaled by the corresponding eigenvalues, while the loading functions were normalized by dividing them by the corresponding eigenvalues. Determination of the number of stochastic trends was based on the information criterion ((ref)), which resulted in two common stochastic trends. The first stochastic trend reveals an upward trajectory, and its loading functions depict a temperature profile characteristic of Australia, with higher temperatures observed at the beginning and end of the year. This stochastic trend, consistent with findings in several other studies such as NSS23, signifies a stochastic trend in the mean temperature, suggesting the presence of global warming. The second stochastic trend exhibits a substantial negative value in 1943 and deviates from zero during the period from 1973 to 1993, suggesting significant temperature fluctuations within those years. Its loading functions display diverse patterns across different weather stations, highlighting their ability to capture station-specific intra-year temperature dynamics.

High-frequency stock prices of S&P 500 index constituents

We next applied functional PANIC to intraday log-prices of S&P 500 stocks. We selected the time period from 3 January 2023 to 1 November 2023, containing 209 trading days after removing a half trading day on 3 July 2023. The sample included $N = 209$ stocks. We adopted the 5-min frequency rather than 1-min frequency in data collection to minimize the impact of microstructure noise effects. Since all stocks trade from 9:30 a.m. to 4:00 p.m., 79 measurements were available per day. Asynchronous missing observations were interpolated by the linear algorithm of Rpackageforecast. The discrete data were converted to a continuous function using Rpackagefda's algorithm, and the resulting curves denoted by ${Z_{1t}(u), Z_{2t}(u),\cdots, Z_{Nt}}(u)$, with $N=209$ and where the index $u$ lies in the time interval between 9:30 a.m. and 4:00 p.m.

The scree plot in Figure (ref) shows the first 50 sample eigenvalues in log-scale. The first three eigenvalues are relatively large, leading to the selection of three stochastic trends based on the proposed information criterion. The cointegrating rank determined by BIC is zero, signifying the absence of cointegration among the estimated stochastic trends. In fact, the lack of a cointegrating relation aligns with the expectation of an efficient market, where hedging opportunities arising from such a relation should not exist.

Further, increments of the three estimated stochastic trends were regressed on the Fama-French five factors\footnote{Data are collected from \url{https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/Data_Library/f-f_factors.html}.}, i.e., market ($r_m-r_f$), size (SMB), value (HML), profitability (RMW), and investment (CMA). The results are reported in Table (ref). For the first stochastic trend, the market factor is significant, whereas for the second stochastic trend, both the market and size factors are significant at the 10% level. No significant factors are identified for the third stochastic trend. The relatively low $R^2$ values indicate that the Fama-French factors may not fully explain the movements of the three common stochastic trends.

Conclusion

\setcounter{equation}{0}

The emergence and growth of vast cross section and time series datasets has substantially increased interest in the development of high-dimensional methods in econometrics. This paper contributes to this growing body of literature by introducing a general dual functional factor model for large-scale nonstationary curve time series. The approach involves the construction of a high-dimensional factor model for the observed curve time series that allows both factors and factor loadings to lie in function spaces with a low-dimensional factor model structure obtained by way of sieve approximation. An important feature of this framework is that both the dimension and time series length diverge to infinity. For the case of full-rank integrated factor curves and stationary functional idiosyncratic components, we employ functional PCA methodology to estimate the common stochastic trends and functional factor loadings and establish mean square convergence and asymptotic distribution theory. An easy-to-implement information criterion is proposed to consistently select the number of common stochastic trends. A functional PANIC methodology is introduced to handle the more general setting with cointegrated factors and possibly nonstationary functional idiosyncratic components. The simulation results reveal that functional PCA outperforms functional PANIC when factors are full-rank integrated and functional idiosyncratic components are stationary, whereas functional PANIC is more reliable when the integrated factors are rank-reduced and functional idiosyncratic components are nonstationary. Two empirical case studies are provided from climatological and financial data, each demonstrating the existence of common stochastic trends for these high dimensional curve time series.