text:
\section{Introduction}

Functional data analysis (FDA) has received growing attention with the increasing capability to store data 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 \cite{RS02, RS05} and \cite{RH17} provide a comprehensive account of the methodology and applications of the FDA; other relevant monographs include \cite{FV06} and \cite{HK12}. More recent advances in this field can be found in many survey papers \citep[see, e.g.,][]{C14,FGG17, GV16, RGS+17, WCM16}. 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:
\begin{equation*}
Y_{ij} = \mathcal{X}_i(u_j) + \eta_{ij}, \quad i=1,\dots, n,\ j = 1,\dots, p.
\end{equation*}
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 \citep[e.g.][]{WJ95}, local polynomial smoothing \citep[e.g.][]{FG96}, and spline smoothing \citep[e.g.][]{W90, E99, GS94}. With pre-smoothed functions, estimates, such as mean and covariance functions, can be further obtained. More recent studies on functional smoothing approaches include \cite{CY11, YL13}, and \cite{ZW16}. 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
\begin{equation*}
Y_{ij} =  \sum_{k=1}^K c_{ik}\phi_k(u_j) + \eta_{ij}, \quad  i=1,\dots, n,\ j = 1,\dots, p.
\end{equation*}

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 these, factor models are widely used \citep[e.g.][]{FFL08, LYB11}.

We propose using a factor model for the measurement error term. Without further information on the measurement error, the 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.
\begin{equation*}
\eta_{ij} = \bm{a}_j^\top\bm{f}_i + \epsilon_{ij}, \quad i = 1,\dots, n, \  j = 1,\dots, p,
\end{equation*}
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:
\begin{equation*}
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.
\end{equation*}

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~\ref{se:nonparam}, we illustrate the use of spline smoothing approaches. In Section~\ref{se:6.7}, the nonparametric smoothing model is applied to simulated data.

This paper motivates 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~\ref{se:2}.
\begin{enumerate}
\item
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.
\item
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.
\item
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.
\end{enumerate}

Since the latent factors are unobserved, we propose an iterative approach to simultaneously estimate the smooth function and factors. 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~\ref{se:2}. In Section~\ref{se: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~\ref{se:4}.
In Section~\ref{se:6}, we conduct Monte-Carlo simulations on the proposed model under different settings. A real data example is given in Section~\ref{se:7}, and conclusions are drawn in Section~\ref{se:8}.
Due to the page limit, we include some important model extensions and data analysis and the proofs for the theorems in the Supplement.

\section{Motivation}\label{se:2}

We introduce three examples to motivate the proposed model. In these cases, the smoothing model is inadequate to capture the raw data's signal information. In the first example, when a large measurement error exists, the residuals after smoothing are large with 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 smoothing residuals is needed.

\subsection{Functional data with measurement error}\label{se:2.1}

Figure~\ref{fig:2} 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. Although there are some perturbations in the temperature plot, 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~\ref{se:7.1}.

\begin{figure}[!htb]
\center
\includegraphics[width = 8.4cm]{CanadianWeather.pdf}
\includegraphics[width = 8.4cm]{CanadianWeather1.pdf}
\caption{Average daily temperature and log precipitation in 35 Canadian weather stations averaged over the year 1960 to 1994.}\label{fig:2}
\end{figure}

\subsection{Misidentification of the basis function}\label{se:2.2}

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~\ref{fig:9a}. 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 unaware 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~\ref{fig:9b}. 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 generating process and further analysis can be found in Section~\ref{se:6.5}.

\begin{figure}[!htb]
\center
\subcaptionbox{Raw data \label{fig:9a}}
{\includegraphics[width = 8.4cm]{wrong_bas.pdf}}
\subcaptionbox{Residuals\label{fig:9b}}
{\includegraphics[width = 8.4cm]{wrong_bas_res1.pdf}}
\caption{A simulated sample of functional data with changing basis functions.}\label{fig:9}
\end{figure}

\subsection{Functional data with step jumps in the mean level}\label{se:2.3}

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~\ref{fig:1a}. 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~\ref{se:6.6}. The residuals after smoothing are presented in Figure~\ref{fig:1b}. 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~\ref{se:6.6} 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.

\begin{figure}[!htb]
\center
\subcaptionbox{Raw data\label{fig:1a}}
{\includegraphics[width = 8.4cm]{Jump_raw.pdf}}
\subcaptionbox{Residuals\label{fig:1b}}
{\includegraphics[width = 8.4cm]{Jump_res1.pdf}}
\caption{A simulated sample of functional data with step jump.}\label{fig:1}
\end{figure}

\section{Model Specification and Estimation}\label{se:3}

This section formally states the proposed model in Section~\ref{se:3.1} and provides the estimation method in Section~\ref{se: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.

\subsection{Factor-augmented smoothing model}\label{se:3.1}

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.
\begin{equation*}
\mathcal{X}_i(u) = \sum_{k=1}^K c_{ik} \phi_k(u), \quad u \in \mathcal{I},
\end{equation*}
where $\{\phi_k(u), k=1, \dots, K\}$ 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
\begin{align*}
Y_{ij} &= \sum_{k=1}^K c_{ik}\phi_k(u_j) + \eta_{ij},\\
\eta_{ij} &= \bm{a}_j^\top\bm{f}_i  + \epsilon_{ij}, \quad i = 1,\dots, n, \  j = 1,\dots, p,
\end{align*}
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.

\begin{condition}\label{id}
We require
\begin{enumerate}
\item[(i)] $\{\mathcal{X}_i(u_j): i=1, \ldots, n; j=1, \ldots, p\}$ are independent of $\{\eta_{ij}: i=1,\dots, n;\ j=1,\dots, p\}$;
\item[(ii)]
$p^{-1}\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$;

$ n^{-1}\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$.
\end{enumerate}
\end{condition}

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): k=1, 2, \ldots, K\}$ as known, and the number of basis functions $K$ could be either fixed or going to infinity. There are various choices for the basis functions in empirical data analysis, and the decision can be 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.
The number of basis functions $K$ controls the smoothness of the predicted functions. As $K$ increases, the estimator variance increases but the bias decay. Cross-validation, for instance,  can be used for decidng $K$ \citep{WW75}. In this article, our model is estimated with a roughness penalty approach explained in the next section. The smoothness of the functional component is switched from being determined by $K$ to being determined by the tuning parameter of the penalty term.  In practice, we can use a relatively large $K$ and select the tuning parameter carefully. In Remark~\ref{re:4}, we discuss the common methods for selecting the tuning parameter.

