Segmentation of High-Dimensional Matrix-Variate Time Series

Abstract

In this chapter, we introduce a new segmentation method for high-dimensional matrix-variate time series. Specifically, we look for linear transformations to segment the matrix into many small sub-matrices for which each of them is uncorrelated with the others both contemporaneously and serially, thus they can be analyzed separately, which will greatly reduce the number of parameters to be estimated in terms of modeling. To overcome the identification issue, we propose a two-step and more structured procedure to segment the rows and columns separately. When the dimension is large in relation to the sample size, we assume the transformation matrices are sparse and use threshold estimators for the (auto) covariance matrices. Unlike principal component analysis (PCA) for independent data, we cannot guarantee that the required linear transformation exists. When it does not, the proposed method provides an approximate segmentation, which may be useful for forecasting. The proposed method is illustrated with simulated data examples.

1. Introduction

Modern scientific studies often gather data under combinations of multiple factors. For example, neuroimaging experiments record brain activity at multiple spatial locations, at multiple time points, and under a variety of experimental stimuli. Studies of social networks record social links of a variety of types from multiple initiators of social activity to multiple receivers of the activity. Data such as these are naturally represented not as lists or tables of numbers, but as multi-indexed arrays or tensors. As many types of such data are collected over time, it is natural to view them as tensor-valued time series. The matrix-valued time series is a sequence of second-order random tensors. For example, financial and economic studies often collect data from different countries with a number of economic indicators (e.g. GDP growth, unemployment rate, etc.) every quarter. Therefore, it is important and interesting to develop appropriate statistical methods to analyze such types of data. The most common approach to modeling such data is to stack the matrix into a large vector, and then apply the standard multivariate methods. However, such an approach will ignore the matrix structure of the data, which can lead to inefficient use of data and important patterns in the data being overlooked. For example, [1] pointed out that after vectorizing the matrices the resulting vectors have a Kronecker structure. Ignoring this structure means that a much larger number of parameters need to be estimated. Therefore, it is urgent to find an effective way to reduce the number of parameters, especially when the dimension is large.

When the data are independent and identically distributed (i.i.d.), [2] introduced dimension folding sufficient reduction for conditional mean functions, [3] proposed a dimension folding method for data with matrix-valued predictors, [4, 5, 6] extended the generalized linear models to the matrix- and tensor-valued predictors for analyzing image data. Ding and Dennis [7] studied the matrix-variate regression with a matrix-variate response. An incomplete list of publications also includes [8, 9, 10, 11, 12]. With temporal dependence, the matrix-valued time series has not been well studied in the literature; [13] handled this kind of data in signal and image processing, and [14, 15] proposed two versions of factor models for matrix-valued time series, which maintain and utilize the matrix structure to achieve the dimension reduction.

In this chapter, we extend the PCA approach of [16] to matrix-variate time series without stacking the matrix into a vector, and the structure can be preserved. Our goal is as follows: Let ${\mathbf{Y}}_t=\left(y_{\mathit{ij}}^t\right)$ be a $p\times q$ matrix-variate time series, that is. there are $\mathit{pq}$ recorded values at each time, for example, $p$ individuals and over $q$ indices or variables. We assume ${\mathbf{Y}}_t$ can be represented as

$${\mathbf{Y}}_t={\mathbf{BW}}_t{\mathbf{A}}^{\mathrm{T}},$$E1

where $\mathbf{B}\in R^{p\times p}$, $\mathbf{A}\in R^{q\times q}$, and ${\mathbf{W}}_t$ is a latent $p\times q$ matrix in which the rows are divided into $p_1\left(\le p\right)$ groups, and there are no correlations across different groups at all time lags, and the columns are divided into $q_1\left(\le q\right)$ groups, and there are no correlations across different groups at all time lags either. With such a decomposition, we only need to model the small sub-matrices in $\mathbf{W}$ separately, and we can achieve substantial dimension reduction. As $\mathbf{B}$, ${\mathbf{W}}_t$, and $\mathbf{A}$ are all latent ones and the identification is a big issue. For example, even when ${\mathbf{W}}_t$ is observable, $\left(\mathbf{A},\mathbf{B}\right)$ can be replaced by $\left(\mathbf{A}/c,c\mathbf{B}\right)$ for any nonzero constant $c$ without changing the relationship of (1).

Instead of estimating them simultaneously, we propose in this paper a two-step and more structured approach; first, we seek a column transformation, that is. We transform linearly the columns of ${\mathbf{Y}}_t$ into $q$ new variables, and ideally, those $q$ new variables form $q_1$ uncorrelated groups with $q_1\le q$. The second step applies the same segmentation method to the $p$ rows of the obtained ones in the first step, and the transformation of the rows will not alter the uncorrelatedness of the column groups in the first step. In the end, this new matrix can be divided into several smaller sub-matrices, and those sub-matrices are uncorrelated with each other both contemporaneously and serially. Our method is a building block for modeling tensor-valued time series, and it turns out that all tensor-valued time series can be rearranged as a matrix time series by matricization, see [17].

2. Methodology

2.1 Setting and method

Let ${\mathbf{Y}}_t=\left({\mathbf{y}}_1^t,\dots ,{\mathbf{y}}_q^t\right)$ be an observable $p\times q$ matrix-valued time series with ${\mathbf{y}}_i^t\in R^p$. We assume ${\mathbf{Y}}_t$ admits a latent segmentation structure:

$${\mathbf{Y}}_t={\mathbf{X}}_t{\mathbf{A}}^{\mathrm{T}},$$E2

where ${\mathbf{X}}_t$ is an unobservable $p\times q$ matrix-valued time series in which the $q$ columns can be classified into $q_1$ ($>1$) groups, and any two groups are contemporaneously and serially uncorrelated, and $\mathbf{A}\in R^{q\times q}$ is an unknown constant matrix. Before we proceed further, we give the definitions of row- and column-covariance matrix between two random matrices.

