Stationary subspace analysis of nonstationary covariance processes: eigenstructure description and testing

TL;DR

This paper develops eigenstructure-based methods to identify and test stationary subspaces within nonstationary covariance processes, providing a computationally efficient alternative to optimization techniques.

Contribution

It characterizes stationary subspaces via eigenvalues and eigenvectors, and introduces statistical tests and estimation methods for their dimensions.

Findings

Eigenstructure characterization of stationary subspaces.

Proposed eigenstructure-based estimation methods.

Validated methods on simulated and real data.

Abstract

Stationary subspace analysis (SSA) searches for linear combinations of the components of nonstationary vector time series that are stationary. These linear combinations and their number defne an associated stationary subspace and its dimension. SSA is studied here for zero mean nonstationary covariance processes. We characterize stationary subspaces and their dimensions in terms of eigenvalues and eigenvectors of certain symmetric matrices. This characterization is then used to derive formal statistical tests for estimating dimensions of stationary subspaces. Eigenstructure-based techniques are also proposed to estimate stationary subspaces, without relying on previously used computationally intensive optimization-based methods. Finally, the introduced methodologies are examined on simulated and real data.

Figures

Tables

Table 1. Table 1 : The 4 cluster centers, cluster sizes (in percentage), and the proportions of u 𝑢 u ’s for the 4 biggest clusters in Figure 2 . The selected (1,1)-local stationary subspace identified as the purple point in Figure 2 corresponds to ( − 0.114 , 0.778 , 0.618 ) 0.114 0.778 0.618 (-0.114,0.778,0.618) .

Cluster center	Cluster size (in %)	Proportion of $u$’s
$(0.661,0.750,-0.003)$	19.79	79.16
$(0.029,-0.641,0.767)$	23.95	91.67
$(-0.114,0.778,0.618)$	23.95	95.83
$(0.938,-0.346,-0.018)$	17.70	70.85

Table 2. Table 2 : Models 1-4: The canonical angle based measure D 3 ​ ( B ^ 1 ) subscript 𝐷 3 subscript ^ 𝐵 1 D_{3}(\widehat{B}_{1}) , between the population and estimated subspace. The empirical standard errors of the respective measures are provided in brackets.

	Model 1		Model 2		Model 3		Model 4
$T$	DSSA	VC	DSSA	VC	DSSA	VC	DSSA	VC
200	0.3277	0.3489	0.7548	0.4963	0.3633	0.3489	0.7453	0.6580
	(0.1429)	(0.1479)	(0.0993)	(0.2348)	(0.1693)	(0.1401)	(0.0914)	(0.1761)
500	0.2976	0.2878	0.6927	0.4145	0.3575	0.3301	0.6866	0.6196
	(0.1399)	(0.1288)	(0.1002)	(0.2078)	(0.1657)	(0.1329)	(0.1348)	(0.1931)
1000	0.2757	0.2968	0.6049	0.3393	0.3241	0.2649	0.6806	0.5655
	(0.1306)	(0.1301)	(0.1070)	(0.1953)	(0.1453)	(0.1228)	(0.1175)	(0.2004)
2000	0.2678	0.2374	0.5094	0.2751	0.2461	0.1931	0.5932	0.4447
	(0.1323)	(0.1311)	(0.1015)	(0.1905)	(0.1347)	(0.1202)	(0.1311)	(0.2261)

Table 3. Table 3 : Left: Out-of-sample classification accuracy (in %) for the 3 subjects S3, S5 and S8 for the two indicated methods with p = 5 𝑝 5 p=5 . Right: Out-of-sample classification accuracy (in %) for 3 subjects S3, S5 and S8 corresponding to d 𝑑 d = 7, 9, 11, and 13 for the VC method with p = 22 𝑝 22 p=22 .

$d$		S3	S5	S8
$1$	DSSA	56.25	49.03	51.11
	VC	50.69	52.08	46.15
$2$	DSSA	54.86	54.16	56.45
	VC	48.61	55.56	45.05
$3$	DSSA	59.02	59.33	64.39
	VC	47.22	57.63	54.54
$4$	DSSA	56.25	62.50	66.28
	VC	65.97	64.58	59.44

Equations

$$X_{t}=MU_{t}=\left[\begin{array}[]{c|c}M_{s}&M_{n}\end{array}\right]\;\begin{bmatrix}U_{t}^{s}\\ U_{t}^{n}\end{bmatrix},$$

$$X_{t}=A\Big{(}\frac{t}{T}\Big{)}Y_{t},\quad t=1,2,\ldots,T,$$

$$B_1^{\prime }{A}^2(u)B_1=\Sigma ,\forall u\in (0,1),$$

$$M(u)=A^2(u)-{\overline{A}}^2.$$

$$B_1^{\prime }M(u)B_1=0,\forall u\in (0,1).$$

$$C_1^{\prime }M{C}_1=0,$$

$$d_0=d_0(M),d_{+}=d_{+}(M),d_{-}=d_{-}(M)$$

$${\lambda }_i\le {\mu }_i\le {\lambda }_{n-d_1+i},i=1,\dots ,d_1,$$

$${\lambda }_i\le 0\le {\lambda }_{n-d_1+i},i=1,\dots ,d_1.$$

