# Factor Modelling for High-Dimensional Functional Time Series

Shaojun Guo  Institute of Statistics and Big Data, Renmin University of China,
P.R. China  Xinghao Qiao  Department of Statistics, London School of
Economics, U.K.  Qingsong Wang  Institute of Statistics and Big Data, Renmin
University of China, P.R. China

###### Abstract

Many economic and scientific problems involve the analysis of high-dimensional
functional time series, where the number of functional variables ($p$)
diverges as the number of serially dependent observations ($n$) increases. In
this paper, we present a novel functional factor model for high-dimensional
functional time series that maintains and makes use of the functional and
dynamic structure to achieve great dimension reduction and find the latent
factor structure. To estimate the number of functional factors and the factor
loadings, we propose a fully functional estimation procedure based on an
eigenanalysis for a nonnegative definite matrix. Our proposal involves a
weight matrix to improve the estimation efficiency and tackle the issue of
heterogeneity, the rationality of which is illustrated by formulating the
estimation from a novel regression perspective. Asymptotic properties of the
proposed method are studied when $p$ diverges at some polynomial rate as $n$
increases. To provide a parsimonious model and enhance interpretability for
near-zero factor loadings, we impose sparsity assumptions on the factor
loading space and then develop a regularized estimation procedure with
theoretical guarantees when $p$ grows exponentially fast relative to $n.$
Finally, we demonstrate that our proposed estimators significantly outperform
the competing methods through both simulations and applications to a U.K.
temperature dataset and a Japanese mortality dataset.

Keywords:  Autocovariance function; Eigenanalysis; Functional thresholding;
Sparse principal component analysis; Strength of factors, Weight matrix.

##  1 Introduction

Functional time series, which refers to a sequential collection of curves
observed over time exhibiting serial dependence, has recently received a great
deal of attention. Despite progress being made in this field, existing
literature has focused on the estimation based on a single or fixed number of
functional time series, see, e.g., Bosq (2000); Hörmann and Kokoszka (2010);
Bathia et al. (2010); Panaretos and Tavakoli (2013); Hörmann et al. (2015);
Aue et al. (2015); Li et al. (2020); Chen, Guo and Qiao (2020) and among many
others.

With the rapid development in technology, datasets containing a large
collection of functional time series are becoming increasingly available in
various applications. Examples include intraday return density trajectories
Bathia et al. (2010) for hundreds of stocks, annual temperature curves
collected at a number of stations, daily energy consumption curves for
thousands of households, age-specific mortality rates in different prefectures
Gao et al. (2019), to list a few. Such data, which are referred to as high-
dimensional functional time series, take the form of ${\bf
Y}_{t}(\cdot)=\\{Y_{t1}(\cdot),\dots,Y_{tp}(\cdot)\\}^{{\mathrm{\scriptscriptstyle
T}}}$ ($t=1,\dots,n$) defined on a compact interval ${\cal U},$ where the
number of random functions ($p$) is comparable to, or even larger than the
number of serially dependent observations ($n$). Under such scenario, not only
$p$ is large, but each $Y_{tj}(\cdot)$ is a functional object with serial
dependence across observations, posing extra challenges. When jointly
modelling the entire curve dynamics, it is of great interest and importance to
explore the latent dynamic structure from a dimension-reduction viewpoint,
while preserving the functional and time series structure. This motivates us
to develop a functional factor model based on ${\bf Y}_{1}(\cdot),\dots,{\bf
Y}_{n}(\cdot).$

In this paper, we deal with the factor modelling for high-dimensional
functional time series, in which each functional datum ${\bf Y}_{t}(\cdot)$
arises as the sum of two unobservable components, one common ${\bf
Z}_{t}(\cdot)$ and one idiosyncratic $\boldsymbol{\varepsilon}_{t}(\cdot).$
The co-movement of $p$ functional time series is assumed to be driven by a
small number of functional factors and their inherited dynamic structures in
${\bf Z}_{t}(\cdot)$. Specifically, our proposed factor model admits the
representation:

${\bf Y}_{t}(\cdot)={\bf
Z}_{t}(\cdot)+\boldsymbol{\varepsilon}_{t}(\cdot)={\bf A}{\bf
X}_{t}(\cdot)+\boldsymbol{\varepsilon}_{t}(\cdot),~{}~{}t=1,\dots,n,$ (1)

where ${\bf
X}_{t}(\cdot)=\\{X_{t1}(\cdot),\dots,X_{tr}(\cdot)\\}^{{\mathrm{\scriptscriptstyle
T}}}$ is a set of latent functional factor time series with (unknown) number
of factors $r<p,$ ${\bf A}=(A_{ij})_{p\times r}$ is the factor loading matrix
and the idiosyncratic component $\boldsymbol{\varepsilon}_{t}(\cdot)$ is a
sequence of white noise process.

For multivariate scalar time series $\\{{\bf Y}_{t}\\},$ the factor model (1)
degenerates to the classical dynamic factor model ${\bf Y}_{t}={\bf A}{\bf
X}_{t}+\boldsymbol{\varepsilon}_{t},$ where two different types of model
assumptions are considered in econometrics and statistics literature. One type
of models assumes that each common factor has non-trivial contribution to
$\\{{\bf Y}_{t}\\}$ while the idiosyncratic noise
$\\{\boldsymbol{\varepsilon}_{t}\\}$ is allowed to have weak cross-
correlations and serial correlations Chamberlain and Rothschild (1983); Forni
et al. (2000); Bai and Ng (2002); Stock and Watson (2012); Fan et al. (2013,
2016). However, the rigorous definition of the common and idiosyncratic
components in those factor models can only be established asymptotically when
the dimension of ${\bf Y}_{t}$ tends to infinity. Another type of models
assumes that the common factors accommodate all dynamics of $\\{{\bf
Y}_{t}\\},$ thus making the idiosyncratic component white noise with no serial
correlations but allowing substantial contemporary cross-correlation among it;
see Pena and Box (1987); Pan and Yao (2008); Lam et al. (2011); Lam and Yao
(2012) and extensions to matrix-valued time series Wang et al. (2019); Chen,
Tsay and Chen (2020) and tensor-valued time series Chen et al. (2021). For
functional time series, the factor modelling remains less touched in the
literature. Hays et al. (2012) and Kokoszka et al. (2015) considered
functional dynamic factor models for univariate functional time series. Gao et
al. (2019) adopted factor modelling techniques to predict high-dimensional
functional time series by fitting a factor model to estimated scores obtained
via dynamic functional principal component analysis Hörmann et al. (2015).
Following the framework in Bai and Ng (2002), a recent manuscript of Tavakoli
et al. (2021) develops a high-dimensional functional factor model with scalar
factors and functional factor loadings, which substantially differs from our
model setup in (1).

Our paper is along the line of Lam and Yao (2012) in the model development.
The adopted white noise assumption is not the most ideal for the purpose to
extract those asymptotically identifiable factors, however, since the white
noise exhibits no serial correlations, the decomposition in (1) is unique,
thus largely simplifying the tasks of model identification and inference.
Moreover, unlike Bai and Ng (2002), we do not place any assumption on the
covariance structure of the idiosyncratic noise. We propose to estimate both
the number of factors and the factor loadings in terms of an eigenanalysis for
a nonnegative definite matrix, which is the double integral and sum of
weighted quadratic forms in autocovariance functions of observed curves from
different time lags. In practical applications, it is common that many
estimated factor loadings are close to zero especially when $p$ is large in
relation to $n,$ This inspires us to add certain sparsity constraints on the
factor loading space defined in terms of the factor loading matrix ${\bf A}.$
Since functional sparsity patterns in autocovariance functions of ${\bf
Y}_{t}(\cdot)$ inherits from the imposed sparsity structures in ${\bf A},$ we
first propose Hilbert–Schmidt-norm based functional thresholding of entries of
the sample auotocovariance functions, the consistencies of which are
guaranteed under $\log p/n\rightarrow 0.$ Furthermore, to recover sparse
factor loading space and produce a parsimonious model with interpretable
factors, we perform sparse principal component analysis Vu and Lei (2013)
rather than standard eigenanalysis for our formed nonnegative definite matrix.

This paper makes useful contributions at multiple fronts. On the method side,
our proposal has four main advantages. First, it involves a weight matrix to
improve the estimation efficiency and tackle the issue of heterogeneity. We
illustrate the rationality of its presence by formulating the estimation
procedure based on a novel regression viewpoint. Our empirical results also
demonstrate the uniform superiority of our method compared to the unweighted
competitors. Furthermore, our strategy of incorporating the weight matrix is
general and can be applied broadly to a range of dimension reduction problems,
e.g., Bathia et al. (2010), Lam and Yao (2012) and Chang et al. (2018).
Second, our estimation relies on the sample autocovariance functions of ${\bf
Y}_{t}(\cdot).$ This not only gets rid of the impact from the noise
$\boldsymbol{\varepsilon}_{t}(\cdot),$ but also makes the good use of the
serial correlation information, which is most relevant in the time series
modelling. Third, the double integral takes advantage of the functional nature
of the data by gathering the autocovariance information at each $(u,v)\in{\cal
U}^{2}$ and then integrating over ${\cal U}^{2}.$ By contrast, the estimation
in Tavakoli et al. (2021) is based on the inner product between pairs of
observed curves and hence ignores the covariance information. Moreover, our
fully functional procedure avoids the information loss incurred by the
dimension reduction step in Gao et al. (2019). Fourth, aided by the enforced
sparsity, we enhance the model interpretability and enlarge the dimension of
the feature space that our model can handle. Our regularized estimation
procedure adopts techniques of functional thresholding and sparse principal
component analysis, each of which leads to the estimation consistency and
sparsity recovery at its own step.

On the theory side, we investigate convergence properties of relevant
estimated terms under high-dimensional scaling, where $p$ grows polynomially
fast relative to $n.$ With additional sparsity constraints and tail
assumptions, we further provide convergence analysis of the regularized
estimation in an ultra-high-dimensional regime, where $p$ diverges at an
exponential rate as $n$ increases. Empirically, we demonstrate through an
extensive set of simulations that our estimators substantially outperform the
competitors. Moreover, we integrate our factor modelling framework into the
challenging task of predicting high-dimensional functional time series and use
a real data example to illustrate the superiority of our strategy over
existing prediction methods.

The paper is set out as follows. In Section 2, we develop a fully functional
procedure to estimate the proposed factor model and explain it from a
regression perspective. The theoretical properties of the proposed estimators
are investigated in Section 3. Under further sparsity assumptions, we develop
the regularized estimation procedure and present the associated convergence
analysis in Section 4. The finite-sample performance of our methods is
illustrated through extensive simulation studies and two real data examples in
Sections 5 and 6, respectively. Section 7 concludes the paper by discussing
two extensions. We relegate all technical proofs to the Supplementary
Material.

Notation. For any matrix ${\bf B},$ we let $\|{\bf
B}\|_{\min}=\lambda_{\min}^{1/2}({\bf B}^{{\mathrm{\scriptscriptstyle T}}}{\bf
B})$ and $\|{\bf B}\|=\lambda_{\max}^{1/2}({\bf
B}^{{\mathrm{\scriptscriptstyle T}}}{\bf B}).$ Let $L_{2}(\mathcal{U})$ denote
the Hilbert space of squared integrable functions defined on ${\cal U}$ and
$\otimes$ denotes the Kronecker product. For $f\in L_{2}({\cal U})$, we define
its $\ell_{2}$-norm by $\|f\|=\\{\int f(u)^{2}\mathrm{d}u\\}^{1/2}$. For a
Hilbert space $\mathbb{H}\subset L_{2}({\cal U})$ and a matrix of bivariate
functions ${\bf K}=\big{(}K_{ij}(\cdot,\cdot)\big{)}_{p\times q},$ we define
the functional version of matrix $\ell_{\infty}$ norm by $\|{\bf K}\|_{{\cal
S},\infty}=\max_{1\leq i\leq p}\sum_{j=1}^{q}\|K_{ij}\|_{{\cal S}},$ where,
for each $K_{ij}\in\mathbb{H}\otimes\mathbb{H},$ we denote its Hilbert–Schmidt
norm by $\|K_{ij}\|_{{\cal S}}=\\{\int\int
K_{ij}(u,v)^{2}\,\mathrm{d}u\mathrm{d}v\\}^{1/2}.$ For two positive sequences
$\\{a_{n}\\}$ and $\\{b_{n}\\}$, we write $a_{n}\lesssim b_{n}$ or
$a_{n}=O(b_{n})$ or $b_{n}\gtrsim a_{n}$ if there exist a positive constant
$c$ such that $a_{n}/b_{n}\leq c$. We write $a_{n}\asymp b_{n}$ if and only if
$a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$ hold simultaneously.

##  2 Methodology

###  2.1 Model setup and estimation procedure

We assume that the $p$-vector of functional time series ${\bf
Y}_{1}(\cdot),\dots,{\bf Y}_{n}(\cdot)$ satisfies the functional factor model
in (1), where each ${\bf Y}_{t}(\cdot)$ is decomposed as the sum of two parts:
a dynamic part driven by a $r$-vector of latent factor process ${\bf
X}_{t}(\cdot)$ and a static component of white noise process
$\boldsymbol{\varepsilon}_{t}(\cdot)=\\{\varepsilon_{t1}(\cdot),\dots,\varepsilon_{tp}(\cdot)\\}^{{\mathrm{\scriptscriptstyle
T}}}$ satisfying $E\\{\boldsymbol{\varepsilon}_{t}(u)\\}=\boldsymbol{0}$ and
$\text{Cov}\\{\boldsymbol{\varepsilon}_{t+k}(u),\boldsymbol{\varepsilon}_{t}(v)\\}=\boldsymbol{0}$
for any $u,v\in{\cal U}$ and $k\neq 0.$ This formulation ensures that all
linear dynamic structure of ${\bf Y}_{t}(\cdot)$ are determined by ${\bf
X}_{t}(\cdot),$ while all white noise elements are absorbed into
$\boldsymbol{\varepsilon}_{t}(\cdot).$ Note that $\\{{\bf
X}_{t}(\cdot)\\}_{t=1}^{n}$ and
$\\{\boldsymbol{\varepsilon}_{t}(\cdot)\\}_{t=1}^{n}$ are unobservable. We
further assume that the rank of the factor loading matrix ${\bf A}$ is $r$. If
this full rank condition is violated, model (1) can be expressed in terms of a
lower-dimensional functional factor. When $r$ is much smaller than $p,$ we
achieve an effective dimension reduction.

######  Remark 1.

In functional time series modelling, one crucial question is how to
characterize the functional and dynamical structure. Our factor model (1)
assumes that both are inherited from common factors ${\bf X}_{t}(\cdot)$. The
recent manuscript of Tavakoli et al. (2021) provides another formulation of
factor modeling for multivariate functional time series, treating the factor
loading matrix as a functional object and the common factors as lower-
dimensional scalar time series. Compared with (1), their model leads to
excessive dimension reduction (from $p\times\infty$ to $r$) under the high-
dimensional setting. We claim that our estimation framework can be applied to
fit their model. See the detailed discussion in Section 2.3.

Similar to Lam and Yao (2012), there is an identifiable issue among ${\bf
X}_{t}(\cdot)$ and ${\bf A},$ since model (1) remains unchanged if $\\{{\bf
A},{\bf X}_{t}(\cdot)\\}$ is replaced by $\\{{\bf
A}\boldsymbol{\Gamma},\boldsymbol{\Gamma}^{-1}{\bf X}_{t}(\cdot)\\}$ for any
invertible matrix $\boldsymbol{\Gamma}\in{\mathbb{R}}^{r\times r}.$ We may
then assume that the columns of ${\bf A}$ are orthonormal, that is ${\bf
A}^{{\mathrm{\scriptscriptstyle T}}}{\bf A}={\bf I}_{r}$ (an $r\times r$
identity matrix). With this constraint, ${\bf X}_{t}(\cdot)$ and ${\bf A}$
still can not be determined uniquely, but the linear space spanned by the
columns of ${\bf A}$ (denoted by $\mathcal{C}({\bf A})$) can. Hence we will
focus on the estimation of the factor loading space ${\mathcal{C}}({\bf A}).$

We next develop a fully functional procedure to estimate $\mathcal{C}({\bf
A})$ and the number of factors $r.$ Under the assumption that $\\{{\bf
X}_{t}(\cdot)\\}$ is weakly stationary, we define the following lag-$k$
autocovariance functions:

$\displaystyle\boldsymbol{\Sigma}^{(k)}_{yy}(u,v)$
$\displaystyle=\text{Cov}\\{{\bf Y}_{t+k}(u),{\bf
Y}_{t}(v)\\},~{}~{}\boldsymbol{\Sigma}^{(k)}_{xx}(u,v)=\text{Cov}\\{{\bf
X}_{t+k}(u),{\bf X}_{t}(v)\\},$ (2)
$\displaystyle\boldsymbol{\Sigma}^{(k)}_{x\varepsilon}(u,v)$
$\displaystyle=\text{Cov}\\{{\bf
X}_{t+k}(u),\boldsymbol{\varepsilon}_{t}(v)\\}.$

We further assume that the future white noise components are uncorrelated with
the common factors up to the present. It then follows from (1) and (2) that

$\boldsymbol{\Sigma}^{(k)}_{yy}(u,v)={\bf
A}\boldsymbol{\Sigma}^{(k)}_{xx}(u,v){\bf A}^{{\mathrm{\scriptscriptstyle
T}}}+{\bf A}\boldsymbol{\Sigma}^{(k)}_{x\varepsilon}(u,v),~{}~{}u,v\in{\cal
U},~{}k\geq 1.$ (3)

Given a prescribed and fixed integer $k_{0}\geq 1,$ we define a nonnegative
definite matrix

${\bf M}=\sum_{k=1}^{k_{0}}\int_{{\cal U}}\int_{{\cal
U}}\boldsymbol{\Sigma}^{(k)}_{yy}(u,v){\bf
W}(u,v)\boldsymbol{\Sigma}^{(k)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle
T}}\,\mathrm{d}u\mathrm{d}v,$ (4)

