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.
98,154 characters
Bootstrap Inference in Nonlinear Panel Data Models with Interactive Fixed Effects
\title{
Bootstrap Inference in Nonlinear Panel Data Models with Interactive Fixed Effects\thanks{
We are grateful to participants at the KU Leuven Applied Micro PhD Workshop (2024), the Netherlands Econometric Study Group (2025), the International Association for Applied Econometrics Conference (2025), the International Panel Data Conference (2025), and the Asian Summer School in Econometrics and Statistics (2025).
We thank Andrei Zeleneev and Weisheng Zhang for sharing the codes from \citet{zeleneev2026tractable}. Haoyuan Xu acknowledges funding from the Research Foundation--Flanders (FWO) fellowship no.~1179925N. Wei Miao acknowledges financial support from the China Scholarship Council (grant no.~202306130036). Jad Beyhum acknowledges financial support from the Research Fund KU Leuven (grant STG/23/014).
}
}
\author{
Haoyuan Xu$^\dag$,
Wei Miao$^\ddag$,
Geert Dhaene$^\dag$,
Jad Beyhum$^\dag$
\\[6pt]
$^\dag$Department of Economics, KU Leuven\\
$^\ddag$ORSTAT, KU Leuven
}
\date{\today}
\maketitle
\vspace{-2em}
\begin{abstract}
\noindent
The maximum likelihood estimator in nonlinear panel data models with interactive fixed effects is biased. Several bias correction methods, such as analytical and jackknife approaches, have been proposed to enable valid inference. This paper shows that the parametric bootstrap also enables valid inference in such models. In particular, we show that the parametric bootstrap replicates the asymptotic distribution of the maximum likelihood estimator. Therefore, it yields asymptotically unbiased estimates and confidence sets with asymptotically correct coverage. We also propose a transformation-based bootstrap confidence interval that delivers improved finite-sample performance. Simulation results support the theoretical findings. Finally, we apply the proposed method to examine technological and product market spillover effects on firms' innovation behavior.
\end{abstract}
\textbf{Keywords:} Panel data, interactive fixed effects, incidental parameter bias, bootstrap inference, skewness correction.
\clearpage
\section{Introduction}
The importance of unobserved heterogeneity for economic analysis and the design of effective public policies has long been recognized by economists and policymakers. Interactive fixed effects models provide a flexible and powerful framework for capturing such heterogeneity. Unlike classical one-way fixed-effect models \citep{chamberlain1984panel}, which rely on agent-specific scalar parameters to absorb unobserved heterogeneity, interactive fixed-effect models allow for multidimensional latent components, thereby accommodating substantially richer forms of unobserved heterogeneity. Moreover, these models naturally incorporate aggregate time-varying shocks through a factor structure, a feature that is particularly important in macroeconomic and financial applications \citep{stock2002forecasting}.
Since the seminal contributions of \citet{pesaran2006estimation} and \citet{Bai2009}, a large body of research has developed the theoretical foundations of linear panel data models with interactive fixed effects. These advances have, in turn, facilitated the widespread empirical adoption of such models. Nevertheless, linear specifications are often inadequate in applied work. In many empirical settings, the outcome variable is discrete, rendering linear models potentially misspecified. Recent work, such as \citet{mao2025adaptive}, demonstrates that nonlinear factor models can provide a more appropriate framework in a wide range of empirical applications.
In contrast to the extensive literature on linear factor models \citep{Bai2009,MoonWeidner2015,MoonWeidner2017}, the theoretical foundations of nonlinear panel data models with interactive fixed effects are less developed. A notable recent contribution is \citet{ChenFVWeidner2021}, who study the maximum likelihood estimator (MLE) in nonlinear panel models with interactive fixed effects. Their approach treats the individual-specific and time-specific effects as fixed parameters and jointly estimates these nuisance parameters together with the common parameters of interest. This fixed-effect framework is attractive because it imposes no restrictions on the joint distribution of the unobserved effects and the covariates, thereby reducing the risk of model misspecification.
\citet{ChenFVWeidner2021} derive the asymptotic distribution of the MLE of the common parameters and show that it is asymptotically biased due to the incidental parameter problem. Importantly, this source of bias differs fundamentally from that encountered in linear factor models estimated by least squares, where bias in the slope parameters arises from cross-sectional and time-series dependence in the error terms \citep{Bai2009}. To address the incidental parameter bias, they propose analytical and jackknife corrections that eliminate the leading bias term and enable asymptotically valid inference. Relatedly, \citet{GaoLiuPengYan2023} study the MLE in binary choice models with interactive fixed effects and heterogeneous slope parameters.
Although analytical and jackknife bias corrections are asymptotically valid, their finite-sample performance may be limited in empirically relevant settings with moderate sample sizes and short time dimensions, as commonly encountered in microeconomic panel datasets. Analytical bias corrections rely on estimating the asymptotic bias and subtracting the estimated bias from the original estimator. However, as shown in \citet{ChenFVWeidner2021}, the asymptotic bias is a highly intricate expression involving up to third-order derivatives of the likelihood function, which may be numerically unstable. The split-panel jackknife method is comparatively straightforward to implement but can suffer from efficiency losses in finite samples \citep{hahn2024efficient}. Moreover, analytical and jackknife corrections rely on normal approximations to construct confidence intervals, which may be inaccurate.
Recently, \citet{HigginsJochmans2024} show that the parametric bootstrap delivers asymptotically valid confidence intervals without explicit bias correction in one-way fixed-effect nonlinear panel data models, and that bootstrap-based inference can exhibit superior finite-sample performance when the time dimension is short. In this paper, we extend the parametric bootstrap approach to nonlinear panel data models with interactive fixed effects. Unlike analytical bias-correction methods, the parametric bootstrap avoids the explicit calculation of the asymptotic bias by approximating it through parametric bootstrap simulations. An important feature of the parametric bootstrap, emphasized by \citet{HigginsJochmans2024}, is that it yields asymptotically valid inference even without preliminary bias correction. Confidence intervals are constructed directly from the quantiles of the bootstrap distribution rather than relying on normal approximations, which can improve accuracy in finite samples. We show that these properties of the parametric bootstrap also hold in nonlinear panel data models with interactive fixed effects.
Simulations corroborate our theoretical results. Furthermore, simulations indicate that bootstrap-based confidence intervals tend to be conservative in finite samples, with coverage probabilities often exceeding the nominal level and thus reflecting a systematic overestimation of uncertainty. This finding is in line with \citet{HigginsJochmans2024}, who document similar behavior in one-way fixed-effect models. They also note that iterating the bootstrap procedure can improve coverage accuracy, but at a substantial computational cost.
To mitigate this issue, we follow the bootstrap literature on skewness correction \citep{hall_removal_1992} and apply a monotone transformation to the estimator prior to constructing confidence intervals. The transformation reduces skewness in the bootstrap distribution, and intervals are then obtained from the transformed scale and mapped back. By the bootstrap delta method \citep{van1998asymptotic}, this procedure preserves asymptotic validity while improving finite-sample performance. Simulations show that the resulting intervals are substantially shorter and achieve more accurate coverage.
All bias correction methods considered in this literature rely on the MLE, i.e., the global maximizer of the loglikelihood function, being obtainable. However, in practice, the loglikelihood function in nonlinear factor models is generally non-concave, implying that standard optimization algorithms may converge to a local instead of the global maximizer. \citet{ChenFVWeidner2021} propose an EM-type algorithm, which guarantees convergence only to a local maximizer. They suggest using multiple starting values to increase the probability of reaching the global maximizer; nevertheless, this strategy does not ensure global maximization and can substantially increase the computational burden. More recently, \citet{zeleneev2026tractable} propose a computationally efficient two-step estimator based on nuclear-norm penalization, which is shown to attain the global maximizer with probability one and is asymptotically equivalent to the MLE. In this paper, we adopt their two-step estimator in conjunction with the bootstrap procedure, thereby ensuring global maximization at low computational cost and preserving asymptotically valid inference.
\\[0.5em]
\noindent \textbf{Related literature}
\noindent
This paper contributes to the literature on large-$T$ bias correction in nonlinear panel data models with fixed effects. Since the seminal work of \citet{HahnNewey2004}, many studies have aimed to correct the first-order incidental parameter bias in a large-$T$ framework. This includes research on one-way fixed effects panel models \citep{fernandez2009fixed,HahnKuersteiner2011,DhaeneJochmans2015,arellano2016likelihood,HigginsJochmans2024}, additive two-way fixed effects panel models \citep{FVWeidner2016}, and network models \citep{Graham2017,dzemski2019empirical,yan2019statistical,HUGHES2026106130}; see \citet{FVWeidner2018} for a review. More recently, some studies have begun to address second-order bias \citep{DHAENE2021227, schumann2023second}
as well as higher-order bias correction \citep{dhaene2023approximate,bonhomme2024neyman}.
We also contribute to the literature on fixed-effect estimation in models with interactive fixed effects. Linear panel models with interactive fixed effects have been extensively studied, including the determination of the number of factors \citep{bai2002determining, onatski2009testing, ahn2013eigenvalue}, the estimation of common parameters \citep{Bai2009, MoonWeidner2015, MoonWeidner2017, BeyhumGautier2019, BeyhumGautier2022}, inference in the presence of weak factors \citep{Onatski2012, BaiNg2023, ChoiYuan2025, jiang2025biascorrectionfactoraugmentedregression, ArmstrongWeidnerZeleneev2025}, and extensions to network models \citep{sassi2024linear}; see \citet{BaiWang2016} for a comprehensive review. The theoretical development for nonlinear factor models remains comparatively limited. Notable contributions in this area include \citet{WANG2022180}, who studies the fixed-effect MLE in nonlinear pure factor models; \citet{GaoLiuPengYan2023}, who analyze binary choice models with interactive fixed effects and heterogeneous slope parameters and propose an information criterion for selecting the number of factors; and \citet{ChenFVWeidner2021}, who address bias correction in general nonlinear factor models and introduce an eigenvalue ratio test for factor number selection. More recently, \citet{zeleneev2026tractable} and \citet{yao2025low} propose a nuclear-norm penalized estimator that provides computational guarantees, facilitating reliable estimation.
\citet{boneva2017discrete} and \citet{chen2025common} generalize the CCE estimator to nonlinear settings.
Our work is also closely related to bootstrap-based bias correction and inference methods in the panel data literature. \citet{GONCALVES2015407} study bootstrap inference for linear dynamic panel models with individual fixed effects. \citet{KimSun2016} and \citet{HigginsJochmans2024} adopt parametric bootstrap methods for nonlinear panel models with one-way fixed effects. \citet{gonccalves2011moving} and \citet{higgins2025inference} introduce moving-block bootstrap procedures for dynamic panel data models. \citet{LI2024105684} employ bootstrap methods for inference for treatment effects in interactive fixed-effect panel models. \citet{cavaliere2024bootstrap} develop a general framework for bootstrap inference in the presence of bias, showing that even when the bias cannot be consistently estimated, appropriately designed bootstrap methods can provide valid inference. Finally, our work is also related to the literature that uses transformation-based methods to improve bootstrap confidence intervals \citep{efron1987better,konishi_normalizing_1991,hall_removal_1992}.
\section{Methodology}
\subsection{Model setup}
We consider nonlinear panel data models with interactive fixed effects. The model specification is
\begin{equation}\label{model 2.1}
Y_{it}|X_{it},\beta_{0},\alpha_0,\gamma_0\sim f\left(\cdot|X_{it}^{\prime}\beta_0+\alpha_{i0}^{\prime}\gamma_{0t}\right)
\end{equation}
for $i=1,\ldots,N$ and $t=1,\ldots,T$, where we observe
$\{Y,X\}=\{Y_{it},X_{it}\}_{i=1\ldots,N;t=1,\ldots,T}$, $Y_{it}$ is the scalar response variable, $X_{it}$ is a $d_x$-dimensional covariate vector, the function $f$ is a known density with respect to some dominating measure, $\beta_0$ is the common parameter vector, and $\alpha_{i0}$ and $\gamma_{0t}$ are unobserved $d_f$-dimensional individual and time effects that appear in the model through a factor structure. We treat the number of factors $d_f$ as known.\footnote{In practice, we can estimate $d_f$ consistently; see \citet{ChenFVWeidner2021} and \citet{GaoLiuPengYan2023}.}
We aim to estimate $\beta_0$, treating
$\alpha_{0} = (\alpha_{1,0}, \ldots, \alpha_{N,0})$ and
$\gamma_{0} = (\gamma_{1,0}, \ldots, \gamma_{T,0})$ as nuisance parameter matrices.
We make no assumptions on the joint distribution of $(X,\alpha_{0},\gamma_{0})$.
This flexibility is important since it reduces the likelihood
of model misspecification. Unlike the individual fixed-effect model \citep{HahnNewey2004, HahnKuersteiner2011}, which only allows the individual effect
to be a time-invariant scalar, the interactive fixed-effect model
accommodates multiplicative fixed effects.
The additive two-way fixed-effect model \citep{FVWeidner2016}, which includes time effects
to capture aggregate shocks---a common feature in economic applications---can, in fact, be viewed as a special case of model \eqref{model 2.1}.
Interactive fixed effects allow for a much richer structure
of unobserved heterogeneity \citep{Freyberger2018}.
\subsection{Interactive fixed-effect estimator}
Let $\mathcal{B}\subset \mathbb{R}^{d_x}$, $\mathcal{A}\subset \mathbb{R}^{d_f}$, and $\mathcal{G}\subset \mathbb{R}^{d_f}$ denote the parameter spaces of $\beta_0$, $\alpha_{i0}$, and $\gamma_{t0}$, respectively, and let $\Theta = \mathcal{B} \times \mathcal{A}^N \times \mathcal{G}^T$ denote the full parameter space of $\theta_0 = (\beta_0, \alpha_0, \gamma_0)$.
The MLE of $\theta_0$ is
\begin{equation}\label{mle problem}
\hat\theta=
(\hat{\beta},\hat{\alpha},\hat{\gamma})
=\underset{(\beta, \alpha, \gamma) \in \mathcal{B} \times \mathcal{A}^N \times \mathcal{G}^T}{\arg\max}\,\mathcal{L}(\beta, \alpha, \gamma; Y, X),
\end{equation}
where $\mathcal{L}(\beta, \alpha, \gamma; Y, X)$ is the loglikelihood function,
\begin{equation}
\mathcal{L}(\beta, \alpha, \gamma; Y, X)
= \sum_{i=1}^N \sum_{t=1}^T
\log f\!\left(Y_{it} \mid X_{it}^{\prime}\beta + \alpha_{i}^{\prime}\gamma_{t}\right).
\label{eq:loglikelihood}
\end{equation}
Since the dimension of the nuisance parameter matrices $\alpha_0$ and $\gamma_0$ increase with $N$ and $T$,
the estimation of $\beta_0$ depends on a large number of nuisance parameter estimates, which can induce substantial bias in $\hat\beta$. Even as $N,T \to \infty$, this bias persists in the limiting distribution of $\hat{\beta}$ and thus affects inference. Our main contribution is to show that the parametric bootstrap eliminates the asymptotic bias and delivers asymptotically valid confidence intervals.
As in the linear factor model \citep{Bai2009}, the unobserved effects $\alpha_{i0}$ and $\gamma_{t0}$ can only be identified up to a rotational indeterminacy since, for any nonsingular $d_f\times d_f$ matrix $A$, we have
$\alpha_{i0}'\gamma_{t0}=(A'\alpha_{i0})'(A^{-1}\gamma_{t0})$ for all $i$ and $t$.
Different normalizations do not affect the estimation of $\beta_0$ nor the estimation of average partial effects (APEs). For further discussion on normalization in factor models, see
\citet{BaiWang2016}.
We are also interested in APEs. For a chosen function $\mu_{it}(\beta,\alpha_{i},\gamma_{t})$ that characterizes the causal effect of interest to the researcher, we define the APE as
\begin{align}\label{APE}
\varDelta(\theta_0) = \frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\mathbb{E}_{\theta_0}\!\left[\mu_{it}(\beta_0,\alpha_{i0},\gamma_{t0})\right],
\end{align}
where $\mathbb{E}_{\theta_0}$ denotes the conditional expectation, given the true values of the unobserved effects.\footnote{We focus on conditional APEs, given the realizations of the unobserved effects. Studying unconditional APEs requires additional assumptions on the distribution of $\{\alpha_{i0}\}_{i=1,\ldots,N}$ and $\{\gamma_{t0}\}_{t=1,\ldots,T}$, such as i.i.d.\ sampling or weak dependence. For further discussion, see \citet{FVWeidner2016}.}
For example, if $Y_{it}$ is binary and we consider the effect of a binary covariate $X_{it,k}$ on the probability that $Y_{it}=1$, we set $\mu_{it}(\beta,\alpha_{i},\gamma_{t})= f(1|Z_{it}+\beta_k)-f(1|Z_{it})$, where
$Z_{it}=X_{it,-k}^{\prime}\beta_{-k}+\alpha_i'\gamma_t$, and $X_{it,-k}$ and $\beta_{-k}$ denote $X_{it}$ and $\beta$, respectively, with their $k$-th element removed.
If $Y_{it}$ is binary and $X_{it,k}$ is continuous, the average marginal effect of interest is obtained by setting $\mu_{it}(\beta,\alpha_{i},\gamma_{t}) = \beta_k\, f(1|X_{it}^{\prime}\beta+\alpha_{i}^{\prime}\gamma_{t})$.
A natural estimator of $\varDelta(\theta_0)$ is the MLE,
\begin{align}\label{plugin-ape}
\varDelta(\hat{\theta})
= \frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}
\mu_{it}(\hat{\beta},\hat{\alpha}_{i},\hat{\gamma}_{t}),
\end{align}
which inherits the bias of $(\hat{\beta},\hat{\alpha},\hat{\gamma})$.
\citet{ChenFVWeidner2021} study the MLE of model~\eqref{model 2.1} and propose a PCA-based algorithm to compute the estimator. They also provide analytical and jackknife corrections for the incidental parameter bias of
$\hat\beta$ and $\varDelta(\hat{\theta})$, thereby restoring asymptotically valid inference.
However, when $T$ is small, the analytical and jackknife bias corrections
may not always perform well; see our simulation results in Section 5. Recently, \citet{HigginsJochmans2024}
showed that the parametric bootstrap corrects the first-order
bias and provides asymptotically valid confidence intervals for $\beta_0$ in nonlinear
panel models with individual effects. Their parametric bootstrap confidence intervals
do not rely on normal approximations, and perform better compared to analytical corrections.
In this paper, we extend this idea to interactive fixed-effect models
and show in simulations that the parametric bootstrap delivers better small-$T$ performance in nonlinear interactive fixed-effect panel models compared with existing bias-correction methods.
We further propose an improved bootstrap confidence interval based on a monotone transformation of the estimator.
\subsection{Bootstrap inference}
In this section, we describe the parametric bootstrap and informally explain why it yields a bias-corrected estimate of $\beta_0$ and asymptotically correct confidence intervals for $\beta_0$.
We begin with the data-generating process (DGP) used to generate
$B$ bootstrap samples, each denoted by \(Y^*= \{Y^*_{it}\}_{i=1\ldots,N;t=1,\ldots,T} \).
Under the assumption that the covariates $X$ are strictly exogenous, the parametric bootstrap DGP is obtained
from model \eqref{model 2.1} by plugging in the MLE \( \hat{\theta} \) and keeping \( X \) fixed:
\begin{equation}\label{bootstrap DGP}
Y^*_{it}|X_{it},\hat{\beta},\hat{\alpha},\hat{\gamma}\sim f\left(\cdot|X_{it}^{\prime}\hat{\beta}+\hat{\alpha_{i}}^{\prime}\hat{\gamma_{t}}\right).
\end{equation}
The bootstrap DGP in \eqref{bootstrap DGP} is the same as in \citet{HigginsJochmans2024} except that we have interactive fixed effects instead of one-way individual effects.
For each bootstrap dataset \( \{Y^*, X\} \), with $Y^*$ independently drawn according to (\ref{bootstrap DGP}), we compute the MLE
\[
\hat{\theta}^*
=(\hat\beta^*, \hat{\alpha}^*, \hat{\gamma}^*) =
\underset{(\beta, \alpha, \gamma) \in \mathcal{B} \times \mathcal{A}^N \times \mathcal{G}^T}{\arg\max}
\mathcal{L}( \beta, \alpha, \gamma;Y^*, X).
\]
Since \( \hat{\theta} \) is a consistent estimator---in the sense that, as $N,T\to\infty$, all of $\hat{\beta},\hat{\alpha}_i,\hat{\gamma}_t$, for every $i$ and $t$, are consistent for $\beta_0,\alpha_{i0},\gamma_{t0}$---\( \hat{\theta} \)
should be close to \( \theta_0 \) with high probability when $N$ and $T$ are sufficiently large. Hence, under certain smoothness conditions, the parametric bootstrap DGP in \eqref{bootstrap DGP} will also be close to the true DGP, that is, to the DGP in \eqref{model 2.1}, which has $\theta=\theta_0$. Furthermore, the bootstrap datasets \( \{Y^*, X\} \) and the corresponding MLEs generated by these two very similar DGPs are expected to exhibit similar properties.
Under suitable regularity conditions, $\hat{\beta}^*$ admits an expansion analogous to that of $\hat{\beta}$. Moreover, the first-order term in the expansion of $\hat{\beta}^*$ converges in probability to the corresponding term in the expansion of $\hat{\beta}$. This implies that the bias
$\mathbb{E}_{\theta_0}(\hat{\beta}) - {\beta_0}$ can be accurately appoximated by
$\mathbb{E}^*(\hat{\beta}^*) - \hat{\beta}$,
where $\mathbb{E}^*$ denotes the expectation under the parametric bootstrap DGP \eqref{bootstrap DGP}.
In turn, $\mathbb{E}^*(\hat{\beta}^*) - \hat{\beta}$ can be approximated arbitrarily well by averaging over a large number, $B$, of parametric bootstrap replications. This average can be used to eliminate the first-order bias of $\hat{\beta}$, yielding the parametric bootstrap bias-corrected estimator
\begin{equation*}
\hat{\beta}_{\rm boot}=
\hat{\beta} - \left(\frac{1}{B}\sum_{b=1}^B \hat{\beta}_b^* - \hat{\beta} \right)
= 2 \hat\beta - \frac{1}{B}\sum_{b=1}^B \hat{\beta}_b^*,
\end{equation*}
where $\hat{\beta}_b^*$ is $\hat{\beta}^*$ computed from the $b$-th bootstrap dataset. Since the first-order bias is eliminated,
\begin{equation*}
\mathbb{E}_{\theta_0}\!(
\hat{\beta}_{\rm boot}
) - \beta_0
=
o\!\left(
\frac{1}{N} + \frac{1}{T}
\right).
\end{equation*}
Note that, like other bias correction methods, the parametric bootstrap bias correction only provides asymptotic guarantees of bias reduction. With finite $N$ and $T$, it may amplify higher-order bias terms and, therefore, amplify the bias. In Section~4, we examine the finite-sample performance of the parametric bootstrap through simulations.
Another important observation is that the parametric bootstrap replicates the distribution of \( \hat{\beta} \) with a vanishingly small error.
For every \( a< b \in \mathbb{R}^{d_x} \),
\begin{equation*}
\mathbb{P}_{\theta_0}(a \leq \sqrt{NT}(\hat{\beta}^* - \hat{\beta}) \leq b) - \mathbb{P}_{\theta_0}(a \leq \sqrt{NT}(\hat{\beta} - \beta_0) \leq b) \to 0,
\end{equation*}
where $\mathbb{P}_{\theta_0}(\cdot)$ denotes probabilities taken under the DGP \eqref{model 2.1}.
As shown in \citet{HigginsJochmans2024}, this yields asymptotically valid confidence intervals for $\beta_0$ without the need for bias correction.
To illustrate, we consider the simple case where $d_x= 1$. In case $d_x > 1$, a valid confidence interval for any linear combination $c^\prime \beta_0$ can be constructed similarly.
Let $F^*$ denote the conditional distribution function of $\hat{\beta}^* - \hat{\beta}$ under DGP~\eqref{bootstrap DGP} given the original sample; $F^*$ can be estimated to arbitrary precision by the empirical distribution of the parametric bootstrap replicates $ \hat{\beta}_b^*$ by setting $B$ large enough. In the next section, we formally show that \( F^* \) weakly converges to the limit distribution of \( \hat{\beta} -\beta_0\). Similarly, the quantile function corresponding to \( F^* \), defined as
\[
Q^*(\alpha) := \inf\{a: F^*(a) \geq \alpha\}, \qquad \alpha\in\mathbb{R},
\]
also converges to the quantile function of the limiting distribution of \( \hat{\beta}-\beta_0 \). Hence, we can simply use \( Q^*(\alpha) \) as the estimator of the \( \alpha \)-quantile of the distribution of \( \hat{\beta} - \beta_0 \), leading to an asymptotically valid confidence interval: for every $\alpha\in(0,1)$,
\begin{equation}\label{CI}
\mathbb{P}_{\theta_0}\left(\beta_0 \in \left[ \hat{\beta} - Q^*(1-\alpha/2), \hat{\beta} - Q^*(\alpha/2) \right]\right) \to 1-\alpha,
\end{equation}
where the interval in brackets is a $(1-\alpha)$-level confidence interval for a scalar \( \beta_0 \). The interval can also be written as
\begin{equation*}
\left[ 2\hat{\beta} - Q^*_{\hat\beta^*}(1-\alpha/2), 2\hat{\beta} - Q^*_{\hat\beta^*}(\alpha/2) \right],
\end{equation*}
where $Q^*_{\hat\beta^*}$ is the conditional quantile function of $\hat\beta^*$ given the original sample.
Other methods to construct confidence intervals often rely on the asymptotic normality of \( \hat{\beta} \). They require a consistent estimator, \( \hat{V} \), of the asymptotic variance $V$ of \( \hat{\beta} \),
and an asymptotically unbiased estimator $\hat{\beta}^{bc}$ with the same asymptotic variance $V$. This gives,
for example, \( [\hat{\beta}^{bc} - 1.96 \sqrt{\hat{V}}, \hat{\beta}^{bc} + 1.96 \sqrt{\hat{V}}] \) as a 95\% confidence interval for a scalar $\beta_0$. The parametric bootstrap
does not require a bias-corrected estimator of $\beta_0$ nor a consistent estimator of $V$. Instead, it simply uses the quantile function of \( \hat{\beta}^* - \hat{\beta} \) and it does not rely on the Gaussian approximation to the distribution
of $\hat{\beta}$.
While the parametric bootstrap delivers asymptotically valid inference, its finite-sample performance may be affected by skewness in the distribution of $\hat{\beta}$, especially when the time dimension is small. As discussed above, the bootstrap procedure reproduces and may even amplify this skewness, which can result in asymmetric confidence intervals with overly wide bounds and thus lead to conservative inference.
To mitigate this issue, we propose to apply a monotone transformation to the bootstrap estimator prior to constructing confidence intervals. Let $\varphi(\cdot)$ be a strictly increasing transformation. We then construct the transformed bootstrap confidence interval as
\begin{equation}\label{transform CI}
\left[
\varphi^{-1}\!\left(2\varphi(\hat{\beta}) - Q^*_{\varphi(\hat{\beta}^*)}(1-\alpha/2)\right), \;
\varphi^{-1}\!\left(2\varphi(\hat{\beta}) - Q^*_{\varphi(\hat{\beta}^*)}(\alpha/2)\right)
\right],
\end{equation}
where $Q^*_{\varphi(\hat{\beta}^*)}$ is the conditional quantile function of $\varphi(\hat{\beta}^*)$ given the original sample.
Since $\varphi(\cdot)$ is strictly increasing, this transformation preserves the ordering of the estimator and thus maintains asymptotic validity. At the same time, by reducing skewness in the bootstrap distribution, it yields more symmetric and tighter confidence intervals, thereby improving finite-sample performance.
An alternative approach to account for skewness is the bias-corrected and accelerated (BCa) bootstrap of \citet{efron1987better}. However, this method is not well suited to our setting. The BCa procedure adjusts confidence intervals based on an estimate of bias and skewness inferred from the bootstrap distribution. When the estimator exhibits non-negligible asymptotic bias, as in our case, the bootstrap distribution inherits this bias and is effectively shifted relative to the true sampling distribution.
As a result, the BCa correction may confound bias with skewness and mischaracterize the shape of the distribution. For example, if the asymptotic bias shifts the distribution to the right, the BCa method may interpret this shift as evidence of skewness and consequently adjust the confidence interval in the opposite direction, leading to an excessive leftward shift. This can distort coverage and produce misleading inference. In contrast, our transformation-based approach directly addresses the asymmetry of the distribution while leaving its location unaffected, thereby providing a more reliable correction in the presence of skewness.
By analogy to (\ref{CI}) and (\ref{transform CI}), we can also construct a parametric bootstrap confidence interval for $\varDelta({\theta_0})$. By plugging $\hat{\theta}^*$ into the APE, we obtain
$$\ \varDelta(\hat{\theta}^*) = \frac{1}{NT} \sum_{i=1}^N \sum_{t=1}^T \mu(X_{it},\hat{\beta}^*,\hat{\alpha}_i^*,\hat{\gamma}_t^*),$$
which allows us to construct a bootstrap confidence interval for
$\varDelta({\theta_0})$, since the distribution of \( \sqrt{NT}(\varDelta(\hat{\theta}^*) - \varDelta(\hat{\theta})) \) approximates that of \( \sqrt{NT}(\varDelta(\hat{\theta}) - \varDelta(\theta_0)) \) as \( N,T \) becomes large.
\section{Asymptotic theory}
In this section, we study the asymptotic properties of the parametric bootstrap in an asymptotic regime where $N$ and $T$ go to infinity with $T/N \to \kappa \in (0,\infty)$. For simplicity, we focus on the original bootstrap estimator. The results for the transformed bootstrap can be obtained by combining bootstrap consistency for the transformed estimator, established via the bootstrap delta method\footnote{This requires that $\varphi$ be continuously differentiable at the true parameter $\beta_0$ with $\varphi'(\beta_0)\neq 0$, so that the delta method and its bootstrap analogue apply; see, e.g., Chapter 23 of \citet{van1998asymptotic}.}, with the fact that, under monotonicity, the coverage of the confidence interval for $\beta_0$ is equivalent to that for $\varphi(\beta_0)$.
Our proofs closely follow those of \citet{KimSun2016} and \citet{HigginsJochmans2024}. We first extend the asymptotic expansion for the general M-estimator established in \citet{FVWeidner2016} to ensure that it holds not only at the true parameter value but also uniformly in a neighborhood of the true value. Achieving this requires that the assumptions in \citet{FVWeidner2016} hold uniformly in a small neighbourhood of $\theta_{0}$. Similar uniformity assumptions are imposed in \citet{HigginsJochmans2024} for the one-way fixed-effect model.
Then we verify that nonlinear panel models with interactive fixed effects satisfy these high-level conditions, thereby establishing the validity of uniform expansions for the MLEs of $\beta_0$ and $\varDelta(\theta_0)$. These uniformity results, combined with the consistency of the MLE and a lemma from \citet{Andrews2005}, show the validity of the parametric bootstrap to approximate the distribution of the MLE.
We treat the unobserved effects as fixed parameters. Alternatively, our analysis can be interpreted as being conditional on the realization of the unobserved effects, as in \citet{FVWeidner2016} and \citet{ChenFVWeidner2021}, provided that the assumptions are modified by conditioning on $(\alpha_0, \gamma_0)$. See Remark~1 of \citet{HahnKuersteiner2011} for more discussion.
Since we generate the bootstrap sample $\{Y^{*}, X\}$ based on model~\eqref{bootstrap DGP}, the covariates are kept fixed when constructing the bootstrap sample. Throughout the paper, we let $\mathbb{E}_{\theta}$ and $\mathbb{P}_{\theta}$ denote the expectation and probability operators under the bootstrap data-generating process.\footnote{Strictly speaking, since the joint distribution of the covariates $X$ and the unobserved effects $\alpha_0,\gamma_0$ is unrestricted, the bootstrap data-generating process is indexed by $(\theta,\alpha_0,\gamma_0)$. Accordingly, the corresponding expectation and probability operators should be written as $\mathbb{E}_{\theta,\alpha_0,\gamma_0}$ and $\mathbb{P}_{\theta,\alpha_0,\gamma_0}$. For notational simplicity, we suppress the dependence on $\alpha_0$ and $\gamma_0$.}
We make the following assumptions. They are similar to, but slightly stronger than, Assumption~1 in \citet{ChenFVWeidner2021}, as some of our assumptions are required to hold uniformly over a neighborhood $\Theta_0$ of the true parameter value.
\begin{assumption}[]\label{assu:parameterspace}
The parameter space $\Theta$ is a compact set. Furthermore, there
exists a constant $\varepsilon > 0$ and an open neighborhood $\Theta_0 \subset \Theta$ that contains $\theta_0$ and satisfies $d(\theta, \theta_0) < \varepsilon$
for every $\theta \in \Theta_0$, where
\begin{align}\label{metric}
d(\theta, \theta_0)
= \big\Vert
\left(
\beta^{\prime},\, \mathrm{vec}(\alpha)^{\prime},\, \mathrm{vec}(\gamma)^{\prime}
\right)^{\prime}
-
\left(
\beta_0^{\prime},\, \mathrm{vec}(\alpha_0)^{\prime},\, \mathrm{vec}(\gamma_0)^{\prime}
\right)^{\prime}
\big\Vert
\end{align}
and $\Vert \cdot \Vert$ denotes the Euclidean norm.
\end{assumption}
The following assumption involves a mixing condition, which we use to bound the moments. Let $\mathcal{A}_t^i$ and $\mathcal{B}_t^i$ be the $\sigma$-algebras generated by $\{ X_{i1}, Y_{i1}. \ldots, X_{it}, Y_{it} \}$
and $\{ X_{it}, Y_{it}, \allowbreak X_{i,t+1}, Y_{i,t+1} , \ldots \}$, respectively, and define the
$\alpha$-mixing coefficients
\[
\alpha_i(m,\theta) := \sup_t \sup_{A \in \mathcal{A}_t^i,\, B \in \mathcal{B}_{t+m}^i}
\left| \mathbb{P}_\theta(A \cap B) - \mathbb{P}_\theta(A)\mathbb{P}_\theta(B) \right|.
\]
\begin{assumption}[Mixing]\label{assu:correct}
(i) Conditional on \( X \), the data $Y=\{Y_{it}\}_{i=1\ldots,N;t=1,\ldots,T}$
are generated from model \eqref{model 2.1}. (ii) For each $t$, the data $\{X_{it},Y_{it}\}_{i=1,\ldots,N}$ is independent across $i$, and
there exists a constant $\mu>0$ such that, for each $i$, $\{X_{it},Y_{it}\}_{t=1,\ldots,T}$ is $\alpha$-mixing with mixing coefficients $\alpha_i(m,\theta)$ satisfying $\sup_{\theta\in\Theta_0}\max_i \alpha_i(m,\theta)=O(m^{-\mu})$ as $m\to\infty$.
\end{assumption}
We assume that the covariates are strictly exogenous.\footnote{However, our approach could be extended to include lagged dependent variables affecting $Y_{it}$. More general feedback mechanisms---such as cases where $Y_{it}$ affects $X_{is}$ for $s>t$---are excluded. For further discussion, see \citet{ChenFVWeidner2021}.} The conditional independence in Assumption \ref{assu:correct}(i) is commonly adopted in the fixed-effects panel data literature \citep{HahnNewey2004, Bonhomme2012, FVWeidner2016, ChenFVWeidner2021}, as it implies that the conditional distribution of $Y_{it}$ does not depend on $\alpha$, $\gamma$, or $X$ beyond what is captured by the specified model.\footnote{Recently, several papers studying average partial effects in nonlinear panel data models with individual fixed effects relax this assumption \citep{GrahamPowell2012, liu2024identification, BotosaruMuris2024}.}
The mixing condition, also assumed in \citet{FVWeidner2016}, is used to bound the moments and covariances, and to apply the WLLN and CLT. Moreover, we require the mixing coefficients to decay uniformly over $\theta \in \Theta_0$ to uniformly bound the moments, which is necessary to establish uniform convergence of the bias term and uniform asymptotic normality. The uniform decay condition is also adopted in \citet{HigginsJochmans2024}.
\begin{assumption}[Strong factors]\label{ass:strong factor}
There exist positive definite matrices $\Sigma_{\alpha}$ and $\Sigma_{\gamma}$ such that
$\frac{1}{N}\sum_{i=1}^{N} \alpha_{i0}\alpha_{i0}^{\prime} \;\to\; \Sigma_{\alpha}$ and
$\frac{1}{T}\sum_{t=1}^{T} \gamma_{t0}\gamma_{t0}^{\prime} \;\to\; \Sigma_{\gamma}$.
\end{assumption}
Assumption~\ref{ass:strong factor}, commonly referred to as the strong factors assumption, is widely used in panel interactive fixed effects models \citep{Bai2009, MoonWeidner2015, MoonWeidner2017, ChenFVWeidner2021, GaoLiuPengYan2023}.\footnote{For discussion of weak factors, see \citet{Onatski2012, BeyhumGautier2019, BeyhumGautier2022, BaiNg2023, ChoiYuan2025, jiang2025biascorrectionfactoraugmentedregression, ArmstrongWeidnerZeleneev2025}.} In our setting, this assumption is required to establish a uniform convergence rate of the MLE, which is needed to apply the asymptotic expansion in Appendix~\ref{Appendix C}, a uniform version of the expansion in \citet{FVWeidner2016}. A necessary step is to ensure that the matrices $\frac{1}{N}\sum_{i=1}^{N}\alpha_{i}\alpha_{i}^{\prime}$ and $\frac{1}{T}\sum_{t=1}^{T}\gamma_{t}\gamma_{t}^{\prime}$ converge uniformly, over $\theta \in \Theta_{0}$, to positive definite limits. Given the fact that $\Theta_{0}$ is an $\epsilon$-neighborhood of $\theta_0$, the differences
$\frac{1}{N}\sum_{i=1}^{N}\alpha_{i}\alpha_{i}^{\prime}-\frac{1}{N}\sum_{i=1}^{N}\alpha_{i0}\alpha_{i0}^{\prime}$ and
$\frac{1}{N}\sum_{i=1}^{N}\gamma_{i}\gamma_{i}^{\prime}-\frac{1}{N}\sum_{i=1}^{N}\gamma_{i0}\gamma_{i0}^{\prime}$
converge to zero for every $\theta\in\Theta_0$. Therefore, Assumption~\ref{ass:strong factor} is sufficient to establish the required uniform consistency result.
For a formal proof, see Lemma~\ref{lemma:A1}.
Let $X_k$ be the $N \times T$ matrix elements $X_{it,k}$ ($i = 1, \ldots, N$; $t = 1, \ldots, T$).
For any matrix $A$, define the residual-maker matrix
$\mathcal{M}_A = \mathbb{I} - A (A^{\prime} A)^{\dagger} A^{\prime}$
where ${}^\dagger$ denotes the Moore-Penrose inverse and $\mathbb{I}$ is the identity matrix of appropriate dimension. Let $\operatorname{Tr}(\cdot)$ denote the trace of a matrix.
\begin{assumption}[Generalized non-collinearity]\label{Generalized non-collinearity}
There exists a constant $c>0$ such that, for every $\theta\in\Theta_0$ and every $d_x \times T$ matrix $\widetilde\gamma$,
the $d_x \times d_x$ matrix $D(\alpha, \widetilde\gamma)$ with elements
\[
D_{k_1 k_2}(\alpha, \widetilde\gamma)
= (NT)^{-1} \operatorname{Tr}\!\left(
\mathcal{M}_{\alpha'} X_{k_1} \mathcal{M}_{\widetilde\gamma'} X_{k_2}^{\prime}
\right),
\quad k_1, k_2 \in \{1, \ldots, d_x\},
\]
satisfies
$\inf_{\theta \in \Theta_0} \lambda_{\min}\!\left(D(\alpha, \widetilde\gamma)\right) > c$
wpa1, where $\lambda_{\min}(\cdot)$ denotes the smallest eigenvalue.
\end{assumption}
Assumption~\ref{Generalized non-collinearity} requires that the unobserved factors and factor loadings cannot fully capture the variation of any covariate, and that this condition holds not only for the true factors, but also for all $\alpha$ in a neighborhood of $\alpha_0$. It ensures that $\beta_0$ is point-identified.
Let $Z_{it}(\theta)=X_{it}^{\prime}\beta + \alpha_i^{\prime}\gamma_t$ and define
$\ell_{it}(z) = \log f\!(Y_{it}|Z_{it}(\theta)=z)$.
Let $\partial_{z^q} \ell_{it}(z)$ denote the $q$-th derivative of $\ell_{it}(z)$.
When evaluating these functions at $z=Z_{it}(\theta_0)$, we omit the argument for simplicity.
For instance, we write $\partial_{z^{q}} \ell_{it}$ instead of
$\partial_{z^{q}} \ell_{it}(Z_{it}(\theta_0))$.
For $z \in \mathbb{R}^q$ and $\delta > 0$, let
$\mathcal{B}(\delta, z) \subset \mathbb{R}^q$ denote the open ball of radius
$\delta$ centered at $z$, taken with respect to the topology induced by the
Euclidean metric on $\mathbb{R}^q$.
\begin{comment}
Denote the conditional log-likelihood function by
\[
\ell_{it}(Z_{it}): = \log f\!\left(Y_{it}|X_{it}^{\prime}\beta + \alpha_i^{\prime}\gamma_t\right),
\]
where $Z_{it}=X_{it}^{\prime}\beta + \alpha_i^{\prime}\gamma_t$ and let $\partial_{z^q} \ell_{it}(z)$ denote its $q$-th derivative with respect to $z$.
When these functions are evaluated at the true parameter values
$(\beta_{0}, \alpha_{0}, \gamma_{0})$,
e.g, $X_{it}^{\prime}\beta_{0} + \alpha_{i0}^{\prime}\gamma_{0t}$,
we omit the argument $z_{it}$ for simplicity.
For instance, we write $\partial_{z^{q}} \ell_{it}$ instead of
$\partial_{z^{q}} \ell_{it}(X_{it}^{\prime}\beta_{0} + \alpha_{i0}^{\prime}\gamma_{0t})$.
\end{comment}
\begin{assumption}[Smoothness and moments]\label{assu:smooth}
There exists $\delta > 0$ such that for every $i,t$ and every $\theta \in \Theta_0$, the following conditions hold almost surely:
\begin{itemize}
\item[(i)] The function $\ell_{it}(\cdot)$ is four times continuously differentiable on
$\mathcal{B}(\delta,Z_{it}(\theta))$.
\vspace{0.4em}
\item[(ii)] The derivatives of $\ell_{it}(\cdot)$ up to fourth order are
uniformly bounded in absolute value on $\mathcal{B}(\delta,Z_{it}(\theta))$ by a random envelope
$L_{it}(Z_{it}(\theta),\theta) > 0$ satisfying
\[
\sup_{\theta \in \Theta_0}
\max_{i,t}
\mathbb{E}_{\theta}\!\left[L_{it}(Z_{it}(\theta),\theta)^{\,8+\nu} \right] < c,
\]
for some constants $\nu>0$ and $c>0$, uniformly in $N,T$.
\item[(iii)] The regressors $X_{it}$ are uniformly bounded in absolute value over $i,t,N,T$.
\end{itemize}
\end{assumption}
\begin{assumption}[Concave loglikelihood]\label{assu:convexity}
The function $\ell_{it}(z)$ is strictly concave in $z$ a.s.
Furthermore, there exist constants $\underline{b}>0$, $\overline{b}>0$, and $\delta>0$ such that, for every $\theta \in \Theta_0$,
\[
\underline{b} \leq \min_{i,t}
\{-\partial_{z^2} \ell_{it}(z)\}
\leq \max_{i,t}
\{-\partial_{z^2} \ell_{it}(z)\} \leq \overline{b} \quad \text{a.s.},
\]
for all $z \in \mathcal{B}(\delta, Z_{it}(\theta))$ and all $N,T$.
\end{assumption}
In Assumption \ref{assu:smooth}, the dominating function \(L_{it}(Z_{it},\theta) \) may depend on \( \theta\). Unlike \citet{FVWeidner2016}, who only require boundedness at
the true value \( \theta_0 \), we impose the condition uniformly over $ \Theta_0$. The boundedness condition on $X_{it}$ is also imposed in \citet{ChenFVWeidner2021} and \citet{GaoLiuPengYan2023}.\footnote{In \citet{GaoLiuPengYan2023}, it is required that
\( Z_{it0} = X_{it}^{\prime}\beta_0 + \alpha_{i0}^{\prime}\gamma_{t0} \)
be uniformly bounded over \( i,t \) wpa1.}
The uniform boundedness condition is necessary for establishing uniform asymptotic expansions and uniform asymptotic normality.
Assumption \ref{assu:convexity} implies point identification for any \( \theta \in \Theta_0 \) that might be the true value, as shown in \citet{FVWeidner2016}. Additionally, we require that the Hessian is uniformly bounded away from zero and infinity, which is necessary for uniform consistency and a uniform expansion of the MLE. Importantly, this assumption is satisfied for many popular nonlinear models, e.g, logit and probit models.
Given Assumptions \ref{assu:parameterspace}--\ref{assu:convexity}, we show in Appendix~\ref{Appendix B} that the following expansion holds uniformly in $\theta\in\Theta_0$:
\begin{equation}\label{eq:uniform_expansion}
\sqrt{NT}(\hat{\beta} - \beta)
= \overline{W}_{NT}^{-1}(\theta)\bigl(U_{NT}(\theta) + \overline{B}_{NT}(\theta)\bigr)
+ r_{NT}(\theta).
\end{equation}
Here, the leading term $U_{NT}(\theta)$ has mean zero and variance
$\overline{W}_{NT}(\theta)$.
As $N,T\to\infty$, and provided that the required limits exist, the term $U_{NT}(\theta)$
determines the asymptotic variance, whereas the term $\overline{B}_{NT}(\theta)$ governs the asymptotic bias. The remainder term $r_{NT}(\theta)$ is asymptotically negligible uniformly in $\theta \in \Theta_0$. For brevity, we omit the explicit expressions for $U_{NT}(\theta)$, $\overline{B}_{NT}(\theta)$, and $\overline{W}_{NT}(\theta)$, which are involved; detailed formulas are provided in Appendix~\ref{Appendix B}.
To establish uniform asymptotic normality, we require that the first term in the above expansion converges in distribution to a normal distribution and that the bias term converges to a constant uniformly. However, these quantities depend on the realization of the unobserved effects $\alpha$ and $\gamma$. For instance, if the time effect $\gamma_t$ is a nonstationary time series, the bias term $\overline{B}_{NT}(\theta)$ may not converge. Moreover, to ensure uniform asymptotic normality, the approximation error must decay uniformly for every $\theta \in \Theta_0$. Formally, we state the following assumption:
\begin{assumption}
\label{assu:limit}
For every $\theta \in \Theta_0$, the limits
$\overline{W}_\infty(\theta)=\lim_{N,T \to \infty} \overline{W}_{NT}(\theta)$ and
$\overline{B}_\infty(\theta)=\lim_{N,T \to \infty} \overline{B}_{NT}(\theta)$ exist, and
$\overline{W}_\infty(\theta)$ is positive definite. Furthermore, the converge of
$\overline{W}_{NT}(\theta)$ and
$\overline{B}_{NT}(\theta)$ to their limits is uniform over $\Theta_0$.
\end{assumption}
\citet{FVWeidner2016} assume that
$\overline{W}_{NT}(\theta_{0})$ and $\overline{B}_{NT}(\theta_{0})$
converge in probability to nonrandom limits as $N,T\to\infty$, which implies the
asymptotic normality of $\hat{\beta}$ unconditionally w.r.t.\ the unobserved effects.
Our assumption is stronger in that it requires the convergence to
hold for any realization of unobserved effects (almost surely), and uniformly
over $\Theta_{0}$.
However, it also allows the limits to depend on $\theta$. A similar assumption is implicitly imposed in \citet{HigginsJochmans2024}.
Given Assumptions \ref{assu:parameterspace}--\ref{assu:limit},
we show that bootstrap consistency holds for any realization of the unobserved effects.
This result provides the theoretical foundation for the parametric bootstrap.
The corresponding result for one-way fixed-effect models is given in
\citet{HigginsJochmans2024}.
\begin{theorem}[Uniform asymptotic normality]
\label{theorem3.1} Suppose Assumptions \ref{assu:parameterspace}--\ref{assu:limit} hold. Then, for every $a\in\mathbb{R}^{d_x}$,
\[
\sup_{\theta\in\Theta_0}|\mathbb{P}_{\theta}(\sqrt{NT}(\hat{\beta}-\beta)\leq a)-G_{\theta}(a)| \to 0,
\]
where $G_{\theta}$ is the distribution function corresponding to
$\mathcal{N}\left(\overline{B}_{\infty}(\theta),\overline{W}_\infty(\theta)^{-1}\right)$.
\end{theorem}
The next result establishes that the parametric bootstrap approximates the distribution of the MLE, ensuring the asymptotic validity of confidence intervals constructed using the quantiles of the bootstrap distribution.
\begin{theorem}[Bootstrap consistency for the MLE]
\label{theorem3.2}
Suppose Assumptions \ref{assu:parameterspace}--\ref{assu:limit} hold. Then
\[
\mathbb{P}_{\theta_0}\left[\sup_{a\in\mathbb{R}^{d_x}}
\left|\mathbb{P}_{\hat{\theta}}(\sqrt{NT}(\hat{\beta}^{*}-\hat{\beta})\leq a)-\mathbb{P}_{\theta_{0}}(\sqrt{NT}(\hat{\beta}-\beta_{0})\leq a)\right|>\epsilon\right] \to 0
\]
for all $\epsilon>0$.
\end{theorem}
The following corollary establishes the asymptotic validity of the transformed bootstrap confidence interval. It follows from the bootstrap delta method \citep{van1998asymptotic} that the bootstrap consistency result extends to smooth transformations of the estimator. A more general result can be found in the proof of Theorem~2 in \citet{higgins2025inference}.
\begin{corollary}[Bootstrap consistency for transformed MLE]
\label{corollary3.3}
Suppose Assumptions \ref{assu:parameterspace}--\ref{assu:limit} hold. Let $c \in \mathbb{R}^{d_x}$ be a fixed vector and define $\phi = \varphi(c^\prime \beta)$, where $\varphi:\mathbb{R}\to\mathbb{R}$ is continuously differentiable at $c^\prime \beta_0$ with $\varphi'(c^\prime \beta_0)\neq 0$. Then
\[
\mathbb{P}_{\theta_0}\left[\sup_{a\in\mathbb{R}}
\left|
\mathbb{P}_{\hat{\theta}}\bigl(\sqrt{NT}(\varphi(c^\prime \hat{\beta}^*)-\varphi(c^\prime \hat{\beta}))\le a\bigr)
-
\mathbb{P}_{\theta_0}\bigl(\sqrt{NT}(\varphi(c^\prime \hat{\beta})-\varphi(c^\prime \beta_0))\le a\bigr)
\right|>\epsilon\right]\to 0
\]
for all $\epsilon>0$.
\end{corollary}
In our simulations, we further restrict $\varphi(\cdot)$ to be strictly increasing and continuously differentiable.
We make the following assumptions to establish the validity of the parametric bootstrap for the APEs defined in \eqref{APE}.
\begin{assumption}[]\label{assu:APEs interactive}
The unobserved effects enter through a factor structure:
\[
\mu_{it}(\beta, \alpha_i, \gamma_t)=\mu^f_{it}(\beta, \pi_{it}), \qquad \pi_{it}=\alpha_i^\prime\gamma_t,
\]
for some function $\mu^f_{it}$.
\end{assumption}
\begin{assumption}[]\label{assu:APEs smooth}
There exists $\delta>0$ such that for every $i,t, N, T$ and every $\theta\in\Theta_0$, (i) \( \mu^f_{it}(\cdot, \cdot) \) is four times continuously differentiable over \( \mathcal{B}(\delta;\beta, \pi_{it}) \) a.s.; (ii) the partial derivatives of \( \mu^f_{it}(\cdot, \cdot) \), up to the fourth order, are uniformly bounded in absolute value over $\mathcal{B}(\delta;\beta, \pi_{it}) $.
\begin{comment}
{\color{red}{OLD: a.s.\ by some function \( M_{it}(D_{it},X_{it},\Theta_0) > 0 \) satisfying
\[
\sup_{\Theta_0 \in \Theta_0} \max_{(i,j)\in\mathcal{D}} \mathbb{E}_{\Theta_0}\left[
M_{it}(D_{it},X_{it},\Theta_0)^{8+\nu'}\right]<c'
\]
for some constants $\nu'>0$ and $c'>0$.}
}
\end{comment}
\end{assumption}
Given Assumptions \ref{assu:parameterspace}--\ref{assu:APEs smooth}, the following uniform asymptotic expansion holds for APEs:\footnote{See Theorem \ref{Thm:C2} for expressions of the terms in the expansion.}
\[
\sqrt{NT}\left(\varDelta(\hat{\theta}) - \varDelta(\theta)\right)
= U_{NT}^{\varDelta}(\theta) + \overline{B}_{NT}^{\varDelta}(\theta) + r_{NT}^{\varDelta}(\theta).
\]
Assuming that limits exist, as $N,T\to\infty$, the leading term $U_{NT}^{\varDelta}(\theta)$ contributes to the asymptotic variance, the second term $\overline{B}_{NT}^{\varDelta}(\theta)$ captures the asymptotic bias, and the remainder term $r_{NT}^{\varDelta}(\theta)$ is asymptotically negligible uniformly over $\Theta_0$.
To establish uniform asymptotic normality of the MLE of APEs and bootstrap consistency, we also require a uniform convergence condition. Let the conditional variance of the leading term be denoted by
$\overline{W}_{NT}^{\varDelta}(\theta)$.
\begin{assumption}
\label{assu:limit:ape}
For every $\theta \in \Theta_0$, the limits
$\overline{W}^\varDelta_{\infty}(\theta) =\lim_{N,T \to \infty} \overline{W}^\varDelta_{NT}(\theta)$ and
$\overline{B}^\varDelta_{\infty}(\theta) =\lim_{N,T \to \infty} \overline{B}^\varDelta_{NT}(\theta)$
exist, and $\overline{W}^\varDelta_{\infty}(\theta)>0$. Furthermore, the convergence of
$\overline{W}^\varDelta_{NT}(\theta)$ and $\overline{B}^\varDelta_{NT}(\theta)$ to their limits is uniform over $\Theta_0$.
\end{assumption}
\begin{comment}
{\color{red} The next assumption is problematic.
\begin{assumption}[]\label{assu: APEs identification}
There exist positive constants $\underline{b}'$ and $\overline{b}'$ such that
for every $(i,j)\in\mathcal{D}$ and every
\( \Theta_0 \in \Theta_0 \), there exists \(\mathcal{B}_{\varepsilon}(\beta_1,\pi_{ij1}) \) such that
for all \( (\beta,\pi)\in\mathcal{B}_{\varepsilon}(\beta_1,\pi_{ij1}) \),
\[
\underline{b}' \leq
\inf_{\Theta_0 \in \Theta_0} \left[\mathbb{E}_{\Theta_0}(\mu_{it}^2)-\mathbb{E}_{\Theta_0}(\mu_{it})^2\right]
\leq \sup_{\Theta_0 \in \Theta_0} \left[\mathbb{E}_{\Theta_0}(\mu_{it}^2)-\mathbb{E}_{\Theta_0}(\mu_{it})^2\right]
\leq \overline{b}' \quad \text{a.s.},
\]
where $\mu_{it}=\mu(X_{it},\beta,\pi)$.
\end{assumption}
}
\end{comment}
We can now show bootstrap consistency for APEs.
\begin{theorem}[Uniform asymptotic normality for APEs]
\label{theorem3.3} Suppose Assumptions \ref{assu:parameterspace}--\ref{assu:limit:ape} hold. Then, for every $a\in\mathbb{R}$,
\[
\sup_{\theta\in\Theta_0}|\mathbb{P}_{\theta}(\sqrt{NT}(\varDelta(\hat{\theta}) - \varDelta(\theta))\leq a)-G^\varDelta_{\theta}(a)| \to 0,
\]
where $G^\varDelta_{\theta}$ is the distribution function corresponding to
$\mathcal{N}(\overline{B}^\varDelta_{\infty}(\theta),\overline{W}^\varDelta_{\infty}(\theta)^{-1})$.
\end{theorem}
\begin{theorem}[Bootstrap consistency for APEs]\label{theorem3.4}
Suppose Assumptions \ref{assu:parameterspace}--\ref{assu:limit:ape} hold. Then
\[
\mathbb{P}_{\theta_0}\left[\sup_{a\in \mathbb{R}}\left|\mathbb{P}_{\hat{\theta}}(\sqrt{NT}(\varDelta(\hat{\theta}^{*})-\varDelta(\hat{\theta}))\leq a)-\mathbb{P}_{\theta_0}
(\sqrt{NT}(\varDelta(\hat{\theta})-\varDelta(\theta_{0}))\leq a)
\right|>\epsilon\right] \to 0
\]
for all $\epsilon>0$.
\end{theorem}
We focus on APEs of the form~\eqref{APE},
which means that we consider the average partial effect evaluated at the realization of the unobserved effects, rather than averaged over the marginal distributions of the unobserved effects (i.e., the unconditional APEs).
There are two reasons for this choice.
First, studying the properties of unconditional APEs requires additional assumptions on the distributions of $\{\alpha_i\}_{i=1,\ldots,N}$ and $\{\gamma_t\}_{t=1,\ldots,T}$,
for example, that $\{\alpha_i\}$ are i.i.d.\ across $i$ and $\{\gamma_t\}$ form a stationary time series.
Second, as discussed in \citet{FVWeidner2016}, although the plug-in estimator $\hat{\Delta}$ is still consistent for the unconditional APEs, its convergence rate is slower than $\sqrt{NT}$.
To see this, let $\mathbb{E}(\Delta)$ denote the unconditional APE, where the expectation is taken over the marginal distributions of the unobserved effects.
The difference between the plug-in estimator and the unconditional APE can be decomposed as
\begin{align} \label{eq:decomposition}
{\Delta(\hat\theta)}-\mathbb{E}(\Delta)
= \underbrace{{\Delta(\hat\theta)}-\Delta(\theta)}_{\text{parameter estimation error}}
+ \underbrace{\Delta(\theta)-\mathbb{E}(\Delta)}_{\text{sample mean error}},
\end{align}
where the first term arises from parameter estimation error and the second term reflects the difference between the sample mean and the population mean.
As discussed in \citet{FVWeidner2016}, under mild weak-dependence assumptions on $\{\alpha_i\}_{i=1,\ldots,N}$ and $\{\gamma_t\}_{t=1,\ldots,T}$, the variance of the second term dominates that of the first term. Hence, in this case, there is no asymptotic bias, but the convergence rate becomes slower than $\sqrt{NT}$.\footnote{\citet{FVWeidner2016} show that if $\{\alpha_i\}_{i=1,\ldots,N}$ and $\{\gamma_t\}_{t=1,\ldots,T}$ are both independent sequences and $\alpha_i$ and $\gamma_t$ are independent for all $i,t$, then the convergence rate of the second term in \eqref{eq:decomposition} is $\sqrt{NT/(N+T-1)}$, see Remark 4 in \citet{FVWeidner2016} for further details.}
Our bootstrap procedure focuses on addressing the asymptotic bias in the first term, although the resulting bias-corrected estimator may still provide improved finite-sample performance for the estimation of unconditional APEs.
For further details, see Remark~4 in \citet{FVWeidner2016}.
\section{Implementation}
Since the loglikelihood function in \eqref{eq:loglikelihood} is not strictly concave due to the presence of interactive fixed effects, it may be difficult to compute the MLE, i.e., the global maximizer of the loglikelihood. \citet{ChenFVWeidner2021} suggest using multiple initial values to increase the probability of reaching the global maximizer. However, this strategy does not guarantee that the global maximizer is reached, and it may substantially increase the computational burden, particularly when combined with bootstrapping.
To address the computational challenge, \citet{zeleneev2026tractable} extend the nuclear-norm penalized estimator proposed by \citet{moon2026nuclear} to nonlinear factor models. Their approach replaces the non-convex rank constraint arising from the interactive fixed effects with a nuclear-norm penalty, thereby relaxing the original optimization problem to a strictly convex, computationally tractable problem. \citet{zeleneev2026tractable} establish the consistency of the resulting penalized estimator and further show that using this estimator as the initial value in a gradient descent algorithm leads to convergence to the global maximizer of the loglikelihood, making it asymptotically equivalent to the MLE.
Given the computational efficiency and asymptotic equivalence of the two-step estimator proposed by \citet{zeleneev2026tractable}, we adopt this estimator as a computationally convenient proxy for the MLE in the subsequent analysis. Formally, we first solve the following nuclear-norm penalized optimization problem:
\begin{equation}\label{nuc problem}
(\hat{\beta}_{\text{nuc}}, \hat{\Sigma}_{\text{nuc}})
= \underset{(\beta, \Sigma) \in \mathcal{B} \times \mathbb{R}^{N \times T}}{\arg\max}
\left\{
\frac{1}{NT}\mathcal{L}(\beta, \Sigma)
+ \frac{\varphi}{\sqrt{NT}} \left\lVert \Sigma \right\rVert_{\text{nuc}}
\right\}
,
\end{equation}
where $\mathcal{L}(\beta,\Sigma)
= \sum_{i=1}^{N}\sum_{t=1}^{T}
\!\log f\big(Y_{it}\mid X_{it}^{\prime}\beta+\Sigma_{it}\big)$
is the loglikelihood function without low-rank constraint on $\Sigma$, $\Sigma$ is an $N \times T$ matrix collecting the unobserved effects for each $(i,t)$,
$\lVert\cdot\rVert_{\text{nuc}}$ is the nuclear norm,
and $\varphi$ is a tuning parameter.
Section 4.3 of \citet{zeleneev2026tractable} proposes a data-dependent procedure for selecting the tuning parameter $\varphi$. Specifically, they first compute a two-way fixed-effect estimator, obtaining $\tilde{\beta}$ and the scalar effects $\tilde{\alpha}_1,\ldots,\tilde{\alpha}_N$, and $\tilde{\gamma}_1,\ldots,\tilde{\gamma}_T$, and form an initial guess
\begin{equation}\label{choicevarphi}
\tilde{\varphi}
= 1.05 \left\| \partial_{\Sigma} \mathcal{L}\bigl(\tilde{\beta}, \tilde{\Sigma}\bigr) \right\|_{\mathrm{op}},
\end{equation}
where $\tilde{\Sigma}_{it}=\widetilde{\alpha_i}+\widetilde{\gamma_t}$ and $\partial_{\Sigma} \mathcal{L}(\beta,\Sigma)\in\mathbb{R}^{N\times T}$ is the matrix of partial derivatives of the objective function with respect to $\Sigma$, and $\|\cdot\|_{\mathrm{op}}$ denotes the operator norm (i.e., the largest singular value).
Solving \eqref{nuc problem} with \(\varphi=\tilde{\varphi}\) yields \(\tilde{\beta}_{\mathrm{nuc}}\) and an initial estimate \(\tilde{\Sigma}_{\mathrm{nuc}}^{\mathrm{init}}\). The singular value decomposition of \(\tilde{\Sigma}_{\mathrm{nuc}}^{\mathrm{init}}\) is then computed, and the first \(d_f\) left and right singular vectors are retained to construct initial estimated factor loadings \(\tilde{\alpha}_{\mathrm{nuc}}^{\mathrm{init}}\in\mathbb{R}^{d_f\times N}\) and factors \(\tilde{\gamma}_{\mathrm{nuc}}^{\mathrm{init}}\in\mathbb{R}^{d_f\times T}\). The tuning parameter is subsequently updated as
\begin{equation}
\label{eq:update}
\hat{\varphi}
=
1.05
\left\|
\partial_{\Sigma}\mathcal{L}\bigl(\tilde{\beta}_{\mathrm{nuc}},\tilde{\Sigma}_{\mathrm{nuc}}\bigr)
\right\|_{\mathrm{op}},
\end{equation}
where the entries of \(\tilde{\Sigma}_{\mathrm{nuc}}\) are given by
$(\tilde{\Sigma}_{\mathrm{nuc}})_{it}
=
(\tilde{\alpha}_{\mathrm{nuc},i}^{\mathrm{init}})^{\prime}(\tilde{\gamma}_{\mathrm{nuc},t}^{\mathrm{init}}).$
In an unreported simulation, we found that the choice in \eqref{choicevarphi} tends to induce excessive shrinkage in the singular values of the initial estimate \(\tilde{\Sigma}_{\mathrm{nuc}}^{\mathrm{init}}\). To mitigate this issue, we instead scale down the initial choice and set
\[
\tilde{\varphi}
= 0.5 \left\| \partial_{\Sigma} \mathcal{L}\bigl(\tilde{\beta}, \tilde{\Sigma}\bigr) \right\|_{\mathrm{op}}
\]
as the initial tuning parameter. We then update the tuning parameter as in \citet{zeleneev2026tractable}, that is, using \eqref{eq:update}, thereby maintaining consistency with their theoretical framework. The final choice of the tuning parameter is then $\phi=\hat{\varphi}$, to be used in \eqref{nuc problem} to obtain the nuclear-norm
penalized estimator
$(\hat{\beta}_{\text{nuc}}, \hat{\Sigma}_{\text{nuc}})$.
Additional computational details can be found in Sections 4.1 and 4.3 of \citet{zeleneev2026tractable}.
Next, we apply a gradient descent algorithm initialized at the nuclear-norm penalized estimator. Specifically, we first recover the factor structure by computing the singular value decomposition of $\hat{\Sigma}_{\text{nuc}}$, which yields the estimators $\hat{\alpha}_{\text{nuc}}$ and $\hat{\gamma}_{\text{nuc}}$. The gradient descent iterations are then carried out starting from $(\hat{\beta}_{\text{nuc}}, \hat{\alpha}_{\text{nuc}}, \hat{\gamma}_{\text{nuc}})$.
The iterative updates are given by
\begin{align}
\beta^{(k+1)} &= \beta^{(k)} - S_\beta \, \partial_\beta \mathcal{L}(\beta^{(k)}, \alpha^{(k)}, \gamma^{(k)}), \notag \\
\alpha^{(k+1)} &= \alpha^{(k)} - S_\alpha \, \partial_\alpha \mathcal{L}(\beta^{(k)}, \alpha^{(k)}, \gamma^{(k)}), \label{eq:gradient-update} \\
\gamma^{(k+1)} &= \gamma^{(k)} - S_\gamma \, \partial_\gamma \mathcal{L}(\beta^{(k)}, \alpha^{(k)}, \gamma^{(k)}), \notag
\end{align}
where $S_\beta$, $S_\alpha$, and $S_\gamma$ denote the step sizes for the corresponding parameter blocks;
see \citet{zeleneev2026tractable} for details about the step sizes. We repeat this algorithm until convergence of all three parameters. Since the resulting two-step estimator is asymptotically equivalent to the MLE, we denote it simply by $\hat\theta$.\footnote{Theoretical guarantees for the convergence of this algorithm are provided in Theorem 5 of \citet{zeleneev2026tractable}. We compute $\hat{\theta}$ using the R package \texttt{NNRPanel} developed by \citet{zeleneev2026tractable}. For more details, see \url{https://github.com/wszhang-econ/NNRPanel}.}
In our theoretical and simulation results, we treat the number of factors, $d_f$, as known and fixed. In applications, however, the number of factors can be estimated using the eigenvalue-ratio test \citep{ahn2013eigenvalue}; see also
\citet{ChenFVWeidner2021}, \citet{GaoLiuPengYan2023}, and \citet{zeleneev2026tractable}.
The following algorithm summarizes our bootstrap inference procedure for $\beta_0$.
\begin{algorithm}[H]
\caption{Bootstrap inference for $\beta_0$}
\label{alg:bootstrap}
\begin{algorithmic}[1]
\State Compute the two-step MLE $\hat{\theta}$ based on the original sample.
\State Substitute the estimated $\hat{\theta}$ into model~\eqref{eq:loglikelihood} to generate $B$ parametric bootstrap samples.
\For{$b = 1, \ldots, B$}
\State Compute the two-step MLE $\hat{\theta}^{(b)}$ based on the $b$-th parametric bootstrap sample.
\EndFor
\State Compute the bias-corrected estimator
\vspace{-0.5cm}
\[
\hat{\beta}_{\text{BC-Mean}} = 2\hat{\beta} - \operatorname{mean}\{\hat{\beta}^{(1)},\ldots,\hat{\beta}^{(B)}\}, \quad
\hat{\beta}_{\text{BC-Median}} = 2\hat{\beta} - \operatorname{median}\{\hat{\beta}^{(1)},\ldots,\hat{\beta}^{(B)}\}.
\]
\State \textbf{(Standard bootstrap CI)} For a scalar $\beta_0$, given a confidence level $1-\alpha$, compute
$Q^*_{\hat\beta^*}(\alpha/2)$
and
$Q^*_{\hat\beta^*}(1-\alpha/2)$,
the quantiles of $\{\hat{\beta}^{(1)},\ldots,\hat{\beta}^{(B)}\}$. The confidence interval is
\[
[ 2\hat{\beta} - Q^*_{\hat\beta^*}(1-\alpha/2),\;
2\hat{\beta} - Q^*_{\hat\beta^*}(\alpha/2)].
\]
\State \textbf{(Transformed bootstrap CI)} Let $\varphi(\cdot)$ be a strictly increasing and continuously differentiable transformation. Compute the transformed bootstrap estimates $\{\varphi(\hat{\beta}^{(1)}), \ldots, \varphi(\hat{\beta}^{(B)})\}$ and the corresponding empirical quantile function $Q^*_{\varphi(\hat{\beta}^*)}$. The transformed bootstrap confidence interval is
\[
\left[
\varphi^{-1}\!\left(2\varphi(\hat{\beta}) - Q^*_{\varphi(\hat{\beta}^*)}(1-\alpha/2)\right), \;
\varphi^{-1}\!\left(2\varphi(\hat{\beta}) - Q^*_{\varphi(\hat{\beta}^*)}(\alpha/2)\right)
\right].
\]
\end{algorithmic}
\end{algorithm}
\section{Simulations}
Following \citet{ChenFVWeidner2021}, the two bias-corrected methods can also be applied to the two-step MLE estimator in \citet{zeleneev2026tractable}. One is a plug-in analytical correction based on the explicit form of the first-order bias, and another one is known as the split-panel jackknife method based on the idea of \citet{DhaeneJochmans2015}. We compare the finite-sample performance of our bootstrap bias correction with the two methods.
We focus on the following static model with interactive fixed effects:
\begin{equation} \label{MC}
D_{it} = \mathbf{1}\!\left\{ X_{it}^{\prime}\beta_0 + \alpha_{i0}^{\prime}\gamma_{t0} - \varepsilon_{it} > 0 \right\},
\quad i = 1,\dots,N,\; t = 1,\dots,T,
\end{equation}
where $\mathbf{1}\{\cdot\}$ is an indicator function and $\varepsilon_{it}$ follows a standard logistic distribution or a standard normal distribution, corresponding to the logit and the probit model. We let $N=30$, and consider $T \in \{20,30,40\}$. There is only $K=1$ covariate, and the number of factors $d_f = 2$, and the true coefficient is $\beta_0 = 0.5$. For all estimators, we fix $d_f = 2$. The following scenarios are included:
\paragraph{Scenario 1 (Covariates independent of factors):} For each $i$ and $t$, draw factor loadings and factors $\alpha_i,\gamma_t \sim \mathcal{N}(0,I_{d_f})$ and independent regressor factors $\alpha_i^{(x)},\gamma_t^{(x)} \sim \mathcal{N}(0,I_{d_f})$ where $I_{d_f}$ is the identity matrix of dimension $d_f$, and generate covariates as $X_{itk}=\alpha_i^{(x)\prime}\gamma_t^{(x)}+\nu_{itk}$ with $\nu_{itk}\sim\mathcal{N}(0,1)$.
\paragraph{Scenario 2 (Covariates correlated with factors):} For each $i$ and $t$, draw factor loadings and factors $\alpha_i,\gamma_t \sim \mathcal{N}(0,I_{d_f})$, and generate covariates as $X_{itk}=0.3 \times \alpha_i^{\prime}\gamma_t+\nu_{itk}$ with $\nu_{itk}\sim\mathcal{N}(0,1)$.
\paragraph{Scenario 3 (Outliers in factor loadings):} For each $i$ and $t$, draw factor loadings and factors $\alpha_i,\gamma_t \sim \mathcal{N}(0,I_{d_f})$, fix a random $10\%$ of $\alpha_i = 3$ as outliers, draw independent regressor factors $\alpha_i^{(x)},\gamma_t^{(x)} \sim \mathcal{N}(0,I_{d_f})$, and generate covariates as $X_{itk}=\alpha_i^{(x)\prime}\gamma_t^{(x)}+\nu_{itk}$ with $\nu_{itk}\sim\mathcal{N}(0,1)$.
\paragraph{Scenario 4 (Sparse spikes in factors):} For each $i$ and $t$, draw factor loadings and factors $\alpha_i,\gamma_t \sim \mathcal{N}(0,I_{d_f})$, multiply a random 5\% of $\gamma_t$ by 10 to create sparse spikes, draw independent regressor factors $\alpha_i^{(x)},\gamma_t^{(x)} \sim \mathcal{N}(0,I_{d_f})$, and generate covariates as $X_{itk}=\alpha_i^{(x)\prime}\gamma_t^{(x)}+\nu_{itk}$ with $\nu_{itk}\sim\mathcal{N}(0,1)$.
We repeat the simulation $1000$ times and use $399$ bootstrap samples in each for the bootstrap method. Tables \ref{sim res logit} and \ref{sim res probit} report the relative bias (defined as the bias divided by the true parameter value), standard deviation, $95\%$ coverage rate, and the rejection probabilities under the null hypotheses: $H_0: \beta_0 = 0.7$ for all scenarios. The simulation results are for the logit and probit models, respectively. We compare the uncorrected maximum likelihood estimator, the split panel jackknife, the analytical bias correction, and the bootstrap method, where Boot-Mean uses the mean of bootstrap estimates as the corrected term, while Boot-Median uses the median.
In most scenarios, the bootstrap-based methods deliver the most effective bias correction, exhibiting uniformly smaller biases than competing approaches. In Scenario 2 of both the logit and probit models, only the bootstrap methods achieve coverage probabilities close to the nominal $95\%$ level, whereas both the split-panel jackknife and the analytical bias correction perform less well in terms of bias reduction and coverage accuracy. As the sample size increases, the gap between the analytical bias correction and the bootstrap methods gradually narrows. The split-panel jackknife also improves with sample size but continues to require relatively large samples to achieve comparable bias reduction. Overall, the results suggest that bootstrap-based bias correction is particularly well suited to relatively small panel data settings, while remaining competitive with analytical bias correction in large samples. Moreover, the bootstrap approach performs well across a range of scenarios, especially when covariates are correlated with factors and loadings.
As discussed above, we find that the bootstrap confidence intervals tend to be overly conservative in finite samples, leading to coverage rates exceeding the nominal confidence levels. This pattern is particularly pronounced when the time dimension is small (e.g., \(T=20\)), as illustrated in Tables~\ref{sim res logit}--\ref{sim res probit}. In both models, the coverage rates of the parametric bootstrap confidence intervals are uniformly above $95\%$ across all scenarios, often substantially so, providing clear evidence of systematic over-coverage.
To address this issue, we consider three monotone transformations to reduce skewness in the bootstrap distribution: the log transformation, the Box--Cox transformation \citep{box1964analysis}, and the Yeo--Johnson transformation \citep{yeo2000new}.
The log transformation provides a simple adjustment for strictly positive parameters and is particularly effective in mitigating right skewness by compressing the upper tail of the distribution.
The Box--Cox transformation introduces a tuning parameter that controls the strength and direction of the transformation; in our implementation, we select this parameter in a data-driven manner by minimizing the skewness of the transformed bootstrap estimates.\footnote{
Specifically, we set $\lambda$ equal to
\[
\hat{\lambda}
=
\arg\min_{\lambda \in [-2,2]}
\left|
\frac{
\frac{1}{B}\sum_{b=1}^{B}
\left(\varphi_{\lambda}(\hat{\beta}^{*(b)})-\bar{\varphi}_{\lambda}\right)^3
}{
\left[
\frac{1}{B}\sum_{b=1}^{B}
\left(\varphi_{\lambda}(\hat{\beta}^{*(b)})-\bar{\varphi}_{\lambda}\right)^2
\right]^{3/2}
}
\right|,
\]
where $\hat{\beta}^{*(b)}$, $b=1,\ldots,B$, are the bootstrap estimates, $\varphi_{\lambda}(\cdot)$ is the Box--Cox transformation, and $\bar{\varphi}_{\lambda} = \frac{1}{B}\sum_{b=1}^{B} \varphi_{\lambda}(\hat{\beta}^{*(b)})$.
}
The Yeo--Johnson transformation extends this approach to accommodate both positive and negative values, while retaining similar flexibility in adjusting skewness.
The results in Tables~\ref{tab:logit_coverage_all}--\ref{tab:probit_coverage_all} show that, relative to the parametric bootstrap, the transformation-based parametric bootstrap substantially reduces the confidence interval length while bringing the coverage rate closer to the nominal level, particularly in Scenarios 1, 3, and 4, reflecting more efficient inference. This improvement is driven by their ability to mitigate skewness, resulting in more balanced lower and upper miss rates.
An exception is Scenario 2, where all methods perform relatively poorly. This is mainly due to the limited effectiveness of the bootstrap-based bias correction, which leads to a discrepancy between the bootstrap and true finite-sample distributions of the estimator. As a result, the bootstrap fails to accurately approximate the sampling distribution, reducing the effectiveness of the skewness correction and leading to distortions in coverage.
By contrast, the untransformed parametric bootstrap, which is typically conservative, exhibits relatively better coverage in this case.
The log transformation often yields the shortest intervals, as it compresses right-skewed distributions. Since the MLE is predominantly right-skewed in our simulations, this effect is particularly pronounced. However, its one-sided nature makes it less robust in more complex settings. By contrast, the Yeo-Johnson transformation provides the most reliable overall performance, offering a favorable balance between coverage accuracy and interval length, especially in more asymmetric scenarios.
Overall, these results demonstrate that appropriate monotone transformations can substantially improve the finite-sample performance of bootstrap confidence intervals by reducing conservativeness and yielding tighter, more informative inference.
\begin{table}[H]
\centering
\caption{Simulation results for the logit model, $N=30$}
\label{sim res logit}
\medskip
\scalebox{0.75}{
\begin{tabular}{llrccc|llrccc}
\hline
\multicolumn{6}{c|}{\textit{Scenario 1: covariates independent of factors}} &
\multicolumn{6}{c}{\textit{Scenario 2: covariates correlated with factors}} \\
$T$ & Method & Bias\,\, & SD & Coverage & Rejection &
$T$ & Method & Bias\,\, & SD & Coverage & Rejection \\
\hline
20 & MLE & 0.326 & 0.122 & 0.615 & 0.148 & 20 & MLE & 0.668 & 0.189 & 0.371 & 0.208 \\
& SplitPJ & $-0.467$ & 0.193 & 0.414 & 0.895 & & SplitPJ & $-0.522$ & 0.339 & 0.494 & 0.732 \\
& Analytical & 0.155 & 0.103 & 0.892 & 0.308 & & Analytical & 0.499 & 0.168 & 0.589 & 0.094 \\
& Boot-Mean & $-0.081$ & 0.083 & 0.985 & 0.483 & & Boot-Mean & 0.171 & 0.141 & 0.995 & 0.056 \\
& Boot-Median & $-0.061$ & 0.085 & 0.985 & 0.483 & & Boot-Median & 0.193 & 0.143 & 0.995 & 0.056 \\
\hline
30 & MLE & 0.218 & 0.085 & 0.675 & 0.309 & 30 & MLE & 0.481 & 0.135 & 0.421 & 0.132 \\
& SplitPJ & $-0.228$ & 0.121 & 0.595 & 0.921 & & SplitPJ & $-0.264$ & 0.214 & 0.575 & 0.724 \\
& Analytical & 0.095 & 0.074 & 0.916 & 0.590 & & Analytical & 0.361 & 0.123 & 0.595 & 0.099 \\
& Boot-Mean & $-0.040$ & 0.066 & 0.979 & 0.756 & & Boot-Mean & 0.177 & 0.111 & 0.966 & 0.136 \\
& Boot-Median & $-0.031$ & 0.066 & 0.979 & 0.756 & & Boot-Median & 0.186 & 0.111 & 0.966 & 0.136 \\
\hline
40 & MLE & 0.177 & 0.066 & 0.679 & 0.472 & 40 & MLE & 0.355 & 0.107 & 0.486 & 0.114 \\
& SplitPJ & $-0.150$ & 0.087 & 0.668 & 0.952 & & SplitPJ & $-0.229$ & 0.164 & 0.594 & 0.822 \\
& Analytical & 0.071 & 0.058 & 0.936 & 0.765 & & Analytical & 0.259 & 0.100 & 0.670 & 0.164 \\
& Boot-Mean & $-0.024$ & 0.053 & 0.977 & 0.882 & & Boot-Mean & 0.128 & 0.092 & 0.952 & 0.265 \\
& Boot-Median & $-0.019$ & 0.054 & 0.977 & 0.882 & & Boot-Median & 0.134 & 0.092 & 0.952 & 0.265 \\
\hline
\multicolumn{6}{c|}{\textit{Scenario 3: outliers in factor loadings}} &
\multicolumn{6}{c}{\textit{Scenario 4: sparse spikes in factors}} \\
$T$ & Method & Bias\,\, & SD & Coverage & Rejection &
$T$ & Method & Bias\,\, & SD & Coverage & Rejection \\
\hline
20 & MLE & 0.314 & 0.119 & 0.661 & 0.129 & 20 & MLE & 0.315 & 0.122 & 0.636 & 0.160 \\
& SplitPJ & $-0.458$ & 0.197 & 0.432 & 0.881 & & SplitPJ & $-0.447$ & 0.192 & 0.423 & 0.895 \\
& Analytical & 0.141 & 0.100 & 0.921 & 0.289 & & Analytical & 0.145 & 0.102 & 0.907 & 0.307 \\
& Boot-Mean & $-0.097$ & 0.080 & 0.979 & 0.473 & & Boot-Mean & $-0.086$ & 0.082 & 0.988 & 0.486 \\
& Boot-Median & $-0.077$ & 0.081 & 0.979 & 0.473 & & Boot-Median & $-0.066$ & 0.083 & 0.988 & 0.486 \\
\hline
30 & MLE & 0.213 & 0.085 & 0.701 & 0.305 & 30 & MLE & 0.215 & 0.085 & 0.697 & 0.303 \\
& SplitPJ & $-0.242$ & 0.116 & 0.567 & 0.921 & & SplitPJ & $-0.241$ & 0.117 & 0.559 & 0.932 \\
& Analytical & 0.085 & 0.074 & 0.926 & 0.576 & & Analytical & 0.090 & 0.075 & 0.922 & 0.579 \\
& Boot-Mean & $-0.049$ & 0.065 & 0.981 & 0.745 & & Boot-Mean & $-0.042$ & 0.066 & 0.981 & 0.749 \\
& Boot-Median & $-0.040$ & 0.066 & 0.981 & 0.745 & & Boot-Median & $-0.033$ & 0.066 & 0.981 & 0.749 \\
\hline
40 & MLE & 0.174 & 0.068 & 0.711 & 0.460 & 40 & MLE & 0.166 & 0.069 & 0.730 & 0.512 \\
& SplitPJ & $-0.152$ & 0.083 & 0.696 & 0.965 & & SplitPJ & $-0.162$ & 0.087 & 0.644 & 0.964 \\
& Analytical & 0.063 & 0.060 & 0.929 & 0.788 & & Analytical & 0.058 & 0.061 & 0.937 & 0.798 \\
& Boot-Mean & $-0.033$ & 0.054 & 0.974 & 0.892 & & Boot-Mean & $-0.034$ & 0.055 & 0.979 & 0.889 \\
& Boot-Median & $-0.027$ & 0.055 & 0.974 & 0.892 & & Boot-Median & $-0.028$ & 0.056 & 0.979 & 0.889 \\
\hline
\end{tabular}
}
\medskip
\begin{minipage}{0.95\textwidth}
\footnotesize
\textit{Notes:} model given in \eqref{MC} with one covariate and $\beta_0=0.5$; $1000$ Monte Carlo replications; $399$ bootstrap replications; coverage rates of nominal $95\%$ confidence intervals; rejection probabilities of $H_0:\beta_0=0.7$ at nominal $5\%$ level.
\end{minipage}
\end{table}
\begin{table}[H]
\centering
\caption{Simulation results for the probit model, $N=30$}
\label{sim res probit}
\medskip
\scalebox{0.75}{
\begin{tabular}{llrccc|llrccc}
\hline
\multicolumn{6}{c|}{\textit{Scenario 1: covariates independent of factors}} &
\multicolumn{6}{c}{\textit{Scenario 2: covariates correlated with factors}} \\
$T$ & Method & Bias\,\, & SD & Coverage & Rejection &
$T$ & Method & Bias\,\, & SD & Coverage & Rejection \\
\hline
20 & MLE & 0.374 & 0.108 & 0.409 & 0.139 & 20 & MLE & 0.472 & 0.144 & 0.457 & 0.129 \\
& SplitPJ & $-0.617$ & 0.242 & 0.249 & 0.959 & & SplitPJ & $-0.651$ & 0.259 & 0.341 & 0.911 \\
& Analytical & 0.152 & 0.082 & 0.911 & 0.384 & & Analytical & 0.283 & 0.124 & 0.744 & 0.130 \\
& Boot-Mean & $-0.184$ & 0.058 & 0.981 & 0.719 & & Boot-Mean & $-0.047$ & 0.090 & 0.995 & 0.308 \\
& Boot-Median & $-0.148$ & 0.058 & 0.981 & 0.719 & & Boot-Median & $-0.016$ & 0.093 & 0.995 & 0.308 \\
\hline
30 & MLE & 0.264 & 0.070 & 0.418 & 0.294 & 30 & MLE & 0.315 & 0.092 & 0.501 & 0.144 \\
& SplitPJ & $-0.264$ & 0.100 & 0.436 & 0.989 & & SplitPJ & $-0.299$ & 0.132 & 0.515 & 0.926 \\
& Analytical & 0.100 & 0.057 & 0.902 & 0.704 & & Analytical & 0.173 & 0.082 & 0.818 & 0.347 \\
& Boot-Mean & $-0.078$ & 0.045 & 0.982 & 0.909 & & Boot-Mean & 0.003 & 0.069 & 0.995 & 0.553 \\
& Boot-Median & $-0.064$ & 0.046 & 0.982 & 0.909 & & Boot-Median & 0.015 & 0.070 & 0.995 & 0.553 \\
\hline
40 & MLE & 0.207 & 0.053 & 0.439 & 0.549 & 40 & MLE & 0.253 & 0.074 & 0.503 & 0.258 \\
& SplitPJ & $-0.169$ & 0.071 & 0.553 & 0.998 & & SplitPJ & $-0.194$ & 0.094 & 0.614 & 0.961 \\
& Analytical & 0.069 & 0.045 & 0.927 & 0.910 & & Analytical & 0.128 & 0.066 & 0.850 & 0.559 \\
& Boot-Mean & $-0.053$ & 0.039 & 0.979 & 0.984 & & Boot-Mean & 0.013 & 0.059 & 0.992 & 0.763 \\
& Boot-Median & $-0.046$ & 0.039 & 0.979 & 0.984 & & Boot-Median & 0.020 & 0.060 & 0.992 & 0.763 \\
\hline
\multicolumn{6}{c|}{\textit{Scenario 3: outliers in factor loadings}} &
\multicolumn{6}{c}{\textit{Scenario 4: sparse spikes in factors}} \\
$T$ & Method & Bias\,\, & SD & Coverage & Rejection &
$T$ & Method & Bias\,\, & SD & Coverage & Rejection \\
\hline
20 & MLE & 0.381 & 0.109 & 0.439 & 0.110 & 20 & MLE & 0.354 & 0.103 & 0.479 & 0.124 \\
& SplitPJ & $-0.643$ & 0.259 & 0.255 & 0.961 & & SplitPJ & $-0.598$ & 0.220 & 0.276 & 0.966 \\
& Analytical & 0.149 & 0.085 & 0.917 & 0.376 & & Analytical & 0.134 & 0.081 & 0.938 & 0.419 \\
& Boot-Mean & $-0.207$ & 0.061 & 0.984 & 0.698 & & Boot-Mean & $-0.197$ & 0.057 & 0.983 & 0.735 \\
& Boot-Median & $-0.165$ & 0.060 & 0.984 & 0.698 & & Boot-Median & $-0.159$ & 0.057 & 0.983 & 0.735 \\
\hline
30 & MLE & 0.266 & 0.075 & 0.458 & 0.276 & 30 & MLE & 0.249 & 0.072 & 0.465 & 0.323 \\
& SplitPJ & $-0.272$ & 0.103 & 0.446 & 0.981 & & SplitPJ & $-0.269$ & 0.101 & 0.440 & 0.985 \\
& Analytical & 0.098 & 0.061 & 0.912 & 0.679 & & Analytical & 0.087 & 0.058 & 0.929 & 0.726 \\
& Boot-Mean & $-0.089$ & 0.048 & 0.980 & 0.889 & & Boot-Mean & $-0.088$ & 0.047 & 0.967 & 0.903 \\
& Boot-Median & $-0.074$ & 0.048 & 0.980 & 0.889 & & Boot-Median & $-0.074$ & 0.048 & 0.967 & 0.903 \\
\hline
40 & MLE & 0.215 & 0.057 & 0.431 & 0.480 & 40 & MLE & 0.213 & 0.057 & 0.445 & 0.499 \\
& SplitPJ & $-0.171$ & 0.075 & 0.555 & 0.994 & & SplitPJ & $-0.158$ & 0.067 & 0.611 & 0.997 \\
& Analytical & 0.069 & 0.047 & 0.943 & 0.889 & & Analytical & 0.070 & 0.048 & 0.930 & 0.891 \\
& Boot-Mean & $-0.058$ & 0.040 & 0.975 & 0.978 & & Boot-Mean & $-0.053$ & 0.041 & 0.981 & 0.971 \\
& Boot-Median & $-0.049$ & 0.040 & 0.975 & 0.978 & & Boot-Median & $-0.045$ & 0.041 & 0.981 & 0.971 \\
\hline
\end{tabular}
}
\medskip
\begin{minipage}{0.95\textwidth}
\footnotesize
\textit{Notes:} model given in \eqref{MC} with one covariate and $\beta_0=0.5$; $1000$ Monte Carlo replications; $399$ bootstrap replications; coverage rates of nominal $95\%$ confidence intervals; rejection probabilities of $H_0:\beta_0=0.7$ at nominal $5\%$ level.
\end{minipage}
\end{table}
\begin{table}[H]
\centering
\caption{Comparison of confidence intervals, logit model, $N=30$}
\label{tab:logit_coverage_all}
\medskip
\scalebox{0.75}{
\begin{tabular}{llcccc|llcccc}
\hline
\multicolumn{6}{c|}{\textit{Scenario 1: covariates independent of factors}} &
\multicolumn{6}{c}{\textit{Scenario 2: covariates correlated with factors}} \\
$T$ & Method & Coverage & Length & LMR & UMR &
$T$ & Method & Coverage & Length & LMR & UMR \\
\hline
20 & Boot & 0.985 & 0.558 & 0.000 & 0.015 & 20 & Boot & 0.995 & 0.783 & 0.002 & 0.003 \\
& Log & 0.941 & 0.342 & 0.048 & 0.011 & & Log & 0.636 & 0.508 & 0.364 & 0.000 \\
& Box-Cox & 0.933 & 0.369 & 0.054 & 0.013 & & Box-Cox & 0.796 & 0.578 & 0.204 & 0.000 \\
& Yeo-Johnson & 0.960 & 0.384 & 0.027 & 0.013 & & Yeo-Johnson & 0.826 & 0.594 & 0.172 & 0.002 \\
\hline
30 & Boot & 0.979 & 0.366 & 0.000 & 0.021 & 30 & Boot & 0.966 & 0.523 & 0.034 & 0.000 \\
& Log & 0.929 & 0.257 & 0.054 & 0.017 & & Log & 0.627 & 0.384 & 0.373 & 0.000 \\
& Box-Cox & 0.934 & 0.281 & 0.045 & 0.021 & & Box-Cox & 0.783 & 0.439 & 0.217 & 0.000 \\
& Yeo-Johnson & 0.945 & 0.287 & 0.034 & 0.021 & & Yeo-Johnson & 0.801 & 0.444 & 0.199 & 0.000 \\
\hline
40 & Boot & 0.977 & 0.288 & 0.000 & 0.023 & 40 & Boot & 0.952 & 0.412 & 0.045 & 0.003 \\
& Log & 0.952 & 0.215 & 0.031 & 0.017 & & Log & 0.692 & 0.319 & 0.308 & 0.000 \\
& Box-Cox & 0.949 & 0.235 & 0.031 & 0.020 & & Box-Cox & 0.835 & 0.364 & 0.163 & 0.002 \\
& Yeo-Johnson & 0.952 & 0.238 & 0.028 & 0.020 & & Yeo-Johnson & 0.843 & 0.366 & 0.155 & 0.002 \\
\hline
\multicolumn{6}{c|}{\textit{Scenario 3: outliers in factor loadings}} &
\multicolumn{6}{c}{\textit{Scenario 4: sparse spikes in factors}} \\
$T$ & Method & Coverage & Length & LMR & UMR &
$T$ & Method & Coverage & Length & LMR & UMR \\
\hline
20 & Boot & 0.979 & 0.578 & 0.000 & 0.021 & 20 & Boot & 0.988 & 0.566 & 0.000 & 0.012 \\
& Log & 0.957 & 0.354 & 0.032 & 0.011 & & Log & 0.950 & 0.349 & 0.043 & 0.007 \\
& Box-Cox & 0.942 & 0.382 & 0.042 & 0.016 & & Box-Cox & 0.950 & 0.375 & 0.038 & 0.012 \\
& Yeo-Johnson & 0.956 & 0.399 & 0.028 & 0.016 & & Yeo-Johnson & 0.967 & 0.391 & 0.021 & 0.012 \\
\hline
30 & Boot & 0.981 & 0.378 & 0.000 & 0.019 & 30 & Boot & 0.981 & 0.369 & 0.000 & 0.019 \\
& Log & 0.942 & 0.264 & 0.043 & 0.015 & & Log & 0.938 & 0.260 & 0.049 & 0.013 \\
& Box-Cox & 0.953 & 0.289 & 0.030 & 0.017 & & Box-Cox & 0.943 & 0.283 & 0.039 & 0.018 \\
& Yeo-Johnson & 0.959 & 0.296 & 0.023 & 0.018 & & Yeo-Johnson & 0.948 & 0.289 & 0.033 & 0.019 \\
\hline
40 & Boot & 0.974 & 0.297 & 0.000 & 0.026 & 40 & Boot & 0.979 & 0.292 & 0.001 & 0.020 \\
& Log & 0.938 & 0.220 & 0.038 & 0.024 & & Log & 0.945 & 0.218 & 0.040 & 0.015 \\
& Box-Cox & 0.952 & 0.240 & 0.023 & 0.025 & & Box-Cox & 0.957 & 0.238 & 0.025 & 0.018 \\
& Yeo-Johnson & 0.955 & 0.243 & 0.020 & 0.025 & & Yeo-Johnson & 0.960 & 0.241 & 0.021 & 0.019 \\
\hline
\end{tabular}
}
\medskip
\begin{minipage}{0.95\textwidth}
\footnotesize
\textit{Notes:} model given in \eqref{MC} with one covariate and $\beta_0=0.5$; $1000$ Monte Carlo replications; $399$ bootstrap replications; coverage rates of nominal $95\%$ confidence intervals; Boot: parametric bootstrap; Log/Box-Cox/Yeo-Johnson: transformed parametric bootstrap; LMR: lower miss rate; UMR upper miss rate.
\end{minipage}
\end{table}
\begin{table}[H]
\centering
\caption{Comparison of confidence intervals, probit model, $N=30$}
\label{tab:probit_coverage_all}
\medskip
\scalebox{0.75}{
\begin{tabular}{llcccc|llcccc}
\hline
\multicolumn{6}{c|}{\textit{Scenario 1: covariates independent of factors}} &
\multicolumn{6}{c}{\textit{Scenario 2: covariates correlated with factors}} \\
$T$ & Method & Coverage & Length & LMR & UMR &
$T$ & Method & Coverage & Length & LMR & UMR \\
\hline
20 & Boot & 0.981 & 0.567 & 0.000 & 0.019 & 20 & Boot & 0.995 & 0.704 & 0.004 & 0.001 \\
& Log & 0.929 & 0.349 & 0.060 & 0.011 & & Log & 0.757 & 0.463 & 0.243 & 0.000 \\
& Box-Cox & 0.904 & 0.375 & 0.080 & 0.016 & & Box-Cox & 0.858 & 0.518 & 0.142 & 0.000 \\
& Yeo-Johnson & 0.934 & 0.392 & 0.052 & 0.014 & & Yeo-Johnson & 0.887 & 0.533 & 0.110 & 0.003 \\
\hline
30 & Boot & 0.982 & 0.377 & 0.000 & 0.018 & 30 & Boot & 0.995 & 0.464 & 0.005 & 0.000 \\
& Log & 0.930 & 0.261 & 0.057 & 0.013 & & Log & 0.828 & 0.333 & 0.172 & 0.000 \\
& Box-Cox & 0.912 & 0.283 & 0.071 & 0.017 & & Box-Cox & 0.899 & 0.373 & 0.101 & 0.000 \\
& Yeo-Johnson & 0.935 & 0.293 & 0.050 & 0.015 & & Yeo-Johnson & 0.915 & 0.381 & 0.084 & 0.001 \\
\hline
40 & Boot & 0.979 & 0.297 & 0.000 & 0.021 & 40 & Boot & 0.992 & 0.365 & 0.008 & 0.000 \\
& Log & 0.934 & 0.217 & 0.047 & 0.019 & & Log & 0.857 & 0.279 & 0.143 & 0.000 \\
& Box-Cox & 0.920 & 0.234 & 0.060 & 0.020 & & Box-Cox & 0.913 & 0.311 & 0.087 & 0.000 \\
& Yeo-Johnson & 0.937 & 0.242 & 0.043 & 0.020 & & Yeo-Johnson & 0.926 & 0.317 & 0.073 & 0.001 \\
\hline
\multicolumn{6}{c|}{\textit{Scenario 3: outliers in factor loadings}} &
\multicolumn{6}{c}{\textit{Scenario 4: sparse spikes in factors}} \\
$T$ & Method & Coverage & Length & LMR & UMR &
$T$ & Method & Coverage & Length & LMR & UMR \\
\hline
20 & Boot & 0.984 & 0.573 & 0.000 & 0.016 & 20 & Boot & 0.983 & 0.548 & 0.000 & 0.017 \\
& Log & 0.942 & 0.355 & 0.047 & 0.011 & & Log & 0.944 & 0.338 & 0.044 & 0.012 \\
& Box-Cox & 0.919 & 0.382 & 0.067 & 0.014 & & Box-Cox & 0.918 & 0.364 & 0.070 & 0.012 \\
& Yeo-Johnson & 0.944 & 0.397 & 0.040 & 0.016 & & Yeo-Johnson & 0.944 & 0.378 & 0.043 & 0.013 \\
\hline
30 & Boot & 0.980 & 0.378 & 0.000 & 0.020 & 30 & Boot & 0.967 & 0.357 & 0.001 & 0.032 \\
& Log & 0.933 & 0.265 & 0.053 & 0.014 & & Log & 0.934 & 0.252 & 0.051 & 0.015 \\
& Box-Cox & 0.915 & 0.289 & 0.068 & 0.017 & & Box-Cox & 0.916 & 0.273 & 0.066 & 0.018 \\
& Yeo-Johnson & 0.937 & 0.296 & 0.046 & 0.017 & & Yeo-Johnson & 0.936 & 0.279 & 0.045 & 0.019 \\
\hline
40 & Boot & 0.975 & 0.297 & 0.000 & 0.025 & 40 & Boot & 0.981 & 0.292 & 0.000 & 0.019 \\
& Log & 0.929 & 0.220 & 0.048 & 0.023 & & Log & 0.932 & 0.218 & 0.050 & 0.018 \\
& Box-Cox & 0.913 & 0.240 & 0.061 & 0.026 & & Box-Cox & 0.918 & 0.238 & 0.058 & 0.024 \\
& Yeo-Johnson & 0.932 & 0.243 & 0.043 & 0.025 & & Yeo-Johnson & 0.936 & 0.241 & 0.040 & 0.024 \\
\hline
\end{tabular}
}
\medskip
\begin{minipage}{0.95\textwidth}
\footnotesize
\textit{Notes:} model given in \eqref{MC} with one covariate and $\beta_0=0.5$; $1000$ Monte Carlo replications; $399$ bootstrap replications; coverage rates of nominal $95\%$ confidence intervals; Boot: parametric bootstrap; Log/Box-Cox/Yeo-Johnson: transformed parametric bootstrap; LMR: lower miss rate; UMR upper miss rate.
\end{minipage}
\end{table}
\section{Empirical application: technology spillovers in the presence of latent heterogeneity}
In this section, we revisit the dataset originally from \citet{bloom} and studied by \citet{burdaapplication}, who use a Bayesian panel probit model to estimate a patent equation using firm-level data on research and development ($R\&D$).\footnote{The dataset is publicly available at \url{http://qed.econ.queensu.ca/jae/datasets/burda002/}.} Their study explores a long-standing question in the $R\&D$ literature concerning the relationship between innovation, proxied by firms' patenting activity, and spillover effects arising from strategic interactions across firms. These interactions arise in two main spaces: technology space, defined by the similarity of firms' underlying technologies and associated with beneficial knowledge spillovers, and product market space, defined by the degree of competition in overlapping product markets and associated with business-stealing effects. Furthermore, as emphasized by \citet{burdaapplication}, the innovation and patenting behavior are also likely to be affected by unobserved heterogeneity at the firm and time level.
In \citet{burdaapplication}, they report results from different models considering unobserved heterogeneity, such as the Bayesian panel probit model with two latent effects, the fixed probit model with time dummies, and the random effects probit model with time dummies. However, they do not consider an interactive fixed effect structure in the models. Hence, we employ our estimator to this empirical application.
The original dataset consists of an unbalanced panel of $729$ U.S. firms observed over the period $1981-2001$. To construct a balanced panel, we restrict attention to firms with complete observations throughout the sample period. In addition, firms that do not register any patents over the entire time span are excluded. The resulting sample comprises $328$ firms observed annually from $1981$ to $2001$. The dependent variable is a binary indicator equal to one if a firm files at least one patent in a given year. The key explanatory variables capture two types of spillover effects. Technological spillovers (SpillTech) are constructed using information on the distribution of patents across technological classes. Specifically, technological proximity between firms is measured using the Jaffe distance \citep{jaffe1986}, calculated as the uncentered correlation (cosine similarity) of their patent shares across different technology classes. The corresponding spillover variable is defined as the weighted sum of other firms' $R\&D$ stocks, where weights reflect technological proximity; see details in \citet{burdaapplication}. Product market spillovers (SpillSIC) are constructed analogously using the distribution of firm sales across industries. Product market proximity is measured as the uncentered correlation of firms' sales shares across industry classifications, and the associated spillover variable is defined as the proximity-weighted sum of other firms' $R\&D$ stocks. In addition, we include firm-level controls for the stock of $R\&D$ ($R\&D$ stock) and firm sales (Sales), both constructed from accounting data. All continuous independent variables are expressed in logarithms, lagged by one period, and all regressions include a dummy where lagged $R\&D$ stock is zero.
Prior to estimation, all continuous independent variables are standardized by subtracting their sample means and dividing by their respective standard deviations. We report the results of the following three estimators: the two-step maximum likelihood estimator proposed by \citet{zeleneev2026tractable} (IFE), a plug-in analytical bias correction (IFE-Analytical), and a median bootstrap-based bias correction (IFE-Bootstrap). We use the eigenvalue-ratio test \citep{ahn2013eigenvalue}, slightly modified as in \citet{zeleneev2026tractable}, to determine the number of factors, which is $1$ in this empirical application. In the bootstrap procedure, we draw $399$ bootstrap samples and treat the number of factors as known, which is $1$.
Economic theory provides clear predictions regarding the effects of $R\&D$ spillovers on patenting activity. In the absence of endogenous patenting behavior, product market spillovers, which captures market rivalry, are expected to have no effect on innovation outcomes, implying a coefficient close to zero. In contrast, technological spillovers are expected to have a positive effect, as knowledge diffusion enhances firms' innovative capacity. These predictions serve as a theoretical benchmark for evaluating the empirical results presented in Table \ref{app}.
\begin{table}[H]
\centering
\caption{Parameter Estimates and Standard Errors}
\medskip
\label{app}
\small
\begin{tabular}{llll}
\toprule
Covariate & IFE & IFE-Analytical & IFE-Bootstrap \\
\midrule
$\log(\text{SpillTech})$ & 0.139*** & 0.222*** & 0.132** \\
& (0.044) & (0.044) & (0.047) \\
\addlinespace
$\log(\text{SpillSIC})$ & 0.065 & 0.060 & 0.045 \\
& (0.039) & (0.039) & (0.044) \\
\addlinespace
$\log(\text{R\&D Stock})$ & 0.269*** & 0.255*** & 0.287*** \\
& (0.074) & (0.074) & (0.082) \\
\addlinespace
$\log(\text{Sales})$ & 0.440*** & 0.420*** & 0.378*** \\
& (0.044) & (0.044) & (0.049) \\
\bottomrule
\multicolumn{4}{p{0.62\linewidth}}{\textit{Notes:} Standard errors are reported in parentheses. Standard errors for the IFE-Bootstrap estimator are obtained as the standard deviation across $399$ bootstrap replications. Statistical significance is denoted by ** and *** at the 5\% and 1\% levels, respectively. All explanatory variables are lagged by one period. All specifications include a dummy variable indicating observations with zero lagged R\&D stock. Additionally, we apply the Yeo-Johnson transformation to adjust the significance levels of the bootstrap estimates, where the significance level of $\log(\text{SpillTech})$ changes to $1\%$ while all other covariates remain unchanged.}
\end{tabular}
\end{table}
Consistent with theory, technological spillovers (SpillTech) are positive and statistically significant across all estimation methods, indicating that knowledge diffusion among technologically similar firms increases patenting activity. The estimated magnitude, however, varies with the estimation method: the analytical bias correction yields a noticeably larger coefficient, whereas the bootstrap correction brings it back close to the baseline IFE estimate. This pattern suggests that the analytical correction may over-adjust upward, while the bootstrap provides a more conservative estimate. In contrast, product market spillovers are small and statistically insignificant across all estimators, and this finding is unaffected by either bias-correction approach. This reinforces the interpretation that competitive pressures in product markets do not exert a direct effect on firms' patenting decisions in this setting. The coefficients on firm fundamentals ($R\&D$ stock and sales) from all estimation methods are large and positive, where the bootstrap correction suggests a smaller estimate of sales.
Finally, although \citet{burdaapplication} apply the proposed Bayesian panel probit model on the full unbalanced panel, their application results are broadly consistent with ours, estimated on a sub-balanced panel, in both sign and statistical significance.
Overall, the bootstrap results confirm the robustness of the main conclusions: technological spillovers matter for innovation, product market spillovers do not, and the qualitative inference remains stable across different estimation methods.
\section{Conclusion}
In this paper, we propose a novel bias-correction method for nonlinear panel data models with interactive fixed effects, based on a parametric bootstrap approach. The presence of the incidental parameters problem renders the maximum likelihood estimator biased. We show that, under suitable regularity conditions, the parametric bootstrap replicates the asymptotic distribution of the maximum likelihood estimator. As a result, the bootstrap can be used to perform both bias correction and statistical inference in a unified framework.
Compared with existing bias-correction methods, such as the split-panel jackknife and plug-in analytical approaches, the proposed method does not require an explicit derivation of the bias term and exhibits superior finite-sample performance in panels with relatively small dimensions, as evidenced by our simulation results. In addition, we apply three monotone transformations to reduce skewness in the bootstrap distribution, thereby improving the previously conservative coverage rates.
Finally, we illustrate the empirical relevance of the proposed method through an application to firm-level innovation behavior. The results provide support for key theoretical predictions and highlight the practical usefulness of the method in applied research.
\section*{Supplementary material} \label{supp m}
\paragraph{Replication files:} The replication code for both the simulations and the empirical application is available at \url{https://github.com/Wei-M-Wei/Factor-Bootstrap-replication}.
\bibliographystyle{apalike}
\bibliography{reference2}