proc iml;

/***************************
 *                         *  
 *   Level Names           *
 *                         *
 *      INPUT              * 
 *    a -- nx1 Factor A    * 
 *                         * 
 *       RETURN            * 
 *  lev -- factor levels   * 
 *                         *  
 ***************************/ 
   start levnam(a);
   tmp=unique(a);
   if type(tmp) = 'N' then lev=trim(left(char(tmp)));
   else lev=trim(tmp);
   return(lev);
   finish levnam;
   

   
/***************************
 *                         *  
 *                         *    
 *  Interaction            *
 *                         * 
 *      INPUT              *
 *   a -- nx1 Factor A     *
 *   b -- nx1 Factor B     *
*                         * 
 *     RETURN              *
 *  ab -- nxab AB levels   *
 *                         * 
 ***************************/

   start inter(a,b);
   n=nrow(a);
   if type(a) = 'N' then loca=left(char(a));
   else loca=a;
   if type(b) = 'N' then locb=left(char(b));
   else locb=b;
   ab=concat(trim(loca),j(n,1," "),trim(locb));
   return(ab);
   finish inter;  

/***************************
 *                         *    
 *     Solution Module     *    
 *                         *    
 *        Input            *
 *    y -- data vector     * 
 *    x -- design Matrix   *   
 *                         *
 *       Output            *
 * beta -- Beta hat        *
 *  cov -- xpx ginv*sigma^2*   
 *  mse -- mse             *   
 *  ssm -- ssm             *   
 *  dfm -- dfm             *     
 *                         *
 ***************************/
   start solut(y,x,beta,cov,mse,ssm,dfm);
   n=nrow(y);
   xpxi=sweep(x`*x);
   beta=xpxi*x`*y;
   ssm=y`*x*beta;
   sse=y`*y-ssm;
   dfm=trace(xpxi*x`*x);
   dfe=n-dfm;
   mse=sse/dfe;
   cov=xpxi*mse;
   finish solut;

   store module=(inter solut levnam);
quit;