/* 

   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; %LET N=40000;
proc iml; 
X = J(20,1)||{14 9 20 25 13 8 12 10 10 5 30 26 18 11 19 22 25 19 8 15}`;
DAT=SHAPE(0,&N,3);
y = SHAPE(0,20,1); NAMES={A,B,RMSE};
DO II= 1 TO &N;
 
do I=1 to 20; 
    Y[I,1] = &A + &B*X[I,2] + &sig*normal(1827689); 
END;

b=INV(x`*x)*x`*y; sse= y`*y-B`*x`*y; 
DAT[II, ] = B`||SSE; 
END;
PART1=DAT[1:5, ]; PRINT PART1; 
CREATE BETAS FROM DAT[COLNAME=NAMES]; 
APPEND FROM DAT; 
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  9.35 seconds (9.34 CPU seconds)
------------------------------------------------------------*/