Approximate Operator Inversion for Average Effects in Nonlinear Panel Models

Approximate Operator Inversion for Average Effects in Nonlinear Panel Models

\thispagestyle{empty} \setcounter{page}{1}

{\bf Keywords:} { Panel data, discrete choice, average effects, incidental parameters, ill-posed inverse problem}

{\bf JEL classification code:} {C14, C23, C25}

Introduction

Nonlinear panel data models with fixed effects are central to empirical work in economics. In the standard semiparametric setting, the researcher specifies for each individual unit $i=1,\ldots,n$, a conditional probability (or density) $f(Y_i|X_i,A_i;\theta_0)$ of outcomes $Y_i = (Y_{i1}, \ldots, Y_{iT})$ given the observed covariates $X_i = (X_{i1}, \ldots, X_{iT})$, the fixed effects $A_i$ that capture unobserved heterogeneity, and the model parameter $\theta_0$. The conditional distribution of fixed effects $\pi_0(A_i|X_i)$, on the other hand, is left completely unspecified. Two objects of typical interest are the common model parameters $\theta_0$ and average effects of the form $\mu_0 = \mathbb{E}[\mu(X_i, A_i, \theta_0)]$, such as average marginal effects or average treatment effects. The presence of unobserved fixed effects leads to the well-known incidental parameter problem neyman1948consistent, which complicates estimation of both $\theta_0$ and $\mu_0$. The literature frames the study of this problem along two asymptotic regimes, each with its own characteristic issues.

Clearly, the large-$T$ literature operates on the level of $A_i$ and not $\pi_0(A_i|X_i)$, sidestepping the fundamental inversion problem in (ref). This poses an inherent limitation. Indeed, the large-$T$ literature almost exclusively focuses on the correction of the leading $O(1/T)$ bias. While this is entirely justified for truly large values of $T$, in practice the remaining bias terms can potentially matter for moderate values of $T$, especially in more complicated models. At a practical level, complete correction of bias remains a virtually hopeless task. At a theoretical level, the focus on $A_i$ rather than $\pi_0(A_i|X_i)$ limits the scope of traditional bias-correction methods.

Motivated by these observations, we propose the novel approximate operator inversion (AOI) method. This method can be used to conduct inference on both $\theta_0$ and $\mu_0$ under large-$n$ large-$T$ asymptotics, though our focus in this paper will be on $\mu_0$. The AOI estimator operates directly on the level of $\pi_0(A_i|X_i)$ rather than of $A_i$. Specifically, our approach treats the estimation of $\mu_0$ as an operator inversion problem and it approximately inverts the likelihood-induced mapping from the fixed-effect distribution to the outcome distribution given in (ref). By shifting the focus to the fundamental inversion problem, our approach bridges large-$T$ analysis with the fixed-$T$ literature.

Crucially, although distinct from traditional bias-correction methods, the AOI estimator admits an iterated bias-correction interpretation and, in a well-defined sense, delivers infinite-order bias correction. This advances the state of the art in the large-$T$ panel literature, where existing methods typically correct bias only up to a fixed, often first, order.

The bias of AOI depends on two approximation components: how well the average-effect function $\mu(X, \cdot, \theta_0)$ is approximated by a chosen basis, and how well the unknown fixed-effect distribution is approximated by another basis. The overall bias is governed by the product of these two approximation errors, a property we call rate double robustness. Rapid convergence in either component is enough for fast overall bias decay. Under regularity conditions, the bias decays at an exponential rate in $T$, i.e., $O(\rho^T)$ for some $0 < \rho < 1$. This exponential decay property is a major improvement over the polynomial rates $O(1/T^s)$ obtained in standard bias-correction methods.

Operating on $A_i$ rather than $\pi_0(A_i|X_i)$ can result in more than the loss of infinite-order bias-correction. We highlight this by studying a third asymptotic scenario that interpolates between the large-$T$ and fixed-$T$ settings.

The traditional large-$T$ asymptotics implicitly assumes that all fixed effects are consistently estimable as long as $T\to\infty$. However, and crucially, consistent estimation of $A_i$ is not a matter of whether $T$ tends to infinity but of whether the Fisher information in data on $A_i$ grows without bound as $T\to\infty$. Depending on the model and the covariate design, $A_i$ can fail to be consistently estimable, rendering existing bias-correction approaches non-operational. Our study of this scenario highlights that the extant bias-correction literature operates not under large-$T$ asymptotics but a subset of it where $A_i$ is consistently estimable. We believe that our paper is the first one to make this distinction. Importantly, under regularity conditions, AOI is inherently immune to this issue.

The ideas in this paper apply to inference on both $\theta_0$ and $\mu_0$. From Section (ref) onwards, however, we focus on estimation of $\mu_0$, which presents the greater challenge. Even in models where $\theta_0$ can be estimated consistently for fixed $T$ (for example, using conditional likelihood methods), point identification of $\mu_0$ typically fails for fixed $T$ because the distribution of $A_i$ is only partially identified from short panels. This motivates our focus on $\mu_0$ and the development of the AOI estimator for this purpose.

\paragraph{Related literature.} AOI builds on two prior contributions. The first is functional differencing bonhomme2012functional, which constructs moment conditions that identify $\theta_0$ while remaining exactly free of $A_i$, without any large-$T$ approximation. Such exact moment conditions exist only in special models and for specific choices of $\mu_0$, and do not extend to average effects in general dano2023transition, aguirregabiria2024identification. bonhomme2017panel adopt a similar inverse-problem perspective on average-effect estimation, but restrict attention to settings where $\mu_0$ is point-identified at fixed $T$, ruling out the discrete choice models considered here.

The second is approximate functional differencing dhaene2023approximate, which replaces exact moment conditions with approximate ones whose error vanishes as $T \to \infty$, broadening the scope to settings where exact moment conditions are unavailable. Both functional differencing and approximate functional differencing are general methods applicable to a broad class of nonlinear panel models. AOI adapts and extends this approach to the estimation of $\mu_0$, going beyond dhaene2023approximate by developing the operator inversion perspective, establishing formal results on bias decay and asymptotic normality, and providing feasible inference procedures.

The bounds literature offers a complementary approach, deriving sharp partial identification regions under minimal assumptions honore2006bounds, chernozhukov2013average, davezies2021identification, dobronyi2021identification, pakel2023bounds, botosaru2024adversarial. In modern applications where $n$ is large and $T$ is only moderately large, these bounds can be very narrow, so that $\mu_0$ is nearly point-identified in practice. AOI and bounds approaches are therefore complementary. Bounds provide robust partial identification guarantees, while AOI provides point estimates and inference under large-$T$ approximations.

More broadly, AOI contributes to the large literature on bias-correction in nonlinear panel models hahn2004jackknife, arellanobonhomme2009intlike, dhaene2015split, higgins2024bootstrap, bonhomme2024neyman. These methods rely on a preliminary consistent estimator of $A_i$. AOI does not require consistent estimation of $A_i$, making it applicable in a broader class of settings and delivering better finite-sample performance when $T$ is small.

\paragraph{Roadmap.} The paper is organized as follows. Section (ref) introduces the model, defines the target average effect, and presents a motivating example. Section (ref) formulates the estimation of average effects as an inversion problem and shows how AOI approximately solves it. Section (ref) defines the AOI estimator and interprets it as an infinitely iterated bias-correction. Section (ref) establishes the large-sample properties of the estimator, including the rate double robustness result, asymptotic normality, and feasible inference. Section (ref) presents numerical results for a random-coefficient logit model. Section (ref) discusses inference on the common parameter $\theta_0$. Section (ref) concludes. All proofs are collected in the Appendix.

