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


# *************** 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.4)
 summary kids
 genr u=uniform()
 genr educ=8*(u<=.1)+12*(u>.1)*(u<=.3)+16*(u>.3)*(u<=.6)+21*(u>.6)
 genr eps = normal()
 genr y_star = 0.07* educ -1.0*kids + eps
 genr work = (0.5<y_star)*(y_star<=2.5)+2*(2.5<y_star)
 
 freq kids
 freq educ
 freq work

 # estimation, testing, and fit
 probit work educ kids 

 # linear probability model gives bad results
 ols work const educ kids --robust

 # For Wald test for joint significance of all slopes, omit does not work
## ? probit work educ kids --quiet
## ? omit educ kids --wald
## No independent variables left after omissions

## Error executing script: halting
## > omit educ kids --wald

 # we can still compute the likelihood-ratio test
 # estimating unrestricted model and storing loglikelihood
 probit work educ kids --quiet
 scalar lur= $lnl

 # estimating restricted model and storing loglikelihood
 probit work const --quiet
 scalar lr= $lnl

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

 printf " LR = %.8g\n p-value = %.8g\n", LR,pval

 # marginal effects of having an additional kid
 probit work educ kids --quiet
 genr beta=$coeff[1:2]
 genr alpha=$coeff[3:4]
 series kids0=kids
 matrix x0={educ,kids0}
 series kids1=kids0+1
 matrix x1={educ,kids1}
 series x1b = x1*beta
 series x0b = x0*beta

 series Mg_kid0 = (cdf(N,x0b-alpha[1])-cdf(N,x1b-alpha[1]))
 series Mg_kid1 = (cdf(N,x1b-alpha[1])-cdf(N,x0b-alpha[1])) \
                 -(cdf(N,x1b-alpha[2])-cdf(N,x0b-alpha[2]))
 series Mg_kid2 = (cdf(N,x1b-alpha[2])-cdf(N,x0b-alpha[2]))

 summary Mg_kid* --simple

 # marginal effects of one extra year of education: calculus approximation 
 genr beta=$coeff[1:2]
 genr alpha=$coeff[3:4]
 matrix X={educ,kids} 
 series Xb=X*beta
 scalar meanXb=mean(Xb)
 scalar Mg_educ0=-pdf(N,meanXb-alpha[1])*$coeff(educ)
 scalar Mg_educ1= (pdf(N,meanXb-alpha[1])-pdf(N,meanXb-alpha[2]))*$coeff(educ)
 scalar Mg_educ2= pdf(N,meanXb-alpha[2])*$coeff(educ)
 
 printf " Mg_educ0 = %.8g\n Mg_educ1 = %.8g\n Mg_educ2 = %.8g\n", Mg_educ0, Mg_educ1, Mg_educ2

 # marginal effects on work
 scalar EMg_kid=mean(Mg_kid1)+2*mean(Mg_kid2)
 scalar EMg_educ=Mg_educ1+2*Mg_educ2

 printf " EMg_kid = %.8g\n EMg_educ = %.8g\n", EMg_kid, EMg_educ
 
quit