Moses, Noah and Joseph Effects in Coupled L\'evy Processes

We study a method for detecting the origins of anomalous diffusion, when it is observed in an ensemble of times-series, generated experimentally or numerically, without having knowledge about the exact underlying dynamics. The reasons for anomalous diffusive scaling of the mean-squared displacement are decomposed into three root causes: increment correlations are expressed by the"Joseph effect"[Mandelbrot 1968], fat-tails of the increment probability density lead to a"Noah effect"[Mandelbrot 1968], and non-stationarity, to the"Moses effect"[Chen et al. 2017]. After appropriate rescaling, based on the quantification of these effects, the increment distribution converges at increasing times to a time-invariant asymptotic shape. For different processes, this asymptotic limit can be an equilibrium state, an infinite-invariant, or an infinite-covariant density. We use numerical methods of time-series analysis to quantify the three effects in a model of a non-linearly coupled L\'evy walk, compare our results to theoretical predictions, and discuss the generality of the method.


I. INTRODUCTION
Diffusive processes that scale anomalously with time, such that the Mean-Squared Displacement (MSD) of the expanding particle packet is and the Hurst exponent H = 1/2, are widely observed. This behavior is found both in theoretical models as well as in many experiments, for reviews see e.g. [3][4][5][6]. Ofcourse, if we know the exact underlying process responsible for the dynamics, H can be determined exactly and the various features of the system that lead to the deviation from the standard linear scaling of the MSD, expected by the Gaussian Central Limit Theorem (CLT), can be understood. However, when anomalous diffusive scaling is detected in measurements it is not always clear what is responsible for the observed behavior of the system. Imagine, for example, that we obtain an ensemble of data-series describing intra-day trades in financial markets [2,7,8], or experimental data obtained from observation of molecules diffusing inside cells, e.g., [9][10][11][12]. Here, the proper characterization of the exact root causes of this phenomenon is very important, since it can have implications on how we understand the underlying functioning of the system. If, for example, we observe that the MSD grows faster than linearly with time, is this due to temporal correlations in the data that cause random large fluctuations to be followed by similar or even greater ones? Is it the result of fat-tailed increment distributions, or is it because there is an actual trend of inflation in the system? Our analysis below allow us to give answer to these questions, despite the fact that we cannot completely restore the underlying process just from the data. To make this more precise: Consider a continuoustime stochastic process x(t ) defined in the time interval t ∈ [0, t]. We can choose a number Q of observation windows of duration ∆ = t/Q, and then, represent this process by a discrete time-series composed of consecutive increments, starting at times 0, ∆, 2∆, . . . , (Q − 1)∆ . The increments are δx 1 , δx 2 , . . . , δ Q = x(∆) − x(0), x(2∆) − x(∆), . . . , x(t) − x(t − ∆) . According to the Gaussian CLT, in the limit of large Q, if the increments are independent, identically distributed (IID) random variables chosen from a distribution with finite variance, then the MSD will grow linearly with Q and thus with time. Each of the three ways that the CLT can be violated corresponds to a constitutive effect that can produce anomalous scaling [2]. For processes with stationary increments, where the probability distribution of δx j is independent of time, anomalous diffusive scaling can occur because of long-time increment correlations. This is called the Joseph effect [1,2,13]. A paradigmatic process that exhibits this effect is fractional Brownian motion [2,14]. Another cause of anomalous scaling may be that the increment distribution is fat-tailed, so that its second moment is not a constant (i.e., it is divergent, or changes in time). This is the Noah effect [1,2,13]. A Lévy flight process where the increments are power-law distributed, independent random variables [4,15], but with infinite variance, is one example of a model with this effect. When the increment distribution is non-stationary, anomalous diffusive scaling can also arise due to the Moses effect [2,13]. A paradigmatic model in this case is scaled Brownian motion [16][17][18]. Each of the three effects can appear individually in a system, or in various combinations. Importantly, the three effects can be interconnected with each-other. For example, in [13], it was shown that statistical aging in the process can be associated not only with a Moses, but Noah effect. Among other things, this manuscript will extend our understanding of the coupling of the Moses and Noah effects. The quantification of the three constitutive effects and the relation between them is given in Sec. II.
In this manuscript, we investigate these three constitutive effects in a well studied stochastic process called coupled Lévy walk [19]. This model is known to have a rich spectrum of statistical behaviors, found by the tuning of a few well defined handles. We explore the emergence of the three effects in different parameter regimes of the model using simulations and methods of time-series analysis of single Lévy walk trajectories, and compare our findings with analytical results based on the well developed theory for this process. This example shows that the analysis based on the three constitutive effects is a useful tool that can be applied to study other systems as well (see discussion).
In a two-state Lévy walk [15,19,20], a particle starts at x = 0 at time t = 0 and then moves in independent steps. Each step has a random duration τ , chosen from a Probability Density Function (PDF) of the form at long τ , where c, γ > 0 are constants. During each step, the particle travels at a constant velocity V , whose magnitude |V | is a deterministic function of τ , but whose direction is chosen randomly to be either toward the right, along the positivex−axis (+), or left (−) along the negative axis. The probability of the direction being to the right or to the left is equal, so the motion is unbiased and the velocity has a symmetric PDF φ(V ). At time t = t, the process stops. Up to this point, the particle has made N − 1 complete steps, and one, final "partial" step of duration The properties of the final step have been shown to have a dramatic affect on the overall behavior of the system [19], as the velocity V N during this step does not necessarily have to be distributed like all its predecessors, see e.g., [20]. For more on this point, see Sec. III. The number of steps in the process N ∈ [1, ∞), in the time interval [0, t] is random, and the particle's position at time t is given by the sum Table I summarizes the main notations we use throughout the paper, by order of their appearance in the main text. The structure of the manuscript is as follows: In Sec. II, we define the three exponents that quantify the Moses, Noah and Joseph effects. We discuss the relation between them, and their role in determining the scaling shape of the increment PDF. In Sec. III, we extend the details on the Lévy walk model. In Sec. IV, we provide a summary of our main results, obtained from time-series analysis of numerical simulations, and a brief comparison of these results with the theoretical predictions. In Sec. V we obtain analytic results for the Moses and Noah effects, and in Sec. VI for the Joseph. We generalize the model in Sec. VII, and the discussion is provided in Sec. VIII.  The complete decomposition of the origin of anomalous diffusion presented in the introduction, was originally derived for discrete-time processes [2]. In this case, the process starts at ξ 0 = 0, at n = 0, and evolves in discrete jumps n = 1...N with duration ∆, until time t = N ∆. The particle's position after n steps is denoted ξ n∆ . The Moses effect is quantified by the exponent M , given by the median of the sum, of the absolute value of the time-series increments [2] m t/∆ n=1 |δξ n | ≡ m t/∆ n=1 |ξ n∆ − ξ (n−1)∆ | ∝ t M +1/2 . Here, M = 1/2 yields a linear relation which is similar to normal diffusion. The Noah effect is defined by the scaling of the median of the sum of square-increments, and quantified by the Latent exponent L: m Here again, normal diffusion leads to linear scaling, where M = L = 1/2. If there is no Moses effect, namely M = 1/2, the deviation from this scaling is quantified only by the exponent L. Finally, the Joseph exponent can be defined via the sum over the auto-correlation function [13] ∆ ∆ =0 δξ n δξ n+∆ / (δξ n ) 2 ∝∆ 2J−1 , where 0 ≤ J ≤ 1. Here, starting from an arbitrary time point n, we sum over a discrete lag time ∆ , up to∆ ∼ O(t/10) (∆ is not related to ∆, defined above), and the scaling shape is valid when∆, t 1. When J > 1/2, the correlations decay very slowly with∆, which leads to a divergent sum when∆ → ∞, and superdiffusion (see discussion on "long-ranged correlations" e.g., in [21]). When J ≤ 1/2, the correlation function decays at least as fast as 1/∆, which may lead either to normal diffusion, or in some particular cases to sub-diffusion, see Appen. A.
For a process x(t) in continuous time, we divide the time series into Q non-overlapping observation windows of duration ∆ = t/Q as mentioned in the introduction, and define the average velocity in each time interval v(t ) ≡ |δx j |/∆, where δx j = x(j∆) − x (j − 1)∆ , and (j − 1)∆ < t < j∆. Fig. 1 illustrates the decomposition of a continuous-time random trajectory, into a time-series of N increments of equal duration ∆ t. Now, we can re-write the definition of the Moses effect in terms of the ensemble-time averaged absolute-velocity (when ∆ t) We use here the ensemble mean, instead of the median, since it is a more convenient property to study analytically and numerically, hence we assume by this definition that this mean does not diverge. In the same spirit, the Noah effect is defined via the ensemble-time average of the squared velocity, when ∆ t where the Latent exponent The upper bound on L is true because v 2 ≤ |v| 2 .
Intuitively it means that a tuning of the parameter that leads to a Noah effect beyond L = 1, would automatically increase the scaling exponent of the first moment and therefore lead to aging and a Moses effect, instead of Noah. The lower bound exists because fat tails of the increment distribution, which are described by a Noah effect, can never lead to a slowing down of the process.
In this work, we will assume that also |v| ∝ t M −1/2 and v 2 ∝ t 2L+2M −2 . We address the relation between our definitions and the original time-averaged definitions of these effects, which were derived when the ensemble means could be divergent, below (Sec. II B). Since the ensemble and time averaging procedures are commutative, if we know the first we can immediately obtain the latter via |v| → |v| = (1/t) t 0 |v|(t ) dt which yields = (Const./t) we can find v 2 in a similar way from its ensemble mean, this yields Note that Eq. (7) introduces additional limits on the possible values of M and L, for processes with finite |v| and v 2 , since the ratio between the time and ensemble averages here has to be positive. These limits are consistent with our results for the Lévy walk model, in Sec. IV. We define the Joseph exponent also in the spirit of the discrete case, via the scaling of the integral  III. At the total measurement time t, the last step is incomplete. Two red dash-dot lines mark the start and end points of one completed Lévy walk step, whose duration τ was selected from the PDF Eq. (19), and the step-velocity is V ∼ τ ν−1 . Here γ = 0.52, ν = 0.5. As explained in Sec. II, the trajectory is decomposed into a series of consecutive increments n = 1, 2, ..., of equal duration ∆, the start and end points of one such increment are marked e.g. by two green dash-dot lines. The size of the average velocity |v| in that increment is also presented.
again,∆ should not be confused with ∆, which is the time duration from which we defined v. In this manuscript we will only focus on the case where see Appen. A for more explanation. Taking the derivative of the integral with respect to∆, the autocorrelation function is where f < insures that the autocorreletation function is regularized at∆ → 0. Note that in data analysis there are several known methods to obtain the Joseph exponent without directly calculating the autocorrelation function. These methods have various advantages and disadvantages in practice, see Sec. VI and Appen. B and C.
We note that by dividing δx j by ∆, and defining the three effects via the mean increment velocity v, we did not limit the generality of the definitions at all. The reason is that we did not at this point take the limit ∆ → 0, hence we do not require the instantaneous velocity to be defined. In any process, one can discuss average velocities and increments of a finite-time duration interchangeably.

