# Renormalization group approach to connect discrete- and continuous-time descriptions of Gaussian processes

Federica Ferretti, Victor Chardès, Thierry Mora, Aleksandra Walczak, Irene Giardina

► To cite this version:

Federica Ferretti, Victor Chardès, Thierry Mora, Aleksandra Walczak, Irene Giardina. Renormalization group approach to connect discrete- and continuous-time descriptions of Gaussian processes. Physical Review E, 2022, 105 (4), pp.044133. 10.1103/PhysRevE.105.044133. hal-03854599

HAL Id: hal-03854599

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

Submitted on 23 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Renormalization group approach to connect discrete- and continuous-time
descriptions of Gaussian processes

Federica Ferretti<sup>1,2</sup> Victor Chardès,<sup>3</sup> Thierry Mora,<sup>3</sup> Aleksandra M. Walczak,<sup>3</sup> and Irene Giardina<sup>1,2,4</sup>

<sup>1</sup>Dipartimento di Fisica, Università Sapienza, 00185 Rome, Italy

<sup>2</sup>Istituto Sistemi Complessi, Consiglio Nazionale delle Ricerche, UOS Sapienza, 00185 Rome, Italy

<sup>3</sup>Laboratoire de Physique de l'École Normale supérieure (PSL University), CNRS,
Sorbonne Université, Université de Paris, 75005 Paris, France

<sup>4</sup>INFN, Unità di Roma 1, 00185 Rome, Italy

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

(Received 16 January 2021; revised 29 May 2021; accepted 17 February 2022; published 22 April 2022)

Discretization of continuous stochastic processes is needed to numerically simulate them or to infer models from experimental time series. However, depending on the nature of the process, the same discretization scheme may perform very differently for the two tasks, if it is not accurate enough. Exact discretizations, which work equally well at any scale, are characterized by the property of invariance under coarse-graining. Motivated by this observation, we build an explicit renormalization group (RG) approach for Gaussian time series generated by autoregressive models. We show that the RG fixed points correspond to discretizations of linear SDEs, and only come in the form of first order Markov processes or non-Markovian ones. This fact provides an alternative explanation of why standard delay-vector embedding procedures fail in reconstructing partially observed noise-driven systems. We also suggest a possible effective Markovian discretization for the inference of partially observed underdamped equilibrium processes based on the exploitation of the Einstein relation.

DOI: 10.1103/PhysRevE.105.044133

I. INTRODUCTION

How to properly represent continuous-time stochastic pro-
cesses by discrete-time descriptions is a central problem in
applied science. Models used in all areas of physics frequently
take the form of stochastic differential equations (SDEs), yet
the experimental observation of any real process produces
discrete sequences of data points. Moreover, numerical inte-
gration of any SDE requires to define transition probabilities
over finite time steps. Hence, finding a good discretization
scheme is key for both numerical integration and parametric
inference of continuous-time processes.

Since exact solutions are not available for arbitrary pro-
cesses, the general strategy is to resort to a Taylor-Itô
expansion of the integrated process. The lowest order ex-
pansion is the Euler-Maruyama scheme, which is (strongly)
convergent as τ¹/², being τ the discretization time step [1].
This discretization scheme has been widely used in the lit-
erature, due to its simplicity and intuitive interpretation (it
amounts to estimating velocities as finite differences). How-
ever, it has also been observed that the Euler-Maruyama
approximation cannot be employed in the derivation of
parametric inference methods for second or higher-order pro-
cesses, as it leads to the extraction of inconsistent parameter

estimators [2–10]. On the contrary, adopting higher-order
discretization schemes, consistent inference algorithms (both
Bayesian and non-Bayesian) can be designed.

The failure of the Euler-Maruyama discretization (Euler
for short) can be explained in the following terms: to learn
second or higher-order stochastic models, the dynamical in-
formation we need to extract from the fluctuations is of
higher order (in τ) than the accuracy with which those are
reconstructed by the Euler scheme. Nonetheless, the correla-
tion functions that we reconstruct from an Euler simulation
faithfully reproduce those of the integrated continuous-time
process, at arbitrary order—provided that the simulation
time step (τ<sub>sim</sub>) is sufficiently small, compared to the typical
time scales of the system (τ<sub>obs</sub>). While this separation of scales
is possible in numerical simulations, state space inference
requires to discretize the process over the same time scales as
those over which we observe the correlations. The lack of such
separation of scales explains the bad performance of the Euler
discretization in parametric inference tasks, as opposed to
numerical integration, and suggests that as we “zoom out” our
lens for the observation of the process, the dynamics recovers
the correct statistical properties that were originally missing.

This property has motivated the introduction of data
augmentation techniques in parametric inference, like
Refs. [11,12], which consist in introducing and marginaliz-
ing additional intermediate states between pairs of observed
points. It also reminds us of what happens in statistical field
theory in the context of the renormalization group (RG). The
bare theory, which provides an effective description of the
system at the microscopic scale, might be incomplete; how-
ever, when we apply the RG procedure to get an effective

<footnote>Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.</footnote>
theory at larger scales, this generates the missing terms and yields a model that correctly describes the large scale statistical fluctuations. More precisely, the augmentation strategy of Refs. [11,12] or the “zoom out” operation described above amount to temporal coarse-graining, which is reminiscent of the coarse-graining implemented in Real Space RG, with time steps playing the role of lattice spacings in the usual setting.

In this paper, we identify the formal framework to develop this analogy with the RG approach, and we exploit it to recognize different classes of discrete-time stochastic processes based on their property of invariance under the RG transformation. The condition of invariance under RG amounts to requiring that, when the discretization interval $\tau$ is small enough, integrating the process over a single step $2\tau$ or through a combination of two steps on intervals $\tau$ should provide the same result.

We restrict our analysis to stationary Gaussian processes, for which an explicit RG map can be easily derived. In the space of these discrete-time models, Euler discretizations of linear SDEs are a subset. The study of the RG map reveals an interesting structure, where the fixed points point out classes of “natural” physical processes. They include full observations of continuous Markov processes, or partial observations of higher-dimensional processes inheriting a non-Markovian structure. No intermediate situation can be obtained. Higher-order Markov processes, whose transition probabilities depend on two or more previous observations, are not RG fixed points. This result underscores the lack of finite-dimensional delay vector embeddings for stochastic systems.

The paper is organized in the following way: in Sec. II we consider as a starting point the linear damped Langevin equation, and illustrate in this simple case how an RG procedure can be defined and applied to discretizations of the continuous process. We derive the RG map, identify the fixed points, and give their interpretation. We also discuss the possibility of building effective Markovian embeddings and their limitations, an issue relevant for inference purposes. In Sec. III we extend the previous analysis to arbitrary higher-order linear processes, and deduce more general conclusions. Finally, in Sec. IV we summarize and discuss our results.

## II. LINEAR SECOND-ORDER PROCESSES

### A. A simple case

The shortcomings of the Euler scheme mentioned in the Introduction are manifested for processes of second or higher order. By nth order processes we mean that they are described by nth order SDEs or, equivalently, that they are obtained from the partial observation of n-dimensional stochastic processes, whose structure is such that the noise is transferred from the hidden coordinates to the single observed degree of freedom (hypoelliptic diffusions). Partial observations of such noise-driven systems break the Markovianity, introducing temporal noise correlations and memory effects in the description of the observed dynamics [13,14]. These features are not captured by the Euler discretization scheme, which has the peculiar property of reducing to a Markovian discrete-time model when applied to this kind of systems.

The simplest example, for $n=2$, is a linear damped Langevin equation in which only the positional coordinates are directly measured. Let us consider

$$dx = vdt, \quad (1)$$

$$dv = -\eta vdt - \kappa xdt + \sigma dW, \quad (2)$$

