Inference on common trends in functional time series

-1cm Inference on common trends in functional time series

Introduction

We consider statistical inference on unit roots and common stochastic trends for time series taking values in a Hilbert space of arbitrary dimension or a subspace of possibly unknown dimension. An important first step in the analysis of such time series is the determination of the dimension of the nonstationary subspace, i.e., the subspace in which the time series behaves like a unit root process (see Section (ref)). This dimension is the number of common stochastic trends. Our objectives are to test hypotheses on the dimension of the nonstationary subspace as well as hypotheses on the stationary and nonstationary subspaces themselves.

Since we want to examine nonstationary time series in a possibly unknown-dimensional subspace of a Hilbert space $\mathcal H$, our tests need to be statistically valid regardless of (i) whether the dimension of the space in which the time series takes values is finite or not and (ii) whether the dimension is known in advance or not. That is, the tests need to be (at least asymptotically) invariant to the dimensionality of the time series, and due to this property they will be called Asymptotically Dimension Invariant (ADI) tests.

This is especially relevant in the recent literature on functional time series, where it is both empirically and theoretically supported that nonstationarity tends to be driven by a finite-dimensional process Chang2016,BS2018,Franchi2017b,LRS2020,LRS2020nonst,NSS. Even when a function-valued random element can be accommodated in a finite-dimensional space, and thus can be represented by a finite number of basis functions, its dimension is generally unknown and large. In such cases, most existing cointegration rank tests are not applicable since they require either finite (preferably small) dimensionality or a priori information on the dimensionality of the time series; this is true even for recently developed cointegration tests for high-dimensional time series Onatski2018,BG2022,BG2023.

Our testing procedures have wide applicability in practice. For example, they can be used with (i) cointegrated vector time series of finite dimension SW1988,Johansen1991, (ii) high-dimensional factor models with a finite number of nonstationary factors NelsonSiegel87,PENA2004291,PENA20061237, (iii) cointegrated curve-valued (or function-valued) time series Chang2016,NSS, and (iv) nonstationary dynamic functional factor models martinez2020nonparametric.

\phantomsectionOur approach to developing ADI tests utilizes a dimension-reduction method that projects a high-dimensional or functional time series onto a fixed $\mathrm{K}$-dimensional subspace. This projection is constructed such that, asymptotically, the $\mathrm{K}$-dimensional projected component includes all stochastic trends, while the residual component (of unknown or possibly infinite dimension) contains no information regarding these trends. By constructing test statistics solely from the projected series, we effectively bypass interference from the residual component, ensuring that the resulting inference is inherently ADI. The analysis proceeds by examining generalized eigenvalues of two variance operators associated with the projected time series. Thus, the proposed tests will be classified as variance ratio-type tests, and some are generalizations of existing tests (Remark (ref) and Appendix (ref)). Variance ratio tests have desirable properties for the study of nonstationary time series. First, they can avoid estimation of the long-run variance, which implies a consistency property not shared by other tests muller2007,muller2008. Second, the limiting behavior of the tests does not depend on a parametric assumption, such as a vector autoregression. As shown in the finite-dimensional case, finite-sample properties of parametric cointegration rank tests, such as those proposed by SW1988,Ahn1990,Johansen1991,Bewley1995,Ahn1997, depend crucially on the model specification; see Toda1995,Haug1996,Bewley1998. \phantomsectionThird, many parametric assumptions, e.g.\ autoregressive structures, are not preserved under projection, and tests relying on such assumptions cannot generally be ADI (Remark (ref)). In these regards, nonparametric variance ratio tests are appealing.

Our theoretical results can be summarized by the following four points. First, we provide limit theory for general variance ratio-type statistics based on partial summation and/or differencing, and this limit theory is applied to obtain our ADI inferential methods. Second, we apply our ADI tests sequentially to determine the dimension of the nonstationary subspace via either a top-down, bottom-up, or hybrid approach. The bottom-up does not require the choice of an initial hypothesis, but the top-down has better finite-sample properties, while the hybrid combines these advantages. Third, we consider also a direct estimator of the dimension of the nonstationary subspace based on ratios of generalized eigenvalues. This is similar to estimation of the dimension of the “dominant subspace” in LRS2020,LRS2020nonst and may be considered complementary to their estimator for curve-valued time series and to those in Zhang2018, Zhang_et_al2019, and Franchi2023 for vector-valued time series. Fourth, we consider hypothesis testing on the nonstationary or stationary subspaces.

Comparisons with existing work

First, we view the time series of interest, $\{X_t\}_{t\geq 1}$, as a sequence in a known Hilbert space $\mathcal H$, but allow for the possibility that $\{X_t\}_{t \geq 1}$ takes values only in a subspace of $\mathcal H$. In the latter case, the variance operator of $X_t$ is singular on $\mathcal H$ and may allow only finitely many nonzero eigenvalues even in an infinite-dimensional setting. This violates the common assumption in the literature of having sufficiently many (often infinitely many) nonzero eigenvalues. Of course, the eigenstructure of such a population variance operator is not known in practice, and assuming that infinitely many eigenvalues are nonzero may not only be unrealistic but is also not testable. For example, each functional observation of interest could be constructed by only a finite number of discrete and regularly spaced points on a grid of the domain, and such empirical examples can easily be found in the literature; see e.g., NSS (age-specific employment rates), LRS2020nonst (US treasury yield curves, in the working paper version), and Section (ref). This poses an obstacle to the use of existing methods, and this is precisely where our ADI inference methods have a distinctive advantage over existing ones. As will be detailed, our methods are designed to be asymptotically valid under a more general scenario where the (long-run) variance operator of the stationary part of $X_t$ permits only a few nonzero eigenvalues (we also provide analysis on the exact number of nonzero eigenvalues required for the proposed tests). In comparison, existing methods rule out this possibility by assumption (including our own recent work, NSS). Thus, our methodology can accommodate more general and realistic functional time series.

Second, the common requirement of infinitely many nonzero eigenvalues for existing methods distinguishes them from cointegration rank tests developed in a finite-dimensional setting, where any variance operator necessarily permits only finitely many nonzero eigenvalues. We show a natural connection between our ADI tests and well-known tests for cointegration (or stationarity) developed in a finite-dimensional setting, and we find that our ADI tests generalize those tests. Furthermore, we demonstrate that, in specific cases, our ADI tests reduce to some recently developed tests for nonstationary functional time series. These results imply that some well-known tests, developed in both conventional and functional setups, can be understood as special cases of our ADI tests, thereby bridging the gap between them.

Third, the advantages of our methods lie not only in their applicability to a more general setting but also in their practical usefulness for real datasets. Many existing methods require a reasonable conjecture on the maximum possible number of stochastic trends Chang2016,LRS2020nonst,NSS. Although some rule-of-thumb choices can be used in practice, such as those based on the scree plot, a more formal procedure may be preferable. The current study complements this aspect by providing a testing procedure based on a generalized version of the functional KPSS test Horvath2014 that can be used in practice to construct an upper bound.

Fourth, we provide statistical inference on the space spanned by stochastic trends based on our ADI tests. Specifically, we can test whether a particular subspace of interest is contained in the nonstationary subspace, spans the nonstationary subspace, or is contained in the stationary subspace. We illustrate the empirical usefulness of these tests in Section (ref) in the context of the well-known NelsonSiegel87 model.

Compared with our earlier work in NSS, all four preceding points apply (see Appendix (ref)). Also, in NSS we only considered a particular variance ratio statistic (VR(2,1) under more restrictive assumptions than in Section (ref) below), whereas in this paper we consider a family of statistics as well as eigenvalue ratio estimators.

Compared with recent work on cointegration testing in high-dimensional vector autoregressive (VAR) processes Onatski2018,BG2022,BG2023, our setup is fundamentally different. They allow only a finite number of cointegrating vectors, whereas our setup requires that the number of stochastic trends is finite and hence is more suitable when nonstationarity is driven by a small number of factors. Interestingly, under the assumption of compact autoregressive operators widely adopted in the literature, a function-valued VAR implies that the number of stochastic trends is finite (Remark (ref)). Moreover, our setup does not rely on a VAR structure, and it covers intrinsically infinite-dimensional time series which theirs does not. Finally, our tests can be naturally applied to determine the number of stochastic trends while their tests are for the existence of cointegration and cannot immediately be applied as cointegration rank tests. In general, therefore, their tests serve a different purpose and may be viewed as complementary; see also Remark (ref).

\phantomsectionIn concurrent work, LiLiPhillips2025 investigated common stochastic trends within a more complex framework involving a growing number of curve-valued time series. Such data, often referred to as high-dimensional functional time series (e.g., Tang2025), can be understood as an integration of high-dimensional time series analysis and functional data analysis. They proposed an approximate factor model for such data and developed inferential methods, including a consistent estimator for the low-dimensional stochastic trends that drive the nonstationarity of the high-dimensional functional series. While the present setup does not accommodate their framework, the exploration of a possible extension of our methodology to their setting may be an interesting avenue for future research.

The remainder is organized as follows. Section (ref) introduces I(1) time series in Hilbert space, and Section (ref) presents our ADI variance ratio tests. Section (ref) discusses determination of the dimension of the nonstationary subspace via sequential testing or via eigenvalue ratio estimation. In Section (ref) we discuss testing hypotheses about the subspaces. Section (ref) presents Monte Carlo simulations, and Section (ref) presents two empirical applications to the term structure of interest rates and to labor market indices. Finally, Section (ref) concludes. The appendix contains mathematical notation, details, and proofs, and a supplementary appendix includes a practical implementation guide, additional discussion, simulations, and proofs.

I(1) time series and stochastic trends in Hilbert space

\phantomsectionWe consider a cointegrated linear I(1) process $X_t$ taking values in a separable Hilbert space, $\mathcal H$, equipped with an inner product $\langle v_1, v_2 \rangle$ and the induced norm $\|v_1\| = \langle v_1,v_1 \rangle^{1/2}$ for $v_1,v_2\in \mathcal H$, as defined by Beare2017. For stationary time series in $\mathcal H$, the real-valued sequence $\langle X_t,v \rangle$ must be stationary for all possible choices of $v \in \mathcal H$ (see Proposition 3.1 and its proof in Beare2017). For example, if $\mathcal H$ is the standard $L^2[a,b]$ Hilbert space and hence $X_t(r)$ is a function defined on $[a,b]$, then $\langle X_t,v \rangle = \int_a^b X_t(r)v(r)dr$ may be viewed as a continuous linear combination of $X_t(r)$. For many empirical applications involving economic functional time series, stationarity is too restrictive Chang2016,SEO2019,NSS. Thus, we consider a sequence of $X_t$ that allows $\langle X_t, v\rangle$ to be either I(1) nonstationary or I(0) stationary. It is then important to characterize under what choices of $v\in \mathcal H$ that $\langle X_t, v\rangle$ becomes stationary or nonstationary, which leads to notions of stationary and nonstationary subspaces. We next discuss this in more detail.

