# Power-Law Distributions and Lévy-Stable Intermittent Fluctuations in Stochastic Systems of Many Autocatalytic Elements

## I Introduction

The origins of power-law distributions as well as their conceptual implications have been an active topic of research in recent years. Power laws are intrinsically related to the emergence of macroscopic features which are scale invariant within some bounds and distinct from the microscopic elementary degrees of freedom. Often, these features are insensitive to the details of the microscopic structures. Well known examples of power law distributions include the energy distribution between scales in turbulence , the distribution of earthquake magnitudes, the diameter distribution of craters and asteroids, the distribution of city populations , the distributions of income and of wealth , the size-distribution of business firms  and the distribution of the frequency of appearance of words in texts. The fact that multiplicative dynamics tends to generate power-law distributions was intuitively invoked long ago  but the limitations in computer simulation power kept the models under the constraints imposed by the applicability of analytical treatment. More recently, a broader class of models has been studied combining computer simulations with theoretical analysis within the Microscopic Representation paradigm proposed in Ref. . In particular, it was shown  that power laws appear in a variety of dynamical processes and are maintained even under highly non-stationary conditions.

In this paper we consider a generic model of stochastic dynamics with many degrees of freedom $`w_i(t)`$, $`i=1,\mathrm{},N`$. The time evolution of the $`w_i`$’s is described by an asynchronous update mechanism in which at each time step one variable is chosen randomly and is multiplied by a factor $`\lambda `$ taken from a predefined distribution. In addition, there is a global coupling constraint which does not allow the $`w_i`$’s to fall below the lower cutoff given by $`c\overline{w}`$, where $`\overline{w}`$ is the momentary average of the $`w_i`$’s and $`0<c<1`$ is a constant. The dynamic variables $`w_i`$ are found to exhibit a power-law distribution of the form $`p(w)w^{1\alpha }`$. The exponent $`\alpha `$ is found to be insensitive to the distribution $`\mathrm{\Pi }(\lambda )`$ of the random factor $`\lambda `$. However, $`\alpha `$ is non-universal, and increases monotonically as a function of $`c`$. In the limit $`c=0`$ (where the $`w_i`$’s become decoupled) $`\alpha =0`$ for any finite $`N`$. However, in the ”thermodynamic” limit $`N=\mathrm{}`$, $`\alpha 1`$ for any positive $`c`$. Thus the two limits do not commute. This is important for applications since typically in empirical systems $`\alpha 1`$  unlike the case of the free multiplicative random walk which predicts a log-normal distribution  corresponding to $`\alpha =0`$ .

The time evolution of $`\overline{w}(t)`$ exhibits intermittent fluctuations parametrized by a truncated Lévy-stable distribution with the same index $`\alpha `$. This intricate relation between the distribution of the $`w_i`$’s at a given time and the temporal fluctuations of their average is examined and its relevance to empirical systems is discussed. Our model indicates that in certain cases the scaling exponent may be insensitive to the distribution of the multiplicative (random) factor $`\lambda `$ and depends only on the ”lower bound” features which control the smallest values of the elementary variables. The relation between the limiting conditions and the power law exponent is to be applied in each particular case and it constitutes a strong instrument in identifying and validating the relevant degrees of freedom responsible for the emergence of scaling.

The present paper proposes to consolidate by numerical simulations the control one has on a specific model and help in this way its further application to additional systems. The paper is organized as follows. In Sec. II we present the model. Simulations and results are reported in Sec. III, followed by a discussion in Sec. IV and a summary in Sec. V.

## II The Model

### II.1 Formal Definition

The model describes the evolution in discrete time of $`N`$ dynamic variables $`w_i(t)`$, $`i=1,\mathrm{},N`$. At each time step $`t`$, an integer $`i`$ is chosen randomly in the range $`1iN`$, which is the index of the dynamic variable $`w_i`$ to be updated at that time step. A random multiplicative factor $`\lambda (t)`$ is then drawn from a given distribution $`\mathrm{\Pi }(\lambda )`$, which is independent of $`i`$ and $`t`$ and satisfies $`_\lambda \mathrm{\Pi }(\lambda )𝑑\lambda =1`$. This can be, for example, a uniform distribution in the range $`\lambda _{min}\lambda \lambda _{max}`$, where $`\lambda _{min}`$ and $`\lambda _{max}`$ are predefined limits. The system is then updated according to the following stochastic time evolution equation

