# Dealing with Logs and Zeros in Regression Models

Christophe Bellégo CREST (UMR 9194), ENSAE, Institut Polytechnique de Paris, 5
Avenue Henry Le Chatelier, 91120 Palaiseau, France (e-mail:
christophe.bellego@ensae.fr). David Benatia HEC Montréal, Département
d’Économie Appliquée, 3000 Chemin de la Côte-Sainte-Catherine, Montréal, QC
H3T 2A7, Canada (Corresponding author, e-mail: david.benatia@hec.ca). Louis-
Daniel Pape CREST (UMR 9194), École Polytechnique, Institut Polytechnique de
Paris, 5 Avenue Henry Le Chatelier, 91120 Palaiseau, France (e-mail:
louis.pape@ensae.fr).

###### Abstract

Log-linear models are prevalent in empirical research. Yet, how to handle
zeros in the dependent variable remains an unsettled issue. This article
clarifies it and addresses the “log of zero” by developing a new family of
estimators called iterated Ordinary Least Squares (iOLS). This family nests
standard approaches such as log-linear and Poisson regressions, offers several
computational advantages, and corresponds to the correct way to perform the
popular $\log(Y+1)$ transformation. We extend it to the endogenous regressor
setting (i2SLS) and overcome other common issues with Poisson models, such as
controlling for many fixed-effects. We also develop specification tests to
help researchers select between alternative estimators. Finally, our methods
are illustrated through numerical simulations and replications of landmark
publications.

Keywords: Contraction mapping, Elasticity, Gravity model, Iterative estimator,
Log-linear, Selection bias.
JEL: C26, C52, C55.

††The authors thank Ghazala Azmat, Michael Bognanno, Marine Carrasco, Sergio
Correia, Laurent Davezies, Xavier D’Haultfœuille, Philippe Gagnepain, Ben
Groom, James McKinnon, Mathieu Marcoux, Ioana Marinescu, Arnaud Maurel,
Isabelle Méjean, Joris Pinkse, João Santos Silva, Alain Trognon, Michael
Visser, Martin Weidner, Jeffrey Wooldridge, Tom Zylkin, and all participants
to seminars and conferences for insightful comments. STATA packages are
available at https://github.com/ldpape/iOLS_delta.

##  1 Introduction

Having to deal with the (natural) logarithm of a zero in the dependent
variable is a recurring problem. There is, unfortunately, a lack of consensus
regarding the best way to address those zeros in log-linear and log-log
models, as evidenced by the many alternative approaches used in recent leading
publications. Besides, these existing approaches may not always be suited to
the data and may suffer from a range of possible issues.111In this paper, we
focus on the log-linear model and address the minor differences of the log-log
model as an extension in Appendix B.4.

This paper not only discusses the existing methods and their limitations for
addressing the “log of zero”, it also develops a new family of estimators and
a model selection procedure which overcome these limitations. Our estimators
are simple iterative extensions of ordinary least squares (OLS) and two-stage
least-squares (2SLS). They are consistent, asymptotically normal,
computationally simple, and can accommodate many fixed-effects along with
instrumental variables. Our approach is general in the sense that it nests a
continuum of models, including the log-linear model and Poisson regression as
special cases, and can be interpreted as the correct way to perform the
$\log(Y_{i}+1)$ transformation. We then develop a specification test aimed at
selecting the most suitable approach to address the log of zero. This test
consists in verifying the external validity of each model with respect to the
observed pattern of zeros in the data. We believe this new family of
estimators, combined with the specification tests, provides a natural
framework for handling zeros in the dependent variable.

Figure 1: Prevalence of the Log of Zero in the AER (2016-2020)

To document the ubiquity of logs and zeros in regression models, and the lack
of consensus thereof, we have collected information from three sources: a
review of empirical publications, a survey conducted during seminars, and the
records from an online research forum. First, we have reviewed all articles
published in the American Economic Review (AER) between 2016 and 2020. Figure
1 summarizes our findings. We find that log-linear and log-log models are
among the most frequent specifications used in empirical research. Nearly 40%
of empirical papers used a log-specification and 36% of these faced the
problem of the log of zero. Several solutions are employed in practice. In
most publications, the authors chose to keep the zero observations and opted
to either (1) add a positive discretionary value to the dependent variable
(48%), (2) use Poisson-type estimators (35%), or (3) apply the inverse
hyperbolic sine (IHS) transformation (15%). Discarding non-positive
observations occurred in 31% of publications. We also note that in around 20%
of cases, the authors compared several methods in order to gauge the
robustness of their results.222This excludes cases where the authors decided
to use a linear specification by fault of having to use such a fix. See Table
D.2 in the Appendix for additional details and information regarding data
collection. Second, we have also conducted a survey in three online seminars
in economic departments asking researchers “What would you do when facing the
log of zero?”. The distribution of responses mirrors that observed in the AER
publications except for the large stated preference for the use of mixture
models.333Appendix D provides more details about the surveys.  Third, the
question “Log transformation of values that include 0 (zero) for statistical
analyses?” asked in 2014 on ResearchGate, a multidisciplinary research-
oriented social network, has been read 120,000 times by August 2020 and has
received 38 contributions from researchers in medicine, biology, statistics,
engineering and other fields. Again, the prevalence of each suggested solution
is comparable to that in the AER.

This evidence suggests three main issues which must be overcome. First, the
problem of the log of zero is widespread in empirical research and extends
beyond the realm of economics. This paper attempts to solve this problem with
a new flexible and computationally-efficient estimator. Second, the prevalence
of simple ad hoc fixes over theoretically-founded methods suggests the
existence of numerical difficulties in some contexts. Our discussions with
empirical researchers indeed revealed that many opt for the popular fix
approach (i.e. $\log(Y+1)$) for its simplicity and convenience, at the cost of
what they assume is a small bias. Our approach, unlike Poisson models, do not
suffer from the convergence issues identified in Santos Silva and Tenreyro
(2011), or difficulties in handling instrumental variables or many fixed-
effects. Third, the choice of one method over another is almost never
justified in publications. In contrast, this paper focuses on the
identification of model parameters and provides a way to substantiate the
choice of a particular model. For many proposed solutions to the log of zero,
like Poisson regression, identification rests on assuming some particular
moment conditions. We show that these conditions correspond to implicit
assumptions about the conditional probability of having a zero in the data.
Based on this observation, we develop a statistical test for evaluating which
model is consistent with the observed pattern of zeros. Our methodological
contribution hence both extends the set of potential solutions while
simultaneously limiting the researcher’s discretionary power.

To handle these issues, we develop a new class of estimators, called iterated
Ordinary Least Squares (iOLS) designed to address the “log of zero”. Our
approach consists in adding an observation-specific value to the outcome,
instead of a constant, which is scaled using a hyper-parameter. This parameter
controls the underlying moment condition used for estimation and can be user-
selected or data-driven using our model selection approach. The range of
admissible moment conditions form a continuum which limits correspond to the
restrictions used in the log-linear model and the Poisson model, respectively.
The model selection procedure amounts to finding the hyper-parameter value
such that the conditional probability of having a zero observation implied by
the model is consistent with the data.

We study the theoretical properties of iOLS (and i2SLS), including consistency
and asymptotic normality, using the asymptotic theory developed in Dominitz
and Sherman (2005). Our estimator corresponds to the fixed point of an
(asymptotic) contraction mapping, which is solved for by running OLS (or 2SLS)
iteratively.444Our model could also be estimated by non-linear least squares
at greater computational costs.  To fix ideas, let us illustrate our
estimation method in its simplest form. Assume $Y$ and $X$ are the dependent
and independent variables, $\beta$ is the parameter of interest, and
$\delta>0$ is the hyper-parameter. The simplest iOLS procedure consists of the
following steps: 1. Initialize $t=0$ and obtain an initial estimate
$\hat{\beta}_{0}$, e.g. by regressing $\log(Y+1)$ onto $X$ using OLS; 2.
Regress
$\log\left(\frac{Y+\delta\exp(X^{\prime}\hat{\beta}_{t})}{1+\delta}\right)$
onto $X$ using OLS to obtain $\hat{\beta}_{t+1}$; 3. Update $t$ to $t+1$; 4.
Repeat steps 2 and 3 until $\hat{\beta}_{t+1}$ converges.555Our approach does
not use Iterated Reweighted Least-Squares, unlike the logit and Poisson models
estimated as generalized linear models, but it can be extended to include a
weighting matrix like in Generalized Least Squares. Remark that this
transformed dependent variable corresponds to the log of a weighted average of
$Y$ and $\exp(X^{\prime}\beta)$, where weights depend on the hyper-parameter
$\delta$.

This estimator has multiple advantages: (a) it can be estimated by ordinary
least squares, hence is computationally fast and easy to implement, with
possibly high dimensional fixed effects; (b) it extends naturally to the
endogenous setting using iterated 2SLS (i2SLS), also in situations with many
fixed effects; (c) it is amenable to different identifying assumptions by
varying the hyper-parameter; (d) robust standard errors are readily available;
and (e) it does not suffer from highly dispersed response variables or the
numerical convergence issues identified for (additive) Poisson estimators.

Our methodological contributions are illustrated through numerical simulations
and (partial) replications of Michalopoulos and Papaioannou (2013) and Santos
Silva and Tenreyro (2006). The former examine the role of pre-colonial ethnic
institutions on economic development by using the popular fix to address the
log of zero, whereas the latter has popularized the use of Pseudo-Poisson
Maximum Likelihood (PPML) for estimating gravity models in trade. Our approach
yields plausible and justifiable estimates in both replications whereas PPML
estimates are found to be externally inconsistent with the observed patterns
of zeros.

The remaining of the paper is organized as follows. Section 2 clarifies the
log of zero issue and discusses existing practices as well as their limits.
Section 3 develops a new family of solutions. Section 4 presents specification
tests and a data-driven model selection procedure. Numerical simulations are
presented in Section 5. The replication exercise is proposed in Section 6.
Section 7 concludes the paper.666The online Appendix also contains several
useful extensions. Appendix B.1 implements Poisson regression as an iOLS
estimator. Appendix B.2 adapts iOLS to the endogenous setting. Appendix B.3
extends it to deal with negative values in $Y$. Appendix B.4 addresses log-log
specifications. Appendix B.5 develops a fast estimator for high-dimensional
fixed-effects. Appendix B.6 deals with the log of a ratio. Appendix B.7 shows
an alternative iOLS approach equivalent to sample-selection models. Finally,
Appendix B.8 details the testing procedures in the endogenous regressors
setting.

##  2  Existing Practices

Log-linear regressions are used in research for many purposes. These include
using the log transformation because (1) the parameter estimate is an
elasticity or a semi-elasticity;777In a log-log model such as
$\log(y)=\beta\log(x)+\epsilon$, the elasticity of $y$ with respect to $x$ is
given by $\frac{\partial\log(y)}{\partial\log(x)}=\frac{\partial y}{\partial
x}\frac{x}{y}=\beta$.  (2) logs can linearize a theoretical model, e.g. a
Cobb-Douglas production function; (3) logs can make heteroskedasticity vanish
in some settings; (4) the data is sometimes naturally related by a log-linear
relationship; or even (5) it provides a concave transformation for a highly
dispersed outcome.

In this section, we first provide a general setup from which we assume the
data to have been generated from. We then review the five main identified
solutions to the “log of zero”, and discuss their assumptions and limitations.
The most popular fix consists in adding a positive constant to all
observations. A second solution is to discard the non-positive observations
from the sample. A third solution uses transformations of the response
variable, such as IHS, akin to the log function. Fourth, the Poisson model
handles the presence of zeros well in many settings. It is especially popular
in international trade for the estimation of gravity equations (Head and
Mayer, 2014). Finally, mixture models (e.g. Tobit or Heckit) also provide a
solution by modelling the occurrence of non-positive observations as a sample
selection problem.

###  2.1 The setup

Let us consider an iid sample of observations $\\{Y_{i},X_{i}\\}_{i=1}^{n}$,
where $n$ denotes the sample size, generated by the “true” model given by

$Y_{i}=\exp(X_{i}^{\prime}\beta+\varepsilon_{i})\xi_{i},$ (1)

where $\beta$ is a fixed parameter of interest in $\mathbb{R}^{K}$, with
$K\geq 1$, $\varepsilon_{i}$ is a random error, and $\xi_{i}\in\\{0,1\\}$ is a
Bernoulli random error. 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$.

$Y_{i}$ can either be equal to zero, when $\xi_{i}=0$, or take positive
values, when $\xi_{i}=1$. Taking logs on both sides of (1) is allowed only if
$Y_{i}$ (and thus $\xi_{i}$) takes only strictly positive values. Doing so
yields the log-linear model given by

$\log(Y_{i})=X_{i}^{\prime}\beta+\varepsilon_{i}.$ (2)

For parsimony, we will focus on the more compact multiplicative
representation,

$Y_{i}=\exp(X_{i}^{\prime}\beta)U_{i},$ (3)

where $U_{i}=\exp(\varepsilon_{i})\xi_{i}$ is the error term, and the related
additive model

$Y_{i}=\exp(X_{i}^{\prime}\beta)+\epsilon_{i},$ (4)

with $\epsilon_{i}=\exp(X_{i}^{\prime}\beta)(U_{i}-1)$ as the error term. We
will discuss the difference between the later two forms when presenting
Poisson models.

###  2.2 A wide variety of solutions

##### The popular fix: $\mathbf{log(Y+1)}$.

The most popular solution is to add a positive constant $\Delta$ to all
observations $Y_{i}$ so that $\tilde{Y}_{i}=Y_{i}+\Delta>0$ making the log-
transformation feasible for $\tilde{Y}_{i}$. The choice of $\Delta$ is
discretionary and may arbitrarily bias the estimates and their standard
errors. Moreover, the size of the bias will depend on the data at hand,
suggesting that adding the smallest possible constant is not necessarily the
least “harmful” choice.888Winkelmann (2008) discusses this approach and Cohn,
Liu and Wardlaw (2021) survey its application in finance.  We reconcile this
approach with theory in Section 3.

To understand the bias, consider the model specified in (1). Adding $\Delta>0$
and applying the log function yields after rearrangement

$\log(Y_{i}+\Delta)=X_{i}^{\prime}\beta+\log\left(U_{i}+\frac{\Delta}{\exp(X_{i}^{\prime}\beta)}\right)$
(5)

where the error term
$\omega_{i}=\log\left(U_{i}+\frac{\Delta}{\exp(X_{i}^{\prime}\beta)}\right)$
is correlated with $X_{i}$ by construction, even when $U_{i}$ and $X_{i}$ are
statistically independent, and creates an endogeneity bias. Although the
choice of $\Delta$ matters, $\exp(X_{i}^{\prime}\beta)$ can be arbitrarily
close to zero hence leading to possibly large biases. Thus, the “popular fix”
estimator is in general not consistent.999This estimator is consistent under
$E(\omega_{i}|X)=constant$ which implies strong assumptions of the joint
distribution of $U_{i}$ and $X_{i}$. Furthermore, standard errors will be too
small because the zero observations transformed into $\log(\Delta)$
artificially reduce the variance of the transformed dependent variable.

##### Discarding zeros.

The simplest solution is to delete the zero observations and estimate (2)
directly with OLS. Formally, discarding zeros introduces a selection bias
unless $E[\varepsilon|\xi=1,X]$ is a constant. Similarly, one could discard
zeros and estimate (4) assuming $E[\exp(\varepsilon)|\xi=1,X]$ to be constant.
Doing so assumes away any role played by the zeros and has context-dependent
consequences; rendering it inadvisable at least since Young and Young (1975).
At the very least, it will change the scope of the study by narrowing down the
focus to observations for which $Y_{i}>0$.

The economic interpretation of the error term should always be discussed when
making such an assumption. For instance, some empirical studies relying on the
Mincer equation for the purpose of estimating the returns to schooling use the
log wage and discard unemployed individuals. Unemployed agents have unobserved
wage rates which can be labelled as zeros. If $\varepsilon_{i}$ captures the
unobserved ability of individual $i$, it will undoubtedly be correlated with
her employment outcome $\xi_{i}=1$ or $\xi_{i}=0$, hence introducing a sample
selection bias when discarding the zeros.

##### Other transformations.

