In SAS, we can do quasi-Monte Carlo simulation by Halton’s low-discrepancy sequences, which are better than random variables from uniform distribution . The only procedure in SAS 9.2 to produce Halton’s sequence is probably the MDC procedure in SAS/ETS [Ref. 2] but hard to output. However, we can build a user-defined function by PROC FCMP to realize our goals. I showed an example by modifying the codes from Paolo Brandimarte’s Matlab book [Ref. 1] to estimate the area of an arbitrary surface. Overall, Halton’s method would converge with more sampling points, while uniform randomization will cause estimates fluctuating along the real value.

/*******************READ ME*********************************************
* - QUASI-MONTE CARLO SIMULATION BY FUNCTIONAL PROGRAMMING -
*
* VERSION: SAS 9.2(ts2m0), windows 64bit
* DATE: 02apr2011
* AUTHOR: [email protected]
*
****************END OF READ ME*****************************************/

****************(1) MODULE-BUILDING STEP******************;
******(1.1) COMPILE USER-DEFINED FUNCTIONS*************;
proc fcmp outlib = sasuser.qmc.funcs;
/***********************************************************
* FUNCTION: h = halton(n, b)
* INPUT: n = the nth number | b = the base number
* OUTPUT: h = the halton's low-discrepancy number
***********************************************************/
function halton(n, b);
n0 = n;
h = 0;
f = 1 / b;
do while (n0 > 0);
n1 = floor(n0 / b);
r = n0 - n1*b;
h = h + f*r;
f = f / b;
n0 = n1;
end;
return(h);
endsub;

/***********************************************************
* FUNCTION: z = surface1(x, y)
* INPUT: x, y = independent variables
* OUTPUT: z = response variable
***********************************************************/
function surface1(x, y);
pi = constant("PI");
z = exp((-x)*y) * (sin(6*pi*x) + cos(8*pi*y));
return(z);
endsub;
run;

******(1.2) BUILD A 3D PLOTTING MACRO*************;
%macro plot3d(ds = , x = , y = , z = , width = , height = , zangle = );
ods html style = money;
ods graphics / width = &width.px height = &height.px imagefmt = png;
proc template;
define statgraph surfaceplotparm;
begingraph;
layout overlay3d / cube = true rotate = &zangle;
surfaceplotparm x = &x y = &y z = &z
/ surfacetype = fill surfacecolorgradient = &z;
endlayout;
endgraph;
end;
run;

proc sgrender data = &ds
template = surfaceplotparm;
run;
ods graphics off;
ods html close;
%mend;

****************END OF STEP (1)******************;

****************(2) PROGRAMMING STEP**********************;
******(2.1) IMPORT USER-DEFINED FUNCTIONS*************;
option cmplib = (sasuser.qmc);

******(2.2) DISTRIBUTIONS BETWEEN HALTON AND UNIFORM RANDOM NUMBERS*;
data test;
do n = 1 to 100;
do b = 2, 7;
halton_random_number = halton(n, b);
uniform_random_number = ranuni(20110401);
output;
end;
end;
run;

proc transpose data = test out = test1;
by n;
var halton_random_number uniform_random_number;
id b;
run;

******(2.3) GENERATE DATA FOR PLOTTING FUNCTION'S SURFACE****;
data surface;
do x = 0 to 1 by 0.01;
do y = 0 to 1 by 0.01;
z = surface1(x, y);
output;
end;
end;
run;

******(2.4) CONDUCT QUASI-MONTE CARLO SIMULATION***;
data simuds;
do i = 1 to 50;
do n = 1 to 100*i;
do b = 2, 7;
halton_random_number = halton(n, b);
uniform_random_number = ranuni(20110401);
output;
end;
end;
end;
run;

proc transpose data = simuds out = simuds_t;
by i n;
id b;
var halton_random_number uniform_random_number;
run;

data simuds1;
set simuds_t;
x = surface1(_2, _7);
run;

proc sql;
create table simuds2 as
select i, _name_, mean(x) as area label = 'Area by simulation'
from simuds1
group by i, _name_
;quit;

****************END OF STEP (2)******************;

****************(3) VISUALIZATION STEP*******************************;
******(3.1) PLOT DATA BY STEP 2.2*************;
ods html style = money;
proc sgscatter data = test1;
plot _2*_7 / grid group = _name_;
label _2 = 'Sequences with 2 as base'
_7 = 'Sequences with 7 as base'
_name_ = 'Methods';
run;
ods html close;

******(3.2) PLOT DATA BY STEP 2.3*************;
%plot3d(ds = surface, x = x , y = y, z = z, width = 800, height = 800, zangle = 60)

******(3.3) PLOT DATA BY STEP 2.4*************;
ods html style = ocean;
proc sgplot data = simuds2;
series x = i y = area / group = _name_ ;
refline 0.0199 / axis = y label = ('Real area');
xaxis label = 'Random numbers --- ×100';
label _name_ = 'Methods';
run;
ods html close;

****************END OF STEP (3)******************;

****************END OF ALL CODING***************************************;


First, we can use Monte Carlo simulation to evaluate the risk of a portfolio of credit instruments. As for a risk manager, the common question is to answer the probability of the losses for a portfolio over certain time period. To conduct a Monte Carlo simulation, four variables, such as Probability of default(PD), Loss given default(LGD), Exposure at default(EAD), Weight(correlations among individual credit events), are required. Here I simulated a portfolio with 100 obligors. Then I ran a 10000- loop for MC and calculated the loss at given percentiles. The histogram of the portfolio loss is displayed above.

data raw;
format pd lgd percent8.;
w = 0.3;
do i = 1 to 50;
PD = 0.01;
LGD= 0.5;
EAD= 100;
output;
end;
do i = 51 to 100;
PD = 0.05;
LGD= 0.5;
EAD= 100;
output;
end;
label pd = 'Probability of default'
lgd = 'Loss given default'
ead = 'Exposure at default'
w = 'Weight';
run;

libname mem 'c:\tmp' memlib;
data mem.raw;
set raw;
run;

%macro mc(cycle = );
filename junk dummy;
proc printto log=junk; run;
%do i = 1 %to &cycle;
%if %eval(&i) = 1 %then %do;
proc sql;
create table mem.result (sumloss num);
quit;
%end;
data mem._tmp01;
set mem.raw;
retain z;
if _n_ = 1 then z = rannor(0);
d = probit(pd);
a = w*z + (1-w**2)**0.5*rannor(0);
if a < d then loss = lgd*ead;
else loss = 0;
run;
proc sql;
insert into mem.result
select sum(loss) from mem._tmp01;
quit;
%end;
proc univariate data = mem.result noprint;
output out = mem.pct pctlpts = 90 95 99 99.5 pctlpre = loss;
run;
proc printto; run;
%mend;
%mc(cycle = 10000);

References:
1. Paolo Brandimarte. Numerical methods in finance. John Wiley & Sons. 2002.
2. SAS/ETS 9.2 user’s guide. SAS Publishing. 2008