% HW4.M
% Solutions to HW4
% calls logit_mle.m and logit_weighted.m 
% sample results are given at the end
% 03/6/06

load hw4b.dat  %this is an edited version of the data, removing the header info

% Logit MLE
y=hw4b(:,1);
w=hw4b(:,2);
n = length(w);
x = [ones(n,1) w];
[beta_ML,se_ML] = logit_mle(y,x); 
beta_ML
se_ML

% weighted logit
n0 = sum(y==0);
n1 = sum(y==1);
Q0 = 1/4;
Q1 = 3/4;
weights = (Q0/(n0/n)).*(y==0) + (Q1/(n1/n)).*(y==1);
beta_weighted = logit_weighted(y,x,weights);
beta_weighted

% nonparametric bootstrap
B=1000;
nb_draws=NaN*ones(B,2);
constant_vec = ones(n,1);

for b=1:B,

  % sample with replacement from observed data points
  % first calculate a random draw for integers between 1 and n:
  b_index = ceil( n*rand(n,1) );  

  % then use these integers as index of the original data
  x_b = [constant_vec w(b_index)]; 
  y_b = y(b_index);

  % calculate MLE and record
  [beta_nb, se_nb] = logit_mle(y_b, x_b);
  nb_draws(b,:) = beta_nb';

end;


% parametric bootstrap
pb_draws=NaN*ones(B,2);

for b=1:B,

  % sample with replacement from observed W data points
  b_index = ceil( n*rand(n,1) );  
  x_b = [constant_vec w(b_index)]; 
  
  % generate y according to probit model based on MLE estimate
  xbb = x_b*beta_ML;
  Eyx = exp(xbb)./(1+exp(xbb));
  y_b =  rand(n,1)<Eyx;

  % calculate MLE and record
  [beta_pb, se_pb] = logit_mle(y_b, x_b);
  pb_draws(b,:) = beta_pb';

end;

% summarize standard deviations and compare to conventional SEs

disp('Standard Error Summary')
se_ML'
SE_NB = std( ( nb_draws - ones(B,1)*beta_ML' ))
SE_PB = std( ( pb_draws - ones(B,1)*beta_ML' ))

% note the se_ML is standard error for betahat-beta, not \sqrt(n)(betahat-beta)
% so i've calculated the bootstrap se's analogously
% alternatively, could rescale the se_ML to be for \sqrt(n)(betahat-beta)

% RESULTS: 
% beta_ML =
%
%  -3.0416
%   3.4614
%
% se_ML =
%
%  0.91702
%  0.90565
%
% beta_weighted =
%
%  -2.2028
%   3.5054
%
% Standard Error Summary
% se_ML =
%
%  0.91702  0.90565
%
% SE_NB =
%
%  1.0642  1.0629
%
% SE_PB =
%
%  2.7458  2.8938    
%
% bootstrap standard errors are somewhat different
% occasional outliers in bootstrap distribution
% very rarely, the ML fails to converge, but this did not happen here
%
% end hw4.m