Setup and motivating example

This section introduces the formal framework and a motivating example. Section (ref) defines the model and the target average effect. Section (ref) illustrates that average effects can be consistently estimated even when the fixed effects cannot.

Model, average effects, and examples

We observe outcomes $Y_i \in \mathcal{Y}=\{y_{(k)},\ k=1,\dots,n_{\mathcal{Y}}\}$ and covariates $X_i \in {\cal X}\subset\mathbb{R}^{d_x}$ for units $i=1,\ldots,n$. We restrict attention to finite outcome sets $\mathcal{Y}$, although in principle the results can be extended to infinite outcome sets under appropriate conditions. Besides the observed variables $Y_i$ and $X_i$, we allow for latent variables $A_i \in {\cal A}\subset\mathbb{R}^{d_a}$. In what follows, we often drop the unit index $i$. For example, instead of $Y_i$, $X_i$, and $A_i$ we simply write $Y$, $X$, and $A$.

We assume that $(Y_i, X_i, A_i)$, $i=1,\ldots,n$, are independent and identically distributed random vectors. The model specifies the conditional outcome probabilities\footnote{In many panel models of interest, the conditional outcome probabilities also depend on a finite-dimensional common parameter $\theta_0$, so that $\mathbb{P}(Y=y \,|\, X=x, A=\alpha) = f(y \,|\, x, \alpha, \theta_0)$ and $\mathbb{P}(Y=y \,|\, X=x) = \int_{\mathcal{A}} f(y \,|\, x, \alpha, \theta_0) \, \pi_0(\alpha \,|\, x) \, \mathrm{d}\alpha$. When a {$\sqrt{nT}$}-consistent estimator of $\theta_0$ is available, $\theta_0$ plays essentially no role for inference on $\mu_0$. To avoid carrying $\theta_0$ in the notation throughout the paper, we absorb it into the model and write simply $f(y \,|\, x, \alpha)$. Inference on $\theta_0$ itself is discussed in Section (ref).}

where the function $f\left(y \, \big| \, x, \alpha \right)$ is known.

Let $\pi_0(\alpha \,| \,x)$ denote the true conditional probability density function of $A$ given $X=x$. Integrating out $A$ from (ref), the conditional outcome probabilities identified from the data are

We impose no restrictions on $\pi_0(\alpha \,| \,x)$ nor on the marginal distribution of $X$, that is, we have a semiparametric model with unknown nonparametric component $\pi_0(\alpha \,| \,x)$.

In models of the form (ref), the primary objects of interest are often functionals of the unknown conditional density $\pi_0(\alpha\,| \, x)$. In particular, we consider average effects of the form

where $\mu(x,\cdot)\in L^2(\mathcal{A})$ is a known function specifying the effect of interest, and $L^2(\mathcal{A})$ denotes the space of square-integrable functions on $\mathcal{A}$ with respect to Lebesgue measure. For example, in a panel data model, the average marginal effect of the $p$-th regressor in period $t$ on the expected outcome in period $t$ sets $\mu( x,\alpha) = \frac{\partial}{\partial x_{t,p}} \sum_{y \in \mathcal{Y}} y_t f(y \,| \, x, \alpha)$, where $y=(y_1,\ldots,y_T)$ and $y_t$ is the $t$-th component of $y$.

Random-coefficient binary choice model

Consider a static binary choice panel model with $T$ periods and $\mathcal{Y} = \{0,1\}^T$. There is no common parameter $\theta_0$, and the unit-specific parameters are $A_i = (A_{i1}, A_{i2})^\top \in \mathcal{A} \subset \mathbb{R}^2$, where $A_{i1}$ is an additive effect and $A_{i2}$ is a slope. The outcome of unit $i$ in period $t$ is \[ Y_{it} = \mathbf{1}(A_{i1} + X_{it} A_{i2} + U_{it} \geq 0), \] where $X_{it} \in \mathbb{R}$ is a scalar covariate and $U_{it}$ is independently and identically distributed across $i$ and $t$ with distribution function $F$, and independent of $A_i$. The conditional outcome probabilities are $$ f(y \,|\, x, \alpha) = \prod_{t=1}^T [F(\alpha_1 + x_t \alpha_2)]^{y_t} [1 - F(\alpha_1 + x_t \alpha_2)]^{1-y_t}, \qquad y \in \mathcal{Y}, x \in \mathcal{X}, \alpha \in \mathcal{A}, $$ where $x = (x_1, \dots, x_T)^\top \in \mathcal{X} \subset \mathbb{R}^T$ and $\alpha = (\alpha_1, \alpha_2)^\top \in \mathcal{A} \subset \mathbb{R}^2$. The average effect of interest is \[ \mu_0 = \mathbb{E}[F(A_{i1} + A_{i2})], \] the average predicted probability when the covariate is set to one. If, for all $i$, $X_{it} = 1$ for some period $t$, then $F(A_{i1} + A_{i2}) = \mathbb{P}(Y_{it} = 1 \,|\, A_i, X_{it}=1)$, and $\mu_0 = \mathbb{E}[Y_{it}]$ is identified from the data without any parametric assumption on $F$ chernozhukov2013average. When $X_{it} \neq 1$ for all $t$, however, evaluating $F(A_{i1} + A_{i2})$ requires extrapolating from the observed covariate values to the target value one. The parametric form of $F$ then becomes essential, and the difficulty of the problem depends on how much information the observed covariate variation provides about this extrapolation.

Consistent estimation of $A_i = (A_{i1}, A_{i2})$ as $T \to \infty$ requires the covariates to exhibit sufficient variation across periods. A necessary and sufficient condition is

where $\bar{X}_{iT} = T^{-1} \sum_{t=1}^T X_{it}$. When this $L^2$ condition on the covariate variation fails, the Fisher information that the data provide about $A_i$ remains bounded as $T \to \infty$, and $A_i$ cannot be consistently estimated.

Consistent estimation of the average effect $\mu_0$ requires weaker assumptions on the covariate variation. {Under appropriate regularity conditions on the model primitives (bounded covariates and a compact parameter space for $A_i$),} a sufficient condition is

This is an $L^1$ condition on the covariate variation, and it is strictly weaker than (ref): condition (ref) always implies (ref), since $\sum_{t=1}^T |X_{it} - \bar{X}_{iT}| \geq \bigl(\sum_{t=1}^T (X_{it} - \bar{X}_{iT})^2\bigr)^{1/2}$, but the converse does not hold in general.

The gap between conditions (ref) and (ref) defines the territory of regime (C) from the introduction: settings in which $\mu_0$ can be consistently estimated even though $A_i$ cannot. To make these conditions concrete, consider a two-block covariate design where \[ X_{it} = 0 \text{ for } t \leq T/2, \qquad X_{it} = c_T \text{ for } t > T/2, \] with $T$ even and $c_T > 0$. Here $c_T$ is allowed to depend on $T$, reflecting a triangular array asymptotic framework in which the covariate design may change as $T$ grows. In this design, condition (ref) reduces to $Tc_T^2 \to \infty$ and condition (ref) reduces to $Tc_T \to \infty$. With a slight abuse of notation, define the within-unit sums \[ Y_{i1} = \sum_{t=1}^{T/2} Y_{it}, \qquad Y_{i2} = \sum_{t=T/2+1}^{T} Y_{it}. \] These two sums are jointly sufficient statistics for $A_i$ given $X_i$, with $Y_{i1} \sim \mathrm{Bin}(T/2,\, F(\alpha_1))$ and $Y_{i2} \sim \mathrm{Bin}(T/2,\, F(\alpha_1 + c_T\,\alpha_2))$ conditionally on $A_i = \alpha$.

