Method of Moments & M-Estimation

The first systematic use of moment matching appears in Pearson’s 1894 paper “Contributions to the Mathematical Theory of Evolution,” where he fit a mixture of two Normals to the carapace-width measurements of 1000 crabs from the Bay of Naples. The mixture has five parameters — two means, two variances, and a mixing weight — so Pearson computed five sample moments and solved the resulting nonlinear system, by hand, to recover the underlying components. He used the method because no other recipe existed. The likelihood-maximization principle was still nearly thirty years away.

The reason MoM predates the MLE is operational. Until R. A. Fisher’s 1922 paper “On the Mathematical Foundations of Theoretical Statistics,” likelihood was viewed as one of several plausible criteria for fitting a distribution; there was no general argument for preferring it. What Fisher established — that the MLE is consistent, asymptotically Normal, and asymptotically efficient under regularity — required machinery (Fisher information, the Cramér–Rao bound, asymptotic theory) that did not exist in Pearson’s era. So a calculation that could be done won out over one that would have required theory not yet developed. The lesson is durable: even today, when the MLE is intractable, MoM is the first thing to reach for.

For an iid sample $X_1, \ldots, X_n$, the $k$-th sample moment is

The first sample moment $\overline{X^1} = \bar{X}_n$ is the sample mean (Definition 13.2). By the strong law of large numbers, $\overline{X^k} \xrightarrow{a.s.} \mu_k = \mathbb{E}[X^k]$ whenever the population moment $\mu_k$ exists.

Let $X_1, \ldots, X_n$ be iid from a parametric family with $\theta = (\theta_1, \ldots, \theta_k) \in \Theta \subseteq \mathbb{R}^k$, and write $\mu_j(\theta) = \mathbb{E}_\theta[X^j]$ for the $j$-th population moment as a function of $\theta$. Suppose the moment map

is invertible on its image. The method-of-moments estimator $\hat{\theta}_{\text{MoM}}$ is the (unique) solution to the system of $k$ equations

Equivalently, $\hat{\theta}_{\text{MoM}} = g^{-1}\bigl(\overline{X^1}, \ldots, \overline{X^k}\bigr)$.

If the moment map $g: \Theta \to \mathbb{R}^k$ is continuous and one-to-one on $\Theta$, and $(\overline{X^1}, \ldots, \overline{X^k})$ lies in the image $g(\Theta)$, then the MoM estimator exists and is unique. If $g$ is continuously differentiable with a non-singular Jacobian $\nabla g(\theta_0)$ at the true parameter $\theta_0$, then the inverse $g^{-1}$ is differentiable in a neighborhood of $g(\theta_0)$, and $\hat{\theta}_{\text{MoM}}$ is well-defined for sufficiently large $n$ with probability tending to one.

For $X \sim \mathcal{N}(\mu, \sigma^2)$, the first two population moments are

Setting $\overline{X^1} = \hat{\mu}$ and $\overline{X^2} = \hat{\mu}^2 + \hat{\sigma}^2$ and solving:

and

The MoM variance estimator is the biased $1/n$ form (the same form the MLE uses, see Example 14.2), not the Bessel-corrected $1/(n-1)$ form. This is structural: MoM matches the population variance to the raw second-moment-minus-squared-first-moment, whose unbiased correction is a $(n-1)/n$ shift. We will see in §15.7 that for the Normal variance, $\text{ARE}(\hat{\sigma}^2_{\text{MoM}}, \hat{\sigma}^2_{\text{MLE}}) = 1$ — they are the same estimator, both dominated in MSE by the unbiased $S^2_{n-1}$.

For $X \sim \text{Exp}(\lambda)$, $\mathbb{E}[X] = 1/\lambda$. So $\bar{X}_n = 1/\hat{\lambda}_{\text{MoM}}$, giving

This is identical to the MLE (Example 14.7). The reason — to be made precise in Theorem 5 of §15.7 — is that the Exponential is an exponential family in its natural parameter, and for such families the MLE is a moment-matching estimator. The two recipes coincide.