with $W(t)$ a Wiener process. We assume that the coordinate $x$ is observed at a finite sampling rate $\tau^{-1}$, producing an infinite time series denoted as $\{X_n, n \in \mathbb{N}\}$. For this process, under stationary conditions, an exact solution can be computed and used in inference and simulation problems [15,16]; however, we are interested here in understanding how discrete models are related to their continuous counterparts. We therefore discretize Eqs. (1) and (2), i.e., we integrate the continuous equations over a time interval of length $\tau$ and expand at the first order in $\tau$ the resulting integral expressions (Euler). Eliminating the $v$ variable, we obtain a discrete update equation for the $x$ coordinate that has the following structure:

$$X_n = \psi X_{n-1} + \theta X_{n-2} + \mu \epsilon_n, \quad (3)$$

where $\psi = 2 - \eta\tau - \kappa\tau^2$, $\theta = (-1 + \eta\tau)$, $\mu = \sigma\tau^{3/2}$ and $\epsilon_n \sim N(0, 1)$.

Equation (3) is an autoregressive model of order two, denoted AR(2) [17], and it is fully characterized by the conditional probability $P(X_n|X_{n-1}, X_{n-2})$. Thanks to the Markovian structure of the discrete process in Eq. (3), the probability of the associated time series reads

$$P(\{X_n, n \in \mathbb{N}\}) = \prod_{n \ge 2} P(X_n|X_{n-1}, X_{n-2})P(X_0, X_1). \quad (4)$$

Under stationary assumptions, we can move the initial condition arbitrarily far in the past to neglect boundary terms. The resulting $P(\{X_n\})$ can be interpreted as the Boltzmann weight of a configuration of spins on an infinite one-dimensional lattice with first- and second-nearest-neighbor interactions. The analogy will be helpful to derive an explicit RG map.

### B. The RG construction

The failure of the Euler discretization Eq. (3) in inference approaches signals that it misses some important information about fluctuations on the local scale $\tau$. The questions we address are the following: is there a way, starting from Eq. (3), to understand what ingredients are missing? Are there specific constraints on the coefficients of discrete equations like Eq. (3), to ensure compliance with any continuous model?

We start our analysis by noticing that a good discretization of the original continuous process must be so regardless of the precise value of the time interval $\tau$ we consider (as long as it is small enough, compared to the physical time scales of the process). If we use the same scheme to integrate over a single step $2\tau$ or through a combination of two steps on intervals $\tau$, then the result should be the same. As we iterate this argument many times, we compel the structure of the discrete equations to remain unaltered at different scales. We look indeed at the discrete process on scales that are larger and larger than the original one, until the ratio between the characteristic observation time scale and the discretization time step becomes infinite ($\tau_{obs}/\tau_{sim} \to \infty$). In this way, upon redefinition of the time units, we approach the continuum limit. This
<figure><img src="image_3.png" /><figcaption>FIG. 1. Sketch of the RG transformation for a time series. The discretely observed coordinates {$X_n$} play the role of spin variables on a lattice, with nearest neighbor and next-nearest neighbor coupling. The time step $\tau$ is the lattice spacing. Model parameters, $A^0$, are transformed into $\tilde{A}^0$ through coarse-graining, then rescaled to $A^1$ as we map $\tau \to \tau/2$.</figcaption></figure>

procedure allows us to check whether a given discrete scheme is scale invariant, and therefore faithfully describes a reference continuous equation, or—in case invariance is violated—how the approach to the continuum limit occurs.

Let us formalize this idea using the language of the Renormalization Group. We explore the continuum limit through a progressive increase in the number of steps contained in a fixed time window. This is obtained by iterating the two operations that make up the RG: (i) coarse graining and (ii) joint rescaling of the time unit and of the parameters of the model. The procedure is sketched in Fig. 1.

We apply the method to Eq. (3). Given the analogy of Euler time series with linear spin chains, we adopt the strategy of decimation to coarse grain [18–20]: the goal is to get rid of half of the sites (e.g., odd ones) in the sequence generated by Eq. (3), thus deriving effective update equations for the even subseries. To implement this transformation, we take a suitable linear combination of neighboring update equations of the form of Eq. (3):

$$ \mathrm{Eqn}(X_n) + \psi \mathrm{Eqn}(X_{n-1}) - \theta \mathrm{Eqn}(X_{n-2}), \quad (5) $$

which results into an update equation for the subseries:

$$ X_n = \tilde{\psi}X_{n-2} + \tilde{\theta}X_{n-4} + \tilde{r}_n, \quad (6) $$

with $\tilde{\psi} = \psi^2 + 2\theta$, $\tilde{\theta} = -\psi^2$. The structure of Eq. (6) looks similar to the one of Eq. (3) with updated parameters. There is, however, a crucial difference: unlike in the original process, the random increment

$$ \tilde{r}_n = \mu[\epsilon_n + \psi\epsilon_{n-1} - \theta\epsilon_{n-2}] \quad (7) $$

is now correlated across nearest neighbors: $\mathbb{E}[\tilde{r}_n\tilde{r}_{n\pm2}] \neq 0$. This fact is better seen if we appropriately rewrite the random increment as a different linear combination of Gaussian variables:

$$ \tilde{r}_n = \tilde{\mu}\tilde{\epsilon}_n + \tilde{\nu}\tilde{\epsilon}_{n-2}, \quad (8) $$

with $\tilde{\epsilon}_i \sim \mathcal{N}(0, 1)$ new I.I.D. variables, and $\tilde{\mu}, \tilde{\nu}$ satisfying

$$ \mathbb{E}[\tilde{r}_n^2] = \tilde{\mu}^2 + \tilde{\nu}^2 = (1 + \psi^2 + \theta^2)\mu^2, \quad (9) $$

$$ \mathbb{E}[\tilde{r}_n \tilde{r}_{n\pm2}] = \tilde{\mu} \tilde{\nu} = -\theta \mu^2. \quad (10) $$

<figure><img src="image_4.png" /><figcaption>FIG. 2. Coarse-graining by decimation maps the ladder topology, corresponding to our AR(2) model, to a fully connected net. The new couplings are however not independent due to the nontrivial correlation properties of the random contributions.</figcaption></figure>

It can be verified that $\mathbb{E}[\tilde{r}_n^l \tilde{r}_{n\pm2k}^l] = 0$ for $k > 1$. Thus, Eq. (6) becomes

$$ X_n = \tilde{\psi}X_{n-2} + \tilde{\theta}X_{n-4} + \tilde{\nu}\tilde{\epsilon}_{n-2} + \tilde{\mu}\tilde{\epsilon}_n, \quad (11) $$

also known as ARMA(2,1) model [17]. General autoregressive moving-average processes of order $(p, q)$, denoted ARMA($p, q$), are time series generated by update equations of the form

$$ X_n = \sum_{i=1}^{p} \phi_i X_{n-i} + \sum_{i=1}^{q} \nu_i \epsilon_{n-i} + \mu \epsilon_n, \quad (12) $$

with $\nu_i, \mu \in \mathbb{R}$ and $\epsilon_n \sim \mathcal{N}(0, 1)$. I.I.D. The autoregressive (AR) part of the equation corresponds to the contribution from the previous $p$ states of the system; the moving average (MA) part, of order $q$, corresponds to the second sum in the right-hand side of Eq. (12), and is responsible for the time correlation of random increments.

It is evident that, after coarse-graining, the process in Eq. (6) is no longer Markovian, as its probability distribution cannot be factorized into the product of single-step transition probabilities. Recalling the analogy with the spin chain, RG creates, since the first iteration, infinite-range effective couplings, even though we started with just first- and second-nearest-neighbor bonds. This is indeed the effect of simple decimation on the zigzag ladder topology, to which the AR(2) process of Eq. (3) corresponds (Fig. 2). However, these emerging couplings are not independent, and the four parameters in Eq. (11) are sufficient to characterize them.