A.
Relation between M, L, J and H Let v(0) ≡ 0, using Eq. (9) and the Green-Kubo relation [22], the MSD of the process can be written as In Eq. (10), c < is a constant, and in the last step note that since the term ∝ t 2M +2L−1 is subdominant with respect to the other when J > 1/2, we neglected it in the long-time limit. Using Eq. (1), this yields The relation in Eq. (11) was previously shown to hold empirically in a number of models in [2,13]. It was conjectured to be broadly valid, even for systems beyond the case we study here, in particular also when ensemble averages diverge and the Moses and Noah effects are only quantified via their original time-averaged definitions. However, a rigorous derivation in other cases is still needed. For more details see Appen. A.

B. Scaling shapes of the increment distribution
Considering ensemble averages allows us to obtain additional insight about the meaning of the Moses and Noah effects. Assume that |v| and v 2 are not divergent. If both means are taken from a regime in the PDF P t (v) with a single scaling shape, then in this regime P t (v) converges to a time-invariant limit as and Notice that since 1/2 ≤ L ≤ 1, Eq. (6), then 0 ≤ α ≤ 1. The limit function W (z β ) is responsible for the mean of |v| and v 2 via for q = 1, 2.
When M, L are such that both α and β are zero, the increment PDF has a stationary asymptotic (equilibrium) state. Coincidentally this occurs only when M = L = 1/2, which as mentioned means that the time-series satisfies at least two of the conditions of the Gaussian CLT. Curiously, M can also be half if α = β = 0. When P t (v) is non-stationary, we always have a Moses effect. The PDF has a normalized scaling shape, when α = 0 but β = 0, namely L = 1/2, M = 1/2. This is the onset of a "pure" Moses effect. Now, the exponent M tells us how to rescale the PDF in-order to find the invariant limit, Note that usually, based on intuition taken from Gaussian processes, there is a tendency to vaguely associate the Hurst exponent H, with the "self-similarity" property of the process. However in anomalous diffusion that is not necessarily the case; one example is when the MSD is diverging, e.g., in Lévy flight, another example is the case of multifractality [23]. In our case, it is β, not H, that may describe this property, from the point of view of the increment PDF.
The onset of a Noah effect means that v 2 becomes nonintegrable with respect to the scaling function which gives the shape of the P t (v) in the bulk. In the paradigmatic example for this effect, Lévy flight [1], the PDF P t (v) can be e.g., a stationary Lévy distribution l ξ,0 (v), with 0 < ξ < 2, defined as the inverse-Laplace transform of exp(−u ξ ), from u → v [24]. In this case, by definition, there is no Moses effect, and the Noah effect rises since ∞ −∞ v 2 l ξ,0 (v)dv → ∞, though of-course, here it can only be quantified by the original definition of L, namely via the time-average of the squared increments of single timeseries [1]. If the increment PDF would have e.g., the scaling shape P t (v) ∼ t −1/ξ l ξ,0 (v/t 1/ξ ), we would find both a Moses effect, and a Noah effect which is still characterized via the time average.
A more involved scenario that can occur, is when the large fluctuations of the system are reduced such that v 2 is not strictly infinity, but is increasing with time as in Eq. (5), because at its tails the PDF P t (v) is scaled differently in time with respect to the bulk. Now, the definitions in Eqs. (4,5) are valid. The Noah effect will now appear if the function which describes the asymptotic shape of P t (v) at the bulk is fat-tailed (in the sense that its variance is infinite), but the mean v 2 will be given by a second scaling function to which P t (v) convergence at the tails. Such a phenomenon is called multifractality, see e.g., [23]. If it happens that the mean of |v| and v 2 are obtained from different scaling regimes, then again Eq. can be found from the Moses M and Latent L exponents, via α, β Eqs. (12,13), if both |v| and v 2 correspond to the same scaling regimes of the PDF. Note that α, β set the restrictions for M, L in the various regimes, not the other way around.
(12) and Eq. (13) are not valid, but one can use methods such as estimating fractional moments [25][26][27] to find the various scaling shapes of P t (v). If both |v| and v 2 correspond to the second scaling function (that describes the large fluctuations), and are proportional to t M −1/2 and t 2M +2L−2 respectively, then Eq. (12) is valid. But in this case, W (z β ) which denotes these moments might not be normalizable, namely Here, α and the Latent exponent L serves as a measure of how far the increment PDF is from having a normalized limit shape. When α > 0 and β = 0, equivalently L > 1/2 and W (z β ) is an infinite-invariant density, a type of quasiequilibrium state, see e.g., [28][29][30][31][32][33]. The relation in Eq. (15), if observed in data, can in-fact be used to indicate that the underlying process has an infinite-invariant density in this regime, and it was also seen in the Pommeau-Manneville map [34]. If α > 0 and β = 0, or equivalently L > 1/2 and M =[Eq. (15)], the limit shape of the increment PDF is given by an infinite-covariant density, see e.g., [25,26,[35][36][37][38]. Note that, in both the invariant and the covariant case, and also in the case when the mean-absolute and meansquared increments are non-divergent, but they correspond to different scaling regimes of the PDF, a Noah effect cannot appear without a Moses effect. The different cases for M, L and α, β are summarized in Table II.

III. THE LÉVY WALK MODEL
As mentioned in the introduction, in this work we analyse a two-state Lévy walk model. Particularly, here, we consider a continuous range of IID random step velocities, whose distribution is φ(V ). In addition, we assume a nonlinear coupling between the ith step duration and the step velocity, namely where The sign of the step velocity is randomly chosen to be positive or negative with equal probability (the motion is unbiased). The constantc 1 has units of distance/(time) ν , but throughout this manuscript we set c 1 = 1 for convenience. Eq. (16) means that (18) Below, in our numerical simulations, we will use a specific example where the IID random step durations are obtained from the distribution though our results are more general (see the discussion, Sec. VIII). Here, τ 0 > 0 can be as small as we wish, and Θ(·) ≡ 1 when the condition inside the brackets is satisfied and zero otherwise. For any g(τ ) in Eq. (2), from Eqs. (16,18), when |V | < 1 one finds . For our example, from Eq. (19) it follows that the step velocity distribution in the first N − 1 complete steps, when ν < 1, is and a it has a similar shape but with Θ(|V | ≥ τ ν−1 0 ) replacing the original one when ν > 1, hence In this manuscript we focus on the parameter regime where τ is divergent. In various models of non-linearly coupled Lévy walk, some of them are summarized in the review [19], it was shown that in addition to the various scaling exponents, the statistical properties of the process depend strongly on the treatment given to the last, incomplete, step in the sequence. We choose to correspond with the model studied in [31,39,40], where V N is determined from the time interval straddling t [41]. With this choice, all the velocities V i , with i = 1..N are IID, though the duration of the last step is given by Eq. (3). As usual, the displacement at each step (complete and incomplete) is the linear product of the step velocity and its duration. Instantaneous velocity PDF. Akimoto et al. [31], studied the instantaneous velocity PDF P t (v) of the Lévy walker in the process described above, at time t 1 and the regime where 0 < ν < 1. We can apply their results to our analysis, since in this model we can associate v and v via v = lim ∆→0 v, see Sec. V. The following analytic results are brought from that referenced paper. At long but finite times, P t (v) assumes different shapes in two separate ranges of v: Let v c = t ν−1 , then [31] Due to the asymptotic shape of φ(v), when v itself is smaller than unity (regardless of t), P t (v/t ν−1 ) corresponds in this regime to the scaling function ∼ t (ν−1) ρ(ṽ), whereṽ = v/t ν−1 and The scaling function ρ(ṽ) is normalized to unity. On the other hand, at long times P t (v) has a second scaling shape valid in the region v > v c , since in the limit t → ∞ the support of the region v/v c < 1 in Eq. (22) goes to zero, and at v/v c 1, we can expand ] γ as a Taylor series for the small parameter (v/v c ) 1/(ν−1) . This yields, to leading order in time, Which means that asymptotically [31], and φ(v) is in Eq. (20). The time-invariant asymptotic limit given by As such, this function is the infinite-invariant density of the process [31]. Note that when ν > 1, the two regimes of the PDF, Eq. (22) simply switch places, but their functional shape remains the same.

