# Multi-level stochastic refinement for complex time series and fields: A
Data-Driven Approach

M. Sinhuber1, J. Friedrich2,3, R. Grauer2 and M. Wilczek1,∗ 1 Max Planck
Institute for Dynamics and Self-Organization, Am Faßberg 17, D-37077
Göttingen, Germany  2 Theoretische Physik I, Ruhr-Universität Bochum,
Universitätsstr. 150, D-44780 Bochum, Germany  3 Univ. Lyon, ENS de Lyon,
Univ. Claude Bernard, CNRS, Laboratoire de Physique, F-69342, Lyon, France  ∗
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Spatio-temporally extended nonlinear systems often exhibit a remarkable
complexity in space and time. In many cases, extensive datasets of such
systems are difficult to obtain, yet needed for a range of applications. Here,
we present a method to generate synthetic time series or fields that reproduce
statistical multi-scale features of complex systems. The method is based on a
hierarchical refinement employing transition probability density functions
(PDFs) from one scale to another. We address the case in which such PDFs can
be obtained from experimental measurements or simulations and then used to
generate arbitrarily large synthetic datasets. The validity of our approach is
demonstrated at the example of an experimental dataset of high Reynolds number
turbulence.

##  1 Introduction

Many high-dimensional nonlinear systems display a remarkable degree of
complexity in space and time. As a consequence, spatio-temporal records of
such systems often exhibit non-trivial statistical features, including multi-
scale correlations and scale-dependent deviations from Gaussianity. In
principle, a comprehensive statistical characterization of such systems
requires obtention of the full multi-time–multi-point statistics – a
prohibitive task given the dimensionality of typical systems.

Significant simplifications arise if a hierarchy of scales in space and/or
time can be identified, which is statistically characterized by transition
probabilities from larger to smaller scales (or vice versa). If the statistics
on one scale only depend on the previous one, the system is Markovian in
scale. This property has been validated empirically for many complex systems,
including hydrodynamic turbulence [1, 2, 3], solar wind turbulence [4],
neuronal spike trains [5] and currency exchange markets [6, 7]. For scale-
Markovian systems, the transition probabilities contain the full statistical
information of the system, leading to a tremendous reduction of complexity.
Even if the Markov property is weakly violated or a priori unknown, a
Markovianization in scale can still provide a valuable approximation,
containing, for example, non-trivial multi-scale correlations.

In many situations, however, a purely statistical characterization is
insufficient. Rather, a complete spatio-temporal record, i.e. a realization of
a complex system, is needed. For example, machine learning applications
require very large training datasets. Providing such datasets from
experimental measurements, computer simulations or historical records is often
challenging, sometimes even impossible. Even if data is available, it
sometimes lacks resolution for a specific target application due to
experimental or computational constraints, such that refinement is necessary.

In this work, we introduce a data-driven method, which allows generating
statistically well-defined time series or fields with non-trivial features in
a hierarchical fashion. The key idea of this _multi-level stochastic
refinement_ (MLSR) is to generate a realization of a stochastic process
representing the multi-scale features of a generic complex system given the
knowledge of transition probability density functions (PDFs) on all relevant
scales. Our emphasis is on a refinement of coarse time series or fields in
scale, which complements previous approaches that focused on forward-in-time
construction of synthetic time series by utilizing Markov properties in scale
in conjunction with a Fokker-Planck-equation based approach [8, 9]. We here
focus on the case in which the transition PDFs can be sourced from
experimental measurements or simulations.

As one prototypical example of a complex system, we apply our approach to
turbulent flows, which are ubiquitous in nature and engineering applications.
Turbulent flows feature a broad range of dynamically active scales in space
and time. For example, the separation of the large, energy-containing length
scales $L$ and the small, dissipative length scales $\eta$ increases in
hydrodynamic turbulence with the Reynolds number like $L/\eta\sim Re^{3/4}$.
As a consequence, fully resolved simulations are computationally prohibitive
for realistic natural and engineering flows. Large-eddy simulations tackle
this problem by only resolving the energy-containing scales [10, 11]. However,
in many applications, resolving small-scale turbulence is key. Examples
include cloud microphysics [12, 13, 14], turbulent combustion [15, 16], mixing
and transport of atmospheric pollutants [17, 18], the assessment of turbulent
loads on wind turbines [19], as well as the development of active control
strategies for wind farms [20]. Recently, deep learning algorithms for feature
identification and extraction have been gaining attention to deal with these
challenges. For the required large amounts of training data [21, 22],
synthetic turbulent fields, like the one presented in the following, could
make a valuable contribution to alleviate the challenge of acquiring the
databases.

