/*
** Probit MLE Estimation
*/

new ;

@ LOADING DATA @

    load dat[1962,25] = mwpsid82.db ; @ Data on married women from PSID 82 @

@ OPEN OUTPUT FILE @

    output file = probit.out reset ;

@ DEFINE VARIABLES @

    nlf     = dat[.,1]  ; @ dummy variable for non-labor-force  @
    emp     = dat[.,2]  ; @ dummy variable for employed workers @
    wrate   = dat[.,3]  ; @ hourly wage rate @ 
    lrate   = dat[.,4]  ; @ log of (hourly wage rate ($)+1) @
    edu     = dat[.,5]  ; @ years of schooling @
    urb     = dat[.,6]  ; @ dummy variable for urban resident @
    minor   = dat[.,7]  ; @ dummy variable for minority race @
    age     = dat[.,8]  ; @ age @
    tenure  = dat[.,9]  ; @ # of months under current employer @
    expp    = dat[.,10] ; @ years of experience since age 18 @
    regs    = dat[.,11] ; @ dummy variable for South @
    occw    = dat[.,12] ; @ dummy variable for white collar @
    occb    = dat[.,13] ; @ dummy variable for blue collar @
    indumg  = dat[.,14] ; @ dummy variable for manufacturing industry @
    indumn  = dat[.,15] ; @ dummy variable for non-manufacturing industry @
    unionn  = dat[.,16] ; @ dummy variable for union membership @
    unempr  = dat[.,17] ; @ local unemployment rate @
    ofinc   = dat[.,18] ; @ other family income in 1980 ($) @
	lofinc  = dat[.,19] ; @ log of (other family income + 1) @
	kids    = dat[.,20] ; @ # of children of age <= 17 @
	wune80  = dat[.,21] ; @ UNKNOWN @
	hwork   = dat[.,22] ; @ hours of house work per week @
	uspell  = dat[.,23] ; @ # of unemployed weeks in 1980 @
	search  = dat[.,24] ; @ # of weeks looking for job in 1980 @
	kids5   = dat[.,25] ; @ # of children of age <= 5 @ 


@ DEFINE # of OBSERVATIONS @

	 n = rows(dat) ;

@  DEFINE DEPENDENT VARIABLE AND REGRESSORS @

	@ Dependent variable @

	yy = emp ;
	vny = {"emp"};

	@ Regressors @

	xx = ones(n,1)~edu~urb~minor~age~expp~kids~kids5~lofinc ;
	vnx = {"cons", "edu", "urb", "minor", "age", "expp",
           "kids", "kids5", "lofinc" };

/*
**  DO NOT CHANGE FROM HERE
*/

@   DEFINE INITIAL VALUE @

    bb = invpd(xx'xx)*(xx'yy) ;

library optmum;
#include optmum.ext;
optset ;

proc  f(b) ;
local ltt   ;

	ltt = yy.*ln( cdfn(xx*b) ) + (1-yy).*ln( 1-cdfn(xx*b) ) ; 

RETP(  -sumc(ltt)  ) ;
ENDP  ;

b0 = bb ;
__title = "Probit MLE ";
_opgtol = 1e-4;
_opstmth = "newton, helf";
__output = 0 ;

{b,func,grad,retcode} = optmum(&f,b0) ;

@ value of likelihood function @

	logl = -func ;

@ likelihood ratio test @

	rlogl = sumc(yy)*ln(sumc(yy)/n) + (n-sumc(yy))*ln(1-sumc(yy)/n) ;
	lrt   = 2*(logl-rlogl);
	lrdf  = cols(xx)-1 ;
	
@ Covariance matrix by Hessian @

	covhess = invpd(hessp(&f,b));

@ Covariance matrix by BHHH @

	proc  ff(b) ;
	local ltt   ;

		ltt = yy.*ln(cdfn(xx*b)) + (1-yy).*ln(1-cdfn(xx*b));

	RETP(  ltt  ) ;
	ENDP  ;

	gtt = gradp(&ff,b);
	covbhhh = invpd(gtt'gtt) ;

@ Robust covariance matrix @

	covrob = covhess*invpd(covbhhh)*covhess ;

@ Computing Standard errors @

	sehess = sqrt(diag(covhess));
	sebhhh = sqrt(diag(covbhhh));
	serob  = sqrt(diag(covrob)) ;

stathess = b~sehess~(b./sehess);
statbhhh = b~sebhhh~(b./sebhhh);
statrob  = b~serob~(b./serob);
yyhess = vnx~stathess;
yybhhh = vnx~statbhhh;
yyrob  = vnx~statrob;

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

format /rd 10,4 ;
"" ;
"Probit Estimation Result" ;
"------------------------" ;
" dependent variable: " $vny ;
"" ;
" log likelihood:     " logl ;
"" ;
" Pseudo R-square:    " (1-logl/rlogl) ;
"" ;
" LR test, df, p-val: " lrt lrdf cdfchic(lrt,lrdf) ;
"" ;
" Using Hessian ";
"" ;
"variable     coeff.  std. err.   t-st " ;
yyprin = printfm(yyhess,mask,fmt); 
"" ;
" Using BHHH ";
"" ;
"variable     coeff.  std. err.   t-st " ;
yyprin = printfm(yybhhh,mask,fmt); 
"" ;
" Using Robust ";
"" ;
"variable     coeff.  std. err.   t-st " ;
yyprin = printfm(yyrob,mask,fmt);

output off ;