IV. SUMMARY OF OUR MAIN RESULTS
This summary brings the main results of our analysis of Lévy walk trajectories generated by the process described in Sec. III, and the detailed derivations appear below. For further discussion about the generality of the three-effect decomposition, also see below. In our simulations, we generated an ensemble of 10 8 realizations of the process x(t) for different values of γ and ν, and observed the increments δx j of the paths at different times ranging from t = 10 4 to 10 8 . We then measured the ensemble averages of |δx j |, δx 2 j (namely, we used v with observation windows of duration ∆ = 1), as well as x 2 , to calculate the values of M , L and H respectively. To obtain the value of the exponent J, we used a method based on the time-averaged MSD δ 2 , as explained in detail, in Sec. VI and Appen. B. The results of this method correspond to those of a direct measurement of the correlation function, but it is numerically more convenient (see Appen. B).
What the data analysis says: Without relying on prior knowledge on the underlying process, we found that in the range defined by Eqs. (17,21), the Lévy walk data exhibits five separate dynamical phases. These phases are summed-up below and in Fig. 2. The summation formula, Eq. (11) is confirmed in all but the "∞" regime. correlation function does not decay with∆, in Eq. (9). The increments are essentially completely correlated, and the Joseph effect is maximal. There is no freedom left in the increment distribution for any Noah effect to be present. There can be, however, a Moses effect as the increment distribution does "age" with time. The existence of a Moses effect without a Noah effect means that in this regime we expect a single scaling function to describe the regime of t ν−1 P t (v/t ν−1 ) which gives rise to the first and second moments of |v| (it is no-fat tailed). Our numerics shows that this regime extends also to the range 1 < ν < γ/2 + 1 (and γ < 1).
• In regime B, γ < ν < γ/2 + 1/2: In this regime all the three effects contribute to the anomalous diffusion. Here, the Joseph effect is present, but is not maximal, as the auto-correlation function decays as a power-law function of∆. This allows for a Noah effect to now be present too. This means that the scaling function at the bulk of P t (v) is fat-tailed, in the sense that its second moment is divergent. But the mean of |v| remains unchanged from regime A, so it is expected to still be given by the same scaling regime of the increment PDF as before. Accordingly our numerical analysis shows that Eq. (12) is not valid in this regime. The Moses effect occurs here in a similar way as it does in regime A, namely also in this regime, the increment PDF is not time-invariant.
• In regime D, ν < γ/2: H = γ/2, J = 1/2, L = 1 − γ/2, M = γ − 1/2. Here, M, L remain coupled as in region C. Hence we expect the same infiniteinvariant density to be valid in this regime too. Interestingly, now there are no long-range increment correlations and, thus, there is no Joseph effect. At this stage anomalous diffusion occurs due to the non-stationarity of P t (v) and the fat tails of the scaling-shape describing this PDF at the bulk.
• When ν > γ/2 + 1, the MSD is divergent. The scaling relations in Eqs. (5-9) don't hold, and in this regime the decomposition is not valid. We call this the "∞" regime. See Appen. D.
What we know from the model, in comparison with the data analysis: When γ, ν < 1, Eq. (23) and Eq. (24) describe two different ways to obtain a time-invariant scaling-shape of the instantaneous velocity PDF P t (v), the first is valid for small v and the second for large. We can associate this velocity PDF with the distribution of the increment velocity v (see Sec. V). As expected from the numerics, the analytic results presented in Sec. V show that the bulk function and the infiniteinvariant density describe the shape of the increment PDF in regimes A and C, D respectively, in the range of v which is responsible for the various moments. In regime B, |v| , v 2 , (hence |v| , v 2 ) are obtained separately from the two scaling regimes. The fact that the Joseph effect, studied in Sec. VI, is "maximal" in regime A, matches to the fact that the bulk limit-function describing the PDF is thin-tailed, from the same reason that in regime D it is "minimal": if the increments are long-(short-) ranged correlated, their size is more (less) predictable from the first step. Therefore, large fluctuations are less (more) possible.
In regime A, when ν > 1, it is easy to show that one can find similar results for |v| and v 2 as in the case when ν < 1 since as mentioned, the shape of P t (v) is similar to Eq. (22), but with the two regimes for v ≤ 1 and v > 1 switching roles. In addition, here τ ν−1 c becomes a lower, instead of an upper cutoff for the step velocity PDF in Eq. (20). The divergence of the MSD in the "∞" regime, was shown analytically in [39,40], further details in Sec. VII and Appen. D.

