Moses, Noah and Joseph effects in Lévy walks

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 and Wallis 1968 Water Resour. Res. 4 909), fat-tails of the increment probability density lead to a ‘Noah effect’ (Mandelbrot and Wallis 1968 Water Resour. Res. 4 909), and non-stationarity, to the ‘Moses effect’ (Chen et al 2017 Phys. Rev. E 95 042141). 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évy walk, compare our results to theoretical predictions, and discuss the generality of the method.


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, see e.g. [3][4][5][6][7]. Of-course, 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,8,9], or experimental data obtained from observation of molecules diffusing inside cells, e.g. [7,[10][11][12][13][14]. 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 a fat-tailed increment distribution, or is it because there is an actual trend of inflation in the system? Our analysis below allows 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 continuous-time 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,15]. A paradigmatic process that exhibits this effect is fractional Brownian motion [2,16]. Another cause of anomalous scaling may be that the increment distribution is fat-tailed, in the sense that its second moment is divergent. This is the Noah effect [1,2,15]. A Lévy flight process where the increments are power-law distributed, independent random variables [4,17], 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,15]. A paradigmatic model in this case is scaled Brownian motion [18][19][20]. 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 [15], 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 section 2.
In this manuscript, we investigate these three constitutive effects in a well studied stochastic process called coupled Lévy walk [21]. 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 [17,21,22], 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| can be either ±1 (sometimes referred to as 'genuine Lévy walk' [21]), or a deterministic function of τ , but whose direction is chosen randomly to be either toward the right, along the positivê x-axis (+), or left (−) along the negative axis. The latter, generalized model, is the case studied in detail in this manuscript (see also e.g. [23][24][25][26][27]), and the results include also the constant-velocity case. 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 effect on the overall behavior of the system [21], as the velocity V N during this step does not necessarily have to be distributed like all its predecessors, see e.g. [22]. For more on this point, see section 3. 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 1 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 section 2, 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 section 3, we extend the details on the Lévy walk model. In section 4, 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 section 5 we obtain analytic results for the Moses and Noah effects, and in section 6 for the Joseph. We generalize the model in section 7, and the discussion is provided in section 8.

