%% load in the data from workplacesmoking.xlsx
%% column 1 is the dependent variable while
%% columns 2-(k+1) are independent variables
%% column 2 contains the constant
[w,varlist]=xlsread('workplacesmoking.xlsx');

dv=varlist{1,1};
%% get dimension of w
nk1=size(w);

% number of observations
n=nk1(1);

% number of independent variables
% k+1 is the no of columns of w.  there are k
% covariates (including the constant)
k=nk1(2)-1;

%% extract y which is in the 1st column of w
y=w(:,1);


%% extract x which is in columns 2 through k
x=w(:,2:(k+1));

%% generate starting values--which are from
%% an OLS model 

xpxi=inv(x'*x);
beta_start=xpxi*x'*y;


%% start quasi-newton maximum likelihood
%% set up initial values
cc=1e-8;          %%convergence criteria
maxit=200;        %%maximum number of interations

beta=beta_start;
logli=logitloglike(y,x,beta);
grad=logitgrad(y,x,beta);
hess=logithess(y,x,beta);
cov=inv(-hess);

%% open file to print out results
file1=fopen('logit_mle_output2.txt','w');
fprintf(file1,'----------------------------------------------------------------\n');


%% now iterate to obtain solution
for it=1:maxit;
betalast=beta;
loglast=logli;
fprintf(file1,'iteration number %4.f, loglikelihood = %9.4f \n' ,it,logli);

step=1;
beta=betalast+step*cov*grad;
disp(it)
logli=logitloglike(y,x,beta);

j=0;
display(logli)
display(loglast)

%% if the log likelihood has not improved
%% keep cutting step size in half

if logli<loglast;
    for subit=1:10
    j=j+1;
    display(j)
    step=0.5*step;
    beta=betalast+step*cov*grad;
    logli=logitloglike(y,x,beta);
    if logli>=loglast, break, end;
    end;
end;
    
display(j)

if j>=10, break, end;
grad=logitgrad(y,x,beta);
hess=logithess(y,x,beta);
cov=inv(-hess);
criti=grad'*cov*grad
if criti<cc, break,end;
end;

display(j)


if j>=10;
    fprintf(file1,'----------------------------------------------------------------\n');
    fprintf(file1,'warning, after 10 halvings of correction vector, loglikelihood  \n');
    fprintf(file1,'has not improved.  the results are not reliable                 \n');
    fprintf(file1,'----------------------------------------------------------------\n');
end;

if it>=maxit && criti>cc;
    fprintf(file1,'----------------------------------------------------------------\n');
    fprintf(file1,'warning, maximum iterations reached and program has not         \n');
    fprintf(file1,'converged.  the results are not reliable                        \n');
    fprintf(file1,'----------------------------------------------------------------\n');
end;


%% print out final results

ym=sum(y)/n;
dof=n-k;
stderr=sqrt(diag(cov));
zscore=beta./stderr;

pvalue=2*(1-normcdf(zscore,0,1));

fprintf(file1,'----------------------------------------------------------------\n');
fprintf(file1,'MLE estimates of Poisson model     \n');
fprintf(file1,'----------------------------------------------------------------\n');
fprintf(file1,'Dependent variable =  %s \n' ,dv);
fprintf(file1,'Mean of dep variable = %9.4f \n' ,ym);
fprintf(file1,'----------------------------------------------------------------\n');
fprintf(file1,'analysis of MLE procedure              \n');
fprintf(file1,'\n');
fprintf(file1,'observations           = %9.f \n' ,n);
fprintf(file1,'parameters             = %9.f \n' ,k);
fprintf(file1,'degrees of freedom     = %9.f \n' ,dof);
fprintf(file1,'number of iterations   = %9.f \n' ,it);
fprintf(file1,'final convergence crit = %e \n' ,criti);
fprintf(file1,'value of log likeliood = %9.4f \n' ,logli);

c1='Covariate'; c2='Beta'; c3='Std Error'; c4='zscore'; c5='P-value: B=0';
fprintf(file1,'----------------------------------------------------------------\n');
fprintf(file1,'%12s %12s %12s %12s %12s \n', c1,c2,c3,c4,c5);
fprintf(file1,'----------------------------------------------------------------\n');
for i=1:k;
    rowname=varlist{1,i+1};
    fprintf(file1,'%12s %12.6f %12.6f %12.6f %12.6f \n', rowname,beta(i,1),stderr(i,1),zscore(i,1),pvalue(i,1));
    end;
fprintf(file1,'----------------------------------------------------------------\n');
fclose(file1);