where the $p\times p$ weight matrix ${\bf W}(u,v)$ is nonnegative definite. In
Section 2.2, we will suggest specific forms of ${\bf W}(u,v)$ and illustrate
the rationality from a regression perspective. Let
$\widetilde{\boldsymbol{\gamma}}_{j}$ be a unit eigenvector of ${\bf M}$
corresponding to the $j$-th largest eigenvalue of ${\bf M}.$ We claim that the
factor loading space ${\cal C}({\bf A})$ can be recovered by
$\text{span}\\{\widetilde{\boldsymbol{\gamma}}_{1},\dots,\widetilde{\boldsymbol{\gamma}}_{r}\\}$
and $\text{rank}({\bf M})=r.$ To see this, let $\widetilde{\bf A}$ be a
$p\times(p-r)$ orthogonal complement of matrix ${\bf A}$ such that
$\widetilde{\bf A}^{{\mathrm{\scriptscriptstyle T}}}{\bf A}={\bf 0}$ and
$\widetilde{\bf A}^{{\mathrm{\scriptscriptstyle T}}}\widetilde{\bf A}={\bf
I}_{p-r}.$ It then follows from (3) and (4) that ${\bf M}\widetilde{\bf
A}={\bf 0},$ which implies that the columns of $\widetilde{\bf A}$ are
eigenvectors of ${\bf M}$ corresponding to $(p-r)$ zero eigenvalues and hence
${\cal C}({\bf A})$ is spanned by the eigenvectors of ${\bf M}$ corresponding
to $r$ nonzero eigenvalues.

######  Remark 2.

(i) We take the double integral and the sum in the definition of ${\bf M}$ in
(4) to accumulate the information of autocovariance functions as much as
possible from each $(u,v)\in{\cal U}^{2}$ and from different time lags,
whereas fixing at certain $(u,v)$ or time lag may lead to spurious estimation
results.
(ii) Replacing ${\bf W}(u,v)$ in (4) by its sample version $\widehat{\bf
W}(u,v),$ we define

${\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}=\sum_{k=1}^{k_{0}}\int_{{\cal U}}\int_{{\cal
U}}\boldsymbol{\Sigma}^{(k)}_{yy}(u,v)\widehat{\bf
W}(u,v)\boldsymbol{\Sigma}^{(k)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle
T}}\,\mathrm{d}u\mathrm{d}v.$ (5)

Let $\nu_{1}\geq\dots\geq\nu_{p}\geq 0$ be eigenvalues of
$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$ $\textstyle{\bf M}$ with the corresponding unit
eigenvectors $\boldsymbol{\gamma}_{1},\dots,\boldsymbol{\gamma}_{p},$ which
differ from the eigenpairs of ${\bf M}.$ By the same arguments below (4), we
can easily show that ${\cal C}({\bf A})={\cal C}({\bf
K})=\text{span}\\{\boldsymbol{\gamma}_{1},\dots,\boldsymbol{\gamma}_{r}\\}$
with ${\bf
K}=(\boldsymbol{\gamma}_{1},\dots,\boldsymbol{\gamma}_{r})\in{\mathbb{R}}^{p\times
r}$ and $\nu_{1}\geq\dots\geq\nu_{r}>0=\nu_{r+1}=\dots=\nu_{p}.$
(iii) For the choice of the maximum lag $k_{0},$ we tend to select a small
value in practice as the strongest autocorrelations usually appear at small
lags while larger $k_{0}$ will add more noise to estimate ${\bf M}.$

To estimate ${\cal C}({\bf A}),$ we need to carry out an eigenanalysis of the
natural estimator for ${\bf M},$

$\widehat{\bf M}=\sum_{k=1}^{k_{0}}\int_{{\cal U}}\int_{{\cal
U}}\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)\widehat{\bf
W}(u,v)\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)^{{\mathrm{\scriptscriptstyle
T}}}\,\mathrm{d}u\mathrm{d}v,$ (6)

where

$\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)=\frac{1}{n-k}\sum_{t=1}^{n-k}\big{\\{}{\bf
Y}_{t+k}(u)-\bar{\bf Y}(u)\big{\\}}\big{\\{}{\bf Y}_{t}(v)-\bar{\bf
Y}(v)\big{\\}}^{{\mathrm{\scriptscriptstyle T}}},~{}~{}u,v\in{\cal U},$ (7)

is the estimator for ${\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)$ and $\widebar{\bf
Y}(\cdot)=n^{-1}\sum_{t=1}^{n}{\bf Y}_{t}(\cdot).$ Performing eigen-
decomposition of $\widehat{\bf M}$ leads to estimated eigenvalues
$\hat{\nu}_{1}\geq\dots\geq\hat{\nu}_{p}$ and the corresponding estimated
eigenvectors
$\widehat{\boldsymbol{\gamma}}_{1},\dots,\widehat{\boldsymbol{\gamma}}_{p}.$
Let $\widehat{\bf
K}=(\hat{\boldsymbol{\gamma}}_{1},\dots,\hat{\boldsymbol{\gamma}}_{r}).$ Then
${\cal C}(\widehat{\bf
K})=\text{span}\\{\widehat{\boldsymbol{\gamma}}_{1},\dots,\widehat{\boldsymbol{\gamma}}_{r}\\}$
forms the estimate of ${\cal C}({\bf A}).$

We have developed the estimation of model (1) assuming that the number of
factors $r$ is known or can be identified correctly. In practice, $r$ is
unknown and there is vast literature on the topic of determining it, see,
e.g., Bai and Ng (2002); Onatski (2010); Lam and Yao (2012); Ahn and
Horenstein (2013) and Fan et al. (2020), to quote a few. Here we take the
commonly adopted ratio-based estimator Lam and Yao (2012) for $r$ as:

$\hat{r}=\operatornamewithlimits{argmin}_{1\leq j\leq
c_{r}p}\hat{\nu}_{j+1}/\hat{\nu}_{j},$ (8)

where $c_{r}\in(0,1)$ is a prespecified constant. In empirical studies, we
take $c_{r}=0.75$ to avoid the fluctuations due to the ratios of extreme small
values.

###  2.2 Weight matrix

Under the nonnegative-definiteness constraint, we suggest the following weight
matrix:

${\bf W}(u,v)={\bf Q}\big{\\{}{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\boldsymbol{\Sigma}_{yy}^{(0)}(v,v){\bf Q}\big{\\}}^{-1}{\bf
Q}^{\mathrm{\scriptscriptstyle T}},~{}~{}u\in{\cal U},$ (9)

where ${\bf Q}\in{\mathbb{R}}^{p\times q}$ is a full-rank matrix and $q$ is
larger than $r$ but much smaller than $p.$ The entries of ${\bf Q}$ are
independently generated from $N(0,1)$. Then we can obtain the sample version
of ${\bf W}(u,v)$ by

$\widehat{\bf W}(u,v)={\bf Q}\big{\\{}{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf Q}\big{\\}}^{-1}{\bf
Q}^{\mathrm{\scriptscriptstyle T}}.$ (10)

From a regression perspective, we next provide an intuitive explanation why
the weight matrix takes the suggested form in (9). Despite the unobservable
common factors ${\bf X}_{t}(\cdot),$ we can simply view (1) fixing at
$u\in{\cal U}$ as a multiple linear regression with multivariate responses
${\bf Y}_{t}(u)$ and the coefficient matrix ${\bf A}$ to be estimated. To
construct the covariate vector from observed data in a regression setup, we
assume that there exists some orthogonal matrix
$\boldsymbol{\Phi}\in{\mathbb{R}}^{q\times r}$ with
$\boldsymbol{\Phi}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\Phi}={\bf
I}_{r}$ such that $\widetilde{\bf X}_{t}(\cdot)=\boldsymbol{\Phi}{\bf
X}_{t}(\cdot)$ can be represented by the linear combination of past observed
curves with the addition of some random noise:

$\widetilde{\bf X}_{t}(\cdot)={\bf Q}^{{\mathrm{\scriptscriptstyle T}}}{\bf
Y}_{t-k}(\cdot)+{\mathbf{e}}_{tk}(\cdot),~{}~{}k=1,\dots,k_{0},$

where ${\mathbf{e}}_{tk}(\cdot)$ is the white noise process, independent of
${\bf Y}_{t-1}(\cdot),\dots,{\bf Y}_{t-k_{0}}(\cdot).$ Then model (1) can be
rewritten as

${\bf Y}_{t}(\cdot)={\bf B}\widetilde{\bf
X}_{t}(\cdot)+\boldsymbol{\varepsilon}_{t}(\cdot)={\bf B}{\bf
Q}^{{\mathrm{\scriptscriptstyle T}}}{\bf
Y}_{t-k}(\cdot)+\widetilde{\boldsymbol{\varepsilon}}_{tk}(\cdot),~{}~{}t=k+1,\dots,n,$
(11)

where
$\widetilde{\boldsymbol{\varepsilon}}_{tk}(\cdot)=\boldsymbol{\varepsilon}_{t}(\cdot)+{\bf
B}{\mathbf{e}}_{tk}(\cdot)$ and ${\bf B}={\bf
A}\boldsymbol{\Phi}^{{\mathrm{\scriptscriptstyle T}}}$ is a $p\times q$ matrix
satisfying $\text{rank}({\bf B})=r$ and ${\cal C}({\bf B})={\cal C}({\bf A}).$
Now we can treat (11) fixing at $u\in{\cal U}$ as a linear regression model
with observed covariate vectors $\\{{\bf Q}^{\mathrm{\scriptscriptstyle
T}}{{\bf Y}_{t-k}}(u)\\}_{t=k+1}^{n}$ and unknown low-rank coefficient matrix
${\bf B}.$ Based on (11) and the population moment equation
$\text{Cov}\\{\widetilde{\boldsymbol{\varepsilon}}_{tk}(u),{\bf
Q}^{{\mathrm{\scriptscriptstyle T}}}{\bf Y}_{t-k}(v)\\}=0$ for $u,v\in{\cal
U},$ we can solve ${\bf B}$ by

${\bf B}=\text{Cov}\Big{\\{}{\bf Y}_{t}(u),{\bf
Q}^{{\mathrm{\scriptscriptstyle T}}}{\bf
Y}_{t-k}(v)\Big{\\}}\text{Cov}\Big{\\{}{\bf Q}^{\mathrm{\scriptscriptstyle
T}}{\bf Y}_{t-k}(u),{\bf Q}^{\mathrm{\scriptscriptstyle T}}{\bf
Y}_{t-k}(v)\Big{\\}}^{-1}$ (12)

To simplify our subsequent derivation, we assume that ${\bf X}_{t}(\cdot)$
have mean zero and
$\text{Cov}\\{\widetilde{\boldsymbol{\varepsilon}}_{tk}(u)\\}={\bf I}_{p}$ for
all $t,$ $k$ and $u\in{\cal U}.$ Replacing the covariance terms in (12) by
their sample versions, we obtain the estimator for ${\bf B}$ by

$\widehat{\bf B}=\Big{\\{}\sum_{t=k+1}^{n}{\bf Y}_{t}(u){\bf
Y}_{t-k}^{\mathrm{\scriptscriptstyle T}}(v){\bf
Q}\Big{\\}}\Big{\\{}\sum_{t=k+1}^{n}{\bf Q}^{\mathrm{\scriptscriptstyle
T}}{\bf Y}_{t-k}(u){\bf Y}_{t-k}^{\mathrm{\scriptscriptstyle T}}(v){\bf
Q}\Big{\\}}^{-1}.$ (13)

Denote the $j$-th row vector of $\widehat{\bf B}$ by
${\widehat{\mathbf{b}}}_{j}\in{\mathbb{R}}^{q}$ with invertible covariance
matrix $\text{Cov}(\widehat{\mathbf{b}}_{j})\in{\mathbb{R}}^{q\times q}.$ Our
target is to identify the rank of $\widehat{\bf B}$ and to recover the space
spanned by columns of $\widehat{\bf B},$ which can be carried out in terms of
an eigenanalysis for nonnegative-definite matrix $\widehat{\bf
B}\\{\text{Cov}(\widehat{\mathbf{b}}_{j})\\}^{-1}\widehat{\bf
B}^{{\mathrm{\scriptscriptstyle T}}}=\widetilde{\bf B}\widetilde{\bf
B}^{{\mathrm{\scriptscriptstyle T}}}$ with $\widetilde{\bf B}=\widehat{\bf
B}\\{\text{Cov}(\widehat{\mathbf{b}}_{j})\\}^{-1/2}.$ The presence of
$\text{Cov}(\widehat{\mathbf{b}}_{j})$ accounts for the heterogeneous
variance-covariance effect and leads to the scaled low-rank coefficient matrix
$\widetilde{\bf B}$ with its $j$-th row vector $\widetilde{\mathbf{b}}_{j}$
satisfying $\text{Cov}(\widetilde{\mathbf{b}}_{j})={\bf I}_{q},$
$\text{rank}(\widetilde{\bf B})=\text{rank}(\widehat{\bf B})$ and ${\cal
C}(\widetilde{\bf B})={\cal C}(\widehat{\bf B}).$ It follows from (13) that

$\begin{split}\text{Cov}(\widehat{\mathbf{b}}_{j})&=\Big{\\{}\sum_{t=k+1}^{n}{\bf
Q}^{\mathrm{\scriptscriptstyle T}}{\bf Y}_{t-k}(v){\bf
Y}_{t-k}^{\mathrm{\scriptscriptstyle T}}(u){\bf
Q}\Big{\\}}^{-1}\Big{\\{}\sum_{t=k+1}^{n}{\bf Q}^{\mathrm{\scriptscriptstyle
T}}{\bf Y}_{t-k}(v){\bf Y}_{t-k}^{\mathrm{\scriptscriptstyle T}}(v){\bf
Q}\Big{\\}}\\\ &~{}~{}\cdot\Big{\\{}\sum_{t=k+1}^{n}{\bf
Q}^{\mathrm{\scriptscriptstyle T}}{\bf Y}_{t-k}(u){\bf
Y}_{t-k}^{\mathrm{\scriptscriptstyle T}}(v){\bf Q}\Big{\\}}^{-1}.\end{split}$
(14)

Combining (13) and (14) yields that

$\widehat{\bf
B}\big{\\{}\text{Cov}(\widehat{\mathbf{b}}_{j})\big{\\}}^{-1}\widehat{\bf
B}^{\mathrm{\scriptscriptstyle
T}}=\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v){\bf Q}\\{{\bf
Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf Q}\\}^{-1}{\bf
Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle
T}}.$ (15)

We integrate the right side of (15) over $(u,v)\in{\cal U}^{2}$ and sum it
over time lags $k=1,\dots,k_{0},$ thus obtaining (6) with the choice of
$\widehat{\bf W}(u,v)$ in (10).

An alternative choice of the weight matrix is ${\bf W}_{1}(u,v)={\bf I}_{p},$
where the homogeneous weights are assigned. As a consequence, we need to
perform an eigenanalysis on the estimated nonnegative-definite matrix:

$\displaystyle\widehat{\bf M}_{1}=\sum_{k=1}^{k_{0}}\int_{{\cal U}}\int_{{\cal
U}}\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle
T}}\,\mathrm{d}u\mathrm{d}v,$ (16)

which can be further simplified by integrating along the diagonal path
$u=v\in{\cal U}:$

$\displaystyle\widehat{\bf M}_{2}=\sum_{k=1}^{k_{0}}\int_{{\cal
U}}\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,u)\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,u)^{\mathrm{\scriptscriptstyle
T}}\,\mathrm{d}u.$ (17)

######  Remark 3.

It is worth noting that, without the double or single integral, the estimators
in (16) and (17) coincide with the proposed method for multivariate scalar
time series in Lam et al. (2011) and Lam and Yao (2012). Compared with the
weighted estimator $\widehat{\bf M},$ the performance of the unweighted
estimator $\widehat{\bf M}_{1}$ is expected to deteriorate especially for the
heterogeneous case. Moreover, due to the loss of autocovariance information
for $u\neq v$ incurred by the single integral, we expect that $\widehat{\bf
M}_{1}$ outperforms $\widehat{\bf M}_{2}.$ We will compare the sample
performance of $\widehat{\bf M},$ $\widehat{\bf M}_{1}$ and $\widehat{\bf
M}_{2}$ in our simulations.

###  2.3 Relationship to relevant work

Our estimation procedure can be extended easily to fit the factor model for
multivariate scalar time series Lam and Yao (2012) and, with more efforts, to
deal with the factor modelling for multivariate matrix-valued time series Wang
et al. (2019); Chen et al. (2021). Despite the integration step being no
longer needed, we can still incorporate the suggested weight matrix to account
for the heterogeneous effect.

It is also possible to apply our method to estimate the functional factor
model Tavakoli et al. (2021) in the form of

${\bf Y}_{t}(u)={\bf A}(u){\bf
X}_{t}+\boldsymbol{\varepsilon}_{t}(u),~{}~{}u\in{\cal U},~{}t=1,\dots,n,$
(18)

where ${\bf X}_{t}=(X_{t1},\dots,X_{tr})^{{\mathrm{\scriptscriptstyle T}}}$ is
a set of latent factor time series and ${\bf A}(\cdot)=\\{{\bf
A}_{jl}(\cdot)\\}_{p\times r}$ is the unknown functional factor loading
matrix. To comply with our estimation framework, we assume that the
idiosyncratic component $\boldsymbol{\varepsilon}_{t}(\cdot)$ is white noise.
Given an orthonormal basis $\\{\phi_{i}(\cdot)\\}_{i=1}^{\infty},$ we expand
the functional objects in (18) by
$Y_{tj}(u)=\sum_{i=1}^{\infty}\xi_{tji}\phi_{i}(u),$
$A_{jl}(u)=\sum_{i=1}^{\infty}a_{jli}\phi_{i}(u),$
$\varepsilon_{tj}(u)=\sum_{i=1}^{\infty}\epsilon_{tji}\phi_{i}(u)$ and hence
(18) can be rewritten as

