text:
\section{Introduction}

Data recorded in a sequence of matrices are increasingly common in this information age,
due to the demand to solve practical problems from, among others, signal processessing,
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 \cite{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. Consequently
there is a growing interest in analyzing matrix time series,
see, for example, a matrix factor model based on the Tucker decomposition
for tensors \citep{wang2019, chen2019factor,
chen2019constrained, chen2019matrix, han2020, chen2020modeling,han2020number},
a bilinear autoregressive model \citep{chen2020autoregressive,xiao2021}, and
a matrix factor model based on the tensor CP decomposition
\citep{chy2021}.

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 (i.e. VARMA) models. Note that a 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 \citep{basu2015, guo2016high, lin2017,
hall2018, gao2019, zhou2018, ghosh2019, lin2020regularized, han2020sparse},
(ii) factor models of various forms \citep{pena1987identifying, bai2002, forni2005generalized, lam2012}, (iii) various versions of independent
component analysis \citep{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 dimension-reduction method for modelling and forecasting
matrix time series. Our proposal is in the spirit of independent
component analysis. It transform a $p\times q$ matrix series into the new 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 many
submatrix series and those submatrix series are uncorrelated across all time lags.
Hence we achieve an effective dimension-reduction 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: 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 dimension model is typically not enough to
offset the error resulted from estimating the large number of parameters.
This is why forecasting a large number of time series together based on
a joint time series model can be worse than that forecasting each component series
separately by ignoring cross serial correlations completely. See \cite{chang2018}
and also
Section \ref{section:data} below.
The proposed transformation channels 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
effectively in future 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. Empirical evidence also
indicates that the decorrelation improves the inference efficiency \citep{chang2017testing}.

The basic idea of the proposed transformation is similar to the so-called principal
component analysis for time series (PCA-TS) of \cite{chang2018}. Actually one might
be tempted to stack a $p\times q$ matrix series into a $(pq)\times 1$ vector time series,
and apply PCA-TS directly. This requires to search for a $(pq)\times (pq)$ transformation
matrix, which is almost prohibited for $p$ and $q$ of the size a few
thousands. In contrast, the proposed bilinear transformation is facilitated by two
matrices of size, respectively, $p\times p$ and $q\times q$.
Perhaps more importantly 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
\ref{section:data}).

The rest of the paper is organized as follows. The proposed model is presented
in Section \ref{sec2}. We introduce an new normalization under which
the required bilinear transformation consists of two orthogonal transformations.
The orthogonality allows us to accumulate the information from different time lags without cancellation.
The estimation of the proposed 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{section:theories}. More precisely the non-asymptotic error bounds
of the estimated transformation are derived based on some concentration
inequalities (e.g., \cite{rlo2000theorie,merlevede2011bernstein}); further leading to the uniform convergence rates of the estimation.
The consistency of the proposed latent uncorrelated blocks selection method is also presented. Numerical illustration with both simulated and real data is reported in
Section \ref{sec4}, which shows superior performance in forecasting over the methods
without transformation. Once again it reinforces the fact that the across serial
correlations are important and useful information for future prediction for a large
number of time series, and, however, it is necessary to adopt appropriate
dimension-reduction techniques, such as the proposed decorrelation transformation,
in order to use the information effectively \citep{chang2018}.
All technical proofs are relegated to an Appendix.

\section{Settings and methodology} \label{sec2}

\subsection{Model} \label{sec21}

Let $X_t= (X_{t,i,j})$ be a $p\times q$ matrix time series, i.e. there are
$pq$ recorded values at each time from, for example, $p$ individuals
and over $q$ indices or variables.
We assume that $ X_t$ is weakly stationary in the sense that all the first
two moments are finite and time-invariant.
With available observations $ X_1, \cdots,  X_T$,
the ultimate goal is to build a dynamic model for $ X_t$ and to forecast the future
values $ X_{T+1},  X_{T+2}, \cdots.$ With moderately large $p$ and $q$,
any direct attempts based on ARMA framework are unlikely to be successful
due to overparametrization.
We assume that $ X_t$ admits a latent
segmentation structure
\begin{equation} \label{a1}
 X_t = {\utwi{B}} {\utwi{U}}_t {\utwi{A}}^\top,
\end{equation}
where ${\utwi{A}}$ and ${\utwi{B}}$ are unknown
$q\times q$ and $p\times p$, respectively, invertible constant matrices,
and random matrix ${\utwi{U}}_t$ is unobservable and is of a block structure
\begin{equation} \label{b1}
{\utwi{U}}_t = \left(
\begin{array}{cccc}
{\utwi{U}}_{t,1,1} & {\utwi{U}}_{t, 1, 2}& \cdots & {\utwi{U}}_{t, 1, n_c} \\
{\utwi{U}}_{t,2,1} & {\utwi{U}}_{t, 2, 2}& \cdots & {\utwi{U}}_{t, 2, n_c} \\
\vdots & \vdots & \vdots &\vdots\\
{\utwi{U}}_{t,n_r,1} & {\utwi{U}}_{t, n_r, 2}& \cdots & {\utwi{U}}_{t, n_r, n_c}
\end{array}
\right)
\end{equation}
in which ${\utwi{U}}_{t, i, j}$ is a $p_i\times q_i$ matrix with unknown $p_i$ and $q_i$,
$\sum_{i=1}^{n_r}p_i=p$, and $\sum_{j=1}^{n_c}q_j=q$.  Furthermore,
${\utwi{U}}_{t+\tau,i_1,j_1}$ and ${\utwi{U}}_{t,i_2,j_2}$, for any $(i_1, j_1)\neq (i_2, j_2)$,
 are uncorrelated for any integer $\tau$.
Note that the rows and the columns of ${\utwi{U}}_t$ are linear combinations of,
respectively, the rows and the columns of $ X_t$.

One would be tempted to stack all the elements of $ X_t$ into a long vector
and to apply PCA-TS of \cite{chang2018} directly. This requires to search for
the $(pq) \times (pq)$ decorrelation transformation matrix ${\utwi{\mathnormal\Phi}}$  such that
$\hbox{\rm vec}( X_t) = {\utwi{\mathnormal\Phi}} \hbox{\rm vec}({\utwi{U}}_t)$. The proposed model (\ref{a1}) is effectively
to impose a low-dimension structure ${\utwi{\mathnormal\Phi}}={\utwi{A}}\otimes
{\utwi{B}}$, where $\otimes$ denotes the matrix Kronecker product.
Also note that simpler models are more likely to be applicable to smaller
matrix series. For example, the bilinear AR models of
\citet{chen2020autoregressive} with small
orders are more plausible for ${\utwi{U}}_{t,i,j}$ than for $ X_t$, the factor models of
\cite{wang2019} or \cite{chy2021} with small number of factors are more
appropriate for ${\utwi{U}}_{t,i,j}$ than for $ X_t$.

The structure in (\ref{b1}) may be too rigid in the sense that both $n_r$ and $n_c$ may be
small or very small in practice. One can apply again the proposed method to each submatrix
series $U_{t, i,j}$; leading to a finer segmentation with irregularly sized small blocks.

\subsection{Normalization and identification} \label{22}
To simplify the statement, 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}({\utwi{U}}_t) = 0$.

