# Factor-augmented Smoothing Model for Functional Data

Yuan Gao
Research School of Finance, Actuarial Studies and Statistics
Australian National University

Han Lin Shang
Department of Actuarial Studies and Business Analytics
Macquarie University

Yanrong Yang
Research School of Finance, Actuarial Studies and Statistics
Australian National University
Postal address: Research School of Finance, Actuarial Studies and Statistics,
Level 4, Building 26C, Kingsley St, Australian National University, Canberra,
ACT 2601, Australia; Email<EMAIL_ADDRESS>
###### Abstract

We propose modeling raw functional data as a mixture of a smooth function and
a high-dimensional factor component. The conventional approach to retrieving
the smooth function from the raw data is through various smoothing techniques.
However, the smoothing model is not adequate to recover the smooth curve or
capture the data variation in some situations. These include cases where there
is a large amount of measurement error, the smoothing basis functions are
incorrectly identified, or the step jumps in the functional mean levels are
neglected. To address these challenges, a factor-augmented smoothing model is
proposed, and an iterative numerical estimation approach is implemented in
practice. Including the factor model component in the proposed method solves
the aforementioned problems since a few common factors often drive the
variation that cannot be captured by the smoothing model. Asymptotic theorems
are also established to demonstrate the effects of including factor structures
on the smoothing results. Specifically, we show that the smoothing
coefficients projected on the complement space of the factor loading matrix is
asymptotically normal. As a byproduct of independent interest, an estimator
for the population covariance matrix of the raw data is presented based on the
proposed model. Extensive simulation studies illustrate that these factor
adjustments are essential in improving estimation accuracy and avoiding the
curse of dimensionality. The superiority of our model is also shown in
modeling Canadian weather data and Australian temperature data.

Keywords: Basis function misspecification; Functional data smoothing; High-
dimensional factor model; Measurement error; Statistical inference on
covariance estimation

##  1 Introduction

With the increasing capability to store data, functional data analysis (FDA)
has received growing attention over the last 20 years. Functional data are
considered realizations of smooth random objects in graphical representations
of curves, images, and shapes. The monographs of Ramsay & Silverman (2002,
2005) and Ramsay & Hooker (2017) provide a comprehensive account of the
methodology and applications of the FDA; other relevant monographs include
Ferraty & Vieu (2006) and Horváth & Kokoszka (2012). More recent advances in
this field can be found in many survey papers (see, e.g., Cuevas, 2014;
Febrero-Bande et al., 2017; Goia & Vieu, 2016; Reiss et al., 2017; Wang et
al., 2016). One main challenge in the FDA lies in the fact that we cannot
observe functional curves directly, but only discrete points, which are often
contaminated by measurement errors. To model a mixture of functional data and
high-dimensional measurement error, we introduce a factor-augmented smoothing
model (FASM).

We denote a random sample of $n$ functional data as
$\mathcal{X}_{i}(u),i=1,\dots,n$, and $u\in\mathcal{I}\subset\mathbb{R}$,
where $\mathcal{I}$ is a compact interval on the real line $\mathbb{R}$. In
practice, the observed data are discrete points and are often contaminated by
noise or measurement error. We use $Y_{ij}$ to represent the $j$th observation
on the $i$th subject; the observed data can then be expressed as a “signal
plus noise” model:

$Y_{ij}=\mathcal{X}_{i}(u_{j})+\eta_{ij},\quad j=1,\dots p,\ i=1,\ldots,n.$

We use $\mathcal{X}_{i}(u_{j})$ to denote the realization of the $j$th
discrete point on the curve $\mathcal{X}_{i}(\cdot)$, and $\eta_{ij}$ is the
noise or measurement error. We assume that measurement error only occurs where
the measurements are taken; thus, the error
$\bm{\eta}_{i}=(\eta_{i1},\dots,\eta_{ip})$ is a multivariate term of
dimension $p$. Though in practice, the signal function component
$\bm{\mathcal{X}}_{i}=(\mathcal{X}_{i}(u_{1}),\dots,\mathcal{X}_{i}(u_{j}))$
is of the same $p$ dimension, it differs from $\bm{\eta}_{i}$ in nature.
Although functions are potentially infinite-dimensional, we may impose
smoothing assumptions on the functions, which usually implies functions
possess one or more derivatives. This smoothness feature is used to separate
the functions from measurement errors – a functional smoothing procedure.

When the variance of the noise level is a tiny fraction of the variance of the
function, we say the signal-to-noise ratio is high. In this case, classic
smoothing tools apply to functional data, including kernel methods (e.g. Wand
& Jones, 1995), local polynomial smoothing (e.g. Fan & Gijbels, 1996), and
spline smoothing (e.g. Wahba, 1990; Eubank, 1999; Green & Silverman, 1999).
With pre-smoothed functions, estimates, such as mean and covariance functions,
can be further obtained. More recent studies on functional smoothing
approaches include Cai & Yuan (2011); Yao & Li (2013), and Zhang & Wang
(2016). In this article, we apply basis smoothing to the functions
$\mathcal{X}_{i}(u)$; that is, we represent $\mathcal{X}_{i}(u)$ as
$\mathcal{X}_{i}(u)=\sum_{k=1}^{K}c_{ik}\phi_{k}(u)$, where $\\{\phi_{k}(u),\
k=1,\dots,K\\}$ are the basis functions and $\\{c_{ik},\ i=1,\dots,n,\
k=1,\dots,K\\}$ are the smoothing coefficients. The smoothing model then
becomes

$Y_{ij}=\sum_{k=1}^{K}c_{ik}\phi_{k}(u_{j})+\eta_{ij},\quad j=1,\dots p,\
i=1,\ldots,n.$

When the signal-to-noise level is low, smoothing tools may not be adequate in
removing the measurement error and may cause an inefficient estimation of the
smoothing coefficients. Let us take a further look at the measurement error
$\eta_{ij}$. In the FDA, the number of discrete points $p$ on each subject is
often large compared with the sample size $n$. Hence the term $\bm{\eta}_{i}$
is a high-dimensional component. In this case, the observed data are, in fact,
a mixture of functional data and high-dimensional data. The existence of the
large measurement error $\eta_{ij}$ raises the curse of dimensionality
problem, which naturally calls for the application of dimension reduction
models to $\eta_{ij}$. Many studies have been conducted on various dimension
reduction techniques for high-dimensional data; among theses, factor models
are widely used (e.g. Fan et al., 2008; Lam et al., 2011).

We propose using a factor model for the measurement error term. Without
further information on the measurement error, factor model is appropriate
since the estimation of latent factors does not require any observed
variables. The high-dimensional measurement error is assumed to be driven by a
small number of unobserved common factors.

$\eta_{ij}=\bm{a}_{j}^{\top}\bm{f}_{i}+\epsilon_{ij},\quad i=1,\dots,n,\
j=1,\dots,p,$

where $\bm{f}_{i}\in\mathbb{R}^{r}$ are the unobserved factors,
$\bm{a}_{j}\in\mathbb{R}^{r}$ are the unobserved factor loadings, $r$ is the
number of latent factors, and $\epsilon_{ij}$ are idiosyncratic errors with
mean zero. Thus, the observed data $Y_{ij}$ can be written as the sum of two
components:

$Y_{ij}=\sum_{k=1}^{K}c_{ik}\phi_{k}(u_{j})+\bm{a}_{j}^{\top}\bm{f}_{i}+\epsilon_{ij},\quad
i=1,\dots,n,\ j=1,\dots,p.$

This is a basis smoothing model with the factor-augmented form. This proposed
model can be easily modified to adopt nonparametric smoothing methods. In
Section 5, we illustrate the use of spline smoothing approaches. In Section
7.5, the nonparametric smoothing model is applied to simulated data.

In this paper, we motivate the FASM in three considerations, as listed below.
In these three cases, using the proposed model remedies the defects of the
traditional smoothing model. Examples of the following three motivations are
provided in Section 2.

  1. 1.

In traditional smoothing models, the measurement error $\eta_{ij}$ is assumed
to be non-informative and independently and identically distributed (i.i.d.)
in both directions. This is an unrealistic assumption when the measurement
errors contain information. With the factor model applied, we assume that a
small number of unobserved factors can capture the covariance in the
measurement error. This is usually reasonable in practice because a few common
factors often drive the occurrence of systematic measurement error.

  2. 2.

When the smoothing basis functions are incorrectly identified, the smoothing
model will lead to an erroneous coefficient estimate and large residuals. The
proposed model deals with this problem since the unexplained variation
resulting from the basis’s misidentification can be modeled with a small
number of unobserved common factors.

  3. 3.

When there are step jumps in the mean level of the functions, neglecting the
mean shift in smoothing models will result in large residuals at the point
where the jumps occur. The changes in the mean levels of the functions come
from a universal source and can be modeled by common factors.

Since the latent factors are unobserved, we propose an iterative approach to
estimate the smooth function and the factors simultaneously. Principal
component analysis (PCA) is used as a tool in estimating the factor model, and
penalized least squares estimation is applied to construct the estimator for
the smoothing coefficient $c_{ik}$. We establish the asymptotic theories of
the smoothing coefficient estimator, where the consistency of the estimator is
proved. We also provide the asymptotic distribution of the projected estimator
in the orthogonal complement of the space spanned by the factors $\bm{f}_{i}$.
The interplay between the smooth component and the factor model component is
manifested.

In the remainder of this article, we elaborate on the previously mentioned
three motivations in detail, with examples given in Section 2. In Section 3,
the model is formally stated, and the iterative estimation approach is
provided. We discuss the asymptotic properties of the smoothing coefficients
under various assumptions in Section 4. We extend the proposed model to a
nonparametric smoothing approach in Section 5. In Section 6, we consider the
statistical inference aspect of the model and propose a covariance matrix
estimator for the raw data. In Section 7, we conduct Monte-Carlo simulations
on the proposed model under different settings. A few real data examples are
given in Section 8, and conclusions are drawn in Section 9. Last, we provide
proofs of the relevant theorems and lemmas in the Appendix.

##  2 Motivation

We introduce three examples to motivate the proposed model. In these cases,
the smoothing model is not adequate to capture the raw data’s signal
information. In the first example, when large measurement error exists, the
residuals after smoothing are large with some extreme values. In the second
example, when the basis functions are selected incorrectly, part of the
functions’ variation cannot be captured by the smoothing model. In the third
example, when there are step jumps in the functional data, the residuals after
smoothing contain gaps. These examples demonstrate that further modeling of
the residuals is needed.

###  2.1 Functional data with measurement error

Figure 1 shows the rainbow plots of the average daily temperature and log
precipitation at 35 locations in Canada. Due to the nature of the two kinds of
data, it is reasonable to assume that temperature and log precipitation are
functions over time. The two graphs, however, display distinct features. In
the temperature plot, though there are some perturbations, it is relatively
easy to discern each curve’s shape. In the precipitation plot, there is a
tremendous amount of variability in the raw data, such that it is almost
impossible to observe the underlying shape of the curves.

Smooth temperature data can be retrieved without much difficulty using basic
smoothing techniques. The residuals are small, with constant variation. On the
other hand, the residuals after smoothing exhibit a high level of variation
for the precipitation data and even contain some extreme values. Our model
endeavors to further explain the large residuals in similar cases to the
precipitation data; we will show the fitting result in Section 8.

Figure 1: Average daily temperature and log precipitation in 35 Canadian
weather stations averaged over the year 1960 to 1994.

###  2.2 Misidentification of the basis function

It is important to choose the appropriate basis functions in the smoothing
method. In this example, we show the inadequacy of the smoothing model when
the basis functions are misidentified. We generate functional data using basis
functions with changing frequencies. The raw data are shown in Figure 2(a).
Fourier basis functions are used. In the second half of the data, the
frequency of the Fourier basis functions increases, so the data set exhibits
more variation toward the right end. Suppose that we were not aware of the
change in the frequencies in the basis functions, and still used the basis of
the first half of the data for the whole curves. The consequence of
misidentifying the basis functions when a smoothing model is applied can be
observed in Figure 2(b). The residuals are large in the second half. The
smoothing model fails to reduce the residuals; a factor model can be used to
further model the signal hidden in the large residuals. The data generation
process and further analysis can be found in Section 7.6.

(a) Raw data  (b) Residuals

Figure 2: A simulated sample of functional data with changing basis functions.

###  2.3 Functional data with step jumps in the mean level

We provide another example of functional data with step jumps to motivate our
proposed model. Suppose we observed a sample of the raw functional data, as
shown in Figure 3(a). It can be seen that there is a jump at around $u=0.5$.
The jump applies to all the sample data, so this sudden shift is at the mean
level. We will explain how the data are generated in Section 7.7. The
residuals after smoothing are presented in Figure 3(b). The large residuals
around the jump clarify that without measures to deal with the step jumps,
smoothing itself is not enough to model these kinds of data. We show in
Section 7.7 that the proposed model applied to the same data generates smaller
residuals and has less flexibility. This is indeed one of the main goals
because of model selection.

(a) Raw data (b) Residuals

Figure 3: A simulated sample of functional data with step jump.

##  3 Model specification and estimation

In this section, we formally state the proposed model in Section 3.1 and
provide the estimation method in Section 3.2. We first show how the smoothing
coefficient $\bm{c}_{i}$ and the latent factors $\bm{f}_{i}$ are estimated
separately and then introduce an iterative approach to simultaneously find
these estimates.

###  3.1 Model statement

We consider a sample of functional data $\mathcal{X}_{i}(u)$, which takes
values in the space $H:=L^{2}(\mathcal{I})$ of real-valued square integrable
functions on $\mathcal{I}$. The space $H$ is a Hilbert space, equipped with
the inner product $\langle x,y\rangle:=\int x(u)y(u)du$. The function norm is
defined as $\|x\|:=\langle x,x\rangle^{1/2}.$ The functional nature of
$\mathcal{X}_{i}(u)$ allows us to represent it as a linear expansion of a set
of $K$ smooth basis functions.

$\mathcal{X}_{i}(u)=\sum_{k=1}^{K}c_{ik}\phi_{k}(u),\quad u\in\mathcal{I},$

where $\phi_{k}(u)$ is a set of common basis functions and $c_{ik}$ is the
$k$th coefficient for the $i$th curve. Therefore, we can express the full
model as

$\displaystyle Y_{ij}$
$\displaystyle=\sum_{k=1}^{K}c_{ik}\phi_{k}(u_{j})+\eta_{ij},$
$\displaystyle\eta_{ij}$
$\displaystyle=\bm{a}_{j}^{\top}\bm{f}_{i}+\epsilon_{ij},\quad i=1,\dots,n,\
j=1,\dots,p,$

where $\bm{f}_{i}\in\mathbb{R}^{r}$ are the unobserved common factors,
$\bm{a}_{j}\in\mathbb{R}^{r}$ are the unobserved factor loadings and $r$ is
the number of factors. We call this model the FASM. For the model to be
identifiable, we require the following condition.

######  Identification Condition 1.

We require

  1. (i)

$\mathcal{X}_{i}(u_{j})$ is independent of $\eta_{ij}$ for $i=1,\dots,n,\
j=1,\dots,p$, and

  2. (ii)

$\frac{1}{p}\sum_{j=1}^{p}\bm{a}_{j}\bm{a}_{j}^{\top}\overset{p}{\to}\Sigma_{\bm{a}}>0$
for some $r\times r$ matrix $\Sigma_{\bm{a}}$, as $p\rightarrow\infty$;

$\frac{1}{n}\sum_{i=1}^{n}\bm{f}_{i}\bm{f}_{i}^{\top}\overset{p}{\to}\Sigma_{\bm{f}}>0$
for some $r\times r$ matrix $\Sigma_{\bm{f}}$, as $n\rightarrow\infty$.

The first part of the identification condition ensures the signal function
component and the factor model component are independent. The second part
ensures the existence of $r$ factors, each of which makes a non-trivial
contribution to the variance of $\eta_{ij}$, which in turn guarantees the
identifiability between the factors and the error term $\epsilon_{ij}$.

We treat the basis functions $\phi_{k}(u)$ as known, and the number $K$ fixed.
This is, of course, a simplification to accommodate for the theoretical
proofs. In real data analysis, there are various choices for the basis
functions, and the decision can be quite subjective. For example, Fourier
bases are preferred for periodic data, while spline basis systems are most
commonly used for non-periodic data. Other bases include wavelet, polynomial,
and some ad-hoc basis functions.

###  3.2 Estimation

We can write the model for the $i$th object as

$\bm{Y}_{i}=\bm{\Phi}\bm{c}_{i}+\bm{A}\bm{f}_{i}+\bm{\epsilon}_{i},\quad
i=1,\dots,n$ (1)

where

$\displaystyle\bm{Y}_{i}=\begin{bmatrix}Y_{i1}\\\ \vdots\\\
Y_{ip}\end{bmatrix},\quad\bm{c}_{i}=\begin{bmatrix}c_{i1}\\\ \vdots\\\
c_{iK}\end{bmatrix},\quad\bm{\Phi}=\begin{bmatrix}\phi_{1}(u_{1})&\dots&\phi_{K}(u_{1})\\\
\vdots&&\vdots\\\
\phi_{1}(u_{p})&\dots&\phi_{K}(u_{p})\end{bmatrix},\quad\bm{A}=\begin{bmatrix}\bm{a}_{1}^{\top}\\\
\vdots\\\
\bm{a}_{p}^{\top}\end{bmatrix},\quad\bm{\epsilon}=\begin{bmatrix}\epsilon_{i1}\\\
\vdots\\\ \epsilon_{ip}\end{bmatrix}.$

Combining all the objects, we have in matrix form

$\bm{Y}=\bm{\Phi}\bm{C}+\bm{A}\bm{F}^{\top}+\bm{E},$ (2)

where $\bm{Y}$ is $p\times n$ and $\bm{C}=(\bm{c}_{1},\dots,\bm{c}_{n})$ is a
$K\times n$ matrix containing all the smoothing coefficients. The matrix
$\bm{F}=(\bm{f}_{1},\dots,\bm{f}_{n})^{\top}$ is $n\times r$ and
$\bm{E}=(\bm{\epsilon}_{1},\dots,\bm{\epsilon}_{n})$ is $p\times n.$ Since
$\bm{\Phi}$ is assumed to be known, we illustrate how the parameters
$\bm{C},\bm{A}$ and $\bm{f}$ are estimated in the following.

For the latent factor estimation, there is an identification problem such that
$\bm{A}\bm{F}^{\top}=\bm{A}\bm{U}\bm{U}^{-1}\bm{F}^{\top}$ for any $r\times r$
invertible matrix $\bm{U}$. Thus we impose the normalization restriction on
the matrices $\bm{A}$ and $\bm{F}$

$\bm{A}^{\top}\bm{A}/p=\bm{I}_{r},\quad\text{and}\ \bm{F}^{\top}\bm{F}\
\text{is a diagonal matrix}.$ (3)

We propose to implement penalized least squares, where the objective function
is defined as

$\text{SSR}(\bm{c}_{i},\bm{A},\bm{f})=\sum_{i=1}^{n}\left[(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i}-\bm{A}\bm{f}_{i})^{\top}(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i}-\bm{A}\bm{f}_{i})+\alpha\times\text{PEN}_{2}(\mathcal{X}_{i})\right],$

where $\text{PEN}(\mathcal{X}_{i})$ is a penalty term used for regularization,
and $\alpha$ is the tuning parameter controlling the degree of regularization.
The same $\alpha$ is used for all the functional observations $i$. This is a
simplified case, where we assume a similar degree of smoothness for all
curves. The tuning parameter can be chosen by cross-validation or information
criteria. We intend to penalize the “roughness” of the function term. To
quantify the notion of “roughness” in a function, we use the square of the
second derivative. Define the measure of roughness as