V. M & L, AND HOW WE OBTAINED THEM
As explained in Sec. II, in order to obtain M and L, we need to examine the temporal behavior of the ensembletime averages |v| and v 2 , where v is the mean velocity obtained at increments δx, whose duration ∆ is defined independently from step duration of the underlying Lévy walk (namely ∆ = τ ). Choosing ∆ 1, the mean velocity v can be exchanged with the instantaneous velocity v of the random walker at various points in time, and then we can replace v ↔ v in Eqs. (4)(5). Accordingly, this means that we can obtain the exponents of the time series from |v| and |v| 2 , where we now use the following definition for the time average of an observable We note that here one should use a bit of care, since during an increment of duration ∆, the particle might have ended one step of the underlying random walk, and started another, and in this interval of the motion the mean velocity is different from the instantaneous value before/ after the transition. However we assume that if ∆ is small enough, the effect of these occurrences is negligible in the context of the results in this manuscript. This is also confirmed by our numerics.

A. Three regimes for M and L
One can obtain the long-time asymptotic behavior of the ensemble mean of any symmetric observable O(v) in the system, as follows: Given Eqs. (20,22), for the mean of |v|, we get . Similarly, for the mean of v 2 , we get To determine the leading behavior of these two means in the long time limit, note that Eqs. (23,24) create a distinction between two different cases, depending on whether O(v) = |v| or v 2 is integrable with respect to ρ(v), Eq. (23), or it is integrable with respect to the infinite-invariant density I(v), Eq. (24). In the first case, the leading order is obtained by first changing variables: v/t ν−1 →ṽ O(v) = 2t q(ν−1) 1/t ν−1 0 O(ṽ)ρ(ṽ)dṽ + 2t q(ν−1) ∞ 1/t ν−1 O(ṽ)P t (ṽt ν−1 )dṽ, and then in the range t 1 the first term is ≈ 2t q(ν−1) ∞ 0 O(ṽ)ρ(ṽ)dṽ and the second term approaches zero since its support vanishes. So, in this case In the second case, when is integrable with respect to the infinite-invariant density, the contribution to its mean from the region v < v c can be neglected in Eq. (25) in the limit t → ∞, to leading order, hence using Eq. (24) we get Notice that in this case, the temporal scaling of O(v) is similar for all the integrable observables, since it is determined only by the scaling of the infinite-density. For O(v) = |v| and v 2 together, there are three regimes of behavior, included within the range γ, ν < 1: The integrable regime, where both |v| and v 2 are integrable with respect to ρ(v); The middle regime, where only the mean-absolute velocity is integrable; And the nonintegrable regime, where neither observable is integrable (details below). Figs. 3a-c display simulation results for the temporal behaviour of |v| and v 2 in the integrable, the middle and the non-integrable regimes, respectively. The simulations match perfectly at long times with the exact expressions in Eqs. (26,27), which denote both the leading order behavior of |v| and v 2 in time, and the next-to-leading order. They also confirm the approach to the leading order asymptotic results, though this approach is slow. The results for the exponents M and L in the various regimes are shown in the lower two panels of Fig. 4. In Fig. 5, we use the results for these exponents in the three regimes, in order to seek for a time invariant asymptotic shape of P t (v), based on Eqs. (12,13).
B. The integrable regime, 1/2 + γ/2 < ν < 1 In this regime, the leading behavior in time of |v| and v 2 is given by the ∼ t ν−1 and ∼ t 2ν−2 terms in Eq. (26) and Eq. (27), respectively. The second term in both equations gives the next-to-leading order behavior. This result agrees with the calculation based on Eq. (28). Similar to the argument in Eq. (7), the ensemble-time averages |v| ∝ t ν−1 and v 2 ∝ t 2ν−2 , like their corresponding ensemble averages, and since we associate v with v, we now obtain the Latent and Moses exponents using Eqs. (4,5): Since here both means are obtained from the same scaling limit of P t (v), we can now associate ρ(ṽ) in this regime with W (z β ), Eq. (12), and here z β =ṽ = v/t ν−1 , so β = ν − 1 and α = 0, in agreement with Eqs. (13,30). Fig. 5a displays the convergence of simulation results of P t (v) at increasing times, rescaled according to Eq. (12), as function of z β , to the scaling limit Eq. (23). Note that the Moses effect originates from the diverging mean duration of the Lévy walk steps, namely because τ → ∞ in g(τ ), Eq. (19), which leads to statistical aging [42].
C. Middle regime, γ < ν < 1/2 + γ/2 In this regime, v 2 is no longer integrable with respect to the scaling function ρ(v). |v|, however, still is. Therefore, the leading-order behavior of |v| and |v| , remains proportional to ∼ t ν−1 , similar to the previous, integrable region. However since v 2 is now integrable with respect to the infinite-density I(v) instead of ρ(v), its leading behavior is now obtained from Eq. (29). The result is equal to the term ∝ t γ−1 in Eq. (27) (and the second term there is now the next-to-leading order behaviour). Therefore, also v 2 ∼ t γ−1 . Note that in this regime we can obtain the time average of v 2 (t) also using arguments based on infinite-ergodic theory [31]. Using Eqs. (4,5), this yields This regime continuously extends the one introduced in Eq. (30). The First moment can still be described by the scaling shape of the PDF at the bulk. However, since the second moment of this PDF diverges with respect to ρ(v), here we see for the first time the emergence of a Noah effect, in addition to Moses. Since the mean of |v| and v 2 are obtained from two different scaling regimes of P t (v), Eqs. (12)(13)(14) are not valid, and α and β are not defined. Fig. 5b shows that, if we did not know the model, and try to obtain α, β from Eq. (13) in this regime from the data, we would find α = γ − 2ν + 1, β = γ − ν, but with this rescaling, P t (v) does not converge to a time-invariant shape.

