1 INTRODUCTION↩︎

Time series forecasting is essential across many domains where anticipated scenarios guide critical decisions, including energy management, healthcare systems, and financial markets. The task is inherently challenging: series often exhibit complex temporal dynamics, requiring models capable of capturing them.

This paper focuses in particular on forecasting electricity demand, a key challenge for maintaining the balance between supply and demand on the electric grid. As energy production is set according to predicted consumption, reliable forecasts are crucial in keeping the grid operating around the clock. Many state-of-the-art approaches rely on Generalised Additive Models (GAMs), which are a flexible way of representing demand as a sum of smooth functions, applied to pre-processed calendar and weather variables. We refer to [1] for an in-depth review of GAMs and other models for electricity forecasting. However, the selection of smooth functions and covariate pre-processing requires specialist knowledge. For example, due to the thermal inertia of buildings, changes in temperature affect electricity demand with a delay, which can be accounted for by exponential smoothing [2]. Furthermore, the relationship between temperature and demand is nonlinear and can be modelled using splines (see, among others, [3]. We propose a new generic framework to automatically extract these complex dependencies of the target to exogenous variables, regardless of any prior assumptions about the nature of relationships or the impact of the past.

The signature transform, introduced by [4] and further developed in the context of rough path theory by [5], provides a representation of sequential data through iterated integrals referred to as signatures. These objects have been shown to capture key geometric and analytic properties of multi-dimensional paths ([6]). Casting data as continuous paths and computing their signatures has proven to be an effective way of obtaining non-parametric feature sets for time-ordered data in machine learning tasks [7], [8]. Recently, a growing line of work has highlighted the relevance of signatures in modern AI pipelines, ranging from their connection to recurrent neural networks [9] to their integration with transformers for time-series tasks [10].

Building on these ideas, we aim to leverage signatures for discrete-time time series forecasting. To this end, we propose a regression model on sliding-window signatures. Sliding windows are key in our approach as they enable the signatures to focus on the most recent past, discarding distant history as opposed to using exponentially fading memory as in [11]. While [12] already used sliding windows for nowcasting in finance, we propose a novel increment-based model that is both theoretically grounded and straightforward to implement.

The goal of this paper is to present this new forecasting framework that consists of three simple steps: data augmentation, sliding-window signature calculation and ridge regression on the increments of the target time series. We support it by the following main contributions.
Theoretical foundations. We establish universal approximation results for sliding-window signatures in discrete time and prove that these features preserve stationarity when derived from a stationary series, providing a solid theoretical basis for our method.
Efficient algorithm. We develop a sequential algorithm that exploits sliding window overlap for fast and efficient signature calculation.
Empirical validation. Through experiments on synthetic and real-world data, we demonstrate that our method successfully leverages signatures as non-parametric features to capture complex temporal dynamics and outperform baselines with incorporated expert knowledge.

2 SIGNATURES FOR TIME SERIES FORECASTING↩︎

In this section, we present signatures, along with necessary theoretical background, frame them as features and place them in a discrete-time forecasting setting. For a more rigorous mathematical treatment, we refer to the Appendix 7.

2.1 Preliminaries↩︎

We define a continuous path in \(\mathbb{R}^d\) as any continuous mapping \(x : [s, t] \longrightarrow \mathbb{R}^d\). We assume \(d\geq2\).

Intuitively, \(\|x\|_{\mathrm{1}-\mathrm{var}}\) can be seen as the length of the path \(x\). From now on, we consider only paths of finite length, i.e. \(\|x\|_{\mathrm{1}-\mathrm{var}} < \infty\). We denote by \(\otimes\) the tensor product and by \((\mathbb{R}^d)^{\otimes k}\) the \(k\)-th tensor power of \(\mathbb{R}^d\) , with \((\mathbb{R}^d)^{\otimes 0} := \mathbb{R}\). We note that \((\mathbb{R}^d)^{\otimes k}\) is a Hilbert space of dimension \(d^k\).

We endow \(\mathscr{T}(\mathbb{R}^d)\) with the scalar product \(\langle \mathbf{a}, \mathbf{b}\rangle_{\mathscr{T}(\mathbb{R}^d)} = \sum_{k=0}^{\infty} \langle a_k, b_k\rangle_{(\mathbb{R}^d)^{\otimes k}} \,,\) which induces the norm \(\|\mathbf{a}\|_{\mathscr{T}\left(\mathbb{R}^d\right)} = \sqrt{\sum_{k=0}^{\infty}\left\|a_k\right\|_{\left(\mathbb{R}^d\right)^{\otimes k}}^2}\).

2.2 Signatures↩︎

Each tensor \(S^k(x)\) can be written in terms of its elements, referred to as signature coefficients, indexed by multi-indices \((i_1,\ldots,i_k) \in \{1,\ldots,d\}^k\): \[S^{\left(i_1, \ldots, i_k\right)} = \idotsint \displaylimits_{s < u_1<\cdots<u_k < t} \mathrm{d}x^{\left(i_1\right)}_{u_1} \cdots \mathrm{d}x^{\left(i_k\right)}_{u_k}.\]

In practice, instead of dealing with infinite sequences, we only consider the levels up to (truncation) order \(N\), defining the truncated signature as \[S^{\leq N}\left( x_{[s,t]} \right) = \left(1, \;S^1\left( x_{[s,t]} \right), \;\ldots, \;S^N\left( x_{[s,t]} \right) \right).\] Furthermore, we often ignore the tensor structure and view \(S^{\leq N}\) as a collection of all signature coefficients with multi-index of length \(k \leq N\), arranged in a vector of size \(s_d(N)= \sum_{k=0}^{N} d^k = \left(d^{N+1}-1\right) /(d-1)\). We also note that, when the interval is not of importance, we use the abbreviation \(S(x)\).

While the definition may look complex, we now illustrate that in the simple case of linear paths, there is no need for evaluating or numerically approximating the iterated integrals. Instead, we obtain a closed formula for signature coefficients which we could interpret as monomials of path increments.

We now turn to key signature properties leveraged in our framework. It follows from their definition as iterated integrals that signatures are invariant to both translation and time reparametrization (see Proposition [prop:invariances]). Consequently, they do not encode the path’s starting point nor the speed with which it was traversed. A standard approach to reintroduce the latter information is \(\textit{time augmentation}\), i.e. adding time as a path component.

The following proposition and example will prove crucial for efficient signature calculation.

Chen’s identity allows us to extend the result from Example 1 to piecewise linear paths.

The next result ([13] Lemma 5.1) is key for switching from infinite to finite-dimensional objects.

The first equation can be seen as a bound on the information carried by each signature level: the factorial decay of the \(k\)-th level norm means that information decreases as the level increases. Consequently, for sufficiently large \(N\), \(S^{\leq N}(x)\) provides a good approximation of \(S(x)\), justifying truncation in practice. The second equation shows that signatures lie in the separable Hilbert space \(\mathscr{T}(\mathbb{R}^d)\), a property we will exploit when working with probability distributions.

Finally, we present the main motivation for using truncated signatures as features in machine learning tasks.

This result follows directly from the Stone-Weierstrass theorem, showing that signature coefficients serve as the path analogue of monomials. Therefore, they can be leveraged to linearize the problem of learning a more complex (non-linear) function on a compact set of paths.

So far, our treatment of paths and signatures has been deterministic. Turning to the statistical perspective, from now on, we consider random paths \(X : [s, t] \to \mathbb{R}^d\). Indeed, for each \(\omega\) in the sample space \(\Omega\), the realisation \(X(\omega)\) is a path for which we can define the signature \(S\left( x(\omega) \right)\) as before. In other words, \(S(X)\) is a random variable taking values in a separable Hilbert space (Proposition 1).

2.3 Signature Features for Time Series↩︎

We now present the pipeline for using signatures as a feature extraction method for discrete-time time series forecasting.

Let \(\left( Y_t \right)_t\) denote the target time series and \(\left( X_t \right)_t\) denote the covariate time series, with \(Y_t \in \mathbb{R}\) and \(X_t \in \mathbb{R}^d\). We consider predictions for \(Y_t\) of the form \(f(X_{u \le t})\) for some continuous, measurable, possibly non-linear function \(f\). In practice, when forecasting electricity demand, covariates from the distant past, such as temperature from a year ago, have little predictive value. Hence, we assume short-term dependence: the forecast of \(Y_t\) depends only on a finite number of the most recent observations of \(\left( X_t \right)_t\).

To apply the previously established results, discrete data must be first embedded into continuous time to yield a path. The choice of embedding can significantly impact the performance of signature-based methods [14]. We turn to linear interpolation, meaning that the path on the interval \(\left[t, t+1\right]\) is defined by \(u \longmapsto X_{t}+{(u-t)}\left(X_{t+1}-X_{t}\right).\) By concatenating finitely many of these segments, we obtain piecewise linear paths. We note that these are finite-length paths and therefore fall within the theoretical framework established earlier. Moreover, as illustrated in Example 2, they allow for a straightforward iterative computation, consisting in calculating the signatures on each linear segment as in Example 1 and concatenating them inductively through the tensor product (see Remark 7). This procedure is efficiently implemented in the Python package iisignature ([15]), used in this work.

In classical frameworks like  [16], signatures summarize the entire path. In time series, this corresponds to expanding windows, where all past points contribute equally. This is unsuitable for our application, where only recent observations matter. A more natural choice is a sliding window of size \(w\), considering only data from \(t-w\) to \(t\). As shown in Section 3, sliding windows also help preserve stationarity, though they introduce \(w\) as a hyperparameter, which we discuss in Section 5.

Finally, we apply time augmentation within each sliding window: \(\left(0, X_{t-w} \right), \left(1, X_{t-w+1} \right), \ldots, \left( w, X_t \right)\) at time \(t\). This ensures that the resulting paths have at least one monotone coordinate.

Our method consists in the following steps. We first fix a sliding window size \(w\). At each time \(t\), we consider \(\left(0, X_{t-w} \right), \left(1, X_{t-w+1} \right), \ldots, \left( w, X_t \right)\). We linearly interpolate these points to obtain a piecewise linear path, denoted by \(X_{[t-w,t]}\). Having constructed a continuous path, we can calculate its signature \(S\left( X_{[t-w,t]} \right)\) and leverage the theory presented in Section 2. In terms of application, we can now use the truncated signature \(S^{\leq N} \left( X_{[t-w,t]} \right)\) as a feature set for \(\left(X_t\right)_t\) on \([t-w,t]\). We summarize this as a pipeline in Figure 1. The final step of fitting a linear model is motivated by the universal approximation theorem, which has yet to be established for our sliding-window, discrete-time setting. In the following section, we identify conditions under which this result holds and examine the important property of stationarity within our framework.

3 THEORETICAL CONTRIBUTIONS↩︎

We introduce two assumptions on the discrete-time time series of covariates \({(X_t)}_t\).

Assuming the time series is uniformly bounded is not restrictive, as real-world data such as temperature lie within physical limits. Similarly, stationarity of increments is a milder assumption than stationarity of the full series; for example, while temperature exhibits seasonality, its half-hour increments can reasonably be considered stationary (see Appendix 10).

3.1 Universal Approximation Theorem↩︎

We treat the paths on sliding windows as functions \({X}_{[t-w, t]}(\omega) : [0,w] \mapsto \mathbb{R}^{d+1}\), which share a common time domain, allowing us to consider the set of such paths and functions over this set.

Sketch of proof. We first introduce paths on windows translated by the starting point, \({X}_{[t-w, t]} - \left( 0, X_{t-w} \right)\). We consider the set of all realisations of these random paths on all sliding windows. These paths have a monotone coordinate (time augmentation), they all begin at the same point (translation), are uniformly bounded (Assumption 1) and they all are Lipschitz with the same constant (Assumption 1 and piecewise linearity). By constructing a superset of paths with the same properties that is compact in the uniform topology, we can apply Theorem 1 to it, completing the proof. For more detail, see Appendix 8.1.

The main takeaway from this theorem is that we can effectively linearize the task of learning the mapping \(f\) by approximating it with a linear functional on the truncated signature. The result allows for an approximation of a wide range of complicated functions \(f\) (such as GAMs or conditional expectation), without requiring assumptions on their specific form, but given that they are continuous with respect to the input.

3.2 Stationarity of Signatures↩︎

Computing the signature of \((X_t)_t\) over a sliding window produces a new signature series \({\left( S(X_{[t-w,t]}) \right)}_t\). Since these signatures form the features of our forecasting model, we study their stationarity to better characterize the resulting prediction process \((\widehat{Y}_t)_t\).

Sketch of proof. The main argument is that signature coefficients, and consequently truncated signatures, are deterministic, measurable functions of the series’ increments. Applying any measurable, deterministic function to two random vectors with the same distribution preserves their equality in distribution. Therefore, given Assumption 2, the stationarity of truncated signatures holds. The conclusion for signatures follows immediately from the previous result and Remark 3. For full proof, see Appendix 8.2.

4 IMPLEMENTATION↩︎

4.1 Model↩︎

The preceding theoretical considerations yield two key implications for the specification of our model. First, our forecast is a function of the translated paths \(X_{[t-w]} - \left(0,X_{t-w}\right)\). Second, in accordance with Remark 4, under the Assumption 2, our forecast is stationary. It is therefore natural to apply our model on a stationary target variable. These requirements are, however, restrictive and fail to hold in the context of electricity demand. To address this, we introduce an additional hyperparameter: a delay \(D\) in the target variable. The delay \(D\) is selected to be sufficiently large to ensure that at time \(t\) the true value of \(Y\) at time \(t-D\) is observable, but small enough that the increment process \(\Delta_D Y_t := Y_t - Y_{t-D}\) can be treated as stationary. Furthermore, it is far more realistic that the increments \(\Delta_D Y_t\), rather than the raw values \(Y_t\), depend on the translated paths (see Remark 2).

We now present our pipeline applied to \(\Delta_D Y_t\). The training procedure consists in fitting a ridge regression \[\begin{align} \hat{\theta} \in \arg\min\limits_\theta \sum\limits_{t=1}^{t_{train}} \left(\Delta_D Y_t - {\theta}^{\top} S^{\leq N} \left( x_{[t-w,t]} \right)\right)^2 + \lambda \|\theta\|_2^2, \end{align}\] where the regularization constant \(\lambda > 0\) is chosen on the validation set. Our final forecast of the target for any \(t> w\) is then \(\widehat{Y}_t = y_{t-D} + \widehat{\Delta_D Y_t}\), with \(\widehat{\Delta_D Y_t} = {\hat{\theta}}^{\top} S^{\leq N} \left( x_{[t-w,t]} \right)\). See Appendix 9 for the full algorithm in pseudocode.

Several challenges arise when implementing this model. We begin by addressing the choice of hyperparameters: namely \(D, w\) and \(N\). The delay \(D\) depends on the data and is chosen so that the assumptions of stationarity of \(\left( \Delta_D Y_t \right)_t\) and the availability of \(Y_{t-D}\) at time \(t\) hold. When selecting the truncation order \(N\), we highlight that the size of the truncated signature \(s_d(N)\) grows polynomially with the dimension of the covariate series \(d\), but exponentially with \(N\). Consequently, large values of \(N\) (above 10) are typically infeasible. In practice, \(N\) is either chosen on a validation set or fixed to a moderate value. The window size \(w\) can likewise be selected on the validation set, noting that only \(w \geq D\) are considered so that \(\Delta_D Y_t\) may be forecasted as a function of \(X_{[t-w]} - \left(0,X_{t-w}\right)\).

The use of a ridge regression is motivated by the ill-conditioning of the design matrix. High correlations among signature components, amplified by overlapping sliding windows, lead to unstable coefficient estimates. We note that lasso and ridge regression are common when working with signature features  [16], [19]; we adopt ridge, as it performs better in presence of collinearity. Additionally, we omit all signature coefficients that depend solely on time. The \(k\)-th such coefficient is given by a constant \(S^{(1,\ldots,1)}(x_{[t-w,t]}) = \frac{w^k}{k!}\), contributing only redundant intercept terms and worsening the collinearity. We thus redefine \(S^{\leq N}\) to exclude these coefficients.

Finally, the truncated signature’s rapidly growing size (in \(N\) and \(d\)) and the tensor-multiplication involved in its computation make calculating signatures for large and numerous windows computationally expensive. To address this, we propose an algorithm that incrementally updates the signature as the window slides, rather than recomputing it from scratch.

4.2 Efficient Signature Computation↩︎

The key observation is that in our pipeline, there is a large overlap in paths captured within two consecutive sliding windows. To exploit this, we use the algebraic properties of signatures from Proposition 2. Recall that the time-reversal property allows us to remove the contribution of \(X_{t-w}\) by computing the signature of the linear section from \(\left(t-w+1, X_{t-w+1}\right)\) back to \(\left(t-w, X_{t-w}\right)\) (reversing time) and concatenating it on the left of the path using Chen’s identity. Similarly, to incorporate the new data point \(X_{t+1}\), it suffices to compute the signature of the linear section between \(\left(t, X_t\right)\) and \(\left(t+1, X_{t+1}\right)\) and concatenate it on the right, again via Chen’s identity. These sequential updates are possible because the path remains unchanged as the window slides, which is still consistent with the theoretical insights, as highlighted in Remarks 1 and 2.

With the proposed procedure, presented in Algorithm 2 and implemented in Python, the number of tensor products required to compute the signature on a new window drops from \(w\) to just 2. This ensures that, even with a very fine time grid and a large window size, the procedure remains computationally feasible.

5 NUMERICAL EXPERIMENTS↩︎

This section demonstrates the proof of concept of our sliding-window signature approach to time series forecasting. Our focus is on forecasting electricity demand in France, which exhibits a strong nonlinear dependence on recent past values [1]. This time series also depends heavily on its own past values, so it can be useful to consider its available observations as features (see, among other [20]). Our experiments on both synthetic and real data examine the effectiveness of signatures in capturing such dependencies. Both code and data are made open source, see Appendix 10.

5.1 Dataset↩︎

We gather half-hourly electrical demand data from Eco2mix¹ open-source dataset published by RTE, France’s electricity transmission system operator. We combine it with temperature observation data from Météo France ². To obtain a time series of temperatures at the national level, we average the observations across French weather stations. We then linearly interpolate the data to get a unified data set comprising 48 observations per day, covering the period from January 1, 2012, to December 30, 2015. This stable period was chosen in order to avoid the significant fluctuations in electricity demand brought about by the COVID-19 pandemic and the 2022 energy crisis. In what follows, we denote by \(Y_t\) the electricity demand and by \(T_t\) the temperature at any time \(t=1, 2, \dots\).We apply our method for a single feature: \(T_t\), and underline that in practice, it is replaced with temperature forecasts. Furthermore, for any smoothing parameter \(\alpha \in [0,1]\), we can define exponentially smoothed temperature as \[\left\{ \begin{array}{ll} \overline{T}_1^\alpha = T_1 \\ \overline{T}_t^\alpha = (1-\alpha)\overline{T}_{t-1}^\alpha + \alpha T_t^\alpha \,, \,\textrm{ for any } t \geq 2 \,. \end{array} \right.\] As mentioned in the introduction, smoothed temperature values are often used to model thermal inertia and were shown to enhance the quality of forecasts. These variables stem from expert knowledge and will be useful for generating synthetic data and comparing the signature method to relevant baselines.

5.2 Experimental Setup↩︎

In what follows, we represent the data as a two-dimensional path consisting of temperature augmented with rescaled time, i.e., \((t/w, T_t)\), where \(w\) denotes the window size. Time is rescaled so that its increments, and therefore the scale of signature coefficients, remain comparable across different choices of \(w\). The selection of window size is data-driven: since half-hourly demand strongly depends on the time of day, we preserve this structure by considering window sizes in units of full days. For the delay, we fix \(D=2 \text{ days}\), reflecting the availability of observed demand. We use the years 2012 and 2013 to train the model, the year 2014 as a validation set on which we fix the hyperparameters (namely the ridge regularization constant \(\lambda\)) and the year 2015 as a test set.

Given a window size \(w\) and a signature truncation order \(N\), we predict the target time series using the approach described in Subsection 4.1, which we denote by \(\mathrm{RidgeSig}\). We evaluate our method and benchmark models, namely linear regressions performed on various features \(X^1,X^2,\dots\) and denoted \(\mathrm{LR}(X^1,X^2,\dots)\), using the two following standard error metrics: Root Mean Squared Error (RMSE), and Mean Absolute Percentage Error (MAPE) on the test set. We recall that for a target values \(Y_1,\dots Y_T\) and the associated forecasts \(\hat{Y}_1,\dots \hat{Y}_T\), we have: \[\mathrm{RMSE} = \frac{1}{T}\sum_{t=1}^T \big( \hat{Y}_t-Y_t \big)^2, \mathrm{MAPE} = \frac{100}{T}\sum_{t=1}^T\frac{|\hat{Y}_t-Y_t|}{|Y_t|}.\]

5.3 Proof of Concept on Synthetic Data↩︎

5.3.0.1 Data generation.

For a smoothing parameter \(\alpha\) and any time step \(t\), we generate the electricity demand time series \(\big(\widetilde{Y}^\alpha_t\big)_t\) as: \[\widetilde{Y}^\alpha_t = \theta_1 \overline{T}_t^\alpha + \theta_2 \big(\overline{T}_t^\alpha\big)^2 + \varepsilon_t, \quad \varepsilon_t \overset{\mathrm{i.i.d}}{\sim} \mathcal{N}(0,\sigma^2)\,.\] Therefore, it depends non-linearly on recent and present temperatures, with \(\alpha\) controlling the impact of past temperature values - the smaller \(\alpha\), the longer the memory. We refer to Appendix D for details on fitting the model parameters (\(\theta_1\), \(\theta_2\) and \(\sigma\)) and synthetic data fidelity with real demand data.

5.3.0.2 Capturing temporal and non-linear dependencies.

Fixing the smoothing parameter \(\alpha = 0.005\) and the signature truncation order \(N=4\), the optimal window obtained on the test dataset is \(w=9\). Results in Table 1 show that our approach outperforms all the benchmarks, except \(\mathrm{RL}(\overline{T}_t^\alpha, (\overline{T}_t^\alpha)^2)\), which constitutes the best possible result because this model knows exactly what the variables are and how they relate to the target used for data generation. These results suggest that our signature approach captures both temporal dependencies, as it outperforms \(\mathrm{RL}(\overline{T}_t^\alpha)\), and nonlinear effects, as it outperforms \(\mathrm{RL}(T_t, T_t^2)\).

5.3.0.3 Optimal sliding window size.

We now investigate how the optimal window size is affected by the strength of past dependencies in the data. Figure 3 shows the RMSE of our approach with truncation order fixed at \(N=5\), plotted for window lengths ranging from 2 to 32 days, applied to synthetic data generated by two different smoothing parameters \(\alpha = 0.005\) and \(\alpha = 0.05\). This confirms the intuition that with \(\alpha\) decreasing, the optimal window size increases. Capturing this additional structure requires a larger window so that the signature representation has access to a broader history of the covariates.

5.3.0.4 Optimal signature truncation order.

We investigate the impact of the signature truncation order by running our algorithm for \(N=2\) to \(N=7\) and window sizes from \(2\) to \(32\) days. As shown in Figure 4, an order that is too small makes it difficult to detect temporal and non-linear dependencies, which degrades performance Therefore, \(N\) must be large enough. On the other hand, we highlight that the larger the value of N, the longer the computational time and the more unstable the results. This suggests that the truncated order should be between 4 and 6.

5.4 Electricity Demand Forecasting↩︎

Finally, we apply our procedure to real electricity demand data. Since real demand, unlike the synthetic data, exhibits a pronounced weekly cycle, we account for it by changing the forecasting delay to \(D = 7\) days. We add to the previous comparison benchmarks using lagged demand from one week prior \(\left(Y_{t-7\text{ days}}\right)\). We emphasize that we tested various smoothing parameters to tune the benchmarks: \(\alpha = 0.005\) gives the best performance. As suggested by synthetic data, we set the window size to \(w=9\) and the truncated order to \(N=6\). Results of Table 2 demonstrates that our procedure yields the best results across both error metrics. Furthermore, in Figure 5, we compare demand predictions from our model and the best performing baseline from Table 2, focusing on a summer and winter week in 2015. Our model follows more closely observed values, particularly in winter when demand is more sensitive to temperature changes.

6 CONCLUSION↩︎

We established theoretical foundations in a new framework for time series forecasting that relies on sliding-window signatures, proving universality of approximation and stationarity of signatures in the discrete-time setting. We highlighted the importance of aligning paths at a common starting point and including appropriate timestamps as an additional coordinate. Incorporating these insights, we built a model that leverages signatures as non-parametric features for covariate time series. Our sequential algorithm exploits sliding window overlap for efficient signature computation. We demonstrated, both on synthetic and real electricity demand data, that signature features capture nonlinear and temporal dependencies without learned representations or expert knowledge. Our approach does require sensible choices of hyperparameters, particularly window size and target delay, to match the data’s temporal structure and granularity. A limitation lies in the rapid growth of signature size with truncation order and path dimension, though we demonstrated good performance at modest orders (e.g., \(N=6\)). Future directions include extending the framework to discrete-time signature kernels, and adapting our sequential algorithm to an online learning setting.

Acknowledgements↩︎

This research was supported by the European Union’s Horizon Europe research and innovation programme through the Marie Skłodowska-Curie Grant Agreement No.101081674.

The authors would like to thank Prof. Christa Cuchiero for valuable inputs during meetings in Vienna, and Adeline Fermanian for constructive discussions that helped improve this work.

7 THEORETICAL BACKGROUND↩︎

In this section, we present a more rigorous treatment of results concerning paths and path integration, tensor spaces, and signatures. To ensure completeness and readability, some definitions and results from the main body are restated.

7.1 Paths and path integration↩︎

We define a continuous path in \(\mathbb{R}^d\) as any continuous mapping from some interval \([s,t]\) to \(\mathbb{R}^d\), i.e. \(x : [s, t] \longrightarrow \mathbb{R}^d\), \(u \longmapsto x_u = \big(x_u^{(1)}, x_u^{(2)}, \ldots, x_u^{(d)}\big)\). We consider multi-dimensional paths, meaning \(d \geq 2\). Moreover, we restrict ourselves to working with paths of finite length, which we define through the notion of the \(1\)-variation.

Intuitively, \(\|x\|_{\mathrm{1}-\mathrm{var}}\) can be interpreted as the length of the path \(x\), therefore we work with paths such that \(\|x\|_{\mathrm{1}-\mathrm{var}} < \infty\). These paths are also referred to as being of bounded variation. This assumption will later prove important for theoretical guarantees, and it also allows us to define path integration in the \(\textit{Riemann-Stieltjes}\) sense. We refer to Chapters 2 and 3 in [6] for more details.

This definition can be generalized to the case when \(x\) and \(y\) are vector-valued; for \(x, y\) : \([s,t] \rightarrow \mathbb{R}^d\) we have \[\int_s^t y_u \mathrm{~d} x_u=\left(\begin{array}{c} \int_s^t y_u^{(1)} \mathrm{~d} x_u^{(1)} \\ \vdots \\ \int_s^t y_u^{(d)} \mathrm{~d} x_u^{(d)} \end{array}\right).\]

7.2 Tensor spaces↩︎

We denote by \(\otimes\) the tensor product and by \((\mathbb{R}^d)^{\otimes k}\) the \(k\)-th tensor power of \(\mathbb{R}^d\) , with \((\mathbb{R}^d)^{\otimes 0} := \mathbb{R}\). Let \(e_1, \ldots, e_d\) be the canonical basis of \(\mathbb{R}^d\), then \(\left(e_{i_1} \otimes \cdots \otimes e_{i_k}\right)_{1 \leq i_1, \ldots, i_k \leq d}\) is an orthonormal basis of \(\left(\mathbb{R}^d\right)^{\otimes k}\), meaning that any element \(a \in\left(\mathbb{R}^d\right)^{\otimes k}\) can be written as \(a=\sum_{1 \leq i_1, \ldots, i_k \leq d} a^{\left(i_1, \ldots, i_k\right)} e_{i_1} \otimes \cdots \otimes e_{i_k}\), where \(a^{\left(i_1, \ldots, i_k\right)} \in \mathbb{R}\).

Furthermore, \(\left(\mathbb{R}^d\right)^{\otimes k}\) is a Hilbert space of dimension \(d^k\), with the following scalar product and norm: \[\langle a, b\rangle_{\left(\mathbb{R}^d\right)^{\otimes k}}=\sum_{1 \leq i_1, \ldots, i_k \leq d} a^{\left(i_1, \ldots, i_k\right)} b^{\left(i_1, \ldots, i_k\right)}, \quad \|a\|_{\left(\mathbb{R}^d\right)^{\otimes k}} = \sqrt{\sum_{1 \leq i_1, \ldots, i_k \leq d} {\left( a^{\left(i_1, \ldots, i_k\right)} \right)} ^2} \quad .\]

In other words, the extended tensor algebra is a set of infinite sequences of tensors of increasing order: \(a_0\) is a scalar, \(a_1\) a vector, \(a_2\) a matrix, \(a_3\) a “cube" (a third-order tensor), and so on. It can be shown that \(T\left(\left(\mathbb{R}^d\right)\right)\) is a non-commutative algebra under the tensor product \(\otimes\), with the neutral element \(\mathbf{1}=(1,0,0,\ldots)\). Furthermore, it holds that the subset \(T_1\left(\left(\mathbb{R}^n\right)\right):=\left\{\mathbf{a} \mid \mathbf{a} \in T\left(\left(\mathbb{R}^n\right)\right) \text{ with } a_0=1\right\}\) is a Lie group. Finally, the following subspace will be of special interest.

We can endow this space with the following scalar product and associated norm: for any \(\mathbf{a}, \mathbf{b} \in \mathscr{T}\left(\mathbb{R}^d\right)\), \[\langle \mathbf{a}, \mathbf{b}\rangle_{\mathscr{T}\left(\mathbb{R}^d\right)}=\sum_{k=0}^{\infty} \langle a_k, b_k\rangle_{\left(\mathbb{R}^d\right)^{\otimes k}}, \quad \|\mathbf{a}\|_{\mathscr{T}\left(\mathbb{R}^d\right)} = \sqrt{\sum_{k=0}^{\infty}\left\|a_k\right\|_{\left(\mathbb{R}^d\right)^{\otimes k}}^2} . \label{normhilbert}\tag{1}\]

7.3 Signatures↩︎

We can now fully define signatures, in the light of the theory presented earlier.

7.3.1 Definitions↩︎

Using the orthonormal basis for \({\left(\mathbb{R}^d\right)}^{\otimes k}\), we can write \[{S^k \left( x_{[s,t]} \right)} = \sum_{\left(i_1, \ldots, i_k\right) \subset\{1, \ldots, d\}^k} S^{\left(i_1, \ldots, i_k\right)}\left(x_{[s,t]}\right) e_{i_1} \otimes \cdots \otimes e_{i_k},\] where \(S^{\left(i_1, \ldots, i_k\right)}\left(x_{[s,t]}\right)\) is referred to as the signature coefficient of \(x\) along the multi-index \(\left(i_1, \ldots, i_k\right) \subset\{1, \ldots, d\}^k, k \geq 1\) on \([s,t]\) and is given by \[S^{\left(i_1, \ldots, i_k\right)} = \idotsint \displaylimits_{s < u_1<\cdots<u_k < t} \mathrm{d}x^{\left(i_1\right)}_{u_1} \cdots \mathrm{d}x^{\left(i_k\right)}_{u_k}.\]

In practice, we don’t work with infinite sequences of tensors. Firstly, we only consider a finite number of signature levels, defining the truncated signature as \[S^{\leq N}\left( x_{[s,t]} \right) = \left(1, \;S^1\left( x_{[s,t]} \right), \;\ldots, \;S^N\left( x_{[s,t]} \right) \right),\] where we refer to \(N\) as the truncation order. Secondly, in applications, we disregard the underlying tensor structure, and consider the truncated signature as a a collection of all signature coefficients with multi-index of length \(k \leq N\), arranged in a vector \[\left(1, S^{(1)}\left( x_{[s,t]} \right), \ldots, S^{(d)}\left( x_{[s,t]} \right), S^{(1,1)}\left( x_{[s,t]} \right), S^{(1,2)}\left( x_{[s,t]} \right), \ldots, S^{\overbrace{ \left( d,d, \ldots, d \right)}^{N \text{times}}} \left( x_{[s,t]} \right) \right),\] that is of size \(s_d(N)= \sum_{k=0}^{N} d^k = \left(d^{N+1}-1\right) /(d-1)\).

We now illustrate the definitions on the example of a linear path, deriving a closed formula for signature coefficients in this simple case. This will prove an important building block in our methodology.

7.3.2 Properties↩︎

Both of these properties can be derived straight from the definition of signatures as iterated integrals. Invariance to reparametrization is a consequence of the change of variable formula in Riemann-Stieltjes integration, while invariance to translation follows from the fact that \(d\overline{x}_t = dx_t\).

We can rephrase Chen’s relation in terms of signature coefficients: for any multi-index \(\left(i_1, \ldots, i_k\right) \subset\{1, \ldots, d\}^k\), it holds that \[S^{(i_1, \ldots, i_k)}\left((x * y)_{[s,u]}\right) = \sum_{\ell=0}^k S^{(i_1, \ldots, i_\ell)}\left(x_{[s,t]}\right) \cdot S^{(i_{\ell+1}, \ldots, i_k)}\left(y_{[s,t]}\right) \label{chencomp}.\tag{2}\]

This result is crucial for the operationalization of signature computation as it provides a recursive formula for the signature of a concatenation of paths. As noted earlier, it also enables an easy calculation in the case of piecewise linear paths, which we explain in greater detail in the following remark.

Computing the signature of a piecewise linear path therefore reduces to applying the iterative formula stated above, requiring no numerical integration. This procedure is efficiently implemented in the Python package iisignature ([15]).

Going back to the time reversal property, we remark that, by definition, the first element of the signature is set to 1 and therefore \(S\left( x_{[s,t]} \right) \in T_1\left(\left(\mathbb{R}^d\right)\right)\). This result now states that the inverse of the signature of a path \(x\) in the tensor group \(T_1\left(\left(\mathbb{R}^d\right)\right)\) is exactly the signature of \(x\) traversed backwards in time.

8 PROOFS OF THEORETICAL CONTRIBUTIONS↩︎

8.1 Proof of Theorem 2↩︎

8.2 Proof of Theorem 3↩︎

We present a more formal statement of the result with the corresponding proof.

9 IMPLEMENTATION DETAILS↩︎

Our method consists in efficiently computing signatures on sliding windows and fitting a ridge regression model, using the covariate signatures as predictors and the increments of the target series as the response variable. Hyperparameters are tuned on a validation set. The full procedure is detailed in the pseudocode below.

Finally, to obtain the forecast \(\widehat{Y}_t\) for \(t\) outside the train-validation set, it is enough to update the signature as in Algorithm 2 and use the estimated ridge regression coefficients \(\widehat{\theta}\) to compute the forecast \(\widehat{Y}_t = y_{t-D} + \widehat{\theta}^\top S(x_{[t-{w},t]})\).

We note that this sequential procedure can easily be transferred to an online learning setting.

10 EXPERIMENTAL DETAILS↩︎

Technical details. The data and code used in this paper are publicly available in our GitHub repository: https://github.com/ninadrobac/slidesig. All experiments were conducted on a standard laptop (MacBook Air, Apple M2 chip, 16 GB RAM), demonstrating that our implementation can be executed efficiently on standard hardware.

Real data. We use half-hourly electricity demand data from RTE (France’s electricity transmission system operator), aggregated on the national level, and temperature observations from Météo France recorded on a three-hour grid. The temperature series on the national level is obtained by averaging across weather stations and linearly interpolating to match the demand frequency. This yields a unified dataset with 48 observations per day from January 1, 2012, to December 30, 2015, resulting in 70128 data points of temperature and demand. This period is selected to avoid the demand fluctuations associated with the COVID-19 pandemic and the 2022 energy crisis. We denote by \(Y_t\) the electricity demand and by \(T_t\) the temperature at any time \(t \geq 1\).

Synthetic data. Electricity consumption depends on many factors, including calendar variables such as weekdays and public holidays. Since our goal is to isolate the role of temperature, we construct a synthetic dataset in which both the strength of past dependency and the functional form of the relationship can be controlled. To model memory effects, we introduce exponentially smoothed temperature as \[\left\{ \begin{array}{ll} \overline{T}_1^\alpha = T_1 \\ \overline{T}_t^\alpha = (1-\alpha)\overline{T}_{t-1}^\alpha + \alpha T_t^\alpha \,, \,\textrm{ for any } t \geq 2 \,, \end{array} \right.\] where the smoothing parameter \(\alpha \in (0,1]\) controls the influence of past values - the smaller \(\alpha\), the stronger the influence of past values. We then fit a linear regression model with observed consumption \(Y_t\) as the target variable and \(\overline{T}_t\) and \(\overline{T}_t^2\) as covariates: \[\widehat{Y}^\alpha_t = \widehat{\theta_1} \overline{T}_t^\alpha + \widehat{\theta_2} \big(\overline{T}_t^\alpha\big)^2, \label{eq:lin95model95synthetic}\tag{3}\]

thereby imposing a quadratic dependence of simulated demand on smoothed temperature. The final synthetic demand time series \(\big(\widetilde{Y}^\alpha_t\big)_t\) is obtained by adding normally distributed noise to the fitted values \(\widehat{Y}^\alpha_t\) : \[\widetilde{Y}^\alpha_t = \widehat{\theta_1} \overline{T}_t^\alpha + \widehat{\theta_2} \big(\overline{T}_t^\alpha\big)^2 + \varepsilon_t, \quad \varepsilon_t \overset{\mathrm{i.i.d}}{\sim} \mathcal{N}(0,\sigma^2)\,. \label{eq:model95synthetic}\tag{4}\]

In all experiments presented, we fix \(\sigma\) to introduce a moderate level of noise (\(\sigma=1000\) for the first and \(\sigma=500\) for the second and third synthetic data experiment) and we set a fixed random seed when simulating the noise to ensure reproducibility. The synthetic demand series closely follows the main patterns of the observed series, effectively capturing the yearly variations driven by temperature, providing a realistic proxy for evaluation. However, as illustrated in Figure 7, Model 4 does not reproduce finer weekly patterns that depend on calendar variables such as the day of the week or time of day, as it only contains temperature information. This is reflected in the fact that for experiments on synthetic data we fix \(D=2\), and for experiments on real data we choose \(D=7\) to account for weekly patterns. Although we could have adapted Model 4 by adding daily seasonality to better match the data locally, this is simply a proof of concept, so we opted for the simplest possible framework.

We note that by varying the smoothing parameter \(\alpha\), we can generate different versions of the dataset in which past values exert more or less influence on the simulated demand. Over several tested parameters, we mainly focus on \(\alpha = 0.005\) as it offers the best fit of Model 4 to the real data.

Assumptions in practice. As discussed earlier, Assumption 1 is satisfied since temperature values in the dataset remain within natural bounds (between \(-10\,^{\circ}\mathrm{C}\) and \(35\,^{\circ}\mathrm{C}\)). In contrast, verifying Assumption 2 is less straightforward, as it concerns the distributional properties of temperature increments. To assess this, we apply two standard statistical tests — the Augmented Dickey–Fuller (ADF) and Kwiatkowski–Phillips–Schmidt–Shin (KPSS) tests, both of which indicate that the time series of half-hourly temperature increments is weakly stationary at a significance level of \(0.05\). The same analysis was performed for weekly demand increments and two-day increments of simulated demand, yielding satisfactory results.

References↩︎

---

1. https://www.rte-france.com/eco2mix↩︎

2. https://donneespubliques.meteofrance.fr↩︎