# Multifractal Fractional Ornstein-Uhlenbeck Processes

Laurent Chevillard  Univ Lyon, Ens de Lyon, Univ Claude Bernard, CNRS,
Laboratoire de Physique, 46 allée d’Italie F-69342 Lyon, France
<EMAIL_ADDRESS>,  Marc Lagoin  and  Stéphane Gaetan Roux

###### Abstract.

The Ornstein-Uhlenbeck process can be seen as a paradigm of a finite-variance
and statistically stationary rough random walk. Furthermore, it is defined as
the unique solution of a Markovian stochastic dynamics and shares the same
local regularity as the one of the Brownian motion. Based on previous works,
we propose to include in the framework of one of its generalization, the so-
called fractional Ornstein-Uhlenbeck process, some Multifractal corrections,
using a Gaussian Multiplicative Chaos. The aforementioned process, called a
Multifractal fractional Ornstein-Uhlenbeck process, is a statistically
stationary finite-variance process. Its underlying dynamics is non-Markovian,
although non-anticipating and causal. The numerical scheme and theoretical
approach are based on a regularization procedure, that gives a meaning to this
dynamical evolution, which unique solution converges towards a well-behaved
stochastic process.

###### Key words and phrases:

Fractional Brownian Motion, Ornstein-Uhlenbeck processes, Multiplicative Chaos

###### 2010 Mathematics Subject Classification:

60G22, 60H05

##  1\. Introduction

In many situations, for example some physical phenomena governed by thermal
fluctuations such as the position of small particles suspended in a quiescent
fluid, the standard Brownian motion, or Wiener process, $W(t)$ is called to
model and give account of certain quantities of interest. As often required by
the nature of the situation, the phenomenon itself is statistically stationary
and of finite-variance, therefore calling for the Ornstein-Uhlenbeck (OU)
process. It is defined as the unique solution of the stochastic dynamics given
by

(1.1) $\displaystyle dX(t)=-\frac{1}{T}X(t)dt+dW(t),$

where the increment $dW(t)$ of the Wiener process $W$ can be viewed as a
Gaussian white noise, and which meaning is clear as a random distribution. As
a continuous-time stochastic process, $X(t)$ is nowhere differentiable, and we
can say that sample paths of $X$, as those of $W$, are H$\ddot{\text{o}}$lder
continuous of any order strictly less than $H=1/2$, where $H$ is the Hurst
exponent.

Other physical phenomena, for example fluid turbulence [19], or some time
series encountered in finance [21, 5], exhibit in a statistically averaged
sense a different Hurst exponent $H\in]0,1[$, a local regularity that can be
reproduced by fractional Brownian motions [22]. Once again, as required for
instance by the physics of turbulence which suggests a roughness of order
$H\approx 1/3$ [13], such phenomena are eventually statistically stationary
and are of finite variance. In this case, using fractional Ornstein-Uhlenbeck
(fOU) processes [9] to model them seems more appealing. They could be define
as the unique solution of the following dynamics

(1.2) $\displaystyle
dX_{H}(t)=-\frac{1}{T}X_{H}(t)dt+\text{``}dW_{H}(t)\text{''},$

where the random noise $dW_{H}$ entering in Eq. 1.2 is indicated between
quotes to remind that its meaning has to be clarified. Whereas fractional
motions can be written as stochastic integrals in a natural way [27, 24, 31]
for any $H\in]0,1[$, the meaning pathwise of the random noise $dW_{H}(t)$, and
consequently the meaning of the dynamics of $X_{H}$ (Eq. 1.2), is not obvious
when $H<1/2$. As proposed in Ref. [10], a regularization procedure over a
small timescale, $\epsilon>0$, gives both a precise and explicit meaning to
this random measure $dW_{\epsilon,H}(t)$. Accordingly, the unique solution
$X_{\epsilon,H}(t)$ of Eq. 1.2 once forced by such a random noise is shown in
[10] to converge, as $\epsilon\to 0$, towards the canonical fOU process of
Ref. [9].

As we can see, this regularization procedure over $\epsilon$ is necessary,
when $H<1/2$, to define pathwise the dynamics schematically presented in Eq.
1.2. This regularization turns out to be crucial when multifractal
corrections, usually referred as the intermittency phenomenon in the
literature of turbulence [13], are added to the picture while considering a
Gaussian multiplicative chaos [20, 18, 29]. It can be formally obtained as

(1.3) $\displaystyle
M_{\gamma}(t)\;\text{``}=\text{''}\;e^{\gamma\widetilde{X}(t)},$

where the Gaussian field $\widetilde{X}$ entering in Eq. 1.3 is asked to be
logarithmically correlated, and $\gamma\in\mathbb{R}$ is a free parameter.
Once again, the equality given in Eq. 1.3 is put in quotes because the meaning
of the exponential of such a field, necessarily being of infinite variance,
has yet to be clarified. In this article, we will show how to give a precise
meaning to $M_{\gamma}$ (Eq. 1.3) while still relying on a regularization
procedure over $\epsilon$, as it is usually done [29]. Moreover, we will see
that the Gaussian process $\widetilde{X}$ can actually be defined as a fOU
process, solution of the stochastic dynamics provided in Eq. 1.2, for a
vanishing Hurst exponent $H=0$. This approach was already proposed and used in
different setups to model some aspects of the random nature of fluid
turbulence [10, 25, 26, 32, 12]. This method of construction differs from
other ones found in the literature [14, 23, 15]. In particular, it includes
the notion of causality, and its respective stochastic evolution, to this
logarithmically correlated Gaussian field.

Finally, we propose a way to incorpore into the dynamics of the fOU process
(Eq. 1.2) a Multiplicative Chaos (Eq. 1.3), and define respectively the
stochastic evolution

(1.4) $\displaystyle
dX_{H,\gamma}(t)=-\frac{1}{T}X_{H,\gamma}(t)dt+\text{``}dW_{H,\gamma}(t)\text{''},$

which eventually makes sense only up to a regularization scale $\epsilon>0$.
Taking the limit $\epsilon\to 0$, the asymptotic process obtained while
considering the unique solution of Eq. 1.2, will be referred as the
Multifractal fractional Ornstein-Uhlenbeck (MfOU) process. It can be seen as a
generalization of the Multifractal random walk of Ref. [3] to a causal
framework, and for any $H\in]0,1[$ instead of the unique value $H=1/2$.

### Organization of the paper.

In Section 2, we present the general setup of the article, and especially the
regularized version of the stochastic evolution proposed in Eq. 1.4. We
include several remarks on the various ingredients entering in the
construction of this dynamics. Section 3 is devoted to the statement of our
results. In particular, we list in Proposition 3.1 some important marginals of
the MfOU process $X_{H,\gamma}$, including variance and high order moments of
its increments. In Section 4, we propose a method to synthesize a
statistically stationary solution of the dynamics (Eq. 1.4) in a periodical
framework, and compare our theoretical predictions to our numerical
estimations. We gather the proofs of Proposition 3.1 in Section 5.

##  2\. Setup, notations and remarks

######  Definition 2.1.

Consider a Wiener processes $W(t)$, note by $dW(t)$ its increments over $dt$
at time $t$. In the sequel, we will be studying the unique statistically
stationary solution $X_{\epsilon,H,\gamma}(t)$ of the following linear
stochastic equation

(2.1)
$dX_{\epsilon,H,\gamma}(t)=-\frac{1}{T}X_{\epsilon,H,\gamma}(t)dt+dW_{\epsilon,H,\gamma}(t),$

where enters a large timescale $T>0$ and a small regularizing timescale
$\epsilon>0$. For a given Hurst exponent $H\in]0,1[$, and an additional free
parameter $\gamma\in\mathbb{R}$, the dynamical evolution proposed in Eq. 2.1
is maintained by a random noise $dW_{\epsilon,H,\gamma}(t)$ that reads

(2.2)
$dW_{\epsilon,H,\gamma}(t)=\omega_{\epsilon,H,\gamma}(t)dt+\epsilon^{H-\frac{1}{2}}M_{\epsilon,\gamma}(t)dW(t),$

where we have introduced the statistically stationary random field

(2.3)
$\omega_{\epsilon,H,\gamma}(t)=\left(H-\frac{1}{2}\right)\int_{-\infty}^{t}(t-s+\epsilon)^{H-\frac{3}{2}}M_{\epsilon,\gamma}(s)dW(s).$

The random noise $dW_{\epsilon,H,\gamma}$ (Eq. 2.2) and field
$\omega_{\epsilon,H,\gamma}$ (Eq. 2.3) involve respectively a local version of
and a linear operation on the white noise $dW$ weighted by a positive random
field $M_{\epsilon,\gamma}(t)$ given by

(2.4)
$M_{\epsilon,\gamma}(t)=e^{\gamma\widetilde{X}_{\epsilon}(t)-\gamma^{2}\mathbb{E}\left(\widetilde{X}_{\epsilon}^{2}\right)}.$

Finally, the Gaussian process $\widetilde{X}_{\epsilon}$ that enters in the
definition of the positive random weight $M_{\epsilon,\gamma}$ (Eq. 2.4) can
be seen as the unique statistically stationary solution of the dynamics
provided in Eq. 2.1 using $H=0$, $\gamma=0$ and an independent white noise
$d\widetilde{W}(t)$ instead of $dW(t)$.

######  Remark 2.2.

As we will see, the parameter $\gamma$ eventually governs entirely the level
of non Gaussianity and multifractality of the process
$X_{\epsilon,H,\gamma}(t)$. To this regard, for $\gamma\neq 0$, the dynamics
provided in Eq. 2.1 can be seen as a generalization of the one of the
regularized Fractional Ornstein-Uhlenbeck process detailed in Ref. [10]. Using
current notations, the unique statistically stationary solution of Eq. 2.1
with $\gamma=0$ can be conveniently written as

(2.5)
$X_{\epsilon,H,0}(t)=\int_{-\infty}^{t}e^{-\frac{t-s}{T}}dW_{\epsilon,H,0}(s).$

It is a finite variance Gaussian process for any $H\in]0,1[$ and finite
$\epsilon>0$, and it remains so as $\epsilon\to 0$. We note by $X_{H,0}$ the
Gaussian stochastic process which marginals coincide with the limiting values
as $\epsilon\to 0$ of those of $X_{\epsilon,H,0}$, and we call it a fractional
Ornstein-Uhlenbeck (fOU) process [9]. We have

(2.6) $\mathbb{E}\left(X_{H,0}^{2}\right)\equiv\lim_{\epsilon\to
0}\mathbb{E}\left(X_{\epsilon,H,0}^{2}\right)=\frac{T^{2H}\left[\Gamma\left(H+\frac{1}{2}\right)\right]^{2}}{2\sin(\pi
H)},$

where enters the function $\Gamma(z)=\int_{0}^{\infty}x^{z-1}e^{-x}dx$. More
generally, since $X_{\epsilon,H,0}$ is a zero average Gaussian process, it is
fully characterized by its covariance function which can be explicitly derived
and behaves in the limit of vanishing regularizing scale $\epsilon\to 0$ as,
considering without loss of generality $\tau\geqslant 0$,

(2.7)
$\displaystyle\mathbb{E}\left(X_{H,0}(t)X_{H,0}(t+\tau)\right)\equiv\lim_{\epsilon\to
0}\mathbb{E}\left(X_{\epsilon,H,0}(t)X_{\epsilon,H,0}(t+\tau)\right)$ (2.8)
$\displaystyle=\frac{1}{2\sin(\pi
H)}\frac{\left[\Gamma\left(H+\frac{1}{2}\right)\right]^{2}}{\Gamma\left(2H\right)}\int_{0}^{\infty}\left[\frac{1}{2}e^{-\frac{\tau+h}{T}}-\frac{\tau-h}{2|\tau-h|}e^{-\frac{|\tau-h|}{T}}\right]h^{2H-1}dh,$
(2.9)
$\displaystyle=T^{2}\left[\Gamma\left(H+\frac{1}{2}\right)\right]^{2}\int_{\mathbb{R}}e^{2i\pi\omega\tau}\frac{|2\pi\omega|^{1-2H}}{1+4\pi^{2}\omega^{2}T^{2}}d\omega.$

Notice that the expression of the correlation function (Eq. 2.8), derived in
Ref. [10], is especially well suited to study its decay towards 0 at large
arguments, in an alternative way as it is proposed in [9], whereas its
spectral form (Eq. 2.9), that can be derived and justified using more general
expressions of Ref. [10], clarifies that indeed, the asymptotic process
$X_{H,0}$ coincides with the fractional Ornstein-Uhlenbeck process of [9], and
differs from the so-called tempered fractional Brownian motion of Ref. [7].

Finally, as a consequence of Eq. 2.8, it is straightforward to see that
fractional Brownian motions [22] and fractional Ornstein-Uhlenbeck processes
[9] share the same local regularity, as it is quantified by the second moment
of the increments. Defining the increment as
$\delta_{\tau}X_{\epsilon,H,0}(t)\equiv
X_{\epsilon,H,0}(t+\tau)-X_{\epsilon,H,0}(t)$, we have

(2.10) $\displaystyle\mathbb{E}\left(\delta_{\tau}X_{H,0}\right)^{2}$
$\displaystyle\equiv\lim_{\epsilon\to
0}\mathbb{E}\left(\delta_{\tau}X_{\epsilon,H,0}\right)^{2}$
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0}}\frac{1}{\sin(\pi
H)}\frac{\left[\Gamma\left(H+\frac{1}{2}\right)\right]^{2}}{\Gamma\left(2H+1\right)}|\tau|^{2H}.$