The paper is structured as follows: We start with a description of the
algorithm and the theoretical background in section 2. We then demonstrate the
data-driven algorithm at the example of state-of-the-art wind tunnel
measurements of high-Reynolds number turbulence in section 3. In the
conclusions 4 we discuss future applications and generalizations as well as an
outlook on model-driven approaches to the MLSR method.

##  2 Description of the data-driven MLSR based on three-point PDFs

Figure 1: Schematic representation of the data-driven multi-level stochastic
refinement at the example of synthetic time series (mapped to a spatial signal
using Taylor’s hypothesis [23]) of wind tunnel turbulence. The initial
boundary values are set on a scale of $L=32\lambda$ ($\lambda$ denotes the
Taylor scale in hydrodynamic turbulence). The velocity at $L=16\lambda$, i.e.
on half the scale, is then drawn from the the conditional PDF
$p(u_{L/2}|u_{0L},u_{L})$ which previously has been obtained from experimental
data. By iterating this procedure, the synthetic time series (gray) is
obtained, which shares statistical features with the original experimental
data.

Let us consider the situation in which we aim to generate a realization of a
statistically homogeneous, one-dimensional field of length $L$ with a
resolution of $\Delta=L/2^{N}$. Given initial field values on scale $r_{0}=L$,
we adopt a hierarchical refinement on scale $r_{1}=L/2$, then progressing to
scale $r_{2}=L/4$ etc. until scale $r_{N}=\Delta=L/2^{N}$ is reached after $N$
refinement steps.

On each refinement level, let us assume we know the three-point PDF
$f(u_{\mathrm{l}},u_{\mathrm{c}},u_{\mathrm{r}};r_{i})$ from experimental
measurements. Here $u_{\mathrm{l}}$, $u_{\mathrm{c}}$ and $u_{\mathrm{r}}$
denote the field values at the left end, in the center, and at the right end
of the interval, respectively, which are separated by the distance $r_{i}$.
Based on the three-point PDF, we can define the conditional PDF

$p_{i}\left(u_{\mathrm{c}}|u_{\mathrm{l}},u_{\mathrm{r}}\right)\equiv
p\left(u_{\mathrm{c}}|u_{\mathrm{l}},u_{\mathrm{r}},r=r_{i}\right)=\frac{f(u_{\mathrm{l}},u_{\mathrm{c}},u_{\mathrm{r}};r_{i})}{f(u_{\mathrm{l}},u_{\mathrm{r}};2r_{i})}\,.$
(1)

Here, $f(u_{\mathrm{l}},u_{\mathrm{r}};2r_{i})$ denotes the two-point PDF on
scale $2r_{i}$, which can be obtained from the three-point PDF by integrating
out the center point:

$f(u_{\mathrm{l}},u_{\mathrm{r}};2r_{i})=\int\mathrm{d}u_{\mathrm{c}}\,f(u_{\mathrm{l}},u_{\mathrm{c}},u_{\mathrm{r}};r_{i})\,.$
(2)

The conditional PDF (1) can be interpreted as a transition PDF in scale. It
describes the probability of drawing $u_{\mathrm{c}}$, which can be used to
obtain a field which is sampled on scale $r_{i}$, given the two values
$u_{\mathrm{l}}$ and $u_{\mathrm{r}}$ at the end points of the interval of
length $2r_{i}$. In our numerical implementation, we employ a standard
rejection method [24] to draw from this conditional PDF. We thereby construct
a realization sampled on the scale $r_{i}$. The three-point statistics on
scale $r_{i}$ are consistent with the three-point statistics of the underlying
experimental data by construction. Employing this multi-level stochastic
refinement algorithm hierarchically starting from $r_{1}$ and progressing down
to scale $r_{N}$ enforces consistency with the underlying data in each
refinement step.

We illustrate the procedure in figure 1 for a time series of flow velocities
in turbulence. Here we construct a synthetic turbulent velocity signal at a
resolution of $L/2^{9}=\lambda/16$ over an interval of length $L=32\lambda$,
where $\lambda$ denotes the Taylor scale of hydrodynamic turbulence. As
detailed in section 3, the corresponding transition PDFs have been obtained
from wind tunnel measurements reaching down to scales $L/2^{9}=\lambda/16$.
Starting from prescribed velocity values at the end points of the interval
$[0,L]$, the velocity $u_{L/2}$ at the mid-point is drawn from the transition
PDF $p_{1}(u_{\mathrm{c}}|u_{\mathrm{l}}=u_{0L},u_{\mathrm{r}}=u_{L})$. The
two velocities $u_{L/4}$ and $u_{3L/4}$ on the next level of refinement are
drawn from
$p_{2}(u_{\mathrm{c}}|u_{\mathrm{l}}=u_{0L},u_{\mathrm{r}}=u_{L/2})$ and
$p_{2}(u_{\mathrm{c}}|u_{\mathrm{l}}=u_{L/2},u_{\mathrm{r}}=u_{L})$,
respectively. The four velocities necessary to refine the field on scale
$r_{3}=L/8$ are then drawn from the corresponding transition PDF $p_{3}$. By
iterating this procedure down to level $9$, we finally obtain the gray time
series shown in figure 1.

