1 Foundations of estimating equations

1.1 Defining estimating functions

An estimating equation is built from an estimating function \(g(\theta)\) whose expectation under the true data-generating mechanism equals zero when \(\theta\) is the true parameter. The core idea is to replace direct likelihood maximization with the solution of \[ \mathbb{E}\big[g(Z,\theta_0)\big]=0, \] where \(Z\) denotes a generic observation and \(\theta_0\) is the parameter of interest. In practice, \( \mathbb{E}[\cdot]\) is replaced by an empirical average, yielding sample moment conditions that are solved for \(\theta\).

Estimating functions can be derived from known moment identities, from unbiased score-like constructions, or from general principles such as balancing moment restrictions that are expected to hold at the true parameter.

1.2 Moment conditions and unbiasedness

Moment conditions express restrictions on expected values of functions of the data and parameters. Unbiasedness refers to the property that the estimating function has mean zero at the true parameter: \[ \mathbb{E}[g(Z,\theta_0)] = 0. \] This property ensures that the target parameter is a root of the population estimating equation. When the model is correctly specified, such unbiased estimating functions often arise naturally. Under misspecification, unbiasedness may hold only for certain “pseudo-true” parameters (parameters defined by the moment restrictions rather than by an exact likelihood model).

1.3 Identifiability and solvability

Solvability requires that the estimating equation has at least one solution in the parameter space, while identifiability requires that the solution corresponds to a unique meaningful parameter. In population terms, identifiability often corresponds to conditions such as:

  • the population moment map \(\theta \mapsto \mathbb{E}[g(Z,\theta)]\) equals zero only at \(\theta_0\);
  • the mapping is sufficiently regular so that small changes in \(\theta\) alter the moment expectation in a non-degenerate way.

If the equation has multiple roots or near-flat regions, estimation can become unstable, sensitive to initialization, or inconsistent without additional structure.

1.4 Relationship to method of moments and M-estimation

Estimating equations generalize classical method of moments: instead of matching a finite set of moments exactly, one solves sample-based versions of moment equalities. They also relate to M-estimation by viewing the estimator as the solution to a first-order condition. Many M-estimators can be rewritten as estimating equation solutions using a gradient (score) of a criterion function, while estimating equations can be more flexible because they do not require an explicit likelihood or objective function.

More broadly, estimating equations sit within a unified framework where consistency follows from matching expected moment conditions, and asymptotic distributions follow from smoothness and limit theorems.

1.5 Notation, parameter vectors, and constraints

Let \(\theta \in \Theta \subset \mathbb{R}^p\) be a parameter vector, possibly with constraints such as positivity, boundedness, or membership in a manifold (e.g., correlations). The estimating equation is typically written as \[ \frac{1}{n}\sum_{i=1}^n g(Z_i,\theta)=0 \] for independent observations, with corresponding modifications for dependence. Constraints are handled either by choosing a constrained parameterization or by applying numerical methods that respect the feasible region. In constrained settings, uniqueness and asymptotic normality may require additional regularity conditions tied to the geometry of \(\Theta\).

2 Core setup and solution methods

2.1 Formulating the estimating equation

2.1.1 Choice of estimating function g(θ)

Selecting \(g(\theta)\) is central to the procedure. Good estimating functions often satisfy:

  • unbiasedness at the target parameter;
  • differentiability in \(\theta\) to support asymptotic theory;
  • robustness properties relevant to the application (e.g., bounded influence, stable behavior under outliers);
  • computational convenience, since solving may require repeated evaluation.

In generalized contexts, \(g(\theta)\) may combine multiple moment conditions into a vector so that the resulting system can be solved simultaneously.

2.1.2 Sample version and empirical averages

Given data \(\{Z_i\}_{i=1}^n\), the sample estimating equation takes the form \[ \hat{G}_n(\theta)=\frac{1}{n}\sum_{i=1}^n g(Z_i,\theta), \] and the estimator \(\hat{\theta}\) solves \(\hat{G}_n(\hat{\theta})=0\) (or, in overidentified cases, approximately minimizes a distance to zero after weighting). Empirical averages are preferred because they preserve the moment structure while enabling asymptotic arguments using laws of large numbers.

In dependent settings (e.g., repeated measures), the sample analogue may aggregate over clusters or time blocks using a correlation-aware structure or cluster-robust replacements.

2.2 Solving for parameter estimates

2.2.1 Root-finding approaches

When the estimating function is square (dimension of moments equals dimension of parameters), the task is to find a root of \(\hat{G}_n(\theta)\). Common approaches include:

  • Newton–Raphson: using a local linear approximation based on the Jacobian;
  • quasi-Newton methods: updating an approximate Jacobian;
  • secant or derivative-free methods when gradients are unstable or unavailable.

