%% 14.382 Problem Set 2, Question 4
% Paul Schrimpf, March 2008 
%
% Recreates main results of Betrand, Duflo, and Mullainathan (2004).
% Also, applies bias and size corrected estimator of Hausman and
% Kuersteiner (2007) to this setup.
%
% See also SIMULATE, OLSDD, FGLSDD, AGGDATA

%% Load Data
% Data is from the CPS-MORG.  Some pre-processing (mainly merging of
% various years into single file) was done in stata.  
% Here, we read the data from a .csv file and save it as a structure.  
% The structure that we choose to save our data in is meant to be a
% generic diff in diff setup with 2 level clustering.  The fields of data
% are:
%  
% * y = outcome variable, size number of observations by 1
% * x = regressors, size number of observations by number of regressors
% * t = time index, size number of observations by 1
% * s = state index, size number of observations by 1
% * p = policy indicator, size number of observations by 1
%

clear;
profile on;
% Whenever importing data, an easy solution is to use "Import Data" under
% the file menu.  It gives you a GUI to import data and will
% automatically generate code. 
cps = importdata('cpsSample.csv',',',1);
% break data into nicely named variables
for c=1:size(cps.data,2)
  eval([cps.colheaders{c} '= cps.data(:,c);']);
end

% only use 25-55 yr old white females with positive earnings
use = earnwke>0 & age>=25 & age<=55 & white==1 & female==1;
y = log(earnwke(use));
t = year(use)-1978;
s = state(use);
x = [yrsEdu(use) age(use) age(use).^2 age(use).^3 age(use).^4];

% create original data structure
data.y = y;
data.x = x;
data.s = s;
data.t = t;
data.p = zeros(length(y),1);
% clear unneeded variables
clear cps y t s x use;

%% Run Simulations
% Just loops through simulations.  All combinations of N number of states
% with T time periods are run S times.  Each simulation is run with both
% a 0 and 2% policy effect.  Each time, the estimated policy effect, its
% variance, and whether or not the null hypothesis of 0 effect is
% rejected at a 5% level is saved.  For each simulation we estimate the
% effect using OLSDD with and without clustering and using FGLS with and
% without bias correction.  

%N = [6 10 20 50]; % number of states to use
%T = [3 5 10 20];  % number of time periods
%S = 100; % number of simulations
N = 10;
T = 5;
S = 20;
for n=1:length(N)
  for t=1:length(T)
    if (T(t)>N(n)-1) % too few T for N
      continue;
    end
    for s=1:S
      % simulate with no effect
      simd = simulate(data,N(n),T(t),0); 
      [bOLS(1,n,t,s), testOLS(1,n,t,s,:) var(1,n,t,s,:) ] = olsDD(simd);
      asimd = aggData(simd);
      [bGLS(1,n,t,s,:), testGLS(1,n,t,s,:) vGLS(1,n,t,s,:)] = fglsDD(asimd);
      % simulate with 2% effect
      simd = simulate(data,N(n),T(t),0.02); 
      [bOLS(2,n,t,s), testOLS(2,n,t,s,:) var(2,n,t,s,:)] = olsDD(simd);
      asimd = aggData(simd);
      [bGLS(2,n,t,s,:), testGLS(2,n,t,s,:) vGLS(2,n,t,s,:)]=fglsDD(asimd);
      fprintf('%4d of %4d\n',s,S);
    end
  end
end

profile viewer;
%% Print Results
% Prints the results in latex formatted tables.

% save results just in case
save(['sim_' datestr(now)],'bOLS','testOLS','var','bGLS','testGLS','vGLS');

% print results
of = fopen('sim.tex','w');
fprintf(of,['Estimator & N & T & Rejection Rate No Effect & Rejection Rate 2' ...
' % Effect \\\\ \\hline\n']);
estOLS = {'OLS' ...
          'OLS state-year clustered' ...
          'OLS state clustered'};
estGLS = {'FGLS' ...
          'FGLS bias-corrected' ...
          'FGLS bias and size corrected' ...
          'OLS aggregate'};
for n=length(N):(-1):1
  for t=length(T):(-1):1
    if (T(t)>N(n)-1) % too few T for N
      continue;
    end
    for i=1:3;
      fprintf(of,'%s & %d & %d & %7.3g & %7.3g \\\\ \n', ...
              estOLS{i},N(n),T(t),mean(testOLS(1,n,t,:,i)), ...
              mean(testOLS(2,n,t,:,i)));  
    end
    i = 4;
    fprintf(of,'%s & %d & %d & %7.3g & %7.3g \\\\ \n', ...
            estGLS{i},N(n),T(t),mean(testGLS(1,n,t,:,i)), ...
            mean(testGLS(2,n,t,:,i)));  
    for i=1:3
      fprintf(of,'%s & %d & %d & %7.3g & %7.3g \\\\ \n', ...
              estGLS{i},N(n),T(t),mean(testGLS(1,n,t,:,i)), ...
              mean(testGLS(2,n,t,:,i)));  
    end
    fprintf(of,'\\hline \n');
  end
end