$\text{PEN}_{2}(\mathcal{X}_{i})=\int_{\mathcal{I}}\left[D^{2}\mathcal{X}_{i}(s)\right]^{2}ds,$

where $D^{2}\mathcal{X}_{i}$ denotes taking the second derivative of the
function $\mathcal{X}_{i}$, the larger the tuning parameter $\alpha$, the
smoother the estimated functions we obtain. Further, we denote

$\bm{\Phi}(u)=\left[\phi_{1}(u),\dots,\phi_{K}(u)\right]^{\top}.$ (4)

Then

$\mathcal{X}_{i}(u)=\bm{c}_{i}^{\top}\bm{\Phi}(u).$

We can re-express the roughness penalty $\text{PEN}_{2}(\mathcal{X}_{i})$ in
matrix form as the following:

$\displaystyle\text{PEN}_{2}(\mathcal{X}_{i})$
$\displaystyle=\int_{\mathcal{I}}\left[D^{2}\mathcal{X}_{i}(s)\right]^{2}ds$
$\displaystyle=\int_{\mathcal{I}}\left[D^{2}\bm{c}_{i}^{\top}\bm{\Phi}(s)\right]^{2}ds$
$\displaystyle=\int_{\mathcal{I}}\bm{c}_{i}^{\top}D^{2}\bm{\Phi}(s)D^{2}\bm{\Phi}^{\top}(s)\bm{c}_{i}ds$
$\displaystyle=\bm{c}_{i}^{\top}\left[\int_{\mathcal{I}}D^{2}\bm{\Phi}(s)D^{2}\bm{\Phi}^{\top}(s)ds\right]\bm{c}_{i}$
$\displaystyle=\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i},\qquad i=1,\dots,n$

where

$\bm{R}\equiv\int_{\mathcal{I}}D^{2}\bm{\Phi}(s)D^{2}\bm{\Phi}^{\top}(s)ds.$
(5)

Thee matrix $\bm{R}$ is the same for all subjects and the penalty term
$\text{PEN}_{2}(\mathcal{X}_{i})$ differs for each subject only by the
coefficient $\bm{c}_{i}$.

######  Remark 1.

The number of smoothing coefficient $\bm{c}_{i}$ needed increases as the
sample size increases. The inclusion of a penalty term not only penalizes the
”roughness” of the smoothed function but also mitigates the effect of
increasing number of parameters to control the model flexibility.

Thus, the objective function can be written as

$\text{SSR}(\bm{c}_{i},\bm{A},\bm{f})=\sum_{i=1}^{n}\left[(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i}-\bm{A}\bm{f}_{i})^{\top}(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i}-\bm{A}\bm{f}_{i})+\alpha\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i}\right],$

subject to the constraint $\bm{A}^{\top}\bm{A}/p=\bm{I}_{r}$.

We aim to estimate the smoothing coefficient of $\bm{c}_{i}$. We left multiply
a matrix to each term in (1) to project the factor model term onto a zero
matrix. Define the projection matrix

$\bm{M}_{\bm{A}}\equiv\bm{I}_{p}-\bm{A}(\bm{A}^{\top}\bm{A})^{-1}\bm{A}^{\top}=\bm{I}_{p}-\bm{A}\bm{A}^{\top}/p.$
(6)

Then

$\bm{M}_{\bm{A}}\bm{A}\bm{f}_{i}=\left(\bm{I}_{p}-\bm{A}\bm{A}^{\top}/p\right)\bm{A}\bm{f}_{i}=\left(\bm{A}-\bm{AA}^{\top}\bm{A}/p\right)\bm{f}_{i}=\bm{0}.$

So we estimate $\bm{c}_{i}$ from the projected equation

$\bm{M}_{\bm{A}}\bm{Y}_{i}=\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}+\bm{M}_{\bm{A}}\bm{\epsilon}_{i}.$

The projected objective function becomes

$\text{SSR}(\bm{c}_{i},\bm{A})=\sum_{i=1}^{n}\left[(\bm{M}_{\bm{A}}\bm{Y}_{i}-\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i})^{\top}(\bm{M}_{\bm{A}}\bm{Y}_{i}-\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i})+\alpha\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i}\right].$
(7)

By taking the derivative of $\text{SSR}(\bm{c}_{i},\bm{A})$ with respective to
each $\bm{c}_{i}$, we can solve for the estimator $\widehat{\bm{c}}_{i}$.

$\frac{\partial\text{SSR}(\bm{c}_{i}\bm{A})}{\partial\bm{c}_{i}}=(\bm{M}_{\bm{A}}\bm{Y}_{i}-\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i})^{\top}(\bm{M}_{\bm{A}}\bm{\Phi})+2\alpha\bm{c}_{i}^{\top}\bm{R}.$

Setting the derivative to zero and rearranging the terms, we have

$\left(\bm{\Phi}^{\top}\bm{M}_{\bm{A}}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}+\alpha\bm{R}\right)\bm{c}_{i}=\bm{\Phi}^{\top}\bm{M}_{\bm{A}}^{\top}\bm{M}_{\bm{A}}\bm{Y}_{i}.$

Using the fact that

$\bm{M}_{\bm{A}}^{\top}\bm{M}_{\bm{A}}=\left(\bm{I}_{p}-\bm{A}\bm{A}^{\top}/p\right)^{\top}\left(\bm{I}_{p}-\bm{A}\bm{A}^{\top}/p\right)=\bm{M}_{\bm{A}},$

we obtain the least squares estimator for $\bm{c}_{i}$ given $\bm{A}$

$\widehat{\bm{c}}_{i}=\left(\bm{\Phi}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{\bm{A}}\bm{Y}_{i}.$

Next, to estimate $\bm{A}$ and $\bm{f}_{i}$, we focus on the factor model

$\bm{\eta}_{i}=\bm{A}\bm{f}_{i}+\bm{\epsilon}_{i},$

and in matrix form

$\bm{Z}=\bm{A}\bm{F}^{\top}+\bm{E},$

where $\bm{Z}=(\bm{\eta}_{1},\dots,\bm{\eta}_{n})$. In high dimensions, the
unknown factors and loadings are typically estimated by least squares (i.e.,
the principal component analysis; see, e.g., Fan et al. (2008); Onatski
(2012). The least squares objective function is

$\text{tr}\left[(\bm{Z}-\bm{AF}^{\top})(\bm{Z}-\bm{AF}^{\top})^{\top}\right].$
(8)

Minimizing the objective function with respect to $\bm{F}^{\top}$, we have
$\bm{F}^{\top}=(\bm{A}^{\top}\bm{A})^{-1}\bm{A}^{\top}\bm{Z}=\bm{A}^{\top}\bm{Z}/p$
using (3). Substituting in (8), we obtain the objective function

$\displaystyle\text{tr}\left[(\bm{Z}-\bm{AA}^{\top}\bm{Z}/p)(\bm{Z}-\bm{AA}^{\top}\bm{Z}/p)^{\top}\right]$
$\displaystyle=\text{tr}\left(\bm{ZZ}^{\top}-\bm{ZZ}^{\top}\bm{AA}^{\top}/p-\bm{ZZ}^{\top}\bm{AA}^{\top}/p+\bm{AA}^{\top}\bm{ZZ}^{\top}\bm{AA}^{\top}/p^{2}\right)$
$\displaystyle=\text{tr}(\bm{ZZ}^{\top})-\text{tr}(\bm{A}^{\top}\bm{ZZ}^{\top}\bm{A})/p,$

where the last equality uses (3) and that
$\text{tr}(\bm{ZZ}^{\top}\bm{AA}^{\top})=\text{tr}(\bm{A}^{\top}\bm{ZZ}^{\top}\bm{A})$.
Thus, minimizing the objective function is equivalent to maximizing
$\text{tr}(\bm{A}^{\top}\bm{ZZ}^{\top}\bm{A})/p$. The estimator for $\bm{A}$
is obtained by finding the first $r$ eigenvectors corresponding to the $r$
largest eigenvalues of the matrix $\bm{ZZ}^{\top}$ in descending order, where

$\bm{Z}\bm{Z}^{\top}=\sum_{i=1}^{n}\bm{\eta}_{i}\bm{\eta}_{i}^{\top}=\sum_{i=1}^{n}(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i})(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i})^{\top}.$

Therefore, knowing $\bm{c}_{i}$, we solve for $\widehat{\bm{A}}$ using

$\left[\frac{1}{np}\sum_{i=1}^{n}(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i})(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i})^{\top}\right]\widehat{\bm{A}}=\widehat{\bm{A}}\bm{V}_{np},$
(9)

where $\bm{V}_{np}$ is a $r\times r$ diagonal matrix containing the $r$
eigenvalues of the matrix in the square brackets in decreasing order. The
coefficient $\frac{1}{np}$ is used for scaling.

######  Remark 2.

The number of factors $r$ is assumed to be known in this paper. In practice,
$r$ is selected based on some criteria regarding the eigenvalues. There have
been many studies on this topic. Examples include Bai & Ng (2002), where two
model selection criteria functions were proposed; and Onatski (2010), where
the number of factors was estimated using differenced eigenvalues; and Ahn &
Horenstein (2013), where this number was selected based on the ratio of two
adjacent eigenvalues.

It can be seen that $\bm{A}$ is needed to find $\widehat{\bm{c}}_{i}$, and in
turn $\bm{c}_{i}$ is needed to find $\widehat{\bm{A}}$. The final estimator
$(\widehat{\bm{c}}_{i},\widehat{\bm{A}})$ is the solution of the set of
equations

$\displaystyle\begin{cases}\widehat{\bm{c}_{i}}=\left(\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{Y}_{i},\qquad
i=1,\dots,n\\\
\left[\frac{1}{np}\sum_{i=1}^{n}\left(\bm{Y}_{i}-\bm{\Phi}\widehat{\bm{c}}_{i}\right)\left(\bm{Y}_{i}-\bm{\Phi}\widehat{\bm{c}}_{i}\right)^{\top}\right]\widehat{\bm{A}}=\widehat{\bm{A}}\bm{V}_{np}.\end{cases}$
(10)

Since there is no closed-form expression of $\widehat{\bm{A}}$ and
$\widehat{\bm{c}}_{i}$, we propose using numerical iterations to find the
estimates. The details of these iterations are as follows:

  1. 1.

Denote the initial value as $\widehat{\bm{A}}^{(0)}$. Using (10), we obtain
$\widehat{\bm{c}}_{i}^{(0)}=\left(\bm{\Phi}^{\top}\bm{M}_{{\widehat{\bm{A}}}^{(0)}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{{\widehat{\bm{A}}}^{(0)}}\bm{Y}_{i}$.

  2. 2.

With $\widehat{\bm{c}}_{i}^{(t)}$, we substitute into the second equation of
(10), to obtain
$\widehat{\bm{A}}^{(t+1)}=(\widehat{\bm{a}}_{1}^{(t+1)},\dots,\widehat{\bm{a}}_{r}^{(t+1)})^{\top}$,
where $\widehat{\bm{a}}_{j}^{(t+1)}$ is the eigenvector of the matrix
$\frac{1}{np}\sum_{i=1}^{n}(\bm{Y}_{i}-\bm{\Phi}\widehat{\bm{c}}_{i}^{(t+1)})(\bm{Y}_{i}-\bm{\Phi}\widehat{\bm{c}}_{i}^{(t+1)})^{\top}$
corresponding to its $j$th largest eigenvalue.

  3. 3.

With $\widehat{\bm{A}}^{(t+1)}$, we obtain
$\widehat{\bm{c}}_{i}^{(t+1)}=\left(\bm{\Phi}^{\top}\bm{M}_{{\widehat{\bm{A}}}^{(t+1)}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{{\widehat{\bm{A}}}^{(t+1)}}\bm{Y}_{i}$
using (10)

  4. 4.

We then repeat steps 2 and 3 until
$\|\widehat{\bm{c}}_{i}^{(t+1)}-\widehat{\bm{c}}_{i}^{(t)}\|<\delta$, where
$\delta$ is a small positive constant.

Algorithm 1 Iterations for estimating FASM

######  Remark 3.

In this paper, we use $\widehat{\bm{A}}^{(0)}=\bm{0}$. This means we start by
ignoring the factor model component so the initial value for the smoothing
coefficient
$\widehat{\bm{c}}_{i}^{(0)}=\left(\bm{\Phi}^{\top}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{Y}_{i}$,
which is simply the ridge estimator. The convergence of Newton’s numeric
iteration requires the convergence of this estimator, which in turn requires
the factor model component $\eta_{ij}$ to have an expectation of zero. The
stopping criterion only focuses on $\widehat{\bm{c}}_{i}$ because we are
interested in estimating $\eta_{ij}$ as a whole.

######  Remark 4.

Common methods for selecting the shrinkage parameter $\alpha$ include the
Akaike’s Information Criterion (AIC Akaike, 1974), the Bayesian Information
Criterion (BIC Schwarz, 1978), and cross-validation. In this paper, we use the
mean generalized cross-validation (mGCV) method (Golub et al., 1979). We
define, at step $t$,

$\text{mGCV}^{(t)}=\frac{1}{n}\sum_{i=1}^{n}\frac{\text{pSSE}_{i}^{(t)}}{[p-df^{(t)}(\alpha)]^{2}},$
(11)

where $\text{SSE}_{i}^{(t)}$ is the sum of squares residual for the $i$th
object at step $t$ and $df^{(t)}(\alpha)$ is the equivalent degrees of freedom
measure, which can be calculated as

$\textnormal{df}^{(t)}(\alpha)=\textnormal{trace}\left[\bm{\Phi}\left(\bm{\Phi}^{\top}\bm{M}_{{\widehat{\bm{A}}}^{(t)}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{{\widehat{\bm{A}}}^{(t)}}\right].$
(12)

At each step of the iteration, the tuning parameter $\alpha$ is chosen by
minimizing the $mGCV^{(t)}$.

######  Remark 5.

Algorithm 1 is an iteration procedure in which ridge regression and PCA are
iterated. The convergence of this iterative algorithm is studied in Jiang et
al. (2020). For instance, Theorem 2 of Jiang et al. (2020) provides some
sufficient conditions under which the recursive algorithm converges to the
true value or some other values. In particular, when the regressors are
independent of the common factors, or the factors involved in regressors are
weaker than the common factors, this algorithm will converge to the true
parameter.

After we obtain the estimates $\widehat{\bm{A}}$ and $\widehat{\bm{c}}_{i}$,
the estimated coefficient matrix $\widehat{\bm{C}}$ is constructed as
$\widehat{\bm{C}}=\left(\widehat{\bm{c}}_{1},\dots,\widehat{\bm{c}}_{n}\right)$,
and the estimated factor can be obtained by

$\widehat{\bm{F}}^{\top}=\widehat{\bm{A}}^{\top}(\bm{Y}-\bm{\Phi}\widehat{\bm{C}}).$

Finally, the functional component can be estimated by
$\widehat{\mathcal{X}}_{i}(u)=\widehat{\bm{c}}_{i}^{\top}\bm{\Phi}(u)$, where
$\bm{\Phi}(u)$ is defined in (4).

######  Remark 6.

Although we have imposed the constraint in (3) and the identification
condition 1, $\bm{A}$ and $\bm{f}_{i}$ are not uniquely determined, since the
model (1) is unchanged if we replace $\bm{A}$ and $\bm{f}_{i}$ with $\bm{AU}$
and $\bm{U}^{\top}\bm{f}_{i}$ for any orthogonal $r\times r$ matrix $\bm{U}$.
However, the linear space spanned by the columns of $\bm{A}$ is uniquely
defined. Although we are not able to estimate $\bm{A}$, we can still estimate
a rotation of $\bm{A}$, which spans the same space as $\bm{A}$ does. The
matrix $\bm{M}_{\bm{A}}$ defined in (6) is a projecting matrix onto the
orthogonal complement of the linear space spanned by the columns of $\bm{A}$.
It is shown in the next section that the estimator $\bm{M}_{\widehat{\bm{A}}}$
for $\bm{M}_{\bm{A}}$ is consistent.

##  4 Asymptotic theory

In this section, we study the asymptotic properties of the coefficient
estimator $\widehat{\bm{c}}_{i}$ with growing sample size and dimension. We
state the assumptions in Section 4.1 and provide the asymptotic results of
$\widehat{\bm{c}}_{i}$ in Section 4.2.

###  4.1 Assumptions

We use $(\bm{c}_{i}^{0},\bm{A}^{0})$ to denote the true parameters. In this
paper, the norm of a vector or matrix $\bm{U}$ is defined as the Frobenius
norm; that is, $\|\bm{U}\|=[\text{tr}(\bm{U}^{\top}\bm{U})]^{1/2}$. We
introduce the matrix

$\bm{D}_{i}(\bm{A})\equiv\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{f}_{i}^{\top}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{\top}\bm{f}_{i}.$
(13)

This matrix plays an important role in this article. It is used in the proof
of the consistency of $\widehat{\bm{c}}_{i}$, as can be found in Appendix A.
The identifying condition for $\bm{c}^{0}_{i}$ is that $\bm{D}_{i}(\bm{A})$ is
positive definite for all $i$, which is stated in Assumption 3.

First, we state the assumptions.

######  Assumption 1.

$\sup_{u}|\phi_{k}(u)|=O(1),\quad k=1,\dots,K.$

The above assumption declares that the basis functions are bounded in the
norm. This is quite natural as some of the most commonly used basis functions
are bounded; for instance, the Fourier basis, B-spline basis, and wavelet
basis functions (Ramsay & Silverman, 2005).

######  Assumption 2.

$\|\bm{c}_{i}^{0}\|=O(1),\ \text{for all }i.$

Above, we assume the smoothing coefficients $\bm{c}^{0}_{i}$ are bounded
uniformly for all $i$. This assumption is introduced to ensure the uniform
consistency of the estimated coefficients of $\widehat{\bm{c}}^{0}_{i}$.

######  Assumption 3.

Let $\mathcal{A}=\\{\bm{A}:\bm{A}^{\top}\bm{A}/p=\bm{I},\text{ and
}\bm{A}\text{ independent of }\bm{\Phi}\\}$. We assume

$\inf_{\bm{A}\in\mathcal{A}}\bm{D}_{i}(\bm{A})>0.$

This assumption is the identification condition for $\bm{c}^{0}_{i}$. The
usual assumption for the least-squares estimator only contains the first term
on the right-hand side of (13). The second term on the right-hand side of (13)
arises because of the unobservable matrices $\bm{F}$ and $\bm{A}$.

######  Assumption 4.

For some constant $M>0$, $\mathbb{E}\|\bm{a}^{0}_{j}\|^{4}\leq M$,
$j=1,\dots,p$ and $\mathbb{E}\|\bm{f}_{i}\|^{4}\leq M$.

######  Assumption 5.

For some constant $M>0$, the error terms
$\epsilon_{ji},j=1,\dots,p,i=1,\dots,n$ are i.i.d. in both directions, with
$\mathbb{E}(\epsilon_{ji})=0$, $\text{Var}(\epsilon_{ji})=\sigma^{2}$, and
$\mathbb{E}|\epsilon_{ji}|^{8}\leq M$.

######  Assumption 6.

$\epsilon_{ji}$ is independent of $\phi_{s}$, $\bm{f}_{t}$, and
$\bm{a}^{0}_{s}$ for all $j,i,s,t$.

We require that the errors are independent in themselves and also of the
functional term $\phi(u)$ and factor model terms $\bm{f}_{i}$ and
$\bm{a}^{0}_{j}$. To not mask the main contribution of our method, we use a
simplified setting on the error terms to exclude endogeneity. Nevertheless,
with simple but tedious modifications, Assumption 5 can be relaxed, and our
model can be extended to more complicated settings where correlations between
the error term and the factor model term are allowed.

######  Assumption 7.

The tuning parameter $\alpha$ satisfies $\alpha=o(p)$.

This is conventionally assumed in ridge regression (see, e.g., Knight & Fu,
2000) and assures that the estimator’s asymptotic bias is zero.

Before stating the next assumption, we introduce some notations. Let
$\bm{\omega}_{j},\ j=1,\dots,p$ denote the $j$th column of the $K\times p$
matrix $\bm{\Phi}^{\top}\bm{M}_{{\bm{A}}^{0}}$, and let $\psi_{ik}$ denote the
$(i,k)$th element of the matrix $\bm{M}_{\bm{F}}$, where

$\bm{M}_{\bm{F}}\equiv\bm{I}_{n}-\bm{F}\left(\bm{F}^{\top}\bm{F}\right)\bm{F}^{\top}.$
(14)

Then, for any vector $\bm{b}=(b_{1},\dots,b_{n})^{\top}$, we can write

$\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}\bm{b}=\frac{1}{\sqrt{np}}\sum_{i}^{n}\sum_{j}^{p}\bm{\omega}_{j}\epsilon_{ji}\sum_{k}^{n}\psi_{ik}b_{k}\equiv\frac{1}{\sqrt{np}}\sum_{i}^{n}\sum_{j}^{p}\bm{x}_{ij}.$
(15)

In (15), for notation simplicity, we define $\bm{x}_{ij}$ as
$\bm{\omega}_{j}\epsilon_{ji}\sum_{k}^{n}\psi_{ik}b_{k}$. The matrix
$\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}$ is of interest
because it is the main component that contributes to the asymptotic
distribution of the estimators, as shall be seen in the next section.

Let

$\bm{L}_{np}\equiv\frac{\sigma^{2}}{np}\sum_{i}^{n}\sum_{j}^{p}\bm{\omega}_{j}^{\top}\bm{\omega}_{j}\left(\sum_{k}^{n}\psi_{ik}b_{k}\right)^{2}.$
(16)

We make the following assumption.

######  Assumption 8.

We assume there exists a $K\times K$ matrix $\bm{L}$ such that

$\bm{L}\equiv\lim_{n,p\rightarrow\infty}\bm{L}_{np},$ (17)

where $\bm{L}_{np}$ is defined in (16). Let $\nu^{2}$ be the smallest
eigenvalue of the matrix $\bm{L}$ defined in (17), then assume that
$\nu^{2}>0$, and that, for all $\varepsilon>0$,

$\lim_{n,p\rightarrow\infty}\frac{1}{np\nu^{2}}\sum_{i=1}^{n}\sum_{j=1}^{p}\mathbb{E}\left[\left\|\bm{x}_{ij}\right\|^{2}\mathbb{1}\left(\left\|\bm{x}_{ij}\right\|^{2}\geq\varepsilon
np\nu^{2}\right)\right]=0.$ (18)

This assumption is the multivariate Lindeberg condition, which is needed in
constructing the central limit theorem in the next section. This is by no
means a strong condition; for instance, when the factor model component is
ignored, $\bm{\omega}_{j}$ is simply $\bm{\phi}_{j}$, and
$\bm{x}_{ij}=\bm{\phi}_{j}b_{i}\epsilon_{ji}$. Since we assume
$\bm{\phi}_{j}=O(1)$ in Assumption 1, the Lindeberg condition in (18) is met.

###  4.2 Asymptotic properties

As we have mentioned previously, the identification problem of the latent
factor implies that we actually use the estimator $\widehat{\bm{A}}$ to
estimate a rotation of $\bm{A}^{0}$. Based on the objective function (7) in
Section 3, we use a center-adjusted objective function, defined as below.

$S_{np}(\bm{c}_{i},\bm{A})=\frac{1}{np}\sum_{i=1}^{n}\left[(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i})^{\top}\bm{M}_{\bm{A}}(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i})+\alpha\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i}\right]-\frac{1}{np}\sum_{i=1}^{n}\epsilon_{i}^{\top}\bm{M}_{\bm{A}^{0}}\epsilon_{i},$
(19)

where $\bm{M}_{\bm{A}}$ is defined in (6), satisfying
$\bm{A}^{\top}\bm{A}/p=\bm{I}_{r}.$ The second term on the right-hand side of
(19) does not contain the unknown $\bm{A}$ and $\bm{c}_{i}$, so the inclusion
of this term does not affect the optimization result. This term is only used
for center adjusting, so that the resulting objective function has an
expectation zero. We estimate $\bm{c}_{i}^{0}$ and $\bm{A}^{0}$ by

$(\widehat{\bm{c}}_{i},\widehat{\bm{A}})=\operatorname*{\arg\\!\min}_{\bm{c}_{i},\bm{A}}S_{np}(\bm{c}_{i},\bm{A}).$
(20)

In the following, we establish the asymptotic properties for the estimated
coefficient matrix $\widehat{\bm{C}}$. In Theorem 1, the consistency of the
matrix $\widehat{\bm{C}}$ is proved. In Theorem 2, we show the rate of
convergence of $\widehat{\bm{C}}$. Theorem 3 provides the asymptotic
distribution of $\widehat{\bm{C}}$.

Let $\bm{P}_{\bm{U}}=\bm{U}(\bm{U}^{\top}\bm{U})^{-1}\bm{U}^{\top}$ for a
matrix $\bm{U}$.

######  Theorem 1.

Under Assumptions 1 \- 6, as $n,p\rightarrow\infty$, we have the following
statements

  1. (i)

$\frac{1}{\sqrt{n}}\left\|\bm{C}-\widehat{\bm{C}}\right\|\overset{p}{\to}0.$

  2. (ii)

$\left\|\bm{P}_{\widehat{\bm{A}}}-\bm{P}_{\bm{A}^{0}}\right\|\overset{p}{\to}0$.

We start by proving consistency for the vector $\widehat{\bm{c}}_{i}$. This
consistency is uniform for all $i=1,\dots,n$. Therefore, we could combine
$\bm{\bm{c}}_{i}$ for all $i=1,\dots,n$, and have the result for the
coefficient matrix $\widehat{\bm{C}}$ in $(i)$. The matrix $\widehat{\bm{C}}$
is of dimension $K\times n$, where $K$ is fixed and the sample size $n$ goes
to infinity, so there is a $\frac{1}{\sqrt{n}}$ scale in the result of $(i)$.
In the second part of the theorem, note that
$\bm{P}_{\bm{A}}=\bm{I}_{p}-\bm{M}_{\bm{A}}$, where $\bm{M}_{\bm{A}}$ is the
projection matrix onto the orthogonal complement of the linear space spanned
by the columns of $\bm{A}$. Thus, $\bm{P}_{\widehat{\bm{A}}}$ and
$\bm{P}_{\bm{A}^{0}}$ represent the spaces spanned by $\widehat{\bm{A}}$ and
$\bm{A}^{0}$, and we show that they are asymptotically the same in $(ii)$.

Next, we obtain the rate of convergence.

######  Theorem 2.

Under Assumptions 1 \- 6, if $p/n\rightarrow\rho>0$,

$\left\|\frac{\left(\bm{C}^{0}-\widehat{\bm{C}}\right)}{\sqrt{n}}\bm{M}_{\bm{F}}\right\|=O_{p}\left(\frac{1}{\sqrt{p}}\right),$

where $\bm{M}_{\bm{F}}$ is defined in (14).

We study the case when the dimension $p$ and the sample size $n$ are
comparable. We achieve rate $\sqrt{p}$ convergence, considering
$\frac{\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|}{\sqrt{n}}$ on the whole. It
is expected that the rate of convergence for smoothing models depends on the
number of discrete points $p$ observed on each curve.

######  Remark 7.

The asymptotic result in Theorem 2 contains a projection matrix
$\bm{M}_{\bm{F}}$. This matrix projects $\bm{C}^{0}-\widehat{\bm{C}}$ onto the
space orthogonal to the factor matrix $\bm{F}$. This theorem shows the
interplay between $\bm{C}^{0}$ and $\bm{F}$. When $\bm{C}^{0}$ and $\bm{F}$
are orthogonal,
$(\bm{C}^{0}-\widehat{\bm{C}})\bm{M}_{\bm{F}}=\bm{C}^{0}-\widehat{\bm{C}}$,
and we obtain the rate of convergence of $\bm{C}^{0}-\widehat{\bm{C}}$. When
$\bm{C}^{0}$ and $\bm{F}$ are not orthogonal, the inference on $\bm{C}^{0}$
will be affected by the existence of the factor model component.

We further begin to establish the limiting distribution. It is shown in
Appendix A that

$\left\|\sqrt{p}\frac{\left(\bm{C}^{0}-\widehat{\bm{C}}\right)}{\sqrt{n}}\bm{M}_{\bm{F}}\right\|=\left\|\left(\frac{1}{p}\Phi^{\top}\bm{M}_{\bm{A}^{0}}\bm{\Phi}\right)^{-1}\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}\right\|+o_{p}(1).$

The limiting distribution is constructed based on the first term on the right-
hand side. Let $\bm{\omega}_{j}$ denote the $j$th column of the $K\times p$
matrix $\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}$. We then have the following
lemma.

######  Lemma 1.

Under Assumptions 1 \- 7, for any vector $\bm{b}=(b_{1},\dots,b_{n})^{\top}$,

$\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}\bm{b}\overset{d}{\to}\mathcal{N}(\bm{0},\bm{L}),$

