# Learning the Dynamics of Sparsely Observed Interacting Systems

Linus Bleistein  Adeline Fermanian  Anne-Sophie Jannot  Agathe Guilloux

###### Abstract

We address the problem of learning the dynamics of an unknown non-parametric
system linking a target and a feature time series. The feature time series is
measured on a sparse and irregular grid, while we have access to only a few
points of the target time series. Once learned, we can use these dynamics to
predict values of the target from the previous values of the feature time
series. We frame this task as learning the solution map of a controlled
differential equation (CDE). By leveraging the rich theory of signatures, we
are able to cast this non-linear problem as a high-dimensional linear
regression. We provide an oracle bound on the prediction error which exhibits
explicit dependencies on the individual-specific sampling schemes. Our
theoretical results are illustrated by simulations which show that our method
outperforms existing algorithms for recovering the full time series while
being computationally cheap. We conclude by demonstrating its potential on
real-world epidemiological data.

Signature Methods, Time series, Dynamical Systems, Applications to Healthcare.

##  1 Introduction

Time series are ubiquitous in many areas such as finance, economics, robotics,
agriculture, and healthcare. One is typically interested in modelling the
evolution of a target quantity through time, which is known to be affected by
a set of time-evolving features. For example, pollution levels in a city are
driven by quantities such as temperature, pressure, traffic, or economic
activity measured through time. Mathematically, one wishes to model the
evolution of a quantity $y_{t}\in{\mathbb{R}}^{p}$, $p\geq 1$, as a function
of some time evolving features $x_{t}\in{\mathbb{R}}^{d}$, $d\geq 1$, for
$t\in[0,1]$. In other words, the goal is to learn the dynamics that link the
target to the features.

Such an interaction is typically modelled via differential equations, which
are a common choice of model in natural sciences (Zwillinger, 1989). In this
article, we assume that there exists a function
$\mathbf{G}:{\mathbb{R}}^{p}\times{\mathbb{R}}^{d}\to{\mathbb{R}}^{p}$ such
that

$y_{t}=y_{0}+\int_{0}^{t}\mathbf{G}(y_{s},x_{s})ds$ (1)

or equivalently

$dy_{t}=\mathbf{G}(y_{t},x_{t})dt,\quad y_{0}\in\mathbb{R}^{p}.$

The value $y_{t}$ depends on the trajectory of the features time series
$x_{s}$ up to time $t$. Learning the dynamics of the system can be framed as
learning the solution map of (1), i.e., a function $\Psi$ which, given a time
$t$, an initial point $y_{0}\in\mathbb{R}^{p}$, and the history of the path up
to time $t$, denoted by $x_{[0,t]}=(x_{s})_{s\in[0,t]}$, outputs the value of
$y$ at time $t$.

If we know $\Psi$, we gain access to the values of $y$ at any point in time
provided we know the values of $x$ up to this point ; this encompasses many
tasks such as forecasting or interpolating between points of $y$. We
specifically have in mind applications where we have an easy access to $x$ but
a limited one to $y$.

This problem is extremely common in healthcare. For example, in obstetrics,
the lactic acidosis (LA) of the fetus, which is a proxy for fetal distress, is
a quantity of high medical interest for predicting complications in the first
hours after birth. This biomarker cannot be measured during pregnancy but only
at birth because the measurement is highly invasive. Some vitals such as heart
rate and fetal movement are however easy to measure during pregnancy. In this
case, $x$ is the non-invasive measurements made during pregnancy, while $y$ is
the invasive measurement of LA at birth. Predicting the value of $y$ at any
time $t$ (both before and at birth) would allow for early diagnosis.
Similarly, after surgery, patients are often monitored to detect hemorrhage.
While some vitals such as heart rate of saturation are monitored in continuous
time, haemoglobin—which is highly predictive of hemorrhage—is only measured by
blood samples taken a few times a day, which can significantly delay
hemorrhage diagnostic.

#### Irregular data.

In practice, the functions $x$ and $y$ are measured on discrete grids and take
the form of time series. These often present a lot of heterogeneity, both
within and across individuals.

  * $(i)$

For every individual, the time between any two measurements can vary, and thus
individuals may not be recorded on the same grid.

  * $(ii)$

The number of total sampling points might vary between individuals.

  * $(iii)$

Each measurement in time might be corrupted by measurement noise.

Mathematically, we consider $n$ pairs of functions
$\\{(x^{1},y^{1}),\dots,(x^{n},y^{n})\\}$. Each $x^{i}$ deterministically
produces a specific $y^{i}$ through the Ordinary Differential Equation (ODE)
(1). We call $x^{i}$ the feature path and $y^{i}$ the target path. Both
$x^{i}$ and $y^{i}$ are only observed at a finite set of times specific to
every individual. We denote by

$D^{i}=\left(t^{i}_{1},\dots,t^{i}_{k_{i}}\right),\quad i=1,\dots,n,$

the sampling grid of $x^{i}$ and by $\bar{D}^{i}$ the sampling grid of
$y^{i}$. We stress that both the number of sampling times $k_{i}$ and the
sampling times $t^{i}_{1},\dots,t^{i}_{k_{i}}$ themselves are individual
specific, as described in $(i)$ and $(ii)$. Moreover, the observations are
corrupted by additive noise, such that we observe

$X^{i}_{t}=x^{i}_{t}+\xi^{i}_{t}$

for all $t\in D^{i}$, and similarly $Y_{t}^{i}=y^{i}_{t}+\varepsilon^{i}_{t}$
for every $t\in\bar{D}^{i}$, where the $\xi^{i}_{t}$ and $\varepsilon^{i}_{t}$
are sub-gaussian i.i.d. random vectors. Each input may therefore be written as
a matrix $\mathbf{X}^{i}=(X^{i}_{t})_{t\in D^{i}}\in\mathbb{R}^{k_{i}\times
d}$ which we call the feature time series. Similarly, the quantity of interest
is a matrix
$\mathbf{Y}^{i}=(Y^{i}_{t})_{t\in\bar{D}^{i}}\in\mathbb{R}^{m_{i}\times p}$
(where $m_{i}$ is the length of $\bar{D}^{i}$) and is called the target time
series. The grid $\bar{D}^{i}$ is assumed to be a subset of $D^{i}$: in our
setup $y^{i}$ is hard to sample and therefore measured at only a few points
(and sometimes only one) while $x^{i}$ is easy to access and measured at high
frequency. Our goal is to approximate the dynamics linking $x$ and $y$ from
the irregular, heterogeneous, and fuzzy data $\mathbf{X}^{i}$ and
$\mathbf{Y}^{i}$.

Such heterogeneity is difficult to handle by classical machine learning
algorithms such as Long short-term memory networks (LSTM, Hochreiter &
Schmidhuber, 1997) which assume that the data is regularly sampled. Some more
recent approaches (Rubanova et al., 2019; De Brouwer et al., 2019; Kidger et
al., 2020) have adapted these models by introducing continuously evolving
hidden states to account for the irregular spacing between observation times.

We build upon the approach of Neural Controlled Differential Equations (Neural
CDE, Kidger et al., 2020; Morrill et al., 2021b), which have proven to be very
successful for time series classification and online prediction tasks (Morrill
et al., 2021a). The key idea of Neural CDE is that under some fairly mild
assumptions, any general ordinary differential equation of the form (1) can be
rewritten as

$y_{t}=y_{0}+\int_{0}^{t}\mathbf{F}(y_{s})dx_{s},$ (2)

where $\mathbf{F}:{\mathbb{R}}^{p}\to{\mathbb{R}}^{p\times d}$ is a matrix-
valued vector field, such that the right-hand-side of (2) is a matrix-vector
product (see, e.g., Fermanian et al., 2021, Proposition 2, for a proof). The
function $x$ is often called the driver of the CDE. In a Neural CDE setting,
the driver $x$ is a continuous interpolation of the feature times series, $y$
corresponds to a continuously-evolving state, and $\mathbf{F}$ is chosen to be
a neural network. This network is then trained such that the values of
$(y_{t})$ can be used as features for classification or regression tasks.
While Neural CDE have been shown to outperform other architectures with
limited memory usage, their training time is considerable and no statistical
guarantees exist.

#### Model.

We model the interactions between the target and the feature paths through a
CDE of the form (2). This modelling choice encapsulates a broad variety of
settings, since the vector field $\mathbf{F}$ can be any (regular enough)
function. A priori, the solution map $\Psi$ of this CDE is a complex function
of time and the history of $x$ up to $t$; however, by linearizing the model,
we are able to approximate $\Psi$ by a simple scalar product between a
deterministic transformation of the history of $x$, called the signature of
$x$ at order $N\geq 1$ and denoted by $S_{N}(x_{[0,t]})$, and a time
independent matrix $\theta_{N}^{\ast}$. Informally, we have

$\Psi\big{(}x_{[0,t]},t\big{)}\approx
S_{N}\big{(}x_{[0,t]}\big{)}^{\top}\theta_{N}^{\ast}.$

Two striking features of this linearized model are (i) that $\theta^{*}_{N}$
can be learned on any time horizon $[0,t]$ since it is independant of time,
and (ii) that once it has been learned, the model can be called at any time
$t$.

#### Contributions.

Our contributions are threefold. First, we frame the task of learning the
interactions between two time series as learning the flow of a CDE, which can
be linearized in the signature space. While the connection between CDEs and
signatures is well-known, this is the first time CDEs are used as a
statistical model. We then leverage this linearization to derive statistical
guarantees on the prediction error with an explicit dependence on both
sampling irregularities and the noise affecting measurements. To our
knowledge, this is the first bound of this type for signature-based models,
allowing for better understanding of the dependencies between prediction
performance and sampling roughness. Finally, the resulting algorithm, called
SigLasso, is shown to be computationally cheap and competitive compared to
existing baselines on a wide range of simulated data and a real-world example
of hospitalization growth rate prediction during the Covid pandemic.

#### Related works.

Signatures originated as a prominent tool in stochastic analysis (Chen, 1958;
Lyons et al., 2007; Friz & Victoir, 2010) and have proven to be a powerful
feature extraction method in machine learning in various domains such as
healthcare (Morrill et al., 2020b; Wang et al., 2020), human action
recognition (Yang et al., 2022), or financial modelling (Lyons et al., 2014;
Buehler et al., 2020). Their appealing properties include a capacity to handle
irregular data, to capture dependence between coordinates, and their links
with the theory of CDE. We refer to Lyons & McLeod (2022) for a recent survey
on their use cases. However, the statistical properties of signatures based
algorithms have received little attention so far, with a few notable
exceptions (Papavasiliou & Ladroue, 2011; Lemercier et al., 2021; Fermanian,
2022).

On the other hand, the interplay between dynamical systems and machine
learning has received considerable attention in the recent years. A first line
of work has focused on approximating the solution of ODE and Partial
Differential Equations (PDE) with neural networks (Lagaris et al., 1998; Han
et al., 2018; Zubov et al., 2021) and directly learning dynamical systems
(Long et al., 2018; Fattahi et al., 2019). Recent approaches have been
interested in combining deep learning algorithms with physical knowledge
(Greydanus et al., 2019; Brunton et al., 2020; Willard et al., 2020). Finally,
dynamical systems, seen as continuous versions of neural network
architectures, have also been a great source of inspiration for analysing and
designing machine learning algorithms in the recent years (Chen et al., 2018;
Fermanian et al., 2021; Marion et al., 2022). We refer to Kidger (2022) for an
extensive review.

Our problem also bears close resemblance to frameworks encountered in
sequence-to-sequence learning (Sutskever et al., 2014; Gehring et al., 2017)
and functional regression (Ramsay & Dalzell, 1991; Marx & Eilers, 1999).

#### Overview.

Section 2 introduces the CDE model for interacting systems, the mathematical
context and the learning procedure. Our main theoretical result is presented
in Section 3. We conclude by an empirical study on synthetic and real-world
data in Section 4. All proofs are postponed to the appendix and the code to
reproduce the experiments is available on our Github.

##  2 Model and assumptions

###  2.1 A CDE-based model on the dynamics

We start by describing our assumptions on the feature and target paths, which
are linked by Equation (2). To correctly define the integral of Equation (2),
we need to impose some conditions on the $x^{i}$ and on $\mathbf{F}$. Note
that we consider that the $x^{i}$ are defined on $[0,1]$ but our results
extend easily to any compact time interval $[a,b]$.

######  Assumption 1.

All paths $(x^{i})_{1\leq i\leq n}$ are continuous and there exists $0<L<1$
such that, for all $i=1,\dots,n$,

$\big{\|}x^{i}\big{\|}_{\textnormal{1-var},[0,1]}=\sup_{D}\sum\limits_{k}\big{\|}x^{i}_{t_{k+1}}-x^{i}_{t_{k}}\big{\|}\leq
L,$

where $\left\lVert\cdot\right\rVert$ is the Euclidean norm and the supremum is
taken over all finite dissections $D=\left\\{0=t_{1}<\dots<t_{k}=1\right\\}$.

We write $C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$ for the set of
continuous paths of total variation bounded by $L$. Outside the statistical
context, when referring to general paths, we will drop the superscripts $i$
and simply write $x$ and $y$ to alleviate notations.

We assume that the target path $y$ is the solution of the ODE (1). This
modelization choice means that the evolution of $y$ is governed by a dynamical
system whose dynamics itself are allowed to vary with the current value of the
feature path. Observe that this model can be seen as a generalized form of a
non-autonomous system (Lyons et al., 2007), which we recover by taking
$x_{t}=t$. Since Equation (1) can be rewritten as a CDE, the starting point of
our work is to assume that the true dynamics of the data follow such a CDE, as
stated in the following assumption.

######  Assumption 2.

There exists a smooth vector field
$\mathbf{F}:{\mathbb{R}}^{p}\to{\mathbb{R}}^{p\times d}$ such that, for all
$i=1,\dots,n$, $y^{i}$ is the solution of the CDE (2) driven by $x^{i}$ with
initial condition $y^{i}_{0}=y_{0}\in\mathbb{R}^{p}$ homogeneous amongst
individuals.

By “smooth” we mean that each coordinate of $\mathbf{F}$ is infinitely
differentiable, that is, is $\mathcal{C}^{\infty}$. The vector field
$\mathbf{F}$ and the initial condition $y_{0}$ are common to all individuals,
which can be seen as homogeneity assumptions on our sample. On the other hand,
since every individual $i$ has her own feature path $x^{i}$, the target paths
$y^{i}$ are individual specific. In other words, there exists a solution map
$\Psi$ that depends only on $y_{0}$ and $\mathbf{F}$ and is such that, for any
$t\in[0,1]$, $\Psi(x^{i}_{[0,t]},t)=y^{i}_{t}$.

The vector field $\mathbf{F}$ encapsulates the common physical dynamics
governing the evolution of $y^{i}$, which are affected by the changes in
$x^{i}$. Note that there is no parametric model on $\mathbf{F}$ (although some
strong smoothness requirements will be needed) contrarily to functional or
traditional time series models (Ramsay & Silverman, 2005; Morris, 2015).

###  2.2 Linearizing the CDE with signatures

Before defining the Taylor expansion of the CDE (2), which will allow us to
linearize the solution map $\Psi$, we need to introduce the notion of
signature, which have emerged as a powerful tool to model time series (Levin
et al., 2013; Kidger et al., 2019).

From now on, for any feature path $x$ with values in $\mathbb{R}^{d}$, we
denote by $x^{(j)}$ its $j$th coordinate, for $j=1,\ldots,d$.

######  Definition 2.1.

Let $x\in C^{\textnormal{{1-var}}}_{L}([0,1],\mathbb{R}^{d}).$ Take a word of
length $k$ from the alphabet $\\{1,\dots,d\\}$, that is, an element
$I=(i_{1},\dots,i_{k})\in\\{1,\dots,d\\}^{k}$. For all $t\in[0,1]$, the
signature coefficient associated to this word is the scalar

$\displaystyle S^{I}\big{(}x_{[0,t]}\big{)}$
$\displaystyle=\idotsint\limits_{0<u_{1}<\dots<u_{k}<t}dx^{(i_{1})}_{u_{1}}\cdots
dx^{(i_{k})}_{u_{k}}.$

Signatures encode geometric properties of paths and have numerous appealing
properties as a feature set for time series. In the following, we draw on
their connection with Taylor expansions of CDEs.

We introduce a series of notation on signature coefficients, grouping them by
the length $k$ of the words. For any $k\in\mathbb{N}$ and $t\in[0,1]$, the
signature of order $k$ of $x_{[0,t]}$ is

$\displaystyle\mathbb{X}_{k,[0,t]}=\big{(}S^{I}(x_{[0,t]})\big{)}_{I\in\\{1,\dots,d\\}^{k}}\in\mathbb{R}^{d^{k}},$

where the words $(i_{1},\dots,i_{k})$ are in lexicographic order. We then
denote the full signature by

$\displaystyle
S(x_{[0,t]})=\big{(}1,\mathbb{X}^{1}_{[0,t]},\dots,\mathbb{X}^{k}_{[0,t]},\dots\big{)}$

and, for any $N\geq 1$, the signature truncated at order $N$ by

$\displaystyle
S_{N}(x_{[0,t]})=\big{(}1,\mathbb{X}^{1}_{[0,t]},\dots,\mathbb{X}^{N}_{[0,t]}\big{)}^{\top}.$

The size of the signature truncated at order $N$ grows exponentially with $N$,
since, for $d\geq 2$, it is equal to

$\displaystyle s_{d}(N)=1+d+d^{2}+\dots+d^{N}=\frac{d^{N+1}-1}{d-1}.$

When computing signatures, it is common practice to add time as a coordinate
to the path (Fermanian, 2021), that is, consider the path $(t,x_{t})^{\top}$.
From now on, we assume that the first dimension of $x$ always corresponds to
time, so that $d\geq 2$.

