*Figures for the paper "The Opportunity Costs of Mandatory Military Service: Evidence from a Draft Lottery"


#delimit;

prog drop _all;
set matsize 11000;
set seed 1234;

*figure 2 - Service probability by lottery draw;
use pnr lotterydraw drawlo shy served using "opp-cost-master", clear;
duplicates drop pnr, force;
egen threshold=min(lotterydraw) if drawlo==0, by(shy);
egen tall=max(threshold), by(shy);
gen lotnorm=lotterydraw-tall-1;
egen sprob=mean(served), by(lotnorm);
gen lot100=int((lotnorm-50)/100);
egen s100=mean(served), by(lot100);
egen c100=sum(served), by(lot100);
replace lot100=lot100*100;
tw (scatter s100 lot100 if lot100>-20000 & lot100<20000 & lot100!=0, graphregion(color(white)) legend(off) msymbol(X)) 
   (lpoly s100 lot100  if lot100>-20000 & lot100<0, deg(2) n(50) lcolor(red))  
   (lpoly s100 lot100  if lot100< 20000 & lot100>0, deg(2) n(50) lcolor(red)), 
   ytitle(proportion serving) xtitle(lottery draw minus threshold);
graph save fig2, replace;
graph export fig2.pdf, replace;



*figure 3 - service age by AFQT;
use spaer_year sessan_year pnr qafqt drawlo served serviceage if spaer_year>sessan_year using "opp-cost-master", clear;
duplicates drop pnr, force;
rename drawlo drafted;
su served if drafted==0;
scalar aa=r(mean);
keep if served;
collapse (count) n=served, by(serviceage qafqt);
reshape wide n, i(serviceage) j(qafqt);
graph bar n4 n3 n2 n1 if serviceage<25, over(serviceage) legend(lab(1 "leftmost - lowest") lab(2 "center left - second quartile AFQT") lab(3 "center right - third") lab(4 "rightmost - highest quartile AFQT") order(1 2 3 4)) ytitle("frequency") b1title(service age) graphregion(color(white));
graph save fig3, replace;
graph export fig3.pdf, replace;


set seed 1234;

use age learn spaer_year sessan_year afqt pnr drawlo pnr served feduc meduc etdisp height birthweight fromsingleparentfamily dane outofhomecare if age>=25 & age<=35 & learn!=. & spaer_year>sessan_year & afqt!=. & afqt>=28 using "opp-cost-master", clear;
duplicates drop pnr, force;
rename drawlo drafted;
keep pnr drafted served learn afqt feduc meduc etdisp height birthweight fromsingleparentfamily dane outofhomecare;
save compliers-count, replace;

for X in any afqt feduc meduc etdisp height birthweight: 
  xtile X4=X, n(4)
\ replace X4=(X4==4);
gen notdane=1-dane;
gen un=1;

prog define matt;
  matrix BN=BN\(y11n,y10n,y01n,y00n);
end;



matrix BN=(0,0,0,0);
for X in any un afqt4 feduc4 meduc4 etdisp4 height4 birthweight4 fromsingleparentfamily notdane outofhomecare: 
  su learn if drafted==1 & served==1 & X==1
\ scalar y11n=r(N)
\ su learn if drafted==1 & served==0 & X==1
\ scalar y10n=r(N)
\ su learn if drafted==0 & served==1 & X==1
\ scalar y01n=r(N)
\ su learn if drafted==0 & served==0 & X==1
\ scalar y00n=r(N)
\ matt;
clear;
svmat BN;
save ncomply101, replace;

rename BN1 n11;
rename BN2 n10;
rename BN3 n01;
rename BN4 n00;

gen d0=n01/(n00+n01);
gen d1=n11/(n11+n10);
gen cc=d1-d0;
gen aa=d0;
gen nn=1-d1;
gen cc1=n11/(n00+n10+n01+n11)*(1-aa);
gen cc0=1-aa-nn-cc1;

log using compliers-counts, text replace;
capture all afqt4 feduc4 meduc4 etdisp4 height4 birthweight4 fromsingleparentfamily notdane outofhomecare;
list aa cc1 cc0 nn;
log close;

