Gen_Samp = function(seed, X, TS = 1000){
  
  set.seed(seed)
  
  L_in = matrix(0, 7, 2)
  L_in[1:4,1] = sample(high.loading, 4)
  L_in[1:4,2] = sample(low.loading, 4)
  L_in[5:7,1] = sample(low.loading, 3)
  L_in[5:7,2] = sample(high.loading, 3)
  
  Phi_in = phi0_82_pb[,,sample(1:82,1)]
  
  Root = chol(Phi_in)
  
  L_in2 = L_in %*% t(Root)
  
  p = 7
  m = 2
  
  M_Target = matrix(0,p,m)
  M_Target[1:4,1] = NA
  M_Target[5:7,2] = NA
  
  Rotated = targetQ(L=L_in2, Target=M_Target)
  
  Order.M = L_in
  Input.M = rbind(Rotated$loadings, Rotated$Phi)
  
  Out = Align.Matrix(Order.M, Input.M) # A function of the R package EFAutilities
  
  lambda = Out[1:7,]
  Phi = Out[8:9,] 
  
  comm = diag(lambda %*% Phi %*% t(lambda))
  
  if (any(comm > 0.8)){
    
    ind = which(comm > 0.8)
    
    for (ii in 1:length(ind)){
      c_temp = runif(1, 0.6, 0.8)
      lambda[ind[ii],] = lambda[ind[ii],] * sqrt(c_temp/comm[ind[ii]])
    }
  }
    
  tmp = mvtnorm::rmvnorm(1, (beta[c(1:2)] + beta[c(3:4)]*X), sigma_A)
  AR = diag(c(tmp))
  
  while (any(c(tmp) > 0.95)){
    tmp = mvtnorm::rmvnorm(1, (beta[c(1:2)] + beta[c(3:4)]*X), sigma_A)
    AR = diag(c(tmp))
  }
    
  Psi = Phi - AR%*%Phi%*%t(AR)
  
  f_m = matrix(NA, c(TS + 100), 2)
  f_m[1,] = 0
  
  
  for (j in 2:c(TS + 100)){
    f_m[j,] = t(AR %*% f_m[(j-1),]) + mvtnorm::rmvnorm(1, c(0,0), Psi)
  }
  
  f = f_m[101:c(TS+100),]
  
  Du = diag(1 - diag(lambda %*% Phi %*% t(lambda)))
  
  e = matrix(mvtnorm::rmvnorm(TS, rep(0,7), Du), TS, 7)
  y = t(lambda %*% t(f)) + e
  
  res = list(X = X, lambda = lambda, Phi = Phi, AR = AR, Psi = Psi, f = f_m, y = y)
  
  return(res)
}