$\xi_{tji}=\sum_{l=1}^{r}a_{jli}X_{tl}+\epsilon_{tji},~{}~{}i=1,\dots,\infty.$
(19)

To simplify notation, we assume the same truncated dimension $M$ across $j$
and stack the basis coefficients $\xi_{tji},$ $a_{jli}$ and $\epsilon_{tji}$
for $j=1,\dots,p$ and $i=1,\dots,M$ to $pM$-dimensional vectors
$\boldsymbol{\xi}_{t},{\mathbf{a}}_{l}$ and $\boldsymbol{\epsilon}_{t},$
respectively. As a result, (19) can be equivalently represented as the
following factor model for $pM$-dimensional vector time series
$\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{n}$:

$\boldsymbol{\xi}_{t}=\boldsymbol{\Lambda}{\bf
X}_{t}+\widetilde{\boldsymbol{\epsilon}}_{t},$ (20)

where the factor loading matrix
$\boldsymbol{\Lambda}=({\mathbf{a}}_{1},\dots,{\mathbf{a}}_{r})\in{\mathbb{R}}^{pM\times
r}$ and the noise vector
$\widetilde{\boldsymbol{\epsilon}}_{t}\in{\mathbb{R}}^{pM}$ is decomposed as
the sum of random and truncation errors. Then we can integrate our suggested
weight matrix into the estimation of $r$ and ${\cal C}(\boldsymbol{\Lambda})$
for the factor model (20). Hence the space spanned by columns of ${\bf
A}(\cdot)$ can be recovered accordingly.

##  3 Theoretical properties

In this section, we study asymptotic properties of the proposed method under a
high-dimensional regime, where $p$ and $n$ tend to infinity together and $r$
is fixed. Before presenting the theoretical results, we impose some regularity
conditions.

######  Condition 1.

(i) The latent functional factor process $\\{{\bf X}_{t}(\cdot)\\}$ is weakly
stationary with ${\mathbb{E}}(\|X_{tj}\|^{4})=O(1)$ for $j=1,\dots,r;$ (ii)
There exists at least one $k\in\\{1,\dots,k_{0}\\}$ such that the rank of
$\iint\boldsymbol{\Sigma}^{(k)}_{xx}(u,v)\boldsymbol{\Sigma}^{(k)}_{xx}(u,v)^{{\mathrm{\scriptscriptstyle
T}}}\,\mathrm{d}u\mathrm{d}v$ is $r.$

######  Condition 2.

(i) The idiosyncratic component $\\{\boldsymbol{\varepsilon}_{t}(\cdot)\\}$ is
white noise with $\max_{j}{\mathbb{E}}(\|\varepsilon_{tj}\|^{2})=O(1);$ (ii)
$\inf_{v\in{\cal
U}}\lambda_{\min}\big{\\{}\boldsymbol{\Sigma}^{(0)}_{\varepsilon\varepsilon}(v,v)\big{\\}}$
is bounded away from zero.

######  Condition 3.

There exists some constant $\delta\in[0,1]$ such that $\|{\bf A}\|^{2}\asymp
p^{1-\delta}\asymp\|{\bf A}\|^{2}_{\min}$.

######  Condition 4.

(i) For $k=1,\dots,k_{0},$
$\|\boldsymbol{\Sigma}^{(k)}_{x\varepsilon}\|_{{{\cal
S}},\infty}=o(p^{(1-\delta)/2});$ (ii) $\text{Cov}\\{{\bf
X}_{t}(u),\boldsymbol{\varepsilon}_{t+k}(v)\\}={\bf 0}$ for any $k\geq 0$ and
$(u,v)\in{\cal U}^{2}.$

Conditions 1 and 2 contain some standard finite moment assumptions in
functional data analysis literature. The parameter $\delta$ in Condition 3 can
be viewed as the strength of factors with smaller values yielding stronger
factors. It measures the relative growth rate of the amount of information
contained in ${\bf Z}_{t}(\cdot)$ as the dimension $p$ increases, compared to
that in $\boldsymbol{\varepsilon}_{t}(\cdot).$ When $\delta=0$, Condition 3
corresponds to the pervasiveness assumption in Fan et al. (2013), which means
that all factors are strong. When $\delta>0,$ the factors are termed as weak
factors. We refer to Lam and Yao (2012) for more detailed discussion on the
factor strength. Condition 4(i) requires that the correlation between ${\bf
X}_{t+k}(\cdot)$ and $\boldsymbol{\varepsilon}_{t}(\cdot)$ is not too strong,
while Condition 4(ii) assumes that the future idiosyncratic components are
uncorrelated with the common factors up to the present.

######  Condition 5.

Let
$\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)=\big{\\{}\widehat{\Sigma}^{(k)}_{yy,ij}(u,v)\big{\\}}_{p\times
p}$ and
${\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)=\big{\\{}\Sigma^{(k)}_{yy,ij}(u,v)\big{\\}}_{p\times
p}$ for $u,v\in{\cal U}.$ Then (i) For $k=1,\dots,k_{0},$
$\max_{i,j}{\mathbb{E}}\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\|_{{\cal
S}}\big{\\}}=O(n^{-1/2});$ (ii) $\boldsymbol{\Sigma}^{(0)}_{yy}(u,v)$ is
Lipschitz-continuous over $(u,v)\in{\cal U}^{2};$ (iii)
$\max_{i,j}{\mathbb{E}}\big{\\{}\sup_{(u,v\in{\cal
U}^{2})}\big{|}\widehat{\Sigma}_{yy,ij}^{(0)}(u,v)-\Sigma_{yy,ij}^{(0)}(u,v)\big{|}\big{\\}}=O\big{\\{}n^{-1/2}(\log
n)^{1/2}\big{\\}}.$

There are several sufficient conditions that have been commonly imposed in
functional time series literature when $p$ is fixed and can lead to the result
in Condition 5(i), which further implies the standard root-$n$ rate for
entrywise sample autocovariance functions, i.e.,
$\|\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\|_{{\cal
S}}=O_{P}(n^{-1/2}).$ The key requirement to establish the consistency is to
control the temporal dependence in $\\{Y_{tj}(\cdot)\\}$ for $j=1,\dots,p.$
Examples include strong mixing conditions Bosq (2000); Chen, Guo and Qiao
(2020), cumulant mixing conditions Panaretos and Tavakoli (2013) and
$L^{q}$-$m$-approximability Hörmann and Kokoszka (2010); Hörmann et al.
(2015). With the additional tail assumption and Lipschitz-continuity in
Condition 5(ii), it is not difficult to apply the partition technique that
reduces the problem from supremum over ${\cal U}^{2}$ to the maximum over a
grid of pairs, and hence the uniform convergence rate of $n^{-1/2}(\log
n)^{1/2}$ in Condition 5(iii) can be achieved. See the same uniform rate with
detailed proof under an independent setting in Qiao et al. (2020).

Now we are ready to present theorems about the rates of convergence for
estimators of the factor loading space ${\cal C}({\bf K})={\cal C}({\bf A})$
and the eigenvalues $\\{\nu_{j}\\}_{j=1}^{p}.$ To measure the accuracy in
estimating ${\cal C}({\bf K}),$ we use the metric of the distance between
${\cal C}({\bf K})$ and ${\cal C}(\widehat{\bf K}).$ For two orthogonal
matrices ${\bf K}_{1}$ and ${\bf K}_{2}$ of dimensions $p\times r_{1}$ and
$p\times r_{2},$ respectively, we define

$\mathcal{D}\Big{(}\mathcal{C}({\bf K}_{1}),\mathcal{C}({\bf
K}_{2})\Big{)}=\Big{\\{}1-\frac{1}{\max(r_{1},r_{2})}\mbox{tr}({\bf K}_{1}{\bf
K}_{1}^{{\mathrm{\scriptscriptstyle T}}}{\bf K}_{2}{\bf
K}_{2}^{{\mathrm{\scriptscriptstyle T}}})\Big{\\}}^{1/2}.$

This distance ranges between $0$ and $1.$ It equals $0$ if and only if
$\mathcal{C}({\bf K}_{1})=\mathcal{C}({\bf K}_{2}),$ and $1$ if and only if
${\bf K}_{1}$ and ${\bf K}_{2}$ are orthogonal. See also Pan and Yao (2008)
and Chang et al. (2018).

######  Theorem 1.

Let Conditions 1–5 hold and $p^{\delta}n^{-1/2}=o(1).$ Then as
$p,n\rightarrow\infty,$ it holds that:

$\displaystyle\mathcal{D}\Big{(}\mathcal{C}({\bf K}),\mathcal{C}(\widehat{\bf
K})\Big{)}=O_{P}(p^{\delta}n^{-1/2}).$

######  Remark 4.

Theorem 1 implies that, as the factors become stronger with smaller $\delta,$
we obtain more efficient estimator for the factor loading space. When the
dimension $p$ is fixed, the $\sqrt{n}$ rate is attained. When all factors are
strong (i.e., $\delta=0$), $\sqrt{n}$ rate still retains, since the signal is
as strong as the noise and hence enhancing $p$ will not affect the estimation
efficiency, circumventing the phenomenon of “curse of dimensionality”. For
weak factors (i.e., $\delta>0$), the noise increases faster than the signal,
and the increase in $p$ will result in a slower convergence rate.

######  Theorem 2.

Let Conditions 1–5 hold and $p^{\delta}n^{-1/2}=o(1).$ Then as
$p,n\rightarrow\infty,$ the following assertions hold:
(i) $|\hat{\nu}_{j}-\nu_{j}|=O_{P}(p^{2-\delta}n^{-1/2})$ for $j=1,\dots,r,$
and
$~{}~{}~{}~{}~{}~{}\hat{\nu}_{j}=O_{P}(p^{2}n^{-1})$ for $j=r+1,\dots,p$;
(ii) $\hat{\nu}_{j+1}/\hat{\nu}_{j}\asymp 1$ for $j=1,\dots,r-1$, and
$\hat{\nu}_{r+1}/\hat{\nu}_{r}=O_{P}(p^{2\delta}n^{-1})$.

######  Remark 5.

We observe from part (i) of Theorem 2 that estimators for nonzero eigenvalues
converge at a slower rate than those for zero eigenvalues. Part (ii) of
Theorem 2 implies that the eigen-ratio $\hat{\nu}_{j+1}/\hat{\nu}_{j}$ will
drop steeply at $j=r,$ thus providing partial theoretical support for the
proposed ratio-based estimator $\hat{r}$ in (8). When all factors are strong
(i.e., $\delta=0$), $\hat{\nu}_{r+1}/\hat{\nu}_{r}=O_{P}(n^{-1}),$ suggesting
that $\hat{r}$ may not suffer from the increase in $p.$ In fact, these results
are consistent to those established in Lam and Yao (2012). In a similar
fashion, we can obtain improved rates as $p$ increases under additional
assumptions. Although we do not theoretically pursue such “blessing of
dimensionality” property here, our simulation results lend empirical support
for this phenomenon.

##  4 Sparse factor model

Despite the phenomenon of “curve of dimensionality” being avoided when all
factors are strong, our method does not guarantee a parsimonious and
interpretable model in the presence of weak factors. In real applications, it
is quite common that many estimated factor loadings are close to zero
especially when $p$ is large, see the gene expression study in Carvalho et al.
(2008) and the sea surface air pressure records example in Lam and Yao (2012).
Such phenomenon motivates us to propose a sparse factor model by imposing
sparsity assumptions on the factor loading space ${\cal C}({\bf A})$ of (1).
In this section, we target to develop the regularized estimation under
sparsity constraints, which not only leads to the enhanced intepretability in
practice, but also theoretically enlarges the dimension of the feature space
that our proposed sparse factor model can handle compared with the nonsparse
factor model with weak factors.

We consider two complementary notions of subspace sparsity defined in terms of
the factor loading matrix ${\bf A}$: row sparsity and column sparsity, which
are consistent with the definitions in Vu and Lei (2013). Intuitively, the row
sparsity entails that only a small subset of components in ${\bf
Y}_{t}(\cdot)$ are driven by the common factors ${\bf X}_{t}(\cdot),$ thus
making the row sparse factor loading space generated by a small number of
variables, independent of the choice of the basis. The column sparsity, on the
other hand, corresponds to the case that each common factor $X_{tj}(\cdot)$
has impact on only a small fraction of components of ${\bf Y}_{t}(\cdot)$ and
hence the column sparse factor loading space has orthonormal basis consisting
of sparse vectors. We begin by introducing parameter spaces of “approximately
column sparse” and “approximately row sparse” factor loading matrices
respectively defined in Conditions 6 and 6′ below. Let ${\mathbf{a}}_{i}$ be
the $i$-th row vector of ${\bf A}=(A_{ij})_{p\times r}.$

######  Condition 6 (Column sparsity).

The matrix ${\bf A}$ is from the class

${\cal V}(\tau,c_{1}(p),L)=\Big{\\{}{\bf A}:\max\limits_{1\leq j\leq
r}\sum_{i=1}^{p}|A_{ij}|^{\tau}\leq c_{1}(p),\ \max\limits_{1\leq i\leq
p,1\leq j\leq r}|A_{ij}|\leq L\Big{\\}},\quad\tau\in[0,1).$

######  Condition 6′ (Row sparsity).

The matrix ${\bf A}$ is from the class

${\cal V}^{*}(\tau,c_{1}(p),L)=\Big{\\{}{\bf
A}:\sum_{i=1}^{p}\|{\mathbf{a}}_{i}\|^{\tau}\leq c_{1}(p),\ \max\limits_{1\leq
i\leq p,1\leq j\leq r}|A_{ij}|\leq L\Big{\\}},\quad\tau\in[0,1).$

The parameters $c_{1}(p)$ and $\tau$ together control the column or row
sparsity in the factor loading matrix ${\bf A}.$ In the special case of
$\tau=0,$ ${\cal V}(0,c_{1}(p),L)$ and ${\cal V}^{*}(0,c_{1}(p),L)$ correspond
to the truly sparse situations, in which ${\bf A}$ has at most $c_{1}(p)$
nonzero entries on each column under Condition 6 or $c_{1}(p)$ nonzero row
vectors under Condition 6′. By Hölder’s inequality, it is easy to check that
${\cal V}^{*}(\tau,c_{1}(p),L)\subseteq{\cal V}(\tau,c_{r,\tau}c_{1}(p),L)$,
where $c_{r,\tau}=1$ for $\tau=0$ and $r^{1-\tau/2}$ for $\tau\in(0,1).$
Provided that $r$ is fixed, our subsequent analysis will focus on the class of
column sparse factor loading matrices, ${\cal V}(\tau,c_{1}(p),L).$

Our estimation procedure in Section 2 is developed based on the estimation of
$\boldsymbol{\Sigma}_{yy}^{(k)}$’s. When $p$ grows larger than $n^{1/2},$ it
follows from Condition 5 that the sample autocovariance functions
$\widehat{\boldsymbol{\Sigma}}_{yy}^{(k)}$’s are no longer consistent.
Nevertheless, under the sparse factor model setup, the decomposition of
$\boldsymbol{\Sigma}_{yy}^{(k)}$ in (3) suggests that our enforced column
sparsity in ${\bf A}$ is inherited by the functional sparsity structure in
$\boldsymbol{\Sigma}_{yy}^{(k)},$ thus enabling us to possibly construct the
threhsolding-based estimator for $\boldsymbol{\Sigma}_{yy}^{(k)}$ to ensure
the consistency. Before introducing the estimator, we slightly modify
Condition 4 to accommodate the sparse setting and present a lemma that reveals
the functional sparsity pattern in $\boldsymbol{\Sigma}_{yy}^{(k)}.$

######  Condition 4′.

(i) For $k=1,\dots,k_{0}$, let
$\boldsymbol{\Sigma}^{(k)}_{x\varepsilon}(u,v)=\big{\\{}\Sigma^{(k)}_{x\varepsilon,ij}(u,v)\big{\\}}_{p\times
p}$ for $u,v\in{\cal U}$ and $\max\limits_{1\leq l\leq
r}\sum_{j=1}^{p}\|\Sigma^{(k)}_{x\varepsilon,lj}\|^{\tau}_{{\cal S}}\leq
c_{2}(p);$ (ii) $\text{Cov}\\{{\bf
X}_{t}(u),\boldsymbol{\varepsilon}_{t+k}(v)\\}=0$ for any integer $k\geq 0$
and $u,v\in{\cal U}.$

Condition 4′(i) requires relatively weak correlations between ${\bf
X}_{t+k}(\cdot)$ and $\boldsymbol{\varepsilon}_{t}(\cdot).$ As long as
$c_{2}(p)=o(p^{(1-\delta)/2}),$ Condition 4(i) follows directly from Condition
4′(i). Moreover, Condition 4′(i) relies on the Hilbert–Schmidt norm to
encourage the functional sparsity in
$\boldsymbol{\Sigma}_{x\varepsilon}^{(k)}$, i.e., each common factor
$X_{(t+k)l}(\cdot)$ is only correlated with a few components of
$\boldsymbol{\varepsilon}_{t}(\cdot).$ With imposed sparsity constraints in
Conditions 4′ and 6, it can be inferred from (3) that
$\boldsymbol{\Sigma}^{(k)}_{yy}$ is functional sparse as justified in Lemma 1
below.

######  Lemma 1.

Under Conditions 1, 2, 4′ and 6, it holds that, for $k=1,\dots,k_{0},$

$\displaystyle\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\|\Sigma^{(k)}_{yy,ij}\|^{\tau}_{{\cal S}}\lesssim
c_{1}(p)~{}~{}\text{and}~{}~{}\max\limits_{1\leq i\leq
p}\sum_{j=1}^{p}\|\Sigma^{(k)}_{yy,ij}\|^{\tau}_{{\cal S}}\lesssim
c_{1}(p)+c_{2}(p).$

