The SIMSYSTEM Procedure

This example demonstrates how to use the table of parameters that PROC SIMSYSTEM displays to reconstruct in DATA step code random variates that have the specified moments for each distribution. Starting with the Pearson system, the moment combinations that are specified for PROC SIMSYSTEM in the following code are chosen to span the families of distributions in the Pearson family, plus the normal distribution, as the parameters table (Figure 22) that PROC SIMSYSTEM produces for these moment combinations confirms:

%let Mean   = -4.5 -3.5 -2.5 -1.5 -0.5  0.5        1.5  2.5        3.5  4.5;
%let Stddev =  1    2    3    4    5    6          7    8          9   10;
%let Skew   =  0   -1.5  1.8 -1.3  0   -2.00831604 0.1 -1.04235831 1.7  0;
%let Kurt   =  3    5.7  7.4  3    2.1  9.05       6    5.15       8.6  3.2;
proc simsystem system=pearson;
   moments mean=&Mean StdDev=&Stddev skew=&Skew kurt=&Kurt;
   ods output Parameters=Parm;
run;

Output 22.5.1: Selected Pearson Distributions

The ODS OUTPUT statement in the preceding PROC SIMSYSTEM code saves the parameters for this sequence of distributions in a data set named Parm. The following PROC PRINT code produces the selected contents of this data set shown in Figure 23. As you can see, the exact numeric values of the parameters that describe each distribution in the Notes column are available in the Shap1Val and Shap2Val variables.

proc print data=Parm noobs;
   var Family Notes Shap: Shift Scale;
run;

As discussed in Parameter Values for Pearson Distributions, you can use the values of the Shap1Val and Shap2Val variables in DATA step code with the formulas given in Table 3 to compute random samples that have the specified moments, as shown in the following DATA step code:

%let N = 1000000;
data mycas.Parm; set Parm; run;
proc ds2 sessref=mysess;
data Sim / overwrite=yes;
   dcl int i;
   dcl double y;
   method run();
      set Parm;
      streaminit(15531);
      do i=1 to &N;
         select (trim(left(Family)));
            when ('Type I'   )
            when ('Type I(J)')
            when ('Type I(U)') y =   rand('beta'  ,Shap1Val,Shap2Val);
            when ('Type II'  ) y =   rand('beta'  ,Shap1Val,Shap1Val);
            when ('Type III' ) y =   rand('gamma' ,Shap1Val);
            when ('Type IV'  ) y = .;
            when ('Type V'   ) y = 1/rand('gamma' ,Shap1Val);
            when ('Type VI'  ) y =   rand('f'     ,Shap1Val,Shap2Val);
            when ('Type VII' ) y =   rand('t'     ,Shap1Val);
            when ('Normal'   ) y =   rand('normal');
         end;
         y = Scale*y + Shift;
         output;
      end;
   end;
enddata;
run;
quit;

The following code computes and displays (in Figure 24) the first four moments for each distribution sampled by the preceding PROC DS2 code. These moments match the specified ones reasonably well.

proc mdsummary data=mycas.Sim;
   var y;
   groupby Family Notes Mean StdDev Skewness Kurtosis;
   output out=mycas.Sim_summary;
run;

proc cas;
   session mysess;
   table.alterTable /
      name="Sim_summary"
      columns={
         {name="_Mean_"       rename="sMean"},
         {name="_Std_"        rename="sStdDev"},
         {name="_Skewness_"   rename="sSkewness"},
         {name="_Kurtosis_"   rename="sKurtosis"}
      };
   table.update /
      table="Sim_summary"
      set={{var="sKurtosis" value="sKurtosis + 3"}};
quit;

proc sort data=mycas.Sim_summary out=Sim_summary; by Family; run;
proc print data=Sim_summary noobs;
   var Family Mean sMean StdDev sStdDev Skewness sSkewness Kurtosis sKurtosis;
   format _NUMERIC_ 6.3;
run;

Finally, for the Johnson system, the following code performs steps precisely like those shown earlier for the Pearson system, except that the moment combinations are chosen to represent the families of the Johnson system, as shown in Output 22.5.2 and Figure 25:

%let Mean   =  -2   -1    0            1  2  ;
%let Stddev =   1    2    3            4  5  ;
%let Skew   =  -1.6  2.5 -1.2922847983 0 -0.7;
%let Kurt   =   7.5  8.4  6.1093660061 3  7  ;
proc simsystem system=johnson;
   moments mean=&Mean StdDev=&Stddev skew=&Skew kurt=&Kurt;
   ods output Parameters=Parm;
run;

Output 22.5.2: Selected Johnson Distributions

Turning to the Johnson system, as discussed in Parameter Values for Johnson Distributions, you can use the values of the Gamma and Delta columns of the "Parameters for Johnson Distributions" table in DATA step code with the formulas given in Table 5 to compute random samples that have the specified moments, as shown in the following code:

%let N = 1000000;
data mycas.Parm; set Parm; run;
proc ds2 sessref=mysess;
data Sim / overwrite=yes;
   dcl int i;
   dcl double y z;
   method run();
      set Parm;
      streaminit(15531);
      do i = 1 to &N;
         z = rand('normal');
         select (trim(left(Family)));
            when ('SB(1)')
            when ('SB(2)') y = 1/(1 + exp(-(z - Gamma)/Delta));
            when ('SL'   ) y =        exp(  z         /Delta) ;
            when ('SU'   ) y =       sinh( (z - Gamma)/Delta) ;
            when ('SN'   ) y =              z                 ;
         end;
         y = Scale*y + Shift;
         output;
      end;
   end;
enddata;
run;
quit;

proc mdsummary data=mycas.Sim;
   var y;
   groupby Family Mean StdDev Skewness Kurtosis;
   output out=mycas.Sim_summary;
run;

proc cas;
   session mysess;
   table.alterTable /
      name="Sim_summary"
      columns={
         {name="_Mean_"       rename="sMean"},
         {name="_Std_"        rename="sStdDev"},
         {name="_Skewness_"   rename="sSkewness"},
         {name="_Kurtosis_"   rename="sKurtosis"}
      };
   table.update /
      table="Sim_summary"
      set={{var="sKurtosis" value="sKurtosis + 3"}};
quit;

proc sort data=mycas.Sim_summary out=Sim_summary; by Family; run;
proc print data=Sim_summary noobs;
   var Family Mean sMean StdDev sStdDev Skewness sSkewness Kurtosis sKurtosis;
   format _NUMERIC_ 6.3;
run;

As was the case with the Pearson system shown earlier, the moments for the samples that the PROC DS2 code computes, shown in Figure 26, match the specified moments reasonably well.