Beyond point forecasts: Uncertainty quantification in tensor extrapolation for relational time series data

Abstract

Contemporary data sets are frequently relational in nature. In retail, for example, data sets are more granular than traditional data, often indexing individual products, outlets, or even customers, rather than aggregating them at the group level. Tensor extrapolation aims to forecast relational time series data; it combines tensor decompositions and time series extrapolation. However, previous approaches to tensor extrapolation are restricted to point forecasts and fail to quantify the uncertainty associated with the forecasts. This paper provides prediction intervals within which we expect the true future values to lie with a specified probability. Numerical experiments suggest the effectiveness of the proposed method in terms of forecast performance. Moreover, the approach proves to be up to 77 times faster than the univariate competitor. For this reason, it exerts positive influence on both inclusiveness and the environment.

Keywords:

• Computational cost
• Forecast performance
• Prediction interval
• Relational data
• Tensor decomposition

1 Introduction

In general, variables can be absolute or relational (Wasserman and Faust 1994). The degree of environmental awareness of a person, their gender, education, and age are absolute variables. If, however, a variable focuses on the relationship among entities, it is called a relational variable. Examples of relational variables include the intensity of communications between members of a group, the extent of economic transactions between corporations, and the existence or nonexistence of trade or treaties between countries. Relational variables are often dynamic, i.e., they evolve over time.

Across numerous domains, contemporary data sets are relational by nature (Kitchin 2014; Müller et al. 2016). In retail, for instance, data sets are more granular than traditional data, often capturing information at the level of individual items, outlets, or even users, rather than aggregating them at the group level (George et al. 2016; Fildes, Ma, and Kolassa 2022). As a result, there are different types of dependencies between the variables of interest (e.g., dependencies among stores, dependencies within products). The relational character of the data is, however, often disregarded, particularly in prediction tasks. Instead, univariate extrapolation approaches are applied, which cannot address inter-series dependencies. In addition, available multivariate forecasting methods (e.g., vector autoregressions) are limited to low-dimensional settings and are, hence, not suitable for practical use in large-scale forecasting problems (Hyndman and Athanasopoulos 2018; De Stefani and Bontempi 2021).

Tensor extrapolation is an attempt to forecast relational time series data using multi-linear algebra. The process unfolds as follows: Multi-way data are organized in the form of tensors, i.e., multi-dimensional arrays. Tensor decompositions are then employed to identify periodic patterns in the data. Afterwards, these patterns serve as input for time series methods. Tensor extrapolation originated in the work of Dunlavy, Kolda, and Acar (2011) and Spiegel et al. (2012), initially limited to preselected time series approaches and binary data. Only recently, Schosser (2020, 2022b) laid the groundwork for applications in large-scale forecasting problems.

The purpose of forecasting is to facilitate decision making in the face of uncertainty. To this end, forecasts shall provide an unbiased guess at what is most likely to happen, the so-called point forecast. Increasingly, an estimate of the uncertainty associated with this forecast is called for (Makridakis, Hyndman, and Petropoulos 2020; Makridakis, Spiliotis, and Assimakopoulos 2022b; Makridakis et al. 2022c). So far, tensor extrapolation is restricted to point forecasts, but fails to quantify the respective uncertainty. The paper at hand intends to fill this gap and, hence, contributes to the literature on tensor extrapolation. It aims at prediction intervals within which we expect the true future values to lie with a specified probability. However, in a broader perspective, there is a pressing need to ensure efficient implementation.

Until recently, the forecasting community placed strong (at times, exclusive) emphasis on achieving “state-of-the-art” results. Typically, accuracy (or similar) measures were reported, while any mention of cost or efficiency was omitted (Nikolopoulos and Petropoulos 2018; Petropoulos et al. 2022). Notwithstanding the benefits of improving model performance, this narrow focus ignores the economic, environmental, and social cost of reaching the reported results (Schwartz et al. 2020; Wu et al. 2022). Contemporary forecasting technology is already very expensive to train or execute, which restricts the ability of many researchers to study it, and of practitioners to adopt it. What is more, the computing landscape undergoes a profound transition as improvements in central processing units (CPUs) are slowing down, and an ever-increasing part of applications moves to specialized processors (Thompson and Spanuth 2021). Taken together, rising computational burdens and the switch toward customized, high-end hardware threaten to fragment the field. Another critical aspect is the environmental footprint left in deriving forecasts. Escalating computational experiments greatly increase the amount of energy consumed and, consequently, the extent of natural resources spent.

Meeting the challenges outlined above necessitates a reevaluation of our assessment standards. We advocate for the integration of efficiency as a prevalent criterion, complementing the traditional focus on forecast performance. This strategic shift aligns seemlessly with our commitment to fostering environmentally friendly and inclusive forecasting practices. The average runtime for generating a result serves as a natural proxy for efficiency. All else being equal, a faster method inherently implies a reduction in computational workload.

