General ideas are from [here][1] and [here][2]. #R code 1. You first need to simulate a vector of uncorrelated Gaussian random variables, $\bf Z $ NVariables=5 VariableLen=1000 Z=matrix(rnorm(NVariables*VariableLen), ncol=NVariables) 2. Create covariance matrix $\Sigma$. Let all variables be correlated with neighbor as 0.5. Sigma=matrix(data=0, ncol=NVariables-1,nrow=NVariables-1) diag(Sigma)<-0.5 Sigma=cbind(matrix(data=0,nrow=NVariables-1),Sigma) Sigma=rbind(Sigma,matrix(data=0,ncol=NVariables)) diag(Sigma)<-1 3. Find a square root of $\Sigma$. Cholesky decomposition is common choice C=chol(Sigma) 4. To obtain random variables with given correlation matrix $\Sigma$ multiply $\bf C$ and $\bf Z$ Y=Z%*%C 5. Use inverse CDF method to obtain any distribution You wish. Here it is lognormal Ylog=qlnorm(pnorm(Y)) #Results [![enter image description here][3]][3] **Correlation Matrix** [,1] [,2] [,3] [,4] [,5] [1,] 1.00000000 0.52817152 -0.01887624 -0.07113405 -0.05551355 [2,] 0.52817152 1.00000000 0.49392903 -0.03233261 -0.01504632 [3,] -0.01887624 0.49392903 1.00000000 0.50604908 0.04029076 [4,] -0.07113405 -0.03233261 0.50604908 1.00000000 0.49000229 [5,] -0.05551355 -0.01504632 0.04029076 0.49000229 1.00000000 [1]: https://math.stackexchange.com/questions/446093/generate-correlated-normal-random-variables [2]: http://stats.stackexchange.com/questions/184325/how-does-the-inverse-transform-method-work [3]: https://i.sstatic.net/OKjTf.png