D. The non-integrable regime, ν < γ
In this regime neither the first, nor the second moment of |v| are integrable with respect to ρ(v), Eq. (23). Instead, both the mean velocity, and the mean squared velocity are integrable with respect to the infinite-density. Here, using Eqs. (20,24,??) we get |v| , |v| ∝ t γ−1 , as well as v 2 , v 2 ∝ t γ−1 , so from Eqs. (4,5), we find In this case, we associate W (z β ), Eq. (12), now with the infinite-invariant density I(v), Eq. (24), and α = 2L − 1, β = 0, so z β = v. The Noah effect tells us that the asymptotic shape of the increment PDF is given by a non-normalizable function, and the relation between M and L here also agrees with Eq. (15), as it should. Fig.  5c shows how simulation results of t α P t (v) at converge increasing times to I(v), the infinite invariant density.
As in the other regimes, here the mean duration of the Lévy walk steps in g(τ ), Eq. (19) is divergent, however, the since ν is small, the step velocity decays very quickly with the duration. Therefore the step displacement χ ∼ τ ν , is almost decoupled from τ . This implies too things: First, the MSD of the process now mostly depends on how many steps the walker can have between t = 0 and t, and that is determined only by the value of γ. So the Hurst exponent in this regime depends only on γ. Second, by a hand-waving argument we can see why M and L depend only γ; because if the step displacement depends only on this parameter, the average velocity v in all the time-series increments withing those steps will depend only on this parameter too.