Let $\{ X_t \}_{t \geq 1}$ be a sequence whose first difference, denoted $\Delta X_t$, satisfies

where $\{ \epsilon_t \}_{t \in \mathbb{Z}}$ is an independent and identically distributed (iid) sequence with $\operatorname{\mathbb{E}} [\epsilon_t] = 0$, $\operatorname{\mathbb{E}} [\|\epsilon_t\|^{4}] < \infty$, and positive definite variance $C_\epsilon$. We further assume that $\{\Phi_j \}_{j \geq 0}$ is a sequence of bounded, linear operators satisfying

where $\Vert \cdot\Vert_{\operatorname{op}}$ denotes the operator norm. Then the long-run variance of $\Delta X_t$ is well defined as $\Lambda_{\Delta X} = \Phi(1) C_\epsilon \Phi (1)^\ast \neq 0$; see Beare2017. Any stationary sequence with nonzero long-run variance is I(0). Hence, because $\{\Delta X_t\}_{t\geq1}$ satisfying \hyperref[{eqlinear}]{\tagform@{\ref*{eqlinear}}} and \hyperref[{eqlinear2}]{\tagform@{\ref*{eqlinear2}}} is necessarily stationary and $\Lambda_{\Delta X}\neq 0$, it is I(0) by construction, so that $\{ X_t \}_{t \geq 1}$ is I(1). The summability condition in \hyperref[{eqlinear2}]{\tagform@{\ref*{eqlinear2}}} is common in the unit root literature PS1992.

Let $\mathcal H_{\mathsf{S}} = \ker \Lambda_{\Delta X}$ denote the kernel of $\Lambda_{\Delta X}$ and let $\mathcal H_{\mathsf{N}}=\mathcal H_{\mathsf{S}}^\perp$ be its orthogonal complement. We call $\mathcal H_{\mathsf{S}}$ and $\mathcal H_{\mathsf{N}}$ the stationary (cointegrating) subspace and nonstationary (attractor) subspace, respectively. These names are related to some distinctive properties possessed by elements in those spaces.

First, it follows from Beare2017 that $\langle X_t , v \rangle$ is stationary for $v \in \mathcal H_\mathsf{S}$ and nonstationary for $v \notin \mathcal H_\mathsf{S} = \mathcal H_{\mathsf{N}}^{\perp}$. In fact, we will allow $X_t$ to take values only in a strict subspace of unknown dimension. For all $v$ that are in the orthogonal complement of the strict subspace it holds that $\langle X_t,v\rangle =0$, and hence also $\langle \Delta X_t,v\rangle =0$ with long-run variance $\langle v, \Lambda_{\Delta X} v \rangle =0$. Because $\Lambda_{\Delta X}$ is self-adjoint and nonnegative, this implies that $v \in \ker \Lambda_{\Delta X}$ and thus $v \in \mathcal H_{\mathsf{S}}$. In general, it therefore follows that if $v \in \mathcal H_{\mathsf{S}}$ then $\{\langle X_t, v\rangle\}_{t \geq 1}$ is either a stationary random sequence or equal to zero almost surely (in which case it is also stationary).

Second, regarding $\mathcal H_{\mathsf{N}}$, we will assume throughout that

If \hyperref[{eqfinite}]{\tagform@{\ref*{eqfinite}}} is true, we say that $\{X_t\}_{t\geq1}$ contains ${\mathbbm{d}_{\mathsf{N}}}$ stochastic trends. To see this in detail, we note that $\{X_t\}_{t\geq 1}$ satisfying \hyperref[{eqlinear}]{\tagform@{\ref*{eqlinear}}} and \hyperref[{eqlinear2}]{\tagform@{\ref*{eqlinear2}}} allows the decomposition

where $\nu_t = \sum_{j=0} ^\infty \tilde \Phi_j \epsilon_{t-j} $ and $\tilde \Phi_j = -\sum_{s=j+1} ^\infty \Phi_s$; see PS1992. From the summability condition in \hyperref[{eqlinear2}]{\tagform@{\ref*{eqlinear2}}}, we also find that $\sum_{j=0}^\infty j \|\tilde \Phi_j\|_{\operatorname{op}}<\infty$. It is known Beare2017 that $\mathcal H_{\mathsf{N}}$ is the closure of $\operatorname{ran} \Phi(1)$. Thus, under \hyperref[{eqbn}]{\tagform@{\ref*{eqbn}}} we have

Since the random walk component in \hyperref[{eqbn}]{\tagform@{\ref*{eqbn}}} clearly takes values in $\mathcal H_{\mathsf{N}}$ and $\dim (\mathcal H_{\mathsf{N}}) = {\mathbbm{d}_{\mathsf{N}}}$, we say that $\{X_t\}_{t\geq1}$ contains ${\mathbbm{d}_{\mathsf{N}}}$ stochastic trends.

Since $\dim (\mathcal H_{\mathsf{N}}) = \dim (\operatorname{ran}\Phi(1)) = \operatorname{rank} (\Phi(1))$, the projection onto $\mathcal H_{\mathsf{N}}$, denoted $P_{\mathsf{N}}$, is a finite-rank operator of rank ${\mathbbm{d}_{\mathsf{N}}}$ under \hyperref[{eqfinite}]{\tagform@{\ref*{eqfinite}}}, regardless of whether $\dim (\mathcal H)$ is finite or not. However, the rank of the projection onto $\mathcal H_{\mathsf{S}}$, denoted $P_{\mathsf{S}} = I - P_\mathsf{N}$ with $I$ denoting the identity operator on $\mathcal H$, depends on $\dim (\mathcal H)$. Specifically, the rank is $p-{\mathbbm{d}_{\mathsf{N}}}$ if $\mathcal H = \mathbb R^p$, while it is $\infty$ if $\mathcal H$ is infinite-dimensional. Importantly, even in the latter case, if $X_t$ only takes values in a strict subspace of dimension $p_X$, then the (long-run) variance of $\{P_{\mathsf{S}}X_t\}_{t \geq 1}$ is not injective on $\mathcal H_{\mathsf{S}}$ and allows at most $p_X-{\mathbbm{d}_{\mathsf{N}}}$ nonzero eigenvalues, which may be only finitely many.

ADI variance ratio tests

We let $\mathcal H$ be a separable Hilbert space and let $\{X_t\}_{t\geq 1}$ take values in $\mathcal H$ or any subspace. $\mathcal H$ can be an infinite-dimensional space of square-integrable functions or sequences depending on whether $X_t$ is a random function or a random vector, but we do not exclude the possibility that $\mathcal H$ is finite-dimensional. We assume that $\{X_t\}_{t\geq 1}$ satisfies the following assumption.

A crucial input to estimation and inference on unit roots and cointegration Chang2016,NSS,seo2024functional is ${\mathbbm{d}_{\mathsf{N}}}=\dim(\mathcal H_{\mathsf{N}})$, i.e., the number of (linearly independent) stochastic trends. To determine this quantity, we consider the testing problem

for $\mathbbm{d}_0 \in \mathbbm{D}_0 = \{1,2,\ldots,\mathbbm{d}_{\max}\}$, where $\mathbbm{d}_{\max}<\infty$ is an upper bound on the number of stochastic trends so that $\mathbbm{D}_0$ is a finite set. In fact, $\mathbbm{d}_{\max}$ can be reasonably chosen from data in such a way that $\mathbbm{d}_{\max} \geq {\mathbbm{d}_{\mathsf{N}}}$, but $\mathbbm{d}_{\max}$ does not greatly exceed ${\mathbbm{d}_{\mathsf{N}}}$ (Remark (ref), Section (ref)).

As mentioned above, our tests can accommodate the case where $\{X_t\}_{t\geq1}$ only takes values in a strict subspace of possibly unknown dimension $p_X$, and the tests do not require knowledge of $p_X$. For this reason, they are called Asymptotically Dimension Invariant (ADI) tests. However, it should be noted that this dimension invariance property does not mean that our tests are not affected by $p_X$ at all. Rather, it means that we do not explicitly require a priori information on $p_X$ in order to implement the tests.

Our approach to developing ADI tests is based on a dimension-reduction technique that preserves stochastic trends. We call this slack extraction of stochastic trends. To elaborate on this, suppose that $\{X_t\}_{t\geq1}$ is a $\mathcal H$-valued cointegrated I(1) process as described in Section (ref) and only takes values in a possibly strict subspace of dimension $p_X \leq p_{\mathcal H} = \dim(\mathcal H) \leq \infty$, where $p_X$ may be much larger than ${\mathbbm{d}_{\mathsf{N}}}$. In addition, assume that we can find a finite integer $\mathrm{K} \in [ {\mathbbm{d}_{\mathsf{N}}} , p_X ]$ and an orthonormal set $\{ f_j\}_{j=1}^{\mathrm{K}}$ in $\mathcal H$ such that $\operatorname{span}\{ f_j\}_{j=1}^{{\mathbbm{d}_{\mathsf{N}}}}=\mathcal H_{\mathsf{N}}$ and $\{\langle X_t, f_j \rangle\}_{t\geq 1}$ is a stationary process for $j={\mathbbm{d}_{\mathsf{N}}} + 1 , \ldots , \mathrm{K}$. Define $P_{\mathrm{K}} = \sum_{j=1} ^{\mathrm{K}} f_j \otimes f_j$, where $\otimes$ denotes the tensor product and thus $f_j \otimes f_j (\cdot) = \langle f_j, \cdot \rangle f_j$; that is, $P_{\mathrm{K}}$ is the unique orthogonal projection onto $\mathcal H_{\mathrm{K}}=\operatorname{span}\{ f_j\}_{j=1}^{\mathrm{K}}$. In this case, we have $\mathcal H = \mathcal H_{\mathrm{K}} \oplus \mathcal H_{\mathrm{K}}^\perp$ and $ \mathcal H_{\mathrm{K}} = \mathcal H_{\mathsf{N}} \oplus \mathcal H_{\mathrm{K}}^{\mathsf{S}}$, where $\mathcal H_{\mathrm{K}}^{\mathsf{S}} = \operatorname{span}\{ f_j\}_{j={\mathbbm{d}_{\mathsf{N}}}+1}^{\mathrm{K}}$. Then, from the fact that $\{\langle X_t, f_j \rangle\}_{t\geq 1}$ is an I$(1)$ sequence for $j \leq {\mathbbm{d}_{\mathsf{N}}}$, we may deduce that (i) the projected time series $\{ P_{\mathrm{K}} X_t \}_{t\geq 1}$ contains ${\mathbbm{d}_{\mathsf{N}}}$ stochastic trends and a $(\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}})$-dimensional stationary component, and (ii) the nonstationary subspace of $\{ P_{\mathrm{K}} X_t \}_{t\geq 1}$ is the same as that of $\{ X_t \}_{t\geq 1}$. Due to these properties, we call $P_{\mathrm{K}}$ a slack extractor of stochastic trends. In other words, we do not lose any information about stochastic trends by using the projected time series $\{ P_{\mathrm{K}} X_t \}_{t\geq 1}$ in the statistical analysis and ignoring the residual time series $\{ (I-P_{\mathrm{K}})X_t\}_{t\geq 1}$ (whose dimension generally depends on $p_X$). Therefore, a statistical test, which is not dependent on the residual part $\{ (I-P_{\mathrm{K}}) X_t \}_{t\geq 1}$ but is constructed only from $\{ P_{\mathrm{K}} X_t \}_{t\geq 1}$, naturally becomes ADI. Based on this idea, one may readily develop ADI tests by, for example, extending existing tests developed in a finite-dimensional setting (because $\{ P_{\mathrm{K}} X_t \}_{t\geq 1}$ itself may be viewed as a $\mathrm{K}$-dimensional cointegrated time series, see Remark (ref)).

