/* discrete choice example
   nonlinear least squares and maximum likelihood  */

/* this example uses the data set ecn1_wag from assignment 3 */

closeall;
output file = dchoice.out reset;

/* read data */

open f201 = ecn1_wag for read;
n = rowsf(f201);
z = readr(f201,n);

/* dependent variable is whether earns wage greater than $15 an hour */
highwage = (z[.,1] .>= 15);

/*create x matrix */
x = ones(n,1) ~ z[.,2 3 5 6];  @constant, age, male, hs, and college@

/* linear probability model */

beta_ols = invpd(x'x)*x'highwage;
print "OLS point estimates";
print beta_ols;

/* for curiousity, how many predicted probabilities NOT in the unit interval?  */
pred_p = x*beta_ols;
print "LPM predictions not in unit interval";
print (1-meanc(pred_p .> 0 .and pred_p .< 1));

/* estimate conditional expectation function cap_phi(x*beta)  */

library optmum,maxlik;
optset;

proc nlin_ce(beta);
local phi;
phi = cdfn(x*beta);
retp(  sumc( (highwage - phi)^2 ./ sqrt(phi.*(1-phi)) ));
endp;


beta0= zeros(5,1);
optprt(optmum(&nlin_ce,beta0));


/* probit estimates */

proc probit_llf(beta,z);
retp( highwage .* ln(cdfn(x*beta)) + (1-highwage) .* ln(cdfnc(x*beta)) );
endp;

print "************* PROBIT ESTIMATES *********************/";

maxset;
maxprt(maxlik(x,0,&probit_llf,beta0));

/* logit estimates */

print "************** LOGIT ESTIMATES ********************/";

proc logit_llf(beta,z);
retp( highwage .* ln(exp(x*beta)./(1+exp(x*beta)))
      + (1-highwage) .* ln(1 / (1+exp(x*beta))) );
endp;

maxset;
maxprt(maxlik(x,0,&logit_llf,beta0));



end;