EconBase
← Back to paper

Linear Regression for Panel With Unknown Number of Factors as Interactive Fixed Effects

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

86,743 characters

Linear Regression for Panel with Unknown Number of Factors as Interactive Fixed Effects



\title{\bf Linear Regression for Panel with Unknown Number of Factors as Interactive Fixed Effects\footnote{This is the last working paper version of the paper published in \textit{Econometrica} \textbf{83}(4), 1543--1579, 2015; doi:10.3982/ECTA9382. The copyright to this article is held by the Econometric Society. It may be downloaded, printed and reproduced only for personal or classroom use. We thank the participants of the 2009
Cowles Summer Conference
``Handling Dependence: Temporal, Cross-sectional, and Spatial'' at Yale University,
of the 2012 North American Summer Meeting of the Econometric Society at Northwestern University,
of the 18th International Conference on Panel Data at the Banque de France,
of the 2013 North American Winter Meeting of the Econometric Society in San Diego,
of the 2014 Asia Meeting of Econometric Society in Taipei,
of the 2014 Econometric Study Group Conference in Bristol, and of the econometrics seminars in USC
and Toulouse
for many interesting comments, and we thank Dukpa Kim, Tatsushi Oka, and Alexei Onatski for helpful discussions.
We are also grateful for the comments and suggestions of James Stock,  Elie Tamer, and anonymous referees.
Moon acknowledges financial supports of the NSF via SES 0920903 and the faculty grant award of USC.
Weidner acknowledges support from the Economic and Social Research Council through the ESRC Centre for Microdata Methods and Practice grant RES-589-28-0001.
}}

\author{\setcounter{footnote}{2}
Hyungsik Roger Moon\footnote{
Department of Economics and USC Dornsife INET, University of Southern California,
Los Angeles, CA 90089-0253.
Email: {\tt [email removed]}.
Department of Economics, Yonsei University, Seoul, Korea.
}
\and Martin Weidner\footnote{Corresponding author.
 Department of Economics,
 University College London,
 Gower Street,
 London WC1E~6BT, U.K.,
 and CeMMaP.
 Email: {\tt [email removed]}.
}}


\date{December 2014}

\maketitle


\abstract{
\begin{center}
\vspace{3mm}
\begin{minipage}{0.75\textwidth}
\footnotesize
  In this paper we study the least squares (LS) estimator in a linear panel
  regression model with \emph{unknown} number of factors appearing as interactive fixed effects.
   Assuming that the number of factors used in estimation is larger than the true number of factors in the data
  we establish the limiting distribution of the LS estimator for the regression coefficients, as
   the number of time periods and the number of cross-sectional units jointly go to infinity.
  The main result of the paper is that under certain assumptions the
  limiting distribution of the LS estimator is independent of the number of factors
  used in the estimation, as long as this number is not underestimated.
  The important practical implication
  of this result is that for inference on the regression
  coefficients one does not necessarily need to estimate the number of
  interactive fixed effects consistently.

\end{minipage}
\vspace{3mm}
\end{center}
}

\bigskip

\noindent{\bf Keywords:}
Panel data, interactive fixed effects, factor models,
perturbation theory of linear operators, random matrix theory.\\[4pt]
\noindent{\bf JEL-Classification:} C23, C33


\linespread{1.3}

\newpage

\section{Introduction}

Panel data models typically incorporate individual and time effects
to control for heterogeneity in cross-section and over time.
While often these individual and time effects enter
the model additively, they can also be interacted multiplicatively, thus giving
rise to so called interactive effects, which we also refer to as a
factor structure. The multiplicative form
captures the heterogeneity in the data more flexibly, since it allows for common
time-varying shocks (factors) to affect the cross-sectional units
with individual specific
sensitivities (factor loadings).\footnote{
The conventional additive model can be interpreted as a two factor interactive fixed effects model.}
It is this flexibility that motivated
the discussion of interactive effects in the econometrics literature, e.g.
Holtz-Eakin, Newey and Rosen~\cite*{HoltzEakin-Newey-Rosen1988},
Ahn, Lee and Schmidt \cite*{AhnLeeSchmidt2001,AhnLeeSchmidt2013},
Pesaran~\cite*{Pesaran2006}, Bai~\cite*{Bai2009,Bai2013likelihood},  Zaffaroni \cite*{Zaffaroni2009},
Moon and Weidner~\cite*{MoonWeidner2013}, and Lu and Su \cite*{LuSu2013}.

Let $N$ be the number of cross-sectional units, $T$ be the number of time periods,
$K$ be the number of regressors, and $R^0$ be the true number of interactive fixed effects.
We consider a linear regression model with observed outcomes $Y$, regressors $X_k$,
and unobserved error structure $\varepsilon$, namely
\begin{align}
   Y &= \sum_{k=1}^{K} \, \beta_{k}^{0} \, X_k \, + \, \varepsilon
        \, ,
   &
   \varepsilon &= \lambda^0 \, f^{0 \, \prime} + e  \, ,
   \label{model}
\end{align}
where $Y$, $X_k$, $\varepsilon$ and $e$ are $N\times T$ matrices,
$\lambda^0$ is an $N \times R^0$ matrix, $f^0$ is a $T \times R^0$ matrix,
and the regression parameters $\beta^0_k$ are scalars
 --- the superscript zero indicates the true value of the parameters.
We write $\beta$ for the $K$-vector of regression parameters,
and we denote the components of the different matrices
by $Y_{it}$, $X_{k,it}$, $e_{it}$, $\lambda^0_{ir}$ and $f^0_{tr}$,
where $i = 1, \ldots, N$, $t=1, \ldots, T$, and $r=1,\ldots,R^0$.
It is convenient to introduce the notation
$\beta \cdot X := \sum_{k=1}^{K} \, \beta_{k} \, X_k$.
All matrices, vectors and scalars in this paper are real valued.

We consider the interactive fixed effect specification,
i.e. we treat $\lambda^0$ and $f^0$ as nuisance parameters, which are estimated
jointly with the parameters of interest $\beta$.\footnote{
When we refer to
interactive fixed effects we mean that both factors and factor loadings
are treated as non-random parameters.
Ahn, Lee and Schmidt \cite*{AhnLeeSchmidt2001}
take a hybrid approach in that they treat the factors as non-random,
but the factor loadings as random.
The common correlated effects estimator of
Pesaran~\cite*{Pesaran2006} was introduced in a context, where both the factor loadings
and the factors follow certain probability laws, but it exhibits many  properties
of a fixed effects estimator.}
The advantages of the fixed effects approach are for instance that it is semi-parametric,
since no assumption on the distribution of the interactive effects needs to be made,
and that the regressors can be arbitrarily correlated with the interactive effect
parameters.

We study the least squares (LS) estimator of model \eqref{model}, which
 minimizes the sum of squared residuals
to estimate the unknown parameters $\beta$, $\lambda$ and $f$.\footnote{
The LS estimator
is sometimes called ``concentrated'' least squares estimator in the literature,
and in an earlier version of the paper we referred to it as the ``Gaussian Quasi Maximum Likelihood Estimator'',
since LS estimation is equivalent to maximizing a conditional Gaussian likelihood function.
Note also that for fixed $\beta$ the LS estimator for $\lambda$
and $f$ is simply the principal components estimator.
}
To our knowledge, this estimator was first discussed in Kiefer~\cite*{Kiefer1980}. Under an asymptotic where $N$ and $T$ grow to infinity, the asymptotic properties of the LS estimator were derived in Bai~\cite*{Bai2009} for strictly exogeneous regressors, and extended in Moon and Weidner~\cite*{MoonWeidner2013} to the case of pre-determined regressors.

An important restriction of these papers is that the number of factors $R^0$ is
assumed to be known.
However, in many empirical applications there is no consensus about the exact number of factors in the data or in the relevant economic model.
If $R^0$ is not known beforehand, then it may be estimated consistently,\footnote{See the discussion
in Bai~\cite*{Bai2009supp}, supplemental material, regarding estimation of $R^0$.}
but difficulties in obtaining reliable estimates for the number of factors
are well-documented in the literature (see, e.g., the simulation results in Onatski~\cite*{Onatski2010}, and also our empirical illustration in Section~\ref{sec:Empirical}). Furthermore, in order to use the existing inference results on  $R^0$ one still needs a good preliminary estimator for $\beta$, so that
working out the asymptotic properties of the LS-estimator for $R \geq R^0$ is still useful when taking that route.

We investigate the asymptotic properties of the LS estimator when the true number of factors $R^0$ is unknown and $R$ ($\geq R^0$) number of factors are used in the estimation.\footnote{
For $R<R^0$ the LS estimator can be inconsistent, since
then there are interactive fixed effects in  the model which can be correlated
with the regressors but are not controlled for in the estimation.
We therefore restrict attention to the case $R \geq R^0$.
}
We denote this estimator by $\widehat{\beta}_{R}$.

The main result of the paper, presented in Section~\ref{sec:main}, is that under certain assumptions
the LS estimator  $\widehat{\beta}_{R}$ has the same limiting distribution as $\widehat{\beta}_{R^0}$
for any $R \geq R^0$ under an asymptotic where both $N$ and $T$ become large, while  $R^0$  and $R$ are constant.
This implies that the LS estimator $\widehat{\beta}_{R}$ is asymptotically robust towards inclusion
of  extra interactive effects in the model, and within the LS estimation framework
there is no asymptotic efficiency loss from choosing $R$ larger than $R^0$.
The important empirical implication of our result is that the number of factors $R^0$ need not be known or estimated accurately to apply the LS estimator.

To derive this robustness result, we impose more restrictive conditions than those typically assumed with known $R^0$. These include that the errors $e_{it}$ are independent and identically (iid) normally distributed and that the regressors are composed of a ``low-rank''
strictly stationary component, a ``high-rank'' strictly stationary component, and a ``high-rank'' pre-determined
component.\footnote{
The pre-determined component of the regressors allows for linear feedback of $e_{it}$ into future realizations of $X_{k,it}$.
}
Notice that while some of these restrictions are necessary for our robustness result,
some of them (e.g.~iid normality of $e_{it}$) are imposed for technical reasons, because
in the proof we use certain results from the
theory of random matrices that are currently only available in that case (see the discussion in Section~\ref{sec:AsymptoticSummary}).
In the Monte Carlo simulations in Section~\ref{sec:MC}, we consider DGPs that violate some technical conditions to demonstrate robustness of the result.

Under less restrictive assumptions we provide intermediate results that sequentially lead to the main result in Section~\ref{sec:AsyTheory} and Appendix~\ref{sec:ConvergenceRate} and~\ref{sec:Equivalence}.
In Section~\ref{sec:consistency} we show $\sqrt{\min(N,T)}$-consistency
of the LS estimator $\widehat{\beta}_{R}$ as $N,T \rightarrow \infty$
under very mild regularity condition on $X_{it}$ and $e_{it}$, and without imposing any assumptions
on $\lambda^0$ and $f^0$ apart from $R \geq R^0$. We thus obtain consistency of the LS estimator
not only for unknown number of factors, but also for weak factors,\footnote{
See  Onatski~\cite*{Onatski2010,Onatski2012} and
Chudik, Pesaran and Tosetti~\cite*{ChudikPesaranTosetti2011}
for a discussion of ``strong'' vs. ``weak'' factors in factor models.} which is an important robustness result.

In Section~\ref{sec:expansion} we derive an asymptotic expansion of the LS profile objective function that concentrates out $f$ and $\lambda$,
for the case $R=R^0$. Given that the profile objective function is a sum of eigenvalues of a covariance matrix, its quadratic approximation is challenging because the derivatives of the eigenvalues with respect to $\beta$ are not generally known. We thus cannot use a conventional Taylor expansion,
but instead apply the perturbation theory of linear operators to derive
the approximation.

In Section~\ref{sec:AsymptoticSummary} we provide an example that satisfies the typical assumptions imposed with known $R^0$, so that $\widehat{\beta}_{R^0}$ is $\sqrt{NT}$ consistent,
but we show that $\widehat{\beta}_{R}$ with $R > R^0$ is only $\sqrt{\min(N,T)}$ consistent in that example. This shows that stronger conditions
are required to derive our main result.