\paragraph{Regime (B): $c_T$ fixed.} When $c_T = c$ is a fixed positive constant that does not depend on $T$, both conditions are satisfied and $A_i$ can be consistently estimated as $T \to \infty$. The standard approach would estimate each $A_i$ by maximum likelihood and form the plug-in estimator $\hat{\mu}_{\rm FE} = n^{-1} \sum_{i=1}^n F(\hat{A}_{i1} + \hat{A}_{i2})$. Due to the incidental parameter problem, this estimator has a bias of order $O(1/T)$ hahn2004jackknife. By contrast, the method developed in this paper achieves a bias converging to $0$ at the exponential rate $\rho^T$ for some $0 < \rho < 1$, which is much faster. This illustrates that our approach offers substantial improvements over standard methods even within regime (B).

\paragraph{Regime (C): $c_T = 1/\sqrt{T}$.} When $c_T = 1/\sqrt{T}$, we have $Tc_T^2 = 1$ and $Tc_T = \sqrt{T}$, so (ref) fails but (ref) holds. The fixed effect $A_i$ cannot be consistently estimated, whereas $\mu_0$ can. Standard bias-correction methods, which rely on a preliminary consistent estimator of $A_i$, are not applicable in this setting. Our method remains applicable because it operates at the distributional level and does not require $A_i$ to be consistently estimated.

A consistent estimator of $\mu_0$ is \[ \hat{\mu}_0 = \frac{1}{n} \sum_{i=1}^n m_T(Y_{i1}, Y_{i2}), \] where $m_T: \{0,\ldots,T/2\}^2 \to \mathbb{R}$ is a function of the two within-unit sums that can be chosen to make $\hat{\mu}_0$ consistent for $\mu_0$ as $n, T \to \infty$. For the standard logistic $F$, an explicit such choice is \[ m_T(y_1, y_2) = \sum_{m=0}^{T/2}\sum_{l=0}^{T/2} g_T(t_m, t_l)\, L_m(y_1)\, L_l(y_2), \] where $g_T(p_1, p_2) = F\bigl((1-\sqrt{T})\,\mathrm{logit}(p_1) + \sqrt{T}\,\mathrm{logit}(p_2)\bigr)$, $t_0,\ldots,t_{T/2}$ are Chebyshev nodes on $[\epsilon, 1-\epsilon]$ {for any fixed $\epsilon \in (0,1/2)$}, and $L_0,\ldots,L_{T/2}$ are the associated unbiased polynomial estimators. With this choice, the bias satisfies \[ \sup_{\alpha \in \mathcal{A}} \bigl|\mathbb{E}[m_T(Y_{i1}, Y_{i2}) \mid A_i = \alpha] - F(\alpha_1 + \alpha_2)\bigr| \leq C\, e^{-c_0 \sqrt{T}} \] for constants $C, c_0 > 0$, so the bias decays to zero faster than any polynomial rate in $T$, even though $A_i$ is not consistently estimable. Full details and the proof are given in Appendix (ref).

Average effect estimation as an inversion problem

The average effect $\mu_0 = \mathbb{E}[\mu(X,A)]$ depends on the distribution of the unobserved $A$. We only observe outcomes $Y$, whose distribution is determined by the distribution of $A$ through the model. This section explains how to estimate $\mu_0$ by approximately inverting this relationship. The covariate value $x \in \mathcal{X}$ is fixed throughout and plays no essential role. Where it aids readability, we suppress $x$ from the discussion.

The model $f(y \,|\, x, \alpha)$ tells us the probability of each outcome $y$ given a particular fixed-effect value $\alpha$. But we are not interested in $f$ as a function of individual fixed effects. What matters is the mapping that $f$ induces at the level of distributions: if the fixed effect $A$ has distribution $\pi$ on $\mathcal{A}$, then integrating $f$ against $\pi$ produces the distribution of the outcome $Y$ on $\mathcal{Y}$. Formally, for any distribution $\pi(\alpha \,|\, x)$ of $A$ given $X = x$,

This is a linear map from distributions on $\mathcal{A}$ to distributions on $\mathcal{Y}$. Its input is a density $\pi(\cdot \,|\, x)$ on $\mathcal{A}$, which is an infinite-dimensional object. Its output is a vector of $n_{\mathcal{Y}} = |\mathcal{Y}|$ outcome probabilities, which is finite-dimensional for any given $T$.

The data identify the output of (ref), that is, the outcome probabilities $\mathbb{P}(Y = y \,|\, X = x)${, while} the average effect $\mu_0$ depends on the input, the true distribution $\pi_0(\cdot \,|\, x)$. Estimating $\mu_0$ therefore {corresponds to} inverting the distributional mapping (ref): given the observed output, recover enough about the input distribution to compute $\mu_0 = \mathbb{E}[\mu(X,A)]$.

{Unfortunately,} this is not possible in general. The input is infinite-dimensional and the output is finite-dimensional, so many different input distributions produce the same output. The average effect $\mu_0$, which depends on $\pi_0$ through (ref), is therefore typically not point-identified for finite $T$.

Approximate inversion

Although the full inversion fails, it succeeds when the input distribution is restricted to a finite-dimensional class. Let $\phi_x^{(j)}(\alpha)$, $j=1,\dots,J$, be a collection of basis functions and consider the linear span \[ \Pi_{\phi_x} := \left\{\pi_x:\mathcal{A}\to\mathbb{R}\ \big|\ \pi_x(\alpha)=\sum_{j=1}^{J} \nu_x^{(j)}\,\phi_x^{(j)}(\alpha),\ \nu_x^{(j)}\in\mathbb{R}\right\}. \] If $\pi_0(\cdot\,|\,x)$ lies in $\Pi_{\phi_x}$, the dimensions match: the mapping (ref) becomes a linear system with $n_{\mathcal{Y}}$ equations and $J \leq n_{\mathcal{Y}}$ unknowns. Under the rank condition that the $n_{\mathcal{Y}}\times J$ matrix with entries $\bigl[\int_{\mathcal{A}} f(y_{(k)}\,| \, x,\alpha)\, \phi_x^{(j)}(\alpha)\,\mathrm{d}\alpha\bigr]_{k,j}$ has rank $J$, this system is invertible and the distribution $\pi_0(\cdot\,|\,x)$ can be recovered from the observed outcome probabilities.

The inversion yields an estimating function $m_T(\cdot,x):\mathcal{Y}\to\mathbb{R}$ that satisfies, for any $\pi_0(\cdot\,|\,x) \in \Pi_{\phi_x}$,

This condition holds for every value of $\alpha$, so integrating both sides over $A$ eliminates the fixed effect entirely: $\mathbb{E}[m_T(Y,x)\,|\,X=x] = \mathbb{E}[\mu(x,A)\,|\,X=x]$. The left-hand side depends only on the observed data. When (ref) holds, $\mu_0$ is identified and can be estimated without any knowledge of the distribution $\pi_0$ or of the individual $A_i$.

