Assume we have n observations $x_i$ (i from 1 to n), each from the a normal distribution with mean 0 and some variance component: $X_i \sim N(0, \sigma^2)$. The random variables $X_i$s have some (let's assume known) correlation structure. I.e.: $Corr(X_i, X_j) = \rho_{ij}$.
How should we estimate $\sigma^2$?
How would the correlated structure impact an estimator such as $\hat \sigma^2 = \frac{\sum(x_i^2)}{n}$? Could we devise a better estimator (given the correlation structure is known)?
Let $R = (\rho_{ij})$ be the correlation matrix so that the covariance matrix is $\Sigma = \sigma^2 R.$ Consider $\mathbf x = (x_1,\ldots, x_n)^\prime$ to be a single observation of the $n$-variate Normal distribution with zero mean and $\Sigma$ covariance. Because the log likelihood of $N\ge 1$ independent such observations is, up to an additive constant (depending only on $N,$ $n,$ and $R$) given by
$$\Lambda(\sigma) = \sum_{i=1}^N (-n \log \sigma) - \frac{1}{2\sigma^2} \sum_{i = 1}^N\mathbf{x}_i^\prime\, R^{-1}\,\mathbf{x}_i,$$
it has critical points where $\sigma\to 0,$ $\sigma\to\infty,$ and at any solutions of
$$0 = \frac{\mathrm d}{\mathrm{d}\sigma}\Lambda(\sigma) = -\frac{nN}{\sigma} + \frac{1}{\sigma^3} \sum_{i = 1}^N\mathbf{x}_i^\prime\, R^{-1}\,\mathbf{x}_i.$$
Unless the form $ \sum_{i = 1}^N\mathbf{x}_i^\prime\, R^{-1}\,\mathbf{x}_i$ is zero, there is a unique global maximum at
$$\hat\sigma^2 = \frac{1}{nN} \sum_{i = 1}^N \mathbf{x}_i^\prime\, R^{-1}\,\mathbf{x}_i.$$
This is the Maximum Likelihood estimate. It exists even when $N=1$ (which is the situation posited in the question).
Intuitively, this will be superior to any estimate that ignores the correlations assumed in $R.$ To check, we could compute the Fisher Information matrix for this estimate -- but I will leave that to you, in part because the result shouldn't be convincing for small values of $N$ (where the maximum likelihood asymptotic results might not apply).
To illustrate what actually happens, here are one thousand estimates of $\sigma^2$ from experiments with $n=4$ and $N=1.$ The value of $\sigma$ was set to $1$ throughout. In these experiments, the correlation was always
$$R = \pmatrix{1 &0.3055569 &0.5513377 &0.5100989\\ 0.3055569 &1 &0.1240151 &0.09634469\\ 0.5513377 &0.12401511 &1 &-0.4209064\\ 0.5100989 &0.09634469 &-0.4209064 &1}$$
(as generated randomly at the outset). In the figure the leftmost panel is a histogram of the foregoing Maximum Likelihood estimates; the middle panel is a histogram of estimates using the usual (unbiased) variance estimator; and the right panel is a QQ plot of the two sets of estimates. The slanted line is the line of equality. You can see the usual variance estimator tends to yield more extreme values. It is also biased (due to ignoring the correlation): the mean of the MLEs is 0.986 -- surprisingly close to the true value of $\sigma^2 =1^2 =1$ while the mean of the usual estimates is only 0.791. (I write "surprisingly" because it is well-known the usual maximum likelihood estimator of $\sigma^2,$ where no correlation is involved, has a bias of order $1/(nN),$ which is pretty large in this case.)
You may experiment with the R code that produced these figures by modifying the values of n, sigma, N, n.sim, Rho, and the random number seed 17.
f <- function(x, Rho) { # The MLE of sigma^2 given data `x` and correlation `Rho`
S <- solve(Rho)
sum(apply(x, 1, function(x) x %*% S %*% x)) / length(c(x))
}
n <- 4
sigma <- 1
N <- 1
n.sim <- 1e3
set.seed(17)
#
# Generate a random correlation matrix. Larger values of `d` yield more
# spherical matrices in general.
#
d <- 1
Rho <- cor(matrix(rnorm(n*(n+d)), ncol=n))
(ev <- eigen(Rho, only.values=TRUE)$values)
#
# Run the experiments.
#
library(MASS)
sim <- replicate(n.sim, {
x <- matrix(mvrnorm(N, rep(0,n), sigma^2 * Rho), N)
c(f(x, Rho), var(c(x)))
})
(rowMeans(sim))
#
# Plot the results.
#
par(mfrow=c(1,3))
hist(sim[1,], col=gray(.93), xlab="Estimate",
main=expression(paste("Histogram of Estimates of ", sigma^2)))
abline(v = sigma^2, col="Red", lwd=2)
hist(sim[2,],col=gray(.93), xlab="Estimate",
main=expression(paste("Histogram of Independent Estimates of ", sigma^2)))
abline(v = sigma^2, col="Red", lwd=2)
y1 <- sort(sim[1,])
y2 <- sort(sim[2,])
plot(y1, y2,asp=1, xlab="Correlation-based estimate", ylab="Independent estimate")
abline(0:1, col="Red", lwd=2)
par(mfrow=c(1,1))
You can apply GLS equations, e.g. top of p.5 in these lectures: http://halweb.uc3m.es/esp/Personal/personas/durban/esp/web/notes/gls.pdf
The equation for MSE is $\sigma^2=1/n \sum_{ij} x_i r^{-1}_{ij} x_j$ if you set X to zero and correspond your r to their V
This equation doesn’t work for extreme cases such as correlation 1 between all variables. This case your usual equation will still produce a number