In practice, the choice balances speed, stability, and the difficulty of evaluating derivatives.

2.2.2 Iterative algorithms and numerical stability

Iterative solvers must manage numerical issues such as ill-conditioning, step sizes that lead outside the feasible region, and sensitivity to scaling of parameters. Stability is improved by:

  • using line searches or trust regions;
  • scaling parameters so that different components have comparable magnitudes;
  • monitoring convergence of both the estimating equation norm and the parameter update;
  • ensuring derivative calculations are accurate, especially for complex estimating functions.

When the estimating system is nearly singular, regularization or alternative parameterizations may be required.

2.3 Existence, uniqueness, and initialization

Because sample estimating equations are random, existence and uniqueness are not guaranteed even if the population equation has a unique root. Initialization matters: some systems have multiple local solutions, and numerical methods can converge to different roots depending on starting values. Convergence diagnostics often include:

- checking whether the final estimate yields a small residual \(\|\hat{G}_n(\hat{\theta})\|\);
  • verifying that the solution lies in the intended parameter space;
  • comparing multiple initializations to assess whether the solver finds a consistent root.

2.4 Handling nuisance parameters

Many models contain nuisance components that are not directly of interest but affect moment conditions. Nuisance handling can be done by:

  • joint estimation: augment the parameter vector so the estimating equation system includes nuisance parameters;
  • profiling: solving nuisance equations conditional on \(\theta\);
  • plugging in consistent nuisance estimators obtained from auxiliary procedures.

The asymptotic impact depends on whether nuisance estimation is root-\(n\) consistent and how it interacts with the estimating function. In semiparametric settings, nuisance may be estimated nonparametrically, requiring careful conditions to preserve asymptotic normality of the parameter of interest.

3 Asymptotic theory

3.1 Consistency of solutions

Consistency typically follows from a combination of:

  • a uniform law of large numbers ensuring \(\hat{G}_n(\theta)\) converges to \(G(\theta)=\mathbb{E}[g(Z,\theta)]\);
  • identifiability, ensuring \(G(\theta)=0\) only at \(\theta_0\);
  • regularity conditions that prevent the estimator from drifting to other roots or boundaries.

For approximate solutions (e.g., minimizing a weighted criterion), consistency is tied to the criterion’s limiting behavior and the stability of the minimizer.

3.2 Asymptotic normality

Under standard smoothness and moment conditions, the estimator admits an asymptotic linear expansion: \[ \sqrt{n}(\hat{\theta}-\theta_0)\ \Rightarrow\ \mathcal{N}(0,\Sigma), \] where \(\Sigma\) is derived from the variability of the estimating functions and the sensitivity of the moment map to parameter changes.

3.2.1 Linearization and influence-function style expansions

Asymptotic normality is commonly obtained via a first-order Taylor expansion around \(\theta_0\): \[ 0=\hat{G}_n(\hat{\theta}) \approx \hat{G}_n(\theta_0)+ \dot{G}(\theta_0)(\hat{\theta}-\theta_0), \] leading to \[ \hat{\theta}-\theta_0 \approx -\dot{G}(\theta_0)^{-1}\hat{G}_n(\theta_0), \] where \(\dot{G}(\theta_0)\) denotes the Jacobian (sensitivity matrix). This yields an influence-function-like representation, revealing that the asymptotic distribution is driven by the sum of centered estimating functions.

3.3 Role of the Jacobian (sensitivity matrix)

The Jacobian \(\dot{G}(\theta)\) quantifies how the expectation of the estimating functions changes with the parameter. If \(\dot{G}(\theta_0)\) is singular or nearly singular, the estimator may have inflated variance, slow convergence, or failure of asymptotic normality. Hence, regularity conditions typically require nonsingularity and boundedness of derivatives in a neighborhood of \(\theta_0\).

3.4 Central limit theorem conditions

To derive asymptotic normality, one needs a central limit theorem for the empirical mean of \(g(Z_i,\theta_0)\). For independent observations, this typically requires finite second moments and possibly Lindeberg-type conditions. For dependent observations, one uses mixing or martingale difference conditions, or cluster-level central limit theorems, depending on the dependence structure.

3.5 Rates of convergence and practical implications

When regularity holds, the convergence rate is typically \(n^{-1/2}\). However, in boundary problems, with weak identification, or under nonregular behavior (e.g., parameters constrained to lie on the boundary), the rate and limiting distribution can differ. Practical implications include:

  • caution in interpreting standard errors when the Jacobian is poorly conditioned;
  • additional diagnostics for small samples or highly nonlinear estimating functions;
  • potential need for robust methods or resampling when asymptotic approximations are unreliable.

