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.
104,065 characters
Tractable Estimation of Nonlinear Panels with Interactive Fixed Effects
\title{\textsc{Tractable Estimation of Nonlinear Panels\\ with Interactive Fixed Effects}\thanks{
We thank Aureo de Paula and the participants of the 2nd UCL–CeMMAP–IFS Ph.D. Econometrics Research Day~(2024) and UCL Econometrics Brownbag Seminar for their valuable comments. We thank Martin Weidner for sharing the codes and data from \citet{chen2021nonlinear}, and thank Wei Miao for careful testing of the package and for helpful comments on the implementation.
}}
\author{Andrei \textsc{Zeleneev}\thanks{University College London: \textsf{[email removed]}.} \quad \quad Weisheng \textsc{Zhang}
\thanks{
University College London: \textsf{[email removed]}.}
}
\maketitle
\begin{abstract}
Interactive fixed effects are routinely controlled for in linear panel models. While an analogous fixed effects (FE) estimator for nonlinear models has been available in the literature \citep{chen2021nonlinear}, it sees much more limited use in applied research because its implementation involves solving a high-dimensional non-convex problem. In this paper, we complement the theoretical analysis of \citet{chen2021nonlinear} by providing a new computationally efficient estimator that is asymptotically equivalent to their estimator. Unlike the previously proposed FE estimator, our estimator avoids solving a non-convex high-dimensional optimization problem and can be feasibly computed in large nonlinear panels. Our proposed method involves two steps. In the first step, we convexify the optimization problem using nuclear norm regularization (NNR) and obtain preliminary NNR estimators of the parameters, including the fixed effects. Then, we find the global solution of the original optimization problem using a standard gradient descent method initialized at these preliminary estimates. Thus, in practice, one can simply combine our computationally efficient estimator with the inferential theory provided in \citet{chen2021nonlinear} to construct confidence intervals and perform hypothesis testing; we also provide an R package for empirical implementation.
\end{abstract}
\thispagestyle{empty}
\clearpage
\setcounter{page}{1}
\clearpage
\section{Introduction}\label{sec:intro}
The importance of accounting for interactive unobserved heterogeneity in panel and network models is well recognized. For example, in linear panel models, interactive fixed effects are routinely controlled for using, for example, the seminal approaches of \citet{bai2009panel} or \citet{pesaran2006estimation}. While analogous methods for nonlinear models have been developed in the literature (e.g., \citealp{chen2021nonlinear}), they see much more limited use in empirical research due to their rapidly growing computational complexity or the lack of inferential theory.\footnote{For example, \citet{zeleneev2020identification} proposes a method for estimating network models with (nonparametric) interactive unobserved heterogeneity that does not require solving a high-dimensional nonconvex problem. However, unlike \citet{chen2021nonlinear}, \citet{zeleneev2020identification} focuses on identification and consistent estimation and does not provide inference tools.}
The main goal of this paper is to bridge the gap between the recent theoretical developments by \citet{chen2021nonlinear} and empirical work by providing a new computationally efficient estimator that can be feasibly implemented in a wide range of nonlinear (semiparametric) settings with unobserved effects following a linear factor structure. We demonstrate that our estimator has two important properties. First, unlike the approach of \citet{chen2021nonlinear}, our method does not require solving a high-dimensional non-convex optimization problem, so our estimator can be efficiently computed for large number of units $N$ and time periods $T$. Second, we argue that our estimator is asymptotically equivalent to the fixed effects (FE) estimator of \citet{chen2021nonlinear}. This means that, in practice, one can combine our computationally efficient estimator with the inferential theory provided in \citet{chen2021nonlinear} to construct confidence intervals for various objects of interest including structural parameters and average partial effects.
Our proposed estimation procedure involves the following two steps. In the first step, we obtain preliminary estimators of the parameters of interest, including the loadings and the factors, by solving a convex relaxation of the original (non-convex) optimization problem in \citet{chen2021nonlinear}. Following the literature, we convexify the original problem by replacing the low-rank constraint imposed on the unobserved effects by the factor model with a nuclear norm penalty. Then, we obtain our final estimator by solving the original optimization problem using a standard gradient descent method initialized at the preliminary nuclear norm regularized (NNR) estimator obtained in the first step.
To demonstrate that our final estimator is asymptotically equivalent to the FE estimator of \citet{chen2021nonlinear} defined as the global solution of the original high-dimensional and non-convex optimization problem, we show that the original problem is locally convex in a \emph{shrinking} neighborhood around the true value of the parameters. Importantly, in the general nonlinear setting studied in this paper (with a growing number of factors and loadings as $N,T \rightarrow \infty$), the size of this neighborhood \emph{shrinks} at a certain rate. To establish the desired result, we characterize the rate of convergence of our preliminary NNR estimator, and demonstrate that this rate is \emph{sufficiently fast} to ensure that our NNR estimator, as well as the FE estimator, falls into that shrinking neighborhood with probability approaching one.
The idea of using a preliminary NNR estimator to initialize local optimization in (globally) non-convex problems has been previously explored in the econometrics literature. For example, \citet{moon2018nuclear} originally proposed an analogous two-step approach for estimating linear panel models with interactive fixed effects. In particular, \citet{moon2018nuclear} also demonstrate that their two-step estimator is asymptotically equivalent to the LS estimator of \citet{bai2009panel}. However, extending these ideas and formally establishing an analogous equivalence result in the general nonlinear setting of \citet{chen2021nonlinear} is a non-trivial task involving additional technical challenges.
As highlighted above, the main conceptual and technical difference is that, in the general nonlinear case, the objective function is locally convex only in a shrinking neighborhood of the true parameters value. In particular, unlike in the linear case, one cannot simply profile out the fixed effects using the singular value decomposition, and demonstrate that the profiled objective function (only depending on the common parameters $\beta$) is locally convex. Since, in the nonlinear case, we cannot work with the profiled objective function directly, we establish local convexity of the original objective function by inspecting its hessian taken with respect to all of the parameters including the loadings and the factors. The analysis is further complicated by the fact that the dimension of the parameter space and hence the dimension of the hessian grows with $N,T \rightarrow \infty$. As a result, local convexity of the objective function can only be established in a shrinking neighborhood of the true parameters. Establishing local convexity in that neighborhood and characterizing at which rate it shrinks is a technical innovation of the paper having important practical implications. Specifically, it imposes an additional requirement on the preliminary estimator's rate of convergence: unless the preliminary estimator falls into that shrinking convexity region with probability approaching one, we cannot guarantee that the second step local optimization finds the global solution. In particular, it turns out that the rate obtained by \citet[Theorem~5]{moon2018nuclear} for the NNR estimator in single-index models is \emph{not} sufficiently fast to satisfy this requirement.
To take advantage of the local convexity result described above, we provide a new improved error bound for the NNR estimator in nonlinear models with interactive fixed effects. Following the literature, we derive this result under a version of the restricted strong convexity (RSC) condition. While various variations of the RSC condition are routinely employed for deriving analogous results in low-rank models (e.g., \citealp{moon2018nuclear,chernozhukov2019inference,ma2022detecting}), these conditions are often difficult to verify. Unlike most previous studies, we provide a set of primitive conditions which can be used to verify the RSC condition in a wide range of panel models allowing, in particular, for predetermined covariates.
To make the proposed approach readily applicable, we supplement it with concrete practical implementation details. In particular, we provide specific optimization algorithms, which can be used to efficiently compute the preliminary NNR and the final estimators, and establish their convergence. We also propose data-driven ways of choosing the regularization parameter involved in the first step and determining the unknown number of factors. To further support empirical use, we release an optimized R package, \texttt{NNRPanel}, which provides an efficient implementation of the proposed method.\footnote{Available at: \url{https://github.com/wszhang-econ/NNRPanel}.} We study the finite sample properties of our estimator in a number of numerical experiments and document its excellent performance and computational efficiency even in large panels with $(N,T) = (1000,200)$. Finally, we revisit the empirical application of \citet{chen2021nonlinear}.
\bigskip
This paper contributes to the literature on estimation of panel (and network) models with interactive fixed effects in two important ways.
First, we complement the theoretical analysis of nonlinear panel models provided by \citet{chen2021nonlinear} by proposing a new estimator that is asymptotically equivalent to their FE estimator. Importantly, unlike their FE estimator, our estimator does not involve solving a high-dimensional non-convex optimization problem, making it an attractive, if not the only available, computationally efficient alternative, which can be feasibly implemented even when both $N$ and $T$ are large. Our two-step approach to solving a non-convex optimization problem essentially extends the proposal of \citet{moon2018nuclear} to nonlinear settings. However, as explained above, establishing the asymptotic equivalence between the two-step and FE estimators in nonlinear settings is more nuanced: it involves careful establishing of local convexity of the criterion function in a shrinking neighborhood of the true parameters value, resulting in additional requirements imposed on the preliminary NNR estimator's convergence rate absent in the linear case studied by \citet{moon2018nuclear}.
Second, we also contribute to the literature on nuclear norm regularized estimation of low-rank models by extending the previously available results established by \citet{moon2018nuclear} and \citet{chernozhukov2019inference} for linear panel models to nonlinear settings.
By verifying the RSC condition in a wide class of nonlinear panel models, we improve on the result of \citet[Theorem~5]{moon2018nuclear}, which extends their original analysis to single-index models. Importantly, our analysis allows for predetermined covariates such as the outcome's lags, which are routinely used in panels, whereas the existing studies providing error bounds for the NNR estimator either only consider strictly exogenous covariates or do not verify the RSC condition at all.
The idea of using nuclear norm regularization to turn estimation of low-rank models, such as factor models, into a convex problem has been extensively applied in various settings in statistics and econometrics. In econometrics, its numerous recent applications include estimation of pure factor models \citep{BaiNg2017,bai2019rank}, estimation of linear \citep{moon2018nuclear,chernozhukov2019inference,beyhum2019square,mugnier2025simple} and quantile panel regressions \citep{belloni2019high,wang2022low,feng_2023}, and treatment effect estimation \citep{athey_matrix_2021,fernandez2021low}. Nuclear norm relaxations have also been proved useful in constructing estimation and inference methods robust to weak factors \citep{armstrong2022robust} and missing data \citep{su2025estimation}. Other recent applications of nuclear norm regularization also include, among others, network recovery and community detection (\citealp{alidaee2020recovering} and \citealp{ma2022detecting}), and estimation of panel threshold models and high-dimensional VARs (\citealp{miao2020panel} and \citealp{miao2023high}).
\bigskip
\noindent {\bf Notation} For any vector $u\in \mathbb{R}^n$, its Euclidean norm is denoted as $\|u\| = \left(u'u\right)^{\frac{1}{2}}$. For any matrix $A\in \mathbb{R}^{m\times n}$, we use $A'$ to denote the transpose of $A$, and use $\|A\|_{\mathrm{F}} = \left(\mathrm{trace}( A'A )\right)^{\frac{1}{2}}$ to denote the Frobenius norm. Furthermore, the singular values of $A$ are arranged in non-increasing order: $\psi_1\left(A\right)\geq \psi_2\left(A\right) \geq \ldots \geq \psi_{\min\{m, n\}}\left(A\right) \geq 0$.
The $\ell^2$ operator norm, $\|A\|_{\mathrm{op}} = \psi_1\left(A\right)$, is the maximum singular value of the matrix, and the nuclear norm is the sum of all singular values:
$\|A\|_{\mathrm{nuc}} = \sum_{i=1}^{\min\{m, n\}}\psi_i\left(A\right)$. We also use $\|A\|_{\max} = \max_{i,j} |A_{ij}| $ to denote the element-wise norm.
When $A$ is a square matrix, we use $\sigma_i(A)$ to denote $A$'s $i$-th largest eigenvalue. We also use $\psi_{\max}$, $\psi_{\min}$, $\sigma_{\max}$, $\sigma_{\min}$ to denote the max/min singular values and max/min eigenvalues respectively. For any matrix $A$, define the coprojection matrix as $M_{A}: = \mathbb{I} - A(A'A)^{\dagger }A'$, where $\mathbb{I}$ denotes the identity matrix of appropriate size and the super-script $\dagger$ denotes the Moore-Penrose generalized inverse. Finally, for any two square matrices $A$ and $B$ of the same dimension, we use $A \geq B$ to denote that $A - B$ is positive semi-definite, and $A > B$ to denote that $A - B$ is positive definite. We use the abbreviation wpa1 instead of with probability approaching to 1.
\bigskip
The remainder of the paper is organized as follows. Section~\ref{sec:model} introduces the model, highlights the computational challenges of the FE estimators, and describes the proposed two-step estimator. Section~\ref{sec:theory} presents the asymptotic equivalence between our two-step estimator and the FE estimator. Section~\ref{sec:implementation} provides practical implementation details, including optimization algorithms, data-driven selection of the tuning parameter, and determination of the number of factors. Section~\ref{sec:MC} reports Monte Carlo simulation results and an empirical application to a gravity model of trade. Section~\ref{sec:conclusion} concludes. Extensions, technical results, and proofs are collected in the Appendix.
\section{The Model and Two Step Estimation}\label{sec:model}
\subsection{The Model}
We observe data $\{(Y_{it}, X_{it})\}_{ 1 \leq i \leq N, 1 \leq t\leq T}$, where $Y_{it}$ is a scalar outcome variable and $X_{it} \in \mathbb R^{d_X}$ a vector of covariates. For concreteness, we adopt the standard panel notation with $i$ indexing units and $t$ indexing time periods, but it should be understood that the considered framework applies to general two-way settings. For example, in a directed network $i$ and $t$ could index senders and receivers (e.g., exporters and importers in an international trade network). The covariates $X_{it}$ could be strictly exogenous or predetermined, e.g., our framework also accommodates lagged outcomes as covariates in panels.
We assume that the (conditional) distribution of $Y_{it}$ belongs to a known family of distributions and is determined by the latent index $Y_{it}^*$, i.e., we assume that the (conditional) log-likelihood takes the form
\begin{align}\label{eq:true_model}
\log f (Y_{it}|X_{it},\lambda_{0,i},\gamma_{0,t}) = \ell (Y_{it}|Y^*_{it}), \quad Y^*_{it} = X'_{it}\beta_0 + \lambda_{0, i}' \gamma_{0, t},
\end{align}
where $\ell(\cdot|Y_{it}^*)$ is a known log-likelihood function, and $\beta_0 \in \mathbb R^{d_X}$ is a parameter of interest. Here, $\lambda_{0,i} \in \mathbb R^R$ and $\gamma_{0, t} \in \mathbb R^R$ are unobserved interactive unit and time effects, commonly referred to as loadings and factors. This formulation is substantially more flexible than the routinely employed two-way fixed effects (TWFE) model, $\lambda_{0, i} + \gamma_{0, t}$, because it allows incorporating multidimensional heterogeneous individual responses $\lambda_{0, i}$ to time-varying aggregate shocks $\gamma_{0, t}$.\footnote{\citet{chen2021nonlinear} also argue that the interactive fixed effects model is sufficiently flexible to allow for homophily based on unobservables (as well as for degree heterogeneity) in network settings.} In particular, the TWFE model corresponds to the special case of the interactive fixed effects model with $R = 2$, where $\lambda_{0, i} = (\lambda_{i1}, 1)'$ and $\gamma_{0, t} = (1, \gamma_{ t1})'$.
Following \citet{chen2021nonlinear}, we consider the so-called large $N,T$ asymptotics with $N,T \rightarrow \infty$ whereas we treat both $d_X$ and $R$ as fixed. For now, we will also assume that the number of factors $R$ is known; we will discuss estimation of $R$ in Section~\ref{sec:implementation}. Finally, we do not put additional restrictions on the relationship between the covariates and the unobserved effects, i.e., we adopt the fixed effects approach.
While the single index formulation~\eqref{eq:true_model} is restrictive, it covers a number of important nonlinear models including binary response models such as Probit and Logit, and Poisson regression.
\begin{example}[Binary response model]
Let $Y_{it}$ be a binary outcome, and consider the binary response model $Y_{it} = \boldsymbol{1}(Y_{it}^* - \epsilon_{it} \geqslant 0)$, where $\boldsymbol{1}(\cdot)$ is the indicator function. Let $F(\cdot)$ be the cumulative distribution function (CDF) of $\epsilon_{it}$, for example, the logistic distribution for the Logit model or the standard normal distribution for the Probit model. The conditional distribution of $Y_{it}$ is given by:
\begin{align*}
\mathbb{P}(Y_{it} = y\mid Y^*_{it}) = F(Y^*_{it})^{y} (1-F(Y^*_{it}))^{(1-y)},\quad y\in \{0, 1\}
\end{align*}
\end{example}
\begin{example}[Poisson model]
Let $Y_{it}$ be count data. Conditional on $Y_{it}^*$, $Y_{it}$ conforms Poisson distribution with parameter $\exp(Y^*_{it})$. The conditional probability of $Y_{it}$ is given by:
\begin{align*}
\mathbb{P}(Y_{it} = y\mid Y^*_{it}) = \frac{\exp (-\exp(Y^*_{it}))(\exp(Y^*_{it}))^{y}}{y!}, \quad y = 0,1,2,\ldots
\end{align*}
\end{example}
\subsection{Fixed Effects MLE Estimator and Computational Challenges}\label{ssec:challenges}
\citet{chen2021nonlinear} propose estimating the model by the fixed effects (FE) MLE estimator maximizing the conditional log-likelihood jointly over the common parameters $\beta$, loadings $\{\lambda_{i}\}_{ 1\leq i\leq N}$ and factors $\{\gamma_{ t}\}_{1\leq t\leq T}$. Specifically, the FE estimator $(\hat{\beta}_{\mathrm{FE}}, \hat{\Lambda}_{\mathrm{FE}}, \hat{\Gamma}_{\mathrm{FE}})$ solves
\begin{equation}\label{eq:FE_estimator}
(\hat{\beta}_{\mathrm{FE}}, \hat{\Lambda}_{\mathrm{FE}}, \hat{\Gamma}_{\mathrm{FE}}) \in \operatorname*{argmin}_{ \beta, \Lambda, \Gamma} \underbrace{-\frac{1}{NT} \sum_{i=1}^{N}\sum_{t=1}^{T} \ell(Y_{it}\mid X_{it}'\beta + \lambda_{ i}'\gamma_{ t})}_{\mathcal{L}_{NT}(\beta, \Lambda, \Gamma)},
\end{equation}
where, for notational simplicity, we collect the unobserved effects $\{\lambda_{i}\}_{ 1\leq i\leq N}$ and $\{\gamma_{ t}\}_{1\leq t\leq T}$ into matrices $\Lambda = \left(\lambda_{1}, \lambda_{2},\ldots, \lambda_{N} \right)'\in \mathbb{R}^{N\times R}$ and $\Gamma = \left(\gamma_{1}, \gamma_{2},\ldots, \gamma_{T} \right)'\in \mathbb{R}^{T\times R}$. Note that problem~\eqref{eq:FE_estimator} does not have a unique solution for $ \hat{\Lambda}_{\mathrm{FE}}$ and $\hat{\Gamma}_{\mathrm{FE}}$ and thus requires a normalization. We will abstract from this issue for now and discuss it in more detail in Section~\ref{sec:theory}.
\cite{chen2021nonlinear} showed that the FE estimator of $\beta_0$ is $\sqrt{NT}$-consistent and asymptotically normal, with an asymptotic incidental parameter bias that can be corrected through various bias reduction methods.
However, despite these well-established theoretical properties, implementing the FE estimator remains a significant computational challenge.
The key computational difficulty is the non-convexity of the objective function $\mathcal{L}_{NT}(\beta, \Lambda, \Gamma)$. To better understand this issue, we reformulate the original optimization problem into an alternative but equivalent form. Let $\theta_{it} = \lambda_i'\gamma_t$ and collect $\theta_{it}$ into a matrix $\Theta\in \mathbb{R}^{N\times T}$. Note that since matrices $\Lambda$ and $\Gamma$ have at most rank $R$, the rank of $\Theta = \Lambda \Gamma'$ is also at most $R$. Equivalently, any matrix $\Theta \in \mathbb R^{N \times T}$ such that $\mathrm{rank}(\Theta)\leq R$ can be represented as $\Lambda \Gamma'$ for some $\Lambda \in \mathbb R^{N \times R}$ and $\Gamma \in \mathbb R^{T \times R}$.\footnote{The representation of $\Theta$ with $\mathrm{rank}(\Theta)\leq R$ as $\Theta = \Lambda \Gamma'$ is not unique. If $\Theta = \Lambda \Gamma'$ for some $\Lambda$ and $\Gamma$, we also have $\Theta = \tilde \Lambda \tilde \Gamma'$ for $\tilde \Lambda = \Lambda G'$ and $\tilde \Gamma = \Gamma G^{-1}$ for any invertible matrix $G \in \mathbb R^{R \times R}$. This non-uniqueness manifests itself in the necessity of normalizing $\Lambda$ and $\Gamma$ in problem~\eqref{eq:FE_estimator} in order to ensure uniqueness of $\hat{\Lambda}_{\mathrm{FE}}$ and $\hat{\Gamma}_{\mathrm{FE}}$.} Thus, problem~\eqref{eq:FE_estimator} can be equivalently reformulated as
\begin{align}\label{eq:FE_rank}
(\hat{\beta}_{\mathrm{FE}}, \hat{\Theta}_{\mathrm{FE}}) \in \operatorname*{argmin}_{\beta\in \mathbb{R}^{d_X}, \Theta \in \mathbb{R}^{N\times T} } \underbrace{ -\frac{1}{NT} \sum_{i=1}^{N}\sum_{t = 1}^{T} \ell(Y_{it}\mid X_{it}'\beta + \theta_{it})}_{\mathcal{L}_{NT}(\beta, \Theta )}, \quad \text{s.t. } \mathrm{rank}(\Theta)\leq R.
\end{align}
The non-convexity arises from the rank constraint $\mathrm{rank}(\Theta)\leq R$: the set of matrices satisfying it is not convex since the sum of two rank-$R$ matrices could have a rank up to $2R$.
The high-dimensional parameter space further exacerbates the computational challenges. When dealing with non-convex optimization problems, it is common practice to start the optimization process with multiple initial values and select the solution that minimizes the objective function. This approach is generally considered effective for finding the global minimum with sufficient trials. However, for problems \eqref{eq:FE_estimator} and \eqref{eq:FE_rank} involving $d_X + R (N+T)$ parameters, this approach becomes intractable even for moderate values of $N$ and $T$.
\begin{remark}
\citet{chen2021nonlinear} propose solving optimization problem~\eqref{eq:FE_estimator} using the EM-algorithm of \citet{chen2016estimation} initialized at multiple initial values. Unfortunately, this method does not overcome the computational challenge discussed above because the EM-algorithm of \citet{chen2016estimation} as well as EM-algorithms in general do not have global convergence guarantees in non-convex problems.
\end{remark}
\subsection{Two-Step Estimation}
To overcome the computational challenges faced by the FE estimator, we propose an alternative two-step estimation procedure. Our procedure does not involve solving a non-convex problem and can be efficiently computed even for large values of $N$ and $T$. Importantly, in Section~\ref{sec:theory}, we demonstrate that, under standard regularity conditions, our two-step estimator is asymptotically equivalent to the FE estimator, whose asymptotic properties have been established in \citet{chen2021nonlinear}. This means that, instead of trying to solve the non-convex and high-dimensional optimization problem~\eqref{eq:FE_estimator} directly, one could compute our two-step estimator and then combine it with the asymptotic theory developed by \citet{chen2021nonlinear} to construct confidence intervals for parameters of interest and for other counterfactuals such as average partial effects (APEs).
\bigskip
Our estimation procedure involves the following two steps.
\bigskip
\noindent {\bf Step 1: Nuclear Norm Regularized (NNR) Estimation}
\noindent The goal of the first step is to construct an easily computable preliminary estimator of $(\beta_0, \Lambda_0, \Gamma_0)$ that is sufficiently close to the global minimizer in \eqref{eq:FE_estimator}. To this end, we consider a convex relaxation of problem~\eqref{eq:FE_rank} of the form
\begin{align}\label{eq:nnr_definition}
\left(\hat{\beta}_{\mathrm{nuc}}, \hat{\Theta}_{\mathrm{nuc}}\right) = \operatorname*{argmin}_{\beta \in \mathbb{R}^{d_X}, \Theta\in \mathbb{R}^{N\times T}} \left\{ \mathcal{L}_{NT} \left( \beta, \Theta\right) + \frac{\varphi_{NT}}{\sqrt{NT}} \|\Theta\|_{\mathrm{nuc}}\right\},
\end{align}
where $\|\Theta\|_{\mathrm{nuc}}$ denotes the nuclear norm of matrix $\Theta$, and $\varphi_{NT} > 0$ is a regularization parameter. We will refer to the solution of this problem $(\hat{\beta}_{\mathrm{nuc}}, \hat{\Theta}_{\mathrm{nuc}})$ as the nuclear norm regularized (NNR) estimator.
Since $\|\Theta\|_{\mathrm{nuc}}$ is a convex function of $\Theta$, problem~\eqref{eq:nnr_definition} is convex when $\mathcal{L}_{NT} \left( \beta, \Theta\right)$ is a convex function of $\beta$ and $\Theta$. This condition is satisfied in important nonlinear models such as Logit, Probit, and Poisson models. Thanks to the convexity of problem~\eqref{eq:nnr_definition}, the NNR estimator can be efficiently computed using, for example, a proximal gradient descent method (e.g., \citealp{hastie2015statistical}) even when the parameter space is high-dimensional. We provide a specific optimization algorithm and a data-dependent recommendation for choosing the regularization parameter $\varphi_{NT}$ in Section~\ref{sec:implementation}.
Notice that problem~\eqref{eq:nnr_definition} can be equivalently rewritten as
\begin{align*}
\left(\hat{\beta}_{\mathrm{nuc}}, \hat{\Theta}_{\mathrm{nuc}}\right) = \operatorname*{argmin}_{\beta \in \mathbb{R}^{d_X}, \Theta\in \mathbb{R}^{N\times T}} \mathcal{L}_{NT} \left( \beta, \Theta\right), \quad \text{s.t. } \|\Theta\|_{\mathrm{nuc}} \leqslant C_{\varphi_{NT}}
\end{align*}
for an appropriately chosen $C_{\varphi_{NT}} > 0$ determined by $\varphi_{NT}$. Thus, problem~\eqref{eq:nnr_definition} can be seen as a convexification of problem~\eqref{eq:FE_rank}, where the non-convex rank constraint is replaced by the slightly looser yet convex constraint $\|\Theta\|_{\mathrm{nuc}} \leqslant C_{\varphi_{NT}}$. Analogously to LASSO using the $\ell_1$-regularization to induce sparsity of the solution in a high-dimensional regression, the nuclear norm regularization (i.e., the $\ell_1$-regularization of the singular values of $\Theta$) induces $\hat \Theta_{\mathrm{nuc}}$ to have low rank (i.e., sparsity of its singular values).
Finally, the nuclear norm regularized estimators $(\hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ are obtained through the singular value decomposition of $\hat{\Theta}_{\mathrm{nuc}}$. Specifically, let
$\hat{\Theta}_{\mathrm{nuc}}/\sqrt{NT} = \hat{U}\hat{D}\hat{V}'$,
where $\hat{U}\in \mathbb{R}^{N\times \min\{N, T\}}$ and $\hat{V}\in \mathbb{R}^{T\times \min\{N, T\}}$ are matrices with left and right orthonormal singular vectors of $\hat{\Theta}_{\mathrm{nuc}}$, and $\hat{D}$ is a diagonal matrix with singular values of $\hat{\Theta}_{\mathrm{nuc}}/\sqrt{NT}$ (arranged in non-increasing order) on its diagonal. Let $\hat{U}_{[:, 1:R]}$ and $\hat{V}_{[:, 1:R]}$ denote the matrices containing the first $R$ columns of $\hat{U}$ and $\hat{V}$, respectively, and $\hat{D}_{[1:R, 1:R]}$ denote the upper-left $R\times R$ diagonal block of $\hat{D}$. We compute $(\hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ as follows:
\begin{gather}\label{eq:space_definition}
\hat{\Lambda}_{\mathrm{nuc}} = \sqrt{N} \hat{U}_{[:, 1:R]}\hat{D}^{1/2}_{[1:R, 1:R]}, \quad \hat{\Gamma}_{\mathrm{nuc}} = \sqrt{T} \hat{V}_{[:, 1:R]}\hat{D}^{1/2}_{[1:R, 1:R]}.
\end{gather}
\bigskip
\noindent {\bf Step 2: Local Estimation}
\noindent While, with appropriately chosen $\varphi_{NT}$, the NNR estimator is consistent for the true values $(\beta_0, \Lambda_0, \Gamma_0)$, it suffers from the regularization bias. To improve on the NNR estimator, in the second step, we solve the original optimization problem~\eqref{eq:FE_estimator} using a standard gradient decent method with $(\hat{\beta}_{\mathrm{nuc}}, \hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ as the initial values. While the original problem~\eqref{eq:FE_estimator} is non-convex, availability of the NNR estimator allows us to guarantee that standard local optimization methods initialized at $(\hat{\beta}_{\mathrm{nuc}}, \hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ converge to the global solution $(\hat{\beta}_{\mathrm{FE}}, \hat{\Lambda}_{\mathrm{FE}}, \hat{\Gamma}_{\mathrm{FE}})$. In particular, in Section~\ref{sec:implementation}, we provide a specific gradient decent algorithm and establish its convergence guarantees.\footnote{In principle, instead of using a gradient decent method, it is possible also employ an EM-algorithm (see, e.g., \citealp{chen2016estimation,chen2021nonlinear}) initialized at $(\hat{\beta}_{\mathrm{nuc}}, \hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$.}
Specifically, to demonstrate that our two-step estimator is (asymptotically) equivalent to the FE estimator, in Section~\ref{sec:theory}, we show that, with probability approaching one, (i) the objective function $\mathcal L_{NT} (\beta, \Lambda, \Gamma)$ is strictly convex in a \emph{shrinking} neighborhood around the true values $(\beta_0, \Lambda_0, \Gamma_0)$ after normalization of the factor, and (ii)~the NNR estimator falls into this neighborhood. The technical difficulty here is that, since the dimension of the parameter space grows with $N,T \rightarrow \infty$, the size of the local neighborhood, in which $\mathcal L(\beta,\Lambda,\Gamma)$ remains convex, shrinks at a certain right. To established the desired result, we characterize (an upper bound for) this rate, and show that the NNR estimator has a sufficiently fast rate of convergence to fall into that neighborhood with probability approaching one.
\bigskip
\noindent {\bf Bias Correction and Inference}
\noindent Since our two-step estimator is asymptotically equivalent to the FE estimator, it also follows the same asymptotic distribution previously derived by \citet{chen2021nonlinear}. In particular, the two-step estimator of $\beta_0$ suffers from the incidental parameter bias caused by estimation of a large number of nuisance parameters. For various bias correction methods as well as for the general asymptotic properties of the FE estimator, we refer the reader to \citet{chen2021nonlinear}.
\section{Asymptotic Analysis}\label{sec:theory}
In this section, we present the consistency of the NNR estimator and establish the local convexity for the original optimization problem~\eqref{eq:FE_estimator}. Combining these results, we demonstrate the asymptotic equivalence between our two-step estimator and the FE estimator.
\subsection{Consistency of NNR Estimator}
We establish the consistency of $(\hat{\beta}_{\mathrm{nuc}}, \hat{\Theta}_{\mathrm{nuc}})$ as in~\eqref{eq:nnr_definition}, as well as the associated nuisance estimators $(\hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ defined in~\eqref{eq:space_definition}, under a set of regularity and high-level conditions. The analysis in the main text focuses on strictly exogenous $X_{it}$ and extends to the more general case in Appendix~\ref{sec:extension}.
For notational simplicity, we collect $X_{it, d}$ into covariate matrices $X_{d} \in \mathbb{R}^{N\times T}$ for each $d= 1,2,\ldots, d_X$, and let $X$ be the collection of all covariate matrices $X= \{X_{1}, \ldots, X_{d_X}\}$. Whenever it does not cause confusion, for any $Y^*$ we abbreviate $\ell_{it}(Y^*) = \log \mathbb{P}(Y_{it}\mid Y^* )$ to denote the log-likelihood evaluated at index $Y^*$. We further denote derivatives of the log-likelihood with respect to the index by $\dot{\ell}_{it}(\cdot), \ddot{\ell}_{it}(\cdot), \ldots$. In addition, we use $\mathbb{P}_{X, \Lambda_0, \Gamma_0} = \mathbb{P}(\cdot\mid X, \Lambda_0, \Gamma_0)$ to denote the conditional probability and $\mathbb{E}_{X, \Lambda_0, \Gamma_0} = \mathbb{E}(\cdot\mid X, \Lambda_0, \Gamma_0)$ to denote the conditional expectation. We now introduce the regularity conditions.
\begin{assumption}[Regularity Conditions]\label{assumption:regularity}
Suppose that
\begin{enumerate}[label=(\roman*)]
\item \label{item:sampling} \textbf{(Sampling)} For each $i = 1,2,\ldots, N$ and $t=1,2,\ldots, T$, conditional on $(X, \Lambda_{0}, \Gamma_{0})$, $\{Y_{it}\}_{1\leq i \leq N, 1\leq t \leq T} $ is distributed as in \eqref{eq:true_model} and
is independent across $i$ and $t$.
\item \label{item:compactness} \textbf{(Compactness)} The parameter spaces for $\beta$, $\lambda_i$, and $\gamma_t$ are uniformly bounded for all $i, t, N, T$. In addition, there exists a constant $\rho_X>0$ such that $\max_{d=1,\ldots, d_X}\|X_d\|_{\max} \leq \rho_X$ for all $i, t, N, T$.
\item \label{item:smoothness} \textbf{(Smoothness and convexity)} $-\ell_{it}(\cdot)$ is four times continuously differentiable and strictly convex almost surely. Furthermore, we assume that $0 < b_{\min} \leq -\ddot{\ell}_{it}( X_{it}' \beta + \lambda_i' \gamma_t) \leq b_{\max} <\infty$ almost surely for all $\beta, \lambda_i, \gamma_t$ in the parameter space uniformly over $i, t, N, T$.
\item \label{item:strong_factors} \textbf{(Strong factors)} Assume that $\frac{1}{N} \sum_{i=1}^{N}\lambda_{0, i}\lambda_{0, i}' \stackrel{p}{\longrightarrow} \Sigma_{\lambda}$ and $\frac{1}{T} \sum_{t=1}^{T}\gamma_{0, t}\gamma_{0, t}'\stackrel{p}{\longrightarrow} \Sigma_{\gamma}$, where $ \Sigma_{\lambda} >0$ and $\Sigma_{\gamma} >0$. In addition, the eigenvalues of $\Sigma_{\lambda}\Sigma_{\gamma}$ are distinct.
\item \label{item:X_generalized_nonlinearity} \textbf{(Generalized non-collinearity)} For any $\Gamma\in \mathbb{R}^{T\times R}$, let $M_{\Lambda_0}$ and $M_{\Gamma}$ be coprojection matrices of $\Lambda_0$ and $\Gamma$ respectively. The $d_{X}\times d_{X}$ matrix $D(\Gamma)$ with elements
\begin{align*}
D(\Gamma)_{d_1, d_2} = \frac{1}{NT} \mathrm{Tr} \left(M_{\Lambda_0}X_{d_1} M_{\Gamma}X_{d_2}'\right), \quad d_{1}, d_{2} = 1, \ldots, d_X
\end{align*}
satisfies $\inf_{\Gamma \in \mathbb{R}^{T\times R}} \sigma_{\min} (D(\Gamma) )>0$, wpa1.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:regularity}\ref{item:sampling} concerns the independence of $Y_{it}$ across $i$ and $t$. This aligns with the sampling assumption in \citet{chen2021nonlinear} and is primarily applicable in contexts where $X_{it}$ is strictly exogenous. Furthermore, it is well suited to contexts involving network data, where the ordering of agents is irrelevant. A more general assumption accommodating predetermined covariates is discussed in Appendix~\ref{sec:extension}, without affecting the main results in the paper.
Assumption~\ref{assumption:regularity}\ref{item:compactness} imposes compactness on the parameter spaces for $\beta$, $\Lambda$, and $\Gamma$, as well as the boundedness on covariates. This assumption is widely adopted in the literature to derive concentration bounds (see, for example, \citealp{chernozhukov2019inference}, \citealp{chernozhukov2023inference}, and \citealp{ma2022detecting}). It is worth noting that \citet{fernandez2016individual} and \citet{chen2021nonlinear} do not impose the boundedness of $\lambda_i$, $\gamma_t$, or $X_{it}$ because their analyses focus on the local properties of the loss function. In contrast, our study requires stronger conditions to ensure global control over parameter behavior. The compactness of $\lambda_i$ and $\gamma_t$ could, in principle, be replaced by assuming that $\{\lambda_i\}_{1\leq i\leq N}$ and $\{\gamma_t\}_{1\leq t\leq T}$ are sub-Gaussian sequences that are independent across $i$ and weak dependence across $t$, respectively.
Assumption~\ref{assumption:regularity}\ref{item:smoothness} is commonly adopted in the nonlinear panel regression literature (see \citealp{fernandez2016individual} and \citealp{chen2021nonlinear}) and is satisfied by Logit, Probit, and Poisson models.
Assumption~\ref{assumption:regularity}\ref{item:strong_factors} refers to the strong factor condition and is standard in the factor model literature.\footnote{Developing estimation and inference methods robust to weak factors is an important but highly nontrivial problem, even in linear panels; see \citet{armstrong2022robust}. In this paper, we simply follow the set-up of \citet{chen2021nonlinear} and leave the important problem of allowing for weak factors in nonlinear models for future research.} Additionally, the compactness of the nuisance parameter space ensures that the maximum eigenvalues of $\Sigma_{\lambda}$ and $\Sigma_{\gamma}$ are bounded. We also impose the assumption that $\Sigma_{\lambda}\Sigma_{\gamma}$ has distinct eigenvalues, but it is not necessary. It is introduced purely to simplify the discussion of technical aspects in the main text. In the Appendix, we demonstrate that relaxing this assumption does not affect our main results.
Assumption~\ref{assumption:regularity}\ref{item:X_generalized_nonlinearity} is a generalized non-collinearity condition that rules out covariates that do not display variation in the individual and time dimensions, such as time-invariant or individual-invariant regressors. This assumption is identical to Assumption 1(vii) in \citet{chen2021nonlinear}, and we refer readers to that paper for further discussion.
\bigskip
We now turn to a key condition for establishing the consistency of our NNR estimator, the restricted strong convexity (RSC) condition. Establishing error bounds for NNR estimators requires analyzing the Hessian of the loss function. As in nonlinear models with a fixed number of parameters, this requires the loss function to have a positive definite Hessian, or equivalently, exhibit strong convexity. However, in the context of problem~\eqref{eq:nnr_definition}, it is not possible for the Hessian matrix to be positive definite, as the number of parameters grows with $N$ and $T$ and exceeds the number of observations. Nonetheless, one can impose an RSC condition that guarantees strong convexity over a restricted subspace in which the parameter space is approximately low-rank. This motivates the adoption of RSC as a key assumption for addressing the high-dimensional nature of our estimator.
To elaborate, for any $(\beta, \Theta)$, the second-order Taylor remainder of the loss function $\mathcal{L}(\beta, \Theta)$ around the true parameter $(\beta_0, \Theta_0)$ satisfies
\begin{equation}\label{eqref:second_order_convexity}
\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}(-\ddot{\ell}_{it}(X_{it}'\tilde{\beta} + \tilde{\theta}_{it})) (X_{it}' \Delta_{\beta} + \Delta_{\theta_{it}} )^2 \geq \frac{b_{\min}}{2} \underbrace{\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T} (X_{it}' \Delta_{\beta} + \Delta_{\theta_{it}} )^2}_{\mathcal{E}_{NT} (\Delta_{\beta}, \Delta_{\Theta}) }
\end{equation}
where $(\tilde{\beta}, \tilde{\Theta})$ lies between $(\beta, \Theta)$ and $(\beta_0, \Theta_0)$, $\Delta_{\beta} = \beta - \beta_0$, $\Delta_{\theta_{it}} = \theta_{it} - \theta_{0,it}$, and the inequality follows from the fact that the second-order derivative $-\ddot{\ell}_{it}(\cdot)$ is uniformly bounded below by $b_{\min}$ (Assumption~\ref{assumption:regularity}\ref{item:smoothness}). Thus, it suffices to study the properties of $\mathcal{E}_{NT}(\Delta_{\beta}, \Delta_{\Theta})$ to ensure that the loss function is strongly convex over a restricted parameter space.
The RSC condition was originally introduced by \citet{negahban2012unified} and has been widely applied in matrix estimation problems with low-rank constraints, including matrix completion (\citealp{negahban2012restricted}), reduced-rank regression (\citealp{rohde2011estimation}), and econometric analyses of low-rank models (\citealp{moon2018nuclear, chernozhukov2019inference}). We are now ready to formalize the RSC condition in our setting.
\begin{assumption}[Restricted strong convexity (RSC)]\label{assumption:RSC}
For any $c_0 >0$, define
\begin{gather*}
\mathcal{C}_1 = \left\{(\Delta_{\beta}, \Delta_{\Theta})\in (\mathbb{R}^{d_X}\times \mathbb{R}^{N\times T})\mid \|M_{\Lambda_0}\Delta_{\Theta}M_{\Gamma_0}\|_{\mathrm{nuc}} \leq c_0 (\sqrt{NT}\|\Delta_{\beta}\| + \|\Delta_{\Theta} - M_{\Lambda_0}\Delta_{\Theta}M_{\Gamma_0}\|_{\mathrm{nuc}})\right\} \\
\mathcal{C}_2 = \left\{ (\Delta_{\beta}, \Delta_{\Theta})\in (\mathbb{R}^{d_X}\times \mathbb{R}^{N\times T})\mid \|\Delta_{\beta}\|^2 + \frac{1}{NT} \|\Delta_{\Theta}\|_{\mathrm{F}}^2 \geq \sqrt{\frac{\log (NT)}{NT}}\right\}
\end{gather*}
We assume that here exist constants $\kappa, \eta >0$, independent of $N, T$, such that for any $(\Delta_{\beta}, \Delta_{\Theta})\in \mathcal{C}_1 \cap \mathcal{C}_2$, the following inequality holds:
\begin{align}\label{eq:RSC}
\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}(X_{it}' \Delta_{\beta} + \Delta_{\theta_{it}} )^2 \geq \kappa \left(\|\Delta_{\beta}\|^2 + \frac{1}{NT}\|\Delta_{\Theta}\|_{\mathrm{F}}^2\right) - \eta \frac{N+T}{NT} (\log(NT))^2, \quad \text{wpa1}
\end{align}
\end{assumption}
The first part of the assumption defines the restricted space of interest $\mathcal{C}_1\cap \mathcal{C}_2$.
$\mathcal{C}_1$ can be viewed as an approximately low-rank space. The term $\Delta_{\Theta} - M_{\Lambda_0}\Delta_{\Theta}M_{\Gamma_0}$ on the right-hand side represents the component that can be explained by $\Lambda_0$ and $\Gamma_0$, serving as a low-rank approximation. In contrast, the left-hand side, $M_{\Lambda_0}\Delta_{\Theta}M_{\Gamma_0}$, corresponds to the residual of $\Delta_{\Theta}$ that cannot be explained by $\Lambda_0$ and $\Gamma_0$, interpreted as the low-rank approximation residual. Therefore, $\mathcal{C}_1$ consists of matrices whose low-rank approximation residuals (in terms of nuclear norm) are small compared to their low-rank approximation (along with the estimation error of $\beta$). The set $\mathcal{C}_2$ is introduced to restrict our attention to scenarios of primary interest. Since we can directly obtain bounds, $\|\Delta_{\beta}\|_{2}^2 \leq \sqrt{\frac{\log (NT)}{NT}}$ and $\frac{1}{NT} \|\Delta_{\Theta}\|_F^2 \leq \sqrt{\frac{\log (NT)}{NT}}$ for matrices that do not belong to this space, focusing on $\mathcal{C}_2$ simplifies the analysis without loss of generality.
The second part introduces an inequality within the restricted space $\mathcal{C}_1 \cap \mathcal{C}_2$, with an additional tolerance term $\eta \frac{N+T}{NT} (\log(NT))^2$ to account for the randomness in $X_{it}$.
Combined with~\eqref{eqref:second_order_convexity}, inequality~\eqref{eq:RSC} ensures that the loss function is strongly convex (up to the tolerance term) over the restricted space $\mathcal{C}_1 \cap \mathcal{C}_2$, even in the high-dimensional setting.
\begin{remark*}
Our RSC conditions are similar to those employed in the low-rank estimation literature. For instance, \citet{moon2018nuclear} employs an RSC condition in the context of linear panel models with interactive fixed effects and homogeneous slopes. Similarly, \citet{chernozhukov2019inference} applies the RSC condition to linear panel regression models with interactive fixed effects and heterogeneous slopes. Additionally, \citet{ma2022detecting} utilizes a comparable condition to detect latent community structures in network formation models.
\end{remark*}
While it is straightforward to see that the RSC condition rules out certain types of regressors, such as time-invariant or individual-specific regressors, verifying RSC directly as a high-level condition can be challenging. \citet{moon2018nuclear} provide important insights into the RSC condition by deriving sufficient conditions for RSC in linear panel settings, although these conditions may not be easily verifiable in practice. \citet{chernozhukov2019inference} provides sufficient conditions for verifying RSC in panel data settings, though their results are mainly applicable to strictly exogenous covariates. We extend the results of \citet{chernozhukov2019inference} to allow for predetermined $X_{it}$ and provide a set of easily verifiable primitive sufficient conditions. We believe that, given the fundamental role of the RSC condition in low-rank estimation, providing accessible verification is essential for broadening the applicability of low-rank estimators and enabling applied researchers to adopt them with greater confidence. The following lemma formalizes the verification of the RSC condition.
\begingroup
\begin{lemma}
Under Assumption~\ref{assumption:regularity_pre} and Assumption~\ref{assumption:conditional_independence_RSC} in the Appendix, the RSC condition (Assumption~\ref{assumption:RSC}) is satisfied.
\end{lemma}
\addtocounter{lemma}{-1}
\endgroup
Since the primary focus of this paper is on developing a tractable estimation method, we defer the discussion of the sufficient conditions for RSC to Lemma~\ref{lemma:sufficient_RSC} in Appendix. This allows us to provide a more detailed examination of these conditions without interrupting the flow of the main argument.
\bigskip
For notational simplicity, we impose sign constraint (e.g., fixing the sign of each factor) and normalization constraint and on the nuisance parameters $(\Lambda_0, \Gamma_0)$ such that $\Lambda_0'\Lambda_0/N$ and $\Gamma_0'\Gamma_0/T$ are diagonal, with $\Lambda_0'\Lambda_0/N = \Gamma_0'\Gamma_0/T$. These constraints are consistent with the construction of $(\hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ in~\eqref{eq:space_definition}. The feasibility of these constraints is discussed in the Appendix~\ref{sec:extension}. The following theorem establishes the convergence rates of NNR estimators.
\begin{theorem}\label{thm:consistency}
For any $\alpha>0$ such that $\varphi_{NT} \geq (1+\alpha) \max\{ \|\nabla_{\beta}\mathcal{L}_{NT}\left(\beta_0 , \Theta_0\right)\|_{2}, \sqrt{NT}\|\nabla_{\Theta}\mathcal{L}_{NT}\left(\beta_0 , \Theta_0\right)\|_{\mathrm{op}} \}$, under Assumption~\ref{assumption:regularity}, Assumption~\ref{assumption:RSC}, as $N, T\rightarrow\infty$, there exist constants $c_1, c_2>0$ that do not depend on $N, T$ such that wpa1:
\begin{align*}
\| \hat{\beta}_{\mathrm{nuc}} - \beta_0\| & \leq c_1 \left(\varphi_{NT} + \log (NT)/\sqrt{\min\{N, T\}}\right) \\
\frac{1}{\sqrt{NT}}\|\hat{\Theta}_{\mathrm{nuc}} - \Theta_0\|_{\mathrm{F}} & \leq c_1 \left(\varphi_{NT} + \log (NT)/\sqrt{\min\{N, T\}}\right)
\end{align*}
In addition, wpa1,
\begin{align*}
\frac{1}{\sqrt{N}}\|\hat{\Lambda}_{\mathrm{nuc}} - \Lambda_0 \|_{\mathrm{F}} & \leq c_2 \left(\varphi_{NT} + \log (NT)/\sqrt{\min\{N, T\}} \right) \\
\frac{1}{\sqrt{T}}\|\hat{\Gamma}_{\mathrm{nuc}} - \Gamma_0 \|_{\mathrm{F}} & \leq c_2 \left(\varphi_{NT} + \log (NT)/\sqrt{\min\{N, T\}} \right)
\end{align*}
\end{theorem}
Theorem~\ref{thm:consistency} establishes that, for a sufficiently large tuning parameter $\varphi_{NT}$, the NNR estimator $\hat{\beta}_{\mathrm{nuc}}$ converges to $\beta_0$ at a rate of order $\varphi_{NT} + \log(NT)/\sqrt{\min\{N,T\}}$. This convergence rate coincides with that obtained by \citet{moon2018nuclear} in the linear panel model. The estimator $\hat{\Theta}_{\mathrm{nuc}}$ achieves the same convergence rate under the normalized Frobenius norm, and the associated estimators of latent factors $(\hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ inherit this rate as well. The convergence rate in Theorem~\ref{thm:consistency} depends on the choice of the tuning parameter $\varphi_{NT}$. While $\varphi_{NT}$ must be sufficiently large to guarantee consistency, choosing it overly large is undesirable. In particular, the theorem suggests that an optimal choice of $\varphi_{NT}$ would be of order $\log(NT)/\sqrt{\min\{N,T\}}$. The following corollary formalizes this choice.
\begin{corollary}\label{corollary:consistency}
Under the conditions of Theorem~\ref{thm:consistency}, let $\varphi_{NT} = O\left(\log (NT)/ \sqrt{ \min\{N, T\}}\right)$. Then there exist constants $c_3, c_4>0$ that do not depend on $N, T$ such that wpa1,
\begin{align*}
\| \hat{\beta}_{\mathrm{nuc}} - \beta_0\| & \leq c_3 \log (NT)/\sqrt{\min\{N, T\}} \\
\frac{1}{\sqrt{NT}}\|\hat{\Theta}_{\mathrm{nuc}} - \Theta_0\|_{\mathrm{F}} & \leq c_3 \log (NT)/\sqrt{\min\{N, T\}}
\end{align*}
In addition, wpa1,
\begin{align*}
\frac{1}{\sqrt{N}}\|\hat{\Lambda}_{\mathrm{nuc}} - \Lambda_0 \|_{\mathrm{F}} & \leq c_4 \log (NT)/\sqrt{\min\{N, T\}} \\
\frac{1}{\sqrt{T}}\|\hat{\Gamma}_{\mathrm{nuc}} - \Gamma_0 \|_{\mathrm{F}} & \leq c_4 \log (NT)/\sqrt{\min\{N, T\}}
\end{align*}
\end{corollary}
Corollary~\ref{corollary:consistency} extends previous results for linear panel models to nonlinear settings. The convergence rate of the NNR estimators, $\log(NT)/\sqrt{\min\{N,T\}}$, is comparable to those established in \citet{moon2018nuclear}, \citet{chernozhukov2019inference}, and \citet{su2025estimation} for linear panel models, up to an additional logarithmic factor that provides an upper bound on the estimation error. In nonlinear models, our result coincides with that of \citet{ma2022detecting} in network settings and improves upon \citet[Theorem~5]{moon2018nuclear}, which extends their original analysis to single-index models. This improvement in the convergence rate is crucial for establishing local convexity and will be discussed further below.
\subsection{Local Convexity}
In this subsection, we establish the asymptotic equivalence between our two-step estimator and the FE estimator. This result builds on the consistency of the NNR estimator and the local convexity of the objective function $\mathcal{L}_{NT}(\beta, \Lambda, \Gamma)$.
Since the NNR estimator $(\hat{\beta}_{\mathrm{nuc}}, \hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ is used to initialize the optimization~\eqref{eq:FE_estimator}, its consistency ensures that the initial estimates lie within a shrinking neighborhood of the true parameters as $N, T \rightarrow \infty$. Thus, it suffices to focus on the local properties of $\mathcal{L}_{NT}(\beta, \Lambda, \Gamma)$ within these shrinking neighborhoods. To formalize the argument, let $\{\delta_{NT}\}$ be a sequence of shrinking radii such that $\delta_{NT}\rightarrow 0$ as $N, T\rightarrow \infty$,
and define the shrinking neighborhood around the true parameters $(\beta_0, \Lambda_0, \Gamma_0)$ as follows:
\begin{align}\label{eq:neighborhood}
\mathcal{B}_{\delta_{NT}} = \bigg\{(\beta, \Lambda, \Gamma) \mid & \|\beta-\beta_0\|, \frac{1}{\sqrt{N}}\|\Lambda - {\Lambda}_0 \|_{\mathrm{F}} , \frac{1}{\sqrt{T}}\|\Gamma - {\Gamma}_0 \|_{\mathrm{F}} \leq \delta_{NT} \bigg\}
\end{align}
The neighborhood $\mathcal{B}_{\delta_{NT}}$ consists of parameters whose distances to the true values $(\beta_0, \Lambda_0, \Gamma_0)$ are less than $\delta_{NT}$.
The second-step estimator (or the local estimator), $(\hat{\beta}_{\mathrm{local}}, \hat{\Lambda}_{\mathrm{local}}, \hat{\Gamma}_{\mathrm{local}}) $, solves:
\begin{equation}\label{eq:definition_local}
\begin{aligned}
(\hat{\beta}_{\mathrm{local}}, \hat{\Lambda}_{\mathrm{local}}, \hat{\Gamma}_{\mathrm{local}}) \in \operatorname*{argmin}_{(\beta, \Lambda, \Gamma) \in \mathcal{B}_{\delta_{NT}} } \mathcal{L}_{NT}(\beta, \Lambda, \Gamma)
\end{aligned}
\end{equation}
where the parameter space is restricted to a shrinking neighborhood $\mathcal{B}_{\delta_{NT}}$.
We establish asymptotic equivalence result provided that, with a properly chosen $\delta_{NT}$, the following conditions hold: (1) the NNR estimator falls within the shrinking neighborhood, i.e., $(\hat{\beta}_{\mathrm{nuc}}, \hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})\in \mathcal{B}_{\delta_{NT}}$ wpa1, (2) the FE estimator, as the global minimizer of problem~\eqref{eq:FE_estimator}, also lies within the shrinking neighborhood up to rotation wpa1, and (3) the objective function $\mathcal{L}_{NT}(\beta, \Lambda, \Gamma)$ is strictly convex within the neighborhood $\mathcal{B}_{\delta_{NT}}$. Under these conditions, we can theoretically guarantee that wpa1: (1) using NNR estimator as the initial value leads to the global minimum, and (2) our estimator is asymptotically equivalent to the FE estimator.
Let $\delta_{NT} = \log(NT) \min\{N^{-3/8}, T^{-3/8}\}$. The NNR estimator falls within the shrinking neighborhood wpa1 by Corollary~\ref{corollary:consistency}. The FE estimator also lies within the shrinking neighborhood wpa1, as shown in Lemma 1 in \citet{chen2021nonlinear}. Therefore, the key step in establishing asymptotic equivalence is to establish local convexity of the original objective function $\mathcal{L}_{NT}(\beta, \Lambda, \Gamma)$ within $\mathcal{B}_{\delta_{NT}}$.
It should be noted that for linear panel models, such as those studied in \citet{moon2018nuclear} and \citet{su2025estimation}, the original objective function is locally convex in $\mathcal{B}_{\delta_{NT}}$ whenever $\delta_{NT}=o_P(1)$. In contrast, local convexity is more delicate in our nonlinear setting, which requires the neighborhood to shrink at a faster rate than in the linear cases.
Analogous to the linear case, we need $R^2$ constraints to identify $(\Lambda, \Gamma)$. However, for any $R$-dimensional non-singular matrix $G$, the conditional distribution of $Y_{it}$ remains invariant under the transformations $\lambda_i' \mapsto G \lambda_i'$ and $\gamma_t' \mapsto G^{-1\prime} \gamma_t'$. This invariance allows us to freely choose different normalization methods for different purposes without affecting the inference of $\beta_0$. Although applied researchers can freely choose different normalization methods, certain choices can greatly facilitate theoretical analysis. In the following text, we follow the idea in \citet{chen2021nonlinear} and focus on the restricted parameter set:
\begin{align*}
\Phi_{NT} = \left\{(\Lambda, \Gamma)\mid \hat{\Lambda}_{\mathrm{nuc}}' \Lambda/N = \Gamma' \hat{\Gamma}_{\mathrm{nuc}} /T \right\}
\end{align*}
which consists of $R^2$ constraints and is consistent with the construction in~\eqref{eq:space_definition}.
Unlike the commonly used approach in the literature, where $\Gamma'\Gamma/T = \mathbb{I}$ and $\Lambda'\Lambda/N$ is diagonal, our restricted set consists of linear functions of $(\Lambda, \Gamma)$, which will bring huge convenience in the following theoretical analysis. Particularly, instead of imposing $\Phi_{NT}$ directly, we convert it into a quadratic penalization term, $\| \hat{\Lambda}_{\mathrm{nuc}}' \Lambda/N - \Gamma' \hat{\Gamma}_{\mathrm{nuc}} /T \|_{\mathrm{F}}^2$,
whose Hessian matrix only depends on $(\hat{\Lambda}_{\mathrm{nuc}}, \hat{\Gamma}_{\mathrm{nuc}})$ due to its linear nature. Define the penalized sample log-likelihood:
\begin{equation}\label{eq:local_estimators}
\begin{aligned}
(\hat{\beta}_{\mathrm{local}}, \hat{\Lambda}_{\mathrm{local}}, \hat{\Gamma}_{\mathrm{local}}) = \operatorname*{argmin}_{(\beta, \Lambda, \Gamma) \in \mathcal{B}_{\delta_{NT}} } \left\{\mathcal{L}_{NT}(\beta, \Lambda, \Gamma) + \frac{1}{2} \| \hat{\Lambda}_{\mathrm{nuc}}' \Lambda/N - \Gamma' \hat{\Gamma}_{\mathrm{nuc}} /T \|_{\mathrm{F}}^2 \right\}
\end{aligned}
\end{equation}
which is equivalent to solving~\eqref{eq:definition_local} with normalization constraint $(\Lambda, \Gamma)\in\Phi_{NT} $. The differentiability of the penalty simplifies the analysis of Hessian, providing a more tractable alternative to the hard constraint.
It is worth noting that the penalization term is introduced for theoretical considerations but is unnecessary in practice.
The sample Hessian of the objective function in \eqref{eq:local_estimators} with respect to $(\beta, \Lambda, \Gamma)$ is given by
\begin{align*}
\mathcal{H}_{NT}(\beta, \Lambda, \Gamma) := \nabla^2 \mathcal{L}_{NT}(\beta, \Lambda, \Gamma) + \frac{1}{2} \nabla^2 \| \hat{\Lambda}_{\mathrm{nuc}}' \Lambda/N - \Gamma' \hat{\Gamma}_{\mathrm{nuc}} /T \|_{\mathrm{F}}^2
\end{align*}
It is a $(d_X + R(N+T))$-dimensional square matrix.\footnote{
The sample Hessian is a matrix-valued function of a $d_X$ -dimensional vector $\beta$, an $N \times R$ parameter matrix $\Lambda$, and a $T \times R$ parameter matrix $\Gamma$. These parameters are arranged as follows:
$$(\beta', \text{vec}(\Lambda')', \text{vec}(\Gamma')')'$$ where $\text{vec}(\cdot)$ denotes the vectorization operator, stacking the columns of a matrix into a vector.
} Establishing local convexity is equivalent to showing that the Hessian $\mathcal{H}_{NT}(\beta, \Lambda, \Gamma)$ is positive definite uniformly over $(\beta, \Lambda, \Gamma)\in \mathcal{B}_{\delta_{NT}}$.
Consider the decomposition of the Hessian as follows:
\begin{align*}
\mathcal{H}_{NT}(\beta, \Lambda, \Gamma) = \underbrace{\mathbb{E}_{X, \Lambda_0, \Gamma_0}\mathcal{H}_{NT}(\beta_0, \Lambda_0, \Gamma_0)}_{\text{population Hessian at true parameters}} + \underbrace{\mathcal{H}_{NT}(\beta, \Lambda, \Gamma)- \mathbb{E}_{X, \Lambda_0, \Gamma_0}\mathcal{H}_{NT}(\beta_0, \Lambda_0, \Gamma_0)}_{\text{deviation}}
\end{align*}
The first term represents the population Hessian evaluated at the true parameters $(\beta_0, \Lambda_0, \Gamma_0)$, whose smallest eigenvalue is strictly positive under general conditions. The second term reflects deviations arising from sampling error and the perturbations of $(\beta, \Lambda, \Gamma)$ from the true parameters. By Weyl's theorem, we can establish the local convexity of~\eqref{eq:local_estimators} if the deviation is negligible (in terms of operator norm) compared with $\mathbb{E}_{X, \Lambda_0, \Gamma_0}\mathcal{H}_{NT}(\beta_0, \Lambda_0, \Gamma_0)$.
The positive definiteness of $\mathbb{E}_{X, \Lambda_0, \Gamma_0}\mathcal{H}_{NT}(\beta_0, \Lambda_0, \Gamma_0)$ can be formalized through the following diagonal block structure assumption:
\begin{assumption}[Diagonal structure]\label{assumption:block}
The population Hessian in~\eqref{eq:local_estimators} evaluated at the true parameters, $\mathbb{E}_{X, \Lambda_0, \Gamma_0}\mathcal{H}_{NT}(\beta_0, \Lambda_0, \Gamma_0)$, admits a diagonal block structure, i.e., there exist a constant $C>0$ (does not depend on $N, T$) such that
\begin{align*}
\mathbb{E}_{X, \Lambda_0, \Gamma_0}\mathcal{H}_{NT}(\beta_0, \Lambda_0, \Gamma_0) \geq C \mathrm{diag}\left\{\mathbb{I}_{d_X}, \frac{1}{N}\mathbb{I}_{N R}, \frac{1}{T}\mathbb{I}_{TR}\right\}
\end{align*}
\end{assumption}
This assumption, which is closely related to the asymptotic diagonal structure condition in \citet{chen2021nonlinear}, \citet{wang2022maximum}, and \citet{su2025estimation}, is mild as it imposes conditions only on the population Hessian evaluated at the true parameters. Easily verifiable sufficient conditions for Assumption~\ref{assumption:block} are provided in Lemma~\ref{lemma:sufficient_convexity} in the Appendix. We now state the main theorem, which establishes local convexity and asymptotic equivalence:
\begin{theorem}\label{thm:convexity_strong}
Under Assumption~\ref{assumption:block} and the conditions in Corollary~\ref{corollary:consistency}, suppose $N, T$ have the same order. Let $ \delta_{NT}= \log (NT) \min\{N^{-3/8}, T^{-3/8}\} $, and $\mathcal{B}_{\delta_{NT}}$ be the neighborhood defined in ~\eqref{eq:neighborhood}. The following results hold:
\begin{enumerate}[label=(\roman*)]
\item \label{item:1} the local optimization problem~\eqref{eq:local_estimators} is strictly convex wpa1;
\item \label{item:2} our two-step estimator is asymptotically equivalent to the FE estimator;
\item \label{item:3} the local optimization problem~\eqref{eq:local_estimators} is strongly convex wpa1 uniformly over $\mathcal{B}_{\delta_{NT}}$, i.e., there exists a constant $c_5>0$ independent of $N, T$ such that for any $(\beta, \Lambda, \Gamma) \in \mathcal{B}_{\delta_{NT}}$,
\begin{align*}
\mathcal{H}_{NT}(\beta, \Lambda, \Gamma) > c_5 \mathrm{diag}\left\{\mathbb{I}_{d_X}, \frac{1}{N}\mathbb{I}_{NR}, \frac{1}{T}\mathbb{I}_{TR}\right\}, \quad \text{wpa1}
\end{align*}
\end{enumerate}
\end{theorem}
Theorem~\ref{thm:convexity_strong} demonstrates that initializing the optimization problem with the NNR estimator ensures convergence to the global minimum, i.e., the FE estimator. The asymptotic equivalence between our two-step estimator and the FE estimator implies that we can leverage the computational advantages of our method while directly applying the theoretical tools from \citet{chen2021nonlinear} for inference and constructing confidence intervals.
Theorem~\ref{thm:convexity_strong} extends the post-estimations in \citet{moon2018nuclear} and \citet{su2025estimation} to nonlinear settings. In linear panel models, consistency of the NNR estimator alone is sufficient to guarantee convergence of the second-step optimization to the global minimum. In contrast, in nonlinear settings, establishing local convexity is more delicate and can only be shown within a shrinking neighborhood of the true parameters under our framework. Therefore, unless the preliminary estimator enters the shrinking convexity region wpa1, convergence of the second-step local optimization to the global solution cannot be guaranteed. In particular, we show that the convergence rate obtained by \citet[Theorem~5]{moon2018nuclear} for the NNR estimator in single-index models is not sufficiently fast to satisfy this requirement.
Finally, Theorem~\ref{thm:convexity_strong}\ref{item:3} further establishes strong local convexity of~\eqref{eq:local_estimators}, implying that simple gradient descent can be effectively applied in the second step to find the global minimum, even in high-dimensional settings.
\section{Implementation}\label{sec:implementation}
In this section, we provide the practical steps for implementing our two-step method, including the algorithm and its theoretical guarantees, as well as a data-dependent approach for selecting the tuning parameter and determining the number of factors.
\subsection{Algorithm: NNR Estimator}
We compute the NNR estimator defined in~\eqref{eq:nnr_definition} using proximal gradient descent method following \citet{hastie2015statistical}. Given the k-step estimates $(\beta^{(k)}, \Theta^{(k)})$, the $k+1$-step estimates are updated by solving
\begin{align*}
\beta^{(k+1)}, \Theta^{(k+1)} \in \arg\min_{\beta, \Theta} \Big\{& \mathcal{L}_{NT}(\beta^{(k)}, \Theta^{(k)}) + \langle \nabla_{\beta}\mathcal{L}_{NT}(\beta^{(k)}, \Theta^{(k)}), \beta - \beta^{(k)}\rangle + \langle \nabla_{\Theta}\mathcal{\mathcal{L}}_{NT}(\beta^{(k)}, \Theta^{(k)}), \Theta - \Theta^{(k)} \rangle \\
& + \frac{1}{2s_{\beta}}\|\beta - \beta^{(k)}\|^2 + \frac{1}{2s_{\theta}} \|\Theta - \Theta^{(k)}\|_{\mathrm{F}}^2 + \frac{\varphi_{NT}}{\sqrt{NT}}\|\Theta\|_{\mathrm{nuc}}\Big\}
\end{align*}
where $\nabla_{\beta}\mathcal{L}_{NT}(\cdot, \cdot)$ is the $d_X$-dimensional vector of gradients with respect to $\beta$, $\nabla_{\Theta}\mathcal{L}_{NT}(\cdot, \cdot)\in \mathbb{R}^{N\times T}$ is a matrix of gradients with respect to $\theta_{it}$, $\langle \cdot, \cdot\rangle$ denotes inner product between two vectors or two matrices, and $s_{\beta}, s_{\theta}>0$ are step sizes. Before introducing the algorithm, it is convenient to define the soft-thresholding operator $\mathcal{S}^*_{s_{\theta}\frac{\varphi_{NT}}{\sqrt{NT}}}: \mathbb{R}^{N\times T}\mapsto \mathbb{R}^{N\times T}$ applied to the singular values of an $N \times T$ matrix with threshold value $s_{\theta}\frac{\varphi_{NT}}{\sqrt{NT}}$. Specifically, for any matrix $A\in \mathbb{R}^{N\times T}$ with singular value decomposition $A = U\Sigma V'$, the operator gives
$$\mathcal{S}^*_{s_{\theta}\frac{\varphi_{NT}}{\sqrt{NT}}}(A) = U\mathrm{diag}\left\{\max\left\{\Sigma_{rr} - s_{\theta}\frac{\varphi_{NT}}{\sqrt{NT}}, 0\right\}_{r=1,\ldots, \min\{N, T\}}\right\}V'$$
We now provide the algorithm obtaining the NNR estimators in \eqref{eq:nnr_definition} as follows:
\begin{algorithm}[Proximal gradient descent]\label{algorithm:NNR}
Compute the nuclear norm regularized estimator as follows:
\begin{itemize}[label={}, leftmargin=2cm]
\item[Step 1:] Fix the step sizes $(s_{\beta}, s_{\theta})$. Initialize $\beta^{(0)}$ and $\Theta^{(0)}$. Set $k = 0$.
\item[Step 2:] Let
\begin{equation}\label{eq:algorithm}
\begin{aligned}
\beta^{(k+1)} & = \beta^{(k)} - s_{\beta} \nabla_{\beta}\mathcal{L}_{NT}\left(\beta^{(k)}, \Theta^{(k)}\right) \\
\Theta^{(k+1)} & = \mathcal{S}^*_{s_{\theta}\frac{\varphi_{NT}}{\sqrt{NT}}}\left(\Theta^{(k)} - s_{\theta}\nabla_{\Theta}\mathcal{L}_{NT}\left(\beta^{(k)}, \Theta^{(k)}\right)\right)
\end{aligned}
\end{equation}
and set $k = k+1$.
\item[Step 3:] Repeat Step 2 until convergence.
\end{itemize}
\end{algorithm}
We establish the convergence of Algorithm~\ref{algorithm:NNR} using a proof strategy similar to that of \citet{nesterov2013gradient}:
\begin{theorem}\label{thm:algorithm_convergence}
Under the conditions of Theorem \ref{thm:consistency}, Algorithm \ref{algorithm:NNR} is guaranteed to converge to the global minimizer if $0 < s_{\beta} < \frac{1}{L_{\beta}}$ and $0 < \frac{s_{\theta}}{NT} < \frac{1}{L_{\theta} }$, where
$ L_{\beta} = 2 d_X b_{\max} \rho_X^2$, $L_{\theta} = 2b_{\max}$, and $\rho_X = \max_{1\leq d\leq d_X} \|X_d\|_{\max}$.
\end{theorem}
The theorem establishes the algorithm's convergence to the global minimizer under sufficiently small $(s_{\beta}, s_{\theta}/(NT))$. The step sizes are on different scales: $s_{\beta}$ and $s_{\theta}/(NT)$ should be of the same order, reflecting their distinct influence on the objective function. This difference arises because a change in $\beta$ affects $\ell_{it}$ for all $(i,t)$, whereas a change in $\theta_{it}$ affects only the corresponding $\ell_{it}$.
Since the optimization problem is convex, the algorithm converges to a global minimizer regardless of the choice of initial values.
Researchers may therefore simply initialize with $\beta^{(0)} = 0$ and $\Theta^{(0)} = 0$. Regarding the choice of step sizes, Theorem~\ref{thm:algorithm_convergence} provides theoretical guarantees but offers limited practical guidance, since $b_{\max}$ is typically unknown. In practice, we recommend starting each iteration with step sizes $s_{\beta} = 1$ and $s_{\theta} = NT$. If the objective function increases, the step sizes are iteratively halved until a decrease in the objective function is achieved.
Solving the nuclear norm-regularized regression is usually more computationally demanding than the second step. The computational bottleneck is computing the singular value decomposition on an $N \times T$ matrix at each iteration. Our R package implementation is computationally efficient for matrices with size $N, T = 1000$. For larger $N$ or $T$, one may apply accelerated proximal gradient method (\citealp{nesterov2013gradient}). Detailed discussion of such acceleration methods is beyond the scope of this paper.
\subsection{Algorithm: local estimator}
In the second step, we employ gradient descent to search for a minimizer of~\eqref{eq:local_estimators}. The following algorithm practically minimizes $\mathcal{L}_{NT}(\beta, \Lambda, \Gamma)$ using the NNR estimator as the initial value. Step 3, which addresses rotational invariance, is optional, since the penalization introduced in~\eqref{eq:local_estimators} is used only for theoretical purposes.
\begin{algorithm}[Gradient descent]\label{algorithm:local}
Compute the local estimator (or the FE estimator) as follows:
\begin{itemize}[label={}, leftmargin=2cm]
\item[Step 1:] Fix the step sizes $(s_{\beta}, s_{\lambda}, s_{\gamma})$. Initialize $\beta^{(0)} = \hat{\beta}_{\mathrm{nuc}}$, $\Lambda^{(0)} = \hat{\Lambda}_{\mathrm{nuc}}$, and $\Gamma^{(0)} = \hat{\Gamma}_{\mathrm{nuc}}$. Set $k = 0$.
\item[Step 2:] Let
\begin{equation}\label{eq:algorithm_local}
\begin{aligned}
\beta^{(k+1)} & = \beta^{(k)} - s_{\beta} \nabla_{\beta}\mathcal{\mathcal{L}}_{NT}(\beta^{(k)}, \Lambda^{(k)}, \Gamma^{(k)} ) \\
\Lambda^{(k+1)} & = \Lambda^{(k)} - s_{\lambda}\nabla_{\lambda}\mathcal{L}_{NT}(\beta^{(k)}, \Lambda^{(k)}, \Gamma^{(k)} ) \\
\Gamma^{(k+1)} & = \Gamma^{(k)} - s_{\gamma}\nabla_{\gamma}\mathcal{L}_{NT}(\beta^{(k)}, \Lambda^{(k)}, \Gamma^{(k)} )
\end{aligned}
\end{equation}
and set $k = k+1$.
\item[Step 3:] (Optional) Normalize $\Lambda^{(k+1)}$ and $\Gamma^{(k+1)}$, for example, let $\frac{1}{N}\Lambda^{(k+1)\prime}\Lambda^{(k+1)} = \frac{1}{T}\Gamma^{(k+1)\prime}\Gamma^{(k+1)}$ and diagonal.
\item[Step 4:] Repeat Step 2 and Step 3 until convergence.
\end{itemize}
\end{algorithm}
The convergence of the gradient descent method can be guaranteed by the following Theorem:
\begin{theorem}\label{thm:algorithm_convergence_local}
Under the conditions of Theorem~\ref{thm:convexity_strong}, Algorithm~\ref{algorithm:local} is guaranteed to converge to the global minimizer when $0 < s_{\beta} < \frac{1}{L_{\beta}}$, $0 < \frac{s_{\lambda}}{N} < \frac{1}{L_{\lambda} }$, and $0 < \frac{s_{\gamma}}{T} < \frac{1}{L_{\gamma} }$, where
$ L_{\beta}, L_{\lambda}, L_{\gamma}$ are sufficiently large constants independent of $N, T$.
\end{theorem}
Similar to the result in Theorem~\ref{thm:algorithm_convergence}, the step sizes have different orders---specifically, $s_{\beta}\sim s_{\lambda}/N \sim s_{\gamma}/T $---reflecting their respective influence on the objective function. In practice, in each step, we recommend starting with $s_{\beta} = 1$, $s_{\lambda} = N$, and $s_{\gamma} = T$. If the objective function increases, we iteratively halve the step sizes $(s_{\beta}, s_{\lambda}, s_{\gamma})$ until the objective function decreases.
\begin{remark}
\citet{chen2021nonlinear} propose using an EM algorithm, a Newton-Raphson-type method that theoretically achieves faster convergence through second-order accuracy. However, we recommend using a gradient descent method instead. This is because in nonlinear panel models, computing and inverting a high-dimensional Hessian matrix can be not only computationally expensive but also numerically unstable. Therefore, we adopt a more robust, albeit slower, gradient-based algorithm.
\end{remark}
\subsection{Determining the number of factors $R$ and tuning parameter $\varphi_{NT}$}
We propose the following data-dependent approach for selecting the tuning parameter $\varphi_{NT}$ and determining the number of factors $R$. Recall that Theorem~\ref{thm:consistency} requires that $\varphi_{NT} > \sqrt{NT}\|\nabla_{\Theta}\mathcal{L}_{NT}(\beta_0, \Theta_0)\|_{\mathrm{op}} $ in order to achieve the desired consistency result.\footnote{
$\|\nabla_{\beta}\mathcal{L}_{NT}(\beta_0, \Theta_0)\|_2$ is negligible compared to $\sqrt{NT}\|\nabla_{\Theta}\mathcal{L}_{NT}(\beta_0, \Theta_0)\|_{\mathrm{op}}$ as $N,T\rightarrow \infty$.
} At the same time, Theorem~\ref{thm:consistency} also indicates that selecting an excessively large $\varphi_{NT}$ should be avoided, as it may induce substantial estimation error in the NNR estimator.
Hence, a preferable choice for $\varphi_{NT}$ is one that slightly exceeds $\sqrt{NT}\|\nabla_{\Theta}\mathcal{L}_{NT}(\beta_0, \Theta_0)\|_{\mathrm{op}}$. Since $\beta_0$ and $\Theta_0$ are unknown, we propose a two-step procedure ---an initial guess followed by an update---following the idea of \citet{chernozhukov2019inference} to determine the tuning parameter, together with a procedure to select the number of factors.
\begin{algorithm}[Determining $\varphi_{NT}$ and $R$]\label{algorithm:tuning}
Compute the local estimator (or the FE estimator) as follows:
\begin{itemize}[label={}, leftmargin=2cm]
\item[Step 1:] Solve the nonlinear panel regression with additive fixed effects to obtain $(\tilde{\beta}_{1}, \tilde{\Lambda}_{1}, \tilde{\Gamma}_{1})$:
\begin{align*}
(\tilde{\beta}_{1}, \tilde{\Lambda}_{1}, \tilde{\Gamma}_{1}) \in \operatorname*{argmin}_{\beta\in \mathbb{R}^{d_X}, \Lambda\in \mathbb{R}^{N}, \Gamma\in\mathbb{R}^{T}} -\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T} \ell(Y_{it}\mid X_{it}'\beta + \lambda_i + \gamma_t)
\end{align*}
We then compute the initial guess for the tuning parameter as
$$\tilde{\varphi}_{NT}:= (1 + \alpha)\|\nabla_{\Theta}\mathcal{L}_{NT}(\tilde{\beta}_1, \tilde{\Lambda}_1, \tilde{\Gamma}_1)\|_{\mathrm{op}}$$
where $\alpha$ is a small positive constant (e.g., $\alpha = 0.05$).
\item[Step 2:] Solve the nuclear-norm regularized optimization problem~\eqref{eq:nnr_definition} with tuning parameter $\tilde{\varphi}_{NT}$ to get $(\tilde{\beta}_2, \tilde{\Theta}_2)$ and then compute $(\tilde{\Lambda}_2, \tilde{\Gamma}_2)$ through~\eqref{eq:space_definition}. The updated tuning parameter used in the first step estimation is
\begin{align*}
\varphi_{NT}:= \left(1 +\alpha\right)\|\nabla_{\Theta}\mathcal{L}_{NT}(\tilde{\beta}_2, \tilde{\Lambda}_2, \tilde{\Gamma}_2)\|_{\mathrm{op}}
\end{align*}
\item[Step 3:] To determine the number of factors, consider the singular value sequence of the NNR estimator $\hat{\Theta}_{\mathrm{nuc}}$, i.e., $\psi_{1}(\hat{\Theta}_{\mathrm{nuc}})\geq \psi_{2}(\hat{\Theta}_{\mathrm{nuc}}) \geq \ldots \psi_{\min\{N, T\}}(\hat{\Theta}_{\mathrm{nuc}})$. We determine $\hat{R}$ by eigenvalue-ratio test (\citealp{ahn2013eigenvalue}) with a preset $R_{\max}>0$:
\begin{align*}
\hat{R} = \operatorname*{argmax}_{r = 1,2,\ldots, R_{\max}} \psi_{r}(\hat{\Theta}_{\mathrm{nuc}})/ \psi_{r+1}(\hat{\Theta}_{\mathrm{nuc}})
\end{align*}
\end{itemize}
\end{algorithm}
Algorithm~\ref{algorithm:tuning} is computationally simple to implement.
Step~1 and Step~2, which solves a nonlinear panel regression with additive fixed effects, is a convex optimization problem for which efficient algorithms are readily available. In addition, the procedure avoids the use of cross-validation for tuning-parameter selection, thereby substantially reducing computational cost.
\subsection{Bias correction}
The FE estimator suffers from the incidental parameter problem because the number of fixed effects grows with $N$ and $T$. As a result, $\hat{\beta}_{\mathrm{FE}}$ exhibits a non-negligible bias of the same order as its standard deviation. Bias correction is therefore required to obtain a $\sqrt{NT}$-consistent estimator centered at the true parameter.
\paragraph{Analytical bias correction}
We follow \citet{chen2021nonlinear} to perform analytical bias correction. For any $d = 1,2,\ldots, d_X$, let
\begin{align*}
(\Lambda^*_{d}, \Gamma^*_{d}) \in \operatorname*{argmin}_{\Lambda_{d}\in \mathbb{R}^{N\times R}, \Gamma_{d}\in \mathbb{R}^{T\times R}}\sum_{i=1}^{N}\sum_{t=1}^{T}\mathbb{E}(-\ddot{\ell}^0_{it})\left(\frac{\mathbb{E}(\ddot{\ell}^0_{it}X_{it, d})}{\mathbb{E}(\ddot{\ell}^0_{it})} - \lambda_{d, i}' \gamma_{0, t} - \lambda_{0, i}' \gamma_{d, t}\right)^2
\end{align*}
and define
\begin{align*}
\Xi_{d, it} = \lambda_{d, i}^{*\prime} \gamma_{0, t} + \lambda_{0, i}' \gamma^*_{d, t}, \quad \tilde{X}_{d, it} = X_{d, it} - \Xi_{d, it}
\end{align*}
Here, $\ddot{\ell}_{it}^0: = \partial^2 \log f(\beta_0'X_{it} + \lambda_0'\gamma_0 )/\partial Y^{*2} $ denotes the second derivative of the log-likelihood $\ell_{it}(\cdot)$ with respect to the index, evaluated at the true parameter values. In addition, let $\hat{\Xi}_{it}$ denote the sample analogue of $\Xi_{it}$, and define
\begin{align*}
\widehat{B} & = -\frac{1}{N} \sum_{i=1}^{N} \sum_{t=1}^{T}\widehat{\gamma}'_{\mathrm{local}, t} \left(\sum_{\tau=1}^{T}\widehat{\gamma}_{\mathrm{local}, \tau}\widehat{\gamma}_{\mathrm{local}, \tau}'\widehat{\ddot{\ell}}_{i\tau}\right)^{-1}\widehat{\gamma}_{\mathrm{local}, t} \left( \widehat{\dot{\ell}}_{it}\widehat{\ddot{\ell}}_{it}(X_{it} - \widehat{\Xi}_{it}) + \frac{1}{2}\widehat{\dddot{\ell}}_{it}(X_{it} - \widehat{\Xi}_{it})\right) \\
\widehat{D} & = -\frac{1}{T} \sum_{t=1}^{T} \sum_{i=1}^{N} \widehat{\lambda}'_{\mathrm{local}, i} \left(\sum_{j=1}^{N}\widehat{\lambda}_{\mathrm{local}, j}\widehat{\lambda}_{\mathrm{local}, j}'\widehat{\ddot{\ell}}_{jt}\right)^{-1}\widehat{\lambda}_{\mathrm{local}, i}\left(\widehat{\dot{\ell}}_{it}\widehat{\ddot{\ell}}_{it}(X_{it} - \widehat{\Xi}_{it}) + \frac{1}{2}\widehat{\dddot{\ell}}_{it}(X_{it} - \widehat{\Xi}_{it})\right) \\
\widehat{W} & = -\frac{1}{NT} \sum_{i=1}^{N}\sum_{t=1}^{T} \widehat{\ddot{\ell}}_{it}(X_{it} - \widehat{\Xi}_{it})(X_{it} - \widehat{\Xi}_{it})'
\end{align*}
Here, $\widehat{\dot{\ell}}_{it}$, $\widehat{\ddot{\ell}}_{it}$, $\widehat{\dddot{\ell}}_{it}$ denote the first, second, and third derivatives of the log-likelihood $\ell_{it}(\cdot)$ with respect to the index, evaluated at the local estimator $(\hat{\beta}_{\mathrm{local}}, \hat{\Lambda}_{\mathrm{local}}, \hat{\Gamma}_{\mathrm{local}})$. The analytical bias correction estimator, $\hat{\beta}_{\mathrm{ABC}} := \hat{\beta}_{\mathrm{local}} - \frac{1}{T} \widehat{W}^{-1} \widehat{B} - \frac{1}{N} \widehat{W}^{-1} \widehat{D}$, follows
\begin{align*}
\sqrt{NT}\left(\hat{\beta}_{\mathrm{ABC}} - \beta_0\right)\stackrel{d}{\longrightarrow} N(0, \widehat{W}^{-1})
\end{align*}
\paragraph{Sample-splitting Jackknife} The sample-splitting Jackknife estimator in \citet{chen2021nonlinear} is given by
\begin{align*}
\hat{\beta}_{JBC} := 3 \hat{\beta}_{\mathrm{local}} - \overline{\beta}_{N, T/2} - \overline{\beta}_{N/2, T}
\end{align*}
where $\overline{\beta}_{N, T/2}$ is the average of the estimators in the half-panels $\{(i, t)\mid i=1, \ldots, N, t = 1,\ldots, \lceil T/2 \rceil \}$ and $\{(i, t)\mid i=1, \ldots, N, t = \lceil T/2 \rceil \ + 1,\ldots, T\}$, $\overline{\beta}_{N/2, T}$ is the average of the estimators in the half-panels $\{(i, t)\mid i=1, \ldots, \lceil N/2 \rceil , t = 1,\ldots, T\}$ and $\{(i, t)\mid i=\lceil N/2 \rceil + 1\ldots, N , t = 1,\ldots, T\}$.
\section{Monte Carlo and Empirical Application}\label{sec:MC}
In this section, we conduct Monte Carlo simulations to evaluate the finite-sample performance of our estimator. We also provide an empirical application to the gravity equation.
\subsection{Monte Carlo}
\paragraph{Binary response model} We evaluate the performance of our estimator in binary response models with strictly exogenous covariates. Let the data be generated from the following binary response model with $R = 2$:
\begin{equation}\label{eq:logit_static}
\begin{aligned}
Y_{it} & = \boldsymbol{1}\left(\beta_{1}X_{it} + \lambda_{i}' \gamma_{t} + \epsilon_{Y, it}>0\right) \\
X_{it} & = \lambda_{i}' \gamma_{t} + \lambda_{i}' \iota + \gamma_{ t}' \iota + \lambda_{X, i} \gamma_{X, t} + \epsilon_{X, it}
\end{aligned}
\end{equation}
The error terms $\{\epsilon_{Y, it}\}_{1\leq i\leq N, 1\leq t\leq T}$ are i.i.d. random variables across both $i$ and $t$. We consider two designs: (i) \emph{Probit}, where $\epsilon_{Y, it}$ follows a standard normal distribution $N(0, 1)$; and (ii) \emph{Logit}, where $\epsilon_{Y, it}$ follows a standard logistic distribution.
In both designs, $ \lambda_{i} = (\lambda_{i1}, \lambda_{i2})'$, $ \gamma_{t} = (\gamma_{t1}, \gamma_{t2})'$, $\{ \lambda_{ir}\}_{\{ 1\leq i\leq N, r=1,2\}}$, $\{ \gamma_{tr}\}_{\{ 1\leq t\leq T, r=1,2\}}$, $\{ \lambda_{X, i}\}_{\{ 1\leq i\leq N\}}$, and $\{ \gamma_{X, t}\}_{\{ 1\leq t\leq T\}}$ consist of independent random variables drawn from the standard normal distribution $N(0, 1)$.
In addition, for both designs, $\{\epsilon_{X, it}\}_{1\leq i\leq N, 1\leq t\leq T}$ are i.i.d. random variables across both $i$ and $t$ following a normal distribution $N(0, 4)$. We set $\alpha = 0.05$ to compute the tuning parameter $\varphi_{NT}$ and $R_{\max} = 5$ to select the number of factors. In addition, $\beta_{1} = 0.2$. It is straightforward to verify that conditional on $\Lambda_0$ and $\Gamma_0$, $X_{it}$ is a strictly exogenous covariate.
For each design, we conduct $1000$ Monte Carlo replications to evaluate the finite-sample performance of our estimator across different sample sizes, ranging from $(N,T) = (50,40)$ to $(N,T) = (1000,200)$. We also report alternative estimators for comparison, as well as bias-corrected estimators to examine whether bias corrections based on our two-step estimator mitigate the incidental parameter problem.
The estimators considered in the Monte Carlo experiments include a pooled estimator that ignores individual and time latent factors (POOL); our first-step nuclear norm-regularized estimator (NNR), defined in~\eqref{eq:nnr_definition}; our two-step estimator using the true number of factors $R$ ($\mathrm{FE}^*$); and its analytical and jackknife bias-corrected versions ($\mathrm{ABC}^*$ and $\mathrm{JBC}^*$).
We also report the corresponding two-step estimator using an estimated number of factors $\hat{R}$ ($\mathrm{FE}$), along with its analytical and jackknife bias-corrected counterparts ($\mathrm{ABC}$ and $\mathrm{JBC}$).
We report both the bias and the standard deviation of each estimator. In addition, we provide the average estimated number of factors $\hat{R}$.
Table~\ref{tab:probit_static} reports results for the Probit models. The first column (POOL) presents the pooled regression results, which exhibit substantial bias that does not diminish as the sample size increases. The second column (NNR) reports the performance of the NNR estimator, whose bias decreases slowly toward zero as $N, T \to \infty$, consistent with Theorem~\ref{thm:consistency}.
Table~\ref{tab:probit_static} shows that our proposed two-step estimator using the true number of factors ($R = 2$), reported in column $\mathrm{FE}^*$, substantially improves upon the NNR estimator: its bias is smaller and converges rapidly to zero as $N$ and $T$ increase. Nevertheless, due to the incidental parameter problem, the bias and standard deviation of the FE estimator remain of the same order even for large sample sizes, such as $(N, T) = (1000, 200)$. The table also shows that applying bias correction to our two-step estimator effectively reduces bias. The analytical bias-corrected estimator ($\mathrm{ABC}^*$) significantly reduces bias even in small samples, such as $(N, T) = (50, 40)$. The Jackknife sample-splitting estimator ($\mathrm{JBC}^*$), while less effective when $T = 40$, achieves effective bias reduction for larger $N$ and $T$. When the number of factors is estimated rather than known, the corresponding estimators ($\mathrm{FE}$, $\mathrm{ABC}$, and $\mathrm{JBC}$) continue to perform well. Since we can estimate the number of factors with high accuracy (see the last column $\bar{R}$), the performance of estimators based on $\hat{R}$ is very similar to that obtained using the true $R$. Overall, the two-step estimator delivers strong finite-sample performance and, when combined with bias correction, achieves substantial bias reduction and supports valid inference.
Table~\ref{tab:logit_static} reports results for the Logit models. Overall, the results are very similar to those for the Probit model: our proposed two-step estimator delivers strong finite-sample performance and, when combined with bias correction, achieves substantial bias reduction and supports valid inference. One notable difference arises in small samples. When $(N, T) = (50, 40)$, the analytical bias-corrected ($\mathrm{ABC}^*$ and $\mathrm{ABC}$) estimators, although still substantially reduces bias, perform less well in the Logit case.
\paragraph{More simulations} Additional numerical results related to dynamic models are presented in Appendix~\ref{sec:extension}.
\begin{table}[H]
\centering
\caption{Simulation Results: Probit Model}\label{tab:probit_static}
\begin{tabular}{cccccccccc}
\toprule
& POOL & NNR & $\mathrm{FE}^*$ & $\mathrm{ABC}^*$ & $\mathrm{JBC}^*$ & FE & ABC & JBC & $\bar{R}$ \\
& $( \times 10^{-2})$ & $( \times 10^{-2})$ & $(10^{-2})$ & $(\times 10^{-2})$ & $(\times 10^{-2})$ & $(\times 10^{-2})$ & $(\times 10^{-2})$ & $(\times 10^{-2})$ & \\
\midrule
N = 50, T = 40 & & & & & & & & & \\
BIAS & 2.14 & 3.69 & 3.78 & 0.78 & -3.05 & 3.78 & 0.85 & -2.87 & 1.962 \\
STD & (1.82) & (1.85) & (2.53) & (2.15) & (3.61) & (2.51) & (2.20) & (3.74) & \\
N = 100, T = 40 & & & & & & & & & \\
BIAS & 2.04 & 3.25 & 2.29 & 0.18 & -1.56 & 2.30 & 0.18 & -1.55 & 1.999 \\
STD & (1.46) & (1.32) & (1.59) & (1.42) & (2.13) & (1.59) & (1.42) & (2.13) & \\
N = 200, T = 40 & & & & & & & & & \\
BIAS & 2.06 & 3.08 & 1.79 & 0.05 & -0.55 & 1.79 & 0.05 & -0.55 & 2.000 \\
STD & (1.35) & (1.11) & (1.06) & (0.96) & (1.36) & (1.06) & (0.96) & (1.36) & \\
N = 100, T = 100 & & & & & & & & & \\
BIAS & 2.07 & 2.70 & 1.17 & 0.07 & -0.34 & 1.17 & 0.07 & -0.34 & 2.000 \\
STD & (1.10) & (0.91) & (0.86) & (0.81) & (0.97) & (0.86) & (0.81) & (0.97) & \\
N = 200, T = 100 & & & & & & & & & \\
BIAS & 2.01 & 2.41 & 0.83 & 0.02 & -0.13 & 0.83 & 0.02 & -0.13 & 2.000 \\
STD & (0.94) & (0.68) & (0.60) & (0.58) & (0.65) & (0.60) & (0.58) & (0.65) & \\
N = 200, T = 200 & & & & & & & & & \\
BIAS & 1.98 & 2.14 & 0.53 & 0.00 & -0.07 & 0.53 & 0.00 & -0.07 & 2.000 \\
STD & (0.75) & (0.49) & (0.41) & (0.40) & (0.44) & (0.41) & (0.40) & (0.44) & \\
N = 1000, T = 200 & & & & & & & & & \\
BIAS & 1.98 & 1.93 & 0.31 & 0.00 & -0.01 & 0.31 & 0.00 & -0.01 & 2.000 \\
STD & (0.57) & (0.28) & (0.18) & (0.18) & (0.25) & (0.18) & (0.18) & (0.25) & \\
\bottomrule
\end{tabular}
\vspace{0.5cm}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} Monte Carlo results based on $1000$ replications for the Probit model with strictly exogenous covariates as in ~\eqref{eq:logit_static}. We set $\alpha = 0.05$ to determine the tuning parameter $\varphi_{NT}$. We report bias and standard deviation for pooled estimators (POOL), nuclear norm regularized estimators (NNR). Using the true number of factors, $R=2$, we report bias and standard deviation (measured in units of $\times 10^{-2}$) for our proposed estimators ($\mathrm{FE}^*$), analytical bias-corrected estimators based on $\mathrm{FE}^*$, ($\mathrm{ABC}^*$) and Jackknife bias-corrected estimators based on $\mathrm{FE}^*$, ($\mathrm{JBC}^*$). Using estimated number of factors $\hat{R}$, we report bias and standard deviation (measured in units of $\times 10^{-2}$) for our proposed estimators ($\mathrm{FE}$), analytical bias-corrected estimators based on $\mathrm{FE}$, ($\mathrm{ABC}$) and Jackknife bias-corrected estimators based on $\mathrm{FE}^*$, ($\mathrm{JBC}^*$). Furthermore, we report the average estimated number of factors $\bar{R}$.
\end{minipage}
\end{table}
\begin{table}[H]
\centering
\caption{Simulation Results: Logit Model}\label{tab:logit_static}
\begin{tabular}{cccccccccc}
\toprule
& POOL & NNR & $\mathrm{FE}^*$ & $\mathrm{ABC}^*$ & $\mathrm{JBC}^*$ & FE & ABC & JBC & $\bar{R}$ \\
& $( \times 10^{-2})$ & $( \times 10^{-2})$ & $(10^{-2})$ & $(\times 10^{-2})$ & $(\times 10^{-2})$ & $(\times 10^{-2})$ & $(\times 10^{-2})$ & $(\times 10^{-2})$ & \\
\midrule
N = 50, T = 40 & & & & & & & & & \\
BIAS & 7.51 & 6.89 & 5.40 & 3.64 & -3.32 & 5.45 & 3.70 & -3.21 & 1.802 \\
STD & (2.48) & (2.40) & (3.64) & (3.48) & (6.56) & (3.68) & (3.52) & (6.75) & \\
N = 100, T = 40 & & & & & & & & & \\
BIAS & 7.57 & 6.64 & 3.09 & 1.70 & -2.25 & 3.10 & 1.71 & -2.23 & 1.959 \\
STD & (1.99) & (1.88) & (2.39) & (2.33) & (4.01) & (2.40) & (2.34) & (4.02) & \\
N = 200, T = 40 & & & & & & & & & \\
BIAS & 7.46 & 6.26 & 1.89 & 0.69 & -1.24 & 1.89 & 0.69 & -1.24 & 1.999 \\
STD & (1.67) & (1.50) & (1.56) & (1.50) & (2.33) & (1.56) & (1.50) & (2.33) & \\
N = 100, T = 100 & & & & & & & & & \\
BIAS & 7.35 & 5.82 & 1.23 & 0.32 & -0.90 & 1.23 & 0.32 & -0.90 & 2.000 \\
STD & (1.40) & (1.22) & (1.17) & (1.12) & (1.51) & (1.17) & (1.12) & (1.51) & \\
N = 200, T = 100 & & & & & & & & & \\
BIAS & 7.36 & 5.39 & 0.76 & 0.11 & -0.29 & 0.76 & 0.11 & -0.29 & 2.000 \\
STD & (1.23) & (1.01) & (0.88) & (0.85) & (1.07) &(0.88) & (0.85) & (1.07) & \\
N = 200, T = 200 & & & & & & & & & \\
BIAS & 7.34 & 4.87 & 0.51 & 0.08 & -0.09 & 0.51 & 0.08 & -0.09 & 2.000 \\
STD & (0.95) & (0.74) & (0.59) & (0.57) & (0.68) & (0.59) & (0.57) & (0.68) & \\
N = 1000, T = 200 & & & & & & & & & \\
BIAS & 7.36 & 4.07 & 0.27 & 0.02 & -0.01 & 0.27 & 0.02 & -0.01 & 2.000 \\
STD & (0.70) & (0.45) & (0.27) & (0.26) & (0.34) & (0.27) & (0.26) & (0.34) & \\
\bottomrule
\end{tabular}
\vspace{0.5cm}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} Monte Carlo results based on $1000$ replications for the Logit model with strictly exogenous covariates as in ~\eqref{eq:logit_static}. We set $\alpha = 0.05$ to determine the tuning parameter $\varphi_{NT}$. We report bias and standard deviation (measured in units of $\times 10^{-2}$) for pooled estimators (POOL), nuclear norm regularized estimators (NNR). Using the true number of factors, $R=2$, we report bias and standard deviation for our proposed estimators ($\mathrm{FE}^*$), analytical bias-corrected estimators based on $\mathrm{FE}^*$, ($\mathrm{ABC}^*$) and Jackknife bias-corrected estimators based on $\mathrm{FE}^*$, ($\mathrm{JBC}^*$). Using estimated number of factors $\hat{R}$, we report bias and standard deviation (measured in units of $\times 10^{-2}$) for our proposed estimators ($\mathrm{FE}$), analytical bias-corrected estimators based on $\mathrm{FE}$, ($\mathrm{ABC}$) and Jackknife bias-corrected estimators based on $\mathrm{FE}^*$, ($\mathrm{JBC}^*$). Furthermore, we report the average estimated number of factors $\bar{R}$.
\end{minipage}
\end{table}
\subsection{Empirical application}
We revisit the empirical analysis of the gravity equation in trade, which employs a count data model to examine the determinants of bilateral trade flows. This setting is also used as the empirical application in \citet{chen2021nonlinear}. The data, originally from \citet{helpman2008estimating}, include bilateral trade flows and other relevant variables for $N = 157$ countries. In our application, we focus on the year $1986$. Given its network structure, we index the countries as $i = 1, 2, \ldots, 157$ (exporters) and $j = 1, 2, \ldots, 157$ (importers). The sample size (the number of distinct country pairs) is therefore $157 \times 156 = 24{,}492$.
The outcome variable $Y_{ij}$ represents the volume of trade (in thousands of constant $2000$ US dollars) from country $i$ to country $j$. The covariates $X_{ij}$ include key determinants of bilateral trade flows, such as the logarithm of the distance between the capitals of the two countries (Log distance) and binary indicators for shared borders (Border), legal systems (Legal), common languages (Language), colonial ties (Colony), currency unions (Currency), regional free-trade agreements (FTAs), and religions (Religion). Descriptive statistics are presented in Table~\ref{tab:summary}.
\begin{table}[H]
\centering
\caption{Summary Statistics}
\begin{tabular}{lcc}
\toprule
& Mean & Standard deviation \\
\midrule
\multicolumn{1}{l}{Trade Volume} & 84,542 & 1,082,219 \\
\multicolumn{1}{l}{Log distance} & 4.18 & 0.78 \\
\multicolumn{1}{l}{Border} & 0.02 & 0.13 \\
\multicolumn{1}{l}{Legal} & 0.37 & 0.48 \\
\multicolumn{1}{l}{Language} & 0.29 & 0.45 \\
\multicolumn{1}{l}{Colony} & 0.01 & 0.10 \\
\multicolumn{1}{l}{Currency} & 0.60 & 1.37 \\
\multicolumn{1}{l}{FTA} & 0.01 & 0.08 \\
\multicolumn{1}{l}{Religion} & 0.17 & 0.25 \\
\bottomrule
\end{tabular}
\vspace{0.5cm}
\begin{minipage}{\textwidth}
\centering
\footnotesize
\textbf{Note:} The table is from \cite{helpman2008estimating}.
\end{minipage}
\label{tab:summary}
\end{table}
We consider the estimation of the following Poisson model:
\begin{align*}
\mathbb{E}(Y_{ij}\mid X_{ij}, \lambda_{1, i} , \gamma_{1, j}, \lambda_{2, i}, \gamma_{2, j} ) = \mathrm{Poisson}(\exp\{\beta' X_{ij} + \lambda_{1, i} + \gamma_{1, j} + \lambda'_{2, i}\gamma_{2, j}\})
\end{align*}
Consistent with \citet{chen2021nonlinear}, we explicitly include additive two-way fixed effects $\lambda_{1,i}+\gamma_{1,j}$ in the regression, together with an interactive fixed-effect component $\lambda_{2,i}'\gamma_{2,j}$. Incorporating additive fixed effects requires only minor modifications to the algorithm (the NNR regression and the gradient-descent method) described in Section~\ref{sec:implementation}. In addition, the dimension of the interactive component (i.e., the number of factors in $R = \mathrm{dim}(\lambda_{2,i})$), as well as the tuning parameters, can be selected using the similar approach as in Section~\ref{sec:implementation}. The corresponding optimization problem and algorithm are detailed in Appendix~\ref{sec:extension}, with the implementation provided in our accompanying R package.
In Table~\ref{tab:trade}, we report estimates from our two-step estimator (column $\mathrm{FE}$), an estimator that incorporates only two-way fixed effects $\lambda_{1,i}+\gamma_{1,j}$ (column $\mathrm{TWFE}$), and the estimator proposed in \citet{chen2021nonlinear}\footnote{
\citet{chen2021nonlinear} obtain their estimator using an EM-type algorithm, but do not provide a discussion of how to choose appropriate initial values.
} (column $\mathrm{CFW}$). Using our algorithm, the estimated number of factors is $\hat{R}=2$. For comparability, we therefore report the estimator of \citet{chen2021nonlinear} using $\hat{R}=2$ only. Standard deviations are reported in brackets.
Column $\mathrm{FE}$ is included as a reference, showing that the signs of most coefficients remain stable after introducing additive latent factors into the gravity equation. Specifically, the results indicate that after controlling for additive latent factors, bilateral trade volumes are larger between countries that are geographically closer, share a common border, have the same legal system, colonial ties, a common currency, or are linked by a free-trade agreement. These findings are consistent with theoretical predictions and previous empirical evidence in the gravity literature. The exceptions are the signs before the common religion and the common language. We do not provide a detailed discussion of the impact of latent factors, as this has already been extensively discussed in \citet{chen2021nonlinear}. Instead, we focus on the differences between columns $\mathrm{CFW}$ and $\mathrm{FE}$.
The estimators in columns $\mathrm{CFW}$ and $\mathrm{FE}$ are very similar. Although both estimates of the effect of a common language are negative (neither is statistically significant), both methods yield positive effects for having the same religion (both are statistically significant). The only notable difference between our estimator and that of \citet{chen2021nonlinear} is the estimated effect of a common currency. Our two-step estimator reports a coefficient of $1.25$ with a standard deviation of $0.34$, while \citet{chen2021nonlinear} report an estimate of $1.37$ with a standard deviation of $0.41$. However, this difference is minor, and both methods yield statistically significant estimates at the $5\%$ level. Finally, the log-likelihood for our estimator is $0.6711$, slightly lower than the $0.6714$ reported by \citet{chen2021nonlinear}. This small discrepancy is likely due to finite-sample variation.
\begin{table}[H]
\centering
\caption{Empirical Application: Gravity Equation}
\begin{tabular}{lccc}
\toprule
& TWFE & CFW & FE \\
& R = 0 & R = 2 & R = 2 \\
\midrule
\multicolumn{1}{c}{Log distance} & -0.64 & -0.71 & -0.71 \\
& (0.07) & (0.06) & (0.05) \\
\multicolumn{1}{c}{Border} & 0.71 & 0.32 & 0.32 \\
& (0.16) & (0.05) & (0.06) \\
\multicolumn{1}{c}{Legal} & 0.30 & 0.26 & 0.26 \\
& (0.06) & (0.04) & (0.04) \\
\multicolumn{1}{c}{Language} & -0.17 & -0.02 & -0.02 \\
& (0.10) & (0.06) & (0.06) \\
\multicolumn{1}{c}{Colony} & 0.36 & 0.39 & 0.39 \\
& (0.12) & (0.09) & (0.10) \\
\multicolumn{1}{c}{Currency} & 0.60 & 1.37 & 1.25 \\
&(0.09) & (0.41) & (0.34) \\
\multicolumn{1}{c}{FTA} & 0.25 & 0.17 & 0.17 \\
& (0.13) & (0.07) & (0.06) \\
\multicolumn{1}{c}{Religion} & -0.25 & 0.24 & 0.24 \\
& (0.12) & (0.13) & (0.08) \\
\midrule
Log-Likelihood & -0.44 & 0.67 & 0.67 \\
\bottomrule
\end{tabular}
\vspace{0.5cm}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} We report estimators and their standard deviations (in brackets) for three regression: Column $\mathrm{TWFE}$ represents the estimation for only including additive fixed effects. Column $\mathrm{CFW}$ is the result for $R = 2$ reported in \cite{chen2021nonlinear} using EM algorithm. Column $\mathrm{FE}$ reports our two-step estimator for $R = 2$. All standard deviations are robust to reciprocity. The log-likelihood in Column (2) is $0.6714$, while our log-likelihood is slightly smaller, $0. 6711$.
\end{minipage}
\label{tab:trade}
\end{table}
\section{Conclusion}\label{sec:conclusion}
We develop a computationally efficient two-step estimator for nonlinear panel models with interactive unobserved effects, thereby helping bridge the gap between recent theory and empirical practice. Our procedure first solves a convex relaxation of the log-likelihood function via nuclear-norm regularization to obtain a consistent initial value, and then refines it with gradient descent on the log-likelihood function. We show that the negative log-likelihood function is locally convex in a shrinking neighborhood of the global optimum and establish a convergence rate for the first-step NNR estimator that is fast enough to enter this neighborhood, ensuring that the second step converges to the global minimum. Consequently, our final estimator is asymptotically equivalent to the fixed-effects estimator of \citet{chen2021nonlinear} while avoiding the high-dimensional non-convex optimization burden. This equivalence allows researchers to conduct inference using their theoretical results, making the estimator practically feasible.
\bibliographystyle{aea}
\bibliography{ref}
\clearpage