For $X \sim \text{Poisson}(\lambda)$, $\mathbb{E}[X] = \lambda$. So

again identical to the MLE. The Poisson is also a natural exponential family. The pattern will become explicit in §15.7.

The construction in Definition 2 uses exactly $k$ moment equations to identify $k$ parameters. If we attempted to use more — say, the first three sample moments to estimate two parameters — the system is over-determined and generally has no exact solution. Choosing which two equations to enforce, or how to weight all three, is the move that takes us from MoM to the generalized method of moments (GMM, Hansen 1982). We preview GMM in §15.12. For the rest of this topic, we hold to the exact-identification setting: $k$ moments for $k$ parameters.

For $X \sim \text{Uniform}(0, \theta)$, $\mathbb{E}[X] = \theta/2$. So $\bar{X}_n = \hat{\theta}_{\text{MoM}}/2$, giving

This is a perfectly reasonable estimator — consistent, with asymptotic variance $\theta^2/(3n)$ from the delta method. But it has two pathologies that the MLE does not. First, with positive probability $\hat{\theta}_{\text{MoM}} < X_{(n)}$, where $X_{(n)} = \max_i X_i$ is the largest observation. That is impossible — the parameter $\theta$ must be at least as large as every observed value, so $\hat{\theta} < X_{(n)}$ is infeasible. Second, the MLE for this problem is

with asymptotic variance of order $\theta^2/n^2$ — a full order of magnitude better than MoM’s $\theta^2/(3n)$. The ratio $\text{ARE}(\hat{\theta}_{\text{MoM}}, \hat{\theta}_{\text{MLE}}) = 3/(n+2) \to 0$ as $n \to \infty$.

We do not call this a failure of MoM — MoM gives an answer, and the answer converges to $\theta$. We call it the most graphic illustration of why efficiency matters: when the MLE achieves a faster convergence rate (here $1/n$ vs. $1/\sqrt{n}$), the asymptotic ratio of variances is not a fixed number but a sequence going to zero.

For $X \sim \text{Bernoulli}(p)$, $\mathbb{E}[X] = p$, so

The Bernoulli is yet another natural exponential family, so MoM = MLE. Three for three so far in the exponential-family column.

The Geometric distribution comes in two conventions; we adopt the failures-before-first-success convention with $X \in \{0, 1, 2, \ldots\}$ and $\mathbb{P}(X = k) = (1-p)^k p$. Then $\mathbb{E}[X] = (1-p)/p$, and setting $\bar{X}_n = (1-\hat{p})/\hat{p}$ gives

If instead we use the trials-until-first-success convention with $X \in \{1, 2, \ldots\}$, then $\mathbb{E}[X] = 1/p$ and $\hat{p}_{\text{MoM}} = 1/\bar{X}_n$. Both are exponential families; both coincide with the MLE.

Example 4 hinted at it; the issue is general. The MoM estimate is the inverse moment map applied to sample moments, which can land anywhere in $\mathbb{R}^k$ — not just in $g(\Theta)$. A few realistic ways this bites:

• Negative variance. For a Gamma$(\alpha, \beta)$ sample, $\hat{\alpha}_{\text{MoM}} = \bar{X}^2/S^2$ requires $S^2 > 0$. With probability one this holds for $n \geq 2$ and a non-degenerate sample, but a near-constant sample (small $S^2$) makes the estimate explode.
• Probability outside [0, 1]. For a Beta$(\alpha, \beta)$, the MoM solves $\hat{\alpha} = \bar{X} \cdot c$ where $c = \bar{X}(1-\bar{X})/S^2 - 1$; if $c \le 0$ the formal MoM gives $\hat{\alpha} \le 0$, outside the parameter space.
• Mixing weight outside [0, 1]. Pearson’s two-Normal mixture famously yielded mixing weights below 0 and above 1 on certain datasets — Pearson reported this as a feature of the data (“the model is mis-specified”) rather than a bug of MoM.

The standard fix is to project onto the constraint set, but at that point we have abandoned pure MoM. Better diagnostic: when MoM gives an absurd answer, suspect the model.

