Gabaix 1999 QJEのMathematica

 数学が難しくてあまりちゃんと読まれていないと思われる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"

] 

無理矢理反射させているので，bottomで崖のようになっていることが確認できる．

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)．

$$S_{t+1} = \gamma_{t+1} S_{t} + \epsilon_{t+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.