Why a Covariance Matrix Is Needed
When a method-of-moments procedure estimates several parameters, reporting one standard error at a time leaves out important information: the estimates may vary together. For example, increasing one fitted parameter may require decreasing another to keep the fitted moments near the observed moments. A covariance matrix records both the variability of each estimate and the dependence between estimates.
The calculation has two parts. First, determine how variable the sample moments are. Second, account for how changes in the parameters change the population moments. The Jacobian of the population moment equations connects these two sources of variation.
This formula concerns a just-identified system: the number of moment equations equals the number of parameters. The usual result requires independent, identically distributed observations, finite second moments for the components of \(h(X)\), a correctly specified and locally identifiable model, and a nonsingular Jacobian at the true parameter. The estimator must also be consistent and lie in a region where the moment equations can be differentiated. If these conditions fail—for example, if \(G\) is singular—the inverse-Jacobian formula does not apply as stated.
The Covariance Formula and Its Justification
Under those conditions, the central limit theorem gives an approximate multivariate normal distribution for the sample moments. A first-order expansion of the moment equations then transfers that variation to the parameter estimates.
Thus, for a large sample, the estimated parameter vector is approximately multivariate normal with covariance matrix \[ \operatorname{Cov}(\hat{\theta})\ \approx\ \frac{1}{n}G(\theta)^{-1}\Sigma(\theta)\{G(\theta)^{-1}\}^{\mathsf{T}}. \] The covariance on the right-hand side of the limiting distribution is for the scaled quantity \(\sqrt n(\hat{\theta}-\theta)\). The covariance of \(\hat{\theta}\) itself includes the factor \(1/n\).
Set \(h(X)\) to the functions used in the estimating equations, then calculate \(g(\theta)=E_\theta[h(X)]\) and \(\Sigma=\operatorname{Cov}_\theta(h(X))\).
Form \(G=\partial g/\partial\theta^\mathsf{T}\), keeping the parameter order consistent with the moment order.
Calculate \(G^{-1}\Sigma(G^{-1})^\mathsf{T}\), then divide by \(n\) for the approximate covariance matrix of \(\hat{\theta}\). Substitute parameter estimates if an estimated covariance matrix is required.
To see why the formula holds, write the moment equations as \(g(\hat{\theta})=\bar{h}\), while \(g(\theta)=E[h(X)]\). A first-order Taylor expansion around the true \(\theta\) gives
The sample mean \(\bar{h}\) has approximate covariance \(\Sigma/n\), because it averages \(n\) independent moment vectors. Multiplying a random vector on the left by a matrix \(A\) changes its covariance from \(V\) to \(AVA^\mathsf{T}\). Taking \(A=G^{-1}\) therefore gives \(G^{-1}(\Sigma/n)(G^{-1})^\mathsf{T}\), which is the stated approximation. The Taylor expansion also shows why the Jacobian must be evaluated at the true parameter for the theoretical result, or at a consistent estimate for a practical covariance estimate.
Worked Calculations
Worked Example: Two Raw Moments for a Normal Model
Suppose \(X_1,\ldots,X_n\) are independent observations from \(N(\mu,\sigma^2)\), with unknown mean \(\mu\) and variance \(\sigma^2\). Use the first two raw moments to estimate \(\theta=(\mu,\sigma^2)^\mathsf{T}\). Find the asymptotic covariance matrix of the estimates at \(\mu=2\), \(\sigma^2=3\), and \(n=100\).
Method. Choose \(h(X)=(X,X^2)^\mathsf{T}\). The population equations are \(g_1(\theta)=E[X]=\mu\) and \(g_2(\theta)=E[X^2]=\mu^2+\sigma^2\). The sample equations are \(\bar X=\hat\mu\) and \(n^{-1}\sum X_j^2=\hat\mu^2+\widehat{\sigma^2}\).
The Jacobian, with parameter order \((\mu,\sigma^2)\), is
For a normal random variable, \(E[X^3]=\mu^3+3\mu\sigma^2\) and \(E[X^4]=\mu^4+6\mu^2\sigma^2+3\sigma^4\). Therefore \[ \operatorname{Var}(X)=\sigma^2,\quad \operatorname{Cov}(X,X^2)=E[X^3]-E[X]E[X^2]=2\mu\sigma^2, \] and \[ \operatorname{Var}(X^2)=E[X^4]-\{E[X^2]\}^2=4\mu^2\sigma^2+2\sigma^4. \] At \(\mu=2\), \(\sigma^2=3\), these give \(\Sigma=\begin{pmatrix}3&12\\12&66\end{pmatrix}\): the lower-right entry is \(4(2^2)(3)+2(3^2)=48+18=66\).
Here \(G^{-1}=\begin{pmatrix}1&0\\-4&1\end{pmatrix}\). Transforming the moment covariance gives \[ G^{-1}\Sigma(G^{-1})^\mathsf{T} = \begin{pmatrix}3&0\\0&18\end{pmatrix}. \] For instance, the off-diagonal entry is \(12-4(3)=0\), and the lower-right entry is \(16(3)-8(12)+66=18\). Dividing by \(n=100\), the estimated parameters have asymptotic covariance matrix \[ \operatorname{Cov}(\hat\mu,\widehat{\sigma^2})\ \approx\ \begin{pmatrix}0.03&0\\0&0.18\end{pmatrix}. \] The corresponding asymptotic standard errors are \(\sqrt{0.03}=0.1732\) and \(\sqrt{0.18}=0.4243\), to four decimal places. The zero covariance here is a result of the normal model and these moment choices; it is not a general property of method-of-moments estimates.
Check. The variance estimate based on the second raw moment has asymptotic variance \(2\sigma^4/n=18/100=0.18\), matching the matrix calculation. A common error is to use \(E[X^2]\) as the variance of \(X^2\); variance requires subtracting \(\{E[X^2]\}^2\).
Worked Example: Two Category Probabilities
In a hypothetical portfolio, each claim is classified as either type A, type B, or neither. Assume independent claims with probabilities \(a=0.20\), \(b=0.30\), and \(1-a-b=0.50\), respectively. Estimate \(a\) and \(b\) from \(n=200\) claims and calculate the asymptotic covariance matrix of the two estimates.
Method. For each claim, let \(h(X)=(I_A,I_B)^\mathsf{T}\), where each indicator is 1 if the claim has that type and 0 otherwise. Then \(g(a,b)=(a,b)^\mathsf{T}\), so \(G\) is the identity matrix and the moment estimates are the observed proportions.
An indicator with success probability \(a\) has variance \(a(1-a)=0.20(0.80)=0.16\). The two indicators cannot both equal 1, so \(E[I_A I_B]=0\), giving \(\operatorname{Cov}(I_A,I_B)=0-ab=-0.06\). Similarly, \(\operatorname{Var}(I_B)=0.30(0.70)=0.21\). Hence
The standard errors are \(\sqrt{0.0008}=0.0283\) for \(\hat a\) and \(\sqrt{0.00105}=0.0324\) for \(\hat b\), to four decimal places. The negative covariance is sensible: for a fixed sample size, observing more type A claims leaves fewer opportunities to observe type B claims.
Check. The probabilities sum to one, and the off-diagonal covariance is negative because the categories are mutually exclusive. A common mistake is to treat the category indicators as independent and set their covariance to zero; mutual exclusivity makes that incorrect.
Worked Example: Paired Exponential Measurements
Suppose each of \(n=100\) hypothetical insured units produces two independent measurements, \(X_1\) and \(X_2\). Assume \(X_1\) is exponential with rate \(\lambda_1=2\), \(X_2\) is exponential with rate \(\lambda_2=0.5\), and the two measurements are independent within each unit. Estimate the two rates by matching the sample means to their population means, and find the asymptotic covariance matrix.
Method. Use \(h(X)=(X_1,X_2)^\mathsf{T}\). For an exponential variable with rate \(\lambda\), \(E[X]=1/\lambda\) and \(\operatorname{Var}(X)=1/\lambda^2\). Thus \(g(\lambda_1,\lambda_2)=(1/\lambda_1,1/\lambda_2)^\mathsf{T}\), and solving the moment equations gives \(\hat\lambda_j=1/\bar X_j\) for \(j=1,2\).
The Jacobian and moment covariance matrix are \[ G= \begin{pmatrix} -1/\lambda_1^2&0\\ 0&-1/\lambda_2^2 \end{pmatrix}, \qquad \Sigma= \begin{pmatrix} 1/\lambda_1^2&0\\ 0&1/\lambda_2^2 \end{pmatrix}. \] The zero covariance follows from the assumed independence of the two measurements. Since \(G^{-1}=\operatorname{diag}(-\lambda_1^2,-\lambda_2^2)\), the transformed covariance is \(\operatorname{diag}(\lambda_1^2,\lambda_2^2)\). Substituting the given rates and dividing by 100 gives
The asymptotic standard errors are \(0.20\) and \(0.05\). The reciprocal moment estimator is nonlinear, but the Jacobian accounts for how a change in a sample mean affects the fitted rate.
Check. A larger rate corresponds to a smaller mean. The negative derivatives in \(G\) are therefore expected, while the resulting variances are positive. A common mistake is to omit the reciprocal derivative and use the variance of the sample mean as the variance of the rate estimate.
Reading and Using the Result
The diagonal entries of the covariance matrix estimate the variances of the parameter estimates. The off-diagonal entries estimate their covariances. To obtain a correlation, divide a covariance by the product of the corresponding standard errors. A negative covariance does not imply that either variance is negative; a valid covariance matrix is symmetric and positive semidefinite.
The matrix order matters. If \(G\) is defined as the derivative of the population moment vector with respect to the parameter vector, the transformation is \(G^{-1}\Sigma(G^{-1})^\mathsf{T}\). Reversing the derivative convention changes how the formula must be written. It is safest to state the dimensions and parameter ordering before multiplying matrices.
In applications, \(\theta\) is unknown, so the theoretical covariance matrix is usually estimated by substituting \(\hat\theta\) into \(G(\theta)\) and \(\Sigma(\theta)\). This plug-in calculation is meaningful only when the fitted parameters satisfy the model's conditions and the Jacobian is nonsingular. A nearly singular Jacobian can produce large estimated variances: the chosen moments then change only slightly for some substantial parameter changes, so the data provide weak information about that parameter combination.
Check Your Understanding
For each problem, identify the moment vector, account for the Jacobian where needed, and remember that covariance of the estimates includes division by \(n\).
- For independent observations from \(N(\mu,\sigma^2)\), use the first two raw moments at \(\mu=1\), \(\sigma^2=2\), and \(n=50\). What is the asymptotic covariance matrix of \((\hat\mu,\widehat{\sigma^2})^\mathsf{T}\)? (Answer: \(\begin{pmatrix}0.04&0\\0&0.16\end{pmatrix}\).)
- In a three-category model, the first two category probabilities are \(0.25\) and \(0.15\), and the remaining probability is \(0.60\). For \(n=100\), what is the covariance matrix of the two observed proportions? (Answer: \(\begin{pmatrix}0.001875&-0.000375\\-0.000375&0.001275\end{pmatrix}\).)
- For paired independent exponential measurements with rates \(\lambda_1=1\) and \(\lambda_2=4\), estimate the rates by reciprocals of sample means. For \(n=80\), give the asymptotic covariance matrix. (Answer: \(\begin{pmatrix}0.0125&0\\0&0.2\end{pmatrix}\).)
- Suppose \(g(\theta_1,\theta_2)=(\theta_1+\theta_2,\theta_1-\theta_2)^\mathsf{T}\), with sample moment covariance \(\Sigma=\begin{pmatrix}4&1\\1&9\end{pmatrix}\) and \(n=100\). Find the asymptotic covariance matrix of the two moment estimates. (Answer: \(\begin{pmatrix}0.0375&-0.0125\\-0.0125&0.0275\end{pmatrix}\).)
- In a two-parameter moment system, \(G\) is singular at the true parameter. Can the inverse-Jacobian covariance formula be applied directly, and what condition in the formula fails? (Answer: No; \(G^{-1}\) does not exist, so the stated just-identified covariance formula cannot be used directly.)