Explainability is a crucial aspect of machine learning that allows us to understand and trust the predictions made from a model. We not only predict the next value, but also explain why the model output that particular value.

This is an essential step to trust the output of a model, but also to understand how different features impact the final outcome. With explainability techniques, we can answer questions like:

• how will a discount impact the sales of my product?
• will a holiday increase or decrease traffic to my website?

Until recently, explainability in time series forecasting was not implemented at scale for deep learning models, which often have the reputation of being complex black box models.

Now, with my latest contribution to the neuralforecast library, an open-source Python package to forecast time series data with state-of-the-art deep learning models, explainability methods are available at scale to attribute the impact of exogenous features and past values of a time series in deep learning models.

In this article, we explore the different methods available for explaining deep learning forecasting models, and we complete a small experiment using Python.

Let’s get started!

Learn the latest time series analysis techniques with my free time series cheat sheet in Python! Get the implementation of statistical and deep learning techniques, all in Python and TensorFlow!

Discover explainability methods

Currently, we implemented three explainability methods for deep learning models, which are summarized in the table below.

From the table above, we first notice that all explainers are local. This means that we get explanations for how a certain input impacts a specific forecast step.

Then, we have two explainers that are based on gradient computation: integrated gradients and input X gradient. The other explainer, Shapley value sampling, is based on feature perturbation.

Let’s explore each method in more detail.

Gradient-based methods

As the name suggests, these methods rely on gradient computation to calculate feature attributions. Methods like integrated gradients and input X gradient are part of that category.

Integrated gradients

With this method, we approximate the integral of the gradient of the model from a straight-line path from a baseline to the model’s output [1]. Typically, the baseline of a model is an input of 0 values. 

To approximate the integral, we use a Riemann sum, and it requires between 20 and 300 gradient computations.

This method has many advantages, the main ones being that it is fast and it has the additivity property. This means that the model’s forecast is the sum of feature attributions and baseline predictions.

As such, it is easy to interpret the results, because we see a clear path from the baseline value all the way to the actual model’s output.

Input X gradient

Input X gradient, unlike integrated gradients, requires only a single pass through the model, making it one of the fastest attribution method. 

It basically multiplies each input value by the gradient of the model [2]. This is like a first-order Taylor approximation of how a model’s output would change if the inputs were set to 0. 

However, this method presents some important drawbacks. First, it does not have the additive property, making it harder to interpret the results. Second, it can bias the results when used with the ReLU function, because the gradient may return 0, but a feature may still carry some information.

It’s also best to avoid using with LSTM and GRU models, because the tanh and sigmoid activation functions can return very low gradient values, but the features can have a significant impact.

Perturbation-based methods

While there are many methods based on feature perturbation, only Shapley value sampling is available for now in neuralforecast.

This method approximates Shapley values using Monte Carlo sampling of feature permutations [3]. In other words, we randomly sample different orderings of features, and compute how much the predictions change when the feature is included versus when it’s excluded.

This method has the advantages of actually returning Shapley values (based on the cooperative game theory) and it has the additive property.

However, it is extremely slow. It typically requires hundreds to thousands of model evaluations to obtain a reasonable approximation of Shapley values.

When to use each explainer

As we can see, there are different explainers that function slightly differently and will return different feature attributions. As such, it is important to consider your needs and use case to choose the right explainer.

In the vast majority of situations, using integrated gradients is a good default choice, since it is a fast method with the additive property, making it easy to understand how to go from the baseline value to the actual forecast.

In cases where you want to consider feature interactions in your explanations, then Shapley value sampling should be used. 

Finally, if you really only care about the importance of features relative to one another, or you have a fairly simple deep learning model (a basic MLP), then input X gradient is an okay choice.

Now that we understand the different categories of explainers and having explored each method in more detail, let’s see how we can interpret our forecasts using neuralforecast and the integrated gradients method.

Tutorial: Interpretability with integrated gradients and NHITS

In this hands-on tutorial, we walk through an example to interpret a deep learning model that forecasts the daily visits to my blog.

