TsSHAP: Robust model agnostic feature-based explainability for univariate time series forecasting

In this paper, we propose a model agnostic feature-based explainability method, named TsSHAP, that can explain the prediction of any black-box forecasting model in terms of SHAP value-based feature importance scores (Lundberg and Lee, 2017). The novelties of the proposed method are the following.

• •

  TsSHAP can explain the prediction(s) of any black-box forecaster model from a set of interpretable features defined by the user a priori. It only needs to access the outputs of the fit() and predict() functions of the model being explained¹¹1Following conventional machine learning terminology, fit() and predict() refer to the training and inference functions of the forecasting algorithm respectively.. It always needs to access the predict() function, and fit() is required only for a specific set of forecasters (see § 4.2 for details).

• •

  It amalgamates the notion of model surrogacy, and the strength of the fast TreeSHAP algorithm into a suitable and ready-to-use post-hoc explainability method for time series forecasting.

• •

  The novel TsSHAP methodology uniquely reduces the surrogate forecasting task into a regression problem where the surrogate model learns a mapping between the interpretable feature space and the output of the forecaster model.

• •

  It introduces multiple scopes of explanation, viz. local, semi-local, and global, which can be useful to obtain useful explanations in several real-world scenarios.

Let $f(t):\mathbb{Z}\to\mathbb{R}^{1}$ represent a latent univariate time series for any discrete time index $t\in\mathbb{Z}$. We observe a sequence of historical noisy (and potentially missing) values $y(t)$ for $t\in[1,\dots,T]$ such that in expectation $\mathbb{E}[y(t)]=\mathbb{E}[f(t)]$. For example, in the retail domain $y(t)$ could represent the daily sales of a product and $f(t)$ the true latent demand for the product.

The task of time series forecasting is to estimate $f(t)$ for all $t>T$ based on the observed historical time series $y(t)$ for $t\in[1,\dots,T]$. The time series forecast is typically done for a fixed number of periods in the future, referred to as the forecast horizon, $H$. We denote

as the projected forecast over time horizon $H$ based on the historical observed time series $y(1),...,y(T)$ of length $T$.

Explainability is the degree to which a human can understand the cause of a decision (or prediction) made by a prediction model (Miller, 2019). An explanation can be based on an instance, a set of features, or (for a time series) inherent components (e.g., trend, seasonality etc.). In this paper, our focus is feature-based explanation.

A feature-based explanation $\phi\in\mathbb{R}^{d}$ is a $d$-dimensional vector where each element corresponds to the importance /contribution score for the corresponding feature. Several feature-based explanation methods have been proposed in the supervised classification and regression literature (see § 3). We will adopt SHAP for explanation because of it’s intuitive nature (Lundberg and Lee, 2017) and state-of-the-art performance (Lundberg et al., 2020).

TsSHAP is model agnostic. Since we do not have access to the actual internals of the forecaster model, we first learn a surrogate model $g$.

The surrogate model $g$ produces a point-wise forecast approximation of the forecaster. While the original forecaster learns to predict $y(T+h)$ based on $y(1),\ldots,y(T)$ the surrogate model is trained to predict the forecasts $\hat{f}(T+h)$ made by the forecaster. Essentially we want to mimic the forecaster using a surrogate model since we aim to explain the forecasts made by the forecaster and not the original time series. Figure 1 illustrates this concept with an example time series from a real-world dataset.

To generate the training data for the surrogate model, we use backtesting. For a given univariate time-series, we produce a sequence of train and test partitions using an expanding window strategy (see Figure 2). The expanding window partitioning strategy uses more and more training data while keeping the test window size fixed to the forecast horizon. The forecaster is trained (fit() is accessed) on the train partition and evaluated (predict() is accessed) on the test split. All the test split predictions are concatenated to get the backtested historical forecasts for each step of the forecast horizon.

Note that we need to access the fit() function only for a set of classical forecasters like SARIMA or exponential smoothing, where we train the forecaster for every expanding window. For most of the machine learning and deep learning forecaster models, we do not need to access their fit() functions because a single model is generally trained (beforehand) on the entire data.

