/* 

   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 a; 
input X @@; 
  do rep= 1 to 4000; 
        
           Y = &A + &B*X + &sig*normal(1827689); 
		   file print; if rep=1 then put Y= X= ; 

   output; end; 
cards; 
14 9 20 25 13 8 12 10 10 5 30 26 18 11 19 22 25 19 8 15
proc sort data=a; by rep; 
proc print data=a; where rep=1; 
  title "Data generated as Y = &a + &b.X + &sig e "; ** try removing the dot **; 

proc reg noprint outest=betas; 
  model Y = X; by rep; 

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 intercept X _RMSE_; 
proc chart data=betas ; vbar Intercept/ levels=40; 
/*
proc g3d; scatter Intercept*_RMSE_=X/noneedle shape="balloon" size=.5 color="black";
*/ 
run; 

data notindep; 
do i=1 to 1000; Z=normal(1827655); Chisq=Z*Z; 
output; end;  title "Uncorrelated but not Independent";
proc gchart; vbar Z Chisq; 
proc corr; var Z Chisq; 
proc gplot; plot Chisq*Z; 
symbol1 v=dot h=.5; 
run; quit;

/*------------------------------------
   PC results:
     4000 regressions in 3.11 seconds (3.06 seconds CPU time)
-------------------------------------*/