<figure><img src="image_1.png" /></figure>

HAL
open science

# Sliding-Window Signatures for Time Series: Application to Electricity Demand Forecasting

Nina Drobac, Margaux Brégère, Joseph de Vilmarest, Olivier Wintenberger

<figure><img src="image_2.png" /></figure>

## To cite this version:

Nina Drobac, Margaux Brégère, Joseph de Vilmarest, Olivier Wintenberger. Sliding-Window Signatures for Time Series: Application to Electricity Demand Forecasting. 2025. hal-05308724

HAL Id: hal-05308724

https://hal.science/hal-05308724v1

Preprint submitted on 13 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sliding-Window Signatures for Time Series: Application to Electricity Demand Forecasting

Nina Drobac

Sorbonne Université

Margaux Brégère

Sorbonne Université
EDF R&D

Joseph de Vilmarest

Viking Conseil

Olivier Wintenberger

Sorbonne Université

## Abstract

Nonlinear and delayed effects of covariates often render time series forecasting challenging. To this end, we propose a novel forecasting framework based on ridge regression with signature features calculated on sliding windows. These features capture complex temporal dynamics without relying on learned or hand-crafted representations. Focusing on the discrete-time setting, we establish theoretical guarantees, namely universality of approximation and stationarity of signatures. We introduce an efficient sequential algorithm for computing signatures on sliding windows. The method is evaluated on both synthetic and real electricity demand data. Results show that signature features effectively encode temporal and nonlinear dependencies, yielding accurate forecasts competitive with those based on expert knowledge.

## 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 Antoniadis et al. (2024) 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 (Goude et al., 2013). Furthermore, the relationship between temperature and demand is nonlinear and can be modelled using splines (see, among others, Nedellec et al., 2014. 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 Chen (1958) and further developed in the context of rough path theory by Lyons et al. (2007), 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 multidimensional paths (Friz and Victoir, 2010). 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 (Chevyrev and Kormilitzin, 2025; Fermanian, 2021b). Recently, a growing line of work has highlighted the relevance of signatures in modern AI pipelines, ranging from their connection to recurrent neural networks (Fermanian et al., 2021) to their integration with transformers for time-series tasks (Moreno-Pino et al., 2025).

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 Jaber and Sotnikov (2025). While Cohen et al. (2023) 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 A.

## 2.1 Preliminaries

We define a continuous path in $\mathbb{R}^d$ as any continuous mapping $x : [s, t] \to \mathbb{R}^d$. We assume $d \ge 2$.

**Definition 1.** Let $x : [s, t] \to \mathbb{R}^d$ be a continuous path. The 1-variation of $x$ is defined by

$$ \|x\|_{1-\text{var}} = \sup_{(t_0, \dots, t_k) \in P} \sum_{i=1}^{k} \|x_{t_i} - x_{t_{i-1}}\|, $$

where $P = \{(t_0, \dots, t_k) \mid k \ge 0, s = t_0 < \dots < t_k = t\}$ denotes the set of all finite partitions of $[s, t]$.