Up to now, we have discussed how to generate a field of length $L$ and
resolution $L/2^{N}$. If the scale $L$ is large enough such that field values
separated by $L$ are statistically independent, arbitrarily large fields can
be generated by drawing from the single-point PDF $f(u)$ on scale $L$.
Statistical correlations on scales smaller than $L$ are then subsequently
introduced by MLSR. Regarding the small-scale resolution, the latter is only
limited by the smallest experimentally accessible scale, from which the
transition PDFs are estimated.

###  2.1 Multi-point single-time statistics and relation to Markov processes
in scale

The MLSR discussed in the previous section allows us to generate a synthetic
field of extent $L$ and resolution $\Delta=L/2^{N}$ sampled on $M+1$ points
where $M=2^{N}$. Let us denote the field values at these points sequentially
by $u_{0}\dots u_{M}$ with the corresponding multi-point PDF
$f(u_{0},u_{1},u_{2},u_{3},\dots,u_{M})$. Owing to the structure of the MLSR
algorithm, an ensemble of such synthetically generated fields obeys the multi-
point statistics

$f(u_{0},u_{1},u_{2},u_{3},\dots,u_{M})=f(u_{0},u_{M})\,\prod_{i=1}^{N}\prod_{j=1}^{i}p_{i}\big{(}u_{(2j-1)M/2^{i}}\big{|}u_{(2j-2)M/2^{i}},u_{2jM/2^{i}}\big{)}\,.$
(3)

Consistent with the hierarchical construction of the random fields, this
expression for the multi-point statistics shows explicitly that random
variables can be systematically integrated out level by level. For example, by
integrating out all random variables which have been introduced in the $N$th
refinement step resulting in a field resolution of $\Delta$, the joint PDF
$f(u_{0},u_{2},u_{4},\dots,u_{M})$ corresponding to an ensemble of fields with
a resolution of $2\Delta$ is obtained.

The hierarchical structure of the multi-point statistics (3), in which the
probability of the field value at the center point of an interval is
determined by the values at the two end points, indicates that MLSR generates
a Markov process in scale. As we show by direct calculation in A, Markovian
properties also follow for related statistical quantities such as field
increments. As a further note, if the original time series obeys Markov
properties of the form
$p\left(u_{1}|u_{0},u_{2},u_{3}\right)=p\left(u_{1}|u_{0},u_{2}\right)$, then
equation (3) is an exact representation of the multi-point statistics. For
turbulent flows, this property is an ongoing topic of investigation [25].

A prototypical example for non-Gaussian, yet Markovian, systems are time
series generated by nonlinear Langevin equations. For such time series, the
MLSR algorithm is capable of generating synthetic equivalents that reproduce
the correct $N$-point statistics precisely. As a proof on concept, we
demonstrate this for the two-point statistics in B. However, even if Markovian
properties are mildly violated, the synthetic time series generated by MLSR
still resemble the original statistics closely as illustrated in section 3.

##  3 Application: Turbulence Data

Previous works (e.g. [1, 2, 3]) have empirically established a Markov property
in scale for hydrodynamic turbulence, which lends support to the MLSR approach
proposed here to capture essential statistical features of turbulent fields.
To illustrate the capabilities of the MLSR method for real-world, close-to-
Markovian complex systems, we utilize experimental turbulence data from the
Variable Density Turbulence Tunnel VDTT at the Max Planck Institute for
Dynamics and Self-Organization [26]. This pressurizable wind tunnel uses
sulfur hexafluoride at densities up to 15bar to achieve high Reynolds numbers
due to the low kinematic viscosity of the pressurized working gas. The data we
used to construct the MLSR method is a subset of the data presented in [27].
Turbulence was generated using a passive grid of rectangular grid bars with a
mesh spacing of 18cm and was measured about 40 mesh sizes downstream of the
grid using a nano-scale thermal anemometry probe (NSTAP) developed at
Princeton University [28, 29]. The NSTAP was a 60$\mu$m long hot-wire that
measured one-dimensional velocity signals and was small enough to resolve the
smaller length scales in the VDTT. The data was obtained at a pressure of
8.5bar of sulfur hexafluoride, resulting in a Taylor-scale Reynolds number of
$R_{\lambda}=1030$. The time series was sampled at a rate of 60kHz with noise
being removed at 15kHz using an 8th-order Butterworth filter and contained
$\mathcal{O}\left(10^{10}\right)$ samples of the turbulent velocity,
corresponding to $\mathcal{O}\left(10^{5}\right)$ large eddy turnover times.
Based on the mean speed $U=4.2$m/s and the root mean square velocity
fluctuation $u_{\mathrm{rms}}=0.13$m/s, the turbulence intensity
$u_{\mathrm{rms}}/U$ was 3%. We invoked Taylor’s hypothesis [23] to convert
temporal information into spatial information. The Kolmogorov length scale,
$\eta$, was 34$\mu$m, the integral length scale, $L_{\mathrm{int}}$, was
12.6cm and the Taylor scale, $\lambda$, was 2.2mm.

