Various approaches to explaining predictions of black box models have been proposed, including model-agnostic techniques that measure feature importance (or e ect) by presenting modi ed test instances to the underlying black-box model. These modi cations rely on choosing feature values from the complete range of observed values. However, when applying machine learning algorithms to the task of forecasting from multivariate time-series, it is suggested that the temporal aspect should be taken into account when analyzing the feature e ect. A modi cation of individual conditional expectation (ICE) plots is proposed, called ICE-T plots, which displays the prediction change for temporally ordered feature values. Results are presented from a case study on predicting the Swedish gross domestic product (GDP) based on a comprehensive set of indicator and prognostic variables. The e ect of calculating feature e ect with and without temporal constraints is demonstrated, as well as the impact of transformations and forecast horizons on what features are found to have a large e ect, and the use of ICE-T plots as a complement to ICE plots.
Machine learning for time series analysis has received signi cant attention over
the years, in particular classi cation of (univariate or multivariate) series, see
e.g., [
In a recent study [
In this study, we will consider the task of forecasting the GDP of Sweden,
a country not included in IMF's original study, and investigate the application
of techniques for explaining predictions, by analyzing the impact each feature
has on predictions of the underlying (black-box) model. In the standard
machine learning setting, where a model is trained from a set of examples that is
assumed to be independently sampled according to some xed but unknown
distribution, the e ect (or importance) of a feature with respect to the model, may
be estimated by measuring how the predictive performance (or the predictions)
change when modifying the values of the features for a set of examples. Such
a procedure was proposed in [
In the next section, we describe the various approaches to explaining predictions by feature e ects in more detail, including the novel ICE-T plot. In Section 3, we describe the considered prediction task and dataset together with the experimental setup and present ndings from the empirical investigation. Finally, in Section 5, we summarize the main conclusions and outline directions for future research. 2
In the following, we will assume that we have an ordered set X 2 RN D of N objects with D features, and an ordered set of labels Y 2 RN . Let Xi;j denote the element at row i and column j in the matrix X, and Yi denote the ith element of the vector Y . Let Xi;: = (Xi;1; : : : ; Xi;D) denote the ith row (object) of X, Xa:b denote the sequence of rows (Xa;:; : : : ; Xb;:), X:;i = (X1;i; : : : ; XN;i) denote the ith column of X, and Ya:b denote the sequence of elements (Ya; : : : ; Yb). Let u(x = (x1; : : : ; xD); i; v) = x0 = (x1; : : : ; xi 1; v; xi+1; : : : ; xD), i.e., given an object x 2 RD, the function returns an updated object x0 2 RD where the ith feature value xi has been replaced by the value v.
We will moreover assume that we have an underlying (black-box) model M , such that given an object x 2 RD, it returns a (predicted) label y^ = M (x) 2 R.
Given an object x 2 RD, a set of values V and a model M , the feature e ect F E on feature i, is de ned as:
F E(x; M; i; V ) = M (x)
Pv2V M (u(x; i; v)) jV j (1)
The above function hence calculates the di erence between the original prediction for the object x and the average prediction from updating the object on feature i with values from V .
In case V = X:;i for some random sample X drawn independently from the same distribution as the training set that was used to construct M , then F E provides an estimate of the expected prediction change relative to this distribution (and feature). However, as we are here primarily interested in understanding what e ect di erent features have on a speci c prediction, rather than providing an unbiased estimate of the expected change for unseen examples, we will, as commonly done, allow the feature e ect to be estimated with respect to the training (or any other) set. It should be noted that the feature e ect is in contrast to variable (or feature) importance, not de ned in relation to the labels (Y ), but only considers the average change of the prediction.
When producing PDP [
ICE(x; M; i; X) = f(v; F E(x; M; i; fvg)) : v 2 X:;ig (2)
When plotted, with the feature values on the x-axis and the feature e ect on the y-axis, the points are normally connected by lines, e ectively providing an interpolation of the predictions of the underlying model between each pair of consecutive feature values.
For multivariate time series data, we typically have timestamps for the objects in X, i.e., T = (t1; : : : ; tN ). Assuming these time points to be unique, we de ne an ICE-T plot by the following set of points:
ICE-T (x; M; i; X; T ) = f(th; F E(x; M; i; fXh;ig)) : h = 1; : : : ; N g (3)
In contrast to the ICE plot, the ICE-T plot hence allows for visualizing the feature e ect over time, where the x-axis represents the time at which feature values have been observed, rather than speci c feature values. Since the actual feature values are not included in an ICE-T plot, it should be considered to give a complementary view to the ICE plot. 3
In this section, we rst describe the prediction task that is considered in this study. We then describe how the empirical investigation has been designed, including the data preparation, before presenting the ndings from the empirical investigation. 3.1
GDP, or gross domestic product, is a measure of the total economic activity
taking place on an economic territory which leads to output meeting the nal
demands of the economy. GDP is an aggregate in the national accounts, an
accounting system meant to summarise and describe the country's economic
activities and development. There are in principle three ways to compute GDP:
1) the production approach, which is the sum of all value added from produced
goods and services, 2) the expenditure approach, which is the sum of all
expenditures made for consuming the output of the economy or adding to wealth, or
3) the income approach, which is the sum of all incomes earned by producing
goods and services [
Swedish GDP is compiled and published by Statistics Sweden (SCB). GDP
and other national accounts are made in accordance with a common European
standard [
In addition to parts and aggregates of the national accounts, there are also other economic features that should have predictive power for estimates of GDP. For this study, we have chosen the most prominent features in the economic outlook in forecasts from the Swedish National Financial Management Authority. Examples of such features are unemployment gures, in ation and interest rates. Apart from Swedish domestic economic features, there are also exchange rates and foreign interest rates as those are relevant for the export-oriented economy of Sweden.
There are also a small number of economic indicators used as features. These can be backward- or forward-looking (GDP in itself is often used as a backwardlooking indicator of the general health of the economy). Backward-looking indicators that are used in this model are for example reported vacancies from employers, redundancy notices and number of newly purchased cars. Indicators like these are often used as they convey trends in actual, not calculated, economic activity. They are used both in forecasting and in business reporting in the general news.
Forward-looking features are in general used to capture intent and predictions
for the future through surveys. Answers from such surveys are then summarised
to an index which can be tracked over time. An example of such a study
incorporated in this model is the Swedish Economic Tendency Survey conducted by the
Swedish National Institute for Economic Research. The results from the survey
is used to construct an indicator. The Economic Tendency Indicator is based on
monthly surveys of households and rms and consequently captures the
sentiment among these players in the Swedish economy. The indicator is based on the
information contained in the con dence indicators for manufacturing, services,
construction, retail and consumers [
All features used in the model have values on a quarterly basis. The values for some of the features are available on a monthly or even daily basis, but they are either summed or averaged to a quarterly value in the model. Data is available for most features from 1993 up to the last quarter of 2018, and features with missing values in this time span have been excluded. The complete multivariate time-series considered in this study contains 104 objects with 68 features, in addition to time points and outcome (GDP). 3.2
The experiment has been designed to emulate a realistic scenario, where only
data available up to the point in time in which a prediction is made, may be used
for generating (and explaining predictions of) a model to predict the outcome
(GDP) at a speci ed later time point. We here consider two time frames;
predicting the GDP one quarter of a year ahead and four quarters of a year ahead,
respectively. Moreover, in addition to considering the actual GDP as the target
variable, we will consider predicting the percentage change relative to the time
point at which the prediction is made, and investigate how this transformation
a ects both predictive performance and feature e ects. In addition to using the
features described in the previous section, we will also, as is common in
timeseries forecasting, use the (recent) previous outcomes as features; in this study
we will use ve features to represent the outcome of the current (the time point
at which the prediction is made) and four preceding quarters, which we refer to
as lagged values. It should be noted that this means that some of the lagged
values will be missing for the rst four objects in the time series (as we do not know
what the earlier outcome was), and rather than handling missing values, e.g., by
imputation, we simply exclude these objects from the training set. Furthermore,
if we at the current time point th want to make a prediction for time point th+a,
i.e., a quarters ahead, we cannot assume that the outcome for the time points
th+1; : : : ; th+a 1 are known at the time of prediction. Consequently, we have
also excluded the corresponding a 1 objects, at time points th a+1; : : : ; th 1,
from the training set. Note that the latter only a ects the second scenario in
which we are making predictions four quarters ahead, i.e., resulting in that three
preceding objects are excluded from the training set. Rather than searching for
the optimal predictive model with careful hyperparameter tuning, we have here
opted to use the standard GradientBoostingRegressor as implemented in [
Given the complete multi-variate time series X 2 RN D, time points T = (t1; : : : ; tN ) and outcomes O 2 RN , a model to make a prediction at the time point th 2 Tb:N a for the outcome at time point th+a is generated from the objects Xb:h a and labels Y = Ob+a:h, where b is the number of lagged values and a is the number of time steps ahead for which a prediction is made. As stated above, we will consider b = 5 together with a = 1 and a = 4, respectively. Moreover, to handle the so-called cold-start problem, we will only make predictions for the last 64 time points for which the outcome a quarters ahead are known.
We will investigate the feature e ect for a subset of the features over time. We have considered two options for aggregating the feature e ect at time-point th; averaging the (absolute) feature e ect using feature values in the entire training set, i.e., Xb:h a, or averaging using feature values in a time window of size w, i.e., Xh a w+1:h a. In this study, we will consider w = 12, i.e., values from the 12 most recent objects in the training set are used. 4
In Fig. 1 and Fig. 2, the predicted vs. actual outcomes in million Swedish Krona (MSEK), are plotted when making predictions one and four quarters ahead, respectively, and when using the original (blue dashed lines) and transformed targets (orange dashed lines), where the latter concerns the percentage change, which here is projected back to the original scale to allow for a direct comparison. When predicting one quarter ahead, the model using transformed targets clearly outperforms the model using non-transformed targets, with the former obtaining a root mean-squared error (RMSE) of 20448 MSEK and a Pearson correlation coe cient of 0.981, while the latter has an RMSE of 37942 MSEK and a correlation coe cient of 0.936. The performance di erence between the two methods is less clear when predicting four quarters ahead; using original targets leads to an RMSE of 36353 MSEK and a correlation coe cient of 0.954, while for the transformed targets the RMSE is 32918 MSEK with a correlation coe cient of 0.951. Interestingly, the task of predicting one year ahead is slightly more easy than predicting one quarter ahead when using the original targets. In terms of RMSE, the use of the transformation (percentage change) is clearly e ective independently of the forecast horizon, while the correlation coe cient is hardly a ected by this transformation when predicting a year ahead.
In order to analyze the di erences between the models using the original and the transformed targets, we can take a look at the aggregated (absolute) feature e ect, averaged over all predictions and for each prediction calculating the feature e ect with respect to all previously observed feature values. In Fig. 3 and Fig. 4, the average absolute e ect is plotted for the top 20 features (listed in descending order) when making predictions one quarter ahead, without and with the transformation, respectively. Note that although the scales di er, since the former model predicts the actual GDP (in MSEK) while the latter predicts the percentage change, we can still compare them by their relative impact. For the former, one of the lagged variables (nbnpmpf-4), which represents the outcome four quarters before the time of prediction, is dominating, while the e ect is distributed di erently for the latter, although four of the lagged variables appear among the top ve.
Fig. 5 shows that when not having applied the transformation, the picture for predicting four quarters ahead is similar to when making predictions one quarter ahead; two of the lagged variables appear among the top ve, although two di erent ones. However, when having applied the transformation, the picture changes quite drastically when predicting four quarters ahead, as shown by Fig. 6. Here, the highest ranked lagged variable (nbnpmpf-4) is ranked behind eight other features, which indicates that the importance of using lagged variables decreases when having applied the (percentage change) transformation and considering a more distant forecast horizon. However, the graphs displaying aggregated feature e ect do not show how it varies over time and can hence not be related to individual predictions. In Fig. 7, the feature e ect for the nine top ranked features (according to Fig. 6), is displayed over time, where the feature e ect is, as above, calculated with respect to all feature values that have been observed prior to each prediction. One may clearly see that the relative sizes of the e ects vary over time. When calculating the feature e ect in relation to all previously observed feature values, we do not know whether the impact is due to deviations to values observed a long time ago or more recently. One way of measuring feature e ect in relation to more recent observations is to include only the latest feature values in the calculation (limiting the set V in Eq. 1). In Fig. 8, the set of feature values to include is limited to the 12 most recent observations. When focusing the analysis to the most recent feature values, one may observe some dramatic changes to using all observed feature values, e.g., the e ect of the Swedish Krona to US Dollar exchange rate (sek_usd) starts to increase signi cantly before 2010, when using the time constraint, while the impact of nexf (Swedish export) during 2014 disappears.
When focusing on a speci c prediction, we may take a look at an ICE plot (Eq. 2). In Fig. 9 and Fig. 10, ICE plots are shown for the predictions at January 1, 2010 and July 1, 2014, respectively, using the above features for which the values have been min-max-normalized to allow for displaying multiple features in one plot. The graphs show for both predictions that they are higher than the ones obtained with lower values of sek_usd and higher values of eurefi (Euro re nancing rate) and sv3mssvx (three-month rate in Sweden). The second prediction is lower than what is output by the model for all but high values values of nexf (Swedish export), while it is higher than what would be output for high values for the exchange rate of the Swedish Krona to Japanese yen (sek_jpy).
In an ICE plot, we can see how a prediction would be a ected by replacing some speci c feature value with all possible (previously observed) values. However, some feature values may be extreme and occur very infrequently. In addition, the values may also not have appeared for a long time, and hence may be of little relevance when reasoning about the current prediction. To allow for reasoning about the feature e ect in a temporal context, we also take a look at the proposed ICE-T plots (Eq. 3), to study how the feature e ect varies when H. Bostrom, P. Hoglund, S-O. Junker, A-S Oberg and M. Sparr selecting values in the order they have been observed. As a complement to the ICE plots, we present the corresponding ICE-T plots for the predictions at January 1, 2010 and July 1, 2014 in Fig. 11 and Fig. 12, respectively. For the rst plot, we can see that the feature sv3mssvx (three-month rate in Sweden) only has an impact if considering feature values that occurred more than a decade before the time of prediction. Moreover, the e ect of sek_usd is (very) high, only if considering relatively recent values, and this feature has hardly any e ect if considering values observed earlier. In the second plot, one may observe that the prediction is lower than what is obtained when using values for nexf (Swedish export) that were observed more than 12 quarters earlier, which explains why the e ect of this feature was not visible for the same date when calculating feature e ect using the time constraint, as shown in Fig. 8. We have investigated ways to explain multivariate time-series forecasting models, by measuring and presenting the feature e ect. In addition to calculating the aggregated feature e ect, measured with or without a time constraint, we have presented an approach to visualize the feature e ect on individual predictions, called the ICE-T plot, which complements the ICE plot by showing the feature e ect over time. We have presented an application of these techniques to the task of predicting the Swedish GDP. In addition to demonstrating the use of the techniques, the empirical investigation has also highlighted the impact of target variable transformation and forecast horizon on the feature e ect.
The work has focused on explaining predictions of a black-box model by analyzing the feature e ect, i.e., how the output of the model changes when changing the input. For retrospective analysis, all the techniques considered in this work can be straightforwardly adapted to instead calculate feature importance, i.e., how the predictive performance, such as absolute error, is a ected by changing the input. However, in a real prediction scenario, one does not have access to the correct target (for the prediction at hand), and hence the latter option is not available.
A natural extension of the current work is to consider the e ect of changing multiple features simultaneously. In contrast to PDP and ICE plots, which need one axis per feature included in the combination, the ICE-T plots would remain two-dimensional, since the time points, rather than feature values, are used to align the predictions.
In the current study, we have only considered one speci c type of blackbox model (generated by gradient boosting), and a direction for future work is to study how the choice of black-box model (including various hyperparameter settings) a ects what features are considered to have an impact. Finally, there are several alternatives to using feature e ect to explain predictions and the application and adaptation of these techniques to the speci c requirements that multivariate time-series forecasting models provide could be a fruitful area of research.
The study was funded by Swedish Governmental Agency for Innovation Systems (grant no. 2019-02252). HB was partly funded also by the Swedish Foundation for Strategic Research (grant no. BD15-0006).