Covariance of Multiple Functions of Estimated Parameters

Why Joint Uncertainty Matters

A fitted model often supplies more than one quantity of interest. An actuary might report two functions of the same estimated severity parameters, such as a product and a ratio, or several transformed parameters used in a later calculation. These estimates are generally dependent because they are built from the same data. Calculating a standard error for each function separately does not describe that dependence; the covariance between the functions is also needed.

The multivariate delta method provides a first-order approximation to the full covariance matrix. Its central calculation is compact: form a Jacobian from the gradients of the functions, then multiply the parameter covariance matrix on the left and right by that Jacobian and its transpose. The details that matter most are using the correct gradient for each function, evaluating it at the estimated parameter vector, and retaining every covariance term.

Definition: Let \(\boldsymbol{\theta}=(\theta_1,\ldots,\theta_p)^T\) be a parameter vector and let \(\hat{\boldsymbol{\theta}}\) estimate it. Suppose \(\hat{\boldsymbol{\theta}}\) is approximately unbiased with covariance matrix \(V\), and let \(g_1,\ldots,g_q\) be differentiable functions of the parameters. The Jacobian \(J\) is the \(q\)-by-\(p\) matrix whose \(j\)th row is the gradient of \(g_j\). The multivariate delta-method approximation is

\operatorname{Cov}\!\left( \begin{pmatrix} g_1(\hat{\boldsymbol{\theta}})\\ \vdots\\ g_q(\hat{\boldsymbol{\theta}}) \end{pmatrix} \right) \approx J V J^T.

The gradients and \(J\) are evaluated at the parameter value, usually the estimate \(\hat{\boldsymbol{\theta}}\). If \(V\) is the covariance of the estimator, use it directly; if the available matrix is an asymptotic covariance per observation, divide by the sample size first.

Why the Jacobian Formula Holds

The approximation follows from a first-order Taylor expansion. Write \(g(\boldsymbol{\theta})\) for the vector containing all \(q\) functions. When the functions are differentiable near the true parameter and the estimator is consistent and approximately normal, the expansion around the true parameter is

g(\hat{\boldsymbol{\theta}}) = g(\boldsymbol{\theta}) + J(\boldsymbol{\theta})(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}) + \boldsymbol{r},

where the remainder \(\boldsymbol{r}\) is of smaller order than the first-order estimation error under the usual large-sample conditions. The constant vector \(g(\boldsymbol{\theta})\) has no effect on covariance. For the linear term, covariance transforms according to

\operatorname{Cov}\!\left(J(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta})\right) = J\operatorname{Cov}(\hat{\boldsymbol{\theta}})J^T = JVJ^T.

Ignoring the smaller-order remainder gives the delta-method result. If \(g\) is linear, the remainder is zero and the covariance transformation is exact. For nonlinear functions it is an approximation, usually most useful when the estimator is concentrated near the true parameter and the gradients do not change sharply over the range of estimation uncertainty.

In a regular model, if \(\hat{\boldsymbol{\theta}}\) is approximately normal with covariance \(\Sigma/n\), then \(V\) in the formula is \(\Sigma/n\). In practice, one commonly substitutes an estimated covariance matrix, such as one obtained from an information matrix, and evaluates the Jacobian at \(\hat{\boldsymbol{\theta}}\). The formula does not repair a poor normal approximation, a non-differentiable function at the estimate, or a parameter estimate near a boundary where local linearization may be unreliable.

List the target functions.
Write each desired quantity as \(g_j(\boldsymbol{\theta})\), keeping their order fixed.

Differentiate each function.
Find its gradient with respect to every component of \(\boldsymbol{\theta}\), then evaluate at \(\hat{\boldsymbol{\theta}}\).

Transform the covariance.
Compute \(JVJ^T\), preserving the off-diagonal entries of \(V\) and all cross terms.

Interpret the result.
The diagonal entries are approximate variances; their square roots are standard errors. The off-diagonal entries are approximate covariances between the transformed estimates.

