set more off
capture log close
clear

// ***************** paths *******************************
local misenda = "/home/ricmora/AAOFICIN/CURSOS/MADE/materiales/slides"

set seed 27031967
set obs 5000
matrix mu = (0,0)
matrix sigma = (1,0.5\0.5,1)

//  ************ variables creation ***********************
drawnorm eps_m eps_h, means(mu) cov(sigma)	// errors
gen kids = uniform()<0.6
sum kids
gen u=uniform()
gen educ=8*(u<=.15)+12*(u>.15)*(u<=.6)+16*(u>.6)*(u<=.9)+21*(u>.9)
tab educ
gen U_m = 0.3 + 0.05* educ + 0.5*kids + eps_m
gen U_h = 0.8 - 0.02* educ + 2*kids + eps_h
gen work=U_m>=U_h
gen educkids=100*kids+educ
bysort educkids: sum work



// estimation
probit work educ kids
mfx
dprobit work educ kids

// predicting probabilities
predict p_hat,p
gen work_hat=p_hat>.5
tab work work_hat

// testing
test educ kids
test (educ=0.07) (kids=-1.5)
test 2*educ = -kids


// marginal effects of having a kid
gen kids_old=kids
replace kids=1
predict p1,p
replace kids=0
predict p0,p
gen Mg_kid = p1 - p0
bysort educ: sum Mg_kid
replace kids=kids_old
drop kids_old

// marginal effects of one extra year of education (individual calculus approximation)
predict xb_hat,xb
gen Mg_educ_cal=normalden(xb_hat)*_b[educ]	    // this is the individual's marginal effect
su Mg_educ_cal

// marginal effects of one extra year of education  (at the mean values)
rename educ educ_old
rename kids kids_old
egen educ = mean(educ_old) if e(sample)
egen kids = mean(kids_old) if e(sample)
predict xb_hat_avg,xb
gen Mg_educ_avg=normalden(xb_hat_avg)*_b[educ]	    // this is the marginal effect on averages
su Mg_educ_avg
replace educ=educ_old
replace kids=kids_old
 
// logit
logit work educ kids
mfx
exit
 // marginal effects of one extra year of education
 genr beta_logit=$coeff
 series educ0=educ
 matrix x0={const,educ0,kids}
 series educ1=educ+1
 matrix x1={const,educ1,kids}
 series x1b = x1*beta_logit
 series x0b = x0*beta_logit
 series Mg_educ_logit = cdf(N,x1b)-cdf(N,x0b)
 summary Mg_educ_logit --by=kids --simple 

quit
 

exit


// probit example
capture program drop prosim
program prosim, rclass
	version 10.0
	syntax [, obs(integer 1) gamma0(real 0) gamma1(real 0) meanZ(real 0) stddevZ(real 1) stddev(real 1)]
	drop _all
	set obs `obs'
	generate Z = invnormal(uniform())*`stddevZ'+`meanZ'
	generate u = invnormal(uniform())*`stddev'
	generate y = `gamma0'+`gamma1'* Z+ u >=0
	probit y Z
	return scalar t1 = _coef[Z]/_se[Z]
end
set seed 0
simulate t1 = r(t1), reps(1000): prosim, obs(30) gamma0(1) gamma1(0.5) stddevZ(4)
sum t1
histogram t1, kdensity plot(function stdnorm = normalden(x,0,1), ra(-4 4) lpattern(dash))



// This generates the dataset for the Probit Example in slides