Definition 1 Let ${\mathbf{U}}_t\in R^{s_1\times r_1}$ and ${\mathbf{V}}_t\in R^{s_2\times r_2}$. If $r_1=r_2=r$, the covariance matrix over the columns between ${\mathbf{U}}_t$ and ${\mathbf{V}}_t$ is defined as

$${\mathrm{Cov}}_c\left({\mathbf{U}}_t,{\mathbf{V}}_t\right)≔\frac{1}{r}E\left({\mathbf{U}}_t-E{\mathbf{U}}_t\right){\left({\mathbf{V}}_t-E{\mathbf{V}}_t\right)}^{\mathrm{T}},$$E3

and if $s_1=s_2=s$, the covariance matrix over the rows is defined as

$${\mathrm{Cov}}_r\left({\mathbf{U}}_t,{\mathbf{V}}_t\right)≔\frac{1}{s}E{\left({\mathbf{U}}_t-E{\mathbf{U}}_t\right)}^{\mathrm{T}}\left({\mathbf{V}}_t-E{\mathbf{V}}_t\right)={\mathrm{Cov}}_c\left({\mathbf{U}}_t^{\mathrm{T}},{\mathbf{V}}_t^{\mathrm{T}}\right).$$E4

The variances ${\mathrm{Var}}_c\left({\mathbf{U}}_t\right)$ and ${\mathrm{Var}}_r\left({\mathbf{U}}_t\right)$ can be defined in a similar way. In particular, when $r=1$ or $s=1$, (3) or (4) reduces to the traditional case for two random vectors.

In model (2), we assume ${\mathbf{Y}}_t$ and ${\mathbf{X}}_t$ are both weakly stationary in the sense that the means and the autocovariances do not vary with respect to time for any fixed $\left(p,q\right)$. The stationarity of ${\mathbf{Y}}_t$ can be inherited from ${\mathbf{X}}_t$ through (2), and a sufficient condition for this is to assume $\mathrm{Vec}\left({\mathbf{X}}_t\right)$ and $\mathrm{Vec}\left({\mathbf{Y}}_t\right)$ are stationary, where $\mathrm{Vec}\left(\cdot \right)$ is the vectorization of a matrix. Denote the segmentation of ${\mathbf{X}}_t$ by

$${\mathbf{X}}_t=\left({\mathbf{x}}_1^t,\dots ,{\mathbf{x}}_q^t\right)=\left({\mathbf{X}}_t^{\left(1\right)},\dots ,{\mathbf{X}}_t^{\left(q_1\right)}\right)$$E5

with ${\mathrm{Cov}}_r\left({\mathbf{X}}_t^{\left(i\right)},{\mathbf{X}}_s^{\left(j\right)}\right)=0$ for all $t,s$ and $i\ne j$. Therefore, all the autocovariances of ${\mathbf{X}}_t^{\mathrm{T}}$ are of the same block-diagonal structure with $q_1$ blocks, and ${\mathbf{X}}_t^{\left(1\right)},\dots ,{\mathbf{X}}_t^{\left(q_1\right)}$ can be modeled or forecasted separately as far as their linear dynamic structure is concerned.

Now, we spell out how to find the segmentation transformation under (2) and (5). Without loss of generality, we assume

$${\mathrm{Var}}_r\left({\mathbf{Y}}_t\right)={\mathbf{I}}_q\text{and}{\mathrm{Var}}_r\left({\mathbf{X}}_t\right)={\mathbf{I}}_q.$$E6

The first equation in (6) is implied by replacing ${\mathbf{Y}}_t$ by ${\mathbf{Y}}_t{\widehat{\mathbf{S}}}_{y,0}^{-1/2}$, where ${\widehat{\mathbf{S}}}_{y,0}$ is a consistent estimator of ${\mathrm{Var}}_r\left({\mathbf{Y}}_t\right)$. The second equation is to conceptually replace ${\mathbf{X}}_t$ by ${\mathbf{X}}_t{\widehat{\mathbf{S}}}_{x,0}^{-1/2}$, where ${\widehat{\mathbf{S}}}_{x,0}^{-1/2}$ is a consistent estimator for ${\mathrm{Var}}_r\left({\mathbf{X}}_t\right)$, and it will not alter the fact that there are no correlations across different groups. As both $\mathbf{A}$ and ${\mathbf{X}}_t$ are unobservable, (6) implies that we can view ${\widehat{\mathbf{S}}}_{y,0}^{-1/2}\mathbf{A}{\widehat{\mathbf{S}}}_{x,0}^{1/2}$ as $\mathbf{A}$. As a consequence of (6), the transformation matrix $\mathbf{A}$ in (2) is orthogonal. Let $l_j$ be the number of columns of ${\mathbf{X}}_t^{\left(j\right)}$ with $l_1+\cdots +l_{q_1}=q$. Write $\mathbf{A}=\left({\mathbf{A}}_1,\dots ,{\mathbf{A}}_{q_1}\right)$, where ${\mathbf{A}}_j\in R^{q\times l_j}$. It follows from (2) and (5) that

$${\mathbf{X}}_t^{\left(j\right)}={\mathbf{Y}}_t{\mathbf{A}}_j,j=1,\dots ,q_1.$$E7