######  Remark 2.3.

Let us now comment on the Gaussian field $\widetilde{X}_{\epsilon}$ entering
in the construction of the regularized Multifractal fractional Ornstein-
Uhlenbeck process defined in 2.1. As already mentioned,
$\widetilde{X}_{\epsilon}$ can be viewed as a regularized fOU process, as it
is recalled in remark 2.2, with the boundary value $H=0$. More precisely, it
can be written as the following stochastic integral

(2.11)
$\widetilde{X}_{\epsilon}(t)=\int_{-\infty}^{t}e^{-\frac{t-s}{T}}d\widetilde{W}_{\epsilon,0,0}(s),$

where

(2.12)
$d\widetilde{W}_{\epsilon,0,0}(t)=-\frac{dt}{2}\int_{-\infty}^{t}(t-s+\epsilon)^{-\frac{3}{2}}d\widetilde{W}(s)+\epsilon^{-\frac{1}{2}}d\widetilde{W}(t).$

As we will see, the Gaussian process $\widetilde{X}_{\epsilon}$ eventually
shares similar statistical properties as those of the processes constructed in
Refs. [14, 23] in the asymptotic limit of vanishing $\epsilon$, although we
recall that it is moreover defined as a stationary solution of a well-posed
stochastic evolution.

In this case, the finiteness of variance is not guaranteed as the regularizing
scale $\epsilon$ goes to zero, and indeed, it is shown in Refs. [10, 26, 1]
that instead,

(2.13)
$\mathbb{E}\left(\widetilde{X}_{\epsilon}^{2}\right)\equiv\mathbb{E}\left(X_{\epsilon,0,0}^{2}\right)\mathrel{\mathop{\kern
0.0pt\sim}\limits_{\epsilon\to 0}}\ln\left(\frac{1}{\epsilon}\right).$

Whereas the variance of $\widetilde{X}_{\epsilon}$ (Eq. 2.13) diverges
logarithmically with the regularizing parameter $\epsilon$, its covariance
remains bounded, and we can write, as a direct application of Eqs. 2.8 and
2.9, for any $\tau>0$,

(2.14)
$\displaystyle\mathbb{E}\left(\widetilde{X}(t)\widetilde{X}(t+\tau)\right)$
$\displaystyle\equiv\lim_{\epsilon\to
0}\mathbb{E}\left(X_{\epsilon,0,0}(t)X_{\epsilon,0,0}(t+\tau)\right)$ (2.15)
$\displaystyle=\int_{0}^{\infty}\left[\frac{1}{2}e^{-\frac{\tau+h}{T}}-\frac{\tau-h}{2|\tau-h|}e^{-\frac{|\tau-h|}{T}}\right]h^{-1}dh,$
(2.16) $\displaystyle=\pi
T^{2}\int_{\mathbb{R}}e^{2i\pi\omega\tau}\frac{|2\pi\omega|}{1+4\pi^{2}\omega^{2}T^{2}}d\omega.$

Starting from Eq. 2.15, we recover that, once the limit $\epsilon\to 0$ has
been taken, the variance of $\widetilde{X}$ is infinite while noticing that

(2.17)
$\mathbb{E}\left(\widetilde{X}(t)\widetilde{X}(t+\tau)\right)\mathrel{\mathop{\kern
0.0pt\sim}\limits_{\tau\to 0}}\ln\left(\frac{1}{|\tau|}\right).$

Because of the logarithmic behavior of the covariance function of
$\widetilde{X}$ for small arguments (Eq. 2.17) and its decay towards 0 at
large arguments $\tau\gtrapprox T$, we will find convenient, and useful for
subsequent calculations, to write in a formal way

(2.18)
$\mathbb{E}\left(\widetilde{X}(t)\widetilde{X}(t+\tau)\right)=\ln_{+}\left(\frac{T}{|\tau|}\right)+g(\tau),$

with $\ln_{+}(s)=\max(\ln s,0)$ and $g$ a continuous, even and bounded
function of its argument, as it is done in Ref. [29]. We could get an exact
expression for the bounded function $g$ entering in Eq. 2.18 using in
particular Eq. 2.15, although, as we will see, it will turn out to be
sufficient to specify its value at the origin to understand with precision the
statistical properties of the MfOU process at small scales. We find

(2.19) $g(0)=\int_{0}^{\infty}\ln(h)e^{-h}dh\approx-0.577216,$

which is known as minus the Euler-Mascheroni constant.

######  Remark 2.4.

The last ingredient entering in the construction of the MfOU is the positive
random weight $M_{\epsilon,\gamma}$ defined in Eq. 2.4. It ressembles a
regularized version of a Gaussian Multiplicative Chaos, as it is called in the
probability theory literature [29], although, the way it is renormalized makes
it converging locally as $\epsilon\to 0$ to a random distribution, which is
zero-average and unit-variance. In the best situations, where subsequent
integrals involve a deterministic function $f$ that ensures dominated
convergence, we will write, making use of the expression of the covariance of
$\widetilde{X}$ (Eq. 2.18),

(2.20) $\displaystyle\lim_{\epsilon\to
0}\;\mathbb{E}\;\bigg{(}\int_{\mathbb{R}}f(t)M_{\epsilon,\gamma}(t)$
$\displaystyle dW(t)\bigg{)}^{2n}=$
$\displaystyle\frac{(2n)!}{2^{n}n!}T^{2\gamma^{2}n(n-1)}\int_{\mathbb{R}^{n}}\prod_{i=1}^{n}f^{2}(t_{i})\prod_{i<j=1}^{n}\frac{e^{4\gamma^{2}g(t_{i}-t_{j})}}{|t_{i}-t_{j}|_{+}^{4\gamma^{2}}}\prod_{i=1}^{n}dt_{i},$

where $|.|_{+}=\exp(\ln_{+}|.|)$ and for a finite range of values of
$\gamma^{2}$, depending on the regularity of $f$ in the vicinity of the origin
and on the requirement of the finiteness of this moment. As we will see, the
very precise behavior of the random noise $M_{\epsilon,\gamma}(t)dW(t)$ when
integrated against test functions, such as it is depicted in Eq. 2.20, is
eventually responsible of Multifractal corrections to the monofractal Gaussian
process fOU.

Let us mention that defining a Multiplicative Chaos (Eq. 2.4) in such a causal
framework is far from anecdotal. In particular, causality is crucial when
modeling several aspects of fluid turbulence [10, 26]. Making use of Itô’s
lemma, the regularization procedure allows to derive a stochastic evolution
for the random weight $M_{\epsilon,\gamma}$ (Eq. 2.4) itself, and we obtain

(2.21)
$\displaystyle\frac{dM_{\epsilon,\gamma}(t)}{M_{\epsilon,\gamma}(t)}=-\frac{1}{T}\left[\ln
M_{\epsilon,\gamma}(t)+\gamma^{2}\mathbb{E}\left(\widetilde{X}_{\epsilon}^{2}\right)+\gamma\widetilde{\omega}_{\epsilon,0,0}(t)+\frac{\gamma^{2}}{2\epsilon}\right]dt+\gamma\epsilon^{-\frac{1}{2}}d\widetilde{W}(t),$

where the causal field $\widetilde{\omega}_{\epsilon,0,0}(t)$ corresponds to
the expression of Eq. 2.3, with $H=0$, $\gamma=0$, and using the increment of
the independent Wiener process $\widetilde{W}$ instead of $W$. It is
interesting to observe that each of the terms entering in the RHS of Eq. 2.21
diverges pathwise as $\epsilon\to 0$, whereas the variance of the stationary
solution remains bounded, with $\mathbb{E}M^{2}_{\epsilon,\gamma}=1$. Notice
also that the dynamical evolution of $M_{\epsilon,\gamma}$ (Eq. 2.21) can be
seen as a stochastic Volterra system [16].

Finally, let us remark that the random weight $M_{\epsilon,\gamma}$ (Eq. 2.4),
which is a log-normal random field since it is taken as the exponential of a
Gaussian process, can be seen as a particular case of the more general class
of log-infinitely divisible measures [6, 30, 4, 8, 28, 29].

##  3\. Statement of the results

######  Proposition 3.1.

The unique solution of the stochastic differential equation Eq. 2.1 converges
at large time towards a statistically stationary process
$X_{\epsilon,H,\gamma}(t)$ that we can write

(3.1)
$X_{\epsilon,H,\gamma}(t)=\int_{-\infty}^{t}e^{-\frac{t-s}{T}}dW_{\epsilon,H,\gamma}(s),$

