#delimit;
set more off;

hetprob voted2004 yrseduc age GWBtherm, het(IraqCost);

matrix params = e(b);
matrix CV = e(V);
drawnorm b1-b5, means(params) cov(CV) double;

/* predicted probabilities w/ 1000 draws */
/* hypothetical case: yrseduc=12, age=40, GWBtherm=50 */
/* vary IraqCost from 0 to 1 */


/* IraqCost = 0 */

quietly generate pvIC0 = normprob((b1*12 + b2*40 + b3*50 + b4)/exp(b5*0));

/* IraqCost = 1 */

quietly generate pvIC1 = normprob((b1*12 + b2*40 + b3*50 + b4)/exp(b5*1));

quietly generate diffI10 = pvIC1 - pvIC0;

summarize pvIC0 pvIC1 diffI10;

drop pvIC0 pvIC1 diffI10;

/* create a graph across yrseduc */

generate pvoteIC0 = .;
generate pvoteIC1 = .;
generate IC0lo = .;
generate IC0hi = .;
generate IC1lo = .;
generate IC1hi = .;
generate gryred = .;

local i=1;

while `i' < 18 {;

local a=1;

quietly generate pvIC0d = .;
quietly generate pvIC1d = .;

  while `a' < 1001{;

  quietly replace pvIC0d = normprob((b1[`a']*`i' + b2[`a']*40 + 
  b3[`a']*50 + b4[`a'])/exp(b5[`a']*0)) if _n==`a';

  quietly replace pvIC1d = normprob((b1[`a']*`i' + b2[`a']*40 + 
  b3[`a']*50 + b4[`a'])/exp(b5[`a']*1)) if _n==`a';

  local a = `a' + 1;
  };

quietly egen mpvIC0 = mean(pvIC0d);
quietly replace pvoteIC0 = mpvIC0 if _n==`i';
quietly egen mpvIC1 = mean(pvIC1d);
quietly replace pvoteIC1 = mpvIC1 if _n==`i';
_pctile pvIC0d, p(2.5, 97.5);
quietly replace IC0lo = r(r1) if _n==`i';
quietly replace IC0hi = r(r2) if _n==`i';
_pctile pvIC1d, p(2.5, 97.5);
quietly replace IC1lo = r(r1) if _n==`i';
quietly replace IC1hi = r(r2) if _n==`i';

quietly replace gryred = `i' if _n== `i';

quietly drop pvIC0d pvIC1d mpvIC0 mpvIC1;

local i = `i' + 1;

};

label var pvoteIC0 "Iraq Cost = 0";
label var pvoteIC1 "Iraq Cost = 1";
label var IC0lo "IC = 0 95% lower limit";
label var IC0hi "IC = 0 95% upper limit";
label var IC1lo "IC = 1 95% lower limit";
label var IC1hi "IC = 1 95% upper limit";
label var gryred "Years of education";

graph twoway line pvoteIC0 pvoteIC1 IC0lo IC1lo IC0hi IC1hi gryred, lwidth(medthick medthick thin thin thin thin) 
lcolor(black black black black black black) lpattern(solid dash dot dot dot dot) ytitle("Probability of voting");

drop b1 b2 b3 b4 b5 pvoteIC0 pvoteIC1 IC0lo IC0hi IC1lo IC1hi gryred;