However, similar to that in [16], $\mathbf{A}$ and ${\mathbf{X}}_t$ are not uniquely identified in (2), even with additional assumption in (6). For example, let ${\mathbf{H}}_j$ be any $l_j\times l_j$ orthogonal matrix and $\mathbf{H}=diag\left({\mathbf{H}}_1,\dots ,{\mathbf{H}}_{q_1}\right)$. Then $\left(\mathbf{A},{\mathbf{X}}_t\right)$ in (2) can be replaced by $\left(\mathbf{AH},{\mathbf{X}}_t\mathbf{H}\right)$_, while (5) still holds. In fact, only $\mathcal{M}\left({\mathbf{A}}_1\right),\dots ,\mathcal{M}\left({\mathbf{A}}_{q_1}\right)$ are uniquely defined by (2), where $\mathcal{M}\left({\mathbf{A}}_j\right)$ denotes the linear space spanned by the columns of ${\mathbf{A}}_j$. As a result, ${\mathbf{Y}}_t{\mathtt{\Gamma }}_j$ can be taken as ${\mathbf{X}}_t^{\left(j\right)}$ for any $q\times l_j$ matrix ${\mathtt{\Gamma }}_j$ as long as ${\mathtt{\Gamma }}_j^{\mathrm{T}}{\mathtt{\Gamma }}_j={\mathbf{I}}_{l_j}$ and $\mathcal{M}\left({\mathtt{\Gamma }}_j\right)=\mathcal{M}\left({\mathbf{A}}_j\right)$. Thus, to estimate $\mathbf{A}=\left({\mathbf{A}}_1,\dots ,{\mathbf{A}}_{q_1}\right)$, it is sufficient to estimate the linear spaces $\mathcal{M}\left({\mathbf{A}}_1\right),\dots ,\mathcal{M}\left({\mathbf{A}}_{q_1}\right)$.

To discover the latent segmentation, we introduce some notation first. We denote ${\mathbf{y}}_{i:}^t$ and ${\mathbf{x}}_{i:}^t$ the row vectors of ${\mathbf{Y}}_t$ and ${\mathbf{X}}_t$, respectively. For any integer $k$, let ${\mathtt{\Sigma }}_y\left(k\right)={\mathrm{Cov}}_r\left({\mathbf{Y}}_{t+k},{\mathbf{Y}}_t\right)$, ${\mathtt{\Sigma }}_x\left(k\right)={\mathrm{Cov}}_r\left({\mathbf{X}}_{t+k},{\mathbf{X}}_t\right)$, ${\mathtt{\Sigma }}_{y,i,j}\left(k\right)=\mathrm{Cov}\left({\mathbf{y}}_{i:}^{t+k},{\mathbf{y}}_{j:}^t\right)$, and ${\mathtt{\Sigma }}_{x,i,j}\left(k\right)=\mathrm{Cov}\left({\mathbf{x}}_{i:}^{t+k},{\mathbf{x}}_{j:}^t\right)$. By (6), we have ${\mathtt{\Sigma }}_y\left(0\right)={\mathtt{\Sigma }}_x\left(0\right)={\mathbf{I}}_q$. For a pre-specified integer $k_0$, define

$${\mathbf{W}}_y=\sum _{k=0}^{k_0}{\mathtt{\Sigma }}_y\left(k\right){\mathtt{\Sigma }}_y{\left(k\right)}^{\mathrm{T}}={\mathbf{I}}_q+\sum _{k=1}^{k_0}{\mathtt{\Sigma }}_y\left(k\right){\mathtt{\Sigma }}_y{\left(k\right)}^{\mathrm{T}}$$E8

and

$${\mathbf{W}}_x=\sum _{k=0}^{k_0}{\mathtt{\Sigma }}_x\left(k\right){\mathtt{\Sigma }}_x{\left(k\right)}^{\mathrm{T}}={\mathbf{I}}_q+\sum _{k=1}^{k_0}{\mathtt{\Sigma }}_x\left(k\right){\mathtt{\Sigma }}_x{\left(k\right)}^{\mathrm{T}}.$$E9

It follows from (2) and (5) that both ${\mathtt{\Sigma }}_x\left(k\right)$ and ${\mathbf{W}}_x$ are block-diagonal and

$${\mathbf{W}}_y={\mathbf{AW}}_x{\mathbf{A}}^{\mathrm{T}}.$$E10

Note that both ${\mathbf{W}}_y$ and ${\mathbf{W}}_x$ are positive definite matrices; therefore, we have the following decomposition

$${\mathbf{W}}_x{\mathtt{\Gamma }}_x={\mathtt{\Gamma }}_x\mathbf{D},$$E11

where ${\mathtt{\Gamma }}_x$ is a $q\times q$ orthogonal matrix with the columns being the orthonormal eigenvectors of ${\mathbf{W}}_x$, and $\mathbf{D}$ is a diagonal matrix with the corresponding eigenvalues as the elements on the main diagonal. By (10) and (11), ${\mathbf{W}}_y\mathbf{A}{\mathtt{\Gamma }}_x=\mathbf{A}{\mathtt{\Gamma }}_x\mathbf{D}$ and hence the columns of ${\mathtt{\Gamma }}_y≔\mathbf{A}{\mathtt{\Gamma }}_x$ are the orthonormal eigenvectors of ${\mathbf{W}}_y$. Consequently,

$${\mathbf{Y}}_t{\mathtt{\Gamma }}_y={\mathbf{Y}}_t\mathbf{A}{\mathtt{\Gamma }}_x={\mathbf{X}}_t{\mathtt{\Gamma }}_x,$$E12

where the last equality follows from (2). Let

$${\mathbf{W}}_x=diag\left({\mathbf{W}}_{x,1},\dots ,{\mathbf{W}}_{x,q_1}\right),$$E13

