# Standard error of estimated covariance

Let $$X_1,...,X_n$$ and $$Y_1,...,Y_n$$ be two independent random samples from $$\mathcal{N}(\mu, \sigma^2)$$ where both $$\mu$$ and $$\sigma$$ are unknown parameters.

I estimate their covariance using: $$\hat{\operatorname{cov}}(X, Y) = \operatorname{E}{\big[(X_i - \operatorname{E}[X])(Y_i - \operatorname{E}[Y])\big]}$$

with replacing $$\operatorname{E}[X]$$ and $$\operatorname{E}[Y]$$ by the according sample mean.

How do i calculate the standard error of $$\hat{\operatorname{cov}}(X, Y)$$?

Edit: The comment by Adam L. Taylor is valid. If both variables are known to be independent, there would be no need to estimate their covariance; so this assumption should be dropped.

• This question may be of interest to you as well: stats.stackexchange.com/questions/287144/… Commented Nov 19, 2020 at 16:09
• Surely you don't mean that $X$ and $Y$ are independent. If that were the case, the perfect estimator of covariance is simply the function that always returns 0, and the standard deviation of this estimator is 0. If you know a priori that $X$ and $Y$ are independent, you already know their covariance! Commented Nov 20, 2020 at 4:12
• Good point @AdamL.Taylor, however OP is not asking for the variance of the 'perfect' estimator, but of a specific one. Commented Nov 22, 2020 at 13:56
• So I think you want to say that you have $n$ samples from a bivariate normal distribution, with common unknown mean $\mu_x=\mu_y=\mu$, common unknown variance $\sigma_x^2=\sigma_y^2=\sigma^2$, and unknown covariance $\rho\sigma^2$. Commented Nov 24, 2020 at 21:47
• @AdamL.Taylor: To be clear, i just have one drawn sample from a bivariate normal distribution, having $n$ observation for both $X$ and $Y$. My best estimate for the covariance is the one (and only) observation i have and would like to know its standard error. Commented Nov 27, 2020 at 11:43

The independence of $$X$$ and $$Y$$ makes this problem straightforward. To make the notation easier, assume $$\mu=0$$. Then

$$Cov (X,Y)=S_{XY}=E[(X-\mu_X)(Y-\mu_Y)]=E[XY]$$,

and the estimator $$\hat S_{XY} = \frac{1}{n}\sum_{i=1}^n x_iy_i$$ has expectation zero, so $$Var (\hat S_{XY}) = E[\hat S_{XY}^2]$$.

\begin{align} Var(\hat S_{XY}) &= E[\hat S_{XY}^2] \\ &= E\left[\left(\frac{1}{n}\sum_{i=1}^n x_iy_i\right)^2\right] \\ &= \frac{1}{n^2}E\left[\sum_{i=1}^n x_i^2y_i^2+2\sum_{i

So the standard error of $$\hat S_{XY} = \sqrt{Var(\hat S_{XY})}=\sigma^2/\sqrt{n}$$.

That's an interesting one :-)

Now if I understood your question right, then the trick is to think in terms of functions and not focus just on COV. Though, the following things are initially important:

• It's normal distributed.
• Covariance is (just) a function as any other.

If you need the standard deviation for the result of the COV function, you automatically assume that:

• The elements $$X_i, Y_i$$ might have standard deviations associated to these measurement values.
• The standard deviation $$E[X], E[Y]$$ is already known and can be computed by mean($$X$$),mean($$Y$$). Do you have any reason not to trust it by assuming a different value?

Now given this information, you can use the standard method for computing the resulting error-estimate of an function with Gaussian Propagation of Uncertainty. Important is the limitation; It works just for Normal-distributed variables. The variance $$\sigma_y^2$$ of an variable $$y$$ which consist of other uncertain variables $$x$$ and their corresponding variances $$\sigma_x^2$$, such as

$$y = x_1 + x_2 + ...+ x_n$$

$$\sigma_{y}^2 = \sigma_{x_1}^2+\sigma_{x_2}^2+...+\sigma_{x_n}^2$$

can be computed as (matrix notation):

$$\sigma_y^2 = \mathbf{A\Sigma}_{xx}\mathbf{A}^\mathrm{T}$$.

where $$\mathbf{A}$$ is the Jacobian matrix and $$\mathbf{\Sigma}_{xx}$$ is the variance-covariance matrix for the values $$X_i,Y_i$$ corresponding to the function. On the diagonal you need to place the variances for $$X_i,Y_i$$ the off-diagonal values are covariances between them (you might want to assume 0 for them). Please keep in mind, this is an general solution for non-linear functions and uses just one (the first) linearization term. It is fast and usually the way to go in productive applications but might have approximation errors compared to a pure analytical solutions.

Another option is to do a small Monte-Carlo simulation. In order to achieve this you can sample $$X_i,Y_i$$ with their expected uncertainty and compute their covariance. Now if you do it several (thousand) times, you get a fair estimate for the resulting error. Here is a pseudo-code for OCTAVE/MATLAB:

% Clean stuff before start to avoid variable conflicts
clc
clear all

% These are the values
X = [ 1 2 3 4 5 ].';
Y = [ 5 4 3 2 1 ].';

% How many tries do you want to have
n_samples = 10000;

% prepare the resulting error
cov_res = zeros( n_samples , 1 );

% loop the computation through n_samples
for i = 1 : n_samples

% generate random distributed noise, 1 sigma [-0.1:0.1]
x_error_sample = 0.1 * randn( size( X , 1 ) , 1 );
y_error_sample = 0.1 * randn( size( Y , 1 ) , 1 );

% Compute the covariance matrix for X and Y
cov_i = cov( X + x_error_sample , Y + y_error_sample );

% Pick only the covariance
cov_res( i ) = cov_i( 1 , 2 );
end

% covariance estimator can be chosen by your own metric (e.g. mean,median,...)
mean( cov_res )

% The error of this estimation can be chosen by your own metric (e.g. std,rms,var,...)
std( cov_res )


This approach might also be used for any distribution for $$X$$ and $$Y$$, just replace the term randn with your choice.

Regards

In addition to @abstrusiosity's analytic solution, you can apply the bootstrap here. This has the advantage of working even when the two samples aren't independent (the true covariance isn't $$0$$)

library(tidyverse)

# Simulate data
m = 0
s = 2
n = 100
X = rnorm(n, m, s)
Y = rnorm(n, m, s)

(expected_se = (s**2) / sqrt(n)) # @abstrusiosity's solution
# [1] 0.2828427

nboot = 200
bootstrap_cov = map_dbl(1:nboot, function(i){
rx = sample(X, n, replace = T)
ry = sample(X, n, replace = T)
cov(rx, ry)
})
sd(bootstrap_cov)
# [1] 0.270266


You can also test that the bootstrap estimate converges to the analytic solution:

sim_bootstrap = function(i){
# Simulate fresh data
X = rnorm(n, m, s)
Y = rnorm(n, m, s)
bootstrap_cov = map_dbl(1:nboot, function(i){
# Do bootstrap sample
rx = sample(X, n, replace = T)
ry = sample(X, n, replace = T)
cov(rx, ry)
})
sd(bootstrap_cov)
}
results = map_dbl(1:200, sim_bootstrap)
mean(results)
# [1] 0.2817101

qplot(results) + geom_vline(xintercept=expected_se, color='red') +
labs(x='Standard Error', caption='Red line shows analytic SE')