When expanding tensor extrapolation, we not only seek to improve forecast performance but also emphasize efficiency. The derivation of prediction intervals frequently emerges as one of the most computationally demanding elements in the forecasting workflow, particularly where simulations are conducted. In recognition of this, we take a step back, embracing a heuristic approach to articulate the uncertainty in forecasts. This decision proves fruitful: Our calculations reveal that tensor extrapolation consistently outperforms the univariate baseline while incurring far lower computational cost.

The remainder of the paper progresses as follows. Section 2 briefly reviews related work. Section 3 presents the proposed method. It provides background material on tensor analysis, introduces the state of the art of tensor extrapolation, and explains our extensions. Section 4 describes the experiments conducted. This includes the decision on the performance metrics, the technical realization, and the application to empirical and synthetic data. The results are discussed in Section 5. Finally, Section 6 concludes the paper.

2 Related work

For an introduction to time series analysis, we suggest Box and Jenkins (1976). A more contemporary overview of the forecasting literature is given by Hyndman and Athanasopoulos (2018). The following areas are highly significant to our problem: tensor autoregression, forecasting competitions, and meta-learning for time series forecasting.

2.1 Tensor autoregression

During the last few years, tensor-based techniques have received growing attention within the statistics community (Bi et al. 2021). For time series data, we can roughly distinguish three categories of approaches. In the first category, tensor analysis is used as a building block in deep learning models, specifically to accelerate and stabilize the training of neural networks (Yu et al. 2017; Deng et al. 2024). The second category involves transforming low-dimensional time series data into higher-dimensional data through a process known as tensorization (Shi et al. 2020; Barsbey and Cemgil 2023). The third category pertains to relational data, meaning data that can naturally be represented as a tensor (Schosser 2022b). A complete survey of the highly dynamic literature is beyond the scope of this paper. Instead, we give examples of work in the second and third class. In comparison with traditional forecasting techniques, autoregressive methods are of particular interest. Hill, Li, and Schneider (2021) recommend a regression model for univariate time series with known seasonal periodicity. Within their framework, past observations of the dependent variable are stacked in tensor format. In contrast, Yu, Rao, and Dhillon (2016) develop a model for relational time series data. Unlike tensor extrapolation, where tensor factorizations are combined with time series extrapolation methods, Yu, Rao, and Dhillon (2016) introduce an autoregressive regularizer and thus impose specific temporal structure on the factorization products. Both Hill, Li, and Schneider (2021) and Yu, Rao, and Dhillon (2016) focus on point forecasts. Tensor extrapolation employs a general class of state-space models, offering a more comprehensive setting compared to autoregressive approaches.

2.2 Forecasting competitions

Competitions (or “challenges”, “common tasks”) are omnipresent in machine learning, statistics, and related fields (Feuerverger, He, and Khatri 2012; Donoho 2017). They lower barriers to entry by providing clearly defined tasks and the resources needed to accomplish them. Moreover, they create research communities with shared goals and assumptions and, perhaps most important, offer proof of incremental scientific progress (Liberman 2010). The so-called M competitions are the most influential and frequently cited competitions in the field of time series forecasting (Hyndman 2020). The focus of the recent M5 competition is on a retail sales forecasting application, examining both accuracy and uncertainty estimation (Makridakis, Spiliotis, and Assimakopoulos 2022b; Makridakis et al. 2022c). The data provided comprise grouped, correlated time series, organized in a hierarchical structure. Depending on the level of aggregation, different forecasting methods stand out. For the top and middle parts of the hierarchy, machine learning-based approaches offer considerable gains in performance over pure statistical methods. At the lower end of the hierarchy, especially at the product-store level, capturing time series patterns becomes more challenging due to the presence of irregularities. This is where traditional methods come into play. We demonstrate the adequacy of tensor extrapolation at a comparable level of granularity.

2.3 Meta-learning for time series forecasting

In relevant literature, there are several classes of forecasting models such as exponential smoothing, Autoregressive Integrated Moving Average (ARIMA), or deep neural networks. Within individual model families, automation of model selection and tuning of forecasting algorithms has been studied in detail (Hyndman and Khandakar 2008; Salinas et al. 2020). Meta-learning for time series forecasting is a natural progression of these efforts. It aims to efficiently train, optimize, and choose the best-performing model among various classes of predictive models for a given data set. Meta-learning platforms have recently been proposed by Gastinger et al. (2021), Ma and Fildes (2021), and Shah et al. (2021). However, these approaches do not feature routines specifically designed for relational data. We illustrate that accounting for the relational character of the data pays off, and suggest tensor extrapolation as an additional building block in meta-learning frameworks.

3 Proposed method

This section covers background information on tensor analysis. In addition, it shows the current state of tensor extrapolation and explains our extensions.

3.1 Tensor analysis