$$S^{\prime }MS=\text{mbox}diag({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p),$$

$$s_i^{\prime }{s}_j={\delta }_{ij},s_i^{\prime }M{s}_j={\delta }_{ij}{\lambda }_i,$$

$$s_{0,i}^{\prime }M{s}_{0,i}=0,i=1,\dots ,d_0.$$

$$s_{+,i}^{\prime }M{s}_{+,i}$$

$$s_{-,i}^{\prime }M{s}_{-,i}$$

$$({\alpha }_i^{1/2}{s}_{+,i}+{\beta }_i^{1/2}{s}_{-,i}{)}^{\prime }M({\alpha }_i^{1/2}{s}_{+,i}+{\beta }_i^{1/2}{s}_{-,i})={\alpha }_i(s_{+,i}^{\prime }M{s}_{+,i})+{\beta }_i(s_{-,i}^{\prime }M{s}_{-,i})=0,$$

$$s_{\pm ,i}=\frac{{\alpha }_i^{1/2}{s}_{+,i}+{\beta }_i^{1/2}{s}_{-,i}}{||{\alpha }_i^{1/2}{s}_{+,i}+{\beta }_i^{1/2}{s}_{-,i}|{|}_2},i=1,\dots ,d_{\pm },$$

$$C_1=(s_{0,1}\dots s_{0,d_0}{s}_{\pm ,1}\dots s_{\pm ,d_{\pm }}).$$

$$w=w_0+w_{\pm },$$

$$w_{\pm }=\sum _{i=1}^{d_{+}}{\alpha }_{+,i}{s}_{+,i}+\sum _{i=1}^{d_{-}}{\alpha }_{-,i}{s}_{-,i},$$

$$\sum _{i=1}^{d_{+}}{\alpha }_{+,i}^2{\lambda }_{+,i}+\sum _{i=1}^{d_{-}}{\alpha }_{-,i}^2{\lambda }_{-,i}=0.$$

$$0=w^{\prime }Mw=\sum _{i=1}^{d_0}{\alpha }_{0,i}^2{s}_{0,i}^{\prime }M{s}_{0,i}+\sum _{i=1}^{d_{+}}{\alpha }_{+,i}^2{s}_{+,i}^{\prime }M{s}_{+,i}+\sum _{i=1}^{d_{-}}{\alpha }_{-,i}^2{s}_{-,i}^{\prime }M{s}_{-,i}=\sum _{i=1}^{d_{+}}{\alpha }_{+,i}^2{\lambda }_{+,i}+\sum _{i=1}^{d_{-}}{\alpha }_{-,i}^2{\lambda }_{-,i}.$$

$$\sum _{i=1}^{d_{+}}{\alpha }_{1,+,i}{\alpha }_{2,+,i}{\lambda }_{+,i}+\sum _{i=1}^{d_{-}}{\alpha }_{1,-,i}{\alpha }_{2,-,i}{\lambda }_{-,i}=0.$$

$$\mathcal{P}(M)={\mathcal{N}}_0(M)\oplus {\mathcal{N}}_{\pm }(M),$$

$$d=d(u)=d(M(u)),{\mathcal{B}}_1={\mathcal{B}}_1(u)=\mathcal{P}(M(u)),\forall u\in (0,1).$$

$$d=d(u)=d(M(u))=d_0(u)+\min (d_{+}(u),d_{-}(u)),\forall u\in (0,1),$$

$$d_0(u)=d_0(M(u)),d_{+}(u)=d_{+}(M(u)),d_{-}(u)=d_{-}(M(u))$$

$$d(M(u))=d_0(u)+\min (d_{+}(u),d_{-}(u))\equiv d^{\ast },u\in (0,1),$$

$${\mathcal{B}}_1(u)=\mathcal{P}(M(u))\equiv {\mathcal{B}}_1^{\ast },u\in (0,1),$$

$$A^2(u)=\text{mbox}diag(2+\sin (2\pi u),3-\sin (2\pi u),1+\sin (2\pi u)).$$

$$M(u)=A^2(u)-{\overline{A}}^2=\sin (2\pi u)\cdot \text{mbox}diag(1,-1,1).$$

$${\mathcal{B}}_1(u)=\mathcal{P}(M(u))=\text{mbox}lin\{{\alpha }_{-,1}{s}_{-,1}+{\alpha }_{+,1}{s}_{+,1}+{\alpha }_{+,2}{s}_{+,2}\}$$

Full text

1 Introduction

The goal of this work is to provide new fundamental insights into the so-called stationary subspace analysis (SSA), a technique for finding linear combinations of components of a multivariate time series that are stationary. More precisely, consider a setup where the observed $p$-vector nonstationary time series $X_{t}$ is a linear transformation of a $d$-vector stationary series $U_{t}^{s}$ and a $(p-d)$-vector nonstationary series $U_{t}^{n}$ through

[TABLE]

and $M$ is an unknown $p\times p$ (invertible) mixing matrix, $M_{s}$ and $M_{n}$ are $p\times d$ and $p\times(p-d)$ matrices, respectively. It is further assumed that no linear transformation of $U_{t}^{n}$ is stationary. Given the data $X_{1},X_{2},\ldots,X_{T}$, SSA seeks to find the demixing matrix $B=(M^{-1})^{\prime}$ so that $B^{\prime}X_{t}=U_{t}$ is naturally partitioned into its stationary and nonstationary sources. The space spanned by the first $d$ columns of $B$ is referred to as a stationary subspace and $d$ as its dimension.

SSA was introduced and studied by von Bünau et al. (2009), with applications to analyzing EEG data from neuroscience experiments. In that work, the observed vector time series is assumed to be independent across time and the notion of stationarity is with respect to the first two moments, that is, the mean and lag-0 covariance are required to be time invariant. The demixing matrix in SSA is found in the spirit of ANOVA by dividing the observed time series data into $N$ segments and minimizing a Kullback-Leibler (KL) divergence between Gaussian distributions measuring differences in the means and covariances across these segments. A sequential likelihood ratio test is used in von Bünau et al. (2009) and Blythe et al. (2012) to determine the dimension of the stationary subspace $d$ under the additional assumption of normality of the data. The frequency domain or dependent SSA (DSSA) in Sundararajan and Pourahmadi (2018) avoids dividing the data into segments and relies on the approximate uncorrelatedness of the discrete Fourier transform (DFT) of a second-order stationary time series at Fourier frequencies. The sum of the Frobenius norms of the estimated covariances of the DFTs at the first few lags is used as a discrepancy measure and the demixing matrix is obtained by optimizing this measure. Then, a sequential test of second-order stationarity is used to determine $d$ and the consistency of the estimated $d$ is discussed using the asymptotic distribution of the test statistic under the alternative hypothesis of local stationarity of the time series (Dahlhaus (1997, 2012)).

Overall, the research thus far suggests the need for a better mathematical formulation and understanding of the problem. A reformulation and, in particular, an optimization-free solution seems to be the key in finding a transparent and interpretable solution of the SSA problem. In this work, we shall tackle these issues for a special but general case of (1.1), namely, that of zero mean vectors $U_{t}$, assuming $U_{t}=C(t)Y_{t}$ with i.i.d. zero mean vectors $Y_{t}$. We shall further write this model formulation as

[TABLE]

where $T$ is the sample size and the time dependence is brought into $A(\frac{t}{T})=M\cdot C(t)$. Additional assumptions on the matrix-valued function $A:(0,1)\rightarrow\mathbb{R}$ and the i.i.d. vectors $Y_{t}$ can be found below in Section 2. The modification (1.2) of (1.1) takes the heterogeneity out of $U_{t}$ and places it into the deterministic matrix-valued function $A(\cdot)$. The nonstationary covariance process (1.2) will be said to follow a varying covariance (VC) model. The assumption of zero mean in (1.2) is made for several reasons. For one, all previous works that identify EEG data analysis as the motivating application involve this assumption. We are currently working in parallel on analogous approaches to the SSA problem for time-varying means (Düker et al. (2019)), and will possibly look at the combined model in the future. In the latter regard, we should also note that dealing with SSA for varying covariances is seemingly much more involved than for varying means. Indeed, as seen from this work, the SSA for the VC model has a surprisingly rich structure.

Our contributions to the SSA for the VC model (1.2) are as follows. First, by using basic ideas from linear algebra, we provide an interpretation of a stationary subspace and its dimension $d$ in terms of eigenvalues and eigenvectors of certain symmetric matrices. This interpretation, in fact, is given at two levels: “local” or for fixed $u\in(0,1)$, and “global” or for $(0,1)=\{u:\;u\in(0,1)\}$, where $u$ is thought here as a variable of $A(u)$ replaced by $t/T$ in (1.2). Second, in the context of the obtained interpretation, we develop formal statistical tests for both “local” and “global” dimensions of stationary subspaces. Together with the algebraic interpretation, these tests are the key theoretical contributions of this work. Third, by leveraging the new interpretation of stationary subspaces, we provide more direct and algebraic ways to construct them. These are shown to outperform the computationally more expensive optimization-based solutions of the previous SSA approaches in a number of simulation settings. We should also note that the proposed dimension tests in Section 4 assume the existence of a common stationary subspace and dimension across the “local” levels (see Section 3.3 for more details); testing for the latter remains an open problem. Furthermore, this work concerns the asymptotics under $T\rightarrow\infty$ with a fixed $p$. Fourth, we revisit an SSA application to EEG data from a Brain-Computer Interface (BCI) experiment and provide additional insights by using the proposed methods.

The outline of the paper is as follows. Section 2 reintroduces more formally the VC model and its stationary subspace and dimension. Section 3 gives an eigenstructure-based characterization of a stationary subspace and its dimension. Section 4 considers statistical tests for the dimension of a stationary subspace at both the “local” and “global” levels. Section 5 discusses estimation methods for stationary subspaces that are based on algebraic constructs and do not involve iterative and computationally heavy optimization methods. Sections 6 and 7 illustrate the proposed methodologies using simulated and real data. Section 8 concludes.

2 Model of interest and its stationary subspace

We focus throughout this work on the varying covariance (VC) model (1.2), where $A:(0,1)\rightarrow\mathbb{R}^{p\times p}$ is a smoothly varying matrix-valued function and $Y_{t}$ are i.i.d. random vectors with i.i.d. entries, $\mathbb{E}(Y_{t})=0$ and $\mathbb{E}(Y_{t}Y_{t}^{\prime})=I_{p}$. Further technical assumptions can be found below.

Definition 2.1**.**

If $d$ is the largest integer in $\{0,1,2,\ldots,p\}$ for which there is a $p\times d$ matrix $B_{1}$ such that

[TABLE]

where $A^{2}(u)=A(u)A(u)^{\prime}$ , $\Sigma$ does not depend on $u$ and $B_{1}^{\prime}B_{1}=I_{d}$, then the space ${\cal B}_{1}$ spanned by the columns of $B_{1}$ will be called a (second-order) stationary subspace of dimension $d$ of the model (1.2).

Note that (2.1) states effectively that the covariance matrix of $B_{1}^{\prime}X_{t}$ does not depend on $t$. It can be reformulated as follows: Let $\overline{A}^{2}=\int_{0}^{1}A^{2}(u)du$ and define a $p\times p$ symmetric matrix $M(u)$ as

[TABLE]

Then, the condition (2.1) is equivalent to

[TABLE]

Indeed, (2.1) implies (2.3) after integrating (2.1) over $u\in(0,1)$, noting that $\int_{0}^{1}\Sigma du=\Sigma$ and subtracting the two sides of the resulting relation $B_{1}^{\prime}\overline{A}^{2}B_{1}=\Sigma$ from those of (2.1). Similarly (2.3) implies (2.1) with $\Sigma=B_{1}^{\prime}\overline{A}^{2}B_{1}$. The matrix $M(u)$ will play a central role henceforth.

Our approach to finding a matrix $B_{1}$ and the corresponding stationary subspace $\mathcal{B}_{1}$ and dimension $d$ is based on the relation (2.3) for a fixed $u$, that is, $B_{1}^{\prime}M(u)B_{1}=0$ for a fixed $u$ and $B_{1}=B_{1}(u)$ of dimension $d=d(u)$. As shown in the next section, for a fixed $u$, the matrix $B_{1}(u)$ and its dimension $d(u)$ can be characterized using the eigenstructure of the matrix $M(u)$. When it comes to $M(u)$, we shall be using the terminology of the following definition.

Definition 2.2**.**

Let $B_{1}=B_{1}(u)$ be a matrix with $d=d(u)$ columns that satisfies (2.3) for a fixed $u\in(0,1)$. The space ${\cal B}_{1}(u)$ spanned by the columns of $B_{1}(u)$ will be called a local stationary subspace of local dimension $d(u)$. The respective quantities in Definition 2.1 will be referred to as a global stationary subspace and a global dimension.

Relationships between local and global stationary subspaces and their dimensions are discussed in Section 3.3.

3 Matrix pseudo nullity and pseudo null space

In this section we describe an eigenstructure-based characterization of a stationary subspace and its dimension. We first define the notions of pseudo null space and psuedo nullity for symmetric matrices and then later, in Section 3.3, associate them to stationary subspaces and their dimensions. In view of the relation (2.3) and Definition 2.2, we start with the following definition.

Definition 3.1**.**

Let $M$ be any $p\times p$ symmetric matrix. A pseudo nullity of $M$, denoted by $d(M)$, is defined as the largest positive number $d_{1}$ such that

[TABLE]

for a $p\times d_{1}$ matrix $C_{1}$ with $C_{1}^{\prime}C_{1}=I_{d_{1}}$. A pseudo null space of $M$, denoted as ${\cal P}(M)$, is defined as the linear span of the $d_{1}$ columns of the matrix $C_{1}$ in (3.1). A column of $C_{1}$, that is, a column vector $s$ such that $s^{\prime}Ms=0$ will be called a pseudo eigenvector.

If $M$ is positive semi-definite, note that its pseudo nullity is its nullity (i.e. the number of zero eigenvalues of $M$) and its pseudo null space is its null space; thus, the psuedo- terminology is used to draw attention to the contrast between these two cases. Otherwise, we should caution the reader against drawing other parallels between the two contexts. For example, if $s_{1}$ and $s_{2}$ are two pseudo eigenvectors (which can be either orthogonal or non-orthogonal), note that there is a priori no reason to have $s_{1}^{\prime}Ms_{2}=0$ and hence $C_{1}^{\prime}MC_{1}=0$ with $C_{1}=(s_{1}\,s_{2})$. In particular, if for e.g. $d(M)=2$, and $s_{1}$ and $s_{2}$ are orthogonal, the linear space spanned by $s_{1}$ and $s_{2}$ is not necessarily a pseudo null space.

Another word of caution is that ${\cal P}(M)$ is not unique in general. This, perhaps surprising, fact will be explained below. By writing ${\cal P}(M)$, we mean one of the pseudo null spaces.

3.1 Characterization of pseudo nullity

We characterize here the pseudo nullity $d(M)$ of a symmetric matrix $M$ in terms of its inertia. More precisely, let

[TABLE]

be the number of zero, positive and negative eigenvalues of $M$, respectively. Part of the proof of the following result is constructive and will be used to construct pseudo null spaces (stationary subspaces) in Section 5.

Proposition 3.1**.**

*Let $M$ be a symmetric matrix. Then, $d(M)=d_{0}+\min(d_{+},d_{-})$. *

Proof.

By the Poincare Separation Theorem (e.g. Magnus and Neudecker (1999)),

[TABLE]

where $\mu_{1}\leq\ldots\leq\mu_{d_{1}}$ are the ordered eigenvalues of $C_{1}^{\prime}MC_{1}$ for any $p\times d_{1}$ matrix $C_{1}$ such that $C_{1}^{\prime}C_{1}=I_{d_{1}}$. Taking $C_{1}$ as in (3.1) with $d_{1}=d(M)$, we have $C_{1}^{\prime}MC_{1}=0$, and hence

[TABLE]

This shows that $d(M)\leq d_{0}+\min(n_{+},n_{-})$.

To prove the reverse inequality $d_{0}+\min(n_{+},n_{-})\leq d(M)$, let $d_{*}=d_{0}+\min(d_{+},d_{-})$. The rest of the proof constructs a $p\times d_{*}$ matrix $C_{1}$ such that $C_{1}^{\prime}MC_{1}=0$ and $C_{1}^{\prime}C_{1}=I_{d_{*}}$, which yields the desired inequality. By the Schur decomposition (e.g. Magnus and Neudecker (1999)),

[TABLE]

where the columns $s_{1},\ldots,s_{p}$ of $S$ are orthonormal and represent the eigenvectors associated with the negative, positive and zero eigenvalues of $M$. We now need to separate the eigenvectors into those associated with the eigenvalues $\lambda_{1},\ldots,\lambda_{p}$ of $M$. Let $s_{0,i}$ be the eigenvectors associated with the zero eigenvalues $\lambda_{0,i}$, $i=1,\ldots,d_{0}$, $s_{+,i}$ be the eigenvectors associated with the positive eigenvalues $\lambda_{+,i}$, $i=1,\ldots,d_{+}$, and $s_{-,i}$ be the eigenvectors associated with the negative eigenvalues $\lambda_{-,i}$, $i=1,\ldots,d_{-}$. Note by (3.3) that

[TABLE]

where $\delta_{ij}=1$ if $i=j$ and [math] otherwise.

For the zero eigenvalues and the corresponding eigenvectors, we have $Ms_{0,i}=0\cdot s_{0,i}$ and hence

[TABLE]

Similarly, for the positive and negative eigenvalues, we have

[TABLE]

Let $d_{\pm}=\min(d_{+},d_{-})$. Then, for some $\alpha_{i}\in(0,1)$ and $\beta_{i}=1-\alpha_{i}$, we have

[TABLE]

for $i=1,\ldots,d_{\pm}$. Set

[TABLE]

and define a $(d_{0}+d_{\pm})\times p=d_{*}\times p$ matrix $C_{1}$ as

[TABLE]

By using (3.4), (3.5) and (3.6), we have $C_{1}^{\prime}MC_{1}=0$. Since $s_{i}$ are orthonormal and in view of (3.7), we also have $C_{1}^{\prime}C_{1}=I_{d_{*}}$. This shows that $d_{0}+d_{\pm}\leq d(M)$ and concludes the proof. ∎

3.2 Characterization of pseudo null space

A closer examination of the proof of Proposition 3.1 also suggests further insights into the structure of a pseudo null space ${\cal P}(M)$. The next auxiliary result describes an element of ${\cal P}(M)$. As in the proof of Proposition 3.1, we let $s_{i}$ ($s_{0,i}$, $s_{+,i}$ and $s_{-,i}$, resp.) be the orthonormal eigenvectors (associated with the zero, positive and negative eigenvalues, resp.) of $M$. The corresponding eigenvalues are denoted $\lambda_{i}$ ($\lambda_{0,i}=0$, $\lambda_{+,i}$, $\lambda_{-,i}$ resp.). We also let ${\cal N}_{0}(M)$ denote the null space of the matrix $M$, that is, the linear space spanned by the eigenvectors $s_{0,i}$, $i=1,\ldots,d_{0}$.

Lemma 3.1**.**

Let $M$ be a symmetric matrix and $w$ be a pseudo eigenvector. Then,

[TABLE]

where $w_{0}\in{\cal N}_{0}(M)$ and

[TABLE]

where $\alpha_{-,i},\alpha_{+,i}\in\mathbb{R}$ are such that

[TABLE]

Conversely, a vector $w$ expressed through (3.8)–(3.10) is a pseudo eigenvector.

Proof.

Any vector $w$ can be expressed as a linear combination of the basis vectors $s_{i}$ as in (3.8)–(3.9). The relation (3.10) follows since

[TABLE]

The converse statement follows similarly. ∎

The next lemma clarifies when two pseudo eigenvectors belong to a pseudo null space. See also the discussion following Definition 3.1.

Lemma 3.2**.**

Let $M$ be a symmetric matrix and $w_{k}$, $k=1,2$, be its two pseudo eigenvectors expressed through (3.9)–(3.10) with the coefficients $\alpha_{k,+,i}$ and $\alpha_{k,-,i}$, $k=1,2$. Then, $w_{k}$, $k=1,2$, belong to the same pseudo null space if and only if

[TABLE]

Proof.

The result follows as in the proof of Lemma 3.1 above from requiring that $w_{1}^{\prime}Mw_{2}=0$. ∎

The next result characterizes a pseudo null space ${\cal P}(M)$.

Proposition 3.2**.**

Let $M$ be a symmetric matrix and ${\cal P}(M)$ be its pseudo null space. Then,

[TABLE]

where ${\cal N}_{\pm}(M)$ is a linear space spanned by $\min(d_{-},d_{+})$ orthogonal pseudo eigenvectors expressed as (3.9)–(3.10) and satisfying (3.11) pairwise. Conversely, the right-hand side of (3.12) defines a pseudo null space. Moreover, a pseudo null space is not unique in general.

Proof.

The first two statements follow from Lemmas 3.1–3.2 and Proposition 3.1. The last statement is illustrated in Example 3.1 given below. ∎

3.3 Implications for stationary subspace and its dimension

In view of Definitions 2.1 and 2.2 and their notation, the global stationary subspace $\mathcal{B}_{1}$ and its dimension $d$ are given by:

[TABLE]

In view of Proposition 3.1, the first relation in (3.13) can now be reformulated as

[TABLE]

where

[TABLE]

are the numbers of zero, positive and negative eigenvalues of $M(u)$, respectively.

On the other hand, if

[TABLE]

for some $d^{*}$, this does not necessarily mean that the dimension $d$ of the stationary subspace is $d^{*}$. The latter is because (3.16) does not guarantee that

[TABLE]

for some ${\cal B}_{1}^{*}$, playing the role of a stationary subspace. But if (3.17) holds, then (3.16) does imply that $d^{*}$ is the dimension $d$ of the stationary subspace. The statistical tests developed in the subsequent sections will, in fact, be for testing (3.16) and hence will lead to the dimension $d$ assuming (3.17). How testing can be done for (3.17) remains an open question, though we shall also suggest new ways to estimate ${\cal B}_{1}$, based on developments in this section.

Finally, we illustrate the observations made above through a simple but instructive example. See also a subsequent remark.

Example 3.1**.**

Consider a VC model with

[TABLE]

Then, $\overline{A}^{2}=\mbox{diag}(2,3,1)$ and

[TABLE]

For fixed $u\neq 1/2$, the eigenvalues of $M(u)$ are $\sin(2\pi u)\cdot 1$ (of multiplicity 2) and $\sin(2\pi u)\cdot(-1)$. For further illustration, suppose $u\in(0,1/2)$, so that $\sin(2\pi u)>0$. Then, $\lambda_{+,1}=\lambda_{+,2}=\sin(2\pi u)\cdot 1$, $\lambda_{-,1}=\sin(2\pi u)\cdot(-1)$ and $d_{+}=2$, $d_{-}=1$, $d_{0}=0$, by using the notation in Section 3.2 with $M=M(u)$. By Proposition 3.1, $d(u)=d(M(u))=0+\min(1,2)=1$. The corresponding eigenvectors are $s_{+,1}=(1\ 0\ 0)^{\prime}$, $s_{+,2}=(0\ 0\ 1)^{\prime}$ and $s_{-,1}=(0\ 1\ 0)^{\prime}$. By Proposition 3.2, a local stationary subspace or a pseudo null space of $M(u)$ can be expressed as

[TABLE]

such that

[TABLE]

and $\alpha_{-,1}^{2}+\alpha_{+,1}^{2}+\alpha_{+,2}^{2}=1$ for the norm to be $1$, where “lin” indicates a linear span. The latter two expressions yield $\alpha_{-,1}=(\alpha_{+,1}^{2}+\alpha_{+,2}^{2})^{1/2}$ (after choosing a positive sign for the square root) and $\alpha_{+,1}^{2}+\alpha_{+,2}^{2}=1/2$. This further yields $\alpha:=\alpha_{+,1}\in[1/\sqrt{2},-1/\sqrt{2}]$, $\alpha_{+,2}=\pm(1/2-\alpha^{2})^{1/2}$ and $\alpha_{-,1}=1/\sqrt{2}$. Thus, a pseudo null space can also be expressed as

[TABLE]

where $\alpha\in[1/\sqrt{2},-1/\sqrt{2}]$. Note that these spaces (vectors) are generally different for different $\alpha$’s. For example, for $\alpha=0$,

[TABLE]

and for $\alpha=1/\sqrt{2}$,

[TABLE]

Remark 3.1**.**

The fact that a pseudo null space in Example 3.1 is not unique should not be surprising from the following perspective. Let $X_{t}=(X_{1,t},X_{2,t},X_{3,t})^{\prime}$ be a $3$-vector process following the VC model with (3.18). The pseudo eigenvectors $w_{1}=(0\ 1\ 1)^{\prime}$ in (3.21) and $w_{2}=(1\ 1\ 0)^{\prime}$ in (3.22) can be checked easily to be such that $w_{1}^{\prime}X_{t}$ and $w_{2}^{\prime}X_{t}$ are stationary.

It is also interesting and important to note here that the stationary dimension for this model is not $2$. For example, note that while $w_{1}^{\prime}X_{t}$ and $w_{2}^{\prime}X_{t}$ are indeed stationary, the $2$-vector process $(w_{1}^{\prime}X_{t},w_{2}^{\prime}X_{t})^{\prime}$ is not stationary. Indeed, this is the case since e.g. $\mathbb{E}(w_{1}^{\prime}X_{t})(w_{2}^{\prime}X_{t})=\mathbb{E}(X_{2,t}+X_{3,t})(X_{1,t}+X_{2,t})=\mathbb{E}X_{2,t}^{2}=3-\sin(2\pi t/T)$.

4 Inference of stationary subspace dimension

Here we discuss the statistical tests for the dimension of a stationary subspace at both the “local” and “global” levels along with their asymptotic properties.

According to (3.14)–(3.15), the dimension of a stationary subspace of the VC model is the local dimension $d(u)$ or pseudo nullity $d(M(u))$ of the matrix $M(u)$ in (2.2), assuming it is the same across $u$, which can further be expressed in terms of $d_{0}(u)$, $d_{+}(u)$ and $d_{-}(u)$. We are interested here in the statistical testing for $d_{0}(u)$, $d_{+}(u)$, $d_{-}(u)$ and hence also for $d(u)=d(M(u))$. We consider below two types of tests: local (that is, for fixed $u$) and global (that is, for an interval of $u$). Though it is the global test which is most relevant for (3.14), we consider the local test since its statistic forms the basis of the global test and also because of independent interest.

For statistical inference, we obviously need an estimator of the matrix $M(u)$. It is set naturally as

[TABLE]

where $K_{h}(u)=h^{-1}K(h^{-1}u)$, $K(\cdot)$ is a kernel function and $h$ denotes the bandwidth. A kernel is a symmetric function which integrates to $1$, with further regularity assumptions possibly made as well. In simulations and data application, we work with the triangle kernel $K(u)=1-|u|$ if $|u|<1$ and [math], otherwise.

4.1 Local dimension test

In this section, $u$ is assumed to be fixed. A pseudo nullity $d(M(u))$ of a symmetric matrix $M(u)$ could be tested by adapting the matrix rank tests found in Donald et al. (2007). For this, we need an asymptotic normality result for the estimator $\widehat{M}(u)$, which is stated next. The proof can be found in Appendix A.1.

Proposition 4.1**.**

Suppose that the assumptions for the proposition stated in Appendix A.1 hold, in particular, $T\to\infty$, $h\to 0$ so that $Th\to\infty$ and $Th^{3}\to 0$. Then, we have

[TABLE]

where $\|K\|_{2}^{2}=\int_{\mathbb{R}}K(v)^{2}dv$, $\mu_{4}=\mathbb{E}(Y_{i,t}^{2}-1)^{2}=\mathbb{E}Y_{i,t}^{4}-1$ and ${\cal Z}_{p}$ is a symmetric $p\times p$ matrix having independent normal entries with variance $1$ on the diagonal and variance $1/\mu_{4}$ off the diagonal. Moreover, $\widehat{A}^{k}(u)\to_{p}A^{k}(u)$, $k=\pm 1,\pm 2$.

Note that the asymptotic normality result (4.2) can be expressed as

[TABLE]

where

[TABLE]

Furthermore, by Proposition 4.1, we have

[TABLE]

with $\widehat{F}(u)=\widehat{A}(u)^{-1}$. In principle, $\mu_{4}$ in $a_{T}$ needs to be estimated as well. But for simplicity and better readability, we shall assume that $\mu_{4}$ is known. Construction of consistent estimators $\widehat{\mu}_{4}$ is discussed in the supplementary technical appendix (Sundararajan et al. (2019)), and would be sufficient for the results below to hold assuming that $\mu_{4}$ is estimated through $\widehat{\mu}_{4}$. Furthermore, if one is willing to assume Gaussianity of $Y_{t}$, note that $\mu_{4}=2$ can also be used.

The convergence (4.3) with (4.5) is the setting for the matrix rank tests considered in Donald et al. (2007), with one small difference that plays little role as noted in the proof of Corollary 4.1 below. We explain below how the results of Donald et al. (2007) can be adapted to test for $d_{0}(u)$, $d_{+}(u)$, $d_{-}(u)$ and hence also for $d(M(u))$. But first it is convenient to gather some of the notation to be used on a number of occasions below. In view of (4.3), note that the matrix $F(u)=A(u)^{-1}$ plays the role of standardization. In this sense, the focus should not be so much on $M(u)$ but rather on $F(u)M(u)F(u)^{\prime}$. Note that $d_{0}(u)=d_{0}(M(u)))=d_{0}(F(u)M(u)F(u)^{\prime})$ and also $d_{\pm}(M(u))=d_{\pm}(F(u)M(u)F(u)^{\prime})$. We let