For a model with $k$ parameters $\theta = (\theta_1, \ldots, \theta_k) \in \Theta \subseteq \mathbb{R}^k$, the moment map $g: \Theta \to \mathbb{R}^k$ is the vector of population moments,

Its Jacobian at $\theta$ is the $k \times k$ matrix

The MoM estimator is $\hat{\theta}_{\text{MoM}} = g^{-1}(\overline{X^1}, \ldots, \overline{X^k})$, well-defined when $\nabla g(\theta_0)$ is invertible at the true parameter.

Let $X_1, \ldots, X_n$ be iid $\text{Gamma}(\alpha_0, \beta_0)$ with $\alpha_0, \beta_0 > 0$, and let $\hat{\alpha}_{\text{MoM}} = \bar{X}_n^2/S^2_n$, $\hat{\beta}_{\text{MoM}} = \bar{X}_n/S^2_n$. Then

and

where $\Sigma_{\text{moments}}$ is the limiting covariance of $\sqrt{n}\,(\bar{X}_n - \mu_1, \overline{X^2} - \mu_2)$ given by the multivariate CLT, and $\nabla g$ is the Jacobian of the moment map.

Restating from the derivation above: for $X_1, \ldots, X_n$ iid $\text{Gamma}(\alpha, \beta)$,

where $S^2_n = (1/n)\sum_i (X_i - \bar{X}_n)^2$ is the biased (1/$n$) sample variance — the same form used by the MLE for the Normal variance, and the form that makes the MoM and the maximum-likelihood approach align as cleanly as possible. With $n = 200$ samples from $\text{Gamma}(3, 2)$, the typical estimate sits within ±0.1 of $\alpha = 3$ on a single replicate. The MLE for $\alpha$ is more accurate by a factor of $1/\sqrt{0.68} \approx 1.21$ in standard error (the ARE = 0.68 of §15.7); whether that matters depends on the application.

For $X \sim \text{Beta}(\alpha, \beta)$, the mean and variance are

A short manipulation yields

with $S^2_n$ the same biased sample variance. The MLE for the Beta has no closed form — the score equations involve digamma functions — so the gap in computational cost is genuine. We will see in §15.8 that fitting a Beta in scipy or R uses the MoM estimate as the warm start for the MLE iteration.

Let $X_1, \ldots, X_n$ be iid from a model with parameter $\theta_0 \in \Theta \subseteq \mathbb{R}^k$. Suppose:

1. The first $k$ population moments $\mu_1(\theta), \ldots, \mu_k(\theta)$ exist and are continuous in $\theta$.
2. The moment map $g(\theta) = (\mu_1(\theta), \ldots, \mu_k(\theta))$ is one-to-one on $\Theta$, and $g^{-1}$ is continuous on a neighborhood of $g(\theta_0)$.

Then $\hat{\theta}_{\text{MoM}} \xrightarrow{a.s.} \theta_0$ as $n \to \infty$.

Under the conditions of Theorem 3, additionally suppose:

3. $\mathbb{E}_{\theta_0}[X^{2k}] < \infty$ (finite higher moments — needed for the CLT to apply to sample moments).
4. The moment map $g$ is continuously differentiable in a neighborhood of $\theta_0$ with non-singular Jacobian $\nabla g(\theta_0)$.

Then

where $\Sigma(\theta_0)$ is the $k \times k$ covariance matrix of $(X, X^2, \ldots, X^k)$ under $\theta_0$.

The four assumptions of Theorem 4 are exactly the ingredients the proof needs. (1)–(2) are the identification assumptions: the moments exist and the moment map distinguishes parameters. (3) is the Lyapunov-style higher-moment assumption that lets the CLT apply componentwise to the iid moments $X_i^j$. (4) is the smoothness assumption that lets the delta method propagate the limit through $g^{-1}$. When any of these fails — most notably (1) for Cauchy and (3) for heavy-tailed Pareto — MoM does not just have a worse asymptotic variance; it has no asymptotic variance, because the underlying CLT does not apply. We see those failure modes in §15.6.