Of course, a slack extractor $P_{\mathrm{K}}$ satisfying the required conditions is not observable, but it can be replaced with a suitable estimator $\widehat{P}_{\mathrm{K}}$. For now, we impose some high-level conditions on $\widehat{P}_{\mathrm{K}}$ and its rank $\mathrm{K}$ (the latter is specified by the practitioner). Define the operator

where $1_{\{\cdot \}}$ is the indicator function, and $a_R$ is a parameter that will be naturally specified (in Assumption (ref) below) for each of the proposed tests. Note that $\Lambda_{\mathsf{S},R}$ is the long-run variance (if $a_R > 0$) or variance (if $a_R = 0$) of the stationary component of $\{X_{t}\}_{t\geq 1}$.

\newcounter{foo} \setcounter{foo}{\value{assumption}}

Assumption (ref)(i) implies that we need an appropriate choice of $\mathrm{K}$ depending on the value of ${\mathbbm{d}_{\mathsf{N}}}$. As will be discussed in Section (ref), see in particular Remarks (ref), (ref), and the subsequent discussion, it is possible to obtain a reasonable upper bound $\mathbbm{d}_{\max}$ for ${\mathbbm{d}_{\mathsf{N}}}$. By choosing $\mathrm{K}$ so that $\mathrm{K} \geq \mathbbm{d}_{\max}$, the condition on $\mathrm{K}$ required by Assumption (ref)(i) does not raise any issues in the practical implementation of our tests. For a given choice of $\mathrm{K}$, feasible choices of $\widehat{P}_{\mathrm{K}}$ satisfying Assumption (ref)(ii)--(iii) are discussed in Section (ref) under mild low-level conditions.

A consequence of Remark (ref) is that, even if we may have a reasonable parametric assumption such as an AR structure for the original time series $\{X_t\}_{t\geq 1}$, such an assumption is not generally preserved under the projection $\widehat{P}_{\mathrm{K}}$. Therefore, when we use the projected time series $\{\widehat{P}_{\mathrm{K}} X_t\}_{t\geq 1}$ for our statistical analysis to obtain ADI tests, these need to be developed without parametric assumptions. To achieve this goal, we consider tests based on various sample (long-run) variance operators that can be computed from the projected time series without any parametric assumptions.

Thus, define $X_{d,t}$ by

so that $d$ denotes the integration order of $X_{d,t}$. Further define the unnormalized sample long-run variance operator of $\{X_{d,t}\}_{t\geq 1}$, denoted $\widehat{\Lambda}_d (h,\mathrm{k})$, as

where $\mathrm{k}(\cdot)$ is a kernel function and $h$ is the associated bandwidth parameter. From \hyperref[{eqslrv}]{\tagform@{\ref*{eqslrv}}} and the fact that $\widehat P_{\mathrm{K}}$ is self-adjoint, it is readily found that $\widehat P_{\mathrm{K}} \widehat{\Lambda}_d (h,\mathrm{k}) \widehat P_{\mathrm{K}}$ is the unnormalized sample (long-run) variance operator of the projected time series $\{\widehat P_{\mathrm{K}} X_{d,t}\}_{t\geq 1}$. As a special case, $\widehat P_{\mathrm{K}} \widehat{\Lambda}_d (0,\mathrm{k}) \widehat P_{\mathrm{K}}$ is the unnormalized sample variance operator of $\{\widehat P_{\mathrm{K}} X_{d,t}\}_{t\geq 1}$.

Our tests are based on the generalized eigenvalue problem

where $\widehat{\Lambda}_{d_m,m} = \widehat{\Lambda}_{d_m} (h_m,\mathrm{k}_\mathit{m})$ for $m \in \{ L,R \}$ (left and right) and choices of integration orders, $d_L$ and $d_R$. Here, $\mu_j$ are the ordered (smallest to largest) eigenvalues of the operators $\widehat P_{\mathrm{K}}\widehat{\Lambda}_{d_L,L}\widehat P_{\mathrm{K}}$ and $\widehat P_{\mathrm{K}}\widehat{\Lambda}_{d_R,R}\widehat P_{\mathrm{K}}$, and $\nu_j$ are the corresponding eigenvectors. These can easily be computed as detailed in Section (ref) in the supplement. Because $\widehat P_{\mathrm{K}}\widehat{\Lambda}_{d_m,m} \widehat P_{\mathrm{K}}$ are sample (long-run) variance operators, we call \hyperref[{eqgev}]{\tagform@{\ref*{eqgev}}} a generalized variance ratio (VR) eigenvalue problem, and any test based on it is called a VR($d_L, d_R$)-based test. We assume $d_L > d_R$ without loss of generality because $d_L <d_R$ is obtained by simply redefining $\mu_j$ as its inverse.

In the sequel, $\mathrm{k}_\mathit{m}(\cdot)$ and $h_m$ are assumed to satisfy the following assumption.

Assumption (ref)(i) is adopted from phillips1995fully. The requirements on $\mathrm{k}_\mathit{m}(\cdot)$ in Assumption (ref)(i) are not restrictive in practice and are satisfied by many widely used kernel functions, including the Epanechnikov, Parzen, Tukey-Hanning, and quartic kernels. The functional form of $h_m$ depending on $T$ in Assumption (ref)(ii) is quite standard in practice.

Tests based on VR($d_L$,1) with $d_L \geq 2$

In this section we consider tests based on the VR($d_L$,1) eigenvalue problem in \hyperref[{eqgev}]{\tagform@{\ref*{eqgev}}} with $d_L \geq 2$ and $d_R = 1$. We first obtain the asymptotic properties of the eigenvalues of this problem.

To this end, we let $\{W_{1,q}(r)\}_{r\in[0,1]} $ denote a $q$-dimensional standard Brownian motion and recursively define the $(d-1)$-fold integrated Brownian motion $W_{d,q}(r) = \int_0^r W_{d-1,q} (u)du$ for $d \geq 2$. For any matrix or compact operator $A$, we let $\lambda_j \{ A \}$ be the $j$-th smallest eigenvalue of $A$, i.e.\ $\lambda_1\{A\}\leq \lambda_2\{A\} \leq \ldots$. Given a kernel $\mathrm{k}_\mathit{m}$ satisfying Assumption (ref), we define

Moreover, when there is no risk of confusion, we write $\int f$ to denote $\int_0^1 f(u)du$.

Note that the normalization factor \hyperref[{nT21}]{\tagform@{\ref*{nT21}}} of the eigenvalues generally depends on the choice of bandwidth parameters $h_m$ associated with the kernel $\mathrm{k}_\mathit{m} ( \cdot )$ for $m \in \{L,R\}$. However, in the case where $a_L=a_R=0$ (i.e., $h_L=h_R=0$), $\widehat{P}_{\mathrm{K}} \widehat{\Lambda}_{d_L,L} \widehat{P}_{\mathrm{K}}$ and $\widehat{P}_{\mathrm{K}} \widehat{\Lambda}_{1,R} \widehat{P}_{\mathrm{K}}$ reduce to the variances of $\{ \widehat{P}_{\mathrm{K}}X_{d_L,t}\}_{t\geq 1}$ and $\{ \widehat{P}_{\mathrm{K}}X_{1,t}\}_{t\geq 1}$, and the convergence rate is simply $n_T^{-1} = T^{2-2d_L}$.

Theorem (ref) shows that the eigenvalues in \hyperref[{eqgev}]{\tagform@{\ref*{eqgev}}} have distinct asymptotic properties depending on whether $j$ is greater than the dimension of $\mathcal H_{\mathsf{N}}$ or not, as motivated in Remark (ref). That is, if we consider the vector $(\widetilde\mu_1 ,\ldots, \widetilde\mu_{\mathrm{K}})$, the first ${\mathbbm{d}_{\mathsf{N}}}$ elements converge jointly in distribution to the vector $(\tilde{\lambda}_1 , \ldots, \tilde{\lambda}_{{\mathbbm{d}_{\mathsf{N}}}})$, while the last $\mathrm{K} - {\mathbbm{d}_{\mathsf{N}}}$ elements are divergent in probability. If $\mathrm{K} \geq \mathbbm{d}_0$, for any continuous map $\mathcal{F} : \mathbb{R}^{\mathbbm{d}_0} \to \mathbb{R}$ we have

If $\mathcal{F}$ additionally satisfies $\mathcal{F}(\{ \widetilde{\mu}_j\}_{j=1}^{\mathbbm{d}_0}) \rightarrow_{\mathrm{P}} \infty$ under $H_1$ of \hyperref[{eqhypo}]{\tagform@{\ref*{eqhypo}}}, we can consistently test the hypothesis of interest in an obvious way; hereafter such an $\mathcal{F}$ is called a proper test functional. Among many possible choices of proper test functionals, we focus on the two most common choices, $\mathcal{F}_{\max}$ and $\mathcal{F}_{\operatorname{trace}}$, which are defined by

The following corollary delivers consistent tests of \hyperref[{eqhypo}]{\tagform@{\ref*{eqhypo}}} based on VR($d_L$,1).

The test statistics in Corollary (ref) are only functions of $\{\widehat{P}_{\mathrm{K}}X_t\}_{t\geq1}$, which may be understood as a $\mathrm{K}$-dimensional vector-valued time series (Remark (ref)), and the associated limiting distributions are functionals of $\mathbbm{d}_0$-dimensional standard Brownian motion and do not depend on $p_X$, $\mathrm{K}$, or any nuisance parameters. These results imply that the proposed tests have the ADI property and that critical values depend only on $\mathbbm{d}_0$ and the test functional.

Tests based on VR($d_L$,0) with $d_L \geq 1$