To get a closed-form RG transformation, we need to ensure that further iterations of the decimation procedure do not keep introducing novel higher-order terms. Luckily, the ARMA(2,1) structure is stable, as an example of a more general result that we discuss in the following (Sec. III). If we call $\mathcal{A} = (\psi, \theta, \mu, \nu)$ the set of parameters of a generic ARMA(2,1) model, and $\tilde{\mathcal{A}} = (\tilde{\psi}, \tilde{\theta}, \tilde{\mu}, \tilde{\nu})$ the parameters of the coarse grained model obtained after decimation, then applying RG yields a well-defined map from $\mathcal{A}$ to $\tilde{\mathcal{A}}$.

The second operation that completes the RG iteration is rescaling the time step, $2\tau \to \tau$, and reabsorbing this change of units through a redefinition of the parameters. The parameters of our models are dimensionless, yet their dependency on $\tau$ is what determines how to connect any discrete process to its continuous-time counterpart. We express $\psi, \theta$ as asymptotic
<table><caption>TABLE I. Fixed-point solutions of the RG recurrence relations up to third order in $\tau$. We find four manifolds of fixed points, corresponding to four types of processes, parametrized by the arbitrary constants $u, s, z$, and $b$. In addition to the reported ones, there are diverging fixed points.</caption><thead><tr><th rowspan="2">Model</th><th colspan="8">AR coefficients</th><th colspan="8">MA coefficients</th></tr><tr><th>ψ<sub>0</sub><sup>*</sup></th><th>θ<sub>0</sub><sup>*</sup></th><th>ψ<sub>1</sub><sup>*</sup></th><th>θ<sub>1</sub><sup>*</sup></th><th>ψ<sub>2</sub><sup>*</sup></th><th>θ<sub>2</sub><sup>*</sup></th><th>ψ<sub>3</sub><sup>*</sup></th><th>θ<sub>3</sub><sup>*</sup></th><th>α<sub>0</sub><sup>*</sup></th><th>β<sub>0</sub><sup>*</sup></th><th>α<sub>1</sub><sup>*</sup></th><th>β<sub>1</sub><sup>*</sup></th><th>α<sub>2</sub><sup>*</sup></th><th>β<sub>2</sub><sup>*</sup></th><th>α<sub>3</sub><sup>*</sup></th><th>β<sub>3</sub><sup>*</sup></th></tr></thead><tbody><tr><td>A MA(0)</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>s</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>B AR(1)</td><td>1</td><td>0</td><td>u</td><td>0</td><td>u<sup>2</sup>/2</td><td>0</td><td>u<sup>3</sup>/6</td><td>0</td><td>0</td><td>0</td><td>s</td><td>0</td><td>us</td><td>0</td><td>2u<sup>2</sup>s/3</td><td>0</td></tr><tr><td>C ARMA(2,1)</td><td>-1</td><td>-1</td><td>u</td><td>2u</td><td>-u<sup>2</sup>/2</td><td>-2u<sup>2</sup></td><td>u<sup>3</sup>/6</td><td>(2u)<sup>3</sup>/6</td><td>0</td><td>0</td><td>4s</td><td>s</td><td>-8us</td><td>-2us</td><td>32u<sup>2</sup>s/3</td><td>13u<sup>2</sup>s/6</td></tr><tr><td>D ARMA(2,1)</td><td>2</td><td>-1</td><td>u</td><td>-u</td><td>z</td><td>-u<sup>2</sup>/2</td><td>u(6z - u<sup>2</sup>)/12</td><td>-u<sup>3</sup>/6</td><td>0</td><td>0</td><td>-2s</td><td>s</td><td>-2us</td><td>us</td><td>4b - (2z + 3u<sup>2</sup>)s</td><td>b</td></tr></tbody></table>

power series of $\tau$, $\psi(\tau) = \sum_k \psi_k \tau^k$ and $\theta(\tau) = \sum_k \theta_k \tau^k$, and work, up to the desired order, with recursive relations for the coefficients of the series expansion, $\psi_k$ and $\theta_k$. These coefficients are now dimensional and get rescaled with the time unit. The same idea can be applied to $\nu$ and $\mu$ by expanding them in powers of $\tau^{1/2}$. It is then convenient to reparametrize the noise amplitudes as

$$ \alpha := \mathbb{E}r_n^2 = \mu^2 + \nu^2; \quad \beta := \mathbb{E}r_n r_{n+1} = \mu\nu, \quad (13) $$

since their asymptotic series expansion involves integer powers of $\tau$: $\alpha(\tau) = \sum_k \alpha_k \tau^k$, $\beta(\tau) = \sum_k \beta_k \tau^k$.

The physical dimension of each coefficient $A_k \in \{\psi_k, \theta_k, \alpha_k, \beta_k\}$ is now set by the order of the corresponding term in the series expansion. Each of them gets rescaled, after coarse graining, as $A_k^{l+1} = 2^{-k} \tilde{A}_k^l$, where $l$ is an index counting the RG iterations. This operation defines the RG map as a set of recursive equations in an infinite-dimensional parameter space, made up of the Taylor coefficients parametrizing the four functions $\psi^l(\tau)$, $\theta^l(\tau)$, $\alpha^l(\tau)$, $\beta^l(\tau)$:

$$ \psi_k^{l+1} = 2^{-k} \left[ 2\theta_k^l + \sum_{i=0}^{k} \psi_i^l \psi_{k-i}^l \right], \quad (14) $$

$$ \theta_k^{l+1} = -2^{-k} \sum_{i=0}^{k} \theta_i^l \theta_{k-i}^l. \quad (15) $$

Similar equations (reported in Appendix A) for the coefficients $\alpha_k$ and $\beta_k$ are obtained from the following definition of the coarse-grained parameters $\tilde{\alpha}^l(\tau)$ and $\tilde{\beta}^l(\tau)$:

$$ \tilde{\alpha}^l = [1 + (\psi^l)^2 + (\theta^l)^2]\alpha^l + 2\psi^l(1 - \theta^l)\beta^l, \quad (16) $$

$$ \tilde{\beta}^l = \psi^l(1 - \theta^l)\beta^l - \theta^l\alpha^l. \quad (17) $$

Notice that since the recursion equations at order $k$ only involve lower orders, they can be solved recursively over $k$, and can also be truncated to an arbitrary order while retaining a closed form.

### C. Fixed points

Our interest is in finding the fixed points of the recursion relations defined above. RG fixed points capture indeed the scale-invariant behavior of a system. Table I reports the solutions for the parameters $A^*$ of the ARMA(2,1) model, showing their coefficients $A_k^*$ up to order $k=3$. There are four manifolds of fixed points in this space, corresponding to different classes of stochastic processes. Their projection onto

the plane of leading order AR coefficients, $\psi_0$ and $\theta_0$, is shown in Fig. 3.

The first class of fixed points A corresponds to sequences of independent random variables ($\psi = \theta = \nu = 0$). Fixed points B are AR(1) processes, i.e., processes where only two time points are involved and there is no memory in the noise term: they can be interpreted as discretizations of linear first-order SDEs of the form

$$ dx = uxdt + \sqrt{sdW}. \quad (18) $$

By induction on the fixed-point equations for the $A_k$ coefficients ($k \in \mathbb{N}$), it is possible to prove that the fixed-point parameters of the B process are: $\theta^*(\tau) = \beta^*(\tau) = 0$, $\psi^*(\tau) = e^{u\tau}$, $\alpha^*(\tau) = \frac{s}{2}(e^{2u\tau} - 1)$.

In addition, there are two other fixed-point manifolds corresponding to ARMA(2,1) models, denoted by C and D. C is not a continuous process, but evolves through finite jumps. The AR coefficients, which appear in the deterministic part of the equation, read $\psi(\tau) = -e^{-u\tau}$ and $\theta(\tau) = -e^{-2u\tau}$. The process has a three-branched phase diagram, and evolves in time jumping from one branch to the next one—approaching the origin or moving away from it depending on the sign of $u$. This structure is invariant under the RG transformation we defined, but it is precisely determined by the details of the