Lemma 1 shows that the functional sparsity patterns in columns/rows of
$\boldsymbol{\Sigma}_{yy}^{(k)}$ are determined by parameters $c_{1}(p)$ and
$c_{2}(p)$ with smaller values yielding functional sparser
$\boldsymbol{\Sigma}_{yy}^{(k)}.$ To obtain a functional sparse estimator for
$\boldsymbol{\Sigma}_{yy}^{(k)},$ we apply the hard functional thresholding
rule, which combines functional versions of hard thresholding and shrinkage
based on the Hilbert-Schmidt norm of functions, on entries of the sample
autocovariance function $\widehat{\boldsymbol{\Sigma}}_{yy}^{(k)}.$ Then the
functional thresholding estimator is constructed as

$\mathcal{T}_{\eta_{k}}(\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy})(u,v)=\
\Big{[}\widehat{\Sigma}^{(k)}_{yy,ij}(u,v)I\big{\\{}\|\widehat{\Sigma}^{(k)}_{yy,ij}\big{\|}_{{\cal
S}}\geq\eta_{k}\big{\\}}\Big{]}_{p\times p},\quad u,v\in{\cal U},$ (21)

where $I(\cdot)$ is the indicator function and $\eta_{k}\geq 0$ is the
thresholding parameter. Under mild regularity conditions, it follows from
Lemma B3 in the Supplementary Material that

$\max_{1\leq i,j\leq
p}\|\widehat{\Sigma}^{(k)}_{yy,ij}-\Sigma^{(k)}_{yy,ij}\|_{{\cal
S}}=O_{P}\left({\cal M}_{y}\sqrt{\frac{\log p}{n}}\right),$ (22)

where ${\cal M}_{y}$ is the functional stability measure Guo and Qiao (2021);
Fang et al. (2021) defined in (A.1) in Appendix A. The rate under the
functional version of $\ell_{\infty}$ norm in (22) plays a crucial role in our
theoretical analysis under an ultra-high-dimensional regime and, in
particular, suggests us to set the thresholding level as $\eta_{k}\asymp{\cal
M}_{y}(n^{-1}\log p)^{1/2}.$ Replacing
$\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}$ in (6) with ${\cal
T}_{\eta_{k}}(\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}),$ we obtain the
corresponding estimator for ${\bf M}$ by

$\widetilde{\bf
M}=\sum_{k=1}^{k_{0}}\int\int\mathcal{T}_{\eta_{k}}(\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy})(u,v)\widehat{\bf
W}(u,v)\mathcal{T}_{\eta_{k}}(\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy})(u,v)^{{\mathrm{\scriptscriptstyle
T}}}\,\mathrm{d}u\mathrm{d}v.$ (23)

To recover the column sparsity structure in ${\cal C}({\bf A}),$ we perform
sparse principal component analysis (PCA) Vu and Lei (2013) on $\widetilde{\bf
M}$ rather than the standard eigenanalysis for $\widehat{\bf M}$ in Section
2.1. For matrices ${\bf A}_{1}$ and ${\bf A}_{2}$ with the same dimension, let
$\langle{\bf A}_{1},{\bf A}_{2}\rangle:=\text{trace}({\bf
A}_{1}^{\mathrm{\scriptscriptstyle T}}{\bf A}_{2}).$ We define $\widetilde{\bf
K}\in\mathbb{R}^{p\times r}$ as a solution to the following constraint
optimization problem:

$\widetilde{\bf K}=\underset{{\bf K}=(K_{jl})_{p\times
r}}{\arg\max}~{}~{}\langle\widetilde{{\bf M}},{\bf K}{\bf
K}^{{\mathrm{\scriptscriptstyle T}}}\rangle~{}~{}\text{subject to}~{}~{}{\bf
K}^{\mathrm{\scriptscriptstyle T}}{\bf K}={\bf I}_{r},~{}\max\limits_{1\leq
l\leq r}\sum_{j=1}^{p}|K_{jl}|^{\tau}\leq C_{\tau},$ (24)

where $C_{\tau}>0$ is a regularization parameter. Alternatively, to estimate
the row sparse factor loading space, we can substitute the second constraint
in (24) by
$\sum_{i=1}^{p}\|{\mathbf{k}}_{i}\|^{\tau}\leq\widetilde{C}_{\tau},$ where
${\mathbf{k}}_{i}$ denotes the $i$-th row vector of ${\bf K}$ and
$\widetilde{C}_{\tau}>0$ is a regularization parameter. It is worth noting
that, without the sparsity constraint in (24), the optimization problem
degenerates to the ordinary PCA. Despite being challenging to solve (24) due
to the non-convex constraint, some efficient and computationally tractable
algorithms have been developed, see, e.g., under the truly sparse case
($\tau=0$), the combinatorial approaches Moghaddam et al. (2006); d’Aspremont
et al. (2008); Mackey (2009), the semi-definite relaxation d’Aspremont et al.
(2007) and its variants, and the random-projection-based method Gataric et al.
(2020). We refer to Zou and Xue (2018) for an extensive review on recent
developments for sparse PCA.

We now present the asymptotic analysis of ${\cal C}(\widetilde{\bf K})$ in the
following theorem.

######  Theorem 3.

Let Conditions 1, 2, 4′, 6 and 7–8 in Appendix A hold and ${\cal
M}_{y}^{2}\log pn^{-1}=o(1).$ Then as $p,n\rightarrow\infty,$ it holds that:

$\nu_{r}\mathcal{D}\Big{(}\mathcal{C}({\bf K}),\mathcal{C}({\widetilde{\bf
K}})\Big{)}=O_{P}\left[c_{1}(p)\big{\\{}c_{1}(p)+c_{2}(p)\big{\\}}\mathcal{M}_{y}^{1-\tau}\Big{(}\frac{\log
p}{n}\Big{)}^{\frac{1-\tau}{2}}\right].$

######  Remark 6.

(i) The convergence rate of ${\mathcal{D}}\big{(}\mathcal{C}({\bf
K}),\mathcal{C}({\widetilde{\bf K}})\big{)}$ is governed by both
dimensionality parameters $\\{n,p,c_{1}(p),c_{2}(p)\\}$ and internal
parameters (${\cal M}_{y},\nu_{r},\tau$). It is easy to see that the rate is
better when $\nu_{r}$ is large and $\\{c_{1}(p),c_{2}(p),{\cal M}_{y},\tau\\}$
are small.
(ii) Under the truly sparse case ($\tau=0$) with $c_{1}(p)\gtrsim c_{2}(p),$
and $|A_{ij}|\asymp\gamma$ for $(i,j)$ such that $A_{ij}\neq 0,$ it follows
from the framework in Section 3 that $c_{1}(p)\gamma^{2}\asymp p^{1-\delta}$
under Condition 3, which together with Theorem 3 imply that

$\nu_{r}\mathcal{D}\big{(}\mathcal{C}({\bf K}),\mathcal{C}({\widetilde{\bf
K}})\big{)}=O_{P}\big{\\{}p^{2-2\delta}\gamma^{-4}{\cal M}_{y}(\log
p)^{1/2}n^{-1/2}\big{\\}}.$ (25)

By comparison, Theorem 1 and (S.3) of the Supplementary Material lead to
$\nu_{r}{\mathcal{D}}\big{(}\mathcal{C}({\bf K}),\mathcal{C}({\widetilde{\bf
K}})\big{)}=O_{P}\big{(}p^{2-\delta}n^{-1/2}\big{)},$ which is slower than the
rate in (25) for larger values of $\gamma$ or $\delta$ (i.e., smaller values
of $c_{1}(p)$ provided that $\gamma\asymp 1$). Hence, when the magnitudes of
nonzero entries in ${\bf A}$ become larger or the factors are weaker in the
sense of Condition 3 (i.e., ${\bf A}$ is sparser in the sense of Condition 6),
our regularized estimation benefits more from the imposed sparsity and enjoys
faster convergence rate than the ordinary method in Section 2.

##  5 Simulation studies

###  5.1 Setup

We illustrate the finite-sample performance of our proposed methods through
extensive simulations. Sections 5.2 and 5.3 consider scenarios where the
factor loading matrix ${\bf A}$ is ordinary and sparse, respectively.

In each simulated scenario, we generate $p$-vector of functional time series
by

${\bf Y}_{t}(\cdot)=\kappa_{0}{\bf A}\widetilde{\bf
X}_{t}(\cdot)+\boldsymbol{\varepsilon}_{t}(\cdot),~{}~{}t=1,\dots,n,$ (26)

where the parameter $\kappa_{0}>0$ controls the strength of common factors
${\bf X}_{t}(\cdot)=\kappa_{0}\widetilde{\bf X}_{t}(\cdot)$ and the entries of
${\bf A}\in\mathbb{R}^{p\times r}$ are sampled from
$\text{Uniform}[-\sqrt{3}p^{-\delta/2},\sqrt{3}p^{-\delta/2}]$ with
$\delta\in[0,1].$ Hence Condition 3 is satisfied, in which $\delta=0$ (or
$\delta>0$) corresponds to the case of strong (or weak) factors. To mimic the
infinite-dimensionality of functional data, we generate each scaled latent
factor by $\widetilde{X}_{tl}(\cdot)=\sum_{i=1}^{50}\xi_{tli}\phi_{i}(\cdot)$
for $l=1,\dots,r$ over ${\cal U}=[0,1],$ where
$\\{\phi_{i}(\cdot)\\}_{i=1}^{50}$ is a $50$-dimensional Fourier basis
function and the basis coefficients
$\boldsymbol{\xi}_{ti}=(\xi_{t1i},\dots,\xi_{tri})^{{\mathrm{\scriptscriptstyle
T}}}$ are generated from a vector autoregressive model,
$\boldsymbol{\xi}_{ti}={\bf
V}\boldsymbol{\xi}_{(t-1)i}+\boldsymbol{\epsilon}_{i}$ with ${\bf
V}=(\rho^{|l-l^{\prime}|+1})_{1\leq l,l^{\prime}\leq r}$ and the innovation
$\boldsymbol{\epsilon}_{ti}=(\epsilon_{t1i},\dots,\epsilon_{tri})^{{\mathrm{\scriptscriptstyle
T}}}$ consisting of independent $N(0,i^{-2\iota})$ components. We set
$\rho=0.45$ and $\iota=0.75.$ For the idiosyncratic component
$\boldsymbol{\varepsilon}_{t}(\cdot),$ we consider the following three
scenarios.

  * Scenario 1.

For each $j=1,\dots,p,$ we generate
$\varepsilon_{tj}(\cdot)=\sum_{i=1}^{20}2^{-(i+1)}Z_{tji}\phi_{i}(\cdot),$
where $Z_{tji}$’s are independent standard normal.

  * Scenario 2.

We generate $\boldsymbol{\varepsilon}_{t}(\cdot)={\bf
H}\tilde{\boldsymbol{\varepsilon}}_{t}(\cdot),$ where ${\bf
H}=5^{-1}\text{diag}(h_{1},\dots,h_{p}),$ $h_{j}$’s are sampled uniformly from
$\\{1,\dots,10\\}$ and each $\tilde{\varepsilon}_{tj}(\cdot)$ is generated in
the same way as $\varepsilon_{tj}(\cdot)$ in Scenario 1.

  * Scenario 3.

We fix $\kappa_{0}=1$ and generate
$\varepsilon_{tj}(\cdot)=\kappa_{1}\varepsilon_{t0}(\cdot)+\tilde{\varepsilon}_{tj}(\cdot),$
where $\kappa_{1}>0,$
$\varepsilon_{t0}(\cdot)=\sum_{i=1}^{4}\widetilde{Z}_{ti}\phi_{i}(\cdot)$ and
each $\widetilde{Z}_{ti}$ is sampled independently from $N(0,1).$

In Scenario 1, each $\varepsilon_{tj}(\cdot)$ is white noise with the
identical covariance function across $j,$ while Scenario 2 corresponds to the
heterogeneous case with different covariance functions. Scenario 3 consists of
$(r+1)$ actual common factors, among which the additional factor
$\varepsilon_{t0}(\cdot)$ is independent of other factors in ${\bf
X}_{t}(\cdot),$ and its signal strength is determined by the parameter
$\kappa_{1}.$ The underlying factor loading matrix is $({\bf A},{\bf 1}_{p}),$
where ${\bf 1}_{p}$ denotes the $p$-vector of ones. For simplicity, we still
denote the factor loading matrix by ${\bf A}.$

###  5.2 Ordinary case

We compare our proposed method based on $\widehat{\bf M}$ with two competing
methods based on $\widehat{\bf M}_{1}$ and $\widehat{\bf M}_{2}.$ Since our
simulation results suggest that $\hat{r}$ and $\mathcal{C}(\widehat{\bf K})$
are insensitive to the choice of $q$ and the maximum lag $k_{0},$ we set
$q=12$ and $k_{0}=4$ in our simulations. In each setting, we generate $n=100$
serially dependent observations of $p=50,100,200$ functional variables based
on $r=4$ functional factors. We ran each simulation $100$ times. The sample
performance of three approaches is examined in terms of their abilities of
correctly identifying the number of factors and the estimation accuracy in
recovering the factor loading space, respectively measured by the relative
frequency estimate for $P(\hat{r}=r)$ and $\mathcal{D}\big{(}\mathcal{C}({\bf
K}),\mathcal{C}(\widehat{\bf K})\big{)}$ using the correct $r.$ For each of
three comparison methods, Figures 1 and 2 plot average relative frequencies
$\hat{r}=r$ as the factor strength increases when the factors are strong
(i.e., $\delta=0$) and weak with $\delta=0.5$, respectively. Figures 3 and 4
plot the corresponding average estimation errors for $\mathcal{C}({\bf K}).$

(a) $p=50$

(b) $p=50$

(c) $p=50$

(d) $p=100$

(e) $p=100$

(f) $p=100$

(g) $p=200$

(h) $p=200$

(i) $p=200$

Figure 1: Scenario 1 (top row), Scenario 2 (middle row) and Scenario 3 (bottom
row) for $p=50,$ $100$ and $200$ when the factors are strong: Plots of average
relative frequency estimates for $P(\hat{r}=r)$ against the factor strength
($\kappa_{0}$ or $\kappa_{1}$) for three methods based on $\widehat{\bf M}$
(black solid), $\widehat{\bf M}_{1}$ (red dashed) and $\widehat{\bf M}_{2}$
(cyan dash dotted).

(a) $p=50$

(b) $p=50$

(c) $p=50$

(d) $p=100$

(e) $p=100$

(f) $p=100$

(g) $p=200$

(h) $p=200$

(i) $p=200$

Figure 2: Scenario 1 (top row), Scenario 2 (middle row) and Scenario 3 (bottom
row) for $p=50,$ $100,$ $200$ when the factors are strong: Plots of average
$\mathcal{D}\big{(}\mathcal{C}({\bf K}),\mathcal{C}(\widehat{\bf K})\big{)}$
against $\kappa_{0}$ or $\kappa_{1}$ for three methods based on $\widehat{\bf
M}$ (black solid), $\widehat{\bf M}_{1}$ (red dashed) and $\widehat{\bf
M}_{2}$ (cyan dash dotted).

Several conclusions can be drawn from Figures 1–4. First, our proposed method
based on $\widehat{\bf M}$ provides highly significant improvements in
accuracies for the identification of $r$ and the recovery of ${\cal C}({\bf
K})$ over the competing methods in all scenarios we consider. The improvement
of our method involving the weight matrix is larger for the heterogeneous case
and even more substantial for the case when the factors are weak. Between two
competitors, the method based on $\widehat{\bf M}_{1}$ outperforms that based
on $\widehat{\bf M}_{2},$ providing empirical evidence for Remark 3. Second,
the estimation for $r$ and ${\cal C}({\bf K})$ performs better as the strength
of factors increases (i.e., $\kappa_{0}$ or $\kappa_{1}$ increases or $\delta$
decreases), which is in line with our theoretical results in Section 3. In
Figure 1, our method makes the sharpest progress as the factors become
stronger, while the two comparison methods require much higher strength of
factors to compete. Third, when the factors are strong, we observe the
phenomenon of “blessing of dimensionality” in Figures 1 and 2 in the sense
that the estimation improves as $p$ increases from $50$ to $200.$ The
improvement is due to the increase of the information from added components on
the factors. Under Scenario 3 with an extra factor independent of the others,
the information and noise on the factors increase simultaneously as $p$ grows,
and hence the estimation does not necessarily improve. Fourth, under the weak
factor setting, many entries of ${\bf A}$ are quite close to zero. As $p$
enlarges, the enhanced information on the factors is accompanied with the
increase of noise. E.g., in Figure 3, we observe that enlarging $p$ does not
necessarily lead to the increased estimation accuracy.

(a) $p=50$

(b) $p=50$

(c) $p=100$

(d) $p=100$

(e) $p=200$

(f) $p=200$

Figure 3: Scenario 1 (top row) and Scenario 2 (bottom row) for $p=50,$ $100$
and $200$ with correct $r$ when the factors are weak: Plots of relative
frequency estimates for $P(\hat{r}=r)$ against $\kappa_{0}$ for three methods
based on $\widehat{\bf M}$ (black solid), $\widehat{\bf M}_{1}$ (red dashed)
and $\widehat{\bf M}_{2}$ (cyan dash dotted).

(a) $p=50$

(b) $p=50$

(c) $p=100$

(d) $p=100$

(e) $p=200$

(f) $p=200$

Figure 4: Scenario 1 (top row) and Scenario 2 (bottom row) for $p=50,$ $100$
and $200$ with correct $r$ when the factors are weak: Plots of average
$\mathcal{D}\big{(}\mathcal{C}({\bf K}),\mathcal{C}(\widehat{\bf K})\big{)}$
against $\kappa_{0}$ for three methods based on $\widehat{\bf M}$ (black
solid), $\widehat{\bf M}_{1}$ (red dashed) and $\widehat{\bf M}_{2}$ (cyan
dash dotted).

###  5.3 Sparse Case

In this section, we conduct simulations to evaluate the performance of the
functional-thresholding-and-sparse-PCA-based approach (TSPCA) developed in
Section 4 to estimate the sparse factor model. Two kinds of subspace sparsity
constraints are imposed on ${\bf A}.$ Specifically, for the row sparsity, we
randomly select $80\%$ rows of ${\bf A}$ to be zero vectors, while for the
column sparsity, we randomly set $80\%$ elements within each column of ${\bf
A}$ as zeros. We focus on Scenario 1 with $\kappa_{0}=1$ and generate nonzero
entries of ${\bf A}$ from $\text{Uniform}[-\sqrt{3},\sqrt{3}].$