$`w_i(t+1)`$ $`=`$ $`\lambda (t)w_i(t)`$   
$`w_j(t+1)`$ $`=`$ $`w_j(t),j=1,\mathrm{},N;ji.`$ (1)   
This is an asynchronous update mechanism. The average value of the system components at time t is given by

$$\overline{w}(t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}w_i(t).$$
(2)   
The term on the right hand side of Eq. (1) describes the effect of auto-catalysis at the individual level. In addition to the update rule of Eq. (1), the value of the updated variable $`w_i(t+1)`$ is constrained to be larger or equal to some lower bound which is proportional to the momentary average value of the $`w_i`$’s according to

$$w_i(t+1)c\overline{w}(t)$$
(3)   
where $`0c<1`$ is a constant factor. This constraint is imposed immediately after step (1) by setting

$$w_i(t+1)\mathrm{max}\{w_i(t+1),c\overline{w}(t)\},$$
(4)   
where $`\overline{w}(t)`$, evaluated just before the application of Eq. (1), is used. This constraint describes the effect of auto-catalysis at the community level.

### II.2 Main Features

Our model is characterized by a fixed (conserved) number of dynamic variables $`N`$, while the sum of their values is not conserved. The conservation of the number of dynamic variables, which is enforced through the lower cutoff constraint is essential since otherwise the system dwindles over time. The non-conservation of the sum of the values of the dynamic variables is important as well. It allows to perform the multiplicative updating on a single variable at a time with no explicit binary interactions since a gain in $`w_i`$ does not require a corresponding immediate loss by other $`w_j`$’s. In fact, the interactions between the dynamic variables are implied only in the step of Eq. (4) in which the lower cutoff is imposed. The dynamic rule (1) can be described by a master equation for the probability distribution $`p(w)`$ of the form

$$p(w,t+1)p(w,t)=\frac{1}{N}\left[_\lambda \mathrm{\Pi }(\lambda )p(w/\lambda ,t)𝑑\lambda p(w,t)\right],$$
(5)   
where the $`1/N`$ factor takes into account the fact that only one of the $`w_i`$’s is updated in each time step. This description applies for the bulk of the distribution of the $`w_i`$’s but not in the vicinity of the lower cutoff where the step of Eq. (4) which is not taken into account by Eq. (5) may be dominant.

For the following analysis it is convenient to normalize the $`w_j`$’s according to

$$w_j(t)w_j(t)/\overline{w}(t),j=1,\mathrm{},N.$$
(6)   
As a result, the new average $`\overline{w}(t)`$ is normalized to

$$\overline{w}(t)=_c^Nwp(w,t)𝑑w=1,$$
(7)   
while $`_iw_i(t)=N\overline{w}=N`$. Performing this normalization step after each iteration removes the non stationary part of the distribution and amounts statistically to an overall multiplicative factor. This (time dependent) factor which represents a global inflation rate can be recorded at each step. It is convenient to represent the dynamics (5) on the logarithmic scale. In terms of the new variables

$$W_i=\mathrm{ln}w_i,$$
(8)   
Eq. (1) defines a random walk with steps of random size $`\mathrm{ln}\lambda `$:

$$W_i(t+1)=W_i(t)+\mathrm{ln}\lambda .$$
(9)   
The corresponding probability distribution $`P(W)`$ becomes

$$P(W)=e^Wp(e^W).$$
(10)   
In terms of $`P`$ and $`W`$, the master equation (5) becomes:

