What is the maximum likelihood estimate of the covariance of bivariate normal data when mean and variance are known?

Suppose we have a random sample from a bivariate normal distribution which has zeroes as means and ones as variances, so the only unknown parameter is the covariance. What is the MLE of the covariance? I know it should be something like $\frac{1}{n} \sum_{j=1}^{n}x_j y_j$ but how do we know this?

-
As a starter, don't you think it's a bit unclever to estimate the means with $\bar{x}$ and $\bar{y}$ when in fact we know that they are 0 and 0? – Wolfgang Aug 14 '14 at 23:37
Very unclever, fixed it. Still don't see how this can easily follow. It's analogous to the sample variance but why is it the MLE (unless it's not and I made another mistake) – Stacy Aug 14 '14 at 23:55
Have you deleted $\frac{1}{n}\sum_{i=1}^n(x_i-\bar x)(y_i - \bar y)$ ? Taking this formula does not mean you consider $\bar x$ and $\bar y$ as the estimates of the means. – Stéphane Laurent Aug 15 '14 at 14:12
@StéphaneLaurent Yes, in the initial post, the formula was given as you have written it. – Wolfgang Aug 15 '14 at 17:25

The estimator for the correlation coefficient (which in the case of a bivariate standard normal equals the covariance)

$$\tilde r = \frac 1n\sum_{i=1}^nx_iy_i$$

is the Method-of-Moments estimator, the sample covariance. Let's see if it coincides with the maximum likelihood estimator, $\hat \rho$.

The joint density of a bivariate standard normal with correlation coefficient $\rho$ is

$$f(x,y) = \frac{1}{2 \pi \sqrt{1-\rho^2}} \exp\left\{-\frac{x^2 +y^2 -2\rho xy}{2(1-\rho^2)}\right\}$$

and so the log-likelihood of an i.i.d. sample of size $n$ is

$$\ln L = -n\ln(2\pi) -\frac n2\ln(1-\rho^2) - \frac 1{2(1-\rho^2)}\sum_{i=1}^n(x_i^2 +y_i^2 -2\rho x_iy_i)$$

(here the i.i.d assumption is with respect to each draw from the two-dimensional population of course)

Taking the derivative with respect to $\rho$ and setting it equal to zero gives a 3d-degree polynomial in $\rho$:

$$\hat \rho: n\hat \rho^3 -\left(\sum_{i=1}^nx_iy_i\right)\hat\rho^2 -\left(1- \frac 1n\sum_{i=1}^n(x_i^2 +y_i^2) \right)n\hat \rho - \sum_{i=1}^nx_iy_i =0$$

That the calculations are correct can be verified if one takes the expected value of the derivative evaluated at the true coefficient $\rho$ -it will equal zero.

For compactness, write $(1/n)\sum_{i=1}^n(x_i^2 +y_i^2) = (1/n)S_2$, which is the sum of the sample variances of $X$ and $Y$. If we divide the 1st-derivative expression by $n$ the MoM estimator will appear, specificaly

$$\hat \rho: \hat \rho^3 -\tilde r \hat \rho^2 + \big[(1/n)S_2-1\big]\hat \rho -\tilde r=0$$

$$\Rightarrow \hat \rho\Big(\hat \rho^2 -\tilde r \hat \rho + \big[(1/n)S_2-1\big] \Big) = \tilde r$$

Doing the algebra, it is not difficult to conclude that we will obtain $\hat \rho = \tilde r$ if,and only if, $(1/n)S_2 =2$, i.e. only if it so happens that the sum of sample variances equals the sum of true variances. So in general

$$\hat \rho \neq \tilde r$$

So what happens here? Somebody wiser will explain it, for the moment, let's try a simulation: I generated an i.i.d. sample of two standard normals with correlation coefficient $\rho=0.6$. The sample size was $n=1.000$. The sample values were

$$\sum_{i=1}^nx_iy_i = 522.05,\;\;S_2 = 1913.28$$

The Method-of-Moments estimator gives us

$$\tilde r = \frac {522.05}{1000} = 0.522$$

What happens with the log-likelihood? Visually, we have

Numerically, we have

\begin{array}{| r | r | r |} \hline \hline ρ&\text{1st deriv}&\text{lnL}\\ \hline 0.5&-70.92&-783.65\\ 0.51&-59.41&-782.47\\ 0.52&-47.7&-781.48\\ 0.53&-35.78&-780.68\\ 0.54&-23.64&-780.1\\ 0.55&-11.29&-779.75\\ 0.56&1.29&-779.64\\ 0.57&14.1&-779.81\\ 0.58&27.15&-780.27\\ 0.59&40.44&-781.05\\ 0.6&53.98&-782.18\\ \hline \end{array}

and we see that the log-likelihood has a maximum a tad before $\rho=0.56$ where also the 1st derivative becomes zero $(\hat \rho = 0.558985)$. No surprises for the values of $\rho$ not shown. Also, the 1st derivative has no other root.

So this simulation accords with the result that the maximum likelihood estimator does not equal the method of moments estimator (which is the sample covariance between the two r.v.'s).

But it appears that "everybody" is saying that it should... so somebody should come up with an explanation.

UPDATE

A reference that proves that the MLE is the Method-of-Moments estimator: Anderson, T. W., & Olkin, I. (1985). Maximum-likelihood estimation of the parameters of a multivariate normal distribution. Linear algebra and its applications, 70, 147-171.
Does it matter that here all means and variances are free to vary and not fixed?

... Probably yes, because @guy's comment in another (now deleted) answer says that, with given mean and variance parameters, the bivariate normal becomes a member of the curved exponential family (and so some results and properties change)... which appears to be the only way that can reconcile the two results.

-
This is a little surprising, but after some reflection it should be expected. The problem can be rephrased as estimating the regression coefficient $\rho$ in the model $Y = \rho X + \epsilon$ where $\epsilon \sim \mathcal N(0, \sqrt{1 - \rho^2}^2)$. This isn't a linear model, so there is no reason to expect the MLE to be a simple dot product. The same logic shows (I think!) that if we only know $\mbox{Var}(X)$ then the MLE is $x'y / x'x$, and $x'y / y'y$ if we only know $\mbox{Var}(Y)$. If we happen to know neither, we get your MOM estimator. – guy Aug 15 '14 at 4:32
@guy: Very interesting. I think these arguments, if slightly expanded, fully deserve to be posted as a separate answer! – amoeba Aug 15 '14 at 10:17
@guy I don't think this formulation is equivalent, because, the log-likelihood in the regression set up contains the square $\epsilon^2=(y-\rho x)^2 = y^2 -2\rho xy + \rho^2 x^2$. The $\rho^2$ coefficient attached to $x^2$ is not present in the bivariate density formulation. – Alecos Papadopoulos Aug 15 '14 at 11:28
My guess is $\frac{1}{n}\sum_{i=1}^n(x_i-\bar x)(y_i - \bar y)$. Imagine $n=2$ and $y_1=y_2$, then a $0$ estimate is expected. – Stéphane Laurent Aug 15 '14 at 14:09
hmmmmm.... no I don't think so cause the estimate should lie between $-1$ and $1$ – Stéphane Laurent Aug 15 '14 at 14:37

Under the stated conditions ($\mu_X = \mu_Y = 0$ and $\sigma_X = \sigma_Y = 1$), the likelihood function for a random sample of size $n$ is $$L(\rho\; |\; X, Y) = \frac{1}{(2\pi[1-\rho^2])^{n/2}}\exp \left[-\frac{1}{2(1-\rho^2)}(X'X - 2\rho X'Y + Y'Y)\right].$$

Now find the log-likelihood and take the derivative with respect to $\rho$. Next, set it equal to 0, solving for $\hat{\rho}$. You should of course do some appropriate test to show what you've found is in fact a global maximum.

-