Estimating equations and asymptotic linearity

Why estimators solve equations, how they fluctuate, and where sandwich standard errors come from.

Before reading: Basic calculus, expectations, and Large-sample reasoning for methods research. Its sample-mean identity motivates the general equation here. The later residualization example uses Conditional expectation and L2 geometry and can be revisited after that page.

1 The one idea

An estimator is often introduced as a formula: a mean, a regression coefficient, a maximum likelihood estimate. For theory, a more useful description is the equation it solves.

Let W_1,\ldots,W_n be i.i.d. observation records: one record might contain a person’s outcome and covariates. Here W denotes the whole record, while Y denotes an outcome within it. Let \theta be a candidate parameter value and \psi(W;\theta) its estimating function, with the same dimension as \theta. Here \psi names an equation contribution, rather than the population mean. Write P_nf=n^{-1}\sum_i f(W_i) for a sample average and Pf=E[f(W)] for its population expectation.

An estimating equation chooses a value \hat\theta that balances the average contribution at zero:

\begin{aligned} P_n\psi(W;\hat\theta)&= \frac{1}{n}\sum_{i=1}^n \psi(W_i;\hat\theta)\\ &=0. \end{aligned}

The population version is

\begin{aligned} P\psi(W;\theta_0)&=E[\psi(W;\theta_0)]\\ &=0. \end{aligned}

The population equation defines the value \theta_0 that this procedure aims to estimate. To connect the two equations, we need consistency: the selected sample root must approach that population root in probability. An LLN is part of this argument, but convergence of the equation at each fixed parameter value is not automatically convergence of a root selected from the data. The root must be identified and the equation controlled near it.

Once consistency is established, a local Taylor expansion explains how random equation noise moves the estimated root. This gives the recurring sequence: identify the population root, show the sample root approaches it, linearize, and justify uncertainty. Means, ordinary least squares, likelihood estimators, and many other procedures fit this pattern.

2 Familiar estimators as equations

Use \mu=E[Y] for the population mean. If a is a candidate mean, its sample equation is

\begin{aligned} P_n(Y-a)&=\bar Y-a\\ &=0. \end{aligned}

The solution is \hat\mu=\bar Y. The population equation is E[Y-a]=\mu-a=0, with unique root a=\mu. The LLN gives \bar Y\xrightarrow{p}\mu, directly proving consistency in this example. Positive and negative deviations balance at the mean because they sum to zero there.

The least-squares coefficient in a regression of Y on a covariate vector X solves the normal equations below. Include a constant component in X when fitting an intercept.

P_n\left[X\{Y-X^\top\beta\}\right]=0.

With suitable finite moments, the population equation is E[XY]-E[XX^\top]\beta=0. If E[XX^\top] is nonsingular, its solution is the projection coefficient \beta^*=E[XX^\top]^{-1}E[XY]. This remains a population functional when the conditional mean is nonlinear. An LLN for the sample moments and continuous mapping through the matrix inverse establish consistency for that coefficient. Whether it answers the scientific question is a separate issue, developed in What is your estimator actually estimating?.

An estimator defined by a minimization problem is often an estimating-equation estimator too. If

\hat\theta = \mathop{\mathrm{argmin}}_\theta P_n\ell(W;\theta),

and \ell is differentiable, an interior, unconstrained minimizer satisfies

P_n\frac{\partial}{\partial\theta}\ell(W;\hat\theta)=0.

These optimization-based procedures are called M-estimators. Estimators defined directly by a zero of P_n\psi are often called Z-estimators. A zero derivative is a necessary condition for the stated interior minimum; it need not distinguish a minimum from another stationary point. Boundary constraints and nonsmooth losses require the corresponding optimization conditions.

Figure 1: For these five observations, P_n(Y-\mu)=3-\mu. Its zero is also the minimum of P_n(Y-\mu)^2=(\mu-3)^2+2. The derivative of that loss is -2P_n(Y-\mu): minimization and zero-solving describe the same estimator.

3 From the sample equation to a linear approximation

