- Email: [email protected]
Production, depreciation and the size distribution of firms
Physica A 387 (2008) 3209–3217 www.elsevier.com/locate/physa
Production, depreciation and the size distribution of firms Qi Ma, Yongwang Chen, Hui Tong, Zengru Di ∗ Department of Systems Science, School of Management, Beijing Normal University, Beijing 100875, China Department of Systems Science, Center for Complexity Research, Beijing Normal University, Beijing 100875, China Received 3 July 2007; received in revised form 13 October 2007 Available online 17 January 2008
Abstract Many empirical researches indicate that firm size distributions in different industries or countries exhibit some similar characters. Among them the fact that many firm size distributions obey power-law especially for the upper end has been mostly discussed. Here we present an agent-based model to describe the evolution of manufacturing firms. Some basic economic behaviors are taken into account, which are production with decreasing marginal returns, preferential allocation of investments, and stochastic depreciation. The model gives a steady size distribution of firms which obey power-law. The effect of parameters on the power exponent is analyzed. The theoretical results are given based on both the Fokker–Planck equation and the Kesten process. They are well consistent with the numerical results. c 2008 Elsevier B.V. All rights reserved.
PACS: 89.65.Gh; 89.75.Da Keywords: Size distribution of firms; Power-law; Production and depreciation; Econophysics
1. Introduction Power-law distributions are widely observed in nature, especially in complex systems, such as the magnitude of earthquakes, moon craters, city populations, the frequency of use of words in many human languages, etc. It draws a great deal of attention from scientists, and a lot of mechanisms to generate power-law distribution have been discussed [1]. In the research of economic systems, power-law distribution is also a typical topic for empirical and theoretical studies. Actually, the Pareto law was first introduced by economists to describe country income distribution [2]. Recently, with the development of econophysics, much more evidence has been found for the powerlaw distribution, from the returns of financial market, welfare, personal income, to the size of firms and so on. In this paper, we concentrate on the mechanism of power-law distribution of firm sizes. In empirical studies, the size of a business firm is usually measured by the sales, the number of employees, the capital employed or the total assets. Firm size distribution is known to be extremely skewed, which has been generally described by a log-normal distribution since Gibrat [3], and the upper tail has been described by the Pareto or Zipf ∗ Corresponding author at: Department of Systems Science, School of Management, Beijing Normal University, Beijing 100875, China. Tel.: +86 10 58807876; fax: +86 10 58800141. E-mail address: [email protected] (Z. Di).
c 2008 Elsevier B.V. All rights reserved. 0378-4371/$ - see front matter doi:10.1016/j.physa.2008.01.064
3210
Q. Ma et al. / Physica A 387 (2008) 3209–3217
law [4,5]. Recent researches on firm size distribution in detail reveal that the firm size distribution can be fitted by power-law precisely with US, Japanese or European firms’ data [6–10]. Furthermore, the power exponents of the cumulative distribution are all around 1 (ranging from 0.7 to 1.2). Since the Pareto–Zipf law of firm size exhibits some similar characters with that of some physical systems which contain numerous interacting units, this interesting property attracts many scholars working on the perspectives of physics in accordance with these empirical facts [11–14]. Amaral et al. have studied firm size growth dynamics since 1997 [15] and developed a stochastic model, which is a complex system with large number of interacting units, and each unit has its own internal structures and independent growth dynamics. The model can fit well with the empirical result [16,17]. Takayasu’s research was based on firm size dynamics analogy to the aggregation-annihilation reaction model [18]. Axtell has insisted on using complexity approach to deal with the problem, and believed that agent-based modelling together with evolutionary dynamics should be helpful to understand the formation of powerlaw [19]. Besides, some other models have also been presented, which are based on the competition dynamics [20], the information transition, herd behaviors [21], and the proportional growth for the firms’ sizes [22]. On the base of multi-agent model, taken account of production and depreciation as firms’ basic economic behaviors, here we present a multi-agent model to describe the evolution of manufacturing firms. The system includes many individual manufacturing firms. Their sizes are characterized by capital stocks [7]. They give aggregate products by the Cobb–Douglas production function. Then each firm gets preferential returns proportional to its capital stock. And the new capital stock of each firm results both from this investment and a stochastic depreciation. According to the computer simulation, we got a power-law distribution of firm sizes. Moreover, We analyzed the effect of parameters on the power exponents and found that power exponents of the distribution do change along with the parameters. Then the model is simplified to a conserved system without influencing the final firm size distribution. Meanwhile, we analyzed the conserved model theoretically based on Fokker–Planck equation and Kesten process. The final distribution obeys power-law, which is well consistent with the numerical analysis. And the numerical power exponent of the firm size distribution changes the same way with the theoretical value when the parameters get changed. The paper is organized as follows. Section 2 describes the production-depreciation model and shows the final steady distribution of firm sizes by computer simulations. The relations of power exponent with the attractor and depreciation rate are also given. Then in Section 3, the model is simplified into a conserved one, i.e. the total capital stock of the system keeps constant with the evolution of firms. The Fokker–Planck equation for the evolution of firm size distribution is constructed. The final steady distribution is given, which obeys power-law for the upper end and the exponent is determined by the attractor and depreciation rate. The Kesten process gives similar results and they are all consistent well with the computer simulations. Some concluding remarks are given in Section 4. 2. The model and simulation results 2.1. The production and stochastic depreciation of firms In this paper, we consider the firms as separate agents of certain sizes. And the size is characterized by the capital a firm has. Considering a system consisting of N firms, firm i has the initial capital k0i , and k0i is a positive random PN i number. At time t, the total capital in the system is kt = i=1 kt . As the growth of capital is associated with the current amount of capital, so in short term, when the labor force is fixed, the capital has diminishing marginal product. Accordingly, we use the Cobb–Douglas function to present the gross products, namely yt = At kte . We let e = 0.4 in our following discussion, and At is a random variable denotes the production coefficient which including the effect of the labor and other realistic elements, like technology, input of other resources, and etc for simplicity. We assume that some fraction of the product is invested to increase the capital stock, and the fraction varies at time, namely 1kt = βt yt , where βt is a random variable. Because the production coefficient At is also a random variable, here for simplicity, we consider At has already carried the information about βt , as a result the capital increment can be simply written as 1kt = At kte . Then the kernel element that affects the final distribution should be how the incremental capital is allocated. Since firms with larger size contribute more to the gross product, we assume that the total incremental capital is allocated in proportion of one firm’s capital to the total capital. But in real economy, the whole income of each firm is not entirely decided by the capital it holds. So we introduce attractor ai to denote elements uncorrelated with the capital of firm i but influence its income. Thus, firm i’s capital is increased by preferential allocation with attractor ai , that is
Q. Ma et al. / Physica A 387 (2008) 3209–3217
3211
Fig. 1. (a) Time evolution of the total capital. (b) Steady state of firm size distribution P(k). Here time T = 30 000 and the distribution is the average of 10 simulations. At is uniformly distributed in [400, 600], e = 0.4, b = 0.2.
P j 1kti = 1kt (kti + ai )/ Nj=1 (kt + a j ). Obviously, the attractor describes other factors that affect the development of a firm besides its capital stock. It just works on the allocation of total investment. And also we consider the depreciation of the capital to make our model more close to the real economy. Here, we suppose the depreciation rate is δti , which is a random variable between 0 and b. Actually, this assumption for depreciation is far beyond the situation of natural depreciation of capital. It has taken the correlation and competition of firms in to account. Then for firm i, its growth is determined by: i kt+1
=
kti (1 − δti ) + 1kt (kti
+ ai )
X N
j
(kt + a j )
(1)
j=1
where, 1kt =
At kte
= At
N X
!e kti
.
(2)
i=1
2.2. Computer simulation results 2.2.1. Steady distribution and corresponding statistics We performed computer simulations with fixed number of firms N . Initially, the capital k0i of firm i is a random number between 0 and 10. The attractor ai is equal to 0.2 for all firms. Then we run the model to show the evolution of system. The results of the total capital stock of all firms and the final distribution are shown in Fig. 1(a) and (b). Fig. 1(a) shows that the total capital of all firms rises rapidly at first and get into steady state soon. Then the total capital does not increase any more and it fluctuates around a fixed value. This result indicates that the total capital in steady state is approximately conserved in average. According to Fig. 1(b), we can see the distribution of the firm sizes is basically a power-law, and with this set of parameters, the power exponent (it is the same as the exponent of cumulative distribution function P(k)) is approximately −1.1. Moreover, it is observed that the firm sizes distribution departs from power-law in the lower tail (where the firm sizes are small). And that is in accord with present empirical research. We noticed both the qualitative and quantitative features of the stable distributions did not change when the production coefficient At was a constant number equal to the mean value of the random variable, which means At can be replaced by the constant number 500 (mean value of the uniformly distributed At in [400, 600]), and will not affect the final results of the simulations. Therefore, we consider At as a constant number A in the following analysis and simulations for simplicity. Then, we adjust the parameters to find out the relations between parameters and the final distributions.
3212
Q. Ma et al. / Physica A 387 (2008) 3209–3217
Fig. 2. Power exponents of the density distribution versus values of the attractor. N = 10 000, A = 500, e = 0.4, b = 0.6, T = 30 000. We run the simulation 10 times and the displayed values are the averaged result.
2.2.2. Parameter effect on the power exponent (1) Effects of Attractors. Both the distribution and value of attractors may affect the final distribution of firm sizes. As to the distribution of attractors, we have tried several common distributions (uniform distribution, normal distribution, etc.). It has been found that the distribution of attractors has only little effect on the results. With the same mean value but different distributions of attractors (even with different deviations), the distributions of firm sizes are basically identical. So in the following discussion, we keep attractor as a same constant for every firm. When itPcomes to study how the value of attractor influences the distribution of firm sizes, we use another quantity N α = N a/ i=1 kti , the proportion of total attractors to total amount of capital in steady state, to reflect the influence of the attractors’ value on the firm sizes distribution. From the simulations, we can find that the final distribution of firm sizes sensitively depends on α. When the values of α lie between 0.001 and 0.1, the final stable firm sizes distributions are approximately power-law. The interval (0.001, 0.1) is related to the maximal depreciation rate b. When the attractors are too small or large, the firm sizes distribution is not close to power-law any more. It approaches exponential distribution when the attractors are large. While the attractors are extremely small, one firm may hold almost all the capital in the system. And other firms’ capital is basically close to zero. Fig. 2 shows how α affects the power exponent of the stable density function. When the maximal depreciation rate b is fixed, the exponent gets smaller as α increases in the effective range. (2) Effects of Depreciation rate. Just like the attractor, the depreciation rate also has two possible ways to affect the stable distribution of firm sizes. Different from the effects of attractors, we can figure out from the simulation that the stochastic properties of depreciation is crucial to the final distribution of firm sizes. Accordingly, we use uniform distribution as the distribution of depreciation rates for simplification in our numerical simulations, which means the depreciation rate δti of firm i is a random variable taken uniformly from the interval [0, b] at any time t. The character of uniform distribution makes that once the maximal depreciation rate b is determined, the mean value and variance is fixed correspondingly. So the maximal depreciation rate b can be used to reflect how the values of depreciation affect the stable firm sizes distributions. Fig. 3 shows the relation between the maximal depreciation rate b and the power exponent of the stable distribution. Generally, the power exponent and the maximal depreciation rate b are positive correlated. Namely when b is increased, the corresponding exponent increases as well. But this result cannot explain how the depreciation rate affects the final firm sizes distribution exactly. Because when other parameters are fixed, firm sizes will change along with the maximal depreciation b. And then α, the proportion of attractors to total capital, will change as well. Since we have already known that the stable firm sizes distribution is sensitively related to α, we cannot distinguish which factor is the essential reason that affects the final distribution at this time. We will solve the problem in the theoretical analysis of the simplified model.
Q. Ma et al. / Physica A 387 (2008) 3209–3217
3213
Fig. 3. Power exponents of the density distribution versus the maximal depreciation rate b. N = 10 000, A = 500, e = 0.4, α = 0.02, T = 30 000. We run the simulation 10 times and the displayed values are the averaged result.
3. Simplified model and analytical results 3.1. The model with conserved capital stock According to our initial model given inP last section, at time t, the total incremental capital in the system is N 1kt = Akte , and the total depreciation is i=1 δti kti . Because the depreciation is not correlated with the capital stock, when the number of firms N is large enough, the total amount depreciation at step t can be expressed as the Pof N mean value of depreciation rate multiplies the total capital, namely i=1 δti kti = δkt (δ denotes the mean value of δti ). e So the total size of all firms will be stable where δkt = Akt , that means when the model achieves equilibrium, there exists that: 1kt = Akte = δkt .
(3)
This will result in a final total capital stock k. According to the former simulations, we know that the distribution of attractors does not affect the stable distribution of firm sizes. So we assume that every firm’s attractor is all a. Thus when the model achieves equilibrium, there exists N X
(kti + a) = k + N a = k(1 + α)
(4)
i=1
here t is the time step of system evolution. k denotes the value of kt when the system is in equilibrium. Using Eqs. (3) and (4) in Eq. (1), one has: i kt+1 = kti [1 − δti + δ/(1 + α)] + aδ/(1 + α)
(5)
where δti is a random variable. From above discussion, for the final steady state, we have simplified the model into a conserved one and have used mean-field approximation to get Eq. (5). To make sure the simplified model and Eq. (5) have kept all the key characters of the original model, we run the simulation with the same set of parameters and then compared the results. For the conserved model, we modified the original model by removing the production function and let the total incremental capital equal to the total depreciation. This makes the system be precisely conserved. The comparison is shown in Fig. 4. Under the same set of parameters, the firm size distribution of the original model Eq. (1), the precisely conserved model and the Eq. (5) are given almost the same final distribution. From the comparison, we can deduce that the simplified models (both the conserved model and the Eq. (5)) have kept the kernel characters of the original model. The fluctuation of the total size basically has no influence on the final distribution. Therefore, the system in our model could be treated as a conserved one, which means the total capital is a constant value.
3214
Q. Ma et al. / Physica A 387 (2008) 3209–3217
Fig. 4. Comparison of the original model, simplified conserved model and Eq. (5). Parameters: N = 10 000, e = 0.4, b = 0.6 are the same for all models. The original model: same attractor ai = 0.2 for all firms; Simplified conserved model: every firm has the same initial firm size 23.4 at the first step (the value 23.4 was calculated according to the total firm size of the original model when it is in the stable state), ai = 0.2 for all firms; Eq. (5): a = 0.2, α = 0.008 54 (the value also came from the stable total capital of the original model).
3.2. Theoretical analysis and the results 3.2.1. Theoretical analysis based on the Fokker–Planck equation Neglecting the firm index i, Eq. (5) can be simply expressed as [23] kt+1 = kt θ + c + εt kt
(6)
with ¯ θ = 1 − α δ/(1 + α) = 1 − αb/(2 + 2α), m = b/2, D = b2 /12, ¯ + α) c = a¯ δ/(1
(7)
where εt is related with the stochastic depreciation. It is a random variable uniformly distributed in [−m, m] with variance D. c is a positive small constant term. From Eq. (6) we obtain kt+1 − kt = kt (θ − 1) + c + εt kt . Therefore, k˙ = k(θ − 1) + c + εt kt .
(8)
Let us assume that the probability density function of k is w(k, t) at time t, and p(k − 1k, 1k, τ ) represents the probability density of the particle moving 1k from k − 1k in time interval τ . For any k and τ , there is Z p(k, 1k, τ )d1k = 1. (9) The probability transfer function of the system is Z ω(k, t + τ ) = p(k − 1k, 1k, τ )ω(k − 1k, t)d1k.
(10)
Q. Ma et al. / Physica A 387 (2008) 3209–3217
3215
When the time interval τ is very small, we can expand ω(k, t + τ ) up to second order in terms of Taylor expansion and have ∂ω(k, t) τ + ◦(τ 2 ), ∂t ∂ p(k, 1k, τ ) 1 ∂ 2 p(k, 1k, τ ) p(k − 1k, 1k, τ ) = p(k, 1k, τ ) − 1k + (1k)2 + · · · , ∂k 2 ∂k 2 ∂ω(k, t) 1 ∂ 2 ω(k, t) (1k)2 + · · · . ω(k − 1k, t) = ω(k, t) − 1k + ∂k 2 ∂k 2 ω(k, t + τ ) = ω(k, t) +
Substituting the above expression for the corresponding terms in Eq. (10), using Eq. (9), and ignoring high order terms, we have ∂ω(k, t) 1 ∂ 1 1 ∂2 =− [ω(k, t)h1ki] + [ω(k, t)h(1k)2 i]. ∂t τ ∂k 2 τ ∂k 2
(11)
Refer to Eq. (8), when the time interval τ is small enough, k could be treated as a constant, then we can obtain Z τ ˙ = [k(θ − 1) + c]τ + k 1k = kτ εt dt. 0
Then the first and second order moment of 1k are Z τ εt dt = [k(θ − 1) + c]τ h1ki = [k(θ − 1) + c]τ + k
(12)
0
h(1k)2 i =
*
τ
Z [k(θ − 1) + c]τ + k
εt dt
2 +
= [k(θ − 1) + c]2 τ 2 + k 2 Dτ.
(13)
0
Using Eqs. (12) and (13) in Eq. (11), and for the final steady distribution (t → ∞) and τ → 0, we obtain 0=−
∂ 1 ∂2 [ωk 2 D]. {ω[k(θ − 1) + c]} + ∂k 2 ∂k 2
It gives (1 − θ )ω + [k(1 − θ ) − c]ω0 +
1 2 00 Dk ω + 2k Dω0 + Dω = 0. 2
Because θ is a constant number less than 1, we ignore c when k is sufficiently large and get 1 2 00 Dk ω + (1 − θ + 2D)kω0 + (1 − θ + D)ω = 0 2
(14)
Eq. (14) is a Cauchy–Euler equation, its general solution is ω = C1 k −2−2(1−θ)/D + C2 k −1 . Where ω is a probability density function. is 0, namely
R +∞ N
ωdk is supposed to be limited, so the coefficient of the second term
ω = Ck −2−2(1−θ)/D .
(15)
Using Eq. (7) in Eq. (15) we can get the power exponent Exponent = −2 −
12α . (1 + α)/b
(16)
3216
Q. Ma et al. / Physica A 387 (2008) 3209–3217
Fig. 5. Comparisons of analytical solutions given by the Fokker–Planck equation and Kesten process with simulations given by the original and conserved models and the numerical solution of Eq. (5). N = 10 000, A = 500, e = 0.4; For (a): b = 0.6; For (b): α = 0.02.
3.2.2. Theoretical analysis based on the Kesten process For the conserved model, δ is the only random variable in Eq. (5). Then we can treat Eq. (5) as a Kesten Process [24, 25], which is generally written as: x(t + 1) = x(t)a(t) + b(t)
(17)
where both a(t) and b(t) are positive random numbers. Providing hlog a(t)i < 0 and other conditions, x(t) obeys power-law distributed of the form [24]: p(x(t)) ∼ x(t)−(1+µ)
(18)
with µ given by ha(t)µ i = 1. ¯ In our model, for Eq. (5), we can find out that a(t) is resulted from the term [1 − δti + δ/(1 + α)] and we have i ¯ ¯ E[1 − δt + δ/(1 + α)] < 1. b(t) is related to the term a¯ δ/(1 + α), which is a constant. Those make the Kesten process — Eq. (5) has a solution. As the a(t) in our model is uniform distributed, the average value is θ = 1 − αb/(2 + 2α), and a(t) fluctuates in the interval [θ −b/2, θ +b/2]. According to ha(t)µ i = 1, we know that the parameter µ satisfies the relation: (θ + b/2)µ+1 − (θ − b/2)µ+1 = b(µ + 1).
(19)
Therefore, when the system is in stable state, the power exponent of the firm size distribution is: Exponent = −1 − µ
(20)
where µ could be given by numerical solutions of Eq. (19). 3.3. Comparison of theoretical and numerical results From Fig. 5 we can observe that the analytical solutions are generally the same as the simulation results. But there still exist some tolerable differences between them. As the numerical power exponents are universally larger than the analytical ones, it indicates the existence of system error. The system error may be caused by the simulation we use or the approximation in the process of solving the Langevin and Kesten equation. Also we can find from Fig. 5 that when α is small or b increases, the error was comparatively less. In the simulations, when α is small or b is large, the stable distribution is more close to power-law. 4. Conclusion and discussion This article focused on the power-law distribution of firm sizes. Considering the elementary economic behaviors of firms, we presented an evolutionary model of manufacturing firms based on the multi-agent model and obtained a
Q. Ma et al. / Physica A 387 (2008) 3209–3217
3217
power-law distribution. The result of our model is consistent with lots of empirical researches. According to our study, we believe that random factors given by the stochastic depreciation and preferential allocation are the crucial causes of the power-law distribution. In our model, preferential allocation is represented by the mechanism that a larger firm can get more capital from the total incremental capital. This is a universal mechanism in both nature and social life. Due to this, our model can explain not only the power-law distribution of firm sizes but also power-laws arising in other fields. In steady state, the value of power exponents is more sensitively dependent on the rescaled attractor than the maximal depreciation rate. The distribution of attractors has nearly no effect on the stable distribution of firm sizes. But for the depreciation rate, both its mean value and variance can affect the final size distribution of firms. We applied some methods of statistical physics in our study. With respect to the Langevin equation and by a series of simplifications, we worked out the analytical solution of the model and it generally accorded with the computer simulation. We also applied the Kesten process and obtained the corresponding analytical solution. It is nearly the same as the result of the Fokker–Planck equation and consistent with the simulation as well. To sum up, we find that random factors and preferential allocation are the two main feasible explanations for the power-law distribution in our model. As the two terms appear a great deal in many fields of nature or social life, they can be widely used for explaining power-law distributions in these fields. Acknowledgments We are thankful to an anonymous referee for the useful suggestions. This research was supported partially by the National Science Foundation of China under Grant No. 70371072 and No. 70601002. References [1] [2] [3] [4] [5] [6] [7] [8]
[9] [10] [11] [12] [13] [14] [15]
[16] [17] [18] [19] [20] [21] [22] [23] [24] [25]
M.E.J. Newman, Contemp. Phys. 46 (2005) 323–351. V. Pareto, Cours d’Economie Politique, Librairie Droz, Geneva, 1987. R. Gibrat, Les InT`egali´es Economiques, Sirey, Paris, 1931. Y. Ijiri, H.A. Simon, Skew Distributions and the Size of Business Firms, North-Holland, Amsterdam, 1977. R. Lucas, Bell J. Econom. 9 (1978) 508. M.H.R. Stanley, S.V. Buldyrev, S. Havlin, R.N. Mantegna, M. Salinger, Econom. Lett. 49 (1995) 453–457. E. Gaffeo, M. Gallegati, A. Palestrini, Physica A 324 (2003) 117–123; C. Di Guilmi, M. Gallegati, P. Ormerodc, Physica A 334 (2004) 267–273. Y. Fujiwara, C. Di Guilmi, H. Aoyama, M. Gallegati, W. Souma, Physica A 335 (2004) 197–216; Y. Fujiwara, Physica A 337 (2004) 219–230; Y. Fujiwara, H. Aoyama, C. Di Guilmi, W. Souma, M. Gallegati, Physica A 344 (2004) 112–116. J.J. Ramsden, Gy. Kiss-Haypl, Physica A 277 (2000) 220–227. R. Axtell, Science 293 (2001) 1819. R.N. Mantegna, H.E. Stanley, An Introduction to Econophysics, Cambridge University Press, Cambridge, 2000. P.F. Peretto, European Economic Rev. 43 (1999) 1747–1773. F. Hashemi, J. Evol. Econ. 10 (2000) 507–521. S.G. Winter, Y.M. Kaniovski, G. Dosi, J. Evol. Econ. 13 (2003) 355–383. L.A.N. Amaral, S.V. Buldyrev, H. Leschhorn, P. Maass, M.A. Salinger, H.E. Stanley, M.H.R. Stanley, J. Phys. I France 7 (1997) 521–633; S.V. Buldyrev, L.A.N. Amaral, S. Havlin, H. Leschhorn, P. Maass, M.A. Salinger, H.E. Stanley, M.H.R. Stanley, J. Phys. I France 7 (1997) 635–650. L.A.N. Amaral, S.V. Buldyrev, S. Havlin, M.A. Salinger, H.E. Stanley, Phys. Rev. Lett. 80 (1998) 1385; Y. Lee, L.A.N. Amaral, D. Canning, M. Meyer, H.E. Stanley, Phys. Rev. Lett. 81 (1998) 3275. L.A.N. Amaral, P. Gopikrishnan, V. Plerou, H.E. Stanley, Physica A 299 (2001) 127–136. H. Takayasu, K. Okuyama, Fractals 6 (1998) 67–79. R. Axtell, CSED Working Paper No. 3, Brookings Institution, 2001. H.M. Gupta, J.R. Campanha, Physica A 323 (2003) 626–634. D. Zheng, G.J. Rodgersa, P.M. Huic, Physica A 310 (2002) 480–486. G. De Fabritiis, F. Pammolli, M. Riccaboni, Physica A 324 (2003) 38–44. D. Sornette, R. Cont, J. Phys. I France 7 (1997) 431–444. H. Kesten, Acta Math. 131 (1973) 207–248. H. Takayasu, A.-H. Sato, M. Takayasu, Phys. Rev. Lett. 79 (1997) 966–969.
Production, depreciation and the size distribution of firms Qi Ma, Yongwang Chen, Hui Tong, Zengru Di ∗ Department of Systems Science, School of Management, Beijing Normal University, Beijing 100875, China Department of Systems Science, Center for Complexity Research, Beijing Normal University, Beijing 100875, China Received 3 July 2007; received in revised form 13 October 2007 Available online 17 January 2008
Abstract Many empirical researches indicate that firm size distributions in different industries or countries exhibit some similar characters. Among them the fact that many firm size distributions obey power-law especially for the upper end has been mostly discussed. Here we present an agent-based model to describe the evolution of manufacturing firms. Some basic economic behaviors are taken into account, which are production with decreasing marginal returns, preferential allocation of investments, and stochastic depreciation. The model gives a steady size distribution of firms which obey power-law. The effect of parameters on the power exponent is analyzed. The theoretical results are given based on both the Fokker–Planck equation and the Kesten process. They are well consistent with the numerical results. c 2008 Elsevier B.V. All rights reserved.
PACS: 89.65.Gh; 89.75.Da Keywords: Size distribution of firms; Power-law; Production and depreciation; Econophysics
1. Introduction Power-law distributions are widely observed in nature, especially in complex systems, such as the magnitude of earthquakes, moon craters, city populations, the frequency of use of words in many human languages, etc. It draws a great deal of attention from scientists, and a lot of mechanisms to generate power-law distribution have been discussed [1]. In the research of economic systems, power-law distribution is also a typical topic for empirical and theoretical studies. Actually, the Pareto law was first introduced by economists to describe country income distribution [2]. Recently, with the development of econophysics, much more evidence has been found for the powerlaw distribution, from the returns of financial market, welfare, personal income, to the size of firms and so on. In this paper, we concentrate on the mechanism of power-law distribution of firm sizes. In empirical studies, the size of a business firm is usually measured by the sales, the number of employees, the capital employed or the total assets. Firm size distribution is known to be extremely skewed, which has been generally described by a log-normal distribution since Gibrat [3], and the upper tail has been described by the Pareto or Zipf ∗ Corresponding author at: Department of Systems Science, School of Management, Beijing Normal University, Beijing 100875, China. Tel.: +86 10 58807876; fax: +86 10 58800141. E-mail address: [email protected] (Z. Di).
c 2008 Elsevier B.V. All rights reserved. 0378-4371/$ - see front matter doi:10.1016/j.physa.2008.01.064
3210
Q. Ma et al. / Physica A 387 (2008) 3209–3217
law [4,5]. Recent researches on firm size distribution in detail reveal that the firm size distribution can be fitted by power-law precisely with US, Japanese or European firms’ data [6–10]. Furthermore, the power exponents of the cumulative distribution are all around 1 (ranging from 0.7 to 1.2). Since the Pareto–Zipf law of firm size exhibits some similar characters with that of some physical systems which contain numerous interacting units, this interesting property attracts many scholars working on the perspectives of physics in accordance with these empirical facts [11–14]. Amaral et al. have studied firm size growth dynamics since 1997 [15] and developed a stochastic model, which is a complex system with large number of interacting units, and each unit has its own internal structures and independent growth dynamics. The model can fit well with the empirical result [16,17]. Takayasu’s research was based on firm size dynamics analogy to the aggregation-annihilation reaction model [18]. Axtell has insisted on using complexity approach to deal with the problem, and believed that agent-based modelling together with evolutionary dynamics should be helpful to understand the formation of powerlaw [19]. Besides, some other models have also been presented, which are based on the competition dynamics [20], the information transition, herd behaviors [21], and the proportional growth for the firms’ sizes [22]. On the base of multi-agent model, taken account of production and depreciation as firms’ basic economic behaviors, here we present a multi-agent model to describe the evolution of manufacturing firms. The system includes many individual manufacturing firms. Their sizes are characterized by capital stocks [7]. They give aggregate products by the Cobb–Douglas production function. Then each firm gets preferential returns proportional to its capital stock. And the new capital stock of each firm results both from this investment and a stochastic depreciation. According to the computer simulation, we got a power-law distribution of firm sizes. Moreover, We analyzed the effect of parameters on the power exponents and found that power exponents of the distribution do change along with the parameters. Then the model is simplified to a conserved system without influencing the final firm size distribution. Meanwhile, we analyzed the conserved model theoretically based on Fokker–Planck equation and Kesten process. The final distribution obeys power-law, which is well consistent with the numerical analysis. And the numerical power exponent of the firm size distribution changes the same way with the theoretical value when the parameters get changed. The paper is organized as follows. Section 2 describes the production-depreciation model and shows the final steady distribution of firm sizes by computer simulations. The relations of power exponent with the attractor and depreciation rate are also given. Then in Section 3, the model is simplified into a conserved one, i.e. the total capital stock of the system keeps constant with the evolution of firms. The Fokker–Planck equation for the evolution of firm size distribution is constructed. The final steady distribution is given, which obeys power-law for the upper end and the exponent is determined by the attractor and depreciation rate. The Kesten process gives similar results and they are all consistent well with the computer simulations. Some concluding remarks are given in Section 4. 2. The model and simulation results 2.1. The production and stochastic depreciation of firms In this paper, we consider the firms as separate agents of certain sizes. And the size is characterized by the capital a firm has. Considering a system consisting of N firms, firm i has the initial capital k0i , and k0i is a positive random PN i number. At time t, the total capital in the system is kt = i=1 kt . As the growth of capital is associated with the current amount of capital, so in short term, when the labor force is fixed, the capital has diminishing marginal product. Accordingly, we use the Cobb–Douglas function to present the gross products, namely yt = At kte . We let e = 0.4 in our following discussion, and At is a random variable denotes the production coefficient which including the effect of the labor and other realistic elements, like technology, input of other resources, and etc for simplicity. We assume that some fraction of the product is invested to increase the capital stock, and the fraction varies at time, namely 1kt = βt yt , where βt is a random variable. Because the production coefficient At is also a random variable, here for simplicity, we consider At has already carried the information about βt , as a result the capital increment can be simply written as 1kt = At kte . Then the kernel element that affects the final distribution should be how the incremental capital is allocated. Since firms with larger size contribute more to the gross product, we assume that the total incremental capital is allocated in proportion of one firm’s capital to the total capital. But in real economy, the whole income of each firm is not entirely decided by the capital it holds. So we introduce attractor ai to denote elements uncorrelated with the capital of firm i but influence its income. Thus, firm i’s capital is increased by preferential allocation with attractor ai , that is
Q. Ma et al. / Physica A 387 (2008) 3209–3217
3211
Fig. 1. (a) Time evolution of the total capital. (b) Steady state of firm size distribution P(k). Here time T = 30 000 and the distribution is the average of 10 simulations. At is uniformly distributed in [400, 600], e = 0.4, b = 0.2.
P j 1kti = 1kt (kti + ai )/ Nj=1 (kt + a j ). Obviously, the attractor describes other factors that affect the development of a firm besides its capital stock. It just works on the allocation of total investment. And also we consider the depreciation of the capital to make our model more close to the real economy. Here, we suppose the depreciation rate is δti , which is a random variable between 0 and b. Actually, this assumption for depreciation is far beyond the situation of natural depreciation of capital. It has taken the correlation and competition of firms in to account. Then for firm i, its growth is determined by: i kt+1
=
kti (1 − δti ) + 1kt (kti
+ ai )
X N
j
(kt + a j )
(1)
j=1
where, 1kt =
At kte
= At
N X
!e kti
.
(2)
i=1
2.2. Computer simulation results 2.2.1. Steady distribution and corresponding statistics We performed computer simulations with fixed number of firms N . Initially, the capital k0i of firm i is a random number between 0 and 10. The attractor ai is equal to 0.2 for all firms. Then we run the model to show the evolution of system. The results of the total capital stock of all firms and the final distribution are shown in Fig. 1(a) and (b). Fig. 1(a) shows that the total capital of all firms rises rapidly at first and get into steady state soon. Then the total capital does not increase any more and it fluctuates around a fixed value. This result indicates that the total capital in steady state is approximately conserved in average. According to Fig. 1(b), we can see the distribution of the firm sizes is basically a power-law, and with this set of parameters, the power exponent (it is the same as the exponent of cumulative distribution function P(k)) is approximately −1.1. Moreover, it is observed that the firm sizes distribution departs from power-law in the lower tail (where the firm sizes are small). And that is in accord with present empirical research. We noticed both the qualitative and quantitative features of the stable distributions did not change when the production coefficient At was a constant number equal to the mean value of the random variable, which means At can be replaced by the constant number 500 (mean value of the uniformly distributed At in [400, 600]), and will not affect the final results of the simulations. Therefore, we consider At as a constant number A in the following analysis and simulations for simplicity. Then, we adjust the parameters to find out the relations between parameters and the final distributions.
3212
Q. Ma et al. / Physica A 387 (2008) 3209–3217
Fig. 2. Power exponents of the density distribution versus values of the attractor. N = 10 000, A = 500, e = 0.4, b = 0.6, T = 30 000. We run the simulation 10 times and the displayed values are the averaged result.
2.2.2. Parameter effect on the power exponent (1) Effects of Attractors. Both the distribution and value of attractors may affect the final distribution of firm sizes. As to the distribution of attractors, we have tried several common distributions (uniform distribution, normal distribution, etc.). It has been found that the distribution of attractors has only little effect on the results. With the same mean value but different distributions of attractors (even with different deviations), the distributions of firm sizes are basically identical. So in the following discussion, we keep attractor as a same constant for every firm. When itPcomes to study how the value of attractor influences the distribution of firm sizes, we use another quantity N α = N a/ i=1 kti , the proportion of total attractors to total amount of capital in steady state, to reflect the influence of the attractors’ value on the firm sizes distribution. From the simulations, we can find that the final distribution of firm sizes sensitively depends on α. When the values of α lie between 0.001 and 0.1, the final stable firm sizes distributions are approximately power-law. The interval (0.001, 0.1) is related to the maximal depreciation rate b. When the attractors are too small or large, the firm sizes distribution is not close to power-law any more. It approaches exponential distribution when the attractors are large. While the attractors are extremely small, one firm may hold almost all the capital in the system. And other firms’ capital is basically close to zero. Fig. 2 shows how α affects the power exponent of the stable density function. When the maximal depreciation rate b is fixed, the exponent gets smaller as α increases in the effective range. (2) Effects of Depreciation rate. Just like the attractor, the depreciation rate also has two possible ways to affect the stable distribution of firm sizes. Different from the effects of attractors, we can figure out from the simulation that the stochastic properties of depreciation is crucial to the final distribution of firm sizes. Accordingly, we use uniform distribution as the distribution of depreciation rates for simplification in our numerical simulations, which means the depreciation rate δti of firm i is a random variable taken uniformly from the interval [0, b] at any time t. The character of uniform distribution makes that once the maximal depreciation rate b is determined, the mean value and variance is fixed correspondingly. So the maximal depreciation rate b can be used to reflect how the values of depreciation affect the stable firm sizes distributions. Fig. 3 shows the relation between the maximal depreciation rate b and the power exponent of the stable distribution. Generally, the power exponent and the maximal depreciation rate b are positive correlated. Namely when b is increased, the corresponding exponent increases as well. But this result cannot explain how the depreciation rate affects the final firm sizes distribution exactly. Because when other parameters are fixed, firm sizes will change along with the maximal depreciation b. And then α, the proportion of attractors to total capital, will change as well. Since we have already known that the stable firm sizes distribution is sensitively related to α, we cannot distinguish which factor is the essential reason that affects the final distribution at this time. We will solve the problem in the theoretical analysis of the simplified model.
Q. Ma et al. / Physica A 387 (2008) 3209–3217
3213
Fig. 3. Power exponents of the density distribution versus the maximal depreciation rate b. N = 10 000, A = 500, e = 0.4, α = 0.02, T = 30 000. We run the simulation 10 times and the displayed values are the averaged result.
3. Simplified model and analytical results 3.1. The model with conserved capital stock According to our initial model given inP last section, at time t, the total incremental capital in the system is N 1kt = Akte , and the total depreciation is i=1 δti kti . Because the depreciation is not correlated with the capital stock, when the number of firms N is large enough, the total amount depreciation at step t can be expressed as the Pof N mean value of depreciation rate multiplies the total capital, namely i=1 δti kti = δkt (δ denotes the mean value of δti ). e So the total size of all firms will be stable where δkt = Akt , that means when the model achieves equilibrium, there exists that: 1kt = Akte = δkt .
(3)
This will result in a final total capital stock k. According to the former simulations, we know that the distribution of attractors does not affect the stable distribution of firm sizes. So we assume that every firm’s attractor is all a. Thus when the model achieves equilibrium, there exists N X
(kti + a) = k + N a = k(1 + α)
(4)
i=1
here t is the time step of system evolution. k denotes the value of kt when the system is in equilibrium. Using Eqs. (3) and (4) in Eq. (1), one has: i kt+1 = kti [1 − δti + δ/(1 + α)] + aδ/(1 + α)
(5)
where δti is a random variable. From above discussion, for the final steady state, we have simplified the model into a conserved one and have used mean-field approximation to get Eq. (5). To make sure the simplified model and Eq. (5) have kept all the key characters of the original model, we run the simulation with the same set of parameters and then compared the results. For the conserved model, we modified the original model by removing the production function and let the total incremental capital equal to the total depreciation. This makes the system be precisely conserved. The comparison is shown in Fig. 4. Under the same set of parameters, the firm size distribution of the original model Eq. (1), the precisely conserved model and the Eq. (5) are given almost the same final distribution. From the comparison, we can deduce that the simplified models (both the conserved model and the Eq. (5)) have kept the kernel characters of the original model. The fluctuation of the total size basically has no influence on the final distribution. Therefore, the system in our model could be treated as a conserved one, which means the total capital is a constant value.
3214
Q. Ma et al. / Physica A 387 (2008) 3209–3217
Fig. 4. Comparison of the original model, simplified conserved model and Eq. (5). Parameters: N = 10 000, e = 0.4, b = 0.6 are the same for all models. The original model: same attractor ai = 0.2 for all firms; Simplified conserved model: every firm has the same initial firm size 23.4 at the first step (the value 23.4 was calculated according to the total firm size of the original model when it is in the stable state), ai = 0.2 for all firms; Eq. (5): a = 0.2, α = 0.008 54 (the value also came from the stable total capital of the original model).
3.2. Theoretical analysis and the results 3.2.1. Theoretical analysis based on the Fokker–Planck equation Neglecting the firm index i, Eq. (5) can be simply expressed as [23] kt+1 = kt θ + c + εt kt
(6)
with ¯ θ = 1 − α δ/(1 + α) = 1 − αb/(2 + 2α), m = b/2, D = b2 /12, ¯ + α) c = a¯ δ/(1
(7)
where εt is related with the stochastic depreciation. It is a random variable uniformly distributed in [−m, m] with variance D. c is a positive small constant term. From Eq. (6) we obtain kt+1 − kt = kt (θ − 1) + c + εt kt . Therefore, k˙ = k(θ − 1) + c + εt kt .
(8)
Let us assume that the probability density function of k is w(k, t) at time t, and p(k − 1k, 1k, τ ) represents the probability density of the particle moving 1k from k − 1k in time interval τ . For any k and τ , there is Z p(k, 1k, τ )d1k = 1. (9) The probability transfer function of the system is Z ω(k, t + τ ) = p(k − 1k, 1k, τ )ω(k − 1k, t)d1k.
(10)
Q. Ma et al. / Physica A 387 (2008) 3209–3217
3215
When the time interval τ is very small, we can expand ω(k, t + τ ) up to second order in terms of Taylor expansion and have ∂ω(k, t) τ + ◦(τ 2 ), ∂t ∂ p(k, 1k, τ ) 1 ∂ 2 p(k, 1k, τ ) p(k − 1k, 1k, τ ) = p(k, 1k, τ ) − 1k + (1k)2 + · · · , ∂k 2 ∂k 2 ∂ω(k, t) 1 ∂ 2 ω(k, t) (1k)2 + · · · . ω(k − 1k, t) = ω(k, t) − 1k + ∂k 2 ∂k 2 ω(k, t + τ ) = ω(k, t) +
Substituting the above expression for the corresponding terms in Eq. (10), using Eq. (9), and ignoring high order terms, we have ∂ω(k, t) 1 ∂ 1 1 ∂2 =− [ω(k, t)h1ki] + [ω(k, t)h(1k)2 i]. ∂t τ ∂k 2 τ ∂k 2
(11)
Refer to Eq. (8), when the time interval τ is small enough, k could be treated as a constant, then we can obtain Z τ ˙ = [k(θ − 1) + c]τ + k 1k = kτ εt dt. 0
Then the first and second order moment of 1k are Z τ εt dt = [k(θ − 1) + c]τ h1ki = [k(θ − 1) + c]τ + k
(12)
0
h(1k)2 i =
*
τ
Z [k(θ − 1) + c]τ + k
εt dt
2 +
= [k(θ − 1) + c]2 τ 2 + k 2 Dτ.
(13)
0
Using Eqs. (12) and (13) in Eq. (11), and for the final steady distribution (t → ∞) and τ → 0, we obtain 0=−
∂ 1 ∂2 [ωk 2 D]. {ω[k(θ − 1) + c]} + ∂k 2 ∂k 2
It gives (1 − θ )ω + [k(1 − θ ) − c]ω0 +
1 2 00 Dk ω + 2k Dω0 + Dω = 0. 2
Because θ is a constant number less than 1, we ignore c when k is sufficiently large and get 1 2 00 Dk ω + (1 − θ + 2D)kω0 + (1 − θ + D)ω = 0 2
(14)
Eq. (14) is a Cauchy–Euler equation, its general solution is ω = C1 k −2−2(1−θ)/D + C2 k −1 . Where ω is a probability density function. is 0, namely
R +∞ N
ωdk is supposed to be limited, so the coefficient of the second term
ω = Ck −2−2(1−θ)/D .
(15)
Using Eq. (7) in Eq. (15) we can get the power exponent Exponent = −2 −
12α . (1 + α)/b
(16)
3216
Q. Ma et al. / Physica A 387 (2008) 3209–3217
Fig. 5. Comparisons of analytical solutions given by the Fokker–Planck equation and Kesten process with simulations given by the original and conserved models and the numerical solution of Eq. (5). N = 10 000, A = 500, e = 0.4; For (a): b = 0.6; For (b): α = 0.02.
3.2.2. Theoretical analysis based on the Kesten process For the conserved model, δ is the only random variable in Eq. (5). Then we can treat Eq. (5) as a Kesten Process [24, 25], which is generally written as: x(t + 1) = x(t)a(t) + b(t)
(17)
where both a(t) and b(t) are positive random numbers. Providing hlog a(t)i < 0 and other conditions, x(t) obeys power-law distributed of the form [24]: p(x(t)) ∼ x(t)−(1+µ)
(18)
with µ given by ha(t)µ i = 1. ¯ In our model, for Eq. (5), we can find out that a(t) is resulted from the term [1 − δti + δ/(1 + α)] and we have i ¯ ¯ E[1 − δt + δ/(1 + α)] < 1. b(t) is related to the term a¯ δ/(1 + α), which is a constant. Those make the Kesten process — Eq. (5) has a solution. As the a(t) in our model is uniform distributed, the average value is θ = 1 − αb/(2 + 2α), and a(t) fluctuates in the interval [θ −b/2, θ +b/2]. According to ha(t)µ i = 1, we know that the parameter µ satisfies the relation: (θ + b/2)µ+1 − (θ − b/2)µ+1 = b(µ + 1).
(19)
Therefore, when the system is in stable state, the power exponent of the firm size distribution is: Exponent = −1 − µ
(20)
where µ could be given by numerical solutions of Eq. (19). 3.3. Comparison of theoretical and numerical results From Fig. 5 we can observe that the analytical solutions are generally the same as the simulation results. But there still exist some tolerable differences between them. As the numerical power exponents are universally larger than the analytical ones, it indicates the existence of system error. The system error may be caused by the simulation we use or the approximation in the process of solving the Langevin and Kesten equation. Also we can find from Fig. 5 that when α is small or b increases, the error was comparatively less. In the simulations, when α is small or b is large, the stable distribution is more close to power-law. 4. Conclusion and discussion This article focused on the power-law distribution of firm sizes. Considering the elementary economic behaviors of firms, we presented an evolutionary model of manufacturing firms based on the multi-agent model and obtained a
Q. Ma et al. / Physica A 387 (2008) 3209–3217
3217
power-law distribution. The result of our model is consistent with lots of empirical researches. According to our study, we believe that random factors given by the stochastic depreciation and preferential allocation are the crucial causes of the power-law distribution. In our model, preferential allocation is represented by the mechanism that a larger firm can get more capital from the total incremental capital. This is a universal mechanism in both nature and social life. Due to this, our model can explain not only the power-law distribution of firm sizes but also power-laws arising in other fields. In steady state, the value of power exponents is more sensitively dependent on the rescaled attractor than the maximal depreciation rate. The distribution of attractors has nearly no effect on the stable distribution of firm sizes. But for the depreciation rate, both its mean value and variance can affect the final size distribution of firms. We applied some methods of statistical physics in our study. With respect to the Langevin equation and by a series of simplifications, we worked out the analytical solution of the model and it generally accorded with the computer simulation. We also applied the Kesten process and obtained the corresponding analytical solution. It is nearly the same as the result of the Fokker–Planck equation and consistent with the simulation as well. To sum up, we find that random factors and preferential allocation are the two main feasible explanations for the power-law distribution in our model. As the two terms appear a great deal in many fields of nature or social life, they can be widely used for explaining power-law distributions in these fields. Acknowledgments We are thankful to an anonymous referee for the useful suggestions. This research was supported partially by the National Science Foundation of China under Grant No. 70371072 and No. 70601002. References [1] [2] [3] [4] [5] [6] [7] [8]
[9] [10] [11] [12] [13] [14] [15]
[16] [17] [18] [19] [20] [21] [22] [23] [24] [25]
M.E.J. Newman, Contemp. Phys. 46 (2005) 323–351. V. Pareto, Cours d’Economie Politique, Librairie Droz, Geneva, 1987. R. Gibrat, Les InT`egali´es Economiques, Sirey, Paris, 1931. Y. Ijiri, H.A. Simon, Skew Distributions and the Size of Business Firms, North-Holland, Amsterdam, 1977. R. Lucas, Bell J. Econom. 9 (1978) 508. M.H.R. Stanley, S.V. Buldyrev, S. Havlin, R.N. Mantegna, M. Salinger, Econom. Lett. 49 (1995) 453–457. E. Gaffeo, M. Gallegati, A. Palestrini, Physica A 324 (2003) 117–123; C. Di Guilmi, M. Gallegati, P. Ormerodc, Physica A 334 (2004) 267–273. Y. Fujiwara, C. Di Guilmi, H. Aoyama, M. Gallegati, W. Souma, Physica A 335 (2004) 197–216; Y. Fujiwara, Physica A 337 (2004) 219–230; Y. Fujiwara, H. Aoyama, C. Di Guilmi, W. Souma, M. Gallegati, Physica A 344 (2004) 112–116. J.J. Ramsden, Gy. Kiss-Haypl, Physica A 277 (2000) 220–227. R. Axtell, Science 293 (2001) 1819. R.N. Mantegna, H.E. Stanley, An Introduction to Econophysics, Cambridge University Press, Cambridge, 2000. P.F. Peretto, European Economic Rev. 43 (1999) 1747–1773. F. Hashemi, J. Evol. Econ. 10 (2000) 507–521. S.G. Winter, Y.M. Kaniovski, G. Dosi, J. Evol. Econ. 13 (2003) 355–383. L.A.N. Amaral, S.V. Buldyrev, H. Leschhorn, P. Maass, M.A. Salinger, H.E. Stanley, M.H.R. Stanley, J. Phys. I France 7 (1997) 521–633; S.V. Buldyrev, L.A.N. Amaral, S. Havlin, H. Leschhorn, P. Maass, M.A. Salinger, H.E. Stanley, M.H.R. Stanley, J. Phys. I France 7 (1997) 635–650. L.A.N. Amaral, S.V. Buldyrev, S. Havlin, M.A. Salinger, H.E. Stanley, Phys. Rev. Lett. 80 (1998) 1385; Y. Lee, L.A.N. Amaral, D. Canning, M. Meyer, H.E. Stanley, Phys. Rev. Lett. 81 (1998) 3275. L.A.N. Amaral, P. Gopikrishnan, V. Plerou, H.E. Stanley, Physica A 299 (2001) 127–136. H. Takayasu, K. Okuyama, Fractals 6 (1998) 67–79. R. Axtell, CSED Working Paper No. 3, Brookings Institution, 2001. H.M. Gupta, J.R. Campanha, Physica A 323 (2003) 626–634. D. Zheng, G.J. Rodgersa, P.M. Huic, Physica A 310 (2002) 480–486. G. De Fabritiis, F. Pammolli, M. Riccaboni, Physica A 324 (2003) 38–44. D. Sornette, R. Cont, J. Phys. I France 7 (1997) 431–444. H. Kesten, Acta Math. 131 (1973) 207–248. H. Takayasu, A.-H. Sato, M. Takayasu, Phys. Rev. Lett. 79 (1997) 966–969.