The terms ${\utwi{A}}, {\utwi{B}}$ and ${\utwi{U}}_t$ on the RHS of (\ref{a1}) are not uniquely defined,
and any ${\utwi{A}}, {\utwi{B}}$ and ${\utwi{U}}_t$ satisfying the required condition will serve the
dimension-reduction purpose.
Therefore we can take the advantage from this lack of uniqueness to identify the `ideal'
${\utwi{A}}$ and ${\utwi{B}}$ to improve the estimation effectiveness. More precisely, by applying a
new normalization, we identify ${\utwi{A}}$ and ${\utwi{B}}$ to be orthogonal,
which facilitates the accumulation of the information from different lags
without cancellation. To this end, we first
list some basic facts as follows.
\begin{enumerate}
\item
 $({\utwi{B}}, {\utwi{U}}_t, {\utwi{A}}^\top)$ in (\ref{a1}) can be replace
by $({\utwi{B}} {\utwi{D}}_1, {\utwi{D}}_1^{-1} {\utwi{U}}_t {\utwi{D}}_2^{-1}, {\utwi{D}}_2 {\utwi{A}}^\top)$ for any
invertible block diagonal matrices ${\utwi{D}}_1$ and ${\utwi{D}}_2$ with the block sizes,
respectively, $(p_1, \cdots, p_{n_r})$ and $(q_1, \cdots, q_{n_c})$.
\item
The $q$ columns of ${\utwi{B}} {\utwi{U}}_t$ can be divided into $n_c$ uncorrelated blocks of
sizes $(q_1, \cdots, q_{n_c})$, i.e. the columns of ${\utwi{B}} {\utwi{U}}_t$ resemble the same pattern
of the uncorrelated blocks
as those of ${\utwi{U}}_t$. Thus
${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(u)} \equiv \mathbb{E}({\utwi{U}}_t^\top{\utwi{B}}^\top {\utwi{B}} {\utwi{U}}_t)/p$ is a block diagonal matrix
with the block sizes $(q_1, \cdots, q_{n_c})$. In the same vein,
${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(u)} = \mathbb{E}({\utwi{U}}_t{\utwi{A}} {\utwi{A}}^\top {\utwi{U}}_t^\top)/q$ is a block diagonal matrix with
the block sizes $(p_1, \cdots, p_{n_r})$.
\item Put ${\utwi{B}}=({\utwi{B}}_1, \cdots, {\utwi{B}}_{n_r})$ and ${\utwi{A}}= ({\utwi{A}}_1, \cdots, {\utwi{A}}_{n_c})$,
where ${\utwi{B}}_i$ is $p_i\times q$ and ${\utwi{A}}_i$ is $q_i\times  p$. Then linear spaces
${\cal M}({\utwi{B}}_1), \cdots, {\cal M}({\utwi{B}}_{n_r})$ and ${\cal M}({\utwi{A}}_1), \cdots,
{\cal M}({\utwi{A}}_{n_c})$ are uniquely defined by (\ref{a1}), though ${\utwi{B}}$ and ${\utwi{A}}$ are not,
where ${\cal M}({\utwi{G}})$ denotes the linear space spanned by the columns of matrix ${\utwi{G}}$.
Note that for any invertible $q\times q$ matrix ${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1$ and $p\times
p$ matrix ${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2$,
\[
{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1 {\utwi{A}} = ({\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1 {\utwi{A}}_1, \cdots, {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1 {\utwi{A}}_{n_c}), \qquad {\cal M}({\utwi{A}}_i) =
{\cal M}({\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1 {\utwi{A}}_i) \;\;\; i=1, \cdots, n_c,
\]
\[
{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2 {\utwi{B}} = ({\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2 {\utwi{B}}_1, \cdots, {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2 {\utwi{B}}_{n_r}), \qquad {\cal M}({\utwi{B}}_i) =
{\cal M}({\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2 {\utwi{B}}_i) \;\;\; i=1, \cdots, n_r.
\]
\end{enumerate}
Put
${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)} = \mathbb{E}( X_t^\top X_t)/p$ and
${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)} = \mathbb{E}( X_t X_t^\top)/q$.
Proposition \ref{prop1} below indicates that if we replace $ X_t$  in
(\ref{a1}) by a normalized
version $\{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-1/2}  X_t \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-1/2}$, we can treat
both ${\utwi{A}}$ and ${\utwi{B}}$ in (\ref{a1}) 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{prop2} in the appendix).

\begin{proposition} \label{prop1}
Let both ${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}$ and ${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}$
be invertible. Then model (\ref{a1}) can be written as
\begin{equation} \label{an1}
 X_t = \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{1/2} {\utwi{B}}_* {\utwi{U}}_t^* {\utwi{A}}_*^\top \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{1/2},
\end{equation}
where $ {\utwi{U}}_t^*$ admits the same decomposition as (\ref{b1}) with
$p\times q$ uncorrelated blocks, and ${\utwi{A}}_*$ and $ {\utwi{B}}_*$ are, respectively,
$q\times q$ and $p\times p$
orthogonal matrices. More precisely,
\begin{equation} \label{an2}
{\utwi{U}}_t^* = \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(u)}\}^{-1/2} {\utwi{U}}_t \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(u)}\}^{-1/2}, \;\;
{\utwi{A}}_* = \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-1/2} {\utwi{A}} \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(u)}\}^{1/2}, \;\;
{\utwi{B}}_* =\{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-1/2} {\utwi{B}} \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(u)}\}^{1/2}.
\end{equation}
\end{proposition}

The proof of Proposition \ref{prop1} is almost trivial.
First note that both ${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(u)}, \, {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(u)}$ are invertible.
Then (\ref{an1}) follows from
(\ref{a1}) and (\ref{an2}) directly. The orthogonality of, e.g., ${\utwi{B}}_*$ follows
from equality $E[\{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-1/2}  X_t  X_t^\top \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-1/2}]/q={\utwi{I}}_p$, (\ref{an1}) and (\ref{an2}).  Also note that $\{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(u)}\}^{-1/2}$
and $\{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(u)}\}^{-1/2}$ are of the same block diagonal structure as, respectively,
${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(u)}$ and ${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(u)}$. This implies that ${\utwi{U}}_t^*$
admits the decomposition in the form of (\ref{b1}).

Note that $({\utwi{B}}_*, {\utwi{U}}_t^*, {\utwi{A}}_*^\top)$ in (\ref{an1}) are (still) not uniquely defined,
as they can be replaced by $({\utwi{B}}_* {\utwi{D}}_1, {\utwi{D}}_1^\top{\utwi{U}}_t^* {\utwi{D}}_2^\top, {\utwi{D}}_2{\utwi{A}}_*^\top)$ for
any orthogonal block diagonal matrices ${\utwi{D}}_1$ and ${\utwi{D}}_2$ with appropriate block sizes.
Proposition \ref{prop2} below shows that we can take
the orthonormal eigenvectors of two properly defined positive definite matrices
as the columns of ${\utwi{A}}_*$ and ${\utwi{B}}_*$.  With
${\utwi{A}}_*$ and ${\utwi{B}}_*$ specified, the segmented ${\utwi{U}}_t^*$ can be solved from (\ref{an1})
directly.  Let
\begin{eqnarray*}
{\utwi{V}}^{(1)}_{\tau,i,j}
= \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-1/2}  \mathbb{E}({ X}_{t+\tau}^\top {\utwi{E}}_{i,j}{ X_t})
\{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-1/2}, \quad
{\utwi{V}}^{(2)}_{\tau,i,j}
= \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-1/2}  \mathbb{E}({ X}_{t+\tau} {\utwi{E}}_{i,j}{ X_t}^\top)
\{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-1/2},
\end{eqnarray*}
where ${\utwi{E}}_{i,j} $
is the unit matrix with
1 at position $(i,j)$ and 0 elsewhere, and ${\utwi{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$, put
\begin{equation} \label{aa9}
{\utwi{W}}^{(1)} = \sum_{\tau=-\tau_0}^{\tau_0} \sum_{i=1}^p\sum_{j=1}^p
\frac{{\utwi{V}}^{(1)}_{\tau,i,j}\big({\utwi{V}}^{(1)}_{\tau,i,j}\big)^\top}{p^2}, \qquad
{\utwi{W}}^{(2)} = \sum_{\tau=-\tau_0}^{\tau_0} \sum_{i=1}^q\sum_{j=1}^q
\frac{{\utwi{V}}^{(2)}_{\tau,i,j}\big({\utwi{V}}^{(2)}_{\tau,i,j}\big)^\top}{q^2}.
\end{equation}

\begin{proposition} \label{prop2}
Let both ${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}$ and ${\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}$ be invertible, and all the eigenvalues of
${\utwi{W}}^{(i)}$ be distinct, $i=1, 2$. Then the $q$ orthonormal eigenvectors of ${\utwi{W}}^{(1)}$
can be taken as the columns of ${\utwi{A}}_*$, and the $p$ orthonormal eigenvectors of ${\utwi{W}}^{(2)}$
can be taken as the columns of ${\utwi{B}}_*$.
\end{proposition}

The proposition above does not give any indication on how to arrange the columns
of ${\utwi{A}}_*$ and ${\utwi{B}}_*$,
which should be ordered according to the latent uncorrelated block structure (\ref{b1}).
We address this issue in Section \ref{section:estimation} below.

\begin{remark}\label{rmk:trans}
In case that ${\utwi{W}}^{(i)}$ has tied eigenvalues and the corresponding eigenvectors
across different blocks, Proposition \ref{prop2} no longer holds.
To avoid the tied eigenvalues, we may use different
values of $\tau_0$, or even different forms of ${\utwi{W}}^{(i)}$ if necessary.
For example, we may replace ${\utwi{W}}^{(1)}$ by
\begin{eqnarray*}
{\utwi{W}}^{(1,f)} = \sum_{\tau=-\tau_0}^{\tau_0} \sum_{i=1}^p\sum_{j=1}^p
\frac{f\left( {\utwi{V}}^{(1)}_{\tau,i,j} {\utwi{V}}^{(1)\top}_{\tau,i,j} \right)}{p^2}.
\end{eqnarray*}
where $f({\utwi{V}})$ is a function of a symmetric matrix ${\utwi{V}}$ in which
$f({\utwi{V}})={\utwi{\mathnormal\Gamma}}{\utwi{D}}^{(f)}{\utwi{\mathnormal\Gamma}}^\top$, ${\utwi{V}}={\utwi{\mathnormal\Gamma}}{\utwi{D}}{\utwi{\mathnormal\Gamma}}^\top$ is
the eigen-decomposition for ${\utwi{V}}$, and
${\utwi{D}}^{(f)}=\hbox{\rm diag}(f(d_1),\ldots,f(d_q))$ is the diagonal-element-wise
transformation of the diagonal matrix ${\utwi{D}}$ of the eigenvalues.
\end{remark}

\begin{remark}\label{rmk:tau0}
In practice, we need to specify integer $\tau_0$. In principle any
$\tau_0$ can be used to estimate ${\utwi{A}}_*$ and ${\utwi{B}}_*$,
as  ${\utwi{V}}^{(k)}_{\tau,i,j}\big({\utwi{V}}^{(k)}_{\tau,i,j}\big)^\top$ is non-negative definite,
there is no information cancellation from different time lags.
Since autocorrelation is typically at its strongest at small lags, a
relatively small $\tau_0$, such as $\tau_0=5$, is often sufficient
\citep{lam2011,chang2018,wang2019,chen2019factor}.
\end{remark}

\subsection{Estimation} \label{section:estimation}

The estimation for ${\utwi{A}}_*$ and ${\utwi{B}}_*$ is based on the eigenanalysis of the sample
versions of matrices defined in (\ref{aa9}). To this end, let
\begin{equation} \label{est-0}
\widehat{{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}}_1^{(x)}=\frac{1}{Tp}\sum_{t=1}^T X_t^\top X_t, \qquad
\widehat{{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}}_2^{(x)}=\frac{1}{Tq}\sum_{t=1}^T X_t X_t^\top,
\end{equation}
\begin{equation} \label{est-1}
\widehat{\utwi{V}}^{(1)}_{\tau,i,j}
= \{ \widehat{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-\frac{1}{2}} \sum_{t=1}^{T-\tau}
{{ X}_{t+\tau}^\top {\utwi{E}}_{i,j}{ X_t}
\over T-\tau}
\{ \widehat{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-\frac{1}{2}}, \;\;
\widehat{\utwi{V}}^{(2)}_{\tau,i,j}
= \{ \widehat{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-\frac{1}{2}} \sum_{t=1}^{T-\tau} {{ X}_{t+\tau} {\utwi{E}}_{i,j}{ X_t}^\top \over T-\tau}
\{ \widehat{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-\frac{1}{2}},
\end{equation}
\begin{equation} \label{est-2}
\widehat{\utwi{W}}^{(1)} = {1 \over p^2} \sum_{\tau=-\tau_0}^{\tau_0} \sum_{i=1}^p\sum_{j=1}^p
{\widehat{\utwi{V}}^{(1)}_{\tau,i,j}\big(\widehat{\utwi{V}}^{(1)}_{\tau,i,j}\big)^\top}, \qquad
\widehat{\utwi{W}}^{(2)} = {1 \over q^2}\sum_{\tau=-\tau_0}^{\tau_0} \sum_{i=1}^q\sum_{j=1}^q
{\widehat{\utwi{V}}^{(2)}_{\tau,i,j}\big(\widehat{\utwi{V}}^{(2)}_{\tau,i,j}\big)^\top}.
\end{equation}
Performing the eigenanalysis for $\widehat{\utwi{W}}^{(1)}$, and arranging the
resulting $q$ orthonormal eigenvectors by the algorithm below, we take the
ordered eigenvectors as the columns of $\widehat {\utwi{A}}_*$. The estimator $\widehat {\utwi{B}}_*$ is obtained
in the same manner via the eigenanalysis for $\widehat{\utwi{W}}^{(2)}$. Then by
(\ref{an1}), we obtain the transformed matrix series
\begin{equation} \label{est-3}
\widehat {\utwi{U}}_t^* = \widehat {\utwi{B}}_*^\top\{ \widehat{{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}}_2^{(x)}\}^{-1/2}  X_t
\{\widehat{{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}}_1^{(x)}\}^{-1/2} \widehat {\utwi{A}}_*.
\end{equation}
Note that the estimation for ${\utwi{A}}_*$ and that for ${\utwi{B}}_*$ are carried out separately.
They do not interfere with each other.

Now we present an algorithm to arrange the order of the columns for $\widehat {\utwi{A}}_*$.
By (\ref{an1}), the columns of
${\utwi{Z}}_t \equiv  X_t \{ {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-1/2} {\utwi{A}}_*$
are divided into the $n_c$ uncorrelated blocks. Define
\begin{equation} \label{a12}
\widehat {\utwi{Z}} \equiv  (\widehat z_{t,i,j}) =  X_t \{ \widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-1/2} (\widehat {\utwi{\mathnormal\gamma}}_1,
\cdots, \widehat {\utwi{\mathnormal\gamma}}_q),
\end{equation}
where $\widehat {\utwi{\mathnormal\gamma}}_1, \cdots, \widehat {\utwi{\mathnormal\gamma}}_q$ are the $q$ orthonormal eigenvectors
of $\widehat {\utwi{W}}^{(1)}$.  We divide the columns of $\widehat {\utwi{Z}} $ into some uncorrelated
blocks according to the maximum pairwise correlations between the columns,
i.e. the correlation between the $k$-th and $\ell$-th column is measured by
\begin{align} \label{a13}
\widehat \rho_{k,\ell} & = \max_{1\le i,j \le p, \, |\tau|\le \tau_1}
\big| \widehat{\rm Corr} (\widehat z_{t+\tau, i,k}, \, \widehat z_{t,j, \ell} ) \big|\\
& = \max_{1\le i,j \le p, \, |\tau|\le \tau_1}
\frac{\big| \widehat{ \gamma}_k^\top \widehat{ V}^{(1)}_{\tau,i,j} \widehat{ \gamma}_\ell\big|}
{\big\{\widehat{ \gamma}_k^\top \widehat{ V}^{(1)}_{0,i,i} \widehat{ \gamma}_k \,
        \widehat{ \gamma}_\ell^\top\widehat{ V}^{(1)}_{0,j,j}
\widehat{ \gamma}_\ell\big\}^{1/2}}, \quad 1 \le k < \ell \le q, \notag
\end{align}
The second equality above follows from \eqref{est-1}, \eqref{est-2} and \eqref{a12}.
In the above expression, $\tau_1 \ge 1$ is a sufficiently large integer, for example,
between 10 and 20, in the spirit of the rule of thumb of \cite{box1970}.
Instead of conducting a multiple test for all $q_0=q(q-1)/2$ pair
columns, we propose a ratio-based statistic to single out those
significantly correlated pairs. To this end, we 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)}$. Define
\begin{align}\label{a14}
\widehat r = \arg\max_{1 \le j \le 0.75 q_0} \widehat \rho_{(j)}/ \widehat \rho_{(j+1)}.
\end{align}
We take the $\widehat r$ pair columns corresponding to $\widehat \rho_{(1)}, \cdots,
\widehat \rho_{(\widehat r)}$ as correlated pairs.
The intuition behind this approach is as follows. Suppose among in total
$q_0$ pair columns of ${\utwi{Z}}_t$ there are $r$ cross correlated pairs only.
Arrange the true maximum absolute cross correlations in the descending order:
$\rho_{(1)} \ge \cdots \ge  \rho_{(q_0)}$. Then $\rho_{(r)}>0$, $\rho_{(q_0)}=0$,
and the ratio $\rho_{(j)}/ \rho_{(j+1)}$ takes value $\infty$ at $j=r$.
See \cite{chang2018} for the application of the same idea in TS-PCA.
Similar ideas
have also been used in determining the number of factors in \cite{lam2012}, \cite{ahn2013} and \cite{han2020number}.
Instead of using 0.75 in
\eqref{a14}, we may choose $c_r\in(0,1)$, to avoid the
fluctuations due to the ratios of extremely small values.

Once we have identified $\widehat r$ correlated pair columns of $\widehat{\utwi{Z}}_t$,
the grouping can
be performed as follows:
\begin{quote}
(i) Start with $q$ groups with one column in each group.\\
(ii) Merge two groups together if two correlated columns are split over
the two groups.\\
(iii) Repeat Step (ii) above until all the correlated columns are
within one group. The number of the final groups is the estimated value for $n_c$.
\end{quote}

The finite sample performance of the above algorithm
 can be improved by prewhitening each column time series
of $\widehat {\utwi{Z}}$. This makes $\widehat \rho_{k,\ell}$, for different $(k, \ell)$,
more comparable. See Remark 2(iii) of \cite{chang2018}. In practice the
prewhitening can be carried out by fitting each column time series an VAR
model with the order between 0 and 5 determined by AIC. The resulting residual series is taken as a prewhitened series.

The ordering for columns of $\widehat {\utwi{B}}_*$ is arranged, in the same manner
as above, by examining the pairwise correlations among the columns of
\[
 X_t^\top \{ \widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-1/2} (\widehat {\utwi{\mathnormal\gamma}}_1,
\cdots, \widehat {\utwi{\mathnormal\gamma}}_p),
\]
where $\widehat {\utwi{\mathnormal\gamma}}_1,
\cdots, \widehat {\utwi{\mathnormal\gamma}}_p$ are now the $p$ orthonormal eigenvectors of $\widehat {\utwi{W}}^{(2)}$.

\section{Theoretical properties}\label{section:theories}

To gain more appreciation of the methodology, we will show the consistency of the proposed detection method of uncorrelated components.

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 ${\utwi{V}}=(V_{ij})$, let $\|{\utwi{V}}\|_{op}$ denote the spectrum norm,
and $\|{\utwi{V}}\|_{\max} = \max_{i,j} |V_{ij}|$. To facilitate the theoretical analysis, we impose the following assumptions.

\begin{assumption} \label{asmp:mixing}
Assume exponentially decay coefficients of the strong $\alpha$-mixing condition,
\begin{eqnarray*}
\alpha(k) \le \exp\big( - c_0 k^{r_1} \big)
\end{eqnarray*}
for some constant $c_0>0$ and $0<r_1\le 1$, where
\begin{equation}\label{def:mixing}
\alpha(k) = \sup_t\Big\{\Big|\mathbb{P}(A\cap B) - \mathbb{P}(A)\mathbb{P}(B)\Big|:
A\in \sigma( X_s, s\le t), B\in \sigma( X_s, s\ge t+k)\Big\}.
\end{equation}
\end{assumption}

Here for any random variable/vector/matrix $ X$, $\sigma( X)$ is understood to be the $\sigma$-field generated by $ X$. Assumption \ref{asmp:mixing} allows a very general class of time series models, including causal ARMA processes with continuously distributed innovations; see also \cite{tong1990non, bradley2005, tsay2005analysis, fan2008nonlinear, rosenblatt2012markov, tsay2018nonlinear}, among others. The restriction $r_1\le 1$ is introduced only for presentation convenience.

As $|\mathbb{P}(A\cap B) - \mathbb{P}(A)\mathbb{P}(B)|\le 1/4 \le e^{-1}$ for all event pairs $\{A,B\}$,
\begin{eqnarray*}
&& \sup_{t}\Big\{\Big|\mathbb{P}(A\cap B) - \mathbb{P}(A)\mathbb{P}(B)\Big|:
A\in \sigma( X_{s+\tau}\otimes  X_s, s\le t,\tau\le\tau_0),
B\in \sigma( X_{s+\tau}\otimes  X_s, s\ge t+k)\Big\}
\cr &\le & \sup_{t}\Big\{\Big|\mathbb{P}(A\cap B) - \mathbb{P}(A)\mathbb{P}(B)\Big|:
A\in \sigma( X_s, s\le t+\tau_0),
B\in \sigma( X_s, s\ge t+k)\Big\}
\cr &\le & \exp\Big[- \max\{1,c_0(k-\tau_0)^{r_1}\}\Big]
\cr &\le & \exp\Big[ - c_0(1+c_0^{1/r_1}\tau_0)^{-r_1}k^{r_1}\Big],
\end{eqnarray*}
where $\otimes$ denotes the Kronecker product for matrices.
Thus, the quadratic process $( X_{t+\tau}\otimes  X_t, 0\le \tau\le\tau_0)$
also satisfies the $\alpha$-mixing condition with
$c_0$ in Assumption \ref{asmp:mixing} replaced by the smaller $c_0(1+c_0^{1/r_1}\tau_0)^{-r_1}$.

\begin{assumption} \label{asmp:eigenvalue}
There exist certain finite constants $c_*, c^*$, such that
\begin{equation}
\sup_{i\le p, \|{\mathbf{u}}\|_2=1} \textsf{var}\Big( X_{t,i,\cdot}^\top {\mathbf{u}}\Big)^2 \le c^*,\quad
\inf_{\|{\mathbf{u}}\|_2=1}\mathbb{E} \| X_t{\mathbf{u}} \|_2^2/p\ \ge\ c_*.
\end{equation}
\end{assumption}
Note that $\sup_{ \|{\mathbf{u}}\|_2=1}\textsf{var}\Big( X_{t,i,\cdot}^\top {\mathbf{u}}\Big)^2= \|\Sigma^{(x)}_{0,i,i} \|_{op}$ and $\mathbb{E} \| X_t{\mathbf{u}} \|_2^2/p=u^\top {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}u$. Assumption \ref{asmp:eigenvalue} on the eigenvalues is a common assumption in the high dimensional setting, for instance, \cite{bickel2008covariance,bickel2008regularized, xia2017}.
It rules out the possibility of latent factor structure in $ X_t$. As the ultimate goal is forecasting, for ultra high dimensional matrix time series, we suggest to de-correlate the estimated latent processes after doing factor analysis first.

\begin{assumption}\label{asmp:tail}
For any $x>0$,
\begin{equation*}
\max_t\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) ,
\end{equation*}
for some constant $c_1,c_2>0$ and $0<r_2\le 2$.
\end{assumption}

Assumption \ref{asmp:tail} 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.

The following theorem provides the non-asymptotic bounds for $\widehat{ \Sigma}^{(1)}_{\tau,i,j}:=\sum_{t=1}^{T-\tau} X_{t+\tau}^\top{\utwi{E}}_{ij} X_t/(T-\tau)$, $\widehat{ V}^{(1)}_{\tau,i,j}$ and $\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1)}$, $1\le i,j\le q$ under both exponential decay $\alpha$ mixing condition and exponential tail condition of $ X_t$.

\begin{theorem}\label{th-1new}
Suppose Assumptions \ref{asmp:mixing}, \ref{asmp:eigenvalue}, \ref{asmp:tail} 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{eqnarray}\label{th-1-1new}
\Big\|\widehat{ \Sigma}^{(1)}_{\tau,i,j} - {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}^{(1)}_{\tau,i,j}\Big\|_{op}
\le q\Big\|\widehat{ \Sigma}^{(1)}_{\tau,i,j} - {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}^{(1)}_{\tau,i,j}\Big\|_{\max}
\le \eta^*_{T,p,q}, \;\; \forall 1\le \tau\le \tau_0, 1\le i,j\le p
\end{eqnarray}
in an event $\Omega_T$ with probability at least $1 -\eps_T$, where
\begin{eqnarray*}
\eta^*_{T,p,q} = C_1q\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*}
for some positive constant $C_1$ depending on $r_1,r_2$. Moreover, if $\eta^*_{T,p,q} \le c_*/2$, then
\begin{eqnarray}\label{th-1-2new}
&& \big\|\widehat{ V}^{(1)}_{\tau,i,j} - {\utwi{V}}^{(1)}_{\tau,i,j}\big\|_{op} \le \eta^{**}_{T,p,q}, \qquad
\big\|\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1)}  - {\utwi{W}}^{(1)} \big\|_{op} \le \eta^{**}_{T,p,q},
\end{eqnarray}
in the same event $\Omega_T$ with probability at least $1-\epsilon_T$, where
$\eta^{**}_{T,p,q} = C_2\eta^*_{T,p,q}$ for a constant $C_2$ depending on $c_*,c^*$ only.
\end{theorem}

When the decays of the $\alpha$ mixing coefficients and the tail probabilities of $ X_t$ are slower than Assumptions \ref{asmp:mixing} and \ref{asmp:tail}, we impose the following Assumptions \ref{asmp:mixing2} and \ref{asmp:tail2}. These conditions ensure the Fuk-Nagaev-type inequalities for $\alpha$-mixing processes; see also \cite{rlo2000theorie}, \cite{liu2013probability, wu2016performance} and \cite{zhang2017gaussian}.

\begin{assumption} \label{asmp:mixing2}
Assume algebraically decay coefficients of the strong $\alpha$-mixing condition
\begin{eqnarray*}
\alpha(k) \le c_0 k^{-r_1},
\end{eqnarray*}
where $\alpha(k)$ is defined in \eqref{def:mixing}, $c_0>0$ and $r_1>1$.
\end{assumption}

\begin{assumption}\label{asmp:tail2}
For any $x>0$,
\begin{equation*}
\max_t\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} ,
\end{equation*}
for some constant $c_1>0$ and $r_2>1$.
\end{assumption}

Theorem 2 below presents the uniform convergence rate for $\widehat{ \Sigma}^{(1)}_{\tau,i,j}$, $\widehat{ V}^{(1)}_{\tau,i,j}$ and $\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1)}$, $1\le i,j\le q$. It is intuitively clear that the rate is much slower than that in Theorem \ref{th-1new}.

\begin{theorem}\label{th-1-poly}
Suppose Assumptions \ref{asmp:eigenvalue}, \ref{asmp:mixing2}, \ref{asmp:tail2} 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, \eqref{th-1-1new} holds in an event $\Omega_T$ with probability at least $1-\epsilon_T$, where
\begin{eqnarray*}
\eta^*_{T,p,q} = C_1q\left( \frac{(Tpq/\epsilon_T)^{1/\beta_3}}{T} + \sqrt{\frac{\log(pq/\epsilon_T)}{T}} \right) ,
\end{eqnarray*}
for some positive constant $C_1$ depending on $r_1,r_2$. Similarly, if $\eta^*_{T,p,q} \le c_*/2$, then \eqref{th-1-2new} holds in the same event $\Omega_T$.
\end{theorem}

Write the eigenvalue decomposition of ${\utwi{W}}^{(1)}$ as
\begin{eqnarray*}
{\utwi{W}}^{(1)} = {\utwi{\mathnormal\Gamma}}^{(1)}{\utwi{D}}\big({\utwi{\mathnormal\Gamma}}^{(1)}\big)^\top,
\end{eqnarray*}
where ${\utwi{D}}=\hbox{\rm diag}(\lambda_1,\ldots,\lambda_q)$ gives the eigenvalues $\lambda_1\ge
\cdots \ge \lambda_q$ and
${\utwi{\mathnormal\Gamma}}^{(1)} = ({\utwi{\mathnormal\gamma}}^{(1)}_1,\ldots,{\utwi{\mathnormal\gamma}}^{(1)}_q)$ gives the corresponding eigenvectors.
For any group $G\subset \{1,\ldots,q\}$, define the eigengap as
\begin{eqnarray}\label{eq:deltag}
\Delta_G = \min\Big\{|\lambda_j-\lambda_k|: j\in G, k\not\in G\Big\}.
\end{eqnarray}
Let ${\utwi{\mathnormal\Gamma}}_{G}^{(1)}=({\utwi{\mathnormal\gamma}}_j^{(1)})_{j\in G}$. We say that the group $G$ is {\it identifiable} if $\Delta_G>0$ as the column space of ${\utwi{\mathnormal\Gamma}}^{(1)}_G$ would then be the span of $|G|$ eigenvectors for all versions of ${\utwi{\mathnormal\Gamma}}^{(1)}$.
Let $\widehat{ \Gamma}^{(1)} = (\widehat{ \gamma}^{(1)}_1,\ldots,\widehat{ \gamma}^{(1)}_q)$ be the
eigenmatrix of $\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1)}$ in the eigenvalue decomposition
\begin{eqnarray*}
\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1)} = \widehat{ \Gamma}^{(1)}\widehat{ D}} \def\bDtil{{\widetilde { D}}\big(\widehat{ \Gamma}^{(1)}\big)^\top,
\end{eqnarray*}
where $\widehat{ D}} \def\bDtil{{\widetilde { D}}=\hbox{\rm diag}(\hlam_1,\ldots,\hlam_q)$ gives the corresponding eigenvalues
with $\hlam_1\ge \cdots \ge \hlam_q$. The following theorem is an immediate
consequence of Theorem \ref{th-1new}.

\begin{theorem}\label{th-2}
Suppose conditions of Theorem \ref{th-1new} or \ref{th-1-poly} hold and $\eta^{**}_{T,p,q} < \Delta_G/2$
for a certain $G\subset \{1,\ldots,q\}$ with the $\eta^{**}_{T,p,q}$ in (\ref{th-1-2new}).
Let ${\utwi{P}}_{G} = {\utwi{\mathnormal\Gamma}}^{(1)}_{G}\Big({\utwi{\mathnormal\Gamma}}^{(1)}_{G}\Big)^\top$ be the orthogonal
projection to the corresponding eigenspace, and
$\widehat{ P}_G = \widehat{ \Gamma}^{(1)}_G\Big(\widehat{ \Gamma}^{(1)}_G\Big)^\top$
the sample version. Then, in the event $\Omega_T$ given in Theorem \ref{th-1new} with probability at least $1-\epsilon_T$,
\begin{eqnarray}\label{th-2-1}
\Big\| \widehat{ P}_G - {\utwi{P}}_G\Big\|_{op} \le \eta^{**}_{T,p,G},\
\hbox{with }\
\eta^{**}_{T,p,G} = \frac{\eta^{**}_{T,p,q}}{\Delta_G - \eta^{**}_{T,p,q}} < 1.
\end{eqnarray}
In particular, if $2\eta^{**}_{T,p,q} \le \lambda_j-\lambda_{j+1}$, then in the event $\Omega_T$
\begin{eqnarray*}
\Big\| \widehat{ P}_{\{1,\ldots,j\}} - {\utwi{P}}_{\{1,\ldots,j\}}\Big\|_{op}
\le \frac{\eta^{**}_{T,p,q}}{\lambda_j-\lambda_{j+1} - \eta^{**}_{T,p,q}} < 1.
\end{eqnarray*}
\end{theorem}

By the structure assumption in (\ref{b1}), there exists a group structure $G_c=\{ G_1,...,G_{n_c}\}$, such that the groups
\begin{eqnarray}\label{gamma-decomp}
 X_t\Big({\Sigma}}\def\hbSigma{{\widehat{\bSigma}}^{(1)}_1\Big)^{-1/2}{\utwi{\mathnormal\Gamma}}^{(1)}_{G_i},\ i = 1,\ldots,n_c,
\end{eqnarray}
are uncorrelated across all time lags. Let $G$ be one of the $G_i$ or a union of several $G_i$. It follows from Theorem \ref{th-2} that when the group $G$ is identifiable with a sufficiently large eigengap $\Delta_G$, it would be consistent with the eigenvalue decomposition of the estimator $\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1)}$, even in the presence of tied eigenvalues within $G$ or within $G^c$.

The decomposition (\ref{gamma-decomp}) is not unique. An alternative decomposition
\begin{eqnarray}\label{alt-decomp}
 X_t\Big({\Sigma}}\def\hbSigma{{\widehat{\bSigma}}^{(x)}_1\Big)^{-1/2}{\utwi{\mathnormal\Pi}}_{H_i}, \ i=1,\ldots, n_c,
\end{eqnarray}
based on an orthonormal matrix ${\utwi{\mathnormal\Pi}}=({\utwi{\mathnormal\gamma}}_1,\ldots,{\utwi{\mathnormal\gamma}}_q)$ and a partition $\{H_1,\ldots,H_{n_c}\}$
of $\{1,\ldots,q\}$ is equivalent to (\ref{gamma-decomp}) if and only if
\begin{eqnarray*}
{\utwi{\mathnormal\Pi}}_{H_i}{\utwi{\mathnormal\Pi}}_{H_i}^\top  = {\utwi{\mathnormal\Gamma}}^{(1)}_{G_i}\big({\utwi{\mathnormal\Gamma}}^{(1)}_{G_i}\big)^\top,\ i=1,\ldots, n_c.
\end{eqnarray*}
This includes rotation of the basis vectors within the group and also permutation of indices of the basis vectors. We say that a group $G_{i_0}$ in (\ref{gamma-decomp}) is {\it inseparable} at level $\rho^*>0$ if
for every equivalent decomposition (\ref{alt-decomp}),
there exists a connected graph $(H_{i_0},E_{i_0})$ with edges $E_{i_0}$ such that
\begin{eqnarray}\label{inseparable}
\rho_{k,\ell}= \max_{i,j,\tau}
\frac{\big|\big\langle {\utwi{\mathnormal\gamma}}_k, {\utwi{V}}^{(1)}_{\tau,i,j} {\utwi{\mathnormal\gamma}}_\ell\big\rangle\big|}
{\big\{\big\langle{\utwi{\mathnormal\gamma}}_k, {\utwi{V}}^{(1)}_{0,i,i} {\utwi{\mathnormal\gamma}}_k\big\rangle
         \big\langle {\utwi{\mathnormal\gamma}}_\ell,{\utwi{V}}^{(1)}_{0,j,j} {\utwi{\mathnormal\gamma}}_\ell\big\rangle \big\}^{1/2}}
         \ge \rho^*,\ \forall (k,\ell)\in E_{i_0}.
\end{eqnarray}
Under a proper identifiability condition on the eigengap $\Delta_{G_{i_0}}$ (e.g., $\Delta_{G_{i_0}}>2\eta_{T,p,q}^{**}$),
such inseparable group can be detected
with the test
\begin{eqnarray}\label{test}
\hrho^{(1)}_{k,\ell}
= \max_{i,j,\tau}
\frac{\big|\big\langle \widehat{ \gamma}^{(1)}_k, \widehat{ V}^{(1)}_{\tau,i,j} \widehat{ \gamma}^{(1)}_\ell\big\rangle\big|}
{\big\{\big\langle\widehat{ \gamma}^{(1)}_k, \widehat{ V}^{(1)}_{0,i,i} \widehat{ \gamma}^{(1)}_k\big\rangle
         \big\langle \widehat{ \gamma}^{(1)}_\ell,\widehat{ V}^{(1)}_{0,j,j} \widehat{ \gamma}^{(1)}_\ell\big\rangle \big\}^{1/2}}
        \ge \rho_0
\end{eqnarray}
with suitable $\rho_0 = \rho_{0,T} >0$,
where $\widehat{ \gamma}^{(1)}_j$ are the eigenvectors of $\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1)}$.
This is done by building a graph $(\{1,\ldots,q\},E)$ with edges
$E = \{(k,\ell): \hrho^{(1)}_{k,\ell} \ge\rho_0\}$.
This graph is split into connected subgraphs $(\Ghat_j, \Ehat_j)$, $j=1,\ldots,\ghat_1$,
which are not connected with each other.
We say that a group $G_{i_0}$ in (\ref{gamma-decomp}) is {\it detected} if
\begin{eqnarray}\label{detect}
\Ghat_j=G_{i_0},
\end{eqnarray}
for some $\ 1\le j\le\ghat_1$. The following theorem spells out the regularity conditions required for (\ref{detect})
in our analysis.

\begin{theorem}\label{th-3}
Suppose the conditions of Theorem \ref{th-2} hold with $G=G_{i_0}$
and that the minimum eigenvalue of ${\utwi{V}}^{(1)}_{0,j,j}$ is
not smaller than a constant $\phi_*>0$, for all $1\le j\le q$.
Suppose $(2\sqrt{2}c^*\eta^{**}_{T,p,G_{i_0}}/c_{*}) \vee \eta^{**}_{T,p,q} < \phi_*$ with the $c^*$ and $c_{*}$ defined in Assumption \ref{asmp:eigenvalue}, $\eta^{**}_{T,p,q}$ in (\ref{th-1-2new}) and $\eta^{**}_{T,p,G}$ in (\ref{th-2-1}). Define
\begin{align*}
&\rho^{+}_{G_{i_0}}=\frac{\rho^* - \big(\eta^{**}_{T,p,q} + 2\sqrt{2} c^*\eta^{**}_{T,p,G_{i_0}}/c_{*}\big)}
{(1 + 2\sqrt{2} c^*\eta^{**}_{T,p,G_{i_0}}/(c_{*}\phi_*))(1 + \eta^{**}_{T,p,q}/\phi_*)}, \\
&\rho^{-}_{G_{i_0}}= \frac{ \eta^{**}_{T,p,q} + 2\sqrt{2} c^*\eta^{**}_{T,p,G_{i_0}}/c_{*} }
{(1 - 2\sqrt{2} c^*\eta^{**}_{T,p,G_{i_0}}/(c_{*}\phi_*))(1 - \eta^{**}_{T,p,q}/\phi_*)}.
\end{align*}
Then, the following holds in the event $\Omega_T$,
\begin{eqnarray*}
\min_{k,\ell \in G_{i_0}} \hrho^{(1)}_{k,\ell} \ge \rho^{+}_{G_{i_0}},\quad
\max_{k\in G_{i_0},\ell\not\in G_{i_0}} \hrho^{(1)}_{k,\ell} \le \rho^{-}_{G_{i_0}}.
\end{eqnarray*}
In particular, when (\ref{inseparable}) holds and $\rho^{+}_{G_{i_0}} > \rho_0 > \rho^{-}_{G_{i_0}}$,
the group $G_{i_0}$ in (\ref{gamma-decomp}) is detected by the test (\ref{test})
in the sense of (\ref{detect}) with probability at least $1-\eps_T$.
\end{theorem}

\begin{remark}
The decomposition can be further refined with
\begin{eqnarray*}
\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1,f)} &=& \sum_{\tau=0}^{\tau_0} \sum_{i=1}^p\sum_{j=1}^p
\frac{f\big(\widehat{ V}^{(1)}_{\tau,i,j} \widehat{ V}^{(1)\top}_{\tau,i,j}\big) }{ p^2},
\cr \widehat{ \Gamma}^{(1,f)} &=& \hbox{the eigenmatrix of }\widehat{ W}} \def\bWtil{{\widetilde { W}}^{(1,f)},
\end{eqnarray*}
where $f({\utwi{V}})$ be a function of the $q\times q$ matrix inputs ${\utwi{V}}$, as defined in Remark \ref{rmk:trans}. Then all the above theorems still hold.
\end{remark}

\begin{remark}
Theorem \ref{th-3} establishes the consistency of the maximum cross correlation method in Section \ref{section:estimation}, whereas the property of its sample version is not presented in \cite{chang2018}. Our theoretical results require that the minimum eigen-gap satisfies $\Delta_G>2\eta^{**}_{T,p,q}$, where the estimation error bound $\eta^{**}_{T,p,q}$ is defined in \eqref{th-1-2new} and $\Delta_G$ is defined in \eqref{eq:deltag}.
\end{remark}

\section{Numerical properties}\label{sec4}

\subsection{Simulation} \label{section:simulation}

We illustrate the proposed decorrelation method with three simulated examples.
We always set $\tau_0=5$ in (\ref{est-2}), and $\tau_1 =15$ in (\ref{a13}).
The prewhitening, as stated at the end of Section \ref{section:estimation}, is
always applied in determining the orders of the columns of $\widehat {\utwi{A}}_*$ and $\widehat {\utwi{B}}_*$.
As the true transformation
matrices ${\utwi{A}}_*$ and ${\utwi{B}}_*$ defined in (\ref{an2}) cannot be computed
easily even for simulated examples, we use the following proxies instead
\[
{\utwi{A}}_*= \{ \widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}\}^{-1/2} {\utwi{A}} \{ \widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(u)} \}^{1/2}, \qquad
{\utwi{B}}_*= \{ \widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}\}^{-1/2} {\utwi{B}} \{ \widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(u)} \}^{1/2},
\]
where $\widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(x)}, \; \widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(x)}$ are given in (\ref{est-0}), and
\[
\widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_1^{(u)}= {1 \over T p} \sum_{t=1}^T {\utwi{U}}_t^\top {\utwi{B}}^\top {\utwi{B}} {\utwi{U}}_t, \qquad
\widehat {\Sigma}}\def\hbSigma{{\widehat{\bSigma}}_2^{(u)}= {1 \over T q} \sum_{t=1}^T {\utwi{U}}_t {\utwi{A}} {\utwi{A}}^\top {\utwi{U}}_t^\top.
\]
To abuse the notation, we still use ${\utwi{A}}_*, \; {\utwi{B}}_*$ to denote their proxies in
this section since the true ${\utwi{A}}_*, \; {\utwi{B}}_*$ never occur in our computation.
For each setting, we replicate the simulation 1000 times.

\begin{example}
We consider model (\ref{a1}) in which
all the component series of
${\utwi{U}}_t$ are independent and
ARMA(1, 2) of the form:
\begin{equation} \label{e1}
z_t = b z_{t-1} + \epsilon_t + a_1 \epsilon_{t-1} + a_2 \epsilon_{t-2},
\end{equation}
where $b$ is drawn randomly 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 ${\utwi{A}}$ and ${\utwi{B}}$
are drawn independently from $U(-1, \, 1)$. For this example, we assume
that we know $p_i=q_i=1$, $n_r=p$
and $n_c=q$ in (\ref{b1}). The focus is to investigate the impact of $p$ and $q$
on the estimation for ${\utwi{B}}$ and ${\utwi{A}}$.

Performing the eigenanalysis for $\widehat {\utwi{W}}^{(1)}$ and $\widehat {\utwi{W}}^{(2)}$ defined in (\ref{est-2}), we take the resulting two sets of orthonormal eigenvectors as the columns of
$\widehat {\utwi{A}}_*$ and $\widehat {\utwi{B}}_*$ respectively.
To measure the estimation error, we define
\begin{align*}
D(\widehat {\utwi{A}}_*, \, {\utwi{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{\utwi{A}}_{*}^{\top}
{\utwi{A}}_{*}$. Note that $D(\widehat {\utwi{A}}_*, \, {\utwi{A}}_{*})$ is always between 0 and 1,
and $D(\widehat{\utwi{A}}_*, \, {\utwi{A}}_{*})=0$ if $\widehat{\utwi{A}}_*$ is
a column permutation of ${\utwi{A}}_{*}$. $D(\widehat {\utwi{B}}_*, \, {\utwi{B}}_*)$ is defined in
the same manner.

\begin{table}[htbp]\caption{Example 1 -- Means and standard errors (SD) of $D(\widehat{\utwi{A}}_*, \,
{\utwi{A}}_{*})$ and $D(\widehat{\utwi{B}}_*, \, {\utwi{B}}_{*})$
in a simulation with 1000 replications.} \label{tab1}
\small
\begin{center}
\begin{tabular}{r|c|cc|cc||c|cc|cc|}
 &  & \multicolumn{2}{c|}{$D(\widehat{\utwi{A}}_*, \,
{\utwi{A}}_{*})$} &
\multicolumn{2}{c||}{$D(\widehat{\utwi{B}}_*, \, {\utwi{B}}_{*})$}
&  & \multicolumn{2}{c|}{$D(\widehat{\utwi{A}}_*, \,
{\utwi{A}}_{*})$} &
\multicolumn{2}{c|}{$D(\widehat{\utwi{B}}_*, \, {\utwi{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{\utwi{A}}_*, \, {\utwi{A}}_{*})$ and $D(\widehat{\utwi{B}}_*, \, {\utwi{B}}_{*})$ over the
1000 replications are reported in Table \ref{tab1}.
As expected, the estimation errors decrease as $T$ increases. Furthermore the error
in estimating ${\utwi{A}}_*$ increases as $q$ increases, and that in estimating ${\utwi{B}}_*$
increases as $p$ increases. Also
noticeable is the fact that the quality of the estimation for ${\utwi{A}}_*$ depends
on $q$ only and that for ${\utwi{B}}_*$ depends on $p$ only, as $D(\widehat{\utwi{A}}_*, \,
{\utwi{A}}_{*})$ is about the same for $(q,p)=(4, 4)$ and $(q, p)=(4,8)$,
and $D(\widehat{\utwi{B}}_*, \, {\utwi{B}}_{*})$ is about the same for $(q, p)=(4,8)$
and $(q, p)=(8,8)$. When $(q,p)=(32, 32)$, each ${\utwi{A}}$ and ${\utwi{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 ${\utwi{A}}_*$ is not accurate enough as $q=100$, while the
estimation for ${\utwi{B}}_*$ remains as good as that for $(q,p)=(16,16)$, and is clearly
better than that for $(q,p)=(32,32)$.
\end{example}

\begin{example}
To examine the performance of the proposed method for
identifying uncorrelated blocks, we consider now a only column-segmentation model
\[
 X_t =  {\utwi{U}}_t {\utwi{A}}^\top,
\]
where $ X_t, {\utwi{U}}_t$ are $p\times q$ matrix time series. Note that adding the row
transformation matrix ${\utwi{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 ${\utwi{A}}_*$ and that for ${\utwi{B}}_*$ are performed separately. See Section
\ref{section:estimation}.

The transformation matrix
${\utwi{A}}$ in the above model is generated in the same way as in Example 1 above. All the element time series
of ${\utwi{U}}_t$, except those in columns 2, 3 and 5, are simulated
independently from ARMA(1,2) model (\ref{e1}).
Denoted by ${\utwi{U}}_{t, j}$ the $j$-th column of ${\utwi{U}}_t$,
we let
\[
{\utwi{U}}_{t,2}={\utwi{U}}_{t+1,1}, \quad {\utwi{U}}_{t,3}={\utwi{U}}_{t+2,1}, \quad
 {\utwi{U}}_{t,5}={\utwi{U}}_{t+1,4}.
\]
Hence the first three columns of ${\utwi{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 {\utwi{W}}^{(1)}$
defined as (\ref{est-2}). We then apply the algorithm in Section \ref{section:estimation}
to arrange the orders of the $q$ resulting orthonormal eigenvectors to form
estimator $\widehat {\utwi{A}}_*$, for which the number of the
correlated pairs of columns is determined by (\ref{a14}).
Table \ref{tab2} 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{tab2} 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
${\utwi{A}}_{*} = ({\utwi{A}}_{1}, \cdots, {\utwi{A}}_{n_c}) $ and
$\widehat {\utwi{A}}_* = (\widehat {\utwi{A}}_{1}, \cdots, \widehat {\utwi{A}}_{n_c})$. We measure
the estimation error by
\begin{align*}
D_1 (\widehat {\utwi{A}}_*, {\utwi{A}}_{*}) = \frac{1}{n_c} \sum_{j=1}^{n_c}
\big\{ 1 - \frac{1}{{\rm rank}({\utwi{A}}_{j})} {\rm trace}\big({\utwi{A}}_{j} {\utwi{A}}_{j}^{\top} \widehat{\utwi{A}}_{j} \widehat{\utwi{A}}_{j}^{\top}\big) \big\}.
\end{align*}
Then $\{{\rm rank}({\utwi{A}}_{j})\}^{-1} {\rm trace}\big({\utwi{A}}_{j}
{\utwi{A}}_{j}^{\top} \widehat{\utwi{A}}_{j} \widehat{\utwi{A}}_{j}^{\top}\big)$ is always between 0 and
1. Furthermore it is equal to 1 if ${\cal M}({\utwi{A}}_{j}) = {\cal M}(\widehat {\utwi{A}}_{j})$,
and 0 if the two spaces are perpendicular with each other; see, for
example, \cite{stewart1990} and \cite{pan2008}. The boxplots of $D_1 (\widehat
{\utwi{A}}_{*}, {\utwi{A}}_{*}) $ presented
in Figure 1 indicate that the estimation errors decrease when $T$ increases, and
the errors with large $q$ are greater than those with small $q$.

\begin{table}[htbp]
\caption{Example 2 -- Relative frequencies correct and incorrect
segmentation in a simulation with 1000 replications.
} \label{tab2}
\small
\begin{center}
\begin{tabular}{ll|cccc}
$(q,p)$&  & $T=100$ & $T=500$ & $T=1000$& $T=5000$\\
\hline
\hline
$(6,3)$ & Correct & 0.369& 0.659& 0.721 & 0.837 \\
        & Merging
        & 0.274& 0.237 & 0.220& 0.142 \\
        & Splitting
        & 0.299& 0.056 & 0.031& 0.008\\
        & Other & 0.058 & 0.048 & 0.028 & 0.013 \\

\hline
$(6,6)$ & Correct & 0.371& 0.639& 0.706 & 0.815   \\
        & Merging & 0.207 & 0.218 & 0.219& 0.156 \\
        & Splitting & 0.208 & 0.071 & 0.030& 0.011 \\
        & Other & 0.214 & 0.072 & 0.045 & 0.018 \\

\hline
$(10,6)$ & Correct & 0.103 & 0.323 & 0.369 & 0.542 \\
        & Merging & 0.087 & 0.148& 0.227& 0.248 \\
        & Splitting & 0.270 & 0.141 & 0.086&0.014 \\
        & Other & 0.540 & 0.388 & 0.318 & 0.196 \\

\hline
$(10,10)$ & Correct & 0.116 & 0.299 & 0.363& 0.495  \\
        & Merging & 0.091 & 0.195& 0.218&    0.269  \\
        & Splitting & 0.249 & 0.116&0.067&   0.019  \\
        & Other & 0.544 & 0.390 & 0.352 & 0.217 \\

\end{tabular}
\end{center}
\end{table}

\begin{figure}[ht]
\begin{center}\includegraphics[scale=0.7]{fig1}\end{center}
\caption{Example 2 -- boxplots of $D_1(\widehat {\utwi{A}}_*, {\utwi{A}}_{*})$ in a simulation
with 1000 replications.} \label{fig1}
\end{figure}
\end{example}

\begin{example}
In this example we examine the gain in post-sample forecasting due to the decorrelation
transformation.
Now in model (\ref{a1}) all the elements of ${\utwi{A}}$ and ${\utwi{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{table:example3} shows the configuration of ${\utwi{U}}_t$ for $(p,q)=(8,8)$.

\begin{table}\label{table:example3}
\small
\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 ${\utwi{U}}_t$ for the case of $(p,q)=(8,8)$}
\end{table}

 Each univariate block of ${\utwi{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 \cite{chen2020autoregressive}, i.e.
\begin{align}\label{model:ex3}
{\utwi{U}}_{t,i,j}=\lambda_{i,j}{\utwi{\mathnormal\Phi}}_{1,i,j} {\utwi{U}}_{t-1,i,j} {\utwi{\mathnormal\Phi}}_{2,i,j}^\top+{\utwi{E}}_{t,i,j}, \quad
i,j =1, 2,
\end{align}
where ${\utwi{\mathnormal\Phi}}_{1,i,j}$ and
${\utwi{\mathnormal\Phi}}_{2,i,j}$ are the coefficient matrices with unit spectral norm,
all elements of ${\utwi{E}}_{t,i,j}$ are independent and $N(0,1)$. Matrices ${\utwi{\mathnormal\Phi}}_{1,i,j}$ and
${\utwi{\mathnormal\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
scaler coefficient $\lambda_{i,j}$ controls
the level of auto-correlation and its values for the four blocks
are marked in Table~\ref{table:example3}.

\end{example}

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+i$ observations for identifying the uncorrelated blocks
and for computing  $\widehat {\utwi{A}}_*, \widehat {\utwi{B}}_*$ and
$\widehat {\utwi{U}}_t^*$ defined in (\ref{est-3}). Depending on its shape/size,
we fit each identified block in $\widehat {\utwi{U}}_t^*$ with an MAR(1), VAR(1) or AR(1) model.
The fitted models are used to predict ${\utwi{U}}_{T+i+1}^*$.
The predicted value for $ X_{T+i+1}$ is then obtained by
\begin{align}\label{model2:ex3}
\widehat X_{T+i+1}= \big(\widehat{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}^{(x)}_2\big)^{1/2}
\widehat{\utwi{B}}_*\widehat{\utwi{U}}_{T+i+1}^* \widehat{\utwi{A}}_*^\top \big(\widehat{\Sigma}}\def\hbSigma{{\widehat{\bSigma}}^{(x)}_1\big)^{1/2}, \quad i=0, 1,
\cdots, M-1.
\end{align}
We then compute the mean squared error:
\begin{equation} \label{rmse}
{\rm MSE}=\frac{1}{pqM}\sum_{s=T+1}^{T+M}\sum_{i=1}^p\sum_{j=1}^q
(\widehat X_{s,i,j}-\mu_{s,i,j})^2,
\end{equation}
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 $ X_{s,i,j}$.
We report overall the mean and the standard deviation of MSE over 1000 replications.
In addition we also report
the mean and the standard deviation over the replications with the completely
correct segmentation of the uncorrelated blocks in ${\utwi{U}}_t^*$, over those with
nearly correct segmentation (i.e. either $\widehat n_r = n_r -1$ or
$\widehat n_c = n_c -1$), and over all the other cases.

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, and (ii)
MAR(1) model of \cite{chen2020autoregressive}, i.e. MAR(1) is fitted directly
to $ X_t$.
In addition, we also include two oracle models:
(O1) the true segmentation blocks are used to model $\widehat{\utwi{U}}_t^*$; and
(O2) True $\Sigma^{(x)}_1,\Sigma^{(x)}_2$ and true ${\utwi{A}}_*,{\utwi{B}}_*$ are used
and the true segmentation blocks are used to model $\widehat{\utwi{U}}_t^*$.

The mean and standard deviations (in bracket) of the MSEs are listed in Table \ref{tab3}.
It is clear that the forecasting based on the proposed decorrelation transformation
is more accurate than those based on MAR(1) and VAR(1) 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).

Table \ref{tab3} 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.

\begin{table}[htbp]
\caption{Example 3 -- One-step ahead post-sample forecasting: means and
standard deviations (in bracket) of RMSEs and the proportions of correct
segmentation, near correct segmentation and others in a simulation with
1000 replications.}\label{tab3}
\small{
\begin{center}
\begin{tabular}{l|cccccccc}
 & \multicolumn{8}{c}{RMSEs and the proportions of segmentations} \\
\hline
 & \multicolumn{2}{c}{$T=500$} & \multicolumn{2}{c}{$T=1000$} & \multicolumn{2}{c}{$T=2000$} & \multicolumn{2}{c}{$T=5000$} \\
 \hline \hline
$(p,q)$= $(6,6)$ \\
\hline
 MAR(1) & 0.571(0.199) & & 0.562(0.196) & & 0.550(0.194) & & 0.556(0.204) & \\
 VAR(1) & 0.421(0.083) & & 0.293(0.082) & & 0.208(0.059) & & 0.130(0.039) & \\
 segmentation & 0.361(0.165) & & 0.253(0.132) & & 0.192(0.104) & & 0.132(0.076)  \\ \hline
 model O1  & 0.337(0.136) & & 0.242(0.096) & & 0.185(0.096) & & 0.132(0.088)    \\
 model O2  & 0.106(0.032) & & 0.073(0.021) & & 0.052(0.014) & & 0.033(0.010) \\
 \hline \hline
$(p,q)$= $(6,8)$ \\
\hline
 MAR(1) & 0.674(0.196) & & 0.664(0.192) & & 0.672(0.201) & & 0.673(0.202) & \\
 VAR(1) & 0.540(0.133) & & 0.385(0.095) & & 0.277(0.069) & & 0.186(0.043) & \\
 segmentation & 0.465(0.175) & & 0.330(0.152) & & 0.250(0.122) & & 0.174(0.107)  \\ \hline
 model O1  & 0.457(0.153) & & 0.322(0.137) & & 0.264(0.137) & & 0.193(0.130)    \\
 model O2  & 0.130(0.034) & & 0.091(0.024) & & 0.063(0.016) & & 0.041(0.011) \\
\hline \hline
$(p,q)$= $(8,8)$ \\
\hline
 MAR(1) & 0.685(0.178) & & 0.674(0.182) & & 0.681(0.194) & & 0.668(0.186) & \\
 VAR(1) & 0.721(0.137) & & 0.517(0.113) & & 0.388(0.077) & & 0.245(0.050) & \\
 segmentation & 0.566(0.167) & & 0.441(0.143) & & 0.342(0.121) & & 0.234(0.106)  \\ \hline
 model O1  & 0.513(0.144) & & 0.415(0.136) & & 0.330(0.120) & & 0.248(0.118)    \\
 model O2  & 0.152(0.034) & & 0.108(0.024) & & 0.075(0.017) & & 0.047(0.010) \\
 \hline \hline
$(p,q)$= $(10,10)$ \\
\hline
 MAR(1) & 0.857(0.176) & & 0.846(0.178) & & 0.840(0.172) & & 0.838(0.162) & \\
 VAR(1) & 1.280(0.171) & & 0.892(0.135) & & 0.622(0.095) & & 0.386(0.062) & \\
 segmentation & 0.747(0.184) & & 0.619(0.155) & & 0.514(0.136) & & 0.369(0.126)  \\ \hline
 model O1  & 0.700(0.150) & & 0.594(0.149) & & 0.510(0.133) & & 0.385(0.129)    \\
 model O2  & 0.192(0.033) & & 0.136(0.024) & & 0.095(0.016) & & 0.061(0.010) \\
\end{tabular}
\end{center}
}
\end{table}

\subsection{Real data analysis} \label{section:data}

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 used
in this excise is January 1, 2015 -- December 31, 2016 for a total
of $T=731$ days.
Figure
\ref{figure:data} displays the time series plot of the $17\times 6=102$ series after removing
the seasonal means.
The readings of different pollutants at different locations are crossly correlated. Figure \ref{figure:row} shows the cross correlogram of five randomly selected stations for PM$_{10}$, and Figure \ref{figure:col} displays the cross correlogram of all six pollutants for a station in Beijing.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.55,page=1]{data_demean}
\caption{Time series plot of 17 by 6 air pollution data}
\label{figure:data}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5,page=1]{ccf_row2}
\caption{Cross correlogram of logarithmic daily PM$_{10}$ readings at five randomly selected monitoring stations in Beijing, Tianjin and Hebei}
\label{figure:row}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5,page=1]{ccf_col2}
\caption{Cross correlogram of logarithmic daily readings of all six pollutants for a station in Beijing}
\label{figure:col}
\end{figure}

We first deduct seasonal mean of the $17\times 6$ matrix time series.
Setting $\tau_0=5$ in (\ref{est-2}), 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{a13}), 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 all the others
of size 1.

To check the post-sample forecasting performance, we calculate the rolling
one-step, and two-step ahead forecasts for each of the daily readings in the last
three months in 2016 (in total 92 days).
The methods included in the comparison are 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/shape,
MAR(1) for the original $17\times 6$ matrix series, VAR(1) for the vectorized
original series, and univariate AR(1) for each of the $17\times 6$ original series.
For the segmentation approach, we fix the transformation obtained using
data upto September 30, 2016.
In addition to the identified segmentation with 15$\times$4 blocks,
we also compute the forecasts based on 14$\times$3 blocks, and 16$\times$5 blocks.
The former is obtained by merger two single row blocks together and two single column
blocks together. The latter is obtained by splitting a block with two rows into two
and splitting the block with 3 columns into two blocks.
This will reveal the impact of `imperfect' segmentation on forecasting performance.

