Analyzing Collaborative Forecast and Response Networks

Partial Correlation refers to the correlation between two random variables ($X$ and $Y$) when the effect of a set of external variables ($Z$) are removed from both of them. We will denote it as $C(X,Y|Z)$, where $C$ refers to correlation. Pearson’s product-moment correlation coefficient is used as a measure of correlation. See Baba et al (2004) and Fisher (1924) for more information.

A good insight about partial correlation comes from a very simple example. Let the random variable $X$ represent ice cream sales over time and $Y$ sunscreen sales over time. It is reasonable to expect a high correlation between these two variables. However, this high correlation does not necessarily mean that one of these events is causing the other. In fact, if a variable expressing the season of the year ($Z$) is introduced, a most likely result will be that $C(X,Y|Z)$ is very small. This means that both $X$ and $Y$ (ice cream and sunscreen sales) are simply reacting to a third variable (summer). When the effect of summer is accounted for, $X$ and $Y$ are not related any more.

While the theoretical partial correlation refers to the theoretical relationship between random variables, the empirical partial correlation is a number estimated using observed samples from these random variables. Empirical partial correlation is an estimator for the theoretical partial correlation.

One way to find the partial correlation between $X$ and $Y$, given $Z$, is to use linear regression:

These regressions allow factoring $Z$ out of both $X$ and $Y$. The correlation between the residuals $\varepsilon_{X|Z}$ and $\varepsilon_{Y|Z}$ can then be used to find the empirical partial correlation: $C(X,Y|Z)=C(\varepsilon_{X|Z},\varepsilon_{Y|Z})$.

An equivalent approach to the partial correlation is to compute the empirical covariance matrix $S$ of the random variable set $[X,Y,Z]$. If $S$ is invertible, then the entry $S^{-1}_{1,1}$ (The 1,2 entry of the inverse covariance matrix) is proportional to the partial correlation between $X$ and $Y$, given $Z$. In the literature, the inverse of the covariance matrix is called the precision matrix, or concentration matrix, and is often denoted with $\theta$.

In a random variable set of arbitrary length, the entry $\theta_{ij}$ of the precision matrix is proportional to the partial correlation between variables $i$ and $j$ when the effect of all other variables in this system are removed from both them:

Partial correlation is a very powerful tool in isolating the interaction of two variables from the effects of the larger system they are in. Given a data set of many variables, the precision matrix of this system provides insights into the isolated interactions of any two variables. The larger entries in this matrix point to strong partial correlation, while smaller entries mean little or no correlation. A very important question at this point is to decide whether any small entry refers to weak but existing partial correlation, or if it is just a noise artifact. We approach this using sparse precision matrix methods. These assume there are relatively few important underlying partial correlations, with the rest being practically zero. The goal is to identify and estimate the non-zero partial correlations. See Dempster (1972) and Cox and Wermuth (1996) for more information.

When the variables in such a system come from a multi-variate Gaussian distribution, the problem of finding a sparse precision matrix in such a system is solved through a Gaussian Graphical Model (GGM). This is a well-researched area in the literature. See Yuan and Lin (2007) and the references therein for more information. Among other uses, they are extensively used to explore the interactions between genes in genetic studies. (Dobra et al (2004) ). These models are built to understand the interactions in a network of variables, where these variables are assumed to be following a multi-variate Gaussian distribution $N(\mu,\Sigma)$. A GGM network is represented by an undirected graph $G=(V,E)$, where $V$ is the set of vertices, and $E$ is the set of edges between the vertices. Each vertex corresponds to a random variable in this network. The edges describe the conditional independence between the random variables. An edge between two variables $X$ and $Y$ is absent whenever $X$ and $Y$ are essentially independent when conditioned on the other variables in the network.

For variables $X$,$Y$,$Z$ following a multi-variate Gaussian distribution, the conditional independence of $X$ and $Y$ given $Z$ is equivalent to $C(X,Y|Z)$ being zero. Therefore finding a sparse precision matrix on this network leads to understanding the dependencies between the nodes.