To obtain the conditional PDFs
$p\left(u_{\mathrm{c}}|u_{\mathrm{l}},u_{\mathrm{r}},r\right)$, we first
subtracted the mean of the signal velocity and computed the three-point PDF
$f(u_{\mathrm{l}},u_{\mathrm{c}},u_{\mathrm{r}};r)$ and the two-point PDF
$f(u_{\mathrm{l}},u_{\mathrm{r}};2r)$, using conventional binning with 100
equidistant bins in all PDF dimensions. Given the substantial length of the
velocity time series, these multi-point PDFs are well-converged. In situations
where one potentially only has access to shorter time series, the estimates of
the PDFs could be improved using e.g. kernel density estimation [30] at a
higher computational cost. To remove spurious noise stemming from divisions of
bin elements with insufficient statistics when computing conditional
statistics, we empirically set all bins in the three-point PDF that contain
less than 10 samples of the turbulent velocity to zero.

Figure 2: Cuts through the conditional PDF
$p(u_{\mathrm{c}}|u_{\mathrm{l}},u_{\mathrm{r}},r)$ for
$u_{\mathrm{r}}/u_{\mathrm{rms}}=0$ and for separations $r$ of 0.5$\lambda$,
4$\lambda$, 32$\lambda$. With decreasing separation, the shape of conditional
PDFs qualitatively changes, showing a more dominant constriction for
$u_{\mathrm{l}}=u_{\mathrm{r}}$. Bottom right: One-dimensional cuts through
the conditional PDFs for $u_{\mathrm{r}}/u_{\mathrm{rms}}=0$ and fixed
$u_{\mathrm{l}}$ indicated by the line color. For large differences
$\left|u_{\mathrm{l}}-u_{\mathrm{r}}\right|$, the PDFs are almost Gaussian,
while significant non-Gaussianity emerges with decreasing difference
$\left|u_{\mathrm{r}}-u_{\mathrm{l}}\right|$. Irrespective of the difference
$\left|u_{\mathrm{l}}-u_{\mathrm{r}}\right|$, all PDFs are centered around the
mean $\left(u_{\mathrm{l}}+u_{\mathrm{r}}\right)/2$.

Two-dimensional cuts through these conditional PDFs computed from the NSTAP
data for $u_{\mathrm{r}}/u_{\mathrm{rms}}=0$ are shown in figure 2. At small
scales $r$, the conditional PDF is noticeably constricted around
$u_{\mathrm{l}}/u_{\mathrm{rms}}=u_{\mathrm{r}}/u_{\mathrm{rms}}=0$ and
substantially widens with increasing difference
$\left|u_{\mathrm{l}}-u_{\mathrm{r}}\right|$. As to be expected, the
constriction around
$u_{\mathrm{l}}/u_{\mathrm{rms}}=u_{\mathrm{r}}/u_{\mathrm{rms}}=0$ weakens
when the separation becomes comparable to the integral length scale such that
the shape of the conditional PDFs, which are centered around
$\left(u_{\mathrm{l}}+u_{\mathrm{r}}\right)/2$, varies less with the
difference $\left|u_{\mathrm{l}}-u_{\mathrm{r}}\right|$. One-dimensional cuts
through the conditional PDFs in figure 2 at a scale of $r=\lambda$ are shown
in the bottom right panel. For small differences
$\left|u_{\mathrm{l}}-u_{\mathrm{r}}\right|$ the cuts through the conditional
PDFs show a noticeably non-Gaussian behavior, while for large differences
$\left|u_{\mathrm{l}}-u_{\mathrm{r}}\right|$, the PDFs become virtually
Gaussian. This emphasizes that the conditional PDFs not only widen with
increasing differences $\left|u_{\mathrm{l}}-u_{\mathrm{r}}\right|$, but
rather change shape and therefore result in a qualitative change of the
underlying dynamics determining $u_{\mathrm{c}}$.

