/* GET_DAT.G: return observed (y,X) matrices for SUR model
       	     John Rust, Yale University, April 1997 */


     proc (5)=get_dat;

         local b0,r0,x,y,i,n,p,k;

	 /* Code below generates (y,x) for artificial data
	    where we know the true parameters generating the 
	    true beta and true omega, denoted b0 and r0 respectively
	    below 

         let b0=.5 1 -.5 .8;
         let r0[2,2]=1 .2 .2 1;
	 trueomeg=invpd(r0);
	 r0=chol(r0);

	 n=1000;
	 p=2;
	 k=4;
	 x=zeros(p*n,k);
	 y=zeros(p*n,1);

	 i=1; do until i > n;

	   x[p*(i-1)+1,1:2]=1~rndn(1,1);
	   x[p*(i-1)+2,3:4]=1~2*rndn(1,1);

           y[p*(i-1)+1:p*(i-1)+2,.]=x[p*(i-1)+1:p*(i-1)+2,.]*b0+r0*rndn(p,1);

         i=i+1; endo; */


	 /* Code below loads p securities from STOCKDAT data set
	    and computes (y,x) for the CAPM problem in Question 1 of
	    the midterm.  */

	 local f1,nms,mkt,ty;

         open f1=stockdat;
	 n=rowsf(f1);
	 x=readr(f1,n);
	 nms=getname("stockdat");
	 mkt=x[.,loc("VALAV",nms)]; /* load the value-weighted average daily
				       returns as the market regressor */

	 "Loading STOCKDAT data set for CAPM model using VALAV and";
	 "ARCO, REXNORD and USR securities ";

	 ty=x[.,loc("ARCO",nms)~loc("REXNORD",nms)~loc("USR",nms)];

         p=3;
	 k=6;
	 x=zeros(p*n,k);
	 y=zeros(p*n,1);

	 i=1; do until i > n;

	   j=1; do until j > p;
   	      x[p*(i-1)+j,1+2*(j-1):2+2*(j-1)]=1~mkt[i];
              y[p*(i-1)+j,.]=ty[i,j];
           j=j+1; endo;

         i=i+1; endo;  

	 /* Code below loads ALL 64 securities from STOCKDAT data set
	    and computes (y,x) for the CAPM problem in Question 1 of
	    the midterm. 

	 local f1,nms,mkt,ty;

         open f1=stockdat;
	 n=rowsf(f1);
	 x=readr(f1,n);
	 nms=getname("stockdat");
	 mkt=x[.,loc("VALAV",nms)]; /* load the value-weighted average daily
				       returns as the market regressor */

	 "Loading STOCKDAT data set for CAPM model using VALAV and";
	 "ALL 63 remaining securities ";

	 ty=x[.,1:22]~x[.,23:63]~x[.,65];

         p=cols(ty);
	 k=2*p;
	 x=zeros(p*n,k);
	 y=zeros(p*n,1);

	 i=1; do until i > n;

	   j=1; do until j > p;
   	      x[p*(i-1)+j,1+2*(j-1):2+2*(j-1)]=1~mkt[i];
              y[p*(i-1)+j,.]=ty[i,j];
           j=j+1; endo;

         i=i+1; endo;  */

     retp(x,y,n,p,k);

     endp;