Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
98,154 characters · 11 sections · 124 citation commands
Bootstrap Inference in Nonlinear Panel Data Models with Interactive Fixed Effects
Keywords: Panel data, interactive fixed effects, incidental parameter bias, bootstrap inference, skewness correction.
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 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 stock2002forecasting.
Since the seminal contributions of pesaran2006estimation and 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 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 Bai2009,MoonWeidner2015,MoonWeidner2017, the theoretical foundations of nonlinear panel data models with interactive fixed effects are less developed. A notable recent contribution is 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.
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 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, 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 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 hahn2024efficient. Moreover, analytical and jackknife corrections rely on normal approximations to construct confidence intervals, which may be inaccurate.
Recently, 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 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 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 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 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. 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, 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] Related literature
This paper contributes to the literature on large-$T$ bias correction in nonlinear panel data models with fixed effects. Since the seminal work of 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 fernandez2009fixed,HahnKuersteiner2011,DhaeneJochmans2015,arellano2016likelihood,HigginsJochmans2024, additive two-way fixed effects panel models FVWeidner2016, and network models Graham2017,dzemski2019empirical,yan2019statistical,HUGHES2026106130; see FVWeidner2018 for a review. More recently, some studies have begun to address second-order bias DHAENE2021227, schumann2023second as well as higher-order bias correction 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 bai2002determining, onatski2009testing, ahn2013eigenvalue, the estimation of common parameters Bai2009, MoonWeidner2015, MoonWeidner2017, BeyhumGautier2019, BeyhumGautier2022, inference in the presence of weak factors Onatski2012, BaiNg2023, ChoiYuan2025, jiang2025biascorrectionfactoraugmentedregression, ArmstrongWeidnerZeleneev2025, and extensions to network models sassi2024linear; see BaiWang2016 for a comprehensive review. The theoretical development for nonlinear factor models remains comparatively limited. Notable contributions in this area include WANG2022180, who studies the fixed-effect MLE in nonlinear pure factor models; 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 ChenFVWeidner2021, who address bias correction in general nonlinear factor models and introduce an eigenvalue ratio test for factor number selection. More recently, zeleneev2026tractable and yao2025low propose a nuclear-norm penalized estimator that provides computational guarantees, facilitating reliable estimation. boneva2017discrete and 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. GONCALVES2015407 study bootstrap inference for linear dynamic panel models with individual fixed effects. KimSun2016 and HigginsJochmans2024 adopt parametric bootstrap methods for nonlinear panel models with one-way fixed effects. gonccalves2011moving and higgins2025inference introduce moving-block bootstrap procedures for dynamic panel data models. LI2024105684 employ bootstrap methods for inference for treatment effects in interactive fixed-effect panel models. 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 efron1987better,konishi_normalizing_1991,hall_removal_1992.
We consider nonlinear panel data models with interactive fixed effects. The model specification is
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 ChenFVWeidner2021 and 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 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 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 (ref). Interactive fixed effects allow for a much richer structure of unobserved heterogeneity Freyberger2018.
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
where $\mathcal{L}(\beta, \alpha, \gamma; Y, X)$ is the loglikelihood function,
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 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 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
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 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,
which inherits the bias of $(\hat{\beta},\hat{\alpha},\hat{\gamma})$.
ChenFVWeidner2021 study the MLE of model (ref) 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, 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.
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 (ref) by plugging in the MLE \( \hat{\theta} \) and keeping \( X \) fixed:
The bootstrap DGP in (ref) is the same as in 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)), 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 (ref) will also be close to the true DGP, that is, to the DGP in (ref), 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 (ref). 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
where $\hat{\beta}_b^*$ is $\hat{\beta}^*$ computed from the $b$-th bootstrap dataset. Since the first-order bias is eliminated,
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} \),
where $\mathbb{P}_{\theta_0}(\cdot)$ denotes probabilities taken under the DGP (ref). As shown in 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 (ref) 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)$,
where the interval in brackets is a $(1-\alpha)$-level confidence interval for a scalar \( \beta_0 \). The interval can also be written as
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
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 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)) and ((ref)), 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.
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 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 KimSun2016 and HigginsJochmans2024. We first extend the asymptotic expansion for the general M-estimator established in 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 FVWeidner2016 hold uniformly in a small neighbourhood of $\theta_{0}$. Similar uniformity assumptions are imposed in 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 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 FVWeidner2016 and ChenFVWeidner2021, provided that the assumptions are modified by conditioning on $(\alpha_0, \gamma_0)$. See Remark 1 of HahnKuersteiner2011 for more discussion.
Since we generate the bootstrap sample $\{Y^{*}, X\}$ based on model (ref), 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 ChenFVWeidner2021, as some of our assumptions are required to hold uniformly over a neighborhood $\Theta_0$ of the true parameter value.
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|. \]
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 ChenFVWeidner2021.} The conditional independence in Assumption (ref)(i) is commonly adopted in the fixed-effects panel data literature 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 GrahamPowell2012, liu2024identification, BotosaruMuris2024.} The mixing condition, also assumed in 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 HigginsJochmans2024.
Assumption (ref), commonly referred to as the strong factors assumption, is widely used in panel interactive fixed effects models Bai2009, MoonWeidner2015, MoonWeidner2017, ChenFVWeidner2021, GaoLiuPengYan2023.\footnote{For discussion of weak factors, see 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), a uniform version of the expansion in 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) is sufficient to establish the required uniform consistency result. For a formal proof, see Lemma (ref).
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.
Assumption (ref) 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$.
In Assumption (ref), the dominating function \(L_{it}(Z_{it},\theta) \) may depend on \( \theta\). Unlike 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 ChenFVWeidner2021 and GaoLiuPengYan2023.\footnote{In 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) implies point identification for any \( \theta \in \Theta_0 \) that might be the true value, as shown in 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)--(ref), we show in Appendix (ref) that the following expansion holds uniformly in $\theta\in\Theta_0$:
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).
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:
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 HigginsJochmans2024.
Given Assumptions (ref)--(ref), 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 HigginsJochmans2024.
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.
The following corollary establishes the asymptotic validity of the transformed bootstrap confidence interval. It follows from the bootstrap delta method 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 higgins2025inference.
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 (ref).
Given Assumptions (ref)--(ref), the following uniform asymptotic expansion holds for APEs:\footnote{See Theorem (ref) 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)$.
We can now show bootstrap consistency for APEs.
We focus on APEs of the form (ref), 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 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
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 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{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 (ref) is $\sqrt{NT/(N+T-1)}$, see Remark 4 in 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 FVWeidner2016.
Since the loglikelihood function in (ref) 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. 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, zeleneev2026tractable extend the nuclear-norm penalized estimator proposed by 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. 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 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:
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 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
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 (ref) 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
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 (ref) 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 zeleneev2026tractable, that is, using (ref), thereby maintaining consistency with their theoretical framework. The final choice of the tuning parameter is then $\phi=\hat{\varphi}$, to be used in (ref) 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 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
where $S_\beta$, $S_\alpha$, and $S_\gamma$ denote the step sizes for the corresponding parameter blocks; see 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 zeleneev2026tractable. We compute $\hat{\theta}$ using the R package NNRPanel developed by 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 ahn2013eigenvalue; see also ChenFVWeidner2021, GaoLiuPengYan2023, and zeleneev2026tractable.
The following algorithm summarizes our bootstrap inference procedure for $\beta_0$.
Following ChenFVWeidner2021, the two bias-corrected methods can also be applied to the two-step MLE estimator in 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 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:
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) and (ref) 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)--(ref). 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 box1964analysis, and the Yeo--Johnson transformation 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)--(ref) 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.
In this section, we revisit the dataset originally from bloom and studied by 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 burdaapplication, the innovation and patenting behavior are also likely to be affected by unobserved heterogeneity at the firm and time level.
In 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 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 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 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 ahn2013eigenvalue, slightly modified as in 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).
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 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.
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.
\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}.