An alternative approach relies on log-like transformations applicable to non-
positive values. The most popular is the IHS or related transformations
(MacKinnon and Magee, 1990; Burbidge, Magee and Robb, 1988; Ravallion, 2017).
It consists in transforming $Y_{i}$ into $\tilde{Y}_{i}=\log(\theta
Y_{i}+\sqrt{\theta^{2}Y_{i}^{2}+1})/\theta$ and estimating
$\tilde{Y}_{i}=X_{i}^{\prime}\beta+\omega_{i}$ by OLS. If the underlying model
writes in log, then this transformation will likely yield biased
estimates.101010Considering model (1), obtaining consistent estimates requires
a moment condition like $E(\log(\theta U_{i}+\frac{\sqrt{\theta
Y_{i}^{2}+1}}{\exp(X_{i}^{\prime}\beta})|X)=0$, which may be difficult to
justify since $\frac{\sqrt{\theta Y_{i}^{2}+1}}{\exp(X_{i}^{\prime}\beta)}$ is
a non-linear function of $X$. Nearly all economic applications set $\theta$ to
1 such that $\tilde{Y}$ tends toward $\log(2Y)$ for large values of $Y$.

This transformation essentially consists in adding a positive _observation-
specific_ value to the response variable before applying the log function. Its
similarity with the $\log$ function may lead to treating them interchangeably.
However, for small values of $Y_{i}$, these transformations can behave
differently. Besides, as shown in Bellemare and Wichman (2020), the
interpretation of the coefficients is not trivial and the underlying
elasticity is potentially biased or undefined.111111The authors show that in
$\tilde{Y}_{i}=X_{i}^{\prime}\beta+\epsilon_{i}$, the elasticity
$\hat{\zeta}_{yx}=\hat{\beta}x\frac{\sqrt{y^{2}+1}}{y}$ is a function of $x$,
$y$, or is not defined for $y=0$. $\beta$ is an elasticity only if $x=1$ and
$y$ is large. It is hence satisfactory in contexts where applying a concave
transformation is the main objective but suffers from similar limitations than
the “popular fix” discussed above.

##### Poisson models.

The model presented in (1) is non-linear in variables and parameters. The
parameters are identified and non-linear estimators, such as non-linear least
squares or Iterated Reweighted Least-Squares (IRLS), yield consistent
estimates of $\beta$ under the strict exogeneity restriction
$E(U_{i}|X_{i})=1$ which implies the unconditional moments121212Choosing the
“best” unconditional moments, or rather picking the optimal instruments, from
a conditional moment restriction is beyond the scope of this paper.

$E(X_{i}(Y_{i}-\exp(X_{i}^{\prime}\beta)))=0.$ (6)

The empirical counterpart yields a solution equivalent to that of the Pseudo
log-likelihood of the Poisson model (Gourieroux, Monfort and Trognon, 1984).
This approach is computationally efficient because it is a well-defined
concave problem. This approach has been popularized by Santos Silva and
Tenreyro (2006) for gravity models and is generally referred to as PPML. It is
based on the additive representation of the model in (4) assuming
$E((U_{i}-1)\exp(X_{i}^{\prime}\beta)|X_{i})=0$, which is equivalent to
$E(U_{i}|X_{i})=1$ in absence of endogenous regressors but leads to a
different objective criterion.

Nevertheless, Poisson models have several shortcomings: (1) Existence of a
solution is not guaranteed (Santos Silva and Tenreyro, 2010);131313The authors
have a dedicated website with helpful resources about PPML
(https://personal.lse.ac.uk/tenreyro/lgw.html). (2) They can exhibit
convergence issues (Santos Silva and Tenreyro, 2011); (3) They can be very
imprecise if the log-scale error is heavy-tailed (Manning and Mullahy, 2001);
(4) They can be difficult to estimate with many fixed-effects (Correia,
Guimarães and Zylkin, 2019); (5) instrumental variables require stronger
assumptions and may dramatically increase computational complexity
(Wooldridge, 2015);141414Non-linear IV estimators require strict exogeneity
between the errors and instruments unlike linear estimators. and (6) they
suffer from an incidental parameter problem in the standard errors when
including two-way fixed effects (Fernández-Val and Weidner, 2016), and in the
point estimates when controlling for three-way fixed-effects. (Weidner and
Zylkin, 2021).

Another important consideration about Poisson models is the role played by
heteroskedasticity. Santos Silva and Tenreyro (2006) conclude that “under
heteroskedasticity, the parameters of the log-linearized models estimated by
OLS lead to biased estimates of the true elasticities” and suggest using PPML
as a solution. This result is not as general as it may seem. It is only valid
under the Poisson restriction $E[U|X]=1$. To show this, we perform a Taylor
expansion of the log-scale error $\varepsilon_{i}$ around $0$ to obtain

$\varepsilon_{i}=\log(1)+(\exp(\varepsilon_{i})-1)-\frac{(\exp(\varepsilon_{i})-1)^{2}}{2}+\frac{(\exp(\varepsilon_{i})-1)^{3}}{3}+...$
(7)

This expression implies that, if $E[\varepsilon_{i}|X]=0$ holds, the
identifying assumption of the Poisson model $E[\exp(\varepsilon_{i})|X]=1$
also holds _only if_ one makes additional assumptions about the higher-order
centered moments of $\exp(\varepsilon_{i})$, i.e. its variance, skewness, etc.
The reverse conclusion hence applies: under heteroskedasticity, the additional
assumptions about the higher-order moments of $\exp(\varepsilon_{i})$ are
unlikely to hold, therefore the Poisson model leads to biased estimates of the
true elasticities. The presence of heteroskedasticity alone hence does not
invalidate one approach over the other because it depends on the underlying
exogeneity restriction, which is unverifiable.

This is not to say that heteroskedasticity is irrelevant. Indeed, Manning and
Mullahy (2001) show that the parameters of the log-linear model, even though
well-identified under $E[\varepsilon_{i}|X]=0$, implies a conditional mean
function $E(Y|X)$ which can be biased under heteroskedasticity, if the
function is not appropriately retransformed, e.g. by using the smearing
estimator (Duan, 1983) or a gaussian approximation (Manning, 1998). In
opposition, Poisson models do not require ex-post retransformation to deliver
an elasticity estimate of the mean response, provided
$E[\exp(\varepsilon_{i})|X]=1$ holds true.

Finally, there are some reasons to prefer the multiplicative Poisson model
over its additive counterpart. The latter (PPML) provides a computationally
efficient solution in many cases but suffers from two significant drawbacks.
First, the estimator minimizes the sum of squared _absolute_ errors
$Y_{i}-\exp(X_{i}^{\prime}\beta)$. Doing so implies that observations with
large values of $Y$ are given much more weight than those with smaller
values.151515They yield different estimates because the additive model uses
the unconditional moment restrictions
$E(X_{i}\exp(X_{i}^{\prime}\beta)(U_{i}-1))=0$ instead of
$E(X_{i}(U_{i}-1))=0$ for the multiplicative model. This is similar to
discarding zero observations in many settings. As such, it cannot be
considered a satisfactory approach to dealing with zeros in the general
setting. This issue is illustrated in our applications.

Second, the additive formulation may be “awkward to motivate directly”
(Mullahy, 1997). This is particularly evident when considering an omitted
variable problem. In such a case, the additive model assumes this omitted
variable to have an additive effect whereas all explanatory variables have
multiplicative effects. The moment conditions for identification in the
additive case hence differ from the usual setting. Finally, elasticities
estimated from the additive model depends on both the exogenous and omitted
variables, unlike those obtained in the multiplicative case.161616These
arguments have been made by Mullahy (1997) and Jeffrey Wooldridge
(https://www.statalist.org/forums/forum/general-stata-
discussion/general/1536854-multiplicative-versus-additive-iv-poisson).

##### Mixture models and Heckman’s correction.

Mixture models consist in modeling the selection explicitly, $Y_{i}=0$ or
$Y_{i}>0$, using a latent variable approach under chosen distributional
assumptions. This approach is rarely used to address the log of zeros but has
been relied upon in the context of gravity equations (Eaton and Tamura, 1994).
The Heckman’s (“Heckit”) correction (Heckman, 1979) applies as follows. In the
setting provided by model (1), it assumes that $\xi_{i}=1$ if
$X_{i}^{\prime}\gamma+\nu_{i}>0$, and $\xi_{i}=0$ otherwise.
$X_{i}^{\prime}\gamma+\nu_{i}$ is hence referred to as the “selection
equation”. The key identifying restriction is that $\varepsilon_{i}$ and
$\nu_{i}$ are bivariate normal, so that $E[\varepsilon_{i}|U_{i}>0,X]$ admits
the closed-form expression
$E[\varepsilon_{i}|\nu_{i}>-X_{i}^{\prime}\gamma,X]=\lambda\frac{\phi(-X_{i}^{\prime}\gamma)}{\Phi(X_{i}^{\prime}\gamma)}$,
for $\phi(\cdot)$ and $\Phi(\cdot)$ denoting the Gaussian probability density
and distribution functions, respectively. $\lambda$ and $\gamma$ are estimable
parameters. Estimation takes two steps. First, a probit model of $Y_{i}>0$
conditional on $X_{i}$ yields $\hat{\gamma}$. Second, the log-linear
regression with an additional term, as specified by

$\log(Y_{i})=X_{i}^{\prime}\beta+\lambda\frac{\phi(-X_{i}^{\prime}\hat{\gamma})}{\Phi(X_{i}^{\prime}\hat{\gamma})}+e_{i},$
(8)

is estimated by OLS to obtain $\beta$ and $\lambda$. The relevance of the
correction term can be tested using a t-test to check whether $\hat{\lambda}$
is different from zero.171717We show how this approach relates to ours in
Appendix B.7. Note that, however, this approach is heavily dependent on the
distributional assumption in absence of instrumental variables in the
selection equation.

##  3  Iterated Ordinary Least Squares (iOLS)

In this section, we develop a new approach which somehow reconciles the
popular fix with econometric theory. This new approach yields a family of
estimators requiring only OLS, and more importantly which embeds both the log-
linear model and Poisson models as special cases. For clarity, we first show
how our estimation procedure arises in the context of the log of zeros.
Second, we present the algorithm. Third, we derive its asymptotic properties,
and detail how minor modifications can accommodate alternative exogeneity
conditions. Finally, we develop a variety of important extensions in the
Appendices.

###  3.1 Fixing the popular fix (iOLSδ)

We let $\Delta_{i}$ vary across observations such that $Y_{i}+\Delta_{i}>0$.
From (5), we have

$\log(Y_{i}+\Delta_{i})=X_{i}^{\prime}\beta+\log\left(U_{i}+\frac{\Delta_{i}}{\exp(X_{i}^{\prime}\beta)}\right).$
(9)

Letting $\Delta_{i}=\delta\exp{(X_{i}^{\prime}\beta)}$, for some $\delta>0$,
this equation becomes

$\log(Y_{i}+\delta\exp(X_{i}^{\prime}\beta))=X_{i}^{\prime}\beta+\upsilon_{i}.$
(10)

where the new error term
$\upsilon_{i}=\log\left(\delta+U_{i}\right)>\log(\delta)$ is unlikely to be
mean-zero. Before proceeding to the more general setting, let us first assume
the following strict exogeneity restriction

$E[\upsilon_{i}|X]=\log(\delta+1).$ (11)

It means that the conditional mean of $\log(\delta+U_{i})$ is independent of
covariates. Under this assumption, (10) can be rewritten into

$\log\left(\frac{Y_{i}+\delta\exp(X_{i}^{\prime}\beta)}{\delta+1}\right)=X_{i}^{\prime}\beta+\overline{\upsilon}_{i},$
(12)

where $\overline{\upsilon}_{i}=\log(\delta+U_{i})-\log(\delta+1)$ is a mean-
zero error term, and $\beta$ could be estimated using a non-linear estimator
or our iterative least-squares estimator to solve

$\min_{\beta}\sum_{i=1}^{n}\left(\log(Y_{i}+\delta\exp(X_{i}^{\prime}\beta))-\log(\delta+1)-X_{i}^{\prime}\beta\right)^{2}.$
(13)

Remark that this transformation consists in taking the log of a weighted
average of $Y$ and its conditional mean, where $\delta$ controls the weights.

##### From log-linear to Poisson.

Interestingly, $\delta$ acts as a model selection parameter within a continuum
of models which includes two important special cases: the log-linear model (as
$\delta\rightarrow 0$) and the multiplicative Poisson model (as
$\delta\rightarrow\infty$), along with a range of intermediate models.

On the one hand, the relation with log-linear models in absence of zeros is
trivial. Let $\xi_{i}=1$ for all $i=1,...,n$, the condition in
$\eqref{eq:newmodel exog1}$ is then equivalent to $E[\varepsilon_{i}|X]=0$ for
$\delta=0$ which is the exogeneity restriction used by the log-linear model
(2), and so the objective in (13) becomes

$\min_{\beta}\sum_{i=1}^{n}(\log(Y_{i})-X_{i}^{\prime}\beta)^{2}.$ (14)

To show this relation when there are zeros, let us develop the error term as

$\log(\delta+U_{i})=(1-\xi_{i})\log(\delta)+\xi_{i}\log(\delta+\exp(\varepsilon_{i})).$
(15)

Taking a Taylor expansion of the second term around $\exp(\varepsilon_{i})$,
substituting in the above expression, and rearranging yield

$\log(\delta+U_{i})=(1-\xi_{i})\log(\delta)+\xi_{i}\varepsilon_{i}+o(\delta).$
(16)

Imposing condition $\eqref{eq:newmodel exog1}$ for $\delta$ small enough is
approximately equivalent to assuming that

$Pr(\xi_{i}=0|X)\log(\delta)+Pr(\xi_{i}=1|X)E[\varepsilon_{i}|X,\xi_{i}=1]=0$
(17)

does not depend on $X$. Therefore, for $\delta$ small enough, our approach
corresponds to the log-linear model where the log-scale errors for zero
observations satisfy the condition
$E[\varepsilon_{i}|X,\xi_{i}=0]=\log(\delta)$.181818One can also show that it
gives a reasonable approximation of the log-linear model as long as
$E[\varepsilon_{i}|X,\xi_{i}=1]$ is proportional to $Pr(\xi_{i}=1|X)^{-1}$, as
illustrated in Section 5’s DGP 2.

On the other hand, taking a Taylor expansion of $\log(\delta+U_{i})$ around
$U_{i}=1$ gives

$E[\log(\delta+U)|X]=\log(1+\delta)+\frac{1}{1+\delta}E[U-1|X]+o(\delta^{-2}).$
(18)

The terms in $o(\delta^{-2})$ will become negligible faster than
$E[U-1|X]/(1+\delta)$ as $\delta$ increases. Thus, for $\delta$ large enough,
the above exogeneity condition is approximately equivalent to $E[U|X]=1$ and
(13) goes to

$\min_{\beta}\sum_{i=1}^{n}\left(Y_{i}\exp(-X_{i}^{\prime}\beta)-1\right)^{2},$
(19)

which corresponds to the optimization problem solved by the (multiplicative)
Poisson model specified in (3). Section B.1 shows that both multiplicative and
additive Poisson models are nested within our general framework and not just
as limiting cases.

###  3.2 Estimation by iOLS

Let us make the following assumptions about the covariates $X_{i}$ and error
term $U_{i}$.

######  Assumption 1 (Covariates)

$X$ has full column rank and $E(X_{i}X_{i}^{\prime})<\infty$.

######  Assumption 2 (Errors)

The aggregate error $U_{i}=\exp(\varepsilon_{i})\xi_{i}$ is independently and
identically distributed.

In the general setting, our approach consists in transforming the response
variable as

$\tilde{Y}_{i}(\beta,\delta)=\log(Y_{i}+\delta\exp(X_{i}^{\prime}\beta))-c(\beta,\delta),$
(20)

where $c(\beta,\delta)$ is a known function of $\beta$ and $\delta$, possibly
different from $\log(\delta+1)$, so as to obtain a (seemingly) linear model:

$\tilde{Y}_{i}(\beta,\delta)=X_{i}^{\prime}\beta+\overline{\upsilon}_{i},$
(21)

with $E[\overline{\upsilon}_{i}]=0$. We refer to this model as iOLSδ, because
it depends on the choice of the parameter $\delta$, which will be further
discussed together with exogeneity restrictions. The moment condition
$E[X_{i}\overline{\upsilon}_{i}]=0$ yields

$\beta=E[X_{i}X_{i}^{\prime}]^{-1}E\left[X_{i}\tilde{Y}_{i}(\beta,\delta)\right],$
(22)

which characterizes $\beta$ as the solution of a fixed-point problem. Based on
this insight, we propose the following iterative least-squares estimator.
Thereafter, we drop the dependence of $\tilde{Y}_{i}$ on $\delta$ to alleviate
notations.

######  Algorithm 1 (iOLS estimator)

The iOLS estimator is defined as follows:

  1. 1.

Initialize $t=0$ and let $\hat{\beta}_{0}$ be an initial estimate, e.g.
obtained with the “popular fix” estimator
$\hat{\beta}^{PF}=[X^{\prime}X]^{-1}X^{\prime}\log(Y+\Delta)\in\mathbb{R}^{K}$,
for some $\Delta>0$;

  2. 2.

Transform the dependent variable into $\tilde{Y}(\hat{\beta}_{t})$ using (20);

  3. 3.

Compute the OLS estimate
$\hat{\beta}_{t+1}=[X^{\prime}X]^{-1}X^{\prime}\tilde{Y}(\hat{\beta}_{t})$,
and update $t$ to $t+1$;

  4. 4.

Iterate steps 2 and 3 until $\hat{\beta}_{t}$ converges.

Remark also that $X^{\prime}X$ needs only be inverted once, making this
approach computationally efficient. In addition, the resulting estimate varies
continuously with $\delta$ which allows using warm starting to solve for a
large range of $\delta$’s, if desired.

###  3.3 Identification and asymptotic properties

The simplified version of iOLS presented earlier assumed
$E[\upsilon_{i}|X]=\log(\delta+1)$. We now turn to the more general setting by
considering the assumption below.

###### Assumption 3

The error term $\upsilon_{i}$ satisfies the weak exogenous restriction
$E[X_{i}^{\prime}(\upsilon_{i}-c(\delta,\beta))]=0$ where the value
$c(\delta,\beta)$ is unknown. In addition, let $E[U_{i}]=1$.

Like in the Poisson model, the assumption that $E[U_{i}]$ is known allows
identifying the intercept term. It also affects $c(\delta,\beta)$ which is
determined by (unconditional) higher-order moments of $U_{i}$. To see this,
consider the Taylor expansion of $\log(\delta+U_{i})$ around $\log(1+\delta)$
to obtain

$c(\delta,\beta)=\log(1+\delta)-\frac{1}{2(1+\delta)^{2}}E[(U_{i}-1)^{2}]+\frac{1}{3(1+\delta)^{3}}E[(U_{i}-1)^{3}]+o(\delta),$
(23)

where the second and third terms are respectively the variance and third
centered moment of $U_{i}$. The first centered moment is zero under Assumption
3. If the infinite weighted sum of higher moments of $U_{i}$ equals zero then
we are back to the simple setting where $c(\delta,\beta)=\log(1+\delta)$.
Hereafter, we omit the dependence of $c$ on $\delta$ and $\beta$ for
notational convenience.

##### Identification.

Demeaning the error term is required to identify all parameters. Let us assume
the exogeneity condition $E[X_{i}\overline{\upsilon}_{i}]=0$, where
$\overline{\upsilon}_{i}=\upsilon_{i}-c$ denotes the centered error term of
the linearized model. This condition yields the set of $k+1$ equations

$\displaystyle
E\left[X_{i}\left(\log(Y_{i}+\delta\exp(X_{i}^{\prime}\beta))-c\right)\right]=E\left[X_{i}X_{i}^{\prime}\right]\beta,$
(24)

with $k+2$ unknowns. This system identifies $\beta$ only if $c$ is known.
Fortunately, the model in (1) provides the additional restriction necessary
for identification. Let us write
$X_{i}^{\prime}\beta=\beta^{1}+X_{i}^{r^{\prime}}\beta^{r}$, where $\beta^{1}$
is the constant term and the other term represents the non-deterministic part.
We rewrite (1) into

$Y_{i}=\exp(\beta^{1}+X_{i}^{r^{\prime}}\beta^{r})U_{i}=\exp(\beta^{1})\exp(X_{i}^{r^{\prime}}\beta^{r})U_{i}.$
(25)

Rearranging, taking expectations and applying the log function gives the
following expression for the intercept given the other parameters

$\displaystyle\beta^{1}_{\beta}=\log(E[Y_{i}\exp(-X_{i}^{r^{\prime}}\beta^{r})]).$
(26)

Therefore, the parameters are identified and $c$ can be written as191919In our
practical implementation, we solve the identification problem by using the
consistent estimator defined for any $\phi$ as
$\hat{c}(\phi)=\frac{1}{n}\sum_{i=1}^{n}\log(Y_{i}+\delta\exp(\tilde{\phi}^{1}_{\phi}+X_{i}^{r^{\prime}}\phi^{r}))-\frac{1}{n}\sum_{i=1}^{n}(\tilde{\phi}^{1}_{\phi}+X_{i}^{r^{\prime}}\phi^{r})$,
where the constant parameter estimate is replaced by the estimator
$\tilde{\phi}^{1}_{\phi}=\log(n^{-1}\sum_{i=1}^{n}Y_{i}\exp(-X_{i}^{r^{\prime}}\phi^{r}))$.

$\displaystyle c(\delta,\beta)=$ $\displaystyle
E[\log(Y_{i}+\delta\exp(\beta^{1}_{\beta}+X_{i}^{r^{\prime}}\beta^{r}))-\beta^{1}_{\beta}-X_{i}^{r^{\prime}}\beta^{r}].$
(27)

Remark that the estimation of $c$ from the unconditional restriction $E(U)=1$
bears similarities with the smearing estimator (Duan, 1983). Indeed, our
approach integrates this idea directly in the estimation procedure. Doing so
allows avoiding the issue that $\exp(X^{\prime}\hat{\beta})$ is not a
consistent estimate of $E(Y|X)$, which arises if $\beta$ is estimated with OLS
in the log-linear model (Manning, 1998). Therefore, the iOLS estimator of
$\beta$ does correspond to the semi-elasticity (or elasticity if $X$ is
replaced by $\log(X)$) of $Y$ with respect to $X$ around the mean.

##### Asymptotic properties.

We establish the asymptotic properties of iOLSδ in the following theorem.

######  Theorem 1 (Consistency and Normality of iOLSδ)

Under Assumptions 1, 2, 3, and suitable regularity conditions, the iOLS
estimator using $c(\delta,\beta)$ defined in (27) is consistent and achieves
the parametric rate of convergence $n^{-1/2}$ for any $\delta\in(0,+\infty)$.
Formally, we have $n^{1/2}|\hat{\beta}_{t(n)}-\beta|=O_{p}(1)$ as $n\to\infty$
for any $t(n)\geq-\frac{1}{2}\log(n)/\log(\kappa)$, where $\kappa\in[0,1)$ is
the modulus of the associated contraction mapping from $\mathbb{R}^{K}$ to
$\mathbb{R}^{K}$. In addition, iOLSδ is asymptotically normally distributed
such that
$\sqrt{n}\left(\hat{\beta}_{t(n)}-\beta\right)\overset{d}{\to}\mathcal{N}(0,\Omega)$,
as $n\to\infty$, where the covariance can be consistently estimated using

$\hat{\Omega}=\left(\frac{1}{n}X^{\prime}(I-W)X\right)^{-1}\hat{\Sigma}_{0}\left(\frac{1}{n}X^{\prime}(I-W)X\right)^{-1},$
(28)

where $W$ is a diagonal weighting matrix with diagonal elements
$\frac{\delta}{\delta+U_{i}}\in[0,1)$, and $\hat{\Sigma}_{0}$ is a (robust)
consistent estimator of the covariance of $X_{i}(\log(1+U_{i})-c)$ across
observations. $\Omega$ hence corresponds to the asymptotic covariance of the
OLS estimator in the last iteration up to reweighting using $W$.

This asymptotic result guarantees root-n consistent estimates and, for any
fixed $n$, the 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. $\kappa$ depends on the DGP and is decreasing with
$\delta$.

Note that very small values of $\delta$ might prevent the algorithm from
converging in finite time. This occurs because a very small $\delta$ may imply
a $\kappa$ very close to 1, hence a very slow numerical convergence. However,
our algorithm checks this condition and recommends a larger $\delta$ whenever
necessary.202020Our algorithm includes a safety check to ensure that $\kappa$
is sufficiently smaller than 1. In practice, it is estimated as the median
across estimates obtained at each iteration by
$\hat{\kappa}_{t+1}=|\beta_{t+1}-\beta_{t}|/|\beta_{t}-\beta_{t-1}|$.

Remark that the initial estimate $\hat{\beta}_{0}$ can be an important
determinant of the number of iterations before convergence. More iterations
will be needed as it is further away from the fixed-point. However, this
initial estimate only needs to be within the parameter space and does not need
to be consistent for the above theorem to hold.

The asymptotic distributions of iOLSδ and of OLS in the last iteration (once
the estimator has converged) are similar. Although the standard errors of the
latter are incorrect for iOLSδ, a reweighting of the corresponding covariance
matrix using simple algebra is sufficient and allows to use any HAC-robust
covariance estimator.212121A simple approximation of the standard errors for
iOLSδ consists in multiplying those of the last step OLS by a factor
$1+\delta$.

####  3.3.1 Poisson regression as iOLS

We have already shown that iOLSδ approaches the multiplicative Poisson model
for an arbitrarily large $\delta$. Nevertheless, it is possible to use the
iOLS framework to (exactly) estimate both the multiplicative and additive
Poisson models. Doing so can be useful in some settings, either to avoid
possible convergence issues or for speeding up computations, as will be
illustrated in simulations with fixed-effects and endogenous regressors. We
refer to those estimators as iOLSMP and iOLSAP, respectively, and replace
Assumption 3 by

######  Assumption 4 (Poisson condition)

The error term $U_{i}$ admits the weak exogeneity restriction
$E[X_{i}^{\prime}(U_{i}-1)]=0$, which implies $E[U_{i}]=1$.

Under this new assumption, the iOLS procedure yields the multiplicative
Poisson estimates if (27) is changed to

$c_{i}(\delta,\beta)=\log(\delta+Y_{i}\exp(-X_{i}^{\prime}\beta))-\frac{1}{1+\delta}(Y_{i}\exp(-X_{i}^{\prime}\beta)-1),$
(29)

and the PPML estimates if changed to

$c_{i}(\delta,\beta)=\log(\delta+Y_{i}\exp(-X_{i}^{\prime}\beta))-\frac{1}{1+\delta}(Y_{i}-\exp(X_{i}^{\prime}\beta)).$
(30)

These estimators are studied in Appendix B.1, where Theorem 3 shows their
consistency and asymptotic normality. These results show that our approach is
flexible with respect to the choice of both the identifying restriction and
objective criterion without significant consequences in large samples, except
for minor modifications to the covariance matrix. Other extensions, including
i2SLS and high-dimensional fixed effects, are presented in Appendix B.

##  4 Specification testing and model selection

Empirical researchers facing the log of zero usually compare several
estimators to gauge the sensitivity of their results. Yet, each estimator is
only valid under specific unverifiable identifying assumptions. The latter can
nonetheless be systematically investigated by evaluating the model’s external
validity with respect to the observed patterns of zeros.

The tests developed in this section offer an opportunity for an ex-post
evaluation of the identifying restrictions used for moment-based estimators.
Our tests are specification tests used to evaluate the validity of conditional
moment restrictions, like $E(U_{i}|X_{i})=1$ for Poisson models.222222Santos
Silva, Tenreyro and Windmeijer (2015) proposed a radically different approach
based on non-nested hypothesis tests (Davidson and MacKinnon, 1981) which
consist in testing two competing models against each other. They are, as such,
similar to the RESET test of Ramsey (1969) for linear regression and its
application for Poisson models by Wooldridge (1997) and Santos Silva and
Tenreyro (2006). Our approach is, however, tailored to addressing the log of
zero by focusing on the conditional probability of observing $Y=0$. It also
provides a much more powerful test of the conditional restrictions in this
context as will be shown in the simulations.

A common limit of these tests, including RESET tests, is their focus on the
conditional moment restrictions (e.g. $E(U_{i}|X_{i})=1$) rather than the
unconditional restrictions (e.g. $E((U_{i}-1)X_{i})=0$). The former is a
_sufficient_ condition whereas the latter is a _necessary_ condition for
consistency. We argue that statistical evidence against a sufficient condition
is still valuable information about the associated model. The main issue is,
however, that a rejection of the sufficient condition is not evidence against
the necessary condition, which remains unverifiable. Bearing these limits in
mind, we proceed to present our methods.

###  4.1 Specification testing

Let us first look at the implicit assumption made by Poisson models about the
zeros.232323Appendix B.8 details how these tests can be implemented in the
endogenous setting. A related approach will be applied for other restrictions
including iOLSδ. Noting that a zero can only be observed if $U_{i}=0$, the
Poisson restriction $E(U_{i}|X_{i})=1$ can be decomposed into

$E[U_{i}|X_{i}]=E[U_{i}|X_{i},U_{i}>0]Pr(U_{i}>0|X_{i})=E(U_{i}),$ (31)

since $E[U_{i}|X_{i},U_{i}=0]=0$. There are only two possibilities for this
condition to hold true. First, both $E[U_{i}|X_{i},U_{i}>0]$ and
$Pr(U_{i}>0|X_{i})$ vary with $X_{i}$ in such a way that the condition holds.
It happens for example if $U_{i}$ is conditionally Poisson. Second, this
condition also holds if, instead, $E[U_{i}|X_{i},U_{i}>0]$ and
$Pr(U_{i}>0|X_{i})$ are constant. The former is an exogeneity restriction
between $X_{i}$ and $U_{i}$, conditional of the error being positive, which
assumes away any selection bias. The latter implies that discarding zeros or
not before estimation is irrelevant for identification.

This equation characterizes the implicit relation between zeros and non-zero
observations given by

$E[U_{i}|X_{i},U_{i}>0]=\frac{E(U_{i}),}{Pr(U_{i}>0|X_{i})}.$ (32)

It means that the conditional error term for non-zero observations is
inversely proportional to the conditional probability of having a non-zero
observation.242424Heckman’s correction model enforces a comparable conditional
moment restriction:
$E[\log(U_{i})|X_{i},U_{i}>0]=\frac{\lambda\phi(-X_{i}^{\prime}\gamma)}{Pr(U_{i}>0|X_{i})}$.
More generally, moment-based methods typically make implicit assumptions about
selection, whereas sample-selection models enforce explicit restrictions. We
propose to investigate whether this implication matches what is observed in
the data. To do so, we develop a test to assess whether the residuals implied
by the chosen model satisfy this relationship where the conditional
probability is estimated outside the model.252525Testing whether the
probability of observing a zero depends on any $X_{i}$ is done by
investigating the statistical significance of the $X_{i}$’s in a conditional
probability model like the logit. This approach amounts to evaluating the null
hypothesis

$H_{0}:E[U_{i}|X_{i},U_{i}>0]=\frac{1}{Pr(U_{i}>0|X_{i})}.$ (33)

Under this null hypothesis, the error term satisfies the following equation

$U_{i}=\lambda Pr(U_{i}>0|X_{i})^{-1}+\nu_{i}$ (34)

with $\lambda=1$ and $E[\nu_{i}|U_{i}>0,X_{i}]=0$. Therefore, one could
evaluate $H_{0}$ by testing $\lambda=1$ after regressing $U_{i}$ onto
$E[U]Pr(U_{i}>0|X_{i})^{-1}$ for non-zero observations. This test would not be
feasible in practice because none of those variables are observable. Instead,
we propose the following approach:

  1. 1.

Estimate of $Pr(U_{i}>0|X_{i})$ denoted $\hat{P}(X)$.

  2. 2.

Compute Poisson estimates $\hat{\beta}$, e.g. using PPML.

  3. 3.

Recover the residuals $\hat{U_{i}}=Y_{i}\exp(-X^{\prime}\hat{\beta})$;

  4. 4.

Define $\hat{W}_{i}=\hat{P}(X_{i})^{-1}$ and estimate the regression
$\hat{U_{i}}=\lambda\hat{W}_{i}+\nu_{i}$, with the $n_{1}>0$ non-zero
observations to obtain
$\hat{\lambda}=\left(n_{1}^{-1}\sum_{i=1}^{n_{1}}\hat{W}_{i}^{2}\right)^{-1}\left(\sum_{j=1}^{n_{1}}\hat{W}_{j}\hat{U}_{j}\right)$.

  5. 5.

Reject or not $H_{0}$ using a test statistic, such as the t-stat

$t=\frac{\hat{\lambda}-1}{\hat{\sigma}_{\lambda}},$ (35)

where $\hat{\sigma}_{\lambda}$ denotes an estimate of the standard deviation
of $\hat{\lambda}$.

The estimation errors from the first and second steps cannot be neglected. We
study the asymptotic properties of the test statistic (35) in the next theorem
using the following assumptions.

######  Assumption 5 (Consistent first-step estimators)

First-step estimators are consistent under $H_{0}$, so that
$\hat{P}(X_{i})\overset{p}{\to}Pr(U_{i}>0|X_{i})$ and
$\hat{\beta}\overset{p}{\to}\beta$ with $n\to\infty$.

######  Assumption 6 (Conditional independence of first-step errors)

Estimation errors from first-step estimates, denoted
$\varepsilon_{i}^{WU}=\hat{W}_{i}\hat{U}_{i}-W_{i}U_{i}$, and the error term
in the second-step, denoted $W_{i}\nu_{i}$, are independent conditionally on
$U_{i}>0$.

######  Theorem 2 (Asymptotic properties of the test statistic)

Under $H_{0}$, Assumptions 1, 5, 6 and some regularity conditions, we have

$\frac{\sqrt{n_{1}}(\hat{\lambda}-1)}{\hat{\sigma}_{\lambda}}\overset{d}{\to}N(0,1),$
(36)

as $n_{1}\to\infty$, and $\hat{\sigma}_{\lambda}^{2}$ is a consistent
estimator of
$\sigma_{\lambda}^{2}=E(W_{i}^{2})^{-2}E(W_{i}^{2}\nu_{i}^{2})+E(W_{i}^{2})^{-2}E(\varepsilon_{i}^{WU^{2}})$,
the asymptotic variance of $\hat{\lambda}$. In addition, the t-stat diverges
under any alternative as long as $E[\nu_{i}|U_{i}>0,X_{i}]\neq 0$.

The practical implementation of this test requires two elements. First,
computing the t-stat requires to estimate $\sigma_{\lambda}^{2}$ which must
account for the additional noise introduced by first-step estimates.262626The
main difficulty in deriving a closed-form expression for
$\hat{\sigma}_{\lambda}$ is to account for the correlation between
$\hat{P}(X_{i})$ and $\hat{\beta}^{iOLS_{MP}}$ if separately estimated. We
calculate the t-stat after estimating $\hat{\lambda}$’s standard error with
the pairs bootstrap. Second, the test necessitates a consistent estimate of
the conditional probability function $P(U>0|\cdot)$. Specifying a parametric
model, like the logit, probit or even (ex-post bounded) linear probability
model, provides a simple option. However, the misspecification of
$P(U>0|\cdot)$ may distort the test’s size and performance. A nonparametric or
semiparametric estimate of the conditional probability should hence provide a
preferable option. Nevertheless, nonparametric and semiparametric estimates
are subject to small-sample bias due to regularization which can lead to
similar distortions. In our implementation, we use a k-nearest neighbors (kNN)
algorithm (Hastie, Tibshirani and Friedman, 2009). Although consistent,
nonparametric estimators exhibit poor small-sample behaviors at the support’s
boundaries (Sricharan, Raich and Hero, 2010). We address this issue by
truncating (trimming) observations associated with predicted probabilities
outside the 5% and 95%. Another issue related to nonparametric estimation is
the choice of the regularization parameter, i.e. the number of neighbors in
kNN. Selecting a small number of neighbors will lead to non-smooth conditional
probability functions, hence a noisy estimate of $P(X)$. It is hence
preferable to have a reasonably large number of neighbors. The optimal choice
of this parameter for inference is beyond the scope of this paper.

The same logic applies to iOLSδ. The null hypothesis is changed to

$H_{0}:E[\log(\delta+U_{i})|X_{i},U_{i}>0]=\frac{c(\delta,\beta^{iOLS})-\log(\delta)(1-Pr(U_{i}>0|X_{i}))}{Pr(U_{i}>0|X_{i})},$
(37)

and the corresponding regression given by
$\log(\delta+\hat{U_{i}})-\log(\delta)=\lambda W_{i}+\nu_{i},$ for strictly
positive errors only, where
$\hat{U_{i}}=Y_{i}\exp(-X^{\prime}\hat{\beta}^{iOLS_{\delta}})$ and
$W_{i}=(\hat{c}-\log(\delta))\hat{P}(X_{i})^{-1}$ based on $\hat{c}$ obtained
from iOLSδ. The rest of the testing procedure is unchanged.

Remark that a similar procedure could be applied to check for which value of
$\Delta$ the popular fix estimator in (5) is coherent with the observed
patterns of zeros. Unfortunately, choosing $\Delta$ using such procedure does
not address the endogeneity bias introduced by the $\log(\Delta+Y)$
transformation.

###  4.2 Data-driven model selection

The wide variety of available methods for addressing the log of zero forces
researchers to use a selection procedure. We propose to do so in a data-driven
way based on the above results.

First, iOLSδ for any $\delta\in(0,\infty)$ is based on the moment condition
$E[\log(\delta+U)|X]=c(\delta)$ which depends on the choice of $\delta$. If
the “true” DGP belongs to this family of models, then there exists a $\delta$
for which $H_{0}$ will not be rejected as $n\to\infty$. Similarly, the data
may better fit a competing estimator such as PPML or IHS. The tests might
however not be sufficiently discriminatory to pinpoint a single model in small
samples.

We propose to select the model for which $\hat{\lambda}$ is the closest to 1,
among all models (iOLSδ, PPML, and iOLSMP) which are not rejected by the
tests. This heuristic is aimed at selecting the model with the least deviation
between the implied and observed patterns of zeros, i.e. the “most plausible”
model. This selection rule is asymptotically valid because only the true model
will be associated with $\lambda=1$ as $n\to\infty$. In the event that all the
considered models are rejected by the tests, we believe that mixture models
should be preferred.

##  5 Simulations

Let us specify the model as
$Y_{i}=\exp(\beta_{0}+\beta_{1}X_{1i}+\beta_{2}X_{2i})U_{i}$, where
$\beta_{0}=\beta_{1}=\beta_{2}=1$, $U_{i}=\exp(\varepsilon_{i})\xi_{i}$ with
$\xi_{i}=0$ or 1, and
$Pr(\xi_{i}=0|X_{i})=P(X_{i})=\frac{1}{1+\exp(\gamma_{0}+\gamma_{1}X_{1i}+\gamma_{2}X_{2i})}$,
with $\gamma_{0}=-0.4$, $\gamma_{1}=0.4$ and $\gamma_{2}=-0.4$. We consider
two DGPs:272727Appendix C presents additional simulations, including under the
exact $iOLS_{\delta}$-restriction and with endogenous regressors and fixed-
effects.

  * •

DGP 1: $E[U_{i}|X_{i}]=1$ (Poisson). This DGP is aimed at comparing our
approach to Poisson models. We assume that $(X_{1i},X_{2i})^{\prime}$ is
bivariate normal with mean 1, variance $\sigma_{X1}^{2}=\sigma_{X2}^{2}=1$ and
covariance $\sigma_{X1X2}=-0.3$. We further assume that $\varepsilon_{i}$
conditional on $\xi_{i}=1$ is Gaussian with mean $-\log(P(X_{i}))-1/2$ and
variance 1 so that $\exp(\varepsilon_{i})$ is log-normal with conditional mean
$1/P(X_{i})$. Recall that $E(\exp(\varepsilon_{i})|X_{i},\xi_{i}=1)=1/P(X)$ is
implied by the (conditional) Poisson restriction. Note that the error $U_{i}$
is heteroskedastic in this DGP.

  * •

DGP 2: $E[\varepsilon_{i}|X]=0$ (Log-linear). This DGP is useful to illustrate
the flexibility of iOLSδ. We assume that $(X_{1i},X_{2i})^{\prime}$ is like in
DGP 1. We further assume that $\varepsilon$ is gaussian with mean
$\xi_{i}P(X)^{-1}-(1-\xi_{i})(1-P(X))^{-1}$, so that
$E[\varepsilon_{i}|X,\xi_{i}=1]=1/P(X)$. Remark that this condition is the
log-scale counterpart of $E(\exp(\varepsilon_{i})|X,\xi_{i}=1)=1/P(X)$ implied
by Poisson models. The variance parameter is fixed at 0.5. Note that $E(U)$ is
different from 1 and the error $U_{i}$ is also heteroskedastic in this DGP.

Remark that none of these DGPs assume ideal conditions for iOLSδ. We simulate
10,000 times each DGP, for two sample sizes ($n=1,000$ and $n=10,000$), and
report the mean and standard deviations for:282828To save space the results
for OLS after discarding zero observations, the IHS transform, and iOLSAP –
which is equivalent to PPML – are in Table C.1. (1) iOLSδ (best) corresponds
to the “oracle estimator”, i.e. where $\delta$ is chosen so that the mean-
squared error of $\hat{\beta}_{1}$ and $\hat{\beta}_{2}$ is minimized; (2)
iOLSδ with automatic selection of $\delta$, where the conditional probability
is specified as a logit model;292929Misspecifying the probability function
does not affect much the results (Figure C.1). (3) iOLSMP, the multiplicative
Poisson model; (5) PPML, the additive Poisson model estimated using IRLS; and
(6) the PF with $\Delta=1$.303030All simulations were performed in MATLAB
2021b with a 3.6GHz 10-Core Intel Core i9 processor and 32 GB 2667 MHz DDR4
memory.

##### Bias and variance.

The mean iOLSδ estimates along with 95% confidence intervals are shown as
functions of the hyper-parameter $\delta$ in Figure 2 for both DGPs. Figure
2(a) confirms that iOLSδ delivers correct estimates under the Poisson DGP
provided that $\delta$ is sufficiently large (DGP 1). Figure 2(b) illustrates
that the best value for $\delta$ is small in the log-linear model (DGP 2).

(a) DGP1 (Poisson)

(b) DGP2 (Log-linear)

Figure 2: iOLS estimates as functions of $\delta$ (n=10,000)

Table 1 reports the results for DGP 1 (Poisson). We first observe that PPML
and iOLSMP are consistent. iOLSMP yields more precise estimates in comparison.
Second, iOLSδ delivers good estimates with a small bias and a variance
comparable to that of iOLSMP. The automatic selection of $\delta$ introduces
some noise in the smaller sample but this effect becomes negligible in the
larger sample. Its bias is however slightly larger than for the best iOLSδ
estimates in both samples. Note that the standard errors of the (best) iOLSδ
estimate illustrates the root-n consistency result: multiplying the sample
size by 10 yields $\sqrt{10}\approx 3$ times smaller standard errors. Finally,
the PF exhibits large biases although relatively small errors.

Table 1:  Simulations: DGP 1 (Poisson)

|  | $n=1000$ |  |  | $n=10,000$ |
---|---|---|---|---|---|---
Estim. | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$
$iOLS_{\delta}$ (best) | $0.99$ | $1.02$ | $0.99$ | $1.00$ | $1.01$ | $0.99$
| $(0.16)$ | $(0.09)$ | $(0.09)$ | $(0.05)$ | $(0.03)$ | $(0.03)$
$iOLS_{\delta}$ (auto) | $1.04$ | $1.13$ | $0.87$ | $1.00$ | $1.03$ | $0.97$
| $(0.26)$ | $(0.20)$ | $(0.20)$ | $(0.05)$ | $(0.03)$ | $(0.03)$
$iOLS_{MP}$ | $0.99$ | $1.01$ | $0.99$ | $1.00$ | $1.00$ | $1.00$
| $(0.16)$ | $(0.09)$ | $(0.09)$ | $(0.05)$ | $(0.03)$ | $(0.03)$
$PPML$ | $1.01$ | $0.99$ | $0.98$ | $1.00$ | $1.00$ | $1.00$
| $(0.46)$ | $(0.16)$ | $(0.20)$ | $(0.18)$ | $(0.06)$ | $(0.08)$
| $(0.11)$ | $(0.07)$ | $(0.07)$ | $(0.03)$ | $(0.02)$ | $(0.02)$
$PF$ | $0.65$ | $0.61$ | $0.15$ | $0.65$ | $0.61$ | $0.15$
| $(0.09)$ | $(0.06)$ | $(0.06)$ | $(0.03)$ | $(0.02)$ | $(0.02)$

  * •

Notes: This table shows the mean and standard errors (in parentheses) of
parameter estimates across 10,000 simulations based on DGP1.

Table 2 reports the results for DGP 2 (Log-linear). All estimators have large
biases except iOLSδ. Poisson estimators perform poorly in this context. In
particular, PPML diverges to extremely large values in many
simulations.313131The same result is observed in Table C.1.

Table 2:  Simulations: DGP 2 (Log-linear)

|  | $n=1000$ |  |  | $n=10,000$ |
---|---|---|---|---|---|---
Estim. | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$
$iOLS_{\delta}$ (best) | $4.52$ | $1.15$ | $0.85$ | $5.36$ | $1.01$ | $0.99$
| $(0.81)$ | $(0.22)$ | $(0.22)$ | $(1.16)$ | $(0.17)$ | $(0.17)$
$iOLS_{\delta}$ (auto) | $4.45$ | $1.06$ | $0.94$ | $5.64$ | $1.05$ | $0.95$
| $(1.34)$ | $(0.17)$ | $(0.17)$ | $(1.68)$ | $(0.11)$ | $(0.10)$
$iOLS_{MP}$ | $3.20$ | $0.44$ | $1.57$ | $3.42$ | $0.30$ | $1.70$
| $(0.36)$ | $(0.19)$ | $(0.20)$ | $(0.51)$ | $(0.22)$ | $(0.23)$
| $(4.50)$ | $(1.17)$ | $(1.32)$ | $(10.21)$ | $(1.95)$ | $(2.25)$
$PPML$ | $-Inf$ | $-Inf$ | $Inf$ | $Inf$ | $-Inf$ | $Inf$
| $(Inf)$ | $(Inf)$ | $(Inf)$ | $(Inf)$ | $(Inf)$ | $(Inf)$
$PF$ | $1.43$ | $0.67$ | $0.14$ | $1.43$ | $0.67$ | $0.14$
| $(0.15)$ | $(0.10)$ | $(0.10)$ | $(0.05)$ | $(0.03)$ | $(0.03)$

  * •

Notes: This table shows the mean and standard errors (in parentheses) of
parameter estimates across 10,000 simulations based on DGP 2.

##### Computational gains with endogenous regressors and many fixed-effects.

We briefly describe the results of some additional simulations detailed in
Appendix C. DGP 1 (Poisson) is modified so that both regressors are now
correlated with the error term. In absence of fixed-effects, the i2SLS
estimators perform similarly to the non-linear IV estimator (Mullahy, 1997) in
terms of bias and variance, and provide the best options. The control function
approach of Wooldridge (1997) delivers biased estimates because the endogenous
regressors are non-linear functions of the error term in this design. Additive
Poisson estimators are not consistent because the unobserved heterogeneity
does not enter $Y$ additively.

i2SLS estimators bring important computational gains in settings with fixed-
effects. We consider $T=100$ periods and $n/T$ individuals to model two-way
fixed effects. For $n=1,000$, i.e. 111 regressors, $i2SLS_{MP}$ is computed in
about 0.08 second whereas the non-linear IV estimator in Mullahy (1997) takes
45 seconds. For $n=10,000$, i.e. 201 regressors, $i2SLS_{MP}$ is computed in
about 0.8 second whereas the non-linear IV estimator takes above 500 seconds.

##### Automatic selection and specification tests.

These simulations are also useful to study our testing procedures and
automatic selection rule. The conditional probabilities to have a zero are
logistic in all DGPs. In what follows, we mainly focus on the correct
parametric specification to compute the conditional probability of observing
zero values (logit). We also report and discuss some results when using a
nonparametric approach (kNN). Figure C.1 in Appendix C compares the results
when misspecifying the conditional probability function as a Probit or a
linear probability model. The differences remain small in all cases.

First, we illustrate the automatic model selection in Figure 3. The plain
lines show $\hat{\lambda}$ as a function of $\log(\delta)$, and the dashed
lines represent the (pointwise) 95% confidence intervals across simulations.
The selection rule is to pick $\delta$ such that the parameter $\lambda$ is as
close to 1 as possible. For the Poisson DGP, Figure 3(a) confirms that
$\hat{\lambda}$ tends to be closer to 1 for larger values of $\delta$. For the
log-linear DGP, Figure 3(b) corroborates that $\hat{\lambda}$ tends to be
closer to 1 for small values of $\delta$.

(a) DGP1 (Poisson)

(b) DGP3 (Log-linear)

Figure 3: iOLS’s $\hat{\lambda}$ as functions of $\delta$ (n=10,000)

Second, we investigate the performance of the inference procedure based on the
pairs bootstrap. We perform 1000 simulations of each DGP with $n=10,000$.
Table 3 reports the mean estimates of $\lambda$, its standards errors across
simulations (in parentheses), and the empirical rejection rate of each test
for a nominal size of 5% and 300 bootstrap draws. We show the results for
various $iOLS_{\delta}$ illustrating the results in Figures 2 and 3,
$iOLS_{MP}$ and $PPML$. All tests, except $RESET$, evaluates the null
hypothesis that the underlying model is correct using the approach laid out in
Section 4.1. Finally, the last column shows rejection rates for Ramsey’s RESET
test for PPML including 3 polynomials terms (Santos Silva and Tenreyro, 2006;
Ramsey, 1969).323232https://personal.lse.ac.uk/tenreyro/reset.do.

Table 3:  Simulations: Specification testing (Logit)

DGP |  |  iOLSδ | MP | PPML | RESET
---|---|---|---|---|---
| $\delta$ | 0.01 | 0.15 | 1.9 | 24 | 86 |  |  |
1 | $\lambda$ | $0.98$ | $0.98$ | $0.99$ | $1.00$ | $1.01$ | $1.00$ | $1.01$ |
(se) | $(0.01)$ | $(0.01)$ | $(0.01)$ | $(0.01)$ | $(0.01)$ | $(0.01)$ | $(0.08)$ |
Rej% | $88.8$ | $72.3$ | $26.1$ | $5.3$ | $15.7$ | $5.7$ | $8.4$ | $7.3$
3 | $\lambda$ | $1.01$ | $1.04$ | $1.10$ | $1.22$ | $1.29$ | $1.43$ | $Inf$ |
(se) | $(0.11)$ | $(0.15)$ | $(0.20)$ | $(0.26)$ | $(0.29)$ | $(0.37)$ | $(Inf)$ |
Rej | $5.6$ | $66.6$ | $94.2$ | $99.7$ | $100.0$ | $100.0$ | $0.2$ | $5.9$

  * •

Notes: This table shows the relative rejection frequency of each null
hypothesis for 1,000 simulations for $n=10,000$ and DGP 1 to 3.

In DGP 1, the tests for iOLSδ=24 and iOLSMP reject the null hypothesis only
5.3% and 5.7% of the time which is close to the nominal test size. The test
for PPML and the RESET test are slightly oversized with 8.4% and 7.3%
rejection rates. The other models are rejected with a larger frequency. In DGP
2, iOLSδ is rejected only 5.6% of the time for $\delta=0.01$, all other models
are vastly rejected. PPML estimates behave poorly leading to a non-informative
test under the alternative. Finally, findings reveal that the RESET test lacks
power against the considered alternative. Similar results are obtained using
KNN (Table C.2).

##  6 Application

###  6.1 Michalopoulos and Papaioannou (2013)

We now revisit the analysis of Michalopoulos and Papaioannou (2013), a leading
empirical study where the log of zero had to be addressed. They examine the
relationship between pre-colonial political centralization and contemporary
development in African countries. The latter is proxied using light density at
night at the regional level and used as the response variable through the
“popular fix”: $\log(\text{Y}_{i}+0.01)$. The effect of political
centralization is measured by the coefficient associated with Murdock’s 1967
index of jurisdictional hierarchy.333333Ranging between 0 and 4, it provides
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. They control
for population density, location, and geography, as well as country fixed
effects.343434We focus on columns (4) of their Table 3 which include all
controls.

We first search for the optimal value of $\delta$ using the test and model
selection procedure developed in Section 4. We estimate the conditional
probability function353535We discard the fixed-effects from the probability
model to prevent over-fitting and avoid computational difficulties. The
results are similar with ethnicity by country fixed effects. using the logit
and kNN algorithm respectively.363636For kNN, we use 100 neighbors and trim
observations associated with predicted probabilities below the 5th and above
the 95th (empirical) percentiles. Figure 4(a) displays the results, for a grid
of candidate $\delta$.373737The grid starts at $\exp(-7)\approx 0.001$ with
increments of $\exp(0.5)$ up to $\exp(7)\approx 1000$. The choice of $\delta$
matters in practice. As shown in Figure 4(b), the estimate for jurisdictional
hierarchy, denoted by $\beta$, ranges from negative to positive, as well as
non-significant to significant at the 5% level when using standard errors
clustered at the ethnicity by country level. The selected $\delta$ is large
for both probability models and lead to the same estimates of $\beta$.

(a) $\hat{\lambda}$ as a function of $\delta$

(b) $\hat{\beta}$ as a function of $\delta$

Figure 4: Jurisdictional hierarchy ($\beta$) with iOLSδ and $\lambda$ for
various candidate $\delta$

Table 4 reports competing estimates of the effect of jurisdictional hierarchy
on economic development. Both the author’s specification (PF), multiplicative
Poisson (iOLSMP), and iterated ordinary least squares (iOLSδ) deliver a
positive and statistically significant estimate. In contrast, PPML reports a
negative and non-significant point-estimate. The $\lambda$’s associated with
the author’s specification and PPML are far from one. Following our model
selection procedure, it means that both exhibit poor external validity with
respect to the probability of observing a zero. Conversely, $\lambda$ is close
to one for both iOLSδ≈1000 and iOLSMP, which are statistically the same.

Therefore, the econometric evidence is in favor of a positive and significant
effect of political centralization on economic development, as found by the
authors.383838In subsequent work, Michalopoulos and Papaioannou (2014) study
the link between contemporary political institutions in Africa and economic
development. They find an absence of statistical significance using both the
popular fix, OLS in level, and PPML (Table 6 in their Appendix). This
replication exercise shows that the popular fix, although fundamentally
flawed, can at times lead to qualitatively correct results. However, it still
delivers artificially lower standard errors by treating zero observations
similarly to non-zero observations. We also note that relying solely on PPML
can yield misleading results.

Table 4:  The effect of political centralization

| PF | PPML |  iOLSMP |  iOLSδ=exp(7)
---|---|---|---|---
Jurisdictional Hierarchy ($\beta$) |  0.177∗∗∗ | -0.110 |  0.149∗∗ |  0.149∗∗
| (0.0473) | (0.0812) | (0.0634) | (0.0733)
Pop. & Loc. & Geo. | Yes | Yes | Yes | Yes
Country Fixed-Effects | Yes | Yes | Yes | Yes
$\lambda$ \- kNN | -0.163 | 2.399 | 1.052 | 1.052
| (0.029) | (3.768) | (0.034) | (0.035)
p-value | 0.000 | 0.371 | 0.126 | 0.134
$N$ | 682 | 682 | 682 | 682
∗ $p<0.05$, ∗∗ $p<0.01$, ∗∗∗ $p<0.001$

  * •

Notes: This table displays the coefficient associated with jurisdictional
hierarchy, standard errors clustered at the country level in parenthesis.

###  6.2 Santos Silva and Tenreyro (2006)

We now revisit the gravity model of Santos Silva and Tenreyro (2006). Their
Table 3 reports models of bilateral trade for data covering 136 countries in
1990. The authors advocate for PPML over a log-linear model arguing that the
latter is biased in presence of heteroskedastic errors. For importer (I) and
exporter (X) countries, they control for log(Distance), along with dummies for
contiguity, shared language, colonial ties, and two variables which, when
summed, measure the impact of free trade agreements (FTA) on trade
flows.393939They also control for log(GDP) and log(GDP per capita), along with
a dummies for access to the oceans, remoteness (which measures the access to
other trading partners). The estimates for the associated coefficients are
available from Table D.5 in the Appendix. We denote the latter sum as $\beta$
and focus on its sensitivity to alternative identifying assumptions.

As in the previous illustration, we search for the optimal $\delta$ and
associated $\lambda$, based on both a logit and a kNN estimator. As shown in
Figure 5(a), The latter provides a series of $\lambda$ which do not reject the
proportionality condition $\lambda=1$, and leads to select $\delta=13$ with
$\lambda=1.004$. The measure of trade agreements $\beta$ evolves from non-
significant for small $\delta$’s to statistically significant for larger
values. In general, those estimates are fairly different from those obtained
using PPML.

(a) $\hat{\lambda}$ as a function of $\delta$

(b) $\hat{\beta}$ as a function of $\delta$

Figure 5: Free Trade Agreement ($\beta$) with iOLSδ and $\lambda$ for various
candidate $\delta$

Table 5 reports estimates of the effect of Free Trade Agreements for the
different models.404040Table D.5 in Appendix 6 reports the complete results.
Whilst the popular fix (with $\Delta=1$) gives a very large effect of trade
agreements on trade flows, PPML’s estimate is close to zero and non-
significant. For both, $\lambda$ is far from one, suggesting a lack of
external validity of the model with respect to the pattern of zeros in the
data. In contrast, iOLSMP and iOLSδ deliver economically comparable results
with $\lambda$ close to one. These latter estimates are also close to what the
authors obtain using the Tobit approach.

These results illustrate two important effects. On the one hand, the popular
fix give too much weight to the zeros by treating them like “normal
observations” which artificially inflates the effect of Free Trade Agreements
and reduces its standard errors. On the other hand, PPML under-estimates the
importance of Free Trade Agreements by putting too much weight on observations
associated with large trade values and too little on smaller values,
especially the zeros. In contrast, our approach considers a wide range of
solutions and allows one to focus on the ones which make sense with respect to
the observed zeros.

Table 5:  The effect of Free Trade Agreements

| PF | PPML |  iOLSMP |  iOLSδ=13
---|---|---|---|---
Free Trade Agreement ($\beta$) |  2.028∗∗∗ | 0.074 |  1.135∗∗ |  1.056∗∗∗
| (0.131) | (0.177) | (0.513) | (0.349)
$\lambda$ \- kNN | 0.027 | 0.449 | 1.020 | 1.004
| (0.017) | (0.060) | (0.072) | (0.054)
p-value | 0.000 | 0.000 | 0.777 | 0.934
$N$ | 18,360 | 18,360 | 18,360 | 18,360
∗ $p<0.05$, ∗∗ $p<0.01$, ∗∗∗ $p<0.001$

  * •

Notes: This table displays the coefficient associated with being part of a
free trade agreement, standard errors based on 300 pairs bootstrap are in
parenthesis.

This application provides a clear illustration of the general result that PPML
estimates tend to be driven by observations with large values in $Y$ while
neglecting those with small values. In Table 6, we report the estimated effect
of trade agreements based on alternative models after discarding observations
based on the quantiles of the dependent variable. As expected, the estimates
barely change for PPML, even when keeping only the top 10% of trade flows.
Comparable results hold for the other explanatory variables (Table D.5). In
contrast, all other methods show a decreasing effect when focusing on
observations with larger $Y$. This suggests that Free Trade Agreements have
relatively large effects between countries with small trade flows but not for
the largest trading nations. In conclusion, both the popular fix and PPML may
yield misleading results.

Table 6:  Trade Inequality and Estimates of Free Trade Agreements ($\beta$)

| PF | PPML |  iOLSMP |  iOLSδ=13 | Observations
---|---|---|---|---|---
Full Sample | 2.028 | 0.074 | 1.135 | 1.048 | 18360
| (0.131) | (0.177) | (0.513) | (0.349) |
No Zeros | 0.325 | 0.040 | 0.207 | 0.141 | 9613
| (0.109) | (0.179) | (0.216) | (0.151) |
Top 25% | 0.155 | 0.064 | -0.036 | -0.038 | 4589
| (0.079) | (0.174) | (0.118) | (0.106) |
Top 10% | -0.045 | 0.076 | -0.248 | -0.237 | 1836
| (0.088) | (0.172) | (0.106) | (0.103) |

  * •

Notes: This table displays the coefficient associated with being part of a
free trade agreement, based on different subsamples of the data.

##  7 Conclusion

This paper brought multiple contributions to address a common yet unresolved
issue faced in empirical research: the log of zero. Our estimation procedure
has several advantages, including: 1) computational simplicity, 2) a natural
extension to instrumental variables, 3) robustness to the inclusion of many
fixed effects, and 4) their flexibility to exogeneity restrictions. Our
testing procedures allows verifying the underlying exogeneity restrictions
imposed on the occurrences of zeros. Our replications of leading publications
have shown how these tests can guide empirical research.