where $\bm{L}$ is defined in (17).

This lemma paves the way for the next theorem on asymptotic normality.

######  Theorem 3.

Under Assumptions 1 \- 7, if $p/n\rightarrow\rho>0$, we have for any vector
$\bm{b}\in\mathbb{R}^{n}$

$\sqrt{p}\left(\frac{\bm{C}^{0}-\widehat{\bm{C}}}{\sqrt{n}}\right)\bm{M}_{\bm{F}}\bm{b}\overset{d}{\to}\mathcal{N}\left(\bm{0},\bm{Q}\left(\bm{A}^{0}\right)^{-1}\bm{L}\bm{Q}\left(\bm{A}^{0}\right)^{-1}\right),$

where $\bm{M}_{\bm{F}}$ is defined in Theorem 2, $\bm{L}$ is defined in (17),
and

$\bm{Q}\left(\bm{A}^{0}\right)\equiv\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{\Phi}.$

The vector $\bm{b}$ comes from the same vector in Lemma 1. The asymptotic bias
is zero since we assume no serial or cross-sectional correlation in the error
terms. This is a simplified setting, which can be extended to allow for weak
correlations in errors in both directions. In that case, the asymptotic
distribution will include a non-zero bias term.

######  Remark 8.

Theorem 3 shows that the asymptotic distribution of the coefficient matrix
$\widehat{\bm{C}}$ relies on the unobserved factor loading matrix
$\bm{A}^{0}$. Although we are unable to consistently estimate $\bm{A}^{0}$
with $\widehat{\bm{A}}$, what we need is in fact the projected matrix
$\bm{M}_{\bm{A}^{0}}$, which can be estimated by $\bm{M}_{\widehat{\bm{A}}}$.
We are able to find estimators for $\bm{Q}$ and $\bm{L}$ based on
$\bm{M}_{\widehat{\bm{A}}}$,

$\displaystyle\widehat{\bm{Q}}$
$\displaystyle=\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}$
$\displaystyle\widehat{\bm{L}}$
$\displaystyle=\frac{\sigma^{2}}{np}\sum_{i}^{n}\sum_{j}^{p}\widehat{\bm{\omega}}_{j}^{\top}\widehat{\bm{\omega}}_{j}\left(\sum_{k}^{n}\widehat{\psi}_{ik}b_{k}\right)^{2},$

where $\widehat{\bm{\omega}}_{j}$ is the $j$th column of the $K\times p$
matrix $\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}$ and $\widehat{\psi}_{ik}$
is the $(i,k)$th element in the matrix $\bm{M}_{\widehat{\bm{F}}}$.

##  5 Extending the model to nonparametric smoothing

In the FASM, we use basis smoothing where we assume the basis functions
$\phi_{k}(u)$ are known, and we show the asymptotic properties in the previous
section. However, in practice, the basis functions are usually unknown, and
nonparametric smoothing techniques are frequently used. In this section, we
extend the proposed model to a smoothing spline.

In spline smoothing, a spline basis is used to model functions. We consider
smoothing splines, where regularized regression is performed, and the knots
are placed on all the observed discrete points. The most commonly considered
basis is the cubic smoothing splines, where the order is 4. The number of
basis functions equals the number of interior knots plus the order. Thus, we
use $p+2$ spline basis functions. Denote the $k$th basis function as
$\psi_{k}(u)$, $k=1,\dots,p+2$. Let the $p\times(p+2)$ matrix $\bm{\Psi}$
denote the basis matrix, where the $(j,k)$th element is $\psi_{k}(u_{j})$. The
objective function can be written as

$\text{SSR}(\bm{w}_{i},\bm{A},\bm{f})=\sum_{i=1}^{n}\left[(\bm{Y}_{i}-\bm{\Psi}\bm{w}_{i}-\bm{A}\bm{f}_{i})^{\top}(\bm{Y}_{i}-\bm{\Psi}\bm{w}_{i}-\bm{A}\bm{f}_{i})+\alpha\bm{w}_{i}^{\top}\bm{U}\bm{w}_{i}\right],$

where $\bm{w}_{i}$ is the vector of smoothing coefficient and the matrix
$\bm{U}=\int_{\mathcal{I}}D^{2}\bm{\Psi}(s)D^{2}\bm{\Psi}^{\top}(s)ds$ is
similarly defined as the matrix $\bm{R}$ in (5). The estimators are the
solution of the equation system