The means and the standard deviations of MSE (see (\ref{rmse}))
of one-step and two-step ahead post-sample forecasts for the readings
in the last 92 days are listed in Table \ref{tab:pollution}.
The prediction based on the decorrelation transformation (even using imperfect segmentations) is clearly  more accurate than those without transformation.
Using the imperfect segmentation deteriorates the performance only slightly,
indicating that a partial (instead of total) decorrelation still leads to
significant gain in prediction. 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 dimension-reduction techniques such as the proposed decorrelation
transformation.

To further evaluate the forecasting stability, we calculate the rolling MSPE for one-step ahead post-sample forecasts based on, respectively, the segmentation, MAR(1), VAR(1) and univariate AR models. We plot the weekly averages of the rolling MSPE in Figure \ref{figure:weakly_mean}. It shows that the forecast based on the proposed segmentation almost always outperforms the other three methods.

\begin{table}[htbp]
\caption{One-step and two-step ahead post-sample forecasting for the air pollution data: means and standard deviations (in bracket) of MSPE} \label{tab:pollution}
\ \\
\centering
\begin{tabular}{lcccccccc}\hline
Method & One-step forecast & Two-step forecast \\\hline
MAR(1) & 0.351 (0.208) & 0.516 (0.303) \\
VAR(1) & 0.400 (0.272) & 0.717 (0.406) \\
Univariate AR & 0.346 (0.205) & 0.501 (0.286) \\ \hline
Segmentation with 15$\times$4 blocks & 0.315 (0.182) & 0.435 (0.246) \\
Segmentation with 14$\times$3 blocks & 0.319 (0.185) & 0.446 (0.256) \\
Segmentation with 16$\times$5 blocks & 0.318 (0.183) & 0.443 (0.247) \\
\hline
\end{tabular}
\end{table}

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.6,page=1]{weakly_mean}
\caption{Weakly averaged one-step forecast errors of four selected methods}
\label{figure:weakly_mean}
\end{figure}