Hopefully, empirical researchers are now better equipped to address the log of
zero and justify their chosen method. The methodology developed in this paper
should help find a consensus among practitioners about the best practice to
address the log of zero. There are also many possible extensions, like
regularized models, which we leave for future research.

## References

  * (1)
  * Bellemare and Wichman (2020) Bellemare, M. F., and C. J. Wichman. 2020. “Elasticities and the Inverse Hyperbolic Sine Transformation.” Oxford Bulletin of Economics and Statistics, 82(1): 50–61.
  * Burbidge, Magee and Robb (1988) Burbidge, J. B., L. Magee, and A. L. Robb. 1988. “Alternative Transformations to Handle Extreme Values of the Dependent Variable.” Journal of the American Statistical Association, 83(401): 123–127.
  * Cohn, Liu and Wardlaw (2021) Cohn, Jonathan B., Zack Liu, and Malcolm Wardlaw. 2021. “Count Data in Finance.” Available at SSRN.
  * Correia, Guimarães and Zylkin (2019) Correia, S., P. Guimarães, and T. Zylkin. 2019. “ppmlhdfe: Fast Poisson Estimation with High-Dimensional Fixed Effects.” arXiv e-prints.
  * Davidson and MacKinnon (1981) Davidson, Russell, and James G. MacKinnon. 1981. “Several Tests for Model Specification in the Presence of Alternative Hypotheses.” Econometrica: Journal of the Econometric Society, 781–793.
  * Dominitz and Sherman (2005) Dominitz, Jeff, and Robert P. Sherman. 2005. “Some Convergence Theory for Iterative Estimation Procedures with an Application to Semiparametric Estimation.” Econometric Theory, 21(4): 838–863.
  * Duan (1983) Duan, Naihua. 1983. “Smearing Estimate: a Nonparametric Retransformation Method.” Journal of the American Statistical Association, 78(383): 605–610.
  * Eaton and Tamura (1994) Eaton, Jonathan, and Akiko Tamura. 1994. “Bilateralism and Regionalism in Japanese and U.S. Trade and Direct Foreign Investment Patterns.” Journal of the Japanese and International Economies, 8(4): 478 – 510.
  * Fernández-Val and Weidner (2016) Fernández-Val, Iván, and Martin Weidner. 2016. “Individual and Time Effects in Nonlinear Panel Models with Large N, T.” Journal of Econometrics, 192(1): 291–312.
  * Gourieroux, Monfort and Trognon (1984) Gourieroux, Christian, Alain Monfort, and Alain Trognon. 1984. “Pseudo Maximum Likelihood Methods: Applications to Poisson Models.” Econometrica, 52(3): 701–20.
  * Hastie, Tibshirani and Friedman (2009) Hastie, Trevor, Robert Tibshirani, and Jerome Friedman. 2009. “Prototype Methods and Nearest-Neighbors.” The Elements of Statistical Learning: Data Mining, Inference, and Prediction, 459–483. New York, NY:Springer New York.
  * Head and Mayer (2014) Head, Keith, and Thierry Mayer. 2014. “Gravity Equations: Workhorse,Toolkit, and Cookbook.” In Handbook of International Economics. Vol. 4 of Handbook of International Economics, 131–195. Elsevier.
  * Heckman (1979) Heckman, James J. 1979. “Sample Selection Bias as a Specification Error.” Econometrica: Journal of the Econometric Society, 153–161.
  * MacKinnon and Magee (1990) MacKinnon, James G, and Lonnie Magee. 1990. “Transforming the Dependent Variable in Regression Models.” International Economic Review, 315–339.
  * Manning (1998) Manning, Willard G. 1998. “The Logged Dependent Variable, Heteroscedasticity, and the Retransformation Problem.” Journal of Health Economics, 17(3): 283–295.
  * Manning and Mullahy (2001) Manning, Willard G, and John Mullahy. 2001. “Estimating Log Models: to Transform or not to Transform?” Journal of Health Economics, 20(4): 461–494.
  * Michalopoulos and Papaioannou (2013) Michalopoulos, Stelios, and Elias Papaioannou. 2013. “Pre-Colonial Ethnic Institutions and Contemporary African Development.” Econometrica, 81(1): 113–152.
  * Michalopoulos and Papaioannou (2014) Michalopoulos, Stelios, and Elias Papaioannou. 2014. “National Institutions and Subnational Development in Africa.” The Quarterly Journal of Economics, 129(1): 151–213.
  * Mullahy (1997) Mullahy, John. 1997. “Instrumental-Variable Estimation of Count Data models: Applications to Models of Cigarette Smoking Behavior.” Review of Economics and Statistics, 79(4): 586–593.
  * Ramsey (1969) Ramsey, James Bernard. 1969. “Tests for Specification Errors in Classical Linear Least-Squares Regression Analysis.” Journal of the Royal Statistical Society: Series B (Methodological), 31(2): 350–371.
  * Ravallion (2017) Ravallion, M. 2017. “A Concave Log-Like Transformation Allowing Non-Positive Values.” Economics Letters, 161: 130 – 132.
  * Santos Silva and Tenreyro (2006) Santos Silva, João M.C., and Silvana Tenreyro. 2006. “The Log of Gravity.” The Review of Economics and Statistics, 88(4): 641–658.
  * Santos Silva and Tenreyro (2010) Santos Silva, João M.C., and Silvana Tenreyro. 2010. “On the Existence of the Maximum Likelihood Estimates in Poisson Regression.” Economics Letters, 107(2): 310 – 312.
  * Santos Silva and Tenreyro (2011) Santos Silva, João M.C., and Silvana Tenreyro. 2011. “Poisson: Some Convergence Issues.” The Stata Journal, 11(2): 207 – 212.
  * Santos Silva, Tenreyro and Windmeijer (2015) Santos Silva, João M.C., Silvana Tenreyro, and Frank Windmeijer. 2015. “Testing Competing Models for Non-negative Data with Many Zeros.” Journal of Econometric Methods, 4(1): 29 – 46.
  * Sricharan, Raich and Hero (2010) Sricharan, Kumar, Raviv Raich, and Alfred O Hero. 2010. “Boundary Compensated k-NN Graphs.” 2010 IEEE International Workshop on Machine Learning for Signal Processing, 277–282.
  * Weidner and Zylkin (2021) Weidner, Martin, and Thomas Zylkin. 2021. “Bias and Consistency in Three-Way Gravity Models.” Journal of International Economics, 132: 103513.
  * Whitney and McFadden (1994) Whitney, K. Newey, and Daniel McFadden. 1994. “Large Sample Estimation and Hypothesis Testing.” Handbook of Econometrics: Chapter 36, 2111–2245.
  * Winkelmann (2008) Winkelmann, Rainer. 2008. Econometric Analysis of Count Data. Springer Science & Business Media.
  * Wooldridge (1997) Wooldridge, Jeffrey M. 1997. “Quasi-Likelihood Methods for Count Data.” Handbook of Applied Econometrics, 2: 352–406.
  * Wooldridge (2015) Wooldridge, Jeffrey M. 2015. “Control Function Methods in Applied Econometrics.” Journal of Human Resources, 50(2): 420–445.
  * Young and Young (1975) Young, K. H., and L. Y. Young. 1975. “Estimation of Regressions Involving Logarithmic Transformation of Zero Values in the Dependent Variable.” The American Statistician, 29(3): 118–120.