First consider a scalar parameter and a consistent sample root. At the population root, P\psi(W;\theta_0)=0, but the sample average P_n\psi(W;\theta_0) usually is not zero. This discrepancy is equation noise: a new sample produces a different discrepancy even though \theta_0 stays fixed.

Changing the candidate parameter changes the average contribution. A first-order expansion gives

\begin{aligned} 0&=P_n\psi(W;\hat\theta)\\ &\approx P_n\psi(W;\theta_0)+A(\hat\theta-\theta_0), \end{aligned}

where

A = E\left[\frac{\partial}{\partial\theta}\psi(W;\theta_0)\right].

Rearranging gives

\hat\theta-\theta_0 \approx -A^{-1}P_n\psi(W;\theta_0).

The multiplier -A^{-1} translates a vertical equation error into a horizontal parameter error. A shallow population equation has small |A|, so the same equation noise can move its root farther. The minus sign moves the candidate in the direction that offsets the noise.

To make this an asymptotic linear expansion, the approximation must have a controlled remainder:

\begin{aligned} \hat\theta-\theta_0&=P_nD_0+R_n,\\ D_0(W)&=-A^{-1}\psi(W;\theta_0),\\ R_n&=o_p(n^{-1/2}). \end{aligned}

The leading term is an average of fixed, mean-zero observation contributions, since PD_0=-A^{-1}P\psi(\theta_0)=0. The remainder becomes negligible after root-n scaling. Under the conditions proved below, the estimator’s influence function is

\operatorname{IF}(W)=-A^{-1}\psi(W;\theta_0).

The Influence functions page develops influence functions for arbitrary smooth functionals. Here we see one of the main ways they arise from an estimator.

For the sample mean, \psi(Y;a)=Y-a and A=-1. Therefore

\begin{aligned} \hat\mu-\mu&=P_n(Y-\mu),\\ \operatorname{IF}(Y)&=Y-\mu. \end{aligned}

The remainder is exactly zero because the sample equation is already linear in the candidate mean. Nonlinear equations need the additional remainder argument.

Figure 2: The population equation is 4-\theta^2 and has root \theta_0=2. The sample equation is 4.6-\theta^2 and has exact root \sqrt{4.6}\approx2.145. Its tangent at 2 is 0.6-4(\theta-2), giving the approximate root 2.15. The magnified panel makes the small Taylor remainder visible. More generally, the vertical fluctuation P_n\psi(\theta_0) becomes a horizontal shift through multiplication by -A^{-1}; applying this multiplier to each observation’s estimating function gives its influence contribution.

4 Standard errors and the sandwich formula

Assume i.i.d. observations and finite, positive equation variance B. Because the population contribution has mean zero, its second moment equals its variance. The CLT gives

\begin{aligned} \sqrt n\,P_n\psi(W;\theta_0)&\xrightarrow{d}N(0,B),\\ B&=E[\psi(W;\theta_0)^2]. \end{aligned}

Combining this with the linear expansion gives

\sqrt n(\hat\theta-\theta_0) \mathrel{\overset{d}{\longrightarrow}} N\left(0,\frac{B}{A^2}\right).

Here B/A^2 is the limiting variance of the scaled error \sqrt n(\hat\theta-\theta_0). The estimator’s first-order variance is B/(nA^2), so its SE is approximately \sqrt{B/n}/|A|. This division by n distinguishes the limiting contribution variance from the uncertainty of an estimate based on n people.

For a vector parameter, the first-order estimator variance becomes

\frac{1}{n}A^{-1}B(A^{-1})^\top,

where A=E[\partial\psi/\partial\theta^\top] is the Jacobian and B=E[\psi\psi^\top] is the covariance matrix of contributions at \theta_0. The transpose on the right is needed because, for a fixed matrix L, covariance transforms as \operatorname{Var}(LX)=L\operatorname{Var}(X)L^\top. This is the sandwich variance: the equation covariance lies between two derivative-inverse factors.

For a scalar equation, compute

