# Asymptotic equivalence of Principal Components and Quasi Maximum Likelihood
estimators in Large Approximate Factor Models

###### Abstract

This paper investigates the properties of Quasi Maximum Likelihood estimation
of an approximate factor model for an $n$-dimensional vector of stationary
time series. We prove that the factor loadings estimated by Quasi Maximum
Likelihood are asymptotically equivalent, as $n\to\infty$, to those estimated
via Principal Components. Both estimators are, in turn, also asymptotically
equivalent, as $n\to\infty$, to the unfeasible Ordinary Least Squares
estimator we would have if the factors were observed. We also show that the
usual sandwich form of the asymptotic covariance matrix of the Quasi Maximum
Likelihood estimator is asymptotically equivalent to the simpler asymptotic
covariance matrix of the unfeasible Ordinary Least Squares. All these results
hold in the general case in which the idiosyncratic components are cross-
sectionally heteroskedastic, as well as serially and cross-sectionally weakly
correlated. The intuition behind these results is that as $n\to\infty$ the
factors can be considered as observed, thus showing that factor models enjoy a
blessing of dimensionality.

Keywords: Approximate Dynamic Factor Model; Principal Component Analysis;
Quasi Maximum Likelihood.

Matteo Barigozzi∗

This version: August 28, 2024

$\ast$$\ast$footnotetext:  Università di Bologna<EMAIL_ADDRESS>
##  1 Introduction

Factor models are one of the major dimension reduction techniques used to
analyze large panels of time series. Some of their most successful
applications are, among many others, in finance (Chamberlain and Rothschild,
1983, Connor et al., 2006, Aït-Sahalia and Xiu, 2017, Kim and Fan, 2019), and
macroeconomics (Stock and Watson, 2002b, Bernanke et al., 2005, Forni et al.,
2005, De Mol et al., 2008, Giannone et al., 2008).

Let us assume to observe an $n$-dimensional zero-mean stochastic process over
$T$ periods: $\\{x_{it},\,i=1,\ldots,n,\,t=1,\ldots,T\\}$, such that

$\displaystyle x_{it}$
$\displaystyle=\bm{\lambda}_{i}^{\prime}\mathbf{F}_{t}+\xi_{it},\quad
i=1,\ldots,n,\quad t=1,\ldots,T,$ (1)

where $\bm{\lambda}_{i}=(\lambda_{i1}\cdots\lambda_{ir})^{\prime}$ and
$\mathbf{F}_{t}=(F_{1t}\cdots F_{rt})^{\prime}$ are $r$-dimensional unobserved
vectors, called loadings and factors, respectively, and with $r<\min(n,T)$. We
also call $\xi_{it}$ the idiosyncratic component and
$\chi_{it}=\bm{\lambda}_{i}^{\prime}\mathbf{F}_{t}$ the common component of
the $i$th observed variable at time $t$. Both the factors and the
idiosyncratic components are allowed to be serially correlated. Furthermore,
the idiosyncratic components are also allowed to be cross-sectionally
correlated. We call such model an approximate factor model. This is the class
of factor models we consider in this paper and studied, e.g., by Bai (2003).
It is a restricted version of the generalized dynamic factor model originally
proposed by Forni et al. (2000), where the factors are loaded with lags (the
loadings are linear filters) and not just contemporaneously. Another popular
alternative, not considered in this paper, is the factor model studied, e.g.,
by Lam and Yao (2012), where the idiosyncratic components are instead assumed
to be serially uncorrelated.

There are two main ways to estimate the factor loadings in (1). First, by
Principal Component (PC) analysis (Stock and Watson, 2002a; Bai, 2003; Fan et
al., 2013), and, second, by Quasi Maximum Likelihood (QML) estimation (Bai and
Li, 2012, 2016). In both cases, we can first estimate the loadings and then
estimate the factors by linear projection, possibly weighted by the
idiosyncratic variances, of the observables onto the estimated loadings. The
PC estimator of the loadings is easily implementable since it just requires to
compute the $r$ leading eigenvectors and eigenvalues of the sample covariance
matrix of the observables. The QML estimator of the loadings does not have a
closed form solution, so numerical maximization is required. It is typically
defined as the maximizer of a mis-specified log-likelihood where the
idiosyncratic components are treated as if they were uncorrelated even if the
true ones are correlated.

Although QML is the classical way to estimate a factor model, dating back more
than fifty years ago (Lawley and Maxwell, 1971; Barigozzi, 2024, and
references therein), in the recent years, PC analysis has gained popularity,
given its non-parametric nature and ease of implementation. Nevertheless, QML
estimation still retains an important role in factor analysis since for many
reasons. For example, it allows to easily impose constraints on the loadings
(see, e.g.,the applications in Coroneo et al., 2016; Delle Chiaie et al.,
2021), and it fully addresses idiosyncratic cross-sectional
heteroskedasticity, while the PC method does not (see the comments in Bai and
Li, 2012, 2016).

In this paper, we compare the PC and QML estimators of the loadings and we
show that, under a minimal unified set of standard assumptions and identifying
constraints, the two estimators are asymptotically equivalent as $n\to\infty$
(see Theorem 3). It is well known that both estimators of the loadings are
$\min(n,\sqrt{T})$-consistent as $n,T\to\infty$ and are also asymptotically
normal if we assume $n^{-1}\sqrt{T}\to 0$. Standard references for these
results are, among others, Bai (2003, Theorem 2) for PC, and Bai and Li (2016,
Theorem 1) for QML. Proving the asymptotic properties of the PC estimator is
long but straightforward. However, things are much more complicated for the
QML estimator, essentially because no closed form expression exists for this
estimator. Moreover, the existing proofs for the two estimators do not make
use of the same assumptions nor of the same identification constraints needed
to uniquely identify the loadings.

Our main result has some important implications. First, if $n^{-1}\sqrt{T}\to
0$, as $n,T\to\infty$, then the QML estimator has the same asymptotic
distribution as the unfeasible OLS estimator we would obtain if the factors
were observed. Therefore, we are able to prove asymptotic normality of the QML
estimator in an indirect and easy way (see Theorem 4). Our approach provides a
much simpler alternative to the approach used by Bai and Li (2016, Theorem 1),
and it does not require to study also the properties of the QML estimators of
the idiosyncratic variances. Second, it nests the case of spherical
idiosyncratic components (Tipping and Bishop, 1999), in which case the
equivalence of PC and QML estimator is often quoted, but, to the best of our
knowledge, it has never been formally proved, at least under the standard
minimal set of assumptions used in this paper. Related to these two points, we
stress that our result is, in fact, more general since it holds without the
need of considering QML estimation based on a mis-specified log-likelihood
with a diagonal idiosyncratic covariance. Third, and last, our result paves
the way towards studying the asymptotic properties of the QML estimator of a
factor model where we also explicitly model the dynamics of the factors, and
for which only partial results exist (Doz et al., 2012, Bai and Li, 2016).

The theoretical analysis of approximate factor models requires a double
asymptotic framework. Indeed, only if $n,T\to\infty$ we can consistently
estimate the eigenvectors of the covariance matrix of the data which are
needed for PC estimation, and we can control for the mis-specifications
introduced in the log-likelihood when considering QML estimation. This means
that QML estimation of approximate factor models does not fall into the
framework of classical QML estimation, where $n$ is fixed.

In a second contribution, we reconcile the asymptotic distribution of the QML
estimator of the loadings, derived in this paper, with the results of
classical likelihood-based inference. In particular, we show that, when
evaluated in the true value of the parameters or in the value of the QML
estimator, the first and second derivatives of the factor model log-likelihood
we maximize are asymptotically, as $n,T\to\infty$, equivalent to the first and
second derivatives of the log-likelihood we would have if the factors were
observed (see Theorem 5). This result also allows us to compute simple
estimators of the asymptotic covariance matrix of the estimated loadings
without the need of computing the derivatives of the log-likelihood which have
rather long and complex expressions.

The paper is organized as follows. In Section 2 we present the model and all
assumptions. In Section 3 we reconsider the PC estimator of the loadings and
in Theorem 1 we derive its asymptotic properties using a new approach
equivalent to the one proposed by Bai (2003) but which is more convenient for
the present work. In Section 4 we review the existing results on the QML
estimator of the loadings by Bai and Li (2016). In Section 5 we present our
first contribution in Theorem 3 where we prove that PC and QML estimators of
the loadings are asymptotically equivalent. In Section 6 we present our second
contribution in Theorem 5, where we prove the asymptotically equivalence of
the first and second derivatives of the factor model log-likelihood and the
log-likelihood of a model with observed factors. In Section 7 we briefly
discuss estimation of factors. In Section 8 we provide simulation results
confirming our theoretical results. In Section 9 we conclude. In the
Supplementary Material we prove all other main and auxiliary theoretical
results.

##  2 Model and Assumptions

Given the $n$-dimensional vector $\mathbf{x}_{t}=(x_{1t}\cdots
x_{nt})^{\prime}$, we can also write model (1) in vector notation:

$\mathbf{x}_{t}=\bm{\Lambda}\mathbf{F}_{t}+\bm{\xi}_{t},\quad t=1,\ldots,T,$

where $\bm{\xi}_{t}=(\xi_{1t}\cdots\xi_{nt})^{\prime}$ is the $n$-dimensional
vector of idiosyncratic components and
$\bm{\Lambda}=(\bm{\lambda}_{1}\cdots\bm{\lambda}_{n})^{\prime}$ is the
$n\times r$ matrix of factor loadings. We call
$\bm{\chi}_{t}=\bm{\Lambda}\mathbf{F}_{t}$ the vector of common components.

Moreover, we can collected all observations into the $T\times n$ matrix
$\bm{X}=(\mathbf{x}_{1}\cdots\mathbf{x}_{T})^{\prime}$, and we can write model
(1) also in matrix notation:

$\bm{X}=\bm{F}\bm{\Lambda}^{\prime}+\bm{\Xi},$ (2)

where $\bm{\Xi}=(\bm{\xi}_{1}\cdots\bm{\xi}_{T})^{\prime}$ is a $T\times n$
matrix of idiosyncratic components, and
$\bm{F}=(\mathbf{F}_{1}\cdots\mathbf{F}_{T})^{\prime}$ is the $T\times r$
matrix of factors.

The following assumptions are similar to those made by Bai (2003) for PC
estimation while slightly differ from those in Bai and Li (2016) for QML
estimation. Throughout, we highlight the main differences or similarities.

We start by characterizing the common component by means of the following
assumption.

######  Assumption 1 (common component).

  1. (a)

$\lim_{n\to\infty}\|n^{-1}\bm{\Lambda}^{\prime}\bm{\Lambda}-\bm{\Sigma}_{\Lambda}\|=0$,
where $\bm{\Sigma}_{\Lambda}$ is $r\times r$ positive definite, and, for all
$i\in\mathbb{N}$, $\|\bm{\lambda}_{i}\|\leq M_{\Lambda}$ for some finite
positive real $M_{\Lambda}$ independent of $i$.

  2. (b)

For all $t\in\mathbb{Z}$, $\mathbb{E}[\mathbf{F}_{t}]=\mathbf{0}_{r}$ and
$\bm{\Gamma}^{F}=\mathbb{E}[\mathbf{F}_{t}\mathbf{F}_{t}^{\prime}]$ is
$r\times r$ positive definite and $\|\bm{\Gamma}^{F}\|\leq M_{F}$ for some
finite positive real $M_{F}$ independent of $t$.

  3. (c)

(i) For all $t\in\mathbb{Z}$, $\mathbb{E}[\|\mathbf{F}_{t}\|^{4}]\leq K_{F}$
for some finite positive real $K_{F}$ independent of $t$;
(ii) for all $i,j=1,\ldots,r$, all $s=1,\ldots,T$, and all $T\in\mathbb{N}$,
$\mathbb{E}\left[\left|\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\left\\{F_{is}F_{jt}-\mathbb{E}[F_{is}F_{jt}]\right\\}\right|^{2}\right]\leq
C_{F}$ for some finite positive real $C_{F}$ independent of $i$, $j$, $s$, and
$T$.

  4. (d)

There exists an integer $N$ such that for all $n>N$, $r$ is a finite positive
integer, independent of $n$.

Part (a) is standard (Bai, 2003, Assumption B). It implies that,
asymptotically, as $n\to\infty$, the loadings matrix has asymptotically
maximum column rank $r$, and that for any given $n\in\mathbb{N}$, each factor
has a finite contribution to each of the $n$ observed series (upper bound on
$\|\bm{\lambda}_{i}\|$). A similar requirement is in Bai and Li (2016,
Assumption B). We only consider non-random factor loadings for simplicity and
in agreement with classical factor analysis where the loadings are the
parameters of the model (Lawley and Maxwell, 1971, see, e.g.,).

Part (b) assumes that the factors have zero mean and have a finite full-rank
covariance matrix $\bm{\Gamma}^{F}$, so they are non-degenerate. In part (c-i)
we assume finite 4th order moments of the factors. These are standard
requirements (Bai, 2003, Assumption A). Part (c-ii) is very general, it
implies that the sample covariance matrix of the factors is a
$\sqrt{T}$-consistent estimator of its population counterpart
$\bm{\Gamma}^{F}$, which has full rank because of part (b). It is immediate to
see that it is equivalent to asking for 4th order summable cross-cumulants
which a necessary and sufficient condition for consistent estimation (Hannan,
1970, pp. 209-211). This approach is high-level in that it does not make any
specific assumption on the dynamics of $\\{\mathbf{F}_{t}\\}$. Obviously this
implies the usual assumption of convergence in probability made in this
literature:
$\mathrm{P}\text{-}\lim_{T\to\infty}\|T^{-1}{\bm{F}^{\prime}\bm{F}}-\bm{\Gamma}^{F}\|=0$
(Bai, 2003, Assumption A). Nothing would change in our proofs if we directly
assumed this latter condition instead of part (c-ii). Notice that Bai and Li
(2016, Assumption A) treat the factors as being deterministic, but essentially
make the same assumption as our parts (b) and (c).

Part (d) implies the existence of a finite number of factors. In particular,
the number of common factors, $r$, is identified only as $n\to\infty$. Here
$N$ is the minimum number of series we need to be able to identify $r$ so that
$r\leq N$. Hereafter, when we say “for all $n\in\mathbb{N}$” we always mean
that $n>N$ so that $r$ can be identified. In practice, we must always work
with $n$ such that $r<n$. Moreover, because PC estimation is based on
eigenvalues of an $n\times n$ matrix estimated using $T$ observations, then we
must also have samples of size $T$ such that $r<T$. Therefore, sometimes it is
directly assumed that $r<\min(n,T)$.

To characterize the idiosyncratic component, we make the following
assumptions.

######  Assumption 2 (idiosyncratic component).

  1. (a)

For all $i\in\mathbb{N}$ and all $t\in\mathbb{Z}$, $\mathbb{E}[\xi_{it}]=0$
and $\sigma_{i}^{2}=\mathbb{E}[\xi_{it}^{2}]$ is such that
$C_{\xi}\leq\sigma_{i}^{2}\leq C_{\xi}^{\prime}$ for some finite positive
reals $C_{\xi}$ and $C_{\xi}^{\prime}$ independent of $i$ and $t$.

  2. (b)

For all $i,j\in\mathbb{N}$, all $t\in\mathbb{Z}$, and all $k\in\mathbb{Z}$,
$|\mathbb{E}[\xi_{it}\xi_{j,t-k}]|\leq\rho^{|k|}M_{ij}$, where $\rho$ and
$M_{ij}$ are finite positive reals independent of $t$ such that $0\leq\rho<1$,
$\sum_{j=1,j\neq i}^{n}M_{ij}\leq M_{\xi}$, and $\sum_{i=1,i\neq
j}^{n}M_{ij}\leq M_{\xi}$ for some finite positive real $M_{\xi}$ independent
of $i$, $j$, and $n$.

  3. (c)

(i) For all $i=1,\ldots,n$, all $t=1,\ldots,T$, and all $n,T\in\mathbb{N}$,
$\mathbb{E}[\xi_{it}^{4}]\leq Q_{\xi}$ for some finite positive real $Q_{\xi}$
independent of $i$ and $t$;
(ii) for all $j=1,\ldots,n$, all $s=1,\ldots,T$, and all $n,T\in\mathbb{N}$,
$\mathbb{E}\left[\left|\frac{1}{\sqrt{nT}}\sum_{i=1}^{n}\sum_{t=1}^{T}\left\\{\xi_{is}\xi_{jt}-\mathbb{E}[\xi_{is}\xi_{jt}]\right\\}\right|^{2}\right]\leq
K_{\xi}$ for some finite positive real $K_{\xi}$ independent of $j$, $s$, $n$,
and $T$.

By part (a), we have that the idiosyncratic components have zero mean. This,
jointly with Assumption 1(b) by which
$\mathbb{E}[\mathbf{F}_{t}]=\mathbf{0}_{r}$, implies that we are implicitly
assuming that each observed series has zero mean, i.e., $\mathbb{E}[x_{it}]=0$
(this is without loss of generality), and strictly positive variance, i.e.,
$\mathbb{V}\text{ar}(x_{it})>0$, for all $i\in\mathbb{N}$ (Bai and Li, 2016,
Assumption C.2). For the case of non-zero mean see Remark 2.

Part (b) has a twofold purpose. First, it limits the degree of serial
correlation of the idiosyncratic components by imposing geometric decay of
their autocovariances. Second, it also limits the degree of cross-sectional
correlation between idiosyncratic components, which is usually assumed in
approximate factor models. In particular, by setting $k=0$, it follows also
that for all $i\in\mathbb{N}$, $\sigma_{i}^{2}\leq M_{\xi}$. Thus, all
idiosyncratic components have finite variance. This jointly with Assumptions
1(a) and 1(b) implies that each observed time series has finite variance,
i.e., $\mathbb{V}\text{ar}(x_{it})<\infty$, for all $i\in\mathbb{N}$. In Lemma
1 we show that part (b) implies all usual conditions on second order moments
typically found in the literature (Bai, 2003, Assumptions C.2, C.3, C.4, and
E, and Bai and Li, 2016, Assumptions C.3, C.4, and E.1).

Part (c-i) assumes finite 4th order moments of the idiosyncratic components.
Part (c-ii) gives summability conditions across the cross-section and time
dimensions for the 4th order cumulants of $\\{\xi_{it}\\}$. Jointly with
Assumptions 1(a) and 2(a) it nests the requirement in Bai and Li (2016,
Assumption E.2). It implies that the sample (auto)covariances between
$\\{\xi_{it}\\}$ and $\\{\xi_{jt}\\}$ are $\sqrt{T}$-consistent estimators of
their population counterparts. In particular, by choosing $s=t$ we see that we
can consistently estimate the $(i,j)$th entry of the idiosyncratic covariance
matrix $\bm{\Gamma}^{\xi}$. Notice that, contrary to the existing literature
(Bai, 2003, Assumptions C.1 and C.5, and Bai and Li, 2016, Assumptions C.1 and
C.5), there is no need to ask for finite 8th order moments and cumulants.

We then make a series of identifying assumptions.

######  Assumption 3 (independence).

The processes $\\{\xi_{it},\,i\in\mathbb{N},\,t\in\mathbb{Z}\\}$ and
$\\{F_{jt},\,j=1,\ldots,r,\,t\in\mathbb{Z}\\}$ are mutually independent.

