Multicollinearity and the solutions

In his book, Rudolf Freund described a confounding phenomenon while fitting a linear regression. Given a small data set below, there are three variables - dependent variable(y) and independent variables(x1 and x2). Using x2 to fit y alone, the estimated parameter of x2 f is positive that is 0.78. Then using x1 and x2 together to fit y, the parameter of x2 becomes -1.29, which is hard to explain since clearly x2 and y has a positive correlation.

data raw;
input y x1 x2;
cards;
2 0 2
3 2 6
2 2 7
7 2 5
6 4 9
8 4 8
10 4 7
7 6 10
8 6 11
12 6 9
11 8 15
14 8 13
;;;
run;

ods graphics on / border = off;
proc sgplot data = raw;
   reg x = x2 y = y;
   reg x = x2 y = y / group = x1 datalable = x1;
run;

The reason is that x1 and x2 have strong correlation each other. Diagnostics are well when using x2 to fit y. However, counting x1 and x2 together into the regression model causes multicollinearity, and therefore demonstrates severe heteroskedasticity and a skewed distribution of the residuals, which violates the assumptions for OLS regressions. Shown in the top scatter plot, 0.78 is the slope of the regression line by y ~ x2 (the longest straight line), while -1.29 is actually the slope of the partial regression lines by y ~ x2|x1 (four short segments).

proc reg data = raw;
   model y = x2;
   ods select parameterestimates diagnosticspanel;
run;

proc reg data = raw;
   model y = x1 x2;
   ods select parameterestimates diagnosticspanel;
run;

Solutions:

1. Drop a variable
Standing alone, x1 seems like a better predictor (higher R-square and lower MSE) than x2. The easiest way to remove this multicollinearity is to keep only x1 in the model.

proc reg data = raw;
   model y = x1;
   ods select parameterestimates diagnosticspanel;
run;

2. Principle component regression
If we want to keep both variables to avoid information loss, principle component regression is a good option. PCA would transform the correlated variables to the orthogonal factors. In this case, the 1st eigenvector explains 97.77% of the total variance, which is fairly enough for the following regression. SAS's PLS procedure can also perform the principle component regression.

proc princomp data = raw out = pca;
   ods select screeplot corr eigenvalues eigenvectors;
   var x1 x2;
run;

proc reg data = pca;
   model y = prin1;
   ods select parameterestimates diagnosticspanel;
run;

3. Repeated measure

We may apply the mixed model and compute the R matrix using x1 as covariate.

proc mixed data = raw plots(only) = ResidualPanel method = ml;
   model y = x2 / s;
   repeated / type=ar(1) subject=x1 r;
run;

Make a frequency function in SAS/IML

Aggregation is probably the most popular operation in the data world. R comes with a handy table() function. Usually in SAS, the FREQ procedure would deal with this job. It will be great if SAS/IML has an equivalent function. I just created a user-defined function or module for such a purpose. Since it contains a DO loop, the efficiency is not very ideal -- always 10 times slower than PROC FREQ for a simulated data set of one million records.

/* 1 - Use IML for simulation and aggregation */
proc iml;
   start freq(invec);
      x = t(unique(invec));
      y = repeat(x, 1, 2);
      do i = 1 to nrow(x);
         y[i, 2] = ncol(loc(invec=y[i, 1]));
      end;
      return(y);
   finish;
   store module = freq;
quit;

proc iml;
   load module = freq;
   /* Simulate a vector with 1 million values */
   test = abs(floor(rannor(1:1e6)*100));
   t0 = time();
      result = freq(test);
   timer = time() - t0;
   print timer;   
   /* Output the result matrix as SAS data set */
   create a var{"level" "frequency"};
      append from result;
   close a;
quit;

proc sgplot data = a;
   series x = level y = frequency;
run;

/* 2 - Use PROC FREQ for simulation and aggregation */
data test;
   do i = 1 to 1e6;
      test = abs(floor(rannor(1234)*100));
      output;
   end;
   drop i;
run;

options fullstimer;
proc freq data = test noprint;
   table test / out = b;
run;

Top 10 reasons for a modeler to learn PROC FCMP

10.  New financial functions
34 pre-compiled financial functions were shipped with SAS for free, which can be called in DATA step.

9. A management GUI
There is a nice-looking Java-powered GUI to manage the user-defined functions: SAS FCmp function editor.

8. Dynamic parameters
Parameters from other statistical procedures in SAS can be passed to a function. We can use a macro to wrap PROC FCMP to create tons of data-driven functions (temporary or permanent).

7. Migration from VBA
Excel dominates every desktop. However, using VBA to process large data is just a pain, such as calculating a transition matrix for probability of default. PROC FCMP is built on the logic of VBA -- easy to port any Excel macro/functions to SAS.

6. Backtesting
PROC FCMP loads data into memory. For backtesting task, it is much more efficient than hard-disk based DATA step.

5. Function/subroutine compiler
A lot of functions/subroutines SAS doesn’t have can be built, such as Binomial tree function or an asset price simulation subroutine.

4.  Encapsulation
To combine a lot of macros/nested macros for a complicated system, encapsulating the macro variables is a serious concern – PROC FCMP may be a rescue.