We now consider VR($d_L$,0)-based tests, which are established from the eigenvalue problem \hyperref[{eqgev}]{\tagform@{\ref*{eqgev}}} with $d_L \geq 1$ and $d_R =0$. As with VR($d_L$,1)-based tests, the tests in this section are based on the asymptotic properties of the sample operators, $\widehat{P}_{\mathrm{K}}\widehat{\Lambda}_{d_L,L}\widehat{P}_{\mathrm{K}}$ and $\widehat{P}_{\mathrm{K}}\widehat{\Lambda}_{0,R}\widehat{P}_{\mathrm{K}}$ for $d_L \geq 1$. The important difference relative to tests based on VR($d_L$,1) is that VR($d_L$,0)-based tests generally require $a_R> 0$; that is, $\widehat{\Lambda}_{0,R}$ must be the sample long-run variance of $X_{0,t}$. This is in contrast with the case in Section (ref), where $a_R = 0$ can be chosen; see Remarks (ref) and (ref).

The asymptotic properties of the VR($d_L$,0) eigenvalues with $d_L \geq 1$ are as follows.

The next corollary delivers consistent tests of \hyperref[{eqhypo}]{\tagform@{\ref*{eqhypo}}} based on VR($d_L$,0) with $d_L \geq 1$.

Inverse VR tests

We showed in Theorem (ref) that the first ${\mathbbm{d}_{\mathsf{N}}}$ eigenvalues $\{\mu_j\}_{j=1}^{{\mathbbm{d}_{\mathsf{N}}}}$ from the VR($d_L$,1) problem with $d_L \geq 2$, properly normalized, converge jointly to a well-defined limit. As we now show, it is also possible to obtain the limiting behavior of the remaining $\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}$ eigenvalues, $\{\mu_j\}_{j={\mathbbm{d}_{\mathsf{N}}} +1}^{\mathrm{K}}$. This enables us to examine the inverse testing problem,

for some $\mathbbm{d}_0 \geq 0$. Unlike in the testing problem \hyperref[{eqhypo}]{\tagform@{\ref*{eqhypo}}}, we do not require any prior information on a reasonable upper bound for ${\mathbbm{d}_{\mathsf{N}}}$.

Related to testing \hyperref[{eqhypo2}]{\tagform@{\ref*{eqhypo2}}}, we will replace Assumption (ref) with the following high-level assumption, which allows the possibility that $\mathrm{K} < {\mathbbm{d}_{\mathsf{N}}}$.

\newcounter{foobar} \setcounter{foobar}{\value{assumption}} \setcounter{assumption}{\value{foo}}

\setcounter{assumption}{\value{foobar}}

Theorem (ref) complements Theorem (ref), where the asymptotic properties of the first ${\mathbbm{d}_{\mathsf{N}}}$ eigenvalues are presented under similar assumptions. Specifically, the asymptotic results given in Theorem (ref) are for ${\mathbbm{d}_{\mathsf{N}}} \geq 0$ under the requirement $a_R>0$, which differs from Theorem (ref), where the limiting distributions are given under the existence of unit roots ensuring ${\mathbbm{d}_{\mathsf{N}}} \geq 1$ but allowing $a_R = 0$. This means that Theorem (ref) is not simply a byproduct of Theorem (ref) that can be obtained with slight modifications, but in fact has its own theoretical and practical justification as we now illustrate.

When $\mathrm{K} > \mathbbm{d}_0$, and given a proper test functional $\mathcal{F}$, the asymptotic results in Theorem (ref) can be used to deliver consistent tests of \hyperref[{eqhypo2}]{\tagform@{\ref*{eqhypo2}}} as follows.

The test functionals described in Corollary (ref) are based on the limiting behavior of the inverse eigenvalues from the VR($d_L$,1) problem, so we call them inverse VR tests. Compared to Corollary (ref) concerning the testing problem \hyperref[{eqhypo}]{\tagform@{\ref*{eqhypo}}}, we allow $\mathbbm{d}_0=0$ in Corollary (ref). This means that we can, for example, test the null of stationarity of the time series $\{X_t\}_{t\geq 1}$ against an alternative of unit root nonstationarity. Hence, this test generalizes some existing, and widely used, KPSS-type tests of stationarity that are essentially obtained when $\mathbbm{d}_0=0$; see Appendix (ref) for a detailed discussion of related procedures.

Estimation of slack extractor

To implement the VR-based tests in practice, a slack extractor $\widehat{P}_{\mathrm{K}}$ satisfying either Assumption (ref) or (ref) is needed. As shown in Sections (ref)--(ref), both the VR and inverse VR tests require Assumption (ref) to hold for $\mathrm{K} \geq {\mathbbm{d}_{\mathsf{N}}}$, and when $\mathrm{K} < {\mathbbm{d}_{\mathsf{N}}}$ the inverse VR tests require instead that Assumption (ref) holds. In this section, we construct such estimators $\widehat{P}_{\mathrm{K}}$ from observations $\{X_t\}_{t\geq 1}$, and we propose low-level conditions on the eigenstructure of the operator $\Lambda_{\mathsf{S},R}$ defined in \hyperref[{eqlongrun}]{\tagform@{\ref*{eqlongrun}}}, under which the estimators satisfy the high-level conditions.

We let $\{\tau_{j}\}_{j \geq 1}$ be the eigenvalues of $\Lambda_{\mathsf{S},R}$, ordered from the largest to the smallest.

Assumption (ref) tells us that the allowable values of $\mathrm{K}$, which is chosen by the practitioner, depend on the number of nonzero eigenvalues of $\Lambda_{\mathsf{S},R}$ defined in \hyperref[{eqlongrun}]{\tagform@{\ref*{eqlongrun}}}. Specifically, it requires that the first $\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}$ eigenvalues of $\Lambda_{\mathsf{S},R}$ are nonzero and that the $(\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}})$-th eigenvalue is different from the next eigenvalue, $\tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}+1}$. The role of the requirement $\tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}} \neq \tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}+1}$ is subtle (and may be relaxed under some additional assumptions). It is employed to make the slack extractor $\widehat{P}_{\mathrm{K}}$ have a certain desirable property in our proof and can be checked using the eigenvalues of the sample counterpart of $\Lambda_{\mathsf{S},R}$; see Assumption (ref) and Theorem (ref).

In a high-dimensional Hilbert space, the number of nonzero eigenvalues of $\Lambda_{\mathsf{S},R}$ is generally very large (or possibly infinite). Particularly, in a typical functional time series setting with $\dim(\mathcal H) = \infty$ (e.g.\ the case considered by NSS), $\Lambda_{\mathsf{S},R}$ generally has infinitely many nonzero eigenvalues, and Assumption (ref) allows $\mathrm{K}$ to be any arbitrary finite integer as long as $\tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}} \neq \tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}+1}$. In a finite-dimensional space, a reasonable upper bound on ${\mathbbm{d}_{\mathsf{N}}}$ will be relevant to choose $\mathrm{K}$ to satisfy Assumption (ref).

We now describe the construction of the estimator of the slack extractor. Let

denote the sample variance (resp.\ long-run variance) estimator of $X_{1,t}$ if $a_R=0$ (resp. $a_R>0$), but allowing a new choice of bandwidth, $h_P$, and kernel, $\mathrm{k}_P(\cdot)$. The slack extractor needs to satisfy a rank condition in Assumption (ref) that depends on $\Lambda_{\mathsf{S},R}$, which is a variance when $a_R=0$ and a long-run variance when $a_R>0$. For this reason, we base our estimator of the slack extractor on $\widehat{\Lambda}_{1,P}$, which depends on $a_R$ in the same way; see also Remark (ref). The advantage is that the estimator $\widehat{P}_{\mathrm{K}}$ described in Theorem (ref) is universal in the sense that it can be applied for any of the VR($d_L,d_R$) tests that we consider.

Note that practitioners can employ the usual functional principal component analysis (FPCA) in $\mathcal H$ to obtain the eigenvectors in Theorem (ref). Thus, Theorem (ref) provides an easy-to-implement way to obtain an estimated projection operator $\widehat{P}_{\mathrm{K}}$ that satisfies the high-level conditions in Assumptions (ref) and (ref).

The condition $\tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}} \neq \tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}+1}$ in Assumption (ref) is generally required for the VR($d_L,d_R$)-based tests and also required for the inverse VR test when $\mathbbm{d}_0$ is close to ${\mathbbm{d}_{\mathsf{N}}}$. Practitioners might be concerned with potential violation of this condition for a chosen $\mathrm{K}$. However, as shown in our proof of Theorem (ref), the ordered eigenvalues $\{\widehat{\tau}_{j}\}_{j \geq {\mathbbm{d}_{\mathsf{N}}}+1}$ of $\widehat{\Lambda}_{1,P}$ converge in probability to $\{ \tau_j \}_{j\geq 1}$. That is, the first ${\mathbbm{d}_{\mathsf{N}}}$ eigenvalues of $\widehat{\Lambda}_{1,P}$ are associated (at least asymptotically) with the nonstationary subspace, and therefore it is the eigenvalues $\widehat\tau_j$ for $j={\mathbbm{d}_{\mathsf{N}}}+1,\ldots$ that correspond to to eigenvalues $\tau_j$ for $j=1,\ldots$ of $\Lambda_{\mathsf{S},R}$.

We may thus avoid violation of the condition $\tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}} \neq \tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}+1}$ by choosing $\mathrm{K}$ such that $\widehat{\tau}_{\mathrm{K}}$ is sufficiently greater than $\widehat{\tau}_{\mathrm{K}+1}$. This approach is pragmatic and also theoretically supported under the following assumption (which is required only when $a_R>0$).

The conditions in Assumption (ref) are stronger than those in the previous sections. In particular, $\{P_{\mathsf{S}}X_t\}_{t\geq 1}$ needs to satisfy a stronger moment condition than required for $\{X_t\}_{t\geq 1}$ in Section (ref), and the coefficients $\{\Phi_j\}_{j\geq 1}$ are required to decay at a faster rate than in Section (ref). The latter can be seen from the fact that $\|\sum_{j=m}^\infty \tilde{\Phi}_j \|_{\operatorname{op}} = \scalebox{0.7}{$\mathcal{O}$} (m^{-\gamma})$ for some $\gamma>1$ implies the summability condition in \hyperref[{eqlinear2}]{\tagform@{\ref*{eqlinear2}}}. Assumption (ref)(iii) imposes a regularity condition on $\mathrm{k}_P(\cdot)$ and a summability condition on the autovariance operators of $\{P_{\mathsf{S}}X_t\}_{t\geq 1}$, and Assumption (ref)(iv) requires that $h_P$ is not too small relative to $T$. Nonetheless, even if Assumption (ref) includes some stronger conditions, we note that they are not essential for our proposed ADI tests but are only needed if practitioners want to check if $\tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}} \neq \tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}+1}$ based on the sample eigenvalues $\widehat\tau_j$. The following theorem establishes the desired result.