We are now ready to define the Taylor expansion of a CDE.

######  Definition 2.2.

Let $N\geq 1$. The Taylor expansion of order $N$ of the solution $y$ of
Equation (2) is defined by

$\displaystyle\overline{y}_{N,t}=y_{0}+\sum\limits_{k=1}^{N}\sum\limits_{I\in\\{1,\dots,d\\}^{k}}S^{I}\left(x_{[0,t]}\right)\times\Phi^{I}_{\mathbf{F}}\left(y_{0}\right),$
(3)

where $\Phi^{I}_{\mathbf{F}}(\cdot)\in{\mathbb{R}}^{p}$ is the differential
product of the vector field $\mathbf{F}$ along $I$, whose definition is
postponed to Appendix A.3.

The differential products $\Phi^{I}_{\mathbf{F}}(y_{0})$ are essentially a
combination of multiplication and summation of derivatives of the different
components of $\mathbf{F}$, evaluated at $y_{0}$.

The Taylor expansion crucially allows to write the solution map as a product
between a time-varying term (the signature of the feature path) and a
constant-over-time term (the differential product). This is similar to a
regular Taylor expansion; we discuss this analogy in Appendix A.4. Note that
Equation (3) may also be written as a matrix-vector product

$\displaystyle\overline{y}_{N,t}{}^{\top}$
$\displaystyle=S_{N}(x_{[0,t]})^{\top}\theta^{\ast}_{N}\in\mathbb{R}^{p},$ (4)

where $\theta^{\ast}_{N}\in{\mathbb{R}}^{s_{d}(N)\times p}$ is a matrix
collecting all differential products $\Phi^{I}_{\mathbf{F}}(y_{0})$ up to
order $N$ and the offset $y_{0}$. Since $\theta^{\ast}_{N}$ depends neither on
$x$ nor on $t$, the Taylor expansion of $y$ at order $N$ is simply a linear
function of the truncated signature.

For this expansion to become exact in the sens of pointwise convergence, we
need quite strong regularity assumption on $\mathbf{F}$. This is the price to
pay for the non-parametric nature of our model.

######  Assumption 3.

The vector field $\mathbf{F}$ has fast decaying derivatives. In other words,
defining

$\displaystyle\Lambda_{k}(\mathbf{F})=\sup\limits_{I\in\\{1,\dots,d\\}^{k}}\left\lVert\Phi^{I}_{\mathbf{F}}\left(y_{0}\right)\right\rVert,$

we assume that $\sum_{k=0}^{\infty}d^{k}\Lambda_{k}(\mathbf{F})/k!<\infty$.

Let $y$ be the solution of the CDE (2), and let
$(\overline{y}_{N,t})_{t\in[0,1]}$ be its $N$-th Taylor expansion. Under
Assumptions 1, 2, and 3, we then have, for any $t\in[0,1]$,

$\displaystyle\left\lVert
y_{t}-\overline{y}_{N,t}\right\rVert\underset{N\to+\infty}{\longrightarrow}0.$
(5)

We refer the reader to Friz & Victoir (2008) and Fermanian et al. (2021,
Proposition 4) for more details on this result.

While we impose relatively strong regularity assumptions on the vector field
$\mathbf{F}$, the assumptions on $x$ are mild (see Assumption 1). Therefore,
our assumptions, while enforcing a fairly high amount of regularity on the
dynamical structure, still accommodate most of the real world data. Also note
that Fermanian et al. (2021) give conditions under which Assumption 3 is
satisfied when $\mathbf{F}$ is a layer of a neural network with smooth
activation functions.

###  2.3 The learning problem

Figure 1: The workflow of our model. Starting from the left, the first panel
describes our modelling hypothesis: the target and feature time series are
linked through an unknown CDE. The second panel shows the observed data. The
third panel shows how every observation of the target is mapped to the
signature of the corresponding path, how to construct the dataset
$\mathbf{S}_{N}^{\mathcal{D}},\mathbf{Y}$ from this data and finally how to
learn the SigLasso estimator.

We go back to the statistical learning problem. Recall that in practice, we do
not observe the continuous paths $\\{(x^{1},y^{1}),\dots,(x^{n},y^{n})\\}$ but
their discretized and fuzzy counterparts
$\\{(\mathbf{X}^{1},\mathbf{Y}^{1}),\dots,(\mathbf{X}^{n},\mathbf{Y}^{n})\\}$
measured on a set of individual grids $\mathcal{D}=\\{D^{1},\dots,D^{n}\\}$
and $\bar{\mathcal{D}}=\\{\bar{D}^{1},\dots,\bar{D}^{n}\\}$. We allow for high
data heterogeneity: two individuals can be sampled at very different
frequencies and at different time-points, and therefore have different
observation grids.

The meshsize of a sampling grid $D$, denoted by $\lvert D\rvert$ is defined as
the largest gap between two successive sampling times, that is,

$\lvert D\rvert=\max\limits_{t_{i}\in D}\,\lvert t_{i+1}-t_{i}\rvert.$

Its cardinality, denoted by $\\#D$, is the number of sampling points in $D$.
We make the following assumption on the sampling procedure.

######  Assumption 4.

For all $i=1,\dots,n$,

$0\in D^{i},\quad\\#D^{i}\geq 2,\quad\text{and}\quad\bar{D}^{i}\subset D^{i}.$

This assumption essentially means that the measurements on all individuals
start at $0$ and that the sampling grid of $\mathbf{Y}^{i}$ is coarser than
the sampling grid of $\mathbf{X}^{i}$. We also let
$\lvert\mathcal{D}\rvert=\max_{i}\lvert D^{i}\rvert$ be the biggest gap
between two successive sampling times within the whole set of individual grids
and $\\#\mathcal{D}=\sum\limits_{i=1}^{n}\\#D^{i}$ the total number of
sampling points of the feature time series. We recall that the random vectors
$\xi_{t}^{i}$, $t\in D^{i}$, are the noises affecting the measurements of
$x^{i}$. The random vectors $\varepsilon^{i}_{t}$ affect the measurements of
$y^{i}$. We end with assumptions on the law of these measurement noises.

######  Assumption 5 (Noise on the feature time series).

The noises $(\xi^{i}_{t})_{i\in\\{1,\dots,n\\},t\in D^{i}}$ are i.i.d.
$v_{\xi}$-subgaussian random vectors.

######  Assumption 6 (Noise on the target time series).

The noises $(\varepsilon^{i}_{t})_{i\in\\{1,\dots,n\\},t\in\bar{D}^{i}}$ are
i.i.d. $v_{\varepsilon}$-subgaussian random vectors and independent from
$(\xi^{i}_{t})_{i\in\\{1,\dots,n\\},t\in D^{i}}$.

The goal of the learning procedure is to learn the solution map $\Psi$, in
order to infer the value of $y$ at any time $t$, given observations of $x$ up
to time $t$. We have seen previously that this problem boils down to
estimating the matrix $\theta_{N}^{\ast}$ via Equation (4). The truncation
order $N\geq 1$ is an hyperparameter and will be selected by cross-validation.

#### Single target measurement.

For the sake of simplicity, we first present the case where $Y^{i}$ is
measured only at the end of the observation period, as in the example of LA
measurements in obstetrics. In this case, for all $i=1,\dots,n$, we have
$m_{i}=1$ such that $M=n$,

$D^{i}=(t^{i}_{1},\dots,t^{i}_{k_{i}}),\quad\text{and}\quad\bar{D}^{i}=(t^{i}_{k_{i}}).$

Then $\mathbf{Y}^{i}\in{\mathbb{R}}^{p}$, and we denote by

$\mathbf{Y}=\begin{bmatrix}\mathbf{Y}^{1},\ldots,\mathbf{Y}^{n}\end{bmatrix}^{\top}\in{\mathbb{R}}^{n\times
p}$

the matrix containing all target measurements. Since we do not have access to
the feature paths $x^{i}$ but only to the discrete measurements
$\mathbf{X}^{i}$, we compute the signature of its linear interpolation
normalized by its total variation, sampled up to final time $t^{i}_{k_{i}}$,
and denote by $\mathbf{S}^{\mathcal{D}}_{N}\in\mathbb{R}^{n\times s_{d}(N)}$
the matrix of stacked signatures. Note that signatures of piecewise linear
functions are fast to compute with packages such as signatory (Kidger & Lyons,
2020) or iisignature (Reizenstein & Graham, 2020). The complexity to compute
the signature truncated at order $N$ of the $i$th feature time series is of
the order $\mathcal{O}(\\#D^{i}d^{N})$. Finally, we approximate
$\theta_{N}^{\ast}\in\mathbb{R}^{s_{d}(N)\times p}$ by solving the
optimisation problem

$\displaystyle\min_{\theta\in\mathbb{R}^{s_{d}(N)\times
p}}\frac{1}{2n}\left\lVert\mathbf{Y}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{\textnormal{F}}^{2}+\Omega(\theta),$
(6)

where $\Omega:\mathbb{R}^{s_{d}(N)\times p}\to{\mathbb{R}}^{+}$ is a
regularization term and $\left\lVert\cdot\right\rVert_{\textnormal{F}}$ is the
Frobenius norm. We have reduced the complex problem of learning the solution
map $\Psi$ to a simple penalized linear regression in the signature space.
This linear model on the signature is close to the one studied by Levin et al.
(2013); Lyons et al. (2014); Fermanian (2022).

#### Multiple target measurements.

We also cover the case when $\\#\bar{D}^{i}>1$, that is, the target is
measured at multiple times for every individual, as in the example of
hemorrhage detection. In this case, we have
$\mathbf{Y}^{i}\in{\mathbb{R}}^{m_{i}\times p}$ and we stack the different
measurement matrices $\mathbf{Y}^{i}$ to obtain a matrix $\mathbf{Y}$ of size
$M\times p$, where $M=m_{1}+\dots+m_{n}$. For any $i=1,\dots,n$ and every
$t\in\bar{D}^{i}$, we predict $Y_{t}^{i}$ using the signature of the linear
interpolation of the normalized $(X^{i}_{0},\dots,X^{i}_{t})$. In this manner,
we will be able to predict $Y_{t}^{i}$ at every point where $X^{i}_{t}$ is
sampled. The exact workflow of our model is described in Figure 1.

##  3 Theoretical guarantees

###  3.1 Mathematical setup

We consider a general multiple target measurements setting. To simplify the
exposure of our results, we consider an univariate target path, i.e., $p=1$.
In this case, the true parameter $\theta_{N}^{*}$ is a vector of size
$s_{d}(N)$ and not a matrix. The general case $p\geq 1$, which our algorithm
handles as running $p$ Lasso regressions in parallel, is considered in the
Appendix B, were we prove all results in this setting. In addition, to lighten
the presentation of the oracle inequality, we also focus in this section on
the case of $\omega$-Lipschitz feature paths, that is, for every $i=1,\dots,n$
and for all $s,t\in[0,1]$, $\left\lVert
x^{i}_{t}-x^{i}_{s}\right\rVert\leq\omega\lvert t-s\rvert$. We stress that our
results are valid for continuous paths of bounded variations. We let
$\mathbf{y}\in\mathbb{R}^{M}$ be the matrix collecting all unobserved values
of the target paths at measurement times such that
$\mathbf{y}=\mathbb{E}\big{(}\mathbf{Y}\big{)},$ where the expectation is
taken over the noises $\varepsilon^{i}_{t}$, and define
$\widehat{\theta}_{N,M}$ as

$\displaystyle\widehat{\theta}_{N,M}\in\operatorname*{arg\,min}\limits_{\theta\in\mathbb{R}^{s_{d}(N)}}\frac{1}{2M}\left\lVert\mathbf{Y}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{2}^{2}+\Omega(\theta).$
(7)

For $\delta\in(0,1)$, we define the set

$\displaystyle
A_{\xi}(\delta)=\Big{\\{}\max\limits\left\lVert\xi^{i}_{t}\right\rVert\leq
v_{\xi}\sqrt{d}+v_{\xi}\sqrt{\frac{1}{C}\log\frac{\\#\mathcal{D}}{\delta}}\Big{\\}}$
(8)

where the maximum is taken on all $i=1,\dots,n$ and $t\in D^{i}$, and $C$ is a
universal constant. This set is of probability greater than $1-\delta$ under
Assumption 5 (see Appendix A.5). Similarly, for $k\geq 0$ and
$\bar{\delta}\in(0,1)$, let

$C_{k}(\bar{\delta})=\sqrt{v_{\varepsilon}\log(2Nd^{k}/\bar{\delta})}$

and define

$\displaystyle A_{\varepsilon}(\bar{\delta})$
$\displaystyle=\bigcap_{k=0}^{N}\Big{\\{}\big{\|}\bm{\varepsilon}^{\top}\mathbf{S}^{\mathcal{D}}_{\cdot,[k]}\big{\|}_{\infty}\leq\frac{M^{\frac{1}{2}}C_{k}(\bar{\delta})}{k!}\Big{\\}},$

where $\mathbf{S}^{\mathcal{D}}_{\cdot,[k]}$ is the sub-matrix of size
$M\times d^{k}$ of $\mathbf{S}^{\mathcal{D}}_{N}$ associated to the signature
coefficients of order $k$, and $\bm{\varepsilon}\in{\mathbb{R}}^{M}$ is a
vector of i.i.d. noise terms satisfying Assumption 6 (see Appendix B.2). Under
Assumptions 5 and 6, $A_{\xi}(\delta)\cap A_{\varepsilon}(\bar{\delta})$ is of
probability at least $(1-\delta)(1-\bar{\delta})$. Let

$\Omega(\theta)=\sum\limits_{k=0}^{N}\frac{C_{k}(\bar{\delta})}{k!\sqrt{M}}\big{\|}\theta_{[k]}\big{\|}_{1},$
(9)

where $\theta_{[k]}$ is the subvector of size $d^{k}$ that collects all
elements of $\theta$ associated to words of size $k$ (see Appendix B.2). This
penalization can be implemented by rescaling the feature matrix
$\mathbf{S}_{N}^{\mathcal{D}}$ and solving a standard $\ell_{1}$-penalized
regression problem (see Appendix C.1). Our result extends to more general
penalties by adapting existing techniques from Chesneau & Hebiri (2008) for
the group-lasso, or Lederer et al. (2019) for the hierarchical lasso.

###  3.2 Main results

The error made when learning $\theta_{N}^{\ast}$ by $\widehat{\theta}_{N,M}$
comes from three different sources.

  1. 1.

Truncating the signature used in the regression at depth $N\geq 1$ results in
a truncation bias.

  2. 2.

Discretization of the feature path and the noise affecting each measurement
point induce a discretization error. In particular, there is a trade-off
between sampling frequency and variance of the noise.

  3. 3.

The measurement error on $\mathbf{y}^{i}$ and the finite-sample setting induce
a classical estimation error.

Table 1: Performance of SigLasso, GRU and Neural CDE in different simulation settings, averaged over 10 iterations. In every setting, $n=50$, $\\#\bar{D}^{i}=5$ for all $i=1,\dots,n$ (and therefore $M=250$). |  $L_{2}$ error | MSE on last point
---|---|---
Setting | SigLasso | GRU | Neural CDE | SigLasso | GRU | Neural CDE
Well-specified | 0.13 $\pm$ 0.07 |  1.05 $\pm$ 0.42 |  0.61 $\pm$ 0.38 | 0.73 $\pm$ 0.56 |  3.32 $\pm$ 1.60 |  1.46 $\pm$ 1.20
Ill-specified | 0.15 $\pm$ 0.02 |  0.24 $\pm$ 0.11 |  0.29 $\pm$ 0.15 | 0.09 $\pm$ 0.05 |  0.19 $\pm$ 0.09 |  0.22 $\pm$ 0.15
OU | 0.01 $\pm$ 0.02 |  0.05 $\pm$ 0.06 |  0.17 $\pm$ 0.12 |  0.018 $\pm$ 0.025 |  0.014 $\pm$ 0.020 | 0.013 $\pm$ 0.016
Tumor growth | 0.16 $\pm$ 0.02 |  0.66 $\pm$ 0.09 |  5.29 $\pm$ 1.38 | 0.35 $\pm$ 0.12 |  2.00 $\pm$ 0.38 |  8.76 $\pm$ 9.26

The following lemmas bound each of those errors. We first bound the variance
of the estimator with arguments borrowed from Bickel et al. (2009).

######  Lemma 3.1.

Under Assumptions 1 and 2, on the set $A_{\varepsilon}(\bar{\delta})\cap
A_{\xi}(\delta)$, the prediction error

$\frac{1}{2M}\big{\|}\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{N,M}\big{\|}_{2}^{2}$

is bounded above by

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta^{\ast}_{N}\right\rVert_{2}^{2}+\frac{2C_{N}(\bar{\delta})}{\sqrt{M}}\sum_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})}{k!}.$

See Appendix B.2 for a proof. This inequality decomposes the error into a bias
and a variance term. We denote the signature matrix of the unobserved paths
$x^{i}$ by $\mathbf{S}_{N}\in\mathbb{R}^{M\times s_{d}(N)}$.

As for the bias term, notice that one can write

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta^{\ast}_{N}\right\rVert_{2}^{2}\leq\frac{1}{M}\underbrace{\left\lVert\mathbf{y}-\mathbf{S}_{N}\theta^{\ast}_{N}\right\rVert_{2}^{2}}_{\text{Truncation
bias}}+\frac{1}{M}\underbrace{\left\lVert\mathbf{S}_{N}\theta^{\ast}_{N}-\mathbf{S}_{N}^{\mathcal{D}}\theta^{\ast}_{N}\right\rVert_{2}^{2}}_{\text{Discretization
error}}.$

Bounding each of these terms corresponds to, respectively, Lemmas 3.2 and 3.3.
We stress that the truncation bias is of a different nature than the
discretization error since it depends on a choice of hyperparameter while the
latter is inherent to the data at hand.

######  Lemma 3.2.

Under Assumptions 1 and 2, for any $N\geq 1$,

$\displaystyle\frac{1}{M}\left\lVert\mathbf{y}-\mathbf{S}_{N}\theta^{*}_{N}\right\rVert^{2}_{2}\leq\Bigg{(}\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\Bigg{)}^{2}.$

Under Assumption 3, the right-hand-side decays exponentially fast with $N$.
This lemma is an immediate consequence of Fermanian et al. (2021) (see
Appendix B.3). We now turn to the error induced by the discretization of the
feature path.

######  Lemma 3.3.

Under Assumptions 1, 4, and 5, on the set $A_{\xi}(\delta)$, one has

$\frac{1}{M}\left\lVert(\mathbf{S}_{N}-\mathbf{S}^{\mathcal{D}}_{N})\theta^{\ast}_{N}\right\rVert^{2}_{2}\leq
C_{\mathcal{D},N}(\delta)\sum\limits_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})^{2}}{k!^{2}},$