\bibliographystyle{apalike}

redpajama_set_name: RedPajamaArXiv


As a leading global manufacturer of crushing and milling equipment, we offer advanced, Posts Related to how to make gold concentrator.
Gold Crushing Machine, Tags Gold Concentrator Machine Used Gold Concentrators Gold Mining Machines. Compare . 1-3 ton/hour stone gold separator machine.
Gold concentrators equipment, Gravity Concentrators. GOLD, CONCENTATORS. The plant uses a plate feed hopper that feeds the ore to a jaw crusher.


text:
import { ExcelXMLTemplate } from '@ag-grid-community/core';
declare const row: ExcelXMLTemplate;
export default row;

redpajama_set_name: RedPajamaGithub


text:
Q: SQL - count without group by? I need to use two ids for a join I thought I could count a column and add it as a column as I can with a sum but I get an error about having to group by / having.  An example of what I want...
Initial table...
 Global ID    Local ID   Name       Role
 100          1          Andy       Manager
 100          2          Andy       Manager
 100          1          John       Co-Manager
 200          1          Andy       Co-Manager
 200          2          John       Manager
 200          2          Mike       Manager

I then want to add a column that counts the number of Manager in each group / local pairing...
 Global ID    Local ID   Name       Role         Manager Count
 100          1          Andy       Manager      1
 100          2          Andy       Manager      1
 100          1          John       Co-Manager   0
 200          1          Andy       Co-Manager   0
 200          2          John       Manager      2
 200          2          Mike       Manager      2