\subsection{Penalized estimation approach}\label{se:3.2}

We can write the model for the $i$th object as
\begin{equation}\label{eq:14}
\bm{Y}_i = \bm{\Phi}\bm{c}_i + \bm{A}\bm{f}_i + \bm{\epsilon}_i,\quad i=1,\dots, n
\end{equation}
where
\begin{align*}
\bm{Y}_i = \begin{bmatrix}
Y_{i1}\\ \vdots \\ Y_{ip}
\end{bmatrix},\
\bm{c}_i = \begin{bmatrix}
c_{i1}\\\vdots\\c_{iK}
\end{bmatrix},\
\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}, \
\bm{A} = \begin{bmatrix}
\bm{a}_1^\top\\ \vdots\\ \bm{a}_p^\top
\end{bmatrix}, \
\bm{\epsilon} = \begin{bmatrix}
\epsilon_{i1}\\ \vdots\\ \epsilon_{ip}
\end{bmatrix}.
\end{align*}
Combining all the objects, we have in matrix form
\begin{equation}\label{eq:16}
\bm{Y} = \bm{\Phi}\bm{C} + \bm{A} \bm{F}^\top + \bm{E},
\end{equation}
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}$
\begin{equation}\label{eq:50}
\bm{A}^\top \bm{A}/p = \bm{I}_r, \quad \text{and} \ \bm{F}^\top\bm{F}\ \text{is a diagonal matrix}.
\end{equation}

We propose to implement penalized least squares, where the objective function is defined as
\begin{equation*}
\text{SSR}(\bm{c}_i, \bm{A}, \bm{f}) = \frac{1}{np}\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 \text{PEN}_2(\mathcal{X}_i) \right],
\end{equation*}
where $\text{PEN}_2(\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
\begin{equation*}
\text{PEN}_2(\mathcal{X}_i) = \int_{\mathcal{I}} \left[D^2\mathcal{X}_i(s)\right]^2 ds,
\end{equation*}
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
\begin{equation}\label{eq:2}
\bm{\Phi}(u) = \left[\phi_1(u), \dots,\phi_K(u)\right]^\top.
\end{equation}
Then
\begin{equation*}
\mathcal{X}_i(u) = \bm{c}_i^\top\bm{\Phi}(u).
\end{equation*}
We can re-express the roughness penalty $\text{PEN}_2(\mathcal{X}_i)$ in a matrix form as follows:
\begin{align*}
\text{PEN}_2(\mathcal{X}_i) &= \int_{\mathcal{I}} \left[D^2\mathcal{X}_i(s)\right]^2 ds\\
&= \int_{\mathcal{I}} \left[D^2\bm{c}_i^\top\bm{\Phi}(s)\right]^2 ds\\
&= \int_{\mathcal{I}} \bm{c}_i^\top D^2\bm{\Phi}(s)D^2\bm{\Phi}^\top (s) \bm{c}_i ds\\
&= \bm{c}_i^\top \left[\int_{\mathcal{I}} D^2\bm{\Phi}(s)D^2\bm{\Phi}^\top (s)ds\right] \bm{c}_i\\
&= \bm{c}_i^\top \bm{R} \bm{c}_i, \qquad i = 1,\dots, n
\end{align*}
where
\begin{equation}\label{eq:55}
\bm{R} \equiv \int_{\mathcal{I}} D^2\bm{\Phi}(s)D^2\bm{\Phi}^\top (s)ds.
\end{equation}
The 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$.

\begin{remark}
The number of smoothing coefficient $\bm{c}_i$ increases as the sample size increases. The inclusion of a penalty term penalizes the ``roughness" of the smoothed function and mitigates the effect of the increasing number of parameters to control the model flexibility.
\end{remark}

Thus, the objective function can be written as
\begin{equation*}
\text{SSR}(\bm{c}_i, \bm{A}, \bm{f}) = \frac{1}{np}\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],
\end{equation*}
subject to the constraint $\bm{A}^\top\bm{A}/p = \bm{I}_r$.

We aim to estimate the smoothing coefficient $\bm{c}_i$. We left multiply a matrix to each term in~\eqref{eq:14} to project the factor model term onto a zero matrix. Define the projection matrix
\begin{equation}\label{eq:51}
\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.
\end{equation}
Then
\begin{equation*}
\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}.
\end{equation*}
So we estimate $\bm{c}_i$ from the projected equation
\begin{equation*}
\bm{M}_{\bm{A}}\bm{Y}_i = \bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_i + \bm{M}_{\bm{A}}\bm{\epsilon}_i.
\end{equation*}