where the integration is performed over the random measure
$dW_{\epsilon,H,\gamma}$ defined in Eq. 2.2. In this statistically steady
range, it is a zero-average process for any $H$ and $\epsilon>0$. Note by
$X_{H,\gamma}(t)$ the process which marginals coincide with the limiting
behavior of those of $X_{\epsilon,H,\gamma}(t)$ as $\epsilon\to 0$, and call
it a Multifractal fractional Ornstein-Uhlenbeck (MfOU) process. For any
$H\in]0,1[$, independently of $\gamma$, the covariance of $X_{H,\gamma}$
coincides with the one of the underlying fractional Ornstein-Uhlenbeck (fOU)
process $X_{H,0}$. In particular, its variance is finite (Eq. 2.6) and its
covariance can be derived in an exact fashion (Eq. 2.8).

More generally, considering the integer $n\geqslant 1$, and for any
$H\in]0,1[$ and $\gamma^{2}<\min(1/4,H/(n-1))$, the MfOU process has a finite
moment of order $2n$, that is

(3.2) $\mathbb{E}\left(X_{H,\gamma}^{2n}\right)<\infty,$

and we can obtain their exact expression as multiple integrals (see Section 5
devoted to proofs, and in particular Eqs. 5.30, 5.31 and 5.37). Furthermore,
the MfOU process exhibits at small scales a multifractal behavior, as it can
be characterized in a simple fashion by the behavior at small scales of the
moments its increments. We have, for an integer $n\geqslant 1$, $H\in]0,1[$
and $\gamma^{2}<\min(1/4,H/(n-1))$,

(3.3)
$\displaystyle\mathbb{E}\left(\delta_{\tau}X_{H,\gamma}\right)^{2n}\mathrel{\mathop{\kern
0.0pt\sim}\limits_{\tau\to
0^{+}}}c_{H,\gamma,2n}\frac{(2n)!}{2^{n}n!}\,\left(\frac{\tau}{T}\right)^{2nH-2n(n-1)\gamma^{2}},$

where the multiplicative constant $c_{H,\gamma,2n}$ reads

(3.4) $\displaystyle c_{H,\gamma,2n}$
$\displaystyle=T^{2nH}e^{2n(n-1)\gamma^{2}g(0)}\times$
$\displaystyle\int_{\\{u_{i}\\}_{1\leqslant i\leqslant
n}\in\mathbb{R}^{n}}\prod_{i=1}^{n}\left[(1-u_{i})^{H-\frac{1}{2}}1_{u_{i}\leqslant
1}-(-u_{i})^{H-\frac{1}{2}}1_{u_{i}\leqslant
0}\right]^{2}\prod_{i<j=1}^{n}\frac{1}{|u_{i}-u_{j}|^{4\gamma^{2}}}\prod_{i=1}^{n}du_{i},$

with $g(0)$ given in Eq. 2.19.

######  Remark 3.2.

Because of the independence of the multiplicative chaos $M_{\epsilon,\gamma}$
and the Wiener process $W$ entering in the definition of the random measure
$dW_{\epsilon,H,\gamma}$ (Eq. 2.2), the MfOU process $X_{H,\gamma}$ has a
symmetric probability distribution. In consequence, all odd moments of
$X_{H,\gamma}$ and its increments vanish.

######  Remark 3.3.

Let us examine the values of the factors $c_{H,\gamma,2n}$ in simple cases.
First, their expression (Eq. 3.4) is understood when $n=1$ as being

(3.5) $\displaystyle c_{H,\gamma,2}=c_{H,0,2}$
$\displaystyle=T^{2H}\int_{u\in\mathbb{R}}\left[(1-u)^{H-\frac{1}{2}}1_{u\leqslant
1}-(-u)^{H-\frac{1}{2}}1_{u\leqslant 0}\right]^{2}du$
$\displaystyle=\frac{T^{2H}}{\sin(\pi
H)}\frac{\left[\Gamma\left(H+\frac{1}{2}\right)\right]^{2}}{\Gamma\left(2H+1\right)},$

which is independent of $\gamma$, and coincides with the one of the fOU
process (its statistical properties are recalled in Remark 2.2), as it is
claimed in Proposition 3.1. Also, it is straightforward to get these factors
for the non-intermittent $\gamma=0$ case since the MfOU process corresponds to
the fOU process, and is thus Gaussian. It is then clear, looking at the
expression provided in Eq. 3.4, that

(3.6) $\displaystyle c_{H,0,2n}=\left(c_{H,0,2}\right)^{n},$

as it is expected for a Gaussian process. The last situation for which the
expression of these factors $c_{H,\gamma,2n}$ (Eq. 3.4) gets simple is when
$H=1/2$. In this case, we obtain

(3.7) $\displaystyle c_{1/2,\gamma,2n}$
$\displaystyle=T^{n}e^{2n(n-1)\gamma^{2}g(0)}\int_{\\{u_{i}\\}_{1\leqslant
i\leqslant
n}\in[0,1]^{n}}\prod_{i<j=1}^{n}\frac{1}{|u_{i}-u_{j}|^{4\gamma^{2}}}\prod_{i=1}^{n}du_{i},$

that can be explicitly computed for the fourth-order moment (i.e. $n=2$) as

(3.8) $\displaystyle c_{1/2,\gamma,4}$
$\displaystyle=T^{2}e^{4\gamma^{2}g(0)}\int_{\\{u_{1},u_{2}\\}\in[0,1]^{2}}\frac{1}{|u_{1}-u_{2}|^{4\gamma^{2}}}du_{1}du_{2}$
$\displaystyle=T^{2}e^{4\gamma^{2}g(0)}\frac{1}{(1-4\gamma^{2})(1-2\gamma^{2})},$

as it was already derived in Ref. [32]. Exact values of $c_{H,\gamma,2n}$ for
any appropriate $H$ and $\gamma$ are not known, but can be evaluated while
performing a numerical estimation of the integrals entering in Eq. 3.4.

######  Remark 3.4.

Notice that the characterization of the sample path properties of
$X_{H,\gamma}$ and its multifractal nature, as it is announced in Proposition
3.1, using only its even moments of integer orders, and those of its
increments, is rather simple and incomplete. A more complete analysis of the
statistical properties of the MfOU process, in the spirit of Ref. [11],
remains to be done.

##  4\. Numerical synthesis and statistical analysis

###  4.1. Numerical synthesis

The purpose of this Section is to propose a simple and efficient numerical
algorithm that generates a trajectory of the statistically stationary MfOU
process $X_{\epsilon,H,\gamma}(t)$ (Eq. 3.1), using similar ideas as they are
developed in [32]. As it is considered in the theoretical framework developed
in Proposition 3.1, trajectories are given for a finite $\epsilon>0$, and for
$H\in]0,1[$ and a certain value of the parameter $\gamma$, the values of $H$
and $\gamma$ being possibly conditioned on the existence of moments of a
certain order as $\epsilon\to 0$. Alternatively, we could solve the underlying
dynamics itself (Eq. 2.1) and thus access the transient towards the
statistically regime using a time-marching discretization, as it is proposed
in [10, 26]. In any case, we need to come up with a numerical approximation of
the underlying random measure $dW_{\epsilon,H,\gamma}(t)$ (Eq. 2.2) which can
be seen as a convolution of some random measure, that we will focus on in the
sequel, with a causal and deterministic kernel, making the overall dynamics
(Eq. 2.1) non-Markovian. It is then tempting to use the Discrete Fourier
transform (DFT) to design a numerical approximation
$\widehat{dW}_{\epsilon,H,\gamma}[t]$ of the continuous measure
$dW_{\epsilon,H,\gamma}(t)$, and incidentally, work in a framework for
$T_{\text{\tiny{tot}}}$-periodic functions.

The causal and deterministic kernel involved in the definition of
$dW_{\epsilon,H,\gamma}(t)$ reads

(4.1)
$h_{\epsilon,H}(t)=\mathfrak{h}_{\epsilon,H}(t)+\epsilon^{H-\frac{1}{2}}\delta(t),$

where

(4.2)
$\mathfrak{h}_{\epsilon,H}(t)=\left(H-\frac{1}{2}\right)(t+\epsilon)^{H-\frac{3}{2}}1_{t\geqslant
0}.$

In this continuous setting, the function $h_{\epsilon,H}$ (Eq. 4.1) behaves
differently depending of the value of the Hurst exponent $H$, as it is
classically interpreted when considering fractional Gaussian noises [22]. For
$H<1/2$, the function $h_{\epsilon,H}$ is negative outside the origin, and is
integrable over $t\in\mathbb{R}$, such that its integral vanishes. On the
opposite, when $H>1/2$, its integral diverges, whereas it is equal to unity
for $H=1/2$. Notice also that, in the vicinity of the origin, $h_{\epsilon,H}$
is distributional, with a weight that depends on $\epsilon$.

From a numerical point of view, over a finite set
$t\in[-T_{\text{\tiny{tot}}}/2,T_{\text{\tiny{tot}}}/2]$, we need to make a
compromise between enforcing integral properties of $h_{\epsilon,H}$ (Eq.
4.1), especially for $H<1/2$, and respecting the strength of its
distributional nature in the vicinity of the origin. The most natural way to
build a numerical estimation $\widehat{h}_{\epsilon,H}$ of the function
$h_{\epsilon,H}$ (Eq. 4.1) is obtained in the Fourier space, such that

(4.3)
$\displaystyle\widehat{h}_{\epsilon,H}[t]=\text{DFT}^{-1}\left\\{\text{DFT}\\{\mathfrak{h}_{\epsilon,H}\\}[\omega]+\epsilon^{H-\frac{1}{2}}\right\\}[t],$

where the function $\mathfrak{h}_{\epsilon,H}$ is defined in Eq. 4.2, properly
periodized over $[-T_{\text{\tiny{tot}}}/2,T_{\text{\tiny{tot}}}/2]$, and the
frequency $\omega$ belonging to the discrete set
$[-1/(2T_{\text{\tiny{tot}}}),1/(2T_{\text{\tiny{tot}}})]$. It is clear that
integral properties of $h_{\epsilon,H}$ (Eq. 4.1) are not exactly fulfilled by
the proposed estimation $\widehat{h}_{\epsilon,H}$. In particular, the
integral over $[-T_{\text{\tiny{tot}}}/2,T_{\text{\tiny{tot}}}/2]$ is equal,
up to numerical errors regarding the finiteness of the time step, to
$(T_{\text{\tiny{tot}}}/2+\epsilon)^{H-1/2}$, which eventually goes towards
the expected limits as $T_{\text{\tiny{tot}}}\to\infty$. In the same time, the
estimation reproduces accurately the distributional nature of the function
$h_{\epsilon,H}$ close to the origin.

An important numerical observation made in Ref. [32] shows that constraining
the estimation $\widehat{h}_{\epsilon,H}$ to fulfill the integral properties
of $h_{\epsilon,H}$ is crucial for the development of logarithmic correlations
for the particular case $H=0$ (Eq. 2.17), and brings numerical stability.
Recall that for this boundary case $H=0$, very different behaviors are
expected. In particular, whereas the variance of the fOU process (Eq. 2.5),
when $H\in]0,1[$, remains bounded as $\epsilon$ goes to zero (Eq. 2.6), it is
expected to diverge when $H=0$ (Eq. 2.11) in a logarithmic fashion (Eq. 2.13).
Because of this observation and for these reasons, we treat the $H=0$ case
separately and propose the following estimator $\widehat{h}_{\epsilon,0}$ of
the function $h_{\epsilon,0}$ as

(4.4)
$\displaystyle\widehat{h}_{\epsilon,0}[t]=\text{DFT}^{-1}\left\\{\text{DFT}\\{\mathfrak{h}_{\epsilon,0}\\}[\omega]-\text{DFT}\\{\mathfrak{h}_{\epsilon,0}\\}[0]\right\\}[t].$

It is now clear that the vanishing integral properties of $h_{\epsilon,0}$ is
exactly fulfilled by its estimation $\widehat{h}_{\epsilon,0}$ (Eq. 4.4),
whereas an unavoidable error is made on the strength of the Dirac distribution
near the origin.

Given these precisions, call $N$ the number of collocation points and $\Delta
t=T_{\text{\tiny{tot}}}/N$ the time step, we propose the following numerical
estimation $\widehat{\widetilde{X}}_{\epsilon}[t]$ of the regularized fOU
process of vanishing Hurst exponent $\widetilde{X}_{\epsilon}(t)$ (Eq. 2.11)
as

(4.5) $\displaystyle\widehat{\widetilde{X}}_{\epsilon}[t]=\Delta
t\;\text{DFT}^{-1}\left\\{\text{DFT}\left\\{e^{-t/T}1_{t\geqslant
0}\right\\}\text{DFT}\left\\{\widehat{h}_{\epsilon,0}\right\\}[\omega]\text{DFT}\left\\{\widehat{d\widetilde{W}}\right\\}\right\\}[t],$

where $\widehat{h}_{\epsilon,0}$ is given in Eq. 4.4, and
$\widehat{d\widetilde{W}}$ a collection of $N$ independent instances of a zero
average Gaussian random variable of variance $\Delta t$. To prevent from
additional aliasing errors, consider long trajectories such that
$T_{\text{\tiny{tot}}}\gg T$, where $T$ is the characteristic large scale of
the OU kernel entering in Eq. 4.5. Accordingly, an estimation
$\widehat{M}_{\epsilon,\gamma}[t]$ of the random positive weight
$M_{\epsilon,\gamma}(t)$ (Eq. 2.4) is obtained while exponentiating the
estimation $\widehat{\widetilde{X}}_{\epsilon}$ (Eq. 4.5), such that

(4.6)
$\widehat{M}_{\epsilon,\gamma}[t]=e^{\gamma\widehat{\widetilde{X}}_{\epsilon}[t]-\gamma^{2}\widehat{\mathbb{E}}\left(\left.\widehat{\widetilde{X}}_{\epsilon}\right.^{2}\right)},$

where enter the empirical expectation $\widehat{\mathbb{E}}$ of the square of
the discrete process $\widehat{\widetilde{X}}_{\epsilon}$. We could have used
alternatively the exact expression of this expectation, that can be easily
derived even for a given finite $\epsilon>0$ (see Eq. A13 of Ref. [32]), which
eventually gives undistinguishable results using the parameters of subsequent
simulations.

Finally, the estimator $\widehat{X}_{\epsilon,H,\gamma}[t]$ of the
statistically stationary MfOU process $X_{\epsilon,H,\gamma}(t)$ (Eq. 3.1)
reads

(4.7) $\displaystyle\widehat{X}_{\epsilon,H,\gamma}[t]=\Delta
t\;\text{DFT}^{-1}\left\\{\text{DFT}\left\\{e^{-t/T}1_{t\geqslant
0}\right\\}\text{DFT}\left\\{\widehat{h}_{\epsilon,H}\right\\}[\omega]\text{DFT}\left\\{\widehat{M}_{\epsilon,\gamma}\widehat{dW}\right\\}\right\\}[t],$

where $\widehat{h}_{\epsilon,H}$ is given in Eq. 4.3,
$\widehat{M}_{\epsilon,\gamma}$ in Eq. 4.6 and $\widehat{dW}$ a collection of
$N$ independent instances of a zero average Gaussian random variable of
variance $\Delta t$, independent of those of $\widehat{d\widetilde{W}}$.

###  4.2. Statistical analysis

We generate numerical instances of the discrete random process
$\widehat{X}_{\epsilon,H,\gamma}$ (Eq. 4.7) taking without loss of generality
$T_{\text{\tiny{tot}}}=1$. We consider $N=2^{30}$, and three different values
for $H$, namely $H=1/3,\;1/2$ and $2/3$, and three different values for
$\gamma$, such that $\gamma^{2}=0,\;0.02$ and $0.04$. For each couple of
values $(H,\gamma^{2})$, we attribute the values
$T=T_{\text{\tiny{tot}}}/2^{10}$ for the characteristic large scale of the OU
kernel entering in Eq. 4.7, and $\epsilon=4\Delta t$ concerning the
regularizing scale, and generate ten independent trajectories.

Figure 4.1. An instance of a MfOU trajectory $\widehat{X}_{\epsilon,H,\gamma}$
(Eq. 4.7), for increasing values of the Hurst parameter $H$ (from left to
right: $H=1/3$, $1/2$ and $2/3$), and increasing values of the intermittent
parameter (from top to bottom, $\gamma^{2}=0$ and $\gamma^{2}=0.04$). See the
text for values of additional numerical parameters.

Figure 4.1 shows an instance of the random process obtained
$\widehat{X}_{\epsilon,H,\gamma}$ (Eq. 4.7) over several time scales $T$ and
for different values of the couple $(H,\gamma)$, with the same instance of the
underlying independent white noises $\widehat{dW}$ and
$\widehat{d\widetilde{W}}$. As $H$ increases, from left to right, the signal
becomes less and less rough. Going from the Gaussian, i.e. non-intermittent,
situation using $\gamma^{2}=0$ (top line), to its multifractal version using
$\gamma^{2}=0.04$ (bottom line), we can see the appearance of spotty
fluctuations, associated to events of larger time variation, which are
reminiscent of the fluctuations of the random weight
$\widehat{M}_{\epsilon,\gamma}$ (Eq. 4.6).

To investigate and quantify further the multiscale nature of the MfOU process
$X_{\epsilon,H,\gamma}$, we present a statistical analysis of some moments of
its increments over a given scale $\tau$, and compare to the predictions made
in Proposition 3.1 of Section 3. More sophisticated analysis could be invoked
to quantify the expected Multifractal phenomenon [33, 2, 17], although we
choose to keep the analysis simple and focus on the power-law behavior with
$\tau$ of these moments. To do so, we compute the empirical moments of the
increments
$\delta_{\tau}\widehat{X}_{\epsilon,H,\gamma}[t]=\widehat{X}_{\epsilon,H,\gamma}[t+\tau]-\widehat{X}_{\epsilon,H,\gamma}[t]$
of our estimation $\widehat{X}_{\epsilon,H,\gamma}$ (Eq. 4.7), and average
over ten independent trajectories. In particular, we display in Fig. 4.2 the
results of our statistical analysis using the so-called second order structure
function $S_{2}(\tau)$, obtained taking $n=2$ in the expression

(4.8)
$S_{n}(\tau)\equiv\widehat{\mathbb{E}}\left[\left(\delta_{\tau}\widehat{X}_{\epsilon,H,\gamma}\right)^{n}\right]\equiv\frac{\Delta
t}{T_{\text{\tiny{tot}}}-\tau}\sum_{i=1}^{N-\tau/\Delta
t}\left(\delta_{\tau}\widehat{X}_{\epsilon,H,\gamma}[t_{i}]\right)^{n},$

and the respective renormalized flatness factor

(4.9) $F(\tau)=\frac{S_{4}(\tau)}{3S_{2}^{2}(\tau)},$

which coincides with unity for Gaussian processes.

Figure 4.2. Second-order moment $S_{2}(\tau)$ (Eq. 4.8) (top line) and
flatness $F(\tau)$ (Eq. 4.9) (bottom line) of the increments of the MFOU
process for different values of $H$ (from left to right $H=1/3$, $1/2$ and
$2/3$) and different values of $\gamma^{2}$ ($\gamma^{2}=0$ ($\star$),
$\gamma^{2}=0.02$ ($\triangledown$) and $\gamma^{2}=0.04$ ($\square$), in a
logarithmic representation. The red lines correspond to the theoretical
predictions concerning the asymptotic regime $\epsilon\to 0$ and given in Eq.
3.3. We have superimposed vertical dashed lines to show the characteristic
scales $\epsilon$ and $T$ used for the simulation.

We display in Fig. 4.2 (top) in a logarithmic fashion our statistical
estimation of the second-order structure function $S_{2}(\tau)$ (Eq. 4.8),
after having checked that all represented scales have converged from a
statistical point of view. As it is expected, and stated in Proposition 3.1,
second-order statistical properties are not impacted by multifractal
corrections, and are thus independent of the parameter $\gamma$: different
symbols of Fig. 4.2, representing different values for $\gamma$, are gathered
on a same curve. We superimpose, using a red solid line the power-law
predicted in Proposition 3.1 and provided in Eqs. 3.3 and 3.4, as far as
second-order $n=2$ is concerned (see Eq. 3.5 for an explicit expression at
second-order). We can see that the prediction reproduces accurately the decay
of the variance of increments as the scale $\tau$ decreases, for both scaling
and amplitude (no additional fitting parameter is involved), in the inertial
range $\epsilon\ll\tau\ll T$, as it is called in the literature of fluid
turbulence [13]. For scales below $\epsilon$, we recover the scaling of the
underlying OU process, that predicts $S_{2}(\tau)$ proportional to $\tau$.
Above $T$, $S_{2}(\tau)$ saturates towards two times the variance of the
process (Eq. 2.6), as expected from a statistically stationary process.
Actually, the whole range of available scales can be predicted and
successfully compared to our estimations using the expression provided in Eq.
2.8, properly generalized to a finite $\epsilon>0$, as it can be derived using
the material developed in Refs. [10, 32].

Let us now focus on the increments flatness $F(\tau)$ (Eq. 4.9) that we
display in Fig. 4.2 (bottom). We can observe, as it is expected from a
Gaussian process, $F(\tau)$ gets independent on scale $\tau$ when $\gamma=0$
(represented with $\star$), coinciding with $F(\tau)=1$, consistently with the
predictions listed in Proposition 3.1. When $\gamma$ differs from zero (using
$\triangledown$ for $\gamma^{2}=0.02$ and $\square$ for $\gamma^{2}=0.04$), we
observe a power-law behavior in the inertial range. Amplitude and exponent of
this power-law are fully captured by our asymptotic predictions (Eq. 3.3),
obtained in the limit $\epsilon\to 0$ and for small scales, that read

(4.10) $F(\tau)\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0^{+}}}\frac{c_{H,\gamma,4}}{c_{H,\gamma,2}^{2}}\left(\frac{\tau}{T}\right)^{-4\gamma^{2}},$