Figure 3: Left column (top to bottom): Mean of the conditional PDF
$p\left(u_{\mathrm{c}}|u_{\mathrm{l}},u_{\mathrm{r}},r\right)$ for separations
$\lambda/4$, 2$\lambda$ and 128$\lambda$ as well as one-dimensional diagonal
cuts through the conditional mean for $u_{\mathrm{l}}=u_{\mathrm{r}}$. The
conditional mean appears to be a linear function of the mean value
$\left(u_{\mathrm{l}}+u_{\mathrm{r}}\right)/2$, with a slope that approaches
one for decreasing separation. Right column (top to bottom): Conditional
variance as well as one-dimensional diagonal cuts for
$u_{\mathrm{l}}=u_{\mathrm{r}}$. In contrast to the conditional mean, the
variance is a slightly asymmetric, non-trivial function of $u_{\mathrm{l}}$
and $u_{\mathrm{r}}$.

This effect can also be seen in the mean and variance of the conditional PDF
in figure 3. While the conditional mean shows a nearly linear behavior across
all separations with a slope that depends on the scale $r$, the variance is
non-trivially depends on the separation and the difference
$\left|u_{\mathrm{l}}-u_{\mathrm{r}}\right|$. However, the deviations from
symmetric behavior are relatively small.

Using these PDFs and the method described in section 2, we can construct a
multi-level refined time series reproducing the key properties of a real
turbulent time series. To initialize the method, we started with a turbulent
time series sampled at $r_{0}=32\lambda$ and utilized a 9-level MLSR process
to reach a final separation $r_{9}=\lambda/16$. Figure 4 shows a direct
comparison of the resulting synthetic time series with an experimentally
obtained turbulent time series sampled at the same separation
$r_{9}=\lambda/16$. Both time series show a visually indistinguishable
behavior, the intermittent character of turbulence is qualitatively well
reproduced by the MLSR time series.

Figure 4: Comparison of time series from wind tunnel data and one obtained
from a 9-level MLSR procedure.

As a quantitative benchmark for the MLSR algorithm, we computed the velocity
increment PDFs $f\left(v;r\right)$, where
$v=u\left(x+r\right)-u\left(x\right)$, at each refinement scale for the
respective separation as well as the velocity increment PDFs for the
experimental turbulence dataset that was originally used to compute the
conditional PDFs. The result is shown in figure 5. At each refinement scale,
the refined time series reproduces the original turbulence PDFs including the
increasingly heavy tails towards smaller scales. As detailed in [31], the
Markovian properties of turbulence start to deteriorate at scales smaller than
$\lambda$, the Taylor scale. While deviations are small, effects of this
violation can supposedly be seen in the tails of the increment PDFs for
smaller scales in figure 5. This is in contrast to the findings of the exact
reproduction for fully Markovian systems in B. It is, however, remarkable how
closely the MLSR-generated time series quantitatively resemble actual fluid
turbulence highlighting its general applicability for systems without ideal
Markov properties.

Figure 5: PDFs of the velocity increments at various separations for each
refinement scale. The curves are staggered for better visibility. The black
line corresponds to PDFs from the original turbulence dataset, the colored
lines to PDFs obtained from the MLSR fields.

##  4 Conclusions and Outlook

We have presented a hierarchical stochastic refinement method that allows to
generate synthetic fields and time series. In this data-driven MLSR, the
synthetic fields (time series) are reconstructed from multi-point (multi-time)
statistics of _real input data_. The method is able to generate synthetic data
which exhibits two-point statistics in good agreement with the original data.
We demonstrated the data-driven MLSR at the example of turbulent wind tunnel
data as well as for nonlinear Langevin equations.

Due to the low requirements of the method, i.e. the determination of the
three-point PDF for different refinement levels, as well as due to the fast
production of surrogate data, the data-driven MLSR appears to be a promising
method for a broad range of applications. Synthetic fields generated by MLSR,
for example, might be used as input to simulate power output fluctuations of
wind parks, the fatigue loads on single wind turbines [32], to test the
stability of power systems [33], as well as to simulate the propagation of
cosmic rays [34, 35, 36].

In certain situations, sufficient statistics to resolve the three-point PDFs
required for the data-driven MLSR method might not be available and
alternative approaches are necessary. An extension of the approach to the
generation of synthetic data based on phenomenological models of turbulence
[37, 38] without the need for measuring three-point PDFs is subject of ongoing
work.

## Acknowledgments

