diary 
% 1. read in the data
load dur_2006.dat
dur=dur_2006;
format short g;
% 1. number of observations with zero duration
sum(dur(:,1)==0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dur(:,1)=dur(:,1).*(dur(:,1)>0)+(dur(:,1)==0);   % adding one for zero durations
censor=dur(:,2);
age=dur(:,3);
educ=dur(:,4);
white=dur(:,5);
locrat=dur(:,6);
dur=dur(:,1);

'summary statistics'
sdur=[min(dur),max(dur),mean(dur),std(dur)]
scensor=[min(censor),max(censor),mean(censor)]
sage=[min(age),max(age),mean(age)]
seduc=[min(educ),max(educ),mean(educ)]
slocrat=[min(locrat),max(locrat),mean(locrat)]
swhite=[min(white),max(white),mean(white)]
'number of observations'
n=length(dur)

'log likelihood function at minus four'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ll = (-4)*sum(1-censor)-exp(-4)*sum(dur)
ll0=[-4,loglik(-4,dur,censor)]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'plot log likelihood function'
lambda=(-10:0.1:-4)';
lllambda=0*lambda;
for i=1:length(lambda);
     lllambda(i,1)=loglik(lambda(i,1),dur,censor);
end,
plot(lambda,lllambda,'-')   % . , - , *, 
%axis([-10 -4 -40000 0])
xlabel('lambda')
ylabel('log likelihood function')
title('Figure 1: Log Likelihood Function') 

'maximum likelihood estimate'
lambdaml=log(sum(1-censor)/sum(dur))
'first derivative at -4:'
ader=fder(-4,dur,censor) %analytical der
c1=0.001;
numder1=(loglik(-4+c1,dur,censor)-loglik(-4,dur,censor))/c1;
c2=0.000001;
numder2=(loglik(-4+c2,dur,censor)-loglik(-4,dur,censor))/c2;
[ader,numder1,numder2]
'second derivative at -4'
ader=sder(-4,dur,censor);
c1=0.001;
numder1=(fder(-4+c1,dur,censor)-fder(-4,dur,censor))/c1;
c2=0.000001;
numder2=(fder(-4+c2,dur,censor)-fder(-4,dur,censor))/c2;
[ader,numder1,numder2]

'newton-raphson with -4 as starting value'
lambda0=-4;
ni=15;
betai=zeros(ni,4); % do 15 iterations
for i=1:15,
   der=fder(lambda0,dur,censor);
   secder=sder(lambda0,dur,censor);
   lambda1=lambda0-inv(secder)*der;
   betai(i,:)=[lambda0,der,secder,(lambda1-lambdaml)/(lambda0-lambdaml)];
   lambda0=lambda1;
end,
betai

'quadratic fit'
lambdal=-8;
lambdah=-4;
lambdam=-6;
betai=zeros(14,5);
lambdanold=lambdam;
for i=1:14;
     ql=loglik(lambdal,dur,censor);
     qm=loglik(lambdam,dur,censor);
     qh=loglik(lambdah,dur,censor);
     dlh=lambdal-lambdah;
     dlm=lambdal-lambdam;
     dmh=lambdam-lambdah;
    dllhh=lambdal*lambdal-lambdah*lambdah;
    dllmm=lambdal*lambdal-lambdam*lambdam;
     dmmhh=lambdam*lambdam-lambdah*lambdah;
     lambdan=0.5*(ql*dmmhh-dllhh*qm+dllmm*qh)/(dmh*ql-dlh*qm+dlm*qh);
     dif=(lambdan-lambdaml)/(lambdanold-lambdaml);
     qn=loglik(lambdan,dur,censor);
     betai(i,:)=[lambdal,lambdam,lambdah,lambdan,dif];
     [lambdal,lambdam,lambdah,lambdan,dif]
    [ql,qm,qh,qn]
     rr=[ql,qm,qh,qn;lambdal,lambdam,lambdah,lambdan]';
     [s,l]=sort(rr(:,1));
     rr=rr(l,:);
     rr=rr(2:4,:);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %sort
     [s,l]=sort(rr(:,2));
     rr=rr(l,:);
     lambdal=rr(1,2);
     lambdam=rr(2,2);
     lambdah=rr(3,2);
end,   
betai


lambdal=-8;

lambdah=-4;

betai=zeros(20,3);

c1=1-2/(1+sqrt(5))

c2=2/(1+sqrt(5))

for i=1:20,

    dif=lambdah-lambdal;

    lambdam1=lambdal+c1*dif;

    lambdam2=lambdal+c2*dif;

    qm1=loglik(lambdam1,dur,censor);

    qm2=loglik(lambdam2,dur,censor);

    [lambdam1,lambdam2];

    [qm1,qm2];

    if qm1<qm2,

       res=[lambdal,lambdah,(0.5*(lambdah+lambdam1)-lambdaml)/(0.5*(lambdah+lambdal)-lambdaml)];

       lambdal=lambdam1;

       else

       res=[lambdal,lambdah,(0.5*(lambdam2+lambdal)-lambdaml)/(0.5*(lambdah+lambdal)-lambdaml)];

       lambdah=lambdam2;

       end,

    betai(i,:)=res;

    end,

    

betai    

    

    

    



       

diary off