3. Monte Carlo simulation
The MCMC procedure, the designated procedure for Monte Carlo simulation in SAS, fully supports the syntax of PROC FCMP with its BEGINCNST and ENDCNST statements.

2. Deal with matrix
As part of SAS/BASE, PROC FCMP, allows all kinds of matrix operations everywhere with a SAS system.

1. Low level computation
The combination of PROC FCMP and PROC PROTO together wraps C functions in SAS. It is interesting to compare the efficiency between a C function and a SAS function.

Rolling regressions for backtesting

Market always generates huge volume time series data with millions of records. Running regressions to obtain the coefficients in a rolling time window is common for many backtesing jobs. In SAS, writing a macro based on the GLM procedures such as PROC REG is not an efficient option. We can imagine the situations: calling PROC REG thousands of times in a big loop would easily petrify any system.

The better way is to go down to the bottom to re-implement the OLS clichés: inverse, transpose and multiply the vectors and matrices. We can do it in either PROC IML, DATA step array or PROC FCMP. For such attempts PROC IML is really powerful but needs extra license. DATA step array would require very high data manipulation skills, since it is not designed for matrix operations. PROC FCMP, a part of SAS/BASE, seems like a portable solution for SAS 9.1 or later. To test this method, I simulated a two-asset portfolio with 100k records, and under a 1000-obs long rolling window, eventually ran 99,001 regressions. The time cost was just 10 seconds on an old laptop. Overall, the speed is quite satisfying.

/* 1 -- Simulate a fund and russell2000 index */
data simuds;
   _beta0 = 0.2;
   _beta1 = 0.6;
   _mse = 0.8;
   format day date9.;
   do day = today() - 20000 to today();
     russell2000 = rannor(1234);
     myfund = _beta0 + _beta1*russell2000 + _mse*rannor(3421);
     output;
   end;
   drop _:;
run;
 
proc print data = simuds(obs = 100);
run;
 
/* 2 -- Decide length of rolling window */
%macro rollreg(data = , wlength = , benchmark = , fund = , out = , rfree = );
   data _1;
      set &data nobs = nobs;
      &benchmark = &benchmark - &rfree;
      &fund = &fund - &rfree;
      call symput('nloop', nobs - &wlength + 1);
      call symput('nobs', nobs);
   run;
   %put &nloop &nobs;
   /* 3 -- Manipulate matrices */
   proc fcmp;
      /* Allocate spaces for matrices */
      array input[&nobs, 2] / nosym;
      array y[&wlength] / nosym;
      array ytrans[1, &wlength] / nosym;
      array xtrans[2, &wlength] / nosym;
      array x[&wlength, 2] / nosym;
      array sscp[2, 2] / nosym;
      array sscpinv[2, 2] / nosym;
      array beta[2] / nosym;
      array result[&nloop, 4] / nosym;
      array beta_xtrans_y[1] / nosym;
      array xtrans_y[2, 1] / nosym;
      array ytrans_y[1] / nosym;
      
      /* Input simulation dataset */
      rc1 = read_array("_1", input, "&benchmark", "&fund");
 
      /* Calculate OLS regression coefficients and r square */
      do j = 1 to &nloop;
         ytotal = 0;
         do i = 1 to &wlength;  
            xtrans[2, i] = input[i+j-1, 1];
            xtrans[1, i] = 1;
            y[i] = input[i+j-1, 2];
          ytotal = ytotal + y[i];
         end;  
         call transpose(y, ytrans);
         call mult(ytrans, y, ytrans_y);
         call transpose(xtrans, x);
         call mult(xtrans, x, sscp);
         call inv(sscp, sscpinv);
         call mult(xtrans, y, xtrans_y);
         call mult(sscpinv, xtrans_y, beta);
         call mult(beta, xtrans_y, beta_xtrans_y);
         nymeansquare = ytotal**2 / &wlength;
         result[j, 1] = beta[1];
         result[j, 2] = beta[2];
         result[j, 3] = (beta_xtrans_y[1]-nymeansquare) / (ytrans_y[1]-nymeansquare);
         result[j, 4] = j;
      end;
 
      /* Output resulting matrix as dataset */
      rc2 = write_array("&out", result, 'beta0', 'beta1', 'rsquare', 'day');
      if rc1 + rc2 > 0 then put 'ERROR: I/O error';
      else put 'NOTE: I/O was successful';
    quit;
%mend;
 
%rollreg(data = simuds, wlength = 50, benchmark = russell2000, fund = myfund, out = result, rfree = 0.005);
 
ods graphics on / noborder;
/* 4 -- Visualize result */
proc sgplot data = result;
   needle x = day y = beta1;
   yaxis label = 'Systematic Risk';
run;
 
proc sgplot data = result;
   needle x = day y = beta0;
   yaxis label = "Jensen's alpha";
run;
 
proc sgplot data = result;
   needle x = day y = rsquare;
   yaxis label = "Coefficient of Determination";
run;

Test A-IRB approach for credit rating