where $C_{\mathcal{D},N}(\delta)$ is equal to

$\displaystyle(2eLN!)^{2}\Bigg{(}\omega\lvert\mathcal{D}\rvert+v_{\xi}\sqrt{d}+v_{\xi}\sqrt{\frac{1}{C}\log\frac{\\#\mathcal{D}}{\delta}}\Bigg{)}^{2}$

and $C$ is a universal constant.

This lemma relies on a fine analysis of the distance between two signature
layers. The dependence of $C_{\mathcal{D},N}(\delta)$ on sampling mechanisms
and noise is of particular interest. First, the term $\omega\lvert D\rvert$
refers to the longest time between two observations amongst individuals. Not
sampling an individual during a long period of time causes a loss in
information, which is bounded by the Lipschitz control $\omega$ of the feature
path.

The remaining part of $C_{\mathcal{D},N}(\delta)$ is a consequence of the
noises $\xi^{i}_{t}$ affecting the measurement points of the feature time
series and does not vanish with $n$. The emergence of such a bias is a well
studied phenomenon in errors-in-variable models, and cannot be corrected
without precise knowledge of the noise’s variance (Loh & Wainwright, 2011).
Note that if we assume that the feature paths are measured without noise
(i.e., $v_{\xi}=0$), the right-hand-side collapses to a discretization error
term proportional to $\lvert\mathcal{D}\rvert^{2}$ and tends to $0$ as $\lvert
D\rvert\to 0$. The proof of this Lemma is given in Appendix B.5.

From Lemmas 3.1 to 3.3 and the definitions of $A_{\xi}$ and $A_{\varepsilon}$
finally we get the following oracle inequality. The proof is given in Appendix
B.6.

######  Theorem 3.4 (An oracle inequality for learning with signatures).

Under Assumptions 1, 2, 4, 5, and 6, with probability at least
$(1-\delta)(1-\bar{\delta})$, the prediction error is bounded above by

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{N,M}\right\rVert_{2}^{2}$
$\displaystyle\leq\Bigg{(}\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\Bigg{)}^{2}$
(10)
$\displaystyle+C_{\mathcal{D},N}(\delta)\sum\limits_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})^{2}}{k!^{2}}$
(11)
$\displaystyle+\frac{2C_{N}(\bar{\delta})}{\sqrt{M}}\sum_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})}{k!}.$
(12)

All terms depend on the regularity of the vector field $\mathbf{F}$ via the
constants $\Lambda_{k}(\mathbf{F})$: the bigger these constants, the faster
the vector field $\mathbf{F}$ may vary, making the CDE harder to predict. The
convergence speed in $1/\sqrt{M}$ is classical. Also note that our bound is
non-asymptotic and is valid for any $M\geq 1$.

The dependence of the bound on $N$ is highly non-trivial and requires an in-
depth analysis of the regularity of $\mathbf{F}$ in order to bound
$\Lambda_{k}(\mathbf{F})$, which is out of the scope of this paper. The
asymptotic behaviour of this oracle inequality is discussed in Appendix B.7.

##  4 Experiments

We study the performance of SigLasso obtained by solving the optimization
problem (7), where $\Omega(\theta)$ is defined by Equation (9). All details
are given in Appendix C.

###  4.1 Simulations

We consider several settings of data generation. First, in the well-specified
setting, the data is generated from a model with regular feature paths $x$
(piecewise polynomials) and target paths $y$ wich are solutions to the CDE
$dy_{t}=\text{tanh}(Ay_{t})dx_{t}$, where $A$ is a randomly drawn matrix. In
the ill-specified setting, the target $y$ is equal to
$y_{t}=\log\|\sum_{h=1}^{10}x_{t-h}\|$ for any $t\in[0,1]$. In the third
setting, called OU setting, the feature paths are realizations of Brownian
motions and the target paths are Ornstein-Uhlenbeck processes (Borodin &
Salminen, 2012) driven by the feature paths. The last setting corresponds to
the tumor growth model from Simeoni et al. (2004). The feature path represents
the concentration of a treatment drug, generated as the squared value of the
smooth paths used in the well-specified setting, and the target path $y$, the
weight of the tumor, is governed by a system of differential equations given
in Appendix C.6.

Figure 2: $L_{2}$ reconstruction error of SigLasso, GRU and Neural CDE in the
well-specified setting, for varying number of target samples.

We compare SigLasso to a GRU and Neural CDE (Kidger et al., 2020). We measure
the performance of the models with two metrics on a test set: the mean squared
error for predicting the last observation point of the target paths and the
$L_{2}$ error for predicting the full path on a fine grid.

The results are shown in Table 1 and Figure 2. In Figure 2 we consider the
well-specified setting and vary the number of sampling points of the target
paths between 1 and 20. In Table 1 it is fixed to 5 but the simulation
settings change. Sampling of both the target and the feature time series is
highly irregular. SigLasso outperforms Neural CDE and GRU models in
generalizing from a few learning points of the target to its full trajectory
in all settings. An additional byproduct of SigLasso’s simple form is its
training speed: its is approximately 10 times faster than GRU and 100 times
faster than Neural CDE, including cross-validation to select $N$ and
regularization strength (see Appendix C.8).

###  4.2 Forecasting the growth rate of hospitalizations in France during the
Covid-19 pandemic

Figure 3: Interpolation (left of dotted line) and prediction (right of dotted
line) of HGR in region Île de France for SigLasso. The lighter the blue, the
smaller the horizon $h$.

Forecasting hospitalizations in real time during the Covid-19 pandemic is a
notably difficult task. In this experiment, we train our model to learn the
dynamics linking population data related to mobility, vaccination, and
weather, and the hospitalization growth rate (HGR) in each of the 9
metropolitan regions of France based on the data of Paireau et al. (2022). The
feature time series is specific to each region and $12$-dimensional. We
consider prediction horizons $h=1,\dots,14$, meaning that we predict the
hospital saturation at time $t$ using the history of the feature time series
up to $t-h$. Using our notations, we have for each region $d=12$, $p=1$,
$n=1$. The target is sampled every day during the training period left of the
dotted line in Figure 3. The model is then asked to predict the HGR on the
days right of the dotted line.

Both GRU and SigLasso learn smooth and precise dynamics, which generalize well
above the learning horizon and yield similar prediction performance (see
Appendix C.7). Neural CDE and the original method proposed by Paireau et al.
(2022) perform poorly. Figure 3 shows an example of reconstruction and
prediction of HGR obtained with SigLasso.

##  5 Conclusion

We have introduced a novel CDE-based model for interacting systems. Drawing on
the theory of signatures, we derive an oracle bound that depends explicitly on
the roughness of the data sampling. We illustrate the high performance of our
approach on synthetic and real-world data.

The obtained theoretical guarantees rely on strong regularity assumptions on
the vector field $\mathbf{F}$. The exact approximation properties of this
class of vector field are a very interesting direction for future work.
Considering other penalties that take into account the underlying structure of
$\theta_{N}^{*}$ would also be an interesting extension of our work.

