The SIMSYSTEM Procedure

The UNIVARIATE procedure in Base SAS software can fit the parameters of all three distribution families in the Johnson system. This example demonstrates how to use PROC SIMSYSTEM in conjunction with PROC UNIVARIATE to study these fits.

Suppose you want to work with data from a Johnson SB distribution that has true moments as given in Table 8.

The following statements use PROC SIMSYSTEM to generate 10,000 observations of such data and to store them in a data table named SB:

proc simsystem system=Johnson n=10000 seed=12345;
   moments mean=5 stddev=2 skew=1 kurt=4.5;
   output out=mycas.SB;
run;

The SIMSYSTEM procedure displays the parameters of this distribution in the table shown in Figure 15 and the observed moments of the 10,000 sampled values in the table shown in Figure 16.

The following statements use PROC UNIVARIATE with the SB option in the HISTOGRAM statement to fit the parameters of a Johnson SB distribution to these data. The default method of fitting an SB distribution uses the percentiles of the data, but in this case, because the moments of the distribution are of particular interest, you use the method of moments to fit the distribution.

proc univariate data=mycas.SB;
   var variate;
   histogram / sb(theta=est sigma=est fitmethod=moments) noplot;
run;

The part of the PROC UNIVARIATE output that shows the fitted parameters for the SB distribution is shown in Figure 17.

As expected from a fit that uses the method of moments, the values of the moments in Figure 17 match those for the sample in Figure 16. The caveat is that PROC SIMSYSTEM works with the direct or raw kurtosis, whereas PROC UNIVARIATE prints the excess kurtosis, which is the raw value minus 3. The values of the four fitted parameters of the distribution match up reasonably well with the true values that they correspond to in the PROC SIMSYSTEM output, with PROC UNIVARIATE’s Theta corresponding to PROC SIMSYSTEM’s Shift and Sigma corresponding to Scale. Delta and Gamma denote the same parameters in both procedures.

It is not at all surprising that, although the sample moments in Figure 16 are fairly close to the specified moments, they do not match very precisely, and the sample skewness and kurtosis are a good deal less precise than the mean and standard deviation. You can explore how far off they are likely to be by using the following statements to simulate 10,000 replicates of this 10,000-observation sample, by using the NREP=10000 option in the PROC SIMSYSTEM code discussed earlier. What you want to look at in this case is a summary analysis of the 10,000 skewness and kurtosis values for each sample, so you specify the MOMENTREPS option to include all the sample moments in the "Distribution Moments" table, instead of just the default averages, and you also use the ODS SELECT NONE statement to suppress printing this table, instead using the ODS OUTPUT statement to direct it to a data set named Moments. Finally, PROC KDE is used to create a contour plot of the distribution of skewness and kurtosis values.

ods select none;
proc simsystem system=johnson n=10000 seed=12345 nrep=10000 momentreps;
   moments mean=5 stddev=2 skew=1 kurt=4.5;
   ods output Moments=Moments;
run;
ods select all;

proc kde data=Moments;
   bivar SampleSkewness SampleKurtosis / plots=contour;
run;

Output 22.1.1: Observed Distribution of Sample Skewness and Kurtosis

The resulting density is shown in Output 22.1.1. Perhaps the most notable feature is that, even though the replicated samples are quite large, containing 10,000 observations each, there is considerable variability in the sample skewness and kurtosis values. This explains the magnitude of the default value for the SIMPLETOL= option. When this option comes into play, PROC SIMSYSTEM does not return a sample that has precisely the expected skewness and kurtosis that you specify, but the difference is undetectable in samples of any reasonable size.