Interactive Effects Panel Data Models with General Factors and Regressors
††Previous versions of the paper were presented at seminars at Aarhus
University, Lund University and Michigan State University. Westerlund would
like to thank seminar participants, and in particular Richard Baillie,
Nicholas Brown, David Edgerton, Yousef Kaddoura, Yana Petrova, Simon Reese,
Tim Vogelsang, Jeffrey Wooldridge, and Morten Ørregaard Nielsen. Westerlund
thanks the Knut and Alice Wallenberg Foundation for financial support through
a Wallenberg Academy Fellowship. Peng and Westerlund also acknowledge the
Australian Research Council Discovery Grants Program for its financial support
under grant number DP210100476. Correspondence: Department of Econometrics and
Business Statistics, Monash University, Australia. Email<EMAIL_ADDRESS>
Bin Peng⋆, Liangju Su‡, Joakim Westerlund♯ and Yanrong Yang†

⋆Monash University, ‡Tsinghua University, ♯Lund University, †Australian
National University

This paper considers a model with general regressors and unobservable factors.
An estimator based on iterated principal components is proposed, which is
shown to be not only asymptotically normal and oracle efficient, but under
certain conditions also free of the otherwise so common asymptotic incidental
parameters bias. Interestingly, the conditions required to achieve
unbiasedness become weaker the stronger the trends in the factors, and if the
trending is strong enough unbiasedness comes at no cost at all. In particular,
the approach does not require any knowledge of how many factors there are, or
whether they are deterministic or stochastic. The order of integration of the
factors is also treated as unknown, as is the order of integration of the
regressors, which means that there is no need to pre-test for unit roots, or
to decide on which deterministic terms to include in the model.

Keywords: Panel data, non-stationarity, principal components, interactive
effects.

JEL classification: C14, C23, C51

##  1 Introduction

The use of panel data in regression analysis has attracted considerable
attention in the empirical literature in economics and elsewhere. A major
reason for this is the ability to deal with the presence of unobserved
heterogeneity, and the problem that this causes when said heterogeneity is
correlated with the regressors. The simplest approach is to assume that the
unobserved heterogeneity is made up of additive individual- and time-specific
constants, or “fixed effects”. Such effects can add hugely to the fit of the
model, and are typically found to be much more important than the regressors
predicted by economic theory. However, in many applications fixed effects are
unlikely to be enough to capture the unobserved heterogeneity (see, for
example, Lemmon et al., 2008, and DeAngelo and Roll, 2015), and this has in
turn spurred much interest in so-called “interactive effects” models, in which
the individual and time effects enter in an multiplicative way. One of the
most common approaches to such interactive effects models by far is the
principal components (PC) approach of Bai (2009). In fact, it is so common
that it has given rise to a separate strand of literature (see Ando and Bai,
2017, Bai and Li, 2014, Li et al., 2016, Moon and Weidner, 2015, to mention a
few).111While popular, PC is not the only approach that can be used to
estimate interactive effects models. One alternative is the so-called “common
correlated effects” (CCE) estimator of Pesaran (2006). The problem with this
approach is that it requires that the number of time effects, or factors, is
bounded by the number of observables, and that the regressors load on the same
factors as the dependent variable, which is not necessary in PC. There is also
the generalized method of moments (GMM) approach of Ahn et al. (2013).
However, this approach supposes that the number of time periods is “small”,
and is not suitable for the type of large panel data sets that we have in
mind. The present paper aims to contribute to this strand, and it does so in
at least three ways.

The first contribution of the paper is to consider a very general data
generating process (DGP) that includes most of the specifications considered
previously in the literature as special cases. The only requirement is that
suitably normalized sample second moment matrices of the factors and
regressors have positive definite limits. This is noteworthy because the
existing literature is almost exclusively based on the assumption that both
the factors and regressors are stationary. The only exceptions known to us are
Bai et al. (2009), Dong et al. (2021) and Huang et al. (2021), but they assume
instead that either the regressors or the factors and the regressors are pure
unit root processes, which is also not realistic. Indeed, regressors and
factors of different order of magnitude are likely to be the rule rather than
the exception, especially in economic and financial data, due to differences
in persistence over time.

The unrestricted DGP is important in itself but also because it can be
accommodated without requiring any knowledge thereof. Hence, not only do we
treat the factors and their number as unknown, but we also do not require any
knowledge of the order of magnitude of both factors and regressors. An
important implication of this is that there is no need to distinguish between
deterministic and stochastic factors, or stationary and non-stationary
factors. In the existing literature, deterministic factors are often treated
as known, and are projected out prior to the application of PC (see, for
example, Moon and Weidner, 2015). The problem here is that there is typically
great uncertainty over which deterministic terms to include, which raises the
issue of model misspecification. The fact that in the present paper
deterministic terms are treated as additional factors means that the problem
of deciding on which terms to include does not arise. Similarly, while the
regressors can be tested for unit roots, and the estimation can be made
conditional on the test outcome, this raises the issue of pre-testing bias. In
the present paper we do not require any knowledge about the order of
integration of the regressors, which means that there is no need for any pre-
testing.

Equally as important as the general model formulation and its empirical appeal
is the extension of the existing econometric theory, which has not yet
ventured much outside the stationary or pure unit root environments. This is
our second major contribution. The main difficulty here is not the
unrestricted specification of the factors and regressors per se, but rather
that the order of magnitude of the factors may differ. In particular, the
problem is that the nonlinearity of the PC estimator distorts the signal
coming from the factors, just as it does in estimation of nonlinear regression
models with mixtures of integrated regressors (see, for example, Park and
Phillips, 2000). This is true if both the number and order of the factors are
known, and the problem does not become any simpler when these quantities are
treated as unknown, as they are here. An additional problem that then arises
is that existing studies on the selection of the number of factors all require
that the data are stationary (see, for example, Bai and Ng, 2002, and Ahn and
Horenstein, 2013), and it is not obvious how one should go about this when the
order of magnitude of the factor is unknown.

Intuitively, the factors whose order is largest should dominate the PC
estimator. This motivates the use of an iterative estimation procedure in
which the factors and their number are estimated in order according to their
magnitude with relatively larger factors being estimated first. We begin by
prescribing a large number of factors, and estimate the resulting model by PC.
The estimated factors only capture the most dominating factors whose order of
magnitude is largest. In spite of this, we can show that the estimator is
consistent, albeit at a relatively low rate of convergence. The rate is,
however, high enough to ensure that the number of dominating factors can be
consistently estimated using a version of the eigenvalue ratio approach of Lam
and Yao (2012), and Ahn and Horenstein (2013). We then apply PC conditional on
the first-step factor estimates, and estimate the second most dominating
factors. This procedure continues until we cannot identify any more factors,
and we can show that both the estimated factors and their number are
consistent. Because of the iterative fashion in which the factors are
estimated, we refer to the new estimation procedure as “iterative PC” (IPC),
which is shown to be asymptotically normal and “oracle efficient”.

Our third contribution is to point to a “blessing” of trending factors. The
blessing occurs if the magnitude of the factors is sufficiently large, in
which case the otherwise so common asymptotic bias of the PC approach can be
completely eliminated without imposing any additional restrictions on the
cross-sectional and time series dependencies of the regression errors. This is
noteworthy, because the conclusion made in the previous literature suggests
that in order to eliminate the asymptotic bias, the errors have to be
independent.

The reminder of the paper is organized as follows. We begin by describing the
model that we will be considering and the proposed IPC approach that we will
use to estimate it. This is done in Section 2. Section 3 presents the formal
assumptions and our main asymptotic results, whose small-sample accuracy is
evaluated using Monte Carlo simulations in Section 4. Section 5 reports the
results obtained from an empirical application to the long-run relationship
between US house prices and income. For the sake of space, we present an extra
empirical study on the returns to scale in the US banking industry in the
supplementary appendix. Section 6 concludes. All proofs are relegated to the
online appendix.

##  2 Model and estimation procedure

Consider the panel data variable $y_{i,t}$, observable for $i=1,\ldots,N$
cross-sectional units and $t=1,\ldots,T$ time periods. The model of this
variable that we will be considering is given by

$\displaystyle
y_{i,t}=\mathbf{x}_{i,t}^{\prime}\boldsymbol{\beta}^{0}+\boldsymbol{\gamma}_{i}^{0\prime}\mathbf{f}_{t}^{0}+\varepsilon_{i,t},$
(2.1)

where $\mathbf{x}_{i,t}=(x_{1,i,t},\ldots,x_{d_{x},i,t})^{\prime}$ is a
$d_{x}\times 1$ vector of regressors,
$\mathbf{f}_{t}^{0}=(f_{1,t}^{0},\ldots,f_{d_{f},t}^{0})^{\prime}$ is a
$d_{f}\times 1$ vector of unobservable common factors with
$\boldsymbol{\gamma}_{i}^{0}=(\gamma_{1,i}^{0},\ldots,\gamma_{d_{f},i}^{0})^{\prime}$
being a conformable vector of factor loadings, and $\varepsilon_{i,t}$ is an
idiosyncratic error term. Moreover, the factors are divided into groups
according to their order of magnitude. There are $G$ groups of size
$d_{1},\ldots,d_{G}$, which means that $d_{1}+\cdots+d_{G}=d_{f}$. Because the
grouping is unknown, we may without loss of generality assume that the factors
are ordered, such that the first $d_{1}$ factors have the highest order of
magnitude, the next $d_{2}$ factors have the second highest order, and so on.
Hence, if we denote by $\mathbf{f}_{g,t}^{0}$ and
$\boldsymbol{\gamma}_{g,i}^{0}$ the $d_{g}\times 1$ vectors of factors and
loadings associated with group $g$, respectively, then
$\boldsymbol{\gamma}_{i}^{0\prime}\mathbf{f}_{t}^{0}=\sum_{g=1}^{G}\boldsymbol{\gamma}_{g,i}^{0\prime}\mathbf{f}_{g,t}^{0}$,
where
$\mathbf{f}_{t}^{0}=(\mathbf{f}_{1,t}^{0\prime},\ldots,\mathbf{f}_{G,t}^{0\prime})^{\prime}=(f_{1,t}^{0},\ldots,f_{d_{f},t}^{0})^{\prime}$
and
$\boldsymbol{\gamma}_{i}^{0}=(\boldsymbol{\gamma}_{1,i}^{0\prime},\ldots,\boldsymbol{\gamma}_{G,i}^{0\prime})^{\prime}=(\gamma_{1,i}^{0},\ldots,\gamma_{d_{f},i}^{0})^{\prime}$.
If $d_{f}=0$, then $G=0$ and
$\boldsymbol{\gamma}_{i}^{0\prime}\mathbf{f}_{t}^{0}=0$.

######  Remark 2.1.

The (implicit) assumption of a common slope coefficient
$\boldsymbol{\beta}^{0}$ is standard in the literature (see, for example, Bai,
2009, and Moon and Weidner, 2015). We speculate that our results can be
extended to accommodate heterogeneous slope coefficients as in, for example,
Castagnetti et al. (2015), Castagnetti et al. (2019), Li et al. (2020), and
Huang et al. (2021).

It is useful to write (2.1) on stacked vector form. Let us therefore introduce
the $T\times 1$ vectors $\mathbf{y}_{i}=(y_{i,1},\ldots,y_{i,T})^{\prime}$ and
$\boldsymbol{\varepsilon}_{i}=(\varepsilon_{i,1},\ldots,\varepsilon_{i,T})^{\prime}$,
the $T\times d_{f}$ matrix
$\mathbf{F}^{0}=(\mathbf{f}_{1}^{0},\ldots,\mathbf{f}_{T}^{0})^{\prime}$, and
the $T\times d_{x}$ matrix
$\mathbf{X}_{i}=(\mathbf{x}_{i,1},\ldots,\mathbf{x}_{i,T})^{\prime}$.
Analogous to $\mathbf{f}_{t}^{0}$, $\mathbf{F}^{0}$ is further partitioned as
$\mathbf{F}^{0}=(\mathbf{F}_{1}^{0},\ldots,\mathbf{F}_{G}^{0})^{\prime}$,
where the $T\times d_{g}$ matrix
$\mathbf{F}^{0}_{g}=(\mathbf{f}_{g,1}^{0},\ldots,\mathbf{f}_{g,T}^{0})^{\prime}$
contains the stacked factor observations for group $g$. In this notation,
(2.1) can be written as

$\displaystyle\mathbf{y}_{i}=\mathbf{X}_{i}\boldsymbol{\beta}^{0}+\mathbf{F}^{0}\boldsymbol{\gamma}_{i}^{0}+\boldsymbol{\varepsilon}_{i}.$
(2.2)

The goal of this paper is to infer $\boldsymbol{\beta}^{0}$. As mentioned in
Section 1, however, because of the generality of the model being considered,
the main difficulty in the estimation process is how to control for
$\mathbf{F}^{0}$. Our proposed estimation procedure consists of three steps.
We first initialize the estimation procedure by applying the PC estimator of
Bai (2009). However, because the first group of factors dominates all the
other groups in terms of order of magnitude, the first-step PC factor
estimator will only be estimating (the space spanned by) $\mathbf{F}_{1}^{0}$.
The second step of the procedure therefore involves iteratively applying PC
conditional on previous factor estimates to estimate all subsequent groups of
factors; hence, the “I” in IPC. In the third and final step, we estimate
$\boldsymbol{\beta}^{0}$ conditional on the second-step IPC estimator of
$\mathbf{F}^{0}$ and the first-step PC estimator of $\boldsymbol{\beta}^{0}$.

######  Step 1 (Initial estimation).

The (concentrated) ordinary least squares (OLS) objective function that we
consider is the same as in Bai (2009). It is given by

$\displaystyle\mathrm{SSR}(\boldsymbol{\beta},\mathbf{F})=\sum_{i=1}^{N}(\mathbf{y}_{i}-\mathbf{X}_{i}\boldsymbol{\beta})^{\prime}\mathbf{M}_{F}(\mathbf{y}_{i}-\mathbf{X}_{i}\boldsymbol{\beta}),$
(2.3)

where $\mathbf{F}$ is a $T\times d_{max}$ matrix satisfying
$T^{-\delta}\mathbf{F}^{\prime}\mathbf{F}=\mathbf{I}_{d_{max}}$,
$d_{max}\in\mathbb{N}\geq d_{f}$ and $\delta\in[0,\infty)$ are user-specified
numbers, and
$\mathbf{M}_{A}=\mathbf{I}_{T}-\mathbf{A}(\mathbf{A}^{\prime}\mathbf{A})^{-1}\mathbf{A}^{\prime}=\mathbf{I}_{T}-\mathbf{P}_{A}$
for any $T$-rowed full column rank matrix $\mathbf{A}$. As we explain in
Remark 2.2, the IPC estimator of $\boldsymbol{\beta}^{0}$ is invariant to the
choice of $\delta$ and the need to select $d_{max}$ is standard. The initial
estimator that we consider is the minimizer of
$\mathrm{SSR}(\boldsymbol{\beta},\mathbf{F})$:

$\displaystyle(\widehat{\boldsymbol{\beta}}_{0},\widehat{\mathbf{F}}_{0})=\operatorname*{\arg\\!\min}_{(\boldsymbol{\beta},\mathbf{F})\in\mathbb{D}}\,\mathrm{SSR}(\boldsymbol{\beta},\mathbf{F}),$
(2.4)

where $\mathbb{D}=\mathbb{R}^{d_{x}}\times\mathbb{D}_{F}$ and
$\mathbb{D}_{F}=\\{\mathbf{F}:T^{-\delta}\mathbf{F}^{\prime}\mathbf{F}=\mathbf{I}_{d_{max}}\\}$.
It is useful to note that $\widehat{\boldsymbol{\beta}}_{0}$ satisfies
$\widehat{\boldsymbol{\beta}}_{0}=\widehat{\boldsymbol{\beta}}(\widehat{\mathbf{F}}_{0})$,
where

$\displaystyle\widehat{\boldsymbol{\beta}}(\mathbf{F})=\left(\sum_{i=1}^{N}\mathbf{X}_{i}^{\prime}\mathbf{M}_{F}\mathbf{X}_{i}\right)^{-1}\sum_{i=1}^{N}\mathbf{X}_{i}^{\prime}\mathbf{M}_{F}\mathbf{y}_{i}.$
(2.5)

######  Step 2 (Iterative estimation of factors).

As already pointed out, the factors are estimated in order according to
magnitude. Therefore, $\widehat{\mathbf{F}}_{0}$ is estimating
$\mathbf{F}_{1}^{0}$. Since $d_{1}\leq d_{f}\leq d_{max}$, in general the
dimension of $\widehat{\mathbf{F}}_{0}$ will be larger than that of
$\mathbf{F}_{1}^{0}$. We therefore begin this step of the estimation procedure
by estimating $d_{1}$, and for this purpose we employ a version of the ratio
of eigenvalue-based estimator considered by, for example, Lam and Yao (2012),
and Ahn and Horenstein (2013), which is given by

$\displaystyle\widehat{d}_{1}=\operatorname*{\arg\\!\min}_{0\leq d\leq
d_{max}}\left\\{\frac{\widehat{\lambda}_{1,d+1}}{\widehat{\lambda}_{1,d}}\cdot\mathbb{I}\left(\frac{\widehat{\lambda}_{1,d}}{\widehat{\lambda}_{1,0}}\geq\tau_{N}\right)+\mathbb{I}\left(\frac{\widehat{\lambda}_{1,d}}{\widehat{\lambda}_{1,0}}<\tau_{N}\right)\right\\},$
(2.6)

where $\mathbb{I}(A)$ is the indicator function for the event $A$ taking the
value one if $A$ is true and zero otherwise,
$\tau_{N}=1/\ln\,(\max\\{\widehat{\lambda}_{1,0},N\\})$,
$\widehat{\lambda}_{1,0}=N^{-1}\sum_{i=1}^{N}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0})^{\prime}\mathbf{M}_{\widehat{F}_{0}}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0})$
and $\widehat{\lambda}_{1,1}\geq\cdots\geq\widehat{\lambda}_{1,d_{max}}$ are
the $d_{max}$ largest eigenvalues of the following $T\times T$ matrix:

$\displaystyle\widehat{\boldsymbol{\Sigma}}_{1}=\frac{1}{N}\sum_{i=1}^{N}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0})(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0})^{\prime}.$
(2.7)

The threshold $\tau_{N}$, the “mock” eigenvalue $\widehat{\lambda}_{1,0}$, and
the indicator function are there to ensure that the estimator is consistent.
The need for these will be explained later. Given $\widehat{d}_{1}$, we update
the estimate of $\mathbf{F}_{1}^{0}$ by setting $\widehat{\mathbf{F}}_{1}$
equal to $T^{\delta/2}$ times the eigenvectors associated with
$\widehat{\lambda}_{1,1},\ldots,\widehat{\lambda}_{1,\widehat{d}_{1}}$, and
estimate $\boldsymbol{\gamma}_{1,i}^{0}$ by
$\widehat{\boldsymbol{\gamma}}_{1,i}=T^{-\delta}\widehat{\mathbf{F}}_{1}^{\prime}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0})$.

The estimation of $\mathbf{F}_{2}^{0},\ldots,\mathbf{F}_{G}^{0}$ is analogous
to that of $\mathbf{F}_{1}^{0}$. The main difference is that we have to
condition on all previous estimates. Let us therefore use
$\widehat{\mathbf{F}}_{-g}=(\widehat{\mathbf{F}}_{1},\ldots,\widehat{\mathbf{F}}_{g-1})$
and
$\widehat{\boldsymbol{\gamma}}_{-g,i}=(\widehat{\boldsymbol{\gamma}}_{1,i}^{\prime},\ldots,\widehat{\boldsymbol{\gamma}}_{g-1,i}^{\prime})^{\prime}$
to denote the matrices containing the previously estimated factors and
loadings, respectively, when estimating group $g>1$. The estimator of $d_{g}$
is given quite naturally by

