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.
 A: 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:


*

*Jointly resample from observed rows of data, allowing for replications and holding the sample size fixed.

*Re-estimate covariance in the resampled data.

*Repeat 1-2 for a sufficient number of iterations.

*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.
A: 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.
A: 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$.]
A: 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