$\displaystyle\begin{cases}\widehat{\bm{w}}_{i}=\left(\bm{\Psi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Psi}+\alpha\bm{U}^{\top}\right)^{-1}\bm{\Psi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{Y}_{i},\quad
i=1,\dots,n\\\
\left[\frac{1}{np}\sum_{i=1}^{n}(\bm{Y}_{i}-\bm{\Psi}\widehat{\bm{c}}_{i})(\bm{Y}_{i}-\bm{\Psi}\widehat{\bm{w}}_{i})^{\top}\right]\widehat{\bm{A}}=\widehat{\bm{A}}\bm{V}_{np},\end{cases}$

which can be solved by iteration approach in Section 3.2. The function
estimator is
$\widehat{\mathcal{X}}_{i}(u)=\widehat{\bm{w}}^{\top}_{i}\bm{\Psi}(u)$, where
$\bm{\Psi}(u)$ is a vector containing all the $p+2$ basis functions
$\psi_{k}(u)$.

It could be seen that the model is almost the same as the model proposed in
Section 3. The difference lies in the dimensions of the matrices $\bm{\Psi}$
and $\bm{\Phi}$. In the parametric model, we assume the basis functions are
known, and the number of basis functions is fixed, so $\bm{\Phi}$ is a
$p\times K$ matrix. While in smoothing spline modeling, the matrix $\bm{\Psi}$
is of dimension $p\times(p+2)$. The number of basis functions $p+2$ goes to
infinity. This does not render our model estimation infeasible, however,
because we have included a penalty term.

##  6 Statistical inference on covariance matrix estimation

Having presented the model estimation approach and the estimators’ asymptotic
properties, we now consider statistical inference with FASM. Our model serves
as a dimension reduction technique and avoids the curse of dimensionality
problem, rendering making inferences from the model convenient.

Covariance estimation is fundamental in both FDA and high-dimensional data
analysis. In these areas, data are of high dimensions, which brings many
challenges. In the FDA, the number of discrete points on each curve is often
larger than the number of curves. Similarly, the dimension $p$ of high-
dimensional data is typical of the same order or larger than the sample size
$n$. In this case, the traditional sample covariance estimator no longer
works. Dimension reduction by imposing some structure on the data is one of
the main ways to solve this problem (see, e.g., Wong et al., 2003; Bickel &
Levina, 2008; Fan et al., 2008). By reducing the data dimension with a
smoothing model and a factor model in FASM, we propose an alternative
covariance matrix estimator.

We consider the covariance matrix of the observed high-dimensional data
$\bm{Y}_{i}$. Let

$\bm{\Sigma}_{Y}\equiv\text{cov}(\bm{Y}).$

Based on the FASM where

$\bm{Y}_{i}=\bm{\Phi}\bm{c}_{i}+\bm{A}\bm{f}_{i}+\bm{\epsilon}_{i},$

we obtain

$\bm{\Sigma}_{\bm{Y}}=\bm{\Phi}\bm{\Sigma}_{\bm{c}}\bm{\Phi}^{\top}+\bm{A}\bm{\Sigma}_{\bm{f}}\bm{A}^{\top}+\bm{\Sigma}_{\bm{\epsilon}},$
(21)

where $\bm{\Sigma}_{\bm{c}}$ and $\bm{\Sigma}_{\bm{f}}$ are covariance
matrices of the vectors $\bm{c}$ and $\bm{F}$ respectively and
$\bm{\Sigma}_{\epsilon}$ denotes the error variance structure and is a
diagonal matrix under Assumption 4. Based on the above equation, we have an
estimator

$\widehat{\bm{\Sigma}}_{\bm{Y}}=\bm{\Phi}\widehat{\bm{\Sigma}}_{\bm{c}}\bm{\Phi}^{\top}+\widehat{\bm{A}}\widehat{\bm{\Sigma}}_{\bm{f}}\widehat{\bm{A}}^{\top}+\widehat{\bm{\Sigma}}_{\bm{\epsilon}},$
(22)

where $\widehat{\bm{\Sigma}}_{\bm{c}}$ and $\widehat{\bm{\Sigma}}_{\bm{f}}$
can be calculated by

$\displaystyle\widehat{\bm{\Sigma}}_{\bm{c}}$
$\displaystyle=\frac{1}{n-1}\bm{C}\bm{C}^{\top}-\frac{1}{n(n-1)}\bm{C}\bm{1}\bm{1}^{\top}\bm{C}^{\top}$
$\displaystyle\widehat{\bm{\Sigma}}_{\bm{f}}$
$\displaystyle=\frac{1}{n-1}\bm{F}\bm{F}^{\top}-\frac{1}{n(n-1)}\bm{F}\bm{1}\bm{1}^{\top}\bm{F}^{\top},$

where $\bm{1}$s are vectors containing ones, the dimensions of which depend on
the matrices multiplied before and after the vectors. The diagonal error
covariance matrix $\bm{\Sigma}_{\bm{\epsilon}}$ is estimated by

$\widehat{\bm{\Sigma}}_{\bm{\epsilon}}=\text{diag}\left(n^{-1}\widehat{\bm{E}}\widehat{\bm{E}}^{\top}\right),$

where $\widehat{\bm{E}}$ is the residual matrix calculated as
$\widehat{\bm{E}}=\bm{Y}-\bm{\Phi}\widehat{\bm{C}}-\widehat{\bm{A}}\widehat{\bm{F}}^{\top}$.

######  Remark 9.

In functional data analysis where the functional signal is the focus, the
estimation of the covariance function
$\bm{\Phi}\bm{\Sigma}_{\bm{c}}\bm{\Phi}^{\top}$ is of main interest. In this
paper, we study the covariance structure of the mixture of functional data and
high-dimensional data. This type of covariance estimator based on factor
models has also been used in previous literature. For example, Fan et al.
(2008) employed a multi-factor model where the factors are assumed observable,
while Fan et al. (2011) considered an extension to approximate factor models
where cross-sectional correlation is allowed in the error terms.

We compare the finite sample performance using mean squared error (MSE) of the
proposed covariance estimator with the ordinary sample covariance estimator.
When the factor structure is ignored, the sample covariance estimator is
expected to have a larger variance than our estimator. The advantage of the
proposed estimator is shown in Section 7.4.

##  7 Simulation studies

In this section, we use simulated data to illustrate the superiority of the
proposed model. The FASM is compared with the smoothing model in Section 7.1
to 7.3. In Section 7.4, we compare the finite sample performance of the
covariance matrix estimator introduced in Section 6 with the ordinary sample
covariance estimator. In Section 7.6, we show how the FASM performs when
applied to functional data with step jumps.

###  7.1 Data generation

We generate simulated data $Y_{ij}$, where $i=1,\dots,n$ and $j=1,\dots,p$
from the following model:

$\displaystyle Y_{ij}$
$\displaystyle=\mathcal{X}_{i}(u_{j})+\eta_{ij}+\epsilon_{ji}$
$\displaystyle=\sum_{k=1}^{13}c_{ik}\phi_{k}(u_{j})+\sum_{k=1}^{4}\lambda_{ik}F_{kj}+\epsilon_{ji},$

where $\phi_{k}(u)$ are chosen as B-spline basis functions of order 4 and the
smoothing coefficients $c_{ik}$ are generated from $\mathcal{N}(0,1.5^{2})$.
The factors $F_{kj}$ follow $\mathcal{N}(0,0.5^{2})$ and the factor loadings
$(\lambda_{i1},\lambda_{i2},\lambda_{i3},\lambda_{i4})^{\top}\sim\mathcal{N}(\bm{\mu},\bm{\Sigma})$,
where $\bm{\Sigma}$ is a 4 by 4 covariance matrix. We set the multivariate
mean term $\bm{\mu}=\bm{0}$ and variance $\bm{\Sigma}=\sigma^{2}\bm{I}_{4}$.
We adjust the value of $\sigma^{2}$ to control the signal-to-noise ratio. When
$\sigma^{2}$ is large, the signal-to-noise level is low, and when $\sigma^{2}$
is small, the signal-to-noise level is high. The random error terms
$\epsilon_{ji}$ follow $\mathcal{N}(0,0.5^{2})$.

###  7.2 Estimation

The numeric iteration procedure for finding
$(\widehat{\bm{c}}_{i},\widehat{\bm{A}},\widehat{\bm{f}})$ is introduced in
Section 3. We compare the FASM with the smoothing model, where the factor
model component is ignored. The smoothing model can be expressed as:

$\bm{Y}_{i}=\bm{\Phi}\bm{c}_{i}+\bm{\epsilon}_{i},$

where the coefficient estimator is calculated as:

$\widehat{\bm{c}_{i}}=\left(\bm{\Phi}^{\top}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{Y}_{i},\qquad
i=1,\dots,n.$

The tuning parameter $\alpha$ is also chosen using mGCV defined in (11).

###  7.3 Results

We repeat the simulation setup 100 times and obtain the estimated smooth
function
$\widehat{\mathcal{X}}_{i}(u)=\widehat{\bm{c}}_{i}^{\top}\bm{\Phi}(u)$. The
averaged mean squared error (aMSE) for function estimation is calculated as

$\text{aMSE}=\frac{1}{np}\sum_{i=1}^{n}\sum_{j=1}^{p}\left[\mathcal{X}_{i}(u_{j})-\widehat{\mathcal{X}}_{i}(u_{j})\right]^{2}.$

The results are reported in Table LABEL:tab:1. With the same sample size $n$,
increasing the number of points $p$ on the curve decreases the estimation
error. However, with the same value for $p$, increasing the sample size does
not decrease the estimation error. This is consistent with the convergence
rate stated in Section 4, where the estimator converges with the rate related
to $p$. When $\sigma$ is large, such that the signal-to-noise ratio is high,
the FASM performs better than the smoothing model.

Table 1: The aMSE of the function estimates with different samples sizes and dimensions. The adjustment of $\sigma^{2}$ value is used to control the signal-to-noise ratio. |  |  aMSE
---|---|---
Dimension |  $\sigma$ |  FASM |  Smoothing model
$n=20,p=51$ |  $\sigma=0.5$ |  $0.1669$ |  $\bm{0.1177}$
|  $\sigma=0.75$ |  $0.1774$ |  $\bm{0.1729}$
|  $\sigma=1$ |  $\bm{0.2151}$ |  $0.2424$
$n=20,p=101$ |  $\sigma=0.5$ |  $0.0719$ |  $\bm{0.0661}$
|  $\sigma=0.75$ |  $\bm{0.0893}$ |  $0.0979$
|  $\sigma=1$ |  $\bm{0.1163}$ |  $0.1384$
$n=50,p=51$ |  $\sigma=0.5$ |  $0.1311$ |  $\bm{0.1207}$
|  $\sigma=0.75$ |  $\bm{0.1522}$ |  $0.1787$
|  $\sigma=1$ |  $\bm{0.1943}$ |  $0.2518$
$n=100,p=101$ |  $\sigma=0.5$ |  $\bm{0.0593}$ |  $0.0674$
|  $\sigma=0.75$ |  $\bm{0.0794}$ |  $0.0989$
|  $\sigma=1$ |  $\bm{0.1051}$ |  $0.1385$

###  7.4 Covariance matrix estimation

This section shows the finite sample performance of the covariance estimator
defined in (22). We also calculate the regular sample covariance estimator
$\widehat{\bm{\Sigma}}^{*}_{Y}$ using

$\widehat{\bm{\Sigma}}_{\bm{Y}}^{*}=\frac{1}{n-1}\left(\bm{Y}-\mkern
1.5mu\overline{\mkern-1.5mu\bm{Y}\mkern-1.5mu}\mkern
1.5mu\right)\left(\bm{Y}-\mkern
1.5mu\overline{\mkern-1.5mu\bm{Y}\mkern-1.5mu}\mkern 1.5mu\right)^{\top},$

where the $p\times n$ matrix $\mkern
1.5mu\overline{\mkern-1.5mu\bm{Y}\mkern-1.5mu}\mkern 1.5mu$ is the sample mean
matrix whose $j$th row elements are $\frac{1}{n}\sum_{i=1}^{n}Y_{ij}$.

Both estimators are compared with the population covariance matrix, which is
calculated using (21). We calculate the estimation errors under the Frobenius
norm as

$\text{MSE}=\frac{1}{p}\left\|\widehat{\bm{\Sigma}}_{\bm{Y}}-\bm{\Sigma}_{\bm{Y}}\right\|^{2}.$

We show the MSE results in Table LABEL:tab:4. It can be seen that the FASM
produces smaller MSE values in almost all cases.

Table 2: The MSE of the two covariance estimators with different sample sizes and dimensions. The adjustment of $\sigma^{2}$ value is used to control the signal-to-noise ratio. |  |  MSE
---|---|---
Dimension |  $\sigma$ |  FASM |  Sample covariance
$n=20,p=51$ |  $\sigma=0.5$ |  $\bm{0.090}$ |  $0.143$
|  $\sigma=0.75$ |  $\bm{0.128}$ |  $0.198$
|  $\sigma=1$ |  $\bm{0.218}$ |  $0.326$
$n=20,p=101$ |  $\sigma=0.5$ |  $\bm{0.089}$ |  $0.145$
|  $\sigma=0.75$ |  $\bm{0.118}$ |  $0.197$
|  $\sigma=1$ |  $\bm{0.211}$ |  $0.333$
$n=50,p=51$ |  $\sigma=0.5$ |  $\bm{0.041}$ |  $0.058$
|  $\sigma=0.75$ |  $\bm{0.059}$ |  $0.078$
|  $\sigma=1$ |  $\bm{0.117}$ |  $0.122$
$n=100,p=101$ |  $\sigma=0.5$ |  $\bm{0.019}$ |  $0.027$
|  $\sigma=0.75$ |  $\bm{0.031}$ |  $0.038$
|  $\sigma=1$ |  $0.066$ |  $\bm{0.062}$

###  7.5 Nonparametric smoothing model

In this section, we apply the factor-augmented nonparametric smoothing model
introduced in Section 5 to simulated data and compare the results with using
nonparametric smoothing models without the factor component.

We generate simulated data $Y_{ij}$, where $i=1,\dots,n$ and $j=1,\dots,p$
from the following model:

$\displaystyle Y_{ij}$
$\displaystyle=\mathcal{X}_{i}(u_{j})+\eta_{ij}+\epsilon_{ji}$
$\displaystyle=\sum_{k=1}^{9}c_{ik}\phi_{k}(u_{j})+\sum_{k=1}^{4}\lambda_{ik}F_{kj}+\epsilon_{ji},$

where $\phi_{k}(u)$ are Fourier basis functions and the smoothing coefficients
$c_{ik}$ are generated from $\mathcal{N}(0,1.5^{2})$. The factors $F_{kj}$
follow $\mathcal{N}(0,0.5^{2})$ and the factor loadings
$(\lambda_{i1},\lambda_{i2},\lambda_{i3},\lambda_{i4})^{\top}\sim\mathcal{N}(\bm{\mu},\bm{\Sigma})$,
where $\bm{\Sigma}$ is a 4 by 4 covariance matrix. The random error terms
$\epsilon_{ji}$ follow $\mathcal{N}(0,0.5^{2})$. We set the multivariate mean
term $\bm{\mu}=\bm{0}$ and variance $\bm{\Sigma}=\sigma^{2}\bm{I}_{4}$. We
adjust the value of $\sigma^{2}$ to control the signal-to-noise ratio. When
$\sigma^{2}$ is large, the signal-to-noise level is low, and when $\sigma^{2}$
is small, the signal-to-noise level is high.

#### Smoothing spline

We use order 4 B-spline basis with knots at every data point. With data of
dimension $p$, we use $p+2$ basis functions. The tuning parameter $\lambda$ is
selected by the mean generalized cross-validation (11) at each iteration step.
The covariance estimate is calculated using (22). Table LABEL:tab:4 presents
the results.

Table 3: Using smoothing splines: The aMSE of the estimated functions and the MSE of the two covariance estimators with different sample sizes and dimensions. The adjustment of $\sigma^{2}$ value is used to control the signal-to-noise ratio. |  | aMSE | MSE
---|---|---|---
Dimension | $\sigma$ | FASM | Smoothing | FASM | Sample covariance
$n=20,p=51$ | $\sigma=0.5$ | $0.154$ | $\bm{0.142}$ | $\bm{1.535}$ | $1.882$
| $\sigma=0.75$ | $\bm{0.212}$ | $0.219$ | $\bm{1.485}$ | $1.931$
| $\sigma=1$ | $\bm{0.293}$ | $0.317$ | $\bm{1.649}$ | $2.336$
$n=20,p=101$ | $\sigma=0.5$ | $\bm{0.076}$ | $\bm{0.076}$ | $\bm{1.479}$ | $1.814$
| $\sigma=0.75$ | $\bm{0.112}$ | $0.121$ | $\bm{1.593}$ | $2.086$
| $\sigma=1$ | $\bm{0.159}$ | $0.175$ | $\bm{1.587}$ | $2.268$
$n=50,p=51$ | $\sigma=0.5$ | $\bm{0.138}$ | $0.142$ | $\bm{0.612}$ | $0.713$
| $\sigma=0.75$ | $\bm{0.192}$ | $0.213$ | $\bm{0.630}$ | $0.765$
| $\sigma=1$ | $\bm{0.270}$ | $0.313$ | $\bm{0.690}$ | $0.869$
$n=100,p=101$ | $\sigma=0.5$ | $\bm{0.070}$ | $0.075$ | $\bm{0.298}$ | $0.349$
| $\sigma=0.75$ | $\bm{0.104}$ | $0.119$ | $\bm{0.309}$ | $0.377$
| $\sigma=1$ | $\bm{0.153}$ | $0.177$ | $\bm{0.352}$ | $0.429$

###  7.6 Misidentification of the basis function

We elaborate on the example presented in Section 2.2. We generate data from

$Y_{ij}=\sum_{k=1}^{7}c_{ik}\phi_{k}(u_{j})+\epsilon_{ji},\quad i=1,\dots,n,\
j=1,\dots,p,$ (23)

where $\\{\phi_{k}(u),k=1,\dots,7\\}$ are a set of Fourier basis functions.
The first Fourier basis function $\phi_{1}(u)$ is the constant function; the
remainder are sine and cosine pairs with integer multiples of the base period.
We generate the Fourier functions with doubled frequencies in the second half
to simulate the change in the basis functions. In particular, when
$u\in[0,0.5],\ \phi_{k}(u)=2\sin(k\pi u),$ for $k=2,4,6,$ and
$\phi_{k}(u)=2\cos[(k-1)\pi u],$ for $k=3,5,7$, and when $u\in(0.5,1],\
\phi_{k}(u)=2\sin(2k\pi u),$ for $k=2,4,6,$ and $\phi_{k}(u)=2\cos[2(k-1)\pi
u],$ for $k=3,5,7$. The coefficients $c_{ik}$ are generated from the normal
distribution with mean 0 and variance $0.5^{2}$. The error terms are also
drawn from the normal distribution with mean 0 and variance $0.5^{2}$. The
generated $Y_{ij}$ are shown in Figure 2(a). It can be seen that the data
exhibit more variation in the second half of the interval.

Suppose we were unaware of the change in the frequencies of the basis
functions and used the bases in the first half to fit the data on the whole
interval. The smoothing model residuals, shown in Figure 2(b), are large in
the second half. When the frequency of the basis functions is misidentified, a
smoothing model with the wrong set of bases is inadequate. We conduct
principal component analysis on the smoothing residuals; the eigenvalues in
descending order are shown in Figure 4(a). The residuals preserve a spiked
structure, where six common factors can explain most of the variation.

We also apply FASM to the same data with the wrong set of basis functions.
According to the eigenvalue scree plot, we retain six factors in the model
($r=6$). The resulting residuals are shown in Figure 4(b). The large residuals
in the second part of Figure 2 (b) are removed. When the basis functions are
misidentified, the FASM serves as a remedy.

(a) Spikes of the smoothing residuals (b) Residuals of FASM

Figure 4: Applying FASM to the data generated by (23).

###  7.7 Functional data with step jumps

We study the case where the functional data exhibit a dramatic change in the
mean level within a small window. We generate data from the following model

$Y_{ij}=\mu(u_{j})+\sum_{k=1}^{7}c_{ik}\phi_{k}(u_{j})+\epsilon_{ji},\quad
i=1,\dots,n,\quad j=1,\dots,p,$

where the basis functions $\phi_{k}(u)$ are order 4 B-spline bases. The
coefficients $c_{ik}$ come from $\mathcal{N}(0,1.5^{2})$ and the error terms
from $\mathcal{N}(0,0.5^{2})$. The mean function $\mu(u)$ is generated by a
linear combination of 25 B-spline basis functions. Figure 5 shows an example
of the mean function- there is a sharp increase in the mean function at around
$u=0.5$.

The change in the mean level happens at $u=0.5$, and $\delta$ denotes the
amount of change. Figure 3 is generated using $\delta=2$. Figure 6 compares
the residuals from the smoothing model and the FASM. With the smoothing model,
the residuals around the jump are large. In contrast, our model explains the
large residuals around the structural break very well. In the aspect of model
selection, we consider the trade-off between model fit and model flexibility.
We first define a notion of flexibility for a fitted model with the degrees of
freedom. We use the same concept as in most textbooks that the degrees of
freedom measure the number of parameters estimated from the data required to
define the model. The degrees of freedom for the smoothing model are
calculated by (12) of the last step of convergence. The degrees of freedom for
the FASM is

$\text{df}=\text{trace}\left[\bm{\Phi}(\bm{\Phi}^{\top}\bm{M}_{{\widehat{\bm{A}}}^{(t)}}\bm{\Phi}+\alpha\bm{R})^{-1}\bm{\Phi}^{\top}\bm{M}_{{\widehat{\bm{A}}}^{(t)}}\right]+r,$

where $r$ is the number of factors retained in the fitted model, the larger
the degrees of freedom, the more flexible the fitted models are. To quantify
the model fitting, we use

$\text{RMSE}=\sqrt{\frac{1}{np}\sum_{i=1}^{n}\sum_{j=1}^{p}\left(Y_{ij}-\widehat{Y}_{ij}\right)^{2}},$

where
$\widehat{Y}_{ij}=\sum_{k=1}^{K}\widehat{c}_{ik}\phi_{k}(u_{j})+\widehat{\eta}_{ij}$.
In Table 4, we show the simulation results by changing the value of the mean
shift $\delta$. The RMSE of the FASM is always smaller than the compared
model. The degrees of freedom when $\delta=1$ are similar. When $\delta$
increases, the degrees of freedom is smaller for the proposed model.
Therefore, we achieve better fit but less flexibility with the FASM.

Figure 5: The mean function $\mu(u)$.

(a) Residuals from applying smoothing model (b) Residuals from FASM

Figure 6: Residual plots of the two models. Table 4: The trade-off between model fitting and flexibility. | RMSE | DF
---|---|---
| Smoothing | Proposed | Smoothing | Proposed
$\delta=1$ | $0.2045$ | $\bm{0.1631}$ | $\bm{10.68}$ | $11.15$
$\delta=2$ | $0.2063$ | $\bm{0.1640}$ | $17.59$ | $\bm{11.03}$
$\delta=3$ | $0.3308$ | $\bm{0.1647}$ | $14.23$ | $\bm{10.94}$

##  8 Application to climatology

In this section, we apply the FASM to two real data sets. In Section 8.1, we
compare Canadian yearly temperature and precipitation data and demonstrate the
advantages of the FASM when the measurement error is large. In Section 8.2, we
analyze Australian daily temperature data and display the necessity of
including the factor model because of the spike structure of the data.

###  8.1 Canadian weather data

In Section 2.1, we introduced Canadian weather data. Raw observations of daily
temperature and precipitation data are presented in Figure 1. Since the true
basis functions are unknown, we apply the FASM with nonparametric smoothing
splines introduced in Section 5 to these two datasets.

We use order 4 B-spline basis functions with knots at every data point. Thus,
when the number of data points is 365, we use 367 basis functions. The number
of factors $r$ is chosen with the scree plot showing the fraction of variation
explained. For temperature data, we presumed the measurement error is small.
The resulting smoothed curves are shown in Figure 7. Compared with using the
smoothing model introduced in Section 7.2, the FASM generates similar results.
This meets our expectation that our model should work the same as a simple
smoothing model when measurement error does not exist.

In Section 2.1, we suspect large measurement errors are contained in the raw
log precipitation data. We apply the two models to the log precipitation data;
the resulting smoothed curves are presented in Figure 8. The plot on the right
shows smoother curves, especially at the drop in the blue curve (the
’Victoria’ Station) at around day 200. Looking at the residual plots in Figure
9, our model mainly explains some extreme residuals left out from solely
applying the smoothing model. As in Section 6, we also compare the RMSE and
degrees of freedom of the two fitted models; they are 0.1933 and 14.41 for the
smoothing model and 0.1659 and 12.71 respectively for the proposed model.
Thus, in terms of model selection, our model performs better across both model
fit and model simplicity.

(a) Smoothed temperature curves from basis smoothing with the penalty (b)
Smoothed temperature curves from the FASM

Figure 7: Comparison between the smoothed curves.

(a) Smoothed log precipitation curves using the smoothing model (b) Smoothed
log precipitation curves using the FASM

Figure 8: Comparison between the smoothed curves.

(a) Residuals from basis smoothing with the penalty (b) Residuals from the
FASM

Figure 9: Comparison between the residuals.

###  8.2 Australian temperature data

In this section, we consider Friday temperature data at Adelaide airport. We
choose Adelaide because it tends to have the hottest temperature among
Australia’s big cities. Data from other weekdays exhibit similar features and
are not shown here. The data are measured every half an hour from the year
1997 to 2007. The sample size $n$ is 508, and the number of discrete data
points $p$ from each curve is 48. The plot of the raw data can be found in
Figure 10(a). It can be seen that the data are quite noisy, with extreme
values in some of the curves due to large measurement errors.

We use the B-spline basis functions of order 4 with knots at every data point.
A penalized smoothing model is fitted to the data, with the tuning parameter
selected to minimize the mGCV value. The residuals are shown in Figure 10(b).
As can be seen, the smoothing model fails to capture the extreme values in the
data.

We check the residuals’ spikiness in Figure 10(b) by conducting a principal
component analysis. The eigenvalues in descending order are shown in Figure
10(c). The first few eigenvalues are significantly larger than the rest. This
means the residuals contain information captured by just a few factors, which
calls for a further dimension reduction model on the residuals.

As a comparison, the FASM is also applied to the data. The tuning parameter
for the smoothing part is selected based on mGCV at each step of the
iteration. The number of factors retained in the factor model component is
five. The residuals are shown in Figure 10(d). The extreme values are almost
all removed from the remaining residuals.

(a) Raw temperature data (b) Residuals of the smoothing model (c) Eigenvalues
of the residuals in the plot (b) (d) Residuals of the FASM

Figure 10: Half-hourly temperature data on Friday at Adelaide airport.

##  9 Conclusion and future work

In this paper, we propose a factor-augmented smoothing model for functional
data. We study raw functional data, which is a mixture of functional curves
and high-dimensional errors. When measurement error is informative, a
smoothing model alone is inadequate to capture data variation and recover the
signal functional component. The proposed model incorporates a factor
structure into the smoothing model to further explain the large residuals. We
propose a numerical iteration approach to simultaneously obtain estimates in
the smoothing model and the factor model. The asymptotic distribution of the
estimators is given with proof. Our model also serves as a dimension reduction
method on functional and high-dimensional mixture data, easing the path to
making inferences. We provide an example of the construction of a covariance
estimator for the raw data. Further, we show that the model can be applied in
situations where there is misidentification in the data structure, two
examples of which are the misspeficications of smoothing basis functions and
the neglect of the step jumps in the mean level of the functions. The
advantages of the proposed model are demonstrated in extensive simulation
studies. We also show how our model performs via applications to Canadian
weather data and Australian temperature data.

The proposed model is a good start point for modeling complex data structures.
The data we deal with are a mixture of smooth functional curves and high-
dimensional measurement error. The factor model component can be regarded as a
”boosting” component that improves model accuracy. Extending from this idea,
the model can be applied to other data structures. One example is the data
that contain change points. Change point is a popular problem in many
statistics and econometric topics and has been extensively studies in the
multivariate setting. Previous literature on change point in functional data
include Berkes et al. (2009); Hörmann & Kokoszka (2010) and Hörmann et al.
(2015). It is shown in simulation examples that our model can be used for
modeling functional data with change point in the cross-sectional direction.
The model can be modified to account for change point also in the sample
direction. Further research can be conducted along this line.

## References

  * (1)
  * Ahn & Horenstein (2013) Ahn, S. C. & Horenstein, A. R. (2013), ‘Eigenvalue ratio test for the number of factors’, Econometrica 81(3), 1203–1227.
  * Akaike (1974) Akaike, H. (1974), ‘A new look at the statistical model identification’, IEEE Transactions on Automatic Control 19(6), 716–723.
  * Bai & Ng (2002) Bai, J. & Ng, S. (2002), ‘Determining the number of factors in approximate factor models’, Econometrica 70(1), 191–221.
  * Berkes et al. (2009) Berkes, I., Gabrys, R., Horváth, L. & Kokoszka, P. (2009), ‘Detecting changes in the mean of functional observations’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71(5), 927–946.
  * Bickel & Levina (2008) Bickel, P. J. & Levina, E. (2008), ‘Regularized estimation of large covariance matrices’, The Annals of Statistics 36(1), 199–227.
  * Cai & Yuan (2011) Cai, T. T. & Yuan, M. (2011), ‘Optimal estimation of the mean function based on discretely sampled functional data: Phase transition’, The Annals of Statistics 39(5), 2330–2355.
  * Cuevas (2014) Cuevas, A. (2014), ‘A partial overview of the theory of statistics with functional data’, Journal of Statistical Planning and Inference 147, 1–23.
  * Eubank (1999) Eubank, R. L. (1999), Nonparametric Regression and Spline Smoothing, 2nd edn, Marcel Dekker, New York.
  * Fan et al. (2008) Fan, J., Fan, Y. & Lv, J. (2008), ‘High dimensional covariance matrix estimation using a factor model’, Journal of Econometrics 147(1), 186–197.
  * Fan & Gijbels (1996) Fan, J. & Gijbels, I. (1996), Local Polynomial Modelling and Its Applications, Chapman & Hall, London.
  * Fan et al. (2011) Fan, J., Liao, Y. & Mincheva, M. (2011), ‘High dimensional covariance matrix estimation in approximate factor models’, The Annals of Statistics 39(6), 3320.
  * Febrero-Bande et al. (2017) Febrero-Bande, M., Galeano, P. & González-Manteiga, W. (2017), ‘Functional principal component regression and functional partial least-squares regression: An overview and a comparative study’, International Statistical Review 85(1), 61–83.
  * Ferraty & Vieu (2006) Ferraty, F. & Vieu, P. (2006), Nonparametric Functional Data Analysis: Theory and Practice, Springer Science & Business Media, New York.
  * Goia & Vieu (2016) Goia, A. & Vieu, P. (2016), ‘An introduction to recent advances in high/infinite dimensional statistics’, Journal of Multivariate Analysis 146, 1–6.
  * Golub et al. (1979) Golub, G. H., Heath, M. & Wahba, G. (1979), ‘Generalized cross-validation as a method for choosing a good ridge parameter’, Technometrics 21(2), 215–223.
  * Green & Silverman (1999) Green, P. J. & Silverman, B. W. (1999), Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach, Chapman & Hall, London.
  * Hörmann et al. (2015) Hörmann, S., Kidziński, L. & Hallin, M. (2015), ‘Dynamic functional principal components’, Journal of the Royal Statistical Society, Statistical Methodology, Series B 77(2), 319–348.
  * Hörmann & Kokoszka (2010) Hörmann, S. & Kokoszka, P. (2010), ‘Weakly dependent functional data’, The Annals of Statistics 38(3), 1845–1884.
  * Horváth & Kokoszka (2012) Horváth, L. & Kokoszka, P. (2012), Inference for functional data with applications, Vol. 200, Springer Science & Business Media, New York.
  * Jiang et al. (2020) Jiang, B., Yang, Y., Gao, J. & Hsiao, C. (2020), ‘Recursive estimation in large panel data models: Theory and practice’, Journal of Econometrics .
  * Knight & Fu (2000) Knight, K. & Fu, W. (2000), ‘Asymptotics for lasso-type estimators’, The Annals of Statistics 28(5), 1356–1378.
  * Lam et al. (2011) Lam, C., Yao, Q. & Bathia, N. (2011), ‘Estimation of latent factors for high-dimensional time series’, Biometrika 98(4), 901–918.
  * Onatski (2010) Onatski, A. (2010), ‘Determining the number of factors from empirical distribution of eigenvalues’, The Review of Economics and Statistics 92(4), 1004–1016.
  * Onatski (2012) Onatski, A. (2012), ‘Asymptotics of the principal components estimator of large factor models with weakly influential factors’, Journal of Econometrics 168(2), 244–258.
  * Ramsay & Hooker (2017) Ramsay, J. O. & Hooker, G. (2017), Dynamic Data Analysis: Modeling Data with Differential Equations, Springer, New York.
  * Ramsay & Silverman (2002) Ramsay, J. O. & Silverman, B. W. (2002), Applied Functional Data Analysis, Springer, New York.
  * Ramsay & Silverman (2005) Ramsay, J. O. & Silverman, B. W. (2005), Functional Data Analysis, Springer, New York.
  * Reiss et al. (2017) Reiss, P. T., Goldsmith, J., Shang, H. L. & Ogden, R. T. (2017), ‘Methods for scalar-on-function regression’, International Statistical Review 85(2), 228–249.
  * Schwarz (1978) Schwarz, G. (1978), ‘Estimating the dimension of a model’, The Annals of Statistics 6(2), 461–464.
  * Wahba (1990) Wahba, G. (1990), Spline models for observational data, Vol. 59, SIAM.
  * Wand & Jones (1995) Wand, M. P. & Jones, C. M. (1995), Kernel Smoothing, Chapman & Hall, Boca Raton, FL.
  * Wang et al. (2016) Wang, J.-L., Chiou, J.-M. & Müller, H.-G. (2016), ‘Functional data analysis’, Annual Review of Statistics and Its Application 3, 257–295.
  * Wong et al. (2003) Wong, F., Carter, C. K. & Kohn, R. (2003), ‘Efficient estimation of covariance selection models’, Biometrika 90(4), 809–830.
  * Yao & Li (2013) Yao, W. & Li, R. (2013), ‘New local estimation procedure for a non-parametric regression function for longitudinal data’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75(1), 123–138.
  * Zhang & Wang (2016) Zhang, X. & Wang, J. L. (2016), ‘From sparse to dense functional data and beyond’, The Annals of Statistics 44(5), 2281–2321.

Supplement to “Factor-augmented Smoothing Model for Functional Data” Yuan
GaoPostal address: Research School of Finance, Actuarial Studies and
Statistics, Level 4, Building 26C, Kingsley St, Australian National
University, Canberra, ACT 2601, Australia; Email<EMAIL_ADDRESS>
Research School of Finance, Actuarial Studies and Statistics
Australian National University

Han Lin Shang
Department of Actuarial Studies and Business Analytics
Macquarie University

Yanrong Yang
Research School of Finance, Actuarial Studies and Statistics
Australian National University

This section contains the proofs for the theorems in the main article. In
Appendix A, we provide the proofs for the theorems in Section 4. In Appendix
B, we include the results of a proposition and its proof. In Appendix C, the
lemmas used for the proofs in Appendix A and B are stated as well as their
proofs.

### Appendix A

Theorem 2 is the main result of the asymptotic theories, and the proof of it
is lengthy. Thus, we include in the following the outlines for the proof
before we show the details.

#### Outlines for proof of Theorem 2

In Theorem 2, we find the order of convergence of the estimated coefficient
matrix $\widehat{\bm{C}}$. The difference between $\widehat{\bm{C}}$ and
$\bm{C}^{0}$ could be written into three terms:

$\frac{1}{p}\left(\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\alpha\bm{R}\right)(\widehat{\bm{C}}-\bm{C}^{0})=\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}\bm{F}^{\top}+\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}.$
(24)