where the coefficients $c_{H,\gamma,2n}$ are provided in Eq. 3.4. Whereas we
can get an explicit expression when $H=1/2$, as it is recalled in Eq. 3.8, we
can get their numerical value performing a double numerical integration. Doing
so, we observe that the coefficient in front of the expected power-law (Eq.
4.10) depends weakly on $H$ (data not shown), and remains close to the exact
value obtained when $H=1/2$ (Eq. 3.8) for the range of values that are
explored in Fig. 4.2. We superimpose these predictions using red lines, and
observe that they collapse with our statistical estimation of the flatness in
the inertial range $\epsilon\ll\tau\ll T$. For scales below $\tau$, the
flatness saturates towards a plateau that depends in particular on the precise
value of $\epsilon$. For scales above $T$, flatness also saturates, towards a
value that depends on chosen parameters. Overall, we can compare with great
success our asymptotic theoretical predictions (Eq. 3.3) with a statistical
analysis of simulated trajectories.

Figure 4.3. Probability density function of increments, estimated while
computing histograms, for different scales $\tau$. Shown scales $\tau$ are
logarithmically spaced between $\epsilon/4$ and $4T$. For each scale, three
curves are displayed corresponding to different values of the Hurst parameter
$H$ ($H=1/3,1/2$ and $2/3$), and are totally indistinguishable.

Finally, for the sake of completeness, we represent in Fig. 4.3 the
Probability Density Functions (PDFs) of increments, estimated as histograms,
for different scales $\tau$ and for different values of $H$ and $\gamma$, as
they are considered in Figs. 4.1 and 4.2. For clarity, curves at different
scales $\tau$ have been vertically shifted by an arbitraly factor, and are
rescaled such that they are all of unit-variance. Firstly, we can observe that
their shape depends indeed very weakly on the value of $H$, as we noticed
while analyzing the values of the flatnesses (Fig. 4.2): we can barely
distinguish them as $H$ varies. Also, we see that they remain Gaussian when
$\gamma=0$ (left colomn), whereas they develop longer and longer tails as the
scale decreases when $\gamma\neq 0$. Furthermore, all PDFs are symmetrical, as
it is expected and recalled in Proposition 3.1.

##  5\. Proofs

Proving the statistical properties of the MfOU process $X_{H,\gamma}$, as they
are stated in Proposition 3.1, boils down to computing the moments of a random
variable defined by a linear operation on the random measure
$dW_{\epsilon,H,\gamma}$ (Eq. 2.2), and then studying their limit when
$\epsilon\to 0$. That is why we consider in the sequel the random variable

(5.1) $Y_{\epsilon,H,\gamma}=\int_{\mathbb{R}}f(s)dW_{\epsilon,H,\gamma}(s),$

