{'fraction_non_alphanumeric': 0.03354954148959964, 'fraction_numerical': 0.012935212107656751, 'mean_word_length': 5.005596597268861, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Multi-Phase Approach for Product Hierarchy Forecasting in Supply Chain Management: Application to MonarchFx Inc

Hierarchical time series demands exist in many industries and are often associated with the product, time frame, or geographic aggregations. Traditionally, these hierarchies have been forecasted using top-down, bottom-up, or middle-out approaches. The question we aim to answer is how to utilize child-level forecasts to improve parent-level forecasts in a hierarchical supply chain. Improved forecasts can be used to considerably reduce logistics costs, especially in e-commerce. We propose a novel multi-phase hierarchical (MPH) approach. Our method involves forecasting each series in the hierarchy independently using machine learning models, then combining all forecasts to allow a second phase model estimation at the parent level. Sales data from MonarchFx Inc. (a logistics solutions provider) is used to evaluate our approach and compare it to bottom-up and top-down methods. Our results demonstrate an 82-90% improvement in forecast accuracy using the proposed approach. Using the proposed method, supply chain planners can derive more accurate forecasting models to exploit the benefit of multivariate data.


Introduction
The efficient movement of goods in a supply chain depends on the ability to accurately forecast product demands. Oftentimes, these forecasts must be produced within a hierarchical structure which may represent geographic regions, product families, (Hyndman et al., 2011), or time periods (Athanasopoulos et al., 2017). The value of hierarchical forecasting is that it can provide decision support information to different stakeholders across various organizational functions and managerial levels (Fliedner and Mabert, 1992). For instance, hierarchical forecasts can be used to improve market positioning, inventory planning, facility layouts, or increased efficiency of operational logistics and transportation networks, leading to increased customer satisfaction and lower costs. Muir (1979) explained how hierarchical forecasting can increase overall forecast accuracy, noting that combining data from two or more homogeneous items can produce a stabilizing effect.
Two dominant approaches exist in the hierarchical forecasting literature: top-down and bottom-up. In the top-down approach, a forecast is initially created at an aggregated level, then disaggregated to lower levels of the hierarchy (Boylan, 2010). A common disaggregation approach involves proration (Fliedner, 1999;Strijbosch et al., 2007) in which the aggregate demand forecast is multiplied by the ratio of corresponding demand to aggregate demand, resulting in an estimate for the next lower-level in the hierarchy. In the bottom-up approach, the steps are reversed. The lowest level of the hierarchy is forecasted first (i.e. SKU level), then aggregated to estimate higher levels in the hierarchy (Hyndman et al., 2011). A third approach called middle-out combines aspects of top-down and bottom-up. In middle-out, the forecast is performed at a middle level of the hierarchy, then aggregated up and disaggregated down to estimate the forecasts for other levels of the hierarchy.
With respect to top down forecasting, Gross and Sohl (1990) argued that two simple disaggregation techniques can be effective; "average historical proportions" and "proportions of the historical averages" (Athanasopoulos et al., 2 2009). In the "average historical proportions" approach, the share of each lower level time series of the aggregated series is calculated across all periods, i.e. a linear average share is used. In the"proportions of the historical averages" approach, a volume weighted share across all time periods is employed. The authors also mention that for the "average historical proportions" approach, one is not required to only use the historical proportions, but can utilize the forecasted proportions instead. Promising results were derived using this approach and it is offered in some forecasting software (Boylan, 2010).
In practice, there may be multiple features in the input data (e.g. date, time, holidays, seasonal discounts, etc.) that can be leveraged to improve forecast accuracy within supply chains. To the best of our knowledge, most of the research in the supply chain hierarchical forecasting literature is univariate.
We found no documented multivariate hierarchical forecasting models that employ lower level forecasts as features in parent level modeling. In this research, we employ multiple features (in contrast to univariate time series data) and child level (SKUs) and parent level (brand) forecasts in a hierarchical supply chain model to improve forecast accuracy at the parent level in the hierarchy.
We utilize Machine Learning (ML) techniques including Multi-Layer Perceptron (MLP), Random Forest (RF), Gradient Boosting (GB), and Extreme Gradient Boosting (XGB) to build competing forecasting models. The rest of the paper is organized as follows: In section 2, we briefly review the various existing hierarchical forecasting methods and the aggregation approaches in use. In section 3, we present the details of our proposed Multi-Phase Hierarchical forecasting approach (MPH). We then describe numerical experiments that demonstrate the performance of MPH using sales data from MonarchFx Inc. (a logistics solutions provider) that is representative of a mid-tier supply chain customer in section 4. We summarize our conclusions and discuss practical aspects of our work in section 5. 3

Literature Review
The performance of top-down and bottom-up forecasting approaches in the literature are mixed (Syntetos et al., 2016). Some authors found top-down approaches to be superior (Barnea and Lakonishok, 1980;Gross and Sohl, 1990;Fliedner, 1999), while others found bottom-up methods to be more accurate. (Dangerfield and Morris, 1992;Gordon et al., 1997). These conflicting results occur because the performance of each approach depends on the nature of demand for the products involved. To illustrate, consider a three-level product hierarchy, with product sales at the lowest level, group sales at the middle level, and category sales at the top level. Since group sales are determined by the sum of product sales (given the additive nature of the hierarchy), and the sum of group sales determines category sales, the underlying demand process is transformed at different levels of the hierarchy. When aggregating, a significant loss of information can occur, which tends to render bottom-up forecasting more favorable. Conversely, in the top-down approach, benefits can occur due to random noise cancellation (Fliedner, 1999). Because the performance of each approach depends on the demand generation process within the data, a wide range of conflicting results appears in the literature. Thus, depending on the demand process and parameter settings, one approach may perform better than the other in different contexts (Widiarta et al., 2007(Widiarta et al., , 2009).
An early study comparing both top-down and bottom-up approaches was conducted by Grunfeld and Griliches (1960), in which they found the top-down approach more accurate, with the explanation that disaggregated data is more susceptible to error. Fogarty and Hoffmann (1983) and Narasimhan et al. (1995) derived similar conclusions in their work. Conversely, the loss of information in a top-down approach was considered substantial in Orcutt et al. (1968) and leading to the conclusion that the bottom-up approach is superior. In Shlifer and Wolff (1979), the authors identified conditions on the hierarchy's structure and forecast horizon, under which they concluded that the bottom-up approach is favorable. The robustness and bias of both approaches were investigated in Schwarzkopf et al. (1988). The authors concluded that the bottom-up approach is more favorable unless there exist unreliable or missing data at the bottom of the hierarchy.
A significant characteristic of the underlying demand process involves the dependencies between demand produced at each level, which can be a reason for the performance differences between top-down and bottom-up approaches (Chen and Boylan, 2007). Sbrana and Silvestrini (2013) summarizes the arguments that are often made against top-down approaches. First, he states that the a high (or low) variance in one level in a hierarchy may be indicative of high (or low) variance at other levels. In such cases, allocating measures of variance from higher levels to lower levels in a hierarchy may yield better results. Second, since different products may be classified in different segments, the aggregation of data will lead to a loss of information, making the top-down approach less appealing.
On the other hand, there are examples in the supply chain forecasting literature where the authors favor the top-down approach. In Boylan (2010), the author found that aggregated data can lead to more accurate sales forecasts when dealing with change policies (e.g. change in pack sizes), compared to individual level forecasts. In such cases, common disaggregation techniques ("average historical proportions" and "proportions of the historical averages") may not be useful, and judgmental estimates are required in disaggregation methods to handle such changes in policy.
One method to overcome these drawbacks involves analysis of the conditions in which each approach produces superior forecasting accuracy outcomes.
In Widiarta et al. (2008), the top-down and bottom-up approaches are compared in the context of production planning. Their goal was to estimate requirements at the SKU level. The aggregate demand series were assumed to have correlated sub-aggregate components, each of which were assumed to follow a first order univariate moving average (stationary) process correlated over time. They con- They concluded that when all sub-aggregate components of the time series follow a first-order univariate moving average process with identical coefficients of the serial correlation term, the relative performance of both top-down and bottom-up approaches are similar. Additionally, the different coefficients of the serial correlation term among sub-aggregate components were examined in a simulation study. The result was that the differences in the performance are relatively insignificant when there are small or moderate correlations between the sub-aggregate components. Sbrana and Silvestrini (2013) found that when moving average parameters are not identical, the performance of top-down and bottom-up approaches is similar.
More recently, Rostami-Tabar et al. (2015) analyzed theoretically and by means of simulation (using theoretically generated data) the relative performance of top-down and bottom-up forecasting methods for both aggregate and SKU level demand. The latter was assumed to follow a non-stationary ARIMA (0,1,1) demand process and exponential smoothing (which is optimal for this demand process). An important finding was that the forecast accuracy improvements achieved by bottom-up and top-down methods for non-stationary demands are higher than those associated with stationary cases. The theoretical findings were validated through empirical analysis on data from a European superstore.
A limitation observed in this work is that the generation of forecasts is dominated by the time series at a single level of aggregation (the point at which forecasts are created). To overcome this issue, a regression-based approach was introduced by Hyndman et al. (2011). In their approach, they estimated the time series at multiple hierarchy levels and then optimized this combination using linear regression. This approach sought to derive the benefits of an ensemble of bottom-up and top-down approaches, employing a linear combination of both. Their method demonstrated a significant improvement in forecast accuracy compared to the traditional approaches. This improvement was believed to be a function of employing a combination of forecasts that reduced the variance 6 of forecast error (Timmermann, 2006;Barrow and Kourentzes, 2016). Hyndman et al. (2011) conclude that their proposed combination method is optimal, and compared to all combination forecasts, leads to the least variance. Their work is inspired by earlier research in economics focusing on revising measurements of macro-economic indicators (Zellner and Tobias, 2000;Espasa et al., 2002;Hubrich, 2005). Other research focuses on using different sources to combine forecasts, e.g. utilizing different available information provided by human experts (Budescu and Chen, 2014;Lamberson and Page, 2012). Additionally, the combination of forecasts may reduce model specification and estimation uncertainty (Kourentzes et al., 2014). In a later work, Hyndman et al. (2014) demonstrate the extendibility of their combination approach for hierarchical forecasting to non-hierarchical time series, and time series with partial hierarchical structure. They also proposed a solution to solve the scalability problem that existed in their previous paper Hyndman et al. (2011). They use a linear model structure for a more efficient coefficient estimation.
In Pennings and van Dalen (2017), the authors utilize all the series in a hierarchy in contrast to a top-down or bottom-up approach. They then incorporate a Kalman filter and state space model to comprehend the dependencies between products (e.g. product substitution, product complementarity). Using a multi-variate state space, one is able to estimate the hierarchical time series efficiently using a Kalman filter as a prediction error decomposition tool (Durbin and Koopman, 2012). In this manner, multiple methods for forecasting hierarchical time series exist (Hyndman and Khandakar, 2008;Snyder et al., 2012).
In their approach, forecasts for the aggregate level is derived by summing the forecast of product sales at the base level. The Kalman filter is then used to track the forecast error of individual series at each level of the hierarchy back to the associated states. In this manner, the forecast leverages the information from all series. The authors conclude that their approach is superior to the traditional top-down and bottom-up approach since they incorporate information from all levels of the hierarchy.
Our work builds on the research by Hyndman et al. (2011), and Pennings and van Dalen (2017) (discussed previously), in which they combine information at all levels of hierarchy to improve forecasting accuracy. However, these authors only employ univariate data as their input, and do not leverage multiple features. Our main contribution in this paper is to propose a novel approach which 1) utilizes forecasts at lower levels to improve forecasts at higher levels, 2) uses multivariate data at each level of the hierarchy instead of univariate data, which is more commonly seen in the literature, and 3) leverages machine learning models. The latter component is, to the best of our knowledge, a novel application in the supply chain forecasting literature. To achieve our goal, we propose an MPH approach which is discussed in the following section.

