%diary c:\teaching\04_spring\econ_241a\prob3\prob3_04feb2.dia


% 1. read in the data
load duration.dat
dur=duration(:,1);
duration=duration(dur>=80,:);
[length(duration),length(dur)]

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)]
size(dur)


% 2a.  log likelihood function with censoring and truncation
beta=[1;0;0;0;0];
Y=dur;
N=length(Y);
X=[ones(N,1),age,educ,white,locrat];
[N,K]=size(X);
LL1=loglik(beta,Y,X,censor)

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

% 2d. numerical first derivatives
fd_num=zeros(K,1);
cc=0.001^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]    
 


% 2i. DFP
'DFP ALGORITHM'
[N,K]=size(X);
beta_k=[1;zeros(K-1,1)];
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,
      
beta=beta';
sd=sder(beta,Y,X,censor);
var=inv(sd);
se=sqrt(diag(var));
beta_se=[beta,se]
pause
LL_unr=loglik(beta,Y,X,censor)

var00=var(4:5,4:5)
beta0=beta(4:5)
WALD=beta0'*inv(var00)*beta0


X=[ones(N,1),age,educ];
% 2i. DFP
'DFP ALGORITHM'
[N,K]=size(X);
beta_k=[1;zeros(K-1,1)];
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,
      
beta=beta';
sd=sder(beta,Y,X,censor);
var=inv(sd);
se=sqrt(diag(var));
beta_se=[beta,se]
LL_r=loglik(beta,Y,X,censor)

LR=-2*(LL_unr-LL_r)

beta=[beta;0;0];
X=[ones(N,1),age,educ,white,locrat];
[N,K]=size(X);
fd=fder_full(beta,Y,X,censor);
on=ones(N,1);
LM=on'*fd*inv(fd'*fd)*fd'*on
on'*fd
inv(fd'*fd)

diary off