CSSM Procedure

This example shows how you can use the CSSM procedure to derive the sampling distribution of a complicated function of the forecasts of a response variable by using simulation smoothing. (For more information, see the section Simulation Smoothing.) The example uses the well-known airline passenger series, given as Series G in (Box and Jenkins 1976), and the equally well-known model of it, the Airline model. This series is a monthly series that consists of the number of airline passengers who traveled during the years 1949 to 1960. The data set sashelp.air contains this series (air) along with the date variable (date). The Airline model postulates an ARIMA (autoregressive, integrated, moving average) model for the logarithm of the air series (logAir). That is, the modeling is done in the log scale. However, often the results of the modeling are needed in the original scale. For example, you might want to obtain a model-based confidence interval for the total number (sum of 12 months) of airline passengers in some future year. The following steps summarize how to do this by using simulation smoothing:

The remainder of this example shows the SAS code for carrying out these steps and the results that they produce. The following DATA step creates a data set, mylib.air, from sashelp.air such that it has the log-transformed passenger series (logAir); logAir is set to missing during 1960 so that PROC CSSM produces forecasts and later the simulations for that period:

data mylib.air;
    set sashelp.air;
    if year(date) < 1960 then logAir = log(air);
    else logAir = .;
run;

The following statements fit the Airline model (ARIMA(0,1,1)(0,1,1)12) to the logAir series and save the model in the BLOB table mylib.airModel:

proc cssm data=mylib.air;
    id date interval=month;
    deplag airDiffs(logAir) logAir(lags=(1 12 13) coeff=(1 1 -1));
    trend ma(arma(q=1 sq=1 s=12));
    model logAir = airDiffs ma;
    simsmooth outmodel=mylib.airModel;
run;

The following statements use the fitted model that is saved in the mylib.airModel table and the input data table mylib.air to generate the desired simulations:

proc cssm data=mylib.air;
    simsmooth inmodel=mylib.airModel seed=1 nsim=1000
        outpred=mylib.airPred;
run;

Output 14.21.1 shows the summary of the simulation smoothing process.

Because of the option NSIM=1000, PROC CSSM produces 1,000 draws of predictions of logAir for the 12 months of 1960. The simulated predictions are output to a table, mylib.airPred, for postprocessing. In total, mylib.airPred contains 1,001 sets of values (12 rows for each set that correspond to the months of 1960):

The table mylib.airPred has three columns: SEED contains the seed values, Date contains the observation dates, and SimSmooth_logAir contains the actual predictions or their simulations. The following statements create a data set, yearly, that contains the 1,000 draws of TotalAir1960 = .

data airPred;
    set mylib.airPred;
    where seed > -1;
    simAir = exp(simsmooth_logAir);
run;
proc sort data=airPred;
     by seed date;
run;
proc means data=airPred noprint;
    by seed;
    var simAir;
    output out=yearly sum=TotalAir1960;
run;

The following statements compute the actual total number of passengers in 1960, which turns out to be 5,714:

proc means data=sashelp.air sum;
    where year(date)=1960;
    var air;
run;

By using the 1,000 draws of TotalAir1960 in the yearly table, you can easily obtain a variety of simulation-based measures of sampling distribution of TotalAir1960. For example, Output 14.21.2 shows the mean, median, and a few percentiles of the sampling distribution of TotalAir1960, which is generated by the following statements:

proc summary data=yearly print p1 q1 mean median q3 p99;
    var TotalAir1960;
run;

The table shows that the mean (5876.4) and the median (5862.9) of the sampling distribution are quite close to the actual total (5714) observed in 1960. Moreover, the interval (5133, 6658) can serve as a 98% confidence band for the total number of passengers in 1960. Finally, the following statements generate the graph shown in Figure 25, which displays the sampling distribution of TotalAir1960:

proc sgplot data=yearly;
    title Sampling Distribution of TotalAir1960;
    histogram TotalAir1960;
    density TotalAir1960 / legendlabel="Normal" name="Normal";
    refline 5714.00 / axis=x lineattrs=GraphFit2
        legendlabel="Actual Total in 1960" name="Coeff";
    discreteLegend "Normal" "Coeff" / across=1 location=inside;
run;

Figure 25: Histogram of TotalAir1960