new;
 output file = ma_ml.out reset ;
library optmum;

N = 50;
let beta0 = .5;
beta1=3 ;
sigma = 1;
x=seqa(1,1,N)/N.*rndn(N,1) ; /* x is the regressor, could be anything, but if not centered around 0, 
the estimates will not be well identified. As for OLS you would like ``variation in the regressor" 
which here mean that we want to observe a good spread between 0s and 1s */

Nsim =3 ;
"Number of simulations:" Nsim ;
/* create a matrix to hold the estimated values */

results = zeros(Nsim,2) ;


/*starting values*/
let sv = 1  1 ;

estim=0 ;
do while estim < Nsim ;
estim=estim+1 ; 

u = sigma*rndn(N,1);

/*generate an a probit */
z0=beta0+beta1*x+u ;
z1=(z0.>0) ; 



/*__output=0;*/
{b_ml,f,g,retcode} = optmum(&logl,sv);
__output=1 ;
results[estim,.]=b_ml' ;
endo ;

/*results of last estimation*/
"estimated parameters:"
"beta0:" b_ml[1] ;
"true beta0 is:" beta0;
"beta1:" b_ml[2] ;
"true beta1:" beta1 ;


"calculate Hessian";
h=hessp(&logl,b_ml); /*this is minus the hessian because we minimize -L*/
"variance matrix";
vmat=inv(h);
vmat;
"" ;
"mean (value, t):" b_ml[1] b_ml[1]/sqrt(vmat[1,1]);
"MA-term (value, t):" b_ml[2] b_ml[2]/sqrt(vmat[2,2]);


/* averages and std dev of estimates */
"mean estimates" meanc(results);
"std of estimates" stdc(results) ;

" " ;

proc logl(b);
  local r,L,muvec,sigma2,beta0,beta1,i;
beta0=b[1] ;
beta1=b[2] ;


muvec=beta0+beta1*x ;


L=0 ;
i=0 ;
do while i<N ;
  i=i+1 ;
  if z1[i]==1 ; 
    L = L + ln(cdfn(muvec[i]));
  else ;
     L = L + ln(1-cdfn(muvec[i]));
  endif ;
endo ;   
  retp(-L); /*optmum minimizes*/
endp;
output off ;