where $f$ is an appropriate deterministic test function, say bounded and
compactly supported for the sake of presentation. For the time being, assume
furthermore that $f$ is differentiable over $\mathbb{R}$, note by $f^{\prime}$
its derivative, an assumption that we will revisit and adapt when we will
treat the particular case of the MfOU process $X_{H,\gamma}$ in Section 5.3.
We recall here for convenience the expression of the random measure
$dW_{\epsilon,H,\gamma}$ that enters in the definition of the random variable
$Y_{\epsilon,H,\gamma}$ (Eq. 5.1), as it was done in Definition 2.1. We have

(5.2)
$dW_{\epsilon,H,\gamma}(s)=\omega_{\epsilon,H,\gamma}(s)ds+\epsilon^{H-\frac{1}{2}}M_{\epsilon,\gamma}(s)dW(s),$

where

(5.3)
$\omega_{\epsilon,H,\gamma}(s)=\left(H-\frac{1}{2}\right)\int_{-\infty}^{s}(s-u+\epsilon)^{H-\frac{3}{2}}M_{\epsilon,\gamma}(u)dW(u),$

and $M_{\epsilon,\gamma}$ a regularized version of a Multiplicative Chaos (see
Remark 2.4), normalized such that $\mathbb{E}M^{2}_{\epsilon,\gamma}=1$.
Recall also that $M_{\epsilon,\gamma}$ (Eq. 2.4) is chosen independently of
$dW$. It is then straightforward to see that $Y_{\epsilon,H,\gamma}$ is a
zero-average random variable, and has more generally vanishing odd moments. In
the sequel, we will find convenient to express the random measure
$dW_{\epsilon,H,\gamma}$ (Eq. 5.2) as a convolution product of the random
measure $M_{\epsilon,\gamma}dW$, i.e.

(5.4)
$dW_{\epsilon,H,\gamma}(s)=ds\int_{\mathbb{R}}h_{\epsilon,H}(s-u)M_{\epsilon,\gamma}(u)dW(u),$

with the deterministic kernel

(5.5)
$h_{\epsilon,H}(s)=\left(H-\frac{1}{2}\right)(s+\epsilon)^{H-\frac{3}{2}}1_{s\geqslant
0}+\epsilon^{H-\frac{1}{2}}\delta(s),$

as it was already considered in Section 4 devoted to numerics (see Eq. 4.1).

###  5.1. Second-order statistical properties

To prepare for subsequent calculations, first consider second-order
statistical properties, and compute the variance of $Y_{\epsilon,H,\gamma}$,
in a similar way, although more general, as it has been done in Ref. [10]. We
have

(5.6) $\displaystyle\mathbb{E}\left[Y_{\epsilon,H,\gamma}^{2}\right]$
$\displaystyle=\int_{\mathbb{R}^{2}}f(s_{1})f(s_{2})\mathbb{E}\left[dW_{\epsilon,H,\gamma}(s_{1})dW_{\epsilon,H,\gamma}(s_{2})\right]$
$\displaystyle=\int_{\mathbb{R}^{3}}f(s_{1})f(s_{2})h_{\epsilon,H}(s_{1}-u)h_{\epsilon,H}(s_{2}-u)\mathbb{E}\left[M^{2}_{\epsilon,\gamma}(u)\right]ds_{1}ds_{2}du$
$\displaystyle=\int_{\mathbb{R}}\mathcal{L}_{\epsilon,H}(u)du,$

where we have used that $M_{\epsilon,\gamma}$ is independent of $dW$ and
normalized such that $\mathbb{E}(M_{\epsilon,\gamma}^{2})=1$ (see Remark 2.4),
and introduced the quantity

(5.7)
$\displaystyle\mathcal{L}_{\epsilon,H}(u)=\int_{\mathbb{R}^{2}}f(s_{1})f(s_{2})h_{\epsilon,H}(s_{1}-u)h_{\epsilon,H}(s_{1}-u)ds_{1}ds_{2}.$

This shows that there are no multifractal corrections on the second-order
statistical properties of the random variable $Y_{\epsilon,H,\gamma}$, which
coincide with those of $Y_{\epsilon,H,0}$. More explicitly, the deterministic
function $\mathcal{L}_{\epsilon,H}$ (Eq. 5.7) reads

(5.8)
$\displaystyle\mathcal{L}_{\epsilon,H}(u)=\left(H-\frac{1}{2}\right)^{2}\int_{\mathbb{R}^{2}}f(s_{1})f(s_{2})(s_{1}-u+\epsilon)^{H-3/2}(s_{2}-u+\epsilon)^{H-3/2}1_{s_{1}\geqslant
u}1_{s_{2}\geqslant u}ds_{1}ds_{2}$
$\displaystyle+2\epsilon^{H-1/2}\left(H-\frac{1}{2}\right)\int_{\mathbb{R}}f(s_{1})f(u)(s_{1}-u+\epsilon)^{H-3/2}1_{s_{1}\geqslant
u}ds_{1}+\epsilon^{2H-1}f^{2}(u).$

