/*
** Tobit MLE Estimation
*/

new ;

@ LOADING DATA @

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

@ OPEN OUTPUT FILE @

    output file = tobit.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
*/

@   DEFINE INITIAL VALUE @

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

@	dummy variable for yy @

	v  = 0 ;
	yd = dummy(yy,v);
	yd = yd[.,2] ;	

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

library optmum;
#include optmum.ext;
optset ;

proc  f(b) ;
local ltt, q ;

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

retp(  -sumc(ltt)  ) ;
endp ;

proc gradien(b) ;
local q, i, gge, g1, g2, gap, gac, poi,e ;

	q  = rows(b)    ;
	gge= zeros(q,1) ;

i = 1; do while i <= n;

	gap  = pdfn(xx[i,.]*b[1:q-1]/b[q]) ;
	gac  = cdfn(xx[i,.]*b[1:q-1]/b[q]) ;
	poi  = xx[i,.]*b[1:q-1]/b[q]       ;
	e    = yy[i] -xx[i,.]*b[1:q-1]     ;
	g1 = yd[i]*e/b[q]^2*xx[i,.]'
              + (1-yd[i])*( -gap/b[q]/(1-gac) )*xx[i,.]' ;
	
	g2 = yd[i]*( -1/b[q] + 1/b[q]^3*e^2  )
			  + (1-yd[i])*( poi*gap/b[q]/(1-gac) ) ;
	
	gge = gge + (g1|g2) ;

i = i + 1 ; endo ;

retp( -gge' ) ; endp ;


proc hessi(b) ;
local q, i, he, h11, h12, h21, h22, gap, gac, poi ;
	q  = rows(b)    ;
	he = zeros(q,q) ;

i = 1; do while i <= n;
	gap  = pdfn(xx[i,.]*b[1:q-1]/b[q]) ;
	gac  = cdfn(xx[i,.]*b[1:q-1]/b[q]) ;
	poi  = xx[i,.]*b[1:q-1]/b[q]       ;

	h11 = yd[i]*(-1/b[q]^2)*xx[i,.]'xx[i,.]
              + (1-yd[i])*( poi*gap/(1-gac)/b[q]^2
                            - gap^2/(1-gac)^2/b[q]^2 )*xx[i,.]'xx[i,.] ;
	
	h12 = yd[i]*( -2/b[q]^3*(yy[i]-xx[i,.]*b[1:q-1])*xx[i,.]' )
			  + (1-yd[i])*( (-poi^2*gap+gap)/b[q]^2/(1-gac)
                            + gap^2*poi/b[q]^2/(1-gac)^2 )*xx[i,.]'   ;
	h21 = h12' ;
	h22 = yd[i]*( 1/b[q]^2 - 3/b[q]^4*(yy[i]-xx[i,.]*b[1:q-1])^2 )
			  + (1-yd[i])*( (-2*poi*gap+poi^3*gap)/b[q]^2/(1-gac)
							- poi^2*gac^2/b[q]^2/(1-gac)^2 ) ;
	
	he = he + ((h11~h12)|(h21~h22)) ;

i = i + 1 ; endo ;

retp( -he ) ; endp ;

b0 = bb ;
__title = "Tobit MLE ";
_opgtol = 1e-6;
_opstmth = "bfgs,half";
__output = 0 ;
_opgdprc = &gradien; 
_ophsprc = &hessi ;  

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

@ value of likelihood function @

	logl = -func ;

@ Covariance matrix by Hessian @

covhess = invpd(hessi(b)) ; 

@ Covariance matrix by BHHH @

	proc  ff(b) ;
	local ltt, q;

	q   = rows(b);
	ltt = yd.*( -.5*ln(2*pi) -.5*ln(b[q]^2) -(.5/b[q]^2).*(yy-xx*b[1:q-1])^2 )
         + (1-yd).*ln( 1- cdfn( xx*b[1:q-1])/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 ;
"" ;
"Tobit Estimation Result" ;
"------------------------" ;
" dependent variable: " $vny ;
"" ;
" 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 ;