This assumption obviously implies that the factors and, therefore, the common
components are independent of the idiosyncratic components at all leads and
lags and across all units. This is compatible for example with a macroeconomic
interpretation of factor models, according to which the factors driving the
common component are independent of the idiosyncratic components representing
measurement errors or local dynamics. This assumption is made for simplicity
and could be easily relaxed (Bai, 2003, Assumption D).

Let
$\bm{\Gamma}^{\chi}=\mathbb{E}[\bm{\chi}_{t}\bm{\chi}_{t}^{\prime}]=\bm{\Lambda}\bm{\Gamma}^{F}\bm{\Lambda}^{\prime}$
with $r$ largest eigenvalues $\mu_{j}^{\chi}$, $j=1,\ldots,r$, sorted in
decreasing order. In Lemma 1(iv) we prove that, for all $j=1,\ldots,r$,

$\underline{C}_{j}\\!\leq\lim\inf_{n\to\infty}\frac{\mu_{j}^{\chi}}{n}\leq\lim\sup_{n\to\infty}\frac{\mu_{j}^{\chi}}{n}\leq\\!\overline{C}_{j},$
(3)

where $\underline{C}_{j}$ and $\overline{C}_{j}$ are finite positive reals.
This means we consider only strong factors, i.e., fully pervasive. Furthermore
in Lemma 1(v) we prove that the largest eigenvalue of $\bm{\Gamma}^{\xi}$ is
such that

$\sup_{n\in\mathbb{N}}\mu_{1}^{\xi}\leq C_{\xi}^{\prime}+M_{\xi},$ (4)

where $C_{\xi}^{\prime}$ and $M_{\xi}$ are defined in Assumptions 2(a) and
2(b), respectively. Because of Weyl’s inequality (see, e.g., Merikoski and
Kumar, 2004, Theorem 1), conditions (3)-(4) and Assumption 3 imply the
eigengap in the eigenvalues $\mu_{j}^{x}$, $j=1,\ldots,n$, of
$\bm{\Gamma}^{x}=\mathbb{E}[\mathbf{x}_{t}\mathbf{x}_{t}^{\prime}]=\bm{\Gamma}^{\chi}+\bm{\Gamma}^{\xi}$
(see Lemma 1(vi)):

$\displaystyle\underline{C}_{r}\leq\lim\inf_{n\to\infty}\frac{\mu_{r}^{x}}{n}\leq\lim\sup_{n\to\infty}\frac{\mu_{r}^{x}}{n}\leq\overline{C}_{r}\;\;\text{
and }\;\;\sup_{n\in\mathbb{N}}\mu_{r+1}^{x}\leq C^{\prime}_{\xi}+M_{\xi}.$

This property allows us to identify, asymptotically, as $n\to\infty$, the
number of factors $r$ (see, e.g., Bai and Ng, 2002, Onatski, 2010, Trapani,
2018, among many others). And, therefore, as $n\to\infty$ we also identify the
common and idiosyncratic components.

######  Assumption 4 (distinct eigenvalues).

For all $n\in\mathbb{N}$ and all $i=1,\ldots,r-1$,
$\mu_{i}^{\chi}>\mu_{i+1}^{\chi}$.

This assumption is needed in order to identify the eigenvectors in PC
estimation. Note that it implies that $\overline{C}_{j}<\underline{C}_{j-1}$,
$j=2,\ldots,r$, in (3) and that the $r$ eigenvalues of
$\bm{\Sigma}_{\Lambda}\bm{\Gamma}^{F}$ are distinct (Bai, 2003, Assumption G).
Indeed, these coincide with the non-zero eigenvalues of
$\lim_{n\to\infty}n^{-1}{\bm{\Gamma}^{\chi}}$ which are given by
$\lim_{n\to\infty}n^{-1}{\mathbf{M}^{\chi}}=(\bm{\Sigma}_{\Lambda})^{1/2}\bm{\Gamma}^{F}(\bm{\Sigma}_{\Lambda})^{1/2}$.

The factors and the loadings can be identified by means of the following
assumption.

######  Assumption 5 (Identification).

(a) For all $n\in\mathbb{N}$, $n^{-1}{\bm{\Lambda}^{\prime}\bm{\Lambda}}$
diagonal.  (b) For all $T\in\mathbb{N}$,
$T^{-1}{\bm{F}^{\prime}\bm{F}}=\mathbf{I}_{r}$.

This assumption a standard requirement in PC based exploratory factor analysis
(see, e.g, Bai and Ng, 2013, identification constraints PC1). It has some
important and usueful implications in terms of identification. In particular,
in Proposition B.1 we prove that, for all $n\in\mathbb{N}$,
$n^{-1}{\bm{\Lambda}^{\prime}\bm{\Lambda}}=n^{-1}{\mathbf{M}^{\chi}}$ and
$\bm{\Lambda}=\mathbf{V}^{\chi}(\mathbf{M}^{\chi})^{1/2}$, where
$\mathbf{M}^{\chi}$ is the $r\times r$ diagonal matrix of the nonzero
eigenvalues of $\bm{\Gamma}^{\chi}$ (sorted in decreasing order) and
$\mathbf{V}^{\chi}$ is the $n\times r$ matrix having as columns the
corresponding normalized eigenvectors. Notice that part (a) differs from Bai
and Li (2016, Condition IC3), where it is assumed that
$n^{-1}{\bm{\Lambda}^{\prime}(\bm{\Sigma}^{\xi})^{-1}\bm{\Lambda}}$ is
diagonal, for all $n\in\mathbb{N}$, with $\bm{\Sigma}^{\xi}$ being the
diagonal matrix with entries the idiosyncratic variances $\sigma_{i}^{2}$,
$i=1,\ldots,n$. This is a common requirement in QML based exploratory factor
analysis. One of the aims of this paper is to reconsider QML estimation under
Assumption 5(a). Part (b) is instead also assumed in Bai and Li (2016,
Condition IC3) but for deterministic factors.

Clearly since Assumption 5 concerns only the product
$\bm{\Lambda}^{\prime}\bm{\Lambda}$, it allows us to identify the loadings
only up to right multiplication by a diagonal matrix with entries $\pm 1$,
i.e., the columns of $\bm{\Lambda}$ are identified only up to a sign. We can
fix such sign by means of the following assumption.

######  Assumption 6 (Global identification).

For all $j=1,\ldots,r$, one of the two following conditions holds:  (a)
$\lambda_{1j}>0$; or  (b) $F_{j1}>0$.

Since loadings are considered as deterministic and are estimated as
eigenvectors, part (a) is more natural and can be easily imposed just by
setting the sign of the eigenvectors of the sample covariance matrix
$\widehat{\bm{\Gamma}}^{x}$ accordingly. As a consequence of Assumption 6 both
the loadings and the factors are globally identified (see also the comments in
Bai and Ng, 2013, Remark 1).

Finally, in order to derive the asymptotic distribution of the considered
estimators of the loadings it is common to assume the following assumption
(Bai, 2003, Assumption F.4, and Bai and Li, 2016, Assumption F.1).

######  Assumption 7 (Central limit theorem).

For all $i\in\mathbb{N}$, as $T\to\infty$,

$\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}\to_{d}\mathcal{N}\left(\mathbf{0}_{r},\bm{\Phi}_{i}\right),$

where
$\bm{\Phi}_{i}=\lim_{T\to\infty}T^{-1}\sum_{t,s=1}^{T}\mathbb{E}[\mathbf{F}_{t}\mathbf{F}_{s}^{\prime}\xi_{it}\xi_{is}]$.

There are many ways to derive this condition from more primitive assumptions,
for example, by assuming strong mixing factors and idiosyncratic components
both with finite $4+\epsilon$ moments.

######  Remark 1.

As discussed above, all assumptions for PC estimation by Bai (2003) are
implied or equivalent to ours. Regarding QML estimation and the assumptions in
Bai and Li (2016), we do not need their Assumption D, which requires the QML
estimators of the idiosyncratic variances, $\sigma_{i}^{2}$, $i=1,\ldots,n$,
to be finite and strictly positive, and we do not require the moment condition
E.3, which is needed only for estimation of $\sigma_{i}^{2}$ and it is not
necessary to prove our results. Assumptions E.4, E.5, E.6, F.2, and F.3 in Bai
and Li (2016) are needed only for studying the behavior of the estimated
factors so are not needed here. All other assumptions in Bai and Li (2016) are
equivalent or nested into ours, with the crucial exception of their
identifying condition on the loadings which, as noted above, differs from
ours.

######  Remark 2.

Allowing for data with non-zero mean simply amounts to adding a constant term,
$\alpha_{i}\neq 0$, to (1), so that
$x_{it}=\alpha_{i}+\bm{\lambda}_{i}^{\prime}\mathbf{F}_{t}+\xi_{it}$. All
estimators described in the following retain the same properties even in this
case. Indeed, estimation of such model by PCs simply requires to center data
first, i.e., to work with $x_{it}-\bar{x}_{i}$ where
$\bar{x}_{i}=T^{-1}\sum_{t=1}^{T}x_{it}$ is clearly a consistent estimator of
$\alpha_{i}$. As for QML estimation, it is straightforward to see that
$\bar{x}_{i}$ is precisely the QML estimator of $\alpha_{i}$ and thus it is
enough to work with the likelihood of the centered data (Bai and Li, 2012,
Section 2).

##  3 Principal Component Analysis

The PC estimators of $\bm{\Lambda}$ and $\bm{F}$ are the solutions of the
following minimization:

$\displaystyle\min_{\underline{\bm{\Lambda}},\underline{\bm{F}}}\frac{1}{nT}\text{tr}\left\\{\left(\bm{X}-\underline{\bm{F}}\,\underline{\bm{\Lambda}}^{\prime}\right)\left(\bm{X}-\underline{\bm{F}}\,\underline{\bm{\Lambda}}^{\prime}\right)^{\prime}\right\\}=\min_{\underline{\bm{\Lambda}},\underline{\bm{F}}}\frac{1}{nT}\text{tr}\left\\{\left(\bm{X}-\underline{\bm{F}}\,\underline{\bm{\Lambda}}^{\prime}\right)^{\prime}\left(\bm{X}-\underline{\bm{F}}\,\underline{\bm{\Lambda}}^{\prime}\right)\right\\},$
(5)

where $\underline{\bm{\Lambda}}$ and $\underline{\bm{F}}$ indicate generic
values of the loadings and the factors, respectively, and it is intended that
they also satisfy Assumptions 1, 4, and 5. Given the identifying constraints
in Assumption 5, the solution to (5) can be found in two steps. There are two
equivalent ways to do that:  (A) based on the $n\times n$ matrix
$\bm{X}^{\prime}\bm{X}$ solve first for loadings and then get the factors by
projecting $\bm{X}$ onto the estimated loadings;  (B) based on the $T\times T$
matrix $\bm{X}\bm{X}^{\prime}$ solve first for the factors and then get the
loadings by projecting $\bm{X}$ onto the estimated factors.

Although the majority of the literature on factor models considers approach B,
thus estimating the factors as normalized eigenvectors, and it derives the
theory accordingly (see, e.g., Bai, 2003), in the rest of the paper we follow
approach A. The main reason for this choice is that approach A does not
require to estimate the factors first, which is convenient given that our
focus is on estimation of the loadings. Notice that approach A is also the
classical one (see, e.g., Lawley and Maxwell, 1971, Chapter 4, Mardia et al.,
1979, Chapter 9.3, Jolliffe, 2002, Chapter 7.2). Nevertheless, which approach
to choose is just a matter of taste and it has no theoretical or practical
implications. Indeed, numerically both approaches give the same results and
all the following theory can be equivalently derived under approach B.
Incidentally, by choosing approach A, we also contribute to PC literature with
new proofs, alternative to those in Bai (2003).

More in detail, consider the $n\times n$ sample covariance matrix (recall that
$\mathbb{E}[\mathbf{x}_{t}]=\mathbf{0}_{n}$ by assumption)
$\widehat{\bm{\Gamma}}^{x}=T^{-1}{\bm{X}^{\prime}\bm{X}}$, having its $r$
largest eigenvalues collected in the $r\times r$ diagonal matrix
$\widehat{\mathbf{M}}^{x}$ (sorted in descending order) with corresponding
normalized eigenvectors as columns of the $n\times r$ matrix
$\widehat{\mathbf{V}}^{x}$. For any given $\underline{\bm{\Lambda}}$ the
solution of (5) for $\underline{\bm{F}}$ is just the linear projection
$\underline{\bm{F}}=\bm{X}\underline{\bm{\Lambda}}(\underline{\bm{\Lambda}}^{\prime}\underline{\bm{\Lambda}})^{-1}$.
By substituting this expression in (5) we have

$\displaystyle\min_{\underline{\bm{\Lambda}}}\frac{1}{nT}\text{tr}\left\\{\bm{X}\left(\mathbf{I}_{n}-\underline{\bm{\Lambda}}(\underline{\bm{\Lambda}}^{\prime}\underline{\bm{\Lambda}})^{-1}\underline{\bm{\Lambda}}^{\prime}\right)\bm{X}^{\prime}\right\\}$
$\displaystyle=\max_{\underline{\bm{\Lambda}}}\frac{1}{nT}\text{tr}\left\\{\bm{X}^{\prime}\bm{X}\underline{\bm{\Lambda}}(\underline{\bm{\Lambda}}^{\prime}\underline{\bm{\Lambda}})^{-1}\underline{\bm{\Lambda}}^{\prime}\right\\}$
$\displaystyle=\max_{\underline{\bm{\Lambda}}}\frac{1}{n}\text{tr}\left\\{(\underline{\bm{\Lambda}}^{\prime}\underline{\bm{\Lambda}})^{-1/2}\underline{\bm{\Lambda}}^{\prime}\frac{\bm{X}^{\prime}\bm{X}}{T}\underline{\bm{\Lambda}}(\underline{\bm{\Lambda}}^{\prime}\underline{\bm{\Lambda}})^{-1/2}\right\\}$
(6)

Now since by construction each column of
$\underline{\bm{\Lambda}}(\underline{\bm{\Lambda}}^{\prime}\underline{\bm{\Lambda}})^{-1/2}$
is normalized (since we assumed
${\underline{\bm{\Lambda}}^{\prime}\underline{\bm{\Lambda}}}$ to be diagonal),
then the above maximizaton, once solved, should return the $r$ largest
eigenvalues of $\widehat{\bm{\Gamma}}^{x}$ divided by $n$, i.e., it must give
$n^{-1}{\widehat{\mathbf{M}}^{x}}$. In other words, our estimator
$\widehat{\bm{\Lambda}}$ must be such that
$\widehat{\bm{\Lambda}}(\widehat{\bm{\Lambda}}^{\prime}\widehat{\bm{\Lambda}})^{-1/2}$
is the matrix of normalized eigenvectors corresponding the $r$ largest
eigenvalues of $(nT)^{-1}{\bm{X}^{\prime}\bm{X}}$, i.e., such that:

$(\widehat{\bm{\Lambda}}^{\prime}\widehat{\bm{\Lambda}})^{-1/2}\widehat{\bm{\Lambda}}^{\prime}\frac{\bm{X}^{\prime}\bm{X}}{nT}\widehat{\bm{\Lambda}}(\widehat{\bm{\Lambda}}^{\prime}\widehat{\bm{\Lambda}})^{-1/2}=\frac{\widehat{\mathbf{M}}^{x}}{n},$
(7)

but also, by definition of eigenvectors,

$\widehat{\mathbf{V}}^{x\prime}\frac{\bm{X}^{\prime}\bm{X}}{nT}\widehat{\mathbf{V}}^{x}=\frac{\widehat{\mathbf{M}}^{x}}{n}.$
(8)

Therefore, since we must have $\bm{\Lambda}^{\prime}\bm{\Lambda}$ diagonal
with distinct entries by Assumption 5(a), from (7) and (8),

$\widehat{\bm{\Lambda}}=\widehat{\mathbf{V}}^{x}(\widehat{\mathbf{M}}^{x})^{1/2},$
(9)

which is such that
$n^{-1}{\widehat{\bm{\Lambda}}^{\prime}\widehat{\bm{\Lambda}}}=n^{-1}{\widehat{\mathbf{M}}^{x}}$
is diagonal. The factors are then estimated as the linear projections:
$\widehat{\bm{F}}=\bm{X}\widehat{\bm{\Lambda}}(\widehat{\bm{\Lambda}}^{\prime}\widehat{\bm{\Lambda}})^{-1}=\bm{X}\widehat{\mathbf{V}}^{x}(\widehat{\mathbf{M}}^{x})^{-1/2}$,
which are the normalized PCs of $\bm{X}$, such that
$T^{-1}{\widehat{\bm{F}}^{\prime}\widehat{\bm{F}}}=\mathbf{I}_{r}$. The
latter, however, are not needed in the following.

Letting $\widehat{\bm{\lambda}}_{i}^{\prime}$ be the $i$th row of the PC
estimator of the loadings in (9), we have the following asymptotic results.

######  Theorem 1.

Under Assumptions 1 through 7, as $n,T\to\infty$,

  1. (a)

$\min(n,\sqrt{T})\|\widehat{\bm{\lambda}}_{i}-{\bm{\lambda}}_{i}\|=O_{\mathrm{P}}(1)$,
for any given $i=1,\ldots,n$, and, if $n^{-1}\sqrt{T}\to 0$,
$\sqrt{T}(\widehat{\bm{\lambda}}_{i}-{\bm{\lambda}}_{i})\to_{d}\mathcal{N}\left(\mathbf{0}_{r},\bm{\Phi}_{i}\right),$
where $\bm{\Phi}_{i}$ is defined in Assumption 7;

  2. (b)

$\min(n,\sqrt{T})\|n^{-1/2}({\widehat{\bm{\Lambda}}-\bm{\Lambda}})\|=O_{\mathrm{P}}(1)$.

An important consequence of Theorem 1 is that the PC estimator is
asymptotically equivalent to the unfeasible OLS estimator, which we would
obtain if the factors were observed, and denoted as $\bm{\Lambda}^{\text{\tiny
OLS}}$, with $i$th row given by ${\bm{\lambda}}_{i}^{\text{\tiny OLS}\prime}$.

######  Corollary 1.

Under Assumptions 1 through 7, as $n,T\to\infty$,  (a)
$\min(n,\sqrt{nT})\|\widehat{\bm{\lambda}}_{i}-{\bm{\lambda}}_{i}^{\text{\tiny
OLS}}\|=O_{\mathrm{P}}\left(1\right)$, for any given $i=1,\ldots,n$;  (b)
$\min(n,\sqrt{nT})\|n^{-1/2}({\widehat{\bm{\Lambda}}-{\bm{\Lambda}}^{\text{\tiny
OLS}}})\|=O_{\mathrm{P}}(1)$.

##  4 Quasi Maximum Likelihood

