/*GAUSS CODE TO GENERATE PREDICTED PROBABILITIES & CONF. INTERVALS*/
new;
library pgraph;

/*SELECTION EQUATION*/ 
/*Coefficients from selection equation (gammas)*/
gammas={4.773204, 0.0493991, 0.0139687, -0.2505272, 0.2440559, -0.4550545, 0.0058177, 0.0000483, 0.0003156, -0.3368851, -0.1112283, 0.0448615, -0.2227449, -0.0789882, -3.822542};

/*Selection equation variance-covariance matrix*/
vcselect={0.195390000	0.004188000	-0.000031000	-0.005970000	0.034829000	-0.042573000	0.000021000	-0.000000180	-0.000239000	0.034058000	0.005108000	0.003583000	-0.016048000	-0.005954000	-0.085818000,
0.004188000	0.000173000	0.000002000	0.000305000	0.000451000	-0.001704000	-0.000003300	0.000000006	-0.000005000	0.001609000	0.000131000	0.000065000	-0.000344000	-0.000122000	-0.004337000,
-0.000031000	0.000002000	0.000006400	0.000017000	-0.000053000	0.000080000	0.000000300	-0.000000017	-0.000000150	-0.000455000	0.000006400	-0.000001400	-0.000017000	-0.000015000	-0.000099000,
-0.005970000	0.000305000	0.000017000	0.021319000	-0.005478000	-0.003681000	-0.000030000	-0.000000280	0.000126000	0.012342000	0.001630000	0.000105000	-0.000880000	0.000375000	-0.007919000,
0.034829000	0.000451000	-0.000053000	-0.005478000	0.032800000	-0.007388000	0.000071000	-0.000000440	-0.000010000	-0.010378000	0.004159000	0.000407000	-0.002244000	-0.001644000	0.002524000,
-0.042573000	-0.001704000	0.000080000	-0.003681000	-0.007388000	0.042789000	-0.000064000	-0.000000076	0.000329000	-0.007288000	-0.001455000	-0.000793000	0.002493000	0.000505000	0.038218000,
0.000021000	-0.000003300	0.000000300	-0.000030000	0.000071000	-0.000064000	0.000005000	-0.000000018	0.000003700	-0.000692000	-0.000003300	0.000004400	0.000004100	0.000000490	-0.000252000,
-0.000000180	0.000000006	-0.000000017	-0.000000280	-0.000000440	-0.000000076	-0.000000018	0.000000001	0.000000010	0.000001400	-0.000000088	-0.000000042	-0.000000010	0.000000070	0.000001000,
-0.000239000	-0.000005000	-0.000000150	0.000126000	-0.000010000	0.000329000	0.000003700	0.000000010	0.000124000	-0.001348000	0.000038000	0.000003600	0.000015000	0.000003900	-0.000848000,
0.034058000	0.001609000	-0.000455000	0.012342000	-0.010378000	-0.007288000	-0.000692000	0.000001400	-0.001348000	0.324147000	-0.026902000	0.001381000	-0.004762000	-0.001427000	-0.106899000,
0.005108000	0.000131000	0.000006400	0.001630000	0.004159000	-0.001455000	-0.000003300	-0.000000088	0.000038000	-0.026902000	0.011888000	0.000046000	-0.000459000	-0.000308000	0.038513000,
0.003583000	0.000065000	-0.000001400	0.000105000	0.000407000	-0.000793000	0.000004400	-0.000000042	0.000003600	0.001381000	0.000046000	0.000208000	-0.000698000	-0.000306000	-0.000852000,
-0.016048000	-0.000344000	-0.000017000	-0.000880000	-0.002244000	0.002493000	0.000004100	-0.000000010	0.000015000	-0.004762000	-0.000459000	-0.000698000	0.003519000	0.001085000	-0.002719000,
-0.005954000	-0.000122000	-0.000015000	0.000375000	-0.001644000	0.000505000	0.000000490	0.000000070	0.000003900	-0.001427000	-0.000308000	-0.000306000	0.001085000	0.000777000	0.001143000,
-0.085818000	-0.004337000	-0.000099000	-0.007919000	0.002524000	0.038218000	-0.000252000	0.000001000	-0.000848000	-0.106899000	0.038513000	-0.000852000	-0.002719000	0.001143000	0.369030000};