Multi-way arrays, or tensors, are multi-dimensional collections of numbers. The dimensions are known as ways, orders, or modes of a tensor. Using this terminology, scalars, vectors, and matrices can be interpreted as zero-order, first-order, and second-order tensors, respectively. Tensors of order three and higher are called higher-order tensors (Cichocki et al. 2009; Kolda and Bader 2009). In the simplest high-dimensional case, the tensor can be thought of as a “data cube” (Papalexakis, Faloutsos, and Sidiropoulos 2016; Rabanser, Shchur, and Günnemann 2017). This case should be formalized in the following: Let $I,J,K\in \mathbb{N}$ denote index upper bounds, i.e., the number of entities in the modes of interest; a third-order tensor is denoted by $\underset{¯}{\mathbf{X}}\in {\mathbb{R}}^{I\times J\times K}$. Here, x_ijk describes the entry in the ith row, jth column, and kth tube of $\underset{¯}{\mathbf{X}}$. The modes of $\underset{¯}{\mathbf{X}}$ are referred to as mode A, mode B, and mode C, respectively. Throughout this section, we will restrict ourselves to third-order tensors for reasons of simplicity. Nevertheless, the concepts introduced naturally extend to tensors of order four and higher. It should be noted that most of the notation and terminology used in our deliberations is borrowed from Kiers (2000). A brief summary of the tensor notation is presented in Table 1.

Table 1 Tensor notation summary.

Symbol	Description
p	scalar (lowercase letter)
a	vector (boldface lowercase letter)
ai	element i of vector a
A	matrix (boldface capital letter)
${\mathbf{a}}_r$	rth column of matrix A
air	element ir of matrix A
$\underset{¯}{\mathbf{X}}$	higher-order tensor (underlined boldface capital letter)
xijk	element ijk of tensor $\underset{¯}{\mathbf{X}}$

The so-called Candecomp/Parafac (CP) decomposition is one of the most common tensor factorizations. It decomposes a tensor as a set of factor matrices. The model was independently proposed by Hitchcock (1927), Carroll and Chang (1970), and Harshman (1970). Optimization entails $$\underset{a_{\mathit{\text{ir}}},b_{\mathit{\text{jr}}},c_{\mathit{\text{kr}}}}{\mathrm{\text{min}}}\sqrt{\sum _{i=1}^I\sum _{j=1}^J\sum _{k=1}^K{\left(x_{\mathit{\text{ijk}}}-\sum _{r=1}^R{a}_{\mathit{\text{ir}}}{b}_{\mathit{\text{jr}}}{c}_{\mathit{\text{kr}}}\right)}^2}.$$

Here, a_ir , b_jr , and c_kr denote elements of the component matrices A (for mode A), B (for mode B), and C (for mode C) of sizes $(I\times R),(J\times R)$, and $(K\times R)$, respectively. The number of components R represents the only hyperparameter to be specified. The optimization problem stated above assumes that the random fluctuation in the tensor data follows a Gaussian distribution. Alternative models can be found in Chi and Kolda (2012) and Papalexakis, Faloutsos, and Sidiropoulos (2016).

3.2 Tensor extrapolation: literature and extensions

Tensor extrapolation seeks to forecast relational time series data. That means, relations for T time steps are given as input and the goal is to predict the relations at future times T + 1, T + 2, $\dots$, T + L. Without loss of generality, the scope of this presentation is limited to the case where the snapshots of relations can be represented in the form of a matrix. Here, tensors provide a straightforward way to integrate the temporal dimension. Consequently, a third-order data array $\underset{¯}{\mathbf{X}}$ of size ($I\times J\times T$) is given, with time being modeled in mode C. Extrapolation entails computing the tensor $\underset{¯}{\widehat{\mathbf{X}}}$ of size ($I\times J\times L$) that includes point estimates concerning future relations. The approach involves the use of tensor decomposition and exponential smoothing. It proceeds as follows: The multi-way data array is decomposed applying Candecomp/Parafac (CP) factorization. Every one of the component matrices gives information about the respective mode, i.e., detects latent structure in the data. For further processing, the “time” component matrix C of size (T × R) is converted into a set of column vectors ${\mathbf{c}}_r$. These vectors capture different periodic patterns, e.g., seasons and trends. The periodic patterns discovered are taken as input for exponential smoothing techniques. In doing so, forecasts are obtained and arranged as columns ${\widehat{\mathbf{c}}}_r$ of the new “time” component matrix $\widehat{\mathbf{C}}$ of size (L × R). Eventually, the tensor $\underset{¯}{\widehat{\mathbf{X}}}$ can be calculated; it contains estimates of future relations.