[TABLE]

be the ordered eigenvalues of $F(u)M(u)F(u)^{\prime}$ and $\widehat{F}(u)\widehat{M}(u)\widehat{F}(u)^{\prime}$, respectively. We have

[TABLE]

Let also

[TABLE]

be the ordered eigenvalues of $(F(u)M(u)F(u)^{\prime})^{2}$ and $(\widehat{F}(u)\widehat{M}(u)\widehat{F}(u)^{\prime})^{2}$, respectively. We have $(\widehat{\gamma}_{i}(u))^{2}=\widehat{\gamma}_{2,j(i)}(u)$ for some $j(i)$, and a similar expression with the hats and also

[TABLE]

In particular, $d_{0}(u)=d_{0}(F(u)M(u)F(u)^{\prime})=d_{0}((F(u)M(u)F(u)^{\prime})^{2})$. By combining these observations, we have

[TABLE]

In the proofs for Section 4.2, we shall also use eigenspaces associated with the eigenvalues above but these will not be discussed here.

We first look at the inference about $d_{0}(u)$. Since $d_{0}(u)$ is equal to $p-\mbox{rk}(M(u))$, where $\mbox{rk}(M(u))$ is the rank of $M(u)$, we can test directly for $d_{0}(u)$ by using the MINCHI2 test of Donald et al. (2007). By the discussion above, $\mbox{rk}\{M(u)\}=\mbox{rk}\{(F(u)M(u)F(u)^{\prime})^{2}\}$. It is then natural to consider for $r=0,\ldots,p$, the test statistic

[TABLE]

Its asymptotics is described in the following result.

Corollary 4.1**.**

Suppose that the assumptions of Proposition 4.1 hold. Then, under $H_{0}:d_{0}(u)=r$,

[TABLE]

and under $H_{1}:d_{0}(u)<r,\;\widehat{\xi}_{r}(u)\to_{p}+\infty$.

Proof.

The result follows from Theorem 4.7 in Donald et al. (2007). Note that the variance in the off-diagonal of the matrix ${\cal Z}_{p}$ in (4.3) plays no role in the derivation of the theorem, since the argument is based on a trace and hence only on the variance on the diagonal. ∎

The corollary can be used to test for $d_{0}(u)$ in a standard way sequentially, namely, testing for $H_{0}:d_{0}(u)=r$ starting with $r=p$ and subsequently decreasing $r$ by $1$ till the null hypothesis is not rejected. Let $\widehat{d}_{0}(u)$ be the resulting estimator, which is consistent for $d_{0}(u)$ under a suitable choice of critical values in the sequential testing as the following result states.

Corollary 4.2**.**

Suppose that the assumptions of Proposition 4.1 hold. Let $\widehat{d}_{0}(u)$ be the estimator of $d_{0}(u)$ defined above through the sequential testing procedure when using a significance level $\alpha=\alpha_{T}$ in testing. If $\alpha_{T}\to 0$ and $(-\log\alpha_{T})/Th\to 0$, then $\widehat{d}_{0}(u)\to_{p}d_{0}(u)$.

Proof.

The result can be proved as e.g. Theorem 4.3 in Fortuna (2008). ∎

We now turn to inference of $d_{+}(u)$ and $d_{-}(u)$. By the discussion surrounding equations (4.6)–(4.9), these quantities are also $d_{\pm}(F(u)M(u)F(u)^{\prime})$. In view of the relation (4.10), since $\widehat{d}_{0}(u)$ estimates $d_{0}(u)$ consistently and the sample eigenvalues converge to their population analogues, it is natural to set

[TABLE]

Corollary 4.3**.**

Under the assumptions of Corollary 4.2, we have $\widehat{d}_{+}(u)\to_{p}d_{+}(u)$ and $\widehat{d}_{-}(u)\to_{p}d_{-}(u)$.

Proof.

The stated result follows from the following two observations. First, by Corollary 4.2, $\widehat{d}_{0}(u)=d_{0}(u)$ a.s. for large enough $T$ (depending on $\omega$ in “a.s.”). Second, since $(\widehat{F}(u)\widehat{M}(u)\widehat{F}(u)^{\prime})^{k}\to_{p}(F(u)M(u)F(u)^{\prime})^{k}$ for $k=1,2$, we have $\widehat{\gamma}_{2,j}(u)\to_{p}\gamma_{2,j}(u)$, $j=1,\ldots,p$, and $\widehat{\gamma}_{i}(u)\to_{p}\gamma_{i}(u)$, $i=1,\ldots,p$. Thus, in the limit, the right-hand side of (4.13) becomes that of (4.10), and hence also its left-hand side, that is, $d_{\pm}(u)$. ∎

Finally, a natural estimator for $d(u)=d(M(u))$ is

[TABLE]

Corollaries 4.2–4.3 imply the consistency of this estimator, which is the main result of this section.

Theorem 4.1**.**

Under the assumptions of Corollary 4.2, we have $\widehat{d}(M(u))\to_{p}d(M(u))$.

4.2 Global dimension test

In the previous section, we considered testing about $d_{0}(u)$, $d_{+}(u)$, $d_{-}(u)$ and $d(M(u))$ for fixed $u$. The relation (3.14) of interest, however, concerns these quantities for all $u\in(0,1)$. We would thus like to be able to make inference “globally,” that is, for an interval of $u$. Since our approach to $d(M(u))$ goes through $d_{0}(u)$, $d_{+}(u)$ and $d_{-}(u)$, we shall make inference about these quantities first. To deal with the possibility that these quantities might differ across $u\in(0,1)$, we shall work under the assumption that

[TABLE]

and develop a global dimension test about $d_{0}$, $d_{+}$ and $d_{-}$ in (4.15) for fixed ${\cal H}$. In practice, we shall apply the developed test over refined dyadic partitions, first for $(0,1)$, then for $(0,1/2)$ and $(1/2,1)$, then for $(0,1/4),(1/4,1/2),(1/2,3/4)$ and $(3/4,1)$ etc. What to expect under this splitting of the interval $(0,1)$ is discussed in Section 4.3, and will be illustrated in Sections 6 and 7.

We first focus on inference about $d_{0}$ in (4.15). Our global test will be based on the quantity $\int_{\cal H}\widehat{\xi}_{r}(u)du$, where $\widehat{\xi}_{r}(u)$ is the statistic (4.11) used in the local dimension test. Its asymptotics are described in the next result, which also defines the global statistic $\widehat{\xi}_{r}$. Recall the notation $\mu_{4}=\mathbb{E}(Y_{i,t}-1)^{4}$ used in Proposition 4.1; see also the discussion following (4.5). We also let $|{\cal H}|$ be the length of the interval ${\cal H}$ and $\overline{K}(u)=\int_{\mathbb{R}}K(v)K(u-v)dv$ be the so-called convolution kernel.

Proposition 4.2**.**

Suppose that the assumptions for the proposition stated in Appendix A.2 hold. Then, under $H_{0}:d_{0}(u)\equiv r$ for all $u\in{\cal H}$,

[TABLE]

and under $H_{1}:d_{0}(u)<r$ for some $u\in{\cal H},\;\widehat{\xi}_{r}\to_{p}\infty$.

The proof of the proposition can be found in Appendix A.2, and follows the approach taken in Donald et al. (2011).

As for the local tests (see the discussion around Corollary 4.2), Proposition 4.2 can be used to test for $d_{0}$ sequentially, namely, testing for $H_{0}:d_{0}=r$ starting with $r=p$ and subsequently decreasing $r$ by $1$ till the null hypothesis is not rejected. Let $\widehat{d}_{0}$ be the resulting estimator.

Corollary 4.4**.**

Suppose that the assumptions of Proposition 4.2 hold. Let $\widehat{d}_{0}$ be the estimator of $d_{0}$ defined above through the sequential testing procedure when using a significance level $\alpha=\alpha_{T}$ in testing. If $\alpha_{T}\to 0$ and $(-\log\alpha_{T})/Th\to 0$, then $\widehat{d}_{0}\to_{p}d_{0}$.

Proof.

The result can be proved as in e.g. Theorem 4.3 in Fortuna (2008), by noting that the two conditions provided in the statement of the theorem are in fact equivalent. ∎

We now turn to inference about $d_{+}$ and $d_{-}$ in (4.15). Recall the definition of the eigenvalues $\widehat{\gamma}_{i}(u)$ in (4.6), whose squares are the eigenvalues $\widehat{\gamma}_{2,j}(u)$ in (4.8) entering the test statistics $\widehat{\xi}_{r}(u)$ and $\widehat{\xi}_{r}$. From the discussions surrounding (4.6)–(4.9), the $\widehat{d}_{0}$ consecutive eigenvalues $\widehat{\gamma}_{i}(u)$ can be thought to be associated with the $d_{0}$ zero eigenvalues of $F(u)M_{u}F(u)^{\prime}$. If we can estimate the starting index for these consecutive eigenvalues, we could then deduce the numbers $d_{\pm}$ of positive and negative eigenvalues of $F(u)M(u)F(u)^{\prime}$. In the case $\widehat{d}_{0}\geq 1$, the above suggests to consider

