CSIS Discussion Paper No. 150

The Size Distribution of ‘Cities’ Delineated with a
Network Theory-based Method and Smartphone GPS
Data∗
Shota Fujishima†

Naoya Fujiwara‡

Yuki Akiyama§ Ryosuke Shibasaki§

Hodaka Kaneda¶

October 10, 2017

Abstract
We apply a network theory-based method to delineate ‘cities’ independent of administrative boundaries in Japan. We divide the country into cells of 1×1 km2 over whicha random walker moves. Its switching probabilities are specified by GPS-based people
flow data. We find the groups of the cells that describe the long-run behavior of therandom walker in the most eﬀective manner. The combination of two lognormal distributions yields a better fit of the city size distribution than the distribution with aPareto upper tail. A jump diﬀusion process is a stochastic process of the city population underlying such a distribution.
Keywords: City size distribution; Power law; Mixture of distributions; Community
detection; GPS data; Jump diﬀusion process
JEL Codes: R12; C46; C55∗

We are grateful to Hiroyasu Inoue, Takaaki Aoki, the participants at the European Meeting of the 7th
Urban Economics Association, and the seminar audience at the University of Tokyo and Kagawa University for their helpful comments and discussions. Ristu Sakuramachi helped us with processing the data
and running the computer programs. Shota Fujishima acknowledges financial support from the Kikawada
Foundation. Any remaining errors are our own.
†
School of Management, Tokyo University of Science, 1-11-2 Fujimi, Chiyoda, Tokyo 102-0071, Japan.
Phone/Fax: +81 3 3288 2563, Email: sfujishima@rs.tus.ac.jp‡
Institute of Industrial Science, the University of Tokyo.
§
Center for Spatial Information Science, the University of Tokyo.
¶
ZENRIN DataCom Co., Ltd.

1

1Introduction

Economic activities are often geographically localized, where the spatial extent of economic
activities is, for instance, measured by people’s commuting flows. However, a spatial extent
that is self-contained in terms of people’s commuting flows is generally not guaranteed tocoincide with preexisting administrative units such as city and county. This has motivatedresearchers and practitioners in regional science to delineate “metropolitan areas.”
The U.S. Oﬃce of Management and Budget defines metropolitan statistical areas (MSAs),
which we can obtain from the U.S. Census Bureau.1 MSAs are constructed by merging
counties that have strong economic and social ties. While MSAs are collections of legallybounded entities, Rozenfeld et al. (2011) consider the micro scale by considering cells of agiven size. By using gridded population data, they delineate metropolitan areas with theircity clustering algorithm in which they recursively grow the cluster by adding populated cells
within a prescribed distance of the cluster.
In general, algorithmic methods such as that mentioned above iteratively merge geographical entities according to some criteria. These criteria depend on free parameters that
researchers need to decide.2 For example, Rozenfeld et al. (2011) must specify the prescribed
distance, which decides the candidate neighboring cells to merge, to start their algorithm.
Moreover, algorithmic methods do not explicitly construct an objective function to be optimized, which makes the underlying model unclear. This is particularly relevant when wewant to relate the delineation of metropolitan areas to the decision problems of economic
agents.3
In this study, we invoke insights from community detection in network theory.4 A network
is generally a collection of nodes and links that connect nodes. A community is then defined
as a collection of nodes that have mutually strong relationships. Whether nodes have strong
ties is decided by the link structure. For example, we can readily think of an urban economy
as a network in which each geographical unit serves as a node and two nodes are linked if
there is a flow of people between them. Thus, the communities in such a network would
naturally correspond to metropolitan areas.
While various methods have been proposed in the literature on community detection, we
use the map equation method developed by Rosvall and Bergstrom (2008). In our model,
a random walker moves over geographical units. In particular, we consider the simplest
possible model in which the switching probability between two units is proportional to the1

See https://www.census.gov/programs-surveys/metro-micro.html.
In addition, the final result may be aﬀected by the initial condition of the algorithm.
3
We return to this point in the Conclusion.
4
See Fortunato (2010) and Fortunato and Hric (2016) for overviews of this literature.
2

2

volume of the flow of people between them. Motivated by information theory, the optimalcommunity structure minimizes the average code length necessary to describe the long-run
behavior of the random walker. How clustering works to reduce the average code length isin the same spirit as the real address. Specifically, grouping geographical units into clusters
allows us to assign the same name to multiple geographical units as long as they belong
to diﬀerent clusters. In fact, we see many Washington counties in the United States. This
reduces the average code length. At the same time, however, clustering incurs a cost because
we must assign names to the clusters. This tradeoﬀ gives rise to the optimal clustering, whereminimizing the “average” code length makes us assign a shorter code to those geographicalunits frequently visited by the random walker. Our task is formulated as a standard optimization problem because it follows that minimizing the average code length boils down tominimizing the entropy of the distribution of the random walker’s long-run visit frequencies
over nodes. As long as we accept the model, this method does not involve a free parameterthat aﬀects the resulting community structure. Moreover, we can readily capture the hierarchical structure of communities by making the entropy function hierarchical (Rosvall andBergstrom, 2011). This approach thus enables us to look at metropolitan areas at various
scales.5
We use GPS data on the flow of people in Japan taken from smartphone apps. Since
we have high-resolution location information for each trip, we are free to choose the scale ofgeographical unit. Extracting about half-a-million commuting trips across Japan, we divide
the country into cells of 1 km by 1 km and then identify the origin and destination cells foreach trip. This yields the volumes of commuting flow for each pair of cells, which are used
to compute the switching rates of the random walker.
Contrary to Rozenfeld et al. (2011) who use population data, we use only people flow
data, which describe the relationships among geographical units.6 If we rely on population
data that describe the characteristics of each geographical unit, we typically need additional
information. For example, it might become necessary to specify how the economic ties
between geographical units decay over distance. In fact, Farmer and Fotheringham (2011),
who use another popular method of community detection to delineate metropolitan areas in
Ireland, make such an assumption.7 On the contrary, people flow data already incorporate
5