The projected objective function becomes
\begin{equation}\label{eq:22}
\text{SSR}(\bm{c}_i, \bm{A}) = \frac{1}{np}\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].
\end{equation}
By taking the derivative of $\text{SSR}(\bm{c}_i, \bm{A})$ with respective to each $\bm{c}_i$ and setting it to zero, we can solve for the estimator $\widehat{\bm{c}}_i$.
\begin{equation*}
\frac{\partial \text{SSR}(\bm{c}_i\bm{A})}{\partial \bm{c}_i} = -\frac{1}{np}(\bm{M}_{\bm{A}}\bm{Y}_i - \bm{M}_{\bm{A}}\bm{\Phi}\bm{c}_i)^\top (\bm{M}_{\bm{A}}\bm{\Phi}) + \frac{1}{np}\alpha \bm{c}_i^\top \bm{R}.
\end{equation*}
Setting the derivative to zero and rearranging the terms, we have
\begin{equation*}
\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.
\end{equation*}
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}$
\begin{equation*}
\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.
\end{equation*}

Next, to estimate $\bm{A}$ and $\bm{f}_i$, we focus on the factor model
\begin{equation*}
\bm{\eta}_i = \bm{A}\bm{f}_i + \bm{\epsilon}_i,
\end{equation*}
and in matrix form
\begin{equation*}
\bm{Z} = \bm{A}\bm{F}^\top + \bm{E},
\end{equation*}
where $\bm{Z} = (\bm{\eta}_1, \dots, \bm{\eta}_n)$. In high-dimensional cases, the unknown factors and loadings are typically estimated by least squares (i.e., the principal component analysis; see, e.g., \cite{FFL08, O12}. The least squares objective function is
\begin{equation}\label{eq:49}
\text{tr}\left[(\bm{Z} - \bm{AF}^\top)(\bm{Z} - \bm{AF}^\top)^\top\right].
\end{equation}
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~\eqref{eq:50}. Substituting in~\eqref{eq:49}, we obtain the objective function
\begin{align*}
&\text{tr}\left[(\bm{Z} - \bm{AA}^\top\bm{Z}/p)(\bm{Z} - \bm{AA}^\top\bm{Z}/p)^\top\right] \\
&=
\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)\\
&= \text{tr}(\bm{ZZ}^\top) - \text{tr}(\bm{A}^\top\bm{ZZ}^\top\bm{A})/p,
\end{align*}
where the last equality uses~\eqref{eq:50} 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 $\frac{1}{p}\bm{ZZ}^\top$ in descending order, where
\begin{equation*}
\frac{1}{p}\bm{Z}\bm{Z}^\top = \frac{1}{p}\sum_{i=1}^n\bm{\eta}_i\bm{\eta}_i^\top = \frac{1}{p}\sum_{i=1}^n(\bm{Y}_i-\bm{\Phi}\bm{c}_i)(\bm{Y}_i-\bm{\Phi}\bm{c}_i)^\top.
\end{equation*}
Therefore, knowing $\bm{c}_i$, we solve for $\widehat{\bm{A}}$ using
\begin{equation}\label{eq:9}
\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},
\end{equation}
where $\bm{V}_{np}$ is an $r\times r$ diagonal matrix containing the $r$ eigenvalues of the matrix in the square brackets in decreasing order. The additional coefficient $1/n$ is used for scaling.

\begin{remark}
The number of factors $r$ is selected based on some criteria regarding the eigenvalues. There have been many studies on this topic. Examples include \cite{BN02}, where two model selection criteria functions were proposed; \cite{O10}, where the number of factors was estimated using differenced adjacent eigenvalues; and \cite{AH13}, where they select the number based on the ratio of two adjacent eigenvalues. A modified ratio criterion of \cite{AH13} is proposed in Section~\ref{se:6.2}. It is worth noting that, simulations and empirical applications in this paper illustrate that our model has good performances under various common selection criteria provided in the literature above.
\end{remark}

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
\begin{align}\label{eq:1}
\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}
\end{align}
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:

\begin{singlespacing}
\small
\begin{algorithm}[H]
\SetAlgoLined
\begin{enumerate}
\item
Denote the initial value as $\widehat{\bm{A}}^{(0)}$. Using~\eqref{eq:1}, 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$.
\item
With $\widehat{\bm{c}}_i^{(t)}$, we substitute into the second equation of~\eqref{eq:1}, 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 $(np)^{-1}\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.
\item
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~\eqref{eq:1}.
\item
We then repeat steps 2 and 3 until $\|\widehat{\bm{c}}_i^{(t+1)} - \widehat{\bm{c}}_i^{(t)}\| < \delta$, where $\delta$ is a prescribed small positive value.
\end{enumerate}
 \caption{Iterations for estimating FASM}\label{ag}
\end{algorithm}

\end{singlespacing}

\normalsize
\begin{remark}
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.
\end{remark}

\begin{remark}\label{re:4}
Common methods for selecting the shrinkage parameter $\alpha$ include the Akaike's Information Criterion \citep[AIC][]{A74}, the Bayesian Information Criterion \citep[BIC][]{S78}, and cross-validation. In this paper, we use the mean generalized cross-validation (mGCV) method \citep{GHW79}. We define, at step $t$,
\begin{equation}\label{eq:53}
\text{mGCV}^{(t)} = \frac{1}{n}\sum_{i=1}^n \frac{\text{pSSE}_i^{(t)}}{[p-df^{(t)}(\alpha)]^2},
\end{equation}
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
\begin{equation}\label{eq:8}
\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].
\end{equation}
At each step of the iteration, the tuning parameter $\alpha$ is chosen by minimizing the $mGCV^{(t)}$.
\end{remark}

\begin{remark}
Algorithm~\ref{ag} is an iteration procedure in which ridge regression and PCA are iterated. The convergence of this iterative algorithm is studied in \cite{JY20}. For instance, Theorem 2 of \cite{JY20} provides some sufficient conditions under which the recursive algorithm converges to the true value or some other values. In particular, this algorithm will converge to the true parameter when the regressors are independent of the common factors, or the factors involved in regressors are weaker than the common ones.
\end{remark}

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
\begin{equation*}
\widehat{\bm{F}}^\top = \widehat{\bm{A}}^\top (\bm{Y}-\bm{\Phi}\widehat{\bm{C}}).
\end{equation*}
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~\eqref{eq:2}.

