/* dur_ex_2007_2.txt*/

/* maximum likelihood estimation of the duration distribution of the 
   right-censored, length-biased search data from the CPS */

/* assume that the population distribution of unemployment spells is Weibull
  with pdf
   pdf(t) = gamma/t  * (t/alpha)^gamma  * exp(-(t/alpha)^gamma )
   S(t) = exp(-(t/alpha)^gamma )
   mean(t) = alpha * gamma_function( 1 + 1/gamma )
*/

/* read in the data */

load x[2907,4] = cps2002.txt;

t = selif(x[.,2],x[.,2] .> 0);

g_gam = 1;
g_alpha = 1;

fn w_survivor(t,gam,alpha) = exp( - (t/alpha)*gam );
fn w_mean(gam,alpha) = alpha * gamma(1+1/gam);

fn w_lbrc_den(t,gam,alpha) = w_survivor(t,gam,alpha) / w_mean(gam,alpha) ;

library maxlik;
maxset;

proc llf(b,x);
local gam,alpha;
gam = exp(b[1]);
alpha = exp(b[2]);
retp( ln(w_lbrc_den(t,gam,alpha) ) );
endp;
 
b0 = {0,0};

{bmle,fval,grad,hessian,retcode}=maxprt(maxlik(t,1,&llf,b0));

/* convert back to parameters of interest - invariance property of maximum likelihood estimators */

gamma_hat = exp(bmle[1]);
alpha_hat = exp(bmle[2]);

print "mle of gamma" gamma_hat;
print "mle of alpha" alpha_hat;


/* use eqsolve to find the estimates */

/* express FOC using numerical derivatives */

llf1 = 0; llf2 = 0; llf3 = 0; gam_eps = 0; alpha_eps = 0;

proc foc(b);
local i,epsilon,gam,alpha,eq1,eq2;
epsilon = .00001;
gam = exp(b[1]);
gam_eps = exp(b[1] + epsilon);
alpha = exp(b[2]);
alpha_eps = exp(b[2]+epsilon);
llf1 = ln(w_lbrc_den(t,gam,alpha));
llf2 = ln(w_lbrc_den(t,gam_eps,alpha));
llf3 = ln(w_lbrc_den(t,gam,alpha_eps));
eq1 = sumc(llf2-llf1) / epsilon ;
eq2 = sumc(llf3-llf1) / epsilon ;
retp( eq1 | eq2 );
endp;

b0 = {0,0};

{b_eqsolve,retcode}=eqsolve(&foc,b0);

/* print solution from eqSolve, and translate back into "natural parameters" */
print "point estimates from eqsolve";
print b_eqsolve;
print "translated back into parameters of interest";
print exp(b_eqsolve);


end;