The Exponential$(\lambda)$ MoM is $\hat{\lambda}_{\text{MoM}} = 1/\bar{X}_n$. By Theorem 4 with $k = 1$ and the delta method on $h(m) = 1/m$,

Computing: $h'(m) = -1/m^2$, so $h'(1/\lambda_0) = -\lambda_0^2$; and $\operatorname{Var}(X) = 1/\lambda_0^2$. Hence

The Cramér–Rao lower bound for the Exponential rate (Example 13.10) is $\lambda_0^2/n$ — exactly what we just computed. So the Exponential MoM is asymptotically efficient, achieving the CRLB. This is the §15.7 statement that the MoM matches the MLE for natural exponential families, made concrete.

A distribution has finite $k$-th moment if $\mathbb{E}[|X|^k] < \infty$. When $\mathbb{E}[|X|^k] = \infty$, the $k$-th moment does not exist (in either the abstract or the practical sense), the SLLN does not apply to $\overline{X^k}$, and the MoM equation $\overline{X^k} = \mu_k(\hat{\theta})$ has no meaningful right-hand side. The MoM estimator is then undefined, not merely badly behaved.

The Cauchy distribution with location $\theta$ and scale 1 has density $f(x; \theta) = 1/[\pi(1 + (x - \theta)^2)]$. Its tails are $O(1/x^2)$, so

i.e., $\mathbb{E}[|X|] = \infty$ and the population mean does not exist. The MoM equation $\hat{\theta} = \bar{X}_n$ has a left-hand side that does not converge as $n \to \infty$: the sample mean of Cauchy data exhibits wild excursions at every scale — large jumps appear with positive probability at every $n$, so the sequence $\bar{X}_n$ has no limit (almost surely it does not converge).

The MLE, by contrast, is consistent. The log-likelihood

is twice continuously differentiable, has a unique maximum (the MLE solves a polynomial equation in $\theta$ of degree $2n - 1$ when expanded), and the Fisher information $I(\theta) = 1/2$ is finite. The MLE is asymptotically Normal with variance $2/n$. Cauchy is the cleanest example of MoM failing while MLE succeeds: the population mean does not exist, but the population is fully identified by the location of its peak — and the likelihood sees the peak.

The Pareto distribution with shape $\alpha$ and scale $x_m$ has density $f(x; \alpha, x_m) = \alpha x_m^\alpha / x^{\alpha+1}$ for $x \geq x_m$. Its first moment is

which converges only when $\alpha > 1$, in which case $\mathbb{E}[X] = \alpha x_m / (\alpha - 1)$. For $\alpha \leq 1$, $\mathbb{E}[X] = \infty$, and the MoM moment-matching equation $\bar{X}_n = \alpha x_m / (\alpha - 1)$ has no consistent population side. Empirically, $\bar{X}_n$ for Pareto$(0.8)$ data does not converge but instead grows roughly like $n^{1/\alpha - 1}$ — heavy-tailed extremes pull it up at every sample size.

The MLE for the Pareto $\alpha$ is closed-form: $\hat{\alpha}_{\text{MLE}} = n / \sum_i \log(X_i / x_m)$, and is consistent with finite asymptotic variance, regardless of whether $\alpha > 1$. So once again — even more cleanly than Cauchy — the MoM fails because $\mathbb{E}[X]$ does not exist, while the MLE works because the log-transformed data has finite moments.

Consider the mixture model $X \sim \tfrac{1}{2}\mathcal{N}(0, 1) + \tfrac{1}{2}\mathcal{N}(0, \sigma^2)$ — equal mixing weights, common mean zero, two variance components. The first three population moments are

Only the second moment depends on $\sigma^2$. The first moment is identically zero (by symmetry), and the third moment is identically zero (by symmetry of each Normal component). So the moment map $g(\sigma^2) = (\mu_1, \mu_2, \mu_3) = (0, (1 + \sigma^2)/2, 0)$ is invertible only if we use the second moment — and it has dimension 1 in moment space, embedded in 3-dimensional moment space. With three moments and one parameter, MoM is over-determined; with one moment and one parameter, MoM works trivially.