<figure><img src="image_5.png" /><figcaption>FIG. 3. Fixed-point manifolds of the RG map, projected on the plane of leading autoregressive coefficients $\psi_0, \theta_0$. The interior of the triangle is the basin of attraction of fixed point A. The equal sides are the basin of attraction of point B. The basins of points C and D are contained in the basis of the triangle (plus vertex on the top for D). The Euler AR(2) process has coordinates $(2, -1)$ in this plane. Color-shaded trajectories represent the solution of recurrence relations from three sample initial conditions. Arrows show the direction of (discrete) moves and should not be interpreted as continuous flow lines.</figcaption></figure>
decimation procedure. If the coarse graining was implemented differently, for instance by trimming two points out of three, then the same type of process would not be a fixed point.

Finally, the D class represents the discretization of a partially observed two-dimensional process of the form:

$$dx = vdt + \sigma_x dW_x, \quad (19)$$

$$dv = -\eta v dt - \kappa x dt + \sigma dW_v, \quad (20)$$

of which Eq. (1) is a particular case [21]. The $(\psi, \theta)$ coefficients of fixed point D reconstruct a second derivative at leading order, and the effect of linear drift at first order (with $u = -\eta$ in Table I). The variable $s = -\sigma_x^2$ in Table I encodes noise added to the $x$ variable, and is zero in inertial models like Eq. (1).

Inertial models are of special interest. In this case, since $s=0$, noise contributions are determined by third-order coefficients. If we start with the Euler model Eq. (3) as an initial condition for the RG recurrence relations, whose noise coefficients are $\alpha_k^0 = \sigma^2 \delta_{k,3}$, $\beta_k^0 = 0$, then the associated RG flow, on the plane of MA parameters, reads

$$\alpha_3^l = \sigma^2 \left(\frac{2}{3} + \frac{1}{3}4^{-l}\right), \quad \beta_3^l = \sigma^2 \left(\frac{1}{6} - \frac{1}{6}4^{-l}\right). \quad (21)$$

Upon renormalization, the Euler time series Eq. (3) falls into class D, with $b = \sigma^2/6$ and $s = 0$. While the initial model is strongly convergent as $\tau^{1/2}$, the resulting fixed point is at least convergent as $\tau^{3/2}$. The asymptotic values $(\alpha_3^{\infty}, \beta_3^{\infty})$ are those we would obtain if we applied a higher-order discretization scheme like [22] to Eq. (1) in the first place. A more exhaustive analysis of fixed point D and of its parameters can be found in Appendix B.

To summarize the analysis performed so far, we have shown that the Euler discretization of the inertial equation we considered is not invariant upon RG transformations. This encapsulates the fact that Euler does not appropriately describe the statistical properties of the continuous model on the original discretization scale. We also showed that, as RG progresses, Euler leads to effective discrete models that are, on the contrary, invariant. This convergence occurs exponentially fast—see Eq. (21)—and explains why numerical simulations performed with the Euler scheme work. Finally, we note that higher-order discretizations of Eqs. (1) and (2) already possess an ARMA(2,1) structure and are therefore invariant upon RG transformations. This explains why they provide good algorithms for local (in time) inference approaches, in contrast to Euler.

### D. Effective Markov description

Although we discussed that Euler-based maximum likelihood inference approaches for stochastic dynamical systems are inconsistent, numerical and analytical evidence [2–5,23] hints at a possible effective Markov discretization of the form of Eq. (3) for generalized Langevin equations at equilibrium. Such discretization neglects noise correlations, but employs a rescaled damping coefficient $\eta' = (2/3)\eta$. In other words, the information missed by the Euler scheme at the discretization scale $\tau$ can be effectively reabsorbed into a single rescaling factor.

For linear processes, this behavior can be interpreted as following from the Einstein relation, which imposes a well-

known relationship between the drag coefficient and the noise amplitude (through temperature and mass) [24,25]—see Appendix D. In our case, this relation involves the parameters $\psi_1$ and $\alpha$. If we select the fixed-point value of the diagonal entry of the noise covariance matrix, $\alpha_3^* = (2/3)\sigma^2$, instead of the Euler value $\alpha_3^0 = \sigma^2$, and neglect off-diagonal elements, then we find that we must also rescale $\eta$ by $2/3$ by imposing the Einstein relation. In Appendix D we show how this effective Euler-like discretization can be derived in the simple case of a linear Langevin equation like Eq. (1). Of course, this effective Markov description can only be used for maximum likelihood inference at the scale of $\tau$. The effective discretization gives, by construction, consistent maximum likelihood estimators, but cannot be iterated in simulations to reproduce the process on longer scales. Because of the convergence of the Euler-Maruyama integration algorithm, the simulated process would be the underdamped Langevin model with a damping coefficient equal to $\eta' = (2/3)\eta$.

In the context of continuous processes, the idea of building an effective Markov model that reproduces the statistics of a non-Markov linear model (only at the stationary level) has been recently investigated in Ref. [26].

### III. GENERALIZATION TO ARMA(p, q) PROCESSES

The RG procedure we detailed for ARMA(2,1) models can be generalized to arbitrary ARMA($p, q$) processes, defined by Eq. (12). Decimation of the time series can be done by combining $p$ neighboring equations of the form of Eq. (12), in a similar way as in Eq. (5) (see Appendix C). This decimation step yields a new update equation with an ARMA($p, \tilde{q}$) structure, where

$$\tilde{q} = \left\lfloor \frac{p+q}{2} \right\rfloor, \quad (22)$$

and $\lfloor \cdot \rfloor$ denotes the rounding down operation. The condition of invariance under RG imposes that fixed points satisfy $\tilde{q} = q$, implying $q = p$ or $q = p - 1$.

This result leads to two important observations. First, it shows that purely autoregressive models of order $p \ge 2$, AR($p$)=ARMA($p, 0$) (i.e., models with no memory in the noise), cannot be stable points and thus cannot be exact discretizations of stochastic differential equations of second or higher order. As a consequence, for partially observed processes ($p \ge 2$) a Markov description of the dynamics ($q = 0$) is impossible. Second, the sharp selection of $q$ reveals that longer memory than $p$ is irrelevant in the RG sense.

Our finding is related to the nonexistence of exact delay-vector embeddings for noise-driven systems. Embedding approaches consist of stacking a finite number of subsequent points to define a new dynamical variable $\mathcal{X} = (X_{n-p+1}, \dots, X_n)$—known as delay vector [27,28]—and assuming that it follows a Markov dynamics. An equivalent approach is to estimate the derivatives of the observed stochastic process through differences of subsequent measurements: the Euler discretization is a simple example of this procedure. In both cases, the Markov dynamics of the embedded process is described by AR($p$) models. In conclusion, although the delay-vector method is standard in deterministic contexts, it cannot be directly extended to
stochastic processes, as partial observation sets strong limitations to phase space reconstruction for stochastic dynamical systems [29–33].

## IV. CONCLUSION

In this work we investigated the connection between discrete and continuous descriptions of stochastic Gaussian processes. Using an RG approach, we studied how the properties of given discretization schemes change as we change the observation scale of the process, from the original discretization step, up to the continuum limit $\tau/\tau_{obs} \to 0$. We started our analysis focusing on the simple case of a linear damped Langevin equation, and then generalized the procedure to higher-order processes. In this broader framework, we defined an RG map on the space of ARMA(p, q) models, i.e., a class of generative models for Gaussian processes, and identified exact discretizations of continuous-time stochastic processes through its fixed points.

The issue relates to the embedding problem and is especially relevant for the design of parametric inference methods. A possible suggestion that can be drawn from our analysis is to abandon the Markov setting in favor of descriptions with correlated noise, by introducing an additional noise delay for any new coordinate in the delay vector. The RG construction shows indeed that time correlations in the noise terms are spontaneously generated to match the original dimension of the partially observed dynamical system. Consistently, inference algorithms that take into account the right noise correlations since the beginning correctly estimate the parameters of the underlying continuous model, without the need of augmentation techniques [2–6].

