Source

mulligan rubinstein bounds / plotBands.m

Full commit
clear; 
close all;
%% settings
outname='june2012';
medname = 'jun2012medR';
savefigs = true;
font = 'cmr17';
fontsize = 24;
ylo = -2.25;
yhi = 1.5;
groups = {'All','Single','Low education','High education', ...
          'Single and high education'};
groupAbbrev = {'all','single','loed','hied','singleHied'};
colnames = {'QR','Low 1970s','High 1970s','Low 1990s','High 1990s', ...
            'Low 1970s with stochastic dominance', ... 
            'High 1970s with stochastic dominance', ...
            'Low 1990s with stochastic dominance', ...
            'High 1990s with stochastic dominance', ...
            'Low 1970s with median restriction', ... 
            'High 1970s with median restriction', ...
            'Low 1990s with median restriction', ...
            'High 1990s with median restriction'};
yearAbbrev = {'7','9','7sd','9sd','7mr','9mr'};

%% load results
res =csvread(sprintf('%s.csv',outname))
medr= csvread(sprintf('%s.csv',medname));
tau = res(2:end,1)
assert(isequal(tau,medr(:,1)));
res = [res [zeros(1,size(medr,2)-1)*NaN; medr(:,2:end)]];
%outname = 'group';
tmp=csvread(sprintf('bs/%s%d.csv',outname,1));
nboot = size(ls(sprintf('bs/%s*.csv',outname)),1);
%assert(nboot==size(ls(sprintf('bs/%s*.csv',medname)),1), ...
%       ['different number of bootstrap rep for main results and median ' ...
%        'restriction.']);
mtmp = csvread(sprintf('bs/%s%d.csv',medname,1));
assert(size(mtmp,1)==size(tmp,1)-1);
tmp = [tmp [zeros(1,size(mtmp,2)-1)*NaN; mtmp(:,2:end)]];
bres = zeros([size(tmp) nboot]);
skip = zeros(nboot,1);
S = zeros([size(bres,1)-1 size(bres,2) nboot]);
oh = res(1,:);
for b=1:nboot
  clear tmp mtmp;
  try 
    tmp=csvread(sprintf('bs/%s%d.csv',outname,b));
    mtmp = csvread(sprintf('bs/%s%d.csv',medname,b));
    assert(size(mtmp,1)==size(tmp,1)-1);
    tmp = [tmp [zeros(1,size(mtmp,2)-1)*NaN; mtmp(:,2:end)]];
    assert(size(tmp,1)==size(bres,1) && size(tmp,2)==size(bres,2));
  catch
    fprintf(['Error while reading bs/%s%d.csv:\n' ...
             '  %s\n  Skipping this file.\n'],outname,b,lasterr());
    skip(b) = 1;
  end  
  if (skip(b)==0)
    bres(:,:,b) = tmp;
    S(:,:,b) = bres(2:end,:,b)-res(2:end,:);
    Soh(:,b) = bres(1,:,b) - oh;
  end
end
bres = bres(:,skip==0);
S = S(:,:,skip==0);
Soh = Soh(:,skip==0);
nboot = size(bres,2);
fprintf('\nSuccessfully read %d bootstrap results, encountered %d errors\n',nboot,sum(skip));
sub=2:10;
res = res(sub,:);
tau = res(:,1);
ng = numel(groups);
ny = numel(yearAbbrev);

% create directory for figures if needed
if (savefigs) 
  if (~(exist(sprintf('fig%s',outname),'dir')))
    mkdir(sprintf('fig%s',outname));
  end
end
set (0, "defaultaxesfontname", font);
set (0, "defaultaxesfontsize", fontsize) ;
set (0, "defaulttextfontname", font);
set (0, "defaulttextfontsize", fontsize) ;
if (savefigs)
  set (0, "defaultlinelinewidth", 8);
  set (0, "defaultlinemarkersize",12);