\begin{remark}
Although we have imposed the constraint in \eqref{eq:50} and the identification condition~\ref{id}, $\bm{A}$ and $\bm{f}_i$ are not uniquely determined, since the model~\eqref{eq:14} 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~\eqref{eq:51} 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.
\end{remark}

\section{Asymptotic Theory}\label{se:4}

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~\ref{se:4.1} and provide the asymptotic results of $\widehat{\bm{c}}_i$ in Section~\ref{se:4.2}.

\subsection{Assumptions}\label{se:4.1}
We use $(\bm{c}_i^0, \bm{A}^0)$ to denote the true parameters. In this paper, we use $L^2$ norm.  The norm of a vector $\bm{x}$ is defined as $\|\bm{x}\| = \sqrt{\bm{x}^\top \bm{x}}$ and the norm of a matrix $\bm{U}$ is defined as $\|\bm{U}\| = \sqrt{\lambda_{\max}(\bm{U}^\top\bm{U})}$, where $\lambda_{\max}(\cdot)$ represents the largest eigenvalue of a matrix. We introduce the matrix
\begin{equation}\label{eq:15}
\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.
\end{equation}
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  \ref{app:a}~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~\ref{as:3}.

First, we state the assumptions.

\begin{assumption}\label{as:1}
\begin{equation*}
\frac{1}{p}\left\|\bm{\Phi}^\top\bm{\Phi}\right\| = O(1), \ \ \left|\left|\textbf{R}\right|\right|=O(1), \ \ K = o(\min(n,p)), \ \ p=o\left(n^2\right).
\end{equation*}
\end{assumption}
In this assumption, the number of basis functions $K$ could be either fixed or tending to infinity.

\begin{assumption}\label{as:2}
\begin{equation*}
\|\bm{c}_i^0\| = O_p(1), \ \text{for all }i=1, 2, \ldots, n.
\end{equation*}
\end{assumption}
This assumption is introduced to ensure the consistency of the estimated coefficients $\widehat{\bm{c}}^0_i$. Notice that the dimension of $\bm{c}_i^{0}$ is $K$, which may go to infinity, but its norm is restricted to be of constant order. Thus this assumption meanwhile controls the flexibility of the model.
\begin{assumption}\label{as: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
\begin{equation*}
\inf_{\bm{A}\in \mathcal{A}} \bm{D}_i(\bm{A}) > 0.
\end{equation*}
\end{assumption}

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~\eqref{eq:15}. The second term on the right-hand side of~\eqref{eq:15} arises because of the unobservable matrices $\bm{F}$ and $\bm{A}$.

\begin{assumption}\label{as:4}
For some constant $M>0$,
\begin{enumerate}
\item
$\mathbb{E}\|\bm{a}^0_j\|^4 \le M$, $j = 1,\dots, p$ and $\frac{1}{p}\bm{A}^\top\bm{A} \overset{p}{\longrightarrow}\bm{\Sigma}_{\bm{a}} >0$ for some $r\times r$ matrix $\bm{\Sigma}_{\bm{a}}$ as $p\rightarrow \infty$.
\item
$\mathbb{E}\|\bm{f}_i\|^4 \le M$ and $\frac{1}{n}\bm{F}\bm{F}^\top \overset{p}{\longrightarrow}\bm{\Sigma}_{\bm{f}} >0$ for some $r\times r$ matrix $\bm{\Sigma}_{\bm{f}}$ as $n\rightarrow \infty$.
\end{enumerate}
\end{assumption}

\begin{assumption}\label{as: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 \le M$.
\end{assumption}

\begin{assumption}\label{as:6}
$\epsilon_{ji}$ is independent of $\phi_{s}$, $\bm{f}_t$, and $\bm{a}^0_s$ for all $j,i,s,t$.
\end{assumption}
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~\ref{as: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.

\begin{assumption}\label{as:7}
The tuning parameter $\alpha$ satisfies $\alpha = o(p)$.
\end{assumption}
This is conventionally assumed in ridge regression \citep[see, e.g.,][]{KF00} 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
\begin{equation}\label{eq:54}
\bm{M}_{\bm{F}} \equiv \bm{I}_n - \bm{F}\left(\bm{F}^\top\bm{F}\right)\bm{F}^\top.
\end{equation}

Then, for any vector  $\bm{b} = (b_1, \dots, b_n)^\top$, we can write
\begin{equation}\label{eq:44}
\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}.
\end{equation}
In~\eqref{eq:44}, 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
\begin{equation}\label{eq:48}
\bm{L}_{np} \equiv \frac{\sigma^2}{np}\sum_i^n\sum_j^p \bm{\omega}_j \bm{\omega}^{\top}_j\left(\sum_k^n \psi_{ik}b_k\right)^2.
\end{equation}
We make the following assumption.

\begin{assumption}\label{as:8}
When $K$ is fixed, we assume there exists a $K\times K$ matrix $\bm{L}$ such that
\begin{equation}\label{eq:45}
\bm{L} \equiv \lim_{n,p\rightarrow \infty}\bm{L}_{np},
\end{equation}
where $\bm{L}_{np}$ is defined in~\eqref{eq:48}. Let $\nu^2$ be the smallest eigenvalue of the matrix $\bm{L}$ defined in~\eqref{eq:45}, then assume that $\nu^2 >0$, and that, for all $\varepsilon >0$,
\begin{equation}\label{eq:56}
\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{I}\left(\left\|\bm{x}_{ij}\right\|^{2} \geq \varepsilon np \nu^{2}\right)\right]=0,
\end{equation}
where $\mathbb{I}(\cdot)$ is the indicator function, defined by $\mathbb{I}(A) = 1$ if $A$ occurs and $\mathbb{I}(A) = 0$ otherwise.

\end{assumption}
This assumption is the multivariate Lindeberg condition, which is needed in constructing the central limit theorem when $K$ is fixed. 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}_jb_i\epsilon_{ji}$. Since we assume $\bm{\phi}_j = O(1)$ in Assumption~\ref{as:1}, the Lindeberg condition in \eqref{eq:56} is met.

