/*-------------------------------------------------------------------------
                                                                         
  OLS.PRC                                                                
  Purpose: Computes Basic OLS Estimates                                  
  Written by Hyeongwoo Kim (Jan 31, 2003)
  Modified on May 18, 2010
                                
-------------------------------------------------------------------------
                                                                         
  Format: {b,rsd,vcb,se,tr,pv,rsq,arq} = ols(y,x,c);
                                                                         
  Input : y    (nX1) Vector of Dependent variables                        
          x    (nXk) Matrix of Independent Variables (w/o constant)      
          c=0  No Constant                                                
           =1  Include Constant                                          
           =2  Include Constant and Linear Time Trend                    
                                                                         
  Output: b    (kx1) OLS Coefficient Estimates                           
          rsd  (nx1) OLS Residuals                                       
          vcb  (kxk) Variance-Covariance Matrix of OLS Coefficients      
          se   (kx1) Standard Errors of OLS Coefficients                 
          tr   (kx1) t-ratios                                            
          pv   (kx1) p-values                                            
          rsq  (1x1) R-square   
          art  (1x1) Adjusted R-square                                         
                                                                         
-------------------------------------------------------------------------*/


proc(8) = ols(y,x,c);
    local n,k,df,b,rsd,ssq,vcb,se,tr,pv,rsq,arq;
    n   = rows(x);
    if c eq 1; x = ones(rows(x),1)~x; endif;
    if c eq 2; x = ones(rows(x),1)~seqa(1,1,rows(x))~x; endif;
    k   = cols(x);
    df  = n - k;
    b   = invpd(x'x)*x'y;			
    rsd = y - x*b;		
    ssq = (rsd'rsd)/df;	
    vcb = ssq*inv(x'x);
    se  = sqrt(diag(vcb));
    tr  = b./se;
    pv  = 2*cdftc(abs(tr),df);
    rsq = 1 - ((rsd'rsd)/((y-meanc(y))'(y-meanc(y))));
    arq = 1 - (1-rsq)*(n-1)/(n-1-k);
    retp(b,rsd,vcb,se,tr,pv,rsq,arq);
endp;