Appendix

######  Proof 1 (Proof of Theorem 1)

Recall that the parameter $\beta\in\mathbb{R}^{K}$ is characterized by the
fixed-point equation

$\beta=E[X_{i}X_{i}^{\prime}]^{-1}E\left[X_{i}\tilde{Y}_{i}(\beta)\right],$
(38)

where
$\tilde{Y}_{i}(\beta)=\log(Y_{i}+\delta\exp(X_{i}^{\prime}\beta))-c(\beta,\delta)$
is the transformed dependent variable. To simplify exposition, we focus on
$\delta=1$ in our calculations. The mapping from $\mathbb{R}^{K}$ to
$\mathbb{R}^{K}$ which characterizes the parameter is hence defined
$\forall\phi\in\mathbb{R}^{K}$ as

$M(\phi)=E[X_{i}X_{i}^{\prime}]^{-1}E\left[X_{i}\tilde{Y}_{i}(\phi)\right].$
(39)

The sample counterpart of this mapping is given by

$\hat{M}_{n}(\phi)=\left[X^{\prime}X\right]^{-1}X^{\prime}\hat{\tilde{Y}}_{i}(\phi),$
(40)

where
$\hat{\tilde{Y}}_{i}(\phi)=\log(Y_{i}+\exp(X_{i}^{\prime}\phi))-\hat{c}(\phi)$,
with
$\hat{c}(\phi)=\frac{1}{n}\sum_{i=1}^{n}\log(Y_{i}+\exp(\hat{\phi}_{1}(\phi)-\phi_{1}+X_{i}^{\prime}\phi))-\log(\frac{1}{n}\sum_{i=1}^{n}(\hat{\phi}_{1}(\phi)-\phi_{1}+X_{i}^{\prime}\phi))$
for
$\hat{\phi}_{1}(\phi)=\log(n^{-1}\sum_{i=1}^{n}Y_{i}\exp(-X_{i}\phi+\phi_{1}))$

