18  IV in Poisson with Fixed Effects

Poisson regression specifies the conditional mean as an exponential function. As a quasi-MLE it needs that mean to be correct, not the full count distribution. If a regressor is endogenous, the IV strategy has to respect the nonlinear mean — and this is where intuition carried over from 2SLS starts to mislead. This chapter works through what goes wrong, what to do instead, and how much the choice costs.

Code 
using Panelest
using Distributions
using DataFrames, DataFramesMeta
using StatsModels
import StatsAPI: coeftable
using Random, StatsBase, LinearAlgebra
using Printf

18.1 The problem, and three routes to it

Throughout, \(y_1\) is a count, \(y_2\) is the endogenous explanatory variable, and \(\mathbf{z} = (\mathbf{z}_1, \mathbf{z}_2)\) collects the exogenous variables — \(\mathbf{z}_1\) appearing in the structural equation and \(\mathbf{z}_2\) excluded from it. The structural conditional mean is exponential,

\[E(y_1 \mid \mathbf{z}_1, y_2, c_1) = \exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + c_1), \tag{18.1}\]

with \(c_1\) an unobservable correlated with \(y_2\). Everything below is about how to remove that correlation.

18.1.1 A note on the notation: the subscript indexes equations, not variables

This chapter follows Wooldridge (2010) (sec. 18.5), because it quotes his equations (18.37)–(18.45) directly and the symbols need to match the source.

The subscript is the source of most of the confusion, so it is worth saying what it means: it indexes the equation, not the variable. \(y_1\) is the variable explained by equation 1, the structural equation; \(y_2\) is the variable explained by equation 2, the reduced form. The same convention governs \(c_1\), \(v_2\), \(\boldsymbol{\delta}_1\) and \(\boldsymbol{\pi}_2\) — each subscript says which equation the object belongs to. This is the simultaneous-equations convention, and it scales: a second endogenous regressor is \(y_3\), with no renaming.

One deliberate deviation: Wooldridge writes the coefficient on \(y_2\) as \(\gamma_1\); this chapter calls it \(\alpha_1\) throughout. Quoted equations are left exactly as he writes them.

18.1.2 The three estimators

Three estimators recur, and it is worth naming them by what they do rather than by letters, since the names then carry the argument:

Name First stage What enters the second stage Family
Linear-residual CF OLS projection of \(y_2\) on \(\mathbf{z}\) the residual \(\hat v_2\), inside \(\exp(\cdot)\) control function
Generalized-residual CF probit or logit for \(y_2\) the generalized residual \(\hat r_2\), inside \(\exp(\cdot)\) control function
Multiplicative-moment GMM any model for \(\Pr(y_2=1 \mid \mathbf{z})\) the fitted value \(\hat p_2\), as an instrument GMM / IV

Earlier drafts of this chapter called these Approach A, B and C. The letters are gone. Linear-residual CF was A, generalized-residual CF was B, and multiplicative-moment GMM was C. The third is not a control function at all: it does not modify the exponential mean, it changes the moment conditions. Grouping all three as “approaches” obscured that, and the distinction matters for when each is valid.

The plan: the theory first, in one place; then the three compared side by side; then simulations, continuous case and binary case.

18.2 Theory

18.2.1 The control function idea

18.2.1.1 A linear second stage: the CF is 2SLS

Set the exponential mean aside and take the linear analogue of Equation 18.1,