We are grateful to Eberhard Bodenschatz for the opportunity of using the data
acquired during the doctoral studies of MS at the Max Planck Institute for
Dynamics and Self-Organization for this project as well as to Gregory P.
Bewley for participating in the data acquisition and for fruitful discussions.
We thank Joachim Peinke for valuable discussions.

MW is supported by the Max Planck Society. JF acknowledges funding from the
Humboldt Foundation within a Feodor-Lynen fellowship and also benefitted from
the financial support of the Project IDEXLYON of the University of Lyon in the
framework of the French program “Programme Investissements d’Avenir”
(ANR-16-IDEX-0005). MS acknowledges support from the Deutsche
Forschungsgemeinschaft under grant no. 396632606.

##  Appendix A Relation of MLSR to Markov processes in scale

To elucidate the relation of the MLSR approach to Markov processes in scale in
terms of increments, let us consider one particular “branch” of random
variables in (3) which connects the largest scales to the smallest scales.
Without loss of generality, this can be achieved by setting $j=1$ in (3) which
results in the reduced multi-point PDF

$f(u_{0},u_{1},u_{2},u_{4},u_{8},u_{16},\dots,u_{M})=f(u_{0},u_{M})\,\prod_{i=1}^{N}p_{i}\big{(}u_{M/2^{i}}\big{|}u_{0},u_{M/2^{(i-1)}}\big{)}\,.$
(4)

The joint PDF of increments $v_{i}=u_{M/2^{i}}-u_{0}$ on scales $r_{i}$
together with $u_{0}$ are obtained from the reduced multi-point PDF (4)
according to

$\displaystyle f(v_{0},v_{1},v_{2},v_{3},\dots,v_{N},u_{0})$ (5)
$\displaystyle=$
$\displaystyle\prod_{i=0}^{N}\int\mathrm{d}u_{M/2^{i}}\,\delta\left(u_{M/2^{i}}-u_{0}-v_{i}\right)f(u_{0},u_{1},u_{2},u_{4},u_{8},u_{16},\dots,u_{M})\,.$

After inserting (4) into (5) the integration can be carried out and we obtain

$f(v_{0},v_{1},v_{2},v_{3},\dots,v_{N},u_{0})=f(u_{0},v_{0}+u_{0})\,\prod_{i=1}^{N}p_{i}\big{(}v_{i}+u_{0}\big{|}u_{0},v_{i-1}+u_{0}\big{)}\,.$
(6)

Similar to the multi-point PDFs, the increment variables can be integrated out
scale by scale. For example, the joint PDF which includes all increments but
the smallest,

$f(v_{0},v_{1},v_{2},v_{3},\dots,v_{N-1},u_{0})=f(u_{0},v_{0}+u_{0})\,\prod_{i=1}^{N-1}p_{i}\big{(}v_{i}+u_{0}\big{|}u_{0},v_{i-1}+u_{0}\big{)}\,,$
(7)

can be straightforwardly obtained by integrating over $v_{N}$. As a
consequence, the conditional increment PDF fulfills a Markov property, which
can be seen by noticing that (7) and (8) only differ by the transition PDF
$p_{N}$. As a result we obtain

$p(v_{N}|v_{0},v_{1},\dots,v_{N-1},u_{0})=\frac{f(v_{0},v_{1},v_{2},v_{3},\dots,v_{N},u_{0})}{f(v_{0},v_{1},v_{2},v_{3},\dots,v_{N-1},u_{0})}=p_{N}\big{(}v_{N}+u_{0}\big{|}u_{0},v_{N-1}+u_{0}\big{)}\,.$
(8)

This relation shows that only the transition PDF $p_{N}$ of the $N$th
refinement level is needed to fully specify increment PDF conditional on all
other increments and $u_{0}$. In other words, the conditional PDF
$p(v_{N}|v_{0},v_{1},\dots,v_{N-1},u_{0})$ indeed is independent of
$v_{0},\dots,v_{N-2}$ such that the generalized Markov property111The Markov
property is generalized in the sense that the conditional PDFs maintain a
dependence on $u_{0}$. The necessity of including $u_{0}$ has also been
discussed in a related approach presented in [3].

$p(v_{N}|v_{0},v_{1},\dots,v_{N-1},u_{0})=p(v_{N}|v_{N-1},u_{0})$ (9)

holds. Of course, this argument holds on all larger scales as well, showing
that MLSR generates a Markov process in scale.

##  Appendix B Proof-of-Concept: Nonlinear Langevin Equations