The term
$\frac{1}{p}\left(\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\alpha\bm{R}\right)$
is $O_{p}(1)$. The first term on the right-hand side of (19) comes from the
penalty, and the order can be found easily from Assumption 6. The third term
contains the random error matrix $\bm{E}$, and the order can be found using
the result in Lemma 10. The second term is the most complicated one, and we
show in the following proof that it could be further broken down into eight
terms. We find the order of each of the eight terms using the lemmas in
Appendix C. Most of the terms can be shown to be
$o_{p}\left(\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)$ and thus can be
omitted. Combining the remaining terms, we arrive at the result

$\displaystyle\left(\widehat{\bm{C}}-\bm{C}^{0}\right)\bm{M}_{\bm{F}}=$
$\displaystyle\bm{Q}^{-1}\left(\widehat{\bm{A}}\right)\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\bm{Q}^{-1}\left(\widehat{\bm{A}}\right)\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}\bm{M}_{\bm{F}}$
$\displaystyle+O_{p}\left(\frac{1}{\min(n,p)}\right)+O_{p}\left(\frac{\sqrt{n}}{p\sqrt{p}}\right)+O_{p}\left(\frac{1}{\sqrt{np}}\right),$
(25)

where matrix $\bm{Q}$ and $\bm{M}_{\bm{F}}$ are

$\bm{Q}\left(\bm{\widehat{\bm{A}}}\right)=\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}\qquad\bm{M}_{\bm{F}}=\bm{I}_{n}-\bm{F}\left(\bm{F}^{\top}\bm{F}\right)^{-1}\bm{F}^{\top}.$

The first term on the right-hand side of (Outlines for proof of Theorem 2) is
$O_{p}(1)$ using the assumption on the tuning parameter $\alpha$. We also show
the second term is $O_{p}(1)$ using results from the lemmas. When $n$ and $p$
are of the same order, we are able to show the norm of
$\left(\widehat{\bm{C}}-\bm{C}^{0}\right)$ projected on the matrix
$\bm{M}_{\bm{F}}$ is $O_{p}(1)$.

Next begins the formal proofs.

#### Proof of Theorem 1

###### Proof.

The concentrated objective function defined in Section 4.2 is

$S_{np}(\bm{c}_{i},\bm{A})=\frac{1}{np}\sum_{i=1}^{n}\left[(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i})^{\top}\bm{M}_{\bm{A}}(\bm{Y}_{i}-\bm{\Phi}\bm{c}_{i})+\alpha\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i}\right]-\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}^{0}}\bm{\epsilon}_{i}.$

Assume $\bm{c}_{i}^{0}=\bm{0}$ for simplicity without loss of generality. From
$\bm{Y_{i}}=\bm{\Phi}\bm{c}_{i}^{0}+\bm{A}^{0}\bm{f}_{i}+\bm{\epsilon}_{i}=\bm{A}^{0}\bm{f}_{i}+\bm{\epsilon}_{i}$,
we have

$\displaystyle S_{np}(\bm{c}_{i},\bm{A})=$
$\displaystyle\frac{1}{np}\sum_{i=1}^{n}\left[(\bm{A}^{0}\bm{f}_{i}+\bm{\epsilon}_{i}-\bm{\Phi}\bm{c}_{i})^{\top}\bm{M}_{\bm{A}}(\bm{A}^{0}\bm{f}_{i}+\bm{\epsilon}_{i}-\bm{\Phi}\bm{c}_{i})+\alpha\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i}\right]-\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}^{0}}\bm{\epsilon}_{i}$
$\displaystyle=$
$\displaystyle\frac{1}{np}\sum_{i=1}^{n}\bm{f}_{i}^{\top}{\bm{A}^{0}}^{\top}\bm{M}_{\bm{A}}\bm{A}^{0}\bm{f}_{i}+\frac{1}{np}\sum_{i=1}^{n}\bm{c}_{i}^{\top}{\bm{\Phi}}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}-\frac{2}{np}\sum_{i=1}^{n}\bm{f}_{i}^{\top}{\bm{A}^{0}}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}$
$\displaystyle+\frac{2}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{A}^{0}\bm{f}_{i}-\frac{2}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}+\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}(\bm{M}_{\bm{A}}-\bm{M}_{\bm{A}^{0}})\bm{\epsilon}_{i}+\frac{\alpha}{np}\sum_{i=1}^{n}\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i}.$

Denote the first three terms in the above equation as

$\widetilde{S}_{np}(\bm{c}_{i},\bm{A})=\frac{1}{np}\sum_{i=1}^{n}\bm{f}_{i}^{\top}{\bm{A}^{0}}^{\top}\bm{M}_{\bm{A}}\bm{A}^{0}\bm{f}_{i}+\frac{1}{np}\sum_{i=1}^{n}\bm{c}_{i}^{\top}{\bm{\Phi}}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}-\frac{2}{np}\sum_{i=1}^{n}\bm{f}_{i}^{\top}{\bm{A}^{0}}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}.$

Then by Lemma 4,

$S_{np}(\bm{c}_{i},\bm{A})=\widetilde{S}_{np}(\bm{c}_{i},\bm{A})+o_{p}(1).$

It is easy to see that
$\widetilde{S}_{np}(\bm{c}_{i}^{0}=\bm{0},\bm{A}^{0}\bm{H})=0$ for any
$r\times r$ invertible $\bm{H}$, because
$\bm{M}_{\bm{A}^{0}\bm{H}}=\bm{M}_{\bm{A}^{0}}$ and
$\bm{M}_{\bm{A}^{0}}\bm{A}^{0}=\bm{0}$.

Here we define two matrix operations before further transformations on
$\tilde{S}_{np}(\bm{c}_{i},\bm{A})$. For an $m\times n$ matrix $\bm{U}$ and a
$p\times q$ matrix $\bm{V}$, the vectorization of $\bm{U}$ is defined as

$\text{vec}(\bm{U})\equiv(u_{1,1},\dots,u_{m,1},u_{1,2},\dots,u_{m,2},u_{1,n},\dots,u_{m,n})^{\top},$

and the Kronecker product $\bm{U}\otimes\bm{V}$ is the $pm\times qn$ block
matrix defined as

$\displaystyle\bm{U}\otimes\bm{V}\equiv\begin{bmatrix}u_{1,1}\bm{V}&\dots&u_{1,n}\bm{V}\\\
\vdots&&\vdots\\\ u_{m,1}\bm{V}&\dots&u_{m,n}\bm{V},\end{bmatrix}$

where $u_{ij}$ represents the element on the $i$th row and $j$th column of
matrix $\bm{U}$.

Next we can further write $\tilde{S}_{n,p}(\bm{c}_{i},\bm{A})$ as

$\displaystyle\widetilde{S}_{np}(\bm{c}_{i},\bm{A})$
$\displaystyle=\text{vec}(\bm{M}_{\bm{A}}\bm{A}^{0})^{\top}\left(\frac{\bm{F}^{\top}\bm{F}}{np}\otimes\bm{I}_{p}\right)\text{vec}(\bm{M}_{\bm{A}}\bm{A}^{0})+\frac{1}{n}\sum_{i=1}^{n}\bm{c}_{i}^{\top}\left(\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\right)\bm{c}_{i}$
$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}2\bm{c}_{i}^{\top}\left(\frac{1}{p}\bm{f}_{i}\otimes\bm{M}_{\bm{A}}\bm{\Phi}\right)\text{vec}(\bm{M}_{\bm{A}}\bm{A}^{0}).$

If we denote

$\displaystyle\bm{P}$
$\displaystyle=\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\bm{A}}\bm{\Phi},$
$\displaystyle\bm{W}$
$\displaystyle=\frac{\bm{F}^{\top}\bm{F}}{np}\otimes\bm{I}_{p},$
$\displaystyle\bm{V}_{i}$
$\displaystyle=\frac{1}{p}\bm{f}_{i}\otimes\bm{M}_{\bm{A}}\bm{\Phi},$

and $\bm{\gamma}=\text{vec}(\bm{M}_{\bm{A}}\bm{A}^{0})$, then we can write

