#============================================================
#function to execute the simulation
#arguments:
#lambda: matrix of factor loadings of dimension (items x factors)
#phi: covariance matrix of factors
#theta: covariance matrix of measurement errors
#epsilon: vector of thresholds
#n: the sample size.
#iter: the number of iterations.
montecarlo = function(lambda, phi, theta, epsilon,sample,iter) {

require(MASS) #Load library MASS which generates multivariate normal distributions.

 #create an empty file to store names of datasets for MPLUS to analyze
 write.table(c(), file = "filelist.txt", append = F, quote=F, row.names=F, col.names=F)
sigma <- lambda %*% phi %*% t( lambda) + theta #obtain implied covariance matrix

for (i in 1:iter) { #loop through iterations

#create a single dataset
dataset <- mvrnorm(sample, mu = rep(0,nrow(sigma)),sigma)
#dichotomize continuous indicators
epsilon <- matrix(epsilon, nrow(dataset),ncol(dataset),byrow=T)
dataset <- (dataset>=epsilon)*1

#save data
filename <- paste("data",i,".csv",sep="") #create file name
write.table(filename, file = "filelist.txt", append = T, quote=F, row.names=F, col.names=F)
write.table(dataset,file=filename,sep=",",quote=F,row.names=F, col.names=F)

  } } #close loop through iterations and close function