The deeper issue surfaces when we try to estimate more parameters in the mixture — say, both variances and the mixing weight. With two distinct components having equal means, the moment system collapses: the same set of population moments is consistent with a continuum of parameter combinations. The moment map is not invertible. MoM fails because the model is not identifiable from moments alone. The MLE faces the same issue — but the MLE has access to the full likelihood, which uses the shape of each observation’s contribution and can identify the components even when their moments coincide.

The pattern in Examples 10–12 is general. When MoM fails, it fails visibly: an undefined limit, an estimate outside the parameter space, or a system that simply has no unique solution. This is, perversely, a strength — MoM tells you when it has failed, often more loudly than the MLE does. If $\hat{\sigma}^2_{\text{MoM}} < 0$ or $\hat{p}_{\text{MoM}} > 1$ on real data, the appropriate response is not to project the estimate onto the constraint set but to ask whether the model is right.

For two consistent, asymptotically Normal estimators $\hat{\theta}_1$ and $\hat{\theta}_2$ of the same parameter $\theta$, with asymptotic variances $\sigma^2_1(\theta)/n$ and $\sigma^2_2(\theta)/n$, the asymptotic relative efficiency of $\hat{\theta}_1$ relative to $\hat{\theta}_2$ is

When this ratio is less than 1, $\hat{\theta}_1$ is asymptotically less efficient than $\hat{\theta}_2$. ARE = 1 means the two estimators are asymptotically equivalent. (Definition 13.13 of Topic 13.)

Let $\{f(x; \eta) : \eta \in \mathcal{E}\}$ be a $k$-parameter exponential family in natural parameter form,

with sufficient statistic $T(x) = (T_1(x), \ldots, T_k(x))$ and log-partition function $A(\eta)$. The MLE $\hat{\eta}_{\text{MLE}}$ satisfies

This is exactly the moment-matching equation: $\nabla A(\eta) = \mathbb{E}_\eta[T(X)]$ is the population mean of the sufficient statistic, so the MLE equates the model’s expected sufficient statistic to the sample average — which is the MoM equation when the sufficient statistics are the moments. In particular, for distributions whose sufficient statistic is $T(x) = x$ (Exponential, Bernoulli, Poisson), the MLE is the moment-matching estimator $\hat{\theta}_{\text{MoM}} = \bar{X}_n$, so $\operatorname{ARE}(\text{MoM}, \text{MLE}) = 1$.

For $X_1, \ldots, X_n$ iid $\text{Gamma}(\alpha, \beta)$, with the MoM estimator $\hat{\alpha}_{\text{MoM}} = \bar{X}_n^2/S^2_n$ and the MLE $\hat{\alpha}_{\text{MLE}}$ from the score equation, the asymptotic variances are

and (from the inverse Fisher-information matrix at $(\alpha, \beta)$)

where $\psi'(\alpha)$ is the trigamma function. Therefore

At $\alpha = 3$: the trigamma function evaluates to $\psi'(3) = \pi^2/6 - 1 - 1/4 \approx 0.394934$, so $\alpha \psi'(\alpha) - 1 \approx 0.184802$. Then

The Monte Carlo verification at $n = 200$ over 1000 replications gives an empirical $\operatorname{Var}(\hat{\alpha}_{\text{MoM}}) / \operatorname{Var}(\hat{\alpha}_{\text{MLE}}) \approx 0.69$ — consistent with the theoretical 0.676 to within Monte-Carlo noise. So the MLE extracts about 1/0.68 ≈ 1.47 times as much information from each observation as MoM does — equivalently, MoM needs 47% more data to match MLE precision.

The Exponential family has $T(x) = x$, so by Theorem 5, the MLE $\hat{\lambda}_{\text{MLE}} = 1/\bar{X}_n$ coincides with the MoM. Asymptotic variance for both is $\lambda^2$ (per Example 9), and $\operatorname{ARE}(\hat{\lambda}_{\text{MoM}}, \hat{\lambda}_{\text{MLE}}) = 1$ exactly. This is not approximate equality — the two estimators are the same function of the data. They cannot have different asymptotic variances.