Our proof follows Dominitz and Sherman (2005), hereafter denoted DS, who
develop a convergence theory for iterative estimators. Following DS, the
convergence of iOLS requires that $M(\cdot)$ and $\hat{M}_{n}(\cdot)$ be
contraction mappings, asymptotically.414141The reader is referred to DS for a
formal definition of an asymptotic contraction mapping.

In order to show the convergence result
$n^{1/2}|\hat{\beta}_{t(n)}-\beta|=O_{p}(1)$ as $n\to\infty$ by applying
Theorem 1 in DS, we need to show that the following conditions hold:

  1. (i)

$\left\\{\hat{M}_{n}(\cdot):n\geq 1,\omega\in\mathcal{S}\right\\}$ is an
asymptotic contraction mapping on $(B_{0},E_{K})$, where $\mathcal{S}$ is a
sample space, $E_{K}$ is the Euclidean metric on $\mathcal{R}^{K}$ and $B_{0}$
is the closed ball centered at $\beta_{0}$ of radius
$|\hat{\beta}_{0}-\beta|$;424242Note that DS’s condition (i) is about
$M(\cdot)$ and not $\hat{M}_{n}(\cdot)$, with $M(\cdot)$ being non-stochastic
in our setting. However those conditions imply each other under conditions
(iii) and (iv) by applying their Lemma 3 with trivial modifications.

  2. (ii)

$n^{1/2}|\beta_{t(n)}-\beta|=O_{p}(1)$ as $n\to\infty$;

  3. (iii)

$n^{1/2}\sup_{\phi\in B_{0}}|\hat{M}_{n}(\phi)-M(\phi)|=O_{p}(1)$ as
$n\to\infty$; and

  4. (iv)

$\sup_{\phi\in B_{0}}||\hat{V}_{n}(\phi)-V(\phi)||=o_{p}(1)$ as $n\to\infty$.

##### Regularity conditions.

Our proofs rely on the sufficient regularity conditions listed below, in
particular for showing the uniform convergence in conditions (iii) and (iv):
(1) $E\left[X_{i}\right]<\infty$; (2) $V\left[X_{i}\right]<\infty$; (3)
$E\left[X_{i}\log(Y_{i}+\exp(X_{i}^{\prime}\phi))\right]<\infty,\text{
}\forall\phi\in B_{0}$; (4)
$V\left[X_{i}\log(Y_{i}+\exp(X_{i}^{\prime}\phi))\right]<\infty,\text{
}\forall\phi\in B_{0}$; (5) $c(\phi)<\infty,\text{ }\forall\phi\in B_{0}$; (6)
$V\left[\hat{c}(\phi)\right]<\infty,\text{ }\forall\phi\in B_{0}$; (7)
$E\left[X_{i}\frac{\exp(X_{i}^{\prime}\phi)}{Y_{i}+\exp(X_{i}^{\prime}\phi)}X_{i}^{\prime}\right]<\infty;\text{
}\forall\phi\in B_{0}$; (8)
$V\left[X_{i}\frac{\exp(X_{i}^{\prime}\phi)}{Y_{i}+\exp(X_{i}^{\prime}\phi)}X_{i}^{\prime}\right]<\infty,\text{
}\forall\phi\in B_{0}$; (9) $\nabla_{\phi}c(\phi)<\infty,\text{
}\forall\phi\in B_{0}$; and (10)
$V\left[\nabla_{\phi}\hat{c}(\phi)\right]<\infty,\text{ }\forall\phi\in
B_{0}$.

##### Condition (i).

Let us adapt the proof of Lemma 5 in DS as follows. The first step is to
consider that $X$ is prewhitened so that $X^{\prime}X=nI_{k}$. This assumption
is useful to establish the local contraction mapping property without loss of
generality. From a multivariate Taylor expansion argument, DS show that
condition (i) boils down to showing that the largest eigenvalue of
$\nabla_{\phi}\hat{M}_{n}(\beta)=\hat{V}_{n}(\beta)$ is strictly less than
unity as $n\to\infty$. Note that we have

$\displaystyle\hat{V}_{n}(\phi)$
$\displaystyle=[X^{\prime}X]^{-1}X^{\prime}\nabla_{\phi}\hat{\tilde{Y}}(\phi)$
(41) $\displaystyle=n^{-1}X^{\prime}\nabla_{\phi}\hat{\tilde{Y}}(\phi),$

where the second equality uses prewhitening and
$\nabla_{\phi}\hat{\tilde{Y}}_{i}(\phi)$ has element $(i,k)$ defined as

$\displaystyle\left[\nabla_{\phi}\hat{\tilde{Y}}(\phi)\right]_{i,k}$
$\displaystyle=\frac{\exp(X_{i}^{\prime}\phi)X_{ki}}{Y_{i}+\exp(X_{i}^{\prime}\phi)}-\frac{\partial\hat{c}(\phi)}{\partial\phi_{k}}.$
(42)

Let us denote $X_{1i}=1$, for all $i$ as the constant. By prewhitening, we
have $\sum_{j=1}^{n}X_{1j}=n$ and $\sum_{j=1}^{n}X_{kj}=0$ for $k>1$.