VI. J, AND HOW WE OBTAINED IT
The Joseph exponent depends on the shape of the autocorrelation function. However, this quantity is difficult to obtain for many systems, analytically and numerically. In practice, the Joseph effect is often quantified by designated methods, such as the so-called rescaled range statistic (R/S) [43], wavelet decomposition [44] or detrended fluctuations analysis [45]. Additional information on the correspondence between our definition of J and the latter method is given in appendix C.
Here, for the Lévy process, we use a measure which is easier to handle analytically; the ensemble averaged time-averaged MSD δ 2 , defined as Note that Eq. (33) should not be confused with v 2 in Eq. (5), since in the latter the increments are strictly non-overlapping, whereas in δ 2 they are. This quantity is related to the auto-correlation function, via [22] when t ∆. The scaling of this function for different types of auto-correlations is discussed in Appen. B, where we also show the correspondence between δ 2 , the autocorrelation function, and our definition in Eq. (9). In all the cases considered in the appendix (even for J ≤ 1/2), the asymptotic scaling is Our model is described by type (II) in the appendix. This means that the Joseph exponent is given by the scaling of δ 2 with the lag time ∆ (note, that this observation was already made in [13], however, due to a typo it was read 't' instead of '∆').
To obtain the Joseph exponent in various regimes, we use the results of the calculation of the ensemble-time averaged MSD, obtained for this model in Ref. [39]. Particularly, in that Ref., the scaling of δ 2 with respect to time and ∆ was calculated for a general shape of g(τ ), with an asymptotic fall-off as in Eq. (19), at large τ , and it was shown to not depend on the exact behavior at small τ s. Given this knowledge, the time-averaged MSD (for 0 < {γ, ν} < 1) for the Lévy walk model we study here, has the following scaling [39] Using Eq. (36) and Eq. (35), we find that Note that since γ < 1, the mean step duration in all these regimes diverges. This implies that the random walker essentially walks in the same direction for almost all of the time t, regardless how long it is. In turn, this means the process is correlated in the whole parameter regime that we study. But when ν < γ/2, the average correlations decay rather quickly with ∆ because the step velocity changes only very little with the step duration, hence in this regime we do not see a Joseph effect (the difference between the mean velocity at increments belonging to the same steps of the Lévy walk, versus increments of other steps, is small). The onset of the effect is above the line ν = γ/2. It is maximal when J = 1, at γ + 1 < 2ν.  15), hence we expect to find that the increment PDF approaches the shape of a non-normalizable infinite-invariant density. This is confirmed by the solid mustard line, that represents Eq. (24). The insets show the same results, but in semi-log plots. Fig. 4c shows a phase diagram summarizing the different regimes of the Joseph effect, shown in Eq. (37). Fig.  4d shows the different regimes of the Hurst exponents which results from the combined affect of the various effects leading to the anomalous diffusion, and calculated using Eq. (11). Our simulation results for several arbitrary samples of values of ν and γ in these regimes agree with the analytic expectation.