4 Variance estimation and inference

4.1 Sandwich (robust) covariance estimator

A widely used covariance estimator for estimating equations is the sandwich form: \[ \widehat{\mathrm{Var}}(\hat{\theta})=\dot{G}_n(\hat{\theta})^{-1}\,\hat{S}_n(\hat{\theta})\,\dot{G}_n(\hat{\theta})^{-\top}, \] where \(\dot{G}_n(\hat{\theta})\) estimates the Jacobian and \(\hat{S}_n(\hat{\theta})\) estimates the variability of \(g(Z,\theta)\). This approach is robust to certain forms of model misspecification because it does not require the full specification that would be needed for likelihood-based Fisher information formulas.

4.2 Estimating components in the sandwich form

The sensitivity matrix \(\dot{G}_n(\hat{\theta})\) can be computed by analytical derivatives, numerical differentiation, or automatic differentiation. The middle matrix \(\hat{S}_n(\hat{\theta})\) is often estimated as the empirical covariance of \(g(Z_i,\hat{\theta})\): \[ \hat{S}_n(\hat{\theta})=\frac{1}{n}\sum_{i=1}^n (g(Z_i,\hat{\theta})-\bar{g}(\hat{\theta}))(g(Z_i,\hat{\theta})-\bar{g}(\hat{\theta}))^\top, \] with modifications for dependence (e.g., cluster aggregates) so that the covariance reflects the sampling dependence structure. Numerical stability can be improved by ensuring proper centering and by avoiding poorly conditioned covariance matrices through regularization when necessary.

4.3 Wald-type tests and confidence intervals

With an estimate of the covariance matrix, inference for a scalar or linear contrast \(a^\top \hat{\theta}\) uses Wald methods. A typical \(100(1-\alpha)\%\) confidence interval for \(a^\top\theta_0\) is \[ a^\top\hat{\theta}\ \pm\ z_{1-\alpha/2}\sqrt{a^\top \widehat{\mathrm{Var}}(\hat{\theta}) a}. \] Wald tests evaluate whether a constraint \(H_0: a^\top\theta_0=0\) holds, based on the standardized contrast. The quality of these procedures depends on the validity of the asymptotic approximation and the quality of the covariance estimator.

4.4 Score-type and likelihood-ratio analogues

Although estimating equations do not rely on likelihood directly, analogous inference procedures exist. Score-type methods may be implemented via generalized score or estimating-function-based tests, where one evaluates standardized moment residuals under the null and uses covariance estimates to calibrate the statistic. Likelihood-ratio analogues can be approximated in frameworks such as empirical likelihood or through quasi-likelihood constructions, but the correspondence is less direct than in pure likelihood theory.

4.5 Small-sample adjustments

In finite samples, sandwich estimators can be biased due to the estimation of \(\dot{G}\) and the covariance structure. Adjustments include:

  • degrees-of-freedom corrections in covariance estimation;
  • using \(m\)-out-of-\(n\) bootstrap or subsampling to better approximate sampling variability;
  • employing bias-reduced variance formulas or penalized corrections in special frameworks (notably for empirical likelihood).

The appropriateness of adjustments depends on the dependence structure, sample size, and the degree of overidentification.

5 Special classes and connections

5.1 Generalized linear models as estimating equations

Generalized linear models can be expressed as estimating equations derived from quasi-score functions. In canonical settings, the estimating function often involves residuals weighted by the derivative of the mean with respect to the linear predictor. Solving these equations yields the familiar iteratively reweighted least squares algorithm. Even when full likelihood assumptions are relaxed, the estimating equation formulation remains useful for producing consistent estimates under weaker conditions.

5.2 Quasi-likelihood and GEE (generalized estimating equations)

For correlated data such as longitudinal measurements, generalized estimating equations extend the estimating-equation approach by incorporating a working correlation structure. The estimating function typically reflects a marginal mean model, while the covariance of the estimating functions is adjusted using either the working correlation or robust “empirical” covariance. This yields consistency for regression coefficients under broad conditions, with efficiency improved when the working correlation is close to the truth.

5.3 Generalized method of moments (GMM)

In generalized method of moments, one may have more moment conditions than parameters (overidentification). Estimation is performed by minimizing a weighted quadratic form of the sample moments, producing an estimate that depends on the choice of weighting matrix. Under standard regularity and correct moment conditions, GMM can achieve asymptotic efficiency when the weighting matrix is chosen as the inverse of the long-run covariance of the moment conditions. Estimating equations appear as special cases when the system is exactly identified.