$\displaystyle\frac{\partial\hat{c}(\phi)}{\partial\phi_{k}}=n^{-1}\sum_{i=1}^{n}\frac{\exp(X_{i}^{r^{\prime}}\phi^{r}+\hat{\phi}^{1})(\frac{\partial\hat{\phi}^{1}}{\partial\phi_{k}}+X_{ki})}{Y_{i}+\exp(X_{i}^{r^{\prime}}\phi^{r}+\hat{\phi}^{1})}-n^{-1}\sum_{i=1}^{n}(\frac{\partial\hat{\phi}^{1}}{\partial\phi_{k}}+X_{ki}),$
(43)

for $k>1$ and $\frac{\partial\hat{c}(\phi)}{\partial\phi_{1}}=0$. This
expression simplifies when evaluated at $\phi=\beta$, as shown by

$\displaystyle\frac{\partial\hat{c}(\beta)}{\partial\phi_{k}}=n^{-1}\sum_{i=1}^{n}\frac{X_{ki}}{1+U_{i}}+O_{p}(1),$
(44)

for $k>1$ because
$\hat{\phi}^{1}(\beta)=\log(n^{-1}\sum_{i=1}Y_{i}\exp(-X_{i}^{r^{\prime}}\beta^{r}))=\beta_{1}+\log(n^{-1}\sum_{i=1}U_{i})$,
where $\log(n^{-1}\sum_{i=1}U_{i})=O_{p}(1)$ by iid assumption and
$E[U_{i}]=1$, and $n^{-1}\sum_{i=1}^{n}X_{ki}=0$ by prewhitening. Thus, we
have $\frac{\partial\hat{\phi}_{1}(\beta)}{\partial\phi_{k}}=0$.

Therefore, each element $(k,l)$ of $\hat{V}_{n}(\beta)$ writes

$\displaystyle\left[\hat{V}_{n}(\beta)\right]_{k,l}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\frac{X_{ki}X_{li}}{1+U_{i}}-n^{-1}\sum_{i=1}^{n}X_{ki}n^{-1}\sum_{j=1}^{n}\frac{X_{lj}}{1+U_{j}},$
(45)

for $l>1$ and

$\displaystyle\left[\hat{V}_{n}(\beta)\right]_{k,l}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\frac{X_{ki}}{1+U_{i}},$ (46)

for $l=1$. Remark that for $k=1,\forall l>1$ we have
$\left[V_{n}(\beta)\right]_{1,l}=0$, and for $k=1,l=1$, we have
$\left[\hat{V}_{n}(\beta)\right]_{1,1}=n^{-1}\sum_{i=1}^{n}\frac{1}{1+U_{i}}<1$.
Therefore, the eigenvalue associated with the constant term is strictly below
$1$, and proving the convergence amounts to showing that the largest
eigenvalue of the $(K-1)\times(K-1)$ lower right-hand submatrix of
$\hat{V}_{n}(\beta)$ is strictly less than unity. All elements $(k,l)$ for
$k,l>1$ of this matrix writes

$\displaystyle\left[\hat{V}_{n}(\beta)\right]_{k,l}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\frac{X_{ki}X_{li}}{1+U_{i}}.$ (47)

because of prewhitening. We can write this in matrix form as

$\displaystyle\left[\hat{V}_{n}(\beta)\right]_{k,l>1}$
$\displaystyle=n^{-1}X^{\prime}WX,$ (48)

where $W$ is a diagonal matrix with elements $(i,i)$ acting as weights given
by $\frac{1}{1+U_{i}}\in(0,1]$. Note that those weights become
$\frac{\delta}{\delta+U_{i}}\in[0,1)$ for $\delta\neq 1$. We can thus write
$W=W^{1/2}W^{1/2}$, and rewrite the submatrix of interest as the quadratic
form

$\displaystyle\left[\hat{V}_{n}(\beta)\right]_{k,l>1}$
$\displaystyle=n^{-1}X^{\prime}W^{1/2}W^{1/2}X.$ (49)

Consequently, this matrix is nonnegative definite and so must have all
nonnegative eigenvalues. We can alternatively write the weight matrix
$W=I_{n}-D$, where $D$ is also a diagonal matrix with elements
$\frac{U_{i}}{1+U_{i}}\in[0,1)$, or more generally
$\frac{U_{i}}{\delta+U_{i}}\in[0,1)$. Therefore, we have the alternative
expression

$\displaystyle\left[\hat{V}_{n}(\beta)\right]_{k,l>1}$
$\displaystyle=n^{-1}X^{\prime}(I_{n}-D)X=I_{K-1}-n^{-1}X^{\prime}D^{1/2}D^{1/2}X,$
(50)

where the second term is also a quadratic form. It follows that as
$n\to\infty$, the maximum eigenvalue is equal to

$\displaystyle\max_{|a|=1}a^{\prime}\left[\hat{V}_{n}(\beta)\right]_{k,l>1}a=\max_{|a|=1}1-a^{\prime}X^{\prime}D^{1/2}D^{1/2}Xa.$
(51)

Assuming the data distribution is non-degenerate,
$a^{\prime}X^{\prime}D^{1/2}D^{1/2}Xa$ is positive and bounded away from zero
for all unit vectors $a\in\mathcal{R}^{K-1}$. Thus, as $n\to\infty$, the
maximum eigenvalue of $\hat{V}_{n}(\beta)$ is strictly less than unity. This
proves the result.

##### Condition (ii).

We want to show that $M(\cdot)$ is a contraction mapping with fixed-point
$\beta$. Following DS, a sufficient condition to satisfy (ii) for contraction
mappings exhibiting linear convergence is
$t(n)\geq-\frac{1}{2}\log(n)/\log(\kappa)$, where $\kappa\in[0,1)$ is the
modulus of the contraction $M(\cdot)$, which can be estimated as the mean or
median of
$\hat{\kappa}=|\hat{\beta}_{t+1}-\hat{\beta}_{t}|/|\hat{\beta}_{t}-\hat{\beta}_{t-1}|$
across several iterations. We must hence show that $M(\cdot)$ is a contraction
mapping converging linearly to $\beta$.

First, let us show that $\beta$ is a fixed-point of $M$. We have

$M(\beta)=E[X_{i}X_{i}^{\prime}]^{-1}E\left[X_{i}\tilde{Y}_{i}(\beta)\right],$
(52)

from which substituting $\tilde{Y}_{i}(\beta)$ yields

$M(\beta)=E[X_{i}X_{i}^{\prime}]^{-1}E\left[X_{i}(X_{i}^{\prime}\beta+\overline{\upsilon}_{i})\right].$
(53)

Rearranging and making use of Assumption 3 gives

$M(\beta)=\beta.$ (54)

Let us now show that $M$ is a contraction mapping exhibiting linear
convergence. Letting $\beta_{t}$ be the parameter after $t$ iterations, we
have

$\beta_{t+1}-\beta=M(\beta_{t})-\beta=M(\beta_{t})-M(\beta),$ (55)

because $M(\beta)=\beta$ by definition. By the mean value theorem, there is a
$b_{t}$ between $\beta_{t}$ and $\beta$ satisfying
$\beta_{t+1}-\beta=M(\beta_{t})-M(\beta)=(\beta_{t}-\beta)V(b_{t}),$ where
$V(\cdot)$ denotes the gradient. Let $e_{t}=||\beta_{t}-\beta||$, where
$||\cdot||$ denotes the sup norm, and thus $e_{t+1}=e_{t}||V(c_{t})||_{op}$,
with $||\cdot||_{op}$ denoting the operator version of the sup norm. A
standard algebra result and the symmetry of the matrix $V(\beta)$, composed of
$K\times L$ elements $E[\delta X_{ki}X_{li}/(\delta+U_{i})]$, imply that
$||V(\beta)||_{op}$ is bounded by the largest eigenvalue of $V(\beta)$. Using
similar derivations than for condition (i), or by applying the limit as
$n\to\infty$, we have that $||V(\beta)||_{op}<1$. Therefore, by the continuity
of $V(\cdot)$ there is a small neighborhood around $\beta$ for which

$||V(\beta)||_{op}<\frac{\kappa+1}{2}<1.$ (56)

If $\beta_{t}$ lies in this neighborhood, then so does $c_{t}$. Therefore, we
have $||e_{t+1}||\leq\frac{\kappa+1}{2}||e_{t}||$, and
$\lim_{t\to\infty}\frac{||e_{t+1}||}{||e_{t}||}=\lim_{t\to\infty}||V(c_{t})||_{op}=||V(\beta)||_{op}=\kappa<1,$
which provides the desired result.

##### Condition (iii).

We now want to show that $\hat{M}_{n}$ converges uniformly to $M$, i.e.
$n^{1/2}\sup_{\phi\in B_{0}}|\hat{M}_{n}(\phi)-M(\phi)|=O_{p}(1)$ as
$n\to\infty$.

For any $\phi\in B_{0}$, recall that
$\hat{M}_{n}(\phi)=X^{\prime}X^{-1}X^{\prime}\hat{\tilde{Y}}_{i}(\phi)$. Under
the iid assumption and assuming $E[X_{i}X_{i}^{\prime}]<\infty$, applying the
weak law of large numbers and Slutsky’s theorem yield
$n^{-1}X^{\prime}X^{-1}\overset{p}{\to}E[X_{i}X_{i}^{\prime}]^{-1}$ and
$\hat{c}(\phi)\overset{p}{\to}c(\phi)$ as $n\to\infty$, and thus
$n^{-1}X^{\prime}\hat{\tilde{Y}}_{i}(\phi)\overset{p}{\to}E\left[X_{i}\tilde{Y}_{i}(\phi)\right]$
as $n\to\infty$. Therefore, $\hat{M}_{n}(\phi)\overset{p}{\to}M(\phi)$ as
$n\to\infty$ and the Lindeberg-Levy’s central limit theorem gives
$|\hat{M}_{n}(\phi)-M(\phi)|=O_{p}(n^{-1/2})$ for any $\phi\in B_{0}$. To show
uniform convergence, let us recall that $B_{0}$ is a closed ball in a
Euclidean space and so is compact. We obtain the following inequality

$\displaystyle|\hat{M}_{n}(\phi)-M(\phi)|\leq$
$\displaystyle|n^{-1}\sum_{i=1}^{n}X_{i}\log(Y_{i}+\exp(X_{i}^{\prime}\phi))-E[X_{i}\log(Y_{i}+\exp(X_{i}^{\prime}\phi))]|$
(57) $\displaystyle+|n^{-1}\sum_{i=1}^{n}X_{i}\hat{c}(\phi)-E[X_{i}]c(\phi)|$
$\displaystyle\leq$
$\displaystyle|n^{-1}\sum_{i=1}^{n}X_{i}\log(Y_{i}+\exp(X_{i}^{\prime}\phi_{l}))-E[X_{i}\log(Y_{i}+\exp(X_{i}^{\prime}\phi_{l}))]|$
$\displaystyle+|n^{-1}\sum_{i=1}^{n}X_{i}-E[X_{i}]||\hat{c}(\phi_{u})|+|E[X_{i}]||\hat{c}(\phi_{u})-c(\phi_{u})|$

where the first inequality follows from prewhitening and the triangular
inequality; the second inequality follows from the compactness of $B_{0}$, by
which there exist $\phi_{u}$ and $\phi_{l}$ in the parameter set such that the
inequality holds, and by the triangular inequality. All three terms on the
right-hand-side (RHS) are finite, and consist in averages of zero-mean iid
random variables with finite first and second moments by assumption, and thus
have order $O_{p}(n^{-1/2})$. We deduce the uniform convergence result from
the continuity of $\hat{M}_{n}(\phi)$ and $M(\phi)$ in $\phi$ by applying
Lemma 2.4 in Whitney and McFadden (1994).

##### Condition (iv).

Let us use the derivations obtained earlier and similar arguments than for
condition (iii). We have that
$\nabla_{\phi}\hat{c}(\phi)\overset{p}{\to}\nabla_{\phi}c(\phi)$ and thus
$\hat{V}_{n}(\phi)\overset{p}{\to}V(\phi)$ as $n\to\infty$. Therefore, the
condition $||\hat{V}_{n}(\phi)-V(\phi)||=o_{p}(1)$ holds. Uniform convergence
follows from similar derivations to obtain

$\displaystyle||\hat{V}_{n}(\phi)-V(\phi)||\leq$
$\displaystyle||n^{-1}\sum_{i=1}^{n}X_{i}\frac{\exp(X_{i}^{\prime}\phi_{l})}{Y_{i}+\exp(X_{i}^{\prime}\phi_{l})}X_{i}^{\prime}-E[X_{i}\frac{\exp(X_{i}^{\prime}\phi_{l})}{Y_{i}+\exp(X_{i}^{\prime}\phi_{l})}X_{i}^{\prime}]||$
(58)
$\displaystyle+|n^{-1}\sum_{i=1}^{n}X_{i}-E[X_{i}]|\cdot||\nabla_{\phi}\hat{c}(\phi_{u})||+|E[X_{i}]|\cdot||\nabla_{\phi}\hat{c}(\phi_{u})-\nabla_{\phi}c(\phi_{u})||,$

where all three terms on the RHS are finite and have finite first and second
moments by assumption. All conditions being satisfied, we apply Theorem 1 in
DS to obtain the convergence result
$n^{1/2}|\hat{\beta}_{t(n)}-\beta|=O_{p}(1)$ as $n\to\infty$.

######  Proof 2 (Proof of Theorem 1: Normality)

We now make use of Theorem 4 in DS to derive the asymptotic distribution of
iOLS. All conditions have been verified in the previous results except that
$\sqrt{n}(\hat{M}_{n}(\beta)-\beta)\overset{d}{\to}Z$ as $n\to\infty$, where
$Z$ is a limit distribution. Note that we have

$\hat{c}(\beta)=n^{-1}\sum_{i=1}^{n}\log(n^{-1}\sum_{j=1}^{n}U_{j}+U_{i})-\log(n^{-1}\sum_{j=1}^{n}U_{j})\overset{p}{\to}E[\log(1+U_{i})]=c,$
(59)

as $n\to\infty$, and
$\hat{\tilde{Y}}_{i}(\beta)=\log(1+U_{i})+X_{i}^{\prime}\beta-\hat{c}(\beta),$
so that

$\sqrt{n}[X^{\prime}X]^{-1}X^{\prime}\hat{\tilde{Y}}_{i}(\beta)=\sqrt{n}\left(\beta+[X^{\prime}X]^{-1}X^{\prime}(\log(1+U)-\hat{c}(\beta))\right).$
(60)

Under the iid assumption and the exogeneity condition
$E[X_{i}\log(1+U_{i})]=c$, the Lindeberg-Levy’s central limit theorem yields

$\sqrt{n}\left(\hat{M}_{n}(\beta)-\beta\right)\overset{d}{\to}\mathcal{N}(0,\Sigma),$
(61)

as $n\to\infty$, where $\Sigma$ is the asymptotic covariance matrix. Remark
that it is the asymptotic covariance of the OLS estimator of the regression of
$\hat{\tilde{Y}}(\beta)$ onto $X$. Heteroskedasticity-robust estimators and
alike apply exactly as in the standard OLS setting. However, the iOLS
estimator has a slightly different asymptotic distribution. Theorem 4 of DS
2005 gives
$\sqrt{n}\left(\hat{\beta}_{i(n)}-\beta\right)\overset{d}{\to}\mathcal{N}(0,\Omega^{-1}),$
as $n\to\infty$, where $\Omega=(I_{k}-V(\beta))^{-1}\Sigma(I_{K}-V(\beta))$
and the gradient $\nabla_{\phi}M(\beta)=V(\beta)$ is defined as

$V(\beta)=E[X_{i}X_{i}^{\prime}]^{-1}E[\frac{X_{i}X_{i}^{\prime}}{1+U_{i}}],$
(62)

of which each element is strictly below $1$. Therefore sandwich-type
covariance estimators are changed from the classical expression
$\hat{\Sigma}=(\frac{1}{n}X^{\prime}X)^{-1}\hat{\Sigma}_{0}(\frac{1}{n}X^{\prime}X)^{-1}$
to

$\tilde{\Sigma}=(\frac{1}{n}X^{\prime}(I-W)X)^{-1}\hat{\Sigma}_{0}(\frac{1}{n}X^{\prime}(I-W)X)^{-1},$
(63)

where $W$ is a diagonal weighting matrix with diagonal element
$\frac{1}{1+U_{i}}$, and $\hat{\Sigma}_{0}$ is an estimator of the covariance
of $X_{i}^{\prime}(\log(1+U_{i})-c)$ across observations. For another
$\delta\neq 1$, we would have the weights
$\frac{\delta}{\delta+U_{i}}\in[0,1)$. In layman’s terms, the “meat” of HAC-
robust estimators is unchanged but the “bread” is modified. As before, the
weights become $\frac{\delta}{\delta+U_{i}}$ when $\delta\neq 1$.

######  Proof 3 (Proof of Theorem 2: Specification tests)

Let us start by assuming some regularity conditions to guarantee existence:
$E({W}_{i}^{2}|U_{i}>0)<\infty$, $E(W_{i}^{2}\nu_{i}^{2}|U_{i}>0)<\infty$,
$E(W_{i}^{2}U_{i}^{2}|U_{i}>0)<\infty$,
$E(\varepsilon_{i}^{WU^{2}}|U_{i}>0)<\infty$. We decompose $\hat{\lambda}-1$
into

$\displaystyle\hat{\lambda}-1=$
$\displaystyle\tilde{\lambda}-1+\left(n_{1}^{-1}\sum_{i=1}^{n_{1}}W_{i}^{2}\right)^{-1}\left(\sum_{j=1}^{n_{1}}\hat{W}_{j}\hat{U}_{j}-\sum_{j=1}^{n_{1}}W_{j}U_{j}\right)$
(64)
$\displaystyle+\left[\left(n_{1}^{-1}\sum_{i=1}^{n_{1}}\hat{W}_{i}^{2}\right)^{-1}-\left(n_{1}^{-1}\sum_{i=1}^{n_{1}}W_{i}^{2}\right)^{-1}\right]\left(\sum_{j=1}^{n_{1}}\hat{W}_{j}\hat{U}_{j}\right),$

where
$\tilde{\lambda}=\left(n_{1}^{-1}\sum_{i=1}^{n_{1}}W_{i}^{2}\right)^{-1}\left(\sum_{j=1}^{n_{1}}W_{j}U_{j}\right)$
is the OLS estimator in absence of first-step errors.

