diary c:\teaching\04_spring\econ_241a\prob2\prob2_04jan27b.dia
%path(path,'c:\teaching\fall97\ec231a\prob1')

% 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
dur=dur.*(dur>0)+(dur==0);   % adding one for zero durations
sdur=[min(dur),max(dur),mean(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)];


% 2a.  log likelihood function
beta=[1;0;0;0;0];
Y=dur;
N=length(Y);
X=[ones(N,1),age,educ,white,locrat];
[N,K]=size(X);
LL0=loglik(beta,Y,X,zeros(N,1))

% 2b. log likelihood function with censoring
LL1=loglik(beta,Y,X,censor)

% 2c. analytic first derivatives
fd=fder(beta,Y,X,censor)

% 2d. numerical first derivatives
fd_num=zeros(K,1);
cc=0.001;
IK=eye(K);
for k=1:K,
    LLc=loglik(beta+cc*IK(:,k),Y,X,censor);
    fd_num(k,1)=(LLc-LL1)/cc;
    end,
[fd,fd_num]    
 
% 2d. numerical first derivatives
fd_num=zeros(K,1);
cc=cc^2;
IK=eye(K);
for k=1:K,
    LLc=loglik(beta+cc*IK(:,k),Y,X,censor);
    fd_num(k,1)=(LLc-LL1)/cc;
    end,
[fd,fd_num]    
 
% 2d. numerical first derivatives
fd_num=zeros(K,1);
cc=cc^2;
IK=eye(K);
for k=1:K,
    LLc=loglik(beta+cc*IK(:,k),Y,X,censor);
    fd_num(k,1)=(LLc-LL1)/cc;
    end,
[fd,fd_num]    

cc=0.001;
sd=sder(beta,Y,X,censor);
sd_num=zeros(k,k);
fd=fder(beta,Y,X,censor);
for k1=1:K,
    fd0=fd(k1,1)
        for k2=1:K,
        fd1=fder(beta+cc*IK(:,k2),Y,X,censor);
        sd_num(k1,k2)=(fd1(k1,1)-fd0)/cc;
        end,
    end,
sd
sd_num
sd-sd_num


% 2h. newton-raphson [1;0;0;0;0]
'NEWTON-RAHPSON'
beta_old=[1;0;0;0;0];
betai=[];
der=[];
stopp=0;
iter=0;
while stopp==0,
   iter=iter+1
   betai=[betai;[iter,beta_old']];
   sd=sder(beta_old,Y,X,censor);
   fd=fder(beta_old,Y,X,censor);
   beta=beta_old'
   first_derivative=fd'
   der=[der;fd'];
   dir=inv(sd)*fd;
   dir=dir/(1+sqrt(dir'*dir));
   beta_new=beta_old-dir;
   stopp=(sum(abs(fd))<0.001)+(iter>20);
   beta_old=beta_new;
end,
betai;
der;

% 2i. DFP
'DFP ALGORITHM'
[N,K]=size(X);
beta_k=[1;0;0;0;0];
A_k=eye(K);
fd_k=fder(beta_k,Y,X,censor);
stopp=0;
iter=0;
while stopp==0,
      iter=iter+1
      %obj_k=loglik(beta_k,Y,X,censor);
      d_k=-A_k*fd_k;
      d_k=d_k/(1+sqrt(d_k'*d_k))
      lambda_l=0;
      lambda_h=1;
      % find lambda_l, lambda_h so that optimum is in between
      q_l=loglik(beta_k+lambda_l*d_k,Y,X,censor);
      q_h=loglik(beta_k+lambda_h*d_k,Y,X,censor);
      [q_l,q_h];
      if q_h<q_l,
         lambda_hh=lambda_l+3*(lambda_h-lambda_l);
         q_hh=loglik(beta_k+lambda_hh*d_k,Y,X,censor);
         while q_hh<q_h,
               [lambda_l,lambda_h,lambda_hh];
               [q_l,q_h,q_hh];
               lambda_l=lambda_h;
               q_l=q_h;
               lambda_h=lambda_hh;    
               q_h=q_hh;
               lambda_hh=lambda_l+3*(lambda_h-lambda_l);
               q_hh=loglik(beta_k+lambda_hh*d_k,Y,X,censor);       
               [lambda_l,lambda_h,lambda_hh];
               [q_l,q_h,q_hh];
               end     
               [lambda_l,lambda_h,lambda_hh];
               [q_l,q_h,q_hh];
          q_h=q_hh;
          lambda_h=lambda_hh;
          end
   
      [lambda_l,lambda_h];
      [q_l,q_h];
      q_h-q_l;
      % find optimal lambda in range [lambda_l,lambda_h]
      cc=1-2/(1+sqrt(5));
      delta=lambda_h-lambda_l;
      lambda_m1=lambda_l+delta*cc;
      lambda_m2=lambda_h-delta*cc;
      q_m1=loglik(beta_k+lambda_m1*d_k,Y,X,censor);
      q_m2=loglik(beta_k+lambda_m2*d_k,Y,X,censor);
      i=0;
      out_gs=[i,lambda_l,lambda_h,lambda_h-lambda_l];
      [q_l,q_m1,q_m2,q_h];
      [lambda_l,lambda_m1,lambda_m2,lambda_h];
      while (lambda_h-lambda_l)>0.0000001,
            i=i+1; 
            if q_m2>q_m1,
               lambda_h=lambda_m2;
               q_h=q_m2;
               lambda_m2=lambda_m1;
               q_m2=q_m1;
               delta=lambda_h-lambda_l;
               lambda_m1=lambda_l+delta*cc;
               q_m1=loglik(beta_k+lambda_m1*d_k,Y,X,censor);
               else        
               lambda_l=lambda_m1;
               q_l=q_m1;
               lambda_m1=lambda_m2;
               q_m1=q_m2;
               delta=lambda_h-lambda_l;
               lambda_m2=lambda_h-delta*cc;
               q_m2=loglik(beta_k+lambda_m2*d_k,Y,X,censor);
               end,
            out=[i,lambda_l,lambda_h,lambda_h-lambda_l];
            out_gs=[out_gs;out];
            end
      lambda_k=(lambda_l+lambda_h)/2
      % update other parameters DFP
      p_k=lambda_k*d_k
      beta_k1=beta_k+p_k;
      fd_k1=fder(beta_k1,Y,X,censor);
      q_k=fd_k1-fd_k
      A_k1=A_k+(p_k*(p_k'))/(p_k'*q_k)-A_k*q_k*(q_k')*A_k/(q_k'*A_k*q_k)
      stopp=(sum(abs(fd_k1))<0.001)+(iter>25);
      beta_k=beta_k1;
      fd_k=fd_k1;
      A_k=A_k1;
      [iter,lambda_k,sum(abs(fd_k))];
      beta=beta_k'
      end,