Implementing TSPCA requires choosing the thresholding parameter $\eta_{k}$ in
(21) and the columnwise number of nonzero elements $C_{0}$ in (24). To select
the optimal $\hat{\eta}_{k}$ for each $k=1,\dots,k_{0},$ we implement a
$G$-fold cross-validation (CV) approach Cai and Liu (2011). To be specific, we
first sequentially divide the set $\\{1,\dots,n\\}$ into $G$ blockwise groups
${\cal D}_{1},\dots,{\cal D}_{G}$ of approximately equal size. We then treat
the $g$-th group as a validation set, compute the sample lag-$k$
autocovariance functions $\widehat{\boldsymbol{\Sigma}}_{yy}^{(k),(g)}(u,v)$
and $\widehat{\boldsymbol{\Sigma}}_{yy}^{(k),(-g)}(u,v)$ for $u,v\in{\cal U}$
based on the validation set and the remaining $G-1$ groups, respectively, and
repeat the above procedure $G$ times. We finally select $\widehat{\eta}_{k}$
by minimizing
$R_{k}(\eta_{k})={G}^{-1}\sum_{g=1}^{G}\big{\|}{\mathcal{T}}_{\eta_{k}}\big{(}\widehat{\boldsymbol{\Sigma}}^{(k),(-g)}_{yy}\big{)}-\widehat{\boldsymbol{\Sigma}}^{(k),(g)}_{yy}\big{\|}_{{\cal
S},\text{F}}^{2}.$ We adopt a similar CV method to select the optimal
$\widehat{C}_{0}.$ Given the $g$-th group as a validation set, we obtain the
solution $\widetilde{\bf K}^{(-g)}$ to the constrained optimization problem in
(24) with $\tau=0$ based on the remaining $G-1$ groups. To solve this problem,
we apply sparse PCA Moghaddam et al. (2006); Mackey (2009) to $\widetilde{\bf
M}^{(-g)},$ which is formed by (23) using
${\mathcal{T}}_{\widehat{\eta}_{k}}\big{(}\widehat{\boldsymbol{\Sigma}}^{(k),(-g)}_{yy}\big{)}$
for $k=1,\dots,k_{0}.$ We also obtain $\widehat{\bf K}^{(g)}$ by carrying out
an eigenanalysis for $\widehat{\bf M}^{(g)}$ based on the validation set. The
above procedure is repeated $G$ times and $\widehat{C}_{0}$ is selected by
minimizing $D(C_{0})={G}^{-1}\sum_{g=1}^{G}{\cal D}\big{(}{\cal
C}(\widetilde{\bf K}^{(-g)}),{\cal C}(\widehat{\bf K}^{(g)})\big{)}.$ Although
the time break created by each leave-out validation set possibly jeopardize
the autocovariance structure on the remaining $G-1$ groups via $k_{0}$ mis-
utilized lagged terms, its effect on the estimation of
$\boldsymbol{\Sigma}_{yy}^{(k)(-g)}$’s is negligible especially when $n$ is
sufficiently large. Hence our proposed CV approach does not place a practical
constraint.

Table 1: The mean and standard error (in parentheses) of the distance between
the estimated and true sparse factor loading spaces over $100$ simulation
runs. All entries have been multiplied by $10^{2}$ for formatting reasons.

Sparsity | $p$ | PCA | TSPCA | Sparsity | $p$ | PCA | TSPCA
---|---|---|---|---|---|---|---
| $100$ | $3.86(0.53)$ | $1.71(0.31)$ |  | $100$ | $4.27(0.83)$ | $3.55(1.09)$
Row | $200$ | $3.79(0.56)$ | $1.74(0.26)$ | Column | $200$ | $4.25(0.65)$ | $3.69(0.59)$
| $400$ | $3.73(0.66)$ | $1.76(0.29)$ |  | $400$ | $4.03(0.74)$ | $3.17(0.88)$

We compare TSPCA with the ordinary PCA-based approach by performing an
eigenanlysis on $\widehat{\bf M}$. The estimation quality is measured in terms
of the distance between the estimated and true sparse factor loading spaces
using the correct $r.$ We report the numerical summaries in Table 1 for
$p=100,200,400.$ It is obvious to see that TSPCA uniformly improves the PCA-
based estimation, demonstrating the advantage of our regularized estimation
procedure to fit sparse factor models. The improvement for the row sparsity
case is more significant than that for the column sparsity case. Compared with
the functional sparsity patterns in Lemma 1 under the column sparsity
constraint, it follows from the decomposition of
$\boldsymbol{\Sigma}_{yy}^{(k)}$ in (3) that the row sparsity in ${\bf A}$ can
lead to functional sparser $\boldsymbol{\Sigma}_{yy}^{(k)}$’s. In this sense,
TSPCA benefits more from the row sparsity structure, thus resulting in
enhanced improvement.

##  6 Real data analysis

###  6.1 U.K. temperature data

Our first dataset, available at
https://www.metoffice.gov.uk/research/climate/maps-and-data/historic-station-
data, consists of monthly average temperature collected at $p=22$ measuring
stations in the U.K. from 1959 to 2020 ($n=62$). Let $Y_{ij}(u_{k})$
($i=1,\dots,n,j=1,\dots,p,k=1,\dots,12$) denotes the average temperature
during month $u_{k}=k\in{\cal U}=[1,12]$ of year $1958+i$ at the $j$-th
measuring station. The observed temperature curves are smoothed based on
$10$-dimensional Fourier basis to capture the periodic structures over the
annual cycle. The smoothed curve series exhibits very small autocorrelations
beyond $k=2,$ so we use $k_{0}=2$ in computing $\widehat{\bf M}.$ Although the
ratio-based estimator suggests $\hat{r}=1,$ we illustrate further results
using $\hat{r}=3.$

Figure 5 displays spatial heatmaps of estimated factor loading matrix and the
rotated matrix using the varimax procedure to maximize the sum of the
variances of the squared loadings. Some interesting patterns can be observed
from those heatmaps. Compared with the original results in Figure 5, the
varimax rotation brings the factor loading matrix closer to a “simple
structure”, where (i) each component in ${\bf Y}_{t}(\cdot)$ has a high
loading on one specific factor but near-zero loadings on other factors and
(ii) each factor has high impact on a few components of ${\bf Y}_{t}(\cdot)$
with high loadings on this factor while the remaining variables have near-zero
loadings on this factor. The varimax rotation leads to larger color variations
in heatmaps for better interpretation. For example, it is apparent that the
first factor via varimax rotation influences the dynamics in the southeast,
while the original factor has uniform impact in all locations. Hence we focus
on the interpretation of remaining factors after varimax rotation.
Specifically, the second factor mainly impacts the dynamics of the northern
region. The third factor can roughly be viewed as the main driving force for
the dynamics in the middle north.

(a) 1st factor (original)

(b) 1st factor (rotated)

(c) 2nd factor (original)

(d) 2nd factor (rotated)

(e) 3rd factor (original)

(f) 3rd factor (rotated)

Figure 5: Spatial heatmaps based on estimated factor loading matrix (top row)
and varimax-rotated loading matrix (bottom row) of 22 U.K. locations.

###  6.2 Japanese mortality data

Our second dataset contains age-specific and gender-specific mortality rates
for $p=47$ Japanese prefectures from 1975 to 2017 ($n=43$). This dataset was
also analyzed in Gao et al. (2019). Due to those sparse observations at old
ages, we focus on data for ages below 96. We apply a log transformation to
mortality rates and denote by $Y_{ij}(u_{k})$
($i=1,\dots,n,j=1,\dots,p,k=1,\dots,96$) the log mortality rate of people aged
$u_{k}=k-1\in{\cal U}=[0,95]$ living in the $j$-th prefecture during the year
$1974+i.$ We then perform smoothing for observed curves and replace the
missing values via smoothing splines. The estimation of model (1) is done
using $k_{0}=2.$ The ratio-based estimator for $r$ is $\hat{r}=1,$ though we
use $\hat{r}=2$ for the illustration. To enhance interpretability for
identified factors under a high-dimensional $p>n$ regime, we also implement
TSPCA to estimate the sparse factor loading matrix.

Figures 6 and 7 in the Supplementary Material show spatial heatmaps of
varimax-rotated loading matrix and sparse loading matrices with different
sparsity levels for Japanese females and males, respectively. For TSPCA, we
implement the CV method to select the optimal $\hat{\eta}_{k}$’s and set the
columnwise sparsity in estimated loading matrix to 27/47 or 32/47 for better
visualization of the results. A more systematic method to determine the
sparsity level, e.g., via a significance testing, needs to be developed.
Compared with the original results, both varimax rotation and TSPCA lead to
enhanced interpretability in the sense of recovering the factor loading matrix
with “simple structure”. Among two competitors, TSPCA tends to reduce some
near-zero loadings to exactly zero, thus providing a parsimonious model with
more interpretable results than varimax rotation. In particular, we can
observe clear regional patterns from those heatmaps for TSPCA, where each of
two factors influences the dynamics of some complementary regions of Japan in
all settings. Take the results under the sparsity 27/47 in Figure 6 as an
example, the first factor serves as the main driving force in regions of
Hokkaido, Kanto, Chubu and Kansal, while Tohoku and Shikoku are heavily loaded
regions on the second factor.

(a) 1st factor (rotated)

(b) 2nd factor (rotated)

(c) 1st factor (27 zeros)

(d) 2nd factor (27 zeros)

(e) 1st factor (32 zeros)

(f) 2nd factor (32 zeros)

Figure 6: Spatial heatmaps based on varimax-rotated loadings (left column) and
sparse loadings with 27 zeros (middle column) and 32 zeros (right column) of
47 prefectures on two factors for Japanese females.

To further illustrate the developed methodology, we set upon the task of
predicting high-dimensional functional time series ${\bf Y}_{t}(\cdot).$
Specifically, we incorporate the factor model framework into the $h$-step-
ahead prediction (denoted as FMP) consisting of three steps below.

  1. 1.

Apply our proposed method to estimate model (1) based on past observations
$\\{{\bf Y}_{t}(\cdot)\\}_{t=1}^{T},$ thus obtaining the estimated factor
loading matrix $\widehat{\bf A}$ and number of factors $\hat{r}.$

  2. 2.

Compute $\hat{r}$ estimated factors by $\widehat{\bf
X}_{t}(\cdot)=\widehat{\bf A}^{\mathrm{\scriptscriptstyle T}}{\bf
Y}_{t}(\cdot).$ For each $k\in\\{1,\dots,\hat{r}\\},$ predict
$\widehat{X}_{(t+h)k}(\cdot)$ based on past fitted values for the $k$-th
common factor, $\\{\widehat{X}_{tk}(\cdot)\\}_{t=1}^{T}.$

  3. 3.

The $h$-step head predict for ${\bf Y}_{T+h}(\cdot)$ is $\widehat{\bf
Y}_{T+h}(\cdot)=\widehat{\bf A}\widehat{\bf X}_{T+h}(\cdot).$

The second step is on the prediction of univariate functional time series. For
each common factor, we obtain the $h$-step-ahead prediction based on the best
fitted functional ARMA model Klepsch et al. (2017) according to the BIC
criterion. We then develop a sparse version of FMP (denoted as SFMP) by
performing TSPCA in the first step to estimate sparse $\widehat{\bf A}.$ For
comparison, we also implement the method of Gao et al. (2019) (denoted as GSY)
to jointly predict $p$-vector of functional time series ${\bf Y}_{t}(\cdot)$,
as well as the univariate prediction method Aue et al. (2015) (denoted as UP)
on each $Y_{tj}(\cdot)$ separately.

To evaluate the predictive accuracy, we use the expanding window approach. The
data is divided into a training set and a test set consisting of the first
$n_{1}$ and the last $n_{2}$ observations, respectively. For any integer
$h>0,$ we implement each fitting method on the training set, obtain $h$-step-
ahead prediction on the test data based on the fitted model, increase the
training size by one and repeat this procedure $n_{2}-h+1$ times to compute
the $h$-step-ahead mean absolute prediction error (MAPE) and mean squared
prediction error (MSPE) by

$\text{MAPE}(h)=\frac{1}{p\times(n_{2}-h+1)\times
96}\sum_{j=1}^{p}\sum_{t=n_{1}+h}^{n}\sum_{k=1}^{96}\Big{|}\widehat{Y}_{tj}(u_{k})-Y_{tj}(u_{k})\Big{|},$
$\text{MSPE}(h)=\frac{1}{p\times(n_{2}-h+1)\times
96}\sum_{j=1}^{p}\sum_{t=n_{1}+h}^{n}\sum_{k=1}^{96}\Big{\\{}\widehat{Y}_{tj}(u_{k})-Y_{tj}(u_{k})\Big{\\}}^{2}.$

The resulting MAPE and MSPE values are summarized in Figure 2. It is obvious
that both FMP and SFMP significantly outperform two competitors in all
settings, while SFMP is slightly inferior to FMP possibly due to the bias
introduced by the enforced sparsity.

Table 2: MAPEs and MSPEs of FMP, SFMP and two competing methods for Japanese
male and female mortality rates. All entries have been multiplied by $10$ for
formatting reasons.

|  | MAPE | MSPE
---|---|---|---
Male | $h$ | FMP | SFMP | GSY | UP | FMP | SFMP | GSY | UP
1 | 1.33 | 1.36 | 1.75 | 1.42 | 0.48 | 0.49 | 0.84 | 0.59
2 | 1.37 | 1.43 | 1.88 | 1.47 | 0.51 | 0.52 | 0.92 | 0.63
3 | 1.43 | 1.49 | 1.87 | 1.53 | 0.54 | 0.55 | 0.86 | 0.67
Female | $h$ | FMP | SFMP | GSY | UP | FMP | SFMP | GSY | UP
1 | 1.40 | 1.44 | 1.77 | 1.52 | 0.61 | 0.64 | 0.87 | 0.82
2 | 1.45 | 1.51 | 1.84 | 1.56 | 0.66 | 0.69 | 0.94 | 0.99
3 | 1.50 | 1.56 | 1.77 | 1.66 | 0.71 | 0.74 | 0.89 | 1.14

##  7 Discussion

We conclude the paper by discussing two possible extensions. First, our
estimation procedure is developed under the assumption that the idiosyncratic
component is a white noise sequence. It is also interesting to estimate the
factor model (1) following the framework in Fan et al. (2013), where the
idiosyncratic noise is allowed to exhibit serial correlations. Under the
uncorrelatedness assumption between ${\bf X}_{t}(\cdot)$ and
$\boldsymbol{\varepsilon}_{t}(\cdot),$ we can decompose the integrated
covariance function as

$\int\int\boldsymbol{\Sigma}^{(0)}_{yy}(u,v)\,\mathrm{d}u\mathrm{d}v={\bf
A}\Big{\\{}\int\int\boldsymbol{\Sigma}^{(0)}_{xx}(u,v)\,\mathrm{d}u\mathrm{d}v\Big{\\}}{\bf
A}^{\mathrm{\scriptscriptstyle
T}}+\int\int\boldsymbol{\Sigma}^{(0)}_{\varepsilon\varepsilon}(u,v)\,\mathrm{d}u\mathrm{d}v.$
(27)

By imposing some eigenvalue conditions on two terms in the right-side of (27)
similar to those in Fan et al. (2013), the common factors are asymptotically
identifiable and hence the factor loading space can be recovered by carrying
out an eigenanalysis for
$\int\int\widehat{\boldsymbol{\Sigma}}^{(0)}_{yy}(u,v)\,\mathrm{d}u\mathrm{d}v.$

Second, it is conceivable to integrate the high-dimensional time series
clustering approach Zhang et al. (2021) into our proposed factor model
framework (1) for the purpose of clustering a large number of functional time
series into different and unknown clusters such that curve dynamics within
each cluster share similar dynamic structure. Specifically, we can model
$p$-vector of functional time series ${\bf Y}_{t}(\cdot)$ based on a
functional latent factor structure, where the dynamic structures are
characterized by the latent common and cluster-specific functional factors
with each factor represented in the form of (1). Whereas the common factors
are strong, the cluster-specific factors are weak and only affect ${\bf
Y}_{t}(\cdot)$ in a specific cluster. The clustering is then based on factor
loadings on all weak factors.

## Appendix

##  Appendix A Additional regularity conditions

To study theoretical properties in an ultra-high-dimensional regime, we need
to rely on non-asymptotic error bounds on $\widehat{\Sigma}_{yy,ij}^{(k)}$’s
and will use the functional stability measure of $\\{{\bf
Y}_{t}(\cdot)\\}_{t\in\mathbb{Z}}$ Guo and Qiao (2021); Fang et al. (2021).

Before presenting regularity conditions, we first solidify some notation. We
denote the $p$-fold Cartesian product of $\mathbb{H}$ by
$\mathbb{H}^{p}=\mathbb{H}\times\cdots\times\mathbb{H}.$ For
${\mathbf{w}}=(w_{1},\dots,w_{g})^{{\mathrm{\scriptscriptstyle
T}}},{\mathbf{g}}=(g_{1},\dots,g_{p})^{{\mathrm{\scriptscriptstyle
T}}}\in\mathbb{H}^{p},$ we define
$\langle{\mathbf{w}},{\mathbf{g}}\rangle=\sum_{j=1}^{p}w_{j}g_{j}.$ For an
integral operator ${\bf K}:\mathbb{H}^{p}\rightarrow\mathbb{H}^{p}$ induced
from the kernel function ${\bf K}=(K_{ij})_{p\times p},$ ${\bf
K}({\mathbf{w}})(u)=\big{(}\sum_{j=1}^{p}\langle
K_{1j}(u,\cdot),w_{j}(\cdot)\rangle,\dots,\sum_{j=1}^{p}\langle
K_{pj}(u,\cdot),w_{j}(\cdot)\rangle\big{)}^{{\mathrm{\scriptscriptstyle T}}}$
for ${\mathbf{w}}\in\mathbb{H}^{p}.$ To simplify notation, we use ${\bf K}$ to
denote both the kernel function and the operator. We then define Hilbert space
valued sub-Gaussian random variable, which corresponds to an infinite-
dimensional analog of the sub-Gaussian random vector in $\mathbb{R}^{p}.$

