Marginal Ricians with correlated errors

First, what's a good way to generate errors with a desired correlation structure? Let's start with the following:

> d <- as.matrix(dist(1:4, diag=T, upper=T))
> s <- 0.3^d
> a <- chol(s)
> m <- matrix(nrow=10000,ncol=4)
> for (i in 1:10000) { x <- rnorm(4); m[i,] <- a %*% x}
> var(m)
           [,1]       [,2]      [,3]       [,4]
[1,] 1.09839578 0.30672920 0.1062927 0.03032918
[2,] 0.30672920 1.01055057 0.3031668 0.09055224
[3,] 0.10629268 0.30316682 1.0044863 0.27050059
[4,] 0.03032918 0.09055224 0.2705006 0.88942696
> s
      1    2    3     4
1 1.000 0.30 0.09 0.027
2 0.300 1.00 0.30 0.090
3 0.090 0.30 1.00 0.300
4 0.027 0.09 0.30 1.000

Is it surprising to find the sample variance such a relatively poor approximation to the true variance even with 10000 observations? (I guess the answer is "no!")

Now, let's use this to see what happens to a mixture of marginally Rician distributions as the correlation gets cranked up.

##Produce elements that are marginally Rician distributed but for
##which the noise between elements is correlated with the
##structure cor(X_i, X_j) = cor(Y_i, Y_j) =  rho^abs(i-j)
##but {X_i} independent of {Y_i}, i = 1:length(mu) and
##result[i] = sqrt(X_i^2 + Y_i^2).
##The value is a matrix of draws from the given distribution,
##each row a single draw so that the ith column corresponds
##to mu[i]
corRice <- function(mu, n, rho, sigma) {
  k <- length(mu)
  d <- as.matrix(dist(1:k, diag=TRUE, upper=TRUE))
  r <- (rho ^ d)
  a <- sigma * chol(r)
  result <- matrix(nrow=n, ncol=k)
  for (i in 1:n) {
    x <- (a %*% rnorm(k)) + mu
    y <- a %*% rnorm(k)
    result[i, ] <- sqrt(x^2 + y^2)
  }
  result
}

> r0 <- corRice(rep(c(0, 50, 100), each=50), n=10000, rho=0, sigma=10)
> r1 <- corRice(rep(c(0, 50, 100), each=50), n=10000, rho=0.9, sigma=10)
> r2 <- corRice(rep(c(0, 50, 100), each=50), n=10000, rho=0.99, sigma=10)
> plot(quantile(r0, probs=(1:999)/1000), quantile(r1, probs=(1:999)/1000), xlab="quantiles with rho=0", ylab="quantiles with rho=0.9")
> plot(quantile(r0, probs=(1:999)/1000), quantile(r2, probs=(1:999)/1000), xlab="quantiles with rho=0", ylab="quantiles with rho=0.99")



Hmmmmmmmm.