In Appendix~\ref{sec:ConvergenceRate}
we show faster than $\sqrt{\min(N,T)}$-convergence of  $\widehat{\beta}_{R}$ under
assumptions that are less restrictive than those employed for the main result, in particular allowing
for either cross-sectional or time-serial correlation of the errors $e_{it}$. In Appendix~\ref{sec:Equivalence} we provide
an alternative version of our main result of asymptotic equivalence of  $\widehat{\beta}_{R^0}$
and  $\widehat{\beta}_{R}$, $R \geq R^0$, which is derived under high-level assumptions.


In Section~\ref{sec:Empirical} we follow Kim and Oka~\cite*{KimOka2014} in employing the interactive fixed effects
specification to study the effect of US divorce law reforms on divorce rates. This empirical example illustrates that
the estimates for the coefficient $\beta$ indeed become insensitive to the choice of $R$, once $R$ is chosen sufficiently
large, as expected from our theoretical results.

Section~\ref{sec:MC} contains Monte Carlo simulation results for a static panel model.
For the simulations we consider a DGP that violates the iid normality restriction of the error term. The simulation results confirm our main result of the paper even with a relatively small sample size (e.g. $N=100$, $T=10$) and non-iid-normal errors. In the supplementary appendix, we report the Monte Carlo simulation results of an AR(1) panel model. It also confirms the robustness result in large samples, but in finite samples it shows more inefficiency than the static case. In general, one should expect some finite sample inefficiency from overestimating the number of factors when the sample size is small or the number of overfitted factors is large.

A few words on notation.
The transpose of a matrix $A$ is denoted by $A'$.
For a column vectors $v$ its Euclidean norm is
defined by $\| v \| = \sqrt{v^{\prime}v}$ .
For an $m\times n$ matrix $A$
the Frobenius or Hilbert Schmidt norm is $\| A \|_{HS} = \sqrt{{\rm Tr}
(AA^{\prime})}$, and the operator or spectral norm is $\| A \| = \max_{0 \neq v \in
\mathbb{R}^n} \, \frac{ \| A v \|} {\| v\|}$.
Furthermore, we use $P_A = A
(A^{\prime}A)^{\dagger} A'$ and $M_A = \mathbbm{1} - A (A^{\prime}A)^{\dagger} A'$,
where $\mathbbm{1}$ is the $m\times m$ identity matrix,
and $(A^{\prime}A)^{\dagger}$ denotes some generalized inverse, in case
$A$ is not of full column rank. For
square matrices $B$, $C$, we use $B>C$ (or $B\geq C$) to indicate that $B-C$ is positive (semi) definite. We use ``wpa1'' for ``with probability approaching one''.








\section{Identification of $\beta^0, \lambda^0 f^{0 \prime}$, and $R^0$}
\label{sec:model_id}

In this section we provide a set of conditions under which the regression coefficient $\beta^0$, the interactive fixed effects $\lambda^0 f^{0 \prime}$, and the number of factors $R^0$ are determined uniquely by the data.
Here, and throughout the whole paper, we treat $\lambda$ and $f$ as non-random
parameters, i.e. all stochastics in the following
are implicitly conditional on $\lambda$ and $f$.
Let $x_k={\rm vec}(X_k)$, the $NT$-vectorization of $X_k$,
and let $x=(x_1,\ldots,x_K)$, which is an
$NT \times K$ matrix.