######  Definition 1.

Let $W(\cdot)$ be a random variable in $\mathbb{H}$ and
$\Sigma_{W}:\mathbb{H}\rightarrow\mathbb{H}$ be the covariance operator of
$W(\cdot)$. Then $W(\cdot)$ is a sub-Gaussian process if there exists a
constant $\alpha>0$ such that for all $w\in\mathbb{H},$
$\mathbb{E}[\exp\\{\langle
w,W-\mathbb{E}(W)\rangle\\}]\leq\exp\\{2^{-1}\alpha^{2}\langle
w,\Sigma_{W}(w)\rangle\\}.$

######  Condition 7.

(i) $\\{{\bf Y}_{t}(\cdot)\\}_{t\in\mathbb{Z}}$ is a sequence of multivariate
functional linear processes with sub-Gaussian errors, i.e., ${\bf
Y}_{t}(\cdot)=\sum_{l=0}^{\infty}{\bf B}_{l}(\boldsymbol{\varepsilon}_{t-l})$
for any $t\in\mathbb{Z},$ where ${\bf B}_{l}=(B_{l,jk})_{p\times p}$ with each
$B_{l,jk}\in\mathbb{H}\otimes\mathbb{H},$
$\boldsymbol{\varepsilon}_{t}=(\varepsilon_{t1},\dots,\varepsilon_{tp})^{{\mathrm{\scriptscriptstyle
T}}}\in\mathbb{H}^{p}$ and the components in
$\\{\boldsymbol{\varepsilon}_{t}\\}_{t\in\mathbb{Z}}$ are independent mean-
zero sub-Gaussian processes satisfying Definition 1; (ii) The coefficient
functions satisfy $\sum_{l=0}^{\infty}\|B_{l,ij}\|_{{\cal S},\infty}=O(1);$
(iii) $\max_{j}\int_{{\cal
U}}\text{Cov}\\{\varepsilon_{tj}(u),\varepsilon_{tj}(u)\\}du=O(1).$

######  Condition 8.

For $\\{{\bf Y}_{t}(\cdot)\\}_{t\in\mathbb{Z}},$ its spectral density function
at frequency $\theta\in[-\pi,\pi],$
${\mathbf{f}}_{y,\theta}=(2\pi)^{-1}\sum_{k\in\mathbb{Z}}\boldsymbol{\Sigma}_{yy}^{(k)}e^{-ik\theta},$
exists, and its functional stability measure defined in (A.1) is finite,

${\cal
M}_{y}=2\pi\cdot\underset{\theta\in[-\pi,\pi],{\mathbf{w}}\in\mathbb{H}_{0}^{p}}{\text{esssup}}\frac{\langle{\mathbf{w}},{\mathbf{f}}_{y,\theta}({\mathbf{w}})\rangle}{\langle{\mathbf{w}},\boldsymbol{\Sigma}_{yy}^{(0)}({\mathbf{w}})\rangle}<\infty,$
(A.1)

where
$\mathbb{H}_{0}^{p}=\big{\\{}{\mathbf{w}}\in\mathbb{H}:\langle{\mathbf{w}},\boldsymbol{\Sigma}_{yy}^{(0)}({\mathbf{w}})\rangle\in(0,\infty)\big{\\}}.$

The multivariate functional linear process in Condition 7(i) generalizes the
multivariate linear process to the functional domain. Under this setting, the
functional stability measure ${\cal M}_{y}$ in (A.1) can be expressed based on
${\mathbf{f}}_{y,\theta}$ and $\boldsymbol{\Sigma}_{yy}^{(0)},$ both of which
have explicit expressions Fang et al. (2021). Under Conditions 7 and 8, the
rate in (22) can be implied from Lemma B3 in the Supplementary Material. In
general, we can relax Conditions 7(ii) and (iii) by allowing both
$\sum_{l=0}^{\infty}\|B_{l,ij}\|_{{\cal S},\infty}$ and $\max_{j}\int_{{\cal
U}}\text{Cov}\\{\varepsilon_{tj}(u),\varepsilon_{tj}(u)\\}du$ to grow at very
slow rates as $p$ diverges, then the rate in (22) will depend on these two
terms.

## References

  * (1)
  * Ahn and Horenstein (2013) Ahn, S. C. and Horenstein, A. R. (2013).  Eigenvalue ratio test for the number of factors, Econometrica 81: 1203–1227.
  * Aue et al. (2015) Aue, A., Norinho, D. D. and Hörmann, S. (2015).  On the prediction of stationary functional time series, Journal of the American Statistical Association 110: 378–392.
  * Bai and Ng (2002) Bai, J. and Ng, S. (2002).  Determining the number of factors in approximate factor models, Econometrica 70: 191–221.
  * Bathia et al. (2010) Bathia, N., Yao, Q. and Ziegelmann, F. (2010).  Identifying the finite dimensionality of curve time series, The Annals of Statistics 38: 3352–3386.
  * Bosq (2000) Bosq, D. (2000).  Linear Processes in Function Spaces, Springer, New York.
  * Cai and Liu (2011) Cai, T. and Liu, W. (2011).  Adaptive thresholding for sparse covariance matrix estimation, Journal of the American Statistical Association 106: 672–684.
  * Carvalho et al. (2008) Carvalho, C. M., Chang, J., Lucas, J. E., Nevins, J. R., Wang, Q. and West, M. (2008).  High-dimensional sparse factor modeling: applications in gene expression genomics, Journal of the American Statistical Association 103: 1438–1456.
  * Chamberlain and Rothschild (1983) Chamberlain, G. and Rothschild, M. (1983).  Arbitrage, factor structure, and mean-variance analysis on large asset markets, Econometrica 51: 1281–1304.
  * Chang et al. (2018) Chang, J., Guo, B. and Yao, Q. (2018).  Principal component analysis for second-order stationary vector time series, The Annals of Statistics 46: 2094–2124.
  * Chen, Guo and Qiao (2020) Chen, C., Guo, S. and Qiao, X. (2020).  Functional linear regression: dependence and error contamination, Journal of Business & Economic Statistics, in press .
  * Chen, Tsay and Chen (2020) Chen, E. Y., Tsay, R. S. and Chen, R. (2020).  Constrained factor models for high-dimensional matrix-variate time series, Journal of the American Statistical Association 115: 775–793.
  * Chen et al. (2021) Chen, R., Yang, D. and Zhang, C.-H. (2021).  Factor models for high-dimensional tensor time series, Journal of the American Statistical Association, in press .
  * d’Aspremont et al. (2008) d’Aspremont, A., Bach, F. and El Ghaoui, L. (2008).  Optimal solutions for sparse principal component analysis., Journal of Machine Learning Research 9: 1269–1294.
  * d’Aspremont et al. (2007) d’Aspremont, A., El Ghaoui, L., Jordan, M. I. and Lanckriet, G. R. (2007).  A direct formulation for sparse pca using semidefinite programming, SIAM review 49(3): 434–448.
  * Fan et al. (2020) Fan, J., Guo, J. and Zheng, S. (2020).  Estimating number of factors by adjusted eigenvalues thresholding, Journal of the American Statistical Association, in press .
  * Fan et al. (2013) Fan, J., Liao, Y. and Mincheva, M. (2013).  Large covariance estimation by thresholding principal orthogonal complements, Journal of the Royal Statistical Society: Series B 75: 603–680.
  * Fan et al. (2016) Fan, J., Liao, Y. and Wang, W. (2016).  Projected principal component analysis in factor models, The Annals of Statistics 44: 219 – 254.
  * Fang et al. (2021) Fang, Q., Guo, S. and Qiao, X. (2021).  Finite sample theory for high-dimensional functional/scalar time series with applications, Electronic Journal of Statistics, in press .
  * Forni et al. (2000) Forni, M., Hallin, M., Lippi, M. and Reichlin, L. (2000).  The generalized dynamic-factor model: Identification and estimation, Review of Economics and statistics 82: 540–554.
  * Gao et al. (2019) Gao, Y., Shang, H. L. and Yang, Y. (2019).  High-dimensional functional time series forecasting: An application to age-specific mortality rates, Journal of Multivariate Analysis 170: 232–243.
  * Gataric et al. (2020) Gataric, M., Wang, T. and Samworth, R. J. (2020).  Sparse principal component analysis via axis-aligned random projections, Journal of the Royal Statistical Society: Series B 82: 329–359.
  * Guo and Qiao (2021) Guo, S. and Qiao, X. (2021).  On consistency and sparsity for high-dimensional functional time series with application to autoregressions, arXiv:2003.11462v2 .
  * Hays et al. (2012) Hays, S., Shen, H. and Huang, J. Z. (2012).  Functional dynamic factor models with application to yield curve forecasting, The Annals of Applied Statistics 6: 870–894.
  * Hörmann et al. (2015) Hörmann, S., Kidzinski, L. and Hallin, M. (2015).  Dynamic functional principal components, Journal of the Royal Statistical Society: Series B 77: 319–348.
  * Hörmann and Kokoszka (2010) Hörmann, S. and Kokoszka, P. (2010).  Weakly dependent functional data, The Annals of Statistics 38: 1845–1884.
  * Klepsch et al. (2017) Klepsch, J., Kluppelberg, C. and Wei, T. (2017).  Prediction of functional ARMA processes with an application to traffic data, Econometrics and Statistics 1: 128–149.
  * Kokoszka et al. (2015) Kokoszka, P., Miao, H. and Zhang, X. (2015).  Functional dynamic factor model for intraday price curves, Journal of Financial Econometrics 13: 456–477.
  * Lam and Yao (2012) Lam, C. and Yao, Q. (2012).  Factor modelling for high-dimensional time series: inference for the number of factors, The Annals of Statistics 40: 694–726.
  * Lam et al. (2011) Lam, C., Yao, Q. and Bathia, N. (2011).  Estimation of latent factors for high-dimensional time series, Biometrika 98: 901–918.
  * Li et al. (2020) Li, D., Robinson, P. M. and Shang, H. L. (2020).  Long-range dependent curve time series, Journal of the American Statistical Association 115: 957–971.
  * Mackey (2009) Mackey, L. (2009).  Deflation methods for sparse PCA, Advances in Neural Information Processing Systems 21: 1017–1024.
  * Moghaddam et al. (2006) Moghaddam, B., Weiss, Y. and Avidan, S. (2006).  Generalized spectral bounds for sparse LDA, Proceedings of the 23rd international conference on Machine learning, pp. 641–648.
  * Onatski (2010) Onatski, A. (2010).  Determining the number of factors from empirical distribution of eigenvalues, The Review of Economics and Statistics 92: 1004–1016.
  * Pan and Yao (2008) Pan, J. and Yao, Q. (2008).  Modelling multiple time series via common factors, Biometrika 95: 365–379.
  * Panaretos and Tavakoli (2013) Panaretos, V. and Tavakoli, S. (2013).  Fourier analysis of stationary time series in function space, The Annals of Statistics 41: 568–603.
  * Pena and Box (1987) Pena, D. and Box, G. E. (1987).  Identifying a simplifying structure in time series, Journal of the American statistical Association 82: 836–843.
  * Qiao et al. (2020) Qiao, X., Qian, C., James, G. M. and Guo, S. (2020).  Doubly functional graphical models in high dimensions, Biometrika 107: 415–431.
  * Stock and Watson (2012) Stock, J. and Watson, M. (2012).  Dynamic Factor Models, Oxford University Press, Oxford.
  * Tavakoli et al. (2021) Tavakoli, S., Nisol, G. and Hallin, M. (2021).  Factor models for high-dimensional functional time series, arXiv:1905.10325 .
  * Vu and Lei (2013) Vu, V. Q. and Lei, J. (2013).  Minimax sparse principal subspace estimation in high dimensions, The Annals of Statistics 41: 2905–2947.
  * Wang et al. (2019) Wang, D., Liu, X. and Chen, R. (2019).  Factor models for matrix-valued high-dimensional time series, Journal of Econometrics 208: 231–248.
  * Zhang et al. (2021) Zhang, B., Pan, G., Yao, Q. and Zhou, W. (2021).  Factor modelling for clustering high-dimensional time series, arXiv:2101.01908 .
  * Zou and Xue (2018) Zou, H. and Xue, L. (2018).  A selective overview of sparse principal component analysis, Proceedings of the IEEE 106: 1311–1320.

Supplementary material to “Factor Modelling for High-Dimensional Functional
Time Series”

Shaojun Guo, Xinghao Qiao and Qingsong Wang

This supplementary material contains all technical proofs in Section B and
additional empirical results in Section C.

##  Appendix B Technical proofs

In addition to the notation defined in the main paper, we summarize here more
notation to be used throughout the supplementary material. For a matrix ${\bf
A}=(A_{ij})_{p\times q},$ we denote its Frobenius norm, matrix $\ell_{1}$ norm
and matrix $\ell_{\infty}$ norm by $\|{\bf
A}\|_{\text{F}}=(\sum_{i=1}^{p}\sum_{j=1}^{q}A_{ij}^{2})^{1/2},$ $\|{\bf
A}\|_{1}=\max_{1\leq j\leq q}\sum_{i=1}^{p}|A_{ij}|$ and $\|{\bf
A}\|_{\infty}=\max_{1\leq i\leq p}\sum_{j=1}^{q}|A_{ij}|,$ respectively. For a
matrix of bivariate functions ${\bf K}=(K_{ij}(\cdot,\cdot))_{p\times q}$ with
each $K_{ij}\in\mathbb{H}\otimes\mathbb{H}$, we define functional versions of
the Frobenius norm and matrix $\ell_{1}$-norm by $\|{\bf K}\|_{{\cal
S},\text{F}}=(\sum_{i=1}^{p}\sum_{j=1}^{q}\|K_{ij}\|_{{\cal S}}^{2})^{1/2}$
and $\|{\bf K}\|_{{\cal S},1}=\max_{1\leq j\leq
q}\sum_{i=1}^{p}\|K_{ij}\|_{{\cal S}}$, respectively.

###  B.1 Proof of Theorem 1

Recall that

$\widehat{\bf
M}=\sum_{k=1}^{k_{0}}\iint\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)\widehat{\bf
W}(u,v)\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle
T}}\,{\rm d}u{\rm d}v,$ ${\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}=\sum_{k=1}^{k_{0}}\iint\boldsymbol{\Sigma}^{(k)}_{yy}(u,v)\widehat{\bf
W}(u,v)\boldsymbol{\Sigma}^{(k)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle
T}}\,{\rm d}u{\rm d}v,$

where $\widehat{\bf W}(u,v)={\bf Q}\big{\\{}{\bf
Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf Q}\big{\\}}^{-1}{\bf
Q}^{\mathrm{\scriptscriptstyle T}}.$ Noting that $\widehat{\bf W}(u,v)$ is
scale-invariant with respect to ${\bf Q},$ we assume here that the entries of
${\bf Q}=(Q_{ij})_{p\times q}$ are independent $N(0,p^{-1})$ instead of
$N(0,1)$ to ensure that the entries of ${\bf Q}^{\mathrm{\scriptscriptstyle
T}}{\widehat{\boldsymbol{\Sigma}}}_{yy}^{(0)}(v,v){\bf Q}$ are $O_{P}(1)$
instead of $O_{P}(p^{2}).$ Provided that the maximum lag $k_{0}$ is fixed,
without loss of generality we only need to consider the case $k_{0}=1.$ We
organize our proof in the following four steps.

Step 1. We will show that ${\bf Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf Q}$ is positive definite
with high probability uniformly over ${\cal U},$ i.e., there exists some
positive constant $C>0$ such that, with probability tending to one,

$\inf_{v\in{\cal U}}\lambda_{\min}\big{\\{}{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf Q}\big{\\}}>C.$ (S.1)

Observe that ${\bf Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf Q}={\bf
Q}^{\mathrm{\scriptscriptstyle T}}\boldsymbol{\Sigma}_{yy}^{(0)}(v,v){\bf
Q}+{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\left\\{\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v)-\boldsymbol{\Sigma}_{yy}^{(0)}(v,v)\right\\}{\bf
Q}.$ Then,

$\left|\lambda_{k}\Big{\\{}{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf
Q}\Big{\\}}-\lambda_{k}\Big{\\{}{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\boldsymbol{\Sigma}_{yy}^{(0)}(v,v){\bf Q}\Big{\\}}\right|\leq\left\|{\bf
Q}^{\mathrm{\scriptscriptstyle
T}}\left\\{\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v)-\boldsymbol{\Sigma}_{yy}^{(0)}(v,v)\right\\}{\bf
Q}\right\|_{\text{F}},$ (S.2)

where $\lambda_{k}({\bf A})$ is the $k$-th largest eigenvalue of a symmetric
matrix ${\bf A}.$ By the Hanson–Wright inequality and union bound of
probability, we obtain that there exists some constant $c>0$ such that for any
matrix ${\bf A}\in\mathbb{R}^{p\times p}$ and $t>0,$

$\displaystyle\mathbb{P}\left\\{\frac{\big{\|}{\bf
Q}^{\mathrm{\scriptscriptstyle T}}{\bf A}{\bf Q}-{\mathbb{E}}{\bf
Q}^{\mathrm{\scriptscriptstyle T}}{\bf A}{\bf Q}\big{\|}_{\text{F}}}{\|{\bf
A}\|_{\text{F}}}\geq tqp^{-1}\right\\}\leq 2q^{2}\exp\\{-c\min(t^{2},t)\\}.$

It follows from the concentration inequality above with ${\bf
A}=\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v)-{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v),$
the independence between ${\bf Q}$ and
$\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v),$ the Lipschitz-continuity of
$\boldsymbol{\Sigma}_{yy}^{(0)}(v,v)$ and the uniform convergence rate in
Condition 5 that