/*Take 1000 random multivariate normal draws for simulations*/
selectsim=rndmn(gammas,vcselect,1000);                  

/*Mean values of selection equation covariates, non-signatories (w^Ns)*/
wN={0.00 38.94177 22.23774 0.6068884 0.2042755 0.6437055 66.35726 1216.923 3.953453 0.1062854 -3.204889  21.46081 9.425178 1.150831 1.00};

/*NON-SIGNATORIES' OUTCOME EQUATION*/
/*Non-signatory outcome equation coeffcients (beta^Ns)*/ 
betaN={0.1532, -0.0111, 0.6284, -0.0111, 0.2764, -0.0358, 2.8077, 0.5649, -1.8870};
/*Non-signatory outcome equation variance-covariance matrix*/  	
vcnonsigoutcome={0.005856  -0.000038  0.002710  0.000077  0.001175  0.000043  0.000587  0.000569  -0.020084,
      -0.000038   0.000038 -0.000434 -0.000005  0.000016 -0.000022 -0.000212 -0.000219	0.000551,
       0.002710	-0.000434  0.078791  0.000162  0.012219  0.000240  0.001008  0.000882	-0.028606,
       0.000077	-0.000005  0.000162  0.000069  0.000150 -0.000005 -0.000074 -0.000051	-0.000627,
       0.001175	 0.000016  0.012219  0.000150  0.012232 -0.000253 -0.000727  0.000629	-0.012508,
       0.000043	-0.000022  0.000240 -0.000005 -0.000253  0.000364  0.002708  0.002332	-0.002769,
       0.000587	-0.000212  0.001008 -0.000074 -0.000727  0.002708  0.034749  0.028685	-0.032804,
       0.000569	-0.000219  0.000882 -0.000051  0.000629  0.002332  0.028685  0.061392	-0.031227,
      -0.020084	 0.000551 -0.028606 -0.000627 -0.012508 -0.002769 -0.032804 -0.031227	 0.110457};

/*Take 1000 random multivariate normal draws for simulations*/
nonsigsim=rndmn(betaN,vcnonsigoutcome,1000);                   
/*Values of outcome equation covariates, non-signatories (x^Ns) (mean values for all variables except lastrest, lr0, lr1 [I vary them below] and constant [set at 1])*/
/*0 years since last restriction*/
xNzero={3.136 -4.442 .106 3.953 .607 0 1 0 1};    
/*1 year since last restriction*/
xNone={3.136 -4.442 .106 3.953 .607 1 0 1 1};
/*2 years since last restriction*/    
xNtwo={3.136 -4.442 .106 3.953 .607 2 0 0 1};
/*3 years since last restriction*/
xNthree={3.136 -4.442 .106 3.953 .607 3 0 0 1};
/*4 years since last restriction*/
xNfour={3.136 -4.442 .106 3.953 .607 4 0 0 1};
/*5 years since last restriction*/
xNfive={3.136 -4.442 .106 3.953 .607 5 0 0 1};
/*6 years since last restriction*/
xNsix={3.136 -4.442 .106 3.953 .607 6 0 0 1};

/*Multiply w^Ns by selection equation simulation*/
selectionpred=wN*selectsim';
/*Multiply x^Ns (holding at means; varying values of lastrest [and lr0, lr1]) by non-signatory outcome equation simulations*/
nonsigpredzero=xNzero*nonsigsim';
nonsigpredone=xNone*nonsigsim';
nonsigpredtwo=xNtwo*nonsigsim';
nonsigpredthree=xNthree*nonsigsim';
nonsigpredfour=xNfour*nonsigsim';
nonsigpredfive=xNfive*nonsigsim';
nonsigpredsix=xNsix*nonsigsim';

