# Full text of "Entropic Analysis of non-Stationary Sequences"

## See other formats

Entropic Analysis of non-Stationary Sequences (N O o (N 5h . < (N 43 O CO I C O o > m (N (N O <N O -a c o o X S3 Michele Virgilio 1 , and Paolo Grigolini 1 ' 2,3 . 1 Dipartimento di Fisica dell'Universitd di Pisa and INFM, Via Buonarrotti, 56127 Pisa, Italy 2 Center for Nonlinear Science, University of North Texas, P.O. Box 311427, Denton, Texas 76203-1427 3 ' Istituto di Biofisica CNR, Area della Ricerca di Pisa, Via Alfieri 1, San Cataldo 56010 Ghezzano-Pisa, Italy (February 1, 2008) The aim of this paper is to shed light on the analysis of non-stationary time series by means of the method of diffusion entropy. For this purpose, we first study the case when in- finitely many time series, as different realizations of the same dynamic process, are available, so as to adopt the Gibbs en- semble perspective. We solve the problem of establishing un- der which conditions scaling emerges from within this perspec- tive. Then, we study the more challenging problem of creating a diffusion process from only one single (non-stationary) time series. The conversion of this single sequence into many dif- fusional trajectories is equivalent to creating a non-Gibbsian ensemble. However, adopting a probabilistic approach to eval- uate the contribution of any system of this non-Gibbsian en- semble, and using for it the theoretical Gibbsian prescription of the earlier case, we find a recipe that fits accurately the numerical results. With the help of this recipe we show that non-stationary time series produce either anomalous scaling with ordinary statistics or ordinary scaling with anomalous statistics. From this recipe we also derive an attractive way to explain the entropy time evolution, as resulting from two distinct uncertainty sources, the missing information on tra- jectory initial condition and the randomness induced lack of control of trajectory time evolution. 89.75.D ,05.70.L,02.30.L ,95.75.W I. INTRODUCTION One of the paradigms of the Science of Complexity is that complex systems, of biological, physical, economical and sociological interest, are characterized by long-range correlations. Thus, a time series obtained by recording data of this kind is expected to reflect these long-range correlations. How to detect these correlations in prac- tice? One the pioneer papers in this field [[lj rests on the study of diffusion process generated by the time series, whose numbers are imagined to be diffusion generating fluctuations. Thus, the detection of these correlations can be made in an indirect way by observing the de- viation of the resulting diffusion process from ordinary Brownian diffusion. To set this important issue on a solid ground, let us refer ourselves to the scaling property of the probabil- ity distribution function (pdf) p(x,t). This means that we are assuming that the time series is transformed into many independent, one-dimensional, diffusion trajecto- ries x(t), and that p(x, t) is the probability density of finding the trajectory in position i at a given time t. Scaling means that p(x,t) = ±F(^). (1) The theoretical foundation of this property is given by the Central Limit Theorem (CLT) j| and by the Gen- eralized Central Limit Theorem (GCLT) ||. In short, these two theorems refer to the case when the variable x is the sum of many uncorrelated variables This is the case when the time series is a sequence of independent fluctuations occurring at the integer times i, and x is the sum of the fluctuations ranging from i = j to i = j + I, with I playing the role of an integer time that in the lim- iting case I 3> 1 is identified with the continuous variable t. Let us now focus our attention on the fluctuations £j, and for simplicity, let us make the assumption that their probability distribution has a vanishing mean value. Let us consider, first, the case when the second moment of distribution is finite. In this case, the CLT applies. The validity of the CLT is broken when the second moment is divergent. However, in this case, we meet the wonderful opportunity of replacing the CLT with the GCLT. In the former case S = ^, and F(y) is a Gaussian function. In the latter case, 5 > h and F(y) is a Levy function, whose asymptotic form is lim\ y \^ 00 F{y) = l/\y\ 1+1 / s . A third case exists, the popular case of Fractional Brownian Mo- tion (FBM), introduced by Mandelbrot M. In this case the fluctuations are correlated, S range from to 1, and F{y) is a Gaussian function again. The dynamic foun- dation of this third case is not yet as solid as that of the earlier two cases, and conjectures have been made as to a possible connection between FBM and non-stationary properties [pi. As to the important problem of detecting scaling from the statistical analysis of real time series, the most popu- lar approach is that developed by the authors of Ref . M . These authors developed a technique of scaling detec- tion based on the analysis of the variance of the diffusion process generated by the time series under study. This technique of analysis was termed Detrended Fluctuation Analysis (DFA) , and has been widely applied since then. In the last few years another method was proposed by the authors of Ref. || . These authors introduced a tech- nique of scaling detection, based on the evaluation of the 1 entropy of the experimental pdf p(x,t), called Diffusion Entropy Analysis (DEA) . We note that in the continuous time limit the Shannon entropy of the diffusion process reads /+oo dxp(x,t)log(p(x,t)). (2) -oo By plugging Eq.(^) into Eq.(|^) and using a change of integration variable, we get S(t) = A + Slogt, (3) where A = - [ dyF(y)log(F(y)). (4) J — oo It is evident that this method can be used to determine 5 without the limitations affecting the DEA. In fact, the DEA works, in principle, even when F(y) does not have a finite second moment, while the DEA, being based on the variance calculation, cannot be applied in the case of Levy statistics. It is important to stress that both DFA and DEA rest on the creation of a Gibbs ensemble from a single se- quence. Thus, to ensure that the systems of this ensem- ble are identical copies of the same system, it is neces- sary to assume a stationary condition. Unfortunately, many complex processes seem to be incompatible with the stationary condition. An interesting example of non- stationary condition is given, for instance, by the teen birth phenomenon studied in Ref. ||. The time series under study in the paper of Ref. || concern the number of babies born per day to teens in Texas. Since fertil- ity has an obvious annual periodicity, and the rate itself of sexual intercourses is holidays dependent, and conse- quently, is affected, too, by annual periodicities, these kinds of time series are not stationary. The fluctuations of interest for statistical analysis take place around a time dependent mean value dictated by both social and bio- logical deterministic processes. If the annual periodicities are not properly taken into account, the social process under study might seem to be characterized by memory effects of surprisingly high intensity. Another case of ex- ternal modulation of fluctuation intensities is given by the solar flares , where the rate of flares is modulated by the 11-years solar cycle. The main purpose of this paper is to establish which is the consequence of addressing the study of a non- stationary time series with a method originally designed for stationary time series. The paper of Ref. || addresses this important issue using DFA. This paper rests on the use of DEA. II. GIBBS ENSEMBLE We notice that the statistical analysis of time series which reflect rules changing with time meets with a big conceptual difficulty. This has to do with the fact that a statistical treatment in physics is usually done with the picture of Gibbs in mind. We imagine that the system under study can ideally be repeated a virtually infinite number of times. The copies are exactly identical to the original, and one copy differs from another only for dif- ferent initial condition. In the case where the system under study undergoes the influence of an environment, expressed under the form of a fluctuation, there is not even the need to use different initial conditions. Moving from one copy to another, the time evolution of the envi- ronmental fluctuation changes, in spite of the constraint that the statistical property of the environmental fluctu- ation remains unchanged. The statistical treatment con- sists of averaging on this set of virtually infinite copies. If the system under study is, partially or totally, random, an even slight difference on initial conditions will create deep differences at later times. The same condition ap- plies to the case of a stochastic system, even if the initial conditions can be chosen to be exactly identical in this case. The observation of the resulting trajectory spread- ing makes it possible to determine in practice the corre- sponding pdf. In the specific case of a time series, we do not have infinitely many identical copies of the same system. However, if the time series is stationary, we can create these copies with the method of the moving win- dow. We move a window of size / along the sequence and all the values spanned by the window located at a given position are regarded as values, random, partially random, or stochastic, generating a diffusion trajectory of time length equal to I. This cannot be done in the non-stationary case, since moving the window from an old position, with respect to sequence reference frame, to a new one, is equivalent to considering a copy that it is not equivalent to the original copy. We are aware of this difficulty. Nevertheless, in Section 3 we shall study non- stationary time series pretending they are a stationary, and we shall try to see if the resulting effects can give us information on the non-stationary nature of the process. To deepen our understanding of those effects, it is conve- nient to study first artificial non-stationary time series. We create a virtually infinite number of copies and we adopt the Gibbs point of view. This is the subject of this section. The statistical analysis of time series of the earlier work (see, for instance, Ref. |9| as well as Ref. Q) is based on turning a time series into a diffusion process. The time series {^} is a sequence of real numbers that are inter- preted as fluctuations. The sum of / of these fluctuations is called x(l) and can be considered as the position of a walker at time I. If the sequence is large enough we can 2 consider values of I so large as to interpret / as a contin- uous time variable that we denote by a different symbol, t, to emphasize the adoption of the continuous time ap- proximation. This means that we study the equation of motion dx ... (5) To mimic the effect of using a non-stationary time se- ries, we write the fluctuation £(t) as follows (G) where g(t) and h(t) are regular functions of time and £°(t) a memoryless stochastic process defined by <e(ti)e(t2)>=s(\h-t 2 \). (7) The function g(t) serves the purpose of making non sta- tionary the noisy portion of the time series under study and the term h(t) is inspired to the teen birth phe- nomenon ||. As mentioned in Section 1, fertility has annual periodicity, and its influence has been mimicked by a term h(t), with a sinuisoidal form. Let us make the change of variables y(t) _ x{t)-j t h{t')dt> g{t) (8) and let us consider the stochastic process described by (9) §=s(t)£ (0) (t) + fc(t). As a consequence of this change of variable, Eq.(|] ) with defined by Eq.ffl), becomes (10) The equation for the probability distribution p(x, t), cor- responding to this stochastic equation, is well known to be dp(y, t) at g(t) dy d 2 y (11) Let us look for a solution of Eq. flll|) with the form: 2 p(y,t) = [2^ (t)} V2 (12) By plugging Eq. (|12|) into Eq. (|ll|) we obtain, after some easy algebra, d<rj(t) dt 2 9 ^-a*(t)-2D = 0. git) v (13) This means that Eq. (12) is a solution of Eq.([LT]) pro- vided that the variance (t^ (t) is a solution of Eq.(|13[) . It is straigthforward to prove that 2D ~gW) g 2 (t')dt'- (14) We have selected a solution fitting the condition a y (0) = 0, if g(0) 0, which corresponds to the way we carry out the statistical analysis of time series ||[)| . This is so be- cause all diffusion trajectories derived from the adoption of mobile windows are assumed to start from x = at t = 0. Due to Eq.(||) this means y = at t = 0. It is straightforward to prove that the probability distribution p(x,t) is then given by p{x,t) 1 m [2^(i)]i/2 with and al{t)=g\tyjt) x (t) h{t')dt' (15) (16) (17) We are now ready to answer the question as to whether the scaling condition of Eq. ([!]) is compatible with an out of equilibrium condition. It is evident that the conditions to fulfill arc: and a z x {t) = 2D / g\t')dt' = at Jo (x-x (t)) 2 _ {x-f h{t')dt'f 2.1 (18) M( ¥ ) (19) 2a%{t) 2at 2S By differentiating Eq.(p"8"|) with respect to time, we get g (t) = (^h*-^. This leads us to the form g(t) = UP, with 1 < (3 < oo. (20) (21) (22) By expanding the square appearing in Eq. (|19[), we find for the function hit) the following form (23) h(t) = ct Thus we conclude that the forms of Eq. (|2lj) and Eq. (|23|) are compatible with the scaling condition of Eq. ([!]) with the scaling parameter S given by 5 = (3 + 1 (24) We note that the linear case, with (3=1, yields 5 = 3/2, the value obtained as inertial effects in Refs. pfl|-^2||. 3 III. SINGLE SEQUENCES We are now in a position to address the intriguing problem of non-stationary time series when only one time series is available. We do not limit ourselves to applying the DE. We also interpret theoretically the numerical re- sult stemming from the DE. It is the proper time here to concisely review how a single sequence is converted into many trajectories. Let us imagine that we have the se- quence From a theoretical point of view, the ideal condition to use would be that the sequence is infinite. Unfortunately, this condition is not fulfilled by real se- quences, and, for obvious computational limitations, is not fulfilled by artificial sequences either. Thus, we ex- plicitly take this condition into account, and we denote the sequence finite size by l seq . We assume that l seq , al- though finite, is large enough as to ensure a fair numerical evaluation of the statistical properties under study. . We denote the maximum temporal length of the trajecto- ries equal to l m ax- Thus, to make predictions concerning diffusion occurring at time t = l max we can use up to heq — Unax + 1 distinct trajectories. These trajectories are not independent of one another. In fact, the discrete time derivative of the s-th trajectory is defined by xW(Hl)-xW(t)=6+f (25) This means the s-th trajectory and the (s + l)-th trajec- tory, both of lenght I, will have I — I numbers in common, and that only if \s — s'\ > l — l the two trajectories are to- tally independent of one another. It is evident that if we adopted the criterion of non-overlapping windows, as the authors of Ref. [0 do, we would get the benefit of using trajectories totally independent of one another without meeting the risk of introducing correlations that do not exist. On the other hand, the number of trajectories to use would be much smaller, and the resulting statis- tics would be poorer. Thus, we decided to check if the method of overlapping windows yields spurious effects or not. The adoption of artificial sequences makes it possi- ble for us to settle this important issue. In fact, we did the calculation using the method of overlapping windows and we got a given result. We did again the same calcu- lation with the same number of trajectories, these being totally independent the ones of the others. This means that we adopt the same stochastic prescription, but we run it several times (one run for each trajectory) so as to create really independent trajectories. The new result turned out to coincide with the earlier, thereby leading us to conclude that the method of overlapping windows can be safely applied. This is in fact the method that we adopt throughout. The key formula for the theoretical treatment of the numerical results of this Section is given by p(M) = Eh-~"+V"(M) , (26) The meaning of this formula is the following. We make the conjecture that the probabilistic contribution of the s-th trajectory to the pdf resulting from this method of numerical observation is the same as that stemming from an arbitrarily large number of trajectories corresponding to the same physical condition. We do not have a rigorous proof for this theoretical expression, but the extremely good agreement between the theoretical prediction itself and the numerical observation. On the other hand, once this theoretical prediction is accepted, the physical inter- pretation of the numerical results become easier, and as we shall see, we can use this formula to account for two distinct sources of entropy increase, the trajectory ran- domness and the uncertainty on the initial conditions. A. Linear increase of noise intensity Let us consider first the case where the noise intensity increases linearly with time, H(t)=^\t){a + bt). (27) The sequence corresponding to the s-th window, becomes ^ s \t) = ^°\t + s)(a + b S + bt). (28) The s-th diffusion trajectory obeys the equation of mo- tion iW =£<f>)(t + s)(a + b8 + bt), (29) which is equivalent to x^(t) = £W(t)(a + bs + bt), (30) due to the assumption that the noise ^°- > is memoryless. This means that we can refer ourselves to Eq.(^) with g(t) =a + bs + bt, h(t) = 0. (31) We plug g(t) of Eq.@ into Eq.((l|). This can be easily integrated and yields 2D t) = [(a + bs + bt) 3 -(a + bs) 3 }. (32) Eq.(p6|) becomes 1 e ~^ih) (2tt) i -' 2 N a x (s,t) where N = l seq - l max + I. (34) For N — > oo, the sum depends on values of s so large as to make a x (s,t) very close to 2D(a + bs) 2 t. This makes 4 it is easy to assess that p(x, t) becomes the sum of distri- butions rescaling with 5 = 1/2. Thus also p(x, t) rescales with 6 = 1/2. This is confirmed by the numerical result illustrated in Fig. 1. It is remarkable that the shape of the resulting distri- bution (see Fig. 2) looks quite different from that of a Gaussian function, as made even more evident from the adoption of the logarithmic scale for the ordinate axis (see Fig. 3). The resulting distribution is very well ap- proximated by the fitting function Pfit M = M!) e -fe(*)M. (35) It is remarkable that an exponential form of this kind has been found by Metzler and Klafter |ll| . Apparently, the case studied by Metzler and Klafter is quite different from the one illustrated in Figs. 2 and 3. The authors of Ref. |H| find that the exponential structure is associated to the subdiffusional scaling 5 < 1/2. Here the scaling, as shown in Fig. 3, is the normal scaling 5=1/2. However, the similarity might not be accidental. The case dis- cussed by Metzler and Klafter refers to a random walk with long rests and a waiting time, between a jump and the next, given by a distribution whose asymptotic limit is the same as that of V(t) = ((i - 1) JT/J- (T + T y ' (36) Sequences of time with this distribution has been stud- ied recently from the point of view of entropy production |jl4| . It has been shown, in agreement with the ear- lier work of Gaspard and Wang jl5), that the entropy production per unit of time is proportional to tf l ~ 2 . This means that the rate of entropy production is not constant and that for t — > oo it tends to vanish. This is a form of non-stationarity that is apparently different from the one studied in this paper. However, in both cases a scaling distribution emerges, and in both cases a striking devia- tion from the Gaussian form is exhibited. We have also to mention that if we adopt the e- entropy prescription p6[ to analyze our non-stationary time series, this form of entropy, being proportional to the square of the noise intensity, tends to increase with time. Thus, from the point of view of entropy production both sequences are non-stationary, and we wonder if the similarity between the exponential distribution of Fig. 3 and the exponen- tial distributions of Metzler and Klafter is more than accidental. Adapting the approach of Section 3 A to this case, we get . sin2uj(t + s) — sin2ujs . .„„. a 2 x (s,t)=a 2 D(t+ ). (38) It is evident that in the asymptotic time limit the linear term becomes more important than the harmonic term, which, consequently, can be neglected, thereby forcing the contributions to the distribution with different s, see Eq. (Eq), to become identical to p(x,t) 1 1/2 C ' (39) with a 2 (t) = a 2 Dt, (40) which yields for the DE the following analytical expres- sion S(t) = ±logt+^log(2ireDa 2 ). (41) Fig. 4 compares this analytical prediction to the numer- ical result and shows that, as expected, the theoretical prediction departs from the numerical result in the time region t < l/u>. C. Time series with a sinuisodal deterministic part Let us study the case where g(t) — a and h(t) = bcosuit. In this case, according to the theory of Section 2, for the contribution corresponding to the s-th window, we must specify the form of xo(t) as well as that of a x (t). With the same method of Section 3 A we get: and al(s,t) = 2Da 2 t. x (s,t) = / h(t' + s)d£. (42) (43) Using h(t) — bcosujt, with a straigtforward integration, we write x (s, t) = I(t)sin(u>s + (44) with B. Time series with noise intensity sinuisoidally dependent on time Let us consider the case where g(t) = acos(u)t), h(t) = 0. (37) = _[2(1 - cosuot)] 1 / 2 . (45) The key formula of Eq. (BO) in this case takes the form p(x,t) 1 N(2ira 2 t) 1 / 2 J 1 e ds (46) 5 > 2 In the lone-time limit t >> ■> 2 n , the contribution of xo(s,t) becomes negligible and we recover Eq. (J3Sj) with cr 2 (i) = 2Dt. However, in this case the regime of transition to scaling is much more interesting than that of Section 3 B. At early times we see a process of entropy increase that might be easily mistaken for a scaling regime with anomalous scaling d > 0.5. This is caused by the deterministic process h(t). The adoption of the method of moving windows implies the adoption of many trajectories, all being the solution of the harmonic equation jj^x + uj 2 x = 0, with different initial condi- tion. This means that in this case the entropy increase does not depend on the random nature of the time se- ries, but rather on the uncertainty on initial conditions originated by the moving window method. It is worth remarking that all these trajectories are given the same initial position, x = 0, but different velocities. The un- certainly on this initial condition makes the entropy in- crease. However, all these trajectories at the end of the period T — 2it /lu must go back to the initial position, x = 0, thereby forcing the entropy, too, to regress to the initial vanishing value. This regression is complete if only the deterministic component is present (a = 0). This be- havior is repeated infinitely many times (see Fig. 5). It is also remarkable that the presence of an even very small stochastic component provokes an abrupt transition to a completely different behavior, where the strength of the oscillations, which are a signature of the presence of a deterministic term, tends very slowly to zero. Fig. 5 must be considered as being one of the main results of this paper. We have found that also in the case of a single sequence the analysis made by means of the DEA yields a scaling larger than that of normal diffu- sion, as an effect of the non-stationary bias, yielding the contribution Xo(t). The agreement between Eq. (p6|) and the numerical result is so good as to give a compelling support to Eq. ( f46| ) and, thus, to Eq.(p6|), from which Eq.([46|) is derived. Eq. (f46|) , on the other hand, serves very well the purpose of illustrating the joint action of the two sources of entropy increase, the trajectory ran- domness and the uncertainty on initial conditions. that reached in earlier work |5J. However, in the non- stationary case considered in this paper, it is possible to go beyond the ballistic limit within which the FBM is confined. Thus, we also establish a contact with the earlier research work of Refs. |lC|-|l^]. These results, as interesting as they are, served the main purpose of developing a proper perspective to an- alyze the diffusion process stemming from a single, non stationary time series. This is a challenging problem that, to the best of our knowledge, has been earlier addressed only in the work of Ref . S . We address this issue using the DE as a method of scaling detection, and we find that the non-stationary condition generates either anomalous scaling and ordinary statistics or ordinary scaling and anomalous statistics. We also prove that the presence of a non-stationary bias yields a process of entropy in- crease that depends on the uncertainty on initial condi- tion, as well as on trajectory randomness. This is a le- gitimate source of entropy increase that , however, might lead to the wrong conclusion that the dynamic process described by the time series under study has memory. There are other interesting aspects, such as the devia- tion from Gaussian statistics and its possible connection with the earlier work of Metzler and Klafter , on one side, and Wang ans Gaspard p6| , on the other. This is left as a subject for further investigation. IV. CONCLUSIONS The first interesting result of this paper is that concern- ing the dynamic derivation of a diffusion process overlap- ping the FBM. The FBM is a very popular generaliza- tion of ordinary Brownian motion, defined by Eq. (|l]), with 8 = 1/2. the distribution F being Gaussian. The FBM keeps the Gaussian form of F and extend S to all possible values of 6 within the interval [0, 1]. With very simple arguments we have found that this generalized form of Brownian diffusion can be realized by releasing the stationary assumption, a conclusion reminiscent of G 10 100 time FIG. 1. S(t) as a function of time t, for the sequence of section III A with a = 1, 6 = 2 • 1(T 5 , N = 10 e (□). The dashed line is the theoretical prediction obtained from eq. @. (+) is the entropy for a = 1, b = 0, N = 10 6 i.e. for the ordinary brownian diffusion. 4000 FIG. 2. p(x,t) for the set of figure gV] calculated at t = 5 • 10 3 (solid line) . The dashed line stand for the theoret- ical prediction. The dotted line is a normalizated gaussian. time FIG. 4. S(t) for the sequence of section III B with a = 1 (solid line). The dashed line is the theoretical prediction ob- tained in the asymptotic limit from eq. fl39). 10 100 lime FIG. 5. S(t) for the set of section III C with o = b = 1, T = 2ir/ui = 60 (solid line). The dashed line is the theoretical prediction obtained from eq. (^) . The dotted line is obtained setting a = 1, b = 0; the dotted-dashed one setting a = 0, 6=1. 0.001 FIG. 3. Theoretical prediction (solid line) and gaussian shape (dotted line) for p(x,t) of figure The scale for the ordinate axis is logarithmic. [1] C.-K.Peng, S.V. Buldyrev, S. Havlin, M. Simons, H.E. Stanley, and A.L. Goldenberger, Phs. Rev. E 49, 1685 (1994). [2] S.-K. Ma, Statistical Mechanics (World Scientific, Singa- pore, 1985). [3] B.V. Gnedenko, A.N. Kolmogorov, Limit Distributions for Sum of Independent Random Variables (Addison Wesley, Reading 1954). [4] B.B. Mandelbrot, The Fractal Geometry of Nature (Free- man, New York, 1983). [5] P. Allegrini, M. Buiatti, P. Grigolini and B.J. West, Phys. Rev. 57, 4558 (1998). [6] N. Scafetta, P. Hamilton, P. Grigolini, Fractals 9, 193 (2001). [7] P. Grigolini, D. Leddon, N. Scafetta, Phys. Rev. E 65, 046203 (2002). 7 [8] K. Hu, P.Ch. Ivanov, Z. Chen, P. Carpena, and H.E. Stanley, Phys. Rev. 64, 011114 (2001). [9] P. Allegrini, M. Barbi, P. Grigolini and B. J. West, Phys. Rev. E 52, 5281 (1995). [10] J. Masoliver, Phys. Rev. A 45, 706 (1992). [11] J. Masoliver, Phys. Rev. E 48, 121 (1993). [12] J. Masoliver and K.-G. Wang, Phys. Rev. E 51, 2987 (1995). [13] R. Metzler, J. Klafter, Phys. Rep. 339, 1 (2000). [14] M. Ignaccolo, P. Grigolini, A. Rosa, Phys. Rev. E 64, 026210. [15] P. Gaspard and X.J. Wang, Proc. Natl. Acad. Sci. U.S.A. 85, 4591 (1988). [16] X.-J. Wang, P. Gaspard, Phys. Rev. A 46, R3000 (1992). 8