Bayesian estimation of the self-similarity exponent of the Nile River ﬂuctuation

. The aim of this paper is to estimate the Hurst parameter of Fractional Gaussian Noise (FGN) using Bayesian inference. We propose an estimation technique that takes into account the full correlation structure of this process. Instead of using the integrated time series and then applying an estimator for its Hurst exponent, we propose to use the noise signal directly. As an application we analyze the time series of the Nile River, where we ﬁnd a posterior distribution which is compatible with previous ﬁndings. In addition, our technique provides natural error bars for the Hurst exponent


Introduction
Many geophysical systems exhibit non-trivial multi-scale correlation structures.In particular fractional Brownian motion and fractional Gaussian noise are often found to explain quite well the heterogeneity and multi-scale properties of geophysical time series in particular those from hydrology.In addition this kind of model is capable to discriminate between long and short term dependency.This difference has observable consequences since for instance a long memory process might explain the patterns observed in sedimentary deposits in river run-off areas (see e.g.Millen and Beard, 2003).For this reason it is important to devise suitable techniques to estimate the underlying self-similarity exponent, known as the Hurst exponent.Moreover it is important to provide methods to quantify the uncertainty of the so obtained measurements.We choose in this paper a Bayesian approach, which provides in a natural way both, a mean to Correspondence to: S. Benmehdi (sabah benmehdi@yahoo.fr)estimate the quantity of interest as well as a way to assess the uncertainty of its value.
Fractional Gaussian noise (FGN) is a Gaussian stochastic process {G H t ,t 0}, that can formally be viewed as the derivative of a fractional Brownian motion (FBM) B H t , t ≥ 0. This is the only Gaussian process with stationary increments which is self-similar and of zero mean.The parameter H characterizes its behavior under rescaling Here denotes equality of distributions.The first and second moments fully characterize this process The trajectories are continuous functions.However often time series of real data represent much noisier features, which however may have hidden self-similar correlation structures with sometimes long-range dependencies.See Fig. 7 for an example showing the level of Nile River from the years 622 to 1284 AD.For this consider the discrete process of increments over fixed time interval t = 1 The choice of time interval t = 1 is no loss of generality, since any other time interval would lead to the same process up to some multiplicative amplitude factor.This process is the coarse grained fractional Gaussian noise process FGN.
It is a stationary process, which is characterized by its auto correlation function.Using the correlation of the FBM we see that Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
For H = 1/2 we have the independent identically distributed Gaussian variables and uncorrelated time series after one single time step.For H = 1/2 however, as the time lag k gets large, the correlation function decays asymptotically like For this reason, H quantifies the correlation at large times.For H < 1/2 correlation is negative, whereas for H > 1/2 the correlation is positive.See Fig. 1 for an example of FGN.Various methods to estimate the self-similarity exponent H from a time series w n of observations have been proposed.A standard method consist in replacing w n by the associated random walk and to apply techniques to estimate the Hurst exponent of fractional Brownian motion.Here the following types of methods are proposed in the literature.The methods can be classified as temporal, spectral and time-scale methods.
The temporal methods are e.g.aggregated variance method (Taqqu et al., 1995), detrended fluctuation analysis (Goldberger et al., 2000;Peng et al., 1994), and range scaled analysis (Mandelbrot and Wallis, 1969).From the group of spectral methods, we would like to mention the log-periodogram method (Geweke and Porter-Hudak, 1983), the modified periodogram method, and the Whittle estimator (Whittle, 1953).The third group is in the wavelet domain, which includes the wavelet maximum likelihood WML estimator (McCoy and Walden, 1996) and the Abry-Veitch Daubechies wavelet-based estimator (Abry and Veitch, 1998).
In this paper we want to introduce a Bayesian estimator for the Hurst exponent H of FGN.This is an extension of the method proposed in (Makarava et al., 2011), where we introduced a Bayesian estimator for the Hurst exponent of fractional Brownian motion.One big advantage of our approach compared to others will be, that we take fully into account the correlation structure present in fractional Gaussian noise.Moreover, we will be able to produce error bars for all quantities, that we estimate.