I tried just joining the two tables on the global / local ID and then adding a column for count of the Role.  The problem is that I get an error about grouping / having by but I don't want to group anything.  I just want to add the column and still have the same number of rows.  Any way around this?
FYI - for the last two rows, the last column has 2 because John and Mike are on the same group / local ID

A: It looks like your definition of a group is a unique combination of [global_id, local_id]. If that is the case you do want to group by those two values, and do a count, where the role is Manager. But because you want other columns from the original table, you must do that count within an inline view, and then join back to the original table, like so:
select t.*, v.mgr_in_grp
  from tbl t
  left join (select global_id, local_id, count(*) as mgr_in_grp
               from tbl
              where role = 'Manager'
              group by global_id, local_id) v
    on t.global_id = v.global_id
   and t.local_id = v.local_id

Fiddle: http://sqlfiddle.com/#!2/fb3ace/2/0
Notice that there is a difference on row 3, as compared to your expected output.
John, Co-manager at global_id and local_id 100 and 1, respectively, belongs to a pair that has a manager in that combination (Andy, your first row). So the count appears should be 1.

redpajama_set_name: RedPajamaStackExchange


Had the chance to shoot the fellas from Interstellar Cave Dweller at a recent performance at Iowa Brewing Company in Cedar Rapids, Iowa. They always put on a great show! Go check them out if you have a chance!


