Just two moments ! A cautionary note against use of high-order moments in multifractal models in hydrology

The need of understanding and modelling the space–time variability of natural processes in hydrological sciences produced a large body of literature over the last thirty years. In this context, a multifractal framework provides parsimonious models which can be applied to a widescale range of hydrological processes, and are based on the empirical detection of some patterns in observational data, i.e. a scale invariant mechanism repeating scale after scale. Hence, multifractal analyses heavily rely on available data series and their statistical processing. In such analyses, high order moments are often estimated and used in model identification and fitting as if they were reliable. This paper warns practitioners against the blind use in geophysical time series analyses of classical statistics, which is based upon independent samples typically following distributions of exponential type. Indeed, the study of natural processes reveals scaling behaviours in state (departure from exponential distribution tails) and in time (departure from independence), thus implying dramatic increase of bias and uncertainty in statistical estimation. Surprisingly, all these differences are commonly unaccounted for in most multifractal analyses of hydrological processes, which may result in inappropriate modelling, wrong inferences and false claims about the properties of the processes studied. Using theoretical reasoning and Monte Carlo simulations, we find that the reliability of multifractal methods that use high order moments ( > 3) is questionable. In particular, we suggest that, because of estimation problems, the use of moments of order higher than two should be avoided, either in justifying or fitting models. Nonetheless, in most problems the first two moments provide enough information for the most important characteristics of the distribution.


Introduction
A simple way to understand the extreme variability of several geophysical processes over a practically important range of scales is offered by the idea that the same type of elementary process acts at each relevant scale.According to this idea, the part resembles the whole as quantified by so-called "scaling laws".Scaling behaviours are typically represented as power laws of some statistical properties, and they are applicable either on the entire domain of the variable of interest or asymptotically.If this random variable represents the state of a system, then we have the scaling in state, which refers to marginal distributional properties.This is to distinguish from another type of scaling, which deals with time-related random variables: the scaling in time, which refers to the dependence structure of a process.Likewise, scaling in space is derived by extending the scaling in time in higher dimensions and substituting space for time (e.g.Koutsoyiannis et al., 2011).
The scaling behaviour widely observed in the natural world (e.g.Newman, 2005) has often been interpreted as a tendency, driven by the dynamics of a physical system, to increase the inherent order of the system (self-organized criticality): this is often triggered by random fluctuations that are amplified by positive feedback (Bak et al., 1987).In another view, the power laws are a necessity implied by the asymptotic behaviour of either the survival and autocovariance function, describing, respectively, the marginal and joint distributional properties of the stochastic process which models the physical system.The main question is whether the two functions decay following an exponential (fast decay) or a power-type law (slow decay).We assume the latter Published by Copernicus Publications on behalf of the European Geosciences Union.

