Temporal Fusion Transformers for Interpretable Multi-horizon Time Series Forecasting

Multi-horizon forecasting problems often contain a complex mix of inputs -- including static (i.e. time-invariant) covariates, known future inputs, and other exogenous time series that are only observed historically -- without any prior information on how they interact with the target. While several deep learning models have been proposed for multi-step prediction, they typically comprise black-box models which do not account for the full range of inputs present in common scenarios. In this paper, we introduce the Temporal Fusion Transformer (TFT) -- a novel attention-based architecture which combines high-performance multi-horizon forecasting with interpretable insights into temporal dynamics. To learn temporal relationships at different scales, the TFT utilizes recurrent layers for local processing and interpretable self-attention layers for learning long-term dependencies. The TFT also uses specialized components for the judicious selection of relevant features and a series of gating layers to suppress unnecessary components, enabling high performance in a wide range of regimes. On a variety of real-world datasets, we demonstrate significant performance improvements over existing benchmarks, and showcase three practical interpretability use-cases of TFT.


INTRODUCTION
Multi-horizon forecasting, i.e the prediction of variables-of-interest at multiple future time steps, is a crucial aspect of machine learning for time series data. In contrast to one-step-ahead predictions, multi-horizon forecasting provides decision makers access to estimates across the entire path, allowing them to optimize their course of action at multiple steps in future. One common aspect of major forecasting scenarios is the availability of different data sourcesincluding known information about the future (e.g. holiday dates), other exogenous time series, and static metadata -without any prior knowledge on how they interact. As such, the identification of key drivers of predictions can be important for decision makers, providing additional insights into temporal dynamics. For instance, static (i.e. time-invariant) covariates often play a key role -such as * Completed as part of internship with Google Cloud AI Research in healthcare where genetic information can determine the expression of a disease [34]. Given the numerous real-world applications of multi-horizon forecasting, e.g. in retail [5,10], medicine [26,42] and economics [6], improvements in existing methods bear much significance for practitioners in many domains.
Deep neural networks have increasingly been used in multihorizon time series forecasting, demonstrating strong performance improvements over traditional time-series models [1,28,31]. While many architectures have focused on recurrent neural network designs [15,31,39], recent improvements have considered the use of attention-based methods to enhance the selection of relevant time-steps in the past [13] -including Transformer-based models in [25]. However, these methods often fail to consider all different types of inputs commonly present in multi-horizon prediction problems, either assuming that all exogenous inputs are known into the future [15,25,31]-a common problem with autoregressive models -or neglecting important static covariates [39] -which are simply concatenated with other time-dependent features at each step. With many improvements in time-series models resulting from the alignment of architectures with unique data characteristics [23,30], similar performance gains could also be reaped by designing networks with suitable inductive biases for multi-horizon forecasting.
Most multi-horizon prediction architectures are 'black-box' models, where forecasts are controlled by complex nonlinear interactions between many parameters, that render explainability challenging. In turn, poor interpretability can make it difficult for model builders to improve the model quality, for business decision makers to trust a model's outputs or for customers in understanding the outcomes of a product -due to the lack of insights into what is driving its forecast. Moreover, commonly-used methods for interpretability in deep neural networks can be further limited in time series settings. Conventional post-hoc explainability methods (e.g. LIME [32] and SHAP [27]), for example, typically do not consider the time ordering of input features -with surrogate models independently constructed for each data-point or with features assumed to be independent of others (including those at neighboring time steps). This can potentially lead to poor quality explanations for time series data, where dependencies between time steps are typically significant. In addition, attention-based models, such as the Transformer architecture [36], only provide insights on the relevant time-steps in their conventional form, but not into important features.
In this paper, we propose the Temporal Fusion Transformer (TFT) -an attention-based architecture which combines high performance multi-horizon forecasting with interpretable insights. For performance improvements over state-of-the-art benchmarks, we introduce several novel adjustments to align the architecture with the full range of potential inputs and temporal relationships common to multi-horizon forecasting -specifically incorporating (1) static covariate encoders which encode context vectors for use in other parts of the network, (2) gating mechanisms throughout and sample-dependent variable selection to minimize the contributions of irrelevant inputs, (3) a sequence-to-sequence layer to locally process known and observed inputs, and 4) a temporal self-attention decoder to learn any long-term dependencies present within the dataset. The use of specialized components also facilitates interpretability, for which we propose three use cases: to identify (i) globally-important variables for the prediction problem, (ii) persistent temporal patterns, and (iii) significant events. On real-world data, we show how these methods can be practically applied and the insights they bring.