## References

  * Bickel et al. (2009) Bickel, P. J., Ritov, Y., and Tsybakov, A. B.  Simultaneous analysis of lasso and dantzig selector.  _The Annals of Statistics_ , 37(4):1705–1732, 2009.
  * Borodin & Salminen (2012) Borodin, A. N. and Salminen, P.  _Handbook of Brownian Motion-Facts and Formulae_.  Birkhäuser, Basel, 2012.
  * Brunton et al. (2020) Brunton, S. L., Noack, B. R., and Koumoutsakos, P.  Machine learning for fluid mechanics.  _Annual review of fluid mechanics_ , 52:477–508, 2020.
  * Buehler et al. (2020) Buehler, H., Horvath, B., Lyons, T., Arribas, I. P., and Wood, B.  A data-driven market simulator for small data environments.  _arXiv preprint arXiv:2006.14498_ , 2020.
  * Chen (1958) Chen, K.-T.  Integration of paths—a faithful representation of paths by non-commutative formal power series.  _Transactions of the American Mathematical Society_ , 89:395–407, 1958.
  * Chen et al. (2018) Chen, R. T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K.  Neural ordinary differential equations.  In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 31, pp. 6572–6583. Curran Associates, Inc., 2018.
  * Chesneau & Hebiri (2008) Chesneau, C. and Hebiri, M.  Some theoretical results on the grouped variables lasso.  _Mathematical Methods of Statistics_ , 17(4):317–326, 2008.
  * De Brouwer et al. (2019) De Brouwer, E., Simm, J., Arany, A., and Moreau, Y.  GRU-ODE-Bayes: Continuous modeling of sporadically-observed time series.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32, pp. 7379–7390. Curran Associates, Inc., 2019.
  * Fattahi et al. (2019) Fattahi, S., Matni, N., and Sojoudi, S.  Learning sparse dynamical systems from a single sample trajectory.  In _2019 IEEE 58th Conference on Decision and Control (CDC)_ , pp. 2682–2689. IEEE, 2019.
  * Fermanian (2021) Fermanian, A.  Embedding and learning with signatures.  _Computational Statistics & Data Analysis_, 157:107148, 2021.
  * Fermanian (2022) Fermanian, A.  Functional linear regression with truncated signatures.  _Journal of Multivariate Analysis_ , 192:105031, 2022.
  * Fermanian et al. (2021) Fermanian, A., Marion, P., Vert, J.-P., and Biau, G.  Framing RNN as a kernel method: A neural ODE approach.  In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), _Advances in Neural Information Processing Systems_ , volume 34, pp. 3121–3134. Curran Associates, Inc., 2021.
  * Friz & Victoir (2008) Friz, P. and Victoir, N.  Euler estimates for rough differential equations.  _Journal of Differential Equations_ , 244:388–412, 2008\.
  * Friz & Victoir (2010) Friz, P. K. and Victoir, N. B.  _Multidimensional Stochastic Processes as Rough Paths: Theory and Applications_ , volume 120 of _Cambridge Studies in Advanced Mathematics_.  Cambridge University Press, Cambridge, 2010.
  * Gehring et al. (2017) Gehring, J., Auli, M., Grangier, D., Yarats, D., and Dauphin, Y. N.  Convolutional sequence to sequence learning.  In Precup, D. and Teh, Y. W. (eds.), _Proceedings of the 34th International Conference on Machine Learning_ , volume 70, pp. 1243–1252. PMLR, 2017.
  * Greydanus et al. (2019) Greydanus, S., Dzamba, M., and Yosinski, J.  Hamiltonian neural networks.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Han et al. (2018) Han, J., Jentzen, A., and E, W.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ , 115(34):8505–8510, 2018.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long short-term memory.  _Neural Computation_ , 9:1735–1780, 1997.
  * Kidger (2022) Kidger, P.  On neural differential equations.  _arXiv preprint arXiv:2202.02435_ , 2022.
  * Kidger & Lyons (2020) Kidger, P. and Lyons, T.  Signatory: differentiable computations of the signature and logsignature transforms, on both cpu and gpu.  In _International Conference on Learning Representations_ , 2020.
  * Kidger et al. (2019) Kidger, P., Bonnier, P., Perez Arribas, I., Salvi, C., and Lyons, T.  Deep signature transforms.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32, pp. 3099–3109. Curran Associates, Inc., 2019.
  * Kidger et al. (2020) Kidger, P., Morrill, J., Foster, J., and Lyons, T.  Neural controlled differential equations for irregular time series.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M. F., and Lin, H. (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 6696–6707. Curran Associates, Inc., 2020.
  * Lagaris et al. (1998) Lagaris, I. E., Likas, A., and Fotiadis, D. I.  Artificial neural networks for solving ordinary and partial differential equations.  _IEEE transactions on neural networks_ , 9(5):987–1000, 1998.
  * Lederer et al. (2019) Lederer, J., Yu, L., and Gaynanova, I.  Oracle inequalities for high-dimensional prediction.  _Bernoulli_ , 25(2):1225–1255, 2019.
  * Lemercier et al. (2021) Lemercier, M., Salvi, C., Damoulas, T., Bonilla, E., and Lyons, T.  Distribution regression for sequential data.  In Banerjee, A. and Fukumizu, K. (eds.), _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130, pp. 3754–3762. PMLR, 2021.
  * Levin et al. (2013) Levin, D., Lyons, T., and Ni, H.  Learning from the past, predicting the statistics for the future, learning an evolving system.  _arXiv preprint arXiv:1309.0260_ , 2013.
  * Loh & Wainwright (2011) Loh, P.-l. and Wainwright, M. J.  High-dimensional regression with noisy and missing data: Provable guarantees with non-convexity.  In Shawe-Taylor, J., Zemel, R., Bartlett, P., Pereira, F., and Weinberger, K. (eds.), _Advances in Neural Information Processing Systems_ , volume 24. Curran Associates, Inc., 2011.
  * Long et al. (2018) Long, Z., Lu, Y., Ma, X., and Dong, B.  PDE-net: Learning PDEs from data.  In Dy, J. and Krause, A. (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80, pp. 3208–3216. PMLR, 2018.
  * Lyons & McLeod (2022) Lyons, T. and McLeod, A. D.  Signature methods in machine learning.  _arXiv preprint arXiv:2206.14674_ , 2022.
  * Lyons et al. (2007) Lyons, T., Caruana, M., and Lévy, T.  _Differential Equations driven by Rough Paths_ , volume 1908 of _Lecture Notes in Mathematics_.  Springer, Berlin, 2007.
  * Lyons et al. (2014) Lyons, T., Ni, H., and Oberhauser, H.  A feature set for streams and an application to high-frequency financial tick data.  In _Proceedings of the 2014 International Conference on Big Data Science and Computing_ , pp. 5. ACM, 2014.
  * Marion et al. (2022) Marion, P., Fermanian, A., Biau, G., and Vert, J.-P.  Scaling resnets in the large-depth regime.  _arXiv preprint arXiv:2206.06929_ , 2022.
  * Marx & Eilers (1999) Marx, B. D. and Eilers, P. H.  Generalized linear regression on sampled signals and curves: a P-spline approach.  _Technometrics_ , 41:1–13, 1999.
  * Morrill et al. (2020a) Morrill, J., Fermanian, A., Kidger, P., and Lyons, T.  A generalised signature method for multivariate time series feature extraction.  _arXiv preprint arXiv:2006.00873_ , 2020a.
  * Morrill et al. (2021a) Morrill, J., Kidger, P., Yang, L., and Lyons, T.  Neural controlled differential equations for online prediction tasks.  _arXiv preprint arXiv:2106.11028_ , 2021a.
  * Morrill et al. (2021b) Morrill, J., Salvi, C., Kidger, P., and Foster, J.  Neural rough differential equations for long time series.  In Meila, M. and Zhang, T. (eds.), _Proceedings of the 38th International Conference on Machine Learning_ , volume 139, pp. 7829–7838. PMLR, 2021b.
  * Morrill et al. (2020b) Morrill, J. H., Kormilitzin, A., Nevado-Holgado, A. J., Swaminathan, S., Howison, S. D., and Lyons, T. J.  Utilization of the signature method to identify the early onset of sepsis from multivariate physiological time series in critical care monitoring.  _Critical Care Medicine_ , 48:e976–e981, 2020b.
  * Morris (2015) Morris, J. S.  Functional regression.  _Annual Review of Statistics and Its Application_ , 2:321–359, 2015.
  * Paireau et al. (2022) Paireau, J., Andronico, A., Hozé, N., Layan, M., Crepey, P., Roumagnac, A., Lavielle, M., Boëlle, P.-Y., and Cauchemez, S.  An ensemble model based on early predictors to forecast covid-19 health care demand in france.  _Proceedings of the National Academy of Sciences_ , 119(18):e2103302119, 2022.
  * Papavasiliou & Ladroue (2011) Papavasiliou, A. and Ladroue, C.  Parameter estimation for rough differential equations.  _The Annals of Statistics_ , 39(4):2047–2073, 2011.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Ramsay & Dalzell (1991) Ramsay, J. O. and Dalzell, C.  Some tools for functional data analysis.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 53:539–561, 1991.
  * Ramsay & Silverman (2005) Ramsay, J. O. and Silverman, B. W.  _Functional Data Analysis. 2nd Edition._ Springer, New York, 2005.
  * Reizenstein & Graham (2020) Reizenstein, J. F. and Graham, B.  Algorithm 1004: The iisignature library: Efficient calculation of iterated-integral signatures and log signatures.  _ACM Transactions on Mathematical Software_ , 46(1):1–21, 2020.
  * Rubanova et al. (2019) Rubanova, Y., Chen, R. T. Q., and Duvenaud, D. K.  Latent ordinary differential equations for irregularly-sampled time series.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32, pp. 5320–5330. Curran Associates, Inc., 2019.
  * Simeoni et al. (2004) Simeoni, M., Magni, P., Cammia, C., De Nicolao, G., Croci, V., Pesenti, E., Germani, M., Poggesi, I., and Rocchetti, M.  Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents.  _Cancer research_ , 64(3):1094–1101, 2004.
  * Sutskever et al. (2014) Sutskever, I., Vinyals, O., and Le, Q. V.  Sequence to sequence learning with neural networks.  In Ghahramani, Z., Welling, M., Cortes, C., Lawrence, N., and Weinberger, K. (eds.), _Advances in Neural Information Processing Systems_ , volume 27. Curran Associates, Inc., 2014.
  * Vershynin (2010) Vershynin, R.  Introduction to the non-asymptotic analysis of random matrices.  _arXiv preprint arXiv:1011.3027_ , 2010.
  * Vershynin (2018) Vershynin, R.  _High-dimensional probability: An introduction with applications in data science_ , volume 47.  Cambridge university press, 2018.
  * Wang et al. (2020) Wang, B., Wu, Y., Taylor, N., Lyons, T., Liakata, M., Nevado-Holgado, A. J., and Saunders, K. E.  Learning to detect bipolar disorder and borderline personality disorder with language and speech in non-clinical interviews.  _Interspeech 2020_ , pp. 437–441, 2020.
  * Willard et al. (2020) Willard, J., Jia, X., Xu, S., Steinbach, M., and Kumar, V.  Integrating physics-based modeling with machine learning: A survey.  _arXiv preprint arXiv:2003.04919_ , 2020.
  * Yang et al. (2022) Yang, W., Lyons, T., Ni, H., Schmid, C., and Jin, L.  Developing the path signature methodology and its application to landmark-based human action recognition.  In Yin, G. and Zariphopoulou, T. (eds.), _Stochastic Analysis, Filtering, and Stochastic Optimization: A Commemorative Volume to Honor Mark H.A. Davis’s Contributions_ , pp. 431–464. Springer, 2022.
  * Zubov et al. (2021) Zubov, K., McCarthy, Z., Ma, Y., Calisto, F., Pagliarino, V., Azeglio, S., Bottero, L., Luján, E., Sulzer, V., Bharambe, A., et al.  NeuralPDE: Automating physics-informed neural networks (PINNs) with error approximations.  _arXiv preprint arXiv:2107.09443_ , 2021.
  * Zwillinger (1989) Zwillinger, D.  _Handbook of differential equations_.  Academic Press, Inc., London, 1989.

Supplementary Material

##  Appendix A Mathematical details

###  A.1 The Riemann-Stieltjes integral

We fitst recall two key properties on the Riemann-Stieljes integral. For a
general presentation of the Riemann-Stieltjes integral, we refer to Friz &
Victoir (2010).

######  Proposition A.1.

Let $x\in C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$ and
$y:[0,1]\rightarrow\mathbb{R}^{d}$ be a continuous path. Then

$\displaystyle\left\lVert\int_{0}^{t}y_{s}dx_{s}\right\rVert\leq\left\lVert
y\right\rVert_{\infty,[0,t]}\left\lVert
x\right\rVert_{\textnormal{1-var},[0,t]}$

We refer the reader to Friz & Victoir (2010, Proposition 2.2) for a proof.

######  Proposition A.2 (Integration by parts).

Let $x,y\in C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$. Then

$\displaystyle\int_{s}^{t}y_{u}dx_{u}+\int_{s}^{t}x_{u}dy_{u}=y_{t}x_{t}-y_{s}x_{s}$

See Friz & Victoir (2010, Proposition 2.4) for a proof.

###  A.2 The truncated tensor algebra

This section introduces notations and definitions on the space in which
signatures are defined, namely, the tensor algebra. While for the exposition
of our main results, the truncated signature of a path $x\in
C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$ at depth $N\geq 1$ can be
assimilated to an element of $\mathbb{R}^{s_{d}(N)}$, it is often useful to
place ourselves in the tensor algebra to obtain finer bounds or technical
results.

Let $x\in C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$ be a path of
bounded variation. For a word
$I=(i_{1},\dots,i_{k})\in\left\\{1,\dots,d\right\\}^{k}$ of size $k$, the
signature coefficient $S^{I}(x_{[0,1]})$ can be seen as an element of the
$k$-th tensor product of $\mathbb{R}^{d}$ with itself, denoted by
$\big{(}\mathbb{R}^{d}\big{)}^{\otimes k}$. For instance, the coefficients of
order $k=1$ can be written as a vector and the coefficients of order $k=2$ as
a matrix, and so on, i.e.,

$\displaystyle\mathbb{X}_{[0,1]}^{1}=\begin{bmatrix}\int_{0}^{1}dx_{s}^{(1)}\\\
\vdots\\\
\int_{0}^{1}dx_{s}^{(d)}\end{bmatrix}\quad\text{and}\quad\mathbb{X}_{[0,1]}^{2}=\begin{bmatrix}\int_{0}^{1}dx_{s}^{(1)}dx_{s}^{(1)}&\dots&\int_{0}^{1}dx_{s}^{(1)}dx_{s}^{(d)}\\\
\vdots&&\vdots\\\
\int_{0}^{1}dx_{s}^{(d)}dx_{s}^{(1)}&\dots&\int_{0}^{1}dx_{s}^{(d)}dx_{s}^{(d)}\end{bmatrix}.$

We now define a norm on $(\mathbb{R}^{d})^{\otimes k}$. Let
$a\in(\mathbb{R}^{d})^{\otimes k}$ and $(e_{1},\dots,e_{d})$ be the canonical
basis of $\mathbb{R}^{d}$. Then $(e_{i_{1}}\otimes\dots\otimes
e_{i_{k}})_{(i_{1},\dots,i_{k})\in\\{1,\dots,d\\}^{k}}$ is a basis of
$(\mathbb{R}^{d})^{\otimes k}$. We can thus write $a$ as
$a=(a^{I})_{I\in\left\\{1,\dots,d\right\\}^{k}}$. For every $k\geq 0$, the
vector space $(\mathbb{R}^{d})^{\otimes k}$ is naturally endowed with the norm

$\left\lVert a\right\rVert_{(\mathbb{R}^{d})^{\otimes
k}}^{2}=\sum\limits_{I\in\left\\{1,\dots,d\right\\}^{k}}\big{(}a^{I}\big{)}^{2}.$

Remark that this norm satisfies for any $x\in(\mathbb{R}^{d})^{\otimes k}$ and
$y\in(\mathbb{R}^{d})^{\otimes m}$,

$\displaystyle\left\lVert x\otimes
y\right\rVert_{(\mathbb{R}^{d})^{\otimes(k+m)}}=\left\lVert
x\right\rVert_{(\mathbb{R}^{d})^{\otimes k}}\left\lVert
y\right\rVert_{(\mathbb{R}^{d})^{\otimes m}}.$ (13)

We refer to Fermanian et al. (2021) for further deails. The signature
truncated at depth $N\geq 1$ collects elements from
$\mathbb{R},\big{(}\mathbb{R}^{d}\big{)}^{\otimes
2},\dots,\big{(}\mathbb{R}^{d}\big{)}^{\otimes N}$. It can thus be seen as an
element of the truncated tensor algebra

$T_{N}(\mathbb{R}^{d})=\mathbb{R}\oplus\big{(}\mathbb{R}^{d}\big{)}^{\otimes
2}\oplus\dots\oplus\big{(}\mathbb{R}^{d}\big{)}^{\otimes N}.$

Let $a=(a_{0},\dots,a_{N})\in T_{N}(\mathbb{R}^{d})$, where every
$a_{k}\in(\mathbb{R}^{d})^{\otimes k}$. We define the norm

$\displaystyle\left\lVert
a\right\rVert_{T_{N}(\mathbb{R}^{d})}=\Big{(}\sum\limits_{k=0}^{N}\left\lVert
a_{k}\right\rVert_{(\mathbb{R}^{d})^{\otimes k}}^{2}\Big{)}^{1/2}.$

To clarify, if we consider the truncated signature of $x$ at depth $N\geq 1$,
which is an element of $T_{N}(\mathbb{R}^{d})$, then

$\displaystyle\left\lVert
S_{N}(x_{[0,t]})\right\rVert_{T_{N}(\mathbb{R}^{d})}=\Big{(}\sum\limits_{k=0}^{N}\left\lVert\mathbb{X}_{[0,t]}^{k}\right\rVert^{2}_{(\mathbb{R}^{d})^{\otimes
k}}\Big{)}^{1/2}=\Big{(}\sum\limits_{k=0}^{N}\quad\sum\limits_{I\in\left\\{1,\dots,d\right\\}^{k}}S^{I}(x_{[0,t]})^{2}\Big{)}^{1/2}.$

Note that this norm is exactly equivalent to the Euclidian norm of
$\mathbb{R}^{s_{d}(N)}$, which is the space we consider in the exposition of
our main results for the sake of simplicity.

We are now ready to define the tensor product on the truncated tensor algebra.
For two elements $a=(a_{0},\dots,a_{N})$ and $b=(b_{0},\dots,b_{N})$ both in
$T_{N}(\mathbb{R}^{d})$, we define

$a\otimes b=(c_{0},\dots,c_{j},\dots,c_{N}),\quad\text{ where
}c_{j}=\sum\limits_{k=0}^{j}a_{k}\otimes b_{j-k}.$

For any $k=0,\dots,N$, we let
$\pi_{k}:T_{N}(\mathbb{R}^{d})\to(\mathbb{R}^{d})^{\otimes k}$ be the
canonical projection of $T_{N}(\mathbb{R}^{d})$ onto
$(\mathbb{R}^{d})^{\otimes k}$. More precisely, for every
$a=(a_{0},\dots,a_{N})\in T_{N}(\mathbb{R}^{d})$,

$\pi_{k}(a)=a_{k}$

We also define the canonical projection $\Pi_{k}:T_{N}(\mathbb{R}^{d})\to
T_{k}(\mathbb{R}^{d})$ defined by

$\Pi_{k}(a)=\left(a_{0},\dots,a_{k}\right).$

###  A.3 The differential product

We first define the differential product, in order to give a precise statement
of Assumption 3.

######  Definition A.3.

Let $F,G:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$ be two smooth vector fields,
i.e., each of their components is $\mathcal{C}^{\infty}$. Denote by $J(\cdot)$
the Jacobian matrix. The differential product $F\star
G:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$ is the smooth vector field defined
for any $h\in\mathbb{R}^{p}$

$\displaystyle(F\star G)(h)=\sum\limits_{j=1}^{e}\frac{\partial G}{\partial
h_{j}}(h)F_{j}(h)=J(G)(h)F(h).$

The differential product is not associative. We therefore use the convention
to evaluate it from right to left, that is,

$\displaystyle F^{1}\star F^{2}\star F^{3}=F^{1}\star\left(F^{2}\star
F^{3}\right).$

Let $F:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p\times d}$ be a smooth vector
field. We write $F^{1},\dots,F^{d}$ the columns of $\mathbf{F}$. Every
$F^{i}$, for $i=1,\dots,d$, can thus be seen as a map from $\mathbb{R}^{p}$ to
$\mathbb{R}^{p}$. Recall that $y_{0}\in\mathbb{R}^{p}$ is the initial
condition of the CDE defined in Assumption 2. Let
$I=(i_{1},\dots,i_{k})\in\left\\{1,\dots,d\right\\}^{k}$. We now define

$\displaystyle\Phi^{I}_{\mathbf{F}}(y_{0})=\left(F^{i_{1}}\star\dots\star
F^{i_{k}}\right)(y_{0})\in\mathbb{R}^{p}.$

We refer to Fermanian et al. (2021) for greater details on the differential
product. We now define for all $k\geq 1$

$\displaystyle\Lambda_{k}(\mathbf{F})=\sup\limits_{1\leq i_{1},\dots,i_{k}\leq
d}\left\lVert\Phi^{I}_{\mathbf{F}}(y_{0})\right\rVert\in\mathbb{R},$ (14)

and use the convention $\Lambda_{0}(\mathbf{F})=\left\lVert
y_{0}\right\rVert$. Remark that Assumption 3 implies that

$\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\underset{N\to+\infty}{\longrightarrow}0.$

As an immediate consequence, the truncation bias

$\Bigg{(}\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\Bigg{)}^{2}$

introduced in Lemma 3.2, vanishes as $N$ grows.

###  A.4 Analogy with the Taylor extension

The definition of the Taylor expansion of a CDE exposed in the previous
subsection is technical. However, it can simply be thought of as a
generalization of the the classical Taylor expansion of a
$\mathcal{C}^{\infty}$ function $f:\mathbb{R}\rightarrow\mathbb{R}$. Recall
that in this case, the Taylor expansion at $0$ evaluated at $t\in\mathbb{R}$
of order $N\in\mathbb{N}$ writes as a power series

$\displaystyle f(t)\approx
f(0)+\frac{f^{\prime}(0)}{1!}t+\dots+\frac{f^{(N)}(0)}{N!}t^{N}.$ (15)

Every element of this power series is a product of two terms: the derivatives
of $f$ encode some information about the regularity of $f$ at the initial
point $0$ and do not depend on $t$, while the polynomial terms $t^{k}$ allow
this linearized form to evolve with time $t$. Remark that these polynomial
terms do not depend on $f$.

Similarly, Equation (3) is also a sum of products of two terms. On the one
hand, the evolving nature of the system, instead of being handled by the
polynomial terms $t^{k}$, are now captured by the signature coefficients
$S^{I}(x_{[0,t]})$. As the polynomial terms, they do not depend on
$\mathbf{F}$. On the other hand, the information about the initial value of
the system at time $t=0$ and the dynamics of $\mathbf{F}$ are summarized by
the differential product $\Phi^{I}_{\mathbf{F}}\left(y_{0}\right)$, which play
the same role as the successive derivatives in Equation (15). To capture the
multivariate nature of the paths, the Taylor expansion is summed over multi-
indexes, or words, $I=(i_{1},\dots,i_{k})\in\\{1,\dots,d\\}^{k}$ of size $k$
for $k\in\mathbb{N}$.

###  A.5 Properties of subgaussian random vectors

We start with the definition of a subgaussian random variable, see Vershynin
(2010) for more details.

######  Definition A.4.

A real-valued random variable $X$ is said to be $\sigma^{2}$-subgaussian if
for all $t>0$

$\mathbb{P}(X>t)\leq\exp(-t^{2}/\sigma^{2}),$

or, equivalently, if for all $t\in\mathbb{R}$

$\mathbb{E}(e^{tX})\leq\exp(-ct^{2}\sigma^{2}),$

where $c$ is an universal constant. A random vector $Z$ is subgaussian if, for
any vector $c$ of norm $1$, $\langle Z,c\rangle$ is subgaussian.

The norm of a sequence of $d$ subgaussian random variables concentrates around
$\sqrt{d}$, as stated by the following lemma.

######  Lemma A.5.

Let $X_{1},\dots,X_{n}$ be a sequence of i.i.d. $\sigma^{2}$-subgaussian
random variables. Let $X=(X_{1},\dots,X_{d})\in\mathbb{R}^{d}$. There exists a
universal constant $C$ such that for all $t>0$

$\mathbb{P}(\left\lVert X\right\rVert_{2}\geq
t+\sigma\sqrt{d})\leq\exp(-Ct^{2}/\sigma^{2}).$

###### Proof.

We refer to Vershynin (2018, Theorem 3.1.1) for a proof. ∎

We can use this lemma to bound the maximum of $n$ sequences of $d$ subgaussian
random variables with high probability.

######  Lemma A.6.

Let $X_{1},\dots,X_{n}$ be a sequence of i.i.d. $\sigma^{2}$-subgaussian
random variables, such that for all $i=1,\dots,n$,
$X_{i}=(X_{i1},\dots,X_{id})$. Then there exists a universal constant such
that for all $\delta\in(0,1)$

$\mathbb{P}\Big{(}\max\limits_{i=1,\dots,n}\left\lVert
X_{i}\right\rVert\leq\sigma\sqrt{d}+\sigma\sqrt{\frac{1}{C}\log(n/\delta)}\Big{)}\geq
1-\delta.$

###### Proof.

Using Lemma A.5 and a union bound, we have

$\displaystyle\mathbb{P}\Big{(}\max\limits_{i=1,\dots,n}\left\lVert
X_{i}\right\rVert\geq\sigma\sqrt{d}+\sigma\sqrt{\frac{1}{C}\log(n/\delta)}\Big{)}$
$\displaystyle=\mathbb{P}\Big{(}\bigcup\limits_{i=1}^{n}\Big{\\{}\left\lVert
X_{i}\right\rVert_{2}\geq\sigma\sqrt{d}+\sigma\sqrt{\frac{1}{C}\log(n/\delta)}\Big{\\}}\Big{)}$
$\displaystyle\quad\leq\sum\limits_{i=1}^{n}\mathbb{P}\Big{(}\left\lVert
X_{i}\right\rVert_{2}\geq\sigma\sqrt{d}+\sigma\sqrt{\frac{1}{C}\log(n/\delta)}\Big{)}$
$\displaystyle\quad\leq\delta,$

which yields the desired inequality. ∎

Notice that the universal constant is identical between both lemmas. As a
consequence of this last lemma, under Assumption 5, the set

$\displaystyle A_{\xi}(\delta)=\Big{\\{}\max\limits_{i=1,\dots,n,t\in
D^{i}}\left\lVert\xi^{i}_{t}\right\rVert\leq
v_{\xi}\sqrt{d}+v_{\xi}\sqrt{\frac{1}{C}\log(\\#\mathcal{D}/\delta)}\Big{\\}}$
(16)

where $\\#\mathcal{D}=\sum\limits_{i=1}^{n}\\#D^{i}$ is of probability at
least $1-\delta$.

We also need the following lemma.

######  Lemma A.7.

Let $X_{1},\dots,X_{n}$ be a sequence of i.i.d. $\sigma^{2}$-subgaussian
random variables. Let $Z_{1},\dots,Z_{n}$ be random variables such that for
all $i=1,\dots,n$, $|Z_{i}|\leq\alpha$ almost surely. Then
$\sum\limits_{i=1}^{n}X_{i}Z_{i}$ is $n\sigma^{2}\alpha^{2}$-subgaussian.

###### Proof.

We use the characterization of subgaussian random variables by their
characteristic function. For all $t>0$,

$\displaystyle\mathbb{E}\bigg{[}e^{t\sum_{i=1}^{n}X_{i}Z_{i}}\bigg{]}=\mathbb{E}\bigg{[}\prod_{i=1}^{n}\mathbb{E}\big{[}e^{tX_{i}Z_{i}}\,|\,Z_{1},\dots,Z_{n}\big{]}\bigg{]}\leq\mathbb{E}\bigg{[}\prod_{i=1}^{n}\mathbb{E}\big{[}e^{tX_{i}\alpha}\big{]}\bigg{]}\leq\mathbb{E}\big{[}e^{ct^{2}n\alpha^{2}\sigma^{2}}\big{]}.$

This finally yields that

$\displaystyle\mathbb{E}\bigg{[}e^{t\sum
X_{i}Z_{i}}\bigg{]}\leq\mathbb{E}\big{[}e^{ct^{2}n\alpha^{2}\sigma^{2}}\big{]},$

which concludes the proof. ∎

##  Appendix B Proofs

###  B.1 Preliminary notations

Let $(E,\left\lVert\cdot\right\rVert_{E})$ be a normed vector space and
$x:[0,1]\to E$. The supremum norm of $x$ is defined for all $t\in[0,1]$ as

$\left\lVert x\right\rVert_{\infty,[0,t]}=\sup_{s\in[0,t]}\left\lVert
x_{s}\right\rVert_{E}.$

When referring to the total variation $\left\lVert
x\right\rVert_{\textnormal{1-var},[0,1]}$ of a path $x:[0,1]\to\mathbb{R}^{d}$
over the whole domain, depending on the mathematical context, we will
sometimes drop the time subscript and simply write $\left\lVert
x\right\rVert_{\textnormal{1-var}}$.

When referring to a matrix $A=(A_{ij})\in\mathbb{R}^{n\times p}$, we define
classicaly the infinite and Frobenius norms by

$\left\lVert A\right\rVert_{\infty}=\max_{\begin{subarray}{c}i=1,\dots,n\\\
j=1,\dots,p\end{subarray}}\lvert A_{ij}\rvert\quad\text{and}\quad\left\lVert
A\right\rVert_{F}=\sqrt{\sum\limits_{\begin{subarray}{c}i=1,\dots,n\\\
j=1,\dots,p\end{subarray}}\lvert A_{ij}\rvert^{2}}.$

We now introduce some notations to take advantage of the structure of
$\theta^{*}_{N}$. The true parameter of the Taylor expansion of the model CDE,
defined in Equation (4), can be written in block notation as

$\displaystyle\theta^{*}_{N}=\left[\begin{array}[]{ccc}\begin{matrix}\theta^{*}_{[0],1}&\cdots&\theta^{*}_{[0],p}\\\
\hline\cr&&\\\ \theta^{*}_{[1],1}&\cdots&\theta^{*}_{[1],p}\\\ &&\\\
\hline\cr\\\ &&\\\ &&\\\ \theta^{*}_{[2],1}&\cdots&\theta^{*}_{[2],p}\\\ &&\\\
&&\\\ \hline\cr&&\\\ &\vdots&\\\ &&\\\ \hline\cr&&\\\ &&\\\ &&\\\
\theta^{*}_{[N],1}&\cdots&\theta^{*}_{[N],p}\\\ &&\\\ &&\\\ &&\\\
\end{matrix}\end{array}\right]\in\mathbb{R}^{s_{d}(N)\times
p},\quad\text{where}\quad\theta^{\ast}_{[k],\ell}\in{\mathbb{R}}^{d^{k}\times
1},k=0,\dots,N,\,\ell=1,\dots,p.$ (18)

Every column of $\theta_{N}^{*}$ corresponds to a dimension of the target,
while blocks of lines correspond to signatures layers. Thus for every
$k=0,\dots,N$ and $\ell=1,\dots,p$, $\theta^{*}_{[k],\ell}$ is a column vector
of size $d^{k}$.

Similarly, for a general $\theta\in\mathbb{R}^{s_{d}(N)\times p}$ and the
SigLasso estimator $\widehat{\theta}_{N,M}$, we will refere to the blocks
forming these matrices as respectively $\theta_{[k],\ell}$ and
$\widehat{\theta}_{[k],\ell}$, for $k=0,\dots,N$ and $\ell=1,\dots p$.

Likewise, the signature feature matrix
$\mathbf{S}_{N}^{\mathcal{D}}\in\mathbb{R}^{M\times s_{d}(N)}$ can be written
in block notation as

$\mathbf{S}_{N}^{\mathcal{D}}=\begin{bmatrix}\,1&\vline&\mathbf{S}^{\mathcal{D}}_{\cdot,[1]}&\vline&\mathbf{S}^{\mathcal{D}}_{\cdot,[2]}&\vline&\cdots&\vline&\mathbf{S}^{\mathcal{D}}_{\cdot,[N]}\end{bmatrix}=\setcounter{MaxMatrixCols}{16}\begin{bmatrix}\,1&\vline&\mathbf{S}^{\mathcal{D}}_{1,[1]}&\vline&&\mathbf{S}^{\mathcal{D}}_{1,[2]}&&\vline&&\vline&&&\mathbf{S}^{\mathcal{D}}_{1,[N]}&&&\\\
\,\vdots&\vline&\vdots&\vline&&\vdots&&\vline&\cdots&\vline&&&\vdots&&&\\\
\,1&\vline&\mathbf{S}^{\mathcal{D}}_{n,[1]}&\vline&&\mathbf{S}^{\mathcal{D}}_{n,[2]}&&\vline&&\vline&&&\mathbf{S}^{\mathcal{D}}_{n,[N]}&&\end{bmatrix},$

where for any $k=1,\dots,N$,
$\mathbf{S}^{\mathcal{D}}_{\cdot,[k]}\in{\mathbb{R}}^{M\times d^{k}}$ and, for
every individual $i=1,\dots,n$,
$\mathbf{S}^{\mathcal{D}}_{i,[k]}\in{\mathbb{R}}^{m_{i}\times d^{k}}$ (recall
that $m_{i}$ is the number of measurements of the target path $y^{i}$). More
precisely, given her target sampling grid
$\bar{D}^{i}=(\bar{t}_{1}^{i},\dots,\bar{t}^{i}_{m_{i}})$, the individual-
specific signature block of depth $k$ is equal to

$\mathbf{S}^{\mathcal{D}}_{i,[k]}=\begin{bmatrix}1&S^{(1)}(X^{i}_{[0,\bar{t}^{i}_{1}]})&\cdots&S^{(d,\dots,d)}(X^{i}_{[0,\bar{t}^{i}_{1}]})\\\
\vdots&\vdots&&\vdots\\\
1&S^{(1)}(X^{i}_{[0,\bar{t}^{i}_{m_{i}}]})&\cdots&S^{(d,\dots,d)}(X^{i}_{[0,\bar{t}^{i}_{m_{i}}]})\\\
\end{bmatrix},$

where the path $t\to X^{i}_{t}$ is a linear interpolation of the observed time
series $\mathbf{X}^{i}$. The same notations will be used for the true
signature feature matrix $\mathbf{S}_{N}$. We use the bracket notation
$[\cdot]$ both in $\theta^{*}_{N}$ and $\mathbf{S}_{N}^{\mathcal{D}}$ to
emphasise that both the columns of the feature matrix and the lines of learned
parameter correspond to words of the alphabet $\\{1,\dots,d\\}$.

The unobserved matrix of true values of the target writes as

$\displaystyle\mathbf{y}=\begin{bmatrix}\mathbf{y}^{1}\\\ \vdots\\\
\mathbf{y}^{n}\end{bmatrix}=\begin{bmatrix}\mathbf{y}^{1}_{1}&\cdots&\mathbf{y}^{1}_{p}\\\
\vdots&&\vdots\\\
\mathbf{y}^{n}_{1}&\cdots&\mathbf{y}^{n}_{p}\end{bmatrix}=\begin{bmatrix}y^{1}_{1,\bar{t}^{1}_{1}}&\cdots&y^{1}_{p,\bar{t}^{1}_{1}}\\\
\vdots&&\vdots\\\
y^{1}_{1,\bar{t}^{1}_{m_{1}}}&\cdots&y^{1}_{p,\bar{t}^{1}_{m_{1}}}\\\
\vdots&&\vdots\\\
y^{n}_{1,\bar{t}^{n}_{1}}&\cdots&y^{n}_{p,\bar{t}^{n}_{1}}\\\
\vdots&&\vdots\\\
y^{n}_{1,\bar{t}^{n}_{m_{n}}}&\cdots&y^{n}_{p,\bar{t}^{n}_{m_{n}}}\end{bmatrix}\in\mathbb{R}^{M\times
p}$ (19)

and the measurement matrix $\mathbf{Y}\in\mathbb{R}^{M\times p}$ can be
written in a similar fashion as

$\displaystyle\mathbf{Y}=\begin{bmatrix}\mathbf{Y}^{1}\\\ \vdots\\\
\mathbf{Y}^{n}\end{bmatrix}=\begin{bmatrix}\mathbf{Y}^{1}_{1}&\cdots&\mathbf{Y}^{1}_{p}\\\
\vdots&&\vdots\\\
\mathbf{Y}^{n}_{1}&\cdots&\mathbf{Y}^{n}_{p}\end{bmatrix}=\begin{bmatrix}Y^{1}_{1,\bar{t}^{1}_{1}}&\cdots&Y^{1}_{p,\bar{t}^{1}_{1}}\\\
\vdots&&\vdots\\\
Y^{1}_{1,\bar{t}^{1}_{m_{1}}}&\cdots&Y^{1}_{p,\bar{t}^{1}_{m_{1}}}\\\
\vdots&&\vdots\\\
Y^{n}_{1,\bar{t}^{n}_{1}}&\cdots&Y^{n}_{p,\bar{t}^{n}_{1}}\\\
\vdots&&\vdots\\\
Y^{n}_{1,\bar{t}^{n}_{m_{n}}}&\cdots&Y^{n}_{p,\bar{t}^{n}_{m_{n}}}\end{bmatrix}=\begin{bmatrix}y^{1}_{1,\bar{t}^{1}_{1}}+\varepsilon^{1}_{1,\bar{t}^{1}_{1}}&\cdots&y^{1}_{p,\bar{t}^{1}_{1}}+\varepsilon^{1}_{p,\bar{t}^{1}_{1}}\\\
\vdots&&\vdots\\\
y^{1}_{1,\bar{t}^{1}_{m_{1}}}+\varepsilon^{1}_{1,\bar{t}^{1}_{m_{1}}}&\cdots&y^{1}_{p,\bar{t}^{1}_{m_{1}}}+\varepsilon^{1}_{p,\bar{t}^{1}_{m_{1}}}\\\
\vdots&&\vdots\\\
y^{n}_{1,\bar{t}^{n}_{1}}+\varepsilon^{n}_{1,\bar{t}^{n}_{1}}&\cdots&y^{n}_{p,\bar{t}^{n}_{1}}+\varepsilon^{n}_{p,\bar{t}^{n}_{1}}\\\
\vdots&&\vdots\\\
y^{n}_{1,\bar{t}^{n}_{m_{n}}}+\varepsilon^{n}_{1,\bar{t}^{n}_{m_{n}}}&\cdots&y^{n}_{p,\bar{t}^{n}_{m_{n}}}+\varepsilon^{n}_{p,\bar{t}^{n}_{m_{n}}}\end{bmatrix}$
(20)

###  B.2 Proof of Lemma 3.1

Using the definition of $\Lambda_{k}$, we get the following proposition which
allows to obtain an explicit dependence of the oracle bound on the regularity
of $\mathbf{F}$.

######  Proposition B.1.

Let $\theta^{\ast}_{N}$ be defined as in Equation (4). Then

$\left\lVert\theta^{*}_{N}\right\rVert_{\textnormal{F}}^{2}\leq\sum\limits_{k=0}^{N}d^{k}\Lambda_{k}(\mathbf{F})^{2},$

and, for all $k=0,\dots,N$ and $\ell=1,\dots,p$,

$\big{\|}\theta^{*}_{[k],\ell}\big{\|}_{1}\leq d^{k}\Lambda_{k}(\mathbf{F}).$

###### Proof.

By definition,

$\left\lVert\theta^{*}_{N}\right\rVert_{\textnormal{F}}^{2}=\sum\limits_{k=0}^{N}\quad\sum\limits_{1\leq
i_{1},\dots,i_{k}\leq d}\left\lVert F^{i_{1}}\star\dots\star
F^{i_{k}}(y_{0})\right\rVert_{2}^{2}.$

Since for all $(i_{1},\dots,i_{k})\in\left\\{1,\dots,d\right\\}^{k}$,

$\left\lVert F^{i_{1}}\star\dots\star
F^{i_{k}}(y_{0})\right\rVert_{2}^{2}\leq\Lambda_{k}(\mathbf{F})^{2},$

we get

$\displaystyle\sum\limits_{k=0}^{N}\quad\sum\limits_{1\leq
i_{1},\dots,i_{k}\leq d}\left\lVert F^{i_{1}}\star\dots\star
F^{i_{k}}(y_{0})\right\rVert_{2}^{2}$
$\displaystyle\leq\sum\limits_{k=0}^{N}\quad\sum\limits_{1\leq
i_{1},\dots,i_{k}\leq d}\Lambda_{k}(\mathbf{F})^{2}$
$\displaystyle\leq\sum\limits_{k=0}^{N}d^{k}\Lambda_{k}(\mathbf{F})^{2}.$

We now turn to the second inequality. For $k=0$, the inequality holds by
definition. For $k=1,\dots,N$ and $\ell=1,\dots,p$, by definition of the
$\ell_{1}$ norm,

$\displaystyle\left\lVert\theta^{*}_{[k],\cdot}\right\rVert_{1}=\sum\limits_{1\leq
i_{1},\dots,i_{k}\leq
d}\left\lVert\Phi^{I}_{\mathbf{F}}(y_{0})\right\rVert_{1},$

This yields

$\left\lVert\theta^{*}_{[k],\cdot}\right\rVert_{1}\leq
d^{k}\Lambda_{k}(\mathbf{F})$

and thus

$\left\lVert\theta^{*}_{[k],\ell}\right\rVert_{1}\leq
d^{k}\Lambda_{k}(\mathbf{F})$

for $\ell=1,\dots,p$. ∎

The following lemma is needed to leverage classical proof techniques to bound
the prediction error of the Lasso estimator.

######  Lemma B.2.

Let $x\in C^{\textnormal{1-var}}_{L}([0,1],\mathbb{R}^{d})$. Then
conditionally on $A_{\xi}(\delta)$, for a given signature layer $k\geq 1$, the
maximum among all signature coefficients and individuals is bounded from
above, that is

$\left\lVert\mathbf{S}^{\mathcal{D}}_{\cdot,[k]}\right\rVert_{\infty}\leq\frac{1}{k!}.$

###### Proof.

It is well known (see, e.g., Fermanian, 2022, Proposition 3) that if
$\mathbb{X}^{k}$ is the signature of a path $x\in
C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$, then

$\displaystyle\left\lVert\mathbb{X}^{k}\right\rVert_{(\mathbb{R}^{d})^{\otimes
k}}\leq\frac{\left\lVert x\right\rVert_{\textnormal{1-var}}^{k}}{k!}.$

As a consequence, for every word $I$ of size $k$, one gets

$\displaystyle\big{|}S^{I}(x))\big{|}\leq\frac{\left\lVert
x\right\rVert_{\textnormal{1-var}}^{k}}{k!}.$

The matrix $\mathbf{S}_{N}^{\mathcal{D}}$ is constructed by taking signatures
of linear interpolations of the $\mathbf{X}^{i}$s normalized by their total
variation. It therefore contains only signatures of paths of total variation
bounded by 1. Taking the maximum on $I\in\\{1,\dots,d\\}^{k}$ and individuals
$i=1,\dots,n$, we get

$\left\lVert\mathbf{S}^{\mathcal{D}}_{\cdot,[k]}\right\rVert_{\infty}\leq\frac{1}{k!}.$

∎

This final inequality being stated, we can now go back to the proof of Lemma
3.1. We prove it in full generality for $p\geq 1$. In this proof, we make
extensive use of the notations introduced in Subsection B.1 and refer the
reader to it if a notation is unclear.

###### Proof.

In all the proof, we place ourselves on the set $A_{\xi}(\delta)$ defined by
Equation (8), which ensures that the matrix $\mathbf{S}_{N}^{\mathcal{D}}$,
seen as a random quantity, is well defined. Recall that we have two sources of
randomness: the feature noises $\xi^{i}_{t}$ on the $\mathbf{X}^{i}$s and the
target noises $\varepsilon^{i}_{t}$ on the $\mathbf{Y}^{i}$s. The feature
noises appear only in $\mathbf{S}_{N}^{\mathcal{D}}$ and make it a random
quantity. For $\mathbf{S}_{N}^{\mathcal{D}}$ to be well-defined, we then need
the total variation of the linear interpolation of the feature time series
$\mathbf{X}^{i}$ to be finite. This holds on the set $A_{\xi}(\delta)$ since
all noises are then bounded.

Recall that we have defined $\widehat{\theta}_{N,M}$ as

$\widehat{\theta}_{N,M}\in\operatorname*{arg\,min}\limits_{\theta\in\mathbb{R}^{s_{d}(N)\times
p}}\frac{1}{2M}\left\lVert\mathbf{Y}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{\textnormal{F}}^{2}+\Omega(\theta).$

Note that

$\frac{1}{2M}\left\lVert\mathbf{Y}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{\textnormal{F}}^{2}+\Omega(\theta)=\sum_{\ell=1}^{p}\frac{1}{2M}\left\lVert\mathbf{Y}_{\ell}-\mathbf{S}_{N}^{\mathcal{D}}\theta_{[\cdot],\ell}\right\rVert_{2}^{2}+\Omega(\theta_{[\cdot],\ell})$

where $\mathbf{Y}_{\ell}\in\mathbb{R}^{M}$ is the $\ell$-th column of the
target measurement matrix defined in Equation (20). The quantity
$\theta_{[\cdot],\ell}\in\mathbb{R}^{s_{d}(N)}$ is the $\ell$-th column of the
parameter matrix defined in Equation (18).

By definition, for any $\theta\in{\mathbb{R}}^{s_{d}(N)}$, we have

$\displaystyle\left\lVert\mathbf{Y}_{\ell}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{[\cdot],\ell}\right\rVert_{2}^{2}\leq\left\lVert\mathbf{Y}_{\ell}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{2}^{2}+\Omega(\theta)-\Omega(\widehat{\theta}_{[\cdot],\ell}).$

Moreover, letting
$\bm{\varepsilon}_{\ell}=(\varepsilon^{1}_{\ell,\bar{t}^{1}_{1}},\dots,\varepsilon^{n}_{\ell,\bar{t}^{n}_{m_{n}}})^{\top}\in\mathbb{R}^{M}$
be a vector of i.i.d. noises (see Equation (20)), we have
$\mathbf{Y}_{\ell}=\mathbf{y}_{\ell}+\bm{\varepsilon}_{\ell}$. The Pythagorean
theorem then yields for any $\theta\in{\mathbb{R}}^{s_{d}(N)}$,

$\displaystyle\left\lVert\mathbf{Y}_{\ell}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{2}^{2}=\left\lVert\mathbf{y}_{\ell}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{2}^{2}+\left\lVert\bm{\varepsilon}_{\ell}\right\rVert^{2}+2\langle\bm{\varepsilon}_{\ell},y_{\ell}-\mathbf{S}_{N}^{\mathcal{D}}\theta\rangle.$

Putting these two equations together, we obtain

$\frac{1}{2M}\left\lVert\mathbf{y}_{\ell}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{[\cdot],\ell}\right\rVert_{2}^{2}\leq\frac{1}{2M}\left\lVert\mathbf{y}_{\ell}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{2}^{2}+\frac{1}{M}\langle\bm{\varepsilon}_{\ell},\mathbf{S}_{N}^{\mathcal{D}}(\widehat{\theta}_{[\cdot],\ell}-\theta)\rangle+\Omega(\theta)-\Omega(\widehat{\theta}_{[\cdot],\ell}).$
(21)

We now work at each layer of the signature matrix
$\mathbf{S}_{N}^{\mathcal{D}}$. Towards that end, we rewrite

$\mathbf{S}_{N}^{\mathcal{D}}\big{(}\widehat{\theta}_{[\cdot],\ell}-\theta\big{)}=\sum_{k=0}^{N}\mathbf{S}_{\cdot,[k]}^{\mathcal{D}}\big{(}\widehat{\theta}_{[k],\ell}-\theta_{[k]}\big{)},$

and bound

$\big{\langle}\bm{\varepsilon}_{\ell},\mathbf{S}_{N}^{\mathcal{D}}(\widehat{\theta}_{[\cdot],\ell}-\theta_{\ell})\big{\rangle}=\sum_{k=0}^{N}\big{\langle}\bm{\varepsilon}_{\ell},\mathbf{S}_{\cdot,[k]}^{\mathcal{D}}(\widehat{\theta}_{[k],\ell)}-\theta_{[k]})\big{\rangle}\leq\sum_{k=0}^{N}\|\bm{\varepsilon}_{\ell}^{\top}\mathbf{S}_{\cdot,[k]}^{\mathcal{D}}\|_{\infty}\|\widehat{\theta}_{[k],\ell}-\theta_{[k]}\|_{1}$

by $\ell_{1}-\ell_{\infty}$ norms duality. We fix $k$ and study the term
$\|\bm{\varepsilon}_{\ell}^{\top}\mathbf{S}_{\cdot,[k]}^{\mathcal{D}}\|_{\infty}$.
Lemma B.2 ensures that each of the words of the signature layer of depth $k$
is bounded by $1/k!$. As a consequence, by Lemma A.7, under Assumption 6,
every element of the vector
$\bm{\varepsilon}_{\ell}^{\top}\mathbf{S}_{\cdot,[k]}^{\mathcal{D}}$ is
$v_{\varepsilon}M/k!^{2}$-subgaussian. It follows that, for any real number
$\mu>0$,

$\mathbb{P}\Big{(}\|\bm{\varepsilon}_{\ell}^{\top}\mathbf{S}_{\cdot,[k]}^{\mathcal{D}}\|_{\infty}>\mu\Big{)}\leq
2d^{k}\exp\Big{(}-\frac{(k!)^{2}\mu^{2}}{v_{\varepsilon}M}\Big{)}.$

We furthermore place ourselves on $A_{\varepsilon}(\bar{\delta})$ defined by

$\displaystyle
A_{\varepsilon}(\bar{\delta})=\bigcap_{\ell=1}^{p}\bigcap_{k=0}^{N}\Big{\\{}\|\bm{\varepsilon}_{\ell}^{\top}\mathbf{S}_{\cdot,[k]}^{\mathcal{D}}\|_{\infty}\leq\frac{1}{k!}\sqrt{v_{\varepsilon}M\log(2pNd^{k}/\bar{\delta})}\Big{\\}}.$

We have just seen that, under Assumption 6 (and still conditionally on
$A_{\xi}(\delta)$), one has $\mathbb{P}(A_{\varepsilon}(\bar{\delta}))\geq
1-\bar{\delta}$. Putting together all terms in Equation (21) and plugging the
definition of $\Omega$ given in Equation (9), we obtain that, on the set
$A_{\varepsilon}(\bar{\delta})\cap A_{\xi}(\delta)$, for all
$\theta\in\mathbb{R}^{s_{d}(N)\times p}$,

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{N,M}\right\rVert_{\textnormal{F}}^{2}$
$\displaystyle\leq\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{\textnormal{F}}^{2}$
$\displaystyle\quad+\sum_{\ell=1}^{p}\sum_{k=0}^{N}\Big{(}\frac{1}{M}\|\bm{\varepsilon}_{\ell}^{\top}\mathbf{S}_{\cdot,[k]}^{\mathcal{D}}\|_{\infty}\|\widehat{\theta}_{[k],\ell}-\theta_{[k],\ell}\|_{1}+\frac{C_{k}(\bar{\delta})}{k!\sqrt{M}}\big{(}\|\theta_{[k],\ell}\|_{1}-\|\widehat{\theta}_{[k],\ell}\|_{1}\big{)}\Big{)}$
$\displaystyle\leq\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{\textnormal{F}}^{2}$
$\displaystyle\quad+\sum_{\ell=1}^{p}\sum_{k=0}^{N}\frac{1}{k!\sqrt{M}}\sqrt{v_{\varepsilon}\log(2pNd^{k}/\bar{\delta})}\big{(}\|\widehat{\theta}_{[k],\ell}-\theta_{[k],\ell}\|_{1}+\|\theta_{[k],\ell}\|_{1}-\|\widehat{\theta}_{[k],\ell}\|_{1}\big{)}.$

Choosing $\theta=\theta_{N}^{\ast}$, by the triangular inequality,

$\displaystyle\|\widehat{\theta}_{[k],\ell}-\theta^{\ast}_{[k],\ell}\|_{1}+\|\theta^{\ast}_{[k],\ell}\|_{1}-\|\widehat{\theta}_{[k],\ell}\|_{1}\leq
2\|\theta^{\ast}_{[k],\ell}\|_{1},$

which finally gives us

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{N,M}\right\rVert_{\textnormal{F}}^{2}$
$\displaystyle\leq\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta_{N}^{\ast}\right\rVert_{\textnormal{F}}^{2}+\frac{2}{\sqrt{M}}\sqrt{v_{\varepsilon}\log(2pNd^{N}/\bar{\delta})}\sum_{\ell=1}^{p}\sum_{k=0}^{N}\frac{\|\theta^{*}_{[k],\ell}\|_{1}}{k!}$
$\displaystyle\leq\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta_{N}^{\ast}\right\rVert_{\textnormal{F}}^{2}+\frac{2p}{\sqrt{M}}\sqrt{v_{\varepsilon}\log(2pNd^{N}/\bar{\delta})}\sum_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})}{k!}$
$\displaystyle=\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta_{N}^{\ast}\right\rVert_{\textnormal{F}}^{2}+\frac{2pC_{N}(\bar{\delta})}{\sqrt{M}}\sum_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})}{k!},$

