# Latent Factor Analysis in Short Panels

Alain-Philippe Fortin1,2, Patrick Gagliardini3,2, Olivier Scaillet1,2
1University of Geneva, 2Swiss Finance Institute, 3Università della Svizzera
italiana. Acknowledgements: We are grateful to A. Onatski for his very
insightful discussion (Onatski (2023)) of our paper Fortin, Gagliardini,
Scaillet (2023) at the 14th Annual SoFiE Conference in Cambridge, which
prompted us to exploit factor analysis for estimation in short panels. We
thank D. Amengual, L. Barras, S. Bonhomme, M. Caner, F. Carlini, I. Chaieb, F.
Ghezzi, A. Horenstein, G. Imbens, S. Kim, F. Kleibergen, H. Langlois, T.
Magnac, S. Ng, E. Ossola, Y. Potiron, F. Trojani, participants at (EC)^2 2022,
QFFE 2023, SoFiE 2023, NASM 2023, IPDC 2023, FinEML, 22${}^{\text{\\`{e}me}}$
Journée d’Économétrie, HKUST workshop, SFI research days, GIGS conference, and
seminars at UNIMIB, UNIBE, UNIGE, Warwick, Bristol, QMUL, CUHK, and
Luxembourg. The first and third author also acknowledge financial support by
the Swiss National Science Foundation (grants UN11140 and 100018_215573).

Abstract

We develop inferential tools for latent factor analysis in short panels. The
pseudo maximum likelihood setting under a large cross-sectional dimension $n$
and a fixed time series dimension $T$ relies on a diagonal $T\times T$
covariance matrix of the errors without imposing sphericity nor Gaussianity.
We outline the asymptotic distributions of the latent factor and error
covariance estimates as well as of an asymptotically uniformly most powerful
invariant (AUMPI) test for the number of factors based on the likelihood ratio
statistic. We derive the AUMPI characterization from inequalities ensuring the
monotone likelihood ratio property for positive definite quadratic forms in
normal variables. An empirical application to a large panel of monthly U.S.
stock returns separates month after month systematic and idiosyncratic risks
in short subperiods of bear vs. bull market based on the selected number of
factors. We observe an uptrend in the paths of total and idiosyncratic
volatilities while the systematic risk explains a large part of the cross-
sectional total variance in bear markets but is not driven by a single factor.
Rank tests show that observed factors struggle spanning latent factors with a
discrepancy between the dimensions of the two factor spaces decreasing over
time.

Keywords: Latent factor analysis, uniformly most powerful invariant test,
panel data, large $n$ and fixed $T$ asymptotics, equity returns.
JEL codes: C12, C23, C38, C58, G12.

##  1 Introduction

Latent variable models have been used for a long time in econometrics (Aigner
et al. (1984)). Here, we study large cross-sectional latent factor models with
small time dimension. Two common methods for estimation of latent factor
spaces are principal component analysis (PCA) and factor analysis (FA), see
Anderson (2003) Chapters 11 and 14. They cover multiple applications in
finance and economics as well as in social sciences in general. They are often
used in exploratory analysis of data. Omitted latent factors are also called
interactive fixed effects in the panel literature (Pesaran (2006), Bai (2009),
Moon and Weidner (2015), Freyberger (2018)). In recent work, Fortin,
Gagliardini and Scaillet (FGS, 2023) show how we can use PCA to conduct
inference on the number of factors in such models without making Gaussian
assumptions. Their methodology relies on sphericity of the idiosyncratic
variances since this restriction is both necessary and sufficient for
consistency of latent factor estimates with small $T$ (Theorem 4 of Bai
(2003)). FGS provides a discussion of the (in)consistency of the PCA estimator
with fixed T from the vantage point of the well-known incidental parameter
problem of the panel data literature (Neyman and Scott(1948); see Lancaster
(2000) for a review). In PCA, sphericity allows to identify the number $k$ of
factors from the $k$ first eigenvalue spacings being larger than zero, and
being zero the subsequent ones. On the contrary, the FA strategy does not
exploit eigenvalue spacings and does not require sphericity. However,
inference with small $T$ up to now mostly relies on (often restrictive)
assumptions such as Gaussian variables (with a notable exception by Anderson
and Amemiya (1988)) and error homoskedasticity across sample units. Those are
untenable assumptions in our application with stock returns. The strong
assumption of sphericity might also fail to hold in some samples. If it
happens, our Monte Carlo experiments with non-Gaussian errors show that the
eigenvalue spacing test of FGS exhibit size distortions of over 80 percentage
points for a nominal size of 5%.111That massive over-rejection translates into
an average estimated number of factors often above 10 instead of 2 (true
value) in our simulation design when $T=24$ (a sample size close to $T=20$ in
our empirics). We get similar Monte Carlo results for a constrained Likelihood
Ratio test when sphericity does not hold (see Section 4.3). Our Monte Carlo
experiments also reveal that the classical chi-square test of the FA theory
obtained under cross-sectionally homoskedastic Gaussian errors (Anderson
(1963)) suffer from massive over-rejection by around 80 percentage points for
a nominal size of 5% when $T=24$.

A central and practical issue in applied work with latent factors is to
determine the number of factors. For models with unobservable (latent) factors
only, Connor and Korajczyk (1993) are the first to develop a test for the
number of factors for large balanced panels of individual stock returns in
time-invariant models under covariance stationarity and homoskedasticity.
Unobservable factors are estimated by the method of asymptotic principal
components developed by Connor and Korajczyk (1986) (see also Stock and Watson
(2002)). For heteroskedastic settings, the recent literature on large balanced
panels with static factors has extended the toolkit available to researchers.
A first strand of that literature focuses on consistent estimation procedures
for the number of factors. Bai and Ng (2002) introduce a penalized least-
squares strategy to estimate the number of factors, at least one. Ando and Bai
(2015) extend that approach when explanatory variables are present in the
linear specification (see Bai (2009) for homogeneous regression coefficients).
Onatski (2010) looks at the behavior of differences in adjacent eigenvalues to
determine the number of factors when $n$ and $T$ are both large and
comparable. Ahn and Horenstein (2013) opt for a similar strategy based on
eigenvalue ratios. Caner and Han (2014) propose an estimator with a group
bridge penalization to determine the number of unobservable factors. Based on
the framework of Gagliardini, Ossola and Scaillet (2016), Gagliardini, Ossola
and Scaillet (2019) build a simple diagnostic criterion for approximate factor
structure in large panel datasets. Given observable factors, the criterion
checks whether the errors are weakly cross-sectionally correlated, or share
one or more unobservable common factors (interactive effects), and selects
their number; see Gagliardini, Ossola and Scaillet (2020) for a survey of
estimation of large dimensional conditional factor models in finance. A second
strand of that literature develops inference procedures for hypotheses on the
number of latent factors. Onatski (2009) deploys a characterization of the
largest eigenvalues of a Wishart-distributed covariance matrix with large
dimensions in terms of the Tracy-Widom Law. To get a Wishart distribution,
Onatski (2009) assumes either Gaussian errors, or $T$ much larger than $n$.
Kapetanios (2010) uses subsampling to estimate the limit distribution of the
adjacent eigenvalues.

This paper puts forward methodological and empirical contributions that
complement the above literature. (i) On the methodological side, we extend the
inferential tools of FA to non-Gaussian and non-i.i.d. settings. First, we
characterize the asymptotic distribution of FA estimators obtained under a
pseudo maximum likelihood approach where the time-series dimension is held
fixed while the cross-sectional dimension diverges. Hence, the asymptotic
analysis targets short panels, and allows for cross-sectionally
heteroskedastic and weakly dependent errors. Cochrane (2005, p. 226) argues in
favour of the development of appropriate large-$n$ small-$T$ tools for
evaluating asset pricing models, a problem only partially addressed in
finance. In a short panel setting, Zaffaroni (2019) considers inference for
latent factors in conditional linear asset pricing models under sphericity
based on PCA, including estimation of the number of factors.222Raponi, Robotti
and Zaffaroni (2020) develop tests of beta-pricing models and a two-pass
methodology to estimate the ex-post risk premia (Shanken (1992)) associated to
observable factors (see Kleibergen and Zhan (2023) for robust-identification
inference based a continuous updating generalized method of moments). Kim and
Skoulakis (2018) deals with the error-in-variable problem of the two-pass
methodology with small $T$ by regression-calibration under sphericity and a
block-dependence structure. The small $T$ setting mitigates concerns for panel
unbalancedness (outside the straightforward missing-at-random mechanism) and
corresponds to a locally time-invariant factor structure accommodating
globally time-dependent features of general forms. It is also appealing to
macroeconomic data observed quarterly. For the sake of space, we put part of
the theory, namely inference for FA estimates, in the Online Appendix (OA). We
refer to Bai and Li (2016) for inference when $n$ and $T$ are both large (see
Bai and Li (2012) for the cross-sectional independent case). Second, we build
on our new theoretical results for FA to develop testing procedures for the
number of latent factors in a short panel which rely on neither sphericity nor
Gaussianity nor cross-sectional independence, thereby extending tests based on
eigenvalues, as in Onatski (2009), to small $T$, and as in FGS, to non-
spherical errors, thanks to an FA device. Here, we deliver asymptotic feasible
distributional theory even if FA residuals are not consistent estimates of the
true errors under fixed $T$. We further derive the Asymptotically Uniformly
Most Powerful Invariant (AUMPI) property of the FA likelihood ratio (LR) test
statistic in the non-Gaussian case under inequality restrictions on the DGP
parameters, and cover inference with weak factors. The AUMPI property is rare
and sought-after in testing procedures (see Engle (1984) for a discussion and
Romano, Shaikh and Wolf (2010) for a survey of optimality approaches in
testing problems), and often holds only under restrictive assumptions such as
Gaussianity. We show that the AUMPI property can hold even if the asymptotic
distribution of the LR statistic is driven by a weighted average of
independent chi-square variates in our context instead of a single chi-square
variate. We achieve that by providing novel sufficient conditions so that the
ratio of the density under local alternative hypotheses and the null
hypothesis satisfies the Monotone Likelihood Ratio (MLR) property. Hence, the
FA theory gathered in the below consists in a body of new results including
their proofs, and differs completely from PCA theory. (ii) On the empirical
side, we apply our FA methodology to panels of monthly U.S. stock returns with
large cross-sectional and small time-series dimensions, and investigate how
the number of driving factors changes over time and particular periods.
Furthermore, month after month, we provide a novel separation333Such a
decomposition with PCA estimates is invalid without sphericity because of
inconsistency of $\hat{F}$. of the risk coming from the systematic part and
the risk coming from the idiosyncratic part of returns in short subperiods of
bear vs. bull market based on the selected number of factors. We observe an
uptrend in the estimated paths of total and idiosyncratic volatilities (see
also Campbell et al. (2023)) while the systematic risk explains a large part
of the cross-sectional total variance in bear markets but is not driven by a
single latent factor. We also investigate whether standard observed factors
span the estimated latent factors using rank tests in order to suit our fixed
$T$ setting. Observed factors struggle spanning latent factors with a
discrepancy between the dimensions of the two factor spaces decreasing over
time.

The outline of the paper is as follows. In Section 2, we consider a linear
latent factor model and introduce test statistics on the number of latent
factors based on FA. Section 3 presents a feasible asymptotic distributional
theory for inference in short panels under a block-dependence structure to
allow for weak dependence in the cross-section. Section 4 discusses three
special cases, i.e., Gaussian errors (yielding the classical chi-square test
of FA theory in the cross-sectionally homoskedastic case), settings where the
asymptotic distribution under Gaussian errors still holds for the test
statistics, and spherical errors. Section 5 is dedicated to local asymptotic
power and AUMPI tests. We provide our empirical application in Section 6 and
our concluding remarks in Section 7. Appendices A and B gather the regularity
assumptions and proofs of the main theoretical results. Appendix C gives a
Monte Carlo assessment of size and power and selection procedure for the
number of factors for the LR test. We place all omitted proofs and additional
analyses in Appendices D-F in Online Appendix (OA). Appendix G collects the
maximum value of $k$ as a function of $T$. Besides, we gather all explicit
formulas not listed in the core text but useful for coding in an online
“Supplementary Materials for Coding” (SMC) attached to the replication files.
We also put there other numerical checks and additional Monte Carlo results to
assess the impact of non-sphericity for the eigenvalue spacing test of FGS and
the constrained LR test, as well as of using the classical chi-square test on
the size and selection procedure for the number of factors.

##  2 Latent factor model and likelihood ratio test

We consider the linear Factor Analysis (FA) model (e.g. Anderson (2003)):

$y_{i}=\mu+F\beta_{i}+\varepsilon_{i},\qquad i=1,...,n,$ (1)

where $y_{i}=(y_{i,1},...,y_{i,T})^{\prime}$ and
$\varepsilon_{i}=(\varepsilon_{i,1},...,\varepsilon_{i,T})^{\prime}$ are
$T$-dimensional vectors of observed data and unobserved error terms for
individual $i$. The $k$-dimensional vectors
$\beta_{i}=(\beta_{i,1},...,\beta_{i,k})^{\prime}$ are latent individual
effects, while $\mu$ and $F$ are a $T\times 1$ vector and a $T\times k$ matrix
of unknown parameters. The number of latent factors $k$ is an unknown integer
smaller than $T$. In matrix notation, model (1) reads $Y=\mu
1_{n}^{\prime}+F\beta^{\prime}+\varepsilon$, where $Y$ and $\varepsilon$ are
$T\times n$ matrices, $\beta$ is the $n\times k$ matrix with rows
$\beta_{i}^{\prime}$, and $1_{n}$ is a $n$-dimensional vector of ones.

###### Assumption 1

The $T\times T$ matrix
$V_{\varepsilon}=\underset{n\rightarrow\infty}{\lim}E[\frac{1}{n}\varepsilon\varepsilon^{\prime}]$
is diagonal.

Matrix $V_{\varepsilon}$ is the limit cross-sectional average of the -
possibly heterogeneous - errors’ unconditional variance-covariance matrices.
The diagonality condition in Assumption 1 is standard in FA (in the more
restrictive formulation involving i.i.d. data).

In our empirics with a large cross-sectional panel of returns for $n$ assets
over a short time span with $T$ periods, vectors $y_{i}$ and $\varepsilon_{i}$
stack the monthly returns and the idiosyncratic errors of stock $i$. Any row
vector $f_{t}^{\prime}:=(f_{t,1},...,f_{t,k})$ of matrix $F$ yields the latent
factor values in a given month $t$, and vector $\beta_{i}$ collects the factor
loadings of stock $i$. In this finance application, we assume the No-Arbitrage
(NA) principle to hold, so that the entries $\mu_{t}$ in the intercept vector
in Equation (1) account for the (possibly time-varying) risk-free rate and
(possibly non-zero) cross-sectional mean of stock betas.444 Under NA, the
intercept term in the asset return model
$y_{i}=\mu_{i}+\tilde{F}\tilde{\beta}_{i}+\varepsilon_{i}$ is
$\mu_{i}=r_{f}+1_{T}\nu^{\prime}\tilde{\beta}_{i}$, where $r_{f}$ is the
$T$-dimensional vector whose entries collect the (possibly time-varying) risk-
free rates, $\nu=(\nu_{1},...,\nu_{k})^{\prime}$ is a $k$-dimensional vector
of parameters, and $1_{T}$ is a $T$-dimensional vector of ones (see e.g.
Gagliardini, Ossola and Scaillet (2016)). We can absorb term
$1_{T}\nu^{\prime}\tilde{\beta}_{i}$ into the systematic part to get
$y_{i}=r_{f}+F\tilde{\beta}_{i}+\varepsilon_{i}$ with
$F=\tilde{F}+1_{T}\nu^{\prime}$. It holds irrespective of the latent factors
being tradable or not. If the factors are tradable, we further have $\nu=0$
from the NA restriction. Akin to standard formulation of FA, we recenter the
latent effects by subtracting their mean
$\tilde{\mu}_{\tilde{\beta}}=\frac{1}{n}\sum_{i=1}^{n}\tilde{\beta}_{i}$, to
get model (1) with $\beta_{i}=\tilde{\beta}_{i}-\tilde{\mu}_{\tilde{\beta}}$
and $\mu=r_{f}+F\tilde{\mu}_{\tilde{\beta}}$. Thus, the linear FA model (1)
yields $y_{i,t}=\mu_{t}+f_{t}^{\prime}\beta_{i}+\varepsilon_{i,t}$, that is
the standard formulation in asset pricing. We cover the Capital Asset Pricing
Model (CAPM) when the single latent factor is the excess return of the market
portfolio. Assumption 1 allows for serial dependence in idiosyncratic errors
in the form of martingale difference sequences, like individual GARCH and
Stochastic Volatility (SV) processes, as well as weak cross-sectional
dependence (see Assumption 2 below). It also accommodates common time-varying
components in idiosyncratic volatilities by allowing different entries along
the diagonal of $V_{\varepsilon}$; see Renault, Van Der Heijden and Werker
(2023) for arbitrage pricing in such settings.555When there is a common random
component in idiosyncratic volatilities, we have
$V_{\varepsilon}=\underset{n\rightarrow\infty}{plim}\frac{1}{n}\varepsilon\varepsilon^{\prime}$
by a suitable version of the Law of Large Number (LLN) conditional on the
sigma-field generated by this common component. With fixed $T$, we treat the
sample realizations of the common component in idiosyncratic volatilities as
unknown time fixed effects (the diagonal elements of matrix
$V_{\varepsilon}$), which yields time heterogeneous distributions for the
errors. It is how the unconditional expectation in Assumption 1 has to be
understood.

This paper focuses mainly on testing hypotheses on the number of latent
factors $k$ when $T$ is fixed and $n\rightarrow\infty$. The fixed $T$
perspective makes FA especially well-suited for applications with short
panels. Indeed, we work conditionally on the realizations of the latent
factors $F$ and treat their values as parameters to estimate. In comparison
with the standard small $n$ and large $T$ framework in traditional asset
pricing (e.g. Shanken (1992) with observable factors), here factors and
loadings are interchanged in the sense that the $\beta_{i}$ and $F$ play the
roles of the “factors" and the “factor loadings" in FA. We depart from
classical FA since the $\beta_{i}$ are not considered as random effects with a
Gaussian distribution but rather as fixed effects, namely incidental
parameters.666Chamberlain (1992) studies semiparametrically efficient
estimation in panel models with fixed effects and short $T$ using moment
restrictions from instrumental variables. Our approach does not rely on
availability of valid instruments.  Moreover, in Assumption 1, we neither
assume Gaussianity nor we impose sphericity of the covariance matrix of the
error terms. Besides we accommodate weak cross-sectional dependence and ARCH
effects in idiosyncratic errors (see Section 3). Hence, the FA estimators
defined below correspond to maximizers of a Gaussian pseudo likelihood. By-
products of our analysis are the feasible asymptotic distributions of FA
estimators of $F$ and $V_{\varepsilon}$ in more general settings than in the
available literature (e.g. Anderson and Amemiya (1988)), which we present in
Appendix E.

The test statistic we consider for conducting inference on the number of
latent factors $k$ is function of the elements of the symmetric matrix

$\hat{S}=\hat{V}_{\varepsilon}^{-1/2}M_{\hat{F},\hat{V}_{\varepsilon}}(\hat{V}_{y}-\hat{V}_{\varepsilon})M_{\hat{F},\hat{V}_{\varepsilon}}^{\prime}\hat{V}_{\varepsilon}^{-1/2},$
(2)

where $\hat{V}_{y}=\frac{1}{n}\tilde{Y}\tilde{Y}^{\prime}$ is the sample
(cross-sectional) variance matrix (the $n$ columns of $\tilde{Y}$ are
$y_{i}-\bar{y}$ and $\bar{y}=\frac{1}{n}\sum_{i=1}^{n}y_{i}$ is the vector of
cross-sectional means),
$M_{F,V}:=I_{T}-F(F^{\prime}V^{-1}F)^{-1}F^{\prime}V^{-1}$ is the Generalized
Least Squares (GLS) projection matrix orthogonal to $F$ for variance $V$, and
$\hat{F}$ and $\hat{V}_{\varepsilon}$ are the FA estimators computed under the
assumption that there are $k$ latent factors. In the following, we use the
same notation for the matrix-to-vector diag operator and the vector-to-matrix
diag operator. Hence, $diag(A)$ for a matrix $A$ denotes the vector in which
we stack the diagonal elements of matrix $A$, and $diag(a)$ for a vector $a$
denotes a diagonal matrix with the elements of $a$ on the diagonal. From
Anderson (2003) Chapter 14, the FA estimators $\hat{F}$,
$\hat{V}_{\varepsilon}$ maximize a Gaussian pseudo likelihood (Appendix E.1)
and meet the first order conditions:777The normalization in (FA2) applies for
$\hat{\gamma}_{j}\geq 0$, which holds with probability approaching $1$.
Otherwise, the first-order conditions of the FA estimators hold with
$\hat{\gamma}_{j}$ replaced by its positive part; see Anderson (2003) for
similar positivity constraints.

(FA1) $diag(\hat{V}_{y})=diag(\hat{F}\hat{F}^{\prime}+\hat{V}_{\varepsilon})$,
and
(FA2) $\hat{F}$ is the $T\times k$ matrix of eigenvectors of
$\hat{V}_{y}\hat{V}_{\varepsilon}^{-1}$ associated to the $k$ largest
eigenvalues

$1+\hat{\gamma}_{j}$, $j=1,...,k$, normalized such that
$\hat{F}^{\prime}\hat{V}_{\varepsilon}^{-1}\hat{F}=diag(\hat{\gamma}_{1},...,\hat{\gamma}_{k})$.

The number of degrees of freedom is $df=\frac{1}{2}((T-k)^{2}-T-k)$.888Integer
$df$ is the number of different elements in data matrix $\hat{V}_{y}$, i.e.,
$\frac{1}{2}T(T+1)$, plus the number of normalization constraints
$\frac{1}{2}k(k-1)$ in equations $F^{\prime}V_{\varepsilon}^{-1}F=diag$, minus
the number of unknown parameters $(k+1)T$ (Anderson (2003)). In OA Table 3, we
list the largest admissible number $k$ of latent factors as a function of $T$
such that $df\geq 0$. Generalizations to non-diagonal $V_{\varepsilon}$ are
possible (Jöreskog (1970)), but stronger restrictions are required to maintain
nonnegative degrees of freedom. Besides, numerical optimisation becomes less
tractable. It is required that $df\geq 0$ for estimation, and we need $df>0$
to test the null hypothesis of $k$ latent factors (see Proposition 2 (a)
below).

Statistic $\hat{S}$ in Equation (2) checks if the difference between the
sample variance-covariance $\hat{V}_{y}$ and diagonal matrix
$\hat{V}_{\varepsilon}$ is a symmetric matrix of reduced rank $k$, with range
spanned by the range of $\hat{F}$. The probability limit of $\hat{S}$ is nil
under the null hypothesis of $k$ latent factors. We get further insights from
the next result.

###### Proposition 1

Under Assumption 1, (a) the eigenvalues of matrix $\hat{S}$ are:
$\hat{\gamma}_{j}$, for $j=k+1,...,T$, and $0$, with multiplicity $k$, where
$1+\hat{\gamma}_{j}$ for $j=k+1,...,T$ are the $T-k$ smallest eigenvalues of
$\hat{V}_{y}\hat{V}_{\varepsilon}^{-1}$, (b) the squared Frobenius norm is
$\|\hat{S}\|^{2}=\sum_{j=k+1}^{T}\hat{\gamma}_{j}^{2}$, (c) $diag(\hat{S})=0$,
and (d) we get

$\hat{S}=\hat{V}_{\varepsilon}^{-1/2}\left(\frac{1}{n}\hat{\varepsilon}\hat{\varepsilon}^{\prime}\right)\hat{V}_{\varepsilon}^{-1/2}-\hat{V}_{\varepsilon}^{-1/2}M_{\hat{F},\hat{V}_{\varepsilon}}\hat{V}_{\varepsilon}^{1/2},$
(3)

where $\hat{\varepsilon}=M_{\hat{F},\hat{V}_{\varepsilon}}\tilde{Y}$ is the
$T\times n$ matrix of GLS residuals.

From Proposition 1 (d), we can interpret matrix $\hat{S}$ in terms of scaled
cross-sectional averages of squared and cross-products of GLS residuals. In
(3), we subtract
$\hat{V}_{\varepsilon}^{-1/2}M_{\hat{F},\hat{V}_{\varepsilon}}\hat{V}_{\varepsilon}^{1/2}$
and not the identity because the residuals are orthogonal to $\hat{F}$ by
construction. From Proposition 1 (c), the diagonal elements of matrix
$\hat{S}$ vanish. Those elements are not informative for inference on the
number of factors, and can be ignored when constructing the test statistics.
This finding is natural because we expect that only the out-of-diagonal
elements of $\frac{1}{n}\hat{\varepsilon}\hat{\varepsilon}^{\prime}$, i.e.,
the cross-sectional averages of cross-products of residuals for two different
dates, are useful to check for omitted factors.999The test in Connor and
Korajczyk (1993) is built on cross-sectional averages of squared residuals,
akin to diagonal terms of $\hat{S}$, but obtained by PCA instead of FA.
However, their test statistic involves the difference of such cross-sectional
averages for two consecutive dates, and relies on error sphericity.

We now introduce the classical FA likelihood ratio (LR) statistic to test the
null hypothesis $H_{0}(k)$ of $k$ latent factors:

$\displaystyle{LR}(k):=-n\sum_{j=k+1}^{T}\log(1+\hat{\gamma}_{j}).$ (4)

The LR statistic in (4) only uses the information contained in the eigenvalues
of matrix $\hat{S}$. Moreover,
$LR(k)=\frac{n}{2}\|\hat{S}\|^{2}+o_{p}(1)$.101010Alternatively, we can use
the squared norm statistic $\mathscr{T}(k)=n\|\hat{S}\|^{2}$. In the below, we
focus on the LR test statistic since both deliver similar results in Monte
Carlo experiments and in our empirics. It comes from a second-order expansion
of the $\log$ function and the use of $\sum_{j=k+1}^{T}\hat{\gamma}_{j}=0$,
from Proposition 1 (a) and (c) (see also Anderson (2003)), and
$\sqrt{n}\hat{\gamma}_{j}=O_{p}(1)$, for $j=k+1,...,T$, from Propositions 1
and 2. Next, we establish feasibility of the asymptotic distribution of the LR
statistic with $n\rightarrow\infty$ and $T$ fixed under a block-dependence
structure.

##  3 Feasible asymptotic distributional theory

We start by defining the normalization for the latent factor matrix
$F=[F_{1}:\cdots:F_{k}]$. Following classical FA, we set $\bar{\beta}=0$,
$\tilde{V}_{\beta}=I_{k}$, and
$F^{\prime}\tilde{V}_{\varepsilon}^{-1}F=diag(\gamma_{1},...,\gamma_{k})$,
where $\bar{\beta}=\frac{1}{n}\sum_{i=1}^{n}\beta_{i}$,
$\tilde{V}_{\beta}=\frac{1}{n}\beta^{\prime}\beta$, and
$\tilde{V}_{\varepsilon}=E[\frac{1}{n}\varepsilon\varepsilon^{\prime}]$. This
normalization of the factor values is sample dependent, i.e., $F=F_{(n)}$. We
skip index $n$ for the purpose of easing notation. Then, under our
assumptions, we have $V_{y}:=\underset{n\rightarrow\infty}{\text{plim
}}\hat{V}_{y}=FF^{\prime}+V_{\varepsilon}$, with
$F^{\prime}V_{\varepsilon}^{-1}F=diag(\gamma_{1},...,\gamma_{k})$.111111The
standardizations $\bar{\beta}=0$ and $\tilde{V}_{\beta}=I_{k}$ of the factor
loadings wash out the incidental parameter problem (Neyman and Scott (1948);
see Lancaster (2000) for a review) since the individual loadings do not appear
in $V_{y}$. It explains why we are able to get consistent estimators $\hat{F}$
and $\hat{V}_{\varepsilon}$ for large $n$ and fixed $T$. In particular, we
have $V_{y}V_{\varepsilon}^{-1}F_{j}=(1+\gamma_{j})F_{j}$, i.e., the columns
of $F$ are the eigenvectors of matrix $V_{y}V_{\varepsilon}^{-1}$ associated
with eigenvalues $1+\gamma_{j}$, $j=1,...,k$.121212The remaining eigenvalue is
equal to $1$ with multiplicity $T-k$. We have
$F_{j}=\sqrt{\gamma_{j}}V_{\varepsilon}^{1/2}U_{j}$, where the $U_{j}$ are the
orthonormal eigenvectors of
$V_{\varepsilon}^{-1/2}V_{y}V_{\varepsilon}^{-1/2}$ for the $k$ largest
eigenvalues $1+\gamma_{j}$.

We use a block-dependence structure to allow for weak cross-sectional
dependence in errors.

###### Assumption 2

(a) The errors are such that $\varepsilon=V_{\varepsilon}^{1/2}W\Sigma^{1/2}$,
where $W=[w_{1}:\cdots:w_{n}]$ is a $T\times n$ random matrix of standardized
errors terms $w_{i,t}$ that are independent across $i$ and uncorrelated across
$t$, and $\Sigma=(\sigma_{i,j})$ is a positive-definite symmetric $n\times n$
matrix, such that
$\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum_{i=1}^{n}\sigma_{ii}=1$.
(b) Matrix $\Sigma$ is block diagonal with $J_{n}$ blocks of size
$b_{m,n}=B_{m,n}n$, for $m=1,...,J_{n}$, where $J_{n}\rightarrow\infty$ as
$n\rightarrow\infty$, and $I_{m}$ denotes the set of indices in block $m$. (c)
There exist constants $\delta\in[0,1]$ and $C>0$ such that $\underset{i\in
I_{m}}{\max}\sum_{j\in I_{m}}|\sigma_{i,j}|\leq Cb_{m,n}^{\delta}$. (d) The
block sizes $b_{m,n}$ and block number $J_{n}$ are such that
$n^{2\delta}\sum_{m=1}^{J_{n}}B_{m,n}^{2(1+\delta)}=o(1)$.

The block-dependence structure as in Assumption 2 is satisfied, for instance,
when there are unobserved industry-specific factors independent among
industries and over time, as in Ang, Liu, and Schwarz (2020). In empirical
applications, blocks in $\Sigma$ can match industrial sectors (Fan, Furger,
and Xiu (2016), Gagliardini, Ossola, and Scaillet (2016)). As already
remarked, the diagonal elements of $V_{\varepsilon}$ are the sample
realizations of the common component driving the variance of the error terms
at times $t=1,...,T$; see e.g. Barigozzi and Hallin (2016), Renault, Van Der
Heijden and Werker (2023) for theory and empirical evidence pointing to
variance factors. A sphericity assumption cannot accommodate such a common
time-varying component. Assumption 2 (a) is coherent with Assumption 1.
Indeed,
$\tilde{V}_{\varepsilon}=V_{\varepsilon}^{1/2}\frac{1}{n}\sum_{i,j=1}^{n}\sigma_{i,j}E[w_{i}w_{j}^{\prime}]V_{\varepsilon}^{1/2}=\frac{1}{n}\sum_{i=1}^{n}\sigma_{ii}V_{\varepsilon}$
is diagonal. Hence, $\tilde{V}_{\varepsilon}$ is a scalar multiple of
$V_{\varepsilon}$, and converges to $V_{\varepsilon}$ under the normalization
$\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum_{i=1}^{n}\sigma_{ii}=1$.
That normalization is without loss of generality by rescaling of the
parameters. Assumption 2 (c) builds on Bickel and Levina (2008), and
$\delta<1$ holds under sparsity, vanishing correlations or mixing dependence
within blocks. With blocks of equal size, Assumption 2 (d) holds for
$J_{n}=n^{\bar{\alpha}}$ and $\bar{\alpha}>\frac{2\delta}{2\delta+1}$. Having
$\delta<1$ helps relaxing this condition on block granularity, however it is
not strictly necessary because we allow value $\delta=1$.

In the proof of Proposition 2 below, we establish an asymptotic expansion for
the LR test statistic under the null hypothesis of $k$ latent factors, by
deriving an asymptotic expansion for $n\|\hat{S}\|^{2}$. For a
$(T-k)\times(T-k)$ symmetric matrix $Z=(z_{i,j})$ and
$p=\frac{1}{2}(T-k)(T-k+1)$, let us define the $p$-dimensional vector
$vech(Z)=\left(\frac{1}{\sqrt{2}}z_{11},...,\frac{1}{\sqrt{2}}z_{T-k,T-k},\\{z_{i,j}\\}_{i<j}\right)^{\prime}$,
where the pairs of indices $(i,j)$ with $i<j$ are ranked as
$(1,2),(1,3),...,(1,T-k),(2,3),...,(T-k-1,T-k)$.131313This definition of the
half-vectorization operator for symmetric matrices differs from the usual one
by the ordering of the elements, and the rescaling of the diagonal elements.
It is more convenient for our purposes (see proof of Lemma 10). For instance,
it holds $\frac{1}{2}\|A\|^{2}=vech(A)^{\prime}vech(A)$, for a symmetric
matrix $A$. Moreover, let us define the $p\times T$ matrix
$\boldsymbol{X}=\left[vech(G^{\prime}E_{1,1}G)~{}:~{}\cdots~{}:~{}vech(G^{\prime}E_{T,T}G)\right]$,
where $E_{t,t}$ denotes the $T\times T$ matrix with entry 1 in position
$(t,t)$ and $0$ elsewhere, and $G$ is a $T\times(T-k)$ matrix such that
$F^{\prime}V_{\varepsilon}^{-1}G=0$ and
$G^{\prime}V_{\varepsilon}^{-1}G=I_{T-k}$. Matrix $G$ is unique up to post-
multiplication by an orthogonal matrix, and the columns of $\boldsymbol{X}$
span the linear space $\\{vech(G^{\prime}DG):D\text{ diagonal}\\}$. Then, the
asymptotic expansion of $n\|\hat{S}\|^{2}$ under $H_{0}(k)$ is:

$\displaystyle\frac{n}{2}\|\hat{S}\|^{2}=vech(Z_{n}^{*})^{\prime}M_{\boldsymbol{X}}vech(Z_{n}^{*})+o_{p}(1),$
(5)

where $Z^{*}_{n}:=G^{\prime}V_{\varepsilon}^{-1}Z_{n}V_{\varepsilon}^{-1}G$,
$Z_{n}:=\sqrt{n}\left(\frac{1}{n}\varepsilon\varepsilon^{\prime}-\tilde{V}_{\varepsilon}\right)$,
and matrix
$M_{\boldsymbol{X}}:=I_{p}-\boldsymbol{X}(\boldsymbol{X}^{\prime}\boldsymbol{X})^{-1}\boldsymbol{X}^{\prime}$
is idempotent of rank $p-T=df$.141414We can extend results like (5) to test
statistics that are generic functions of the eigenvalues of matrix $\hat{S}$
by using the Weyl’s inequalities (see e.g. Bernstein (2009)), and develop test
statistics along the lines of FGS. The full-rank condition for matrix
$\boldsymbol{X}$ corresponds to the local identification condition in
Assumption A.5 ((see Lemma 6)), analogously as in linear regression. The
expansion in (5) does not depend on the diagonal elements of $Z_{n}$, since
$vech(G^{\prime}V_{\varepsilon}^{-1}diag(Z_{n})V_{\varepsilon}^{-1}G)$ is
spanned by the columns of $\boldsymbol{X}$, and thus is annihilated by the
projection matrix $M_{\boldsymbol{X}}$. Such an irrelevance of the diagonal
elements of $Z_{n}$ is an implication of Proposition 1 (c).

We now establish the distributional convergence $Z^{*}_{n}\Rightarrow Z^{*}$
as $n\rightarrow\infty$ and $T$ is fixed, where $Z^{*}$ is a Gaussian
symmetric matrix variate. By the block structure in Assumption 2 (b), we can
write $Z^{*}_{n}$ as a sum of independent zero-mean terms:
$Z^{*}_{n}=\frac{1}{\sqrt{n}}\sum_{m=1}^{J_{n}}z^{*}_{m,n},$ where the
variables in the triangular array $z^{*}_{m,n}=\sum_{i\in
I_{m}}G^{\prime}V_{\varepsilon}^{-1}\left(\varepsilon_{i}\varepsilon_{i}^{\prime}-E[\varepsilon_{i}\varepsilon_{i}^{\prime}]\right)V_{\varepsilon}^{-1}G$
are independent across $m$ and such that $E[z^{*}_{m,n}]=0$. In Appendix B, we
invoke the CLT for independent heterogeneous variables to
$vech(Z^{*}_{n})=\frac{1}{\sqrt{n}}\sum_{m=1}^{J_{n}}vech(z^{*}_{m,n})$ and
use Assumptions 2 (c) and (d) to check the Liapunov condition. We get
$Z^{*}_{n}\Rightarrow Z^{*}$, where $vech(Z^{*})\sim N(0,\Omega_{Z^{*}})$ and
$\Omega_{Z^{*}}=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum_{m=1}^{J_{n}}V[vech(z^{*}_{m,n})]$.
Then, the asymptotic expansion in (5) yields the asymptotic distribution for
the LR test statistic defined in (4) under the null hypothesis, and is given
in Proposition 2 (a) below.

In Proposition 2 (b) hereafter, we use
$\hat{G}=\hat{V}_{\varepsilon}^{1/2}\hat{Q}$, where $\hat{Q}$ is the
$T\times(T-k)$ matrix of standardized eigenvectors of $\hat{S}$ corresponding
to its $T-k$ non-zero eigenvalues ranked in decreasing order (see Proposition
1 (a)),151515 $\hat{Q}$ is the $T\times(T-k)$ matrix of standardized
eigenvectors of
$\hat{V}^{-1/2}_{\varepsilon}\hat{V}_{y}\hat{V}^{-1/2}_{\varepsilon}$
corresponding to its $T-k$ smallest eigenvalues ranked in decreasing order.
and $\hat{\boldsymbol{X}}$ the matrix obtained by replacing $G$ with $\hat{G}$
in $\boldsymbol{X}$.

###### Proposition 2

Let Assumptions 1-2 and A.1-A.7 hold. As $n\rightarrow\infty$ and $T$ is
fixed, under the null hypothesis $H_{0}(k)$ of $k$ latent factors, (a)
$LR(k)\Rightarrow\sum_{j=1}^{df}\mu_{j}\chi^{2}_{j}(1)$, where the
$\chi^{2}_{j}(1)$ are independent chi-square variables with one degree of
freedom, and the $\mu_{j}$ are the $df$ non-zero eigenvalues of matrix
$M_{\boldsymbol{X}}\Omega_{Z^{*}}M_{\boldsymbol{X}}$. (b)
$\hat{\mu}_{j}\overset{p}{\rightarrow}\mu_{j}$, where $\hat{\mu}_{j}$,
$j=1,...,df$, are the non-zero eigenvalues of matrix
$M_{\hat{\boldsymbol{X}}}\hat{\Omega}_{Z^{*}}M_{\hat{\boldsymbol{X}}}$, and
$\hat{\Omega}_{Z^{*}}=\frac{1}{n}\sum_{m=1}^{J_{n}}vech(\hat{z}^{*}_{m,n})vech(\hat{z}^{*}_{m,n})^{\prime}$
with $\hat{z}^{*}_{m,n}=\sum_{i\in
I_{m}}\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}\hat{\varepsilon}_{i}\hat{\varepsilon}_{i}^{\prime}\hat{V}_{\varepsilon}^{-1}\hat{G}$
and $\hat{\varepsilon}_{i}=M_{\hat{F},\hat{V}_{\varepsilon}}(y_{i}-\bar{y})$.
Under the alternative hypothesis $H_{1}(k)$ of more than $k$ latent factors,
(c) $LR(k)\geq Cn$, w.p.a. $1$ for a constant $C>0$, and
$\hat{\mu}_{j}=O_{p}(n\sum_{m=1}^{J_{n}}B_{m,n}^{2})=o_{p}(n)$.

Proposition 2 (a) shows that we depart from classical FA theory since we have
convergence to a weighted sum of chi-square variates instead of a single chi-
square variate (see the discussion in the next section). In Proposition 2 (b),
matrix $\hat{G}$ is consistent for $G$ up to a rotation. The eigenvalues
$\hat{\mu}_{j}$ are unaffected by such rotation because eigenvalues are
invariant under pre and post multiplication by an orthogonal matrix and its
transpose.161616In Proposition 9 in Appendix E.6, we study how
$\boldsymbol{X}$ and $M_{\boldsymbol{X}}\Omega_{Z^{*}}M_{\boldsymbol{X}}$ are
transformed under different choices for the rotation of $G$. The eigenvalues
$\mu_{j}$ are invariant to such rotation as expected. With fixed $T$, the GLS
residuals $\hat{\varepsilon}_{i}$ are asymptotically close to
$M_{F,V_{\varepsilon}}\varepsilon_{i}$ and not to the true errors
$\varepsilon_{i}$. However, it does not impede the consistency of the
eigenvalues $\hat{\mu}_{j}$, since
$G^{\prime}V_{\varepsilon}^{-1}M_{F,V_{\varepsilon}}=G^{\prime}V_{\varepsilon}^{-1}$.
When we apply the CLT, centering is done implicitly since
$M_{\boldsymbol{X}}vech(G^{\prime}V_{\varepsilon}^{-1}\varepsilon_{i}\varepsilon_{i}^{\prime}V_{\varepsilon}^{-1}G)=M_{\boldsymbol{X}}vech(G^{\prime}V_{\varepsilon}^{-1}(\varepsilon_{i}\varepsilon_{i}^{\prime}-E[\varepsilon_{i}\varepsilon_{i}^{\prime}])V_{\varepsilon}^{-1}G)$.
Indeed
$vech(G^{\prime}V_{\varepsilon}^{-1}E[\varepsilon_{i}\varepsilon_{i}^{\prime}]V_{\varepsilon}^{-1}G)$
is spanned by the columns of $\boldsymbol{X}$ since
$E[\varepsilon_{i}\varepsilon_{i}^{\prime}]=\sigma_{ii}V_{\varepsilon}$ is
diagonal. We can consistently estimate the critical values of the statistics
by simulating a large number of draws from
$\sum_{j=1}^{df}\hat{\mu}_{j}\chi^{2}_{j}(1)$. Hence, even if FA residuals are
not consistent estimates of the true errors in short panels, we are still able
to supply a feasible asymptotic distributional theory for our empirical
applications under a block-dependence structure. Proposition 2 (c) gives test
consistency against global alternatives.

##  4 Impact of assuming Gaussian errors and sphericity

In this section, we particularize the general distributional results of
Proposition 2 to three important cases, namely Gaussian errors, settings where
the asymptotic distribution under Gaussian errors still holds for the LR test
statistics (up to scaling), and spherical errors. We further discuss the
impact of wrongly assuming Gaussian errors or sphericity when conducting
inference on the number of latent factors.

###  4.1 Gaussian errors

Let us consider the case where the errors
$\varepsilon_{i}\overset{ind}{\sim}N(0,\sigma_{ii}V_{\varepsilon})$ are
independent Gaussian vectors. From classical FA theory, we expect that the
statistic $LR(k)$ admits asymptotically a chi-square distribution with $df$
degrees of freedom in the cross-sectionally homoskedastic case, i.e.,
$\sigma_{ii}=1$ for all assets $i$. We cannot expect that this distributional
result applies to the Gaussian framework in full generality, since - even in
such a case - our setting corresponds to a pseudo model (because the
$\sigma_{ii}$ may be heterogeneous across $i$, and the $\beta_{i}$ are treated
as fixed effects, namely incidental parameters, instead of Gaussian random
effects). Under the normality assumption for the error terms, we have
${\varepsilon}^{*}_{i}:=G^{\prime}V_{\varepsilon}^{-1}\varepsilon_{i}\overset{ind}{\sim}N(0,\sigma_{ii}I_{T-k})$.
Thus, by the Liapunov CLT, the distributional limit of
$\frac{1}{\sqrt{q}}Z_{n}^{*}=\sqrt{n/q}\left(\frac{1}{n}{\varepsilon^{*}}(\varepsilon^{*})^{\prime}-\frac{1}{n}E\left[{\varepsilon^{*}}(\varepsilon^{*})^{\prime}\right]\right)$
is in the Gaussian Orthogonal Ensemble (GOE) for dimension $T-k$ (see e.g. Tao
(2012)), i.e., $\frac{1}{\sqrt{q}}vech(Z^{*})\sim N(0,I_{p})$, where
$q:=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum_{i=1}^{n}\sigma_{ii}^{2}$.
Then, from (5) we get $LR(k)\Rightarrow q\chi^{2}(df)$, i.e., we get
convergence to a scaled chi-square variate $q\chi^{2}(df)$. In the cross-
sectionally homoskedastic case, we have $q=1$ yielding the classical
$\chi^{2}(df)$ result. On the contrary, cross-sectional heterogeneity in the
unconditional idiosyncratic variances yields $q>1$ and a deviation from
classical FA theory even in the Gaussian case. In the Gaussian case,
unobserved heterogeneity across asset idiosyncratic variances leads to an
oversized LR test if we use critical values from the chi-square table without
proper scaling. In our Monte Carlo experiments with non-Gaussian errors, we
get size distortion by around 80 percentage points when $T=24$ for a nominal
size of 5%, and thus an average estimated number of factors over 3 instead of
2.

###  4.2 Validity of the scaled asymptotic chi-square test

In this subsection, we investigate sufficient conditions for the validity of
the convergence of the LR statistic to a scaled chi-square variate
$\chi^{2}(df)$, but in special cases beyond Gaussianity of errors. For this
purpose, let us recall that the asymptotic expansion of $LR(k)$ in (5) only
involves the out-of-diagonal elements of $Z_{n}$. Under independent Gaussian
errors (Section 4.1), by the Liapunov CLT, we have $Z_{n}\Rightarrow Z$, where
$Z_{t,s}\sim N(0,qV_{\varepsilon,tt}V_{\varepsilon,ss})$, for $t>s$, mutually
independent, where the $V_{\varepsilon,tt}$ are the diagonal elements of
matrix $V_{\varepsilon}$. We deduce that any setting featuring the same joint
asymptotic distribution for the out-of-diagonal elements of random matrix
$Z_{n}$ leads to an asymptotic distribution for the LR statistic similar to
the Gaussian case.

###### Proposition 3

Let Assumptions 1-2, A.1-A.6 hold with (a)
$\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum_{i=1}^{n}E[\varepsilon_{i,t}\varepsilon_{i,s}\varepsilon_{i,r}\varepsilon_{i,p}]=qV_{\varepsilon,tt}V_{\varepsilon,ss}$,
when $t=r>s=p$, for a constant $q>0$, and $=0$ in all other cases with $t>s$
and $r>p$, and (b) let
$\kappa=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum_{m=1}^{J_{n}}\sum_{i\neq
j\in I_{m}}\sigma_{ij}^{2}$ as in Assumption A.3 (b). Then,
$LR(k)\Rightarrow\bar{q}\chi^{2}(df)$ under $H_{0}(k)$ for
$\bar{q}:=q+\kappa$.

Conditions (a) and (b) in Proposition 3 generalize the correctness of the
scaled chi-square test beyond Gaussianity and error independence across time
and assets. Under Assumption 2, Condition (a) is satisfied if the standardized
error terms $w_{i,t}$ are conditionally homoskedastic martingale difference
sequences. However, Condition (a) excludes empirically relevant cases such as
ARCH processes for $w_{i,t}$, because, in that case,
$\frac{1}{V_{\varepsilon,tt}V_{\varepsilon,ss}}E[\varepsilon_{i,t}^{2}\varepsilon_{i,s}^{2}]$
depends on lag $t-s$. Hence, serial correlation in squared idiosyncratic
errors is responsible for the deviation of the LR test from the scaled chi-
square asymptotic distribution. This setting is covered by the general results
in Proposition 2. Anderson and Amemiya (1988) establish the asymptotic
distribution of FA estimates assuming that the error terms are i.i.d. across
sample units and deploy an assumption that is analogue to Condition (a) above
in their Corollary 2. The i.i.d. assumption in our case implies
$\sigma_{ii}=1$ for all $i$, which results in a cross-sectionally
homoskedastic setting.171717If the $\sigma_{ii}$ are treated as i.i.d. random
effects independent of the errors, and we exclude cross-sectional correlation
of errors to simplify, we recover the i.i.d. condition of the data. However,
the random $\sigma_{ii}$ yields a stochastic common factor across time that
breaks the condition in Corollary 2 of Anderson and Amemiya (1988). That
setting is irrealistic in our application, as it would imply that the
idiosyncratic variance is the same for all assets. Our results show that
establishing the asymptotic distribution of the test statistics, especially
the AUMPI property of LR test (see Section 5), in a general setting with non-
Gaussian errors, heterogeneous idiosyncratic variances and ARCH effects, is
challenging, but still possible.

###  4.3 Spherical errors

If errors are spherical, i.e., matrix $V_{\varepsilon}=\bar{\sigma}^{2}I_{T}$
is a multiple of the identity with unknown parameter $\bar{\sigma}^{2}>0$,
then the asymptotic distribution of $LR(k)$ corresponds to a special case of
Proposition 2 (a). If sphericity is imposed in the estimation procedure, i.e.,
$\hat{V}_{\varepsilon}$ becomes
$\hat{V}_{\varepsilon,c}=\hat{\sigma}^{2}I_{T}$, the constrained FA estimator
boils down to the Principal Component Analysis (PCA) estimator; see Anderson
and Rubin (1956) Section 7.3. Let $\hat{F}_{c}$ denote the FA estimator under
sphericity constraint. Then, $\hat{F}_{c}$ is the matrix of eigenvectors of
matrix $\hat{V}_{y}$ standardized such that
$\hat{F}_{c}^{\prime}\hat{F}_{c}=diag(\hat{\delta}_{1}-\hat{\sigma}^{2},...,\hat{\delta}_{k}-\hat{\sigma}^{2})$,
and $\hat{\sigma}^{2}=\frac{1}{T-k}\sum_{j=k+1}^{T}\hat{\delta}_{j}$, where
$\hat{\delta}_{j}$ denotes the $j$th largest eigenvalue of matrix
$\hat{V}_{y}$. The matrix $\hat{S}$ under sphericity constraint becomes
$\hat{S}_{c}=\frac{1}{\hat{\sigma}^{2}}M_{\hat{F}_{c}}(\hat{V}_{y}-\hat{\sigma}^{2}I_{T})M_{\hat{F}_{c}}=\frac{1}{\hat{\sigma}^{2}}\left(\frac{1}{n}\hat{\varepsilon}_{c}\hat{\varepsilon}_{c}^{\prime}\right)-M_{\hat{F}_{c}},$
where
$M_{\hat{F}_{c}}=I_{T}-\hat{F}_{c}(\hat{F}_{c}^{\prime}\hat{F}_{c})^{-1}\hat{F}^{\prime}$
and $\hat{\varepsilon}_{c}=M_{\hat{F}_{c}}\tilde{Y}$ is the matrix of OLS
residuals. Then, the constrained LR statistic becomes
$LR_{c}(k):=-n\sum_{j=k+1}^{T}\log(1+\frac{\hat{\delta}_{j}-\hat{\sigma}^{2}}{\hat{\sigma}^{2}})$
and reduces to the LR statistic invoqued by Onatski (2023) in his discussion
of FGS.181818We have $\hat{\gamma}_{j}=\hat{\delta}_{j}/\hat{\sigma}^{2}-1$,
and
$LR_{c}(k)=\frac{n}{2\hat{\sigma}^{4}}\sum_{j=k+1}^{T}(\hat{\delta}_{j}^{2}-\hat{\sigma}^{2})^{2}+o_{p}(1)$,
i.e., the constrained LR statistic is asymptotically equivalent to the sum of
squared deviations of the $T-k$ smallest eigenvalues from their mean. Besides,
by similar results, we have that the eigenvalue spacing statistic
$\sqrt{n}\mathscr{S}(k):=\hat{\gamma}_{k+1}-\hat{\gamma}_{T}$ corresponds to
the statistic considered in FGS divided by $\hat{\sigma}^{2}$, and its
asymptotic distribution under sphericity coincides with that obtained by FGS.
Under spherical errors, Onatski (2023) shows
$LR_{c}(k)\Rightarrow\frac{1}{2}(Tr[(Z^{*})^{2}]-\frac{1}{T-k}[Tr(Z^{*})]^{2})=vech(Z^{*})^{\prime}M_{\boldsymbol{x}}vech(Z^{*})$,
where
$M_{\boldsymbol{x}}=I_{p}-\boldsymbol{x}(\boldsymbol{x}^{\prime}\boldsymbol{x})^{-1}\boldsymbol{x}^{\prime}$
with $\boldsymbol{x}=vech(I_{T-k})$,
$Z^{*}=\frac{1}{\bar{\sigma}^{4}}G^{\prime}ZG$, $Z$ is the distributional
limit of $Z_{n}$, and $G$ is a $T\times(T-k)$ matrix such that $F^{\prime}G=0$
and $G^{\prime}G=\bar{\sigma}^{2}I_{T-k}$. Hence the asymptotic distribution
of the constrained LR statistic under sphericity is
$\sum_{j=1}^{p-1}\mu_{j}\chi^{2}(1)$, where the $\chi^{2}(1)$ are independent
and the $\mu_{j}$ are the non-zero eigenvalues of matrix
$M_{\boldsymbol{x}}\Omega_{Z^{*}}M_{\boldsymbol{x}}$, with
$\Omega_{Z^{*}}=V[vech(Z^{*})]$. Under Gaussian errors, it simplifies to
$LR_{c}(k)\Rightarrow q\chi^{2}(p-1)$ (see Section 4.1). In a PCA setting when
sphericity fails to hold, a massive over-rejection by over 80 percentage
points for a nominal size of 5% is again observed in our Monte Carlo
experiments. It is true both for the constrained LR test and for the
eigenvalue spacing test of FGS. We get an average estimated number of factors
often above 10 instead of 2 when $T=24$.

##  5 Local asymptotic power

In this section, we study the asymptotic power of the test statistics against
local alternatives in which we have $k$ (strong) factors plus a weak factor.
Specifically, under $H_{1,loc}(k)$, we have $\sqrt{n}\gamma_{k+1}\rightarrow
c_{k+1}$ as $n\rightarrow\infty$, with $c_{k+1}>0$. The (drifting) DGP is
$Y=\mu
1_{n}^{\prime}+F\beta^{\prime}+F_{k+1}\beta_{loc}^{\prime}+\varepsilon$, where
$\beta_{loc}$ is the loading vector for the $(k+1)$th factor, and the factor
vector is normalized such that $F_{k+1}=\sqrt{\gamma_{k+1}}\rho_{k+1}$ with
$\rho_{k+1}^{\prime}V_{\varepsilon}^{-1}\rho_{k+1}=1$ and
$F^{\prime}V_{\varepsilon}^{-1}\rho_{k+1}=0$. Thus, we can write
$\rho_{k+1}=G\xi_{k+1}$ for a $T-k$ dimensional vector $\xi_{k+1}$ with unit
norm. Scalar $c_{k+1}$ and vector $\xi_{k+1}$ yield the (normalized) strength
and the direction of the local alternative.

###  5.1 Asymptotic distributions under local alternatives

We derive an asymptotic expansion of $n\|\hat{S}\|^{2}$ under $H_{1,loc}(k)$
using similar arguments as in the proof of Proposition 2 (a) (see the proof of
Proposition 4 in Appendix B for the derivation):

$\displaystyle\frac{n}{2}\|\hat{S}\|^{2}=vech(Z_{n,loc}^{*})^{\prime}M_{\boldsymbol{X}}vech(Z_{n,loc}^{*})+o_{p}(1),$
(6)

where $Z_{n,loc}^{*}=Z^{*}_{n}+c_{k+1}\xi_{k+1}\xi_{k+1}^{\prime}$. From the
CLT, we have $Z_{n,loc}^{*}\Rightarrow Z_{loc}^{*}$ where
$Z_{loc}^{*}=Z^{*}+c_{k+1}\xi_{k+1}\xi_{k+1}^{\prime}$. Matrix variate
$Z^{*}_{loc}$ is a non-central symmetric Gaussian matrix. The non-zero mean
depends in general on both $c_{k+1}$ and $\xi_{k+1}$, while the variances and
covariances of the elements of $Z^{*}_{loc}$ are the same as those of $Z^{*}$.
Then, we deduce from (6) that the asymptotic distribution of the $LR(k)$
statistic under the local alternative hypothesis is a weighted average of $df$
mutually independent non-central chi-square variables, with non-centrality
parameters depending on
$vech(\Delta):=M_{\boldsymbol{X}}vech(c_{k+1}\xi_{k+1}\xi^{\prime}_{k+1})$.

###### Proposition 4

Let Assumptions 1-2, A.1-A.6 hold. Under the local alternative hypothesis
$H_{1,loc}(k)$, we have as $n\rightarrow\infty$ and $T$ is fixed,
$LR(k)\Rightarrow\sum_{j=1}^{df}\mu_{j}\chi^{2}(1,\lambda_{j}^{2})$, where
$\lambda_{j}^{2}=\mu_{j}^{-1}[v_{j}^{\prime}vech(\Delta)]^{2}$, and the
$\mu_{j}$ and $v_{j}$ are the non-zero eigenvalues and the associated
standardized eigenvectors of matrix
$M_{\boldsymbol{X}}\Omega_{Z^{*}}M_{\boldsymbol{X}}$.

The non-centrality term $vech(\Delta)$ is in charge of the asymptotic local
power of the statistics. When this vector is null, the asymptotic local power
is zero. Indeed, for some local alternatives the $(k+1)$th weak factor can be
absorbed in the diagonal variance matrix $V_{\varepsilon}$ of the error terms.
More precisely, in Appendix E.3 ii), we show that
$V_{y}+\frac{c_{k+1}}{\sqrt{n}}\rho_{k+1}\rho_{k+1}^{\prime}=F^{*}(F^{*})^{\prime}+V_{\varepsilon}^{*}+\frac{1}{\sqrt{n}}G\Delta
G^{\prime}+o(1/\sqrt{n})$ for some $T\times k$ matrix $F^{*}$ and diagonal
matrix $V_{\varepsilon}^{*}$, which yields asymptotically a $k$-factor model
when $\Delta=0$. We have
$\lambda^{2}:=\sum_{j=1}^{df}\mu_{j}\lambda_{j}^{2}=\frac{1}{2}\|\Delta\|^{2}$,
i.e., the half squared Frobenius norm of the matrix measuring local distance
from the $k$-factor specification. It follows that the asymptotic local power
of the LR statistic is non null as long as $\lambda^{2}>0$, i.e., it has non-
trivial asymptotic power against any proper local alternative. In our Monte
Carlo experiments in Appendix C, we find that the LR statistic has size close
to the nominal value, and power against global as well as local alternatives
with time dimension as small as $T=6$.

Under the normality of errors, or more generally the conditions of Proposition
3, using that matrix $\frac{1}{\sqrt{\bar{q}}}Z^{*}$ is in the GOE for
dimension $T-k$, i.e., $vech(Z^{*})\sim N(0,\bar{q}I_{p})$, we have
$LR(k)\Rightarrow\bar{q}\chi^{2}(df,\lambda^{2}/\bar{q})$ from (6). The local
power is a function solely of the squared Euclidean norm of the vector
$vech(\Delta)$ measuring local distance from the $k$-factor specification,
divided by $\bar{q}$.

###  5.2 AUMPI tests

In this subsection, we investigate asymptotic local optimality of the LR
statistic for testing hypotheses on the number of latent factors. In our
framework with composite null and alternative hypotheses and multi-dimensional
parameter, we cannot expect in general to establish Uniformly Most Powerful
(UMP) tests. Instead, we can establish an optimality property by restricting
the class of tests to invariant tests (e.g. Lehmann and Romano (2005)). We
focus on statistics with test functions $\phi$ written on the elements of
matrix $\hat{S}$. To eliminate the asymptotic redundancy in the elements of
$\hat{S}$, we actually consider the test class
$\mathscr{C}=\left\\{\phi~{}:~{}\phi=\phi(\hat{W})\right\\}$ with
$\hat{W}:=\sqrt{n}\hat{\boldsymbol{D}}^{\prime}vech(\hat{S}^{*})$, where
$\hat{S}^{*}=\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1/2}\hat{S}\hat{V}_{\varepsilon}^{-1/2}\hat{G}$
and $\hat{\boldsymbol{D}}$ is a $p\times df$ full-rank matrix, such that
$M_{\hat{\boldsymbol{X}}}=\hat{\boldsymbol{D}}\hat{\boldsymbol{D}}^{\prime}$
and $\hat{\boldsymbol{D}}^{\prime}\hat{\boldsymbol{D}}=I_{df}$. We have that
$\hat{S}^{*}=diag(\hat{\gamma}_{k+1},\dots,\hat{\gamma}_{T})$, since $\hat{S}$
has the eigendecomposition
$\hat{V}_{\varepsilon}^{-1/2}\hat{G}diag(\hat{\gamma}_{k+1},\dots,\hat{\gamma}_{T})\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1/2}$
from Proposition 1 (a). The diagonal matrix $\hat{S}^{*}$ contains the
information in $\hat{S}$ beyond orthogonality of its rows and columns to
$\hat{V}_{\varepsilon}^{-1/2}\hat{F}$. The vector $\hat{W}$ contains the
information in $\sqrt{n}vech(\hat{S}^{*})$ beyond orthogonality to
$\hat{\boldsymbol{X}}$.191919From Proposition 1 (c), we have
$0=diag(\hat{S})=2\hat{V}^{-1}_{\varepsilon}\hat{\boldsymbol{X}}^{\prime}vech(\hat{S}^{*})$.
Therefore, $vech(\hat{S}^{*})$ lies in the orthogonal complement of the range
of $\hat{\boldsymbol{X}}$ in sample.

Matrices $\hat{G}$ and $\hat{\boldsymbol{D}}$ are both consistent up to post-
multiplication by an orthogonal matrix. This point yields a group of
orthogonal transformations under which we require the test statistics to be
invariant.202020Here, we do not deal with invariance to data transformations
but rather with invariance to parameterization of $\hat{G}$ and
$\hat{\boldsymbol{D}}$. However, if we consider tests based on the elements of
vector $\hat{W}$, this difference is immaterial. In Appendix E.6, we show that
the maximal invariant under this group is provided by
$\hat{W}^{\prime}\hat{W}=nvech(\hat{S}^{*})^{\prime}M_{\hat{\boldsymbol{X}}}vech(\hat{S}^{*})$.
Since $\sqrt{n}vech(\hat{S}^{*})$ belongs to the range of matrix
$M_{\hat{\boldsymbol{X}}}$, we have
$\hat{W}^{\prime}\hat{W}=\frac{n}{2}\|\hat{S}^{*}\|^{2}=\frac{n}{2}\|\hat{S}\|^{2}$.
Therefore, the invariant tests are functions of the squared norm of $\hat{S}$,
which is asymptotically equivalent to the LR statistic (up to the factor
$1/2$).

In the Gaussian case, or more generally under the conditions of Proposition 3,
the LR statistic follows asymptotically a scaled non-central chi-square
distribution with $df$ degrees of freedom and non-centrality parameter
$\lambda^{2}/\bar{q}=\sum_{j=1}^{df}\lambda_{j}^{2}$ as shown in the previous
subsection. Thus, we can simplify the null and alternative hypotheses of our
testing problem asymptotically and locally to a one-sided test with null
hypothesis $H_{0}(k):\lambda^{2}=0$ vs. alternative hypothesis
$H_{1,loc}(k):\lambda^{2}>0$. The scaling constant $q>0$ plays no role in the
power analysis. It means that the LR test is an AUMPI test (Lehmann and Romano
(2005) Chapters 3 and 13). Indeed, the density $g(z;df,\lambda^{2})$ of the
$\chi^{2}(df,\lambda^{2})$ distribution is Totally Positive of order 2 ($TP2$)
in $z$ and $\lambda^{2}$ (Eaton (1987) Example A.1 p. 468); see Miravete
(2011) for a review of applications of TP2 in economics. A density, which is
$TP2$ in $z$ and $\lambda^{2}$, has the Monotone Likelihood Ratio (MLR)
property (Eaton (1987) p. 467). Since $g(z;df,\lambda^{2})/g(z;df,0)$ is an
increasing function in $z$, it gives the AUMPI property.

In the general case with $df>1$, when neither Gaussianity nor the conditions
of Proposition 3 apply, we cannot use the same reasoning, since the density
$f(z;\lambda_{1},...,\lambda_{df})$ of
$\sum_{j=1}^{df}\mu_{j}\chi^{2}(1,\lambda_{j}^{2})$, with $\mu_{j}>0$,
$j=1,...,df$, is not a function of
$\lambda^{2}=\sum_{j=1}^{df}\mu_{j}\lambda_{j}^{2}$ only, and thus cannot be
$TP2$ in $z$ and $\lambda^{2}$. Instead, our strategy to get a new AUMPI
result uses a power series representation of the density of
$\sum_{j=1}^{df}\mu_{j}\chi^{2}(1,\lambda_{j}^{2})$ in terms of central chi-
square densities from Kotz, Johnson, and Boyd (1967). Under the sufficient
condition (7) in Proposition 5, the density ratio
$\frac{f(z;\lambda_{1},...,\lambda_{df})}{f(z;0,...,0)}$ is monotone
increasing in $z$.

###### Proposition 5

Let Assumptions 1-2, A.1-A.6 hold. (a) Let us assume that, for any DGP in the
subset $\bar{H}_{1,loc}(k)\subset H_{1,loc}(k)$ of the local alternative
hypothesis, we have for any integer $m\geq 3$:

$\displaystyle\sum_{j>l\geq
0,j+l=m}\frac{(j-l)\Gamma(\frac{df}{2})^{2}}{\Gamma(\frac{df}{2}+j)\Gamma(\frac{df}{2}+l)}[c_{j}(\lambda_{1},...,\lambda_{df})c_{l}(0,...,0)-c_{l}(\lambda_{1},...,\lambda_{df})c_{j}(0,...,0)]\geq
0,\quad$ (7)

where $\Gamma(\cdot)$ is the Gamma function,
$c_{j}(\lambda_{1},...,\lambda_{df}):=E[Q(\lambda_{1},...,\lambda_{df})^{j}]/j!$
for
$Q(\lambda_{1},...,\lambda_{df})=\frac{1}{2}\sum_{j=1}^{df}(\sqrt{\nu_{j}}X_{j}+\sqrt{1-\nu_{j}}\lambda_{j})^{2}$,
$\nu_{j}=1-\frac{1}{\mu_{j}}\mu_{1}$ with the $\mu_{j}$ ranked in increasing
order, and $X_{j}\sim N(0,1)$ are mutually independent. Then, the statistic
$LR(k)$ yields an AUMPI test against $\bar{H}_{1,loc}(k)$. (b) Suppose that
either $\lambda_{1}^{2}+(1-\nu_{2})\lambda_{2}^{2}\geq\nu_{2}$ and
$(1-\nu_{2})\lambda_{2}^{2}\geq\frac{1}{2}\nu_{2}$ when $df=2$, or

$1\\{i=0\\}\lambda_{1}^{2}+\sum_{j=2}^{df-1}\rho_{j}^{i}(1-\nu_{j})\lambda_{j}^{2}+(1-\nu_{df})\lambda_{df}^{2}\geq\frac{\nu_{df}}{i+1}\left(df-2-\sum_{j=2}^{df-1}\rho_{j}^{i+1}\right),$
(8)

for all $i\geq 0$, where $\rho_{j}:=\frac{\nu_{j}}{\nu_{df}}$, when $df\geq
3$. Then, Inequalities (7) hold for any $m\geq 3$.

Conditions (7) involve polynomial inequalities in the parameters $\lambda_{j}$
of the alternative hypothesis, and parameters $\nu_{j}$ of the weights of the
non-central chi-square distributions, $j=1,...,df$. It is challenging to
establish an explicit characterization of the $\lambda_{j}$ and $\nu_{j}$
equivalent to Inequalities (7), unless $df=1$.212121Inequalities (7) with
$df=1$ are easily proved to hold. In such a case, we can use the asymptotic
distribution of a scaled chi-square variable and its MLR property. By
deploying a novel characterization of the
$c_{j}(\lambda_{1},...,\lambda_{df})$ in terms of a recurrence relation (Lemma
2), we establish explicit sufficient conditions in part b) of Proposition 5.
Inequalities (8) are linear in the $\lambda_{j}^{2}$, and define a non-empty
convex domain in the $(\lambda_{1}^{2},..,\lambda_{df}^{2})$ space, that does
not contain the origin $\lambda_{1}=...=\lambda_{df}=0$ (unless the DGP is
such that $\nu_{2}=...=\nu_{df}$, in which case the RHS of (8) is nil for all
$i$ and thus any $\lambda_{j}^{2}$ meet the inequalities). Proposition 5 (b)
implies that, for a given set of values of $df$, the MLR property holds if
$\lambda_{j}\geq\underline{\lambda}$ for all $j$, uniformly for
$\nu_{j}\leq\bar{\nu}$, where $\underline{\lambda}>0$ is a constant that
depends on $\bar{\nu}<1$. Vanishing values of the $\nu_{j}$ correspond to
homogenous weights $\mu_{j}$, i.e., the scaled non-central chi-square
distribution with $df$ degrees of freedom. Hence, the AUMPI property in
Proposition 5 holds in neighborhoods of DGPs that match the conditions of
Proposition 3 (e.g. Gaussian errors) for alternative hypotheses that are
sufficiently separated from the null hypothesis. Besides, Proposition 5 shows
that the Gaussian case is not the only design delivering an AUMPI test.
Further, in the SMC, we establish a new analytical representation of the
coefficients $c_{k}(\lambda_{1},...,\lambda_{df})$ in terms of matrix product
iterations. That analytical representation allows us to check numerically the
validity of Inequalities (7) for given $df$, $\lambda_{j}$, $\nu_{j}$, and
$m=1,....,M$, for a large bound $M$ (see Appendix F). In Appendix F, when
Inequalities (8) are met, we always conclude to the MLR property in the
numerical checks as predicted by the theory of Proposition 5. There, we also
provide numerical evidence that the domain of validity of the MLR property is
relevant for our empirical application. The sufficient conditions (7) and (8)
in Proposition 5 yielding the monotone property of density ratios have
potentially broad application outside the current setting to show AUMPI
properties since other test statistics share an asymptotic distribution
characterized by a positive definite quadratic form in normal vectors.

##  6 Empirical application

In this section, we test hypotheses about the number of latent factors driving
stock returns in short subperiods of the Center for Research in Securities
Prices (CRSP) panel. Then, we decompose the cross-sectional variance into
systematic and idiosyncratic components. We also check whether there is
spanning between the estimated latent factors and standard observed factors.

###  6.1 Testing for the number of latent factors

We consider monthly returns of U.S. common stocks trading on the NYSE, AMEX or
NASDAQ between January 1963 and December 2021, and having a non-missing
Standard Industrial Classification (SIC) code. We partition subperiods into
bull and bear market phases according to the classification methodology of
Lunde and Timmermann (2004).222222We fix their parameter values
$\lambda_{1}=\lambda_{2}=0.2$ for the classification based on the nominal
S&P500 index. Bear periods are close to NBER recessions. We implement the
tests using a rolling window of $T=20$ months, moving forward 12 months each
time (adjacent windows overlap by 8 months), thereby ensuring that we can test
up to $14$ latent factors in each subperiod. The size of the cross-section $n$
ranges from $1768$ to $6142$, and the median is $3680$. We only consider
stocks with available returns over the whole subperiod, so that our panels are
balanced. In each subperiod, we sequentially test $H_{0}(k)$ vs $H_{1}(k)$,
for $k=0,\dots,k_{max}$, where $k_{max}=14$ is the largest nonnegative integer
such that $df>0$ (see Table 3 in OA). We compute the variance-covariance
estimator $\hat{\Omega}_{\bar{Z}^{*}}$ using a block structure implied by the
partitioning of stocks by the first two digits of their SIC code. The number
of blocks ranges from 61 to 87 over the sample, and the number of stocks per
block ranges from 1 to 641. The median number of blocks is 76 and the median
number of stocks per block is 21. We display the p-values of the statistic
$LR(k)$ over time for each subperiod in the upper panel of Figure 1, stopping
at the smallest $k$ such that $H_{0}(k)$ is not rejected at level
$\alpha_{n}=10/n_{max}$, where $n_{max}$ is the largest cross-sectional sample
size over all subperiods, so that $\alpha_{n}=0.16\%$ in our data. If no such
$k$ is found then p-values are displayed up to $k_{max}$. The $n$-dependent
size adjustment controls for the over-rejection problem induced by sequential
testing (see Section 6.2 below). Overall, the results point to a higher number
of latent factors during bear market phases compared to bull market phases and
a decrease of the number of factors over time.232323We also investigate
stability of the factor structure by dividing each window of 20 months into
two overlapping subperiods of 16 months (overlap of 12 months) and by
estimating canonical correlations between the betas in each subperiod (see
SMC). We find that the fraction of the number of latent factors which are
common factors is 1 in 70% of the windows. The fraction is between 0.8 and 1
in 25% of the windows. It is between 0.5 and 0.8 in the remaining periods. It
remains true for the three-month recession periods 1987/09-1987/11 and
2020/01-2020/03, which represent only a fraction of their respective
subperiods, although there are "bull" market periods finding a similar number
of latent factors. In particular, our results based on a fixed $T$ and large
$n$ approach contradict the comprehension of a single factor model during
market downturns due to estimated correlations between equities approaching
$1$. It is consistent with the presence of risk factors, such as tail risk or
liquidity risk, only showing in distress periods. A rise in the estimated $k$
often happens towards the end of the recession periods. It is consistent with
the methodology of Lunde and Timmermann (2004) being early in detecting bear
periods (early warning system). The average estimated number of factors is
around $7$, close to the 4 to 6 factors found by PCA in Bai and Ng (2006) on
large time spans of individual stocks.242424With fixed $T$, the selection
procedure of Zaffaroni (2019), being by construction more conservative than a
(multiple) testing procedure (see the discussion on p. 508 of Gagliardini,
Scaillet, and Ossola (2019)), yields a smaller number of factors. Imposing
cross-sectional independence (resp., Gaussianity and cross-sectional
independence) for the LR test gives most of the time an increase by 1 or 2 (1
or 3). We have an average increase of 3 under sphericity.

###  6.2 Decomposing the cross-sectional variance

Building on the results in Pötscher (1983), we can obtain a consistent
estimator of the number of latent factors in each subperiod by allowing the
asymptotic size $\alpha$ go to zero as $n\to\infty$ in the sequential testing
procedure. We let $\hat{k}$ be defined as the smallest nonnegative integer $k$
satisfying $\text{pval}(k)>\alpha_{n}$, where $\text{pval}(k)$ is the p-value
from testing $H_{0}(k)$, and $\alpha_{n}$ is a sequence in $[0,1]$ with
$\alpha_{n}\to 0$. In practice, we take $\alpha_{n}=10/n_{max}$.252525This
choice satisfies the theoretical rule $\log\alpha_{n}/n\to 0$ given in
Pötscher (1983). If no such $k$ is found after sequentially testing
$H_{0}(k)$, for $k=0,\dots,k_{max}$ at level $\alpha_{n}$, then we take
$\hat{k}=k_{max}+1$. The Monte Carlo results show that such a selection
procedure works well. We use the estimate $\hat{k}$ at each subperiod to
decompose the path of the cross-sectional variance of stock returns into its
systematic and idiosyncratic parts:
$\hat{V}_{y,tt}=\hat{F}_{t}^{\prime}\hat{F}_{t}+\hat{V}_{\varepsilon,tt}$,
where $\hat{F}$ and $\hat{V}_{\varepsilon}$ are the FA estimates obtained by
extracting $\hat{k}$ latent factors. The condition (FA1) ensures that the
decomposition holds for any $t$. Such a decomposition is invariant to the
choice of normalization for the latent factors. If we look at time averages on
a subperiod, we get the decomposition
$\overline{\hat{V}}_{y}=\overline{\hat{F}^{\prime}\hat{F}}+\overline{\hat{V}}_{\varepsilon}$,
where the overline indicates averaging
$\hat{V}_{y,tt}=\hat{F}_{t}^{\prime}\hat{F}_{t}+\hat{V}_{\varepsilon,tt}$ on
$t$. In the lower panels of Figure 1, the blue dots correspond to the square
root of those quantities for the volatilities, while the ratios
$\hat{R}^{2}=\overline{\hat{F}^{\prime}\hat{F}}/\overline{\hat{V}}_{y}$ and
$\hat{R}^{2}$ under a single-factor model in the two last panels give measures
of goodness-of-fit.262626We do not plot the whole paths date $t$ by date $t$,
but only averages, for readability. If we sum over time instead of averaging
the estimated variances, we get a quantity similar to an integrated volatility
(see e.g. Barndorff-Nielsen and Shephard (2002), Andersen, Bollerslev,
Diebold, and Labys (2003), and references in Aït-Sahalia and Jacod (2014)),
and $\hat{R}^{2}$ is the ratio of such quantities. We can observe an uptrend
in total and idiosyncratic volatilities, while the systematic volatility
appears to remain stable over time even if the number of factors has overall
decreased over time.272727Cross-sectional independence increases estimated
systematic risk in average by 0.6% and decreases estimated idiosyncratic
volatility in average by 0.5%, so that estimated $R^{2}$ is inflated in
average by 4% per month. We get the same magnitude under cross-sectional
independence and Gaussianity. Under sphericity, PCA estimates give an increase
of 1%, a decrease of 1.3%, and an increase of 12.8% for the same quantities.
As a result, $\hat{R}^{2}$ is lower on average after the year 2000, indicating
a noisier environment. During the 2007-2008 financial crisis, we can observe a
rise in systematic volatility, causing $\hat{R}^{2}$ to reach 59% during that
period. In bear markets, $\hat{R}^{2}$ is often higher. It means that over a
bear subperiod, the systematic risk explains a large part of the cross-
sectional total variance even if it is not driven by a single factor as
reported in Section 6.1. The lowest panel in Figure 1 also signals that
$\hat{R}^{2}$ under the constraint of a single-factor model can be way below
the one given by the multifactor model. It also means that the idiosyncratic
volatility is overestimated if we use a single latent factor only. The plots
of the equal-weighted market and firm volatilities used as measures of total
and idiosyncratic volatility from a single observed factor (CAPM)
decomposition in Campbell et al. (2023) show similar patterns as our panels in
Figure 1.282828As in Campbell et al. (2023), we have also made the estimation
on value-weighted returns and we confirm that the results are qualitatively
similar. Section 4 of Campbell et al. (2023) discusses economic forces (firm
fundamentals and investor sentiments) driving the observed time-series
variation in average idiosyncratic volatility.

###  6.3 Spanning with observed factors

As discussed in Bai and Ng (2006), we get economic interpretation of latent
factors with observed factors when we have spanning between the latent factors
and the observed factors to be used as proxies in asset pricing (Shanken
(1992)). When $n$ and $T$ are large, Bai and Ng (2006) exploit the asymptotic
normality of the empirical canonical correlations between the two sets of
factors to investigate spanning under a symmetric role of the two sets. When
$T$ is fixed, we suggest the following strategy based on testing for the rank
of a matrix. Let us consider $k^{O}\geq k$ empirical factors that are excess
returns of portfolios,292929If $k^{O}<k$, empirical factors cannot span the
latent space by construction. The condition $k^{O}\geq k$ eases discussion but
is not needed for the rank tests. and let $\hat{F}^{O}$ denote the $T\times
k^{O}$ matrix of their values with row $t$ given by the transpose of
$\hat{f}^{O}_{t}=\frac{1}{n}\sum_{i=1}^{n}(y_{i,t}-r_{f,t})z_{i,t}$, where
$\frac{1}{n}z_{i,t}$ is a $k^{O}\times 1$ vector of time-varying portfolio
weights (long or short positions) based on stocks characteristics. Let matrix
$F^{O}$ with rows
$f_{t}^{O}=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum_{i=1}^{n}E[(y_{i,t}-r_{f,t})z_{i,t}]$
be the corresponding large-$n$ population limit. The notation $\hat{F}^{O}$
makes clear that the sample average of weighted excess returns is an estimate
of the population values $F^{O}$. We need to take this into account in the
asymptotic analysis of the rank test statistics when $n\rightarrow\infty$.
From the factor model under NA,
$y_{i,t}=r_{f,t}+f_{t}^{\prime}\tilde{\beta}_{i}+\varepsilon_{i,t}$ (see
footnote 4), and assuming cross-sectional non-correlation of idiosyncratic
errors and portfolio weights, we get $F^{O}=F\Phi^{\prime}$, where
$\Phi=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum_{i=1}^{n}E[z_{i,t}]\tilde{\beta}_{i}^{\prime}$
is assumed independent of $t$, $t=1,...,T$. Hence, the range of $F^{O}$ is a
subset of the range of $F$, namely the latent factors span the observed
factors (in the population limit sense) by construction. Moreover,
$\text{Rank}(F^{O})\leq k$. We can test the null hypothesis that $F$ and
$F^{O}$ span the same linear spaces, namely matrices $F$ and $F^{0}$ have the
same range. Such a null hypothesis is equivalent to the rank condition:
$\text{Rank}(F^{O})=k$.

We build on the rank testing literature; see e.g. Cragg and Donald (1996),
Robin and Smith (RS, 2000), Kleibergen and Paap (KP, 2006), Al-Sadoon
(2017).303030Ahn, Horenstein and Wang (2018) use that technology in a
fixed-$n$ large-$T$ setting, and find that ranks of beta matrices estimated
from either portfolios, or individual stocks, excess returns are often
substantially smaller than the (potentially large) number $k^{O}$ of observed
factors. The explanation in large economies is that the portfolio beta
matrices coincide with $\Phi$, and thus they cannot have a rank above the
(potentially small) number $k$ of latent factors. We use in particular the RS
and KP statistics. For those tests, the null hypothesis is that a given matrix
has a reduced rank $r$ against the alternative hypothesis that the rank is
greater than $r$. Hence, to test for spanning by the empirical factors, we
consider the null hypothesis $H_{0,sp}(r):\text{Rank}(F^{O})=r$ against the
alternative hypothesis $H_{1,sp}(r):\text{Rank}(F^{O})>r$, for any integer
$r<k$.313131Spanning holds if we can reject $H_{0,sp}(r)$ for any $r<k$. We
use the asymptotic expansion
$\hat{F}^{O}=\tilde{F}^{O}+\frac{1}{\sqrt{n}}\Psi_{F^{O},n}$, where
$\tilde{F}^{O}=F\Phi_{n}^{\prime}$ with
$\Phi_{n}=\frac{1}{n}\sum_{i=1}^{n}E[z_{i,t}]\tilde{\beta}_{i}^{\prime}$, and
the rows of matrix $\Psi_{F^{O},n}$ are given by
$\Psi_{F,n,t}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\eta_{i,t}f_{t}+\varepsilon_{i,t}z_{i,t})$
with $\eta_{i,t}:=(z_{i,t}-E[z_{i,t}])\tilde{\beta}_{i}^{\prime}$. Under
$H_{0,sp}(r)$, we assume that $\Phi_{n}$ has the same null space as $\Phi$, in
particular $\Phi_{n}$ has rank $r$, for $n$ large enough.323232Under
$\text{Rank}(F)=k$, we have $\text{Rank}(F^{O})=\text{Rank}(\Phi)$. Hence,
under $H_{0,sp}(r)$, matrix $\Phi$ has reduced rank $r$.  We assume the CLT
$vec(\Psi_{F^{O},n})\Rightarrow N(0,\Omega_{\Psi})$. Further, we use the
Singular Value Decomposition (SVD) of matrix
$\hat{F}^{O}=\hat{U}\hat{\mathcal{S}}\hat{V}^{\prime}$. Then, the RS and KP
statistics are the quadratic forms
$\mathscr{S}_{RS}=nvec(\hat{\mathcal{S}}_{22})^{\prime}vec(\hat{\mathcal{S}}_{22})$
and
$\mathscr{S}_{KP}=nvec(\hat{\mathcal{S}}_{22})^{\prime}\hat{\Omega}^{-1}_{\mathcal{S}}vec(\hat{\mathcal{S}}_{22})$,
where $\hat{\mathcal{S}}_{22}$ is the lower-right $(T-r)\times(k^{O}-r)$ block
of matrix $\hat{\mathcal{S}}$. Here,
$\hat{\Omega}_{\mathcal{S}}=(\hat{V}_{k^{O}-r}\otimes\hat{U}_{T-r})^{\prime}\hat{\Omega}_{\Psi}(\hat{V}_{k^{O}-r}\otimes\hat{U}_{T-r})$,
where $\hat{U}_{T-r}$ and $\hat{V}_{k^{O}-r}$ are the $T\times(T-r)$ and
$k^{O}\times(k^{O}-r)$ matrices in the block forms
$\hat{U}=[\hat{U}_{r}:\hat{U}_{T-r}]$ and
$\hat{V}=[\hat{V}_{r}:\hat{V}_{k^{O}-r}]$. In the SMC, we design a consistent
estimator $\hat{\Omega}_{\Psi}$ of $\Omega_{\Psi}$ building on a block
structure for the characteristics akin to Assumption 2 and a stationarity
condition. The definitions of the test statistics $\mathscr{S}_{RS}$ and
$\mathscr{S}_{KP}$ are equivalent to those in the original RS and KP papers.
The asymptotic distributions under $H_{0,sp}(r)$ are
$\mathscr{S}_{RS}\Rightarrow\sum_{j=1}^{(T-r)(k^{O}-r)}\delta^{\mathcal{S}}_{j}\chi^{2}_{j}(1)$
and $\mathscr{S}_{KP}\Rightarrow\chi^{2}[(T-r)(k^{O}-r)]$, where the
$\delta^{\mathcal{S}}_{j}$ are the eigenvalues of matrix
$\Omega_{\mathcal{S}}=(V_{k^{O}-r}\otimes
U_{T-r})^{\prime}\Omega_{\Psi}(V_{k^{O}-r}\otimes U_{T-r})$, assumed non-
singular.

We build the empirical matrix $\hat{F}^{O}$ with the time-varying portfolio
weights of the Fama-French five-factor model (Fama and French (2015)) plus the
momentum factor (Carhart (1997)), i.e., $k^{O}=6$. In the two panels of Figure
2, we can observe that the rank tests point most of the time at a low reduced
rank $r$ either 1 or 2, with only occasionally 3 or 4, for the matrix
$\hat{F}^{O}$. Observed factors struggle spanning latent factors since their
associated linear space is of a dimension smaller than the one of the latent
factor space. The discrepancy between the dimensions of the two factor spaces
has decreased over time. According to the KP statistic, the rank deficiency of
$\hat{F}^{O}$ is often less pronounced in bear markets indicating less
redundancy between the observed factors.

##  7 Concluding remarks

In this paper, we develop a new theory of latent Factor Analysis in short
panels beyond the Gaussian and i.i.d. cases. We establish the AUMPI property
of the LR statistic for testing hypotheses on the number of latent factors.
Our results for short subperiods of the CRSP panel of US stock returns
contradict the comprehension of a single factor during market downturns. In
bear markets, systematic risk driven by a latent multifactor structure
explains a large part of the cross-sectional variance, and is not spanned by
traditional empirical factors with a discrepancy between the dimensions of the
two factor spaces decreasing over time. The estimated paths of total and
idiosyncratic volatilities month after month feature an uptrend through time.

## References

Ahn, S., and Horenstein, A., 2013. Eigenvalue ratio test for the number of
factors. Econometrica 81 (3), 1203-1227.

Ahn, S., Horenstein, A., and Wang, N., 2018. Beta matrix and common factors in
stock returns. Journal of Financial and Quantitative Analysis 53 (3),
1417-1440.

Aigner, D., Hsiao, C., Kapteyn A., and Wansbeek, T., 1984. Latent variable
models in econometrics, in Handbook of Econometrics, Volume II, Z. Griliches
and M.D. Intriligator Eds., 1321-1393.

Aït-Sahalia, Y., and Jacod, J., 2014. High-frequency financial econometrics.
Princeton University Press.

Al-Sadoon, M., 2017. A unifying theory of tests of rank. Journal of
Econometrics 199 (1), 49-62.

Ang, A., Liu, J., and Schwarz, K., 2020. Using stocks or portfolios in tests
of factor models. Journal of Financial and Quantitative Analysis 55 (3),
709-750.

Andersen, T., Bollerslev, T., Diebold, F., and Labys, P., 2003. Modeling and
forecasting realized volatility. Econometrica 71 (2), 579-625.

Anderson, T. W., 2003. An introduction to multivariate statistical analysis.
Wiley.

Anderson, T. W., 1963. Asymptotic theory for Principal Components Analysis.
Annals of Mathematical Statistics 34, 122-148.

Anderson, T. W., and Rubin, H., 1956. Statistical inference in factor
analysis. Proceedings of the Third Berkeley Symposium in Mathematical
Statistics and Probability 5, 11-150.

Anderson, T. W. and Amemiya, Y., 1988. The asymptotic normal distribution of
estimators in factor analysis under general conditions. Annals of Statistics
16 (2), 759-771.

Ando, T., and Bai, J., 2015. Asset pricing with a general multifactor
structure. Journal of Financial Econometrics 13 (3), 556-604.

Bai, J., 2003. Inferential theory for factor models of large dimensions.
Econometrica 71 (1), 135-171.

Bai, J., 2009. Panel data models with interactive effects. Econometrica 77
(4), 1229-1279.

Bai, J., and Li, K., 2012. Statistical analysis of factor models of high
dimension. Annals of Statistics 40 (1), 436-465.

Bai, J., and Li, K., 2016. Maximum likelihood estimation and inference for
approximate factor models of high dimension. Review of Economics and
Statistics 98 (2), 298-309.

Bai, J., and Ng, S., 2002. Determining the number of factors in approximate
factor models. Econometrica 70 (1), 191-221.

Bai, J., and Ng, S., 2006. Evaluating latent and observed factors in
macroeconomics and finance. Journal of Econometrics 131 (1), 507-537.

Barigozzi, M., and Hallin, M., 2016. Generalized dynamic factor models and
volatilities: recovering the market volatility shocks. Econometrics Journal 19
(1), 33-60.

Barndorff-Nielsen, O., and Shephard, N., 2002. Econometric analysis of
realized volatility and its use in estimating stochastic volatility models.
Journal of the Royal Statistical Society, Series B, 64 (2), 253-280.

Bell, E. T., 1934. Exponential polynomials. Annals of Mathematics 35 (2),
258-277.

Bernstein, D., 2009. Matrix mathematics: Theory, facts and formulas. Princeton
University Press.

Bickel, P. J., and Levina, E., 2008. Covariance regularization by
thresholding. Annals of Statistics 36 (6), 2577-2604.

Campbell, J., Lettau, M., Malkiel, B., and Xu, Y., 2023. Idiosyncratic equity
risk two decades later. Critical Finance Review, 12, 203-223.

Caner, M., and Han, X., 2014. Selecting the correct number of factors in
approximate factor models: the large panel case with group bridge estimator.
Journal of Business and Economic Statistics 32 (3), 359-374.

Carhart, M., 1997. On persistence in mutual fund performance. Journal of
Finance. 52 (1), 57-82.

Chamberlain, G., 1992. Efficiency bounds for semi-parametric regression.
Econometrica 60 (3), 567-596.

Cochrane, J., 2005. Asset pricing. Princeton University Press.

Connor, G., and Korajczyk, R., 1986. Performance measurement with the
arbitrage pricing theory: A new framework for analysis. Journal of Financial
Economics 15 (3), 373-394.

Connor, G., and Korajczyk, R., 1993. A test for the number of factors in an
approximate factor model. Journal of Finance 48 (4), 1263-1291.

Cragg, J., and Donald, S., 1996. On the asymptotic properties of LDU-based
tests of the rank of a matrix. Journal of the American Statistical Association
91 (435), 1301-1309.

Eaton, M., 1987. Multivariate statistics. A vector space approach. Institute
of Mathematical Statistics Lecture notes-Monograph series, Vol. 53.

Engle, R., 1984. Wald, likelihood ratio, and Lagrange multiplier tests in
econometrics, in Handbook of Econometrics, Volume II, Z. Griliches and M.D.
Intriligator Eds., 775-826.

Fama, E., and French, K., 2015. A five-factor asset pricing model. Journal of
Financial Economics 116 (1), 1-22.

Fan, J., Furger, A., and Xiu, D., 2016. Incorporating global industrial
classification standard into portfolio allocation: A simple factor-based large
covariance matrix estimator with high-frequency data. Journal of Business and
Economic Statistics 34 (4), 489-503.

Fortin, A.-P., Gagliardini, P., and Scaillet, O., 2023. Eigenvalue tests for
the number of latent factors in short panels. Journal of Financial
Econometrics, nbad024.

Freyberger, J., 2018. Non-parametric panel data models with interactive fixed
effects. Review of Economic Studies 85 (3), 1824-1851.

Gagliardini, P., Ossola, E., and Scaillet, O., 2016. Time-varying risk premium
in large cross-sectional equity datasets. Econometrica 84 (3), 985-1046.

Gagliardini, P., Ossola, E., and Scaillet, O., 2019. A diagnostic criterion
for approximate factor structure. Journal of Econometrics 21 (2), 503-521.

Gagliardini, P., Ossola, E., and Scaillet, O., 2020. Estimation of large
dimensional conditional factor models in finance, in Handbook of Econometrics,
Volume 7A, S. Durlauf, L. Hansen, J. Heckman, and R. Matzkin Eds., 219-282.

Jöreskog, K., 1970. A general method for analysis of covariance structures.
Biometrika 57 (2), 239-251.

Kapetanios, G., 2010. A testing procedure for determining the number of
factors in approximate factor models with large datasets. Journal of Business
and Economic Statistics 28 (3), 397-409.

Kim, S., and Skoulakis, G., 2018. Ex-post risk premia estimation and asset
pricing tests using large cross-sections: the regression-calibration approach.
Journal of Econometrics 204 (2), 159-188.

Kleibergen, F., and Paap, R., 2006. Generalized reduced rank tests using the
singular value decomposition. Journal of Econometrics 133 (1), 97-126.

Kleibergen, F., and Zhan, Z., 2023. Identification-robust inference for risk
premia in short panels. Working paper, University of Amsterdam.

Kotz, S., Johnson, N., and Boyd, D., 1967. Series representations of
distributions of quadratic forms in Normal variables II. Non-central case.
Annals of Mathematical Statistics 38 (3), 838-848.

Lancaster, T., 2000. The incidental parameter problem since 1948. Journal of
Econometrics 95 (2), 391-413.

Lehmann, E., and Romano, D., 2005. Testing statistical hypotheses, Springer
Texts in Statistics.

Lunde, A., and Timmermann, A., 2004. Duration dependence in stock prices: An
analysis of bull and bear markets. Journal of Business and Economic Statistics
22 (3), 253-273.

Magnus, J., and Neudecker, H., 2007. Matrix differential calculus, with
applications in statistics and econometrics. Wiley.

Miravete, E., 2011. Convolution and composition of totally positive random
variables in economics. Journal of Mathematical Economics 47 (4), 479-490.

Moon, H.R., and Weidner, M., 2015. Linear regression for panel with unknown
number of factors as interactive fixed effects. Econometrica 83 (4),
1543-1579.

Neyman, J., and Scott, E., 1948. Consistent estimation from partially
consistent observations. Econometrica 16(1), 1-32.

Onatski, A., 2009. Testing hypotheses about the number of factors in large
factor models. Econometrica 77 (5), 1447-1479.

Onatski, A., 2010. Determining the number of factors from empirical
distribution of eigenvalues. Review of Economics and Statistics 92 (4),
1004-1016.

Onatski, A. 2023. Comment on “Eigenvalue tests for the number of latent
factors in short panels" by A.-P. Fortin, P. Gagliardini and O. Scaillet,
Journal of Financial Econometrics, nbad28.

Pesaran, M.H., 2006. Estimation and inference in large heterogeneous panels
with a multifactor error structure. Econometrica 74 (4), 967-1012.

Pötscher, B., 1983. Order estimation in ARMA-models by Lagrangian multiplier
tests, Annals of Statistics 11 (3), 872-885.

Raponi, V., Robotti, C., and Zaffaroni, P., 2020. Testing beta-pricing models
using large cross-sections. Review of Financial Studies 33 (6), 2796-2842.

Renault, E., Van Der Heijden, T., and Werker, B., 2023. Arbitrage pricing
theory for idiosyncratic variance factors. Journal of Financial Econometrics
21 (5), 1403-1442.

Robin, J.-M., and Smith, R., 2000. Tests of rank. Econometric Theory 16 (2),
151-175.

Romano, J., Shaikh, A., and Wolf, M., 2010. Hypothesis testing in
econometrics. Annual Review of Economics 2, 75-104.

Shanken, J. 1992. On the estimation of beta pricing models. Review of
Financial Studies 5 (1), 1-33.

Stock, J., and Watson, M., 2002. Forecasting using principal components from a
large number of predictors. Journal of the American Statistical Association 97
(460), 1167-1179.

Tao, T. 2012. Topics in random matrix theory. Graduate Studies in Mathematics,
Volume 132, American Mathematical Society.

White, H., 1982. Maximum likelihood estimation of misspecified models.
Econometrica 50 (1), 1-25.

Zaffaroni, P., 2019. Factor models for conditional asset pricing. Working
Paper.

Appendix

## A Regularity assumptions

In this appendix, we list and comment the additional assumptions used to
derive the large sample properties of the estimators and test statistics. We
often denote by $C>0$ a generic constant, and we use $\delta_{j}(A)$ to denote
the $j$th largest eigenvalue of a symmetric matrix $A$. Set $\Theta$ is a
compact subset of
$\\{\theta=(vec(F)^{\prime},diag(V_{\varepsilon})^{\prime})^{\prime}\in\mathbb{R}^{r}~{}:~{}V_{\varepsilon}~{}\text{is
diagonal and positive definite},~{}F^{\prime}V_{\varepsilon}^{-1}F~{}\text{is
diagonal, with}$$\text{diagonal elements ranked in decreasing order}\\}$ with
$r=(T+1)k$, and function
$L_{0}(\theta)=-\frac{1}{2}\log|\Sigma(\theta)|-\frac{1}{2}Tr\left(V_{y}\Sigma(\theta)^{-1}\right)$
is the population FA criterion, where
$\Sigma(\theta)=FF^{\prime}+V_{\varepsilon}$ and
$V_{y}=\underset{n\rightarrow\infty}{\text{plim}}~{}\hat{V}_{y}$. Further,
$\theta_{0}=(vec(F_{0})^{\prime},diag(V_{\varepsilon}^{0})^{\prime})^{\prime}$
denotes the vector of true parameter values under $H_{0}(k)$ and is an
interior point of set $\Theta$.

###### Assumption A.1

The non-zero eigenvalues of $V_{y}V_{\varepsilon}^{-1}-I_{T}$ are distinct,
i.e., $\gamma_{1}>...>\gamma_{k}>0$.

###### Assumption A.2

The loadings are normalized such that
$\bar{\beta}=\frac{1}{n}\sum_{i=1}^{n}\beta_{i}=0$ and
$\tilde{V}_{\beta}:=\frac{1}{n}\sum_{i=1}^{n}\beta_{i}\beta_{i}^{\prime}\break=I_{k}$,
for any $n$. Moreover, $|\beta_{i}|\leq C$, for all $i$.

###### Assumption A.3

We have $E[w_{i,t}^{8}]\leq C$ and $|\sigma_{i,j}|\leq C$, for all $i,j,t$.

###### Assumption A.4

Under the null hypothesis $H_{0}(k)$, we have:
$\Sigma(\theta)=\Sigma(\theta_{0})$, $\theta\in\Theta$ $\Rightarrow$
$\theta=\theta_{0}$, up to sign changes in the columns of $F$.

###### Assumption A.5

Matrix $M_{F_{0},{V}^{0}_{\varepsilon}}^{\odot 2}$ is non-singular.

###### Assumption A.6

(a) The $\frac{T(T+1)}{2}\times\frac{T(T+1)}{2}$ symmetric matrix
$D=\underset{n\rightarrow\infty}{\lim}D_{n}$ exists, where
$D_{n}=\frac{1}{n}\sum_{i=1}^{n}\sigma_{ii}^{2}V[vech(w_{i}w_{i}^{\prime})]$.
(b) We have
$\delta_{T(T+1)/2}\left(V[vech(w_{i}w_{i}^{\prime})]\right)\geq\underline{c}$,
for all $i\in\bar{S}$, where $\bar{S}\subset\\{1,...,n\\}$ with
$\frac{1}{n}\sum_{i=1}^{n}1_{i\in\bar{S}}\geq 1-\frac{1}{2\bar{C}}$, for
constants $\bar{C},\bar{c}>0$, such that $\sigma_{ii}\leq\bar{C}$. (c) We have
$\underset{n\rightarrow\infty}{\lim}\kappa_{n}=\kappa$ for a constant
$\kappa\geq 0$, where
$\kappa_{n}:=\frac{1}{n}\sum_{m=1}^{J_{n}}\left(\sum_{i\neq j\in
I_{m}}\sigma_{ij}^{2}\right)$.

###### Assumption A.7

Under the alternative hypothesis $H_{1}(k)$, (a) function $L_{0}(\theta)$ has
a unique maximizer
$\theta^{*}=(vec(F^{*})^{\prime},diag(V_{\varepsilon}^{*})^{\prime})^{\prime}$
over $\Theta$, and (b) we have $V_{y}\neq
F^{*}(F^{*})^{\prime}+V^{*}_{\varepsilon}$.

###### Assumption A.8

Matrix
$Q_{\beta}:=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\beta^{\prime}\Sigma\beta$
is positive definite.

Assumption A.1 removes the rotational indeterminacy (up to sign) that occurs
when identifying the columns of matrix $F$ as eigenvectors of
$V_{y}V^{-1}_{\varepsilon}$. Assumptions A.2 and A.3 require uniform bounds on
factor loadings as well as on covariances and higher-order moments of the
idiosyncratic errors. Assumption A.4 implies global identification in the FA
model (see Lemma 4). Assumptions A.2-A.4 yield consistency of FA estimators
(see proof of Lemma 5). Assumption A.5 is the local identification condition
in the FA model (see Lemma 6). We use Assumption A.6 together with Assumption
A.3 to invoke a CLT based on a multivariate Lyapunov condition (see proof of
Lemma 1) to establish the asymptotic distribution of the LR statistic. To ease
the verification of the Lyapunov condition, we bound a fourth-order moment of
squared errors, which explains why we require finite eight-order moments in
Assumption A.3. We could relax this condition at the expense of a more
sophisticated proof of Lemma 1. The mild Assumption A.6 (b) requires that the
smallest eigenvalue of $V[vech(w_{i}w_{i}^{\prime})]$ is bounded away from $0$
for all assets $i$ up to a small fraction. In Assumption A.6 (c), in order to
have $\kappa_{n}$ bounded, we need either mixing dependence in idiosyncratic
errors within blocks, i.e., $|\sigma_{i,j}|\leq C\rho^{|i-j|}$ for $i,j\in
I_{m}$ and $0\leq\rho<1$, or vanishing correlations, i.e., $|\sigma_{i,j}|\leq
Cb_{m,n}^{-\bar{s}}$ for all $i\neq j\in I_{m}$ and a constant $\bar{s}\geq
1/2$, with blocks of equal size. In Assumption A.7, part (a) defines the
pseudo-true parameter value (White (1982)) under the alternative hypothesis,
and part (b) is used to establish the consistency of the LR test under global
alternative hypotheses (see proof of Proposition 2). Finally, Assumption A.8
is used to apply a Lyapunov CLT (see proof of Lemma 7) when deriving the
asymptotic normality of the FA estimators.

## B Proofs of Propositions 1-5

Proof of Proposition 1: Let $\hat{U}$ be the $T\times k$ matrix whose
orthonormal columns are the eigenvectors for the $k$ largest eigenvalues of
matrix $\hat{V}_{\varepsilon}^{-1/2}\hat{V}_{y}\hat{V}_{\varepsilon}^{-1/2}$.
Those eigenvalues are $1+\hat{\gamma}_{j}$, $j=1,...,k$, while it holds
$\hat{F}=\hat{V}_{\varepsilon}^{1/2}\hat{U}\hat{\Gamma}^{1/2}$, where
$\hat{\Gamma}=diag(\hat{\gamma}_{1},...,\hat{\gamma}_{k})$. We have
$I_{T}-\hat{U}\hat{U}^{\prime}=I_{T}-\hat{V}_{\varepsilon}^{-1/2}\hat{F}\hat{\Gamma}^{-1}\hat{F}^{\prime}\hat{V}_{\varepsilon}^{-1/2}$
$=I_{T}-\hat{V}_{\varepsilon}^{-1/2}\hat{F}(\hat{F}^{\prime}\hat{V}_{\varepsilon}^{-1}\hat{F})^{-1}\hat{F}^{\prime}\hat{V}_{\varepsilon}^{-1/2}$
$=\hat{V}_{\varepsilon}^{-1/2}M_{\hat{F},\hat{V}_{\varepsilon}}\hat{V}_{\varepsilon}^{1/2}=\hat{V}_{\varepsilon}^{1/2}M_{\hat{F},\hat{V}_{\varepsilon}}^{\prime}\hat{V}_{\varepsilon}^{-1/2}.$
Thus,
$\hat{S}=(I_{T}-\hat{U}\hat{U}^{\prime})\left(\hat{V}_{\varepsilon}^{-1/2}\hat{V}_{y}\hat{V}_{\varepsilon}^{-1/2}-I_{T}\right)(I_{T}-\hat{U}\hat{U}^{\prime}).$
By the spectral decomposition of
$\hat{V}_{\varepsilon}^{-1/2}\hat{V}_{y}\hat{V}_{\varepsilon}^{-1/2}$, we get
$(I_{T}-\hat{U}\hat{U}^{\prime})\left(\hat{V}_{\varepsilon}^{-1/2}\hat{V}_{y}\hat{V}_{\varepsilon}^{-1/2}-I_{T}\right)$
$(I_{T}-\hat{U}\hat{U}^{\prime})=\sum_{j=k+1}^{T}\hat{\gamma}_{j}\hat{P}_{j}$,
where the $\hat{P}_{j}$ are the orthogonal projection matrices onto the
eigenspaces for the $T-k$ smallest eigenvalues. Then, Part (a) follows. Part
(b) is a consequence of the squared Frobenius norm of a symmetric matrix being
equal to the sum of its squared eigenvalues. For part (c), let
$P_{\hat{F},{\hat{V}}_{\varepsilon}}=I_{T}-M_{\hat{F},\hat{V}_{\varepsilon}}$
and note that
$\hat{F}\hat{F}^{\prime}=P_{\hat{F},{\hat{V}}_{\varepsilon}}(\hat{V}_{y}-\hat{V}_{\varepsilon})+(\hat{V}_{y}-\hat{V}_{\varepsilon})P_{\hat{F},{\hat{V}}_{\varepsilon}}^{\prime}-P_{\hat{F},{\hat{V}}_{\varepsilon}}(\hat{V}_{y}-\hat{V}_{\varepsilon})P_{\hat{F},{\hat{V}}_{\varepsilon}}^{\prime}$
$=\hat{V}_{y}-\hat{V}_{\varepsilon}-M_{\hat{F},\hat{V}_{\varepsilon}}(\hat{V}_{y}-\hat{V}_{\varepsilon})M_{\hat{F},\hat{V}_{\varepsilon}}^{\prime},$
where the first equality is because the three terms on the RHS are all equal
to $\hat{F}\hat{F}^{\prime}$ by (FA2). The conclusion follows from (FA1) and
$\hat{V}_{\varepsilon}$ being diagonal. Finally, part (d) follows since
$\frac{1}{n}\hat{\varepsilon}\hat{\varepsilon}^{\prime}=M_{\hat{F},\hat{V}_{\varepsilon}}\hat{V}_{y}M_{\hat{F},\hat{V}_{\varepsilon}}^{\prime}$,
$\hat{V}_{\varepsilon}M_{\hat{F},\hat{V}_{\varepsilon}}^{\prime}=M_{\hat{F},\hat{V}_{\varepsilon}}\hat{V}_{\varepsilon}$,
and $M_{\hat{F},\hat{V}_{\varepsilon}}$ is idempotent, which implies
$\hat{V}_{\varepsilon}^{-1/2}M_{\hat{F},\hat{V}_{\varepsilon}}\hat{V}_{\varepsilon}M_{\hat{F},\hat{V}_{\varepsilon}}^{\prime}\hat{V}_{\varepsilon}^{-1/2}=\hat{V}_{\varepsilon}^{-1/2}M_{\hat{F},\hat{V}_{\varepsilon}}\hat{V}_{\varepsilon}^{1/2}$.

Proof of Proposition 2: (a) Let us write
$\hat{G}=(\hat{g}_{t,i})=[\hat{g}_{1}:\cdots:\hat{g}_{T-k}]$. By definition of
the $vech$ operator, we have
$vech(\hat{G}^{\prime}E_{t,t}\hat{G})=\left[\frac{1}{\sqrt{2}}\hat{g}^{2}_{t,1}~{}:~{}\cdots~{}:~{}\frac{1}{\sqrt{2}}\hat{g}^{2}_{t,T-k}~{}:~{}\\{\hat{g}_{t,i}\hat{g}_{t,j}\\}_{i<j}\right]^{\prime}$.
Therefore, we have
$\hat{\boldsymbol{X}}^{\prime}=\left[\frac{1}{\sqrt{2}}\hat{g}_{1}\odot\hat{g}_{1}~{}:~{}\cdots~{}:~{}\frac{1}{\sqrt{2}}\hat{g}_{T-k}\odot\hat{g}_{T-k}~{}:~{}\\{\hat{g}_{i}\odot\hat{g}_{j}\\}_{i<j}\right]$,
where $\odot$ denotes the Hadamard product (i.e., element-wise matrix
product). It follows that, for any $(T-k)\times(T-k)$ symmetric matrix
$A=(a_{i,j})$,
$diag(\hat{G}A\hat{G}^{\prime})=\sum_{i=1}^{T-k}a_{i,i}diag(\hat{g}_{i}\hat{g}^{\prime}_{i})+2\sum_{i<j}a_{i,j}diag(\hat{g}_{i}\hat{g}^{\prime}_{j})=\sum_{i=1}^{T-k}a_{i,i}(\hat{g}_{i}\odot\hat{g}_{i})+2\sum_{i<j}a_{i,j}(\hat{g}_{i}\odot\hat{g}_{j})=2\hat{\boldsymbol{X}}^{\prime}vech(A)$.

Let us now show the asymptotic expansion (5). Since
$\hat{G}\hat{G}^{\prime}\hat{V}^{-1}_{\varepsilon}=M_{\hat{F},\hat{V}_{\varepsilon}}$,
we have
$\hat{S}=\hat{V}^{-1/2}_{\varepsilon}\hat{G}\hat{S}^{*}\hat{G}^{\prime}\hat{V}^{-1/2}_{\varepsilon}$,
where
$\hat{S}^{*}=\hat{G}^{\prime}\hat{V}^{-1}_{\varepsilon}(\hat{V}_{y}-\hat{V}_{\varepsilon})\hat{V}^{-1}_{\varepsilon}\hat{G}$.
Therefore, from Proposition 1 (c) we get
$0=diag(\hat{S})=\hat{V}^{-1}_{\varepsilon}diag(\hat{G}\hat{S}^{*}\hat{G}^{\prime})=2\hat{V}^{-1}_{\varepsilon}\hat{\boldsymbol{X}}^{\prime}vech(\hat{S}^{*})$,
i.e., the condition $diag(\hat{S})=0$ is equivalent to $vech(\hat{S}^{*})$
being in the orthogonal complement of the range of $\hat{\boldsymbol{X}}$. It
follows from the local identification assumption A.5 that
$vech(\hat{S}^{*})=M_{\hat{\boldsymbol{X}}}vech(\hat{S}^{*})$.333333Assumption
A.5 is equivalent to $\boldsymbol{X}$ having full column rank by Lemma 6 in
Appendix E.4. Besides, from Proposition 9 in Appendix E.6 and the fact that
$\hat{G}\hat{O}=G+o_{p}(1)$ for some rotation matrix $\hat{O}$ (see below), we
have $\hat{\mathscr{R}}\hat{\boldsymbol{X}}=\boldsymbol{X}+o_{p}(1)$ from some
orthogonal matrix $\hat{\mathscr{R}}$ . Hence, $\hat{\boldsymbol{X}}$ is
invertible with probability approaching $1$, and $M_{\hat{\boldsymbol{X}}}$ is
well-defined.  Next, we have

$M_{\hat{\boldsymbol{X}}}vech(\hat{S}^{*})=M_{\hat{\boldsymbol{X}}}vech(\hat{G}^{\prime}\hat{V}^{-1}_{\varepsilon}(\hat{V}_{y}-\hat{V}_{\varepsilon})\hat{V}^{-1}_{\varepsilon}\hat{G})=M_{\hat{\boldsymbol{X}}}vech(\hat{G}^{\prime}\hat{V}^{-1}_{\varepsilon}(\hat{V}_{y}-\tilde{V}_{\varepsilon})\hat{V}^{-1}_{\varepsilon}\hat{G}),$
(B.1)

because the kernel of $M_{\hat{\boldsymbol{X}}}$ is
$\\{vech(\hat{G}^{\prime}D\hat{G}):D\text{ diagonal}\\}$. Besides, we have the
expansion
$\sqrt{n}vech(\hat{G}^{\prime}\hat{V}^{-1}_{\varepsilon}(\hat{V}_{y}-\tilde{V}_{\varepsilon})\hat{V}^{-1}_{\varepsilon}\hat{G})=vech(\hat{Z}^{*}_{n})+o_{p}(1)$,
where
$\hat{Z}^{*}_{n}=\hat{G}^{\prime}\hat{V}^{-1}_{\varepsilon}Z_{n}\hat{V}^{-1}_{\varepsilon}\hat{G}$.
It is because
$\hat{V}_{y}-\tilde{V}_{\varepsilon}=FF^{\prime}+\frac{1}{\sqrt{n}}\Psi_{y}+o_{p}(\frac{1}{\sqrt{n}})$,
with $\Psi_{y}=\frac{1}{\sqrt{n}}(\varepsilon\beta
F^{\prime}+F\beta^{\prime}\varepsilon^{\prime})+\sqrt{n}\left(\frac{1}{n}\varepsilon\varepsilon^{\prime}-\tilde{V}_{\varepsilon}\right)$
(see Equation (E.2) and Lemma 5 in Appendix E.2), and
$\hat{G}^{\prime}\hat{V}^{-1}_{\varepsilon}F=\hat{G}^{\prime}\hat{V}^{-1}_{\varepsilon}M_{\hat{F},\hat{V}_{\varepsilon}}F=O_{p}(\frac{1}{\sqrt{n}})$
by the root-$n$ consistency of FA estimators (see Appendix E.5.1). Using
$\|\hat{S}\|^{2}=\|\hat{S}^{*}\|^{2}$, it follows that

$\frac{n}{2}\|\hat{S}\|^{2}=nvech(\hat{S}^{*})^{\prime}vech(\hat{S}^{*})=nvech(\hat{S}^{*})^{\prime}M_{\hat{\boldsymbol{X}}}vech(\hat{S}^{*})=vech(\hat{Z}^{*}_{n})^{\prime}M_{\hat{\boldsymbol{X}}}vech(\hat{Z}^{*}_{n})+o_{p}(1).$
(B.2)

From $M_{\hat{F},\hat{V}_{\varepsilon}}=M_{F,V_{\varepsilon}}+o_{p}(1)$, we
have $\hat{G}\hat{O}=G+o_{p}(1)$ for some (possibly data-dependent)
$(T-k)\times(T-k)$ orthogonal matrix $\hat{O}$. Since
$vech(\hat{Z}^{*}_{n})^{\prime}M_{\hat{\boldsymbol{X}}}vech(\hat{Z}^{*}_{n})$
is invariant to post-multiplication of $\hat{G}$ by an orthogonal matrix (see
Proposition 9 in Appendix E.6), from (B.2) we get
$\frac{n}{2}\|\hat{S}\|=vech(Z^{*}_{n})^{\prime}M_{\boldsymbol{X}}vech(Z^{*}_{n})+o_{p}(1)$,
which is the asymptotic expansion (5).

Let us now establish the asymptotic normality of $vech(Z^{*}_{n})$. For any
integer $m$, we let $A_{m}$ denote the unique $m^{2}\times\frac{m(m+1)}{2}$
matrix satisfying $vec(S)=A_{m}vech(S)$ for any $m\times m$ symmetric matrix
$S$.343434The explicit form for $A_{m}$ is $A_{m}=\left[\sqrt{2}(e_{1}\otimes
e_{1}):\cdots:\sqrt{2}(e_{m}\otimes e_{m}):\\{e_{i}\otimes e_{j}+e_{j}\otimes
e_{i}\\}_{i<j}\right]$, with $e_{i}$ being the $i$th unit vector of dimension
$m$. Matrix $A_{m}$ satisfies $A_{m}^{\prime}A_{m}=2I_{\frac{m(m+1)}{2}}$,
$A_{m}A_{m}^{\prime}=I_{m^{2}}+K_{m,m}$, and $K_{m,m}A_{m}=A_{m}$, where
$K_{m,m}$ is the commutation matrix (see also Magnus, Neudecker (2007) Theorem
12 in Chapter 2.8). Then, we have
$vech(Z^{*}_{n})=\boldsymbol{R}^{\prime}vech(\mathscr{Z}_{n})$, where
$\mathscr{Z}_{n}=V_{\varepsilon}^{-1/2}Z_{n}V_{\varepsilon}^{-1/2}$,
$\boldsymbol{R}=\frac{1}{2}A^{\prime}_{T}(Q\otimes Q)A_{T-k}$, and
$Q=V_{\varepsilon}^{-1/2}G$. Matrix $\boldsymbol{R}$ satisfies
$\boldsymbol{R}^{\prime}\boldsymbol{R}=I_{p}$. The next lemma establishes the
asymptotic normality of $vech(\mathscr{Z}_{n})$.

###### Lemma 1

(a) Under Assumptions 1-2, A.2, A.6 (a)-(b), we have
$\Omega_{n}^{-1/2}vech(\mathscr{Z}_{n})\Rightarrow N(0,I_{\frac{T(T+1)}{2}})$
as $n\rightarrow\infty$ and $T$ is fixed, where
$\Omega_{n}=D_{n}+\kappa_{n}I_{\frac{T(T+1)}{2}}$, and
$\kappa_{n}=\frac{1}{n}\sum_{m=1}^{J_{n}}\left(\sum_{i\neq j\in
I_{m}}\sigma_{ij}^{2}\right)$. If additionally Assumption A.6 (c) holds, then
$vech(\mathscr{Z}_{n})\Rightarrow N(0,\Omega)$, with $\Omega:=D+\kappa
I_{\frac{T(T+1)}{2}}$.

Lemma 1 yields the asymptotic normality of $vech(Z^{*}_{n})$, namely
$vech(Z^{*}_{n})\Rightarrow N(0,\Omega_{Z^{*}})$, with
$\Omega_{Z^{*}}=\boldsymbol{R}^{\prime}\Omega\boldsymbol{R}$. Part (a) then
follows from expansion (5) and the standard result on the distribution of
idempotent quadratic forms of Gaussian vectors.

(b) We have $\hat{z}^{*}_{m,n}=\sum_{i\in
I_{m}}\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}\left(\tilde{y}_{i}\tilde{y}_{i}^{\prime}\right)\hat{V}_{\varepsilon}^{-1}\hat{G}$
with $\tilde{y}_{i}=y_{i}-\bar{y}$, since
$\hat{\varepsilon}_{i}=M_{\hat{F},\hat{V}_{\varepsilon}}\tilde{y}_{i}$ and
$\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}M_{\hat{F},\hat{V}_{\varepsilon}}=\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}$.
We get $\hat{z}^{*}_{m,n}=\sum_{i\in
I_{m}}\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}\left(\tilde{\varepsilon}_{i}\tilde{\varepsilon}_{i}^{\prime}\right)\hat{V}_{\varepsilon}^{-1}\hat{G}+\sum_{i\in
I_{m}}\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}\left(F\beta_{i}\beta_{i}^{\prime}F^{\prime}\right)\hat{V}_{\varepsilon}^{-1}\hat{G}\break+\sum_{i\in
I_{m}}\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}\left(F\beta_{i}\tilde{\varepsilon}_{i}^{\prime}+\tilde{\varepsilon}_{i}\beta_{i}^{\prime}F^{\prime}\right)\hat{V}_{\varepsilon}^{-1}\hat{G}=:\tilde{z}^{*}_{m,n}+z^{*}_{m,n,1}+z^{*}_{m,n,2}$,
where $\tilde{\varepsilon}_{i}=\varepsilon_{i}-\bar{\varepsilon}$ by using
$\tilde{y}_{i}\tilde{y}_{i}^{\prime}=\tilde{\varepsilon}_{i}\tilde{\varepsilon}_{i}^{\prime}+F\beta_{i}\beta_{i}^{\prime}F^{\prime}+F\beta_{i}\tilde{\varepsilon}_{i}^{\prime}+\tilde{\varepsilon}_{i}\beta_{i}^{\prime}F^{\prime}$.
Then, we can decompose $\hat{\Omega}_{Z^{*}}$ into a sum of a leading term and
other terms, which are asymptotically negligible, so that
$\hat{\Omega}_{Z^{*}}=\tilde{\Omega}_{Z^{*}}+o_{p}(1)$, with
$\tilde{\Omega}_{Z^{*}}=\frac{1}{n}\sum_{m=1}^{J_{n}}vech(\bar{z}^{*}_{m,n})vech(\bar{z}^{*}_{m,n})^{\prime}$,
with $\bar{z}^{*}_{m,n}$ defined as $\tilde{z}^{*}_{m,n}$ after replacing
$\tilde{\varepsilon}_{i}$ with $\varepsilon_{i}$. Let us now show that
$M_{\hat{\boldsymbol{X}}}\tilde{\Omega}_{Z^{*}}M_{\hat{\boldsymbol{X}}}=M_{\boldsymbol{X}}\Omega_{Z^{*}}M_{\boldsymbol{X}}+o_{p}(1)$
up to pre- and post-multiplication by a rotation matrix and its inverse. We
have
$M_{\hat{\boldsymbol{X}}}vech\left(\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}\left(\varepsilon_{i}\varepsilon_{i}^{\prime}\right)\hat{V}_{\varepsilon}^{-1}\hat{G}\right)=M_{\hat{\boldsymbol{X}}}vech\left(\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}\left(\varepsilon_{i}\varepsilon_{i}^{\prime}-\sigma_{ii}V_{\varepsilon}\right)\hat{V}_{\varepsilon}^{-1}\hat{G}\right)$,
because of the kernel of $M_{\hat{\boldsymbol{X}}}$ . Moreover, from the
properties of matrix $A_{m}$ introduced in part (a), we have:
$vech\left(\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}(\varepsilon_{i}\varepsilon_{i}^{\prime}-\sigma_{ii}V_{\varepsilon})\hat{V}_{\varepsilon}^{-1}\hat{G}\right)\break=vech\left(\tilde{Q}^{\prime}(e_{i}e_{i}^{\prime}-\sigma_{ii}I_{T})\tilde{Q}\right)=\hat{\boldsymbol{R}}^{\prime}vech(e_{i}e_{i}^{\prime}-\sigma_{ii}I_{T}),$
where $e_{i}=V_{\varepsilon}^{-1/2}\varepsilon_{i}$, and
$\hat{\boldsymbol{R}}:=\frac{1}{2}A_{T}^{\prime}(\tilde{Q}\otimes\tilde{Q})A_{T-k}$
with $\tilde{Q}=V_{\varepsilon}^{1/2}\hat{V}_{\varepsilon}^{-1}\hat{G}$. We
get
$M_{\hat{\boldsymbol{X}}}vech(\bar{z}^{*}_{m,n})={M}_{\hat{\boldsymbol{X}}}\hat{\boldsymbol{R}}^{\prime}vech(\zeta_{m,n})$,
where $\zeta_{m,n}:=\sum_{i\in I_{m}}(e_{i}e_{i}^{\prime}-\sigma_{ii}I_{T})$.
Besides,
$vech(\mathscr{Z}_{n})=\frac{1}{\sqrt{n}}\sum_{m=1}^{J_{n}}vech(\zeta_{m,n})$.
Then,
$M_{\hat{\boldsymbol{X}}}\tilde{\Omega}_{Z^{*}}M_{\hat{\boldsymbol{X}}}\break={M}_{\hat{\boldsymbol{X}}}\hat{\boldsymbol{R}}^{\prime}\tilde{\Omega}_{n}\hat{\boldsymbol{R}}{M}_{\hat{\boldsymbol{X}}}$
for
$\tilde{\Omega}_{n}:=\frac{1}{n}\sum_{m=1}^{J_{n}}vech(\zeta_{m,n})vech(\zeta_{m,n})^{\prime}$.
Further, $E[\tilde{\Omega}_{n}]=V[vech(\mathscr{Z}_{n})]\break=\Omega_{n}$.
Moreover, $\tilde{\Omega}_{n}-E[\tilde{\Omega}_{n}]=o_{p}(1)$, by using
$vec(\tilde{\Omega}_{n})=\frac{1}{n}\sum_{m=1}^{J_{n}}vech(\zeta_{m,n})\otimes
vech(\zeta_{m,n})$ and $\|V[vec(\tilde{\Omega}_{n})]\|\leq
C\frac{1}{n^{2}}\sum_{m=1}^{J_{n}}E\left[\|vech(\zeta_{m,n})\|^{4}\right]=o(1)$,
where the latter bound is shown in the proof of Lemma 1 using Assumption 2
(d). Additionally, by Assumption A.6, we have $\Omega_{n}=\Omega+o(1)$. Thus,
$\tilde{\Omega}_{n}=\Omega+o_{p}(1)$. Now, from the proof of part (a) we have
$\hat{G}\hat{O}=G+o_{p}(1)$ for some $(T-k)\times(T-k)$ orthogonal matrix
$\hat{O}$. Then, by Proposition 9 (e) in Appendix E.6, we have
$\hat{\boldsymbol{R}}{M}_{\hat{\boldsymbol{X}}}\hat{\mathscr{R}}^{-1}=\boldsymbol{R}M_{\boldsymbol{X}}+o_{p}(1)$,
for a $p$ dimensional orthogonal matrix
$\hat{\mathscr{R}}\equiv\mathscr{R}(\hat{O})$. We conclude that
$\hat{\mathscr{R}}M_{\hat{\boldsymbol{X}}}\tilde{\Omega}_{Z^{*}}M_{\hat{\boldsymbol{X}}}\hat{\mathscr{R}}^{-1}$
is a consistent estimator of
$M_{\boldsymbol{X}}{\Omega}_{Z^{*}}M_{\boldsymbol{X}}$ as $n\rightarrow\infty$
and $T$ is fixed. Part (b) then follows from the continuity of eigenvalues for
symmetric matrices, and their invariance under pre- and post-multiplication by
an orthogonal matrix and its transpose.

(c) Under $H_{1}(k)$ and Assumption A.7 (a), we have
$\hat{F}\overset{p}{\rightarrow}F^{*}$ and
$\hat{V}_{\varepsilon}\overset{p}{\rightarrow}V_{\varepsilon}^{*}$. Then,
$\hat{S}\overset{p}{\rightarrow}S^{*}$ with
$S^{*}=(V_{\varepsilon}^{*})^{-1/2}M_{F^{*},V_{\varepsilon}^{*}}(V_{y}-V_{\varepsilon}^{*})M_{F^{*},V_{\varepsilon}^{*}}^{\prime}(V_{\varepsilon}^{*})^{-1/2}\neq
0$. Indeed, if $S^{*}$ were the null matrix, then we would have
$M_{F^{*},V_{\varepsilon}^{*}}(V_{y}-V_{\varepsilon}^{*})M_{F^{*},V_{\varepsilon}^{*}}^{\prime}=0$,
which implies
$V_{y}-V_{\varepsilon}^{*}=P_{F^{*},V_{\varepsilon}^{*}}(V_{y}-V_{\varepsilon}^{*})+(V_{y}-V_{\varepsilon}^{*})P^{\prime}_{F^{*},V_{\varepsilon}^{*}}-P_{F^{*},V_{\varepsilon}^{*}}(V_{y}-V_{\varepsilon}^{*})P^{\prime}_{F^{*},V_{\varepsilon}^{*}}$,
with $P_{F^{*},V_{\varepsilon}^{*}}=I_{T}-M_{F^{*},V_{\varepsilon}^{*}}$. From
the probability limits of Equation (FA2) for pseudo values, we have
$P_{F^{*},V_{\varepsilon}^{*}}(V_{y}-V_{\varepsilon}^{*})=(V_{y}-V_{\varepsilon}^{*})P^{\prime}_{F^{*},V_{\varepsilon}^{*}}=P_{F^{*},V_{\varepsilon}^{*}}(V_{y}-V_{\varepsilon}^{*})P^{\prime}_{F^{*},V_{\varepsilon}^{*}}=F^{*}(F^{*})^{\prime}$
(see proof of Proposition 1 (c)). Thus
$V_{y}=F^{*}(F^{*})^{\prime}+V_{\varepsilon}^{*}$, in contradiction with
Assumption A.7 (b). Thus, $n\|\hat{S}\|^{2}\geq Cn$, w.p.a. 1, for a constant
$C>0$. Moreover, using
$vech(\hat{z}^{*}_{m,n})=vech(\hat{G}^{\prime}\hat{V}_{\varepsilon}^{-1}(\sum_{i\in
I_{m}}\tilde{y}_{i}\tilde{y}_{i}^{\prime})\hat{V}_{\varepsilon}^{-1}\hat{G})$
and the conditions on $\Theta$, we get $\|vech(\hat{z}^{*}_{m,n})\|\leq
C\sum_{i\in I_{m}}\|\tilde{y}_{i}\|^{2}$. Then, from Assumptions A.2 and A.3,
$E[\|M_{\hat{\boldsymbol{X}}}\hat{\Omega}_{Z^{*}}M_{\hat{\boldsymbol{X}}}\|]\leq
C\frac{1}{n}\sum_{m=1}^{J_{n}}b_{m,n}^{2}=O(n\sum_{m=1}^{J_{n}}B_{m,n}^{2})$.
Moreover, $\sum_{m=1}^{J_{n}}B_{m,n}^{2}=o(1)$. Indeed, Assumption 2 (d)
implies $B_{m,n}\leq cn^{-\frac{\delta}{\delta+1}}$ uniformly in $m$, for any
$c>0$ and $n$ large enough, and hence
$\sum_{m=1}^{J_{n}}B_{m,n}^{2}=cn^{-\frac{\delta}{\delta+1}}\sum_{m=1}^{J_{n}}B_{m,n}\leq
c$, for any $c>0$ and $n$ large. Part (c) follows from the Lipschitz
continuity of eigenvalues for symmetric matrices.

Proof of Proposition 3: We have
$LR(k)=\frac{n}{2}\|\hat{S}\|^{2}+o_{p}(1)=vech(Z_{n}^{*})^{\prime}M_{\boldsymbol{X}}vech(Z_{n}^{*})+o_{p}(1)$
from expansion $(\ref{asy:exp:stats})$. Moreover, the kernel of matrix
$M_{\boldsymbol{X}}$ implies that
$M_{\boldsymbol{X}}vech(Z_{n}^{*})=A(F,V_{\varepsilon})z^{AD}_{n}$, where
vector $z^{AD}_{n}$ stacks the $T(T-1)/2$ above-diagonal elements of matrix
$Z_{n}$ and $A(F,V_{\varepsilon})$ is a deterministic matrix whose elements
only depend on $F,V_{\varepsilon}$. From Conditions (a) and (b) of Proposition
3, and Lemma 1, we have $z^{AD}_{n}\Rightarrow N(0,\Omega_{z})$, where the
diagonal matrix $\Omega_{z}$ is the same as if the errors were independent
normally distributed - up to replacing $q$ with ${q}+\kappa$.

Proof of Proposition 4: Let us first get the asymptotic expansion of
$\hat{V}_{y}-\tilde{V}_{\varepsilon}=\frac{1}{n}\tilde{Y}\tilde{Y}^{\prime}-\tilde{V}_{\varepsilon}$.
With the drifting DGP $Y=\mu
1_{n}^{\prime}+F\beta^{\prime}+F_{k+1}\beta_{loc}^{\prime}+\varepsilon$, and
using $\bar{\beta}=0$, $\bar{\beta}_{loc}=0$,
$\frac{1}{n}[\beta:\beta_{loc}]^{\prime}[\beta:\beta_{loc}]=I_{k+1}$ and Lemma
5 (a) in Appendix D, we get
$\hat{V}_{y}-\tilde{V}_{\varepsilon}=FF^{\prime}+\frac{1}{\sqrt{n}}\Psi_{y,loc}+R_{y}$,
where

$\Psi_{y,loc}=c_{k+1}\rho_{k+1}\rho_{k+1}^{\prime}+\frac{1}{\sqrt{n}}(\varepsilon\beta
F^{\prime}+F\beta^{\prime}\varepsilon^{\prime})+\sqrt{n}\left(\frac{1}{n}\varepsilon\varepsilon^{\prime}-\tilde{V}_{\varepsilon}\right),$
(B.3)

and
$R_{y}=\frac{1}{n}(\varepsilon\beta_{loc}F_{k+1}^{\prime}+F_{k+1}\beta_{loc}^{\prime}\varepsilon^{\prime})+[F_{k+1}F_{k+1}^{\prime}-n^{-1/2}c_{k+1}\rho_{k+1}\rho_{k+1}^{\prime}]+o_{p}(\frac{1}{\sqrt{n}})$.
Using $F_{k+1}=\sqrt{\gamma_{k+1}}\rho_{k+1}$ and
$\sqrt{n}\gamma_{k+1}=c_{k+1}+o(1)$, we get $R_{y}=o_{p}(1/\sqrt{n})$.
Subsituting the expansion for $\hat{V}_{y}-\tilde{V}_{\varepsilon}$ into
(B.1), and repeating the arguments leading to expansion
$(\ref{asy:exp:stats})$ yields expansion $(\ref{asy:exp:stats:alternative})$.
From Lemma 1, we get $vech(Z_{n,loc}^{*})\Rightarrow
N(c_{k+1}vech(\xi_{k+1}\xi_{k+1}^{\prime}),\Omega_{Z^{*}})$ as
$n\rightarrow\infty$. The result then follows from the standard result on the
distribution of idempotent quadratic forms of non-central Gaussian vectors.

Proof of Proposition 5: The proof of part (a) is in three steps. (i) The
testing problem asymptotically simplifies to the null hypothesis
$H_{0}:\lambda_{1}=...=\lambda_{df}=0$ vs. the alternative hypothesis
$H_{1}:\exists\lambda_{j}>0,j=1,...,df$. Let us define
$\boldsymbol{\lambda}_{0}=(0,...,0)^{\prime}$ for the null hypothesis and pick
a given vector
$\boldsymbol{\lambda}_{1}=(\lambda_{1},...,\lambda_{df})^{\prime}$ in the
alternative hypothesis, and consider the test of $\boldsymbol{\lambda}_{0}$
versus $\boldsymbol{\lambda}_{1}$ (simple hypothesis). By Neyman-Pearson
Lemma, the most powerful test for $\boldsymbol{\lambda}_{0}$ versus
$\boldsymbol{\lambda}_{1}$ rejects the null hypothesis when
$f(z;\lambda_{1},...,\lambda_{df})/f(z;0,...,0)$ is large, i.e., the test
function is
$\phi(z)=\boldsymbol{1}\left\\{\frac{f(z;\lambda_{1},...,\lambda_{df})}{f(z;0,...,0)}\geq
C\right\\}$ for a constant $C>0$ set to ensure the correct asymptotic size.

(ii) Let us now show that the density ratio
$\frac{f(z;\lambda_{1},...,\lambda_{df})}{f(z;0,...,0)}$ is an increasing
function of $z$. To show this, we can rely on an expansion of the density of
$\sum_{j=1}^{df}\mu_{j}\chi^{2}(1,\lambda_{j}^{2})$ in terms of central chi-
square densities (Kotz, Johnson, and Boyd (1967) Equations (144) and (151)):

$\displaystyle
f(z;\lambda_{1},...,\lambda_{df})=\sum_{k=0}^{\infty}\bar{c}_{k}(\lambda_{1},...,\lambda_{df})g(z;df+2k,0),$
(B.4)

where the coefficients
$\displaystyle\bar{c}_{k}(\lambda_{1},...,\lambda_{df})=Ae^{-\sum_{j=1}^{df}\lambda_{j}^{2}/2}E[Q(\lambda_{1},...,\lambda_{df})^{k}]/k!$
involve moments of the quadratic form $\displaystyle
Q(\lambda_{1},...,\lambda_{df})=(1/2)\sum_{j=1}^{df}\left(\nu_{j}^{1/2}X_{j}+\lambda_{j}(1-\nu_{j})^{1/2}\right)^{2}$
of the mutually independent variables $X_{j}\sim N(0,1)$,
$A=\prod_{j=1}^{df}\mu_{j}^{-1/2}$, and
$\nu_{j}=1-\frac{1}{\mu_{j}}\min_{\ell}\mu_{\ell}$. Without loss of generality
for checking the monotonicity, we have rescaled the density so that
$\min_{j}\mu_{j}=1$. Then, from (B.4), we get the ratio:
$\frac{f(z;\lambda_{1},...,\lambda_{df})}{f(z;0,...,0)}=\frac{\sum_{k=0}^{\infty}\bar{c}_{k}(\lambda_{1},...,\lambda_{df})g(z;df+2k,0)}{\sum_{k=0}^{\infty}\bar{c}_{k}(0,...,0)g(z;df+2k,0)}.$
By dividing both the numerator and the denominator by the central chi-square
density $g(z;df,0)$, we get
$\frac{f(z;\lambda_{1},...,\lambda_{df})}{f(z;0,...,0)}=e^{-\sum_{j=1}^{df}\lambda_{j}^{2}/2}\frac{\sum_{k=0}^{\infty}c_{k}(\lambda_{1},...,\lambda_{df})\psi_{k}(z)}{\sum_{k=0}^{\infty}c_{k}(0,...,0)\psi_{k}(z)}=:e^{-\sum_{j=1}^{df}\lambda_{j}^{2}/2}\Psi(z;\lambda_{1},...,\lambda_{df}),$
where $\psi_{k}(z):=$
$g(z;df+2k,0)/g(z;df,0)=\frac{\Gamma(\frac{df}{2})}{2^{k}\Gamma(\frac{df}{2}+k)}z^{k}$
is the ratio of central chi-square distributions with $df+2k$ and $df$ degrees
of freedom, and
$c_{k}(\lambda_{1},...,\lambda_{df})=E[Q(\lambda_{1},...,\lambda_{df})^{k}]/k!$.
We use the short notation
$c_{k}(\lambda):=c_{k}(\lambda_{1},...,\lambda_{df})$ and
$c_{k}(0):=c_{k}(0,...,0)$. The factor $e^{-\sum_{j=1}^{df}\lambda_{j}^{2}/2}$
does not impact on the monotonicity of the density ratio. We take the
derivative of $\Psi(z;\lambda_{1},...,\lambda_{df})$ with respect to argument
$z$ and get
$\partial_{z}\Psi(z;\lambda_{1},...,\lambda_{df})=\frac{\left(\sum_{k=1}^{\infty}c_{k}(\lambda)\psi_{k}^{\prime}(z)\right)\left(1+\sum_{k=1}^{\infty}c_{k}(0)\psi_{k}(z)\right)}{\left(\sum_{k=0}^{\infty}c_{k}(0)\psi_{k}(z)\right)^{2}}-\frac{\left(1+\sum_{k=1}^{\infty}c_{k}(\lambda)\psi_{k}(z)\right)\left(\sum_{k=1}^{\infty}c_{k}(0)\psi_{k}^{\prime}(z)\right)}{\left(\sum_{k=0}^{\infty}c_{k}(0)\psi_{k}(z)\right)^{2}}.$
The sign is given by the difference of the numerators, which is
$\sum_{k=1}^{\infty}[c_{k}(\lambda)-c_{k}(0)]\psi_{k}^{\prime}(z)+\sum_{k,l=1,k\neq
l}^{\infty}c_{k}(\lambda)c_{l}(0)[\psi_{k}^{\prime}(z)\psi_{l}(z)-\psi_{k}(z)\psi_{l}^{\prime}(z)]=\sum_{k=1}^{\infty}[c_{k}(\lambda)-c_{k}(0)]\psi_{k}^{\prime}(z)+\sum_{k,l=1,k>l}^{\infty}[c_{k}(\lambda)c_{l}(0)-c_{l}(\lambda)c_{k}(0)][\psi_{k}^{\prime}(z)\psi_{l}(z)-\psi_{k}(z)\psi_{l}^{\prime}(z)].$
We use
$\psi_{k}^{\prime}(z)=\frac{\Gamma(\frac{d}{2})k}{2^{k}\Gamma(\frac{d}{2}+k)}z^{k-1}$
and
$\psi_{k}^{\prime}(z)\psi_{l}(z)-\psi_{k}(z)\psi_{l}^{\prime}(z)=(k-l)\frac{\Gamma(\frac{d}{2})^{2}}{2^{k+l}\Gamma(\frac{d}{2}+k)\Gamma(\frac{d}{2}+l)}z^{k+l-1}$
for $k>l$ and $z\geq 0$. The difference of the numerators in the derivative of
the density ratio becomes:
$\frac{1}{2}\frac{\Gamma(\frac{d}{2})}{\Gamma(\frac{d}{2}+1)}[c_{1}(\lambda)-c_{1}(0)]+\frac{1}{2^{2}}\frac{2\Gamma(\frac{d}{2})}{\Gamma(\frac{d}{2}+2)}[c_{2}(\lambda)-c_{2}(0)]z+\sum_{m=3}^{\infty}\frac{1}{2^{m}}\left(m\frac{\Gamma(\frac{d}{2})}{\Gamma(\frac{d}{2}+m)}[c_{m}(\lambda)-c_{m}(0)]\right.$
$\left.+\sum_{k>l\geq
1,k+l=m}\frac{(k-l)\Gamma(\frac{d}{2})^{2}}{\Gamma(\frac{d}{2}+k)\Gamma(\frac{d}{2}+l)}[c_{k}(\lambda)c_{l}(0)-c_{l}(\lambda)c_{k}(0)]\right)z^{m-1}=\sum_{m=1}^{\infty}\frac{1}{2^{m}}\kappa_{m}z^{m-1},$
with $\kappa_{m}:=\sum_{k>l\geq
0,k+l=m}(k-l)\frac{\Gamma(\frac{d}{2})^{2}}{\Gamma(\frac{d}{2}+k)\Gamma(\frac{d}{2}+l)}$
$[c_{k}(\lambda)c_{l}(0)-c_{l}(\lambda)c_{k}(0)]$. A direct calculation shows
that $\kappa_{1},\kappa_{2}\geq 0$. Hence, a sufficient condition for
monotonicity of the density ratio is $\kappa_{m}\geq 0$, for all $m\geq 3$,
i.e., Inequalities (7). Thus, the test rejects for large values of the
argument, i.e., $\phi(z)=\boldsymbol{1}\\{z\geq\bar{C}\\}$, where the constant
$\bar{C}$ is determined by fixing the asymptotic size under the null
hypothesis.

(iii) Since the test function $\phi$ does not depend on
$\boldsymbol{\lambda}_{1}$, it is AUMPI in the class of hypothesis tests based
on the LR statistic (or the squared norm statistic). It yields part (a).

Let us now turn to the proof of part (b). From the definition of the
$\kappa_{m}$ coefficients written as $\kappa_{m}=\sum_{j>l\geq
0,j+l=m}\frac{(j-l)\Gamma(\frac{df}{2})^{2}}{\Gamma(\frac{df}{2}+j)\Gamma(\frac{df}{2}+l)}c_{j}(0)c_{l}(0)[\frac{c_{j}(\lambda)}{c_{j}(0)}-\frac{c_{l}(\lambda)}{c_{l}(0)}]$,
it is sufficient to get $\kappa_{m}\geq 0$, for all $m$, that sequence
$\frac{c_{j}(\lambda)}{c_{j}(0)}$, for $j=0,1,...$, is increasing. To prove
that, we link the coefficients $c_{j}(\lambda)$ to the complete exponential
Bell’s polynomials (Bell (1934)) and establish the following recurrence.

###### Lemma 2

We have
$c_{l+1}(\lambda)=\frac{1}{l+1}\sum_{i=0}^{l}\left(\frac{1}{2}\sum_{j=1}^{df}\nu_{j}^{i}\left[\nu_{j}+(i+1)(1-\nu_{j})\lambda_{j}^{2}\right]\right)c_{l-i}(\lambda)$,
for $l\geq 0$.

We use
$\frac{c_{l}(\lambda)}{c_{l}(0)}=\frac{\tilde{c}_{l}(\lambda)}{\gamma_{l}}$,
where we obtain the sequences $\gamma_{l}:=c_{l}(0)\nu_{df}^{-l}$ and
$\tilde{c}_{l}(\lambda):=c_{l}(\lambda)\nu_{df}^{-l}$ by standardization with
$\nu_{df}^{-l}$. From Lemma 2, we have
$\gamma_{l+1}=\frac{1}{l+1}\sum_{i=0}^{l}\frac{1}{2}\left(1+\sum_{j=2}^{df-1}\rho_{j}^{i+1}\right)\gamma_{l-i}$
with $\gamma_{0}=1$, and
$\tilde{c}_{l+1}(\lambda)=\frac{1}{l+1}\sum_{i=0}^{l}\left(\frac{1}{2}\sum_{j=1}^{df}\rho_{j}^{i}\left[\rho_{j}+\frac{i+1}{\nu_{df}}(1-\nu_{j})\lambda_{j}^{2}\right]\right)\tilde{c}_{l-i}(\lambda)$
with $\tilde{c}_{0}(\lambda)=1$ (note that $\rho_{1}=0$ and $\rho_{df}=1$). To
prove that sequence$\frac{\tilde{c}_{l}(\lambda)}{\gamma_{l}}$ is increasing,
the next lemma provides a sufficient condition from "separation" of the
coefficients that define the recursive relations.

###### Lemma 3

Let $(a_{i})$ be a real sequence, and let
$b_{i}=\frac{1}{2}\left(1+\sum_{j=2}^{df-1}\rho_{j}^{i}\right)$, for $i\geq
1$, where $0\leq\rho_{j}\leq 1$. Let sequences $(g_{l})$ and $(c_{l})$ be
defined recursively by
$g_{l+1}=\frac{1}{l}(b_{1}g_{l}+b_{2}g_{l-1}+...+b_{l})$ and
$c_{l+1}=\frac{1}{l}(a_{1}c_{l}+a_{2}c_{l-1}+...+a_{l})$, with
$g_{1}=c_{1}=1$. Suppose that $a_{i}\geq\max\\{\frac{df-1}{2},1\\}$, for all
$i$ (separation condition). Then, sequence $(\frac{c_{l}}{g_{l}})$ is
increasing.

We apply Lemma 3 to sequences $\tilde{c}_{l}(\lambda)$ and $\gamma_{l}$. We
detail the case $df\geq 3$ (for $df=2$ the analysis is simpler). The
separation condition
$\frac{1}{2}\sum_{j=1}^{df}\rho_{j}^{i}\left[\rho_{j}+\frac{i+1}{\nu_{df}}(1-\nu_{j})\lambda_{j}^{2}\right]\geq\frac{df-1}{2}$,
for $i=0$, yields
$\lambda_{1}^{2}+\sum_{j=2}^{df}(1-\nu_{j})\lambda_{j}^{2}\geq\nu_{df}\left(df-2-\sum_{j=2}^{df-1}\rho_{j}\right)$,
and, for $i\geq 1$, it yields
$\sum_{j=2}^{df-1}\rho_{j}^{i}(1-\nu_{j})\lambda_{j}^{2}+(1-\nu_{df})\lambda_{df}^{2}\geq\frac{\nu_{df}}{i+1}\left(df-2-\sum_{j=2}^{df-1}\rho_{j}^{i+1}\right)$.
Inequalities (8) follow.

## C Monte Carlo experiments

This appendix gives a Monte Carlo assessment of size and power and selection
procedure for the number of factors for the LR test under non-gausian errors.
Let us start with a description of the DGP we use in our simulations. In the
DGP, the betas are $\beta_{i}~{}\overset{i.i.d.}{\sim}~{}N(0,I_{k})$, with
$k=3$, and the matrix of factor values is
$F=V_{\varepsilon}^{1/2}U\Gamma^{1/2}$, where
$U=\tilde{F}(\tilde{F}^{\prime}\tilde{F})^{-1/2}$ and $vec(\tilde{F})\sim
N(0,I_{Tk})$. We generate the diagonal elements of
$V_{\varepsilon}=diag(h_{1},...,h_{T})$ through a common time-varying
component in idiosyncratic volatilities (Renault, Van Der Heijden and Werker
(2023)) via the ARCH $h_{t}=0.6+0.5h_{t-1}z_{t-1}^{2}$, with $z_{t}\sim
IIN(0,1)$. This common component induces a deviation from spherical errors.
The diagonal matrix $\Gamma=Tdiag(3,2,n^{-\bar{\kappa}})$ yields
$\frac{1}{T}F^{\prime}V_{\varepsilon}^{-1}F=diag(3,2,n^{-\bar{\kappa}})$,
i.e., the "signal-to-noise" ratios equal $3$, $2$ and $n^{-\bar{\kappa}}$ for
the three factors. We take $\bar{\kappa}=\infty$ to study the size of $LR(2)$.
To study the power of $LR(2)$, we take $\bar{\kappa}=0$ to get a global
alternative and $\bar{\kappa}=1/2$ to get a local alternative (weak factor).
We generate the idiosyncratic errors by
$\varepsilon_{i,t}=h_{t}^{1/2}h_{i,t}^{1/2}z_{i,t}$, where
$h_{i,t}=c_{i}+\alpha_{i}h_{i,t-1}z_{i,t-1}^{2}$, with $z_{i,t}\sim IIN(0,1)$
mutually independent of $z_{t}$. We use the constraint
$c_{i}=\sigma_{ii}(1-\alpha_{i})$ with uniform draws for the idiosyncratic
variances
$V[\varepsilon_{i,t}]=\sigma_{ii}~{}\overset{i.i.d.}{\sim}~{}U[1,4]$, so that
$V[\varepsilon_{i,t}/h_{t}^{1/2}]=\frac{c_{i}}{1-\alpha_{i}}=\sigma_{ii}$.
Such a setting allows for cross-sectional heterogeneity in the variances of
the scaled $\varepsilon_{i,t}/h_{t}^{1/2}$. The ARCH parameters are uniform
draws $\alpha_{i}~{}\overset{i.i.d.}{\sim}~{}U[0.2,0.5]$ with an upper
boundary of the interval ensuring existence of fourth-order moments. We
generate $5,000$ panels of returns of size $n\times T$ for each of the $100$
draws of the $T\times k$ factor matrix $F$ and common ARCH process $h_{t}$,
$t=1,...,T$, in order to keep the factor values constant within repetitions,
but also to study the potential heterogeneity of size and power results across
different factor paths. The factor betas $\beta_{i}$, idiosyncratic variances
$\sigma_{ii}$, and individual ARCH parameters $\alpha_{i}$ are the same across
all repetitions in all designs of the section. We use three different cross-
sectional sizes $n=500,1000,5000$, and three values of time-series dimension
$T=6,12,24$. The variance matrix $\hat{\Omega}_{\bar{Z}^{*}}$ is computed
using the parametric structure of Lemma 8. We get the $T-1$ estimated
parameters by least squares, as detailed in OA Section E.5.3 i). The
$p$-values are computed over $5,000$ draws.

We provide the size and power results in % in Table 1. Size of $LR(2)$ is
close to its nominal level $5\%$, with size distortions smaller than $1\%$,
except for the case $T=24$ and $n=500$. The impact of the factor values on
size is small for $T$ above $6$. The labels global power and local power refer
to $\bar{\kappa}=0$ and $\bar{\kappa}=1/2$, and power computation is not size
adjusted. The global power is equal to 100%, while the local power ranges from
80% to 85% for $T=6$, and is equal to 100% for $T=12$ and $T=24$. The
approximate constancy of local power w.r.t. $n$, for large $n$, is coherent
with theory implying convergence to asymptotic local power. In the last panel
of Table 1, we provide the average of the estimated number $\hat{k}_{LR}$ of
factors, obtained by sequential testing with $LR(k)$, for $k=0,\dots,k_{max}$,
with $k_{max}=2,7,17$ for $T=6,12,24$ (see Table 3 of OA). We follow the
procedure described in Section 6.2, with size $\alpha_{n}=10/n$. If we reject
for all $k=0,\dots,k_{max}$, then the estimated number of factors is set to
$\hat{k}_{LR}=k_{max}+1$. For all sample sizes $T=6,12,24$, the average
estimated number of factors is very close to the true number 2. We can
conclude that our selection procedure for the number of factors works well in
our simulations.

Figure 1: The upper panel displays the p-values for the statistic $LR(k)$ for
the subperiods from January 1963 to December 2021, stopping at the smallest
$k$ such that $H_{0}(k)$ is not rejected at level $\alpha_{n}=10/n_{max}$. If
no such $k$ is found then p-values are displayed up to $k_{max}$. We use
rolling windows of $T=20$ months moving forward by $12$ months each time. The
first bar of p-values covers the whole 20 months. Other bars cover the last 12
months of the 20 months subperiod. We flag bear market phases with grey shaded
vertical bars. The five lower panels display
$\left(\overline{\hat{V}}_{y}\right)^{1/2}$ for total cross-sectional
volatility, $\left(\overline{\hat{F}^{\prime}\hat{F}}\right)^{1/2}$ for
systematic volatility, $\left(\overline{\hat{V}}_{\varepsilon}\right)^{1/2}$
for idiosyncratic volatility, as well as $\hat{R}^{2}$ and $\hat{R}^{2}$ under
a single-factor model.

Figure 2: The upper and lower panels display the p-values for the RS and KP
statistics for the subperiods from January 1963 to December 2021, for the rank
test of the null hypothesis $H_{0,sp}(r)$ that $F^{O}$ has rank $r$ against
the alternative hypothesis of rank larger than $r$, for any integer $r\leq
k-1$. The empirical matrix $\hat{F}^{O}$ is computed with the time-varying
portfolio weights of the Fama-French five-factor model plus momentum. We stop
at the smallest $r$ such that $H_{0,sp}(r)$ is not rejected at level
$\alpha_{n}=10/n$. If no such $r$ is found then p-values are displayed up to
$k-1$. The red horizontal segments give $\hat{k}-1$, i.e., the estimated
number of latent factors obtained from Figure 1 minus $1$. We flag bear market
phases with grey shaded vertical bars, and use the same rolling windows as in
Figure 1.

| Size (%) | Global Power (%) | Local Power (%) | $\hat{k}_{LR}$
---|---|---|---|---
$T$ | 6 | 12 | 24 | 6 | 12 | 24 | 6 | 12 | 24 | 6 | 12 | 24
$n=500$ | 6.0 | 5.2 | 6.7 | 100 | 100 | 100 | 80 | 100 | 100 | 2.0 | 2.0 | 2.1
| (2.8) | (0.3) | (0.4) | (0.1) | (0.0) | (0.0) | (20.5) | (0.0) | (0.0) | (0.1) | (0.1) | (0.2)
$n=1000$ | 5.6 | 4.9 | 5.5 | 100 | 100 | 100 | 81 | 100 | 100 | 2.0 | 2.0 | 2.0
| (2.3) | (0.3) | (0.3) | (0.0) | (0.0) | (0.0) | (21.1) | (0.0) | (0.0) | (0.0) | (0.0) | (0.1)
$n=5000$ | 5.3 | 5.0 | 4.9 | 100 | 100 | 100 | 85 | 100 | 100 | 2.0 | 2.0 | 2.0
| (0.9) | (0.3) | (0.3) | (0.0) | (0.0) | (0.0) | (20.4) | (0.0) | (0.0) | (0.0) | (0.0) | (0.1)
Table 1: For each sample size combination $(n,T)$, we provide the average size
and power in % for the statistic $LR(2)$ (first three panels), and the average
of the estimated number $\hat{k}_{LR}$ of factors obtained by sequential
testing (last panel). Nominal size is $5\%$ for the first three panels, and
$\alpha_{n}=10/n$ for the last panel. Global power refers to the global
alternative $\bar{\kappa}=0$, and local power refers to the local alternative
$\bar{\kappa}=0.5$. In parentheses, we report the standard deviations for
size, power, and $\hat{k}_{LR}$ across $100$ different draws of the factor
path.

ONLINE APPENDIX
Latent Factor Analysis in Short Panels
Alain-Philippe Fortin, Patrick Gagliardini, and Olivier Scaillet