$\displaystyle\widetilde{S}_{np}(\bm{c}_{i},\bm{A})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left[\bm{c}_{i}^{\top}\bm{P}\bm{c}_{i}+\bm{\gamma}^{\top}\bm{W}\bm{\gamma}-2\bm{c}_{i}^{\top}\bm{V}_{i}^{\top}\bm{\gamma}\right]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left[\bm{c}_{i}^{\top}\left(\bm{P}-\bm{V}_{i}^{\top}\bm{W}^{-1}\bm{V}_{i}\right)\bm{c}_{i}+(\bm{\gamma}^{\top}-\bm{c}_{i}^{\top}\bm{V}_{i}^{\top}\bm{W}^{-1})\bm{W}(\bm{\gamma}^{\top}-\bm{W}^{-1}\bm{V}_{i}\bm{c}_{i})\right]$
$\displaystyle\equiv\frac{1}{n}\sum_{i=1}^{n}\left[\bm{c}_{i}^{\top}\bm{D}_{i}(\bm{A})\bm{c}_{i}+\bm{\theta}_{i}^{\top}\bm{W}\bm{\theta}_{i}\right].$

In the last equation,

$\displaystyle\bm{D}_{i}(\bm{A})$
$\displaystyle\equiv\bm{P}-\bm{V}_{i}^{\top}\bm{W}^{-1}\bm{V}_{i}$
$\displaystyle=\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{f}_{i}^{\top}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\bm{f}_{i},$
$\displaystyle\bm{\theta}_{i}$
$\displaystyle\equiv\bm{\gamma}^{\top}-\bm{W}^{-1}\bm{V}_{i}\bm{c}_{i}.$

By Assumption 2 and 3, the matrices $\bm{D}_{i}$ and $\bm{W}$ are positive
definite for each $i$. Thus we have $\widetilde{S}_{np}(\bm{c}_{i},\bm{A})\geq
0$. In addition, if either $\bm{c}_{i}\neq\bm{c}_{i}^{0}$ or
$\bm{A}\neq\bm{A}^{0}\bm{H}$, then $\widetilde{S}_{np}(\bm{c}_{i},\bm{A})>0$.
Thus $\widetilde{S}_{np}(\bm{c}_{i},\bm{A})$ achieves its unique minimum at
$(\bm{c}_{i}^{0},\bm{A}^{0})$. Thus we have

$\widehat{\bm{c}}_{i}-\bm{c}_{i}^{0}=o_{p}(1),\qquad i=1,\dots,n.$

Next, we show $\widehat{\bm{c}}_{i}$ is consistent uniformly in $i$.

We can write

$\displaystyle
S_{np}(\bm{c}_{i},\bm{A})-\widetilde{S}_{np}(\bm{c}_{i},\bm{A})=$
$\displaystyle\frac{2}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{A}^{0}\bm{f}_{i}-\frac{2}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}$
$\displaystyle+\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}(\bm{M}_{\bm{A}}-\bm{M}_{\bm{A}^{0}})\bm{\epsilon}_{i}+\frac{\alpha}{np}\sum_{i=1}^{n}\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i}.$

Using Taylor’s expansion at $\bm{c}_{i}$,

$\displaystyle
S_{np}(\bm{c}_{i},\bm{A})-\widetilde{S}_{np}(\bm{c}_{i},\bm{A})=$
$\displaystyle\frac{2}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{A}^{0}\bm{f}_{i}-\frac{2}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}^{0}$
$\displaystyle+\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}(\bm{M}_{\bm{A}}-\bm{M}_{\bm{A}^{0}})\bm{\epsilon}_{i}+\frac{\alpha}{np}\sum_{i=1}^{n}{\bm{c}_{i}^{0}}^{\top}\bm{R}\bm{c}_{i}^{0}$
$\displaystyle+\left(-\frac{2}{np}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}+\frac{2\alpha}{np}{\bm{c}_{i}^{0}}^{\top}\bm{R}\right)(\bm{c}_{i}-\bm{c}_{i}^{0})+\Delta,$

where $\Delta$ denotes the small order terms. Then we have

$\displaystyle\left(-\frac{2}{np}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}+\frac{2\alpha}{np}{\bm{c}_{i}^{0}}^{\top}\bm{R}\right)\left(\bm{c}_{i}-\bm{c}_{i}^{0}\right)=$
$\displaystyle
S_{np}(\bm{c}_{i},\bm{A})-\widetilde{S}_{np}(\bm{c}_{i},\bm{A})-\frac{2}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{A}^{0}\bm{f}_{i}$
$\displaystyle+\frac{2}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}^{0}+\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}(\bm{M}_{\bm{A}}-\bm{M}_{\bm{A}^{0}})\bm{\epsilon}_{i}$
$\displaystyle-\frac{\alpha}{np}\sum_{i=1}^{n}{\bm{c}_{i}^{0}}^{\top}\bm{R}\bm{c}_{i}^{0}+\Delta$
(26)

In the above equation, the right-hand side is $o_{p}(1)$ uniformly in $i$.
This is
because$S_{np}(\bm{c}_{i},\bm{A})-\widetilde{S}_{np}(\bm{c}_{i},\bm{A})=o_{p}(1)$
and $S_{np}(\bm{c}_{i},\bm{A})$ and $\widetilde{S}_{np}(\bm{c}_{i},\bm{A})$
both consist of summations over $i$. Furthermore, all other terms on the
right-hand side are $o_{p}(1)$ as proved in Lemma 4 and all contain summations
over $i$. On the left-hand side of (Proof.),
$-\frac{2}{np}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}$ is $o_{p}(1)$
uniformly because
$\mathbb{E}\left(\frac{1}{np}\left\|\bm{\epsilon}^{\top}_{i}\bm{M}_{\bm{A}}\bm{c}_{i}\right\|\right)=o(1)$
as shown in Lemma 4 $(ii)$. Moreover, the term
$\frac{2\alpha}{np}{\bm{c}_{i}^{0}}^{\top}\bm{R}$ is also $o_{p}(1)$ uniformly
using Lemma 4 $(iv)$ and that we assume $\bm{c}_{i}$ are bounded uniformly in
Assumption 2. This leads us to the result that

$\widehat{\bm{c}}_{i}-\bm{c}_{i}^{0}=o_{p}(1),\quad\text{uniformly for all
}i=1,\dots,n$

Combining the $i$, we have

$\frac{\left\|\widehat{\bm{C}}-\bm{C}^{0}\right\|}{\sqrt{n}}=o_{p}(1).$

To prove part $(ii)$, note that the centred objective function satisfies
$S_{np}(\bm{c}_{i}^{0}=\bm{0},\bm{A}^{0})=0$ and, by definition in (20), we
have $S_{np}(\widehat{\bm{c}}_{i},\widehat{\bm{A}})\leq 0$. Therefore,

$0\geq
S_{np}(\widehat{\bm{c}}_{i},\widehat{\bm{A}})=\widetilde{S}_{np}(\widehat{\bm{c}}_{i},\widehat{\bm{A}})+o_{p}(1).$

Combined with $\widetilde{S}_{np}(\widehat{\bm{c}}_{i},\widehat{\bm{A}})\geq
0$, it must be true that

$\widetilde{S}_{np}(\widehat{\bm{c}}_{i},\widehat{\bm{A}})=o_{p}(1).$

This implies that

$\frac{1}{np}\sum_{i=1}^{n}\bm{F}_{i}^{\top}{\bm{A}^{0}}^{\top}\bm{M}_{\bm{A}}\bm{A}^{0}\bm{F}_{i}=\text{tr}\left[\frac{{\bm{A}^{0}}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}}{p}\frac{\bm{F}^{\top}\bm{F}}{n}\right]=o_{p}(1).$

Since $\bm{F}^{\top}\bm{F}/n=O_{p}(1)$, it must be true that

$\frac{{\bm{A}^{0}}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}}{p}=\frac{{\bm{A}^{0}}^{\top}\bm{A}^{0}}{p}-\frac{{\bm{A}^{0}}^{\top}\widehat{\bm{A}}}{p}\frac{\widehat{\bm{A}}^{\top}\bm{A}^{0}}{p}=o_{p}(1).$
(27)

By Assumption 4, ${\bm{A}^{0}}^{\top}\bm{A}^{0}/p$ is invertible. Thus
${\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p$ is also invertible. Next,

$\|\bm{P}_{\widehat{\bm{A}}}-\bm{P}_{\bm{A}^{0}}\|^{2}=\text{tr}[(\bm{P}_{\widehat{\bm{A}}}-\bm{P}_{\bm{A}^{0}})^{2}]=2\text{tr}(\bm{I}_{r}-\widehat{\bm{A}}^{\top}\bm{P}_{\bm{A}^{0}}\widehat{\bm{A}}/p).$

But (27) implies
$\widehat{\bm{A}}^{\top}\bm{P}_{\bm{A}^{0}}\widehat{\bm{A}}/p\rightarrow\bm{I}_{r}$,
which means $\|\bm{P}_{\widehat{\bm{A}}}-\bm{P}_{\bm{A}}\|\rightarrow 0$. ∎

#### Proof of Theorem 2

###### Proof.

Writing the first equation in (10) in matrix notation, we have

$\widehat{\bm{C}}=\left(\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{Y}.$
(28)

Substitute $\bm{Y}=\bm{\Phi}\bm{C}^{0}+\bm{A}^{0}\bm{f}^{\top}+\bm{E}$ into
(28) and subtract the matrix $\bm{C}^{0}$ on both sides, we get

$\displaystyle\widehat{\bm{C}}-\bm{C}^{0}=\left[\left(\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}-\bm{I}_{K}\right]\bm{C}^{0}$
$\displaystyle+\left(\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}\bm{F}^{\top}$
$\displaystyle+\left(\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\alpha\bm{R}\right)^{-1}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E},$

or

$\frac{1}{p}\left(\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\alpha\bm{R}\right)\left(\widehat{\bm{C}}-\bm{C}^{0}\right)=\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}\bm{F}^{\top}+\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}$
(29)

We first look at the second term on the right-hand side of (29). Recall that
$\bm{M}_{\widehat{\bm{A}}}=\bm{I}_{p}-\widehat{\bm{A}}\widehat{\bm{A}}^{\top}/p$.
We have $\bm{M}_{\widehat{\bm{A}}}\widehat{\bm{A}}=\bm{0}$. Thus

$\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}=\bm{M}_{\widehat{\bm{A}}}\left(\bm{A}^{0}-\widehat{\bm{A}}\bm{H}^{-1}+\widehat{\bm{A}}\bm{H}^{-1}\right)=\bm{M}_{\widehat{\bm{A}}}\left(\bm{A}^{0}-\widehat{\bm{A}}\bm{H}^{-1}\right),$

where $\bm{H}$ is defined in (43). Using (47), it follows that

$\displaystyle\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}\bm{F}^{\top}$
$\displaystyle=-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\left(I1+\dots+I8\right)\left(\frac{{\bm{A}^{0}}^{\top}\widehat{\bm{A}}}{p}\right)^{-1}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\bm{F}^{\top}$
(30) $\displaystyle\equiv J1+\dots+J8.$

In the following, we calculate the order for each from $J1$ to $J8$. Note that
$I1$ to $I8$ are defined in (Proof.). Before we begin, for simplicity, denote

$\bm{G}\equiv\left(\frac{{\bm{A}^{0}}^{\top}\widehat{\bm{A}}}{p}\right)^{-1}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}.$
(31)

We prove in Lemma 5 that $\bm{G}=O_{p}(1)$. We also use the fact that
$\left\|\bm{M}_{\widehat{\bm{A}}}\right\|=O_{p}(1)$. Now

$J1=-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\left(I1\right)\bm{G}\bm{F}^{\top}.$
(32)

Since
$I1=O_{p}\left(\frac{\sqrt{p}}{n}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|^{2}\right)$,
using the result from Lemma 2 $(i)$, the term $J1$ is bounded in norm by
$O_{p}\left(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|^{2}\right)$.
Thus it is also
$o_{p}\left(\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)$.

$\displaystyle J2$
$\displaystyle=-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\left(I2\right)\left(\frac{{\bm{A}^{0}}^{\top}\widehat{\bm{A}}}{p}\right)^{-1}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\bm{F}^{\top}$
$\displaystyle=\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}\left(\widehat{\bm{C}}-\bm{C}^{0}\right)\bm{F}\left(\bm{F}^{\top}\bm{F}\right)^{-1}\bm{F}^{\top}.$
(33)

For the term $J2$, since it is not a small order term, we keep it as what it
is.

Now consider

$\displaystyle J3$
$\displaystyle=-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\left(I3\right)\left(\frac{{\bm{A}^{0}}^{\top}\widehat{\bm{A}}}{p}\right)^{-1}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\bm{F}^{\top}$
$\displaystyle=\frac{1}{np^{2}}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}\left(\widehat{\bm{C}}-\bm{C}^{0}\right)\bm{E}^{\top}\widehat{\bm{A}}\bm{G}\bm{F}^{\top}.$
(34)

We take
$\bm{E}^{\top}\widehat{\bm{A}}=\bm{E}^{\top}\left(\widehat{\bm{A}}-\bm{A}^{0}\bm{H}\right)+\bm{E}^{\top}\bm{A}^{0}\bm{H}$,
where the order of each term can be found in Lemma 6 $(i)$. Again using the
result of Lemma 2 $(i)$ and $(iii)$, it can be shown that $J3$ is
$o_{p}\left(\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)$.

Next

$\displaystyle\|J4\|$
$\displaystyle=\left\|-\frac{1}{p}\bm{\Phi}\bm{M}_{\widehat{\bm{A}}}I4\left(\frac{{\bm{A}^{0}}^{\top}\widehat{\bm{A}}}{p}\right)^{-1}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\bm{F}^{\top}\right\|$
$\displaystyle=O_{p}\left(\frac{\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}}{\sqrt{p}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right).$
(35)

Using Proposition 1, we have
$\frac{1}{\sqrt{p}}\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}=\bm{M}_{\widehat{\bm{A}}}\frac{1}{\sqrt{p}}\left(\bm{A}^{0}-\widehat{\bm{A}}\bm{H}^{-1}\right)=o_{p}\left(1\right)$,
Thus, $\|J4\|=o_{p}\left(\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)$.

It can also be proven that
$\|J5\|=o_{p}\left(\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)$.

Then we consider

$\displaystyle J6$
$\displaystyle=-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}I6\bm{G}\bm{F}^{\top}$
$\displaystyle=-\frac{1}{np^{2}}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{A}^{0}\bm{F}^{\top}\bm{E}^{\top}\widehat{\bm{A}}\bm{G}\bm{F}^{\top}$
$\displaystyle=-\frac{1}{np^{2}}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\left(\bm{A}^{0}-\widehat{\bm{A}}\bm{H}^{-1}\right)\bm{F}^{\top}\bm{E}^{\top}\widehat{\bm{A}}\bm{G}\bm{F}^{\top},$

where the last equation comes from
$\bm{M}_{\widehat{\bm{A}}}\widehat{\bm{A}}\bm{H}^{-1}=\bm{0}$. Now

$\displaystyle\left\|\bm{F}^{\top}\bm{E}^{\top}\widehat{\bm{A}}\right\|$
$\displaystyle=\left\|\bm{E}^{\top}\widehat{\bm{A}}\right\|$
$\displaystyle\leq\left\|\
\bm{E}^{\top}\left(\widehat{\bm{A}}-\bm{A}^{0}\bm{H}\right)\right\|+\left\|\bm{E}^{\top}\bm{A}^{0}\bm{H}\right\|$
$\displaystyle=O_{p}\left(\frac{p}{\min\left(\sqrt{n},\sqrt{p}\right)}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\sqrt{n}\right)+O_{p}\left(\frac{p}{\sqrt{n}}\right)+O_{p}\left(\sqrt{np}\right)$
$\displaystyle=O_{p}\left(\frac{p}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{p}{\sqrt{n}}\right)+O_{p}\left(\sqrt{np}\right),$

using Lemma 6. Thus,

$\displaystyle\|J6\|$
$\displaystyle\leq\left\|\frac{1}{np^{2}}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\left(\bm{A}^{0}-\widehat{\bm{A}}\bm{H}^{-1}\right)\right\|\left\|\bm{F}^{\top}\bm{E}^{\top}\widehat{\bm{A}}\right\|\left\|\bm{G}\right\|\left\|\bm{F}^{\top}\right\|$
$\displaystyle=-\frac{1}{np^{2}}\times\left[O_{p}\left(\frac{p}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{p}{\min\left(n,p\right)}\right)\right]$
$\displaystyle\times\left[O_{p}\left(\frac{p}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{p}{\sqrt{n}}\right)+O_{p}\left(\sqrt{np}\right)\right]\times
O_{p}\left(\sqrt{n}\right)$
$\displaystyle=o_{p}\left(\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{1}{n^{2}}\right)+O_{p}\left(\frac{1}{n\sqrt{p}}\right)+O_{p}\left(\frac{1}{np}\right)+O_{p}\left(\frac{1}{p\sqrt{p}}\right),$
(36)

where Proposition 1 is used in the first equation, and the second equation is
a result of the calculation on the orders. Next

$\displaystyle J7$
$\displaystyle=-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}I7\bm{G}\bm{F}^{\top}$
$\displaystyle=-\frac{1}{np}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}\bm{F}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\bm{F}^{\top}.$
(37)

This term is not a small order term, so we keep it as what it is. And lastly,
the proof of order for the term $J8$ is too long, so we show in Lemma 11 that

$J8=o_{p}\left(\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{1}{\min\left(n,p\right)}\right)+O_{p}\left(\frac{\sqrt{n}}{\sqrt{p}}\frac{1}{\min\left(n,p\right)}\right).$
(38)

Collecting terms from $J1$ to $J8$, we can write (29) as

$\left(\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+\frac{1}{p}\alpha\bm{R}\right)\left(\widehat{\bm{C}}-\bm{C}^{0}\right)=\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+J1+\dots+J8+\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}.$

Combining the results we have found for $J1,J3,J4,J5,J6$ and $J8$ in Proof.,
Proof., Proof. and 38,

$\displaystyle\left(\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+o_{p}(1)\right)\left(\widehat{\bm{C}}-\bm{C}^{0}\right)-J2=$
$\displaystyle\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}+J7$
$\displaystyle+O_{p}\left(\frac{1}{\min(n,p)}\right)+O_{p}\left(\frac{\sqrt{n}}{p\sqrt{p}}\right)+O_{p}\left(\frac{1}{\sqrt{np}}\right).$
(39)

Substitute $J2$ and $J7$ from Proof. and Proof. into (Proof.), we have

$\displaystyle\left(\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}+o_{p}(1)\right)\left(\widehat{\bm{C}}-\bm{C}^{0}\right)-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}\left(\widehat{\bm{C}}-\bm{C}^{0}\right)\bm{F}\left(\bm{F}^{\top}\bm{F}\right)^{-1}\bm{F}^{\top}$
$\displaystyle=$
$\displaystyle\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}-\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}\bm{F}\left(\bm{F}^{\top}\bm{F}\right)^{-1}\bm{F}^{\top}$
$\displaystyle+O_{p}\left(\frac{1}{\min(n,p)}\right)+O_{p}\left(\frac{\sqrt{n}}{p\sqrt{p}}\right)+O_{p}\left(\frac{1}{\sqrt{np}}\right).$
(40)

We combine the two terms on the left-hand side of (Proof.) and also combine
the second and third term on the right-hand side of (Proof.), then we get

$\displaystyle\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}(\widehat{\bm{C}}-\bm{C}^{0})\left(\bm{I}_{n}-\bm{F}\left(\bm{F}^{\top}\bm{F}\right)^{-1}\bm{F}^{\top}\right)=$
$\displaystyle\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}\left(\bm{I}_{n}-\bm{F}(\bm{F}^{\top}\bm{F})^{-1}\bm{F}^{\top}\right)$
$\displaystyle+O_{p}\left(\frac{1}{\min(n,p)}\right)+O_{p}\left(\frac{\sqrt{n}}{p\sqrt{p}}\right)+O_{p}\left(\frac{1}{\sqrt{np}}\right).$