Multi-Phase Hierarchical Forecasting Approach
Our goal is to find a small loss value, l(.), in the parent level of the hierarchy, to optimize: where ω is the matrix of the weights, x i is the vector of the inputs from the i th instance, y i is the dependent variable, e.g., demand (sales) values, and θ(.) is the output function defined by the forecasting model. The well-known Mean Absolute Error (MAE) is being used as our loss function, l(.), in which, the average of differences between the actual demands and estimated demand is calculated.
To achieve a higher level of accuracy in the parent level of the hierarchy, we utilize an MPH approach. In the first phase, we forecast at both child level (SKU level) and parent level (brand) demands using several machine learning approaches. Then, for each individual time series, we select the most promising forecast method in terms of MAE. MAE is calculated based on a cross-validation technique. In the second phase, we aggregate the forecasts at the child level and parent level and use them as an input for the multi-feature forecasting approach.

Overview of Forecasting Methods
Conventional parametric forecasting techniques include ARIMA, GARCH, and TRANSFER models (Box et al., 2015;Shumway and Stoffer, 2011). Moreover, Taylor (2000) forecasts Supply chain forecasting is a field which generally consists of very noisy data, thus it is important to control for noise and learn the true underlying demand patterns which are likely to be repeated in the future. The universal approximators discussed earlier have two desirable features which make them suitable for the supply chain forecasting problem, while dealing with noise. The first is that they are capability of learning any arbitrary function, while the second feature is the capability to control the learning process.
Since we want to exploit additional information provided by multiple input features, we are faced with a multi-dimensional input data vector. The traditional parametric forecasting models such as ARIMA are not able to integrate multi-dimensional inputs, thus we exploit the ability of universal approximators to take multi-dimensional inputs and utilize them in our forecasting model. The details of the MPH approach are explained in the following subsection.