Nonetheless, we also showed that, up to a parameter rescaling, effective discretizations based on second-order Markov models—AR(2)—can still be used for the inference of damped equilibrium processes. This observation raises the question of whether there exist alternative (e.g., variational) RG transformations exist that can map AR(2) models into other AR(2) processes, allowing for an effective memoryless description of Gaussian processes.

How to extend these results to nonlinear processes also remains an open question.

## ACKNOWLEDGMENTS

We thank A. Vulpiani, M. Baldovin, A. Cavagna for helpful conversations, and A. C. Costa for useful comments and discussions. This work was partially supported by the ERC Consolidator Grant No. 724208, ERC Advanced Grant No. 785932 and by the Italian Ministry of Foreign Affairs and International Cooperation through the Adinmat project.

## APPENDIX A

In this Appendix we look in greater details at the fixed points of the recursive equations. Our starting point are Eqs. (14)–(17). Since the recursive equations for $\psi, \theta$ are independent of those for $\alpha$ and $\beta$, it is convenient to start with the coefficients of the autoregressive (AR) part of the model, and later focus on the stochastic (MA) contribution.

## 1. AR coefficients

From Eqs. (14) and (15), at each order $k$ of the series expansion, we obtain a two-dimensional map from which the fixed points for the considered coefficients, $(\psi_k^*, \theta_k^*)$, can be extracted. Notice that higher-order recurrence relations are only coupled to lower order ones, so one can solve them iteratively. Here we are just interested in the study of the fixed points.

The recursive relations at the leading order ($k=0$) read

$$ \psi_0^{l+1} = (\psi_0^l)^2 + 2\theta_0^l, \quad (A1) $$

$$ \theta_0^{l+1} = -(\theta_0^l)^2. \quad (A2) $$

As it is evident from Eqs. (14) and (15), only at this order the system is nonlinear; at any subsequent order the map is linear. There are four fixed points for the nonlinear recurrence relations Eqs. (A1) and (A2):

A. $(\psi_0^*, \theta_0^*) = (0, 0)$. This is a process with null autoregressive part. Moving average contributions are specified by the recurrence relations for the MA coefficients.

B. $(\psi_0^*, \theta_0^*) = (1, 0)$. This fixed point corresponds to a first order process of the form $X_n = X_{n-1} + r_n$, i.e., an ARMA(1,q) process; the moving average order q will be determined by the recurrence relations for $\alpha_k$ and $\beta_k$.

C. $(\psi_0^*, \theta_0^*) = (-1, -1)$. This is an ARMA(2,q) process, of the form $X_n = -X_{n-1} - X_{n-2} + r_n$. Here again, the random contribution will be specified by the recursion relations for $\alpha_k$ and $\beta_k$.

D. $(\psi_0^*, \theta_0^*) = (2, -1)$. This model is a specific case of an ARMA(2,q) model, known as ARIMA(1,1,q) model [17]. Thanks to the specific value assumed by the coefficients of the AR part, one can indeed rewrite the process as $(1-L)^2 X_n = r_n$, where L is the lag operator: $LX_n = X_{n-1}$. So $(1-L)$ is the discrete differencing operator and $\{(1-L)X_n\} = \{\bar{V}_n\}$ is the reconstructed velocity series. We leave once again the value of q unspecified for the moment, as it is determined by the analysis of the recurrence relations for the MA coefficients.

Notice that only the fixed point A is an asymptotically stable point, whose basin of attraction is the interior of the triangle in Fig. 3. The other points are unstable, at least in some directions.

We are especially interested in the class of discrete-time models represented by point D, since, at the leading order, the model can be considered as a discretized version of a second-order SDE of the kind $\ddot{x} = \xi$, with $\xi(t)$ a Gaussian noise (white or colored, depending on the MA coefficients). Higher-order contributions can modify the coefficients $\psi$ and $\theta$ in front of $X_{n-1}$ and $X_{n-2}$, but they will not affect the interpretation of such a process as a discretization of a second-order SDE. The only difference will be in the addition of position- or velocity-dependent linear forces.

Proceeding to the next order, $k=1$, the study of the fixed points reveals that their number and their nature (i.e., the kind of dynamical models they correspond to) is left unchanged. The AR coefficients of the A fixed point are still constrained to zero $(\psi_1 = \theta_1 = 0)$ while a new free parameter ($u$ in Table I) appears at order $k=1$ for the other fixed points. The new parameter accommodates the arbitrariness of what can be
interpreted as a linear force in the overdamped case B, or a damping force in case D.

Let us now examine what happens at order $k = 2$. Given the recurrence relations:

$$ \psi_2^{l+1} = \frac{1}{4}[2\psi_0^l\psi_2^l + (\psi_1^l)^2 + 2\theta_2^l], \quad (A3) $$

$$ \theta_2^{l+1} = -\frac{1}{4}[2\theta_0^l\theta_2^l + (\theta_1^l)^2], \quad (A4) $$

the fixed points, at this order, are expanded as follows:

$$ \text{A. } (\psi_0^*, \theta_0^*, \psi_1^*, \theta_1^*, \psi_2^*, \theta_2^*) = (0, 0, 0, 0, 0, 0). $$

$$ \text{B. } (\psi_0^*, \theta_0^*, \psi_1^*, \theta_1^*, \psi_2^*, \theta_2^*) = (1, 0, u, 0, u^2/2, 0). $$

$$ \text{C. } (\psi_0^*, \theta_0^*, \psi_1^*, \theta_1^*, \psi_2^*, \theta_2^*) = (-1, -1, u, 2u, -u^2/2, -2u^2). $$

$$ \text{D. } (\psi_0^*, \theta_0^*, \psi_1^*, \theta_1^*, \psi_2^*, \theta_2^*) = (2, -1, u, -u, z, -u^2/2). $$

Compared to the previous fixed points, we have here an additional arbitrariness on $\psi_2^*$. The novel free parameter $z$ accounts for an $x$-dependent force in the second-order continuous-time SDE, of which the fixed-point model D can be considered a discretization. Indeed the discretization

picture above keeps holding, even if the ARIMA(1,1,q) structure is lost in favor of an ARMA(2,q) one, when $\psi_2^* \neq u^2/2$. For points B and C, on the contrary, no new free parameter appears and the values of $\psi_2$ and $\theta_2$ are fixed by the values of lower order coefficients. In these cases, the dimension of the fixed-point manifold does not change by moving from $k = 1$ to $k = 2$.

Of course the study of the fixed points can be developed to any desired order, but the discussion we made so far is already sufficient to characterize their autoregressive nature. Indeed, it can be proven by induction that the AR order of the fixed points is left unchanged at subsequent orders in the expansion.

## 2. MA coefficients

Recurrence relations for the MA coefficients $\mu$ and $\nu$ or, equivalently, for the parameters $\alpha$ and $\beta$, depend on those of the AR coefficients as shown in Eqs. (16) and (17). The great advantage of working with $\alpha$ and $\beta$ is that, when the AR coefficients $\psi$ and $\theta$ are fixed, their recurrence relations are linear:

$$ \tilde{\alpha}_k^l = (1 + \psi_0^2 + \theta_0^2)\alpha_k^l + 2\beta_k^l \psi_0(1 - \theta_0) + \sum_{i=0}^{k-1} \left[ \alpha_i^l \sum_{j=0}^{k-i} (\psi_j \psi_{k-i-j} + \theta_j \theta_{k-i-j}) + 2\beta_i^l \psi_{k-i}(1 - \theta_0) - 2\beta_i^l \sum_{j=0}^{k-i-1} \psi_j \theta_{k-i-j} \right]. \quad (A5) $$