Let
$\bm{Q}(\widehat{\bm{A}})\equiv\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{\Phi}$,
and
$\bm{M}_{\bm{F}}\equiv\bm{I}_{n}-\bm{F}(\bm{F}^{\top}\bm{F})^{-1}\bm{F}^{\top}$.
Left multiplying $\bm{Q}^{-1}(\widehat{\bm{A}})$ to both sides of the equation
above, we have

$\displaystyle(\widehat{\bm{C}}-\bm{C}^{0})\bm{M}_{\bm{F}}=$
$\displaystyle\bm{Q}(\widehat{\bm{A}})^{-1}\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\bm{Q}(\widehat{\bm{A}})^{-1}\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}\bm{M}_{\bm{F}}$
$\displaystyle+O_{p}\left(\frac{1}{\min(n,p)}\right)+O_{p}\left(\frac{\sqrt{n}}{p\sqrt{p}}\right)+O_{p}\left(\frac{1}{\sqrt{np}}\right)$
$\displaystyle=$
$\displaystyle\bm{Q}^{-1}(\bm{A}^{0})\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\bm{Q}^{-1}(\bm{A}^{0})\frac{1}{p}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}$
$\displaystyle+O_{p}\left(\frac{1}{\min(n,p)}\right)+O_{p}\left(\frac{\sqrt{n}}{p\sqrt{p}}\right)+O_{p}\left(\frac{1}{\sqrt{np}}\right),$

where in the last equation, we substitute $\bm{Q}(\widehat{\bm{A}})$ with
$\bm{Q}(\bm{A}^{0})$ using Lemma 8 and substitute
$\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\widehat{\bm{A}}}\bm{E}$ with
$\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}$ using Lemma 10.
Note that
$\sqrt{p}O_{p}\left(\frac{\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|^{2}}{n}\right)$
in the result of Lemma 10 is dominated by
$\sqrt{p}\frac{\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|}{\sqrt{n}}$. Next by
multiplying a scale of $\frac{\sqrt{p}}{\sqrt{n}}$,

$\displaystyle\frac{\sqrt{p}}{\sqrt{n}}(\widehat{\bm{C}}-\bm{C}^{0})\bm{M}_{\bm{F}}=$
$\displaystyle\bm{Q}(\bm{A}^{0})^{-1}\frac{1}{p}\alpha\bm{R}\bm{C}^{0}+\bm{Q}(\bm{A}^{0})^{-1}\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}$
$\displaystyle+\frac{\sqrt{p}}{\sqrt{n}}\times
O_{p}\left(\frac{1}{\min(n,p)}\right)+\frac{\sqrt{p}}{\sqrt{n}}\times
O_{p}\left(\frac{\sqrt{n}}{p\sqrt{p}}\right)+O_{p}\left(\frac{1}{\sqrt{np}}\right)$
$\displaystyle=$ $\displaystyle O_{p}(1),$ (41)

when $n$ and $p$ are of the same order, that is $p/n\rightarrow\rho>0$. ∎

#### Proof of Theorem 3

From (Proof.), we have when $p/n\rightarrow\rho>0$,

$\frac{\sqrt{p}}{\sqrt{n}}(\widehat{\bm{C}}-\bm{C}^{0})\bm{M}_{\bm{F}}=\bm{Q}(\bm{A}^{0})^{-1}\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}+o_{p}(1).$

Using Lemma 1 we have, for any vector $\bm{b}=(b_{1},\dots,b_{n})^{\top}$,

$\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}\bm{b}\overset{d}{\to}\mathcal{N}(0,\bm{L}),$
(42)

where $\bm{L}$ is defined in (17).

Multiplying the constant matrix $\bm{Q}(\bm{A}^{0})^{-1}$ to (42), we have the
result

$\bm{Q}(\bm{A}^{0})^{-1}\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{M}_{\bm{F}}\bm{b}\overset{d}{\to}\mathcal{N}\left(0,\bm{Q}(\bm{A}^{0})^{-1}\bm{L}\bm{Q}(\bm{A}^{0})^{-1}\right).$

The theorem is thus proved.

### Appendix B

In this section, we provide the proposition used in Appendix A, along with its
proof.

######  Proposition 1.

Under Assumptions 1 to 4, we have the following statements:

  1. (i)

The matrix $\bm{V}_{np}$ defined in (9) is invertible and
$\bm{V}_{np}\overset{p}{\to}\bm{V}$, where the $r\times r$ matrix $\bm{V}$ is
a diagonal matrix consisting of the eigenvalues of
$\bm{\Sigma}_{\bm{F}}\bm{\Sigma}_{\bm{A}}$;

  2. (ii)

Let

$\bm{H}=(\bm{F}^{\top}\bm{F}/n)^{-1}({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)^{-1}\bm{V}_{np},$
(43)

then $\bm{H}$ is $r\times r$ invertible matrix and

$\frac{1}{p}\|\widehat{\bm{A}}-\bm{A}^{0}\bm{H}\|^{2}=O_{p}\left(\frac{1}{n}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|^{2}\right)+O_{p}\left(\frac{1}{\min(n,p)}\right).$

###### Proof.

Write the second equation in (10) in a matrix form, we have

$\frac{1}{np}(\bm{Y}-\bm{\Phi}\widehat{\bm{C}})(\bm{Y}-\bm{\Phi}\widehat{\bm{C}})^{\top}\widehat{\bm{A}}=\widehat{\bm{A}}\bm{V}_{np}.$

By (2), we also have

$\bm{Y}-\bm{\Phi}\widehat{\bm{C}}=\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})+\bm{A}^{0}\bm{F}^{\top}+\bm{E}.$
(44)

Plugging it in (44) and by expanding terms, we obtain

$\displaystyle\widehat{\bm{A}}\bm{V}_{np}=$
$\displaystyle\frac{1}{np}\left[\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})+\bm{A}^{0}\bm{F}^{\top}+\bm{E}\right]\left[\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})+\bm{A}^{0}\bm{F}^{\top}+\bm{E}\right]^{\top}\widehat{\bm{A}}$
$\displaystyle=$
$\displaystyle\frac{1}{np}\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})(\bm{C}^{0}-\widehat{\bm{C}})^{\top}\bm{\Phi}^{\top}\widehat{\bm{A}}+\frac{1}{np}\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})\bm{F}{\bm{A}^{0}}^{\top}\widehat{\bm{A}}$
$\displaystyle+\frac{1}{np}\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})\bm{E}^{\top}\widehat{\bm{A}}+\frac{1}{np}\bm{A}^{0}\bm{F}^{\top}(\bm{C}^{0}-\widehat{\bm{C}})^{\top}\bm{\Phi}^{\top}\widehat{\bm{A}},$
$\displaystyle+\frac{1}{np}\bm{E}(\bm{C}^{0}-\widehat{\bm{C}})^{\top}\bm{\Phi}^{\top}\widehat{\bm{A}}+\frac{1}{np}\bm{A}^{0}\bm{F}^{\top}\bm{E}^{\top}\widehat{\bm{A}}$
$\displaystyle+\frac{1}{np}\bm{E}\bm{F}{\bm{A}^{0}}^{\top}\widehat{\bm{A}}+\frac{1}{np}\bm{EE}^{\top}\widehat{\bm{A}}$
$\displaystyle+\frac{1}{np}\bm{A}^{0}\bm{F}^{\top}\bm{F}{\bm{A}^{0}}^{\top}\widehat{\bm{A}}$
$\displaystyle\equiv$ $\displaystyle I1+\dots+I9.$ (45)

The above can be rewritten as