The GGM problem refers to estimating the existence or non-existence of the edges of the graph $G$, given a random sample from this network. Directly calculating the partial correlations from the samples will not give the desired results: because the input is a random sample, zero partial correlation will almost never be observed. However, values close to zero will indicate that the corresponding vertices are most probably independent, and the small non-zero partial correlation is just noise. Therefore, the challenge of GGM problem is to identify which entries of the precision matrix are actually zero with a certain level of confidence. This problem has been named the covariance selection problem, or the model selection problem in the Gaussian Graphical Network.

For any GGM method, a parameter $\lambda$ is needed to control the tightness of the sparse model selection. A higher $\lambda$ will result in a network where many edges are pushed to zero. As $\lambda$ is lowered, the model selection will be more relaxed, and links with less certain conditional dependencies will appear. An appropriate $\lambda$ should be determined based on the data and the needs of the analysis.

The literature contains many methods proposed to estimate the precision matrix given a desired tightness for model selection. In this work, the graphical lasso method proposed by Friedman, Hastie and Tibshirani in 2008 is used. This method is shown to be robust to data deviations, and it is computationally very fast.

The graphical lasso method solves the following problem for any given Gaussian graph $G$ and tightness parameter $\lambda$:

Where $S$ denotes the empirical covariance matrix. In this formulation, the first two components represent the log likelihood function for the precision matrix $\theta$, and the last component is the linear shrinkage penalty to obtain sparseness.

Other approaches to the sparse GGM problem exist in literature. Some recent examples are explained in Banerjee et al (2008), Dobra and West (2004) and Meinshausen and Bühlmann (2006).

Some other studies in literature seek to relax the normality assumption of GGM’s and therefore a provide more general framework. See Xue and Zou (2012) for a nonparanormal model, and Dobra and Lenkoski (2004) for a copula representation.

The general GGM model does not consider a time dimension. The vertices do not represent events that happen in a sequence. To find the conditional dependence/independence between two vertices $X$ and $Y$, GGM controls for all the other random variables in the network. Another consequence of the lack of a time-line is that the conditional dependence discovered between $X$ and $Y$ does not immediately provide information on causality. We do not know if one variable may have caused the other; we only know there is dependence.

An information flow network like the subject of this paper, however, has an inherent time dimension. All events in the collaborative forecasting network happen in a sequence. This fact has two consequences.

First, when searching for the partial correlation between two events $X$ and $Y$, it is not appropriate to condition them on all the other events. The appropriate action is to condition them only on the events of the past. The events that happened after both $X$ and $Y$ should not influence their dependence.

For deeper insight, suppose that this principal is not followed in finding the dependencies of this network. as an illustrative example, also suppose that the information flow network contains three consecutive events $X$, $Y$ and $Z$, and there is an information signal that entered the system with X, and then flowed to $Y$ and then $Z$. When investigating the flow between $X$ and $Y$, this information should be represented as a link between these two nodes. To apply the classical GGM methods to this network, the partial correlation between $X$ and $Y$ conditioned on $Z$ is calculated. Since $Z$ contains future information, this would result in inappropriate conditioning, where $X$ and $Y$ may appear to be conditionally independent.

Secondly, the time dimension allows us to construct our graph so that when a conditional dependence is discovered between $X$ and $Y$, the information contained in the earlier event ($X$) is propagated to the latter ($Y$), and not the other way around. Therefore we use a directed graph rather than the classical undirected graph of a GGM.

To accommodate the time dimension of our forecasting network, we first develop some appropriate notation, and then propose a scheme called Expanding Window GGM, or EW GGM.

Recall that $V$ denotes the set of nodes in this network. Let:

where $v_{n}$ refers to the $n^{th}$ event in the network. Notice that the indices indicate the sequence of events. Therefore $V_{n}$ refers to the set of first $n$ events. Furthermore, let $C^{n}$ denote the partial correlation matrix of the event set $V_{n}$.

In the information flow network, the partial correlation between any two nodes $i,j$ conditioned on all the events that happened before $i,j$ or between $i$ and $j$, is of interest. Let the matrix of all desired partial correlations be denoted by $C^{\prime}$. The $i,j^{th}$ entry of $C^{\prime}$ is:

Then:

Therefore:

where $C^{i}_{i}$ is the column $i$ of the partial correlation matrix $C^{i}$, as defined above.