5.4 Semiparametric estimating equations

Semiparametric settings include models where parts of the data-generating mechanism remain unspecified or are modeled flexibly. Estimating equations are constructed to target parameters in the presence of nuisance functions. Robustness often comes from using moment conditions designed to be orthogonal to nuisance estimation errors, enabling valid inference even when nuisance components are learned with complex methods. Asymptotic analysis then relies on conditions ensuring nuisance estimation does not dominate the asymptotic variability of the parameter of interest.

Empirical likelihood replaces likelihood with a nonparametric likelihood constructed from moment constraints. Estimation proceeds by maximizing an empirical likelihood subject to constraints that the sample moments match their theoretical expectations. This yields inference tools that can mimic likelihood-based procedures while remaining tied to moment restrictions. Empirical likelihood is connected to estimating equations through the duality between constraint satisfaction and Lagrange multipliers, though the exact form of resulting estimators and tests depends on the specific implementation.

6 Dependence, clustering, and repeated measures

6.1 Independent vs correlated observations

Standard estimating-equation theory assumes independent contributions, but many applications involve dependence. Dependence affects the variance of the estimating functions and can invalidate naive covariance formulas that assume independence. To correct this, variance estimation must reflect the actual sampling structure, either by using robust covariance estimators or by modeling the dependence in the estimating function design.

6.2 Cluster-robust estimating equations

When observations cluster (e.g., multiple measurements per subject), cluster-robust methods treat within-cluster dependence as potentially arbitrary while assuming independence across clusters. The estimating function contributions are aggregated at the cluster level, and the covariance is estimated using cluster sums. This approach yields valid large-sample inference for parameters under appropriate conditions even when the within-cluster correlation is misspecified.

6.3 Working correlation in GEE

In GEE, a working correlation structure is used to improve efficiency. Even if the working correlation is incorrect, regression coefficients can remain consistent when the mean model is correct, provided robust covariance is employed. When the working correlation matches the true dependence pattern, variance estimates become smaller and inference more powerful.

6.4 Robustness to misspecification

The estimating-equation framework supports robust inference through sandwich covariance and related adjustments. Robustness is not universal: it depends on what is misspecified (mean model, variance structure, dependence model) and on how moment conditions are constructed. Typically, the most reliable results are obtained when the estimating function is unbiased for the target parameter under the data-generating mechanism, even if auxiliary assumptions are violated.

6.5 Efficiency considerations

Efficiency depends on the choice of estimating function and weighting. Optimal or near-optimal choices can reduce estimator variance, but they may require stronger modeling assumptions. Robust procedures trade some efficiency for protection against misspecification, often leading to larger standard errors but improved credibility of inference across scenarios.

7 Practical workflow and diagnostics

7.1 Model checking for moment conditions

Because estimation is based on moment restrictions, diagnostics often focus on whether moment conditions appear plausible in the data. Common strategies include:

  • examining empirical moment residuals at the fitted parameter;
  • checking whether sample moments behave stably across subgroups or time blocks;
  • using plots of estimating-function components versus covariates to detect systematic deviations.

Moment checks are not the same as likelihood goodness-of-fit, but they provide targeted evidence about whether the assumptions driving the estimating equations are reasonable.

7.2 Sensitivity to outliers and leverage

Estimating equations may be sensitive if the estimating function grows quickly with certain covariate values or outcomes. Outliers can distort the sample moments and shift roots. Diagnostics include:

  • influence analysis (approximating how \(\hat{\theta}\) changes when an observation or cluster is removed);
  • leverage inspection using norms of covariate vectors or fitted residual magnitudes;
  • robust variants of \(g(\theta)\) that downweight extreme observations.

7.3 Tuning, convergence, and stopping criteria

Numerical solvers require termination rules, such as:

  • small norm of \(\hat{G}_n(\theta)\);
  • small change in \(\theta\) between iterations;
  • bounded gradient or step size thresholds.

Tuning may also involve selecting weighting matrices (e.g., in GMM), regularization parameters (in regularized estimating equations), or constraints handling methods. Poor tuning can lead to slow convergence or convergence to undesired roots.

7.4 Missing data and estimating equations

Missing data can be addressed within the estimating-equation framework by using:

  • complete-case estimating equations, with potential loss of efficiency;
  • inverse-probability weighting when missingness is modeled through selection probabilities;
  • doubly robust constructions that combine outcome modeling and propensity modeling within moment conditions.

Validity hinges on missingness assumptions and on whether the estimating function remains unbiased for the target parameter under the missingness mechanism being used.

7.5 Diagnostic plots and residual-based summaries