\[y_1 = \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \varepsilon_1, \qquad E(\mathbf{z}'\varepsilon_1) = 0,\]

which is what defines \(\varepsilon_1\). Two-stage least squares fits

\[ \begin{aligned} \text{first stage:} &\quad y_2 = \mathbf{z}\boldsymbol{\pi}_2 + v_2, \\ \text{second stage:} &\quad y_1 = \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 \hat y_2 + (\varepsilon_1 + \alpha_1 \hat v_2), \end{aligned} \]

with \(\hat y_2 = \mathbf{z}\hat{\boldsymbol{\pi}}_2\) and \(\hat v_2 = y_2 - \hat y_2\). The second-stage error is not \(\varepsilon_1\): replacing \(y_2\) by \(\hat y_2\) moves \(\alpha_1\hat v_2\) into it. Add \(\hat v_2\) to an OLS second stage instead, and you recover the 2SLS coefficient exactly. The algebra:

  1. First stage. \(y_2 = \mathbf{z}\boldsymbol{\pi}_2 + v_2\), by OLS. The residual satisfies \(E(\mathbf{z}'\hat v_2) = 0\) by construction.
  2. Error decomposition. Project \(\varepsilon_1\) on \(v_2\): \(\varepsilon_1 = \rho v_2 + \eta_1\) with \(\rho = \mathrm{Cov}(\varepsilon_1, v_2)/\mathrm{Var}(v_2)\). Then \(E(v_2\eta_1) = 0\) by construction of the projection, and \(E(\mathbf{z}'\eta_1) = E(\mathbf{z}'\varepsilon_1) - \rho E(\mathbf{z}'v_2) = 0\).
  3. Augmented second stage. \(y_1 = \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho\hat v_2 + \eta_1\).

OLS on step 3 is the control function estimator, and Frisch–Waugh–Lovell makes it numerically identical to 2SLS. No assumption beyond \(E(\mathbf{z}'\varepsilon_1) = 0\) is required: step 2 is a projection, so both of its orthogonality conditions are consequences of that moment condition and the first stage, not additions to them.

18.2.1.2 A nonlinear second stage: the CF is a genuine assumption

With an exponential mean, that safety net disappears. A direct IV/GMM approach solves moment conditions in \(\mathbf{z}\) without touching the exponential mean, and the form of those moments matters:

  • Additive moments \(E[\mathbf{z}'(y_1 - \exp(\mathbf{x}\boldsymbol{\beta}))] = 0\) assume an additive error orthogonal to \(\mathbf{z}\).
  • Multiplicative moments \(E[\mathbf{z}'(y_1 e^{-\mathbf{x}\boldsymbol{\beta}} - 1)] = 0\) (Mullahy 1997) assume a multiplicative error with mean independent of \(\mathbf{z}\).

When the unobservable enters inside the exponential — \(y_1\) has mean \(\exp(\mathbf{x}\boldsymbol{\beta} + \rho v)\), as in both DGPs of this chapter — the error is multiplicative, the additive moments are misspecified, and additive-moment IV is inconsistent. The multiplicative form is the robust default, and both are coded below so the difference is visible.

The control function does something different. It imposes the stronger conditional mean assumption (Wooldridge 2010, sec. 18.5; 2014), Wooldridge’s eq. (18.40):

\[E[y_1 \mid \mathbf{z}_1,\, y_2,\, v_2] = \exp\!\bigl(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho v_2\bigr). \tag{18.2}\]

The unobservable enters the mean linearly inside \(\exp(\cdot)\). Under this, substituting \(\hat v_2\) for \(v_2\) and running Poisson QMLE is consistent for \(\alpha_1\). Without it, nothing rescues the procedure.

Property Linear (OLS/2SLS) Nonlinear (Poisson CF)
Extra assumption beyond IV exogeneity None — FWL applies CF mean assumption required
CF vs. 2SLS/GMM numerically Identical Different — FWL fails
SE from second stage alone Correct Understated — bootstrap needed
Test of endogeneity \(t\)-stat on \(\hat{v}_2\) (DWH) Same — \(t\)-stat on \(\hat{v}_2\)

Where the decomposition gets applied. The decomposition \(\varepsilon_1 = \rho v_2 + \eta_1\) is the same in both cases. What changes is where it enters: additively outside any transformation in the linear model, inside \(\exp(\cdot)\) in Poisson, so that the quasi-MLE score can absorb it. If \(\rho \neq 0\) and \(\hat v_2\) is omitted from the exponential mean, a multiplicative \(e^{\rho v_2}\) term correlated with \(y_2\) survives, biasing \(\hat\alpha_1\) in the direction of \(\rho\).

18.2.2 A continuous \(y_2\): the linear reduced form is enough

For continuous \(y_2\) the assumptions are the natural ones. Wooldridge’s eqs. (18.38)–(18.39) are

\[y_2 = \mathbf{z}\boldsymbol{\pi}_2 + v_2, \qquad c_1 = \rho_1 v_2 + e_1, \tag{18.3}\]

with \(v_2\) independent of \(\mathbf{z}\) and \(e_1\) independent of \((\mathbf{z}, v_2)\). Those two independence conditions deliver Equation 18.2 directly, and the estimator follows: regress, take residuals, add them to the index.

Note which condition is doing the work. It is independence, not uncorrelatedness — and Wooldridge flags this at the point he states the assumption, before discreteness enters the discussion at all (Wooldridge 2010, 743):

“(We could relax the independence assumptions to some degree, but we cannot just assume that \(v_2\) is uncorrelated with \(\mathbf{z}\) and that \(e_1\) is uncorrelated with \(v_2\).)”

The distinction is invisible for continuous \(y_2\) — a linear model with a normal additive error satisfies both — which is exactly why it is easy to carry the procedure across to a binary \(y_2\) without noticing that it has stopped holding.

18.2.2.1 Why independence, when the linear model needs only uncorrelatedness?

Because of the exponential. Substituting Equation 18.3 into Equation 18.1 gives

\[E(y_1 \mid \mathbf{z}, y_2, v_2) = \exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2)\cdot\exp(\rho_1 v_2)\cdot E[\exp(e_1) \mid \mathbf{z}, y_2, v_2], \tag{18.4}\]

and Equation 18.2 follows only if that final factor is a constant — which is why Wooldridge normalises \(E[\exp(e_1)] = 1\) rather than \(E(e_1) = 0\).

That is a condition on the moment generating function of \(e_1\), that is, on every moment at once, which is essentially full independence. Compare the linear model, where the requirement is \(E[\mathbf{z}'\varepsilon_1] = 0\) — a statement about one moment, delivered free by the OLS projection.

This is the same fact as “CF \(=\) 2SLS in a linear model, CF \(\neq\) 2SLS in Poisson,” seen from the assumption side rather than the algebra side. \(\exp(\cdot)\) is what turns a moment condition into a distributional one, and everything difficult in this chapter follows from it.

18.2.2.2 Absorb or average in the first stage: the answer is the same

In a panel the linear first stage can handle the unit effects two ways, and it makes no difference which. Regress \(y_{it2}\) on \(\mathbf{z}_{it}\) and unit dummies and you get the fixed-effects residual \(\hat{v}_{it2}^{FE}\). Regress it on \(\mathbf{z}_{it}\) and the unit means \(\bar{\mathbf{z}}_i\) — the Mundlak device — and you get \(\hat{v}_{it2}^{M}\). Put either into the second-stage FE-Poisson and \(\hat{\alpha}_1\) and \(\hat{\rho}_1\) come out numerically identical (Lin and Wooldridge 2019, app. A.1–A.2).

The reason is that Mundlak’s regression reproduces the within slope, so the two residuals differ only by a unit-specific constant, \(\hat{v}_{it2}^{FE} = \hat{v}_{it2}^{M} + g_i\). In the second-stage index that constant is absorbed:

\[\alpha_1 y_{it2} + \rho_1 \hat{v}_{it2}^{FE} + c_i = \alpha_1 y_{it2} + \rho_1 \hat{v}_{it2}^{M} + (c_i + \rho_1 g_i),\]

and \(c_i\) is a free parameter, so shifting it costs nothing. The model is reparameterised, not changed.

Two conditions on this. The averaging dimension in the first stage has to be the dimension the second stage absorbs. And it is a fact about linear first stages: once \(y_2\) is binary the first stage is nonlinear, the residual no longer shifts by a constant, and choosing between dummies and Mundlak means becomes a real decision — which is the subject of the next section.

18.2.3 A binary \(y_2\): why the linear reduced form fails

OLS always produces a residual uncorrelated with \(\mathbf{z}\). That is a property of the population linear projection and it holds whatever \(y_2\) looks like. But the control-function derivation needs \(v_2\) independent of \(\mathbf{z}\), which is a substantive assumption about the reduced form and not something OLS supplies.

For binary \(y_2\) it cannot hold. Conditional on \(\mathbf{z}\), the residual \(v_2 = y_2 - \mathbf{z}\boldsymbol{\pi}_2\) takes exactly two values,

\[v_2 \mid \mathbf{z} = \begin{cases} 1 - \mathbf{z}\boldsymbol{\pi}_2 & \text{with probability } p(\mathbf{z}) \\ -\mathbf{z}\boldsymbol{\pi}_2 & \text{with probability } 1 - p(\mathbf{z}),\end{cases} \tag{18.5}\]

so both the two support points and the probabilities of landing on them are functions of \(\mathbf{z}\). Its variance, \(\mathrm{Var}(v_2 \mid \mathbf{z}) = p(\mathbf{z})[1 - p(\mathbf{z})]\), is a function of \(\mathbf{z}\) too. Independence fails by construction, not by approximation.

Wooldridge states this outright, at (Wooldridge 2010, 746):

“The assumption that \(e_1\) in (18.39) is independent of \(v_2\) and \(\mathbf{z}\) rules out most cases where \(y_2\) has some discreteness, such as with binary, count, or corner responses. (In particular, the independence assumption likewise fails unless \(v_2\) is independent of \(\mathbf{z}\), and it would be very unusual for a discrete variable to be expressible as a linear equation with additive error independent of \(\mathbf{z}\).)”

That sentence is compressed, so it is worth unpacking. It names one assumption, but the failure runs through three links:

  1. What the correction needs. From Equation 18.4, \(e_1\) must be independent of \((\mathbf{z}, y_2, v_2)\), so that \(E[\exp(e_1)\mid\cdot]\) is constant.
  2. That silently requires \(v_2 \perp \mathbf{z}\). \(e_1\) is defined as what is left of \(c_1\) after projecting on \(v_2\), so if \(v_2\) depends on \(\mathbf{z}\), \(e_1\) inherits the dependence. This is what “the independence assumption likewise fails unless” refers to.
  3. And \(v_2 \perp \mathbf{z}\) cannot hold for discrete \(y_2\). This is the step he leaves as “very unusual”, and it is in fact airtight. Independence means the conditional distribution of \(v_2\) given \(\mathbf{z}\) is the same for every \(\mathbf{z}\). But by Equation 18.5 the support is \(\{-\mathbf{z}\boldsymbol{\pi}_2,\ 1-\mathbf{z}\boldsymbol{\pi}_2\}\), and both points slide as \(\mathbf{z}\) changes. A distribution whose support moves with \(\mathbf{z}\) cannot be independent of it — unless \(\boldsymbol{\pi}_2 = \mathbf{0}\), in which case there is no first stage at all. The same argument covers count and corner responses: the mass point sits at \(-\mathbf{z}\boldsymbol{\pi}_2\) and travels.

So the sentence reads: the correction needs \(e_1 \perp (v_2, \mathbf{z})\); that needs \(v_2 \perp \mathbf{z}\); and \(v_2 \perp \mathbf{z}\) cannot hold when \(y_2\) is discrete, because the residual’s support is pinned to \(\mathbf{z}\boldsymbol{\pi}_2\).

Note what does not break. \(E[\mathbf{z}'v_2] = 0\) still holds exactly — the projection is fine, and OLS has done nothing wrong. It is independence that fails, and only the exponential mean cares about the difference.

18.2.3.1 The subtlety: discreteness of the residual is not the problem

It is tempting to summarise this as “the residual is binary, so it cannot be a control function.” That argument proves too much, because the generalized residual of the next section is two-valued given \(\mathbf{z}\) as well.

The difference is not the discreteness of the derived residual. It is which object independence is asserted about. The linear-residual CF asserts it of the observed \(v_2\), where both the support points and their probabilities move with \(\mathbf{z}\). Terza’s setup asserts it of the continuous latent probit error, \(v_2 \mid \mathbf{z} \sim \mathrm{Normal}(0,1)\) — which is the assumption every probit already makes, and costs nothing extra.

So the fix is not to find a continuous control function. It is to relocate the independence assumption to a latent variable where it is plausible.

18.2.3.2 What survives: the test, but not the correction

The linear-residual construction is not worthless for binary \(y_2\). It remains fully valid as a Hausman-type exogeneity test, and Wooldridge is explicit about why (Wooldridge 2010, 744):

“For the purposes of getting a limiting chi-square distribution, it does not matter where the linear combination \(\hat v_2\) comes from. In other words, under the null hypothesis none of the assumptions we made about \((c_1, v_2)\) need to hold: \(v_2\) need not be independent of \(\mathbf{z}\), and \(e_1\) in equation (18.39) need not be independent of \(v_2\). Therefore, as a test, this procedure is very robust, and it can be applied when \(y_2\) contains binary, count, or other discrete variables. Unfortunately, if \(y_2\) is endogenous, the correction does not work without something like the assumptions made previously.”

That is the cleanest statement of the position: under \(H_0: \rho = 0\) there is no endogeneity left to approximate, so the test needs no reduced-form assumption. Once \(\rho \neq 0\), the correction needs one, and for binary \(y_2\) it is unavailable.

18.2.4 Terza’s exact correction

Terza (1998), exposited in Wooldridge (2010) (sec. 18.5) and Imbens and Wooldridge (2007) (sec. 3.3, eqs. 3.32–3.36), replaces the linear reduced form with a threshold model — Wooldridge’s eq. (18.43):

\[y_2 = \mathbb{1}[\mathbf{z}\boldsymbol{\pi}_2 + v_2 \geq 0], \qquad v_2 \mid \mathbf{z} \sim \mathrm{Normal}(0,1), \tag{18.6}\]

and keeps the independence assumption, now about the latent \(v_2\). The derivation then runs in two steps. With \((c_1, v_2)\) jointly normal, mean zero and independent of \(\mathbf{z}\), write \(c_1 = \rho_1 v_2 + e_1\) where \(e_1 \mid \mathbf{z}, v_2 \sim \mathrm{Normal}(0, \tau_1^2 - \rho_1^2)\). Conditioning on the latent error gives eq. (18.44),

\[E(y_1 \mid \mathbf{z}, v_2) = E[\exp(e_1)]\exp(\mathbf{x}_1\boldsymbol{\beta}_1 + \rho_1 v_2) = \exp\!\bigl((\tau_1^2 - \rho_1^2)/2 + \mathbf{x}_1\boldsymbol{\beta}_1\bigr)\exp(\rho_1 v_2), \tag{18.7}\]

using the lognormal mean. But \(v_2\) is not observed — only \(y_2\) is. Averaging over the truncated normal implied by \(y_2\) gives eq. (18.45), and evaluating that expectation yields the exact result:

\[E(y_1 \mid \mathbf{z}, y_2) = \exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2)\cdot h(y_2, \mathbf{z}\boldsymbol{\pi}_2, \rho) \tag{18.8}\] \[h(y_2, \mathbf{z}\boldsymbol{\pi}_2, \rho) = \exp(\rho^2/2)\left\{y_2\,\frac{\Phi(\rho + \mathbf{z}\boldsymbol{\pi}_2)}{\Phi(\mathbf{z}\boldsymbol{\pi}_2)} + (1-y_2)\,\frac{1-\Phi(\rho + \mathbf{z}\boldsymbol{\pi}_2)}{1-\Phi(\mathbf{z}\boldsymbol{\pi}_2)}\right\} \tag{18.9}\]

This is estimated by plugging in the first-stage probit’s \(\hat{\boldsymbol{\pi}}_2\) and jointly QMLE-fitting \((\boldsymbol\delta_1, \alpha_1, \rho)\) against this mean — a genuinely nonlinear second stage, not a linear regressor addition.

\(\rho\) here is a coefficient, not a correlation. It is the coefficient in the linear projection of \(c_1\) on \(v_2\), i.e. \(E[\exp(c_1) \mid v_2] = \exp(\rho v_2)\). With \(v_2\) standard normal the two differ by a factor of \(\sigma_{c_1}\), since \(\rho = \mathrm{Corr}(c_1, v_2)\,\sigma_{c_1}\). The distinction is not cosmetic: this chapter’s binary DGP sets \(c_1 = 0.2\,e_2\) with \(v_2 = e_2\), so the correlation is exactly \(1\) while the \(\rho\) governing the approximation below is \(0.2\). Reading \(\rho\) as a correlation is also what \(h(\cdot)\) forbids — it is \(\rho\) in the coefficient sense that makes \(\exp(\rho^2/2)\) the lognormal correction and puts \(\rho\) on the same scale as the probit index.

18.2.4.1 Terza’s FIML, and Stata’s etpoisson

Terza (1998) offers three estimators for this model, and it is worth knowing which is which. The first is full-information maximum likelihood: write the Poisson density for \(y_1\) and the probit for \(y_2\), let the two errors be jointly normal, and integrate the unobservable out. This is what Stata’s etpoisson fits, by Gauss–Hermite quadrature. The second is the two-step above — probit, then the exact mean by QMLE — which Terza calls a two-stage method of moments and describes as the relatively robust partially parametric option, since it needs the conditional mean and nothing else. The third sits between them, a nonlinear weighted least squares for the fully parametric case.

The distinction matters more in a panel than in a cross-section, for two reasons.

The first is what FIML assumes. Plain Poisson QMLE is consistent for the mean under any conditional distribution, because its first-order condition is the mean moment condition \(\sum \mathbf{x}(y - \exp(\mathbf{x}\boldsymbol\beta)) = 0\). That robustness is the reason fepois and ppmlhdfe are safe on count data that is nothing like Poisson. FIML’s score instead carries the Poisson density inside the integral, so getting the mean right is not enough, and vce(robust) repairs inference rather than the estimator.

The second is the fixed effects. etpoisson has no absorb() option and no panel version, so unit effects must enter as dummies — and here the Poisson rescue does not apply. The closed-form concentration \(\exp(\hat\alpha_i) = \sum_i y / \sum_i \exp(\mathbf{x}\boldsymbol\beta)\) works because the unit effect factors out of the Poisson likelihood; integrating over a normal unobservable destroys that separability, so the unit effects become ordinary nonlinear parameters and the incidental-parameter problem returns in full. With a few observations per unit that is fatal, and non-convergence is the usual symptom.

Replacing the dummies with unit means makes it estimable again, and in simulation it behaves reasonably: close to unbiased when its own distributional assumption holds, and off by a few percent when the count is overdispersed or zero-inflated. But it then rests on the Mundlak projection and joint normality and the Poisson shape, where the control function needs only the first two, and it converges unreliably — in our own runs it failed on five to twelve percent of samples, and the failures were the samples with the most dispersion, so discarding them is a selection. It also identifies the treatment effect from within- and between-unit variation, where absorbing the dummies uses within-unit variation only. The efficiency it buys is not worth that stack.

18.2.5 The generalized-residual CF as a first-order approximation

Taylor-expand \(h\) in \(\rho\) around \(\rho = 0\). Since \(\tfrac{d}{d\rho}\exp(\rho^2/2)\big|_{\rho=0} = 0\), only the \(\Phi\)-ratios contribute to first order:

\[\left.\frac{\partial}{\partial\rho}\frac{\Phi(\rho+\mathbf{z}\boldsymbol\pi_2)}{\Phi(\mathbf{z}\boldsymbol\pi_2)}\right|_{\rho=0} = \frac{\phi(\mathbf{z}\boldsymbol\pi_2)}{\Phi(\mathbf{z}\boldsymbol\pi_2)} = \lambda(\mathbf{z}\boldsymbol\pi_2), \qquad \left.\frac{\partial}{\partial\rho}\frac{1-\Phi(\rho+\mathbf{z}\boldsymbol\pi_2)}{1-\Phi(\mathbf{z}\boldsymbol\pi_2)}\right|_{\rho=0} = -\lambda(-\mathbf{z}\boldsymbol\pi_2) \tag{18.10}\]

where \(\lambda(\cdot) = \phi(\cdot)/\Phi(\cdot)\) is the inverse Mills ratio — exactly the gr2 computed via pdf(nd, η) / cdf(nd, η) below. Since \(h(y_2, \mathbf{z}\boldsymbol\pi_2, 0) = 1\),

\[h(y_2,\mathbf{z}\boldsymbol\pi_2,\rho) \approx 1 + \rho\cdot\underbrace{\bigl[y_2\lambda(\mathbf{z}\boldsymbol\pi_2) - (1-y_2)\lambda(-\mathbf{z}\boldsymbol\pi_2)\bigr]}_{=\ \hat r_2,\ \text{the generalized residual}} \tag{18.11}\]

and \(\exp(\mathbf{x}_1\boldsymbol\beta_1)\cdot h(\cdot) \approx \exp(\mathbf{x}_1\boldsymbol\beta_1)(1+\rho\hat r_2) \approx \exp(\mathbf{x}_1\boldsymbol\beta_1 + \rho\hat r_2)\) for small \(\rho\hat r_2\) — which is exactly “add \(\hat r_2\) linearly inside \(\exp(\cdot)\).”

So the generalized-residual CF is the first-order-in-\(\rho\) local approximation to Terza’s exact correction, with error \(O(\rho^2)\). A direct large-\(n\) check against the R companion chapter’s version of this DGP gives probability limits of \(0.839\) for the linear-residual CF and \(0.822\) for the generalized-residual CF against a true \(0.8\) — both biased, the latter roughly half the former, matching the theoretical ranking, and neither exact.

18.2.5.1 Why use the approximation rather than the exact correction?

Not because the exact correction is infeasible. It is worth being clear about this, because \(h(\cdot)\) looks forbidding: the \(\Phi\)-ratios and the \(\exp(\rho^2/2)\) term make Equation 18.8 genuinely nonlinear in \(\rho\), which suggests custom optimisation and no obvious way to combine it with the absorption that makes FE-Poisson tractable.

That reading is wrong, for a simple reason. Hold \(\rho\) fixed. Then \(\log h(y_2, \mathbf{z}\hat{\boldsymbol\pi}_2, \rho)\) is a variable — a known function of the data and the first-stage index — so the mean becomes

\[E(y_1 \mid \mathbf{z}, y_2, \alpha_i) = \exp\bigl(\mathbf{z}_1\boldsymbol\delta_1 + \alpha_1 y_2 + \alpha_i + \log h\bigr),\]

which is an ordinary Poisson mean with an offset. Any Poisson routine that takes one will fit it with the fixed effects absorbed exactly as usual — ppmlhdfe ..., absorb(id) offset(lnh) in Stata, fepois(..., offset = ~lnh) in R. Because \(h\) enters multiplicatively, the concentration that rescues FE-Poisson still goes through — \(\exp(\hat\alpha_i) = \sum_i y_1 / \sum_i \exp(\mathbf{x}\boldsymbol\beta)h\) — so there is no incidental-parameter problem. Only \(\rho\) itself is nonlinear, and it is one scalar: profile it on a grid and keep the value maximising the Poisson quasi-likelihood. This is Terza’s two-step, and it is what Wooldridge (2010) (sec. 18.5) prescribes for the cross-section — “we can use NLS or a quasi-MLE, such as the Poisson or gamma QMLE” — with the fixed effects added.

So the honest reasons for the approximation are convenience and the test, not feasibility. It runs as one command with no profile and no bootstrap; the coefficient on \(\hat r_2\) doubles as a Hausman-type \(t\)-test sitting next to the coefficient of interest, which is in fact how Wooldridge introduces it — as a variable-addition test of \(\rho_1 = 0\), noting that its \(\hat\rho_1\) “is not the estimate we obtain from Terza’s two-step method.” Joint nonlinear estimation is also not free of trouble: for the closely analogous binary-\(y_1\) case Wooldridge notes the iterations can be hard to converge, with \(\hat\rho_1\) drifting toward \(\pm 1\) (Wooldridge 2010, sec. 15.7.3), though a bounded profile grid makes that visible rather than fatal.

What matters is how large \(\rho\) is, and Equation 18.11 gives the rule: the error is \(O(\rho^2)\). In a Monte Carlo at 2,000 units with four observations each, the approximation is indistinguishable from the exact fit at \(|\rho| = 0.2\) — the value this chapter’s binary DGP uses — and the exact version is unbiased throughout.

The size of the error is not symmetric in \(\rho\), which is worth knowing before worrying about it. At \(\rho = +0.5\) the approximation overstates the effect by about nine percent and at \(\rho = +0.8\) by about thirty; at \(\rho = -0.5\) the bias is invisible and at \(\rho = -0.8\) it is about six percent downward. The asymmetry comes from \(h\) itself: a positive \(\rho\) scales the treated mean by a factor bounded by \(1/\Phi(\mathbf{z}\boldsymbol\pi_2)\), which grows without bound in the tail, while a negative \(\rho\) shrinks it smoothly toward zero. The \(O(\rho^2)\) rate describes the behaviour near zero; it does not put the same coefficient on both sides.

So: report the approximation when \(\hat\rho_1\) comes back small or negative, and fit the exact version alongside it when it is large and positive. Two cautions on the exact version. It is exact only under the joint normality of Equation 18.6 — when the unobservable depends on \(v_2\) nonlinearly, fitting \(h\) can do worse than truncating it, because the profile absorbs curvature \(h\) cannot represent. And its standard errors need a cluster bootstrap, since the first stage is plugged in and \(\rho\) is profiled.

18.2.5.2 Probit or logit — and why logit is simpler than it looks

Nothing in the local-approximation argument is specific to the normal. The generalized residual is defined for any correctly specified binary-response GLM (Gourieroux et al. 1987), and for a logit it collapses. With \(\Lambda(\cdot)\) the logistic CDF and \(\lambda(\eta) = \Lambda(\eta)[1-\Lambda(\eta)]\) its density,

\[\hat r_2 = y_2\frac{\lambda(\hat\eta)}{\Lambda(\hat\eta)} - (1-y_2)\frac{\lambda(\hat\eta)}{1-\Lambda(\hat\eta)} = y_2[1-\Lambda(\hat\eta)] - (1-y_2)\Lambda(\hat\eta) = y_2 - \Lambda(\hat\eta) = y_2 - \hat p_2, \tag{18.12}\]

i.e. the logit generalized residual is just the plain response residual \(y_2 - \hat p_2\). No Mills-ratio machinery: a logit first stage, then y2 .- predict(fit), is the generalized residual — not a shortcut to it.

Two limits. Terza’s exact correction is probit-specific: \(h(\cdot)\) comes from the truncated moments of the bivariate normal, and no analogue for a logistic \(v_2\) appears in Terza (1998), Wooldridge (2010) or Imbens and Wooldridge (2007). And the link-robustness that the GMM route enjoys — the fitted probability is a valid instrument even if the link is wrong, since any function of \(\mathbf{z}\) is a valid instrument (Imbens and Wooldridge 2007) — belongs to using the fitted value as an instrument, not as a residual. It does not carry over to either control function, where the second-stage conditional mean still has to be right.

18.2.6 A different family: multiplicative-moment GMM

The third route abandons the control function. Any function of the exogenous \(\mathbf{z}\) is a valid instrument, so the fitted probability \(\hat p_2\) can instrument \(y_2\) directly in a GMM problem that leaves the exponential mean untouched. Probit, logit and linear-probability fitted values are all equally valid for consistency; the choice affects only instrument strength. The moments are Mullahy’s multiplicative ones,

\[E[\mathbf{z}'(y_1 e^{-\mathbf{x}\boldsymbol{\beta}} - 1)] = 0, \tag{18.13}\]

which is the right choice when the unobserved heterogeneity sits inside \(\exp(\cdot)\), as it does for counts.

Two consequences follow, and they are the whole tradeoff. The validity of this route does not depend on the first-stage model being right — \(\hat p_2\) need only be a function of \(\mathbf{z}\). But it does not accommodate high-dimensional fixed effects: absorbing them would require within-group demeaning of the regressor and instrument matrices before solving nonlinear moment conditions, which is not straightforward. Stata’s ivpoisson similarly drops absorb().

18.3 The three estimators compared

Two questions separate these estimators, and everything else they share. Is it valid when \(y_2\) is binary, and can it absorb high-dimensional fixed effects?

Linear-residual CF Generalized-residual CF Multiplicative-moment GMM
How it works OLS residual \(\hat v_2\) in the index probit/logit generalized residual \(\hat r_2\) in the index fitted \(\hat p_2\) as an instrument
Valid for binary \(y_2\)? No — reduced-form assumption cannot hold Approximately — error \(O(\rho^2)\), derived Yes — consistent
First stage must be correct? moot yes no
Takes high-dimensional FE? yes yes no

No estimator gets both. The generalized-residual CF comes closest, buying fixed effects at the price of an approximation whose error is at least derived and shrinks with \(\rho\). The GMM route is the only genuinely consistent one, and pays for it by giving up fixed effects — and, as the simulations show, by roughly tripling the variance. The linear-residual CF is the fallback: for binary \(y_2\) it has no valid reduced form behind it, though it remains a fully valid exogeneity test and is accurate when \(\rho\) is small.

Common to all three: the second-stage exponential mean must be correctly specified, and the two control functions need bootstrapped standard errors because their residual is generated.

18.4 Simulation 1 — continuous endogenous variable

18.4.1 Data generating process

We simulate a panel where time (continuous, e.g. minutes on site) is endogenous: it is correlated with an unobserved factor e that also affects visits. The excluded instrument is phone (binary, e.g. phone-access indicator). Fixed effects are ad (20 groups) and female.

Code 
Random.seed!(42)
n_ad = 20; obs_per = 250; n = n_ad * obs_per

fe_ad  = randn(n_ad) .* 0.5
fe_grp = repeat(1:n_ad, inner = obs_per)
female = Int.(rand(n) .< 0.5)
phone  = Int.(rand(n) .< 0.4)
frfam  = rand(n)

# Common error drives endogeneity
e = randn(n)

# Continuous endogenous variable (linear first stage holds exactly)
time = 1.5 .* phone .+ 0.5 .* frfam .+ fe_ad[fe_grp] .+ e

# Poisson outcome — true coefficient on time = 0.8
true_α = 0.8
λ = exp.(0.5 .+ true_α .* time .+ 0.4 .* frfam .+
         fe_ad[fe_grp] .+ 0.3 .* female .+ 0.5 .* e)
visits = [rand(Poisson(λi)) for λi in λ]

df1 = DataFrame(visits=visits, time=time, phone=phone,
                frfam=frfam, female=female, ad=fe_grp)

println("n = $n   groups = $n_ad   visits mean = $(round(mean(visits),digits=2))")
n = 5000   groups = 20   visits mean = 15.93

18.4.2 Naive Poisson — the endogeneity bias

Ignoring the endogeneity of time gives a biased estimate of α₁:

Code 
naive1 = fepois(df1, @formula(visits ~ time + frfam + fe(ad) + fe(female)))
println("Naive (biased):")
show_regression_html(naive1; title="Naive Poisson — time is endogenous")
Naive (biased):

Naive Poisson — time is endogenous

Estimate Std. Error t value Pr(>|t|)
time 1.150 0.003 392.827 < 2e-16 ***
frfam 0.291 0.013 22.698 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

18.4.3 The linear-residual control function

Applying the framework above, the two-step procedure is:

  1. First stage — linear projection of time on phone, frfam, and fixed effects via feols. The residual \(\hat{v}_2 = \text{time} - \hat{\text{time}}\) satisfies \(E(\mathbf{z}'\hat{v}_2)=0\) by construction. A linear first stage is appropriate here because time is continuous — the condition that fails for binary \(y_2\) holds in this design.
  2. Second stage — Poisson QMLE with \(\hat{v}_2\) added as a regressor inside the exponential mean. Equation 18.2 holds exactly here because the DGP specifies \(0.5\,e\) linearly inside \(\exp(\cdot)\).

Note that unlike 2SLS in a linear model, \(\hat{v}_2\) does not replace \(y_2\) — both enter together. This is what makes the control function different from IV/GMM Poisson, confirmed numerically in Simulation 2.

Code 
# Step 1: linear projection
fs1 = feols(df1, @formula(time ~ phone + frfam + fe(ad) + fe(female)))
df1.v2hat = fs1.residuals

# Step 2: Poisson CF
m1_cf = fepois(df1, @formula(visits ~ time + v2hat + frfam + fe(ad) + fe(female)))
println("CF estimate (should recover α₁ = $(true_α)):")
show_regression_html(m1_cf; title="CF: linear projection first stage")
CF estimate (should recover α₁ = 0.8):

CF: linear projection first stage

Estimate Std. Error t value Pr(>|t|)
time 0.810 0.005 158.862 < 2e-16 ***
v2hat 0.489 0.006 76.972 < 2e-16 ***
frfam 0.405 0.013 31.664 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The coefficient on time recovers the true \(\alpha_1 = 0.8\). The v2hat coefficient tests exogeneity — its \(t\)-statistic is the Hausman-type test.

18.4.4 Bootstrap standard errors

The second-stage standard errors are incorrect because the first-stage estimation error propagates. Julia’s Threads.@threads parallelises the cluster bootstrap across all CPU cores:

Code 
function bootstrap_resample(df::DataFrame, B::Int; seed::Int=42)
    Random.seed!(seed)
    clusters = unique(df.ad)
    K = length(clusters)
    draws = [sample(clusters, K, replace=true) for _ in 1:B]
    return draws
end

function build_boot_df(df::DataFrame, drawn_clusters::Vector)
    frames = DataFrame[]
    for (id, c) in enumerate(drawn_clusters)
        sub = copy(df[df.ad .== c, :])
        sub.ad_boot = fill(id, nrow(sub))
        push!(frames, sub)
    end
    return vcat(frames...)
end

# Part 1 bootstrap: continuous time
function bootstrap_cf_continuous(df::DataFrame, B::Int=500; seed::Int=42)
    draws   = bootstrap_resample(df, B; seed=seed)
    results = zeros(B)
    Threads.@threads for b in 1:B
        df_b = build_boot_df(df, draws[b])
        fs   = feols(df_b,  @formula(time ~ phone + frfam + fe(ad_boot) + fe(female)))
        df_b.v2hat = fs.residuals
        ss   = fepois(df_b, @formula(visits ~ time + v2hat + frfam + fe(ad_boot) + fe(female)))
        idx  = findfirst(==("time"), coefnames(ss))
        results[b] = coef(ss)[idx]
    end
    return results
end

boot1 = bootstrap_cf_continuous(df1, 500)
naive_se1 = stderror(m1_cf)[findfirst(==("time"), coefnames(m1_cf))]
coef1     = coef(m1_cf)[findfirst(==("time"), coefnames(m1_cf))]

@printf("time coefficient:  %.4f  (true %.1f)\n", coef1, true_α)
@printf("Naive 2nd-stage SE: %.4f\n", naive_se1)
@printf("Bootstrap SE:       %.4f\n", std(boot1))
@printf("Bootstrap 95%% CI: [%.4f, %.4f]\n",
        quantile(boot1, 0.025), quantile(boot1, 0.975))
time coefficient:  0.8103  (true 0.8)
Naive 2nd-stage SE: 0.0051
Bootstrap SE:       0.0112
Bootstrap 95% CI: [0.7868, 0.8312]

The bootstrapped SE is larger than the naive second-stage SE, confirming that ignoring the first-stage uncertainty understates uncertainty.

18.5 Simulation 2 — binary endogenous variable

18.5.1 Data generating process

A clean probit DGP, with a fresh structural error e2 independent of Simulation 1. The endogeneity coefficient ρ = 0.2 is kept small so the first-order approximation stays tight and the estimators can be compared without that error dominating.

Code 
Random.seed!(123)

# Fresh structural error for Part 2 (independent of Part 1's e)
e2 = randn(n)

# Binary endogenous variable: clean probit latent index
# phone is the excluded instrument; no FE in latent index for clarity
time_hi = Int.(1.5 .* phone .+ 0.4 .* frfam .+ e2 .>= 0)

# Poisson outcome — true coefficient 0.8, endogeneity via ρ*e2 (ρ=0.2)
# Small ρ ensures the approximate CF assumption holds closely
ρ2 = 0.2
λ2 = exp.(0.5 .+ true_α .* time_hi .+ 0.4 .* frfam .+
          fe_ad[fe_grp] .+ 0.3 .* female .+ ρ2 .* e2)
visits2 = [rand(Poisson(λ2i)) for λ2i in λ2]

df2 = DataFrame(visits=visits2, time_hi=time_hi, phone=phone,
                frfam=frfam, female=female, ad=fe_grp, e2=e2)

println("time_hi mean:      $(round(mean(time_hi), digits=3))")
println("Cor(time_hi, e2):  $(round(cor(time_hi, e2), digits=3))  (endogeneity)")
println("visits mean:       $(round(mean(visits2), digits=2))")
time_hi mean:      0.727
Cor(time_hi, e2):  0.613  (endogeneity)
visits mean:       4.78

18.5.2 Linear-residual CF

The estimator with no valid reduced-form assumption here, included to see how far wrong it goes.

Code 
# Step 1: linear projection of binary time_hi
fs_a = feols(df2, @formula(time_hi ~ phone + frfam + fe(ad) + fe(female)))
df2.v2hat = fs_a.residuals

# Step 2: Poisson CF
m_a = fepois(df2, @formula(visits ~ time_hi + v2hat + frfam + fe(ad) + fe(female)))
println("Linear-residual CF:")
show_regression_html(m_a; title="Linear-residual CF (OLS residual)")
Linear-residual CF:

Linear-residual CF (OLS residual)

Estimate Std. Error t value Pr(>|t|)
time_hi 0.813 0.036 22.719 < 2e-16 ***
v2hat 0.330 0.035 9.396 < 2e-16 ***
frfam 0.413 0.023 17.921 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

18.5.3 Generalized-residual CF

Code 
# Step 1: probit first stage
fs_b = feprobit(df2, @formula(time_hi ~ phone + frfam + fe(ad) + fe(female)))

# Generalized residual (Gourieroux et al.): φ(η)/Φ(η) if y=1 (inverse Mills ratio),
# -φ(η)/(1-Φ(η)) if y=0 (the reverse Mills ratio).
η  = fs_b.eta
nd = Normal()
df2.gr2 = @. ifelse(df2.time_hi == 1,
    pdf(nd, η) / cdf(nd, η),
    -pdf(nd, η) / (1.0 - cdf(nd, η)))

# Step 2: Poisson CF with generalized residual
m_b = fepois(df2, @formula(visits ~ time_hi + gr2 + frfam + fe(ad) + fe(female)))
println("Generalized-residual CF (probit):")
show_regression_html(m_b; title="Generalized-residual CF (probit)")
Generalized-residual CF (probit):

Generalized-residual CF (probit)

Estimate Std. Error t value Pr(>|t|)
time_hi 0.764 0.039 19.549 < 2e-16 ***
gr2 0.230 0.024 9.721 < 2e-16 ***
frfam 0.415 0.023 17.969 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

18.5.4 Probit or logit

Confirming Equation 18.12 numerically — the logit generalized residual is time_hi .- phat:

Code 
fs_b_logit = felogit(df2, @formula(time_hi ~ phone + frfam + fe(ad) + fe(female)))
phat_logit = fs_b_logit.mu
df2.gr2_logit = df2.time_hi .- phat_logit    # algebraically = generalized residual

m_b_logit = fepois(df2, @formula(visits ~ time_hi + gr2_logit + frfam + fe(ad) + fe(female)))
println("Generalized-residual CF with a logit first stage:")
show_regression_html(m_b_logit; title="Generalized-residual CF (logit)")

coef_b_probit = coef(m_b)[findfirst(==("time_hi"), coefnames(m_b))]
coef_b_logit  = coef(m_b_logit)[findfirst(==("time_hi"), coefnames(m_b_logit))]
@printf("Probit first stage: alpha1 = %.4f\n", coef_b_probit)
@printf("Logit  first stage: alpha1 = %.4f\n", coef_b_logit)
Generalized-residual CF with a logit first stage:

Generalized-residual CF (logit)

Estimate Std. Error t value Pr(>|t|)
time_hi 0.800 0.036 22.145 < 2e-16 ***
gr2_logit 0.346 0.036 9.698 < 2e-16 ***
frfam 0.413 0.023 17.878 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1		 
Probit first stage: alpha1 = 0.7636
Logit  first stage: alpha1 = 0.8000

18.5.5 The dummies in the first stage need a long panel

Both first stages above put the group dummies inside the probit or logit, which is safe here because each group has 250 observations. Applied panel work usually has the opposite shape: many units, a handful of periods each. A logit with unit dummies is badly biased in that case, and the size of the bias depends on observations per unit, not on how many units there are. Against a true slope of 1 it returns about 2.0 at two observations per unit, 1.37 at four, 1.15 at six, 1.04 at twenty-five, and essentially the truth by a hundred. A panel of a few thousand units with four or five observations each sits in the worst part of that range.

The problem is specific to the nonlinear first stage. Poisson with unit dummies is consistent with the number of periods fixed, and so is a linear first stage.

The fix is the Mundlak device: drop the unit dummies from the first stage and add unit means of the instrument and the time-varying controls instead. The instrument’s own mean has to be among them, since the unit effect is being modelled as a function of those means and the instrument is one of the regressors it must be purged of. Conditional logit is also consistent with the number of periods fixed, but it eliminates the unit effects rather than estimating them, so it leaves no fitted index from which to build the generalized residual; the Mundlak version keeps the index. The construction is worked through in IV Poisson with many fixed effects.

Lin and Wooldridge (2019) treat this case directly for the exponential model with unobserved effects, and land in the same place. Their Procedure 3 adds a first-stage residual to FE Poisson and tests \(\rho_1 = 0\); for a binary endogenous variable they write that one can “use as a control function the generalized residuals … where \(\mathbf{w}_{it} = (1, \mathbf{z}_{it}, \bar{\mathbf{z}}_i)\)” — the time-varying exogenous variables and their unit means, the instrument’s mean included. They add the caveat in the same breath: the Mundlak device “might work reasonably well, but neither might be flexible enough,” and they “leave investigations into the quality of CF approximations in discrete cases to future research.” So the recommendation below is theirs as well as ours, and the open question they name is the one worth remembering: how good the approximation is when \(y_2\) is discrete has not been settled, and Equation 18.11 says the error grows with \(\rho\).

18.5.6 Multiplicative-moment GMM

We implement IV/GMM Poisson directly, keeping time_hi as the regressor and using probit fitted probabilities as the instrument. Both moment forms are coded: the additive moments are what “IV Poisson” often means in practice (e.g. Stata’s ivpoisson gmm default), but because this chapter’s DGP puts the unobservable inside \(\exp(\cdot)\), only the multiplicative moments of Mullahy (1997) are consistent here.

Code 
"""
    ivpois(y, X, Z) -> (β, se)

IV/GMM Poisson (just-identified). X = [endog, exog], Z = [instrument, exog].
Damped IRLS: intercept initialised from log(mean(y)), half-steps for stability.
"""
function ivpois(y::Vector{Float64}, X::Matrix{Float64}, Z::Matrix{Float64};
                maxiter=500, tol=1e-8, step=0.5)
    p  = size(X, 2)
    β  = zeros(p)
    β[p] = log(mean(y))    # intercept from marginal mean — prevents divergence

    for _ in 1:maxiter
        μ = exp.(clamp.(X * β, -30.0, 30.0))
        H = Z' * (μ .* X)
        g = Z' * (y .- μ)
        Δ = try H \ g catch; pinv(H) * g end
        β .+= step .* Δ
        maximum(abs.(step .* Δ)) < tol && break
    end

    μ = exp.(clamp.(X * β, -30.0, 30.0))
    H = Z' * (μ .* X)
    B = Z' * ((y .- μ).^2 .* Z)       # robust sandwich meat
    V = try inv(H) * B * inv(H)' catch; pinv(H) * B * pinv(H)' end
    return β, sqrt.(max.(0.0, diag(V)))
end

"""
    ivpois_mult(y, X, Z) -> β

IV/GMM Poisson with MULTIPLICATIVE (Mullahy 1997) moments
Z'(y .* exp(-Xβ) .- 1) = 0 — consistent when the unobservable enters
inside the exponential mean (multiplicative error).
"""
function ivpois_mult(y::Vector{Float64}, X::Matrix{Float64}, Z::Matrix{Float64};
                     maxiter=1000, tol=1e-10, step=0.5)
    p  = size(X, 2)
    β  = zeros(p)
    β[p] = log(mean(y))
    for _ in 1:maxiter
        w = y .* exp.(clamp.(-X * β, -30.0, 30.0))
        g = Z' * (w .- 1)
        H = Z' * (w .* X)
        Δ = try H \ g catch; pinv(H) * g end
        β .+= step .* Δ
        maximum(abs.(step .* Δ)) < tol && break
    end
    return β
end

# Probit (without FE, to avoid collinearity with instrument matrix)
fs_c  = feprobit(df2, @formula(time_hi ~ phone + frfam))
phat  = fs_c.mu

# Regressor matrix X: [time_hi, frfam, intercept]
# Instrument matrix Z: [phat,    frfam, intercept]  ← phat replaces time_hi
X_c = hcat(Float64.(df2.time_hi), Float64.(df2.frfam), ones(n))
Z_c = hcat(Float64.(phat),        Float64.(df2.frfam), ones(n))

β_c, se_c = ivpois(Float64.(df2.visits), X_c, Z_c)
β_cm      = ivpois_mult(Float64.(df2.visits), X_c, Z_c)
@printf("\nMultiplicative-moment GMM (no FE, for illustration)\n")
@printf("  additive moments:       time_hi = %7.4f   SE = %7.4f\n",
        β_c[1], se_c[1])
@printf("  multiplicative moments: time_hi = %7.4f   (true %.1f)\n",
        β_cm[1], true_α)

Multiplicative-moment GMM (no FE, for illustration)
  additive moments:       time_hi =  0.5273   SE =  0.0706
  multiplicative moments: time_hi =  0.7260   (true 0.8)

The additive-moment estimate is visibly biased — the DGP’s unobservable \(\rho\,e_2\) sits inside \(\exp(\cdot)\), so the additive moments are misspecified (a 30-replication check gives an additive-moment mean of about \(0.61\) against a multiplicative-moment mean of about \(0.79\), with true value \(0.8\)). The multiplicative moments repair it.

18.5.7 Control function and GMM are not equivalent in Poisson

In a linear model, using \(\hat{y}_2\) as an instrument (2SLS) is numerically identical to adding \(\hat{v}_2 = y_2 - \hat{y}_2\) as a regressor, by Frisch–Waugh–Lovell. In Poisson this equivalence breaks, because the control function puts \(\hat{v}_2\) inside \(\exp(\cdot)\) while the GMM route uses \(\hat{y}_2\) in the moment conditions without modifying the exponential mean.

Code 
# Use linear fitted values as instruments (GMM route, linear first stage)
ŷ2  = Float64.(df2.time_hi) .- Float64.(fs_a.residuals)
X_l = hcat(Float64.(df2.time_hi), Float64.(df2.frfam), ones(n))
Z_l = hcat(ŷ2,                    Float64.(df2.frfam), ones(n))

β_l, _ = ivpois(Float64.(df2.visits), X_l, Z_l)

coef_a = coef(m_a)[findfirst(==("time_hi"), coefnames(m_a))]
@printf("Linear-residual CF (residual as regressor, with FE): α̂₁ = %.4f\n", coef_a)
@printf("GMM, additive moments (linear fitted values, no FE): α̂₁ = %.4f\n", β_l[1])
@printf("GMM, additive moments (probit fitted values, no FE): α̂₁ = %.4f\n", β_c[1])
@printf("GMM, multiplicative moments (no FE):                 α̂₁ = %.4f\n", β_cm[1])
@printf("True α₁ = %.1f\n", true_α)
Linear-residual CF (residual as regressor, with FE): α̂₁ = 0.8128
GMM, additive moments (linear fitted values, no FE): α̂₁ = 0.6194
GMM, additive moments (probit fitted values, no FE): α̂₁ = 0.5273
GMM, multiplicative moments (no FE):                 α̂₁ = 0.7260
True α₁ = 0.8

The estimates differ despite using the same first-stage linear projection, confirming that 2SLS \(\neq\) CF for Poisson — and, within the GMM route, that the moment form matters as much as the instrument.

18.5.8 Bootstrap for the binary case

Code 
# Part 2 bootstrap: binary time_hi
function bootstrap_cf_binary(df::DataFrame, B::Int=500; seed::Int=42)
    draws   = bootstrap_resample(df, B; seed=seed)
    results = zeros(B)
    Threads.@threads for b in 1:B
        df_b = build_boot_df(df, draws[b])
        fs   = feols(df_b,  @formula(time_hi ~ phone + frfam + fe(ad_boot) + fe(female)))
        df_b.v2hat = fs.residuals
        ss   = fepois(df_b, @formula(visits ~ time_hi + v2hat + frfam + fe(ad_boot) + fe(female)))
        idx  = findfirst(==("time_hi"), coefnames(ss))
        results[b] = coef(ss)[idx]
    end
    return results
end

boot2 = bootstrap_cf_binary(df2, 500)
naive_se2 = stderror(m_a)[findfirst(==("time_hi"), coefnames(m_a))]

@printf("Linear-residual CF bootstrap (binary EEV, 500 reps, cluster by ad)\n")
@printf("  Point estimate:     %.4f  (true %.1f)\n", coef_a, true_α)
@printf("  Naive 2nd-stage SE: %.4f\n", naive_se2)
@printf("  Bootstrap SE:       %.4f\n", std(boot2))
@printf("  Bootstrap 95%% CI:  [%.4f, %.4f]\n",
        quantile(boot2, 0.025), quantile(boot2, 0.975))
Linear-residual CF bootstrap (binary EEV, 500 reps, cluster by ad)
  Point estimate:     0.8128  (true 0.8)
  Naive 2nd-stage SE: 0.0358
  Bootstrap SE:       0.0343
  Bootstrap 95% CI:  [0.7380, 0.8766]

18.5.9 All methods side by side

Code 
coef_naive2 = coef(
    fepois(df2, @formula(visits ~ time_hi + frfam + fe(ad) + fe(female)))
)[findfirst(==("time_hi"),
    coefnames(fepois(df2, @formula(visits ~ time_hi + frfam + fe(ad) + fe(female)))))]
coef_b = coef(m_b)[findfirst(==("time_hi"), coefnames(m_b))]

methods = ["True value", "Naive Poisson (binary y₂)",
           "Linear-residual CF (with FE)", "Generalized-residual CF (with FE)",
           "GMM, additive moments (no FE)",
           "GMM, multiplicative moments (no FE)"]
ests    = [true_α, coef_naive2, coef_a, coef_b, β_c[1], β_cm[1]]

println("\n── Binary EEV: coefficient on time_hi ─────────────────────")
for (m, e) in zip(methods, ests)
    @printf("  %-42s  α̂₁ = %6.4f\n", m, e)
end

── Binary EEV: coefficient on time_hi ─────────────────────
  True value                                  α̂₁ = 0.8000
  Naive Poisson (binary y₂)                   α̂₁ = 1.0852
  Linear-residual CF (with FE)                α̂₁ = 0.8128
  Generalized-residual CF (with FE)           α̂₁ = 0.7636
  GMM, additive moments (no FE)               α̂₁ = 0.5273
  GMM, multiplicative moments (no FE)         α̂₁ = 0.7260

The two control functions recover the truth closely — with small \(\rho = 0.2\) the approximate CF mean assumption is mild (the linear-residual CF has no valid first-stage assumption for binary \(y_2\); the generalized-residual CF is correct only to first order in \(\rho\), per the derivation above). Additive-moment IV Poisson does not: the unobservable enters inside \(\exp(\cdot)\), so its moments are misspecified (Mullahy 1997).

The multiplicative-moment version repairs it, and in one specific sense does more than repair it: the ad and female effects it omits are independent multiplicative heterogeneity, absorbed by the intercept under the multiplicative moments (\(\beta_0 \to \beta_0 + \log E[e^{c}]\)), so dropping the fixed effects costs no slope bias at all. Mullahy’s moments are consistent for this DGP, where the control functions are only approximately unbiased — \(0.8015\) at \(n=10^6\) (sd \(0.0013\)), against plims of \(0.839\) and \(0.822\).

What consistency does not buy is precision, and that is why the table above can look as though it says the opposite. Without the fixed effects the multiplicative-moment estimator is roughly three times as noisy as either control function (sd \(0.068\) against \(0.024\) and \(0.026\), over 10 seeds at this chapter’s \(n=5000\)), so a single draw can leave the only consistent estimator of the three further from \(0.8\) than either CF estimator. At this sample size its RMSE is in fact the worst of the three (\(0.065\), against \(0.039\) and \(0.031\)); the ranking reverses as \(n\) grows. Read the comparison as one draw from that tradeoff, not as a verdict.

18.6 Practical guidance

  • With fixed effects and a binary \(y_2\), use the generalized-residual CF. Among the estimators that drop into FE-Poisson unmodified it is the theoretically motivated choice, and its error is derived, not assumed: \(O(\rho^2)\), shrinking as the endogeneity shrinks. Put the unit dummies in its first stage only when each unit has many observations; with the short panels usual in applied work, use the Mundlak means instead (§ above).

  • Don’t reach for etpoisson in a panel. It fits Terza’s FIML, which is the exact model, but it needs the Poisson distribution rather than just the mean, it has no absorb(), and unit dummies inside its quadrature likelihood bring back the incidental-parameter problem that FE-Poisson is normally free of (§ above).

  • Multiplicative-moment GMM trades bias for variance — decide which binds. Mullahy’s moments are consistent here, not merely accurate to first order, whenever the multiplicative heterogeneity is independent of \(\mathbf{z}\). But consistency is not accuracy at a given sample size: over 10 seeds at \(n=5000\) it is nearly unbiased (mean \(0.7985\)) and yet almost three times as noisy (sd \(0.068\) against \(0.026\)), giving it the worst RMSE of the three. Prefer it once \(n\) is large enough that variance is not the binding constraint — and note that it does not accommodate high-dimensional fixed effects either way.

  • Use the multiplicative moments, not the additive ones. When the unobservable sits inside \(\exp(\cdot)\) — the natural case for counts — the additive moments are misspecified, and the simulation above shows the damage.

  • The linear-residual CF is a fallback, not a first choice. It lands close to the truth at \(\rho = 0.2\) here, but for binary \(y_2\) it has no valid reduced-form assumption behind it. Do not lean on it when \(\rho\) is not known to be small — and note that \(\rho\) is estimated by the coefficient on the control function, so this is checkable rather than a matter of faith.

  • It remains fully valid as an exogeneity test. Under \(H_0: \rho = 0\) none of the reduced-form assumptions are needed (Wooldridge 2010, 744). Use it to test; do not use it to correct.

  • Test, but do not let the test pick the model. The obvious move is to run the \(t\)-test on the control function and, if it comes back insignificant, drop the residual and report the plain FE-Poisson. That is a pretest estimator and it is worse than either model it selects between. Guggenberger (2010) shows for the linear case that a Hausman pretest followed by OLS-or-2SLS, with the second-stage \(t\)-test read at its nominal level, has asymptotic size one. The mechanism is local endogeneity: let \(\rho\) drift to zero at rate \(n^{-1/2}\) and the pretest usually fails to reject, while the bias of the estimator it falls back on is the same order as that estimator’s standard error. No such theorem exists for the Poisson CF, but the mechanism needs only a low-powered pretest and a fallback whose bias does not vanish faster than its standard error, and both are here. Report the naive and the control-function columns side by side.

  • Probit or logit both work for the generalized-residual CF. Only Terza’s exact correction needs normality.

  • Always bootstrap. In every control-function specification above the second-stage standard errors are wrong, because the residual is generated.

18.7 Why Julia is faster

The parallel bootstrap above uses Threads.@threads which runs each of the 500 replications simultaneously across all CPU cores. A single replication calls feols + fepois, both implemented with the fast within-demeaning algorithm from Panelest.jl. On a machine with 8 cores this is roughly 6–8× faster than Stata’s serial bootstrap command, and the speed advantage grows with replication count and sample size.