The problem of finding a sparse approximation of the $C^{\prime}$ matrix is a difficult one. Note that, a simpler problem of finding a sparse precision matrix for the classical GGM problem was solved in a computationally efficient manner only recently. In this paper, instead of working towards an optimal solution, we will propose a heuristic approach that provides a robust approximation to the desired $C^{\prime}$ matrix. We will also take advantage of the fast solution that is available for the classical GGM problem.

The Expanding Window GGM method relies upon providing sparse approximations for the columns $C^{i}_{i}$ through the sparse precision matrices $\theta^{i}$’s, calculated using the graphical lasso method for all $V_{i}$’s in the network. An illustration of the expanding windows of this idea is given in Figure 2.

Letting $S_{k}$ denote the empirical covariance matrix of the event set $V_{k}$, this heuristic method can be summarized as follows:

In this system, each event is a combination of information that is available from the past periods, and the new information obtained in that period. The EW GGM method of the previous section provides a way to determine which past events contributed information to any given event. The next step is to determine the proportion of these contributions. Although the magnitudes of partial correlations expressed in $\hat{C}^{\prime}$ give an indication of the size of the contributions, another step is needed for a precise decomposition.

Going back to the notation of Section 4, the events consist of forecasts (F), responses (R) and a shipment event (S). We assume that the information in each event can be decomposed linearly.

Figure 3 presents the intended design of the rolling horizon forecast collaboration system. Every forecast and response should include the signals that entered the system before it. Furthermore, we will also investigate the Markovian property. This corresponds to the situation where the effect of the past on $F_{i}$ is fully captured by the events $F_{i+1}$ and $R_{i+1}$, and any previous event does not carry additional information. Therefore, the effect of, say, $F_{i+2}$ on $F_{i}$ should be close to zero when the effects of $R_{i+2}$, $F_{i+1}$ and $R_{i+1}$ are accounted for through conditioning.

Collaborative forecasting is intended to be Markovian, so it is interesting to investigate this property for the actual forecasting data. The EW GGM method of the previous section provides the actual links that are observed in this network. Based on the linear decomposition assumption, the following equations hold:

For lag $N$:

For lag $i=1…(N-1)$:

Where $F_{i:N}=[F_{i},\ldots,F_{N}]$.

For S:

The coefficients in this system are calculated using regression.

Recall that the empirical partial correlation matrix of the network has zero entries for events that are partially independent. In the decomposition step, the coefficients that correspond to a partially independent pair of events are set to zero. In other words, the information in each event is decomposed using only the past events from which an incoming arrow to this event exists. Therefore the EW GGM method serves as a model selection tool in decomposition. Note that the sparsity of the estimated partial correlation matrix is controlled by the parameter $\lambda$. Higher values of $\lambda$ will result in tighter model selection in this step.

The non-zero $\alpha$ and $\beta$’s obtained through regression correspond to the magnitude of the influence of past events, while the $\varepsilon$’s correspond to the new information and noise.

Standardizing the input data is a common practice in statistics. It prevents scaling issues from dominating the results, and allows the correct detection of correlations. Another common practice is to test for the normality of input data, and apply a monotone transformation to the data if its distribution is noticeably different from normal. The transformation allows the analysts to take advantage of the attractive properties of the normal distribution. Moreover, hypothesis testing (checking whether an observed relationship is actually there or not) can be reliably done with appropriately transformed input.

Our method assumes Gaussian input data. Unfortunately, there tends to be strong skewness in collaborative data. Therefore we use recommend use of the Box-Cox transformation of Box and Cox (1964):

The proper $\gamma$ is found through exploration. In our data set, a value of $\gamma=-0.5$ seems to achieve reasonably Gaussian data for some representative parts. The plots of empirical distributions obtained using this value of $\gamma$ for one example is shown in Figure 4.

After normalization, as the second step, we standardize the transformed input. Notice that the equations given in the previous subsection do not contain a constant term, because all columns are scaled to have zero mean.

We propose a new visualization technique to represent the decomposition of the information flow network. Figure 5 shows an example forecasting network illustrated in this fashion.

