How to simulate n Uniform[0,1] variables with a specified correlation? For a simulation study, I need to generate $n$ Uniform$[0,1]$ random variables with population correlation $\rho$. I'm not aware of any simple ways to do this. 
I have considered sampling from a multivariate normal with mean 0 and exchangeable correlation matrix $ \Sigma= \left[ {\begin{array}{cccc}
   1 & \rho &\dots  & \rho \\
   \rho & 1 & \dots  & \rho \\
   \vdots & \vdots & \vdots & 1
  \end{array} } \right] \\ $ 
then using an inverse normal transformation, but the correlation will not be $\rho$.
 A: Your method works.
Take a standard bivariate normal random variable $(X_1,X_2)$ with correlation $\tau$.  Let $\Phi$ be the standard Normal distribution function.  The copula $(U_1,U_2)$ given by $U_i = \Phi(X_i)$ has some correlation $\rho=f(\tau)$.  Although this function $f$ cannot be determined analytically (as far as I can tell), numerical integration with a full suite of values $\tau$ spanning the interval $[-1,1]$ suggests its inverse $f^{-1}$, which determines what $\tau$ needs to be to yield a correlation $\rho$, can be approximated to about $0.000001$ throughout the interval by the function
$$\hat f^{-1}(\rho) =  \frac{4118}{3163}  \sin(\rho/3) + \frac{3183}{3149}  \sin(2  \rho/3) - \frac{145}{2391}  \sin(3  \rho/3).$$
Tests with the following R implementation support the claimed accuracy.  (Actually, it would take a large simulation--greater than $10^{12}$ observations--to detect the error of approximation.)
Edit
Let $f_\tau$ be the standard bivariate Normal density,
$$f_\tau(x_1,x_2) = \frac{1}{2\pi \sqrt{1-\tau^2}} \exp\left(-\frac{1}{2(1-\tau^2)}(x_1^2+x_2^2-2\tau x_1x_2)\right).$$
  Then the correlation of the $U_i$ is a function of the moments
$$m_\tau(i,j) = E[U_1^i U_2^j] = E[\Phi(X_1)^i\Phi(X_2)^j] = \iint f_\tau(x_1,x_2) \Phi(x_1)^i \Phi(x_2)^j dx_1 dx_2;$$
$$f(\tau)=\operatorname{Cor}(U_1,U_2) = \frac{m_\tau(1,1) - m_\tau(1,0)^2}{m_\tau(2,0)-m_\tau(1,0)^2}.$$
Most of these have easy analytical expressions based on univariate moments of the uniform distribution, because they do not depend on $\rho$: $m_\rho(1,0)=1/2,$ $m_\rho(2,0)=1/3.$  The one I integrated numerically is $m_\rho(1,1)$.  I enforced the symmetry $f(-\tau)=-f(\tau)$ by numerically computing both $f(\pm \tau)$ and averaging their absolute values.

unif.to.norm <- function(rho) {
  4118/3163 * sin(rho/3) + 3183/3149 * sin(2 * rho/3) - 145/2391 * sin(rho)
}
# curve(unif.to.norm(rho)-rho, 0, 1, xname="rho") # Compare to a linear function
library(MASS) # mvrnorm
n <- 1e5      # observations to simulate
d <- 5        # dimensions
rho <- 0.6    # Must lie in [-1/(d-1), 1]
#
# Compute the correlation matrix for the multivariate Normal distribution.
#
rho.norm <- unif.to.norm(rho)
Sigma <- matrix(rho.norm, d, d) + diag(rep(1-rho.norm, d))
#
# Simulate data.
#
x <- pnorm(mvrnorm(n, rep(0, d), Sigma))
#
# Display the data.
# 
# pairs(x)               # The scatterplots
# par(mfrow=c(1, d))
# apply(x, 2, hist)      # The histograms (will be uniform)
# par(mfrow=c(1,1))
Sigma.hat <- cor(x)
prec <- ceiling(log10(n)/2)
signif(Sigma.hat, prec)  # Should almost match `rho` in all off-diagonal entries
signif(mean(Sigma.hat[lower.tri(Sigma.hat)]), prec+1)