end
table = fopen(sprintf('table_%s.tex',outname),'w');
%% make plots and tables
for g=1:ng
  for c=1:ny    
    % indices of results in res and bres
    lols = 1 + (g-1)*2 + c;
    lh = 1 + 2*ng + (g-1)*2 + c;
    %if (c<=2) 
    lqr(c) = 1+(g-1)*2+2-mod(c,2);    
    lh = lh - 2*(floor((c-1)/2));
    lols = lols - 2*(floor((c-1)/2));
    %else 
    %lqr(c) = 1+(g-1)*2+c-2;
    %lh = lh - 2;
    %lols = lols - 2;
    %end    
    %if (c<=4)
    %  llo(c) = 1+ng*2+(c-1)*2*ng+g;
    %  lhi(c) = 1+ng*2+(c-1)*2*ng+ng+g;
    %else
      llo(c) = 1+ng*2*5+(c-5)*2*ng+g;
      lhi(c) = 1+ng*2*5+(c-5)*2*ng+g+ng;
      %end
    
    %display([g c lqr(c) llo(c) lhi(c)]);
    fprintf(table,'%10s & %6.3f & %6.3f & %6.3f & %6.3f & %6.3f\n', ...
            yearAbbrev{c}, oh([lols lh]), res(5,[lqr(c) llo(c) lhi(c)]));
    se = sqrt([mean((Soh(lols,:)).^2) ...
               mean((Soh(lh,:)).^2) ...
               mean((squeeze(S(5,lqr(c),:))).^2) ...
               mean((squeeze(S(5,llo(c),:))).^2) ...
               mean((squeeze(S(5,lhi(c),:))).^2) ]);
    fprintf(table,'%10s & (%0.3f) & (%0.3f) & (%0.3f) & (%0.3f) & (%0.3f)\n', ...
            '', se);
    %% plot results
    figure;
    plot(tau,res(:,lqr(c)),'k');
    line(tau,res(:,llo(c)),'Color',[1 0 0],'marker','^');
    line(tau,res(:,lhi(c)),'Color',[0 0 1],'marker','v');
    % pointwise bands
    lo = quantile(S(:,llo(c),:),0.05,3);
    lo(lo>0) = -lo(lo>0);
    hi = quantile(S(:,lhi(c),:),0.95,3);
    hi(hi<0) = -hi(hi<0);
    % adjust to be uniform
    p = @(f) mean(all(squeeze(S(:,llo(c),:))>=f*lo*ones(1,nboot) ... %,1)); 
                      & ...
                      squeeze(S(:,lhi(c),:))<=f*hi*ones(1,nboot),1));
    flo = 1;
    while(p(flo)<0.9) 
      flo = flo+0.05;
    end
    %flo = fminbnd(@(f) abs(p(f)-0.9),1,4)
    fhi = flo;
    %p = @(f) mean(all(squeeze(S(:,lhi(c),:))<=f*hi*ones(1,nboot),1));
    %fhi = 1;
    %while(p(fhi)<0.95)
    %  fhi = fhi+0.1;
    %end
    %fhi = fminbnd(@(f) abs(p(f)-0.95),1,4)
    line(tau,flo*lo+res(:,llo(c)),'Color',[1 0.5 0.5],'LineStyle','--'); %,'marker','^')
    line(tau,fhi*hi+res(:,lhi(c)),'Color',[0.5 0.5 1],'LineStyle','--'); %,'marker','v')
    %line(tau,zeros(size(tau)),'Color','k','LineStyle',':')
    grid('on');
    axis([0.1 0.9 ylo yhi]);
    xlabel('quantile');
    ylabel('log wage gap');
    legend('QR','Lower','Upper','location','northwest');
    set(gca,'FontName',font);
    if (savefigs)
      print('-depsc2',sprintf('fig%s/%s%s.eps',outname,...
                              groupAbbrev{g},yearAbbrev{c}));    
      close;
    else
      title([groups{g} yearAbbrev{c}]);
    end
  end % for c=1:ny
  %% change in gap
  for sd=0:2:4
    dqr = res(:,lqr(2+sd))-res(:,lqr(1+sd));
    dhi = res(:,lhi(2+sd))-res(:,llo(1+sd));
    dlo = res(:,llo(2+sd))-res(:,lhi(1+sd));
    dS(:,1,:) = S(:,llo(2+sd),:) - S(:,lhi(1+sd),:);
    dS(:,2,:) = S(:,lhi(2+sd),:) - S(:,llo(1+sd),:);
    % pointwise bands
    lo = quantile(dS(:,1,:),0.05,3);
    lo(lo>0) = -lo(lo>0);
    hi = quantile(dS(:,2,:),0.95,3);
    hi(hi<0) = -hi(hi<0);
    % adjust to be uniform
    p = @(f) mean(all(squeeze(dS(:,1,:))>=f*lo*ones(1,nboot) ...
                      & ...
                      squeeze(dS(:,2,:))<=f*hi*ones(1,nboot),1));
    flo = 1;
    while(p(flo)<0.9) 
      flo = flo+0.05;
    end
    %flo = fminbnd(@(f) abs(p(f)-0.9),1,4)
    fhi = flo;
    figure;
    plot(tau,dqr,'k');
    line(tau,dlo,'Color',[1 0 0],'marker','^');
    line(tau,dhi,'Color',[0 0 1],'marker','v');
    line(tau,flo*lo+dlo,'Color',[1 0.5 0.5],'LineStyle','--'); %,'marker','^')
    line(tau,fhi*hi+dhi,'Color',[0.5 0.5 1],'LineStyle','--'); %,'marker','v')
    grid('on');
    axis([0.1 0.9 -2.25 2.25]);
    xlabel('quantile');
    ylabel('log wage gap');
    legend('QR','Lower','Upper','location','northwest');
    set(gca,'FontName',font);
    if (savefigs)
      if (sd==2) 
        nameend='ChangeSD';
      elseif (sd==0)
        nameend='Change';
      elseif (sd==4)
        nameend='ChangeMR';
      end
      print('-depsc2',sprintf('fig%s/%s%s.eps',outname,...
                              groupAbbrev{g},nameend));    
    else
      title([groups{g} ' change']);
    end
  end % for sd=1:2
end
fclose(table);
if (savefigs)
  close all;
end