\begin{aligned} \widehat A&=P_n\dot\psi(W;\hat\theta),\\ \widehat B&=P_n\{\psi(W;\hat\theta)^2\},\\ \widehat{\operatorname{SE}}(\hat\theta) &=\frac{\sqrt{\widehat B}}{\sqrt n\,|\widehat A|}. \end{aligned}

The dot denotes a derivative with respect to the parameter. Consistency of \widehat A and \widehat B and A\ne0 let continuous mapping and Slutsky justify these substitutions. Because the evaluation point \hat\theta is random, the required LLN/continuity controls must hold locally, not only at a fixed parameter. For the mean, \widehat A=-1 and \widehat B=n^{-1}\sum_i(Y_i-\bar Y)^2; the resulting SE is asymptotically equivalent to s_Y/\sqrt n, with the familiar finite-sample divisor difference.

The sandwich is the general variance calculation for a justified smooth equation. Its use still depends on the sampling unit, moments, root consistency, and remainder conditions.

Figure 3: For a scalar equation, \operatorname{SE}(\hat\theta)\approx\sqrt{B/n}/|A|. Relative to the blue case, either quadrupling B or halving |A| doubles the standard error and gives the orange curve. The areas under both densities are one: the lower peak reflects a wider distribution, not less total probability. This comparison changes A or B separately; merely rescaling \psi changes both and leaves B/A^2 unchanged.

5 Example: a regression slope under heteroscedasticity

Suppose Y=\alpha+\beta X+\varepsilon with E[\varepsilon\mid X]=0, while the variance of \varepsilon depends on X. This is heteroscedasticity. Under suitable moment and rank conditions, ordinary least squares is consistent for the intercept and slope, but the usual constant-variance SE can misdescribe its uncertainty.

For the covariate vector C=(1,X)^\top, the estimating function is C(Y-C^\top\gamma), where \gamma=(\alpha,\beta)^\top. Thus A=-E[CC^\top] and B=E[CC^\top\varepsilon^2]. The dependence of squared errors on X belongs inside this expectation. The code fits \alpha=1 and \beta=0.8, with error SD 0.3+|X|.

set.seed(1)
n <- 250
x <- rnorm(n)
eps <- rnorm(n, sd = 0.3 + abs(x))
y <- 1 + 0.8 * x + eps
fit <- lm(y ~ x)

# Manual HC0 sandwich standard error for the slope.
X <- model.matrix(fit)
u <- residuals(fit)
bread <- solve(crossprod(X))
meat <- t(X) %*% diag(u^2) %*% X
vcov_hc0 <- bread %*% meat %*% bread

c(estimate = coef(fit)["x"],
  usual_SE = summary(fit)$coef["x", "Std. Error"],
  sandwich_SE = sqrt(vcov_hc0["x", "x"]))

 estimate.x    usual_SE sandwich_SE 
 0.90177237  0.08591487  0.13819053 

Both SE calculations aim to describe uncertainty in the same fitted slope. They use different assumptions about B: the usual calculation replaces the individual squared residuals with one common variance, while the sandwich retains their joint variation with the covariates. In the code, X is the design matrix with rows C_i^\top, bread is (X^\top X)^{-1}, and meat is \sum_i C_iC_i^\top\hat\varepsilon_i^2. These use sums rather than averages; their factors already include the required sample-size scaling.

This HC0 example assumes independent observations and the stated conditional mean. More generally, a sandwich can support inference for a projection coefficient under misspecification when its own moment and expansion conditions hold. It does not change that coefficient into another scientific target. Clustered data require contributions and variance calculations at the independent-unit level.

Figure 4: Left: a simulated sample from the guide’s error model. Right: population squared slope-score contributions E[(X\varepsilon)^2\mid X]=X^2(0.3+|X|)^2, compared with the homoscedastic substitution X^2E[\varepsilon^2] for standard-normal X. The sandwich “meat” retains the joint variation of leverage and error size through X_i^2\hat\varepsilon_i^2. The next simulation checks the resulting confidence-interval coverage.

5.1 Check repeated samples

