set seed 547
# This generates the dataset for the Probit Example in slides


# ***************** paths *******************************
#string misenda = "/media/Elements/AAOFICIN/CURSOS/MICCUA/materiales/Slides18_ProbitEstimation_QM"


# *************** Before the loop *********************
nulldata 5000
matrix SIG = {1,0.5;0.5,1}
matrix A = cholesky(SIG)

#****** open a loop, to be repeated R=500 times **********
#loop 5 --progressive --quiet

 #************* variables creation ***********************
 genr kids = uniform()<0.6
 summary kids
 genr u=uniform()
 genr educ=8*(u<=.15)+12*(u>.15)*(u<=.6)+16*(u>.6)*(u<=.9)+21*(u>.9)
 genr u1 = normal()
 genr u2 = normal()
 genr eps_m = A[1,1]*u1+A[1,2]*u2
 genr eps_h = A[2,1]*u1+A[2,2]*u2
 genr U_m = 0.3 + 0.05* educ + 0.5*kids + eps_m
 genr U_h = 0.8 - 0.02* educ + 2*kids + eps_h
 genr work=U_m>=U_h
 genr educkids=100*kids+educ
 
 summary work --by=educkids --simple 

 # estimation, testing, and fit
 probit work const educ kids
 omit educ kids --wald
 
 probit work const
 add educ kids --quiet
 
 
 outfile --write null
 # estimating unrestricted model and storing loglikelihood
 probit work const educ kids --quiet
 scalar lur= $lnl

 # estimating restricted model and storing loglikelihood
 genr x=educ-2*kids
 probit work const x --quiet
 scalar lr= $lnl

 # computing the LR statistic and p-value
 scalar LR=2*(lur-lr)
 scalar pval = pvalue(X, 1, LR)

 # printout
 outfile --close
 printf "\nLikelihood Ratio test\nH0: 2*beta_educ+beta_kids=0\nLR %.8g   p-value %.8g,\n", LR, pval

 # back to the unrestricted model
 probit work const educ kids --quiet
 genr yhat=$yhat
 summary work yhat --by=educkids --simple
 
 # marginal effects of having a kid
 genr beta=$coeff
 series kids0=0
 matrix x0={const,educ,kids0}
 series kids1=1
 matrix x1={const,educ,kids1}
 series x1b = x1*beta
 series x0b = x0*beta
 series Mg_kid = cnorm(x1b)-cdf(N,x0b)
 summary Mg_kid --by=educ --simple 
 summary Mg_kid --simple
 
 # marginal effects of one extra year of education
 genr beta=$coeff
 series educ0=educ
 matrix x0={const,educ0,kids}
 series educ1=educ+1
 matrix x1={const,educ1,kids}
 series x1b = x1*beta
 series x0b = x0*beta
 series Mg_educ = cdf(N,x1b)-cdf(N,x0b)
 summary Mg_educ --by=kids --simple 
 summary Mg_educ --simple


 # marginal effects of one extra year of education: calculus approximation 
 genr beta=$coeff
 matrix X={const,educ,kids} 
 series Xb=X*beta

 genr meanXb=mean(Xb)
 series Mg_educ_slope=pdf(N,meanXb)*$coeff(educ)	# this is the slope in gretl output
 series Mg_educ_cal=pdf(N,Xb)*$coeff(educ)	    # this is the individual's marginal effect
 summary Mg_educ_slope Mg_educ_cal --by=kids --simple
 

 # linear probability model
 ols work const educ kids --robust
 genr yhat_linear=$yhat
 summary work yhat yhat_linear --by=educkids --simple
 
 
 # logit
 logit work const educ kids
 # 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