From Theorem (ref) we deduce that, for any $\eta_T$ satisfying $\eta_T M_T^{-1/2} \to 0$ and $\eta_T \to \infty$,

This result can be applied in practice as an informal way to check whether the condition $\tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}} \neq \tau_{\mathrm{K}-{\mathbbm{d}_{\mathsf{N}}}+1}$ in Assumption (ref) is likely to hold. Since we do not need to impose Assumption (ref) in Theorem (ref) when $\mathrm{K} \leq {\mathbbm{d}_{\mathsf{N}}}$, we do not need to consider $\mathrm{K} \leq {\mathbbm{d}_{\mathsf{N}}}$ in Theorem (ref).

Inclusion of deterministic components

So far we have assumed that $\{X_t\}_{t \geq 1}$ has no deterministic component. However, in practice it is quite common that observed time series contain a deterministic component. In this section, we focus on the case with a nonzero intercept and a linear trend since those seem to be most relevant in practice.

Suppose that the observed time series $\{Y_t\}_{t\geq 1}$ is

where $\{X_t\}_{t\geq 1}$ satisfies Assumption (ref). We let $U_t^{(1)}$ and $U_t^{(2)}$ be defined by

where $ \bar{Y} = T^{-1}\sum_{t=1}^T Y_t$ and $\bar{t} = T^{-1}\sum_{t=1}^T t = (T+1)/2$. Thus, $U_t^{(1)}$ and $U_t^{(2)}$ are, respectively, the mean-adjusted and the trend-adjusted residuals of $Y_t$ kokoszka2016kpss.

Our ADI variance ratio tests can be adjusted to accommodate deterministic terms by replacing the sample operators computed from $\{X_t\}_{t\geq1}$ with the corresponding operators computed from $\{U_t^{(1)}\}_{t\geq1}$ or $\{U_t^{(2)}\}_{t\geq1}$. When the model for $Y_t$ is given by \hyperref[{equnobs}]{\tagform@{\ref*{equnobs}}} with $\zeta_2 =0$ then we use $U_t^{(1)}$, and otherwise we use $U_t^{(2)}$.

This adjustment to accommodate deterministic terms produces some obvious changes in the asymptotic results given in Theorems (ref) and (ref) as in the existing literature (e.g., NSS, Section 3.4; seo2024functional, Section S3). For example, the extension for the VR$(d_L,1)$ test can easily be done as in Theorem 3 of NSS. More specifically, in Theorems (ref) and (ref), $W_{1,{\mathbbm{d}_{\mathsf{N}}}}$ needs to be replaced by ${\mathbbm{d}_{\mathsf{N}}}$-dimensional demeaned standard Brownian motion (if $\zeta_2=0$) or detrended standard Brownian motion (if $\zeta_2\neq 0$), and the definition of $W_{d_L,{\mathbbm{d}_{\mathsf{N}}}}$ needs to change accordingly. Similar changes have to be made in Theorem (ref), but where, as may be deduced from our proof of Theorem (ref), $B_{\mathrm{K} - {\mathbbm{d}_{\mathsf{N}}}}$ needs to be replaced with the standard Brownian bridge (if $\zeta_2=0$) or the second-level standard Brownian bridge (if $\zeta_2\neq 0$), while $W_{d_L,{\mathbbm{d}_{\mathsf{N}}}}$ is the $(d_L-1)$-fold integrated demeaned (if $\zeta_2=0$) or detrended (if $\zeta_2\neq 0$) Brownian motion. The asymptotic results given in Theorems (ref) and (ref) do not require any changes once $\widehat{\Lambda}_{1,P}$ is constructed from the relevant residuals (our proofs of Theorems (ref) and (ref) contain related discussions). Following these changes to the theorems, it is quite obvious to make relevant changes to the corollaries given in the previous sections, and hence the details are omitted.

Determination of dimension of nonstationary subspace

Practitioners may mostly be interested in determining the dimension of $\mathcal H_{\mathsf{N}}$ rather than testing a specific hypothesis on the dimension. In this section we first propose sequential testing procedures for doing so, and then we consider an eigenvalue ratio estimation approach.

Estimation via sequential testing

First, we propose a “top-down” (TD) sequential procedure to estimate ${\mathbbm{d}_{\mathsf{N}}}$. Here we test \hyperref[{eqhypo}]{\tagform@{\ref*{eqhypo}}} sequentially (using any of the tests given in Corollaries (ref) and (ref)) for $\mathbbm{d}_0 = \mathbbm{d}_{\max}, \mathbbm{d}_{\max}-1, \ldots$, where $\mathbbm{d}_{\max}$ is an upper bound on the dimension of $\mathcal H_{\mathsf{N}}$ that can reasonably be chosen to be (slightly) greater than ${\mathbbm{d}_{\mathsf{N}}}$ in practice; e.g., Remark (ref). We let $\widehat\mathbbm{d}_{\rm TD}$ denote the value of $\mathbbm{d}_0$ under the first non-rejected null hypothesis; if the null hypothesis is rejected for every $\mathbbm{d}_0 = \mathbbm{d}_{\max}, \mathbbm{d}_{\max}-1, \ldots , 1$, then $\widehat\mathbbm{d}_{\rm TD} =0$.

Second, we consider a “bottom-up” (BU) procedure, where we test \hyperref[{eqhypo2}]{\tagform@{\ref*{eqhypo2}}} sequentially for $\mathbbm{d}_0 = 0 , 1 , \ldots$ using any of the inverse VR tests in Corollary (ref). We then let $\widehat\mathbbm{d}_{\rm BU}$ denote the first non-rejected null hypothesis. This procedure is attractive because no $\mathbbm{d}_{\max}$ needs to be specified, and the condition $\mathbbm{d}_{\max} \geq {\mathbbm{d}_{\mathsf{N}}}$ is not required; see also Remark (ref).

The sequential procedures can be applied with a different $\mathrm{K}$ and $\widehat{P}_{\mathrm{K}}$ for each value of $\mathbbm{d}_0$, i.e.\ in each step of the algorithm. Then $\mathrm{K}$ should be understood as a function of $\mathbbm{d}_0$, denoted $\mathrm{K} = \mathrm{K}(\mathbbm{d}_0)$, although we mostly suppress the dependence on $\mathbbm{d}_0$ to simplify the notation. In practice, $\mathrm{K}$ can simply be set to $\mathrm{K} = \mathbbm{d}_{\max} + m$ or $\mathrm{K} = \mathrm{K}(\mathbbm{d}_0) = \mathbbm{d}_0 + m$ for some integer $m \geq 1$ as in NSS; see Remark (ref) for details.

The consistency results for both procedures are given in the following theorem.

The proof of Theorem (ref) follows directly from earlier results and is therefore omitted. For both procedures we note that if $\alpha \to 0$ as $T \to \infty$ then the procedures are consistent; that is $\mathrm{P}(\widehat\mathbbm{d}_{\rm TD} = {\mathbbm{d}_{\mathsf{N}}})~\to~1$ and $\mathrm{P}(\widehat\mathbbm{d}_{\rm BU} = {\mathbbm{d}_{\mathsf{N}}})~\to~1$.

As discussed in Remarks (ref) and (ref), the conditions required for the TD procedure are stronger than those required for the BU procedure. The BU procedure has the additional advantage that it can be implemented without specifying an upper bound $\mathbbm{d}_{\max}$ on ${\mathbbm{d}_{\mathsf{N}}}$ and, specifically, the choices of $\mathrm{K}$ for the BU procedure to consistently estimate ${\mathbbm{d}_{\mathsf{N}}}$ (with shrinking significance levels as discussed above) depend only on $\mathbbm{d}_0$ (not ${\mathbbm{d}_{\mathsf{N}}}$). However, Monte Carlo simulations have shown that the TD procedure has better finite-sample properties. Thus, we combine the above two sequential procedures into an “up-down” (UD) hybrid procedure that enjoys the advantages of both the TD and BU procedures. Specifically, we first apply the BU procedure to obtain $\widehat\mathbbm{d}_{\rm BU}$ and use this to specify a data-driven value of the upper bound as $\mathbbm{d}_{\max} = \widehat\mathbbm{d}_{\rm BU} +m_\mathbbm{d}$ for an integer $m_\mathbbm{d} \geq 0$. We then apply the TD procedure initiated at this upper bound to obtain the UD hybrid estimate, $\widehat\mathbbm{d}_{\rm UD}$. Note that the upper bound obtained in the first (BU) part of the UD procedure can also be used to choose the same $\mathrm{K} = \mathbbm{d}_{\max} + m$ for an integer $m \geq 0$ in all steps in the second (TD) part of the UD procedure. The finite-sample properties of this UD hybrid procedure, along with those of the other procedures, will be investigated by simulations in Section (ref).

Eigenvalue ratio estimator

Recent work concerning vector-valued time series provide ways to consistently estimate the number of stochastic trends using the fact that the eigenvalues of the sample autocovariance matrices have different stochastic orders if the time series contains components of different integration orders Zhang2018,Zhang_et_al2019,Chang2021. In a functional time series setting, the eigenvalue ratio estimator proposed by LRS2020,LRS2020nonst can be used to estimate the dimension of the component of the highest integration order. A slight modification of their estimator was discussed in Remark (ref).

Alternatively, based on the asymptotic results in Section (ref), we may use the generalized eigenvalue problem \hyperref[{eqgev}]{\tagform@{\ref*{eqgev}}} to construct an eigenvalue ratio estimator for ${\mathbbm{d}_{\mathsf{N}}}$. Let $\mu_j$ be an eigenvalue of \hyperref[{eqgev}]{\tagform@{\ref*{eqgev}}}. As we prove in Theorem (ref) below, we can use Theorems (ref)--(ref) and arguments from their proofs to show that $\mu_j^{-1}\mu_{j+1} = \mathcal{O}_{\mathrm{P}} (1)$ if $j < {\mathbbm{d}_{\mathsf{N}}}$ (resp.\ $j > {\mathbbm{d}_{\mathsf{N}}}$) since both $n_T\mu_j$ and $n_T\mu_{j+1}$ converge (resp.\ diverge at the same rate), while $\mu_j^{-1}\mu_{j+1} \rightarrow_{\mathrm{P}} \infty$ if $j={\mathbbm{d}_{\mathsf{N}}}$. Hence, the ratio $\mu_j^{-1}\mu_{j+1}$ will be maximized at $j={\mathbbm{d}_{\mathsf{N}}}$. This motivates the estimator

where $\mathbbm{d}_{\max}$ is a pre-specified integer indicating an upper bound on ${\mathbbm{d}_{\mathsf{N}}}$. We note that these eigenvalue ratio-type estimators are only applicable when ${\mathbbm{d}_{\mathsf{N}}} \geq 1$; hence this needs to be examined first in practical implementation of these estimators. A convenient way to do this is to test \hyperref[{eqhypo2}]{\tagform@{\ref*{eqhypo2}}} with $\mathbbm{d}_0=0$ using the tests in Corollary (ref).