MPH Algorithm
Phase I: • Step 0: Choose forecasting model types which support multi-feature inputs, e.g. MLP, RF, GB, XGB, etc. Suppose we have chosen N forecasting approaches. Set i = 1. • Step 1: Use the i th forecasting approach to forecast parent level demand (brand demand) and child level demands (SKU demand). • Step 2: Optimize the hyperparameters of the i th forecasting method using a search approach, e.g. Bayesian optimization method, grid search, successive halving. • Step 3: Set i = i + 1. If i = N + 1, go to step 4. Otherwise, go to step 1. • Step 4: Using the outputs of the previous steps, record the best forecasting approach and the associated outputs for demands at all levels.
Phase II: • Step 5: Append the recorded outputs of step 4 to input data of the parent level, as additional features. • Step 6: Repeat steps 1 and 2 once more, using the new input data with additional features. The only modification is to only forecast the parent level.

•
Step 7: Choose the best forecasting output among the forecasting methods used in step 6. Hence, one needs to find an approach to optimize the hyperparameters of the forecasting models, which is illustrated below:

Hyperparameter Optimization
There are several approaches in the literature that address hyperparameter optimization in machine learning (Maclaurin et al., 2015;Feurer et al., 2015;Li et al., 2017;Bergstra et al., 2013). In this paper, we use the hyperOpt algorithm proposed by Bergstra et al. (2013), and combine it with the successive halving approach (Jamieson and Talwalkar, 2016) to obtain a more efficient search. In the following a summary of the HyperOpt algorithm is provided and the details of the proposed hyperparameter optimization algorithm is explained.