where ${\mathbf{W}}_{x,j}$ is an $l_j\times l_j$ positive definite matrix, and the eigenvalues of ${\mathbf{W}}_{x,j}$ are also the eigenvalues of ${\mathbf{W}}_x$. Suppose that ${\mathbf{W}}_{x,i}$ and ${\mathbf{W}}_{x,j}$ do not share the same eigenvalues for any $i\ne j$. Then if we line up the eigenvalues of ${\mathbf{W}}_x$ (i.e. the eigenvalues of ${\mathbf{W}}_{x,1},\dots ,{\mathbf{W}}_{x,q_1}$ combining together) in the main diagonal of $\mathbf{D}$ according to the order of the blocks in ${\mathbf{W}}_x$, ${\mathtt{\Gamma }}_x$ must be a block-diagonal orthogonal matrix of the same shape as ${\mathbf{W}}_x$; see Proposition 1 (i). However, the order of the eigenvalues is latent, and any ${\mathtt{\Gamma }}_x$ defined by (11) is nevertheless a column-permutation (i.e. a matrix consisting of the same column vectors but arranged in a different order) of such a block-diagonal orthogonal matrix; see Proposition 1(ii). By Proposition 1(i), write ${\mathtt{\Gamma }}_x=\mathit{diag}\left({\mathtt{\Gamma }}_{x,1},\dots ,{\mathtt{\Gamma }}_{x,q_1}\right)$, it follows from (5) and (12) that

$${\mathbf{Y}}_t{\mathtt{\Gamma }}_y={\mathbf{X}}_t{\mathtt{\Gamma }}_x=\left({\mathbf{X}}_t^{\left(1\right)}{\mathtt{\Gamma }}_{x,1},\dots ,{\mathbf{X}}_t^{\left(q_1\right)}{\mathtt{\Gamma }}_{x,q_1}\right),$$E14

Hence, ${\mathbf{X}}_t{\mathtt{\Gamma }}_x$ does not alter the fact that there are no correlations between different groups in ${\mathbf{X}}_t$, and ${\mathtt{\Gamma }}_y$ can be regarded as $\mathbf{A}$ so long as the eigenvalues of ${\mathbf{W}}_x$ are ordered appropriately. As they are all latent, ${\mathbf{Y}}_t{\mathtt{\Gamma }}_y$ can be taken as a permutation of ${\mathbf{X}}_t$, and ${\mathtt{\Gamma }}_y$ can be viewed as a column-permutation of $\mathbf{A}$; see the discussion below (7). This leads to the following three-step estimation for $\mathbf{A}$ and ${\mathbf{X}}_t$:

Step 1. Let${\widehat{\mathbf{S}}}_{y,0}$be a consistent estimator for${\mathrm{Var}}_r\left({\mathbf{Y}}_t\right)$. Replace${\mathbf{Y}}_t$by${\mathbf{Y}}_t{\widehat{\mathbf{S}}}_{y,0}^{-1/2}$.

Step 2. Let$\widehat{\mathbf{S}}$be a consistent estimator for${\mathbf{W}}_y$. Calculate an$q\times q$orthogonal matrix${\widehat{\mathtt{\Gamma }}}_y$with columns being the orthonormal eigenvectors of$\widehat{\mathbf{S}}$.

Step 3. The columns of$\widehat{\mathbf{A}}=\left({\widehat{\mathbf{A}}}_1,\dots ,{\widehat{\mathbf{A}}}_{q_1}\right)$are a permutation of the columns of${\widehat{\mathtt{\Gamma }}}_y$such that${\widehat{\mathbf{X}}}_t={\mathbf{Y}}_t\widehat{\mathbf{A}}$is segmented into$q_1$uncorrelated sub-matrix series${\widehat{\mathbf{X}}}_t^{\left(j\right)}={\mathbf{Y}}_t{\widehat{\mathbf{A}}}_j$,$j=1,\dots ,q_1$.

In Steps 1 and 2, the estimators ${\widehat{\mathbf{S}}}_{y,0}$ and $\widehat{\mathbf{S}}$ should be consistent and will be constructed under various scenarios in Section 3 below. The permutation in Step 3 can be carried out by grouping the columns of ${\widehat{\mathbf{Z}}}_t≔{\mathbf{Y}}_t{\widehat{\mathtt{\Gamma }}}_y$.

We now state a proposition that demonstrates the assertion after (13); the proof is similar to Proposition 1 in [16], and we therefore omit it.

Proposition 1: (i) The orthogonal matrix ${\mathtt{\Gamma }}_x$ in (11) can be taken as a block-diagonal orthogonal matrix with the same block structure as ${\mathbf{W}}_x$. (ii) An orthogonal matrix ${\mathtt{\Gamma }}_x$ satisfied (11) if and only if its columns are a permutation of the columns of a block-diagonal orthogonal matrix described in (i), provided that any two different blocks ${\mathbf{W}}_{x,i}$ and ${\mathbf{W}}_{x,j}$ do not share the same eigenvalues.

From Proposition 1, we can see that the proposed method will not be able to separate ${\mathbf{X}}_t^{\left(i\right)}$ and ${\mathbf{X}}_t^{\left(j\right)}$ if ${\mathbf{W}}_{x,i}$ and ${\mathbf{W}}_{x,j}$ share one or more common eigenvalues. But, it does not rule out the possibility that each block ${\mathbf{W}}_{x,j}$ may have multiple eigenvalues.

2.2 Permutation

2.2.1 Permutation rule