Tensor extrapolation was originally developed for binary data, indicating the presence or absence of relations of interest. Dunlavy, Kolda, and Acar (2011) use the temporal profiles computed by CP as a basis for the so-called Holt-Winters Method, i.e., exponential smoothing with additive trend and additive seasonality. Spiegel et al. (2012) apply CP decomposition in connection with simple exponential smoothing, meaning exponential smoothing without trend and seasonality. In both articles, the model parameters are hard-coded. Later on, Schosser (2022b) breaks ground for data-driven model selection and estimation in large-scale forecasting problems. First, he draws on an automatic forecasting procedure based on a general class of state-space models encompassing all standard exponential smoothing methods (Hyndman et al. 2002; Hyndman and Khandakar 2008). Model selection and optimization of both smoothing parameters and initial state variables are performed “individually” (Fildes 1989; Petropoulos et al. 2014). This means that they are done based on each individual periodic pattern contained in a column of the matrix C. Second, he investigates how to use the methodology with real-valued data. This highlights the importance of data preprocessing. Schosser (2022a) transfers the framework into a setting characterized by missing values and thus increases the practical relevance of tensor extrapolation. Note that there are adaptations of tensor extrapolation that go beyond the scope of pure time series forecasting. Particular mention should be made of approaches that use additional information from distinct sources, for example in the form of regularized (Bi et al. 2022) or coupled decompositions (Araujo, Ribeiro, and Faloutsos 2017).

To date, tensor extrapolation is used to generate point forecasts. There is no quantification of the uncertainty associated with the predictions. The paper at hand provides an appropriate extension and, consequently, contributes to the literature on the forecasting of relational data. In line with the goal of computational efficiency, we take a heuristic approach. For a single time series, the lower and upper endpoints of the prediction interval are supposed to be given by $${\widehat{x}}_{\mathit{\text{ijl}}}^{\downarrow }={\widehat{x}}_{\mathit{\text{ijl}}}-0.5\lambda \sqrt{l}{\sigma }_{\mathit{\text{ij}}·}^e$$

and $${\widehat{x}}_{\mathit{\text{ijl}}}^{\uparrow }={\widehat{x}}_{\mathit{\text{ijl}}}+0.5\lambda \sqrt{l}{\sigma }_{\mathit{\text{ij}}·}^e,$$ respectively. Here, the point forecast ${\widehat{x}}_{\mathit{\text{ijl}}}$ marks the middle of the range of possible values. The interval width is proportional to ${\sigma }_{\mathit{\text{ij}}·}^e$, i.e., the standard deviation of residual errors in the estimation sample. A characteristic common to prediction intervals is that they increase in width as the forecast horizon increases. The further ahead the forecast, the more uncertainty is associated with it, and, hence, the wider the prediction interval. Our framework ensures this effect; the interval width grows at the square root of the temporal distance $\sqrt{l}$ ($l\in 1,\dots ,L$). The multiplier λ depends on the desired coverage probability. In this paper, we derive 80% and 95% prediction intervals, although any percentage may be used. Inspired by the quantiles of a standardized Gaussian distribution, we specify λ as 1.28 and 1.96, respectively. Since the CP decomposition employed is based on the assumption of Gaussian errors, this imposes no further restrictions on the scope of application. Both the lower and upper interval bounds are organized in tensor form. Algorithm 1 displays our modifications to tensor extrapolation. Numerical experiments suggest the effectiveness of our proposed method. We use both real-world and synthetic data to demonstrate the superiority of the procedure over an established baseline in terms of forecast performance and computational burden.

Algorithm 1.

Tensor extrapolation for large-scale data including point and interval forecasts

Table

Input $\underset{¯}{\mathbf{X}}$, R, λ
Output $\underset{¯}{\widehat{\mathbf{X}}},{\underset{¯}{\widehat{\mathbf{X}}}}^{\downarrow },{\underset{¯}{\widehat{\mathbf{X}}}}^{\uparrow }$
1:	Preprocessing / transformation: ${\underset{¯}{\mathbf{X}}}_{\mathit{\text{trans}}}$ $\leftarrow$ $\underset{¯}{\mathbf{X}}$
2:	CP decomposition: A, B, C $\leftarrow$ ${\underset{¯}{\mathbf{X}}}_{\mathit{\text{trans}}}$, R
3:	Separate column vectors: ${\mathbf{c}}_1,{\mathbf{c}}_2,\dots ,{\mathbf{c}}_R$ $\leftarrow$ C
4:	for each $r\in 1,\dots ,R$ do
5:	Automatic exponential smoothing: ${\mathbf{c}}_r^{\mathit{\text{fit}}},{\widehat{\mathbf{c}}}_r$ $\leftarrow$ ${\mathbf{c}}_r$
6:	end for
7:	Merge column vectors: ${\mathbf{C}}^{\mathit{\text{fit}}}$ $\leftarrow$ ${\mathbf{c}}_1^{\mathit{\text{fit}}},{\mathbf{c}}_2^{\mathit{\text{fit}}},\dots ,{\mathbf{c}}_R^{\mathit{\text{fit}}}$; $\widehat{\mathbf{C}}$ $\leftarrow$ ${\widehat{\mathbf{c}}}_1,{\widehat{\mathbf{c}}}_2,\dots ,{\widehat{\mathbf{c}}}_R$
8:	Calculate tensors: ${\underset{¯}{\mathbf{X}}}_{\mathit{\text{trans}}}^{\mathit{\text{fit}}}$ $\leftarrow$ A, B, ${\mathbf{C}}^{\mathit{\text{fit}}}$; ${\underset{¯}{\widehat{\mathbf{X}}}}_{\mathit{\text{trans}}}$ $\leftarrow$ A, B, $\widehat{\mathbf{C}}$
9:	Back-transformation: ${\underset{¯}{\mathbf{X}}}^{\mathit{\text{fit}}}$ $\leftarrow$ ${\underset{¯}{\mathbf{X}}}_{\mathit{\text{trans}}}^{\mathit{\text{fit}}}$; $\underset{¯}{\widehat{\mathbf{X}}}$ $\leftarrow$ ${\underset{¯}{\widehat{\mathbf{X}}}}_{\mathit{\text{trans}}}$
10:	Calculate prediction errors based on residuals: $\underset{¯}{\mathbf{E}}$ $\leftarrow$ $\underset{¯}{\mathbf{X}},{\underset{¯}{\mathbf{X}}}^{\mathit{\text{fit}}}$
11:	Calculate interval limits: ${\underset{¯}{\widehat{\mathbf{X}}}}^{\downarrow },{\underset{¯}{\widehat{\mathbf{X}}}}^{\uparrow }$ $\leftarrow$ $\underset{¯}{\widehat{\mathbf{X}}},\underset{¯}{\mathbf{E}}$, λ