/*generate bivariate normal cdf of (-w^Ngamma, x^Nbeta^N, -rho^N),  at various values of lastrest (and lr0, lr1)*/
nonsigbinormzero=cdfbvn((-1*selectionpred), nonsigpredzero, .334);                         
nonsigbinormone=cdfbvn((-1*selectionpred), nonsigpredone, .334);            
nonsigbinormtwo=cdfbvn((-1*selectionpred), nonsigpredtwo, .334);  
nonsigbinormthree=cdfbvn((-1*selectionpred), nonsigpredthree, .334);
nonsigbinormfour=cdfbvn((-1*selectionpred), nonsigpredfour, .334);   
nonsigbinormfive=cdfbvn((-1*selectionpred), nonsigpredfive, .334); 
nonsigbinormsix=cdfbvn((-1*selectionpred), nonsigpredsix, .334);       

/*generate cumulative standard normal distribution of "non-selection" equation (-w^Ngamma)*/
norm=cdfn(-1*(selectionpred));

/*Divide cdfbvn by cdfn at each value of lastrest*/
nonsigprobzero= nonsigbinormzero./norm; 
nonsigprobone= nonsigbinormone./norm; 
nonsigprobtwo= nonsigbinormtwo./norm; 
nonsigprobthree= nonsigbinormthree./norm; 
nonsigprobfour= nonsigbinormfour./norm; 
nonsigprobfive= nonsigbinormfive./norm; 
nonsigprobsix= nonsigbinormsix./norm; 

/*Invert matrices*/
nonsigprobzero2=nonsigprobzero';
nonsigprobone2=nonsigprobone';
nonsigprobtwo2=nonsigprobtwo';
nonsigprobthree2=nonsigprobthree';
nonsigprobfour2=nonsigprobfour';
nonsigprobfive2=nonsigprobfive';
nonsigprobsix2=nonsigprobsix';

/*Sort*/
nonsigprobzero2=sortc(nonsigprobzero2,1);
nonsigprobone2=sortc(nonsigprobone2,1);
nonsigprobtwo2=sortc(nonsigprobtwo2,1);
nonsigprobthree2=sortc(nonsigprobthree2,1);
nonsigprobfour2=sortc(nonsigprobfour2,1);
nonsigprobfive2=sortc(nonsigprobfive2,1);
nonsigprobsix2=sortc(nonsigprobsix2,1);

/*Generate mean probability of restrictions at different values of lastrest (and lr0, lr1), holding other variables at means*/
nonsigmeanzero=meanc(nonsigprobzero2);
nonsigmeanone=meanc(nonsigprobone2);
nonsigmeantwo=meanc(nonsigprobtwo2);
nonsigmeanthree=meanc(nonsigprobthree2);
nonsigmeanfour=meanc(nonsigprobfour2);
nonsigmeanfive=meanc(nonsigprobfive2);
nonsigmeansix=meanc(nonsigprobsix2);

/*Generate lower-bound confidence intervals at different values of lastrest (and lr0, lr1)*/
nonsiglbzero=nonsigprobzero2[25,.];
nonsiglbone=nonsigprobone2[25,.];
nonsiglbtwo=nonsigprobtwo2[25,.];
nonsiglbthree=nonsigprobthree2[25,.];
nonsiglbfour=nonsigprobfour2[25,.];
nonsiglbfive=nonsigprobfive2[25,.];
nonsiglbsix=nonsigprobsix2[25,.];

/*Generate upper-bound confidence intervals at different values of lastrest (and lr0, lr1)*/
nonsigubzero=nonsigprobzero2[975,.];
nonsigubone=nonsigprobone2[975,.];
nonsigubtwo=nonsigprobtwo2[975,.];
nonsigubthree=nonsigprobthree2[975,.];
nonsigubfour=nonsigprobfour2[975,.];
nonsigubfive=nonsigprobfive2[975,.];
nonsigubsix=nonsigprobsix2[975,.];