RELATED WORKS
Deep Learning Models for Multi-horizon Forecasting. In line with traditional methods for multi-horizon forecasting [29,35], recent deep learning methods can be categorized into iterated approaches using autoregressive models [15,25,31] or direct methods using sequence-to-sequence models [13,39].
Iterated approaches utilize one-step-ahead prediction models, with multi-step predictions obtained by recursively feeding predictions into future inputs. For instance, approaches with Long Short-term Memory (LSTM) [20] networks have been considered -such as Deep AR models [15] which use 3 stacked LSTM layers to generate parameters of one-step-ahead Gaussian predictive distributions. Deep State-Space Models (DSSM) [31] adopt a similar approach, utilizing LSTMs to generate parameters of a predefined linear state-space model with predictive distributions produced via Kalman filtering -with extensions for multivariate time series data in [38]. More recently, Transformer-based architectures have been explored in [25], which propose the use of convolutional layers for local processing, and a sparse attention mechanism to increase the size of the receptive field during forecasting. Despite their simplicity, iterative methods rely on the assumption that the values of all variables excluding the target are known at forecast timesuch that only the target needs to be recursively fed into future inputs. However, in many practical scenarios, numerous useful time-varying inputs exist. As they are unknown in advance, their straightforward use is limited for iterative approaches. TFTs, on the other hand, explicitly account for the diversity of inputs -naturally handling static covariates and known/unknown time-varying inputs.
Direct methods are trained to explicitly generate forecasts for multiple predefined horizons at each time step. Their architectures typically rely on sequence-to-sequence models, using LSTM encoders to summarize historical inputs, and a variety of methods to generate future predictions. The Multi-horizon Quantile Recurrent Forecaster (MQRNN) [39], for example, utilizes LSTM or convolutional encoders to generate context vectors, which are feed into multi-layer perceptrons (MLPs) for each horizon. In [13] a multi-modal attention mechanism is used with LSTM encoders to construct context vectors for a bi-directional LSTM decoder. Despite performance gains over LSTM-based iterative methods, interpretability is still not straightforward with standard direct methods. In contrast, we show 3 use-cases for interpreting attention patterns in TFTs to produce general insights about temporal dynamics, while maintaining state-of-the-art performance on a variety of datasets.
Time Series Interpretability with Attention Weights. Attention mechanisms are used in translation [36], image classification [37] or tabular learning [3] to identify salient portions for a specific example -using the magnitude of attention weights to determine the importance of different locations. Recent papers have also proposed the use of attention-based mechanisms for time-series interpretability [1,7,25], with both LSTM-based [33] and transformer-based [25] models. However, the importance of static covariates -which may be applicable across all time-steps -may be lost with temporal importance, as these methods typically blend variables at each input. TFT alleviates this by using separate encoder-decoder attention for static features at each time step, on top of the self-attention to determine the contribution time-varying inputs.
Instance-wise Variable Importance with Deep Neural Networks. Instance-wise variable importance can be obtained with post-hoc explanation methods [27,32,40] and inherently-intepretable models [7,18]. Post-hoc explanation methods, such as LIME [32], SHAP [27] and RL-LIM [40], are applied on pre-trained black-box models. They are often based on distilling into a surrogate interpretable model, or decomposing into feature attributions. These methods are not designed to take into account the time-wise ordering of inputs -i.e. they ignore sequential dependencies between the input -making it challenging to directly apply them to complex time series data. Inherently-interpretable model design approaches build components for feature selection directly into the architecture itself. For time-series forecasting specifically, they are based on explicitly quantifying time-dependent variable contributions. For example, Interpretable Multi-Variable LSTMs [18] partition the hidden state such that each variable contributes uniquely to its own memory segment, and weights memory segments to determine variable contributions. Methods combining temporal importance and variable selection have also been considered in [7], which computes a single contribution coefficient based on attention weights from each. However, in addition to modelling only one-step-ahead prediction problems, existing methods also focus on instance-specific interpretations of attention weights -without providing insights into general temporal dynamics. In contrast, the use-cases demonstrated in Section 7 demonstrate the capability of TFT in analyzing global temporal relationships to build insights about the data as a whole.