The restriction $\pi_0(\cdot\,|\,x)\in\Pi_{\phi_x}$ is, of course, too strong to maintain as an assumption. AOI uses $\Pi_{\phi_x}$ as an approximation device. Even when $\pi_0(\cdot\,|\,x)$ does not lie in $\Pi_{\phi_x}$, there is a unique element $\pi_x^*\in\Pi_{\phi_x}$ whose image under (ref) matches the observed outcome probabilities, that is, $\pi_x^*$ uniquely solves \[ \mathbb{P}(Y=y\,| \, X=x)=\int_{\mathcal{A}} f(y\,| \, x,\alpha)\, \pi_x^*(\alpha)\,\mathrm{d}\alpha, \qquad y\in\mathcal{Y}. \] AOI estimates $\mathbb{E}[\mu(x,A)]$ by $\mathbb{E}_{\pi_x^*}[\mu(x,A)]$, the average-effect functional evaluated at $\pi_x^*$ instead of $\pi_0(\cdot\,|\,x)$. Since $\pi_x^*$ only approximates $\pi_0(\cdot\,|\,x)$, the exact condition (ref) no longer holds. Over all distributions $\pi_0$, we can only achieve the approximation

with an error that depends on how well $\Pi_{\phi_x}$ captures $\pi_0(\cdot\,|\,x)$.

The approximation error in (ref) can be made small because the dimension of $\Pi_{\phi_x}$ can grow with $T$. This is possible because the number of distinct outcome values $n_{\mathcal{Y}} = |\mathcal{Y}|$ itself grows with $T$.\footnote{For instance, in a binary-outcome panel model with $T$ periods, $n_{\mathcal{Y}}=2^T$.} As $T$ increases, there are more outcome probabilities available to pin down the input distribution, so the approximating space $\Pi_{\phi_x}$ can be made richer while the rank condition continues to hold. Under appropriate smoothness conditions on $\mu$ and $\pi_0$, this means the approximation error in (ref) vanishes rapidly as $T$ grows. The formal rates are established in Section (ref).

A central feature of the procedure described above is that it operates entirely at the level of distributions. The distributional mapping (ref) takes as input the distribution of $A$ and produces as output the distribution of $Y$. Inverting this mapping recovers information about the distribution of $A$, which is all that is needed to compute $\mu_0 = \mathbb{E}[\mu(X,A)]$. At no point does the method construct an estimator of the fixed effect $A_i$ for any individual unit $i$. All that is needed is the aggregate distribution of outcomes, not the ability to trace outcomes back to individual fixed-effect values.

This is a fundamental difference from standard bias-correction methods, which first estimate each $A_i$ and then correct the resulting bias. Those methods require $A_i$ to be consistently estimable, which places them in regime (B) of the classification in Section (ref). Because AOI bypasses individual fixed-effect estimation, it applies equally in regimes (B) and (C). Whether or not $A_i$ can be consistently estimated is irrelevant to the procedure. This is the reason why AOI remains valid in settings where existing methods fail.

The approach is related to functional differencing bonhomme2012functional, which constructs exact moment conditions of the form (ref) that hold for any distribution $\pi_0$ dano2023transition, aguirregabiria2024identification. Exact moment conditions of this kind exist only in specific models and for specific choices of $\mu_0$. AOI works with the approximate moment conditions (ref) instead, which are broadly available, at the cost of a bias that vanishes as $T$ grows.

The AOI estimator

This section presents the AOI estimator. Section (ref) introduces the prior and posterior densities that serve as building blocks. Section (ref) defines the transition matrix $Q(x)$, with entries given by posterior predictive probabilities, and its pseudoinverse. Section (ref) presents the estimator itself and characterizes the sieve space it implicitly uses. Section (ref) shows that the estimator arises as the limit of an iterated bias correction.

Prior, posterior, and predictive densities

The AOI estimator is built from a user-chosen prior density $\pi_{\rm prior}(\alpha\,| \, x)$ for $A$ given $X=x$. This prior is not a belief about $A$ in the Bayesian sense; rather, it is a computational device that determines the function space $\Pi_{\phi_x}$ used to approximate the unknown true density $\pi_0(\cdot\,| \, x)$. As we show in Section (ref), the choice of prior pins down functions that span $\Pi_{\phi_x}$, but the AOI estimator remains consistent as $T\to\infty$ for any prior satisfying the conditions below.

We require that $\pi_{\rm prior}(\alpha\,| \, x)$ integrates to one and satisfies the following positivity condition:

The prior need not depend on $x$: choosing $\pi_{\rm prior}(\alpha\,| \, x) = \pi_{\rm prior}(\alpha)$ simplifies computation, but $x$-dependent priors are permitted. The prior $\pi_{\rm prior}$ may differ from the true $\pi_0$.

Given $\pi_{\rm prior}$, the posterior density of $A$ conditional on $Y=y$ and $X=x$ follows from Bayes' rule:

where

is the prior predictive probability of outcome $y$. We assume that the prior is chosen such that

Condition (ref) ensures that the posterior (ref) is well-defined for every possible outcome $y\in\mathcal{Y}$. It is automatically satisfied when $\pi_{\rm prior}$ satisfies (ref) and $\mathcal{A}$ has positive Lebesgue measure.

The transition matrix and its pseudoinverse

Given $x \in {\cal X}$, the posterior predictive probability of a “future" outcome $\widetilde y\in {\cal Y}$ after having observed $y \in {\cal Y}$ is

Collecting these probabilities into a matrix, let $Q(x)$ be the $n_{\cal Y} \times n_{\cal Y}$ matrix with entries $Q_{k,\ell}(x) = Q(y_{(k)} \, | \, y_{(\ell)},x)$. Thus $Q(x)$ is a transition matrix: its $(k,\ell)$-entry is the probability that an independent replicate of $Y$ equals $y_{(k)}$, given that the original observation was $y_{(\ell)}$ and that the prior $\pi_{\rm prior}$ is used to form beliefs about $A$. We have the following lemma from dhaene2023approximate.

Let $\lambda_1(x) \ge \ldots \ge \lambda_{n_{\mathcal Y}}(x)$ denote the eigenvalues of $Q(x)$, ordered in descending order, and let $U(x)$ be the $n_{\mathcal Y}\times n_{\mathcal Y}$ matrix whose columns are the corresponding right eigenvectors. Define $ \Lambda(x) := \operatorname{diag}\bigl(\lambda_k(x)\bigr)_{k=1}^{n_{\mathcal Y}}. $ By Lemma (ref), $\lambda_k(x)\in[0,1]$ for all $k$, and $$ Q(x) = U(x)\,\Lambda(x)\,U^{-1}(x). $$ Our definition of the AOI estimator involves the inverse $Q(x)^{-1}$, or, when $Q(x)$ is singular, its Drazin inverse, \[ Q(x)^D := U(x)\,\Lambda(x)^D\,U(x)^{-1}, \] where $ \Lambda(x)^D$ is obtained from $ \Lambda(x)$ by inverting the nonzero eigenvalues $\lambda_k(x)$ and setting the zero eigenvalues to zero. If $Q(x)$ is nonsingular, then $Q(x)^D=Q(x)^{-1}$. In general, when $Q(x)$ is singular, the Drazin inverse differs from the Moore-Penrose pseudoinverse, unless $Q(x)$ is symmetric, which is usually not the case in our setup.

Definition of the estimator

Suppose we are interested in estimating the average effect $\mu_0=\mathbb{E}[\mu(X,A)]$. The AOI estimator is

where the estimating function $w^{(\infty)}$ is defined as

The estimating function has a transparent structure. For each hypothetical outcome $\tilde{y}$, the integral computes the posterior mean of $\mu(x,A)$. These posterior means are then reweighted by the entries of the Drazin inverse $Q(x)^D$, which corrects for the distortion introduced by using the posterior rather than the true conditional distribution of $A$.

The estimator uses the approximate moment condition \[ \mathbb{E}\!\left[w^{(\infty)}(Y,X) -\mu(X,A)\right]\approx 0, \] where the approximation becomes exact under the conditions given below.