\begin{IDassumption}[\bf Assumptions for Identification]
   There exists a non-negative integer $R$ such that
   \label{ass:id}
   \begin{itemize}
      \item[(i)] The second moments of $X_{it}$ and
                 $e_{it}$ exist for all $i$, $t$.

      \item[(ii)]  $\mathbbm{E}(e_{it})=0$,
     $\mathbbm{E}(X_{it} e_{it})=0$,
   for all $i$, $t$.

      \item[(iii)] $\mathbbm{E}[x' (M_F \otimes M_{\lambda^0}) x]
     >0 $, for all $F \in \mathbbm{R}^{T \times R}$.

     \item[(iv)]  $R \geq R^0 := {\rm rank}(\lambda^0 f^{0 \prime})$.

   \end{itemize}
\end{IDassumption}

\begin{theorem}[\bf Identification]
   \label{th:id}
     Suppose that  the Assumptions~\ref{ass:id} are satisfied. Then, $\beta^0$, $\lambda^0 f^{0 \prime}$, and $R^0$ are identified.\footnote{
Here, identification means that
$\beta^0$ and $\lambda^0 f^{0 \prime}$
 can be uniquely recovered from the distribution of
$(Y,X)$ conditional on those parameters.
Identification of the number of factors follows
since $R^0 = {\rm rank}( \lambda^0 f^{0 \prime} )$.
    The factor loadings
and factors $\lambda^0$ and $f^{0}$ are not separately
identified without further normalization restrictions,
but the product $\lambda^0 f^{0 \prime}$ is identified.
}
\end{theorem}

Assumption~\ref{ass:id}$(i)$ imposes existence of second moments.
Assumption~\ref{ass:id}$(ii)$ is an exogeneity condition,
which demands that $x_{it}$ and $e_{it}$ are not correlated contemporaneously,
but allows for pre-determined regressors like lagged dependent variables.
Assumption~\ref{ass:id}$(iv)$ imposes that the true number of factors
$R^0 := {\rm rank}(\lambda^0 f^{0 \prime})$ is bounded by a non-negative integer
$R$, which cannot be too large (e.g. the trivial bound $R=\min(N,T)$ is not possible),
since otherwise Assumption~\ref{ass:id}$(iii)$
cannot be satisfied.

Assumption~\ref{ass:id}$(iii)$ is a non-collinearity condition,
which demands that the regressors have significant variation across $i$ and over $t$ after projecting out all variation that
can be explained by the factor loadings $\lambda^0$ and by arbitrary
factors $F \in \mathbbm{R}^{T \times R}$.
This generalizes the within variation assumption in the conventional panel regression with time invariant individual fixed effects, which in our notation reads
$\mathbbm{E}[x' (M_{1_T} \otimes \mathbbm{1}_N) x] >0$.\footnote{
The conventional panel regression with additive individual fixed effects and time effects requires
a non-collinearity condition of the form $\mathbbm{E}[x' (M_{1_T} \otimes M_{1_N}) x] >0$.}
This conventional fixed effect assumption rules out time-invariant regressors.
Similarly, Assumption~\ref{ass:id}$(iii)$ rules out more general ``low-rank regressors'',\footnote{
We do not consider such ``low-rank regressors'' in this paper.
Note also that Assumption A in Bai~\cite*{Bai2009} is the sample version of our Assumption~\ref{ass:id}$(iii)$.}
see our discussion of Assumption~\ref{ass:NC} below.

\section{Main Result}
\label{sec:main}

The estimator we investigate in this paper is the least squares (LS) estimator,
which for a given choice of $R$ reads\footnote{
The optimal $\widehat \Lambda_R$ and $\widehat F_R$ in \eqref{estimator} are not unique,
since the objective function is invariant under right-multiplication of $\Lambda$
with a non-degenerate $R \times R$ matrix $S$, and simultaneous right-multiplication
of $F$ with $(S^{-1})'$. However, the column spaces of $\widehat \Lambda_R$ and $\widehat F_R$
are uniquely determined.
}
\begin{align}
   \left( \widehat \beta_{R},\,\widehat \Lambda_R,\,\widehat F_R \right) &\in
     \;
  \operatorname*{argmin}_{ \left\{ \beta \in \mathbbm{R}^K, \; \Lambda \in \mathbbm{R}^{N\times R}, \;
                F \in \mathbbm{R}^{T\times R} \right\}  } \;
    \left\| Y \, - \, \beta \cdot X \, - \, \Lambda \, F'
    \right\|^2_{HS}    \; ,
   \label{estimator}
\end{align}
where $\|.\|_{HS}$ refers to the Hilbert Schmidt norm, also called Frobenius norm.
The objective function $ \left\| Y \, - \, \beta \cdot X \, - \, \Lambda \, F'
    \right\|^2_{HS}$ is simply the sum of squared residuals.
The estimator for $\beta^0$ can equivalently be defined by minimizing
the profile objective function that concentrates out the $R$ factors and the $R$ factor loadings, namely
\begin{align}
   \widehat \beta_R &= \operatorname*{argmin}_{\beta \in \mathbbm{R}^K} \, {\cal L}_{NT}^R(\beta) \; ,
   \label{DefLS estimator}
\end{align}
with\footnote{
The profile objective function ${\cal L}_{NT}^R(\beta)$ need not be convex in $\beta$
and can have multiple local minima.
Depending on the dimension of $\beta$
one should either perform an initial grid search or try multiple starting values
for the optimization when calculating the global minimum $\widehat \beta_R$ numerically. See also
Section~\ref{app:Numerics} of the supplementary material.
 }
\begin{align}
   {\cal L}_{NT}^R(\beta) &= \min_{\left \{ \Lambda \in \mathbbm{R}^{N\times R}, \;
                F \in \mathbbm{R}^{T\times R} \right\}  } \;
                \frac 1 {NT}  \;
    \left\| Y \, - \, \beta \cdot X \, - \, \Lambda \, F'
    \right\|^2_{HS}
  \nonumber \\
      &= \min_{F \in \mathbbm{R}^{T\times R}} \; \frac 1 {NT}  \;
          {\rm Tr}\left[ \left(Y - \beta \cdot X\right)
                         M_F
                        \left(Y - \beta \cdot X\right)' \right]
 \nonumber \\
      &= \; \frac 1 {NT}  \; \sum_{r=R+1}^{T}
            \mu_r\left[ \left(Y - \beta \cdot X \right)'
                     \left(Y - \beta \cdot X \right) \right] \; ,
   \label{LSobjective}
\end{align}
where, $\mu_r(.)$ is the $r$'th largest eigenvalue of the matrix argument.
Here, we first concentrated out $\Lambda$ by use of its own first order condition.
The resulting optimization problem for $F$ is a principal components
problem, so that the the optimal $F$ is given by the $R$ largest principal components
of the $T \times T$ matrix $\left(Y - \beta \cdot X\right)'\left(Y - \beta \cdot X\right)$. At the optimum the projector $M_F$ therefore exactly projects out
the $R$ largest eigenvalues
of this matrix, which gives rise to the final formulation of the
profile objective function as the sum over its $T-R$ smallest eigenvalues.\footnote{
This last formulation of ${\cal L}_{NT}^R(\beta)$ is very convenient
since it does not involve any explicit optimization over nuisance parameters.
Numerical calculation of eigenvalues is very fast,
so that the numerical evaluation of ${\cal L}_{NT}^R(\beta)$ is unproblematic
for moderately large values of $T$.
Since the model is symmetric under $N \leftrightarrow T$,
$\Lambda \leftrightarrow F$, $Y \leftrightarrow Y'$, $X_k \leftrightarrow X_k'$ there also exists
a dual formulation of ${\cal L}_{NT}^R(\beta)$ that involves solving an eigenvalue
problem for an $N \times N$ matrix.}
We write ${\cal L}_{NT}^0(\beta)$ for ${\cal L}_{NT}^{R^0}(\beta)$, the profile objective function obtained for the
true number of factors.
Notice that in \eqref{DefLS estimator} the parameter set for $\beta$ is the whole Euclidean space $\mathbbm{R}^K$ and we do not restrict the parameter set to be compact.

\begin{SFassumption}[\bf Strong Factor Assumption]~
\label{ass:SF}
\begin{itemize}
   \item[(i)] $0 < \operatorname*{plim}_{N,T \rightarrow \infty}
                               \frac 1 N \, \lambda^{0\prime} \lambda^0
                             < \infty$. \quad
   \item[(ii)] $0 <  \operatorname*{plim}_{N,T \rightarrow \infty}
                               \frac 1 T \, f^{0\prime} f^0  < \infty$.
\end{itemize}
\end{SFassumption}

\begin{NCassumption}[\bf Non-Collinearity of $X_k$]~
  \label{ass:NC}
  Consider linear combinations
  $\alpha \cdot X :=\sum_{k=1}^K \alpha_k X_k$ of the regressors $X_k$ with
  $K$-vector $\alpha$ such that $\|\alpha\|=1$.
  We assume that there exists a constant $b>0$ such that
  \begin{align*}
     \min_{\{\alpha \in \mathbb{R}^{K}, \, \|\alpha\|=1\}}   \,
     \sum_{r=R+R^0+1}^{T} \,
      \mu_{r} \left[ \frac{ (\alpha \cdot X)' (\alpha \cdot X)} {NT}\right] \;  &\geq b \; , \qquad
     \text{wpa1.}
  \end{align*}
\end{NCassumption}


\begin{LLassumption}[\bf Low Level Conditions for Main Result]  $\phantom{a}$
\label{ass:LL}
\begin{itemize}
      \item[(i)] {\bf Decomposition of Regressors:}
      $X_k = \overline X_k  + \widetilde X_k^{\rm str} +  \widetilde X_k^{\rm weak}$,
      for $k=1,\ldots,K$, where
     $\overline X_k$, $\widetilde X_k^{\rm str}$ and $\widetilde X_k^{weak}$ are $N \times T$ matrices, and
     \begin{itemize}

     \item[(i.a)] {\bf Low-Rank (strictly exogenous) Part of Regressors:} ${\rm rank}(\overline X_k)$ is bounded as $N,T \rightarrow \infty$,
     and  $\frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T \overline X_{k,it}^2  = {\cal O}_P(1)$.

    \item[(i.b)] {\bf High-Rank (strictly exogenous) Part of Regressors:}
       $\| \widetilde X_k^{\rm str} \| = {\cal O}_P(N^{3/4})$, as can be justified e.g. by Lemma~\ref{lemma:SpectralNormBound}
       in the appendix.

     \item[(i.c)] {\bf Weakly Exogenous Part of Regressors:}
     $\widetilde X_{k,it}^{\rm weak} = \sum_{\tau=1}^{t-1} \gamma_\tau e_{i,t-\tau}$,
     where the real valued coefficients $\gamma_\tau$ satisfy
     $\sum_{\tau=1}^\infty | \gamma_\tau | < \infty$.

    \item[(i.d)]  {\bf Bounded Moments:}            We assume that
                   $\mathbbm{E} \left| X_{k,it}    \right|^{2}$,
                  $\mathbbm{E} \left|(M_{\lambda^0} X_k M_{f^0})_{it}
                               \right|^{26}$,
                 $\mathbbm{E} \left|(M_{\lambda^0} X_k)_{it}
                               \right|^{8}$
                 and
                 $\mathbbm{E} \left|(X_k M_{f^0})_{it}
                               \right|^{8}$
                  are bounded uniformly over $k$, $i$, $j$, $N$ and $T$.

    \end{itemize}

      \item[(ii)] {\bf Errors are iid Normal:}
      The error matrix $e$ is independent of
                 $\lambda^0$, $f^0$,
                 $\overline X_k$, and $\widetilde X_k^{\rm str}$, $k=1,\ldots,K$,
                 and its elements $e_{it}$ are
                 independent and identically distributed as ${\cal N}(0,\sigma^2)$ across $i$ and over $t$.

     \item[(iii)] {\bf Number of Factors not Underestimated:}
     $R \geq R^0 := {\rm rank}(\lambda^0 f^{0 \prime})$.
   \end{itemize}
\end{LLassumption}






\paragraph{Remarks}


\begin{itemize}
  \item[(i)] Assumption~\ref{ass:SF} imposes that the factor $f^0$ and the factor loading $\lambda^0$ are strong. The strong factor assumption is regularly imposed in the literature on large $N$ and $T$ factor models, including Bai and Ng~\cite*{BaiNg2002},
Stock and Watson~\cite*{StockWatson2002} and Bai~\cite*{Bai2009}.

  \item[(ii)] Assumption~\ref{ass:NC} demands that there exists significant sampling variation in the regressors after concentrating out $R+R^0$ factors (or factor
  loadings). It is a sample version of the identification Assumption~\ref{ass:id}$(iii)$, and it is essentially equivalent to Assumption~A of Bai~\cite*{Bai2009},
  but avoids mentioning the unobserved loadings $\lambda^0$.\footnote{By dropping the expected value from
  Assumption~\ref{ass:id}$(iii)$ and replacing the zero lower bound by
  a positive constant
   one obtains
  $\inf_{F} \left[ x' (M_F \otimes M_{\lambda^0}) x /NT \right] \geq b >0$, wpa1,
  which is equivalent to Assumption A of Bai~\cite*{Bai2009},
  and can also be rewritten as
 $\min_{\| \alpha \|=1} \inf_{F}
 {\rm Tr}\left[ M_{\lambda^0} (\alpha \cdot X)' M_{F} (\alpha \cdot X)/NT \right] \geq b$.
  A slightly stronger version of the Assumption, which avoids mentioning
  the unobserved factor loading $\lambda^0$,
  reads
  $\min_{\| \alpha \|=1} \inf_{F} \inf_{\lambda}
 {\rm Tr}\left[ M_{\lambda} (\alpha \cdot X)' M_{F} (\alpha \cdot X)/NT \right] \geq b$,
where $F \in \mathbbm{R}^{T \times R}$ and $\lambda \in \mathbbm{R}^{N \times R^0}$,
and this slightly stronger version is equivalent to Assumption~\ref{ass:NC}.
}

  \item[(iii)] Assumption~\ref{ass:NC} is violated if there exists a linear combination $\alpha \cdot X$ of the regressors
with $\alpha \neq 0$ and ${\rm rank}(\alpha \cdot X) \leq R+R^0$, i.e. the assumption
rules out ``low-rank regressors'' like time invariant regressors
or cross-sectionally invariant regressors. These low-rank regressors
require a special treatment in the interactive fixed effect model,
see Bai~\cite*{Bai2009} and Moon and Weidner~\cite*{MoonWeidner2013},
and we do not consider them in the present paper.
If one is not interested explicitly
in their regression coefficients, then one can always
eliminate the low-rank regressors by an appropriate projection of the data, e.g.
subtraction of the time (or cross-sectional) means from the data eliminates all
time-invariant (or cross-sectionally invariant) regressors,
see Section~\ref{sec:Empirical} for an example of this.

 \item[(iv)] The norm restriction in Assumption~\ref{ass:LL}$(i.b)$ is a high level assumption.
 It is satisfied as long as $\widetilde X_{k,it}^{\rm str}$ is mean zero and
  weakly correlated across $i$ and over $t$,  for details see Appendix~\ref{app:SpectralNorm} and
  Lemma~\ref{lemma:SpectralNormBound} there.

  \item[(v)] Assumption~\ref{ass:LL}$(i)$ imposes that each regressor consists of three parts: (a) a strictly exogenous low rank component , (b) a strictly exogenous component satisfying a norm restriction, and (c) a weakly exogenous component
that follows a linear process with innovation given by the lagged error term $e_{it}$.
For example, if $X_{k,it} \, \sim \, iid \, {\cal N}( \mu_k, \sigma_k^2)$, independent of $e$, then we have
$\overline X_{k,it}=\mu_k$, $\widetilde X_{k,it}^{\rm str} \, \sim \, iid \, {\cal N}( 0, \sigma_k^2)$ and
$\widetilde X_k^{\rm weak} = 0$.
Assumption~\ref{ass:LL}$(i)$ is also satisfied for a stationary panel VAR with interactive fixed effects as in Holtz-Eakin, Newey and Rosen~\cite*{HoltzEakin-Newey-Rosen1988}. A special case of this is a dynamic panel regression with fixed effects, where $Y_{it} = \beta Y_{i,t-1} + \lambda^{0 \prime}_i f^0_t + e_{it}$, with $| \beta |<1$ and ``infinite history''. In this case, we have
$X_{it} = Y_{i,t-1} = \overline X_{it}  + \widetilde X_{it}^{\rm str} +  \widetilde X_{it}^{\rm weak}$, where
$\overline X_{it} =  \lambda^{0 \prime}_i  \sum_{\tau=1}^\infty \beta^{\tau-1} f^0_{t-\tau}$,
 $\widetilde X_{it}^{\rm str} = \sum_{\tau=t}^\infty \beta^{\tau-1} e_{i,t-\tau}$,
 and $\widetilde X_{it}^{\rm weak} = \sum_{\tau=0}^{t-1}  \beta^{\tau-1} e_{i,t-\tau}$.

  \item[(vi)] Assumption~\ref{ass:LL}$(i)$ is more restrictive than Assumption 5 in Moon and Weidner~\cite*{MoonWeidner2013}, where $R^0$ is assumed to be known. However, it is more general than the restriction on the regressors in Pesaran~\cite*{Pesaran2006}, where  -- in our notation -- the decomposition
  $X_k = \overline X_k  + \widetilde X_k^{\rm str}$ is  imposed,
  but the lower rank component $\overline X_k$ needs to satisfy further
  assumptions,
   and the weakly exogenous component   $\widetilde X_k^{\rm weak} $ is not considered.
Bai~\cite*{Bai2009} requires no such decomposition, but
imposes strict exogeneity of the regressors.

    \item[(vii)]    Among the conditions in Assumption~\ref{ass:LL},
    the iid normality condition in Assumption~\ref{ass:LL}$(ii)$ may be the most restrictive.
    In Appendix~\ref{sec:Equivalence} we provide an alternative version of Theorem~\ref{th:MAIN}
    that imposes more general high-level conditions.
    Verifying those high-level conditions requires results on the eigenvalues and eigenvectors
    of random covariance matrices, which can be verified for iid normal errors by using
    known results from the random matrix theory literature,
    see Section~\ref{sec:AsymptoticSummary} for more details.
    We believe, however, that those high-level conditions and thus our main result hold more generally,
    and we explore non-normal and serially correlated errors in our Monte Carlo simulations below.
\end{itemize}



\begin{theorem}[\bf Main Result]
\label{th:MAIN}
Let Assumption~\ref{ass:SF}, \ref{ass:NC}
and \ref{ass:LL} hold and
consider a limit $N,T \rightarrow \infty$ with $N/T \rightarrow \kappa^2$, $0<\kappa<\infty$. Then we have
\begin{align*}
    \sqrt{NT}\big(\widehat \beta_{R} - \beta^0\big)
    =\sqrt{NT}\big(\widehat \beta_{R^0} - \beta^0\big) + o_P(1).
\end{align*}
\end{theorem}

Theorem~\ref{th:MAIN} follows from Theorem~\ref{th:LimitingDistribution} and  Lemma~\ref{lemma:JustifyEV} in the appendix, whose proof is given in the supplementary material.
The theorem guarantees that the asymptotic distribution of $\widehat \beta_{R}$, $R \geq R^0$,
is identical to that of $\widehat \beta_{R^0}$ in \eqref{AsyDistribution} below.

The limiting distribution of $\sqrt{NT}\big(\widehat \beta_{R^0} - \beta^0\big)$ with known $R^0$ is available in the existing literature. According to Bai~\cite*{Bai2009} and Moon and Weidner~\cite*{MoonWeidner2013},
\begin{align}
    \sqrt{NT}\big(\widehat \beta_{R^0} - \beta^0\big) \;  \Rightarrow \;
    {\cal N}\left( -  \kappa \, \operatorname*{plim}  W^{-1} B ,  \;
    \sigma^2  \operatorname*{plim}  W^{-1} \right) ,
    \label{AsyDistribution}
\end{align}
where $W$ is the $K \times K$ matrix with elements
   $W_{k_1 k_2} =  \frac 1 {NT} {\rm Tr}( M_{\lambda^0} X_{k_1} M_{f^0} X_{k_2}' )$ and
 $B$ is the $K$-vector with elements
$B_k = \frac 1 N {\rm Tr}[ P_{f^0} \mathbbm{E}(e'  X_k)]$.\footnote{
The asymptotic distribution in \eqref{AsyDistribution} can also be derived from Corollary~\ref{cor:LimitR0} below
under more general conditions than in Assumption~\ref{ass:LL} (see
Moon and Weidner~\cite*{MoonWeidner2013} for details).
Here we have used the homoscedasticity of $e_{it}$ to simplify the structure of the asymptotic variance
and bias. Bai~\cite*{Bai2009} finds further asymptotic bias
in $\widehat \beta_{R^0}$ due to heteroscedasticity and correlation in $e_{it}$,
which in our asymptotic result is ruled out by Assumption~\ref{ass:LL}$(ii)$,
but is studied in our Monte Carlo simulations below.
Moon and Weidner~\cite*{MoonWeidner2013} work out the additional asymptotic
bias in $\widehat \beta_{R^0}$ due to pre-determined regressors, which
is allowed for in Theorem~\ref{th:MAIN}.
}

The result \eqref{AsyDistribution} holds under the assumptions of Theorem~\ref{th:MAIN}
and also assuming that $ \operatorname*{plim}  W^{-1} B$ and   $\operatorname*{plim}  W^{-1}$ exist, where
$\operatorname*{plim}$ refers to the probability limit as $N,T \rightarrow \infty$.
Note that Assumption~\ref{ass:NC} guarantees that $W$ is invertible asymptotically.
The asymptotic bias in \eqref{AsyDistribution} is an incidental parameter bias due to pre-determined regressors and is equal to zero for
strictly exogenous regressors (for which $\mathbbm{E}(e'  X_k)=0$); it generalizes the well-known
Nickell~\cite*{Nickell1981} bias of the within-group estimator for dynamic panel models.

Estimators for $\sigma^2$, $W$ and $B$ are given by\footnote{The first factor in $ \widehat \sigma^2$ reflects the degree of freedom correction from
estimating $\Lambda$, $F$ and $\beta$, but could simply be chosen as $1/NT$ for the purpose of consistency.
Note also that  $P_{\widehat F_R,t \tau} = {\cal O}_P(1/T)$, which explains why no $1/T$  factor is required in
the definition of $\widehat B_{R,k}$.}
\begin{align*}
     \widehat \sigma_R^2 &= \frac 1 {(N-R)(T-R)-K} \sum_{i=1}^N \sum_{t=1}^T
       \left( \widehat e_{R,it}  \right)^2 ,
    &
     \widehat W_{R,k_1 k_2} &= \frac 1 {NT}
      {\rm Tr}\left( M_{\widehat \Lambda_R} X_{k_1} M_{\widehat F_R} X_{k_2}'  \right) ,
    \\
     \widehat B_{R,k} &=
     \sum_{t=1}^T \sum_{\tau=t+1}^{t+M} P_{\widehat F_R,t \tau}  \left[ \frac 1 N  \sum_{i=1}^N  \widehat e_{R,it}  X_{k,i \tau} \right] ,
\end{align*}
where $\widehat e_{R,it}$ denotes the $(i,t)^{th}$ element of
$\widehat e_R = Y - \widehat \beta_R \cdot X- \widehat \Lambda_R \widehat F_R'$,
and $P_{\widehat F_R,t \tau}$ denotes the $(t,\tau)^{th}$ element of
$P_{\widehat F_R} = \mathbbm{1}_T - M_{\widehat F_R} =
\widehat F_R (\widehat F_R' \widehat F_R)^\dagger \widehat F_R' $,
and $M \in \{1,2,3,\ldots\}$ is a bandwidth parameter that also depends on the sample size $N,T$.
Let $\widehat W_R$ and $\widehat B_R$ be the matrix and vector with elements
$ \widehat W_{R,k_1 k_2}$ and $\widehat B_{R,k}$, respectively.

The next theorem establishes the consistency of these estimators.
Let $\lambda^{\rm red} \in \mathbbm{R}^{N \times (R-R^0)}$ and $f^{\rm red}  \in \mathbbm{R}^{T \times (R-R^0)}$ be the leading $R-R^0$ principal components obtained from the $N \times T$ matrix
$M_{\lambda^0} e M_{f^0}$, i.e. $\lambda^{\rm red}$ and $f^{\rm red}$ minimize the
objective function $\left\| M_{\lambda^0} e M_{f^0} -  \lambda^{\rm red} \, f^{\rm red \prime}  \right\|^2_{HS}$,
analogous to $\widehat \Lambda_R$ and $\widehat F_R$ defined in \eqref{estimator}.\footnote{The superscript
``red'' stands for redundant, because it turns out that $\lambda^{\rm red}$ and $f^{\rm red}$ are
asymptotically close to the $R-R^0$ redundant principal components that are estimated in~\eqref{estimator}.
}
\begin{theorem}[\bf Consistency of Bias and Variance Estimators]~
\label{th:Estimators}
\begin{itemize}
\item[(i)]
Let the conditions of Theorem~\ref{th:MAIN} hold.
Then we have $\left\| P_{\widehat F_R} - P_{[f^0,f^{\rm red}]} \right\| = o_p(1)$,
$\left\| P_{\widehat \Lambda_R} - P_{[\lambda^0,\lambda^{\rm red}]} \right\| = o_p(1)$,
$\widehat \sigma_R^2 = \sigma^2 + o_P(1)$,
and
$\widehat  W_R =  W + o_P(1)$.

\item[(ii)]
In addition, let $X_{k,\cdot t} = (X_{k,1t},...,X_{k,Nt})^{\prime}$, and assume that (1) $\gamma_\tau$ in Assumption~\ref{ass:LL}(i.c) satisfies $|\gamma_\tau| < c \tau^{-d}$
for some $c>0$ and $d>1$,
(2) $\| \lambda^0_i \|$ and $\| f^0_t \|$
are uniformly bounded over $i,t$ and $N,T$,
(3) $\max_t \| X_{k,\cdot t} \| = {\cal O}_P(\sqrt{N} \log N)$,\footnote{
The high-level assumption $\max_t \| X_{k,\cdot t} \| = {\cal O}_P(\sqrt{N} \log N)$
can be shown to be satisfied for the regressor component $\widetilde X^{\rm weak}_{k,it}$ above, and can be justified for
the other regressor components e.g. by assuming that
$\overline X_{k}$ and $\widetilde X_k^{\rm str}$ are uniformly bounded.
}
and (4) the bandwidth $M \rightarrow \infty$ such that
$M  (\log T)^2 T^{-1/6} \rightarrow 0$. Then, we have $\widehat B_R =  B + o_P(1)$.

\end{itemize}

\end{theorem}

Combining Theorems~\ref{th:MAIN} and~\ref{th:Estimators} and the asymptotic distribution in \eqref{AsyDistribution}
allows inference on $\beta$, for $R \geq R^0$. In particular, the bias corrected estimator
$\widehat \beta^{\rm BC}_{R} = \widehat \beta_{R} + \frac 1 T \widehat W_R^{-1} \widehat B_R$
satisfies\footnote
{Instead of estimating the bias analytically one can use the result that the bias
is of order $T^{-1}$ and perform split panel bias correction as in Dhaene and Jochmans~\cite*{DhaeneJochmans2010},
which instead of the conditions of Theorem~\ref{th:Estimators}(ii) only
requires some stationary condition over time.}
\[
\sqrt{NT}\big(\widehat \beta^{\rm BC}_{R} - \beta^0\big)  \Rightarrow
    {\cal N}( 0 , \sigma^2 W^{-1} ).
\]

\paragraph{Heuristic Discussion of the Main Result}~
\\
Intuitively, the inclusion of
unnecessary factors in the LS estimation is similar to the
inclusion of irrelevant regressors in an OLS regression.
In the OLS case it is well known that if those irrelevant extra regressors are uncorrelated with the regressors
of interest, then they have no effect on the asymptotic distribution of the regression coefficients of interest.
It is therefore natural to expect that if the extra estimated factors in $\widehat F_R$ are asymptotically uncorrelated
with the regressors, then the result of Theorem~\ref{th:MAIN} should hold.
To explore this,
remember that  $\widehat F_R$ is given by the
first $R$ principal components of the matrix
$(Y-\widehat \beta_R \cdot X)' (Y-\widehat \beta_R \cdot X)$, and write
\begin{align*}
  Y-\widehat \beta_R \cdot X
   &= \lambda^0 f^{0 \prime} + e - (\widehat \beta_R - \beta^0) \cdot X.
\end{align*}
The strong factor assumption and the consistency of $\widehat \beta_R$ guarantee
 that the first
$R^0$ principal components of $(Y-\widehat \beta_R \cdot X)' (Y-\widehat \beta_R \cdot X)$
are close to $f^0$ asymptotically, i.e. the true factors are correctly picked up by
the principal component estimator.
 The additional $R-R^0$ principal  components that are
estimated for $R>R^0$ cannot pick up anymore true factors and are thus
mostly determined by the remaining term
$e - (\widehat \beta_R - \beta^0) \cdot X$. The key question for the properties of
the extra estimated factors, and thus of $\widehat \beta_R$, is therefore
whether the principal components obtained from $e - (\widehat \beta_R - \beta^0) \cdot X$
are dominated by $e$ or by $(\widehat \beta_R - \beta^0) \cdot X$.
Only if they are dominated by $e$ can we expect the extra factors in $\widehat F_R$ to be uncorrelated with $X$
and thus the result in Theorem~\ref{th:MAIN} to hold.  The result on $P_{\widehat F_R}$ in
Theorem~\ref{th:Estimators} shows that the additional estimated factors are indeed close to $f^{\rm red}$,
i.e. are mostly determined by $e$, but this result is far from obvious a priori, as the following discussion shows.

Under our assumptions we have $\| e\| = {\cal O}_P(\sqrt{N})$ and $\| X_k \| = {\cal O}_P( \sqrt{NT} )$
as $N$ and $T$ grow at the same rate.
Thus, if the convergence rate of $\widehat \beta_R$ is faster than $\sqrt{N}$, i.e.
$\| \widehat \beta_R - \beta^0 \| = o_P(\sqrt{N})$, then we have
$\| e \|  \gg \left\|  (\widehat \beta_R - \beta^0) \cdot X \right\|$ asymptotically, and we expect
the extra $\widehat F_R$ to be dominated by $e$. A crucial step in the derivation
of Theorem~\ref{th:MAIN} is therefore to show faster than $\sqrt{N}$ convergence of $\widehat \beta_R$.
Conversely, we expect counter examples to the main result to be such that the convergence rate
of the estimator $\widehat \beta_R$ is not faster than $\sqrt{N}$, and we provide such a counter example
-- which, however, violates Assumptions~\ref{ass:LL} -- in Section~\ref{sec:AsymptoticSummary} below.
Whether the intuition about ``inclusion of irrelevant regressors'' carries over to the
``inclusion of irrelevant factors'' thus crucially depends on the convergence rate of $\widehat \beta_R$.

\section{Asymptotic Theory and Discussion}
\label{sec:AsyTheory}

Here we introduce key intermediate results for the proof of the main
Theorem~\ref{th:MAIN} stated above. These intermediate results may be useful independently of the main result, e.g. Moon and Weidner~\cite*{MoonWeidner2013} and Moon, Shum, and Weidner~\cite*{MoonShumWeidner2014}
crucially use the results established in Section~\ref{sec:expansion} for the case of known $R=R^0$.
The assumptions introduced below are all implied by the low-level Assumptions~\ref{ass:LL}
above, see to Lemma~\ref{lemma:JustifyEV} in the appendix.


\subsection{Consistency of $\widehat \beta_R$}
\label{sec:consistency}

Here we present a consistency result for $\widehat \beta_R$ under an arbitrary
asymptotic $N,T \rightarrow \infty$, i.e. without the assumption that $N$ and $T$ grow at the same rate,
which is imposed everywhere else in the paper. In addition to Assumption~\ref{ass:NC} we require
the following high level assumptions to obtain the result.

\begin{SNassumption}[\bf Spectral Norm of $X_k$ and $e$]~
\label{ass:SN}
\begin{itemize}
  \item [(i)] $\|X_k\| = {\cal O}_P(\sqrt{NT})$, \quad $k=1,\ldots,K$.
  \item [(ii)] $\| e \| = {\cal O}_P(\sqrt{\max(N,T)})$.
\end{itemize}
\end{SNassumption}


\begin{EXassumption}[\bf Weak Exogeneity of $X_k$]~
  \label{ass:EX}
  $\frac 1 {\sqrt{NT}} {\rm Tr}(X_k e^{\prime}) = {\cal O}_P(1)$,
              \quad $k=1,\ldots,K$.
\end{EXassumption}

\begin{theorem}
   \label{th:consistency}
      Let Assumptions \ref{ass:SN}, \ref{ass:EX} and \ref{ass:NC} be satisfied
      and let $R \geq R^0$. For $N,T \rightarrow \infty$ we then have
   $\sqrt{\min(N,T)} \left( \widehat \beta_R - \beta^0 \right) = {\cal O}_P(1)$.
\end{theorem}


\paragraph{Remarks}

\begin{itemize}
  \item[(i)] One can justify Assumption~\ref{ass:SN}$(i)$ by use of the norm inequality
$\|X_k\| \leq \|X_k\|_{HS}$ and the fact that
$\|X_k\|^2_{HS} = \sum_{i,t} X_{k,it}^2 = {\cal O}_P(NT)$, where the last step
follows e.g. if $X_{k,it}$ has a uniformly bounded second moment.
  \item[(ii)] Assumption \ref{ass:SN}$(ii)$ is a condition on the largest
eigenvalue of the random covariance matrix $e'e$, which is
often studied in the literature on random matrix theory, e.g.
Geman \cite*{Geman1980}, Bai, Silverstein, Yin \cite*{BaiSilvYin1988},
Yin, Bai, and Krishnaiah \cite*{BaiKrishYin1988}, Silverstein \cite*{Silverstein1989}.
The results in Latala \cite*{Latala2005} show that
$\| e \| = {\cal O}_P(\sqrt{\max(N,T)})$ if $e$ has independent entries with
mean zero and uniformly bounded fourth moment. Weak dependence of the
entries $e_{it}$ across $i$ and over $t$ is also permissible, see Appendix~\ref{app:SpectralNorm}

    \item[(iii)]
    Assumption \ref{ass:EX} requires exogeneity of the regressors
$X_k$, allowing for pre-determined regressors, and some weak dependence of $X_{k,it}e_{it}$ across $i$ and over $t$.\footnote{
Note that $\frac 1 {\sqrt{NT}} {\rm Tr}(X_k e^{\prime})= \frac 1 {\sqrt{NT}} \sum_i \sum_t X_{k,it}e_{it}$.
}

  \item [(iv)] The theorem
imposes no restriction at all on $f^0$ and $\lambda^0$, apart from the
condition $R \geq {\rm rank}(\lambda^0 f^{0 \prime})$.\footnote{This is the main reason why we use a slightly different
non-collinearity Assumption~\ref{ass:NC}, which avoids mentioning $\lambda^0$,
compared to Bai~\cite*{Bai2009}.}
  In particular, the strong factor Assumption~\ref{ass:SF} is not imposed here, i.e. consistency
  of $\widehat \beta_R$ holds independently of whether the factors
  are strong, weak, or not present at all. This is an  important robustness result, which
  is new in the literature.

  \item[(v)] Under an asymptotic where $N$ and $T$ grow at the same rate, which is imposed everywhere
  else in the paper, Theorem~\ref{th:consistency} shows $\sqrt{N}$ (or equivalently $\sqrt{T}$)
consistency of the estimator $\widehat \beta_R$.  To prove the consistency, we do not use the argument of the standard consistency proof for an extremum estimator which is to apply a uniform law of large numbers to the sample objective function to find the limit function that is uniquely minimized at the true parameter. Deriving the uniform limit of the objective function ${\cal L}_{NT}^{R^0}(\beta)$ is difficult. In the proof that is available in the supplementary appendix, we find a lower bound of the objective function ${\cal L}_{NT}^{R^0}(\beta)$ that is quadratic in $\beta - \beta^0$ asymptotically and establish the desired consistency, extending
the consistency proof in Bai~\cite*{Bai2009}.

\item[(vi)] $\sqrt{N}$ consistency of $\widehat \beta_R$
 implies that the residuals $Y - \widehat \beta_R \cdot X$ will be asymptotically close to
$\lambda^0 f^{0 \prime} + e$.\footnote{In the sense that
$ \| (Y - \widehat \beta_R \cdot X) - (\lambda^0 f^{0 \prime} + e) \| = \| (\widehat \beta_R - \beta) \cdot X \| = {\cal O}_P(\sqrt{N})$.}
 This allows consistent estimation
of $R^0$ under a strong factor Assumption~\ref{ass:SF} by employing the known techniques  on factor
models without regressors (by applying, e.g.~,
Bai and Ng~\cite*{BaiNg2002} to  $Y - \widehat \beta_R \cdot X$), as also
discussed in Bai~\cite*{Bai2009supp}.\footnote{
Bai~\cite*{Bai2009supp} does not prove the required
consistency and convergence rate of $\widehat \beta_R$, for $R > R^0$.}

\item[(vii)]
Having  a consistent estimator for $R^0$, say $\widehat{R}$, one can calculate $\widehat \beta_{\widehat{R}}$, which will be
asymptotically equal to $\widehat \beta_{R^0}$.
In practice, however, the finite sample properties of the estimator $\widehat \beta_{\widehat{R}}$ crucially depend on the finite sample properties of $\widehat{R}$. Many recent papers have documented difficulties in obtaining reliable estimates for $R^0$ at finite sample (see, e.g., the simulation results of Onatski~\cite*{Onatski2010} and Ahn and Horenstein~\cite*{AhnHorenstein2013}), and those difficulties
are also illustrated by our empirical example in Section~\ref{sec:Empirical}.

\end{itemize}





\subsection{Quadratic Approximation of ${\cal L}_{NT}^0(\beta) (:= {\cal L}_{NT}^{R^0}(\beta))$}
\label{sec:expansion}




To derive the limiting distribution of $\widehat \beta_R$,
we study the asymptotic properties of the
profile objective function ${\cal L}_{NT}^R(\beta)$ around $\beta^0$.
The expression in \eqref{LSobjective} cannot easily be discussed by
analytic means, since  no explicit formula for the
eigenvalues of a matrix is available.
In particular, a standard Taylor expansion
of ${\cal L}^R_{NT}(\beta)$ around $\beta^0$ cannot easily be derived.
Here, we consider the case of known $R=R^0$
and we perform a joint expansion
of  the corresponding profile objective function ${\cal L}^0_{NT}(\beta)$
 in the regression parameters $\beta$ and in the idiosyncratic error terms $e$. To perform this joint
expansion we apply the perturbation theory of linear operators (e.g., Kato~\cite*{Kato}).
We thereby obtain an approximate quadratic expansion of ${\cal L}_{NT}^0(\beta)$ in $\beta$,
which can be used to derive the first order asymptotic theory of
the LS estimator $\widehat \beta_{R^0}$, see Appendix~\ref{app:ExpansionDiscussion} for details.
In addition to the $K \times K$ matrix $W$ already defined in Section~\ref{sec:main} we now also define
\begin{align}
        C^{(1)}_k &=
         \frac 1 {\sqrt{NT}} \, {\rm Tr}( M_{\lambda^0} \, X_k \,
                  M_{f^0} \, e^{\prime} ) \; ,  \notag \\
        C^{(2)}_k  &=
        - \, \frac 1 {\sqrt{NT}} \, \bigg[
       {\rm Tr}\left(e M_{f^0} \, e' \, M_{\lambda^0} \, X_k \,
              f^0 \, (f^{0\prime}f^0)^{-1} \, (\lambda^{0\prime}\lambda^0)^{-1} \, \lambda^{0\prime} \right)
    \nonumber \\ & \qquad \qquad \quad
       +{\rm Tr}\left(e^{\prime}M_{\lambda^0} \, e \, M_{f^0} \, X^{\prime}_k \,
              \lambda^0 \, (\lambda^{0\prime}\lambda^0)^{-1} \, (f^{0\prime}f^0)^{-1} \, f^{0\prime} \right)
    \nonumber \\ & \qquad \qquad \quad
       +{\rm Tr}\left(e^{\prime}M_{\lambda^0} \, X_k \, M_{f^0} \, e^{\prime}
                \, \lambda^0 \, (\lambda^{0\prime}\lambda^0)^{-1} \, (f^{0\prime}f^0)^{-1} \, f^{0\prime} \right)
                        \bigg]  \; .
\end{align}
Let $C^{(1)}$ and $C^{(2)}$ be the $K$-vectors with elements $C^{(1)}_k$ and $C^{(2)}_k$, respectively.



\begin{theorem}
   \label{th:expansion}
   Let Assumptions  \ref{ass:SF} and \ref{ass:SN} be satisfied.
   Suppose that $ N,T \rightarrow \infty$ with $ N/T \rightarrow \nolinebreak[4] \kappa^2$,
   $0<\kappa<\infty$.
   Then we have
   \begin{align*}
      {\cal L}_{NT}^{0}(\beta) &= {\cal L}_{NT}^{0}(\beta^0)
                           -  \ft 2 {\sqrt{NT}}   \left(\beta-\beta^0 \right)'
                           \left( C^{(1)} + C^{(2)} \right)
                    + \left(\beta-\beta^0 \right)'  W  \left(\beta-\beta^0 \right)
                          + {\cal L}_{NT}^{0,{\rm rem}}(\beta)  ,
   \end{align*}
   where the remainder term ${\cal L}_{NT}^{0,{\rm rem}}(\beta)$
    satisfies for any sequence $c_{NT}\rightarrow 0$
   \begin{align*}
     \sup_{\{\beta :\left\| \beta -\beta^{0} \right\| \leq c_{NT}\}} \frac{
        \left| {\cal L}_{NT}^{0,{\rm rem}}(\beta) \right|  } { \left( 1 + \sqrt{NT} \, \left\| \beta -\beta^{0} \right\| \right)^2 } = o_{p}\left( \frac 1 {NT} \right) .
    \end{align*}
\end{theorem}

The bound on remainder\footnote{
The expansion in Theorem~\ref{th:expansion} contains a term that is linear  in $\beta$
and linear in $e$ ($C^{(1)}$ term), a term that is linear in $\beta$  and quadratic in $e$ ($C^{(2)}$ term),
and a term that is quadratic in $\beta$ ($W$ term). All higher order terms of the expansion
are contained in the remainder term ${\cal L}_{NT}^{0,{\rm rem}}(\beta)$.}
in  Theorem~\ref{th:expansion} is such that it has no effect on the first order asymptotic theory of
$\widehat \beta_{R^0}$, as stated in the following corollary (see also  Andrews~\cite*{Andrews1999}).

\begin{corollary}
   \label{cor:LimitR0}
   Let Assumptions \ref{ass:SF}, \ref{ass:SN}, \ref{ass:EX} and \ref{ass:NC}
   be satisfied.  In the limit $N,T \rightarrow \infty$ with
   $N/T \rightarrow \kappa^2$, $0< \kappa < \infty$, we then have
                  $\sqrt{NT}\left(\widehat \beta_{R^0} - \beta^0\right)
                   = W^{-1} \left( C^{(1)} + C^{(2)} \right) + o_P \left(1 + \| C^{(1)} \| \right)$.
   If we furthermore assume that $C^{(1)} = {\cal O}_P(1)$, then we obtain
    \[\sqrt{NT}\left(\widehat \beta_{R^0} - \beta^0\right)
                   = W^{-1} \left( C^{(1)} + C^{(2)} \right) + o_P(1)
                   = {\cal O}_P(1).\]
\end{corollary}

Note that our assumptions already guarantee
$C^{(2)}={\cal O}_P(1)$ and that $W$ is invertible with $W^{-1}={\cal O}_P(1)$, so this need not be
explicitly assumed in Corollary~\ref{cor:LimitR0}.


\paragraph{Remarks}

\begin{itemize}
  \item[(i)] More details on the expansion of  ${\cal L}_{NT}^{0}(\beta)$ are provided in
Appendix~\ref{app:ExpansionDiscussion} and the formal proofs can be found in
in Section~\ref{app:expansion1} of the supplementary appendix.

  \item[(ii)] Corollary~\ref{cor:LimitR0} allows to replicate the results in
Bai~\cite*{Bai2009}
and Moon and Weidner~\cite*{MoonWeidner2013}
on the asymptotic distribution of $\widehat \beta_{R^0}$,
including the result in formula~\eqref{AsyDistribution} above.\footnote{
Let $\rho$, $D(.)$, $D_0$, $D_Z$, $B_0$ and $C_0$ be the notation used in Assumption~A
and Theorem~3 of Bai~\cite*{Bai2009}, and let Bai's assumptions be satisfied.
Then, our $\kappa$, $W$, $C^{(1)}$ and $C^{(2)}$ satisfy
$\kappa=\rho^{-1/2}$, $W= D(f^0) \rightarrow_p D>0$, $C^{(1)} \rightarrow_d  {\cal N}(0,D_Z)$
and $W^{-1} C^{(2)} \rightarrow_p \rho^{1/2} B_0 + \rho^{-1/2} C_0$.
Corollary~\ref{cor:LimitR0} can therefore be used to replicate Theorem~3 in Bai~\cite*{Bai2009}.
For more details and extensions of this we refer to Moon and Weidner~\cite*{MoonWeidner2013}.
}
The assumptions of the corollary do not restrict the regressors
to be strictly exogenous and do not impose Assumption~\ref{ass:LL}.


  \item[(iii)] If one weakens
Assumption \ref{ass:SN}$(ii)$ to $\|e\|=o_P(N^{2/3})$,
then Theorem \ref{th:expansion} still continues to hold.
If $C^{(2)}={\cal O}_P(1)$,
then Corollary~\ref{cor:LimitR0} also holds under this weaker condition on $\|e\|$.
\end{itemize}

\subsection{Remarks on Deriving the Convergence Rate and Asymptotic Distribution of  $\widehat \beta_R$ for $R>R^0$.}
\label{sec:AsymptoticSummary}


\subsubsection*{An example that motivates stronger restrictions}

The results in Bai~\cite*{Bai2009} and
Corollary~\ref{cor:LimitR0} above show that under appropriate assumptions
the estimator
$\widehat \beta_R$ is $\sqrt{NT}$-consistent for $R=R^0$.
For $R>R^0$
we know from Theorem~\ref{th:consistency} that $\widehat \beta_R$ is $\sqrt{N}$ consistent
as $N$ and $T$ grow at the same rate, but we have not shown faster than $\sqrt{N}$ converge of
$\widehat \beta_R$ for $R>R^0$, yet, which according to the heuristic discussion at the end of
Section~\ref{sec:main} is a very important intermediate step to obtain our main result.\footnote{
One reason why  $\widehat \beta_R$ might only converge at $\sqrt{N}$ rate, but not faster, are
weak factors (both for $R>R^0$ and for $R=R^0$). A weak factor
(see e.g. Onatski~\cite*{Onatski2010,Onatski2012} and Chudik, Pesaran and Tosetti~\cite*{ChudikPesaranTosetti2011})
might not be picked up at all
or might only be estimated very inaccurately by the principal components estimator $\widehat F_R$,
in which case that factor is not properly accounted for in the LS estimation procedure. If this happens
and the weak factor is correlated with the regressors, then there is some uncorrected weak
endogeneity problem, and $\widehat \beta_R$ will only converge at $\sqrt{N}$ rate.
We do not consider the issue of weak factors any further in this paper.
}
However, one might not obtain a faster than $\sqrt{N}$ convergence rate of $\widehat \beta_R$
for $R>R^0$ without imposing further restrictions, as the following example shows.

\begin{example}
     \label{example:negative}
      Let $R^0=0$ (no true factors) and $K=1$ (one regressor).
      The true model reads $Y_{it} = \beta^0 X_{it} + e_{it}$, and we consider the following data generating process (DGP)
  \begin{align*}
   X_{it} &= a \widetilde X_{it} + \lambda_{x,i} f_{x,t}, &
   e &= \left( \mathbbm{1}_{N} + c \, \frac{\lambda _{x}\lambda _{x}^{\prime }}{N} \right)
   u \left( \mathbbm{1}_{T}+c \, \frac{f_{x}f_{x}^{\prime }}{T}\right) ,
\end{align*}
where $e$ and $u$ are $N \times T$ matrices with entries $e_{it}$ and $u_{it}$, respectively,
and $\lambda _{x}$ is an $N$-vector with entries $\lambda_{x,i}$, and $f_x$ is a $T$-vector with entries $f_{x,t}$.
Let $\widetilde X_{it}$ and $u_{it}$ be mutually independent
iid standard normally distributed random variables.
Let $\lambda_{x,i} \in {\cal B}$ and $ f_{x,t} \in {\cal B}$ be
non-random sequences with bounded range ${\cal B} \subset \mathbbm{R}$ such that
$\frac 1 N \sum_{i=1}^N  \lambda^2_{x,i} \rightarrow 1$
and $\frac 1 T \sum_{t=1}^T  f^2_{x,t}  \rightarrow 1$ asymptotically.\footnote{
    We could also allow $\lambda_x$ and $f_x$ to be random (but independent of $e$
    and $\widetilde X$) and we could let the range of ${\cal B}$ be unbounded.
    We only assume non-random $\lambda_x$ and $f_x$ to guarantee that the DGP
    satisfies Assumption~D of  Bai~\cite*{Bai2009}, namely that
    $X$ and $e$ are independent (otherwise we only have mean-independence, i.e. $\mathbbm{E}(e|X)=0$).
    Similarly, we only assume bounded ${\cal B}$ to satisfy
    the restrictions on $e_{it}$ imposed in
    Assumption~C of Bai~\cite*{Bai2009}.
}
Consider $N,T \rightarrow \infty$ such that
$N/T \rightarrow \kappa^2$, $0< \kappa < \infty$,
and let $0<a <(1/2)^{2/3} \min(\kappa^2,\kappa^{-2})$ and $c \geq    \frac{ (2+\sqrt{2})  \left( 1+\kappa \right)
 (1+ \sqrt{3} a^{-1/4} )}   { \min(1,\kappa) [1/2 - a^{3/2} \max(\kappa,\kappa^{-1})]} $.\footnote{
 The bounds on the constants $a$ and $c$ imposed here are sufficient, but not necessary
 for the result of no faster than $\sqrt{N}$ convergence of $\widehat \beta_1$.
Simulation evidence suggests that this result holds for a much larger range of $a$, $c$ values.
}
For this DGP one can show that
$\widehat \beta_1$, the LS-estimator with $R=1>R^0$, only converges at a rate of $\sqrt{N}$ to $\beta^0$, but not faster.
\end{example}

The proof of the last statement is provided in the supplementary material.
The DGP in this example satisfies
      all the assumptions imposed in Corollary~\ref{cor:LimitR0} to derive the
  limiting distribution of the LS-estimator for $R=R^0$, including $\sqrt{NT}$-consistency of
  $\widehat \beta_R$ for $R=R^0$ (=0 in this example).
    It also satisfies all the regularity conditions imposed in
   Bai~\cite*{Bai2009}.\footnote{See Section~\ref{app:CheckBai} in the supplementary material  for details.}
The aspect that is special about this DGP is
  that $\lambda_x$ and $f_x$ feature both in $X_{it}$
  and in the second moment structure
   of $e_{it}$.
   The heuristic discussion at the end of Section~\ref{sec:main} provides some intuition why this
   can be problematic, because the leading principal components obtained from only the error matrix $e$
   will have a strong sample correlation with $X_{it}$ for this DGP.

\subsubsection*{Faster than $\sqrt{N}$ convergence of $\widehat \beta_R$}

 In Appendix~\ref{sec:ConvergenceRate}, we summarize our results on faster than $\sqrt{N}$ convergence of $\widehat \beta_R$ for $R \geq R^0$.  The above example shows that this requires more restrictive
 assumptions than those imposed for the analysis of the case $R=R^0$ above,
 but the assumptions that we impose for this intermediate results are
 still significantly weaker than the Assumption~\ref{ass:LL} required for our main result above,
   in particular either cross-sectional correlation or time-serial correlation of
   $e_{it}$ are still allowed.

In that appendix we also provide one set of assumptions (Assumption~\ref{ass:DX-2})
   for faster than $\sqrt{N}$ convergence such that
     no additional conditions on $e$ are required, but where the regressors are restricted
   to essentially be  lagged dependent variables in an AR(p) model with factors.

\subsubsection*{On the role of the iid normality of $e_{it}$}

We establish the asymptotic equivalence of $\widehat \beta_R$ and $\widehat \beta_{R^0}$ in Theorem \ref{th:MAIN}
by showing that the LS objective function $\mathcal{L}_{NT}^{R}(\beta)$ can, up to a constant,
 be uniformly well approximated
 by $\mathcal{L}_{NT}^{0}(\beta)$ in shrinking neighborhoods around the true parameter. For this, we need not only the faster than $\sqrt{N}$ convergence rate of $\widehat \beta_R$, but also require the Assumption~\ref{ass:EV} in Appendix~\ref{sec:Equivalence}. This is a high-level assumption on the eigenvalues and eigenvectors of the random covariance matrices
$E E'$ and $E' E$, where $E=M_{\lambda^0} e M_{f^0}$. The assumption essentially requires
the eigenvalues of those matrices to be sufficiently separated from each other,
as well as the eigenvectors of those matrices to be
sufficiently uncorrelated with the regressors $X_k$, and with $e P_{f^0}$ and $P_{\lambda^0} e$.

We use the iid normality of $e_{it}$ to verify those high-level conditions in
Section~\ref{ass:SufficiencyLL} of the supplementary appendix.
There are three reasons why we can currently only
verify those conditions for iid normal errors:


\begin{itemize}
    \item[(i)] The random matrix theory literature studies the eigenvalues and eigenvectors of random
    covariance matrices of the form $e e'$ and $e' e$, while we have to deal with the additional projectors $M_{\lambda^0}$
    and $M_{f^0}$ in the random covariance matrices.
    These additional projections stem from integrating out the true factors and factor loadings of the model.
If the error distribution is $iid$ normal, and independent from $\lambda^0$ and $f^0$, then these projections
are unproblematic, since the distribution of $e$ is rotationally invariant from the left and right in that case,
so that the projections are mathematically equivalent to a reduction of the sample size
by $R^0$ in both panel dimensions.

    \item[(ii)]
    In the iid normal case one can furthermore use the invariance of the distribution of $e$ under orthonormal
    rotations from the left and from the right to also fully
    characterize the distribution of the eigenvectors of $E E'$ and $E E'$.\footnote{Rotational invariance implies
    that the distribution of the normalized eigenvectors
    is given by the Haar measure of a rotation group manifold.}
      The conjecture in the random matrix theory literature is that
the limiting distribution of the eigenvectors of a random covariance matrix
is ``distribution free'',
i.e. is independent of the particular distribution of $e_{it}$, see, e.g., Silverstein~\cite*{Silverstein1990} and
Bai~\cite*{bai1999review}. However, we are not currently aware of a formulation
and corresponding proof of this conjecture that is sufficient for our purposes, i.e. that would allow us to
verify our high-level Assumption~\ref{ass:EV} more generally.

    \item[(iii)]  We also require certain properties of the eigenvalues of $E E'$ and $E E'$.
    Eigenvalues are studied more intensely than eigenvectors
in the random matrix theory literature,
and it is well-known that the properly normalized
empirical distribution of the eigenvalues
(the so called empirical spectral distribution)
of an $iid$ sample covariance matrix converges
to the Mar{\v{c}}enko-Pastur-law (Mar{\v{c}}enko and Pastur~\cite*{MarcenkoPastur1967})
for asymptotics where $N$ and $T$ grow at the same rate.
This result does not require normality, and results on the limiting spectral distribution are also
known for non-iid matrices. However,
to check our high-level Assumption~\ref{ass:EV}
we also need results on the
convergence rate of the empirical spectral
distribution to its limit law, which is an ongoing research subject in the
literature, e.g. Bai~\cite*{Bai1993}, Bai, Miao and Yao \cite*{BaiMiaoYao2004},
G{\"o}tze and Tikhomirov \cite*{GotzeTikhomirov2010}, and we are currently only aware of
results on this convergence rate for the case of either iid or iid normal errors.
To verify the high-level assumption we furthermore
use a result from
Johnstone~\cite*{Johnstone2001} and Soshnikov~\cite*{Soshnikov2002}
that shows that the properly normalized few
largest eigenvalues of  $E E'$ and $E E'$ converge to the Tracy-Widom law, and to our knowledge this result
is not established for error distributions that are not iid normal.
\end{itemize}

In spite of these severe mathematical challenges,
we believe that in principle our high-level Assumption~\ref{ass:EV} could be verified for more
general error distributions, implying that our main result of asymptotic equivalence of $\widehat \beta_R$
and $\widehat \beta_{R^0}$ holds more generally. This is also supported by our Monte Carlo
simulations, where we explore non-independent and non-normal error distributions.



\section{Empirical Illustration}
\label{sec:Empirical}

As an illustrative empirical example, we estimate the dynamic effects of
unilateral divorce law reforms on the state-wise divorce rates in the US. The impact of the
divorce law reform has been studied by many researches (e.g., Allen~\cite*{Allen1992},
Peters~\cite*{Peters1986,Peters1992}, Gray~\cite*{Gray1998}, Friedberg~\cite*{Friedberg1998},
Wolfers~\cite*{Wolfers2006}, and Kim and Oka~\cite*{KimOka2014}). In this section we revisit this topic,
extending Wolfers~\cite*{Wolfers2006} and Kim and Oka~\cite*{KimOka2014} by controlling
for interactive fixed effects and also a lagged dependent variable.

Let $Y_{it}$ denote the number of divorces per 1000 people in state $i$
at time $t$, and let $D_{i}$ denote the year in which state $i$ introduced the unilateral
divorce law, i.e. before year $D_i$ state $i$ had a consent divorce law, while from $D_i$ onwards state $i$
had  a unilateral ``no-fault'' divorce law, which loweres the barrier for divorce.
The goal is to estimate the dynamic effects of this law change on the divorce rate. The empirical model we estimate is
\begin{equation}
Y_{it}=\beta_{0} \, Y_{i,t-1} + \sum_{k=1}^{8} \beta _{k}  X_{k,it}
+\alpha_{i}+ \gamma_i \, t + \delta_i \, t^{2} + \mu_t +\lambda _{i}^{\prime }f_{t}+e_{it},
\label{m:empirical.ex}
\end{equation}
where we follow Wolfers~\cite*{Wolfers2006} in defining the regressors as bi-annual dummies:
\begin{eqnarray*}
X_{k,it} &=&\mathbbm{1}\{D_{i}+2(k-1)\leq t\leq D_{i}+2k-1\},
\quad \text{for} \; \; k=1,...,7, \\
X_{8,it} &=&\mathbbm{1}\{D_{i}+2(k-1)\leq t\}.
\end{eqnarray*}
The dummy variable and quadratic trend specification $\alpha_{i}+ \gamma_i \, t + \delta_i \, t^{2} + \mu_t$
is also used in Friedberg~\cite*{Friedberg1998} and Wolfers~\cite*{Wolfers2006}. The additional
interactive fixed effects $\lambda _{i}^{\prime }f_{t}$ were added in Kim and Oka~\cite*{KimOka2014}
to control for additional unobserved heterogeneity in the divorce rate, e.g. due to social, cultural
or demographic factors. We extend the specification further by
adding a lagged dependent variable $Y_{i,t-1}$ to control for state dependence of the divorce rate,
but we also report results without $Y_{i,t-1}$ below.
We use the dataset of Kim and Oka~\cite*{KimOka2014},\footnote{
The data is available from {\tt http://qed.econ.queensu.ca/jae/2014-v29.2/kim-oka/}}
which is a balanced panel of $N=48$ states over $T=33$ years,
leaving $T=32$ time periods if the lagged dependent variable is included.

For estimation we first eliminate $\alpha_{i}$, $\gamma_i$, $\delta_i$ and $\mu_t$ from the model by projecting the
outcome variable and all regressors accordingly, e.g. $\widetilde Y = M_{1_N} Y M_{(1_T, {\bf t}, {\bf t}^2)}$,
where $1_N$ and $1_T$ are $N$- and $T$-vectors, respectively, with all entries equal to one,
and ${\bf t}$ and ${\bf t}^2$ are $T$-vectors with entries $t$ and $t^2$, respectively.
The model after projection reads
$\widetilde Y_{it}=\beta_{0} \, \widetilde Y_{i,t-1} + \sum_{k=1}^{8} \beta _{k}  \widetilde X_{k,it}
  +\widetilde \lambda _{i}^{\prime } \widetilde f_{t} + \widetilde e_{it}$, which is exactly the model we have studied
  so far in this paper.\footnote{To construct $\widetilde Y_{i,t-1}$ we first apply the lag-operator and then apply the
  projections $M_{1_N}$ and $M_{(1_T, {\bf t}, {\bf t}^2)}$.}
   We use the LS estimator described above to estimate this model.
The projection reduces the effective sample size to $N=48-1=47$ and $T=32-3=29$, which
should be accounted for when calculating standard errors, e.g. in the formula for
$\widehat \sigma_R^2$ above (degree of freedom correction). Our theoretical results are still applicable.\footnote{
If $e_{it}$ is iid normal, then $\widetilde e_{it}$ is not, but one can apply appropriate orthogonal rotations in $N$- and $T$-space such that $\widetilde e_{it}$ becomes iid normal again, although with sample size reduced to $N=47$ and $T=29$.
The rotation has no effect on the LS estimator, i.e. it does not matter whether we work in the original or the rotated
frame.
}

We need to decide on a number of factors $R$ when implementing the LS estimator. As already mentioned in
the last remark  in Section~\ref{sec:consistency} above, we can can apply known techniques from the literature
on factor models without regressors to obtain a consistent estimator of $R^0$. To do so we choose a maximum
number of factors of $R_{\max}=9$ to obtain the preliminary estimate $\widehat \beta_{R_{\max}}$
and then calculate the residuals
$\widehat u_{it} = \widetilde Y_{it} - \widehat \beta_{R_{\max},0} \, \widetilde Y_{i,t-1}
  - \sum_{k=1}^{8} \widehat \beta _{R_{\max},k}  \widetilde X_{k,it}$.
We then apply the IC, PC and BIC3 criteria of Bai and Ng~\cite*{BaiNg2002},\footnote{
Following Onatski~\cite*{Onatski2010} and Ahn and Horenstein~\cite*{AhnHorenstein2013}
we report only BIC3 among the AIC and BIC criteria of Bai and Ng~\cite*{BaiNg2002}.
}
the criterion described in Onatski~\cite*{Onatski2010}, and the ER and GR criteria of
Ahn and Horenstein~\cite*{AhnHorenstein2013} to $\widehat u$.\footnote{
To include $R=0$ as a possible outcome for the Ahn and Horenstein (2013) criterion, we use the mock
eigenvalue used in their simulations.}
Most of these criteria also require specification
of $R_{\max}$, and we continue to use $R_{\max} = 9$.
The corresponding estimation results for $R$ are presented in Table~\ref{table:facor.number.est}.
In addition, we also report the log scree plot, i.e. the sorted eigenvalues of $\widehat u' \widehat u$
in Figure 1.

The log scree plot already shows that it is not obvious how to decompose the eigenvalue spectrum
into a few larger eigenvalues stemming from factors and the remaining smaller eigenvalues stemming from
the idiosyncratic error term.\footnote{The first largest eigenvalue is 2.2 times larger than the second eigenvalue,
 the second is 1.6 times larger than the third, the third is 1.9 times larger than fourth. So the largest view eigenvalues
 are larger than the remaining ones, and the strong factor assumption might not be completely inappropriate here.
 However, deciding on a cutoff between factor and non-factor eigenvalues is difficult.}
 This problem is also reflected in the very different estimates for $R$ that one obtains from the various criteria.
 It might appear that IC1, IC3, PC1, PC2 and PC3 all agree on $\widehat R=9$, but this is simply $\widehat R=R_{\max}$,
 and if we choose $R_{\max}=10$, then all these criteria deliver $\widehat R=10$, so this should not be considered
 a reliable estimate.

\begin{table}[tb!]
\begin{minipage}{\textwidth}

\begin{minipage}[b]{0.47\textwidth}
\begin{center}
\begin{tabular}{c@{\;}c|c@{\;}c|c@{\;}c}
Criterion: & $\widehat{R}$ & Criterion: & $\widehat{R}$ & Criterion: & $\widehat{R}$
\\
\hline \hline
IC1: & $9$ & PC1: & $9$ & Onatski: & 1 \\
IC2: & $7$ & PC2: & $9$ & ER: & 1 \\
IC3: & $9$ & PC3: & $9$ & GR: & 3 \\
BIC3: & 6 &  &  &  &
\end{tabular}
\end{center}
\caption{\label{table:facor.number.est} \footnotesize
Estimated number of factors in the residuals $\widehat u$,
using different criteria for estimation and $R_{\max}=9$. The  IC, PC and BIC criteria
are described in Bai and Ng~(2002), the ER and GR criteria are from Ahn and Horenstein (2013),
and we also use the criterion of Onatski~(2010).
}

\end{minipage}
\hfill
\begin{minipage}[b]{0.47\textwidth}
\begin{center}
\includegraphics[width=0.9\textwidth]{./TablesFigures/LogScreePlot.eps}
\end{center}
\vspace{-0.3cm}
Figure~1: {\footnotesize
Log scree plot.
The natural logarithm of the sorted eigenvalues (corresponding to the principal components, or factors)
of $\widehat u' \widehat u$ are plotted.
}

\end{minipage}
\end{minipage}
\end{table}

\begin{table}[tb!]
\centering
\includegraphics[width=\textwidth]{./TablesFigures/Beta-Estimates-Lag.eps}
\caption{\label{tab:EstimateBetaLag}\footnotesize
Dynamic effects of divorce law reform. We report
bias corrected LS-estimates for the regression coefficients in model~\eqref{m:empirical.ex}. Each column corresponds to
a different number of factors $R \in \{0,1,\ldots,9\}$ used in the estimation. t-values are reported in parenthesis.}
\end{table}

\begin{table}[htb!]
\centering
\includegraphics[width=\textwidth]{./TablesFigures/Beta-Estimates-No-Lag.eps}
\caption{\label{tab:EstimateBetaNoLag}\footnotesize
Same as Table~\ref{tab:EstimateBetaLag}, but without including the lagged dependent variable into the model.}
\end{table}

On the other hand, our asymptotic theory suggests, that the exact choice of $R$ in the estimation of
$\widehat \beta_R$ should not matter too much, as long as $R$ is chosen large enough to cover all relevant factors.
Table~\ref{tab:EstimateBetaLag} contains the estimation results for the bias corrected $\widehat \beta_R$ for $R \in \{0,1,\ldots,9\}$.
Table~\ref{tab:EstimateBetaNoLag} contains estimates if the lagged dependent variable is not included into the model.\footnote{
The result for $R=7$  in Table~\ref{tab:EstimateBetaNoLag} should be equal to column (6) in Table III
of Kim and Oka~\cite*{KimOka2014}. The discrepancy is explained
by a coding error in their bias computation.
Note also that the result for $R=0$ in Table~\ref{tab:EstimateBetaNoLag} does not match the one in Wolfers~\cite*{Wolfers2006},
because he uses WLS with state population weights, while we use OLS for simplicity.
Kim and Oka~\cite*{KimOka2014} estimate both WLS and OLS and find that the difference between the resulting
estimates becomes insignificant, once a sufficient number of interactive fixed effects is controlled for.}
For all reported estimates we perform bias correction and standard error estimation as described in
Bai~\cite*{Bai2009} and Moon and Weidner~\cite*{MoonWeidner2013}.\footnote{We correct for the biases
due to heterscedasticity in both panel dimensions worked out in Bai~\cite*{Bai2009},
 as well as for the dynamic bias worked out in Moon and Weidner~\cite*{MoonWeidner2013}.
 For the latter we use the formula for $ \widehat B_{R,k}$ above, with bandwidth $M=2$.
For the standard error estimation we allow for heterscedasticity in both panel dimensions, also
following Bai~\cite*{Bai2009} and Moon and Weidner~\cite*{MoonWeidner2013}.
The bias and standard error formulas in those paper assume $R=R^0$ known, but we strongly expect that those formulas
are robust towards $R>R^0$, as partly justified by Theorem~\ref{th:Estimators} above.
For the model without lagged dependent variable we also allow for serial correlation in $e_{it}$ when estimating the
bias and standard deviation of $\widehat \beta_R$.
}

When ignoring the lagged dependent variable coefficient,
one finds that in both Table~\ref{tab:EstimateBetaLag} and Table~\ref{tab:EstimateBetaNoLag}
the estimation results for $\widehat \beta_R$ and the corresponding t-values
are quite sensitive to changes in $R$ for very small values of $R$, but become
much more stable as $R$ increases, and actually do not change too much anymore from roughly $R=2$ onwards.
These findings are very well in line with our asymptotic theory, and the dynamic effect of divorce law reform
that we find are also similar to the findings in Wolfers~\cite*{Wolfers2006}
and Kim and Oka~\cite*{KimOka2014}. The effect of the law reform on the divorce rates initially increases over time,
is certainly significant in year 3-4 after the reform, and declines and becomes insignificant afterwards.\footnote{
The magnitude of the estimates is smaller than those in Wolfers~\cite*{Wolfers2006}, i.e. controlling
for unobserved factors reduced the effect size, as already pointed out by Kim and Oka~\cite*{KimOka2014}.}

In contrast,
the estimated coefficient on the lagged dependent variable in Table~\ref{tab:EstimateBetaLag} is quite large and
highly significant for small values of $R$, but decreases steadily with $R$, until it gets close to zero and
insignificant for $R \geq 8$. A plausible interpretation of this finding is that the model that includes the lagged dependent
variable is misspecified, and that the estimated value of $\beta_0$ for small values of $R$ does not correspond
to a true state dependence of $Y_{it}$, but simply reflects the time-serial correlation of the error process
being picked up by the autoregressive model.
According to this interpretation,
once we include more
and more factors into the model we control for more and more serial dependence of the unobserved error term,
thus uncovering the true insignificance of $\beta_0$ in the estimates for $R \geq 8$.

This empirical example shows that instead of relying on a single estimate $\widehat R$ for the number of factors
and reporting the corresponding $\widehat \beta_{\widehat R}$ it can be very informative to calculate
$\widehat \beta_{R}$ for multiple values of $R$. Whether the estimated coefficients become stable for sufficiently
large $R$ values, as our asymptotic theory suggests, is a useful robustness check for the model. When reporting the final results, then, it is better,
within a reasonable range, to choose
an $R$ that is too large than one that is too small.

We also perform a Monte Carlo
simulation that is tailored towards the empirical application.
For this we use the static model without lagged dependent variable.
To generate $Y_{it}$ in equation \eqref{m:empirical.ex} with $\beta_0=0$
we use the observed regressors $X_{k,it}$, as described above, and
as true parameters we use the
$\beta$ (bias corrected), $\alpha_i$, $\gamma_i$, $\delta_i$, $\mu_t$, $\lambda_i$ and $f_t$ obtained from the estimation with $R=4$
(i.e. $\beta_k$ as reported in the $R=4$ column of Table~\ref{tab:EstimateBetaNoLag}). We generate $e_{it}$
from an MA(1) model with $t(5)$ distributed
innovations. Note that this error distribution violates the
assumption~\ref{ass:LL}$(ii)$.

In this ``empirical Monte Carlo'' we have
$N=48$, $T=33$ and true number of factors $R^0=4$.
We find that the
bias corrected estimates for $\beta_k$, $k=1,\ldots,8$,
are essentially unbiased when $R \geq R^0$ factors are used in the estimation, but for $R<R^0$ the coefficient estimates are often biased.
For $\beta_k$, $k \geq 3$, there are only small changes in the
 standard deviation of the estimator between $R=4$ and $R=9$,
  but for $\beta_k$, $k=1,2$, we observe standard deviation inflation
of up to $25 \%$ between $R=4$ and $R=9$.
Given the relatively small sample size
the difference between $R=9$ and $R^0=4$ is relatively large, and some
finite sample inefficiency is not too surprising.
The detailed results are available in the supplementary appendix.




\section{Monte Carlo Simulations}
\label{sec:MC}

In addition to the ``empirical Monte Carlo'' discussed above
we now investigate the finite sample properties of $\widehat \beta_{R}$
and $\widehat \beta_{R}^{\rm BC}$  further. In the simulations in this section we use
a generated regressor $X_{it}$ that is correlated with the interactive fixed effects.
The serial correlation of the error term $e_{it}$ together with the
data generating process (DGP) for $X_{it}$, $\lambda_i$ and $f_t$ are such that
the naive LS estimator has an asymptotic bias. This allows to verify
whether the bias is essentially unchanged for $R>R^0$ and whether
bias correction works well for $R>R^0$ in finite sample. We also study
various combinations of $N$ and $T$.


The model is a static panel model with one regressor ($K=1$), two factors ($R^0=2$),
and the following DGP:
\begin{align}
   Y_{it} &= \beta^0 X_{it} + \sum_{r=1}^2 \lambda_{ir} f_{tr} + e_{it} ,
   \nonumber \\
   X_{it} &= 1 + \widetilde X_{it}
               + \sum_{r=1}^2 (\lambda_{ir}+\chi_{ir}) (f_{tr} + f_{t-1,r} ) ,
   \nonumber \\
   e_{it} &= \frac 1 {\sqrt{2}} (v_{it} +v_{i,t-1} )  .
   \label{DGP-Static}
\end{align}
The random variables
$\widetilde X_{it}$, $\lambda_{ir}$, $f_{tr}$, $\chi_{ir}$ and $v_{it}$
are mutually independent; with $\widetilde X_{it}$ and $f_{tr} \, \sim \, iid  \, {\cal N}(0,1)$;
$\lambda_{ir}$ and $\chi_{ir} \, \sim iid \, {\cal N}(1,1)$; and
$v_{it} \, \sim \, iid \, t(5)$, i.e.~$v_{it}$ has a Student's t-distribution  with 5 degrees of freedom.


Note that this model satisfies Assumptions~\ref{ass:SF},~\ref{ass:NC}, and ~\ref{ass:LL}(i), but not ~\ref{ass:LL}(ii).
The error term $e_{it}$ is \emph{not} distributed as $iid$ normal. The time series of $e_{it}$ follows an MA(1)
process with innovations distributed as $t(5)$.

We choose $\beta^0=1$, and
use $10,000$ repetitions in our simulation.
The true number of factors is chosen to be $R^0=2$.
For each draw of $Y$ and $X$
we compute the LS estimator $\widehat \beta_R$ according to equation \eqref{estimator}
for different values of $R$, namely $R \in \{0,1,2,3,4,5\}$.

\begin{table}[tb!]
\centering
\includegraphics[width=14cm]{./TablesFigures/table_bias_sd_static.eps}
\caption{\label{tab:MC-Static-1}\footnotesize
For different combinations of sample sizes $N$ and $T$ we report
the bias and standard deviation of the estimator $\widehat \beta_R$, for $R=0,1,\ldots,5$,  based on
simulations with $10,000$ repetition of design \eqref{DGP-Static}, where the true number of
factors is $R^0=2$.}
\end{table}

Table~\ref{tab:MC-Static-1} reports bias and standard deviation of the estimator $\widehat \beta_R$
for different combinations of $R$, $N$ and $T$. For $R<R^0=2$ the model is misspecified and $\widehat \beta_R$
turns out to be severely biased.
There is also bias in $\widehat \beta_R$ for $R \geq R^0$, due to time-serial correlation of $e_{it}$.
This bias was worked out in Bai~\cite*{Bai2009}, and bias correction is also discussed there.

\begin{table}[tb!]
\centering
\includegraphics[width=14cm]{./TablesFigures/table_quantile_static.eps}
\caption{\label{tab:MC-Static-2}\footnotesize
Quantiles of the distribution of $\sqrt{NT}( \widehat \beta_R - \beta^0 )$
are reported for
$N=T=100$ and $N=T=300$, with $R=2,3,4,5$,
 based on
simulations with $10,000$ repetition of design \eqref{DGP-Static}, where the true number of
factors is $R^0=2$.
 }
\end{table}

Table~\ref{tab:MC-Static-2} reports various quantiles of
the distribution of $\sqrt{NT}( \widehat \beta_R - \beta^0 )$
for $N=T=100$ and $N=T=300$, and different values of $R \geq R^0$.
From these tables, we see that as $N,T$ increases the distribution of $\widehat \beta_R$ gets closer to that of $\widehat\beta_{R^0}$.

\begin{table}[tb!]
\centering
\includegraphics[width=14cm]{./TablesFigures/table_size_static.eps}
\caption{\label{tab:MC-Static-3}\footnotesize
The empirical size of a t-test with $5 \%$ nominal size is reported for different combinations of $N$, $T$ and $R$,
based on $10,000$ repetition of design \eqref{DGP-Static}. A bias corrected estimator $\widehat \beta^{\rm BC}_R$ is used to calculate the
test statistics, and we allow for heteroscedasticity and time-serial correlation when estimating bias and standard
deviation. Results for $R=0,1$ are not reported since those have
size=1 due to misspecification.  }
\end{table}

Table~\ref{tab:MC-Static-3} reports the size of a t-test with nominal size equal to $5 \%$ for $R \geq R^0$.
We use the results in Bai~\cite*{Bai2009}
to correct for the leading $1/N$ (not actually present in our DGP)
and $1/T$ (present in our DGP) biases in $\widehat \beta_R$ before calculating the t-test statistics, allowing for heteroscedsticity in both panel dimensions and for time-serial correlation when estimating the bias and standard
deviation of $\widehat \beta_R$. The finite sample size distortions are mostly due to residual bias after bias correction,
but also partly due to some finite sample downward bias in the standard error estimates. The size distortions increase
with $R$, but for all values of $R \geq R^0$ in Table~\ref{tab:MC-Static-3} the size distortions decrease rapidly as $T$ increases.

Monte Carlo Simulation results for an AR(1) model with factors can be found in Section~\ref{app:MC}
of the supplementary material. Those additional simulations show that the finite sample
properties (e.g. for $T=30$) of $\widehat \beta_{R^0}$ and $\widehat \beta_{R}$, $R>R^0$,
can be quite different, but those differences vanish as $T$ becomes large, as predicted
by our asymptotic theory.
In general, we always expect some finite sample inefficiency from overestimating the number
of factors.








\section{Conclusions}
\label{sec:conclusions}

We show that under certain assumptions
the limiting distribution of the LS estimator of a linear panel regression with interactive fixed effects
does not change when we include redundant factors in the estimation. The implication of this is
that one can use an upper bound of the number of factors $R$ in the estimation without
asymptotic efficiency loss. However, some finite sample efficiency loss from overestimating $R$
is likely, so that $R$ should not be chosen too large in actual applications.
We  impose $iid$ normality of the regression errors to derive the asymptotic result, because we require certain
results on the eigenvalues and eigenvectors of random covariance matrices that are only known in that case.
We expect that progress in the literature on large dimensional
random covariance matrices will allow verification of our high-level assumptions
under more general error distributions,
and our Monte Carlo simulations suggest that the result also holds for non-normal
and correlated errors.
We also provide multiple intermediate asymptotic results under more general conditions.


\begin{appendix}