Here, we use the NHITS model to generate forecasts, because it supports exogenous features, allowing us to measure the impact of external regressors on our target variable.

Before diving into the code, I want to clarify a few aspects about the current implementation in neuralforecast:

• The functionality is in beta, and it is subject to change in the future.
• We currently only support univariate models. See here to see which models are univariate (channel independence) or multivariate (channel dependence).
• While the explainers can attribute all three types of exogenous features (static, historical and future), not all models support all types of exogenous features. Again, see here to see which features are supported for each model. 

As always, the entire source code for this section is fully reproducible and available on GitHub.

Initial setup

First, we import the required packages for this experiment.

import datetime
import numpy as np
import pandas as pd
import shap
import torch
import matplotlib.pyplot as plt

from neuralforecast.core import NeuralForecast
from neuralforecast.models import NHITS
from neuralforecast.losses.pytorch import MAE, MSE

Notice that we import shap, but we only use it for its plotting capabilities. The actual computation of feature attributions is done with captum, which is a library for model interpretability for PyTorch models. 

It is very important to install captum in order to access the explainability functionality of neuralforecast.

Once that this is done, we can load the data.

DATA_URL = "https://raw.githubusercontent.com/marcopeix/time-series-analysis/refs/heads/master/data/medium_views_published_holidays.csv"

df = pd.read_csv(DATA_URL, parse_dates=["ds"])

This dataset compiles the daily number of visitors to my blog from January 1st 2020 to October 12th 2023. I compiled this data myself, and I also included a flag to indicate a US holiday and when a new article is published.

From the figure above, we can see that when a new article is published (marked by a red dot), the number of visitors tends to peak on the same day or a day later. This makes sense because new content tends to attract more readers.

We also notice a weekly seasonality in the data, where the number of visitors is lower during the weekend and higher during the week.

With all of that in mind, let’s train a model to get forecasts and interpret them.

Forecasting with NHITS

To reproduce the following results, we first set the seeds.

np.random.seed(42);
torch.manual_seed(42);

Then, we split the data into a training and test set. We keep the last 28 time steps for the test period, and the rest is for training. 

In the last 28 steps, we have a date where a new article was published, and another instance where it was a holiday, so we will be able to see how those features impact the final model forecast on those particular days.

test = df.tail(28)
train = df.drop(test.index)

Then, we can initialize a NHITS model and train it on our data. Note that both features in this dataset are future exogenous features, because we know their future values with certainty. Ultimately, I decide when to publish an article, and we know in advance the future dates of holidays.

models = [
    NHITS(
        h=28,
        input_size=56,
        futr_exog_list=["published", "is_holiday"],
        max_steps=1000,
        loss=MSE(),
        valid_loss=MAE(),
        early_stop_patience_steps=3,
        scaler_type="robust",
        accelerator="cpu",
    ),
]

nf = NeuralForecast(
    models=models, 
    freq="D", 
)

nf.fit(df=train, val_size=28)

In the code block above, I set accelerator="cpu", because it’s a little quirk for using captum on a Mac. If you are on a Windows or Linux machine, you can omit that parameter to leverage the GPU if it’s available.

Once the model is trained, we can then make predictions and compute feature attributions using the explain method of NeuralForecast.

futr_df = test.drop(columns=["y"])

fcsts_df, explanations = nf.explain(
    futr_df=futr_df, 
    explainer="IntegratedGradients"
)

In the code block above, notice that we pass the future values of our exogenous features over the forecast horizon in futr_df. Again, this is reasonable, because we safely know in advance when an article will be published, and when the next holiday will be.

Also, we specify explainer="IntegratedGradients". We could also specify ShapleyValueSampling or InputXGradient.

And just like that, we have stored the baseline predictions and feature attributions in the explanations variable.

We can also plot the predictions against the actual values in the image below.

While the model captured the seasonality correctly, it failed to predict the big peak in visitors after the article was published.

Now, let’s use the information stored in explanations to interpret our model’s forecast.