where the last inequality comes from Proposition B.1. To conclude the proof,
we just need to compute the probability of the set $A_{\xi}(\delta)\cap
A_{\varepsilon}(\bar{\delta})$. It is an immediate consequence of Lemma A.6
that $\mathbb{P}(A_{\xi}(\delta))\geq 1-\delta$, and we have seen that
$\mathbb{P}(A_{\varepsilon}(\bar{\delta})|A_{\xi}(\delta))\geq
1-\bar{\delta}$, which yields that

$\mathbb{P}(A_{\xi}(\delta)\cap
A_{\varepsilon}(\bar{\delta}))\geq(1-\bar{\delta})(1-\delta).$

∎

###  B.3 Proof of Lemma 3.2

This proof relies on bouding the remainder of the Taylor expansion of the CDE.

###### Proof.

For every $i=1,\dots,n$ and a given point $t_{i}\in\bar{D}^{i}$, one has,
using the upper bound of the approximation error of a CDE by its Taylor
expansion provided by Fermanian et al. (2021, Proposition 4)

$\left\lVert
y^{i}_{t_{i}}-S_{N}(x^{i}_{[0,t_{i}]})\theta_{N}^{*}\right\rVert\leq\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}.$

This immediately gives

$\displaystyle\frac{1}{M}\left\lVert\mathbf{y}-\mathbf{S}_{N}\theta_{N}^{*}\right\rVert^{2}_{\textnormal{F}}=\frac{1}{M}\sum\limits_{i=1}^{n}\sum\limits_{t_{i}\in\bar{D}^{i}}\left\lVert
y^{i}_{t_{i}}-S_{N}(x^{i}_{[0,t_{i}]})\theta_{N}^{*}\right\rVert^{2}$
$\displaystyle\leq\frac{1}{M}\sum\limits_{i=1}^{M}\Big{(}\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\Big{)}^{2}=\Big{(}\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\Big{)}^{2},$

