Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Simultaneous Decorrelation of Matrix Time Series
\if11
{
} \fi
\if01
{
center[center omitted — 89 chars of source]
} \fi
abstractWe propose a contemporaneous bilinear transformation for a $p\times q$ matrix time series to alleviate the difficulties in modeling and forecasting matrix time series when $p$ and/or $q$ are large.
The resulting transformed matrix assumes a block structure consisting of several small matrices, and those small matrix series are uncorrelated across all times. Hence an overall parsimonious model is achieved by modelling each of those small matrix series separately without the loss of information on the linear dynamics. Such a parsimonious model often has better forecasting performance, even when the underlying true dynamics deviates from the assumed uncorrelated block structure after transformation.
The
uniform convergence rates of the estimated transformation are derived, which vindicate an important virtue of the proposed bilinear transformation, i.e. it is technically equivalent to the decorrelation of a vector time series of dimension max$(p,q)$ instead of $p\times q$. The proposed method is illustrated numerically via both simulated and real data examples.
{\it Keywords:}
Decorrelation transformation;
Eigenanalysis;
Matrix time series;
Forecasting;
Uniform convergence rates.
Introduction
Let $X_t= (X_{t,i,j})$ be a $p\times q$ matrix time series, i.e. there are $p\times q$ recorded values at each time from, for example, $p$ individuals and over $q$ indices or variables. Data recorded in this form are increasingly common in this information age, due to the demand to solve practical problems from, among others, signal processing, medical imaging, social networks, IT communication, genetic linkages, industry production and distribution, economic activities and financial markets. Extensive developments on statistical inference for matrix data under i.i.d. settings can be found in negahban2011estimation,rohde2011,xia2019statistical, and the references therein. Many matrix sequences are recorded over time, exhibiting significant serial dependence which is valuable for modelling and future prediction.
The surge of development in analyzing matrix time series includes bilinear autoregressive models Hoff2015,chen2021autoregressive,xiao2021,
factor models based on Tucker's decomposition for tensors wang2019, chen2022factor, chen2020constrained, han2020, chen2020modeling,han2022rank, and factor models based on the tensor CP decomposition han2021CP,chy2021,han2022tensor.
A common feature in the aforementioned approaches is dimension-reduction, as the so-called `curse-of-dimensionality' is more pronounced in modelling time series than i.i.d. observations, which is exemplified by the limited practical usefulness of vector ARMA (VARMA) models. Note that an unregularized VAR(1) model with dimension $p$ involves at least $p^2$ parameters.
Hence finding an effective way to reduce the number of parameters is of fundamental importance in modelling and forecasting high dimensional time series. In the context of vector time series, most available approaches may be divided into three categories: (i) regularized VAR or VARMA methods incorporating LASSO or alternative penalties basu2015, guo2016high, lin2017, gao2019, zhou2018, ghosh2019, han2020high, han2020sparse, han2023high, (ii) factor models of various forms pena1987identifying, bai2002, forni2005generalized, lam2012, (iii) various versions of independent component analysis back1997first, tiao1989model, huang2014refined, matteson2011dynamic, chang2018. Note that the literature in each of the three categories is large, it is impossible to list all the relevant references here.
In this paper we propose a new parsimonious approach for analyzing matrix time series,
which is in the spirit of independent component analysis. It transforms a $p\times q$ matrix series into the new matrix series of the same size by a contemporaneous bilinear transformation (i.e. the values at different time lags are not mixed together). The new matrix series is divided into several submatrix series and those submatrix series are uncorrelated across all time lags. Hence an overall parsimonious model can be developed
as those submatrix series can be modelled separately without the loss of information on the overall linear dynamics. This is particularly advantageous for forecasting future values. In general cross-serial correlations among different component series are valuable for future prediction.
However the gain from incorporating those cross-serial correlations directly in a moderately high dimensional model is typically not enough to
offset the error resulted from estimating the additional large number of parameters. This is why forecasting a large number of time series together based on a joint time series model can often be worse than that forecasting each component series separately by ignoring cross-serial correlations completely.
The proposed transformation alleviates the problem by channeling all cross-serial correlations into the transformed submatrix series and those subseries are uncorrelated with each other across all time lags. Hence the relevant information can be used more effectively in forecasting within each small models. Note that the transformation is one-to-one, therefore the good forecasting performance of the transformed series can be easily transformed back to the forecasting for the original matrix series.
Our empirical study in Section (ref) also demonstrates that, even when the underlying model does not exactly follow the assumed uncorrelated block structure, the decorrelation transformation often produces superior out-of-sample prediction, as the transformation enhances the within block autocorrelations, and the cross-correlations between the blocks are weaker or too weak to be practically useful.
The basic idea of the proposed transformation is similar to the so-called principal component analysis for time series (TS-PCA) of chang2018. Actually one might be tempted to stack a $p\times q$ matrix series into a $(pq)\times 1$ vector time series, and apply TS-PCA directly. This requires to search for a $(pq)\times (pq)$ transformation matrix. In contrast, the proposed bilinear transformation is facilitated by two matrices of size, respectively, $p\times p$ and $q\times q$. Indeed our asymptotic analysis indicates that the proposed new bilinear transformation is technically equivalent to a TS-PCA transformation with dimension max$(p,q)$ instead of $(p\times q)$; see Remark (ref) in Section (ref) below. Furthermore the bilinear transformation does not mix rows and columns together, as they may represent radically different features in the data. See, e.g., the example in Section (ref) for which the bilinear transformation also outperforms the vectorized TS-PCA approach in a post-sample forecasting.
The rest of the paper is organized as follows. The targeted bilinear decorrelation transformation is characterized in Section (ref). We introduce a new normalization under which the required transformation consists of two orthogonal transformations. The orthogonality allows us to accumulate the information from different time lags without cancellation. The estimation of the bilinear transformation boils down to the
eigenanalysis of two positive-definite matrices of sizes $p\times p$ and $q\times q$ respectively. Hence the computation can be carried out on a laptop/PC with $p$ and $q$ equal to a few thousands. Theoretical properties of the proposed estimation are presented in Section
(ref). The non-asymptotic error bounds of the estimated transformation are derived based on some concentration
inequalities (e.g., rlo2000theorie,merlevede2011bernstein); further leading to the uniform convergence rates of the estimation.
Numerical illustration with both simulated and real data is reported in Section (ref), which shows superior performance in forecasting over the methods without transformation and TS-PCA. Once again it reinforces the fact that the cross-serial correlations are important and useful information for future prediction for a large number of time series, and, however, it is necessary to adopt the proposed decorrelation transformation (or other dimension-reduction techniques) in order to use the information effectively. All technical proofs are relegated to a supplementary.
Methodology
Decorrelation transformations
Again, let $X_t= (X_{t,i,j})$ be a $p\times q$ matrix time series.
We assume that $X_t$ is weakly stationary in the sense that all the first two moments are finite and time-invariant. Our goal is to seek for a bilinear transformation such that the transformed matrix series admits the following segmentation structure:
equation[equation omitted — 326 chars of source]
where $A$ and $B$ are unknown, respectively, $q\times q$ and $p\times p$ invertible constant matrices, and random matrix $U_t$ is unobservable,
$U_{t, i, j}$ is a $p_i\times q_j$ matrix with unknown $p_i$ and $q_j$, $\sum_{i=1}^{n_r}p_i=p$, and $\sum_{j=1}^{n_c}q_j=q$. Furthermore,
Cov$\{ \hbox{\rm vec}({U}_{t+\tau,i_1,j_1}), \hbox{\rm vec}({U}_{t,i_2,j_2})\} =0$ for any $(i_1, j_1)\neq (i_2, j_2)$ and any integer $\tau$,
i.e. all those submatrix series are uncorrelated with each other across all time lags.
remark(i) The decorrelation bilinear transformation ((ref)) is in the same spirit as TS-PCA of chang2018 which transforms linearly a vector time series to a new vector time series of the same dimension but segmented into several subvector series, and those subvectors are uncorrelated across all time lags. Thus as far as the linear dynamics is concerned, one can model each of those subvector time series separately. It leads to appreciable improvement in future forecasting, as the transformed process encapsulates all the cross-serial correlations into the auto-correlations of those uncorrelated subvector processes.
(ii) In the matrix time series setting, one would be tempted to stack all the elements of $ X_t$ into a long vector and to apply TS-PCA of chang2018 directly, though it destroys the original matrix structure.
This requires to search for the $(pq) \times (pq)$ decorrelation transformation matrix $\Phi$ such that $\hbox{\rm vec}( X_t) = \Phi \hbox{\rm vec}(U_t)$. The proposed model ((ref)) is to impose a low-dimensional structure $\Phi=A\otimes B$, where $\otimes$ denotes the matrix Kronecker product, such that the technical difficulty is reduced to that of estimating two transformation matrices of size $p\times p$ and $q\times q$ respectively, with significant faster convergence rate; see Remark (ref) in Section (ref) below.
The empirical results in Section (ref) also demonstrate the benefit of maintaining the matrix structure in out-sample prediction performance.
(iii) The segmentation structure in ((ref)) may be too rigid and the division of the submatrices may not be that regular. In practice a small number of row blocks and/or column blocks may be resulted from applying the estimation method proposed in Section (ref) below. Then applying the same method again to each submatrix series $U_{t, i,j}$ may lead to a finer segmentation with irregularly sized small blocks.
(iv) Similar to TS-PCA, the desired segmentation structure ((ref)) may not exist {in real applications}. Then the estimated bilinear transformation, by the method in Section (ref), leads to an approximate segmentation which encapsulates most significant correlations across different component series into the segmented submatrices while ignoring some small correlations. With enhanced auto-correlations, those submatrix processes become more predictable while those ignored small correlations are typically practically negligible. Consequently the improvement in future prediction
still prevails in spite of the lack of an exact segmentation structure. See chang2018, and also a real data example in Section (ref) below in which the three versions of segmentation (therefore, at least two of them are approximations) uniformly outperform the various methods without the transformation. A simulation study in Section (ref) also shows that when the underlying true model deviates from the desired segmentation structure by a moderate amount, the approximated segmentation model produces superior out-sample prediction performance than those without the transformation.
(v) The use of bilinear form in (ref) is common in matrix and tensor data analysis. For example, it is used in regression with matrix-type covariates basu2012blr,zhou2013tensor,wang2019symmetric, matrix autoregressive model chen2021autoregressive,Hoff2015 and matrix factor model wang2019, chen2020constrained, chen2022factor.
(vi) For dealing with time series with structural changes, it is possible to allow $A$ and $B$ to be time varying, which is technically more demanding and will be explored elsewhere.
Normalization and identification
To simplify the statements, let $\mathbb{E}( X_t) = 0$. This amounts to centre the data by the sample mean first, which will not affect the asymptotic properties of the estimation under the stationarity assumption. Thus $\mathbb{E}({U}_t) = 0$.
The terms ${A}, {B}$ and ${U}_t$ on the RHS of ((ref)) are not uniquely defined, and any ${A}, {B}$ and ${U}_t$ satisfying the required condition will serve the decorrelation purpose. Therefore we can take the advantage from this lack of uniqueness to identify the `ideal' ${A}$ and ${B}$ to improve the estimation effectiveness. More precisely, by applying a new normalization, we identify ${A}$ and ${B}$ to be orthogonal, which facilitates the accumulation of the information from different time lags without cancellation. To this end, we first list some basic facts as follows.
enumerate[(i)]
• $(B, U_t, A^\top)$ in ((ref)) can be replaced by $(B D_1, D_1^{-1} U_t D_2^{-1}, D_2 A^\top)$ for any invertible block diagonal matrices $D_1$ and $D_2$ with the block sizes, respectively, $(p_1, \cdots, p_{n_r})$ and $(q_1, \cdots, q_{n_c})$.
In particular, $D_1$ and $D_2$ can be the permutation matrices which permute the columns and rows within each block. In addition, $D_1$ and $D_2$ can be the block matrices that permute $n_r$ row blocks and $n_c$ column blocks respectively. Note both $(p_1, \cdots, p_{n_r})$ and $(q_1, \cdots, q_{n_c})$ are defined by ((ref)) only up to any permutation.
• The $q$ columns of $B U_t$ can be divided into $n_c$ uncorrelated blocks of sizes $(q_1, \cdots, q_{n_c})$, i.e. the columns of ${B} {U}_t$ resemble the same pattern of the uncorrelated blocks as those of ${U}_t$. Thus ${\Sigma}_1^{(u)} \equiv \mathbb{E}(U_t^\top B^\top B U_t)/p$ is a block diagonal matrix with the block sizes $(q_1, \cdots, q_{n_c})$. In the same vein, ${\Sigma}_2^{(u)} = \mathbb{E}(U_t A^\top A U_t^\top)/q$ is a block diagonal matrix with the block sizes $(p_1, \cdots, p_{n_r})$.
• Put $B=(B_1, \cdots, B_{n_r})$ and $A= (A_1, \cdots, A_{n_c})$, where $B_i$ is $p\times p_i$ and $A_i$ is $q\times q_i$. Then linear spaces ${\cal M}(B_1), \cdots, {\cal M}(B_{n_r})$ and ${\cal M}(A_1), \cdots, {\cal M}(A_{n_c})$ are uniquely defined by ((ref)) up to any permutation, though $B$ and $A$ are not, where ${\cal M}(G)$ denotes the linear space spanned by the columns of matrix $G$.
Note that for any $q\times q$ positive definite matrix $\Sigma_1$ and $p\times p$ positive definite matrix $\Sigma_2$,
\[
{\Sigma}_1^{-1/2} {A} = ({\Sigma}_1^{-1/2} {A}_1, \cdots, {\Sigma}_1^{-1/2} {A}_{n_c}), \qquad {\cal M}({A}_i) =
{\Sigma}_1^{1/2}{\cal M}({\Sigma}_1^{-1/2} {A}_i) \;\;\; i=1, \cdots, n_c,
\]
\[
{\Sigma}_2^{-1/2} {B} = ({\Sigma}_2^{-1/2} {B}_1, \cdots, {\Sigma}_2^{-1/2} {B}_{n_r}), \qquad {\cal M}({B}_i) =
{\Sigma}_2^{1/2}{\cal M}({\Sigma}_2^{-1/2} {B}_i) \;\;\; i=1, \cdots, n_r.
\]
Put
${\Sigma}_1^{(x)} = \mathbb{E}( X_t^\top X_t)/p$ and ${\Sigma}_2^{(x)} = \mathbb{E}( X_t X_t^\top)/q$. Proposition (ref) below indicates that if we replace $ X_t$ in ((ref)) by a normalized version $\{ {\Sigma}_2^{(x)}\}^{-1/2} X_t \{ {\Sigma}_1^{(x)}\}^{-1/2}$, we can treat both ${A}$ and ${B}$ in ((ref)) as orthogonal matrices. This orthogonality plays a key role in combining together the information from different time lags in our estimation (see the proof of Proposition (ref) in the appendix).
propositionLet both ${\Sigma}_1^{(x)}$ and ${\Sigma}_2^{(x)}$ be invertible. Then it holds that
\begin{equation}
X_t = \{ {\Sigma}_2^{(x)}\}^{1/2} {B}_* {U}_t^* {A}_*^\top \{ {\Sigma}_1^{(x)}\}^{1/2},
\end{equation}
where $ {U}_t^*$ admits the same segmentation structure as ${U}_t$ in ((ref)), and ${A}_*$ and $ {B}_*$ are, respectively, $q\times q$ and $p\times p$ orthogonal matrices. More precisely,
\begin{equation}
{U}_t^* = \{ {\Sigma}_2^{(u)}\}^{-1/2} {U}_t \{ {\Sigma}_1^{(u)}\}^{-1/2},
{A}_* = \{ {\Sigma}_1^{(x)}\}^{-1/2} {A} \{ {\Sigma}_1^{(u)}\}^{1/2},
{B}_* =\{ {\Sigma}_2^{(x)}\}^{-1/2} {B} \{ {\Sigma}_2^{(u)}\}^{1/2}.
\end{equation}
The proof of Proposition (ref) is almost trivial. First note that both ${\Sigma}_1^{(u)}, \, {\Sigma}_2^{(u)}$ are invertible. Then ((ref)) follows from ((ref)) and ((ref)) directly. The orthogonality of, e.g., ${B}_*$ follows from equality $\mathbb{E}[\{ {\Sigma}_2^{(x)}\}^{-1/2} X_t X_t^\top \{ {\Sigma}_2^{(x)}\}^{-1/2}]/q={I}_p$, ((ref)) and ((ref)). Also note that $\{ {\Sigma}_2^{(u)}\}^{-1/2}$ and $\{ {\Sigma}_1^{(u)}\}^{-1/2}$ are of the same block diagonal structure as, respectively, ${\Sigma}_2^{(u)}$ and ${\Sigma}_1^{(u)}$. This implies that ${U}_t^*$ admits the same segmentation structure as ${U}_t$.
Write $A_* = (A_{*1}, \cdots,A_{*n_c}), B_* = (B_{*1},\cdots,B_{*n_r})$, where $A_{*j}$ has $q_j$ columns and $B_{*i}$ has $p_i$ columns. Then $U_{t,i,j}^*=A_{*j} \{ \Sigma_2^{(x)}\}^{-1/2} X_t \{ {\Sigma}_1^{(x)}\}^{-1/2} B_{*i}$.
Note that $({B}_*, {U}_t^*, {A}_*^\top)$ in ((ref)) are (still) not uniquely defined, similar to the property (i) above.
In fact, only the linear spaces ${\cal M}(B_{*1}), \cdots, {\cal M}(B_{*n_r})$ and ${\cal M}(A_{*1}), \cdots, {\cal M}(A_{*n_c})$ are uniquely defined.
Proposition (ref) below shows that we can take the orthonormal eigenvectors of two properly defined positive definite matrices as the columns of ${A}_*$ and ${B}_*$. With ${A}_*$ and ${B}_*$ specified, the segmented ${U}_t^*$ can be solved from ((ref)) directly. Let
eqnarray*[eqnarray* omitted — 282 chars of source]
where ${E}_{i,j} $ is the unit matrix with 1 at position $(i,j)$ and 0 elsewhere, and ${E}_{i,j} $ is $p\times p$ in the first equation, and $q\times q$ in the second equation. For a prespecified integer $\tau_0 \ge 1$, let
equation[equation omitted — 304 chars of source]
propositionLet both ${\Sigma}_1^{(x)}$ and ${\Sigma}_2^{(x)}$ be invertible, and all the eigenvalues of ${W}^{(i)}$ be distinct, $i=1, 2$. Also let $\tau_0\ge1$. Then the $q$ orthonormal eigenvectors of ${W}^{(1)}$ can be taken as the columns of ${A}_*$, and the $p$ orthonormal eigenvectors of ${W}^{(2)}$ can be taken as the columns of ${B}_*$.
The proposition above does not give any indication on how to arrange the columns of ${A}_*$ and ${B}_*$, which should be ordered according to the latent uncorrelated block structure ((ref)). We address this issue in Section (ref) below.
remarkThe representation (ref) suffers from the indeterminacy due to the fact that the column blocks and row blocks can be permuted, and the columns and rows within each block can be rotated, see property (i) above. However this indeterminacy does not have material impact on identifying the desired structure ((ref)), as any one of such representations will serve the purpose. The developed asymptotic theory guarantees that the proposed estimator converges to a representation which fulfills the conditions imposed on ((ref)).
remarkIn case that $W^{(1)}, W^{(2)}$ have tied eigenvalues and the corresponding eigenvectors across different blocks, Proposition (ref) no longer holds. Then the blocks sharing the tied eigenvalues cannot be separated. To avoid the tied eigenvalues, we may use different values of $\tau_0$, or different form of $W^{(i)}$. For example, we may replace $W^{(1)}$ by
\begin{eqnarray*}
W^{(1,f)} = \sum_{\tau=-\tau_0}^{\tau_0} \sum_{i=1}^p\sum_{j=1}^p
\frac{f\left( V^{(1)}_{\tau,i,j} V^{(1)\top}_{\tau,i,j} \right)}{p^2}.
\end{eqnarray*}
where $f(V)$ is a function of a symmetric matrix $V$ in which $f(V)=\Gamma D^{(f)}\Gamma^\top$, $V=\Gamma D\Gamma^\top$ is the eigen-decomposition for $V$, and $D^{(f)}=\hbox{\rm diag}(f(d_1),\ldots,f(d_q))$ is the diagonal-element-wise transformation of the diagonal matrix $D$ of the eigenvalues.
In fact the condition that all eigenvalues of $W^{(i)}$ are different can be relaxed. For example, we only require any two eigenvalues of $W^{(i)}$ corresponding two different blocks to be different while the eigenvalues corresponding to the same black can be the same. See also a similar condition in the eigen-gap $\Delta$ in (ref) in Section (ref) below.
remarkIn practice, we need to specify $\tau_0$ in ((ref)). In principle any $\tau_0\geq 1$ can be used for the purpose of segmentation. A larger $\tau_0$ may capture more lag-dependence, but may also risk adding more `noise' when the dependence decays fast. Since autocorrelation is typically at its strongest at small lags, a relatively small $\tau_0$, such as $\tau_0\le 5$, is often sufficient lam2011,chang2018,wang2019,chen2022factor.
Estimation
The estimation for ${A}_*$ and ${B}_*$, defined in (ref), is based on the eigenanalysis of the sample versions of matrices defined in ((ref)). To this end, let
equation[equation omitted — 168 chars of source]
equation[equation omitted — 387 chars of source]
equation[equation omitted — 361 chars of source]
Performing the eigenanalysis for $\widehat{W}^{(1)}$, and arranging the order of the resulting $q$ orthonormal eigenvectors $\widehat \gamma_1, \cdots, \widehat \gamma_q$ by the algorithm below, we take the
re-ordered eigenvectors as the columns of $\widehat {A}_*$. The estimator $\widehat {B}_*$ is obtained in the same manner via the eigenanalysis for $\widehat{W}^{(2)}$. Then by ((ref)), we obtain the transformed matrix series
equation[equation omitted — 165 chars of source]
Note that the estimation for ${A}_*$ and that for ${B}_*$ are carried out separately. They do not interfere with each other.
Now we present an algorithm to determine the order of the columns for $\widehat {A}_*$. By ((ref)), the columns of $ Z_t \equiv X_t \{ {\Sigma}_1^{(x)}\}^{-1/2} {A}_*$ are divided into the $n_c$ uncorrelated blocks. Define
equation[equation omitted — 168 chars of source]
We divide the columns of $\widehat {Z}_t $ into uncorrelated blocks according to the pairwise maximum cross correlations between the columns. Specifically, the maximum cross correlation between the $k$-th and $\ell$-th columns is defined as
align[align omitted — 501 chars of source]
The second equality above follows from (ref), (ref) and (ref).
In the above expression, $\tau_1 \ge 1$ is a user-defined tuning parameter (See Remark (ref) below).
To determine all significantly correlated pairs of variable, rearrange $\widehat\rho_{k,\ell}$, $1\le k < \ell \le q$, in the descending order: $\widehat \rho_{(1)} \ge \cdots \ge \widehat \rho_{(q_0)}$ and define
align[align omitted — 132 chars of source]
where $\delta_T>0$ is a small constant.
We take the $\widehat r$ pair of columns corresponding to $\widehat \rho_{(1)}, \cdots, \widehat \rho_{(\widehat r)}$ as correlated pairs, and treat the rest as uncorrelated pairs. The intuition is that $\rho_{(r)}/\rho_{(r+1)}$ is $\infty$ if $\rho_{(r)}> 0$ but $\rho_{(r+1)}=0$. The use of $\delta_T$ is to smooth out the large variation of $\widehat \rho_{(j)}/\widehat \rho_{(j+1)}$ when both $\rho_{(j)}$ and $\rho_{(j+1)}$ are small. Similar ideas have also been used in determining the number of factors in lam2012, ahn2013 and han2022rank.
With the $\widehat r$ identified significantly correlated pairs, a connection graph is built, with the column indexes as the vertices and the edges between the identified correlated column pairs. The number of unconnected sub-graphs is the estimated number of blocks $\widehat{n}_c$, and each of maximum connected sub-graphs forms the estimated groups $\widehat{G}_i$, $i=1,\ldots \hat{n}_c$. The corresponding $\widehat{n}_c$ groups of $\widehat \gamma_1, \cdots, \widehat \gamma_q$ are taken as the columns of
$\widehat A_{*i}$, $i=1,\ldots, \hat{n}_c$.
Algorithmically,
start with $q$ groups with one column in each group; then iteratively check all pairs of groups and merge two groups together if there exists at least one pair of columns (one in each group) are significantly correlated; and stop when no groups can be merged.
The finite sample performance of the above algorithm can be improved by prewhitening each column time series of $\widehat {Z}_t$. This makes $\widehat \rho_{k,\ell}$, for different $(k, \ell)$, more comparable. See Remark 2(iii) of chang2018. In practice the prewhitening can be carried out by fitting each column time series a VAR model with the order between 0 and 5 determined by AIC. The resulting residual series is taken as a prewhitened series.
remark(i) The ratio estimator in ((ref)) picks the $\widehat r$ most correlated pairs of columns, and ignores the other small correlations in constructing the segmentation structure. Theorem (ref) in Section (ref) below shows that the partition of $\widehat A_*$ into $\{\widehat A_{*1}, \cdots, \widehat A_{*\widehat n_c}\}$, determined by the above algorithm, provides a consistent estimation for the column segmentation of ${A}_*$.
(ii) There are two tuning parameters used in the procedure, the maximum lag $\tau_0$ used in constructing $W^{(1)}$ and $W^{(2)}$ in (ref) and the maximum lag $\tau_1$ used in measuring the dependency between two columns (rows) in (ref). Lag $\tau_1$ is usually a sufficiently large integer, for example, between 10 and 20, in the spirit of the rule of thumb of box1970. Different values of $\tau_0$ and $\tau_1$
may, or may not, lead to different segmentation.
However the impact on, for example, future prediction is minimum, as the most information on linear dynamics is encapsulated in the most correlated pairs, as indicated by numerical examples in the Appendix.
The ordering for columns of $\widehat {B}_*$ is arranged, in the same manner as above, by examining the pairwise correlations among the columns of
\[
X_t^\top \{ \widehat {\Sigma}_2^{(x)}\}^{-1/2} (\widehat {\gamma}_1^{(2)},
\cdots, \widehat {\gamma}_p^{(2)}),
\]
where $\widehat {\gamma}_1^{(2)}, \cdots, \widehat {\gamma}_p^{(2)}$ are now the $p$ orthonormal eigenvectors of $\widehat {W}^{(2)}$.
Theoretical properties
To gain more appreciation of the methodology, we will show the consistency of the proposed detection method of uncorrelated components. We mainly focus on the estimation error and the ordering of $\widehat A_*$, as those for $\widehat B_*$ are similar.
Denote by $ X_{t,i,\cdot}$ the $i$-th row of $ X_t$, and $X_{t,i,j}$ the $(i,j)$-th element of $ X_t$.
For matrix ${V}=(V_{ij})$, let $\|{V}\|_{\rm op}$ denote the spectrum norm,
and $\|{V}\|_{\max} = \max_{i,j} |V_{ij}|$. We introduce some regularity conditions first.
assumptionAssume exponentially decay coefficients of the strong $\alpha$-mixing condition,
$\alpha(k) \le \exp\big( - c_0 k^{r_1} \big)$
for some constant $c_0>0$ and $0<r_1\le 1$, where
\begin{equation}
\alpha(k) = \sup_t\Big\{\Big|{\mathbb{P}}({\cal E}_1\cap {\cal E}_2) - {\mathbb{P}}({\cal E}_1){\mathbb{P}}({\cal E}_2)\Big|:
{\cal E}_1\in \sigma( X_s, s\le t), {\cal E}_2\in \sigma( X_s, s\ge t+k)\Big\}.
\end{equation}
Here for any random variable/vector/matrix $ X$, $\sigma( X)$ is understood to be the $\sigma$-field generated by $ X$. Assumption (ref) allows a very general class of time series models, including causal ARMA processes with continuously distributed innovations; see bradley2005, and also Section 2.6 of fan2008nonlinear. The restriction $r_1\le 1$ is introduced only for simple presentation.
assumptionAssume $\mathbb{E} X_t=0$. There exist certain finite constants $c_*, c^*$, such that
\begin{equation}
\sup_{i\le p, \|{u}\|_2=1} \hbox{\rm var}\Big( X_{t,i,\cdot} {u}\Big) \le c^*,\quad
\inf_{\|{u}\|_2=1}\mathbb{E} \| X_t{u} \|_2^2/p\ \ge\ c_*.
\end{equation}
Note that $\sup_{ \|{u}\|_2=1}\hbox{\rm var}\Big( X_{t,i,\cdot}^\top {u}\Big)= \|\mathbb{E} X_{t,i\cdot}^\top X_{t,i\cdot} \|_{\rm op}$ and $\mathbb{E} \| X_t{u} \|_2^2/p=u^\top {\Sigma}_1^{(x)}u$. Assumption (ref) on the eigenvalues is a common assumption in the high dimensional setting, for instance, bickel2008covariance,bickel2008regularized, xia2017.
assumptionFor any $x>0$,
$\max_{1\le i\le p,1\le j\le q} {\mathbb{P}} \left(|X_{t,i,j}|\ge x \right) \le c_1 \exp\left(-c_2 x^{r_2} \right) ,$
for some constant $c_1,c_2>0$ and $0<r_2\le 2$.
Assumption (ref) requires that the tail probability of each individual series of $ X_t$ decay exponentially fast. In particular, when $r_2=2$, each $X_{t,i,j}$ is sub-Gaussian.
Write the eigenvalue-eigenvector decomposition of ${W}^{(1)}$ as
eqnarray*[eqnarray* omitted — 76 chars of source]
where ${D}=\hbox{\rm diag}(\lambda_1,\ldots,\lambda_q)$
with $\lambda_1\ge \cdots \ge \lambda_q$, and ${\Gamma}^{(1)} = ({\gamma}_1,\ldots,{\gamma}_q)$. By the structure assumption in ((ref)),
the index set $\{1, \cdots, q\}$ is partitioned
into $n_c$ subsets
$ G_1,...,G_{n_c}$ such that {the columns of different sub-matrices}
$ X_t\Big({\Sigma}^{(1)}_1\Big)^{-1/2}{\Gamma}^{(1)}_{G_i},\ i = 1,\ldots,n_c,$
are uncorrelated with each other across all time lags.
Define eigen-gap
align[align omitted — 130 chars of source]
The following theorem provides the non-asymptotic bounds for the estimators of ${V}^{(1)}_{\tau,i,j}$, ${W}^{(1)}$ {and $ A_{*j}$}, $1\le i,j\le q$ under both exponential decay $\alpha$ mixing condition and exponential tail condition of $ X_t$.
theoremSuppose Assumptions (ref), (ref), (ref) hold with constants $c_*,c^*, r_1,r_2$, and $\tau_0$ is a finite constant. Let $1/\beta_1=1/r_1+2/r_2$ and $1/\beta_2=1/r_1+1/r_2$. Then,
\begin{align}
&\big\|\widehat{ V}^{(1)}_{\tau,i,j} - {V}^{(1)}_{\tau,i,j}\big\|_{\rm op} \le C_1 \eta_{T,p,q}, \;\; \forall 1\le \tau\le \tau_0, 1\le i,j\le p, \\
& \big\|\widehat{ W}^{(1)} - {W}^{(1)} \big\|_{\rm op} \le C_1 \eta_{T,p,q},
\end{align}
in an event $\Omega_T$ with probability at least $1 -\epsilon_T$, where
\begin{eqnarray}
\eta_{T,p,q} = q\left(\sqrt{\frac{\log(pq/\epsilon_T)}{T}} + \frac{[\log(Tpq/\epsilon_T)]^{1/\beta_1}} {T} + \frac{[\log(Tpq/\epsilon_T)]^{2/\beta_2}} {T^2}\right),
\end{eqnarray}
$C_1$ is a constant depending on $c_*,c^*,r_1,r_2$ only. Moreover, there exists $\widetilde A_* \equiv (\widetilde A_{*1}, \cdots,
\widetilde A_{*n_c}) $ of which the columns are a permutation
of $(\widehat \gamma_1, \cdots, \widehat\gamma_q)$ such that
\begin{eqnarray}
\big\|\widetilde A_{*j} \widetilde A_{*j}^\top - A_{*j} A_{*j}^\top \big\|_{\rm op} \le C_2 \eta_{T,p,q} \quad {\rm for}
\quad 1\le j\le n_c
\end{eqnarray}
in the same event $\Omega_T$ provided $C_1 \eta_{T,p,q} \le \Delta/2$, where $C_2$ is a constant depending on $c_*,c^*,r_1,r_2$ only.
remark(i) Let $P_G$ be the projection operator onto the column space of $G$, which can be written as $P_G=G (G^\top G)^{-1} G^\top$.
Under the conditions of Theorem (ref), we can obtain {$\| P_{\widetilde A_{*j}\otimes \widetilde B_{*i}} - P_{A_{*j}\otimes B_{*i}}\|_{\rm op}=O_{{\mathbb{P}}} (\max\{p,q\} [(\log(pq)/T)^{1/2} + \log(Tpq)^{1/\beta_1}/T])$ for each sub-group $1\le i\le n_r, 1\le j\le n_c$}, where the columns of $\widetilde B_*$ are a permutation of the
orthonormal eigenvectors of ${\widehat W}_{(2)}$. If we stack all the elements of $X_t$ into a long vector such that $\hbox{\rm vec}( X_t) = {\Phi} \hbox{\rm vec}({U}_t)$ and apply TS-PCA of chang2018 directly, the estimation of the decorrelation transformation matrix would satisfy {$\| P_{\widetilde\Phi_k} - P_{\Phi_k}\|_{\rm op}=O_{{\mathbb{P}}} (pq[(\log(pq)/T)^{1/2} + \log(Tpq)^{1/\beta_1}/T])$ for $1\le k\le n_r n_c$ and $\Phi=(\Phi_1,...,\Phi_{n_r n_c})$}. Obviously, $\| P_{\widetilde A_{*j}\otimes \widetilde B_{*i}} - P_{A_{*j}\otimes B_{*i}}\|_{\rm op}$ has much sharper rate which is equivalent to that for the TS-PCA estimation with dimension max$(p,q)$.
(ii) The consistency of $\widehat W^{(1)}$ requires $\max(p,q)=o(\sqrt{T})$. When $A$ and $B$, or equivalently $W^{(1)}$ and $W^{(2)}$, have certain sparsity structures, this condition can be further relaxed. In TS-PCA, threshold estimator of bickel2008covariance is employed to construct a sparse $\widehat W^{(1)}$, and then the convergence rates are improved to allow the dimension to be much larger than $T$ chang2018.
The similar extension can be established in our setting.
When the decays of the $\alpha$-mixing coefficients and the tail probabilities of $ X_t$ are slower than Assumptions (ref) and (ref), we impose the following Assumptions (ref) and (ref) {instead}. These conditions ensure the Fuk-Nagaev-type inequalities for $\alpha$-mixing processes; see also rlo2000theorie, liu2013probability, wu2016performance and zhang2017gaussian.
assumptionLet the $\alpha$-mixing coefficients satisfy the condition
$\alpha(k) \le c_0 k^{-r_1},$
where $\alpha(k)$ is defined in (ref), $c_0>0$ and $r_1>1$.
assumptionFor any $x>0$,
$\max_{1\le i\le p,1\le j\le q} {\mathbb{P}} \left(|X_{t,i,j}|\ge x \right) \le c_1 x^{-2r_2} ,$
for some constant $c_1>0$ and $r_2>1$.
Theorem (ref) below presents the uniform convergence rate for $\widehat{ V}^{(1)}_{\tau,i,j}$ and $\widehat{ W}^{(1)}$, $1\le i,j\le q$. It is intuitively clear that the rate is much slower than that in Theorem (ref).
theoremSuppose Assumptions (ref), (ref), (ref) hold with constants $c_*,c^*, r_1,r_2$, and $\tau_0$ is a finite constant. Let $\beta_3=r_2(r_1+1)/(r_1+r_2)>1$. Then, (ref) and (ref) hold in an event $\Omega_T$ with probability at least $1-\epsilon_T$, where
\begin{eqnarray}
\eta_{T,p,q} = q\left( \frac{(Tpq/\epsilon_T)^{1/\beta_3}}{T} + \sqrt{\frac{\log(pq/\epsilon_T)}{T}} \right) ,
\end{eqnarray}
and $C_1$ depends on $c_*,c^*,r_1,r_2$ only. Similarly, if $C_1 \eta_{T,p,q} \le \Delta/2$, then there exists $\widetilde A_* \equiv (\widetilde A_{*1}, \cdots,
\widetilde A_{*n_c}) $ of which the columns are a permutation of $(\widehat \gamma_1, \cdots, \widehat\gamma_q)$ such that (ref) holds in the same event $\Omega_T$ for some positive constant $C_2$ depending on $c_*,c^*,r_1,r_2$ only.
Theorem (ref) below implies that the partition
of $\widehat A_*$ into $\{ \widehat A_{1*}, \cdots, \widehat A_{*\widehat n_c}\}$ in Section (ref)
provides a consistent estimation
for the column segmentation of $A_*
=(A_{*1}, \cdots, A_{*n_c})$. To this end, let
eqnarray[eqnarray omitted — 285 chars of source]
See also ((ref)).
By ((ref)), the columns of $ Z_t \equiv X_t \{ {\Sigma}_1^{(x)}\}^{-1/2} {A}_*$ are divided into the $n_c$ uncorrelated blocks.
We assume that those $n_c$ blocks can also be obtained by
the algorithm in Section (ref)
with $\widehat \rho_{k, \ell}$ replaced by $\rho_{k,\ell}\cdot I(\rho_{k,\ell}\ge \rho)$ for some constant $\rho>0$.
Denote by $G_1, \cdots, G_{n_c}$, respectively, the indices of the components of $Z_t$ in each of those
$n_c$ blocks. Denote by $\widehat G_1, \cdots, \widehat G_{\widehat n_c}$, respectively, the indices of the components
of $\widehat Z_t$ in each of the $\widehat n_c$ uncorrelated blocks identified by the algorithm in Section (ref). Re-arrange of the order of
$\widehat G_1, \cdots, \widehat G_{\widehat n_c}$ if necessary. Then the theorem below implies that
${\mathbb{P}}(\hat n_c =n_c) \to 1$ and ${\mathbb{P}}(\widehat G_i = G_i | \hat n_c = n_c) \to 1$ for $1\le i \le n_c$.
theoremSuppose conditions of Theorem (ref) or (ref) hold and $\tau_1$ is a finite constant.
Suppose
\begin{align}
\kappa_1 \eta_{T,p,q}<\Delta, \qquad \kappa_2 \eta_{T,p,q} < \rho^* ,
\end{align}
where $\kappa_1,\kappa_2$ are certain positive constants depending on $c_*,c^*,r_1,r_2$ only, $\eta_{T,p,q}$ is defined in (ref) or (ref) and depends on $\epsilon_T$. Then in an event $\Omega_T$ with probability at least $1 -\epsilon_T$, we have
\[
\widehat n_c=n_c \mbox{\ \ and \ \ }
\widehat G_i=G_i,\quad 1\le i\le n_c.
\]
remarkThe first inequality in (ref) requires that the minimum eigen-gap $\Delta$ between different uncorrelated groups is sufficiently larger than the estimation error $\eta_{T,p,q}$, such that all the groups are identifiable. The second inequality in (ref) ensures that there are no cross-group edges among $G_1, \cdots, G_{n_c}$.
The constants $\kappa_1$ and $\kappa_2$ are specified in the proof of the theorem.
Numerical properties
Simulation
We illustrate the proposed decorrelation method with simulated examples. We always set $\tau_0=5$ in ((ref)), and $\tau_1 =15$ in ((ref)).
The prewhitening, as stated at the end of Section (ref), is always applied in determining the groups of the columns of $\widehat {A}_*$ and $\widehat {B}_*$. As the true transformation matrices ${A}_*$ and ${B}_*$ defined in ((ref)) cannot be computed easily even for simulated examples, we use the following proxies instead
\[
{A}_*= \{ \widehat {\Sigma}_1^{(x)}\}^{-1/2} {A} \{ \widehat {\Sigma}_1^{(u)} \}^{1/2}, \qquad
{B}_*= \{ \widehat {\Sigma}_2^{(x)}\}^{-1/2} {B} \{ \widehat {\Sigma}_2^{(u)} \}^{1/2},
\]
where $\widehat {\Sigma}_1^{(x)}, \; \widehat {\Sigma}_2^{(x)}$ are given in ((ref)), and
\[
\widehat {\Sigma}_1^{(u)}= {1 \over T p} \sum_{t=1}^T {U}_t^\top {B}^\top {B} {U}_t, \qquad
\widehat {\Sigma}_2^{(u)}= {1 \over T q} \sum_{t=1}^T {U}_t {A} {A}^\top {U}_t^\top.
\]
To abuse the notation, we still use ${A}_*, \; {B}_*$ to denote their proxies in this section since the true ${A}_*, \; {B}_*$ never occur in our computation. For each setting, we replicate the simulation 1000 times.
exampleWe consider model ((ref)) in which all the component series of ${U}_t$ are independent and ARMA(1, 2) of the form:
\begin{equation}
z_t = b z_{t-1} + \epsilon_t + a_1 \epsilon_{t-1} + a_2 \epsilon_{t-2},
\end{equation}
where $b$ is drawn independently from the uniform distribution on $(-.98,\,-0.5)\cup (0.5,\, 0.98)$, $a_1, a_2$ are drawn independently from the uniform distribution on $(-.98,\,-0.3)\cup (0.3,\, 0.98)$, and $\epsilon_t$ are independent and $N(0,1)$. The elements of ${A}$ and ${B}$ are drawn independently from $U(-1, \, 1)$. For this example, we assume that we know
$n_r=p$ and $n_c=q$ in ((ref)). The focus is to investigate the impact of $p$ and $q$ on the estimation for ${B}$ and ${A}$.
Performing the eigenanalysis for $\widehat {W}^{(1)}$ and $\widehat {W}^{(2)}$ defined in ((ref)), we take the resulting two sets of orthonormal eigenvectors as the columns of $\widehat {A}_*$ and $\widehat {B}_*$ respectively. To measure the estimation error, we define
\begin{align*}
D(\widehat {A}_*, \, {A}_*) = \frac{1}{2 q (\sqrt{q} -1)} \sum_{j=1}^q \big( \frac{1}{\max_i |d_{i,j}|}
+ \frac{1}{\max_i |d_{j,i}|}
-2 \big),
\end{align*}
where $d_{i,j}$ is the $(i,j)$-th element of matrix $\widehat{A}_{*}^{\top} {A}_{*}$. Note that $D(\widehat {A}_*, \, {A}_{*})$ is always between 0 and 1,
and $D(\widehat{A}_*, \, {A}_{*})=0$ if $\widehat{A}_*$ is a column permutation of ${A}_{*}$.
\begin{table}[ht]\caption{Example 1 -- Means and standard errors (SE) of $D(\widehat{A}_*, \, {A}_{*})$ and $D(\widehat{B}_*, \, {B}_{*})$ in a simulation with 1000 replications.}
\scriptsize
\begin{center}
\begin{tabular}{r|c|cc|cc||c|cc|cc|}
& & \multicolumn{2}{c|}{$D(\widehat{A}_*, \,
{A}_{*})$} &
\multicolumn{2}{c||}{$D(\widehat{B}_*, \, {B}_{*})$}
& & \multicolumn{2}{c|}{$D(\widehat{A}_*, \,
{A}_{*})$} &
\multicolumn{2}{c|}{$D(\widehat{B}_*, \, {B}_{*})$}\\
$T$ & ($q,\; p$) & Mean& SE & Mean & SE & ($q,\; p$) & Mean& SE & Mean & SE\\
\hline
100 & (4, 4) & 0.135 & 0.103 & 0.116 & 0.094 & (4, 8) & 0.124 & 0.096& 0.169 & 0.056 \\
500 & & 0.081 & 0.080 & 0.052 & 0.065 & & 0.074 & 0.078& 0.089& 0.047 \\
1000& & 0.048 & 0.062 & 0.050 & 0.063 & & 0.057 & 0.066& 0.056& 0.037 \\
5000& & 0.021 & 0.042 & 0.020 & 0.042 & & 0.019 & 0.039& 0.032& 0.027 \\
\hline
100 & (8, 8)& 0.172 & 0.057 & 0.172 & 0.055 & (16, 16) & 0.204 & 0.031 & 0.206 & 0.032 \\
500 & & 0.098 & 0.048 & 0.086 & 0.044 & & 0.132 & 0.032 & 0.131 & 0.031 \\
1000& & 0.074 & 0.045 & 0.070 & 0.041 & & 0.097 & 0.027 & 0.093 & 0.026 \\
5000& & 0.034 & 0.029 & 0.035 & 0.030 & & 0.041 & 0.017 & 0.042 & 0.018 \\
\hline
100& (32, 32) & 0.224 & 0.018 & 0.225 & 0.018 & (100, 16)& 0.251 & 0.006 & 0.187 & 0.032 \\
500 & & 0.162 & 0.019 & 0.161 & 0.020 & & 0.212 & 0.008 & 0.127 & 0.031 \\
1000& & 0.127 & 0.018 & 0.125 & 0.018 & & 0.185 & 0.009 & 0.099 & 0.028 \\
5000& & 0.058 & 0.013 & 0.058 & 0.013 & & 0.106 & 0.008 & 0.044 & 0.019 \\
\end{tabular}
\end{center}
\end{table}
We set $T=100, 500, 1000, 5000$, and $p, q = 4, 8, 16, 32, 100$. The means and standard errors of $D(\widehat{A}_*, \, {A}_{*})$ and $D(\widehat{B}_*, \, {B}_{*})$ over the 1000 replications are reported in Table (ref). As expected, the estimation errors decrease as $T$ increases. Furthermore the error in estimating ${A}_*$ increases as $q$ increases, and that in estimating ${B}_*$ increases as $p$ increases. Also noticeable is the fact that the quality of the estimation for ${A}_*$ depends on $q$ only and that for ${B}_*$ depends on $p$ only, as $D(\widehat{A}_*, \, {A}_{*})$ is about the same for $(q,p)=(4, 4)$ and $(q, p)=(4,8)$, and $D(\widehat{B}_*, \, {B}_{*})$ is about the same for $(q, p)=(4,8)$ and $(q, p)=(8,8)$. When $(q,p)=(32, 32)$, each ${A}$ and ${B}$ contains $1024$ unknown parameters. The quality of their estimation with $T=100$ is not ideal, but not substantially worse than that for $(q,p)=(16, 16)$, and the estimation is very accurate with $T=5000$. For $(q,p)=(100, 16)$, matrix time series $ X_t$ contains 1600 component series. The estimation for ${A}_*$ is not accurate enough as $q=100$, while the estimation for ${B}_*$ remains as good as that for $(q,p)=(16,16)$, and is clearly better than that for $(q,p)=(32,32)$.
exampleTo examine the performance of the proposed method for identifying uncorrelated blocks, we consider now a only column-segmentation model
\[
X_t = {U}_t {A}^\top,
\]
where $ X_t, {U}_t$ are $p\times q$ matrix time series. Note that adding the row transformation matrix ${B}$ in the model entails the application of the same method twice without adding any new insights on the performance of the method, as the estimation for ${A}_*$ and that for ${B}_*$ are performed separately. See Section (ref).
The transformation matrix ${A}$ in the above model is generated in the same way as in Example 1 above. All the element time series
of ${U}_t$, except those in columns 2, 3 and 5, are simulated independently from ARMA(1,2) model ((ref)). Denoted by ${U}_{t, j}$ the $j$-th column of ${U}_t$,
we let
\[
{U}_{t,2}={U}_{t+1,1}, \quad {U}_{t,3}={U}_{t+2,1}, \quad
{U}_{t,5}={U}_{t+1,4}.
\]
Hence the first three columns of ${U}_t$ forms one block of size 3, columns 4 and 5 forms a block of size 2, and the rest of the columns are
uncorrelated with each other, and are uncorrelated to the first two blocks.
We set $T=100, 500, 1000, 5000$. For each sample size, we set $(p,q)=(3,6)$, $(6,6)$, $(6, 10)$ and $(10, 10)$. Hence the number of groups is $n_c=3,3,7$ and $7$, respectively. For each setting, we perform the eigenanalysis for $\widehat {W}^{(1)}$ defined as ((ref)). We then apply the algorithm in Section (ref) to arrange the orders of the $q$ resulting orthonormal eigenvectors to form estimator $\widehat {A}_*$, for which the number of the correlated pairs of columns is determined by ((ref)).
Table (ref) reports the relative frequencies of the correct specification of all the $n_c$ uncorrelated blocks. We also report the relative frequencies of two types of wrong segmentation: (a) Merging with $\widehat n_c = n_c -1$, i.e. $n_c-2$ blocks are correctly specified, and the remaining two blocks are put into one; (b) Splitting with $\widehat n_c = n_c+1$, i.e. $n_c -1$ blocks are correctly specified, and the remaining block
incorrectly splits into two. Table (ref) indicates that the relative frequency of the correct specification increases as the sample size $T$ increases, and it decreases as $q$ increases while the performance is much less sensitive to the increase of $p$.
It is noticeable that the probability of the event $\{ \widehat n_c = n_c +1 \}$ is very small especially when $T$ is large. With $q=6$, the probability for the correct specification is not smaller than 64% with $T=500$, and is greater than 70% with $T=1000$. Furthermore the probability for $\widehat n_c = n_c$ or $n_c-1$ is over 92% for $T\ge 500$. Note that when $\widehat n_c = n_c-1$, an effective dimension reduction is still achieved in spite of missing a split.
For the instances with the correct specification, we also calculate the estimation errors for the $n_c$ subspaces. More precisely, put ${A}_{*} = ({A}_{1}, \cdots, {A}_{n_c}) $ and $\widehat {A}_* = (\widehat {A}_{1}, \cdots, \widehat {A}_{n_c})$. We measure the estimation error by
\begin{align*}
D_1 (\widehat {A}_*, {A}_{*}) = \frac{1}{n_c} \sum_{j=1}^{n_c}
\big\{ 1 - \frac{1}{{\rm rank}({A}_{j})} {\rm trace}\big({A}_{j} {A}_{j}^{\top} \widehat{A}_{j} \widehat{A}_{j}^{\top}\big) \big\}.
\end{align*}
Note that $\{{\rm rank}({A}_{j})\}^{-1} {\rm trace}\big({A}_{j} {A}_{j}^{\top} \widehat{A}_{j} \widehat{A}_{j}^{\top}\big)$ is always between 0 and 1. Furthermore it is equal to 1 if ${\cal M}({A}_{j}) = {\cal M}(\widehat {A}_{j})$, and 0 if the two spaces are perpendicular with each other; see, for example, stewart1990 and pan2008. The boxplots of $D_1 (\widehat{A}_{*}, {A}_{*}) $ presented in Figure (ref) indicate that the estimation errors decrease when $T$ increases, and the errors with large $q$ are greater than those with small $q$.
\begin{table}[ht]
\caption{Example 2 -- Relative frequencies correct and incorrect segmentation in a simulation with 1000 replications.
}
\scriptsize
\begin{center}
\begin{tabular}{c|c|cccc||c|cccccccccc}
$T$& $(q,p)$
& Correct & Merging & Splitting & Other & $(q,p)$ & Correct & Merging & Splitting & Other \\\hline
100 & $(6,3)$ & 0.369 & 0.274 & 0.299 & 0.058 & $(6,6)$ & 0.371 & 0.207 & 0.208 & 0.214 \\
500 & & 0.659 & 0.237 & 0.056 & 0.048 & & 0.639 & 0.218 & 0.071 & 0.072 \\
1000 & & 0.721 & 0.220 & 0.031 & 0.028 & & 0.706 & 0.219 & 0.030 & 0.045 \\
5000 & & 0.837 & 0.142 & 0.008 & 0.013 & & 0.815 & 0.156 & 0.011 & 0.018 \\ \hline
100 & $(10,6)$& 0.103 & 0.087 & 0.270 & 0.540 &$(10,10)$& 0.116 & 0.091 & 0.249 & 0.544 \\
500 & & 0.323 & 0.148 & 0.141 & 0.388 & & 0.299 & 0.195 & 0.116 & 0.390 \\
1000 & & 0.369 & 0.227 & 0.086 & 0.318 & & 0.363 & 0.218 & 0.067 & 0.352 \\
5000 & & 0.542 & 0.248 & 0.014 & 0.196 & & 0.495 & 0.269 & 0.019 & 0.217 \\
\end{tabular}
\end{center}
\end{table}
\begin{figure}[ht]
\begin{center}\end{center}
\caption{Example 2 -- boxplots of $D_1(\widehat {A}_*, {A}_{*})$ in a simulation
with 1000 replications.
}
\end{figure}
exampleIn this example we examine the gain in post-sample forecasting due to the decorrelation transformation. Now in model ((ref)) all the elements of ${A}$ and ${B}$ are drawn independently from $U(-1, \, 1)$, with $(p,q)$ equal to $(6,6)$, $(6,8)$, $(8,8)$ and $(10,10)$. For each $(p,q)$ combination, the first three rows (or columns) form one block, the next two rows (or columns) form a second block, and all the other
row (or column) is independent of the rest of the rows (or columns). Table (ref) shows the configuration of ${U}_t$ for $(p,q)=(8,8)$.
\begin{table}
\tiny{
\begin{center}
\begin{tabular}{c|ccc|cc|c|c|c|}
& 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ \hline
1 & \multicolumn{3}{|c|}{\multirow{3}{*}{.81}} & \multicolumn{2}{|c|}{\multirow{3}{*}{.64}} & & & \\
2 & \multicolumn{3}{|c|} & \multicolumn{2}{|c|} & & & \\
3 & \multicolumn{3}{|c|} & \multicolumn{2}{|c|} & & & \\ \hline
4 & \multicolumn{3}{|c|}{\multirow{2}{*}{.25}} & \multicolumn{2}{|c|}{\multirow{2}{*}{.25}} & & & \\
5 & \multicolumn{3}{|c|} & \multicolumn{2}{|c|} & & & \\ \hline
6 & & & & & & & & \\ \hline
7 & & & & & & & & \\ \hline
8 & & & & & & & & \\ \hline
\end{tabular}
\end{center}
\caption{Example 3: Block configuration of ${U}_t$ for the case of $(p,q)=(8,8)$}
}
\end{table}
Each univariate block of ${U}_t$ (e.g. $U_{t,6,6}$) is an AR(1) process with the autoregressive coefficient drawn from the $U[0.1,0.3]$ and
independent and $N(0,1)$ innovations. For each block of size $(3,1)$, $(1,3)$, $(2,1)$ and $(1,2)$, a vector AR(1) model is used. All elements
of the coefficient matrix is drawn independently from $U[0,1]$, then normalized so its singular values are within $[0.1, 0.3]$. The vector innovations consist of independent and N$(0,1)$ random variables. Each of the remaining four blocks follow a matrix AR(1) (i.e. MAR(1)) model of chen2021autoregressive, i.e.
\begin{align}
{U}_{t,i,j}=\lambda_{i,j}{\Phi}_{1,i,j} {U}_{t-1,i,j} {\Phi}_{2,i,j}^\top+{E}_{t,i,j}, \quad
i,j =1, 2,
\end{align}
where ${\Phi}_{1,i,j}$ and ${\Phi}_{2,i,j}$ are the coefficient matrices with unit spectral norm,
all elements of ${E}_{t,i,j}$ are independent and $N(0,1)$. Matrices ${\Phi}_{1,i,j}$ and ${\Phi}_{2,i,j}$ are constructed such that their singular values are all one and the left and right singular matrices are orthonormal, and the scalar coefficient $\lambda_{i,j}$ controls the level of auto-correlation and its values for the four blocks are marked in Table (ref).
We set sample sizes $T=500,1000,2000,5000$, and $(p,q)=(6,6),(8,6),(8,8),(10,10)$. For each setting, we simulate a matrix time series of length $T+M$ with $M=10$. To perform one-step ahead post-sample forecasting, for each $i=0, 1, \cdots, M-1$, we use the first $T$ observations for identifying the uncorrelated blocks and for computing $\widehat {A}_*, \widehat {B}_*$ and $\widehat {U}_t^*$ defined in ((ref)). Depending on its shape/size, we fit each identified block in $\widehat {U}_t^*$ with an MAR(1), VAR(1) or AR(1) model. The fitted models are used to predict ${U}_{T+i+1}^*$. The predicted value for $ X_{T+i+1}$ is then obtained by
align[align omitted — 214 chars of source]
We then compute the mean squared error:
equation[equation omitted — 111 chars of source]
where $\widehat X_{s,i,j}$ denote the one-step ahead predicted value for the $(i,j)$-th element $X_{s,i,j}$ of $ X_s$, and $\mu_{s,i,j} =
\mathbb{E}(X_{s,i,j}| X_{s-1}, X_{s-2}, \cdots)$. We compare $\widehat X_{s,i,j}$ with $\mu_{s,i,j}$, instead of $ X_{s,i,j}$, to remove the impact of the noise in the observed $ X_{s,i,j}$. We report overall the mean and the standard deviation of MSE over 1000 replications.
For the comparison purpose, we also include several other models in our simulation:
(i) VAR(1) model, i.e. $ X_t$ is stacked into a vector of length $pq$ and a VAR(1) is fitted to the vector series directly, (ii) MAR(1) model, i.e. MAR(1) of chen2021autoregressive, is fitted directly to $ X_t$, and (iii) TS-PCA, i.e. all the elements of $ X_t$ are stacked into a long vector and the segmentation method of chang2018 is applied to the vector series.
In addition, we also include two oracle models:
(O1) the true segmentation blocks are used to model $\widehat{U}_t^*$; and
(O2) True $\Sigma^{(x)}_1,\Sigma^{(x)}_2$ and true ${A}_*,{B}_*$ are used and the true segmentation blocks are used to model $\widehat{U}_t^*$.
The mean and standard deviations (in bracket) of the MSEs are listed in Table (ref). It is clear that the forecasting based on the proposed decorrelation transformation is more accurate than those based on MAR(1), VAR(1) and TS-PCA directly. The improvement is often substantial.
The only exception is the case of $(p,q)=(6,6)$ and $T=5000$: the 36-dimension VAR(1) performs marginally better (i.e. the decrease of 0.002 in MSE) than the transformation based method. On the other hand, the improvement of the transformation based forecast over VAR(1) for sample size $T \le 2000$ is more significant. The relative poor performance of MAR(1) is due to the substantial discrepancy between MAR(1) model and the data generating mechanism used in this example. Note that there are $(p+q)^2$ free parameters in VAR(1) while there are merely $(p^2+q^2)$ parameters in MAR(1). With large $T$, VAR(1) is more flexible than MAR(1). TS-PCA performs poorly, due to the {large} error in
estimating the large transformation matrix of size $pq\times pq$.
See also Remark (ref). In addition, TS-PCA also requires to fit VAR models for several large blocks of size 9, 6, 6 and etc.
Table (ref) indicates that the oracle model O1 performs only slightly better than segmentation when sample size $T=500$ or 1000.
Oracle Model O2 only deviates from the true model in terms of the estimation error of the coefficients in the AR models, hence performs the best.
The MSE of O2 is very close to 0, since the only source of error in MSE is the estimation error of the AR models for each blocks of $U_t$, which goes to zero quickly as $T$ increases since these models are of small dimensions. When $T$ is large, the VAR(1), proposed segmentation and O1 all perform about the same, which is worse than O2.
table[table omitted — 2,717 chars of source]
exampleNext, we assess the sensitivity of the proposed matrix segmentation method with respective to the Kronecker product structure ($\hbox{\rm vec}(X_t)=(A\otimes B) \hbox{\rm vec}(U_t)$). We use the same model setting in Example 3 and choose $(p,q)=(8,8), T=500$. Instead of model (ref), we consider $\hbox{\rm vec}(X_t)=(A\otimes B+\Psi) \hbox{\rm vec}(U_t)$, where $\Psi$ is a perturbation matrix with $N$ none-zero elements. We also set each non-zero element of $\Psi$ to be $c \|A\otimes B\|_F/(pq)\times \omega$, $\omega\overset{\rm iid}{\sim}$ Rademacher distribution. Table (ref) shows MSEs and standard deviations of VAR(1) and the proposed matrix segmentation method over 1000 replications, for several combinations of $N$ and $c$. Again, $N$ controls the number of non-zero perturbations and $c$ controls the level of the perturbations, comparing to the average size of the elements in $A\otimes B$. When $c=0$, it becomes the original case in Table (ref). From Table (ref), it is seen that the segmentation is still useful in prediction when the model is deviated from the assumed block structure within certain perturbation. The prediction performance of the VAR model is not affected by the perturbation since it does not assume any block structure. It is seen that the segmentation outperforms the VAR model when $c<0.3$ when $N=10$ and $c<0.2$ when $N=100$. This feature demonstrates the benefit of dimension reduction through segmentation, even when the underlying model does not have the exact Kronecker product structure.
table[table omitted — 974 chars of source]
Real data analysis
We now illustrate the proposed decorrelation method with a real data example. The data concerned is a 17$\times$6 matrix time series consisting of
the logarithmic daily averages of the 6 air pollutant concentration readings (i.e. PM$_{2.5}$, PM$_{10}$, SO$_2$, NO$_2$, CO, O$_3$) from the 17 monitoring stations in Beijing and its surrounding areas (i.e. Tianjin and Hebei Province) in China.
The sampling period is January 1, 2015 -- December 31, 2016 for a total of $T=731$ days.
The readings of different pollutants at different locations are crossly correlated.
We first remove the seasonal mean of the $17\times 6$ matrix time series. Setting $\tau_0=5$ in ((ref)), we apply the proposed bilinear transformation to discover the uncorrelated block structure.
The finding is stable:
with $5\le\tau_1\le 30$ in ((ref)), the transformed matrix series admits $\widehat n_c=4$ uncorrelated column blocks with sizes $3,1,1, 1$ respectively, and $\widehat n_r=15$ uncorrelated row blocks with two blocks of size 2 and the rest of size 1. Overall there are $15\times 4$ blocks, among which there are 2 blocks of size $2\times 3$, 13 blocks of size $1\times 3$, 6 blocks of size $2\times 1$ and 39 blocks of size $1\times 1$.
To check the post-sample forecasting performance, we calculate the rolling one-step, and two-step ahead {out-sample} forecasts for each of the daily readings in the last three months in 2016 (total 92 days). The methods included in the comparison are (i) the forecasting based on the segmentation derived from the proposed decorrelation transformation by fitting an MAR(1), VAR(1) or AR(1) to each identified block according to its size and shape; (ii) MAR(1) for the original $17\times 6$ matrix series; (iii) VAR(1) for the vectorized original series, (iv) univariate AR(1) for each of the $17\times 6$ original series, and (v) the TS-PCA chang2018 for the vectorized original series. The TS-PCA leads to the segmentation consisting of one group of size 3, three groups of size 2, and 93 groups of size 1. {For each group, a VAR(1) model is used for prediction.}
For the two segmentation approaches, we fix the needed transformation obtained using data up to September 30, 2016, and the corresponding segmentation structure. The time series model for each individual block is updated throughout the rolling forecasting period.
In addition to the identified segmentation with 15$\times$4 blocks, we also compute the forecasts based on two other segmentations: one with 14$\times$3 blocks, and another with 16$\times$5 blocks. The former is obtained by merging two most correlated single row blocks in the discovered segmentation into one block and two most correlated single column blocks into one block. The latter is obtained by splitting one of the two blocks with two rows in the discovered segmentation into two single row blocks and splitting the block with 3 columns into two blocks. This will reveal the impact of slightly `wrong' segmentations on forecasting performance.
table[table omitted — 1,145 chars of source]
Let the mean squared forecast error for $s$-th observation be
equation[equation omitted — 112 chars of source]
where $\widehat X_{s,i,j}$ denote the one-step ahead predicted value for the $(i,j)$-th element $X_{s,i,j}$ of $ X_s$. Two-step ahead MSPE$_s$ is defined similarly. The means and the standard deviations of MSPE$_s$ of one-step and two-step ahead post-sample forecasts for the pollution readings in the last 92 days are listed in Table (ref). The prediction based on the bilinear decorrelation transformation (even using `wrong' segmentations) is clearly more accurate than those without transformation as well as that of the TS-PCA. Using the `wrong' segmentation deteriorates the performance only slightly, indicating that a partial (instead of total) decorrelation still leads to significant gain in prediction. Applying TS-PCA to the vectorized series requires to estimate a $102\times 102$ transformation matrix (instead of the $17\times 17$ and $6\times 6$ matrices {by utilizing the matrix time series structure}). It leads to larger estimation errors (see Remark (ref) in Section (ref)). Nevertheless it still provides more accurate predicts than those without transformation.
Also note that fitting each of the original $17\times 6$ time series with a univariate AR(1) separately leads to a better performance
than those from fitting MAR(1) and VAR(1) to the original matrix series jointly. This indicates that the cross-serial correlation is useful and important information for future forecasting. However to make efficient use of the information, it is necessary to adopt some effective independent component analysis or dimension-reduction techniques such as the proposed decorrelation transformation which pushes correlations across different series into the autocorrelations of some transformed series.
Also included in the table are `regret' defined as the difference between MSPE and the in-sample residual variance of the fitted VAR(1) model for the vectorized data (i.e. a vector time series of dimension $17\times 6= 102$), and `adjusted ratio' defined as the ratio of the regret to the regret of the best prediction model (i.e. the segmentation with $15\times 4$ blocks). The regret tries to measure the forecasting error on the predictable signal, removing the impact from unpredictable noise in the model. Since the fitted VAR(1) model uses more than 10,000 parameters,
its sample residual variance underestimates the noise variance in the model, and is taken as a proxy for the latter.
The adjusted ratios shows that the `regret' of TS-PCA is 14.4% worse than the $15\times 4$ segmentation for one-step ahead forecast, and 6.4% worse than the $15\times 4$ segmentation for two-step ahead forecast.
To further evaluate the forecasting stability, Figure (ref) shows weekly average of one-step rolling forecast MSPE within the forecasting period under various models. It is seen that the proposed segmentation outperforms the other three methods in most of the periods in terms of prediction.
figure[figure omitted — 377 chars of source]
In this example, the identified segmentation of $(\widehat n_r, \widehat n_c)= (15,4)$ for the original 17$\times$6 matrix series is more likely to be an approximation {of} the underlying dependence structure rather than a true model. The fact that the two `wrong' segmentation models (though quite close to the discovered one) as well as the aggressive segmentation via vectorized TS-PCA transformation provide comparable, though inferior, post-sample forecasting performance lends further support to the claim that the proposed decorrelation method makes the transformed series more predictable than the original ones, regardless model ((ref)) holds or not. See Remark (ref)(iv) in Section (ref) above.