Story of three exponents: M, L and J
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 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, and it arises if the increment PDF is fat-tailed. Finally, the Joseph exponent can be defined via the sum over the auto-correlation function [15] Δ Δ =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 e.g.Δ ∼ 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 can cause superdiffusion (see discussion on 'long-ranged correlations' e.g. in [28]). 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 appendix 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 Figure 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 1/2 L 1.
In this definition, one can notice that the manifestation of the Noah effect is somewhat different from the case of e.g. a Lévy flight, since the mean of the squared increments is not divergent. In fact, as we explain in detail below, what leads to L = 1/2 in this case, is that the increment PDF has a regime where its shape is 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 equation (19), and the step-velocity is V ∼ τ ν−1 . Here γ = 0.52, ν = 0.5. As explained in section 2, 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.
fat-tailed, but these fat tails have a time-dependent cutoff which is pushed towards ∞ or −∞ as time increseases. The resemblance between this observation, and the source of the Noah effect in its original definition on P.1, is the reason that we can make the association between the two cases and refer to L in equations (5) and (6) throughout this manuscript as the Latent exponent. The upper bound on L, in equation (6), 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 (section 2.2). 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) Note that equation (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 section 4. We define the Joseph exponent also in the spirit of the discrete case, via the scaling of the integral Here 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 appendix A for more explanation. Taking the derivative of the integral with respect toΔ, the autocorrelation function is atΔ 1. For smallΔ, we define f (Δ) ≡ f < (Δ), where f < insures that the autocorrelation 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 section 6 and appendices 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.

Relation between M, L, J and H
Let v(0) ≡ 0, using equation (9) and the Green-Kubo relation [29], the MSD of the process can be written as In equation (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 equation (1), this yields The relation in equation (11) was previously shown to hold empirically in a number of models in [2,15]. 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 appendix A.

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. Let P t (v) be the PDF of finding an increment velocity v at time t, given that the process started at at rest at t = 0. This increment PDF is said to have a single scaling shape, if for any v and t it can be described by a time-independent function W(z β ), such W(z β ) = t β P(v/t β ) and z β = v/t β . In our case, we do not restrict P t (v) to only one such scaling regime, and it can have two different scaling shapes, one in the bulk, and another in the tails. This is a situation not uncommon in anomalous diffusion, which is associated with multifractality, see e.g. [27,[30][31][32][33][34]. If both the mean of |v| and v 2 are taken from the same scaling regime of P t (v), then in this regime and Notice that since 1/2 L 1, equation (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 Table 2. Summary of the different scaling limit of P t (v), that can be found from the Moses M and Latent L exponents, via α, β equations (12) and (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.
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 re-scale the PDF in-order to find the invariant limit, since ). According to equations (12) and (14), 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 [30]. 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 non-integrable 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 symmetric Lévy distribution l ξ,1,0 (v), with 0 < ξ < 2, defined as the inverse-Laplace transform of exp(−|u| ξ ), from u → v [35]. In this case, by definition, there is no Moses effect, and the Noah effect rises since here it can only be quantified by the original definition of L, namely via the time-average of the squared increments of individual time-series [1]. If the increment PDF would have e.g. the scaling shape P t (v) ∼ t −1/ξ l ξ,1,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 equation (5), because at its tails the PDF P t (v) is scaled differently in time with respect to the bulk. Now, the definitions in equations (4) and (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. If it happens that the mean of |v| and v 2 are obtained from different scaling regimes, then again equations (12) and (13) are not valid, but one can use methods such as estimating fractional moments [27,32,34] 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 equation (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 quasi-equilibrium state, see e.g. [36][37][38][39][40][41][42]. The relation in equation (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 observed in the Pommeau-Manneville map [43]. If α > 0 and β = 0, or equivalently L > 1/2 and M = (equation (15)), the limit shape of the increment PDF is given by an infinite-covariant density, see e.g. [27,32,33,[44][45][46][47]. Note that, in both the invariant and the covariant case, and also in the case when the mean-absolute and mean-squared 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 2.

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. Equation (16) means that 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, section 8). 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 equation (2), from equations (16) and (18), . For our example, from equation (19) it follows that the step velocity distribution in the first N − 1 complete steps, when ν < 1, is and it has a similar shape but with Θ(|V| τ ν−1 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 [21], 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 [39,48,49], where V N is determined from the time interval straddling t [50]. With this choice, all the velocities V i , with i = 1, . . . , N are IID, though the duration of the last step is given by equation (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 [39], 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 with v via v = lim Δ→0 v, see section 5. 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: 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 . Here, since in the limit t → ∞ the support of the region v/v c < 1 in equation (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 [39], and φ(v) is in equation (20). The time-invariant asymptotic limit given by I(v) in equation (24) is non-integrable around v = 0, hence it is non-normalizable: As such, this function is the infinite-invariant density of the process [39]. Note that when ν > 1, the two regimes of the PDF, equation (22) simply switch places, but their functional shape remains the same.

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 section 3, 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 section 6 and appendix B. The results of this method correspond to those of a direct measurement of the correlation function, but it is numerically more convenient (see appendix B).
What the data analysis says: without relying on prior knowledge about the underlying process, we found that in the range defined by equations (17) and (21), the Lévy walk data exhibits five separate dynamical phases. These phases are summed-up below and in figure 2. The summation formula, equation (11) is confirmed in all but the '∞' regime.
• In regime A, when γ/2 + 1/2 < ν < γ/2 + 1: H = ν, J = 1, L = 1/2, M = ν − 1/2. Here, the auto-correlation function does not decay withΔ, in equation (9), namely the increments are essentially completely correlated. In this situation, we say that the Joseph effect is maximal, since by definition J can never be bigger than its value here. 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 in the form of t ν−1 P t (v/t ν−1 ) to describe the regime of the PDF which gives rise to the first and second moments of |v| (which is therefore not-fat tailed). Our numerics show that this regime extends also to the range 1 < ν < γ/2 + 1 (and γ < 1). • In regime B, γ < ν < γ/2 + 1/2: H = ν, J = (1 + 2ν − γ)/2, L = 1 − ν + γ/2, M = ν − 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 be present too. Here, the Noah effect means that the scaling shape 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, namely |v| and v 2 correspond to different regimes of P t (v). Accordingly our numerical analysis shows that equation (12) is not valid in this case. 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 C, γ/2 < ν < γ: H = ν, J = (1 + 2ν − γ)/2, L = 1 − γ/2, M = γ − 1/2. Still, all three effects contribute to the anomalous diffusion. Here, just as in regime B, the Joseph effect is present, but is not maximal. In this regime, the Moses and Noah effects are coupled, with the Moses and the Latent exponents obeying equation (15). This suggests that the large fluctuations of the system are described by an infinite-invariant density, equation (12) with α = 1 − γ, β = 0. • 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 infinite-invariant 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 equations (5)-(9) do not hold, and in this regime the decomposition is not valid. We call this the '∞' regime. See appendix D. What we know from the model, in comparison with the data analysis: when γ, ν < 1, equations (23) and (24) describe two different ways to obtain a time-invariant scaling-shape for 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 section 5). As expected from the numerics, the analytic results presented in section 5 show that the bulk function and the infinite-invariant 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 section 6, 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 equation (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 equation (20). The divergence of the MSD in the '∞' regime, was shown analytically in [48,49], further details in section 7 and appendix D.

M & L, and how we obtained them
As explained in section 2, in order to obtain M and L, we need to examine the temporal behavior of the ensemble-time 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 equations (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 f: f = (1/t) t 0 f (t )dt . 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 influence of these occurrences is negligible in the context of the results in this manuscript. This is also confirmed by our numerics.

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 equations (20) and (22), for the mean of |v|, we get . (26) 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 equations (23) and (24) create a distinction between two different cases, depending on whether O(v) = |v| or v 2 is integrable with respect to ρ(v), equation (23), or it is integrable with respect to the infinite-invariant density I(v), equation (24). In the first case, the leading order is obtained by first changing variables: 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 equation (25) in the limit t → ∞, to leading order, hence using equation (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 non-integrable regime, where neither observable is integrable (details below). Figures 3(a)-(c) display simulation results for the temporal behavior 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 equations (26) and (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 figure 4. In figure 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 equations (12) and (13).

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 equations (26) and (27), respectively. The second term in both equations gives the next-to-leading order behavior. This result agrees with the calculation based on equation (28). Similar to the argument in equation (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 equation (4) and (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 β ), equation (12), and here z β =ṽ = v/t ν−1 , so β = ν − 1 and α = 0, in agreement with equations (13) and (30). Figure 5(a) displays the convergence of simulation results of P t (v) at increasing times, rescaled according to equation (12), as function of z β , to the scaling limit equation (23). Note that the Moses effect originates from the diverging mean duration of the Lévy walk steps, namely because τ → ∞ in g(τ ), equation (19), which leads to statistical aging [51].

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 equation (29). The result is equal to the term ∝ t γ−1 in equation (27) (and the second term there is now the next-to-leading order behavior). 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 [39]. Using equations (4) and (5), this yields This regime continuously extends the one introduced in equation (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), equations (12)- (14) are not valid, and α and β are not defined. Figure 5(b) shows that, if we did not know the model, and try to obtain α, β from equation (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.

The non-integrable regime, ν < γ
In this regime neither the first, nor the second moment of |v| are integrable with respect to ρ(v), equation (23). Instead, both the mean velocity, and the mean squared velocity are integrable with respect to the infinite-density. Here, using equations (20), (24) and (29) we get |v| , |v| ∝ t γ−1 , as well as v 2 , v 2 ∝ t γ−1 , so from equations (4) and (5), we find In this case, we associate W(z β ), equation (12), now with the infinite-invariant density I(v), equation (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 equation (15), as it should. Figure 5(c) 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(τ ), equation (19) is divergent, however since ν is small, the displacement in the Levy walk steps χ ∼ τ ν , is almost decoupled from τ . This implies that the MSD of the process now mostly depends on how many steps the walker can have between t = 0 and t, and this in turn is determined only by the value of γ. Consequently, the Hurst exponent in this regime depends only on this parameter too.

J, and how we obtained it
The Joseph exponent depends on the shape of the auto-correlation 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) [52], wavelet decomposition [53] or detrended fluctuations analysis [54]. Additional information on the correspondence between our definition of J and the latter method is given in appendix C.  (12) and (14) and the quantification of the Moses and Noah effects. From equation (14): α = 2L − 1 and β = M + α − 0.5. In the log-log plots (a)-(c), in symbols, we see the rescaled PDF t α+β P t (v), obtained from simulation results of 3 × 10 8 paths, in regimes A, B and D respectively (in regime C the shape of the PDF is behave similar to the last, at increasing times). The measurements were performed at times t = 10 4 (red dots), 10 5 (green dots), 10 6 (blue dots) and 10 7 (black dots).  (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 equation (24). The insets show the same results, but in semi-log plots.
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 equation (33) should not be confused with v 2 in equation (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 [29] when t Δ. The scaling of this function for different types of auto-correlations is discussed in appendix B, where we also show the correspondence between δ 2 , the autocorrelation function, and our definition in equation (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 [15], however, due to a typo it was written '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 reference [48]. Particularly, in that reference, the scaling of δ 2 with respect to time and Δ was calculated for a general shape of g(τ ), with an asymptotic fall-off as in equation (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 [48] Using equations (35) and (36), 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 increment-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ν. Figure 4(c) shows a phase diagram summarizing the different regimes of the Joseph effect, shown in equation (37). Figure 4(d) shows the different regimes of the Hurst exponents which results from the combined effect of the various effects leading to the anomalous diffusion, and calculated using equation (11). Our simulation results for several arbitrary samples of values of ν and γ in these regimes agree with the analytic expectation.

A generalized model
In this section, following [48,49] we extend the model displayed above by introducing a new parameter η. This parameter generalize equation (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 [55,56], 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 [49]. 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 figure 2 does not change, however the regimes themselves may expand or shrink and disappear. Let us look again at the PDF P t (x), of the particles' displacement x at time t. Here, 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 reference [49]. (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 equation (3). After Fourier transform, we getP Letf (k, s) = ∞ 0f (k, t) exp(−st)dt be the Laplace transform off (k, t). In Fourier and Laplace space, from equation (39) we obtainP On comparing equations (40) and (43), we can now obtain the MSD using x 2 (s) = A 0 (s)r 2 (s) + A 2 (s)r 0 (s). (44) Calculating the values on the right-hand side of equation (44), and taking the inverse Laplace transform, one can now derive the MSD. Part of this calculation, performed in [49], 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 [49] x 2 (t) ≈ γ 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 in figure 4(d).
When ν < γ/2 + η, but γ < 1, 2ν > γ, the MSD reads [49] which gives the value of the Hurst exponent as H = ν also similar to figure 4(d). The results shown in figure 4 are for η = 1, whereas equations (43) and (44) are calculated for any value of η. This shows that the power-law dependence of the MSD 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 figure 2 survives, and beyond it we have the non-scaling '∞' regime. When η is very large, regime A expends higher into the realm of ν > 1.

Discussion
Imagine that you get hold of a 'blind' set of data series, containing the positions of an ensemble of random walkers at various times in the interval [0, t]. This data was generated by a Lévy walk model, but you do not have this prior knowledge. Our analysis allows us to uncover the main features of the hidden process that cause its behavior to scale anomalously with time, despite the fact that we do not know what process generated the data. Elucidating the origins of anomalous diffusion observed in experimental data is crucial in order to understand the underlying functioning of the system, and it is studied therefore these days, e.g. using new advanced methods for single-particle tracing [7,13,57,58]. We encourage the verification of our results for example in (but not limited to) future such experiments, in particular e.g. the scaling relation in equations (11) and (12), and consequently its application.
In addition to learning about 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 be at least good candidates to represent the underlying dynamics. These days, there are many studies which use techniques such as machine learning [59][60][61], Bayesian statistics [62] and more, e.g. [63,64], to try to infer the Hurst exponent or 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 analysis of data obtained from single trajectories. This issue is even being studied today as part of a multi-group competition to characterize the properties of anomalous diffusion in data, called the ANDI challenge [65]. Though we cannot fully and uniquely restore the underlying dynamics just by discerning the scaling properties of the process from the data, the characterization of anomalous diffusion using three additional exponents M, L and J, in addition to the Hurst, does bring additional tools which can be helpful for modeling it. In this sense, this decomposition should also be useful for example for the modeling of diffusion in the membranes of living cells done in [13], where a Moses and a Joseph effect seem to have been observed. Another interesting example is found in [66], where the authors observed intercellular transport of insulin granules in eukaryotic cells, and then used information from the time-averaged MSD (Joseph), and the evolution of the absolute-mean of the increments (Moses) to model it. The authors compared two candidate models to describe their dynamics: fractional Brownian motion and continuous-time random walk, and concluded that none is sufficiently good for a full description of the system. They therefore continued by proposing a different, 'hybrid' model based on the previous two [66]. Since the first model leads only a Joseph effect, but the second leads to both Moses and Noah, a full three-effect decomposition here, which takes into account also the inherent relation between them, equation (11), might shed more light on the unified model. Of-course, in any case if one seeks to fully reconstruct the underlying process from the data, a complete knowledge of the entire correlation structure would be required, including two-point, and all the higher order correlations.
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 other processes as well. The results in figures 2 and 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 step duration and 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 also includes other processes such as the Pommeau-Manneville map [15] and ATTM [67], will display similar properties as in the various regimes in figures 2 and 4. These phase diagrams describe their dynamics as well, after change of variables (see e.g. [15]). 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 equation (11). Further generalizations of this summation relation are discussed below, in appendix B. 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 value is subtracted and equation (C1) simplifies to a discrete version of the time-averaged MSD, with non-overlapping windows.
So what is the 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 [72] 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 [73]. 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 [74]). 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.