which concludes the proof. ∎

###  B.4 A layer-wise bound on the signature

We now prove that signature layers are locally Lipschitz mappings. We start
with the following proposition.

######  Proposition B.3.

Let $x\in C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$. Then for all
$t\in[0,1]$, the path $t\mapsto\mathbb{X}^{k}_{[0,t]}$ has $1$-variation

$\displaystyle\left\lVert\mathbb{X}^{k}\right\rVert_{\textnormal{1-var},[0,t]}\leq\frac{L^{k}}{k!}.$

###### Proof.

By definition of the total variation,

$\left\lVert\mathbb{X}^{k}\right\rVert_{\textnormal{1-var},[0,t]}=\sup_{D}\sum\limits_{i=1}^{m}\left\lVert\mathbb{X}^{k}_{[0,t_{i+1}]}-\mathbb{X}^{k}_{[0,t_{i}]}\right\rVert_{(\mathbb{R}^{d})^{\otimes
k}}=\sup_{D}\sum\limits_{i=1}^{m}\left\lVert\mathbb{X}^{k}_{[t_{i},t_{i+1}]}\right\rVert_{(\mathbb{R}^{d})^{\otimes
k}},$

since $\mathbb{X}^{k}_{[0,t]}=\int_{0}^{t}dx_{u_{1}}\otimes\dots\otimes
dx_{u_{k}}$, and where the supremum is taken over finite dissections
$D=\\{0=t_{1},\dots,t_{m}=1\\}$ of $[0,1]$. Notice that the signature layer of
depth $k$ is here written as an element of $(\mathbb{R}^{d})^{\otimes k}$,
which is more convenient for this proof. Then

$\sup_{D}\sum\limits_{i=1}^{m}\left\lVert\mathbb{X}^{k}_{[t_{i},t_{i+1}]}\right\rVert_{(\mathbb{R}^{d})^{\otimes
k}}\leq\sup_{D}\sum\limits_{i=1}^{m}\frac{\left\lVert
x\right\rVert_{\textnormal{1-var},[t_{i},t_{i+1}]}^{k}}{k!}\leq\frac{1}{k!}\sup_{D}\Big{(}\sum\limits_{i=1}^{m}\left\lVert
x\right\rVert_{\textnormal{1-var},[t_{i},t_{i+1}]}\Big{)}^{k}=\frac{1}{k!}\sup_{D}\left\lVert
x\right\rVert_{\textnormal{1-var},[0,1]}^{k}\leq\frac{L^{k}}{k!},$

where the second inequality follows from the multinomial theorem and the last
equality comes from the fact that for all $s<u<t$, $\left\lVert
x\right\rVert_{\textnormal{1-var},[s,u]}+\left\lVert
x\right\rVert_{\textnormal{1-var},[u,t]}=\left\lVert
x\right\rVert_{\textnormal{1-var},[s,t]}$. This ends our proof. ∎

We now state a bound on the difference between the $k$-th layer of the
signatures of two different paths.

######  Theorem B.4.

Let $x,z\in C_{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$. Then for all $k\geq
2$, the difference in supremum norm between the paths
$t\to\mathbb{X}^{k}_{[0,t]}$ and $t\to\mathbb{Z}^{k}_{[0,t]}$ is bounded by

$\left\lVert\mathbb{X}^{k}-\mathbb{Z}^{k}\right\rVert_{\infty,[0,t]}\leq
2L^{k-1}\sum\limits_{j=1}^{k-1}\frac{1}{j!}\left\lVert
x-z\right\rVert_{\infty,[0,t]}\leq 2eL^{k-1}\left\lVert
x-z\right\rVert_{\infty,[0,t]}$

and

$\left\lVert\mathbb{X}^{1}_{[0,t]}-\mathbb{Z}^{1}_{[0,t]}\right\rVert\leq
2\left\lVert x-z\right\rVert_{\infty,[0,t]}.$

###### Proof.

Our proof works by induction. Let $x,z\in
C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$, and for $t\in[0,1]$ denote
by $\mathbb{X}^{k}_{[0,t]}$ resp. $\mathbb{Z}^{k}_{[0,t]}$ the $k$-th layer of
the signature of $x$ resp. $z$. For $k=1$ and $t\in[0,1]$, remark that

$\mathbb{X}^{1}_{[0,t]}-\mathbb{Z}^{1}_{[0,t]}=\int_{0}^{t}d(x_{u}-z_{u})=x_{t}-z_{t}-(x_{0}-z_{0})$

such that

$\left\lVert\mathbb{X}^{1}_{[0,t]}-\mathbb{Z}^{1}_{[0,t]}\right\rVert\leq\left\lVert
x-z\right\rVert_{\infty,[0,t]}+\left\lVert x_{0}-z_{0}\right\rVert\leq
2\left\lVert x-z\right\rVert_{\infty,[0,t]}.$

Consider now $k\geq 2$. We have

$\mathbb{X}^{k}_{[0,t]}-\mathbb{Z}^{k}_{[0,t]}=\int_{0}^{t}\mathbb{X}^{k-1}_{[0,s]}\otimes
dx_{s}-\int_{0}^{t}\mathbb{Z}^{k-1}_{[0,s]}\otimes
dz_{s}=\int_{0}^{t}\mathbb{X}^{k-1}_{[0,s]}\otimes
d(x_{s}-z_{s}+z_{s})-\int_{0}^{t}\mathbb{Z}^{k-1}_{[0,s]}\otimes dz_{s},$

and thus

$\mathbb{X}^{k}_{[0,t]}-\mathbb{Z}^{k}_{[0,t]}=\int_{0}^{t}\mathbb{X}^{k-1}_{[0,s]}\otimes
d(x_{s}-z_{s})+\int_{0}^{t}\big{(}\mathbb{X}^{k-1}_{[0,s]}-\mathbb{Z}^{k-1}_{[0,s]}\big{)}\otimes
dz_{s}.$

We now bound each of these terms separately. First,

$\left\lVert\int_{0}^{t}\big{(}\mathbb{X}^{k-1}_{[0,s]}-\mathbb{Z}^{k-1}_{[0,s]}\big{)}\otimes
dz_{s}\right\rVert_{(\mathbb{R}^{d})^{\otimes
k}}\leq\left\lVert\mathbb{X}^{k-1}-\mathbb{Z}^{k-1}\right\rVert_{\infty,[0,t]}\left\lVert
z\right\rVert_{\textnormal{1-var},[0,t]}\leq\left\lVert\mathbb{X}^{k-1}-\mathbb{Z}^{k-1}\right\rVert_{\infty,[0,t]}L.$

Moving to the first integral, integration by parts yields

$\int_{0}^{t}\mathbb{X}^{k-1}_{[0,s]}\otimes
d(x_{s}-z_{s})=\mathbb{X}^{k}_{[0,t]}\otimes(x_{t}-z_{t})-\mathbb{X}^{k}_{[0,0]}\otimes(x_{0}-z_{0})-\int_{0}^{t}(x_{s}-z_{s})\otimes
d\mathbb{X}^{k-1}_{[0,s]}.$

We stress that Proposition (A.2) applies since the integral over the tensor
product is taken coordinate-wise. Since $\mathbb{X}^{k-1}_{[0,0]}=0$, we are
left with

$\int_{0}^{t}\mathbb{X}^{k-1}_{[0,s]}\otimes
d(x_{s}-z_{s})=\mathbb{X}^{k}_{[0,t]}\otimes(x_{t}-z_{t})-\int_{0}^{t}(x_{s}-z_{s})\otimes
d\mathbb{X}^{k-1}_{[0,s]}.$

Using Lemma B.3 and submultiplicativity of the tensor norms, this can thus be
bounded by

$\displaystyle\left\lVert\int_{0}^{t}\mathbb{X}^{k-1}_{[0,s]}\otimes
d(x_{s}-z_{s})\right\rVert_{(\mathbb{R}^{d})^{\otimes k}}$
$\displaystyle\leq\left\lVert\mathbb{X}^{k-1}_{[0,t]}\right\rVert_{(\mathbb{R}^{d})^{\otimes(k-1)}}\left\lVert
x-z\right\rVert_{\infty,[0,t]}+\left\lVert
x-z\right\rVert_{\infty,[0,t]}\left\lVert\mathbb{X}^{k-1}\right\rVert_{\textnormal{1-var},[0,t]}$
$\displaystyle=\frac{2L^{k-1}}{(k-1)!}\left\lVert
x-z\right\rVert_{\infty,[0,t]}.$

Finally, we are left with

$\displaystyle\left\lVert\mathbb{X}^{k}-\mathbb{Z}^{k}\right\rVert_{\infty,[0,t]}\leq\frac{2L^{k-1}}{(k-1)!}\left\lVert
x-z\right\rVert_{\infty,[0,t]}+\left\lVert\mathbb{X}^{k-1}-\mathbb{Z}^{k-1}\right\rVert_{\infty,[0,t]}L,$

which can be recursively bounded by

$\left\lVert\mathbb{X}^{k}-\mathbb{Z}^{k}\right\rVert_{\infty,[0,t]}\leq
2L^{k-1}\left\lVert
x-z\right\rVert_{\infty,[0,t]}\sum\limits_{j=1}^{k-1}\frac{1}{j!}\leq
2L^{k-1}e\left\lVert x-z\right\rVert_{\infty,[0,t]}.$

∎

Note that this inequality implies that if $z$ is chosen as the linear
interpolation of a discretization of $x$ on a grid $D$, and if the grid gets
finer, all signature layers converge at speed $\left\lVert
x-z\right\rVert_{\infty,[0,t]}$ but the multiplicative constant increases with
depth (if $L\geq 1$). Figure 4 illustrates this phenomenon.

Figure 4: Difference between the signature of a continuous path $x$ and the
signature of its discretized and noisy counterpart $X$, without noise on the
discretization points (left), with noise of variance $v_{\xi}=0.08^{2}$
(middle) and with noise of variance $v_{\xi}=0.5^{2}$. For every number of
sampling points, we average the distance between the two signature over $50$
randomly chosen discretizations of the interval $[0,1]$. The discretized path
is generated as in the well-specified setting (see Appendix C.4).

###  B.5 Proof of Lemma 3.3

First, recall that for a generic path $x:[0,1]\rightarrow\mathbb{R}^{d}$, a
modulus of continuity is a continuous function $\omega_{x}:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0}$ vanishing at $0$ such that for all
$t,s\in[0,1]$