The following experiment repeats that data-generating process. A nominal 95% confidence interval built with the sandwich standard error should cover more reliably than the homoscedastic interval.

OneRep <- function(n = 250) {
  x <- rnorm(n)
  y <- 1 + 0.8 * x + rnorm(n, sd = 0.3 + abs(x))
  fit <- lm(y ~ x)
  X <- model.matrix(fit); u <- residuals(fit)
  V <- solve(crossprod(X)) %*% (t(X) %*% diag(u^2) %*% X) %*% solve(crossprod(X))
  b <- coef(fit)["x"]
  c(usual = abs(b - 0.8) <= 1.96 * summary(fit)$coef["x", "Std. Error"],
    sandwich = abs(b - 0.8) <= 1.96 * sqrt(V["x", "x"]))
}
colMeans(t(replicate(1000, OneRep())))

   usual.x sandwich.x 
     0.784      0.946 

Each replication generates a new study of 250 independent people, refits the slope, and checks whether each interval contains the fixed truth 0.8. Averaging the indicators estimates repeated-sampling coverage. With 1,000 replications and coverage near 0.95, its Monte Carlo SE is about \sqrt{0.95(0.05)/1000}\approx0.007. A single simulated proportion therefore has its own numerical uncertainty.

The simulation evaluates this generator. A proof additionally verifies the expansion, CLT, and variance consistency for the estimator under stated conditions. Simulation studies for methods research explains how to build a fuller evaluation.

6 Identification comes before asymptotics

The expansion only helps when the population equation has a meaningful, unique solution in the stated parameter space. For example, E[Y-\theta^2]=\mu-\theta^2 has two roots when \mu>0 if both signs of \theta are allowed. Restricting the target to the positive square root gives a unique root. A numerical solver must respect that choice.

For the mean, the proposed target is not defined as a finite quantity if E[Y] does not exist. A distinct identification failure occurs when different parameter values give the same observable distribution or population equation. For a regression slope, it fails if X has no variation. For the partially linear model from the Conditional expectation and L2 geometry page,

\begin{aligned} \beta&= \frac{E[\widetilde X\widetilde Y]}{E[\widetilde X^2]},\\ \widetilde X&=X-E[X\mid Z], \end{aligned}

identification requires E[\widetilde X^2]>0. If X is almost fully predicted by Z, the denominator is small: the coefficient may be identified in theory but estimated with large uncertainty. This is the semiparametric version of weak information.

It is useful to separate three claims:

1. The population equation defines the scientific target or a clearly named approximation target.
2. The target is identified by the observed data distribution.
3. The sample solution is sufficiently close to the population solution for the linear expansion to hold.

These claims are often bundled together in a theorem. Read them separately.

Figure 5: For \psi=\widetilde X(\widetilde Y-\beta\widetilde X), the population slope is A=-v and the equation variance at the truth is B=v. Hence \operatorname{SE}(\hat\beta)\approx1/\sqrt{nv}: 0.10 and 0.50 at n=100. The orange equations use the population derivative and illustrative intercept perturbations of one negative standard deviation, -\sqrt{v/n}; they show a local sample approximation rather than simulated data. Noise shrinks as v falls, but sensitivity shrinks faster. At v=0, both equations are identically zero and no unique \beta is identified.

7 A smooth estimating-equation proof

Here is the argument behind the earlier approximation. Allow an approximate numerical root, with residual r_n=P_n\psi(W;\hat\theta)=o_p(n^{-1/2}). An exact root has r_n=0; for an approximate root, the numerical error must be smaller than sampling fluctuation on this scale.

Assume P\psi(W;\theta_0)=0 and \hat\theta\xrightarrow{p}\theta_0. For a scalar target, suppose the equation is differentiable near \theta_0, its population derivative is continuous there, and

\begin{aligned} A&=P\dot\psi(W;\theta_0)\\ &\ne0,\\ 0&<B\\ &=P\psi(W;\theta_0)^2\\ &<\infty. \end{aligned}

Also require a uniform derivative LLN on a neighborhood \mathcal N of \theta_0:

\sup_{\theta\in\mathcal N} |(P_n-P)\dot\psi(W;\theta)|\xrightarrow{p}0.

“Uniform” means the sample and population derivatives are close throughout that neighborhood, including at a random intermediate value selected using the data. A fixed-parameter LLN alone does not establish this.

These are conditions to verify for the problem. Consistency itself can follow from a uniform LLN for the equation and a separated unique population root, with an appropriate parameter set and root-selection rule. Separation means values outside any fixed neighborhood of \theta_0 cannot have population equation values arbitrarily close to zero. The LLN then prevents a well-chosen approximate sample root from remaining outside that neighborhood.

A Taylor expansion uses an intermediate value \tilde\theta between the consistent estimate and target:

\begin{aligned} P_n\psi(\hat\theta)&=P_n\psi(\theta_0)+\widetilde A_n(\hat\theta-\theta_0),\\ \widetilde A_n&=P_n\dot\psi(\tilde\theta)\\ &\xrightarrow{p}A. \end{aligned}

Consistency keeps \tilde\theta near \theta_0 with probability tending to one. Uniform derivative control and population continuity therefore give \widetilde A_n\xrightarrow{p}A\ne0. Rearranging the equality shows both sources of error:

\hat\theta-\theta_0=-\widetilde A_n^{-1}P_n\psi(\theta_0) +\widetilde A_n^{-1}r_n.

The CLT makes P_n\psi(\theta_0)=O_p(n^{-1/2}). Continuous mapping gives \widetilde A_n^{-1}-A^{-1}=o_p(1), so replacing the intermediate inverse by A^{-1} changes the expression by o_p(n^{-1/2}). The numerical-root term has that same negligible order because r_n=o_p(n^{-1/2}) and the inverse is O_p(1). Hence

\sqrt n(\hat\theta-\theta_0)=\sqrt nP_nD_0+o_p(1).

The fixed-function CLT handles the leading average; Slutsky removes the negligible remainder. Thus the influence function is -A^{-1}\psi(W;\theta_0) and the limiting variance is B/A^2.

The intermediate \widetilde A_n establishes the proof and need not be observable. For computation, use \widehat A=P_n\dot\psi(W;\hat\theta) and \widehat B=P_n\psi(W;\hat\theta)^2 from the sandwich section. Their consistency needs the corresponding local LLN and continuity arguments. These complete the inference proof after root consistency and linearization.

For a vector equation, use the Jacobian integrated along the segment from target to estimate, with convergence to a nonsingular matrix A. A single intermediate point for every component is not generally a valid multivariate mean-value argument.

The smooth-root checkpoint SP-01 supplies a concrete example with \hat\theta=\sqrt{\bar Y}. Its exact factorization makes the negligible-remainder argument especially transparent. Quantiles and local expansions is a later transfer that handles a nonsmooth sample equation with a different local argument.

8 What the theorem hides

The Taylor expansion uses conditions that can fail.

• The observations must have enough independence for an average to obey a central limit theorem. Repeated observations, time series, and related samples need a different variance calculation or a different estimating equation.
• A must be nonzero in one dimension, or invertible in several dimensions. Low equation sensitivity relative to its random fluctuation can produce unstable estimates. Interpret the derivative together with B: multiplying the whole equation by a constant changes both without changing the estimator or its variance.
• The estimating function must be smooth enough near the truth. Quantiles, thresholds, and selection procedures need extra work because they are not ordinary differentiable equations.
• The remainder in the Taylor expansion must be smaller than n^{-1/2}. In semiparametric work that remainder includes error from estimating nuisance functions.

The point is not to memorize a list of technical assumptions. It is to know what a theorem is buying: a valid first-order approximation, and therefore a standard error and confidence interval.

9 Where this connects

• The Conditional expectation and L2 geometry page gives a population residualization identity. Its sample analogue becomes an estimating equation.
• The Likelihood, scores, and parametric efficiency page is an important special case: a maximum likelihood estimator solves the score equation.
• The Influence functions page generalizes the linear expansion from a chosen estimating equation to smooth functionals of a distribution.
• The Orthogonality and cross-fitting page studies the remainder when \psi contains estimated nuisance functions.

