Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
174,926 characters · 25 sections · 80 citation commands
Dealing with Logs and Zeros in Regression Models
\begingroup \footnote{The authors thank Ghazala Azmat, Michael Bognanno, Marine Carrasco, Sergio Correia, Laurent Davezies, Xavier D'Haultf{\oe}uille, Philippe Gagnepain, Christian Gourieroux, Koen Jochmans, James McKinnon, Mathieu Marcoux, Ioana Marinescu, Arnaud Maurel, Isabelle M\'ejean, Joris Pinkse, Alain Trognon, Michael Visser, Martin Weidner, Tom Zylkin, and all participants to seminars and conferences for insightful comments. STATA & R packages under development are available at \url{https://github.com/ldpape/iOLS-i2SLS}.} \addtocounter{footnote}{-1} \endgroup
\thispagestyle{empty}
The log-transformation of the outcome variable is widely used in linear regression models. Its main appeal is that coefficients represent the average treatment effect (ATE) in percentage terms, defined as $ ATE_{\%} = E[(Y(1)-Y(0))/Y(0)]$ for a binary treatment. Yet, this approach faces three fundamental limitations. First, the log is undefined when the outcome takes zero values, and no alternative monotone transformation consistently identifies $ ATE_{\%}$ chen2024logs. Second, in the binary regressor case, the OLS estimator of the log-linear model is biased for $ATE_{\%}$ under heterogeneous treatment effects or heteroskedasticity Kennedy1981. Third, the estimand does not, in general, coincide with the structural parameter of an exponential conditional mean model manning2001estimating.
In this paper, we develop a new estimator that targets a distinct parameter that preserves the percentage-change interpretation while addressing all three limitations. This parameter, the normalized ATE, is defined as $ ATE_{norm} = E[Y(1)-Y(0)]/E[Y(0)]$ for a binary treatment, which measures the average effect of treatment relative to the mean untreated outcome.
Our estimator, iterated OLS (iOLS), is the theoretically justified analogue of the ad-hoc $\log(1+Y)$ transformation commonly used to accommodate zeros. Whereas this popular fix imposes an arbitrary constant shift, iOLS replaces this shift with a data-driven adjustment that is updated iteratively through OLS until convergence. This procedure yields an estimator that is nearly as simple to implement as the popular fix but is consistent for the structural parameter of the exponential mean model. More generally, we show that iOLS solves the same objective function as Gamma Pseudo-Maximum Likelihood (GPML), via a sequence of fixed-point updates. Our key methodological contribution is to reinterpret the exponential regression model as a simple iterative linear regression model, thereby making available the full machinery of linear models, such as partialling out fixed effects and using instrumental variables.\footnote{Our framework naturally extends to any power variance function barlev1986 beyond GPML, i.e. the class of Tweedie distributions with a log link including PPML and the compound Poisson-Gamma. A systematic study of their properties is left for future research.}
Building on the asymptotic theory of ds2005, we establish $\sqrt{n}$-consistency and asymptotic normality of iOLS (and i2SLS) under standard assumptions. Our procedure offers several advantages over existing approaches: (i) it can be implemented by ordinary least squares with a single $X'X$ matrix inversion, making it computationally fast and directly comparable to log-linear regression; (ii) it extends naturally to endogenous regressors through iterated two-stage least squares (i2SLS); (iii) it avoids the incidental-parameter bias and the computational burden of high-dimensional fixed effects using Frisch-Waugh-Lovell projections; and (iv) it is globally convergent, with convergence guaranteed even in finite samples. In contrast, generalized linear models (GLMs) for exponential mean specifications are typically estimated by Iteratively Reweighted Least Squares (IRLS), a Hessian-based iterative method that is only locally convergent and does not allow partialling out fixed effects or incorporating instrumental variables.
To document the ubiquity of log specifications, the presence of zeros, and the lack of consensus on how to handle them, we collected evidence from three sources: (i) a review of empirical publications, (ii) a survey conducted in seminars, and (iii) an online research forum. First, we reviewed all articles published in the American Economic Review (AER) between 2016 and 2020. Figure (ref) summarizes our findings. Log-linear and log-log models rank among the most frequent specifications: nearly 40% of empirical papers employed a log specification, and 36% of these faced the problem of the log of zero. Several solutions were used in practice. Most authors retained zero observations and either (1) added a positive discretionary constant (48%), (2) employed Poisson Pseudo Maximum Likelihood (PPML) or GPML (35%), or (3) applied the inverse hyperbolic sine (15%). Discarding non-positive observations occurred in 31% of publications, and in about 20% of cases authors compared multiple methods to assess robustness.\footnote{This excludes cases where a linear specification was used as a default. See Table (ref) in the Appendix for details on data collection.} Second, we conducted a survey during three online seminars in economics departments, asking researchers: “What would you do when facing the log of zero?” The distribution of responses closely mirrored the AER evidence, except for a stronger stated preference for mixture models.\footnote{Appendix (ref) provides details about the survey design and responses.} Third, the issue extends well beyond economics. A ResearchGate thread titled “Log transformation of values that include 0 (zero) for statistical analyses?”, posted in 2014, had been viewed over 120,000 times and had attracted 38 contributions from researchers in medicine, biology, statistics, and engineering by August 2020. The suggested solutions were again broadly consistent with those observed in economics.
Taken together, the evidence shows that the log transformation and the log-of-zero problem are pervasive in empirical research and not confined to economics. Moreover, the widespread reliance on simple ad hoc fixes over theoretically grounded methods points either to numerical difficulties in implementation or to a reluctance to adopt alternative approaches.
Researchers have long attempted to adapt the log-linear model to accommodate zeros. The most common fix is to add a positive constant before taking logs. This method is simple but arbitrary: the choice of constant affects both estimates and standard errors, and the resulting bias depends on the data at hand winkelmann2008econometric,cohn2021count. Alternative monotone transformations, most prominently the inverse hyperbolic sine, have been proposed mackinnon1990transforming,burbidge_et_al_1988,ravallion_2018. These avoid discarding zeros but effectively add an observation-specific constant before logging, leading to coefficients with nonstandard or undefined interpretations bellemare_wichman_2018. While these transformations are widespread in applied work, they fail to identify $ATE_{\%}$, as shown by chen2024logs. A second strategy is to drop zero outcomes and estimate the log model on the remaining sample. This practice may however introduce a selection bias young_young_1975. The literature has occasionally turned to more structural fixes, such as mixture models or Heckman-type corrections eaton_tamura,heckman1979sample, which explicitly model the zero outcomes or address the selection bias. Yet these approaches hinge on restrictive distributional assumptions and are not widely applied in contemporary practice.
Even in the absence of zeros, important limitations arise when targeting $ATE_{\%}$. Under an exponential conditional mean specification, log-linear regression does not consistently recover the structural parameter whenever errors are heteroskedastic manning1998logged,manning2001estimating. The problem stems from Jensen's inequality implies that the conditional mean of the outcome cannot in general be recovered from its conditional log-mean, so OLS on logs delivers biased estimates of the structural parameter, as is well documented in the trade and labor literatures santos_silva_tenreyro_2006,Blackburn2007. Similar concerns apply to reduced-form interpretations: when the regressor of interest is binary, the log-linear model recovers only an approximation to $ATE_{\%}$. With homogeneous effects and homoskedastic errors, simple retransformation can deliver the correct parameter, but in the presence of heteroskedasticity more elaborate corrections are required Kennedy1981. Moreover, as we show in this paper, when treatment effects are heterogeneous the log specification yields a systematically downward-biased estimate of the true average percentage effect.
To address these concerns, a natural and theoretically justified approach is to focus on $ATE_{norm}$. As emphasized by chen2024logs, this parameter is scale-invariant and preserves the interpretation of treatment effects in percentage terms. Moreover, we show that it is intimately connected with the widely used log(Y+1) transformation through our iOLS procedure. Estimation of $ATE_{norm}$ typically proceeds via pseudo-maximum likelihood (PML) methods gourieroux_etal_1984, most notably PPML and GPML,\footnote{A more systematic comparison betweem PPML and GPML is provided in Section (ref).} both special cases of generalized linear models with an exponential conditional mean. These estimators, however, face two important limitations: an incidental-parameter bias in the presence of multi-way fixed effects, and the difficult extension to settings with endogenous regressors. These limitations may in part explain their limited adoption in empirical practice.
First, the dimensionality of the fixed effects can make estimation computationally burdensome. Second, and more fundamentally, the incidental parameter problem induces an asymptotic bias in the point estimates and their estimated standard errors whenever the number of fixed effects increases with the sample while at least one dimension of the fixed effects remains fixed. The resulting bias is typically of order $O(1/T)$, where $T$ denotes the length of the short panel dimension, and does not vanish asymptotically as $N \to \infty$ with $T$ fixed fernandez2016individual.
For PPML with two-way fixed effects, fernandez2016individual show that the coefficient estimates are still consistent in short panels, but the incidental parameter problem distorts the asymptotic distribution, leading to biased standard errors. They develop a jackknife correction that removes the leading bias, thereby restoring valid inference. For two-way gravity models, jochmans2017two develops a transformation to condition out the nuisance parameters analytically, thereby eliminating the incidental parameter bias and delivering valid inference without recourse to jackknife or analytical corrections. weidner2021bias extend fernandez2016individual to multi-way fixed effects and provide explicit bias expansions for PPML. In the canonical three-way gravity specification, where the number of countries grows and the associated exporter-time and importer-time effects increase only at the square-root rate relative to the number of observations, they show that PPML coefficients remain consistent when $N \to \infty$ and one panel dimension is fixed, but the limiting distribution is not centered at the true parameter and conventional sandwich standard errors are likewise biased. They also provide analytical corrections that substantially improve finite-sample inference, and demonstrate that PPML is unique among a broad class of PML estimators in being generally consistent under these asymptotics. More generally, when the number of fixed effects grows at the same rate as the number of observations, rather than at the square-root rate implied by trade panels, even PPML is inconsistent in multi-way settings.
Our iOLS estimator addresses both the computational burden and the incidental parameter problem in general fixed-effects settings by eliminating the need to estimate nuisance parameters directly, thereby extending the validity of GPML to arbitrary fixed-effects structures without recourse to case-specific analytical corrections.
A further limitation of PPML and GPML arises in the presence of endogenous regressors. These estimators do not easily accommodate instrumental variables in general settings. Existing approaches have relied either on nonlinear least squares implementations mullahy1997instrumental or on GMM-type estimators santos_mulhally, both of which substantially increase computational complexity, or on control-function methods wooldridge_cf, which require stronger assumptions. Endogeneity is also the context in which the differences between PPML and GPML become most fundamental. PPML rests on an additive error formulation, making the estimator “awkward to motivate directly” mullahy1997instrumental. This is particularly evident with omitted variables: the additive specification implicitly treats the omitted factor as having an additive effect, whereas all included regressors enter multiplicatively. Consequently, the implied moment conditions differ from those in the multiplicative case, so that a given set of instruments will not generally be orthogonal to both error specifications santos_mulhally. Moreover, elasticities in the additive formulation depend on both observed and omitted variables, unlike in the multiplicative case.
Recent work has sought to address these limitations in restricted settings. Extending the approach of jochmans2017two, jochmans2022instrumental propose an IV estimator for exponential regression models with two-way fixed effects. Their method adapts classical IV moment conditions to the pseudo-demeaned structure, yielding consistent and bias-corrected estimates even when regressors are correlated with unobserved heterogeneity. These advances complement earlier works, showing how endogeneity can be handled in particular designs. In contrast, our iOLS framework delivers a general solution: by exploiting the linear structure of each iteration, instrumental variables can be incorporated through iterated 2SLS, thereby extending the validity of GPML to settings with arbitrary fixed effects and endogenous regressors, while retaining the computational simplicity of linear regression.
The small-sample properties of iOLS and i2SLS are examined through numerical simulations adapted from weidner2021bias and jochmans2022instrumental in the context of gravity models. To demonstrate the broader applicability of our approach, we revisit three influential studies outside trade. The first is michalopoulos_etal_2013, who investigate the long-run consequences of pre-colonial political centralization for contemporary African development. Their analysis employs the popular transformation $\log(Y+\Delta)$ to handle zeros. The second is azoulay, who examine how the premature death of a star biomedical scientist reshapes the trajectory of her field using PPML in a three-way fixed-effects panel design. The third is inequality, who study the relationship between innovation, income inequality, and social mobility in the United States, using log-linear regressions with endogenous regressors.
The remainder of the paper is organized as follows. Section (ref) presents our methodology and asymptotic results. Numerical simulations are in Section (ref). Section (ref) contains the three applications. Section (ref) provides practical guidance to empirical researchers. Section (ref) concludes the paper.
We first present the parameter of interest. Consider an i.i.d. sample of $n$ observations $\{Y_i, X_i\}_{i=1}^n$. For now, let us assume that the variable $X_i \in \{0,1\}$ is a binary treatment and adopt the potential outcomes framework by expressing the observed outcome as $Y_i = X_i \, Y_i(1) + (1 - X_i) \, Y_i(0)$, where $Y_i(1)$ and $Y_i(0)$ denote the potential outcomes for unit $i$ under treatment and control, respectively.
The average treatment effect in percentage terms is defined as
If $Y_i>0$ for all $i=1,\ldots,n$, a common approach is to estimate this target parameter using the log-linear regression
with parameters $\alpha,\beta \in \mathbb{R}$ and an i.i.d. error term $\varepsilon_i$ satisfying the standard exogeneity assumption $E[\varepsilon_i \mid X_i] = 0$. Under these conditions, the coefficient $\beta$ measures the average treatment effect in log terms expressed as
Unfortunately, $ATE_{\%}$ suffers from three important limitations: it is not identified with zeros, biased in the binary regressor case under heterogeneity or heteroskedasticity, and generally distinct from the structural parameters of an exponential conditional mean. These limitations motivate our emphasis on the normalized ATE and the exponential mean framework.\footnote{Further discussion of these issues, along with practical guidance on how to choose between alternative target parameters and estimation methods, is provided in Section (ref).} Define
which preserves the percentage interpretation, now relative to the control mean.\footnote{In a difference-in-differences setting, the identified parameter would be the ATT relative to the average for the treatment group in the pre-treatment period.} This parameter is consistently estimable from the exponential conditional mean model
where $U_i = \exp(\varepsilon_i)$ is an i.i.d. error term satisfying $E[U_i \mid X_i] = 1$. This specification can be viewed as the exponential transformation of (ref), but under a different exogeneity restriction.\footnote{Remark that $E[U_i \mid X_i] = 1$ is not equivalent to $E[\varepsilon_i \mid X_i] = 0$ except under independence between $\varepsilon_i$ and $X_i$ or homoskedastic log-normal errors $\varepsilon_i$.} It delivers the conditional mean function $E[Y_i \mid X_i] = \exp(\alpha+X_i\beta)$, which implies $ATE_{norm} = \exp(\beta)-1$ in the binary case.
In the continuous case (i.e., $X_i \in \mathbb{R}$), $ATE_\%$ corresponds to the average semi-elasticity whereas $ATE_{norm}$ is the semi-elasticity around the average.\footnote{They are elasticity parameters if $X_i$ is log-transformed.} Formally, under $E[\varepsilon_i \mid X_i] = 0$, we have \[ ATE_{\%} = E\left[\frac{\partial Y_i/\partial X_i}{Y_i}\right] = E\left[\frac{\partial E[\log(Y_i) \mid X_i]}{\partial X_i}\right] = \beta, \] whereas under $E[U_i \mid X_i] = 1$, it follows that \[ ATE_{norm} = E\left[\frac{E[\partial Y_i/\partial X_i \mid X_i]}{E[Y_i \mid X_i]}\right] = E\left[\frac{\partial E[Y_i \mid X_i]}{\partial X_i}\frac{1}{E[Y_i \mid X_i]}\right] = \beta. \]
To estimate the normalized ATE, we consider the multiplicative representation of the exponential conditional mean model, given by
where $U_i\geq 0$ is an i.i.d. error satisfying $E[U_i|X_i]=1$, $\beta$ is a fixed parameter of interest in $\mathbb{R}^{K}$, with $K\geq 1$. Let $X$ denote the $n\times K$ matrix comprised of the $K$-dimensional column vector $X_i$ with elements $X_{ki}$, for $1\leq k \leq K$.
In contrast, many researchers focus on the alternative additive representation $Y_i = \exp(X_i'\beta) + \epsilon_i$ with $\epsilon_i = \exp(X_i'\beta)(U_i-1)$ as the error term, or use alternative transformations to approximate the log-linear regression.
Before proceeding to the general algorithm and our theoretical results, we provide a pedagogical overview of our approach. It involves three main ingredients: 1. a log-linear transformation akin to $\log(1+Y)$, 2. OLS estimates, and 3. a fixed-point equation.
The popular fix used in empirical settings is to add a small positive constant $\Delta$ to all observations $Y_i$ so that $\tilde{Y}_i = Y_i+\Delta>0$ can be log-transformed. We can reconcile this practice with the exponential conditional mean model by letting $\Delta_i$ vary across observations such that $Y_i+\Delta_i>0$. From (ref), we have
Letting $\Delta_i = \delta \exp{(X_i'\beta ) } $, for some arbitrary $\delta>0$, and subtracting $\log(1+\delta)$ on both sides yield the pseudo-linear model\footnote{Although the underlying model is not linear in $Y_i$, the transformation recovers a linear structure on the right-hand-side.}
where the error term $\overline{\upsilon}_i = \log \left( \frac{U_i + \delta}{1+\delta} \right)$ is approximately conditionally mean-independent of $X_i$ under $E[U_i|X_i] = 1$. To see this, take the Taylor expansion of $\log(U_i+\delta)$ around $U_i=1$ to obtain
Therefore, under the conditional mean restriction $E[U_i|X_i]=1$, the leading term vanishes and the residual bias is of order $O(\delta^{-2})$. Hence, the pseudo-linear representation in (ref) satisfies the exogeneity condition asymptotically.
The classical approach to estimating the transformed model (ref) would be to solve the non-linear least-squares problem
for fixed $\delta$. As $\delta \to \infty$, the minimization problem in (ref) converges uniformly to
which corresponds to the non-linear least-squares problem associated with the unconditional moment restriction $E[X_i(U_i-1)] = 0$, i.e. the GPML score equations for a log-link exponential mean model with a gamma distribution gourieroux_etal_1984.
The key advantage of transformation (ref) is that $\beta$ can be estimated using only OLS, instead of typical gradient- or Hessian-based algorithms. Indeed, the OLS estimator applied to (ref) is given by
where the vector $\tilde{Y}(\hat{\beta})$ has components $\tilde{Y}_i = \log\left(\frac{Y_{i} + \delta\,\exp(X_i'\hat{\beta})}{1+\delta}\right)$. Since $\hat{\beta}$ appears on both sides of (ref), the solution is characterized by a fixed-point equation, so OLS must simply be applied repeatedly until convergence.
To formalize this, define the mapping $T:\mathbb{R}^k \to \mathbb{R}^k$, $ T(\hat{\beta}) = (X'X)^{-1}X'\tilde{Y}(\hat{\beta})$, so we have $\hat{\beta} = T(\hat{\beta})$. It is intuitive to see that $T$ is a contraction. To that end, consider the derivative of $\tilde{Y}_i(\hat{\beta})$ with respect to $\hat{\beta}$. Differentiating yields \[ \frac{d}{d\hat{\beta}}\tilde{Y}_i(\hat{\beta}) =\frac{\delta\,\exp(X_i'\hat{\beta})}{Y_i + \delta\,\exp(X_i'\hat{\beta})}\,X_i \leq X_i, \] so the norm of the derivative of $\tilde{Y}_i(\hat{\beta})$ with respect to $\hat{\beta}$ is bounded by the norm of $X_i$. The Jacobian is therefore $ T'(\hat{\beta}) = (X'X)^{-1}X'D(\hat{\beta})X$ using the diagonal matrix $D(\hat{\beta})=\operatorname{diag}\left(\frac{\delta\exp(X_1'\hat{\beta})}{Y_1+\delta\exp(X_1'\hat{\beta})};\ldots;\frac{\delta\exp(X_n'\hat{\beta})}{Y_n+\delta\exp(X_n'\hat{\beta})}\right)$. Under a mild overlap on positive outcomes condition,\footnote{There is overlap on positive outcomes when regressors display some variability on observations with \(Y>0\), so the iteration does not place all weight on boundary cases with \(Y=0\). Intuitively, this excludes separation-type patterns, e.g. a binary regressor that equals \(1\) only when \(Y=0\).} the operator norm of $T'(\hat{\beta})$ is strictly below $1$ so the mapping $T$ is a contraction and the fixed-point iteration $\hat{\beta}^{(t+1)} = T(\hat{\beta}^{(t)}) $ converges to the unique fixed point $\hat{\beta}$ by the Banach fixed-point theorem. Importantly, this iterative algorithm is globally convergent since this operator norm is below 1 for any $\hat{\beta}$.
Therefore, for any fixed $\delta>0$ and initial value $\hat{\beta}^{(0)}$, we can repeatedly use OLS to obtain a unique fixed-point that corresponds to an approximate solution $\hat{\beta}_{\delta}$.
Finally, we can obtain a consistent estimator of $\beta$ in a computationally efficient way by solving a sequence of fixed-point equations for an increasing set of $\delta$ values. Remark that $(X'X)^{-1}$ needs only be computed once for the full sequence. For each new value of $\delta$, we use the previously computed solution as a warm start. This strategy is justified because (i) numerical convergence requires a smaller number of fixed-point iterations for smaller $\delta$, and (ii) the solution of the fixed-point equation depends continuously on $\delta$; that is, small changes in $\delta$ induce only small changes in the solution, thereby reducing the number of iterations required for convergence.
This continuity property can be seen as follows. By construction, the transformed response $\tilde{Y}_i$ is a continuous function of $(\beta;\delta)$. Let us write the fixed-point equation as \[ F(\beta,\delta)=\beta-(X'X)^{-1}X'\tilde{Y}(\beta;\delta), \] with \(F(\hat{\beta}(\delta),\delta)=0\). Total differentiation with respect to \(\delta\) yields $\frac{\partial F}{\partial \beta}\frac{d\hat{\beta}(\delta)}{d\delta}+\frac{\partial F}{\partial \delta}=0$. Since the contraction property guarantees that the Jacobian \(\frac{\partial F}{\partial \beta}\) is invertible, we have $\frac{d\hat{\beta}(\delta)}{d\delta} = -\left[\frac{\partial F}{\partial \beta}\right]^{-1}\frac{\partial F}{\partial \delta}$. All the functions involved being continuously differentiable, \(\hat{\beta}(\delta)\) is continuously differentiable with respect to \(\delta\).
In practice, our algorithm approaches the limiting case by solving the fixed-point problem for a sequence of increasing values of $\delta$. The continuity argument above implies that the sequence of solutions $\hat{\beta}(\delta)$ converges smoothly to the limiting solution as $\delta \to \infty$. Hence, setting $\delta$ arbitrarily large yields an approximation error of order $O(\delta^{-2})$. To eliminate this residual bias exactly, we use a modified fixed-point equation that enforces the unconditional moment condition associated with GPML in a final step. This ensures that the final estimator is both computationally efficient and free of approximation bias.
Let us make the following assumptions about $Y_i$, $X_i$ and the error term $U_i$.
These conditions are standard in econometrics and ensure the invertibility of $X'X$, the existence of probability limits, and a finite asymptotic covariance matrix for the estimator.
Assumption (ref) guarantees a causal interpretation of the parameters, although consistency alone only requires $E[X_i(U_i-1)]=0$.
Assumption (ref) requires that, in every regressor direction \(v\), a fixed fraction $\alpha$ of the sample variability \(v'X_i\) is contributed by observations with $Y_i>0$. This rules out boundary patterns where regressors vary only when $Y_i=0$. In particular, for a binary regressor $X\in\{0,1\}$, the condition excludes the degenerate case where $X_i=1$ only when $Y_i=0$.
In the general setting, our approach consists of transforming the response variable into
where $c_i(\beta,\delta)$ was set to $\log(1+\delta)$ in the preceding section to illustrate our approach in a simple setting. In practice, it is preferable to enforce the mean-zero condition of the transformed error $E[\overline{\upsilon}_i]=0$ in order to identify the intercept term in the pseudo-linear model:
giving the transformed error term $\overline{\upsilon}_i=\log(U_i + \delta) - c_i(\delta,\beta)$.
For any fixed $\delta<\infty$, we can simply choose $c_i(\beta,\delta) = E[\log(U_i+\delta)]$,\footnote{ In practice, one must use an empirical analog $\hat c(\delta,\beta)$ obtained as follows. Let us write $X_i'\beta = \alpha+X_{i1}'\beta_1$, with $\beta = (\alpha,\beta_1')'$ where $\alpha$ is the constant term and $X_{i1}'\beta_1$ represents the non-deterministic part. From $E[U_i]=1$ it is easy to obtain $\alpha(\beta_1) = \log(E[Y_i\exp(-X_{i1}'\beta_1)])$ then to characterize the desired term as a function of $\beta_1$: $c(\delta,\beta) = E[\log( Y_{i} + \delta\exp(\alpha(\beta_1)+X_{i1}'\beta_1) ) -\alpha(\beta_1) - X_{i1}'\beta_1]$, finally replace by the empirical counterparts.} so that under Assumption (ref), we have $E[\overline{\upsilon}_i]=0$ and $E[\overline{\upsilon}_i|X_i] = O(\delta^{-2})$. The solution of the fixed-point equation $\hat{\beta} = X'X^{-1}X' \tilde{Y}(\hat\beta,\delta)$ is thus a biased estimator of $\beta$, denoted $\beta_{\delta} = \beta + O(\delta^{-2})$.
Although the transformation (ref) is not defined for $\delta = \infty$, we can still enforce this limiting case using
leading to the alternative pseudo-linear model
where the tuning parameter $\rho>0$ governs the convergence of the fixed-point algorithm but, unlike $\delta$, does not induce an approximation bias.\footnote{Note that one could use $c_i^{\infty}(\rho,\beta) = \log(\rho+U_i)-\frac{1}{1+\rho}(U_i-1)\exp(X_i'\beta)$ to retrieve the PPML solution rather than GPML. Consistency and asymptotic normality continue to hold.} The fixed-point mapping associated with (ref) is computationally more demanding and, unlike the mapping derived from (ref), it is not globally convergent. As a result, convergence is only guaranteed when the initialization is sufficiently close to the true parameter.
To overcome this difficulty, we adopt a two-phase strategy: the first phase solves a sequence of approximate problems for increasing values of $\delta$, where the log transformation ensures global and fast convergence to increasingly accurate approximations. The second phase then applies the limiting transformation in (ref), starting from the Phase-1 solution, to eliminate the residual $O(\delta^{-2})$ bias. This design ensures both uniform global convergence and consistency.
We establish the asymptotic properties of this algorithm in the following theorem.
This theorem guarantees $\sqrt{n}$-consistent estimates and, for any fixed $n$, each iterative process converges after a finite number of iterations: $t(n)\geq - \frac{1}{2}\log(n)/\log(\kappa)$, where $\kappa\in[0,1)$ is the modulus of the associated contraction mapping. The numerical convergence will hence be slower for larger sample sizes $n$ and modulus $\kappa$ closer to 1, where $\kappa$ depends on the DGP and is affected by the choice of the sequence for $\delta$ and stability parameter $\rho$.
The algorithm is globally convergent in its first phase: for any finite $n$ and any initialization, it converges to a biased estimate that can be made arbitrarily close (as $\delta\to\infty$) to the consistent solution obtained in the second phase. Although the second phase is only locally convergent, this is not problematic because the first phase can approach the consistent solution.
In terms of computation, the first phase solves a sequence of inexpensive subproblems, converging more rapidly for smaller $\delta$ and more slowly as $\delta$ increases. The second phase requires choosing $\rho$ large enough to ensure $\kappa<1$ but not so large as to unduly slow convergence. In practice, we select $\rho$ by evaluating an increasing grid and estimating $\kappa$ from the first few iterations as the median of $\hat{\kappa}_{t+1} = |\beta_{t+1}-\beta_{t}|/|\beta_{t}-\beta_{t-1}|$. If $\hat{\kappa}>1$, we increment $\rho$ until the contraction condition $\hat{\kappa}<1$ is satisfied.
This theorem also provides the asymptotic distribution of our estimated parameters. It corresponds to the distribution of OLS estimates in the last iteration (once the estimator has converged) after a simple reweighting of the corresponding covariance matrix. We can therefore use any robust covariance estimator.
This approach extends naturally to the endogenous regressor setting. Let us assume $E[U_i-1|X_i] \neq 0$ and define $Z$ as a $n\times L$ matrix with $L\geq K$ instrumental variables. We denote $P_z$ as the projection matrix $Z(Z'Z)^{-1}Z'$.
We make the following standard assumptions about the rank conditions, finiteness of higher-order moments, relevance of instruments, compactness of the parameter space, and exogeneity conditions of the instrumental variables.
Assumption (ref) requires that, in every instrumented regressor direction \(v\), a fixed fraction $\alpha$ of the sample variability is contributed by observations with $Y_i>0$. Again, this rules out separation patterns where the instrumented part of $X$ varies only when $Y_i=0$.
We establish the asymptotic properties of this estimator in the next theorem.
This result shows that the advantages of 2SLS carry over to our iterative setting.
Nonlinear regressions with multi-way fixed effects (FE) face two related challenges: a high-dimensional set of nuisance parameters and an IPP when some FE dimensions are short. We tackle both using Frisch-Waugh-Lovell projections to partial out FE indicators. We develop three complementary algorithms tailored to the empirical design: (i) when every FE dimension is large, iOLS-HDFE handles the high dimensionality but would suffer an IPP if any FE dimension were short; (ii) when FE dimensions are short, iOLS-FD removes them ex ante and iterates entirely in the residualized space, thereby avoiding the IPP; and (iii) with a mix of large and short FE dimensions, iOLS-FD-HDFE estimates the large FE dimensions and partials out the short ones.
First, let us partition the regressor matrix as $X = [X_{0}, X_{1}]$, where $X_{0}$ comprises all fixed-effect dummies and $X_{1}$ contains the remaining covariates (including the intercept). Let us denote the aggregate fixed-effect term by $\Lambda = X_0'\beta_{0}$, and the linear projections $P_0$ (onto $span(X_0)$) and its orthogonal complement $M_0$.
This estimator is consistent and asymptotically normal only if every fixed-effect (FE) dimension grow sufficiently fast with the sample size, and some modified overlap condition. If any FE dimension remains fixed, an incidental-parameter problem (IPP) emerges.\footnote{A formal proof of this estimator's asymptotic properties parallels that of Algorithm (ref) and is hence omitted.} The issue can be observed from the estimating equations.
The transformed regression used to update $\beta$ is
while the fixed effects are updated using the projection onto $span(X_0)$:
with $P_0\Lambda = 0$. The term $P_0\overline{\upsilon}$ is the sample group means of the transformed error. Under Assumption (ref), these group means converge to zero whenever each FE group size grows to infinity with the sample, so (ref) identifies $\Lambda$ asymptotically and the noise carried into the $\beta$-update through $\hat\Lambda$ disappears. However, if at least one FE dimension is short, the associated group means do not converge to zero, and neither does $P_0\overline{\upsilon}$, leaving asymptotic bias in $\hat\Lambda$. This noise contaminates the $\beta$-update and produces an IPP: even as $n\to\infty$, the limiting distribution of $\hat\beta$ is biased. This is the same mechanism behind the IPP in PPML: with a short FE dimension, the projection step cannot “average out” the FE component.
We now introduce an estimator that removes the short fixed-effects dimensions ex-ante and thereby sidesteps the IPP. The idea is simple: apply the projection \(M_0\) to both the covariates and the transformed outcome, and then iterate OLS entirely in the residualized space. No fixed effect is ever estimated, so noise from short FE groups cannot feed back into \(\hat\beta\). Theorem (ref) applies with minor modifications, assuming the following overlap condition.
Assumption (ref) requires that a non-trivial share of the within-group sample variation of the residualized regressors is carried by units with $Y_i>0$. This rules out within-group separation patterns.
The resulting estimator is based on the transformed model
where the projector $M_0$ is computed once and reused throughout, so the presence of high-dimensional fixed effects has very limited impact on the computational complexity of the iterative procedure. The asymptotic properties are established in Theorem (ref).
The properties of iOLS established in Theorem (ref) hence carry over to iOLS-FD, under the added assumption $E[M_{0}X_{1}\mid X_{0}]=0$. This additional orthogonality condition ensures that, once additive fixed effects are partialled out from $X_{1}$, the residual variation in $X_{1}$ is uncorrelated with interactions of those same fixed effects. This condition is automatically satisfied in a one-way fixed-effects model, since $X_0'\beta_0$ is constant within each cell. With multi-way fixed effects, however, the nonlinear transformation $\exp(X_0'\beta_0)$ implicitly generates cross-products of fixed effects that vary within cell and cannot be captured by the linear projection. For instance, if unit and time dummies both enter the model, then a within-transformation of the error will itself contain products of unit and time effects. The resulting projection error is hence proportional to the within-cell variance of $\exp(X_0'\beta_0)$. This bias is hence negligible only if $\exp(X_0'\beta_0)$ spans a narrow range of values within cell.
Should \(X_{1}\) also exhibit genuine latent factors, one may replace the standard within-operator \(M_0\) with a factor-purged projector, following bai2009panel or phillips2024simple, to partial out multiple fixed-effect dimensions.\footnote{For two-way fixed-effects, the best rank-1 approximation of the double-demeaned matrix can quickly be obtained via singular-value decomposition. For larger dimensions, an alternating least-squares routine must be used.} Although we do not pursue this extension here, our algorithm is readily adapted to such settings.
More generally, we view the simplifying assumption \(E[M_{0}X_{1}\mid X_{0}]=0\) as natural whenever the outcome \(Y_{i}\) does not call for an interactive fixed-effects specification. For simplicity, all simulations and applications presented in this paper remove all fixed-effect dimensions using this algorithm, thereby abstracting from the possibility of remaining interaction terms.
In the presence of a mix of large and short fixed-effects dimensions, we can combine both approaches. We can estimate the large dimensions of the fixed effects (e.g., time effects) using Algorithm (ref) while applying Algorithm (ref) only to the small dimensions (e.g., unit effects).\footnote{With multiple small dimensions, again a factor-purged projector will be needed. } In that setup, the transformed error no longer contains cross-terms between units and periods and therefore retains the same exogeneity properties as in the one-way case.\footnote{For instance, in a panel with $T=2$ and $n=100$, one should estimate the time dummies and remove the unit dummies, which would otherwise be poorly estimated given only two time periods.} This leads to the following procedure.
Let us partition the regressor matrix as $X = [X_{0}, X_{1}, X_2]$, where $X_{0}$ comprises the fixed-effect dummies for the large dimensions, $X_{1}$ contains the covariates (including the intercept), and $X_2$ contains the the fixed-effect dummies for the small dimensions to be partialled-out from $\Lambda$. We use the linear projections $P_0$ (onto $span(X_0)$) and $\widetilde P_2$ (onto $span(M_0X_2)$) and their orthogonal complements $M_0$ and $\widetilde M_2$. We define the residualized aggregate fixed-effect term by $\tilde{\Lambda} = \widetilde M_2X_0'\beta_{0}$.
The resulting estimator is based on the transformed model.
and the observation that $\widetilde M_2\Lambda$ can be recovered using the same approach as for iOLS-HDFE noting that $P_0\widetilde M_2 = P_0$ by construction, so again we need $P_0 \overline{\upsilon}_i$ to vanish asymptotically with the sample size.
To ensure global convergence of the doubly residualized map, the overlap requirement used for iOLS-FD must be strengthened to Assumption (ref). For iOLS-FD, it suffices that there is positive overlap on $C$, i.e. a uniform slack away from $\tau_i=1$ along directions spanned by $M_0 X_1$. By contrast, with the aggregate fixed-effects term the Jacobian contains the resolvent $(I-P_0T)^{-1}$ which can amplify weight concentration within FE groups.
Assumption (ref) says that usable variation on the positive-$Y$ part of the data dominates any fixed-effects-induced concentration: specifically, $\eta_C>\kappa_{\mathrm{FE}}$. The first bound ensures that, along any direction spanned by the doubly residualized regressors, the iteration weights never place all mass on boundary observations with $Y_i=0$; there is a uniform slack $\eta_C$ away from that boundary. The second bound limits how much these boundary weights can concentrate within fixed-effects groups; $\kappa_{\mathrm{FE}}$ captures the worst such concentration toward $Y_i=0$ within any group. Consequently, when some FE groups contain many zeros, the share of variation contributed by units with $Y_i>0$ must be correspondingly larger. In particular, groups composed entirely of zeros are ruled out.
The asymptotic properties are established in Theorem (ref).
Therefore, the properties of iOLS established in Theorem (ref) hence carry over to iOLS-FD-HDFE, provided FE groups in $X_0$ do not contain too many zero observations.
We assess the performance of our estimators in a series of simulation experiments adapted from fernandez2016individual, weidner2021bias, jochmans2017two, and jochmans2022instrumental. We consider simulations with and without zero observations in the dependent variable, with and without fixed-effects, with and without endogenous regressors, for both count data and continuous outcomes.
We use the Monte Carlo framework of weidner2021bias. Let $Y_{ijt} = \mu_{ijt}\,U_{ijt}$, where: (i) $\mu_{ijt}=\exp\bigl(\alpha_{it} \;+\;\gamma_{jt}\;+\;\eta_{ij}\;+\;X_{ijt}\,\beta\bigr)$; (ii) $U_{ijt}\sim\text{LogNormal}\bigl(\mu=1,\;\sigma^2=\tfrac1{16}\bigr)$; (iii) \(\beta = 1\); (iv) the model-relevant fixed effects \(\alpha\), \(\gamma\), and \(\eta\) are each \(\mathcal{N}(0,1/16)\); (v) $X_{ijt} = X_{ijt-1}/2 + \alpha_{it} + \gamma_{jt} + v_{ijt}$ with $v_{ijt}\sim\mathcal{N}(0,1/16)$; and (vi) the error \(U_{ijt}\) follows one of the following DGP:\footnote{Serial correlation within pairs is specified as \[\mathrm{Cov}\bigl[U_{ijs},U_{ijt}\bigr] = \exp\!\Bigl[ 0.3^{\lvert s-t\rvert} \,\sqrt{\ln\bigl(1+\sigma_{ijs}^{2}\bigr)} \,\sqrt{\ln\bigl(1+\sigma_{ijt}^{2}\bigr)} \Bigr] - 1.\]}
Under DGP I, the conditional variance remains constant, similarly to a Gaussian model with i.i.d.\ errors. In DGP II, the variance tracks the conditional mean, as in a Poisson setting. DGP III is the “log-homoskedastic” case in which variance is proportional to the square of the mean, yielding homoskedastic errors in the log-linear model. Finally, DGP IV imposes a variance structure that combines the features of DGP II and III and adds dispersion as a function of \(X_{ijt}\). None of these DGPs have zero values in the outcome variable.
We run 5,000 simulations for $T=2$ and $N=20$, $50$, and $100$, for a total of $N^2\times T$ observations. We report summary statistics in Table (ref) as well as kernel density estimates in Figure (ref). In line with weidner2021bias, the average bias for PPML is generally larger for DGPs I and IV than II and III. This bias is decreasing in $N$ but remain large relative to the estimated standard error. For example in DGP I, the ratio of average bias to standard deviation of $\hat{\beta}^{PPML}$ is $.07/.129 = .54$ for $T=2,N=20$, $.028/.052 = .54$ for $T=2,N=50$, and $.016/.026 = .61$ for $T=2,N=100$. In contrast, the results confirm the viability of the iOLS estimator across all DGPs, which perform similarly to their Jackknife estimator.
For our second simulation design, we take inspiration from the two-way fixed-effects framework of jochmans2017two, and extend it along the lines of jochmans2022instrumental to allow for endogenous regressors.
Let $Y_{it} = \mu_{it}\,U_{it}$, where: (i) $\displaystyle \mu_{it} = \exp\bigl(X_{it}'\beta + \alpha_i + \delta_t \bigr)$; (ii) $U_{it}$ is a mean-one multiplicative error whose distribution and exogeneity restriction varies by DGP; (iii) $\beta = (1,1,0)'$; (iv) the regressors are generated as $Z_{0it} = 1$, and \[ (Z_{1it},Z_{2it})' \sim N\bigl(\mu_Z + \frac{\alpha_i + \delta_t}{2}\times \mathbf{1}(\alpha_i \geq \delta_t),\Sigma_Z\bigr), \quad \mu_Z =
,\; \Sigma_Z =
, \] with $X_{it} = Z_{it}$ in the exogenous regressor cases, and $X_{kit} = Z_{kit} + 0.5\varepsilon_{it} + 0.5\varepsilon_{it}^2$ for $k=1,2$ in the endogenous regressor case; (v) $\varepsilon_{it}\sim N\bigl(-\sigma_\varepsilon/2,\;\sigma_\varepsilon^2\bigr)$ with $\sigma_\varepsilon = 0.3 + 0.03\,Z_{1it}^2 + 0.03\,Z_{2it}^2$ so $E(\exp(\varepsilon)|Z)=1$; (v) $U_{it} = \exp(\varepsilon_{it})\nu_{it}$, where the error term $\nu_{it}$ is generated differently across each of the seven DGPs:
\[
\]
DGPs I to III produce count data. The variance is larger in negative binomial models compared to the Poisson model, and even more so in the NegBin 5 case. DGP IV (LogNormal 1) is a log-normal design with homoskedastic errors and no zero observations. In contrast, DGP V (LogNormal 0) has zero observations and heteroskedastic errors introduced using a Bernoulli error $\nu_{it} \sim \text{B}(1,P(Z_{it}))$, where the success probability is specified as a logistic function so we have the conditional mean of the error equal to $E(\nu_{it}|Z_{it})=1/P(Z_{it})$ and $E(Y_{it}|Z_{it})=\mu_{it}$ continues to hold. DGPs VI and VII generate produce a mix of zero observations and continuous observations, inspired by santos_silva_2011.
We consider three variants of each of the 7 DGPs:
We simulate 10,000 times each of the 21 DGPs for a sample size of $1,000$ observations with $N=100$ and $T=10$. We report the mean and standard deviations in parentheses for a selection of estimators.\footnote{All simulations were performed in MATLAB 2021b with a 3.6GHz 10-Core Intel Core i9 processor and 32 GB 2667 MHz DDR4 memory.}
For the baseline DGPs, Table (ref) reports summary statistics of: (i) iOLS estimates, which are numerically equivalent to GPML estimates;\footnote{GPML is generally estimated using a Newton-Raphson algorithm or IRLS. However, these two algorithms deliver biased estimates with fixed-effects and do not accommodate instrumental variables. They also do not converge in many simulations when FEs are added.} (ii) OLS estimates of the log-linear model after discarding zero values in the dependent variable; (iii) Popular Fix (PF) estimate using the $\log(Y+1)$ transformation; and (iv) PPML estimates.
We find that iOLS and PPML provide unbiased estimates and have similar precisions across all designs. In comparison, OLS and the Popular Fix (PF) yield biased estimates in most cases. OLS is unbiased only for Lognormal 1, in which there are no zero observations and errors are homoskedastic. Note that OLS and PF estimates of $\beta_2$ are unbiased in all designs, except in Lognormal 0 in which the probability $Pr(Y=0|X)$ is a function of the irrelevant regressor $X_2$. This result illustrates the selection bias problem arising from the presence of zeros.
Table (ref) shows the results for the simulations with fixed-effects. We report summary statistics for iOLS, PPML, GMM-A and GMM-B which are two estimators designed for two-way fixed-effect gravity models developed in jochmans2017two.\footnote{These estimators aim at addressing the IPP bias in standard errors exhibited by PPML in two-way gravity models. Also, PPML can be slow to compute in large samples as the number of parameters to estimate grows with the sample size unlike the others.} iOLS behave very closely to GMM-A in terms of bias and precision. PPML fails to converge in around 2% of simulations.\footnote{This non-convergence rate ranges from about 0.05% in Lognormal1 to about 9% for Lognormal 0, Negbin 1 and Mixture 1.} Discarding these outliers, PPML still appears to be less precise than the other methods.
Finally, Table (ref) reports the results for the simulations with fixed-effects and endogenous regressors. We report i2SLS, the DIFF estimator from jochmans2022instrumental, the PPML-IV estimator from santos_mulhally, and the NL-IV estimator from mullahy1997instrumental.\footnote{In addition to their DIFF estimator, jochmans2022instrumental also develops asymptotic-bias corrections for both PPML-IV and NL-IV, which we do not use here. Note that all three estimators are designed for two-way exponential regression model on a balanced panel data set.}
We find that i2SLS and DIFF yield unbiased estimates with similar precisions. In comparison, PPML-IV and NL-IV fail to converge or deliver outlying estimates in around 0.5% and 4% of simulations, respectively. Neglecting these outliers, both estimators exhibit some bias and are relatively less precise than the other two methods. In terms of computational efficiency, i2SLS with fixed-effects takes about 0.3 second, with linear time complexity in $n$.
We demonstrate the generality of our estimators through three applications beyond gravity models.
We revisit the analysis of michalopoulos_etal_2013, a leading empirical study on the long-run relationship between pre-colonial political centralization and contemporary African development. The authors show that regions historically governed by more centralized ethnic institutions, such as chiefdoms and kingdoms with hierarchical layers of authority, are in contemporary Africa significantly more developed than areas that were organized as stateless societies. Their empirical strategy relates satellite-based measures of regional economic activity to Murdock's index of jurisdictional hierarchy, a proxy for this dimension of institutional complexity, using log-linear regressions with the popular fix transformation.
Night light density at the regional level, denoted $Y_i$, is used as a proxy for economic development, and transformed into $\log(\text{Y}_{i}+0.01)$ to “use all observations and [...] minimize the problem of outliers”. The effect of political centralization is measured by the coefficient associated with jurisdictional hierarchy, denoted $H_i$.\footnote{The index takes discrete values between 0 and 4, it corresponds to the number of jurisdictions above the local level for each ethnicity as reported in 1967. A large number indicates a more centralized political organization.} The cross-sectional unit is ethnicity-by-country.
The main results of the original article are reported in Table III Panel (A), using the popular fix with all observations, and Panel (B), focusing on the intensive margin, i.e. with non-zero observations. We replicate only the columns (4) of each panel from their Table III, which use all available control variables. The corresponding exponential mean model is written as
where $X_i$ includes controls for population density, location, geography, as well as country fixed effects. Although the authors do not claim causality, the causal interpretation of $\beta$ as being an $ATE_{norm}$ of political centralization on economic development would hinge on the strict exogeneity assumption $E(U_i-1|H_i,X_i) = 0$. Under this assumption, both PPML and iOLS can be used to estimate $ATE_{norm}$. However, they use different unconditional moment restrictions and may deliver different point estimates in any given sample. In particular, PPML estimates tend to be sensitive to very large values of the dependent variable.
Table (ref) reports competing estimates of $\beta$ including: (1) the author's choice of the popular fix $\log(Y+0.01)$, (2) the popular fix with a larger constant $\log(Y+1)$, (3) the popular fix with a smaller constant $\log(Y+10^{-5})$ (4) PPML, (5) OLS estimates of the log-linear model after discarding zeros, and (iv) iOLS. To illustrate the sensitivity to large values, we report the estimation results for different samples: (a) the full sample of 682 observations ($\bar{Y}=0.36$, Var$(Y)=2.33$), (b) the full sample except the maximum observation ($Y_{max}=25.14$), (c) the subsample below the 98th percentile ($Y_{98\%}=3.98$); along with the variance-mean ratio (VMR). We then repeat (a)-(c) after dropping zeros to focus on the intensive margin. Values in bold correspond to results in the original paper.
The point estimates derived from the $\log(Y + \Delta)$ transformations differ in magnitude and statistical significance depending on both the choice of shift parameter ($\Delta = 10^{-5}, 0.01, 1$) and the subsample analyzed. Compared to the authors' selected value of $\Delta = 0.01$, employing a smaller shift parameter slightly affects the magnitude of the point estimates but noticeably reduces their statistical significance across all subsamples. Conversely, adopting a larger shift parameter significantly enhances statistical significance while simultaneously shrinking the point estimates. A key implication is that the arbitrary selection of the shift parameter can substantially alter both the estimated coefficients and their statistical significance.
Nevertheless, in this application, the popular fix transformation chosen by the authors yields results comparable to those from the OLS estimation of the log-linear model for subsamples without zeros, and aligns closely with iOLS estimates across all subsamples. Therefore, using this transformation has minimal implications for the paper's main results.
At the intensive margin, the estimated $ATE_{\%}$ is 0.149 and statistically significant at the 5% level, whereas iOLS yields $ATE_{norm} = 0.109$, which is not statistically significant. The iOLS estimate becomes significant and closer to the OLS estimate once the maximum value of $Y$ is excluded. This greater sensitivity of iOLS compared to OLS is unsurprising, since trimming zeros or extreme values mechanically affects the denominator in the computation of $ATE_{norm}$. Nonetheless, both approaches deliver consistent qualitative conclusions about the relationship of interest. The main takeaway is that empirical estimates of $ATE_{\%}$ and $ATE_{norm}$ can be fairly similar. The advantage of $ATE_{norm}$ is that it can also be estimated on the full sample, thereby showing that the main effect persists even when zeros are included, though it may remain sensitive to outliers.
Finally, PPML estimates are of opposite sign relative to the other estimators and are not statistically significant, whether estimated on the full sample or after excluding zeros. The variance-mean ratio exceeds 6 in both cases and is heavily influenced by extreme observations. In fact, removing the single largest observation reverses the sign of the PPML point estimate, and trimming the top 2% of values yields estimates that are closer to those from other methods.
azoulay show that the premature death of a star scientist reshapes her research subfield: collaborators experience a sharp and persistent decline in productivity and funding, whereas non-collaborators increase their output and secure more grants.
Their study considers 452 biomedical researchers who died unexpectedly between 1975 and 2006, across 34,218 subfields. To identify causal effects, they implement a staggered difference-in-differences design:\footnote{In revisiting this application, we abstract from potential concerns related to staggered difference-in-differences designs chaisemartin_haultfoeuille_2020,callaway2021difference,BellegoBenatiaDortetBernadet2024.} each star scientist's death is linked to a specific subfield, which is then paired with a control subfield unaffected by such a death using a matching procedure. The analysis is conducted at the subfield-year level.
Treatment is defined as belonging to a subfield affected by a star scientist's death ($\mathbbm{1}(\text{Treated}_i)=1$) in the years following the event ($\mathbbm{1}(\text{After-Death}_t)=1$). The main specification can be written as:
where $Y_{it}$ denotes the outcome for subfield $i$ in year $t$. The coefficient of interest, $\tau$, captures the differential post-death change in treated subfields relative to controls. The term $\theta$ measures any general shift occurring after a scientist's death that affects both groups. Fixed effects include subfield ($\alpha_i$), subfield age in year $t$ ($\alpha_{it}^{\text{age}}$), and calendar year ($\alpha_t$).
The outcomes considered are the number of publications and the number of NIH grants in each subfield-year, with further disaggregation by whether outputs are linked to the deceased scientist's collaborators or to non-collaborators. A notable feature of the data is the prevalence of zeros: 31.4% of subfield-years record no publications, and 56.7% record no NIH grants.
Their main findings are reproduced in Table (ref) under the row “PPML” of Panel A, reported in bold. They find that following a star researcher's death, the overall amount of publications in the subfield increases by 5.23% but the number of publications of her collaborators fall by 20.7%.\footnote{In the text, we report effects as $\exp(\beta)-1$, using estimates from Table (ref), consistent with a binary treatment variable.} In contrast, non-collaborators publish 7.8% more. Similarly, the star researcher's subfield does not observe a statistically significant change in the number of provisioned NIH grants. They nonetheless find that subsequent to the star's death, collaborators earned 23.2% fewer grants when non-collaborators see their funding rise by 11.6%.
In Panel B. and C., we decompose the effect of a star researcher's death into the extensive and intensive margins, respectively. Results from a linear probability model (LPM) show that collaborators are more likely to exit the subfield entirely, as reflected by publishing zero papers and receiving no NIH grants after the death. By contrast, there is no statistically significant intensive-margin effect for collaborators, suggesting that the negative impact of a star's death on their output is primarily driven by subfield exit rather than reduced productivity among those who remain. For non-collaborators, the opposite pattern emerges: those who are already active increase both their publication output and their grant success.
The choice of estimator affects the magnitude of these effects. Compared with the authors' preferred specification (PPML), the widely used log-linear transformation (“LOG 1+” and “LOG”) produces estimates of smaller magnitude, though qualitatively similar in sign. By contrast, iOLS generally yields larger magnitudes than PPML, with standard errors of comparable size. It is worth noting that in this setting, because the treatment variable is binary, coefficients from log-linear regressions are negatively biased in the presence of heterogeneous treatment effects. In addition, PPML suffers from an incidental parameter bias in both point estimates and standard errors in this three-way fixed-effects setting, which is not addressed here.
The key takeaway is that, in settings with many zeros, complementing total effects with extensive- and intensive-margin estimates provides a clearer picture of the relationship of interest.
inequality study the relationship between innovation, top income inequality, and social mobility. Their main findings are that innovation raises the income share of the top 1%, has little impact on broad inequality, and enhances social mobility when driven by entrants.
Their analysis is based on a balanced panel of 51 U.S. states observed annually from 1976 to 2011, which records top income inequality-measured by the income share of the top 1% of earners in each state-year, and several indicators of innovation intensity constructed from patent data.\footnote{There are six different measures of innovation defined at the year by state level : (1) Patents measures the number of patents per capita, (2) Cit5 measures the number of patents per capita weighted by citations within 5 years of the patent application, (3) Claims measures the number of patents per capita weighted by the patents' breadth, (4) Generality measures the number of patents per capita weighted by herfindalh index of the technological classes that cite the patent, and (5-6) Top5 and Top1 measure the number of patents per capita among, respectively, the top 5% and top 1% most cited patents of the year.}
Their baseline specification (reported in their Table 3) is
where $\log(Y_{it})$ is the log share of adjusted gross income earned by the top 1% in state $i$ and year $t$, $\log(\text{Innov}_{it-2})$ is the log of one of the innovation measures lagged by two years, $X_{it}$ denotes state-year controls, and $\alpha_i$ and $\delta_t$ are state and year fixed effects.
In this setting, endogeneity may arise because inequality itself can influence innovation: higher top incomes may foster lobbying or the erection of barriers to entry, thereby altering innovation dynamics and biasing OLS estimates. To address this concern, the authors extend their specification by adding additional controls and instrumenting for innovation intensity with time-varying state representation on the U.S. Senate Appropriations Committee, which directs federal research funding to universities and stimulates local patenting, while being plausibly exogenous to contemporaneous state-level inequality.
Although the dependent variable is always strictly positive, this application is useful for two reasons: (i) to illustrate our i2SLS estimator, and (ii) to discuss the treatment of regressors that may take the value zero. To handle the latter, the authors correctly estimate
where the indicator function separates observations with positive and zero innovation intensity. However, they only instrument $\mathbbm{1}_{\{\text{Innov}_{it-2} > 0\}}\cdot \log(\text{Innov}_{it-2})$ and treat $\mathbbm{1}_{\{\text{Innov}_{it-2} = 0\}}$ as exogenous. This is problematic if the incidence of zero innovation is correlated with unobserved determinants of inequality, in which case the estimate of $\tau$ may be biased unless the zero-innovation indicator is also instrumented.
We reproduce their main findings in Table (ref), where Panel A reports results under covariate exogeneity and Panel B presents 2SLS and i2SLS estimates. OLS indicates a positive and statistically significant association between the income share of the top 1% and the various innovation measures (LOG in Panel A.). PPML yields coefficients of slightly smaller magnitude with comparable standard errors, though one measure loses statistical significance. By contrast, iOLS produces estimates close to zero and statistically insignificant for three of the six innovation measures, all of which are not focused on top-performing patents.
Turning to Panel B, the instrumental-variable estimates from i2SLS retain the same sign and have standard errors comparable to the 2SLS estimates (LOG), but their magnitudes are smaller, rendering the effects statistically insignificant across all specifications.
From a theoretical perspective, it is not clear which results should be preferred, yet the discrepancy itself is informative. The two estimands capture different objects: $ATE_{\%}$ averages individual semi-elasticities, while $ATE_{norm}$ corresponds to the semi-elasticity of the conditional mean outcome. They coincide when treatment effects are proportional to baseline outcomes, but they generally differ under heterogeneity. A useful rule of thumb is that the sign of the gap depends on how marginal effects correlate with the baseline outcome: if larger effects occur in states with low initial inequality, then $ATE_{\%}$ will tend to exceed $ATE_{norm}$; if they occur in states with high inequality, the reverse will hold. In this application, relatively stronger effects in low-inequality states could explain why the i2SLS estimates appear smaller than those from 2SLS.
Building on our literature review, the formal results introduced in Section (ref), and practical lessons from our applications in Section (ref), this section offers intuitive guidance for applied researchers choosing among alternative estimands and estimation strategies. In particular, we focus on settings where the object of interest is a (semi-)elasticity, a treatment effect expressed in percentage terms, or the conditional mean of a non-negative outcome. The discussion highlights the main considerations that should guide empirical practice in these contexts.
The first step is to define the target parameter. While $ATE_{\%}$ is often a natural choice, it has important limitations. By contrast, $ATE_{norm}$ is always well defined and should generally be reported alongside $ATE_{\%}$. It is particularly relevant when treatment effects are heterogeneous or when the conditional mean function $E(Y_i|\cdot)$ is the primary object of interest.
When outcomes can take the value zero, $ATE_{\%}$ is not identified unless attention is restricted to the intensive margin under suitable exogeneity conditions. By contrast, $ATE_{norm}$ is well defined on the full sample and captures the total effect, while restricting estimation to positive outcomes isolates the intensive margin. This decomposition clarifies whether the overall impact operates primarily through the extensive or the intensive margin. The extensive margin may also be of independent interest and can be estimated separately, for example with a linear probability model.
When the treatment variable is binary, $ATE_{\%}$ estimates may be biased in the presence of heterogeneous treatment effects or heteroskedasticity, as detailed below. In such cases, $ATE_{norm}$ offers a useful alternative. Similarly, when the objective is to estimate parameters of an exponential conditional mean model derived from theory, e.g. in gravity models or wage determination models, only $ATE_{norm}$ corresponds to the structural parameters of interest.
Ultimately, the choice of a target parameter hinges on plausible exogeneity assumptions. Causal interpretation requires strict exogeneity between the error term and regressors (or valid instruments), and these conditional moment restrictions provide not only identification but also an opportunity for specification testing DovononGospodinov2024.
Once the target parameter is defined, the next step is to choose an estimation method. If $ATE_{\%}$ is the parameter of interest, it is typically estimated using OLS on a log-linear regression model. Heteroskedasticity and heterogeneous effects are the two main sources of concern.
\paragraph{Structural parameters.} Suppose a theoretical model implies that the conditional mean takes the form $E[Y_i|X_i] = \exp(X_i'\beta)$. In this case, OLS on $\log Y_i$ does not consistently estimate $\beta$ in the presence of heteroskedasticity. The reason is that the exponential mean specification requires the exogeneity restriction $E[U_i|X_i]=1$, whereas log-linear regression relies instead on $E[\log(U_i)|X_i]=0$. These two conditions coincide only under homoskedastic log-normal errors (or under stronger independence assumptions more generally).
This problem can also be understood as a consequence of Jensen's inequality: $\exp(E[\log Y_i|X_i]) \neq E[Y_i|X_i]$ whenever errors are heteroskedastic on the log scale. As manning2001estimating emphasize, the log-linear model therefore yields biased estimates of $E[Y_i|X_i]$, except in restrictive cases. manning1998logged show that corrections are possible under log-normality by modeling the variance, but extending such adjustments to more general error distributions is substantially more difficult.
Evidence from both applied and theoretical work reinforces this implication. In gravity models, santos_silva_tenreyro_2006 show that OLS on $\log Y_i$ delivers biased estimates of $\beta$ whenever heteroskedasticity is present. Similarly, Blackburn2007 highlight the same issue in the context of wage equations. The central message is that OLS on logs does not consistently recover the structural parameter $\beta$ under the exponential mean specification.
\paragraph{Binary regressors.} Even for reduced-form parameters, similar issues arise when the regressor of interest is binary. In this case, the log-linear regression identifies only an approximation to $ATE_{\%}$. With homogeneous effects and homoskedastic log-normal errors, it is possible to recover $ATE_{\%}$ using the retransformation $\exp(\beta)-1$. When errors are heteroskedastic, however, more elaborate corrections are required, typically involving estimates of the variance function Kennedy1981. This bias diminishes asymptotically as the sample size increases, and is therefore frequently disregarded in empirical applications.
A further problem arises in the presence of heterogeneous treatment effects, even if errors are homoskedastic. The true average percentage effect is
whereas the log-linear regression (or comparing log outcomes between treated and control units) yields an estimate of $\widetilde{ATE}_{\%} = E[\log(Y_i(1)/Y_i(0))] $. Applying the standard transformation yields
By Jensen's inequality, $\exp(E[\log T_i]) \leq E[T_i]$ for any positive random variable $T_i$. Hence $ATE_{\%} \geq \widetilde{TATE}_{\%}\geq \widetilde{ATE}_{\%}$. Thus, with binary regressors and heterogeneous effects, the log-linear approach identifies a downward-biased approximation of $ATE_{\%}$ even in randomized experiments. Formally, a second-order expansion shows that the leading bias is proportional to the variance of the heterogeneous effects. By contrast, $ATE_{norm}$ remains robust to heterogeneity and heteroskedasticity and is therefore a more reliable target parameter.
Several methods are available to estimate $ATE_{norm}$, each with limitations arising from small-sample bias, overdispersion, fixed effects, endogeneity, or sensitivity to outliers.
Our main contribution is to develop alternative algorithms for GPML that accommodate instrumental variables and high-dimensional fixed effects, thereby overcoming key shortcomings of earlier approaches. In applied work, however, the dominant estimator remains PPML, largely due to its widespread use in gravity models.
In what follows, we compare the respective strengths and weaknesses of PPML and GPML. Other generalized linear models with additional parameters, such as negative binomial models, are also available but fall outside the scope of this discussion.
\paragraph{Small-sample bias.} Monte Carlo evidence suggests that the PPML estimator is less prone to small-sample bias than the GPML estimator, although both are consistent under the conditional mean restriction $E(U_i|X_i)=1$ santos_silva_2011. From a theoretical perspective, the small-sample bias of generalized linear model estimators arises from higher-order terms in the expansion of the score equations. CordeiroMcCullagh1991 show that the maximum likelihood estimator of $\beta$ in GLMs generally suffers from an $O(n^{-1})$ bias, the magnitude of which depends on the variance function and the choice of link.
In the case of GPML, where the variance is quadratic in the mean and the link is non-canonical for that variance structure, the first-order bias in $\hat{\beta}$ does not cancel and can be appreciable in finite samples. By contrast, in PPML the log link is canonical and the Poisson likelihood has fixed dispersion. In this canonical case, the $O(n^{-1})$ bias term cancels at first order, leaving only higher-order terms. This explains why PPML estimates of $\beta$ are generally less prone to small-sample bias.
\paragraph{Overdispersion.} Overdispersion in GLMs does not bias coefficient estimates, but it does invalidate conventional standard errors if uncorrected. PPML assumes equality between the conditional variance and the conditional mean of the dependent variable. Violations of this assumption play the same role as heteroskedasticity in OLS: the estimator remains consistent, but it is no longer efficient, and robust standard errors must be used for valid inference. As emphasized by manning2001estimating, the associated efficiency loss can be substantial when log-scale errors are heavy-tailed, in which case alternative estimators may deliver important gains in precision.
Tests for overdispersion are available CameronTrivedi1990, and if it is found to be important, one may prefer specifications that explicitly allow for it, such as GPML which has a quadratic variance-mean relationship, or negative binomial models with an additional dispersion parameter.
\paragraph{Fixed-effects.} As discussed in the introduction, fixed effects in GLMs give rise to an incidental parameter problem. For PPML, the severity depends on dimensionality: with two-way fixed effects, coefficients remain consistent but standard errors are biased fernandez2016individual; with three-way fixed effects, even point estimates are asymptotically biased in general weidner2021bias. In the context of gravity models, recent work has made PPML particularly well suited: jochmans2017two develop a pseudo-demeaned representation that eliminates nuisance parameters in two-way panels, while weidner2021bias derive analytical corrections for three-way panels. Outside these specific trade settings, however, both PPML and GPML remain exposed to the incidental parameter problem. By contrast, our iOLS procedure circumvents it altogether: fixed effects are partialled out at each iteration using Frisch-Waugh-Lovell projections, rendering GPML valid in general multi-way environments without the need for case-specific corrections.
\paragraph{Endogeneity.} Endogeneity poses distinct challenges for PPML and GPML. Because PPML relies on an additive error specification, the implied moment conditions differ from those in the multiplicative case, so that a given instrument will not generally be valid across both formulations. Recent work has shown how instrumental variables can be accommodated in two-way gravity models jochmans2022instrumental. Outside such tailored designs, however, we recommend i2SLS, which is a general and computationally simple way to extend GPML to arbitrary fixed-effects structures with endogenous regressors.
\paragraph{Robustness.} A final consideration in the choice between PPML and GPML estimators concerns their sensitivity to influential observations. In GLMs, the moment conditions differ in a way that affects robustness. For PPML, the score equation $\sum (Y_i - \exp(X_i'\beta))X_i=0$ assigns relatively larger weight to observations with high $Y_i$, since reducing large discrepancies between actual and fitted values requires a proportionally greater adjustment of the coefficient estimates. Formally, the influence function is proportional to this score and is unbounded, which means that outliers in either the dependent or explanatory variables can exert disproportionate influence on the estimates CantoniRonchetti2001. Although the same non-robustness holds for GPML, its moment conditions distribute weights more evenly across observations and are therefore less dominated by very large outcomes. This contrast reflects differences in the underlying variance functions.
Neither estimator is formally robust. In practice, applied researchers often rely on diagnostic checks or sensitivity analyses, such as trimming extreme observations, to assess whether results are driven by outliers. Robust GLM theory also provides formal methods, including bounded quasi-likelihood methods, which mitigate the effect of extreme observations KunschStefanskiCarroll1989.
This paper contributes to the long-standing discussion on how to deal with logs and zeros in regression models. We show that no single solution is universally appropriate: the choice of estimator must depend on the empirical context and the parameter of interest. Our proposed estimator, iOLS, provides a theoretically grounded analogue to the popular $\log(1+Y)$ transformation, combining computational simplicity with consistency for the normalized average treatment effect, while accommodating both fixed effects and endogenous regressors.
Future work may build on this framework in several directions. Iterated linearization can be adapted to other GLMs with log links, to nonlinear difference-in-differences settings, or to regularized estimators such as the LASSO. Further progress may also come from specification testing and model selection procedures tailored to the log-of-zero problem.