A new stochastic process to model Heart Rate series during exhaustive run and an estimator of its fractality parameter

A ne w stochas tic proces s to model Hear t Rate series during e xhaustiv e run and an estima tor of its fractality parameter Jean-Marc Bardet Université P aris 1, SAMOS-MA TISSE-CES, 90 rue de T olbiac, 75013 P ari s, F rance. Véronique Bill at Université d’Evry , LEPHE, E.A. 3872 Genop ole, Boule vard F . Mitterrand, Evry Cedex, Fr a nce. Imen Kammoun Université P aris 1, SAMOS-MA TISSE-CES, 90 rue de T olbiac, 75013 P ari s, F rance. Summary . In order to interpr et and explain the physiological signal behaviors, it can be in- teresting to find some constants among the fluctuati ons of these da ta durin g all the effor t or duri ng different stages of the race (wh ich can be detecte d using a change poin ts d etection method). Sev eral recent papers hav e proposed the long -range depende nce (Hurst) parameter as such a constan t. Howe ver , the ir resul ts i nduce two main problems. F irstly , DF A method is usually appli ed for estimating thi s paramete r . Clearly , su ch a method d oes not provide the most efficie nt estimator and moreover it is not at all robust ev en i n the ca se of smooth tren ds. Secondly , this method ofte n gives e stimated Hurst parameters larger th an 1 , which is the larger possible v alue for long me mor y statio nar y processes. In this ar ticle we pro pose solutio ns for both these problems and we define a new model allowing such estimated paramete rs. On the one h and, a wav e let-ba se estimator is ap plie d to data. Such an estimator p rovides opti mal conv e rgence rates in a semiparametr ic context an d can b e used f o r smoothly t rended pro- cesses. On the oth er hand, a new semiparametric model so-called locally fractional Gaussian noise i s introduced a nd is characterized by a so-called parameter wh ich can be la rger t han 1 . Such semiparametr ic process is tested to be relev a nt for modeling HR data in the t hree characteristi c phase s of the race. It al so sh ows an ev ol ution o f the local fractality pa rameter duri ng the race confir ming the re sults obta ined by P e ng et a l. (1995) i n thei r stud y rega rding Hurst parameter of HR time series dur ing the ex ercise for healthy a dults (where the estimated parameter is close to t hat ob ser ved in the race begin ning) and h ear t failure ad ults (wher e the estimated parameter is close to that obser ved in the e nd of race). So , this ev oluti on, whi ch can not b e o bser ved with DF A method, may be associated wit h fatigue app ear ing d uri ng the last phase of the marathon. K eywords: W av e let an alysis; Detrended fluctuatio n a nalysis; F ractional Gaussian noise; Self- similar ity; Hurst parameter; Long -range depen dence processes; Hear t r ate time seri es 1. Introduction The on ten t of this artile w as motiv ated b y a general study of ph ysiologial signals of run- ners reorded during endurane raes as marathons. More preisely , after dieren t signal pro edures for "leaning" data, one onsiders the time series resulting of the ev olution of E-mail: bardet @univ-pa ris1.f r 2 Bardet et al. heart rate (HR) data during the rae. The follo wing gure pro vides sev eral examples of su h data (reorded during Marathon of P aris). F or ea h runner, the p erio ds (in ms) b et w een the suessiv e pulsations (see Fig. 1 ) are reorded. The HR signal in n um b er of b eats p er min ute (bpm) is then dedued (the HR a v erage for the whole sample is of 162 bpm). This pap er, fo uses on the mo deling and the estimation of relev an t parameters  hara- terizing these instan taneous heart rate signals of athletes reorded during the marathon. W e ha v e  hosen to fo us in an exp onen t that an b e alled "F ratal", whi h indiates the lo al regularit y of the path and the dep endeny b et w een data. In ertain stationary ases, this parameter is lose to the Hurst parameter, dened for long range dep enden t (LRD) pro esses. 0.5 1 1.5 2 x 10 4 320 340 360 380 400 420 HR(Msec.) Ath.1 0.5 1 1.5 2 x 10 4 2.2 2.4 2.6 2.8 3 HR(Hertz) Ath.1 0.5 1 1.5 2 x 10 4 140 150 160 170 180 HR(BPM) Ath.1 0.5 1 1.5 2 2.5 x 10 4 140 150 160 170 180 HR(BPM) Ath.2 0.5 1 1.5 2 2.5 x 10 4 110 120 130 140 150 160 170 HR(BPM) Ath.3 0.5 1 1.5 2 2.5 x 10 4 130 140 150 160 170 180 HR(BPM) Ath.4 Fig. 1. Heat rate sig nals of At hlete 1 in ms, Her tz and BPM (u p), of Athlete s 2, 3 an d 4 in BPM (down) The LRD b eha vior is often seen on v arious data. This phenomenon w as dev elop ed in man y elds b eginning with h ydrology (Hurst, 1951), teleomm uniation, biome hanis and reen tly in eonom y and nane. A mathematial and signal pro essing metho ds ha v e also noted the presene of LRD in time series desribing the utuations o v er time of ph ysiologial signals (Goldb erger (2001 ), Goldb erger et al. (2002 )). Indeed, n umerous authors ha v e studied heartb eat time series (see for instane P eng et al. (1993 ), (1995) or Absil et al. (1999 )) and the mo del prop osed to t these data is a trended long memory pro ess with an estimated Hurst parameter lose to 1 (and sometimes more than 1 ). In this artile, three impro v emen ts ha v e b een prop osed to su h a mo del: (a) Data are stepp ed in three dieren t stages whi h are deteted using a  hange p oin ts detetion metho d (see for instane La vielle (1999 )) dev elop ed in Setion 2. The main idea of the detetion's metho d is to onsider that the signal distribution dep ends on a v etor of unkno wn  harateristi parameters onstituted b y the mean and the v ariane. The dieren t stages (b eginning, middle and end of the rae) and therefore the dieren t v etors of parameters, whi h  hange at t w o unkno wn instan ts, are estimated. This rst step is imp ortan t sine the mo del dened b elo w ts w ell for ea h sub-series but not at all for the whole HR data. (b) The parameter H whi h is v ery in teresting for in terpreting and explaining the ph ysio- logial signal b eha viors, is estimating using t w o metho ds the DF A metho d and w a v elet A new stochastic process to m odel HR seri es and an estimator of its fractality parameter 3 analysis in Setion 3. The DF A whi h is a v ersion, for time series with trend, of the aggregated v ariane metho d w as studied in details in Bardet and Kammoun (2008 ). In partiular, the asymptoti prop erties of the DF A funtion and the dedued estimator of H are stud- ied in the ase of frational Gaussian noise and extended to a general lass of sta- tionary semiparametri long-range dep enden t pro esses with or without trend. W e ha v e sho wn that DF A metho d is not as eien t to estimate the Hurst parameter of stationary long memory pro esses, than other metho ds su h as log-p erio dogram (see for instane, Moulines and Soulier (2003) or w a v elet analysis (see Abry et al. (1998), or Moulines et al. (2007)), whi h pro vide the optimal on v ergene rate (in sense of the minimax riterium). Moreo v er, if the DF A metho d is not at all robust in the ase of p olynomial trends, this is not su h a ase of the w a v elet analysis metho d. Finally , a go o dness-of-t test an b e dedue from the w a v elet analysis metho d. Despite the p opularit y of DF A metho d in n umerous pap ers onerning su h ph ysiologial signals (see for instane, Absil et al. (1999), Iv ano v et al. (2001), P eng et al. (1993), (1995)), it is therefore learly more in teresting to use the w a v elet based estimator in view of estimate a fratal exp onen t of HR data. As a rst mo del, the usual frational Gaussian noise (F GN) is then prop osed for mo deling HR data. In su h a on text, the w a v elet based estimator pro vides t w o re- sults. Firstly , the estimated parameter often exeeds the v alue 1 , whi h is the largest p ossible v alue for a F GN. Seondly , ev en for the 3 dieren t stages of the rae, the go o dness-of-t test is alw a ys rejeted. () In Setion 4, w e prop ose to mo del HR data, during ea h stage, with a generalization of frational Gaussian noise, alled lo ally frational Gaussian noise. Su h stationary pro ess is built from a parameter alled lo al fratalit y parameter whi h is a kind of Hurst parameter in a restrited band frequeny (that ma y tak e v alues in R and not only in (0 , 1) as usual Hurst parameter). The estimation of lo al fratalit y parameter and also the onstrution of go o dness-of-t test an b e made with w a v elet analysis. W e also sho w the relev ane of mo del and an ev olution of the parameter during the rae, whi h onrms results obtained b y other authors in their study regarding the distinguish of health y from pathologi data (see P eng et al. (1995)). 2. Abrupt change detection During eort, one or more phases an b e observ ed in reorded HR series, whi h ev olv e and  hange dieren tly from an athlete to another: the transition step, reorded b et w een the rae b eginning and the stage of HR rea hed during the eort, the main stage during the exerise and an arriv al phase un til the rae end. So, after leaning the HR data (implying that only 9 athlete HR data are no w onsidered) an automati detetion of  hanges is applied to HR time series utting its in dieren t rae phases - b eginning, middle and end. The  hange p oin t detetion metho d used here is dev elop ed b y La vielle (see for instane La vielle (1999 )). The main idea is to onsider that the signal distribution dep ends on a v etor of unkno wn  harateristi parameters in ea h stage. The dieren t stages and there- fore the dieren t v etors of parameters,  hange at unkno wn instan ts. F or instane and it will b e our  hoie,  hanges in mean and v ariane an b e deteted. Applied to the data, the  hange p oin t detetion metho d along these t w o phenomena distinguishes b eginning, middle 4 Bardet et al. and end of rae. So, it ma y b e p ossible to en visage a pieewise stationarit y i.e. that the pro ess is almost stationary on xed time in terv als and it remains the mo deling of this stationary omp onen t. General principle of the method of change detection Assume that a sample of a time series ( Y ( i ) , i = 1 , . . . , n ) is observ ed. Assume also that it exists τ = ( τ 1 , τ 2 , . . . , τ K − 1 ) with 0 = τ 0 < τ 1 < τ 2 < ... < τ K − 1 < n = τ K and su h that for ea h j ∈ { 1 , 2 , . . . , K } , the distribution la w of Y ( i ) is dep ending on a parameter θ j ∈ Θ ⊂ R d (with d ∈ N ) for all τ j − 1 < i ≤ τ j . Therefore, K is the n um b er of segmen ts to b e dedued starting from the series and τ = ( τ 1 , τ 2 , . . . , τ K − 1 ) is the ordered  hange instan ts. No w, dene a on trast funtion U θ  Y ( τ j + 1) , Y ( τ j + 2) , . . . , Y ( τ j +1 )  , of θ ∈ R d applied on ea h v etor  Y ( τ j + 1) , Y ( τ j + 2) , . . . , Y ( τ j +1 )  for all j ∈ { 0 , 2 , . . . , K − 1 } . A general example of su h a on trast funtion is U θ  Y ( τ j + 1) , Y ( τ j + 2) , . . . , Y ( τ j +1 )  = − 2 log L θ  Y ( τ j + 1) , Y ( τ j + 2) , . . . , Y ( τ j +1 )  , where L θ is the lik eliho o d. Then, for all j ∈ { 0 , 2 , . . . , K − 1 } , dene: b θ j = Argmin θ ∈ Θ U θ  Y ( τ j + 1) , Y ( τ j + 2) , . . . , Y ( τ j +1 )  . No w, set: b G ( τ 1 , . . . , τ K − 1 ) = K − 1 X j =0 U b θ j  Y ( τ j + 1) , Y ( τ j + 2) , . . . , Y ( τ j +1 )  As a onsequene, an estimator ( b τ 1 , . . . , b τ K − 1 ) an b e dened as: ( b τ 1 , . . . , b τ K − 1 ) = Argmin 0 <τ 1 <τ 2 <...<τ K − 1 0 . As a onsequene, b y minimizing V in τ 1 , . . . , τ K − 1 , K , an estimator b K is obtained whi h v aries with the p enalization parameter β . F or HR data, the  hoie of p en ( K ) w as K . Let b G K = b G ( b τ 1 , . . . , b τ K − 1 ) , for K = K 1 , . . . , K M AX w e dene β i = b G K i − b G K i +1 K i +1 − K i and l i = β i − β i +1 with i ≥ 1 . Then the retained K is the greatest v alue of K i su h that l i >> l j for j > i . Applied to the whole set HR data, the n um b er of abrupt  hanges is estimated at 4 or 3. Three phases w ere seleted to b e studied, whi h are lo ated in the b eginning of the rae, in the middle and in the end (see for example Fig. 2). Ho w ev er for ertain reorded signals the rst or the last phase an not b e distinguished probably for measuremen ts reasons. 0 0.5 1 1.5 2 2.5 x 10 4 130 140 150 160 170 180 190 200 Fig. 2. The estimated configuratio n of change s in a HR time ser ies of an athlet e 6 Bardet et al. In order to un v eil if a  hange of b eha vior of HR series w as happ ened during these three deteted phases of the marathon, w e prop ose a new mo del for these sub-series  haraterized b y a parameter H . Sev eral ommon estimators of this parameter, so-alled saling b eha vior exp onen ts, onsist in p erforming a linear regression t of a sale-dep enden t quan tit y v ersus the sale in a logarithmi represen tation. This inludes the Detrended Flutuation Analysis (DF A) metho d P eng et al. (1994 ) and the w a v elet analysis metho d Abry et al. (2003 ). 3. A fractional Gaussian noise for modeling HR series and the estimation of the Hurst parameter In this setion, a rst mo del, the frational Gaussian noise, is prop osed for mo deling HR data. After a statisti study , one  ho oses to estimate the Hurst parameter with a w a v elet based estimator instead of the DF A metho d (whi h is ommonly used in ph ysiologi pap ers despite its w eak p erformanes). Moreo v er, a test built from the w a v elet based metho d sho ws the badness-of-t of this mo del to the data. 3.1. A first model: the fractional Gaussian noi se When w e observ e en tire or partial (during the three phases) HR time series, w e remark that it exhibits a ertain p ersistene and the related orrelations dea ys v ery slo wly with time what  haraterizes tra jetories of a long memory Gaussian noise. Also, the distribution of data reorded during the phases leads as to susp et a Gaussian b eha vior in these data. Of ourse this is only an assumption and w e an  he k it with tests onsidered for long range dep enden t pro esses. But in our ase w e will try to test whether a Gaussian pro ess ould mo del these data. Moreo v er, the aggregated signals (see for example Fig. 3) presen t a ertain regularit y v ery lose to that of frational Bro wnian motion sim ulated tra jetories with a parameter lose to 1 (Fig. 5 ). So, frational Gaussian noise ould b e an appropriated mo del to HR series. 0 0.5 1 1.5 2 2.5 x 10 4 −14000 −12000 −10000 −8000 −6000 −4000 −2000 0 2000 4000 6000 1 1.1 1.2 1.3 1.4 1.5 x 10 4 −13000 −12000 −11000 −10000 −9000 −8000 −7000 −6000 −5000 −4000 −3000 1.2 1.25 1.3 x 10 4 −1.25 −1.2 −1.15 −1.1 −1.05 x 10 4 Fig. 3. T he self-similar ity of the aggregate d HR sign als (representation of the aggregate d HR fluctu- ations at 3 different time resolu tions) A new stochastic process to m odel HR seri es and an estimator of its fractality parameter 7 The follo wing Figure 4 presen ts a omparison b et w een the graphs of HR data during a stage (deteted previously) and a frational Gaussian noise (F GN in the sequel) with parameter H = 0 . 99 (see the denition ab o v e). Before using statistial to ols for testing the similarities of b oth these graphs, let us remind some elemen ts onerning the F GN. 0 2000 4000 6000 8000 10000 12000 14000 16000 −2 −1 0 1 2 3 x 10 −5 FGN 0 2000 4000 6000 8000 10000 12000 14000 16000 150 155 160 165 170 175 180 HR (BPM) Fig. 4. Compa rison of HR data in the middle of r ace (Ath4) and generated FGN(H=0.99) trajector ies The F GN is one of the most famous example of stationary long range dep enden t (LRD in the sequel) pro ess. The LRD phenomenon w as observ ed in man y elds inluding teleomm uniation, h ydrology , biome hani, eonom y ... A stationary seond order pro ess Y = { Y ( k ) , k ∈ N } is said to b e a LRD pro ess if: X k ∈ N | r Y ( k ) | = ∞ with r Y ( k ) = E  Y (0) Y ( k )  . Th us Y ( k ) is dep ending on Y (0) ev en if k is a v ery large lag. Another w a y for writing the LRD prop ert y is the follo wing: r Y ( k ) ∼ k 2 H − 2 L ( k ) , as k → ∞ , with L ( k ) a slo wly v arying funtion ( i.e. ∀ t > 0 , L ( xt ) /L ( x ) → 1 when x → ∞ ) and the Hurst parameter H ∈ ( 1 2 , 1 ) . The LRD is losely related to the self-similarit y onept. A pro ess X = { X ( t ) , t ≥ 0 } is so alled a self-similar pro ess with self-similarit y exp onen t H , if ∀ c > 0 :  X ( ct )  t L = c H  X ( t )  t . No w, if w e onsider the aggregated pro ess { X ( t ) , t ≥ 0 } dened b y X ( k ) = P k i =1 Y ( i ) with a LRD pro ess Y , then under w eak onditions (for instane Y is a Gaussian or a ausal linear pro ess), it an b e pro v ed that, roughly sp eaking, for k → ∞ , the distribution of { X ( t ) , t ≥ k } is a self-similar distribution (see Doukhan et al. , 2003, for more details). The F GN is an example of a LRD Gaussian pro ess. More preisely , Y H = { Y H ( k ) , k ∈ N } is a F GN, when r Y H ( k ) = σ 2 2 ( | k + 1 | 2 H − 2 | k | 2 H + | k − 1 | 2 H ) ∀ k ∈ N , 8 Bardet et al. with H ∈ (0 , 1) and σ 2 > 0 . As a onsequene, for H ∈ ( 1 2 , 1 ) , a usual T a ylor form ula implies r Y H ( k ) ∼ σ 2 H (2 H − 1 ) k 2 H − 2 , when k → ∞ . F or a zero-mean F GN, the orresp onding aggregated pro ess, denoted here X H , is so-alled the frational Bro wnian motion (FBM) and X H is a self-similar Gaussian pro ess with self-similar parameter H and therefore satises, V ar ( X H ( k )) = σ 2 | k | 2 H ∀ k ∈ N (it an b e ev en pro v ed that X H is the only Gaussian self-similar pro ess with stationary inremen ts). It is ob vious that Y H ( k ) = X H ( k ) − X H ( k − 1 ) , the sequene of the inremen ts of a FBM, is a F GN. 0 500 1000 −1 −0.5 0 0.5 1 FGN 0 500 1000 −3 −2 −1 0 1 2 FBM H=0.2 0 500 1000 −0.1 −0.05 0 0.05 0.1 0.15 0 500 1000 0 0.5 1 1.5 2 H=0.5 0 500 1000 −0.01 −0.005 0 0.005 0.01 0.015 0 500 1000 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 H=0.8 Fig. 5. Generated F GN trajecto rie s and correspond ing aggregate d series (FBM) for H = 0 . 2 < 0 . 5 anti-p ersistant noise (left), H = 0 . 5 white noise (center) and H = 0 . 8 > 0 . 5 LRD process (r ight) Sev eral generated tra jetories of F GN and orresp onding FBM are presen ted in Fig. 5 for dieren t v alues of H . F or testing if a HR path an b e suitably mo del b y a F GN, a rst step onsists in estimating H . Here w e  hose to use t w o estimators (but there exist man y else, see for instane Doukhan et al. , 2003) that are kno wn to b e un hanged to the presene of a p ossible trend. 3.2. T wo estimators of the Hurst parameter: DF A and wa velet based estimators F or estimating H , a frequen tly used metho d in the ase of ph ysiologial data pro essing is the Detrended Flutuation Analysis (DF A). The DF A metho d w as in tro dued b y P eng et al. P eng et al. (1994 ). The DF A metho d is a v ersion for trended time series of the metho d of aggregated v ariane used for long-memory stationary pro ess. It onsists briey on: (a) Aggregated the pro ess and divided it in to windo ws with xed length, (b) Detrended the pro ess from a linear regression in ea h windo ws, () Computed the standard deviation of the residual errors (the DF A funtion) for all data, A new stochastic process to m odel HR seri es and an estimator of its fractality parameter 9 (d) Estimated the o eien t of the p o w er la w from a log-log regression of the DF A funtion on the length of the  hosen windo w (see Fig. 6). After the rst stage, the pro ess is supp osed to b eha v e lik e a self-similar pro ess with stationary inremen ts added with a trend (see previously). The seond stage is supp osed to remo v e the trend. Finally , the third and fourth stages are the same than those of the aggregated metho d (for zero-mean stationary pro ess). An example of the DF A metho d applied to a path of a F GN with dieren t v alues of H is sho wn in Fig. 6. −3 −2 −1 0 1 2 −14 −12 −10 −8 −6 −4 −2 0 2 4 log. of chosen frequencies (Hz) log. of IN 1.5 2 2.5 3 3.5 −3 −2.5 −2 −1.5 −1 −0.5 log10(ni) log10(F(ni)) Fig. 6. Results of the DF A metho d a nd wav e let anal ysis ap plie d to a pa th of a discretized FGN for different v alues of H = 0 . 2 , 0 . 4 , 0 . 5 , 0 . 7 , 0 . 8 , with N = 10000 In Bardet and Kammoun (2008 ), the asymptoti prop erties of the DF A funtion in ase of a F GN path ( Y (1 ) , . . . , Y ( N )) are studied. In su h a ase the estimator b H DF A on v erges to H with a non-optimal on v ergene rate ( N 1 / 3 instead of N 1 / 2 rea hed for instane b y max- im um lik eliho o d estimator). An extension of these results for a general lass of stationary Gaussian LRD pro esses is also established. In this semiparametri frame, w e ha v e sho wn that the estimator b H DF A on v erges to H with an optimal on v ergene rate (follo wing the minimax riteria) when an optimal length of windo ws is kno wn. The pro essing of exp erimen tal data, and in partiular ph ysiologial data, exhibits a ma jor problem that is the nonstationarit y of the signal. Hu et al. (2001) ha v e studied dieren t t yp es of nonstationarit y asso iated with examples of trends and dedued their eet on an added noise and the kind of omp etition who exists b et w een this t w o signals. They ha v e also explained (2002) the eets of three other t yp es of nonstationarit y , whi h are often enoun tered in real data. In Bardet and Kammoun (2008 ), w e pro v ed that b H DF A do es not on v erge to H when a p olynomial trend (with degree greater or equal to 1 ) or a piee- wise onstan t trend is added to a LRD pro ess: the DF A metho d is learly a non robust estimation of the Hurst parameter in ase of trend. F or impro ving this estimation at least for p olynomial trended LRD pro ess, a w a v elet based estimator is no w onsidered. This metho d has b een in tro dued b y Flandrin (1992) and w as 10 Bardet et al. dev elop ed b y Abry et al. (2003) and Bardet et al. (2000). In W esfreid et al. (2005), a m ultifratal analysis of HR time series is presen ted for trying to un v eil their saling la w b eha vior using the W a v elet T ransform Mo dulus Maxima (WTMM) metho d. Let ψ : R → R a funtion so-alled the mother w a v elet. Let ( a, b ) ∈ R ∗ + × R and de- note λ = ( a, b ) . Then dene the family of funtions ψ λ b y ψ λ ( t ) = 1 √ a ψ  t a − b  P arameters a and b are so-alled the sale and the shift of the w a v elet transform. Let us underline that w e onsider a on tin uous w a v elet transform. Let d Z ( a, b ) b e the w a v elet o eien t of the pro ess Z = { Z ( t ) , t ∈ R } for the sale a and the shift b , with d Z ( a, b ) = 1 √ a Z R ψ ( t a − b ) Z ( t ) dt = < ψ λ , Z > L 2 ( R ) . F or a time series instead of a on tin uous time pro ess, a Riemann sum an replae the previous in tegral for pro viding a disretized w a v elet o eien t e Z ( a, b ) . The funtion ψ is supp osed to b e a funtion su h that it exists M ∈ N ∗ satisfying , Z R t m ψ ( t ) dt = 0 for all m ∈ { 0 , 1 , . . . , M } . (2) Therefore, ψ has its M rst v anishing momen ts. The w a v elet based metho d an b e applied to LRD or self-similar pro esses for resp etiv ely estimating the Hurst or the self-similarit y parameter. This metho d is based on the follo wing prop erties: for Z a stationary LRD pro ess or a self-similar pro ess ha ving stationary inremen ts, for all a > 0 , ( d Z ( a, b )) b ∈ R is a zero-mean stationary pro ess and • If Z is a stationary LRD pro ess, E ( d 2 Z ( a, b )) = V ar ( d Z ( a, b )) ∼ C ( ψ, H ) a 2 H − 1 when a → ∞ • If Z is a self-similar pro ess ha ving stationary inremen ts, E ( d 2 Z ( a, b )) = V ar ( d Z ( a, b )) ∼ K ( ψ , H ) a 2 H +1 for all a > 0 with C ( ψ , H ) and K ( ψ , H ) t w o p ositiv e onstan ts dep ending only on ψ and H (those results are pro v ed in Flandrin, 1992, Abry et al. , 1998). Therefore, in b oth these ases, the v ariane of w a v elet o eien ts is a p o w er la w of a , and a log-log regression pro vides an estimator of H . F rom a path ( Z (1) , . . . , Z ( N )) , the estimator will b e dedued from the log-log regression of the "natural" sample v ariane of disretized w a v elet o eien ts, i.e. , S N ( a ) = 1 [ N /a ] [ N/a ] X i =1 e 2 Z ( a, i ) . (3) A graph (log a i , lo g S N ( a i )) 1 ≤ i ≤ ℓ is dra wn from a priori family of sales and the slop e of the least square regression line pro vides the estimator b H W AV of H . In the semiparametri frame of a general lass of stationary Gaussian LRD pro esses (more general than the previous A new stochastic process to model HR series and an estimator of its fractality parameter 11 T able 1. Compar ison of th e two samples o f estimation s o f H with 100 realizati ons of fGn path (N=1000 0) wit h DF A and wav e lets method s (The p-val was ded uced from compar ison of th e mean of each sample with theore tical one at the 5% lev e l) H f Gn | B ias b H DF A | | B ias b H W AV | p-val D F A p-val W AV √ M S E D F A √ M S E W AV 0.50 0.0064 0.0071 0.0152 0.0983 0.0271 0.0433 0.60 0.0092 0.0009 0.0017 0.8289 0.0304 0.0405 0.70 0.0141 0.0015 10 − 5 0.7342 0.0347 0.0436 0.80 0.0125 0.0050 0.0002 0.1978 0.0349 0.0391 0.90 0.0179 0.0062 2 · 10 − 6 0.1030 0.0407 0.0448 semiparametri on text required for b H DF A ), it w as established b y Moulines et al. (2007) that the estimator b H W AV on v erges to H with an optimal on v ergene rate (follo wing the minimax riteria) when an optimal length of windo ws is kno wn. Th us, theoretial asymp- toti b eha viors of b H DF A and b H W AV are omparable for F GN and a semiparametrial lass of LRD Gaussian pro esses (ho w ev er more general for b H W AV ). This is not true an y more when a p olynomial trended LRD (or self-similar) pro esses is onsidered. Indeed, Abry et al. (1998) remark ed that ev ery degree M p olynomial trend is without eets on b H W AV sine ψ has its M rst v anishing momen ts. Therefore, the larger M , the more robust b H W AV is. Finally , Bardet (2002) established a  hi-squared go o dness-of-t test for a path of FBM (therefore for aggregated F GN) using w a v elet analysis. This test is based on a (p enalized) distane b et w een the p oin ts (log a i , lo g S N ( a i )) 1 ≤ i ≤ ℓ and a pseudo-generalized least square regression line (here the sales a i are seleted to b eha v e as N 1 / 3 ). In the T able 1 app ear the dieren t estimations of H omputed from the DF A and w a v elet analysis metho ds for 100 realizations of F GN paths with N = 100 00 . W e  ho ose for these sim ulations the onrete pro edure of w a v elet analysis dev elop ed b y Abry et al. (2003) (a Daub e hies w a v elet is  hosen and a Mallat's fast p yramidal algorithm is used to ompute w a v elet o eien ts). In one hand, the w a v elets metho d app ear sligh tly more eetiv e than DF A metho d onsidering the p-v alue whi h is v ery lo w for the sample of the DF A estima- tions ompared to w a v elet analysis estimations. This is essen tially due to the estimator bias whi h is more imp ortan t in the ase of DF A. In the other hand, if w e onsider the ro ot of MSE whi h is the sum of the squared bias and the v ariane, the DF A estimator seems to b e sligh tly more eetiv e. Note that for F GN pro esses (without trend), the Whittle maxim um lik eliho o d estimator of H giv es a "b etter" results (see T aqqu et al. , 1999). 3.3. Applicatio n of both the estimators to HR data Both these estimators of H an also b e applied to the HR time series of the 9 athletes. The follo wing gures Fig. 7 and Fig. 8 exhibit examples of appliations of b oth the estimation metho d to HR data. F or ea h athlete, it w as rst done to the whole time series, and then to the dieren t phases of the rae (as it w as obtained from the detetion of abrupt  hanges, see Setion 2). The estimation results of H , for the dieren t signals observ ed during the three phases of the rae, are reapitulated in the T able 2 using w a v elets metho d and in T able 3 using DF A metho d. 12 Bardet et al. 0 0.5 1 1.5 2 2.5 3 x 10 4 −1.5 −1 −0.5 0 0.5 1 x 10 4 cumsum(HR) (a) 2 2.5 3 3.5 4 1.5 2 2.5 3 3.5 log10(ni) log10(F(ni)) (b) 2 2.5 3 3.5 4 1 1.5 2 2.5 3 3.5 log10(ni) log10(F(ni)) (b) Fig. 7. T wo first steps of the DF A metho d appl ied to a HR serie s (up) and results of the DF A metho d appli ed to HR serie s f o r two different a thlet es (down) 4 5 6 7 8 9 4 6 8 10 12 14 16 18 Logarithm of IN Log. of chosen scales 4 5 6 7 8 9 4 6 8 10 12 14 16 18 Logarithm of IN Log. of chosen scales Fig. 8. The lo g-log graph of th e varian ce of wav e let coef ficients relat ing to the HR seri es obser ved duri ng the race and in the end of race (Ath2) T w o main problems result from these dieren t estimations. First, b H DF A and b H W AV are often larger than 1 . Ho w ev er, the F GN is only dened for H ∈ (0 , 1) . F or dening a pro ess allo wing H > 1 , three main assumptions of F GN ha v e to b e  hanged: (a) the assumption that the pro ess is a stationary pro ess; (b) the assumption that the pro ess is a Gaussian pro ess; () the assumption that only t w o parameters ( H and σ 2 ) are suien t to dene the pro ess. A new stochastic process to model HR series and an estimator of its fractality parameter 13 T able 2. Estimated H with wavelets methods for HR se- ries of different athletes P hases HR series Beginning Middle Rae end A th1 0.8931 1.1268 1.1064 1.2773 A th2 1.1174 0.7871 1.0916 0.8472 A th3 1.0208 1.0315 1.1797 - A th4 0.9273 - 1.0407 0.7925 A th5 1.0986 1.3110 1.0113 1.3952 A th6 1.0769 1.5020 1.1597 1.3673 A th7 1.0654 1.4237 1.1766 1.0151 A th8 0.9568 1.6600 0.9699 1.1948 A th9 0.9379 1.5791 0.9877 0.7263 In the sequel (see b elo w), a new mo del is prop osed. Both the rst assumptions are still satised and the third one is replaed b y a semiparametri assumption. The seond problem is implied b y the results of the go o dness-of-t test (for w a v elet analysis metho d). Indeed, this test is nev er aepted as w ell for the whole time series as for the par- tial times series. An explanation of su h a phenomenon an b e dedued from Figure 8 : for the w a v elet analysis, the p oin ts (log a i , lo g S N ( a i )) 1 ≤ i ≤ ℓ are learly lined for a i ≤ a m , but not exatly lined for a i ≥ a m . Th us the HR time series seems to nearly b eha v e lik e a F GN for "small" sales (or high frequenies), but not for "large" sales (or small frequenies). A pro ess follo wing this onlusion an not b e the b etter t of HR time series... Remark: this last onlusion leads also to a lear adv an tage of w a v elet based o v er DF A estimator. Indeed, the DF A algorithm measures only one exp onen t  haraterizing the en tire signal. Then, this metho d orresp onds rather to the study of "monofratal" signals su h as F GN. A t the on trary , the w a v elet metho d pro vides the graph (log a i , lo g S N ( a i )) 1 ≤ i ≤ ℓ whi h an b e v ery in teresting for analyze the m ultifratal b eha vior of data (see also Billat et al. , 2005). 4. A new model f o r modeling HR data: a locally fractional Gaussian noise 4.1. The locally fractional Gaussian noise In Bardet and Bertrand (2007), a generalization of the FBM, so-alled the ( M K ) -m ultisale FBM, w as in tro dued. The ( M 0 ) -FBM is a FBM with self-similarit y parameter H 0 . Roughly sp eaking, the ( M K ) -FBM has the same harmonizable represen tation (and therefore quite the same b eha vior as the FBM) than a FBM with self-similarit y parameter H i for frequenies | ξ | ∈ [ ω i , ω i +1 [ for all i = 0 , . . . , K ( K ∈ N ) . F or instane, a ( M 1 ) -FBM b eha v es as a FBM with self-similarit y parameter H 0 for small frequenies and as a FBM with self- similarit y parameter H 1 for high frequenies. Su h a mo del w as fruitfully used for mo deling biome hanial signals (p osition of the en ter of pressure on a fore platform during quiet p ostural stane measured at a frequeny of 100 Hz for the one min ute p erio d). Here, Fig. 9 suggests than a tted mo del for aggregated HR data should b eha v e lik e a FBM with self-similarit y parameter H for lo w frequenies and dieren tly for high frequenies 14 Bardet et al. −2 −1 0 1 2 4 6 8 10 12 14 16 18 20 Logarithm of the chosen frequencies (Hz) Logarithm of IN −4 −2 0 2 5 10 15 20 25 Logarithm of the chosen frequencies (Hz) Logarithm of IN Fig. 9. The lo g-log graph of th e varian ce of wav e let coef ficients relat ing to the HR seri es obser ved duri ng the arrival phase (Ath6 ) with a frequency band of [0.01 12](rig ht) and of [0.2 4](left ). (and not neessary lik e a FBM). Th us dene a lo ally frational Bro wnian motion X ρ = { X ρ ( t ) , t ∈ R } as the pro ess su h that: X ρ ( t ) = Z R e itξ − 1 ρ ( ξ ) c W ( dξ ) where the funtion ρ : R → [0 , ∞ ) is an ev en on tin uous funtion su h that: • ρ ( ξ ) = 1 σ | ξ | H +1 / 2 for | ξ | ∈ [ ω 0 , ω 1 ] with H ∈ R , σ > 0 and 0 < ω 0 < ω 1 • Z R  1 ∧ | ξ | 2  1 ρ 2 ( ξ ) dξ < ∞ . and W ( dξ ) is a Bro wnian measure and c W ( dξ ) its F ourier transform in the distribution meaning. Cramér and Leadb etter (1967) pro v ed the existene of su h Gaussian pro ess with stationary inremen ts. The main adv an tages of su h pro ess ompared to usual FBM are the follo wing: (a) X ρ "b eha v es" lik e a FBM only for lo al band of frequenies; (b) In this band, the parameter H is not restrited to b e in (0 , 1 ) : it is in R . F rom this denition, one dedues a p ossible mo del for HR data: Y ρ ( t ) = X ρ ( t + 1) − X ρ ( t ) = 2 · R e  Z R e itξ sin( ξ / 2) ρ ( ξ ) c W ( dξ )  for t ∈ R . Note that Y ρ = { Y ρ ( t ) , t ∈ R } is a stationary Gaussian pro ess and the funtion 2 sin( ξ / 2) ρ − 1 ( ξ ) is so-alled the sp etral densit y of Y ρ . A new stochastic process to model HR series and an estimator of its fractality parameter 15 Let ∆ N → 0 and N ∆ N → ∞ when N → ∞ . The w a v elet based estimator an pro vide a on v ergen t estimation of H when a path ( Y ρ (∆ N ) , Y ρ (2∆ N ) , . . . , Y ρ ( N ∆ N )) and therefore a path ( X ρ (∆ N ) , . . . , X ρ ( N ∆ N )) is observ ed. Indeed, onsider a "mother" w a v elet ψ su h that ψ : R 7→ R is a C ∞ funtion satisfying : • for all s ≥ 0 , Z R | t s ψ ( t ) | dt < ∞ ; • its F ourier transform b ψ ( ξ ) is an ev en funtion ompatly supp orted on [ − β , − α ] ∪ [ α, β ] with 0 < α < β . Then, using results of Bardet and Bertrand (2007), for all a > 0 su h that [ α a , β a ] ⊂ [ ω 0 , ω 1 ] , i.e. a ∈ [ β ω 1 , α ω 0 ] , ( d X ρ ( a, b )) b ∈ R is a stationary Gaussian pro ess and E  d 2 X ρ ( a, . )  = K ( ψ , H, σ ) · a 2 H +1 , with K ( ψ , H , σ ) > 0 only dep ending on ψ , H and σ . Ho w ev er this prop ert y is  he k ed if and only if the funtion ψ is  hosen su h that: β α < ω 1 ω 0 . Moreo v er, for a ∈ [ β ω 1 , α ω 0 ] , the sample v ariane S N ( a ) dened in (3) and omputed from a path ( X ρ (∆ N ) , . . . , X ρ ( N ∆ N )) on v erges to E  d 2 X ρ ( a, . )  and satises a en tral limit theo- rem with on v ergene rate √ N ∆ N . Th us, with xed sales ( a 1 , . . . , a ℓ ) ∈ [ β ω 1 , α ω 0 ] ℓ , a log- log-regression of ( a i , S N ( a i )) 1 ≤ i ≤ ℓ pro vides an estimation of H (and a en tral limit theorem with on v ergene rate N ∆ N satised b y b H W AV an also b e established). As previously , w e onsider also  hi-squared go o dness-of-t test based on the w a v elet analysis and dened as a w eigh ted distane b et w een p oin ts (log( a i ) , log ( S N ( a i ))) 1 ≤ i ≤ ℓ and a pseudo-generalized regression line. Remark: The main problem with these estimator and test is the lo alization of the suitable frequeny band [ ω 0 , ω 1 ] ( ω 0 and ω 1 are assumed to b e unkno wn parameters). A solution onsists in seleting a v ery large band of sales and determining then graphially the "most" linear part of the set of p oin ts (log( a i ) , log ( S N ( a i ))) 1 ≤ i ≤ ℓ . Another p ossible w a y ma y b e to ompute an adaptiv e estimator of this band using a quadrati riterion (follo wing a similar pro edure than in Bardet and Bertrand, 2007). Here, lik e 9 dieren t paths of HR data are observ ed, a ommon frequeny band [ ω 0 , ω 1 ] an b e graphially obtained and used for whole HR data (see ab o v e). 4.2. Applicatio n to HR data First, one onsiders that a HR time series ( Y (1) , . . . , Y ( n )) an b e written ( Y ρ (∆ n ) , Y ρ (2∆ n ) , . . . , Y ρ ( n ∆ n )) , Y ρ = { Y ρ ( t ) , t ∈ R } a pro ess dened as previously . Seondly , the w a v elet 16 Bardet et al. T able 3. Estimated b H , with DF A and wavelets me thods, for HR ser ies of dif ferent at hlete s ( ∗ ) Th e serie s for which the test i s rejected . Co mpari son of the two samp les ( b H D F A ) 1 ,..., 9 and ( b H W AV ) 1 ,..., 9 for whole and par tial series (p-value) HR series Rae b eginning During the rae End of rae b H D F A b H W AV b H D F A b H W AV b H D F A b H W AV b H D F A b H W AV A th1 0.928 1 . 288 ∗ 1.032 1 . 192 1.060 1 . 214 0.429 1 . 400 A th2 1.095 1 . 268 ∗ 0.905 0 . 973 1.126 1 . 108 1.240 1 . 452 ∗ A th3 1.163 1 . 048 0.553 0 . 898 1.130 1 . 172 - - A th4 1.193 0 . 916 ∗ - - 1.098 1 . 249 ∗ 1.172 1 . 260 A th5 1.239 1 . 110 1.267 1 . 117 ∗ 1.133 1 . 205 1.273 1 . 348 ∗ A th6 1.247 1 . 084 ∗ 1.237 1 . 106 1.091 1 . 172 1.436 1 . 338 A th7 1.155 1 . 095 0.850 1 . 295 1.182 1 . 186 ∗ 1.129 1 . 209 A th8 1.258 1 . 011 1.304 1 . 128 ∗ 0.995 1 . 134 1.122 1 . 247 A th9 1.243 1 . 429 ∗ 0.820 1 . 019 1.127 1 . 535 ∗ 1.250 1 . 238 ∗ p-value 0 . 6414 0 . 3723 0 . 0225 0 . 1260 F-stat 0 . 23 0 . 85 6 . 38 2 . 65 analysis is applied to the 9 (whole or partial) HR time series (the  hosen "mother" w a v elet is a kind of Lemarié-Mey er w a v elet su h that β = 2 α ). Using rst a v ery large band of sales for all HR time series (for example [0.01, 12℄ in Fig. 10 ), one estimation of frequeny band is dedued: [ ω 0 , ω 1 ] = [0 . 2 , 4 ] is the  hosen frequeny band for the whole and partial signals. −6 −4 −2 0 2 5 10 15 20 25 Logarithm of the chosen frequencies (Hz) Logarithm of IN −3 −2 −1 0 1 6 8 10 12 14 16 18 20 Logarithm of the chosen frequencies (Hz) Logarithm of IN Fig. 10. The log-l og g raph of the variance of wavelet coefficient s relating to the HR ser ies ob ser ved in the middle of the ex ercise (Ath5) The estimation results of H, for the dieren t signals observ ed during the three phases of the rae, are reapitulated in the T able 3 . Both DF A and w a v elet analysis metho ds pro vide estimations of Hurst exp onen t whi h re- et the p ossible mo deling of HR data with long range dep endene time series. A new stochastic process to model HR series and an estimator of its fractality parameter 17 W e also note that with a p-v alue of 0 . 64 , b oth the samples ( b H DF A ) 1 ,..., 9 and ( b H W AV ) 1 ,..., 9 obtained from all HR time series are signian tly lose. The same omparison an also b e done when the three  harateristi stages of the rae (b e- ginning, middle and end of the rae) are distinguished. The result is dieren t. Indeed, the orresp onding p-v alues b et w een ( b H DF A ) 1 ,..., 9 and ( b H W AV ) 1 ,..., 9 are signiativ ely dieren t in the middle part of the rae (and relativ ely dieren t in the stage of rae end). In spite of v alues relating to the estimator of H for all the athletes in the dieren t phases whi h are relativ ely large, the DF A has sometimes tendeny to under estimating this param- eter lik e in the rae b eginning (A th3) and the end of rae (A th1). Indeed, these v alue are learly due to a ertain trend supp orts b y the fat that data p oin ts in log-log plot (Fig. 11 ) ha v e not a straigh t line form, and w e ha v e pro v ed in Bardet and Kammoun (2008 ) that the DF A metho d is not robust in the ase of trended long range dep enden t pro ess. Ho w ev er in b oth the ases, the w a v elets metho d is more eetiv e sine it remo v es suien tly this kind of trend. 500 1000 1500 2000 2500 130 140 150 160 170 180 190 HR (BPM) 1 1.5 2 2.5 3 0 0.5 1 1.5 2 log10(ni) log10(F(ni)) 500 1000 1500 2000 100 120 140 160 HR (BPM) 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 log10(ni) log10(F(ni)) Fig. 11. The results of the DF A method appli ed to records f or race beg inni ng (Ath 3) ( left) and for end of race (Ath1) (right) F or HR data and when the go o dness-of-t test is aepted, the w a v elet metho d sho ws a fratal parameter H lose to 1 . A ording to the dieren t studies (using DF A metho d) ab out ph ysiologi time series for distinguishing health y from pathologi data sets (see Goldb erger et al. (2002 ), P eng et al. (1995 ), P eng et al. (1993 )), an exp onen t H ≃ 1 indiate a health y ardia HR time series. Indeed, for the study onerning a 24 hours reorded in terb eat time series during the exerise for health y adults and heart failure adults, the follo wing results are obtained: for health y sub jets, H = 1 . 01 ± 0 . 16 , for the group of heart failure sub jets H = 1 . 24 ± 0 . 22 . During the dieren t stages of the marathon rae, a small inrease of the fratal param- eter H is observ ed esp eially at the end of raes. This b eha vior and this ev olution ma y b e asso iated with fatigue app earing during the last phase of the marathon. This ev olution an not b e observ ed with DF A metho d. Indeed, in one hand, when w e observ e the three 9 -samples of w a v elet estimators (related to the 3 phases of the rae), the p-v alue (see Fig. 18 Bardet et al. 12 ) indiates a signian tly dierene due preisely to this ev olution of the fratal param- eter. On the other hand, a large p-v alue (0.85) is obtained for the same test using DF A estimation. b H DF A b H W AV p-value 0.8570 0.0158 F-stat 0.16 5.27 Race beginning During the race End of race 1 1.1 1.2 1.3 1.4 1.5 H WAV Boxplot of the three samples of estimations by applying wavelet analysis Fig. 12. Compa rison of the three samples constituti ng by estimati ons in the beg innin g of race, during the race and then in the race end by the DF A a nd wav elet methods The represen tation giv en b y Fig. 12 , highligh t a dierene in the b eha viors of HR series in the b eginning of the rae and in the end of rae. Indeed, the disp ersions in the rst and last sample are more imp ortan t than in the middle of rae and it seems that ea h athlete starts and nished the rae at his o wn rh ythm but in the middle athletes seems to ha v e the same rate. 5. Conclusion As indiated in the b eginning of the last setion, our main goal is to see whether the heart rate time series during the rae ha v e sp ei prop erties that of saling la w b eha vior. The w a v elet analysis and the DF A metho ds are applied to 9 HR time series during the whole and also the dieren t three phases of the rae (b eginning, middle and end of rae) obtained b y an automati pro edure. Ev en if their results are not exatly the same, b oth metho ds pro vide Hurst exp onen ts whi h reet the p ossible mo deling of HR data b y a LRD time series. Ho w ev er, in Bardet and Kammoun (2008 ), ev en if the DF A estimator of Hurst pa- rameter is pro v ed to b e on v ergen t with a reasonable on v ergene rate for LRD stationary Gaussian pro esses, it is not at all a robust metho d in ase of trend. The w a v elet based metho d pro vides a more preise and robust estimator of the Hurst parameter. Th us, the results obtained from this w a v elet estimator seem to b e more v alid. Moreo v er, a  hi-squared go o dness-of-t test an also b e dedued from this metho d. It seems to sho w that a lassial LRD stationary Gaussian pro ess is not exatly a suitable mo del for HR data. Graphs obtained with w a v elet analysis also sho w that a lo ally fra- tional Gaussian noise, a semiparametri pro ess dened in Setion 3 ould b e more relev an t to mo del these data. A  hi-squared test onrms the go o dness-of-t of su h a mo del. Th us, using the w a v elet estimation of a fratal parameter in a sp ei frequeny band, one obtains a onlusion relativ ely lose to those obtained b y other studies (onlusion whi h an not b e deteted with DF A metho d): these fratal parameters inrease through the rae phases, what ma y b e explained with fatigue app earing during the last phase of the marathon. Th us A new stochastic process to model HR series and an estimator of its fractality parameter 19 this fratal parameter ma y b e a relev an t fator to detet a  hange during a long-distane rae. Finally , for the 9 athletes and as the test is v alidated with signiane lev el around 0.65, w e an estimate b H beginning at 1 . 1 , the b H middle at 1 . 2 and b H end at 1 . 3 with a larger ondene in terv al at the b eginning and the end of the rae. This b eha vior ould bring a new w a y of understanding what is happ ening during a rae. References Abry , P ., Flandrin, P ., T aqqu, M. and V eit h, D. (2003) Self-similarit y and long-range de- p endene through the w a v elets lens. The ory and appli ations of L ong-R ange Dep enden e , 527-556, Birkhäuser. Abry , P ., V eit h, D. and Flandrin, P . (1998) Long-range dep enden t: revisiting aggregation with w a v elets. J. Time Ser. A nal. , 253-266. Absil, P .A., Sepul hre, R., Bilge, A. and Gérard, P . (1999) Nonlinear analysis of ardia rh ythm utuations using DF A metho d. J. Physi a A: Statisti al me hanis and its ap- pli ations , 235-244. Bagger, M., P etersen, P .H. and P edersen, P .K. (2003) Biologial V ariation in v ariables asso iated with exerise training. Int. J. Sp orts Me d , 433-440. Bai, J. (1998) Least squares estimation of a shift in linear pro esses. J. of Time Series A nal. , 5, 453-472. Bai, J. and P erron, P . (1998) Estimating and testing linear mo dels with m ultiple strutural  hanges. E onometri a , 66, 47-78. Bardet, J.M. (2002) Statistial study of the w a v elet analysis of frational Bro wnian motion. IEEE T r ans. Inform. The ory. , 48, 991-999. Bardet, J.M., Bertrand, P . (2007) Iden tiation of the m ultisale frational Bro wnian mo- tion with biome hanial appliations. Journal of Time Series A nalysis , 28, 1-52. Bardet, J.M. and Kammoun, I. (2008) Asymptoti Prop erties of the Detrended Flutua- tion Analysis of Long Range Dep endene Pro esses. T o app ear in IEEE T rans. Inform. Theory . Bardet, J.M., Lang, G., Moulines, E. and Soulier, P . (2000) W a v elet estimator of long-range dep enden t pro esses. Statisti al Infer en e for Sto hasti Pr o  esses , 3, 85-99. Bardet, J.M., Lang, G., Opp enheim, G., Philipp e, A. and T aqqu, M.S. (2003) Generators of long-range dep enden t pro esses : A surv ey . L ong-r ange Dep enden e : The ory and Appli ations, Birkhauser . Billat, V., W esfreid, E., Cottin, F., Kapfer, C. and K oralztein, J.P . (2003) F ratal analysis of sp eed and ph ysiologial osillations in long- and middle-distane running: Eet of train- ing. International Journal of Computer Sien e in Sp ort , 16-30. h ttp://www.billat.net/ . 20 Bardet et al. Flandrin, P . (1992) W a v elet analysis and syn thesis of frational Bro wnian motion. IEEE T r ans. on Inform. The ory , 38, 910-917. Goldb erger, A.L. (2001) Heartb eats, hormones, and health : is v ariabilit y the spie of life? A m J R espir Crit Car e Me d , 163, 1289-1290. Goldb erger, A.L., Amaral, L.A.N., Hausdor, J.M., Iv ano v, P .C., P eng, C.K. and Stanley , H.E. (2002) F ratal Dynamis in Ph ysiology: Alterations with Disease and Aging. PNAS , 99, 2466-2472. La vielle, M. (1999) Detetion of m ultiple  hanges in a sequene of random v ariables. Sto h. Pr o  ess Appl , 79-102. h ttp://www.math.u-psud.fr/ ∼ la vielle/ . La vielle, M. and Moulines, E. (2000) Least-squares estimation of an unkno wn n um b er of shifts in a time series. J. of Time Series A nal. , 21, 33-59. Moulines, E., Roue, F. and T aqqu, M.S. (2007), On the sp etral densit y of the w a v elet o eien ts of long memory time series with appliation to the log-regression estimation of the memory parameter. J. Time Ser. A nal. , 28, 155-187. Moulines, E. and Soulier, P . (2003). Semiparametri sp etral estimation for frational pro esses. In P . Doukhan, G. Opp enheim and M.S. T aqqu editors. The ory and appli ations of long-r ange dep enden e , 251301, Birkhäuser, Boston. P eng, C.K., Buldyrev, S.V., Ha vlin, S., Simons, M., Stanley , H.E. and Goldb erger, A.L. (1994) Mosai organization of DNA n uleotides. Physi al R eview E , 49, 1685-1689. P eng, C.K., Ha vlin, S., Stanley , H.E. and Goldb erger, A.L. (1995) Quan tiation of saling exp onen ts and rosso v er phenomena in nonstationary heartb eat time series. Chaos , 5, 82. P eng, C.K., Mietus, J., Hausdor, J., Ha vlin, S., Stanley H.E. and Goldb erger A.L. Long- Range An tiorrelations and Non-Gaussian Beha vior of the Heartb eat. Phys. R ev. L ett. , 70, 1343-1346. Pikkujämsä, S.M., Mäkik allio, T.H., Airaksinen, K.E. and Huikuri, H.V. (2001) Deter- minan ts and in terindividual v ariation of R-R in terv al dynamis in health y middle-aged sub jets. A m. J. Physiol. He art Cir . Physiol. , 280, 1400-1406. Samoro dnitsky , G. and T aqqu, M.S. (1994) Stable non-Gaussian Random Pro esses. Chap- man and Hal l . Sto ev, S., T aqqu, M.S., P ark, C. and Marron, J.S. (2005) On the w a v elet sp etrum di- agnosti for Hurst parameter estimation in the analysis of In ternet tra. Elsevier , 49, 423-445. T aqqu, M.S. and T ev ero vsky , V. (1997) Robustness of whittle t yp e estimators for time series with long range dep endene. J. Sto hasti Mo del , 13, 723-757. T aqqu, M.S., T ev ero vsky , V. and Willinger, W. (1995) Estimators for long-range dep en- dene: an empirial study . F r atals , 3, 785-798. W esfreid, E., Billat, V. and Mey er, Y. (2005) Multifratal analysis of heartb eat time series in h uman raes. Appl. Comput. Harmon. A nal. , 18, 329-335.