Next, we introduce a similar assumption to Assumption~\ref{as:8}, used for the central limit theorem when $K$ goes to infinity. Let $\widetilde{\bm{\omega}}_j, \ j =  1, \dots, p$ denote the $j$th column of the $K\times p$ matrix $\left(\frac{1}{p}\Phi^\top \bm{M}_{\bm{A}^0}\bm{\Phi}\right)^{-1}\bm{\Phi}^\top\bm{M}_{{\bm{A}}^0}$. Denote
\begin{eqnarray}
&\boldsymbol{\gamma}^{\top}\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}}\textbf{b} \nonumber\\
=& \frac{1}{\sqrt{np}}\sum_i^n\sum_j^p\boldsymbol{\gamma}^{\top}\widetilde{\bm{\omega}}_j\epsilon_{ji}\sum_k^n\psi_{ik}b_k
\equiv \frac{1}{\sqrt{np}}\sum_i^n\sum_j^p\widetilde{x}_{ij}.
\end{eqnarray}

\begin{assumption}\label{as:9}
Assume that there exists a constant $\widetilde{L}>0$ such that
\begin{eqnarray}\label{2022(8)}
\widetilde{L}\equiv\lim_{n,p,K\rightarrow\infty}\widetilde{L}_{np},
\end{eqnarray}
where $\widetilde{L}_{np}$ is defined as
\begin{eqnarray}
\widetilde{L}_{np}\equiv \frac{\sigma^2}{np}\sum_i^n\sum_j^p \left(\boldsymbol{\gamma}^{\top}\widetilde{\bm{\omega}}_j\right)^2\left(\sum_k^n \psi_{ik}b_k\right)^2.
\end{eqnarray}
Moreover, for any small value $\varepsilon>0$, we assume
\begin{equation}\label{eq:566}
\lim _{n, p, K\rightarrow \infty} \frac{1}{np \nu^{2}}\sum_{i=1}^{n}\sum_{j=1}^p \mathbb{E}\left[\left|\widetilde{x}_{ij}\right|^{2} \mathbb{I}\left(\left|\widetilde{x}_{ij}\right|^{2} \geq \varepsilon np \nu^{2}\right)\right]=0,
\end{equation}
where $\mathbb{I}(\cdot)$ is the indicator function, defined by $\mathbb{I}(A) = 1$ if $A$ occurs and $\mathbb{I}(A) = 0$ otherwise.

\end{assumption}

\subsection{Asymptotic properties}\label{se:4.2}

As we have mentioned previously, the identification problem of the latent factor implies that we use the estimator $\widehat{\bm{A}}$ to estimate a rotation of $\bm{A}^0$. Based on the objective function~\eqref{eq:22} in Section~\ref{se:3}, we use a center-adjusted objective function, defined as below.
\begin{equation}\label{eq:23}
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,
\end{equation}
where $\bm{M}_{\bm{A}}$ is defined in \eqref{eq:51}, satisfying $\bm{A}^\top \bm{A}/p = \bm{I}_r.$ The second term on the right-hand side of~\eqref{eq:23} 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
\begin{equation}\label{eq:36}
(\widehat{\bm{c}}_i ,\widehat{\bm{A}}) = \argmin_{\bm{c}_i, \bm{A}} S_{np} (\bm{c}_i, \bm{A}).
\end{equation}

In the following, we establish the asymptotic properties for the estimated coefficient matrix $\widehat{\bm{C}}$. In Theorem~\ref{th:1}, the consistency of the matrix $\widehat{\bm{C}}$ is proved. In Theorem~\ref{th:2}, we show the rate of convergence of $\widehat{\bm{C}}$. Theorem~\ref{th: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}$.
\begin{theorem}\label{th:1}
Under Assumptions~\ref{as:1} - \ref{as:6}, as $n, p \rightarrow \infty$, we have the following statements
\begin{enumerate}
\item[(i)]
$\frac{1}{\sqrt{n}}\left\|\bm{C}-\widehat{\bm{C}}\right\|  \overset{p}{\to} 0.$
\item[(ii)]
$\left\|\bm{P}_{\widehat{\bm{A}}}-\bm{P}_{\bm{A}^0}\right\|  \overset{p}{\to} 0$.
\end{enumerate}
\end{theorem}
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 $(\sqrt{n})^{-1}$ 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.
\begin{theorem}\label{th:2}
Under Assumptions~\ref{as:1} - \ref{as:6}, if $p/n\rightarrow \rho >0$ as $n, p\rightarrow\infty$,
\begin{equation*}
\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),
\end{equation*}
where $\bm{M}_{\bm{F}}$ is defined in~\eqref{eq:54}.
\end{theorem}

We study the case when the dimension $p$ and the sample size $n$ are comparable. We achieve rate $\sqrt{p}$ convergence, considering $(\sqrt{n})^{-1}\left\|\bm{C}^0-\widehat{\bm{C}}\right\|$ 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.

\begin{remark}
The asymptotic result in Theorem~\ref{th: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.
\end{remark}

We further begin to establish the limiting distribution. It is shown in \ref{app:a}~Appendix A that
\begin{equation*}
\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).
\end{equation*}
The limiting distribution is constructed based on the first term on the right-hand side.