Define
$\bm{\mathcal{X}}=\text{vec}(\bm{X}^{\prime})=(\mathbf{x}_{1}^{\prime}\cdots\mathbf{x}_{T}^{\prime})^{\prime}$
and
$\bm{\mathcal{Z}}=\text{vec}(\bm{\Xi}^{\prime})=(\bm{\xi}_{1}^{\prime}\cdots\bm{\xi}_{T}^{\prime})^{\prime}$
as the $nT$-dimensional vectors of observations and idiosyncratic components,
$\bm{\mathcal{F}}=\text{vec}(\bm{F}^{\prime})=(\mathbf{F}_{1}^{\prime}\cdots\mathbf{F}_{T}^{\prime})^{\prime}$
as the $rT$-dimensional vector of factors, and
$\bm{\mathfrak{L}}=\mathbf{I}_{T}\otimes\bm{\Lambda}$ as the $nT\times rT$
matrix containing all factor loadings replicated $T$ times. Then, by
vectorizing the transposed of (2) we have:

$\displaystyle\bm{\mathcal{X}}=\bm{\mathfrak{L}}\bm{\mathcal{F}}+\bm{\mathcal{Z}}.$
(10)

Let $\bm{\Omega}^{x}=\mathbb{E}[\bm{\mathcal{X}}\bm{\mathcal{X}}^{\prime}]$
and $\bm{\Omega}^{\xi}=\mathbb{E}[\bm{\mathcal{Z}}\bm{\mathcal{Z}}^{\prime}]$,
be the $nT\times nT$ covariance matrices of the $nT$-dimensional vectors of
data and idiosyncratic components, respectively. Let also
$\bm{\Omega}^{F}=\mathbb{E}[\bm{\mathcal{F}}\bm{\mathcal{F}}^{\prime}]$ be the
$rT\times rT$ covariance matrix of the $rT$-dimensional factor vector. Then,
because of Assumption 3,

$\bm{\Omega}^{x}=\bm{\mathfrak{L}}\,\bm{\Omega}^{F}\bm{\mathfrak{L}}^{\prime}+\bm{\Omega}^{\xi}.$

In principle, the parameters that need to be estimated are then given by the
vector

$\bm{\varphi}=(\mathrm{vec}(\bm{\Lambda})^{\prime},\mathrm{vech}(\bm{\Omega}^{\xi})^{\prime},\mathrm{vech}(\bm{\Omega}^{F})^{\prime})^{\prime},$

and the Gaussian quasi-log-likelihood computed in a generic value of the
parameters, denoted as $\underline{\bm{\varphi}}$, is given by (omitting the
constant term for simplicity)

$\displaystyle\ell^{*}(\bm{\mathcal{X}};\underline{\bm{\varphi}})$
$\displaystyle=-\frac{1}{2}\log\det\left(\underline{\bm{\Omega}}^{x}\right)-\frac{1}{2}\bm{\mathcal{X}}^{\prime}\left(\underline{\bm{\Omega}}^{x}\right)^{-1}\bm{\mathcal{X}}.$
(11)

In general, maximization of (11) is an unfeasible task since the parameters
vector ${\bm{\varphi}}$ to be estimated has $\simeq(nT)^{2}$ elements. The
common practice is then to consider simpler mis-specified log-likelihoods
which depend on fewer parameters (Bai and Li, 2012, 2016).

First of all, consistently with the fact that we do not assume any parametric
model describing the dynamics of the factors, we can consider a simpler log-
likelihood with $\bm{\Omega}^{F}=\mathbf{I}_{nr}$, where we imposed also the
identification constraint $\bm{\Gamma}^{F}=\mathbf{I}_{r}$ implied by
Assumption 5(b). A second simplification consists in considering a mis-
specified log-likelihood of an approximate factor model where also the
idiosyncratic components are treated as serially uncorrelated. As a result of
these two mis-specifications the log-likelihood (11) is reduced to:

$\ell(\bm{\mathcal{X}};\underline{\bm{\varphi}})=-\frac{T}{2}\log\det(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Gamma}}^{\xi})-\frac{1}{2}\sum_{t=1}^{T}\mathbf{x}_{t}^{\prime}(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Gamma}}^{\xi})^{-1}\mathbf{x}_{t}.$
(12)

The parameters to be estimated are then reduced:
$\bm{\varphi}=(\mathrm{vec}(\bm{\Lambda})^{\prime},\mathrm{vech}(\bm{\Gamma}^{\xi})^{\prime})^{\prime}$.
Nevertheless, estimation of $\bm{\varphi}$ by means of maximization of (12)
seems still hopeless since in general ${\bm{\Gamma}}^{\xi}$ contains
$n(n+1)/2$ distinct elements.

Some further mis-specification of the log-likelihood (12), based on
regularizing $\bm{\Gamma}^{\xi}$, is then usually introduced in order to
reduce the number of parameters to be estimated. A possibility in this sense
is explored, for example, by Bai and Liao (2016) who propose to maximize (12)
subject to an $\ell_{1}$ penalty imposed on the off-diagonal entries of
$\bm{\Gamma}^{\xi}$. This approach forces sparsity, thus reducing the number
of parameters to be estimated, but at the same time it makes estimation
dependent on the chosen penalization level, which affects also the rate of
consistency.

An even simpler approach consists in estimating only the diagonal entries of
$\bm{\Gamma}^{\xi}$. Specifically, letting
$\bm{\Sigma}^{\xi}=\text{dg}(\sigma_{1}^{2}\cdots\sigma_{n}^{2})$ be the
diagonal matrix with entries the diagonal entries of $\bm{\Gamma}^{\xi}$, we
can focus on maximization of the further mis-specified log-likelihood:

$\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})=-\frac{T}{2}\log\det(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Sigma}}^{\xi})-\frac{1}{2}\sum_{t=1}^{T}\mathbf{x}_{t}^{\prime}(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Sigma}}^{\xi})^{-1}\mathbf{x}_{t}.$
(13)

The parameters to be estimated are then reduced to:
$\bm{\varphi}=(\mathrm{vec}(\bm{\Lambda})^{\prime},\sigma_{1}^{2},\cdots,\sigma_{n}^{2})^{\prime}$,
which are just $nr+n$. This is a feasible task given that we have $nT$ data
points.

The log-likelihood (13) is the one considered in classical factor analysis,
where, however, $n$ is assumed to be fixed and small, and, moreover, the true
idiosyncratic covariance matrix is assumed to be diagonal, i.e., the
likelihood is not mis-specified (see, e.g., Lawley and Maxwell, 1971, Chapter
2, and Rubin and Thayer, 1982). In the high-dimensional case, i.e., when we
allow $n\to\infty$, maximization of (13) has been studied by Bai and Li (2012,
2016) under a variety of possible identifying constraints. In particular,
while Bai and Li (2012) consider the case of no idiosyncratic serial or cross-
correlation, i.e., they assume
$\bm{\Omega}^{\xi}=\mathbf{I}_{T}\otimes\bm{\Sigma}^{\xi}$, which is diagonal,
thus considering (13) as a correctly specified log-likelihood, Bai and Li
(2016) allow instead for idiosyncratic serial and cross-sectional
correlations, as we do in this paper, thus treating (13) as a mis-specified
log-likelihood.

In the rest of this section we briefly review the properties of the QML
estimator of the loadings. We refer to Barigozzi (2024) for a full review of
QML estimation of factor models. We start with the simplest case in which we
consider a further mis-specification of the log-likelihood (13), where we
treat the idiosyncratic components as if they were homoskedastic, thus in the
log-likelihood we replace $\bm{\Sigma}^{\xi}$ with an even simpler covariance
matrix ${\sigma}^{2}\mathbf{I}_{n}$, with ${\sigma}^{2}>0$ and finite. In this
case Tipping and Bishop (1999) prove that the QML estimator of the loadings
matrix $\bm{\Lambda}$, has a closed form given by

$\widehat{\bm{\Lambda}}^{\text{\tiny
QML,E}_{0}}=\widehat{\mathbf{V}}^{x}\left(\widehat{\mathbf{M}}^{x}-\widehat{\sigma}^{2\text{\tiny
QML,E}_{0}}\mathbf{I}_{r}\right)^{1/2}=\widehat{\mathbf{V}}^{x}\left(\widehat{\mathbf{M}}^{x}-\left\\{\frac{1}{n-r}\sum_{j=r+1}^{n}\widehat{\mu}^{x}_{j}\right\\}\mathbf{I}_{r}\right)^{1/2},$
(14)

with $\widehat{\sigma}^{2\text{\tiny QML,E}_{0}}$ being the QML estimator of
$\sigma^{2}$. Intuitively, since under our assumptions we should have
$\widehat{\mathbf{M}}^{x}=O_{\mathrm{P}}(n)$, as $n\to\infty$, then
$\widehat{\bm{\Lambda}}^{\text{\tiny QML,E}_{0}}$ seems to coincide
asymptotically with the PC estimator given in (9). This is a well known fact
and it is often quoted in the literature (see, e.g., Doz et al., 2012): in the
case of spherical idiosyncratic components the PC and QML estimators are
asymptotically equivalent. However, to the best of our knowledge no formal
proof exists, at least under the present set of assumptions. In fact, the
proof would essentially require to prove that the estimated $(r+1)$th largest
eigenvalue is such that $\widehat{\mu}^{x}_{r+1}=O_{\mathrm{P}}(1)$, which is
not an easy task in a high-dimensional setting, because, although we know that
under our assumptions ${\mu}^{x}_{r+1}=O(1)$ (see Lemma 1(v)), in general
$\widehat{\mu}^{x}_{r+1}$ is not a consistent estimator of ${\mu}^{x}_{r+1}$
(see, e.g., Trapani, 2018, Lemma 2.2).

Things become more complicated if we allow for heteroskedasticity. Indeed, no
closed form solution exists for the QML estimator maximizing the mis-specified
log-likelihood (13) and numerical maximization is required instead (see, e.g.,
Bai and Li, 2012, Section 8, for a proposed algorithm). Still it is possible
to derive its asymptotic properties. Let us denote as
$\widehat{\bm{\lambda}}_{i}^{\text{\tiny QML,E}}$, $i=1,\ldots,n$, the QML
estimator of the $i$th row of $\bm{\Lambda}$. In the classical fixed $n$ case,
such estimator retains the classical propertied of the QML estimators, so it
is $\sqrt{T}$-consistent and asymptotically normal. However, since the first
and second derivatives of (13) are very complex, the asymptotic covariance
matrix has also a very complicated form, which in turn makes its estimation
not at all easy (see, e.g., Anderson and Rubin, 1956, Theorem 12.3, Amemiya et
al., 1987, Theorem 2F, and Anderson and Amemiya, 1988, Theorems 1, 2, and 3).
If, instead, we study the properties of the QML estimator when allowing also
for $n\to\infty$, things become, perhaps surprisingly, simpler. This is shown
in the following theorem proved by Bai and Li (2016, Theorem 1).

######  Theorem 2.

Assume $n^{-1}{\bm{\Lambda}^{\prime}(\bm{\Sigma}^{\xi})^{-1}\bm{\Lambda}}$ to
be diagonal for all $n\in\mathbb{N}$ and
$T^{-1}{\bm{F}^{\prime}\bm{F}}=\mathbf{I}_{r}$ for all $T\in\mathbb{N}$. Then,
under Assumptions 1, 2, 3, and 7, if $n^{-1}\sqrt{T}\to 0$, as $n,T\to\infty$,
for any given $i=1,\ldots,n$,
$\sqrt{T}(\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML,E}}-\bm{\lambda}_{i})\to_{d}\mathcal{N}\left(\mathbf{0}_{r},\bm{\Phi}_{i}\right),$
where $\bm{\Phi}_{i}$ is defined in Assumption 7.

The proof of this result is based on asymptotic expansions of a set of
conditions derived from first order conditions computed for the log-likelihood
(13). It is a very long proof, aimed at showing that
$\widehat{\bm{\lambda}}_{i}^{\text{\tiny QML,E}}$ is asymptotically equivalent
to the unfeasible OLS we would obtain if we knew the factors, which in turn
has a known asymptotic distribution (see the online supplement of Bai and Li,
2016). The proof requires a series of technical assumptions different from,
but nested into, ours (see Remark 1). However, as noticed already in Section
2, it is important to stress that Theorem 2 holds under an identification
constraint for the loadings which differs from our Assumption 5(a). At the end
of the next section we prove that this theorem holds also under our
assumptions (see Theorem 4 and Corollary 2).

######  Remark 3.

Notice that the fact that we consider a mis-specified log-likelihood with
serially and cross-sectionally uncorrelated idiosyncratic components, does not
affect consistency of the estimated loadings, but only their asymptotic
covariance. And in particular, notice that what matters for the asymptotic
covariance is just the idiosyncratic serial correlation. Indeed, if the
idiosyncratic components were cross-sectionally uncorrelated but serially
correlated the asymptotic covariance in Theorem 2 would be the same, while if
they were serially uncorrelated the asymptotic covariance would be
$\sigma_{i}^{2}\mathbf{I}_{r}$ (due to the identification
$\bm{\Gamma}^{F}=\mathbf{I}_{r}$), regardless of the presence or not of cross-
sectional correlation. This result, which is a special case of Theorem 2, does
not require any constraint between the rates of divergence of $n$ and $T$ and
it is proved by Bai and Li (2012, Theorem 5.2).

######  Remark 4.

The case of QML estimation when $\bm{\Omega}^{F}$ depends explicitly on
additional parameters capturing the autocorrelations in the factors is
considered in Doz et al. (2012). However, in that case QML estimation requires
the use of the EM algorithm jointly with the Kalman smoother and the results
of this paper do not directly apply unless we first prove convergence of the
EM to the same QML estimator considered here.

##  5 The PC and QML estimators are asymptotically equivalent

Given the discussion in the previous section, we might argue that the PC and
QML estimators have the same asymptotic properties. However, Theorems 1 and 2
are derived under a similar but different sets of identifying assumptions and
the two results cannot be directly compared. Even in the spherical case the
proof seems to be not so easy due to the unknown properties of the smallest
$N-r$ sample eigenvalues.

It is then natural to ask the following question. Can we prove in a simple way
that the PC and QML estimators are asymptotically equivalent under the same
set of assumptions given in this paper, which are the standard PC assumptions?
Moreover, if we prove such equivalence and since the PC estimator of the
loadings does not depend on the idiosyncratic covariance matrix, it is natural
to ask also the following additional question. Can we prove the asymptotic
equivalence of the two estimators when considering QML based on the log-
likelihood (12) of an approximate factor model, i.e., without constraining the
idiosyncratic covariance to be diagonal or even homoskedastic? In this section
we answer both questions.

Denote as $\widehat{\bm{\Lambda}}^{\text{\tiny QML}}$ the QML estimator of the
loadings matrix maximizing the mis-specified log-likleihood (12), where the
idiosyncratic components are treated as serially uncorrelated, but their
covariance matrix is unrestricted so it is correctly specified. Let also
$\widehat{\bm{\lambda}}_{i}^{\text{\tiny QML}\prime}$, $i=1,\ldots,n$, be the
$i$th row of $\widehat{\bm{\Lambda}}^{\text{\tiny QML}}$. Then, we state our
main result.

######  Theorem 3.

Under Assumptions 1 through 6 and assuming also that $\bm{\Gamma}^{\xi}$ is
positive definite, as $n,T\to\infty$,  (a)
$n\|n^{-1/2}({\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}-\widehat{\bm{\Lambda}}})\|=O_{\mathrm{P}}(1)$;  (b)
$n\|\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML}}-\widehat{\bm{\lambda}}_{i}\|=O_{\mathrm{P}}(1)$, for any given
$i=1,\ldots,n$.

Consistency and asymptotic normality of the QML estimator of the loadings
maximizing the log-likelihood (12) immediately follow.

######  Theorem 4.

Under Assumptions 1 through 6 and assuming also that $\bm{\Gamma}^{\xi}$ is
positive definite, as $n,T\to\infty$:

  1. (a)

$\min(n,\sqrt{T})\|\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML}}-{\bm{\lambda}}_{i}\|=O_{\mathrm{P}}(1)$, for any given $i=1,\ldots,n$,,
and, if $n^{-1}\sqrt{T}\to 0$
$\sqrt{T}(\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML}}-\bm{\lambda}_{i})\to_{d}\mathcal{N}\left(\mathbf{0}_{r},\bm{\Phi}_{i}\right),$
where $\bm{\Phi}_{i}$ is defined in Assumption 7;

  2. (b)

$\min(n,\sqrt{T})\|n^{-1/2}({\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}-\bm{\Lambda}})\|=O_{\mathrm{P}}(1)$.

The results in Theorems 3 and 4 hold when considering the log-likelihood (12)
which depends on any generic idiosyncratic covariance matrix provided that it
is positive definite and Assumption 2 is satisfied. Therefore, from Theorem 3
it follows that if we replace in the log-likelihood the full idiosyncratic
covariance with a diagonal covariance matrix as in Bai and Li (2012, 2016) or
if we also impose homoskedasticity as in Tipping and Bishop (1999), we still
get QML estimators that are asymptotically equivalent to the PC estimator.

######  Corollary 2.

Under Assumptions 1 through 6, as $n,T\to\infty$,  (a)
$n\|\widehat{\bm{\lambda}}_{i}-\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML,E}}\|=O_{\mathrm{P}}(1)$;  (b)
$n\|\widehat{\bm{\lambda}}_{i}-\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML,E}_{0}}\|=O_{\mathrm{P}}(1)$.

From these results and Theorem 4, we directly have a proof of Theorem 2 by Bai
and Li (2016), which now holds under our identifying constraint in Assumption
5(a), and we also prove the often quoted statement that the PC and QML
estimators are asymptotically equivalent under sphericity of idiosyncratic
components.

######  Remark 5.

The appealing feature of the log-likelihood (12), is that it does not depend
on the factors, and, under the identifying constraint
$T^{-1}{\bm{F}^{\prime}\bm{F}}=\mathbf{I}_{r}$ for all $T\in\mathbb{N}$, it
does not depend on the second moments of the factors either. Still, if we
follow the classical approach to consider the vector of factors
$\bm{\mathcal{F}}$ as an $rT$-dimensional sequence of deterministic constants
(Lawley and Maxwell, 1971), then a log-likelihood alternative to (13) can be
considered, namely:

$\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}},\underline{\bm{\mathcal{F}}})=-\frac{T}{2}\log\det(\underline{\bm{\Gamma}}^{\xi})-\frac{1}{2}\sum_{t=1}^{T}(\mathbf{x}_{t}-\underline{\bm{\Lambda}}\,\underline{\mathbf{F}}_{t})^{\prime}(\underline{\bm{\Gamma}}^{\xi})^{-1}(\mathbf{x}_{t}-\underline{\bm{\Lambda}}\,\underline{\mathbf{F}}_{t}).$
(15)

It is straightforward to see that, for given factors, the loadings maximizing
(15) are given by their OLS estimator, while, for given loadings, the factors
maximizing (15) are given by their GLS estimator. Hence, full maximization of
(15) requires knowing, or estimating, the factors too. Although seemingly
simpler this approach presents at least two major drawbacks with respect to
the approach followed in the proof of Theorem 3 (see also the comments in Bai
and Li, 2012, p.440, and Anderson and Rubin, 1956, pp.129-130).