According to the discussion in Section 2.1, $\widehat{\mathbf{A}}$ is a permutation of the columns of ${\widehat{\mathtt{\Gamma }}}_y$_, and the permutation can be carried out by grouping the columns of ${\widehat{\mathbf{Z}}}_t≔{\mathbf{Y}}_t{\widehat{\mathtt{\Gamma }}}_y$ into $q_1$ groups, where $q_1$ and the number of columns $l_j$ ($1\le j\le q_1$) are unknown. Let ${\widehat{\mathbf{Z}}}_t=\left({\widehat{\mathbf{z}}}_1^t,\dots ,{\widehat{\mathbf{z}}}_q^t\right)$, ${\mathbf{Z}}_t={\mathbf{Y}}_t{\mathrm{\Gamma }}_y=\left({\mathbf{z}}_1^t,\dots ,{\mathbf{z}}_q^t\right)$, and ${\mathtt{\Gamma }}_{i,j}\left(h\right)$ denote the covariance matrix between two series ${\widehat{\mathbf{z}}}_i^t$ and ${\widehat{\mathbf{z}}}_j^t$ at lag $h$, that is. ${\mathtt{\Gamma }}_{i,j}\left(h\right)=\text{Corr}\left({\widehat{\mathbf{z}}}_i^{t+h},{\widehat{\mathbf{z}}}_j^t\right)$. We say that ${\widehat{\mathbf{z}}}_i^t$ and ${\widehat{\mathbf{z}}}_j^t$ are connected if the multiple null hypothesis

$$H_0:{\mathtt{\Gamma }}_{i,j}\left(h\right)=\mathbf{0}\text{for}\mathrm{any}h=0,\pm 1,\pm 2,\dots ,\pm m$$E15

Is rejected, where $m\ge 1$ is a prescribed integer. We should mention that the true ${\mathtt{\Gamma }}_{i,j}\left(h\right)$ is not known since ${\widehat{\mathbf{z}}}_i^t$ is also one estimator for ${\mathbf{z}}_i^t$, but it will be asymptotically equivalent to $\text{Corr}\left({\mathbf{z}}_i^{t+h},{\mathbf{z}}_j^t\right)$ so long as ${\widehat{\mathtt{\Gamma }}}_y$ is consistent to ${\mathtt{\Gamma }}_y$. Given the structure of ${\mathbf{W}}_x$, this can be done under some regularity conditions, and therefore we also denote the true ${\mathtt{\Gamma }}_{i,j}\left(h\right)=\text{Corr}\left({\mathbf{z}}_i^{t+h},{\mathbf{z}}_j^t\right)$, and the estimator ${\widehat{\mathtt{\Gamma }}}_{i,j}\left(h\right)=\widehat{\text{Corr}}\left({\widehat{\mathbf{z}}}_i^{t+h},{\widehat{\mathbf{z}}}_j^t\right)$. The permutation in Step 3 in Section 2.1 can be performed as follows.

1. Start with the$q$groups with each group containing one column of${\widehat{\mathbf{Z}}}_t$only.

2. Combine two groups together if one connected pair is found.

3. Repeat Step ii above until all connected pairs are within one group.

We introduce below one way to identify the connected pairs of the transformed matrix ${\widehat{\mathbf{Z}}}_t$.

2.2.2 Maximum cross-correlation method

Similar to [16], one natural way to test hypothesis $H_0$ in (15) is to use the maximum cross-correlation over all elements of ${\mathtt{\Gamma }}_{i,j}\left(h\right)$ and all the lags between $-m$ to $m$:

$${\widehat{L}}_n\left(i,j\right)=\underset{\mid h\mid \le m}{max}{\left|{\widehat{\mathtt{\Gamma }}}_{i,j}\left(h\right)\right|}_{\infty },$$E16

where ${\widehat{\mathtt{\Gamma }}}_{i,j}\left(h\right)$ is a sample correlation matrix between ${\widehat{\mathbf{z}}}_i^t$ and ${\widehat{\mathbf{z}}}_j^t$ at lag $h$ when the dimension $p$ and $q$ are fixed, and it is a block-wisely thresholded sample correlation matrix when $p$ and $q$ are moderately high, and it will be constructed in the same way as that in [18]. We would reject $H_0$ for the pair $\left({\widehat{\mathbf{z}}}_i^t,{\widehat{\mathbf{z}}}_j^t\right)$ if ${\widehat{L}}_n\left(i,j\right)$ is greater than an appropriate threshold value.

For the $q_0=q\left(q-1\right)/2$ pairs of ${\widehat{\mathbf{Z}}}_t$, we propose a ratio-based method to those pairs for which $H_0$ will be rejected. We rearrange the $q_0$ pairs obtained ${\widehat{L}}_n\left(i,j\right)$ ‘s in descending order: ${\widehat{L}}_1\ge \cdots \ge {\widehat{L}}_{q_0}$. Define

$$\widehat{d}=arg\underset{1\le j<c_0{q}_0}{max}{\widehat{L}}_j/{\widehat{L}}_{j+1},$$E17

where $c_0\in \left(0,1\right)$ is a prescribed constant. Similar ideas can be found in [16, 19].

3. Numerical results

3.1 Simulation

In this section, we illustrate the finite sample properties of the proposed methodology using simulated data. We only study the performance of the column transformation in (2) since the row transformation in the second stage is essentially the same. As the estimated, $\widehat{\mathbf{A}}$ is an orthogonal matrix for the normalized model in which ${\mathrm{Var}}_r\left({\mathbf{Y}}_t\right)={\mathrm{Var}}_r\left({\mathbf{X}}_t\right)={\mathbf{I}}_q$. We should use ${\widehat{\mathbf{S}}}_{y,0}^{-1/2}\mathbf{A}{\widehat{\mathbf{S}}}_{x,0}^{1/2}$ instead of $\mathbf{A}$ in computing estimation error defined as

$$D\left(\mathcal{M}\left({\mathbf{H}}_1\right),\mathcal{M}\left({\mathbf{H}}_2\right)\right)=\sqrt{1-\frac{1}{min\left(r_1,r_2\right)}\mathrm{tr}\left({\mathbf{P}}_1{\mathbf{P}}_2\right)},$$E18