$\left\lVert x_{t}-x_{s}\right\rVert\leq\omega_{x}(\lvert t-s\rvert).$

Also recall that by Heine’s theorem, we can define such a modulus of
continuity for every continuous mapping $[0,1]$ to $\mathbb{R}^{d}$.

We start by giving a general lemma that bounds the difference between the
signature layers of a path and its discretized version. Its proof is based on
the results of the previous section.

######  Lemma B.5.

Let $x\in C_{L}^{\textnormal{1-var}}([0,1],\mathbb{R}^{d})$ and
$\omega_{x}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ its modulus of
continuity. Let $x^{D}:[0,1]\rightarrow\mathbb{R}^{d}$ be the path obtained by
linear interpolation of the discretization of $x$ on a grid $D$ corrupted by
additive noise $\xi$. Let $\mathbb{X}^{k}_{[0,t]}$ and
$\mathbb{X}^{k,D}_{[0,t]}$ be their respective $k$-th layers of signature.
Then for all $t\in[0,1]$ and $k\geq 2$

$\left\lVert\mathbb{X}^{k}-\mathbb{X}^{k,D}\right\rVert_{\infty,[0,1]}\leq
2L^{k-1}\sum\limits_{j=1}^{k-1}\frac{1}{j!}\Big{(}\max\limits_{0\leq
s\leq\lvert D\rvert}\omega_{x}(s)+\max\limits_{t\in
D}\left\lVert\xi_{t}\right\rVert\Big{)},$

and for $k=1$

$\left\lVert\mathbb{X}^{1}-\mathbb{X}^{1,D}\right\rVert_{\infty,[0,1]}\leq
2\Big{(}\max\limits_{0\leq s\leq\lvert D\rvert}\omega_{x}(s)+\max\limits_{t\in
D}\left\lVert\xi_{t}\right\rVert\Big{)}.$

###### Proof.

Theorem B.4 yields for $k\geq 2$

$\left\lVert\mathbb{X}^{k}-\mathbb{X}^{k,D}\right\rVert_{\infty,[0,1]}\leq
2L^{k-1}\sum\limits_{j=1}^{k-1}\frac{1}{j!}\left\lVert
x-x^{D}\right\rVert_{\infty,[0,t]}$

Now, remark that

$\displaystyle\left\lVert x-x^{D}\right\rVert_{\infty,[0,1]}\leq\left\lVert
x-\tilde{x}\right\rVert_{\infty,[0,1]}+\max\limits_{t\in
D}\left\lVert\xi_{t}\right\rVert$ (22)

from the triangular inequality, where $\tilde{x}$ is the piecewise linear path
obtained by linear interpolation of $x_{0},x_{t_{1}},\dots,x_{t_{j}}$. Now,
since the paths $x$ and $\tilde{x}$ coincide on $0,t_{1},\dots,t_{j}$, we have

$\displaystyle\left\lVert
x-\tilde{x}\right\rVert_{\infty,[0,1]}=\max\limits_{i=0,\dots,j-1}\left\lVert
x-\tilde{x}\right\rVert_{\infty,[t_{i},t_{i+1}]}\leq\max\limits_{i=0,\dots,j-1}\omega_{x}\big{(}\lvert
t_{i+1}-t_{i}\rvert\big{)}=\max\limits_{0\leq s\leq\lvert
D\rvert}\omega_{x}(s).$

This gets us

$\displaystyle\left\lVert\mathbb{X}^{k}-\mathbb{X}^{k,D}\right\rVert_{\infty,[0,1]}\leq
2L^{k-1}\sum\limits_{j=1}^{k-1}\frac{1}{j!}\Big{(}\max\limits_{0\leq
s\leq\lvert D\rvert}\omega_{x}(s)+\max\limits_{t\in
D}\left\lVert\xi_{t}\right\rVert\Big{)}.$

For the case $k=1$, we immediately get

$\left\lVert\mathbb{X}^{1}-\mathbb{X}^{1,D}\right\rVert_{\infty,[0,1]}\leq
2\Big{(}\max\limits_{0\leq s\leq\lvert D\rvert}\omega_{x}(s)+\max\limits_{t\in
D}\left\lVert\xi_{t}\right\rVert\Big{)}$

using the same technique as above. ∎

This result is illustrated in Figure 4. One can notice that as predicted by
our theoretical bounds, the convergence of signature of high order happens at
the same rate than the convergence of signatures of lower order. However, the
multiplicative constant controlling the tightness of the bound increases with
$N$, leading to a slower convergence when $N$ increases. Strong noise hinders
the convergence of the signature of the discretized path since in this case,
the noise’s variance is independent of the number of sampling points : adding
more sampling points means adding more noise. There thus are two trade-offs
when learning with signatures. A first trade-off is between sampling frequency
and order: with paths sampled at low resolution, one should prefer lower order
signatures, which trade model complexity against precise features. A second
trade-off is between sampling and noise: if the feature time series are very
noisy, the precision of the features increases up to a certain point, past
which noise prevails.

With this result in hand, we can now prove Lemma 3.3.

###### Proof.

We restrict ourselves to the $\omega$-Lipschitz case. In this case,
$\omega_{x}(s)=\omega s$, and Lemma B.5 becomes

$\displaystyle\left\lVert\mathbb{X}^{k}-\mathbb{X}^{k,D}\right\rVert_{\infty,[0,1]}\leq
2L^{k-1}\sum\limits_{j=0}^{k-1}\frac{1}{j!}\Big{(}\omega\lvert
D\rvert+\max\limits_{t\in D}\left\lVert\xi_{t}\right\rVert\Big{)}\leq
2eL^{k-1}\Big{(}\omega\lvert D\rvert+\max\limits_{t\in
D}\left\lVert\xi_{t}\right\rVert\Big{)}$

This results holds for a single path on the time interval $[0,t]$. Now moving
to the feature matrices, we have

$\displaystyle\frac{1}{M}\left\lVert(\mathbf{S}_{N}-\mathbf{S}^{\mathcal{D}}_{N})\theta^{\ast}_{N}\right\rVert^{2}_{F}$
$\displaystyle\leq\frac{1}{M}\sum\limits_{i=1}^{n}\sum\limits_{k=0}^{N}\left\lVert(\mathbf{S}_{i,[k]}-\mathbf{S}^{\mathcal{D}}_{i,[k]})\theta^{*}_{[k],\cdot}\right\rVert_{F}^{2}$
$\displaystyle\leq\frac{1}{M}\sum\limits_{i=1}^{n}\sum\limits_{t\in\bar{D}^{i}}\sum\limits_{k=0}^{N}d^{k}\Lambda_{k}(\mathbf{F})^{2}\big{(}2eL^{k-1}\big{(}\omega\lvert
D^{i}\rvert+\max\limits_{t^{\prime}\in D,t^{\prime}\leq
t}\left\lVert\xi_{t^{\prime}}\right\rVert\big{)}\Big{)}^{2}$
$\displaystyle\leq 4e^{2}\big{(}\omega\lvert
D\rvert+\max\limits_{t^{\prime}\in D,t^{\prime}\leq
t}\left\lVert\xi_{t^{\prime}}\right\rVert\big{)}\Big{)}^{2}L^{2}\sum\limits_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})^{2}}{k!^{2}}\times
k!^{2}$ $\displaystyle\leq 4e^{2}N!^{2}\big{(}\omega\lvert
D\rvert+\max\limits_{t^{\prime}\in D,t^{\prime}\leq
t}\left\lVert\xi_{t^{\prime}}\right\rVert\big{)}\Big{)}^{2}L^{2}\sum\limits_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})^{2}}{k!^{2}}.$

On the set $A_{\xi}(\delta)$, one has

$\displaystyle\max\limits_{i=1,\dots,n,t\in
D^{i}}\left\lVert\xi^{i}_{t}\right\rVert\leq
v_{\xi}\sqrt{d}+v_{\xi}\sqrt{C^{-1}\log(\delta^{-1}\\#\mathcal{D})}.$ (23)

Writing

$C_{\mathcal{D},N}(\delta)=4e^{2}L^{2}N!^{2}\Big{(}\omega\lvert
D\rvert+v_{\xi}\sqrt{d}+v_{\xi}\sqrt{C^{-1}\log(\delta^{-1}\\#\mathcal{D})}\Big{)}^{2}$

One finally gets with probability $1-\delta$ that

$\displaystyle\frac{1}{M}\left\lVert(\mathbf{S}_{N}-\mathbf{S}^{\mathcal{D}}_{N})\theta^{\ast}_{N}\right\rVert^{2}_{F}\leq
C_{\mathcal{D},N}(\delta)\sum\limits_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})^{2}}{k!^{2}}.$

∎

###  B.6 Proof of the main Theorem

We finally combine all Lemmas to obtain the desired oracle bound.

###### Proof.

First, we have from Lemma 3.1 that on $A_{\varepsilon}(\bar{\delta})$,

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{N,M}\right\rVert_{\textnormal{F}}^{2}\leq\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta^{\ast}_{N}\right\rVert_{\textnormal{F}}^{2}+\frac{2pC_{N}(\bar{\delta})}{\sqrt{M}}\sum_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})}{k!}.$

The first term of the right-hand side of this inequality is bounded by

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta^{\ast}_{N}\right\rVert_{\textnormal{F}}^{2}\leq\frac{1}{M}\left\lVert\mathbf{y}-\mathbf{S}_{N}\theta^{\ast}_{N}\right\rVert_{\textnormal{F}}^{2}$
$\displaystyle+\frac{1}{M}\left\lVert\mathbf{S}_{N}\theta^{\ast}_{N}-\mathbf{S}_{N}^{\mathcal{D}}\theta^{\ast}_{N}\right\rVert_{\textnormal{F}}^{2}.$

By Lemma 3.2 and Lemma 3.3, this can in turn be bounded on
$A_{\varepsilon}(\bar{\delta})\cap A_{\xi}(\delta)$ by

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\theta^{\ast}_{N}\right\rVert_{\textnormal{F}}^{2}\leq\Bigg{(}\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\Bigg{)}^{2}+C_{\mathcal{D},N}(\delta)\sum\limits_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})^{2}}{k!^{2}}$

Combining all the pieces, this finally gives us, on
$A_{\varepsilon}(\bar{\delta})\cap A_{\xi}(\delta)$,

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{N,M}\right\rVert_{\textnormal{F}}^{2}$
$\displaystyle\leq\Bigg{(}\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\Bigg{)}^{2}$
$\displaystyle+C_{\mathcal{D},N}(\delta)\sum\limits_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})^{2}}{k!^{2}}$
$\displaystyle+\frac{2pC_{N}(\bar{\delta})}{\sqrt{M}}\sum_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})}{k!}.$

∎

###  B.7 Asymptotics

We briefly discuss the asymptotic behaviour of the upper bound of the oracle
inequality.

#### Sampling coarseness $\mathcal{D}$.

In the noiseless case, $C_{\mathcal{D},N}\to 0$ as $\lvert\mathcal{D}\rvert\to
0$. This means that the discretization bias of our estimator vanishes as the
sampling gets finer. However, in the noisy case - that is, when $v_{\xi}>0$ \-
our estimator is durably biased. This is consistent with the literature on
learning with fussy features (Loh & Wainwright, 2011).

#### Truncation depth $N$.

A natural question is whether the bias of our estimator vanishes as
$N\to\infty$. If we have perfect sampling, i.e. the limit case where
$\mathcal{D}=0$ and $v_{\xi}=0$, our bound on the prediction error becomes on
$A_{\varepsilon}(\bar{\delta})$

$\displaystyle\frac{1}{2M}\left\lVert\mathbf{y}-\mathbf{S}_{N}^{\mathcal{D}}\widehat{\theta}_{N,M}\right\rVert_{\textnormal{F}}^{2}$
$\displaystyle\leq\Bigg{(}\frac{d^{N+1}\Lambda_{N+1}(\mathbf{F})}{(N+1)!}\Bigg{)}^{2}$
$\displaystyle+\frac{2pC_{N}(\bar{\delta})}{\sqrt{M}}\sum_{k=0}^{N}\frac{d^{k}\Lambda_{k}(\mathbf{F})}{k!}.$

The first term of this bound vanishes as an immediate consequence of
Assumption 3, while the second term is a statistical error term that behaves
like $\frac{\sqrt{\log(Nd^{N})}}{\sqrt{M}}$. In order to obtain an
assymptotical convergence, we thus need that $N\log(dN)=o(M)$.

In the more realistic setting where $\lvert\mathcal{D}\rvert>0$, the
discretization bias behaves like $L^{N-1}N!\lvert\mathcal{D}\rvert$. It is
thus sufficient to assume that $\lvert D\rvert=o(1/N!)$. If $v_{\xi}>0$, our
estimator is durably biased due to the measurement noise, and this bias
increases with $N\to\infty$. This is due to a ”propagation of chaos”
phenomenon: the difference between the unobserved feature path and the
interpolated feature time series is amplified by taking the successive
iterated integrals that define the signature. This advocates for using simple,
low-order signature models in the presence of noise, as the gain in precision
obtained when taking higher $N$ and reducing the truncation bias will at some
point be lost because of the amplified noise.

#### Dimension $p$ of the target path.

Our oracle bound only depends on $p$ through the statistical error term. This
term is proportional to $p\sqrt{\log p}$, which is expected in multitask
regression.

#### Dimension $d$ of the feature path.

Our oracle bound exhibits multiple dependencies in $d$. First, the truncation
bias grows polynomially with $d$. Similarly, the discretization bias also
depends polynomially on $d$. Finally, the statistical error term is
proportional to $\log d$ times a polynomial term.

##  Appendix C Algorithms, experiments, and supplementary results

###  C.1 Implementation details

Recall that the SigLasso estimator $\widehat{\theta}_{N,M}$ is defined as

$\displaystyle\widehat{\theta}_{N,M}\in\operatorname*{arg\,min}\limits_{\theta\in\mathbb{R}^{s_{d}(N)\times
p}}\frac{1}{2M}\left\lVert\mathbf{Y}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{\textnormal{F}}^{2}+\Omega(\theta),$

where

$\displaystyle\Omega(\theta)=\sum\limits_{k=0}^{N}\frac{M^{\frac{1}{2}}C_{k}(\bar{\delta})}{k!}\left\lVert\theta_{[k],\cdot}\right\rVert_{1},$