First, by iterating, between OLS and GLS we can think of finding a solution.
This approach is similar to the one recently considered by Zadrozny (2023),
but the convergence properties of such algorithm are not proved nor discussed,
so it is unclear how this estimator is related to the maximizer of the full
likelihood (12). Moreover, this strategy would require a positive definite
estimator of $\bm{\Gamma}^{\xi}$ in order to compute the GLS estimator of the
factors, a hard task when $n$ is large. This, in general, requires again mis-
specifying or regularizing the log-likelihood (15), e.g., by replacing
$\bm{\Gamma}^{\xi}$ with the diagonal $\bm{\Sigma}^{\xi}$. The asymptotic
properties of the estimator of the loadings will then explicitly depend on the
properties of an estimator of the idiosyncratic covariance, or at least of its
diagonal elements. As shown in Theorem 4, this is not the case for the QML
estimator maximizing the log-likelihood (12).

Second, if the factors are treated as random variables, as, e.g., in the
popular Factor Augmented VAR models (Bernanke et al., 2005), then they cannot
be considered as constant parameters, which means that (15) is not the full
log-likelihood of the data $\bm{\mathcal{X}}$ but it is just the conditional
log-likelihood of $\bm{\mathcal{X}}$ given the factors, so, in principle, not
all information is used to estimate the loadings. Indeed, if the factors are
random variables then the log-likelihood (13) is decomposed as

$\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})=\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})+\ell_{\text{\tiny
E}}(\bm{\mathcal{F}};\underline{\bm{\varphi}})-\ell_{\text{\tiny
E}}(\bm{\mathcal{F}}|\bm{\mathcal{X}};\underline{\bm{\varphi}}),$ (16)

where $\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})$ coincides
with (15) but it does not coincide with $\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})$ anymore.

This last point has both a theoretical and an applied implication. From the
theory point of view, to directly show from (16) that the QML estimator of the
loadings maximizing $\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})$ is asymptotically equivalent
to the unfeasible OLS maximizing $\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})$ would require
showing that $\ell_{\text{\tiny
E}}(\bm{\mathcal{F}};\underline{\bm{\varphi}})$ and $\ell_{\text{\tiny
E}}(\bm{\mathcal{F}}|\bm{\mathcal{X}};\underline{\bm{\varphi}})$ are
asymptotically negligible. This is the argument sketched by Breitung and
Tenhofen (2011), but to make it formal is not an easy task. Consider the
simplest case in which the factors are treated as serially uncorrelated, then,
while $\ell_{\text{\tiny E}}(\bm{\mathcal{F}};\underline{\bm{\varphi}})$ does
not depend on any parameter and can be discarded, the expression of
$\ell_{\text{\tiny
E}}(\bm{\mathcal{F}}|\bm{\mathcal{X}};\underline{\bm{\varphi}})$ will depend
on the conditional moments (mean and covariance) of the factors given
$\bm{\mathcal{X}}$. These in turn have simple expressions only if we are
willing to assume Gaussianity, in which case the conditional mean is just a
linear projection. Otherwise computation of those moments is not
straightforward. The proof of Theorem 3 relies instead only on
$\ell_{\text{\tiny E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})$, so, as
noticed above, it does not require knowing the factors or their conditional
moments.

Finally, from a practical point of view, we could use the right-hand-side of
(16) to compute the QML estimator by means of the Expectation Maximization
(EM) algorithm, which simplifies our task since it allows us to discard
$\ell_{\text{\tiny
E}}(\bm{\mathcal{F}}|\bm{\mathcal{X}};\underline{\bm{\varphi}})$ (Wu, 1983).
However, once again we would still have to compute the conditional moments of
the factors when in the E-step we need compute the conditional expectation of
$\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})+\ell_{\text{\tiny
E}}(\bm{\mathcal{F}};\underline{\bm{\varphi}})$ given $\bm{\mathcal{X}}$. This
introduces a correction in the estimation of the loadings and the resulting
estimator is not given by a simple OLS anymore.

##  6 Asymptotic covariance matrices of the QML and PC estimators

In this section, we focus on the mis-specified log-likelihood (13), which is
commonly used in empirical work (Bai and Li, 2012, 2016). For such log-
likelihood, denote the Fisher information and the population Hessian matrices
for $\bm{\lambda}_{i}$, $i=1,\ldots,n$, as:

$\displaystyle\bm{\mathcal{I}}_{i}(\bm{\mathcal{X}};\bm{\varphi})$
$\displaystyle=\lim_{T\to\infty}\mathbb{E}\left[\left(\frac{1}{\sqrt{T}}\left.\frac{\partial\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}}\right|_{\underline{\bm{\varphi}}={{\bm{\varphi}}}}\right)\left(\frac{1}{\sqrt{T}}\left.\frac{\partial\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}}\right|_{\underline{\bm{\varphi}}={{\bm{\varphi}}}}\right)\,\right],$
$\displaystyle\bm{\mathcal{H}}_{i}(\bm{\mathcal{X}};\bm{\varphi})$
$\displaystyle=\lim_{T\to\infty}\mathbb{E}\left[\frac{1}{T}\left.\frac{\partial^{2}\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}\partial\underline{\bm{\lambda}}_{i}}\right|_{\underline{\bm{\varphi}}={{\bm{\varphi}}}}\right].$

From QML theory the asymptotic covariance of the QML estimator should be given
by:

$\text{A$\\!\mathbb{V}$ar}_{0}(\sqrt{T}\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML,E}})=\\{\bm{\mathcal{H}}_{i}(\bm{\mathcal{X}};\bm{\varphi})\\}^{-1}\bm{\mathcal{I}}_{i}(\bm{\mathcal{X}};\bm{\varphi})\\{\bm{\mathcal{H}}_{i}(\bm{\mathcal{X}};\bm{\varphi})\\}^{-1}.$
(17)

This would be the matrix to estimate if we were to conduct QML based inference
on the loadings.

In general, estimation of (17) is very difficult given the complex expressions
of the Hessian and Fisher information matrices (see also (A.15) and (A.42) in
the Supplementary Material). Moreover, from Theorem 4 we know that, in fact,
if $n^{-1}\sqrt{T}\to 0$, as $n,T\to\infty$, then the asymptotic covariance of
the QML estimator is

$\text{A$\\!\mathbb{V}$ar}_{1}(\sqrt{T}\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML,E}})=\bm{\Phi}_{i}=\lim_{T\to\infty}\frac{1}{T}\sum_{s,t=1}^{T}\mathbb{E}\left[\mathbf{F}_{t}\mathbf{F}_{s}^{\prime}\xi_{it}\xi_{is}\right],$
(18)

where we used the definition of $\bm{\Phi}_{i}$ in Assumption 7. So what is
the relation between the asymptotic covariances (17) and (18)?

First, notice that (18) coincides with the asymptotic covariance of the
unfeasible OLS estimator, which, in turn, is the QML estimator maximizing the
log-likelihood for an exact factor model conditional on observing the factors,
i.e.,

$\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})=-\frac{T}{2}\log\det(\underline{\bm{\Sigma}}^{\xi})-\frac{1}{2}\sum_{t=1}^{T}(\mathbf{x}_{t}-\underline{\bm{\Lambda}}\mathbf{F}_{t})^{\prime}(\underline{\bm{\Sigma}}^{\xi})^{-1}(\mathbf{x}_{t}-\underline{\bm{\Lambda}}\mathbf{F}_{t}).$
(19)

It is easy to compute the Fisher information and the population Hessian
matrices of (19) for $\bm{\lambda}_{i}$: (see also (A.37) and (A.47) in the
Supplementary Material)

$\displaystyle\bm{\mathcal{I}}_{i}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\bm{\varphi})$
$\displaystyle=\lim_{T\to\infty}\mathbb{E}\left[\left(\frac{1}{\sqrt{T}}\left.\frac{\partial\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}}\right|_{\underline{\bm{\varphi}}={{\bm{\varphi}}}}\right)\left(\frac{1}{\sqrt{T}}\left.\frac{\partial\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}}\right|_{\underline{\bm{\varphi}}={{\bm{\varphi}}}}\right)\,\right]=\frac{1}{\sigma_{i}^{4}}\bm{\Phi}_{i},$
$\displaystyle\bm{\mathcal{H}}_{i}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\bm{\varphi})$
$\displaystyle=\lim_{T\to\infty}\mathbb{E}\left[\frac{1}{T}\left.\frac{\partial^{2}\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}\partial\underline{\bm{\lambda}}_{i}}\right|_{\underline{\bm{\varphi}}={{\bm{\varphi}}}}\right]=-\frac{1}{\sigma_{i}^{2}}\mathbf{I}_{r},$

where we imposed orthonormality of the factors as required by Assumption 5(b).
Clearly,

$\bm{\Phi}_{i}=\left\\{\bm{\mathcal{H}}_{i}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\bm{\varphi})\right\\}^{-1}\bm{\mathcal{I}}_{i}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\bm{\varphi})\left\\{\bm{\mathcal{H}}_{i}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\bm{\varphi})\right\\}^{-1}.$
(20)

It follows that for the asymptotic covariances (17) and (18) to be equivalent
it must be that the Fisher information and Hessian matrices of the full log-
likelihood $\ell_{\text{\tiny E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})$
in (12) and of the conditional log-likelihood $\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})$ in (19) are
asymptotically equivalent when computed in the true values of the parameters
$\bm{\varphi}$. This is proved in the following theorem.

######  Theorem 5.

Under Assumptions 1 through 6, as $n,T\to\infty$, for any given
$i=1,\ldots,n$,

  1. (a)

$\frac{1}{\sqrt{T}}\left\|\left.\frac{\partial\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}}\right|_{\underline{\bm{\varphi}}={\bm{\varphi}}}-\left.\frac{\partial\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}}\right|_{\underline{\bm{\varphi}}={\bm{\varphi}}}\right\|=O_{\mathrm{P}}\left(\max\left(\frac{1}{\sqrt{n}},\frac{\sqrt{T}}{n}\right)\right);$

  2. (b)

$\frac{1}{T}\left\|\left.\frac{\partial^{2}\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}\partial\underline{\bm{\lambda}}_{i}}\right|_{\underline{\bm{\varphi}}={\bm{\varphi}}}-\left.\frac{\partial^{2}\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}\partial\underline{\bm{\lambda}}_{i}}\right|_{\underline{\bm{\varphi}}={\bm{\varphi}}}\right\|=O_{\mathrm{P}}\left(\max\left(\frac{1}{n},\frac{1}{\sqrt{nT}}\right)\right);$

  3. (c)

$\frac{1}{\sqrt{T}}\left\|\left.\frac{\partial\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}}\right|_{\underline{\bm{\varphi}}=\widehat{\bm{\varphi}}^{\text{\tiny
QML,E}}}-\left.\frac{\partial\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}}\right|_{\underline{\bm{\varphi}}=\widehat{\bm{\varphi}}^{\text{\tiny
QML,E}}}\right\|=O_{\mathrm{P}}\left(\max\left(\frac{1}{\sqrt{n}},\frac{\sqrt{T}}{n}\right)\right);$

  4. (d)

$\frac{1}{T}\left\|\left.\frac{\partial^{2}\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}\partial\underline{\bm{\lambda}}_{i}}\right|_{\underline{\bm{\varphi}}=\widehat{\bm{\varphi}}^{\text{\tiny
QML,E}}}-\left.\frac{\partial^{2}\ell_{\text{\tiny
E}}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\lambda}}_{i}^{\prime}\partial\underline{\bm{\lambda}}_{i}}\right|_{\underline{\bm{\varphi}}=\widehat{\bm{\varphi}}^{\text{\tiny
QML,E}}}\right\|=O_{\mathrm{P}}\left(\max\left(\frac{1}{n},\frac{1}{\sqrt{nT}}\right)\right).$

From Theorem 5 it is clear that if $n^{-1}\sqrt{T}\to 0$, as $n,T\to\infty$,
then we can consistently estimate
$\text{A$\\!\mathbb{V}$ar}_{0}(\sqrt{T}\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML,E}})$ in (17) by means of any consistent estimator of
$\text{A$\\!\mathbb{V}$ar}_{1}(\sqrt{T}\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML,E}})=\bm{\Phi}_{i}$, as the classical HAC estimator:

$\displaystyle\widehat{\bm{\Phi}}_{i}=\left(\frac{1}{T}\sum_{t=1}^{T}\widehat{\mathbf{F}}_{t}\widehat{\mathbf{F}}_{t}^{\prime}\widehat{\xi}_{it}^{2}\right)$
$\displaystyle+\left\\{\sum_{k=1}^{m_{T}}\left(1-\frac{k}{M_{T}+1}\right)\left[\left(\frac{1}{T}\sum_{t=k+1}^{T}\widehat{\mathbf{F}}_{t}\widehat{\mathbf{F}}_{t-k}^{\prime}\widehat{\xi}_{it}\widehat{\xi}_{it-k}\right)\right.\right.$
$\displaystyle+\left.\left.\left(\frac{1}{T}\sum_{t=k+1}^{T}\widehat{\mathbf{F}}_{t}\widehat{\mathbf{F}}_{t-k}^{\prime}\widehat{\xi}_{it}\widehat{\xi}_{it-k}\right)^{\prime}\right]\right\\},$

where
$\widehat{\mathbf{F}}_{t}=(\widehat{\mathbf{M}}^{x})^{-1/2}\widehat{\mathbf{V}}^{x\prime}\mathbf{x}_{t}$
and
$\widehat{\xi}_{it}=x_{it}-\widehat{\bm{\lambda}}_{i}^{\prime}\widehat{\mathbf{F}}_{t}$
are the PC estimators of the factors and idiosyncratic components,
respectively. Consistency of this estimator is proved by Bai (2003, Theorem
6).

######  Remark 6.

The sandwich form of the asymptotic covariances (17) and (20) comes from the
fact that in the log-likelihood we treat the idiosyncratic components as
serially uncorrelated, while, in fact, they might be autocorrelated. Indeed,
if we assumed $\mathbb{E}[\xi_{it}\xi_{is}]=0$ for all $s,t=1,\ldots,T$ with
$s\neq t$ and all $i=1,\ldots,T$, then
$\bm{\Phi}_{i}=\sigma_{i}^{2}\mathbf{I}_{r}$ (see also Remark 3). It follows
that:
$\bm{\mathcal{I}}_{i}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\bm{\varphi})={\sigma_{i}^{-2}}\mathbf{I}_{r}=-\bm{\mathcal{H}}_{i}(\bm{\mathcal{X}}|\bm{\mathcal{F}};\bm{\varphi})$.
And, by virtue of Theorem 5,
$\bm{\mathcal{I}}_{i}(\bm{\mathcal{X}};\bm{\varphi})$ must coincide
asymptotically with $-\bm{\mathcal{H}}_{i}(\bm{\mathcal{X}};\bm{\varphi})$,
thus, from (17), we have
$\text{A$\\!\mathbb{V}$ar}_{0}(\sqrt{T}\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML,E}})=\\{\bm{\mathcal{I}}_{i}(\bm{\mathcal{X}};\bm{\varphi})\\}^{-1}$.

##  7 Estimation of the factors

Up to this point we did not say anything about estimating the factors. This is
because, first, the factors are typically estimated once we have estimated the
loadings, and, second, there is not a clear definition of what a QML estimator
of the factors is.

If we treat the factors as random variables, then they are not parameters and
they do not have a QML estimator. We know that their optimal (in mean-squared
sense) estimator is their conditional mean given $\bm{\mathcal{X}}$, which,
under Gaussianity can be estimated as the linear projection. If, consistently
with the discussion in Section 4 about QML estimation, we mis-specify the
second order structure of the factors, by replacing $\bm{\Omega}^{F}$ with
$\mathbf{I}_{rT}$, and of the idiosyncratic components, by replacing
$\bm{\Omega}^{\xi}$ with $\bm{\Sigma}^{\xi}$, then such linear projection is
given by:

$\mathbf{F}_{t}^{\text{\tiny
LP}}=\bm{\Lambda}^{\prime}\left(\bm{\Lambda}\bm{\Lambda}^{\prime}+\bm{\Sigma}^{\xi}\right)^{-1}\mathbf{x}_{t}=\left(\bm{\Lambda}^{\prime}(\bm{\Sigma}^{\xi})^{-1}\bm{\Lambda}+\mathbf{I}_{r}\right)^{-1}\bm{\Lambda}^{\prime}(\bm{\Sigma}^{\xi})^{-1}\mathbf{x}_{t},\quad
t=1,\ldots,T,$ (21)

where we used Woodbury formula. Alternatively, it is common to use the OLS or
GLS estimators:

$\mathbf{F}_{t}^{\text{\tiny
OLS}}=\left(\bm{\Lambda}^{\prime}\bm{\Lambda}\right)^{-1}\bm{\Lambda}^{\prime}\mathbf{x}_{t},\;\text{
or }\;\mathbf{F}_{t}^{\text{\tiny
GLS}}=\left(\bm{\Lambda}^{\prime}(\bm{\Sigma}^{\xi})^{-1}\bm{\Lambda}\right)^{-1}\bm{\Lambda}^{\prime}(\bm{\Sigma}^{\xi})^{-1}\mathbf{x}_{t},\quad
t=1,\ldots,T.$ (22)

If instead we treat the factors as constant parameters then, as discussed in
Remark 5, we can see the GLS estimator as the QML estimator maximizing the
joint log-likelihood (15) of factors and data.

The three estimators defined above are unfeasible unless we first compute
estimates of the parameters. By virtue of our result in Theorem 3 there is
asymptotically no difference if we use $\widehat{\bm{\Lambda}}$ or
$\widehat{\bm{\Lambda}}^{\text{\tiny QML}}$, and an estimator of
$\bm{\Sigma}^{\xi}$ can easily be computed either from the residuals of PC
estimation or by QML as suggested by Bai and Li (2012). Once we use these
estimated parameters in (21) and (22), we have the estimators
$\widehat{\mathbf{F}}_{t}^{\text{\tiny LP}}$,
$\widehat{\mathbf{F}}_{t}^{\text{\tiny OLS}}$, and
$\widehat{\mathbf{F}}_{t}^{\text{\tiny GLS}}$. Notice that, as shown in
Section 3, $\widehat{\mathbf{F}}_{t}^{\text{\tiny OLS}}$ is nothing else but
the PC estimator of the factors. The GLS has also been studied by Breitung and
Tenhofen (2011) and Choi (2012) when using different estimators of
$\bm{\Sigma}^{\xi}$.

By construction, the OLS and GLS in (22) are both less efficient than the
linear projection in (21). Moreover, the GLS is always more efficient than the
OLS if we could use an estimator of the full idiosyncratic covariance matrix,
but, since this is in general unfeasible and we typically estiamate only its
diagonal as in (22), then we can just conjecture that the more sparse is the
true covariance matrix the more likely the GLS is to be more efficient.

Finally, these estimators are all $\min(\sqrt{n},T)$-consistent. For the OLS
we refer to Bai (2003, Theorem 1). For the GLS we refer to Bai and Li (2016,
Theorem 2). Moreover, it is straightforward to see that, by Lemmas 1(v) and
2(i), we have $\|\widehat{\mathbf{F}}_{t}^{\text{\tiny
GLS}}-\widehat{\mathbf{F}}_{t}^{\text{\tiny LP}}\|=O_{\mathrm{P}}(n^{-1})$.

######  Remark 7.

