数学が難しくてあまりちゃんと読まれていないと思われるGabaixのランダム成長理論．成長率が都市規模に依存しない (Gibrat) とき，都市規模分布がZipfになるというメッセージはある程度知られていると思われるが，Mathematicaでシミュレーションしてみよう．
1. reflected GBM
都市規模S_{it}が幾何ブラウン運動に従う
$$dS_{it}/S_{it} = \mu dt + \sigma dB_{it},$$
が，外生的な下限S_{min}で反射するとする．下限がないと分布が拡散し続け，定常分布を持たなくなってしまう．
パラメーターを設定する．ここがいまいちだとZipf係数は-1とけっこう乖離する．
mu = 0.001(* 期待成長率. zipfになるのはmu=0 *)
sigma = Sqrt[0.1] (* SD. 少し大きめにして動きを激しくする. 下限にヒットして跳ね返す必要．実証はQJEの脚註17 *)
Smin = 1 (* 反射壁, 最小人口 *)
nCities = 1024 (* 都市数 *)
tMax = 4096 (* シミュレーション期間 *)
dt = 1 (* time step *)
ブラウン運動をMapで書く．
nextSize[s_List] :=
Map[Max[Smin, #] &,
s*Exp[(mu - sigma^2/2) dt + sigma Sqrt[dt]*
RandomVariate[NormalDistribution[0, 1], nCities]]]
あとは走らせるだけ．
initialSizes = ConstantArray[1.1, nCities] (* 初期の人口. Sminの近くに設定. *);
SeedRandom["gabaix"] (* mathematicaは文字列でも乱数のseedを設定できる *)
history = NestList[nextSize, initialSizes, tMax];
finalSizes = Last@history;
finalSizesNorm = finalSizes / Mean@finalSizes (* 基準化 *);
プロットは
sortedSizes = ReverseSort@finalSizesNorm;
ranks = Range@nCities;
ListLogLogPlot[Transpose@{sortedSizes, ranks - 1/2},
AxesLabel -> {"log Pop", "log Rank (-1/2)"},
PlotLabel -> "Zipf Plot at T=" <> ToString@tMax,
PlotStyle -> PointSize[Medium], PlotFit -> "Linear"
]
2. heterogeneous cities
命題2は，muやsigmaが異なっていても，それぞれGibrat則を満たしていれば集計したときZipf則が成り立つ，というものである．ここではsigmaが大きい群と小さい群の2グループに分けてみよう．
(* 2グループ *)
nCities1 = 500; nCities2 = 500;
sigma1 = 1; sigma2 = 1;
mu1 = 0.05; mu2 = 0.20;
(*シミュレーション関数 (幾何ブラウン運動+反映壁) *)
simulateZipf[nCities_, tMax_, sigma_, sMin_] :=
Module[{mu = -sigma^2/2, nextSize, initialSizes},
nextSize[s_List] :=
Map[Max[sMin, #] &,
s*Exp[mu + sigma*RandomVariate[NormalDistribution[0, 1], nCities]]];
initialSizes = ConstantArray[sMin, nCities];
Last[NestList[nextSize, initialSizes, tMax]]]
SeedRandom["eeckhout"] (* seed *)
(*グループ1 *)
group1 = simulateZipf[nCities1, tMax, mu1, sigma1];
(*グループ2*)
group2 = simulateZipf[nCities2, tMax, mu2, sigma2];
(*合計（集計）*)
aggregate = Join[group1, group2];
あとはプロット
zipfData[list_] := Transpose[{Range[Length@list] - 1/2, ReverseSort@list}]
ListLogLogPlot[{zipfData[group1], zipfData[group2],
zipfData[aggregate]},
PlotLegends -> {"Group 1 (sigma=0.05)", "Group 2 (sigma=0.20)",
"Aggregated"}, AxesLabel -> {"log Rank", "log Pop"},
PlotStyle -> {Automatic, Automatic, {Black, Thick}},
PlotLabel -> "Aggregation of Heterogeneous Cities"]
3. Kesten Processes
Zipf法則自体はもっとシンプルなKesten過程でもいける (Appendix 1)．
パラメーターの設定
(*gamma は期待値が1に近いがわずかに下回るのが定常性を得る上で重要 *)
gammasigma = 0.1
gammamu = - gammasigma^2/1.25
epsilonConst = 0.01; (*加法的ショック（最小規模を支える）*)
nCities = 1024 (* 都市数 *)
tMax = 4096 (* シミュレーション期間 *)
gammaDist = LogNormalDistribution[gammamu, gammasigma];
Mean@gammaDist (* 期待値0.997004 *)
遷移を定義する:
kestenStep[s_List] := RandomVariate[gammaDist, nCities]*s + epsilonConst
シミュレーションの実行
initialSizes = ConstantArray[1.1, nCities];
SeedRandom["pareto"]
history = NestList[kestenStep, initialSizes, tMax];
finalSizes = Last@history;
finalSizesNorm = finalSizes / Mean@finalSizes;
同様にプロットすると，
毎期epsilonを加えてなめらかにbottomが形成されるので崖はなくなるし，収束も早めだが，上の例はconcaveっぽい感じで係数-1.066ぐらい (SE=.022なので-1は棄却)．
Reference
Gabaix, Xavier. "Zipf's law for cities: an explanation." The Quarterly journal of Economics 114.3 (1999): 739-767.