\begin{theorem}\label{th:3}
In addition to Assumptions~\ref{as:1} -~\ref{as:8}, we assume that $K$ is fixed and, as $n, p\rightarrow\infty$, there exists a positive definite matrix $\widetilde{\textbf{Q}}$ such that
\begin{equation*}
\bm{Q}\left(\bm{A}^0\right) \equiv \frac{1}{p}\bm{\Phi}^\top \bm{M}_{\bm{A}^0} \bm{\Phi}\stackrel{p}{\longrightarrow}\widetilde{\textbf{Q}}.
\end{equation*}
Then we have, for any vector $\bm{b} \in \mathbb{R}^n$
\begin{equation*}
\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}, \widetilde{\textbf{Q}}^{-1}\bm{L}\widetilde{\textbf{Q}}^{-1}\right),
\end{equation*}
where $\bm{M}_{\bm{F}}$ is defined in Theorem~\ref{th:2}, $\bm{L}$ is defined in~\eqref{eq:45}. \\
Moreover if the limit of $\bm{\Phi}\widetilde{\textbf{Q}}^{-1}\bm{L}\widetilde{\textbf{Q}}^{-1}\bm{\Phi}^\top$ exists, we have
\begin{equation*}
\sqrt{p}\bm{\Phi}\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}, \lim_{n,p\rightarrow\infty}\bm{\Phi}\widetilde{\textbf{Q}}^{-1}\bm{L}\widetilde{\textbf{Q}}^{-1}\bm{\Phi}^\top\right),
\end{equation*}
\begin{equation*}
\sqrt{\frac{p}{n}}\left(\bm{X}-\widehat{\bm{X}}\right)\bm{M}_{\bm{F}}\bm{b} \overset{d}{\to} \mathcal{N}\left(\bm{0}, \lim_{n,p\rightarrow\infty}\bm{\Phi}\widetilde{\textbf{Q}}^{-1}\bm{L}\widetilde{\textbf{Q}}^{-1}\bm{\Phi}^\top\right).
\end{equation*}
\end{theorem}
The vector $\bm{b}$ comes from the same vector in Lemma~\ref{le:4}. The asymptotic bias is zero since we assume no serial or cross-sectional correlation in the error terms. This simplified setting 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.

\begin{remark}
Theorem~\ref{th: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}}}$,
\begin{align*}
\widehat{\bm{Q}} &= \frac{1}{p}\bm{\Phi}^\top \bm{M}_{\widehat{\bm{A}}} \bm{\Phi}\\
\widehat{\bm{L}} &= \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,
\end{align*}
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}}}$.
\end{remark}

When $K$ goes to infinity, the random vector $\left(\frac{\bm{C}^0-\widehat{\bm{C}}}{\sqrt{n}}\right)\bm{M}_{\bm{F}}\bm{b}$ considered in Theorem \ref{th:3} is a high-dimensional random vector. Then we establish the asymptotic distribution for some linear combination of this random vector in the following theorem.
\begin{theorem}\label{th:4}
Suppose Assumptions ~\ref{as:1} -~\ref{as:7} and Assumption~\ref{as:9} are satisfied. As $n, p, K$ tend to infinity, we have for any vectors $\boldsymbol{\gamma}\in \mathbb{R}^{K}$ and $\bm{b} \in \mathbb{R}^n$,
\begin{equation}
\sqrt{p}\boldsymbol{\gamma}^{\top}\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}, \widetilde{L}\right),
\end{equation}
where $\widetilde{L}$ is defined in (\ref{2022(8)}).
\end{theorem}

\section{Simulation Studies}\label{se:6}

In this section, we use simulated data to illustrate the superiority of the proposed model. We compare FASM with three other modes. The first one is the basis smoothing model with a penalty (Bsmooth). We use the same basis functions as FASM and the L2 penalty. So the Bsmooth model is only different from FASM in that it does not consider the augmented factor component. The second one is the local linear smoothing model (Localin). Cross-validation is used to choose the bandwidth parameter. The third one is the principal component analysis with conditional expectation (PACE). This method is proposed in \cite{YMW05} and is mainly used for sparse data, where very few data points are observed on each functional curve.

\subsection{Data generating process}\label{se:6.1}

\subsubsection*{Setting 1}

