The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
70,161 characters
Debiased Inference for Dynamic Nonlinear Panels with Multi-dimensional Heterogeneities
\title
{Debiased Inference for Dynamic Nonlinear Panels with Multi-dimensional Heterogeneities}
\date{\today}
\author{
Xuan Leng
\thanks
{Department of Statistics and Data Science at School of Economics, Xiamen University,
Xiamen (361000), China. Email: [email removed].}
\and Jiaming Mao
\thanks{Wang Yanan Institute for Studies in Economics, Xiamen University,
Xiamen (361000), China. Email: [email removed]
.}
\and Yutao Sun
\thanks
{Institute for Advanced Economic Research, Dongbei University of Finance and Economics,
Dalian (116025), China. Email: [email removed].}
\thanks{Corresponding Author.}
}
\maketitle
\begin{abstract}
We introduce a generic class of dynamic nonlinear heterogeneous parameter
models that incorporate individual and time fixed effects in both the
intercept and slope. These models are subject to the incidental parameter
problem, in that the limiting distribution of the point estimator is not
centered at zero, and that test statistics do not follow their standard
asymptotic distributions as in the absence of the fixed effects. To address
the problem, we develop an analytical bias correction procedure to construct a
bias-corrected likelihood. The resulting estimator follows an asymptotic
normal distribution with mean zero. Moreover, likelihood-based test
statistics---including likelihood-ratio, Lagrange-multiplier, and Wald
tests---follow the limiting chi-squared distribution under the null
hypothesis. Simulations demonstrate the effectiveness of the proposed
correction method, and an empirical application on the labor force
participation of single mothers underscores its practical importance.
\begin{description}
\item[\textbf{Keywords:}] incidental parameter problem; bias correction; fixed
effects; panel data models
\item[\textbf{JEL Classification}:] C23
\end{description}
\end{abstract}
\section{Introduction}
\label{Section.Introduction}
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 \emph{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 \citep{bc2007}. Cross-sectional response
heterogeneity reflects individual-specific factors such as preferences,
productivity, and access to resources \citep{heckman2001}, whereas
time-varying response heterogeneity captures changes in behavior driven by
evolving circumstances---including policy reforms, business cycles, and
technological progress \citep{ow2021}. When such \emph{response
heterogeneities} are correlated with covariates, failing to account for them
can lead to biased estimates and invalid inference \citep{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 \citep{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.
\label{moremotivation}In this paper, we introduce a general class of
\emph{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 \citep{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 \citep{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 \citep{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{Result.EstLogitPreview} 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.
\begin{center}
\begin{figure}[H]
\caption
{Comparison of Uncorrected and Corrected Maximum Likelihood Estimators}
\label{Result.EstLogitPreview}
\includegraphics[width=0.35\linewidth]{BoxPlot-Param2}
\qquad
\includegraphics[width=0.35\linewidth]{BoxPlot-Param3}
\begin{flushleft}
\begin{justify}
\begin{footnotesize}
\noindent\textit{Notes}: Model. $Y_{it}=
\mathds{1}
\{\rho Y_{it-1}+(\beta+\alpha_{1,i}+\gamma_{1,t})Z_{it}+\alpha_{2,i}
+\gamma_{2,t}+\varepsilon_{it}>0\}$ where $
\mathds{1}
\{\cdot\}$ denotes the indicator function, $\varepsilon_{it}$ is standard
logistic independent of the exogenous regressor $Z_{it}$, and the true values
$\rho_{0}=\beta_{0}=0.5$. $\widehat{\beta}$ and $\widehat{\rho}$ are the
uncorrected estimators while $\widehat{\beta}_{L}$ and $\widehat{\rho}_{L}$
are the corrected estimators. $1000$ replications. The detailed
data-generating process is given in Section \ref{Section.Simulation}.
\end{footnotesize}
\end{justify}
\end{flushleft}
\end{figure}
\end{center}
\emph{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 \citet{c1980}, \citet{ab1991}, and \citet{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. \citet{hn2004} and
\citet{hk2011} propose \textquotedblleft parameter-based" corrections that
remove the bias from the likelihood estimator; \citet{w2002} and
\citet{llw2003} propose \textquotedblleft score-based" methods that modify the
profiled score; while \textquotedblleft likelihood-based" approaches such as
\citet{bh2009} and \citet{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 \citep{dj2015}, the
bootstrap \citep{ks2016,bkss2020,hj2024}, and integrated-likelihood methods \citep{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. \citet{mw2015a} develop a
parameter-based analytical correction for dynamic linear models with
interactive effects. For nonlinear models, \citet{fw2016} propose
parameter-based techniques for additive fixed effects, while \citet{cfw2014}
study interactive effects in static models. Alternatively, \citet{ko2018}
provide a likelihood-based approach for static models with an arbitrary but
known fixed-effect structure. Other model-specific contributions include
\citet{b2009} and \citet{c2017}. For a comprehensive overview of these
developments, see \citet{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
\citep{robertson1992, ps1995} and time-varying coefficients
\citep{robinson1989,sw1996}, most often in linear or static contexts
\citep{h2014, pesaran2015}. Recent research in dynamic nonlinear panels has
advanced these ideas further. \citet{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 \citet{fglv2025} for a panel distribution
regression with individual-heterogeneous coefficient. \citet{cfhn2013} derive
partial identification with uniform inference for nonseparable panels under
time-homogeneity. \citet{bc2014} provide mixture-based point identification
for dynamic binary models with maximal cross-sectional heterogeneity, while
\citet{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. \citet{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. \citet{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.
\emph{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 \cite{ah2016} and
\cite{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{Section.IPP}
explains our settings, the identification restriction, and the incidental
parameter problem of DN-HP models. Section \ref{Section.Asymptotics} describes
our bias correction procedure and provides relevant statistical properties.
Section \ref{Section.Simulation} presents simulation studies to demonstrate
the performance of the method. Section \ref{Section.Empirical} applies our
method to study single mothers' labor force participation. Finally in Section
\ref{Section.Conclusion}, 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.
\emph{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$.
\section{Bias in Heterogeneous Parameter Models}
\label{Section.IPP}
\subsection{Model and Estimation}
\label{ME}
In this section, we explain our models and the estimation procedure, leaving
two motivating examples in Section
\ref{Section.TechDetails.ModelAndEstimation} 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
\begin{equation}
Y_{it}|X_{it},\ldots,X_{i1},\phi^{0}\sim f(y|X_{it},\theta_{0}+\alpha_{i}
^{0}+\gamma_{t}^{0}), \label{Eq.IPP.DGP}
\end{equation}
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{Eq.IPP.DGP}) 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{Eq.IPP.DGP}), 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
\begin{equation}
\ell(\theta,\phi):=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}l_{i,t}(\theta
,\phi),\qquad l_{i,t}(\theta,\phi):=\log f(Y_{it}|X_{it},\theta+\alpha
_{i}+\gamma_{t}). \label{Eq.IPP.OriginalLike}
\end{equation}
Due to the additive structure $\theta+\alpha_{i}+\gamma_{t}$, $\theta$ is not
identified because the likelihood function is invariant to transformations
\begin{equation}
\alpha_{i}\mapsto\alpha_{i}+c_{1},\qquad\gamma_{t}\mapsto\gamma_{t}
+c_{2},\qquad\theta\mapsto\theta-c_{1}-c_{2} \label{Eq.IPP.RotInd}
\end{equation}
for all pairs of constant $(c_{1},c_{2})$ and all $(i,t)$. To resolve this, we
follow \cite{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{Eq.IPP.RotInd}) without explicitly using the Lagrangian.
\subsection{Incidental Parameter Problem}
\label{Section.IPP.IPPBias}
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{Section.TechDetails.BiasExpansion} 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 \cite{ps2006}, we continue our discussion by
defining, for a given $\theta$ and each pair of $(N,T)$,
\begin{equation}
\widehat{\psi}(\theta):=\arg\max_{\psi}l(\theta,\psi),\qquad\psi(\theta
):=\arg\max_{\psi}\mathbb{E}[l(\theta,\psi)]. \label{Eq.IPP.PsiEstimator}
\end{equation}
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., \cite{hn2004} and
\cite{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)$:
\begin{equation}
l(\theta,\widehat{\psi}(\theta))=l(\theta,\psi(\theta))+s^{\prime}
(\theta)[\widehat{\psi}(\theta)-\psi(\theta)]+\frac{1}{2}[\widehat{\psi
}(\theta)-\psi(\theta)]^{\prime}H(\theta)[\widehat{\psi}(\theta)-\psi
(\theta)]+r^{l}(\theta), \label{Eq.IPP.LikeExpansion}
\end{equation}
where $r^{l}(\theta)$ is a remainder, $\widehat{\psi}(\theta)-\psi(\theta)$
represents the estimation errors of the fixed effects,
\begin{align*}
s(\theta) & :=s(\theta,\psi(\theta)),\qquad s(\theta,\psi):=\partial_{\psi
}l(\theta,\psi),\\
H\left( \theta\right) & :=H(\theta,\psi(\theta)),\qquad H(\theta
,\psi):=\partial_{\psi\psi^{\prime}}l(\theta,\psi)
\end{align*}
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
\begin{align}
0 & =s\left( \theta\right) +H\left( \theta\right) [\widehat{\psi}
(\theta)-\psi(\theta)]+r^{s}(\theta),\nonumber\\
\widehat{\psi}(\theta)-\psi(\theta) & =-[H\left( \theta\right)
]^{-1}s\left( \theta\right) -[H\left( \theta\right) ]^{-1}r^{s}(\theta),
\label{Eq.IPP.PsiExpansion}
\end{align}
where $r^{s}(\theta)$ is a remainder. Combining (\ref{Eq.IPP.PsiExpansion})
and (\ref{Eq.IPP.LikeExpansion}), and taking expectation, we obtain
\begin{align}
\sqrt{NT}\mathbb{E}[l(\theta,\widehat{\psi}(\theta))-l(\theta,\psi(\theta))]
& =-\sqrt{NT}B(\theta)+\sqrt{NT}\mathbb{E}r(\theta),\nonumber\\
B(\theta) & :=\frac{1}{2}\mathbb{E\{}s^{\prime}(\theta)[H(\theta
)]^{-1}s(\theta)\}, \label{Eq.IPP.BiasExpansion}
\end{align}
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.
\section{Bias Correction and Asymptotic Theory}
\label{Section.Asymptotics}
The bias expansion of $\widehat{l}(\theta):=l(\theta,\widehat{\psi}(\theta))$
in Equation (\ref{Eq.IPP.BiasExpansion}) 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{Section.TechDetails.BiasTerms} of the supplementary appendix i) more
technical details; ii) formal statements and discussions about our secondary
results (Equations \ref{Eq.Asymptotics.LikeExpansion},
\ref{Eq.Asymptotic.NormalTheta}, \ref{Eq.Asymptotic.Chi2Tests}, etc.); iii)
useful remarks; etc. In addition, we explain the relation of our method to
\cite{bh2009} and \cite{ah2016} in Section \ref{Section.AH2016} 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
\begin{equation}
\mathbb{E}l(\theta)=\mathbb{E}\widehat{l}(\theta)+B_{\alpha}(\theta
)+B_{\gamma}(\theta)+o(1/\sqrt{NT}), \label{Eq.Asymptotics.LikeExpansion}
\end{equation}
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{\label{Footnote.AsymBiasInKappa}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., \cite{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
\begin{equation}
\lbrack\mathbb{E}H(\theta)]^{-1}=:\left(
\begin{array}
[c]{cc}
H_{\alpha\alpha}^{\mathcal{\ast}}(\theta) & H_{\alpha\gamma}^{\mathcal{\ast}
}(\theta)\\
H_{\gamma\alpha}^{\mathcal{\ast}}(\theta) & H_{\gamma\gamma}^{\mathcal{\ast}
}(\theta)
\end{array}
\right) ,\qquad s(\theta)-\mathbb{E}s\left( \theta\right) =:\left(
\begin{array}
[c]{c}
\widetilde{s}_{\alpha}(\theta)\\
\widetilde{s}_{\gamma}(\theta)
\end{array}
\right) . \label{Eq.Asymptotics.Partitions}
\end{equation}
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
\begin{align}
L(\theta) & :=\widehat{l}(\theta)+\widehat{B}_{\alpha}(\theta)+\widehat{B}
_{\gamma}(\theta),\label{Eq.Asymptotic.CorrectedLike}\\
\widehat{B}_{\alpha}(\theta) & :=\frac{1}{2}\operatorname*{trace}
[\widehat{S}_{\alpha\alpha}(\theta)\widehat{H}_{\alpha\alpha}^{\mathcal{\ast}
}(\theta)],\qquad\widehat{B}_{\gamma}(\theta):=\frac{1}{2}
\operatorname*{trace}\mathbb{[}\widehat{S}_{\gamma\gamma}(\theta
)\widehat{H}_{\gamma\gamma}^{\ast}(\theta)].\nonumber
\end{align}
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
\begin{equation}
\lbrack\widehat{H}(\theta)]^{-1}=:\left(
\begin{array}
[c]{cc}
\widehat{H}_{\alpha\alpha}^{\mathcal{\ast}}(\theta) & \widehat{H}
_{\alpha\gamma}^{\mathcal{\ast}}(\theta)\\
\widehat{H}_{\gamma\alpha}^{\mathcal{\ast}}(\theta) & \widehat{H}
_{\gamma\gamma}^{\mathcal{\ast}}(\theta)
\end{array}
\right) \label{Eq.Asymptotic.Hessian}
\end{equation}
in the same manner as in (\ref{Eq.Asymptotics.Partitions}). 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:
\begin{align}
\widehat{S}_{\alpha\alpha}(\theta) & :=D_{1}\widehat{\mathcal{S}}
_{\alpha\alpha}(\theta)D_{1}^{\prime},\label{Eq.Asymptotics.SaaHat}\\
\widehat{\mathcal{S}}_{\alpha\alpha}(\theta) & :=\frac{1}{N^{2}T^{2}}
\sum_{t=1}^{T}\sum_{s=1}^{T}
\mathds{1}
\{\left\vert t-s\right\vert \leq\tau\}\operatorname*{diag}\{\widehat{s}
_{1,t}^{\alpha}(\theta)[\widehat{s}_{1,s}^{\alpha}(\theta)]^{\prime}
,\ldots,\widehat{s}_{N,t}^{\alpha}(\theta)[\widehat{s}_{N,s}^{\alpha}
(\theta)]^{\prime}\},\nonumber\\
\widehat{s}_{i,t}^{\alpha}(\theta) & :=s_{i,t}^{\alpha}(\theta
,\widehat{\phi}(\theta))-\frac{1}{T}\sum_{t=1}^{T}s_{i,t}^{\alpha}
(\theta,\widehat{\phi}(\theta)),\nonumber
\end{align}
where $\widehat{\phi}(\theta):=D^{\prime}\widehat{\psi}\left( \theta\right)
$ is the \textquotedblleft unparameterized" counterpart of $\widehat{\psi
}\left( \theta\right) $; and
\begin{align}
\widehat{S}_{\gamma\gamma}(\theta) & :=D_{2}\widehat{\mathcal{S}}
_{\gamma\gamma}(\theta)D_{2}^{\prime},\label{Eq.Asymptotics.SggHat}\\
\widehat{\mathcal{S}}_{\gamma\gamma}(\theta) & :=\frac{1}{N^{2}T^{2}}
\sum_{i=1}^{N}\widehat{s}_{i}^{\gamma}(\theta)[\widehat{s}_{i}^{\gamma}
(\theta)]^{\prime},\qquad\widehat{s}_{i}^{\gamma}(\theta):=\{[\widehat{s}
_{i,1}^{\gamma}(\theta)]^{\prime},\ldots,[\widehat{s}_{i,T}^{\gamma}
(\theta)]^{\prime}\}^{\prime},\nonumber\\
\widehat{s}_{i,t}^{\gamma}(\theta) & :=s_{i,t}^{\gamma}(\theta
,\widehat{\phi}(\theta))-\frac{1}{N}\sum_{i=1}^{N}s_{i,t}^{\gamma}
(\theta,\widehat{\phi}(\theta)).\nonumber
\end{align}
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.
\begin{remark}
[efficient computation of bias terms]\label{Remark.Computation1}The terms
$\widehat{S}_{\alpha\alpha}(\theta)$ and $\widehat{S}_{\gamma\gamma}(\theta)$
are constructed using the derivatives of the original (i.e., \textquotedblleft
unparameterized") likelihood. This is valid because the reparameterization
produces $\partial_{\psi}l(\theta,\psi)=D\partial_{\phi}\ell(\theta,\phi)$ and
$\partial_{\psi\psi^{\prime}}l(\theta,\psi)=D\partial_{\phi\phi^{\prime}}
\ell(\theta,\phi)D^{\prime}$ holding true for every $\phi=D^{\prime}\psi$ and
every $\theta$. Generally, it may be easier to construct $\widehat{S}
_{\alpha\alpha}(\theta)$ and $\widehat{S}_{\gamma\gamma}(\theta)$ from the
original likelihood, because analytical expressions for $\partial_{\phi}
\ell(\theta,\phi)$ and $\partial_{\phi\phi^{\prime}}\ell(\theta,\phi)$ are
well-known for many frequently used models (e.g., the probit, the logit, and
the Poisson).
In addition, if the model is equipped with linear indices, these derivatives
may be calculated efficiently using the chain rule. We have $\partial_{\phi
}l_{i,t}(\theta,\phi)=\partial_{\pi}l_{i,t}(\pi_{i,t}\left( \theta
,\phi\right) )\partial_{\phi}\pi_{i,t}\left( \theta,\phi\right) $ by
viewing $l_{i,t}(\theta,\phi):=l_{i,t}(\pi_{i,t}\left( \theta,\phi\right)
)$, where $\pi_{i,t}\left( \theta,\phi\right) :=X_{it}^{\prime}
(\theta+\alpha_{i}+\gamma_{t})$ is the linear index. Here $\partial_{\phi}
\pi_{i,t}\left( \theta,\phi\right) $ is high-dimensional but can be
calculated easily, because $\pi_{i,t}\left( \theta,\phi\right) $ is only
linear. The more complex component $\partial_{\pi}l_{i,t}(\cdot)$ is only
scalars and, therefore, can be calculated afforably, even with numerical
differentiation. The complexity of this only grows with the number of linear
indices, but not the dimension of $\phi$.
\end{remark}
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{Theorem.CorrectedLike} below.
\begin{asu}
\
\begin{itemize}[leftmargin=*,itemindent=-2em]
\label{Assumption.Expansion}
\item[]\refstepcounter{subassumption}(\textit{\roman{subassumption}})\
\label{Assumption.Expansion.AsySeq}Suppose $N/T\rightarrow\kappa$ for some
$0<\kappa<\infty$ as $N,T\rightarrow\infty$.
\item[]\refstepcounter{subassumption}(\textit{\roman{subassumption}})\
\label{Assumption.Expansion.Smooth}Let $\Theta$ be a compact subset of
$\mathbb{R}^{K}$ with $\theta_{0}$ in its interior and$\ \Phi$ be a compact
subset of $\mathbb{R}^{(N+T)K}$. For each $\theta\in\Theta$, $\Phi$ contains
both $\widehat{\phi}(\theta)$ and $\phi(\theta)$ in its interior.
$l_{i,t}(\theta,\phi)$ is three-time continuously differentiable w.r.t. to
$\phi\in\Phi$. There exists a function $g(w_{it})$, $w_{it}=(Y_{it}
,X_{it}^{\prime})^{\prime}$, independent of $\theta$ and $\phi$, such that
\[
\sup_{\theta\in\Theta}\sup_{\phi\in\Phi}\left\vert \partial_{\phi_{1}
\cdots\phi_{S}}l_{i,t}(\theta,\phi)\right\vert <g(w_{it}),
\]
where $\phi_{s}$ represents any element of $\phi$ and $S\in\{0,1,2,3\}$ ($S=0$
is understood as not taking derivatives), and
\[
\sup_{N,T}\max_{1\leq i\leq N,1\leq t\leq T}\mathbb{E}_{\phi}\{[g(w_{it}
)]^{\eta}\}<\infty
\]
for some $\eta>2,$ where $\mathbb{E}_{\phi}$ denotes the conditional
expectation w.r.t. the joint distribution of $w_{it}$, given the heterogeneous
effects $\phi^{0}$.
\item[]\refstepcounter{subassumption}(\textit{\roman{subassumption}})\
\label{Assumption.Expansion.Normalize}For each $k=1,\ldots,K$, $\sum_{i=1}
^{N}\alpha_{k,i}^{0}=\sum_{t=1}^{T}\gamma_{k,t}^{0}=0$, where $\alpha
_{k,i}^{0}$ and $\gamma_{k,t}^{0}$ are the $k$-th component of $\alpha_{i}
^{0}$ and $\gamma_{t}^{0}$, respectively.
\item[]\refstepcounter{subassumption}(\textit{\roman{subassumption}})\
\label{Assumption.Expansion.Ident} For each $\theta\in\Theta$, $\mathbb{P}
[l(\theta,\psi)\neq l(\theta,\psi(\theta))]>0$ for every $\psi$ such that
$D^{\prime}\psi\in\Phi$ and $\psi\neq\psi(\theta)$.
\item[]\refstepcounter{subassumption}(\textit{\roman{subassumption}})\
\label{Assumption.Expansion.Data}Conditional on $\phi^{0}$, $\{w_{it}\}$ are
independent across $i$ and conditionally strong mixing with mixing coefficient
\[
a_{i}(m):=\sup_{t\geq1}\sup_{A\in\mathcal{A}_{it},B\in\mathcal{B}_{it+m}
}\left\vert \mathbb{P}(A\cap B|\phi^{0})-\mathbb{P}(A|\phi^{0})\mathbb{P}
(B|\phi^{0})\right\vert
\]
such that, for $\eta>2$ as above,
\[
\sup_{N}\max_{1\leq i\leq N}\sum_{m=0}^{\infty}[a_{i}(m)]^{1-2/\eta}<\infty,
\]
where $\mathcal{A}_{it}$ and $\mathcal{B}_{it}$ are the $\sigma$-algebras
generated by $\{w_{is}:1\leq s\leq t\}$ and $\{w_{is}:t\leq s\leq T\}$, respectively.
\item[]\refstepcounter{subassumption}(\textit{\roman{subassumption}})\
\label{Assumption.Expansion.Hessian}The matrix $\sqrt{NT}\mathbb{E}
[H(\theta)]$ has eigenvalues $\lambda_{p}(\theta)$ for $p=1,\ldots,(N+T-2)K$
which satisfy
\[
\sup_{N>N_{0},T>T_{0}}\max_{1\leq p\leq(N+T-2)K}\sup_{\theta\in\Theta}
\lambda_{p}(\theta)<0
\]
a.s. for some integers $N_{0}$ and $T_{0}$.
\end{itemize}
\end{asu}
\begin{theorem}
\label{Theorem.CorrectedLike}Under Assumption \ref{Assumption.Expansion} and
for some $\tau\rightarrow\infty$ such that $\tau/T\rightarrow0$, we have
\[
L(\theta)-l(\theta)=o_{\mathbb{P}}(1/\sqrt{NT}),
\]
uniformly over $\theta\in\Theta$ as $N,T\rightarrow\infty$, where $L(\theta)$
is the corrected likelihood defined in Equation
(\ref{Eq.Asymptotic.CorrectedLike}).
\end{theorem}
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,
\begin{equation}
\sqrt{NT}[-\mathbb{E}\triangledown_{\theta\theta^{\prime}}l(\theta_{0}
)]^{1/2}(\widehat{\theta}_{L}-\theta_{0})\overset{\mathbb{D}}{\longrightarrow
}\mathcal{N}(0,\mathbb{I}_{K}), \label{Eq.Asymptotic.NormalTheta}
\end{equation}
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{Eq.Asymptotic.NormalTheta}) makes use of the information matrix equality
to simplify the presentation. We refer the reader to, e.g., \cite{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
\begin{align*}
\widehat{\xi}_{LR} & :=-2NT\{L(\widehat{\theta}_{R})-L(\widehat{\theta}
_{L})\},\\[4pt]
\widehat{\xi}_{LM} & :=-NT\{\triangledown_{\theta^{\prime}}L(\widehat{\theta
}_{R})[\triangledown_{\theta\theta^{\prime}}L(\widehat{\theta}_{R}
)]^{-1}\triangledown_{\theta}L(\widehat{\theta}_{R})\},\\[4pt]
\widehat{\xi}_{Wald} & :=-NT\{R^{\prime}(\widehat{\theta}_{L}
)[J(\widehat{\theta}_{L})[\triangledown_{\theta\theta^{\prime}}
L(\widehat{\theta}_{L})]^{-1}J^{\prime}(\widehat{\theta}_{L})]^{-1}
R(\widehat{\theta}_{L})\},
\end{align*}
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
\begin{equation}
\widehat{\xi}_{LR},\widehat{\xi}_{LM},\widehat{\xi}_{Wald}\overset{\mathbb{D}
}{\longrightarrow}\chi^{2}(r), \label{Eq.Asymptotic.Chi2Tests}
\end{equation}
under $H_{0}$ as $N,T\rightarrow\infty$, where $\chi^{2}(r)$ is the $\chi^{2}
$-distribution with degrees of freedom $r$.
\section{Simulation}
\label{Section.Simulation}
In this section, we present a simulation study. We consider dynamic binary
response models specified according to the data-generating process (Design 1):
\begin{equation}
Y_{it}=
\mathds{1}
\{\rho_{0}Y_{it-1}+(\beta_{0}+\alpha_{1,i}^{0}+\gamma_{1,t}^{0})Z_{it}
+\alpha_{2,i}^{0}+\gamma_{2,t}^{0}+\varepsilon_{it}>0\},
\label{Equation.Simulation.BinResp}
\end{equation}
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{Equation.Simulation.BinResp}), 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 \cite{cfw2018}; and the
bootstrap-corrected estimators $\widehat{\rho}_{B}$ and $\widehat{\beta}_{B}$
of \cite{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{Table.Simulation.Design1Est} and
\ref{Table.Simulation.Design1LR}, our findings are as follows.
\begin{enumerate}
\item The MLEs $\widehat{\beta}$ and $\widehat{\rho}$ may exhibit significant
bias, especially when the sample size is small. Our procedure reduces the bias
considerably without inflating the RMSE, even at $\left( N,T\right) =\left(
30,30\right) $.
\item Our corrected estimators show smaller biases than the jackknife and the
bootstrap for $\beta$ when the sample size is small. The biases of all
candidate estimators, except the MLEs, become similar as the sample size
increases. For $\left( N,T\right) =\left( 30,30\right) $, the jackknife
estimators may show larger biases than the MLEs for structural parameters
associated with two-way heterogeneities. We discuss the possible cause in
Section \ref{Section.Jackknife} of the supplementary appendix.
\item The LR test based on $\widehat{l}\left( \theta\right) $ has severe
size-distortions even with $\left( N,T\right) =\left( 90,90\right) $. On
the contrary, the LR test based on our $L\left( \theta\right) $ is able to
deliver an empirical size close to the nominal level of $5\%$ as the sample
size increases, maintaining relatively strong powers. For small sample sizes,
$L\left( \theta\right) $ reduces the type-I error risk considerably. The
empirical sizes from the bootstrap are close to the nominal level with small
sample sizes. However, the powers are relatively low, especially near the true
null hypothesis.
\end{enumerate}
\begin{table}[H]
\begin{centering}
\caption{Comparisons of Uncorrected and Corrected Logit Estimates, Design 1}
\label{Table.Simulation.Design1Est}
\setstretch{1.2}
\footnotesize
\begin{tabular*}
{\linewidth}[c]{@{\extracolsep{\fill}}llrrcrrcrr}\hline\hline
$(N,T)$ & & \multicolumn{2}{c}{$(30,30)$} & & \multicolumn{2}{c}{$(60,60)$}
& & \multicolumn{2}{c}{$(90,90)$}\\
& & Bias & RMSE & \multicolumn{1}{r}{} & Bias & RMSE & \multicolumn{1}{r}{} &
Bias & RMSE\\\hline
& & \multicolumn{8}{c}{Static Logit Model}\\
$\widehat{\beta}$ & & $0.165$ & $0.197$ & \multicolumn{1}{r}{} & $0.073$ &
$0.086$ & \multicolumn{1}{r}{} & $0.045$ & $0.053$\\
$\widehat{\beta}_{L}$ & & $0.055$ & $0.106$ & \multicolumn{1}{r}{} & $0.015$
& $0.044$ & \multicolumn{1}{r}{} & $0.006$ & $0.027$\\
$\widehat{\beta}_{J}$ & & $-0.776$ & $1.716$ & \multicolumn{1}{r}{} &
$-0.022$ & $0.049$ & \multicolumn{1}{r}{} & $-0.008$ & $0.028$\\
$\widehat{\beta}_{B}$ & & $-0.069$ & $0.098$ & \multicolumn{1}{r}{} &
$-0.008$ & $0.040$ & \multicolumn{1}{r}{} & $-0.004$ & $0.026$\\\hline
& \multicolumn{1}{c}{} & \multicolumn{8}{c}{Dynamic Logit Model}\\
$\widehat{\beta}$ & & $0.170$ & $0.210$ & \multicolumn{1}{r}{} & $0.072$ &
$0.087$ & \multicolumn{1}{r}{} & $0.046$ & $0.054$\\
$\widehat{\beta}_{L}^{(1)}$ & & $0.061$ & $0.120$ & \multicolumn{1}{r}{} &
$0.016$ & $0.046$ & \multicolumn{1}{r}{} & $0.007$ & $0.027$\\
$\widehat{\beta}_{L}^{(2)}$ & & $0.063$ & $0.122$ & \multicolumn{1}{r}{} &
$0.017$ & $0.046$ & \multicolumn{1}{r}{} & $0.007$ & $0.027$\\
$\widehat{\beta}_{J}$ & & $-1.014$ & $2.120$ & \multicolumn{1}{r}{} &
$-0.029$ & $0.111$ & \multicolumn{1}{r}{} & $-0.009$ & $0.029$\\
$\widehat{\beta}_{B}$ & & $-0.079$ & $0.110$ & \multicolumn{1}{r}{} &
$-0.010$ & $0.042$ & \multicolumn{1}{r}{} & $-0.004$ & $0.026$\\
$\widehat{\rho}$ & & $-0.109$ & $0.209$ & \multicolumn{1}{r}{} & $-0.051$ &
$0.097$ & \multicolumn{1}{r}{} & $-0.030$ & $0.059$\\
$\widehat{\rho}_{L}^{(1)}$ & & $-0.047$ & $0.165$ & \multicolumn{1}{r}{} &
$-0.018$ & $0.079$ & \multicolumn{1}{r}{} & $-0.007$ & $0.049$\\
$\widehat{\rho}_{L}^{(2)}$ & & $-0.042$ & $0.168$ & \multicolumn{1}{r}{} &
$-0.013$ & $0.079$ & \multicolumn{1}{r}{} & $-0.003$ & $0.049$\\
$\widehat{\rho}_{J}$ & & $0.010$ & $0.177$ & \multicolumn{1}{r}{} & $0.002$ &
$0.081$ & \multicolumn{1}{r}{} & $0.003$ & $0.051$\\
$\widehat{\rho}_{B}$ & & $0.009$ & $0.153$ & \multicolumn{1}{r}{} & $0.000$ &
$0.077$ & \multicolumn{1}{r}{} & $0.003$ & $0.049$\\\hline\hline
\end{tabular*}
\begin{footnotesize}
\justify
\end{footnotesize}
\end{centering}
\end{table}
\begin{table}[H]
\begin{centering}
\caption{Empirical Size and Power from LR Test for Logit Model, Design 1}
\label{Table.Simulation.Design1LR}
\setstretch{1.2}
\scriptsize
\begin{tabular*}
{\linewidth}[c]{@{\extracolsep{\fill}}cccccccccccc}\hline\hline
& \multicolumn{3}{c}{Size} & \multicolumn{4}{c}{Power at $\delta$
(Analytical)} & \multicolumn{4}{c}{Power at $\delta$ (Bootstrap)}\\
$(N,T)$ & Uncorrected & Analytical & Bootstrap & $-0.2$ & $-0.1$ & $0.1$ &
$0.2$ & $-0.2$ & $-0.1$ & $0.1$ & $0.2$\\\hline
& \multicolumn{11}{c}{Static Logit Model}\\
$(30,30)$ & $38$ & $13$ & $0$ & $88$ & $48$ & $12$ & $45$ & $53$ & $12$ & $0$
& $0$\\
$(60,60)$ & $40$ & $8$ & $2$ & $100$ & $85$ & $59$ & $99$ & $99$ & $74$ & $0$
& $32$\\
$(90,90)$ & $38$ & $5$ & $3$ & $100$ & $98$ & $95$ & $100$ & $100$ & $98$ &
$9$ & $97$\\\hline
& \multicolumn{11}{c}{Dynamic Logit Model}\\
$(30,30)$ & $37$ & $13$ & $1$ & $84$ & $42$ & $19$ & $55$ & $42$ & $8$ & $1$ &
$4$\\
$(60,60)$ & $36$ & $10$ & $3$ & $100$ & $79$ & $64$ & $99$ & $99$ & $65$ & $7$
& $68$\\
$(90,90)$ & $33$ & $6$ & $3$ & $100$ & $98$ & $95$ & $100$ & $100$ & $96$ &
$34$ & $99$\\\hline\hline
\end{tabular*}
\begin{footnotesize}
\justify
\end{footnotesize}
\end{centering}
\end{table}
In Section \ref{Section.AdditionalSim.Design1} 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{Section.AdditionalSim.Design1}. In Table
\ref{Table.Simulation.Design1UnequalNT}, 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{Section.AdditionalSim.Poisson} of the supplementary appendix
for models with heterogeneous autoregressive coefficient. Next, in Section
\ref{Section.AdditionalSim.Design2} 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{Section.AdditionalSim.FVW2016} of the supplementary
appendix, we simulate a dynamic panel probit model with additive two-way fixed
effects in the intercept, under Design 1 of \cite{fw2016}. Finally, in Section
\ref{Section.APE} of the supplementary appendix, we present some additional
discussion and simulation results regarding the average partial effects.
\section{Empirical Analysis}
\label{Section.Empirical}
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
\citep{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 \citep{mincer_labor_1962}, education
\citep{heath_causes_2016}, childcare costs \citep{connelly_effect_1992}, the
cost of home technology \citep{greenwood_technology_2016}, and culture norms
\citep{fernandez_cultural_2013}. Among these variables, the number of children
consistently emerges as one of the most important determinants of female labor
supply \citep{nakamura_econometrics_1992}. This is unsurprising since women
continue to bear a disproportionate share of childrearing responsibilities \citep{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
\citep{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
\citet{dj2015}, we focus on the \textit{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{sec:APP} 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{fig:lfpct}--\ref{fig:lfpflow}).
Table \ref{tab:DS} 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:
\begin{align}
Y_{it} & =\mathds{1}\!\left\{ \rho Y_{i,t-1}+\beta X_{it}+c_{i}+\tau
_{t}+\varepsilon_{it}>0\right\} ,\label{eq:em01}\\
Y_{it} & =\mathds{1}\!\left\{ (\rho+\zeta_{i}+\eta_{t})Y_{i,t-1}
+(\beta+\alpha_{i}+\gamma_{t})X_{it}+c_{i}+\tau_{t}+\varepsilon_{it}
>0\right\} , \label{eq:em02}
\end{align}
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 \eqref{eq:em01} is the dynamic homogeneous-slope
model, which includes two-way fixed effects in the intercept. Model
\eqref{eq:em02} 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 \eqref{eq:em02} nests Model
\eqref{eq:em01} as a special case when $\rho_{it}$ and $\beta_{it}$ are
constant across $i$ and $t$.
As discussed in Section \ref{Section.Introduction}, 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 \eqref{eq:em02} 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{tab:M1} presents the estimation results.
Panel A reports results for Model \eqref{eq:em01} 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 \eqref{eq:em02}, 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{tab:M1} 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{fig:tauA}, we plot the distributions of the estimated
individual effects $\alpha_{i}$ and time effects $\gamma_{t}$ from Model
\eqref{eq:em02}, with corresponding summary statistics reported in Table
\ref{tab:M2d} (supplementary appendix).\footnote{\label{footnote.oy2020}We
acknowledge that the estimated densities of $\alpha_{i}$ and $\gamma_{t}$ may
themselves be affected by the incidental-parameter problem; see
\citet{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{fig:tauB} 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.
\setlength{\tabcolsep}{10pt}
\begin{table}[H]
\centering\begin{threeparttable}
\caption{Single Mother Labor Force Participation: Estimation Results \protect
\label{tab:M1}}
\begin{tabular}{V{\linewidth}V{\linewidth}V{\linewidth}V{\linewidth
}V{\linewidth}}
\hline\hline& \multicolumn{2}{c}{Number of Children} & \multicolumn{2}
{c}{Lagged Participation}\tabularnewline\hline
& Uncorrected & Corrected & Uncorrected & Corrected\tabularnewline
\hline\multicolumn{5}{V{\linewidth}}{\vspace{1bp}
\emph{Panel A. Homogeneous Slope}
\ }\tabularnewline
Model parameter & -0.923 & -0.851 & 0.589 & 0.540\tabularnewline
& (0.232) & (0.228) & (0.108) & (0.107)\tabularnewline\hline\multicolumn
{5}{V{\linewidth}}{\vspace{1bp}
\emph{Panel B. Heterogeneous Slope}
\ }\tabularnewline
Model parameter & 0.704 & 0.436 & 0.191 & 0.253\tabularnewline
& (0.141) & (0.144) & (0.106) & (0.105)\tabularnewline\hline\hline
\end{tabular}
\vspace{1mm}
\begin{tablenotes}
\footnotesize\item\emph{Notes:}
Standard errors in parentheses. Panel A reports the estimation results of the dynamic homogeneous-slope model. Panel B reports the estimation results of the dynamic heterogeneous-slope model. The ``uncorrected'' and the ``corrected'' columns show respectively the MLE and the bias-corrected parameter estimates. Data source: PSID 1987--1997.
\end{tablenotes}
\end{threeparttable}
\end{table}
\begin{center}
\begin{figure}[H]
\caption
{Estimated Heterogeneous Effects of Children on Single Mother Labor Force Participation}
\label{figure2}
\centering\begin{subfigure}[b]{0.48\textwidth}
\centering\includegraphics[width=\textwidth]{taudist_bw.eps}
\caption{\small Distribution of Estimated Effects}
\label{fig:tauA}
\end{subfigure}
\hfill\begin{subfigure}[b]{0.48\textwidth}
\centering\includegraphics[width=\textwidth]{tautrend_bw.eps}
\caption{\small Temporal Evolution of Estimated Effects}
\label{fig:tauB}
\end{subfigure}
\vspace{1mm}
\begin{flushleft}
\begin{justify}
\begin{footnotesize}
\noindent\textit{Notes}
: Panel (a) shows the distributions of the estimated individual effects $\alpha
_i$ (top) and time effects $\gamma
_t$ (bottom) from the dynamic heterogeneous-slope model. Panel (b) plots the median (solid), interquartile range (25th\textendash
{}75th, dashed), and decile range (10th\textendash{}
90th, dotted) of the slope coefficients $\beta_{it}$ over time. Each $\beta
_{it}$ is constructed as the sum of the bias-corrected average effect $\beta
$ and the estimated individual and time effects $\alpha_i$ and $\gamma_t$.
\end{footnotesize}
\end{justify}
\end{flushleft}
\end{figure}
\end{center}
\section{Final Remarks}
\label{Section.Conclusion}
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.
\section*{Acknowledgement}
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.