Simulation of a generalized asset exchange model with economic growth and wealth distribution
Kang K. L. Liu, N. Lubbers, W. Klein, J. Tobochnik, B. M. Boghosian, Harvey Gould
SSimulation of a generalized asset exchange model with economicgrowth and wealth distribution
Kang K. L. Liu, N. Lubbers,
2, 1
W. Klein,
1, 3
J.Tobochnik,
4, 1
B. M. Boghosian, and Harvey Gould
1, 6, ∗ Department of Physics, Boston University, Boston, Massachusetts 02215 Los Alamos National Laboratory, Los Alamos, New Mexico 87545 Center for Computational Science,Boston University, Boston, Massachusetts 02215 Department of Physics, Kalamazoo College, Kalamazoo, Michigan 49006 Department of Mathematics, Tufts University, Medford, Massachusetts 02155 Department of Physics, Clark University, Worcester, Massachusetts 01610 (Dated: February 3, 2021) a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b bstract The agent-based Yard-Sale model of wealth inequality is generalized to incorporate exponentialeconomic growth and its distribution. The distribution of economic growth is nonuniform andis determined by the wealth of each agent and a parameter λ . Our numerical results indicatethat the model has a critical point at λ = 1 between a phase for λ < λ ≥ λ <
1. Our estimates of various criticalexponents are consistent with a mean-field theory (see following paper). The exponents do notobey the usual scaling laws unless a combination of parameters that we refer to as the Ginzburgparameter is held fixed as the transition is approached. The model illustrates that both poorerand richer agents benefit from economic growth if its distribution does not favor the richer agentstoo strongly. This work and the accompanying theory paper contribute to understanding whetherthe methods of equilibrium statistical mechanics can be applied to economic systems.
I. INTRODUCTION AND THE GED MODEL
Although economies are complex systems that are difficult to understand [1], the consid-eration of simple models can provide insight if the questions are of a statistical nature andabout the economy as a whole. One such question is whether economic growth can benefitall members of society [1–3]. Another question is to what degree is wealth accumulation anatural consequence of the way that wealth is exchanged and distributed. Whether thesequestions and others can be treated using methods appropriate to equilibrium systems isnot settled [4].In this paper we approach these questions by simulating an agent-based model that incor-porates wealth exchange, economic growth, and its distribution. Agent-based asset-exchangemodels have been useful for understanding the statics and dynamics of wealth distributionsand exhibit a rich phenomenology [5–12]. The agent-based asset-exchange model we gener-alize belongs to a class of wealth exchange models that has been of considerable interest tophysicists and has provided insight into how economies function [13–23]. In these models ∗ [email protected] f of thewealth of the poorer agent is transferred to the winning agent. The latter is determinedat random with probability 1/2. If an equal amount of wealth is initially assigned to the N agents, the result is that after many exchanges, the wealth is concentrated among in-creasingly fewer agents, culminating, in the limit of infinitely many wealth exchanges and N → ∞ , to a single agent holding a fraction of the total wealth [5, 9, 14, 15, 17, 21–23, 25].To investigate the effect of economic growth on the wealth distribution, we generalize theYard-Sale model so that the wealth µW ( t ) is added to the system after N exchanges, where W ( t ) is the total wealth and the time t is defined such that one unit of time corresponds to N exchanges; W ( t ) grows exponentially with the rate µ >
0. The motivation for this typeof growth is the expected annual increase in the gross domestic product.To distribute the increased total wealth due to growth, ∆ W ( t ) = µW ( t ), to the individualagents, we introduce the distribution parameter λ ≥ i according to their wealth w i ( t ) at time t as∆ w i ( t ) = ∆ W ( t ) w λi ( t ) (cid:80) Ni =1 w λi ( t ) . (1)The form of Eq. (1) implies that as λ → − , the allocation of the added wealth is weightedmore toward agents with greater wealth. (The limit λ → λ -dependent distribution mechanism in Eq. (1) as the Growth,Exchange, and Distribution (GED) model.The motivation for the form of the wealth distribution in Eq. (1) is that in practice, notall agents benefit in the same way from economic growth, and that agents with more assetsand resources are able to take more advantage of the growth of the economy. We argue inthe Appendix that the allocation of growth according to Eq. (1) is consistent with economicdata.The distribution parameter λ , exchange parameter f , growth parameter µ , and the num-ber of agents N determine the wealth distribution in the model. Our primary results can begrouped into two categories – the implications for economic systems and the implications3or our understanding of the statistical mechanics of systems that are near-mean-field. Wefind that there is a phase transition at λ = 1 such that for λ <
1, all agents benefit fromeconomic growth, there is economic mobility, and the system is in thermodynamic equilib-rium. In contrast, for λ ≥
1, the system is non-stationary, there is no economic mobility, andthere is wealth condensation as in the original Yard-Sale model. In the context of statisticalmechanics we note that we can define an energy that satisfies the Boltzmann distribution.The existence of the latter is consistent with the assumption that that the system is not justin a steady state, but is in thermodynamic equilibrium for λ < λ <
1. InSec. III we show that the GED model is effectively ergodic and that there is economicmobility for λ <
1. In Sec. IV we find that we can define an energy that satisfies theBoltzmann distribution. We characterize the phase transition at λ = 1 in Sec. V andestimate the critical exponents β and γ associated with the order parameter and its variance,respectively. The consequences of the system being describable by a mean-field theory arediscussed in Sec. VI, where we estimate the critical exponent α associated with the energyand the specific heat. In Sec. VII we show that there is critical slowing down. We summarizeand discuss our results in Sec. VIII. In the appendix we argue from economic data that ourmethod of distributing the growth is a reasonable zeroth order approximation. II. STEADY STATE WEALTH DISTRIBUTION
Because the total wealth increases exponentially, it is convenient to introduce the rescaledwealth of an agent, (cid:101) w i ( t ) = NW ( t ) w i ( t ) , (2)and consider the rescaled wealth distribution of the N agents. That is, after the increasedwealth due to growth is distributed to the agents, their wealth is scaled so that the totalrescaled wealth equals N , the initial total wealth. In the following all references to the wealthof the agents will be to their rescaled wealth, and we will omit the tilde for simplicity.To go from time t to time t + 1, we do N wealth exchanges and then apply Eq. (1) tocompute the wealth to be added to each agent due to economic growth. After the addedwealth is distributed to all agents, we rescale the wealth of each agent according to Eq. (2)4 data for figure1 4:17:13 PM 11/18/18 H100H40016006400 l og ( w ) rank FIG. 1. The time dependence of the rescaled wealth distribution as a function of the rank of N = 2500 agents for λ = 0 at t = 100 ( • ), t = 400 ( (cid:78) ), t = 1600 ( (cid:3) ), and t = 6400 ( × ). For t > N . so that the total rescaled wealth is N . The simulations in this section are for N = 2500, f = 0 . µ = 0 . N , f , and µ .We show in Fig. 1 the time dependence of the rescaled wealth distribution for λ = 0,starting from the initial condition w i ( t = 0) = 1 for all i . The wealth disparity betweenricher and poorer agents initially increases until a steady state is established. Once a steadystate is reached, the rescaled wealth distribution remains fixed, and the wealth in every rankincreases as e µt .According to Eq. (1), the growth allocation is weighted more toward the richer agents as λ →
1, thus leading to a less equal steady state rescaled wealth distribution (see Fig. 2).The time to reach a steady state increases as λ approaches 1, and as for λ = 0, the wealth ofall agents increases exponentially. We find that for 0 ≤ λ <
1, “a rising tide lifts all boats”5 abridged l og ( w ) rank FIG. 2. The rescaled wealth distribution versus the rank of N = 2500 agents for λ = 0 . • ), λ = 0 . (cid:78) ), λ = 0 . (cid:4) ), and λ = 0 . (cid:7) ) at t = 10 after a steady state has been reached. Thewealth distribution becomes less equal as λ → − . and all agents benefit from economic growth.The time-dependence of the rescaled wealth distribution is shown for λ = 1 in Fig. 3. Asteady state is not reached, and the slope of the rescaled wealth distribution increases withtime, corresponding to the accumulation of wealth by fewer and fewer agents until eventuallya single agent gains almost all the wealth. Similar results are found for λ >
1. The time fora single agent to dominate decreases as λ increases for λ > λ = 1 is a special case for which the increase in an agent’s wealth due to growthis proportional to the agent’s wealth. Hence, for λ = 1, the ratio of the rescaled wealth ofany two agents does not change after the distribution of the growth in wealth according toEq. (1). Consequently, aside from the exponential growth of the total wealth, the entiredynamics of the system is driven only by the wealth exchange mechanism, and the evolutionof the wealth distribution for λ = 1 is identical to its evolution in the original Yard-Salemodel; that is, the model with no economic growth.6 rankeddataforfig3 t500t10000t50000 l og ( w ) rank t = 500t = 10000t = 50000 FIG. 3. The rescaled wealth distribution at t = 500 ( • ), t = 10000 ( (cid:4) ), and t = 50000 ( (cid:7) ) for λ = 1. In contrast to the behavior for λ <
1, the slope of the rescaled wealth distribution decreaseswith time and a steady state is not reached. The time for a single agent to gain almost all of thewealth is a decreasing function of λ for λ > III. EFFECTIVE ERGODICITY AND ECONOMIC MOBILITY
The numerical results in Sec. II indicate that the GED model exhibits distinct behaviorfor λ < λ ≥
1. That is, for λ <
1, all agents benefit from economic growth, whereasfor λ ≥ λ <
1, but is not ergodic for λ ≥
1. In Sec. III B we find that theagents have nonzero economic mobility for λ <
1, but have zero mobility for λ ≥
1. Unlikein Sec. II, we randomly assign the wealth of each agent at t = 0 from a uniform distributionand then rescale the agents’ wealth so that the initial total wealth equals N .7 . Effective ergodicity To determine whether the system is effectively ergodic, we define the (rescaled) wealthmetric as [26] Ω( t ) = 1 N N (cid:88) i =1 (cid:2) w i ( t ) − (cid:104) w ( t ) (cid:105) (cid:3) , (3)where w i ( t ) = 1 t (cid:90) t w i ( t (cid:48) ) dt (cid:48) (4)and (cid:104) w ( t ) (cid:105) = 1 N N (cid:88) i =1 w i ( t ) . (5)The metric Ω( t ) in Eq. (3) is a measure of how the time averaged rescaled wealth of eachagent approaches the rescaled wealth averaged over all agents. If the system is effectivelyergodic, Ω( t ) ∝ /t [26]. Effective ergodicity is a necessary, but not sufficient condition forergodicity.The linear time-dependence of Ω(0) / Ω( t ) shown in Fig. 4(a) for λ < λ <
1. In contrast, Ω(0) / Ω( t ) for λ = 1 does not increaselinearly with t [see Fig. 4(b)], and hence the system is not ergodic. Similar results are foundfor λ > B. Economic mobility
In a system with economic mobility, poorer agents can become wealthier and richer agentscan become poorer. In contrast, in a system with very low economic mobility, agents rarelychange their rank and the rich become richer.To determine the mobility, we rank the agents according to their wealth at various timesand compute the correlation function C ( t ) of the ranks of the agents once a steady state hasbeen reached for λ <
1. The rank correlation function C i ( t ) of agent i is defined as C i ( t ) = (cid:104) R i ( t ) R i (0) (cid:105) − (cid:104) R i (cid:105) (cid:104) R i (cid:105) − (cid:104) R i (cid:105) , (6)where R i ( t ) is the rank of agent i at time t and (cid:104) R i (cid:105) = N/
2. The corresponding quantityfor λ ≥
1, where a steady state is not reached, is the Pearson correlation function given by8
Data 69 lam=0.9lam=0.5 Ω ( ) / Ω ( t ) t Data 68 lam = 1 Ω ( ) / Ω ( t ) t FIG. 4. (a) The linear time-dependence of the inverse rescaled wealth metric indicates that thesystem is effectively ergodic for λ = 0 . λ = 0 .
9. The corresponding inverse slopesare 3777 and 44300 for λ = 0 . λ = 0 .
9, respectively. (b) The inverse wealth metric for λ = 1 . N = 5000, f = 0 .
01, and µ = 0 . coefficient [27] C i ( t ) = (cid:2) R i ( t ) − (cid:104) R ( t ) (cid:105) (cid:3)(cid:2) R i (0) − (cid:104) R (0) (cid:105) (cid:3)(cid:113)(cid:2)(cid:0) R i ( t ) − (cid:104) R i ( t ) (cid:105) (cid:1) (cid:3)(cid:2)(cid:0) R i (0) − (cid:104) R i (0) (cid:105) (cid:1) (cid:3) , (7)The correlation function averaged over all agents is C ( t ) = (1 /N ) (cid:80) i C i ( t ). As can be seenfrom Fig. 5(a), C ( t ) → t → ∞ for λ <
1, which indicates that the rank of an agent as t → ∞ is not correlated with its rank at t = 0. Hence, the agents have a nonzero economicmobility for λ <
1. In contrast, in Fig. 5(b) we see that C ( t ) approaches a constant for λ ≥
1, indicating that the ranks are strongly correlated at different times, and there is noeconomic mobility.
IV. EQUILIBRIUM, NOT JUST STEADY STATE
We have seen that the GED model approaches a steady state and is effectively ergodic for λ <
1. In the following, we will show that a reasonable definition of the total energy yieldsan energy distribution that is consistent with the Boltzmann distribution. The existence ofthe latter is consistent with the idea that that the system is not just in a steady state, butis in thermodynamic equilibrium for λ <
1. 9
Data 46 log C(t)y = -0.00066049 - 2.6543e-6x R= 0.99945 l og C ( t ) t Data 49
C(t)C(t)1.1 C ( t ) t FIG. 5. (a) The rank correlation decays exponentially (red curve) for λ = 0 .
9, indicating that themobility is nonzero. (b) The Pearson correlation function for λ = 1 . • ) and λ = 1 . (cid:3) ). For λ ≥ C ( t ) remains nonzero even as t → ∞ , which indicates that the rank of an agent remainscorrelated and there is no mobility ( N = 2500, f = 0 .
1, and µ = 0 . As discussed in Ref. [28] (following paper), a mean-field theory treatment of the GEDmodel yields a quantity that can be interpreted as the total energy of the system E = N (cid:88) i =1 (1 − w i ) . (8)Note that the energy is zero if all agents have the same wealth ( w i = 1). Equation (8)yields a mean energy that is extensive, that is, (cid:104) E (cid:105) ∝ N for fixed values of λ , f , and µ .For example, for λ = 0 . f = 0 .
01, and µ = 0 .
1, we find that (cid:104) E N =4000 (cid:105) / (cid:104) E N =1000 (cid:105) =336 . / . . ≈ P ( E ) is shown in Fig. 6(a) for N = 5000, λ = 0 . f = 0 .
01, and µ = 0 .
1. As expected, P ( E ) is fit well by a Gaussian. Fits of P ( E ) to a Gaussian becomeless robust as λ → N , f , and µ . We will discuss the behavior of P ( E ) for largervalues of λ in Sec. VI.If the system is in thermodynamic equilibrium, we expect that the probability density P ( E, β ) to be proportional to g ( E ) e − βE , where β is an effective inverse temperature thatdepends on the parameters λ , f , µ , and g ( E ) is the density of states, which is independentof λ , f , and µ and hence independent of β . Because g ( E ) is independent of β , the ratio P ( E, β ) /P ( E, β ) is an exponential proportional to exp[ − ( β − β )] E if the system is char-10 -3 -3 -3 -3
22 23 24 25 26 27 28
P(E)2005637019
P(E) P ( E ) E y = m1*exp(-m3*(x-m2)^ 2)ErrorValue 3.0684e-60.0085316m1 0.0001942424.764m2 0.0018992.2857m3 NA2.3875e-7Chisq NA0.99997R -2 -1
22 23 24 25 26 27 ratio of probabilities ratioy = 3.5357e+11 * e^(-1.0769x) R= 0.85056 P ( E , f = . ) / P ( E , f = . ) E FIG. 6. (a) The energy probability density P ( E ) for N = 5000, λ = 0 . f = 0 .
01, and µ = 0 . P ( E ) ∝ exp( E − E ) /σ E ) with E = 24 . σ E = 0 .
44[ (red) curve]. (b) The ratio P ( E, f = 0 . /P ( E, f = 0 . e − ∆ βE with ∆ β ≈ .
08 (red curve). acterized by the Boltzmann distribution with the energy given by Eq. (8). The range ofvalues of E over which this ratio is nonzero and finite is limited by the overlap of the twoprobabilities, which becomes smaller as N is increased. In Fig. 6(b) we see that the ratio P ( E, β ) /P ( E, β ) is consistent with the Boltzmann distribution e − ( β − β ) E with β ∝ f − , f = 0 . f = 0 .
01, with a larger value of f (for fixed values of λ and µ ) corre-sponding to a smaller value of β and hence a higher value of the effective temperature. Thedependence of an effective temperature β − on f is consistent with the association of f withthe presence of noise in the system [28]. Similar results hold for two similar values of λ for fixed f and µ . The existence of an energy, and the observation that its probability isproportional to the Boltzmann distribution implies that the system can be considered to bein thermodynamic equilibrium for λ = 0 .
8. 11 . CHARACTERIZATION OF THE PHASE TRANSITIONA. Fixed number of agents
The numerical results discussed in this section are for N = 5000, f = 0 .
01, and µ = 0 . ( N × exchanges) after a transient time of 10 . Themajor source of uncertainty in our estimations of the values of the various critical exponentsis the choice of the range of values to be retained in the least squares fits.Because we will characterize the approach to the phase transition at λ = 1 in terms ofpower laws, it is convenient to define the quantity (cid:15) ≡ − λ. (9)To characterize the phase transition at (cid:15) = 0, we need to identify an order parameter.A common measure of income or wealth inequality is the Gini coefficient G n [29]. If allagents have the same wealth, G n = 0. In contrast, if one agent has all the wealth, G n = 1,corresponding to the maximum degree of inequality. These characteristics would appear tomake 1 − G n a good choice for the order parameter. However, for reasons that are discussedin Ref. [28], the fluctuations of the Gini coefficient are zero in the limit N → ∞ , and hencethe susceptibility, which would be associated with the variance of G n , would be zero, making G n not an appropriate choice of the order parameter. The order parameter and the value of β . Another possible choice of the order parameteris the fraction of the wealth held by all the agents except the richest agent, that is, φ = N − w max N , (10)where w max is the wealth of the richest agent. For λ < w max (cid:28) N and dependsonly weakly on N for given values of λ , f , and µ . For example, w max = 2 . N = 1000and w max = 2 . N = 4000, λ = 0 . f = 0 .
01, and µ = 0 .
1. The weak dependence of w max on N and also on λ implies that φ defined in Eq. (10) approaches one as N → ∞ for λ < (cid:15) . Hence, the value of the critical exponent β associatedwith the order parameter is β = 0. The susceptibility and the value of γ . The (cid:15) -dependence of the variance of w max is shown inFig. 7(a) and can be fit to a power law to give an effective exponent for the susceptibility closeto one for (cid:15) ≥ × − ; fits for smaller values of (cid:15) give an effective exponent approximately12 -4 -3 -2 -1 -4 -3 -2 -1 varwmaxN5000data var wmaxfar fitnnear fit v a r i a n ce o f t h e w ea lt h o f t h e r i c h e s t a g e n t λ slope = -1slope = -1.7 -1 -4 -3 -2 C2N5000 C2 v a r i a n ce o f t h e w ea lt h d i s t r i bu ti on λ y = m1*M0^(-1) ErrorValue 1.5631e-60.00063521m1 NA0.14011Chisq NA0.99964R FIG. 7. (a) The (cid:15) -dependence of w max , the variance of the wealth of the richest agent, showscrossover behavior. Fits of the variance to a power law for (cid:15) ≥ × − give an effective exponentclose to one (red curve); fits for smaller values of (cid:15) give an effective exponent of approximately1.7 (blue curve). (b) The (cid:15) -dependence of C , the variance of the wealth of all the agents, showsless curvature than the (cid:15) -dependence of the variance of w max , and gives an effective exponent closeto 1.0 (red curve) if fits are made for (cid:15) < . . , . (cid:15) that are included in the least squares fits. equal to 1.7. Both estimates indicate that the variance of the wealth of the richest agentdiverges strongly, and hence we choose the order parameter to be as given in Eq. (10).More consistent results for the susceptibility can be found from C , the variance of thedistribution of wealth of all the agents, and not just the richest agent [see Fig. 7(b)]. Wesee that there is less curvature in the plot of log C versus log (cid:15) , and we find an effectiveexponent close to one if fits are made for (cid:15) < . C are included for largervalues of (cid:15) , the effective exponent from the least squares fits is in the [0 . , . N C : χ = N C , (11)and conclude that the critical exponent associated with the susceptibility is γ ≈ The mean energy, heat capacity, and the value of α . Although the total rescaled wealthis a constant, the energy depends on the way the wealth is distributed. We use Eq. (8) todetermine the mean energy by averaging the quantity (cid:80) Ni =1 (1 − w i ) over many realizations .13 -3 -2 -1 -5 -4 -3 -2 -1 E t o t a l e n e r gy λ y = m1*x^(-1) ErrorValue 1.0199e-60.00064679m1 NA0.061092Chisq NA0.99987R -1 -4 -3 -2 -1 varE var sum v a r i a n ce o f t h e e n e r gy λ y = m1*M0^(-2) ErrorValue 7.8833e-50.0037189m1 NA9.7497e+10Chisq NA0.96279R FIG. 8. (a) The divergent (cid:15) -dependence of the mean total energy (cid:104) E (cid:105) for fixed N is consistent withthe power law (cid:15) − (red curve). Least squares fit of (cid:104) E (cid:105) yield values of the divergence in the range[0 . , . (cid:15) -dependence of the heat capacity C is consistent with the power law (cid:15) − (redcurve). Least squares fits of C yield values of the effective exponent in the range [1 . , . Our results for the mean energy (cid:104) E (cid:105) are shown in Fig. 8(a). We see that the (cid:15) -dependence of (cid:104) E (cid:105) is consistent with (cid:15) − α as (cid:15) → α ≈
2. This divergence of (cid:104) E (cid:105) is inconsistent withequilibrium statistical mechanics, which requires that the total energy be finite for finitevalues of N .We define the heat capacity C to be the variance of the total energy defined in Eq. (8).The (cid:15) -dependence of of C is consistent with (cid:15) − α with the exponent α = 2 as shown inFig. 8(b). Our determination of the value of α depends on the range of values of (cid:15) that areincluded in the least squares fits and are in the range [1 . , . α , β , and γ , determined as (cid:15) is varied for a given value of N , are consistent with β = 0 , γ = 1 , and α = 2 . (12)These numerical values are inconsistent with the usual scaling law α + 2 β + γ = 2 . (fixed N ) (13)However, the value of α determined from the power law behavior of the heat capacity isconsistent with the λ -dependence of the mean total energy, that is, C = ∂ (cid:104) E (cid:105) /∂λ .14 I. MEAN-FIELD THEORY AND CONSTANT GINZBURG PARAMETER
Although it is natural to determine the critical behavior of the GED model as the criticalpoint is approached for a fixed number of agents, our numerical results for the criticalexponents are not consistent with the scaling law, Eq. (13), nor consistent with equilibriumstatistical mechanics because the mean energy (cid:104) E (cid:105) diverges as the critical point is approachedeven for finite N . Simulations for fixed λ , f , and µ also show that (cid:104) E (cid:105) and the heat capacity C are proportional to N so that the divergent behavior of (cid:104) E (cid:105) is not removed by first takingthe limit N → ∞ before taking the limit (cid:15) → β = 0 , γ = 1 , and α = 1 . (14)The predicted mean-field values of the exponents in Eq. (14) are consistent with Eq. (13).The mean-field theory [28] also predicts that the mean energy approaches a constant as (cid:15) → N weneed to account for the fact f and µ are rates and that they depend on the definition of theunit of time (see Sec. VII). It is shown in the mean-field theory of Ref. [28] that we need torescale f and µ as, f = f /N and µ = µ /N, (15)to achieve a consistent thermodynamic description of the GED model.Another condition for the applicability of the mean-field treatment of the GED model isthat the Ginzburg parameter G , defined as G ≡ µ f N (cid:15), (16)be much greater than one and be held constant as (cid:15) → .60 x 10 -3 -3 -3 -3 -3 -3 -3 -3 sumw2Gdata sum w^2/Ny = 0.0046516 + 0.0017344x R= 0.99566 E / N λ FIG. 9. The (cid:15) -dependence of (cid:104) E (cid:105) /N , the mean energy per agent, for G = 10 is consistent withthe linear (cid:15) -dependence a + a (cid:15) as λ →
1, with a ≈ .
005 and a ≈ .
002 (red curve). fixed as the critical point is approached. In contrast, we find in the following that the resultsfor β and γ do not depend on keeping G fixed as has been found for other fully connectedmodels [30–32].We emphasize that if µ and f are not rescaled in the simulations, the energy per agentwould diverge as N → ∞ even for fixed Ginzburg parameter.To compare the mean-field theory predictions to the simulations, we choose f = 0 . µ = 0 .
1, and N = 5000, with f = ( N /N ) f and µ = ( N /N ) µ . For a particular choice ofthe value of λ , we determine the value of N needed to keep the value of G in Eq. (16) fixedat G = 10 . Our simulations are for 0 . ≤ λ ≤ .
998 and 5 × ≤ N ≤ × .Because the results for (cid:104) E (cid:105) and C depend on keeping G fixed, we first discuss their (cid:15) dependence. In Fig. 9 we see that (cid:104) E (cid:105) /N approaches a constant as (cid:15) →
0, in contrast to itsdivergent (cid:15) behavior for fixed N . Because (cid:104) E (cid:105) ∼ (cid:15) − α , we see that simulations for constantGinzburg parameter are consistent with α = 1.Our numerical results for (cid:104) E (cid:105) for fixed N are consistent with the relation (cid:104) E (cid:105) ∼ Ng , (17)where g ≡ µ(cid:15)f . (18)16quations (17) and (18) imply that (cid:104) E (cid:105) exhibits different behaviors for fixed N and for fixed G . For fixed N we find (cid:104) E (cid:105) ∼ Ng = N f µ(cid:15) ∝ N(cid:15) (fixed N ) (19)Equation (19) implies that for fixed values of λ , f , and µ , (cid:104) E (cid:105) is proportional to N forfixed λ and diverges as (cid:15) − for a given value of N ; both behaviors are consistent with oursimulations.For fixed Ginzburg parameter we use Eqs. (15)–(18) to write (cid:104) E (cid:105) as (cid:104) E (cid:105) = N f N µ (cid:15) = f µ (cid:15) = NG (fixed Ginzburg parameter) . (20)Equation (20) is predicted by mean-field theory [28]. If simulations are done at constantGinzburg parameter, then (cid:104) E (cid:105) /N is predicted to approach a constant, consistent with theresults of our simulations.The (cid:15) -dependence of C , the variance of the total energy, for constant G is shown inFig. 10. We find that C ∼ (cid:15) − α , with α ≈
1, consistent with the prediction of mean-fieldtheory [28].Our numerical results for the variance of the total energy are consistent with the relation C ∼ Ng . (21)Equation (21) implies that for fixed NC ∼ N(cid:15) (fixed N ) . (22)In this case α = 2 and C is proportional to N for fixed values of λ , f , and µ . Because C isproportional to N , it is tempting to interpret it as a measure of the heat capacity. As wewill see the interpretation is more subtle.For fixed Ginzburg parameter, the interpretation of C is different. From Eq. (21) we have C ∼ N f N µ (cid:15) = NG ∼ G (cid:15) (fixed Ginzburg parameter) , (23)where N ∼ (cid:15) − for fixed G . The dependence of C on G and N is predicted in Ref. [28].Note that C in Eq. (23) is independent of N , and hence might be interpreted as a specificheat. However, the relation of C for fixed N to that for fixed G is not the usual one because17 -1 -3 -2 -1 specificheat var sum v a r i a n ce o f e n e r gy λ y = m1*M0^(-1) ErrorValue 0.0011320.033248m1 NA168.54Chisq NA0.91058R FIG. 10. The (cid:15) -dependence of C , the variance of the total energy, for constant Ginzburg parameteris consistent with the power law (cid:15) − α , with α = 1 (red line). A least square fit gives an exponent of ≈ . the number of agents N in the Ginzburg parameter is not independent of (cid:15) as is usually thecase.The (cid:15) − dependence of the variance of the total energy near (cid:15) = 0 implies that themean energy must include a logarithmic dependence on λ . For example, the form, (cid:104) E (cid:105) ∼ a + a (cid:15) + a L / ln (cid:15) , implies that C = ∂ (cid:104) E (cid:105) /∂λ scales as (cid:15) − with a logarithmic correction.Mean-field theory is incapable of finding logarithmic corrections, and our data for (cid:104) E (cid:105) is notsufficiently accurate to detect the presence of logarithmic factors.Although our numerical results for α depend on whether N or G is held fixed as λ → β and γ are independentof the nature of the approach to the phase transition. We find that w max (cid:28) N and isindependent of λ for fixed G and hence φ = 1 for λ <
1, implying that β = 0 as was foundfor fixed N . The divergence of the susceptibility χ = N C shown in Fig. 11(b) is consistentwith χ ∼ (cid:15) − γ with γ = 1. 18 -3 -2 -1 Data 177
NC2 N C λ y = m1*(M0)^(-1) ErrorValue 0.00128644.6513m1 NA437.36Chisq NA0.99999R FIG. 11. The divergence of
N C is consistent with χ ∼ (cid:15) − γ with γ = 1 (red curve). VII. CRITICAL SLOWING DOWN
We find that various time scales increase rapidly as (cid:15) →
0, which limits how close thesimulations can approach the transition. One time scale of interest is τ ra , the mean lifetimeof the richest agent. We expect that if the mobility of the agents is nonzero, then the richestagent at t = 0 will no longer be the richest after some time has elapsed. We define τ ra asthe mean time that a particular agent remains the richest and assume that τ ra is a simplemeasure of the decorrelation time of the wealth of individual agents.Another time scale of interest is the mixing time τ m associated with the time-dependenceof the wealth metric; τ m is related to the inverse slope of the wealth metric asΩ(0) / Ω( t ) = t/τ m . (24)We also computed the time-displaced energy autocorrelation function given by C E ( t ) = (cid:104) E ( t ) E (0) (cid:105) − (cid:104) E (cid:105) (cid:104) E (cid:105) − (cid:104) E (cid:105) , (25)where E ( t ) is the value of the energy of the system at time t . We find that C E ( t ) relaxes19 lifetimeG lifetimey = 0.60694 * x^(-0.90501) R= 0.99068 li f e ti m e o f t h e r i c h e s t a g e n t − λ bottom curve is linear fit energydecorrelationtime energyenergymixing m i x i ng a nd e n e r gy d ec o rr e l a ti on ti m e s λ mixing timedecorrelation time FIG. 12. (a) The (cid:15) -dependence for constant Ginzburg parameter ( G = 10 ) of τ ra , the lifetime ofthe richest agent, diverges as τ ra ∼ (cid:15) − . (red curve). (b) The (cid:15) -dependence of the mixing time τ m (top red curve) and τ E (bottom blue curve) diverge as ∼ (cid:15) − . Least squares fits give an exponentin the range [1 . , . exponentially, and hence we can extract the energy decorrelation time τ E . Our simulationresults for these various times are for constant Ginzburg parameter ( G = 10 ).The simulation results for the (cid:15) -dependence of τ ra are shown in Fig. 12(a) and are con-sistent with the power law (cid:15) − . To obtain accurate results for the wealth metric Ω( t ), weaveraged Ω( t ) over ten origins and found that the linear dependence of Ω(0) / Ω( t ) holds overa wide range of t , thus yielding robust values of τ m . The exponential dependence of C E ( t )holds for t (cid:46) τ E , yielding some uncertainty in the fitted values of τ E . The (cid:15) -dependence of τ m and τ E are shown in Fig. 12(b). We see that both τ m and τ E increase rapidly as (cid:15) → (cid:15) -dependence is consistent with (cid:15) − . Although the mean-field theory [28]makes no direct predictions for the mixing time, the (cid:15) -dependence of τ m can be understoodby noting that the metric measures the time it takes for the average wealth of each agentto equal the global average. Because the time it takes for the richest agent to cease beingthe richest and for another to assume that role diverges as approximately (cid:15) − , the time fora system of N agents to mix is N (cid:15) − . Because N ∝ (cid:15) − for fixed Ginzburg parameter [seeEq. (16)], we find that τ m ∼ (cid:15) − in agreement with the simulations.The mean-field theory of Ref. [28] predicts that the decorrelation time for fixed Ginzburg20arameter scales as τ mf ∼ µ (cid:15) . (26)The reason for the apparent discrepancy between the (cid:15) -dependence predicted by Eq. (26) andour simulation results for τ m and τ E is that the time unit in the simulations corresponds to N exchanges, during which one agent exchanges wealth with only one agent on the average.In contrast, the applicability of mean-field theory requires that in one time unit one agentexchanges wealth with N agents on the average. Hence, the simulation and mean-field timeunits differ by a factor of N or 1 /(cid:15) for constant Ginzburg parameter.The divergent behavior of τ m and τ E are examples of critical slowing down, which isassociated with a cooperative effect and is not a property of a single agent. In contrast, τ ra is a property of a single agent rather than of the system as a whole and becomes independentof (cid:15) if we define the time as required by the applicability of mean-field theory.Although the results of our simulations are consistent with the (cid:15) -dependence in Eq. (26)for constant Ginzburg parameter, Eq. (26) also predicts that τ E is independent of the valuesof f and N . As discussed in Ref. [28], the derivation of Eq. (26) neglects the effects of bothadditive and multiplicative noise. The weak dependence of τ m and τ E on N and f , as wellas their dependence on µ , is discussed in Ref. [28]. VIII. DISCUSSION
We have generalized the Yard-Sale model to incorporate economic growth and its distri-bution according to the wealth of the agents as determined by Eq. (1) and the parameter λ . Our numerical results suggest that there are two phases. For λ < λ ≥ t → ∞ , there is condensation ofa finite fraction of the system’s wealth in a vanishingly small number of agents. In addition,the system is not ergodic and shares some of the characteristics of the geometric randomwalk which is also not ergodic and cannot be treated by equilibrium methods [4, 34].It is remarkable that it is possible to define a thermodynamic energy for a system thatinvolves wealth and has no obvious energy analogue. The interpretation of the energy andits variance is subtle and thermodynamic consistency is achieved only if the mean-field limit21s taken appropriately.We showed in Sec. IV that P ( E ), the probability density of the energy of the system, iswell fit by a Gaussian function for N = 5000 and λ = 0 .
8. Simulations for N = 5000 andvalues of λ much closer to one show departures from a Gaussian, even though the wealthfluctuation metric still indicates that the system is effectively ergodic. The deviation of P ( E ) from a Gaussian for fixed N (and fixed µ and f ) is due to the fact that G decreases as λ → λ much closer to one. However, although the Ginzburg parameter is large and fixed, theimportance of the multiplicative noise increases as λ →
1, and eventually the effect of themultiplicative noise can no longer be ignored [28].Although our numerical values of the various exponents are consistent with the mean-fieldtheory of Ref. citekleinmf, their estimated numerical values must be viewed with cautionbecause they are obtained by extrapolation over a limited range of λ and for finite valuesof G and N . Much larger values of G would be needed to obtain more accurate numericalresults for λ closer to one.The transition at λ = 1 is from a system in thermodynamic equilibrium for λ < λ ≥
1. In Ref. [28] it is shown that theevolution of the model for λ ≥ λ has been changed from λ = 0 . λ = 1 .
1. We see thatthe wealth of the richest agent initially increases exponentially as predicted by the theory.Because every agent can trade with every other agent, the GED model and similar agent-based models can be considred to be fully connected. As we have seen, fully connectedmodels, such as the fully connected Ising model [31, 32], can be treated consistently bymean-field theory near the critical point only if the Ginzburg parameter is much greaterthan one and is held constant as the critical point is approached. If the Ginzburg parameteris not held constant, anomalous results for the energy and the specific heat are found.Our results have possible important economic implications. As λ is increased, the benefitsof growth are weighted more toward the wealthy, and wealth inequality increases. Neverthe-less, as long as λ <
1, the wealth of all ranks grows at the same rate once a steady state isreached, and all agents benefit from economic growth. However, if the benefits of growth are22 superquenchN20000 richestfit w ea lt h o f r i c h e s t a g e n t t FIG. 13. The wealth of the richest agent after a instantaneous change of λ from λ = 0 . λ = 1 . e − t/τ q (red curve) with τ q ≈
137 for t (cid:46) N = 20000, f = 0 .
01, and µ = 0 . skewed too much toward the wealthy ( λ ≥ λ = 1 the model reduces to thegeometric random walk with resulting wealth condensation [4], as we show in the followingpaper [28].There is some question whether economic systems can be treated as being in equilibriumor even exhibit effective ergodicity [4, 10, 34]. Our results suggest that ergodicity andequilibrium may depend on various system parameters. Because parameters such as λ and µ are not temporal constants in real economies, our results also suggest that the applicabilityof equilibrium methods may be situational and vary with time. ACKNOWLEDGMENTS
We thank Ole Peters, John Ogren, Alan Gabel, Timothy Khouw, and Louis Colonna-Romano for useful conversations. 23 ppendix: Comparison to some economic data
Modeling the economy of a country as large and diverse as the United States by com-pressing economic growth and transactions into three parameters, λ , f , and µ is a grosssimplification. The assumption that these parameters are independent of time also is unre-alistic. In the following, we analyze the growth data [35] and wealth distribution data [36]compiled by Karen Smith and published by the Urban Institute. Our analysis suggests thatthe assumptions that the distribution of growth can be modeled as in Eq. (1) and that theparameters are independent of time is a reasonable zeroth order approximation to the distri-bution of wealth in the real economy. We will discuss the exchange term and its relevanceto the real economy in the following paper [28].From the growth rate of the gross domestic product shown in constant dollars in Ref. [35],we note that the temporal fluctuations of the (inflation adjusted) growth rate of the grossdomestic product exhibits large swings that appear to be damped as a function of time. Thedecline in the growth caused by the great recession starting in 2008 is an example of a largefluctuation, but the growth rate has remained close to the mean rate of roughly 3% over thelast 30 years, thus implying that µ ≈ . λ as λ = log (cid:32) W r ( t ) − W r ( t ) W r ( t ) (cid:33) , (A.1)where W r ( t ) is the wealth of people of economic rank (percentile) r at time t .We estimated λ for the 50th, 90th and 95th wealth percentiles in the intervals 1983–1989, 1995 –1998 and 2013–2016 (see Table I). Although the values of λ are not constantfor different time intervals and percentiles, they vary by only a few percentage points as afunction of percentile. They vary more as a function of time, with the 50% percentile havingthe greatest variation. The change of λ appears to decrease for later times, consistent withthe damping of the variation of the growth. The variation of λ is more pronounced for evenlower rankings. However, because the wealth of the lower rankings is considerably smaller,the variation of the value of λ has less effect on the wealth of the poor.We conclude from the growth and wealth distribution data [35, 36] that the distributionof economic growth assumed in the GED model is a reasonable zeroth order approximation,24 ercentile 1983–1989 1995–1998 2013–201695% 0.85 0.90 0.9090% 0.84 0.89 0.9050% 0.75 0.90 0.84TABLE I. The calculated values of λ for the percentiles and time intervals indicated using economicdata from Refs. [35] and [36]. particularly for the upper half of the wealth ranks of the United States over the past 30years. Of course, there is much that it is not included in the model, such as the effects ofwars, famines, storms, and recessions, which are not obtainable from the simple GED model.However, it appears from the data that the model is a reasonable approximation over timescales of the order of decades and yields insights into the importance of how economic growthis distributed. [1] A. Chakraborti, I. M. Toke, M. Patriarca, and F. Abergel, “Econophysics review: II. Agent-based models,” Quant. Finance , 1013 (2011).[2] John F. Kennedy, “Remarks in Heber Springs Arkansas at the dedication of the Greers FerryDam,” The American Presidency Project, 3 October 1963, .[3] E. Gudreis, “Unequal America,” Harvard Magazine, July–August (2008).[4] O. Peters, “Optimal leverage from non-ergodicity,” Quant. Finance , 1593 (2011).[5] J. Angle, “The surplus theory of social stratification and the size distribution of personalwealth,” Social Forces , 293 (1986).[6] J. Angle, “Deriving the size distribution of personal wealth from the rich get richer, the poorget poorer,” J. Math. Sociology , 27 (1993).[7] A. Chakraborti and B. K. Chakrabarti, “Statistical mechanics of money: how saving propen-sity affects its distribution,” Eur. Phys. J. B , 167 (2000).[8] A. Dr¨agulescu and V. M. Yakovenko, “Statistical mechanics of money,” Eur. Phys. J. B ,723 (2000).[9] C. F. Moukarzel, S. Goncalves, J. R. Iglesias, M. Rodriguez-Achach, and R. Huerta- uintanilla, “Wealth condensation in a multiplicative random asset exchange model,” Eur.Phys. J.-Spec. Top. , 75 (2007).[10] V. M. Yakovenko and J. B. Rosser Jr., “Colloquium: Statistical mechanics of money, wealth,and income,” Rev. Mod. Phys. , 1703 (2009).[11] M. Patriarca and A. Chakraborti, “Kinetic exchange models: From molecular physics to socialscience,” Am. J. Phys. , 618 (2013).[12] J.-P. Bouchaud and M. M´ezard, “Wealth condensation in a simple model of economy,” PhysicaA , 536 (2000).[13] Z. Burda, D. Johnston, J. Jurkiewicz, M. Kami´nski, M. A. Nowak, G. Papp, and I. Zahed,“Wealth condensation in Pareto macroeconomies,” Phys. Rev. E , 026102 (2002).[14] B. M. Boghosian, “Kinetics of wealth and the Pareto law,” Phys. Rev. E , 042804 (2014).[15] B. M. Boghosian, “Fokker-Planck description of wealth dynamics and the origin of Pareto’sLaw,” Int. J. Mod. Phys. C , 1441008 (2014).[16] B. M. Boghosian, M. Johnson, and J. A. Marcq, “An H theorem for Boltzmann’s equation forthe yard-sale model of asset exchange,” J. Stat. Phys. , 1339 (2015).[17] Christophe Chorro, “A simple probabilistic approach of the Yard-Sale model,” Statistics andProbability Letters , 35 (2016).[18] B. M. Boghosian, A. Devitt-Lee, M. Johnson, J. Li, J. A. Marcq, and H. Wang, “Oligarchyas a phase transition: The effect of wealth-attained advantage in a Fokker-Planck descriptionof asset exchange,” Physica A , 15 (2017). “Wealth-attained advantage” is implementedby biasing the coin flip for the exchange of wealth in favor of the wealthier agent. In contrast,we implement a form of wealth-attained advantage by an uneven redistribution of wealth dueto growth.[19] A Devitt-Lee, H. Wang, J. Ii, and B. Boghosian, “A nonstandard description of wealth con-centration in large-scale economies,” Siam J. Appl. Math , 996 (2018).[20] J. Li, B. M. Boghosian, and C. Li, “The affine wealth model: An agent-based model of assetexchange that allows for negative-wealth agents and its empirical validation,” Physica A ,423 (2019).[21] B. M. Boghosian, “The inescapable casino,” Sci. Amer. (11), 72 (2019).[22] P. L. Krapivsky and S. Redner, “Wealth distributions in asset exchange models,” Science andCulture , 424 (2010).
23] S. Ispolatov, P. L. Krapivsky, and S. Redner, “Wealth distributions in asset exchange models,”J. Eur. Phys. B , 267 (1998).[24] A. Chakraborti, “Distributions of money in model markets of economy,” Int J. Mod. Phys. C , 1315 (2002).[25] B. Hayes, “Follow the money,” Am. Sci. , 400 (2002).[26] D. Thirumalai and R. D. Mountain, “Activated dynamics, loss of ergodicity, and transport insupercooled liquids,” Phys. Rev. A , 4574 (1990) and Phys. Rev. E , 479 (1993).[27] J. L. Rodgers and W. A. Nicewander, “Thirteen ways to look at the correlation coefficient,”Amer. Stat. , 59 (1988).[28] W. Klein, N. Lubbers, K. K. L. Liu, and H. Gould, “Mean-field theory of an asset exchangemodel with economic growth and wealth distribution,” following manuscript.[29] For a discussion of the Gini coefficient, see, for example,