We generate simulated data $Y_{ij}$, where $i = 1,\dots, n$ and $j = 1,\dots, p$ from the following model:
\begin{align*}
Y_{ij} &= \mathcal{X}_i(u_j) + \eta_{ij} + \epsilon_{ji}\\
&= \sum_{k=1}^{13} c_{ik} \phi_k(u_j) + \sum_{k=1}^4 f_{ik}a_{kj} + \epsilon_{ji},
\end{align*}
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 loadings $a_{kj}$ follow $\mathcal{N}(0, 0.6^2)$ and the factors $(f_{i1}, f_{i2}, f_{i3}, f_{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)$.

\subsubsection*{Setting 2}

We generate simulated data $Y_{ij}$, where $i = 1,\dots, n$ and $j = 1,\dots, p$ from the following model:
\begin{align*}
Y_{ij} &= \mathcal{X}_i(u_j) + \eta_{ij} + \epsilon_{ji}\\
&= \sum_{k=1}^{5} c_{ik} \phi_k(u_j) + \sum_{k=1}^4 f_{ik}a_{kj} + \epsilon_{ji},
\end{align*}
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, \gamma_k^2)$, where $\gamma_1 = 3,\gamma_2 = 2.5, \gamma_3= 2, \gamma_4= 1.5, \gamma_5= 1$. The factors loadings $a_{kj}$ follow $\mathcal{N}(0, 0.8^2)$ and the factors $(f_{i1}, f_{i2}, f_{i3}, f_{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)$.

\subsubsection*{Setting 3}

We generate simulated data as in Setting 1, except we use $K = 21$ B-spline basis functions. We use $n = 40$ and $p  = 50$.  We show how the model works when the number of basis functions is large compared to $n$ and $p$.

\subsection{Results}\label{se:6.3}

We repeat the simulation setup 500 times and obtain the estimated smooth function $\widehat{\mathcal{X}}_i(u) = \widehat{\bm{c}}_i^\top \bm{\Phi}(u)$.  The mean integrated squared error (MISE) for function estimation is calculated as
\begin{equation}\label{eq:100}
\text{MISE} = \frac{1}{n}\sum_{i=1}^n \int\left[ \mathcal{X}_i(u)-\widehat{\mathcal{X}}_i(u)\right]^2du.
\end{equation}

The results for Setting 1 and 2 are reported in Table~\ref{tab:1} and~\ref{tab:2}. 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~\ref{se:4}, where the estimator converges with the rate related to $p$.

When $\sigma$ is large, i.e., the signal-to-noise ratio is high, the FASM performs better than the smoothing model. When $\sigma = 0$, the actual data are generated without the factor model component. The performance of the two methods is identical. This implies that the proposed model is robust in that it improves the estimation when there is measurement error, and it is no worse than the ordinary smoothing model when there is no measurement error. This should be credited to the selection criterion used in \eqref{eq:58}. It always selects zero factors when there are no common factors.

In Table~\ref{tab:3}, we show the MISE result from Setting 3. Our proposed model also shows its supiority when $K$ is large compared to $n$ and $p$. This setting corresponds to allowing the number of basis functions $K$ to go to infinity in Section~\ref{se:4}.

\begin{singlespace}

\begin{table}[!htb]
\caption{{\bf Setting 1}: The MISE 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.}
\label{tab:1}\par
\centering
{\begin{tabular}{|cclcccc|} \hline
&&&
    			\multicolumn{4}{c|}{MISE} \\
  Sample Size & Dimension & Size of $\eta$ & FASM & Bsmooth &  Localin & PACE\\\hline
  $n = 20$ & $ p = 51$&   $\sigma = 0$& $\bm{0.058}$ & $\bm{0.058}$ & $0.075$ & $0.238$ \\
 & & $\sigma = 0.5$& $\bm{0.131}$ & $\bm{0.131}$ & $0.149$ & $0.587$ \\
&&$\sigma = 0.75$ & $\bm{0.204}$ & $0.208$ & $0.226$ & $1.026$ \\
&& $\sigma = 1$& $\bm{0.297}$ & $0.318$ & $0.336$ & $1.676$ \\ \hline
  $n = 20$ & $p = 101$ &   $\sigma = 0$&$\bm{0.031}$ & $\bm{0.031}$ & $0.043$ & $0.238$ \\
&&$\sigma = 0.5$ & $\bm{0.076}$ & $\bm{0.076}$ & $0.089$ & $0.595$ \\
&&$\sigma = 0.75$  & $\bm{0.115}$ & $0.123$ & $0.136$ & $1.020$ \\
&& $\sigma = 1$& $\bm{0.159}$ & $0.187$ & $0.201$ & $1.686$ \\ \hline
  $n = 50$ & $p = 51$&   $\sigma = 0$ & $\bm{0.058}$ & $\bm{0.058}$ & $0.076$ & $0.195$ \\
&& $\sigma = 0.5$ & $\bm{0.131}$ & $0.135$ & $0.153$ & $0.543$ \\
&&$\sigma = 0.75$ & $\bm{0.186}$ & $0.214$ & $0.233$ & $1.013$ \\
&& $\sigma = 1$& $\bm{0.248}$ & $0.314$ & $0.331$ & $1.668$ \\\hline
  $n = 100$ & $ p = 101$ &   $\sigma = 0$& $\bm{0.031}$ & $\bm{0.031}$ & $0.044$ & $0.170$ \\
 && $\sigma = 0.5$ & $\bm{0.062}$ & $0.074$ & $0.089$ & $0.463$ \\
&&$\sigma = 0.75$  & $\bm{0.090}$ & $0.118$ & $0.132$ & $0.908$ \\
&& $\sigma = 1$ & $\bm{0.133}$ & $0.181$ & $0.194$ & $1.570$ \\ \hline
\end{tabular}}
\end{table}

\end{singlespace}

\begin{singlespace}

\begin{table}[!htb]
\caption{{\bf Setting 2}: The MISE 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.}
\label{tab:2}\par
\centering
{\begin{tabular}{|cclcccc|} \hline
&&&
    			\multicolumn{4}{c|}{MISE} \\
  Sample Size & Dimension & Size of $\eta$ & FASM & Bsmooth &  Localin & PACE\\\hline
  $n = 20$ & $ p = 51$&   $\sigma = 0$& $\bm{0.024}$ & $\bm{0.024}$ & $0.070$ & $0.116$ \\
 & & $\sigma = 0.5$ & $0.086$ & $\bm{0.085}$ & $0.195$ & $0.621$ \\
&&$\sigma = 0.75$ & $\bm{0.160}$ & $\bm{0.160}$ & $0.340$ & $1.431$ \\
&& $\sigma = 1$ & $\bm{0.272}$ & $0.280$ & $0.538$ & $2.685$ \\  \hline
  $n = 20$ & $p = 101$ &   $\sigma = 0$  & $\bm{0.012}$ & $\bm{0.012}$ & $0.040$ & $0.089$ \\
&&$\sigma = 0.5$ & $0.047$ & $\bm{0.046}$ & $0.117$ & $0.660$ \\
&&$\sigma = 0.75$  & $\bm{0.079}$ & $\bm{0.079}$ & $0.192$ & $1.434$ \\
&& $\sigma = 1$& $\bm{0.135}$ & $0.138$ & $0.301$ & $2.549$ \\  \hline
  $n = 50$ & $p = 51$&   $\sigma = 0$ & $\bm{0.025}$ & $\bm{0.025}$ & $0.070$ & $0.076$ \\
&& $\sigma = 0.5$  & $0.087$ & $\bm{0.086}$ & $0.202$ & $0.584$ \\
&&$\sigma = 0.75$ & $\bm{0.162}$ & $0.169$ & $0.348$ & $1.435$ \\
&& $\sigma = 1$ & $\bm{0.250}$ & $0.269$ & $0.522$ & $2.557$ \\ \hline
  $n = 100$ & $ p = 101$ &   $\sigma = 0$ & $\bm{0.012}$ & $\bm{0.012}$ & $0.041$ & $0.033$ \\
 && $\sigma = 0.5$ & $\bm{0.045}$ & $\bm{0.045}$ & $0.114$ & $0.518$ \\
&&$\sigma = 0.75$ & $\bm{0.085}$ & $\bm{0.085}$ & $0.195$ & $1.337$ \\
&& $\sigma = 1$ & $\bm{0.138}$ & $0.139$ & $0.299$ & $2.618$ \\ \hline
\end{tabular}}
\end{table}

\end{singlespace}

\begin{singlespace}

\begin{table}[!htb]
\caption{{\bf Setting 3}: The MISE of the function estimates when $K$ is large compared to $n$ and $p$. The adjustment of $\sigma^2$ value is used to control the signal-to-noise ratio.}
\label{tab:3}\par
\centering
{\begin{tabular}{|cclcccc|} \hline
&&&
    			\multicolumn{4}{c|}{MISE} \\
  Sample Size & Dimension & Size of $\eta$ & FASM & Bsmooth &  Localin & PACE\\\hline
  $n = 50$ & $ p = 40$ &   $\sigma = 0$ & $\bm{0.097}$ & $\bm{0.097}$ & $0.121$ & $0.225$ \\
 && $\sigma = 0.5$ & $\bm{0.206}$ & $0.207$ & $0.235$ & $0.590$ \\
&&$\sigma = 0.75$ & $\bm{0.291}$ & $0.309$ & $0.335$ & $1.044$ \\
&& $\sigma = 1$& $\bm{0.390}$ & $0.457$ & $0.466$ & $1.747$ \\  \hline
\end{tabular}}
\end{table}

\end{singlespace}

\section{Application to Climatology}\label{se:7}

We apply the aforementioned FASM, Bsmooth, Localin and PACE methods to some real data sets. This section looks at the Friday temperature data at Adelaide airport. We also provide the analysis of Canadian yearly temperature and precipitation data in the Supplement.

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 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~\ref{fig:7a}. 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.  We show the residuals from the four models in Figure~\ref{fig:8}.  Apart from FASM, all the other models produce large residuals with some extreme values. To further look at the residual struture from the Bsmooth model, we conduct a principal component analysis on the residuals. In Figure~\ref{fig:7c}, the screeplot of the eigenvalues are presented. The residuals from Bsmooth model exhibit a ``spike" structure, where the first few eigenvalues are significantly larger than the rest. This means the residuals contain information that can be captured by just a few factors, which calls for a further dimension reduction model on the residuals.

\begin{figure}[!htb]
\centering
\subcaptionbox{Raw temperature data\label{fig:7a}}
{\includegraphics[width = 8.41cm]{Friday.pdf}}
\subcaptionbox{Eigenvalues of the residuals in from Bsmooth \label{fig:7c}}
{\includegraphics[width = 8.41cm]{friday_spike.pdf}}
\caption{Half-hourly temperature data on Friday at Adelaide airport.}\label{fig:7}
\end{figure}

\begin{figure}[!htb]
\centering
\subcaptionbox{Residuals of FASM\label{fig:8a}}
{\includegraphics[width = 8.4cm]{friday_fasm_resid.pdf}}
\subcaptionbox{Residuals of the Bsmooth\label{fig:8b}}
{\includegraphics[width = 8.4cm]{friday_basis_resid.pdf}}
\subcaptionbox{Residuals of Localin\label{fig:8c}}
{\includegraphics[width = 8.4cm]{friday_locpol_resid.pdf}}
\subcaptionbox{Residuals of PACE\label{fig:8d}}
{\includegraphics[width = 8.4cm]{friday_pace_resid.pdf}}
\caption{Half-hourly temperature data on Friday at Adelaide airport.}\label{fig:8}
\end{figure}

\section{Conclusion and Future Work}\label{se:8}

In this paper, we propose a factor-augmented smoothing model for functional data. We study raw functional data, 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 obtain estimates in the smoothing and factor models simultaneously. 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 constructing 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 misspecification 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 Australian temperature data.

The proposed model is a good starting point for modeling complex data structures. We deal with 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. The change point is a popular problem in many statistics and econometric topics and has been extensively studied in the multivariate setting. Previous literature on change points in functional data include \cite{BG09, HK10} and \cite{HK15}. It is shown in simulation examples that our model can be used for modeling functional data with change points in the cross-sectional direction. The model can be modified to account for change points also in the sample direction. Further research can be conducted along this line.

\section*{Supplementary Materials}

The Supplement includes our extension of the model to non-parametric settings, introducing a covariance estimator and some extensive simulated and real data analysis. It also contains the proofs for the theorems in the main article. Due to the page limit, many important results are placed in the Supplement of this paper.

\bibhang=1.7pc
\bibsep=2pt
\fontsize{12}{14pt plus.8pt minus .6pt}\selectfont
\renewcommand\bibname{\large \bf References}
\expandafter\ifx\csname
natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL}\fi

\newpage
\bibliographystyle{chicago}

redpajama_set_name: RedPajamaArXiv