/*
**	OLS Program Using Proc
*/

@	Data loading   @

@ LOADING DATA @

    load dat[500,17] = wage22.txt ; 

@  Selecting observations @

   @ dat = selif(dat,dat[.,11] .== 1) ; @
     dat = delif(dat,dat[.,11] .== 1) ; 

@ open OUTPUT FILE @

    output file = new_ols.out reset ;

@ DEFINE VARIABLES @

    wage    = dat[.,1]  ; 
    hours   = dat[.,2]  ; 
    iq      = dat[.,3]  ;
    kww     = dat[.,4]  ; 
    educ    = dat[.,5]  ; 
    exper   = dat[.,6]  ; 
    tenure  = dat[.,7]  ; 
    age     = dat[.,8]  ; 
    married = dat[.,9]  ; 
    black   = dat[.,10] ; 
    south   = dat[.,11] ; 
    urban   = dat[.,12] ; 
    sibs    = dat[.,13] ; 
    brthord = dat[.,14] ; 
    meduc   = dat[.,15] ; 
    feduc   = dat[.,16] ; 
    lwage   = dat[.,17] ; 

@ Dependent variable @

	yy = lwage ;
	vny = {"lwage"};

@ Regressors @

    xx = ones(rows(yy),1)~educ~exper~(exper^2)  ; 
	vnx = {"cons", "edu", "exper", "exper^2" } ;

@   Do not change below  @
@	OLS using yy and xx  @

proc (6) = reg(a,b) ; @ a=dep. var; b=regressors; 6 = # of outcomes @
local c, e, s2, v, t, k, cc, dd, ts, pv, sst, sse, r2 ;

t  = rows(a);
k  = cols(b);
c  = invpd(b'b)*(b'a);
e  = a - b*c ;
s2 = (e'e)/(t-k);
v  = s2*invpd(b'b);
ts = (c./sqrt(diag(v)))    ;
pv = 2*cdftc(abs(ts),t-k)  ;
cc = c~sqrt(diag(v))~ts~pv ;
sst   = a'a - t*meanc(a)^2 ;
sse   = e'e;
r2    = 1 - sse/sst;

retp(t,k,c,v,r2,cc);
endp;

{t,k,c,v,r2,cc} = reg(yy,xx);
@ t = # of observations @
@ k = # of regressors   @
@ c = vector of estimated coefficients  @
@ v = covariance matrix @
@ r2 = R-squared @
@ cc = matrix of c, se, t-st, and p-val @

econ = vnx~cc;


@ F-Test for H_o: b_3 = b_4 @

    lr = zeros(2,k);
        lr[1,3] = 1;
        lr[2,4] = 1;
    sr = zeros(2,1);

    ft   = (lr*c-sr)'invpd(lr*v*lr')*(lr*c-sr)/2;
    p_ft = cdffc(ft,rows(lr),t-k);
    

@ Printing out OLS results @

output file = ols.out reset;

let mask[1,5] = 0 1 1 1 1;
let fmt[5,3] =
    "-*.*s" 8 8
    "*.*lf" 10 4
    "*.*lf" 10 4
    "*.*lf" 10 4
    "*.*lf" 10 4;

format /rd 11,3 ;
"" ;
"OLS Regression Result" ;
"------------------------" ;
" dependent variable: " $vny ;
"" ;
" R-Squares         : " r2 ;
"" ;
" # of Observations : " t  ;
"" ;
"variable    coeff.   std. err.   t-st      p-val" ;
yyprin = printfm(econ,mask,fmt); 
"" ;
"F-statistic for H_o: b_3 = b_4  and p-val";
"            " ft p_ft;
""; 

output off ;