$\displaystyle\widehat{d}_{g}=\operatorname*{\arg\\!\min}_{0\leq d\leq
d_{max}}\left\\{\frac{\widehat{\lambda}_{g,d+1}}{\widehat{\lambda}_{g,d}}\cdot\mathbb{I}\left(\frac{\widehat{\lambda}_{g,d}}{\widehat{\lambda}_{g,0}}\geq\tau_{N}\right)+\mathbb{I}\left(\frac{\widehat{\lambda}_{g,d}}{\widehat{\lambda}_{g,0}}<\tau_{N}\right)\right\\},$
(2.8)

where
$\widehat{\lambda}_{g,0}=N^{-1}\sum_{i=1}^{N}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0})^{\prime}\mathbf{M}_{\widehat{F}_{-g}}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0})$
and $\widehat{\lambda}_{g,1}\geq\cdots\geq\widehat{\lambda}_{g,d_{max}}$ are
the $d_{max}$ largest eigenvalues of

$\displaystyle\widehat{\boldsymbol{\Sigma}}_{g}=\frac{1}{N}\sum_{i=1}^{N}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0}-\widehat{\mathbf{F}}_{-g}\widehat{\boldsymbol{\gamma}}_{-g,i})(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0}-\widehat{\mathbf{F}}_{-g}\widehat{\boldsymbol{\gamma}}_{-g,i})^{\prime}.$
(2.9)

The resulting estimator $\widehat{\mathbf{F}}_{g}$ of $\mathbf{F}_{g}^{0}$ is
given by the eigenvectors associated with
$\widehat{\lambda}_{g,1},\ldots,\widehat{\lambda}_{g,\widehat{d}_{g}}$ and
$\widehat{\boldsymbol{\gamma}}_{g,i}=T^{-\delta}\widehat{\mathbf{F}}_{g}^{\prime}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}}_{0}-\widehat{\mathbf{F}}_{-g}\widehat{\boldsymbol{\gamma}}_{-g,i})$.
New groups of factors are estimated until $\widehat{d}_{g}=0$. At this point,
we set $\widehat{G}=g-1$ and define
$\widehat{\mathbf{F}}=(\widehat{\mathbf{F}}_{1},\ldots,\widehat{\mathbf{F}}_{\widehat{G}})$.
This is the IPC estimator of $\mathbf{F}^{0}$.

######  Step 3 (Estimation of $\boldsymbol{}\beta^{0}$).

Given $\widehat{\mathbf{F}}$, we compute
$\widehat{\boldsymbol{\beta}}_{1}=\widehat{\boldsymbol{\beta}}(\widehat{\mathbf{F}})$
using (2.5). The IPC-based estimator of $\boldsymbol{\beta}^{0}$ is given by

$\displaystyle\widehat{\boldsymbol{\beta}}=\widehat{\boldsymbol{\beta}}_{0}+\left(\sum_{i=1}^{N}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i}\right)^{-1}\sum_{i=1}^{N}\mathbf{X}_{i}^{\prime}\mathbf{M}_{\widehat{F}}\mathbf{X}_{i}(\widehat{\boldsymbol{\beta}}_{1}-\widehat{\boldsymbol{\beta}}_{0}),$
(2.10)

where
$\widehat{\mathbf{Z}}_{i}=\mathbf{M}_{\widehat{F}}\mathbf{X}_{i}-\sum_{j=1}^{N}\mathbf{M}_{\widehat{F}}\mathbf{X}_{j}\widehat{a}_{ij}$
with
$\widehat{a}_{ij}=\widehat{\boldsymbol{\gamma}}_{i}^{\prime}(\widehat{\boldsymbol{\Gamma}}^{\prime}\widehat{\boldsymbol{\Gamma}})^{-1}\widehat{\boldsymbol{\gamma}}_{j}$,
$\widehat{\boldsymbol{\gamma}}_{i}=(\widehat{\boldsymbol{\gamma}}_{1,i}^{\prime},\ldots,\widehat{\boldsymbol{\gamma}}_{\widehat{G},i}^{\prime})^{\prime}$
and
$\widehat{\boldsymbol{\Gamma}}=(\widehat{\boldsymbol{\gamma}}_{1},\ldots,\widehat{\boldsymbol{\gamma}}_{N})^{\prime}$.

######  Remark 2.2.

The reason for why in Steps 1 and 2 $\delta$ can be set arbitrarily is that
the only use of this parameter is to normalize $\widehat{\mathbf{F}}$ and
$T^{-\delta/2}\|\widehat{\mathbf{F}}_{0}\|=\sqrt{d_{max}}$ for any
$\delta\in[0,\infty)$. This result is nice from both applied and theoretical
points of view. In the bulk of the previous literature, the appropriate value
of $\delta$ to use depends on whether $\mathbf{F}^{0}$ is stationary or unit
root non-stationary (see, for example, Bai, 2004). The assumed knowledge of
$\delta$ is therefore tantamount to assuming that the order of integration of
$\mathbf{F}^{0}$ is known, which is again not a requirement here. The need to
specify a maximum $d_{max}$ for the number of factors is standard in the
literature (see, for example, Bai and Ng, 2002). In Sections 4 and 5, we
explain how $d_{max}$ is set in the Monte Carlo and empirical studies,
respectively.

######  Remark 2.3.

Many papers use information criteria to select the number of factors (see, for
example, Bai and Ng, 2002). We do not. The reason is that the eigenvalue ratio
is self-normalizing, which makes it possible to handle factors that are of
different order of magnitude. The general idea behind the eigenvalue ratio
approach in Step 2 is similar to the one used in ocular inspection of “scree
plots”; one looks for the point at which the ordered eigenvalues drop
substantially, and set the number of estimated factors accordingly. The most
natural way to mimic this decision rule is to minimize
$\widehat{\lambda}_{g,d+1}/\widehat{\lambda}_{g,d}$ over $1\leq d\leq
d_{max}$. However, this raises two issues. One is that since
$\widehat{\lambda}_{g,d+1}/\widehat{\lambda}_{g,d}$ is not defined for $d=0$,
we cannot have $d_{f}=0$, and we want to be able to entertain the possibility
that there are no factors. The use of the mock eigenvalue
$\widehat{\lambda}_{g,0}$ allows us to do just that. The other problem is that
under the conditions of this paper (laid out in Section 3) the limiting
behavior of $\widehat{\lambda}_{g,d+1}/\widehat{\lambda}_{g,d}$ when $d>d_{g}$
is unknown. Lam and Yao (2012) discuss this issue at length. They conjecture
that $\widehat{\lambda}_{g,d+1}/\widehat{\lambda}_{g,d}\asymp 1$, where
$a\asymp b$ means that $a=O_{P}(b)$ and $b=O_{P}(a)$, but there is no proof.
The challenge when studying
$\widehat{\lambda}_{g,d+1}/\widehat{\lambda}_{g,d}$ for $d>d_{g}$ is to bound
this ratio from below, as both of the numerator and denominator converge to
zero for a certain choice of $\delta$ at the same rate. The use of the
indicator function allows us to circumvent this problem. The idea is to look
at $\widehat{\lambda}_{g,d}$ only. If this eigenvalue is “small”, we take it
as a sign of $d>d_{g}$ and set
$\widehat{\lambda}_{g,d+1}/\widehat{\lambda}_{g,d}$ to one. However, because
the order of $\mathbf{f}_{t}^{0}$ is assumed to be unknown, we cannot look at
$\widehat{\lambda}_{g,d}$ directly but rather we look at
$\widehat{\lambda}_{g,d}/\widehat{\lambda}_{g,0}$, which in contrast to
$\widehat{\lambda}_{g,d}$ is self-normalizing. Barigozzi and Trapani (2021)
also look at ratios of eigenvalues as a way to determine the number of factors
when their order is unknown. However, they only consider factors that are
either trending linearly, unit root non-stationary or stationary, and the
implementation of their procedure is quite involved. We therefore do not
consider it here.

######  Remark 2.4.

Intuition suggests to take $\widehat{\boldsymbol{\beta}}_{1}$, the OLS
estimator conditional on $\widehat{\mathbf{F}}$, as the final estimator of
$\boldsymbol{\beta}^{0}$ in Step 3. Interestingly, while consistent, because
of the stepwise estimation of the factors, the asymptotic distribution of
$\widehat{\boldsymbol{\beta}}_{1}$ is generally not (mixed) normal and
nuisance parameter-free. In fact, as we show later in Lemma 3.1, even
$\widehat{\boldsymbol{\beta}}_{0}$ is consistent but again with non-normal
asymptotic distribution. In Section 4, we use Monte Carlo simulations to shed
some light on the consistency and non-normality of
$\widehat{\boldsymbol{\beta}}_{0}$ and $\widehat{\boldsymbol{\beta}}_{1}$.

##  3 Assumptions and asymptotic results

Before we state our assumptions and asymptotic results, we introduce some
notation. Specifically, if $\mathbf{A}$ is a matrix,
$\lambda_{min}(\mathbf{A})$ and $\lambda_{max}(\mathbf{A})$ signify its
smallest and largest eigenvalues, respectively, $\text{tr}\,\mathbf{A}$
signifies its trace, and
$\|\mathbf{A}\|=\sqrt{\text{tr}\,\mathbf{A}^{\prime}\mathbf{A}}$ and
$\|\mathbf{A}\|_{2}=\sqrt{\lambda_{max}(\mathbf{A})}$ signify its Frobenius
and spectral norms, respectively. We write $\mathbf{A}>0$ to signify that
$\mathbf{A}$ is positive definite. If $\mathbf{B}$ is also a matrix, then
$\mathrm{diag}(\mathbf{A},\mathbf{B})$ denotes the block-diagonal matrix that
takes $\mathbf{A}$ ($\mathbf{B}$) as the upper left (lower right) block. The
symbols $\to_{D}$, $\to_{P}$ and $MN(\cdot,\cdot)$ signify convergence in
distribution, convergence in probability and a mixed normal distribution,
respectively. We use $N,\,T\to\infty$ to indicate that the limit has been
taken while passing both $N$ and $T$ to infinity. We use w.p.a.1 to denote
with probability approaching one.

Assumption 1 is concerned with the order of magnitude of $\mathbf{f}_{t}^{0}$
and $\mathbf{x}_{i,t}$. It is therefore key. The assumption is stated in terms
of the required moment conditions rather than primitive assumptions on
$\mathbf{f}_{t}^{0}$ and $\mathbf{x}_{i,t}$. The reason is that we would like
to be agnostic about these variables.

######  Assumption 1 (Moments).

  * (a)

There exists a matrix $\boldsymbol{\Sigma}_{X}$ such that
$E\|(NT)^{-1}\sum_{i=1}^{N}\mathbf{D}_{T}\mathbf{X}_{i}^{\prime}\mathbf{M}_{F^{0}}\mathbf{X}_{i}\mathbf{D}_{T}-\boldsymbol{\Sigma}_{X}\|^{2}=o(1)$,
where
$\mathbf{D}_{T}=\mathrm{diag}(T^{-\kappa_{1}/2},\ldots,T^{-\kappa_{d_{x}}/2})$
with $0\leq\kappa_{j}<\infty$ for $j=1,\ldots,d_{x}$,
$E\|\boldsymbol{\Sigma}_{X}\|^{2}<\infty$, and
$0<\lambda_{min}(\boldsymbol{\Sigma}_{X})\leq\lambda_{max}(\boldsymbol{\Sigma}_{X})<\infty$
w.p.a.1.

  * (b)

$\|(NT)^{-1}\sum_{i=1}^{N}\mathbf{D}_{T}\mathbf{X}_{i}^{\prime}\boldsymbol{\varepsilon}_{i}\|=O_{P}(1/\min\\{\sqrt{N},\sqrt{T}\\})$
and $\|\boldsymbol{\varepsilon}\|_{2}=O_{P}(\max\\{\sqrt{N},\sqrt{T}\\})$ with
$\boldsymbol{\varepsilon}=(\boldsymbol{\varepsilon}_{1},\ldots,\boldsymbol{\varepsilon}_{N})$.

  * (c)

There exists a matrix $\boldsymbol{\Sigma}_{F^{0}}$ such that
$E\|\mathbf{C}_{T}\mathbf{F}^{0\prime}\mathbf{F}^{0}\mathbf{C}_{T}-\boldsymbol{\Sigma}_{F^{0}}\|^{2}=o(1)$,
where
$\mathbf{C}_{T}=\mathrm{diag}(T^{-\nu_{1}/2}\mathbf{I}_{d_{1}},\ldots,T^{-\nu_{G}/2}\mathbf{I}_{d_{G}})$
with $\nu_{1}>\cdots>\nu_{G}>1/2$,
$E\|\boldsymbol{\Sigma}_{F^{0}}\|^{2}<\infty$ and
$0<\lambda_{min}(\boldsymbol{\Sigma}_{F^{0}})\leq\lambda_{max}(\boldsymbol{\Sigma}_{F^{0}})<\infty$
w.p.a.1.

  * (d)

There exists a matrix $\boldsymbol{\Sigma}_{\Gamma^{0}}$ such that
$\|N^{-1}\boldsymbol{\Gamma}^{0\prime}\boldsymbol{\Gamma}^{0}-\boldsymbol{\Sigma}_{\Gamma^{0}}\|=o_{P}(1)$
and $\max_{i\geq 1}E\|\boldsymbol{\gamma}_{i}^{0}\|^{4}<\infty$, where
$\boldsymbol{\Gamma}^{0}=(\boldsymbol{\gamma}_{1}^{0},\ldots,\boldsymbol{\gamma}_{N}^{0})^{\prime}$
is $N\times d_{f}$ and
$0<\lambda_{min}(\boldsymbol{\Sigma}_{\Gamma^{0}})\leq\lambda_{max}(\boldsymbol{\Sigma}_{\Gamma^{0}})<\infty$.

Assumption 1 (a) is very general in that it imposes almost no restrictions on
the type of trending behaviour that $\mathbf{x}_{i,t}$ may have. The trending
can be deterministic but it can also be stochastic, as in the presence of unit
roots. Either way, the degree of the trending is not restricted in any way,
provided that it is finite. The only requirement is that
$0<\lambda_{min}(\boldsymbol{\Sigma}_{X})\leq\lambda_{max}(\boldsymbol{\Sigma}_{X})<\infty$
w.p.a.1, which means that the elements of $\mathbf{x}_{i,t}$ cannot be
asymptotically collinear. Note that $\boldsymbol{\Sigma}_{X}$ is not required
to be a constant matrix, as this would rule out regressors that are
stochastically integrated. This is very different when compared to the bulk of
the previous PC-based interactive effects literature in which
$\mathbf{x}_{i,t}$ is assumed to be either stationary (see, for example, Bai,
2009, and Moon and Weidner, 2015), such that
$\kappa_{1}=\cdots=\kappa_{d_{x}}=0$, or unit-root non-stationary (see Bai et
al., 2009, and Dong et al., 2021), such that
$\kappa_{1}=\cdots=\kappa_{d_{x}}=1$. Assumption 1 (a) is analogous to the
work of, for example, Dong and Linton (2018), wherein time series of different
order of magnitude are considered.

The first requirement of Assumption 1 (b) is quite mild and holds if a central
limit theorem in only one of the two panel dimensions applies to the
normalized sum of
$\mathbf{D}_{T}\mathbf{X}_{i}^{\prime}\boldsymbol{\varepsilon}_{i}$. The
second requirement is quite common in the literature, and is expected to hold
as long as $\varepsilon_{i,t}$ has zero mean, and weak serial and cross-
sectional correlation (see Moon and Weidner, 2015, for a discussion).

Assumption 1 (c) is similar to (a) in that it leaves the trending behavior of
the factors essentially unrestricted. A majority of previous PC-based works
assume that $T^{-1}\mathbf{F}^{0\prime}\mathbf{F}^{0}$ converges to positive
definite matrix (see, for example, Bai, 2009, and Moon and Weidner, 2015).
Notable exceptions include Bai (2004), Bai et al. (2009), and Choi (2017), in
which $\mathbf{f}_{t}^{0}$ is assumed to follow pure unit root process, and
Bai and Ng (2004), who allow for a mix of stationary and unit root factors.
The only study that comes close to ours in terms of the generality of the
factors is that of Westerlund (2018). However, he assumes that
$\mathbf{x}_{i,t}$ has a factor structure that loads on the same set of
factors as $y_{i,t}$, which is not required here. Also, unlike Westerlund
(2018), we do not require $\nu_{G}\geq 1$, but also allow $1/2<\nu_{G}<1$,
such that $T^{-1}\mathbf{F}_{G}^{0\prime}\mathbf{F}_{G}^{0}$ is asymptotically
singular. The signal coming from $\mathbf{F}_{G}^{0}$ is therefore even weaker
than under stationarity, and we will therefore refer to these factors as
“signal-weak”, as opposed to “weak” in usual sense (see, for example, Chudik
et al., 2011 and Uematsu and Yamagata (2020), for discussions). One example of
such signal-weak factors is when $\mathbf{f}_{G,t}^{0}$ is stationary and
sparse.

Assumption 1 (d) is standard and ensures that each factor has a nontrivial
contribution to the variance of $y_{i,t}$.

######  Assumption 2 (Identification).

$\inf_{\mathbf{F}\in\mathbb{D}_{F}}\lambda_{min}(\mathbf{B}(\mathbf{F}))\geq
c_{0}>0$ for all $N$ and $T$, where

$\displaystyle\mathbf{B}(\mathbf{F})=\frac{1}{NT}\sum_{i=1}^{N}\mathbf{D}_{T}\mathbf{Z}_{i}(\mathbf{F})^{\prime}\mathbf{Z}_{i}(\mathbf{F})\mathbf{D}_{T}$
(3.1)

with
$\mathbf{Z}_{i}(\mathbf{F})=\mathbf{M}_{F}\mathbf{X}_{i}-\sum_{j=1}^{N}\mathbf{M}_{F}\mathbf{X}_{j}a_{ij}$
and
$a_{ij}=\boldsymbol{\gamma}_{i}^{0\prime}(\boldsymbol{\Gamma}^{0\prime}\boldsymbol{\Gamma}^{0})^{-1}\boldsymbol{\gamma}_{j}^{0}$.

The requirement that the smallest eigenvalue of $\mathbf{B}(\mathbf{F})$
should be positive is equivalent to requiring that $\mathbf{B}(\mathbf{F})$ be
positive definite uniformly in $\mathbf{F}$, which is a non-collinearity
condition that rules out low-rank regressors (see Moon and Weidner, 2015, for
a discussion). It demands that the regressors in $\mathbf{X}_{i}$ have enough
variation after projecting out all variation that can be explained by
$\mathbf{F}^{0}$ and $\boldsymbol{\gamma}_{i}^{0}$. Time-invariant regressors
are therefore not allowed, which is similar to the conventional cross-section
fixed effects-only OLS condition. Assumption 2 is a high-level condition, just
like Assumption 1 (see, for example, Bai, 2009, and Moon and Weidner, 2015,
for similar high-level conditions). The reason for formulating it in this way
is again that we want to avoid making highly specific and potentially invalid
assumptions about the DGP of $\mathbf{x}_{i,t}$, $\mathbf{f}_{t}^{0}$ and
$\boldsymbol{\gamma}_{i}^{0}$. If all the regressors are stationary,
Assumption 2 reduces to Assumption A in Bai (2009).

Assumptions 1 and 2 are enough to ensure that the initial estimator
$\widehat{\boldsymbol{\beta}}_{0}$ of $\boldsymbol{\beta}^{0}$ is consistent.

######  Lemma 3.1 (Consistency of $\widehat{\boldsymbol{}\beta}_{0}$).

Under Assumptions 1 and 2, as $N,T\to\infty$,

$\displaystyle\min\\{\sqrt{N},\sqrt{T}\\}\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}_{0}-\boldsymbol{\beta}^{0})=O_{P}(1).$
(3.2)