For $H\in]1/2,1[$, we can see that, in the limit $\epsilon\to 0$, all
integrals entering in $\mathcal{L}_{\epsilon,H}$ (Eq. 5.8) involve singular
kernels that are all locally integrable, thus only the first term of the RHS
survives. Call this limit the function $\mathcal{L}_{H}(u)$. It is a
continuous and bounded function, which is furthermore compactly supported if
it is assumed so for the test function $f$. In this case, we can write

(5.9) $\displaystyle\mathbb{E}\left[Y_{H>1/2,\gamma}^{2}\right]$
$\displaystyle\equiv\lim_{\epsilon\to
0}\mathbb{E}\left[Y_{\epsilon,H,\gamma}^{2}\right]\equiv\int_{\mathbb{R}}\mathcal{L}_{H}(u)du$
$\displaystyle=\left(H-\frac{1}{2}\right)^{2}\int_{\mathbb{R}}\left(\int_{\mathbb{R}}f(s)(s-u)^{H-3/2}1_{s\geqslant
u}ds\right)^{2}du.$

This last expression could be even further simplified, as it was done in Ref.
[10], while introducing the correlation product, although, as we will see,
this procedure is difficult to generalize for higher order moments. Notice
that this expression makes perfect sense since all singularities are locally
integrable, and $f$ is compactly supported.

For $H=1/2$, the situation is even simpler since the process
$\omega_{\epsilon,H,\gamma}$ (Eq. 5.3) vanishes, and we obtain

(5.10)
$\displaystyle\mathbb{E}\left[Y_{H=1/2,\gamma}^{2}\right]=\int_{\mathbb{R}}f^{2}(s)ds.$

Take now $H\in]0,1/2[$. Going back to the expression of
$\mathcal{L}_{\epsilon,H}$ (Eq. 5.8), we can see that the first term in the
RHS has no pointwise limit when $\epsilon\to 0$ since singularities are non
locally integrable, neither do the two other terms. To extract finite and
diverging contributions, split the integral of the first term of the RHS of
Eq. 5.8 in two parts while considering the integration over the sets
$s_{1}\leqslant s_{2}$ and $s_{1}\geqslant s_{2}$, and notice that they give
an equal contribution. Perform then an integration by parts over the dummy
variable $s_{1}$ and obtain

(5.11) $\displaystyle
2\left(H-\frac{1}{2}\right)^{2}\int_{\mathbb{R}^{2}}f(s_{1})f(s_{2})(s_{1}-u+\epsilon)^{H-3/2}(s_{2}-u+\epsilon)^{H-3/2}1_{s_{1}\geqslant
u}1_{s_{2}\geqslant s_{1}}ds_{1}ds_{2}$
$\displaystyle=\mathcal{I}_{\epsilon,H}(u)+\mathcal{J}_{\epsilon,H}(u),$

where we have introduced the quantities

$\displaystyle\mathcal{I}_{\epsilon,H}(u)=\mathcal{I}_{1,H,\epsilon}(u)-\mathcal{I}_{2,H,\epsilon}(u)$

with

(5.12) $\displaystyle\mathcal{I}_{1,H,\epsilon}(u)$
$\displaystyle=2\left(H-\frac{1}{2}\right)\int
f^{2}(s)(s-u+\epsilon)^{2H-2}1_{u\leqslant s}ds$
$\displaystyle=-\epsilon^{2H-1}f^{2}(u)-2\int
f(s)f^{\prime}(s)(s-u+\epsilon)^{2H-1}1_{u\leqslant s}ds,$ (5.13)
$\displaystyle\mathcal{I}_{2,H,\epsilon}(u)$
$\displaystyle=2\left(H-\frac{1}{2}\right)\epsilon^{H-\frac{1}{2}}\int
f(u)f(s)(s-u+\epsilon)^{H-\frac{3}{2}}1_{u\leqslant s}ds,$

and

(5.14)
$\displaystyle\mathcal{J}_{\epsilon,H}(u)=-2\left(H-\frac{1}{2}\right)\int
f^{\prime}(s_{1})f(s_{2})(s_{1}-u+\epsilon)^{H-\frac{1}{2}}(s_{2}-u+\epsilon)^{H-\frac{3}{2}}1_{u\leqslant
s_{1}}1_{s_{1}\leqslant s_{2}}ds_{1}ds_{2}.$

Adding everything up, noticing that $\mathcal{I}_{2,H,\epsilon}(u)$ (Eq. 5.13)
will cancel with the second term of the RHS of Eq. 5.8, and that the
contribution of order $\epsilon^{2H-1}$ of $\mathcal{I}_{1,H,\epsilon}(u)$
(Eq. 5.12) will cancel with the third one, we obtain

(5.15)
$\displaystyle\mathcal{L}_{\epsilon,H}(u)=\mathcal{J}_{\epsilon,H}(u)+\mathcal{K}_{\epsilon,H}(u),$

where $\mathcal{J}_{\epsilon,H}$ is defined in Eq. 5.14 and

(5.16) $\displaystyle\mathcal{K}_{\epsilon,H}(u)$ $\displaystyle=-2\int
f(s)f^{\prime}(s)(s-u+\epsilon)^{2H-1}1_{u\leqslant s}ds$
$\displaystyle=-2\int f(s+u)f^{\prime}(s+u)(s+\epsilon)^{2H-1}1_{s\geqslant
0}ds.$

Whereas the pointwise limit as $\epsilon\to 0$ of the function
$\mathcal{L}_{\epsilon,H}$ as written in Eq. 5.8 is not obvious when $H<1/2$,
it becomes clear once written as in Eq. 5.15. To see this, remark that the
function $\mathcal{K}_{\epsilon,H}$ (Eq. 5.16) is continuous, bounded and its
integral over $u\in\mathbb{R}$ vanishes. Note its pointwise limit as the
function $\mathcal{K}_{H}(u)$. Similarly, while taking the limit $\epsilon\to
0$, the strongest singularity entering in the expression of
$\mathcal{J}_{\epsilon,H}$ (Eq. 5.14) is obtained along the diagonal
$s_{1}=s_{2}=s=u$ and is of order $(s-u)^{2H-2}$, which is locally integrable
in $\mathbb{R}^{2}$ as long as $H>0$. Call then $\mathcal{J}_{H}(u)$ the
pointwise limit of $\mathcal{J}_{\epsilon,H}(u)$.

We finally obtain, after taking the limit $\epsilon\to 0$,

(5.17)
$\displaystyle\mathbb{E}\left[Y_{H<1/2,\gamma}^{2}\right]=\int_{\mathbb{R}}\mathcal{L}_{H}(u)du=\int_{\mathbb{R}}\mathcal{J}_{H}(u)+\mathcal{K}_{H}(u)du=\int_{\mathbb{R}}\mathcal{J}_{H}(u)du$
$\displaystyle=-2\left(H-\frac{1}{2}\right)\int_{\mathbb{R}^{3}}f^{\prime}(s_{1})f(s_{2})(s_{1}-u)^{H-1/2}(s_{2}-u)^{H-3/2}1_{u\leqslant
s_{1}}1_{s_{1}\leqslant s_{2}}duds_{1}ds_{2}.$

###  5.2. Higher-order moments

Let us now develop and generalize our calculations, as they were done to
derive the variance of the random variable $Y_{\epsilon,H,\gamma}$ (Eq. 5.1),
for the moment of order $2n$, with $n$ and integer greater than unity. As
mentioned in Proposition 3.1, it is clear that its probability law is
symmetrical, as a direct consequence of the independence of the white noise
field $dW(t)_{t\in\mathbb{R}}$ and the random positive weight
$M_{\epsilon,\gamma}(t)_{t\in\mathbb{R}}$ (Eq. 2.4). Thus all moments of odd
order vanish.

Start with introducing the following quantity that encodes the correlation
structure of the multiplicative chaos $M_{\epsilon,\gamma}$ (Eq. 2.4) and
write

(5.18) $\displaystyle\mathcal{M}_{\epsilon,\gamma}(\\{u_{i}\\}_{1\leqslant
i\leqslant
n})=\mathbb{E}\left[\prod_{i=1}^{n}M^{2}_{\epsilon,\gamma}(u_{i})\right],$

which converges as $\epsilon\to 0$ towards a locally integrable function in
$\mathbb{R}^{n}$ when $\gamma^{2}<\min[1/4,1/(2(n-1))]$ (see Remark 2.4), and
that reads

(5.19) $\displaystyle\mathcal{M}_{\gamma}(\\{u_{i}\\})\equiv\lim_{\epsilon\to
0}\mathcal{M}_{\epsilon,\gamma}(\\{u_{i}\\})=T^{2n(n-1)\gamma^{2}}\prod_{i<j=1}^{n}\frac{e^{4\gamma^{2}g(u_{i}-u_{j})}}{|u_{i}-u_{j}|_{+}^{4\gamma^{2}}}.$

We have

(5.20)
$\displaystyle\mathbb{E}\left(Y_{\epsilon,H\in]0,1[,\gamma}^{2n}\right)=\int_{\\{s_{i}\\}\in\mathbb{R}^{2n}}\prod_{i=1}^{2n}f(s_{i})\mathbb{E}\left[\prod_{i=1}^{2n}dW_{\epsilon,H,\gamma}(s_{i})\right]$
$\displaystyle=\int_{\\{s_{i},u_{i}\\}\in\mathbb{R}^{4n}}\prod_{i=1}^{2n}f(s_{i})h_{\epsilon,H}(s_{i}-u_{i})\mathbb{E}\left[\prod_{i=1}^{2n}dW(u_{i})\right]\mathbb{E}\left[\prod_{i=1}^{2n}M_{\epsilon,\gamma}(u_{i})\right]\prod_{i=1}^{2n}ds_{i}du_{i}$
$\displaystyle=\frac{(2n)!}{2^{n}n!}\int_{\\{s_{i},u_{i}\\}\in\mathbb{R}^{3n}}\mathcal{M}_{\epsilon,\gamma}(\\{u_{i}\\})\prod_{i=1}^{n}f(s_{2i-1})f(s_{2i})h_{\epsilon,H}(s_{2i-1}-u_{i})h_{\epsilon,H}(s_{2i}-u_{i})ds_{2i-1}ds_{2i}du_{i}$
$\displaystyle=\frac{(2n)!}{2^{n}n!}\int_{\\{u_{i}\\}\in\mathbb{R}^{n}}\mathcal{M}_{\epsilon,\gamma}(\\{u_{i}\\})\prod_{i=1}^{n}\mathcal{L}_{\epsilon,H}(u_{i})du_{i},$

where the function $\mathcal{L}_{\epsilon,H}(u)$ is defined in Eq. 5.7. Recall
from Section 5.1 that the function $\mathcal{L}_{\epsilon,H}(u)$ converges
pointwise as $\epsilon\to 0$ towards a function $\mathcal{L}_{H}(u)$ which
expression depends on the value of the parameter $H$. In particular, we
obtained

(5.21) $\displaystyle\mathcal{L}_{H}(u)$ $\displaystyle\mathrel{\mathop{\kern
0.0pt=}\limits^{H>1/2}}\left(H-\frac{1}{2}\right)^{2}\left(\int_{\mathbb{R}}f(s)(s-u)^{H-3/2}1_{s\geqslant
u}ds\right)^{2}$ (5.22) $\displaystyle\mathrel{\mathop{\kern
0.0pt=}\limits^{H=1/2}}f^{2}(u)$ (5.23) $\displaystyle\mathrel{\mathop{\kern
0.0pt=}\limits^{H<1/2}}\mathcal{J}_{H}(u)+\mathcal{K}_{H}(u),$

with

(5.24) $\displaystyle\mathcal{J}_{H}(u)$
$\displaystyle=-2\left(H-\frac{1}{2}\right)\int_{\mathbb{R}^{2}}f^{\prime}(s_{1})f(s_{2})(s_{1}-u)^{H-\frac{1}{2}}(s_{2}-u)^{H-\frac{3}{2}}1_{u\leqslant
s_{1}}1_{s_{1}\leqslant s_{2}}ds_{1}ds_{2},$ (5.25)
$\displaystyle\mathcal{K}_{H}(u)$
$\displaystyle=-2\int_{\mathbb{R}}f(s+u)f^{\prime}(s+u)s^{2H-1}1_{s\geqslant
0}ds.$

Remark that for any $H\in]0,1[$, the functions $\mathcal{L}_{H}$ and $f$ share
similar supports, which ensures that a product on them will be integrable at
large arguments. So justifying the existence of the limit

(5.26) $\displaystyle\mathbb{E}\left(Y_{H,\gamma}^{2n}\right)$
$\displaystyle\equiv\lim_{\epsilon\to
0}\mathbb{E}\left(Y_{\epsilon,H,\gamma}^{2n}\right)$
$\displaystyle=\frac{(2n)!}{2^{n}n!}\int_{\\{u_{i}\\}\in\mathbb{R}^{n}}\mathcal{M}_{\gamma}(\\{u_{i}\\})\prod_{i=1}^{n}\mathcal{L}_{H}(u_{i})du_{i},$

where the pointwise limit $\mathcal{M}_{\gamma}$ of
$\mathcal{M}_{\epsilon,\gamma}$ is provided in Eq. 5.19, and in particular
determining the range of parameters $H$ and $\gamma$ for which this limit
makes sense, requires more assumptions on the test function $f$ than only
continuity and boundedness. Instead of treating this question from a general
point of view, we propose in the following Section to focus on the test
function suggested by the MfOU process $X_{H,\gamma}$, and to justify its
statistical properties as they are stated in Proposition 3.1.

###  5.3. Statistical properties of the MfOU process

####  5.3.1. Variance and higher order moments of the process

Let us now go back to our initial problem, that is obtaining the marginals of
the process $X_{\epsilon,H,\gamma}$ (Eq. 3.1), and their limit as the
regularizing scale $\epsilon$ goes to zero. We make use of the statistical
stationarity of this process, and consider expectations at the time $t=0$.
Concerning the moments of the process, consider the particular test function

(5.27) $f(s)=e^{\frac{s}{T}}1_{s\leqslant 0},$

such that indeed $X_{\epsilon,H,\gamma}(t)$ (Eq. 3.1) at the given time $t=0$
can be written as

(5.28)
$X_{\epsilon,H,\gamma}(0)=\int_{\mathbb{R}}f(s)dW_{\epsilon,H,\gamma}(s).$

As we have seen in Section 5.2, for vanishing regularizing scale $\epsilon\to
0$, its moment of order $2n$ (Eq. 5.26) brings into play the function
$\mathcal{L}_{H}(u)$ which explicit expression depends of the value of $H$
(Eqs. 5.21, 5.22 and 5.23). We have

(5.29) $\displaystyle\mathcal{L}_{H}(u)$ $\displaystyle\mathrel{\mathop{\kern
0.0pt=}\limits^{H>1/2}}1_{u\leqslant
0}\left(H-\frac{1}{2}\right)^{2}\left(\int_{s=u}^{0}e^{\frac{s}{T}}(s-u)^{H-3/2}ds\right)^{2}$
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{u\to
0^{-}}}(-u)^{2H-1},$

