/* with detailed outputs from simulation size 5 */

option ls=76 nocenter;

%let seed      = 71853048;
%let var_cage  = 5.03;
%let var_cc    = -6.24;
%let var_resid = 31.56; 
%let nsim      = 5;

proc iml;
   cov = &var_cc * j(3,3) + &var_resid * i(3);
   covh = half(cov)`;
   print cov covh;
run;
				
data noise;
   retain i j;
   do sim = 1 to ≁
      i = 1;
      j = 2;
      do cage = 1 to 6;
         u_cage = sqrt(&var_cage) * rannor(&seed);
         do condition = i,j;
            e1 = rannor(&seed);
            e2 = rannor(&seed);
            e3 = rannor(&seed);
            diet = "normal  ";
            noise = u_cage + 5.0318983*e1;
            output;
            diet = "restrict";
            noise = u_cage - 1.240089*e1 + 4.8766977*e2;
            output;
            diet = "suppleme";
            noise = u_cage - 1.240089*e1 - 1.594895*e2 + 4.6085237*e3;
            output;	
         end;
         j = j + 1;
         if (j = 5) then do;
            i = i + 1;
            j = i + 1;
         end;
      end;
   end;
run;

PROC PRINT DATA=noise;
TITLE 'The output from noise simulation';
run;