Lemma 3.1 establishes that $\widehat{\boldsymbol{\beta}}_{0}$ is consistent
for $\boldsymbol{\beta}^{0}$ and that the rate of convergence is
$\|\mathbf{D}_{T}\|/\min\\{\sqrt{N},\sqrt{T}\\}=\max\\{T^{-\kappa_{1}/2},\ldots,T^{-\kappa_{d_{x}}/2}\\}/\min\\{\sqrt{N},\sqrt{T}\\}$.
To put this into perspective, suppose that $\mathbf{x}_{i,t}$ is stationary,
such that $\kappa_{1}=\cdots=\kappa_{d_{x}}=0$. In this case,
$\mathbf{D}_{T}=\mathbf{I}_{d_{x}}$ and the rate of convergence is given by
$1/\min\\{\sqrt{N},\sqrt{T}\\}$, which is the slowest of the regular rates in
pure time series and cross-section regressions. Still, the rate is fast enough
for the estimation of the number of factors. This brings us to Step 2 of the
estimation procedure. In order to be able to show that
$\widehat{d}_{1},\ldots,\widehat{d}_{G+1}$ and $\widehat{\mathbf{F}}$ are
consistent, however, we need more structure.

######  Assumption 3 (Errors).

  * (a)

$E(\varepsilon_{i,t})=0$ and
$E(\boldsymbol{\varepsilon}_{i}\boldsymbol{\varepsilon}_{i}^{\prime})=\boldsymbol{\Sigma}_{\varepsilon,i}$.

  * (b)

Let
$\boldsymbol{\varepsilon}_{t}=(\varepsilon_{1,t},\ldots,\varepsilon_{N,t})^{\prime}$
in this assumption only. $\\{\boldsymbol{\varepsilon}_{t}:\ t\geq 1\\}$ is
strictly stationary and $\alpha$-mixing such that $\max_{i\geq
1}E|\varepsilon_{i,1}|^{4+\mu}<\infty$ for some $\mu>0$ and the mixing
coefficient
$\alpha(t)=\sup_{A\in\mathcal{F}_{-\infty}^{0},B\in\mathcal{F}_{t}^{\infty}}|P(A)P(B)-P(AB)|$
satisfies $\sum_{t=1}^{\infty}\alpha(t)^{\mu/(2+\mu)}<\infty$, where
$\mathcal{F}_{-\infty}^{t}$ and $\mathcal{F}_{t}^{\infty}$ are the sigma-
algebras generated by $\\{\boldsymbol{\varepsilon}_{s}:s\leq 0\\}$ and
$\\{\boldsymbol{\varepsilon}_{s}:s\geq t\\}$, respectively.

  * (c)

Suppose that
$\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T}\sum_{s=1}^{T}|E\varepsilon_{i,t}\varepsilon_{j,s}|=O(NT)$
and $\sum_{i=1}^{N}\sum_{j\neq i}|\sigma_{\varepsilon,ij}|=O(N)$, where
$\sigma_{\varepsilon,ij}=E(\varepsilon_{i,t}\varepsilon_{j,t})$.

  * (d)

$\varepsilon_{i,t}$ is independent of $\boldsymbol{\gamma}_{j}^{0}$,
$\mathbf{f}_{s}^{0}$ and $\mathbf{x}_{j,s}$ for all $i$, $j$, $t$ and $s$.

Assumption 3 is similar to Assumptions C and D in Bai (2009). Assumption 3 (b)
and (c) ensure that the serial and cross-sectional dependencies of
$\varepsilon_{i,t}$ are at most weak. Assumption 3 (d) requires that the
regressors are exogenous, which rules out the presence of lagged dependent
variables in $\mathbf{x}_{i,t}$. Assumption 3 (d) can be relaxed by instead
requiring that $\varepsilon_{i,t}$ satisfies a martingale type condition
similar to Assumption 3 in Dong et al. (2021). However, since this will add to
the already lengthy derivations and heavy notation, in the present paper we
maintain Assumption 3.

If $\nu_{G}\geq 1$, Assumptions 1–3 are enough to ensure that
$\widehat{d}_{1}$ is consistent. If, however, $\nu_{G}<1$, so that some of the
factors are signal-weak, then we also need Assumption 4.

######  Assumption 4 (Signal-weak factors).

If $\nu_{G}<1$, then $T/N^{2}\to c_{1}\in[0,\infty)$.

As already pointed out, $\mathbf{F}_{1}^{0}$ dominates all other factors and
is thus easiest to estimate. It is therefore not surprising to find that
$\widehat{d}_{1}$ is consistent for $d_{1}$. Lemma 3.2 states this result
formally.

######  Lemma 3.2 (Consistency of $\widehat{d}_{1}$).

Under Assumptions 1–4, as $N,\,T\to\infty$,

$\displaystyle P(\widehat{d}_{1}=d_{1})\to 1.$ (3.3)

Hence, $\widehat{d}_{1}$ is consistent. In order to ensure that also
$\widehat{d}_{2},\ldots,\widehat{d}_{G+1}$ are consistent, the groups have to
be “distinct”. The next assumption formalizes this requirement. The assumption
is stated in terms of $\mathbf{F}_{g}^{0}$ and $\boldsymbol{\Gamma}_{g}^{0}$,
where $\mathbf{F}_{g}^{0}$ is as before and
$\boldsymbol{\Gamma}_{g}^{0}=(\boldsymbol{\gamma}_{g,1}^{0},\ldots,\boldsymbol{\gamma}_{g,N}^{0})^{\prime}$
is the $N\times d_{g}$ matrix of stacked factor loadings for group $g$.

######  Assumption 5 (Orthogonal groups).

$\max_{g\neq
h}\|\boldsymbol{\Gamma}_{g}^{0\prime}\boldsymbol{\Gamma}_{h}^{0}\|=O_{P}(N^{p})$
and $\max_{g\neq
h}\|\mathbf{F}_{g}^{0\prime}\mathbf{F}_{h}^{0}\|=O_{P}(T^{q})$, where
$g,h=1,\ldots,G$, $G>1$, $p<1$, $q<(\nu_{G}+\nu_{G-1})/2$ and $\nu_{G-1}\geq
1$.

The condition that $\nu_{G-1}\geq 1$ means that we only allow for one group of
weak factors. This can be seen as a form of normalization and is not
particularly restrictive. Let us therefore instead consider $\max_{g\neq
h}\|\mathbf{F}_{g}^{0\prime}\mathbf{F}_{h}^{0}\|=O_{P}(T^{q})$, which is less
restrictive that the exact orthogonality condition typically required in
papers on grouped factor structures (see, for example, Ando and Bai, 2017). As
is well known,
$\boldsymbol{\gamma}_{i}^{0\prime}\mathbf{f}_{t}^{0}=\boldsymbol{\gamma}_{i}^{0\prime}\mathbf{W}^{-1}\mathbf{W}\mathbf{f}_{t}^{0}$
for any positive definite matrix $\mathbf{W}$. Now set
$\mathbf{W}=(\mathbf{F}^{0\prime}\mathbf{F}^{0}\mathbf{C}_{T}^{2})^{-1/2}$.
This implies

$\mathbf{C}_{T}\mathbf{W}\mathbf{F}^{0\prime}\mathbf{F}^{0}\mathbf{W}^{\prime}\mathbf{C}_{T}=\mathbf{I}_{d_{f}}$
(3.4)

which means that the renormalized factors are exactly orthogonal, and hence
that Assumption 5 is satisfied for $\mathbf{F}^{0}\mathbf{W}^{\prime}$ with
$q=-\infty$. This is one way in which $\max_{g\neq
h}\|\mathbf{F}_{g}^{0\prime}\mathbf{F}_{h}^{0}\|=O_{P}(T^{q})$ can be
justified. Another way is through orthogonal basis functions on $[0,1]$. For
example, if $\mathbf{f}_{t}^{0}=(1,t)^{\prime}$ and
$\mathbf{C}_{T}=\mathrm{diag}(T^{-1/2},T^{-3/2})$, then
$\sqrt{T}\mathbf{C}_{T}\mathbf{f}_{t}^{0}\rightarrow\mathbf{f}^{0}(w)=(1,w)^{\prime}$
and
$\mathbf{C}_{T}\mathbf{F}^{0\prime}\mathbf{F}^{0}\mathbf{C}_{T}\rightarrow\int_{0}^{1}\mathbf{f}^{0}(w)\mathbf{f}^{0}(w)^{\prime}dw$
as $T\rightarrow\infty$ with $w\in[0,1]$. Because they are square integrable,
the elements of $\mathbf{f}^{0}(w)$ can be represented in terms of orthogonal
basis functions (see, for example, Dong and Linton, 2018). This is true not
only for the simple example given here but for most commonly used
deterministic trend terms. The loading condition, $\max_{g\neq
h}\|\boldsymbol{\Gamma}_{g}^{0\prime}\boldsymbol{\Gamma}_{h}^{0}\|=O_{P}(N^{p})$,
can be rationalized in the same way as the factor condition. Another
possibility is that
$\boldsymbol{\Gamma}_{1}^{0},\ldots,\boldsymbol{\Gamma}_{G}^{0}$ are
independent and at most one of them has non-zero mean. Independence is often
assumed and we therefore do not justify it here (see, for example, Chudik et
al., 2011, and Pesaran, 2006). In order to justify the zero mean assumption,
suppose for simplicity that $G=2$, that
$d_{1}=\boldsymbol{\gamma}_{1,i}^{0}=1$ and that
$\boldsymbol{\gamma}_{2,i}^{0}=\boldsymbol{\gamma}_{2}^{0}+\boldsymbol{\eta}_{i}$
with $E(\boldsymbol{\eta}_{i})=\mathbf{0}_{d_{2}\times 1}$. Hence,
$\boldsymbol{\gamma}_{i}^{0}=(\boldsymbol{\gamma}_{1,i}^{0},\boldsymbol{\gamma}_{2,i}^{0\prime})^{\prime}=(1,\boldsymbol{\gamma}_{2}^{0\prime}+\boldsymbol{\eta}_{i}^{\prime})^{\prime}$
and
$\mathbf{f}_{t}^{0}=(\mathbf{f}_{1,t}^{0},\mathbf{f}_{2,t}^{0\prime})^{\prime}$,
which in turn implies that

$\boldsymbol{\gamma}_{i}^{0\prime}\mathbf{f}_{t}^{0}=\mathbf{f}_{1,t}^{0}+(\boldsymbol{\gamma}_{2}^{0\prime}+\boldsymbol{\eta}_{i}^{\prime})\mathbf{f}_{2,t}^{0}=(\mathbf{f}_{1,t}^{0}+\boldsymbol{\gamma}_{2}^{0\prime}\mathbf{f}_{2,t}^{0})+\boldsymbol{\eta}_{i}^{\prime}\mathbf{f}_{2,t}^{0}.$
(3.5)

The common component with $\boldsymbol{\gamma}_{i}^{0}$ as loading and
$\mathbf{f}_{t}^{0}$ as factor can therefore be written equivalently as one in
which
$(\mathbf{f}_{1,t}^{0}+\boldsymbol{\gamma}_{2}^{0\prime}\mathbf{f}_{2,t}^{0},\mathbf{f}_{2,t}^{0\prime})^{\prime}$
is the factor and $(1,\boldsymbol{\eta}_{i}^{\prime})^{\prime}$ is the
loading.

######  Lemma 3.3 (Consistency of $\widehat{d}_{2},\ldots,\widehat{d}_{G+1}$
and $\widehat{\mathbf{F}}$).

Suppose that Assumptions 1–5 are satisfied. Then, the following results hold
as $N,\,T\to\infty$:

  * $\mathrm{(a)}$

$P(\widehat{d}_{g}=d_{g})\to 1$ for $g=2,\ldots,G+1$, where $G>1$ and
$d_{G+1}=0$;

  * $\mathrm{(b)}$

$\|\mathbf{P}_{\widehat{F}}-\mathbf{P}_{F^{0}}\|=o_{P}(1)$.

Together with Lemma 3.2, Lemma 3.3 (a) implies that
$\widehat{d}_{1},\ldots,\widehat{d}_{G+1}$ are all consistent. The consistency
of $\widehat{d}_{1},\ldots,\widehat{d}_{G}$ is important for obvious reasons.
The consistency of $\widehat{d}_{G+1}$ ensures that the stopping rule in Step
2 of the estimation procedure is asymptotically valid, which in turn implies
that $\widehat{G}$ is consistent. Hence, under the conditions of Lemma 3.3, we
have

$\displaystyle P(\widehat{G}=G)\to 1$ (3.6)

as $N,\,T\to\infty$. The consistency of $\widehat{G}$ further implies that
$d_{f}$ can be consistently estimated using
$\widehat{d}_{f}=\widehat{d}_{1}+\cdots+\widehat{d}_{\widehat{G}}$.

As we alluded to in the discussion of Assumption 5 above, $\mathbf{F}^{0}$ and
$\boldsymbol{\gamma}_{i}^{0}$ are only identified up to a rotation matrix
(see, for example, Bai and Ng, 2002). However, we cannot claim that
$\widehat{\mathbf{F}}$ is rotationally consistent for $\mathbf{F}^{0}$, as the
number of rows of both objects is growing with $T$. We therefore have to
resort to alternative consistency concepts. This is where Lemma 3.3 (b) comes
in. It shows that the spaces spanned by $\widehat{\mathbf{F}}$ and
$\mathbf{F}^{0}$ are asymptotically the same. This establishes that all the
Step 2 estimates are consistent. We therefore move on to investigate the
asymptotic properties of the IPC estimator $\widehat{\boldsymbol{\beta}}$ of
$\boldsymbol{\beta}^{0}$ obtained in Step 3.

######  Assumption 6 (Rates).

  * (a)

$N/T^{\nu_{G}}\to\rho_{1}\in[0,\infty)$;

  * (b)

$T^{2-\nu_{G}}/N\to\rho_{2}\in[0,\infty)$.

It is important to note that if all the factors in $\mathbf{f}_{t}^{0}$ are
stationary such that $G=1$ and $\nu_{G}=\nu_{1}=1$, then Assumption 6 (a) and
(b) require that $N/T\to\rho_{1}=1/\rho_{2}\in(0,\infty)$, which is the
condition considered in Bai (2009), and Moon and Weidner (2015). Assumption 6
(a) and (b) are therefore more general than the condition considered in these
other papers.

######  Assumption 7 (Asymptotic normality).

$\displaystyle\frac{1}{\sqrt{NT}}\sum_{i=1}^{N}\mathbf{D}_{T}\mathbf{Z}_{i}(\mathbf{F}^{0})^{\prime}\boldsymbol{\varepsilon}_{i}\to_{D}MN(\mathbf{0}_{d_{x}\times
1},\boldsymbol{\Omega})$ (3.7)

as $N,\,T\to\infty$, where

$\displaystyle\boldsymbol{\Omega}=\operatorname*{p\\!\lim}_{N,T\to\infty}\frac{1}{NT}\sum_{i=1}^{N}\sum_{j=1}^{N}\mathbf{D}_{T}E[\mathbf{Z}_{i}(\mathbf{F}^{0})^{\prime}\boldsymbol{\varepsilon}_{i}\boldsymbol{\varepsilon}_{j}\mathbf{Z}_{j}(\mathbf{F}^{0})|\mathcal{C}]\mathbf{D}_{T},$
(3.8)

with $\mathcal{C}$ being the sigma-algebra generated by $\mathbf{F}^{0}$.

Assumption 7 is a central limit theorem that is analogous to Assumption E of
Bai (2009). The reason for requiring that the asymptotic distribution is mixed
normal as opposed to normal is that by doing so we can accommodate
stochastically integrated factors (see, for example, Bai et al., 2009, and
Huang et al., 2021). In the absence of such integrated factors, the mixed
normal becomes normal. Either way, Assumption 7 ensures that standard normal
and chi-squared inference based on $\widehat{\boldsymbol{\beta}}$ is possible.

######  Theorem 3.1 (Asymptotic distribution of
$\widehat{\boldsymbol{\beta}}$).

Under Assumptions 1–7, as $N,\,T\to\infty$,

$\displaystyle\sqrt{NT}\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})\to_{D}MN(\mathbf{B}_{0}^{-1}(\sqrt{\rho_{1}}\mathbf{A}_{1}+\sqrt{\rho_{2}}\mathbf{A}_{2}),\mathbf{B}_{0}^{-1}\boldsymbol{\Omega}\mathbf{B}_{0}^{-1})$
(3.9)

where

$\displaystyle\mathbf{B}_{0}$
$\displaystyle=\operatorname*{p\\!\lim}_{N,T\to\infty}E[\mathbf{B}(\mathbf{F}^{0})|\mathcal{C}],$
(3.10) $\displaystyle\mathbf{A}_{1}$
$\displaystyle=-\operatorname*{p\\!\lim}_{N,T\to\infty}\frac{1}{T^{(1-\nu_{G})/2}}\sum_{i=1}^{N}\mathbf{D}_{T}E[\mathbf{X}_{i}^{\prime}\mathbf{M}_{F^{0}}\boldsymbol{\Sigma}_{\varepsilon}\mathbf{F}_{G}^{0}(\mathbf{F}_{G}^{0\prime}\mathbf{F}_{G}^{0})^{-1}(\boldsymbol{\Gamma}_{G}^{0\prime}\boldsymbol{\Gamma}_{G}^{0})^{-1}\boldsymbol{\gamma}_{G,i}^{0}|\mathcal{C}],$
(3.11) $\displaystyle\mathbf{A}_{2}$
$\displaystyle=-\operatorname*{p\\!\lim}_{N,T\to\infty}\frac{1}{T^{(3-\nu_{G})/2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\mathbf{D}_{T}E[\mathbf{Z}_{i}(0)^{\prime}\mathbf{F}_{G}^{0}(\mathbf{F}_{G}^{0\prime}\mathbf{F}_{G}^{0})^{-1}(\boldsymbol{\Gamma}_{G}^{0\prime}\boldsymbol{\Gamma}_{G}^{0})^{-1}\boldsymbol{\gamma}_{G,j}^{0}\boldsymbol{\varepsilon}_{j}^{\prime}\boldsymbol{\varepsilon}_{i}|\mathcal{C}],$
(3.12)

with
$\boldsymbol{\Sigma}_{\varepsilon}=N^{-1}\sum_{i=1}^{N}\boldsymbol{\Sigma}_{\varepsilon,i}$
and $\mathbf{Z}_{i}(0)=\mathbf{X}_{i}-\sum_{j=1}^{N}\mathbf{X}_{j}a_{ij}$.

According to Theorem 3.1, the asymptotic bias is driven by the factors and
loadings of group $G$, which is intuitive as the factors of this group are
smallest in order of magnitude. They therefore dominate the asymptotic bias.
By bounding $\nu_{G}$ from below Assumption 6 ensures that
$\mathbf{B}_{0}^{-1}(\sqrt{\rho_{1}}\mathbf{A}_{1}+\sqrt{\rho_{2}}\mathbf{A}_{2})$
is not diverging. It is important to note that $\rho_{1}$ and $\rho_{2}$ may
both be zero, which will be the case if $\nu_{G}>1$ and $T/N$ converges to a
constant. In general, the larger $\nu_{G}$ is, the less restrictive the
condition on $T/N$ has to be for $\rho_{1}$ and $\rho_{2}$ to be zero. For
example, if $\nu_{G}=2$, then we only require $N/T^{2}\to 0$. In this sense,
trending factors are a blessing. In order to put this into perspective,
suppose again that all the factors are stationary, such that $G=1$,
$\nu_{G}=\nu_{1}=1$ and $\rho_{1}=1/\rho_{2}\in(0,\infty)$, and that the
regressors are stationary, too, such that $\mathbf{D}_{T}=\mathbf{I}_{d_{x}}$.
In this case, the bias in Theorem 3.1 reduces to
$\mathbf{B}_{0}^{-1}(\sqrt{\rho_{1}}\mathbf{A}_{1}+\rho_{1}^{-1/2}\mathbf{A}_{2})$,
which is identical to the bias reported in Theorem 3 of Bai (2009). In fact,
it is not difficult to show that under these conditions,
$\sqrt{NT}(\widehat{\boldsymbol{\beta}}-\widehat{\boldsymbol{\beta}}_{0})=o_{P}(1)$,
so that the IPC estimator is asymptotically equivalent to the first step PC
estimator. The point about the bias is that while
$\sqrt{\rho_{1}}\mathbf{A}_{1}$ can be made arbitrarily small (large) by just
taking $\rho_{1}$ to zero (infinity), this will make
$\rho_{1}^{-1/2}\mathbf{A}_{2}$ divergent (negligible). It follows that unless
$\rho_{1}\in(0,\infty)$ the bias will diverge, which in turn means that there
is no way to make the bias disappear by just manipulating $\rho_{1}$, which in
practical terms means restricting $T/N$. By allowing $\nu_{G}>1$, we break the
inverse relationship between $\rho_{1}$ and $\rho_{2}$, which means that one
can be zero without for that matter forcing the other to infinity.

######  Remark 3.1.

Bai (2009, Corollary 1) shows that PC can be unbiased even if
$\mathbf{f}_{t}^{0}$ and $\mathbf{x}_{i,t}$ are stationary. However, this
requires either that $\varepsilon_{i,t}$ is independent and identically
distributed across both $i$ and $t$, or that it is uncorrelated and
homoskedastic in $t$ ($i$) with $T/N\to 0$ ($N/T\to 0$).

######  Remark 3.2.

Note that while biased, $\widehat{\boldsymbol{\beta}}$ is still consistent at
the best achievable rate. Again, if $\mathbf{x}_{i,t}$ is stationary,
$\mathbf{D}_{T}=\mathbf{I}_{d_{x}}$ and the rate of convergence is given by
$1/\sqrt{NT}$, which is the same as in Bai (2009). This is in contrast to
Lemma 3.1 and the relatively slow rate of convergence reported there. The
reason for this difference is that unlike Theorem 3.1, which requires that
Assumptions 1–7 all hold, Lemma 3.1 only requires Assumptions 1 and 2, and
under these very relaxed conditions the Theorem 3.1 rate is not attainable.
If, however, the conditions of Theorem 3.1 are met, then we can show that both
$\widehat{\boldsymbol{\beta}}_{0}$ and $\widehat{\boldsymbol{\beta}}_{1}$ are
consistent at the same rate as $\widehat{\boldsymbol{\beta}}$ (see the
appendix for a formal proof). However, unlike $\widehat{\boldsymbol{\beta}}$,
they are generally not mixed normal, which we verify using Monte Carlo
simulations in Section 4.

######  Corollary 3.1 (Unbiased asymptotic distribution).

Suppose that the conditions of Theorem 3.1 are met and that
$\rho_{1}=\rho_{2}=0$. Then, as $N,\,T\to\infty$,

$\displaystyle\sqrt{NT}\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})\to_{D}MN(\mathbf{0}_{d_{x}\times
1},\mathbf{B}_{0}^{-1}\boldsymbol{\Omega}\mathbf{B}_{0}^{-1}).$ (3.13)