If we explicitly model the dynamics of $\mathbf{F}_{t}$ then the expression of
$\widehat{\mathbf{F}}_{t}^{\text{\tiny LP}}$ in (21) is replaced by the Kalman
smoother (Doz et al., 2011). This, in fact, can be shown to be asymptotically
equivalent, as $n\to\infty$, to the GLS estimator in (22) (Bai and Li, 2016;
Ruiz and Poncela, 2022).

######  Remark 8.

In the case of deterministic factors, we could also write the mis-specified
log-likelihood (13) of an exact factor model as:

$\ell_{\text{\tiny
E}}(\bm{\mathcal{X}};\underline{\bm{\varphi}})=-\frac{1}{2}\sum_{i=1}^{n}\log\det(\underline{\bm{F}}\,\underline{\bm{\lambda}}_{i}\underline{\bm{\lambda}}_{i}^{\prime}\underline{\bm{F}}^{\prime}+\underline{\sigma}_{i}^{2}\mathbf{I}_{T})-\frac{1}{2}\sum_{i=1}^{n}\bm{x}_{i}^{\prime}(\underline{\bm{F}}\,\underline{\bm{\lambda}}_{i}\underline{\bm{\lambda}}_{i}^{\prime}\underline{\bm{F}}^{\prime}+\underline{\sigma}_{i}^{2}\mathbf{I}_{T})^{-1}\bm{x}_{i},$
(23)

where $\bm{x}_{i}=(x_{i1}\cdots x_{iT})^{\prime}$, thereby exchanging the role
of $n$ and $T$. Then, we can conjecture that the QML estimator of the factors
maximizing (23) will be asymptotically equivalent, this time as $T\to\infty$,
to their PC estimator. This is approach is also considered in Fortin et al.
(2023). However, since, as noted above, the PC estimator of the factors is
asymptotically equivalent to the OLS, it is not the most efficient estimator
because it does not account for possible cross-sectional heteroskedasticity of
idiosyncratic components.

##  8 Monte Carlo study

Throughout, we let $n\in\\{20,50,100,200\\}$, $T=100$, and $r=2$, and, for all
$i=1,\ldots,n$ and $t=1,\ldots,T$, we simulate the data according to

$\displaystyle x_{it}$
$\displaystyle={\bm{\ell}}_{i}^{\prime}{\bm{f}}_{t}+\phi_{i}\xi_{it},\quad{\bm{f}}_{t}={\bm{A}}{\bm{f}}_{t-1}+\mathbf{u}_{t},\quad\xi_{it}=\alpha_{i}\xi_{it-1}+e_{it},$

where $\bm{\ell}_{i}$ and $\bm{f}_{t}$ are $r$-dimensional vectors.
Specifically,  • $\bm{\ell}_{i}$has entries
${\ell}_{ij}{\sim}iid\,\mathcal{N}(1,1)$, $i=1,\ldots,N$, $j=1,\ldots,r$;  •
${\bm{A}}=0.9\check{\bm{A}}\|{\bm{A}}\|^{-1}$, where $\check{\bm{A}}$ is
$r\times r$ with entries $\check{a}_{jj}{\sim}iid\,\mathcal{U}[0.5,0.8]$ for
all $j$, and $\check{a}_{jk}{\sim}iid\,\mathcal{U}[0,0.3]$ if $j\neq k$;  •
$u_{jt}{\sim}(0,1)$, $j=1,\ldots,r$, $t=1,\ldots,t$, with either a Gaussian or
an Asymmetric Laplace distribution and with
$\mathbb{C}\mathrm{ov}(u_{it},u_{jt})=0$ if $i\neq j$, and
$\mathbb{C}\mathrm{ov}(u_{it},u_{jt-k})=0$ for all $i,j$ and if $k\neq 0$;  •
$e_{it}{\sim}(0,\sigma_{ei}^{2})$, with either a Gaussian or an Asymmetric
Laplace distribution and with $\sigma_{ei}^{2}\sim\mathcal{U}[0.5,1.5]$ for
all $i$,
$\mathbb{C}\mathrm{ov}(e_{it},e_{jt})=\tau^{|i-j|}\mathbb{I}(|i-j|\leq 10)$
with $\tau\in\\{0,0.5\\}$ if $i\neq j$, and
$\mathbb{C}\mathrm{ov}(e_{it},e_{jt-k})=0$ for all $i,j$ and if $k\neq 0$;  •
$\delta_{i}{\sim}iid\,\mathcal{U}(0,\delta)$, and $\delta\in\\{0,0.5\\}$;  •
$\phi_{i}=\\{{\theta_{i}(\sum_{t=1}^{T}\chi_{it}^{2})/(\sum_{t=1}^{T}\xi_{it}^{2})}\\}^{1/2}$,
and $\theta_{i}{\sim}iid\,\mathcal{U}(0.25,0.5)$.  In the case of the
Asymmetric Laplace distribution, all the innovations have location 0,
asymmetry index $\kappa$, with $\kappa\sim\mathcal{U}(.9,1.1)$ and scale index
$\lambda=\sqrt{{1+\kappa^{4}}/{\kappa^{2}}}$, so that the variance is 1. The
parameters $\tau$ and $\delta$ control the degrees of cross-sectional and
serial idiosyncratic correlation in the idiosyncratic components. The the
noise-to-signal ratio for series $i$ is given by $\theta_{i}$.

Finally, in order to satisfy Assumptions 5(a) and 5(b) we proceed as follows.
Given the common components is generated as
$\chi_{it}={\bm{\ell}}_{i}^{\prime}{\bm{f}}_{t}$, let
$\bm{\chi}_{t}=(\chi_{1t}\cdots\chi_{nt})^{\prime}$ and compute
$\check{\bm{\Gamma}}^{\chi}=T^{-1}\sum_{t=1}^{T}\bm{\chi}_{t}\bm{\chi}_{t}^{\prime}$.
Collect its $r$ non-zero eigenvalues into the $r\times r$ diagonal matrix
$\check{\mathbf{M}}^{\chi}$ and the corresponding normalized eigenvectors as
the columns of the $n\times r$ matrix $\check{\mathbf{V}}^{\chi}$, with sign
fixed such that it has non-negative entries in the first row. The loadings are
then simulated as
$\bm{\Lambda}=\check{\mathbf{V}}^{\chi}(\check{\mathbf{M}}^{\chi})^{1/2}$ and
the factors as
$\mathbf{F}_{t}=(\check{\mathbf{M}}^{\chi})^{-1/2}\check{\mathbf{V}}^{\chi\prime}\bm{\chi}_{t}$.

We simulate the model described above $B=500$ times, and at each replication
we estimate the loadings via PC and QML, where the latter is defined as the
maximizer of the log-likelihood (13) and it is computed numerically in the way
proposed by Bai and Li (2012, 2016). In this way, we obtain two estimates of
the loadings matrix $\widehat{\bm{\Lambda}}^{(b)}$ and
$\widehat{\bm{\Lambda}}^{\text{\tiny QML,E}(b)}$, respectively,
$b=1,\ldots,B$. We also compute the unfeasible OLS estimator
${\bm{\Lambda}}^{\text{\tiny OLS}(b)}$ by regressing $\mathbf{x}_{t}$ onto the
simulated factors.

Table 1: MSEs Gaussian innovations
---
$n$ | $T$ | $\tau$ | $\delta$ | $\text{MSE}_{1}^{\text{\tiny OLS}}$ | $\text{MSE}_{1}^{\text{\tiny PC}}$ | $\text{MSE}_{1}^{\text{\tiny QML}}$ | $\text{MSE}_{2}^{\text{\tiny OLS}}$ | $\text{MSE}_{2}^{\text{\tiny PC}}$ | $\text{MSE}_{2}^{\text{\tiny QML}}$
20 | 100 | 0 | 0 | 0.0103 | 0.0123 | 0.0116 | 0.0099 | 0.0153 | 0.0118
|  |  |  | (0.0032) | (0.0060) | (0.0053) | (0.0034) | (0.0069) | (0.0053)
50 | 100 | 0 | 0 | 0.0102 | 0.0109 | 0.0108 | 0.0102 | 0.0116 | 0.0107
|  |  |  | (0.0021) | (0.0023) | (0.0023) | (0.0022) | (0.0025) | (0.0023)
100 | 100 | 0 | 0 | 0.0100 | 0.0103 | 0.0103 | 0.0100 | 0.0105 | 0.0104
|  |  |  | (0.0015) | (0.0016) | (0.0015) | (0.0015) | (0.0016) | (0.0016)
200 | 100 | 0 | 0 | 0.0101 | 0.0102 | 0.0102 | 0.0101 | 0.0104 | 0.0103
|  |  |  | (0.0011) | (0.0011) | (0.0011) | (0.0011) | (0.0011) | (0.0011)
20 | 100 | 0.5 | 0.5 | 0.0180 | 0.0239 | 0.0230 | 0.0134 | 0.0270 | 0.0276
|  |  |  | (0.0081) | (0.0115) | (0.0110) | (0.0055) | (0.0188) | (0.0217)
50 | 100 | 0.5 | 0.5 | 0.0183 | 0.0201 | 0.0199 | 0.0135 | 0.0166 | 0.0160
|  |  |  | (0.0054) | (0.0060) | (0.0060) | (0.0037) | (0.0050) | (0.0047)
100 | 100 | 0.5 | 0.5 | 0.0182 | 0.0190 | 0.0189 | 0.0134 | 0.0146 | 0.0145
|  |  |  | (0.0038) | (0.0041) | (0.0041) | (0.0027) | (0.0031) | (0.0030)
200 | 100 | 0.5 | 0.5 | 0.0184 | 0.0187 | 0.0187 | 0.0135 | 0.0141 | 0.0141
|  |  |  | (0.0027) | (0.0028) | (0.0028) | (0.0021) | (0.0022) | (0.0022)
Asymmetric Laplace innovations
$n$ | $T$ | $\tau$ | $\delta$ | $\text{MSE}_{1}^{\text{\tiny OLS}}$ | $\text{MSE}_{1}^{\text{\tiny PC}}$ | $\text{MSE}_{1}^{\text{\tiny QML}}$ | $\text{MSE}_{2}^{\text{\tiny OLS}}$ | $\text{MSE}_{2}^{\text{\tiny PC}}$ | $\text{MSE}_{2}^{\text{\tiny QML}}$
20 | 100 | 0 | 0 | 0.0401 | 0.0499 | 0.0440 | 0.0408 | 0.0654 | 0.0389
|  |  |  | (0.0136) | (0.0211) | (0.0180) | (0.0144) | (0.0254) | (0.0145)
50 | 100 | 0 | 0 | 0.0410 | 0.0442 | 0.0423 | 0.0407 | 0.0480 | 0.0390
|  |  |  | (0.0091) | (0.0102) | (0.0097) | (0.0090) | (0.0112) | (0.0087)
100 | 100 | 0 | 0 | 0.0412 | 0.0429 | 0.0418 | 0.0410 | 0.0439 | 0.0391
|  |  |  | (0.0068) | (0.0073) | (0.0071) | (0.0068) | (0.0077) | (0.0066)
200 | 100 | 0 | 0 | 0.0409 | 0.0420 | 0.0435 | 0.0411 | 0.0421 | 0.0421
|  |  |  | (0.0046) | (0.0050) | (0.0078) | (0.0051) | (0.0054) | (0.0075)
20 | 100 | 0.5 | 0.5 | 0.0755 | 0.1169 | 0.1016 | 0.0572 | 0.1257 | 0.1050
|  |  |  | (0.0353) | (0.1469) | (0.1405) | (0.0270) | (0.1577) | (0.1378)
50 | 100 | 0.5 | 0.5 | 0.0745 | 0.0864 | 0.0797 | 0.0549 | 0.0726 | 0.0587
|  |  |  | (0.0252) | (0.0292) | (0.0271) | (0.0179) | (0.0264) | (0.0191)
100 | 100 | 0.5 | 0.5 | 0.0742 | 0.0814 | 0.0765 | 0.0566 | 0.0640 | 0.0555
|  |  |  | (0.0183) | (0.0216) | (0.0194) | (0.0136) | (0.0162) | (0.0128)
200 | 100 | 0.5 | 0.5 | 0.0760 | 0.0800 | 0.0774 | 0.0548 | 0.0594 | 0.0542
|  |  |  | (0.0160) | (0.0178) | (0.0169) | (0.0106) | (0.0128) | (0.0108)
Table 2: Comparison between PC and QML estimators. Gaussian innovations
---
$n$ | $T$ | $\tau$ | $\delta$ | $D_{1}$ | $D_{2}$ | ${\text{MSE}_{1}^{\text{\tiny REL}}}$ | ${\text{MSE}_{2}^{\text{\tiny REL}}}$
20 | 100 | 0 | 0 |  5.06$\times 10^{-4}$ |  2.79$\times 10^{-3}$ | 1.06 | 1.30
|  |  |  | (1.72$\times 10^{-4}$) | (2.61$\times 10^{-3}$) |  |
50 | 100 | 0 | 0 |  6.07$\times 10^{-5}$ |  4.50$\times 10^{-4}$ | 1.01 | 1.08
|  |  |  | (1.59$\times 10^{-5}$) | (2.01$\times 10^{-4}$) |  |
100 | 100 | 0 | 0 |  1.23$\times 10^{-5}$ |  6.28$\times 10^{-5}$ | 1.00 | 1.01
|  |  |  | (3.54$\times 10^{-6}$) | (3.00$\times 10^{-5}$) |  |
200 | 100 | 0 | 0 |  4.54$\times 10^{-6}$ |  2.05$\times 10^{-5}$ | 1.00 | 1.00
|  |  |  | (1.02$\times 10^{-6}$) | (6.96$\times 10^{-6}$) |  |
20 | 100 | 0.5 | 0.5 |  5.72$\times 10^{-4}$ |  5.86$\times 10^{-3}$ | 1.04 | 0.98
|  |  |  | (3.27$\times 10^{-4}$) |  (1.15$\times 10^{-2}$) |  |
50 | 100 | 0.5 | 0.5 |  7.55$\times 10^{-5}$ |  5.63$\times 10^{-4}$ | 1.01 | 1.04
|  |  |  | (3.65$\times 10^{-5}$) |  (4.88$\times 10^{-4}$) |  |
100 | 100 | 0.5 | 0.5 |  1.77$\times 10^{-5}$ |  7.93$\times 10^{-5}$ | 1.00 | 1.01
|  |  |  | (7.26$\times 10^{-6}$) |  (5.36$\times 10^{-5}$) |  |
200 | 100 | 0.5 | 0.5 |  6.06$\times 10^{-6}$ |  2.34$\times 10^{-5}$ | 1.00 | 1.00
|  |  |  | (1.95$\times 10^{-6}$) |  (1.02$\times 10^{-5}$) |  |
Asymmetric Laplace innovations
$n$ | $T$ | $\tau$ | $\delta$ | $D_{1}$ | $D_{2}$ | ${\text{MSE}_{1}^{\text{\tiny REL}}}$ | ${\text{MSE}_{2}^{\text{\tiny REL}}}$
20 | 100 | 0 | 0 |  4.35$\times 10^{-3}$ |  2.09$\times 10^{-2}$ | 1.13 | 1.68
|  |  |  | (1.67$\times 10^{-3}$) |  (1.68$\times 10^{-2}$) |  |
50 | 100 | 0 | 0 |  1.13$\times 10^{-3}$ |  5.67$\times 10^{-3}$ | 1.05 | 1.23
|  |  |  | (3.33$\times 10^{-4}$) |  (3.19$\times 10^{-3}$) |  |
100 | 100 | 0 | 0 |  6.23$\times 10^{-4}$ |  2.64$\times 10^{-3}$ | 1.03 | 1.12
|  |  |  | (1.48$\times 10^{-4}$) |  (1.28$\times 10^{-3}$) |  |
200 | 100 | 0 | 0 |  2.20$\times 10^{-3}$ |  2.07$\times 10^{-3}$ | 0.97 | 1.00
|  |  |  | (3.81$\times 10^{-3}$) |  (2.41$\times 10^{-3}$) |  |
20 | 100 | 0.5 | 0.5 |  9.93$\times 10^{-3}$ |  4.12$\times 10^{-2}$ | 1.15 | 1.20
|  |  |  | (6.77$\times 10^{-3}$) |  (5.68$\times 10^{-2}$) |  |
50 | 100 | 0.5 | 0.5 |  4.68$\times 10^{-3}$ |  9.16$\times 10^{-3}$ | 1.08 | 1.24
|  |  |  | (2.07$\times 10^{-3}$) |  (9.07$\times 10^{-3}$) |  |
100 | 100 | 0.5 | 0.5 |  3.13$\times 10^{-3}$ |  4.56$\times 10^{-3}$ | 1.06 | 1.15
|  |  |  | (1.05$\times 10^{-3}$) |  (2.92$\times 10^{-3}$) |  |
200 | 100 | 0.5 | 0.5 |  2.63$\times 10^{-3}$ |  3.41$\times 10^{-3}$ | 1.03 | 1.10
|  |  |  | (1.78$\times 10^{-3}$) |  (2.73$\times 10^{-3}$) |  |

In Table 1 we report the Mean-Squared-Error (MSE) for each column,
$j=1,\ldots,r$, of the considered estimators, averaged over the $B$
replications (with standard deviations in parenthesis):

$\displaystyle\text{MSE}_{j}^{\text{\tiny
OLS}}=\frac{1}{B}\sum_{b=1}^{B}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left({\bm{\Lambda}}^{\text{\tiny
OLS}(b)}_{ij}-\bm{\Lambda}_{ij}\right)^{2}\right\\},$
$\displaystyle\text{MSE}_{j}^{\text{\tiny
PC}}=\frac{1}{B}\sum_{b=1}^{B}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{\bm{\Lambda}}^{(b)}_{ij}-\bm{\Lambda}_{ij}\right)^{2}\right\\},\quad\text{MSE}_{j}^{\text{\tiny
QML}}=\frac{1}{B}\sum_{b=1}^{B}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}(b)}_{ij}-\bm{\Lambda}_{ij}\right)^{2}\right\\}.$

Results show that QML and PC estimator have similar MSEs and both improve as
$n$ increases to the point that when $n=200$ their MSEs are comparable with
the one of the unfeasible OLS, i.e., the estimators behave as if the factors
were observed. In most cases the QML estimator has a smaller MSE even when the
true distribution is not Gaussian.

Finally, in Table 2 for each column of the loadings, $j=1,\ldots,r$, we report
the distance between the QML and PC estimators measured as (with standard
deviations in parenthesis):

$\displaystyle\text{D}_{j}=\frac{1}{B}\sum_{b=1}^{B}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}(b)}_{ij}-\widehat{\bm{\Lambda}}_{ij}^{(b)}\right)^{2}\right\\}.$

And we also report the relative MSE of the PC estimator with respect to the
MSE of the QML estimator: $\text{MSE}_{j}^{\text{\tiny
REL}}={\text{MSE}_{j}^{\text{\tiny PC}}}/{\text{MSE}_{j}^{\text{\tiny QML}}}$.
Results clearly show that as $n$ grows the PC and QML estimators become almost
indistinguishable.

##  9 Concluding remarks