[TABLE]

that is, the quantities involving the sums of the $\widehat{d}_{0}$ consecutive eigenvalues $\widehat{\gamma}_{i}(u)$, and to estimate this starting index through

[TABLE]

(Whenever $r$ does not match the starting index, a larger value of $\widehat{\zeta}_{r}$ is expected, since it will be driven by $\widehat{\gamma}_{i}(u)$ associated with the positive or negative eigenvalues of $F(u)M(u)F(u)^{\prime}$.) With the estimated index $\widehat{r}$, it is then natural to set further

[TABLE]

If $\widehat{d}_{0}=0$, the preceding argument does not apply and, in fact, the quantity (4.17) is not even defined. In this case, we suggest to consider

[TABLE]

and to set

[TABLE]

The idea behind this definition and further motivation for using (4.17)–(4.19) can be found in the proof of the following corollary.

Corollary 4.5**.**

Under the assumptions of Corollary 4.4, we have $\widehat{d}_{+}\to_{p}d_{+}$ and $\widehat{d}_{-}\to_{p}d_{-}$.

Proof.

The result follows from the following two observations. First, by Corollary 4.4, we may assume that for large enough $T$, $\widehat{d}_{0}=d_{0}$. Second, the eigenvalues entering (4.17) and (4.20) converge to their population counterparts, so that the relations (4.17) and (4.20) become in the limit, respectively, the relations $\zeta_{r}=|\int_{\cal H}(\gamma_{r}(u)+\ldots+\gamma_{r+d_{0}-1}(u))du|$ and $\eta_{r}=|\int_{\cal H}(\gamma_{1}(u)+\ldots+\gamma_{r}(u))du|+|\int_{\cal H}(\gamma_{r+1}(u)+\ldots+\gamma_{p}(u))du|$. The conclusion follows by observing that the population quantities satisfy (4.18), (4.19) and (4.21) with all the hats removed. ∎

Finally, a natural estimator for $d:=d_{0}+\min(d_{+},d_{-})$ is

[TABLE]

Corollaries 4.4–4.5 imply the consistency of this estimator, which is the main result of this section.

Theorem 4.2**.**

Under the assumptions of Corollary 4.4, we have $\widehat{d}\to_{p}d$.

4.3 Global dimension tests under interval splitting

The global pseudo dimension test was developed in Section 4.2 under the assumption (4.15) for a subinterval $\mathcal{H}\subset(0,1)$. When global testing is to be performed on (0,1), we suggest to carry out the introduced global test over subintervals of refined partitions. In this section, we describe how the procedure is carried out and the results can be interpreted.

We shall need the following basic property of the global test statistic $\widehat{\xi}_{r,{\cal H}}$ in (4.16). Suppose $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ are two disjoint intervals such that

[TABLE]

Then, since $\int_{\mathcal{H}}=\int_{\mathcal{H}_{1}}+\int_{\mathcal{H}_{2}}$ and $|\mathcal{H}|=|\mathcal{H}_{1}|+|\mathcal{H}_{2}|$, it follows from the definition (4.16) of $\widehat{\xi}_{r,{\cal H}}$ that

[TABLE]

In particular, in view of Proposition 4.2, under $d_{0}(u)\equiv r,\;u\in\mathcal{H}_{i}$,

[TABLE]

where $Z_{i}$’s can be considered independent because the VC model involves independent variables across $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$. The relations (4.24) and (4.25) then imply

[TABLE]

which is consistent with what is expected under $d_{0}(u)\equiv r,\;u\in\mathcal{H}=\mathcal{H}_{1}+\mathcal{H}_{2}$ by Proposition 4.2. Such considerations will allow having some consistency over refined partitions in the sense described below. We first consider the case of $d_{0}(u)$, and then discuss those of $d_{\pm}(u)$ and $d(u)=d(M(u))$.

Thus, let

[TABLE]

form refined partitions of $(0,1]$ for $k\geq 0$. Let $\widehat{d}_{0,i}^{(k)}$ be the global estimator of $d_{0}$ over $\mathcal{H}_{i}^{(k)}$ by using the test statistic $\widehat{\xi}_{r,\mathcal{H}_{i}^{(k)}}$. In view of (4.24)–(4.26) and for the sake of consistency, when estimating $d_{0}$ over finer partitions, we suggest to adjust the normal critical value for comparing $\widehat{\xi}_{r,\mathcal{H}_{i}^{(k)}}$. More precisely, if $c_{\alpha}^{(0)}=c_{\alpha}$ is a normal critical value used at the level $k=0$, we use the critical value $c_{\alpha}^{(k)}=2^{-k/2}c_{\alpha}$ at level $k$.

As a consequence of the choice of the critical values, there are only the following three possibilities for estimators $\widehat{d}_{0,i}^{(k)}$ when refining a partition:

• (P1)

$\widehat{d}_{0,i}^{(k)}=r,\;\widehat{d}_{0,2i-1}^{(k+1)}=r,\;\widehat{d}_{0,2i}^{(k+1)}=r$ ;

• (P2)

$\widehat{d}_{0,i}^{(k)}=r$ and either $\widehat{d}_{0,2i-1}^{(k+1)}=r,\;\widehat{d}_{0,2i}^{(k+1)}<r$ or $\widehat{d}_{0,2i-1}^{(k+1)}<r,\;\widehat{d}_{0,2i}^{(k+1)}=r$ ;

• (P3)

$\widehat{d}_{0,i}^{(k)}=r$ and either $\widehat{d}_{0,2i-1}^{(k+1)}>r,\;\widehat{d}_{0,2i-1}^{(k+1)}>\widehat{d}_{0,2i}^{(k+1)}=r$ or $\widehat{d}_{0,2i}^{(k+1)}>r,\;\widehat{d}_{0,2i}^{(k+1)}>\widehat{d}_{0,2i-1}^{(k+1)}$.

Indeed, let us explain the first two of these possibilities, and also indicate a case which is not one of the possibilities.

The possibility (P1) arises in the following scenario: one has $\widehat{d}_{0,i}^{(k)}=r$ when $\widehat{\xi}_{j,\mathcal{H}_{i}^{(k)}}>2^{-k/2}c_{\alpha}$ for $j=p,p-1,\ldots,r+1$ and $\widehat{\xi}_{r,\mathcal{H}_{i}^{(k)}}\leq 2^{-k/2}c_{\alpha}$. Similarly, $\widehat{d}_{0,2i-1}^{(k+1)}=r$ and $\widehat{d}_{0,2i}^{(k+1)}=r$ when $\widehat{\xi}_{j,\mathcal{H}_{2i-1}^{(k+1)}},\widehat{\xi}_{r,\mathcal{H}_{2i}^{(k+1)}}>2^{-(k+1)/2}c_{\alpha}$ for $j=p,p-1,\ldots,r+1$ and $\widehat{\xi}_{r,\mathcal{H}_{2i-1}^{(k+1)}},\widehat{\xi}_{r,\mathcal{H}_{2i}^{(k+1)}}\leq 2^{-(k+1)/2}c_{\alpha}$. This is consistent with the special case of (4.24),

[TABLE]

in the sense that the relationships of the two summands of (4.28) to the respective critical values is the same as that for their sum.

The possibility (P2), on the other hand, corresponds to the scenario when $\widehat{\xi}_{i,\mathcal{H}_{2i-1}^{(k+1)}},\widehat{\xi}_{j,\mathcal{H}_{2i}^{(k+1)}}>2^{-(k+1)/2}c_{\alpha}$ for $j=p,p-1,\ldots,r+1$, but then either $\widehat{\xi}_{r,\mathcal{H}_{2i-1}^{(k+1)}}\leq 2^{-(k+1)/2}c_{\alpha}$ and $\widehat{\xi}_{r,\mathcal{H}_{2i}^{(k+1)}}>2^{-(k+1)/2}c_{\alpha}$ or $\widehat{\xi}_{r,\mathcal{H}_{2i-1}^{(k+1)}}>2^{-(k+1)/2}c_{\alpha}$ and $\widehat{\xi}_{r,\mathcal{H}_{2i}^{(k+1)}}\leq 2^{-(k+1)/2}c_{\alpha}$. The case that is excluded from the possibilities listed above, is $\widehat{d}_{0,2i-1}^{(k+1)}<r$ and $\widehat{d}_{0,2i}^{(k+1)}<r$ (and $\widehat{d}_{0,i}^{(k)}=r$) which would happen if $\widehat{\xi}_{r,\mathcal{H}_{2i-1}^{(k+1)}},\widehat{\xi}_{r,\mathcal{H}_{2i}^{(k+1)}}>2^{-(k+1)/2}c_{\alpha}$ but is impossible in view of (4.28). In our experiments with simulated and real data, the possibilities (P1) and (P2) seem to be the most common.

Figure 1, left plot, presents global testing results under splitting for one realization of a VC model (Model 5 from Section 6.1). In this plot, the estimate of $\widehat{d}_{0,1}^{(0)}$ for $(0,1]$ is presented at the point 1/2 of the $x$-axis which is the midpoint of $(0,1]$. The estimator $\widehat{d}_{0,1}^{(1)}$ and $\widehat{d}_{0,2}^{(1)}$ are presented at points 1/4 and 3/4 respectively, which are the midpoints of the intervals $(0,1/2]$ and $(1/2,1]$. The presentation is continued in the same way, till the level $k=2$ is reached and the estimates $\widehat{d}_{0,1}^{(2)},\widehat{d}_{0,2}^{(2)},\widehat{d}_{0,3}^{(2)},\widehat{d}_{0,4}^{(2)}$ are presented at the finest considered level.

The estimators $\widehat{d}_{0,i}^{(k)}$ lead to the corresponding estimates $\widehat{d}_{\pm,i}^{(k)}$ and $\widehat{d}_{i}^{(k)}$. The latter could be presented similarly as in the described left plot of Figure 1. This is illustrated in the right plot of Figure 1.

5 Estimation of stationary subspace

We turn here to the estimation of a stationary subspace of the VC model, which is related to pseudo null spaces, given in Definition 3.1, of matrices $M(u)$ through (3.13). We shall not provide here any formal statistical tests related to a stationary subspace but rather make a number of related comments, inspired by the developments in Section 3.2. In Section 5.1, we introduce a particular class of local stationary subspaces. In Section 5.2, a graph-based method is provided that aims at forming clusters of the many local stationary subspaces introduced in Section 5.1. Finally, a technique to select one subspace out of the clusters is provided and this serves as our estimate of a stationary subspace.

5.1 (1,1)-local stationary subspaces

The results of Section 3.2 show that pseudo null spaces have a rich structure and are typically not unique, for a given matrix. Motivated by the proof of Proposition 3.1, we shall restrict our attention to their special cases, given in the following definition. The definition and subsequent developments use the notation of Section 3, namely that of $M$, $d_{0}$, $d_{\pm}$, $\lambda_{i}$, $\lambda_{0,i}$, $\lambda_{\pm,i}$, $s_{i}$, $s_{0,i}$, $s_{\pm,i}$.

Definition 5.1**.**

Let $M$ be a symmetric matrix and suppose $\min(d_{+},d_{-})\geq 1$. A $(1,1)$-pseudo null space of $M$ is defined as

[TABLE]

for a fixed $\alpha_{i},\beta_{i}$, where $p(i)\in\{1,\ldots,d_{+}\}$ are different across $i$, $n(i)\in\{1,\ldots,d_{-}\}$ are different across $i$ and

[TABLE]

When $M=M(u)$ with $M(u)$ as in (2.2), a (1,1)-pseudo null space $\mathcal{P}(M(u))$ will be called a (1,1)-local stationary subspace and denoted as $\mathcal{B}_{(1,1)}(u)$.

The fact that ${\cal P}_{(1,1)}(M)$ defines a pseudo null space for $M$ follows as in the proof of Proposition 3.1.

Remark 5.1**.**

The total number of $(1,1)$-pseudo null spaces of $M$ is

[TABLE]

where the first two terms account for the selection of eigenvectors associated with the positive and negative eigenvalues, and the last term for pairing them off. Depending on the values of $d_{\pm}$, the total number (5.3) can be quite large: e.g. with $d_{+}=3$ and $d_{-}=5$, the number is $60$.

Remark 5.2**.**

The prefix “$(1,1)$-” in Definition 5.1 refers to the fact that a pseudo eigenvector of a pseudo null space is constructed by taking one (1) eigenvector $s_{+,p(i)}$ associated with the positive eigenvalues and one (1) eigenvector $s_{-,n(i)}$ associated with the negative eigenvalues. More elaborate constructions are possible as well, for example, by taking two (2) eigenvectors associated with the positive eigenvalues and one (1) eigenvector with the negative ones, as in Example 3.1, which could be called a $(2,1)$-pseudo null space. In this work though, we shall consider only $(1,1)$-pseudo null spaces.

If $\widehat{M}$ estimates $M$ with the analogous estimators $\widehat{d}_{0}$, $\widehat{d}_{\pm}$, $\widehat{\lambda}_{i}$, $\widehat{\lambda}_{0,i}$, $\widehat{\lambda}_{\pm,i}$, $\widehat{s}_{i}$, $\widehat{s}_{0,i}$, $\widehat{s}_{\pm,i}$ of the respective quantities, we would similarly like to have the sample counterparts of the $(1,1)$-pseudo null space in (5.1). Definition 5.1 may not, however, extend directly to the sample quantities since $\widehat{d}_{\pm}$ may not necessarily represent the actual number of positive/negative eigenvalues $\widehat{\lambda}_{\pm,i}$ and hence the condition (5.2) may not be satisfied. To deal with this possibility, we define the sample counterparts $\widehat{\cal P}_{(1,1)}(M)$ in the same way as in (5.1) by using the quantities with the hats, except that $\widehat{d}_{+}$ and $\widehat{d}_{-}$ are replaced by $\widetilde{d}_{+}$ and $\widetilde{d}_{-}$, which are defined as

[TABLE]

and

[TABLE]

where we assumed that the eigenvalues $\widehat{\lambda}_{i}$ appear in the non-decreasing order. That is, we define a sample $(1,1)$-pseudo null space as

[TABLE]

where $p(i)\in\{1,\ldots,\widetilde{d}_{+}\}$ are different across $i$, $n(i)\in\{1,\ldots,\widetilde{n}_{+}\}$ are different across $i$ and

[TABLE]

Replacing $\widehat{M}$ above by $\widehat{M}(u)$ from (4.1) for $u\in(0,1)$, estimators $\widehat{\mathcal{B}}_{(1,1)}(u)=\mathcal{P}_{(1,1)}(\widehat{M}(u))$ can be defined for (1,1)-local stationary subspaces. Techniques to cluster these (1,1)-local stationary subspaces is discussed next.