VII. A GENERALIZED MODEL
In this section, following [39,40] we extend the model displayed above by introducing a new parameter η. This parameter generalize Eq. (16) by modifying the relation between the ith step velocity V i , its duration τ i and the actual time in motion t , as follows Some values of η correspond to special cases: The Lévy walk we studied above corresponds to η = 1, when η = ν we get a Drude-like model [46,47], and when η → 0 or η → ∞ we approach either a jump-then-wait type of coupled continuous-time random walk, or a wait-then-jump model, respectively [40]. As we will now show, modifying this parameter changes the onset of the "∞" regime. Our simulation results suggest that when η is within the open range (0, ∞), the behavior of all the effects in regimes A,B,C,D in Fig. 2 does not change, however the regimes themselves may expand or shrink and disappear.
Let's look again at the PDF P t (x), of the particles' displacement x at time t. Here, (39) where A(x , t ) is the joint probability density to land on x between x and x + dx in a complete step ending at t < t, and r(x − x |t − t ) is the conditional probability density of the displacement in the last, incomplete step given the duration of the walk is t. The following calculation of the MSD for this model is adapted from Ref. [40]. Letf (k, t) = ∞ −∞ f (x, t) exp(−ikx)dx, be the Fourier transform of some function f (x, t), from x → k. Eqs. (40,41) and Eq. (42) below, represent the characteristic functions of the probability densities P, r and A, respectively. All the functions except for P lack normalization on unity; the zero terms of the expansions are denoted as r 0 (t) = r(x|t)dx = 1 and A 0 (t) = A(x, t)dx = 1. Let x 2 (t) ≡ x 2 (t) be MSD at time t, A 2 (τ ) = χ 2 A(χ, τ )dχ is the marginal second moment of displacement χ in a single complete step of duration τ and r 2 (τ * ) ≡ χ * 2 r(χ * |τ * )dχ * is the MSD of the displacement χ * in the last, incomplete step. The duration τ * of the latter is defined in Eq. (3). After Fourier transform, we get Letf (k, s) = ∞ 0f (k, t) exp(−st)dt be the Laplace transform off (k, t). In Fourier and Laplace space, from Eq. (39) we obtain On comparing Eq. (40) and Eq. (43), we can now obtain the MSD using Calculating the values on the right-hand side of Eq. (44), and taking the inverse Laplace transform, one can now derive the MSD. Part of this calculation, performed in [40], was to obtain the second marginal moment of the function r(x|t): From here, we can see that r 2 , and therefore also x 2 , can only obtain a finite value when γ > 2(ν − η). This explains the crossover to the "∞" regime, which occurs when ν > γ/2 + η. When ν < γ/2 + η and γ < 1, 2ν < γ, the MSD is [40] x where B(a, b) is the Beta-function. This is dominated by the second term, since 2ν < γ, and therefore x 2 (t) ∝ t γ , which gives the value of the Hurst Exponent as H = γ/2, similar to what is seen seen in Fig. 4d. When ν < γ/2+η, but γ < 1, 2ν > γ, the MSD reads [40] x which gives the value of the Hurst Exponent as H = ν also similar to Fig. 4(d). The results shown in Fig. 4 are for η = 1, whereas Eqs. (43,44) are calculated for any value of η. This shows that the power-law dependence of the mean squared displacement is independent of the exponent η and the particular value of η only enters in the prefactors, when x 2 is finite. When η is very small, only regime D in Fig. 2 survives, and beyond it we have the non-scaling "∞" regime. When η is very large, regime A expends higher into the realm of ν > 1.

