Sunday, May 29, 2011

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

2 comments:

  1. So good i like your Blogspot and i have blog to talking about Perte de Poids and a lot of thing related to Perte de Poids as aide à la perte de poids, comment maigrir en se musclant, perte de poids inexpliquée, perte de poids rapide, perte de poids anormale, forum perte de poids, programme perte de poids, perdre du poids and thanks a lot again admin ,,, Weight-loss-steps.info

    ReplyDelete
  2. can someone explain what is the weight
    thank you

    ReplyDelete