Basel II or Basel III framework would allow qualified financial institutions to apply their own rating systems for credit risks, such as advanced internal ratings-based approach (A-IRB). The equations of required capital play a big role in it. With full-fledged facility of SAS, an automatic selection system can be built by its’ macros, SQL syntax and function compiler. I just ran a simple experiment based on a standard seven-level grading system, toward a data set with a varying PD and the fixed EAD, LGD and maturity.

First 1000 rating structures were simulated. SAS's seeds are used to memorize the structures (Liang has an interesting post regarding the "seed trap" problem).

%macro structgen(data = , var = , out = );
   proc sql noprint;
      create table &out(y num, seed num, step num);
   quit;
   %do i = 1 %to 1000;
      data _tmp1;
         set &data;
         retain y 0;
         y = y + ranuni(&i)*(&var-y);
         seed = &i;
         step = _n_ ;
      run;
      proc sql noprint;
         insert into &out
         select y, seed, step 
         from _tmp1;
      quit;
   %end;
%mend;
%structgen(data = maxbound, var = maxvalue, out = final);

proc sgplot data = final;
   series x = step y = y / group = seed;
   yaxis label = 'Percentage'; xaxis grid;
run;

Second, for each of the rating structure, GINI and required capital were calculated. The min value of required capital turned out to be 0.069, while the max value of GINI is 0.749.

/*Global graphic options*/

ods graphics on / border = off height =4in;
ods html image_dpi=300;

/*Part 2 - Test A-IRB approach for credit rating*/
data maxbound;
   input @1 maxvalue ;
   maxvalue = maxvalue / 100;
   cards;   
0
0.2
0.5
2
10
20
25
;;;
run;

data p;
   input @1 pd ;
   cards;