In what follows, we omit the conditioning of expectations are conditioned on
$U_{i}>0$. If first-step estimators are consistent, and by
$E({W}_{i}^{2})<\infty$, which implies $E(W_{i}U_{i}|U_{i}>0)<\infty$ under
$H_{0}$, the weak law of large numbers for iid observations implies that
$n_{1}^{-1}\sum_{i=1}^{n_{1}}\hat{W}_{i}^{2}\overset{p}{\to}E({W}_{i}^{2})$,
$n_{1}^{-1}\sum_{i=1}^{n_{1}}\hat{W}_{i}\hat{U}_{i}\overset{p}{\to}E(W_{i}U_{i})$
as $n_{1}\to\infty$,
$n_{1}^{-1}\sum_{i=1}^{n_{1}}W_{i}^{2}\overset{p}{\to}E({W}_{i}^{2})$, and
$n_{1}^{-1}\sum_{i=1}^{n_{1}}W_{i}U_{i}\overset{p}{\to}E(W_{i}U_{i})$ as
$n_{1}\to\infty$. Moreover, $H_{0}$ implies $E[\nu_{i}|U_{i}>0,X_{i}]=0$ which
in turn yields $\tilde{\lambda}\overset{p}{\to}1$ as $n_{1}\to\infty$, leading
to the consistency result: $\hat{\lambda}\overset{p}{\to}1$.

Let us now characterize the asymptotic distribution of $\tilde{\lambda}$. We
decompose $\hat{W}_{i}\hat{U}_{i}=W_{i}U_{i}+\varepsilon_{i}^{WU}$, where the
estimation error $\varepsilon_{i}^{WU}$ is $o_{p}(1)$. Applying the Lindeberg-
Levy’s central limit theorem given the iid assumption and finite second-order
moments $E(W_{i}^{2}\nu_{i}^{2})<\infty$, $E(\varepsilon_{j}^{WU^{2}})<\infty$
imply that the first two terms in (64) are asymptotically normal as
$n_{1}\to\infty$. In particular, the firm term satisfies

$\sqrt{n}\left(\tilde{\lambda}-1\right)\overset{d}{\to}N(0,\sigma_{\lambda_{0}}^{2}),$
(65)

as $n_{1}\to\infty$ with
$\sigma_{\lambda_{0}}^{2}=E(W_{i}^{2})^{-2}E(W_{i}^{2}\nu_{i}^{2})$. The
second term in (64) satisfies

$\left(n_{1}^{-1}\sum_{i=1}^{n_{1}}W_{i}^{2}\right)^{-1}\sqrt{n}\sum_{j=1}^{n_{1}}\varepsilon_{j}^{WU}\overset{d}{\to}N(0,\sigma_{\lambda_{1}}^{2}),$
(66)

as $n_{1}\to\infty$ with
$\sigma_{\lambda_{1}}^{2}=E(W_{i}^{2})^{-2}E(\varepsilon_{j}^{WU^{2}})$.
Finally, the third term satisfies

$\left[\left(n_{1}^{-1}\sum_{i=1}^{n_{1}}\hat{W}_{i}^{2}\right)^{-1}-\left(n_{1}^{-1}\sum_{i=1}^{n_{1}}W_{i}^{2}\right)^{-1}\right]\sqrt{n}\sum_{j=1}^{n_{1}}\hat{W}_{j}\hat{U}_{j}\overset{p}{\to}0,$
(67)

as $n_{1}\to\infty$, because
$\sqrt{n}\sum_{j=1}^{n_{1}}\hat{W}_{j}\hat{U}_{j}$ is stochastically bounded
under $E(W_{i}^{2}U_{i}^{2}|U_{i}>0)<\infty$ and the term in square bracket
converges to zero in probability with $n_{1}$. Therefore, assuming
$W_{i}\nu_{i}$ and $\varepsilon_{i}^{WU}$ to be independent random variables,
conditionally on $U>0$, we have

$\sqrt{n}\left(\hat{\lambda}-1\right)\overset{d}{\to}N(0,\sigma_{\lambda}^{2}),$
(68)

where
$\sigma_{\lambda}^{2}=E(W_{i}^{2})^{-2}E(W_{i}^{2}\nu_{i}^{2})+E(W_{i}^{2})^{-2}E(\varepsilon_{j}^{WU^{2}})$,
and the test statistic (35) is asymptotically standard normal provided that
$\hat{\sigma}_{\lambda}$ is a consistent estimator of $\sigma_{\lambda}$.

Online Appendix

##  Appendix A Mathematical Appendix

###  A.1 Theorem 3: Poisson models estimated with iOLS

######  Proof 4 (Proof of Theorem 3: iOLSMP)

This proof is similar to that of Theorem 1, with small modifications.
Regularity conditions must be modified accordingly.

Let us now consider

$\displaystyle\hat{V}_{n}(\phi)$
$\displaystyle=[X^{\prime}X]^{-1}X^{\prime}\nabla_{\phi}\hat{\tilde{Y}}(\phi)$
(69) $\displaystyle=n^{-1}X^{\prime}\nabla_{\phi}\hat{\tilde{Y}}(\phi),$

where $\nabla_{\phi}\hat{\tilde{Y}}_{i}(\phi)$ has element $(i,k)$ defined as

$\displaystyle\left[\nabla_{\phi}\hat{\tilde{Y}}(\phi)\right]_{i,k}$
$\displaystyle=\frac{\delta\exp(X_{i}^{\prime}\phi)X_{ki}}{Y_{i}+\delta\exp(X_{i}^{\prime}\phi)}+\frac{\partial\hat{U}_{i}(\phi)}{\partial\phi_{k}}\left(\frac{1}{1+\delta}-\frac{1}{\hat{U}_{i}(\phi)+\delta}\right).$
(70)

This expression simplifies, when evaluated at $\phi=\beta$, to

$\displaystyle\left[\nabla_{\beta}\hat{\tilde{Y}}(\beta)\right]_{i,k}$
$\displaystyle=X_{ki}\left(1-\frac{U_{i}}{1+\delta}\right),$ (71)

which yields

$\displaystyle\left[\hat{V}_{n}(\beta)\right]_{k,l}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}X_{ki}X_{li}\left(1-\frac{U_{i}}{1+\delta}\right).$
(72)

Following the same reasoning as in the previous theorem, a sufficient
condition for convergence is that $\frac{U_{i}}{1+\delta}$ is between $0$ and
$1$ for all $i$. Therefore, the choice of $\delta$ will affect both the speed
of convergence and whether the estimator converges at all. An efficient
strategy for choosing $\delta$ is to start at a relatively small value and
increment it if convergence fails – which can be checked by estimating
$\kappa$ as explained above.

The proof of asymptotic normality is also unchanged, except that now the
diagonal weighting matrix $W$ in

$\tilde{\Sigma}=(\frac{1}{n}X^{\prime}(I-W)X)^{-1}\hat{\Sigma}_{0}(\frac{1}{n}X^{\prime}(I-W)X)^{-1},$
(73)

has element $1-\frac{U_{i}}{1+\delta}$, and $\hat{\Sigma}_{0}$ is an estimator
of the covariance of $X_{i}^{\prime}U_{i}$ across observations.

######  Proof 5 (Proof of Theorem 3: iOLSAP))

This proof follows the same lines, with small modifications to the previous
one. The gradient $\nabla_{\phi}\hat{\tilde{Y}}_{i}(\phi)$ has now element
$(i,k)$ defined as

$\displaystyle\left[\nabla_{\phi}\hat{\tilde{Y}}(\phi)\right]_{i,k}$
$\displaystyle=\frac{\delta\exp(X_{i}^{\prime}\phi)X_{ki}}{Y_{i}+\delta\exp(X_{i}^{\prime}\phi)}-\frac{1}{\hat{U}_{i}(\phi)+\delta}\frac{\partial\hat{U}_{i}(\phi)}{\partial\phi_{k}}+\frac{1}{1+\delta}\frac{\partial(Y_{i}-\exp(X_{i}^{\prime}\phi))}{\partial\phi_{k}}.$
(74)

This expression simplifies, when evaluated at $\phi=\beta$, to

$\displaystyle\left[\nabla_{\beta}\hat{\tilde{Y}}(\beta)\right]_{i,k}$
$\displaystyle=X_{ki}\left(1-\frac{\exp(X_{i}^{\prime}\beta)}{1+\delta}\right),$
(75)

which yields

$\displaystyle\left[\hat{V}_{n}(\beta)\right]_{k,l}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}X_{ki}X_{li}\left(1-\frac{\exp(X_{i}^{\prime}\beta)}{1+\delta}\right).$
(76)

Following the same reasoning as in the previous theorem, a sufficient
condition for convergence is that $\frac{\exp(X_{i}^{\prime}\beta)}{1+\delta}$
is between $0$ and $1$ for all $i$. We suggest using the same trial and error
approach based on estimating $\kappa$.

The proof of asymptotic normality is also unchanged, except that now the
diagonal weighting matrix $W$ in

$\tilde{\Sigma}=(\frac{1}{n}X^{\prime}(I-W)X)^{-1}\hat{\Sigma}_{0}(\frac{1}{n}X^{\prime}(I-W)X)^{-1},$
(77)

has element $1-\frac{\exp(X_{i}^{\prime}\beta)}{1+\delta}$, and
$\hat{\Sigma}_{0}$ is an estimator of the covariance of
$X_{i}^{\prime}\epsilon_{i}$ across observations.

###  A.2 Theorem 4: i2SLS

######  Proof 6 (Proof of Theorem 4 : Instrumental Variables Consistency)

Recall that the parameter $\beta\in\mathbb{R}^{K}$ is characterized by the
fixed-point equation

$\beta^{IV}=E[\breve{X_{i}}\breve{X_{i}}^{\prime}]^{-1}E\left[\breve{X_{i}}\tilde{Y}_{i}(\beta)\right],$
(78)

where $\breve{X}=P^{Z}X$, $P^{Z}=Z(Z^{\prime}Z)^{-1}Z^{{}^{\prime}}$,
$Z\in\mathbb{R}^{M}$ with $M\geq K$, $E(Z_{i}^{\prime}X_{i})$ has rank $K$,
and $\tilde{Y}_{i}(\beta)=\log(Y_{i}+\exp(X_{i}^{\prime}\beta))-c(\beta)$ is
the transformed dependent variable. The mapping from $\mathbb{R}^{K}$ to
$\mathbb{R}^{K}$ which characterizes the parameter is hence defined
$\forall\phi\in\mathbb{R}^{K}$ as

$M^{IV}(\phi)=E[\breve{X_{i}}\breve{X_{i}}^{\prime}]^{-1}E\left[\breve{X_{i}}\tilde{Y}_{i}(\phi)\right].$
(79)

The sample counterpart of this mapping is given by

$\hat{M}_{n}^{IV}(\phi)=[\breve{X_{i}}^{\prime}\breve{X_{i}}]^{-1}\breve{X_{i}}^{\prime}\hat{\tilde{Y}}_{i}(\phi),$
(80)

where $\hat{\tilde{Y}}_{i}(\phi)$ is defined as before.

Our proof is very similar to the one used to show Theorem 1. We do not state
the modified regularity conditions and only focus on showing condition (i)
because the others consist in simple extensions. For condition (i), the first
step is to consider that $Z$ is standardized so that $\breve{X}$ is
prewhitened: $\breve{X}^{\prime}\breve{X}=nI_{k}$. As before, showing
condition (i) boils down to showing that the largest eigenvalue of
$\nabla_{\phi}\hat{M}^{IV}_{n}(\beta)=\hat{V}^{IV}_{n}(\beta)$ is strictly
less that unity as $n\to\infty$. Note that we have

$\displaystyle\hat{V}^{IV}_{n}(\phi)$
$\displaystyle=[\breve{X}^{\prime}\breve{X}]^{-1}\breve{X}^{\prime}\nabla_{\phi}\hat{\tilde{Y}}(\phi)$
(81)
$\displaystyle=n^{-1}\breve{X}^{\prime}\nabla_{\phi}\hat{\tilde{Y}}(\phi),$

where the second equality uses prewhitening on $\breve{X}$. Moreover,
$\nabla_{\phi}\hat{\tilde{Y}}_{i}(\phi)$ has element $(i,k)$ defined as

$\displaystyle\left[\nabla_{\phi}\hat{\tilde{Y}}(\phi)\right]_{i,k}$
$\displaystyle=\frac{\exp(X_{i}^{\prime}\phi)X_{ki}}{Y_{i}+\exp(X_{i}^{\prime}\phi)}-\frac{\partial\hat{c}(\phi)}{\partial\phi_{k}}.$
(82)

Let us denote $X_{1i}=1$ and $Z_{1i}=1$, for all $i$ as the constant. By
prewhitening $\breve{X}$, we have $\sum_{j=1}^{n}\breve{X}_{1j}=n$ and
$\sum_{j=1}^{n}\breve{X}_{kj}=0$ for $k>1$. The derivative of the nuisance
parameter estimate writes

$\displaystyle\frac{\partial\hat{c}(\phi)}{\partial\phi_{k}}=n^{-1}\sum_{i=1}^{n}\frac{\exp(X_{i}^{r^{\prime}}\phi^{r}+\hat{\phi}^{1})(\frac{\partial\hat{\phi}^{1}}{\partial\phi_{k}}+X_{ki})}{Y_{i}+\exp(X_{i}^{r^{\prime}}\phi^{r}+\hat{\phi}^{1})}-n^{-1}\sum_{i=1}^{n}(\frac{\partial\hat{\phi}^{1}}{\partial\phi_{k}}+X_{ki}),$
(83)

for $k>1$ and $\frac{\partial\hat{c}(\phi)}{\partial\phi_{1}}=0$. As before,
this expression simplifies when evaluated at $\phi=\beta$, as shown by

$\displaystyle\frac{\partial\hat{c}(\beta)}{\partial\phi_{k}}=$ $\displaystyle
n^{-1}\sum_{i=1}^{n}\frac{X_{ki}}{1+U_{i}}-n^{-1}\sum_{i=1}^{n}X_{ki}+O_{p}(1)$
(84) $\displaystyle=$ $\displaystyle
n^{-1}\sum_{i=1}^{n}\frac{X_{ki}U_{i}}{1+U_{i}}+O_{p}(1),$

for $k>1$ because
$\hat{\phi}^{1}(\beta)=\log(n^{-1}\sum_{i=1}Y_{i}\exp(-X_{i}^{r^{\prime}}\beta^{r}))=\beta_{1}+\log(n^{-1}\sum_{i=1}U_{i})$,
where $\log(n^{-1}\sum_{i=1}U_{i})=O_{p}(1)$ by iid assumption and
$E[U_{i}]=1$.

Therefore, each element $(k,l)$ of $\hat{V}^{IV}_{n}(\beta)$ writes

$\displaystyle\left[\hat{V}^{IV}_{n}(\beta)\right]_{k,l}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\frac{\breve{X}_{ki}X_{li}}{1+U_{i}}-(n^{-1}\sum_{i=1}^{n}\breve{X}_{ki})(n^{-1}\sum_{j=1}^{n}\frac{X_{lj}U_{j}}{1+U_{j}}),$
(85)

for $l>1$ and

$\displaystyle\left[\hat{V}^{IV}_{n}(\beta)\right]_{k,l}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\frac{\breve{X}_{ki}}{1+U_{i}},$ (86)

for $l=1$. Remark that for $k=1,\forall l>1$ we have
$\left[V^{IV}_{n}(\beta)\right]_{1,l}=n^{-1}\sum_{i=1}^{n}\frac{X_{li}}{1+U_{i}}$,
and for $k=1,l=1$, we have
$\left[\hat{V}^{IV}_{n}(\beta)\right]_{1,1}=n^{-1}\sum_{i=1}^{n}\frac{1}{1+U_{i}}<1$.
Therefore, all elements $(k,l)$ for $k,l\geq 1$ of this matrix writes

$\displaystyle\left[\hat{V}^{IV}_{n}(\beta)\right]_{k,l}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\frac{\breve{X}_{ki}X_{li}}{1+U_{i}}.$ (87)

because of prewhitening. We can write this in matrix form as

$\displaystyle\left[\hat{V}^{IV}_{n}(\beta^{IV})\right]$
$\displaystyle=n^{-1}X^{\prime}P_{z}WX,$ (88)

where $W$ is a diagonal matrix with elements $(i,i)$ acting as weights given
by $\frac{1}{1+U_{i}}\in(0,1]$. The projection matrix $P_{z}$ being symmetric
and idempotent, its eigenvalues are equal to either $0$ or $1$. $P_{z}$ is
hence a positive semi-definite matrix. The product $P_{z}W$ is thus a positive
semi-definite matrix because it is the product of two symmetric positive semi-
definite matrices.

Nevertheless $P_{z}W$ is not necessarily symmetric. For any vector
$a\in\mathbb{R}^{K}$, $a^{\prime}X^{\prime}P_{z}WXa$ and
$a^{\prime}X^{\prime}\frac{1}{2}(P_{z}W+W^{\prime}Pz)Xa$ are the same
quadratic forms. We have that $X^{\prime}\frac{1}{2}(P_{z}W+W^{\prime}Pz)X$ is
positive semi-definite matrix and all its eigenvalues are nonnegative and
corresponds to those of $X^{\prime}P_{z}WX$.

We can alternatively write the weight matrix $W=I_{n}-D$, where $D$ is also a
diagonal matrix with elements $\frac{U_{i}}{1+U_{i}}\in[0,1)$. Therefore, we
have the alternative expression

$\displaystyle\left[\hat{V}^{IV}_{n}(\beta)\right]$
$\displaystyle=n^{-1}X^{\prime}P_{z}(I_{n}-D)X$ (89)
$\displaystyle=X^{\prime}P_{z}X-n^{-1}X^{\prime}P_{z}DX$
$\displaystyle=I_{K}-n^{-1}X^{\prime}P_{z}DX,$

where the second equality comes from $P_{z}$ being idempotent, and
prewhitening. It follows that as $n\to\infty$, the maximum eigenvalue is equal
to

$\displaystyle\max_{|a|=1}a^{\prime}\left[\hat{V}^{IV}_{n}(\beta)\right]a=\max_{|a|=1}1-a^{\prime}X^{\prime}\frac{1}{2}(P_{z}D+D^{\prime}P_{z})Xa.$
(90)

Assuming the data distribution is non-degenerate,
$a^{\prime}X^{\prime}\frac{1}{2}(P_{z}D+D^{\prime}P_{z})Xa$ is positive and
bounded away from zero for all unit vectors $a\in\mathcal{R}^{K}$. Thus, as
$n\to\infty$, the maximum eigenvalue of $\hat{V}^{IV}_{n}(\beta)$ is strictly
less than unity. This proves the result. The other conditions follow similar
derivations as for Theorem 1 which complete the proof.