The estimator in \hyperref[{eqeigratio}]{\tagform@{\ref*{eqeigratio}}} is very similar to the eigenvalue ratio estimator proposed by LRS2020,LRS2020nonst, see also Remark (ref), but there are some important differences. First, their estimator is constructed directly from the eigenvalues of a sample variance operator, while we use the generalized eigenvalue problem \hyperref[{eqgev}]{\tagform@{\ref*{eqgev}}}. Second, their estimator requires a tuning parameter, say $\kappa >0$, to deal with estimation error associated with small eigenvalues; specifically, they set $\widehat\mathbbm{d}_{\text{\tiny{LRS}}} = \arg\min_{1\leq j\leq \mathbbm{d}_{\max}} \{ \hat\kappa_{j+1}/\hat\kappa_j \}$, where $\hat\kappa_j$ is the $j$-th largest eigenvalue of the sample variance operator and $\hat\kappa_j$ is regarded as zero if $\hat\kappa_j / \hat\kappa_1$ is smaller than $\kappa$, and they define $0/0=1$. However, our estimator in \hyperref[{eqeigratio}]{\tagform@{\ref*{eqeigratio}}} does not require such a tuning parameter.

Inference on subspaces

In practical applications it may be of interest to test hypotheses about $\mathcal H_{\mathsf{N}}$ or $\mathcal H_{\mathsf{S}}$. For a specified subspace $\mathcal H_0$, consider the following hypotheses:

For example, we may be interested in testing if a specific element $v \in \mathcal H$ is included in $\mathcal H_{\mathsf{N}}$, and hence can be interpreted as one of the common stochastic trends. In that case we can use \hyperref[{hptest1}]{\tagform@{\ref*{hptest1}}} with $\mathcal H_0 = \operatorname{span}\{v\}$. We can also test if a specified subspace, e.g.\ $\mathcal H_0 = \operatorname{span}\{v_1,v_2\}$, contains the entire nonstationary subspace by using \hyperref[{hptest2}]{\tagform@{\ref*{hptest2}}}. Finally, it may be of interest to examine if a specified element or subspace is in $\mathcal H_{\mathsf{S}}$, and then we can use \hyperref[{hptest3}]{\tagform@{\ref*{hptest3}}}. In a finite-dimensional setting, \hyperref[{hptest3}]{\tagform@{\ref*{hptest3}}} corresponds to testing if a vector or matrix is cointegrating.

Let $P_{\mathcal H_0}$ denote the projection onto $\mathcal H_0$, let $p_0 = \dim(\mathcal H_{0})$, and let ${\mathbbm{d}_{\mathsf{N}}}$ be known (or replaced by a consistent estimator from Section (ref)). \phantomsectionThe hypotheses \hyperref[{hptest1}]{\tagform@{\ref*{hptest1}}} and \hyperref[{hptest2}]{\tagform@{\ref*{hptest2}}} can then be tested by investigating the dimension of the nonstationary subspace associated with the residual series $\{ (I-P_{\mathcal H_0}) X_t\}_{t\geq 1}$, and the hypothesis \hyperref[{hptest3}]{\tagform@{\ref*{hptest3}}} can be tested by investigating the dimension of the nonstationary subspace associated with the projected series $\{ P_{\mathcal H_0} X_t\}_{t\geq 1}$. Specifically, if $H_0$ of \hyperref[{hptest1}]{\tagform@{\ref*{hptest1}}} is true, then $\{(I-P_{\mathcal H_0}) X_t\}_{t\geq 1}$ contains ${\mathbbm{d}_{\mathsf{N}}}-p_0$ stochastic trends, if $H_0$ of \hyperref[{hptest2}]{\tagform@{\ref*{hptest2}}} is true then $\{ (I-P_{\mathcal H_0}) X_t\}_{t\geq 1}$ contains zero stochastic trends, and if $H_0$ of \hyperref[{hptest3}]{\tagform@{\ref*{hptest3}}} is true then $\{ P_{\mathcal H_0} X_t\}_{t\geq 1}$ contains zero stochastic trends. Under $H_1$ in \hyperref[{hptest1}]{\textup{\tagform@{\ref*{hptest1}}}}--\hyperref[{hptest3}]{\textup{\tagform@{\ref*{hptest3}}}}, there are more stochastic trends in the relevant time series than under the null. Thus, by analyzing $\{ (I-P_{\mathcal H_0}) X_t\}_{t\geq 1}$ or $\{P_{\mathcal H_0} X_t\}_{t\geq 1}$, \hyperref[{hptest1}]{\textup{\tagform@{\ref*{hptest1}}}}--\hyperref[{hptest3}]{\textup{\tagform@{\ref*{hptest3}}}} may be readily tested by our inverse VR test, assuming the projected time series satisfy the conditions for the inverse VR test.

Of course, the high-level conditions in Assumption (ref), which are applied in Theorem (ref), may be replaced by appropriate low-level conditions as discussed in Theorems (ref) and (ref).

The approach to testing the hypotheses \hyperref[{hptest1}]{\tagform@{\ref*{hptest1}}} and \hyperref[{hptest2}]{\tagform@{\ref*{hptest2}}} based on $\{(I-P_{\mathcal H_0})X_t\}_{t\geq 1}$ is parallel to that in seo2024functional who relied on an FPCA-based test for examining the number of stochastic trends embedded in a cointegrated functional time series. However, seo2024functional requires the condition that the long-run variance of $\{P_{\mathsf{S}} X_t\}_{t\geq 1}$ is positive definite on $\mathcal H_{\mathsf{S}}$, and this may be restrictive in cases where $X_t$ only takes values on an unknown (and possibly finite-dimensional) subspace. This condition is not required under our setup.

It is not possible, in general, to investigate the hypothesis \hyperref[{hptest1}]{\tagform@{\ref*{hptest1}}} (resp.\ \hyperref[{hptest2}]{\tagform@{\ref*{hptest2}}}) by testing if there are $p_0$ (resp.\ ${\mathbbm{d}_{\mathsf{N}}}$) stochastic trends in the projected time series $\{ P_{\mathcal H_0} X_t\}_{t\geq 1}$. For example, let $\{f_j\}_{j\geq 1}$ be an orthonormal basis of $\mathcal H$ and suppose $\mathcal H_{\mathsf{N}} = \operatorname{span}\{f_1\}$, i.e., ${\mathbbm{d}_{\mathsf{N}}}=1$. Consider $\mathcal H_0= \operatorname{span}\{g\}$ for $g = 0.5f_1+0.5f_2$, which can be a null hypothesis in either \hyperref[{hptest1}]{\tagform@{\ref*{hptest1}}} or \hyperref[{hptest2}]{\tagform@{\ref*{hptest2}}}. Both null hypotheses are false since neither $\operatorname{span}\{g\} \subseteq \operatorname{span}\{f_1\}$ nor $\operatorname{span}\{g\} \supseteq \operatorname{span}\{f_1\}$ is true. However, the projected time series $P_{\mathcal H_0} X_t=\langle X_t,g \rangle = \langle X_t,f_1 \rangle + \langle X_t,f_2 \rangle$ is nonstationary and has one stochastic trend. \phantomsectionSimilarly, $H_0$ in \hyperref[{hptest3}]{\tagform@{\ref*{hptest3}}} cannot be investigated by testing if there are ${\mathbbm{d}_{\mathsf{N}}}$ stochastic trends in $\{ (I-P_{\mathcal H_0}) X_t\}_{t\geq 1}$. These impossibilities are discussed in Theorem 4.2 in FranchiParuolo2025, which characterizes the structure of the parameter space under the null and alternative of hypotheses of type \hyperref[{hptest1}]{\tagform@{\ref*{hptest1}}} and \hyperref[{hptest2}]{\textup{\tagform@{\ref*{hptest2}}}}.

Monte Carlo study

We study the finite-sample performance of the proposed ADI tests using a Monte Carlo simulation setup similar to the functional AR(1) model in LRS2020nonst, NSS, and seo2024functional. Let $\{ g_j \}_{j \geq 1}$ be the Fourier basis functions on $[0,1]$. For each ${\mathbbm{d}_{\mathsf{N}}}$, we let $P_{\mathsf{N}} \Delta X_t$ and $P_{\mathsf{S}} X_t$ be generated by the following stationary functional AR(1) models,

where $\{ g_j^N \}_{j=1}^{{\mathbbm{d}_{\mathsf{N}}}}$ are randomly drawn from $\{ g_1 , \ldots , g_{{\mathbbm{d}_{\mathsf{N}}}+3} \}$ (${\mathbbm{d}_{\mathsf{N}}} \leq 8$ in our simulation experiments), $\{ g_j^S \}_{j=1}^{12}$ are randomly drawn from $\{ g_{15}, \ldots , g_{30} \}$, and $\varepsilon_t$ is given by $\sum_{j=1}^{40} (0.9)^{j-1} \theta_{j,t} g_j$ for standard normal random variables $\{ \theta_{j,t} \}_{j \geq 1}$ that are independent across $j$ and $t$. We let $\alpha_j$ and $\beta_j$ be independent uniform random variables on $[-0.8,0.8]$ for each realization of the data generating process (DGP). In practice, it is common to have a nonzero intercept or a linear time trend as in Section (ref). We here consider the former case and add an intercept $\zeta_1$, which is also randomly chosen, to each realization of the DGP. Specifically, $\zeta_1 = \sum_{j=1}^{30} (0.9)^{j-1} \tilde{\theta}_j g_j$, where $\{ \tilde{\theta}_j \}_{j=1}^{30}$ are independent standard normal random variables. Viewed as a time series taking values in the $L^2[0,1]$ Hilbert space, we may obtain the eigenelements of the sample (long-run) variance of $X_t$ using FPCA (see e.g., Chang2016; NSS; LRS2020nonst). To this end, we represent $X_t$, which is assumed to be observed on 200 regularly spaced points of $[0,1]$, using the first 40 Legendre basis functions.

To choose bandwidths, we follow Remarks (ref), (ref), (ref), and (ref). For VR tests with $d_R = 1$ we let $h_P=h_L=h_R=0$ (Remarks (ref) and (ref)). For the remaining tests, we let $h_P=[T^{1/3}]$ because a large bandwidth is supported by Assumption (ref). For VR tests with $d_R = 0$ we let $h_R = h_{\rm opt}$, which is the andrews1991 optimal bandwidth (Remark (ref)), but for the inverse tests the series are nonstationary, so we let $h_R=[T^{1/5}]$, which is the same rate as $h_{\rm opt}$ (Remark (ref)). We consider both $h_L=0$ and $h_L=[ T^{1/3}]$ (Remark (ref)). Throughout, we let $\mathrm{k}_\mathit{m}(\cdot)$ be the Tukey-Hanning kernel for $m \in \{ L,R,P \}$, all test statistics apply the trace test functional $\mathcal{F}_{\operatorname{trace}}$, and the slack extraction uses $\mathrm{K} = \mathbbm{d}_0 +2$ following Remark (ref).

