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.
70,264 characters · 9 sections · 61 citation commands
Debiased Inference for Dynamic Nonlinear Panels with Multi-dimensional Heterogeneities
Panel data are common in empirical research. Panel data models with two-way fixed effects are widely used to control for unobserved heterogeneity that may correlate with covariates. Individual effects capture time-invariant, individual-specific characteristics, while time effects account for shifts that are common across individuals but vary over time. Traditionally, these models address level heterogeneity by incorporating fixed effects only in the intercept, adjusting for baseline outcome differences across individuals and time periods. However, both economic theory and empirical evidence indicate substantial variation in how outcomes respond to covariates across individuals and over time bc2007. Cross-sectional response heterogeneity reflects individual-specific factors such as preferences, productivity, and access to resources heckman2001, whereas time-varying response heterogeneity captures changes in behavior driven by evolving circumstances---including policy reforms, business cycles, and technological progress ow2021. When such response heterogeneities are correlated with covariates, failing to account for them can lead to biased estimates and invalid inference sc2013.
To illustrate, consider the classic problem of estimating the impact of children on a mother's labor force participation. Research indicates that mothers with more children may differ systematically from those with fewer or none br1988, ae1998, kleven2019. First, mothers who choose to have more children might be inherently less inclined to work. Second, mothers with more children may be those whose labor supply is less sensitive to family size, perhaps due to better financial resources or childcare access. Moreover, over time, declining fertility and rising female labor force participation suggest a shifting baseline, while advances in home technology and welfare reforms may lessen the burden of additional children, enhancing labor supply responsiveness. In a panel dataset tracking households over time, we can control for the first type of heterogeneity---level differences in baseline participation---by including fixed effects in the intercept. However, to address the second type of heterogeneity---variations in response to children across individuals and over time---we must include fixed effects in the slope. Given that labor force participation is a binary outcome and typically exhibits serial correlation, this example underscores the need for a dynamic nonlinear panel data model with fixed effects in both the intercept and slope.
In this paper, we introduce a general class of dynamic nonlinear heterogeneous-parameter (DN-HP) models that incorporate individual and time fixed effects in both the intercept and the slope. This flexible framework accommodates multi-dimensional heterogeneities and encompasses many commonly used panel-data models, including both static and dynamic specifications for linear and limited-dependent-variable outcomes. It nests individual-specific slope models that capture heterogeneous responses across individuals, as well as time-varying-coefficient models that allow response parameters to evolve over time due to shifts in the underlying economic or structural environment. Within this class of models, we focus on estimation and inference for average (common) slope coefficients under a large-$N,T$ framework. In many empirical applications, slope coefficients are primary objects of interest because they summarize economically meaningful responsiveness---often in elasticity or semi-elasticity form---and are the quantities most directly reported and compared across studies. For example, in international trade, gravity models interpret coefficients on distance and policy variables as elasticities or percentage effects on bilateral flows using long panels of trading partners silva2006. In innovation economics, nonlinear count models for patenting rely on slope coefficients to measure the responsiveness of innovative activity to market structure or policy incentives aghion2005. In differentiated product demand estimation, discrete-choice models focus on slope coefficients that capture consumers' price sensitivity and valuation of product characteristics, and are often estimated using product--market panel data with a large number of products or markets observed repeatedly over an extended period of time nevo2001.
In such settings, DN-HP models are generally subject to the incidental-parameter problem: the limiting distribution of the estimator is not centered at zero, and standard test statistics (e.g., Wald, LM, and LR) no longer follow their conventional asymptotic $\chi^{2}$ distributions. To address this, we propose an analytical bias-correction procedure that restores valid large-sample inference for both parameter estimates and test statistics. As a preview, Figure (ref) presents simulated boxplots of the estimation bias in a two-way heterogeneous parameter logit model, comparing estimators from our bias correction procedure to the uncorrected ones. The results reveal that our bias correction procedure significantly reduces the estimation bias while maintaining comparable mean squared errors.
Literature Review. The estimation and inference of fixed-effects models in the presence of the incidental-parameter problem have been extensively studied. Early work, such as c1980, ab1991, and l2002, focused on frameworks with only individual effects in the intercept, establishing fixed-$T$ consistency (short panels) for structural parameter estimators in specific models. In recent decades, the availability of long panel datasets has motivated a large--$N,T$ framework, where the cross-sectional size $N$ and the time dimension $T$ grow at similar rates under rectangular asymptotics. In this setting, estimators remain consistent but exhibit a non-negligible bias of order $O\left( 1/T\right) $ because each individual effect is estimated from only $T$ observations. This incidental-parameter bias is especially pronounced in nonlinear or dynamic models. To address it, researchers have developed a variety of bias correction techniques within a maximum likelihood framework. hn2004 and hk2011 propose \textquotedblleft parameter-based" corrections that remove the bias from the likelihood estimator; w2002 and llw2003 propose \textquotedblleft score-based" methods that modify the profiled score; while \textquotedblleft likelihood-based" approaches such as bh2009 and ah2016 adjust the log-likelihood directly. Most of these procedures are analytical, relying on closed-form approximations to the bias. Alternatively, numerical corrections estimate the bias through resampling or re-estimation, including the jackknife dj2015, the bootstrap ks2016,bkss2020,hj2024, and integrated-likelihood methods ab2009.
For two-way fixed effects models, incorporating time effects in the intercept generates an additional bias of order $O\left( 1/N\right) $ on top of the $O\left( 1/T\right) $ term from individual effects. Several studies extend bias-correction methods to this setting. mw2015a develop a parameter-based analytical correction for dynamic linear models with interactive effects. For nonlinear models, fw2016 propose parameter-based techniques for additive fixed effects, while cfw2014 study interactive effects in static models. Alternatively, ko2018 provide a likelihood-based approach for static models with an arbitrary but known fixed-effect structure. Other model-specific contributions include b2009 and c2017. For a comprehensive overview of these developments, see fw2017. Despite this progress, the two-way literature primarily addresses heterogeneity in levels---that is, intercept effects---rather than heterogeneity in responses (slopes). Extending bias-correction methods to settings with slope heterogeneity remains largely unexplored and is the focus of our contribution.
A substantial related literature explores heterogeneity in slope coefficients within panel models, emphasizing individual-specific parameters robertson1992, ps1995 and time-varying coefficients robinson1989,sw1996, most often in linear or static contexts h2014, pesaran2015. Recent research in dynamic nonlinear panels has advanced these ideas further. fl2013 study linear and nonlinear panel data models with individual-heterogeneous coefficients and endogenous regressors, and propose a bias correction method based on the generalized method of moments. See also fglv2025 for a panel distribution regression with individual-heterogeneous coefficient. cfhn2013 derive partial identification with uniform inference for nonseparable panels under time-homogeneity. bc2014 provide mixture-based point identification for dynamic binary models with maximal cross-sectional heterogeneity, while bc2010 analyze short-panel dynamic binary models with heterogeneous transition probabilities and propose a mean-integrated-MSE (MIMSE) estimator that better balances bias--variance trade-offs in small--$T$ settings than analytical bias correction. While these studies focus on identification or fixed--$T$ settings within semiparametric or nonparametric frameworks, we adopt a parametric, likelihood-based approach that achieves point identification and enables bias-corrected inference for dynamic nonlinear models with two-way slope and intercept effects under large--$N,T$ asymptotics.
Two recent studies are closely aligned with our framework. kn2020 examine linear models with two-way slope and intercept effects, employing an iterative “mean-observation OLS”\ estimator to address bias. Their analysis of U.S. agricultural data reveals pronounced regional differences in heat-yield sensitivity, alongside temporal adaptation as farmers adopt new technologies over time. ls2023 study a similar two-way heterogeneous-slope specification and implement a parameter-based jackknife estimator that enables uniform inference for structural parameters. Their cross-country analysis of the Feldstein--Horioka relation uncovers variation across nations and periods, driven by financial integration and evolving policy regimes. Both studies underscore the importance of slope heterogeneity across individuals and time but remain confined to linear settings. Our study advances this literature by extending the framework to a dynamic nonlinear context and developing a likelihood-based analytical bias-correction procedure that ensures valid estimation and inference. To our knowledge, neither our DN-HP model nor its associated bias correction has been previously explored. Together, the model and method provide a unified approach to estimation and inference in dynamic panels with multi-dimensional heterogeneity.
This paper. We employ the maximum likelihood framework and focus on additive multi-dimensional two-way fixed effects, whose number grows with the sample size. We consider an arbitrary DN-HP model whose log-likelihood function is specified up to the unknown parameters. We construct a modified, or corrected, log-likelihood function by adding two bias correction terms to the original one. These two bias terms are analytically derived by combining i) the Taylor expansion of the original log-likelihood function, in terms of the fixed effects, and ii) an asymptotic expansion of the fixed effects themselves. We show that the corrected log-likelihood function does not suffer from the incidental parameter problem, under appropriate regularity conditions. Our method can be viewed as an extension of ah2016 and bh2009 to the two-way DN-HP models.
One key benefit of our procedure is that we achieve bias corrections for the point estimators and the test statistics by modifying a single object, the log-likelihood function. We rigorously show that the estimators obtained by maximizing the corrected log-likelihood function retain a limiting normal distribution with zero mean, consistent with the asymptotic theory of the classical maximum likelihood estimators (MLE) in the absence of incidental parameters. Beyond the estimators, we also show that the likelihood ratio (LR), Lagrange-multiplier (LM), and Wald statistics, derived from the same corrected log-likelihood function, are asymptotically equivalent, sharing the same asymptotic $\chi^{2}$ distribution under the null hypothesis.
We demonstrate the finite-sample performance of our method through Monte Carlo simulations, comparing it to the original likelihood and potential bias correction devices of the bootstrap and the jackknife. We show that our correction procedure reduces the bias significantly, without increasing the root mean squared errors, and restores the test sizes of the LR test to its nominal level. We find that our procedure outperforms the jackknife when it comes to heterogeneous slope coefficients. We also find that the LR test based on our approach has a stronger power near the true hypothesis than the LR test based on the bootstrap. In the empirical application, we apply our approach with a probit model to study the aforementioned question of how the number of children affects a mother's labor force participation, with a focus on single-mother households. Our analysis reveals that this impact varies significantly across individuals and over time. Neglecting such variation can lead to biased estimates and markedly different conclusions, highlighting the need to address both level and response heterogeneities, as well as the value of our framework in empirical research.
The remainder of the paper is organized as follows. Section (ref) explains our settings, the identification restriction, and the incidental parameter problem of DN-HP models. Section (ref) describes our bias correction procedure and provides relevant statistical properties. Section (ref) presents simulation studies to demonstrate the performance of the method. Section (ref) applies our method to study single mothers' labor force participation. Finally in Section (ref), we leave some closing remarks. Proofs, technical details, additional results, elaborated discussions, etc. are provided in the appendix. These items have references starting with letters.
Notation. We denote $\mathbb{I}_{n}$ to be the $n\times n$ identity matrix and $\iota_{n}$ the $n\times1$ vector of ones. $\otimes$ denotes the Kronecker product and $ \mathds{1} \{\cdot\}$ is the indicator function. Next, for a sequence of square matrices $A_{i}$, $i=1,\ldots,n$, $\operatorname*{diag}\{A_{1},\ldots,A_{n}\}$ represents a block diagonal matrix where each $A_{i}$ is the $i$-th diagonal element. For a vector $v$ and a function $f(v)$, we write $\partial_{v}$ and $\partial_{vv^{\prime}}$ to represent the first and second partial derivatives of $f(v)$ with respect to (w.r.t.) $v$.
In this section, we explain our models and the estimation procedure, leaving two motivating examples in Section (ref) of the supplementary appendix. Let $\{(Y_{it},X_{it}^{\prime})^{\prime}:i=1,\ldots,N;t=1,\ldots,T\}$ be a panel data set of $N$ individuals and $T$ time periods, where, $i$ represents the individual and $t$ represents the time period. Here $Y_{it}$ is a scalar response variable, and $X_{it}$ is a $K\times1$ vector of regressors, with $K$ known and fixed. For each individual $i$, the response $Y_{it}$ is generated sequentially over $t$, by
where $\theta_{0}$ is a $K\times1$ vector of (non-random) unknown structural parameters of interests, $\phi^{0}:=(\alpha_{1}^{0\prime},\ldots,\alpha _{N}^{0\prime},\gamma_{1}^{0\prime},\ldots,\gamma_{T}^{0\prime})^{\prime}$ consists of the $K\times1$ random vectors $\alpha_{i}^{0}$ and $\gamma_{t} ^{0}$, and $f(y|\cdot)$ is a density known up to $\theta_{0}$, $\alpha_{i} ^{0}$, and $\gamma_{t}^{0}$. We restrict our attention to the case where $f\left( \cdot\right) $ depends on $\theta_{0}$, $\alpha_{i}^{0}$, and $\gamma_{t}^{0}$ through the additive structure $\theta_{0}+\alpha_{i} ^{0}+\gamma_{t}^{0}$, the \textquotedblleft heterogeneous parameter". The vectors $\alpha_{i}^{0}$ and $\gamma_{t}^{0}$ capture the individual-specific and time-specific heterogeneities, respectively. Generally, $\alpha_{i}^{0}$ and $\gamma_{t}^{0}$ may be correlated with the regressor $X_{it}$. Throughout the paper, we employ the \textquotedblleft fixed-effect" framework to condition on the realizations of $\alpha_{i}^{0}$ and $\gamma_{t}^{0}$ (denoted by the same symbols) in the density, and treat them as unknown nuisance parameters to be estimated together with $\theta_{0}$. Conditioning on $\alpha_{i}^{0}$ and $\gamma_{t}^{0}$, we assume $(Y_{it},X_{it}^{\prime })^{\prime}$ independent across $i$ but serially dependent over $t$. In particular, we allow $X_{it}$ to contain both strictly exogenous components and predetermined components w.r.t. $Y_{it}$. For instance, $X_{it}$ may contain lagged values of $Y_{it}$. In this paper, we focus on the case where all parameters in ((ref)) are heterogeneous. This is without loss of generality, because a \textquotedblleft homogeneous parameter" can be accommodated by setting the relevant component of $\alpha_{i}^{0}+\gamma _{t}^{0}$ to $0$ for all $\left( i,t\right) $.
We are interested in the estimation of and the inference about the structural parameter $\theta_{0}$ in ((ref)), under the presence of the nuisance parameter $\phi^{0}$. Our approach is likelihood-based. Denote $\phi:=(\alpha_{1}^{\prime},\ldots,\alpha_{N}^{\prime},\gamma_{1}^{\prime },\ldots,\gamma_{T}^{\prime})^{\prime}$ and define the log-likelihood function (\textquotedblleft likelihood" hereafter) as
Due to the additive structure $\theta+\alpha_{i}+\gamma_{t}$, $\theta$ is not identified because the likelihood function is invariant to transformations
for all pairs of constant $(c_{1},c_{2})$ and all $(i,t)$. To resolve this, we follow ls2023 to reparameterize the likelihood by setting $\alpha _{N}=-\sum_{i=1}^{N-1}\alpha_{i}$ and $\gamma_{T}=-\sum_{t=1}^{T-1}\gamma_{t} $, thereby imposing the normalization $\sum_{i=1}^{N}\alpha_{i}=0$ and $\sum_{t=1}^{T}\gamma_{t}=0$. Any \textquotedblleft average effect" is absorbed by $\theta$. Let $\psi:=\left( \alpha_{1}^{\prime},\ldots ,\alpha_{N-1}^{\prime},\gamma_{1}^{\prime},\ldots,\gamma_{T-1}^{\prime }\right) ^{\prime}$, the reparameterized likelihood can be constructed from the original $\ell\left( \theta,\phi\right) $ as \[ \ell\left( \theta,D^{\prime}\psi\right) =:l(\theta,\psi),\qquad D:=\operatorname*{diag}\{D_{1},D_{2}\}, \] where $D_{1}:=(\mathbb{I}_{N-1},-\iota_{N-1})\otimes\mathbb{I}_{K}$ and $D_{2}:=(\mathbb{I}_{T-1},-\iota_{T-1})\otimes\mathbb{I}_{K}$. This rules out the indeterminacy in ((ref)) without explicitly using the Lagrangian.
The dynamic nonlinear heterogeneous parameter, DN-HP, models generally suffer from the incidental parameter problem (IPP). We briefly explain this problem here, relegating to Section (ref) of the supplementary appendix i) the exact characterization of the remainders, ii) an illustrative example, etc. We use $\mathbb{E}$ to denote the expectation w.r.t. $\prod_{i=1}^{N}\prod_{t=1}^{T}f(y_{it}|X_{it},\theta_{0}+\alpha _{i}^{0}+\gamma_{t}^{0})$, which is the true density evaluated at $\theta_{0} $, given: i) initial values of the predetermined regressors, ii) all strictly exogenous regressors, and iii) the unobserved effect $\phi^{0}$. Following the profiled likelihood framework of ps2006, we continue our discussion by defining, for a given $\theta$ and each pair of $(N,T)$,
Here $\psi(\theta)$ is referred to as the pseudo-true value and $l(\theta ,\psi(\theta))$ is not subject to the IPP. See, e.g., hn2004 and ah2016. In practice, however, $\psi(\theta)$ is generally infeasible and the plug-in version $l(\theta,\widehat{\psi}(\theta))$ is used instead. Specifically, the MLE of $\theta_{0}$, $\widehat{\theta}:=\arg\max_{\theta }\ell(\theta,\widehat{\psi}(\theta))$. The LR and LM test statistics (for hypothesis about $\theta_{0}$) are also constructed using $l(\theta ,\widehat{\psi}(\theta))$. A fundamental problem is that $l(\theta ,\widehat{\psi}(\theta))$ contains a non-negligible asymptotic bias of order $O(1/\sqrt{NT})$, due to the estimation errors in $\widehat{\psi}(\theta)$. This is referred to as the IPP in the context of the DN-HP model. The consequence is that, $\widehat{\theta}$ and many likelihood-based test statistics deviate from their standard asymptotic distribution, leading to incorrect inferences.
To demonstrate the non-negligible bias in $l(\theta,\widehat{\psi}(\theta))$, consider a Taylor expansion of $l(\theta,\widehat{\psi}(\theta))$ around $\psi(\theta)$:
where $r^{l}(\theta)$ is a remainder, $\widehat{\psi}(\theta)-\psi(\theta)$ represents the estimation errors of the fixed effects,
are the score and Hessian. Here the estimation errors $\widehat{\psi} (\theta)-\psi(\theta)$ can be derived by another Taylor expansion. Seeing $s(\theta,\widehat{\psi}(\theta))=0$ as $\widehat{\psi}(\theta)$ is the maximizer, a Taylor expansion of $s(\theta,\widehat{\psi}(\theta))$ around $\psi(\theta)$ gives
where $r^{s}(\theta)$ is a remainder. Combining ((ref)) and ((ref)), and taking expectation, we obtain
where $r(\theta)$ is the final remainder depending on $r^{l}(\theta)$ and $r^{s}(\theta)$ satisfying $\mathbb{E}r(\theta)=o(1/\sqrt{NT})$, uniformly in $\theta$, as $N,T\rightarrow\infty$ with $N/T\rightarrow\kappa$ for some $0<\kappa<\infty$. This implies that $\sqrt{NT}\mathbb{E}r(\theta)=o(1)$ and is negligible. $B(\theta)$ is the IPP bias arising from the estimation errors of $\widehat{\psi}(\theta)$. As opposite to $\mathbb{E}r(\theta)$, the bias $B(\theta)$ is not negligible, in the sense that $\sqrt{NT}B(\theta )\not \rightarrow 0$ as $N,T\rightarrow\infty$ with $N/T\rightarrow\kappa$. Consequently, $\widehat{\theta}$ and many likelihood-based test statistics inherit the bias from the likelihood, leading to invalid inferences.
The bias expansion of $\widehat{l}(\theta):=l(\theta,\widehat{\psi}(\theta))$ in Equation ((ref)) motivates our approach of correcting the likelihood $l(\theta,\widehat{\psi}(\theta))$ by adding $B(\theta)$. Our corrected likelihood can be viewed as an approximation to the IPP-free infeasible likelihood $l(\theta):=l(\theta,\psi(\theta))$, with approximation error $o_{\mathbb{P}}(1/\sqrt{NT})$. That is, it does not suffer from the IPP to the first order. Intuitively, the estimator of $\theta_{0}$, and the LR, LM, and Wald test statistics, will follow their standard asymptotic distributions when they are derived from the corrected likelihood (instead of $\widehat{l}(\theta)$). In what follows, we explain our correction procedure and give the main result. We leave in Section (ref) of the supplementary appendix i) more technical details; ii) formal statements and discussions about our secondary results (Equations (ref), (ref), (ref), etc.); iii) useful remarks; etc. In addition, we explain the relation of our method to bh2009 and ah2016 in Section (ref) of the supplementary appendix.
In this paper, we decompose $B(\theta)=B_{\alpha}(\theta)+B_{\gamma} (\theta)+o(1/\sqrt{NT})$ to obtain the bias expansion
where \[ B_{\alpha}(\theta):=\frac{1}{2}\operatorname*{trace}[S_{\alpha\alpha} (\theta)H_{\alpha\alpha}^{\mathcal{\ast}}(\theta)],\qquad B_{\gamma} (\theta):=\frac{1}{2}\operatorname*{trace}\mathbb{[}S_{\gamma\gamma} (\theta)H_{\gamma\gamma}^{\ast}(\theta)] \] appear because of the estimation errors of, respectively, the individual-specific heterogeneities $\alpha_{i}^{0}$ and the time-specific heterogeneities $\gamma_{t}^{0}$.\footnote{In addition, we have $B_{\alpha}(\theta)=O\left( T^{-1}\right) $ and $B_{\gamma}(\theta)=O\left( N^{-1}\right) $, indicating that $\sqrt {NT}(\mathbb{E}l(\theta)-\mathbb{E}\widehat{l}(\theta))$ is proportional to $\sqrt{\kappa}+\sqrt{\kappa^{-1}}$ as $N,T\rightarrow\infty$ provided that $N/T\rightarrow\kappa$. This is similar to, e.g., fw2016.} The terms $H_{\alpha\alpha}^{\mathcal{\ast}}(\theta)$, $H_{\gamma\gamma}^{\ast}(\theta )$, $S_{\alpha\alpha}(\theta)$, and $S_{\gamma\gamma}(\theta)$ are defined from partitioning\footnote{Note that $\mathbb{E}s\left( \theta\right) =0$ by definition. We keep $\mathbb{E}s\left( \theta\right) $ here for the convenience of constructing the corrected likelihood.} $[\mathbb{E} H(\theta)]^{-1}$ and $s(\theta)-\mathbb{E}s\left( \theta\right) $ as
Here $H_{\alpha\alpha}^{\mathcal{\ast}}(\theta)$ is $(N-1)K\times(N-1)K$, which corresponds to the inverse expected Hessian w.r.t. $\left( \alpha _{1}^{\prime},\ldots,\alpha_{N-1}^{\prime}\right) ^{\prime}$, and $H_{\gamma\gamma}^{\ast}(\theta)$ is $(T-1)K\times(T-1)K$, corresponding to the inverse expected Hessian w.r.t. $\left( \gamma_{1}^{\prime},\ldots ,\gamma_{T-1}^{\prime}\right) ^{\prime}$. The off-diagonal blocks $H_{\alpha\gamma}^{\mathcal{\ast}}(\theta)$ and $H_{\gamma\alpha }^{\mathcal{\ast}}(\theta)$ are defined accordingly (but are irrelevant to the construction of the corrected likelihood). Similarly, $\widetilde{s}_{\alpha }(\theta)$, the score w.r.t. $\left( \alpha_{1}^{\prime},\ldots,\alpha _{N-1}^{\prime}\right) ^{\prime}$, is $(N-1)K\times1$ and $\widetilde{s} _{\gamma}(\theta)$, the score w.r.t. $\left( \gamma_{1}^{\prime} ,\ldots,\gamma_{T-1}^{\prime}\right) ^{\prime}$ is $(T-1)K\times1$. $\widetilde{s}_{\alpha}(\theta)$ and $\widetilde{s}_{\gamma}(\theta)$ are used to construct \[ S_{\alpha\alpha}(\theta):=\mathbb{E\{}\widetilde{s}_{\alpha}(\theta )\mathbb{[}\widetilde{s}_{\alpha}(\theta)]^{\prime}\},\qquad S_{\gamma\gamma }(\theta):=\mathbb{E\{}\widetilde{s}_{\gamma}(\theta)\mathbb{[}\widetilde{s} _{\gamma}(\theta)]^{\prime}\}, \] which are essentially the covariance matrices of the scores $\widetilde{s} _{\alpha}(\theta)$ and $\widetilde{s}_{\gamma}(\theta)$, respectively.
Estimating $B_{\alpha}(\theta)$ and $B_{\gamma}(\theta)$ by plug-in estimates, the corrected likelihood is established as
Here $\widehat{H}_{\alpha\alpha}^{\mathcal{\ast}}(\theta)$, $\widehat{H} _{\gamma\gamma}^{\ast}(\theta)$, $\widehat{S}_{\alpha\alpha}(\theta)$, and $\widehat{S}_{\gamma\gamma}(\theta)$ are the plug-in versions of the corresponding quantities in $B_{\alpha}(\theta)$ and $B_{\gamma}(\theta )$.\ They are constructed as follows. First, denoting $\widehat{H} (\theta):=H(\theta,\widehat{\psi}(\theta))$, the two Hessian terms $\widehat{H}_{\alpha\alpha}^{\mathcal{\ast}}(\theta)$ and $\widehat{H} _{\gamma\gamma}^{\ast}(\theta)$ are defined from the partition of $[\widehat{H}(\theta)]^{-1}$ as
in the same manner as in ((ref)). Here $\widehat{H}_{\alpha\alpha}^{\mathcal{\ast}}(\theta)$ and $\widehat{H} _{\gamma\gamma}^{\mathcal{\ast}}(\theta)$ estimate $H_{\alpha\alpha }^{\mathcal{\ast}}(\theta)$ and $H_{\gamma\gamma}^{\mathcal{\ast}}(\theta)$, respectively. Second, denote \[ s_{i,t}^{\gamma}(\theta,\phi):=\partial_{\gamma_{t}}l_{i,t}(\theta ,\phi),\qquad s_{i,t}^{\alpha}(\theta,\phi):=\partial_{\alpha_{i}} l_{i,t}(\theta,\phi), \] which are the derivatives of the original likelihood $l_{i,t}(\theta,\phi)$ w.r.t. $\alpha_{i}$ and $\gamma_{t}$, respectively, evaluated at $\phi (\theta)$. $\widehat{S}_{\alpha\alpha}(\theta)$ and $\widehat{S}_{\gamma \gamma}(\theta)$ are defined as:
where $\widehat{\phi}(\theta):=D^{\prime}\widehat{\psi}\left( \theta\right) $ is the \textquotedblleft unparameterized" counterpart of $\widehat{\psi }\left( \theta\right) $; and
The indicator $ \mathds{1} \{\left\vert t-s\right\vert \leq\tau\}$ is a truncation mechanism, with truncation parameter $\tau$, which is common in the relevant literature. In our simulation, we use $\tau=1$ and $2$ and find the difference relatively insignificant.
Intuitively, by adding back $\widehat{B}_{\alpha}(\theta)$ and $\widehat{B} _{\gamma}(\theta)$, $L(\theta)$ serves as an approximation to the infeasible IPP-free likelihood $l(\theta)$ with rate $o_{\mathbb{P}}(1/\sqrt{NT})$. Denote $\phi(\theta):=D^{\prime}\psi(\theta)$, which is the \textquotedblleft unparameterized" counterpart of $\psi(\theta)$. We impose the following assumption and state this result formally in Theorem (ref) below.
Using the corrected likelihood $L(\theta)$, we may obtain the corrected estimator of $\theta_{0}$ as \[ \widehat{\theta}_{L}:=\arg\max_{\theta}L(\theta). \] Intuitively, since $L(\theta)$ approximates $l(\theta)$ with an error negligible relative to $1/\sqrt{NT}$, the estimator $\widehat{\theta}_{L}$ inherits this rate and is consistent and asymptotically normal with mean zero under $N/T\rightarrow\kappa$ as $N,T\rightarrow\infty$. In particular,
where $\mathcal{N}(\mu,\Sigma)$ stands for the normal distribution with mean $\mu$ and covariance matrix $\Sigma$, $\nabla_{\theta\theta^{\prime}}$ denotes the second total derivative w.r.t. to $\theta$. Note that the result in ((ref)) makes use of the information matrix equality to simplify the presentation. We refer the reader to, e.g., w1982 when such an equality does not hold.
For hypothesis testing procedures, we consider a generic null hypothesis $H_{0}:R(\theta_{0})=0$, where $R(\theta)$ is a known $r\times1$ ($r\leq K$) vector-valued non-random function with Jacobian $J(\theta):=\partial _{\theta^{\prime}}R(\theta)$ satisfying $\operatorname*{rank}[J(\theta)]=r$. The LR ($\widehat{\xi}_{LR}$), LM ($\widehat{\xi}_{LM}$), and Wald ($\widehat{\xi}_{LM}$) test statistics can be constructed as
where $\nabla_{\theta}$ denotes the first total derivative w.r.t. to $\theta$, and $\widehat{\theta}_{R}:=\arg\max_{\theta}L(\theta)$ subject to $R(\theta)=0$. The same argument above intuitively indicates that
under $H_{0}$ as $N,T\rightarrow\infty$, where $\chi^{2}(r)$ is the $\chi^{2} $-distribution with degrees of freedom $r$.
In this section, we present a simulation study. We consider dynamic binary response models specified according to the data-generating process (Design 1):
where the regressor vector is $X_{it}=\left( Y_{it-1},Z_{it}\right) ^{\prime}$ ($Z_{it}$ is described below); the parameter of interests is $\theta_{0}=(\rho_{0},\beta_{0})^{\prime}=(0.5,0.5)^{\prime}$; $\varepsilon _{it}$ is standard normal\ (for the probit model) or standard logistic (for the logit model), independent and identically distributed (i.i.d.) over $\left( i,t\right) $ and independent from the regressor $Z_{it}$; for $k=1,2$, $\{\alpha_{k,i}^{0},\gamma_{k,t}^{0}\}\sim\mathcal{N}(0,0.04)$ and are demeaned\footnote{During the estimation, we do not normalize the individual-effect parameters, because the model only contains fixed effects in the intercept.} after being generated; the regressor $Z_{it}\sim \mathcal{N}(\left( \alpha_{1,i}^{0}+\alpha_{2,i}^{0}+\gamma_{1,t}^{0} +\gamma_{2,t}^{0}\right) /2,1)$ i.i.d. over $(i,t)$; and the initial value $Y_{i0}= \mathds{1} \{(\beta_{0}+\alpha_{1,i}^{0})Z_{i0}+\alpha_{2,i}^{0}+\varepsilon_{i0}>0\}$ with $Z_{i0}\sim\mathcal{N}(\left( \alpha_{1,i}^{0}+\alpha_{2,i}^{0}\right) /2,1)$ i.i.d. over $(i,t)$. We also simulate a static version of ((ref)), where everything is the same, except that we set $\rho_{0}=0$ (hence $\theta_{0}=\beta_{0}$) and remove $Y_{it-1}$ from the regressors of estimated model. We consider $\left( N,T\right) \in\left\{ (30,30),(60,60),(90,90)\right\} $ and run $1000$ replications, comparing respectively the MLEs $\widehat{\rho}\,$and$\ \widehat{\beta}$ (of $\rho_{0}\ $and $\beta_{0}$ respectively); the corrected estimators $\widehat{\rho}_{L}^{\left( \tau\right) }\,$and$\ \widehat{\beta} _{L}^{\left( \tau\right) }$ (for dynamic models, setting $\tau\in\left\{ 1,2\right\} $), or $\widehat{\rho}_{L}\,$and$\ \widehat{\beta}_{L}$ (for static models, setting $\tau=0$); the split-panel jackknife estimators $\widehat{\rho}_{J}$ and $\widehat{\beta}_{J}$ of cfw2018; and the bootstrap-corrected estimators $\widehat{\rho}_{B}$ and $\widehat{\beta}_{B}$ of hj2024, with $499$ repetitions. We also compare the empirical sizes and powers, at the $5\%$ level, of the LR tests based on the uncorrected likelihood $\widehat{l}\left( \theta\right) $ (reporting size only), on the corrected likelihood $L\left( \theta\right) $ with $\tau=1$ (for dynamic models) or $\tau=0$ (for static model), and on the bootstrap. For the powers, the null hypotheses are $H_{0}:\theta_{0}=(0.5,0.5)^{\prime}+\delta$ for $\delta\in\left\{ \pm0.2,\pm0.1\right\} $.
Summarizing important insights in Tables (ref) and (ref), our findings are as follows.
In Section (ref) of the supplementary appendix, we report the empirical sizes and powers from the LM and the Wald test, which are similar to the LR here. We also report results from the probit model in Section (ref). In Table (ref), we present simulation results for $\left( N,T\right) =\left( 90,10\right) $ showing the performance of our approach for very small $T$. Second, we present additional simulation results in Section (ref) of the supplementary appendix for models with heterogeneous autoregressive coefficient. Next, in Section (ref) of the supplementary appendix, we present results from an alternative design (Design 2), where the biases of the MLEs are larger. Our bias correction procedure is effective under both Designs 1 and 2. In Section (ref) of the supplementary appendix, we simulate a dynamic panel probit model with additive two-way fixed effects in the intercept, under Design 1 of fw2016. Finally, in Section (ref) of the supplementary appendix, we present some additional discussion and simulation results regarding the average partial effects.
In this section we apply our likelihood-based bias correction method to examining the determinants of the labor force participation (LFP) decision of single mothers. In particular, we look at the impact of the number of children on the decision of the mother to engage in paid employment.
An extensive literature in labor economics has studied the labor supply decisions of married women killingsworth_chapter_1986,ae1998,blau_changes_2007. These studies have uncovered the impacts of a variety of economic variables, including female market wage and husband income mincer_labor_1962, education heath_causes_2016, childcare costs connelly_effect_1992, the cost of home technology greenwood_technology_2016, and culture norms fernandez_cultural_2013. Among these variables, the number of children consistently emerges as one of the most important determinants of female labor supply nakamura_econometrics_1992. This is unsurprising since women continue to bear a disproportionate share of childrearing responsibilities aguero_motherhood_2008.
In contrast to the substantial body of research on the labor supply behavior of married women, the labor supply decisions of single mothers have received relatively limited attention, with only a few exceptions kimmel_child_1998,blundell_female_2016. Single mothers, however, face distinctive challenges when it comes to balancing work and childrearing responsibilities due to the absence of a second earner in the household. Moreover, their employment decisions may have a more pronounced impact on their children's well-being compared to the decisions of married women and should thus be of great importance to economists and policymakers.
To study the labor supply decisions of single mothers, we compile a dataset from waves 20--30 of the Panel Study of Income Dynamics (PSID), which span the period of 1987 to 1997. Our sample includes only single mothers, defined as unmarried female household heads with at least one child. The dependent variable is labor force participation status, defined as whether the individual worked any hours during the interview year. Our main explanatory variable is the number of children under 18 in the household. Following dj2015, we focus on the informative sample of individuals aged 18--60 whose current and lagged LFP status each changed at least once during the period. This restriction ensures within-individual variation in both variables, which is necessary for identification of the dynamic model. After applying this criterion, the estimation sample comprises $N=86$ single mothers observed for $T=10$ consecutive years (1987--1996), with lagged participation status measured from 1986 to 1995. Further details on data construction are provided in Section (ref) of the supplementary appendix. On average, 28 percent of these individuals switched into or out of the labor force in a given year (Figure (ref)--(ref)). Table (ref) reports the summary statistics.
Let $Y_{it}\in\{0,1\}$ denote the labor-force-participation status of single mother $i$ in year $t$, and let $X_{it}$ denote her number of children. To assess the impact of fertility on labor supply, we estimate the following dynamic probit models:
where $\varepsilon_{it}$ follows the standard normal distribution, and $\sum_{i}c_{i}=\sum_{i}\zeta_{i}=\sum_{i}\alpha_{i}=\sum_{t}\eta_{t}=\sum _{t}\gamma_{t}=0$. Model (ref) is the dynamic homogeneous-slope model, which includes two-way fixed effects in the intercept. Model (ref) is the dynamic heterogeneous-slope model, which extends the specification by allowing the slope coefficients---that is, the coefficients on both the lagged dependent variable and the number of children---to vary across individuals and over time: $\rho_{it}=\rho+\zeta_{i}+\eta_{t}$ and $\beta_{it}=\beta+\alpha_{i}+\gamma_{t}$. The parameter $\beta_{it}$ measures the heterogeneous impact of the number of children on a mother's latent propensity to work, while $\rho_{it}$ captures heterogeneous state dependence in labor force participation. Model (ref) nests Model (ref) as a special case when $\rho_{it}$ and $\beta_{it}$ are constant across $i$ and $t$.
As discussed in Section (ref), when the true effect of $X_{it}$ on the outcome, $\beta_{it}$, varies across the population, it is important to account for its potential correlation with $X_{it}$. In practice, such correlation often arises from self-selection. For instance, single mothers who choose to have more children might be those with greater financial resources or better access to childcare, such that their labor supply is less affected by additional children. Alternatively, correlation between $\beta_{it}$ and $X_{it}$ may stem from common time trends: over time, fertility rates declined while concurrent factors such as rising childcare and schooling costs, advances in home technology, and evolving welfare policies (e.g., child tax credits) altered the effect of each additional child on their mother's LFP. In both scenarios, Model (ref) is the appropriate specification. It controls for unobserved heterogeneity in $\beta_{it}$ by incorporating individual and time fixed effects into the slope and allowing them to correlate with $X_{it}$. Controlling for heterogeneous state dependence further improves robustness by allowing the persistence of labor-force participation to differ across individuals and periods, capturing additional sources of dynamic heterogeneity.
Given our small sample size, direct estimation of both models could suffer from severe incidental parameter problems, which we address using our bias correction procedure. Table (ref) presents the estimation results. Panel A reports results for Model (ref) under the homogeneous slope assumption. For each parameter, we provide the MLE and bias-corrected estimates, along with their standard errors. Examining the bias-corrected estimates, we observe that an increase in the number of children appears to reduce a single mother's likelihood of labor force participation, while past LFP status strongly predicts current LFP. However, this relationship shifts dramatically once we account for slope heterogeneity. Panel B presents the results for Model (ref), allowing for heterogeneous slopes. Here, the bias-corrected estimates reveal a positive average impact of additional children on a mother's labor supply: each additional child increases the latent propensity to participate in the labor force by 0.436, which corresponds to an average increase of about 7 percent in the probability of employment. This indicates that, on average, single mothers with more children are more likely to work, likely driven by the increased financial demands of supporting multiple children independently. Additionally, unlike the homogeneous-slope model, the serial correlation in labor-force participation remains statistically significant but diminishes in magnitude, suggesting that what appears to be strong state dependence in the homogeneous model partly reflects unobserved heterogeneity in individuals' persistence of employment rather than uniform structural dynamics. Together, these findings underscore the importance of controlling for unobserved slope heterogeneity, as failing to do so can lead to substantially different conclusions.
In addition to highlighting the difference between homogeneous and heterogeneous slope models, Table (ref) demonstrates the importance of bias-correction in estimating dynamic nonlinear panel data models. Examining the results in Panel B, we observe that the MLE overestimates the impact of the number of children while underestimating the degree of serial correlation in LFP. Bias correction results in a 38 percent lower estimate of the former and a 32 percent higher estimate of the latter.
Finally, in Figure (ref), we plot the distributions of the estimated individual effects $\alpha_{i}$ and time effects $\gamma_{t}$ from Model (ref), with corresponding summary statistics reported in Table (ref) (supplementary appendix).\footnote{We acknowledge that the estimated densities of $\alpha_{i}$ and $\gamma_{t}$ may themselves be affected by the incidental-parameter problem; see oy2020.} These estimates reveal substantial heterogeneity in the effect of children on single mothers' labor force participation both across individuals and over time. Figure (ref) illustrates the temporal evolution of the slope coefficients $\beta_{it}$, constructed as the sum of the bias-corrected average effect $\beta$ and the estimated individual and time effects $\alpha_{i}$ and $\gamma_{t}$. The median of $\beta_{it}$ rises steadily from 1987 to 1994---implying an increasingly positive labor-supply response to additional children---before declining modestly thereafter. The interquartile and decile bands remain wide throughout, highlighting persistent dispersion across individuals. Together, these results reinforce the importance of accounting for both individual and temporal variations in slope coefficients, as the labor supply effects of children can vary significantly based on unobserved factors unique to each mother and time period.
In conclusion, our analysis demonstrates the efficacy of our likelihood-based bias correction method for estimating dynamic nonlinear models of labor supply with multi-dimensional heterogeneities. We show the importance of controlling for unobserved slope heterogeneity when it may correlate with regressors and the value of panel data models that accommodate two-way fixed effects in both the intercept and the slope. Our findings reveal that the number of children has, on average, a positive influence on a single mother's labor force participation, but such effect is highly heterogeneous across individuals and time periods. This heterogeneity underscores the complex nature of labor supply decisions and the unique challenges each single mother faces in balancing work and childrearing responsibilities.
{10pt}
In this paper, we propose a likelihood-based analytical bias correction procedure for a general class of two-way dynamic nonlinear heterogeneous parameter models subject to the incidental parameter problem. We give the analytical form of a corrected likelihood and show that it delivers point estimators that are asymptotically unbiased and test statistics that are asymptotically $\chi^{2}$-distributed. Simulation studies and empirical analyses support our claims.
Several issues deserve further studies. First, we have not rigorously investigated the asymptotic properties of the average partial effects. Admittedly, average partial effect is very important, especially for nonlinear models. It is therefore worthwhile to formally establish its asymptotic theory. Second, there has been various fixed-$T$ consistent solutions for the logit model with fixed effects. However, to the best of our knowledge, we have not seen any such developments for the two-way heterogeneous parameter logit model. It could be useful to push such a development, because the logit model is widely used in many studies. Finally, certain micro panels tend to have a small $T$, in which case the higher-order bias terms may be non-negligible. For this reason, it may be potentially important to develop high-order bias correction approaches for two-way dynamic nonlinear heterogeneous parameter models.
We would like to thank the editor Iv \'a n Fern \'a ndez-Val and two anonymous referees for their valuable comments and suggestions. We would like to express our gratitude to comments and suggestions from Geert Dhaene, Yanqin Fan, Jinyong Hahn, Chen Hsiao, Whitney Newey, and Jun Yu. Our research assistants Yuchi Liu and Xinrui Zhang helped us with various tasks. Yutao Sun is the sole corresponding author of this paper. All correspondence shall be sent to Yutao Sun. Leng's research is partially supported by the National Natural Science Foundation of China (72573136), the Fujian Natural Science Foundation in China (2024J08013), and the Basic Scientific Center Project of National Natural Science Foundation of China (71988101). Mao acknowledges support from the NSFC Basic Science Center Project for Econometric Modeling and Economic Policy Studies (Grant No. 71988101). Sun acknowledges the support from the National Natural Science Foundation of China under Grant Number 72203032.