30.85375387
30.18106395
29.56316569
28.99087964
28.45727718
27.95698093
27.4857216
27.04004492
26.61711101
26.21455288
25.83037412
25.46287322
25.11058664
24.77224523
24.44674027
24.13309679
23.83045233
23.53803976
23.25517337
22.98123742
22.71567665
22.45798844
22.20771616
21.96444353
21.72778996
21.49740641
21.27297191
21.05419056
20.84078894
20.63251378
20.42913001
20.230419
20.03617701
19.84621387
19.66035172
19.47842399
19.30027441
19.12575616
18.95473112
18.78706914
18.62264749
18.46135023
18.30306772
18.1476962
17.99513729
17.84529771
17.69808884
17.55342646
17.41123044
17.27142447
17.13393582
16.9986951
16.86563609
16.73469551
16.60581286
16.47893024
16.25034834
16.02806181
15.81176519
15.60117514
15.39602834
15.19607959
15.00110014
14.81087624
14.62520783
14.44390735
14.26679871
14.09371639
13.92450453
13.75901624
13.59711286
13.4386634
13.28354392
13.13163706
12.98283158
12.83702189
12.69410772
12.55399373
12.4165892
12.28180773
12.14956696
12.01978834
11.89239685
11.76732084
11.6444918
11.52384418
11.40531522
11.28884481
11.1743753
11.06185142
10.95122011
10.8424304
10.73543333
10.63018182
10.52663059
10.42473606
10.32445627
10.22575079
10.12858064
10.03290827
9.938697416
9.845913103
9.754521552
9.664490136
9.575787328
9.488382647
9.40224662
9.317350728
9.233667376
9.151169842
9.069832251
8.989629531
8.910537386
8.83253226
8.75559131
8.679692377
8.604813957
8.530935179
8.458035775
8.386096064
8.315096923
8.245019769
8.175846541
8.107559677
8.040142096
7.973577185
7.907848776
7.842941134
7.778838943
7.715527287
7.65299164
7.59121785
7.530192128
7.469901036
7.410331474
7.351470668
7.293306163
7.235825809
7.179017752
7.122870426
7.067372541
7.012513079
6.958281279
6.904666636
6.851658888
6.799248009
6.747424206
6.696177908
6.645499759
6.595380616
6.545811537
6.496783781
6.448288796
6.40031822
6.352863871
6.305917742
6.259471998
6.213518973
6.16805116
6.12306121
6.078541928
6.034486267
5.990887326
5.947738345
5.9050327
5.862763902
5.820925592
5.779511537
5.738515628
5.697931877
5.657754412
5.617977476
5.578595423
5.539602717
5.500993925
5.462763719
5.424906873
5.387418256
5.350292835
5.313525671
5.277111914
5.241046805
5.205325672
5.169943926
5.134897062
5.100180657
5.065790365
5.03172192
4.997971129
4.964533873
4.931406106
4.898583852
4.866063205
4.833840324
4.801911435
4.77027283
4.738920861
4.707851943
4.677062552
4.646549222
4.616308545
4.586337168
4.556631796
4.527189185
4.498006147
4.469079543
4.440406287
4.411983341
4.383807717
4.355876473
4.328186714
4.300735593
4.273520305
4.246538091
4.219786232
4.193262055
4.166962926
4.140886251
4.115029478
4.089390091
4.063965615
4.038753611
4.013751675
3.988957443
3.964368584
3.9399828
3.915797829
3.891811444
3.868021446
3.844425673
3.821021991
3.797808298
3.774782523
3.751942624
3.729286589
3.706812434
3.684518203
3.662401968
3.640461829
3.618695911
3.597102366
3.575679372
3.554425133
3.533337876
3.512415853
3.491657341
3.47106064
3.450624072
3.430345984
3.410224743
3.390258739
3.370446385
3.350786112
3.331276376
3.311915649
3.292702426
3.273635222
3.254712571
3.235933025
3.217295156
3.198797554
3.180438829
3.162217606
3.14413253
3.126182262
3.108365481
3.090680883
3.073127179
3.055703097
3.038407383
3.021238796
3.004196111
2.98727812
2.970483627
2.953811455
2.937260438
2.920829425
2.904517281
2.888322882
2.87224512
2.8562829
2.840435139
2.824700769
2.809078734
2.79356799
2.778167505
2.762876262
2.747693254
2.732617486
2.717647975
2.70278375
2.688023851
2.673367329
2.658813247
2.644360678
2.630008706
2.615756425
2.601602941
2.587547369
2.573588835
2.559726473
2.54595943
2.532286861
2.518707929
2.50522181
2.491827686
2.47852475
2.465312204
2.452189257
2.439155131
2.426209051
2.413350256
2.40057799
2.387891506
2.375290065
2.362772938
2.350339401
2.337988741
2.32572025
2.31353323
2.301426988
2.28940084
2.277454111
2.265586129
2.253796233
2.242083768
2.230448085
2.218888543
2.207404507
2.195995348
2.184660447
2.173399187
2.16221096
2.151095165
2.140051206
2.129078492
2.118176441
2.107344476
2.096582023
2.085888519
2.075263403
2.06470612
2.054216124
2.04379287
2.033435822
2.023144448
2.012918221
2.002756619
1.992659128
1.982625236
1.972654438
1.962746233
1.952900125
1.943115624
1.933392244
1.923729504
1.914126928
1.904584044
1.895100385
1.885675489
1.876308898
1.867000159
1.857748823
1.848554445
1.839416585
1.830334806
1.821308678
1.812337773
1.803421666
1.794559939
1.785752176
1.776997965
1.768296899
1.759648574
1.75105259
1.742508552
1.734016066
1.725574745
1.717184202
1.708844057
1.700553932
1.692313452
1.684122247
1.675979949
1.667886194
1.659840622
1.651842875
1.643892599
1.635989444
1.628133062
1.620323108
1.612559243
1.604841126
1.597168424
1.589540804
1.581957938
1.5744195
1.566925166
1.559474616
1.552067533
1.544703603
1.537382514
1.532168714
1.526976586
1.521806021
1.516656906
1.511529131
1.506422588
1.501337167
1.496272761
1.491229261
1.486206561
1.481204555
1.476223138
1.471262204
1.46632165
1.461401373
1.456501269
1.451621237
1.446761175
1.441920981
1.437100557
1.432299802
1.427518617
1.422756903
1.418014564
1.413291502
1.408587621
1.403902823
1.399237015
1.394590101
1.389961987
1.385352579
1.380761785
1.376189511
1.371635666
1.367100159
1.362582899
1.358083795
1.353602758
1.349139699
1.344694529
1.340267159
1.335857503
1.331465474
1.327090984
1.322733948
1.318394281
1.314071897
1.309766712
1.305478643
1.301207605
1.296953516
1.292716293
1.288495854
1.284292119
1.280105006
1.275934434
1.271780323
1.267642595
1.263521169
1.259415968
1.255326913
1.251253926
1.247196931
1.24315585
1.239130607
1.235121127
1.231127333
1.227149152
1.223186508
1.219239327
1.215307535
1.21139106
1.207489829
1.203603768
1.199732807
1.195876873
1.192035896
1.188209805
1.184398528
1.180601998
1.176820143
1.173052895
1.169300185
1.165561944
1.161838106
1.158128602
1.154433365
1.150752328
1.147085425
1.14343259
1.139793757
1.136168861
1.132557837
1.128960621
1.125377148
1.121807355
1.118251177
1.114708553
1.111179419
1.107663713
1.104161373
1.100672337
1.097196543
1.093733932
1.090284442
1.086848012
1.083424584
1.080014097
1.076616491
1.073231709
1.069859691
1.066500379
1.063153715
1.059819642
1.056498101
1.053189036
1.049892391
1.046608108
1.043336132
1.040076407
1.036828877
1.033593487
1.030370182
1.027158909
1.023959611
1.020772236
1.01759673
1.014433039
1.01128111
1.008140891
1.005012328
1.001895371
0.998789965
0.995696061
0.992613606
0.98954255
0.986482841
0.983434429
0.980397264
0.977371296
0.974356474
0.97135275
0.968360074
0.965378397
0.96240767
0.959447846
0.956498875
0.95356071
0.950633303
0.947716608
0.944810575
0.94191516
0.939030314
0.936155993
0.933292148
0.930438735
0.927595708
0.924763021
0.921940629
0.919128488
0.916326552
0.913534776
0.910753118
0.907981532
0.905219974
0.902468403
0.899726772
0.896995041
0.894273166
0.891561103
0.888858812
0.886166249
0.883483372
0.880810141
0.878146512
0.875492445
0.872847899
0.870212832
0.867587205
0.864970976
0.862364105
0.859766552
0.857178277
0.854599241
0.852029403
0.849468726
0.846917169
0.844374694
0.841841262
0.839316834
0.836801373
0.834294841
0.831797198
0.829308408
0.826828433
0.824357236
0.82189478
0.819441028
0.816995942
0.814559487
0.812131626
0.809712323
0.807301541
0.804899246
0.8025054
0.80011997
0.797742919
0.795374213
0.793013816
0.790661693
0.788317811
0.785982134
0.783654629
0.781335261
0.779023997
0.776720802
0.774425644
0.772138488
0.769859301
0.767588051
0.765324704
0.763069229
0.760821591
0.758581759
0.7563497
0.754125383
0.751908775
0.749699844
0.74749856
0.74530489
0.743118803
0.740940268
0.738769254
0.73660573
0.734449665
0.73230103
0.730159792
0.728025923
0.725899392
0.723780169
0.721668223
0.719563527
0.717466049
0.715375761
0.713292634
0.711216637
0.709147743
0.707085922
0.705031146
0.702983386
0.700942614
0.698908801
0.696881919
0.694861941
0.692848838
0.690842583
0.688843149
0.686850506
0.684864629
0.68288549
0.680913062
0.678947319
0.676988232
0.675035776
0.673089923
0.671150649
0.669217925
0.667291727
0.665372027
0.663458801
0.661552021
0.659651664
0.657757702
0.65587011
0.653988864
0.652113937
0.650245305
0.648382943
0.646526826
0.64467693
0.642833229
0.640995698
0.639164315
0.637339054
0.635519891
0.633706803
0.631899764
0.630098752
0.628303743
0.626514713
0.624731638
0.622954495
0.621183262
0.619417914
0.617658428
0.615904783
0.614156954
0.612414919
0.610678656
0.608948142
0.607223354
0.60550427
0.603790868
0.602083126
0.600381022
0.598684533
0.596993639
0.595308316
0.593628545
0.591954302
0.590285568
0.58862232
0.586964537
0.585312198
0.583665282
0.582023768
0.580387635
0.578756863
0.57713143
0.575511317
0.573896503
0.572286967
0.570682689
0.569083649
0.567489827
0.565901203
0.564317757
0.562739469
0.561166319
0.559598288
0.558035356
0.556477504
0.554924713
0.553376962
0.551834233
0.550296506
0.548763764
0.547235986
0.545713153
0.544195248
0.542682251
0.541174144
0.539670908
0.538172525
0.536678976
0.535190244
0.533706309
0.532227154
0.53075276
0.529283111
0.527818187
0.526357971
0.524902446
0.523451593
0.522005395
0.520563835
0.519126896
0.517694559
0.516266808
0.514843625
0.513424993
0.512010896
0.510601316
0.509196237
0.507795641
0.506399512
0.505007833
0.503620588
0.502237761
0.500859333
0.49948529
0.498115616
0.496750293
0.495389305
0.494032637
0.492680273
0.491332196
0.489988391
0.488648841
0.487313532
0.485982447
0.484655571
0.483332889
0.482014384
0.480700042
0.479389846
0.478083783
0.476781836
0.475483991
0.474190232
0.472900544
0.471614913
0.470333323
0.46905576
0.467782209
0.466512656
0.465247085
0.463985482
0.462727833
0.461474123
0.460224337
0.458978462
0.457736484
0.456498387
0.455264158
0.454033783
0.452807248
0.451584538
0.450365641
0.449150541
0.447939226
0.446731682
0.445527894
0.44432785
0.443131535
0.441938937
0.440750042
0.439564836
0.438383306
0.437205439
0.436031222
0.434860641
0.433693684
0.432530337
0.431370587
0.430214422
0.429061829
0.427912794
0.426767306
0.425625351
0.424486916
0.423351989
0.422220558
0.42109261
0.419968132
0.418847112
0.417729538
0.416615398
0.415504678
0.414397368
0.413293454
0.412192926
0.41109577
0.410001974
0.408911528
0.407824418
0.406740634
0.405660162
0.404582992
0.403509112
0.40243851
0.401371175
0.400307095
0.399246258
0.398188653
0.397134268
0.396083093
0.395035116
0.393990325
0.39294871
0.391910259
0.390874961
0.389842805
0.38881378
0.387787875
0.386765079
0.385745381
0.38472877
0.383715236
0.382704767
0.381697353
0.380692983
0.379691647
0.378693334
0.377698033
0.376705734
0.375716427
0.3747301
0.373746744
0.372766348
0.371788902
0.370814396
0.369842819
0.368874162
0.367908413
0.366945564
0.365985603
0.365028522
0.364074309
0.363122956
0.362174451
0.361228786
0.36028595
0.359345934
0.358408728
0.357474322
0.356542707
0.355613872
0.354687809
0.353764508
0.352843959
0.351926153
0.35101108
0.350098731
0.349189097
0.348282168
0.347377935
0.346476389
0.345577521
0.344681321
0.34378778
0.342896889
0.34200864
0.341123023
0.340240028
0.339359648
0.338481873
0.337606694
0.336734103
0.335864091
0.334996648
0.334131766
0.333269437
0.332409652
0.331552401
0.330697678
0.329845472
0.328995775
0.328148579
0.327303876
0.326461656
0.325621912
0.324784635
0.323949817
0.32311745
0.322287524
0.321460033
0.320634967
0.319812319
0.318992081
0.318174243
0.317358799
0.31654574
0.315735058
0.314926745
0.314120793
0.313317195
0.312515942
0.311717026
0.31092044
0.310126176
0.309334226
0.308544582
0.307757237
0.306972183
0.306189412
0.305408916
0.304630689
0.303854721
0.303081007
0.302309538
0.301540307
0.300773307
0.300008529
0.299245967
0.298485613
0.297727461
0.296971501
0.296217728
0.295466134
0.294716712
0.293969455
0.293224355
0.292481405
0.291740598
0.291001928
0.290265386
0.289530966
0.288798662
0.288068465
0.287340369
0.286614368
0.285890453
0.285168619
0.284448859
0.283731165
0.28301553
0.282301949
0.281590414
0.280880919
0.280173456
0.27946802
0.278764603
0.278063199
0.277363802
0.276666404
0.275970999
0.275277581
0.274586143
0.273896679
0.273209183
0.272523647
0.271840065
0.271158431
0.270478739
0.269800983
0.269125155
0.268451251
0.267779263
0.267109185
0.266441011
0.265774735
0.26511035
0.264447851
0.263787232
0.263128486
0.262471607
0.26181659
0.261163427
0.260512114
0.259862644
0.259215011
0.25856921
0.257925233
0.257283077
0.256642734
0.256004198
0.255367465
0.254732527
0.25409938
0.253468018
0.252838434
0.252210623
0.25158458
0.250960298
0.250337772
0.249716997
0.249097967
0.248480675
0.247865117
0.247251288
0.24663918
0.24602879
0.245420111
0.244813138
0.244207866
0.243604289
0.243002401
0.242402198
0.241803674
0.241206823
0.24061164
0.240018121
0.239426259
0.238836049
0.238247486
0.237660565
0.237075281
0.236491628
0.235909602
0.235329196
0.234750406
0.234173227
0.233597654
0.233023681
0.232451304
0.231880517
0.231311316
0.230743695
0.230177649
0.229613174
0.229050264
0.228488915
0.227929121
0.227370878
0.226814181
0.226259024
0.225705404
0.225153315
0.224602752
0.224053711
0.223506187
0.222960174
0.222415669
0.221872667
0.221331162
0.22079115
0.220252627
0.219715588
0.219180027
0.218645942
0.218113326
0.217582175
0.217052485
0.216524251
0.21578875
0.21505607
0.214326198
0.213599121
0.212874829
0.212153308
0.211434547
0.210718534
0.210005256
0.209294702
0.20858686
0.207881719
0.207179265
0.206479489
0.205782378
0.205087921
0.204396106
0.203706922
0.203020358
0.202336402
0.201655043
0.200976269
0.200300071
0.199626436
0.198955353
0.198286812
0.197620802
0.196957311
0.19629633
0.195637846
0.194981851
0.194328332
0.193677279
0.193028682
0.19238253
0.191738814
0.191097521
0.190458642
0.189822167
0.189188086
0.188556388
0.187927063
0.187300101
0.186675492
0.186053226
0.185433293
0.184815683
0.184200386
0.183587393
0.182976693
0.182368278
0.181762136
0.18115826
0.180556638
0.179957262
0.179360122
0.178765209
0.178172513
0.177582025
0.176993736
0.176407636
0.175823716
0.175241967
0.174662379
0.174084945
0.173509654
0.172936498
0.172365468
0.171796555
0.171229749
0.170665043
0.170102427
0.169541893
0.168983431
0.168427034
0.167872693
0.167320398
0.166770142
0.166221916
0.165675712
0.165131521
0.164589334
0.164049145
0.163510943
0.162974721
0.16244047
0.161908184
0.161377852
0.160849468
0.160323023
0.159798509
0.159275918
0.158755242
0.158236474
0.157719605
0.157204628
0.156691534
0.156180316
0.155670967
0.155163478
0.154657842
0.154154051
0.153652098
0.153151975
0.152653674
0.152157189
0.151662511
0.151169633
0.150678549
0.15018925
0.149701729
0.149215979
0.148731992
0.148249762
0.147769282
0.147290543
0.14681354
0.146338264
0.145864709
0.145392869
0.144922735
0.144454301
0.14398756
0.143522505
0.14305913
0.142597427
0.14213739
0.141679011
0.141222285
0.140767204
0.140313763
0.139861953
0.139411769
0.138963204
0.138516251
0.138070904
0.137627157
0.137185003
0.136744435
0.136305447
0.135868034
0.135432187
0.134997902
0.134565172
0.13413399
0.133704351
0.133276248
0.132849675
0.132424626
0.132001094
0.131579074
0.13115856
0.130739546
0.130322025
0.129905991
0.129491439
0.129078363
0.128666757
0.128256615
0.127847931
0.127440699
0.127034914
0.126630569
0.126227659
0.125826179
0.125426122
0.125027484
0.124630257
0.124234437
0.123840019
0.123446996
0.123055363
0.122665114
0.122276244
0.121888748
0.12150262
0.121117855
0.120734447
0.120352391
0.119971681
0.119592313
0.119214281
0.118837579
0.118462203
0.118088147
0.117715407
0.117343976
0.116973849
0.116605023
0.11623749
0.115871247
0.115506289
0.115142609
0.114780204
0.114419068
0.114059197
0.113700584
0.113343227
0.112987118
0.112632254
0.11227863
0.111926241
0.111575081
0.111225147
0.110876434
0.110528936
0.110182649
0.109837568
0.109493689
0.109151006
0.108809516
0.108469214
0.108130094
0.107792153
0.107455385
0.107119787
0.106785353
0.106452079
0.106119961
0.105788994
0.105459174
0.105130497
0.104802957
0.10447655
0.104151273
0.10382712
0.103504088
0.103182171
0.102861367
0.102541669
0.102223075
0.10190558
0.101302325
0.100703034
0.100107676
0.099516223
0.098928646
0.098344918
0.097765009
0.097188893
0.096616541
0.096047926
0.095483021
0.094921798
0.094364232
0.093810295
0.093259962
0.092713205
0.09217
0.09163032
0.09109414
0.090561434
0.090032179
0.089506348
0.088983917
0.088464861
0.087949157
0.08743678
0.086927707
0.086421913
0.085919375
0.08542007
0.084923974
0.084431065
0.083941321
0.083454717
0.082971233
0.082490845
0.082013532
0.081539271
0.081068042
0.080599822
0.080134591
0.079672326
0.079213008
0.078756615
0.078303126
0.077852521
0.07740478
0.076959882
0.076517807
0.076078536
0.075642048
0.075208324
0.074777346
0.074349092
0.073923545
0.073500685
0.073080494
0.072662953
0.072248042
0.071835745
0.071426042
0.071018916
0.070614349
0.070212322
0.069812818
0.06941582
0.06902131
0.068629271
0.068239686
0.067852537
0.067467808
0.067085482
0.066705542
0.066327973
0.065952757
0.065579878
0.065209321
0.064841069
0.064475106
0.064111416
0.063749985
0.063390796
0.063033835
0.062679085
0.062326531
0.06197616
0.061627955
0.061281901
0.060937985
0.060596192
0.060256507
0.059918915
0.059583403
0.059249956
0.05891856
0.058589202
0.058261867
0.057936542
0.057613213
0.057291866
0.056972489
0.056655067
0.056339588
0.056026039
0.055714406
0.055404677
0.055096838
0.054790877
0.054486782
0.054184539
0.053884137
0.053585562
0.053288804
0.052993848
0.052700685
0.0524093
0.052119683
0.051831822
0.051545705
0.05126132
0.050978656
0.050697701
0.050418444
0.050140874
0.049864979
0.049590748
0.04931817
0.049047235
0.04877793
0.048510247
0.048244173
0.047979698
0.047716811
0.047455503
0.047195762
0.046937578
0.046680941
0.046425841
0.046172268
0.045920211
0.04566966
0.045420607
0.04517304
0.04492695
0.044682327
0.044439163
0.044197446
0.043957168
0.04371832
0.043480891
0.043244874
0.043010257
0.042777033
0.042545192
0.042314726
0.042085624
0.041857879
0.041631482
0.041406423
0.041182695
0.040960288
0.040739193
0.040519404
0.04030091
0.040083704
0.039867776
0.03965312
0.039439727
0.039227588
0.039016696
0.038807042
;;;
run;

