%diary c:\teaching\04_spring\econ_241a\prob1\prob1_04jan19.dia

% 1. read in the data
load duration.dat
dur=duration(:,1);
censor=duration(:,2);
age=duration(:,3);
educ=duration(:,4);
white=duration(:,5);
locrat=duration(:,6);
clear duration
sdur=[min(dur),max(dur),mean(dur),median(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)]
on=ones(length(dur),1);

% 2. log likelihood function at 5
beta=5;
ll0=loglik(dur,1-on,beta)

% 3. plot log likelihood function
beta=log(mean(dur))
nb=100;
bmax=10;
bmin=0;
bf=zeros(nb,2);
N=length(dur);
for i=1:nb,
    b=bmin+(i-1)*(bmax-bmin)/(nb-1);
    bf(i,1)=b;
    bf(i,2)=loglik(dur,1-on,b);
    end,
%plot(bf(:,1),bf(:,2),'-')
%title('Figure 1: Log Likelihood Function')
%xlabel('beta')
%ylabel('log likelihood function')
%print c:\teaching\04_spring\econ_241a\prob1\fig1

% 4. checking analytic and numerical first derivatives
beta=5;
c1=0.00001;
c2=0.000001;
ll0=loglik(dur,1-on,beta);
ll1=loglik(dur,1-on,beta+c1);
ll2=loglik(dur,1-on,beta+c2);
fder_num1=(ll1-ll0)/c1;
fder_num2=(ll2-ll0)/c2;
fder_ana=loglik_fder(dur,1-on,beta)
[fder_ana,fder_num1,fder_num2]

% 5. checking analytic and numerical second derivatives
beta=5;
c2=0.000001;
ll_fder_0=loglik_fder(dur,1-on,beta);
ll_fder_1=loglik_fder(dur,1-on,beta+c2);
sder_num=(ll_fder_1-ll_fder_0)/c2;
sder_ana=loglik_sder(dur,1-on,beta)
[sder_ana,sder_num]

% 6. Newton-Raphson iterations
% starting value: zero
% 20 iterations
out_nr=[];
beta=0;
crit=0;
max_crit=0.00001;
i=0;
while crit==0,
    i=i+1;
    fd=-loglik_fder(dur,1-on,beta);
    sd=-loglik_sder(dur,1-on,beta);
    dir=-inv(sd)*fd;
    dir=dir/(1+sqrt(dir'*dir));
    beta_new=beta-inv(sd)*fd;
    criter=sum(abs(fd))+sum(abs(beta-beta_new));
    crit=criter<max_crit;
    out=[i,beta,fd,beta_new-beta,dir,criter];
    beta=beta_new;
    out_nr=[out_nr;out];
    end,
out_nr
[out_nr(:,1:2),out_nr(:,4:5)]
[out_nr(4:5,3),out_nr(4:5,6)]
[out_nr(6:11,3),out_nr(6:11,6)]
pause

[nn,mm]=size(out_nr);
out_nr(nn-3:nn,:)

% 7. golden section
beta_l=0;
beta_h=1;
q_l=-loglik(dur,1-on,beta_l);
q_h=-loglik(dur,1-on,beta_h);
[q_l,q_h]
if q_h<q_l,
   beta_hh=beta_l+3*(beta_h-beta_l);
   q_hh=-loglik(dur,1-on,beta_hh);
   while q_hh<q_h,
         [beta_l,beta_h,beta_hh]
         [q_l,q_h,q_hh]
         %pause
         beta_l=beta_h;
         q_l=q_h;
         beta_h=beta_hh;    
         q_h=q_hh;
         beta_hh=beta_l+3*(beta_h-beta_l);
         q_hh=-loglik(dur,1-on,beta_hh);       
         end     
   q_h=q_hh;
   beta_h=beta_hh;
   end
   
[beta_l,beta_h]
[q_l,q_h]

cc=1-2/(1+sqrt(5));
delta=beta_h-beta_l;
beta_m1=beta_l+delta*cc;
beta_m2=beta_h-delta*cc;
q_m1=-loglik(dur,1-on,beta_m1);
q_m2=-loglik(dur,1-on,beta_m2);
i=0;
out_gs=[i,beta_l,beta_h,beta_h-beta_l]
[q_l,q_m1,q_m2,q_h]
[beta_l,beta_m1,beta_m2,beta_h]
while (beta_h-beta_l)>0.0000001,
      i=i+1; 
      if q_m2>q_m1,
         beta_h=beta_m2;
         q_h=q_m2;
         beta_m2=beta_m1;
         q_m2=q_m1;
         delta=beta_h-beta_l;
         beta_m1=beta_l+delta*cc;
         q_m1=-loglik(dur,1-on,beta_m1);
         else        
         beta_l=beta_m1;
         q_l=q_m1;
         beta_m1=beta_m2;
         q_m1=q_m2;
         delta=beta_h-beta_l;
         beta_m2=beta_h-delta*cc;
         q_m2=-loglik(dur,1-on,beta_m2);
         end,
      out=[i,beta_l,beta_h,beta_h-beta_l];
      out_gs=[out_gs;out];
      end

out_gs

%diary off