MULTI-HORIZON FORECASTING
The general problem of multi-horizon forecasting is depicted in Fig. 1. Let there be I unique entities in a given time series dataset -such as different stores for retail forecasting or patients in the medical context. Each entity i is associated with a set of static covariates s i ∈ R m s , as well as inputs χ i,t ∈ R m χ and scalar targets y i,t ∈ R at each time-step t ∈ [0,T i ]. In a general sense, which can only be measured at each step and are unknown beforehand, and known inputs x i,t ∈ R m x which can be predetermined (e.g. the day-of-week at time t).
In many scenarios, the provision for prediction intervals can be useful for optimizing business decisions and risk management, by giving decision makers an indication of likely best and worstcase values that the target can take. As such, we adopt quantile regression to our multi-horizon forecasting setting (e.g. outputting the 10 t h , 50 t h and 90 t h percentiles at each time step). Each quantile forecast takes the form: whereŷ i,t +τ (q, t, τ ) is the predicted q t h sample quantile of the τstep-ahead forecast at time t, and f q (.) is a prediction model. In line with other direct methods, our model simultaneously outputs forecasts for τ max discrete prediction horizons -i.e. τ ∈ {1, . . . , τ max }. We incorporate all past information within a finite look-back window k, using target and known inputs only up till and including the forecast start time t (i.e. y i,t −k :t = y i,t −k , . . . , y i,t ) and known inputs across the entire range (i.e. x i,t −k:t +τ = x i,t −k , . . . , x i,t , . . . , x i,t +τ ). For notational simplicity, we omit the subscript i through the papers unless explicitly required.

MODEL ARCHITECTURE
We design the Temporal Fusion Transformer (TFT) to use canonical components to efficiently build feature representations for each input type (i.e. static, known inputs, observed inputs), enabling it to obtain high forecasting performance on a wide range of problems. The major constituents of the TFT are: (1) Gating Mechanisms -to skip over any unused components of the architecture, providing adaptive depth and network complexity to accommodate a wide range of datasets and scenarios. Gated Linear Units extensively are utilized throughout our architecture, and Gated Residual Network is proposed as a main building block. (2) Variable Selection Networks -to select relevant input variables at each time step. (3) Static Covariate Encoders -to integrate static features into the network, through encoding of context vectors to condition temporal dynamics. (4) Temporal Processing -to learn both long-and short-term temporal relationships, while naturally handling both observed and a priori know time-varying inputs. A sequenceto-sequence layer is employed for local feature processing, whereas long-term dependencies are captured using a novel interpretable multi-head attention block.

Gating Mechanisms
As previously highlighted, the precise relationship between exogenous inputs and targets is often unknown in advance, making it difficult to anticipate which variables are relevant. Moreover, it makes it difficult to determine the extent of non-linear processing required, and there may be instances where simpler models can be beneficial -e.g. when datasets are small or noisy.
To apply non-linear processing only where needed, we introduce the Gated Residual Network (GRN) as a basic building block of TFT, as shown in in Fig. 2. At the lowest level, the GRN takes in a primary input a and an optional context vector c and yields: where ELU is the Exponential Linear Unit activation function [8], is standard layer normalization of [24], and ω is an index used to denote how weights are shared. We adopt component gating layers based on Gated Linear Units (GLUs) [11] to provide the flexibility to suppress any parts of the architecture that are not required for a given dataset. Letting γ ∈ R d mod el be the input, the GLU then takes the form: where σ (.) is the sigmoid activation function,W (.) ∈ R d mod el ×d mod el , b (.) ∈ R d mod el are the weights and biases, ⊙ is the element-wise Hadamard product, and d model is the hidden state size (common across the TFT). GLU allows the TFT to control the extent to which the GRN contributes to the original input a -potentially skipping over the layer entirely if necessary. For instances without a context vector, the GRN simply treats the contex input as zero -i.e. c = 0 in Eq. (4). During training, dropout is applied before the gating layer and layer normalization -i.e. to η 1 in Eq. (3).

Variable Selection Networks
While multiple variables may be available, their relevance and specific contribution to the output are typically unknown. The TFT is designed to provide instance-wise variable selection -through the use of variable selection networks applied to both static covariates and time-dependent covariates. Beyond providing insights into which variables are most significant for the prediction problem, variable selection also allows the TFT to remove any unnecessary noisy inputs which could negatively impact performance. We use entity embeddings [16] for categorical variables and linear transformations for continuous variables, to transform each input variable into a (d model )-dimensional vector -matching the dimensions in subsequent layers for skip connections. In addition, all static, past and future inputs make use of separate variable selection networks as denoted by different colors in Fig. 2. Without loss of generality, we present the variable selection network for historical inputs -noting that those for other inputs take the same form.
Let ξ (j) t ∈ R d mod el denote the transformed input of the j-th vari- being the flattened vector of all historical inputs at time t. Variable selection weights are generated by feeding both Ξ t and an external context vector c s through a GRN, followed by a Softmax layer: where v χ t ∈ R m χ is a vector of variable selection weights, and c s is obtained from a static covariate encoder (see Section 4.3). For static variables, we note that the context vector c s is omitted -given that it already has access to static information. At each time step, an additional layer of non-linear processing is employed by feeding each ξ (j) t through its own GRN: whereξ (j) t is the processed feature vector for variable i. We note that each variable has its own GRN ξ (j) , with weights shared across all time steps t. Processed features are then weighted by their variable selection weights and combined as below: where v (j) χ t is the j-th element of vector v χ t .