Interpret the model’s forecast

In explanations, we stored a dictionary with the attribution score for each type of exogenous features and past values of the series. We also store the baseline predictions.

The first key in the dictionary is the model’s name. In this case, since we only trained a single model, we access the attribution scores using explanations["NHITS"].

The following list specifies what tensors are stored and their shape:

• insample: [batch_size, horizon, n_series, n_output, input_size, 2 (y attribution, mask attribution)]
• futr_exog: [batch_size, horizon, n_series, n_output, input_size+horizon, n_futr_features]
• hist_exog: [batch_size, horizon, n_series, n_output, input_size, n_hist_features]
• stat_exog: [batch_size, horizon, n_series, n_output, n_static_features]
• baseline_predictions: [batch_size, horizon, n_series, n_output]

Let’s break down each element.

Here, insample designates the past values of the series. Then, we have each type of exogenous feature supported in neuralforecast. Finally, we have the baseline predictions. This is important for methods with the additive property, because the sum of the baseline with all the other elements gives the final forecast.

Then, the batch_size designates how many series we are forecasting. In this case, the value is 1.

The horizon is the forecast horizon, so its value is 28.

Next, n_series represents how many series the model is using to forecast. This is related to univariate or multivariate models. When it’s univariate, the value is 1 (as is the case here). Note that for now, only univariate models can compute feature attributions in neuralforecast.

We then have n_output which is related to point forecasts or probabilistic forecasts. Here, we only made point forecasts, so the value is 1. If we trained a probabilistic model to output prediction intervals, then we could have 3 outputs:

1. median forecast
2. upper bound
3. lower bound

Finally, for all types of exogenous features, the last dimension is the number of features. For insample, the last dimension is 2 because it contains the attribution of the past values and the mask. The mask is simply a flag to say if a value is missing or not. We usually sum both attribution scores to get rid of that granularity.

So, in this current scenario, we can access these elements:

• explanations["NHITS"]["insample"] — shape of [1, 28, 1, 1, 56, 2]
• explanations["NHITS"]["futr_exog"] — shape of [1, 28, 1, 1, 84, 2]
• explanations["NHITS"]["baseline_predictions"] — shape of [1, 28, 1, 1]

So, let’s now create a waterfall plot to see how different features and past values affected the model’s final forecast. Here, let’s focus on the date where a new article was published.

holiday_idx = 24
published_idx = 25

batch_idx = 0
horizon_idx = published_idx
output_idx = 0

attributions = []
feature_names = []

# Insample attributions
y_attr = explanations["NHITS"]["insample"][batch_idx, horizon_idx, 0, output_idx, :, 0]
mask_attr = explanations["NHITS"]["insample"][batch_idx, horizon_idx, 0, output_idx, :, 1]
combined_insample = (y_attr + mask_attr).cpu().numpy()
for i, attr in enumerate(combined_insample):
    attributions.append(attr)
    feature_names.append(f"y_lag{i+1}")

# futr_exog attributions
futr_attr = explanations["NHITS"]["futr_exog"][batch_idx, horizon_idx, 0, output_idx]
futr_attr = futr_attr.cpu().numpy()

feature_labels = ["published", "is_holiday"]

for feat_idx in range(futr_attr.shape[1]):
    for t in range(futr_attr.shape[0]):
        attributions.append(futr_attr[t, feat_idx])
        if t < 56:  # Historical values (input_size part)
            feature_names.append(f"{feature_labels[feat_idx]}_hist_t{t+1}")
        else:  # Future values (horizon part)
            feature_names.append(f"{feature_labels[feat_idx]}_futr_h{t-55}")

shap_values = np.array(attributions)

baseline = float(explanations["NHITS"]["baseline_predictions"][batch_idx, horizon_idx, 0, output_idx].cpu())

# Create SHAP Explanation
shap_explanation = shap.Explanation(
    values=shap_values,
    base_values=baseline,
    feature_names=feature_names
)

shap.plots.waterfall(shap_explanation)

