Sample mean and variance independence in the case of correlated observations Are the sample mean and sample variance of correlated normal observations independent? 
The classic theory relies on the independence of observations, but this is not the case. So as an example I tried to simulate 100k samples of 5 elements taken from a stationary, ergodic, normal random process with known mean and covariance, and I evaluated the cross distribution between the sample mean and the sample variance and I got this:
So from this (and others) simulations I'm starting to believe that sample mean and sample variance for a normal process are independent even if the samples are correlated, but I don't know how to prove this. Any help?
 A: The sample mean of a multivariate normal vector $\mathbf{X}=(X_1, X_2, \ldots, X_n)$ is a function of 
$$M = X_1+X_2+\cdots X_n$$
and the sample variance is a function of the residual vector with components
$$Z_i = -X_1 - X_2 - \cdots - X_{i-1} + (n-1)X_i - X_{i+1} - \cdots - X_n,$$
$i=1, 2, \ldots, n$.
Let $\Sigma$ be the covariance matrix of $\mathbf{X}$.  Write $\sigma_i$ for the sum of column (or row) $i$ of $\Sigma$, $\sigma_i = \Sigma_{1i} + \Sigma_{2i} + \cdots + \Sigma_{ni}$, and let $\sigma$ be the sum of all the entries of $\Sigma$.   We may compute
$$\operatorname{Cov}(M, Z_i) = n\sigma_i - \sigma.$$
Because both $M$ and $Z_i$ are linear combinations of multivariate Normal variables, they are jointly Normal, whence they are independent if and only if their covariance is zero.  Consequently $M$ is independent of all the $Z_i$ if and only if
$$n\sigma_1 = n\sigma_2 = \cdots = n\sigma_n = \sigma.$$
In other words, equality of the column sums guarantees independence of the mean and the components of the sample variance, whence it will guarantee independence of the mean and the sample variance itself.
Although the converse is not true--it is possible for $M$ not to be independent of the $Z_i$, yet for $M$ to be independent of the sample mean--this requires exceptional circumstances.  In almost all cases, inequality of the column sums creates a dependence between the sample mean and sample standard deviation.
By definition, in a stationary process the covariances $\Sigma_{ij}$ may depend only on $i-j$.  Although this does not guarantee the column sums are all equal, for large $n$ and a covariance that decays sufficiently rapidly with $|i-j|$, it will approximately be true, because in the limit the column sums are all equal:
$$\sigma_i = \sum_{j=-\infty}^\infty \Sigma_{ji} =\sum_{j=-\infty}^\infty \Sigma_{jk} =\sigma_k.$$
All that is required is the convergence of these sums.

A good way to see the dependence in the scatterplot is to render the points more carefully.  When they are made semitransparent, you can see the underlying density better.  A lowess smooth helps demonstrate a variation in the standard deviation with the mean in this example where $n=8$ and the column sums of $\Sigma$ vary appreciably.

Here is the R code that generated it.
library(MASS)   # mvrnorm()
set.seed(17)
n <- 5e4        # Simulation size
d <- 8          # Dimension
k <- 4          # Size of upper block of Sigma
rho <- 0.99     # Correlation in upper block
mu <- rep(0, d) # Mean
Sigma <- outer(1:d, 1:d, function(i,j) ifelse(i <= k & j <= k, rho^abs(i-j), i==j))
colSums(Sigma)

x <- mvrnorm(n, mu, Sigma)
sim <- t(apply(x, 1, function(y) c(mean(y), sd(y))))

plot(sim, pch=16, cex=0.5, col="#00000008",
     xlab="Mean", ylab="SD")
i <- order(sim[, 1])
lines(sim[i, 1], lowess(sim[i, 2], f=1/20)$y, col="Red", lwd=2)
# g <- cut(sim[, 1], quantile(sim[, 1], seq(0, 1, by=0.025)))
# boxplot(sim[, 2] ~ g)

A: Only if $\Sigma$ is diagonal it appears. (see edit). I'm going to look at the covariance between $[\bar{X}, \ldots, \bar{X}]'$ and $X -  [\bar{X}, \ldots, \bar{X}]'$. It's easier to work with the unsquared deviates rather than the sum of squared deviates.
Let $X \sim \text{Normal}(\mu, \Sigma)$, where $X, \mu \in R^n$. Let $\Sigma$ look like whatever you want it to. Not diagonal, because you don't want independence between samples. Just as long as it's symmetric and positive definite (I'm assuming we're dealing with full rank normals here). 
Then I use bilinearity of $Cov(\cdot, \cdot)$. 
$Cov(1\frac{1}{n}1'X, (I-\frac{1}{n}11')X) = \frac{1}{n}11' \Sigma (I-\frac{1}{n}11') = \frac{1}{n}11' \Sigma-\frac{1}{n}11' \Sigma  \frac{1}{n}11' = \frac{1}{n}11' \Sigma-\frac{1}{n^2}11' \Sigma 11'$
As whuber points out, 0 covariance/correlation only implies independence in the case of normal vectors. This can be shown by writing down the density and seeing it factors.
Edit: whuber is correct that my final conclusion is incorrect. He came up with a more general criterion that guarantees independence (+1). Below I continue on verifying his answer with my notation. Let $e_i = (0,\ldots, 1,\ldots, 0)'$ be the vector with a one in the $i$th spot. Whuber's condition that the row or column sums be equal is equivalent to assuming that $\Sigma e_i = \Sigma e_k$ for $i \neq k$ or $e_i' \Sigma = e_k' \Sigma$. If we re-write $\Sigma 1$ as $ \sum_{i=l}^n \Sigma e_l$ then $\Sigma 1 = n \Sigma e_1$ and $Cov(1\frac{1}{n}1'X, (I-\frac{1}{n}11')X)$ becomes diagonal using the derivation above. Side note: apologies for overloading $\Sigma$.