$\displaystyle\sup_{v\in{\cal U}}\left\|{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\left\\{\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v)-\boldsymbol{\Sigma}_{yy}^{(0)}(v,v)\right\\}{\bf
Q}\right\|_{\text{F}}$ $\displaystyle\leq$ $\displaystyle
O_{P}\Big{\\{}qp^{-1}\\{\log(qn)\\}^{1/2}\sup_{v\in{\cal
U}}\|\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v)-\boldsymbol{\Sigma}_{yy}^{(0)}(v,v)\|_{\text{F}}\Big{\\}}$
$\displaystyle=$ $\displaystyle O_{P}\big{\\{}q\log(qn)n^{-1/2}\big{\\}}.$

Combining the above result with (S.2) and the fact that

$\inf_{v\in{\cal U}}\lambda_{q}\\{{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\boldsymbol{\Sigma}_{yy}^{(0)}(v,v){\bf Q}\\}\geq\inf_{v\in{\cal
U}}\lambda_{q}\\{{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{(0)}(v,v){\bf Q}\\}>c>0$

under Condition 2 implies that (S.1) holds with high probability.

Step 2. The bounds of $r$ largest eigenvalues of  $\textstyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule height=0.0pt,width=9.16669pt}$
$\textstyle{\bf M}$ (i.e., $\nu_{1}\geq\cdots\geq\nu_{r}$) satisfy

$\nu_{j}\asymp p^{2-2\delta},~{}~{}j=1,\cdots,r,$ (S.3)

with probability tending to one.

First, note that

$\nu_{r}\leq\nu_{1}\leq\|{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}\|\leq\iint\|\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\|^{2}\|\widehat{\bf
W}(u,v)\|\,{\rm d}u{\rm d}v.$ (S.4)

Applying the fact $\|{\bf Q}\|\asymp 1+o_{P}(1)$ and (S.1) from Step 1, we
have that

$\sup_{u,v\in{\cal U}}\|\widehat{\bf W}(u,v)\|\leq\|{\bf
Q}\|^{2}\big{[}\inf_{v\in{\cal U}}\lambda_{\min}\\{{\bf
Q}^{T}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf
Q}\\}\big{]}^{-1}\lesssim 1+o_{P}(1).$ (S.5)

It follows from Conditions 3–4, the decomposition (3) and Lemma B2 that

$\displaystyle\iint\big{\|}\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\big{\|}^{2}\,{\rm
d}u{\rm d}v$ $\displaystyle\leq$ $\displaystyle\iint\left(\big{\|}{\bf
A}\boldsymbol{\Sigma}_{xx}^{(1)}(u,v){\bf A}^{{\mathrm{\scriptscriptstyle
T}}}\big{\|}+\big{\|}{\bf
A}\boldsymbol{\Sigma}_{x\varepsilon}^{(1)}(u,v)\big{\|}\right)^{2}\,{\rm
d}u{\rm d}v$ (S.6) $\displaystyle\leq$ $\displaystyle 2\Big{\\{}\big{\|}{\bf
A}\big{\|}^{4}\iint\big{\|}\boldsymbol{\Sigma}_{xx}^{(1)}(u,v)\big{\|}^{2}\,{\rm
d}u{\rm d}v+\big{\|}{\bf
A}\big{\|}^{2}\iint\big{\|}\boldsymbol{\Sigma}_{x\varepsilon}^{(1)}(u,v)\big{\|}^{2}{\rm
d}u{\rm d}v\Big{\\}}$ $\displaystyle=$ $\displaystyle
O\Big{\\{}p^{2-2\delta}+p^{1-\delta}\iint\big{\|}\boldsymbol{\Sigma}_{x\varepsilon}^{(1)}(u,v)\big{\|}^{2}\,{\rm
d}u{\rm d}v\Big{\\}},$

where, by Condition 4,

$\iint\big{\|}\boldsymbol{\Sigma}_{x\varepsilon}^{(1)}(u,v)\big{\|}^{2}\,{\rm
d}u{\rm d}v\leq\big{\|}\boldsymbol{\Sigma}_{x\varepsilon}^{(1)}\big{\|}_{{\cal
S},\text{F}}^{2}=O\Big{\\{}\max\limits_{1\leq i\leq
r}\sum_{j=1}^{p}\big{\|}\Sigma_{x\epsilon,ij}^{(1)}\big{\|}^{2}_{{\cal
S}}\Big{\\}}\leq
O\Big{\\{}\big{\|}\boldsymbol{\Sigma}_{x\varepsilon}^{(1)}\big{\|}_{{\cal
S},\infty}^{2}\Big{\\}}=o(p^{1-\delta}).$ (S.7)

Combining (S.4)–(S.7) yields that $\nu_{r}\lesssim
p^{2-2\delta}\\{1+o_{P}(1)\\}.$

We next give a lower bound on $\nu_{r}.$ By Conditions 1, 3, Lemma 3 of Wang
and Xi (1997) and (S.1) from Step 1, we first obtain that

$\displaystyle I_{1}$
$\displaystyle=\lambda_{r}\left\\{\iint\boldsymbol{\Sigma}^{(1)}_{xx}(u,v){\bf
A}^{{\mathrm{\scriptscriptstyle T}}}\widehat{\bf W}(v,v){\bf
A}\boldsymbol{\Sigma}^{(1)}_{xx}(u,v)^{\mathrm{\scriptscriptstyle T}}\,{\rm
d}u{\rm d}v\right\\}$
$\displaystyle\geq\lambda_{r}\left\\{\iint\boldsymbol{\Sigma}^{(1)}_{xx}(u,v)\boldsymbol{\Sigma}^{(1)}_{xx}(u,v)^{\mathrm{\scriptscriptstyle
T}}\,{\rm d}u{\rm d}v\right\\}\inf_{u,v\in{\cal U}}\lambda_{r}\left\\{{\bf
A}^{{\mathrm{\scriptscriptstyle T}}}\widehat{\bf W}(u,v){\bf A}\right\\}$
$\displaystyle\gtrsim\lambda_{r}({\bf A}^{\mathrm{\scriptscriptstyle T}}{\bf
Q}{\bf Q}^{\mathrm{\scriptscriptstyle T}}{\bf A})\cdot\inf_{v\in{\cal
U}}\lambda_{q-r+1}\left(\big{\\{}{\bf
Q}^{T}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf
Q}\big{\\}}^{-1}\right)$ $\displaystyle\gtrsim p^{1-\delta}\\{1+o_{P}(1)\\}.$

It follows from (S.5), (S.7), Lemma B2 and Condition 3 that

$\displaystyle I_{2}$
$\displaystyle=\left\|\iint\boldsymbol{\Sigma}^{(1)}_{x\varepsilon}(u,v)\widehat{\bf
W}(u,v)\boldsymbol{\Sigma}^{(1)}_{xx}(u,v){\bf A}^{{\mathrm{\scriptscriptstyle
T}}}\,{\rm d}u{\rm d}v\right\|$
$\displaystyle\leq\left\\{\iint\|\boldsymbol{\Sigma}^{(1)}_{x\varepsilon}(u,v))\|^{2}\,{\rm
d}u{\rm d}v\iint\|\boldsymbol{\Sigma}^{(1)}_{xx}(u,v)\|^{2}\,{\rm d}u{\rm
d}v\right\\}^{1/2}\sup_{u,v\in{\cal U}}\|\widehat{\bf W}(u,v)\|\|{\bf
A}\|=o_{P}(p^{1-\delta}),$

By Weyl’s inequality, the decomposition (3) and Condition 3, $\nu_{r}$ can be
bounded from below by

$\displaystyle\nu_{r}$
$\displaystyle=\lambda_{r}\Big{\\{}\iint\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\widehat{\bf
W}(u,v)\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)^{{\mathrm{\scriptscriptstyle
T}}}\,{\rm d}u{\rm d}v\Big{\\}}$ $\displaystyle\geq\lambda_{r}({\bf
A}^{\mathrm{\scriptscriptstyle T}}{\bf
A})\cdot\lambda_{r}\Big{[}\iint\big{\\{}\boldsymbol{\Sigma}^{(1)}_{xx}(u,v){\bf
A}^{{\mathrm{\scriptscriptstyle
T}}}+\boldsymbol{\Sigma}^{(1)}_{x\varepsilon}(u,v)\big{\\}}\widehat{\bf
W}(u,v)\big{\\{}\boldsymbol{\Sigma}^{(1)}_{xx}(u,v){\bf
A}^{{\mathrm{\scriptscriptstyle
T}}}+\boldsymbol{\Sigma}^{(1)}_{x\varepsilon}(u,v)\big{\\}}^{{\mathrm{\scriptscriptstyle
T}}}\,{\rm d}u{\rm d}v\Big{]}$ $\displaystyle\geq\|{\bf
A}\|_{\min}^{2}\cdot\lambda_{r}\Big{\\{}\iint\boldsymbol{\Sigma}^{(1)}_{xx}(u,v){\bf
A}^{{\mathrm{\scriptscriptstyle T}}}\widehat{\bf W}(u,v){\bf
A}\boldsymbol{\Sigma}^{(1)}_{xx}(u,v)^{\mathrm{\scriptscriptstyle T}}\,{\rm
d}u{\rm d}v\Big{\\}}$ $\displaystyle~{}~{}~{}~{}+2\|{\bf
A}\|_{\min}^{2}\cdot\lambda_{r}\Big{\\{}\iint\boldsymbol{\Sigma}^{(1)}_{x\varepsilon}(u,v)\widehat{\bf
W}(u,v)\boldsymbol{\Sigma}^{(1)}_{xx}(u,v){\bf A}^{{\mathrm{\scriptscriptstyle
T}}}\,{\rm d}u{\rm d}v\Big{\\}}$ $\displaystyle\geq\|{\bf
A}\|_{\min}^{2}\cdot(I_{1}-2I_{2})\gtrsim p^{2-2\delta}\\{1+o_{P}(1)\\}.$

Hence we obtain that $\nu_{r}\asymp p^{2-2\delta}$ with high probability.

Step 3. We will show that

$\|\widehat{\bf M}-{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}\|=O_{P}(p^{2-\delta}n^{-1/2}).$ (S.8)

Observe that

$\displaystyle\widehat{\boldsymbol{\Sigma}}^{(1)}_{yy}(u,v)\widehat{\bf
W}(u,v)\widehat{\boldsymbol{\Sigma}}^{(1)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle
T}}-\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\widehat{\bf
W}(u,v)\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle T}}$
$\displaystyle=$
$\displaystyle\big{\\{}\widehat{\boldsymbol{\Sigma}}^{(1)}_{yy}(u,v)-\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\big{\\}}\widehat{\bf
W}(u,v)\big{\\{}\widehat{\boldsymbol{\Sigma}}^{(1)}_{yy}(u,v)-\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\big{\\}}^{\mathrm{\scriptscriptstyle
T}}$
$\displaystyle+\big{\\{}\widehat{\boldsymbol{\Sigma}}^{(1)}_{yy}(u,v)-\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\big{\\}}\widehat{\bf
W}(u,v)\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)^{\mathrm{\scriptscriptstyle T}}$
$\displaystyle+\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\widehat{\bf
W}(u,v)\big{\\{}\widehat{\boldsymbol{\Sigma}}^{(1)}_{yy}(u,v)-\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\big{\\}}^{\mathrm{\scriptscriptstyle
T}}.$

By $\|{\bf Q}\|\asymp 1+o_{P}(1)$, Condition 5 and Markov’s inequality, we
have that

$\displaystyle\iint\left\|\left\\{\widehat{\boldsymbol{\Sigma}}^{(1)}_{yy}(u,v)-\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\right\\}{\bf
Q}\right\|^{2}{\rm d}u{\rm d}v$
$\displaystyle\leq\sum_{i=1}^{p}\sum_{j=1}^{p}\left\|\widehat{\Sigma}^{(1)}_{yy,ij}-\Sigma^{(1)}_{yy,ij}\right\|_{{\cal
S}}^{2}\|{\bf Q}\|^{2}=O_{P}(p^{2}n^{-1}).$

Furthermore, (S.6) and (S.7) imply that

$\displaystyle\iint\left\|\boldsymbol{\Sigma}^{(1)}_{yy}(u,v){\bf
Q}\right\|^{2}{\rm d}u{\rm
d}v\leq\iint\left\|\boldsymbol{\Sigma}^{(1)}_{yy}(u,v)\right\|^{2}\,{\rm
d}u{\rm d}v\,\|{\bf Q}\|^{2}=p^{2-2\delta}\\{1+o_{P}(1)\\}.$

It follows from the above decomposition, two upper bound results,
Cauchy–Schwartz inequality and $\sup_{v\in{\cal
U}}\lambda_{\max}\big{(}\\{{\bf Q}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\Sigma}}_{yy}^{(0)}(v,v){\bf
Q}\\}^{-1}\big{)}=O_{P}(1)$ from Step 1 that (S.8) holds.

Step 4. Combining (S.3) and (S.8) from Steps 2 and 3 with Lemma 3 and the
proof of Theorem 1 in Lam et al. (2011), we have that there exists
$\widetilde{{\bf K}}$ and an orthogonal matrix ${\bf U}$ such that
$\widehat{{\bf K}}=\widetilde{{\bf K}}{\bf U}$ and

$\|\widetilde{\bf K}-{\bf K}\|=O_{P}(p^{\delta}n^{-1/2}).$

It then follows from the orthogonality of ${\bf K}$ that

$\displaystyle\mathcal{D}^{2}\Big{(}\mathcal{C}({\bf
K}),\mathcal{C}(\widehat{\bf K})\Big{)}$
$\displaystyle=r^{-1}\mbox{tr}\Big{(}{\bf I}-{\bf K}{\bf
K}^{{\mathrm{\scriptscriptstyle T}}}\widehat{\bf K}\widehat{{\bf
K}}^{{\mathrm{\scriptscriptstyle T}}}\Big{)}=r^{-1}\mbox{tr}({\bf
K}^{{\mathrm{\scriptscriptstyle T}}}({\bf I}-\widehat{\bf K}\widehat{{\bf
K}}^{{\mathrm{\scriptscriptstyle T}}}){\bf K})$
$\displaystyle=r^{-1}\Big{\\{}\mbox{tr}({\bf K}^{{\mathrm{\scriptscriptstyle
T}}}({\bf I}-\widehat{\bf K}\widehat{{\bf K}}^{{\mathrm{\scriptscriptstyle
T}}}){\bf K})-\mbox{tr}({\bf K}^{{\mathrm{\scriptscriptstyle T}}}({\bf I}-{\bf
K}{\bf K}^{{\mathrm{\scriptscriptstyle T}}}){\bf K})\Big{\\}}$
$\displaystyle=r^{-1}\mbox{tr}\Big{(}{\bf K}^{{\mathrm{\scriptscriptstyle
T}}}({\bf K}{\bf K}^{{\mathrm{\scriptscriptstyle T}}}-\widehat{\bf
K}\widehat{{\bf K}}^{{\mathrm{\scriptscriptstyle T}}}){\bf
K}\Big{)}\leq\Big{\|}{\bf K}^{{\mathrm{\scriptscriptstyle T}}}({\bf K}{\bf
K}^{{\mathrm{\scriptscriptstyle T}}}-\widetilde{\bf K}\widetilde{{\bf
K}}^{{\mathrm{\scriptscriptstyle T}}}){\bf K}\Big{\|},$

where the last term is further bounded by

$\displaystyle\|-{\bf K}^{{\mathrm{\scriptscriptstyle T}}}({\bf
K}-\widetilde{{\bf K}})({\bf K}-\widetilde{{\bf
K}})^{{\mathrm{\scriptscriptstyle T}}}{\bf K}+({\bf K}-\widetilde{{\bf
K}})^{{\mathrm{\scriptscriptstyle T}}}({\bf K}-\widetilde{{\bf K}})\|\leq
2\|\widetilde{{\bf K}}-{\bf K}\|^{2}=O_{P}(p^{2\delta}n^{-1}).$

The result in Theorem 1 follows immediately, which completes the proof.
$\square$

###  B.2 Proof of Theorem 2

Denote $(\nu_{j},\boldsymbol{\gamma}_{j})$ and
$(\hat{\nu}_{j},\widehat{\boldsymbol{\gamma}}_{j})$ by the eigenpairs of
$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$ $\textstyle{\bf M}$ and $\widehat{\bf M}$,
respectively. For simplicity we assume
$\boldsymbol{\gamma}_{j}^{\mathrm{\scriptscriptstyle
T}}\widehat{\boldsymbol{\gamma}}_{j}\geq 0$ here. First note that for
$j=1,2,\dots,r,$

$\displaystyle\hat{\nu}_{j}-\nu_{j}=\widehat{\boldsymbol{\gamma}}_{j}^{\mathrm{\scriptscriptstyle
T}}\widehat{\bf
M}\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j}^{\mathrm{\scriptscriptstyle
T}}{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}\boldsymbol{\gamma}_{j}=I_{1}+I_{2}+I_{3}+I_{4}+I_{5},$

where

$\displaystyle I_{1}$
$\displaystyle=(\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j})^{\mathrm{\scriptscriptstyle
T}}(\widehat{\bf M}-{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}})\widehat{\boldsymbol{\gamma}}_{j},\qquad
I_{2}=(\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j})^{\mathrm{\scriptscriptstyle
T}}{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}(\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j})$
$\displaystyle I_{3}$
$\displaystyle=(\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j})^{\mathrm{\scriptscriptstyle
T}}{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}\boldsymbol{\gamma}_{j},\qquad
I_{4}=\boldsymbol{\gamma}_{j}^{\mathrm{\scriptscriptstyle T}}(\widehat{\bf
M}-{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}})\widehat{\boldsymbol{\gamma}}_{j},\qquad
I_{5}=\boldsymbol{\gamma}_{j}^{\mathrm{\scriptscriptstyle
T}}{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}(\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j}).$

Since
$\|\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j}\|\leq\|\widetilde{{\bf
K}}-{\bf K}\|=O_{P}(p^{\delta}n^{-1/2}),$ $\|\widehat{\bf
M}-{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}\|=O_{P}(p^{2-\delta}n^{-1/2})$ and
$\|{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}\|=O_{P}(p^{2-2\delta})$, we have that