4 Experimental evaluation

Next, we detail the choice of performance measures, the implementation of our approach, and its application to empirical and synthetic data.

4.1 Performance measures

We rely on the following performance measures to compare our predictions with the observations retained: The Mean Absolute Percentage Error (MAPE) represents a scale-free performance measure for point forecasts. It is given by $$\mathit{\text{MAPE}}=\frac{1}{I·J·L}\sum _{i=1}^I\sum _{j=1}^J\sum _{l=1}^L\frac{\left|x_{\mathit{\text{ijl}}}-{\widehat{x}}_{\mathit{\text{ijl}}}\right|}{\left|x_{\mathit{\text{ijl}}}\right|}·100.$$

The literature on time series forecasting often favors the Symmetric Mean Absolute Percentage Error (sMAPE) and the Mean Absolute Scaled Error (MASE) as performance metrics (Hyndman and Koehler 2006). For the data sets used in the following, MAPE has been proven to be a reliable proxy for sMAPE (Schosser 2020, 2022a, 2022b). With MASE, the denominator is given by the average absolute error of a naïve one-step prediction. Since a naïve forecast does not offer a suitable benchmark for data sets with heterogeneous periodic patterns, this may lead to distortions. Therefore, and in the light of our focus on interval prediction, we concentrate on a single measure of the quality of point forecasts.

The Degree of Miscalibration (DM) is used to compare the performance of the methods on how frequently prediction intervals contain the true values of the series. We first compute the actual coverage, i.e., the percentage of true values included in the (nominally 80% or 95%) prediction intervals. From this, we derive DM, understood as the difference between the actual coverage for the forecasting method under study and the nominal probability level. Where the actual coverage falls short of the nominal probability level (negative difference), we denote the interval as overconfident. Where the reverse relationship holds (positive difference), we refer to the interval as underconfident. Miscalibration typically results from inappropriate interval widths or location of intervals.

The Mean Scaled Interval Score (MSIS) is a more elaborate attempt at quantifying the adequacy of prediction intervals. Formally, $$\mathit{\text{MSIS}}=\frac{1}{I·J·L}\sum _{i=1}^I\sum _{j=1}^J\sum _{l=1}^L\frac{I{S}_{\mathit{\text{ijl}}}}{\left|x_{\mathit{\text{ijl}}}\right|}$$

and $$I{S}_{\mathit{\text{ijl}}}=\left({\widehat{x}}_{\mathit{\text{ijl}}}^{\uparrow }-{\widehat{x}}_{\mathit{\text{ijl}}}^{\downarrow }\right)+p\left({\widehat{x}}_{\mathit{\text{ijl}}}^{\downarrow }-x_{\mathit{\text{ijl}}}\right)1_{\{x_{\mathit{\text{ijl}}}<{\widehat{x}}_{\mathit{\text{ijl}}}^{\downarrow }\}}+p\left(x_{\mathit{\text{ijl}}}-{\widehat{x}}_{\mathit{\text{ijl}}}^{\uparrow }\right)1_{\{x_{\mathit{\text{ijl}}}>{\widehat{x}}_{\mathit{\text{ijl}}}^{\uparrow }\}}$$ hold. Here, $1_{\left\{\right\}}$ describes an indicator function that returns the value one if the condition in curly brackets is fulfilled and zero otherwise. A forecast method is rewarded for narrow prediction intervals, and incurs a penalty, the size of which depends on the constant p, if the observation misses the interval (Gneiting and Raftery 2007; Landon and Singpurwalla 2008). Typically, the parameter p is chosen in an ad-hoc manner (Makridakis et al. 2022c); this is why we include a range of values in our analysis. To allow for aggregation over the entire data set, scaling is required. Within the M4 and M5 competitions, the average absolute error of a naïve one-step forecast was used as a scaling factor. For reasons discussed above, we discard this approach. Instead, MSIS is scaled similarly to MAPE.