To our knowledge, this is the first work that delineates metropolitan areas while capturing their hierarchical structure.
6
There is a large literature on detecting industrial agglomeration (see, e.g., Mori and Smith (2015) andreferences therein). In the context of industrial agglomeration, addressing the economic connections among
geographical units requires micro data on, for example, transaction volumes among establishments, which isdiﬃcult to obtain. In fact, agglomeration is usually detected with only data on the spatial distribution of
establishments.
7
See the Appendix for more details.

3the geographical structure because the volume of the flow of people tends to be small as the
distance becomes large.
By adopting detected communities, or metropolitan areas, we study their size distribution, which has received considerable research attention in the urban economics literature(Gabaix and Ioannides, 2004). To simplify the expositions, we hereafter call metropolitan
areas cities. Studying the city size distribution is particularly interesting here because ourcities are not constrained to follow legal boundaries. Moreover, we can examine the size distribution of cities of all sizes. We understand that the main point at issue in the literatureis whether the city size distribution is a single lognormal distribution (Eeckhout, 2004) orwhether its upper tail is Pareto (Giesen and Suedekum, 2014; Ioannides and Skouras, 2013;
Giesen et al., 2010). Our conclusion is aligned with neither of these distributions. In fact, wefind that a combination of two lognormal distributions is a better fit with the data than isa single lognormal distribution or a combination of lognormal and Pareto distributions. We
further show that the stochastic process of the city population behind such a distribution isapproximated by a jump diﬀusion process, which has a long history of application to financesince Merton (1976).
The rest of the paper proceeds as follows. In Section 2, we formally present our method.
In Section 3, we explain the data we use. In Section 4, we visualize the detected communities
and, in Section 5, we study their size distribution. The last section concludes.

2

The Model

We divide the whole of Japan into cells of approximately 1 km by 1 km and consider a modelin which a random walker moves over the cells. Let nij be the number of workers commuting
from cells i to j and mij be the number of workers returning home from cells i to j. The
total number of commuting trips from cells i to j is then Nij = nij + mij . The probabilityof a random walker switching from cells i to j, Pij , is given by
Pij =

Nij
,
Ni

(1)

∑
where Ni = Sj=1 Nij and S is the total number of cells. Because mij = nji , Nij = nij + nji
and thus Ni is the sum of the daytime and nighttime working populations of cell i.
We focus on the largest recurrent class of the Markov chain defined by (1). As we
see in Section 4, this makes up around 95% of all populated cells in our data. Then, we
are interested in the probabilities of a random walker staying in each cell in the long run,which we call the long-run visit frequencies. Because the vector of long-run visit frequencies4p = (p1 , ..., pS ) is an invariant distribution of the Markov chain defined by (1), it satisfies
pP = p,(2)where P = [Pij ]i,j and S is the total number of cells in the largest recurrent class. Because
we focus on the largest recurrent class, p uniquely exists.
We describe the long-run behavior of the random walker in our network by assigningbinary codes, which are enumerations of numbers taking the values of either 0 or 1 suchas ‘01’ and ‘0010’, to each state. Our objective is to code these as eﬀectively as possible.
To make this more precise, let ℓj be the length of the binary code assigned to state j. For
example, the length of the binary code ‘0010’ is 4. Then, we assign codes to states tominimize the following average code length:∑ℓj pj ,

(3)jwhere pj is the long-run probability of state j. The key idea here is to reduce the averagecode length by bringing cells together into several communities because this allows us toassign the same code to diﬀerent cells as long as they belong to diﬀerent communities. We
also observe this type of information saving in real addresses: we see the same street nameeverywhere. However, grouping cells into communities makes it necessary for us to assign
codes to each community. This tradeoﬀ yields the optimal community structure.
We code by dividing the long-run behavior of the random walker into those inside each
community and those across communities. Let {C k } be a community partition such that∪C k ∩ C ℓ = ∅ for any k ̸= ℓ and k C k = S where C k ⊆ {1, 2, ..., S} is the set of cells
that belong to community k. Observe that the number of communities is endogenously
determined here. We describe the long-run behavior of the random walker inside community
k by assigning binary codes to the states of visiting each cell in the community and thestate of exiting the community. Then, the average code length for describing the long-run
behavior of the random walker inside community k is∑qkpi
Lk = ℓke k+ℓ,
ik + qkp + qk
p
k

(4)

i∈C

where ℓi is the length of the binary code assigned to the state of visiting cell i, ℓke is thelength of the binary code assigned to the state of exiting community k, q k is the probability∑of exiting community k, and pk = i∈C k pi . From (1), the probability of exiting community