We now have an original time series and the corresponding backtested forecast time series for each step of the forecast horizon. The goal of the surrogate model is to learn to predict the backtested forecast time series based on the original time series. We construct a surrogate time series forecasting task by reducing it to a supervised regression problem. A standard supervised regression task takes a $d$-dimensional feature vector as input ($\mathbf{x}\in\mathbb{R}^{d}$) and predicts a scalar $y\in\mathbb{R}$. The regressor $y=f(\mathbf{x})$ is learnt based on a labelled training dataset $\left(\mathbf{x}_{i},y_{i}\right)$, for $i=1,\ldots,n$. The input time series sequence does not satisfy classical feature definitions. Instead, we process the original time series to compute various interpretable features.

For each time point $t$, we compute features $\mathbf{x}(t)\in\mathbb{R}^{d}$ based on the original time series values observed so far. The surrogate model uses these $d$ features to predict the backtested forecast for each step in the forecast horizon. We focus only on interpretable features such that the explanations produced on the features are easily comprehensible by the end-user. The set of interpretable features can be defined by the (expert) user as well. As a starting point, we define a set of seven feature families, viz., lag, seasonal lag, rolling window, expanding window, date and time encoding, holidays, and trend-related features. Table 1 and 2 list the features.

Recall that we have $H$ backtested time series forecasts which are to be used as targets to learn the surrogate regressor model. A single regressor is trained for a one-step-ahead forecast horizon and then called recursively to predict multiple steps. Let, $\mathcal{G}(y(1),...,y(T))$ represent the one-step-ahead forecast by the surrogate model trained on the training data. The surrogate produces a one-step-ahead forecast using features based on the time-series till $T$, that is, $y(1),\ldots,y(T)$. The final forecasts from the surrogate model are made recursively as follows,

For $h=2,3,\ldots,H$

In this work, we have used a tree-based ensemble regression model, XGBoost (Chen and Guestrin, 2016)⁴⁴4https://github.com/dmlc/xgboost, as the surrogate model⁵⁵5Other models like CatBoost (Prokhorenkova et al., 2018) and LightGBM (Ke et al., 2017) can also be employed.. The selection of tree-based model is biased by its reasonably accurate forecast capability and, the support of the fast (polynomial time) TreeSHAP algorithm which we rely on for the various TsSHAP explanations as described below. (Lundberg et al., 2020)⁶⁶6https://github.com/slundberg/shap.

We introduce the notion of scope of explanation in the context of time series forecasting. An explanation answers to either a why-question or a what if-scenario (Molnar, 2019). We generalize three scopes of explanations.

1. (1)

   A local explanation explains the forecast made by the forecaster at a certain point in time. For example,

   • •

     Why is the sale forecast for July 1, 2019 much higher than the average sales?

   • •

     What is the impact on the sale forecast for July 1, 2019 if I offer a discount of 60% that week?

2. (2)

   A global explanation explains the forecaster trained on the historical time series.

   • •

     What are the most important attributes the forecaster algorithms rely on to make the forecast?

   • •

     What is the overall impact of discount on the sales forecast?

3. (3)

   A semi-local explanation explains the overall forecast made by a forecaster in a certain time interval.

   • •

     Why is the sales forecast over the next 4 weeks much higher?

   • •

     What is the impact of a particular markdown schedule on a 4 week planning horizon sales forecast?

Notions of local and global explanations exist in the supervised learning paradigms. In addition, we introduce the concept of semi-local explanations in the context of time series forecasting, which returns one (semi-local) explanation aggregated over many successive time steps in the forecast horizon. TsSHAP is able to provide explanations suitable for all the scopes.

For each of the scopes, TsSHAP provides multiple types of explanation.

We experimentally validate the proposed explainer algorithm and the corresponding various explanations on five univariate time series forecasting datasets (with any available external regressors) listed in Table 3. For all datasets, we use a forecast horizon of 10% of the actual data. Availability of data is provided in the Reprodicibility section (§ 9).