As formally established in Theorem (ref), the AOI estimator is exactly unbiased whenever, for all $x \in {\cal X}$, there exists $\nu(\cdot,x):\mathcal{Y}\to\mathbb{R}$ such that

In the notation of Section (ref), condition (ref) means that $\pi_0(\cdot\,| \, x)$ lies in the sieve space $\Pi_{\phi_x}$ with basis functions \[ \phi_x^{(k)}(\alpha) := \pi_{\rm prior}(\alpha\,|\,x)\, f(y_{(k)}\,| \, x,\alpha), \qquad k=1,\dots,n_{\mathcal{Y}}. \] This particular choice of sieve space -- prior-weighted likelihood functions -- has two desirable properties:

Interpretation as iterated bias correction

We now show that $w^{(\infty)}$ arises as the limit of an iterated bias correction, extending the approach of dhaene2023approximate to average effects. The starting point is the plug-in estimating function

which replaces the unknown $\pi_0(\alpha\,| \, x)$ by the posterior $\pi_{\rm post}(\alpha\,| \, y,x)$. This introduces a bias

where $ \widetilde\mu^{(0)}(x,\alpha) := \sum_{y\in \cal Y} w^{(0)}(y,x)\, f(y\,| \, x,\alpha) - \mu(x,\alpha), $ see Appendix (ref) for details. The right-hand side of (ref) has the same form as the average effect $\mu_0 = \mathbb{E}[\mu(X,A)]$, but with $\mu$ replaced by the conditional bias function $\widetilde\mu^{(0)}$. The bias of $w^{(0)}$ is therefore itself an average effect, and can be estimated by the same plug-in construction. Replacing $\pi_0$ by $\pi_{\rm post}$ in $\mathbb{E}[\widetilde\mu^{(0)}(X,A)]$ yields the estimating function \[ b^{(0)}(y,x) := \int_{\cal A} \widetilde\mu^{(0)}(x,\alpha)\, \pi_{\rm post}(\alpha\,|\,y,x)\,\mathrm{d}\alpha, \] which, by the definition of $Q(x)$ in (ref), simplifies to $b^{(0)}(\cdot,x) = Q(x)\,w^{(0)}(\cdot,x) - w^{(0)}(\cdot,x)$. Subtracting this estimated bias from $w^{(0)}$ gives the bias-corrected estimating function \[ w^{(1)}(\cdot,x) = w^{(0)}(\cdot,x) - b^{(0)}(\cdot,x) = \bigl(2\,\mathbb{I}_{n_{\cal Y}} - Q(x)\bigr)\,w^{(0)}(\cdot,x), \] where $\mathbb{I}_{n_{\cal Y}}$ denotes the $n_{\mathcal{Y}}\times n_{\mathcal{Y}}$ identity matrix. This bias correction idea can be iterated: The bias of $w^{(q)}$ takes the analogous form $\mathbb{E}[\widetilde\mu^{(q)}(X,A)]$ with conditional bias function \[ \widetilde\mu^{(q)}(x,\alpha) := \sum_{y\in \cal Y} w^{(q)}(y,x)\, f(y\,|\,x,\alpha) - \mu(x,\alpha), \] which is again an average-effect function. Estimating its expectation by the same plug-in construction produces $b^{(q)}(\cdot,x) = Q(x)\,w^{(q)}(\cdot,x) - w^{(0)}(\cdot,x)$ --- the second term is just the plug-in $\int_{\cal A}\mu(x,\alpha)\pi_{\rm post}(\alpha\,|\,y,x)\mathrm{d}\alpha = w^{(0)}(y,x)$, which does not depend on $q$. The next estimating function is $w^{(q+1)} := w^{(q)} - b^{(q)}$, and rearranging gives the simple recursion

which yields the closed-form expression

The corresponding estimator is $$\widehat \mu^{(q)} := \frac 1 n \sum_{i=1}^n w^{(q)} ( Y_i,X_i).$$ The partial sum $\sum_{r=0}^q [\mathbb{I}_{n_{\cal Y}} - Q(x)]^r$ is a truncated Neumann series for $Q(x)^D$. The following lemma confirms that this series converges (in the appropriate sense) to $Q(x)^D$ as $q\to\infty$, so that $w^{(q)}$ converges to $w^{(\infty)}$.

The convergence holds because the eigenvalues of $\mathbb{I}_{n_{\cal Y}}-Q(x)$ lie in $[0,1)$ for all nonzero eigenvalues of $Q(x)$, so the geometric series converges on the range of $Q(x)$. Eigenvalues equal to zero contribute divergent terms, but these lie in the kernel of $Q(x)$ and are annihilated when multiplied by the posterior means.

Asymptotic theory

This section establishes the large-sample properties of the AOI estimator. We study the estimator $\widehat{\mu}^{(\infty)}$ defined in (ref)--(ref). The analysis proceeds in two steps. Section (ref) characterizes the bias $\mu_*^{(\infty)}-\mu_0$, where $\mu_*^{(\infty)}:=\mathbb{E}[w^{(\infty)}(Y,X)]$, and establishes rates of convergence. Section (ref) provides asymptotic normality and feasible inference under asymptotics where $n,T\to\infty$ jointly.

Bias

Characterization of the bias

The bias of $\widehat{\mu}^{(\infty)}$ admits a clean decomposition in terms of two approximation errors. Fix $x \in\mathcal{X}$ and define the following two subspaces of $L^2(\mathcal{A})$.

\paragraph{Likelihood span.} Let $\mathcal{F}(x)$ denote the subspace of $L^2(\mathcal{A})$ spanned by the likelihood functions $\{f(y\,| \, x,\cdot):y\in\mathcal{Y}\}$, and write \[ \mu(x,\cdot)=\mu_{\mathcal{F}(x)}(x,\cdot) + \mu_{\mathcal{F}(x)^\perp}(x,\cdot), \] where $\mu_{\mathcal{F}(x)}$ and $\mu_{\mathcal{F}(x)^\perp}$ denote the orthogonal projections of $\mu(x,\cdot)$ onto $\mathcal{F}(x)$ and its complement, respectively. The residual $\mu_{\mathcal{F}(x)^\perp}$ measures how well the likelihood basis approximates the average-effect functional.

\paragraph{Prior-weighted likelihood span.} Let $\mathcal{P}(x)$ denote the subspace of $L^2(\mathcal{A})$ spanned by the functions $\phi_x^{(k)}(\alpha)=\pi_{\rm prior}(\alpha\,|\,x) f(y_{(k)}\,| \, x,\alpha)$, $k=1,\dots,n_{\mathcal{Y}}$, and write \[ \pi_0(\cdot\,| \, x)=\pi_{\mathcal{P}(x)}(\cdot\,| \, x) + \pi_{\mathcal{P}(x)^\perp}(\cdot\,| \, x), \] where $\pi_{\mathcal{P}(x)}$ and $\pi_{\mathcal{P}(x)^\perp}$ are the projections of $\pi_0(\cdot\,| \, x)$ onto $\mathcal{P}(x)$ and its complement. The residual $\pi_{\mathcal{P}(x)^\perp}$ measures how well linear combinations of the prior-weighted likelihood functions approximate the true fixed-effect distribution.

The following theorem states our main result on bias.

Discussion of the bias structure

Theorem (ref) decomposes the bias as the inner product of two approximation residuals. We discuss each component and then derive the implied rates.