5

k is
qk =

∑ ∑

pi Pij =i∈C k j ∈C
/ k

∑ ∑

pi

i∈C k j ∈C
/ k

Nij
.
Ni

(5)

On the contrary, we describe the long-run behavior of the random walker across communities by assigning binary codes to the states of visiting each community. Note that because
our Markov chain is stationary,∑ ∑pi Pij =i∈C k j ∈C
/ k

∑ ∑

pj Pji .

(6)

i∈C k j ∈C
/ kHence, q k also represents the probability of visiting community k from another community.
Therefore, the average code length for describing the long-run behavior of the random walker
across communities is
K
∑
qk
L=
ℓk ,
(7)
q
k=1
where ℓk is the length of the binary code assigned to the state of visiting community k, K∑
kis the number of communities, and q = K
k=1 q .
We consider the weighted sum of L and {Lk }, where the weight of L is q whereas theweight of Lk is pk + q k . Therefore, what we try to minimize is given by
qL +

K
∑

(pk + q k )Lk .

(8)

k=1

It might seem a daunting task to find the community structure and coding that minimize
the average code length above. However, Shannon’s source coding theorem simplifies our task.
To state the theorem, we define that the entropy of a probability distribution (π1 , ..., πI ) isH(π1 , ..., πI ) =

I
∑
i=1

πi log

1
.
πi

(9)

Then, the theorem states that it is possible to make the average code length arbitrarily close
to the entropy of the underlying probability distribution.8
In our coding problem, the entropies of the probability distributions in Lk and L are
k
q1
qK
i
H( pkq+qk , { pkp+qk }i∈C k ) and H( q , ..., q ), respectively. Therefore, our task reduces to finding

8

This is a fundamental theorem in information theory. See, e.g., Theorem 5.4.2 of Cover and Thomas(2012).

6

the community partition that minimizes∗

1

K

L (C , ..., C ) = qH

( 1

K
q
, ..., qq
q

)

K
∑

+

(pk + q k )H

(

{

qk
pk +q k,pi
pk +q k

k=1

}

)

i∈C k

.

(10)

Rosvall and Bergstrom (2008) call the objective function above the map equation.
Moreover, this method can be extended to allow for the hierarchical structure of communities, as is done by Rosvall and Bergstrom (2011). Let C kℓ ⊆ C k be a subcommunityof community k, C kℓm ⊆ C kℓ be a subsubcommunity of subcommunity ℓ, and so forth. Our
objective function is then constructed recursively as( 1L∗ ({C k }k ) = qHkq
, ..., qq
q

)

∑

+

L∗ ({C kℓ }ℓ ),

(11)

k

where
L∗ ({C kℓ }ℓ ) = (pk + q k )H(qk,
pk +q k{pi
pk +q k

)

}
i∈C k

+

∑

L∗ ({C kℓm }m ),

(12)

ℓ

and, at the lowest level of the hierarchy,∗

L (Ckℓ···r) = (pkℓ···r+qkℓ···r()H

q kℓ···r,pkℓ···r +q kℓ···r{pipkℓ···r +q kℓ···r}i∈C kℓ···r

)
.

(13)

Note that the depth of the hierarchy is also a choice variable. We exploit the flexibility ofthis method, and detect the hierarchical structure of metropolitan areas.

3

Data

R (Congestion Statistics), a GPS dataset of the flow of people in
We use Konzatsu-Tokei ⃝R refers to the people flow
Japan provided by ZENRIN DataCom Co., Ltd. Konzatsu-Tokei ⃝data collected from the individual location data sent from mobile phones with the Auto-GPS
function, enabled under users’ consent, through the “docomo map navi” service provided by
NTT DOCOMO, INC (Japan’s primary mobile service provider). Those data are processedcollectively and statistically to conceal private information. The number of users, who live
across Japan, is around half-a-million. The information on people’s trips is so detailed thatthe time interval between the acquisitions of location information via GPS is five minutes at
its the shortest interval. Moreover, the data are panel data that record each individual’s trips
every day.9 The data we use were collected for one year from January 1, 2012 to December
R NTT Docomo, INC performed an overall and statistical processing ofIn generating Konzatsu Tokei ⃝,GPS data as per order of ZENRIN DataCom CO., LTD. This applies to all the figures presented in this97