In the appendix, we provide some alternative conditions that ensure
$\mathbf{A}_{1}=\mathbf{A}_{2}=\mathbf{0}_{d_{x}\times 1}$. If $\rho_{1}$,
$\rho_{2}$, $\mathbf{A}_{1}$ and $\mathbf{A}_{2}$ are all different from zero,
one possibility is to use bias correction. In the appendix, we explain how the
Jackknife approach can be used for this purpose.

Theorem 3.1 imposes only minimal conditions on the correlation and
heteroskedasticity of $\varepsilon_{i,t}$, and is in this sense very general.
Such generality is, however, not possible if we also want to ensure consistent
estimation of $\boldsymbol{\Omega}$. Let us therefore assume for a moment that
$E(\varepsilon_{i,t}\varepsilon_{j,s})=0$ for all $(i,t)\neq(j,s)$, so that
$\varepsilon_{i,t}$ is serially and cross-sectionally uncorrelated. In this
case,

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}E[\mathbf{Z}_{i}(\mathbf{F}^{0})^{\prime}\boldsymbol{\varepsilon}_{i}\boldsymbol{\varepsilon}_{j}\mathbf{Z}_{j}(\mathbf{F}^{0})|\mathcal{C}]=\sum_{i=1}^{N}\sigma_{\varepsilon,i}^{2}E[\mathbf{Z}_{i}(\mathbf{F}^{0})^{\prime}\mathbf{Z}_{i}(\mathbf{F}^{0})|\mathcal{C}].$
(3.14)

A natural estimator of this matrix is given by

$\displaystyle\sum_{i=1}^{N}\widehat{\sigma}_{\varepsilon,i}^{2}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i},$
(3.15)

where
$\widehat{\sigma}_{\varepsilon,i}^{2}=T^{-1}\sum_{t=1}^{T}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}})^{\prime}\mathbf{M}_{\widehat{F}}(\mathbf{y}_{i}-\mathbf{X}_{i}\widehat{\boldsymbol{\beta}})$
and $\widehat{\mathbf{Z}}_{i}$ is as in the definition of
$\widehat{\boldsymbol{\beta}}$. It is not difficult to show that under the
conditions of Theorem 3.1,

$\displaystyle\left\|\frac{1}{NT}\sum_{i=1}^{N}\widehat{\sigma}_{\varepsilon,i}^{2}\mathbf{D}_{T}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i}\mathbf{D}_{T}-\boldsymbol{\Omega}\right\|=o_{P}(1).$
(3.16)

Of course, in this paper we do not assume knowledge of the order of the
regressors, which in practice means that the appropriate normalization matrix
$\mathbf{D}_{T}$ to use is unknown. This is not a problem, however, as the
usual Wald and $t$-test statistics are self-normalizing. As an illustration,
consider testing the null hypothesis of
$H_{0}:\mathbf{R}\boldsymbol{\beta}^{0}=\mathbf{r}$, where $\mathbf{R}$ is a
$r_{0}\times d_{x}$ matrix of rank $r_{0}\leq d_{x}$ and $\mathbf{r}$ is a
$r_{0}\times 1$ vector. The Wald test statistic for testing this hypothesis is
given by

$\displaystyle
W_{\widehat{\boldsymbol{\beta}}}=(\mathbf{R}\widehat{\boldsymbol{\beta}}-\mathbf{r})^{\prime}\left[\mathbf{R}\left(\sum_{i=1}^{N}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i}\right)^{-1}\sum_{i=1}^{N}\widehat{\sigma}_{\varepsilon,i}^{2}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i}\left(\sum_{i=1}^{N}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i}\right)^{-1}\mathbf{R}\right]^{-1}(\mathbf{R}\widehat{\boldsymbol{\beta}}-\mathbf{r}),$
(3.17)

which has a limiting chi-squared distribution with $r_{0}$ degrees of freedom
under $H_{0}$, as is clear from

$\displaystyle W_{\widehat{\boldsymbol{\beta}}}$
$\displaystyle=[\mathbf{R}\mathbf{D}_{T}\sqrt{NT}\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})]^{\prime}\Bigg{[}\mathbf{R}\mathbf{D}_{T}\left(\frac{1}{NT}\sum_{i=1}^{N}\mathbf{D}_{T}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i}\mathbf{D}_{T}\right)^{-1}\frac{1}{NT}\sum_{i=1}^{N}\widehat{\sigma}_{\varepsilon,i}^{2}\mathbf{D}_{T}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i}\mathbf{D}_{T}$
$\displaystyle\times\left(\frac{1}{NT}\sum_{i=1}^{N}\mathbf{D}_{T}\widehat{\mathbf{Z}}_{i}^{\prime}\widehat{\mathbf{Z}}_{i}\mathbf{D}_{T}\right)^{-1}\mathbf{D}_{T}\mathbf{R}\Bigg{]}^{-1}\mathbf{R}\mathbf{D}_{T}\sqrt{NT}\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})\to_{D}\chi^{2}(r_{0}).$
(3.18)

In the next section, we use Monte Carlo simulations as a means to evaluate the
accuracy of this last results in small samples.

The above results are for the case when $\varepsilon_{i,t}$ is serially and
cross-sectionally uncorrelated. If $\varepsilon_{i,t}$ is serially and/or
cross-sectionally correlated, we recommend following Bai (2009), who discusses
the issue of consistent covariance matrix estimation at length. The same
arguments can be applied without change in current context.

##  4 Monte Carlo results

This section reports the results obtained from a small-scale Monte Carlo
simulation exercise. One important goal of the simulations is to examine the
finite sample performance of $\widehat{\boldsymbol{\beta}}_{0}$,
$\widehat{\boldsymbol{\beta}}_{1}$, and $\widehat{\boldsymbol{\beta}}$, and
justify the fact that $\widehat{\boldsymbol{\beta}}$ is superior to
$\widehat{\boldsymbol{\beta}}_{0}$ and $\widehat{\boldsymbol{\beta}}_{1}$.

We consider two DGPs, henceforth denoted “DGP 1” and “DGP 2”, which we now
describe.

######  DGP 1 (Artificial data).

In this DGP, $y_{i,t}$ is generated according to a restricted version of (2.1)
that sets $d_{x}=2$, $\boldsymbol{\beta}^{0}=\mathbf{1}_{2\times 1}$,
$\varepsilon_{i,t}\sim N(0,1)$ and $N,\,T\in\\{40,80,160,320\\}$. We further
set $d_{f}=3$ and generate the factor loadings as $\gamma_{1,i}^{0}\sim
N(1,1)$, $\gamma_{2,i}^{0}\sim N(0,1)$ and $\gamma_{3,i}^{0}\sim N(0,1)$. The
factors themselves are generated as $f_{1,t}^{0}=t$, $f_{2,t}^{0}=\mu_{t}$ and
$f_{3,t}^{0}=c_{t}$, where $\mu_{t}=\mu_{t-1}+\xi_{t}$, $\mu_{0}=0$,
$\xi_{t}\sim N(0,1/4)$ and $c_{t}=\sin(8\pi t/T)$. Hence, in this DGP, the
common component is a random walk with drift and cycle. Also,
$d_{1}=d_{2}=d_{3}=1$ and $(\nu_{1},\nu_{2},\nu_{3})=(3,2,1)$. Let us denote
by $x_{j,i,t}$ the $j$-th element of $\mathbf{x}_{i,t}$. The following
specification makes $x_{j,i,t}$ correlated with the common component of
$y_{i,t}$:

$\displaystyle
x_{j,i,t}=\frac{1}{d_{x}}\left(\sum_{j=1}^{d_{f}}|\gamma_{j,i}^{0}|+|\xi_{t}|+|c_{t}|\right)+\left(\frac{t}{4}\right)^{(j-1)/4}+v_{j,i,t},$
(4.1)

where $v_{j,i,t}$ is the $i$-th element of the $N\times 1$ vector
$\mathbf{v}_{j,t}=(v_{j,1,t},\ldots,v_{j,N,t})^{\prime}$, which we generate as

$\displaystyle\mathbf{v}_{j,t}=0.5\mathbf{v}_{j,t-1}+\boldsymbol{\omega}_{j,t},$
(4.2)

where $\boldsymbol{\omega}_{j,t}\sim
N(\mathbf{0},\boldsymbol{\Sigma}_{\omega})$ and $\boldsymbol{\Sigma}_{\omega}$
has $0.5^{|m-n|}$ in row $m$ and column $n$. Thus, $v_{j,i,t}$ is weakly
correlated across both $i$ and $t$.

######  DGP 2 (Calibrated data).

This DGP is calibrated by using the US commercial banks data set studied in
Section A.3 in which $N=466$ and $T=80$. The $d_{x}=5$ regressors are the same
as those included in Section 5 and
$\boldsymbol{\beta}^{0}=(0.4063,0.0199,0.3805,0.0912,0.04726)^{\prime}$ is
calibrated using the empirical estimates taken from Table A.1. The factors are
the estimated ones reported in Figure A.2. Hence, in this DGP, $d_{f}=2$. The
loadings are generated as in DGP 1 with $\gamma_{1,i}^{0}\sim N(1,1)$ and
$\gamma_{2,i}^{0}\sim N(0,1)$.

For each combination of $N$ and $T$, we report the correct selection frequency
for $(\widehat{d}_{1},\ldots,\widehat{d}_{\widehat{G}})$ when seen as an
estimator of $(d_{1},\ldots,d_{G})$ and for $\widehat{d}_{g}$ individually for
each group $g=1,\ldots,G$. Hence, while the former frequency captures the
accuracy of the estimation of both $(d_{1},\ldots,d_{G})$ and $G$, the latter
frequency only captures the accuracy of the estimation of $d_{g}$ for each
$g=1,\ldots,G$. We also report the root mean squared error (RMSE) of
$\widehat{\boldsymbol{\beta}}$ and $\mathbf{P}_{\widehat{F}}$, as measured by
the square root of the average of
$\|\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0}\|^{2}$ and
$\|\mathbf{P}_{\widehat{F}}-\mathbf{P}_{F^{0}}\|^{2}$, respectively, over the
replications. In interest of comparison, RMSEs of
$\widehat{\boldsymbol{\beta}}$ is compared to that of
$\widehat{\boldsymbol{\beta}}_{0}$, $\widehat{\boldsymbol{\beta}}_{1}$ and the
infeasible OLS estimator of $\boldsymbol{\beta}^{0}$ based on taking
$\mathbf{F}^{0}$ as known, $\widehat{\boldsymbol{\beta}}(\mathbf{F}^{0})$.
Some results on the 5% size of the Wald test based on
$\mathbf{R}=\mathbf{I}_{d_{x}}$ and $\mathbf{r}=\boldsymbol{\beta}^{0}$ are
also reported. In particular, we report the size of
$W_{\widehat{\boldsymbol{\beta}}}$, $W_{\widehat{\boldsymbol{\beta}}_{0}}$ and
$W_{\widehat{\boldsymbol{\beta}}_{1}}$, which are computed in an obvious
fashion by replacing $(\widehat{\boldsymbol{\beta}},\widehat{\mathbf{F}})$
with $(\widehat{\boldsymbol{\beta}}_{0},\widehat{\mathbf{F}}_{0})$ and
$(\widehat{\boldsymbol{\beta}}_{1},\widehat{\mathbf{F}})$, respectively, and
$W_{\widehat{\boldsymbol{\beta}}(\boldsymbol{F}^{0})}$, which is calculated in
the same way as $W_{\widehat{\boldsymbol{\beta}}}$ but with
$\mathbf{M}_{F}\mathbf{X}_{i}$ in place of $\widehat{\mathbf{Z}}_{i}$. The
critical values are taken from $\chi^{2}(d_{x})$. As pointed out in Section 2,
the IPC estimator is invariant with respect to $\delta$. The results reported
here are based on $\delta=1$. In time series analysis, rules such as (the
integer part of) $8(T/100)^{1/4}$ are sometimes used for the maximum model
dimension to be considered. In our theory, however, $d_{max}$ is fixed and
some of the arguments used in the proofs fail if $d_{max}$ is allowed to
increase with $N$ and $T$. In this section, we therefore follow the bulk of
the previous literature and set $d_{max}$ to a fixed large number (see, for
example, Ahn and Horenstein, 2013, Bai and Ng, 2002, and Moon and Weidner,
2015). We choose $d_{max}=10$, which led to the same results as some of the
other values we tried. The number of replications is set to 1,000.

We begin by considering the results for DGP 1, which are reported in Tables 1
and 2. In particular, while Table 1 contains the results for the estimated
common component, Table 2 contains the results for the estimated slopes. We
first look at the correct selection frequencies for all and each individual
group of factors reported in Table 1. The results for the individual groups
suggest that the accuracy is generally degreasing in $g$, which is partly
expected because the signal strength of the factors, as measured by $\nu_{g}$,
decreases with increasing values of $g$. The accuracy of
$(\widehat{d}_{1},\ldots,\widehat{d}_{\widehat{G}})$ is almost identical to
that of $\widehat{d}_{G}$, suggesting that the accuracy of $\widehat{G}$ is
driven by the accuracy of the group whose factors has the weakest signal.
Looking next at the RMSE results reported in Table 2 for estimating
$\boldsymbol{\beta}^{0}$, we see that there is a clear improvement as the
sample size increases. The best overall performance is generally obtained when
taking $\mathbf{F}^{0}$ as known, which is in accordance with our priori
expectations. However, the improvement is not very large and it decreases with
increases in $N$ and $T$. The reason for this is the accuracy of the estimated
factors, which according to the RMSE of $\mathbf{P}_{\widehat{F}}$ reported
Table 1 is high and increasing in $N$ and $T$. The second best performance is
obtained by using $\widehat{\boldsymbol{\beta}}$, followed by
$\widehat{\boldsymbol{\beta}}_{1}$, and then
$\widehat{\boldsymbol{\beta}}_{0}$. If our asymptotic theory is correct, while
the size of $W_{\widehat{\boldsymbol{\beta}}}$ and
$W_{\widehat{\boldsymbol{\beta}}(\boldsymbol{F}^{0})}$ should converge to 5%
as the sample size increases, that of $W_{\widehat{\boldsymbol{\beta}}_{0}}$
and $W_{\widehat{\boldsymbol{\beta}}_{1}}$ should not, and this is exactly
what we see in Table 2. Note in particular how the size of
$W_{\widehat{\boldsymbol{\beta}}_{0}}$ and
$W_{\widehat{\boldsymbol{\beta}}_{1}}$ is not only nonconvergent but that it
is in fact increasing in $N$ and $T$.

The RMSEs of $\widehat{\boldsymbol{\beta}}_{0}$,
$\widehat{\boldsymbol{\beta}}_{1}$, $\widehat{\boldsymbol{\beta}}$ and
$\widehat{\boldsymbol{\beta}}(\mathbf{F}^{0})$ in DGP 2 are given by 0.0444,
0.0405, 0.0405 and 0.0391, respectively, and the sizes of the Wald tests
associated with these estimators are given by 0.285, 0.103, 0.078 and 0.061,
respectively. Hence, just as in DGP 1, the best performance is obtained by
using $\widehat{\boldsymbol{\beta}}(\mathbf{F}^{0})$ with
$\widehat{\boldsymbol{\beta}}$, $\widehat{\boldsymbol{\beta}}_{1}$ and
$\widehat{\boldsymbol{\beta}}_{0}$ ending up in second, third and fourth
place, respectively.

##  5 Empirical illustration — house prices and income

Economists have become concerned that recently house prices have grown too
quickly, and that prices are now too high relative to per capita incomes. If
this is correct and there is any truth to the theory on the matter, prices
should stagnate or fall until they are better aligned with income, which in
statistical terms mean that house prices should be cointegrated with income.
The validity of this assumption has important implications for policy, because
a failure could be due to a housing bubble.

In this section, we revisit the real house price data set of Holly et al.
(2010), which comprises data on log real house prices ($p_{i,t}$) and log real
per capita income ($w_{i,t}$) for 49 US states across the 1975–2003 period.
According to theory, $p_{i,t}$ and $w_{i,t}$ should be cointegrated with
cointegrating vector $(1,\,-1)^{\prime}$. The previous empirical evidence of
this prediction has, however, been mixed and far from convincing (see, for
example, Gallin, 2006). Holly et al. (2010) argue that this lack of empirical
support can be attributed in part to a failure to account for cross-sectional
dependence, leading to deceptive conclusions. The authors therefore apply the
“CIPS” panel unit root test of Pesaran (2007), which allow for cross-section
dependence in the form of a common factor. The test is applied both to
$p_{i,t}$ and $w_{i,t}$ separately, and to $p_{i,t}-w_{i,t}$. According to the
results, while the variables are unit root non-stationary, their difference is
not. Holly et al. (2010) also report CCE results suggesting that the estimated
income elasticity is indeed close to one. They therefore conclude that
$p_{i,t}$ and $w_{i,t}$ are cointegrated with cointegrating vector
$(1,\,-1)^{\prime}$, just as predicted by theory.

Figure 1: The variables in the house prices and income illustration.

Our interest in the work of Holly et al. (2010) stems from their preference to
apply the CIPS test, which tests for a unit root in the defactored data. This
means that if $p_{i,t}$ and $w_{i,t}$ are not cointegrated by themselves, but
only when conditioning on unit root common factors, because of the way that
the data are defactored prior to the testing, the unit root null hypothesis is
likely to be rejected by the CIPS test. That is, the test is likely to lead to
the conclusion of cointegration when in fact there is none. In this section,
we use IPC as a means to investigate this possibility.