\paragraph{(i) Approximation error to $\mu(X,\cdot)$.} The factor $\mu_{\mathcal{F}(X)^\perp}(X,\alpha)$ is the component of $\mu(X,\cdot)$ orthogonal to the likelihood span $\mathcal{F}(X)$. It is small when $\mu(X,\cdot)$ is well approximated by linear combinations of the likelihood functions $f(y\,| \, X,\cdot)$, $y\in\mathcal{Y}$. In particular, Theorem (ref) implies exact unbiasedness when $\mu(x,\cdot)\in\mathcal{F}(x)$ for every $x\in\cal X$: the average effect $ \mathbb{E}\left[ \mu(x,A) \right]$ is then a linear combination of the outcome probabilities $\mathbb{P}(Y=y\,|\,X=x),y\in\cal Y$, and is therefore point-identified. The residual norm $\|\mu_{\mathcal{F}(X)^\perp}\|_{L^2(\mathcal{A})}$ can be made small when $\mu(X,\cdot)$ is sufficiently smooth. Rates are discussed in Section (ref).

\paragraph{(ii) Approximation error to $\pi_0(\cdot\,| \, X)$.} The factor $\pi_{\mathcal{P}(X)^\perp}(\cdot\,| \, X)$ is the component of $\pi_0(\cdot\,| \, X)$ orthogonal to the prior-weighted likelihood span $\mathcal{P}(X)$. It vanishes when $\pi_0(\cdot\,| \, x)\in\mathcal{P}(x)$ for all $x\in\cal X$. In that case, the AOI estimator exactly recovers $\pi_0$ by inverting the outcome probabilities $\mathbb{P}(Y=y\,|\,X=x)$, and the bias is zero. The residual norm $\|\pi_{\mathcal{P}(X)^\perp}\|_{L^2(\mathcal{A})}$ will converge to zero when $\pi_0(\cdot\,| \, X)$ is regular enough. Rates are discussed in Section (ref).

\paragraph{(iii) Rate double robustness.} Since the bias is the expectation of the inner product of two residuals in $L^2(\cal A)$, it is governed by residual product. By the Cauchy--Schwarz inequality,

Thus, the bias vanishes exactly when either $\mu(x,\cdot)\in\mathcal{F}(x)$ or $\pi_0(\cdot\,| \, x)\in\mathcal{P}(x)$, for all $x\in\mathcal{X}$ and is small whenever either approximation is good. The product structure is analogous to the double robustness property in semiparametric estimation funk2011doubly,chernozhukov2018double, but here it operates at the level of rates rather than point identification. The following corollary makes this explicit.

The rate $r_{\mu,T}$ depends on the smoothness of $\mu(\cdot, \cdot)$, and $r_{\pi,T}$ on the smoothness of $\pi_0(\cdot\,|\,\cdot)$. Smoothness of either one suffices for consistency; smoothness of both yields faster rates.

Rates of convergence of the approximation errors

We now derive the approximation rates $r_{\mu,T}$ and $r_{\pi,T}$. To keep the exposition concrete, we consider the following setting.

This example isolates the core approximation-theoretic challenge: for each unit $i$, only a single draw $Y_i\sim\mathrm{Bin}(T,p(A_i))$ is available to learn about $A_i$.

\paragraph{(i) Approximation rate to $\pi_0(\cdot)$.}

In Example (ref), the space $\mathcal{P}$ consists of the functions of the form \[ \pi_\nu(\alpha) = \pi_{\rm prior}(\alpha) \sum_{y=0}^T \nu(y)\binom{T}{y}a^{y}(1-a)^{T-y}, \qquad a := p(\alpha), \] where $\nu(\cdot)$ is any function in $N:=\{\nu:{\cal Y}\to \mathbb{R}\}$. The functions $\pi_\nu(\alpha)$ are prior-weighted linear combinations of Bernstein basis polynomials in $a=p(\alpha)$. Assume that ${\cal A}=[\underline{\alpha},\overline{\alpha}]$ is a compact interval, and let $[\underline{a},\overline{a}]=[p(\underline{\alpha}),p(\overline{\alpha})]$. Bounding $\|\pi_0(\cdot)-\pi_\nu(\cdot)\|_{L^2(\mathcal{A})}$ then reduces to bounding

Assuming also that $\pi_{\rm prior}(\cdot)$ is bounded away from zero on $\mathcal{A}$, (ref) is bounded if

is bounded. Since the Bernstein basis polynomials $\binom{T}{y}a^{y}(1-a)^{T-y}, y\in\cal Y$, span the space of polynomials of degree at most $T$, we can invoke standard results from polynomial approximation theory to bound (ref). The rate of convergence depends on the smoothness of the function $a\mapsto \pi_0(p^{-1}(a))/\pi_{\rm prior}(p^{-1}(a))$ on $[\underline{a},\overline{a}]$. In particular, by Theorem 8.1 in devore1993constructive, we have the following lemma.

\paragraph{(ii) Approximation rate to $\mu(\cdot)$.}

A similar approximation argument applies to $\mu(\cdot)$, yielding the following result.

\paragraph{Summary.} Combining Lemmas (ref) and (ref) with Corollary (ref) yields the following conclusion. If both mappings are analytic, the bias decays at rate $O(\rho_\pi^T\rho_\mu^T)$, that is, exponentially in $T$. If only one mapping---say, the one involving $\pi_0$---is analytic while the other has $s$ derivatives, the bias decays at rate $O(\rho_\pi^T T^{-s})$, which remains exponentially fast in $T$.

Inference

We establish asymptotic normality of $\widehat{\mu}^{(\infty)}$ under triangular array asymptotics where $n \to \infty$ and $T = T_n \to \infty$. Define $$ \sigma_T^2 := {\rm Var}\left(w^{(\infty)}(Y,X)\right). $$ We assume that $\sigma^2_T < \infty$ for all $T$; this is implicit in Assumption (ref) below.

Assumption (ref) restates the sampling assumption from Section (ref) for clarity. Assumption (ref) is a nondegeneracy condition ensuring that sampling variability of $w^{(\infty)}(Y,X)$ does not vanish as $T\to\infty$. Assumption (ref) is the Lindeberg condition adapted to the triangular array setting; it controls the tail behavior of the standardized influence function uniformly across the sequence $(n, T_n)$.

When $r_T = \rho^T$ for some $0<\rho <1$, as discussed in Section (ref), the condition $\sqrt{n} \, r_T \to 0$ becomes $\sqrt{n} \, \rho^T \to 0$ or, equivalently, $$ T \gg \frac{\log n}{2|\log \rho|}. $$ Thus, $T$ needs only grow slightly faster than logarithmically in $n$ for valid inference. This is a remarkably mild requirement, and stands in contrast with the $T \gg n^{\frac{1}{2(s+1)}}$ condition typically needed, for $s$-th order bias correction, when the bias decays at a polynomial rate $O(T^{-(s+1)})$.

For feasible inference, the variance $\sigma_T^2$ must be estimated. Define

Consistency of $\widehat{\sigma}_T^2$ in the triangular array setting requires control on higher moments.

Random-coefficient binary logit model: numerical results

We illustrate the AOI estimator numerically in the random-coefficient binary choice model of Section (ref). Section (ref) reports exact ($n=\infty$) bias and asymptotic standard deviation in the two-block design across a range of $T$ and iteration depths $q$. Section (ref) complements this with a Monte Carlo study at $n=500$ in a design with a continuous covariate, where we also assess coverage of the variance estimator.

Two-block covariate design

