Source

mulligan rubinstein bounds / plotMedRBands.m

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

%% load results
res =csvread(sprintf('%s.csv',outname))
tau = res(2:end,1)
tmp=csvread(sprintf('bs/%s%d.csv',outname,1));
nboot = size(ls(sprintf('bs/%s*.csv',outname)),1);
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
  try 
    bres(:,:,b) = csvread(sprintf('bs/%s%d.csv',outname,b));
  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)
    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 director 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", 4);
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+c;    
    else 
      lqr(c) = 1+(g-1)*2+c-2;
      lh = lh - 2;
      lols = lols - 2;
    end    
    llo(c) = 1+ng*2+(c-1)*2*ng+g;
    lhi(c) = 1+ng*2+(c-1)*2*ng+ng+g;
    %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('\tau');
    ylabel('log wage gap');
    legend('QR','Low','High','Low CB','High CB');
    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:2
    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('\tau');
    ylabel('log wage gap');
    legend('QR','Low','High','Low CB','High CB');
    set(gca,'FontName',font);
    if (savefigs)
      if (sd>0) 
        nameend='ChangeSD';
      else
        nameend='Change';
      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);