where ${\mathbf{H}}_i$ is a $p\times r_i$ matrix with rank $\left({\mathbf{H}}_i\right)=r_i$ and ${\mathbf{P}}_i={\mathbf{H}}_i{\left({\mathbf{H}}_{i^{\prime }}{\mathbf{H}}_i\right)}^{-1}{\mathbf{H}}_{i^{\prime }}$ for $i=1,2$. Let ${\mathbf{A}}^{\ast }={\left\{{\mathtt{\Sigma }}_y\left(0\right)\right\}}^{-1/2}\mathbf{A}{\left\{{\mathtt{\Sigma }}_x\left(0\right)\right\}}^{1/2}\equiv \left({\mathbf{A}}_1^{\ast },\dots ,{\mathbf{A}}_{q_1}^{\ast }\right)$, and ${\left\{{\mathtt{\Sigma }}_y\left(0\right)\right\}}^{-1/2}\mathbf{A}=\left({\mathbf{H}}_1,\dots ,{\mathbf{H}}_{q_1}\right)$. Since ${\left\{{\mathtt{\Sigma }}_x\left(0\right)\right\}}^{1/2}$ is a block-diagonal matrix, it holds that $\mathcal{M}\left({\mathbf{H}}_j\right)=\mathcal{M}\left({\mathbf{A}}_j^{\ast }\right)$ for $1\le j\le q_1$. Therefore, we only need to replace $\mathbf{A}$ by ${\widehat{\mathbf{S}}}_{y,0}^{-1/2}\mathbf{A}$.

Since the goal is to specify (via estimation) the $q_1$ linear spaces $\mathcal{M}\left({\mathbf{A}}_j\right),j=1,\dots ,q_1$, simultaneously, we first introduce the concept of a ‘correct’ specification. We call $\widehat{\mathbf{A}}=\left({\widehat{\mathbf{A}}}_1,\dots ,{\widehat{\mathbf{A}}}_{{\widehat{q}}_1}\right)$ a correct specification for $\mathbf{A}$ if (i) ${\widehat{q}}_1=q_1$, and (ii) rank $\left({\widehat{\mathbf{A}}}_j\right)=rank\left({\mathbf{A}}_j\right)$ for $j=1,\dots ,q_1$, after rearranging the order of ${\widehat{\mathbf{A}}}_1,\dots ,{\widehat{\mathbf{A}}}_{{\widehat{q}}_1}$ (we still denote the rearranged sub-matrices as ${\widehat{\mathbf{A}}}_1,\dots ,{\widehat{\mathbf{A}}}_{{\widehat{q}}_1}$ for simplicity in notation). When more than one ${\mathbf{A}}_j$ have the same rank, we pair each of those ${\mathbf{A}}_j$ with the ${\widehat{\mathbf{A}}}_j$ for which

$$D\left(\mathcal{M}\left({\mathbf{A}}_j\right),\mathcal{M}\left({\widehat{\mathbf{A}}}_j\right)\right)=\underset{\mathrm{r}\mathit{ank}\left({\widehat{\mathbf{A}}}_i\right)=\mathrm{r}\mathit{ank}\left({\mathbf{A}}_j\right)}{min}D\left(\mathcal{M}\left({\mathbf{A}}_j\right),\mathcal{M}\left({\widehat{\mathbf{A}}}_i\right)\right).$$E19

Note that a correct specification for $\mathbf{A}$ implies a structurally correct segmentation for ${\mathbf{X}}_t$, which will be abbreviated as ‘correct segmentation’ hereafter. For a correct segmentation, we report the estimation error defined as

$$\overline{D}\left(\widehat{\mathbf{A}},\mathbf{A}\right)=1q_1\sum _{j=1}^{q_1}D\left(\mathcal{M}\left({\mathbf{A}}_j\right),\mathcal{M}\left({\widehat{\mathbf{A}}}_j\right)\right).$$E20

In addition to the correct segmentations, we also report the proportions of the near-complete (NC) segmentations with ${\widehat{q}}_1=q_1-1$ in all the following examples.

Example 1. We consider the model (2) with $p=3$ and $q=6$. The columns of ${\mathbf{X}}_t$ are generated as follows:

$${\mathbf{x}}_i^t={\mathit{\eta }}_{t+i-1}^{\left(1\right)}\left(i=\mathrm{1,2,3}\right),{\mathbf{x}}_i^t={\mathit{\eta }}_{t+i-4}^{\left(2\right)}\left(i=4,5\right),\text{and}{\mathbf{x}}_6^t={\mathit{\eta }}_t^{\left(3\right)},$$E21

where

$${\mathit{\eta }}_t^{\left(j\right)}={\mathtt{\Phi }}^{\left(j\right)}{\mathit{\eta }}_{t-1}^{\left(j\right)}+{\epsilon }_t^{\left(j\right)}-{\mathtt{\Theta }}^{\left(j\right)}{\epsilon }_{t-1}^{\left(j\right)},j=\mathrm{1,2,3}.$$E22