Three reasons, in increasing order of importance.

(1) Closed form. When the moment map inverts in elementary functions — Gamma, Beta, Weibull (numerically), Pareto with $\alpha > 1$ — MoM gives an estimate in two lines of arithmetic. No iteration, no convergence diagnostics, no risk of missing a local maximum.

(2) Warm-start for MLE. The next section makes this precise. When MLE requires Newton-Raphson, the MoM estimate is almost always a good starting point: it is consistent (so it converges to $\theta_0$ at the same rate as the MLE), it is computable, and it lies inside the parameter space whenever the data are not pathological. Newton from MoM converges in 3–5 iterations at typical $n$; Newton from an arbitrary start can take 8–12 or fail to converge at all.

(3) Diagnostic. When MoM and MLE disagree wildly, the model is suspect. The two recipes use different summaries of the data — sample moments for MoM, the full likelihood for MLE — and a large discrepancy means the model is not capturing the moment structure consistently with the likelihood structure. This is a free goodness-of-fit signal, available without additional computation.

The MLE is asymptotically efficient. The MoM is consistent. The two estimates converge to the same true $\theta_0$ at the same $\sqrt{n}$ rate, and the difference between them is $O_p(1/\sqrt{n})$. So the MoM estimate places us, with high probability, within an $O(1/\sqrt{n})$ neighborhood of the MLE — close enough that one or two Newton steps land at the MLE to machine precision. This is the algorithmic reason MoM matters even when we have decided to fit by maximum likelihood.

The conceptual analogue in optimization is multi-resolution: solve a cheap approximate problem to get a good initial point, then refine with the expensive exact procedure. MoM-warm-started MLE is exactly this pattern at the level of statistical estimation.

For a $\text{Gamma}(\alpha = 3, \beta = 2)$ sample of size $n = 100$, the MoM estimate of the shape is typically $\hat{\alpha}_{\text{MoM}} \approx 3.0 \pm 0.3$ on a single replicate. Starting Newton-Raphson on the score equation $\log\beta - \psi(\alpha) + \overline{\log X} = 0$ from this initial value, convergence to machine precision (step size $< 10^{-8}$) requires 4 iterations in the typical case. Starting from $\alpha_0 = 0.5$ takes 9 iterations — Newton has to traverse the steep left tail of the score function. Starting from $\alpha_0 = 2 \times \hat{\alpha}_{\text{MoM}} \approx 6$ takes 10 iterations, with the first step often overshooting. The pattern holds across sample sizes from $n = 50$ to $n = 1000$: MoM-warm-started Newton converges in 3–5 steps; arbitrary-start Newton needs 8–12.

This is not a textbook idea waiting to be implemented. It is the dominant pattern in the actual codebases that fit distributions to data.

• scipy.stats.fit: when a parametric family is specified, scipy first computes a MoM estimate (it calls this “method of moments” internally), then refines via L-BFGS or a Nelder–Mead simplex on the negative log-likelihood. MoM is the explicit default initial point.
• R’s MASS::fitdistr: the same pattern. MoM init, then optim() for the MLE.
• Statsmodels (Python): iteratively reweighted least squares for GLMs, with the OLS estimate (a MoM-like start) as the initial guess.

The dominant alternative — searching for a starting point by random restart or by grid search — is much slower and rarely necessary when MoM is available. The exceptions are models where MoM is undefined (mixtures with non-identifiable components, latent-variable models) or pathological (Cauchy). For everything else, MoM is the warm-start.

Let $X_1, \ldots, X_n$ be iid from a distribution $P$ with parameter $\theta \in \Theta$. Given an estimating function $\psi: \mathcal{X} \times \Theta \to \mathbb{R}^k$ satisfying the moment condition $\mathbb{E}_{\theta_0}[\psi(X; \theta_0)] = 0$ at the true parameter, the M-estimator $\hat{\theta}_n$ is any solution of the estimating equation