run;
data pd;
   set p;
   pd = pd / 100;
run;

%macro structgen(data = , var = , num = , out = );
   proc sql noprint;
      create table &out(y num, seed num, grade num);
   quit;
   %do i = 1 %to #
      data _tmp1;
         set &data;
         retain y 0;
         y = y + ranuni(&i)*(&var-y);
         seed = &i;
         grade = _n_ ;
      run;
      proc sql noprint;
         insert into &out
         select y, seed, grade 
         from _tmp1;
      quit;
   %end;
%mend;
%structgen(data = maxbound, var = maxvalue, num = 1000, out = final);

proc sgplot data = final;
   series x = grade y = y / group = seed;
   yaxis label = 'Percentage'; xaxis grid;
run;

/* Part 1 - THE FUNCTION OF CAPITAL REQUIREMENT UNDER BASEL III*/
proc fcmp outlib = work.myfunclib.finance;
   function reqcap(pd, lgd, m);
      corr = 0.12*(1-exp(-50*pd))/(1-exp(-50)) + 0.24*(1-(1-exp(-50*pd))/(1-exp(-50)));
      mtradj = (0.11852 - 0.05478 * log(pd))**2;
      return( (lgd * probnorm((probit(pd) + corr**0.5 * probit(0.999)) 
            / (1-corr)** 0.5) - pd*lgd) * (1 + (m-2.5)*mtradj) / (1-1.5*mtradj) );
   endsub;
