new;
// PROBIT MODEL USING MLE
// DATA

x=rndn(100,1);
z=rndu(100,1);
ys=x+z+rndn(100,1);
y=ys.>0;


// Maximization of the Log-Likelihood 

b0={0.5, 0.5}; // starting values for the coefficients

{b,f,g,retcode} = qnewton(&llik,b0);

// Var-Cov Matrix via Hessian
vcov_Hes = invpd(hessp(&llik,b)); // this uses GAUSS's own numerical hessian procedure 
se_Hes= sqrt(diag(vcov_Hes));    // to calcuate the matrix of the second derivatives 
                                 // it is not as accurate as Hessian computed analytically

// Var-Cov Matrix via Outer Product of Gradient

Gmat = glik(b); // this uses the first derivative procedure we have below

vcov_OPG=(invpd(Gmat'Gmat));
se_OPG= sqrt(diag(vcov_OPG));


//OUTPUT

titles = "Coeff."~"SERR_Hes"~"SERR_OPG";

// cls;
"MLE ESTIMATION OF THE PROBIT MODEL";
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~";
$titles;
b~se_Hes~se_OPG;
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~";
"Log-Likelihood value";;-f;
"Note: SERR_Hes and SERR_OPG are Standard Errors Computed";
"using Empirial Hessian estimator and OPG estimator respectively";


// PROCEDURES

// Negative of the Log-Likelihood
proc llik(beta);
	local lik;
lik = ln(cdfn(x*beta[1]+z*beta[2])).*y + ln(1-cdfn((x*beta[1]+z*beta[2]))).*(1-y);

retp(-sumc(lik));
endp;

// Gradient matrix
proc glik(beta);
	local g;
g = ((x~z).*pdfn(x*beta[1]+z*beta[2])./cdfn(x*beta[1]+z*beta[2])).*y 
    + (-(x~z).*pdfn(x*beta[1]+z*beta[2])./cdfn(x*beta[1]+z*beta[2])).*(1-y);
	
	retp((g));
endp;

/* Notes :
pdfn() is the normal density function
cdfn() is the normal cumulative distribution function
*/