new;
       cls;
    library pgraph;

       /*========DGP=========*/
       T = 200;
       Ymat = zeros(t,1);
       X1mat = rndn(t,1);
       X2mat = rndn(t,1);
       xmat = X1mat~X2mat;
       b1 = 2;
       b2 = -1;
       b = b1|b2;
       tru_para = b;
    
       itr = 1;
       do until itr > t;
          et = rndn(1,1);
          xb = xmat[itr,.]*b;
          if et < xb; ymat[itr] = 1; endif;
          itr = itr + 1;
       endo;

//       xy(T,ymat);
//stop;

       /*======== Estimation =========*/
       b1_in = b1 + 0.5;
       b2_in = b2 - 0.5;
       prmtr_in = b1_in|b2_in;
       {xout,fout,cout,gout} = qnewton(&lik_f,prmtr_in);
       cov = inv(hessp(&lik_f,xout));
       T_val = xout./sqrt(diag(cov));
       "===============================";
       "    True value      Estimates      t value" ;
       "------------------------------------------------";
       tru_para~xout~t_val;
       "================================";
    
       /*======== Likelihood Function =======*/
       proc lik_f(prmtr);
       local bb,lnL,cuml_den,lnf,itr;
       bb = prmtr[1:2];
       lnL = 0;
    
       itr = 1;
       do until itr > t;
          cuml_den = cdfn(xmat[itr,.]*bb);
          lnf = ymat[itr]*ln(cuml_den) + (1 - ymat[itr])*ln(1-cuml_den);
          lnL = lnL + lnf;
          itr = itr + 1;
       endo;
    
       retp(-lnL);
       endp;
       end;