Often times hunters wonder what a particular buck would score if he were older. We have developed an equation that will predict what a buck will score if the buck followed normal age-related growth patterns. This equation was developed from deer harvested in Mississippi and Texas, so most likely the equation will apply very well to your area. Simply enter the Boone & Crockett score in the box below, select the age of the buck, and the equation will do the rest.
For example, say I was lucky enough to harvest a buck that scored 125 inches and a wildlife biologist examined the jawbone and determined the buck to be 3.5 years old. What would this buck have scored at 4.5 or 5.5 years of age? Enter 125 in the B&C score box, click on the 3.5 button, and click the calculate button. The predicted Boone & Crockett scores of 145 inches at 4.5 years of age, and 157 inches at 5.5 years of age will be shown in the boxes below.


text:
And research published in November 2017 in the journal NeuroImage showed that aerobic exercise may be helpful in improving memory function and maintaining brain health as we age. The study, led by researchers at Australia's National Institute of Complementary Medicine at Western Sydney University, looked at brain scans of 737 people ranging in ages from 26 to 76. The group included a mix of healthy adults, people with Alzheimer's and other cognitive impairments, and people with a clinical diagnosis of mental illness, including depression and schizophrenia. The researchers found that exercises, like riding a stationary bike, running on a treadmill, or walking, slowed down the deterioration of brain size and slowed the effect of age on brain health.
Studies have consistently shown that physical activity can help treat depression, and on the flipside, that low activity levels are a big risk factor for it. The antidepressant effect of exercise seems to be moderated in part through serotonin, the brain chemical that's targeted with some antidepressants, and in part through bone-derived neurotrophic factor (BDNF). And this goes back to the generation of new cells mentioned earlier—exercise, though various mechanisms, seems to make the brain more plastic and more capable of growing new cells.
Exercise has long been correlated with a longer life, but it's only recently started to become clear why this might be. Studies, like a new one in the journal Preventive Medicine which found that exercise is linked to longer caps at the ends of chromosomes, have helped flesh this out a bit more. These caps, called telomeres, naturally shorten as we age, with each cell division. People who live a long time have telomeres that are in better shape than those who don't—but there's a lot we can do to affect the rate at which they shorten over the years. The team behind the new study looked at data from CDC's National Health and Nutrition Examination Survey, and found that for people who exercised regularly, their telomeres were 140 base pairs longer on average than sedentary people's. Which correlates to being years "younger" than their sedentary peers.

