/*
** Truncation MLE Estimation
*/

new ;

@ LOADING DATA @

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

@ OPEN OUTPUT FILE @

    output file = turnc.out reset ;

@ DEFINE VARIABLES @

	dat     = delif(dat,dat[.,2] .== 0) ; @ Delete nonemployed people @ 
    nlf     = dat[.,1]  ; @ dummy variable for non-labor-force  @
    emp     = dat[.,2]  ; @ dummy varrable 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-manufactoring 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 = uspell ;
	vny = {"uspell"};

	@ Regressors @

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

/*
**  DO NOT CHANGE FROM HERE
*/

	vnx1 = {"sigma"};
	vnx = vnx|vnx1 ;

@	Deleting observations with y = 0  @

	zz = yy~xx ;
	zz = delif(zz,zz[.,1] .== 0);
	yy = zz[.,1] ;
	xx = zz[.,2:cols(zz)];

@   DEFINE INITIAL VALUE @

    bb = invpd(xx'xx)*(xx'yy) ;
	ss = sqrt((yy-xx*bb)'(yy-xx*bb)/n) ;
	bb = bb|ss ;

library optmum;
#include optmum.ext;
optset ;

proc  f(b) ;
local ltt, q ;

	q   = rows(b);
	ltt = (-.5*ln(2*pi) -.5*ln(b[q]^2) -(.5/b[q]^2).*(yy-xx*b[1:q-1])^2 )
         - ln( cdfn( xx*b[1:q-1]/abs(b[q]) ) ) ; 

retp(  -sumc(ltt)  ) ;
endp ;

b0 = bb ;
__title = "Truncation MLE ";
_opgtol = 1e-6;
_opstmth = "bfgs,half";
__output = 0 ;

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

@ value of likelihood function @

	logl = -func ;

@ Covariance matrix by Hessian @

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

@ Covariance matrix by BHHH @

	proc  ff(b) ;
	local ltt, q;

	q   = rows(b);
	ltt = (-.5*ln(2*pi) -.5*ln(b[q]^2) -(.5/b[q]^2).*(yy-xx*b[1:q-1])^2 )
         - ln( cdfn( xx*b[1:q-1]/abs(b[q]) ) ) ; 

	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 ;
"" ;
"Truncation MLE Result" ;
"------------------------" ;
" dependent variable: " $vny ;
"" ;
" # of obs:           " rows(yy);
"";
" log likelihood:     " logl ;
"" ;
" 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 ;