Existence and uniqueness require regularity on $\psi$ (e.g., monotonicity in $\theta$, or strict convexity of an associated loss); see van der Vaart Chapter 5 for the formal conditions.

When $\psi(x; \theta) = \nabla_\theta \rho(x; \theta)$ for some scalar loss function $\rho: \mathcal{X} \times \Theta \to \mathbb{R}$, the M-estimator is equivalently the minimizer of the empirical risk:

In this gradient form, the estimating equation $\nabla_\theta \sum_i \rho = 0$ is the first-order condition for the minimization. The label Z-estimator (for “zero” of the estimating equation) is sometimes reserved for the general form even when no underlying loss exists; M-estimator (for “minimum” of an empirical risk) emphasizes the loss-function origin.

For an iid sample from a parametric family $\{f(x; \theta) : \theta \in \Theta\}$, the maximum-likelihood estimator $\hat{\theta}_{\text{MLE}}$ is the M-estimator with estimating function

The moment condition $\mathbb{E}_{\theta_0}[\psi(X; \theta_0)] = 0$ holds because the population score has mean zero (Theorem 13.8). The associated loss is the negative log-density: $\rho(x; \theta) = -\log f(x; \theta)$, with $\hat{\theta}_{\text{MLE}}$ the minimizer of the negative log-likelihood.

For a $k$-parameter model with first $k$ population moments $\mu_j(\theta)$, the method-of-moments estimator $\hat{\theta}_{\text{MoM}}$ is the M-estimator with vector-valued estimating function

The moment condition $\mathbb{E}_{\theta_0}[\psi(X; \theta_0)] = (\mu_1 - \mu_1, \ldots, \mu_k - \mu_k) = \mathbf{0}$ holds at the true parameter by definition. There is no associated loss function in general — MoM is a Z-estimator, not an M-estimator in the strict loss-minimization sense.

For a regression model $Y_i = X_i^\top \beta + \varepsilon_i$ with $\mathbb{E}[\varepsilon | X] = 0$, the ordinary least-squares estimator $\hat{\beta}_{\text{OLS}}$ is the M-estimator with

The moment condition $\mathbb{E}[X(Y - X^\top \beta_0)] = \mathbb{E}[X \varepsilon] = 0$ holds by the regression assumption. The loss is $\rho(x, y; \beta) = \tfrac{1}{2}(y - x^\top \beta)^2$. We develop OLS in detail in Linear Regression; here it suffices to note that OLS, MLE for a Normal-error linear model, falls into the M-estimation framework.

For an M-estimator with smooth estimating function $\psi(x; \theta)$, define

the sensitivity matrix — how strongly the population estimating function responds to changes in $\theta$. And define

the variability matrix — the covariance of $\psi$ at $\theta_0$ (using the moment condition $\mathbb{E}[\psi] = 0$ to drop the mean). When $\theta$ is scalar, both $A$ and $B$ are scalars; for vector $\theta$, they are $k \times k$ matrices.

Sign conventions vary across references. We adopt the convention that $A$ is positive definite in the regular case (so the negative sign appears explicitly). The alternative convention $A = \mathbb{E}[\nabla \psi]$ flips the sign and treats $A$ as negative definite.

Under suitable identifiability and regularity conditions on $\psi$ (e.g., $\theta \mapsto \mathbb{E}_{\theta_0}[\psi(X; \theta)]$ is uniquely zero at $\theta = \theta_0$ and continuous in a neighborhood, plus uniform convergence of $\Psi_n/n$ to its expectation), the M-estimator satisfies

The argument parallels the MLE consistency proof (Theorem 14.3): uniform LLN on $\Psi_n/n$, plus the unique-zero hypothesis, plus a continuity argument.

Under the conditions of Theorem 9 plus continuous differentiability of $\psi$ in $\theta$ near $\theta_0$, finite second moments of $\psi$, and invertibility of $A(\theta_0)$,

where the asymptotic variance is the sandwich form

For scalar $\theta$ this reduces to $V = B/A^2$.