Hierarchical Forecasting with SAS® Visual Forecasting
Overview
This document presents an overview and example of hierarchical forecasting methods using SAS Visual Forecasting procedures and packages. A hierarchical time series is one in which each of the unique series are arranged into a hierarchical structure of nodes using time series attributes such as geographic location, product type, etc. The bottom level contains all the time series with each having a unique combinations of attribute values. Each successively higher level considers one less attribute for delineating the time series and aggregates each set of child nodes from the lower level into a parent node. The terms child and parent are relative, such that any node in any level (except the top level) is a child of a parent in a higher level. Likewise, any node in any level (except the bottom level) is a parent of one or more child nodes. The time series in a parent node is typically aggregated as either the sum or the average of the data from its associated child nodes.
For example, the SASHELP library dataset, pricedata, tracks the product sales delineated by the attributes regionName, productLine, and productName. Figure 1 shows a hierarchy of pricedata where unique values for the regionName attribute are denoted as R1, R2, etc., productLine as L1, L2, etc., and productName as N1, N2, etc.
Figure 1: Hierarchical structure of pricedata data set used in this paper
The motivation for hierarchical forecasting lies in the amount of detail contained in the time series at each level of aggregation. Time series at lower levels of the hierarchy tend to be sparse or noisy, while higher levels may have important details washed out by too much aggregation. Choosing a reconciliation level somewhere in the middle of the hierarchy, where the forecasts are stronger, enables propagation of the forecasts up or down to obtain more reliable forecasts at every level.
Steps for generalized hierarchical forecasting
Define the hierarchical structure and aggregate data for each level
Generate independent statistical forecasts for each node in the hierarchy using a suitable forecasting model for the included series.
Lower level, sparse, and noisy series may work best with ESM (Exponential Smoothing Models)
Higher levels with more aggregated data may work better with ARIMA (AutoRegressive Integrated Moving Average), UCM (Unobserved Component Model), or ESM
Choose the best hierarchy level to be the reconciliation level
Lower levels with little or no aggregation may contain too much noise
Higher levels may wash out details from lower levels as they are aggregated together
Reconcile the statistical forecasts per the hierarchy structure rules. In other words, adjust the forecast values such that the reconciled forecasts add up correctly to match the statistical forecasts at the reconciliation level.
Example:
Define Hierarchy and Generate Statistical Forecasts
Regardless of reconciliation method, the first step in hierarchical forecasting is to aggregate the data into individual time series for each hierarchy node and generate the statistical forecasts. The following code illustrates an example of the three level (plus level zero) hierarchical forecast of the pricedata data set as illustrated in Figure 1.
The first code section sets up a CAS session called mycas, and a libref called mylib, to point to the input and output data. The options statement is required for the table naming convention in this example where each table will be prefixed by “Lx.”, where x is the level number of the hierarchy. For example, the output forecast table for level two will be named L2.outFor. This option is not required for strictly alpha-numeric table names. The pricedata data set is then loaded into CAS under the mylib reference with a DATA step.
cas mycas;
libname mylib cas sessref = mycas;
options VALIDMEMNAME = EXTEND;
data mylib.pricedata;
set sashelp.pricedata;
run;
PROC TSMODEL will be called once for each level of the hierarchy to aggregate the data and generate statistical forecasts using the ATSM (Automatic Time Series Modelling) package. In this particular example, similar code will be used in each call. The following macro, timeDataCode, can be used to avoid rewriting similar code several times, but this is not always the case.
The macro takes one argument, an indicator to only use ESM (Exponential Smoothing Models) when choosing the best forecast method. ESM is more robust and less time consuming than ARIMA and UCM for sparse or noisy time series, so this example specifies to use only ESM on the bottom level, but consider ESM, ARIMA and UCM otherwise. Within the macro, the objects used by the ATSM package are declared first and then each one is configured as follows:
The dataFrame object is initialized and the dependent (sale) and independent (price) variables are added to it.
The diagSpec object is populated with the time series diagnose specifications of each model to be considered.
The diagnose object is initialized with the dataFrame and the diagSpec objects and run to generate candidate models.
The forecast object is then initialized with the diagnose object. The model options are specified and the forecasts are generated.
The collector objects, outFor, outEst, and outStat collect the forecast, model parameter estimates, and model fit statistics, respectively.
*define the script to run in TSMODEL;
%macro timeDataCode (ESM_ONLY = NO);
*declare ATSM objects;
declare object dataFrame(TSDF);
declare object diagnose(DIAGNOSE);
declare object diagSpec(DIAGSPEC);
declare object forecast(FORENG);
declare object outFor(OUTFOR);
declare object outEst(OUTEST);
declare object outStat(OUTSTAT);
*setup dependent and independent variables;
rc = dataFrame.initialize();
rc = dataFrame.addY(sale);
rc = dataFrame.addX(price);
*setup time series diagnose specifications;
rc = diagSpec.open();
rc = diagSpec.setEsm();
%if %upcase(&ESM_ONLY) = NO %then %do;
rc = diagSpec.setArimax();
rc = diagSpec.setUcm();
%end;
rc = diagSpec.close();
*diagnose time series to generate candidate model list;
rc = diagnose.initialize(dataFrame);
rc = diagnose.setSpec(diagSpec);
rc = diagnose.run();
*run model selection and forecast;
rc = forecast.initialize(diagnose);
rc = forecast.setOption('lead', 12, 'holdoutpct', 0.1);
rc = forecast.run();
*collect forecast results;
rc = outFor.collect(forecast);
rc = outEst.collect(forecast);
rc = outStat.collect(forecast);
%mend;
The first call to PROC TSMODEL aggregates the time series and generates statistical forecasts for level zero of the hierarchy. For this call, the output table names are prefixed with “L0.” to indicate the hierarchical level. This naming convention requires enclosing the table name in quotes and appending the letter n to indicate that the entire string should be treated as a name literal. This is also the reason for the options statement discussed previously. In level zero, all the time series are aggregated into one and none of the time series attributes are considered, so no BY statement is required. Note that the sale variable is aggregated as a summation while the price variable is averaged.
/* Aggregate time series data and generate forecasts
at level 0 of hierarchy (no BY variables) */
proc tsmodel data = mylib.pricedata
outSum = mylib.'L0.outSum'n
outobj = (
outFor = mylib.'L0.outFor'n
outEst = mylib.'L0.outEst'n
outStat = mylib.'L0.outStat'n
)
;
id date interval = month;
var sale /acc = sum;
var price/acc = average;
*use ATSM package;
require atsm;
submit;
%timeDataCode(ESM_ONLY = NO);
endsubmit;
run;
The second call to PROC TSMODEL aggregates the time series and generates forecasts for level one of the hierarchy. For this call, the output table names are prefixed with “L1.” Here the time series are delineated by regionName, as specified in the BY statement.
/* Aggregate time series data and generate forecasts
at level 1 of hierarchy (regionName) */
proc tsmodel data = mylib.pricedata
outSum = mylib.'L1.outSum'n
outobj = (
outFor = mylib.'L1.outFor'n
outEst = mylib.'L1.outEst'n
outStat = mylib.'L1.outStat'n
)
;
id date interval = month;
var sale /acc = sum;
var price/acc = average;
by regionName;
*use ATSM package;
require atsm;
submit;
%timeDataCode(ESM_ONLY = NO);
endsubmit;
run;
Then, the third call to PROC TSMODEL aggregates the time series and generates forecasts for level two, where the series are delineated by regionName and productLine attributes. For this call, the output table names are prefixed with “L2.”
/* Aggregate time series data and generate forecasts
at level 2 of hierarchy (regionName and productLine) */
proc tsmodel data = mylib.pricedata
outSum = mylib.'L2.outSum'n
outobj = (
outFor = mylib.'L2.outFor'n
outEst = mylib.'L2.outEst'n
outStat = mylib.'L2.outStat'n
)
;
id date interval = month;
var sale /acc = sum;
var price/acc = average;
by regionName productLine;
*use ATSM package;
require atsm;
submit;
%timeDataCode(ESM_ONLY = NO);
endsubmit;
run;
Finally, PROC TSMODEL is called a fourth time to generate forecasts for level three, where the series are delineated by regionName, productLine, and productName. For this call, the output tables are prefixed with “L3.” Note that the timeDataCode macro argument specifies to consider Exponential Smoothing Models only for this level.
/* Aggregate time series data and generate forecasts at
level 3 of hierarchy (regionName, productLine, and productName) */
proc tsmodel data = mylib.pricedata
outSum = mylib.'L3.outSum'n
outobj = (
outFor = mylib.'L3.outFor'n
outEst = mylib.'L3.outEst'n
outStat = mylib.'L3.outStat'n
)
;
id date interval = month;
var sale /acc = sum;
var price/acc = average;
by regionName productLine productName;
*use ATSM package;
require atsm;
submit;
%timeDataCode(ESM_ONLY = YES);
endsubmit;
run;
Reconcile Statistical Forecasts
A reconciliation level, where the forecasts are believed to be the strongest, is chosen as the starting point for reconciliation. From there, the forecasts may be reconciled down to lower levels (Top Down Reconciliation) or up to higher levels (Bottom Up Reconciliation). In this example, the reconciliation level is chosen as level two, where the time series are delineated by regionName and productLine. From this reconciliation level, both top down and bottom up reconciliation will be demonstrated as illustrated in Figure 2.
Figure 2: Direction of top down and bottom up reconciliation relative to the reconciliation level (level two).
Due to the distributed nature of the data, different coding methods are required for the two reconciliation methods. In top down reconciliation, the lower level time series to be reconciled exist on different workers that do not communicate with each other. For this method, PROC TSRECONCILE manages the communication and reconciliation. In bottom up reconciliation, the starting time series and all attribute values exist on the same worker and can be easily aggregated to a subset of the attributes using PROC TSMODEL. Examples of each method are shown in the following sections.
Top Down Reconciliation
PROC TSRECONCILE is used to reconcile hierarchical forecasts from a higher level to a lower level and requires forecast tables PROC TSMODEL from the parent and child levels as input. In this example, the parent= argument specifies the name of the level two forecast table and the child= argument is the level three forecast. The outfor= argument specifies the name of the output table which contains the complete reconciled forecast for level three after the code is run.
/* Reconcile forecasts from the level 2 of the hierarchy
(regionName and productLine), to the lower levels of the hierarchy (productName) */
proc tsreconcile child = mylib.'L3.outfor'n
parent = mylib.'L2.outfor'n
outfor = mylib.'L3.recfor'n;
by regionName productLine productName;
id date;
run;
Finally, if two or more levels to be reconciled exist below the reconciliation level, then PROC TSRECONCILE is called iteratively to propagate reconciliation down the hierarchy. After the first level is reconciled, each successive iteration will call PROC TSRECONCILE with the previously generated reconciled forecast passed in as the parent and the next level down as the child.
Additional Options in PROC TSRECONCILE
PROC TSRECONCILE also includes additional functionality to specify names for confidence limits, error, and standard error variables in both child and parent tables if they exist. Other options available include specifying aggregation and disaggregation methods, calculation or shifting of confidence limits, handling of zero and negative values, and weighting of loss function to emphasize the effect of forecasts with lower variance during reconciliation. The syntax for PROC TSRECONCILE is shown below.
proc tsreconcile <options>;
by variables;
id variable <options>;
childRoles <options>;
parentRoles <options>;
run;
The following tables summarize the additional options available in PROC TSRECONCILE. Table 1 lists options for the procedure statement itself, while Table 2 lists options for specifying variable names in the childRoles and parentRoles statements.
Table 1: Summary of options for the PROC TSRECONCILE statement
Table 2: Summary of options for childRoles and parentRoles statements in PROC TSRECONCILE
Bottom Up Reconciliation
Reconciliation from the bottom up requires a little extra care to ensure the data is complete. If any lagged or differenced values were considered in the model, then the fitted model will have missing values at the beginning of each time series where the historical data to calculate the lags and/or difs was not available before the beginning of the time series. During aggregation, the missing values will be ignored resulting in artificially low aggregated values. To prevent this, any missing values are replaced with the actual historical values before aggregation begins. Then, the aggregated forecast can be used to reconcile the statistical forecast of the higher level. Additionally, the confidence interval will need to be adjusted to center around the new reconciled forecast. This process is demonstrated in the following code sections.
Replace missing values and aggregate data to next hierarchy level
Instead of using PROC TSMODEL, the following code demonstrates how to use the CAS actions directly to perform bottom up reconciliation. More specifically, the timeData action will be used within PROC CAS to create a new table, replace missing values, and aggregate the forecasts.
It should be noted that the input table, L2.outFor, is formatted such that the name of the variable from the input data set is written to one variable, _NAME_, while the values for the actual historical data and the predicted values are written to new variables named ACTUAL and PREDICT, respectively. During the following reconciliation steps, these names will be used to refer to the time series variables.
The cmpCode statement defines the script to be run by the timeData action, then the input and output specifications are defined within the timeData action. The table statement specifies the input table, L2.outFor, and to group the values by the list of attributes for the next highest level in the hierarchy. In this case, the next level up uses only one attribute, regionName, to delineate the time series. The replacement of missing values is easily handled using the computedVars and computedVarsProgram options in the table statement. This creates a new variable, new_predict, and populates it with the forecast model values or the actual values if the model value is missing.
The arrayOut statement specifies the output table name, L2._aggTable, and the replace existing table option. Then, the series statement gives the variables to include in the output and specifies they are to be aggregated as summations. The code statement runs the script defined at the beginning of the procedure which copies the aggregated variables from new_predict into the predict variable
proc cas;
cmpCode = "do i = 1 to dim(predict);
predict[i] = new_predict[i];
end;
";
timeData.runTimeCode /
table={name='L2.outFor',
groupBy={'regionName'},
computedVars = {'new_predict'},
computedVarsProgram = 'if predict = . then new_predict = actual;
else new_predict = predict;'}
timeId= 'date'
interval='month'
arrayOut={
table={name='L2._aggTable', replace=true}
}
series={
{name='actual', acc='sum'},
{name='predict', acc='sum'},
{name='new_predict', acc='sum'}
}
code= cmpCode;
run;
quit;
Merge parent forecast with aggregate forecast
Now the statistical forecasts in the upper level must be merged with the reconciled forecasts aggregated from the reconciliation level. Proc fedSQL performs this task and renames the predict variable from the lower level to _agg_predict. The resulting table, _tmp_merge, contains all the data from the upper level time series and the reconciled forecast from the lower level. The predict variable is the statistical forecast and _agg_predict is the reconciled forecast.
proc fedSQL sessref=mycas;
create table _tmp_merge {options replace=true} as
select high.*, low.predict as _agg_predict
from "L1.outFor" as high, "L2._aggtable" as low
where high.regionName=low.regionName and high.date=low.date;
quit;
Adjust confidence interval limits
The confidence interval limits are still centered around the statistical forecast and must be shifted to center around the reconciled forecast. In the data step below, the statistical forecast is renamed to _orig_predict and the aggregate forecast is _agg_predict. The difference is added to the values for the confidence interval variables, lower and upper.
data mylib.'L1.recfor'n / sessref=mycas;
set mylib._tmp_merge(rename=(predict=_orig_predict) );
if missing(_agg_predict) then predict = _orig_predict;
else predict = _agg_predict;
*recenter upper/lower limits;
lower = lower + (predict - _orig_predict);
upper = upper + (predict - _orig_predict);
keep regionName date actual predict std lower upper;
run;
Now the bottom up reconciliation is complete and the data table, recfor, contains the reconciled forecasts at the regionName level based on the aggregated forecasts from the reconciliation level. This sequence of replacing missing values, aggregating the lower level forecasts, reconciling to the next higher level, and adjusting the confidence interval limits may be repeated iteratively in a macro for hierarchies with additional levels.
Model Fit Statistics
The model fit statistics for the final reconciled forecast can be obtained with the UTL package. The final code section shows an example of passing the reconciled forecast, recfor, to PROC TSMODEL and using the UTL package to compute and collect the fit statistics in the table recstat. To obtain the fit statistics for different levels, one need only change the by variables and the file name for the input data.
*Generate model fit statistics;
proc tsmodel data=mylib.'L1.recfor'n
outobj = (utlStat = mylib.'L1.recstat'n);
require utl;
id date interval = month;
by regionName;
var actual predict;
submit;
declare object utlStat(UTLSTAT);
rc = utlStat.collect(actual, predict);
endsubmit;
run;
Summary
This document has demonstrated how to generate hierarchical forecasts by reconciling statistical forecasts of multiple time series arranged by time series attributes. PROC TSMODEL was used to generate statistical forecasts for each level of the hierarchical structure and level two was chosen as the reconciliation level. Top down reconciliation was demonstrated using PROC TSRECONCILE and bottom up reconciliation was shown using PROC TSMODEL. Finally, the reconciled forecast model fit statistics were obtained using PROC TSRECONCILE and the UTL package. For more information, see the documentation for the ATSM and UTL packages, and PROC TSRECONCILE in the SAS Viya forecasting documentation.
... View more
