Ok, let's make sure I understand properly here - you are looking to generate a multivariate normal variable N(mu, sigma) with sigma = (1 .5 , .5 1) so you generate two independant normal random vars and use the covar matrix to generate the resulting bivariate normal random variable? So in the end you wind up with a rescaled value X+.5Y in the same row as your X simulation and .5X+Y in the row of your Y simulation for each simulation pair. As Snoopy pointed out, this can probably be done in moreorless a single pass of data which would save you tremendous amounts of time for processing. The bigger question would be do you need to build a generic process that will be able to further handle N-variate random variables or are you looking for a solution tailored to your problem? A simple approach to save time if you don't need an expandable solution is to do everything in one line and hard the multiplication since its easy to achieve when N = 2 dimensions. E.g. data want;      do sim = 1 to 100;           Z_A = rand('norm', 0, 1);           Z_B = rand('norm', 0, 1);           W_A=Z_A+.5*Z_B;           W_B=.5*Z_A+Z_B;           output;      end; run; /* you could use transpose if you really needed the bivariate data over 2 rows */ Assuming you want this to be somewhat generic, here's how I would approach the problem: proc fcmp OUTLIB=WORK.FUNCS.MATRIX;      subroutine mmult(z {*,*}, corr {*,*}, y {*,*});      outargs y;      call mult(z, corr, y);      endsub; run;    options cmplib=work.funcs; %macro multivarnormsim(numsim=, nvar=, corrmat= ); data want;      array z {&nvar.} z1-z&nvar.;   array zz {1, &nvar.} _TEMPORARY_;      array corr {&nvar., &nvar.} _temporary_ (&corrmat);      array w {&nvar.} w1-w&nvar.;   array ww {1, &nvar.} _TEMPORARY_;      do i=1 to &numsim;           do j=1 to dim(z);                z{j} = rand('norm', 0, 1);       zz{1, j}=z{j};           end;           call mmult(zz, corr, ww);     do j=1 to dim(w);    w{j}=ww{1, j};     end;     output;      end;      drop i j; run; %mend; %multivarnormsim(numsim=100, nvar=5, corrmat=     %quote( 1 0 0 0 .5,                                                                                       0 1 .2 .2 0,                                                                                        0 .2 1 0  0,                                                                                        0 .2 0 1 0,                                                                                        .5 0 0 0 1)                               ); /*the nvar*nvar values could've been written in a single line it doesn't matter */ ... View more