# Standard errors for covariance estimate in R

This is a very simple question: how does one get the standard error for the covariance estimate in R? I estimate the covariance using the cov function but there seems to be no place for it to return a standard error on the estimate. I would prefer the derivation so I can implement myself.

• If you're looking for some black-box function to do this for you, please let us know so we can migrate your question to SO. Otherwise, if you're looking for a formula you can implement, this is the right place to post your question. Could you clarify? – whuber Jan 23 '13 at 18:26
• @whuber: i would be happy with either. i prefer the formula that's why i posted here because i can just implement that myself. i am sure there must be somethign already built in to R to do this – Alex Jan 23 '13 at 18:32

In response to whuber's follow-up, I would advocate that an all-purpose black-box approach would be using a non-parametric bootstrap.

The basic pseudocode is:

1. Jointly resample from observed rows of data, allowing for replications and holding the sample size fixed.
2. Re-estimate covariance in the resampled data.
3. Repeat 1-2 for a sufficient number of iterations.
4. Use the simulated values to compute variance estimates or empirical 0.025 and 0.975 quantiles to form confidence intervals.

An example here:

set.seed(1)

x <- seq(-3, 3, length.out=100)

do.one <- function(x) {
y <- rnorm(100, x)
d <- data.frame(x, y)

## bootstrap out
bs.out <- replicate(1000, {
dd <- d[sample(1:100, replace=TRUE), ]
cov(dd)[1, 2]
})

bs.lower <- quantile(bs.out, 0.025)
bs.upper <- quantile(bs.out, 0.975)

## in the absence of random error, y=x so cov(x, y)=var(x)
(bs.lower < var(x)) & (bs.upper > var(x))
}

o <- replicate(1000, do.one(x))
mean(o) ## should be 95% if bs estimates correct CIs


Feel free to try this simulation with any random or non-random distribution of $X$ and functional form of the mean model. I am unsure (though cautiously optimistic) CIs based on bootstrapped covariance estimates give correct 95% coverage.

• thak you very much. this makes a lot of sense to me. would you be able to derive a closed form solution assuming the variables are normal for example? – Alex Jan 23 '13 at 20:58
• That request is a bit irrelevant to this response. The parametric and non-parametric approaches to variance-component estimation are on polar opposites of the spectrum. Another caveat: unlike the t-test, inference on variance components isn't robust to distributional assumptions, so data which are non-normal will have standard errors that are very miscalibrated using the exact normal sample model – AdamO Jan 23 '13 at 22:19

This is not an answer to the original question, but to your request to AdamO. (As far as I'm concerned he's covered the original question.)

I'd make it a comment but I think it's too long

Would you be able to derive a closed form solution assuming the variables are normal for example?

see

http://en.wikipedia.org/wiki/Estimation_of_covariance_matrices#Concluding_steps

and

http://en.wikipedia.org/wiki/Wishart_distribution

The second link gives the variance of the $(i,j)\,$ element of the distribution of the scatter matrix for multivariate normal random variables. From there you can get the variance of the sample covariance and hence the standard error.

Specifically, $\sum _{{i=1}}^{n}(X_{i}-\overline {X})(X_{i}-\overline {X})^{{\mathrm {T}}}\sim W_{p}(\Sigma ,n-1)$ implies

$\text{Var}(\sum _{{i=1}}^{n}(X_{i}-\overline {X})(Y_{i}-\overline {Y}))=(n-1)(\Sigma_{XY}^2+\Sigma_{XX}\Sigma_{YY})$, or

$\text{Var}(\frac{1}{n-1}\sum _{{i=1}}^{n}(X_{i}-\overline {X})(Y_{i}-\overline {Y}))=(n-1)^{-1}(\Sigma_{XY}^2+\Sigma_{XX}\Sigma_{YY})$

Or, for a more general result,

If $S_{XY}=\frac{1}{n}\sum _{{i=1}}^{n}(X_{i}-\overline {X})(Y_{i}-\overline {Y}))$ then these notes by Thomas S. Richardson, here give

$\text{Var}(S_{XY})=\frac{(n−1)^2}{n^3}(μ_{22}−μ_{11}^2)+ \frac{(n−1)}{n^3} (μ_{11}^2 + μ_{20} μ_{02})$

(where $\mu_{rs}=E[(X-\mu_{_X})^r\,(Y-\mu_{_Y})^s]$)

however, wolfies notes in his answer here that this is incorrect. If I haven't made an error, his result corresponds to a flip of sign on the second $\mu_{11}$ term:

$\text{Var}(S_{XY})=\frac{(n−1)^2}{n^3}(μ_{22}-μ_{11}^2)+ \frac{(n−1)}{n^3} ( μ_{20} μ_{02}-μ_{11}^2 )$

Note that correcting this for the $\frac{1}{n-1}$ version is a simple matter of multiplying the above result by $(\frac{n}{n-1})^2$.

IIRC, there's more details in vol. I of Kendall and Stuart

