■
ずっと気になっている(しかしマイナな)問題が解決しないので,そのメモを記録しておきます.
所得分布の生成モデルは,投資行動の反復から対数正規分布を導出する数理モデルである.仮定は次のとおりである(Hamada 2003; 2004; 2016).
- 人々が最初に初期資産
y_0
を持っている. - 資産のうち割合
b
を投資して,確率p
で成功すると,資産が増える.失敗すると投資した資産を失う. - 資産の一定割合を所得として得る.人々が投資行動を
n
回繰り返すと,資産と所得の分布は対数正規分布にしたがう.
このモデルの結論の一つは,導出された対数正規分布のパラメータがモデルの外生的パラメータであるy_0, b, p, n
によって内生的に決まる,というものである.
具体的には,資産の分布は
となる.証明はHamada (2003; 2004; 2016)をみよ.
これは2項分布が正規分布で近似できる,という性質を利用した命題である(確率変数を基準化した場合にはド・モアブル=ラプラスの定理となる).
さて
このモデルから得た命題をシミュレーションによって再現できない,というのが問題だった.
どういうことか説明しよう.
まず次の関数によって,確率p
で資産y_0 b
を増やし,1-p
で資産y_0 b
を失う仮定を表現する.コードは例によって一番楽に書けるMathematica.
gmi[y0_, b_, p_, n_] := Module[{y = y0}, Do[(* repeat n time *) If[Random[] < p, y = y + (y*b), y = y - (y*b)], {n}]; y];
一方で,理論的に導出される定まる確率分布のパラメータを使って,対数正規分布に従う疑似乱数をコンピュータに作らせる.
(* generate samples from a pdf with parameters derived from Hamada model *) lognsamples[y0_, b_, p_, n_, pop_] := Module[{m, s, a, B }, a = Log[(1 + b)/(1 - b)]; B = Log[y0] + n Log[1 - b]; m = B + a n p; s = Sqrt[n p (1 - p) a^2]; RandomVariate[LogNormalDistribution[m, s], {pop}]];
モデルが正しく対数正規分布を導出できているのであれば,ミクロレベルでの資産生成プロセスを集約した結果と,分布からサンプリングした結果は,ほぼ等しくなるはずである.
modelfit[y0_, b_, \[Gamma]_, n_, pop_] := Module[{dat0, dat1}, (* generate samples from a model *) dat0 = Table[gmi[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*)], {pop(* population *)}]; dat1 = lognsamples[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*), pop]; Histogram[{dat0, dat1}, Automatic, "Probability", ChartStyle -> {Red, Blue}]];
ところが,この関数を実行すると,「理論値から得たサンプルデータ」と「疑似乱数で作ったサンプルデータ」がズレる.
ズレの解消の仕方として,疑似乱数でサンプルデータを作るのではなく,確率密度関数を理論値で得た離散的所得階級値間のレンジで積分して,確率密度を少し「まとめて」やれば,うまくいくことは分かっている.
だが,いちいち積分させるのもコードが煩わしい.
先日,ふと気づいた.疑似乱数から生成したデータベクトルの中には同じ値がdouble精度で1つもないのに,理論に基づいて生成した方は,double精度で同じ値がいくつも含まれているのだ.ようするに,理論値のベクトルの方が疑似乱数のベクトルよりも分散が小さいのである.
これは一体,どういうことなのか?
理論値は理論値で間違っていない.疑似乱数ベクトルも指定したパラメータを持つ分布からランダムにサンプルを抽出している.
分布を導出しているはずの数理モデルから生成したサンプルが,なぜPCで発生させた確率分布疑似乱数と分布がずれてしまうのか?
そこで理論値は厳密すぎて,疑似乱数のランダムさを失っていると考えてみる.
確かに関数gmi[y0_, b_, p_, n_]
の中で一様乱数を使っているので投資結果はランダムに変化するのだが,獲得した資本額がある決まった値しかとらない.だから理論値にノイズを足してやればいいのでは? というわけだ.
別の言い方をすれば疑似乱数の方は連続確率変数の実現値を表現しているのに,理論モデルは離散的確率変数の実現値しか表現できていない.だから確率密度関数を区間毎に積分する操作が,実現値の離散化に対応していたのである.よって修正はただノイズ
RandomVariate[NormalDistribution[0, 0.1]]
を加えるだけでよいのでは,と考えた.
修正したコードはこうなる
lognsamples[y0_, b_, p_, n_, pop_] := Module[{m, s, a, B}, a = Log[(1 + b)/(1 - b)]; B = Log[y0] + n Log[1 - b]; m = B + a n p; s = Sqrt[n p (1 - p) a^2]; RandomVariate[LogNormalDistribution[m, s], {pop}]]; gmi[y0_, b_, p_, n_] := Module[{y = y0}, Do[(* repeat n time *) If[Random[] < p, y = y + (y*b), y = y - (y*b)], {n}]; y]; modelfit[y0_, b_, p_, n_, pop_] := Module[{dat0, dat1}, dat0 = Table[gmi[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*)] + RandomVariate[NormalDistribution[0, 1/Sqrt[n]]], {pop(* population *)}]; (* add randome error since the model is too rigorous *) dat1 = lognsamples[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*), pop]; Histogram[{dat0, dat1}, ChartStyle -> {Red, Blue}]];
これで,だいぶマシにはなる.だがパラメータの値によってノイズの分散に微妙に調整が必要なので 完璧ではない.
多分,記録しておかないと数日したら忘れちゃうので,書き留めておきます.