#delimit;
clear;
set more off;
set mem 20m;
set matsize 200;
cap log close;
log using c:\stataproject\micrometrics\hw8\hw8.log,replace;
use C:\stataproject\micrometrics\hw8\labsup.dta,clear;

des;

*1.1;
reg worked morekids edu age nonmomi,robust;

*1.2;
ivreg worked (morekids=samesex) edu age nonmomi,robust;

*1.3;
probit worked morekids edu age nonmomi;
dprobit worked morekids edu age nonmomi;

*1.4;
reg morekids samesex edu age nonmomi;
predict vhat,residual;
probit worked morekids edu age nonmomi vhat;

*1.5;
*Before implementing the IV probit estimation, see whether own written probit works fine.;
*maximum likelihood estimation;
cap program drop myprobit;
program define myprobit;
	version 7;
	args lnf theta;
	quietly replace `lnf' = 
	ln(normprob(`theta')) if $ML_y1==1;
	quietly replace `lnf' = 
	ln(1-normprob(`theta')) if $ML_y1==0;
end;

ml model lf myprobit (worked=morekids edu age nonmomi);
ml maximize;
probit worked morekids edu age nonmomi;
*The results are exactly same;

*Implementing IV probit estimation;
*Since treating y2 as binary variable is not trivial, I assume y2 is conditionally normally distributed.;
*Note that this assumption cannot be true when y2 is binary. Thus this is a short-cut and dirty answer;
*Model;
* y1 = 1(x delta + alpha1 y2 + e > 0);
* y2 = z theta + v;
* (e,v)|x ~ N(Ee=0,Ev=0,var(e)=1,var(v)=tau2^2,corr(e,v)=rho1);
* x=[edu age nonmomi], z=[samesex edu age nonmomi];
* x delta is called as `delta' in the following estimation.;
* z theta is called as `theta' in the following estimation.;
*See Wooldridge's "Econometric Analysis of Cross Section and Panel Data", p476, equation (15.50) for likelihood function.;

*maximum likelihood estimation;
cap program drop ivprobit;
program define ivprobit;
	version 7;
	args lnf delta theta alpha1 tau2 rho1;

	quietly replace `lnf' = 
	ln(normprob((`delta'+`alpha1'*morekids+
	(`rho1'/`tau2')*(morekids -`theta'))/(1-`rho1'^2)^(1/2))) 
	-ln(`tau2') - (1/2)*(morekids-`theta')^2/`tau2'^2 
	if $ML_y1==1;

	quietly replace `lnf' = 
	ln(1-normprob((`delta'+`alpha1'*morekids+
	(`rho1'/`tau2')*(morekids -`theta'))/(1-`rho1'^2)^(1/2))) 
	-ln(`tau2') - (1/2)*(morekids-`theta')^2/`tau2'^2 
	if $ML_y1==0;
end;

ml model lf ivprobit (y1:worked = edu age nonmomi) (y2:samesex edu age nonmomi) (alpha1:) (tau2:) (rho1:);
ml maximize;

*confirmation of the correctness of the second equation;
reg morekids samesex edu age nonmomi;

log close;
exit;