In order to validate the explanations we present results on six different univariate forecasting algorithms listed in Table 4. Only some forecasters (Prophet (Taylor and Letham, 2018) and XGBoost (Chen and Guestrin, 2016)) can explicitly incorporate external regressors.

Robustness of TsSHAP explanations is evaluated by perturbing the original time series (in § 6.4, we will describe the evaluation metrics that will use the perturbed time series). To generate bootstrapped samples (i.e., multiple perturbed versions) of a time series $y(t)$, we employ the block bootstrap algorithm (Bühlmann, 2002; Kreiss and Lahiri, 2012). The block bootstrap algorithm extends Efron’s bootstrap (Efron, 1992) setup for independent observations to time series where the observations ($y(t)$ for different values of $t$) might not be independent.

To preserve the trend-cycle of the time series (Hyndman and Athanasopoulos, 2018), we first decompose $y(t)$ into trend-cycle and residual components using a moving average model of order $m$,

The residual series $y_{\text{residual}}(t)$ is then utilized by the block bootstrap algorithm. It samples contiguous blocks of $y_{\text{residual}}(t)$ at random (with replacement) and concatenates them together to produce a bootstrapped residual,

where, $L_{b}$ is the block-length. The final perturbed series is then obtained by summing the bootstrapped residual to the trend-cycle component of the original series,

We have extended three well-adopted evaluation metrics in the context of time series, viz. faithfulness, sensitivity, and complexity. Faithfulness measures the consistency of an explanation, sensitivity measures the robustness, and complexity measures the sparsity of an explanation (Bhatt et al., 2020).

We first validate the accuracy of the surrogate model forecasts from the explainer via backtesting (temporal cross-validation). Using an expanding window approach, we partition the univariate time series into a sequence of train and test datasets. The forecaster and the explainer are trained on the train partition and then evaluated on the test dataset. We report various error metrics between the forecasts from the original forecaster and the forecasts made by the surrogate model of the explainer, viz. MAE (Mean Absolute Error), RMSE (Root Mean Squared Error), MAPE (Mean Absolute Percentage Error), and MASE (Mean Absolute Scaled Error) (Hyndman and Koehler, 2006). Table 5 (in § 9 Reproducibility) shows the accuracy scores between the forecaster and the surrogate model for the six forecasters and five datasets.

Figure 3 shows the variation in MAPE with increasing size of the datasets and increasing complexity of the models. It is evident that the MAPE increases with an increase in the model complexity. This might be because it is harder for the surrogate model to mimic a complex forecaster (e.g. XGBoost) than a simpler one (e.g., Naive). However, there is no such pattern in the variation of MAPE with increase in the dataset size. This is particularly beneficial since it conveys that the dataset size does not affect the surrogate model training.

We evaluate the feature-based explanations in terms of the (high) faithfulness, (low) sensitivity and (low) complexity metrics defined in § 6.4. In the Reproducibility section, § 9, we provide various metrics for global, semi-local and local explanations for all the forecasters and the datasets in Table 6, 8 and 7. The metrics for the local explanations were averaged over the entire forecast horizon.

Figure 4 provides a graphical summary of Table 6, 8 and 7 by aggregating the metric values across all the datasets⁷⁷7For sensitivity, we do not show mean (std) because its value depends on the absolute value of the data (see § 6.4). Instead, we show median with 1st and 3rd quartiles as the error bars.. We can observe that, for all the forecasters, the faithfulness of semi-local and local explanations are higher than that of global explanations. A closer look revelas that the faithfulness of a semi-local explanation is always higher than the faithfulness of a global explanation. However, when we move to local explanations (from semi-local), there is no definite pattern in change in the faithfulness metric. The complexity of the explanations do not vary much with the scope of explanation, as can be seen in Figure 4(c). This might refer to an implicit relationship between the complexity of the forecast explanation and that of the forecaster. The median sensitivity curves stay (approximately) flat for all the forecasters.