The Mean Scaled Interval Width (MSIW) indicates the degree of uncertainty associated with the forecasts. Wider intervals reflect higher uncertainty. In our case, it is given by $$\mathit{\text{MSIW}}=\frac{1}{I·J·L}\sum _{i=1}^I\sum _{j=1}^J\sum _{l=1}^L\frac{{\widehat{x}}_{\mathit{\text{ijl}}}^{\uparrow }-{\widehat{x}}_{\mathit{\text{ijl}}}^{\downarrow }}{\left|x_{\mathit{\text{ijl}}}\right|}.$$

For MSIW, we apply the same scaling as with MAPE.

4.2 Implementation

For our calculations, we use the programming language Python. In particular, tensor extrapolation employs the functions parafac (library TensorLy; Kossaifi et al. 2019) and ETSModel (library statsmodels; Seabold and Perktold 2010). When parafac is called, the number of components R must be declared. The function ETSModel provides a Python-based implementation of the automatic exponential smoothing algorithm designed by Hyndman et al. (2002) and Hyndman and Khandakar (2008). This algorithm offers a variety of possible models and has been proven to be very powerful in comparisons with other forecasting methods (Januschowski et al. 2020; Gastinger et al. 2021; Karlsson Rosenblad 2021). Further, in relation to more sophisticated forecasting techniques, the requirements on data availability and computational resources are fairly low (De Stefani and Bontempi 2021; Schosser 2022b). As our baseline, we employ univariate extrapolation through ETSModel, producing both point forecasts and prediction intervals. Diverging from tensor extrapolation, which explicitly models the entire data set in tensor shape, this univariate approach forecasts individual time series one by one, effectively partitioning the tensor into vectors. The choice of this baseline is driven by the outcomes of the M5 competition, where traditional statistical methods demonstrated superiority in handling disaggregated data (Makridakis, Spiliotis, and Assimakopoulos 2022b; Makridakis et al. 2022c). To quantify the computational burden associated with the methods in question, we resort to the Python library timeit. The measurements refer to a commodity notebook with AMD Ryzen 5 4500U 6x2.40 GHz processor and 8 GB RAM.

4.3 Application: empirical data

Following Schosser (2022b), our empirical study relies on data provided by the United Nations. Trade data between countries of the world are accessible at the UN Comtrade website https://comtradeplus.un.org. We choose a set of 30 countries, usually large or developed countries with high gross domestic products and trade volumes. As a trade category, we select “export of goods.” We analyze monthly data collected from 2012 to 2016. Since the exports of a country to itself are not defined, we end up with a total of 870 time series, each comprising 60 observations. The data are partitioned into an estimation sample and a hold-out sample. The latter consists of the twelve most recent observations that are used to evaluate the forecast accuracy. Estimation sample and hold-out sample are organized in tensor form. We thus obtain ${\underset{¯}{\mathbf{X}}}^{\mathit{\text{est}}}$ of size $(30\times 30\times 48)$ and ${\underset{¯}{\mathbf{X}}}^{\mathit{\text{hold}}}$ of size $(30\times 30\times 12)$, respectively. The methods under study are implemented, or trained, on the estimation sample. Point forecasts are produced for the whole of the hold-out sample and arranged as tensor $\underset{¯}{\widehat{\mathbf{X}}}$ of size $(30\times 30\times 12)$. In addition, we determine the interval boundaries and organize them into the tensors ${\underset{¯}{\widehat{\mathbf{X}}}}^{\downarrow }$ and ${\underset{¯}{\widehat{\mathbf{X}}}}^{\uparrow }$ (of corresponding size). It is worth mentioning that Hoff (2011) and Ma, Qin, and Xia (2023) examine United Nations trade data within a tensor-based setting. Assuming corresponding multivariate Gaussian distributions, the authors detect correlations among the index sets of the data. However, the information obtained is not utilized for prediction purposes.

Table 2 displays our results. In accordance with Schosser (2022b), tensor extrapolation consistently outperforms univariate exponential smoothing in terms of MAPE, i.e., the quality of point forecasts. However, the findings concerning the associated prediction intervals are also very encouraging. On the basis of MSIS, our proposed approach dominates the baseline. This holds true irrespective of the parameter p. DM remains within reasonable limits. With 80% prediction intervals, there are modest advantages for tensor extrapolation. With 95% prediction intervals, the inverse relationship applies. As measured by MSIW, the proposed method generates slightly wider prediction intervals. Obviously, this drawback is more than compensated by the better location of the intervals.