HyperOpt:
HyperOpt is a module proposed by Bergstra et al. (2013) , and is focused on intelligently searching through the hyperparameter space. One approach is to use the Tree-structured Parzen Estimator (TPE) algorithm (Bergstra et al., 2011), in which the search space is explored in an intelligent way, while the parameter values are narrowed down to the best estimated parameters. In contrary to the Grid Search, in which the hyperparameters must be pre-determined and the increment steps fixed, HyperOpt is an oriented random search and is proven to work efficiently (Bergstra et al., 2013). Hence, it serves as a good candidate to tune and optimize hyperparameters for universal approximators, as adopted in this paper.

Proposed Hyperparameter Optimization Algorithm:
• Initializing the sample space by HyperOpt. Suppose that we choose to start with N parameter settings to search more rigorously among them.
In our modified approach, we use HyperOpt to search intelligently through the search space, and we store the first N parameter settings that are used by HyperOpt. Note that each HyperOpt iteration is only performed on one set of train/test data. Now we use the generated N parameter settings as an input for a more rigorous search by successive halving (Jamieson and Talwalkar, 2016).
• We follow the idea of successive halving proposed by Jamieson and Talwalkar (2016). Using N parameter settings generated by HyperOpt, the well-known K-fold algorithm (Kohavi et al., 1995) is used to evaluate each parameter settings for a fixed amount of time/budget, e.g. T. Then, we select the top-performing half of the parameter settings (N/2), and again, we evaluate them via k-fold for time 2T. This procedure is repeated until the search space is singular or the designated budget is exhausted.
The idea behind the above algorithm is quite intuitive. Initially, HyperOpt is used as the screening procedure on the search region by expending a small In the second phase we add the best performing forecasts as additional features to the input data of the parent (See Figure 2). Next, a parent level forecasting model is re-estimated using the new input and then the hyperparameter optimization process is conducted. After identifying the best parameter settings, we select the best performing forecast as our final forecasting model.