To compute in practice the QML estimator of the loadings there are at least
two main issues. First, in finite samples the QML estimator of the loadings
has no closed form and depends also on the estimator of the idiosyncratic
covariance for which no closed form exists either. This is true even if we use
the the log-likelihood (13) of an exact factor model. Second, the convergence
properties of the various available EM algorithms used to compute the QML
estimator (Rubin and Thayer, 1982; Bai and Li, 2012, 2016; Zadrozny, 2023)
have never been fully investigated. On the one hand, it is easy to prove, that
at each iteration of an EM algorithm the log-likelihood evaluated in the
parameters estimated at that iteration is larger than at the previous
iteration (Wu, 1983), but, on the other hand, no formal proof exists of
convergence of those algorithms to a global maximum of the likelihood, at
least to the best of our knowledge.

The results of this paper offer a possible solution by showing that, if we are
just interested in the factor loadings and we do not need to estimate the
idiosyncratic variances, then we can simply use the PC estimator of the
loadings and of its asymptotic covariance matrix to approximate the
corresponding QML estimator and its asymptotic covariance matrix. Once this is
done, the factors can be estimated via OLS as in PC analysis.

As a consequence, we might think that there is no apparent advantage in
directly computing the QML estimator of the loadings and of the idiosyncratic
variances.

Nevertheless, QML estimation has at least three advantages. First, it allows
us to easily impose restrictions on the parameters of the model. Second,
having also the QML estimator of the idiosyncratic variances allows us to
compute estimators of the factors as the GLS, which are possibly more
efficient. Third, QML estimation, as presented in this paper, is a first step
towards estimating a model where we explicitly model the dynamics of the
factors, something we cannot do with PC analysis. This last point, already
briefly discussed in Remarks 4 and 7, is the subject of our ongoing research
(Barigozzi and Luciani, 2024).

## References

  * Aït-Sahalia and Xiu (2017) Aït-Sahalia, Y. and D. Xiu (2017).  Using principal component analysis to estimate a high dimensional factor model with high-frequency data.  Journal of Econometrics 201, 384–399.
  * Amemiya et al. (1987) Amemiya, Y., W. A. Fuller, and S. G. Pantula (1987).  The asymptotic distributions of some estimators for a factor analysis model.  Journal of Multivariate Analysis 22, 51–64.
  * Anderson and Amemiya (1988) Anderson, T. W. and Y. Amemiya (1988).  The asymptotic normal distribution of estimators in factor analysis under general conditions.  The Annals of Statistics 16, 759–771.
  * Anderson and Rubin (1956) Anderson, T. W. and H. Rubin (1956).  Statistical inference in factor analysis.  In Proceedings of the third Berkeley symposium on mathematical statistics and probability, Volume 5, pp. 111–150.
  * Bai (2003) Bai, J. (2003).  Inferential theory for factor models of large dimensions.  Econometrica 71, 135–171.
  * Bai and Li (2012) Bai, J. and K. Li (2012).  Statistical analysis of factor models of high dimension.  The Annals of Statistics 40, 436–465.
  * Bai and Li (2016) Bai, J. and K. Li (2016).  Maximum likelihood estimation and inference for approximate factor models of high dimension.  The Review of Economics and Statistics 98, 298–309.
  * Bai and Liao (2016) Bai, J. and Y. Liao (2016).  Efficient estimation of approximate factor models via penalized maximum likelihood.  Journal of Econometrics 191, 1–18.
  * Bai and Ng (2002) Bai, J. and S. Ng (2002).  Determining the number of factors in approximate factor models.  Econometrica 70, 191–221.
  * Bai and Ng (2013) Bai, J. and S. Ng (2013).  Principal components estimation and identification of static factors.  Journal of Econometrics 176, 18–29.
  * Barigozzi (2024) Barigozzi, M. (2024).  Quasi maximum likelihood estimation of high-dimensional factor models.  Oxford Research Encyclopedia of Economics and Finance.  forthcoming.
  * Barigozzi and Luciani (2024) Barigozzi, M. and M. Luciani (2024).  Quasi maximum likelihood estimation and inference of large approximate dynamic factor models via the EM algorithm.  Technical Report arXiv:1910.03821.
  * Bernanke et al. (2005) Bernanke, B. S., J. Boivin, and P. S. Eliasz (2005).  Measuring the effects of monetary policy: A Factor-Augmented Vector Autoregressive (FAVAR) approach.  The Quarterly Journal of Economics 120, 387–422.
  * Breitung and Tenhofen (2011) Breitung, J. and J. Tenhofen (2011).  GLS estimation of dynamic factor models.  Journal of the American Statistical Association 106, 1150–1166.
  * Chamberlain and Rothschild (1983) Chamberlain, G. and M. Rothschild (1983).  Arbitrage, factor structure, and mean-variance analysis on large asset markets.  Econometrica 51, 1281–1304.
  * Choi (2012) Choi, I. (2012).  Efficient estimation of factor models.  Econometric Theory 28, 274–308.
  * Connor et al. (2006) Connor, G., R. A. Korajczyk, and O. Linton (2006).  The common and specific components of dynamic volatility.  Journal of Econometrics 132, 231–255.
  * Coroneo et al. (2016) Coroneo, L., D. Giannone, and M. Modugno (2016).  Unspanned macroeconomic factors in the yield curve.  Journal of Business and Economic Statistics 34, 472–485.
  * De Mol et al. (2008) De Mol, C., D. Giannone, and L. Reichlin (2008).  Forecasting using a large number of predictors: Is bayesian shrinkage a valid alternative to principal components?  Journal of Econometrics 146, 318–328.
  * Delle Chiaie et al. (2021) Delle Chiaie, S., L. Ferrara, and D. Giannone (2021).  Common factors of commodity prices.  Journal of Applied Econometrics.  forthcoming.
  * Doz et al. (2011) Doz, C., D. Giannone, and L. Reichlin (2011).  A two-step estimator for large approximate dynamic factor models based on Kalman filtering.  Journal of Econometrics 164, 188–205.
  * Doz et al. (2012) Doz, C., D. Giannone, and L. Reichlin (2012).  A quasi maximum likelihood approach for large approximate dynamic factor models.  The Review of Economics and Statistics 94(4), 1014–1024.
  * Fan et al. (2013) Fan, J., Y. Liao, and M. Mincheva (2013).  Large covariance estimation by thresholding principal orthogonal complements.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75, 603–680.
  * Forni et al. (2000) Forni, M., M. Hallin, M. Lippi, and L. Reichlin (2000).  The Generalized Dynamic Factor Model: Identification and estimation.  The Review of Economics and Statistics 82, 540–554.
  * Forni et al. (2005) Forni, M., M. Hallin, M. Lippi, and L. Reichlin (2005).  The Generalized Dynamic Factor Model: One sided estimation and forecasting.  Journal of the American Statistical Association 100, 830–840.
  * Fortin et al. (2023) Fortin, A.-P., P. Gagliardini, and O. Scaillet (2023).  Latent factor analysis in short panels.  Technical Report arXiv:2306.14004.
  * Giannone et al. (2008) Giannone, D., L. Reichlin, and D. Small (2008).  Nowcasting: The real-time informational content of macroeconomic data.  Journal of Monetary Economics 55, 665–676.
  * Hannan (1970) Hannan, E. J. (1970).  Multiple time series.  John Wiley & Sons.
  * Jolliffe (2002) Jolliffe, I. T. (2002).  Principal component analysis.  Springer.
  * Kim and Fan (2019) Kim, D. and J. Fan (2019).  Factor GARCH-Itô models for high-frequency data with application to large volatility matrix prediction.  Journal of Econometrics 208, 395–417.
  * Lam and Yao (2012) Lam, C. and Q. Yao (2012).  Factor modeling for high-dimensional time series: inference for the number of factors.  The Annals of Statistics 40, 694–726.
  * Lawley and Maxwell (1971) Lawley, D. N. and A. E. Maxwell (1971).  Factor Analysis as a Statistical Method.  Butterworths, London.
  * Mardia et al. (1979) Mardia, K., J. Kent, and J. Bibby (1979).  Multivariate analysis.
  * Merikoski and Kumar (2004) Merikoski, J. K. and R. Kumar (2004).  Inequalities for spreads of matrix sums and products.  Applied Mathematics E-Notes 4, 150–159.
  * Onatski (2010) Onatski, A. (2010).  Determining the number of factors from empirical distribution of eigenvalues.  The Review of Economics and Statistics 92, 1004–1016.
  * Rubin and Thayer (1982) Rubin, D. B. and D. T. Thayer (1982).  EM algorithms for ML factor analysis.  Psychometrika 47, 69–76.
  * Ruiz and Poncela (2022) Ruiz, E. and P. Poncela (2022).  Factor extraction in Dynamic Factor Models: Using Kalman Filter and Principal Components in practice.  Foundations and Trends in Econometrics.  forthcoming.
  * Stock and Watson (2002a) Stock, J. H. and M. W. Watson (2002a).  Forecasting using principal components from a large number of predictors.  Journal of the American Statistical Association 97, 1167–1179.
  * Stock and Watson (2002b) Stock, J. H. and M. W. Watson (2002b).  Macroeconomic forecasting using diffusion indexes.  Journal of Business and Economic Statistics 20, 147–162.
  * Tipping and Bishop (1999) Tipping, M. E. and C. M. Bishop (1999).  Probabilistic principal component analysis.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 61, 611–622.
  * Trapani (2018) Trapani, L. (2018).  A randomized sequential procedure to determine the number of factors.  Journal of the American Statistical Association 113, 1341–1349.
  * Wu (1983) Wu, J. C. F. (1983).  On the convergence properties of the EM algorithm.  The Annals of Statistics 11, 95–103.
  * Yu et al. (2015) Yu, Y., T. Wang, and R. J. Samworth (2015).  A useful variant of the Davis-Kahan theorem for statisticians.  Biometrika 102, 315–323.
  * Zadrozny (2023) Zadrozny, P. A. (2023).  Gaussian maximum likelihood estimation of a static-form factor model using the Expectation-Maximization algorithm.  mimeo, Bureau of Labor Statistics.

##  Appendix A Proof of main results

###  A.1 Proof of Theorem 1

From (7) and (8), and since by (9) we have
$\widehat{\bm{\Lambda}}^{\prime}\widehat{\bm{\Lambda}}=\widehat{\mathbf{M}}^{x}$
and
$\widehat{\mathbf{V}}^{x}=\widehat{\bm{\Lambda}}(\widehat{\mathbf{M}}^{x})^{-1/2}$
which is well defined because of Lemma 8(iv), we get

$\frac{\bm{X}^{\prime}\bm{X}}{nT}\widehat{\bm{\Lambda}}=\widehat{\bm{\Lambda}}\frac{\widehat{\mathbf{M}}^{x}}{n}.$
(A.1)

Then, substituting
$\bm{X}^{\prime}\bm{X}=(\bm{\Lambda}\bm{F}^{\prime}+\bm{\Xi}^{\prime})^{\prime}(\bm{F}\bm{\Lambda}^{\prime}+\bm{\Xi})$
into (A.1)

$\displaystyle\frac{\bm{\Lambda}\bm{F}^{\prime}\bm{F}\bm{\Lambda}^{\prime}\widehat{\bm{\Lambda}}}{nT}+\frac{\bm{\Lambda}\bm{F}^{\prime}\bm{\Xi}\widehat{\bm{\Lambda}}}{nT}+\frac{\bm{\Xi}^{\prime}\bm{F}\bm{\Lambda}^{\prime}\widehat{\bm{\Lambda}}}{nT}+\frac{\bm{\Xi}^{\prime}\bm{\Xi}\widehat{\bm{\Lambda}}}{nT}=\widehat{\bm{\Lambda}}\frac{\widehat{\mathbf{M}}^{x}}{n}.$
(A.2)

Define

$\widehat{\mathbf{H}}=\left(\frac{\bm{F}^{\prime}\bm{F}}{T}\right)\left(\frac{\bm{\Lambda}^{\prime}\widehat{\bm{\Lambda}}}{n}\right)\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}=\left(\frac{\bm{\Lambda}^{\prime}\widehat{\bm{\Lambda}}}{n}\right)\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1},$
(A.3)

by Assumption 5(b). Notice that, as $n,T\to\infty$, $\widehat{\mathbf{H}}$ is
well defined because $\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}$ is
well defined because of Lemma 8(iv). From (A.2) and (A.3)

$\displaystyle\widehat{\bm{\Lambda}}$
$\displaystyle-\bm{\Lambda}\widehat{\mathbf{H}}=\left(\frac{\bm{\Lambda}\bm{F}^{\prime}\bm{\Xi}\widehat{\bm{\Lambda}}}{nT}+\frac{\bm{\Xi}^{\prime}\bm{F}\bm{\Lambda}^{\prime}\widehat{\bm{\Lambda}}}{nT}+\frac{\bm{\Xi}^{\prime}\bm{\Xi}\widehat{\bm{\Lambda}}}{nT}\right)\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}$
$\displaystyle=$
$\displaystyle\,\left(\frac{\bm{\Lambda}\bm{F}^{\prime}\bm{\Xi}{\bm{\Lambda}}}{nT}+\frac{\bm{\Xi}^{\prime}\bm{F}\bm{\Lambda}^{\prime}{\bm{\Lambda}}}{nT}+\frac{\bm{\Xi}^{\prime}\bm{\Xi}{\bm{\Lambda}}}{nT}\right)\bm{\mathcal{H}}\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}\\!\\!\\!\\!\\!+\left(\frac{\bm{\Lambda}\bm{F}^{\prime}\bm{\Xi}}{nT}+\frac{\bm{\Xi}^{\prime}\bm{F}\bm{\Lambda}^{\prime}}{nT}+\frac{\bm{\Xi}^{\prime}\bm{\Xi}}{nT}\right)(\widehat{\bm{\Lambda}}-\bm{\Lambda}\bm{\mathcal{H}})\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}.$
(A.4)

Taking the $i$th row of (A.4)

$\displaystyle\widehat{\bm{\lambda}}_{i}^{\prime}-{\bm{\lambda}}_{i}^{\prime}\widehat{\mathbf{H}}=$
$\displaystyle\,\left(\underbrace{\frac{1}{nT}{\bm{\lambda}}_{i}^{\prime}\sum_{t=1}^{T}\sum_{j=1}^{n}\mathbf{F}_{t}\xi_{jt}{\bm{\lambda}}_{j}^{\prime}}_{\text{(1.a)}}+\underbrace{\frac{1}{nT}\sum_{t=1}^{T}\xi_{it}\mathbf{F}_{t}^{\prime}\sum_{j=1}^{n}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\prime}}_{\text{(1.b)}}+\underbrace{\frac{1}{nT}\sum_{t=1}^{T}\sum_{j=1}^{n}\xi_{it}\xi_{jt}\bm{\lambda}_{j}^{\prime}}_{\text{(1.c)}}\right)\bm{\mathcal{H}}\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}$
$\displaystyle+\left(\underbrace{\frac{1}{nT}{\bm{\lambda}}_{i}^{\prime}\sum_{t=1}^{T}\sum_{j=1}^{n}\mathbf{F}_{t}\xi_{jt}(\widehat{\bm{\lambda}}_{j}^{\prime}-\bm{\lambda}_{j}^{\prime}\bm{\mathcal{H}})}_{\text{(1.d)}}+\underbrace{\frac{1}{nT}\sum_{t=1}^{T}\xi_{it}\mathbf{F}_{t}^{\prime}\sum_{j=1}^{n}\bm{\lambda}_{j}(\widehat{\bm{\lambda}}_{j}^{\prime}-\bm{\lambda}_{j}^{\prime}\bm{\mathcal{H}})}_{\text{(1.e)}}\right.$
$\displaystyle\;\;\left.+\underbrace{\frac{1}{nT}\sum_{t=1}^{T}\sum_{j=1}^{n}\xi_{it}\xi_{jt}(\widehat{\bm{\lambda}}_{j}^{\prime}-\bm{\lambda}_{j}^{\prime}\bm{\mathcal{H}})}_{\text{(1.f)}}\right)\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}.$
(A.5)

From Proposition B.3 we see that, under Assumptions 1 through 4, the terms
(1.a), (1.c), (1.d), (1.e), and (1.f) are all
$o_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right)$. In particular, from (A.5),
for any $i=1,\ldots,n$, we get

$\displaystyle\widehat{\bm{\lambda}}_{i}-\widehat{\mathbf{H}}^{\prime}{\bm{\lambda}}_{i}=\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}\bm{\mathcal{H}}^{\prime}\left(\frac{\bm{\Lambda}^{\prime}\bm{\Lambda}}{n}\right)\left(\frac{1}{T}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}\right)+O_{\mathrm{P}}\left(\max\left(\frac{1}{n},\frac{1}{\sqrt{nT}}\right)\right).$
(A.6)

Consistency follows immediately since the first term in (A.6) is
$O_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right)$ because of Assumption 7 (see
also (B.42) in the proof of Theorem 5) and since
$\left\|\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}\bm{\mathcal{H}}^{\prime}\right\|=O_{\mathrm{P}}(1)$
because of Lemma 8(iv) and 11.

Then, from (A.6), by Proposition B.4(b):

$\displaystyle\sqrt{T}(\widehat{\bm{\lambda}}_{i}-\widehat{\mathbf{H}}^{\prime}{\bm{\lambda}}_{i})$
$\displaystyle=\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}\bm{\mathcal{H}}^{\prime}\left(\frac{\bm{\Lambda}^{\prime}\bm{\Lambda}}{n}\right)\left(\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}\right)+o_{\mathrm{P}}(1)$
$\displaystyle=\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}\left(\frac{\widehat{\bm{\Lambda}}^{\prime}\bm{\Lambda}}{n}\right)\left(\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}\right)+o_{\mathrm{P}}(1).$
(A.7)

Now, from Proposition B.6 we have that, as $n,T\to\infty$,
$\min(\sqrt{n},\sqrt{T})\|\widehat{\mathbf{H}}-\bm{J}\|=o_{\mathrm{P}}(1)$,
and, by using the definition of $\widehat{\mathbf{H}}$ in (A.3), it follows
that (A.1) is equivalent to (note that $\|\bm{\lambda}_{i}\|=O(1)$ by
Assumption 1(a))

$\displaystyle\sqrt{T}(\widehat{\bm{\lambda}}_{i}-\bm{J}\bm{\lambda}_{i})$
$\displaystyle=\sqrt{T}(\widehat{\bm{\lambda}}_{i}-\widehat{\mathbf{H}}^{\prime}\bm{\lambda}_{i})+\sqrt{T}(\widehat{\mathbf{H}}^{\prime}-\bm{J})\bm{\lambda}_{i}+o_{\mathrm{P}}(1)$
$\displaystyle=\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}\left(\frac{\widehat{\bm{\Lambda}}^{\prime}\bm{\Lambda}}{n}\right)\left(\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}\right)+o_{\mathrm{P}}(1)$
$\displaystyle=\widehat{\mathbf{H}}^{\prime}\left(\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}\right)+o_{\mathrm{P}}(1)$
$\displaystyle=\bm{J}\left(\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}\right)+o_{\mathrm{P}}(1)$
$\displaystyle\to_{d}\mathcal{N}\left(\mathbf{0}_{r},\bm{\Phi}_{i}\right),$
(A.8)

where we used Slutsky’s theorem and Assumption 7. Notice that, $\bm{J}$ plays
no role in the covariance since it is diagonal and
$\bm{J}^{2}=\mathbf{I}_{r}$. Because of Assumption 6 the sign indeterminacy on
the left-hand-side can be easily fixed so that $\bm{J}=\mathbf{I}_{r}$. By
substituting $\mathbf{I}_{r}$ in place of $\widehat{\mathbf{H}}$ in (A.6) and
using Assumption 7 it follows also that:

$\left\|\widehat{\bm{\lambda}}_{i}-{\bm{\lambda}}_{i}\right\|=O_{\mathrm{P}}\left(\max\left(\frac{1}{n},\frac{1}{\sqrt{T}}\right)\right).$

This completes the proof of part (a).

For part (b), from (A.4)

$\displaystyle\left\|\frac{\widehat{\bm{\Lambda}}-\bm{\Lambda}\widehat{\mathbf{H}}}{\sqrt{n}}\right\|\leq$
$\displaystyle\,\left(\left\|\frac{\bm{F}^{\prime}\bm{\Xi}}{\sqrt{n}T}\right\|\,\left\|\frac{\bm{\Lambda}}{\sqrt{n}}\right\|^{2}+\left\|\frac{\bm{\Xi}^{\prime}\bm{F}}{\sqrt{n}T}\right\|\,\left\|\frac{\bm{\Lambda}^{\prime}{\bm{\Lambda}}}{n}\right\|+\left\|\frac{\bm{\Xi}^{\prime}\bm{\Xi}\bm{\Lambda}}{n^{3/2}T}\right\|\right)\left\|\bm{\mathcal{H}}\right\|\,\left\|\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}\right\|$
$\displaystyle+\left(\left\|\frac{\bm{F}^{\prime}\bm{\Xi}}{\sqrt{n}T}\right\|\,\left\|\frac{\bm{\Lambda}}{\sqrt{n}}\right\|+\left\|\frac{\bm{\Xi}^{\prime}\bm{F}}{\sqrt{n}T}\right\|\,\left\|\frac{\bm{\Lambda}}{\sqrt{n}}\right\|+\left\|\frac{\bm{\Xi}^{\prime}\bm{\Xi}}{nT}\right\|\right)\left\|\frac{\widehat{\bm{\Lambda}}-\bm{\Lambda}\bm{\mathcal{H}}}{\sqrt{n}}\right\|\left\|\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{-1}\right\|,$

and the proof of part (b) follows from Proposition B.2 and Lemma 2(i), 3,
4(i), 4(ii), 4(iii), 8(iv), 11. This proves part (b) and completes the proof.
$\Box$

###  A.2 Proof of Corollary 1

From (A.6) and the second last line of (A.1) in the proof of Theorem 1 and by
imposing the identification constraint of orthonormal factors in Assumption
5(b) and $\bm{J}=\mathbf{I}_{r}$ by Assumption 6, we have

$\displaystyle\widehat{\bm{\lambda}}_{i}-\bm{\lambda}_{i}$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}+O_{\mathrm{P}}\left(\max\left(\frac{1}{n},\frac{1}{\sqrt{nT}}\right)\right),$
(A.9)

where the rate of the last term comes from (A.6) in the proof of Theorem 1(a).
By definition of OLS, and again imposing Assumption 5(b), we have:

$\displaystyle{\bm{\lambda}}_{i}^{\text{\tiny
OLS}}-\bm{\lambda}_{i}=\frac{1}{T}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}.$
(A.10)

By comparing (A.9) and (A.10) we complete the proof of part (a). Part (b)
follows also directly from Theorem 1(b). This completes the proof. $\Box$

###  A.3 Proof of Theorem 3

In principle, we could try to replicate the proofs by Bai and Li (2016) under
our identifying constraints and using only our assumptions. However, there is
a much simpler and intuitive way to proceed.

Consider the log-likelihood (12). The parameters to be estimated are given by
${\bm{\varphi}}=(\mathrm{vec}({\bm{\Lambda}})^{\prime},\mathrm{vech}({\bm{\Gamma}}^{\xi})^{\prime})^{\prime}$.
Let also $\widehat{\bm{\varphi}}^{\text{\tiny
QML}}=(\mathrm{vec}(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}})^{\prime},\mathrm{vech}(\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}})^{\prime})^{\prime}$ denote the maximizer of (12) and
$\underline{\bm{\varphi}}=(\mathrm{vec}(\underline{\bm{\Lambda}})^{\prime},\mathrm{vech}(\underline{\bm{\Gamma}}^{\xi})^{\prime})^{\prime}$
denote a generic value of the parameters. Whenever we consider
$\underline{\bm{\varphi}}$ it is intended that its elements satisfy
Assumptions 1 through 6.

Then, the elements of $\widehat{\bm{\varphi}}^{\text{\tiny QML}}$ are such
that:

$\widehat{\bm{\Lambda}}^{\text{\tiny QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}=\widehat{\bm{\Gamma}}^{x},$ (A.11)

where $\widehat{\bm{\Gamma}}^{x}=T^{-1}{\bm{X}^{\prime}\bm{X}}$. To see that
(A.11) defines the global maximum of the log-likelihood we proceed in two
steps.

First, notice that the first order conditions derived from the log-likelihood
(12) are satisfied when (A.11) holds:

$\displaystyle\left.\frac{\partial\ell(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\Lambda}}}\right|_{\underline{\bm{\Lambda}}=\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}}$ $\displaystyle=T\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Gamma}}^{x}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}-T\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Lambda}}^{\text{\tiny QML}}$
$\displaystyle=T\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}-T\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Lambda}}^{\text{\tiny QML}}=\mathbf{0}_{r\times
n},$
$\displaystyle\left.\frac{\partial\ell(\bm{\mathcal{X}};\underline{\bm{\varphi}})}{\partial\underline{\bm{\Gamma}}^{\xi}}\right|_{\underline{\bm{\Gamma}}^{\xi}=\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}}$ $\displaystyle=\frac{T}{2}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Gamma}}^{x}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}-\frac{T}{2}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny QML}}\right)^{-1}$
$\displaystyle=\frac{T}{2}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}-\frac{T}{2}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}=\mathbf{0}_{n\times n}.$

Notice also that the conditions given in Bai and Li (2012, Equations
(2.7)-(2.8)), which are derived from the first order conditions above, are
also satisfied. Namely, it holds that:

$\displaystyle\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\left\\{\widehat{\bm{\Gamma}}^{x}-\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}-\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right\\}=\mathbf{0}_{r\times n},$
$\displaystyle\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}=\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Gamma}}^{x}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny QML}}\right)^{-1},$
$\displaystyle\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}=\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Gamma}}^{x}\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)^{-1}\widehat{\bm{\Lambda}}^{\text{\tiny QML}}.$

Second, given the log-likelihood (12), for any $\underline{\bm{\varphi}}$, we
have:

$\displaystyle\ell(\bm{\mathcal{X}};\widehat{\bm{\varphi}}^{\text{\tiny
QML}})$ $\displaystyle-\ell(\bm{\mathcal{X}};\underline{\bm{\varphi}})$
$\displaystyle=-\frac{T}{2}\log\frac{\det\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)}{\det\left(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Gamma}}^{\xi}\right)}-\frac{nT}{2}+\frac{T}{2}\text{tr}\left\\{\widehat{\bm{\Gamma}}^{x}\left(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Gamma}}^{\xi}\right)^{-1}\right\\}$
(A.12)
$\displaystyle=-\frac{T}{2}\log\frac{\det\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)}{\det\left(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Gamma}}^{\xi}\right)}-\frac{nT}{2}+\frac{T}{2}\text{tr}\left\\{\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)\left(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Gamma}}^{\xi}\right)^{-1}\right\\},$

because of (A.11). Now, denote as $\zeta_{j}$, $j=1,\ldots,n$, the $n$ roots
of

$\det\left\\{\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)-\zeta\left(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Gamma}}^{\xi}\right)\right\\}=0,$

which are all real since $\left(\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}+\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)\left(\underline{\bm{\Lambda}}\,\underline{\bm{\Lambda}}^{\prime}+\underline{\bm{\Gamma}}^{\xi}\right)^{-1}$
is a symmetric matrix. Then, (A.12) reads:

$\ell(\bm{\mathcal{X}};\widehat{\bm{\varphi}}^{\text{\tiny
QML}})-\ell(\bm{\mathcal{X}};\underline{\bm{\varphi}})=\frac{T}{2}\sum_{j=1}^{n}\left\\{-\log\zeta_{j}-1+\zeta_{j}\right\\}\geq
0,$ (A.13)

since $x\leq e^{x-1}$ and so $-\log x-1+x\geq 0$. Therefore, from (A.13) we
see that (A.11) defines indeed the global maximum of the log-likleihood (12).

Now, consider the Singular Value Decomposition (SVD) of the true loadings:

$\frac{{\bm{\Lambda}}}{\sqrt{n}}={\mathbf{V}}{\mathbf{D}}{\mathbf{U}},$ (A.14)

where ${\mathbf{V}}$ is $n\times r$ and such that
${\mathbf{V}}^{\prime}{\mathbf{V}}=\mathbf{I}_{r}$ for all $n\in\mathbb{N}$,
${\mathbf{D}}$ is $r\times r$ diagonal with strictly positive entries, and
${\mathbf{U}}$ is $r\times r$ such that
${\mathbf{U}}{\mathbf{U}}^{\prime}={\mathbf{U}}^{\prime}{\mathbf{U}}=\mathbf{I}_{r}$.

We know that under Assumption 5 and 6, $\bm{\Lambda}$ is globally identified.
Let us show that ${\mathbf{V}}$, ${\mathbf{D}}$, and ${\mathbf{U}}$ in (A.14)
are also globally identified. First, notice that, given Assumption 5(a) which
requires $n^{-1}{{\bm{\Lambda}}^{\prime}{\bm{\Lambda}}}$ to be diagonal, in
order to estimate ${\bm{\Lambda}}$ we need to estimate $nr-{r(r-1)}/{2}$
parameters. Then, by looking at the right-hand-side of (A.14) we see that to
estimate ${\mathbf{V}}$ we need to estimate $nr-\frac{r(r+1)}{2}$ parameters
and to estimate ${\mathbf{D}}$ we need to estimate $r$ parameters, thus
${\mathbf{V}}{\mathbf{D}}$ depends on $nr-{r(r+1)}/{2}+r=nr-{r(r-1)}/{2}$
parameters as ${\bm{\Lambda}}$. However, in principle the right-hand-side of
(A.14) depends also on ${\mathbf{U}}$ which in turn requires estimating
${r(r+1)}/2$ parameters more. But if we impose Assumption 5(a) also to the
right-hand-side of (A.14) we have that
$n^{-1}{{\bm{\Lambda}}^{\prime}{\bm{\Lambda}}}={\mathbf{U}}^{\prime}{\mathbf{D}}^{2}{\mathbf{U}}$
has to be diagonal, and since ${\mathbf{D}}$ is diagonal, without loss of
generality we can set ${\mathbf{U}}=\mathbf{J}$, a diagonal $r\times r$ matrix
with entries $\pm 1$.

Furthermore, from Proposition B.1(a) we also see that we must have
$n^{-1}{{\bm{\Lambda}}^{\prime}{\bm{\Lambda}}}=n^{-1}{\mathbf{M}^{\chi}}=\mathbf{D}^{2}$.
Hence,

$\mathbf{D}=\left(\frac{\mathbf{M}^{\chi}}{n}\right)^{1/2},$ (A.15)

and by Lemma 1(iv) the entries of $\mathbf{D}$, denoted as $d_{j}$,
$j=1,\ldots,r$, are such that

$\sqrt{\underline{C}_{j}}\\!\leq\lim\inf_{n\to\infty}d_{j}\leq\lim\sup_{n\to\infty}d_{j}\leq\\!\sqrt{\overline{C}_{j}},$
(A.16)

where $\underline{C}_{j}$ and $\overline{C}_{j}$ are finite positive reals.
Last, from (A.14) and Proposition B.1(b), we must have
$n^{-1/2}{{\bm{\Lambda}}}={\mathbf{V}}{\mathbf{D}}\mathbf{J}=\mathbf{V}^{\chi}\left(n^{-1}{\mathbf{M}^{\chi}}\right)^{1/2}$,
and by (A.15) it follows that

${\mathbf{V}}=\mathbf{V}^{\chi}\mathbf{J}.$ (A.17)

This shows that under Assumption 5(a), the parameters in ${\mathbf{D}}$ are
globally identified while the parameters in $\mathbf{V}$ are uniquely
identified up to a right-multiplication by $\mathbf{J}$, which can be pinned
down by means of Assumption 6, thus achieving global identification of
$\mathbf{V}$ as well. Given this discussion, hereafter, we can directly set
$\mathbf{U}=\mathbf{I}_{r}$.

It is clear that the problem of QML estimation of the loadings can be
rewritten as a problem of QML estimation of their SVD in (A.14), namely of
$\mathbf{V}$ and $\mathbf{D}$. To this end, we introduce also the SVDs of the
QML estimator of the loadings and of a generic value of the loadings:

$\displaystyle\frac{\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}}{\sqrt{n}}=\widehat{\mathbf{V}}^{\text{\tiny
QML}}\widehat{\mathbf{D}}^{\text{\tiny QML}}\widehat{\mathbf{U}}^{\text{\tiny
QML}},\qquad\frac{\underline{\bm{\Lambda}}}{\sqrt{n}}=\underline{\mathbf{V}}\,\underline{\mathbf{D}}\,\underline{\mathbf{U}},$
(A.18)

where $\widehat{\mathbf{V}}^{\text{\tiny QML}}$ and $\underline{\mathbf{V}}$
have the same properties as ${\mathbf{V}}$, $\widehat{\mathbf{D}}^{\text{\tiny
QML}}$ and $\underline{\mathbf{D}}$ have the same properties as
${\mathbf{D}}$, and $\widehat{\mathbf{U}}^{\text{\tiny QML}}$ and
$\underline{\mathbf{U}}$ have the same properties as ${\mathbf{U}}$.

Because we set $\mathbf{U}=\mathbf{I}_{r}$, it follows that we can set
$\widehat{\mathbf{U}}^{\text{\tiny QML}}=\mathbf{I}_{r}$, and we are left with
the task of finding $\widehat{\mathbf{V}}^{\text{\tiny QML}}$ and
$\widehat{\mathbf{D}}^{\text{\tiny QML}}$. From (A.11) and since we must have
$\widehat{\mathbf{V}}^{\text{\tiny
QML}\prime}\widehat{\mathbf{V}}^{\text{\tiny QML}}=\mathbf{I}_{r}$, it follows
that

$\frac{\widehat{\bm{\Gamma}}^{x}}{n}-\widehat{\mathbf{V}}^{\text{\tiny
QML}}\left(\widehat{\mathbf{D}}^{\text{\tiny
QML}}\right)^{2}\widehat{\mathbf{V}}^{\text{\tiny
QML}\prime}-\frac{\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}}{n}=\mathbf{0}_{n\times n},$

which is equivalent to

$\frac{\widehat{\mathbf{V}}^{\text{\tiny
QML}\prime}\widehat{\bm{\Gamma}}^{x}\widehat{\mathbf{V}}^{\text{\tiny
QML}}}{n}-\left(\widehat{\mathbf{D}}^{\text{\tiny
QML}}\right)^{2}-\frac{\widehat{\mathbf{V}}^{\text{\tiny
QML}\prime}\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\widehat{\mathbf{V}}^{\text{\tiny QML}}}{n}=\mathbf{0}_{r\times r}.$
(A.19)

Now, let $\widehat{\bm{v}}_{j}^{\text{\tiny QML}}$, $j=1,\ldots,r$, be the
$j$th column of $\widehat{\mathbf{V}}^{\text{\tiny QML}}$ and let
$\widehat{d}_{j}^{\,\text{\tiny QML}}$, $j=1,\ldots,r$, be the $j$th diagonal
entry of $\widehat{\mathbf{D}}^{\text{\tiny QML}}$. Then, from (A.19), for all
$j=1,\ldots,r$, we have

$\frac{\widehat{\bm{v}}_{j}^{\text{\tiny
QML}\prime}\widehat{\bm{\Gamma}}^{x}\widehat{\bm{v}}_{j}^{\text{\tiny
QML}}}{n}-\left(\widehat{d}_{j}^{\,\text{\tiny
QML}}\right)^{2}-\frac{\widehat{\bm{v}}_{j}^{\text{\tiny
QML}\prime}\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\widehat{\bm{v}}_{j}^{\text{\tiny QML}}}{n}=0.$ (A.20)

Then, trivially, the QML estimators are such that:

$\left(\widehat{\bm{v}}_{j}^{\text{\tiny QML}},\widehat{d}_{j}^{\,\text{\tiny
QML}},\widehat{\bm{\Gamma}}^{\xi,\text{\tiny
QML}}\right)=\arg\\!\\!\\!\\!\min_{\underline{\bm{v}}_{j},\underline{d}_{j},\underline{\bm{\Gamma}}^{\xi}}\left(\frac{\underline{\bm{v}}_{j}^{\prime}\widehat{\bm{\Gamma}}^{x}\underline{\bm{v}}_{j}}{n}-\underline{d}_{j}^{2}-\frac{\underline{\bm{v}}_{j}^{\prime}\underline{\bm{\Gamma}}^{\xi}\underline{\bm{v}}_{j}}{n}\right)^{2}=\arg\\!\\!\\!\\!\min_{\underline{\bm{v}}_{j},\underline{d}_{j},\underline{\bm{\Gamma}}^{\xi}}\mathcal{L}_{0j}\left(\underline{\bm{v}}_{j},\underline{d}_{j},\underline{\bm{\Gamma}}^{\xi}\right),\;\text{say,}$
(A.21)

where $\underline{\bm{v}}_{j}$, $j=1,\ldots,r$, is the $j$th column of
$\underline{\mathbf{V}}$ and $\underline{d}_{j}$, $j=1,\ldots,r$, is the $j$th
diagonal entry of $\underline{\mathbf{D}}$.

Define also the estimators $\widetilde{\bm{v}}_{j}$ and $\widetilde{d}_{j}$
such that:

$\left(\widetilde{\bm{v}}_{j},\widetilde{d}_{j}\right)=\arg\\!\min_{\underline{\bm{v}}_{j},\underline{d}_{j}}\left(\frac{\underline{\bm{v}}_{j}^{\prime}\widehat{\bm{\Gamma}}^{x}\underline{\bm{v}}_{j}}{n}-\underline{d}_{j}^{2}\right)^{2}=\arg\\!\min_{\underline{\bm{v}}_{j},\underline{d}_{j}}\mathcal{L}_{1j}\left(\underline{\bm{v}}_{j},\underline{d}_{j}\right),\;\text{say.}$
(A.22)

Consistently with (A.14) and (A.18), these define an estimator of
$\bm{\Lambda}$ by means of its SVD:

$\frac{\widetilde{\bm{\Lambda}}}{\sqrt{n}}=\widetilde{\mathbf{V}}\widetilde{\mathbf{D}}\widetilde{\mathbf{U}},$
(A.23)

where $\widetilde{\mathbf{V}}$ has columns $\widetilde{\bm{v}}_{j}$,
$j=1,\ldots,r$, and it has the same properties as ${\mathbf{V}}$ (because
$\underline{\mathbf{V}}$ does), and $\widetilde{\mathbf{D}}$ has entries
$\widetilde{d}_{j}$, $j=1,\ldots,r$, and it has the same properties as
${\mathbf{D}}$ (because $\underline{\mathbf{D}}$ does). We set
$\widetilde{\mathbf{U}}=\mathbf{I}_{r}$ consistently with the fact that
$\mathbf{U}=\mathbf{I}_{r}$.