and

$\displaystyle
C_{k}(\bar{\delta})=\sqrt{v_{\varepsilon}\log(2pNd^{k}/\bar{\delta})}$

for $\bar{\delta}\in[0,1]$. Our goal is first to rewrite the penalty
$\Omega(\theta)$ as

$\Omega(\theta)=C\sum\limits_{k=0}^{N}\lambda_{k}\left\lVert\theta_{[k],\cdot}\right\rVert_{1},$

such that training will only require to scale each layer of $\theta$ and to
crossvalidate the multiplicative constant $C$. Since for $k\geq 1$,

$\displaystyle
C_{k}(\bar{\delta})=\sqrt{k}\times\sqrt{v_{\varepsilon}(\bar{\delta}/k+\log(pN)/k+\log
d)}\leq\sqrt{k}\times\sqrt{v_{\varepsilon}(\bar{\delta}+\log(pN)+\log d)},$

we let $\lambda_{k}=\frac{\sqrt{k}}{k!}$.

We now show that the minimization problem with layer-specific penalty can be
written as a standard regression problem with $\ell_{1}$ penalization by
rescaling the feature matrix, that is, multiplying
$\mathbf{S}_{N}^{\mathcal{D}}$ by a well-chosen diagonal matrix. Consider the
$\ell_{1}$-penalized problem

$\displaystyle\min\limits_{\theta\in\mathbb{R}^{s_{d}(N)\times
p}}\frac{1}{2M}\left\lVert\mathbf{Y}-\mathbf{S}_{N}^{\mathcal{D}}\theta\right\rVert_{\textnormal{F}}^{2}+C\sum\limits_{k=0}^{N}\lambda_{k}\left\lVert\theta_{[k],\cdot}\right\rVert_{1},$

where $C>0$ controls the strength of the penalization.

Making the change of variable

$\tilde{\theta}=\text{diag}(1,\underbrace{\lambda_{1},\dots,\lambda_{1}}_{d\text{
repetitions}},\underbrace{\lambda_{2},\dots,\lambda_{2}}_{d^{2}\text{
repetitions}},\dots,\underbrace{\lambda_{k},\dots,\lambda_{k}}_{d^{N}\text{
repetitions}})\theta,$

which is equivalent to

$\theta=\text{diag}(1,1/\lambda_{1},\dots,1/\lambda_{1},\dots,1/\lambda_{k},\dots,1/\lambda_{k})\tilde{\theta},$

and denoting by $W$ this last weight matrix, we get the equivalent
minimization problem

$\displaystyle\min\limits_{\tilde{\theta}\in\mathbb{R}^{s_{d}(N)\times
p}}\frac{1}{2M}\left\lVert\mathbf{Y}-\mathbf{S}_{N}^{\mathcal{D}}W\tilde{\theta}\right\rVert_{\textnormal{F}}^{2}+C\sum\limits_{k=0}^{N}\left\lVert\tilde{\theta}_{[k],\cdot}\right\rVert_{1}.$

We can thus obtain the SigLasso estimator by (i) multiplying the feature
matrix $\mathbf{S}_{N}^{\mathcal{D}}$ by $W$ and solving the associated
$\ell_{1}$-penalized problem (ii) multiplying the obtained solution by $W$.

The Learn-And-Reconstruct algorithm is the generic algorithm used in our work.
It is applicable for a wide variety of tasks such as missing values inference,
trajectory reconstruction, forecasting and many more. It is described in
Algorithm 1.

Algorithm 1 Learn-and-Reconstruct Algorithm. The algorithm infers for every
individual in the test set a reconstructed time series $\hat{Y}^{i}_{t}$.

1\. Learn the dynamics

Input: train dataset of normalized paths
$(\mathbf{X}^{1},\mathbf{Y}^{1}),\dots,(\mathbf{X}^{n},\mathbf{Y}^{n})$
sampled on $(D^{1},\bar{D}^{1}),\dots,(D^{n},\bar{D}^{n})$.

Construct the feature matrix $\mathbf{S}_{N}^{\mathcal{D}}$ and the target
vector $\mathbf{Y}$

for $i=1$ to $n$ do

for $t$ in $\bar{D}^{i}$ do

$\mathbf{S}_{N}^{\mathcal{D}}$ $\leftarrow$ Append $S_{N}(X^{i}_{[0,t]})$

$\mathbf{Y}$ $\leftarrow$ Append $Y^{i}_{t}$

end for

end for

Compute  $\hat{\theta}_{N,M}$ by solving (6) with
$\mathbf{Y},\mathbf{S}_{N}^{\mathcal{D}}$ using coordinate descent.

2\. Reconstruct trajectories

Input: test dataset $\tilde{\mathbf{X}}^{1},\dots,\tilde{\mathbf{X}}^{n}$
sampled on $\tilde{D}^{1},\dots,\tilde{D}^{n}$.

for $i=1$ to $n$ do

for $t$ in $\tilde{D}^{i}$ do

$\hat{Y}^{i}_{t}=\hat{\theta}_{N,M}S_{N}(\tilde{X}^{i}_{[0,t]})$

end for

end for

###  C.2 Assessing feature importance

We two metrics used to asses the importance of the different dimensions of the
feature path.

Given a truncation depth $N$ and a dimension $i\in\left\\{1,\dots,d\right\\}$,
we define its pure feature importance (PFI) as the sum of the norm of the
coefficients (or vectors in the case the target is multivariate) of
$\widehat{\theta}_{N,M}$ that are associated to signatures taken on the words
$I_{1}=(i)$, $I_{2}=(i,i)$, and so forth until $I_{N}=(i,\dots,i)$.
Mathematically,

$PFI(i)=\frac{1}{N}\Big{(}\left\lVert\theta^{I_{1}}\right\rVert_{2}+\left\lVert\theta^{I_{2}}\right\rVert_{2}+\dots+\left\lVert\theta^{I_{N}}\right\rVert_{2}\Big{)}.$

Since signatures also capture interactions between dimensions of the feature
path, we also define the cross feature importance (CFI) as the sum of norms of
the coefficients (or vectors) of $\widehat{\theta}_{N,M}$ that are associated
to signatures coefficients of words of length $\leq N$ in which the letter $i$
appears. Mathematically,

$CFI(i)=\frac{1}{s_{d}(N)-s_{d-1}(N)}\sum\limits_{I\text{ s.t.
}i\,\in\,I}\left\lVert\theta^{I}\right\rVert_{2}.$

For a given truncation depth $N$, note that there are
$s_{d}(N)-s_{d-1}(N)=\sum\limits_{k=0}^{N}d^{k}-\sum\limits_{k=0}^{N}(d-1)^{k}$
terms in the last sum, which justifies our choice of normalization.

###  C.3 Details on model implementation and evaluation

#### SigLasso.

The SigLasso model is implemented using the CVLasso class in scikit-learn
(Pedregosa et al., 2011). This implementation optimises the objective function
using coordinate descent and features automatic cross-validation of the
penalty strength. We use iisignature (Reizenstein & Graham, 2020) to compute
the signature of the feature time series. Every time series is standardized
prior to this through division by its own total variation, as suggested by
Morrill et al. (2020a). The depth of the signature is a hyperparameter chosen
between $2$ and $9$ or $6$ depending on the experiment. An intercept is added.

#### GRU.

The GRU is of width $128$ and systematically trained with $100$ epoches using
a learning rate of $0.001$.

#### Neural CDE.

We use the implementation of Neural CDE provided by torchcde (Kidger et al.,
2020). We use the original vector field described in the documentation of this
package, with the small tweak that we use a smoother non-linearity (tanh
instead of ReLU). We observed that using the rk4 solver instead of dopri5
significantly accelerates the training time of the Neural CDE without
affecting the model’s performances. The learning rate is hand-tuned to either
$0.001$ or $0.0001$ depending on the experiment. We train the model for $100$
epochs and asses its convergence by using a standard stopping criteria.

Metrics. The MSE is computed in a classical fashion. To compute the integrate
MSE, we compute the $L^{2}$ distance between the piecewise constant
interpolations of the true $y_{t}$ and the predicted $\widehat{y}_{t}$.

###  C.4 Details on the well specified model

#### Generation of the training data.

We generate a two-dimensional feature path by interpolating for every
dimension $15$ points in $[0,1]$, each of them being draw randomly for a
normal distribution $\mathcal{N}(0,1)$. The interpolation is done with Hermite
cubic splines with backward differences using the package torchcde (Kidger et
al., 2020). Time is added as a supplementary channel, which is a standard
practice when learning with signatures and Neural CDEs. These paths are then
downsampled by randomly drawing sampling points for the target and the feature
time series specific to every individual. The target path is the solution of a
CDE of the form

$dy_{t}=\sigma\Big{(}Ay_{t}\Big{)}dx_{t}$

where $\sigma:\mathbb{R}^{p}\to\mathbb{R}^{p}$ is the hyperbolic tangent,
$A\in\mathbb{R}^{d\times p}$ is a matrix drawn randomly from
$\mathcal{N}(0,I_{d\times p})$ and $x_{t}$ is the feature path constructed as
above. The solution of this CDE is computed using torchcde.

#### Generation of the test data.

We generate the test data in the same way than the training data. However,
this data is not downsampled as we wish to assess the generalization
capacities of our model—i.e., is our model capable of approximating the
dynamics and extrapolating to continuous feature paths.

###  C.5 Details on Ornstein-Uhlenbeck experiment

We take $(x_{t})_{t\in[0,1]}$ to be a 1-dimensional Brownian motion with
variance $\sigma^{2}=0.1$, and generate $(y_{t})$ as a $1$-dimensional
Ornstein-Uhlenbeck process driven by $(x_{t})$, that is, for all $t\in[0,1]$

$dy_{t}=\theta(\mu-y_{t})dt+dx_{t}.$

Simulation of $(y_{t})$ is done using a standard Euler-Maruyama simulation
scheme. We let $\theta=3$ and $\mu=1$. The training data is then downsampled
as in the well specified experiment.

###  C.6 Details on the tumor growth experiment

We consider the following tumor growth model taken from (Simeoni et al.,
2004). Let $x\in C^{1-var}([0,1],\mathbb{R})$. The weight $y\in
C([0,1],\mathbb{R}^{+})$ under the concentration of a treatment drug $x$ is
governed by the differential system

$\displaystyle
du^{1}_{t}=\Big{[}\Big{(}\lambda_{0}u_{t}^{1}\big{[}1+(\frac{\lambda_{0}}{\lambda_{1}}y_{t})^{\psi}\big{]}\Big{)}^{-1/\psi}-k_{2}x_{t}u_{t}^{1}\Big{]}dt$
$\displaystyle du^{2}_{t}=\Big{[}k_{2}x_{t}u_{t}^{1}-k_{1}u_{t}^{2}\Big{]}dt$
$\displaystyle du^{3}_{t}=\Big{[}k_{1}(u_{t}^{2}-u_{t}^{3})\Big{]}dt$
$\displaystyle du^{4}_{t}=\Big{[}k_{1}(u_{t}^{3}-u_{t}^{4})\Big{]}dt$
$\displaystyle y_{t}=u_{t}^{1}+u_{t}^{2}+u_{t}^{3}+u_{t}^{4}$

with initial condition
$(u_{0}^{1},u_{0}^{2},u_{0}^{3},u_{0}^{4},y_{0})=(2,0,0,0,2)$ and parameters
$(k_{1},k_{2},\lambda_{0},\lambda_{1},\psi)=(10,0.5,0.9,0.7,20)$. The
concentration $(x_{t})$ is chosen to be the squared value of the paths used
for the well-specified experiment. Notice that this system is non-linear
w.r.t. $x$. Indeed, writing $dy_{t}=\mathbf{G}(y_{t},x_{t})dt$, one has, for
$\alpha\in\mathbb{R}$, $G(y_{t},\alpha x_{t})\neq\alpha G(y_{t},x_{t})$. The
training data is then downsampled as in the well specified experiment.

###  C.7 Details on the French Covid experiment

We illustrate the performance of our method and competitors on French Covid
data from 2021-03-31 to 2021-07-07 available on Gitlab. Hospital data was
obtained from the SI-VIC database, the national inpatient surveillance system.

#### Target path.

Following Paireau et al. (2022), we chose the predict the growth rate of
incident hospitalisations in each of the 9 metropolitan regions of France. The
exponential growth rate was computed from raw data using a 2 days rolling
window and then smoothed using local polynomial regression as in Paireau et
al. (2022). Mathematically, our target time series is the $\mathbb{R}$-valued
growth rate, and we fit a different model for every of the 12 regions. It is
displayed for all 12 regions in Figure 6.

#### Feature path.

As in Paireau et al. (2022), we consider a set of 12 time-dependant predictors
of different types summarized in Table 2 and plotted in Figure 5. Both SIDEP
(”Système d’Information de Dépistage Populationnel”) and VAC-SI datasets are
publicly available. The mobility data was obtained from Google. The mobility-
related predictors describe travel trends for different kind of public spaces
such as such as shops and leisure spaces, food stores and pharmacies, parks,
public transport stations, workplaces and residential areas. The
meteorological data was obtained from Météo France.

Predictor | Type | Source | Description
---|---|---|---
prop_pos_symp | Epidemiological | SIDEP databas | proportion of positive tests among symptomatics
P_r | Epidemiological | SIDEP database | growth rate of positive tests
couv-complet | Epidemiological | VAC-SI database | proportion of vaccinated
grocery_and_pharmacy | Mobility | Google | visits to grocery and pharmacy stores
parks | Mobility | Google | visits to parks
transit_stations | Mobility | Google | visits visits to transit stations
workplaces | Mobility | Google | visits to workplaces
residential | Mobility | Google | visits to residential places
IPTCC | Meteorological | Météo France | Index PREDICT of climatic transmissivity
temperature | Meteorological | Météo France | temperature
rel_humidity | Meteorological | Météo France | relative humidity
abs_humidity | Meteorological | Météo France | absolute humidity
Table 2: The set of time-dependant predictors used to predict the hospital
admission growth rate. See Figure 5 for a vizualisation of all different
features across regions through time. Figure 5: The 12 different feature time
series used to forecast the hospitalization growth rate. Every different color
corresponds to a given region of France. Figure 6: Hospitalization growth rate
through time during the full period for the 12 different regions of France.

#### Models.

SigLasso (weighted and unweighted), NCDE and GRU algorithms were trained on
the period from 2021-03-31 to 2021-06-23 and tested on the period from
2021-06-24 to 2021-07-07. We included a history of 10 days at each point and
performed prediction for different horizons ranging from $1$ to $14$. In
others words, at horizon $h$, features values from day $t-h-10$ to day $t-h$
to were used to compute the prediction at time $t$. All feature time series
are normalized to have total variation equal to $1$.

#### Architectural details.

The GRU has width $128$ and is trained for $100$ epochs with a learning rate
of $0.0001$. The NCDE is trained for $30$ epochs with a learning rate of
$0.001$. It has $2$ hidden layers of width $128$, an intermediate
$\textnormal{Tanh}(\cdot)$ non-linearity and a final linear readout. This
architecture is identical to the one proposed in (Kidger, 2022). Penalty
strenght of the SigLasso is crossvalided using the internal implementation
LassoCV of scikit-learn (Pedregosa et al., 2011).

All details, in particular the features used for each individual prediction,
can be found in Paireau et al. (2022).

We refer to the supplementary information file of Paireau et al. (2022) and
our code for more details.

#### Results.

Figure 12 displays the RMSE (on all regions) of NCDE, SigLasso, GRU, and the
Ensemble method for all prediction horizons $h=1,\ldots,14$. Figures 7, 8, 10
and 9 display the obtained interpolation for weighted and unweighted SigLasso,
GRU and NCDE at different horizons (corresponding to different line colors in
winter matplotlib palette). The lighter the blue, the smaller the time
horizon: the lightest curve corresponds to a time horizon equal to $h=1$.
Truth is in red.

Figure 7: Interpolation (left of dotted line) and prediction (right of dotted
line) of hospitalization growth rate for all 12 french regions using weighted
SigLasso. Figure 8: Interpolation (left of dotted line) and prediction (right
of dotted line) of hospitalization growth rate for all 12 french regions using
unweighted SigLasso. Figure 9: Interpolation (left of dotted line) and
prediction (right of dotted line) of hospitalization growth rate for all 12
french regions using NCDE. Figure 10: Interpolation (left of dotted line) and
prediction (right of dotted line) of hospitalization growth rate for all 12
french regions using GRU. Figure 11: Interpolation (left of dotted line) and
prediction (right of dotted line) of hospitalization growth rate for all 12
french regions using ensemble methods (Paireau et al., 2022). Figure 12: RMSE
accross all regions on the training period (left) and the testing period
(right) for the ensemble methode (Paireau et al., 2022), NCDE, GRU, and
SigLasso (weighted and unweighted). See Figure 13 for a zoom-in on GRU and
SigLasso performances. Figure 13: RMSE accross all regions on the training
period (left) and the testing period (right) for GRU and SigLasso (weighted
and unweighted).

###  C.8 Additional results

We give in Table 3 some additional results on the experiments described above.

Table 3: Training time of SigLasso, GRU and Neural CDE in different simulation settings, averaged over 10 iterations. In every setting, $n=50$, $\\#\bar{D}^{i}=5$ for all $i=1,\dots,n$ (and therefore $M=250$). | Training time (s)
---|---
| SigLasso | GRU | Neural CDE
Well-specified | 0.37 $\pm$ 0.23 | 269 $\pm$ 109 | 1754 $\pm$ 587
OU | 0.057 $\pm$ 0.005 | 27 $\pm$ 0.44 | 216 $\pm$ 2.7
Tumor growth | 0.056 $\pm$ 0.007 | 31 $\pm$ 3.5 | 250 $\pm$ 14