F. Lombardo et al.: Just two moments!
to hold in the form of scaling in state (heavy-tailed distributions) and in time (long-term persistence), which have also been verified in geophysical time series (e.g.Markonis and Koutsoyiannis, 2013;Papalexiou et al., 2013).According to this view, scaling behaviours are just manifestations of enhanced uncertainty and are consistent with the principle of maximum entropy (Koutsoyiannis, 2011).The connection of scaling with maximum entropy constitutes also a connection of stochastic representations of natural processes with statistical physics.The emergence of scaling from maximum entropy considerations may thus provide theoretical background in modelling complex natural processes by scaling laws.
In the literature, natural processes showing scaling behaviour are often classified as multifractal systems (i.e.multiscaling) that generalize fractal models, in which a single scaling exponent (the fractal dimension) is enough to describe the system dynamics.For a detailed review on the fundamentals of multifractals, the reader is referred to Schertzer and Lovejoy (2011).
Multifractal models generally provide simple power-law relationships to link the statistical distribution of a stochastic process at different scales of aggregation.All power laws with a particular scaling exponent are equivalent up to constant factors, since each is simply a scaled version of the others.Therefore, the multifractal framework provides parsimonious models to study the variability of several natural processes in geosciences, such as rainfall.Rainfall models of multifractal type have, indeed, for a long time been used to reproduce several statistical properties of actual rainfall fields, including the power-law behaviour of the moments of different orders and spectral densities, rainfall intermittency and extremes (see e.g.Koutsoyiannis and Langousis, 2011, and references therein).However, published results vary widely, calling into question whether rainfall indeed obeys scaling laws, what those laws are, and whether they have some degree of universality (Nykanen and Harris, 2003;Veneziano et al., 2006;Molnar and Burlando, 2008;Molini et al., 2009;Serinaldi, 2010;Verrier et al., 2010Verrier et al., , 2011;;Gires et al., 2012;Veneziano and Lepore, 2012;Papalexiou et al., 2013).In fact, significant deviations of rainfall from multifractal scale invariance have also been pointed out.These deviations include breaks in the power-law behaviour (scaling regimes) of the spectral density (Fraedrich and Larnder, 1993;Olsson, 1995;Verrier et al., 2011;Gires et al., 2012), lack of scaling of the non-rainy intervals in time series (Veneziano and Lepore, 2012;Mascaro et al., 2013), differences in scaling during the intense and moderate phases of rainstorms (Venugopal et al., 2006), and more complex deviations (Veneziano et al., 2006;Marani, 2003).
Multifractal signals generally obey a scale invariance that yields power law behaviours for multi-resolution quantities depending on their scale .These multi-resolution quantities at discrete time steps (j = 1, 2,. . .), denoted by x ( ) j in the following, are local time averages in boxes of size (notice that we use the so-called Dutch convention according to which random variables are underlined; see Hemelrijk (1996), and the additional notational conventions in Koutsoyiannis (2013)).This is the basis of the fixed-size boxcounting approach (see e.g.Mach et al., 1995).
For multifractal processes, one usually observes a powerlaw scaling of the form where E[•] denotes expectation (ensemble average) and K(q) is the moment scaling function, at least in some range of scales and for some range of orders q.Generally, the multifractal behaviour of a physical system is directly characterized by the multiscaling exponents K(q), whose estimation relies on the use of the sample q-order moments at different scales and their linear regressions in log-log diagrams.
A fundamental problem in the multifractal analysis of data sets is to estimate the moment scaling function K(q) from data (Villarini et al., 2007;Veneziano and Furcolo, 2009).Considerable literature has been dealing with estimation problems in the context of so-called scaling multifractal measures for at least three decades (see e.g.Grassberger and Procaccia, 1983;Pawelzik and Schuster, 1987;Schertzer and Lovejoy, 1992;Ashkenazy, 1999;Mandelbrot, 2003;and Neuman, 2010).Interestingly, Mandelbrot (2003) and Neuman (2010) recognize the crucial role played by time dependence in estimating multifractal properties from finite length data.Nonetheless, in this work we remain strictly within the framework of the standard statistical formalism, which is actually a novelty with respect to the literature cited above.In this context, we highlight the problematic estimation of moments for geophysical processes, because the statistical processing of geophysical data series is usually based upon classical statistics.The classical statistical approaches rely on several simplifying assumptions, tacit or explicit, such as independence in time and exponentially decaying distribution tails, which are invalidated in natural processes thus causing bias and uncertainty in statistical estimations.In many studies, it has been a common practice to neglect this problem, which is introduced when the process exhibits dependence in time and is magnified when the distribution function significantly departs from the Gaussian form, which itself is an example of an exceptionally light-tailed distribution.In their pioneering work on statistical hydrology, Wallis et al. (1974) already provided some insight into the sampling properties of moment estimators when varying the marginal probability distribution function of the underlying stochastic process.The main results of the paper agree well with those found here, but its Monte Carlo experiments were carried out under a classical statistical framework assuming independent samples.The purpose of this paper is to explore, at different timescales, the information content in estimates of raw moments of processes exhibiting temporal dependence (see Sect. 2).In order for the true moments to be fully known a priori, we use synthetic examples in a Monte Carlo simulation framework.We explore processes with both normal and non-normal distributions including ones with heavy tails.We show (Sect.3) that, even in quantities whose estimates are in theory unbiased, the dependence and non-normality affect significantly their statistical properties, and sample estimates based on classical statistics are characterized by high bias and uncertainty.

