/* 

   Goal: to simulate 4000 simple linear regressions and study the joint
limiting distributions of a, b, MSE.                           */

* -------------------;
Options notes;        
/*proc datasets;        ;
delete allbetas run;  ;
---------------------*/
%let a=10; %let b=2; %let sig=3; 
Data betas; array XX(20); 
input Xx1-Xx20; 
  do rep= 1 to 4000;
    SSX=0; SSY=0; SXY=0; SX=0; SY=0; 
    do i=1 to 20;  X=XX(i);
           Y = &A + &B*X + &sig*normal(1827689); 
               SSX=SSX+X*X; SSY=SSY+Y*Y; SXY=SXY+X*Y; 
                   SX=SX+X; SY=SY+Y; 
                   file print; if rep=1 then put Y= X= ; 
               if rep=1 then put SSX= SSY= SXY= SX=SY=;
        end;  b=(SXY-SX*SY/20)/(SSX-SX*SX/20);
              a = (SY-b*SX)/20;
                  MSE= (SSY-SY*SY/20)- b*(SXY-SX*SY/20); 
                  RMSE=sqrt(MSE/18); keep a b RMSE;
    output; end; 
cards; 
14 9 20 25 13 8 12 10 10 5 30 26 18 11 19 22 25 19 8 15
 ;
proc print data=betas(obs=10); title "BETAS";

* Do this the first time; 
data allbetas; set betas;

/* Then do this; 
proc append base=allbetas data=betas; 
  ** Dataset allbetas must exist and have right variables **; 
*/
proc print data=allbetas(obs=10);  title "ALLBETAS"; run;

proc corr data=allbetas; var a b RMSE; 
proc chart data=betas ; vbar a / levels=40; 
run; 

/*----------------------------------------------------------
  PC results: 
     40,000 regressions in  1.05 seconds (1.05 CPU seconds)
    140,000 regressions in  3.61 seconds (3.59 CPU seconds)
    250,000 regressions in  6.58 seconds (6.39 )
  1,000,000 regressions in 25.57 seconds (25.40 )
------------------------------------------------------------*/