Figure 2: The estimated factors in the house prices and income illustration.

As in the RTS illustration, we begin by plotting the variables. This is done
in Figure 1. As expected, both variables are highly persistent and the ADF
test provides no evidence against the unit root null. This corroborates the
unit root test results reported by Holly et al. (2010). However, we also see
that the trending behaviours of $p_{i,t}$ and $w_{i,t}$ are very different,
suggesting that their stochastic trends are not the same, which they should be
under cointegration. We also see that the trending behaviour is very similar
across states, which is suggestive of non-stationary common factors. Of
course, the IPC procedure does not require cointegration and it does allow for
very general types of factors. We therefore proceed with the estimation of the
model. Hence, in this illustration, $y_{i,t}=p_{i,t}$ and
$\mathbf{x}_{i,t}=w_{i,t}$. For comparison purposes, the IPC results are
presented together with the results obtained by applying the PC estimator of
Bai (2009), as well as the usual OLS estimator with time and state fixed
effects. The results reported in Table 3. The first thing to note is that the
estimated slopes vary a lot depending on the estimator used. Interestingly,
the point estimates are increasing in the generality of the estimator with
fixed effects OLS (IPC) leading to the lowest (highest) estimate. We therefore
begin by considering the IPC results. The point estimate of 2.1024 is far from
the theoretically predicted value of one, which is also not included in the
reported 95% confidence interval. As for the factors, we estimate
$\widehat{d}_{1}=\widehat{d}_{2}=\widehat{d}_{3}=1$, implying that
$\widehat{d}_{f}=3$. The estimated factors, denoted $\widehat{f}_{1,t}$,
$\widehat{f}_{2,t}$ and $\widehat{f}_{3,t}$, are plotted in Figure 2. As
expected given Figure 1, all three factors are highly persistent, although
with clearly distinct trends. We take this as evidence against cointegration
between $p_{i,t}$ and $w_{i,t}$, since under cointegration the factors should
be stationary.

Because we estimate three time-varying factors, fixed effects OLS is invalid,
as it only allows for a common time effect. Moreover, since the factors come
from three distinct groups, and are not all stationary, PC is invalid, too.
This leaves us with the IPC estimator, which again provides strong evidence
against the theoretically predicted one-to-one cointegrated relationship
between $p_{i,t}$ and $w_{i,t}$. Holly et al. (2010, page 172) conclude that
“[o]ur results support the hypothesis that real house prices have been rising
in line with fundamentals (real incomes), and there seems little evidence of
house price bubbles at the national level.” The results reported here reveal a
completely different picture with housing prices being long run disconnected
with real income.

##  6 Conclusion

The PC approach of Bai (2009) has attracted considerable interest in recent
years, so much so that it has given rise to a separate PC literature. A key
assumption in this literature is that both the unknown factors and regressors
are stationary, which is rarely the case in practice. In the present paper, we
relax this assumption by considering a very general specification in which the
factors and regressors are essentially unrestricted. In spite of this
generality, the proposed IPC estimator can be applied without any input from
the practitioner, except for the maximum number of factors to be considered.
The fact that in IPC there is no need to distinguish between deterministic and
stochastic factors means that the usual problem in applied work of deciding on
which deterministic terms to include in the model does not arise, as these are
estimated along with the other factors of the model. There is also no need to
pre-test the regressors for unit roots, which is otherwise standard practice
when using procedures that do not require data to be stationary. In other
words, the proposed IPC is not only very general but also extremely user-
friendly. It should therefore be a valuable addition to the already existing
menu of techniques for panel regression models with interactive effects.

## References

  * (1)
  * Ahn and Horenstein (2013) Ahn, S. C. and Horenstein, A. R. (2013), ‘Eigenvalue ratio test for the number of factors’, Econometrica 81, 1203–1227.
  * Ahn et al. (2013) Ahn, S., Lee, Y. H. and Schmidt, P. (2013), ‘Panel data models with multiple time-varying individual effects’, Journal of Econometrics 174, 1–14.
  * Ando and Bai (2017) Ando, T. and Bai, J. (2017), ‘Clustering huge number of financial time series: A panel data approach with high-dimensional predictors and factor structures’, Journal of the American Statistical Association 112, 1182–1198.
  * Bai (2004) Bai, J. (2004), ‘Estimating cross-section common stochastic trends in nonstationary panel data’, Journal of Econometrics 122, 137–183.
  * Bai (2009) Bai, J. (2009), ‘Panel data models with interactive fixed effects’, Econometrica 77, 1229–1279.
  * Bai et al. (2009) Bai, J., Kao, C. and Ng, S. (2009), ‘Panel cointegration with global stochastic trends’, Journal of Econometrics 149, 82–99.
  * Bai and Li (2014) Bai, J. and Li, K. (2014), ‘Theory and methods of panel data models with interactive effects’, Annals of Statistics 42, 142–170.
  * Bai and Ng (2002) Bai, J. and Ng, S. (2002), ‘Determining the number of factors in approximate factor models’, Econometrica 70, 191–221.
  * Bai and Ng (2004) Bai, J. and Ng, S. (2004), ‘A Panic attack on unit roots and cointegration’, Econometrica 72, 1127–1177.
  * Barigozzi and Trapani (2021) Barigozzi, M. and Trapani, L. (2021), ‘Testing for common trends in non-stationary large datasets’, Journal of Business and Economic Statistics p. forthcoming.
  * Castagnetti et al. (2015) Castagnetti, C., Rossi, E. and Trapani, L. (2015), ‘Inference on factor structures in heterogeneous panels’, Journal of Econometrics 184, 145 – 157.
  * Castagnetti et al. (2019) Castagnetti, C., Rossi, E. and Trapani, L. (2019), ‘A two-stage estimator for heterogeneous panel models with common factors’, Econometrics and Statistics 11, 63 – 82.
  * Chen et al. (2021) Chen, M., Fernández-Val, I. and Weidner, M. (2021), ‘Nonlinear factor models for network and panel data’, Journal of Econometrics 220, 296–324.
  * Choi (2017) Choi, I. (2017), ‘Efficient estimation of nonstationary factor models’, Journal of Statistical Planning and Inference 183, 18–43.
  * Chudik et al. (2011) Chudik, A., Pesaran, M. H. and Tosetti, E. (2011), ‘Weak and strong cross section dependence and estimation of large panels’, Econometrics Journal 14, 45–90.
  * DeAngelo and Roll (2015) DeAngelo, H. and Roll, R. (2015), ‘How stable are corporate capital structures?’, Journal of Finance 70, 373–418.
  * Dhaene and Jochmans (2015) Dhaene, G. and Jochmans, K. (2015), ‘Split-panel Jackknife Estimation of Fixed-effect Models’, The Review of Economic Studies 82(3), 991–1030.
  * Dong et al. (2021) Dong, C., Gao, J. and Peng, B. (2021), ‘Varying-coefficient panel data models with nonstationarity and partially observed factor structure’, Journal of Business and Economic Statistics 39, 700–711.
  * Dong and Linton (2018) Dong, C. and Linton, O. (2018), ‘Additive nonparametric models with time variable and both stationary and nonstationary regressors’, Journal of Econometrics 207, 212–236.
  * Feng et al. (2017) Feng, G., Gao, J., Peng, B. and Zhang, X. (2017), ‘A varying-coefficient panel data model with fixed effects: Theory and an application to the US commercial banks’, Journal of Econometrics 196, 68–82.
  * Feng and Serletis (2008) Feng, G. and Serletis, A. (2008), ‘Productivity trends in U.S. manufacturing: Evidence from the NQ and AIM cost functions’, Journal of Econometrics 142, 281–311.
  * Feng and Zhang (2012) Feng, G. and Zhang, X. (2012), ‘Productivity and efficiency at large and community banks in the U.S.: A bayesian true random effects stochastic distance frontier analysis’, Journal of Banking and Finance 36, 1883–1895.
  * Fernández-Val and Weidner (2018) Fernández-Val, I. and Weidner, M. (2018), ‘Fixed effects estimation of large-$t$ panel data models’, Annual Review of Economics 10, 109–138.
  * Gallin (2006) Gallin, J. (2006), ‘The long run relationship between house prices and income: Evidence from local housing markets’, Real Estate Economics 34, 417–438.
  * Holly et al. (2010) Holly, S., Pesaran, M. H. and Yamagata, T. (2010), ‘A spatio-temporal model of house prices in the USA’, Journal of Econometrics 158, 160–173.
  * Huang et al. (2021) Huang, W., Jin, S., Phillips, P. C. B. and Su, L. (2021), ‘Nonstationary panel models with latent group structures and cross-section dependence’, Journal of Econometrics 221, 198–222.
  * Jiang et al. (2020) Jiang, B., Yang, Y., Gao, J. and Hsiao, C. (2020), ‘Recursive estimation in large panel data models: Theory and practice’, Journal of Econometrics p. forthcoming.
  * Lam and Yao (2012) Lam, C. and Yao, Q. (2012), ‘Factor modeling for high-dimensional time series: Inference for the number of factors’, Annals of Statistics 40, 694–726.
  * Lam et al. (2011) Lam, C., Yao, Q. and Bathia, N. (2011), ‘Estimation of latent factors for high-dimensional time series’, Biometrika 98, 901–918.
  * Lemmon et al. (2008) Lemmon, M., Roberts, M. and Zender, J. (2008), ‘Back to the beginning: Persistence and the cross-section of corporate capital structure’, Journal of Finance 60, 2575–2619.
  * Li et al. (2016) Li, D., Qian, J. and Su, L. (2016), ‘Panel data models with interactive fixed effects and multiple structural breaks’, Journal of the American Statistical Association 111, 1804–1819.
  * Li et al. (2020) Li, K., Cui, G. and Lu, L. (2020), ‘Efficient estimation of heterogeneous coefficients in panel data models with common shocks’, Journal of Econometrics 216, 327 – 353.
  * Magnus and Neudecker (2007) Magnus, J. R. and Neudecker, H. (2007), Matrix Differential Calculus with Applications in Statistics and Econometrics, third edn, John Wiley & Sons Ltd.
  * Moon and Weidner (2015) Moon, H. R. and Weidner, M. (2015), ‘Linear regression for panel with unknown number of factors as interactive fixed effects’, Econometrica 83, 1543–1579.
  * Park and Phillips (2000) Park, J. Y. and Phillips, P. C. B. (2000), ‘Nonstationary binary choice’, Econometrica 68, 1249–1280.
  * Pesaran (2006) Pesaran, M. H. (2006), ‘Estimation and inference in large heterogeneous panels with a multifactor error structure’, Econometrica 74, 967–1012.
  * Pesaran (2007) Pesaran, M. H. (2007), ‘A simple panel unit root test in the presence of cross section dependence’, Journal of Applied Econometrics 22, 265–312.
  * Stiroh (2000) Stiroh, K. (2000), ‘How did bank holding companies prosper in the 1990s?’, Journal of Banking and Finance 24, 1703–1745.
  * Uematsu and Yamagata (2020) Uematsu, Y. and Yamagata, T. (2020), Estimation of weak factor models.  Available online at: https://ssrn.com/abstract=3374750.
  * Westerlund (2018) Westerlund, J. (2018), ‘CCE in panels with general unknown factors’, Econometrics Journal 21, 264–276.
  * Wheelock and Wilson (2012) Wheelock, D. C. and Wilson, P. W. (2012), ‘Do large banks have lower costs? New estimates of returns to scale for US banks’, Journal of Money, Credit and Banking 44, 171–199.

Table 1: Monte Carlo results for the estimated common component in DGP 1.

|  | Correct selection frequency |  | RMSE
---|---|---|---|---
$N$ | $T$ | $\widehat{d}_{1},...,\widehat{d}_{\widehat{G}}$ | $\widehat{d}_{1}$ | $\widehat{d}_{2}$ | $\widehat{d}_{3}$ |  | $\mathbf{P}_{\widehat{F}}$
40 | 40 | 0.341 | 1.000 | 0.406 | 0.341 |  | 0.9453
| 80 | 0.628 | 1.000 | 0.646 | 0.628 |  | 0.5599
| 160 | 0.835 | 1.000 | 0.837 | 0.835 |  | 0.3715
| 320 | 0.954 | 1.000 | 0.954 | 0.954 |  | 0.3520
80 | 40 | 0.348 | 1.000 | 0.378 | 0.348 |  | 0.9452
| 80 | 0.661 | 1.000 | 0.668 | 0.661 |  | 0.4523
| 160 | 0.867 | 1.000 | 0.867 | 0.867 |  | 0.2634
| 320 | 0.958 | 1.000 | 0.958 | 0.958 |  | 0.2401
160 | 40 | 0.329 | 1.000 | 0.337 | 0.329 |  | 0.6401
| 80 | 0.684 | 1.000 | 0.686 | 0.684 |  | 0.4359
| 160 | 0.864 | 1.000 | 0.864 | 0.864 |  | 0.1761
| 320 | 0.960 | 1.000 | 0.960 | 0.960 |  | 0.1696
320 | 40 | 0.220 | 1.000 | 0.225 | 0.220 |  | 0.4043
| 80 | 0.603 | 1.000 | 0.603 | 0.603 |  | 0.2083
| 160 | 0.882 | 1.000 | 0.882 | 0.882 |  | 0.1326
| 320 | 0.988 | 1.000 | 0.988 | 0.988 |  | 0.1195

  * •

_Notes_ : The correct selection frequencies are for the estimated factor
groups, which is Step 2 of the IPC procedure. The results for
$\widehat{d}_{1},\ldots,\widehat{d}_{\widehat{G}}$ treat both the groups and
their number, $G$, as unknown, while the results for $\widehat{d}_{1}$,
$\widehat{d}_{2}$ and $\widehat{d}_{3}$ take $G$ as given. The reported RMSE
results for $\mathbf{P}_{\widehat{F}}$ refer to the square root of the average
of $\|\mathbf{P}_{\widehat{F}}-\mathbf{P}_{F^{0}}\|^{2}$ across the Monte
Carlo replications.

Table 2: Monte Carlo results for the estimated slopes in DGP 1.

|  | RMSE |  | 5% size
---|---|---|---|---
Estimator |  $N$\$T$ | 40 | 80 | 160 | 320 |  | 40 | 80 | 160 | 320
$\widehat{\boldsymbol{\beta}}_{0}$ | 40 | 0.0573 | 0.0423 | 0.0380 | 0.0410 |  | 0.6800 | 0.7280 | 0.8410 | 0.9160
| 80 | 0.0335 | 0.0234 | 0.0212 | 0.0220 |  | 0.5990 | 0.6550 | 0.7940 | 0.9050
| 160 | 0.0219 | 0.0169 | 0.0158 | 0.0169 |  | 0.5820 | 0.3500 | 0.8680 | 0.9560
| 320 | 0.0148 | 0.0109 | 0.0105 | 0.0142 |  | 0.5680 | 0.6230 | 0.7750 | 0.9040
$\widehat{\boldsymbol{\beta}}_{1}$ | 40 | 0.0435 | 0.0296 | 0.0252 | 0.0262 |  | 0.3330 | 0.4320 | 0.6080 | 0.7520
| 80 | 0.0286 | 0.0177 | 0.0144 | 0.0138 |  | 0.2520 | 0.3130 | 0.4720 | 0.6830
| 160 | 0.0172 | 0.0127 | 0.0104 | 0.0102 |  | 0.2300 | 0.7260 | 0.5610 | 0.7560
| 320 | 0.0113 | 0.0079 | 0.0069 | 0.0084 |  | 0.2330 | 0.2770 | 0.4580 | 0.7060
$\widehat{\boldsymbol{\beta}}$ | 40 | 0.0383 | 0.0212 | 0.0129 | 0.0091 |  | 0.1320 | 0.0950 | 0.0580 | 0.0680
| 80 | 0.0280 | 0.0146 | 0.0092 | 0.0062 |  | 0.1500 | 0.0670 | 0.0650 | 0.0550
| 160 | 0.0164 | 0.0105 | 0.0064 | 0.0043 |  | 0.0960 | 0.0750 | 0.0690 | 0.0430
| 320 | 0.0099 | 0.0064 | 0.0044 | 0.0032 |  | 0.0720 | 0.0480 | 0.0420 | 0.0660
$\widehat{\boldsymbol{\beta}}(\mathbf{F}^{0})$ | 40 | 0.0254 | 0.0171 | 0.0114 | 0.0081 |  | 0.0550 | 0.0610 | 0.0440 | 0.0410
| 80 | 0.0186 | 0.0117 | 0.0084 | 0.0057 |  | 0.0750 | 0.0470 | 0.0560 | 0.0410
| 160 | 0.0127 | 0.0086 | 0.0059 | 0.0041 |  | 0.0560 | 0.0600 | 0.0530 | 0.0490
| 320 | 0.0087 | 0.0059 | 0.0041 | 0.0030 |  | 0.0640 | 0.0420 | 0.0520 | 0.0660

  * •

_Notes_ : The RMSE of $\widehat{\boldsymbol{\beta}}$ refers to the square root
of the average of
$\|\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0}\|^{2}$ across the Monte
Carlo replications. The RMSEs of $\widehat{\boldsymbol{\beta}}_{0}$,
$\widehat{\boldsymbol{\beta}}_{1}$ and
$\widehat{\boldsymbol{\beta}}(\mathbf{F}^{0})$ are constructed in an analogous
fashion, where $\widehat{\boldsymbol{\beta}}_{0}$ is the initial Step 1 IPC
estimator, $\widehat{\boldsymbol{\beta}}_{1}$ is the regular PC estimator
based on the IPC estimator of the factors, $\widehat{\mathbf{F}}$, and
$\widehat{\boldsymbol{\beta}}(\mathbf{F}^{0})$ is the infeasible OLS estimator
based on the true factors. The 5% size results are for the Wald test
associated with each estimator.

Table 3: Empirical results for the house price and income illustration.

Estimator | Point estimate | Wald $p$-value
---|---|---
FE | 0.3453 | 0.0000
PC | 1.1602 | 0.0000
IPC | 2.1024 | 0.0000

  * •

_Notes_ : “FE”, “PC” and “IPC” refer to the two-way fixed effects OLS
estimator, the PC estimator of Bai (2009) and the proposed IPC estimator,
respectively, in a regression of $p_{i,t}$ onto $w_{i,t}$. The reported Wald
$p$-values test the null hypothesis that the relevant coefficient is zero.

Online Appendix to: “Interactive Effects Panel Data Models with General
Factors and Regressors”

Bin Peng⋆, Liangju Su‡, Joakim Westerlund♯ and Yanrong Yang†

⋆Monash University, ‡Tsinghua University, ♯Lund University, †Australian
National University

We begin this appendix by laying out the notation that will be used
throughout. This is done in Section A.1, which in Section A.2 is followed by a
discussion of some conditions that ensure that the asymptotic distribution of
$\widehat{\boldsymbol{\beta}}$ is correctly centered at zero. Section A.3
provides an extra empirical study using the US bank data. Section A.4
describes the outline of the proofs. Sections A.5 and A.6 provide some
auxiliary lemmas and their proofs, respectively. Proofs of the main results
are provided in Section A.7.

##  A.1 Notation