The code block above gives the plot below.

From the figure above, we can see that the top feature with the highest contribution is our published feature, which pushed the forecast upward by 610. This is great to see, because it aligns with our expectations that a new article does indeed bring more visitors to the blog.

Also, from the plot above, we can see that the explainer shows us that the baseline value is 1276, but then after considering the impact of past values and exogenous features, it eventually reached the value of 1871.

This is clearer if we aggregate all time steps together as shown below.

batch_idx = 0
horizon_idx = published_idx
output_idx = 0

baseline = float(explanations["NHITS"]["baseline_predictions"][batch_idx, horizon_idx, output_idx, output_idx].cpu())
insample_sum = float(explanations["NHITS"]["insample"][batch_idx, horizon_idx, output_idx, output_idx, :, :].sum().cpu())

# Separate sums for each future exogenous feature
published_sum = float(explanations["NHITS"]["futr_exog"][batch_idx, horizon_idx, output_idx, output_idx, :, 0].sum().cpu())
is_holiday_sum = float(explanations["NHITS"]["futr_exog"][batch_idx, horizon_idx, output_idx, output_idx, :, 1].sum().cpu())

feature_names = []
shap_values = []

if insample_sum != 0:
    feature_names.append("Historical Y (all lags)")
    shap_values.append(insample_sum)

if published_sum != 0:
    feature_names.append("Published (futr_exog)")
    shap_values.append(published_sum)

if is_holiday_sum != 0:
    feature_names.append("Is Holiday (futr_exog)")
    shap_values.append(is_holiday_sum)

shap_values = np.array(shap_values)

# Create SHAP Explanation
shap_explanation = shap.Explanation(
    values=shap_values,
    base_values=baseline,
    feature_names=feature_names
)

shap.plots.waterfall(shap_explanation)

By summing across all time steps, we get a broader view of the impact of each feature. Here, the past values of the series still has the greatest impact. Then, publishing an article pushed the forecast further up. 

Now, let’s take a look at the date where a holiday occured to see if that feature has an impact on the final forecast.

batch_idx = 0
horizon_idx = holiday_idx
output_idx = 0

attributions = []
feature_names = []

# Insample attributions
y_attr = explanations["NHITS"]["insample"][batch_idx, horizon_idx, 0, output_idx, :, 0]
mask_attr = explanations["NHITS"]["insample"][batch_idx, horizon_idx, 0, output_idx, :, 1]
combined_insample = (y_attr + mask_attr).cpu().numpy()
for i, attr in enumerate(combined_insample):
    attributions.append(attr)
    feature_names.append(f"y_lag{i+1}")

# futr_exog attributions
futr_attr = explanations["NHITS"]["futr_exog"][batch_idx, horizon_idx, 0, output_idx]
futr_attr = futr_attr.cpu().numpy()

feature_labels = ["published", "is_holiday"]

for feat_idx in range(futr_attr.shape[1]):
    for t in range(futr_attr.shape[0]):
        attributions.append(futr_attr[t, feat_idx])
        if t < 56:  # Historical values (input_size part)
            feature_names.append(f"{feature_labels[feat_idx]}_hist_t{t+1}")
        else:  # Future values (horizon part)
            feature_names.append(f"{feature_labels[feat_idx]}_futr_h{t-55}")

shap_values = np.array(attributions)

baseline = float(explanations["NHITS"]["baseline_predictions"][batch_idx, horizon_idx, 0, output_idx].cpu())

# Create SHAP Explanation
shap_explanation = shap.Explanation(
    values=shap_values,
    base_values=baseline,
    feature_names=feature_names
)

shap.plots.waterfall(shap_explanation)

From the plot above, we see the breakdown of the impact of each feature for the forecast on October 9th 2023, which was Columbus Day, a national holiday in the US. 

Interestingly, it seems like knowing it was a holiday did not impact the model. In fact, past values and knowing if an article is published or not remain the most important.

Thus, it seems that holidays may not have that big of an impact as we might think.

Studying the impact of past values