$\widehat{\bm{A}}\bm{V}_{np}-\bm{A}^{0}(\bm{F}^{\top}\bm{F}/n)({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)=I1+\dots+I8.$
(46)

Right multiplying
$(\bm{F}^{\top}\bm{F}/n)^{-1}({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)^{-1}$ on
each side, we obtain

$\widehat{\bm{A}}\left[\bm{V}_{np}({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)^{-1}(\bm{F}^{\top}\bm{F}/n)^{-1}\right]-\bm{A}^{0}=(I1+\dots+I8)({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)^{-1}(\bm{F}^{\top}\bm{F}/n)^{-1}.$
(47)

Note that the matrix in the square brackets is $\bm{H}^{-1}$, but the
invertibility of $\bm{V}_{np}$ hasn’t been proved yet. We can write

$\frac{1}{\sqrt{p}}\left\|\widehat{\bm{A}}\left[\bm{V}_{np}({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)^{-1}(\bm{F}^{\top}\bm{F}/n)^{-1}\right]-\bm{A}^{0}\right\|\leq\frac{1}{\sqrt{p}}(\|I1\|+\dots+\|I8\|)\|\bm{G}\|,$
(48)

where $\bm{G}$ is defined in (31) and $\|\bm{G}\|$ is proved to be $O_{p}(1)$
in Lemma 5. In the following, we find the order for each term on the right-
hand side of (48). We repeatedly use results from Lemma 2, where the orders of
the matrices $\bm{\bm{\Phi}},\bm{A}$ and $\bm{F}$ are given. The first term

$\displaystyle\frac{1}{\sqrt{p}}\|I1\|$
$\displaystyle\leq\frac{1}{\sqrt{p}}\frac{1}{np}\|\bm{\Phi}\|\|(\bm{C}^{0}-\widehat{\bm{C}})(\bm{C}^{0}-\widehat{\bm{C}})^{\top}\|\|\bm{\Phi}^{\top}\|\|\widehat{\bm{A}}\|$
$\displaystyle=O_{p}\left(\frac{1}{n}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|^{2}\right)=o_{p}\left(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right).$

For the second term

$\displaystyle\frac{1}{\sqrt{p}}\|I2\|$
$\displaystyle\leq\frac{1}{\sqrt{p}}\frac{1}{np}\|\bm{\Phi}\|\|(\bm{C}^{0}-\widehat{\bm{C}})\|\|\bm{F}\|\|{\bm{A}^{0}}^{\top}\|\|\widehat{\bm{A}}\|$
$\displaystyle=O_{p}\left(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right).$

The terms $I3$ to $I5$ are all
$O_{p}(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|)$. The
proofs are similar to the proof for $I2$ since they are only a switch in the
order of the matrices. For the sixth term

$\frac{1}{\sqrt{p}}\|I6\|\leq\frac{1}{\sqrt{p}}\frac{1}{np}\|\bm{A}^{0}\|\|\bm{F}^{\top}\bm{E}^{\top}\|\|\widehat{\bm{A}}\|=O_{p}\left(\frac{1}{\sqrt{n}}\right),$

by Lemma 3 $(i)$. Similarly, for the next term

$\frac{1}{\sqrt{p}}\|I7\|\leq\frac{1}{\sqrt{p}}\frac{1}{np}\|\bm{E}\bm{F}\|\|{\bm{A}^{0}}^{\top}\|\|\widehat{\bm{A}}\|=O_{p}\left(\frac{1}{\sqrt{n}}\right).$

For the last term

$\frac{1}{\sqrt{p}}\|I8\|\leq\frac{1}{\sqrt{p}}\frac{1}{np}\|\bm{EE}^{\top}\|\|\widehat{\bm{A}}\|=O_{p}\left(\frac{1}{\sqrt{n}}\right)+\left(\frac{1}{\sqrt{p}}\right),$

where Lemma 3 $(iv)$ is used.

Putting all the above together, we have

$\displaystyle\frac{1}{\sqrt{p}}\left\|\widehat{\bm{A}}\left[\bm{V}_{np}({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)^{-1}(\bm{F}^{\top}\bm{F}/n)^{-1}\right]-\bm{A}^{0}\right\|=$
$\displaystyle
O_{p}\left(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)$
$\displaystyle+O_{p}\left(\frac{1}{\min(\sqrt{n},\sqrt{p})}\right).$ (49)

To show $(i)$, left multiply (46) by $\frac{1}{p}\widehat{\bm{A}}^{\top}$.
Using $\widehat{\bm{A}}^{\top}\widehat{\bm{A}}/p=\bm{I}_{r}$, we have

$\bm{V}_{np}-(\widehat{\bm{A}}^{\top}\bm{A}^{0}/p)(\bm{F}^{\top}\bm{F}/n)({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)=\frac{1}{p}\widehat{\bm{A}}^{\top}(I1+\dots+I8)=o_{p}(1),$

where the last equality is using Lemma 2 $(v)$ and that
$p^{-1/2}(\|I1\|+\dots+\|I8\|)=o_{p}(1)$ from (Proof.). Thus,

$\bm{V}_{np}=(\widehat{\bm{A}}^{\top}\bm{A}^{0}/p)(\bm{F}^{\top}\bm{F}/n)({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)+o_{p}(1).$

We have shown in (27) that $\widehat{\bm{A}}^{\top}\widehat{\bm{A}}^{0}$ is
invertible, thus $\bm{V}_{np}$ is invertible. To obtain the limit of
$\bm{V}_{np}$, left multiply (46) by $\frac{1}{p}{\bm{A}^{0}}^{\top}$ to yield

$({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)\bm{V}_{np}-({\bm{A}^{0}}^{\top}\bm{A}^{0}/p)(\bm{F}^{\top}\bm{F}/n)({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)=o_{p}(1),$

or

$({\bm{A}^{0}}^{\top}\bm{A}^{0}/p)(\bm{F}^{\top}\bm{F}/n)({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)+o_{p}(1)=({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)\bm{V}_{np}$
(50)

because $p^{-1}{\bm{A}^{0}}^{\top}(\|I1\|+\dots\|I8\|)=o_{p}(1)$. Equation
(50) shows that the columns of $({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p)$ are
the eigenvectors of the matrix
$({\bm{A}^{0}}^{\top}\bm{A}^{0}/p)(\bm{F}^{\top}\bm{F}/n)$, and that
$\bm{V}_{np}$ consists of the eigenvalues of the same matrix in the limit.
Thus, $\bm{V}_{np}\overset{p}{\to}\bm{V}$, where the $r\times r$ matrix
$\bm{V}$ is a diagonal matrix consisting of the eigenvalues of
$\bm{\Sigma}_{\bm{F}}\bm{\Sigma}_{\bm{A}}$.

For $(ii)$, since $\bm{V}_{np}$ is invertible, $\bm{H}$ is also invertible we
can write (Proof.) as

$\frac{1}{\sqrt{p}}\left\|\widehat{\bm{A}}\bm{H}^{-1}-\bm{A}^{0}\right\|=O_{p}\left(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{1}{\min(\sqrt{n},\sqrt{p})}\right).$

By right multiplying the matrix $\bm{H}$, we obtain $(ii)$.

∎

### Appendix C

In this section, we state all the lemmas used for previous theorems and
propositions, along with the proofs of the lemmas.

Lemma 1 is stated in Section 4.2. We provide the proof here.

###### Proof.

For any vector $\bm{b}=(b_{1},\dots,b_{n})^{\top}$,

$\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{b}=\frac{1}{\sqrt{np}}\sum_{i}^{n}\sum_{j}^{p}\bm{\omega}_{j}\epsilon_{ji}b_{i}\equiv\frac{1}{\sqrt{np}}\sum_{i}^{n}\sum_{j}^{p}\bm{x}_{ij},$

where $\bm{\omega}_{j}$ is the $j$th column in the matrix
$\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}$. Since we assume $\epsilon_{ji}$ are
i.i.d., the variance of the above quantity is given by

$\text{var}\left(\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{b}\right)=\text{var}\left(\frac{1}{\sqrt{np}}\sum_{i}^{n}\sum_{j}^{p}\bm{x}_{ij}\right)=\frac{1}{np}\sum_{i}^{p}\sum_{j}^{n}b_{j}b_{i}\sigma^{2}E\left(\bm{\omega}_{i}\bm{\omega}_{j}^{\top}\right).$

The Lindeberg condition is assumed to hold in Assumption 7. Thus we have a
central limit theorem result

$\frac{1}{\sqrt{np}}\bm{\Phi}^{\top}\bm{M}_{\bm{A}^{0}}\bm{E}\bm{b}=\frac{1}{\sqrt{np}}\sum_{i}^{n}\sum_{j}^{p}\bm{x}_{ij}\overset{d}{\to}\mathcal{N}(0,\bm{L}),$

where $\bm{L}$ is defined in (17).

∎

######  Lemma 2.

Under Assumptions 1-3, we have

  1. (i)

$\frac{1}{\sqrt{p}}\|\bm{\Phi}\|=O_{p}(1)$

  2. (ii)

$\frac{1}{\sqrt{p}}\|\bm{A}\|=O_{p}(1)$

  3. (iii)

$\frac{1}{\sqrt{n}}\|\bm{F}\|=O_{p}(1)$

  4. (iv)

$\frac{1}{\sqrt{np}}\|\bm{E}\|=O_{p}(1)$

  5. (v)

$\frac{1}{\sqrt{p}}\|\widehat{\bm{A}}\|=O_{p}(1)$

###### Proof.

In Assumption 1, we assume the basis functions $\phi_{k}(u),\ k=1,\dots,K$ are
bounded. The $p\times K$ basis matrix $\bm{\Phi}$ contains discrete
evaluations on the basis functions so each element is $O_{p}(1)$, thus
$\bm{\Phi}$ is of order $\sqrt{p}$. Similarly, using Assumption 3, we have
results $(ii)$ and $(iii)$. Using Assumption 4, we have result $(iv)$. Lastly,
$(v)$ is directly from the restriction
$\widehat{\bm{A}}^{\top}\widehat{\bm{A}}/p=\bm{I}_{r}$. ∎

######  Lemma 3.

Under Assumptions 1 to 5, we have

  1. (i)

$\frac{1}{np}\|\bm{EF}\|^{2}=O_{p}(1)$

  2. (ii)

$\frac{1}{np}\|\bm{E}^{\top}\bm{\Phi}\|^{2}=O_{p}(1)$ and
$\frac{1}{np}\|\bm{E}^{\top}\bm{A}^{0}\|^{2}=O_{p}(1)$

  3. (iii)

$\frac{1}{np}\|\bm{F}^{\top}\bm{E}^{\top}\bm{A}^{0}\|^{2}=O_{p}(1)$ and
$\frac{1}{np}\|\bm{\Phi}^{\top}\bm{E}^{\top}\bm{F}\|^{2}=O_{p}(1)$

  4. (iv)

$\|\bm{E}^{\top}\bm{E}\|^{2}=O_{p}(n^{2}p)+O_{p}(p^{2}n);\\\
\|\bm{EE}^{\top}\|^{2}=O_{p}(n^{2}p)+O_{p}(p^{2}n);\\\
\|\bm{F}^{\top}\bm{E}^{\top}\bm{E}\|^{2}=O_{p}(n^{2}p)+O_{p}(p^{2}n);\\\
\|\bm{\Phi}^{\top}\bm{E}^{\top}\bm{E}\|^{2}=O_{p}(n^{2}p)+O_{p}(p^{2}n);\\\
\|\bm{\Phi}^{\top}\bm{E}^{\top}\bm{E}\bm{A}^{0}\|^{2}=O_{p}(n^{2}p)+O_{p}(p^{2}n);\\\
\|\bm{F}^{\top}\bm{E}^{\top}\bm{E}\bm{F}\|^{2}=O_{p}(n^{2}p)+O_{p}(p^{2}n)$.

###### Proof.

For $(i)$

$\displaystyle\mathbb{E}\left(\frac{1}{np}\|\bm{EF}\|^{2}\right)$
$\displaystyle=\mathbb{E}\left(\frac{1}{np}\sum_{k=1}^{p}\sum_{i=1}^{n}\sum_{j=1}^{n}\epsilon_{ki}\epsilon_{kj}\bm{f}_{i}^{\top}\bm{f}_{j}\right)$
$\displaystyle=\frac{1}{np}\sum_{k=1}^{p}\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbb{E}(\epsilon_{ki}\epsilon_{kj})\mathbb{E}(\bm{f}_{i}^{\top}\bm{f}_{j})=O(1),$

where the second equation uses the independence between $\epsilon_{ki}$ and
$\bm{f}_{j}$ assumed in Assumption 5.

The proof of $(ii)$ and $(iii)$ is similar to $(i$). For $(iv)$,

$\displaystyle\mathbb{E}\left(\|\bm{E}^{\top}\bm{E}\|^{2}\right)$
$\displaystyle=\mathbb{E}\left(\sum_{ij}^{n}\sum_{kl}^{p}\epsilon_{kj}\epsilon_{lj}\epsilon_{ki}\epsilon_{li}\right)$
$\displaystyle=\sum_{i\neq
j}^{n}\sum_{k=l}^{p}\mathbb{E}(\epsilon^{2}_{kj})\mathbb{E}(\epsilon_{ki}^{2})+\sum_{i=j}^{n}\sum_{k\neq
l}^{p}\mathbb{E}(\epsilon^{2}_{kj})\mathbb{E}(\epsilon_{lj}^{2})+\sum_{i=j}^{n}\sum_{k=l}^{p}\mathbb{E}(\epsilon_{kj}^{4})$
$\displaystyle=O(n^{2}p)+O(p^{2}n)+O(np)$ $\displaystyle=O(n^{2}p)+O(p^{2}n),$

where Assumption 4 is used. The proof of $\|\bm{EE}^{\top}\|$ is the same. The
orders of $\|\bm{F}^{\top}\bm{E}^{\top}\bm{E}\|$ and
$\|\bm{F}^{\top}\bm{E}^{\top}\bm{E}\bm{F}\|^{2}$ are the same since

$\mathbb{E}\left(\|\bm{F}^{\top}\bm{E}^{\top}\bm{E}\|^{2}\right)=\mathbb{E}\left(\sum_{ij}^{n}\sum_{kl}^{p}\epsilon_{kj}\epsilon_{lj}\epsilon_{ki}\epsilon_{li}\|\bm{f}_{i}\|^{2}\right),$

and

$\mathbb{E}\left(\|\bm{F}^{\top}\bm{E}^{\top}\bm{E}\bm{F}\|^{2}\right)=\mathbb{E}\left(\sum_{ij}^{n}\sum_{kl}^{p}\epsilon_{kj}\epsilon_{lj}\epsilon_{ki}\epsilon_{li}\|\bm{f}_{i}\|^{4}\right),$

where the order of $\bm{f}_{i}$ is assumed to be $O_{p}(1)$ in Assumption 3.

∎

######  Lemma 4.

Under Assumptions 1-6,

  1. (i)

$\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{A}^{0}\bm{F}_{i}=o_{p}(1)$

  2. (ii)

$\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}=o_{p}(1)$

  3. (iii)

$\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}(\bm{M}_{\bm{A}}-\bm{M}_{\bm{A}^{0}})\bm{\epsilon}_{i}=o_{p}(1)$

  4. (iv)

$\frac{\alpha}{np}\sum_{i=1}^{n}\bm{c}_{i}^{\top}\bm{R}\bm{c}_{i}=o_{p}(1)$

###### Proof.

We prove $(ii)$. First, we have

$\mathbb{E}\left(\left\|\sum_{i=1}^{n}\bm{\epsilon}_{i}\right\|^{2}\right)=\mathbb{E}\left(\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{p}\epsilon_{ik}\epsilon_{jk}\right)=\sum_{i=j}^{n}\sum_{k=1}^{p}\mathbb{E}\left(\epsilon_{ik}^{2}\right)=O(np).$

Since $\bm{M}_{\bm{A}}=\bm{I}_{p}-\bm{A}\bm{A}^{\top}/p$, we have

$\displaystyle\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_{i}=\frac{1}{np}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{\Phi}\bm{c}_{i}-\frac{1}{np^{2}}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\top}\bm{A}{\bm{A}}^{\top}\bm{\Phi}\bm{c}_{i}.$
(51)

The first term on the right of (51) is $o_{p}(1)$ since

$\displaystyle\mathbb{E}\left(\left\|\sum_{i=1}^{n}\bm{\epsilon}^{\top}_{i}\bm{\Phi}\bm{c}_{i}\right\|^{2}\right)$
$\displaystyle=\mathbb{E}\left(\left\|\sum_{j=1}^{p}\sum_{i=1}^{n}\epsilon_{ji}\bm{\phi}_{j}^{\top}\bm{c}_{i}\right\|^{2}\right)$
$\displaystyle=\mathbb{E}\left(\sum_{t}^{p}\sum_{s}^{n}\sum_{j}^{p}\sum_{i}^{n}\epsilon_{ji}\epsilon_{ts}\bm{c}_{i}^{\top}\bm{\phi}_{j}\bm{\phi}_{t}^{\top}\bm{c}_{s}\right)$
$\displaystyle=\sum_{t}^{p}\sum_{s}^{n}\sum_{j}^{p}\sum_{i}^{n}\mathbb{E}(\epsilon_{ji}\epsilon_{ts})\mathbb{E}(\bm{c}_{i}^{\top}\bm{\phi}_{j}\bm{\phi}_{t}^{\top}\bm{c}_{s})$
$\displaystyle=\sum_{j}^{p}\sum_{i}^{n}\sigma^{2}\mathbb{E}(\bm{c}_{i}^{\top}\bm{\phi}_{j}\bm{\phi}_{j}^{\top}\bm{c}_{i})$
$\displaystyle=O(np),$

where the third equation uses Assumption 5; the fourth equation uses the
assumption that $\epsilon_{ji}$ are independent in both directions.

The second term on the right-hand side of (51) is also $o_{p}(1)$ since

$\displaystyle\mathbb{E}\left(\left\|\sum_{i=1}^{n}\bm{\epsilon}^{\top}_{i}\bm{AA}^{\top}\bm{\Phi}\bm{c}_{i}\right\|^{2}\right)$
$\displaystyle=\mathbb{E}\left(\left\|\sum_{j=1}^{p}\sum_{i=1}^{n}\epsilon_{ji}\bm{a}_{j}^{\top}\bm{A}^{\top}\bm{\Phi}\bm{c}_{i}\right\|^{2}\right)$
$\displaystyle=\mathbb{E}\left(\sum_{t}^{p}\sum_{s}^{n}\sum_{j}^{p}\sum_{i}^{n}\epsilon_{ji}\epsilon_{ts}\bm{c}_{i}^{\top}\bm{\Phi}^{\top}\bm{Aa}_{j}\bm{a}_{t}^{\top}\bm{A}^{\top}\bm{\Phi}\bm{c}_{s}\right)$
$\displaystyle=\sum_{t}^{p}\sum_{s}^{n}\sum_{j}^{p}\sum_{i}^{n}\mathbb{E}(\epsilon_{ji}\epsilon_{ts})\mathbb{E}(\bm{c}_{i}^{\top}\bm{\Phi}^{\top}\bm{Aa}_{j}\bm{a}_{t}^{\top}\bm{A}^{\top}\bm{\Phi}\bm{c}_{s})$
$\displaystyle=\sum_{j}^{p}\sum_{i}^{n}\sigma^{2}\bm{c}_{i}^{\top}\mathbb{E}(\bm{\Phi}_{i}^{\top}\bm{Aa}_{j}\bm{a}_{j}^{\top}\bm{A}^{\top}\bm{\Phi})\bm{c}_{i}$
$\displaystyle=O(np^{3}),$

where the third equality uses the independence in Assumption 5 and the last
equality uses the results in Lemma 2, where $\bm{\Phi}$ and $\bm{A}$ are both
$O_{p}(\sqrt{p})$.

The proofs for $(i)$ and $(iii)$ are similar. And $(iv)$ is a direct result
from Assumption 6. ∎

######  Lemma 5.

Under Assumptions 1-5 , we have

$\bm{G}\equiv\left({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p\right)^{-1}\left(\bm{F}^{\top}\bm{F}/n\right)^{-1}=O_{p}(1).$

###### Proof.

The matrix $\bm{F}^{\top}\bm{F}/n$ is positive definite by Assumption 3. We
have shown in the proof of Theorem 1 in (27) that the matrix
${\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p$ is invertible, thus is also positive
definite. Therefore,
$\lambda_{\min}\left({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p\right)>0$, and
$\lambda_{\min}\left(\bm{F}^{\top}\bm{F}/n\right)>0$, where
$\lambda_{\min}(\cdot)$ denotes the smallest eigenvalue of a matrix. So we
have

$\left({\bm{A}^{0}}^{\top}\widehat{\bm{A}}/p\right)^{-1}=O_{p}(1),\quad\left(\bm{F}^{\top}\bm{F}/n\right)^{-1}=O_{p}(1).$

∎

######  Lemma 6.

We have the following

  1. (i)

$\left\|\bm{E}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})\right\|=O_{p}\left(\frac{p}{\min(\sqrt{n},\sqrt{p})}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}(\sqrt{n})+O_{p}\left(\frac{p}{\sqrt{n}}\right).$

  2. (ii)

$\left\|\bm{F}^{\top}\bm{E}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})\right\|=O_{p}\left(\frac{p}{\min(\sqrt{n},\sqrt{p})}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}(\sqrt{n})+O_{p}\left(\frac{p}{\sqrt{n}}\right).$

###### Proof.

For $(i)$, from Proposition 1, we can write

$\displaystyle\|\bm{E}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})\|$
$\displaystyle=\|\bm{E}^{\top}(I1+\dots,I8)\bm{G}\|$
$\displaystyle\leq\|\bm{E}^{\top}I1\bm{G}\|+\dots+\|\bm{E}^{\top}I8\bm{G}\|$
$\displaystyle=\|a1\|+\dots+\|a8\|.$

To find the order for each term, the results from Lemma 2 are repeatedly used
where the order of the matrices $\bm{\Phi},\bm{A},\bm{\widehat{\bm{A}}}$ and
$\bm{F}$ are given.

$\displaystyle\|a1\|$
$\displaystyle=\left\|\bm{E}^{\top}\frac{1}{np}\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})(\bm{C}^{0}-\widehat{\bm{C}})^{\top}\bm{\Phi}^{\top}\widehat{\bm{A}}\bm{G}\right\|$
$\displaystyle\leq\frac{1}{np}\|\bm{E}^{\top}\bm{\Phi}\|\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|^{2}\|\bm{\Phi}\|\|\widehat{\bm{A}}\|\|\bm{G}\|$
$\displaystyle=O_{p}\left(\frac{\sqrt{p}}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|^{2}\right)=o_{p}\left(\sqrt{p}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right),$

where the order of $\|\bm{E}^{\top}\bm{\Phi}\|$ is from Lemma 3 $(ii)$. The
orders of $\|\bm{\Phi}\|$, $\|\widehat{\bm{A}}\|$ and $\|\bm{G}\|$can be found
from Lemmas 2 and 5. Similarly,

$\displaystyle\|a2\|$
$\displaystyle=\left\|\bm{E}^{\top}\frac{1}{n}\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})\bm{F}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\right\|$
$\displaystyle\leq\frac{1}{n}\|\bm{E}^{\top}\bm{\Phi}\|\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\|\bm{F}\|\left\|\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\right\|$
$\displaystyle=O_{p}\left(\sqrt{p}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right).$
$\displaystyle\|a3\|$
$\displaystyle=\left\|\bm{E}^{\top}\frac{1}{np}\bm{\Phi}(\bm{C}^{0}-\widehat{\bm{C}})\bm{E}^{\top}\widehat{\bm{A}}\bm{G}\right\|$
$\displaystyle\leq\frac{1}{np}\|\bm{E}^{\top}\bm{\Phi}\|\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\|\bm{E}^{\top}\|\|\widehat{\bm{A}}\|\|\bm{G}\|$
$\displaystyle=O_{p}\left(\sqrt{p}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right).$
$\displaystyle\|a4\|$
$\displaystyle=\left\|\bm{E}^{\top}\frac{1}{np}\bm{A}^{0}\bm{F}^{\top}(\bm{C}^{0}-\widehat{\bm{C}})^{\top}\bm{\Phi}^{\top}\widehat{\bm{A}}\bm{G}\right\|$
$\displaystyle\leq\frac{1}{np}\|\bm{E}^{\top}\bm{A}^{0}\|\|\bm{F}^{\top}\|\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\|\bm{\Phi}\|\|\widehat{\bm{A}}\|\|\bm{G}\|$
$\displaystyle=O_{p}\left(\sqrt{p}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right),$

where Lemma 3 $(ii)$ is used.

$\displaystyle\|a5\|$
$\displaystyle=\|\bm{E}^{\top}\frac{1}{np}\bm{E}(\bm{C}^{0}-\widehat{\bm{C}})^{\top}\bm{\Phi}^{\top}\widehat{\bm{A}}\bm{G}\|$
$\displaystyle\leq\frac{1}{np}\|\bm{E}^{\top}\bm{E}\|\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\|\bm{\Phi}\|\|\widehat{\bm{A}}\|\|\bm{G}\|$
$\displaystyle=O_{p}\left(\sqrt{p}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{p}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right),$

where Lemma 3 $(iv)$ is used.

$\displaystyle\|a6\|$
$\displaystyle=\|\bm{E}^{\top}\frac{1}{np}\bm{A}^{0}\bm{F}^{\top}\bm{E}^{\top}\widehat{\bm{A}}\bm{G}\|$
$\displaystyle\leq\frac{1}{np}\|\bm{E}^{\top}\bm{A}^{0}\bm{F}^{\top}\bm{E}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})\bm{G}\|+\frac{1}{np}\|\bm{E}^{\top}\bm{A}^{0}\bm{F}^{\top}\bm{E}^{\top}\bm{A}^{0}\bm{HG}\|$
$\displaystyle\leq\frac{1}{np}\|\bm{E}^{\top}\bm{A}^{0}\|\|\bm{F}^{\top}\bm{E}^{\top}\|\|\widehat{\bm{A}}-\bm{A}^{0}\bm{H}\|\|\bm{G}\|+\frac{1}{np}\|\bm{E}^{\top}\bm{A}^{0}\|\|\bm{F}^{\top}\bm{E}^{\top}\bm{A}^{0}\|\|\bm{HG}\|$
$\displaystyle=O_{p}\left(\sqrt{p}(\frac{\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|}{\sqrt{n}}\right)+O_{p}\left(\frac{1}{\text{min}(\sqrt{n},\sqrt{p})}\right)+O_{p}(1)$
$\displaystyle=O_{p}\left(\sqrt{p}(\frac{\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|}{\sqrt{n}}\right)+O_{p}\left(\frac{1}{\text{min}(\sqrt{n},\sqrt{p})}\right),$

where the order of $\widehat{\bm{A}}-\bm{A}^{0}\bm{H}$ is proved in
Proposition 1 and the order of other matrix norms can be found in Lemma 3
$(i),(ii)$ and $(iii)$.

$\displaystyle\|a7\|$
$\displaystyle=\|\bm{E}^{\top}\frac{1}{np}\bm{EF}{\bm{A}^{0}}^{\top}\widehat{\bm{A}}\bm{G}\|$
$\displaystyle=\left\|\frac{1}{n}\bm{E}^{\top}\bm{EF}\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\right\|$
$\displaystyle\leq\frac{1}{n}\left\|\bm{E}^{\top}\bm{EF}\right\|\left\|\left(\frac{\bm{F}^{\top}\bm{F}}{n}\right)^{-1}\right\|$
$\displaystyle=O_{p}\left(\sqrt{p}\right)+O_{p}\left(\frac{p}{\sqrt{n}}\right),$

where Lemma 3 $(iv)$ is used.

$\displaystyle\|a8\|=$
$\displaystyle\frac{1}{np}\left\|\bm{E}^{\top}\bm{EE}^{\top}\widehat{\bm{A}}\bm{G}\right\|$
$\displaystyle\leq$
$\displaystyle\frac{1}{np}\left\|\bm{E}^{\top}\bm{EE}^{\top}\bm{A}^{0}\bm{HG}\right\|+\frac{1}{np}\left\|\bm{E}^{\top}\bm{EE}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})\bm{G}\right\|$
$\displaystyle\leq$
$\displaystyle\frac{1}{np}\|\bm{E}^{\top}\bm{E}\|\|\bm{E}^{\top}\bm{A}^{0}\|\|\bm{H}\|\|\bm{G}\|+\frac{1}{np}\|\bm{E}^{\top}\bm{E}\|\|\bm{E}^{\top}\|\|\widehat{\bm{A}}-\bm{A}^{0}\bm{H}\|\|\bm{G}\|$
$\displaystyle=$
$\displaystyle\frac{1}{np}\left[O_{p}(n\sqrt{p})+O_{p}(p\sqrt{n})\right]O_{p}(\sqrt{np})$
$\displaystyle+\frac{1}{np}\left[O_{p}(n\sqrt{p})+O_{p}(p\sqrt{n})\right]O_{p}(\sqrt{np})\left[O_{p}\left(\frac{\sqrt{p}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|}{\sqrt{n}}\right)+O_{p}\left(\frac{\sqrt{p}}{\min(\sqrt{n},\sqrt{p})}\right)\right]$
$\displaystyle=$ $\displaystyle
O_{p}(\sqrt{n})+O_{p}(\sqrt{p})+O_{p}\left(\frac{p}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\sqrt{p}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\sqrt{n}\right)+O_{p}\left(\frac{p}{\sqrt{n}}\right)$
$\displaystyle=$ $\displaystyle
O_{p}\left(\frac{p}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\sqrt{p}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\sqrt{n}\right)+O_{p}\left(\frac{p}{\sqrt{n}}\right),$

where the order of $\widehat{\bm{A}}-\bm{A}^{0}\bm{H}$ is proved in
Proposition 1 and the order of other matrix norms can be found in Lemma 3
$(i),(ii)$ and $(iii)$.

Combining all the terms, we have

$\|\bm{E}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})\|=O_{p}\left(\frac{p}{\min(\sqrt{n},\sqrt{p})}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}(\sqrt{n})+O_{p}\left(\frac{p}{\sqrt{n}}\right).$

For $(ii)$, multiplying the matrix $\bm{F}^{\top}$ in the front does not
change the order, using the fact that
$\|\bm{F}^{\top}\bm{E}^{\top}\bm{\Phi}\|$ is of the same order as
$\|\bm{E}^{\top}\bm{\Phi}\|$ and that $\|\bm{F}^{\top}\bm{E}^{\top}\bm{E}\|$
and $\|\bm{F}^{\top}\bm{E}^{\top}\bm{EF}\|$ are of the same order as
$\|\bm{E}^{\top}\bm{E}\|$, as proved in Lemma 3. ∎

######  Lemma 7.

Under Assumptions 1-5, we have the following

  1. (i)

$\frac{1}{p}\bm{\Phi}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})=O_{p}\left(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{1}{\min(n,p)}\right)$

  2. (ii)

$\frac{1}{p}{\bm{A}^{0}}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})=O_{p}\left(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{1}{\min(n,p)}\right)$

  3. (iii)

$\frac{1}{p}\widehat{\bm{A}}^{\top}(\widehat{\bm{A}}-\bm{A}^{0}\bm{H})=O_{p}\left(\frac{1}{\sqrt{n}}\left\|\bm{C}^{0}-\widehat{\bm{C}}\right\|\right)+O_{p}\left(\frac{1}{\min(n,p)}\right)$

  4. (iv)