The elements of ${\mathtt{\Phi }}^{\left(j\right)}$ are drawn independently from $U\left(-3,3\right)$_, and then normalized by $0.9\times {\mathtt{\Phi }}^{\left(j\right)}/\parallel {\mathtt{\Phi }}^{\left(j\right)}{\parallel }_2$ so that ${\mathit{\eta }}^{\left(j\right)}$ is stationary, and the elements of ${\mathtt{\Theta }}^{\left(j\right)}$ are drawn independently from $U\left(-1,1\right)$. Meanwhile, the elements of the transformation matrix $\mathbf{A}$ are also drawn independently from $U\left(-3,3\right)$. Thus ${\mathbf{X}}_t$ consists of three independent sub-matrices with, respectively, $3$, $2$, and $1$ columns. In the experiments, we choose $c_0=0.75$ in (17) and $k_0=2$ in (8), and the other $k_0$ ‘s give similar results. The sample sizes $\left(n\right)$ are $100$, $200$, $300$, $400$, $500$, $1000$ and $1500$, and the number of replications is $500$ for each case. The proportions of the correct, incorrect, and near-complete (NC) segmentations are reported in Table 1. From Table 1, we can see that the proposed method improves as the sample size increases. With a dimension of $\mathit{pq}=18$, we can see that the performance is reasonably well even for a small sample size, and the sum of the proportions of the complete and the near complete ones is more than $90\%$ for a small sample size $n=100$, from which we can see that we have achieved sufficient dimension reduction. We next study the estimation errors (20), and the box plots of the errors in the complete segmentations are shown in Figure 1. From Figure 1, we can see that the proposed method improves as the sample size increases.

$n$	100	200	300	400	500	1000	1500
Correct segmentation	0.608	0.668	0.748	0.806	0.864	0.948	0.980
Incorrect segmentation	0.392	0.332	0.252	0.294	0.246	0.052	0.020
NC segmentation with ${\widehat{q}}_1=2$	0.298	0.306	0.250	0.190	0.136	0.052	0.020

Table 1.

The proportions of correct and incomplete segmentations in Example 1. $p=3,q=6$.

Figure 1.

The boxplots of estimation errors $\overline{D}\left(\widehat{\mathit{A}},\mathit{A}\right)$ in Example 1.

As a concrete example, we report the correlogram of ${\mathbf{Y}}_t$ for one replication. We choose $m=10$_, and for $0\le k\le m$ in Figure 2, the correlation between ${\mathbf{y}}_i^t$ and ${\mathbf{y}}_j^t$ are computed by

Figure 2.

Cross correlogram of the ${\mathit{y}}_i^t$ in Example 1, $n=$ 1500.

$$\widehat{\text{Corr}}\left({\mathbf{y}}_i^{t+k},{\mathbf{y}}_j^t\right)={\left|diag{\left\{{\widehat{\mathtt{\Sigma }}}_{y,i,i}\left(0\right)\right\}}^{-1/2}{\widehat{\mathtt{\Sigma }}}_{y,i,j}\left(k\right)diag{\left\{{\widehat{\mathtt{\Sigma }}}_{y,j,j}\left(0\right)\right\}}^{-1/2}\right|}_{\infty }.$$E23

From Figure 2, we can see that each of the columns of ${\mathbf{Y}}_t$ are highly correlated with the others. We then apply our method to ${\mathbf{Y}}_t$_, and Figure 3 depicts the cross correlogram of the transformed matrix ${\mathbf{Z}}_t={\mathbf{Y}}_t{\widehat{\mathtt{\Gamma }}}_y$, and the scales of the $y$-axes are not the same. We can see from Figure 3 that the columns of ${\mathbf{Z}}_t$ can be divided into three groups: $\left\{\mathrm{1,3,6}\right\}$, $\left\{2,4\right\}$, and $\left\{5\right\}$.

Figure 3.

Cross correlogram of the transposed series ${\widehat{\mathit{z}}}_i^t={\mathit{Y}}_t{\widehat{\mathit{v}}}_i$ in Example 1, $n=1500$. The components of ${\widehat{\mathit{Z}}}_t$ can be segmented into 3 groups: {1,3,6}, {2,4}, and {5}.

Example 2. In this example, we slightly increase the dimension as $p=q=6$, and the data-generating processes are the same as those in Example 1. We observe that the performance of the proposed method is not as good as that in Example 1 when the dimension is higher with $\mathit{pq}=36$. Nevertheless, similar conclusions can also be obtained from Table 2 such as the proportion of the complete segmentation increases as the sample size becomes larger, and the sum of the proportions of the complete and the near complete ones is more than $90\%$ even for a small sample size. The box plots of the estimation errors are similar to that in Example 1, and hence we do not report it here. From this example, we can see that when the dimension is higher, the proposed method without thresholding may not work well, and we will illustrate this point in the next example.

$n$	100	200	300	400	500	1000	1500
Correct segmentation	0.504	0.620	0.648	0.670	0.732	0.766	0.780
Incorrect segmentation	0.496	0.380	0.352	0.330	0.268	0.234	0.220
NC segmentation with ${\widehat{q}}_1=2$	0.354	0.320	0.312	0.308	0.264	0.230	0.210

Table 2.

The proportions of correct and incomplete segmentations in Example 2. $p=6,q=6$.

Example 3. We consider model (2) with $p=q=10$_, and the dimension is $\mathit{pq}=100$. The columns of ${\mathbf{X}}_t$ are generated as follows:

$${\mathbf{x}}_i^t={\mathit{\eta }}_{t+i-1}^{\left(1\right)}\left(i=\mathrm{1,2,3,4}\right),{\mathbf{x}}_i^t={\mathit{\eta }}_{t+i-5}^{\left(2\right)}\left(i=\mathrm{5,6,7}\right),{\mathbf{x}}_i^t={\mathit{\eta }}_{t+i-8}^{\left(3\right)}\left(i=8,9\right),\text{and}{\mathbf{x}}_{10}^t={\mathit{\eta }}_t^{\left(4\right)},$$E24

where

$${\mathit{\eta }}_t^{\left(j\right)}={\mathtt{\Phi }}^{\left(j\right)}{\mathit{\eta }}_{t-1}^{\left(j\right)}+{\epsilon }_t^{\left(j\right)}-{\mathtt{\Theta }}^{\left(j\right)}{\epsilon }_{t-1}^{\left(j\right)},j=\mathrm{1,2,3,4},$$E25