31, 2012. Each data entry includes information such as the unique user ID, location (latitude
and longitude), and time stamp. The attributes of users such as age and sex are not available.
Although detailed information on people’s trips is available, the purpose of each trip is
not specified. Moreover, each user’s place of residence and work place, if he or she works,
are not known. Hence, to delineate metropolitan areas by using these data, we need to
extract commuting trips. To this end, we start by detecting the stops in each trip. Let us
represent the trip of a user by a sequence of locations {x1 , x2 , x3 , ...}, where xi is the i-th
location of the user. The location is composed of latitude y, longitude z, and time t. Hence,xi = (yi , zi , ti ). Then, given the threshold values S and T , we define the stop as a set of
locations {xk , xk+1 , xk+2 , ...xk+m } such that the distance between (yi , zi ) and (yj , zj ) is lessthan S for any i, j ∈ {k, k + 1, k + 2, ..., k + m} and tk+m − tk > T . That is, we regard a
set of locations as a stop if all the locations in the set are close to each other and the user
stayed there for a reasonably long time. We set S to 200 meters and T to five minutes.
Given these stops, we then identify the residential and working zones of users. Specifically,
the residential zone is identified as the zone to which the most frequent stop during the night(10pm-6am) belongs, whereas the working zone is identified as the zone to which the mostfrequent stop during daytime hours (9am-5pm) belongs.10
By following the procedure outlined above, we extract around 540,000 commuting trips.
To assess the reliability of our extracted data, we conduct two reliability checks. First, weaggregate the data to calculate the residential densities for each grid of cell size 1 km by
1km and compare those with the residential densities from the Grid Square Statistics of the2010 Population Census. As shown in Figure 1(a), we obtain a correlation coeﬃcient ofaround 0.9393. Second, we compute the trip volumes for each pair of municipal districts
and compare them with those from the 2011 Person Trip Survey conducted in the Chukyo
metropolitan area.11 As shown in Figure 1(b), we obtain a correlation coeﬃcient of around0.9369. These results show that our data are consistent with aggregated public data.

4

Detected Communities

By using the commuting data developed above, we first construct the Markov chain definedby (1) and find its largest recurrent class. Of the 85,607 cells, the largest recurrent class has
paper.
10
We refer to Horanont et al. (2013) for this method.
11
The data are similar to those derived from the National Household Travel
SurveyintheUnitedSates.
Thedata
arepublicly
available
onlineat
http://www.cbr.mlit.go.jp/kikaku/chukyo-pt/persontrip/p01.html (in Japanese).
The Chukyo
metropolitan area covers parts of the Aichi, Gifu, and Mie prefectures. The largest city in this metropolitan
area is Nagoya, which is located between Tokyo and Osaka.

8

0.040

4

0.035

3.0
Trip share from GPS data

Residential density from GPS data

3.5 ×10

2.5
2.0
1.5
The correlation
coefficient = 0.9393

1.0
0.5
0.0
0.0

0.030
0.025
0.020
0.015
The correlation
coefficient = 0.9369

0.010
0.005

0.5

1.5
1.0
2.0
2.5
Residential density from Census

3.0

0.000
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040
Trip share from PT survey

3.5
×10 4

「Konzatsu Tokei(r)」(c)ZENRIN DataCom CO., LTD.

「Konzatsu Tokei(r)」(c)ZENRIN DataCom CO., LTD.

(a) GPS data vs. Census(b) GPS data vs. Person Trip Survey

Figure 1: Reliability checks of GPS data
80,926 cells. We then apply the map equation method to our data limited to the largest
recurrent class.12 As a result, the cells are grouped into three level-1 (i.e., the highest level ofthe hierarchy) communities. As shown in Figure 2(a), the rightmost community is prominent:
it has major cities such as Tokyo, Osaka, and Nagoya and, moreover, the Hokkaido and thewhole of the Tohoku region are included in the community.
However, as the first level is too coarse to obtain insights for urban agglomerations, we
proceed to the lower level. Here, we detect 55 level-2 communities, which are depicted in
Figure 2(b). This level has an intuitively relevant scale for metropolitan areas. For example,
we can think of the rose community in the Kanto region as the Tokyo metropolitan area,
the salmon pink community in the Kansai region as the Osaka metropolitan area, the lime
community in the Tohoku region as the Sendai metropolitan area, and so forth. These
communities are divided into subcommunities (level-3 communities). The total number oflevel-3 communities is 2,048, although 30% are composed of fewer than 10 cells. Figure
3 depicts the level-3 communities in the Tokyo metropolitan area. We detect up to thesixth level for the hierarchy of communities, although the depth of the hierarchy is generally
diﬀerent among parent communities.

5

City Size Distribution

In this section, we study the size distribution of the detected communities in terms of pop∑
ulation, where the population of community k is computed as i∈C k Ni . To simplify the12

To carry out the map equation method, we use the code provided by Daniel Edler and Martin Rosvall
at http://www.mapequation.org/.

9

0

250

500

750

1000 km










6HQGDL0HWURSROLWDQ$UHD

2VDND0HWURSROLWDQ$UHD

7RN\R0HWURSROLWDQ$UHD「Konzatsu Tokei(r)」(c)ZENRIN DataCom CO., LTD.

「Konzatsu Tokei(r)」(c)ZENRIN DataCom CO., LTD.

(a) Level-1 communities(b) Level-2 communities

Figure 2: Maps of communities
0

25

50

75

100 km「Konzatsu Tokei(r)」(c)ZENRIN DataCom CO., LTD.

Figure 3: Map of level-3 communities in Tokyo area
expositions, communities, or metropolitan areas, are hereafter interchangeably called cities.
Studying the city size distribution, which has been a major research topic in the urban
economics literature, is particularly interesting here because our cities are independent ofadministrative boundaries and, by definition, the populations of all cities are available.
In the context of the city size distribution, Zipf ’s law has been regarded as an important
regularity condition for urban economics models (Gabaix and Ioannides, 2004). Let us review