The matrices $\boldsymbol{\Sigma}_{F^{0}}$ and
$\boldsymbol{\Sigma}_{\Gamma^{0}}$ have been defined in Assumption 1. In this
appendix, we use $\boldsymbol{\Sigma}_{F_{g}^{0}}$ and
$\boldsymbol{\Sigma}_{\Gamma_{g}^{0}}$ to denote the sub-matrices of
$\boldsymbol{\Sigma}_{F^{0}}$ and $\boldsymbol{\Sigma}_{\Gamma^{0}}$
corresponding to $T^{-\nu_{g}}\mathbf{F}_{g}^{0\prime}\mathbf{F}_{g}^{0}$ and
$N^{-1}\boldsymbol{\Gamma}_{g}^{0\prime}\boldsymbol{\Gamma}_{g}^{0}$,
respectively, for $g=1,\ldots,G$. We also define
$\mathbf{V}_{g}=\mathrm{diag}(\widehat{\lambda}_{g,1},\ldots,\widehat{\lambda}_{g,d_{max}})$,
where $\widehat{\lambda}_{g,d}$ has been defined in Step 2 of the IPC
estimation procedure. We partition
$\mathbf{F}^{0}=(\mathbf{F}_{1}^{0},\ldots,\mathbf{F}_{g}^{0},\mathbf{F}_{+g}^{0})$
and
$\mathbf{C}_{T}=\mathrm{diag}(T^{-\nu_{1}/2}\mathbf{I}_{d_{1}},\ldots,T^{-\nu_{g}/2}\mathbf{I}_{d_{g}},\mathbf{C}_{+g,T})$,
where $\mathbf{F}_{+g}^{0}=(\mathbf{F}_{g+1}^{0},\ldots,\mathbf{F}_{G}^{0})$
is $T\times(d_{g+1}+\cdots+d_{G})$, and
$\mathbf{C}_{+g,T}=\mathrm{diag}(T^{-\nu_{g+1}/2}\mathbf{I}_{d_{g+1}},\ldots,T^{-\nu_{G}/2}\mathbf{I}_{d_{G}})$
is $(d_{g+1}+\cdots+d_{G})\times(d_{g+1}+\cdots+d_{G})$. We partition
$\widehat{\mathbf{F}}$, $\boldsymbol{\gamma}_{i}^{0}$ and
$\boldsymbol{\Gamma}^{0}$ conformably as
$\widehat{\mathbf{F}}=(\widehat{\mathbf{F}}_{1},\ldots,\widehat{\mathbf{F}}_{g},\widehat{\mathbf{F}}_{+g})$,
$\boldsymbol{\gamma}_{i}^{0}=(\boldsymbol{\gamma}_{1,i}^{0\prime},\ldots,\boldsymbol{\gamma}_{g,i}^{0\prime},\boldsymbol{\gamma}_{+g,i}^{0\prime})^{\prime}$
and
$\boldsymbol{\Gamma}^{0}=(\boldsymbol{\Gamma}_{1}^{0},\ldots,\boldsymbol{\Gamma}_{g}^{0},\boldsymbol{\Gamma}_{+g}^{0})$,
respectively.

We introduce
$\lambda_{g,d}=T^{-\nu_{g}}\mathbf{h}_{g,d}^{0\prime}\mathbf{F}_{g}^{0\prime}\boldsymbol{\Sigma}_{g}^{0}\mathbf{F}_{g}^{0}\mathbf{h}_{g,d}^{0}$,
where
$\boldsymbol{\Sigma}_{g}^{0}=N^{-1}\mathbf{F}_{g}^{0}\boldsymbol{\Gamma}_{g}^{0\prime}\boldsymbol{\Gamma}_{g}^{0}\mathbf{F}_{g}^{0\prime}$
and $\mathbf{h}_{g,d}^{0}$ is the $d$-th column of
$\mathbf{H}_{g}^{0}=N^{-1}T^{(\nu_{g}-\delta)/2}\boldsymbol{\Gamma}_{g}^{0\prime}\boldsymbol{\Gamma}_{g}^{0}\mathbf{F}_{g}^{0\prime}\widehat{\mathbf{F}}_{g}^{0}(\mathbf{V}_{g}^{0})^{-1}$
with $\widehat{\mathbf{F}}_{g}^{0}$ being the $T\times d_{g}$ matrix
consisting of the first $d_{g}$ columns of $\widehat{\mathbf{F}}_{g}$ and
$\mathbf{V}_{g}^{0}$ being the leading $d_{g}\times d_{g}$ principal submatrix
of $\mathbf{V}_{g}$. In other words, $\widehat{\mathbf{F}}_{g}^{0}$ and
$\mathbf{V}_{g}^{0}$ are $\widehat{\mathbf{F}}_{g}$ and $\mathbf{V}_{g}$ based
on treating the number of factors for each group $g$, $d_{g}$, as known. We
also define $\mathbf{H}_{g}=T^{-(\nu_{g}-\delta)/2}\mathbf{H}_{g}^{0}$. In
order to appreciate the implication of the difference in normalization with
respect to $T$, let us consider $\mathbf{H}_{g}^{0}$. By Assumption 1,
$N^{-1}\boldsymbol{\Gamma}_{g}^{0\prime}\boldsymbol{\Gamma}_{g}^{0}$ is
asymptotically of full rank, and hence
$\|N^{-1}\boldsymbol{\Gamma}_{g}^{0\prime}\boldsymbol{\Gamma}_{g}^{0}\|=O_{P}(1)$.
Hence, since

$\displaystyle\|T^{-(\nu_{g}+\delta)/2}\mathbf{F}_{g}^{0\prime}\widehat{\mathbf{F}}_{g}\|^{2}\leq
T^{-\delta}\|\widehat{\mathbf{F}}_{g}\|_{2}^{2}T^{-\nu_{g}}\|\mathbf{F}_{g}^{0}\|^{2}=T^{-\nu_{g}}\|\mathbf{F}_{g}^{0}\|^{2}=O_{P}\left(1\right).$
(A.1)

and $\|(T^{-\nu_{1}}\mathbf{V}_{1}^{0})^{-1}\|=O_{P}(1)$ as explained under
(A.28), we can show that

$\displaystyle\|\mathbf{H}_{g}^{0}\|\leq\|N^{-1}\boldsymbol{\Gamma}_{g}^{0\prime}\boldsymbol{\Gamma}_{g}^{0}\|\|T^{-(\nu_{g}+\delta)/2}\mathbf{F}_{g}^{0\prime}\widehat{\mathbf{F}}_{g}^{0}\|\|(T^{-\nu_{g}}\mathbf{V}_{g}^{0})^{-1}\|=O_{P}(1),$
(A.2)

which in turn implies

$\displaystyle\|\mathbf{H}_{g}\|=T^{-(\nu_{g}-\delta)/2}\|\mathbf{H}_{g}^{0}\|=O_{P}(T^{-(\nu_{g}-\delta)/2}).$
(A.3)

We further use $\widehat{\mathbf{F}}_{g,d}^{0}$ to refer to the $d$-th column
of $\widehat{\mathbf{F}}_{g}^{0}$. In this notation,
$\widehat{\lambda}_{g,d}=T^{-\delta}\widehat{\mathbf{F}}_{g,d}^{0\prime}\widehat{\boldsymbol{\Sigma}}_{g}\widehat{\mathbf{F}}_{g,d}^{0}$
for $d=1,\ldots,d_{g}$.

We also partition $\mathbf{X}_{i}$ as
$\mathbf{X}_{i}=(\mathbf{X}_{1,i},\ldots,\mathbf{X}_{d_{x},i})$ with
$\mathbf{X}_{j,i}$ being the $j$-th column of $\mathbf{X}_{i}$. The $j$-th
column of $\mathbf{X}_{i}\mathbf{D}_{T}$ is therefore given by
$T^{-\kappa_{j}/2}\mathbf{X}_{j,i}$. Moreover, $\mathrm{vec}\,\mathbf{A}$,
$\mathrm{rank}\,\mathbf{A}$, $\mathrm{span}\,\mathbf{A}$ and
$\lambda(\mathbf{A})$ denote the vectorized version, rank, span and
eigenvalues of $\mathbf{A}$, respectively, $a\wedge b=\min\\{a,b\\}$ and
$a\vee b=\max\\{a,b\\}$.

##  A.2 Conditions that ensure asymptotic unbiasedness

In this section, we provide a set of assumptions that ensure that the
asymptotic distribution of
$\sqrt{NT}\mathbf{D}_{T}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})$
given in Theorem 3.1 is free of bias without for that matter requiring that
$\varepsilon_{i,t}$ is serially and cross-sectionally independent. One way to
accomplish this is to assume that $\rho_{1}=\rho_{2}=0$, as in Corollary 3.1.
The assumptions considered here, which are stated in Assumption A.1, can be
seen as alternatives to this last condition. In terms of the notation of
Theorem 3.1, they ensure that
$\mathbf{A}_{1}=\mathbf{A}_{2}=\mathbf{0}_{d_{x}\times 1}$.

######  Assumption A.1 (No asymptotic bias).

One of the following set of conditions is met:

  * (a)

$T^{1-\nu_{g}}E(\mathbf{f}_{g,t}^{0\prime}\mathbf{f}_{g,s}^{0}|\mathbf{x}_{t},\mathbf{x}_{s})=\phi_{ts}$
w.p.a.1 and $\sum_{t=1}^{T}\sum_{s=1}^{T}|\phi_{ts}|=O(T)$, where
$\mathbf{x}_{t}=(\mathbf{x}_{1,t},\ldots,\mathbf{x}_{N,t})^{\prime}$. If
$G\geq 2$, then $q<(\nu_{G}+\nu_{G-1})/2-1/4$ and $\nu_{g-1}-\nu_{g}>1/2$ for
$g=2,\ldots,G$.

  * (b)

$T/N\to c_{4}\in(0,\infty)$ and $\nu_{G}>1$. If $G\geq 2$, then
$q<(\nu_{G}+\nu_{G-1}-1)/2$ and $\nu_{g-1}-\nu_{g}>1/2$.

  * (c)

$T^{1-\nu_{g}-\kappa_{j}}E(\sum_{t=1}^{T}\sum_{s=1}^{T}x_{j,i,t}x_{j,k,s}\mathbf{f}_{g,t}^{0\prime}\mathbf{f}_{g,s}^{0})=O(T^{2-r})$,
where $r<2$, $r+\nu_{G}-1>0$ and $x_{j,i,t}$ is the $j$-th row of
$\mathbf{x}_{i,t}$. If $G\geq 2$, then $\nu_{g-1}-\nu_{g}>1/2$.

Assumption A.1 ensures that the asymptotic distribution of
$\sqrt{NT}\mathbf{D}_{T}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})$
is bias-free. It is, however, not necessary, and (a)–(b) should therefore be
viewed as examples of conditions under which there is no asymptotic
distribution bias. These conditions all have their strengths and weaknesses,
and so their suitability will in general depend on the context. Take as an
example condition (b), which has the advantage of not requiring any more
moment conditions than those that are already in Assumption 1. It does,
however, require that $\nu_{G}>1$, which rules out both signal-weak and
stationary factors ($\nu_{G}\leq 1$). Conditions (a) and (c) are more general
in this regard, but then at the expense of requiring additional moment
conditions. Condition (c) is a functional central limit theorem style moment
condition.

The next corollary to Theorem 3.1 verifies that the asymptotic distribution of
$\sqrt{NT}\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})$
is indeed bias-free under Assumption A.1.

######  Corollary A.1 (Unbiased asymptotic distribution).

Suppose that Assumptions 1–5, 7, and A.1 are met and that $N/T^{\nu_{G}}\to
0$. Then, as $N,\,T\to\infty$,

$\displaystyle\sqrt{NT}\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})\to_{D}MN(\mathbf{0}_{d_{x}\times
1},\mathbf{B}_{0}^{-1}\boldsymbol{\Omega}\mathbf{B}_{0}^{-1}).$ (A.1)

The asymptotic distribution in Corollary A.1 is the same as the one given in
Corollary 3.1.

If $\rho_{1}$, $\rho_{2}$, $\mathbf{A}_{1}$ and $\mathbf{A}_{2}$ are all
different from zero, one possibility is to use bias correction. Dhaene and
Jochmans (2015) were first to bring attention to the relevance of the half-
panel Jackknife approach for bias correction in panel data. These authors
focus on the fixed effects case, but Chen et al. (2021) have shown that the
Jackknife can be used to correct for bias also in models with interactive
effects. In our setting, the value of $\nu_{G}$ is unknown, and therefore the
standard Jackknife is not directly applicable. We thus turn to the hybrid
half-panel Jackknife correction proposed by Fernández-Val and Weidner (2018).
The bias-corrected version of the IPC estimator is given by

$\displaystyle\widehat{\boldsymbol{\beta}}_{\text{bc}}=3\widehat{\boldsymbol{\beta}}-\frac{1}{2}(\widehat{\boldsymbol{\beta}}_{1,N}+\widehat{\boldsymbol{\beta}}_{2,N}+\widehat{\boldsymbol{\beta}}_{1,T}+\widehat{\boldsymbol{\beta}}_{2,T}),$
(A.2)

where $\widehat{\boldsymbol{\beta}}$ is the standard IPC estimator defined in
(2.10), $\widehat{\boldsymbol{\beta}}_{1,N}$ and
$\widehat{\boldsymbol{\beta}}_{2,N}$ are defined in the same way but applied
to cross-sectional units $\\{1,\ldots,\lfloor N/2\rfloor\\}$ and $\\{\lfloor
N/2\rfloor+1,\ldots,N\\}$, respectively, and
$\widehat{\boldsymbol{\beta}}_{1,T}$ and $\widehat{\boldsymbol{\beta}}_{2,T}$
are the IPC estimators applied to odd and even numbered time periods,
respectively. The splitting in even and odd numbered time periods is needed in
order to account for the behaviour of the factors. See Fernández-Val and
Weidner (2018) for a more detailed discussion of the hybrid Jackknife
approach.

##  A.3 Extra empirical study — Returns to scale

In this section, we estimate the returns to scale (RTS) of US banks, an area
that has attracted considerable attention in the empirical production
literature (see, for example, Feng and Serletis, 2008). Let us therefore
denote by $CT_{i,t}$ the observed total cost of bank $i$ during period $t$.
The total cost is a function $C(\cdot)$ of a $J\times 1$ vector of input
prices $\mathbf{p}_{i,t}=(p_{1,i,t},\ldots,p_{J,i,t})^{\prime}$ and a $L\times
1$ vector of output quantities
$\mathbf{q}_{i,t}=(q_{1,i,t},\ldots,q_{L,i,t})^{\prime}$. Hence,
$CT_{i,t}=C(\mathbf{p}_{i,t},\mathbf{q}_{i,t})$, or when expressed in logs,
$\ln CT_{i,t}=\ln C(\mathbf{p}_{i,t},\mathbf{q}_{i,t})$. The RTS is simply the
reciprocal of the sum of the output elasticities, which can be estimated for a
given $C(\cdot)$. In order to ensure that $C(\cdot)$ is linearly homogenous in
prices, however, it is very common to first normalize costs and prices by one
of the prices, $p_{J,i,t}$ say. Hence, $\ln(CT_{i,t}/p_{J,i,t})=\ln
C(\mathbf{p}_{i,t}/p_{J,i,t},\mathbf{q}_{i,t})$, which in turn implies that
the RTS is given by

$\displaystyle\mathrm{RTS}=\left(\sum_{j=1}^{L}\frac{\partial\ln
C(\mathbf{p}_{i,t}/p_{J,i,t},\mathbf{q}_{i,t})}{\partial\ln
q_{j,i,t}}\right)^{-1}.$ (A.1)

In order to estimate this quantity, estimates of the partial derivatives are
needed. The standard way to obtain such estimates in the literature is to
assume that $\ln C(\mathbf{p}_{i,t}/p_{J,i,t},\mathbf{q}_{i,t})$ is linear, to
estimate the resulting linear relationship between $\ln(CT_{i,t}/p_{J,i,t})$,
$\ln(\mathbf{p}_{i,t}/p_{J,i,t})$ and $\ln\mathbf{q}_{i,t}$ by OLS, and to
estimate the partial derivatives by simply replacing the parameters by their
OLS estimates. A number of modelling issues then arise. First, as pointed out
by Feng and Zhang (2012, page 1884), “[t]here are many … types of … unobserved
heterogeneity in the US banking industry, which include, but are not limited
to, location-related corporate income tax rate, property tax rate, and
personal income tax rate.” Hence, there is a need to control for unobserved
heterogeneity. Second, the literature typically assumes that the regressors in
$\ln(\mathbf{p}_{i,t}/p_{J,i,t})$ and $\ln\mathbf{q}_{i,t}$ are stationary,
which, as pointed out by Dong et al. (2021), is highly unlikely to be the case
in practice. There is therefore a need to consider approaches that do not
require the regressors to be stationary. Third, to account for the fact that
costs are typically trending, many researchers fit their models with linear
and quadratic trend terms. Of course, deterministic trends can account for
some trending behaviour, but not all. Moreover, the results tend to be highly
sensitive to how the trending is modelled, suggesting that this is an
important issue (see Feng and Serletis, 2008, for a discussion).

The discussion of the last paragraph suggests that there is a need for an
approach that is general enough to accommodate not only unobserved
heterogeneity but also trending behaviour. The proposed IPC approach fits this
bill and we will therefore use it in this empirical illustration to the RTS of
US commercial banks. The data that we will use for this purpose, which are the
same as in Feng et al. (2017), are quarterly and cover the period 1986–2005,
which means that $T=80$. To avoid the impact of entry and exit, we only
include continuously operating large banks with assets of at least one billion
USD. This gives us a total of $N=466$ banks. The data set contains
observations on three input prices and output quantities. They are the wage
rate of labour ($p_{1,i,t}$), the price of deposits and purchased funds
($p_{2,i,t}$), the price of physical capital ($p_{3,i,t}$), consumer loans
($q_{1,i,t}$), non-consumer loans, which is composed of industrial, commercial
and real estate loans ($q_{2,i,t}$), and securities, including non-loan
financial assets ($q_{3,i,t}$). All outputs are deflated by the GDP deflator.
Following the previous literature (see, for example, Stiroh, 2000), total cost
($CT_{i,t}$) is computed as the sum of total salaries and benefits divided by
the number of full-time employees, the price of deposits and purchased funds
equals total interest expense divided by total deposits and purchased funds,
and the price of capital equals expenses on premises and equipment divided by
premises and fixed assets. Hence, after normalization by $p_{3,i,t}$, there
are six variables included in the model; $\ln(CT_{i,t}/p_{3,i,t})$,
$\ln(p_{1,i,t}/p_{3,i,t})$, $\ln(p_{2,i,t}/p_{3,i,t})$, $\ln q_{1,i,t}$, $\ln
q_{2,i,t}$ and $\ln q_{3,i,t}$.

Figure A.1: The variables in the RTS illustration.

In order to get a feeling for the trending behaviour of the data, in Figure
A.1 we plot all $N$ series for each variable. A few observations are
noteworthy. First, there is strong co-movement among the series. This is true
for all six variables, including the dependent variable,
$\ln(CT_{i,t}/p_{3,i,t})$, which we take as evidence in support of our
interactive effects specification. Second, the variables are highly persistent
and most are trending over time, suggesting that they are non-stationary. This
last observation is supported by a formal augmented Dickey–Fuller (ADF) test
(with intercept and linear trend) applied to each series. The highest
rejection frequency across all $N$ series is obtained for $\ln q_{2,i,t}$ and
is given by 6.7%, which in turn suggests that the evidence against unit root
null hypothesis is weak. In fact, since the ADF test does not account for the
multiplicity of the testing problem, the true proportion of (trend-)stationary
series is likely to be even lower. Third, it is unclear if the trend in
$\ln(CT_{i,t}/p_{3,i,t})$ is the same as those in the regressors, which means
that it is important to allow the factors to be trending, as otherwise any
unaccounted for trend in $\ln(CT_{i,t}/p_{3,i,t})$ will be pushed into the
regression errors.

Figure A.2: The estimated factors in the RTS illustration.

We now go on to discuss the actual estimation results, which are reported in
Table A.1. The estimated model is given by (2.1) with
$y_{i,t}=\ln(CT_{i,t}/p_{3,i,t})$,
$\boldsymbol{\beta}^{0}=(\beta_{1}^{0},\ldots,\beta_{5}^{0})^{\prime}$ and
$\mathbf{x}_{i,t}=(\ln(p_{1,i,t}/p_{3,i,t}),\ln(p_{2,i,t}/p_{3,i,t}),\ln
q_{1,i,t},\ln q_{2,i,t},\ln q_{3,i,t})^{\prime}$. Hence, in this notation,