and ${\mathtt{\Phi }}^{\left(j\right)}$ and ${\mathtt{\Theta }}^{\left(j\right)}$ are generated in the same way as Example 1. To fulfill a sparsity assumption, let

$${\mathbf{B}}_i=\left(\begin{array}{cc}\cos \left({\theta }_i\pi \right)& sin\left({\theta }_i\pi \right)\\ -sin\left({\theta }_i\pi \right)& cos\left({\theta }_i\pi \right)\end{array}\right),\mathbf{A}=diag\left({\mathbf{B}}_1,{\mathbf{B}}_2,{\mathbf{B}}_3,{\mathbf{B}}_4,{\mathbf{B}}_5\right),$$E26

where ${\theta }_1=\pi /5$, ${\theta }_2=\pi /6$, ${\theta }_3=\pi /7$, ${\theta }_4=\pi /8$, ${\theta }_5=\pi /9$. Thus, $\mathbf{A}$ is a sparse orthogonal transformation matrix. Since the covariance of ${\mathbf{X}}_t$ and ${\mathbf{Y}}_t$ are all block-diagonal by the data generating process, hence ${\widehat{\mathbf{S}}}_{y,0}^{-1/2}{\mathbf{AS}}_{x,0}^{1/2}$ is also block-diagonal and hence satisfies Assumption 3. If we apply our methodology to this model without thresholding the covariance matrices, the results are reported in Table 3. From Table 3, we can see that the proportions of the complete segmentations are pretty low for all the sample sizes, and it does not necessarily improve as the sample size increases. Now, we adopt the thresholding technique as that in [18], and the choice of the threshold is computed by a cross-validation method therein. Table 4 presents the proportions of the correct, incorrect, and near-complete segmentations for model (2) with $\mathit{pq}=100$ dimensions. From Table 4, we can see that the thresholding technique works reasonably well for even small sample sizes. Figure 4 reports the boxplots of the estimation errors, from which we can see that the estimation is very accurate, especially when the sample size is large. Figure 5 displays the cross correlogram of ${\mathbf{Y}}_t$ for one instance of the $500$ replications, and the correlations are calculated in the same way as that in Example 1 with thresholded estimators for the (auto) covariance matrices, and Figure 6 gives the correlogram of the transformed data matrix ${\mathbf{Z}}_t={\mathbf{Y}}_t{\widehat{\mathtt{\Gamma }}}_y$. We can see from Figure 5 that all the columns of ${\mathbf{Y}}_t$ are highly correlated, while those of ${\mathbf{Z}}_t$ can be separated into _four groups: $\left\{\mathrm{1,3,4,8}\right\}$, $\left\{\mathrm{2,5,7}\right\}$, $\left\{6,10\right\}$, and $\left\{9\right\}$. The threshold of the correlation is $0.39$ according to the rule in [18], and those pairs will be treated as uncorrelated ones if their maximal correlation is below this threshold level. Together with Tables 3 and 4 and Figures 4 and 5, we can achieve substantial dimension reduction using our proposed method.

$n$	100	200	300	400	500	1000	1500
Correct segmentation	0.148	0.184	0.186	0.170	0.216	0.174	0.164
Incorrect segmentation	0.852	0.816	0.814	0.830	0.784	0.826	0.836
NC segmentation with ${\widehat{q}}_1=3$	0.216	0.310	0.388	0.470	0.462	0.674	0.744

Table 3.

The proportions of correct and incomplete segmentations in Example 3 without thresholding. $p=$ 10, q $=$ 10.

$n$	100	200	300	400	500	1000	1500
Correct segmentation	0.440	0.572	0.654	0.666	0.840	0.910	0.936
Incorrect segmentation	0.560	0.428	0.346	0.334	0.260	0.080	0.064
NC segmentation with ${\widehat{q}}_1=3$	0.490	0.294	0.282	0.286	0.084	0.062	0.042

Table 4.

The proportions of correct and incomplete segmentations in Example 3 by thresholding. $p=$ 10, q $=$ 10.

Figure 4.

The boxplots of estimation errors $\overline{D}\left(\widehat{\mathit{A}},\mathit{A}\right)$ in Example 3.

Figure 5.

Cross correlogram of the ${\mathit{y}}_i^t$ in Example 3, $n=$ 1500.

Figure 6.

Cross correlogram of the transposed series ${\widehat{\mathit{z}}}_i^t={\mathit{Y}}_t{\widehat{\mathit{v}}}_i$ in Example 3, $n=$ 1500. The components of ${\widehat{\mathit{Z}}}_t$ can be segmented into 4 groups: {1,3,4,8}, {2,5,7}, {6,10}, and {9}.

4. Conclusions

In this chapter, we proposed a new procedure to reduce the dimension of matrix-variate time series by looking for linear transformations to segment the matrix into many small sub-matrices, for which each of them is uncorrelated with the others both contemporaneously and serially, thus they can be analyzed separately, which will greatly reduce the number of parameters to be estimated in terms of modeling. We propose a two-step and more structured procedure to segment the rows and columns separately. When the dimensions are large in relation to the sample size, we assume the transformation matrices are sparse and use threshold estimators for the (auto) covariance matrices. Simulation studies suggest that our procedure works well in segmenting matrix-variate series, and it provides another option for researchers and practitioners to choose from the toolbox in matrix-variate time series modeling.

Acknowledgments

This research was supported in part by the Start-up Fund granted by the “100 Talents Program” of Zhejiang University and the National Natural Science Foundation of China (NSFC, 12201558).

Additional information

An earlier version of this chapter was previously published as a preprint in: Gao Z. Segmenting High-dimensional Matrix-valued Time Series via Sequential Transformations [Internet].arXiv.org. 2020 [cited 2023 Sep 4]. Available from:https://arxiv.org/abs/2002.03382.