In Figure 5, the grey horizontal line represents the time-line. Each of the forecast, response and shipment events is shown as a node. The $X$ coordinate of the position of each node indicates the time point that the event happened on the time-line. The vertical grey lines from the nodes can be used to identify when each event happened. The $Y$ coordinates of the nodes are determined such that the nodes are arranged on an ellipse. This allows us to present the edges with least amount of overlap. The forecasts are arranged on the top hemisphere of the ellipse and responses are on the lower hemisphere to enable quick visual inspection of forecast-response interaction. The shipment node, if it exists, sits on the rightmost end of the time-line.

An edge from node $i$ to node $j$ is drawn if the EW GGM method indicates that these events are conditionally dependent. The directions of these edges are always from left to right (from past to more recent). The color and thickness are determined by the decomposition regression. For each node, the color of the incoming edge is blue if the corresponding parameter in that node’s regression is positive, and it is red if the parameter is negative. The thickness of the arrow is proportional to the magnitude of the parameter, indicating which past events exerted large and small influences.

The methods described are applied to the data set coming from the collaborative forecasting practice between HP and a microprocessor provider. The data contains forecasts with up to $12$ lags, and responses with up to $7$ lags.

In this section, we will illustrate the use of this tool on an example unit from this data set. The selected unit is a microprocessor used in North America region. The unit is selected due to its strategic importance and high volume. The observations obtained from this unit are typical among a larger collection that was considered. Thus the similar results for other representative units are presented in supplemental files and not in the main text.

The first decision the user will make is the choice of the model selection parameter. The forecast analytics software tool allows users to interactively change this parameter and review results in real time. For the purposes of this paper, the networks obtained by changing the model selection parameter in the range $[0.7,1]$ are provided in a supplementary document. A quick exploration of this document using the up and down arrow buttons allows the user to observe the effect on the application. In our experience, the parameter range around $0.7-0.8$ gave too little model selection, i.e. way too many very small links were shown while $0.9-1.0$ was too tight, erasing most of the links. Approximately $0.8-0.9$ gave the most reasonable range. In fact, $0.8$ seems to be a reasonable value for most of the units in the data set. The forecast-response network of our example part for this parameter value is in Figure 5.

Recall that we expect to see the Markovian property in this network. If it is there, past influence on each event should be fully captured by the most recent forecast and response. It would visually present itself by short, fat arrows coming in to each event only from the most recent events. In Figure 5, we see that this property is mostly there for within forecasts and within responses.

A striking observation in this graph is the lack of arrows going from response nodes to forecast nodes. This indicates that the buyer is not strongly using the information coming from the supplier when revising the forecasts. The adoption of the forecasts based on the supplier’s signals is an important goal of this setting, so our visualization provides an important diagnostic tool, which gives insight into how well this system is working in practice.

The next observation is, the forecasts are closely linked with each other, and the responses are closely linked with each other, as expected. On the other hand, the only influence between forecasts and responses appear from $F5$, $F4$, $F3$ and $F2$ to the responses $R5$, $R4$, $R3$ and $R2$, respectively. This makes sense, because for this particular microprocessor, the production delay is $4$ periods. Therefore, when deciding how much to produce $4$ periods ahead, the most recent demand forecast that the supplier has access to is $F5$. The supplier does not react to earlier forecasts because they are not recent at the time of production decision. The updates communicated in $F4$, $F3$ and $F2$ are also incorporated to supplier’s response through supplier’s manufacturing flexibility, buffer inventory and re-arrangement of product commitments. However, $F1$ carries a very late update, perhaps too late, so the supplier does not react to it.

We see that the events $F12$, $F11$, $R1$ and $S$ are not linked to any other event. $F12$ and $F11$ are the earliest events. It may be the case that it is yet too early to make meaningful projections for weekly purchase, and thus these nodes do not carry significantly meaningful information. Similarly $R1$ is the response right before the realization of the purchase. It is not meant to influence any forecasts. Knowing this, the responders do not seem to put meaningful information in it. Given that these nodes do not provide meaningful information, and they take effort to produce, it may be beneficial to remove them from the process to save forecasters’ time. Thus our visualization tool provides a large amount of useful information for streamlining the forecast response system.