5.2 Clustering (1,1)-local stationary subspaces

According to (3.13), there is a stationary subspace $\mathcal{B}_{1}$ for a VC model if there is at least one identical stationary subspace for all $u\in(0,1)$. In Section 5.1, we defined (1,1)-local stationary subspaces $\mathcal{B}_{(1,1)}(u)$ whose number can already be quite large for a fixed $u$; see Remark 5.1. A natural possibility in defining a candidate for a stationary subspace $\mathcal{B}_{1}$ is to consider all (1,1)-local stationary subspaces across $u$’s and select a subspace representing their “majority” in some suitable sense. In light of this observation, using clustering seems natural and this approach is pursued here on the estimated (1,1)-local stationary subspaces.

More specifically, we discuss a graph-based approach to clustering the many (1,1)-local stationary subspaces, or equivalently, the many (1,1)-pseudo null spaces using distances that are computed based on canonical angles between spaces. In addition, a method to select one (1,1)-local stationary subspace out of the many is discussed and this is considered as our estimate of a stationary subspace.

As in Section 5.1, we let $\widehat{\mathcal{B}}_{(1,1)}(u)=\mathcal{P}_{(1,1)}(\widehat{M}(u))$ denote a sample (1,1)-local stationary subspace of the matrix $\widehat{M}(u)$ from (4.1). Let $\widehat{\mathcal{B}}_{(1,1)}=\bigcup_{u}\widehat{\mathcal{B}}_{(1,1)}(u)$ be the union of all (1,1)-local stationary subspaces over $u$. In practice, the union $\bigcup_{u}$ is replaced by a set of discrete points $\{u_{1},u_{2},\ldots,u_{n_{u}}\}$ in $(0,1)$.

Consider a graph $G=(V,E)$ where every vertex $v\in V$ corresponds to a (1,1)-local stationary subspace in $\widehat{\cal B}_{(1,1)}$. The adjacency matrix $E$ will be defined in terms of a distance between (1,1)-local stationary subspaces using canonical angles. Let $\mathscr{B}_{1}$ and $\mathscr{B}_{2}$ be two (1,1)-local stationary subspaces in $\widehat{\cal B}_{(1,1)}$ of dimensions $d_{1}$ and $d_{2}$, respectively. Letting $d=\min(d_{1},d_{2})$, the canonical angles computed between these two spaces are given by $\theta_{1}\leq\theta_{2}\leq\ldots\leq\theta_{d}$, where

[TABLE]

The vectors $x_{1},x_{2},\ldots,x_{d}$ and $y_{1},y_{2},\ldots,y_{d}$ are called canonical vectors. We measure the distance between spaces $\mathscr{B}_{1}$ and $\mathscr{B}_{2}$ as $\max_{1\leq i\leq d}\theta_{i},\;$ and define the adjacency matrix $E=(e_{i,j})$ for $i,j=1,2,\ldots,|V|$, as

[TABLE]

where $\theta(v_{i},v_{j})$ is the maximum canonical angle between the (1,1)-local stationary subspaces corresponding to vertices $v_{i}$ and $v_{j}$, and $\theta_{0}$ is a threshold. The choice of $\theta_{0}$ dictates the number of clusters estimated with more being formed for lower values of $\theta_{0}$. In our numerical work, we set $\theta_{0}=20^{\circ}$.

Finally, in order to obtain the clusters of vertices in the graph $G$, we utilize the Walktrap algorithm of Pons and Latapy (2005). Here, a transition probability matrix $P=(p_{ij})$ is constructed with $p_{ij}=\frac{A_{ij}}{d(i)}$ where $A=(A_{ij})$ denotes the adjacency matrix of $G$ and $d(i)$ denotes the degree of vertex $v_{i}$. A random walk process defined on this graph $G$ is based on the powers of the matrix $P$, that is, the probability of moving from vertex $v_{i}$ to $v_{j}$ through a random walk of length $t$ is given by $P^{t}_{ij}$. The closeness of vertices in the graph is measured by these probabilities from the observation that if two vertices $v_{i}$ and $v_{j}$ are in the same cluster, $P^{t}_{ij}$ must be high.

Let $C_{1},C_{2},\ldots,C_{L}$ be the $L$ cluster of vertices produced by the Walktrap algorithm that have a size of at least 3 vertices. We first obtain the centers $\{c_{1},c_{2},\ldots,c_{L}\}$ of these $L$ clusters by computing the sine of the maximum canonical angle,

[TABLE]

The (1,1)-local stationary subspaces corresponding to the $L$ cluster centers are considered as the candidate stationary subspaces returned by our method. Additionally, in order to select a single (1,1)-local stationary subspace (our stationary subspace estimator) out of these $L$ representative subspaces, we assess the “denseness” of each cluster through

[TABLE]

where $d(v)$ denotes the degree of vertex $v$ within cluster $C_{l}$. We then identify the cluster with maximum $T(\cdot)$ among the $L$ clusters and select the final (1,1)-local stationary subspace corresponding to the center of that cluster. That is, we select our stationary subspace estimate as the (1,1)-local stationary subspace corresponding to the cluster center $c_{s}$ in cluster $C_{s}$, where

[TABLE]

Example 5.1**.**

To illustrate the above technique, we consider the VC model (1.2) with $p=3$, i.i.d. $Y_{t}\sim N(0,I_{3})$, $T=1000$ and $A^{2}(u)=\mbox{diag}(2+0.5\sin(2\pi u),3-\sin(2\pi u),1.5+\sin(2\pi u))$. The stationary subspace dimension for this model is $d=1$ and the true (1,1)-local stationary subspaces are given by $(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0)$, $(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$, $(\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}},0)$ and $(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0)$ for all $u$’s. We generated one realization of the series $X_{t}$ based on this model and obtained all the estimated (1,1)-local stationary subspaces across a set of points $\{0.04,0.08,0.12,\ldots,0.96\}$. Figure 2 depicts a 3D plot that includes the population (1,1)-local stationary subspaces (solid circles), the estimated (1,1)-local stationary subspaces (crosses), the cluster centers $c_{l}$ for $l=1,2,3,4$ of the 4 clusters (open circles) and the selected (1,1)-local stationary subspace based on (5.13) (solid square, marked as VC Final). Finally, the estimated stationary subspace from DSSA (Sundararajan and Pourahmadi (2018)) is plotted along with the other spaces (diamond).

In Table 1, we list the cluster centers of the 4 main clusters in Figure 2, the sizes of those 4 clusters and the proportions of $u\in\{0.04,0.08,0.12,\ldots,0.96\}$ with points in the respective clusters. Observe that the (1,1)-local stationary subspace selected as the stationary subspace based on (5.13) lies in the most “dense” cluster and the DSSA stationary subspace also lies in the same cluster. Here, the 4 biggest clusters formed by the method comprise roughly 84% of the total number of (1,1)-local stationary subspaces. The (1,1)-local stationary subspace selected as the stationary subspace and identified as the solid square point in Figure 2 lies in the biggest cluster that contains roughly 24% of the (1,1)-local stationary subspaces. This subspace also lies in the cluster with the highest proportion of $u$’s with points in that cluster.

6 Simulation study

In Section 6.1, we evaluate the empirical performance of the proposed method in estimating the dimension of a stationary subspace for several VC models. In Section 6.2, we assess the ability of the proposed method in estimating a stationary subspace using a few discrepancy measures.

6.1 Dimension estimation comparison

We first consider several VC models (1.2), characterized through $A(u)$, for which the dimensions $d_{0}(u)=d_{0}$, $\min(d_{+}(u),d_{-}(u))=\min(d_{+},d_{-})$ and $d(u)\equiv d$ do not depend on $u$ (and neither do the local stationary subspaces). The model matrices $A(u)$, the respective matrices in (2.2) and the dimensions $p$, $d_{0}$, $\min(d_{+},d_{-})$ are:

Model 1: $p=3$, $d_{0}=0$, $\min(d_{+},d_{-})=1$, $d=1$,

[TABLE]

Model 2: $p=3$, $d_{0}=2$, $\min(d_{+},d_{-})=0$, $d=2$,

[TABLE]

Model 3: $p=4$, $d_{0}=1$, $\min(d_{+},d_{-})=1$, $d=2$,

[TABLE]

$M(u)=\mbox{diag}(0,\sin(2\pi u),1-2u,-\sin(2\pi u))$.

Model 4: $p=4$, $d_{0}=3$, $\min(d_{+},d_{-})=0$, $d=3$,

[TABLE]

Models 1 and 2 are 3-dimensional and have diagonal $A^{2}(u)$. Models 3 and 4 are 4-dimensional with non-diagonal $A^{2}(u)$.

We suppose that the models above are Gaussian with i.i.d. $N(0,I_{p})$ vectors $Y_{t}$. In the simulations we take $T\in\{200,500,1000,2000\}$ as the sample sizes. In estimating $A^{2}(u)$ in (4.1), bandwidth choices $h$ ranging from $T^{-0.1}$ to $T^{-0.5}$ were attempted and the best results were obtained for $h\in(T^{-0.3},T^{-0.4})$. We fix $h=T^{-0.35}$ and present the results for this choice. We focus on testing for $\mathcal{H}=(0,1)$ only and use 100 Monte Carlo replications.

We compare the performance of the proposed dimension estimator $\widehat{d}$ given in (4.22) with the sequential estimation procedure provided in Section 2.2.5 of Sundararajan and Pourahmadi (2018). The method found in the latter work will be referred to as DSSA and the proposed method will be denoted as VC. The estimation results for the two methods and the four considered models are presented through violin plots in Figure 3. Violin plots are intended as a somewhat qualitative visualization of the results – perhaps slightly more informative histograms for the results can be found in Sundararajan et al. (2019). The plots show that estimation improves with increasing sample size $T$, and that the proposed VC method performs better than the competing DSSA method in detecting the true dimension $d$.

We now turn to VC models whose dimensions depend on $u$’s. We consider the following models:

Model 5: $p=3$, $d_{0}(u)=2$ and $d(u)=2$ if $u<0.5$ and $d_{0}(u)=1$ and $d(u)=2$ if $u>0.5$,

[TABLE]

Model 6: $p=3$, $d_{0}(u)=3$ and $d_{u}=3$ if $u<0.5$ and $d_{0}(u)=2$ and $d(u)=2$ if $u>0.5$,

[TABLE]

The entries 2.0901 in Model 5 and 3.125 in Model 6 ensure smoothness of $A^{2}(u)$ at $u=0.5$.

We report the estimation results for the models in Figures 4 and 5. The plots across the two figures are analogous and hence only those in Figure 4 will be explained. The top two plots in Figure 4 concern the estimates of $d_{0}(u)$, and the bottom two plots are those of $d(u)$. The structure of the left plots is similar to that of Figure 1. That is, the estimates over $\mathcal{H}=(0,1)$ are depicted at $u=0.5$, over $\mathcal{H}=(0,0.5)$ at $u=0.25$ and those over $\mathcal{H}=(0.5,1)$ at $u=0.75$. The only difference here is that the results are reported over 100 realizations in the form of violin plots. The plots on the right, on the other hand, present the local estimates at $u$’s indicated on the horizontal axis, according to the methods discussed in Section 4.1.

Several observations are in place regarding Figures 4 and 5. Note that the estimates are generally sensitive to the choice of the subinterval, (0,0.5) or (0.5,1), in the direction of the true dimension value. The estimates of the local dimensions also tend to capture their variability across the different $u$’s.

6.2 Subspace estimation comparison

We compare here the performance of the proposed method in estimating a subspace (from Section 5.2) to DSSA of Sundararajan and Pourahmadi (2018), in terms of three discrepancy measures. The first measure concerns departure from stationarity based on the sizes of the DFT covariances as given in Eq. (9) of Sundararajan and Pourahmadi (2018). More precisely, for an estimated stationary subspace process $Y_{t}=\widehat{B}_{1}^{{}^{\prime}}X_{t}$, we set

[TABLE]

where $||A||_{F}$ denotes the Frobenius norm of a matrix $A\in\mathbb{R}^{d\times d}$, $\Re(\cdot)$ and $\Im(\cdot)$ denote the entrywise real and imaginary parts respectively, and $\widehat{\Gamma}_{r}^{Y}$ is the $d\times d$ lag-$r$ DFT sample autocovariance matrix given by

[TABLE]

with $\omega_{k}=2\pi k/T$ referring to a Fourier frequency, $J_{Y}(\cdot)$ being the discrete Fourier transform (DFT) of the $d$-variate series $Y_{t}$ and $J_{Y}(\cdot)^{\ast}$ denoting the complex conjugate transpose. The number of DFT covariance lags $m$ in (6.1) is fixed to 3.

The second measure is based on the relation (2.3). For any candidate subspace $\widehat{B}_{1}$, we set

[TABLE]

The last measure is based on canonical angles computed between the population and estimated subspaces $\mathcal{B}_{1}$ and $\widehat{\mathcal{B}}_{1}$ that are spanned by the columns of the $p\times d$ matrices $B_{1}$ and $\widehat{B}_{1}$, respectively. As in (5.8), let $\theta_{1}\leq\theta_{2}\leq\ldots\leq\theta_{d}$ be the $d$ canonical angles between the spaces $\mathcal{B}_{1}$ and $\widehat{\mathcal{B}}_{1}$. Then, set

[TABLE]

In the simulations here, we consider the same four models, Models 1–4, as in Section 6.1. We first present estimation results for the measure $D_{3}$ in Table 2. For our method, labeled as VC in the table, the stationary subspace estimate is taken as discussed in Section 5.2. At the population level, the corresponding stationary subspace is also selected using the same technique. In the case when such multiple subspaces are available at the population level, we compute the distance $D_{3}$ to all of them and then take the minimum. For the DSSA method, we take the subspace estimate as defined in Sundararajan and Pourahmadi (2018).

As seen from the table, the VC method generally performs better than DSSA in yielding smaller average distances over multiple realizations. We stress here again that the VC method is computationally much less intensive than DSSA.

We now turn to the measures $D_{1}$ and $D_{2}$. Here, we shall examine the VC and DSSA methods through these measures from a different perspective. The relevant plots can be found in Figure 6 for Model 1, associated with the indicated sample sizes. In Figure 6, top panel, the horizontal solid circles in the plots represent the measure $D_{1}(\widehat{B}_{1})$ for the DSSA estimate $\widehat{B}_{1}$, averaged over multiple realizations. The other two curves correspond to the measure $D_{1}(\widehat{B}_{1}(u))$ computed for the estimates of $\widehat{B}_{1}(u)$ of (1,1)-local stationary subspaces from Section 5, either averaged over multiple realizations (triangles) or with the minimum value taken (squares). The interpretation of the plots of Figure 6, bottom panel, is analogous but for the measure $D_{2}$. The figures show that VC method performs better than DSSA even “locally” for most values of $u$ under measure $D_{2}$ whereas DSSA performs better under measure $D_{1}$. Analogous figures for Models 2–4 can be found in Sundararajan et al. (2019). The results are similar to Figure 6, especially for larger sample sizes.