Intuitively, $\|x\|_{1-\text{var}}$ can be seen as the length of the path $x$. From now on, we consider only paths of finite length, i.e. $\|x\|_{1-\text{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$.

**Definition 2.** We denote by $\mathcal{T}$ the space of square summable sequences of tensors of increasing order:

$$ \mathcal{T}(\mathbb{R}^d) = \left\{ (a_k)_{k \ge 0} \mid a_k \in (\mathbb{R}^d)^{\otimes k}, \sum_{k=0}^{\infty} \|a_k\|_{(\mathbb{R}^d)^{\otimes k}}^2 < \infty \right\}. $$

We endow $\mathcal{T}(\mathbb{R}^d)$ with the scalar product $\langle \mathbf{a}, \mathbf{b} \rangle_{\mathcal{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}\|_{\mathcal{T}(\mathbb{R}^d)} = \sqrt{\sum_{k=0}^{\infty} \|a_k\|_{(\mathbb{R}^d)^{\otimes k}}^2}$.

**Proposition 1.** $(\mathcal{T}(\mathbb{R}^d), \langle \cdot, \cdot \rangle_{\mathcal{T}(\mathbb{R}^d)})$ is a separable Hilbert space.

## 2.2 Signatures

**Definition 3.** The signature of a finite-length path $x$ on $[s, t]$ is defined as an infinite tensor sequence:

$$ S(x_{[s,t]}) = (1, S^1(x_{[s,t]}), \dots, S^k(x_{[s,t]}), \dots), $$

where the k-th element (called level) is given by

$$ S^k(x_{[s,t]}) = \int \dots \int dx_{u_1} \otimes \dots \otimes dx_{u_k} \in (\mathbb{R}^d)^{\otimes k}. $$

Each tensor $S^k(x)$ can be written in terms of its elements, referred to as signature coefficients, indexed by multi-indices $(i_1, \dots, i_k) \in \{1, \dots, d\}^k$:

$$ S^{(i_1, \dots, i_k)} = \int \dots \int dx_{u_1}^{(i_1)} \dots dx_{u_k}^{(i_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^{\le N}(x_{[s,t]}) = (1, S^1(x_{[s,t]}), \dots, S^N(x_{[s,t]})). $$

Furthermore, we often ignore the tensor structure and view $S^{\le N}$ as a collection of all signature coefficients with multi-index of length $k \le N$, arranged in a vector of size $s_d(N) = \sum_{k=0}^{N} d^k = (d^{N+1} - 1)/(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.

**Example 1.** (Signature of a linear path). Let $x : [s, t] \mapsto \mathbb{R}^d$, $u \to \frac{x_t - x_s}{t-s}(u-s) + x_s$ be a $d$-dimensional linear path. It follows by integration that $S^{(i_1, \dots, i_k)}(x_{[s,t]}) = \frac{1}{k!} \prod_{j=1}^k (x_t^{(i_j)} - x_s^{(i_j)})$ and $S^k(x_{[s,t]}) = \frac{1}{k!} (x_t - x_s)^{\otimes k}$. (Detailed calculation in Example 4.)

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 4). 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 time augmentation, i.e. adding time as a path component.

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

**Proposition 2 (Algebraic properties).** Let $x : [s, t] \mapsto \mathbb{R}^d$ and $y : [t, u] \mapsto \mathbb{R}^d$ denote two paths of finite length.

1. (Chen's identity) Let $x * y : [s, u] \mapsto \mathbb{R}^d$ be the concatenation of $x$ and $y$, meaning $(x * y)_v = x_v$ for $v \in [s, t]$ and $(x*y)_v = x_t + y_v - y_t$ for $v \in [t, u]$. Then $S((x * y)_{[s,u]}) = S(x_{[s,t]}) \otimes S(y_{[t,u]}).$

2. (Time reversal) We denote the time-reversal of $x$ as the path $\overleftarrow{x} : [s, t] \mapsto \mathbb{R}^d$ where $\overleftarrow{x}(u) = x_{s+t-u}$. Then $S(x_{[s,t]}) \otimes S(\overleftarrow{x}_{[s,t]}) = (1, 0, 0, \dots)$.

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

**Example 2.** Let $x : [s, t] \to \mathbb{R}^d$ be a piecewise linear path and let $s = u_0 < u_1 < \dots < u_k = t$ be a partition such that $x$ is linear on each $[u_{k-1}, u_k]$. From Chen's identity we have:

$$S(x_{[s,t]}) = S(x_{[s,u_1]}) \otimes \dots \otimes S(x_{[u_{k-1}, t]}),$$
where every $S(x_{[u_{k-1}, u_k]})$ is given in Example 1.

The next result (Lyons, 2014, Lemma 5.1) is key for switching from infinite to finite-dimensional objects.

**Proposition 3.** Let $x : [s, t] \to \mathbb{R}^d$ be a path of finite length. It holds that:

$$\|S^k(x_{[s,t]})\|_{(\mathbb{R}^d)^{\otimes k}} \le \frac{\|x\|_{1-var}^k}{k!},$$
$$\|S(x_{[s,t]})\|_{\mathcal{T}(\mathbb{R}^d)} \le \exp(\|x\|_{1-var}) < \infty.$$

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^{\le 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 $\mathcal{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.

**Theorem 1 (Universal approximation).** Let $K$ be a compact subset of the space of finite-length paths from $[s, t]$ to $\mathbb{R}^d$ and such that for any $x \in K$, $x_s = a$ for

some $a \in \mathbb{R}^d$ and $x$ has at least one monotone coordinate. Let $f : K \to \mathbb{R}$ be continuous. Then, for every $\varepsilon > 0$, there exists $N \ge 1$, $\theta \in \mathbb{R}^{s_d(N)}$, such that, for any $x \in K$,

$$|f(x) - \theta^\top S^{\le N}(x_{[s,t]})| \le \varepsilon.$$

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(x(\omega))$ 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 $(Y_t)_t$ denote the target time series and $(X_t)_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 $(X_t)_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 (Fermanian, 2021a). We turn to linear interpolation, meaning that the path on the interval $[t, t+1]$ is defined by $u \mapsto X_t + (u-t)(X_{t+1} - X_t)$. 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** (Reizenstein and Graham, 2020), used
in this work.

In classical frameworks like Fermanian (2022), sig-
natures 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, con-
sidering only data from t − w to t. As shown in Sec-
tion 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 slid-
ing window: (0, X<sub>t-w</sub>), (1, X<sub>t-w+1</sub>), . . . , (w, X<sub>t</sub>) at
time t. This ensures that the resulting paths have at
least one monotone coordinate.

**Remark 1.** This approach yields the same signature values as if the time coordinate was added over the entire time horizon: $(t-w, X_{t-w}), (t-w+1, X_{t-w+1}), \dots, (t, X_t)$. Consequently, when the window slides, only the new observation needs to be added, without adjusting previous timestamps.

Our method consists in the following steps. We first fix a sliding window size $w$. At each time $t$, we consider $(0, X_{t-w})$, $(1, X_{t-w+1})$, ..., $(w, X_t)$. 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(X_{[t-w,t]})$ and leverage the theory presented in Section 2. In terms of application, we can now use the truncated signature $S^{\le N}(X_{[t-w,t]})$ as a feature set for $(X_t)_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$.

**Assumption 1 (Uniform boundedness).** Suppose that
$|X_t^{(i)}| \le M \text{ a.s. for some } M > 0 \text{ and all } 1 \le i \le d$.

**Assumption 2 (Stationarity of increments).** Suppose that the increments $(X_{t+1}-X_t)_{t \ge 1}$ are strictly stationary.

Assuming the time series is uniformly bounded is not

restrictive, as real-world data such as temperature lie
within physical limits. Similarly, stationarity of in-
crements 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 D).

**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.

**Theorem 2** (Universal approximation on sliding windows). Let $(X_t)_t$ be a $d$-dimensional time series such that Assumption 1 holds. Let $w \in \mathbb{N}$ be a fixed window size and $f: C([0, w], \mathbb{R}^{d+1}) \mapsto \mathbb{R}$ a continuous function given the uniform topology. It holds that for every $\varepsilon > 0$, there exists $N \in \mathbb{N}$, $\theta \in \mathbb{R}^{s_{d+1}(N)}$, such that

$$
\sup_t |f(X_{[t-w,t]} - (0, X_{t-w})) - \theta^\top S^{\le N} (X_{[t-w,t]})| \le \varepsilon.
$$

Sketch of proof. We first introduce paths on windows translated by the starting point, $X_{[t-w,t]} - (0, X_{t-w})$. 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 B.1.

**Remark 2.** *Assuming all paths start at the same point is vital, since signatures cannot encode the path’s starting value due to translation invariance. This is commonly addressed in practice through basepoint augmentation, prepending a zero to each path (Morrill et al., 2021). As it was not suitable for our sliding-window framework, we instead subtract the initial value from each window. The truncated signature remains unchanged under this transformation, $S^{\le N}(X_{[t-w,t]} - (0, X_{t-w})) = S^{\le N}(X_{[t-w,t]})$, again due to translation invariance. This is key for efficient sequential computation, as there is no need to modify the path at each step. The subtraction has mainly a theoretical implication: we consider the function $f$ as acting on translated paths, though in practice the starting point still carries predictive information. We address this explicitly in Section 4.*

The main takeaway from this theorem is that we can
effectively linearize the task of learning the mapping
<figure><img src="image_4.png" /><figcaption>Figure 1: Method Pipeline: leveraging sliding-window signature features to linearize the forecasting task.</figcaption></figure>

$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 $(S(X_{[t-w,t]}))_t$. Since these signatures form the features of our forecasting model, we study their stationarity to better characterize the resulting prediction process $(\hat{Y}_t)_t$.

**Remark 3.** Although signatures are defined as infinite sequences, the usual definition of strict stationarity (via equality in distribution) still applies, since they lie in the separable Hilbert space $\mathcal{T}(\mathbb{R}^d)$, where finite-dimensional marginal distributions determine the law of the process (see, e.g. Chapter 1 of Bosq, 2000).

**Theorem 3.** Let $(X_t)_t$ be a discrete $d$-dimensional time series such that Assumption 2 holds. Then, for any truncation order $N \in \mathbb{N}$, the time series $(S^{\le N}(X_{[t-w,t]}))_t$ is strictly stationary. Furthermore, the time series $(S(X_{[t-w,t]}))_t$ is strictly stationary.

*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 B.2.

**Remark 4.** Let us now suppose that Assumptions 1, 2 hold for our time series of covariates $(X_t)_t$. Combining the two theoretical results, we can now conclude that the prediction process $(\hat{Y}_t)_t$, with $\hat{Y}_t = \theta^\top S^{\le N} (X_{[t-w,t]})$, is also strictly stationary as a linear transformation of a stationary process. This has important implications for choosing a model.

# 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]} - (0, X_{t-w})$. 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

$$ \hat{\theta} \in \arg \min_{\theta} \sum_{t=1}^{t_{\text{train}}} (\Delta_D Y_t - \theta^{\top} S^{\le N} (x_{[t-w,t]}))^2 + \lambda \|\theta\|_2^2, $$

where the regularization constant $\lambda > 0$ is chosen on the validation set. Our final forecast of the target for any $t > w$ is then $\hat{Y}_t = y_{t-D} + \hat{\Delta}_D \hat{Y}_t$, with $\hat{\Delta}_D \hat{Y}_t = \hat{\theta}^{\top} S^{\le N} (x_{[t-w,t]})$. See Appendix C 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 $(\Delta_D Y_t)_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 \ge D$ are considered so that $\Delta_D Y_t$ may be forecasted as a function of $X_{[t-w]} - (0, X_{t-w})$.
The use of a ridge regression is motivated by the ill-
conditioning of the design matrix. High correlations
among signature components, amplified by overlap-
ping sliding windows, lead to unstable coefficient es-
timates. We note that lasso and ridge regression are
common when working with signature features (Guo
et al., 2024; Fermanian, 2022); we adopt ridge, as
it performs better in presence of collinearity. Addi-
tionally, we omit all signature coefficients that depend
solely on time. The k-th such coefficient is given by
a constant $S^{(1,...,1)}(x_{[t-w,t]}) = \frac{w^k}{k!}$, contributing only
redundant intercept terms and worsening the collinear-
ity. We thus redefine $S^{\le 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 incremen-
tally 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 $(t-w+1, X_{t-w+1})$ back to
$(t-w, X_{t-w})$ (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
$(t, X_t)$ and $(t+1, X_{t+1})$ and concatenate it on the
right, again via Chen's identity. These sequential up-
dates 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 1 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 fore-
casting. Our focus is on forecasting electricity demand
in France, which exhibits a strong nonlinear depen-
dence on recent past values (see, e.g. Section 1.2.4
of Antoniadis et al., 2024). This time series also de-

**Algorithm 1:** Signature update on sliding windows

**Input** : observations $\{x_t | t \ge 0\}$, window size $w$, truncation order $N$

**Output:** signatures $S^{\le N}(x_{[t-w,t]})$ for all $t \ge w$
Compute initial signature $S \leftarrow S^{\le N}(x_{[0,w]})$;

1. **for** $t > w$ **do**

2. Compute signature of the oldest segment
   $S_{old} \leftarrow S^{\le N}(x_{[t-w-1,t-w]})$;

3. Update $S$ by overwriting the oldest segment:
   $S \leftarrow S_{old} \otimes S$; // time-reversal update

4. Compute signature of the new segment
   $S_{new} \leftarrow S^{\le N}(x_{[t-1,t]})$;

5. Update $S$ by adding the new segment
   $S \leftarrow S \otimes S_{new}$; // Chen's identity

pends heavily on its own past values, so it can be useful
to consider its available observations as features (see,
among other De Vilmarest and Goude, 2022). Our ex-
periments on both synthetic and real data examine the
effectiveness of signatures in capturing such dependen-
cies. Both code and data are made open source, see
Appendix D.

## 5.1 Dataset

We gather half-hourly electrical demand data from Eco2mix<fnref n="1" /> open-source dataset published by RTE, France's electricity transmission system operator. We combine it with temperature observation data from Météo France <fnref n="2" />. 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

$$ \begin{cases} \bar{T}_1^\alpha = T_1 \\ \bar{T}_t^\alpha = (1-\alpha)\bar{T}_{t-1}^\alpha + \alpha T_t^\alpha, & \text{for any } t \ge 2. \end{cases} $$

As mentioned in the introduction, smoothed temper-
ature values are often used to model thermal inertia
and were shown to enhance the quality of forecasts.

<footnote n="1">¹https://www.rte-france.com/eco2mix</footnote>

<footnote n="2">²https://donneespubliques.meteofrance.fr</footnote>
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$ 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 RidgeSig. We evaluate our method and benchmark models, namely linear regressions performed on various features $X^1, X^2, \dots$ and denoted $\text{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:

$$ \text{RMSE} = \frac{1}{T} \sum_{t=1}^{T} (\hat{Y}_t - Y_t)^2, \text{MAPE} = \frac{100}{T} \sum_{t=1}^{T} \frac{|\hat{Y}_t - Y_t|}{|Y_t|} $$

## 5.3 Proof of Concept on Synthetic Data

**Data generation.** For a smoothing parameter $\alpha$ and any time step $t$, we generate the electricity demand time series $(\tilde{Y}_t^\alpha)_t$ as:

$$ \tilde{Y}_t^\alpha = \theta_1 \bar{T}_t^\alpha + \theta_2 (\bar{T}_t^\alpha)^2 + \varepsilon_t, \quad \varepsilon_t \stackrel{\text{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.

**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 shows that our approach outperforms all the benchmarks, except $\text{RL}(\bar{T}_t^\alpha, (\bar{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 $\text{RL}(\bar{T}_t^\alpha)$, and nonlinear effects, as it outperforms $\text{RL}(T_t, T_t^2)$.

<table><caption>Table 1: RMSE (MW) and MAPE (%) on test data set for synthetic data generated with $\alpha = 0.005$.</caption><thead><tr><th>Model</th><th>LR(T)</th><th>LR(T, T<sup>2</sup>)</th><th>LR(T<sup>α</sup>)</th><th>RidgeSig</th><th>LR(T<sup>α</sup>, (T<sub>t</sub><sup>α</sup>)<sup>2</sup>)</th></tr></thead><tbody><tr><td>RMSE (MW)</td><td>5.435</td><td>4.729</td><td>3.079</td><td>1.637</td><td>995</td></tr><tr><td>MAPE (%)</td><td>8.3</td><td>6.4</td><td>4.9</td><td>2.5</td><td>1.5</td></tr></tbody></table>

**Optimal sliding window size.** We now investigate how the optimal window size is affected by the strength of past dependencies in the data. Figure 2 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.

<figure><img src="image_5.png" /><figcaption>Figure 2: RMSE (MW) on the test set as a function of the sliding window size for synthetic data generated with $\alpha = 0.005$ (top: optimal window = 9 days) and $\alpha = 0.05$ (bottom: optimal window = 3 days).</figcaption></figure>

**Remark 5.** We can show by induction that $\bar{T}_t^\alpha = \sum_{s=0}^{t-1} \alpha(1-\alpha)^s T_{t-s} + (1-\alpha)^t T_1$, so the ratio between the proportion of $T_{t-s}$ and the one of $T_t$ in the smoothing equals $(1-\alpha)^s$. After 3 days, this ratio still equals 0.48 for $\alpha = 0.005$, while it is already $6e^{-4}$ for $\alpha = 0.05$. This suggests that the temperature from
<figure><img src="image_6.png" /><figcaption>Figure 3: RMSE (MW) on the test set for synthetic data ($\alpha$=0.005) as a function of the sliding window size and for different truncation orders.</figcaption></figure>

three days ago continues to have a significant impact on the data generated by the first model, whereas it is already being forgotten by the second model. Figure 2 perfectly illustrates this phenomenon.

**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 3, 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 ($Y_{t-7}$ days). 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 4, 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.

<table><caption>Table 2: Test RMSE and MAPE on real demand.</caption><thead><tr><th>Model</th><th>RMSE (MW)</th><th>MAPE (%)</th></tr></thead><tbody><tr><td>LR (T&#773;, T&#773;<sup>2</sup>)</td><td>6518</td><td>10.7</td></tr><tr><td>LR (T, T<sup>2</sup>, T&#773;, T&#773;<sup>2</sup>)</td><td>6172</td><td>9.9</td></tr><tr><td>LR (Y<sub>t-7</sub> days)</td><td>4149</td><td>5.3</td></tr><tr><td>LR (T, T<sup>2</sup>, T&#773;, T&#773;<sup>2</sup>, Y<sub>t-7</sub> days)</td><td>3714</td><td>5.3</td></tr><tr><td>RidgeSig</td><td>3150</td><td>4.4</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 4: Observed and predicted demand for a winter (top) and summer (bottom) week in the test set.</figcaption></figure>

# 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 non-linear 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.

## References

Antoniadis, A., Cugliari, J., Fasiolo, M., Goude, Y., and Poggi, J.-M. (2024). *Statistical Learning Tools for Electricity Load Forecasting*. Springer.

Bosq, D. (2000). *Linear Processes in Function Spaces: Theory and Applications*, volume 149 of *Lecture Notes in Statistics*. Springer.

Chen, K.-T. (1958). Integration of paths—a faithful representation of paths by non-commutative formal power series. *Transactions of the American Mathematical Society*, 89:395–407.

Chevyrev, I. and Kormilitzin, A. (2025). A primer on the signature method in machine learning.

Cohen, S. N., Lui, S., Malpass, W., Mantoan, G., Nesheim, L., de Paula, A., Reeves, A., Scott, C., Small, E., and Yang, L. (2023). Nowcasting with signature methods. *arXiv preprint arXiv:2305.10256*.

De Vilmarest, J. and Goude, Y. (2022). State-space models for online post-covid electricity load forecasting competition. *IEEE Open Access Journal of Power and Energy*, 9:192–201.

Fermanian, A. (2021a). Embedding and learning with signatures. *Computational Statistics & Data Analysis*, 157:107148.

Fermanian, A. (2021b). *Learning time-dependent data with the signature transform*. Theses, Sorbonne Université.

Fermanian, A. (2022). Functional linear regression with truncated signatures. *Journal of Multivariate Analysis*, 192:105031.

Fermanian, A., Marion, P., Vert, J.-P., and Biau, G. (2021). Framing rnn as a kernel method: A neural ode approach.

Friz, P. K. and Victoir, N. B. (2010). *Multidimensional Stochastic Processes as Rough Paths: Theory and Applications*, volume 120 of *Cambridge Studies in Advanced Mathematics*. Cambridge University Press, Cambridge.

Goude, Y., Nedellec, R., and Kong, N. (2013). Local short and middle term electricity load forecasting

with semi-parametric additive models. *IEEE transactions on smart grid*, 5(1):440–446.

Guo, X., Wang, B., Zhang, R., and Zhao, C. (2024). On consistency of signature using lasso.

Jaber, E. A. and Sotnikov, D. (2025). Exponentially fading memory signature.

Lyons, T. (2014). Rough paths, signatures and the modelling of functions on streams.

Lyons, T., Caruana, M., and Lévy, T. (2007). *Differential Equations Driven by Rough Paths*, volume 1908 of *Lecture Notes in Mathematics*. Springer.

Moreno-Pino, F., Álvaro Arroyo, Waldon, H., Dong, X., and Álvaro Cartea (2025). Rough transformers: Lightweight and continuous time series modelling through signature patching.

Morrill, J., Fermanian, A., Kidger, P., and Lyons, T. (2021). A generalised signature method for multivariate time series feature extraction.

Nedellec, R., Cugliari, J., and Goude, Y. (2014). Gefcom2012: Electric load forecasting and backcasting with semi-parametric models. *International Journal of Forecasting*, 30(2):375–381.

Reizenstein, J. and Graham, B. (2020). Algorithm 1004: The iisignature library: Efficient calculation of iterated-integral signatures and log signatures. *ACM Transactions on Mathematical Software (TOMS)*.
# A 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.

## A.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] \to \mathbb{R}^d$, $u \mapsto x_u = (x_u^{(1)}, x_u^{(2)}, \dots, x_u^{(d)})$. We consider multi-dimensional paths, meaning $d \ge 2$. Moreover, we restrict ourselves to working with paths of finite length, which we define through the notion of the 1-variation.

**Definition.** Let $x : [s, t] \to \mathbb{R}^d$ be a continuous path. The 1-variation of $x$ is defined by

$$ \|x\|_{1-\text{var}} = \sup_{(t_0, \dots, t_k) \in P} \sum_{i=1}^{k} \|x_{t_i} - x_{t_{i-1}}\|, $$

where $P = \{(t_0, \dots, t_k) | k \ge 0, s = t_0 < \dots < t_k = t\}$ denotes the set of all finite partitions of $[s, t]$.

Intuitively, $\|x\|_{1-\text{var}}$ can be interpreted as the length of the path $x$, therefore we work with paths such that $\|x\|_{1-\text{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 Riemann-Stieltjes sense. We refer to Chapters 2 and 3 in Friz and Victoir (2010) for more details.

**Definition.** Let $x, y$ be two continuous paths from $[s, t]$ to $\mathbb{R}$, with $x$ being of finite length. Let $P^n = \{s = t_0^n < t_1^n < \dots < t_{N_n}^n = t\}$ for $n \ge 0$ be a sequence of partitions with vanishing mesh size $|P^n| = \max_{1 \le i \le N_n} |t_i - t_{i-1}|$, meaning $|P^n| \to 0$ as $n \to \infty$. The Riemann-Stieltjes integral of $y$ against $x$, is defined as

$$ \int_s^t y_u \, dx_u := \lim_{n \to \infty} \sum_{i=0}^{N_n-1} y_i^{(n)} (x_{t_{i+1}}^n - x_i^n), $$

where the limit does not depend on the choice of the sequence of partitions $(P^n)_{n \ge 0}$.

This definition can be generalized to the case when $x$ and $y$ are vector-valued; for $x, y : [s, t] \to \mathbb{R}^d$ we have

$$ \int_s^t y_u \, dx_u = \begin{pmatrix} \int_s^t y_u^{(1)} \, dx_u^{(1)} \\ \vdots \\ \int_s^t y_u^{(d)} \, dx_u^{(d)} \end{pmatrix}. $$

**Remark 6.** Take two continuous paths $x : [s, t] \to \mathbb{R}$ and $y : [s, t] \to \mathbb{R}$, where $x$ is continuously differentiable. Then the Riemann-Stieltjes integral of $y$ against $x$ comes down to

$$ \int_s^t y_u \, dx_u = \int_s^t y_u x'_u \, du, $$

where the last integral is the classical Riemann integral and $x'_u = dx_u/du$ denotes differentiation with respect to a single variable.

**Example 3.** Consider the constant path $y_u = 1$, $u \in [s, t]$. It follows that the integral of $y$ against any $x : [s, t] \to \mathbb{R}$ that is continuously differentiable is just the increment of $x$ on $[s, t]$:

$$ \int_s^t dx_u = \int_s^t x'_u \, du = x_t - x_s. $$

## A.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, \dots, e_d$ be the canonical basis of $\mathbb{R}^d$, then $(e_{i_1} \otimes \dots \otimes e_{i_k})_{1 \le i_1, \dots, i_k \le d}$ is an orthonormal basis of $(\mathbb{R}^d)^{\otimes k}$, meaning that any element $a \in (\mathbb{R}^d)^{\otimes k}$ can be written as $a = \sum_{1 \le i_1, \dots, i_k \le d} a^{(i_1, \dots, i_k)} e_{i_1} \otimes \dots \otimes e_{i_k}$, where $a^{(i_1, \dots, i_k)} \in \mathbb{R}$.
Furthermore, $(\mathbb{R}^d)^{\otimes k}$ is a Hilbert space of dimension $d^k$, with the following scalar product and norm:

$$ \langle a, b \rangle_{(\mathbb{R}^d)^{\otimes k}} = \sum_{1 \le i_1, \dots, i_k \le d} a^{(i_1, \dots, i_k)} b^{(i_1, \dots, i_k)}, \quad \|a\|_{(\mathbb{R}^d)^{\otimes k}} = \sqrt{\sum_{1 \le i_1, \dots, i_k \le d} (a^{(i_1, \dots, i_k)})^2} . $$

**Definition 4.** We define the extended tensor algebra $T((\mathbb{R}^d))$ as

$$ T((\mathbb{R}^d)) := \{ (a_0, a_1, \dots, a_k, \dots) \mid a_k \in (\mathbb{R}^d)^{\otimes k}, k \ge 0 \}. $$

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((\mathbb{R}^d))$ is a non-commutative algebra under the tensor product $\otimes$, with the neutral element $1 = (1, 0, 0, \dots)$. Furthermore, it holds that the subset $T_1((\mathbb{R}^n)) := \{\mathbf{a} \mid \mathbf{a} \in T((\mathbb{R}^n))$ with $a_0 = 1\}$ is a Lie group. Finally, the following subspace will be of special interest.

**Definition.** We denote with $\mathcal{T}((\mathbb{R}^d))$ the space of square-summable elements of $T((\mathbb{R}^d))$:

$$ \mathcal{T}((\mathbb{R}^d)) = \left\{ \mathbf{a} \in T((\mathbb{R}^d)) \mid \sum_{k=0}^{\infty} \|a_k\|_{((\mathbb{R}^d)^{\otimes k})}^2 < \infty \right\}. $$

We can endow this space with the following scalar product and associated norm: for any $\mathbf{a}, \mathbf{b} \in \mathcal{T}((\mathbb{R}^d))$,

$$ \langle \mathbf{a}, \mathbf{b} \rangle_{\mathcal{T}((\mathbb{R}^d))} = \sum_{k=0}^{\infty} \langle a_k, b_k \rangle_{((\mathbb{R}^d)^{\otimes k})}, \quad \|\mathbf{a}\|_{\mathcal{T}((\mathbb{R}^d))} = \sqrt{\sum_{k=0}^{\infty} \|a_k\|_{((\mathbb{R}^d)^{\otimes k})}^2}. $$

**Proposition (Proposition 1).** $(\mathcal{T}((\mathbb{R}^d)), \langle \cdot, \cdot \rangle_{\mathcal{T}((\mathbb{R}^d))})$ is a separable Hilbert space.

*Proof. (Completeness.)* The detailed proof that $(\mathcal{T}((\mathbb{R}^d)), \langle \cdot, \cdot \rangle_{\mathcal{T}((\mathbb{R}^d))})$ is a Hilbert space can be found in Fermanian et al. (2021), Appendix A.

(Separability.) Let $e_1, \dots, e_d$ be the canonical basis of $\mathbb{R}^d$, and $B_k = \{e_{i_1} \otimes \dots \otimes e_{i_k} \mid 1 \le i_1, \dots, i_k \le d\}$ denote the associated orthonormal basis of $(\mathbb{R}^d)^{\otimes k}$ that has $d^k$ elements. We now embed all of the elements of $B_k$, $k \ge 0$ in the space $\mathcal{T}((\mathbb{R}^d))$ by setting the other tensor sequence elements to zero, i.e. we define the sets

$$ \mathbf{B}_k = \left\{ \left( 0, 0, \dots, \overbrace{e_{i_1} \otimes \dots \otimes e_{i_k}}^{\text{k-th position}}, 0, \dots \right) \mid 1 \le i_1, \dots, i_k \le d \right\}, $$

with $\mathbf{B}_0 = \{(1, 0, 0, \dots)\}$. Let $\mathbf{B} = \bigcup_{k=0}^{\infty} \mathbf{B}_k$. The first thing to note is that $\mathbf{B}$ is a countable union of finite sets and is therefore countable. Secondly, $\mathbf{B}$ is a Schauder basis for $\mathcal{T}((\mathbb{R}^d))$, meaning that every element $\mathbf{a} \in \mathcal{T}((\mathbb{R}^d))$ can be uniquely represented as $\mathbf{a} = \sum_{n=1}^{\infty} \alpha_n \mathbf{v}_n$, $\mathbf{v}_n \in \mathbf{B}$, where the convergence of the infinite sum is the one of the topology induced by the norm (A.2). Indeed, we have that $\mathbf{a} = \sum_{k=0}^{\infty} (0, 0, \dots, a_k, 0, \dots)$ and $(0, 0, \dots, a_k, 0, \dots) = \sum_{1 \le i_1, \dots, i_k \le d} a^{(i_1, \dots, i_k)} (0, 0, \dots, e_{i_1} \otimes \dots \otimes e_{i_k}, 0, \dots)$, which gives us a unique representation of $\mathbf{a}$ in terms of elements of $\mathbf{B}$. Furthermore, from the fact that $B_k$ is an orthonormal basis for $(\mathbb{R}^d)^{\otimes k}$ and the definition of the scalar product (A.2), it is easy to see that $\mathbf{B}$ is also orthonormal, making it a countable orthonormal basis. Finally, a Hilbert space is separable if and only if it has a countable orthonormal basis, which leads to the conclusion that $(\mathcal{T}((\mathbb{R}^d)), \langle \cdot, \cdot \rangle_{\mathcal{T}((\mathbb{R}^d))})$ is a separable Hilbert space. $\square$

## A.3 Signatures

We can now fully define signatures, in the light of the theory presented earlier.
### A.3.1 Definitions

**Definition.** The signature of a finite-length path $x$ on $[s, t]$ is defined as an infinite tensor sequence:

$$S(x_{[s,t]}) = (1, S^1(x_{[s,t]}), \dots, S^k(x_{[s,t]}), \dots) \in T((\mathbb{R}^d)),$$

where the k-th element (called level) is given by

$$S^k(x_{[s,t]}) = \int \dots \int_{s<u_1<\dots<u_k<t} dx_{u_1} \otimes \dots \otimes dx_{u_k} \in (\mathbb{R}^d)^{\otimes k}.$$

Using the orthonormal basis for $(\mathbb{R}^d)^{\otimes k}$, we can write

$$S^k(x_{[s,t]}) = \sum_{(i_1, \dots, i_k) \subset \{1, \dots, d\}^k} S^{(i_1, \dots, i_k)}(x_{[s,t]}) e_{i_1} \otimes \dots \otimes e_{i_k},$$

where $S^{(i_1, \dots, i_k)}(x_{[s,t]})$ is referred to as the signature coefficient of $x$ along the multi-index $(i_1, \dots, i_k) \subset \{1, \dots, d\}^k, k \ge 1$ on $[s, t]$ and is given by

$$S^{(i_1, \dots, i_k)} = \int \dots \int_{s<u_1<\dots<u_k<t} dx_{u_1}^{(i_1)} \dots dx_{u_k}^{(i_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^{\le N}(x_{[s,t]}) = (1, S^1(x_{[s,t]}), \dots, S^N(x_{[s,t]})),$$

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 collection of all signature coefficients with multi-index of length $k \le N$, arranged in a vector

$$\left( 1, S^{(1)}(x_{[s,t]}), \dots, S^{(d)}(x_{[s,t]}), S^{(1,1)}(x_{[s,t]}), \dots, S^{(1,2)}(x_{[s,t]}), \dots, S^{(\overbrace{\underbrace{\dots, \dots, \dots}_{N \text{ times}}} d, \dots, d)}(x_{[s,t]}) \right),$$

that is of size $s_d(N) = \sum_{k=0}^{N} d^k = (d^{N+1} - 1)/(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.

**Example 4** (Example 1 revisited). Let $x : [s, t] \mapsto \mathbb{R}^d$, $u \to \frac{x_t - x_s}{t-s}(u-s) + x_s$ be a $d$-dimensional linear path.
For any $(i_1, \dots, i_k) \in \{1, \dots, d\}^k$ we have the following:

$$
\begin{align*}
S^{(i_1, \dots, i_k)} (x_{[s,t]}) &= \int \dots \int dx_{u_1}^{(i_1)} \dots dx_{u_k}^{(i_k)} = \int \dots \int \frac{x_t^{(i_1)} - x_s^{(i_1)}}{t-s} du_1 \dots \frac{x_t^{(i_k)} - x_s^{(i_k)}}{t-s} du_k \\
&= \prod_{j=1}^k \left( \frac{x_t^{(i_j)} - x_s^{(i_j)}}{t-s} \right) \cdot \int \dots \int du_1 \dots du_k \\
&= \prod_{j=1}^k \left( \frac{x_t^{(i_j)} - x_s^{(i_j)}}{t-s} \right) \cdot \int_s^t \int_s^{u_k} \dots \left( \int_s^{u_2} du_1 \right) \dots du_{k-1} du_k \\
&= \prod_{j=1}^k \left( \frac{x_t^{(i_j)} - x_s^{(i_j)}}{t-s} \right) \cdot \int_s^t \int_s^{u_k} \dots \left( \int_s^{u_3} (u_2-s) du_2 \right) \dots du_{k-1} du_k \\
&= \prod_{j=1}^k \left( \frac{x_t^{(i_j)} - x_s^{(i_j)}}{t-s} \right) \cdot \int_s^t \int_s^{u_k} \dots \left( \int_s^{u_4} \frac{(u_3-s)^2}{2} du_3 \right) \dots du_{k-1} du_k \\
&= \dots = \prod_{j=1}^k \left( \frac{x_t^{(i_j)} - x_s^{(i_j)}}{t-s} \right) \cdot \frac{(t-s)^k}{k!} = \frac{1}{k!} \prod_{j=1}^k (x_t^{(i_j)} - x_s^{(i_j)}).
\end{align*}
$$

More compactly written, $S^k(x_{[s,t]}) = \frac{1}{k!} (x_t - x_s)^{\otimes k}$.

### A.3.2 Properties

**Proposition 4 (Invariances).** Let $x : [s, t] \mapsto \mathbb{R}^d$ be a path of finite length. The following holds:

1. Let $\tilde{x}_u = x_{\psi(u)}$ be the smooth reparametrization of $x$, where $\psi : [s, t] \to [s, t]$ is a continuously differentiable non-decreasing surjection. Then, for any $[v, w] \subset [s, t]$ we have $S(\tilde{x}_{[v,w]}) = S(x_{[\psi(v),\psi(w)]})$.

2. Let $\bar{x}_u = x_u + a$ be a path obtained by translating $x$ by some $a \in \mathbb{R}^d$. Then $S(\bar{x}_{[s,t]}) = S(x_{[s,t]})$.

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\bar{x}_t = dx_t$.

**Proposition (Algebraic properties).** Let $x : [s, t] \mapsto \mathbb{R}^d$ and $y : [t, u] \mapsto \mathbb{R}^d$ denote two paths of finite length.

1. (Chen's identity) Let $x \ast y : [s, u] \mapsto \mathbb{R}^d$ be the concatenation of $x$ and $y$, meaning $(x \ast y)_v = x_v$ for $v \in [s, t]$ and $(x \ast y)_v = x_t + y_v - y_t$ for $v \in [t, u]$. Then $S((x \ast y)_{[s,u]}) = S(x_{[s,t]}) \otimes S(y_{[t,u]})$.

2. (Time reversal) We denote the time-reversal of $x$ as the path $\overleftarrow{x} : [s, t] \mapsto \mathbb{R}^d$ where $\overleftarrow{x}(u) = x_{s+t-u}$. Then $S(x_{[s,t]}) \otimes S(\overleftarrow{x}_{[s,t]}) = (1, 0, 0, \dots)$.

We can rephrase Chen's relation in terms of signature coefficients: for any multi-index $(i_1, \dots, i_k) \subset \{1, \dots, d\}^k$, it holds that

$$S^{(i_1, \dots, i_k)}((x \ast y)_{[s,u]}) = \sum_{\ell=0}^{k} S^{(i_1, \dots, i_\ell)}(x_{[s,t]}) \cdot S^{(i_{\ell+1}, \dots, i_k)}(y_{[s,t]}). \quad (1)$$

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.

**Remark 7.** Suppose that $x$ is a piecewise linear path such that for $s = u_1 < u_2 < \dots < u_k = t$, $x|_{[u_\ell, u_{\ell+1}]}$, $1 \le \ell \le k$, is linear. In order to obtain its signature, we first need to calculate the signature coefficients for each
linear segment, as in Example 4:

$$S^{(i_1, \dots, i_k)}(x_{[u_\ell, u_{\ell+1}]}) = \frac{1}{k!} \prod_{j=1}^{k} (x_{u_{\ell+1}}^{(i_j)} - x_{u_\ell}^{(i_j)}).$$

We proceed to calculate the signature on the whole time horizon $[s, t]$ by inductively concatenating the linear segments using the Chen's relation (1):

$$ \begin{aligned} S^{(i_1, \dots, i_k)}(x_{[s, u_{\ell+1}]}) &= \sum_{m=0}^{k} S^{(i_1, \dots, i_m)}(x_{[s, u_\ell]}) \cdot S^{(i_{m+1}, \dots, i_k)}(x_{[u_\ell, u_{\ell+1}]}) \\ &= \sum_{m=0}^{k} \left( S^{(i_1, \dots, i_m)}(x_{[s, u_\ell]}) \cdot \frac{1}{(k-m)!} \prod_{j=m+1}^{k} (x_{u_{\ell+1}}^{(i_j)} - x_{u_\ell}^{(i_j)}) \right). \end{aligned} \quad (2) $$

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` (Reizenstein and Graham, 2020).

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(x_{[s,t]}) \in T_1((\mathbb{R}^d))$. This result now states that the inverse of the signature of a path $x$ in the tensor group $T_1((\mathbb{R}^d))$ is exactly the signature of $x$ traversed backwards in time.

# B PROOFS OF THEORETICAL CONTRIBUTIONS

## B.1 Proof of Theorem 2

**Theorem** (Universal approximation on sliding windows). Let $(X_t)_t$ be a $d$-dimensional time series such that Assumption 1 holds. Let $w \in \mathbb{N}$ be a fixed window size and $f: C([0, w], \mathbb{R}^{d+1}) \to \mathbb{R}$ a continuous function given the uniform topology. It holds that for every $\varepsilon > 0$, there exists $N \in \mathbb{N}$, $\theta \in \mathbb{R}^{s_{d+1}(N)}$, such that

$$ \sup_t |f(X_{[t-w,t]} - (0, X_{t-w})) - \theta^\top S^{\le N} (X_{[t-w,t]})| \le \varepsilon. $$

*Proof.* We first introduce $\bar{X}_{[t-w,t]} := X_{[t-w,t]} - (0, X_{t-w})$ as a random path obtained by linear interpolation of $(0,0), (1, X_{t-w+1} - X_{t-w}), \dots, (w, X_t - X_{t-1})$. In other words, we consider paths on windows translated by the starting point. We now define $U = \{\bar{X}_{[t-w,t]}(\omega) \mid \omega \in \Omega; t \in \mathbb{Z}\}$ as the set of all realisations of random paths on all sliding windows. We note that $U$ is a set of piecewise linear functions whose first coordinate (corresponding to time) is strictly monotone and such that $|(\bar{X}_{[t-w,t]}^{(i)}(\omega))| \le \max(2M, w) =: B$ for all $1 \le i \le d$. It follows that they are therefore Lipschitz-continuous with the same Lipschitz constant $L := \max(4M, 1)$.

Let us define a set of paths $K \subset C([0, w], \mathbb{R}^{d+1})$ such that:

1. $\forall x \in K$ has at least one coordinate that is linear, i.e. $\exists j \in \{1, 2, \dots, d\}, x_u^j = a_x u$ and such that there exists $\varepsilon > 0$ such that $a_x \ge \varepsilon \forall x \in K$,

2. $\forall x \in K$ has the same starting point at 0, i.e. $x_0 = 0 \in \mathbb{R}^{d+1}$,

3. the paths are uniformly bounded by $B$, i.e. $|x_u^{(i)}| \le B \forall i \in \{1, 2, \dots, d\}, \forall x \in K$,

4. the paths are Lipschitz with the same constant $L$, i.e. $\|x_v - x_u\|_{\infty} \le L|v - u|, \forall x \in K$.

Functions in $K$ have the same Lipschitz constant, therefore they are also uniformly equicontinuous. Furthermore, it can be shown that the set $K$ is also closed in $C([0, w], \mathbb{R}^{d+1})$. We can now use the Arzelà-Ascoli theorem to conclude that $K$ is compact in $C([0, w], \mathbb{R}^{d+1})$ in the uniform topology. This allows us to apply the Theorem 1 to $K$.
Leveraging Theorem 1 and the fact that $U \subset K$, we obtain:

$$ \sup_U \left| f(\bar{X}_{[t-w,t]}(\omega)) - \theta^\top S^{\le N} (\bar{X}_{[t-w,t]}) \right| \le \sup_K \left| f(x) - \theta^\top S^{\le N} (x_{[0,w]}) \right| \le \varepsilon. $$

We conclude to the desired result using the identity $S^{\le N}(\bar{X}_{[t-w,t]}) = S^{\le N}(X_{[t-w,t]})$ following from the invariance by translation of the signatures; see Proposition 4, 2. □

## B.2 Proof of Theorem 3

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

**Theorem.** Let $(X_t)_{t \ge 1}$ be a discrete $d$-dimensional time series such that the Assumption 2 is satisfied. The following statements hold:

i) For any multi-index $(i_1, \dots, i_k)$ of any length $k \in \mathbb{N}$, the time series $(S^{(i_1, \dots, i_k)}(X_{[t-w,t]}))_t$ is strictly stationary.

ii) For any truncation order $N \in \mathbb{N}$, the time series $(S^{\le N}(X_{[t-w,t]}))_t$ is strictly stationary.

iii) The time series $(S(X_{[t-w,t]}))_t$ is strictly stationary.

*Proof.* i) For any multi-index $I = (i_1, \dots, i_k)$ of length $k \in \mathbb{N}$, we denote with $f^I$ the function from $\mathbb{R}^{w \times d}$ to $\mathbb{R}$ given as:

$$ (X_{t-w+1} - X_{t-w}, \dots, X_t - X_{t-1}) \xrightarrow{f^I} S^I(X_{[t-w,t]}). $$

Examining the iterative formula (2), we see that the $f^I$ (the signature component $S^I(X_{[t-w,t]})$) is the same deterministic, measurable function (as a composition of measurable functions: multiplication, sums and powers) of the increments $(X_{t-w+1} - X_{t-w}, X_{t-w+2} - X_{t-w+1}, \dots, X_t - X_{t-1})$ of the original time series $(X_t)_t$. Furthermore, it does not depend on $t$. We also note that if the $d$-dimensional time series $(X_t - X_{t-1})_t$ is strictly stationary, then it follows straight from the definition of strict stationarity that the $(d \times w)$-dimensional time series $(X_{t-w+1} - X_{t-w}, \dots, X_t - X_{t-1})_t$ is also strictly stationary. If we now remember that for any deterministic, measurable function $g$ and any random vectors $\mathbf{X}, \mathbf{Y}$ such that $\mathbf{X} \stackrel{D}{=} \mathbf{Y}$, we have $g(\mathbf{X}) \stackrel{D}{=} g(\mathbf{Y})$, and combine it with the previous remark for $f^I$, we have the strict stationarity for $(S^I(X_{[t-w,t]}))_t$.

ii) Let us now fix an $N$ and consider the $s_d(N)$-dimensional time series of the signature truncated at level $N$, $(S^{\le N}(X_{[t-w,t]}))_t$. Let $I_1, I_2, \dots, I_{s_d(N)}$ be the multi-indexes of the components of the truncated signature, in alphabetical order. We now have:

$$ \begin{aligned} f^{\le N}(X_{t-w+1} - X_{t-w}, \dots, X_t - X_{t-1}) &:= (f^{I_1}, \dots, f^{I_{s_d(N)}}) (X_{t-w+1} - X_{t-w}, \dots, X_t - X_{t-1}) \\ &= (S^{I_1}(X_{[t-w,t]}), \dots, S^{I_{s_d(N)}}(X_{[t-w,t]})) \\ &= S^{\le N}(X_{[t-w,t]}), \end{aligned} $$

where $f = (f^{I_1}, \dots, f^{I_{s_d(N)}}) : \mathbb{R}^{d \times w} \to \mathbb{R}^{s_d(N)}$ is a deterministic, measurable function of the increments that does not depend on $t$, as noted above. Following the same arguments as before, we can now conclude that $(S^{\le N}(X_{[t-w,t]}))_t$ is strictly stationary.

iii) It is easy to see that the previous proof holds not only for the time series of signatures truncated at level $N$, but also for the time series of any collection of signature components denoted by multi-indexes $(I_1, \dots, I_j)$. Additionally, as discussed in Remark 3, the finite-dimensional distributions of the signature components uniquely determine the law of the entire infinite-dimensional time series of signatures. From these two arguments, it follows that $(S(X_{[t-w,t]}))_t$ is strictly stationary. □
## C 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.

**Algorithm 2:** Fitting a ridge regression on signatures

**Input:** data $\{(x_0, y_0), \dots, (x_n, y_n)\}$, truncation order $N$, delay $D$, window size $w$

1. Split the data into the train and validation and test set by fixing the horizons $w \le t_{train} < t_{valid} \le n$;

2. **for** $\lambda$ **in set of possible regularization parameters** **do**

3. Compute the signatures $S(x_{[t-w,t]})$ for $w \le t \le t_{train}$ using Algorithm 1;

4. Fit a ridge regression model with parameter $\lambda$ on pairs $\{(S(x_{[t-w,t]}), \Delta_D y_t) | w \le t \le t_{train}\}$;

5. Compute forecasts $\hat{Y}_t = y_{t-D} + \hat{\theta}^\top S(x_{[t-w,t]})$ and error metric on the validation set;

6. Identify the best regularisation parameter $\hat{\lambda}$ (smallest error metric on the validation set);

7. Refit the ridge regression with parameter $\hat{\lambda}$ on the combined train and validation set;

**Output:** Ridge coefficients estimator $\hat{\theta}$

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

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

## D 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 \ge 1$.

**Synthetic data.** Electricity consumption depends on many factors, including calendar variables such as week-days 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

$$ \begin{cases} \bar{T}_1^\alpha = T_1 \\ \bar{T}_t^\alpha = (1-\alpha)\bar{T}_{t-1}^\alpha + \alpha T_t^\alpha, & \text{for any } t \ge 2, \end{cases} $$

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 $\bar{T}_t$ and $\bar{T}_t^2$ as covariates:

$$ \hat{Y}_t^\alpha = \hat{\theta}_1 \bar{T}_t^\alpha + \hat{\theta}_2 (\bar{T}_t^\alpha)^2, \quad (3) $$

thereby imposing a quadratic dependence of simulated demand on smoothed temperature. The final synthetic demand time series $(\tilde{Y}_t^\alpha)_t$ is obtained by adding normally distributed noise to the fitted values $\hat{Y}_t^\alpha$:

$$ \tilde{Y}_t^\alpha = \hat{\theta}_1 \bar{T}_t^\alpha + \hat{\theta}_2 (\bar{T}_t^\alpha)^2 + \varepsilon_t, \quad \varepsilon_t \stackrel{\text{i.i.d}}{\sim} \mathcal{N}(0, \sigma^2). \quad (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 5, 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.

<figure><img src="image_8.png" /><figcaption>Figure 5: Comparison of observed and simulated electrical demand. Left: full four-year dataset showing that the synthetic series captures yearly patterns driven by temperature. Right: zoom-in on a weekly segment (Feb 11–17, 2013) illustrating that weekly fluctuations driven by calendar variables are not retained in synthetic demand.</figcaption></figure>

**Assumptions in practice.** As discussed earlier, Assumption 1 is satisfied since temperature values in the dataset remain within natural bounds (between -10°C and 35°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.