We set $F=\Lambda$, where $\Lambda$ denotes the standard logistic distribution function, and retain the two-block design \[ X_{it}=0 \quad \text{for } t\leq T/2, \qquad X_{it}=c_T \quad \text{for } t>T/2, \] with $T$ even. The average effect of interest is $ \mu_0 = \mathbb{E}\big[\Lambda(A_{i,1}+A_{i,2})\big]. $ We set the true distribution of $(A_{i,1},A_{i,2})$ equal to the product of two independent logistic distributions with mean $1$ and scale $\frac{1}{2}$ (so the variance is $\frac{1}{4}\pi^2/3$). The prior used by AOI is deliberately severely misspecified: we set it equal to the product of two independent Gaussian distributions with mean $0$ and variance $4\pi^2/3$. Since the two-dimensional integrals entering $Q$ are not available in closed form, both the true distribution, $\pi_0$, and the prior, $\pi_{\rm prior}$, are approximated numerically by $K=1000$ quantiles in each dimension, each with mass $1/K$. The calculations are carried out in Matlab for $n=\infty$, $T\in\{2,4,6,8,10,20,30\}$, and $q\in\{0,1,2,5,10,50,100,1000,10^6,\infty\}$, and we also report a regularized version obtained by truncating eigenvalues of $Q$ below $\lambda_{\min}=10^{-4}$. Because the estimators are linear in the outcome frequencies, the reported bias is the exact fixed-$n$ bias, and the reported asymptotic standard deviation is the exact fixed-$n$ standard deviation times $\sqrt{n}$.

To keep the presentation compact, Tables (ref)--(ref) report only $T\in\{2,6,20,30\}$ and $q\in\{0,1,2,10,1000,10^6,\infty\}$, together with the regularized AOI estimator. Here $\mathrm{AOI}(q)$ denotes the $q$-th order bias-corrected estimator, while AOI denotes the limit estimator $\mathrm{AOI}(\infty)$. We consider three scenarios of increasing difficulty to estimate $\mu_0$, whose true value, given our choice of $\pi_0$, is $0.8276$.

Scenario 1 sets $c_T=1$, so the target covariate value $x=1$ is directly observed in the second block. This is the point-identified benchmark in Section (ref). Table (ref) shows that $\mathrm{AOI}(q)$ converges rapidly to AOI as $q$ increases, and that AOI has zero bias. By $q=10$, the bias is already numerically negligible for all four reported values of $T$, while the increase in asymptotic standard deviation relative to $q=0$ remains moderate. Up to $T=20$, regularization has essentially no visible effect. At $T=30$, however, regularization reduces the reported asymptotic standard deviation of AOI while leaving the bias essentially unchanged.

Scenario 2 sets $c_T=1/2$. Here the target value $x=1$ is not observed, so estimation of $\mu_0$ requires extrapolation. Relative to $q=0$ and $q=1$, increasing $q$ reduces the bias substantially, but the asymptotic standard deviation can increase sharply. At $T=20$, the bias falls from $-0.077086$ at $q=0$ to $-4.100\times 10^{-4}$ for AOI, while the asymptotic standard deviation rises from $0.207323$ to $296.6$. At $T=30$, the same pattern is even more pronounced. Regularization leaves the bias essentially unchanged but dramatically reduces the asymptotic standard deviation of the AOI estimator.

Scenario 3 sets $c_T=1/\sqrt{T}$, which is the weak-variation design highlighted in Section (ref). This is the conceptually most interesting and challenging case: the slope $A_{i,2}$ is not consistently estimable, but the target average effect $\mu_0$ remains estimable. Table (ref) shows that the AOI sequence continues to reduce the bias toward zero, but the variance explosion is much stronger than in Scenario 2. At $T=20$, AOI has bias $1.1\times 10^{-3}$ and asymptotic standard deviation $4881$; at $T=30$, the reported asymptotic standard deviation of AOI is $1.433\times 10^{8}$. Regularization again stabilizes the computation sharply, reducing the asymptotic standard deviation at $T=30$ to $3.836$, while preserving the qualitative message that AOI still targets the average effect $\mu_0$ in this weak-variation design.

Taken together, the three tables make three points. First, when the target covariate value is observed, AOI reproduces the fixed-$T$ identified benchmark. Second, when estimation requires extrapolation, higher-order bias correction can remove most of the bias, though at a potentially large variance cost. Third, in the weak-variation design, AOI still reduces the bias toward zero, which is a key conceptual advantage of the method, but regularization becomes practically important because the high-order calculations are numerically unstable.

Continuous covariate design

We now turn to a Monte Carlo study with finite $n$ and a continuous covariate. The data-generating process is a binary choice logit model with unobserved heterogeneity and random coefficients, given by

where

and

We consider three average effects of interest. The first one is the average partial effect

which is a standard object of interest for a continuous covariate $X_{it}$. In addition to this, we also consider the average effects given by

and

These correspond to $\mathbb{E}[\Lambda(A_{i,1}+A_{i,2})]$ and $\mathbb{E}[\Lambda(A_{i,1})]$, respectively.

We focus on panels with $N=500$ and $T\in\{2,4,6\}$. This choice of $T$ corresponds to our setting of interest with moderately large panels. All results are based on $K=1000$ replications. We consider three choices of prior functions, given by

which are all misspecified relative to the true distributions of $A_{i,1}$ and $A_{i,2}$. Note that Prior 3 is misspecified as it ignores the dependence between $(A_{i,1},A_{i,2})$ and $X_i$.

We calculate the AOI estimator $w^{(q)}(y,x)$ by numerical integration over a discretised support for $(A_{i,1},A_{i,2})$. In particular, for a given choice of priors, we obtain a grid of $L=99$ points for each of $A_{i,1}$ and $A_{i,2}$, corresponding to equi-distant percentiles on the prior distributions of $A_{i,1}$ and $A_{i,2}$. This yields a grid size of $99 \times 99$. We also regularise $Q$ for numerical stability by clamping all eigenvalues smaller than $10^{-4}$ to $10^{-4}.$

The simulation results for the estimation of the average effects in (ref), (ref) and (ref) are presented in Tables (ref), (ref) and (ref), respectively. Each table presents the average bias across replications, the standard deviation of estimators across replications, the ratio of estimated standard errors to simulation standard deviations (SE/SD), and the 95% coverage rate of the AOI estimator. Several important patterns stand out.

\afterpage{

}

Regarding bias, in all cases the AOI estimator reduces the bias substantially compared to $q=0$, especially at $q=\infty$. Interestingly, the bias of the average partial effect $\mathbb{E} \left[ \partial P(Y_{it}=1|X_{it},A_{i,1},A_{i,2})/ \partial X_{it} \right]$ remains just below 0.01 even for $q=\infty$; see Table (ref). However, it is possible that this average effect is inherently difficult to estimate and requires greater $T$ than considered here. For $\mathbb{E} \left[ P(Y_{it}=1 | X_{it}=1, A_{i,1}, A_{i,2}) \right]$ bias appears to be small relative to its true value of $0.679$, even at $q=0$; see Table (ref). The average effect $\mathbb{E} \left[ P(Y_{it}=1 | X_{it}=0, A_{i,1}, A_{i,2}) \right]$, on the other hand, shows an interesting pattern: depending on the choice of priors, the bias can be quite small or significantly large at $q=0$; see Table (ref). All in all, the results for the uncorrected case of $q=0$ show that the severity of bias also depends on the average effect itself. Nevertheless, at $q=\infty$ and $T=6$ the bias becomes negligibly small in most settings.

The simulation results also reveal that estimator variance increases with $q$. This is not an unexpected reflection of the classical bias-variance trade-off. However, while the estimator standard deviation has a clearly increasing trend with $q$, we do not observe an explosive behaviour. As for the estimation of the standard deviation, across most configurations the variance estimator is quite accurate (as revealed by the SE/SD ratios).