$$ \tilde{\beta}_k^l = \beta_k^l \psi_0(1 - \theta_0) - \alpha_k^l \theta_0 - \sum_{i=0}^{k-1} \alpha_i^l \theta_{k-i} + \sum_{i=0}^{k-1} \beta_i^l \left[ \psi_{k-i}(1 - \theta_0) - \sum_{j=0}^{k-i-1} \psi_j \theta_{k-i-j} \right]. \quad (A6) $$

We analyze the behavior of these parameters for each of the four fixed-point classes we found above.

At the leading order we have a homogeneous system:

$$ \alpha_0^{l+1} = (1 + \psi_0^2 + \theta_0^2)\alpha_0^l + 2\psi_0(1 - \theta_0)\beta_0^l, \quad (A7) $$

$$ \beta_0^{l+1} = \psi_0(1 - \theta_0)\beta_0^l - \theta_0\alpha_0^l. \quad (A8) $$

The single fixed point of such system is the origin $(\alpha_0^*, \beta_0^*) = (0, 0)$, unless the coefficients $\theta_0$ and $\phi_0$ take values that render the fixed-point equations linearly dependent. This happens only for class A, where the condition $(\psi_0^*, \theta_0^*) = (0, 0)$ leaves $\alpha_0^*$ as a free parameter, while $\beta_0^* = 0$. As a result, $q = 0$, and this fixed-point model is just a sequence of I.I.D. Gaussian variables. For autoregressive models (i.e., fixed points B, C, and D), the condition $(\alpha_0^*, \beta_0^*) = (0, 0)$ tells us that $\tau$-independent noise contributions are prohibited. This feature enables us to interpret these models as discretizations of SDEs [34].

Moving to a first-order expansion (the lowest nontrivial one for the last three fixed points), the recursion relations Eqs. (A5) and (A6) take the form of a 2D affine map. For the four fixed points we find

$$ \text{A. } \text{MA}(q): (\alpha_0^*, \beta_0^*, \alpha_1^*, \beta_1^*) = (s, 0, 0, 0), s \in \mathbb{R}. $$

$$ \text{B. } \text{ARMA}(1,q): (\alpha_0^*, \beta_0^*, \alpha_1^*, \beta_1^*) = (0, 0, s, 0), \text{ with } s \text{ a real parameter. As expected, even at this order} $$

$q = 0$ and the process reduces to a simple AR(1) model.

C. ARMA(2,q): $(\alpha_0^*, \beta_0^*, \alpha_1^*, \beta_1^*) = (0, 0, 4s, s)$, $s \in \mathbb{R}$. We have a manifold of fixed points, represented by a line on the plane of first order covariance coefficients.

D. AR(I)MA(2,q): $(\alpha_0^*, \beta_0^*, \alpha_1^*, \beta_1^*) = (0, 0, -2s, s)$. Again, we have a line of fixed-point solutions. The corresponding process can be interpreted as the first order discretization of a partially observed SDE. More details can be found in Appendix B.

At second order, the structure of the MA part of the fixed points remains almost the same: A and B remain memoryless processes ($q = 0$), and no new parameters appear for C and D. Therefore we report directly the result for $k = 3$, the last order where arbitrariness can be introduced to modify the structure of the fixed points:

A. MA(0): $(\alpha_0^*, \beta_0^*, \alpha_1^*, \beta_1^*, \alpha_2^*, \beta_2^*, \alpha_3^*, \beta_3^*) = (s, 0, 0, 0, 0, 0, 0, 0)$.

B. ARMA(1,0): $(\alpha_0^*, \beta_0^*, \alpha_1^*, \beta_1^*, \alpha_2^*, \beta_2^*, \alpha_3^*, \beta_3^*) = (0, 0, s, 0, us, 0, \frac{2}{3}\psi_1^2 s, 0)$.

C. ARMA(2,1): $(\alpha_0^*, \beta_0^*, \alpha_1^*, \beta_1^*, \alpha_2^*, \beta_2^*, \alpha_3^*, \beta_3^*) = (0, 0, 4s, s, -8us, -2us, \frac{32}{3}u^2s, \frac{13}{6}u^2s)$.

D. AR(I)MA(2,1): $(\alpha_0^*, \beta_0^*, \alpha_1^*, \beta_1^*, \alpha_2^*, \beta_2^*, \alpha_3^*, \beta_3^*) = (0, 0, -2s, s, -2us, us, 4b - (2z + 3u^2)s, b)$. For this nontrivial fixed point, the recurrence relations at this order are linearly dependent, admit infinitely many fixed-point solutions. We parametrize them by taking $\beta_3^* = b$.
APPENDIX B

This section is devoted to a more detailed discussion of the fourth fixed point, which, due to its physical meaning, we think deserves a special focus. We have already highlighted that, at leading order, the AR coefficients reproduce a second time derivative through the double differencing operator $\Delta^2 = (1 - L)^2$, with L the lag operator: $LX_n = X_{n-1}$. This fact gives to the model an “integrated process” nature at leading order, also known as ARIMA(1,1,1) [17]. Since the relation

$$ \lim_{\tau \to 0} \frac{\psi - \psi_0^*}{\tau} = - \lim_{\tau \to 0} \frac{\theta - \theta_0^*}{\tau} \quad (B1) $$

also holds, this integrated process structure is kept up to $O(\tau)$. Deviations from it occur at higher order and are due to the presence of linearly x-dependent forces.

The fixed-point model is indeed the consistent discretization of a general class of partially observed continuous-time processes of the form:

$$ dy = Aydt + BdW, \quad (B2) $$

where $y = (x, v)^T$, $W = (W_x, W_v)^T$, with $W_x(t)$ and $W_v(t)$ independent Wiener processes and

$$ A = \begin{pmatrix} -\lambda & 1 \\ -\kappa & -\eta \end{pmatrix}, \quad BB^T = \begin{pmatrix} \sigma_{xx}^2 & \sigma_{xv}^2 \\ \sigma_{vx}^2 & \sigma_{vv}^2 \end{pmatrix}. \quad (B3) $$

Notice that setting $A_{12} = 1$ does not imply a loss of generality: compared to $A_{12} = a \in \mathbb{R} \setminus \{0\}$, it just corresponds to a rescaling of the time unit, which does not alter the process (the only caveat is that $a < 0$ would revert the time direction). The

case $A_{12} = 0$ is not of interest for us, as it would decouple the dynamics of the unobserved degrees of freedom from that of the observed ones. The entries in the drift matrix A must satisfy the stability condition, i.e., assuming that time evolves in the positive direction, ($t \ge 0$), it must be negative semidefinite. Finally, we have $\sigma_{xv}^2 = \sigma_{vx}^2$ for the symmetry of the diffusion matrix.

Partial observation of the process in Eq. (B2) yields a Gaussian process described, at the continuous level, by a generalized Langevin equation (GLE), and, at the discrete level, by an ARMA(2,1) model. The GLE can be obtained by integrating the continuous equation for the v variable and plugging it into the first one in Eq. (B2):

