The SIMSYSTEM Procedure

The value at risk () of an investment is the most that you can expect to lose in a given future time frame, with a given probability . One way to think of this is that if you have cash on hand, then you can be % confident that you will be able to cover a loss in that investment. This example uses PROC SIMSYSTEM to explore the accuracy of a common method of estimating the when the assumptions of the method break down.

A popular way to model the periodic returns of an investment is to use a GARCH (generalized autoregressive conditional heteroscedasticity) model, in which, conditional on past returns, future returns are distributed as independent normal values, with a variance that also depends autoregressively on the past. The following code for simulating and fitting a GARCH process is modeled on the first example in SAS Institute Inc. (2019):

data Index;
   call streaminit(12345);
   do t = 1 to 10000;
      e = normal(2);
      output;
      end;
run;

data Index; set Index;
   retain ReturnRate Var eLag varLag;
   if (_N_ = 1) then do;
      eLag   = 1;
      VarLag = 1;
      end;
   Var = 0.4 + 0.8 * eLag**2  + 0.1 * VarLag;
   ReturnRate = sqrt(Var) * e;
   output;
   eLag   = ReturnRate;
   varLag = Var;
run;

proc autoreg data=Index;
   model ReturnRate = / noint garch=(p=1, q=1) dist=t;
   output out=o1 alphacli=%sysevalf(2*0.05) lcl=lcl;
run;

The OUTPUT statement for PROC AUTOREG computes – for this sequence of returns as the 95% lower confidence limit for the forecasted returns, and if you use PROC MEANS as follows to count up the proportion of returns that are less than this (that is, that have greater loss than this), you see that this indeed accounts for about 5% of them, as shown in Figure 20:

data o1; set o1;
   AtRisk = (ReturnRate < lcl);
proc means data=o1 mean;
   var AtRisk;
run;

The PROC AUTOREG code given earlier for estimating – depends on the assumption that the conditional returns are normally distributed, corresponding to the use of the NORMAL() random number generating function in the first DATA step shown previously. How important is this assumption? The following code uses PROC SIMSYSTEM to simulate the random part of 10,000 conditional returns from each of a systematic grid of mostly nonnormal distributions. A subsequent call of the alterTable action renames the output variables for convenience.

proc simsystem n=10000 system=pearson seed=1 plot(only)=mrmap(skewscale=skewness);
   momentgrid skew=-1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1
              kurt=1.25  1.5
                   1.75  1.8
                   1.85  2.0
                   2.1   2.5
                   3     3.5
                   6    10;
   output out=mycas.GARCHSim_error;
run;

proc cas;
   session mysess;
   table.alterTable /
      name="GARCHSim_error"
      columns={
         {name="iObs"    rename="t" label=""},
         {name="Variate" rename="e" label=""}
      };
quit;

These statements specify negative as well as positive skewness values and kurtosis values that run from appreciably platykurtic (kurtosis less than 3, the kurtosis of the normal distribution) to appreciably leptokurtic (kurtosis greater than 3). The MOMENTGRID statement specifies a distribution for every pair of these values, although some are ignored because of the constraint that the kurtosis for a random distribution must always be greater than or equal to . The SKEWSCALE=SKEWNESS option produces a moment-ratio map that shows both negative and positive skewness values. The placement of this grid of skewness and kurtosis values within the moment-ratio map for the Pearson system is shown in Output 22.4.1.

Output 22.4.1: Systematic Grid of Sampled Errors

The following PROC DS2 code uses these simulated values in the same way that the previous DATA step code used random normal values: to produce many samples of returns. A thread block in the PROC DS2 step processes different samples in parallel. The subsequent code performs the PROC AUTOREG analysis discussed earlier for each sample of 10,000 simulated returns, as indexed by the specified skewness and kurtosis values, and similarly uses PROC SUMMARY to count up the proportion whose loss exceeds the computed .

/*using THREAD to parallelize the set statement*/
proc ds2 sessref=mysess;
thread GARCHSim_thd / overwrite=yes;
   dcl double ReturnRate Var eLag VarLag;
   retain ReturnRate Var eLag VarLag;
   method run();
      set GARCHSim_error;
      by SimIndex;
      if FIRST.SimIndex then do;
         eLag   = 1;
         VarLag = 1;
      end;
      Var = 0.4 + 0.8 * eLag**2 + 0.1 * VarLag;
      ReturnRate = sqrt(Var) * e;
      output;
      eLag   = ReturnRate;
      VarLag = Var;
   end;
endthread;
data GARCHSim / overwrite=yes;
   dcl thread GARCHSim_thd GARCHSim_thd_instance;
   method run();
      set from GARCHSim_thd_instance;
   end;
enddata;
run;
quit;

proc autoreg data=mycas.GARCHSim noprint;
   by notsorted Skewness Kurtosis;
   model ReturnRate = / noint garch=(p=1, q=1);
   output out=mycas.g11 alphacli=%sysevalf(2*0.05) lcl=lcl;
run;

data g11; set mycas.g11;
   AtRisk = (ReturnRate < lcl);
proc summary data=g11 mean;
   class Kurtosis Skewness;
   ways 2;
   var AtRisk;
   output out=AtRisk mean=pAtRisk;
run;

Finally, the following code prints these simulated actual coverage values for the nominal , arranging them in a fashion similar to that of the moment-ratio map in Output 22.4.1. The results are shown in Figure 21.

proc sort data=AtRisk;
   by Kurtosis Skewness;
proc transpose data=AtRisk(rename=(Skewness=_NAME_)) out=AtRisk;
   by Kurtosis;
   var pAtRisk;
proc print data=AtRisk(drop=_NAME_) noobs label;
   var N1 N0D75 N0D5 N0D25 _0 _0D25 _0D5 _0D75 _1;
   id Kurtosis;
   format _: N: 6.3;
   label
      N1    = "-1"
      N0D75 = "-0.75"
      N0D5  = "-0.5"
      N0D25 = "-0.25"
      _0    =  "0"
      _0D25 =  "0.25"
      _0D5  =  "0.5"
      _0D75 =  "0.75"
      _1    =  "1";
run;

In Figure 21, you can see that the simulated coverage is reasonably accurate for distributions near skewness 0 and kurtosis 3 of the normal, but the coverage can be severely liberal or conservative for negatively or positively skewed errors, respectively. This fact motivated Koenker and Zhao (1996) to use quantile regression to construct a distribution-free estimate of VaR that is less sensitive to skewness and kurtosis in the underlying error distribution for the returns, as discussed in Rodriguez and Yao (2017).