$`P(W,t+1)P(W,t)=`$ (11)   
$`{\displaystyle \frac{1}{N}}\left[{\displaystyle _\lambda }\mathrm{\Pi }(\lambda )P(W\mathrm{ln}\lambda ,t)𝑑\lambda P(W,t)\right].`$   
The asymptotic stationary solution, is found to be 

$$P(W)e^{\alpha W}.$$
(12)   
In terms of the original variable $`w_i`$, we get according to Eq. (10) a power law distribution:

$$p(w)=Kw^{1\alpha }.$$
(13)   
The value of the exponent $`\alpha `$ is determined by the normalization condition \[Eq. (7)\] divided by the probability normalization condition $`_c^Np(w,t)𝑑w=1`$ (in order to eliminate the constant factor $`K`$), which yields:

$$N=\frac{\alpha 1}{\alpha }\left[\frac{\left(\frac{c}{N}\right)^\alpha 1}{\left(\frac{c}{N}\right)^\alpha \left(\frac{c}{N}\right)}\right].$$
(14)   
The exponent $`\alpha `$ is given implicitly as a function of $`c`$ and $`N`$ by Eq. (14). We identify two regimes within $`0c<1`$ in which Eq. (14) can be simplified and $`\alpha `$ can be obtained explicitly. For a given $`N`$ and values of $`c`$ in the range $`1/\mathrm{ln}Nc<1`$ one obtains $`\alpha >1`$ as well as $`(c/N)^\alpha c/N1`$. Consequently, in this range, one can neglect the $`(c/N)^\alpha `$ terms in Eq. (14) to obtain to a good approximation

$$N=\frac{\alpha 1}{\alpha }\left[\frac{1}{\left(\frac{c}{N}\right)}\right].$$
(15)   
which gives the explicit, $`N`$-independent solution

$$\alpha \frac{1}{1c}.$$
(16)   
This relation is exact in the ”thermodynamic” limit $`N=\mathrm{}`$. The relation (16) has two remarkable properties: (a) it does not depend on the distribution $`\mathrm{\Pi }(\lambda )`$; (b) it gives rise to $`\alpha `$ values in the experimentally realistic range $`\alpha 1`$.

For finite $`N`$ and values of $`c`$ lower than $`1/\mathrm{ln}N`$ the approximation Eq. (16) breaks down and values $`\alpha <1`$ become possible. However, for any finite N, another approximation holds in the range $`c1/N<1`$. In this range $`(c/N)(c/N)^\alpha 1`$ and therefore one can neglect $`(c/N)^\alpha `$ in the numerator of Eq. (14) and $`c/N`$ in the denominator to obtain:

$$N=\frac{\alpha 1}{\alpha }\left[\frac{1}{\left(\frac{c}{N}\right)^\alpha }\right].$$
(17)   
By taking the logarithm on both sides and neglecting terms of order $`1`$ we obtain

$$\alpha \frac{\mathrm{ln}N}{\mathrm{ln}(N/c)}.$$
(18)   
Note that even for systems in which the lower bound (which is due to some microscopic discretization) given by $`c`$, is orders of magnitude smaller than $`1/N`$, the resulting $`\alpha `$ may differ significantly from the free multiplicative random walk result $`\alpha =0`$. Since $`c`$ enters in the formula (18) for $`\alpha `$ through its logarithm, the system gives away information on its microscopic scale cut-off $`c`$ through the exponent $`\alpha `$ of its macroscopic power law behavior.

One should emphasize that in the region where $`\alpha <1`$ the average $`\overline{w}`$ of the distribution $`p(w)`$ in Eq. (7) is not well defined and in fact one expects in the actual runs very wide macroscopic fluctuations of this mean. These fluctuations are however never infinite because according to the formulae above, as one increases the size of the system $`N`$, the region along the $`c`$ axis where $`\alpha <1`$ shrinks to 0. For $`1<\alpha <2`$ it is only the standard deviation of the distribution $`p(w)`$ which is formally divergent. This gives rise in the actual computer simulations to wide fluctuations of the individual values of $`w_i`$. However, this divergence is kept in check too by the fact that no $`w_i`$ can possibly exceed $`N\overline{w}`$, namely $`p(N\overline{w})=0`$. This amounts to a truncation from above of the power law Eq. (13).