The fact that the event $S$, shipment, is not connected to any other event is an indication of forecast inaccuracy in the system. This may be the most important problem to fix in this system, since forecast inaccuracy against the actuals is the biggest driver of buffer inventory and stock out costs.

Finally, recall from Section 3.3 that the magnitude of any influence is communicated through the thickness and color of arrows. When it is more useful to know the exact values in the underlying decomposition, these will be given using a mouse-over action on the tool. As an example, the decomposition of the $F1$ event in this network is:

This decomposition implies that $50\%$ of the information contained in $F1$ propagated from $F2$, $40\%$ came from $F4$, and $10\%$ is new information and noise.

In this section, we will illustrate the CCC method on a microprocessor part. In the ANANSI tool, the method is implemented such that the user can interactively change the $\alpha$ parameter and see how both the summary plot and the forecast-response scatter plot change. In the paper, we will provide these plots for the $\alpha$ values $0$, $0.5$ and $1$.

In Figure 6, we see the results of applying CCC to this part.

The first row shows the plots with $\alpha=0$. This means the method only seeks correlation structure. On the left panel, we see that a very strong correlation exists in later periods, which indicates close collaboration and coordination between the parties. In earlier periods, although some correlation exists, it is not as strong. This is most probably due to the fact that in the initiation periods, supplier kept a certain level of inventory, but the buyer orders weren’t strong yet, so there wasn’t a strong need for close collaboration.

Similarly, the scatter plot on the right hand side indicates a very strong linear trend, with cooler colors (earlier periods) occupying the lower end of order sizes, and warmer colors (later periods) make up the higher order size periods. An outlier in this plot would indicate a period where the forecast and response values in all horizons were unrelated to each other, suggesting a possible communication problem or a maybe a disagreement. Such a period does not exist in the given example.

Even though some insights can be gained from these graphs, it is clear that both plots show an excessive amount of correlation, due to overfitting.

In the second row, we have $\alpha=0.5$. This option considers both the variance and correlation equally. We see that the summary vectors of forecasts and responses very closely follow their $\alpha=1$ vector summaries, which only considers the variance. This means that the within-sets variability for this data set is much stronger than the correlation effect, causing these vectors to feel the variance more strongly than the correlation. Note that this conclusion can only be made when $\alpha=0.5$ results are compared to the other two results.

In the third row, we have $\alpha=1$, which gives the summaries of forecasts of responses without taking into account the relationship between them. It only considers to retain the information existing within the forecasts (blue line), and separately, the information within the responses (red line).

The summary forecast and response series show a low level of volume in the first half of the periods, and a high and rather volatile forecasting pattern for the second half of the periods (Left panel). This suggests that this part was initiated around periods 1-10 with low volumes (and sufficient supplies), and the consumption picked up in the later periods (where supply predictions were slightly lower for certain periods). The trends of forecasts and responses are close to each other, indicating that a long term strong supply shortage is not in effect. Although these summary vectors only consider the within-set variance, the second half of the periods show a large amount of correlation. This can also be seen in the scatter plot of the right panel. A linear trend in this plot would suggest that responses and forecasts are moving together. Like with the other $\alpha$ values, the blue-green hue dots crowd the lower left area of the scatter plot, while reds are closer to upper right. This tells us that high predictions exist in later periods, while the earlier periods have relatively less volume. Outliers from the linear trend would suggest a big shortage period (this would show up in lower right) or an excessive inventory period (this would show up in upper left). Our example does not possess strong outliers.

If a certain default value for $\alpha$ to be followed is desired, these plots make it clear that a good balance between the variance and correlation can be obtained somewhere between $\alpha=0$ and $\alpha=0.5$, Interactively going through all the possible $\alpha$ values (not given here), we see that a value around $0.1$ provides an acceptable balance.

For future, it may also be beneficial to approach this method with a mechanism design perspective. A certain value of $\alpha$ acceptable to both buyer and supplier can be decided on to follow, and deviations and fluctuations defined over the summaries based on this $\alpha$ can be penalized. The effect of such an agreement on the forecast accuracy and collaboration, and the best values of $\alpha$ to achieve supply chain goals of interest are left to future research.