which shows that $\mathcal{L}_{H}$ is a bounded function that goes to zero at
the origin as fast as a power-law. This behavior in the vicinity of the origin
weakens the singular behavior of the multiplicative chaos, and extends its
range of integrability (Eq. 5.19). Hence, the $2n$-order moment of the MfOU
process is given by, for $H\in]1/2,1[$ and
$\gamma^{2}<\min\left(1/4,H/(n-1)\right)$,

(5.30)
$\displaystyle\mathbb{E}\left(X_{H>1/2,\gamma}^{2n}\right)=\frac{(2n)!}{2^{n}n!}\int_{\\{u_{i}\\}\in\mathbb{R}^{n}}\mathcal{M}_{\gamma}(\\{u_{i}\\})\prod_{i=1}^{n}\mathcal{L}_{H}(u_{i})du_{i},$

where the function $\mathcal{L}_{H}$ is given in Eq. 5.29 and
$\mathcal{M}_{\gamma}$ in Eq. 5.19. The case $H=1/2$ is rather simple since
$\mathcal{L}_{H}(u)=f^{2}(u)=\exp(2u/T)1_{u\leqslant 0}$, and we obtain for
$\gamma^{2}<\min\left(1/4,1/(2(n-1))\right)$

(5.31)
$\displaystyle\mathbb{E}\left(X_{H=1/2,\gamma}^{2n}\right)=\frac{(2n)!}{2^{n}n!}\int_{\\{u_{i}\\}\in\mathbb{R}^{n}}\mathcal{M}_{\gamma}(\\{u_{i}\\})\prod_{i=1}^{n}f^{2}(u_{i})du_{i}.$

Take now $H\in]0,1/2[$, and remark that, using the expression of the test
function $f$ (Eq. 5.27),

(5.32) $\displaystyle f^{\prime}(s)=\frac{1}{T}f(s)-\delta(s),$

which shows that the function $f$ is not differentiable on $\mathbb{R}$ in a
classical sense, as expected from a causal kernel. For this reason, the
expression of $\mathcal{J}_{H}$ (Eq. 5.24) where enters the derivative
$f^{\prime}$ is questionable. Nonetheless, we can argue, going back to the
integration by parts procedure that takes place in the decomposition proposed
in Eq. 5.11, that it makes sense to use the distributional derivative (Eq.
5.32) in the expression of the function $\mathcal{J}_{H}$ (Eq. 5.24) that
enters the expression of $\mathcal{L}_{H}$ (Eq. 5.23). In this case, the Dirac
distribution entering in the expression of $f^{\prime}$ (Eq. 5.32) gives no
contribution, and we obtain

(5.33) $\displaystyle\mathcal{J}_{H}(u)$
$\displaystyle=-\frac{2}{T}\left(H-\frac{1}{2}\right)\int_{\mathbb{R}^{2}}f(s_{1})f(s_{2})(s_{1}-u)^{H-\frac{1}{2}}(s_{2}-u)^{H-\frac{3}{2}}1_{u\leqslant
s_{1}}1_{s_{1}\leqslant s_{2}}ds_{1}ds_{2}$
$\displaystyle=-\frac{2}{T}\left(H-\frac{1}{2}\right)1_{u\leqslant
0}\int_{[u,0]^{2}}e^{\frac{s_{1}+s_{2}}{T}}(s_{1}-u)^{H-\frac{1}{2}}(s_{2}-u)^{H-\frac{3}{2}}1_{s_{1}\leqslant
s_{2}}ds_{1}ds_{2}.$

It is a positive, bounded and integrable function. To see this, we have for
large negative arguments

$\displaystyle\mathcal{J}_{H}(u)\mathrel{\mathop{\kern
0.0pt\sim}\limits_{u\to-\infty}}-\frac{2}{T}\left(H-\frac{1}{2}\right)(-u)^{2H-2}\int_{(\mathbb{R}^{-})^{2}}e^{\frac{s_{1}+s_{2}}{T}}1_{s_{1}\leqslant
s_{2}}ds_{1}ds_{2}=-T\left(H-\frac{1}{2}\right)(-u)^{2H-2},$

which is indeed integrable. To see the behavior of $\mathcal{J}_{H}(u)$ at the
origin, rescale the dummy variables $s_{1}$ and $s_{2}$ of Eq. 5.33 by $-u$
and get

$\displaystyle\mathcal{J}_{H}(u)$ $\displaystyle\mathrel{\mathop{\kern
0.0pt\sim}\limits_{u\to
0^{-}}}-\frac{2}{T}\left(H-\frac{1}{2}\right)(-u)^{2H}\int_{[-1,0]^{2}}(s_{1}+1)^{H-\frac{1}{2}}(s_{2}+1)^{H-\frac{3}{2}}1_{s_{1}\leqslant
s_{2}}ds_{1}ds_{2}$ $\displaystyle\mathrel{\mathop{\kern
0.0pt\sim}\limits_{u\to
0^{-}}}-\frac{1}{T}\frac{H-\frac{1}{2}}{H\left(H+\frac{1}{2}\right)}(-u)^{2H}.$

Whereas the function $\mathcal{J}_{H}(u)$ remains bounded, the function
$\mathcal{K}_{H}(u)$ entering in the expression of $\mathcal{L}_{H}(u)$ (Eq.
5.23) will be singular at the origin, although integrable. To see this, notice
first that we cannot make sense in a simple fashion of the expression of the
function $\mathcal{K}_{H}$ (Eq. 5.25) while injecting the distributional
derivative of $f$ (Eq. 5.32) in it. Instead, go back to the integration by
parts procedure that takes place in the decomposition proposed in Eq. 5.11,
and obtain, using the causal kernel $f$ (Eq. 5.27),

(5.34) $\displaystyle\mathcal{K}_{H}(u)$ $\displaystyle=\lim_{\epsilon\to
0}1_{u\leqslant
0}\left[2\left(H-\frac{1}{2}\right)\int_{u}^{0}e^{2\frac{s}{T}}(s-u+\epsilon)^{2H-2}ds+\epsilon^{2H-1}e^{2\frac{u}{T}}\right]$
$\displaystyle=\lim_{\epsilon\to 0}1_{u\leqslant
0}\left[(-u+\epsilon)^{2H-1}-\frac{2}{T}\int_{u}^{0}e^{2\frac{s}{T}}(s-u+\epsilon)^{2H-1}ds\right]$
$\displaystyle=1_{u\leqslant
0}\left[(-u)^{2H-1}-\frac{2}{T}\int_{u}^{0}e^{2\frac{s}{T}}(s-u)^{2H-1}ds\right].$

This shows that $\mathcal{K}_{H}(u)$ is singular in the vicinity of the
origin, such that

(5.35) $\displaystyle\mathcal{K}_{H}(u)\mathrel{\mathop{\kern
0.0pt\sim}\limits_{u\to 0^{-}}}(-u)^{2H-1},$

since the second term of the third equality of Eq. 5.34 is bounded in $u$, and
actually behaves proportionally to $(-u)^{2H}$ near the origin. To check the
integrability at large arguments, write $\mathcal{K}_{H}(u)$ (Eq. 5.34) as

(5.36) $\displaystyle\mathcal{K}_{H}(u)$ $\displaystyle=1_{u\leqslant
0}\left[\frac{2}{T}\int_{u}^{0}e^{2\frac{s}{T}}\left[(-u)^{2H-1}-(s-u)^{2H-1}\right]ds+e^{2\frac{u}{T}}(-u)^{2H-1}\right]$
$\displaystyle\mathrel{\mathop{\kern
0.0pt\sim}\limits_{u\to-\infty}}-\frac{2}{T}(2H-1)(-u)^{2H-2}\int_{-\infty}^{0}s\;e^{2\frac{s}{T}}ds$
$\displaystyle\mathrel{\mathop{\kern
0.0pt\sim}\limits_{u\to-\infty}}\frac{T}{2}(2H-1)(-u)^{2H-2},$

which is indeed integrable at large arguments. Actually, it can be shown that
its integral over $u\in\mathbb{R}$ vanishes. To see this, before taking the
limit $\epsilon\to 0$ in Eq. 5.34, observe that the integral vanishes for any
$\epsilon>0$, which is also true in the limit.

Hence, the existence of the $2n$-order moment of the MfOU process is governed
by the singular behaviors of $\mathcal{K}_{H}(u)$ (Eq. 5.35) and
$\mathcal{M}_{\gamma}$ (Eq. 5.19), and is given by, for $H\in]0,1/2[$ and
$\gamma^{2}<\min\left(1/4,H/(n-1)\right)$,

(5.37)
$\displaystyle\mathbb{E}\left(X_{H<1/2,\gamma}^{2n}\right)=\frac{(2n)!}{2^{n}n!}\int_{\\{u_{i}\\}\in\mathbb{R}^{n}}\mathcal{M}_{\gamma}(\\{u_{i}\\})\prod_{i=1}^{n}\mathcal{L}_{H}(u_{i})du_{i},$

where the function $\mathcal{L}_{H}$ is given by
$\mathcal{L}_{H}=\mathcal{J}_{H}+\mathcal{K}_{H}$, with $\mathcal{J}_{H}$
given in Eq. 5.33 and $\mathcal{K}_{H}$ in Eq. 5.34.

####  5.3.2. Variance and higher order moments of the increments of the
process

Proceed now with the calculation of the moments of velocity increments, and
determine in particular their behavior at small scales. To do so, consider the
particular test function

(5.38) $\displaystyle f(s)\equiv f_{\tau}(s)$
$\displaystyle=e^{-\frac{\tau-s}{T}}1_{\tau-s\geqslant
0}-e^{-\frac{-s}{T}}1_{-s\geqslant 0}$
$\displaystyle=e^{\frac{s}{T}}\left(e^{-\frac{\tau}{T}}1_{\tau-s\geqslant
0}-1_{-s\geqslant 0}\right).$

such that the increment $\delta_{\tau}X_{\epsilon,H,\gamma}(t)\equiv
X_{\epsilon,H,\gamma}(t+\tau)-X_{\epsilon,H,\gamma}(t)$ of the MfOU process
$X_{\epsilon,H,\gamma}$ (Eq. 3.1) at the given time $t=0$ can be written as

(5.39)
$\delta_{\tau}X_{\epsilon,H,\gamma}(0)=\int_{\mathbb{R}}f(s)dW_{\epsilon,H,\gamma}(s).$

First consider $H>1/2$. In this case, make use of Eq. 5.21 to get

(5.40) $\displaystyle\mathcal{L}_{H}(u)$ $\displaystyle\mathrel{\mathop{\kern
0.0pt=}\limits^{H>1/2}}\left(H-\frac{1}{2}\right)^{2}\left(\int_{\mathbb{R}}e^{\frac{s}{T}}\left(e^{-\frac{\tau}{T}}1_{\tau-s\geqslant
0}-1_{-s\geqslant 0}\right)(s-u)^{H-3/2}1_{s\geqslant u}ds\right)^{2}.$

Rescaling by the positive quantity $\tau$, we get

(5.41) $\displaystyle\mathcal{L}_{H}(\tau u)$
$\displaystyle\mathrel{\mathop{\kern
0.0pt=}\limits^{H>1/2}}\left(H-\frac{1}{2}\right)^{2}\tau^{2H-1}\left(\int_{\mathbb{R}}e^{\frac{\tau
s}{T}}\left(e^{-\frac{\tau}{T}}1_{s\leqslant 1}-1_{s\leqslant
0}\right)(s-u)^{H-3/2}1_{s\geqslant u}ds\right)^{2}$
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0^{+}}}\left(H-\frac{1}{2}\right)^{2}\tau^{2H-1}\left(\int_{\mathbb{R}}\left(1_{s\leqslant
1}-1_{s\leqslant 0}\right)(s-u)^{H-3/2}1_{s\geqslant u}ds\right)^{2}$
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0^{+}}}\tau^{2H-1}\left[(1-u)^{H-1/2}1_{u\leqslant
1}-(-u)^{H-1/2}1_{u\leqslant 0}\right]^{2}.$

Similarly for the term associated to the multiplicative chaos
$\mathcal{M}_{\gamma}(\\{u_{i}\\})$ (Eq. 5.19), we have, once its variables
rescaled by the scale factor $\tau$,

(5.42) $\displaystyle\mathcal{M}_{\gamma}(\\{\tau u_{i}\\})$
$\displaystyle=T^{2n(n-1)\gamma^{2}}\prod_{i<j=1}^{n}\frac{e^{4\gamma^{2}g[\tau(u_{i}-u_{j})]}}{|\tau(u_{i}-u_{j})|_{+}^{4\gamma^{2}}}$
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0^{+}}}\left(\frac{\tau}{T}\right)^{-2n(n-1)\gamma^{2}}e^{2n(n-1)\gamma^{2}g(0)}\prod_{i<j=1}^{n}\frac{1}{|u_{i}-u_{j}|^{4\gamma^{2}}},$

where $g(0)$ is given in Eq. 2.19, such that, for $H\in]1/2,1[$ and
$\gamma^{2}<\min\left(1/4,H/(n-1)\right)$, we have

(5.43)
$\displaystyle\mathbb{E}\left(\delta_{\tau}X_{H>1/2,\gamma}\right)^{2n}$
$\displaystyle=\frac{(2n)!}{2^{n}n!}\int_{\\{u_{i}\\}\in\mathbb{R}^{n}}\mathcal{M}_{\gamma}(\\{u_{i}\\})\prod_{i=1}^{n}\mathcal{L}_{H}(u_{i})du_{i}$
$\displaystyle=\frac{(2n)!}{2^{n}n!}\tau^{n}\int_{\\{u_{i}\\}\in\mathbb{R}^{n}}\mathcal{M}_{\gamma}(\\{\tau
u_{i}\\})\prod_{i=1}^{n}\mathcal{L}_{H}(\tau u_{i})du_{i}$ (5.44)
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0^{+}}}c_{H,\gamma,2n}\frac{(2n)!}{2^{n}n!}\,\left(\frac{\tau}{T}\right)^{2nH-2n(n-1)\gamma^{2}},$

with $c_{H,\gamma,2n}$ provided in Eq. 3.4. This proves the power law
announced in Eq. 3.3. For $H=1/2$, since $\mathcal{L}_{H}$ is rather simple in
this case (Eq. 5.22), we have

(5.45) $\displaystyle\mathcal{L}_{H}(\tau u)$
$\displaystyle\mathrel{\mathop{\kern 0.0pt=}\limits^{H=1/2}}e^{2\frac{\tau
u}{T}}\left(e^{-\frac{\tau}{T}}1_{u\leqslant 1}-1_{u\leqslant 0}\right)^{2}$
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0^{+}}}\left(1_{u\leqslant 1}-1_{u\leqslant 0}\right)^{2},$

which justifies the expression of the constant $c_{1/2,\gamma,2n}$ (Eq. 3.4)
and of the power law announced in Eq. 3.3 for $H=1/2$, with again
$\gamma^{2}<\min\left(1/4,1/(2(n-1))\right)$.