• Last link is dead, here is a copy from wayback machine web.archive.org/web/20130226032010/http://… – mgilbert Oct 28 '15 at 18:37
• Thanks @mgilbert - I've updated the link, pulled in the equation and given credit (since I now quote it). – Glen_b -Reinstate Monica Oct 28 '15 at 21:41
• I'm assuming the n's in the formula refer to the sample size/number of observations, but could someone help with clarifying what the mu's (with subscript 22, 11, 20, and 02) signify? Sorry for being such a dummy :P. – Jonna Oct 20 '16 at 16:14
• No problem asking; I probably should have defined them. Those are called (central) mixed moments. If you have two subscripts it means you're dealing with two variables ($X,Y$ say); then $\mu_{rs}=E[(X-\mu_{_X})^r\,(Y-\mu_{_Y})^s]$. So $\mu_{11}$ is covariance, $\mu_{20}$ is the variance of the first variable, and so on. They extend in the natural way to more than two variables. – Glen_b -Reinstate Monica Oct 20 '16 at 23:25
• Thanks a lot for clarifying @Glen_b! Just to check, μ22 should then be the product of the variances of X and Y? Here's the code I wrote: V <- ((n-1)^2/n^3)*(Xvar*Yvar - covar^2) + ((n-1)/n^3)*(covar^2 - Xvar*Yvar) – Jonna Oct 21 '16 at 11:39

The OP's question does not define what formula R uses as sample covariance estimator. However, following the link in Glen's answer, I assume that R is using:

$$m_{11} = \frac{1}{n} \sum _{i=1}^n \left(X_i-\bar{X}\right) \left(Y_i-\bar{Y}\right)$$

... also known as the $m_{11}$ sample central moment, which can be expressed in power sum notation $s_{r,t}=\sum _{i=1}^n X_i^r Y_i^t$ (using mathStatica here) as :

... which is a familiar alternative notation.

We seek the variance of the estimator i.e. $\text{Var}(m_{1,1})$. Since the variance operator denotes the $2^\text{nd}$ central moment of $m_{1,1}$, we can find the exact symbolic solution (for any distribution whose moments exists) with the mathStatica function:

I would note that the solution so obtained is different to that referenced in the link given by Glen above to a paper. Perhaps they are computing something else?! There is now a long list of published 'moment of moments' papers that have been shown to contain incorrect results by mathStatica, including some results by Fisher himself, and some of the results in Stuart and Ord - see for instance Spot The Error

There is an alternative defn of sample covariance using $\frac{1}{n-1}$ but that is not the one used in the paper referenced by Glen either.

• taking your output as correct, there's a simple sign error in Richardson's version which I will note in my answer. Thanks. – Glen_b -Reinstate Monica Oct 23 '16 at 22:42
• Actually, R's cov function seems to use $\frac{1}{n-1}$; I don't recall whether I even checked that detail before. – Glen_b -Reinstate Monica Oct 23 '16 at 22:51

There is some confusion in the discussion above.

Simple algebra shows that the expression ascribed to Richardson:

Var$(S_{XY})=\frac{(n−1)^2}{n^3}(\mu_{22}−\mu^2_{11})+\frac{(n−1)}{n^3}(\mu^2_{11}+\mu_{20}\mu_{02})$

is identical to that obtained by wolfies using MathStatica.

Both expressions clearly agree on the coefficients of $\mu_{22}$ and $\mu_{20}\mu_{02}$.

For $\mu_{11}^2$, collecting terms in Richardson's expression gives:

$[-(n-1)^2 + (n-1)]/n^3$

$= [(1-n)(n-1) + (n-1)]/n^3$

$= [(1-n) + 1](n-1)/n^3$

$= (2-n)(n-1)/n^3$

$= -(n-2)(n-1)/n^3$

$= -(-2+n)(-1+n)/n^3$

which is the coefficient for $\mu_{11}^2$ obtained by MathStatica.

[The expression provided by Glen as a "correction" to MathStatica's is not equivalent, as can be seen by substituting $n=2$ and comparing coefficients for $\mu_{11}^2$.]

The correct expression may also be derived in a few steps from results in Kendall's Advanced Theory of Statistics, Kendall & Stuart (1987), Fifth Edition, p.441, Example 13.3, where it is stated that:

Var$(k_{11}) = \frac{1}{n}\kappa_{22} + \frac{1}{n-1}\kappa_{20}\kappa_{02} +\frac{1}{n-1}\kappa^2_{11}$.

Simple algebra shows this is equivalent to the above expressions, noting that $k_{11}$ in K+S is the $k$-statistic, which is the unbiased estimator of the population covariance (aka the (1,1) product cumulant $\kappa_{11}$), so $k_{11} = \frac{n}{n-1} S_{XY}$. It is also necessary to note the following relations between cumulants and moments: $\mu_{11}=\kappa_{11}$, $\mu_{20}=\kappa_{20}$, $\mu_{02}=\kappa_{02}$ and

$\mu_{22} = \kappa_{22} + \kappa_{20}\kappa_{02} + 2\kappa_{11}^2$.

See K+S p.105, p.102 following (3.69) and p.87.

[Lastly, Goldberger (1991), A Course in Econometrics, p.108 gives an expression for $V(S_{XY})$ that is incorrect. Specifically it contains a term $2(n-1)(\mu_{20}\mu_{02})/n^3$ that should instead be $(n-1)(\mu_{11}^2+ \mu_{20}\mu_{02})/n^3$.]