/*Print results!*/
print "NON-SIGNATORY PREDICTED PROBABILITIES AS FUNCTION OF YEARS SINCE LAST RESTRICTION:";
print "0" nonsigmeanzero; 
print "1" nonsigmeanone;
print "2" nonsigmeantwo;
print "3" nonsigmeanthree;
print "4" nonsigmeanfour;
print "5" nonsigmeanfive;
print "6" nonsigmeansix;

print "NON-SIGNATORY LOWER BOUND CONF. INTERVALS AS FUNCTION OF YEARS SINCE LAST RESTRICTION:";
print "0" nonsiglbzero;
print "1" nonsiglbone;
print "2" nonsiglbtwo;
print "3" nonsiglbthree;
print "4" nonsiglbfour;
print "5" nonsiglbfive;
print "6" nonsiglbsix;

print "NON-SIGNATORY UPPER BOUND CONF. INTERVALS AS FUNCTION OF YEARS SINCE LAST RESTRICTION:";
print "0" nonsigubzero;
print "1" nonsigubone;
print "2" nonsigubtwo;
print "3" nonsigubthree;
print "4" nonsigubfour;
print "5" nonsigubfive;
print "6" nonsigubsix;

/*SIGNATORIES' OUTCOME EQUATION*/
/*Signatory outcome equation coeffcients (beta^Ss)*/ 
betaS={0.2618, -0.0019, -0.5229, -0.0219, 0.5113, -0.0365, 2.2672, 0.1483, -2.2291};

/*Signatory outcome equation variance-covariance matrix*/
vcsigoutcome={0.006700 0.000021 -0.004415 -0.000042 -0.001897 0.000030 -0.000048 -0.001442 -0.019038,
       0.000021 0.000010	0.000233	-0.000014  0.000056 0.000011  0.000138  0.000052 -0.000126,
      -0.004415 0.000233  0.369003 -0.000102  0.007011 0.000478  0.013359  0.001496 -0.031840,
      -0.000042 -0.000014 -0.000102 0.000119 -0.000045 0.000013 -0.000040  0.000053 -0.000306,
      -0.001897  0.000056	 0.007011 -0.000045 0.014125 0.000245  0.001746  0.002275 -0.004550,
       0.000030	0.000011  0.000478  0.000013 0.000245 0.000250  0.002198  0.001984 -0.002490,
      -0.000048	0.000138	 0.013359 -0.000040 0.001746 0.002198  0.033998  0.024816 -0.028245,
      -0.001442  0.000052  0.001496  0.000053 0.002275 0.001984  0.024816  0.083863 -0.021723,
      -0.019038 -0.000126 -0.031840 -0.000306 -0.004550 -0.002490 -0.028245 -0.021723 0.093614};

/*Take 1000 random multivariate normal draws for simulations*/
sigsim=rndmn(betaS,vcsigoutcome,1000);     

/*Multiply x^Ns (holding at means; varying values of lastrest [and lr0, lr1]) by signatory outcome equation simulations*/
sigpredzero=xNzero*sigsim';
sigpredone=xNone*sigsim';
sigpredtwo=xNtwo*sigsim';
sigpredthree=xNthree*sigsim';
sigpredfour=xNfour*sigsim';
sigpredfive=xNfive*sigsim';
sigpredsix=xNsix*sigsim';

/*generate bivariate normal cdf of (-w^Ngamma, x^Nbeta^S, -rho^S) at various values of lastrest (and lr0, lr1)*/
sigbinormzero=cdfbvn((-1*selectionpred), sigpredzero, .536);                         
sigbinormone=cdfbvn((-1*selectionpred), sigpredone, .536);            
sigbinormtwo=cdfbvn((-1*selectionpred), sigpredtwo, .536);  
sigbinormthree=cdfbvn((-1*selectionpred), sigpredthree, .536);
sigbinormfour=cdfbvn((-1*selectionpred), sigpredfour, .536);   
sigbinormfive=cdfbvn((-1*selectionpred), sigpredfive, .536); 
sigbinormsix=cdfbvn((-1*selectionpred), sigpredsix, .536);       

