nobs <- 500
nsim <- 1000

for (i in 1:nsim) {
  x1 <- rnorm(nobs)
  x2 <- rnorm(nobs)
  
  x.matrix <- cbind(x1,x2)%*%chol(matrix(c(0.085,0.093,0.093,0.907), ncol=2))
  f1 <- x.matrix[,1]
  f2 <- x.matrix[,2]
  
  y <- 1.074*f1 + 0.109*f2 + rnorm(nobs, mean=1.554, sd=sqrt(0.871))
  
  ex1 <- 1*f1 + rnorm(nobs, mean=1.146, sd=sqrt(0.207))
  ex2 <- 1.086*f1 + rnorm(nobs, mean=1.122, sd=sqrt(0.203))
  ex3 <- 1.624*f1 + rnorm(nobs, mean=1.393, sd=sqrt(0.553))
  ex4 <- 1.262*f1 + rnorm(nobs, mean=1.203, sd=sqrt(0.298))
  in1 <- 1*f2 + rnorm(nobs, mean=2.286, sd=sqrt(0.945))
  in2 <- 0.747*f2 + rnorm(nobs, mean=1.956, sd=sqrt(0.860))
  in3 <- 0.888*f2 + rnorm(nobs, mean=1.983, sd=sqrt(0.865))
  in4 <- 0.675*f2 + rnorm(nobs, mean=1.636, sd=sqrt(0.631))
  
  obs.matrix <- cbind(ex1,ex2,ex3,ex4,in1,in2,in3,in4,y)
  print(i)
  write.table(obs.matrix, 
              file=paste("z:/mplus/mvn/mvn",i,".dat",sep=""), 
              row.names=FALSE, col.names=FALSE)
}

write.table(paste("mvn", 1:nsim, ".dat", sep=""),
            file="z:/mplus/mvn/mvnlist.dat",
            row.names=FALSE, col.names=FALSE, quote=FALSE)