Hi, I would like to be able to calculate the CDF of the Bivariate (and multivariate) T distribution. Has anyone got a reliable and accurate way to do this in SAS? I suspect IML would be needed hence posting here. R has a package mvtnorm that can do it, and I think it is based methods in the paper by Genz & Bretz (attached). Rather than try to digest what looks to be very complicated separation algebra in the paper, I hoped someone may have a solution already? A very simple way to do it would be a MC integration algorithm; generating large numbers of MVT values from IML and averaging them. I have tried this for the bivariate T (see embedded SAS program) but I don’t think this would be accurate and quick enough (I compared the accuracy of the process for the MVN in the code with the SAS function PROBBNRM and it isn’t great). Any ideas or solutions would be welcome! Thanks, Tom. Using SAS 9.3 and 9.4. **********************************************;
*** TEST BASIC MONTE CARLO INEGRATION ***;
**********************************************;
ods html close;
ods listing;
*** MACRO TO GENERATE MONTE CARLO INTEGRATED MVT CDF ESTIMATE IN IML ***;
%macro mc_mvt;
proc iml;
call randseed(&seed.);
mvector = {0 0};
smatrix = {1 &rho., &rho. 1};
yvector = randmvt(10000, &df., mvector, smatrix);
ind = (yvector[,1] <= &x1.) & (yvector[,2] <= &x2.);
probt = ind[:,];
call symput ("probt",char(probt));
quit;
run;
%mend;
*** MACRO TO GENERATE MONTE CARLO INTEGRATED MVN CDF ESTIMATE IN IML ***;
%macro mc_mvn;
proc iml;
call randseed(&seed.);
mvector = {0 0};
smatrix = {1 &rho., &rho. 1};
yvector = randnormal(10000, mvector, smatrix);
ind = (yvector[,1] <= &x1.) & (yvector[,2] <= &x2.);
probn = ind[:,];
call symput ("probn",char(probn));
quit;
run;
%mend;
*** FCMP FUNCTION TO RUN MACRO FROM WITHIN DATASTEP ***;
proc fcmp outlib = work.stats.functions;
function cdf_mvt(x1,x2,rho,df,seed);
rc = run_macro("mc_mvt",x1,x2,rho,df,seed,probt);
if rc = 0 then return(probt);
else return(.);
endsub;
quit;
run;
*** FCMP FUNCTION TO RUN MACRO FROM WITHIN DATASTEP ***;
proc fcmp outlib = work.stats.functions;
function cdf_mvn(x1,x2,rho,seed);
rc = run_macro("mc_mvn",x1,x2,rho,seed,probn);
if rc = 0 then return(probn);
else return(.);
endsub;
quit;
run;
*** CALL FUNCTIONS FROM DATASTEP ***;
options cmplib = (work.stats);
data results;
do i = -1 to 1 by 0.25;
do j = -1 to 1 by 0.25;
r=0.8;
v=2;
seeds=12348765;
mc_probt = cdf_mvt(i,j,r,v,seeds);
mc_probn = cdf_mvn(i,j,r,seeds);
sas_probn = probbnrm(i,j,r);
output;
end;
end;
run;
... View more