Numerical Experiment
The MPH forecasting algorithm was implemented on sales data provided by MonarchFx Inc., which consists of 935 days of data for ten Stock Keeping Units (SKUs) and aggregated data which represents total brand sales. This data is representative of one of MonarchFx's mid-tier supply chain customers.
In addition to the historical sales data, the input also contained additional features including: • Promotion: a binary variable indicating if a promotion was present.
• Holiday: a binary variable indicating holiday periods.
• Day of the week: seven dummy variables corresponding to day of the week.
• Date: in the format day/month/year Each of these factors may increase the predictive power of forecasting models, both independently as well as in combinations.
To measure model accuracy, the Mean Absolute Error (MAE) is used: Where y i corresponds to the actual values of sales, and yˆi is the forecasted sales on day i.
The well-known k-fold cross-validation method (Kohavi et al., 1995) with k=5 is used to test each forecasting model and measure the forecast accuracy.
The parameter k refers to the number of groups that the input data will be divided to. We chose this method because it provides a less biased estimate of the model compared to a single train/test split of the data. In this procedure, initially the data is randomly shuffled and is divided into k different groups.
Then, for each group, it is selected as the test dataset and the remaining data is considered as the train set. The forecasting model is trained on the train set and the accuracy is measured on the test set. This procedure is repeated for every group and the average MAE across k train/test splits is reported as the final MAE.
Multi-Layer Perceptron (MLP), Random Forest (RF), Gradient Boosting (GB), and Extreme Gradient Boosting (XGB) are the forecasting models selected for the experiments, due to their popularity in machine learning forecasting literature (Mei et al., 2014;De'Ath, 2007;Cigizoglu, 2004;Deo et al., 2018;Zieba et al., 2016). Using these four forecasting models, the MPH algorithm in conjunction with the hyperparameter optimization method (explained in section 3-2) is implemented on the data and the forecasting error at the parent level is compared to the top-down and bottom-up approach.  The results of phase I are added as additional features to input data at the parent level. As phase II of the algorithm suggests, MLP, RF, GB, and XGB models are estimated again using the new input data.       is that MPH algorithm uses information at both levels of the hierarchy (SKU level and brand level), which helps the algorithm to provide significantly more accurate forecasts.

Conclusions
In this paper, we develop a novel multi-phase hierarchical approach (MPH) for supply chain forecasting using machine learning techniques supporting multifeature input data (e.g. MLP, RF, GB, and XGB). In the proposed two-phase model, the information at the child level is leveraged to improve forecasting accuracy at the parent level, by adding the results of the best forecasting model for each child as additional features at the parent level. The MPH algorithm is implemented on sales data provided by MonarchFx Inc. and the results were compared to a top-down and bottom-up approach. The results demonstrate that a considerable improvement can be achieved by utilizing the MPH algorithm (90% improvement in comparison with top-down approach, and 82% improvement comparing to bottom-up approach). This improvement is possible due to the fact that the MPH algorithm leverages information both at the child level and parent level.
Based on the experience of one of the co-authors who leads the supply chain analytics function at MonarchFx Inc., there are multiple applications possible for our approach. Indeed, the majority of companies employing supply chain forecasting solutions generally apply top-down and bottom-up approaches and use traditional models that only support single feature input data. However, in practice, multiple factors can impact future sales and can be controlled for in this manner to improve forecast accuracy. Using the machine learning forecasting models developed in this paper, supply chain planners can derive more accurate forecasting models to exploit the benefit of multivariate data.
There are multiple possible future extensions to this work. One is to use the MPH algorithm on hierarchies with more than two levels. The other is to utilize the reconciliation techniques used by Hyndman et al. (2011). Another possible path is to characterize the situations under which each of the forecasting models perform best in the child level and parent level. Researchers may also use the forecasting model selection method developed by Taghiyeh et al. (2020) to select the best forecasting model among existing machine learning models.
The noisy optimization method in Taghiyeh and Xu (2016) may also be used to find the parameters for optimal reconciliation for the levels of the hierarchy. To improve the speed of the model, the parallelization method proposed in Rosen et al. (2016) can be utilized. We believe there lies great promise for using these approaches in the future.