######  Proof 7 (Proof of Theorem 4: Instrumental Variables Normality)

We now derive the asymptotic distribution of i2SLS. We must show that
$\sqrt{n}(\hat{M}^{IV}_{n}(\beta)-\beta)\overset{d}{\to}Z$ as $n\to\infty$,
where $Z$ is a limit distribution. As before, we have

$\hat{c}(\beta)\overset{p}{\to}E[\log(1+U_{i})]=c,$ (91)

as $n\to\infty$, and

$\hat{\tilde{Y}}_{i}(\beta)=\log(1+U_{i})+X_{i}^{\prime}\beta-\hat{c}(\beta),$
(92)

so that

$\sqrt{n}[\breve{X}^{\prime}\breve{X}]^{-1}\breve{X}^{\prime}\hat{\tilde{Y}}_{i}(\beta)=\sqrt{n}\left(\beta+[\breve{X}^{\prime}\breve{X}]^{-1}\breve{X}^{\prime}(\log(1+U)-\hat{c}(\beta))\right).$
(93)

Under the iid assumption and the exogeneity condition
$E[\breve{X}_{i}(\log(1+U_{i})-c)]=0$, the Lindeberg-Levy’s central limit
theorem yields

$\sqrt{n}\left(\hat{M}^{IV}_{n}(\beta)-\beta\right)\overset{d}{\to}\mathcal{N}(0,\Sigma),$
(94)

as $n\to\infty$, where $\Sigma$ is the asymptotic covariance matrix. Remark
that it is the asymptotic covariance of the 2SLS estimator of the regression
of $\hat{\tilde{Y}}(\beta)$ onto $X$ using $Z$ as IV. Heteroskedasticity-
robust estimators apply as in the standard setting. However, the i2SLS
estimator has a slightly different asymptotic distribution, because the true
$\beta$ is unknown. Using the same reasoning as for iOLS, we obtain

$\sqrt{n}\left(\hat{\beta}^{IV}_{i(n)}-\beta^{IV}\right)\overset{d}{\to}\mathcal{N}(0,[\Omega^{IV}]^{-1}),$
(95)

as $n\to\infty$, where
$\Omega^{IV}=(I_{k}-V^{IV}(\beta))^{-1}\Sigma(I_{K}-V^{IV}(\beta))^{-1}$ and
the gradient $\nabla_{\phi}M^{IV}(\beta)=V^{IV}(\beta)$ is defined as

$V(\beta)=E[\breve{X}_{i}\breve{X}_{i}^{\prime}]^{-1}E[\frac{\breve{X}_{i}X_{i}^{\prime}}{1+U_{i}}].$
(96)

Therefore sandwich-type covariance estimators are given by

$\tilde{\Sigma}=(\frac{1}{n}X^{\prime}\frac{1}{2}(P_{z}(I-W)+(I-W)P_{z})X)^{-1}\hat{\Sigma}_{0}(\frac{1}{n}X^{\prime}\frac{1}{2}(P_{z}(I-W)+(I-W)P_{z})X)^{-1},$
(97)

where $W$ is a diagonal weighting matrix with diagonal element
$\frac{1}{1+U_{i}}$, and $\hat{\Sigma}_{0}$ is an estimator of the covariance
of $P_{z}X^{\prime}(\log(1+U_{i})-c)$ across observations. Symmetrizing the
weight matrix, as explained in the proof of the preceding theorem, is required
to have a symmetric positive definite matrix, hence invertible.

##  Appendix B Model Extensions

###  B.1 Poisson models as iOLS

##### Multiplicative Poisson.

First, we consider the multiplicative version of the model. It relies on the
identifying assumption $E(U_{i}|X_{i})=1$, but only requires
$E((U_{i}-1)X_{i})=0$ for consistency. To enforce this condition, we can add
$\frac{1}{1+\delta}(U_{i}-1)$ on both sides of (10) and rearrange to obtain

$\log(Y_{i}+\delta\exp(X_{i}^{\prime}\beta))-\left(\log(\delta+U_{i})-\frac{1}{1+\delta}(U_{i}-1)\right)=X_{i}^{\prime}\beta+\frac{1}{1+\delta}(U_{i}-1).$
(98)

with $U_{i}=Y_{i}\exp(-X_{i}^{\prime}\beta)$, the second term on the left-
hand-side can be rewritten into

$c_{i}(\delta,\beta)=\log(\delta+Y_{i}\exp(-X_{i}^{\prime}\beta))-\frac{1}{1+\delta}(Y_{i}\exp(-X_{i}^{\prime}\beta)-1),$
(99)

to obtain a new transformed dependent variable

$\tilde{Y}_{i}(\beta)=\log(Y_{i}+\delta\exp(X_{i}^{\prime}\beta))-c_{i}(\delta,\beta).$
(100)

and associated model

$\tilde{Y}_{i}(\beta)=X_{i}^{\prime}\beta+\eta_{i},$ (101)

where $\eta_{i}=\frac{1}{1+\delta}(U_{i}-1)$ is a mean-zero error term, and is
exogenous to $X_{i}$. This estimator will be referred to as iOLSMP. The choice
of $\delta$ will be discussed shortly.

##### Additive Poisson (PPML).

Similarly, one can enforce the additive representation based on model (4),
which assumes $E[\epsilon_{i}|X_{i}]=0$ with
$\epsilon_{i}=Y_{i}-\exp(X_{i}^{\prime}\beta)$. This assumption is equivalent
to $E[U_{i}|X]=1$ but leads to a different least-squares objective function.
iOLS can be adapted to this setting by adding and substracting
$\frac{1}{1+\delta}(Y_{i}-\exp(X_{i}^{\prime}\beta))$ to (10) and defining

$c_{i}(\delta,\beta)=\log(\delta+Y_{i}\exp(-X_{i}^{\prime}\beta))-\frac{1}{1+\delta}(Y_{i}-\exp(X_{i}^{\prime}\beta)).$
(102)

This estimator, hereafter referred to as iOLSAP, is equivalent to PPML but can
yield numerically different point estimates because the numerical algorithm
differs.

We derive the asymptotic properties of both estimators in the following
theorem.

######  Theorem 3 (Consistency and Normality of iOLSMP and iOLSAP)

Under Assumptions 1, 2, 4, and suitable regularity conditions, the iOLS
estimators using $c(\delta,\beta)$ in (99) and (102) are consistent, achieve
the parametric rate of convergence $n^{-1/2}$, and correspond to the
multiplicative and additive Poisson regression estimates, respectively.
Formally, we have $n^{1/2}|\hat{\beta}_{t(n)}-\beta|=O_{p}(1)$ as $n\to\infty$
for any $t(n)\geq-\frac{1}{2}\log(n)/\log(\kappa)$, where $\kappa\in[0,1)$ is
the modulus of the associated contraction mapping from $\mathbb{R}^{K}$ to
$\mathbb{R}^{K}$. In addition, they are asymptotically normally distributed
such that
$\sqrt{n}\left(\hat{\beta}_{t(n)}-\beta\right)\overset{d}{\to}\mathcal{N}(0,\Omega)$,
as $n\to\infty$, where $\Omega$ is the same than in Theorem 1 except for the
diagonal weighting matrix.

Unlike in the previous case, the parameter $\delta$ does not modify the
relevant moment condition. Its only purpose is now to control the convergence
of the algorithm. The modulus $\kappa$ is a function of $\delta$ with two
important features. First, the algorithm will diverge for too small values of
$\delta$, which ultimately depends on the underlying DGP, because it implies
$\kappa$ above 1. Second, a too large $\delta$ implies $\kappa$ very close to
1, hence a slow convergence as discussed earlier. Therefore, the optimal
$\delta$ is large enough to guarantee convergence but small enough so that
convergence is fast. We address these issues by defining a grid of values for
$\delta$ based on a sufficient condition guaranteeing convergence. We start
with a relatively small value and sequentially increment it if the algorithm
detects divergence.

###  B.2 Instrumental variables

The estimation of causal relationships is central to social sciences. Yet,
doing so is fraught with difficulties. Simultaneity, an omitted variable, or
the presence of measurement errors could result in biased estimates. For
example, if a researcher is interested in estimating the causal effect of the
number of police officers on crime, one may observe that the police is more
often deployed in areas where crime is high and conclude that police causes
more crime.

A popular solution consists on finding an instrumental variable which affects
the outcome only through the endogenous variable. Using variation in the
instrument, one can recover the impact of the main variable of interest on the
outcome through an estimation procedure known as Two Stage Least Squares
(2SLS).

Our iterated solution extends directly to this situation and consists, in
turn, in running 2SLS iteratively. Let us define $Z$ as a $n\times L$ matrix
with $L\geq K$ instrumental variables so that $E[X^{\prime}Z]\neq 0$. We
assume $E(Z^{\prime}Z)<\infty$ and denote $P_{z}$ as the projection matrix
$Z(Z^{\prime}Z)^{-1}Z^{\prime}$. The following algorithm characterizes the
i2SLS estimators.

######  Algorithm 2 (i2SLS estimator)

Let $\hat{\beta}_{0}$ be an initial estimate, for instance the 2SLS “popular
fix” estimator
$\hat{\beta}^{2SPF}=[X^{\prime}P_{z}X]^{-1}X^{\prime}P_{z}\log(Y+\Delta)\in\mathbb{R}^{K}$,
for some $\Delta>0$. the i2SLS estimator is obtained as follows.

  1. 1.

Initialize $t$ at $0$;

  2. 2.

Transform the dependent variable into $\tilde{Y}(\hat{\beta}_{t})$;

  3. 3.

Compute the 2SLS estimate
$\hat{\beta}^{2SLS}_{t+1}=(X^{\prime}P_{z}X)^{-1}(X^{\prime}P_{z}\hat{\tilde{Y}}(\hat{\beta}_{t}))$,
and update $t$ to $t+1$;

  4. 4.

Iterate steps 2 and 3 until $\hat{\beta}^{2SLS}_{t}$ converges.

This iterative estimator converges under some conditions on
$\tilde{Y}(\cdot)$. The same transformations studied earlier apply without
further modifications. We prove the consistency of this estimator using the
following assumptions in the next theorem.

######  Assumption 7 (Covariates)

$X$ and $Z$ have full column rank and $E(X_{i}X_{i}^{\prime})<\infty$ and
$E(Z_{i}Z_{i}^{\prime})<\infty$.

###### Assumption 8

The error term $\upsilon_{i}$ satisfies the weak exogenous restriction
$E[Z_{i}^{\prime}(\upsilon_{i}-c(\delta,\beta))]=0$ where the value
$c(\delta,\beta)$ is unknown. In addition, let $E[U_{i}]=1$.

######  Theorem 4 (Consistency and Asymptotic Normality)

Under Assumptions 2, 7, 8, and suitable regularity conditions, the i2SLS
estimator is consistent and achieves the parametric rate of convergence
$n^{-1/2}$. Formally, we have

$n^{1/2}|\hat{\beta}_{t(n)}^{IV}-\beta|=O_{p}(1)$ (103)

as $n\to\infty$ for any $t(n)\geq-\frac{1}{2}\log(n)/\log(\kappa)$, where
$\kappa\in[0,1)$ is the modulus of the associated contraction mapping from
$\mathbb{R}^{K}$ to $\mathbb{R}^{K}$. In addition, the i2SLS estimator is
asymptotically normally distributed such that

$\sqrt{n}\left(\hat{\beta}_{t(n)}^{IV}-\beta\right)\overset{d}{\to}\mathcal{N}(0,\Omega^{IV}),$
(104)

as $n\to\infty$, where $\Omega^{IV}$, as given in the proof, corresponds to
the asymptotic covariance of the 2SLS estimator in the last iteration up to
minor modifications.

This asymptotic result reveals several desirable properties of our procedure.
First, the i2SLS estimators can be obtained easily using available software.
Second, this iterative procedure makes non-linear instrumental variable
estimation computationally tractable even when many control variables are
included. This is particularly important because current count models are hard
to estimate when using instrumental variables.434343For example, to our
knowledge, there are no packages in Stata which allow one to estimate
instrumental variable count models, as in Mullahy (1997), with many
categorical control variables. Finally, researchers often rely on the control
function approach in non-linear models. This method requires the error in the
second stage to be an additively separable function of the first-stage error
and an independent error term. It also rules out settings where the endogenous
variable is not continuous (Wooldridge, 2015). In contrast, 2SLS (and thus
i2SLS) does not require such assumptions and can leverage the Frisch-Waugh-
Lovell theorem to greatly alleviate computations when many fixed-effects are
included.

Finally, the specification tests developed for iOLS are easily adapted for
situations with endogenous regressors. The main difference is that one must
estimate $Pr(Y>0|Z)$ instead of $Pr(Y>0|X)$. Further details are provided in
Appendix B.8.

###  B.3 Negative values

Our estimator extends to dependent variables taking negative values. However,
one needs to specify a slightly different model. We modify (1) into

$Y_{i}=\alpha+\exp(X_{i}^{\prime}\beta)U_{i},$ (105)

where $\exp(X_{i}^{\prime}\beta)U_{i}\geq 0$ is as before but $\alpha<0$
shifts leftwise the log function’s vertical asymptote at zero towards the
minimum value of $Y$. It is hence fundamentally different from the IHS, which
imposes a S-shape transformation around zero. The iOLS transformation becomes

$\log(Y_{i}-\alpha+\delta\exp(X_{i}^{\prime}\beta))=X_{i}^{\prime}\beta+\log(\delta+U_{i})$
(106)

where $Y_{i}-\alpha\geq 0$. Estimation requires an additional step to find
$\alpha$ before proceeding with the iOLS algorithm to address observations
bunched at the lower bound (instead of zero). In this model, $Y_{i}$ is
bounded below by $\alpha$. The bound can be estimated by taking the first-
order statistic $\hat{\alpha}=\min\limits_{i}Y_{i}$. Here, the convergence
rate of $\hat{\alpha}$ is crucial to determine that of the iOLS estimator. For
instance, if $Y_{i}$ is uniformly distributed, the first-order statistic will
converge at rate $n^{-1}$ to the true lower bound and the convergence result
of the iOLS estimator will remain unaffected. Reversely, slower convergence
rates will prevail if the first-order statistic converges at a rate slower
than $n^{-1/2}$.

###  B.4 Log-log specifications

In many econometric applications, the main parameter of interest is an
elasticity of $Y_{i}$ with respect to some variable $X_{i}$. Elasticities are
often estimated using a log-log specification. However, it is common to have
both dependent and independent variables that are equal to zero for some
observations. Taking the log-transform of either of these variables is
impossible. We propose to address this issue as follows.

Let us consider the following data generating process

$Y_{i}=X_{i}^{\beta}U_{i},$ (107)

with $X_{i}>0$ and $U_{i}\geq 0$. The iOLSδ estimator directly applies using
the transformation

$\log(Y_{i}+\delta X_{i}^{\beta})=\beta\log(X_{i})+\eta_{i},$ (108)

under the exogeneity restriction $E[\log(X_{i})\eta_{i}]=0$, where
$\eta_{i}=\log(\delta+U_{i})-c$ is the mean-zero error term of the linearized
model. The only difference with the log-linear setting is that the regressors
are also in log form.

A potential issue arises when $X_{i}$ can take zero values with positive
probability. For any independent variable, let us rewrite the above
restriction as

$E[\log(X_{i})\eta_{i}|X_{i}>0]Pr(X_{i}>0)+\lim\limits_{\epsilon\to
0}E[\log(\epsilon)\eta_{i}|X_{i}=0]Pr(X_{i}=0)=0,$ (109)

which can be rewritten into

$E[\log(X_{i})\eta_{i}|X_{i}>0]Pr(X_{i}>0)+\lim\limits_{\epsilon\to
0}\log(\epsilon)E[\mathbbm{1}_{(X_{i}=0)}\eta_{i}]Pr(X_{i}=0)=0.$ (110)

A sufficient condition for this equality to hold is to have both
$E[\log(X_{i})\eta_{i}|X_{i}>0]=0$ and $E[\mathbbm{1}_{(X_{i}=0)}\eta_{i}]=0$.
The former is the standard exogeneity condition stated for non-negative values
of $X_{i}$, whereas the latter means that the occurrences of zeros in $X_{i}$
are exogenous to the errors. In the single covariate setting, one can simply
discard observations where $X_{i}=0$ and estimate the model based on the
condition $E[\log(X_{i})\eta_{i}|X_{i}>0]=0$. In the multivariate case, this
approach would lead to discard possibly many observations and greatly dampen
statistical power. Instead, one can make use of both restrictions and
introduce an extra binary variable in the modelas in

$\log(Y_{i}+X_{i}^{\beta})=\beta_{0}\mathbbm{1}_{(X_{i}=0)}+\beta\tilde{X}_{i}+\eta_{i},$
(111)

where $\tilde{X}_{i}=\log(X_{i})$ for $X_{i}>0$ and is equal to $0$
otherwise.For ease of exposition, we have supposed the existence of a single
explanatory variable but this strategy can be used along with an intercept and
other covariates.

###  B.5 High-dimensional fixed-effects

The inclusion of fixed-effects creates some computational issues in non-linear
panel data models. A modified version of the iOLS/i2SLS algorithm can be used
to accommodate many fixed effects by making use of the Frisch-Waugh-Lovell
theorem as follows. Let us decompose the set of regressors $X=[X_{0},X_{1}]$,
where $X_{0}$ are binary variables capturing all fixed-effects and $X_{1}$ the
remaining regressors (including the constant term). Define the projection
matrix $P_{0}=X_{0}(X_{0}^{\prime}X_{0})^{-1}X_{0}^{\prime}$ and denote the
aggregate fixed-effect term by $\Lambda=X_{0}^{\prime}\beta_{0}$.

######  Algorithm 3 (iOLS estimator with many fixed effects)

Let $\hat{\beta}_{0}$, and $\hat{\Lambda}_{0}$ be initial estimates. The iOLS
estimator is defined as the following iterative procedure:

  1. 1.

Initialize $t$ at $0$;

  2. 2.

Transform the dependent variable into
$\tilde{Y}_{iOLS}(\hat{\beta}_{t},\hat{\Lambda}_{t})$, where the term
$X^{\prime}\hat{\beta}_{t}$ is replaced by
$X_{1}^{\prime}\hat{\beta}_{t}+\hat{\Lambda}_{t}$;

  3. 3.

Partial out the transformed dependent variable