# IV in Poisson with Fixed Effects
```{r}
#| include: false
library(fixest)
```
Poisson regression specifies the conditional mean as an exponential function:
$$E(y \mid \mathbf{x}) = \exp(\mathbf{x}\boldsymbol{\beta})$$
As a quasi-maximum likelihood estimator (QMLE), Poisson regression requires only that the conditional mean is correctly specified; the true count distribution does not need to follow a Poisson distribution, and standard errors are easily made robust to arbitrary overdispersion or heteroskedasticity.
However, when an explanatory variable is **endogenous**—correlated with unobserved omitted factors—standard linear two-stage least squares (2SLS) intuition fails. Because the conditional mean is nonlinear, unobserved confounders enter multiplicatively inside the exponential function.
This chapter covers:
1. Why linear 2SLS intuition fails in nonlinear exponential models.
2. The control function approach for continuous vs. binary endogenous variables.
3. Terza's exact correction [@terza-1998] and its profile-likelihood implementation.
4. Generalized residuals as a first-order local approximation, logit response residuals, and the Mundlak device for panel fixed effects [@lin-wooldridge-2019].
5. Why applied econometrics generally recommends the generalized residual approach over the exact Terza correction.
6. Multiplicative-moment GMM [@mullahy-1997] and quasi-differencing for panel fixed effects.
7. Why the exogeneity test should not be used to choose between the naive and control-function models [@guggenberger-2010].
8. Step-by-step implementation demos in R (`fixest`) and Stata (`ppmlhdfe`).
---
## 1. The Model and Setup
Let $y_1$ denote the count outcome, $y_2$ the endogenous explanatory variable, and $\mathbf{z} = (\mathbf{z}_1, \mathbf{z}_2)$ the vector of exogenous variables ($\mathbf{z}_1$ included in the structural equation, $\mathbf{z}_2$ excluded instruments).
The structural conditional mean is:
$$E(y_1 \mid \mathbf{z}_1, y_2, c_1) = \exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + c_1)$$ {#eq-poisson-iv-1}
where $c_1$ is an unobserved omitted variable correlated with $y_2$. Our primary estimand of interest is $\alpha_1$, the structural semi-elasticity of $y_1$ with respect to $y_2$.
### Subscripts Index Equations
Following the simultaneous-equations convention of @wooldridge-2010-crosssection (ch. 18), **subscripts index equations rather than individual variables**:
- **Equation 1 (Structural equation):** outcome $y_1$, parameters $(\boldsymbol{\delta}_1, \alpha_1)$, unobserved error $c_1$.
- **Equation 2 (Reduced-form equation):** outcome $y_2$, coefficients $\boldsymbol{\pi}_2$, reduced-form error $v_2$.
| Role | Linear IV Notation | Nonlinear Notation Here | Description |
|:---|:---:|:---:|:---|
| Outcome (the count) | $Y$ | $y_1$ | Dependent variable in structural equation |
| Endogenous regressor | $W$ or $D$ | $y_2$ | Explanatory variable correlated with $c_1$ |
| Excluded instrument | $Z$ | $\mathbf{z}_2$ | Exogenous instrument excluded from Eq. 1 |
| Included exogenous | $X$ | $\mathbf{z}_1$ | Exogenous covariates included in Eq. 1 |
| Reduced-form error | $V$ | $v_2$ | Error from projecting $y_2$ on $\mathbf{z}$ |
| Structural unobservable | $U$ | $c_1$ | Confounder correlated with $y_2$ |
| Endogeneity parameter | $\rho$ | $\rho_1$ | Projection coefficient: $c_1 = \rho_1 v_2 + e_1$ |
Two standard deviations from Wooldridge's notation are maintained for clarity:
1. Wooldridge writes the coefficient on $y_2$ as $\gamma_1$; we write $\alpha_1$ to align with standard treatment effect notation.
2. Where Wooldridge writes $\mathbf{x}_1\boldsymbol{\beta}_1$, we write $\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2$.
### Fixed Effects in the Exponential Model
Everything below absorbs fixed effects without comment, and it is worth saying once why that is allowed. Fixed-effects Poisson has no incidental parameters problem [@wooldridge-1999-multiplicative; @lin-wooldridge-2019]. With $N \to \infty$ and $T$ fixed, the unit effects enter multiplicatively, factor out of the score, and are conditioned out cleanly. The estimator of $\boldsymbol{\beta}$ is consistent even though the number of unit parameters grows with the sample. Probit and logit have no such property, which is why the first stage of @sec-poisson-iv-gr needs the Mundlak device while the second stage does not.
The software follows from this. In Stata, `xtpoisson, fe` fits the conditional model, but the community standard is `ppmlhdfe` [@correia-guimaraes-zylkin-2020]: it absorbs several dimensions of fixed effects at once, it is fast, it detects separated and singleton observations, and it recovers predicted values with the absorbed effects put back. In R, `fixest::fepois()` does the same, with `fixest::feols()` for the linear first stage.
Stata's `ivpoisson` handles endogenous regressors, but it cannot absorb high-dimensional fixed effects. That is the practical reason the panel literature routes endogeneity through a control function rather than through a packaged IV command: the control function leaves the second stage as an ordinary FE Poisson regression, which `ppmlhdfe` and `fepois` already do well.
---
## 2. Theoretical Foundations: Control Functions vs. GMM
### 2.1 The Linear Second Stage: Control Function *is* 2SLS
To see where the nonlinear complication begins, consider first the linear model:
$$y_1 = \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \varepsilon_1, \qquad E(\mathbf{z}'\varepsilon_1) = \mathbf{0}$$
Two-stage least squares (2SLS) replaces $y_2$ with its linear projection $\hat{y}_2 = \mathbf{z}\hat{\boldsymbol{\pi}}_2$:
$$y_1 = \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 \hat{y}_2 + (\varepsilon_1 + \alpha_1 \hat{v}_2)$$
where $\hat{v}_2 = y_2 - \hat{y}_2$ is the OLS residual.
Alternatively, the **control function (CF)** approach decomposes the structural error $\varepsilon_1$ via a linear projection on the first-stage residual $v_2$:
$$\varepsilon_1 = \rho_1 v_2 + \eta_1$$
where $\rho_1 = \frac{\mathrm{Cov}(\varepsilon_1, v_2)}{\mathrm{Var}(v_2)}$. By the properties of linear projection, $E(v_2 \eta_1) = 0$ and $E(\mathbf{z}'\eta_1) = E(\mathbf{z}'\varepsilon_1) - \rho_1 E(\mathbf{z}'v_2) = \mathbf{0}$.
Plugging this into the structural equation yields the augmented regression:
$$y_1 = \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 \hat{v}_2 + \eta_1$$
By the **Frisch–Waugh–Lovell (FWL) theorem**, OLS estimation of this augmented regression produces an estimate $\hat{\alpha}_1$ that is **numerically identical to 2SLS**. No additional distributional assumptions are required beyond the standard orthogonality condition $E(\mathbf{z}'\varepsilon_1) = \mathbf{0}$.
### 2.2 The Nonlinear Second Stage: Control Function Requires Independence
When the conditional mean is exponential, this automatic numerical equivalence breaks down.
Substituting the error decomposition $c_1 = \rho_1 v_2 + e_1$ into @eq-poisson-iv-1 yields:
$$\begin{aligned}
E(y_1 \mid \mathbf{z}, y_2, v_2) &= E\left[ \exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 v_2 + e_1) \mid \mathbf{z}, y_2, v_2 \right] \\
&= \exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 v_2) \cdot E\left[ \exp(e_1) \mid \mathbf{z}, y_2, v_2 \right]
\end{aligned}$$ {#eq-poisson-iv-2}
For the augmented index $\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 v_2$ to correctly specify the conditional mean in a Poisson second stage, the remaining expectation factor must be **strictly constant**:
$$E\left[ \exp(e_1) \mid \mathbf{z}, y_2, v_2 \right] = \text{constant}$$
This is why @wooldridge-2010-crosssection (sec. 18.5) normalizes $E[\exp(e_1)] = 1$ rather than $E(e_1) = 0$.
This condition is a requirement on the **moment-generating function** of $e_1$—a condition across all moments simultaneously, which requires that $e_1$ is **fully independent** of $(\mathbf{z}, v_2)$, and in turn that the reduced-form error $v_2$ is **independent** of $\mathbf{z}$ ($v_2 \perp \mathbf{z}$).
| Property | Linear Model (OLS / 2SLS) | Nonlinear Model (Poisson CF) |
|:---|:---|:---|
| **Second stage** | $y_1 = \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 \hat{v}_2$ | $y_1 \sim \text{Poisson}(\exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 \hat{v}_2))$ |
| **CF vs. 2SLS/GMM** | **Identical** (by FWL theorem) | **Different** (FWL does not hold) |
| **Required assumption** | Uncorrelatedness: $E(\mathbf{z}'\varepsilon_1) = \mathbf{0}$ | **Full Independence:** $v_2 \perp \mathbf{z}$ and $e_1 \perp (\mathbf{z}, v_2)$ |
| **Second-stage standard errors** | Valid under $H_0: \rho_1 = 0$ | **Invalid** (must bootstrap generated regressor) |
| **Endogeneity test** | $t$-statistic on $\hat{v}_2 = 0$ (DWH test) | $t$-statistic on $\hat{v}_2 = 0$ (Variable addition test) |
---
## 3. Estimators for Endogenous Regressors
### 3.1 Continuous $y_2$: The Linear Reduced Form Holds
When $y_2$ is continuous (e.g., hours worked, prices, expenditures), a linear reduced form with additive errors is natural:
$$y_2 = \mathbf{z}\boldsymbol{\pi}_2 + v_2, \qquad c_1 = \rho_1 v_2 + e_1$$ {#eq-poisson-iv-3}
with $v_2 \perp \mathbf{z}$ and $e_1 \perp (\mathbf{z}, v_2)$. These are two separate restrictions and it is worth keeping them apart. Joint normality of $(c_1, v_2)$ delivers the second: $e_1$ is the projection residual, and for jointly normal variables uncorrelated implies independent. It says nothing about the first. $v_2 \perp \mathbf{z}$ is a restriction on the reduced form — the conditional distribution of $v_2$ must not vary with $\mathbf{z}$ — and a linear equation with a homoskedastic normal additive error satisfies it. Both hold here; neither is implied by the other. The structural conditional mean becomes:
$$E(y_1 \mid \mathbf{z}_1, y_2, v_2) = \exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 v_2)$$ {#eq-poisson-iv-4}
**Estimation Procedure:**
1. **Step 1 (First Stage):** Estimate the linear regression $y_2 = \mathbf{z}\boldsymbol{\pi}_2 + v_2$ by OLS (or with fixed effects absorbed) and compute residuals $\hat{v}_2 = y_2 - \mathbf{z}\hat{\boldsymbol{\pi}}_2$.
2. **Step 2 (Second Stage):** Estimate Poisson QMLE with $y_2$ and $\hat{v}_2$ entering simultaneously:
$$y_1 \sim \text{Poisson}(\exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 \hat{v}_2))$$
3. **Inference:** Standard errors must be bootstrapped (e.g., clustered by panel unit) to account for the first-stage estimation uncertainty of $\hat{v}_2$.
This is @lin-wooldridge-2019's Procedure 3. In a panel, step 1 absorbs the unit effects (`reghdfe` in Stata, `feols` in R) and step 2 absorbs them again (`ppmlhdfe`, `fepois`).
#### The First Stage May Absorb or Average: Both Give the Same Answer
There are two ways to handle unit effects in a linear first stage, and it does not matter which we pick. Regressing $y_{it2}$ on $\mathbf{z}_{it}$ and unit dummies gives the fixed-effects residual $\hat{v}_{it2}^{FE}$. Regressing $y_{it2}$ on $\mathbf{z}_{it}$ and the unit means $\bar{\mathbf{z}}_i$ — the Mundlak device — gives the Mundlak residual $\hat{v}_{it2}^{M}$. Substituting either one into the second-stage FE Poisson gives **numerically identical** $\hat{\alpha}_1$ and $\hat{\rho}_1$ [@lin-wooldridge-2019, app. A.1–A.2].
The reason is short. Mundlak's result is that the coefficient on $\mathbf{z}_{it}$ in the second regression equals the within estimator $\hat{\boldsymbol{\pi}}_2^{FE}$, so the two residuals share a slope and differ only in what is subtracted off unit by unit:
$$\hat{v}_{it2}^{FE} = \hat{v}_{it2}^{M} + g_i$$
for some unit-specific constant $g_i$. Put that into the second-stage index:
$$\mathbf{z}_{1,it}\boldsymbol{\delta}_1 + \alpha_1 y_{it2} + \rho_1 \hat{v}_{it2}^{FE} + c_i = \mathbf{z}_{1,it}\boldsymbol{\delta}_1 + \alpha_1 y_{it2} + \rho_1 \hat{v}_{it2}^{M} + (c_i + \rho_1 g_i)$$
The swap moves $\rho_1 g_i$ into the unit effect, and the unit effect is a free parameter that FE Poisson absorbs without restriction. The model is reparameterized, not changed. Nothing reaches $\hat{\alpha}_1$ or $\hat{\rho}_1$.
This holds only when the first stage averages or absorbs on the same dimension the second stage absorbs. It is also specific to the linear first stage: for a binary $y_2$ the first stage is nonlinear, the residual is not shifted by a constant, and the choice between dummies and Mundlak means becomes a real one (@sec-poisson-iv-gr).
---
### 3.2 Binary $y_2$: Why the Linear Reduced Form Fails
When $y_2 \in \{0, 1\}$ is binary, the linear reduced-form residual $v_2 = y_2 - \mathbf{z}\boldsymbol{\pi}_2$ takes exactly two values conditional on $\mathbf{z}$:
$$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}$$ {#eq-poisson-iv-5}
The support $\{-\mathbf{z}\boldsymbol{\pi}_2, 1 - \mathbf{z}\boldsymbol{\pi}_2\}$, the probabilities, and the conditional variance $\mathrm{Var}(v_2 \mid \mathbf{z}) = p(\mathbf{z})[1 - p(\mathbf{z})]$ are all **direct functions of $\mathbf{z}$**.
As @wooldridge-2010-crosssection [p. 746] emphasizes:
> "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}$.)"
A distribution whose support moves with $\mathbf{z}$ cannot be independent of $\mathbf{z}$ (unless $\boldsymbol{\pi}_2 = \mathbf{0}$, meaning there is no first stage). Therefore, $v_2 \perp \mathbf{z}$ fails mathematically for any discrete endogenous variable.
#### Where the Independence Assumption Belongs
The issue is not that the observed residual is discrete. The generalized residual of Section 3.4 is also two-valued given $\mathbf{z}$.
The critical distinction is **which object independence is asserted about**:
- The linear CF asserts independence of the **observed discrete residual** $v_2$, which is impossible.
- Terza's threshold model asserts independence of the **continuous latent probit error** $v_2 \mid \mathbf{z} \sim N(0, 1)$, which is the standard assumption of probit models.
#### What Survives: The Exogeneity Test
The linear-residual construction remains **fully valid as a Hausman-type exogeneity test** [@wooldridge-2010-crosssection, p. 744]. Under the null hypothesis of exogeneity ($H_0: \rho_1 = 0$), $c_1$ does not enter the model, so none of the joint independence assumptions are needed. A standard $t$-test on $\hat{v}_2 = 0$ in the Poisson second stage provides an asymptotically valid test for endogeneity.
---
### 3.3 Terza's (1998) Exact Correction
To resolve the discreteness problem, @terza-1998 and @wooldridge-2010-crosssection (sec. 18.5) replace the linear probability model with a **latent threshold model**:
$$y_2 = \mathbf{1}[\mathbf{z}\boldsymbol{\pi}_2 + v_2 \ge 0], \qquad v_2 \mid \mathbf{z} \sim N(0, 1)$$ {#eq-poisson-iv-6}
Let $(c_1, v_2)$ be jointly normal, mean zero, and independent of $\mathbf{z}$, with $c_1 = \rho_1 v_2 + e_1$, where $e_1 \mid \mathbf{z}, v_2 \sim N(0, \tau_1^2 - \rho_1^2)$.
Conditioning on the latent error $v_2$ yields:
$$E(y_1 \mid \mathbf{z}, v_2) = \exp\bigl((\tau_1^2 - \rho_1^2)/2 + \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2\bigr)\exp(\rho_1 v_2)$$
Averaging over the truncated normal distribution of $v_2$ conditional on the observed binary outcome $y_2 \in \{0, 1\}$ gives the exact conditional mean:
$$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_1)$$ {#eq-poisson-iv-7}
where:
$$h(y_2, \mathbf{z}\boldsymbol{\pi}_2, \rho_1) = \exp(\rho_1^2 / 2) \left[ y_2 \frac{\Phi(\rho_1 + \mathbf{z}\boldsymbol{\pi}_2)}{\Phi(\mathbf{z}\boldsymbol{\pi}_2)} + (1 - y_2) \frac{1 - \Phi(\rho_1 + \mathbf{z}\boldsymbol{\pi}_2)}{1 - \Phi(\mathbf{z}\boldsymbol{\pi}_2)} \right]$$ {#eq-poisson-iv-8}
Here $\Phi(\cdot)$ denotes the standard normal CDF, $\rho_1 = \mathrm{Cov}(c_1, v_2)$ is a projection coefficient on the scale of the probit index, and $\exp(\rho_1^2 / 2)$ is the lognormal correction factor.
#### Estimation via 1D Profile Likelihood with Offset
Notice that for any fixed value of $\rho_1$, $\log h(y_2, \mathbf{z}\hat{\boldsymbol{\pi}}_2, \rho_1)$ is a known data column. The conditional mean can be written as:
$$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(\rho_1)\bigr)$$ {#eq-poisson-iv-9}
This is a standard Poisson regression with an **offset**.
Because $h(\cdot)$ enters multiplicatively, individual fixed effects factor out in closed form:
$$\exp(\hat{\alpha}_i) = \frac{\sum_t y_{1,it}}{\sum_t \exp(\mathbf{z}_{1,it}\boldsymbol{\delta}_1 + \alpha_1 y_{2,it}) h_{it}(\rho_1)}$$
Fixed effects are absorbed without an incidental parameters problem. The full model is estimated by profiling out $(\boldsymbol{\delta}_1, \alpha_1)$ and optimizing over the single scalar $\rho_1 \in [-2, 2]$ using 1D optimization.
#### Terza's Estimators and FIML
@terza-1998 proposes three estimators:
1. **Full-Information Maximum Likelihood (FIML):** Evaluates the joint likelihood of $y_1$ and $y_2$ by integrating out $v_2$ via Gauss–Hermite quadrature.
2. **Two-Stage Method of Moments / QMLE:** The two-step estimator above (probit followed by QMLE on @eq-poisson-iv-7).
3. **Nonlinear Weighted Least Squares:** For the fully parametric specification.
In panel data, FIML has two major drawbacks:
- It requires the full Poisson density inside the quadrature integral, losing QMLE robustness to overdispersion.
- It cannot absorb fixed effects; entering unit dummies destroys likelihood separability, leading to severe incidental parameters bias and convergence failures.
By contrast, the two-step profile QMLE concentrates out fixed effects cleanly.
---
### 3.4 Generalized-Residual CF and the Mundlak Device {#sec-poisson-iv-gr}
A first-order Taylor expansion of $h(y_2, \mathbf{z}\boldsymbol{\pi}_2, \rho_1)$ around $\rho_1 = 0$ yields:
$$\left. \frac{\partial}{\partial \rho_1} \frac{\Phi(\rho_1 + \mathbf{z}\boldsymbol{\pi}_2)}{\Phi(\mathbf{z}\boldsymbol{\pi}_2)} \right|_{\rho_1 = 0} = \frac{\phi(\mathbf{z}\boldsymbol{\pi}_2)}{\Phi(\mathbf{z}\boldsymbol{\pi}_2)} \equiv \lambda(\mathbf{z}\boldsymbol{\pi}_2)$$ {#eq-poisson-iv-10}
$$\left. \frac{\partial}{\partial \rho_1} \frac{1 - \Phi(\rho_1 + \mathbf{z}\boldsymbol{\pi}_2)}{1 - \Phi(\mathbf{z}\boldsymbol{\pi}_2)} \right|_{\rho_1 = 0} = -\frac{\phi(\mathbf{z}\boldsymbol{\pi}_2)}{1 - \Phi(\mathbf{z}\boldsymbol{\pi}_2)} = -\lambda(-\mathbf{z}\boldsymbol{\pi}_2)$$
Since $h(y_2, \mathbf{z}\boldsymbol{\pi}_2, 0) = 1$ and $\left.\frac{d}{d\rho_1}\exp(\rho_1^2/2)\right|_{\rho_1=0} = 0$, we have:
$$h(y_2, \mathbf{z}\boldsymbol{\pi}_2, \rho_1) \approx 1 + \rho_1 \hat{r}_2$$ {#eq-poisson-iv-11}
where $\hat{r}_2$ is the **generalized residual** [@gourieroux-monfort-renault-trognon-1987]:
$$\hat{r}_2 = y_2 \frac{\phi(\mathbf{z}\hat{\boldsymbol{\pi}}_2)}{\Phi(\mathbf{z}\hat{\boldsymbol{\pi}}_2)} - (1 - y_2)\frac{\phi(\mathbf{z}\hat{\boldsymbol{\pi}}_2)}{1 - \Phi(\mathbf{z}\hat{\boldsymbol{\pi}}_2)}$$
Using $1 + u \approx \exp(u)$ for small $u = \rho_1 \hat{r}_2$, the conditional mean becomes:
$$E(y_1 \mid \mathbf{z}, y_2) \approx \exp(\mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 \hat{r}_2)$$
The generalized-residual CF is a **first-order local approximation to Terza's exact correction with error $O(\rho_1^2)$**.
#### Probit vs. Logit Generalized Residuals
For a **logit** first stage ($\Lambda(\eta) = \frac{1}{1 + e^{-\eta}}$), the generalized residual collapses algebraically:
$$\hat{r}_2 = y_2 \frac{\Lambda'(\eta)}{\Lambda(\eta)} - (1 - y_2)\frac{\Lambda'(\eta)}{1 - \Lambda(\eta)} = y_2 (1 - \hat{p}_2) - (1 - y_2)\hat{p}_2 = y_2 - \hat{p}_2$$ {#eq-poisson-iv-12}
The logit generalized residual is simply the **raw response residual** ($y_2 - \hat{p}_2$).
#### Dealing with Panel Fixed Effects: The Mundlak Device
In short panels ($T$ small, $N \to \infty$), estimating individual unit dummies in a nonlinear first stage (probit or logit) suffers from the **incidental parameters problem**, producing heavily biased first-stage estimates and distorted generalized residuals.
As @lin-wooldridge-2019 show, the standard solution is the **Mundlak device** (Correlated Random Effects):
1. Add the unit-specific time averages of all time-varying exogenous covariates and instruments:
$$\mathbf{w}_{it} = (1, \mathbf{z}_{it}, \bar{\mathbf{z}}_i)$$
2. Estimate a pooled probit or logit model on $\mathbf{w}_{it}$ and compute the generalized residuals $\hat{r}_{2,it}$.
3. Include $\hat{r}_{2,it}$ in a fixed-effects Poisson second stage. The unit dummies are absorbed there, and FE Poisson is consistent with $T$ fixed, so the incidental parameters problem does not reappear.
---
### 3.5 Multiplicative-Moment GMM
Instead of a control function, @mullahy-1997 proposes multiplicative moment conditions:
$$E\left[ \mathbf{z}' \left( y_1 \exp(-\mathbf{z}_1\boldsymbol{\delta}_1 - \alpha_1 y_2) - 1 \right) \right] = \mathbf{0}$$ {#eq-poisson-iv-13}
Under $E[\exp(c_1) \mid \mathbf{z}] = 1$, this condition holds without requiring the first-stage model for $y_2$ to be correctly specified.
#### Panel Fixed Effects via Quasi-Differencing
In panel data with unobserved individual heterogeneity $\eta_i$, $E[y_{it} \mid \mathbf{z}_i, \eta_i] = \eta_i \exp(\mathbf{x}_{it}\boldsymbol{\beta})$. Individual effects cannot be removed by additive demeaning; they must be eliminated via **multiplicative quasi-differencing** [@chamberlain-1992; @wooldridge-1999-multiplicative; @windmeijer-2000]:
$$E\left[ \mathbf{z}_{it}' \left( y_{it} \frac{\exp(\mathbf{x}_{is}\boldsymbol{\beta})}{\exp(\mathbf{x}_{it}\boldsymbol{\beta})} - y_{is} \right) \right] = \mathbf{0}$$ {#eq-poisson-iv-14}
This requires repeated time periods and strict exogeneity of instruments across periods.
---
## 4. Why Recommend Generalized Residuals Instead of Exact Terza?
Given that Terza's correction [@terza-1998] is theoretically exact under bivariate normality, why do applied econometricians and textbooks [@wooldridge-2010-crosssection; @wooldridge-2014-qmle; @lin-wooldridge-2019] frequently recommend the **Generalized-Residual Control Function**?
There are five key reasons:
### 1. Built-in Hausman Exogeneity Test in a Single Step
In the generalized residual model:
$$\log E(y_1 \mid \mathbf{z}, y_2) = \mathbf{z}_1\boldsymbol{\delta}_1 + \alpha_1 y_2 + \rho_1 \hat{r}_2$$
The coefficient on $\hat{r}_2$ estimates $\rho_1$, and its $t$-statistic is a **direct, asymptotically valid variable-addition test for exogeneity** ($H_0: \rho_1 = 0$).
- Under $H_0$, none of the joint normality or reduced-form assumptions are required [@wooldridge-2010-crosssection, p. 744].
- Terza's exact method does not yield a simple linear $t$-statistic directly beside $\hat{\alpha}_1$.
### 2. Standard Software Drop-in
- **Generalized Residual:** Fits directly in standard software without custom numerical optimization.
- **Terza’s Exact 2-Step:** Requires writing a custom 1D profile likelihood wrapper to search over $\rho_1$ with an offset.
### 3. Distributional Flexibility: Beyond the Probit Link
- **Terza’s $h(\cdot)$ is strictly probit-specific:** The formula in @eq-poisson-iv-8 is derived exclusively from truncated moments of a bivariate normal distribution, and no closed-form analogue for logit appears in @terza-1998, @wooldridge-2010-crosssection or @imbens-wooldridge-2007-cf.
- **Generalized residuals apply to any GLM link:** For logit, the generalized residual collapses to the simple **raw response residual** $\hat{r}_2 = y_2 - \hat{p}_2$, requiring zero Mills-ratio machinery.
### 4. Robustness Against Non-Normality (Curvature Misspecification)
- Terza’s correction is **exact only under strict joint bivariate normality**.
- If the true unobservables $(c_1, v_2)$ are non-normal (e.g. skewed or heavy-tailed), $h(\cdot)$ has the wrong functional form. The profile then picks the $\rho_1$ that best fits a curve the data do not follow, and that misspecification passes straight into $\hat{\alpha}_1$. This is misspecification bias, not overfitting — there is only one extra parameter, so there is nothing to overfit with.
- The generalized residual commits to less. It controls for $E[v_2 \mid y_2, \mathbf{z}]$ and stops at first order, so a wrong $h$ cannot drag it. Whether this actually wins under non-normality is an open question; @lin-wooldridge-2019 leave "investigations into the quality of CF approximations in discrete cases to future research.
### 5. The Approximation Error Shrinks with the Endogeneity
The error is $O(\rho_1^2)$, so it vanishes quadratically as $\rho_1 \to 0$ — and $\rho_1$ is not a quantity you have to guess at, because it is estimated as the coefficient on $\hat{r}_2$. The error is also asymmetric: $h$ scales the treated mean by a factor bounded by $1/\Phi(\mathbf{z}\boldsymbol{\pi}_2)$, which is unbounded in the tail for $\rho_1 > 0$ but shrinks smoothly toward zero for $\rho_1 < 0$, so positive confounding is the case to worry about. We give no numerical threshold here: the size of the error at a given $\rho_1$ depends on how far the probit index reaches into the tail, so it is design-specific rather than a property of $\rho_1$ alone.
### Practical Strategy
| Situation | Recommended Estimator | Rationale |
|:---|:---|:---|
| **Routine empirical analysis & testing** | **Generalized-Residual CF** | Fast, standard 1-line command, built-in $t$-test for $\rho_1 = 0$, works with Logit ($y_2 - \hat{p}_2$), handles panel FE via Mundlak. |
| **$\hat{\rho}_1$ small, or negative** | **Generalized-Residual CF alone** | The $O(\rho_1^2)$ error is smallest exactly here, and it is the negative side of the asymmetry. |
| **$\hat{\rho}_1$ large and positive** | **Add the Terza Exact Profiled 2-Step** | This is where the Taylor error bites hardest. Fit the exact model alongside and report both; if they agree the approximation was safe, and if they do not, the exact one is right provided the probit is. |
---
## 5. Overview and Comparison of the Four Estimators
| Estimator | First Stage | Second Stage Adjustment | Consistency (Binary $y_2$) | Panel Fixed Effects | Standard Errors |
|:---|:---:|:---:|:---:|:---:|:---:|
| **Linear CF** | Linear OLS | Add $\hat{v}_2$ inside $\exp(\cdot)$ | Inconsistent | Yes (Poisson with FE) | Cluster Bootstrap |
| **Generalized CF** | Probit / Logit (with Mundlak) | Add $\hat{r}_2$ inside $\exp(\cdot)$ | Approx. ($O(\rho_1^2)$ error) | Yes (Mundlak in 1st stage) | Cluster Bootstrap |
| **Terza Exact 2-Step** | Probit | Profile $\rho_1$ with offset $\log h(\rho_1)$ | **Yes** (under normality) | Yes (Poisson with offset) | Cluster Bootstrap |
| **Multiplicative GMM** | Any ($\hat{p}_2$) | Instrument $y_2$ in moments | **Yes** (distribution-free) | Needs quasi-differencing | Sandwich vcov |
---
## 6. Practical Guidance: An Observable Diagnostic Workflow
In empirical practice, we do not need to guess the confounding parameter $\rho_1$—**it is estimated directly as the coefficient on the generalized residual $\hat{r}_2$ in your second-stage Poisson regression table.**
A standard applied workflow proceeds in four steps:
1. **Estimate the Generalized-Residual CF as the baseline model:**
- For continuous $y_2$: Add the OLS residual $\hat{v}_2$.
- For binary $y_2$: Add the generalized residual $\hat{r}_2$ (or raw response residual $y_2 - \hat{p}_2$ if using logit).
- In short panels ($T$ small): Use the Mundlak device ($\mathbf{w}_{it} = (1, \mathbf{z}_{it}, \bar{\mathbf{z}}_i)$) in the first stage to avoid incidental parameters bias.
2. **Check the $t$-statistic on the residual (Exogeneity Test):**
The $t$-statistic on $\hat{r}_2$ (testing $\hat{\rho}_1 = 0$) is an asymptotically valid Hausman-type exogeneity test. Read it as a description of how much endogeneity the data can detect. Do not read it as a licence to drop the residual — see the pretest warning below.
3. **Read the size and sign of $\hat{\rho}_1$ (Approximation Diagnostic):**
- **If $\hat{\rho}_1$ is small or negative:** the $O(\rho_1^2)$ error is at its smallest, and on the benign side of the asymmetry. The Generalized-Residual CF stands on its own.
- **If $\hat{\rho}_1$ is large and positive:** the Taylor error is upward and grows quadratically. Fit Terza's exact 2-step alongside it and report both. Do not read a threshold off $\hat{\rho}_1$ alone — how much the error matters also depends on how far your probit index reaches into the tail, so the comparison is the diagnostic, not a rule of thumb.
4. **Inference:**
Always compute standard errors using a cluster bootstrap over panel units to account for first-stage estimation uncertainty.
### Do Not Use the Exogeneity Test to Choose the Model
The tempting next move is to let step 2 pick the specification: if the $t$-statistic on $\hat{r}_2$ is insignificant, drop the residual and report the plain FE Poisson; if it is significant, keep the control function. This is a pretest estimator, and it is worse than either of the two models it chooses between.
@guggenberger-2010 works this out for the linear case. A Hausman pretest followed by OLS-or-2SLS, with the second-stage $t$-test reported at its nominal level, has **asymptotic size equal to one**: there are points in the parameter space where the procedure rejects a true null with probability approaching 1. The mechanism is a sequence of local alternatives. Let the endogeneity drift toward zero at rate $n^{-1/2}$. The pretest then has low power and often fails to reject, so the procedure falls back on the estimator that ignores endogeneity — but at that drift rate the bias of that estimator is the same order as its standard error, so the reported confidence interval is centred in the wrong place by about its own width. Taking the supremum of rejection probability over the parameter space picks out exactly these points. @guggenberger-2010-panel obtains the same conclusion for the panel Hausman test between fixed and random effects.
These results are proved for linear models, and we are not aware of an analogous theorem for the Poisson control function. The mechanism does not depend on linearity: it needs only a pretest with low power against local endogeneity and a fallback estimator whose bias does not vanish faster than its standard error, and both are present here. We treat it as a warning that carries over rather than a result that has been established in this setting.
The practical rule is to **report the naive FE Poisson and the control-function FE Poisson side by side** and let the reader see both. The $t$-statistic on $\hat{r}_2$ still belongs in the table; it is evidence about the size of $\rho_1$, which is what step 3 needs. What it must not do is silently delete one of the two columns.
---
## 7. Implementation Demos in R and Stata
Below are concise, copy-pasteable implementation templates in R (`fixest`) and Stata (`ppmlhdfe`).
### 7.1 R Implementation Demo (`fixest`)
```{r demo-r}
#| echo: true
library(fixest)
# 1. Simulate a short panel dataset (50 units, 5 periods)
set.seed(123)
n_i <- 50; n_t <- 5; n <- n_i * n_t
id <- rep(1:n_i, each = n_t)
time <- rep(1:n_t, times = n_i)
z <- rnorm(n) # Excluded instrument
x <- rnorm(n) # Exogenous covariate
c_i <- rnorm(n_i)[id]
e2 <- rnorm(n)
# Endogenous binary treatment and Poisson outcome
y2 <- as.integer(0.8 * z + 0.5 * x + c_i + e2 >= 0)
lambda <- exp(0.5 + 0.8 * y2 + 0.4 * x + c_i + 0.2 * e2)
y1 <- rpois(n, lambda)
df <- data.frame(y1 = y1, y2 = y2, z = z, x = x, id = id, time = time)
# 2. Mundlak Device: calculate unit means
df$mean_z <- ave(df$z, df$id)
df$mean_x <- ave(df$x, df$id)
# 3. Step 1: Probit First Stage & Generalized Residual
fs <- feglm(y2 ~ z + x + mean_z + mean_x, data = df, family = binomial(probit))
eta <- predict(fs, type = "link")
df$r2 <- ifelse(df$y2 == 1, dnorm(eta) / pnorm(eta), -dnorm(eta) / (1 - pnorm(eta)))
# (Note: If using Logit first stage, r2 is simply y2 - predict(fs, type = "response"))
# 4. Step 2: Second Stage Poisson with High-Dimensional Fixed Effects
m_naive <- fepois(y1 ~ y2 + x | id + time, data = df, notes = FALSE)
m_cf <- fepois(y1 ~ y2 + r2 + x | id + time, data = df, notes = FALSE)
etable(m_naive, m_cf, headers = c("Naive FE Poisson", "Generalized Residual CF"))
cat(sprintf("\ntrue alpha1 = %.2f, true rho1 = %.2f\n", 0.8, 0.2))
cat(sprintf("naive alpha1 = %.4f\n", coef(m_naive)["y2"]))
cat(sprintf("generalized CF alpha1 = %.4f (SE %.4f), rho1 = %.4f\n",
coef(m_cf)["y2"], se(m_cf)["y2"], coef(m_cf)["r2"]))
```
With 50 units over 5 periods the standard error is large enough that this single
draw cannot separate the two estimators; it is a template for the syntax, not
evidence about the estimator. The naive coefficient is reported beside it so the
direction of the correction is visible.
### 7.2 Stata Implementation Demo (`ppmlhdfe`)
```stata
* ==============================================================================
* Generalized-Residual Control Function in Stata using ppmlhdfe
* ==============================================================================
* 1. Mundlak Device: Calculate unit-level means of time-varying variables
bysort id: egen mean_z = mean(z)
bysort id: egen mean_x = mean(x)
* 2. Step 1: Probit First Stage on pooled data with Mundlak means
probit y2 z x mean_z mean_x
predict double eta, xb
gen double r2 = cond(y2 == 1, normalden(eta)/normal(eta), -normalden(eta)/(1 - normal(eta)))
* (Alternatively, for a Logit first stage: raw response residual)
* logit y2 z x mean_z mean_x
* predict double phat, pr
* gen double r2 = y2 - phat
* 3. Step 2: Fixed-Effects Poisson Second Stage
ppmlhdfe y1 y2 r2 x, absorb(id time)
* 4. Inference: Panel Cluster Bootstrap Program
capture program drop boot_poisson_cf
program define boot_poisson_cf, rclass
* Mundlak means must be rebuilt on the RESAMPLED clusters (new_id), not
* carried in from the original data, or every replication reuses the
* means of the original panel and the bootstrap understates the variance.
quietly bysort new_id: egen double mz_b = mean(z)
quietly bysort new_id: egen double mx_b = mean(x)
* First stage
quietly probit y2 z x mz_b mx_b
quietly predict double eta_b, xb
quietly gen double r2_b = cond(y2 == 1, normalden(eta_b)/normal(eta_b), -normalden(eta_b)/(1 - normal(eta_b)))
* Second stage
quietly ppmlhdfe y1 y2 r2_b x, absorb(new_id time)
return scalar alpha1 = _b[y2]
return scalar rho1 = _b[r2_b]
drop eta_b r2_b mz_b mx_b
end
bootstrap r(alpha1) r(rho1), reps(500) cluster(id) idcluster(new_id): boot_poisson_cf
```