Before concluding this article, let’s plot the impact of only past values on the forecast. Here, we plot the top five lags for the last 12 forecast steps.

data = explanations["NHITS"]["insample"][0, :-12, 0, 0, :, 0]  # Shape: [12, 56]

# Create 4x3 subplot grid
fig, axes = plt.subplots(4, 3, figsize=(18, 16))
axes = axes.flatten()

for h in range(12):
    horizon_attrs = data[h, :].cpu().numpy()
    
    # Get top 5 lags by absolute value
    top5_indices = np.argsort(np.abs(horizon_attrs))[-5:]
    top5_values = horizon_attrs[top5_indices]
    top5_lags = top5_indices + 1
    
    # Sort by absolute value (most important at top)
    sort_order = np.argsort(np.abs(top5_values))
    
    # Create horizontal bar plot
    colors = ['red' if v < 0 else 'blue' for v in top5_values[sort_order]]
    axes[h].barh(range(5), top5_values[sort_order], color=colors, alpha=0.7)
    axes[h].set_yticks(range(5))
    axes[h].set_yticklabels([f'Lag {lag}' for lag in top5_lags[sort_order]])
    axes[h].set_title(f'Horizon Step {h+17}', fontsize=12, pad=10)
    axes[h].axvline(x=0, color='black', linestyle='-', alpha=0.5, linewidth=0.8)
    axes[h].grid(axis='x', alpha=0.3)
    axes[h].set_xlabel('Attribution Score', fontsize=10)

plt.suptitle('Top 5 Most Important Lags for Each Horizon Step (First Series)\nNHITS Model - IntegratedGradients', 
             fontsize=16, y=0.98)

# Add legend
from matplotlib.patches import Patch
fig.legend(handles=[Patch(facecolor='blue', alpha=0.7, label='Positive'),
                    Patch(facecolor='red', alpha=0.7, label='Negative')],
          loc='upper center', bbox_to_anchor=(0.5, 0.02), ncol=2)

plt.tight_layout()
plt.subplots_adjust(top=0.92, bottom=0.08)
plt.show()

The main thing I want to point out in the plot above is that lag 56 (the most recent time step in the input sequence) appears on all plots for each forecast step.

This is to be expected, especially in series with a strong trend. Generally, forecasting models will give more weight to the most recent history to forecast the next set of values, and the plot above confirms this expectation.

Conclusion

In this article, we have explored different methods to explain deep learning forecasting models like integrated gradients and Shapley value sampling.

We distinguished between gradient-based and perturbation-based methods and learned that gradient-based methods are often preferred due to their speed and accuracy in estimating features attributions.

We then applied that in a scenario where we want to measure the impact of exogenous features in predicitng daily visits to a blog. Using neuralforecast’s intuitive interface, we managed to observe that knowing when an article is published has a large impact on the final forecast, resulting in the model increasing its final prediction.

While the functionality is still in beta and actively being developed in neuralforecast, I believe this is an important step in making our forecasts explainable. Deep learning models do not have to be black box models anymore.

Thanks for reading! I hope that you enjoyed it and that you learned something new!

Learn the latest time series analysis techniques with my free time series cheat sheet in Python! Get the implementation of statistical and deep learning techniques, all in Python and TensorFlow!

Cheers 🍻

References

[1] M. Sundararajan, A. Taly, and Q. Yan, “Axiomatic Attribution for Deep Networks.” Available: https://arxiv.org/pdf/1703.01365

[2] A. Shrikumar, P. Greenside, A. Shcherbina, and A. Kundaje, “Not Just a Black Box: Learning Important Features Through Propagating Activation Differences,” arXiv:1605.01713 [cs], Apr. 2017, Available: https://arxiv.org/abs/1605.01713

[3] J. Castro, D. Gómez, and J. Tejada, “Polynomial calculation of the Shapley value based on sampling,” Computers & Operations Research, vol. 36, no. 5, pp. 1726–1730, May 2009, doi: https://doi.org/10.1016/j.cor.2008.04.004.