The method described in section 2 can, given Markov properties, produces
synthetic time series with correct statistics from any experimental or
numerical time series with sufficient statistics to resolve the three-point
PDFs $f(u_{\mathrm{l}},u_{\mathrm{c}},u_{\mathrm{r}};r_{i})$ at a hierarchy of
scales $r_{i}$. From these three-point PDFs and the corresponding two-point
PDFs at each scale, $f(u_{\mathrm{l}},u_{\mathrm{r}};2r_{i})$, one can then
compute the conditional PDFs,
$p_{i}(u_{\mathrm{c}}|u_{\mathrm{l}},u_{\mathrm{r}})$, through equation (1)
and use those to refine a coarsely sampled time series at scale $r_{0}$ down
to scale $r_{N}$.

Figure 6: Implementation of the MLSR method to the nonlinear Langevin equation
(10) with $\gamma=1$, $c_{1}=5$, $\sigma=1$ and $c_{2}=1$. The top plot shows
a time series of the simulated system, the bottom plot a 11-level refinement
from the MLSR method.

The quality of the refinement depends on the statistical properties of the
original time series as well as on the convergence of the measured PDFs. In
general, the further the original time series deviates from exhibiting
Markovian properties, the less accurate the reconstruction becomes.

As a proof-of-concept test case, let us consider a prototypical nonlinear
Langevin equation of the form

$\mathrm{d}u=-\gamma
u^{c_{1}}\mathrm{d}t+\sqrt{\left(\sigma+c_{2}u\right)^{2}}\mathrm{d}W_{t}.$
(10)

Here, $\mathrm{d}W_{t}$ is the increment of a Wiener process in time, $\gamma$
is a parameter controlling the strength of the (nonlinear) drift term, and
$\sigma$ a diffusion parameter. The parameters $c_{1}$ and $c_{2}$ determine
the degree of nonlinearities of the respective terms and give rise to non-
Gaussian behavior. For $c_{1}=1$ and $c_{2}=0$ this equation reduces to the
standard Ornstein-Uhlenbeck process with linear drift and constant diffusion
that results in Gaussian solutions. The system can be directly simulated using
a standard Euler-Maruyama method [39]. To obtain sufficient data for the
three-point PDFs, we compute a total of $1.2\cdot 10^{10}$ data points with a
time step of $\Delta t=10^{-3}$ for various choices of the parameters. We then
apply the MLSR method to refine initial conditions at scale
$\tau_{0}=2^{11}\Delta t$ down to scale $\tau_{11}=\tau_{0}/2^{11}$. Such a
refinement is shown in figure 6.

Figure 7: Increment PDFs of simulated (black solid lines) Langevin equation
(10) and the respective MLSR reproductions (blue dashed lines). The curves are
staggered for better visibility corresponding to scales $\tau=\tau_{0}$
(bottommost) through $\tau=\tau_{7}$ (topmost). (Left): The parameters are
$\gamma=1$, $c_{1}=1$, $\sigma=1$ and $c_{2}=1$. This choice results in a
heavy-tailed increment distribution at all scales that is correctly reproduced
by the MLSR algorithm. (Right): The parameters are $\gamma=1$, $c_{1}=5$,
$\sigma=1$ and $c_{2}=1$. The sub-Gaussian cores of the increment PDFs
produced by the strongly nonlinear drift as well as the heavy tails are well-
reproduced.

The properties of the MLSR reproduction are virtually indistinguishable from
the original, simulated time series. This can also be tested quantitatively by
the means of the scale-dependent PDFs $f\left(w;\tau\right)$ of the increment
$w=u\left(t+\tau\right)-u\left(t\right)$ at a given scale $\tau$. This is
shown in figures 7 for two sets of parameters.

## References