quit;

options cmplib = (work.myfunclib);

%macro grdselect(lgd = 0.45, ead = 100, m = 2.5, out = want);
   proc sql noprint;
      select mean(pd)/100 into: avgpd from pd
      ;
      select count(pd) into: nobs from pd
      ;
      create table &out(structure char(50), gini num, cr num)
      ;
   quit;   
   %do i = 1 %to 1000;
   data bound;
      set final;
      lowbound = y ;
      where seed = &i;
      run;

      data _myfmt;
         merge bound(rename=(lowbound=start)) 
         bound(rename=(lowbound=end) firstobs=2);
         fmtname = 'range'; type = 'n';
         label = _n_;
         if end = . then end = 1;
      run; 
      proc format cntlin = _myfmt;
      run;

      data _1;
         set pd;
         grade = put(pd, range.);
         cr = reqcap(pd, &lgd, &m);
      run;


      proc sql noprint;
         select mean(pd) into: avgpd from pd
      ;
      select count(pd) into: nobs from pd
         ;
         select lowbound format=percent8.2 into: structure separated by '-' from bound
         ;
         create table _2 as
         select grade, 
         count(grade) / &nobs as ps 'share of portfolio',
         mean(pd) as gpd 'grade pd',
         calculated ps * calculated gpd / &avgpd as esod 'expected share of default',
         reqcap(calculated gpd, &lgd, &m) as gcr 'grade cr'
         from _1 
         group by grade
       ;
   quit;


   proc sql noprint;
      select avg(gcr) into: cr
      from (   select a.*, b.gcr
           from _1 as a left join _3 as b on a.grade = b.grade)
   ;quit;

   proc sort data = _2;
      by descending grade;
   run;
   data _3;
      set _2;
      retain _y 0;
      _y + esod;
      _x = _y - esod;
      label egc = 'Expected GINI values';
      egc = ps * (_x + esod/2);
   run;

   proc sql noprint;
      select (sum(egc)-0.5) / ((&avgpd/2+(1-&avgpd))-0.5) into:gini 
   from _3
   ;
        insert into &out 
        values ("&structure", &gini, &cr)
        ;
     quit;
   %end;