redpajama_set_name: RedPajamaCommonCrawl


"Fluted" interior textured mirrored glass wall, viewed from a distance. Note: Installed glass appearance can vary significantly based on actual lighting, viewing angle and distance. Nothing comes close to the representative accuracy of a full-size mockup.
"Fluted" interior textured mirrored glass wall, viewed from a greater distance. Note: Installed glass appearance can vary significantly based on actual lighting, viewing angle and distance. Nothing comes close to the representative accuracy of a full-size mockup.


text:
Skanska/
Skanska to renovate and expand Brookline High School in Massachusetts, USA, for USD 128 M, about SEK 1.1 billion
Tue, Jan 08, 2019 07:30 CET
Skanska has signed a contract with Brookline High School in Massachusetts, USA. The contract is worth USD 128 M, about SEK 1.1 billion, which will be included in the order bookings for Skanska USA for the fourth quarter of 2018.
Skanska will construct a new, approximately 11,000-square-meter academic building that will include classrooms; event space with dining and food service; a medical clinic; a library; and spaces for general administration and collaboration. In addition, Skanska will add 6,500 square meters to the school's existing STEM Wing that will accommodate new biology and chemistry classrooms and labs; collaboration and maker spaces; a culinary arts kitchen and student restaurant with seating; and a 300-square-meter lobby and Community Wellness Room in the existing Tappan Athletic Complex.
Skanska will also perform substantial upgrades to all the school's existing buildings, landscaping and streetscapes.
Construction began in November 2018, and completion is scheduled for the third quarter of 2021. The project architect is William Rawn Associates.
Skanska is one of the leading construction and development companies in the USA, specializing in building construction, civil infrastructure, public-private partnerships and commercial property developments in select US markets. Skanska USA had sales of SEK 59 billion in 2017 and has approximately 9,500 employees.
Katie Koch, Communications, Skanska USA, tel +1 (917) 438 33 67
Patric Elmén, acting Head of Media Relations, tel +46 (0) 449 31 34
Direct line for media, tel +46 (0)10 448 88 99
This and previous releases can also be found at www.skanska.com
Skanska is one of the world's leading construction and project development companies, focused on selected home markets in the Nordic region, Europe and USA. Supported by global trends in urbanization and demographics, and by being at the forefront in sustainability, Skanska offers competitive solutions for both simple and the most complex assignments, helping to build a sustainable future for customers and communities. The Group has about 40,000 employees. Skanska's sales in 2017 totaled SEK 161 billion.
Entré Lindhagen, Warfvinges väg 25, 112 74 Stockholm, Sweden
http://www.skanska.com
Jacob Birkeland
Head of Media Relations
jacob.birkeland@skanska.se
20190108 Brookline High School
Privacy statement | Cookies | About | Content Policy | Contact
©Cision Privacy statement Cookies About Content Policy Contact


text:
Anime Expo Hosts Fuji TV Producer Yōko Matsuzaki
posted on 2019-06-26 19:00 EDT by Rafael Antonio Pineda
+Ultra panel also announces new project with BONES , Wit Studio staff
Anime Expo announced on Tuesday that it will host Fuji TV general manager Yōko Matsuzaki . Matsuzaki will participate in a Fuji TV panel about its +Ultra programming block on July 6 at 4:30 p.m. Manga artist Yoshiyuki Sadamoto , director Hiro Kaburagi , scriptwriter Ryota Kosawa , Wit Studio president and CEO George Wada , and BONES studio president Masahiko Minami will join Matsuzaki to announce a new project and screen the first episode of Bones' Carole & Tuesday anime.
Matsuzaki has had a long career producing both live-action and anime series for Fuji TV . She was involved in the launch of Fuji TV 's +Ultra and Noitamina programming blocks. Anime Matsuzaki has produced include the long-running Sazae-san and One Piece series, as well as Psycho-Pass , Kabaneri of the Iron Fortress , After the Rain , The Promised Neverland , and Sarazanmai .
This year's Anime Expo will take place from July 4 to 7 at the Los Angeles Convention Center and the surrounding area.
Source: Anime Expo
discuss this in the forum (2 posts) |