10several concepts. A power law is a distribution function of the form Pr(X ≥ x) ∝ a/xη wherea, η > 0. Alternatively, we can say that X follows the Pareto distribution. Then, Zipf’s law
corresponds to the case where η = 1. When visualizing a power law, it is convenient torank the realized values of the variable (population in our context) by their sizes. Let x(i)be the i-th largest city’s population and K be the total number of cities. Under Zipf’s law,xη(i) Pr(X ≥ x(i) ) = a, where, as long as K is suﬃciently large, Pr(X ≥ x(i) ) may be regarded
as the share of cities larger than city i. Thus,
Pr(X ≥ x(i) ) ≈ i/K ⇒ log x(i) ≈ c − η −1 log i,(14)where c = η −1 log(aK). Therefore, under a power law, plotting log-population against logrank should yield the slope of −1/η as long as the sample size is suﬃciently large.13 See
Figure 4 for this type of graph plotted with our data.
Although Zipf’s law has been documented for the city size distributions in many countries,
Eeckhout (2004) points out that it appears to be important because only large cities are
considered. By using population data on “places,” which were newly introduced geographicalunits in the U.S. Census at that time, he claims that the city size distributions of all cities
are best fit by using a lognormal distribution. Observe that we also have the whole sampleof cities (i.e., detected communities).
Some authors, however, claim that a single lognormal distribution is insuﬃcient to describe the city size distribution and that the power law is important for, at least, the upper
tail of the distribution. It seems that this debate has not yet been settled.14 However, even
if the power law is relevant to the upper tail, it is still sound to consider the whole sample
because it is diﬃcult to truncate the sample in a convincing way.15 Then, a sensible suggestion to this debate would be to consider a distribution such that its body is characterized by
the lognormal distribution, whereas its tails are characterized by the Pareto distributions.
In fact, Giesen et al. (2010) and Giesen and Suedekum (2014) consider the double Pareto
lognormal distribution (DPLN) for U.S. city size distributions, concluding that the DPLN
has a better fit with the data than the lognormal distribution.16 This distribution, first
Note that when (14) holds with equality and η = 1, we have x(i) = 1i x(1) because log x(1) = c. Thus, thesize of the second largest city is half that of the largest city, the size of the third largest city is two-thirds
that of the second largest city, and so forth. Therefore, Zipf’s law is also called the rank-size rule.
14
See, for example, Levy (2009), Eeckhout (2009), Giesen et al. (2010), Bee et al. (2013), Ioannides andSkouras (2013), and Giesen and Suedekum (2014).
15
Eeckhout (2009) argues that truncating sample results in biased conclusions.
16
Ioannides and Skouras (2013) also consider a combination of lognormal and Pareto distributions. They
assume that there exists a switching population S̄ such that the distribution is lognormal for S ≤ S̄ and is
Pareto for S ≥ S̄, subject to the regularity condition that the density integrates to one. Their estimation
results also indicate that such a mixed distribution is better than the single lognormal distribution.
13

11

introduced by Reed (2002),17 is obtained as a normal variance mixture by the exponential
distribution (see the Appendix).
Thus, we also consider the lognormal distribution and DPLN. However, although the result that the combination of lognormal and Pareto distributions is better than the lognormal
distribution appears to be reasonable, it might be unfair to require a single lognormal distribution to describe a city size distribution for the whole sample. Therefore, we additionallyconsider a combination of two lognormal distributions. We present the stochastic process ofthe city population behind this type of distribution in the next section.
In sum, the density functions we fit to our data are given as follows. For convenience, we
consider the densities of log-population which is denoted by x:
(
)
1
x−µ
fN (x) = ϕ
,
σ
σ
(
){ (
)
(
)}
αβ
x−ν
x−ν
x−ν
fDP LN (x) =ϕm ατ −+ m βτ +
,
α+β
τ
τ
τ
(
)
(
)
x − µ11x − µ2
1+ (1 − θ) ϕ,fmixN (x) = θ ϕ
σ1
σ1
σ2
σ2

(15)
(16)
(17)

where ϕ is the standard normal density and m is the Mills ratio of the standard normallaw (i.e., m(x) = 1−Φ(x), where Φ is the cumulative standard normal distribution). fN
ϕ(x)corresponds to the case where the city population follows the lognormal distribution, fDP LN
corresponds to the case where the city population follows the DPLN, and fmixN corresponds
to the case where the city population follows the combination of two lognormal distributions.
See the Appendix for the derivation of the density of the DPLN.
We fit the three distributions above to our data by using the maximum likelihood estimation. Table 1 summarizes the estimation results.

17

See also Reed and Jorgensen (2004).

12

Table 1: Estimated parameters of the three distributions fitted to the level-3 communities
NormalNormal + Normal

µ1
σ1
µ2
σ2
θ
α
β
ν
τ

3.3472 (0.041)
1.8549 (0.029)

2.7910 (0.033)
1.0635 (0.024)
7.0222 (0.280)
1.7743 (0.171)
0.8685 (0.014)

Log-Likelihood
AIC
BIC-4,171
8,346
8,357

DPLN

0.5634 (0.013)
8.6315 (0.013)
1.6852 (0.006)
0.5843 (0.002)-3,7787,566
7,594-3,7997,606
7,628Standard errors are in parentheses.