\phantomsectionWe present simulation results for the VR(2,1), VR(2,0), VR(1,0), and inverse VR tests (the latter based on the (2,1) problem). The VR choices (2,1), (1,0), and the inverse VR are motivated by well-known tests from finite-dimensional Euclidean space, see Remark (ref) and Appendix (ref), and the VR(2,0) is in-between. We also performed simulations for VR tests based on larger values of $d_L$, but found that their performances were inferior to the tests reported here (results available upon request).

Table (ref) reports rejection frequencies. We consider true dimension ${\mathbbm{d}_{\mathsf{N}}} \in \{ 0,1,3,5,7 \}$. For the VR tests we consider the null $\mathbbm{d}_0 = {\mathbbm{d}_{\mathsf{N}}}$ and alternatives $\mathbbm{d}_0 = {\mathbbm{d}_{\mathsf{N}}} +1$ and $\mathbbm{d}_0 = {\mathbbm{d}_{\mathsf{N}}} +2$, and for the inverse VR test we consider the null $\mathbbm{d}_0 = {\mathbbm{d}_{\mathsf{N}}}$ and alternatives $\mathbbm{d}_0 = {\mathbbm{d}_{\mathsf{N}}} -1$ and $\mathbbm{d}_0 = {\mathbbm{d}_{\mathsf{N}}} -2$. The results show that the VR(1,0) and VR(2,0) tests are seriously under-sized when $a_L=0$, as expected based on shintani2001simple, especially for the larger values of ${\mathbbm{d}_{\mathsf{N}}}$ considered, and both tests suffer from lower power in those cases as a consequence. Size control is much better for $h_L=[T^{1/3}]$ for the VR(1,0) test, but not for the VR(2,0) test. The VR(2,1) test has very good size properties in all cases considered except the largest value of ${\mathbbm{d}_{\mathsf{N}}}$ with the smallest sample size, and it has excellent power. Finally, we see that the inverse VR test is slightly over-sized in some cases. Its power is nearly as impressive as that of the VR(2,1) test, but of course the alternatives are different for the inverse VR test, and it does not have the natural upper bound $\mathbbm{d}_{\max}$ that the VR tests have.

In Table (ref) we investigate the finite-sample performance of the tests when they are applied to estimate ${\mathbbm{d}_{\mathsf{N}}}$ using one of the sequential procedures described in Section (ref). For comparison, we also consider the eigenvalue ratio estimator $\ddot\mathbbm{d}$ computed from the VR(2,1) problem (see Theorem (ref)), and \citepos{LRS2020nonst} estimator, which is denoted $\widehat\mathbbm{d}_{\text{\tiny{LRS}}}$. For each method and each value of ${\mathbbm{d}_{\mathsf{N}}}$, the table reports the frequencies of correctly estimating the true dimension ${\mathbbm{d}_{\mathsf{N}}}$. Overall, the direct estimators $\ddot\mathbbm{d}$ and $\widehat\mathbbm{d}_{\text{\tiny{LRS}}}$ perform worse than the VR(2,1) and inverse VR sequential testing procedures, and between the two eigenvalue ratio estimators, $\ddot\mathbbm{d}$ tends to perform better than $\widehat\mathbbm{d}_{\text{\tiny{LRS}}}$. The sequential TD procedure based on the VR(2,1) tests performs very well, and the BU procedure (based on the inverse VR test) performs nearly as well. On the other hand, the TD procedures based on the VR(2,0) and VR(1,0) tests perform very poorly for ${\mathbbm{d}_{\mathsf{N}}} \geq 3$, which was expected based on their poor power observed in Table (ref). The choice of $h_L$ seems to significantly affect the performance of the VR(1,0) tests, especially when ${\mathbbm{d}_{\mathsf{N}}}$ is large. Of course, this is not the case for the TD procedure based on the VR(2,1) test, and indeed, this method appears to perform the best overall.

Results for the UD hybrid procedure are also reported in Table (ref). This enjoys the advantage of a data-driven selection of $\mathbbm{d}_{\max} = \widehat\mathbbm{d}_{\rm BU}+5$. This $\mathbbm{d}_{\max}$ is random, but the good performance of the BU procedure, the addition of 5 from $\widehat\mathbbm{d}_{\rm BU}$, and the strong robustness of the TD procedures to the choice of $\mathbbm{d}_{\max}$, imply that in fact the UD procedures enjoy nearly identical frequencies of correct estimation of ${\mathbbm{d}_{\mathsf{N}}}$ as the corresponding TD procedures. Thus, the advantage of the data-driven selection of $\mathbbm{d}_{\max}$ comes without cost, at least for these sample sizes, and consequently the UD procedure based on VR(2,1) is the preferred method in our Monte Carlo simulations.

Finally, recall that the functional observations are constructed from 200 regularly spaced points on $[0,1]$ by smoothing using the first 40 Legendre basis functions. However, as discussed earlier, our methodology can also be applied to the time series of 200-dimensional vectors of discrete realizations of $X_t$ with no additional theoretical modification. We thus repeated the analysis in Table (ref) letting $\{X_t\}_{t\geq 1}$ be viewed as a 200-dimensional vector-valued time series, but the results were nearly identical and are hence omitted.

Empirical applications

In this section we illustrate our methdology with two real-world data examples.

Corporate bond yield curves

