%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 14.662: Problem Set 3 solution (s3)
% Fall 2015
% P. Hull
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all 
close all
cd 'C:\Users\Hull\Dropbox (MIT)\14.662-spring-2015\problem-sets\ps3';

% optimization options
options=struct('Display','iter','Diagnostics','on','MaxFunEvals',10000,'MaxIter',10000,'TolFun',1e-8);
 
% set seed
s = RandStream('mcg16807','Seed',0);
RandStream.setGlobalStream(s)

% import data 
cpsdata = csvread('morg11_lwage.csv');
  lwage = cpsdata(:,5);
  wght  = cpsdata(:,4);
  
% starting values
theta_0 = [1,2,log(0.5),log(0.5),0];

% mle
theta_hat = fminsearch(@(x)maxbivnorm(x,lwage,wght),theta_0,options);
theta_hat_rf = theta_hat;
theta_hat_rf(3) = exp(theta_hat(3));
theta_hat_rf(4) = exp(theta_hat(4));
theta_hat_rf(5) = (exp(theta_hat(5))-1)/(exp(theta_hat(5))+1);
theta_hat_rf %#ok

% plot
[fx,x] = ksdensity(lwage,'npoints',200);
wghtx = ones(1,200);
[temp,fx_hat] = maxbivnorm(theta_hat,x',wghtx');
figure(1)
hold on
scatter(x,fx,'.r')
plot(x,fx_hat,'LineWidth',3)
legend('empirical','fitted')
xlabel('log(wage)')
hold off

%theta_hat_rf = [2.81,1.74,0.647,0.637,-0.03]; %# authors' estimates

% draw sample based on estimates
mu = theta_hat_rf(1,1:2);
sigma = [theta_hat_rf(3)^2 theta_hat_rf(5)*theta_hat_rf(3)*theta_hat_rf(4); ...
         theta_hat_rf(5)*theta_hat_rf(3)*theta_hat_rf(4) theta_hat_rf(4)^2];
potential_ln_w = mvnrnd(mu,sigma,100000);
potential_w = exp(potential_ln_w);
max_w = max(potential_w(:,1),potential_w(:,2));

% estimate parameters
potential_e = (0.75*potential_w).^0.5;
max_e = max(potential_e(:,1),potential_e(:,2));
potential_y = potential_w.*potential_e;
max_y = max(potential_y(:,1),potential_y(:,2));
w_theta = sort(max_w(max_y==potential_y(:,1)==1,:));
w_phi   = sort(max_w(max_y~=potential_y(:,1)==1,:));
e_theta = sort(max_e(max_y==potential_y(:,1)==1,:));
e_phi   = sort(max_e(max_y~=potential_y(:,1)==1,:));
y_Theta = e_theta'*w_theta;
y_Phi = e_phi'*w_phi;
alpha = y_Theta/(y_Theta+y_Phi)             %#ok

% optimal tax data
optimal = load('MTR.mat');
figure(2)
plot(optimal.wage,optimal.mtr,'LineWidth',3)
xlabel('wage')
ylabel('optimal marginal tax rate')

% fit optimal tax schedule to simulated wages. 
% you guys probably have a more efficient way of doing this
 potential_t = nan(size(potential_w,1),2);
 for j = 1:2
   for i = 1:size(potential_w,1)
       for k = 1:size(optimal.wage,1)-1
            if optimal.wage(k)<= potential_w(i,j) && ...
                optimal.wage(k+1)>potential_w(i,j)
            potential_t(i,j)=optimal.mtr(k);
            end
       end
       if optimal.wage(size(optimal.wage,1)) <= potential_w(i,j)
           potential_t(i,j)=optimal.mtr(size(optimal.wage,1));
       end
   end
 end
                     
potential_e = ((1-potential_t).*potential_w).^0.5;
max_e = max(potential_e(:,1),potential_e(:,2));
potential_y = potential_w.*potential_e;
max_y = max(potential_y(:,1),potential_y(:,2));
w_theta = sort(max_w(max_y==potential_y(:,1),:));
w_phi   = sort(max_w(max_y~=potential_y(:,1),:));
e_theta = sort(max_e(max_y==potential_y(:,1),:));
e_phi   = sort(max_e(max_y~=potential_y(:,1),:));

figure(3)
plot(w_theta,e_theta,'LineWidth',3)
grid on
axis([0, 250, 0, 15])
xlabel('wage')
ylabel('optimal effort')

istheta = max_w==potential_w(:,1);
bins = logspace(log10(min(potential_w(:,1))),log10(max(potential_w(:,1))),100);
[~,bin_id] = histc(potential_w(:,1), [bins(1:end-1) Inf]);
bin_id = bin_id+1;
theta_share = accumarray(bin_id, istheta, [length(bins) 1], @mean);
w_theta_iftheta = potential_w(:,1).*istheta;
avg_theta_w = accumarray(bin_id, w_theta_iftheta, [length(bins) 1], @mean);

figure(4)
plot(bins,theta_share,'LineWidth',3)
grid on
xlabel('wage')
ylabel('share')