$\displaystyle\mathrm{RTS}=\frac{1}{\beta_{3}^{0}+\beta_{4}^{0}+\beta_{5}^{0}}$
(A.2)

is the object of interest, which we estimate using

$\displaystyle\widehat{\mathrm{RTS}}=\frac{1}{\widehat{\beta}_{3}+\widehat{\beta}_{4}+\widehat{\beta}_{5}},$
(A.3)

where $\widehat{\beta}_{3}$, $\widehat{\beta}_{4}$ and $\widehat{\beta}_{5}$
are from
$\widehat{\boldsymbol{\beta}}=(\widehat{\beta}_{1},\ldots,\widehat{\beta}_{5})^{\prime}$.
The values of $\delta$ and $d_{max}$ are set as in Section 4. The value of
$\delta$ is again arbitrary and does not require justification. As for
$d_{max}$, 10 factors should be large enough to capture the true number, yet
not “excessively” large (see Ahn and Horenstein, 2013, for a discussion). In
order to enable inference not only for $\boldsymbol{\beta}^{0}$ but also for
the RTS, we follow Feng et al. (2017) and report 95% bootstrap confidence
intervals. The estimated factor groups are given by
$\widehat{d}_{1}=\widehat{d}_{2}=1$, implying that the estimated number of
factors equals $\widehat{d}_{f}=\widehat{d}_{1}+\widehat{d}_{2}=2$. In Figure
A.2, we plot the estimated factors, denoted $\widehat{f}_{1,t}$ and
$\widehat{f}_{2,t}$. We see that while both factor estimates are trending,
$\widehat{f}_{2,t}$ is growing much faster. As a measure of this difference in
the trends, we look at
$\sum_{i=1}^{N}\widehat{\gamma}_{1,i}^{2}/\sum_{i=1}^{N}\widehat{\gamma}_{2,i}^{2}$.
By using the results provided in the appendix, we can show that this ratio
should be $O_{P}(T^{\nu_{1}-\nu_{2}})$, implying
$\ln(\sum_{i=1}^{N}\widehat{\gamma}_{1,i}^{2}/\sum_{i=1}^{N}\widehat{\gamma}_{2,i}^{2})/\ln
T=\nu_{1}-\nu_{2}+o_{P}(1)$. By plugging in the known value of $T$ and
$\sum_{i=1}^{N}\widehat{\gamma}_{1,i}^{2}/\sum_{i=1}^{N}\widehat{\gamma}_{2,i}^{2}=123.4$,
we get
$\ln(\sum_{i=1}^{N}\widehat{\gamma}_{1,i}^{2}/\sum_{i=1}^{N}\widehat{\gamma}_{2,i}^{2})/\ln
T=1.1$, which is thus an estimate of $\nu_{1}-\nu_{2}$. Hence, there is indeed
a substantial difference in the degree of trending of the two factors. But it
is not as large as it would have been if $\widehat{f}_{1,t}$ and
$\widehat{f}_{2,t}$ were made up of a quadratic and a linear trend, say, in
which case $\nu_{1}=5$ and $\nu_{2}=3$, such that $\nu_{1}-\nu_{2}=2$. Hence,
just as expected given Figure A.2, while trending, the behaviour of the
factors is not just a simple deterministic function of time, which of course
casts doubt on much of the existing research.

While our main interest is in the RTS, for completeness in Table A.1, we also
report the estimated elasticities. The first thing to note is that the RTS is
estimated to be significantly larger than one, suggesting that the large US
banks considered here have had increasing returns to scale during the period
of investigation. This is consistent with the results reported by, for
example, Feng et al. (2017) and Wheelock and Wilson (2012). One explanation
for this finding is that the banks have been adopting productivity-enhancing
internet technologies, leading to increased RTS. The confidence intervals are
very narrow not only for the RTS, but also for the elasticities, which are all
significant.

Table A.1: Empirical results for the RTS illustration.

Estimator | Point estimate |  $p$-value
---|---|---
$\widehat{\beta}_{1}$ | 0.4063 | 0.0000
$\widehat{\beta}_{2}$ | 0.0199 | 0.0000
$\widehat{\beta}_{3}$ | 0.3805 | 0.0000
$\widehat{\beta}_{4}$ | 0.0912 | 0.0000
$\widehat{\beta}_{5}$ | 0.4726 | 0.0000
$\widehat{\mathrm{RTS}}$ | 1.0590 | 0.0000

  * •

_Notes_ : $\widehat{\beta}_{1}$, $\widehat{\beta}_{2}$, $\widehat{\beta}_{3}$,
$\widehat{\beta}_{4}$ and $\widehat{\beta}_{5}$ refer to the estimated IPC
slopes of $\ln(p_{i,t}/p_{3,i,t})$, $\ln(p_{2,i,t}/p_{3,i,t})$, $\ln
q_{1,i,t}$, $\ln q_{2,i,t}$ and $\ln q_{3,i,t}$, respectively, when the
dependent variable is $\ln(CT_{i,t}/p_{3,i,t})$. The reported Wald $p$-values
for $\widehat{\beta}_{1}$–$\widehat{\beta}_{5}$ test the null hypothesis that
the relevant coefficient is zero. The $p$-value for RTS test the null
hypothesis that $\text{RTS}=1$.

##  A.4 Outline of the proofs

In this section, we describe the outline of the proofs. We assume throughout
that $G\geq 2$. The proofs for the cases when $G\in\\{0,1\\}$ are much
simpler, and can be obtained by manipulating the proofs for $G\geq 2$.

Lemma A.1 is a widely used result for studying the eigenvalues of large
dimensional matrices (Lam et al., 2011), and is presented here for
convenience.

Our own theoretical development is carried out stepwise, starting with Step 1
of the IPC estimation procedure. Lemma A.2 is very important in this regard.
It presents two order results that are used repeatedly throughout the proofs.
Given Lemma A.2, we are able to establish Lemma 3.1 of the main text, which
provides a lower bound on the rate of convergence of the initial Step 1
estimator $\widehat{\boldsymbol{\beta}}_{0}$.

As a first step towards establishing the consistency of the Step 2 estimators
of $(d_{2},\ldots,d_{G})$, in Lemma A.3 we study limiting behaviour of the
eigenvalues associated with $\widehat{\boldsymbol{\Sigma}}_{1}$ in (2.7). The
consistency of $\widehat{d}_{1}$ reported in Lemma 3.2 of the main text is a
direct consequence of Lemma A.3. Lemmas A.4 and A.5 enable consistent
estimation of $(d_{2},\ldots,d_{G})$, and are thus key in proving Lemma 3.3 of
the main text. For simplicity, in the proofs of Lemmas A.4 and A.5 we treat
$d_{1}$ as known, which according to Lemma 3.2 is justified in large samples.

The final step in our theoretical analysis is to prove the asymptotic mixed
normality of the Step 3 IPC estimator $\widehat{\boldsymbol{\beta}}$, as
stated in Theorem 3.1. In so doing, we treat $(d_{1},\ldots,d_{G})$ as known,
which is again justified in large samples. We begin by proving Lemma A.6,
which is needed in the proof of Theorem 3.1. The proofs of Corollaries 3.1 and
A.1 are immediate consequences of that of Theorem 3.1. We also derive the rate
of convergence of $\widehat{\boldsymbol{\beta}}_{0}$ under the conditions of
Theorem 3.1. This rate is provided as a part of Lemma A.6.

##  A.5 Auxiliary lemmas

######  Lemma A.1.

Suppose that $\mathbf{A}$ and $\mathbf{A}+\mathbf{E}$ are $n\times n$
symmetric matrices and that $\mathbf{Q}=(\mathbf{Q}_{1},\mathbf{Q}_{2})$,
where $\mathbf{Q}_{1}$ is $n\times r$ and $\mathbf{Q}_{2}$ is $n\times(n-r)$,
is an orthogonal matrix such that $\text{span}\,\mathbf{Q}_{1}$ is an
invariant subspace for $\mathbf{A}$. Decompose
$\mathbf{Q}^{\prime}\mathbf{A}\mathbf{Q}$ and
$\mathbf{Q}^{\prime}\mathbf{E}\mathbf{Q}$ as
$\mathbf{Q}^{\prime}\mathbf{A}\mathbf{Q}=\mathrm{diag}(\mathbf{D}_{1},\mathbf{D}_{2})$
and

$\displaystyle\mathbf{Q}^{\prime}\mathbf{E}\mathbf{Q}=\left(\begin{array}[]{cc}\mathbf{E}_{11}&\mathbf{E}_{12}\\\
\mathbf{E}_{21}&\mathbf{E}_{22}\end{array}\right).$ (A.3)

Let
$\mathrm{sep}(\mathbf{D}_{1},\mathbf{D}_{2})=\min_{\lambda_{1}\in\lambda(\mathbf{D}_{1}),\
\lambda_{2}\in\lambda(\mathbf{D}_{2})}|\lambda_{1}-\lambda_{2}|$. If
$\mathrm{sep}(\mathbf{D}_{1},\mathbf{D}_{2})>0$ and
$\|\mathbf{E}\|_{2}\leq\mathrm{sep}(\mathbf{D}_{1},\mathbf{D}_{2})/5$, then
there exists a $(n-r)\times r$ matrix $\mathbf{P}$ with
$\|\mathbf{P}\|_{2}\leq
4\|\mathbf{E}_{21}\|_{2}/\mathrm{sep}(\mathbf{D}_{1},\mathbf{D}_{2})$, such
that the columns of
$\mathbf{Q}_{1}^{0}=(\mathbf{Q}_{1}+\mathbf{Q}_{2}\mathbf{P})(\mathbf{I}_{r}+\mathbf{P}^{\prime}\mathbf{P})^{-1/2}$
define an orthonormal basis for a subspace that is invariant for
$\mathbf{A}+\mathbf{E}$.

######  Lemma A.2.

Under Assumption 1, as $N,\,T\to\infty$,

  * $\mathrm{(a)}$

$\sup_{\mathbf{F}\in\mathbb{D}_{F}}(NT)^{-1}\sum_{i=1}^{N}\boldsymbol{\varepsilon}_{i}^{\prime}\mathbf{P}_{F}\boldsymbol{\varepsilon}_{i}=O_{P}(N^{-1}\vee
T^{-1})$;

  * $\mathrm{(b)}$

$\sup_{\mathbf{F}\in\mathbb{D}_{F}}\|(NT)^{-1}\sum_{i=1}^{N}\mathbf{D}_{T}\mathbf{X}_{i}^{\prime}\mathbf{P}_{F}\boldsymbol{\varepsilon}_{i}\|=O_{P}(N^{-1/2}\vee
T^{-1/2})$.

If Assumption 3 also holds, then

  * $\mathrm{(c)}$

$(NT)^{-1}\|\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^{\prime}\|=O_{P}(N^{-1/2}\vee
T^{-1/2})$ and
$(NT)^{-1}\|\boldsymbol{\varepsilon}^{\prime}\boldsymbol{\varepsilon}\|=O_{P}(N^{-1/2}\vee
T^{-1/2})$;

  * $\mathrm{(d)}$

$\|\boldsymbol{\Gamma}^{0\prime}\boldsymbol{\varepsilon}\|=O_{P}(\sqrt{NT})$.

######  Lemma A.3.

Let Assumptions 1–4 hold. Then, as $N,\,T\to\infty$,

  * $\mathrm{(a)}$

$T^{-\nu_{1}}|\widehat{\lambda}_{1,d}-\lambda_{1,d}|=O_{P}(T^{-(\nu_{1}-\nu_{2})/2})$
for $d=1,\ldots,d_{1}$;

  * $\mathrm{(b)}$

$T^{-\nu_{1}}|\widehat{\lambda}_{1,d}|=O_{P}(T^{-(\nu_{1}-\nu_{2})})$ for
$d=d_{1}+1,\ldots,d_{max}$.

######  Lemma A.4.

As $N,\,T\to\infty$, the following hold under Assumptions 1–5:

  * $\mathrm{(a)}$

$T^{-\delta}\|\mathbf{F}_{2}^{0\prime}\widehat{\mathbf{F}}_{1}\|=O_{P}(T^{-(\delta+\nu_{1}-\nu_{2})/2}+N^{-1/2}T^{-(\delta+\nu_{1}-\nu_{2}-1)/2}+N^{-(1-p)}T^{-(\delta+\nu_{1}-2\nu_{2})/2}+T^{-(\nu_{1}+\delta)/2-q})$;

  * $\mathrm{(b)}$

$\sum_{i=1}^{N}\|\mathbf{F}_{1}^{0\prime}\boldsymbol{\gamma}_{1,i}^{0}-\widehat{\mathbf{F}}_{1}\widehat{\boldsymbol{\gamma}}_{1,i}\|^{2}=O_{P}(N\vee
T+NT^{2q-\nu_{1}}+N^{-(1-2p)}T^{\nu_{2}})$.

######  Lemma A.5.

Let
$\tau_{NT}=N^{-1/2}T^{-(\nu_{2}-1)/2}+T^{-\nu_{2}/2}+T^{-(\nu_{2}-\nu_{3})}+T^{q-(\nu_{1}+\nu_{2})/2}+N^{-(1-p)}$.
As $N,\,T\to\infty$, the following hold under the conditions of Lemma A.4:

  * $\mathrm{(a)}$

$T^{-\nu_{2}}|\widehat{\lambda}_{2,d}-\lambda_{2,d}|=O_{P}(\tau_{NT})$ for
$d=1,\ldots,d_{2}$;

  * $\mathrm{(b)}$

$T^{-\nu_{2}}|\widehat{\lambda}_{2,d}|=O_{P}(\tau_{NT}^{2})$ for
$d=d_{2}+1,\ldots,d_{max}$.

######  Lemma A.6.

Let Assumptions 1–5 hold. In addition, let $NT^{-\nu_{G}}<\infty$, as
$N,\,T\to\infty$,

  * $\mathrm{(a)}$

$\|\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}_{1}-\widehat{\boldsymbol{\beta}}_{0})\|=O_{P}((NT)^{-1/2}\vee\|\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}_{0}-\boldsymbol{\beta}^{0})\|)$;

  * $\mathrm{(b)}$

$\|\mathbf{D}_{T}^{-1}(\widehat{\boldsymbol{\beta}}_{0}-\boldsymbol{\beta}^{0})\|=O_{P}((NT)^{-1/2})$.

##  A.6 Proofs of auxiliary lemmas

Proof of Lemma A.1.

This is Lemma 3 of Lam et al. (2011). The proof is therefore
omitted.$\blacksquare$

Proof of Lemma A.2.

Consider (a). We have

$\displaystyle\sup_{\mathbf{F}\in\mathbb{D}_{F}}\frac{1}{NT}\sum_{i=1}^{N}\boldsymbol{\varepsilon}_{i}^{\prime}\mathbf{P}_{F}\boldsymbol{\varepsilon}_{i}$
$\displaystyle=\sup_{\mathbf{F}\in\mathbb{D}_{F}}(NT)^{-1}\text{tr}(\mathbf{P}_{F}\boldsymbol{\varepsilon}^{\prime}\boldsymbol{\varepsilon})\leq
O(1)\sup_{\mathbf{F}\in\mathbb{D}_{F}}(NT)^{-1}\|\mathbf{P}_{F}\|_{2}\|\boldsymbol{\varepsilon}^{\prime}\boldsymbol{\varepsilon}\|_{2}$
$\displaystyle\leq
O(1)\sup_{\mathbf{F}\in\mathbb{D}_{F}}(NT)^{-1}\|\boldsymbol{\varepsilon}\|_{2}^{2}=O_{P}(N^{-1}\vee
T^{-1}),$ (A.1)

where the first inequality follows from the fact that
$|\text{tr}\,\mathbf{A}|\leq\text{rank}\,\mathbf{A}\,\|\mathbf{A}\|_{2}$, the
second inequality follows from the fact that $\|\mathbf{P}_{{F}}\|_{2}=1$, and
the second equality hold by Assumption 1 (b).

The result in (b) is due to

$\displaystyle\sup_{\mathbf{F}\in\mathbb{D}_{F}}\left\|\frac{1}{NT}\sum_{i=1}^{N}\mathbf{D}_{T}\mathbf{X}_{i}^{\prime}\mathbf{P}_{F}\boldsymbol{\varepsilon}_{i}\right\|$
$\displaystyle\leq\frac{1}{NT}\sum_{j=1}^{d_{x}}\sup_{\mathbf{F}\in\mathbb{D}_{F}}\left|\text{tr}\,(T^{-\kappa_{j}/2}\mathbf{X}_{j}\mathbf{P}_{F}\boldsymbol{\varepsilon}^{\prime})\right|$
$\displaystyle\leq
O(1)\frac{1}{NT}\sum_{j=1}^{d_{x}}\sup_{\mathbf{F}\in\mathbb{D}_{F}}\|T^{-\kappa_{j}/2}\mathbf{X}_{j}\|_{2}\|\mathbf{P}_{F}\|_{2}\|\boldsymbol{\varepsilon}\|_{2}$
$\displaystyle=O(1)(NT)^{-1}O_{P}(\sqrt{NT})O_{P}(\sqrt{N}\vee\sqrt{T})$
$\displaystyle=O_{P}(N^{-1/2}\vee T^{-1/2}),$ (A.2)

where, with a slight abuse of notation and in this proof only,
$\mathbf{X}_{j}=(\mathbf{X}_{j,1},\ldots,\mathbf{X}_{j,N})^{\prime}$, the
second inequality follows from
$|\text{tr}\,\mathbf{A}|\leq\text{rank}\,\mathbf{A}\,\|\mathbf{A}\|_{2}$,
while the first equality is due to Assumption 1.

For (c), we use

$\displaystyle(NT)^{-2}E\|\boldsymbol{\varepsilon}^{\prime}\boldsymbol{\varepsilon}\|^{2}$
$\displaystyle=\frac{1}{(NT)^{2}}\sum_{t=1}^{T}\sum_{s=1}^{T}\left(\sum_{i=1}^{N}E[\varepsilon_{i,t}^{2}\varepsilon_{i,s}^{2}]+\sum_{i=1}^{N}\sum_{j\neq
i}E[\varepsilon_{i,t}\varepsilon_{i,s}\varepsilon_{j,t}\varepsilon_{j,s}]\right)$
$\displaystyle=\frac{1}{(NT)^{2}}\sum_{t=1}^{T}\left(\sum_{i=1}^{N}E(\varepsilon_{i,t}^{4})+\sum_{i=1}^{N}\sum_{j\neq
i}E[(\varepsilon_{i,t}\varepsilon_{j,t}-\sigma_{\varepsilon,ij})^{2}]\right)$
$\displaystyle+\frac{1}{(NT)^{2}}\sum_{t=1}^{T}\sum_{s\neq
t}\left(\sum_{i=1}^{N}E[\varepsilon_{i,t}^{2}\varepsilon_{i,s}^{2}]+\sum_{i=1}^{N}\sum_{j\neq
i}E[(\varepsilon_{i,t}\varepsilon_{j,t}-\sigma_{\varepsilon,ij})(\varepsilon_{i,s}\varepsilon_{j,s}-\sigma_{\varepsilon,ij})]\right)$
$\displaystyle+\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j\neq
i}\sigma_{\varepsilon,ij}^{2}=O(N^{-1})+O(T^{-1}),$ (A.3)

where the third equality follows from using the mixing condition on
$\varepsilon_{i,t}\varepsilon_{j,t}$ across $t$. The above result implies that
$(NT)^{-1}\|\boldsymbol{\varepsilon}^{\prime}\boldsymbol{\varepsilon}\|=O_{P}(N^{-1/2})+O_{P}(T^{-1/2})$
as required for the first result in (c). The second follows from

