global Y X

alldata=csvread('PS1discreteOutcomeTwovar.csv');
Y=alldata(:,1);
vars=alldata(:,2:end);
d=size(Y);
cons=ones(d(1,1),1);
X=[vars cons];

% Maximum Likelihood parameter starting values.  Let's use ones.
beta0=ones(size(X,2),1)

MAXITERS=[];
PARAMTOL=0.00001;
LLTOL=0.00001;
NAMES={'slope1' 'slope2' 'intercept'}; % Labels for independent variables


disp('Start estimation');
tic; % Start stopwatch timer
param=beta0;  % starting values: must be col vector
options=optimset('LargeScale','off','Display','iter','GradObj','off',...
    'MaxFunEvals',10000,'MaxIter',MAXITERS,'TolX',PARAMTOL,'TolFun',LLTOL,'DerivativeCheck','off');
[paramhat,fval,exitflag,output,grad,hessian]=fminunc(@nllLogit,param,options);

disp(' ');
if exitflag == 1
  disp('Convergence achieved.');
elseif exitflag == 2
  disp('Convergence achieved by criterion based on change in parameters.');
elseif exitflag == 3
  disp('Convergence achieved by criterion based on change in log-likelihood value.');
else
    disp('Convergence not achieved.');
    return
end

disp(['Value of the log-likelihood function at convergence: ' num2str(-fval)]);

%Calculate standard errors of parameters
disp('Taking inverse of hessian for standard errors.');
ihess=inv(hessian);
%ihessmean=ihess/NCS; % for long-form logit data
%stderr=sqrt(diag(ihessmean));
stderr=sqrt(diag(ihess));
disp(['The value of grad*inv(hessian)*grad is: ' num2str(grad'*ihess*grad)]);

disp(' ');
disp('ESTIMATION RESULTS');
disp(' ')
disp('              ---------------------------- ');
disp('                Est         SE      t-stat');
for r=1:size(NAMES,2);
    fprintf('%-10s %10.4f %10.4f %10.4f\n', NAMES{1,r}, [paramhat(r,1) stderr(r,1) paramhat(r,1)./stderr(r,1) ]);
end
disp(' ');

disp(' ');
disp('You can access the estimated parameters as variable paramhat,');
disp('the gradient of the negative of the log-likelihood function as variable grad,');
disp('the hessian of the negative of the log-likelihood function as variable hessian,');
disp('and the inverse of the hessian as variable ihess.');
disp('The hessian is calculated by the BFGS updating procedure.');