Moreover, we can notice that a particular forecaster might not achieve the best performance in all the metrics, and this highlights the complementary nature of the metrics defined in § 6.4. For example, in Table 6, for “jeans-sales-weekly” dataset, MovingAverage forecaster achieves the highest faithfulness and the lowest sensitivity, but the lowest complexity is obtained by the Naive forecaster. This type of behaviors is also observed in local and semi-local explanations (Table 7 and Table 8).

A closer look on Figure 4 (and on Table 6, 7 and 8) reveals that the MovingAverage forecaster achieves the average highest faithfulness (0.244 for global, 0.654 for local, and 0.658 for semi-local explanations) across all datasets. The Naive forecaster obtains the average lowest complexity (0.074 for global, 0.086 for local, and 0.088 for semi-local explanations) across all datasets. This might be because the inherent simplicity of the Naive forecaster results in the fractional distribution of the feature importance scores to have low entropy (see § 6.4).

However, the median lowest sensitivity is shared among MovingAverage and SeasonalNaive forecasters across different scopes of explanations. The MovingAverage forecaster has the lowest median sensitivity for global (0.8) and local (12.64) explanations, but the SeasonalNaive attains the lowest value (7.7) for semi-local explanation.

Another noticeable observation from our experiments is that the complex (and possibly non-linear) models (like Prophet and XGBoost) have high complexity, which is expected. However, they might achieve better faithfulness and sensitivity than a simpler model. For example, in Table 6, Prophet attains higher faithfulness than a Naive forecaster, and in Table 8, XGBoost achieves lower sensitivity and higher faithfulness than a Naive forecaster in multiple datasets.

Figure 5, 6, and 7 illustrate a few local, semi-local, and global explanations for the jeans-sales-daily dataset for all different forecasters. Note that the PDP and SDP plots are only shown for a local scope. It can be seen that the explanations reflect the underlying aspects the forecaster is able to learn from the data. From the local, semi-local, and global explanations of Figure 5, we can see that the main contributing factor for the forecast is the seasonal history at $(t-52)$ weeks time, which is expected from a seasonal naive forecaster. The partial dependence plot also depicts a linear relationship between the sales(t-52) and forecaster output. From the local explanation of Figure 6, we can see that the mean value of the past 6 observations is primarily responsible for bringing the forecast down. Other explanations depict the same story. From the local explanation of Figure 7, we can see that, mainly, discount(t) and is_quarters_start are pushing the forecast up, but year and sales_max(t-1, t-6) are pushing the forecast down. The Prophet forecaster supports external regressors, and here, it is picking discount as the primary explanatory factor. The global explanation clearly shows the main contributing features responsible for the forecast across the entire time series dataset. Some other illustrations are provided in § 9.

1. (1)

   jeans-sales-daily: Proprietary dataset.

2. (2)

   jeans-sales-weekly: Proprietary dataset.

3. (3)

   us-unemployment:
   https://www.bls.gov/data/#unemployment

4. (4)

   peyton-manning:
   https://en.wikipedia.org/wiki/Peyton_Manning

5. (5)

   bike-sharing: http://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset

We have provided a summarized observation from our experiments in the main text (§ 7). Here, we demonstrate the detailed results obtained in the experiments.

Figure 8 shows the TsSHAP explanations for XGBoost forecaster on the jeans-sales-weekly data. We can see that the surrogate model accurately approximates the forecaster in this case. Moreover, a closer look into the local explanation says that sales(t-1), last year’s value sales(t-52), and the external regressor discount(t) are primarily bringing the forecast up from the average sales. The global explanation shows that sales(t-1) and discount(t) are the main explanatory factors for the forecast made by the forecaster across the entire time series data. The semi-local explanation depicts similar picture for the entire forecast horizon interval. The specific PDP shows that after a certain threshold, increasing discount is helpful for increasing sales, but it can also saturate at some point. The specific SDP shows sales(t) is positively correlated with sales(t-1).

TsSHAP also provides correct explanations for the Naive forecasters, i.e., it picks the last observation as the most important feature in all the scopes. We skip the illustration for the Naive forecaster due to space limitation.