Although “residuals” depend on the moment formulation, useful summaries often include:

  • component-wise plots of \(g(Z_i,\hat{\theta})\) to locate systematic patterns;
  • residual-type quantities aligned with the chosen moment conditions (e.g., standardized working residuals);
  • sandwich-based diagnostics that track standardized contributions to the score-like components.

These tools help detect heterogeneity, nonlinearity, and dependence misalignment.

8 Extensions and advanced topics

8.1 Regularized estimating equations

When estimating equations are high-dimensional or ill-posed, regularization can stabilize the solution. Approaches include penalizing norms of \(\theta\), adding ridge-type penalties in linearized systems, or using constrained optimization to control variance inflation. Regularization changes both finite-sample behavior and asymptotic theory, requiring modified proofs and potentially different limiting distributions.

8.2 High-dimensional parameter settings

In high-dimensional regimes where the number of parameters grows with sample size, classical asymptotic normality may fail or require stronger conditions. Modern extensions use sparsity assumptions, debiasing strategies, and careful control of empirical process terms. In this context, estimating equations often become part of a two-stage procedure: selection or regularized estimation followed by correction to recover valid inference.

8.3 Empirical process considerations

Asymptotic theory in advanced settings relies on controlling fluctuations of the estimating function over neighborhoods of \(\theta\). This can involve empirical process tools such as uniform convergence, covering numbers, or bracketing arguments. These techniques ensure that sample-based moment maps approximate their population counterparts uniformly, a key step for consistency and for asymptotic linearization.

8.4 Combining multiple estimating equations

One may combine distinct sets of estimating functions capturing different aspects of the model (e.g., separate moment conditions for mean and variance). Combination can be done by stacking moment equations, or by forming weighted sums that prioritize certain conditions. Care is needed to avoid redundancy that can cause near collinearity, and to account for correlations among moment conditions in the variance estimator.

8.5 Optimal weighting and efficiency gains

Efficiency gains often come from choosing weighting matrices that minimize asymptotic variance, especially in GMM or overidentified systems. In practice, optimal weights are unknown and must be estimated, leading to feasible two-step or iterated GMM procedures. When dependence is present, “long-run” covariance estimation plays the role of the optimal weighting matrix, improving efficiency while maintaining robust behavior under misspecification.

9 Examples and worked applications

9.1 Estimating equations for means and variances

A basic example uses estimating functions built from centered moments. For instance, to estimate a mean \(\mu\), one may solve \[ \sum_{i=1}^n (Y_i-\mu)=0. \] To estimate both mean and variance, one can solve a system using both first and second centered moments, for example \[ \sum_{i=1}^n (Y_i-\mu)=0,\quad \sum_{i=1}^n \big((Y_i-\mu)^2-\sigma^2\big)=0. \] This illustrates how multiple moment conditions yield a multivariate estimating equation and how the Jacobian affects standard error computation.

9.2 Estimation under generalized linear models

In generalized linear models, estimating equations typically enforce a mean-structure relationship between predictors and the conditional expectation of the response. The estimating function is constructed from residuals and the derivative of the mean with respect to the linear predictor. Solving the resulting equations reproduces familiar coefficient estimates, while sandwich covariance provides robust standard errors when distributional assumptions (such as variance specification) are imperfect.

9.3 Longitudinal modeling with GEE

For repeated measurements, GEE uses a marginal mean model and constructs estimating functions that account for within-subject dependence through a working correlation. A typical workflow specifies:

  • the mean model for each time point;
  • a working correlation (exchangeable, autoregressive, independence);
  • robust covariance estimation to protect against correlation misspecification.

The resulting coefficients estimate population-average effects, and the working correlation primarily influences efficiency rather than consistency under common assumptions.

9.4 GMM examples with instrumental variables

Instrumental-variable settings provide a canonical GMM example: if \(X\) is endogenous, one can use instruments \(Z\) to impose moment restrictions such as \[ \mathbb{E}[Z(Y-\beta X)]=0. \] The sample moments are formed from the product of instruments and residuals, and the parameter \(\beta\) is obtained by minimizing a weighted criterion. This demonstrates how overidentification and weighting choices affect efficiency, and how robust covariance estimates support valid inference when certain distributional assumptions are relaxed.

9.5 Comparing results with likelihood-based methods

In many models, likelihood-based estimators correspond to particular estimating equations (often involving score functions). Comparing results highlights:

  • when likelihood and estimating-equation approaches agree closely, both the model and the moment conditions are consistent with each other;
  • when they diverge, the divergence can reflect misspecification in variance structure, correlation modeling, or departures from distributional assumptions.

The estimating-equation approach is often favored for its flexibility and for providing robust standard errors through sandwich formulas.