## References

  * [1] Friedrich R and Peinke J 1997 Phys. Rev. Lett. 78 863–866
  * [2] Renner C, Peinke J and Friedrich R 2001 J. Fluid Mech. 433 383–409
  * [3] Stresing R and Peinke J 2010 New J. Phys. 12(10) 103046 http://stacks.iop.org/1367-2630/12/i=10/a=103046
  * [4] Strumik M and Macek W 2008 Phys. Rev. E 78 026414 https://link.aps.org/doi/10.1103/PhysRevE.78.026414
  * [5] Sherry C, Barrow D and Klemm W 1982 Brain Research Bulletin 8(2) 163 – 169 http://www.sciencedirect.com/science/article/pii/0361923082900430
  * [6] Friedrich R, Peinke J and Renner C 2000 Phys. Rev. Lett. 84 5224–5227 https://link.aps.org/doi/10.1103/PhysRevLett.84.5224
  * [7] Bassler K E, McCauley J L and Gunaratne G H 2007 PNAS 104(44) 17287–17290 https://www.pnas.org/content/104/44/17287
  * [8] Nawroth A and Peinke J 2006 PLA 360(2) 234 – 237 http://www.sciencedirect.com/science/article/pii/S0375960106012837
  * [9] Behnken C, Wächter M and Peinke J 2020 Wind Energy Science 5(3) 1211–1223
  * [10] Meneveau C and Katz J 2000 Annu. Rev. Fluid Mech. 32 1–32
  * [11] Sagaut P 2006 Large Eddy Simulation for Incompressible Flows (Springer, Berlin/Heidelberg, Germany)
  * [12] Shaw R A 2003 Annu. Rev. Fluid Mech. 35 183–227
  * [13] Grabowski W W and Wang L P 2013 Annu. Rev. Fluid Mech. 45 293–324
  * [14] Pumir A and Wilkinson M 2016 Annu. Rev. Cond. Ma. P. 7 141–170
  * [15] Meneveau C and Poinsot T 1991 Combust. Flame 86 311 – 332 http://www.sciencedirect.com/science/article/pii/001021809190126V
  * [16] Dopazo C, Cifuentes L, Martin J and Jimenez C 2015 Combust. Flame 162 1729–1736
  * [17] Gifford F A 1982 Atmospheric Dispersion Models for Environmental Pollution Applications (Boston, MA: American Meteorological Society) https://doi.org/10.1007/978-1-935704-23-2_2
  * [18] Vallis G K 2017 Atmospheric and oceanic fluid dynamics (Cambridge University Press, Cambridge, United Kingdom)
  * [19] Mücke T, Kleinhans D and Peinke J 2011 Wind Energy 14(2) 301–316 https://onlinelibrary.wiley.com/doi/abs/10.1002/we.422
  * [20] Spencer M D, Stol K A, Unsworth C P, Cater J E and Norris S E 2013 Wind Energy 16(3) 417–434 https://onlinelibrary.wiley.com/doi/abs/10.1002/we.1501
  * [21] Foresti L, Tuia D, Kanevski M and Pozdnoukhov A 2011 SERRA 25(1) 51–66 https://doi.org/10.1007/s00477-010-0405-0
  * [22] Heinermann J and Kramer O 2016 Renewable Energy 89 671 – 679 http://www.sciencedirect.com/science/article/pii/S0960148115304894
  * [23] Taylor G I 1938 Proc. Roy. Soc. Lond. A 164 476–490
  * [24] von Neumann J 1951 in A S Householder, G E Forsythe and H H Germond, eds, Monte Carlo Method vol 12 of National Bureau of Standards Applied Mathematics Series (Washington, DC: US Government Printing Office) chap 13, pp 36–38
  * [25] Peinke J personal communication
  * [26] Bodenschatz E, Bewley G P, Nobach H, Sinhuber M and Xu H 2014 Rev. Sci. Instrum. 85 093908
  * [27] Sinhuber M, Bewley G P and Bodenschatz E 2017 Phys. Rev. Lett. 119 134502
  * [28] Bailey S C C, Kunkel G J, Hultmark M, Vallikivi M, Hill J P, Meyer K A, Tsay C, Arnold C B and Smits A J 2010 J. Fluid Mech. 663 160–179
  * [29] Vallikivi M, Hultmark M, Bailey S and Smits A 2011 Exp. Fluids 51 1521–1527
  * [30] Gramacki A 2018 Nonparametric Kernel Density Estimation and Its Computational Aspects (Springer International Publishing)
  * [31] Lück S, Renner C, Peinke J and Friedrich R 2006 Phys. Lett. Sect. A Gen. At. Solid State Phys. 359(5) 335–338
  * [32] Lind P G, Herráez I, Wächter M and Peinke J 2014 Energies 7(12) 8279–8293
  * [33] Schmietendorf K, Peinke J, Friedrich R and Kamps O 2014 EUR PHYS J-SPEC TOP 223(12) 2577–2592
  * [34] Schlickeiser R 2015 PoP 22(9) 091502
  * [35] Zweibel E G 2013 Phys. Plasmas 20(5) 055501
  * [36] Reichherzer P, Becker Tjus J, Zweibel E G, Merten L and Pueschel M J 2020 Monthly Notices of the Royal Astronomical Society 498(4) 5051–5064 (Preprint https://academic.oup.com/mnras/article-pdf/498/4/5051/33806251/staa2533.pdf) https://doi.org/10.1093/mnras/staa2533
  * [37] Friedrich J and Grauer R 2020 Atmosphere 11(9) https://www.mdpi.com/2073-4433/11/9/1003
  * [38] Friedrich J, Gallon S, Pumir A and Grauer R 2020 Phys. Rev. Lett. 125(17) 170602
  * [39] Kloeden P and Platen E 1992 Numerical Solution of Stochastic Differential Equations (Springer-Verlag, Berlin)