Monte Carlo and quasi-Monte Carlo simulations for credit risk

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:      hchao8@gmail.com
*
****************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