These estimation results indicate that the DPLN is a better fit to our data than the lognormal
distribution in terms of any of likelihood, AIC, and BIC, which is in line with the findings ofGiesen et al. (2010) and Giesen and Suedekum (2014). However, the results also show thatthe combination of two lognormal distributions is better than the DPLN in terms of anyof likelihood, AIC, and BIC. This finding implies that although the lognormal distribution
is suﬃcient to describe the city size distribution, we must consider a combination of such
distributions. Hence, our results indicate that Pareto distributions are not relevant, even for
the tails.
We can also confirm this result graphically. Figure 4 plots the log-populations against
the log-ranks for the level-2 and level-3 communities. We can hardly see from Figure 4(b)that Zipf’s law, or more generally a power law, holds even for the upper tail. In the upper
left of Figure 4(b), the graph seems to be concave rather than a straight line; moreover,there is a jump between the log-populations of the third and fourth largest cities.18

18

On the contrary, the upper left of Figure 4(a) seems to be convex. This finding implies that the powerlaw does not hold for level-2 communities either.

13

12

10

10
log-population

log-population

12

8
6
4
2

8
6
4
2

0

2
4
log-rank

0

6

0

2

4
6
log-rank

8「Konzatsu Tokei(r)」(c)ZENRIN DataCom CO., LTD.

「Konzatsu Tokei(r)」(c)ZENRIN DataCom CO., LTD.

(a) Level-2 Communities(b) Level-3 Communities

Figure 4: The size distributions of communities in terms of population

5.1

Stochastic process behind the Gaussian mixture

One of the important motivations behind studying the city size distribution is to provide
regularity for the underlying theoretical model because the relevant model should yield thecity size distribution observed in the real world as an equilibrium. Although Eeckhout (2004)
and Giesen and Suedekum (2014) present the stochastic processes of the city population forthe lognormal distribution and DPLN, respectively, we are not aware of work that does so for
a combination of lognormal distributions. In this section, we show that adding a compound
Poisson process, where the jump amplitude is stochastic, to a geometric Brownian motion
yields the desired distribution. The resulting process has been commonly used in financesince the seminal work of Merton (1976).
∑ (t)Let Q(t) = N
i=1 Yi , where N is a Poisson process with rate λ and 1+Yi is iid lognormally
distributed with mean µ and δ 2 . Q is called a compound Poisson process. Then, we assume
that S(t) follows the following process:dS(t) = γS(t)dt + ξS(t)dW (t) + S(t−)dQ(t),(18)

where W is a Wiener process and S(t−) = lims↗t S(s). This is called a jump diﬀusion
process (see, e.g., Shreve, 2004, ch. 11). Given S(0) = s0 , the solution to (18) is given by

14

the Doleans-Dade exponential:{(
S(t) = s0 exp

)
} N∏
(t)
1 2
γ − ξ t + ξW (t)(1 + Yi ).
2
i=1

(19)

Note that because Yi is lognormally distributed, Yi > −1. This assures that S(t) does notjump to a negative value. Let X(t) = log S(t). Then,
(
)
N (t)
∑
1 2
X(t) = X(0) + γ − ξ t + ξW (t) +log(1 + Yi ).
2
i=1

(20)

By discretizing the time with interval ∆, we obtain
(

)
(
)
1 2
X(t) = X(t − ∆) + γ − ξ ∆ + ξ W (t) − W (t − ∆) +
2

N (t)
∑

(
)
log 1 + Yi .

(21)

i=N (t−∆)+1

We consider the probability transition density p(x; ∆|y, k) that satisfies19[]
Pr X(t) ∈ [x, x + dx] X(t − ∆) = y, N (t) − N (t − ∆) = k = p(x; ∆|y, k)dx.

(22)

