Robust Priors in Nonlinear Panel Models with Individual and Time Effects

Robust Priors in Nonlinear Panel Models with Individual and Time Effects

\setcounter{footnote}{0}

Introduction

Nonlinear panel models with additive unobserved individual and time effects are widely used to control for persistent heterogeneity and aggregate shocks. Yet in nonlinear models, fixed effects estimation that treats these unobserved effects as nuisance parameters to be estimated suffers from the incidental parameters problem neyman1948consistent. Estimating $N+T$ nuisance parameters from $NT$ observations produces a first-order bias of order $\max\{N^{-1},T^{-1}\}$ in the other common model parameters, and this bias propagates to structural functions such as average partial effects (APEs).

Likelihood-based bias corrections modify the criterion function to produce first-order debiased estimators of model parameters. Among them, the integrated likelihood approach, which averages the likelihood over fixed effects using a weight function (or prior), is particularly attractive. Because many structural quantities, including APEs, are themselves defined through integration, debiasing the model parameters carries over naturally into bias reduction for partial effects within the same integrated-likelihood framework. Moreover, by an appropriate choice of the prior, the same device can also absorb the higher-order bias from estimating the fixed effects, which makes the subsequent APE correction especially transparent. In one-way panels, this approach is tractable because the likelihood separates at the unit level, allowing stable evaluation and debiasing arellano2009robust,pakel2019bias. With two-way effects, however, separability is lost and the integration becomes $(N+T)$-dimensional yielding the classical Laplace expansions invalid. As a result, despite its conceptual simplicity, a rigorous integrated-likelihood debiasing method for two-way panels has been missing.

This paper develops such a framework under the large-$N$, large-$T$ asymptotic scheme of fernandez2016individual. Our starting point is the maximum integrated-likelihood estimator, obtained by integrating out the $(N+T)$ fixed effects. Since direct evaluation of this integral is infeasible in nonlinear two-way panels, we implement the correction through a bias-reducing adjustment to the objective function, which can be interpreted as a robust prior for the integrated likelihood. Our analysis is entirely frequentist and characterizes the sampling behavior of the resulting estimators.

Justifying this strategy requires a new asymptotic theory. With two-way fixed effects, the integrated likelihood involves a high-dimensional integration, and standard Laplace expansions for posterior moments can fail when the parameter dimension grows with the sample size shun_mccullagh_1995. We resolve this difficulty by developing a target-centered full-exponential Laplace approximation that expands directly around the true parameter values, rather than the sample modes, and isolates the leading bias created by integrating out the fixed effects. This approximation turns bias reduction into a constructive problem: choose a prior whose gradient offsets the leading score bias induced by integrating out the fixed effects.

To validate the approximation, we develop a cumulant-based expansion that exploits sparsity in additive two-way models. For each observation, only few third- and higher-order log-likelihood derivative terms are nonzero. This sparsity limits the tensor contractions that drive the high-order remainder even as the fixed-effects dimension diverges. The resulting expansion yields asymptotic linear representations for posterior means of the full model parameter vector and for the associated APE estimators.

These representations have two practical implications. First, they guide the construction of bias-reducing priors for both model parameters and fixed effects. In particular, the prior removes the leading term of the bias for the model parameters and fixed effects. A key message is that debiasing the fixed effects is essential for debiasing APEs. Second, the same expansion characterizes APE estimators under the same priors. Our results show that, even after applying the bias-reducing priors, APE estimators generally retain an additional first-order component driven by fixed-effect estimation uncertainty. We show that this remaining bias has a simple variance form and provide a closed-form correction.

We also clarify the conceptual relationship between integrated and joint likelihood corrections. It is often assumed that priors for the integrated likelihood and penalties for the joint likelihood are interchangeable, since both can be interpreted as Bayesian priors. We show that integration introduces an additional curvature term, so the bias-reducing prior for the integrated likelihood generally differs from the penalty that debiases the joint likelihood.\footnote{A notable exception arises in dynamic AR($p$) panels: we find the same closed-form object delivers both the integrated-likelihood robust prior and the joint-likelihood penalty.} We provide the corresponding joint-likelihood penalties as a companion result.

Finally, we derive model-specific priors and provide a Guidance for Practitioners for workhorse specifications. While the generic robust priors work broadly, exploiting specific model priors yields improvements in finite-sample performance. With strictly exogenous regressors, exponential-family panels (e.g., binary logit and Poisson) admit closed-form, non-data-dependent priors that are particularly stable in finite samples. When regressors are predetermined, as in dynamic specifications with lagged outcomes, we propose a dedicated prior that improves small-$N$ performance by absorbing the cross-sectional component of the leading bias. When $T$ is also small, a simple add-on adjustment addresses the remaining cross-time component. For dynamic AR($m$) panels, we derive an explicit fixed-$T$ consistent prior based on the roots of the AR polynomial, and we show that the same likelihood-based correction carries over to general AR($m$) specifications. For probit panels, comparable closed-form simplifications are unavailable, but the generic priors remain valid. Beyond these cases, the approach extends directly to richer likelihoods: we implement it for multinomial and ordinal response panels with two-way fixed effects and, to our knowledge, provide the first debiasing results for these models in the presence of two-way effects.\footnote{ Recently, honore2025dynamic develop fixed-$T$ identification and GMM estimation for a dynamic ordered logit with individual fixed effects. Other existing work adopt random effects albarran2019estimation or focus on large-$(N,T)$ debiasing in models with one-way effects carro2014state,fernandez2017evaluating. Our focus is large-$(N,T)$ debiasing under two-way effects. }

On the computation side, we propose a scalable Metropolis--Hastings algorithm with random blocking and self-adaptive proposals. This algorithm efficiently handles the high-dimensional parameter space, making integrated-likelihood inference practical for real-world datasets. We provide a Python package, twowaypanel, which covers the main model classes studied in our simulations and empirical illustration.\footnote{ The twowaypanel package (source code and the replication archive) is available at \url{github.com/zizhongyan/twowaypanel}. Online documentation is maintained at \url{twowaypanel.readthedocs.io}.}

Literature review. This paper relates to the literature on the incidental parameters problem in nonlinear panels. While early work emphasized short panels and focused on fixed-$T$ consistency chamberlain1980analysis,arellano1991some, a complementary line adopts the large-$N,T$ framework in which fixed-effects estimators are consistent but display non-negligible first-order bias hahn2004jackknife,hahn2011bias,dhaene2015split,kim2016bootstrap,alvarez2022robust,higgins2024bootstrap. For nonlinear models with individual and time effects, our paper relies on the rigorous asymptotic framework of fernandez2016individual,chen2021nonlinear; see also fernandez2018fixed for a broader review. Our contribution connects to three main strands of this literature.

The first strand develops bias corrections based on the integrated likelihood in one-way panels. lancaster2002orthogonal and arellano2009robust characterize bias-reducing priors and show how modified likelihood objectives can be interpreted as integrated likelihoods; related perspectives appear in woutersen2002robustness and bester2005bias. pakel2019bias extends robust-prior ideas to models with state dependence and schumann2021integrated proposes a debiased method based on transforming the likelihood function. We deal instead with two-way panels, work with the standard likelihood, and develop corrections for both structural parameters and APEs.

Second, a related literature debiases the joint/profile likelihood via penalties, including the general penalized likelihood principle firth1993bias,bester2005bias, panel-specific likelihood adjustments bester2009penalty,arellano2016likelihood, and second-order refinements dhaene2021second. Our analysis shows that, with two-way effects, the bias-reducing prior for the integrated likelihood generally differs from the penalty that debiases the joint likelihood, and we derive the latter as a companion result.

The third strand addresses models with two-way fixed effects. Existing solutions include the analytical and jackknife bias corrections of fernandez2016individual, the differencing method of charbonneau2017multiple for logit/Poisson models, and the recent likelihood-based corrections for dynamic heterogeneous panels by leng2025debiased. jochmans2019likelihood propose a modified likelihood for two-way models and illustrate an MCMC implementation. While their discussion does not develop a formal large-$N,T$ characterization of the resulting posterior mean, our results provide such a characterization within an integrated-likelihood framework. We formally show that posterior-mean bias reduction is tied to the prior associated with the integrated likelihood, whereas bias reduction for the maximized penalized joint likelihood requires a different correction term. This distinction also underlies our bias-corrected APE procedures for two-way panels.

Outline of the paper. The rest of the paper is organized as follows. Section (ref) introduces the model and the integrated likelihood framework. Section (ref) develops the target-centered approximation and robust priors, and derives large-sample distributions for parameters and APEs. Section (ref) provides model-specific priors and practical guidance. The remaining sections present simulation evidence and an empirical illustration, followed by concluding remarks. Proofs and computational details are collected in the Appendix.

Model and Estimators

Model

Let the observed data be denoted by $\left\{\left(Y_{it},X_{it}\right): i=1,\ldots,N, t=1,\ldots,T\right\}$, where $Y_{it}$ is the outcome variable and $X_{it}$ is a vector of explanatory variables. Observations are independent across $i$ and weakly dependent in $t$. The data are generated sequentially over time according to \[ Y_{it}\mid X_i^t,\alpha,\gamma,\beta \sim f_Y\left(\cdot | X_{it};\ \beta, \alpha_i+\gamma_t\right), \quad i=1,\ldots,N,\ t=1,\ldots,T, \] where $f_Y$ is a known probability function, $X_i^t=(X_{i1},\ldots,X_{iT})$, $\beta$ collects the finite dimensional model parameters of interest, and $(\alpha',\gamma')'=: \phi$ stacks individual and time effects with $\alpha=(\alpha_1,\ldots,\alpha_N)'$ and $\gamma=(\gamma_1,\ldots,\gamma_T)'$. The regressors $X_{it}$ may be strictly exogenous or predetermined with respect to $Y_{it}$. Note that $X_{it}$ can include lags of $Y_{it}$ to allow dynamics.