VIII. DISCUSSION
Imagine that you get hold of a "blind" set of data series, containing the positions of a random walker at various times in the interval [0, t]. This data was generated by such a Lévy walk, but you do not have this prior knowledge. Our analysis allow us to uncover the main features of the hidden process that cause its behavior to scale anomalously with time, despite the fact that we cannot restore the underlying process just from the data. We encourage the verification of our results in future experiments, in particular e.g. the scaling relation in Eq. (11) and Eq. (12), and consequently its application.
In addition to learning on the origins of the anomalous diffusion, one may use the knowledge about the Moses, Noah and Joseph effects in order to try to extrapolate which processes can and cannot at least be good candidates to represent the underlying dynamics. These days, there are many studies which use techniques such as machine learning [48,49], Bayesian statistics [50] and more, e.g. [51], to try to distinguish between various known models such as continuous-time random walk, fractional Brownian motion and others, which lead to anomalous scaling of the MSD, based only on data obtained from single trajectories generated in a simulation. The characterization of anomalous diffusion using three additional exponents M, L and J, in addition to the Hurst, and the various time-series-analysis based methods to obtain them that we presented in this paper, may add additional tools to be used for such purpose.
The Lévy walk model that we studied in this paper, is a prototypical example which shows how the three effects analysis can be used for many processes in continuous time. The results in Fig. 2 and Fig. 4 eventually only depend on two inputs: the shape of the step durations PDF at large τ s, and the coupling between the step durations and the velocity, which can also be translated to the coupling between the duration and the step displacement, since the step-displacement χ = V τ . Therefore, a class of process which can be mapped into a coupled step-duration and step-displacement process, which includes other processes such as the Pommeau-Manneville map [34] and ATTM [? ], will display similar properties as in the various regimes in Fig. 2 and Fig. 4. These phase diagrams describe their dynamics as well, after change of variables (see e.g., [34]). proof is still required in a future work. The following example, shows that we can prove this relation analytically also for a widely useful system where J ≤ 1/2 (in the derivation in Eq. (10) we assumed that J > 1/2), and the correlation function is negative. Consider the ARFIMA(0, d, 0) process [21,53] x(t) = t i=0 X i , in discrete time, whose increments are defined via the transformation (1 −B) d X i = σ 2 η(i), where η i is Gaussian white noise with zero mean and η i η j = δ ij . Here,B n X i = 2! B 2 + .... When −1/2 ≤ d ≤ 0, this process is long-ranged anti-correlated, and the autocorrelation function of the increments of this process is [21]: 2d−1 , so according to the definition in Eq. (9), the Joseph exponent is J = d + 1/2. The MSD is related to the correlation function via Plugging Eq. (A1) into Eq. (A2), we find that at long t and from the leading-order term, using J = d − 1/2, we find x 2 (t) ∝ t 2J , namely H = J. This is a well known result, and note that in the standard ARFIMA process the increment distribution is stationary and thintailed, hence M = L = 1/2 and the summation relation in the section title, and Eq. (11) is fulfilled. Now, consider the related process: Here we introduced the time dependence of the variance of the increments in the same way is in Eq. (5), which means that now the process can have both a Noah and a Moses effect, in addition to Joseph. A calculation in this case, which follows exactly the same lines as the above, will now yield x 2 (t) ∼ t 2L+2M +2J−2 , which again leads to Eq. (11). Further generalizations of this summation relation are discussed below, in Appen. B.
The conditions above are necessary in order to ensure that the correlation function decays with ∆, and at the same time does not blow up with time t. Since in the q = 0 case, the correlation function is equal to the velocity displacement Eq. (5), continuity demands Now, the correlation function for q → 0 can be inserted into Eq. (B2) Integration yields with t 0 t 1 , t 2 Inserting this result into Eq. (B1) only yields an additional pre-factor. Therefore, the EATA MSD scales like (III) The third case occurs for processes with correlations that do scale with the measurement time t, but decay faster with the lag time ∆. Here, we write the correlation function as with a positive function Φ(∆), which decays faster than ∆ −1 , i.e. the integral over the autocorrelation function with respect to ∆ becomes finite for ∆ → ∞. The dependence on H can again be verified using Eq. (10).
Here, we specify the shape to be Φ (∆) = (1 + ∆) 2ϕ−2 with ϕ < 1/2. The calculation is as simple for the relevant cases of exponential decay or the velocity correlation function being (Dirac-) delta-distributed with δ(∆). Using Eqs. (B2) and (B11) we find for t 0 t 1 , t 2 Since ϕ < 1/2, for ∆ → ∞, the linear term is dominant. Accordingly, the scaling of the EATA MSD for short range correlated processes is The Hurst exponent is independent of ϕ and of the exact shape of the distribution as long as it decays sufficiently fast. The Joseph exponent is therefore always J = 1/2. (IV) In all cases discussed so far, the Joseph exponent is J ≥ 1/2. This is also a necessary condition for the derivation in equation (10). However, it is possible, to construct antipersistent processes with a autocorrelation functions, that lead to J < 1/2. Looking at the calculation (B12) we see, that in order to obtain a scaling with H < 1/2, the constant factor after the first integration, which leads to the linear scaling, has to be eliminated. This is possible if the correlation function is not strictly positive. The most prominent example of a continuous process, which can fulfill this condition is fractional Gaussian noise, which is the continuous-time version of the discrete ARFIMA(0, d, 0) studied in Appen. A. This process has a correlation function similar in shape as in case III above, but without the dependence on t, it reads v(t)v(t + ∆) = 1 2 |∆ + 1| 2J − 2|∆| 2J + |∆ − 1| 2J .
In conclusion we want to point out that the exponent J is given by the scaling of the autocorrelation function with ∆ 2J−2 directly for J > 1/2, and by the scaling of the integral over the autocorrelation function with ∆ 2J−1 for all cases. There are several methods that are used in practice in order to quantify long range correlations in discrete measured time series. Examples are R/S statistics, detrended moving averages [55], scaling analysis based on the wavelet transform [44] and detrended fluctuation analysis [45]. In [56] it was shown that the latter three can be expressed in the same framework. In this section we want to discuss whether or not the definition of long range correlations in DFA is different from our definition, Eq. (9) for J.
The squared fluctuation function of DFA is defined as Here p t,q is a polynomial that is fitted to the series x t for segments of lengths s. The squared error (x t − p t,q ) 2 of these fits is then averaged over all the non-overlapping segments of equal length s. The index q is the order of the polynomial. DFA is not sensitive to trends with polynomial shape of order q − 1, i.e. the slow dynamics is filtered by the method. For the definitions used above, data with trends does not yield a meaningful exponent. So we want to concentrate on stationary data. If the detrending order is zero, in each window just the mean So what is difference if detrending is performed with q > 0? Here, for stationary systems, a relation between the fluctuation function and the autocorrelation function was derived in [57]  (C2) L q (0, s) is some sophisticated kernel. Its leading order is linear. So the fluctuation function is a measure of the integral (here discrete) over the autocorrelation function. Thus it measures the Joseph effect for 0 < J < 1.
If a Noah effect is present, i.e. the variance is infinite in theory, the formula can still be used, since a measured time series always has a finite variance [58]. The Moses effect is more complicated. Even though a scaling of the increment distribution as in scaled Brownian motion (or in the parameter range A in the Lévy walk) is not visible in DFA due to the averaging over the segments, the scaling exponent DFA still might differ from J. This is true if v is diffusive as in fractional Brownian motion (DFA results shown in [59]). Here the DFA exponent is > 1 in contrast to J. In fact it is equal to H. So DFA is usually, but not always, a measure of the Joseph effect.
Appendix D: The "∞" regime, ν > γ/2 + 1 As explained in Sec. VII, the analytic results show that when η = 1 and ν > γ/2 + 1, the Hurst exponent diverges. We demonstrate this in Fig. 6a, by showing simulation results for the time derivative of the ratio log( x 2 )/ log(t) for five pairs of ν, γ in the different regions shown in Fig. 4d. These simulation were performed with an ensemble size of 10 7 trajectories. Clearly, in all but the "∞" regime, this ratio is constant. We expect to find a similar behaviour if we fix time and increase the ensemble size. In Fig. 6b-d, we checked the time evolution of log( δ 2 )/ log(∆), log √ v 2 |v| / log(t) + 1/2 and log( |v| )/ log(t) + 1/2 respectively, at increasing times, which yield the Joseph, Moses and Latent exponents when these ratios are constant (again, as seen in Fig. 4). We observe that in panels (b) and (c) of Fig. 6 the above ratios are constant in all the regions except for the "∞" regime. In the latter, the Joseph and the Latent exponents are clearly not well defined, which renders also the definition of M at least irrelevant. A similar behavior will appear also when η = 1.