3.1 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 (3.4) 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, that is, $\lambda _1\{A\}\leq \lambda _2\{A\} \leq \cdots $ . Given a kernel $\mathrm {k}_{\mathit {m}}$ satisfying Assumption 3, we define

(3.5) $$ \begin{align} c_m = \int_{-1}^1 \mathrm{k}_{\mathit{m}} (u) du = 2 \int_0^1 \mathrm{k}_{\mathit{m}} (u) du. \end{align} $$

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

Theorem 3.1. Suppose that Assumptions 1–3 hold and define $\widetilde {\mu }_j = n_T \mu _j$ , where $\{\mu _j\}_{j=1}^{\mathrm {K}}$ are the eigenvalues from (3.4) with $d_L \geq 2$ , $d_R = 1$ , and

(3.6) $$ \begin{align} n_T = T^{2d_L-2} (h_Lc_L + 1_{\{ a_L=0 \}})/(h_Rc_R + 1_{\{ a_R = 0\}}) , \end{align} $$

where $h_m,a_m$ , and $c_m$ are given in Assumption 3(ii) and (3.5). Then

(3.7) $$ \begin{align} \widetilde{\mu}_j \overset{\mathrm{d}}\rightarrow \widetilde{\lambda}_j &= \lambda_j \left\{ \left(\int W_{d_L,{\mathbb{d}_{\mathsf{N}}}} W_{d_L,{\mathbb{d}_{\mathsf{N}}}}'\right)^{-1} \int W_{1,{\mathbb{d}_{\mathsf{N}}}} W_{1,{\mathbb{d}_{\mathsf{N}}}}' \right\}\nonumber\\&\quad \text{ jointly for } j \leq {\mathbb{d}_{\mathsf{N}}}\ (\text{if }{\mathbb{d}_{\mathsf{N}}} \geq 1), \end{align} $$

(3.8) $$ \begin{align} \widetilde{\mu}_j & \overset{\mathrm{P}}\rightarrow \infty \quad \text{for } j \geq {\mathbb{d}_{\mathsf{N}}} +1. \end{align} $$

Note that the normalization factor (3.6) 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 3.1 shows that the eigenvalues in (3.4) have distinct asymptotic properties depending on whether j is greater than the dimension of $\mathcal H_{\mathsf {N}}$ or not, as motivated in Remark 3.4. That is, if we consider the vector $(\widetilde \mu _1 ,\ldots , \widetilde \mu _{\mathrm {K}})$ , the first ${\mathbb {d}_{\mathsf {N}}}$ elements converge jointly in distribution to the vector $(\tilde {\lambda }_1 , \ldots , \tilde {\lambda }_{{\mathbb {d}_{\mathsf {N}}}})$ , while the last $\mathrm {K} - {\mathbb {d}_{\mathsf {N}}}$ elements are divergent in probability. If $\mathrm {K} \geq \mathbb {d}_0$ , for any continuous map $\mathcal {F} : \mathbb {R}^{\mathbb {d}_0} \to \mathbb {R}$ , we have

$$ \begin{align*} \mathcal{F}(\{ \widetilde{\mu}_j\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{d}}\rightarrow \mathcal{F}(\{\widetilde{\lambda}_j\}_{j=1}^{\mathbb{d}_0}) \quad \text{under } H_0 \text{ of }(3.1). \end{align*} $$

If $\mathcal {F}$ additionally satisfies $\mathcal {F}(\{ \widetilde {\mu }_j\}_{j=1}^{\mathbb {d}_0}) \rightarrow _{\mathrm {P}} \infty $ under $H_1$ of (3.1), 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 {\mathrm {trace}}}$ , which are defined by

$$ \begin{align*} \mathcal{F}_{\max}(\{x_j\}_{j=1}^{\mathbb{d}_0}) = \max_{1\leq j\leq \mathbb{d}_0}\{x_j\} \quad\text{and}\quad \mathcal{F}_{\operatorname{\mathrm{trace}}} (\{x_j\}_{j=1}^{\mathbb{d}_0}) = \sum_{j=1}^{\mathbb{d}_0} x_j. \end{align*} $$

The following corollary delivers consistent tests of (3.1) based on VR( $d_L$ ,1).

Corollary 3.1. Consider the setup of Theorem 3.1 with $\mathrm {K} \geq \mathbb {d}_0$ . Under $H_0 : {\mathbb {d}_{\mathsf {N}}} = \mathbb {d}_0$ ,

(3.9) $$ \begin{align} \mathcal{F}_{\max}(\{{\widetilde{\mu}_j}\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{d}}\rightarrow \max_{1\leq j\leq \mathbb{d}_0}\{\lambda_{j}\{ \mathcal A \} \} , \quad \mathcal{F}_{\operatorname{\mathrm{trace}}}(\{ {\widetilde{\mu}_j}\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{d}}\rightarrow \sum_{j=1}^{\mathbb{d}_0}\lambda_{j}\{ \mathcal A \} , \end{align} $$

where $\mathcal A= (\int W_{d_L,\mathbb {d}_0} W_{d_L,\mathbb {d}_0}' )^{-1} \int W_{1,\mathbb {d}_0} W_{1,\mathbb {d}_0}'$ . Under $H_1 : {\mathbb {d}_{\mathsf {N}}} <\mathbb {d}_0$ ,

(3.10) $$ \begin{align} \mathcal{F}_{\max}(\{\widetilde{\mu}_j\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{P}}\rightarrow \infty , \quad \mathcal{F}_{\operatorname{\mathrm{trace}}}(\{\widetilde{\mu}_j\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{P}}\rightarrow \infty. \end{align} $$

The test statistics in Corollary 3.1 are only functions of $\{\widehat {P}_{\mathrm {K}}X_t\}_{t\geq 1}$ , which may be understood as a $\mathrm {K}$ -dimensional vector-valued time series (Remark 3.2), and the associated limiting distributions are functionals of $\mathbb {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 $\mathbb {d}_0$ and the test functional.

3.2 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 (3.4) 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 3.1, where $a_R = 0$ can be chosen (see Remarks 3.5 and 3.6).

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

Theorem 3.2. Suppose that Assumptions 1–3 hold with $a_R>0$ and define $\widetilde {\mu }_j = n_T \mu _j$ , where $\{\mu _j\}_{j=1}^{\mathrm {K}}$ are the eigenvalues from (3.4) with $d_L \geq 1 , d_R=0$ , and $n_T = T^{2d_L-1}(h_Lc_L + 1_{\{a_L=0\}})$ , where $h_L,a_L$ , and $c_L$ are given in Assumption 3(ii) and (3.5). Then

$$ \begin{align*} \widetilde{\mu}_j &\overset{\mathrm{d}}\rightarrow \lambda_j \left\{ \left( \int W_{d_L,{\mathbb{d}_{\mathsf{N}}}} W_{d_L,{\mathbb{d}_{\mathsf{N}}}}' \right)^{-1} \right\} \quad \text{jointly for } j \leq {\mathbb{d}_{\mathsf{N}}}\ (\text{if } {\mathbb{d}_{\mathsf{N}}} \geq 1), \\ \widetilde{\mu}_j &\overset{\mathrm{P}}\rightarrow \infty \quad \text{for } j \geq {\mathbb{d}_{\mathsf{N}}} +1. \end{align*} $$

The next corollary delivers consistent tests of (3.1) based on VR( $d_L$ ,0) with $d_L \geq 1$ .

Corollary 3.2. Consider the setup of Theorem 3.2 with $\mathrm {K} \geq \mathbb {d}_0$ . Under $H_0 : {\mathbb {d}_{\mathsf {N}}} = \mathbb {d}_0$ ,

(3.11) $$ \begin{align} \mathcal{F}_{\max}(\{{\widetilde{\mu}_j}\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{d}}\rightarrow \max_{1\leq j\leq \mathbb{d}_0} \{\lambda_{j}\{ \mathcal A \} \}, \quad \mathcal{F}_{\operatorname{\mathrm{trace}}}(\{ {\widetilde{\mu}_j}\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{d}}\rightarrow \sum_{j=1}^{\mathbb{d}_0}\lambda_j \{ \mathcal A \} , \end{align} $$

where $\mathcal A = ( \int W_{d_L,\mathbb {d}_0} W_{d_L,\mathbb {d}_0}' )^{-1}$ . Under $H_1 : {\mathbb {d}_{\mathsf {N}}} < \mathbb {d}_0$ ,

(3.12) $$ \begin{align} \mathcal{F}_{\max}(\{\widetilde{\mu}_j\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{P}}\rightarrow \infty , \quad \mathcal{F}_{\operatorname{\mathrm{trace}}}(\{\widetilde{\mu}_j\}_{j=1}^{\mathbb{d}_0}) \overset{\mathrm{P}}\rightarrow \infty. \end{align} $$

3.3 Inverse VR Tests

We showed in Theorem 3.1 that the first ${\mathbb {d}_{\mathsf {N}}}$ eigenvalues $\{\mu _j\}_{j=1}^{{\mathbb {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}-{\mathbb {d}_{\mathsf {N}}}$ eigenvalues, $\{\mu _j\}_{j={\mathbb {d}_{\mathsf {N}}} +1}^{\mathrm {K}}$ . This enables us to examine the inverse testing problem,

(3.13) $$ \begin{align} H_0: {\mathbb{d}_{\mathsf{N}}} = \mathbb{d}_0 \quad\text{vs}\quad H_1: {\mathbb{d}_{\mathsf{N}}} \geq \mathbb{d}_0 + 1, \end{align} $$

for some $\mathbb {d}_0 \geq 0$ . Unlike in the testing problem (3.1), we do not require any prior information on a reasonable upper bound for ${\mathbb {d}_{\mathsf {N}}}$ .

Related to testing (3.13), we will replace Assumption 2 with the following high-level assumption, which allows the possibility that $\mathrm {K} < {\mathbb {d}_{\mathsf {N}}}$ .

Theorem 3.3. Suppose that Assumptions 1–3 hold with $a_R>0$ and define $\widetilde {\mu }_j = n_T \mu _j$ , where $\{\mu _j\}_{j=1}^{\mathrm {K}}$ are the eigenvalues from (3.4) with $d_L \geq 2$ , $d_R = 1$ , and $n_T = T^{2d_L -3}( h_Lc_L + 1_{\{a_L=0\}})$ , where $h_L,a_L$ , and $c_L$ are given in Assumption 3(ii) and (3.5). Then

(3.14) $$ \begin{align} (\widetilde{\mu}_j)^{-1} &\overset{\mathrm{d}}\rightarrow \lambda_{\mathrm{K} -j+1} \{ \mathcal A \} \quad \text{jointly for } j = {\mathbb{d}_{\mathsf{N}}} +1 , \ldots , \mathrm{K}\ (\text{if } {\mathbb{d}_{\mathsf{N}}} < \mathrm{K}), \end{align} $$

(3.15) $$ \begin{align} (\widetilde{\mu}_j)^{-1} &\overset{\mathrm{P}}\rightarrow \infty \quad \text{for } j\leq {\mathbb{d}_{\mathsf{N}}}\ (\text{if } {\mathbb{d}_{\mathsf{N}}} \geq 1), \end{align} $$

where $\mathcal A = \int B_{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}} B_{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}}' - 1_{\{ {\mathbb {d}_{\mathsf {N}}} \geq 1 \}} \int B_{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}} W_{d_L,{\mathbb {d}_{\mathsf {N}}}}' (\int W_{d_L,{\mathbb {d}_{\mathsf {N}}}} W_{d_L,{\mathbb {d}_{\mathsf {N}}}}' )^{-1} \int W_{d_L,{\mathbb {d}_{\mathsf {N}}}}B_{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}}'$ and $B_{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}}$ is a $(\mathrm {K}-{\mathbb {d}_{\mathsf {N}}})$ -dimensional standard Brownian motion which is independent of $W_{d_L,{\mathbb {d}_{\mathsf {N}}}}$ .

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

When $\mathrm {K}> \mathbb {d}_0$ , and given a proper test functional $\mathcal {F}$ , the asymptotic results in Theorem 3.3 can be used to deliver consistent tests of (3.13) as follows.

Corollary 3.3. Consider the setup of Theorem 3.3 with $\mathrm {K}> \mathbb {d}_0$ . Under $H_0 : {\mathbb {d}_{\mathsf {N}}}=\mathbb {d}_0$ ,

(3.16) $$ \begin{align} \mathcal{F}_{\max} ( \{ \widetilde{\mu}_j^{-1} \}_{j=\mathbb{d}_0+1}^{\mathrm{K}} ) \overset{\mathrm{d}}\rightarrow \max_{1\leq j\leq \mathrm{K}-\mathbb{d}_0} \{\lambda_j\{ \mathcal A \} \}, \quad \mathcal{F}_{\operatorname{\mathrm{trace}}} ( \{ \widetilde{\mu}_j ^{-1} \}_{j=\mathbb{d}_0+1}^{\mathrm{K}}) \overset{\mathrm{d}}\rightarrow \sum_{j=1}^{\mathrm{K}-\mathbb{d}_0} \lambda_j \{ \mathcal A \}, \end{align} $$

where $\mathcal A = \int B_{\mathrm {K}-\mathbb {d}_0} B_{\mathrm {K}-\mathbb {d}_0}' - 1_{\{ \mathbb {d}_0\geq 1 \}} \int B_{\mathrm {K}-\mathbb {d}_0} W_{d_L,\mathbb {d}_0}' (\int W_{d_L,\mathbb {d}_0} W_{d_L,\mathbb {d}_0}' )^{-1} \int W_{d_L,\mathbb {d}_0} B_{\mathrm {K}-\mathbb {d}_0}'$ . Under $H_1 : {\mathbb {d}_{\mathsf {N}}}>\mathbb {d}_0$ ,

(3.17) $$ \begin{align} \mathcal{F}_{\max} ( \{ \widetilde{\mu}_j ^{-1} \}_{j=\mathbb{d}_0+1}^{\mathrm{K}}) \overset{\mathrm{P}}\rightarrow \infty , \quad \mathcal{F}_{\operatorname{\mathrm{trace}}} ( \{ \widetilde{\mu}_j ^{-1} \}_{j=\mathbb{d}_0+1}^{\mathrm{K}}) \overset{\mathrm{P}}\rightarrow \infty. \end{align} $$

The test functionals described in Corollary 3.3 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 3.1 concerning the testing problem (3.1), we allow $\mathbb {d}_0=0$ in Corollary 3.3. 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 $\mathbb {d}_0=0$ (see Section S.2.3 of the Supplementary Material for a detailed discussion of related procedures).

3.4 Estimation of Slack Extractor

To implement the VR-based tests in practice, a slack extractor $\widehat {P}_{\mathrm {K}}$ satisfying either Assumption 2 or 2’ is needed. As shown in Sections 3.1–3.3, both the VR and inverse VR tests require Assumption 2 to hold for $\mathrm {K} \geq {\mathbb {d}_{\mathsf {N}}}$ , and when $\mathrm {K} < {\mathbb {d}_{\mathsf {N}}}$ the inverse VR tests require instead that Assumption 2’ 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 (3.2), 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 4 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 (3.2). Specifically, it requires that the first $\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}$ eigenvalues of $\Lambda _{\mathsf {S},R}$ are nonzero and that the $(\mathrm {K}-{\mathbb {d}_{\mathsf {N}}})$ -th eigenvalue is different from the next eigenvalue, $\tau _{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}+1}$ . The role of the requirement $\tau _{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}} \neq \tau _{\mathrm {K}-{\mathbb {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 5 and Theorem 3.5).

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 Nielsen et al., Reference Nielsen, Seo and Seong2023), $\Lambda _{\mathsf {S},R}$ generally has infinitely many nonzero eigenvalues, and Assumption 4 allows $\mathrm {K}$ to be any arbitrary finite integer as long as $\tau _{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}} \neq \tau _{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}+1}$ . In a finite-dimensional space, a reasonable upper bound on ${\mathbb {d}_{\mathsf {N}}}$ will be relevant to choose $\mathrm {K}$ to satisfy Assumption 4.

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

(3.18) $$ \begin{align} \widehat{\Lambda}_{1,P} = \widehat{\Lambda}_{1,P} (a_R) = \widehat{\Gamma}_{1,0}+ 1_{\{a_R>0\}}\sum_{s=1}^{T-1} \mathrm{k}_P (s/h_P)(\widehat{\Gamma}_{1,s} + \widehat{\Gamma}_{1,-s}) \end{align} $$

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 2 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 3.11). The advantage is that the estimator $\widehat {P}_{\mathrm {K}}$ described in Theorem 3.4 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 3.4. Thus, Theorem 3.4 provides an easy-to-implement way to obtain an estimated projection operator $\widehat {P}_{\mathrm {K}}$ that satisfies the high-level conditions in Assumptions 2 and 2’.

The condition $\tau _{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}} \neq \tau _{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}+1}$ in Assumption 4 is generally required for the VR( $d_L,d_R$ )-based tests and also required for the inverse VR test when $\mathbb {d}_0$ is close to ${\mathbb {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 3.4, the ordered eigenvalues $\{\widehat {\tau }_{j}\}_{j \geq {\mathbb {d}_{\mathsf {N}}}+1}$ of $\widehat {\Lambda }_{1,P}$ converge in probability to $\{ \tau _j \}_{j\geq 1}$ . That is, the first ${\mathbb {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={\mathbb {d}_{\mathsf {N}}}+1,\ldots $ that correspond to eigenvalues $\tau _j$ for $j=1,\ldots $ of $\Lambda _{\mathsf {S},R}$ .

We may thus avoid violation of the condition $\tau _{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}} \neq \tau _{\mathrm {K}-{\mathbb {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 5 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 2, and the coefficients $\{\Phi _j\}_{j\geq 1}$ are required to decay at a faster rate than in Section 2. The latter can be seen from the fact that $\|\sum _{j=m}^\infty \tilde {\Phi }_j \|_{\operatorname {\mathrm {op}}} = {\mathcal {O}} (m^{-\gamma })$ for some $\gamma>1$ implies the summability condition in (2.2). Assumption 5(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 5(iv) requires that $h_P$ is not too small relative to T. Nonetheless, even if Assumption 5 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}-{\mathbb {d}_{\mathsf {N}}}} \neq \tau _{\mathrm {K}-{\mathbb {d}_{\mathsf {N}}}+1}$ based on the sample eigenvalues $\widehat \tau _j$ . The following theorem establishes the desired result.

Theorem 3.5. Suppose that Assumption 1 holds. If $a_R>0,$ then assume also that $\mathrm {k}_P(\cdot )$ and $h_P$ in (3.18) satisfy Assumption 3 and that Assumption 5 holds. If $\mathrm {K}> {\mathbb {d}_{\mathsf {N}}}$ , then

$$ \begin{align*} (\widehat{\tau}_{\mathrm{K}} - \widehat{\tau}_{\mathrm{K}+1}) - ( \tau_{\mathrm{K}-{\mathbb{d}_{\mathsf{N}}}} - \tau_{\mathrm{K}-{\mathbb{d}_{\mathsf{N}}}+1} ) = \mathcal{O}_{\mathrm{P}} ( M_T^{-1/2}), \end{align*} $$

where $M_T = T/h_P$ if $a_R>0$ and $M_T = T$ if $a_R = 0$ .

From Theorem 3.5, we deduce that, for any $\eta _T$ satisfying $\eta _T M_T^{-1/2} \to 0$ and $\eta _T \to \infty $ ,

$$ \begin{align*} \eta_T (\widehat{\tau}_{\mathrm{K}}-\widehat{\tau}_{\mathrm{K}+1}) &\overset{\mathrm{P}}\rightarrow 0 \phantom{\infty} \quad \text{if} \quad \tau_{\mathrm{K}-{\mathbb{d}_{\mathsf{N}}}} = \tau_{\mathrm{K}-{\mathbb{d}_{\mathsf{N}}}+1}, \\ \eta_T (\widehat{\tau}_{\mathrm{K}}-\widehat{\tau}_{\mathrm{K}+1}) &\overset{\mathrm{P}}\rightarrow \infty \phantom{0} \quad \text{if} \quad \tau_{\mathrm{K}-{\mathbb{d}_{\mathsf{N}}}} \neq \tau_{\mathrm{K}-{\mathbb{d}_{\mathsf{N}}}+1}. \end{align*} $$

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

Remark 3.12. In our proof of Theorem 3.5, Assumption 5 is used to show that, when $a_R>0$ ,

(3.19) $$ \begin{align} \| T^{-1}P_{\mathsf{S}} \widehat{\Lambda}_{1,P} P_{\mathsf{S}} - \Lambda_{\mathsf{S},R} \|_{\operatorname{\mathrm{op}}} = \mathcal{O}_{\mathrm{P}} ((T/h_P)^{-1/2}) \end{align} $$

(see Berkes et al., Reference Berkes, Horváth and Rice2016, Theorems 2.2 and 2.3). The desired result then follows from (3.19). Thus, any other conditions that imply (3.19) can replace Assumption 5 in Theorem 3.5 (see also Lemma 4.1 of Rice and Shang, Reference Rice and Shang2017).

Remark 3.13. A potentially useful result that may be deduced from our proofs of Theorems 3.4 and 3.5 is that the eigenvalues $\{\widehat {\tau }_j\}_{j\geq 1}$ of $\widehat {\Lambda }_{1,P}$ satisfy

(3.20) $$ \begin{align} \frac{\widehat{\tau}_j}{\widehat{\tau}_{j+1}} \overset{\mathrm{P}}\rightarrow \infty \text{ if } j={\mathbb{d}_{\mathsf{N}}} \quad\text{and}\quad \frac{\widehat{\tau}_j}{\widehat{\tau}_{j+1}} =\mathcal{O}_{\mathrm{P}} (1) \text{ if } j\neq{\mathbb{d}_{\mathsf{N}}}. \end{align} $$

This result is analogous to that in Theorem 3.3 of Chang et al. (Reference Chang, Kim and Park2016) for the sample variance operator of $\{X_t\}_{t\geq 1}$ (see also Theorem 3.2 of Li et al., Reference Li, Robinson and Shang2023). Based on (3.20), we may in practice easily obtain a reasonable upper bound on ${\mathbb {d}_{\mathsf {N}}}$ , say $\mathbb {d}_{\max }$ , by computing the ratio given in (3.20) and set $\mathbb {d}_{\max } = \max _{1\leq j \leq \bar {s}} \{{\widehat {\tau }_j}/{\widehat {\tau }_{j+1}}\} + k$ for some large integer $\bar {s}$ and a small nonnegative integer k. Of course, in view of the asymptotic properties of $\widehat {\tau }_j$ , $\max _{1\leq j \leq \bar {s}} \{{\widehat {\tau }_j}/{\widehat {\tau }_{j+1}}\}$ is itself a consistent estimator of ${\mathbb {d}_{\mathsf {N}}}$ . The latter is a slight modification of the eigenvalue ratio estimator of Li et al. (Reference Li, Robinson and Shang2023). However, we found that the finite-sample properties of this estimator tend to be worse than those of a different eigenvalue ratio-based estimator to be discussed in Section 4.2. Nonetheless, either eigenvalue ratio-based estimator could be useful as a simple way to construct a reasonable upper bound on ${\mathbb {d}_{\mathsf {N}}}$ in practice.

3.5 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

(3.21) $$ \begin{align} Y_t = \zeta_1 + \zeta_2 t + X_t, \quad t \geq1, \end{align} $$

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

(3.22) $$ \begin{align} U_t^{(1)} = Y_t - \bar{Y}, \quad U_t^{(2)} = Y_t - \bar{Y} - ( t-\bar{t} ) \frac{\sum_{t=1}^T ( t-\bar{t} ) Y_t}{\sum_{t=1}^T ( t-\bar{t} )^2}, \end{align} $$

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$ (Kokoszka and Young, Reference Kokoszka and Young2016).

Our ADI VR tests can be adjusted to accommodate deterministic terms by replacing the sample operators computed from $\{X_t\}_{t\geq 1}$ with the corresponding operators computed from $\{U_t^{(1)}\}_{t\geq 1}$ or $\{U_t^{(2)}\}_{t\geq 1}$ . When the model for $Y_t$ is given by (3.21) 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 3.1 and 3.2 as in the existing literature (e.g., Nielsen et al., Reference Nielsen, Seo and Seong2023, Sect. 3.4 and Seo, Reference Seo2024, Sect. S3). For example, the extension for the VR $(d_L,1)$ test can easily be done as in Theorem 3 of Nielsen et al. (Reference Nielsen, Seo and Seong2023). More specifically, in Theorems 3.1 and 3.2, $W_{1,{\mathbb {d}_{\mathsf {N}}}}$ needs to be replaced by ${\mathbb {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,{\mathbb {d}_{\mathsf {N}}}}$ needs to change accordingly. Similar changes have to be made in Theorem 3.3, but where, as may be deduced from our proof of Theorem 3.3, $B_{\mathrm {K} - {\mathbb {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,{\mathbb {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 3.4 and 3.5 do not require any changes once $\widehat {\Lambda }_{1,P}$ is constructed from the relevant residuals (our proofs of Theorems 3.4 and 3.5 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.

4.1 Estimation via Sequential Testing

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

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

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

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

The proof of Theorem 4.1 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 {\mathbb {d}}_{\mathrm {TD}} = {\mathbb {d}_{\mathsf {N}}})~\to ~1$ and $\mathrm {P}(\widehat {\mathbb {d}}_{\mathrm {BU}} = {\mathbb {d}_{\mathsf {N}}})~\to ~1$ .

Remark 4.2. To implement the proposed sequential procedures under the low-level conditions discussed in Section 3.4, in each step we require Assumption 4 with $\mathrm {K} = \mathrm {K}(\mathbb {d}_0) \geq \mathbb {d}_0$ for the TD procedure or $\mathrm {K} = \mathrm {K}(\mathbb {d}_0)> \mathbb {d}_0$ for the BU procedure. In either case, for each $\mathbb {d}_0$ , our tests do not require knowledge of the dimension of the (possibly strict sub)space in which $X_t$ takes values, but do require $\Lambda _{\mathsf {S},R}$ to have $\mathrm {K} - {\mathbb {d}_{\mathsf {N}}}$ nonzero eigenvalues. It may be deduced from Theorem 3.4 and Assumption 4 that the TD procedures based on VR( $d_L,d_R$ ) tests require a stronger condition on $\Lambda _{\mathsf {S},R}$ than the BU procedure based on inverse VR tests. To see this, suppose that a common choice of $\mathrm {K}$ (as a function of $\mathbb {d}_0$ ) is used for both of the procedures, such as $\mathrm {K}=\mathbb {d}_0+m$ for some $m \geq 1$ . Because the BU procedure starts from $\mathbb {d}_0 \leq {\mathbb {d}_{\mathsf {N}}}$ and the TD procedure from $\mathbb {d}_0 \geq {\mathbb {d}_{\mathsf {N}}}$ , the number of nonzero eigenvalues that we require for the TD procedure is always larger than that required for the BU procedure; this is illustrated for $m=2$ in Figure 1. Moreover, the number of nonzero eigenvalues of $\Lambda _{\mathsf {S},R}$ required by the TD procedure is always positive when $\mathbb {d}_0> {\mathbb {d}_{\mathsf {N}}}$ , and this number increases as $\mathbb {d}_0$ gets larger. This suggests that it is important to have a reasonable starting point $\mathbb {d}_{\max } (\geq {\mathbb {d}_{\mathsf {N}}} )$ for the TD procedure. Similarly, it should also be noted that, all else equal, $\Lambda _{\mathsf {S},R}$ is required to have more nonzero eigenvalues as $\mathrm {K}$ increases for both procedures, so it is preferable to choose $\mathrm {K}$ not too much larger than $\mathbb {d}_0$ . A similar suggestion for the choice of $\mathrm {K}$ can be found in Nielsen et al. (Reference Nielsen, Seo and Seong2023), but their argument was based on computational advantages in a more restrictive setting than ours. Our finding thus gives a theoretical justification for their choice in our more general setting.

Figure 1

Eigenvalue requirements for the proposed tests, $\mathrm {K}=\mathrm {K}(\mathbb {d}_0)=\mathbb {d}_0+2$ .

Figure 1 Long description

A legend at the top shows a thick red line representing the Required number of nonzero eigenvalues. Below are three horizontal axes.

* The first axis is labeled V R (d sub L, d sub R). It shows a thick red line segment starting at d sub N, passing through d sub 0, and ending at K. The points are ordered from left to right as d sub N, then d sub 0, then K.

* The second axis is labeled Inverse V R with K greater than d sub N. It shows a thick red line segment starting at d sub N and ending at K. On this axis, d sub 0 is located to the left of d sub N, and K is located to the right of d sub N.

* The third axis is labeled Inverse V R with K less than or equal to d sub N. There is no red line segment on this axis. The points are ordered from left to right as d sub 0, then K, then d sub N.

As discussed in Remarks 4.1 and 4.2, 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 $\mathbb {d}_{\max }$ on ${\mathbb {d}_{\mathsf {N}}}$ and, specifically, the choices of $\mathrm {K}$ for the BU procedure to consistently estimate ${\mathbb {d}_{\mathsf {N}}}$ (with shrinking significance levels as discussed above) depend only on $\mathbb {d}_0$ (not ${\mathbb {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 {\mathbb {d}}_{\mathrm {BU}}$ and use this to specify a data-driven value of the upper bound as $\mathbb {d}_{\max } = \widehat {\mathbb {d}}_{\mathrm {BU}} +m_{\mathbb {d}}$ for an integer $m_{\mathbb {d}} \geq 0$ . We then apply the TD procedure initiated at this upper bound to obtain the UD hybrid estimate, $\widehat {\mathbb {d}}_{\mathrm {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} = \mathbb {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 6.

4.2 Eigenvalue Ratio Estimator

Recent work concerning vector-valued time series provides 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 (e.g., Zhang and Chan, Reference Zhang and Chan2018; Zhang et al., Reference Zhang, Robinson and Yao2019; Chang, Cheng, and Yao, Reference Chang, Cheng and Yao2021). In a functional time series setting, the eigenvalue ratio estimator proposed by Li et al. (Reference Li, Robinson and Shang2020, Reference Li, Robinson and Shang2023) 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 3.13.

Alternatively, based on the asymptotic results in Section 3, we may use the generalized eigenvalue problem (3.4) to construct an eigenvalue ratio estimator for ${\mathbb {d}_{\mathsf {N}}}$ . Let $\mu _j$ be an eigenvalue of (3.4). As we prove in Theorem 4.2 below, we can use Theorems 3.1–3.3 and arguments from their proofs to show that $\mu _j^{-1}\mu _{j+1} = \mathcal {O}_{\mathrm {P}} (1)$ if $j < {\mathbb {d}_{\mathsf {N}}}$ (resp. $j> {\mathbb {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={\mathbb {d}_{\mathsf {N}}}$ . Hence, the ratio $\mu _j^{-1}\mu _{j+1}$ will be maximized at $j={\mathbb {d}_{\mathsf {N}}}$ . This motivates the estimator

(4.1) $$ \begin{align} \ddot{\mathbb{d}} = \underset{1\leq j\leq \mathbb{d}_{\max}}{\arg\max} \left\{\frac{\mu_{j+1}}{\mu_{j}}\right\}, \end{align} $$

where $\mathbb {d}_{\max }$ is a pre-specified integer indicating an upper bound on ${\mathbb {d}_{\mathsf {N}}}$ . We note that these eigenvalue ratio-type estimators are only applicable when ${\mathbb {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 (3.13) with $\mathbb {d}_0=0$ using the tests in Corollary 3.3.

The estimator in (4.1) is very similar to the eigenvalue ratio estimator proposed by Li et al. (Reference Li, Robinson and Shang2020, Reference Li, Robinson and Shang2023), see also Remark 3.13, 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 (3.4). Second, their estimator requires a tuning parameter, say $\kappa>0$ , to deal with estimation error associated with small eigenvalues; specifically, they set $\widehat {\mathbb {d}}_{\text { {LRS}}} = \arg \min _{1\leq j\leq \mathbb {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 (4.1) does not require such a tuning parameter.

7.1 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 are from the U.S. Department of Treasury available at https://home.treasury.gov/ and spans January 1984 to December 2018. At each time t, we observe a corporate bond with 200 different maturities ranging from six months to 100 years. This data were studied by Barigozzi and Trapani (Reference Barigozzi and Trapani2022), 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 6.

Figure 2a shows the time series of HQM yield curves (observed at 200 different maturities) and Figure 2b 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 Barigozzi and Trapani (Reference Barigozzi and Trapani2022) for the data span of January 1985 to September 2017. Thus, we apply the trend-adjusted testing procedures in Section 3.5.

Figure 2

HQM yield curves.

Figure 2 Long description

The image consists of two panels labeled a and b.

Panel a, titled Yields at 200 maturities, is a three-dimensional surface plot. The vertical axis represents yield values from 0 to 15. The horizontal depth axis is labeled Maturities, ranging from 0 to 100. The horizontal width axis is labeled Time, with markers for 1991 Q 1, 2001 Q 1, 2010 Q 1, and 2018 Q 1. The surface is color-coded with a gradient from dark purple at low yields to bright yellow-green at high yields. The plot shows a complex, undulating surface that generally trends downward over time with significant volatility and peaks in the early 1990s and mid-2000s.

Panel b, titled Yields at a few different maturities, is a two-dimensional line graph. The vertical axis shows yield values from 5 to 15. The horizontal axis represents Time from 1991 Q 1 to 2018 Q 1. Four lines are plotted:

- A solid blue line representing 3 years maturity, which is the lowest and most volatile curve.

- A dashed orange line representing 6 years maturity, positioned slightly above the 3-year line.

- A solid green line representing 10 years maturity, following a similar trend but higher.

- A dotted red line representing 30 years maturity, which is the highest and smoothest curve.

All four curves show a long-term downward trend with synchronized peaks around 1991, 2000, and 2008, followed by a sharp decline after 2008.

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

Table 3

Empirical results for HQM yield curve data.

Table 3 Long description

The table is divided into two primary sections.

Section 1: Eigenvalue ratio estimates.

The columns are labeled j equals 1 through j equals 7.

Row 1: mu sub j plus 1 over mu sub j values are 2.89, 3.00, 2.02, 1.91, 1.15, 3.25, and 2.15.

Row 2: kappa-hat sub j plus 1 over kappa-hat sub j values are 0.42, 0.14, 0.08, 0.29, 0.00, 1, and 1.

Section 2: Variance ratio test statistics.

The columns are labeled m sub 0 equals 0 through m sub 0 equals 7.

Row 1: Inv.V R values range from 0.94 double-asterisk at m sub 0 equals 0 to 0.03 at m sub 0 equals 7. Significance decreases as m sub 0 increases.

Row 2: V R (2,1) values range from 151.61 at m sub 0 equals 1 to 23,154.51 double-asterisk at m sub 0 equals 7.

Row 3: V R (2,0) values range from 7,007.65 at m sub 0 equals 1 to 6,936,298.27 double-asterisk at m sub 0 equals 7.

Row 4: V R (1,0) values range from 24.25 at m sub 0 equals 1 to 3,825.57 double-asterisk at m sub 0 equals 7.

Footnote: Dagger, asterisk, and double-asterisk indicate significance at the 10 percent, 5 percent, and 1 percent levels, respectively.

Note: $^{\dagger }$ , $^{*}$ , and $^{**}$ indicate significance at 10%, 5%, and 1% level, respectively.

In the second part of Table 3, we report the VR and inverse VR test statistics for a range of $\mathbb {d}_0$ . The inverse VR tests imply that $\widehat {\mathbb {d}}_{\mathrm {BU}} = 3$ (at 1% level) or $\widehat {\mathbb {d}}_{\mathrm {BU}} = 6$ (at 5% level). Given this result, the choice $\mathbb {d}_{\max } = 7$ seems like a reasonable input to the UD procedures. The latter give either $\widehat {\mathbb {d}}_{\mathrm {UD}} = 2$ or $\widehat {\mathbb {d}}_{\mathrm {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 {\mathbb {d}} = 3$ . Overall, the evidence thus supports ${\mathbb {d}_{\mathsf {N}}} = 3$ stochastic trends.

7.2 The Nelson–Siegel Term Structure Model

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

$$ \begin{align*} \gamma_0 &= 1, \quad \gamma_1 = \gamma_1 (\varsigma, \tau) = \varsigma (1-e^{-\tau/\varsigma})/\tau , \\ \gamma_2 &= \gamma_2 (\varsigma, \tau) = \varsigma (1-e^{-\tau/\varsigma})/\tau - e^{-\tau/\varsigma} , \end{align*} $$

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 3. 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 are annualized, we set $\varsigma = 1.37$ following Diebold and Li (Reference Diebold and Li2006).

We first test the hypothesis $H_0$ in (5.3) with $\mathcal H_0 = \operatorname {\mathrm {span}} \{ \gamma _0 , \gamma _1, \gamma _2 \}$ using Theorem 5.1(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 (5.3) 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 (5.2). From our earlier results in Table 3, we set the dimension of the latter to ${\mathbb {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 {\mathrm {span}} \{ \gamma _0 , \gamma _1, \gamma _2 \} \supseteq \mathcal H_{\mathsf {N}}$ in (5.2) is very interesting, and its non-rejection would be the strongest conclusion that we could obtain. However, implementing the test as in Theorem 5.1(ii) with $\mathbb {d}_0 = 0$ and $p_0 = 3$ , the inverse VR statistic is $1.47$ (see Table 4) and rejects the null at 1% level.

Table 4

Test results for $( I - P_{\mathcal H_0})X_t$ with $\mathcal H_0 = \operatorname {\mathrm {span}}\{\gamma _0,\gamma _1 ,\gamma _2 \}$ and $\varsigma =1.37$ .

Table 4 Long description

The table consists of five columns and five rows.

Header Row: The first column is labeled Statistic. The subsequent columns are d sub 0 = 0, d sub 0 = 1, d sub 0 = 2, and d sub 0 = 3.

Row 1: Statistic Inv.V R. Values are 1.47 double asterisk for d sub 0 = 0; 0.27 double asterisk for d sub 0 = 1; 0.24 double asterisk for d sub 0 = 2; and 0.12 asterisk for d sub 0 = 3.

Row 2: Statistic V R (2,1). The d sub 0 = 0 cell is empty. Values are 113.48 for d sub 0 = 1; 629.06 dagger for d sub 0 = 2; and 1,577.24 asterisk for d sub 0 = 3.

Row 3: Statistic V R (2,0). The d sub 0 = 0 cell is empty. Values are 6,842.81 for d sub 0 = 1; 39,832.64 asterisk for d sub 0 = 2; and 123,243.78 double asterisk for d sub 0 = 3.

Row 4: Statistic V R (1,0). The d sub 0 = 0 cell is empty. Values are 72.34 double asterisk for d sub 0 = 1; 167.27 double asterisk for d sub 0 = 2; and 346.79 double asterisk for d sub 0 = 3.

Footnote: Dagger, asterisk, and double asterisk indicate significance at the 10 percent, 5 percent, and 1 percent levels, respectively.

Note: $^{\dagger }$ , $^{*}$ , and $^{**}$ indicate significance at 10%, 5%, and 1% level, respectively.

We proceed to test $H_0$ in (5.1), 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 5.1(i), the inverse VR statistics with $\mathbb {d}_0 = 2$ are $0.24$ , $0.21$ , and $0.22$ for $\mathcal H_0 = \operatorname {\mathrm {span}} \{ \gamma _0 \}$ , $\mathcal H_0 = \operatorname {\mathrm {span}} \{ \gamma _1 \}$ , $\mathcal H_0 = \operatorname {\mathrm {span}} \{ \gamma _2 \}$ , respectively, in each case rejecting the null at 1% level. For two factors, the inverse VR statistics with $\mathbb {d}_0 = 1$ are $0.40$ and $0.38$ for $\mathcal H_0 = \operatorname {\mathrm {span}} \{ \gamma _0 , \gamma _1 \}$ and $\mathcal H_0 = \operatorname {\mathrm {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, that is, that any one factor can be interpreted as one of the common stochastic trends.

Finally, following Remark 5.1, 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 {\mathrm {span}} \{ \gamma _0 , \gamma _1, \gamma _2 \}$ , is of smaller dimension than that of $X_t$ itself. Letting the former be denoted $\mathbb {d}_{\mathcal H_0}$ , we apply the sequential procedures in (5.4) and (5.5). The results are reported in Table 4. The BU procedure yields $\widehat {\mathbb {d}}_{\mathcal H_0 , \mathrm {BU}} = 3$ at 1% level and rejects $\mathbb {d}_0 = 3$ at 5% level. Applying the TD sequential procedure starting from $\mathbb {d}_{\max } = 3$ , we find $\widehat {\mathbb {d}}_{\mathcal H_0 , \mathrm {TD}} = 2$ with the preferred VR(2,1) test at 5% level (but $\widehat {\mathbb {d}}_{\mathcal H_0 , \mathrm {TD}} = 1$ at 10% level with the VR(2,1) test, $\widehat {\mathbb {d}}_{\mathcal H_0 , \mathrm {TD}} = 1$ at 5% level with the VR(2,0) test, and $\widehat {\mathbb {d}}_{\mathcal H_0 , \mathrm {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.

7.3 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 (McCracken and Ng, Reference McCracken and Ng2016). 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 McCracken and Ng (Reference McCracken and Ng2016). We apply log-transformation to the data.Footnote ¹ The time span used in this analysis is January 1990 to December 2019.

In Figure 3, 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^{\mathrm {th}}$ score processes (see Chang et al., Reference Chang, Kim and Park2016 and Section 5 of Nielsen et al., Reference Nielsen, Seo and Seong2023). It is clear from Figure 3 that there is some nonstationary behavior (e.g., $j=1$ ), but this seems to eventually disappear leaving stationary processes (e.g., $j=10$ ).

Figure 3

Time series of $\langle X_t,\widehat {\nu }_j \rangle $ .

Note: Time series of $\langle X_t, \widehat {\nu }_j \rangle $ for some selected eigenvectors $\widehat {\nu }_j$ of $\widehat \Lambda _{0,R}$ .

Figure 3 Long description

A three-panel figure showing time series data. All panels share a common X axis representing time from January 1991 to January 2019, with major ticks at 1991 forward slash 01, 2001 forward slash 01, 2010 forward slash 01, and 2019 forward slash 01.

* Panel a, labeled with the inner product of X sub t and nu-hat sub 1, shows a smooth, oscillating blue line. It features three distinct peaks, with the highest peak occurring around 2010 reaching a value above 2 on the Y axis. The Y axis ranges from below 0 to above 2.

* Panel b, labeled with the inner product of X sub t and nu-hat sub 5, shows a more volatile blue line with higher frequency fluctuations. The Y axis ranges from approximately minus 0.2 to 0.2. There is a notable peak around 2005 and a sharp dip just before 2010.

* Panel c, labeled with the inner product of X sub t and nu-hat sub 10, shows the most volatile data with rapid, high-frequency oscillations around a zero baseline. The Y axis ranges from minus 0.1 to 0.1. The fluctuations appear as dense noise compared to the smoother trends in panel a.

McCracken and Ng (Reference McCracken and Ng2016) suggested that some of these labor market indices may be I(2) variables.Footnote ² To examine the existence of this type of higher-order stochastic trends, we first display in Figure 4 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 7.1, and in this part of the analysis, we use mean-adjusted statistics (Section 3.5). The results are reported in Table 5, where we see that $\ddot {\mathbb {d}}= 2$ , $\widehat {\mathbb {d}}_{\text { {LRS}}} =1$ (but recall that they cannot be zero), whereas $\widehat {\mathbb {d}}_{\mathrm {BU}} = \widehat {\mathbb {d}}_{\mathrm {TD}} = \widehat {\mathbb {d}}_{\mathrm {UD}} = 0$ regardless of the choice of test statistic and level (except inv.VR at 10% level). From Figure 4 and Table 5, it seems unlikely that there are any I(2) stochastic trends.

Table 5

Empirical results for labor indices in first differences

Table 5 Long description

The table is titled Empirical results for labor indices in first differences. It is organized into two main horizontal sections.

Section 1: Eigenvalue ratio estimates.

The columns are labeled j = 1, j = 2, j = 3, j = 4, and j = 5.

* Row 1 (mu sub j+1 / mu sub j): 2.06, 2.77, 2.67, 2.32, 1.88.

* Row 2 (kappa-hat sub j+1 / kappa-hat sub j): 0.57, 0.82, 0.96, 0.62, 0.70.

Section 2: Variance ratio test statistics.

The columns are labeled by the parameter d sub 0 from 0 to 5.

* Row 1 (Inv.VR): d sub 0 = 0 is 0.72 with a dagger symbol. d sub 0 = 1 is 0.07. d sub 0 = 2 is 0.04. d sub 0 = 3 is 0.01. d sub 0 = 4 is 0.01. d sub 0 = 5 is 0.01.

* Row 2 (VR 2,1): d sub 0 = 1 is 423.84**. d sub 0 = 2 is 10,347.35**. d sub 0 = 3 is 39,933.38**. d sub 0 = 4 is 101,141.36**. d sub 0 = 5 is 139,311.88**.

* Row 3 (VR 2,0): d sub 0 = 1 is 17,922.13**. d sub 0 = 2 is 481,545.37**. d sub 0 = 3 is 812,958.94**. d sub 0 = 4 is 2,993,476.52**. d sub 0 = 5 is 5,194,233.83**.

* Row 4 (VR 1,0): d sub 0 = 1 is 117.55**. d sub 0 = 2 is 558.39**. d sub 0 = 3 is 1,291.18**. d sub 0 = 4 is 1,695.46**. d sub 0 = 5 is 3,269.15**.

A footer note indicates that dagger, asterisk, and double asterisk indicate significance at the 10 percent, 5 percent, and 1 percent levels, respectively.

Note: $^{\dagger }$ , $^{*}$ , and $^{**}$ indicate significance at 10%, 5%, and 1% level, respectively.

Figure 4

Time series of $\langle \Delta X_t,\widehat {\nu }_j \rangle $ .

Note: Time series of the first three score functions of $\Delta X_t$ . The estimated eigenvectors $\widehat {\nu }_1,\widehat {\nu }_2$ , and $\widehat {\nu }_3$ are computed from $\widehat \Lambda _{1,R}$ .

Figure 4 Long description

The figure consists of three panels labeled a, b, and c, each displaying a highly volatile blue line graph against a light gray grid.

* Panel a is titled with the inner product of Delta X sub t and nu hat sub 1. The horizontal x-axis represents time, with major ticks at 1991 forward slash 01, 2001 forward slash 01, 2010 forward slash 01, and 2019 forward slash 01. The vertical y-axis ranges from negative 0.2 to 0.2. The data shows high-frequency oscillations around zero, with a notable positive peak around 2008.

* Panel b is titled with the inner product of Delta X sub t and nu hat sub 2. It shares the same x-axis timeline. The y-axis ranges from negative 0.1 to 0.1. The oscillations appear more uniform in amplitude compared to panel a, with frequent spikes above and below the zero baseline throughout the entire period.

* Panel c is titled with the inner product of Delta X sub t and nu hat sub 3. It shares the same x-axis timeline. The y-axis ranges from negative 0.2 to 0.2. This panel shows the highest density of fluctuations, with several sharp downward spikes reaching negative 0.2, particularly around the year 2000.

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 6. For the eigenvalue ratio estimators, we find $\ddot {\mathbb {d}} = \widehat {\mathbb {d}}_{\text { {LRS}}} = 1$ . However, the sequential testing results are quite different. Although the inverse VR test results are reported only for $\mathbb {d}_0=6,\ldots , 11$ , the null hypotheses $\mathbb {d}_0=0,1,\ldots , 5$ are all rejected at 1% level. Thus, we find that $\widehat {\mathbb {d}}_{\mathrm {BU}} = 10$ at 5% level and $\widehat {\mathbb {d}}_{\mathrm {BU}} = 9$ at 1% level, and we use this to set $\mathbb {d}_{\max } = 15$ . We report results for the VR tests for $\mathbb {d}_0=6,7,\ldots ,11$ . The larger values of $\mathbb {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 1. We then find $\widehat {\mathbb {d}}_{\mathrm {UD}} = 8$ based on VR(1,0) at 5% level, $\widehat {\mathbb {d}}_{\mathrm {UD}}=9$ based on VR(2,1) at 5% level or VR(1,0) at 1% level, and $\widehat {\mathbb {d}}_{\mathrm {UD}}=10$ based on VR(2,1) at 1% level. On balance, we conclude that there are likely nine 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 2.2).

Table 6

Empirical results for labor indices in levels

Table 6 Long description

The table is organized into two main sections under the header Eigenvalue ratio estimates.

Section 1: Eigenvalue ratio estimates.

Columns are labeled j equals 1 through j equals 6.

Row 1: Statistic mu sub j plus 1 all over mu sub j. Values are 2.06, 2.02, 1.74, 1.37, 1.34, and 1.42.

Row 2: Statistic kappa-hat sub j plus 1 all over kappa-hat sub j. Values are 0.07, 0.27, 0.34, 0.82, 0.80, and 0.50.

Section 2: Variance ratio test statistics.

Columns are labeled d sub 0 equals 6 through d sub 0 equals 11.

Row 1: Inv.V R. Values are 0.15 double asterisk, 0.11 double asterisk, 0.06 double asterisk, 0.03 asterisk, 0.03, and 0.02.

Row 2: V R (2,1). Values are 2,557.11; 4,531.11; 6,119.56; 10,730.08; 16,691.44 asterisk; and 22,714.94 double asterisk.

Row 3: V R (2,0). Values are 74,227.09; 157,057.93; 235,784.03; 738,700.29; 1,581,152.06; and 2,481,609.80 asterisk.

Row 4: V R (1,0). Values are 169.35; 229.45; 300.84; 564.11 asterisk; 762.54 double asterisk; and 898.93 double asterisk.

Note: Asterisk and double asterisk indicate significance at 5 percent and 1 percent levels respectively.

Note: $^{\dagger }$ , $^{*}$ , and $^{**}$ indicate significance at 10%, 5%, and 1% level, respectively.

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