Worked Examples

Worked Example: Two Linear Functions of Gamma Parameters

A hypothetical actuarial analyst fits a Gamma severity model parameterized by shape \(\alpha\) and scale \(\beta\) to independent loss observations. Suppose the fitted values are \(\hat{\alpha}=2\) and \(\hat{\beta}=1\), and the estimated covariance matrix of \((\hat{\alpha},\hat{\beta})^T\) is

V= \begin{pmatrix} 0.14&0.03\\ 0.03&0.21 \end{pmatrix}.

Find the covariance matrix and standard errors for \(g_1(\alpha,\beta)=\alpha+2\beta\) and \(g_2(\alpha,\beta)=3\alpha+\beta\). The parameter covariance matrix is assumed to be a valid estimated covariance matrix from the fitted model; because these functions are linear, the covariance transformation is exact conditional on that matrix.

Method and Jacobian. The gradients are \(\nabla g_1=(1,2)\) and \(\nabla g_2=(3,1)\), so

J= \begin{pmatrix} 1&2\\ 3&1 \end{pmatrix}.

First variance. The variance of \(g_1\) is

\begin{aligned} \operatorname{Var}(g_1) &\approx (1,2)V(1,2)^T\\ &=0.14+4(0.21)+2(1)(2)(0.03)\\ &=0.14+0.84+0.12=1.10. \end{aligned}

Second variance. Similarly,

\begin{aligned} \operatorname{Var}(g_2) &\approx (3,1)V(3,1)^T\\ &=9(0.14)+0.21+2(3)(1)(0.03)\\ &=1.26+0.21+0.18=1.65. \end{aligned}

Covariance. The cross-covariance is

\begin{aligned} \operatorname{Cov}(g_1,g_2) &\approx (1,2)V(3,1)^T\\ &=1(3)(0.14)+2(1)(0.21) +(1(1)+2(3))(0.03)\\ &=0.42+0.42+0.21=1.05. \end{aligned}

Thus the transformed covariance matrix is \(\begin{pmatrix}1.10&1.05\\1.05&1.65\end{pmatrix}\). The standard errors are \(\sqrt{1.10}=1.0488\) and \(\sqrt{1.65}=1.2845\), rounded to four decimal places. As a check, the determinant is \(1.10(1.65)-1.05^2=0.7125>0\), consistent with a positive-definite covariance matrix. The common mistake avoided here is dropping the original parameter covariance \(0.03\): it contributes to both variances and to the covariance of the functions.

Worked Example: Product and Ratio of Gamma Parameters

For a second hypothetical independent-loss sample, a Gamma severity model uses shape \(\alpha\) and scale \(\beta\). Suppose \((\hat{\alpha},\hat{\beta})=(2,3)\) and the estimated covariance matrix is

V= \begin{pmatrix} 0.04&0.01\\ 0.01&0.09 \end{pmatrix}.

Estimate the covariance matrix for \(g_1(\alpha,\beta)=\alpha\beta\) and \(g_2(\alpha,\beta)=\alpha/\beta\). This example uses the first-order delta method because both functions are nonlinear.

Gradients at the fitted values. The derivatives are \(\nabla g_1=(\beta,\alpha)\) and \(\nabla g_2=(1/\beta,-\alpha/\beta^2)\). At \((2,3)\),

J= \begin{pmatrix} 3&2\\ 1/3&-2/9 \end{pmatrix}.

Variance of the product.

\begin{aligned} \operatorname{Var}(g_1) &\approx 3^2(0.04)+2^2(0.09)+2(3)(2)(0.01)\\ &=0.36+0.36+0.12=0.84. \end{aligned}

Variance of the ratio.

\begin{aligned} \operatorname{Var}(g_2) &\approx \left(\frac13\right)^2(0.04) +\left(-\frac29\right)^2(0.09) +2\left(\frac13\right)\left(-\frac29\right)(0.01)\\ &=0.0044444+0.0044444-0.0014815\\ &=0.0074074=\frac{1}{135}. \end{aligned}