Local average process
Central to the development of robust multifractal models is the concept of "local average" of a stochastic process.Practical interest often revolves around local average or aggregates (temporal or spatial) of random variables, because it is seldom useful or necessary to describe in detail the local pointto-point variation occurring on a microscale in time or space.Even if such information were desired, it may be impossible to obtain: there is a basic trade-off between the accuracy of a measurement and the (time or distance) interval within which the measurement is made (Vanmarcke, 1983).For example, rain gauges (owing to size, inertia, and so on) measure some kind of local average of rainfall depth over time.Moreover, through information processing, "raw data" are often transformed into average or aggregate quantities such as, e.g.sub-hourly averages or daily totals.
Mathematically, let x(t) be a stationary stochastic process in continuous time t with mean µ = E [x], and autocovariance c(τ ) = Cov[x(t), x(t + τ )], where τ is the time lag.Consider now the random process x ( ) j obtained by local averaging x(t) over the window at discrete time steps j (=1, 2,. . .), defined as where n = T / is the number of the sample steps of x ( ) j in the observation period T o , and T = T o is the observation period rounded off to an integer multiple of .The relationship between the processes x(t) and x ( ) j is illustrated in Fig. 1.
The mean of the process x ( ) j is not affected by the averaging operation, i.e. (3) Let us now investigate the climacogram of the process x ( ) j , which is defined to be the variance (or the standard deviation) of the time-averaged process x ( ) j as a function of the timescale of averaging (Koutsoyiannis, 2010).The climacogram of x ( ) j can be calculated from the autocovariance function c(τ ) of the continuous-time process as follows (see e.g.Vanmarcke, 1983, p. 186;Papoulis, 1991, p. 299): which shows that the climacogram γ ( ) generally decreases with and fully characterizes the dependence structure of x(t).The climacogram γ ( ) and the c(τ ) are fully dependent on each other; thus, the latter can be obtained by the former from the inverse transformation (see Koutsoyiannis, 2013, for further details): Thus, the dependence structure of x(t) is represented either by the climacogram γ ( ) or the autocovariance function c(τ ).In addition, the Fourier transform of the latter, the spectral density function s(w), where w is the frequency, is of common use.Selection of an analytical model for c(τ ) or s(w) is usually based on the quality of fit in the range of observed (observable) values of τ and w which, for reasons mentioned above, does not include the "microscale" (τ → 0 or w → ∞) or in general the asymptotic behaviour.However, asymptotic stochastic properties of the processes are crucial for the quantification of future uncertainty, as well as for planning and design purposes (Montesarchio et al., 2009;Russo et al., 2006).Any model choice does, of course, imply an assumption about the nature of random variation asymptotically.Therefore, we may want this assumption (although fundamentally unverifiable) to be theoretically supported.
In this context, Koutsoyiannis (2011) connected statistical physics (the extremal entropy production concept, in particular) with stochastic representations of natural processes, which are otherwise solely data-driven.He demonstrated that extremization of entropy production of stochastic representations of natural systems, performed at asymptotic times (zero or infinity) results in the Hurst-Kolmogorov process (HKp), else known as fractional Gaussian noise (Mandelbrot and Van Ness, 1968).
HKp can be defined in continuous time by the following autocovariance function (Koutsoyiannis, 2013): which shows that autocovariance is a power function of lag τ ; consequently, it can be shown that the spectral density function s(w) is also a power law of the frequency w with exponent 1-2H .The three nominal parameters of the HKp are ν, α and H : the units of α and ν are [τ ] and [x] 2 , respectively, while H , the so-called Hurst coefficient, is dimensionless.Substituting Eq. ( 6) in ( 4) we obtain the climacogram of the process x Thus, the variance of x is a power law of the averaging time with exponent 2H − 2, precisely the same as that of c(τ ).
The climacogram contains the same information as the autocovariance function c(τ ) or the power spectrum s(w), because they are transformations of one another.Its relationship with the latter is given by (Koutsoyiannis, 2013) It has been observed that, when there is temporal dependence in the process of interest, the classical statistical estimation of the climacogram involves bias (Koutsoyiannis and Montanari, 2007), which is obviously transferred to transformations thereof, e.g.c(τ ) or s(w).The bias in the climacogram estimation can be determined analytically and included in the estimation itself (Koutsoyiannis, 2013).However, in the next section we show how the problems of uncertainty in statistical estimation may be extremely remarkable when using other uncontrollable quantities (e.g.high order moments) to justify or calibrate stochastic models.