$\displaystyle(NT)^{-2}E\|\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^{\prime}\|^{2}$
$\displaystyle=\frac{1}{N^{2}T^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T}\sum_{s=1}^{T}E[\varepsilon_{i,t}\varepsilon_{j,t}\varepsilon_{i,s}\varepsilon_{j,s}]$
$\displaystyle=\frac{1}{(NT)^{2}}\sum_{t=1}^{T}\sum_{s=1}^{T}\left(\sum_{i=1}^{N}E[\varepsilon_{i,t}^{2}\varepsilon_{i,s}^{2}]+\sum_{i=1}^{N}\sum_{j\neq
i}E[\varepsilon_{i,t}\varepsilon_{i,s}\varepsilon_{j,t}\varepsilon_{j,s})\right)$
$\displaystyle=O(N^{-1})+O(T^{-1}),$ (A.4)

where the last step follows by the same arguments used to establish the first
result of (c).

It remains to prove (d), which is a direct consequence of Assumptions 1 and 3,
as seen from

$\displaystyle
E\|\boldsymbol{\Gamma}^{0\prime}\boldsymbol{\varepsilon}\|^{2}=\sum_{t=1}^{T}\sum_{i=1}^{N}\sum_{j=1}^{N}E[\boldsymbol{\gamma}_{i}^{0\prime}\boldsymbol{\gamma}_{j}^{0}\varepsilon_{i,t}\varepsilon_{j,t}]\leq
O(T)\sum_{i=1}^{N}\sum_{j=1}^{N}|\sigma_{\varepsilon,ij}|=O(NT).$ (A.5)

This establishes (d) and hence the proof of the lemma is
complete.$\blacksquare$

Proof of Lemma A.3.

Consider (a). As in Appendix A.1, decompose
$\mathbf{F}^{0}=(\mathbf{F}_{1}^{0},\mathbf{F}_{+1}^{0})$ and
$\mathbf{C}_{T}=\mathrm{diag}(T^{-\nu_{1}/2}\mathbf{I}_{d_{1}},\mathbf{C}_{+1,T})$,
where $\mathbf{F}_{+1}^{0}=(\mathbf{F}_{2}^{0},\ldots,\mathbf{F}_{G}^{0})$ is
$T\times(d_{f}-d_{1})$ and
$\mathbf{C}_{+1,T}=\mathrm{diag}(T^{-\nu_{2}/2}\mathbf{I}_{d_{2}},\ldots,T^{-\nu_{G}/2}\mathbf{I}_{d_{G}})$
is $(d_{f}-d_{1})\times(d_{f}-d_{1})$. We partition $\widehat{\mathbf{F}}$,
$\boldsymbol{\gamma}_{i}^{0}$ and $\boldsymbol{\Gamma}^{0}$ conformably as
$\widehat{\mathbf{F}}=(\widehat{\mathbf{F}}_{1},\widehat{\mathbf{F}}_{+1})$,
$\boldsymbol{\gamma}_{i}^{0}=(\boldsymbol{\gamma}_{1,i}^{0\prime},\boldsymbol{\gamma}_{+1,i}^{0\prime})^{\prime}$
and
$\boldsymbol{\Gamma}^{0}=(\boldsymbol{\Gamma}_{1}^{0},\boldsymbol{\Gamma}_{+1}^{0})$,
respectively.

By the definition of the eigenvectors and eigenvalues,
$\widehat{\boldsymbol{\Sigma}}_{1}\widehat{\mathbf{F}}_{1}=\widehat{\mathbf{F}}_{1}\mathbf{V}_{1}$.
By using this and the definition of $\widehat{\boldsymbol{\Sigma}}_{1}$,

$\displaystyle T^{-(\nu_{1}+\delta)/2}\widehat{\mathbf{F}}_{1}\mathbf{V}_{1}$
$\displaystyle=\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathbf{X}_{i}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})^{\prime}\mathbf{X}_{i}^{\prime}\widehat{\mathbf{F}}_{1}$
$\displaystyle+\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathbf{X}_{i}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\boldsymbol{\gamma}_{1,i}^{0\prime}\mathbf{F}_{1}^{0\prime}\widehat{\mathbf{F}}_{1}+\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathbf{F}_{1}^{0}\boldsymbol{\gamma}_{1,i}^{0}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})^{\prime}\mathbf{X}_{i}^{\prime}\widehat{\mathbf{F}}_{1}$
$\displaystyle+\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathbf{X}_{i}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})^{\prime}\widehat{\mathbf{F}}_{1}$
$\displaystyle+\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})^{\prime}\mathbf{X}_{i}^{\prime}\widehat{\mathbf{F}}_{1}$
$\displaystyle+\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})^{\prime}\widehat{\mathbf{F}}_{1}$
$\displaystyle+\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathbf{F}_{1}^{0}\boldsymbol{\gamma}_{1,i}^{0}(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})^{\prime}\widehat{\mathbf{F}}_{1}+\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})\boldsymbol{\gamma}_{1,i}^{0\prime}\mathbf{F}_{1}^{0\prime}\widehat{\mathbf{F}}_{1}$
$\displaystyle+\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathbf{F}_{1}^{0}\boldsymbol{\gamma}_{1,i}^{0}\boldsymbol{\gamma}_{1,i}^{0\prime}\mathbf{F}_{1}^{0\prime}\widehat{\mathbf{F}}_{1}$
$\displaystyle=\sum_{j=1}^{9}\mathbf{J}_{j},$ (A.6)

with implicit definitions of $\mathbf{J}_{1},\ldots,\mathbf{J}_{9}$. Note that

$\displaystyle\mathbf{J}_{9}=\mathbf{F}_{1}^{0}(N^{-1}\boldsymbol{\Gamma}_{1}^{0\prime}\boldsymbol{\Gamma}_{1}^{0})(T^{-(\nu_{1}+\delta)/2}\mathbf{F}_{1}^{0\prime}\widehat{\mathbf{F}}_{1}).$
(A.7)

Hence, moving this term over to the left-hand side, the above expression for
$T^{-(\nu_{1}+\delta)/2}\widehat{\mathbf{F}}_{1}\mathbf{V}_{1}$ becomes

$\displaystyle
T^{-(\nu_{1}+\delta)/2}\widehat{\mathbf{F}}_{1}\mathbf{V}_{1}-\mathbf{J}_{9}=\sum_{j=1}^{8}\mathbf{J}_{j}.$
(A.8)

We now evaluate each of the terms on the right-hand side.

Because $T^{-\delta}\|\widehat{\mathbf{F}}_{1}\|^{2}=d_{max}$ and
$(NT)^{-1}\sum_{i=1}^{N}\|\mathbf{X}_{i}\mathbf{D}_{T}\|^{2}=O_{p}(1)$ by
Assumption 1, the order of $\mathbf{J}_{1}$ is given by

$\displaystyle T^{-\delta/2}\|\mathbf{J}_{1}\|$
$\displaystyle\leq\left\|\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathbf{X}_{i}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})^{\prime}\mathbf{X}_{i}^{\prime}\right\|T^{-\delta/2}\|\widehat{\mathbf{F}}_{1}\|$
$\displaystyle\leq
O(1)\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\|\mathbf{X}_{i}\mathbf{D}_{T}\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|^{2}$
$\displaystyle=O_{P}(T^{1-(\nu_{1}+\delta)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|^{2}).$
(A.9)

Moreover, since

$\displaystyle\frac{1}{NT^{\nu_{1}}}\sum_{i=1}^{N}\|\mathbf{F}_{1}^{0}\boldsymbol{\gamma}_{1,i}^{0}\|^{2}\leq\frac{1}{N}\sum_{i=1}^{N}\|\boldsymbol{\gamma}_{1,i}^{0}\|^{2}T^{-\nu_{1}}\|\mathbf{F}_{1}^{0}\|^{2}=O_{P}(1)$
(A.10)

by Assumption 1, we can show that

$\displaystyle T^{-\delta/2}\|\mathbf{J}_{2}\|$ $\displaystyle\leq
O(1)\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\|\mathbf{X}_{i}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\boldsymbol{\gamma}_{1,i}^{0\prime}\mathbf{F}_{1}^{0\prime}\|$
$\displaystyle\leq
O(1)\left(\frac{1}{NT^{\delta}}\sum_{i=1}^{N}\|\mathbf{X}_{i}\mathbf{D}_{T}\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|^{2}\right)^{1/2}\left(\frac{1}{NT^{\nu_{1}}}\sum_{i=1}^{N}\|\mathbf{F}_{1}^{0}\boldsymbol{\gamma}_{1,i}^{0}\|^{2}\right)^{1/2}$
$\displaystyle=O_{P}(T^{(1-\delta)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|),$
(A.11)

and by exactly the same arguments,

$\displaystyle
T^{-\delta/2}\|\mathbf{J}_{3}\|=O_{P}(T^{(1-\delta)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|).$
(A.12)

For $\mathbf{J}_{4}$, we use

$\displaystyle\frac{1}{NT^{(\nu_{1}+\delta)/2}}\left\|\sum_{i=1}^{N}\mathbf{X}_{i}\mathbf{D}_{T}\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0})^{\prime}\right\|$
$\displaystyle\leq\left(\frac{1}{NT^{\delta}}\sum_{i=1}^{N}\|\mathbf{X}_{i}\mathbf{D}_{T}\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|^{2}\right)^{1/2}\left(\frac{1}{NT^{\nu_{1}}}\sum_{i=1}^{N}\|\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}\|^{2}\right)^{1/2}$
$\displaystyle\leq
O_{P}(T^{(1-\delta)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|)O_{P}(T^{-\nu_{1}/2}\|\mathbf{C}_{+1,T}^{-1}\|)$
$\displaystyle=O_{P}(T^{(1-\delta-\nu_{1}+\nu_{2})/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|),$
(A.13)

where the last equality makes use of the fact that
$\|\mathbf{C}_{+1,T}^{-1}\|=O(T^{\nu_{2}/2})$, as $\nu_{2}>\cdots>\nu_{G}$ by
Assumption 1. We can further show that

$\displaystyle\frac{1}{NT^{(\nu_{1}+\delta)/2}}\left\|\sum_{i=1}^{N}\mathbf{X}_{i}\mathbf{D}_{T}\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\boldsymbol{\varepsilon}_{i}^{\prime}\right\|$
$\displaystyle=\left\|\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathrm{vec}\,[\mathbf{X}_{i}\mathbf{D}_{T}\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\boldsymbol{\varepsilon}_{i}^{\prime}]\right\|$
$\displaystyle\leq\left\|\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}(\boldsymbol{\varepsilon}_{i}\otimes\mathbf{X}_{i}\mathbf{D}_{T})\right\|\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|$
$\displaystyle=O_{P}(N^{-1/2}T^{(2-\nu_{1}-\delta)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|),$
(A.14)

where the last equality holds, because by Assumption 3 and
$(\mathrm{tr}\,\mathbf{A}^{\prime}\mathbf{B})^{2}\leq(\mathrm{tr}\,\mathbf{A}^{\prime}\mathbf{A})(\mathrm{tr}\,\mathbf{B}^{\prime}\mathbf{B})$,
we have

$\displaystyle
E\left\|\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}(\boldsymbol{\varepsilon}_{i}\otimes\mathbf{X}_{i}\mathbf{D}_{T})\right\|^{2}$
$\displaystyle=\frac{1}{N^{2}T^{\nu_{1}+\delta}}\sum_{t=1}^{T}\sum_{s=1}^{T}\sum_{i=1}^{N}\sum_{j=1}^{N}E(\mathbf{x}_{j,t}^{\prime}\mathbf{D}_{T}^{2}\mathbf{x}_{i,t}\varepsilon_{i,s}\varepsilon_{j,s})\leq\frac{1}{N^{2}T^{\nu_{1}+\delta}}\sum_{i=1}^{N}\sum_{j=1}^{N}E|\mathrm{tr}(\mathbf{D}_{T}\mathbf{X}_{j}^{\prime}\mathbf{X}_{i}\mathbf{D}_{T})|\cdot|\sigma_{\varepsilon,ij}|$
$\displaystyle\leq\frac{1}{N^{2}T^{\nu_{1}+\delta-2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\sqrt{T^{-1}E\|\mathbf{D}_{T}\mathbf{X}_{j}\|^{2}}\sqrt{T^{-1}E\|\mathbf{X}_{i}\mathbf{D}_{T}\|^{2}}|\sigma_{\varepsilon,ij}|$
$\displaystyle=O(1)\frac{1}{N^{2}T^{\nu_{1}+\delta-2}}\sum_{i=1}^{N}\sum_{j=1}^{N}|\sigma_{\varepsilon,ij}|=O(N^{-1}T^{2-\nu_{1}-\delta}).$
(A.15)

Hence, by adding the results,

$\displaystyle T^{-\delta/2}\|\mathbf{J}_{4}\|$ $\displaystyle\leq
O(1)\frac{1}{NT^{(\nu_{1}+\delta)/2}}\left\|\sum_{i=1}^{N}\mathbf{X}_{i}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})^{\prime}\right\|$
$\displaystyle=O_{P}(T^{(1-\delta-\nu_{1}+\nu_{2})/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|+N^{-1/2}T^{(2-\nu_{1}-\delta)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|)$
$\displaystyle=o_{P}(T^{-(1-\delta)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|),$
(A.16)

where we have used Assumptions 1 and 4 ($T/N^{2}=O(1)$ under $\nu_{G}<1$) to
show that $T^{(1-\delta-\nu_{1}+\nu_{2})/2}$ and
$N^{-1/2}T^{(2-\nu_{1}-\delta)/2}$ are $o(1)$. The same steps can be used to
show that

$\displaystyle T^{-\delta/2}\|\mathbf{J}_{5}\|\leq
o_{P}(T^{-(1-\delta)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|).$
(A.17)

For $\mathbf{J}_{6}$, we use

$\displaystyle T^{-\delta/2}\|\mathbf{J}_{6}\|$ $\displaystyle\leq
O(1)\left\|\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})^{\prime}\right\|$
$\displaystyle\leq
O(1)\cdot(N^{-1}T^{-(\nu_{1}+\delta)/2}\|\mathbf{F}_{+1}^{0}\boldsymbol{\Gamma}_{+1}^{0\prime}\boldsymbol{\Gamma}_{+1}^{0}\mathbf{F}_{+1}^{0\prime}\|+N^{-1}T^{-(\nu_{1}+\delta)/2}\|\boldsymbol{\varepsilon}^{\prime}\boldsymbol{\varepsilon}\|$
$\displaystyle+2N^{-1}T^{-(\nu_{1}+\delta)/2}\|\mathbf{F}_{+1}^{0}\boldsymbol{\Gamma}_{+1}^{0\prime}\boldsymbol{\varepsilon}\|).$
(A.18)

By Assumption 1,

$\displaystyle
N^{-1}T^{-(\nu_{1}+\delta)/2}\|\mathbf{F}_{+1}^{0}\boldsymbol{\Gamma}_{+1}^{0\prime}\boldsymbol{\Gamma}_{+1}^{0}\mathbf{F}_{+1}^{0\prime}\|=O_{P}(T^{\nu_{2}-(\nu_{1}+\delta)/2}).$
(A.19)

Another application of Assumption 1 and Lemma A.2 gives

$\displaystyle
N^{-1}T^{-(\nu_{1}+\delta)/2}\|\boldsymbol{\varepsilon}^{\prime}\boldsymbol{\varepsilon}\|$
$\displaystyle=O_{P}(T^{1-(\nu_{1}+\delta)/2}(N^{-1/2}\vee T^{-1/2})),$ (A.20)
$\displaystyle
N^{-1}T^{-(\nu_{1}+\delta)/2}\|\mathbf{F}_{+1}^{0}\boldsymbol{\Gamma}_{+1}^{0\prime}\boldsymbol{\varepsilon}\|$
$\displaystyle\leq
N^{-1/2}T^{-(\nu_{1}+\delta-\nu_{2}-1)/2}T^{-\nu_{2}/2}\|\mathbf{F}_{+1}^{0}\|(NT)^{-1/2}\|\boldsymbol{\Gamma}_{+1}^{0\prime}\boldsymbol{\varepsilon}\|$
$\displaystyle=O_{P}(N^{-1/2}T^{-(\nu_{1}+\delta-\nu_{2}-1)/2}).$ (A.21)

These results can be inserted into the expression for
$T^{-\delta/2}\|\mathbf{J}_{6}\|$, giving

$\displaystyle
T^{-\delta/2}\|\mathbf{J}_{6}\|=O_{P}(T^{\nu_{2}-(\nu_{1}+\delta)/2}+N^{-1/2}T^{1-(\nu_{1}+\delta)/2}).$
(A.22)

Next up is $\mathbf{J}_{7}$. By using Assumptions 1 and 4, and Lemma A.2, and
the arguments use in evaluating $\mathbf{J}_{6}$,

$\displaystyle T^{-\delta/2}\|\mathbf{J}_{7}\|$ $\displaystyle\leq
O(1)\left\|\frac{1}{NT^{(\nu_{1}+\delta)/2}}\sum_{i=1}^{N}\mathbf{F}_{1}^{0}\boldsymbol{\gamma}_{1,i}^{0}(\mathbf{F}_{+1}^{0}\boldsymbol{\gamma}_{+1,i}^{0}+\boldsymbol{\varepsilon}_{i})^{\prime}\right\|$
$\displaystyle\leq
O(1)N^{-1}T^{-(\nu_{1}+\delta)/2}\|\mathbf{F}_{1}^{0}\boldsymbol{\Gamma}_{1}^{0\prime}\boldsymbol{\Gamma}_{+1}^{0}\mathbf{F}_{+1}^{0\prime}\|+O(1)N^{-1}T^{-(\nu_{1}+\delta)/2}\|\mathbf{F}_{1}^{0}\boldsymbol{\Gamma}_{1}^{0\prime}\boldsymbol{\varepsilon}\|$
$\displaystyle=O_{P}(T^{(\nu_{2}-\delta)/2}+N^{-1/2}T^{(1-\delta)/2})=O_{P}(T^{(\nu_{2}-\delta)/2}),$
(A.23)

and we can similarly show that

$\displaystyle T^{-\delta/2}\|\mathbf{J}_{8}\|=O_{P}(T^{(\nu_{2}-\delta)/2}).$
(A.24)

By putting everything together, (A.8) becomes

$\displaystyle
T^{-\delta/2}\|T^{-(\nu_{1}+\delta)/2}\widehat{\mathbf{F}}_{1}\mathbf{V}_{1}-\mathbf{J}_{9}\|$
$\displaystyle=O_{P}(T^{-(\delta-1)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|)$
$\displaystyle+O_{P}(N^{-1/2}T^{1-(\nu_{1}+\delta)/2})+O_{P}(T^{-(\delta-\nu_{2})/2}).$
(A.25)

We now left multiply (A.8) by
$T^{-(\nu_{1}+\delta)/2}\widehat{\mathbf{F}}_{1}^{\prime}$ to obtain that

$\displaystyle
T^{-\nu_{1}}\mathbf{V}_{1}-(T^{-(\nu_{1}+\delta)/2}\widehat{\mathbf{F}}_{1}^{\prime}\mathbf{F}_{1}^{0})(N^{-1}\boldsymbol{\Gamma}_{1}^{0\prime}\boldsymbol{\Gamma}_{1}^{0})(T^{-(\nu_{1}+\delta)/2}\mathbf{F}_{1}^{0\prime}\widehat{\mathbf{F}}_{1})$
$\displaystyle=T^{-(\nu_{1}-\delta)/2}O_{P}(T^{-(\delta-1)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|+N^{-1/2}T^{1-(\nu_{1}+\delta)/2}+T^{-(\delta-\nu_{2})/2})$
$\displaystyle=O_{P}(T^{-(\nu_{1}-1)/2}\|\mathbf{D}_{T}^{-1}(\boldsymbol{\beta}^{0}-\widehat{\boldsymbol{\beta}}_{0})\|)+O_{P}(T^{-(\nu_{1}-\nu_{2})/2})$
$\displaystyle=O_{P}(T^{-(\nu_{1}-\nu_{2})/2}),$ (A.26)

where the third equality follows from Assumption 4 and Lemma 3.1. This implies
that $\mathbf{V}_{1}$ is at most of rank $d_{1}$. Similarly, we can left-
multiply (A.8) by $T^{-\nu_{1}}\mathbf{F}_{1}^{0\prime}$ to obtain