Table 2 Forecasting performance based on empirical data in terms of MAPE, DM, MSIS, and MSIW.

		Tensor extrapolation
Measure	ETSModel	R = 4	R = 5	R = 6	R = 7	R = 8
MAPE	272.7493	266.4299	266.7067	264.6828	195.7938	193.7586
80% pred. int.
DM	0.0603	–0.0017	–0.0043	–0.0077	–0.0207	–0.0116
MSIS (p = 5)	12.0713	11.5867	11.5985	11.5315	7.3062	7.2776
MSIS (p = 10)	21.9902	20.7811	20.8060	20.6866	11.8007	11.7394
MSIS (p = 40)	81.5042	75.9479	76.0508	75.6170	38.7677	38.5098
MSIW	2.1523	2.3922	2.3911	2.3764	2.8117	2.8159
95% pred. int.
DM	–0.0003	–0.0488	–0.0499	–0.0499	–0.0552	–0.0475
MSIS (p = 5)	11.7629	11.3958	11.3988	11.3254	6.7231	6.7076
MSIS (p = 10)	20.2342	19.1285	19.1363	19.0119	9.1407	9.1033
MSIS (p = 40)	71.0619	65.5247	65.5610	65.1310	23.6463	23.4776
MSIW	3.2916	3.6631	3.6614	3.6389	4.3054	4.3118

The estimation sample includes the years 2012–2015, the hold-out sample covers the year 2016. The results obtained from tensor extrapolation are listed with respect to the selected hyperparameter, the number of components R. The constant p affects the penalty in case the observation misses the prediction interval.

As a next step, we delve into the computational efforts associated with the methods employed. In case of 100 executions, the average runtime for tensor extrapolation with R = 4 components is a mere 0.19 seconds. This calculation includes 0.12 seconds attributed to the tensor decomposition process. With R = 4 components, the automatic exponential smoothing algorithm is applied on an equivalent number of time series. In conjunction with the heuristic derivation of prediction intervals, this operation consumes an additional 0.07 seconds. In stark contrast, the baseline approach exhibits a significantly higher average runtime of 13.40 seconds. This value stems from 870 calls to the function ETSModel, one for each time series within the data set. It is crucial to note that the application of the automatic exponential smoothing algorithm, encompassing the selection and estimation of an appropriate exponential smoothing method, results in highly volatile runtimes contingent on the characteristics of the respective time series. In essence, tensor extrapolation substantially reduces the dimensions of the forecasting problem, rendering it 70 times faster.

4.4 Application: synthetic data

Regardless of the promising results of the previous subsection, our approach should be applied to more than just one data set. An obvious use case would be the data set underlying the M5 competition. However, this data record exhibits two features that preclude the use of tensor extrapolation for the moment. First, apart from time series data, it provides exogenous (or explanatory) variables, such as holidays, product prices, promotions, and social events. Second, the competition involved series that display intermittency, meaning sporadic demand, including many zeros (Syntetos, Boylan, and Croston 2005). Appropriate extensions of our framework to these problems are the subject of future work.

For further testing, it is essential to base the assessment on time series with both diverse and controllable characteristics. The relevant literature increasingly calls for the use of synthetic data (Kang, Hyndman, and Li 2020; Makridakis et al. 2022a) or already resorts to such data (Hill, Li, and Schneider 2021; Shah et al. 2021; Das et al. 2024). With this in mind, we create a comprehensive data set closely related to the field of retailing. Relational data appear, for example, in companies that sustain relationships with a large number of business partners or users and offer a multitude of products or services. Therefore, our simulated data should be understood as sales volume by user, product, and time. We presume that the data can be modeled by R = 4 components and strive for a tensor of size $(I\times J\times T)=(160\times 120\times 100)$. In this way, we arrive at a total of 19,200 time series with 100 observations each.

When constructing our synthetic data set, we adhere to the procedure introduced by Dunlavy, Kolda, and Acar (2011) and adapted by Schosser (2022a, 2022b). Apart from missing values, the resulting data correspond to those used by Schosser (2022a). To start with, we generate the component matrices. In this regard, we do not encounter specific limitations. Notably, the CP model does not impose orthogonality constraints on the component matrices (Dunlavy, Kolda, and Acar 2011; Papalexakis, Faloutsos, and Sidiropoulos 2016). The matrices A and B of size $(I\times R)$ and $(J\times R)$, respectively, can be considered as “entity participation” matrices. That is to say, column ${\mathbf{a}}_r$ (resp. ${\mathbf{b}}_r$) is the vector of participation levels of all the entities in component r. Matrix A demonstrates the participation of the users: Here, the users are expected to respond in a different way to a specific product-time combination. We measure the magnitude of this reaction according to the density function of a Gaussian distribution (cf. Fig. 1(a)). Matrix B shows the product participation: Groups of products can be distinguished that respond similarly to the same time and user effect combination (cf. Fig. 1(b)). The columns of matrix C of size $(T\times R)$ record a number of periodic patterns. We create a sinusoidal pattern, a trend, a structural break, and another break (cf. Fig. 1(c)). Aggregation of the matrices A, B, and C yields a noise-free version of the tensor. Eventually, we add Gaussian noise to any entry. The standard deviation of the error term is taken to equal half the mean of the noise-free tensor entries. We split our simulated data into an estimation sample and a hold-out sample. The latter is made up of the 20 most recent observations. Once more, the forecasts are produced for the whole of the hold-out sample.