Finally---and most importantly---in almost all cases, for $q=\infty$ (and even for many large but finite $q$ settings) the coverage rates are very close to the nominal coverage rate of 95%. The coverage rates at $q=0$, on the other hand, are in stark contrast to this result and often fall below 0.5. This confirms the validity of the asymptotic distribution as $q\to\infty$ even for very small values of $T$ in this complicated model, and strongly supports the validity of our approach in terms of inference.

Inference on common parameters

The main body of this paper focuses on estimation of $\mu_0$ in models where the conditional outcome probabilities $f(y\,|\,x,\alpha)$ do not depend on any common parameter. As noted before, many panel models of interest include a finite-dimensional common parameter $\theta_0$, so that the outcome probabilities take the form $f(y\,|\,x,\alpha,\theta_0)$. This section discusses two issues that arise in that setting: how inference on $\mu_0$ is affected by the presence of $\theta_0$, and how $\theta_0$ itself can be estimated.

Inference on $\mu_0$ in models with $\theta_0$

All results in this paper carry over directly to models with a common parameter $\theta_0$, provided a consistent estimator $\widehat{\theta}$ of $\theta_0$ is available. Given $\widehat{\theta}$, one simply evaluates the AOI estimating function at $\widehat{\theta}$, that is, $\widehat{\mu}^{(\infty)} = n^{-1} \sum_{i=1}^n w^{(\infty)}(Y_i, X_i, \widehat{\theta})$, where $w^{(\infty)}(y,x,\theta)$ is the estimating function from Section (ref) applied to the model $f(y\,|\,x,\alpha,\theta)$.

Consistent estimators of $\theta_0$ are available in a wide range of nonlinear panel models. For essentially every type of discrete outcome variable (binary, count data, ordered choice, multinomial choice), there exist model specifications that allow point identification and $\sqrt{n}$-consistent estimation of $\theta_0$ even at fixed $T$. In static models, this is typically achieved through conditional likelihood methods that exploit the existence of a sufficient statistic for $A_i$, as in exponential-family models rasch1961general, andersen1970asymptotic, chamberlain1980analysis. In dynamic models, appropriate specifications similarly allow estimation of $\theta_0$ via generalized method of moments honore2020dynamic. More generally, the functional differencing method of bonhomme2012functional provides a unifying framework for point estimation of $\theta_0$ in both static and dynamic panel models. These methods are well established in the literature and widely implemented in statistical software. In models where $\theta_0$ is not point-identified at fixed $T$, the approximate functional differencing (AFD) method of dhaene2023approximate can be used to obtain consistent estimators of $\theta_0$ whose bias vanishes rapidly as $T$ grows.

When $\theta_0$ is estimated, the estimation error in $\widehat{\theta}$ may need to be accounted for when conducting inference on $\mu_0$. If $\widehat{\theta}$ is $\sqrt{nT}$-consistent, as is the case for many conditional likelihood and AFD estimators, then $\sqrt{n}(\widehat{\theta} - \theta_0) = O_P(T^{-1/2})$, and the contribution of the estimation error in $\widehat{\theta}$ to the asymptotic distribution of $\widehat{\mu}^{(\infty)}$ vanishes as $T \to \infty$. In this case, no formal adjustment is required, and the asymptotic theory of Section (ref) applies directly. For good finite-sample performance, however, it may still be advisable to account for the estimation error in $\widehat{\theta}$.

More generally, if $\widehat{\theta}$ is $\sqrt{n}$-consistent and asymptotically linear with influence function $\psi(Y,X,\theta_0)$, and if $\theta \mapsto w^{(\infty)}(y,x,\theta)$ is sufficiently smooth, then the delta method gives $$ \sqrt{n}(\widehat{\mu}^{(\infty)} - \widetilde{\mu}^{(\infty)}) \approx \sqrt{n}\, G(\theta_0)^\top (\widehat{\theta} - \theta_0), $$ where $G(\theta_0) := \mathbb{E}[\nabla_\theta w^{(\infty)}(Y,X,\theta)]_{\theta=\theta_0}$ and $\widetilde{\mu}^{(\infty)}$ denotes the infeasible estimator evaluated at the true $\theta_0$. The adjusted asymptotic variance becomes $$ \widetilde{\sigma}_T^2 = \sigma_T^2 + 2\, G(\theta_0)^\top \mathrm{Cov}(\psi, w^{(\infty)}) + G(\theta_0)^\top \mathrm{Var}(\psi)\, G(\theta_0), $$ and the analogue of Theorem (ref) holds with $\sigma_T^2$ replaced by $\widetilde{\sigma}_T^2$.

Estimation of $\theta_0$

{The ideas developed in this paper also apply to inference on $\theta_0$ itself. For most of the paper, we have suppressed $\theta_0$ from the notation. Reinstating $\theta_0$, the AOI construction produces estimating functions $w(Y_i, X_i, \theta_0)$ that satisfy $\mu_0 \approx \mathbb{E}[w(Y_i, X_i, \theta_0)]$ for large $T$, with the bias decaying exponentially under our regularity conditions. Estimation of $\theta_0$ calls for moment functions of the same family but targeting zero rather than $\mu_0$, that is, functions $m(Y_i, X_i, \theta_0)$ such that $0 \approx \mathbb{E}[m(Y_i, X_i, \theta_0)]$, which can be combined into a GMM estimator for $\theta_0$.}

{The construction of such moment functions follows the same logic as the construction of $w$ in this paper, starting from the score of the MLE for $\theta$ and applying the iterated bias correction machinery to the resulting average-effect expression --- this construction is detailed, without asymptotic theory, in dhaene2023approximate. At $\theta = \theta_0$, the resulting moment functions are themselves average effects in the sense of the present paper, so the bias and variance results of Section (ref) apply directly to $\mathbb{E}[m(Y_i, X_i, \theta_0)]$ as $n, T \to \infty$. Two ingredients are not addressed in that paper: identification of $\theta_0$ from the moment conditions, and well-behavedness of the Jacobian $\nabla_\theta \mathbb{E}[m(Y_i, X_i, \theta)]\big|_{\theta = \theta_0}$. Once these are established, standard GMM cross-sectional asymptotics combined with the bias control developed here yields the asymptotic distribution of the GMM estimator for $\theta_0$. Controlling the bias at $\theta_0$ is the technically demanding step, and it is precisely what the AOI theory of our paper provides.}

Conclusion

This paper develops the approximate operator inversion method (AOI) for estimating average effects in nonlinear panel data models with fixed effects. The central idea is to recast the estimation problem as an inversion of the distributional mapping from the fixed-effect distribution to the outcome distribution. This mapping goes from an infinite-dimensional space to a finite-dimensional space, so it cannot be inverted exactly, but the approximation improves as $T$ grows because the outcome space becomes richer. The resulting estimator can be understood as the limit of infinitely iterated large $T$ bias corrections.

Two properties of AOI are worth highlighting. First, the bias has a product structure (rate double robustness), decaying at the product of the approximation rates for the average-effect function and the fixed-effect distribution. Under analyticity conditions, this yields exponential bias decay in $T$, so that only $T \gg \log n$ is needed for valid inference. Second, the method operates entirely at the distributional level and never estimates individual fixed effects. This makes it applicable in what we call regime (C) in the introduction, where the fixed effects cannot be consistently estimated, a setting not covered or discussed by any existing papers.

Acknowledgment

This paper benefited from the use of generative AI tools to assist with language editing and \LaTeX formatting; all output was carefully reviewed by the authors. All substantive content, results, and any remaining errors are the authors' responsibility.