Definition of the model
Often the data comes with an unknown offset.Instead of removing it prior to the analysis, we will incorporate it into the model.So we consider discrete time series of the form where G H n is the FGN and H ∈]0,1[ is the Hurst exponent, λ > 0 and λ ∈ R is the amplitude, and β ∈ R is the offset.Suppose we are given the observations X n , n = 1,...,N, the problem we want to address is how to obtain estimators for all involved parameters.For this we propose to use Bayes theorem, that reads and P(β,λ,H ) is some prior information we might have about the parameters.In the absence of such information, we suggest to use uninformative priors, P(β,λ,H ) λ −1 where as usual, we use the Jeffreys prior for the scale parameter λ.
For fixed parameters, the observations X = [X 1 ,...,X N ] t are multivariate Gaussian random variables with mean value and covariance given by with given by Eq. ( 1) and F t = [1,...,1] is the vector with N components, where [.] t denotes the transpose.The likelihood function for the parameters λ,β,H can now be written as From Bayes theorem, we now obtain the following posterior density of our parameters The normalization constant C is chosen to make the posterior probability density integrate to one.Then the numerator of the exponent can be written as the following quadratic polynomial in β  From this we see, that the posterior can be written in the following "Gaussian" form Here f µ,σ 2 is the one dimensional Gaussian probability density function with mean µ and variance σ 2 .Note that H and therefore, the residuum R 2 and the mode β * depend on Equation ( 6) describes the full posterior information that we have about all the parameters jointly.In case, we are interested in only one of them, we may treat the other parameters as completely unknown and integrate them out, to obtain the marginal distributions.Integrating over λ we obtain

P(β,H
The integral over the offset yields The integral over H reads dH. Integrating Eq. ( 12) over λ yields the posterior distribution of H that we can infer from the observations In practice, this integral has to be performed numerically.
Taking the position of the maximum of this posterior density, we obtain the maximum posterior estimate of the scaling exponent H .In the same way we may obtain posterior estimates of λ and of β and From these expressions we also may produce point estimators.For instance we may set This in turn yields an estimator for the offset β = β * ( Ĥ ). (17)

Application to synthetic data
In this section, we want to show, how to apply the analysis to synthetic data.In the first part, these are data of fractional Gaussian noise.In the second part, we will investigate synthetically generated Rosenblatt processes as examples of non-Gaussian H -self-similar processes.

Fractional Gaussian noise
For the investigation of fractional Gaussian noise, we generate a random sample Several methods have been proposed to numerically generate such a realization of FGN, and we refer for more details to (Dieker, 2004;Kang, 2008) and the reference therein.Our simulations are based on the use of the Cholesky decomposition of H = L t L where we then apply L t e, where e = [e 1 ,...,e N ] t is a random vector of N independent standard Gaussian random variables.In Fig. 2 we have shown this data.
We then have computed the Bayesian posterior distribution of λ and H together with one of their marginal distributions (see Fig. 6a).As you can see, posterior information is well localized around the data generating parameters.However, since the posterior itself is computed from a random realization its maximum does not coincide with the true value but randomly fluctuates around it.
Next, we use the wavelet based joint estimator by Veitch and Abry (1999) and the periodogram estimator by Robinson (1995) to compare with the point estimator we propose in this work (Eq.16).We perform 20 000 realizations by Monte-Carlo simulations for the fixed Hurst exponent H = 0.7 and with length N = 128.For each of them we produce the maximum posterior estimator Ĥ of our method, wavelet based joint estimator and periodogram method estimator.We implement the wavelet based joint estimator as is proposed by Veitch and Abry (2002).Figure 3 shows the comparison between the methods and in Table 1 the intervals containing ≥ 90 % of the distribution are presented.It is clearly depicted that the Bayesian method outperforms the wavelet and periodogram methods.
To make the validation test of the proposed method we quantify the bias of the maximum posterior estimator Ĥ , E( Ĥ ) − H, as a function of the number of data points.For that we generate 500 realizations of fractional Gaussian noise starting from one single observation point.In Fig. 4 it is shown that even starting from 20 data points, the bias quickly decays.

Rosenblatt process
In this section the Rosenblatt process as an example of a Hself-similar process, whose finite-dimensional distributions are non-Gaussian, is considered.The Rosenblatt process (see Rosenblatt, 1961) is the H -self-similar process with Hurst exponent H ∈ (1/2,1) and stationary increments.It can be written in explicit form as a stochastic integral where B(y),y ∈ R is pure Brownian motion, a(H ) is a positive normalization constant chosen such that E(Z(1) 2 ) = 1, and x + = max{x,0}.We implement Matlab code from Abry and Pipiras (2006) for the wavelet-based synthesis of the Rosenblatt process and derive the increments of the generated processes.In Fig. 5 we show the example of the increment process for the Rosenblatt process with Hurst exponent H = 0.7 of the length N = 129.The obtained averaged posterior densities for point estimator derived with Bayesian method for 1000 realizations is also shown.
This data set was studied before by several authors to estimate the Hurst exponent.For example in Liu et al. (2009) the modified Periodogram Estimators (MPE) is used and a comparison with the following methods is presented: the Periodogram Estimator (PE), the Whittle estimator, the wavelet maximum likelihood (WML) estimator and estimators based on the associated FBM.
We have computed the posterior distribution according to Eq. ( 6).For reasons of visualization, we have computed the marginals over H , λ and β, which are then functions of the remaining two variables only.The results are visible in Fig. 6b-d).We thus find the following posterior estimates (maximum of posterior) of our parameters together with their 90 % confidence intervals (Table 2).
In Table 3 we give a summary of the known results from (Liu et al., 2009), and how they compare to this study.We see that the values obtained by the other methods are in the confidence interval of our result, except for the result obtained by the periodogram estimator (PE).In addition, our method provides estimations of all parameters involved in the model and we also obtain error bars for each of them.

Conclusions
In this paper we have proposed a Bayesian estimation technique of the Hurst parameter for the fractional Gaussian noise process.We have considered a slightly more general model, where in addition we have taken the offset and the amplitude as parameters.This technique yields, besides the point estimations of the model's parameters, confidence intervals that enclose a given percentage of the posterior uncertainties.The method was tested successfully on synthetic realizations of fractional Gaussian noise processes.In addition, we performed tests on synthetic realizations of the Rosenblatt process as an example for a non-Gaussian H -self-similar process.Our method turned out to give good results also for Rosenblatt processes.But we can not generalize this to any non-Gaussian self-similar processes.When applied to the historical time series of the Nile River, our results are compatible with previous findings.However in addition we are able to provide error bars for all estimates.
Edited by: J. Kurths Reviewed by: two anonymous referees

Fig. 1 .
Fig. 1. (Top) Simulation of fractional Gaussian noise for N = 200 and H = 0.05.(Botton) Simulation of fractional Gaussian noise for N = 200 and H = 0.95.The difference in correlation structure is clearly visible.

Fig. 3 .
Fig. 3.The distribution of the estimator of the Hurst exponent for H = 0.7,N = 128 with 20 000 realizations by wavelet based estimator (solid line), periodogram method (dotted line) and Bayestian method (dashed line).

Fig. 4 .
Fig. 4. The validation test for Bayesian estimation with 500 realizations for the Hurst exponent H = 0.7.

Fig. 5 .
Fig. 5. (Top) Simulation of Rosenblatt process for N = 129 and H = 0.7.(Bottom) The averaged of the obtained posterior densities for point estimator derived with Bayesian method for 1000 such realizations.

Fig. 6 .
Fig. 6.Normalized two dimensional marginal posterior densities.The maxima indicate the most likely estimates in: (a) H − λ plain for signal X n = 1.5G 0.3 n + 8.0, n = 1,...,N = 100.Note that the posterior is well localized around the point corresponding to the true value (white spot); (b) H − β plain for Nile River; (c) λ − β plain for Nile River; (d) H − λ plain for Nile River; on the axis, the one dimensional projections of the posterior densities are depicted.The white contour-line encloses 90 % of posterior probability.It therefore quantifies the posterior uncertainty of the parameters together with their posterior dependency.

Fig. 7 .
Fig. 7. (Top) The time series of minimal water levels of Nile River near Cairo.(Bottom) The integrated time series.

Table 1 .
The interval for Ĥ with p ≥ 90% for synthetic data with H = 0.7 and N = 128.

Table 2 .
Confidence intervals on the estimated parameters.

Table 3 .
Estimation of the Hurst exponent for different methods.
The interval containing ≥ 90 % of the distribution here is [0.596,0.838].Benmehdi et al.: Bayesian estimation of the self-similarity exponent of the Nile River fluctuation Gauge near Cairo.The data set we analyse is publicly available at StatLib archive: http://lib.stat.cmu.edu/S/beran.In Fig.7we have plotted the time series together with its first discrete integral that is defined by www.nonlin-processes-geophys.net/18/441/2011/ Nonlin.Processes Geophys., 18, 441-446, 2011S.