Fig. 1 Generation of synthetic data.

Table 3 shows our results. With regard to the accuracy of point forecasts, tensor extrapolation easily beats the baseline. The related forecast intervals are, however, also highly convincing. In terms of miscalibration, our proposed method dominates. What is more, it generates slightly narrower prediction intervals. Finally, tensor extrapolation exhibits better interval scores in almost all parameter settings.

Table 3 Forecasting performance based on synthetic data in terms of MAPE, DM, MSIS, and MSIW.

		Tensor extrapolation
Measure	ETSModel	R = 3	R = 4	R = 5
MAPE	229.2167	185.4018	179.8527	181.2276
80% pred. int.
DM	0.0848	0.0438	0.0469	0.0468
MSIS (p = 5)	7.8454	7.2656	7.0958	7.1177
MSIS (p = 10)	8.7375	8.2508	7.9383	7.9805
MSIS (p = 40)	14.0899	14.1620	12.9931	13.1573
MSIW	6.9534	6.2804	6.2533	6.2550
95% pred. int.
DM	0.0292	–0.0016	–0.0007	–0.0010
MSIS (p = 5)	10.7367	9.9298	9.8693	9.8826
MSIS (p = 10)	10.8392	10.2426	10.1631	10.1873
MSIS (p = 40)	11.4539	12.1198	11.9261	12.0154
MSIW	10.6343	9.6169	9.5754	9.5779

The estimation sample includes the first 80 observations, and the hold-out sample covers the next 20 observations. The results obtained from tensor extrapolation are listed with respect to the selected hyperparameter, the number of components R. The constant p affects the penalty in case the observation misses the prediction interval.

Given the variety and volume of our synthetic data, the computational cost are more pronounced. Based on 100 executions, the average runtime of tensor extrapolation with R = 4 components equals 4.34 seconds. Specifically, 4.20 seconds are dedicated to tensor decomposition. Operating the automatic exponential smoothing algorithm on R = 4 time series, along with the calculation of prediction intervals, adds 0.14 seconds to the process. The competitor takes on average 334.34 seconds. This figure is derived from 19,200 calls to the function ETSModel, each corresponding to a time series within the data set. The method advocated is, thus, 77 times faster.

5 Discussion

The findings of our study highlight the superior performance of tensor extrapolation, excelling in both the accuracy of point forecasts and the estimation of associated uncertainty. Moreover, our method demonstrates a noteworthy speed advantage over the baseline, resulting in a significantly reduced computational burden. This advantage is expected to persist across various setups, indicating that the relative differences in computational time will endure regardless of the system used. Therefore, tensor extrapolation disrupts the conventional tradeoff between forecast performance and computational efficiency. Essentially, even when the objective is to quantify uncertainty, accounting for the relational nature of the data proves beneficial. This approach enables us to capture intrinsic structural information and various types of dependencies, thereby enhancing the overall effectiveness of the forecasting workflow.

In selecting our data records, we emphasized diverse and controllable features. Nevertheless, one aspect we did not address was the length of the time series to be predicted—the number of observations at our disposal for training the methods. Generally, the availability of time series data depends on the specifics of the domain and application. In retail, for instance, product life cycles are becoming shorter and shorter, limiting the access to historical sales figures to a specific period (Fildes, Ma, and Kolassa 2022). As demonstrated, tensor extrapolation is well-equipped for such scenarios (see also Schosser 2022b). However, there is no reason to question the suitability of our method for longer time series; the capacity of CP decomposition to identify relevant underlying patterns remains unaffected.

6 Conclusions

In spite of the opportunities emerging from today’s “big data,” the relational character of many time series is largely neglected in forecasting tasks. Recently, tensor extrapolation has been shown to be effective in forecasting large-scale relational data (Schosser 2022a, 2022b). However, the results so far are limited to point forecasts. The paper at hand gives an indication of the uncertainty in forecasts. It offers prediction intervals that prove to be highly accurate even though they lack a theoretical underpinning. The proposed approach is characterized by exceptionally low computational cost. The latter is increasingly important for the acceptance and successful implementation of a method, especially where thousands, or even millions of forecasts and estimates of uncertainty are needed on a regular basis. Taken together, the findings of this paper provide a compelling argument in favor of tensor extrapolation.