prog define matt2;
  matrix BN=BN\(`1',y11n,y10n,y01n,y00n);
end;


use compliers-count, clear;
xtile pafqt=afqt, n(100);
matrix BN=(0,0,0,0,0);
for Y in num 1/80:
  su learn if drafted==1 & served==1 & pafqt>=Y & pafqt<=Y+20
\ scalar y11n=r(N)
\ su learn if drafted==1 & served==0 & pafqt>=Y & pafqt<=Y+20
\ scalar y10n=r(N)
\ su learn if drafted==0 & served==1 & pafqt>=Y & pafqt<=Y+20
\ scalar y01n=r(N)
\ su learn if drafted==0 & served==0 & pafqt>=Y & pafqt<=Y+20
\ scalar y00n=r(N)
\ matt2 Y;
clear;
svmat BN;
save ncomply80, replace;

use ncomply80, clear;
drop if BN1==0;
rename BN1 pafqt;
replace pafqt=pafqt+10;
label variable pafqt "AFQT centile";

rename BN2 n11;
rename BN3 n10;
rename BN4 n01;
rename BN5 n00;

gen d0=n01/(n00+n01);
gen d1=n11/(n11+n10);
gen cc=d1-d0;
gen aa=d0;
gen nn=1-d1;
gen cc1=n11/(n00+n10+n01+n11)*(1-aa);
gen cc0=1-aa-nn-cc1;

twoway line cc0 cc1 aa nn pafqt, lpattern(solid shortdash longdash dash_dot) xlabel(10[20]90) ylabel(0.18[0.02]0.28) legend(lab(1 "compliers, serve=0") lab(2 "compliers, serve=1") lab(3 "always takers") lab(4 "never takers") order(1 2 4 3)) ytitle("sample proportion") xtitle("AFQT score centile +/-10") graphregion(color(white));
graph save figa1, replace;
graph export figa1.pdf, replace;



set seed 1234;
prog drop _all;


*Figure A2 - understanding OLS-IV differences by AFQT;

use age learn pnr drawlo served afqt qafqt if age>=25 & age<=35 & learn!=. using "opp-cost-master", clear;
rename drawlo drafted;
keep pnr drafted served learn afqt qafqt;
save compliers, replace;
by pnr, sort: gen nvals=_n==1;
count if nvals;
scalar npnr=r(N);

tab drafted served, su(learn);
tab drafted, su(served);



prog define compliertypep;
  use if pafqt>=`1' & pafqt<=`1'+20 using compliersboot, clear;
  su served if drafted==0;
  scalar d0=r(mean);
  su served if drafted==1;
  scalar d1=r(mean);
  scalar cc=d1-d0;
  scalar aa=d0;
  scalar nn=1-d1;
  su learn if drafted==1 & served==1;
  scalar y11=r(mean);
  su learn if drafted==1 & served==0;
  scalar y10=r(mean);
  su learn if drafted==0 & served==1;
  scalar y01=r(mean);
  su learn if drafted==0 & served==0;
  scalar y00=r(mean);
  su learn if served==1;
  scalar yx1=r(mean);
  su learn if served==0;
  scalar yx0=r(mean);
  su learn if drafted==1;
  scalar y1x=r(mean);
  su learn if drafted==0;
  scalar y0x=r(mean);
  scalar yc1=(d1*y11-d0*y01)/(d1-d0);
  scalar yc0=((1-d0)*y00-(1-d1)*y10)/(d1-d0);
  scalar ols=yx1-yx0;
  scalar iv=(y1x-y0x)/(d1-d0);
  scalar ynn=y10;
  scalar yaa=y01;
  scalar yc1aa=yc1-yaa;
  scalar yc0nn=yc0-ynn;
  for X in any y11 y10 y01 y00 yx1 yx0 y1x y0x ols iv yc1 yc0 ynn yaa yc1aa yc0nn: display "X " X;
  matrix B=B\(y11,y10,y01,y00,yx1,yx0,y1x,y0x,ols,iv,yc1,yc0,ynn,yaa,yc1aa,yc0nn);
end;


prog define times80;
  matrix B=(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
  for X in num 1/80: compliertypep X;
  clear;
  svmat B;
end;


for Y in num 1/100:
  display Y
\ use compliers, clear
\ bsample npnr, cluster(pnr)
\ drop if afqt==. | afqt<28
\ xtile pafqt=afqt, n(100)
\ save compliersboot, replace
\ times80
\ save ddY, replace;


for Y in num 1/100:
  use ddY, clear
\ gen nn=_n
\ drop if nn==1
\ replace nn=nn-1
\ gen brep=Y
\ save dddY, replace;
for Y in num 1/99: append using dddY;
save dddall, replace;


use dddall, clear;
for Y in num 1/16: 
  egen p05BY=pctile(BY), p(5) by(nn)
\ egen p95BY=pctile(BY), p(95) by(nn)
\ egen meaBY=mean(BY), by(nn);

keep if brep==1;
replace nn=nn+10;
label variable nn "AFQT score centile";

twoway rarea p05B16 p95B16 nn, color(gs8) 
  || rarea p05B15 p95B15 nn, color(gs10) xlabel(10[20]90) ylabel(-0.15[0.05]0.1, format(%5.2f)) ytitle("") 
  || line meaB16 nn, lpattern(solid) 
  || line meaB15 nn, lpattern(dash) legend(lab(3 "serve=0 (comply minus never)") lab(4 "serve=1 (comply minus always)") order(3 4) cols(1)) graphregion(color(white));
graph save b1516, replace;
twoway rarea p05B9 p95B9 nn, color(gs8) 
  || rarea p05B10 p95B10 nn, color(gs10) xlabel(10[20]90) ylabel(-0.15[0.05]0.1, format(%5.2f)) ytitle("log earnings difference") 
  || line meaB9 nn, lpattern(solid) 
  || line meaB10 nn, lpattern(dash) legend(lab(3 "OLS") lab(4 "IV") order(3 4) cols(1)) graphregion(color(white));
graph save b0910, replace;
graph combine b0910.gph b1516.gph,  graphregion(color(white));
graph save figa2, replace;
graph export figa2.pdf, replace;



prog drop _all;

*figure A3 - effect of military service on earnings by age;

prog define keepmat;
  matrix C=e(b);
  matrix V=e(V);
  scalar v=sqrt(V[1,1]);
  matrix B=B\(C[1,1],v);
end;


prog define matdat;
  clear;
  svmat B;
  drop if _n==1;
  gen nn=_n;
  gen page=nn+17;
  gen B1hi=B1+2*B2;
  gen B1lo=B1-2*B2;
  keep page B1*;
  matrix B=(0,0);
end;


prog define byearnageafqt;
use pnr afqt using "opp-cost-master", clear;
duplicates drop pnr, force;
for X in any afqt: xtile qqX=X, n(4) \ tab qqX;
compress;
keep pnr qq*;
save quartilesonly, replace;
joinby pnr using "opp-cost-master";
keep learn qqafqt served drawlo height height2 afqt afqt2 yob mob year shy sessan_year age;

*drop if (spaer_year==(year-1))|(spaer_year==year);
keep if learn!=.;
keep if qqafqt>=`1' & qqafqt<=`2';
matrix B=(0,0);
for A in num 19/35:
  xi: ivreg2 learn (served=drawlo) height* afqt* i.yob i.mob i.year i.shy i.sessan_year if age==A
  \ keepmat;
matdat;
label variable page "age";
save matdat_servedq`1'`2', replace;
end;


for Y in num 1 1 4 \ Z in num 4 1 4: byearnageafqt Y Z; 

use matdat_servedq14, clear;
renpfix B1 AB1;
joinby page using matdat_servedq11;
renpfix B1 BB1;
joinby page using matdat_servedq44;
renpfix B1 CB1;
label variable page "age";
replace page=page+1;

twoway rarea AB1hi AB1lo page, color(gs10) xlabel(19[2]35) ytitle("Service coefficient") ylabel(, format(%3.1f))|| line AB1 page, legend(off) graphregion(color(white));
graph save AB, replace;
twoway rarea CB1hi CB1lo page, color(gs10) xlabel(19[2]35) ylabel(, format(%3.1f)) || line CB1 page, legend(off) graphregion(color(white));
graph save CB, replace;
graph combine AB.gph CB.gph, ycom graphregion(color(white));
graph save figa3, replace;
graph export figa3.pdf, replace;







