Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
97,595 characters · 19 sections · 13 citation commands
-0.5cmLow-Rank Estimation of Nonlinear Panel Data Models
\onehalfspacing Unobserved heterogeneity is of broad interest in both reduced-form and structural work in economics and other social sciences. Empirical data often include individual units sampled over time from diverse backgrounds, with factors unobserved by econometricians. Accommodating this heterogeneity in a flexible, yet parsimonious manner is challenging but essential in practice.
Panel data, which involve observations on individual units over time, provide a valuable framework to model latent structures within low-dimensional manifolds. This is often achieved by incorporating unobserved individual and time effects into the model, controlling for unobserved covariates that remain either time or cross-sectionally invariant. A fixed effects approach imposes no distributional assumptions on these unobserved effects, allowing them to be arbitrarily related to observed covariates. A notable example is two-way fixed effects estimation for difference-in-differences, which has become a leading method in applied economics.
The additive structure, however, assumes that the influence of unobserved factors on individual units remains constant over time, failing to capture more complex dynamics. Alternatively, fixed effects can interact multiplicatively, giving rise to interactive fixed effects (IFEs) or factor structures. This multiplicative form provides a more flexible representation of heterogeneity, as it allows common time-varying shocks (factors) to affect cross-sectional units with individual-specific sensitivities (factor loadings). For example, they can account for aggregate shocks with heterogeneous impacts on agents in macroeconomic models or capture multidimensional individual heterogeneity with time-varying effects in microeconomic models. This flexibility motivated the discussion of interactive effects in the econometrics literature. \citet*{bai_panel_2009} and \citet*{moon_linear_2015,moon_dynamic_2017} study the linear regression model, treating individual and time effects as nuisance parameters estimated by least squares. \citet*{chen_nonlinear_2021} and \citet*{wang_maximum_2022} extend the least squares estimator to nonlinear panel data models with convex log-likelihood functions.
However, many important and widely used economic models, such as random coefficient logit panel models, do not have convex objective functions. The problem becomes even more challenging when interactive fixed effects are introduced to capture unobserved heterogeneity. Interactive fixed effect panel models with non-convex objective functions remain largely unexplored in the literature. The current approaches in the literature are not applicable as they rely on closed-form solutions or global convexity.
The present paper proposes a new method for econometric estimation and inference in panel models with interactive fixed effects. The main estimation procedure consists of two steps. In the first step, we obtain preliminary consistent estimators for all coefficients via a nuclear-norm regularization (NNR) procedure. The second step uses an iterative procedure to conduct post-nuclear-norm regularization estimation for the parameters of interest, treating our consistent regularized estimator as the initial value. We demonstrate a faster convergence of our iterative estimator compared to the preliminary estimator. Furthermore, we derive the asymptotic distribution of our iterative estimator.
In this context, I make two contributions. First, our general M-estimation framework accommodates models defined by potentially non-convex objective functions with respect to the unknown parameters of interest, while the results in \citet*{bai_panel_2009} and \citet*{chen_nonlinear_2021} rely on closed-form solutions or global convexity. Our theoretical findings, therefore, broaden the applicability of IFEs in econometrics, complementing the existing toolbox for applied economists when a factor structure is central to the analysis. Additionally, when the objective function is convex in the index, the proposed estimator formulates and solves a convex optimization problem, avoiding the issue of multiple local minima. Moreover, it does not require prior knowledge of the number of factors, and we also propose a consistent estimator for the number of factors.
The second contribution is to establish the large sample properties of the second-step iterative estimator. Specifically, this estimator exhibits a contraction mapping property, and thus converges to the true parameter value in probability at a rapid rate. Demonstrating the numerical convergence properties is technically challenging, as it requires controlling both the asymptotic bias and variance of the iterative estimator at each iteration. The iterative procedure in our second-step estimation is conceptually similar to that of \citet*{moon_nuclear_2019} and \citet*{hong_profile_2023}, but we do not have closed-form solutions as they do in linear models. We show that the estimator obtained through the second step exhibits the same asymptotic distribution as found in \citet*{chen_nonlinear_2021} in the case of convex objective functions.
This paper relates to three branches of the literature. First, it adds to the extensive literature on panel data models with IFEs. \citet*{bai_panel_2009} and \citet*{moon_linear_2015,moon_dynamic_2017} propose estimators based on quasi-maximum likelihood estimation and principal component analysis. Much previous work has focused on linear models, while nonlinear models have recently attracted attention. \citet*{chen_nonlinear_2021} extend \citet*{fernandez-val_individual_2016}'s results to nonlinear panel data models with IFEs. \citet*{chen_quantile_2021} study quantile factor models, and \citet*{gao_binary_2023-1,ando_bayesian_2022} study binary panel choice models. Relatedly, \citet*{boneva_discrete-choice_2017} and \citet*{chen_common_2025} extend the common correlated effects (CCE) framework of \citet*{pesaran_estimation_2006} to nonlinear settings. We refer to \citet*{fernandez-val_fixed_2018} for a recent survey. We complement and extend the literature by generalizing previous results to a broad range of M-estimators.
Second, our work contributes to the burgeoning literature on nuclear norm penalization, a technique that has gained popularity in estimating low-rank matrices in statistics and econometrics. This approach has been explored in various works such as \citet*{agarwal_noisy_2012,agarwal_causal_2021,athey_matrix_2021,belloni_high_2023,beyhum_square-root_2019,candes_exact_2009,chernozhukov_inference_2019-1,chernozhukov_inference_2023,fan_generalized_2017,feng_nuclear_2023,hong_profile_2023,huang_low-rank_2018,miao_high-dimensional_2022,ma_detecting_2021,negahban_estimation_2011,negahban_unified_2012}, among others. Most of these works focus on establishing estimation error bounds, specifically in the Frobenius norm, for the NNR estimators. The convergence rate of our estimator is consistent with the existing results of the mean and quantile regressions, indicating that NNR can be successfully extended to more complex models without compromising performance.
More recent studies have advanced statistical inference for NNR-based estimators. For instance, linear models have been studied by \citet*{armstrong_robust_2023, chernozhukov_inference_2023, miao_high-dimensional_2022, moon_nuclear_2019, hong_profile_2023}, among others. Two main debiasing strategies have emerged to address shrinkage bias. The first, inspired by the debiased-Lasso framework in the statistics literature, typically requires relatively few iterations. \citet*{chernozhukov_inference_2023} combine a rotation-based debiasing step with sample splitting, while \citet*{choi_inference_2024} show that sample splitting is not essential. \citet*{armstrong_robust_2023} apply minimax linear estimation theory to construct robust debiased estimators. The second line of research—where our paper contributes—employs iterative bias-correction procedures, as exemplified by \citet*{hong_profile_2023}, \citet*{miao_high-dimensional_2022}, and \citet*{moon_nuclear_2019}. Two very recent papers developed independently and in parallel with ours examine nonlinear panel data models with interactive fixed effects: \citet*{ChenMiaoSu2025} derive the asymptotic distribution of factor and loading estimators in logistic panel models without covariates, and \citet*{Zeleneev_tractable_2025} study single-index models under global convexity of the transformation function. In contrast, our results do not require global convexity, which allows us to accommodate models such as random-coefficients specifications.
Lastly, this paper is also related to the literature on estimating the parameters of high-dimensional models using $\ell_1$ regularization. Rather than listing all relevant papers, we refer the interested reader to the comprehensive textbook by \citet*{Hastie_lasso_2015} and the recent survey by \citet*{belloni_high-dimensional_2018}, and instead focus on a few key references. Specifically, \citet*{chetverikov_selecting_2025} develop a method called bootstrapping after cross-validation for selecting the penalty parameter in $\ell_1$-penalized M-estimators in high dimensions. \citet*{stadler_l1-penalization_2010} derive an oracle inequality for the $\ell_1$-penalized approach to estimating mixture regression models, while \citet*{beyhum_high-dimensional_2024} study high-dimensional nonconvex Lasso-type M-estimators and establish their convergence rates. However, none of these studies consider panel data settings or nuclear norm regularization.
An empirical illustration of our procedure is also provided, where we use multiple firm-level financial data sets to study US non-financial corporation joint financing decisions between debt issuance and equity repurchase, revisiting ma_nonfinancial_2019. The results show that our estimator produces economically sensible outcomes that align with financial intuition. For example, firms are more likely to issue equity and retire debt when the cost of debt is high and the cost of equity is low, consistent with the view that corporations act as cross-market arbitrageurs in their own securities. In addition, we document substantial heterogeneity in firm sensitivities to market valuation measures.
The remainder of the paper is organized as follows. Section (ref) introduces the class of IFE models. Section (ref) proposes a general M-estimation framework and presents the main estimation procedure. In Section (ref), we examine the asymptotic properties of the estimators, including the consistency and rate of convergence of the initial estimator and the rank estimator, as well as the convergence and asymptotic distribution of the second-step estimator. Section (ref) outlines the implementation algorithms and provides Monte Carlo results. Section (ref) contains the empirical applications. Finally, Section (ref) concludes. All proofs and additional results are provided in the Appendix.
\paragraph*{Notation.}
For a natural number $m\in\mathbb{N}$, we introduce the notation $[m]=\{1,\ldots,m\}$. For a vector $v=\left(v_{1},\ldots,v_{p}\right)^{\prime}\in\mathbb{R}^{p}$, we define its $\ell_{1}$-norm as $\|v\|_{1}=\sum_{j=1}^{p}\left|v_{j}\right|$, and its $\ell_{2}$-norm norm as $\|v\|=\sqrt{\sum_{j=1}^{p}v_{j}^{2}}$. For a matrix $A=\left(a_{it}\right)_{i,t}\in\mathbb{R}^{N\times T}$, we define its entry-wise $\ell_{1}$-norm as $\|A\|_{1}=\sum_{i=1}^{N}\sum_{t=1}^{T}\left|a_{it}\right|$, its Frobenius norm as $\|A\|_{F}=\sqrt{\sum_{i=1}^{N}\sum_{t=1}^{T}a_{it}^{2}}$, its infinity norm as $\|A\|_{\infty}=\max\left\{ \left|a_{it}\right|:i\in[N],t\in[T]\right\} $, its nuclear norm as $\|A\|_{*}=\operatorname{trace}\left(\sqrt{A^{\prime}A}\right)$, its spectral norm as $\|A\|=\sup_{x:\|x\|=1}\sqrt{x^{\prime}A^{\prime}Ax}$, and its rank by $\operatorname{rank}(A)$. Moreover, for a real matrix, we use $\sigma_s(\cdot)$ to denote its $s$-th largest singular value. For a real and symmetric matrix, we use $\mu_s(\cdot), \mu_{\max }(\cdot)$ and $\mu_{\min }(\cdot)$ to denote its $s$-th largest eigenvalue, the largest and smallest eigenvalues, respectively.
For a real number $a$, let $\operatorname{sgn}(a)=1$ if $a \geq 0$ and $\operatorname{sgn}(a)=-1$ if $a<0$. For a square matrix $A$ whose $j$-th diagonal element is denoted as $A_{j j}$, define $\operatorname{sgn}(A)$ as a diagonal matrix whose $j$-th diagonal element is equal to $\operatorname{sgn}\left(A_{j j}\right)$. We also define the operations $\vee$ and $\wedge$ as $a\vee b=\max\{a,b\}$ and $a\wedge b=\min\{a,b\}$. Finally, for sequences $\left\{ a_{m}\right\} _{m=1}^{\infty}$ and $\left\{ b_{m}\right\} _{m=1}^{\infty}$ we use $a_{m}\lesssim b_{m}$ as the shorthand for the inequality $a_{m}\leq\overline{c}b_{m}$ for some finite positive $\overline{c}$ independent of $a_{m}$ and $b_m$ for sufficiently large $m$. $a_{m}\asymp b_{m}$ means that $a_{m}\lesssim b_{m}$ and $b_{m}\lesssim a_{m}$.
Let the data observations be denoted as $\left\{ \left(Y_{it},X_{it}\right):i\in\left[N\right],t\in\left[T\right]\right\} $, where $Y_{it}\in\mathcal{Y}\subset\mathbb{R}$ is the outcome variable and $X_{it}\in\mathcal{X}\subset\mathbb{R}^{d_x}$ is a vector of exogenous covariates for a fixed $d_x$. The indices $i$ and $t$ denote individuals and time periods, respectively. In matrix notation, let $Y$ and $X_{j}$ as $N\times T$ matrices representing the outcome and the $j$-th covariate for $j=1,...,d_x$. Let $X=\left(X_{j}\right)_{j\in\left[d_x\right]}$ collect all covariates. The support of $\left(Y_{it},X_{it}\right)$ is given by $\mathcal{Y}\times\mathcal{X}$.
We assume that for each individual $i$ and time $t$, the outcome $Y_{it}$ is allowed to depend on both the observed covariates $X_{it}$ and the latent interactive effects given by individual-specific loadings $\lambda_{i}\in\mathbb{R}^{r}$ and time factors $f_{t}\in\mathbb{R}^{r}$, where $r$ is a fixed rank. The effects $\lambda_{i}$ and $f_t$, although unobserved by the econometrician, may confound the effect of $X_{it}$ on $Y_{it}$. Following the fixed effects approach, we treat the realizations of $\left\{\lambda_i\right\}_{i=1}^N$ and $\left\{f_t\right\}_{t=1}^T$ as unrestricted parameters to be estimated.
Formally, we consider a class of nonlinear panel data models in which the true values of the common parameter vector $\theta \in \mathbb{R}^p$ and the $N\times r$ and $T\times r$ matrices of fixed effects, $\Lambda=\left(\lambda_{1},...,\lambda_{N}\right)^{\prime}$ and $F=\left(f_{1},...,f_T\right)^{\prime}$, are defined by the solution to the population optimization problem,\footnote{The model is identified only up to a normalization; see Remark (ref). A formal discussion of the identification assumptions is provided in Assumption (ref).}
where $\ell:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ is a known loss function, which may be nonconvex with respect to $\theta$ and $\lambda_i'f_t$, and where $W_{it} = (Y_{it}, X_{it})$. The sets $\Theta \subset \mathbb{R}^p$, $\Phi_{\lambda} \subset \mathbb{R}$, and $\Phi_f \subset \mathbb{R}$ are parameter spaces. The expectation $\mathbb{E}[\cdot]$ is taken with respect to the distribution of the data, conditional on the fixed realizations of unobserved individual and time effects for all $N$ and $T$. It is also referred to as a factor model with factor loadings $\lambda_{i}$ and common factors $f_{t}$, and we will use the terms “factor” and “interactive fixed effect” synonymously. The conventional additive structure is a special case of the factor structure with $r=2$, $\lambda_{i}=\left(\lambda_{i1},1\right)^{\prime}$, and $f_{t}=\left(1,f_{1t}\right)^{\prime}$.
Let $\pi_{0,it} = \lambda_{0i}^{\prime} f_{0t}$ denote the true value of $\pi_{it}$ that generates the data, where $\pi_{0,it} \in \Phi \subset \mathbb{R}$ for each $i$ and $t$. The matrix $\Pi_0 \in \Phi^{N \times T}$ collects all the fixed effects. Similarly, the matrix $\Pi = (\pi_{it})$ is treated as a parameter to be estimated. Thus, we can rewrite Model (ref) as
Some examples of models that fall within this framework and the scope of our methodology are as follows.
In this subsection, we describe our two-step estimation procedure for the model: the initial estimator is based on nuclear norm regularization (NNR). It is followed by a second step iterative estimator. As we will show in Sections (ref) and (ref), under a set of regularity conditions, the first step estimation provides consistent estimators for $\theta_0$, $\Lambda_0$ and $F_0$. The second step uses the estimators from the first step as the initial values, and iteratively updates estimators to enhance their properties. The asymptotic normality of the iterative estimator is established in Section (ref).
We define the loss function as
Motivated by the literature on nuclear norm regularization, we propose estimating $\left(\theta_{0},\Pi_{0}\right)$ by minimizing the following penalized criterion function:
where $\nu$ is a tuning parameter,\footnote{The tuning parameter $\nu$ depends on $N$ and $T$, but we omit the subscript for simplicity.} and $\|\cdot\|_{*}$ denotes the nuclear norm, which regularizes the singular values of the matrix $\Pi$. The estimation problem at hand is inherently high-dimensional, as it involves estimating a total of $p+r\left(N+T\right)$ parameters, where the number of parameters grows linearly with $N$ and $T$. Specifically, $p$ represents the number of parameters associated with $\theta$, and $r(N + T)$ corresponds to those related to the low-rank matrix $\Pi$.
Our estimator can be viewed as a generalization of Lasso, with a key difference being that we impose sparsity not directly on individual parameters but on the singular values of the matrix $\Pi$. This regularization ensures that the number of non-zero singular values of $\Pi$ is relatively sparse compared to $N$ and $T$. Since we assume that the rank $r$ of the matrix $\Pi_0$ is fixed and much smaller than both $N$ and $T$, the nuclear norm regularization reduces the complexity by encouraging a low-rank structure in $\Pi$, similar to how Lasso enforces sparsity in individual coefficients.
Similar to Lasso, (ref) can be seen as a convex relaxation of the following problem:
where the rank $r$ is predetermined, and $\theta$, $\lambda_{i}$ and $f_{t}$ are jointly estimated, as studied by \citet*{bai_panel_2009} and \citet*{chen_nonlinear_2021}, among others. While the formulation in (ref) seems appealing, it presents two challenges.
First, the rank constraint makes (ref) a non-convex optimization problem, even when the loss function is convex. In contrast, our proposed estimator, defined in (ref), avoids directly penalizing or constraining the rank of the estimated interactive fixed effects matrix. Instead, it seeks a $\hat{\Pi}$ with a small nuclear norm, serving as a convex relaxation of $\eqref{ref:Model1}$. Just as $\ell_1$ minimization is the tightest convex relaxation of the combinatorial $\ell_0$ minimization problem, nuclear norm minimization offers the tightest convex relaxation of the NP-hard rank minimization problem \citep*{candes_matrix_2010}.
Second, when the loss function $\mathcal{L}_{NT}$ is non-convex, as in random coefficients logit panel models, there are no guarantees of convergence for the estimators defined in (ref). Global convexity plays a crucial role in the literature on interactive fixed effects. For instance, \citet*{chen_nonlinear_2021} rely on global convexity to establish the consistency of $\theta$ and to bound the remainder terms in their stochastic expansions. In contrast, our penalized estimator accommodates a broader class of important economic models without global convexity.
Note that the NNR-based initial estimation does not require the number of factors $r$ to be known beforehand. When $r$ is unknown, we propose estimating $r$ using singular value thresholding (SVT) as follows:
Alternatively, $r$ can be estimated using other methods, such as the panel information criteria (IC) and panel criteria (PC) methods of \citet*{bai_determining_2002}, or the eigenvalue ratio (ER) and growth ratio (GR) methods of \citet*{ahn_eigenvalue_2013}, among others. These methods remain valid as long as $\hat{\theta}$ is consistent. In Section (ref), we first establish the consistency of $\hat{\theta}$ and $\hat{\Pi}$, then we prove that $\hat{r}=r$ with probability approaching one (w.p.a.1). Consequently, for the second-step estimation, we assume $r$ is known.
It is widely recognized that NNR estimators are subject to shrinkage bias, which complicates statistical inference. Theorem (ref) establishes the convergence rate of the initial estimator, which, although consistent, is slower than the $\sqrt{NT}$ rate. To address this issue, we propose a post nuclear-norm-regularization estimator to improve the convergence rate.
For simplicity, we focus primarily on maximum likelihood estimation. We assume that the outcome is generated by \[ Y_{it}\mid\ensuremath{X_{it};\theta,\ensuremath{\lambda_{i},f_{t}\sim g\left(\cdot|X_{it};\theta,\lambda_{i},f_{t}\right)}} \] where $g\left(\cdot\right)$ is a known probability density function with respect to some dominating measure. The log-likelihood function is then given by \[ \ell\left(W_{it};\theta, \lambda_{i}^{\prime}f_{t}\right)=-\log g\left(Y_{it}\mid X_{it};\theta,\lambda_{i},f_{t}\right) \]
Since the common factors $F$ and factor loadings $\Lambda$ cannot be separately identified without imposing normalization (see Remark (ref)), we impose the same normalization on both $F$ and $\Lambda$ as in \citet*{bai_panel_2009}, without loss of generality. Specifically, we normalize such that $F^{\prime}F/T=\mathbb{I}_{r}$ and $\Lambda^{\prime}\Lambda/N$ is diagonal with non-increasing diagonal elements. It is important to note that the ultimate asymptotic results remain unchanged regardless of the specific normalization chosen for $F$ and $\Lambda$.
To refine the initial estimators, we employ a localized iterative procedure. Starting with the initial estimators $\left(\hat{\theta},\hat{\Lambda},\hat{F}\right)$, the procedure iteratively updates these estimators by searching within a local region of radius $d_{NT}=c\log\left(N\land T\right)\gamma_{NT}$, where $c$ is a positive constant, and $\gamma_{NT}$ represents the convergence rate of the first-step estimator, as specified in Corollary (ref) below. The procedure is defined as follows:
In this section, we examine the asymptotic properties of both the initial and iterative estimators. First, we study the consistency and convergence rate of the initial regularized estimator. We then establish the convergence and asymptotic distribution of the second-step iterative estimator. We consider an asymptotic framework where $N$ and $T$ tend to infinity jointly. To keep the exposition simple, we focus on the single-index model with interactive fixed effects as the leading case. Similar results can be derived for other models with multiple indices.
Recall the definition of $\mathcal{L}_{NT}$ in (ref). We further define \[ \mathcal{\bar{L}}\left(\theta,\Pi\right)=\mathbb{E}\left[\mathcal{L}_{NT}\left(\theta,\Pi\right)\right],\quad\mathcal{\tilde{L}}_{NT}\left(\theta,\Pi\right)=\mathcal{L}_{NT}\left(\theta,\Pi\right)-\mathbb{E}\left[\mathcal{L}_{NT}\left(\theta,\Pi\right)\right] \] and the excess risk function as \[ \ensuremath{\mathcal{E}}\left(\theta,\Pi\right)=\mathcal{\bar{L}}\left(\theta,\Pi\right)-\mathcal{\bar{L}}\left(\theta_{0},\Pi_{0}\right) \] Let $\underline{c}$ and $\overline{c}$ denote generic positive constants that may vary in their occurrences.
Assumption (ref) imposes basic regularity on the data generating process. Assumption (ref)(ii) limits inter-temporal dependence within the data. Specifically, $\mu$ controls the strength of this temporal dependence. We adopt the assumption of exponentially beta-mixing for simplicity, but similar results can be achieved with polynomially beta-mixing data. If the observations $\left\{ \left(W_{it}\right)\right\} _{i\in[N],t\in[T]}$ are i.i.d. across $i$ and over $t$, our theoretical results will hold without imposing Assumption (ref)(ii). Assumption (ref)(iii) imposes constraints on the relative rates at which $N$ and $T$ tend to infinity.\footnote{For the first-step estimator, it is sufficient that $T^{-1}\log N = o(1)$, which requires only that $T$ increase sufficiently fast relative to $\log N$. This weaker condition guarantees the convergence rate of the first-step estimates and is implied by the joint asymptotic regime $N/T \to \kappa^{2}$ adopted here.} It is the standard assumption in the large-$T$ panel data literature \citep*{hahn_asymptotically_2002,bai_panel_2009,chen_nonlinear_2021}. Assumption (ref)(iv) is included for convenience, but could be relaxed to allow covariates with sufficiently light tails. For instance, if the covariates are standard Gaussian, $c_x$ can grow like $\sqrt{\log(pNT)}$ with probability approaching one. This avenue is not pursued here to maintain a clear and concise exposition. Additionally, we focus only on the case of fixed $r$; the analysis of approximately low-rank $\Pi_0$, analogous to approximate sparsity in the Lasso literature, is beyond the scope of this paper.
Assumption (ref) imposes a Lipschitz continuity condition on the loss function with respect to the index. In the specific case where the loss function is globally Lipschitz with a constant independent of $y$ or $x$, this assumption holds trivially.
We show that the first-step regularized estimator is consistent, which is the starting point for establishing its rate of convergence. This section states the high-level conditions required for consistency, while Proposition (ref) and Appendix (ref) present primitive sufficient conditions.
This is an identification assumption that imposes restrictions on the shape of the excess risk function. When $N$ and $T$ are fixed, the condition is satisfied if the loss function is continuous and $(\theta_{0},\Pi_{0})$ is its unique minimizer. In the context of panel data, this condition needs to be satisfied uniformly in $N$ and $T$.
Thanks to the penalty term in (ref), we can show that with probability approaching one, $\left(\hat{\theta},\hat{\Pi}\right)$ belongs to a restricted set \[ \ensuremath{\mathcal{B}=\left\{ \left(\theta,\Pi\right)\in\Theta\times\Phi^{N\times T}:\left\Vert \Pi\right\Vert _{*}\leq c_{\mathcal{L}}\nu^{-1}+\left\Vert \Pi_{0}\right\Vert _{*}\right\} } \] where $c_{\mathcal{L}}=\mathcal{\bar{L}}\left(\theta_{0},\Pi_{0}\right)+1$. It is formally established in the Appendix. This result is important because it allows in the mathematical development to restrict the attention to a smaller set $\mathcal{B}$ included in the parameter space.
Assumption (ref) assumes uniform convergence. In the low-dimensional context, a similar condition is usually required on a compact set which does not depend on $N$ and $T$. The main difference here is that the radius of $\mathcal{B}$ grows with $N$ and $T$. We provide sufficient conditions for Assumption (ref) to hold in Proposition (ref).
The result depends on the condition that $\nu\left\Vert \Pi_{0}\right\Vert _{*}=o(1)$, which implies that the added penalty term has a negligible effect on the objective function when evaluated at the true values. This condition on $\nu$ is satisfied in the subsequent theorems and corollaries.
Similar to other high-dimensional settings such as those in negahban_unified_2012, we impose an invertibility condition involving another restricted set, denoted as $\mathcal{A}$, which we describe next and use to derive our main results. Before stating the next condition, we first introduce some notation.
Let $\Pi_{0}=UDV^{\prime}$ represent the singular value decomposition (SVD) of $\Pi_{0}$, where $U$ and $V$ are the matrices of singular vectors corresponding to all singular values. In particular, let $U_{0}\in\mathbb{R}^{N\times r}$ and $V_{0}\in\mathbb{R}^{T\times r}$ denote the columns of $U$ and $V$ associated with the non-zero singular values. For an $N\times T$ matrix $\Delta$, we define the operators
Essentially, $\mathcal{P}\left(\cdot\right)$ can be thought of as a projection onto the subspace spanned by the columns of $U_{0}$ and $V_{0}$, constituting the “ low-rank” space of $\Pi_{0}$. Similarly, $\mathcal{M}\left(\cdot\right)$ is the projection onto the space orthogonal to this low-rank space.
Next, we define the restricted set as
It can be shown that $\left(\hat{\theta}-\theta_{0},\hat{\Pi}-\Pi_{0}\right)$ lies in this “cone” under certain conditions. Namely, the set contains matrices $\Pi$ that are close to $\Pi_{0}$, in the sense that the part that cannot be explained by $\lambda_{0i}$ and $f_{0t}$ is small in terms of nuclear norm.
Assumption (ref) is based on Restricted Strong Convexity (RSC), which relaxes the definition of strong convexity by only needing strong convexity in certain directions or over a subset of the ambient space \citep*{negahban_unified_2012,wainwright_high-dimensional_2019}. This assumption represents a version of a widely used condition in the matrix completion literature, although verifying it typically requires imposing additional structure on the parameter space. See also the discussions in \citet*{moon_nuclear_2019} and \citet*{miao_high-dimensional_2022}, following their Assumptions 1 and 2, respectively, for the linear case. A detailed discussion for the single-index model is provided in Appendix (ref).
Finally, to control the empirical risk associated with the estimation errors, we define the set \[ \mathcal{V}=\left\{ (\theta,\Pi)\in\mathcal{B}:\left\Vert \delta\right\Vert ^{2}+\frac{1}{NT}\left\Vert \Delta\right\Vert _{F}^{2}\leq c_{l}\right\}. \] We also introduce the norm $\rho(\cdot,\cdot)$, defined as \[ \rho(\delta,\Delta)=\left[\left\Vert \delta\right\Vert ^{2}+\frac{1}{NT}\left\Vert \Delta\right\Vert _{*}^{2}\right]^{1/2}. \] According to Lemma (ref), it suffices to analyze the convergence rate of the empirical risk within the restricted set $\mathcal{V}$.
The proof of Proposition (ref) is based on the empirical process theory and does not rely on the differentiability of the loss function. It derives sufficient conditions under which Assumption (ref) holds and provides a stochastic bound that will be used for the subsequent theorem. Specifically, we have $c_{\varepsilon,NT}=O\left(\sqrt{\log\left(NT\right)/\left(N\wedge T\right)}\right)$, since, under our exponential $\beta$-mixing setting, the dependence length satisfies $c_T=O(\log (N T))$. The scaling sequence $\psi_{N T}$ is allowed to diverge slowly, at a rate exceeding $O(\log \left(\log (N T)\right))$, so that the probabilistic bound in Proposition (ref) holds uniformly over the local parameter space $\mathcal{V}$.
We now present our first main result for estimating $\left(\theta_{0},\Pi_{0}\right)$.
Theorem (ref) establishes the convergence rates of the estimation errors for $\hat{\theta}$ and $\hat{\Pi}$ in the $\ell_{2}$ norm. Here, $\theta$ is a low-dimensional parameter vector, while $\Pi$ is a high-dimensional matrix with $NT$ elements; hence, the Frobenius norm of $\Pi$ is normalized by $1/\sqrt{NT}$ to ensure comparability. For simplicity, consider the case of i.i.d. panel data across $i$ over $t$. When the regularization parameter is set to $\nu \asymp \frac{\sqrt{N}\lor\sqrt{T}}{NT}$, the convergence rate of our estimator under the Euclidean norm is of the order of $\sqrt{1 / (N \land T)}$, with all other factors ignored. These results are consistent with previous work on penalized mean and quantile regression models for panel data \citep*{athey_matrix_2021, feng_nuclear_2023, moon_nuclear_2019, belloni_high_2023}, while the present paper develops a unified M-estimation framework that extends these results to a broader class of models. Finally, note that Theorem (ref) does not require differentiability of the loss function; however, subsequent asymptotic normality results rely on smoothness conditions to establish limiting distributions.
Next, we impose an assumption on the common factors and factor loadings.
Assumption (ref) formalizes the strong factor condition commonly imposed in the literature; see, for example, \citet*{bai_panel_2009} and \citet*{chen_nonlinear_2021}. The requirement that the singular values $\sigma_1, \ldots, \sigma_r$ are distinct is analogous to Assumption G in \citet*{bai_inferential_2003}, and closely related to Assumption 4.2 in \citet*{chernozhukov_inference_2023}. In particular, Theorem (ref) does not depend on Assumption (ref).\footnote{\citet*{armstrong_robust_2023} use NNR and study robust inference for weak factors in the linear panel model.}
The following corollary establishes the consistency of $\hat{r}$ and the mean squared convergence rates of $\hat{\Lambda}$ and $\hat{F}$.
Nuclear-norm regularization does not require prior knowledge or specification of the number of factors $r$. Corollary (ref) indicates that $\hat{r}$ is a consistent estimator of $r$. Thus, in what follows, we assume that the number of factors has been correctly selected. Without loss of generality, we further assume that $\hat{\mathrm{S}}=\mathbb{I}_r$ to simplify the notation.
Here, we use the shorthand notation $\ell_{it}\left(\theta,\pi\right)=\ell\left(W_{it};\theta, \pi\right)$ for convenience, where $\theta\in\Theta$ and $\pi\in\Phi$. Additionally, we omit the function arguments when they are evaluated at the true parameter values $\left(\theta_{0},\pi_{0,it}\right)$, e.g., $\ell_{it}=\ell\left(W_{it};\theta_0, \pi_{0,it}\right)$.
To study the asymptotic behavior of the iterative estimator $\hat{\theta}^{\left(m+1\right)}$, it is necessary to introduce additional quantities that characterize how the parameter $\theta$ interacts with the incidental components $\pi_{it}$.
Let $\Xi_{it}$ denote a $p$-dimensional vector defined by the following population weighted least squares projection for each component of $\mathbb{E}\left[\partial_{\theta_{k}\pi}\ell_{it}\right]$ onto the space spanned by the incidental parameters, under a metric given by $\mathbb{E}\left[\partial_{\pi^{2}}\ell_{it}\right]$. Specifically, $\Xi_{it,k}=\lambda_{i,k}^{*\prime}f_{0t}+\lambda_{0i}^{\prime}f_{t,k}^{*}$, where $\left(\lambda_{i,k}^{*},f_{t,k}^{*}\right)$ is defined as the solution to the following optimization problem, \[ \left(\lambda_{i,k}^{*},f_{t,k}^{*}\right)\in\underset{\lambda_{i,k},f_{t,k}}{\operatorname{argmin}}\sum_{i,t}\mathbb{E}\left[\partial_{\pi^{2}}\ell_{it}\right]\left(\frac{\mathbb{E}\left[\partial_{\theta_{k}\pi}\ell_{it}\right]}{\mathbb{E}\left[\partial_{\pi^{2}}\ell_{it}\right]}-\lambda_{i,k}^{\prime}f_{0t}-\lambda_{0i}^{\prime}f_{t,k}\right)^{2} \] In addition, we define the operator $D_{\theta\pi^{q}}\ell_{it}=\partial_{\theta\pi^{q}}\ell_{it}-\Xi_{it}\partial_{\pi^{q+1}}\ell_{it}$ for $q=0,1,2$. Intuitively, these operators remove the component of $\partial_{\theta\pi^q}\ell_{it}$ that can be explained by the individual and time fixed effects. Furthermore, define a $p \times p$ matrix $$ \bar{W}_{NT}:=\frac{1}{NT}\sum_{i=1}^N \sum_{t=1}^T \mathbb{E}\left[\partial_{\theta \theta} \ell_{i t}-\partial_{\pi^2} \ell_{i t} \Xi_{i t} \Xi_{i t}^{\prime}\right] $$ which will be used to characterize the information matrix for $\theta$.
We now introduce the regularity conditions that are required for the asymptotic results.
Assumption (ref) strengthens the conditions used for Theorem (ref).
Assumption (ref)(i) and (ii) require the loss function to exhibit sufficient smoothness and local strong convexity. Assumption (ref)(iii) is a generalized noncollinearity condition, similar to those commonly imposed in the factor literature, such as Assumption A in \citet*{bai_panel_2009} and Assumption 1(vi) in \citet*{chen_nonlinear_2021}). In the linear case, Assumption (ref)(iii) simplifies to requiring $ \sum_{i=1}^N \sum_{t=1}^T \mathbb{E}\left[\left(X_{i t}-\Xi_{i t}\right)\left(X_{i t}-\Xi_{i t}\right)^{\prime}\right]$ to be positive definite. Intuitively, this condition ensures that the covariates exhibit sufficient variation across individuals and over time, thus guaranteeing the identification of $\theta$.
Theorem (ref) examines the numerical convergence properties of $\hat{\theta}^{\left(m+1\right)}$. Specifically, Theorem (ref)(i) guarantees the convergence of the iterative procedure, while Theorem (ref)(ii) indicates that after $O\left(\log\left(NT\right)\right)$ iterations, the iterative estimator $\hat{\theta}^{\left(m+1\right)}$ achieves the desired convergence rate for inference. It converges to the true value $\theta_{0}$ much faster than $\gamma_{NT}$, which is the rate at which the initial estimator $\hat{\theta}$ converges in probability. The upper bound of our contraction parameter, specifically the spectral norm of $\bar{C}^{(0)}$, is crucial in the iterative procedure. The closer $\|\bar{C}^{\left(0\right)}\|$ is to $0$, the faster is the contraction and the numerical convergence of the iterative procedure. Conversely, if $\|\bar{C}^{\left(0\right)}\|$ is close to $1$, the numerical convergence of $\hat{\theta}^{\left(m+1\right)}$ is slow.
The following theorem establishes the asymptotic distribution of the second-step estimator $\hat{\theta}^{\left(m+1\right)}$.
Theorem (ref) indicates that $\hat{\theta}^{\left(m+1\right)}$ contains two asymptotic bias terms associated with $\frac{1}{T}\bar{B}_{\infty}$ and $\frac{1}{N}\bar{D}_{\infty}$, respectively. For convex objective functions, the estimator is asymptotically equivalent to the one obtained from (ref), as studied by \citet*{chen_nonlinear_2021}.
This section describes methods for removing the asymptotic bias of the second‐step estimator. We first outline the analytical bias‐correction procedure and then discuss alternative approaches, including the split‐panel jackknife and the bootstrap.
The analytical correction is constructed using sample analogs of the expressions in Theorem (ref), replacing the true values of $(\theta,\pi)$ with their second‐step estimates. Both analytical bias correction and variance estimation require consistent estimators of the quantities $\bar{B}{\infty}$, $\bar{D}{\infty}$, and $\bar{W}_{\infty}$ defined in Theorem (ref). Let $\hat{B}$, $\hat{D}$, and $\hat{W}$ denote the corresponding sample analogs, obtained by substituting sample averages for expectations and replacing the true parameters with their second‐step estimates. For example, $$\hat{W}=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\partial_{\theta\theta}\hat{\ell}_{it}-\partial_{\pi^{2}}\hat{\ell}_{it}\hat{\Xi}_{it}\hat{\Xi}_{it}^{\prime}$$ where $\partial_{\theta\theta}\hat{\ell}_{it}=\partial_{\theta\theta}\ell_{it}\left(W_{it};\hat{\theta},\hat{\lambda_{i}}^{\prime}\hat{f_{t}}\right)$, $\partial_{\pi^{2}}\hat{\ell}_{it}=\partial_{\pi^{2}}\ell_{it}\left(W_{it};\hat{\theta},\hat{\lambda_{i}}^{\prime}\hat{f_{t}}\right)$, and $\hat{\Xi}_{it}$ is a $p$-dimensional vector with elements $\hat{\Xi}_{it,k}=\lambda_{i,k}^{\#\prime}\hat{f}_{t}+\hat{\lambda}_{i}^{\prime}f_{t,k}^{\#}$ where the pair $\left(\lambda_{i,k}^{\#},f_{t,k}^{\#}\right)$ is defined by $$\left(\lambda_{i,k}^{\#},f_{t,k}^{\#}\right)\in\underset{\lambda_{i,k},f_{t,k}}{\operatorname{argmin}}\sum_{i,t}\left(\partial_{\pi^{2}}\hat{\ell}_{it}\right)\left(\frac{\partial_{\theta_{k}\pi}\hat{\ell}_{it}}{\partial_{\pi^{2}}\hat{\ell}_{it}}-\lambda_{i,k}^{\prime}\hat{f}_{t}-\hat{\lambda}_{i}^{\prime}f_{t,k}\right)^{2}$$ Once those sample analogs are constructed, then the analytic bias correction of $\hat{\theta}$ reads $$\hat{\theta}_{ABC}=\hat{\theta}^{\left(m+1\right)}-\frac{1}{T}\hat{W}^{-1}\hat{B}-\frac{1}{N}\hat{W}^{-1}\hat{D}$$
Although the analytical correction offers a convenient closed‐form adjustment, it requires estimating high‐order derivatives and moment matrices that may be sensitive to model specification and sample size. In practice, alternative approaches—such as the split‐panel jackknife and the bootstrap—can remove the leading bias terms without explicit analytical derivations.
Following \citet*{dhaene_split-panel_2015} and \citet*{chen_nonlinear_2021}, bias correction can also be implemented using the split‐panel jackknife method. Let $\hat{\theta}_{N / 2, T}^1$ and $\hat{\theta}_{N / 2, T}^2$ be the second-step estimators based on the subsamples $\{(i, t): i=1, \ldots, N / 2 ; t=1, \ldots, T\}$ and $\{(i, t): i=N / 2+1, \ldots, N ; t=1, \ldots, T\}$, respectively. Similarly, let $\hat{\theta}_{N, T / 2}^1$ and $\hat{\theta}_{N, T / 2}^2$ be the estimators obtained from the subsamples $\{(i, t): i=1, \ldots, N ; t=1, \ldots, T / 2\}$ and $\{(i, t): i=1, \ldots, N ; t=T / 2+1, \ldots, T\}$, respectively. The jackknife bias‐corrected estimator is \[ \hat{\theta}_{\text{JBC}} = 3 \hat{\theta} - \frac{1}{2}\left(\hat{\theta}_{N / 2, T}^1+\hat{\theta}_{N / 2, T}^2\right) - \frac{1}{2}\left(\hat{\theta}_{N, T / 2}^1+\hat{\theta}_{N, T / 2}^2\right). \] Under suitable homogeneity and stationarity conditions ensuring that the asymptotic biases of all estimators converge to the same limit, \[ \sqrt{N T}\left(\hat{\theta}_{\text{JBC}}-\theta_0\right) \xrightarrow{d} \mathcal{N}\left(0,\bar{W}_{\infty}^{-1} \right) \] Other bias‐correction methods include the leave‐one‐out jackknife of \citet*{hughes_estimating_2022} and the bootstrap procedures of \citet*{higgins_bootstrap_2024} and \citet*{sun_k-step_2010}. However, these approaches are not directly applicable to models with interactive fixed effects. Extending them, particularly adapting the bootstrap of \citet*{higgins_bootstrap_2024} to account for time-specific or interactive effects, remains an interesting avenue for future research.
In this section, we assess the finite sample performance of our proposed approach using Monte Carlo simulations. Section (ref) outlines the main steps of the first-step estimation algorithm. Section (ref) discusses the choice of the tuning parameter $\nu$. Section (ref) presents the simulation designs and results.
This section outlines the computational procedures for the proposed estimation algorithms. We first present the single-index panel model and its implementation; the binary logit specification in Example (ref) is treated as a special case. We then describe the random coefficient logit model in Example (ref). Additional details and formal descriptions can be found in Section (ref).
\paragraph{Single-Index Panel Model} Algorithm (ref) summarizes the first-stage estimation procedure for the class of single-index panel models, which includes the binary logit specification in Example (ref) as a special case. We introduce a slack variable $Z_{\Pi}$ to separate the low-rank interactive effects from the linear index component. This reformulation allows the optimization problem in (ref) to be expressed in an equivalent and computationally convenient form:
To solve the minimization problem ((ref)), we use an Alternating Direction Method of Multipliers (ADMM) algorithm, which relies on the following augmented Lagrangian
Here, $U_{p}$ and $U_{v}$ are the scaled dual variables corresponding to the linear constraints $V=\sum_{j=1}^{p}X_{j}\theta_{j}+Z_{\Pi}$ and $Z_{\Pi}-\Pi=0$, and $\eta>0$ is the penalty parameter for constraint violations. Due to the separability of the parameters in $\mathscr{L}$, the ADMM algorithm proceeds by iteratively minimizing the augmented Lagrangian in blocks with respect to the original variables, in this case $\left(V,\theta,\Pi\right)$ and $Z_{\Pi}$, and then then updating the dual variables $U_{p}$ and $U_{v}$.
As shown in Line (ref) of Algorithm (ref), $\Pi$ is updated via singular value thresholding (SVT). Specifically, recall that the singular value decomposition (SVD) of a matrix $A$ is \[ A=U_{A}\Sigma_{A}V_{A}^{\prime}\in\mathbb{R}^{N\times T},\quad\Sigma_{A}=\operatorname{diag}\left(\left\{ \sigma_{s}\right\} _{s\in\left[N\land T\right]}\right) \] For $\iota>0$, the soft-thresholding operator $S_{\iota}\left(\cdot\right)$ is defined as \[ S_{\iota}\left(A\right)\coloneqq U_{A}S_{\iota}\left(\Sigma_{A}\right)V_{A}^{\prime},\quad S_{\iota}\left(\Sigma_{A}\right)=\operatorname{diag}\left(\left\{ \sigma_{s}-\iota\right\} _{+}\right) \] where $t_{+}$ is the positive part of $t$, namely, $t_{+}=\max\left(0,t\right)$. See Theorem 2.1 in \citet*{cai_singular_2008} for details.
\paragraph{Random Coefficient Logit Panel}
Algorithm (ref) outlines the first-stage estimation algorithm for the binary random coefficient logit panel model in Example (ref) using a majorization-maximization (MM) approach. For simplicity, we assume that all coefficients $\beta$ are random, but incorporating fixed coefficients, such as $\beta_{j,it} = \beta_j$ for some $j$, is straightforward. Our algorithm extends the work of train_discrete_2009 and james_mm_2017 by introducing a novel surrogate function to handle the penalty term. Specifically, to update $\Pi$ at step $k+1$, we use the surrogate function $Q\left(\Pi|\theta^{\left(k\right)},\Pi^{(k)}\right)$, defined as $$
$$ where the closed-form solution is given by $$ \Pi^{(k+1)} \gets S_{NT\nu} \left( \Pi^{(k)} - \left(\sum_{r=1}^R w^{(k)}_{i t r} h\left(Y_{it},X_{it}^{\prime}\beta_{itr}+\pi_{it}^{(k)}\right) \right)_{i,t} \right)$$ with $h(y, v) = \Lambda(v)-y$ and $\Lambda(v) = \frac{1}{1 + e^{-v}}$. This is detailed in Line (ref) of Algorithm (ref).
The selection of the tuning parameter $\nu$ is crucial for the proper implementation of the estimator. The objective is to choose $\nu$ such that it dominates the "score" of our penalized estimator. This ensures a desirable rate of convergence and permits the consistent estimation of the rank of $\Pi_{0}$, even in scenarios where it is not known a priori. Nonetheless, the computation of the "score" is generally infeasible, as it depends on unknown true parameters $\left(\theta_{0},\Pi_{0}\right)$. This challenge becomes more complex in the context of time series data, where the selection of $\nu$ must also account for the degree of temporal dependence.
Several general approaches are commonly used to choose the optimal tuning parameters. As is standard in the statistics and machine learning literature, the most popular data-driven method is cross-validation. Another common approach involves information criteria (IC). Methods like grid search or randomized search are also employed, where a range or a sample of tuning parameters is systematically explored to identify the best value.
In this paper, we select the tuning parameter $\nu$ using a modified information criterion (IC), motivated by prior work on low-rank and grouped panel estimators (e.g., \citet*{bai_determining_2002}; \citet*{belloni_high_2023}; \citet*{su_identifying_2016}) and provides a practical way to balance model complexity and statistical fit.
Specifically, for each candidate value of $\nu$, we compute the corresponding rank estimate $\hat{r}(\nu)$, as defined in Equation (ref),\footnote{The consistency of the rank estimator is established in Corollary (ref).} and evaluate
The first term on the right-hand side of (ref) measures the overall fit to the data, and the second term introduces a penalty proportional to model complexity through the factor $\varrho_{NT}$. We experimented with several alternatives and found that $\varrho_{NT} = \frac{1}{2}\log(N \wedge T)\,\tfrac{N \vee T}{N T}$ works fairly well in Design 1 and so does $\varrho_{NT} = \frac{1}{2}\tfrac{\log\log(\sqrt{NT})}{\sqrt{NT}}$ in Design 2. Accordingly, we adopt these specifications for $\varrho_{NT}$ in the analysis that follows.
\paragraph{Data Generating Process} To evaluate the finite-sample performance of the estimation procedure, we consider two data generating processes (DGPs) that cover models with convex and non-convex objective functions.
Design 1: The data are generated from the following model \[ y_{it}=\mathbb{I}\left\{ X_{it}^{\prime}\theta_{0}+\sum_{s=1}^{2}\lambda_{is}f_{ts}-\varepsilon_{it}\geq0\right\} \] where $\mathbb{I}\left\{ \cdot\right\} $ denotes the indicator function, and the error term $\varepsilon_{it}$ follows a logistic distribution. We set $p=3$ and $r=2$ and take the coefficients $\theta_{0}$ to satisfy $\theta_{0}=\left(1,1,1\right)^{\prime}$. The regressor vector $X_{it}=(x_{it,1},x_{it,2},x_{it,3})^\prime$ consists of three variables, $\lambda_i=(\lambda_{i1},\lambda_{i2})^\prime$, and $f_t=(f_{t1},f_{t2})^\prime$. Specifically, we define $x_{it,k}=\tilde{x}_{it,k}+\rho\left(\lambda{}_{{ik}}^{2}+f_{tk}^{2}\right)$ for $k=1,2$ and $x_{it,3}=\tilde{x}_{it,3}$. The parameter $\rho=0.2$ governs the correlation between the covariates and the fixed effects. Here, $\lambda_i$, $f_t$, and $\tilde{x}_{it}$ are mutually independent across both $i$ and $t$ and $\lambda_{i}\sim \mathcal{N}\left(\ensuremath{\binom{1}{1}},I_{2}\right)$ and $f_{t}\sim \mathcal{N}\left(\ensuremath{\binom{0}{0}},I_{2}\right)$.
We consider two alternative specifications for the distribution of $\tilde{x}_{it}$:
Design 2: We extend Design 1 by allowing random coefficients $\beta_{it}$. With a slight abuse of notation, \[ y_{it}=\mathbb{I}\left\{ X_{it}^{\prime}\beta_{0,it}+\sum_{s=1}^{2}\lambda_{is}f_{ts}-\varepsilon_{it}\geq0\right\} \] The covariates $X_{it}$, factor loadings $\lambda_{is}$, and factors $f_{ts}$ follow the same distributions as in Design 1, and we adopt the same specifications of $\tilde{x}_{it}$. The coefficients $\beta_{0,it}\sim \mathcal{N}\left(\bar{\beta}_{0},\Sigma_{0}\right)$ with $\bar{\beta}_{0}=\left(1,1,1\right)^{\prime}$ and $\Sigma_{0}=0.3 I_{3}$.
\paragraph{Performance} For each model design, we consider different combinations of sample sizes $(N,T) \in \{(100,100),\,(150,100),\,(200,100),\,(150,150),\,(200,150),\,(200,200)\}$ to examine how the first-step and second-step estimators behave as the sample size increases. Each experiment is replicated $100$ times. The ADMM algorithm is used to estimate the model in Design 1 and Design 2, while the MM algorithm is applied to estimate the model in Design 3. Detailed descriptions of the implementation are provided in Section (ref).
To evaluate the performance of my estimator, we report the root mean squared error (RMSE) for both steps, defined as $\operatorname{RMSE}(\hat{\theta})=\sqrt{S^{-1}\sum_{s=1}^{S}|\hat{\theta}^{(s)}-\theta_0|^2}/\|\theta_0\|$, and the average estimated rank, $\mathbb{E}[\hat{r}]=S^{-1}\sum_{s=1}^{S}\hat{r}_s$. We set $S=100$ Monte Carlo replications. These measures summarize the finite-sample accuracy of the estimators and the typical rank selected across simulations.
Table (ref) summarizes the finite-sample performance of the proposed estimators. Several key observations can be made from the results. (i) At each sample size, the second-step estimators exhibit smaller RMSE compared to the first-step estimators; (ii) As the sample size increases, particularly when both $N$ and $T$ grow, the RMSE of the first-step estimators tends to decrease; (iii) Similarly, the RMSE of the second-step estimators generally decreases with increasing sample size; (iv) The rank estimator correctly estimates the rank as the sample size increases, with the average rank stabilizing at $\hat{r} = 2$ for larger samples. These observations confirm the consistency results for both the first-step and the second-step estimators. For example, the RMSE decreases from $3.06\%$ when $N=T=100$ to $1.57\%$ at $N=T=200$, consistent with the theoretical prediction that estimator performance improves with larger sample sizes.
Table (ref) summarizes the results for Design 2, which incorporates random coefficients. The findings parallel those of Design 1: the second-step estimator consistently improves upon the initial estimator, and both bias and RMSE decline with sample size. In this case, we focus on the accuracy of the estimated mean of the random coefficients, which is reliably recovered in larger samples.
In general, the simulations confirm the main theoretical insights. The two-step procedure delivers substantial efficiency gains, the estimators are consistent as both dimensions of the panel grow, and rank selection is reliable in sufficiently large samples.
In this section, we apply our methodology to revisit the joint financing behavior of US non-financial corporations, focusing on how firms adjust debt issuance and equity repurchase decisions in response to relative valuations across capital markets. The objective is to illustrate the empirical applicability of our approach by extending the baseline model of ma_nonfinancial_2019 and to compare its benchmark fixed-effects logit with our proposed specification.
Empirical studies suggest that firms often act as cross-market arbitrageurs, substituting between debt and equity to exploit differences in financing costs. When debt is relatively inexpensive, firms tend to issue debt and repurchase equity, whereas undervalued equity prompts them to issue equity and retire debt. In a recent study, ma_nonfinancial_2019 documents that a substantial share of financing flows originate from non-financial firms that simultaneously issue in one market while repurchasing in another. About 45 percent of quarterly net equity repurchases (by value) occur in periods when firms are concurrently net issuers of debt, and roughly 50 percent of seasoned equity issuance takes place among firms that are net retiring debt. These coordinated patterns indicate that firms actively rebalance their capital structures to exploit valuation differentials across markets, with this behavior most pronounced among large, profitable, and financially unconstrained firms.
To quantify this relationship, the study estimates the probability that a firm simultaneously issues equity and retires debt using a panel logit model with firm fixed effects:
where $\Lambda(\cdot)$ denotes the logistic cumulative distribution function. The binary outcome variable $y_{it}$ equals one if the firm issues equity (negative net equity repurchases) and retires debt (negative net debt issuance) in the same quarter, and zero otherwise. Net equity repurchases are measured as purchases minus sales of common and preferred stock (PRSTKC$-$SSTK), and net debt issuance is defined as long-term debt issuance minus long-term debt reduction (DLTIS$-$DLTR). The vectors $X_{D,it-1}$ and $X_{E,it-1}$ capture debt and equity market valuations measured at the end of the quarter $t-1$, while $Z_{it}$ includes additional firm-level controls that can affect financing behavior.
The model uses three firm-level indicators of market valuation. Two bonds-based measures, the credit spread and the term spread, capture relative valuations in the debt market ($X_{D,it-1}$). The firm-level credit spread is computed as the face-value-weighted average yield differential between the firm’s bonds and the nearest-maturity Treasury, while the term spread is the yield difference between the nearest-maturity Treasury and the three-month Treasury bill.
The valuation of the equity market ($X_{E,it-1}$) is proxied by the firm's specific value-price ratio (V / P), where $V$ denotes the intrinsic value of the equity estimated from the residual income model and $P$ represents the market price of equity dong_overvalued_2012. Control variables include net income, cash holdings, capital expenditures, deviations from target leverage, asset growth, and firm size.\footnote{See ma_nonfinancial_2019 for detailed definitions and data construction.}
Before introducing our model, we first examine whether the data exhibit latent time-varying structure by estimating a two-way fixed-effects (TWFE) specification following \citet*{fernandez-val_individual_2016} and analyzing its residuals, $$\widehat{u}_{TWFE,i t}=y_{i t}-\Lambda\left(\hat{\lambda}_i +\hat{f}_t+\hat{\beta}_{D}^{\prime} X_{D, i, t-1}+\hat{\beta}_{E}^{\prime} X_{E, i, t-1}+\hat{\gamma}^{\prime} Z_{i t}\right)$$ Figure (ref) shows systematic temporal and cross-sectional patterns: firms display persistent deviations across quarters (horizontal streaks), and many firms experience synchronized movements within the same periods (vertical color bands). These residual patterns point to time-varying unobserved heterogeneity.
We begin with a brief overview of the data and the model. Following ma_nonfinancial_2019, this analysis integrates several primary data sources, combining firm-level information on bond prices, equity valuations, and accounting fundamentals with aggregate market indicators. Bond data is obtained primarily from the Trade Reporting and Compliance Engine (TRACE) and supplemented by Datastream and Mergent’s Fixed Income Securities Database (FISD). Equity returns and valuation measures come from CRSP and the Institutional Brokers’ Estimate System (IBES), while balance sheet and cash flow variables are drawn from Compustat.\footnote{Access to TRACE, Compustat, CRSP, and IBES is provided through Wharton Research Data Services (WRDS), while Datastream is accessed via the LESG portal.} The sample covers the period 2003-2024, producing a panel of $N = 212$ firms observed in $T = 88$ quarters.\footnote{Additional details on data sources, variable definitions, and sample construction are provided in Appendix (ref).} Summary statistics for key variables appear in Table (ref).
We extend this framework by incorporating interactive fixed effects and random coefficients, allowing for unobserved common shocks and heterogeneous firm responses to market valuations. Formally, the latent index is given by
where $f_t$ denotes latent factors that capture aggregate market shocks, $\lambda_i$ are the corresponding firm-specific factor loadings, and $\beta_{D,it}$ and $\beta_{E,it}$ represent random coefficients reflecting heterogeneity in sensitivities to debt and equity valuations. Specifically, we assume that $
\sim \mathcal{N}\!\left(
, \, \Sigma_{\beta} \right)$. We employ the same set of market valuation and control variables as in (ref).
We report our estimation results in Table (ref). Model A corresponds to the conventional panel logit specification with individual fixed effects, as defined in (ref), and serves as the benchmark. It is estimated using the conditional maximum likelihood approach following ma_nonfinancial_2019. In contrast, Model B adopts the specification in (ref), which extends the benchmark by allowing for random coefficients on lagged valuation measures and by incorporating interactive fixed effects. The rank of the interactive component is estimated using the procedure described in (ref); the first-step estimation yields an estimated rank of $\hat{r} = 2$. Accordingly, Table (ref) reports the second-step estimates, including the means and standard deviations of the random coefficients.
In general, the results re-confirm the importance of market valuations in shaping joint financing decisions of firms and are consistent with economic intuition: firms are more likely to simultaneously issue equity and retire debt when the cost of debt is high (e.g., elevated bond spreads or expected excess returns) and the cost of equity is low (e.g., a low V/P ratio or low expected excess equity returns).
We find that the signs of the estimated effects are generally robust to the inclusion of latent factors, with variations in their magnitudes across specifications. For instance, the effect of lagged term spread is larger in our model than in the benchmark specification: a one-standard-deviation decrease in term spread multiplies the odds of joint equity issuance and debt retirement by about 1.33 on average across firms, compared with roughly 1.25 in Model A. We also uncover substantial heterogeneity across firms in their sensitivities to valuation measures. In particular, the estimated standard deviation of the random coefficient on the lagged credit spread is 0.363 and statistically significant, indicating meaningful variation in firms' responses to credit market conditions.
We regard the comparison between our results and those from the benchmark model as an illustration that highlights the empirical value of allowing richer forms of unobserved heterogeneity. Firms adjust their financing decisions in response to evolving macrofinancial conditions and to unobserved, time-varying shocks—such as shifts in credit availability, monetary policy, or market liquidity—that jointly influence financing behavior and asset valuations. At the same time, firms differ in the extent to which they adjust their financing in response to observed market valuations. By incorporating interactive fixed effects and random coefficients, the proposed approach flexibly captures both common dynamics and firm-specific heterogeneity, demonstrating its applicability to complex empirical settings.
This paper explores nonlinear panel data models with interactive fixed effects and tackles nonconvex objective functions in large $N$ and $T$ settings. We propose a two-step estimation method using nuclear norm regularization and an iterative process. First, we obtain preliminary estimates of the slope coefficient, factors, and factor loadings using NNR. The second step refines these estimates through an iterative process. We establish the consistency of the preliminary estimator and derive the asymptotic distribution for the second-step estimator. The performance of the estimator is demonstrated through Monte Carlo simulations and an empirical application to a random coefficient binary logit model.
Several extensions remain open for future research. First, it would be valuable to explore whether post-NNR inference can be achieved within a finite number of iterations, extending the work of chernozhukov_inference_2023 and choi_inference_2024 to nonlinear panel settings. Second, while this paper accounts for certain heterogeneity in slope coefficients within the random coefficient logit panel model, allowing for even richer forms of coefficient heterogeneity, as discussed by bonhomme_fixed_2024, presents an interesting avenue for further study. These topics are left for future exploration.
\appendixtitleon \appendixtitletocon