By assumption, W (t) − W (t − ∆) ∼ N (0, ∆) and log(1 + Yi ) ∼ N (µ, δ 2 ). Then, from (21),p(x; ∆|y, k) = √(1ξ 2 ∆ + kδ 2

ϕ

)
)
(
x − y − γ − 12 ξ 2 ∆ − kµ√.
ξ 2 ∆ + kδ 2

(23)

Furthermore, because X(t − ∆) and N (t) − N (t − ∆) are independent, the variable representing the number of jumps is integrated out from the probability transition density as∞p(x, y; ∆) ∑ p(x, y, k; ∆) Pr[N (t) − N (t − ∆) = k]p(x; ∆|y) ==p(y; ∆)
p(y; ∆) Pr[N (t) − N (t − ∆) = k]
k=0

=
=

∞
∑
p(x, y, k; ∆) Pr[N (t) − N (t − ∆) = k]k=0∞∑p(y, k; ∆)p(x; ∆|y, k) Pr[N (t) − N (t − ∆) = k].

(24)
(25)
(26)

k=0

19

We focus on the transition density here because the stationary distribution does not exist, as in the caseof the geometric Brownian motion. However, it is possible to “stabilize” the process so that the stationary
distribution exists. One possible way is to allow for the possibility of people’s death. In such a model, Gabaix
et al. (2016) show in Proposition 8 that if the initial distribution and distribution of jump amplitudes bothhave Pareto tails, the stationary distribution also has a Pareto (upper) tail.

15

Because N is a Poisson process with rate λ,Pr[N (t) − N (t − ∆) = 0] = 1 − λ∆ + o(∆),(27)

Pr[N (t) − N (t − ∆) = 1] = λ∆ + o(∆),(28)

Pr[N (t) − N (t − ∆) ≥ 2] = o(∆).

(29)

Hence, we approximate N (t) − N (t − ∆) by the random variable Z that follows the Bernoulli
distribution with parameter λ∆ for a short time interval ∆. Then, we havep(x; ∆|y) ≈∑p(x; ∆|y, z) Pr(Z = z)(30)

= λ∆p(x; ∆|y, 1) + (1 − λ∆)p(x; ∆|y, 0).

(31)z∈{0,1}

Therefore,
(
(
)
)
1 2
2
2
X(t)|X(t − ∆) = y ∼ λ∆ × N y + γ − ξ ∆ + µ, ξ ∆ + δ
2
(
(
)
)
1 2
2
+ (1 − λ∆) × N y + γ − ξ ∆, ξ ∆ . (32)2Note that our estimation result is consistent with the regularity conditions behind (32).
First, for the approximation through the Bernoulli distribution to be accurate, ∆ must be
suﬃciently small. Thus, one of the two normal distributions must have a suﬃciently small
weight. This means that θ ≥ 12 must be suﬃciently large, and our estimate of θ is around0.87. Second, because µ, δ > 0, the normal distribution having the smaller weight must havethe larger mean and variance. Because the estimates of µ2 and σ2 are larger than those of µ1
and σ1 , respectively, this condition is also met. Therefore, we may conclude that the above
model is relevant to our data.20

Note that if we do not approximate N (t) − N (t − ∆) by using the Bernoulli distribution,
the transition density is an infinite mixture of the normal distributions according to the Poisk
son distribution.
In fact, ( because Pr[N (t) − N (t )− ∆) = k] = (λ∆)
e−λ∆ , p(x; ∆|y) =
k!
(
)
∞
∑
x − y − γ − 12 ξ 2 ∆ − kµ(λ∆)k −λ∆1√√
ϕ
. Yu (2007) proposes a method of approximatek!
ξ 2 ∆ + kδ 2
ξ 2 ∆ + kδ 2
k=0
ing the likelihood function based on this density.
20

16

6

Conclusion

By using the network theory-based map equation method following Rosvall and Bergstrom(2008) and smartphone GPS data, we delineated “cities” in Japan that are independent ofadministrative boundaries and then examined the size distribution of the delineated cities.
Contrary to previous observations in which the power law has been regarded as an important
regularity, we found that mixing lognormal distributions is suﬃcient to describe the city size
distribution and hence the Pareto distribution does not play a role. We argued that in sucha case, the jump diﬀusion process is relevant to the stochastic process of the city population.
Having said that, a fundamental problem in this literature is that how our cities look
like is decided by our definition of city. Hence, it is important to connect the definition ofcity to economics. One possible way of doing so might be to consider a coalition formation
game of local governments (Weese, 2015). One conceptual advantage of the map equation
method is that it is explicitly formulated as an optimization problem. Moreover, it is wellknown that minimizing the entropy, as in the map equation method, is related to finding theequilibrium of logit-type discrete choice models. Hence, there would be potential for doing
this in the map equation method.

Appendix
A1 Comparison with the modularity method
Several methods have been proposed for community detection as summarized in Fortunatoand Hric (2016). Among other things, Farmer and Fotheringham (2011) use the modularity
method of Newman and Girvan (2004) to delineate metropolitan areas. This method searches
for a community partition such that the volume of flows within a community is large, whereasthat of flows between communities is small. The performance of a community partition isevaluated relative to the case where cells are placed completely at random. Specifically, if
cells are placed at random, the expected number of commuting trips from cells i to j isN
Ni Nj
. The modularity is then defined as
N × NNi × Nj = N
)
(
1 ∑
Ni Nj
Q=
δij ,
Nij −
N i,j
N

(33)

where δij = 1 if cells i and j belong to the same community and δij = 0 otherwise. The
modularity method finds the community partition that maximizes the modularity, which does
Ni Nj
not depend on the free parameters. However, because N
only depends on the numbers of17

workers in each cell, it can take a large value even if the two cells are far apart as long as
they host large numbers of workers. Owing to this, Farmer and Fotheringham (2011), who
use commuting flow data in Ireland, needed to discount the flow volumes by geographicaldistance to obtain communities of reasonable sizes. However, this means that the resulting
community structure depends on the discounting method. The map equation method, on
the contrary, does not use information on the population of each cell and hence spatial
discounting is not necessary. In fact, the volume of the flow of people, which is the only
information used by the map equation method, naturally tends to be small as the distance
between origin and destination becomes large. Moreover, the map equation method canaccommodate the hierarchical structure of communities and hence detect various scales.
Hence, we do not have to explicitly take the distance into account.

A2 DPLN
As demonstrated by Eeckhout (2004), conditional on the city age and population in theprevious period, the city size distribution is given by a lognormal distribution if the logarithmof the city population follows a geometric Brownian motion. Specifically, let S(t) be the city
population at time t, and suppose that it obeys the following process:dS(t) = γS(t)dt + ξS(t)dW (t),(34)

where W is a Wiener process. Given S(0) = s0 , it follows from the standard Itô calculusthat the solution to (34) is
1 2
S(t) = s0 e(γ− 2 ξ )t+ξW (t) .
(35)
Let X(t) = log S(t). Then,()
1 2
X(t) = X(0) + γ − ξ t + ξW (t).
2(36)

Suppose that S(0) is lognormally distributed with mean ν and variance τ 2 (i.e., X(0) ∼N (ν, τ 2 )). Because W is a Wiener process (i.e., W (t) ∼ N (0, t)), we have
)
]
(
1 2
2
2
X(t) ∼ N ν + γ − ξ t, τ + ξ t .
2
[

(37)

The distribution above depends on t, which is the city age in our context. Giesen and
Suedekum (2014) assume that the city age is exponentially distributed with rate parameter18

λ. Hence, the density function of t is fT (t) = λe−λt . Then, the unconditional density of X is∫fX (x) =

1
ϕ
σt

(

x − µtσt)fT (t)dt,

(38)

(
)
where ϕ is the standard normal density, µt = ν + γ − 12 ξ 2 t, and σt2 = τ 2 + ξ 2 t. Letm(z) = 1−Φ(z)where Φ is the standard normal cumulative distribution. It is shown thatϕ(z)αβfX (x) =
ϕ
α+β

(

x−ν
τ

){

(

x−ν
m ατ −
τ)(x−ν
+ m βτ +
τ

)}
,

(39)

where α and −β (α, β > 0) are the roots of the following equation:
(
)
ξ2 2
1 2
z + γ − ξ z − λ = 0.
2
2

(40)

References
Bee, M., Riccaboni, M. and Schiavo, S. (2013). The size distribution of us cities: Not
pareto, even in the tail. Economics Letters, 120 (2), 232–237.
Cover, T. M. and Thomas, J. A. (2012). Elements of information theory. John Wiley &Sons.
Eeckhout, J. (2004). Gibrat’s law for (all) cities. The American Economic Review, 94 (5),1429–1451.
— (2009). Gibrat’s law for (all) cities: Reply. The American Economic Review, 99 (4),1676–1683.
Farmer, C. J. and Fotheringham, A. S. (2011). Network-based functional regions.
Environment and Planning A, 43 (11), 2723–2741.
Fortunato, S. (2010). Community detection in graphs. Physics reports, 486 (3), 75–174.
— and Hric, D. (2016). Community detection in networks: A user guide. Physics Reports,
659, 1–44.
Gabaix, X. and Ioannides, Y. M. (2004). The evolution of city size distributions. Handbook of regional and urban economics, 4, 2341–2378.
—, Lasry, J.-M., Lions, P.-L. and Moll, B. (2016). The dynamics of inequality. Econometrica, 84 (6), 2071–2111.
19

Giesen, K. and Suedekum, J. (2014). City age and city size. European Economic Review,71, 193–208.
—, Zimmermann, A. and Suedekum, J. (2010). The size distribution across all cities–
double pareto lognormal strikes. Journal of Urban Economics, 68 (2), 129–137.
Horanont, T., Phithakkitnukoon, S., Leong, T. W., Sekimoto, Y. andShibasaki, R. (2013). Weather eﬀects on the patterns of people’s everyday activities:
a study using gps traces of mobile phone users. PloS one, 8 (12), e81153.
Ioannides, Y. and Skouras, S. (2013). Us city size distribution: Robustly pareto, but
only in the tail. Journal of Urban Economics, 73 (1), 18–29.
Levy, M. (2009). Gibrat’s law for (all) cities: Comment. The American Economic Review,
99 (4), 1672–1675.
Merton, R. C. (1976). Option pricing when underlying stock returns are discontinuous.
Journal of financial economics, 3 (1-2), 125–144.
Mori, T. and Smith, T. E. (2015). On the spatial scale of industrial agglomerations.
Journal of Urban Economics, 89, 1–20.
Newman, M. E. and Girvan, M. (2004). Finding and evaluating community structure in
networks. Physical review E, 69 (2), 026113.
Reed, W. J. (2002). On the rank-size distribution for human settlements. Journal of Regional Science, 42 (1), 1–17.
— and Jorgensen, M. (2004). The double pareto-lognormal distribution - a new parametricmodel for size distributions. Communications in Statistics-Theory and Methods, 33 (8),1733–1753.
Rosvall, M. and Bergstrom, C. T. (2008). Maps of random walks on complex networks
reveal community structure. Proceedings of the National Academy of Sciences, 105 (4),1118–1123.
— and — (2011). Multilevel compression of random walks on networks reveals hierarchicalorganization in large integrated systems. PloS one, 6 (4), e18209.
Rozenfeld, H. D., Rybski, D., Gabaix, X. and Makse, H. A. (2011). The area andpopulation of cities: new insights from a diﬀerent perspective on cities. The American
Economic Review, 101 (5), 2205–2225.
20

Shreve, S. E. (2004). Stochastic calculus for finance II: Continuous-time models, vol. 11.
Springer Science & Business Media.
Weese, E. (2015). Political mergers as coalition formation: an analysis of the heisei municipal amalgamations. Quantitative Economics, 6 (2), 257–307.
Yu, J. (2007). Closed-form likelihood approximation and estimation of jump-diﬀusions with
an application to the realignment risk of the chinese yuan. Journal of Econometrics,
141 (2), 1245–1280.

21