Continuum time limit and stationary states of the minority game

Continuum time limit and stationary states of the Minority Game Matteo Marsili arXiv:cond-mat/0102257v3 [cond-mat.stat-mech] 3 Jul 2001 Istituto Nazionale per la Fisica della Materia (INFM), Trieste-SISSA Unit, Via Beirut 2-4, Trieste 34014, Italy Damien Challet Theoretical Physics, Oxford University, 1 Keble Road, Oxford OX1 3NP, United Kingdom (February 1, 2008) We discuss in detail the derivation of stochastic differential equations for the continuum time limit of the Minority Game. We show that all properties of the Minority Game can be understood by a careful theoretical analysis of such equations. In particular, i) we confirm that the stationary state properties are given by the ground state configurations of a disordered (soft) spin system; ii) we derive the full stationary state distribution; iii) we characterize the dependence on initial conditions in the symmetric phase and iv) we clarify the behavior of the system as a function of the learning rate. This leaves us with a complete and coherent picture of the collective behavior of the Minority Game. Strikingly we find that the temperature like parameter which is introduced in the choice behavior of individual agents turns out to play the role, at the collective level, of the inverse of a thermodynamic temperature. I. INTRODUCTION acterize the average behavior of agents by computing the frequency with which they play their strategies. This step can be translated in the study of the ground state properties of a soft spin disordered Hamiltonian. Secondly we characterize the fluctuations around the average behavior. To do this, we explicitly solve the Fokker-Planck equation associated to the stochastic dynamics. The new results which we derive are: Even under the most demanding definition, the Minority Game (MG) [1,2] definitely qualifies as a complex system. The MG can be regarded as an Ising model for systems of heterogeneous adaptive agents, who interact via a global mechanism that entails competition for limited resource, as found for instance in biology and financial markets. In spite of more than three years of intense research, its rich dynamical behavior is still the subject of investigations, many variations of the basic MG being proposed, each uncovering new surprising regions of phase space. Most importantly, Refs. [3–7] have shown that much theoretical insight can be gained on the behavior of this class of models, using non-equilibrium statistical physics and statistical mechanics of disordered systems. The approach of Refs. [3–6] rests on the assumption that, in a continuum time limit (CTL), the dynamics of the MG can be described by a set of deterministic equations. From these, one derives a function H which is minimized along all trajectories; hence, the stationary state of the system corresponds to the ground state of H, which can be computed exactly by statistical mechanics techniques. This approach has been challenged in Refs. [8,9], which have proposed a stochastic dynamics for the MG, thus leading to some debate in the literature [10,11]. In this paper, our aim is to analyze in detail the derivation of the CTL in order to clarify this issue. We show that a proper derivation of the CTL indeed reconciles the two approaches: the resulting dynamical equations – Eqs. (15-17) below, which are our central result – are indeed stochastic, as suggested in Ref. [8,9], but still the stationary state of the dynamics is described by the minima of the function H, as suggested in Refs. [3,4]. We then confirm the analytic results derived previously. In few words, our analysis follows two main steps: first we char- 1. we derive the full probability distribution in the stationary state. Remarkably we find that the parameter which is introduced as a temperature in the individual choice model, turns out to play the role of the inverse of a global temperature; 2. for α > αc the distribution factorizes over the agents whereas in the symmetric phase (α < αc ) agents play in a correlated way. In the latter case, the correlations contribute to the stochastic force acting on agents. We show how the dependence of global efficiency on individual temperature found in Ref. [8] arises as a consequence of these correlations; 3. we extend the analytic approach of Refs. [3,4] to the α < αc phase and asymmetric initial conditions. The dependence on the initial conditions in this phase, first noticed and discussed in Refs. [3,4], has been more recently studied quantitatively in Refs. [9,7]. We clarify the origin of this behavior and derive analytic solutions in the limit Γ → 0. 4. we show that the stronger is the initial asymmetry in agents evaluation of their strategies, the larger is the efficiency and the more stable is the system against crowd effects [12]. 5. we derive the Hamiltonian of MGs with non-linear payoffs. 1 degrees of freedom. As a consequence also Us,i (t) and hence the probability with which agents chose si (t) are subject to stochastic fluctuations. Our analysis will indeed focus on the characterization of the low-frequency fluctuations of Us,i , by integrating our the high frequency fluctuations of µ(t) and si (t). This will become clearer in the next section. For the time being let it suffice to say that there are two level of fluctuations, that of “fast” variables µ(t) and si (t) and that of “slow” degrees of freedom Us,i (t). The key parameter is the ratio α = P/N [18] and two relevant quantities are This leaves us with a coherent picture of the collective behavior of the Minority Game which is an important reference framework for the study of complex systems of heterogeneous adaptive agents. II. THE MODEL The dynamics of the MG is defined in terms of dynamical variables Us,i (t) in discrete time t = 0, 1, . . . . These are scores, propensities or “attractions” [14] which each agent i = 1, . . . , N attaches to each of his possible choices s = 1, . . . , S. Each agent takes a decision si (t) with eΓi Us,i (t) Prob{si (t) = s} = P Γ U ′ (t) i s ,i s′ e σ 2 = hA2 i, (1) N X µ(t) i=1 Our approach, which follows that of Refs. [3,4], is based on two key observations: each agent updates his scores according to µ(t) A(t) Us,i (t + 1) = Us,i (t) − as,i P . (4) III. THE CONTINUUM TIME LIMIT (2) asi (t),i P 1 X hA|µi2 P µ=1 which measure, respectively, global efficiency and predictability† . Generalizations of the model, where agents account for their market impact [3,4], where deterministic agents – so-called producers – are present [5], or where agents are allowed not to play [20–24], have been proposed. Rather than dealing with the most generic model which would depend on too many parameters, we shall limit our discussion to the plain MG. Furthermore we shall specialize, in the second part of the paper to the case S = 2 which lends itself to a simpler analytic treatment. The analysis carries through in obvious ways to the more general cases discussed in Refs. [4,5,3,24]. where Γi > 0 appears as an “individual inverse temperature”. The original MG corresponds to Γi = ∞ [1], and was generalized later to Γi ≡ Γ < ∞ [8]. The public information variable µ(t) is given to all agents; it belongs to the set of integers (1, . . . , P ), and can either be the binary encoding of last M winning choices [1], or drawn at random from a uniform distribution [15]; we stick to the latter case for sake of simplicity∗ . The acµ(t) tion asi (t),i of each agent depends the on its choice si (t) and on µ(t). The coefficients aµs,i , called strategies, plays the role of quenched disorder: they are randomly drawn signs (Prob{aµs,i = ±1} = 1/2), independently for each i, s and µ. On the basis of the outcome A(t) = H= 2 1. the scaling σ√ ∼ N , at fixed α, suggests that typically A(t) ∼ N . Hence time increments of U√ s,i (t), √ in Eq. (3) are small (i.e. of order N /P ∼ 1/ N ); (3) The idea of this equation is that agents reward [Us,i (t + 1) > Us,i (t)] those strategy which would have predicted the minority sign −A(t)/|A(t)|. The MG was initially proposed with a nonlinear dependence on A(t), i.e. with µ(t) a dynamics Us,i (t + 1) = Us,i (t) − as,i sign A(t). This leads to qualitatively similar results. The extension of our theory to non-linear cases is dealt with in the appendix. We shall not discuss any longer the interpretation of the model, which is discussed at length elsewhere [4,17–19]. The source of randomness are in the choices of µ(t) by Nature and of si (t) by agents. These are fast fluctuating 2. characteristic times of the dynamics are proportional to P . Naively this is because agents need to “test” their strategies against all P values of µ, which requires of order P time steps. More precisely, one can reach this conclusion by measuring relaxation or correlation times and verifying that they indeed grow linearly with P (see Ref. [10]). The second observation implies that one needs to study the dynamics in the rescaled time τ = t/P . This makes † Averages h. . .i stand for time averages in the stationary state of the process. Then h. . . |µi stands for time averages conditional on µ(t) = µ. ∗ Both prescriptions lead to qualitatively similar results for the quantities we study here. See [16] for more details. 2 Now the first term is of order dτ as required for a deterministic term. In addition it remains finite as N → ∞ [25]. The second term is a sum of P dτ random variables Xs,i (t) with zero average. We take dτ fixed and N very large, so that P dτ ≫ 1 and we can use limit theorems. The variables Xs,i (t) are independent from time to time, because both µ(t) and sj (t) are drawn independently at each time. Hence Xs,i (t) for P τ ≤ t < P (τ + dτ ) are independent and identically‡ distributed. For P dτ ≫ 1 we may approximate the second term dWs,i of Eq. (10) by a Gaussian variable with zero average and variance our approach differ from that of Refs. [8,9], where the time is not rescaled. In order to study the dynamics for P, N ≫ 1 at fixed α, we shall focus on a fixed small increment dτ such that P dτ = αN dτ ≫ 1. This means that we take the continuum time limit dτ → 0 only after the thermodynamic limit N → ∞. We focus only on the leading order in N . Furthermore we shall also consider Γi finite and Γi dτ ≪ 1 (5) which means that the limit Γi → ∞ should be taken after the limit dτ → 0. The orders in which these limits are taken, given the agreement with numerical simulation results, does not really matters: as we shall see differences only enter in the finite size corrections. We shall come back later to these issues. Iteration of the dynamics for P dτ time steps, from t = P τ to t = P (τ + dτ ) gives hdWs,i (τ )dWr,j (τ ′ )i = δ(τ − τ ′ ) P2 P (τ +dτ ) X t=P τ = δ(τ − τ ′ )dτ hXs,i (t)Xr,j (t)iπ hXs,i (t)Xr,j (t)iπ P (6) where the δ(τ − τ ′ ) comes from independence of Xs,i (t) and Xr,j (t′ ) for t 6= t′ and the fact that Xs,i (t) are identically distributed in time leads to the expression in the second line. Now: where we have introduced the functions us,i (τ ) = Us,i (P τ ). Let us separate a deterministic (dus,i ) from a stochastic (dWs,i ) term in this equation by replacing hXs,i (t)Xr,j (t)iπ as,i ar,j hA2 iπ as,i hAiπ ar,j hAiπ = − P P P (11) us,i (τ + dτ ) − us,i (τ ) = − 1 P P (τ +dτ )−1 X µ(t) as,i A(t). t=P τ µ(t) as,i A(t) = as,i hAiπ + Xs,i (t). The second term always vanishes for N → ∞ because as,i hAiπ is of order N 0 [25]. In the first term, instead, (7) Here and henceforth, we denote averages over µ by an over-line hA2 |µiπ = N + P 1 X µ R= R , P µ=1 1 P dτ P (τ +dτ )−1 X t=P τ eΓi Us,i (t) P Γ U (t) , i r,i re (8) which is the frequency with which agent i plays strategy s in the time interval P τ ≤ t < P (τ + dτ ). Notice that πs,i (τ ) will themselves be stochastic variables, hence we also define the average on the stationary state as Z 1 τ0 +T dτ h. . .iπ(τ ) (9) h. . .i = lim τ0 ,T →∞ T τ 0 hA2 |µiπ ≈ hA2 i ≡ σ 2 . P (τ +dτ )−1 X (13) ‡ Strictly speaking, si (t) is drawn from the distribution in Eq. (1) and not from πs,i of Eq. (8). However these two distributions differ by a negligible amount as long as the condition (5) holds (see later). us,i (τ + dτ ) − us,i (τ ) = dus,i (τ ) + dWs,i (τ ) 1 P (12) This approximation can be justified naively by observing that the dependence on πs,i of the correlations only involves the global quantity hA2 |µiπ /P for which one expects some sort of self-averaging properties. Numerical results suggest that terms which are ignored by Eq. (13) are negligible for N ≫ 1, but we were unable to prove where the average inside the integral is performed with the probabilities πs,i (τ ). Hence: = −as,i hAiπ dτ + k6=l=1 aµs′ ,k aµr′ ,l πs′ ,k πr′ ,l s′ ,r ′ is of order N which then gives a positive contribution in Eq. (11) for N → ∞. Eq. (12) leads to a stochastic dynamics where the noise covariance depends on the stochastic variables us,i (τ ) themselves. The complications which result from this fact can be avoided if one takes the approximation while h. . .iπ stands for an average over the distributions πs,i (τ ) = N X X Xs,i (t) (10) t=P τ 3 confirm Eq. (17) is valid even for α < αc . Fig. 2 compares the results of numerical simulations of the MG, as a function of Γ, with those of a semi-analytic solution of Eqs. (15-17), to be discussed later. The agreement of the two approaches shows that Eqs. (15-17) are valid even in the symmetric phase (α ≤ αc ) for all values of Γ. The instantaneous probability distribution of s, in the continuum time limit, reads this in general [26] (Eq. (13) holds trivially in the limit Γ → 0). Within this approximation, the correlation, for N ≫ 1, become hdWs,i (τ )dWr,j (τ ′ )i ∼ = σ2 as,i ar,j δ(τ − τ ′ )dτ αN (14) Note that, for r 6= s or j 6= √i, correlations hdWs,i (τ )dWr,j (τ ′ )i ∝ as,i ar,j ∼ 1/ N vanishes as N → ∞. However it is important to keep the offdiagonal terms because they keeps the dynamics of the phase space point |U (t)i = {Us,i (t)}s=1,...,S, i=1,...,N constrained to the linear space spanned by the vectors |aµ i = {aµs,i }s=1,...,S, i=1,...,N which contains the initial condition |U (0)i. The original dynamics of Us,i (t) indeed posses this property. It is important to remark that the approximation Eq. (13) makes our approach a self-consistent theory for σ 2 . We introduce σ 2 as a constant in Eq. (13) which then has to be computed self-consistently from the dynamic equations. eΓi us,i (τ ) πs,i (τ ) = P Γ u (τ ) . i r,i re This and Ito calculus then lead to a dynamic equation for πs,i (τ ). We prefer to exhibit this for Γi = Γ and using the rescaled time t = Γτ : i h dπs,i = −πs,i as,i hAi − ~πi · ~ai hAi dt
√ σ2 Γ + πs,i πs,i − ~πi2 + Γπs,i (ηs,i − ~πi ~ηi ) . (19) αN The first term in the r.h.s. comes from the deterministic part of Eq. (15), the second from the Ito term √ (where we neglected terms proportional to ai,s ai,s′ ∼ 1/ N for s 6= s′ ) and the third from the stochastic part. It is clear that, in the limit Γ → 0 the last two term vanish and the dynamics becomes deterministic. We see then that Γ tunes the strength of stochastic fluctuations in much the same way as temperature does for thermal fluctuations in statistical mechanics. The “individual inverse temperature” Γ should indeed more correctly be interpreted as a learning rate§ . Furthermore it plays the role of a “global temperature”. We shall pursue this discussion in detail below for the case S = 2. At this point, let us comment on the limit Γ → ∞, which is of particular importance, since it corresponds to the original MG. It is clear that in the limit Γ → ∞ the dynamical equations (19) become problematic. The origin of the problem lies in the order in which the limits N → ∞ and Γ → ∞ is performed. Indeed in Eq. (8) Us,i (t) ≃ us,i (τ ) + O(dτ ) for P τ ≤ t < P (τ + dτ ). Therefore, as long as Γdτ ≪ 1 the difference between Eq. (8) and πs,i in Eq. (18) is negligible. In practice, in order to satisfy both Γdτ ≪ 1 and P dτ ≫ 1, one needs Γ ≪ P . When this condition is not satisfied the instantaneous probability Eq. (1) fluctuates very rapidly at each time-step. Eq. (8) averages out these high frequency fluctuations so that, even for Γ = ∞, the distribution πs,i (τ ) of Eq. (8) is not a discontinuous step function of 10 fluctuations X αX 2 σ /N 1 0.1 0.1 1 α 10 FIG. 1. Noise strength averaged over all agents (stars); when it is multiplied by α (diamonds), one recovers σ 2 /N (circles) (P = 32, S = 2, 300P , iterations, Γ = ∞, average over 50 samples). Dashed lines are for eye guidance only Summarizing, the dynamics of us,i is described by a continuum time Langevin equation: dus,i (τ ) = −as,i hAi + ηs,i (τ ) dτ hηs,i (τ )i = 0 (18) (15) (16) 2 σ hηs,i (τ )ηr,j (τ ′ )i ∼ as,i ar,j δ(τ − τ ′ ). = αN (17) § This is an a posteriori learning rate. Indeed 1/Γ is the time the dynamics of the scores needs in order to learn a payoff difference. From a different viewpoint, Γ tunes the randomness of the response of agents. The larger the randomness, the longer it takes to average fluctuations out. Equation (15), given its derivation, has to be interpreted in the Ito sense. The expression for the noise strength by figure 1 where the measure of P is confirmed 2 i/(P N S) in a MG is reported; note that X = i,s h(Xi,s these numerical simulations were done for Γ = ∞ and 4 (22) are one and the same problem∗∗ . In other words fs,i can be computed from the constrained minimization of H as proposed in Refs. [3,4]. Hence the statistical mechanics approach based on the study of the ground state of H is correct. This approach gives the frequency fs,i with which agents play their strategies. We remark once more that H is a function of the stationary state probabilities fs,i . Also note that us,i (τ ), as suggested by Eq. (18). High frequency fluctuations contribute to the functional form of πs,i on us,i which will differ from Eq. (18). Summarizing, when we let Γ → ∞ only after the limit N → ∞ has been taken no problem arises. There is no reason to believe that results change if the order of the limits is interchanged. This expectations, as we shall see, is confirmed by numerical simulations (see Fig. 2): direct numerical simulations of the MG deviate from the prediction of Eqs. (15-17) only for finite size effects which vanish as N → ∞. Eqs. (15-17) are our central result. We shall devote the rest of the paper to discuss their content and to show that all of the observed behavior of the MG can be derived from these equations. H̃{πs,i } = Let us take the average, denoted by h. . .i, of Eq. (15) on the stationary state (SS). Let fs,i = hπs,i i be the frequency with which agent i plays strategy s in the SS. Then we have   X dUs,i = −as,i hAi, hA|µi = fs′ ,j aµs′ ,j vs,i ≡ dτ ′ j,s A. Independence on Γ for α > αc Given the relation between πs,i and Us,i and considering that the long time dynamics of Us,i in the SS is Us,i (τ ) = const+vs,i τ , we have that i) each strategy which is played in the SS by agent i must have the same “velocity” vs,i = vi∗ , and ii) strategies which are not played (i.e. with fs,i = 0) must have vs,i < vi∗ . In other words −as,i hAi ≤ vi∗ , When α > αc the solution to Eq. (22) is unique and H > 0. Hence fs,i does not depends on Γ, neither does H. In addition we shall see that hπs,i πs′ ,j i = hπs,i ihπs′ ,j i = fs,i fs′ ,j min {fs,i ≥0} for i 6= j (23) implying that ∀i, s such that fs,i > 0 ∀i, s such that fs,i = 0. (20) σ2 ≡ N + (21) Consider now the problem of constrained minimization of H in Eq. (4), subject to fs,i ≥ 0 for all s, i and the normalization conditions. Introducing Lagrange multipliers P λi to enforce s fs,i = 1 for all i, this problem reads ( as,i as′ ,j πs,i πs′ ,j i,j=1 s,s′ =1 as a function of the instantaneous probabilities πs,i , is not a Lyapunov function of the dynamics. The dynamical variables πs,i √ (t) are subject to stochastic fluctuations of the order of Γi around their average values fs,i . Only in the limit Γi → 0, when the dynamics becomes deterministic and πs,i → fs,i , the quantity H̃{πs,i } becomes a Lyapunov function. The solution to the minimization of H reveals two qualitatively distinct phases [3,4] which are separated by a phase transition occurring as α → αc . We discuss qualitatively the behavior of the solution for a generic S and leave for the next section a more detailed discussion in the simpler case S = 2. IV. STATIONARY STATE − as,i hAi = vi∗ , N S X X hAi2 − N X i=1 λi 1− S X s=1 fs,i !) . XX i6=j s,r as,i ar,j hπs,i πr,j i does not depend on Γ either. Hence the solution {fs,i } uniquely determines all quantities in the SS, as well as ∗∗ Indeed both problems can be put in the form of a Linear Complementarity problem [27]: (22) Taking derivatives, we find that if fs,i > 0 then as,i hAi + λi = 0 whereas if fs,i = 0 then as,i hAi + λi ≥ 0. These are exactly Eqs. (20,21) where vi∗ = λi . We then conclude that the two problems Eqs. (20,21) and Eq. X j,s′ fs,i " X j,s′ fs,i ≥ 0 as,i as′ ,j fs′ ,j + vi∗ ≥ 0 as,i as′ ,j fs′ ,j + vi∗ # =0 This problem has a solution for all values of vi∗ because of non-negativity of the matrix as,i as′ ,j , see Ref. [27]. 5 the parameters which enter into the dynamics (notice the dependence on σ 2 in Eq. 17). In particular σ 2 does not depend on Γ. V. THE CASE S = 2 We work in this section with the simpler case of S = 2 strategies, labelled by s = ±. We also set Γi = Γ for all i. Following Refs. [3,13] we introduce the variables B. Dependence on Γ and on initial conditions for α < αc ξiµ = For α < αc the solution to the minimization problem is not unique: there is a connected set of points {fs,i } such that H = 0. Let us first discuss the behavior of the system in the limit Γ → 0, where the dynamics becomes deterministic. The dynamics reaches a stationary state {fs,i } which depends on the initial conditions. In order to see this, let us introduce the vector notation |vi = {vs,i , s = 1, . . . , S, i = 1, . . . , N }. Then for all times |u(τ )i is of the form |u(τ )i = |u(0)i + P X µ=1 µ aµ+,i − aµ−,i , 2 Ωµ = N X aµ+,i + aµ−,i i=1 2 . Let us rescale time t = Γτ and introduce the variables yi (t) = Γ U+,i (τ ) − U−,i (τ ) 2 Then, using Eq. (8), the dynamical equations (15-17) become N X dyi = −ξi Ω − ξi ξj tanh(yj ) + ζi dt j=1 µ |a iC (τ ) hζi (t)ζj (t′ )i = where C µ (τ ) are P functions of time. If there are vectors hv| such that hv|aµ i = 0 for all µ, then hv|u(τ )i = hv|u(0)i, i.e. the components of the scores will not change at all along these vectors. As a result the SS will depend on initial conditions |u(0)i. These vectors hv| exist exactly for α < αc [4], because the “dimensionality” of the vectors |u(τ )i is larger than P †† . The picture is made even more complex by the fact that, for α < αc , when Γ is finite Eq. (23) does not hold. Hence σ 2 has a contribution which depends on the stochastic fluctuations around fs,i . The strength of these fluctuations, by Eqs. (17,19), depends on Γ and σ 2 itself. We face, in this case, a self-consistent problem: σ 2 enters as a parameter of the dynamics but should be computed in the stationary state of the dynamics itself. Therefore the solution to this problem and hence σ 2 depends on Γ. The solution {fs,i } to the minimization of H should also be computed self-consistently. As a result, the SS properties acquire a dependence on Γ. The condition (23), which is similar to the clustering property in spin glasses [28], plays then a crucial role. We show below how the condition (23), the dependence on initial conditions and on Γ enter into the detailed solution for S = 2. By similar arguments our conclusion can be generalized to all S > 2. Γσ 2 ξi ξj δ(t − t′ ). αN (24) (25) The Fokker-Planck (FP) equation for the probability distribution P ({yi }, t) under this dynamics reads  N N X ∂P ({yi }, t) X ∂  ξi Ω + ξi ξj tanh(yj )+ = ∂t ∂yi  j=1 i=1  N 1X ∂  ξi ξj P ({yi }, t) (26) β ∂yj  j=1 where we have introduced the parameter β= 2αN . Γσ 2 (27) Multiply Eq. (26) by yi and integrate over all variables. Using integration by parts, assuming that P → 0 fast as yj → ∞, one gets N X ∂ ξi ξj htanh(yj )i . hyi i = −ξi Ω − ∂t j=1 Let us look for solutions with hyi i ∼ vi t and define mi = htanh(yi )i. Hence for t → ∞ we have vi = −ξi Ω − †† In order to compute the dimensionality of the vectors |ui we have to take into account the N normalization conditions and the fact that strategies which are not played (fs,i = 0) should not be counted. So if there are N> variables fs,i > 0, the relevant dimension of the space of |ui is N> − N . Hence vectors hv| orthogonal to all |aµ i exist for N> − N > P , i.e. for α < αc = N> (αc )/N − 1. N X ξi ξj mj . (28) j=1 Now, either vi = 0 and hyi i is finite, or vi 6= 0, which means that yi → ±∞ and mi = sign vi . In the latter case (vi 6= 0) we say that agent i is frozen [13], we call F the set of frozen agents and φ = |F |/N the fraction of frozen agents. 6 As in the general case, the parameters vi for i ∈ F and mi ≡ htanh(yi )i for i 6∈ F are obtained by solving the constrained minimization of #2 " N P X 1 X ξiµ mi . Ωµ + H= P µ=1 i=1 A. α > αc For α > αc the solution of min H is unique, and hence mi depends only on the realization of the disorder. In addition, the number N − |F | ≡ N (1 − φ) of unfrozen agents is less than P and the constraint is ineffective, i.e. Py(0) ≡ 1. The scores |yi of unfrozen agents span a linear space which is embedded in the one spanned by the vectors |ξ µ i. Hence the dependence on initial conditions yi (0) drops out. Therefore the probability distribution of yi factorizes, as in Eq. (29). Then the third term of Eq. (32), which is the only one which depends on β, vanishes identically. We conclude that, for α > αc , σ 2 only depends on mi and is hence independent of Γ as confirmed by numerical simulations. Summarizing, for α > αc one derives a complete solution of the MG by finding first the minimum {mi } of H and then by computing σ 2 , β and the full distribution of yi from Eq. (29). When the solution of min H is unique, i.e. for α > αc , the parameters mi depend only on the realization of disorder {ξiµ , Ωµ }, and their distribution can be computed as in Ref [3]. When the solution is not unique, i.e. for α < αc , we are left with the problem of finding which solution the dynamics selects. Let us suppose that we have solved this problem (we shall come back later to this issue), so that all mi are known. Using the stationary condition Eq. (28), we can write the FP equation for the probability distribution Pu (yi , i 6∈ F) of unfrozen agents. For times so large that all agents in F are indeed frozen (i.e. si (t) = sign vi ) this reads:   X ∂ X ∂Pu 1 ∂ Pu . = ξi ξj tanh(yj ) − mj + ∂t ∂yi β ∂yj B. α ≤ αc : dependence on Γ and crowd effects j6∈F i6∈F This has a solution     X Pu ∝ exp −β [log cosh yj − mj yj ] .   When α < 1 − φ, on the other hand, the solution of min H is not unique. Furthermore the constraint cannot be integrated out and the stationary distribution depends on the initial conditions. Numerical simulations [8] show that σ 2 increases with Γ for α < αc (see Fig. 2). This effect has been related to crowd effects in financial markets [12]. Ref. [6] has shown that crowd effects can be fully understood in the limit α → 0: as Γ exceeds a critical learning rate Γc , the time independent SS becomes unstable and a bifurcation to a period two orbit occurs. Neglecting the stochastic term ζi , Ref. [6] shows that this picture can be extended to α > 0‡‡ . This approach suggests a crossover to a (29) j6∈F Finally we have to impose the constraint that |y(t)i = {yi (t)}N i=1 must lie on the linear space spanned by the vectors |ξ µ i which contains the initial condition |y(0)i. This means that     X Pu ∝ Py(0) exp −β [log cosh yj − mj yj ] (30)   j6∈F where the projector Py(0) is given by # " P P Z ∞ N X Y Y µ µ µ δ yi − yi (0) − c ξi . (31) Py(0) ≡ dc −∞ µ=1 i=1 µ=1 ‡‡ The idea of Ref. [6] is the following: imagine that our system is close to a SS point yi∗ at time t = tk , when µ(t) = 1. Will the system be close to yi∗ when the pattern µ = 1 occurs again the next time t′ = tk+1 ? To see this, let us integrate Eq. (24) from t = tk to tk+1 . In doing this we i) neglect the noise term (i.e. ζi = 0) and ii) assume that tanh yi (t) ≈ tanh yi (tk ) stays constant in the integration time interval. This latter assumption is similar to the recently introduced [7] batch version of the MG, where agents update their strategies every P time-steps. This leads to study a discrete time dynamical system We find it remarkable that Γ, which is introduced as the inverse of an individual “temperature” in the definition of the model, actually turns out to be proportional to β −1 (see Eq. 27) which plays collectively a role quite similar to that of temperature. Using the distribution Eq. (30), we can compute σ2 = H + N X i=1 X i6=j ξi2 (1 − m2i ) + ξi ξj h(tanh yi − mi )(tanh yj − mj )i. (32) " yi (tk+1 ) = yi (tk ) − Γ Ωξi + This depends on β, i.e. on σ 2 itself by virtue of Eq. (27). The stationary state is then the solution of a selfconsistent problem. Let us analyze in detail the solution of this self-consistent problem. N X j=1 # ξi ξj tanh[yj (tk )] (33) where the factor Γ comes because tk+1 −tk is on average equal to Γ. The linear stability of fixed point solutions is analyzed setting yi (tk ) = yi∗ + δyi (k) and computing the eigenvalues of 7 N N 1 X 1 X hsi (t)si (0)i = mi t→∞ N N i=1 i=1 “turbulent” dynamics for Γ > Γc , where 4 √ 2 Γc (α) = (1 + α) (1 − Q) M ≡ lim (34) which measures the overlap of the SS configuration with the initial condition. Symmetric initial conditions are related to M = 0 SS. These are the states we focus on. The solution is derived in two steps: and N 1 X 2 Q= m . N i=1 i 1. find the minimum {mi } of H, with M = 0; Both Q and Γc can be computed exactly in the limit N → ∞ within the statistical mechanics approach [3,6]. This approach however i) does not properly takes into account the stochastic term, ii) does not explain what happens for Γ > Γc and iii) does not explains why such effects occur only for α < αc . 2. compute self-consistently σ 2 . The numerical procedure for solving the problem is the following: given the realization of disorder {ξiµ , Ωµ }, step (1) — finding the minimum {mi } of H — is straightforward. For step (2) we sample the distribution Eq. (30) with the Montecarlo method§§ at inverse temperature β and measure the β dependent contribution of σ 2 in Eq. (32): X Σ(β) = ξi ξj h(tanh yi − mi )(tanh yj − mj )iβ . 4 N=160 MG N=320 MG N=640 MG N=160 Montecarlo N=320 Montecarlo N=640 Montecarlo Theory Γ −> 0 2 i6=j Here h. . .iβ stands for an average over the distribution Eq. (30) with parameter β. Finally we solve the equation   2αN 2 2 σ (Γ) = σ (0) + Σ . (36) Γσ 2 (Γ) 0.8 2 σ (Γ)/Ν 3 Γc 0.6 1 0.4 0 −2 10 10 −1 (35) 10 0 Γ 0 0.1 0.2 10 1 0.3 10 This procedure was carried out for different system sizes and several values of Γ. The results, shown in Fig. 2, agree perfectly well with direct numerical simulations of the MG. Actually Fig. 2 shows that, for Γ ≫ 1, the solution of the self-consistent equation (36) suffers much less of finite size effects than the direct numerical simulations of the MG. Fig. 2 also shows that, even if only approximate, Eq. (34) provides an useful estimate of the point where the crossover occurs. It is possible to compute σ 2 (Γ) to leading order in Γ ≪ 1. The calculation is carried out in the appendix in detail. The result is   1 − Q + α(1 − 3Q) σ2 ∼ 1 − Q 2 1+ Γ + O(Γ ) . (37) = N 2 4α 2 FIG. 2. Global efficiency σ 2 /N versus Γ for α = 0.1 < αc and different system sizes. Lines refer to direct simulations of the MG with N = 160, 320 and 640. Finite size effect for Γ ≫ 1 are evident. Symbols refer instead to the solution of the self-consistent equation (36) for the same system sizes. For both methods and all values of N , σ 2 is averaged over 100 realizations of the disorder. The arrow marks the location of Γc predicted by Eq. (34). In the inset, the theoretical prediction Eq. (37) on the leading behavior of σ 2 /N for Γ ≪ 1 (solid line) is tested against numerical simulations of the MG (points) for the same values of N . The stochastic dynamics derived previously gives detailed answers to all these issues. We first restrict attention to symmetric initial conditions yi (0) = 0 and then discuss the dependence on initial conditions. The choice of yi (0) = 0 ∀i is convenient because it allows one to use this same symmetry to identify the solution {mi } to the minimization of H, independently of Γ. To be more precise, one can introduce a “magnetization” §§ The Montecarlo procedure follows the usual basic steps: i) A move yi → yi + ǫξiµ is proposed, with µ and ǫ drawn at random, ii) the “energy” E{yi } = N X i=1 [log cosh yi − mi yi ] of the new configuration is computed and iii) The move is accepted with a probability equal to min(1, e−β∆E ) where ∆E is the “energy” difference. P the linearized map δyi (k+1) ≃ j Ti,j δyj (k). There is a critical value of Γ above which the solution yi∗ become unstable, which is given by Eq. (34). 8 and χ = β(Q − q)/α is a “spin susceptibility”. There are two possible solutions: one with χ < ∞ finite as η → 0 which describes the α > αc phase. The other has χ ∼ 1/η which diverges as η → 0. This solution describes the α < αc phase. We focus on this second solution, which can be conveniently parameterized by two parameters z0 and ǫ0 . We find  if z ≤ −z0 − ǫ0  −1 z+ǫ0 if −z s0 (z) = 0 − ǫ0 < z < z0 − ǫ0  z0 1 if z ≥ z0 − ǫ0 The inset of Fig. 2 shows that this expression indeed reproduces quite accurately the small Γ behavior of σ 2 . Note finally that Eq. (36) has a finite solution σ 2 (∞) = σ 2 (0) + Σ(0) in the limit Γ → ∞. Furthermore it is easy to understand the origin of the behavior σ 2 /N ∼ 1/α for Eq. (36). Because of the constraint, when ξi ξj is positive (negative) the fluctuations of tanh(yi )− mi are positively (negatively) correlated with tanh(yj ) − mj . If we assume that h[tanh(yi ) − mi ][tanh(yj ) − mj ]i ≃ cξi ξj for some P 2 constant c, we find Σ ≃ c i6=j ξi ξj . This leads easily to Σ/N ≃ c/(4α), which explains the divergence of σ 2 /N as α → 0 for Γ ≫ 1. Indeed Eq. (38) p gives Q(z0 , ǫ0 )R and Eq. (39) which for χ → ∞ reads α(1 + Q) = Dzzs0 (z), then gives α(z0 , ǫ0 ). With ǫ0 6= 0 one finds solutions with a non-zero “magnetization” M = hsi i. This quantity is particularly meaningful, in this context, because it measures the overlap of the behavior of agents in the SS with their a priori preferred strategies C. Selection of different initial conditions in the Replica calculation As discussed above, the stationary state properties of the MG in the symmetric phase depend on the initial conditions. Can the statistical mechanics approach to the MG [3,4] be extended to characterize this dependence for Γ ≪ 1? If this is possible, how do we expect the resulting picture to change when Γ increases? We shall first focus on the first question (i.e. Γ ≪ 1) and then discuss the second. Of course one can introduce the constraint on the distribution of yi in the replica approach in a straightforward manner. This leads however to tedious calculations. We prefer to follow a different approach. In the α < αc phase the minimum of H is degenerate, i.e. H = 0 occurs on a connected set of points. Each of these points corresponds to a different set of initial conditions, as discussed above. In order to select P a particular point with H = 0 we can add a potential η i (si − s∗i )2 /2 to the Hamiltonian H, which will favor the solutions closer to s∗i , and then let the strength η of the potential go to zero. This procedure lifts the degeneracy and gives us the statistical features of the equilibrium close to s∗i . The nature of the stationary state changes as the asymmetry in the initial conditions changes. If we take s∗i = s∗ , the state at s∗ = 0 describes symmetric initial conditions and increasing s∗ > 0 gives asymmetric states. The saddle point equations of the statistical mechanics approach of Ref. [3] can be reduced to two equations: Z ∞ Q= Dzs20 (z) (38) −∞ Z ∞ 1+χ Dzzs0 (z) (39) χ= p α(1 + Q) −∞ N 1 X hsi (t)si (0)i. Dz s0 (z) = lim M≡ t→∞ N −∞ i=1 Z ∞ (41) Note indeed that one can always perform a “gauge” transformation in order to redefine s = +1 as the initially preferred strategy. This amounts to taking yi (0) ≥ 0 for all i. Which SS is reached from a particular initial condition is, of course, a quite complex issue which requires the integration of the dynamics. However, the relation between Q and M derived analytically from Eqs. (38,41) can easily be checked by numerical simulations of the MG. Figure 3 shows that the self-overlap Q and the magnetization M computed in numerical simulations with initial conditions yi (0) = y0 for all i, perfectly match the analytic results. The inset of this figure shows how the final magnetization M and the self-overlap Q in the SS depend on the asymmetry y0 of initial conditions. 1 1 M Q 0.8 0.5 Q 0.6 0 −1 0 10 10 y 0 0.4 1 2 10 10 2 where Dz = √dz e−z /2 , s0 (z) ∈ [−1, 1] is the value of s 2π which minimizes r 1 1+Q 1 2 s − zs + η(1 + χ)(s − s∗ )2 (40) Vz (s) = 2 α 2 0.2 0 0 0.2 0.4 M 9 0.6 0.8 FIG. 3. Relation between Q and M , for α = 0.1, derived from analytic calculation (full line) and from numerical simulations of the MG with different initial conditions y0 (⋄, P = 32, N = 320, Γ = 0.1). The inset shows the dependence of Q and M from the initial condition y0 . D. The maximally magnetized stationary state The maximally magnetized SS (MMSS), obtained in the limit y0 → ∞, is also the one with the largest value of Q, and hence with the smallest value of σ 2 = N (1 − Q)/2. σ 2 /N is plotted against α both for symmetric y0 = 0 initial conditions and for maximally asymmetric ones y0 → ∞ in fig. 5. The inset shows the behavior of Q and M in the MMSS. Let us finally discuss the dependence on Γ for asymmetric initial conditions. Eq. (34) provides a characteristic value of Γ as a function of α and Q. This theoretical prediction is tested against numerical simulations of the MG in Fig. 4: when plotted against Γ/Γc , the curves of σ 2 /N obtained from numerical simulations, approximately collapse one onto the other in the large Γ region. Fig. 4 suggests that Γc in Eq. (34) provides a close lower bound for the onset of saturation to a constant σ 2 for large Γ. We find it remarkable that a formula such as Eq. (34) which is computed in the limit Γ → 0, is able to predict the large Γ behavior. 1.00 1 0.8 0.6 0.4 2 σ /Ν 0 2 Γc= 0.260 Γc= 0.447 Γc= 0.757 Γc=1.382 Γc=2.099 0.00 −2 10 0.1 0.2 α 0.3 −1 10 0 α 10 1 10 FIG. 5. σ 2 /N for the MG with initial conditions y0 = 0 (full line) and y0 → ∞ (dashed line). The inset reports the behavior of M and Q in the y0 → ∞ SS. Note that Q is linear in α. 1 0 −1 10 0 0.50 y0=0 M=Mmax 2 σ /Ν 3 M(α) Q(α) 0.2 0 10 1 Γ/Γc 10 Remarkably we find that σ 2 /N vanishes linearly with α in the MMSS∗∗∗ . This means that, at fixed P , as N increases the fluctuation σ 2 remains constant. This contrast with what happens in the y0 = 0 state, for Γ ≪ Γc , where σ 2 increases linearly with N , and with the case Γ ≫ Γc where σ 2 ∝ N 2 [18,6]. Note also that the lowest curve of Fig. 5 also gives an upper bound to the σ 2 of Nash equilibria (see Refs. [3,4,29]). The MMSS is also the most stable state against crowd effects: if we put Q(α, y0 = ∞) ∼ = 1 − cα, as appropriate for the MMSS we find that Γc ∼ 1/α diverges with α. 2 10 FIG. 4. Global efficiency σ 2 /N versus Γ/Γc for α = 0.1 < αc , N = 160 agents and different initial conditions y0 . The value Γc is computed from Eq. (34) and is shown in the legend. With respect to the dependence on initial conditions, we observe that Γc is an increasing function of Q and hence it increases with the asymmetry y0 of initial conditions. Hence, for a fixed Γ, the fluctuation dependent part Σ of σ 2 decreases with y0 because it is an increasing function of Γ/Γc . This effects adds up to the decrease of the Γ independent part of σ 2 discussed previously. Fig. 4 also shows that the Γ ≫ Γc state is independent of initial conditions y0 . This can naively be understood observing that stochastic fluctuations induce fluctuations δyi which increase with Γ. For Γ ≫ 1 the asymmetry y0 of initial conditions is small compared to stochastic fluctuations δyi and hence the system behaves as if y0 ≈ 0. VI. CONCLUSIONS We have clarified the correct derivation of continuous time dynamics for the MG. This on the one hand reconciles the two current approaches [3,9]. On the other it leads to a complete understanding of the collective behavior of the MG. We confirm that stationary states are characterized by the minimum of a Hamiltonian which ∗∗∗ 10 This results was also found analytically in [7]   P Z √ 1 X ∞ dx −x2 /2 X µ √ e Hg = fs,i as,i + Dx G P µ=1 −∞ 2π i,s measures the predictability of the game. For α > αc we find a complete analytic solution, whereas for α < αc the statistical mechanics approach of Ref. [3] is valid for Γ → 0. It is in principle possible to introduce the new elements discussed here in the approach of Ref. [3] and to derive a full analytic solution. We have indeed derived the first term of the series expansion for Γ ≪ 1, which agrees perfectly with numerical data. The extension of the approach of Ref. [3] involves lengthy calculations and it shall be pursued elsewhere. Finally we note that the results derived in this paper generalize to more complex models. It is worth to remark that the solution to the FP equation is no more factorizable, in general, when agents account for their market impact as in Refs. [3,6,4]. Hence, as long as there are unfrozen agents, we expect that the stationary state depends on Γ. However, when the agents take fully into account their market impact, all of them are frozen and the conclusion that agents converge to Nash equilibria remains valid. We acknowledge constructive discussions with A.C.C. Coolen, A. Engel, J.P. Garrahan, J.A.F. Heimel and D. Sherrington. with g(x) = dG(x) dx and D = σ 2 − H which must be determined selfconsistently. Indeed taking the derivative of Hg w.r.t. fs,i and imposing the constraint fs,i ≥ 0 and normalization, we arrive at exactly the same equations which describe the stationary state of the process. The Hamiltonian for the original MG is derived setting g(x) = sign x, which leads to Hsign   P  1 X 1 −hA|µi2 /D hA|µi hA|µi √ √ √ = . e + erf P µ=1 π D D vi = −as,i hg(A)i, if fs,i > 0 The analysis of stochastic fluctuations can be extended to non-linear cases in a straightforward manner. Again the key point is that the dynamics is constrained to the linear space spanned by the vectors |aµ i. For α > αc we have no dependence on initial conditions. However it is not easy to show in general that the distribution of scores factorizes across agents. This means that there may be a contribution of fluctuations to σ 2 – i.e. Σ > 0 – so we cannot rule out a dependence of σ 2 on Γ. Numerical simulations for g(x) = sign x show that such a dependence, if it exists, is very weak. Anyway even though σ 2 only depend on fs,i , the minimization problem depends on D = σ 2 − H which must then be determined self-consistently. For α < αc the dependence on initial conditions induces a correlation of scores across agents. As a result σ 2 depends on Γ just as in the linear case discussed above. vi > −as,i hg(A)i, if fs,i = 0. APPENDIX B: SMALL Γ EXPANSION APPENDIX A: NON-LINEAR MINORITY GAMES Take a generic dynamics µ(t) Us,i (t + 1) = Us,i (t) − as,i g [A(t)] with g(x) some function. When we carry out the limit to continuous time we find a deterministic term which is proportional to −as,i hg(A)i. The stationary state conditions then read and For any fixed µ, A(t) is well approximated by a Gaussian variable with mean X fs,i aµs,i hA|µi = For Γ ≪ 1 it is appropriate to consider β ≫ 1 and to take zi yi = arc tanh mi + √ β i,s and variance D = σ 2 − H. Here we neglect dependences on µ. Also we treat D as a parameter and neglect its dependence on the stationary state probabilities fs,i . Hence we can write Z ∞  √  2 dx √ e−x /2 g hA|µi + Dx hg(A)|µi = 2π −∞ so that β[log cosh yi − mi yi ] ≃ 21 (1 − m2i )zi2 + O(β −1/2 ). Hence we have to sample a distribution P 2 2 1 P {zi } ∝ e− 2 i (1−mi )zi where zi has the form The stationary state conditions of the dynamics above can again be written as a minimization problem of the functional zi = P X µ=1 11 cµ ξiµ . It is convenient to express everything in terms of the coefficients cµ . Their pdf is derived from that of zi and it reads: P {cµ } ∝ e − 21 P µ,ν cµ T µ,ν cν T µ,ν = , N X i=1 [5] Challet D., M. Marsili and Y.-C. Zhang, Physica A 276, 284 (2000) e-print cond-mat/9909265 [6] M. Marsili and D. Challet, Adv Compl. Sys. 3-I (2001), e-print cond-mat/0004376 e-print cond-mat/0004196 [7] J. A. F. Heimel, A. C. C. Coolen, Phys. Rev. E 63, 056121 (2001), e-print cond-mat/0012045. [8] A. Cavagna, J.P. Garrahan, I. Giardina, D. Sherrington, Phys. Rev. Lett. 83, 4429 (1999) e-print condmat/9903415 (1999). [9] J. P. Garrahan, E. Moro, D. Sherrington, Phys. Rev. E 62, R9 (2000) [10] D. Challet, M. Marsili and R. Zecchina, Phys. Rev. Lett. 85, 5008 (2000), e-print cond-mat/0004308 [11] A. Cavagna, J. P. Garrahan, I. Giardina, D. Sherrington, Phys. Rev. Lett. 85, 5009 (2000), e-print condmat/0005134 [12] M. Hart, P. Jefferies, P. Hui, N. F. Johnson, e-print condmat/0008385 [13] D. Challet and M. Marsili, Phys. Rev. E 60, 6271 (1999), e-print cond-mat/9904071 [14] The learning process belongs to the class of models described in C. Camerer and T.-H. Ho, Econometrica 67, 827 (1999) with φ = κ = δ = 1. [15] A. Cavagna, Phys. Rev. E 59, R3783 (1999) [16] D. Challet and M. Marsili, Phys Rev. E 62, 1862 (2000), [17] Y.-C. Zhang, Evolving models of financial markets, Europhys. News 29, 51 (1998). [18] Savit R., Manuca R., and Riolo R., Phys. Rev. Lett. 82(10), 2203 (1999), e-print adap-org/9712006. [19] D. Challet, A. Chessa, M. Marsili, Y.-C. Zhang, Quant. Fin. I-1 (2001), e-print cond-mat/0011042 [20] F. Slanina and Y.-C. Zhang, Physica A 272, 257 (1999), e-print cond-mat/9906249 [21] P. Jefferies, M.L. Hart, P.M. Hui, N.F. Johnson, Int. J. Th and Appl. Fin. 3-3 (2000), e-print cond-mat/9910072. [22] T. S. Lo, P. M. Hui, N. F. Johnson, e-print condmat/0008387 [23] J.-P. Bouchaud, I. Giardina, M. Mezard, e-print condmat/0012156 [24] D. Challet, M. Marsili and Y.-C. Zhang, Physica A 294, 514 (2001), e-print cond-mat/0101326 √ [25] Indeed aµs,i hA|µiπ is of order N but its sign fluctuates. When we average over P ∼ N different µ, we get a quantity of order N 0 . We are implicitly assuming that hAi ≃ 0 which definitely holds for large times. [26] We thank J. A. F. Heimel for pointing out this issue. [27] Katta G. Murty, Linear Complementarity, Linear and Nonlinear Programming, Heldermann Verlag, (1988) [28] M. Mezard, G. Parisi, M. A. Virasoro, Spin glass theory and beyond World Scientific (1987). [29] A. De Martino and M. Marsili (2000), e-print condmat/0007397 (1 − m2i )ξiµ ξiν . From this we find hcµ cν i = [T −1 ]µ,ν . Now we split the term Σ(β) in two contributions, " Σ(β) = h N X i=1 N X i=1 #2 ξi (tanh yi − mi ) i+   ξi2 m2i − (tanh yi )2 and work them out separately. For the first we use N X 1 X µ,ν ν ξiµ (tanh yi − mi ) = √ T c β ν i=1 so that " h N X i=1 #2 ξi (tanh yi − mi ) i= = 1 X µ,ν µ,γ ν γ T T hc c i βP µ,ν,γ TrT ∼ 1 − Q N. = βP 2β Within the approximation (1 − 3m2i ) ≈ (1 − 3Q) we are able to derive a closed expression also for the second term: N X N   1X 2 ξi2 m2i −(tanh yi )2 = ξi (1−m2i )(1−3m2i )hzi2 i β i=1 i=1 ≈ α 1 − 3Q N. β 2 Hence we find   1 − Q 1 − 3Q 1 ∼ Σ(β) = + α + O(β −2 ). 2 2 β This and Eq. (27) lead to Eq. (37). [1] Challet D. and Zhang Y.-C., Physica A 246, 407 (1997) e-print adap-org/9708006 [2] See the Minority Game’s web page on http://www.unifr.ch/econophysics. [3] D. Challet, M. Marsili and R. Zecchina, Phys. Rev. Lett. 84, 1824 (2000), e-print cond-mat/9904392 [4] M. Marsili, D. Challet and R. Zecchina, Physica A 280, 522 (2000), e-print cond-mat/9908480 12