Covariance. First calculate

V \begin{pmatrix}1/3\\-2/9\end{pmatrix} = \begin{pmatrix}1/90\\-1/60\end{pmatrix}.

Therefore,

\operatorname{Cov}(g_1,g_2) \approx (3,2) \begin{pmatrix}1/90\\-1/60\end{pmatrix} =\frac{3}{90}-\frac{2}{60}=0.

The approximate covariance matrix is \(\begin{pmatrix}0.84&0\\0&1/135\end{pmatrix}\); the standard errors are \(0.9165\) and \(0.0861\), rounded to four decimal places. The zero off-diagonal entry means the two transformed estimates have zero covariance to first order under this calculation. It does not establish that the nonlinear estimates are exactly independent. The calculation is checked by retaining the negative cross term in the ratio variance; omitting it would give a different variance.

Worked Example: Two Functions of a Shifted-Lognormal Model

Suppose independent hypothetical loss observations are fitted with a shifted-lognormal severity model \(X=\tau+\exp(Z)\), where \(Z\sim N(\mu,\sigma^2)\), \(\sigma>0\), and \(\tau\) is the shift. Let \((\hat{\mu},\hat{\sigma},\hat{\tau})=(1,2,1)\), with estimated covariance matrix

V= \begin{pmatrix} 0.04&0.01&0\\ 0.01&0.09&0.02\\ 0&0.02&0.16 \end{pmatrix}.

For \(g_1(\mu,\sigma,\tau)=\mu\sigma+\tau\) and \(g_2(\mu,\sigma,\tau)=\mu+\tau/\sigma\), calculate their approximate covariance matrix. Assume the fitted model is regular and the displayed matrix is the covariance estimate for the parameter estimates.

The gradients are \((\sigma,\mu,1)\) and \((1,-\tau/\sigma^2,1/\sigma)\). At the fitted values, they are \((2,1,1)\) and \((1,-1/4,1/2)\), respectively. Thus \(J=\begin{pmatrix}2&1&1\\1&-1/4&1/2\end{pmatrix}\).

First variance. The diagonal contributions are \(2^2(0.04)+1^2(0.09)+1^2(0.16)=0.41\). The off-diagonal contributions are \(2(2)(1)(0.01)=0.04\) and \(2(1)(1)(0.02)=0.04\). Hence \(\operatorname{Var}(g_1)\approx0.41+0.04+0.04=0.49\).

Second variance. Its diagonal contributions are \(1^2(0.04)+(-1/4)^2(0.09)+(1/2)^2(0.16)=0.04+0.005625+0.04=0.085625\). The \((\mu,\sigma)\) and \((\sigma,\tau)\) contributions are \(2(1)(-1/4)(0.01)=-0.005\) and \(2(-1/4)(1/2)(0.02)=-0.005\); the \((\mu,\tau)\) contribution is zero. Therefore \(\operatorname{Var}(g_2)\approx0.075625\).

Cross-covariance. Multiplying the parameter covariance matrix by the second gradient gives

V \begin{pmatrix}1\\-1/4\\1/2\end{pmatrix} = \begin{pmatrix}0.0375\\-0.0025\\0.075\end{pmatrix}.

Taking the inner product with the first gradient yields \(2(0.0375)-0.0025+0.075=0.1475\). The result is \(\begin{pmatrix}0.49&0.1475\\0.1475&0.075625\end{pmatrix}\), with standard errors \(0.7000\) and \(0.2750\). The positive determinant, \(0.49(0.075625)-0.1475^2=0.0153\), provides a numerical check. This example also shows why a three-parameter problem requires accounting for each parameter pair, not just the three diagonal variances.

Reading and Checking the Result