7 Application to BCI and EEG data

Brain-Computer Interface (BCI) aims at connecting the human brain and the computer in a non-invasive manner. During the BCI study used here, individuals are asked to imagine movements with their left and right hands and these are referred to as trials. The trials are interspersed with break periods. The multivariate EEG brain signal is recorded during the entire course of the experiment and the objective is to associate the movements imagined by the individuals with the corresponding EEG signal. Note that the EEG signal is recorded through different channels (locations on the scalp) and each channel constitutes a component of the multivariate signal.

Several works including Ombao et al. (2005), von Bünau et al. (2009), Sundararajan and Pourahmadi (2018) and Sundararajan et al. (2017) have treated the multivariate EEG signal from such experiments as a nonstationary time series. The advantages of finding a stationary subspace of the observed EEG signal are discussed in the above mentioned works. The key observation made was that nonstationary sources in the brain signal were associated with variations in the mental state that are unrelated to the motor imagery task. Hence, relying solely on stationary components was seen to improve classification accuracy.

We study the classification performance of the proposed VC method using the BCI Competition IV 444See http://www.bbci.de/competition/iv/ dataset II in Naeem et al. (2006). It consists of EEG signals from 9 subjects performing 4 different motor imagery tasks: 1-left hand, 2-right hand, 3-feet and 4-tongue. We analyze the EEG signals only from classes 1 and 2 and treat the problem as a two-group classification. The continuous signal was sampled at discrete time steps at the sampling rate of 250 Hz where 1 second corresponds to 250 observations on the digital signal scale. The signal was recorded through 22 electrodes over the course of the experiment and the signal was band-pass filtered as in Lotte and Guan (2011). The experiment involved 144 trials for each subject wherein 72 trials belonged to Class 1 (left hand) and the other 72 to Class 2 (right hand). Every trial is followed by an adequate resting period for the subject before the start of the next one. In each trial, we use the observations between 0.5 seconds to 2.5 seconds after the cue instructing the subject to perform the motor imagery task. More precisely, for trial $j$ where $j=1,2,\ldots,144$, this interval comprises of 500 observations on the digital signal scale, denoted by $X_{t}^{(j)}$, $t=1,2,\ldots,500$.

We first restrict attention to 5 EEG electrodes and treat the input signal to have dimension $p=5$. These are 5 locations that can be viewed as representatives of the different regions on the brain, namely, Frontal (Fz), Pre-Frontal (Pz) and Cortical (C3, C4, Cz).555See http://www.bbci.de/competition/iv/desc_ 2a.pdf for additional details on the dataset. We use the VC method to obtain a $d$ dimensional stationary process where $d<p$. For every trial $j=1,2,\ldots,144$, we have $\{X_{t}^{(j)}\}\in\mathbb{R}^{p}$ as the observed multivariate signal and $\{Y_{t}^{(j)}\}\in\mathbb{R}^{d}$ as the stationary transformation.

We first report on the estimated pseudo dimension $d$ for the observed signals for the 9 subjects in this study labeled S1, S2,$\ldots$, S9. For subjects S3, S5 and S8, the percentage of times the candidate dimensions ($d=0,1,2,3,4$) were estimated by the 2 competing methods, DSSA or VC, out of the 144 trials is provided in the left panel of Figure 7. This plot also includes the estimates of $d$ over the first and second halves of the data denoted as VC (First) and VC (Second). It is noted that DSSA always provides a lower estimate of $d$ than the VC method. Similar plots for all 9 subjects can be found Sundararajan et al. (2019).

Given the $d$-variate stationary processes $Y_{t}^{(j)}$, we aim to find differences between the two classes (0 and 1) based on the covariance structure. For a given subject, this is achieved by computing the average spectral density matrices for the two classes over the Fourier frequencies:

[TABLE]

where $g_{j}(\omega_{k})$ is the estimated $d\times d$ spectral matrix for trial $j$ using observations $\{Y_{t}^{(j)}\}$, $n_{i}=72$, for $i=0,1$ and $\omega_{k}=\frac{2\pi k}{500}$, $k=1,2,\ldots,500$, are the Fourier frequencies.

In order to train a classifier, for every trial $j=1,2,\ldots,144$, a distance vector $p_{j,AB}=(p_{0,j,AB},p_{1,j,AB})$ is computed, where

[TABLE]

and $||\cdot||_{F}$ is the Frobenius norm of a matrix. It measures the distance to the center of each of the two classes. This distance measure serves as our two-dimensional feature vector to be used in constructing a logistic regression classifier and assessing its out-of-sample classification accuracy.

Table 3 (left) shows the out-of-sample classification accuracies for three subjects corresponding to $d=1,2,3,4$. A similar table containing results for all 9 subjects can be found in Sundararajan et al. (2019). The accuracy rates reflect a comparable average performance in the two competing methods. Finally, the accuracy rate increases as the pseudo dimension $d$ increases from 1 to 4, a phenomenon witnessed and discussed in von Bünau et al. (2009) and Sundararajan and Pourahmadi (2018). In dealing with brain signals from healthy individuals in this experiment, the nonstationarity is believed to be caused by artifacts such as as fatigue, physical movement, blinking. Hence, more stationary sources means greater elimination of nonstationary sources in the signal that are unrelated to the experimental task at hand.

The results in Figure 7 (right) and Table 3 (right) are analogous to those on their respective left panels but taking the input signal to have dimension $p=22$, that is, without restricting attention to 5 EEG electrodes. For the out-of-sample classification accuracies, the values $d=7,9,11,13$ are considered. The results in Table 3 (right) indicate a better performance of VC in comparison to results on the same dataset (with $p=22$) in Sundararajan and Pourahmadi (2018).

It must be noted that the results in Table 3 is based on an estimated stationary subspace for a pre-fixed dimension $d$. The estimates of the dimensions of the stationary subspaces from (4.22), namely $\widehat{d}_{0}$, $\widehat{d}_{+}$ and $\widehat{d}_{-}$, could potentially be different across different trials. For each trial, we obtained these estimates and compared them with the fixed $d$ and investigated results for several different combinations of $\widehat{d}_{0}$, $\widehat{d}_{+}$ and $\widehat{d}_{-}$ that results in $d=\widehat{d}_{0}+\min(\widehat{d}_{+},\widehat{d}_{-})$. Having witnessed very similar results for the various combinations, we only present, in Table 3, the case wherein $d=\widehat{d_{0}}$.

8 Concluding remarks

Our goal in this work is to (i) study existence of linear combinations of components of a multivariate nonstationary process which are stationary, and (ii) find the number of such stationary linear combinations. The true nature of the problem and richness of its solution present themselves naturally when the general dependence setup in Sundararajan and Pourahmadi (2018) is abandoned in favor of heterogenous independent observations. In this simplified setup (von Bünau et al. (2009)), solution of the problem reduces to the study of inertia or signs of the eigenvalues and the corresponding eigenvectors of certain symmetric time-varying matrices constructed from varying covariance or heterogeneity of the vector observations. This enable us to provide a direct linear-algebraic method to construct stationary subspaces which outperforms the earlier computationally more expensive optimization-based SSA solutions.

Several directions related to this work could be explored in the future. The developed framework involving pseudo spaces and dimensions is general enough to apply to zero mean locally stationary processes, when working with their time-frequency spectra. Another possibility is to try to combine the models of this work and Düker et al. (2019), so that both mean and the covariance are allowed to vary in time. Yet another direction is to explore connections to multivariate stochastic volatility models that are concerned with modeling the changing covariance across time.

Appendix A Proofs

The appendix concerns the technical aspects of this work, including various assumptions and proofs of some of the stated results.

A.1 Assumptions and proof of Proposition 4.1

We use the following assumptions for Proposition 4.1, labeled according to the quantities they concern.

• (Y1)

$Y_{t}$, $t=1,\ldots,T$, are i.i.d. random vectors with i.i.d. entries, $\mathbb{E}(Y_{t})=0$ and $\mathbb{E}(Y_{t}Y_{t}^{\prime})=I_{p}$.

• (Y2)

The entries of $Y_{t}$ have finite absolute moments of order $4+\epsilon$ for some $\epsilon>0$.

• (K)

The kernel function $K$ is even (i.e. $K(u)=K(-u)$, $u\in\mathbb{R}$), has bounded support $(-S,S)$, where it is positive, and is continuously differentiable on $(0,S)$. Furthermore, $\int_{\mathbb{R}}K(u)du=1$.

• (A1)

The matrix $A^{2}(u)$ is positive definite for $u\in(0,1)$.

• (A2)

The matrix $A^{2}(u)$ is continuously differentiable for $u\in(0,1)$.

We note that (A1) and (A2) imply the continuous differentiability of $A(u)$, $A(u)^{-1}$ and $A(u)^{-2}$ for $u\in(0,1)$.

Proof of Proposition 4.1: In view of (2.2) and (4.1), write

[TABLE]

where

[TABLE]

We will show that $(Th)^{1/2}S_{j}\to 0$ for $j=2,4$, and $\to_{p}0$ for $j=3$, and that $S_{1}$ yields the normal limit in (4.2). This shall prove the convergence (4.2). In what follows, $C$, $C^{\prime}$ and so on, denote generic positive constants that can change from line to line.

Since the derivatives of the entries of $A(u)$ (and $A^{2}(u)$) are assumed continuous and hence bounded on $(0,1)$, it follows that $|S_{4}|\leq C/T$ and hence indeed $(Th)^{1/2}S_{4}\to 0$, since $T\to\infty$ and $h\to 0$. Since by assumptions, $Y_{t}$ are i.i.d. and their entries have finite fourth moments, and $A(u)$ is bounded on $(0,1)$, we have

[TABLE]

Thus, $S_{3}=O_{p}(T^{-1/2})$ and $(Th)^{1/2}S_{3}=O_{p}(h^{1/2})=o_{p}(1)$. For $S_{2}$, write $S_{2}=S_{2,1}+S_{2,2}$, where

[TABLE]

Since the entries of $A^{2}(u)$ have bounded first derivatives, we have

[TABLE]

where $G_{1}(u)=|u|K(u)$, $G_{1,h}(u)=h^{-1}G_{1}(h^{-1}u)$ and we used Lemma A.1 below in the last step. Thus, $(Th)^{1/2}S_{2,1}\to 0$ since $Th^{3}\to 0$ by assumption. Similarly, $S_{2,2}=O((Th)^{-1})$.

Finally, we show that $S_{1}$ gives the desired asymptotic limit. For this, we consider $A(u)^{-1}S_{1}A(u)^{\prime-1}$ and express it as

[TABLE]

where

[TABLE]

We will argue that $S_{1,1}$, $S_{1,2}$ and $S_{1,3}$ are asymptotically negligible for the limit, and that $S_{1,4}$ gives the desired normal limit.

For $S_{1,1}$, let $K_{2}(u)=K^{2}(u)/\|K\|_{2}^{2}$ be a kernel function obtained from the square of $K$ and $K_{2,h}(u)=h^{-1}K_{2}(h^{-1}u)$. As in dealing with $S_{2}$ above, note that $\mathbb{E}\|X_{t}X_{t}^{\prime}-A^{2}(t/T)\|_{2}^{2}$ is bounded by a constant uniformly over $t=1,\ldots,T$. Then,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $G_{2}(u)=u^{2}K_{2}(u)$, $G_{2,h}(u)=h^{-1}G_{2}(h^{-1}u)$ and we used the assumed smoothness of $A(u)^{-1}$. By using Lemma A.1 below, we have $S_{1,1}=O_{p}((h/T)^{1/2})$ and hence $(Th)^{1/2}S_{1,1}=O_{p}(h)=o_{p}(1)$. One can show similarly that $S_{1,2}=O_{p}(1/T^{1/2})$ and $S_{1,3}=O_{p}(1/T^{1/2})$, and hence that $(Th)^{1/2}S_{1,2}$ and $(Th)^{1/2}S_{1,3}$ are negligible as well. We next analyze the last term $S_{1,4}$.

To establish the asymptotic normality of $(Th)^{1/2}S_{1,4}$, it is enough to check the Lyapunov condition and to make sure that the variances converge to the desired quantities. For the Lyapunov condition, note that, with $\delta>0$,

[TABLE]

[TABLE]

where $K_{2+\delta}(u)=K(u)^{2+\delta}/\|K\|_{2+\delta}^{2+\delta}$. The last bound converges to [math] since $Th\to\infty$ and since the sum in the bound converges to $\int_{\mathbb{R}}K_{2+\delta}(u)du$ by Lemma A.1 below. For the convergence of the variances, note that the different entries of the symmetric matrix $Y_{t}Y_{t}^{\prime}-I_{p}$ are independent, which translates into them being uncorrelated in the limit, as stated in the proposition. We thus need to consider the variances of

[TABLE]

for the diagonal and off-diagonal asymptotic variances, respectively. For $S_{\scriptsize\rm diag}$, we have

[TABLE]

where again, the kernel $K_{2}(u)=K^{2}(u)/\|K\|_{2}^{2}$. By Lemma A.1 below, $\mathbb{E}S_{\scriptsize\rm diag}^{2}\to\mu_{4}\|K\|_{2}^{2}$, which matches the diagonal variance in (4.2). For $S_{\scriptsize\rm off}$, we have similarly

[TABLE]

This concludes the proof of the convergence (4.2).

Finally, the consistency of $\widehat{A}^{2}(u)$ follows from the proved convergence to [math] of the terms $S_{1}$ and $S_{2}$ above. This implies the consistency of $\widehat{A}^{-2}(u)$ since $A^{2}(u)$ is assumed to be positive-definite. Since the square root operation is a continuous functional, we also get the consistency of $\widehat{A}(u)$ and $\widehat{A}^{-1}(u)$.   $\Box$

Part $(i)$ of the following auxiliary lemma was used in the proof above. Its part $(ii)$ will be used in Appendices A.2 and A.3 below.

Lemma A.1**.**

Let $G:\mathbb{R}\to\mathbb{R}$ be even ($G(-v)=G(v)$, $v\in\mathbb{R}$), have bounded support $(-S,S)$ and be continuously differentiable on $(0,S)$. Let also $G_{h}(v)=\frac{1}{h}G(\frac{v}{h})$, $v\in\mathbb{R}$, with $h>0$ such that $h\to 0$ and $Th\to\infty$. Then:

$(i)$* for $u\in(0,1)$,*