%mend;
%grdselect;


data want1;
   set want;
   seed = _n_;
run;

proc sgplot data = want1;
   series x = seed y = cr;
   series x = seed y = gini / y2axis;
run;

Conclusion:
1. The optimal rating structure depends on multiple factors, including accuracy ratio and required capital which I used in this test. Business decision would make significant difference, while an integrated evaluation index may need more modeling work.
2. To save time from the high I/O operations involved with lots of loops, such as simulation, backtesting, stress testing, etc., try to decrease the steps to write data to the hard disk. Memory-intensive PROC SQL or PROC IML sometimes is a more efficient way than the DATA step. Also avoid %IF-%THEN statement inside of a macro loop (some interesting discussions similar to this topic: Rick's post here and my post here).

Mahalanobis distances on a heat map

I just learned Mahalanobis distance from Rick’s blog post yesterday, and realized its significance in detecting outliers. One of SAS’s online documents shows how to use PCA method to find Mahalanobis distances. And in SAS 9.3, the popular heat map becomes available

SAS’s classic help dataset SASHELP.CLASS has weight, height, age and some other information for 19 teenagers. I calculated the pair-wise Mahalanobis distances according to their age, weight and height, and showed those distances on a heat map. It seems that it is helpful to tell how similar two teenagers are to each other.

/* 1 -- Find pairwise Mahalanobis distances */
proc princomp data=sashelp.class std out=_1 noprint;
     var age weight height;
run;

proc distance data=_1 out=_2;
   var interval(prin:);
   id name;
run;

/* 2 -- Restructrue data */
data _3(where=(missing(distance)=0));
   set _2;
   array a[*] _numeric_;
   do i = 1 to dim(a);
      x = name;
      y = vlabel(a[i]);
      distance = a[i];
      output;
   end;
   keep x y distance;
run;

data _4;
   set _3 _3(rename=(x=y y=x));
run;

/* 3 -- Draw Mahalanobis distances on a heat map */
proc template;
  define statgraph heatmapparm;
    begingraph;
      layout overlay / xaxisopts=(label=" ") yaxisopts=(label=" ");
        heatmapparm x = x y = y colorresponse = distance / name = "heatmap";
        continuouslegend "heatmap" / orient = vertical location = outside;
      endlayout;
    endgraph;
  end;
run;

ods html style = money;
proc sgrender data=_4 template=heatmapparm;
run;

Valentine's Day

Happy Valentine's Day!

data one;
   do t = 1 to 3*constant("pi") by 0.05;
      x = 16*sin(t)**3;
      y = 13*cos(t) - 5*cos(2*t) - 2*cos(3*t) - cos(4*t);
      output;
   end;
run;

data two;   
   set one;
   if _n_ = 70 then label = "Valentine's Day"; 
run;

ods graphics on / width=6in height= 5in;
proc sgplot data = two;
   series x = x y = y /lineattrs=(color=red thickness=5) 
      datalabel = label datalabelattrs=(color=red 
      family="garamond" style=italic size=45 weight= bold);
run;