Take now $H\in]0,1/2[$ and write the derivative of the kernel $f$ (Eq. 5.38)
as

(5.46) $\displaystyle f^{\prime}(s)$
$\displaystyle=\frac{1}{T}f(s)-\delta(\tau-s)+\delta(s).$

which makes sense as a distribution. Once again, it can be shown that we can
inject safely this distribution into the expression of $\mathcal{J}_{H}(u)$
(Eq. 5.24), and obtain, while noticing that the Dirac function centered on
$\tau$ entering in Eq. 5.46 has no contribution,

(5.47) $\displaystyle\mathcal{J}_{H}(u)$
$\displaystyle=-\frac{2}{T}\left(H-\frac{1}{2}\right)\int
f(s_{1})f(s_{2})(s_{1}-u)^{H-\frac{1}{2}}(s_{2}-u)^{H-\frac{3}{2}}1_{u\leqslant
s_{1}}1_{s_{1}\leqslant s_{2}}ds_{1}ds_{2}$
$\displaystyle-2\left(H-\frac{1}{2}\right)(-u)^{H-\frac{1}{2}}1_{u\leqslant
0}\int f(s_{2})(s_{2}-u)^{H-\frac{3}{2}}1_{0\leqslant s_{2}}ds_{2}.$

Rescaling by the positive quantity $\tau$, we get

(5.48) $\displaystyle\mathcal{J}_{H}(\tau u)$
$\displaystyle=-\frac{2}{T}\left(H-\frac{1}{2}\right)\tau^{2H}\int f(\tau
s_{1})f(\tau
s_{2})(s_{1}-u)^{H-\frac{1}{2}}(s_{2}-u)^{H-\frac{3}{2}}1_{u\leqslant
s_{1}}1_{s_{1}\leqslant s_{2}}ds_{1}ds_{2}$
$\displaystyle-2\left(H-\frac{1}{2}\right)\tau^{2H-1}(-u)^{H-\frac{1}{2}}1_{u\leqslant
0}\int f(\tau s_{2})(s_{2}-u)^{H-\frac{3}{2}}1_{0\leqslant s_{2}}ds_{2}$
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0^{+}}}-2\left(H-\frac{1}{2}\right)\tau^{2H-1}(-u)^{H-\frac{1}{2}}1_{u\leqslant
0}\int(1_{s_{2}\leqslant 1}-1_{s_{2}\leqslant
0})(s_{2}-u)^{H-\frac{3}{2}}1_{0\leqslant s_{2}}ds_{2}$
$\displaystyle\mathrel{\mathop{\kern 0.0pt\sim}\limits_{\tau\to
0^{+}}}-2\tau^{2H-1}(-u)^{H-\frac{1}{2}}1_{u\leqslant
0}\left[(1-u)^{H-\frac{1}{2}}-(-u)^{H-\frac{1}{2}}\right].$

As encountered when deriving the variance of the process, we cannot make sense
of the expression of $\mathcal{K}_{H}(u)$ (Eq. 5.25) when the derivative of
$f$ (Eq. 5.46) is distributional. Instead, write $\mathcal{K}_{H}$ as

(5.49) $\displaystyle\mathcal{K}_{H}(u)\equiv\lim_{\epsilon\to
0}\left[2\left(H-\frac{1}{2}\right)\int_{\mathbb{R}}f^{2}(s)(s-u+\epsilon)^{2H-2}1_{s\geqslant
u}ds+\epsilon^{2H-1}f^{2}(u)\right]$
$\displaystyle=(\tau-u)^{2H-1}1_{u\leqslant\tau}+(-u)^{2H-1}\left(1-2e^{-\frac{\tau}{T}}\right)1_{u\leqslant
0}-\frac{2}{T}\int_{\mathbb{R}}f^{2}(s)(s-u)^{2H-1}1_{s\geqslant u}ds.$

Once rescaled by $\tau>0$, we obtain

(5.50) $\displaystyle\mathcal{K}_{H}(\tau u)$
$\displaystyle=\tau^{2H-1}\left[(1-u)^{2H-1}1_{u\leqslant
1}+(-u)^{2H-1}\left(1-2e^{-\frac{\tau}{T}}\right)1_{u\leqslant 0}\right]$
$\displaystyle-\tau^{2H}\frac{2}{T}\int_{\mathbb{R}}f^{2}(\tau
s)(s-u)^{2H-1}1_{s\geqslant u}ds$ $\displaystyle\mathrel{\mathop{\kern
0.0pt\sim}\limits_{\tau\to 0^{+}}}\tau^{2H-1}\left[(1-u)^{2H-1}1_{u\leqslant
1}-(-u)^{2H-1}1_{u\leqslant 0}\right].$

Finally, summing up the equivalents of $\mathcal{J}_{H}(\tau u)$ (Eq. 5.48)
and $\mathcal{K}_{H}(\tau u)$ (Eq. 5.50) as $\tau\to 0$, we get

$\displaystyle\mathcal{L}_{H}(\tau u)\mathrel{\mathop{\kern
0.0pt\sim}\limits_{\tau\to
0^{+}}^{H<1/2}}\tau^{2H-1}\left[(1-u)^{H-1/2}1_{u\leqslant
1}-(-u)^{H-1/2}1_{u\leqslant 0}\right]^{2},$

which justifies the expression of the constant $c_{H,\gamma,2n}$ (Eq. 3.4) and
of the power law announced in Eq. 3.3 for $H<1/2$, with again
$\gamma^{2}<\min\left(1/4,H/(n-1)\right)$.

Acknowledgements. L.C. is partially supported by the Simons Foundation Award
ID: 651475.

## References

  * [1] G. B. Apolinário and L. Moriconi.  Shot noise multifractal model for turbulent pseudo-dissipation.  Journal of Statistical Mechanics: Theory and Experiment, 2020(7):073208, jul 2020.
  * [2] A. Arneodo, B. Audit, P. Kestener, and S. Roux.  Wavelet-based multifractal analysis.  Scholarpedia, 3(3):4103, 2008.
  * [3] E. Bacry, J. Delour, and J.-F. Muzy.  Multifractal random walk.  Phys. Rev. E, 64:026103, 2001.
  * [4] E. Bacry and J.-F. Muzy.  Log-infinitely divisible multifractal processes.  Comm. Math. Phys., 236(3):449–475, 2003.
  * [5] O. Barndorff-Nielsen and N. Shephard.  Non-gaussian ornstein–uhlenbeck-based models and some of their uses in financial economics.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(2):167–241, 2001.
  * [6] J. Barral and B. B. Mandelbrot.  Multifractal products of cylindrical pulses.  Prob. Th. Rel. Fields, 124(3):409–430, 2002.
  * [7] B. C. Boniece, G. Didier, and F. Sabzikar.  Tempered fractional brownian motion: Wavelet estimation, modeling and testing.  Applied and Computational Harmonic Analysis, 2019.
  * [8] P. Chainais, R. Riedi, and P. Abry.  On non-scale-invariant infinitely divisible cascades.  IEEE transactions on Information Theory, 51(3):1063–1083, 2005\.
  * [9] P. Cheridito, H. Kawaguchi, and M. Maejima.  Fractional ornstein-uhlenbeck processes.  Electron. J. Probab., 8:no. 3, 1–14, 2003.
  * [10] L. Chevillard.  Regularized fractional ornstein-uhlenbeck processes and their relevance to the modeling of fluid turbulence.  Physical review E, 96(3):033111, 2017.
  * [11] L. Chevillard, C. Garban, R. Rhodes, and V. Vargas.  On a skewed and multifractal unidimensional random field, as a probabilistic representation of kolmogorov’s views on turbulence.  Annales Henri Poincaré, 20(11):3693–3741, Nov 2019.
  * [12] J.-B. Durrive, P. Lesaffre, and K. Ferrière.  Magnetic fields from multiplicative chaos.  Monthly Notices of the Royal Astronomical Society, 496(3):3015–3034, 06 2020.
  * [13] U. Frisch.  Turbulence, The Legacy of A.N. Kolmogorov.  Cambridge University Press, Cambridge, 1995.
  * [14] Y. V. Fyodorov, B. A. Khoruzhenko, and N. J. Simm.  Fractional brownian motion with hurst index $h=0$ and the gaussian unitary ensemble.  Ann. Probab., 44(4):2980–3031, 07 2016.
  * [15] P. Hager and E. Neuman.  The multiplicative chaos of $h=0$ fractional brownian fields.  arXiv preprint arXiv:2008.01385, 2020.
  * [16] A. Jacquier and A. Pannier.  Large and moderate deviations for stochastic volterra systems.  arXiv preprint arXiv:2004.10571, 2020.
  * [17] S. Jaffard, S. Seuret, H. Wendt, R. Leonarduzzi, S. Roux, and P. Abry.  Multivariate multifractal analysis.  Applied and Computational Harmonic Analysis, 46(3):653–663, 2019\.
  * [18] J.-P. Kahane.  Sur le chaos multiplicatif.  Ann. Sci. Math. Québec, 9:105, 1985.
  * [19] A. N. Kolmogorov.  The local structure of turbulence in a incompressible viscous fluid for very large Reynolds number.  Dokl. Akad. Nauk SSSR, 30:299, 1941.
  * [20] B. B. Mandelbrot.  Possible refinement of the lognormal hypothesis concerning the distribution of energy dissipation in intermittent turbulence.  In M. Rosenblatt and C. Van Atta, editors, Statistical Models and Turbulence, volume 12 of Lecture Notes in Physics, pages 333–351. Springer Berlin Heidelberg, 1972.
  * [21] B. B. Mandelbrot.  Fractals and scaling in finance: Discontinuity, concentration, risk. Selecta volume E.  Springer Science & Business Media, 2013.
  * [22] B. B. Mandelbrot and J. W. Van Ness.  Fractional Brownian motion, fractional noises and applications.  SIAM Reviews, 10:422, 1968.
  * [23] E. Neuman and M. Rosenbaum.  Fractional brownian motion with zero hurst parameter: a rough volatility viewpoint.  Electron. Commun. Probab., 23:12 pp., 2018.
  * [24] D. Nualart.  Stochastic integration with respect to fractional brownian motion and applications.  Contemporary Mathematics, 336:3–40, 2003.
  * [25] R. M. Pereira, C. Garban, and L. Chevillard.  A dissipative random velocity field for fully developed fluid turbulence.  J. Fluid Mech., 794:369–408, 5 2016.
  * [26] R. M. Pereira, L. Moriconi, and L. Chevillard.  A multifractal model for the velocity gradient dynamics in turbulent flows.  J. Fluid Mech., 839:430–467, 2018.
  * [27] V. Pipiras and M. Taqqu.  Integration questions related to fractional brownian motion.  Probability Theory and Related Fields, 118(2):251–291, 2000.
  * [28] R. Rhodes, J. Sohier, and V. Vargas.  Levy multiplicative chaos and star scale invariant random measures.  The Annals of Probability, 42(2):689–724, 2014.
  * [29] R. Rhodes and V. Vargas.  Gaussian multiplicative chaos and applications: A review.  Probability Surveys, 11:315, 2014.
  * [30] F. Schmitt and D. Marsan.  Stochastic equations generating continuous multiplicative cascades.  Eur. Phys. J. B, 20(1):3–6, 2001.
  * [31] C. A. Tudor and F. G. Viens.  Statistical aspects of the fractional stochastic calculus.  The Annals of Statistics, 35(3):1183–1212, 2007.
  * [32] B. Viggiano, J. Friedrich, R. Volk, M. Bourgoin, R. B. Cal, and L. Chevillard.  Modelling lagrangian velocity and acceleration in turbulent flows as infinitely differentiable stochastic processes.  Journal of Fluid Mechanics, 900:A27, 2020.
  * [33] H. Wendt, P. Abry, and S. Jaffard.  Bootstrap for empirical multifractal analysis.  IEEE signal processing magazine, 24(4):38–48, 2007.