[TABLE]

$(ii)$* The bound also hold uniformly for $u\in{\cal H}$, where ${\cal H}$ is a closed subinterval in $(0,1)$.*

Proof.

The result $(i)$ follows from: for large enough $T$,

[TABLE]

[TABLE]

since, for example, $Tu-1>STh$ for large enough $T$. The last sum above converges to $\int_{-S}^{S}G(v)dv$ no slower than the given rate since $Th\to\infty$ and $G$ is continuously differentiable. The part $(ii)$ follows in the same way since $Tu-1>STh$ for large enough $T$ uniformly for $u\in{\cal H}$. ∎

A.2 Assumptions and proof of Proposition 4.2

The proof of Proposition 4.2 follows the path taken by Donald et al. (2011); see, in particular, its technical appendix online. We shall try to minimize repetitions by indicating only the key needed assertions. Some of the developments will be somewhat simpler since many of the considered matrices are symmetric. But we shall also need several new auxiliary results to account for the key difference from Donald et al. (2011) in that smoothing through a kernel is carried out here in time $t$ while this was over the values of a random variable in Donald et al. (2011). Some of the auxiliary results for the proof of Proposition 4.2 can be found in Appendix A.3 and the reader interested in proofs might want to look at them first, before going through the arguments in this section.

We assume implicitly as stated in Section 4.2 that ${\cal H}$ is a closed subinterval of $(0,1)$. We shall also need to strengthen Assumption (A2) in Appendix A.1 to:

• (A2g)

The entries of the matrix $A^{2}(u)$ are real analytic for $u\in{\cal H}$.

According to one possible definition, a function $f$ is real analytic if its Taylor series converges to the function $f$ in a neighborhood of each point. As noted for similar assumptions G4, G5 in Donald et al. (2011), analyticity is assumed to have smoothness of the eigenvectors of analytic matrices involving $A^{2}(u)$. It is well known that smoothness of a matrix is not sufficient to have smooth eigenvectors (see, e.g., Kato (1976), Bunse-Gerstner et al. (1991)). Alternatively, the smoothness of the eigenvectors of interest can be assumed.

A number of comments concerning the notation are also in place. To simplify the notation, we shall drop the dependence on $u$, and write $\widehat{A}$, $A$, $\widehat{M}$, $M$, $\widehat{\xi}_{r}$, $\widehat{\gamma}_{2,i}$, etc. instead of $\widehat{A}(u)$, $A(u)$, $\widehat{M}(u)$, $M(u)$, $\widehat{\xi}_{r}(u)$, $\widehat{\gamma}_{2,i}(u)$, etc. Similarly, $\sup$ will denote $\sup_{u\in{\cal H}}$ and $\xi=O_{p,{\scriptsize\rm sup}}(b_{T})$ will stand for $\sup|\xi|=O_{p}(b_{T})$. As throughout this work, $G_{h}(u)$ will stand for a scaled kernel function $h^{-1}G(h^{-1}u)$. The kernel function $G$ will be normalized to integrated to $1$ when it is important. We shall use both $K$ in Assumption (K) of Appendix A.1 and other functions related to $K$, which will be denoted as

[TABLE]

where $C$ is such that $\|K^{(q)}_{p}\|_{1}=1$. Note that $K^{(0)}_{1}(u)=K(u)$. When $q=0$, we shall simply write $K_{p}(u)$.

Proof of Proposition 4.2: We shall first prove (4.16). The first step consists of showing that the ordered eigenvalues $0\leq\widehat{\gamma}_{2,1}\leq\ldots\leq\widehat{\gamma}_{2,r}$ entering $\widehat{\xi}_{r}$ can be replaced in the asymptotic limit by the ordered eigenvalues $0\leq\widehat{\lambda}_{2,1}\leq\ldots\leq\widehat{\lambda}_{2,r}$ of the $r\times r$ matrix

[TABLE]

where $D$ is a $p\times r$ matrix described below. The key here is that the matrix (A.3) is $r\times r$ so that the sum of its $r$ eigenvalues is just its trace, which is amenable to easier manipulations. The matrix $D$ will also play an important role of standardization.

The $p\times r$ matrix $D$ and another $p\times(p-r)$ matrix $\widetilde{D}$ enter into a $p\times p$ matrix $D_{0}=(\widetilde{D},D)$ characterized as follows: $D_{0}$ consists of the “eigenvectors” associated with the eigenvalues $\gamma_{2,i}$ through the characteristic equation

[TABLE]

that satisfy

[TABLE]

Said differently, $F^{-1}D_{0}$ consists of the eigenvectors of $(FMF)^{2}$, that is,

[TABLE]

Since we deal with squared symmetric matrices, $F^{-1}D_{0}$ also consists of the eigenvectors of $FMF$, whose squared eigenvalues are $\gamma_{2,i}$, so that

[TABLE]

where $(\gamma_{i^{\prime}})^{2}=\gamma_{2,i}$ and the prime indicates that the order is not necessarily that of the increasing order in $\gamma_{1}\leq\ldots\leq\gamma_{p}$. This fact will be used below. Note also that $\gamma_{2,r}=\ldots=\gamma_{2,1}=0$ by assumption and hence that $D^{\prime}M=0$.

The next result will justify the replacement of the eigenvalues $\widehat{\gamma}_{2,i}$.

Lemma A.2**.**

With the above notation and under Assumptions (A.1), (A.2g), we have for $i=1,\ldots,r$,

[TABLE]

Proof: Let $|B|$ denote the determinant of a matrix $B$. As on p. 173 of Robin and Smith (2000), we have

[TABLE]

[TABLE]

By using the relation $|(B_{11}\ B_{12};\ B_{21}\ B_{22})|=|B_{11}|\cdot|B_{22}-B_{21}B_{11}^{-1}B_{12}|$, we have further that

[TABLE]

where

[TABLE]

By Proposition A.1 and the smoothness of $\widetilde{D}$ by Proposition A.2, note that

[TABLE]

where the second equality follows from (A.6) and (A.5). This shows that, asymptotically, $|\widehat{S}|>0$. Hence, in view of (A.9), we may suppose without loss of generality that $a_{T}^{2}\widehat{\gamma}_{2,i}$ are the eigenvalues of the matrix $\widehat{W}\widehat{V}$. The matrix $\widehat{V}$ is symmetric and its eigenvalues are positive asymptotically since $\widehat{V}\to_{p}D^{\prime}F^{-2}D$. Then, $\widehat{V}$ may be assumed to be positive definite, and $a_{T}^{2}\widehat{\gamma}_{2,i}$ be taken as eigenvalues of $\widehat{V}^{1/2}\widehat{W}\widehat{V}^{1/2}$. Since the matrix is symmetric, by applying the Wielandt-Hoffman theorem (e.g. Golub and Van Loan (2012)), we get that

[TABLE]

By using Proposition A.1 and the fact that the square root is a continuous operation on positive definite matrces, $\widehat{V}^{1/2}=(D^{\prime}F^{-2}D)^{1/2}+O_{p,{\scriptsize\rm sup}}((Th/\ln T)^{-1/2})=I_{r}+O_{p,{\scriptsize\rm sup}}((Th/\ln T)^{-1/2})$. Similarly, for $\widehat{W}$, we have $\widehat{W}=a_{T}^{2}D^{\prime}\widehat{M}D^{2}\widehat{M}D+O_{p,{\scriptsize\rm sup}}((Th/\ln T)^{-1/2})=a_{T}^{2}D^{\prime}(\widehat{M}-M)D^{2}(\widehat{M}-M)D+O_{p,{\scriptsize\rm sup}}((Th/\ln T)^{-1/2})$, where the first equality follows from the relation

[TABLE]

The latter relation holds by the following argument. Note that it is equivalent to

[TABLE]

and in view of (A.6), follows from $FM\widetilde{D}\mbox{\rm diag}(\gamma_{2,p},\ldots,\gamma_{2,r+1})^{-1}\widetilde{D}^{\prime}MF=F^{-1}\widetilde{D}^{2}F^{-1}$ or

[TABLE]

which is a consequence of (A.7).   $\Box$

By Lemma A.2, instead of working with $\int\widehat{\xi}_{r}du$, we can focus instead on

[TABLE]

and $\int\widehat{L}_{r}du$. Write $\widehat{M}-M=S_{1}+S_{2}-S_{3}-S_{4}$ as in (A.1) of the proof of Proposition 4.1. Then,

[TABLE]

where

[TABLE]

and $(\pm 1)$ in (A.11) accounts for the signs of $S_{j},S_{k}$ in the decomposition (A.1). The next two lemmas concern the asymptotics of $\widehat{G}_{r,jk}$.

Lemma A.3**.**

Under Assumptions (Y1), (Y2), (K), (A1), (A2g), and

[TABLE]

we have

[TABLE]

As noted following Assumption (A2g), this assumption is needed to apply Proposition A.2 to have the smoothness of $D=D(u)$.

Proof: In view of the definition of $S_{1}$ following (A.1), we can write

[TABLE]

[TABLE]

[TABLE]

We shall argue first that $D=D(u)$ in $I_{1}$ and $I_{2}$ above can be replaced by $D(t/T)$, and shall denote the respective terms by $\widetilde{I}_{1}$ and $\widetilde{I}_{2}$. The relation (A.5) will then be used to simplify $\widetilde{I}_{1}$ and $\widetilde{I}_{2}$. Finally, we will show that $\widetilde{I}_{1}$ produces the centering in (A.13) and $\widetilde{I}_{2}$ yields the asymptotic normality in (A.13).

To see why $I_{1}$ can be replaced by $\widetilde{I}_{1}$, that is, $D=D(u)$ be replaced by $D(t/T)$, consider one of the terms in the difference between $I_{1}$ and $\widetilde{I}_{1}$, namely,

[TABLE]

(The other terms in the difference can be dealt with similarly.) Then, by using the smoothness of $D(u)$ by Proposition A.2, the expression (4.4) for $a_{T}$ and Lemma A.1, $(ii)$, we have

[TABLE]

[TABLE]

[TABLE]

for a random constant $C(w)$. The term $R_{1}$ would then not affect the asymptotics (A.13) since $h/\sqrt{h}\to 0$.

For $I_{2}$ and $\widetilde{I}_{2}$, consider similarly a term from their difference given by

[TABLE]

[TABLE]

We shall use Proposition A.3 to obtain a convergence rate for $R_{2}$. To apply the proposition, a number of (conditional) expectations involving $b_{T,t_{1},t_{2}}(X_{t_{1}}X_{t_{1}}^{\prime},X_{t_{2}}X_{t_{2}}^{\prime})$ need to be evaluated. Note that $\mathbb{E}b_{T,t_{1},t_{2}}(X_{t_{1}}X_{t_{1}}^{\prime},X_{t_{2}}X_{t_{2}}^{\prime})=0$ and $\mathbb{E}(b_{T,t_{1},t_{2}}(X_{t_{1}}X_{t_{1}}^{\prime},X_{t_{2}}X_{t_{2}}^{\prime})|X_{t_{1}}X_{t_{1}}^{\prime})=0$ and similarly when conditioning on $X_{t_{2}}X_{t_{2}}^{\prime}$. We thus only need to consider $\mathbb{E}(b_{T,t_{1},t_{2}}(X_{t_{1}}X_{t_{1}}^{\prime},X_{t_{2}}X_{t_{2}}^{\prime}))^{2}$. By using the generalized Minkowski’s inequality and the smoothness of $D(u)$ by Proposition A.2, note that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we used the definiton of $a_{T}$. It follows from Lemma A.5 that

[TABLE]

Hence, by Proposition A.3, $R_{2}$ is of the order $O_{p}(h^{3/2})$ and hence does not affect the asymptotics (A.13) since $h^{3/2}/\sqrt{h}\to 0$.

We can thus replace $I_{1}$ and $I_{2}$ by $\widetilde{I}_{1}$ and $\widetilde{I}_{2}$, respectively, which in view of (A.5) can be expressed as

[TABLE]

[TABLE]

where $\widetilde{Y}_{t}=D^{{}^{\prime}}(t/T)X_{t}$.

For $\widetilde{I}_{1}$, write it as

[TABLE]

[TABLE]

[TABLE]

It can be checked that

[TABLE]

Hence, by Lemma A.1, $(ii)$,

[TABLE]

For $\widetilde{I}_{1,1}$, setting $\xi_{t}=\sum_{i,j=1}^{r}((\widetilde{Y}_{t}\widetilde{Y}_{t}^{\prime}-I_{r})_{ij}^{2}-\mathbb{E}(\widetilde{Y}_{t}\widetilde{Y}_{t}^{\prime}-I_{r})_{ij}^{2})$, $p=1+\epsilon/4$ and using the von Bahr-Essen inequality (von Bahr and Esseen (1965)), we obtain that

[TABLE]

and hence $\widetilde{I}_{1,1}=O_{p}(T^{1/p-1})$. By assumption, $Th^{3/2}\to\infty$ and $T^{1-1/p}h^{1/2}=T^{\epsilon/(4+\epsilon)}h^{1/2}\to\infty$. This shows that the error term in $\widetilde{I}_{1,1}$ and $\widetilde{I}_{1,2}$ do not affect the asymptotics (A.13), and that $\widetilde{I}_{1,2}$ produces the desired asymptotic mean.

For $\widetilde{I}_{2}$, consider

[TABLE]

where

[TABLE]

We will argue that $\widetilde{I}_{2}/\sqrt{h}$ is asymptotically normal with the desired limiting variance. By Proposition 3.2 in de Jong (1987), it is enough to show that

$\mbox{Var}(\frac{\widetilde{I}_{2}}{\sqrt{h}})\to\frac{\|\overline{K}\|_{2}^{2}}{\|K\|_{2}^{4}}\frac{2(r\mu_{4}^{2}+2r(r-1))}{\mu_{4}^{2}}|{\cal H}|$; 2. 2.

$G_{T,i}=o(1)$, $i=1,2,4$, where

[TABLE]

The first point above follows from

[TABLE]

the observation that

[TABLE]

and Lemma A.5. For the second point above, $G_{T,1}$ can be bounded up to a constant by

[TABLE]

For example, the first term in the sum of $G_{T,2}$ can be bounded up to a constant by

[TABLE]

For example, the first term in the sum of $G_{T,4}$ can be bounded up to a constant by

[TABLE]

[TABLE]

The rate $o(1)$ for each of these bounds follows from Lemma A.5. This completes the proof of the lemma.   $\Box$

Lemma A.4**.**

Under Assumptions (Y1), (Y2), (K), (A1), (A2g), and

[TABLE]

we have, for $(j,k)\neq(1,1)$,

[TABLE]