Use the core theorems and stochastic-order definitions in Large-sample reasoning for methods research to identify each step of the proof above.

10 Exercises

Checkpoints are completion targets for this page. Additional practice below is optional. Each checkpoint gives a completion standard, a hint, and a solution outline.

Checkpoint · SP-01 — Complete a smooth estimating-equation argument

Task type: Derivation and asymptotic reasoning.

Let Y_i>0 be i.i.d., with E[Y]=\mu>0 and 0<\operatorname{Var}(Y)=\sigma^2<\infty. The target is \theta_0=\sqrt\mu, and \hat\theta=\sqrt{\bar Y} solves P_n(Y-\theta^2)=0.

1. Prove consistency using the law of large numbers.
2. Use the identity (\hat\theta-\theta_0)(\hat\theta+\theta_0)=\bar Y-\mu to derive an expansion with an o_p(n^{-1/2}) remainder.
3. Identify A, the influence function, and the limiting distribution of \sqrt n(\hat\theta-\theta_0).
4. Give a consistent standard-error estimate. Explain which steps use independence, finite variance, and \mu>0.

Ready to continue when: you can connect the population root to consistency, a negligible remainder, a CLT, and an estimated standard error.

Transfer: Apply the same reasoning to a smooth function g(\bar Y). What derivative multiplies Y-\mu, and what would require rechecking if \mu approached a boundary?

Hint

The denominator \hat\theta+\theta_0 converges to 2\theta_0>0. Combine that with \bar Y-\mu=O_p(n^{-1/2}).

Solution outline — open after trying

The LLN and continuity give \hat\theta\to_p\theta_0. Then

\begin{aligned} \hat\theta-\theta_0&=\frac{\bar Y-\mu}{2\theta_0}+R_n,\\ R_n&=(\bar Y-\mu)D_n,\\ D_n&=\frac1{\hat\theta+\theta_0}-\frac1{2\theta_0}. \end{aligned}

Here D_n=o_p(1) and \bar Y-\mu=O_p(n^{-1/2}) by the CLT. Their product R_n is o_p(n^{-1/2}). Here A=-2\theta_0, so \phi(Y)=(Y-\mu)/(2\theta_0) and the limiting variance is \sigma^2/(4\mu). A standard error is s_Y/(2\hat\theta\sqrt n). The i.i.d. assumption and finite variance support the LLN, CLT, and variance consistency; positivity keeps the root away from the derivative-zero boundary. Slutsky’s theorem justifies replacing unknown quantities with consistent estimates.

Additional practice (optional)

These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.

1. Write the estimating equation for the sample variance with divisor n, treating both the mean and variance as unknown parameters. Derive its sandwich variance.
2. Show that the least-squares normal equation identifies the coefficient \beta^\ast = E[XX^\top]^{-1}E[XY] when E[XX^\top] is invertible. What does it target if the conditional mean of Y is nonlinear?
3. Repeat the coverage simulation with n=50, 250, and 1000. Report the Monte Carlo standard error for each estimated coverage probability.
4. In the partially linear model, let X=Z exactly. Which step in the identification argument fails? Simulate X=Z+\sigma U for several small values of \sigma and describe what happens to the standard error of \hat\beta.

11 Reading

• van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge. Chapters 5 and 7 give the general Z-estimation argument and M-estimation examples.
• Huber, P. J. (1967). The behavior of maximum likelihood estimates under nonstandard conditions. In Proceedings of the Fifth Berkeley Symposium. An early source for the general sandwich result.
• White, H. (1982). Maximum likelihood estimation of misspecified models. Econometrica, 50(1), 1 to 25. The classic treatment of pseudo-true targets and robust variance under misspecification.
• Stefanski, L. A. and Boos, D. D. (2002). The calculus of M-estimation. The American Statistician, 56(1), 29 to 38. A compact and practical introduction.