Static Covariate Encoders
To build complex representations of static metadata, we use four separate GRN encoders to produce different context vectors. These are then wired into various locations in the temporal fusion decoder (Section 4.5) where static variables play an important role in processing. Specifically, this includes contexts for 1) temporal variable selection (c s ), 2) local processing of temporal features (c c , c h ), and 3) enriching of temporal features with static information (c e ). As an example, taking ζ to be the output of the static variable selection network, contexts for temporal variable selection would be encoded according to c s = GRN c s (ζ ).

Interpretable Multi-Head Attention
To learn long-term relationships across different time steps, TFT employs a self-attention mechanism. In a broad sense, attention mechanisms scale values V ∈ R N ×d V based on relationships between keys K ∈ R N ×d at t n and queries Q ∈ R N ×d at t n as below: where A() is a normalization function. A common choice is scaled dot-product attention [36]: In the canonical form used in the Transformer [36], multi-head attention uses different heads to attend to different representation subspaces, with each head applying the mechanism of Eq. (9): where W Given that different values are used in each head, analyzing attention weights alone would not be indicative of a particular feature's overall importance. As such, we modify multi-head attention to share values in each head, and employ additive aggregation of all heads at the output: where W V ∈ R d mod el ×d V are value weights shared across all heads, and W H ∈ R d at t n ×d mod el is used for final linear mapping. From Eq. (15), we see that each head is able to learn different temporal patterns, while attending to a common set of input features -which can be interpreted as a simple ensemble over attention weights into combined matrixÃ(Q, K) in Eq. (14). Compared to A(Q, K) in Eq. (10), we can see thatÃ(Q, K) yields an increased representation capacity in an efficient way.

Temporal Fusion Decoder
The temporal fusion decoder uses the series of layers described below to learn temporal relationships present in the dataset: 4.5.1 Locality Enhancement with Sequence-to-Sequence Layer. Points of significance in time series data are often identified in relation to its surrounding values -such as anomalies, change-points or cyclical patterns. Leveraging local context, through the construction of features that utilize pattern information on top of point-wise values, can thus lead to performance improvements in attention-based architectures, as also highlighted in [25]. For instance, [25] adopt a single convolutional layer for locality enhancement -extracting local patterns using the same filter across all time. However, this might not be suitable for cases when observed inputs exist, due to the differing number of past and future inputs. As such, we propose the application of a sequence-to-sequence model to naturally handle these differences -feedingξ t −k :t into the encoder and ξ t +1:t +τ max into the decoder. This then generates a set of uniform temporal features which serve as inputs into the temporal fusion decoder itself -denoted by ϕ(t, n) ∈ {ϕ(t, −k), . . . , ϕ(t, τ max )} with n being a position index. For comparability with commonly-used sequence-to-sequence baselines, we consider the use of an LSTM encoder-decoder model -although other models can potentially be adopted as well. This also serves as a replacement for standard positional encoding, providing an appropriate inductive bias for the time ordering of the inputs. Moreover, to allow static metadata to influence local processing, we use the c c , c h context vectors from the static covariate encoders to initialize the cell state and hidden state respectively for the first LSTM in the layer. We also employ a gated skip connection over this layer: where n ∈ [−k, τ max ] is a position index.

Static Enrichment Layer.
As static covariates often have a significant influence on the temporal dynamics (e.g. genetic information on disease risk), we introduce a static enrichment layer that enhances temporal features with static metadata. For a given position index n, static enrichment takes the form: where the weights of GRN ϕ are shared across the entire layer, and c e is a context vector from a static covariate encoder.

Temporal Self-Attention Layer.
Following static enrichment, we next apply self-attention to the temporal features. All staticenriched temporal features are first grouped into a single matrix -i.e. Θ(t) = [θ (t, −k), . . . , θ (t, τ )] T -and interpretable multi-head attention (see Section 4.4) is applied at each forecast time (with N = τ max + k + 1): to yield where m H is the number of heads. Decoder masking [25,36] is applied to the multi-head attention layer to ensure that each temporal dimension can only attend to features preceding it. Besides preserving causal information flow via masking, the selfattention layer allows the TFT to pick up long-range dependencies that may be challenging for RNN-based architectures to learn. Following the self-attention layer, an additional gating layer is also applied to facilitate training: 4.5.4 Position-wise Feed-forward Layer. Lastly, we apply an additional non-linear processing to the outputs of the self-attention layer. Similar to the static enrichment layer, this makes use of a series of GRNs: where the weights of GRN ψ are shared across the entire layer. As per Fig. 2, we also apply a gated residual connection which skips over the entire transformer block, providing a direct path to the sequence-to-sequence layer -yielding a simpler model if additional complexity is not required, as shown below:

Quantile Outputs
In line with previous work [39], the TFT also generates prediction intervals on top of point forecasts. This is achieved by the simultaneous prediction of various percentiles (e.g. 10 t h , 50 t h and 90 th ) at each time step. Quantile forecasts are generated using linear transformation of the output from the temporal fusion decoder: where W q ∈ R 1×d , b q ∈ R are linear coefficients for the specified quantile q. We note that forecasts are only generated for horizons in the future -i.e. τ ∈ {1, . . . , τ max }.

TRAINING PROCEDURE
As per [39], the TFT is trained by jointly minimizing the quantile loss terms summed across all quantile outputs: QL(y,ŷ, q) = q(y −ŷ) where Ω is the domain of training data containing M samples, W represents the weights of the TFT, Q = {0.1, 0.5, 0.9} is the set of output quantiles, and (.) + = max(0, .). For out-of-sample testing, we evaluate the normalized quantile losses across the entire forecasting horizon -focusing on P50 and P90 risk for consistency with previous work [15,25,31]: whereΩ is the domain of test samples. For additional information, full details on hyperparameter optimization and training can be found in Appendix B.

PERFORMANCE EVALUATION 6.1 Datasets
We choose datasets to reflect commonly observed characteristics across a wide range of challenging multi-horizon forecasting problems. To establish a baseline and position with respect to prior academic work, we first evaluate performance on the Electricity and Traffic datasets used in [15,25,31] -which focus on simpler univariate time series containing known inputs only alongside the target. Next, the Retail dataset helps us benchmark the model using the full range of complex inputs observed in multi-horizon prediction applications (see Section 3) -including static metadata and observed time-varying inputs. Finally, to evaluate robustness to over-fitting on smaller noisy datasets, we consider the financial application of volatility forecasting -using a dataset much smaller than others. Broad descriptions of each dataset can be found below: Electricity. The UCI Electricity Load Diagrams Dataset contains the electricity consumption of 370 customers -aggregated on an hourly level as in [41]. In accordance with [15], we use the past week of data (i.e. 168 hours) to forecast consumption over the next day (i.e. 24 hours).
Traffic. The UCI PEM-SF Traffic Dataset describes the occupancy rate (with y t ∈ [0, 1])for 440 San Francisco Bay Area freewaysas in [41]. This is also aggregated on an hourly level as per the electricity dataset, with the same look back window and forecast horizon.
Retail. Favorita Grocery Sales Dataset from the Kaggle competition [14], that combines metadata for different products and the stores, along with other exogenous time-varying inputs sampled at the daily level. We forecast log product sales 30 days into the future, using 90 days of historical information.
Volatility. The OMI realized library [19] contains daily realized volatility values of 31 stock indices computed from intraday data, along with their daily returns. For our experiments, we consider forecasts over the next week (i.e. 5 business days) using information over the past year (i.e. 252 business days).
For each dataset, we partition all time series into 3 sections -a training set for network calibration, a validation set for hyperparameter optimisation, and a hold-out test set for performance evaluation. Full details on the feature engineering steps and train/test splits are provided for each dataset in Appendix A.

Benchmarks
We extensively compare the TFT to a wide range of machine learning models for multi-horizon forecasting, based on the categories described in Section 2. Hyperparameter optimization is conducted using random search over a pre-defined search space, using the same number of iterations across all benchmarks for a give dataset. Additional details on benchmark model training are also included in Appendix B for reference.
Direct methods. As the TFT falls within this class of multi-horizon models, we primarily focus comparisons on deep learning methods  which directly generate prediction at future horizons. This specifically includes 1) simple sequence-to-sequence models with global contexts (Seq2Seq), and 2) the Multi-horizon Quantile Recurrent Forecaster (MQRNN) -both of which are described in [39].
Iterative methods. To position with respect to the rich body of work on iterative models, we evaluate the TFT using the same setup as [15] for the Electricity and Traffic datasets. This extends the results from [25] for 1) DeepAR models [15], 2) Deep State Space Models (DSSM) [31], and 3) the Transformer-based architecture of [25] with local convolutional processing -which refer to as ConvTrans. For more complex datasets, we focus on the ConvTrans model given its strong outperformance over other iterative models in prior work. As models in this category require knowledge of all inputs in the future to generate predictions, we accommodate this for complex datasets by imputing unknown inputs with their last available value.