In fact, the proof of lemma establishes sharper rates of convergence of $\widehat{G}_{r,jk}$ to 0. The rate given in the lemma is what is needed to conclude the convergence (4.16). Indeed, the latter now follows immediately from the arguments above and, in particular, Lemmas A.3 and A.4 .

Proof: We consider only the cases $(j,k)=(2,2)$, $(3,3)$, $(4,4)$, $(1,2)$, $(1,3)$ and $(1,4)$. The other mixed cases can be dealt with similarly.

For $(j,k)=(2,2)$, we have $\int\widehat{L}_{r,22}du=a_{T}^{2}\int\mbox{tr}\{D^{\prime}S_{2}D^{2}S_{2}D\}du$ where $S_{2}$ is defined in (A.2). As in the proof of Proposition 4.1, $S_{2}=O(h+(Th)^{-1})$ uniformly in $u$, where the latter follows by using Lemma A.1, $(ii)$. Then, by using the smoothness of $D=D(u)$ by Proposition A.2, we have $\int\widehat{L}_{r,22}du=O(Th(h+(Th)^{-1})^{2})=O(Th^{3}+(Th)^{-1}+h)$. This is of the order $o(h^{1/2})$ since, in particular, $Th^{5/2}\to 0$ by assumption. For $(j,k)=(3,3)$, we have similarly $\int\widehat{L}_{r,33}du=a_{T}^{2}\int\mbox{tr}\{D^{\prime}S_{3}D^{2}S_{3}D\}du$ where $S_{3}$ is defined in (A.2). The term $S_{3}$ does not depend on $u$ and its rate is $O_{p}(T^{-1/2})$ as shown in the proof of Proposition 4.1. This leads to $\int\widehat{L}_{r,33}du=O_{p}(Th(T^{-1/2})^{2})=O_{p}(h)$, which is again of the order $o_{p}(h^{1/2})$ as desired. For $(j,k)=(4,4)$, $S_{4}$ does not depend on $u$ either and its rate is $O(T^{-1})$ as shown in the proof of Proposition 4.1. This leads to $\int\widehat{L}_{r,44}du=O(Th(T^{-1})^{2})=O(h/T)=o(h^{1/2})$.

For $(j,k)=(1,2)$, we need to consider $\int\widehat{L}_{r,12}du=a_{T}^{2}\int\mbox{tr}\{D^{\prime}S_{1}D^{2}S_{2}D\}du$. After matrix multiplication and taking the trace, a general term in $\int\widehat{L}_{r,12}du$ has the form

[TABLE]

where $\ast$ refers to an index pair that can change from matrix to matrix. Furthermore, in view of (A.2),

[TABLE]

By using the facts that $S_{2}=O(h+(Th)^{-1})$ uniformly in $u$ as above and $\int_{\mathbb{R}}K_{h}(u-t/T)du=1$, it follows that $\mathbb{E}R_{1,2}^{2}=O(a_{T}^{4}T^{-1}(h+(Th)^{-1})^{2})=O(Th^{4}+h^{2}+1/T)$ and hence that $R_{1,2}=o_{p}(h^{1/2})$ since, in particular, $Th^{3}\to 0$. For $(j,k)=(1,3)$, a general term of interest is similarly,

[TABLE]

and hence $\mathbb{E}R_{1,3}^{2}=a_{T}^{4}T^{-2}=h^{4}$. This leads to $R_{1,3}=o_{p}(h^{1/2})$. The case $(j,k)=(1,4)$ can be dealt with similarly by using the fact that $S_{4}=O(T^{-1})$.   $\Box$

Finally, we prove the last statement of Proposition 4.2 concerning the behavior of the test statistic under the alternative. For this, note that by Proposition A.1, $\sup_{u\in{\cal H}}|a_{T}^{-2}\xi_{r}(u)-\sum_{i=1}^{r}\gamma_{2,i}(u)|\to_{p}0$. Furthermore, under the considered alternative $H_{1}$, and by using smoothness of $\gamma_{2,i}(u)$, we have $\sup_{u\in{\cal H}}\sum_{i=1}^{r}\gamma_{2,i}(u)>0$. It follows that $\widehat{\xi}_{r}=(a_{T}^{2}\int_{\cal H}a_{T}^{-2}\xi_{r}(u)du-C_{1})/(C_{2}h^{1/2})$ with constants $C_{1}$, $C_{2}$ behaves asymptotically as $Ca_{T}^{2}/h^{1/2}=C^{\prime}Th^{1/2}\rightarrow\infty$ (for example, since Lemma A.3 assumes $Th^{3/2}\rightarrow\infty$). This concludes the proof of Proposition 4.2.

A.3 Auxiliary technical results for proof of Proposition 4.2

The following auxiliary results were used in the proof of Proposition 4.2. The first result is analogous to the results in e.g. Newey (1994), Lemma B.1. A separate proof though is needed since averaging and smoothing in the estimators here are with respect to time $t$, rather than the values of a random variable as in the aforementioned results.

Proposition A.1**.**

Under Assumptions (Y1), (Y2), (K), (A1), (A2g), and

[TABLE]

we have,

[TABLE]

and

[TABLE]

Proof: We outline the proof in the case $p=1$. Write

[TABLE]

where

[TABLE]

Set $\delta_{T}=(Th/\ln T)^{-1/2}$ or $\delta^{-1}=(Th/\ln T)^{1/2}$. We will show that $\delta_{T}^{-1}\sup_{u\in{\cal H}}|R_{i}(u)|=O_{p}(1)$, $i=1,2,3,4$. This will prove the convergence (A.15) for $k=1$, and the convergence for $k=-1$ will also follow since $A^{2}(u)$ is smooth for $u\in{\cal H}$.

By Lemma A.1, $(ii)$, we have $\sup_{u\in{\cal H}}|R_{2}(u)|=O((Th)^{-1})$ and hence also $\delta_{T}^{-1}\sup_{u\in{\cal H}}|R_{2}(u)|\to 0$ since $Th\to\infty$. Similarly, by the smoothness of $A^{2}(u)$ and Lemma A.1, $(ii)$, applied to $G=K_{1}$,

[TABLE]

and hence $\delta_{T}^{-1}\sup_{u\in{\cal H}}|R_{4}(u)|\to 0$ since $Th^{3}\to 0$. The terms $R_{1}(u)$ and $R_{3}(u)$ are a more delicate to deal with. In their analysis, we will follow the proof of Lemma B.1 in Newey (1994) which, unsurprisingly, involves Bernstein’s inequality.

We will consider the term $R_{1}(u)$ only, since the term $R_{3}(u)$ can be dealt with similarly. For the term $R_{1}(u)$, note that

[TABLE]

where we used the Lipschitz continuity of $K$ and the law of large numbers in the last step with a constant $C^{\prime}(\omega)$ which does not depend on $T,h,u,v$. The difference $\delta_{T}^{-1}|R_{1}(u)-R_{1}(v)|$ can then be made small as long as, for example, $|u-v|\leq T^{-2}$. Now, cover ${\cal H}$ with balls of radius $T^{-2}$ and centered at $u_{j}$, where we can take $j=1,\ldots,cT^{2}$. Then, for large $T$, any arbitrarily small $\epsilon$ and any fixed $\eta>0$,

[TABLE]

[TABLE]

Assume first that $Y_{t}^{2}-1$, $t\in\mathbb{Z}$, are bounded by a constant, almost surely. Then, by Bernstein’s inequality, the last sum in (A.18) is bounded by (up to a constant)

[TABLE]

By Lemma A.1, $(ii)$, $\frac{1}{Th}\sum_{t=1}^{T}K^{2}(u_{j}-t/T)/h)$ behaves as a positive constant uniformly over $u_{j}$. Since $(Th\ln T)^{1/2}/(Th)\to 0$ by assumption, the above bound becomes, for large enough $T$ and constants $C_{1},C_{2}>0$,

[TABLE]

The latter converges to [math] as long as $\eta$ is large enough (so that $C_{2}\eta^{2}-2>0$). When the variables $\widetilde{\epsilon}_{t}:=Y_{t}^{2}-1$ are not bounded almost surely, a standard truncation argument is used. Let $\widetilde{\epsilon}_{t,T}=Y_{t}^{2}-1$ if $|Y_{t}^{2}-1|\leq C_{0}T^{1/p}$ and $\widetilde{\epsilon}_{t,T}=C_{0}T^{1/p}$ otherwise, where $p=2+\epsilon/2$ and $\epsilon$ appears in Assumption (Y2). Let also $\widetilde{R}_{1}(u)=\frac{1}{T}\sum_{t=1}^{T}\widetilde{\epsilon}_{t,T}K_{h}(u-t/T)$. Suppose for simplicity that $\mathbb{E}\widetilde{\epsilon}_{t,T}=0$. Then,

[TABLE]

[TABLE]

which can be made arbitrarily small by taking large enough $C_{0}$. The above argument for $R_{1}(u)$ can now be applied to $\widetilde{R}_{1}(u)$ for this large fixed $C_{0}$, but keeping track of and dealing with the truncation $C_{0}T^{1/p}$ which depends on $T$. More precisely, the bound (A.19) becomes

[TABLE]

and the same argument as above applies as long as $T^{1/p}(Th\ln T)^{1/2}/(Th)\to 0$. The latter follows from the assumption $T^{1-2/p}h/\ln T=T^{\epsilon/(4+\epsilon)}h/\ln T\to\infty$.

The convergence (A.16) follows from (A.15) and the fact that $\widehat{\overline{A}^{2}}=\frac{1}{T}\sum_{t=1}^{T}X_{t}X_{t}^{\prime}$ converges to $\overline{A}^{2}=\int_{0}^{1}A^{2}(u)du$ at a rate faster than $\delta_{T}$. Indeed, by using the notation of the proof of Proposition 4.1, this difference is $S_{3}+S_{4}$ which is of the order $1/T^{1/2}$ as shown in that proof.   $\Box$

The next result allows having smooth eigenvectors associated with the matrix $M$. This is the only place where the analyticity of the matrix $M$ needs to be assumed.

Proposition A.2**.**

Under Assumption (A2g), the matrix $D_{0}$ in the proof of Proposition 4.2 can be chosen analytic.

Proof: By Assumption (A2g), the matrix $FMF$ is analytic. By using the analytic singular value decomposition (e.g. Bunse-Gerstner et al. (1991)), there are $p\times p$ analytic matrices $U$ and $T$ such that $FMF=UTU^{\prime}$, where $T=\mbox{diag}(t_{1},\ldots,t_{p})$ are the singular values and an orthogonal matrix $U$ consists of the eigenvectors of $FMF$. Now take $D_{0}=FU$. Then, $D_{0}$ is analytic and satisfies (A.5).   $\Box$

We also used the following result on several occasions concerning the limiting behavior of a quadratic form

[TABLE]

where $W_{t}$, $t=1,\dots,T$, are i.i.d. random vectors in $\mathbb{R}^{d}$ and $b_{T,t_{1},t_{2}}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$.

Proposition A.3**.**

Let $Q_{T}$ be a quadratic form defined by (A.21) and assume that $\mathbb{E}b_{T,t_{1},t_{2}}(W_{t_{1}},W_{t_{2}})^{2}<\infty$. Then,

[TABLE]

When $b_{T,t_{1},t_{2}}(\cdot)\equiv b_{T}(\cdot)$ and the quadratic form $Q_{T}$ becomes a $U$-statistic, the result restates Lemma C.1 in Fortuna (2008), p. 181. As for that lemma, the proof below uses the arguments of the proof of Lemma 3.1 in Powell et al. (1989).

Proof: Let

[TABLE]

[TABLE]

Then,

[TABLE]

where

[TABLE]

This yields, by using the independence of $W_{t}$’s,

[TABLE]

since $\mathbb{E}c_{T,t_{1},t_{2}}(W_{t_{1}},W_{t_{2}})c_{T,t_{1}^{\prime},t_{2}^{\prime}}(W_{t_{1}^{\prime}},W_{t_{2}^{\prime}})=0$ whenever $(t_{1},t_{2})\neq(t_{1}^{\prime},t_{2}^{\prime})$ (the interested reader should check this, especially when e.g. $t_{1}=t_{1}^{\prime}$, $t_{2}\neq t_{2}^{\prime}$). Since $\mathbb{E}c_{T,t_{1},t_{2}}(W_{t_{1}},W_{t_{2}})^{2}=O(\mathbb{E}b_{T,t_{1},t_{2}}(W_{t_{1}},W_{t_{2}})^{2})$, we obtain that

[TABLE]

which is the last term in the relation (A.22). The first term on the right-hand side of the relation (A.22) is the first sum in the definition of $\widehat{Q}_{T}$. Finally, to deal with the second sum in the definition of $\widehat{Q}_{T}$, note that, for example,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where for the first inequality above we used the generalized Minkowski inequality. This yields the second term on the right-hand side of the relation (A.22). The third term results from a similar argument for the term in the definition of $\widehat{Q}_{T}$, where the conditioning is on $W_{t_{2}}$. This proves the result (A.22).   $\Box$

Finally, the following auxiliary result concerns the behavior of the sum of power integrals of kernel functions, and was used in the proof of Lemma A.3. We also note that the result extends easily to other powers than $2$ and $4$ appearing below, and a product of two different kernel functions. It is formulated only for what is needed in Lemma A.3 for the shortness sake.

Lemma A.5**.**

For two kernel functions $K,\widetilde{K}$ satisfying Assumption (K), let

[TABLE]

Then, as $T\to\infty$, $h\to 0$, $Th\to\infty$, we have

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

Proof: We first show (A.23) which requires a more delicate treatment. Denote the left hand-side of (A.23) by $I_{0}$ and the interval ${\cal H}$ by $[a,b]$. Suppose $\mbox{supp}\{K\}=[-S,S]$ as in Assumption (K). Then,

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and otherwise all the sums for $t_{1},t_{2}$ are over $1\leq t_{1}<t_{2}\leq T$. Note that $I_{0,4}=0$ for small enough $h$, since $Th\to\infty$. For $I_{0,1}$, we have

[TABLE]

[TABLE]

Note that

[TABLE]

and, after changing the summation indices,

[TABLE]

[TABLE]

Thus, $I_{0,1}\to(b-a)\|\overline{K}\|_{2}^{2}/2$. For $I_{0,2}$, we have

[TABLE]

[TABLE]

[TABLE]

One can show similarly that $I_{0,3}\to 0$.

The relations (A.24)–(A.26) are easier to obtain because ${\cal H}$ can be replaced by $\mathbb{R}$ and the involved integrals simplified to convolutions. We will show how this can be argued for (A.25). Denote the sum on the left-hand side of (A.25) by $I_{2}$. Then,

[TABLE]

[TABLE]

[TABLE]

The relation (A.24) and (A.26) can be argued for similarly.   $\Box$