We adopt the fixed effects approach for estimation, treating unobserved individual and time effects as parameters to estimate. The joint log–likelihood is \[ \mathcal L(\beta,\phi) =\tfrac{1}{\sqrt{NT}}\left\{\textstyle\sum_{i=1}^N\sum_{t=1}^T \ln f_Y\left(Y_{it}\,\middle|\,X_{it};\ \beta, \alpha_i+\gamma_t\right) -\frac{b}{2}\,(v'\phi)^2\right\}, \] where $v=(\iota_N',-\iota_T')'$ and $b>0$ is an arbitrary constant. The quadratic term $-\frac{b}{2}(v'\phi)^2$ is included as a normalization device for the two-way effects. Adding a constant to all $\alpha_i$ and subtracting it from all $\gamma_t$ leaves $\alpha_i+\gamma_t$ unchanged, so $\phi$ is unidentified in one location direction. Following fernandez2016individual, this term selects a normalized representative (e.g., $v'\phi=0$) and yields a unique maximizer in $\phi$ without affecting $\beta$.\footnote{Equivalently, one may impose $\sum_i \alpha_i=\sum_t \gamma_t$ or fix one effect such as $\alpha_1=0$ or $\gamma_1=0$).} The factor $1/\sqrt{NT}$ is a normalization that leaves maximizers unchanged and simplifies the stochastic expansions used below. Our main focus is on integrated likelihood corrections. Given a prior $p(\beta,\phi)$, the integrated log–likelihood is \[ \mathcal L^{I}(\beta) =\tfrac{1}{\sqrt{NT}}\ln \left(\textstyle\int \exp\!\big(\sqrt{NT}\,\mathcal L(\beta,\phi)\big)\,p(\beta,\phi) \, \mathrm{d}\phi. \right) \]

In parallel, one may add a penalty directly to the joint likelihood. Writing $p^J(\beta,\phi)$ for a prior (or penalty) used with $\mathcal L(\beta,\phi)$, the penalized objective is $\mathcal L(\beta,\phi)+\sqrt{NT}\,\log p^J(\beta,\phi).$\footnote{Section (ref) will return to the joint–likelihood route, and review the robust penalty used with the joint likelihood and clarify how it differs from the robust prior for the integrated likelihood. In Section (ref), we give model–specific penalty forms that complement the model–specific priors.} From a probabilistic standpoint, both $p(\beta,\phi)$ and $p^J(\beta,\phi)$ can be read as priors.\footnote{ The marginal posterior of $\beta$ under a prior $p(\beta, \phi)$ satisfies $ f(\beta|Y,X) \propto \int f (Y|\beta,\phi;X) p(\beta,\phi) \mathrm{d} \phi $, which is precisely the integrated likelihood weighted by $p(\beta, \phi)$. A Bayesian maximum a posteriori estimator based on the joint likelihood uses $ f(\beta,\phi|Y,X) \propto f(Y|\beta,\phi;X) \times p^{J}(\beta,\phi),$ so $p^J$ plays the role of a prior. Thus both $p$ and $p^J$ enter as priors: maximizing the penalized log-likelihood is a MAP problem, while maximizing $\mathcal{L}^I(\beta)$ corresponds to integrating the prior out. } When the purpose is bias reduction, however, Section (ref) shows that the prior that removes the leading bias in the integrated likelihood generally differs from the penalty for the joint likelihood, due to the unique structure of the panel data model. We therefore keep the notational distinction.

We define the true parameters $(\beta_0,\phi_0) $ as the solution to the conditional population problem $ \max_{\beta,\phi} \mathbb{E}_{\phi}\left[ \mathcal{L}(\beta,\phi)\right], $ where $\mathbb E_{\phi}$ denotes the expectation taken with respect to the distribution of the data conditional on the unobserved effects and initial conditions (including strictly exogenous regressors).

Finally, the quantities of primary economic interest are structural functions, most notably, the average partial effect (APE). Let $\Delta(X_{it},\beta,\alpha_i+\gamma_t)$ be a per observation partial effect that depends on the fixed effects only through $\alpha_i+\gamma_t$. The sample APE at $(\beta,\phi)$ is \[ \Delta(\beta,\phi)=\tfrac{1}{NT}\textstyle\sum_{i=1}^N\sum_{t=1}^T \Delta(X_{it},\beta,\alpha_i+\gamma_t), \] and the population target APE is $\mathbb E[\Delta(\beta_0,\phi_0)]$, where the expectation $\mathbb E$ integrates over both the data and unobserved effects, provided that the expectation exists.

Estimators and likelihood-based corrections

In nonlinear panel models, estimating high-dimensional fixed effects induces an incidental parameter bias of order $\max\{T^{-1},N^{-1}\}$ in model parameter estimators fernandez2016individual. For any fixed $\beta$, we distinguish the sample mode of the fixed effects $ \widehat\phi(\beta)=\arg\max_{\phi}\ \mathcal L(\beta,\phi), $ from the target mode $ \overline\phi(\beta)=\arg\max_{\phi}\ \mathbb E_{\phi}\!\left[\mathcal L(\beta,\phi)\right]. $ The concentrated and target log–likelihoods are then $\mathcal L\!\big(\beta,\widehat\phi(\beta)\big)$ and $\mathcal L\!\big(\beta,\overline\phi(\beta)\big)$, respectively. Under standard conditions, $\widehat\phi(\beta)$ converges to $\overline\phi(\beta)$ as $N,T\to\infty$, but its convergence rate is too slow relative to the sample size $NT$, and induces a bias of order $\max\{T^{-1},N^{-1}\}$ in $\widehat\beta=\arg\max_{\beta}\mathcal L(\beta,\widehat\phi(\beta))$. In contrast, the hypothetical maximizer based on the infeasible target mode, $\arg\max_{\beta}\ \mathcal L\big(\beta,\overline\phi(\beta)\big)$, does not suffer from this first–order bias. These target quantities are our benchmark and guide the bias corrections that follow.

Our main focus is the integrated likelihood. The associated maximum integrated likelihood estimator (MILE) is $\widehat\beta^{I} = \arg\max_{\beta}\ \mathcal L^I(\beta)$. We construct bias‐reducing priors under which $\widehat\beta^{I}$ is free of first–order bias. Direct evaluation of $\mathcal L^{I}(\beta)$ is infeasible with two–way effects because it requires integration over an $(N{+}T)$–dimensional $\phi$.\footnote{ With one–way effects the likelihood factorizes by unit and the integral can be computed at the “atom” level, which admits stable quadrature arellano2009robust,pakel2019bias. Two–way effects couple individuals and time periods and break this factorization, so the integral is high dimensional. } We therefore work with the joint posterior for $(\beta,\phi)$ induced by $p(\beta,\phi)$, and propose to use the posterior means\footnote{ Posterior mode is a possible alternative. However, computing them typically requires an additional nonparametric approximation, and variance estimation is more involved. }

Under the priors developed in Sections (ref)–(ref), $\widehat{\beta}^E$ removes the leading bias without inflating asymptotic variance. In contrast to standard corrections that target $\beta$ only, the same priors eliminate the leading bias in $\widehat{\phi}^E$, which is useful for bias correction of nonlinear functionals such as APEs.

Under the integrated–likelihood setup, we consider two APE estimators. The first estimator is the posterior mean APE\footnote{The estimator $\widehat\Delta^E$ is referred to as the Bayesian fixed effect APE in arellano2009robust. Related posterior-based average effects are studied by bonhomme2022posterior, which has a different focus on posterior conditioning and distributional robustness in the random effects setting.} $$ \widehat\Delta^E =\tfrac{\iint \Delta(\beta,\phi)\ \exp\!\big(\sqrt{NT}\ \mathcal L(\beta,\phi)\big)\,p(\beta,\phi)\,d\phi\,d\beta} {\iint \exp\!\big(\sqrt{NT}\ \mathcal L(\beta,\phi)\big)\,p(\beta,\phi)\,d\phi\,d\beta}. $$ The second is the posterior plug–in APE, $\Delta(\widehat\beta^{E},\widehat\phi^{E})$, which evaluates the APE at the posterior means $(\widehat\beta^{E},\widehat\phi^{E})$. Section (ref) shows that, once both leading biases in estimating $(\beta,\phi)$ are removed via the integrated-likelihood based correction, the remaining APE bias has a simple variance–type form that can be subtracted in closed form.

We approximate the posterior integrals by Markov chain Monte Carlo (MCMC). Appendix (ref) describes a scalable random-blocking, self-adaptive Metropolis--Hastings algorithm for the high-dimensional two-way setting.

For completeness we also consider penalized joint likelihood. The maximum penalized likelihood estimator is $(\widehat{\beta}^P, \widehat{\phi}^P) = \mathop{\rm arg\,max}_{\beta,\phi}\big(\mathcal L(\beta,\phi) + \sqrt{NT}\ln p^J(\beta,\phi)\big).$ One can analogously report the plug–in APE estimator $\Delta(\widehat\beta^{P},\widehat\phi^{P})$. Section (ref) will discuss model‐specific penalties that parallel our bias-reducing priors.

Notation and Assumptions

We write $\mathcal S(\beta,\phi)=\partial_{\phi}\mathcal L(\beta,\phi)$ for the score in the fixed–effects direction and $\mathcal H(\beta,\phi)= - \partial_{\phi\phi'}\mathcal L(\beta,\phi)$ for the corresponding negative Hessian, where the notation $\partial_{x}f$ denotes the partial derivative of $f$ with respect to $x$, and additional subscripts denote higher-order partial derivatives. Bars denote conditional expectations given $\phi$, e.g., $\partial_{\beta}\overline{\mathcal L}=\mathbb E_{\phi}[\partial_{\beta}\mathcal L]$, and tildes denote deviations, e.g., $\partial_{\beta}\widetilde{\mathcal L}=\partial_{\beta}\mathcal L-\partial_{\beta}\overline{\mathcal L}$. When no confusion arises, arguments are omitted and functions are understood to be evaluated at the truth $(\beta_0,\phi_0)$, e.g. $\mathcal{L}= \mathcal{L}(\beta_0,\phi_0)$. Additionally, let $\pi_{it}=\alpha_i+\gamma_t$, we define the per–observation log–likelihood as $$ \ell_{it}(\beta,\pi_{it})=\ln f_Y\left(Y_{it}\,\middle|\,X_{it};\ \beta, \alpha_i+\gamma_t \right). $$

We impose the following conditions, which are identical to those in fernandez2016individual except for the last condition of Assumption (ref) on the prior.\footnote{ Assumption (ref)(vi) is a mild high-level requirement ensuring that the prior acts as a bias-correction term and does not alter the leading stochastic orders of the score and Hessian of the criterion function. }

Characterization of the Robust Priors

This section develops a target-centered Laplace approximation, uses it to construct bias-reducing priors for generic nonlinear panels, and derives the large-sample distributions for parameters and APEs. Section (ref) then specializes these priors to widely used models and provides practical guidance.

Target-Centered Full-Exponential Laplace Approximations

We analyze the integrated likelihood and develop a new Laplace approximation that is pivotal for constructing bias-reducing priors for both the model parameters and the APE. Our expansion is centered at the target quantities, i.e., the true parameter values, rather than at sample modes (MLE). This choice significantly simplifies the derivations and provides a clean asymptotic framework. In contrast to classical Laplace methods tierney1986accurate,tierney1989fully and to arellano2009robust and pakel2019bias in the panel data context, which expand around the MLE, we expand around the truth. No smoothness beyond the regularity in fernandez2016individual is required. The next lemma links the integrated and target likelihoods:

The term $\mathcal{{D}}_{\mathcal{L}^I}\left( \beta,\overline{\phi}(\beta) \right) $coincides with the leading bias in $\mathcal{L}^I(\beta)$ and is $ \mathcal{O}_P(\sqrt{NT})$. Hence, if $\frac{1}{\sqrt{NT}} \left[ \ln p\left(\beta, \overline{\phi}(\beta)\right) - \mathcal{D}_{\mathcal{L}^I}\left( \beta,\overline{\phi}(\beta) \right) \right] = o_P(1)$, then the first-order bias in the integrated log-likelihood vanishes and the maximum integrated likelihood estimator $\widehat{\beta}^{I}$ is free of first-order bias. In the two-way case, direct evaluation of $\mathcal{L}^{I}(\beta)$ becomes intractable, so we work instead with the posterior means $(\widehat{\beta}^E,\widehat{\phi}^E)$. Building on Lemma (ref), we next develop expansions for these posterior means centered at the truth $(\beta_0,\phi_0)$.

Validating this expansion in the two-way setting is challenging. shun_mccullagh_1995 demonstrate that standard Laplace approximations can fail unless the parameter dimension grows slower than the cube root of the sample size. This condition is violated in our context.\footnote{ Under Assumption (ref)(i), the parameter dimension grows at rate $\mathcal{O}((NT)^{1/2})$. This exceeds the rate $o((NT)^{1/3})$ required by shun_mccullagh_1995 for the convergence of remainder terms. } To overcome this, we propose a cumulant-based expansion that exploits the extreme sparsity of higher-order derivatives in two-way additive models. This sparsity restricts tensor contractions and keeps remainder terms under control. We obtain the following asymptotic representations for $(\widehat{\beta}^E,\widehat{\phi}^E)$:

Within the integrated-likelihood framework, both the posterior mean APE $\widehat{\Delta}^E$ and the posterior plug-in APE $\Delta(\widehat{\beta}^E,\widehat{\phi}^E)$ admit comparable asymptotic representations around the true APE $\Delta := \Delta(\beta_0,\phi_0)$:

Robust Priors for Generic Models

We turn to a general formulation of the bias-reducing prior. Our aim is straightforward: choose a prior $p(\beta,\phi)$ so that the posterior means $(\widehat{\beta}^E,\widehat{\phi}^E)$ have no first-order bias. Theorem (ref) shows that this bias is driven by the gap between the gradient of the log–prior and the gradient of the leading bias term in the integrated log–likelihood.

To make this decomposition explicit, let $\Upsilon_{it}^{\pi}\in\mathbb{R}$ and $\Upsilon_{it}^{\beta}\in\mathbb{R}^{\dim\beta}$ denote the observation level contributions to the leading bias of gradients in the fixed-effects and common parameter directions, respectively. In the proof of Theorem (ref) we show that

Here,

To discuss the last term in (ref) jointly, let $\theta\in\{\beta,\pi\}$, where $\partial_{\theta\pi}\ell_{i\tau}$ denotes $\partial_{\beta\pi}\ell_{i\tau}$ when $\theta=\beta$ and $\partial_{\pi^2}\ell_{i\tau}$ when $\theta=\pi$. Then the term $\sum_{\tau=t+1}^T\mathbb E_\phi(\partial_{\theta\pi}\ell_{i\tau}\,\partial_\pi\ell_{it})$ captures cross-time dependencies. It is a spectral expectation that is generally nonzero in dynamic or predetermined settings, and as discussed in Section (ref), it largely determines the fixed-$T$ performance. Following the standard spectral estimation, we replace $\sum_{\tau=t+1}^T\mathbb E_\phi(\partial_{\theta\pi}\ell_{i\tau}\,\partial_\pi\ell_{it})$ by $\sum_{\tau=1}^{\min\{L,T-t\}}T/(T-\tau) \mathbb E_\phi(\partial_{\theta\pi}\ell_{i,t+\tau}\,\partial_\pi\ell_{it})$ for a fixed integer $L$, and apply the degrees-of-freedom adjustment $T/(T-\tau)$ as in hahn2011bias.\footnote{ Here $L$ is the trimming parameter for estimating spectral expectations and we assume $L\rightarrow \infty$ and $L/T\to 0$ asymptotically. }

Representation (ref) characterizes an equivalence class of bias-reducing priors, that is, any log-prior whose gradient matches $\Upsilon$ at the below stated rates will eliminate the leading bias in both $\widehat{\beta}^E$ and $\widehat{\phi}^E$:

Among the priors that satisfy (ref), we report two convenient closed-form choices:

Implementation requires evaluating the expectations $\mathbb{E}_{\phi}(\cdot)$, which are normally data-dependent. Following arellano2016likelihood, the expectations can be computed either in observed form (sample analog) or in closed form from the model. We adopt a mixed strategy: use observed forms by default (consistent for the expectations), and substitute analytical expectations whenever available to improve small-sample performance. The large-$N,T$ validity of the correction is unaffected, and our Monte Carlo evidence confirms that the two approaches coincide asymptotically.

In Section (ref) we show that, for specific models, certain components of $\Upsilon$ (or of their expectations) vanish, which delivers closed-form simplifications of the prior. To our knowledge, these are the first model-specific bias-reducing priors.

Asymptotic Distribution and Inference

The next results summarize our asymptotics under bias-reducing priors. The first result describes the limiting distribution of the posterior means, while the second treats the two APE estimators and isolates a residual variance-type bias term:

Before stating the asymptotic result for APE estimators, it is useful to emphasize that our APE estimators center at the population target $\mathbb E(\Delta)$ rather than at the sample average $\Delta$. As a consequence, the convergence rate need not be the parametric $\sqrt{NT}$ rate, and instead depends on the sampling properties of $(X_{it},\alpha_i,\gamma_t)$, as shown by fernandez2016individual.

To our knowledge, there exists no tractable prior that simultaneously corrects both parameter and APE biases. Therefore, we address the remaining bias term directly. We construct the bias-corrected APE estimators by estimating $\mathcal{B}^{\Delta}$ at the posterior means and performing the following subtraction: $$ \widetilde{{\Delta}}^E = \widehat{\Delta}^E -2{\mathcal{\widehat{B}}}^{\Delta}( \widehat{\beta}^E, \widehat{\phi}^E ), \qquad \widetilde{\Delta}( \widehat{\beta}^E, \widehat{\phi}^E ) = \Delta( \widehat{\beta}^E, \widehat{\phi}^E ) -{\mathcal{\widehat{B}}}^{\Delta}( \widehat{\beta}^E, \widehat{\phi}^E). $$

For inference on model parameters and APEs, one may either (i) compute analytical asymptotic variances (i.e., $\overline{W}^{-1}_{\infty}$, $\overline{\mathcal{H}}^{-1}_{\infty}$ and $\overline{V}_{\infty}^{\delta}$), or (ii) use posterior draws. The latter is justified by the Bernstein–von Mises theorem: although our priors are improper, Bernstein–von Mises applies under local positivity and continuity of the prior near the true parameter van2000asymptotic, and implies asymptotic equivalence. Our assumption ensures these conditions, so both routes lead to the same first-order inference.

Bias-reducing Penalties for Joint Likelihood

Together with jochmans2019likelihood, the connection between integrated- and joint-likelihood corrections is direct. The penalty term is given by Prior GE minus $\frac{1}{2}\ln{\det\left( \overline{\mathcal{H}}\left( \beta,\phi \right) \right)}$.\footnote{ To see the link, the expansion of the concentrated log-likelihood around the target mode $\overline{\phi}(\beta)$ gives the leading bias $-\tfrac{1}{2}\,\mathrm{tr}\!\left\{\overline{\mathcal H}^{-1}(\beta,\overline{\phi}(\beta))\,\sqrt{NT}\,\mathcal S(\beta,\overline{\phi}(\beta))\,\mathcal S(\beta,\overline{\phi}(\beta))^{\prime}\right\}.$ When moving to the integrated log-likelihood, the standard Laplace expansion about the sample mode $\widehat{\phi}(\beta)$ (cf. tierney1989fully) contributes an additional term $\frac{1}{2}\ln{\det\left(\cdot \right)}$. Hence the joint–likelihood penalty is obtained by subtracting this term: $\mathrm{Penalty GE}=\mathrm{Prior GE}-\frac{1}{2}\ln{\det\left(\cdot\right)}$. } Specifically, we obtain detailed penalty formulas for generic nonlinear panels: \footnote{Proof sketch. Start from the leading bias of the concentrated log-likelihood in jochmans2019likelihood, and apply the same decomposition steps as in the proof of Theorem (ref) to collapse this matrix expression into scalar form, which yields Penalty GEs in closed form. A step-by-step derivation note is available upon request.}

jochmans2019likelihood illustrate an MCMC scheme that targets the joint likelihood combined with Penalty GE, and report the corresponding posterior mean of $\beta$. Our theory distinguishes the roles of Penalty GE and Prior GE. Penalty GE is designed to remove the leading bias of the maximized penalized joint- likelihood estimator. For MCMC-based posterior means, in contrast, Corollary (ref) shows that bias reduction is achieved when the likelihood is paired with Prior GE.

Given $(\widehat{\beta},\widehat{\phi})$ from the penalized joint likelihood, the APE is computed in plug-in form and can be bias-corrected in the same way: $ \widetilde{\Delta}( \widehat{\beta}, \widehat{\phi} ) = \Delta( \widehat{\beta}, \widehat{\phi} ) -{\mathcal{\widehat{B}}}^{\Delta}( \widehat{\beta}, \widehat{\phi} ). $

Model-Specific Bias-reducing Priors

We begin with two workhorse specifications in applied microeconometrics: binary logit and Poisson counts. Both models are one-parameter exponential-family models.\footnote{ Examples of one-parameter exponential family models include binary logit, Poisson, and normal regression with known variance, among others. Each observation is indexed by a single natural parameter and under canonical links, the variance is pinned down by the mean.} Precisely because of this structure, the generic priors simplify to non–data-dependent formulas in the settings with strict exogenous covariates. We then carefully cover logit and Poisson models with predetermined covariates, and other widely used specifications such as dynamic AR($m$) panels, multinomial logit, and probit, allowing for either strict exogeneity or predetermined covariates.

Our aims are twofold. First, where closed-form simplifications exist, we provide explicit priors that improve small sample behavior and are easy to compute. Second, where such simplifications are unavailable, we indicate when the generic priors from Section (ref) are the appropriate default and when score- or estimator-based corrections may be preferable.

One-Parameter Exponential Family Panels

We describe one-parameter exponential family panel models with additive fixed effects, following mccullagh2019generalized. Let the conditional density of $Y_{it}$ belong to the exponential family with natural index $X_{it}^{\prime}\beta+\pi$, and write the log-likelihood in its canonical form

for known functions $\mathcal{C}_1$ and $\mathcal{C}_2$. Then, $\partial_{\pi}\ell_{it} = Y_{it} - \partial_{\pi}\left[\mathcal{C}_1(X_{it}'\beta_0 + \pi_{it,0})\right]$, $\partial_{\pi^2}\ell_{it} = - \partial_{\pi^2}\left[\mathcal{C}_1(X_{it}'\beta_0 + \pi_{it,0})\right]$, $\partial_{\pi^3}\ell_{it} = - \partial_{\pi^3}\left[\mathcal{C}_1(X_{it}'\beta_0 + \pi_{it,0})\right]$, and the mixed derivatives satisfy the canonical identities $\partial_{\beta\pi}\ell_{it} = -\partial_{\pi^2}\ell_{it} X_{it}$, and $\partial_{\beta\pi^2}\ell_{it} = -\partial_{\pi^3}\ell_{it} X_{it}$. Hence $\partial_{\pi^{2}}\ell_{it}$ and its higher $\pi$-mixed derivatives do not depend on the data $Y_{it}$.

Exponential family: strict exogeneity.

When $X_{it}$ is strictly exogenous, the expressions for $\Upsilon^{\theta}_{it}$ in (ref) collapse to $ \Upsilon_{it}^{\theta} = \frac{\partial_{\theta\pi^{2}}\ell_{it} }{\sum_{\tau = 1}^{T}{\partial_{\pi^{2}}\ell_{i\tau} }} + \frac{\partial_{\theta\pi^{2}}\ell_{it} }{\sum_{j = 1}^{N}{\partial_{\pi^{2}}\ell_{jt} }} , $ for $\theta\in\{\beta, \pi\}$. Because $\partial_{\pi^{2}}\ell_{it}$ and its higher derivatives are $Y$-free, observed and expected forms coincide. Integrating $\Upsilon^{\theta}$ yields a non–data-dependent prior:

(Prior SE: exponential family with strictly exogenous regressors)

This prior exactly cancels the leading bias in the integrated log-likelihood for one-parameter exponential families with strict exogenous covariates. While the generic priors remain valid, Prior SE is preferable in finite samples because it avoids estimated expectations.

Exponential family: predetermined regressors.

Using exponential-family properties, we have that $\mathbb{E}_{\phi}(\partial_{\theta\pi}\ell_{it}\,\partial_{\pi}\ell_{it})=0$. With predetermined variable in $X_{it}$, the cross-time term in (ref) survives: $$ \Upsilon_{it}^{\theta} = \tfrac{\mathbb{E}_{\phi}(\partial_{\theta\pi^{2}}\ell_{it})}{\sum_{\tau = 1}^{T}{ \mathbb{E}_{\phi}( \partial_{\pi^{2}}\ell_{i\tau} ) }} +\tfrac{\mathbb{E}_{\phi}(\partial_{\theta\pi^{2}}\ell_{it})}{\sum_{j = 1}^{N}{\mathbb{E}_{\phi}( \partial_{\pi^{2}}\ell_{jt} )}} + \tfrac{\sum_{\tau=1}^{\min\{L,T-t\}}T/(T-\tau) \mathbb E_\phi(\partial_{\theta\pi}\ell_{i,t+\tau}\,\partial_\pi\ell_{it})} {\sum_{\tau = 1}^{T}{ \mathbb{E}_{\phi}(\partial_{\pi^{2}}\ell_{i\tau} )}}, $$ for $\theta\in\{\beta, \pi\}$. An exact closed-form primitive for the differential system (ref) is unavailable. Nevertheless, the following simplified prior is bias-reducing:

(Prior PE: exponential family models with predetermined regressors)

Relative to the generic priors, Prior PE improves small-$N$ performance and and its finite-sample accuracy does not depend on the individual dimension. Its remaining finite-sample error is driven by the cross-time dependence captured by the trimmed lead--lag term. To see this, consider $\partial_{\alpha_i}\ln p^{PE}$. A direct calculation gives the remaining mismatch between $\partial_{\alpha_i}\ln p^{PE}$ and the target $\textstyle\sum_{t=1}^T \Upsilon^\pi_{it} $:

Both remainder terms involve the trimmed lead--lag moments, i.e., $\mathbb{E}_{\phi}(\partial_{\pi}\ell_{i,t+\tau}\partial_{\beta\pi}\ell_{it})$ and $\mathbb{E}_{\phi}(\partial_{\pi}\ell_{i,t+\tau}\partial_{\pi}\ell_{it})$. By the martingale difference property, these expectations are zero at the truth. In practice, however, they must be estimated by plug-in sample analogs. As $T$ grows, these moments are consistently estimated, and the mismatch vanishes. When $T$ is small, they can be noisy to estimate, so a finite-$T$ discrepancy may remain even under Prior PE.

When $T$ is small and one seeks additional accuracy, a quick add-on correction is available. This alternative is to estimate $(\widehat{\beta}^E,\widehat{\phi}^E)$ under Prior SE and adjust for the cross-time component $$\mathcal{B}_i^{(\theta,1)} = \tfrac{\textstyle\sum_{t = 1}^{T}{\sum_{\tau=1}^{\min\{L,T-t\}}\frac{T}{T-\tau}\mathbb{E}_{\phi}(\partial_{\beta\pi}\ell_{i,t+\tau}\partial_{\pi}\ell_{it})}}{\sum_{t = 1}^{T}{ \mathbb{E}_{\phi}(\partial_{\pi^{2}}\ell_{it}) }} $$ to obtain $$ \widetilde{\beta}^E = \widehat{\beta}^E - \tfrac{1}{\sqrt{NT}} \widehat{{W}}^{-1} \left( \textstyle\sum_{i=1}^{N}\widehat{\mathcal{B}}_i^{(\beta,1)} + [\partial_{\beta\phi'}\widehat{{\mathcal{L}}}]\widehat{{\mathcal{H}}}^{-1} \widehat{\mathcal{B}}^{(\pi,1)} \right), \quad \widetilde{\phi}^E = \widehat{\phi}^E(\widetilde{\beta}^E) - \widehat{{\mathcal{H}}}^{-1}\widehat{\mathcal{B}}^{(\pi,1)}, $$ where ${\mathcal{B}}^{(\pi,1)}=({\mathcal{B}}_1^{(\pi,1)},\ldots,{\mathcal{B}}_N^{(\pi,1)},0_T')'$, and the notation $\widehat{f}=f(\widehat{\beta}^E,\widehat{\phi}^E)$ denotes the plug-in fixed effects estimator.

Poisson panels: strict exogeneity.

After the general analysis of exponential families, we single out the Poisson model as a worked example because we obtain an additional result. In this case, $\mathcal{C}_1\left(X_{it}'\beta+\pi\right) = \exp\left(X_{it}'\beta+\pi\right)$. With $\omega_{it}=\exp\left(X_{it}'\beta_0 + \pi_{it,0}\right) $, we have $\partial_{\pi}\ell_{it}=Y_{it}-\omega_{it}$, $\partial_{\pi^2}\ell_{it}=\partial_{\pi^3}\ell_{it}=-\omega_{it}$, $\partial_{\beta}\ell_{it}=\partial_{\pi}\ell_{it}X_{it}$, and $\partial_{\beta\pi}\ell_{it}=\partial_{\beta\pi^2}\ell_{it}=\partial_{\pi^2}\ell_{it}X_{it}$. When all covariates are strictly exogenous, substituting into Theorem (ref), the bias of $\widehat{\beta}^E$ depends on $ \mathcal{B}^{\beta} + [\partial_{\beta\phi'}\overline{\mathcal{L}}]\overline{\mathcal{H}}^{-1} \mathcal{B}^{\phi} = \frac{\partial_{\beta}\ln p}{\sqrt{NT}} - \frac{1}{N} \sum_{i} \frac{\sum_{t}\omega_{it}X_{it} (\partial_{\alpha_i} \ln p) }{\sum_{t}{ \omega_{it} }} - \frac{1}{T} \sum_{t} \frac{\sum_{i}\omega_{it}X_{it} (\partial_{\gamma_t} \ln p) }{\sum_{i}{ \omega_{it} }} + o_P(1); $ and leading bias of $\widehat{\phi}^E$ depends on $ \mathcal{B}^{\phi_k} = \tfrac{\partial_{\phi_k} \ln p - 2 }{\sqrt{NT}} + o_P((NT)^{-1/2}), $ for $k=1,\ldots,\dim\phi$.

If the goal is only to debias $\widehat{\beta}^E$, a flat prior already removes the leading bias. This finding echoes fernandez2016individual, who show that in the Poisson model with strict exogenous $X_{it}$, the fixed-effects estimator of $\beta$ is asymptotically unbiased. However, our derivations indicates that a flat prior does not remove the second-order bias in $\widehat{\phi}^E$, and as a consequence, complicates the APE debiasing (Corollary (ref)). In contrast, Prior SE eliminates the leading bias in both $\widehat{\beta}^E$ and $\widehat{\phi}^E$ and, for the APE correction, leaves only the variance-driven term in Corollary (ref) to subtract.

Dynamic AR($m$) Panels

Dynamic linear panels provide a transparent benchmark for incidental-parameter bias. The classic “Nickell bias” arises under fixed-$T$ asymptotics nickell1981bias, while explicit bias expansions and higher-order corrections are available under large-$N$, large-$T$ asymptotics kiviet1995bias, hahnkuersteiner2002,hahn2006reducing. Related likelihood-based adjustments for autoregressive panels with fixed effects are studied in dhaenejochmans2016 and alvarez2022robust. This subsection shows that, under our bias-reducing prior framework, a simple closed-form prior is available for Gaussian AR($m$) panels even in the presence of two-way additive effects.

We allow unrestricted time effects $\{\gamma_t\}$ to absorb aggregate shocks. Such time effects are generally incompatible with strict stationarity of the level process $\{Y_{it}\}$ (unless $\gamma_t$ is constant and thus absorbed by normalization). Accordingly, our derivation does not impose stationarity on $\{Y_{it}\}$. Instead, we impose stability of the AR polynomial on the innovation-driven component obtained after removing $\alpha_i+\gamma_t$, which ensures an MA($\infty$) representation and delivers the moment calculations below.\footnote{Throughout, the stability restriction is imposed on $\widetilde Y_{it}=Y_{it}-\alpha_i-\gamma_t$. In estimation, $\gamma_t$ is still included (or profiled out) to absorb common shocks, and it will not enter the leading bias terms for common parameters.}

Consider $ Y_{it} = \mu_{1,0} Y_{i,t-1}+\cdots+\mu_{m,0} Y_{i,t-m} + \alpha_{i,0} + \gamma_{t,0} + \varepsilon_{it}, $ with $ (\varepsilon_{i1},\ldots,\varepsilon_{iT}) \sim\mathcal{N}(0,\sigma_{0}^{2}I_{T}) $ conditional on the true fixed effects and initial conditions. Let $X_{it}=(Y_{i,t-1},\ldots,Y_{i,t-m})^{\prime}$ and $\beta=(\mu^{\prime},\sigma)^{\prime}$, $\mu=(\mu_1,\ldots,\mu_m)'$. The per-observation log-likelihood is given by $ \ell_{it}(\beta,\pi) = -\frac{1}{2}\ln(2\pi)-\frac{1}{2}\ln(\sigma^2)-\frac{1}{2}\frac{(Y_{it}-X_{it}'\mu-\pi)^2}{\sigma^2}. $ Then $\partial_{\pi}\ell_{it} = \varepsilon_{it}\sigma^{-2}$, $\partial_{\pi^2}\ell_{it} = -\sigma^{-2}$, $\partial_{\mu\pi}\ell_{it} = \partial_{\pi^2}\ell_{it}X_{it}$, $\partial_{\sigma\pi}\ell_{it} = -2\sigma^{-1}\partial_{\pi}\ell_{it}$, $\partial_{\sigma\pi^2}\ell_{it} = -2\sigma^{-1}\partial_{\pi^2}\ell_{it}$, and $\partial_{\pi^3}\ell_{it} =\partial_{\mu\pi^2}\ell_{it} = 0$.

Assume the AR polynomial $\Phi(z)=1-\sum_{k=1}^{m}\mu_{k}z^{k}$ is stable. We now work with the innovation driven component $\widetilde Y_{it}=Y_{it}-\alpha_i-\gamma_t$. Let $\widetilde Y_{it}=\sum_{h\ge 0}\psi_h \varepsilon_{i,t-h}$ be the MA($\infty$) representation with $\psi_0=1$ and $\psi_h=\sum_{l=1}^{\min(m,h)}\mu_l\,\psi_{h-l}$ for $h\ge 1$ hamilton2020time. The required cross-moments involve the innovation-driven component only. The fixed effects contribute a deterministic term that is orthogonal to the innovations and therefore drops out from the leading bias expressions.

Standard moments imply, for $k=1,\ldots,m$, $ \sum_{t,\tau=1}^T \mathbb E(\widetilde Y_{i,\tau-k}\varepsilon_{it}) = \sigma_\varepsilon^2\sum_{h=0}^{T-k-1}(T-k-h)\psi_h, $ $ \sum_{t,\tau=1}^T \mathbb E(\varepsilon_{i\tau}\varepsilon_{it}) = T\,\sigma_\varepsilon^2, $ $ \sum_{t<\tau}^T \mathbb E(\varepsilon_{i\tau}\varepsilon_{jt})=0 (\forall i\neq j). $ Substituting into (ref) gives $\Upsilon_{it}^{\pi}=0$ and $\Upsilon_{it}^{\sigma}=0$, while for $k=1,\ldots,m$, we obtain $ \sum_i\sum_t\Upsilon_{it}^{\mu_k} = \sum_{h=0}^{T-k-1}(T-k-h)\psi_h . $ Hence, for this Gaussian AR($m$) specification, the two-way additive effects do not enter the leading bias terms that determine the bias-reducing prior. Equivalently, the same primitive applies whether one includes only individual effects or both individual and time effects, because the relevant moments depend on $\mu$ only through the innovation-driven recursion.

(Prior AR($m$): dynamic AR($m$) panels )

where $r_{1},\ldots,r_{m}$ are the roots of the AR polynomial $\Phi(z)=1-\sum_{k=1}^{m}\mu_{k}z^{k}$. \footnote{Proof sketch. Since $\sum_{t=1}^T\frac{z^t}{t}\sum_{k=1}^{m}r_k^t = - \sum_{k=1}^{m}\ln(1-r_k z)= - \ln \Phi(z)$, the log prior can be written as $\ln p = -[z^0]\{A(z)z^{-m}\ln \Phi(z)\}$, with $A(z)= \sum_{t=0}^{T-1} (T-t) z^t$ and $[z^0]\{\cdot\}$ denoting the constant term in the Laurent expansion (ahlfors1979complex Ch. 4) around $z=0$. For $k=1,\ldots,m$, using $\frac{\partial}{\partial\mu_k}\ln\Phi(z)=z^k/\Phi(z)$ and $1/\Phi(z)=\sum_{t\ge0}^{\infty}\psi_h z^t$ gives that $\partial_{\mu_k}\ln p = [z^0]\left\{A(z) \frac{\partial}{\partial\mu_k}[z^{-k}\ln\Phi(z)]\right\} = [z^0]\{A(z)/\Phi(z)\} =\sum_{t=0}^{T-1}(T-t)\psi_{t-k} = \sum_{h=0}^{T-k-1}(T-k-h)\psi_h =\sum_i\Upsilon_{\alpha,i}^{\mu_k}+\sum_t\Upsilon_{\gamma,t}^{\mu_k}$, which is exactly the desired expression. }

Although Prior AR($m$) is expressed in terms of the characteristic roots, computing those roots via convolution or fast Fourier transformation is unnecessary. The power sums $\sum_{k=1}^{m} r_k^{t}$ can be generated recursively computed using Newton’s identities knuth1997art, which is numerically stable for moderate $m$. This recursion yields a simple and fast implementation for AR models. For AR(2) panel model, two roots can be obtained explicitly from the quadratic formula, yielding a closed-form representation. For AR(1), the single root is simply $r_1=\mu_1$, so we have:

(Prior AR(1): dynamic AR(1) panel models)

Note that Prior AR(1) coincides with lancaster2002orthogonal,arellano2009robust in the one-way case. The invariance to two-way additive effects is specific to this Gaussian linear AR($m$) benchmark. In nonlinear panels, unrestricted time effects can alter the leading bias terms and may require a modified asymptotic expansion fernandez2016individual.

Given the prior, a convenient point estimator is the penalized fixed-effects Gaussian MLE, $ (\hat\mu,\hat\sigma,\hat\alpha,\hat\gamma) \in \arg\max_{\mu,\sigma,\alpha,\gamma} \left\{ \sum_{i=1}^N\sum_{t=m+1}^T \ell_{it}(\beta,\alpha_i+\gamma_t) +\ln p(\mu) \right\}. $ Profiling out $(\alpha,\gamma)$ yields the usual two-way within transformation. Let $\ddot Y_{it}$ and $\ddot X_{it}$ denote the double-demeaned outcome and regressors over $i$ and $t$. Then the bias corrected $\hat\mu$ solves the penalized normal equations: $$ \big(\textstyle\sum_{i,t}\ddot X_{it}\ddot X_{it}'\big)\hat\mu = \sum_{i,t}\ddot X_{it}\ddot Y_{it} +\hat\sigma^{2}\,\partial_{\mu}\ln p(\hat\mu), $$ with $\hat\sigma^{2}$ given by the profiled Gaussian variance (e.g., $\hat\sigma^{2}=(NT)^{-1}\sum_{i,t}\hat\varepsilon_{it}^{2}$). Moreover, the gradient admits the closed form $\partial_{\mu_k}\ln p(\mu)=\sum_{h=0}^{T-k-1}(T-k-h)\psi_h$, where ${\psi_h}$ are computed by the usual AR recursion.

Multinomial Logit Panels

Multinomial logit is a standard specification for polytomous panel outcomes. Building on the one-way fixed effects literature chamberlain1980analysis,wooldridge2010econometric, we consider two-way fixed effects. Our approach extends seamlessly. We assume that $Y_{it}\in\{1,\ldots,J\}$ is a polytomous categorical variable, with category $1$ taken as the baseline. For each non-baseline category $k=2,\ldots,J$, let $ \pi_{it,k}=\alpha_{i,k}+\gamma_{t,k}, $ where $\alpha_{i,k}$ and $\gamma_{t,k}$ are the individual and time effects for category $k$, relative to the baseline. Stack these objects as $\beta=(\beta_2',\ldots,\beta_J')',$ $\alpha_i=(\alpha_{i,2},\ldots,\alpha_{i,J})',$ $\gamma_t=(\gamma_{t,2},\ldots,\gamma_{t,J})',$ and $ \pi_{it}=(\pi_{it,2},\ldots,\pi_{it,J})'=\alpha_i+\gamma_t. $ Accordingly, the fixed effect parameter is now vector valued. The choice probabilities are given by $$ \Pr\left( Y_{it} = k | X_{it},\pi_{it},\beta \right) =

$$ and per-observation log-likelihood is $\ell_{it}(\beta,\pi) = \sum_{k=1}^{J} 1\left\{Y_{it}=k\right\}\ln\Pr\left(Y_{it}=k|X_{it},\pi_{it},\beta\right)$.

Here $\beta=(\beta_2',\ldots,\beta_J')'$ and $\pi_{it}=(\pi_{it,2},\ldots,\pi_{it,J})'$ are stacked vectors. The key difference is that $\pi_{it}$ is now a $(J-1)$-dimensional vector to account for each categorical choice. Consequently, all scalar second derivative objects are replaced by their matrix analogues. In particular, $\partial_{\pi^{2}}\ell_{it}$ becomes the Hessian block $\partial_{\pi\pi'}\ell_{it}$ and cross-time products use $\partial_{\pi}\ell_{i\tau}\,\partial_{\pi'}\ell_{it}$. With this notation upgrade, the Prior SE/PE constructions carry over:

(Prior SML: multinomial logit model with strict exogenous regressors)

(Prior PML: multinomial logit model with predetermined regressors)

Probit Panels

Let $\Phi$ denote the standard normal cdf, $\Phi_{it}=\Phi(X_{it}^{\prime}\beta+\pi)$, and $H_{it}=\partial\Phi_{it}/(\Phi_{it}-\Phi_{it}^2)$. The log-likelihood is $\ell_{it}(\beta,\pi)=Y_{it}\ln\Phi_{it}+(1-Y_{it})\ln\{1-\Phi_{it}\},$ so that $\partial_{\pi}\ell_{it} = H_{it} (Y_{it}-\Phi_{it} )$, $\partial_{\pi^2}\ell_{it} = - H_{it} \partial\Phi_{it} + \frac{\partial H_{it}}{H_{it}} \partial_{\pi}\ell_{it} $, $\partial_{\pi^3}\ell_{it} = - H_{it} \partial^2\Phi_{it} - 2\partial H_{it} \partial\Phi_{it} + \partial^2H_{it} (Y_{it}-\Phi_{it} ) $, $\partial_{\beta\pi}\ell_{it} = - \partial_{\pi^2}\ell_{it} X_{it}$, and $\partial_{\beta\pi^2}\ell_{it} = - \partial_{\pi^3}\ell_{it} X_{it}$. Also, $\mathbb{E}_{\phi}(\partial_{\pi^2}\ell_{it}\,\partial_{\pi}\ell_{it})=\mathbb{E}_{\phi}(\partial H_{it} \partial \Phi_{it})$. Substituting these into (ref) gives $$ \Upsilon_{it}^{\pi} = \tfrac{\mathbb{E}_{\phi} (H_{it}\partial^2\Phi_{it} +\partial H_{it} \partial \Phi_{it} ) }{\sum_{\tau}{ \mathbb{E}_{\phi} (H_{i\tau}\partial\Phi_{i\tau}) }} + \tfrac{ \mathbb{E}_{\phi} (H_{it}\partial^2\Phi_{it} +\partial H_{it} \partial \Phi_{it} ) }{\sum_{j}{ \mathbb{E}_{\phi} (H_{jt}\partial\Phi_{jt} )}} + \tfrac{{\sum_{\tau>t} \mathbb{E}_{\phi} (H_{i\tau}\partial\Phi_{i\tau}\partial_{\pi}\ell_{it}})}{\sum_{\tau}{\mathbb{E}_{\phi}( H_{i\tau}\partial\Phi_{i\tau} )}} . $$

Carrying through the algebra shows that an explicit non–data-dependent prior that precisely cancels the bias is not available, even in the strict exogenous case. The generic priors (GE-1/GE-2) are therefore the default for probit panels. With small $N$ or small $T$, score or estimator-based corrections may yield tighter finite-sample performance.

Guidance for Practitioners: Choosing a Prior and a Correction

As a default, we recommend the robust priors (GE-1/GE-2) given in Section (ref). They remove the leading incidental parameter bias for a wide range of nonlinear panels. The only data dependence arises through terms of the form $\mathbb{E}_{\phi}(\cdot)$ and primarily reflects the cross-time dependence inherent in panel data context. In practice, whenever analytical expectations are available, they typically deliver better small-sample stability. When no closed form is available, the observed (sample-analog) form provides a convenient alternative, though it may be less stable in small samples.

With strictly exogenous regressors, one-parameter exponential family panels (logit, Poisson, and multinomial logit) work best with the non-data-dependent priors SE/SML. These priors match the relevant bias exactly and usually deliver tighter finite-sample performance. With predetermined regressors, the dedicated priors PE/PML are recommended. The remaining time-driven component vanishes as $T$ grows. When $T$ is small, a simple improvement is to estimate under SE/\textit{SML} and then apply the add-on correction $\mathcal{B}^{(\theta,1)}$ to account for the cross-time component.

For dynamic AR($m$) panels, the explicit Prior AR($m$) is best: it is non–data-dependent, fixed-$T$ consistent, and easy to compute via power-sum recursions. Prior AR($m$) itself can be used as the joint-likelihood penalty as well. In probit models, no closed-form model-specific prior is available. The safe choice for probit panel models is GE-1/GE-2. With small $N$ or $T$, one may consider score or estimator corrections.

Table (ref) lists bias-reducing priors for the integrated likelihood. For joint likelihood corrections, as discussed in Section (ref), the only additional component when moving from the joint to the integrated likelihood is the $\tfrac{1}{2}\ln\det\!\big(\overline{\mathcal{H}}(\beta,\phi)\big)$ term. The same logic delivers model-specific penalties by inserting the model’s Hessian into this term. The algebra is omitted for brevity.

Monte Carlo Experiments

This section studies the finite sample performance of our likelihood-based bias corrections. We consider four distinct nonlinear panel models: ordered logit, binary logit, multinomial logit, and probit, considering both static designs with strictly exogenous regressors and dynamic designs with predetermined regressors (including lagged dependent variables). These models span settings where model-specific priors are available (e.g., binary logit) as well as cases where only generic corrections are feasible. To our knowledge, bias correction for ordered and multinomial response panels with two-way fixed effects has not previously been studied. See honore2025dynamic for related work on dynamic ordered logit panels with individual effects.

Data-generating processes

We specify four data generating processes (DGPs) that encompass the models of interest. The design largely follows fernandez2016individual. In all cases, the unobserved individual effects $\alpha_i$ and time effects $\gamma_t$ are drawn independently from $\mathcal{N}(0, 0.25^2)$. The exogenous covariate $Z_{it}$ evolves according to an AR(1) process: $ Z_{it} = 0.5 Z_{i,t-1} + \alpha_i + \gamma_t + v_{it}, $ where $v_{it} \sim \mathcal{N}(0, 0.5),$ and $Z_{i0} \sim \mathcal{N}(0, 1)$.

For binary and ordered designs, outcomes are generated from a latent representation. In the static specification, \[ Y_{it}^* = Z_{it} \beta_Z + \alpha_i + \gamma_t + \varepsilon_{it}, \] and we set $\beta_Z = 1$. In the dynamic case we introduce state dependence via \[ Y_{it}^* = Y_{i,t-1} \beta_Y + Z_{it} \beta_Z + \alpha_i + \gamma_t + \varepsilon_{it}, \] with $\beta_Z = 1$ and $\beta_Y = 0.5$.\footnote{The multinomial logit design uses an alternative-specific index and is described separately below.} The disturbance $\varepsilon_{it}$ is standard logistic in logit-family designs and standard normal in probit. Given $Y_{it}^*$, the observed outcome $Y_{it}$ is constructed as described next:

(1) Binary response models (logit/probit). The outcome is an indicator function of the latent variable: $ Y_{it} = \mathbf{1}\{Y_{it}^* > 0\}. $

(2) Ordered logit models. The outcome takes values in $\{1, 2, 3, 4\}$. We generate the observed responses by discretizing the latent variable $Y_{it}^*$ according to three true cutoffs $\tau_1=-2.5$, $\tau_2=0.5$, and $\tau_3=2.5$: $ Y_{it} = \sum_{k=1}^{4} k \cdot \mathbf{1}\left\{ \tau_{k-1} < Y_{it}^* \leq \tau_{k} \right\}, $ with $\tau_{0} = -\infty$ and $\tau_4 = \infty$. In estimation, the cutoffs are treated as common parameters. Since the index is identified only up to location, we fix $\tau_1$ and estimate $\beta=(\beta_Z,\beta_Y,\tau_2,\tau_3)'$. Note that the cutoffs $\tau_2$ and $\tau_3$ are subject to incidental parameter bias, which our method corrects. In our simulations, the choice of cutoffs ensures that the ordered-response likelihood is well behaved under two-way effects, in the sense of the identification condition in Assumption (ref).\footnote{ For the ordered logit, Assumption (ref)(v) rules out asymptotically individuals or time periods with outcomes concentrated in a single category, and implies local identification (a nonsingular Hessian in a neighborhood) with probability approaching one as $N,T\to\infty$.}

(3) Multinomial logit models. We consider a setting with $J=3$ alternatives, where $Y_{it} \in \{1, 2, 3\}$. The probability of observing outcome $k$ is given by the standard multinomial logit formulation with base category 1: \[ \Pr(Y_{it} = k \mid \cdot) = \exp(V_{it,k}) / \big( 1 + \textstyle\sum_{j=2}^{3} \exp(V_{it,j}) \big), \quad \text{for } k \in \{2, 3\}, \] where the linear index $V_{it,k} = X_{it}'\beta_k + \alpha_{i,k} + \gamma_{t,k}$. Here, $\beta_k$ represents the alternative specific coefficients, and $\alpha_{i,k}, \gamma_{t,k}$ are the alternative specific fixed effects. In the static case, $X_{it} = Z_{it}$; in the dynamic case, $X_{it} = (\mathbf{1}\{Y_{i,t-1}=1\}, Z_{it})'$.

We set $N=45$ and consider $T\in\{15,45\}$,\footnote{Whenever a correction requires estimating cross-time spectral expectations, we set $L=1$ for $T=15$ and $L=2$ for $T=45$, in line with the small-lag choices commonly used in practice.} with 500 replications. Across designs we compare: (i) the uncorrected fixed-effects MLE; (ii) the analytical and jackknife corrections of fernandez2016individual (FW16) when available; and (iii) our bias-reducing estimators, implemented either as posterior mean under the corresponding prior or as a penalized likelihood estimators. We report Monte Carlo bias, standard deviation (SD), the ratio of average standard errors to SD, and empirical 95% coverage. APEs are computed using the plug-in estimator described in Section (ref).\footnote{ All simulation designs and estimators reported in this section can be reproduced using our open-source Python package twowaypanel. In particular, the package implements the four model classes considered here under both static and dynamic specifications, and provides likelihood-based prior/penalty corrections together with inference for parameters and APEs. For binary logit and probit panels, it also includes the analytical bias correction of FW16. }

Monte Carlo results

Ordered logit. We begin with ordered logit panels. Table (ref) reports the results for both static and dynamic specifications. The uncorrected estimator shows severe bias not only in the slope coefficients ($\beta$) but also in the cutoff parameters ($\tau$). Our likelihood-based corrections (GE-1) reduce these biases sharply and restore coverage close to the nominal level, even in the dynamic specification. For average partial effects, the static model shows only modest bias to begin with, whereas the dynamic model displays large distortions. In both cases the corrections deliver well calibrated inference. A key advantage is that the procedure treats $\beta$ and the ancillary cutoffs $\tau$ in a unified way, avoiding derivation of specific bias formulas for each parameter.

Binary logit. Table (ref) reports results for the binary logit model, a one-parameter exponential-family specification. For the static model, we use the model-specific SE correction, which is non-data-dependent and cancels the leading bias in the integrated likelihood. Consistent with this property, prior/penalty SE is nearly unbiased already at $T=15$, with no visible variance inflation. For the dynamic logit, we use the dedicated PE correction, which delivers substantial bias reduction with minimal implementation burden. At $T=15$, prior PE cuts the large bias in $\beta_{Y_{t-1}}$ from $-64.1\%$ to $-13.0\%$. With short $T$, the remaining bias is largely cross-time. Accordingly, the add-on adjustment (SE+AC) further improves performance (to $-5.6\%$ and $-3.7\%$). By $T=45$, residual bias under \textit{PE} is small (about $-4.2\%$), and \textit{SE+AC} remains close to the analytical benchmark. The same ordering holds for APEs. All likelihood-based corrections improve bias and coverage, with \textit{SE+AC} more helpful when $T$ is small.

Multinomial logit. We next consider the multinomial logit model (Table (ref)). This model extends the one-parameter exponential family to a vector-valued setting, increasing the number of nuisance parameters. Despite this added complexity, the matrix versions Prior SML (static) and Prior PML (dynamic) substantially reduce the bias in the coefficients and yield well-calibrated inference. In the dynamic case with $T=15$, the uncorrected estimator for $\beta_Y$ exhibits a large bias of $-60.6\%$. The Prior PML reduces this to $-8.3\%$, while maintaining a standard error close to the uncorrected estimator. This demonstrates the robustness of our approach in handling vector-valued outcomes and higher-dimensional fixed effects parameters.

Binary Probit. In probit panels (Table (ref)), GE-1 reduces bias relative to the uncorrected estimator and improves coverage. While the bias reduction is significant, it is slightly less complete than in the logit cases for small $T$. For example, in the dynamic probit model with $T=15$, the bias in $\beta_{Y_{t-1}}$ is reduced from $-39.8\%$ (Uncorrected) to $-15.3\%$ (Prior GE-1), compared to $-4.6\%$ for FW16. However, as $T$ increases to 45, the performance gap closes, and the Prior GE-1 estimator becomes nearly unbiased. This pattern is consistent with our theoretical predictions of the generic correction for binary probit panel, which relies on large-$N,T$ asymptotics.

Summary. Across all four model classes, the proposed likelihood-based corrections provide a simple and reusable route to debiasing two-way nonlinear panels, including ordered and multinomial response models where alternative corrections were previously unavailable. For standard logit models, the model-specific priors deliver excellent finite-sample performance, while for probit models, the generic priors provide a robust improvement that converges rapidly as the time dimension grows.

Empirical Illustration: Female Labor Force Participation

To illustrate the utility of our method, we revisit the classic relationship between fertility and female labor force participation. This relationship is central to labor economics but is difficult to identify due to the joint determination of fertility decisions and labor supply (see, e.g., angrist1998children). In contrast to earlier one-way applications fernandez2009fixed, we incorporate both individual heterogeneity (e.g., preferences for work) and aggregate time shocks (e.g., business cycles) via two-way fixed effects.

The data come from the Panel Study of Income Dynamics (PSID), covering calendar years 1979--1988 (waves 13--22). Our estimation window is 1980--1988 (9 waves) and the 1979 participation outcome is retained only to initialize the lagged dependent variable in the dynamic specification. The sample is restricted to women aged 18--60 in 1985 who remained married to a husband in the labor force throughout the sample period. The full sample contains 1,461 women observed for 9 years. Because the logit model with two-way additive fixed effects is not identified from individuals whose participation status never changes over the panel, estimation of the index parameters uses the subsample of 664 “switchers” (women who change participation status at least once).\footnote{In switcher subsample, about 16% are black women; average age in 1985 is 35.6; average schooling is 12 years; the mean participation rate is 57%. Average numbers of children per woman are 0.28 (ages 0--2), 0.36 (ages 3--5), and 1.11 (ages 6--17).}

We consider static and dynamic two-way fixed-effects logit models. Let $Y_{it}$ denote labor force participation and write $\beta=(\beta_Y,\beta_Z')'$ and $X_{it}=(Y_{i,t-1},Z_{it}')'$, where $Z_{it}$ includes the numbers of children aged 0--2, 3--5, and 6--17, and log husband's earnings. We estimate $$ \Pr(Y_{it}=1\mid X_{it},\alpha_i,\gamma_t) = F\!\left(X_{it}'\beta+\alpha_i+\gamma_t\right), \qquad F(v)={e^{v}}/{(1+e^{v})}, $$ using the $664$ switchers ($i=1,\ldots,664$) observed over $t=1,\ldots,9$.

Although index parameters are estimated on switchers, we report population-scaled averages over the full panel of $1461\times 9$ person-years. In our data, 121 women never participate and 676 always participate over 1980--1988. For these 797 “stayers”, their fixed effects index $\alpha_i$ diverges. As a result, the logit derivative $F'(\cdot)$ vanishes and so do the partial effects.\footnote{ For a continuous $X_{it,k}$, the per-observation partial effect is $ \Delta_{it,k}(\beta,\phi)=\beta_k\,F'\!\left(X_{it}'\beta+\alpha_i+\gamma_t\right), $ and for a binary regressor $X_{it,k}$, it is defined by $ \Delta_{it,k}(\beta,\phi) = F\!\left(X_{it,-k}'\beta_{-k}+\beta_k+\alpha_i+\gamma_t\right) - F\!\left(X_{it,-k}'\beta_{-k}+\alpha_i+\gamma_t\right). $ When $|\alpha_i|\to\infty$, $\Delta_{it,k}(\beta,\phi)\to 0$. } Accordingly, for APEs it is without loss to average over switchers and scale by the population size $$ \Delta(\beta,\phi):=\tfrac{1}{1461\times 9}\textstyle\sum_{i=1}^{664}\sum_{t=1}^{9} \Delta(X_{it},\beta,\alpha_i+\gamma_t). $$

In the dynamic specification, we follow browning2014dynamic and define one-step transition probabilities in a two-state first order Markov representation. Specifically, conditional on $(Z_{it},\alpha_i+\gamma_t)$:

The marginal dynamic effect is the difference in transition probabilities, $$ \Delta_{M}(Z_{it},\beta,\alpha_i+\gamma_t) :=P_{11}(Z_{it},\beta,\alpha_i+\gamma_t)-P_{01}(Z_{it},\beta,\alpha_i+\gamma_t), $$ which coincides exactly with the partial effect $\Delta(X_{it},\beta,\alpha_i+\gamma_t)$ for the discrete change in the lagged indicator $Y_{i,t-1}$. Hence the reported APE for $Participation_{t-1}$ in Table (ref) can be read as the average marginal dynamic effect (AMDE): $$ \Delta_M(\beta,\phi)=\tfrac{1}{1461\times 9}\textstyle\sum_{i=1}^{664}\sum_{t=1}^{9} \Delta_{M}(Z_{it},\beta,\alpha_i+\gamma_t). $$

browning2014dynamic also report the long-run proportion of unit values implied by a time-homogeneous two state Markov chain. In our two-way fixed effects setting, $\gamma_t$ is unrestricted and the implied transitions need not be time invariant. We therefore evaluate the long-run proportion locally by freezing $(Z_{it},\alpha_i+\gamma_t)$ at each $(i,t)$ and applying their stationary formula: \[ \Delta_{L}(Z_{it},\alpha_i+\gamma_t) := \frac{P_{01}(Z_{it},\alpha_i+\gamma_t)}{1-P_{11}(Z_{it},\alpha_i+\gamma_t)+P_{01}(Z_{it},\alpha_i+\gamma_t)}. \] We summarize this local long-run participation probability by the same population-scaled normalization used for APEs, \[ \Delta_L(\beta,\phi) =\tfrac{1}{1461\times 9}\textstyle\sum_{i=1}^{664}\sum_{t=1}^{9}\Delta_{L,it}(Z_{it},\alpha_i+\gamma_t). \]

Table (ref) reports results for four estimators: (i) the uncorrected fixed effects MLE; (ii) the analytical correction of fernandez2016individual (FW16); (iii) the bias-reducing penalized likelihood estimator (Penalty); and (iv) our posterior-mean estimator (Prior).\footnote{ For the cross-time spectral expectations entering the dynamic corrections we use $L=1$ given the short panel length ($T=9$). } We employ the closed form corrections (SE) for the static specification and (PE) for the dynamic case. For the prior-based approach, we report the posterior plug-in estimates for the APEs. Standard errors are derived from analytical asymptotic formulae for the likelihood based estimators.

The posterior quantities were computed using the algorithm detailed in Appendix (ref).\footnote{As an empirical demonstration, we deliberately run a long MCMC chain so that the reported estimates are numerically stable in this high-dimensional two-way panel setting. Specifically, we run 150,000 iterations, discard the first 50,000 as burn-in, and thin every 20th draw thereafter. The reported estimates are not sensitive to shorter runs.} Convergence diagnostics indicate excellent performance. The geweke1992evaluating test fails to reject the null hypothesis of stationarity for the key dynamic parameters (Table (ref)), and Figure (ref) shows stable traces and fast-decaying autocorrelation, confirming good mixing.

Turning to Table (ref), the uncorrected estimates in both specifications exhibit the inflated bias predicted by the theory. All three corrected estimators move the coefficients in the same direction, which is consistent with bias reduction. In the static model, the estimates from FW16, Penalty SE, and Prior SE are very similar, and the standard errors are almost unchanged. This provides a useful benchmark. In a setting where an analytical correction is already available, the likelihood-based implementations lead to essentially the same empirical conclusions.

The dynamic model is also the case in which the likelihood-based approach goes beyond the existing analytical correction. Using Prior PE and its penalized analogue, the largest changes arise for the state-dependence parameter and its associated APE. The APE for Participation$_{t-1}$, equivalently the average marginal dynamic effect $\Delta_M$ in browning2014dynamic, increases under all bias corrections. This points to stronger persistence once incidental-parameter bias is taken into account.

The likelihood-based correction extends directly to nonlinear dynamic functionals. In particular, Table (ref) reports the local long-run participation probability $\Delta_L$. We do not report an FW16-corrected counterpart because their analytical correction is derived for the index parameters and APEs, whereas $\Delta_L$ is a nonlinear Markov functional and would require a separate derivation. In our framework, once the corrected model parameters are available, $\Delta_L$ follows directly by plug-in together with the same variance-type adjustment. The similarity of the penalty and prior estimates for $\Delta_L$ is also reassuring for the substantive conclusion.

Concluding Remarks

This paper develops a practical integrated likelihood approach to bias reduction in nonlinear panel models with additive two-way fixed effects under large-$N,T$ asymptotics. Because maximizing the integrated likelihood is infeasible in this setting, we implement the correction through posterior means under a carefully constructed prior, while maintaining an entirely frequentist interpretation.

Our main theoretical contribution is to make integrated-likelihood debiasing operational in the two-way regime, where classical Laplace approximations can fail. A target-centered full-exponential Laplace expansion, together with a cumulant-based argument, yields valid asymptotic approximations for posterior moments. This delivers robust priors that are bias-reducing for both common parameters and fixed effects without variance inflation, and supports valid inference for APE after a simple variance-type adjustment. For empirical practice, we derive model-specific priors and companion joint-likelihood penalties for workhorse specifications and provide implementation guidance.

Our approach opens several avenues for future research. The theoretical machinery developed here, particularly the handling of slow-converging nuisance parameters via sparse cumulant expansions, could be adapted to models with interactive fixed effects or network dependencies. Extending the framework to unbalanced panels or to semiparametric settings, and developing corresponding bias-reducing priors, would further expand the set of tools available to empirical researchers.

Figures & Tables

{\fontsize{10pt}{11pt}\selectfont

}

{\fontsize{10pt}{13pt}\selectfont

}

{\fontsize{9.5pt}{10.0pt}\selectfont

}