Results and Discussion
Tables 1 show that the TFT significantly outperforms all benchmarks over the variety of datasets described in Section 6.1. For median forecasts, the TFT yields 7% lower P50 and 9% lower P90 losses on average compared to the next best model -demonstrating the benefits of explicitly aligning the architecture with the general multi-horizon forecasting problem. Further ablation analyses can be found in Appendix C.
Comparing direct and iterative model performance, we observe the importance of accounting for the observed inputs -noting the poorer results of ConvTrans on complex datasets where observed input imputation is required (i.e. Volatility and Retail). Furthermore, the benefits of quantile regression are also observed when targets are not modelled well by conditional Gaussians with directly method outperforming in those scenarios. This can be seen, for example, from the Traffic dataset where target distribution is significantly skewed -with more than 90% of occupancy rates falling between 0 and 0.1, and the remainder distributed evenly until 1.0.

INTERPRETABILITY USE CASES
Having established the performance benefits of our model, we next demonstrate how to analyze components of the TFT to interpret the general relationships it has learned. We demonstrate three interpretability use-cases: 1) examining the importance of each input variable in prediction, 2) visualizing persistent temporal patterns, and 3) identifying any regimes or events that lead to significant changes in temporal dynamics. In contrast with other examples of attention-based interpretability [1,25,33], which zoom in on interesting but instance-specific examples, we note that our methods focus on ways to aggregate the patterns across the entire dataset -allowing us to extract generalizable insights about temporal dynamics.

Analyzing Variable Importance
We first quantify variable importance by analyzing the variable selection weights described in Section 4.2. Concretely, we aggregate selection weights (i.e. v (j) χ t in Eq. (8)) for each variable across our entire test set, recording the 10 th , 50 th and 90 th percentiles of each sampling distribution. Given its wide range of inputs, we present results on the Retail dataset in Table 2 -with the remainder presented in Appendix D.1. Overall, the TFT focuses on only a subset of key inputs that significantly contribute to predictions.

Visualizing Persistent Temporal Patterns
The analysis of persistent temporal patterns is often key to understanding the time-dependent relationships present in a given dataset. For instance, lag models are frequently adopted to study length of time required for an intervention to take effect [12]   as the impact of a government's increase in public expenditure on the resultant growth in Gross National Product [4]. Seasonality models are also commonly used in econometrics to identify periodic patterns in a target-of-interest [21] and measure the length of each cycle. Using the attention weights present in the self-attention layer of the temporal fusion decoder, we present a method below to identify similar persistent patterns -by measuring the contributions of features at fixed lags in the past on forecasts at various horizons. Combining Eq. (14) and (19), we see that the self-attention layer contains a matrix of attention weights at each forecast time t -i.e. A(ϕ(t), ϕ(t)). As such, multi-head attention outputs at each forecast horizon τ (i.e. β(t, τ )) can be described as an attention-weighted sum of lower level features at each position n: where α(t, n, τ ) is the (τ , n)-th element ofÃ(ϕ(t), ϕ(t)), andθ (t, n) is a row ofΘ(t) = Θ(t)W V . Due to decoder masking, we also note that α(t, i, j) = 0, ∀i > j. For each forecast horizon τ , the importance of a previous time point n < τ can hence be determined by analyzing distributions of α(t, n, τ ) across all time steps and entities. We present results for the Traffic dataset below, with similar findings on Electricity and Retail presented in Appendix D.2. Fig. 3 shows the temporal patterns learned by the TFT -with the top chart recording the mean along the 10 t h , 50 t h and 90 th percentiles of the attention weights for one-step-ahead forecasts (i. e. α(t, 1, τ )) over the test set, and the average attention weights for various horizons (i.e. τ ∈ {5, 10, 15, 20}) on the bottom. Based on the regularly-spaced peaks at 24-hour intervals, we can infer that the TFT has learned a strong daily seasonal pattern for predictions -placing the largest attention on the same hour of preceding days.

Identifying Regimes & Significant Events
Apart from persistent patterns, identifying sudden changes in temporal patterns can also be very useful, as temporary shifts can occur due to the presence of significant regimes or events. For instance, regime-switching behavior has been widely documented in financial markets [2], with returns characteristics -such as volatilitybeing observed to change abruptly between regimes. Firstly, for a given entity, we define the average attention pattern per forecast horizon to be: and then constructᾱ (τ ) = [ᾱ(−k, τ ), . . . ,ᾱ(τ max , τ )] T . To compare similarities between attention weight vectors, we use the distance metric proposed by [9]: where ρ(p, q) = j √ p j q j is the Bhattacharya coefficient [22] measuring the overlap between discrete distributions -with p j , q j being elements of probability vectors p, q respectively. For each entity, significant shifts in temporal dynamics are then measured using the distance between attention vectors at each point with the average pattern, aggregated for all horizons as below: where α (t, τ ) = [α(t, −k, τ ), . . . , α(t, τ max , τ )] T . We use the volatility dataset as a test case for regime identification, specifically applying our distance metric to the attention patterns for the S&P 500 index over our training period from 2001 to 2015. Plotting dist(t) against the target (i.e. log realized volatility) in the bottom chart of Fig. 4, significant deviations in attention patterns can be observed around periods of high volatility -corresponding to the peaks observed in dist(t). From the plots, we can see that the TFT appears to alter its behaviour between regimes -placing equal attention across historical inputs when volatility is low, while attending more to sharp trend changes during high volatility periods -suggesting differences in temporal dynamics learned in each.