$|I_{1}|+|I_{2}|=O_{P}(p^{2}n^{-1}),\qquad|I_{3}|+|I_{4}|+|I_{5}|=O_{P}(p^{2-\delta}n^{-1/2}),$

which means that $|\hat{\nu}_{j}-\nu_{j}|=O_{P}(p^{2-\delta}n^{-1/2})$ under
$p^{\delta}n^{-1/2}=o(1).$

To prove the result for $j=r+1,\dots,p$, we introduce

$\widebar{\bf
M}=\sum_{k=1}^{k_{0}}\iint\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}(u,v)\widehat{\bf
W}(v,v)\boldsymbol{\Sigma}^{(k)}_{yy}(u,v)^{{\mathrm{\scriptscriptstyle
T}}}\,\mathrm{d}u\mathrm{d}v.$

For each $j=r+1,\ldots,p,$ consider the following decomposition

$\displaystyle\hat{\nu}_{j}=\widehat{\boldsymbol{\gamma}}_{j}^{\mathrm{\scriptscriptstyle
T}}\widehat{\bf M}\widehat{\boldsymbol{\gamma}}_{j}=I_{6}+I_{7}+I_{8},$

where

$\displaystyle I_{6}$
$\displaystyle=\widehat{\boldsymbol{\gamma}}_{j}^{\mathrm{\scriptscriptstyle
T}}(\widehat{\bf M}-\widebar{{\bf M}}-\widebar{{\bf
M}}^{\mathrm{\scriptscriptstyle T}}+{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}})\widehat{\boldsymbol{\gamma}}_{j},\qquad
I_{7}=2\widehat{\boldsymbol{\gamma}}_{j}^{\mathrm{\scriptscriptstyle
T}}(\widebar{{\bf M}}-{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}})(\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j}),$
$\displaystyle I_{8}$
$\displaystyle=(\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j})^{\mathrm{\scriptscriptstyle
T}}{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}(\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j}).$

Note that, under Condition 5,

$|I_{6}|\leq\|\widehat{\bf M}-\widebar{{\bf M}}-\widebar{{\bf
M}}^{\mathrm{\scriptscriptstyle T}}+{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}\|\lesssim\sum_{k=1}^{k_{0}}\|\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}-\boldsymbol{\Sigma}^{(k)}_{yy}\|_{{\cal
S}}^{2}=O_{P}(p^{2}n^{-1}).$

By the similar arguments to prove Theorem 1, we can show that $\|\widebar{\bf
M}-{\mathchoice{{\ooalign{\hbox{\raise
7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\displaystyle\widehat{\vrule
width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\displaystyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 7.09259pt\hbox{\scalebox{1.0}[-1.0]{\lower
7.09259pt\hbox{$\textstyle\widehat{\vrule width=0.0pt,height=6.83331pt\vrule
height=0.0pt,width=9.16669pt}$}}}}\cr\hbox{$\textstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 6.40926pt\hbox{\scalebox{1.0}[-1.0]{\lower
6.40926pt\hbox{$\scriptstyle\widehat{\vrule width=0.0pt,height=4.78333pt\vrule
height=0.0pt,width=6.41667pt}$}}}}\cr\hbox{$\scriptstyle{\bf
M}$}}}}{{\ooalign{\hbox{\raise 5.9537pt\hbox{\scalebox{1.0}[-1.0]{\lower
5.9537pt\hbox{$\scriptscriptstyle\widehat{\vrule
width=0.0pt,height=3.41666pt\vrule
height=0.0pt,width=4.58334pt}$}}}}\cr\hbox{$\scriptscriptstyle{\bf
M}$}}}}}\|=O_{P}(p^{2-\delta}n^{-1/2})$ and
$\|\widehat{\boldsymbol{\gamma}}_{j}-\boldsymbol{\gamma}_{j}\|=O_{P}(p^{\delta}n^{-1/2}),$
which together imply that $|I_{7}|+|I_{8}|=O_{P}(p^{2}n^{-1}).$ Combining the
above results yields that $\hat{\nu}_{j}=O_{p}(p^{2}n^{-1})$ for
$j=r+1,\dots,p,$ which completes the proof of the first part.

We next turn to prove the second part. It follows from $\nu_{j}\asymp
p^{2-2\delta}$ for $j=1,\ldots,r$ in (S.3) and $p^{\delta}n^{-1/2}=o(1)$ that

$\frac{|\hat{\nu}_{j}-\nu_{j}|}{\nu_{j}}=O_{P}(p^{2-\delta}n^{-1/2}p^{2\delta-2})=O_{P}(p^{\delta}n^{-1/2})=o_{P}(1),~{}~{}j=1,\dots,r.$

As a result, for $j=1,\ldots,r-1,$

$\frac{\hat{\nu}_{j+1}}{\hat{\nu}_{j}}=\frac{\nu_{j+1}}{\nu_{j}}\frac{1+\nu_{j+1}^{-1}(\hat{\nu}_{j+1}-\nu_{j+1})}{1+\nu_{j}^{-1}(\hat{\nu}_{j}-\nu_{j})}=\frac{\nu_{j+1}}{\nu_{j}}\frac{1+o_{P}(1)}{1+o_{P}(1)}\asymp
1,$

and

$\frac{\hat{\nu}_{r+1}}{\hat{\nu}_{r}}=\frac{1}{\nu_{r}}\frac{\hat{\nu}_{r+1}}{1+\nu_{r}^{-1}(\hat{\nu}_{r}-\nu_{r})}=O_{P}(p^{2\delta-2}p^{2}n^{-1})=O_{P}(p^{2\delta}n^{-1}).$

We complete the proof of the second part. $\square$

###  B.3 Proof of Lemma 1

Denote by ${\bf A}=(A_{ij})_{p\times r}$. By the decomposition (3) and the
fact $(x+y)^{\tau}\leq x^{\tau}+y^{\tau}$ for $x,y\geq 0$ and $\tau\in[0,1),$
we have that

$\displaystyle\big{\|}\Sigma^{(k)}_{yy,ij}\big{\|}^{\tau}_{{\cal S}}$
$\displaystyle\leq\Big{\|}\sum_{1\leq m,m^{\prime}\leq
r}A_{im}\Sigma^{(k)}_{xx,mm^{\prime}}A_{jm^{\prime}}\Big{\|}_{{\cal
S}}^{\tau}+\Big{\|}\sum_{m=1}^{r}A_{im}\Sigma^{(k)}_{x\epsilon,mj}\Big{\|}_{{\cal
S}}^{\tau}$ $\displaystyle\leq\sum_{1\leq m,m^{\prime}\leq
r}\big{\|}A_{im}\Sigma^{(k)}_{xx,mm^{\prime}}A_{jm^{\prime}}\big{\|}_{{\cal
S}}^{\tau}+\sum_{m=1}^{r}\big{\|}A_{im}\Sigma^{(k)}_{x\epsilon,mj}\big{\|}_{{\cal
S}}^{\tau}$ $\displaystyle=\sum_{1\leq m,m^{\prime}\leq
r}\big{|}A_{im}A_{jm^{\prime}}\big{|}^{\tau}\big{\|}\Sigma^{(k)}_{xx,mm^{\prime}}\big{\|}_{{\cal
S}}^{\tau}+\sum_{m=1}^{r}\big{|}A_{im}\big{|}^{\tau}\big{\|}\Sigma^{(k)}_{x\epsilon,mj}\big{\|}_{{\cal
S}}^{\tau},$

where, under Conditions 1 and 2, it holds that

$\begin{split}\big{\|}\Sigma^{(k)}_{xx,mm^{\prime}}\big{\|}_{{\cal
S}}^{2}=&\iint[\mathbb{E}X_{(t+k)m}(u)X_{t,m^{\prime}}(v)]^{2}\,{\rm d}u{\rm
d}v\leq\mathbb{E}\|X_{(t+k)m}\|^{2}\,\mathbb{E}\|X_{tm^{\prime}}\|^{2}=O(1),\\\
\max\limits_{1\leq j\leq p}\big{\|}\Sigma^{(k)}_{x\epsilon,mj}\big{\|}_{{\cal
S}}^{2}\leq&\mathbb{E}\|X_{(t+k)m}\|^{2}\max\limits_{1\leq j\leq
p}\mathbb{E}\|\epsilon_{tj}\|^{2}=O(1).\end{split}$ (S.9)

For each $j$, summing over $i$ and applying (S.9) under Condition 6 with fixed
$L$ and $r,$ we obtain that

$\displaystyle\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\big{\|}\Sigma^{(k)}_{yy,ij}\big{\|}^{\tau}_{{\cal S}}$
$\displaystyle\leq\max\limits_{1\leq j\leq p}\sum_{1\leq m,m^{\prime}\leq
r}\sum_{i=1}^{p}\Big{\\{}|A_{im}A_{jm^{\prime}}|^{\tau}\|\Sigma^{(k)}_{xx,mm^{\prime}}\|_{{\cal
S}}^{\tau}\Big{\\}}+\max\limits_{1\leq j\leq
p}\sum_{m=1}^{r}\sum_{i=1}^{p}\Big{\\{}|A_{im}|^{\tau}\|\Sigma^{(k)}_{x\epsilon,mj}\|_{{\cal
S}}^{\tau}\Big{\\}}$ $\displaystyle\leq
c_{1}(p)\Big{\\{}r^{2}\cdot\max\limits_{1\leq j\leq p,\atop 1\leq
m,m^{\prime}\leq
r}|A_{jm^{\prime}}|^{\tau}\|\Sigma^{(k)}_{xx,mm^{\prime}}\|_{{\cal
S}}^{\tau}+r\cdot\max\limits_{1\leq j\leq p,\atop 1\leq m\leq
r}\|\Sigma^{(k)}_{x\epsilon,mj}\|_{{\cal S}}^{\tau}\Big{\\}}\lesssim
c_{1}(p).$

By the similar arguments above and Condition 4′,

$\displaystyle\max\limits_{1\leq i\leq
p}\sum_{j=1}^{p}\big{\|}\Sigma^{(k)}_{yy,ij}\big{\|}^{\tau}_{{\cal S}}$
$\displaystyle\leq\max\limits_{1\leq i\leq p}\sum_{1\leq m,m^{\prime}\leq
r}\sum_{j=1}^{p}\Big{\\{}|A_{im}A_{jm^{\prime}}|^{\tau}\|\Sigma^{(k)}_{xx,mm^{\prime}}\|_{{\cal
S}}^{\tau}\Big{\\}}+\max\limits_{1\leq i\leq
p}\sum_{m=1}^{r}\sum_{j=1}^{p}\Big{\\{}|A_{im}|^{\tau}\|\Sigma^{(k)}_{x\epsilon,mj}\|_{{\cal
S}}^{\tau}\Big{\\}}$ $\displaystyle\leq c_{1}(p)r^{2}\cdot\max_{1\leq i\leq
p,\atop 1\leq m,m^{\prime}\leq
r}|A_{im}|^{\tau}\|\Sigma^{(k)}_{xx,mm^{\prime}}\|_{{\cal
S}}^{\tau}+c_{2}(p)r\cdot\max\limits_{1\leq i\leq p,\atop 1\leq m\leq
r}|A_{im}|^{\tau}$ $\displaystyle\lesssim c_{1}(p)+c_{2}(p).$

Hence we complete the proof of this lemma. $\square$

###  B.4 Proof of Theorem 3

We organize our proof in the following three steps.

Step 1. With the choice of $\eta_{k}\asymp(\log p/n)^{1/2},$ we will show that

$\displaystyle\big{\|}\mathcal{T}_{\eta_{k}}\\{\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}-\boldsymbol{\Sigma}^{(k)}_{yy}\big{\|}_{{\cal
S},1}$
$\displaystyle=O_{P}\left\\{c_{1}(p)\mathcal{M}_{y}^{1-\tau}\Big{(}\frac{\log
p}{n}\Big{)}^{\frac{1-\tau}{2}}\right\\},$ (S.10)
$\displaystyle\big{\|}\mathcal{T}_{\eta_{k}}\\{\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}-\boldsymbol{\Sigma}^{(k)}_{yy}\big{\|}_{{\cal
S},\infty}$
$\displaystyle=O_{P}\left[\big{\\{}c_{1}(p)+c_{2}(p)\big{\\}}\mathcal{M}_{y}^{1-\tau}\Big{(}\frac{\log
p}{n}\Big{)}^{\frac{1-\tau}{2}}\right].$

By the definition of functional matrix $\ell_{1}$ norm $\|\cdot\|_{{\cal
S},1}$ and the triangle inequality,

$\displaystyle\big{\|}\mathcal{T}_{\eta_{k}}\\{\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}-\boldsymbol{\Sigma}^{(k)}_{yy}\big{\|}_{{\cal
S},1}\leq\big{\|}\mathcal{T}_{\eta_{k}}\\{\boldsymbol{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}-\boldsymbol{\Sigma}^{(k)}_{yy}\big{\|}_{{\cal
S},1}+\big{\|}\mathcal{T}_{\eta_{k}}\\{\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}-\mathcal{T}_{\eta_{k}}\\{\boldsymbol{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}\big{\|}_{{\cal
S},1}.$

For the first term, it follows from Lemma 1 that

$\displaystyle\big{\|}\mathcal{T}_{\eta_{k}}\\{\boldsymbol{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}-\boldsymbol{\Sigma}^{(k)}_{yy}\big{\|}_{{\cal
S},1}$ $\displaystyle=\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\|\Sigma^{(k)}_{yy,ij}\|_{{\cal
S}}I\big{\\{}\|\Sigma^{(k)}_{yy,ij}\|_{{\cal S}}<\eta_{k}\big{\\}}$ (S.11)
$\displaystyle<\eta_{k}^{1-\tau}\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\|\Sigma^{(k)}_{yy,ij}\|_{{\cal
S}}^{\tau}\lesssim\eta_{k}^{1-\tau}c_{1}(p).$

For the second term, observe that

$\displaystyle\big{\|}\mathcal{T}_{\eta_{k}}\\{\widehat{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}-\mathcal{T}_{\eta_{k}}\\{\boldsymbol{\boldsymbol{\Sigma}}^{(k)}_{yy}\\}\big{\|}_{{\cal
S},1}$ $\displaystyle=\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\big{\|}\widehat{\Sigma}_{yy,ij}^{(k)}I\\{\|\widehat{\Sigma}_{yy,ij}^{(k)}\|_{{\cal
S}}\geq\eta_{k}\\}-\Sigma_{yy,ij}^{(k)}I\\{\|\Sigma_{yy,ij}^{(k)}\|_{{\cal
S}}\geq\eta_{k}\\}\big{\|}_{{\cal S}}$ $\displaystyle\leq\max\limits_{1\leq
j\leq
p}\sum_{i=1}^{p}\big{\|}\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}I\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}\|_{{\cal
S}}\geq\eta_{k},\|\Sigma_{yy,ij}^{(k)}\|_{{\cal S}}\geq\eta_{k}\big{\\}}$
$\displaystyle~{}~{}~{}~{}+\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\big{\|}\widehat{\Sigma}_{yy,ij}^{(k)}\|_{{\cal
S}}I\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}\geq\eta_{k},\|\Sigma_{yy,ij}^{(k)}\|_{{\cal S}}<\eta_{k}\big{\\}}$
$\displaystyle~{}~{}~{}~{}+\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\big{\|}\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}I\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}\|_{{\cal
S}}<\eta_{k},\|\Sigma_{yy,ij}^{(k)}\|_{{\cal S}}\geq\eta_{k}\big{\\}}$
$\displaystyle:=I_{1}+I_{2}+I_{3}.$

Denote $Z=\max\limits_{1\leq i,j\leq
p}\big{\|}\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}$. We first bound $I_{1}.$ By Lemma 1,

$\displaystyle I_{1}$ $\displaystyle\leq\max\limits_{1\leq i,j\leq
p}\big{\|}\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}\cdot\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}I\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}\|_{{\cal
S}}\geq\eta_{k},\|\Sigma_{yy,ij}^{(k)}\|_{{\cal S}}\geq\eta_{k}\big{\\}}$
$\displaystyle\leq Z\cdot\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}I\big{\\{}\|\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}\geq\eta_{k}\\}\lesssim Z\eta_{k}^{-\tau}c_{1}(p).$

We next bound $I_{2}.$ By the triangle inequality and Lemma 1,

$\displaystyle I_{2}$ $\displaystyle\leq\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\big{\|}\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}I\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}\|_{{\cal
S}}\geq\eta_{k},\|\Sigma_{yy,ij}^{(k)}\|_{{\cal S}}<\eta_{k}\big{\\}}$
$\displaystyle~{}~{}~{}~{}+\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\big{\|}\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}I\big{\\{}\|\Sigma_{yy,ij}^{(k)}\|_{{\cal S}}<\eta_{k}\big{\\}}$
$\displaystyle\lesssim\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\big{\|}\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}I\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}\|_{{\cal
S}}\geq\eta_{k},\|\Sigma_{yy,ij}^{(k)}\|_{{\cal
S}}<\eta_{k}\big{\\}}+\eta_{k}^{1-\tau}c_{1}(p)$
$\displaystyle:=I_{4}+\eta_{k}^{1-\tau}c_{1}(p).$

We take certain $\theta\in(0,1).$ Let $N(t)=\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}I\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\|_{{\cal
S}}\geq t\big{\\}}.$ By the triangle inequality and Lemma 1,

$\displaystyle I_{4}$ $\displaystyle\leq\max\limits_{1\leq j\leq
p}\sum_{i=1}^{p}\big{\|}\widehat{\Sigma}_{yy,ij}^{(k)}-\Sigma_{yy,ij}^{(k)}\big{\|}_{{\cal
S}}I\big{\\{}\|\widehat{\Sigma}_{yy,ij}^{(k)}\|_{{\cal
S}}\geq\eta_{k},\|\Sigma_{yy,ij}^{(k)}\|_{{\cal