Multifractal analysis
Multifractal analysis has been used in several fields in science to characterize various types of data sets, which have been investigated by means of the mathematical basis of multifractal theory.This is the basis for a series of calculations that reveal and explore the multiple scaling rules, if any, from data sets, in order to calibrate multifractal models.From a practical perspective, multifractal analysis is usually based upon the following steps (Lopes and Betrouni, 2009).
-Estimate the sample raw moments of different orders q over a range of aggregation scales .
-Plot the sample q-moments against the scale in a log-log diagram.
-Fit least-squares regression lines (one for each order q) through the data points.
-Estimate the multiscaling exponents K(q) as the slopes of regression lines (see Eq. 1).
The classical estimator of the qth raw moment of the local average process x High moments, i.e. q ≥ 3, mainly depend on the distribution tail of the process of interest.If we assume, for reasons mentioned in Sect. 1, scaling in state, i.e. a power-type (e.g.Pareto) tail, then raw moments are theoretically infinite beyond a certain order q max .However, their numerical estimates from a time series by Eq. ( 9) are always finite, thus resulting in infinite biases from a practical perspective, because the estimate is a finite number while the true value is infinity.Even below q max , where it can be proved that the estimates are unbiased, we show that the estimation of moments can be still problematic.It is easily shown, indeed, that the expected value of the moment estimator equals its theoretical value E[(x which can be used to derive the variance of the moment estimator as follows: This quantity can be assumed as a measure of uncertainty in the estimation of the qth moment of the local average process x ( ) j .Therefore, the estimator m ( ) q is theoretically unbiased (because of Eq. 10) but involves uncertainty (quantified by Eq. 11), which is expected to depend on statistical properties of the instantaneous process x(t) (i.e.marginal and joint distributional properties), the averaging scale , the sample size n, and the moment order q.

Estimation of the mean
The (unbiased) estimator of the common mean µ of the local average process x ( ) j is given by Eq. ( 9) for q =1, m ( ) where T is the largest timescale of averaging multiple of in a given observation period T o (Fig. 1).As a consequence of Eqs. ( 4) and ( 12), the variance of the estimator above can be expressed as follows: Therefore, the estimator m 1 is a function of the dependence structure of the continuous-time (instantaneous) process x(t), and the rounded observation period T .Note that the uncertainty in the estimation of the sample mean is independent of the timescale of averaging while it depends on the observation period T .
Considering now the HKp, the autocovariance function is given by Eq. ( 6).Hence, the climacogram γ (T ) takes the form of Eq. ( 7).In Fig. 2, we show how the temporal dependence (governed by the Hurst coefficient H for the HKp) influences the reliability of moment estimates.For simplicity and without loss of generality, we plot the ratio of Var[m to Var[x ( ) j ] for = 1 against the scale T , which equals the sample size n for = 1.As a consequence of Eqs. ( 13) and ( 7), the ratio is given by Var Notice that large values of H result in a much higher ratio than in the iid case (which is given by 1/n), and the convergence to the iid case is extremely slow (see Fig. 2).In essence, it can be argued that the greater the dependence in time, the harder it is to estimate the moment; in the sense that larger samples are required in order to obtain estimates of similar quality.

Estimation of higher moments
Let us now investigate the behaviour of estimators of higher order moments (q > 1) when the underlying random process exhibits dependence in time and when changing the process marginal distribution; this can be done by Monte Carlo simulation.Specifically, we use the Gaussian distribution and three one-sided distributions whose tails are sub-exponential, i.e. heavier than the former (as observed in several geophysical processes).All synthetic time series are generated in a way to have similar dependence structures based on the HKp, which are therefore governed by the Hurst coefficient H .In this study, we estimate the performance of qth moment estimators for four different common tail types (ordered from heavier to lighter): the Pareto, the log-normal, the Weibull and the Gaussian tails (see e.g.El Adlouni et al., 2008;and Papalexiou et al., 2013).The Pareto is the only power-type distribution, while the remaining three are of exponential type with all their moments finite.Specifically, we use the Pareto type II distribution, defined in [0, ∞), with survival function where β > 0 is the scale parameter, and κ > 0 the shape parameter.The latter, also known as the tail index, controls the asymptotic behaviour of the tail, which is given by x −1/κ ; as the value of κ increases the tail becomes heavier and consequently extreme values occur more frequently.Moreover, the shape parameter κ unequivocally defines the order q max =1/κ beyond which the qth moments are theoretically infinite, i.e.E[(x ( ) j ) q ] = ∞ for q ≥ 1/κ; in our study we assume κ = 0.2, and thus q max = 5.
The log-normal distribution, also defined in [0, ∞), is very commonly used in geosciences and has the survival function www.hydrol-earth-syst-sci.net/18/243/2014/ Hydrol.Earth Syst.Sci., 18, 243-255, 2014 where erfc( x exp −t 2 dt is the complementary error function, β > 0 is the scale parameter, and κ > 0 is the shape parameter that controls the behaviour of the tail (notice some differences from the more typical notational convention in the literature; see Forbes et al. (2011, p. 131) for further details).Despite all its moments being theoretically finite, the log-normal distribution is very similar in shape to a power-type distribution (Pareto), in the sense that the two distributions appear almost indistinguishable from each other for a large portion of their body (Mitzenmacher, 2004).Therefore, log-normal is regarded as a heavy-tailed distribution.
Another widely used distribution is the Weibull distribution, again defined in [0, ∞).Its survival function is a stretched exponential function (obtained by inserting a fractional power law into the exponential function), i.e.
where β > 0 is the scale parameter, and the stretching exponent 0 < κ < 1 (shape parameter) actually modifies the shape of the exponential distribution so as to obtain a heavier tail.
Consequently, the Weibull distribution can be regarded as a generalization of the exponential distribution, which is recovered with κ = 1.The case with κ > 1 (compressed exponential function, i.e. a tail lighter than the exponential one) has less practical importance, with the notable exception of κ = 2, which gives the Rayleigh distribution, closely related to the Gaussian distribution.

Monte Carlo simulation
As the log-normal model has been the most common in multifractal literature, we start our study from this model.For the Monte Carlo simulation we use the model introduced by Lombardo et al. ( 2012), which follows a disaggregation approach.In that respect it resembles the discrete multifractal cascade models, yet it is not affected by uncontrollable nonstationary issues that are typical in these multifractal cascades.The model starts the generation from the coarsest scale and then disaggregates into finer scales applying a specific scale-dependent exponential transformation to the HKp in a way to preserve part of its scaling properties.For the Monte Carlo experiment we generate 30 000 time series with sample size n = 2 10 = 1024, unit mean, standard deviation σ = 1.29 and H = 0.85.Later we will compare with the other models in a different setting, i.e. aggregating rather than disaggregating, using the same statistical properties (note that σ = 1.29 is the standard deviation of the Pareto type II with unit mean and tail index κ = 0.2).
The results of the Monte Carlo simulation experiment are depicted in Figs.3-6.Specifically, Fig. 3 shows the probability distribution of the natural logarithm of the ratio of qth moment estimates to their expected values (i.e. the theoretical values, following Eq.10).It can be noticed that the information content of the sample moments strongly decreases when increasing the order q (i.e. the distribution is less concentrated around 0): only low moments have reasonably low variation, all others vary within several orders of magnitude (notice that the horizontal axis is logarithmic and spans more than 10 orders of magnitude!).Despite the sample raw moment being an unbiased estimator of the true (population) raw moment, the probability distribution of the statistical estimator is very broad and skewed.This is particularly the case for high moments.Note that the averaging scale has negligible influence on the statistical characteristics of low moment estimators, while it slightly regularizes the behaviour of higher moment estimators.
In addition, in Fig. 4 we show the empirical frequency distribution of the sample fifth moment estimated from lognormal time series averaged locally over different timescales .Again here the bias is theoretically zero, but the most probable value of the moment estimate (the mode) is very different from its expected value.For example when = 1 (upper-left panel in Fig. 4), the mode of the distribution of m ( =1) 5 (green line) is almost two orders of magnitude less than the expected value (red line) and the probability of calculating from a unique sample a value equal to the mode is much greater (almost one order of magnitude) than the probability of obtaining the expected value itself.Recall that the expected value of the sample moment equals the true value of the moment, because of unbiasedness, but according to the distributions in Fig. 4 we can hardly expect the moment estimate from a unique sample to be close to this expected value.Increasing the averaging scale reduces the difference between the mean and the mode.Nonetheless, this difference is still remarkable at large scales (see e.g.lower-right panel in Fig. 4).
The large difference between the mode and the expected value of the moment estimators is not the only problem.Another problem is the high estimation uncertainty.In order to illustrate the uncertainty in the moment estimation, Fig. 5 shows semi-logarithmic plots of the prediction intervals of the sample moments, calculated from the Monte Carlo simulations, against the moment order, for various scales .The logarithmic scale on the vertical axis highlights the huge variability of estimates when the order increases.Note that the mean of raw moments (i.e. the true expected value) moves closer to the upper prediction limit for orders q > 3, thus making the use of high moments unreliable.Furthermore, Fig. 6 depicts log-log diagrams of the prediction intervals of the sample moments against the scale of averaging , for various orders q.In addition to the observations made with respect to Fig. 5, Fig. 6 shows that the increase of  the averaging scale has little influence on the variability of the moments, meaning that the sample size reduction is somewhat compensated by the time averaging.Nevertheless, it is clear that larger samples provide better estimates than smaller.For example, Meneveau and Sreenivasan (1991) propose a criterion of statistical convergence for the moments of local average processes, and find that data records of size 10 q may be sufficiently long to ensure statistical convergence for qth order moments.However, this is not immediately straightforward in case of highly correlated data series, as we show in Fig. 2. To further investigate this issue accounting for the criterion of convergence above, in Fig. 7 we show the trend of the interquartile range (IQR) of the prediction intervals for the third (q = 3) moment when increasing the sample size from 2 10 to 2 14 (the ensemble consists of 10 000 log-normal time series for each sample size generated by the model of Lombardo et al., 2012).It can be noticed that the sample size should be increased more than one order of magnitude to obtain roughly a 10 % improvement over the results presented in Fig. 5 for = 1.
In the second part of the Monte Carlo simulation experiment we use a different approach, first generating at the finest scale and then aggregating into coarser scales.In this case we generate 30 000 synthetic time series from the four distributions described in Sect.3.2 above (ordered from heavier to lighter tail type: Pareto, log-normal, Weibull with shape parameter smaller than 1 and Gaussian) with characteristics same as those in the previous experiment.In this case we Fig. 5. Semi-logarithmic plots of the prediction intervals of the sample moments versus the order q for various timescales , where "Q" stands for quantile.investigate how the classical estimators of raw moments behave when varying the tail type of the marginal distribution of the underlying stochastic process.To accomplish this aim, in Fig. 8 we plot on a semi-logarithmic scale the prediction intervals of the sample moments against the moment order (assuming = 1), for the four distributions.It can be seen that the tail type significantly influences the reliability of moment estimators.The heavier the distribution tail, the more uncertain the sample moments are.This is especially the case for high moments, because they depend enormously on the distribution tail and non-normality affects significantly their statistical properties.Analogous considerations apply to aggregated series (i.e.> 1).
It is emphasized that the vertical axes in Fig. 8 span more than 10 orders of magnitude yet the prediction limits do not necessarily bracket the true value of the moment.Particularly for the Pareto distribution the true (population) values of the fifth and sixth moments are infinite while their statistical estimates are finite and the entire graph does not provide any hint that these high moments differ so essentially from the lower ones.Another important conclusion drawn from Fig. 8 is that the prediction limits in the case of the Gaussian distribution are dramatically narrower than in all other cases.As the Gaussian distribution has been dominating in classical statistical applications and perhaps in statistical thinking, this fact may explain why the multifractal applications were Fig. 7. Semi-logarithmic plot of the interquartile range (IQR) (standardized with respect to the IQR for n = 2 10 ) of the prediction intervals for the third moment versus the sample size n for the lognormal series generated by our downscaling model (Lombardo et al., 2012).misled to neglect the huge uncertainty of high moment estimates and its impact on modelling.

Empirical moment scaling function
Since the ultimate aim of a multifractal analysis is to study the scaling of raw moments, we have carried out some additional numerical investigations on the generated samples by simply taking an average slope of linear regressions of sample moments at different scales in log-log diagrams (actually, this is commonly the case when dealing with real world data).Despite being not really crucial to the focus of our work (which aims to answer the question about how many raw moments we can estimate reliably), we believe it is worth exploring the variability in the estimates of the moment scaling function K(q), when using the statistical tools which we cautioned against.To accomplish this purpose, we use the log-normal synthetic series generated by the downscaling model of Lombardo et al. (2012).
In order to estimate an empirical exponent function K(q) describing the scaling of raw moments over a range of timescales, we should define the following non-dimensional quantities commonly used in the literature (e.g. de Lima and Grasman, 1999;Serinaldi, 2010).The scale ratio λ so that λ = 1 for the largest scale of interest max , i.e. λ = max / .In our case, we assume that max = [n/8] = 128, where the sample size n = 1024, so that sample moments can be estimated from at least eight data values, while the generic aggregated scale is bounded in [1,128].Similarly, we form the non-dimensional process ε(λ) dividing the local average of the continuous-time process x(t) by its mean at the largest scale max (or equivalently λ = 1); then where m is the temporal mean of the data series.The scaling behaviour of the process is characterized by the moment scaling function K(q) as follows: If K(q) linearly increases with q, then the process is said to be "simple scaling", otherwise it exhibits a "multiple scaling" behaviour.In Fig. 9, we graphically show how uncertainty in sample moments is reflected in the uncertainty in the estimates of scaling exponents.It can be noticed that the function K(q) shows a non-linear behaviour for the log-normal series, thus suggesting a multifractal behaviour.Analogous considerations apply to the series generated by the other Monte Carlo experiments described in Sect.3.3 above (not reported here).
The prediction intervals in Fig. 9 spread out widely while increasing the moment order q, which is consistent with an enhancement of uncertainty.We clarify that we used the ratios of moment estimates in all calculations to compute ε(λ).Nonetheless, recalling that we assumed the unit ensemble mean µ = E[x(t)] = 1 in all our Monte Carlo experiments, we found (not shown here) the same numerical results if using raw moments without taking any ratios.This is to stress that ratios of moments do not seem to play any significant role in the estimation of multiscaling exponents in our case.
It may be useful to add here some theoretical aspects.The theory of multifractals depends on the fact that raw moments obey power laws as the scale → 0 (or equivalently λ → ∞) (Falconer, 1990;Gneiting and Schlather, 2004), and so it depends on taking limits which cannot be achieved in reality.For most experimental purposes, the multifractal behaviour of a process x(t) is usually found by estimating the gradient of a graph of log(E[(ε(λ)) q ]) against log λ over an "appropriate" range of scales, where empirical points are closely matched by a straight line of slope K(q).Being the latter an asymptotic slope, it is difficult to find the "appropriate" range of scales to estimate K(q), because we could be misled by some artificial slopes which do not indicate the multifractal behaviour of the underlying process (see e.g.Koutsoyiannis, 2013).In addition, we should emphasize that the empirical moment scaling function K(q) varies across scales for ergodic processes.The simple proof for this is given in the Appendix A in the special case of q = 2.
Furthermore, we show in Appendix B that the theoretical moment scaling function K Th (q) for the model by Lombardo et al. (2012) is given by where H is the Hurst coefficient.Based on these findings, the empirical results in Fig. 9 do not seem to agree well with their theoretical counterparts.For example, in our case H = 0.85, for q = 4 the theoretical value should be K Th (q) =1.8, while the estimated mean value is about K(q) = 0.5 in the scale range of our Monte Carlo experiments.Hence, not finding  versus the order q for the log-normal series generated by the downscaling model (Lombardo et al., 2012).
the "appropriate" range of scales, in addition to estimation problems reported in our work, may lead to remarkable underestimation of the moment scaling function.

Conclusions
During recent decades, there has been a growing interest in multifractal analyses especially for the study of hydrological processes, particularly in rainfall modelling.Indeed, the multifractal framework provides parsimonious models to study the variability of several natural processes in geosciences, such as rainfall.Models following this approach require the scaling of the sample moments of different orders q, which is used in model identification and fitting.A common problem with the application of multifractal models, which in some cases may have led to incorrect results, is their disconnection from stochastic methodology and reasoning, and the (un-stated) naïve consideration that statistical estimates represent the true properties of a process.
Using theoretical reasoning and Monte Carlo simulations we find that the reliability of multifractal methods which use high order (> 3) is questionable.In particular, we highlight the problems in inference from time series of geophysical processes.The classical statistical approaches, often used in geophysical modelling, are based upon several simplifying assumptions, tacit or explicit, such as independence in time and exponential distribution tails, which are invalidated in natural processes.Indeed, the study of natural processes reveals scaling behaviours in state (departure from exponential distribution tails) and in time (departure from independence).While the multifractal models are based on these scaling behaviours per se, they may have failed to explore their statistical consequences with respect to the implied dramatic increase of uncertainty.
The following list briefly summarizes the main findings of our analyses.
-As natural processes are characterized by dependence in time, while classical statistics typically assumes independence, much larger samples are required in order to obtain estimates of similar reliability with classical statistics.
-Estimators of high moments whose distribution ranges over several orders of magnitude cannot support inference about a natural behaviour nor fitting of models.
-The most probable value of sample high moments (the mode) can strongly differ (by orders of magnitude) from its expected value (i.e. the true value), thus making the statistical estimate problematic even in the case of unbiasedness.
-The calculation of numerical values of high order moments is misleading as the theoretical moments may tend to infinity for high orders, while the sample estimates are always finite.Even smaller order moments can be very uncertain.
-Even if the generated process is multifractal, the sample estimates of the q moments from a unique sample can provide misleading results.
Hence, we have shown that distribution tails heavier than the exponential one and temporal dependence result in enormously increased uncertainty and/or infinite biases from a practical in raw moments.This paper warns practitioners against the blind use in geophysical time series analyses of classical statistical tools, which neglect dependence and heavy tails in distributions.Ossiander and Waymire (2000) already caution against using high moments in multifractal estimation, but their particular focus is on discrete multiplicative cascade models.Indeed, they demonstrate that the estimators of multiscaling exponents converge almost surely to the structure function of the cascade generators as the sample becomes large for all moment orders within a certain critical interval, whose upper bound is consistent with our results.Ignorance of increased uncertainty and inattentive use of high order moments may result in inappropriate modelling, wrong inferences and false claims about the properties of the processes.Evidently, the first two moments need to be used in all problems as they define the most important characteristics of the distribution, marginal (the first two moments) and joint (the second moment).Even for these two lowest moments it is important to always study their uncertainty and this only can be done in connection with a model fitted for the process of interest (as it is not possible to define uncertainty without specifying a model for the marginal distribution and dependence).The third moment is often useful as a measure of skewness but we should always be aware of its uncertainty; however using the third moment is not the only way to identify and assess the skewness of a distribution.For example, in parameter estimation of three-parameter distributions, it is better to avoid the method of moments and use other fitting methods such as maximum likelihood, L-moments, etc. Moments of order > 3 should be avoided in model identification and fitting because their estimation is problematic.If we have to use them, then it is imperative to specify their uncertainty and involve this uncertainty in any type of modelling and inference.

Fig. 1 .
Fig. 1.Sketch of the local average process x ( ) j obtained by averaging the continuous-time process x(t) locally over intervals of size .

)Fig. 2 .]
Fig. 2. Estimator variance of the mean of the local average process x ( =1) j standardized by the process variance, i.e.

Fig. 3 .
Fig. 3. Empirical cumulative distribution function (ecdf) of the natural logarithm of the ratio of qth moment estimates to their expected values E[(x ( )j) q ] = µ ( ) q when varying .

Fig. 4 .
Fig. 4. The epdf of the sample fifth moment estimated from log-normal time series averaged locally over different timescales .

Fig. 6 .
Fig. 6.Log-log plots of the prediction intervals of the sample moments versus the scale for various orders q.

Fig. 8 .
Fig.8.Semi-logarithmic plots of the prediction intervals of the sample moments versus the order q for various marginal probability distributions, assuming = 1.

Fig. 9 .
Fig. 9. Prediction intervals of the moment scaling function K(q)versus the order q for the log-normal series generated by the downscaling model(Lombardo et al., 2012).