/*Divide cdfbvn by cdfn at each value of lastrest (and lr0, lr1)*/
sigprobzero= sigbinormzero./norm; 
sigprobone= sigbinormone./norm; 
sigprobtwo= sigbinormtwo./norm; 
sigprobthree= sigbinormthree./norm; 
sigprobfour= sigbinormfour./norm; 
sigprobfive= sigbinormfive./norm; 
sigprobsix= sigbinormsix./norm;

/*Invert matrices*/
sigprobzero2=sigprobzero';
sigprobone2=sigprobone';
sigprobtwo2=sigprobtwo';
sigprobthree2=sigprobthree';
sigprobfour2=sigprobfour';
sigprobfive2=sigprobfive';
sigprobsix2=sigprobsix';

/*Sort*/
sigprobzero2=sortc(sigprobzero2,1);
sigprobone2=sortc(sigprobone2,1);
sigprobtwo2=sortc(sigprobtwo2,1);
sigprobthree2=sortc(sigprobthree2,1);
sigprobfour2=sortc(sigprobfour2,1);
sigprobfive2=sortc(sigprobfive2,1);
sigprobsix2=sortc(sigprobsix2,1);

/*Generate mean probability of restrictions at different values of lastrest (and lr0, lr1), holding other variables at means*/
sigmeanzero=meanc(sigprobzero2);
sigmeanone=meanc(sigprobone2);
sigmeantwo=meanc(sigprobtwo2);
sigmeanthree=meanc(sigprobthree2);
sigmeanfour=meanc(sigprobfour2);
sigmeanfive=meanc(sigprobfive2);
sigmeansix=meanc(sigprobsix2);

/*Generate lower-bound confidence intervals at different values of lastrest (and lr0, lr1)*/
siglbzero=sigprobzero2[25,.];
siglbone=sigprobone2[25,.];
siglbtwo=sigprobtwo2[25,.];
siglbthree=sigprobthree2[25,.];
siglbfour=sigprobfour2[25,.];
siglbfive=sigprobfive2[25,.];
siglbsix=sigprobsix2[25,.];

/*Generate upper-bound confidence intervals at different values of lastrest (and lr0, lr1)*/
sigubzero=sigprobzero2[975,.];
sigubone=sigprobone2[975,.];
sigubtwo=sigprobtwo2[975,.];
sigubthree=sigprobthree2[975,.];
sigubfour=sigprobfour2[975,.];
sigubfive=sigprobfive2[975,.];
sigubsix=sigprobsix2[975,.];

/*Print results!*/
print "SIGNATORY PREDICTED PROBABILITIES AS FUNCTION OF YEARS SINCE LAST RESTRICTION:";
print "0" sigmeanzero;
print "1" sigmeanone;
print "2" sigmeantwo;
print "3" sigmeanthree;
print "4" sigmeanfour;
print "5" sigmeanfive;
print "6" sigmeansix;

print "SIGNATORY LOWER BOUND CONF. INTERVALS AS FUNCTION OF YEARS SINCE LAST RESTRICTION:";
print "0" siglbzero; 
print "1" siglbone; 
print "2" siglbtwo; 
print "3" siglbthree; 
print "4" siglbfour; 
print "5" siglbfive; 
print "6" siglbsix;

print "SIGNATORY UPPER BOUND CONF. INTERVALS AS FUNCTION OF YEARS SINCE LAST RESTRICTION:";
print "0" sigubzero;
print "1" sigubone;
print "2" sigubtwo; 
print "3" sigubthree; 
print "4" sigubfour;
print "5" sigubfive;
print "6" sigubsix;
end;