$$ dx(t) = \left\{ e^{-\eta t} v(0) - \kappa \int_{0}^{t} dt' e^{-\eta(t-t')} x(t') + \int_{0}^{t} [B_{21} dW_{x}(t') + B_{22} dW_{v}(t')] dt - \lambda x(t) + B_{11} dW_{x}(t) + B_{12} dW_{v}(t) \right\} dt \quad (B4) $$

This equation contains a parametric dependence on the initial condition $v(0)$, which can be removed during discretization as in Refs. [3,22]:

$$ X_n = \frac{(e^{2A\tau})_{12}}{(e^{A\tau})_{12}} X_{n-1} - \left[ (e^{2A\tau})_{11} - \frac{(e^{2A\tau})_{12}}{(e^{A\tau})_{12}} (e^{A\tau})_{11} \right] X_{n-2} \\ = r_n, \quad (B5) $$

with

$$ r_n = \int_{t_n-2\tau}^{t_n} (e^{A(2\tau-s)})_{11} [B_{11} dW_x(s) + B_{12} dW_v(s)] - \frac{(e^{2A\tau})_{12}}{(e^{A\tau})_{12}} \int_{t_n-2\tau}^{t_n-\tau} (e^{A(\tau-s)})_{11} [B_{11} dW_x(s) + B_{12} dW_v(s)] \\ + \int_{t_n-2\tau}^{t_n} (e^{A(2\tau-s)})_{12} [B_{21} dW_x(s) + B_{22} dW_v(s)] - \frac{(e^{2A\tau})_{12}}{(e^{A\tau})_{12}} \int_{t_n-2\tau}^{t_n-\tau} (e^{A(\tau-s)})_{12} [B_{21} dW_x(s) + B_{22} dW_v(s)]. \quad (B6) $$

Equation (B5) corresponds to an ARMA(2,1) process, since Eq. (B6) implies that $E[r_n r_m] = \alpha \delta_{n,m} + \beta \delta_{n,m \pm 1}$.

It is possible to explicitly work out the calculation to find the mapping from (A, BB$^T$) to $(\psi, \theta, \alpha, \beta)$—and hence find the relation between Eq. (B3) and the variables parametrizing the fixed point D in Table I. Let us start by performing a small $\tau$ expansion for the AR coefficients:

$$ \psi = \frac{(e^{2A\tau})_{12}}{(e^{A\tau})_{12}} \approx 2 - (\eta + \lambda)\tau + \frac{1}{2}\tau^2(-2\kappa + \eta^2 + \lambda^2) + \frac{1}{6}\tau^3(3\kappa(\eta + \lambda) - \eta^3 - \lambda^3), \quad (B7) $$

$$ \theta = (e^{2A\tau})_{11} - \frac{(e^{2A\tau})_{12}}{(e^{A\tau})_{12}} (e^{A\tau})_{11} = -e^{-(\eta+\lambda)\tau} \approx -1 + (\eta + \lambda)\tau - \frac{1}{2}\tau^2(\eta + \lambda)^2 + \frac{1}{6}\tau^3(\eta + \lambda)^3. \quad (B8) $$

Identifying the coefficients $\psi_0, \dots, \psi_3$ and $\theta_0, \dots, \theta_3$, we deduce that the model parameters in Eq. (B3) and the fixed-point parameters of process D in Table I are linked by the following relation:

$$ u = -(\lambda + \eta); \quad z = -\kappa + \frac{1}{2}(\eta^2 + \lambda^2). \quad (B9) $$

The algebra for the derivation of $\alpha$ and $\beta$ is more laborious but not complicated. The results are

$$ \alpha \approx 2\sigma_{xx}^2\tau - 2\sigma_{xx}^2(\eta + \lambda)\tau^2 + \frac{2}{3}[\sigma_{vv}^2 + 2\eta\sigma_{xv}^2 + \sigma_{xx}^2(-\kappa + 3\eta^2 + 3\eta\lambda + 2\lambda^2)]\tau^3, \quad (B10) $$

$$ \beta \approx -\sigma_{xx}^2\tau + (\eta + \lambda)\sigma_{xx}^2\tau^2 + \frac{1}{6}[\sigma_{vv}^2 + 2\eta\sigma_{xv}^2 + \sigma_{xx}^2(2\kappa - 3\eta^2 - 6\eta\lambda - 4\lambda^2)]\tau^3. \quad (B11) $$

Hence, we can deduce

$$ s = -\sigma_{xx}^2, \quad (B12) $$

$$ b = \frac{1}{6}[\sigma_{vv}^2 + 2\eta\sigma_{xv}^2 + \sigma_{xx}^2(2\kappa - 3\eta^2 - 6\eta\lambda - 4\lambda^2)]. \quad (B13) $$
Further manipulation of the sums in Eq. (C3) allows us to find formal expressions for the AR coefficients of the coarse grained process, $\tilde{\phi}_{i=1...p}$:

$$ 
\begin{aligned}
\tilde{\phi}_1 &= 2\phi_2 + \phi_1^2, \\
\tilde{\phi}_i &= 2\phi_{2i} + (-1)^{i+1}\phi_i^2 + 2\phi_{2i-1} \sum_{k=1}^{i-1} (-1)^{k+1}\phi_k \quad \text{for} \quad 2 \le i \le \lfloor \frac{p}{2} \rfloor - 1, \\
\tilde{\phi}_i &= (-1)^{i+1}\phi_i^2 \quad \text{for} \quad \lfloor \frac{p}{2} \rfloor + 1 \le i \le p.
\end{aligned}
\quad (C7) $$

The coarse-grained coefficient $\tilde{\phi}_{\lfloor p/2 \rfloor}$ takes a different form depending on $p$ being even or odd:

$$ p \text{ odd}: \quad \tilde{\phi}_{\lfloor p/2 \rfloor} = 2\phi_{2\lfloor p/2 \rfloor} + (-1)^{1+\lfloor p/2 \rfloor} \phi_{\lfloor p/2 \rfloor}^2, \quad (C8) $$

$$ p \text{ even}: \quad \tilde{\phi}_{\lfloor p/2 \rfloor} = 2\phi_{2\lfloor p/2 \rfloor} + (-1)^{1+\lfloor p/2 \rfloor} \phi_{\lfloor p/2 \rfloor}^2 + 2\phi_{2\lfloor p/2 \rfloor-1} \sum_{k=1}^{\lfloor \frac{p-1}{2} \rfloor} (-1)^{k+1} \phi_k. \quad (C9) $$

Let us now restart from Eq. (C3) and focus on the random term $r_n$. Since linear combinations of Gaussian variables are still Gaussian, one can properly redefine the $\epsilon_m$'s and rearrange the coefficients in front of them to rewrite $\tilde{r}_n = \sum_{i=0}^{\tilde{q}} \tilde{\mu}_i \epsilon_{n-2i}$, where

$$ \tilde{q} = \lfloor \frac{p+q}{2} \rfloor. \quad (C10) $$

We deduce there are only two invariant scenarios for ARMA($p, q$) processes: $q = p$ or $q = p-1$. This fact tells us that partial (discrete) observation of continuous-time processes let memory emerge: each hidden degree of freedom increases by one the order of both the AR part and the MA part of the discrete model, thus introducing color.

## APPENDIX D

In this Appendix we motivate the exploitation of the seemingly universal 2/3 rescaling factor in effective Markov models discussed in the main text, and provide a physical interpretation for it, working out the reference problem of an integrated Ornstein-Uhlenbeck (OU) process.

The integrated OU process is the simplest example of second-order SDE for which the Euler-related inconsistency appears. It is described by Eq. (1) with $\kappa = 0$, where we suppose we can only observe (with infinite accuracy) the inertial degree of freedom, at a sampling rate $\tau^{-1}$. Let us recall the notation for the time series of empirical observations $\{X_n\}_{n \in \mathbb{N}}$, and for the time series of reconstructed velocities $\{\tilde{V}_n\}_{n \in \mathbb{N}}$, where

$$ \tilde{V}_n = \frac{(X_{n+1} - X_n)}{\tau}. \quad (D1) $$

Let us also introduce the series of real velocities $\{V_n\}_{n \in \mathbb{N}}$, corresponding to the one we would obtain if we were able to measure directly the velocity degree of freedom. Because the evolution of the $v$ variable is described by an independent first-order SDE when $\kappa = 0$, the time series $\{V_n\}_{n \in \mathbb{N}}$ is described by an AR(1) process. On the contrary, the evolution of the $x$ variables is non Markovian and expressed, at the continuous level, via a generalized Langevin equation. Consequently, the time series $\{\tilde{V}_n\}_{n \in \mathbb{N}}$ inherits a nonzero MA order, ending in an ARMA(1,1).

Nonetheless, we may ask whether it is possible to describe it with an effective AR(1) process, which would correspond

to an effective AR(2) process for the $\{X_n\}_{n \in \mathbb{N}}$ series. Let us write a putative AR(1) model for the series of reconstructed velocities:

$$ \tilde{V}_n - (1 - \alpha)\tilde{V}_{n-1} = \sigma\epsilon_n, \quad (D2) $$

where $\epsilon_n \sim \mathcal{N}(0, 1)$ I.I.D. and $\alpha, \sigma$ are parameters to fix. The goal is to find an effective memoryless discrete model for $\{\tilde{V}_n\}_{n \in \mathbb{N}}$ that reproduces correctly the sufficient statistics used by Bayesian and non-Bayesian parametric inference approaches. The common characteristic of these approaches is that of being derived from a Taylor-Itô expansion in the small $\tau$ limit. They just exploit local dynamical information to learn the laws governing the evolution of the system, typically carried by the first few elements of the autocovariance of the time series.

In the case of AR(1) models, a sufficient statistics corresponds to the set $S_1 = \{\mathbb{E}[\tilde{V}_n^2], \mathbb{E}[\tilde{V}_n \tilde{V}_{n+1}]\}$, i.e., the autocorrelation functions at equal times and at a distance of one time step [35]. We impose on the observables $S_1$ the two following consistency conditions:

i. The stationary distribution of $\tilde{V}_n$ is the Maxwell Boltzmann distribution at temperature $T$ ($k_B = 1$):

$$ \mathbb{E}[\tilde{V}_n^2] = T. \quad (D3) $$

ii. The relation between the reconstructed acceleration $\bar{A}_n = (\tilde{V}_{n+1} - \tilde{V}_n)/\tau$ and the reconstructed velocity $\tilde{V}_n$ is the one we can exactly compute for the integrated Ornstein-Uhlenbeck process:

$$ \mathbb{E}[\bar{A}_n | \tilde{V}_n] \underset{\tau \to 0}{\sim} -\frac{2}{3} \eta \tilde{V}_n. \quad (D4) $$

A detailed derivation is in Ref. [2]. From this relation we immediately derive the condition to impose on the observables of interest:

$$ \mathbb{E}[\tilde{V}_{n+1} \tilde{V}_n] = (1 - \frac{2}{3} \eta \tau) \mathbb{E}[\tilde{V}_n^2]. \quad (D5) $$

The self-correlation function $\mathbb{E}[\tilde{V}_n \tilde{V}_{n+k}]$ of an AR(1) process of the form of Eq. (D2) is explicitly known:

$$ \mathbb{E}[\tilde{V}_n \tilde{V}_{n+k}] = \frac{(1 - \alpha)^{|k|} \sigma^2}{1 - (1 - \alpha)^2}. \quad (D6) $$

Taking its value at $k = 1$ and using Eq. (D3) in Eq. (D5) yields the expected result $\alpha = (2/3)\eta\tau$. Computing the self-correlation of the reconstructed velocities at $k = 0$ and
imposing equipartition from Eq. (D3), we set the value of σ :

$$
\sigma^2 = T[1 - (1 - \alpha)^2] \underset{\tau \to 0}{\sim} 2T\alpha. \qquad (D7)
$$

This is the celebrated Einstein relation.

In conclusion, it is possible to describe the sequence of
measurements of an integrated OU process as an effective

AR(2) series with a rescaled damping coefficient η' = (2/3)η.
This result is intuitive and could have been grossly derived by
imposing the Einstein relation [which comes from Eq. (D3)
alone], and selecting only the main diagonal of the covariance
matrix of random increments to cancel color, thus implying
σ² = 2/3(2Tητ).

[1] P. E. K. E. Platen and P. E. Kloeden, *Numerical Solution of Stochastic Differential Equations*, 1st ed., Vol. 23 of Stochastic Modelling and Applied Probability (Springer-Verlag, Berlin, 1992).

[2] J. N. Pedersen, L. Li, C. Grădinaru, R. H. Austin, E. C. Cox, and H. Flyvbjerg, Phys. Rev. E **94**, 062401 (2016).

[3] F. Ferretti, V. Chardès, T. Mora, A. M. Walczak, and I. Giardina, Phys. Rev. X **10**, 031018 (2020).

[4] B. Lehle and J. Peinke, Phys. Rev. E **97**, 012113 (2018).

[5] B. Lehle and J. Peinke, Phys. Rev. E **91**, 062113 (2015).

[6] A. Gloter, Scand. J. Stat. **33**, 83 (2006).

[7] A. Samson and M. Thieullen, Stoch. Processes Appl. **122**, 2521 (2012).

[8] S. Ditlevsen and A. Samson, J. Roy. Stat. Soc.: Ser. B (Stat. Methodol.) **81**, 361 (2019).

[9] Q. Clairon and A. Samson, Stat. Inference Stoch. Processes **23**, 105 (2020).

[10] Y. Pokern, A. M. Stuart, and P. Wiberg, J. Roy. Stat. Soc.: Ser. B (Stat. Methodol.) **71**, 49 (2009).

[11] O. Elerian, S. Chib, and N. Shephard, Econometrica **69**, 959 (2001).

[12] B. Eraker, J. Bus. Econ. Stat. **19**, 177 (2001).

[13] R. Zwanzig, *Nonequilibrium Statistical Mechanics* (Oxford University Press, New York, NY, 2001).

[14] M. S. Miguel and J. M. Sancho, J. Stat. Phys. **22**, 605 (1980).

[15] R. E. Kalman, J. Basic Eng. **82**, 35 (1960).

[16] D. T. Gillespie, Phys. Rev. E **54**, 2084 (1996).

[17] P. Brockwell and R. Davis, *Introduction to Time Series and Forecasting* (Springer, Berlin, 2002).

[18] L. P. Kadanoff, Ann. Phys. **100**, 359 (1976).

[19] G. Parisi, *Statistical Field Theory, Frontiers in Physics* (Addison-Wesley, Redwood City, CA, 1988).

[20] E. Brézin, *Introduction to Statistical Field Theory* (Cambridge University Press, Cambridge, UK, 2010).

[21] In fact, the fixed point D includes a more general class of two-dimensional continuous-time processes of the form dy = Ay + BdW, as described in Appendix B. However, all these processes are projected, after marginalization of the velocities, to the same ARMA(2,1) model. What changes is just the mapping from the coefficients of the continuous-time process A and B to those of the ARMA(2,1) model (clearly noninvertible).

[22] R. D. Skeel and J. A. Izaguirre, Mol. Phys. **100**, 3885 (2002).

[23] A. Gloter, Statistics **35**, 225 (2001).

[24] R. Kubo, Rep. Prog. Phys. **29**, 255 (1966).

[25] U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Phys. Rep. **461**, 111 (2008).

[26] G. Gradziuk, G. Torregrosa, and C. P. Broedersz, Phys. Rev. E **105**, 024118 (2021).

[27] F. Takens, *Detecting Strange Attractors in Turbulence* (Springer-Verlag, Berlin, 1981), Vol. 898, pp. 366–381.

[28] M. Casdagli, S. Eubank, J. Farmer, and J. Gibson, Physica D **51**, 52 (1991).

[29] J. Stark, J. Nonlin. Sci. **9**, 255 (1999).

[30] J. Stark, D. S. Broomhead, M. E. Davies, and J. Huke, J. Nonlin. Sci. **13**, 519 (2003).

[31] L. Onsager and S. Machlup, Phys. Rev. **91**, 1505 (1953).

[32] F. Borra and M. Baldovin, Chaos **31**, 023102 (2021).

[33] A. C. Costa, T. Ahamed, D. Jordan, and G. Stephens, Maximally predictive ensemble dynamics from data, arXiv:2105.12811 (2021).

[34] C. Gardiner, *Stochastic Methods* (Springer-Verlag, Berlin, 2009), Vol. 13.

[35] J. Lee and S. Pressé, J. Chem. Phys. **137**, 074103 (2012).