Why Stack the Equations?
A two-step estimator first estimates a parameter \(\eta\), then uses that estimate to estimate a target parameter \(\beta\). The second stage is therefore affected by two sources of sampling variation: its own estimating equation and the estimated first-stage quantity it uses. If the two stages use the same observations, their sampling errors can also be correlated. Stacking both equations makes these contributions visible in one joint asymptotic distribution.
The point is not that two sequential calculations become independent when written together. Rather, stacking expresses their dependence explicitly. The resulting covariance matrix supplies standard errors for either stage and for functions of both stages.
Here \(\eta\) and \(\beta\) may be vectors of different dimensions, and the functions \(g\) and \(h\) are vector-valued with matching dimensions. Write their population values as \(\eta_0\) and \(\beta_0\), where \(E[g(W_i,\eta_0)]=0\) and \(E[h(W_i,\beta_0,\eta_0)]=0\). The first equation identifies the first-stage parameter; the second identifies the target when the first-stage argument is at its population value.
For the basic result below, assume the observations are identically distributed and independent, both population equations have unique local roots, and the estimators are consistent for those roots. Assume the estimating functions are differentiable near the population values, their sample Jacobians converge to their population Jacobians, and the Taylor remainders are negligible at the root-\(n\) scale. Also assume finite second moments and a nonsingular block Jacobian. These are regularity conditions that allow a first-order expansion and a multivariate central limit theorem; they do not themselves guarantee identification or consistency.
The Stacked Linearization
Let \(A\) be the derivative of the first equation with respect to \(\eta\), \(B\) the derivative of the second equation with respect to \(\beta\), and \(C\) the derivative of the second equation with respect to \(\eta\), all evaluated at the population values:
The prime denotes transpose, so each derivative is a matrix with one row per equation and one column per parameter. Since the first equation does not depend on \(\beta\), the population Jacobian of the stacked system is block lower triangular:
Define the stacked observation-level score \(m_i=(g(W_i,\eta_0)',h(W_i,\beta_0,\eta_0)')'\), and let \(\Omega=E[m_i m_i']\). Because the population means of both score blocks are zero, \(\Omega\) is their covariance matrix. In particular, its off-diagonal blocks contain the covariance between the first- and second-stage scores.
A first-order Taylor expansion around \((\eta_0,\beta_0)\) gives the stacked score average plus \(J\) times the parameter error, with a remainder negligible after multiplication by \(\sqrt{n}\).
The scaled average \(n^{-1/2}\sum_{i=1}^n m_i\) converges in distribution to a mean-zero normal vector with covariance \(\Omega\).
Multiplication by the inverse Jacobian gives the joint limit, with the lower block reflecting how first-stage error enters the second equation.
More explicitly, the first row of the linearized system gives the first-stage expansion. Substituting it into the second row gives the second-stage expansion:
where \(g_i=g(W_i,\eta_0)\) and \(h_i=h(W_i,\beta_0,\eta_0)\). The term involving \(C A^{-1}g_i\) is the first-stage contribution to the second-stage influence function. If \(C\) is zero, a small first-stage error has no first-order effect on the second equation. Otherwise it generally matters.
Why the Covariance Formula Has This Form
The expansion of the full stacked equation is \(0=n^{-1}\sum_i m_i+J((\hat{\eta}-\eta_0)',(\hat{\beta}-\beta_0)')'+o_p(n^{-1/2})\). Multiplying by \(\sqrt{n}\) and solving yields minus \(J^{-1}\) times the scaled score average, plus a term converging to zero in probability. The central limit theorem gives the score average covariance \(\Omega\). A linear transformation by \(-J^{-1}\) therefore gives covariance \((-J^{-1})\Omega(-J^{-1})'=J^{-1}\Omega J^{-T}\); the two negative signs cancel.
For the two-stage block structure, direct inversion gives
This inverse produces the influence functions shown above. The full \(\Omega\) is essential: omitting its cross-stage covariance generally changes both the second-stage variance and the covariance between estimates. In practice, the unknown Jacobian and score covariance are replaced by consistent sample estimates. The resulting sandwich estimate approximates the covariance of the scaled estimator; divide it by \(n\) to estimate the covariance of the estimator itself.
Worked Example: A Scalar Stacked Covariance
Suppose a hypothetical two-equation estimator has scalar \(\eta\) and \(\beta\). At the population values, let \(A=2\), \(B=3\), and \(C=1\). Suppose the covariance matrix of the two scores is \(\Omega=\begin{pmatrix}4&1\\1&9\end{pmatrix}\). These are assumed population quantities for illustration, not estimates from a real study. Compute the asymptotic covariance of \(\sqrt{n}((\hat{\eta}-\eta_0),(\hat{\beta}-\beta_0))'\).
The first-stage influence function is \(-g_i/2\). The second-stage influence function is \(-\frac{1}{3}(h_i-\frac{1}{2}g_i)\), since \(C A^{-1}=1/2\). Thus
The last \(-1\) inside the parentheses is the cross term \(-2(1/2)\operatorname{Cov}(g_i,h_i)=-1\). The covariance between the two influence functions is \(\operatorname{Cov}(-g_i/2,-h_i/3+g_i/6)=1/6-1/3=-1/6\). Therefore the joint asymptotic covariance is \(\begin{pmatrix}1&-1/6\\-1/6&1\end{pmatrix}\). The asymptotic covariance of the estimator vector itself is approximately this matrix divided by \(n\).
Worked Calculations from Hypothetical Observations
Worked Example: Plugging In the First-Stage Estimate
Consider the illustrative equations \(g_i=z_i-\eta\) and \(h_i=y_i-\beta-\eta z_i\), using the same four hypothetical paired observations: \((z_i,y_i)=(0,1),(1,2),(2,2),(3,5)\). Assume these represent independent draws from a population for the purpose of the asymptotic calculation. The first-stage equation sets \(\hat{\eta}=\bar{z}\), and the second sets \(\hat{\beta}=\bar{y}-\hat{\eta}\bar{z}\).
The sample means are \(\bar{z}=(0+1+2+3)/4=1.5\) and \(\bar{y}=(1+2+2+5)/4=2.5\). Hence \(\hat{\eta}=1.5\) and \(\hat{\beta}=2.5-(1.5)(1.5)=0.25\). For these equations, the sample Jacobian blocks at the estimates are \(A=-1\), \(B=-1\), and \(C=-\bar{z}=-1.5\). The observation-level estimated influence values are \(\widehat{\mathrm{IF}}_{\eta,i}=z_i-\hat{\eta}\) and \(\widehat{\mathrm{IF}}_{\beta,i}=h_i-C A^{-1}g_i\), where \(g_i=z_i-\hat{\eta}\) and \(h_i=y_i-\hat{\beta}-\hat{\eta}z_i\).
The values are shown below; the influence-function formula is the one derived from the stacked equations, with the estimated Jacobians substituted.
| Observation | \(g_i\) | \(h_i\) | \(\widehat{\mathrm{IF}}_{\eta,i}\) | \(\widehat{\mathrm{IF}}_{\beta,i}\) |
|---|---|---|---|---|
| 1 | -1.5 | 0.75 | -1.5 | 3 |
| 2 | -0.5 | 0.25 | -0.5 | 1 |
| 3 | 0.5 | -1.25 | 0.5 | -2 |
| 4 | 1.5 | 0.25 | 1.5 | -2 |
Here \(C A^{-1}=1.5\), so the second-stage influence values are \(h_i-1.5g_i\). For example, observation 1 gives \(0.75-1.5(-1.5)=3\). The empirical covariance of the influence vectors, using divisor \(n=4\), is
For instance, the off-diagonal entry is \([(-1.5)(3)+(-0.5)(1)+(0.5)(-2)+(1.5)(-2)]/4=-2.25\). In the usual large-sample approximation, the estimated covariance of \((\hat{\eta},\hat{\beta})'\) is this matrix divided by \(n\). With only four observations the approximation may be poor; the calculation illustrates the stacked covariance mechanics, not a claim of reliable small-sample inference.
Worked Example: A Two-Component Second-Stage Target
Suppose a hypothetical system has scalar \(\eta\) and two-component \(\beta\). Let \(A=2\), \(B=\begin{pmatrix}3&0\\0&2\end{pmatrix}\), and \(C=(1,2)'\). Assume the score covariance, ordered as \((g,h_1,h_2)'\), is \(\Omega=\begin{pmatrix}4&1&-2\\1&9&2\\-2&2&16\end{pmatrix}\). Find the joint covariance of the scaled estimator and the variance for the sum \(\beta_1+\beta_2\).
The influence functions are \(\mathrm{IF}_{\eta}=-g/2\), \(\mathrm{IF}_{\beta_1}=-(h_1-g/2)/3=-h_1/3+g/6\), and \(\mathrm{IF}_{\beta_2}=-(h_2-g)/2=-h_2/2+g/2\). Using the stated score variances and covariances gives \(\operatorname{Var}(\mathrm{IF}_{\eta})=4/4=1\), \(\operatorname{Var}(\mathrm{IF}_{\beta_1})=(9+1-1)/9=1\), and \(\operatorname{Var}(\mathrm{IF}_{\beta_2})=(16+4+4)/4=6\).
The cross-covariances are \(\operatorname{Cov}(\mathrm{IF}_{\eta},\mathrm{IF}_{\beta_1})=1/6-4/12=-1/6\) and \(\operatorname{Cov}(\mathrm{IF}_{\eta},\mathrm{IF}_{\beta_2})=-2/4-4/4=-1.5\). For the two second-stage components, \(\operatorname{Cov}(\mathrm{IF}_{\beta_1},\mathrm{IF}_{\beta_2})=(2-1+1+2)/6=2/3\). Thus the joint covariance, in the order \((\eta,\beta_1,\beta_2)\), is \(\begin{pmatrix}1&-1/6&-1.5\\-1/6&1&2/3\\-1.5&2/3&6\end{pmatrix}\). The variance of the scaled sum \(\beta_1+\beta_2\) is \(1+6+2(2/3)=25/3\), or approximately \(8.333\). Its estimator has asymptotic variance approximately \((25/3)/n\). The positive covariance between the second-stage components increases the variance of their sum.
Common Pitfall and Practical Interpretation
A frequent mistake is to calculate the second-stage variance as if \(\hat{\eta}\) were known, using only the second-stage score \(h_i\). This ignores the \(C A^{-1}g_i\) term and can omit both the first-stage variance contribution and covariance between \(g_i\) and \(h_i\). Even if the two score blocks are uncorrelated, first-stage uncertainty can still affect the second stage when \(C\ne0\). Conversely, when the scores are correlated, dropping their covariance can change the answer in either direction.
The stacked approach also clarifies what is being approximated. The covariance \(J^{-1}\Omega J^{-T}\) describes the limiting distribution of the scaled joint estimation error. Its plug-in estimate, divided by \(n\), is used for the covariance of the unscaled estimator. A standard error for one component is the square root of its estimated variance; a standard error for a linear combination uses the full covariance matrix, not just the separate marginal variances.
Check Your Understanding
Use the stacked-equation framework to answer the following questions.
- Why is the upper-right block of the stacked Jacobian zero for the two-step system defined here?
- In the second-stage influence function, what do \(C\) and \(A^{-1}g_i\) represent?
- For scalar \(A=2\), \(B=3\), and \(C=1\), write the second-stage influence function in terms of \(g_i\) and \(h_i\).
- Why can setting the covariance between the two score blocks to zero still leave a first-stage contribution to second-stage variance?
- If a target is \(\beta_1+\beta_2\), which entries of the joint covariance matrix are needed to compute its asymptotic variance?