We apply the proposed tests to examine the number of common stochastic trends in the monthly time series of high quality market corporate bond yield curves (a.k.a.\ the HQM yield curve). The data is from the US Department of Treasury available at \url{https://home.treasury.gov/} and spans Jan.\ 1984 to Dec.\ 2018. At each time $t$, we observe a corporate bond with 200 different maturities ranging from 6 months to 100 years. This data was studied by BT2022, who assumed that this high-dimensional time series is driven by a few factors of deterministic trends and I(1) and I(0) processes, and estimated the number of those factors using their own testing procedure. However, it may be more natural to view this time series as realizations of yield curves that are either finite-dimensional (of high and unknown dimension) or curve-valued. With this point of view, we apply our methodology to determine the number of stochastic trends. To implement our methodology, we represent the HQM yield data by the first 40 Legendre basis functions as in Section (ref).

Figure (ref) shows the time series of HQM yield curves (observed at 200 different maturities) and Figure (ref) shows the time series of yields at fixed maturities of 3, 6, 10, and 30 years. The figures suggest that the HQM yield data includes a linear time trend. The existence of a linear trend is tested and confirmed by BT2022 for the data span of Jan.\ 1985 to Sep.\ 2017. Thus, we apply the trend-adjusted testing procedures in Section (ref).

The empirical results are reported in Table (ref). We implement all procedures as in Section (ref) using $h_L=[T^{1/3}]$ for the VR(2,0), VR(1,0), and inverse VR tests, and using $\mathrm{K} = \mathbbm{d}_0 +2$ for slack extraction. The eigenvalue ratio estimates in the first part of the table are $\ddot\mathbbm{d} = 6$ and $\widehat\mathbbm{d}_{\text{\tiny{LRS}}} = 5$ (using $\kappa=10^{-4}$), although the ratios are close to suggesting $\ddot\mathbbm{d} = 2$ and $\widehat\mathbbm{d}_{\text{\tiny{LRS}}} = 3$.

In the second part of Table (ref) we report the VR and inverse VR test statistics for a range of $\mathbbm{d}_0$. The inverse VR tests imply that $\widehat\mathbbm{d}_{\rm BU} = 3$ (at 1% level) or $\widehat\mathbbm{d}_{\rm BU} = 6$ (at 5% level). Given this result, the choice $\mathbbm{d}_{\max} = 7$ seems like a reasonable input to the UD procedures. The latter give either $\widehat\mathbbm{d}_{\rm UD} = 2$ or $\widehat\mathbbm{d}_{\rm UD} = 3$ depending on the choice of ($d_L,d_R)$ and level. However, the BU and TD procedures (with (2,1) and (2,0)) only agree at 1% level, in which case all three find $\widehat\mathbbm{d} = 3$. Overall, the evidence thus supports ${\mathbbm{d}_{\mathsf{N}}} = 3$ stochastic trends.

The Nelson-Siegel term structure model

Given an estimated dimension, it is possible to examine various hypotheses on the nonstationary subspace as in Section (ref). A popular model for the term structure of interest rates was proposed by NelsonSiegel87; henceforth the N-S model. They represent the term structure of yield curves by a linear combination of three parametric functions,

where $\tau$ is time to maturity, $\varsigma$ is a shape parameter, and $\gamma_0$, $\gamma_1$, and $\gamma_2$ denote level, slope, and curvature factors. Because $\gamma_1$ and $\gamma_2$ are nearly collinear for large $\tau$, we subsequently consider only data for $\tau \in [0,30]$ years. This has almost no impact on the results reported in Table (ref). We now investigate to what extent the N-S model can explain the dominant (nonstationary) variation in the data. To this end, we need to specify the shape parameter, $\varsigma$. Since our data is annualized, we set $\varsigma = 1.37$ following DieboldLi2006.

We first test the hypothesis $H_0$ in \hyperref[{hptest3}]{\tagform@{\ref*{hptest3}}} with $\mathcal H_0 = \operatorname{span} \{ \gamma_0 , \gamma_1, \gamma_2 \}$ using Theorem (ref)(iii). That is, we test whether “fitted values,” $\{ P_{\mathcal H_0}X_t \}_{t\geq 1}$, from the N-S model are stationary. Not surprisingly, this hypothesis is strongly rejected. The same conclusion holds for tests of \hyperref[{hptest3}]{\tagform@{\ref*{hptest3}}} when $\mathcal H_0$ is specified as the span of any one or two of the factors $\gamma_i$, $i=0,1,2$. This confirms that none of the factors are stationary.

Next, we test if the level, slope, and curvature factors span the nonstationary subspace as in \hyperref[{hptest2}]{\tagform@{\ref*{hptest2}}}. From our earlier results in Table (ref), we set the dimension of the latter to ${\mathbbm{d}_{\mathsf{N}}} = 3$ (testing at 1% level). Thus, there is no possibility that the N-S model with only one or two factors can yield stationary residuals. With three factors, it is of interest to test if the N-S model can explain all the nonstationary variation in the yield curve data and, as a result, can yield stationary residuals. In this sense, the hypothesis $H_0 : \mathcal H_0 = \operatorname{span} \{ \gamma_0 , \gamma_1, \gamma_2 \} \supseteq \mathcal H_{\mathsf{N}}$ in \hyperref[{hptest2}]{\tagform@{\ref*{hptest2}}} is very interesting, and its non-rejection would be the strongest conclusion that we could obtain. However, implementing the test as in Theorem (ref)(ii) with $\mathbbm{d}_0 = 0$ and $p_0 = 3$, the inverse VR statistic is $1.47$ (see Table (ref)) and rejects the null at 1% level.

We proceed to test $H_0$ in \hyperref[{hptest1}]{\tagform@{\ref*{hptest1}}}, where $\mathcal H_0$ is specified as various choices of factor(s). That is, we test whether any one factor can reduce the dimension of the nonstationary subspace by one, and also whether any two factors can reduce the dimension by two. Of course, the latter is less likely to be accepted than the former. For one factor, implementing the test as in Theorem (ref)(i), the inverse VR statistics with $\mathbbm{d}_0 = 2$ are $0.24$, $0.21$, and $0.22$ for $\mathcal H_0 = \operatorname{span} \{ \gamma_0 \}$, $\mathcal H_0 = \operatorname{span} \{ \gamma_1 \}$, $\mathcal H_0 = \operatorname{span} \{ \gamma_2 \}$, respectively, in each case rejecting the null at 1% level. For two factors, the inverse VR statistics with $\mathbbm{d}_0 = 1$ are $0.40$ and $0.38$ for $\mathcal H_0 = \operatorname{span} \{ \gamma_0 , \gamma_1 \}$ and $\mathcal H_0 = \operatorname{span} \{ \gamma_1 , \gamma_2 \}$, respectively, in both cases rejecting the null at 1% level. Thus, there is no evidence to suggest that any one factor is in the nonstationary subspace, i.e., that any one factor can be interpreted as one of the common stochastic trends.

Finally, following Remark (ref), we investigate if the nonstationary subspace of the residual series, $\{ (I-P_{\mathcal H_0})X_t \}_{t\geq 1}$ with $\mathcal H_0 = \operatorname{span} \{ \gamma_0 , \gamma_1, \gamma_2 \}$, is of smaller dimension than that of $X_t$ itself. Letting the former be denoted $\mathbbm{d}_{\mathcal H_0}$, we apply the sequential procedures in \hyperref[{sequential1}]{\tagform@{\ref*{sequential1}}} and \hyperref[{sequential2}]{\tagform@{\ref*{sequential2}}}. The results are reported in Table (ref). The BU procedure yields $\widehat\mathbbm{d}_{\mathcal H_0 , {\rm BU}} = 3$ at 1% level and rejects $\mathbbm{d}_0 = 3$ at 5% level. Applying the TD sequential procedure starting from $\mathbbm{d}_{\max} = 3$, we find $\widehat\mathbbm{d}_{\mathcal H_0 , {\rm TD}} = 2$ with the preferred VR(2,1) test at 5% level (but $\widehat\mathbbm{d}_{\mathcal H_0 , {\rm TD}} = 1$ at 10% level with the VR(2,1) test, $\widehat\mathbbm{d}_{\mathcal H_0 , {\rm TD}} = 1$ at 5% level with the VR(2,0) test, and $\widehat\mathbbm{d}_{\mathcal H_0 , {\rm TD}} = 0$ at 1% level with the VR(1,0) test). Thus, we do find some evidence that three factors in the N-S model explain some of the nonstationary variation in the data.

The empirical results for the N-S model are specific to the value of $\varsigma$, at least for the slope and curvature factors. As a robustness check, we repeated the analysis for $\varsigma = 1,2,3,4,5$, and the results were qualitatively the same. In conclusion, the three parametric factors of the N-S model cannot fully explain the nonstationarity in the corporate yield curve data, so the residuals from the model will still exhibit nonstationarity. However, the model can capture some of the nonstationarity, in the sense that the dimension of the nonstationary subspace of the residuals is smaller than that of the original yield curve data.

Labor market indices

In empirical (macro) analysis, high-dimensional observations are often given by a collection of potentially nonstationary variables that are closely related to each other but do not have a natural ordering. In such cases, not only is each observation naturally understood as a realization of a high-dimensional random vector that cannot generally be smoothed to a curve (unlike in our previous example of yield curve data), but also the number of linearly independent I(1) stochastic trends may be substantially smaller than the total number of variables. As discussed, in this case as well, our proposed methodology can be applied to study the number of stochastic trends and/or stationarity without any modifications. To illustrate this versatility of our methodology, we now consider another empirical example analyzing 29 monthly labor market indices provided in the FRED-MD data set McCracken2016. The indices include the civilian labor force, employment/unemployment rate, hourly wage, and other labor-related economic variables; a detailed list can be found in the table labeled “Group 2” in the appendix of McCracken2016. We apply log-transformation to the data.\footnote{The log transformation was suggested by McCracken2016 for some variables. We found that the remaining variables are often very large in scale (e.g., “HWI”) or have only positive values. Thus, those variables are also log-transformed in our testing procedure.} The time span used in this analysis is Jan.\ 1990 to Dec.\ 2019.

In Figure (ref) we display time series of $\langle X_t,\widehat{\nu}_j \rangle$ for some eigenvectors $\widehat{\nu}_j$ that are computed from $\widehat\Lambda_{0,R}$. These are the $j^{\rm th}$ score processes (see Chang2016 and Section 5 of NSS). It is clear from Figure (ref) that there is some nonstationary behavior (e.g., $j=1$), but this seems to eventually disappear leaving stationary processes (e.g., $j=10$).

McCracken2016 suggested that some of these labor market indices may be I(2) variables.\footnote{Specifically, FRED ID variables CES0600000008, CES2000000008, and CES3000000008.} To examine the existence of this type of higher-order stochastic trends, we first display in Figure (ref) the first three score processes, $\{\langle \Delta X_t,\widehat{\nu}_j \rangle\}_{t=1}^T$. Some of these must behave as a unit root process if the original time series is I(2), and that does not appear to be the case. Next, we apply our testing procedures to the first differenced time series, $\{\Delta X_t\}_{t=1}^T$. The bandwidth parameters are set to those applied in Section (ref), and in this part of the analysis, we use mean-adjusted statistics (Section (ref)). The results are reported in Table (ref), where we see that $\ddot\mathbbm{d}= 2$, $\widehat\mathbbm{d}_{\text{\tiny{LRS}}} =1$ (but recall that they cannot be zero), whereas $\widehat\mathbbm{d}_{\rm BU} = \widehat\mathbbm{d}_{\rm TD} = \widehat\mathbbm{d}_{\rm UD} = 0$ regardless of the choice of test statistic and level (except inv.VR at 10% level). From Figure (ref) and Table (ref), it seems unlikely that there are any I(2) stochastic trends.

We therefore proceed to apply our methodology to the level time series, $\{X_t\}_{t=1}^T$, to determine the number of I(1) stochastic trends. We now correct for a mean and linear time trend. The results are presented in Table (ref). For the eigenvalue ratio estimators we find $\ddot\mathbbm{d} = \widehat\mathbbm{d}_{\text{\tiny{LRS}}} = 1$. However, the sequential testing results are quite different. Although the inverse VR test results are reported only for $\mathbbm{d}_0=6,\ldots, 11$, the null hypotheses $\mathbbm{d}_0=0,1,\ldots , 5$ are all rejected at 1% level. Thus, we find that $\widehat\mathbbm{d}_{\rm BU} = 10$ at 5% level and $\widehat\mathbbm{d}_{\rm BU} = 9$ at 1% level, and we use this to set $\mathbbm{d}_{\max} = 15$. We report results for the VR tests for $\mathbbm{d}_0=6,7,\ldots,11$. The larger values of $\mathbbm{d}_0$ are rejected at 1% level for VR(2,1) and VR(1,0), but not for VR(2,0). This is consistent with the simulation findings in Table (ref). We then find $\widehat\mathbbm{d}_{\rm UD} = 8$ based on VR(1,0) at 5% level, $\widehat\mathbbm{d}_{\rm UD}=9$ based on VR(2,1) at 5% level or VR(1,0) at 1% level, and $\widehat\mathbbm{d}_{\rm UD}=10$ based on VR(2,1) at 1% level. On balance, we conclude that there are likely 9 stochastic trends. This means that the number of cointegrating relationships is about 2/3 of the dimension of the original data. Given this result, one should be careful in the application of cointegration rank tests that require a small number of cointegrating relationships (Remark (ref)).

Finally, for both empirical applications we investigate the robustness of $\widehat\mathbbm{d}_{\rm TD}$ to the initial hypothesis, $\mathbbm{d}_{\max}$, by computing $\widehat\mathbbm{d}_{\rm TD}$ using the VR(2,1), VR(2,0), and VR(1,0) tests at 5% level for a range of $\mathbbm{d}_{\max}$. In the yield example, $\widehat\mathbbm{d}_{\rm TD} = 2$ for $\mathbbm{d}_{\max}=2,\ldots ,20$ using either of the tests. In the labor index example in levels, $\widehat\mathbbm{d}_{\rm TD}$ is unchanged for $\mathbbm{d}_{\max} = 9,\ldots ,20$ using VR(2,1) and VR(1,0). The VR(2,0) test is sensitive to $\mathbbm{d}_{\max}$ as expected given the simulation results in Table (ref). The former two tests are remarkably robust to the choice of $\mathbbm{d}_{\max}$.

Conclusion

We have considered statistical inference on unit roots and cointegration for time series taking values in a Hilbert space of an arbitrarily large dimension or a subspace of possibly unknown dimension. This has wide applicability in practice; for example, in cointegrated vector time series of finite dimension, in a high-dimensional factor model that includes a finite number of nonstationary factors, in cointegrated curve-valued (or function-valued) time series, and in nonstationary dynamic functional factor models. We considered mainly determination of the dimension of the nonstationary subspace (number of common stochastic trends), but we also considered hypothesis testing on the stationary and nonstationary subspaces themselves.

We provided limit theory for general variance ratio-type statistics based on partial summation and/or differencing, and we demonstrated how to apply these variance ratio statistics to sequentially test for the dimension of the nonstationary subspace or to test hypotheses on the subspaces. Our theory covers many possible tests, and choosing among them is not trivial. Our preferred test is the VR(2,1) test due to several theoretical and practical advantages. Most importantly, it does not require the choice of any bandwidth parameters as it can be implemented only with estimated sample variance operators and not long-run variance operators. This is not only a practical advantage, but also implies a consistency property not shared by tests that need to estimate a long-run variance muller2007,muller2008. Also, critical values are readily available Breitung2002 by the isomorphism with $\mathbb{R}^\mathrm{K}$. To illustrate our methods, we included a small Monte Carlo simulation study, wherein the VR(2,1) test also showed superior performance, as well as two empirical illustrations to the term structure of interest rates and labor market indices, respectively.