## III Numerical Simulations and Results

Numerical simulations of the stochastic multiplicative process described by Eqs. (1) and (4), confirm the validity of Eq. (13) for a wide range of lower bounds $`c`$. It appears that the exponent $`\alpha `$ is largely independent of the shape of the probability distribution $`\mathrm{\Pi }(\lambda )`$. Fig. 1 shows the distribution of $`w_i`$, $`i=1,\mathrm{},N`$, obtained for $`N=1000`$, $`c=0.3`$, and $`\lambda `$ uniformly distributed in the range $`0.9\lambda 1.1`$. A power law distribution is found for a range of three decades between $`w_{min}=0.0003`$ and $`w_{max}=0.3`$. The slope of the best linear fit within this range is given by $`\alpha =1.4`$, in agreement with Eqs. (14) and (16). On the horizontal axis of this graph the sum of all $`w_i`$’s is normalized to 1 and therefore $`\overline{w}=0.001`$. The exponent $`\alpha `$ as a function of the lower cutoff $`c`$ is shown in Fig. 2. Numerical results are presented for $`N=100`$ (empty dots), $`1000`$ (full dots) and $`5000`$ (squares). The prediction of Eq. (14) is shown for $`N=1000`$ (solid line), which is in good agreement with the numerical results for all values of $`c`$. The approximate expression Eq. (16) is also shown (dashed line). It is observed that for $`N=1000`$ this approximation gradually starts to hold as $`c`$ is increased beyond $`1/\mathrm{ln}(1000)`$, in agreement with the theoretical analysis. In general, for a given $`N`$, $`\alpha `$ is monotonically increasing as a function of $`c`$, starting from $`\alpha =0`$ (which corresponds to $`1/w`$ distribution) at $`c=0`$, where the $`w_i`$’s are uncoupled. It is also observed that as $`N`$ is raised, the value of $`\alpha `$ which corresponds to a given $`c`$ increases monotonically. As a result, the range of validity $`1/\mathrm{ln}Nc<1`$ of the approximation Eq. (16) is extended and the knee adjacent to $`c=0`$ sharpens and becomes a discontinuity for $`N\mathrm{}`$. The range $`0c<1/N`$ in which the approximation of Eq. (18) is valid, shrinks correspondingly.

Let us turn now to the dynamics of the system as a whole. The dynamics of the system involves, according to Eq. (1), a generalized random walk with step sizes distributed according to Eq. (13). Therefore, the stochastic fluctuations of $`\overline{w}(t)`$ after $`\tau `$ time steps:

$$r(\tau )=\frac{\overline{w}(t+\tau )\overline{w}(t)}{\overline{w}(t)}$$
(19)   
are governed  by a truncated Lévy distribution $`L_\alpha (r)`$.

In Fig. 3 we show the distribution of the stochastic fluctuations $`r(\tau )`$ for $`\tau =50`$, which is given by a (truncated) Lévy distribution $`L_\alpha (r)`$. According to Ref. , the peak of the (truncated) Lévy-stable distribution scales with $`\tau `$ as

$$L_\alpha (r=0)\tau ^{1/\alpha }$$
(20)   
where $`\alpha `$ is the index of the Lévy distribution. In Fig. 4 we show the height of the peak P$`(r=0)`$ of Fig. 3 as a function of $`\tau `$. It is found that the slope of the fit in Fig. 4 is $`0.71`$, which following the scaling relation (20) means that the index of the Lévy distribution in Fig. 3 is $`\alpha =1/(0.71)=1.4`$. These results were obtained for the same parameters which gave rise to the power law distribution with $`\alpha =1.4`$ in Fig. 1. Thus, the prediction that the fluctuations of $`\overline{w}`$ in Fig. 3 follow a (truncated) Lévy-stable distribution with an index $`\alpha `$ which equals the exponent $`\alpha `$ of the power-law distribution in Fig. 1, is confirmed.