CONCLUSIONS
We introduce the Temporal Fusion Transformer (TFT) -a novel attention-based deep neural network model for interpretable highperformance multi-horizon time series forecasting. The TFT utilizes specialized components to handle the full range of inputs typically present in multi-horizon forecasting problems (i.e. static covariates, a priori known inputs, and observed inputs). Specifically, these include: 1) sequence-to-sequence and attention based temporal processing components that capture time-varying relationships at different timescales, 2) static covariate encoders that allow the network to condition temporal forecasts on static metadata, 3) gating components that enable skipping over any parts of the network that are unnecessary for a given dataset, 4) variable selection networks that select relevant input features at each time step, and 5) quantile predictions to obtain output intervals across all prediction horizons. Through tests on a series of real-world datasets, we show that the TFT achieves state-of-the-art forecasting performance on both simple datasets that contain only known inputs, and complex datasets which encompass the full range of possible inputs. Finally, we investigate the general relationships learned by the TFT through a series of interpretability use-cases -proposing novel methods to use the TFT to 1) analyze important variables for a given prediction problem, 2) visualize persistent temporal relationships learned (e.g. seasonality), and 3) identify significant regimes present in the dataset. 9

APPENDIX A DATASET DESCRIPTION
Additional details on each dataset can be found below. We provide all the sufficient information on feature pre-processing and train/test splits to ensure reproducibility of our results.
Electricity. Per [15], we use 500k samples taken between 2014-01-01 to 2014-09-01 -using the first 90% for training, and the last 10% as a validation set. Testing is done over the 7 days immediately following the training set -as described in [15,41]. Given the large differences in magnitude between trajectories, we also apply z-score normalization separately to each entity for real-valued inputs. In line with previous work, we consider the electricity usage, dayof-week, hour-of-day and and a time index -i.e. the number of time-steps from the first observation -as real-valued inputs, and treat the entity identifier as a categorical variable.
Traffic. Tests on the Traffic dataset are also kept consistent with previous work, using 500k training samples taken before 2008-06-15 as per [15], and split in the same way as the Electricity dataset. For testing, we use the 7 days immediately following the training set, and z-score normalization was applied across all entities. For inputs, we also take traffic occupancy, day-of-week, hour-of-day and and a time index as real-valued inputs, and the entity identifier as a categorical variable.
Retail. For the retail dataset, we treat each product number-store number pair as a separate entity, with over 135k entities present within the full dataset. For network calibration, the training set is made up of 450k samples taken between 2015-01-01 to 2015-12-01, validation set of 50k samples from the 30 days after the training set, and test set of all entities over the 30-day horizon following the validation set. We use all inputs supplied by the Kaggle competition (full list in the variable importance section of Appendix D.1 -including additional variables for the day-of-week, day-ofmonth, and month. Data is also resampled at regular daily intervals, imputing any missing days using the last available observation. We also include an additional 'open' flag to denote whether data is present on a given day. For holidays, we group national, regional, and local holidays as separate categorical variables. We also apply a log-transform on the sales data, and adopt z-score normalization across all entities. We consider log sales, transactions, oil to be real-valued variables -with the remainder treated as categorical inputs. Volatility. Data is downloaded from 2000-01-03 to 2019-06-28with the training set consisting of data before 2016, the validation set from 2016-2017, and the test set data from 2018 onwards. For the target, we focus on 5-min sub-sampled realized volatility (i.e. the rv5_ss column ), and add daily open-to-close returns as an extra exogenous input. Also, additional variables are included for the day-of-week, day-of-month, week-of-year, and month -along with a 'region' variable for each index (i.e. Americas, Europe or Asia). Finally, a time index is added to denote the number of days from the first day in our training set. We treat all date-related variables (i.e. day-of-week, day-of-month, week-of-year, and month) and the region input as categorical variables. A log transformation is also applied to the target, and all inputs are z-score normalized across all entities.

B TRAINING DETAILS
Hyperparameter optimization is conducted via random search, using 240 iterations for the smaller Volatility dataset, and 60 iterations for the larger Electricity, Traffic and Retail datasets. Full search ranges for all hyperparameters are below, and the optimal TFT hyperparameters can be found in Table 3: To preserve the explainability of interpretable multi-head attention, we adopt only a single stack (i.e. temporal fusion decoder block) for the TFT itself. For ConvTrans, we adopt the same fixed stack size and number of heads used in the original paper [25] setting them to 8 heads and 3 layers respectively. We also used the full attention model of [25], and treated kernel sizes for the CNN local processing layer as a hyperparameter (i.e. kenel size ∈ {1, 3, 6, 9}) -as optimal kernel sizes were observed to be dataset dependent in [25].

C ABLATION ANALYSIS
To quantify the benefits of each of our proposed architectural contribution, we perform an extensive ablation analysis -removing each component from the network as below, and quantifying the percentage increase in loss versus the original architecture: Gating layers. The effects of gated skip connections are tested by replacing each GLU layer (Eq. (5)) with a simple linear layer passed through an ELU activation function.
Static covariate encoders. The importance of specialized static encoders are tested by setting all context vectors to zero -i.e. c s = c e = c c = c h = 0 -and concatenating all transformed static inputs to all time-dependent past and future inputs.
Variable selection networks. The effects of instance-wise variable selection are tested by replacing the softmax outputs of Eq. 8 with a vector of trainable coefficients, and removing the networks generating the variable selection weights. We retain, however, the variable-wise GRNs (i.e. Eq. (7), maintaining the same degree nonlinear processing as before.
Self-attention layers. The benefits of the self-attention layer are quantified by replacing the attention matrix used in the interpretable multi-head attention layer (Eq. 14) with a matrix of trainable parameters W A -i.e.Ã(Q, K) = W A , where W A ∈ R N ×N . This prevents the TFT from attend to different input features at different forecast times, helping us evaluate the importance of instance-wise attention weights.  Table 3: Dataset details and optimal TFT configuration for each dataset.
Sequence-to-sequence layers for local processing. We evaluate the importance of local processing by removing the sequence-tosequence layer of Section 4.5.1 -replacing this with standard positional encoding used in [36].
Ablated networks are trained across for each dataset using the hyperparameters of Table 3, with full results shown in Figure 5. From the charts, the effects on both P50 and P90 losses are found to be similar across all datasets, with all components contributing to performance improvements on the whole. In general, the components responsible for capturing temporal relationships (i.e. local processing and self-attention layers) have the largest impact on performance, with P90 loss increases of > 6% on average and > 20% on select datasets when ablated. Static encoders and variable selection have the next largest impact -increasing P90 losses by more than 2.6% on average and up to 7.9% for specific datasets. Finally, gating layer ablation also significant increases in P90 losses, with a 1.9% increase on average. This is most significantly show on the the volatility dataset (with a 4.1% P90 loss increase), demonstrating the utility of component gating for smaller, noisier datasets.

D ADDITIONAL INTERPRETABILITY RESULTS
On top of the interpretability use cases of Section 7, which highlight our most prominent findings, we also include the remaining results in this section for completeness. Table 4 shows the variable importance scores for the remaining Electricity, Traffic and Volatility datasets. Given that only one static input is present for these datasets, the network allocates full importance for the entity identifier for Electricity and Traffic, as well as for the region input for Volatility. We also observe two general types of import time-dependent inputs -those related to past values of the target as before, and those related to calendar effects. For instance, the hour-of-day plays a significant roles for Electricity and Traffic datasets, echoing the daily seasonality observed in the next section. In the Volatility dataset, the day-of-month is observe to play a significant role in future inputs -potentially reflecting turn-of-month effects [17].

D.2 Persistent Temporal Patterns
Fig . 6 shows the attention weight patterns across all datasets, and extends the results of Section 7.2. We observe that the three datasets exhibit a seasonal pattern, with clear attention spikes at daily intervals observed for the Electricity and Traffic datasets, and a slightly weaker weekly patterns for the Retail dataset. No strong persistent patterns were observed for the Volatility datasets however, with attention weights equally distributed across all positions on average. This resembles a moving average filter at the feature level, andgiven the high degree of randomness associated with the volatility process -could be useful in extracting the trend over the entire period by smoothing out high-frequency noise. Results per dataset shown on the left, and the range across datasets shown on the right. While the precise importance of each is dataset-specific, all components contribute significantly on the whole -with the maximum percentage increase over all datasets ranging from 3.6% to 23.4% for P50 losses, and similarly from 4.1% to 28.4% for P90 losses.
10% 50% 90%  Table 4: Variable importance scores for the Electricity, Traffic and Volatility datasets. The most significant variable of each input category is highlighted in purple. As before, past values of the target play a significant role -being the top 1 or 2 most significant past input across datasets. The role of seasonality can also be seen in Electricity and Traffic, where the past and future values of the hour-of-day is important for forecasts.