Then, it is easily seen that

$\mathcal{L}_{0j}\left(\underline{\bm{v}}_{j},\underline{d}_{j},\underline{\bm{\Gamma}}^{\xi}\right)=\mathcal{L}_{1j}\left(\underline{\bm{v}}_{j},\underline{d}_{j}\right)+\left(\frac{\underline{\bm{v}}_{j}^{\prime}\underline{\bm{\Gamma}}^{\xi}\underline{\bm{v}}_{j}}{n}\right)^{2}-2\left(\frac{\underline{\bm{v}}_{j}^{\prime}\underline{\bm{\Gamma}}^{\xi}\underline{\bm{v}}_{j}}{n}\right)\sqrt{\mathcal{L}_{1j}\left(\underline{\bm{v}}_{j},\underline{d}_{j}\right)}.$
(A.24)

Now, we know that for any generic value of the idiosyncratic covariance matrix
satisfying Assumption 2(b), it holds that

$\frac{\underline{\bm{v}}_{j}^{\prime}\underline{\bm{\Gamma}}^{\xi}\underline{\bm{v}}_{j}}{n}\leq\max_{\bm{w}\,:\,\bm{w}^{\prime}\bm{w}=1}\frac{\bm{w}^{\prime}\underline{\bm{\Gamma}}^{\xi}\bm{w}}{n}=\frac{\underline{\mu}^{\xi}_{1}}{n}\leq\frac{M_{2\xi}}{n},$
(A.25)

since $\underline{\bm{v}}_{j}^{\prime}\underline{\bm{v}}_{j}=1$ and because of
Lemma 1(v) which implies also that $M_{2\xi}$ is independent of $n$. Moreover,
for any generic $\underline{\bm{\Lambda}}$ satisfying Assumption 1(a):

$\sqrt{\mathcal{L}_{1j}\left(\underline{\bm{v}}_{j},\underline{d}_{j}\right)}=\left|\frac{\underline{\bm{v}}_{j}^{\prime}\widehat{\bm{\Gamma}}^{x}\underline{\bm{v}}_{j}}{n}-\underline{d}_{j}^{2}\right|\leq\frac{\underline{\bm{v}}_{j}^{\prime}\widehat{\bm{\Gamma}}^{x}\underline{\bm{v}}_{j}}{n}+\underline{d}_{j}^{2}\leq\frac{\widehat{\mu}_{1}^{x}}{n}+\underline{d}_{j}^{2}=O_{\mathrm{P}}(1),$
(A.26)

since $\underline{d}_{j}=O(1)$ because
$\|n^{-1/2}{\underline{\bm{\Lambda}}}\|=O(1)$ by Lemma 2(i), and
$\widehat{\mu}_{1}^{x}=O_{\mathrm{P}}(n)$ because of Lemma 8(iii). Therefore,
from (A.24), (A.25), and (A.26), we have

$\left|\mathcal{L}_{0j}\left(\underline{\bm{v}}_{j},\underline{d}_{j},\underline{\bm{\Gamma}}^{\xi}\right)-\mathcal{L}_{1j}\left(\underline{\bm{v}}_{j},\underline{d}_{j}\right)\right|=O_{\mathrm{P}}\left(\frac{1}{n}\right),$
(A.27)

which holds for any $\underline{\bm{\Gamma}}^{\xi}$ satisfying Assumption
2(b). By continuity of these loss functions, from (A.27) it follows that their
minima satisfy:

$\displaystyle\left\|\widehat{\bm{v}}_{j}^{\text{\tiny
QML}}-\widetilde{\bm{v}}_{j}\right\|=O_{\mathrm{P}}\left(\frac{1}{n}\right),\qquad\left|\widehat{d}_{j}^{\,\text{\tiny
QML}}-\widetilde{d}_{j}\right|=O_{\mathrm{P}}\left(\frac{1}{n}\right),$ (A.28)

for all $j=1,\ldots,r$.

Let us now find $\widetilde{\bm{v}}_{j}$ and $\widetilde{d}_{j}$. From (A.22),
it is clear that the solutions must be such that:

$\frac{\widetilde{\bm{v}}_{j}^{\prime}\widehat{\bm{\Gamma}}^{x}\widetilde{\bm{v}}_{j}}{n}=\widetilde{d}_{j}^{\,2},$
(A.29)

which means that $\widetilde{d}_{j}^{\,2}$ must be an eigenvalue of
$n^{-1}{\widehat{\bm{\Gamma}}^{x}}$ and $\widetilde{\bm{v}}_{j}$ is the
corresponding normalized eigenvector. Obviously, the solution in (A.31)
defines a global minimum of the loss
$\mathcal{L}_{1}(\underline{\bm{v}}_{j},\underline{d}_{j})$ since
$\mathcal{L}_{1}(\widetilde{\bm{v}}_{j},\widetilde{d}_{j})=0$ while
$\mathcal{L}_{1}(\underline{\bm{v}}_{j},\underline{d}_{j})>0$ for any other
value $\underline{\bm{v}}_{j}\neq\widetilde{\bm{v}}_{j}$ and
$\underline{d}_{j}\neq\widetilde{d}_{j}$.

Now, let us show that indeed it must be that
$\widetilde{d}_{j}^{\,2}=n^{-1}{\widehat{\mu}_{j}^{x}}$, $j=1,\ldots,r$, i.e.,
they have to be the $r$ largest eigenvalues of
$n^{-1}{\widehat{\bm{\Gamma}}^{x}}$. First, by Lemma 6(i) and Weyl’s
inequality, for all $k=1,\ldots,r$, as $n,T\to\infty$,

$\left|\frac{\widehat{\mu}_{k}^{x}}{n}-\frac{{\mu}_{k}^{x}}{n}\right|\leq\left\|\frac{\widehat{\bm{\Gamma}}^{x}}{n}-\frac{{\bm{\Gamma}}^{x}}{n}\right\|=O_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right).$
(A.30)

Therefore, if for any given $j=1,\ldots,r$ we were to choose
$\widetilde{d}_{j}^{\,2}=n^{-1}{\widehat{\mu}_{k}^{x}}$ for, say, $k=r+1$,
then, from (A.30), we would have

$\displaystyle\lim\inf_{n\to\infty}\widetilde{d}_{j}^{\,2}=\lim\inf_{n\to\infty}\frac{\widehat{\mu}_{r+1}^{x}}{n}=\lim\inf_{n\to\infty}\frac{{\mu}_{r+1}^{x}}{n}+O_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right)\geq\lim\inf_{n\to\infty}\frac{{\mu}_{n}^{\xi}}{n}+O_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right)=O_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right),$
$\displaystyle\lim\sup_{n\to\infty}\widetilde{d}_{j}^{\,2}=\lim\sup_{n\to\infty}\frac{\widehat{\mu}_{r+1}^{x}}{n}=\lim\sup_{n\to\infty}\frac{{\mu}_{r+1}^{x}}{n}+O_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right)\leq\lim_{n\to\infty}\frac{M_{\xi}}{n}+O_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right)=O_{\mathrm{P}}\left(\frac{1}{\sqrt{T}}\right),$

since we assumed $\bm{\Gamma}^{\xi}$ to be positive definite, so
$\mu^{\xi}_{n}>0$ for all $n\in\mathbb{N}$, and by Lemma 1(vi). Therefore, as
$n,T\to\infty$, this choice for $\widetilde{d}_{j}$ cannot be a consistent
estimator of $d_{j}$ (nor an approximation of $\widehat{d}_{j}^{\,\text{\tiny
QML}}$), since while $\widetilde{d}_{j}\to 0$, as $n,T\to\infty$, it must be
that $d_{j}>0$ for all $n\in\mathbb{N}$, as required in (A.16).

So from (A.29) and the above reasoning it follows that, for all
$j=1,\ldots,r$,

$\widetilde{d}_{j}^{\,2}=\frac{\widehat{\mu}_{j}^{x}}{n},\qquad\widetilde{\bm{v}}_{j}=\widehat{\bm{v}}_{j}^{x}.$
(A.31)

where $\widehat{\bm{v}}_{j}^{x}$ is the eigenvector of
$n^{-1}{\widehat{\bm{\Gamma}}^{x}}$ corresponding to its $j$th largest
eigenvalue.

Then, from (A.31), first we have

$\displaystyle\left\|\widehat{\mathbf{V}}^{\text{\tiny
QML}}-\widehat{\mathbf{V}}^{x}\right\|=\left\|\widehat{\mathbf{V}}^{\text{\tiny
QML}}-\widetilde{\mathbf{V}}\right\|=O_{\mathrm{P}}\left(\frac{1}{n}\right),$
(A.32) $\displaystyle\left\|\widehat{\mathbf{D}}^{\,\text{\tiny
QML}}-\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{1/2}\right\|=\left\|\widehat{\mathbf{D}}^{\,\text{\tiny
QML}}-\widetilde{\mathbf{D}}\right\|=O_{\mathrm{P}}\left(\frac{1}{n}\right),$
(A.33)

because of (A.28), and, second, we have

$\frac{\widetilde{\bm{\Lambda}}}{\sqrt{n}}=\widehat{\mathbf{V}}^{x}\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{1/2}=\frac{\widehat{\bm{\Lambda}}}{\sqrt{n}},$
(A.34)

which shows that the estimator $\widetilde{\bm{\Lambda}}$ is the PC estimator
defined in (9). Therefore, from (A.32), (A.33), and (A.34), and by using the
SVD of the QML estimator in (A.18), it follows that:

$\displaystyle\left\|\frac{\widehat{\bm{\Lambda}}^{\text{\tiny
QML}}}{\sqrt{n}}-\frac{\widehat{\bm{\Lambda}}}{\sqrt{n}}\right\|=$
$\displaystyle\,\left\|\widehat{\mathbf{V}}^{\text{\tiny
QML}}\widehat{\mathbf{D}}^{\text{\tiny
QML}}-\widehat{\mathbf{V}}^{x}\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{1/2}\right\|$
$\displaystyle\leq$ $\displaystyle\,\left\|\widehat{\mathbf{V}}^{\text{\tiny
QML}}-\widehat{\mathbf{V}}^{x}\right\|\,\left\|\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{1/2}\right\|+\left\|\widehat{\mathbf{D}}^{\text{\tiny
QML}}-\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{1/2}\right\|\,\left\|\widehat{\mathbf{V}}^{x}\right\|$
$\displaystyle+\left\|\widehat{\mathbf{V}}^{\text{\tiny
QML}}-\widehat{\mathbf{V}}^{x}\right\|\,\left\|\widehat{\mathbf{D}}^{\text{\tiny
QML}}-\left(\frac{\widehat{\mathbf{M}}^{x}}{n}\right)^{1/2}\right\|=O_{\mathrm{P}}\left(\frac{1}{n}\right),$

because $\|\widehat{\mathbf{V}}^{x}\|=1$ (it is a matrix of normalized
eigenvectors) and
$\|(n^{-1}{\widehat{\mathbf{M}}^{x}})^{1/2}\|=O_{\mathrm{P}}(1)$ by Lemma
8(iii). This proves part (a).

Finally, consider the $r$ dimensional $i$th rows, $i=1,\ldots,n$, of
${\mathbf{V}}^{\chi}$, $\widehat{\mathbf{V}}^{x}$, $\underline{\mathbf{V}}$,
and $\widehat{\mathbf{V}}^{\text{\tiny QML}}$, denoted as
${\mathbf{v}}_{i}^{\chi\prime}$, $\widehat{\mathbf{v}}_{i}^{x\prime}$,
$\underline{\mathbf{v}}_{i}^{\prime}$ , and
$\widehat{\mathbf{v}}_{i}^{\text{\tiny QML}\prime}$, respectively. From Lemma
7(ii) and 7(iii), we know that
$\sqrt{n}\|{\mathbf{v}}_{i}^{\chi\prime}\|=O(1)$ and
$\sqrt{n}\|\widehat{\mathbf{v}}_{i}^{x\prime}\|=O_{\mathrm{P}}(1)$, which
means that we must have $\sqrt{n}\|\underline{\mathbf{v}}_{i}^{\prime}\|=O(1)$
since any generic parameter we consider is assumed to satisfy the same
assumptions as the corresponding true parameters
(${\mathbf{v}}_{i}^{\chi\prime}$ in this case). Therefore, since the search of
the elements of $\widehat{\mathbf{V}}^{\text{\tiny QML}}$ is made over all
$\underline{\mathbf{V}}$ satisfying the same assumptions as
${\mathbf{V}}^{\chi}$, it follows that we must have also
$\sqrt{n}\|\widehat{\mathbf{v}}_{i}^{\text{\tiny
QML}\prime}\|=O_{\mathrm{P}}(1)$.

To see that this is the case, notice also that because of Assumption 5(a) the
QML estimator must be such that $n^{-1}{\widehat{\bm{\Lambda}}^{\text{\tiny
QML}\prime}\widehat{\bm{\Lambda}}^{\text{\tiny QML}}}$ is diagonal and
positive definite. Thus, we must have
$\|n^{-1/2}{\widehat{\bm{\Lambda}}^{\text{\tiny QML}}}\|=O_{\mathrm{P}}(1)$
(see also Lemma 2(i)). From the SVD in (A.18) when
$\widehat{\mathbf{U}}^{\text{\tiny QML}}=\mathbf{I}_{r}$ and since it must be
that $\widehat{\mathbf{V}}^{\text{\tiny
QML}\prime}\widehat{\mathbf{V}}^{\text{\tiny QML}}=\mathbf{I}_{r}$ for all
$n\in\mathbb{N}$, then it follows that $\|\widehat{\mathbf{V}}^{\text{\tiny
QML}}\|=O_{\mathrm{P}}(1)$ and $\|\widehat{\mathbf{D}}^{\text{\tiny
QML}}\|=O_{\mathrm{P}}(\sqrt{n})$. Moreover, since by Assumption 1(a) the
$i$th row of $\widehat{\bm{\Lambda}}^{\text{\tiny QML}}$ must be such that
$\|\widehat{\bm{\lambda}}_{i}^{\text{\tiny QML}}\|=O_{\mathrm{P}}(1)$ and
$\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML}\prime}=\widehat{\mathbf{v}}_{i}^{\text{\tiny
QML}\prime}\widehat{\mathbf{D}}^{\text{\tiny QML}}$, then it must be that
$\sqrt{n}\|\widehat{\mathbf{v}}_{i}^{\text{\tiny
QML}\prime}\|=O_{\mathrm{P}}(1)$.

From this reasoning and (A.32) it follows that:

$\sqrt{n}\left\|\widehat{\mathbf{v}}_{i}^{\text{\tiny
QML}}-\widehat{\mathbf{v}}_{i}^{x}\right\|=O_{\mathrm{P}}\left(\frac{1}{n}\right).$
(A.35)

Now, by taking the $i$th row of the PC estimator in (A.34) and the SVD of the
QML estimator in (A.18), for any $i=1,\ldots,n$, we get

$\displaystyle\left\|\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML}\prime}-\widehat{\bm{\lambda}}_{i}^{\prime}\right\|=$
$\displaystyle\,\left\|\widehat{\mathbf{v}}_{i}^{\text{\tiny
QML}\prime}\sqrt{n}\widehat{\mathbf{D}}^{\text{\tiny
QML}}-\widehat{\mathbf{v}}_{i}^{x\prime}(\widehat{\mathbf{M}}^{x})^{1/2}\right\|$
$\displaystyle\leq$
$\displaystyle\,\left\|\widehat{\mathbf{v}}_{i}^{\text{\tiny
QML}}-\widehat{\mathbf{v}}_{i}^{x}\right\|\,\left\|\left({\widehat{\mathbf{M}}^{x}}\right)^{1/2}\right\|+\left\|\sqrt{n}\widehat{\mathbf{D}}^{\text{\tiny
QML}}-\left({\widehat{\mathbf{M}}^{x}}\right)^{1/2}\right\|\,\left\|\widehat{\mathbf{v}}_{i}^{x}\right\|$
$\displaystyle+\left\|\widehat{\mathbf{v}}_{i}^{\text{\tiny
QML}}-\widehat{\mathbf{v}}_{i}^{x}\right\|\,\left\|\sqrt{n}\widehat{\mathbf{D}}^{\text{\tiny
QML}}-\left({\widehat{\mathbf{M}}^{x}}\right)^{1/2}\right\|=O_{\mathrm{P}}\left(\frac{1}{n}\right),$

because of (A.33), (A.35), and since
$\|\widehat{\mathbf{v}}_{i}^{x}\|=O_{\mathrm{P}}(n^{-1/2})$ by Lemma 7(ii) and
7(iii), and $\|({\widehat{\mathbf{M}}^{x}})^{1/2}\|=O_{\mathrm{P}}(\sqrt{n})$
by Lemma 8(iii). This proves part (b) and completes the proof. $\Box$

###  A.4 Proof of Theorem 4

From Theorem 3(b) and Corollary 1 it follows that:

$\displaystyle\|\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML}}-{\bm{\lambda}}_{i}\|$
$\displaystyle=\|\widehat{\bm{\lambda}}_{i}-{\bm{\lambda}}_{i}\|+O_{\mathrm{P}}\left(\frac{1}{n}\right)$
$\displaystyle=\|\widehat{\bm{\lambda}}_{i}-{\bm{\lambda}}_{i}^{\text{\tiny
OLS}}\|+\|{\bm{\lambda}}_{i}^{\text{\tiny
OLS}}-{\bm{\lambda}}_{i}\|+O_{\mathrm{P}}\left(\frac{1}{n}\right)$
$\displaystyle=\|{\bm{\lambda}}_{i}^{\text{\tiny
OLS}}-{\bm{\lambda}}_{i}\|+O_{\mathrm{P}}\left(\max\left(\frac{1}{n},\frac{1}{\sqrt{nT}}\right)\right)$
$\displaystyle=O_{\mathrm{P}}\left(\max\left(\frac{1}{n},\frac{1}{\sqrt{nT}},\frac{1}{\sqrt{T}}\right)\right),$
(A.36)

since the unfeasible OLS estimator is $\sqrt{T}$-consistent. Moreover, if
$\sqrt{T}/n\to 0$ as $n,T\to\infty$, by imposing the identification constraint
of orthonormal factors in Assumption 5(b), we have

$\displaystyle\sqrt{T}(\widehat{\bm{\lambda}}_{i}^{\text{\tiny
QML}}-\bm{\lambda}_{i})$
$\displaystyle=\sqrt{T}(\widehat{\bm{\lambda}}_{i}^{\text{\tiny
OLS}}-\bm{\lambda}_{i})+o_{\mathrm{P}}(1)=\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathbf{F}_{t}\xi_{it}+o_{\mathrm{P}}(1)\to_{d}\mathcal{N}\left(\mathbf{0}_{r},\bm{\Phi}_{i}\right).$

by Slutsky’s theorem and Assumption 7. This completes the proof of part (a).
Part (b) follows similarly from Theorem 3(a). This completes the proof. $\Box$

###  A.5 Proof of Corollary 2