## IV Discussion

The model considered in this paper may be relevant to a variety of empirical systems in the physical, biological and social sciences which can be described by a set of interacting dynamic variables which follow a stochastic multiplicative dynamics. Such dynamical processes may play a role in the formation of the mass distribution in the universe where clusters of galaxies accumulate and eventually form super-clusters. In a different context, the growth of cities is basically a multiplicative process governed by the reproduction rate of the local population in addition to mobility between cities.

Enhanced diffusion processes, which can be described by the Lévy-stable distribution have been observed in a variety of nonlinear dynamical systems . Unlike the stochastic model studied here, these systems are governed by deterministic rules. They exhibit intermittent chaotic motion which gives rise to enhanced diffusion.

In population dynamics, the number of individuals in each specie varies stochastically from one season to the next with a multiplicative factor which depends on the local conditions. The lower bound may represent the minimal number of individuals required for the species to survive in the given environment. In this case the number of species may not be strictly a constant, but species that are wiped out may be replaced by others which invade their area. In this context it was found that the number of species of a given size often follows a decreasing power-law distribution as a function of their size (see e.g. Ref. ).

In the economic context of a stock-market system the dynamic variables $`w_i`$, $`i=1,\mathrm{},N`$ may represent the wealth of individual investors. In this case the dynamics represents the increase (or decrease) by a random factor $`\lambda (t)`$ of the wealth $`w_i`$ of the investor $`i`$ between times $`t`$ and $`t+1`$. The lower bound may represent a minimal wealth required in order to participate in stock market trading. In a more general economic model, this lower bound may be related to a basket of basic publicly funded services which every individual receives. In another possible interpretation, the $`w_i`$’s represent the capitalization (total market value) of the firm $`i`$, which may increase (or decrease) by a factor $`\lambda (t)`$ at each time step. In this case the lower bound may represent the minimal requirements for a company stock to be publicly traded.

Studies of the distribution of wealth in the general population revealed a power-law behavior (see e.g. Ref. ). More recently it was shown  that the distribution of individual wealth of the 400 richest people in the United States (Forbes 400) corresponds to a power law with $`\alpha =1.36`$ \[more precisely $`W(n)=Cn^{1/\alpha }`$ where $`W(n)`$ is the wealth of the $`n`$-th richest person on the list\]. Recent analysis of stock market returns, measured over many years found a truncated Lévy distribution $`L_\alpha (r)`$ with the index $`\alpha =1.4`$ for an extended (but finite) range of returns $`r`$ . These results indicate that the property observed in our model, namely that the same value of the index $`\alpha `$ appears both in the power law distribution and in the Lévy-stable distribution of the fluctuations may be of relevance in the economic context. To further explore this possibility it would be interesting to examine whether the distribution of total market values of companies in the stock market exhibits a power law behavior of the form (13) with $`\alpha =1.4`$.

## V Summary

We have studied a generic model of stochastic auto-catalytic dynamics of many degrees of freedom using computer simulations. The model consists of dynamic variables $`w_i`$, $`i=1,\mathrm{},N`$ which are updated randomly one at a time through an autocatalytic process at the individual level. In addition, the variables are coupled through a lower bound constraint which enhances the variables which fall below a fraction of the global average. The model may describe a large variety of systems such as stock markets and city populations. The distribution $`p(w,t)`$ of the system components $`w_i`$ turns out to fulfill a power law distribution of the form $`p(w,t)w^{1\alpha }`$. In the limit $`N=\mathrm{}`$, $`c0`$ one obtains the case often encountered in nature: $`\alpha 1`$. The average $`\overline{w}(t)`$ exhibits intermittent fluctuations following a Lévy-stable distribution with the same index $`\alpha `$. This relation between the distribution of system components and the temporal fluctuations of their average may be relevant to a variety of empirical systems. For example, it may provide a connection between the distribution of wealth/capitalization in a stock market and the distribution of the index fluctuations.