For \(q\) functions, \(JVJ^T\) is a \(q\)-by-\(q\) matrix even when the original parameter vector has a different dimension. The diagonal entries are variances on the scales of the target functions. Taking their square roots gives standard errors on those scales. An off-diagonal entry describes how two estimated functions vary together: a positive value indicates that their first-order estimation errors tend to move in the same direction, while a negative value indicates opposite movement.

The resulting matrix should be symmetric and positive semidefinite, apart from tiny numerical discrepancies caused by rounding. For two functions, a useful check is that both diagonal entries are nonnegative and the squared covariance does not exceed the product of the variances:

\operatorname{Cov}(g_1,g_2)^2 \leq \operatorname{Var}(g_1)\operatorname{Var}(g_2).

This follows because every covariance matrix is positive semidefinite. A violation is a strong sign of an arithmetic error, inconsistent rounding, or an invalid input covariance matrix. It is also important to distinguish covariance from correlation: if a standardized measure is needed, divide the covariance by the product of the two standard errors.

Common Exam Trap

The most frequent error is to treat parameter estimates as independent when the supplied covariance matrix has nonzero off-diagonal entries. For a scalar function \(g(\theta_1,\theta_2)\), the correct first-order variance is

\operatorname{Var}(g) \approx \left(\frac{\partial g}{\partial\theta_1}\right)^2V_{11} + \left(\frac{\partial g}{\partial\theta_2}\right)^2V_{22} + 2\frac{\partial g}{\partial\theta_1} \frac{\partial g}{\partial\theta_2}V_{12}.

For example, if the gradient is \((1,2)\) and \(V_{12}=0.03\), the omitted term is \(2(1)(2)(0.03)=0.12\). Leaving it out changes the variance, even though both diagonal terms have been calculated correctly. For covariance between two different functions, the corresponding cross terms must also be included; the full matrix multiplication is often the safest way to avoid losing one.

The method is aligned with the SOA Exam ASTAM (Advanced Short-Term Actuarial Mathematics) syllabus. In an exam calculation, state the parameter ordering, write the Jacobian in that order, and make clear whether the supplied covariance is for the estimator itself or must first be divided by sample size. Those checks prevent errors that can otherwise pass unnoticed in a plausible-looking standard error.

Check Your Understanding

For each problem, use the multivariate delta method and preserve any nonzero parameter covariances.

A hypothetical two-parameter Gamma severity model has fitted shape and scale estimates \((2,1)\), with estimator covariance matrix \(\begin{pmatrix}0.09&0.02\\0.02&0.04\end{pmatrix}\). For \(g_1=\alpha+\beta\) and \(g_2=2\alpha-\beta\), find the approximate covariance matrix. (Answer: \(\begin{pmatrix}0.17&0.16\\0.16&0.32\end{pmatrix}\))
For another Gamma fit, \((\hat{\alpha},\hat{\beta})=(2,4)\) and \(V=\begin{pmatrix}0.04&0.01\\0.01&0.09\end{pmatrix}\). Find the approximate covariance matrix for \(g_1=\alpha+\beta\) and \(g_2=\alpha\beta\). (Answer: \(\begin{pmatrix}0.15&0.40\\0.40&1.16\end{pmatrix}\))
For Gamma shape and scale estimates \((\hat{\alpha},\hat{\beta})=(2,4)\) with \(V=\begin{pmatrix}0.04&0.01\\0.01&0.09\end{pmatrix}\), find the delta-method covariance matrix for \(g_1=\alpha/\beta\) and \(g_2=\beta/\alpha\). (Answer: \(\begin{pmatrix}0.00328125&-0.013125\\-0.013125&0.0525\end{pmatrix}\))
A hypothetical three-parameter severity model has estimate \((\hat{a},\hat{b},\hat{c})=(2,3,1)\) and covariance matrix \(\operatorname{diag}(0.01,0.04,0.09)\). For \(g_1=ab\) and \(g_2=a+c\), find the approximate covariance matrix. (Answer: \(\begin{pmatrix}0.25&0.03\\0.03&0.10\end{pmatrix}\))

What You'll Learn