確率モデルの数値計算

確率モデルの数値計算

Numerical Simulations of Stochastic Models

物理学専攻 小野清敬

Kiyotaka Ono 1.1

はじめに

 物理のみならず様々な分野で用いられている確率 モデルを経済に適用する。経済を確率的に扱う理由 として、為替相場、株式市場の出来高などは刻々と 変わっており正確に推定、測定することができない。

よって、こうした変化はランダムな過程、つまり確 率過程として扱うほうが自然である。

^∗ [1]

経済や社会の個々の活動主体

(人・会社・国など)

をエージェントと呼ぶ。従来の経済学では主として 同一のエージェントの集まりがその対象とされてき た。しかし、青木

¹

の主張である、「経済システムは 異質的なクラスターの集まりで、平均的経済量を扱 うのではなくゆらぎを取り入れた確率的な扱いが必 要。」から、エージェントをいくつかの「タイプ」に 分類し、異なったタイプのエージェント間の相互作 用を調べ、いくつかの異なったタイプの群れが時間 と共にどういう経過をたどるか解析し、それによる 確率的ゆらぎを分析することを本研究の目的とする。

1.2

ダイナミクス

 将来における事象の起こる確率は現在の状態だけ から決まり、過去の状態には依存しない性質をマル コフ性といい、状態

i

から次の時点で状態

j

に移る 確率

P { X n+1 = j | X n = i, ..., X 1 = i 1 , X 0 = i 0 }

= P { X n+1 = j | X n = i } (1)

が成り立つ。

1

青木正直 カリフォルニア大学ロサンゼルス校

(UCLA)

経 済学部名誉教授

有限または可算の状態空間を持つモデルを連続マ ルコフ・チェーンとよび、その変化は、ある状態

j

から状態

k

への遷移として、時間が連続的に変化す る場合には、

P r(X t+τ = k | X t = j) = w j,k +o(τ), j ̸ = k, j, k ∈ S (2)

と規定できる

(w j,k :遷移率)。

1.3

マスター方程式

 マスター方程式は確率モデルの時間変動を捉える 方程式で、式

(3)

のように表される。これはある時 点である特定の状態をとる確率が、確率流の流入

(右

辺第

1

項)と流出

(右辺第 2

項)によってどのように 変化をするかを示している。

dP _j (t)

dt = ∑

k ̸ =j

[P _k (t)w _kj − w _jk P _j (t)], j ∈ S (3)

2.

モデル

 エージェントのタイプが二つあるクラスターモデル を考える。一つは技術力の優れた

advanced company

のクラスター１、もう一つは

less advanced com- pany

のクラスター２である。各クラスターのサイ ズを

n 1 , n 2

とすると、クラスターサイズの遷移率は

w { (n 1 , n 2 ), (n 1 + 1, n 2 ) } = c 1 n 1 + f 1 (4) w { (n 1 , n 2 ), (n 1 , n 2 + 1) } = c 2 n 2 (5) w { (n 1 , n 2 ), (n 1 − 1, n 2 ) } = d 1 n 1 (6) w { (n ₁ , n ₂ ), (n ₁ , n ₂ − 1) } = d ₂ n ₂ (7) w { (n 1 , n 2 ), (n 1 + 1, n 2 − 1) } = λf 1 d 2 (e + an 1 )n 2

(8)

w { (n ₁ , n ₂ ), (n ₁ − 1, n ₂ + 1) } = bλf ₁ d ₁ n ₁ (9)

となる。また、このモデルのマスター方程式は

∂P (n 1 , n 2 ; t)

∂t = P(n 1 + 1, n 2 ; t)d 1 (n 1 + 1) +P (n 1 , n 2 + 1; t)d 2 (n 2 + 1)

+P (n ₁ − 1, n ₂ ; t) { c ₁ (n ₁ − 1) + f ₁ } +P (n ₁ , n ₂ − 1; t)c ₂ (n ₂ − 1)

+P (n ₁ + 1, n ₂ − 1; t)(n ₁ + 1)(n ₁ + 1)b ₁

+P (n ₁ − 1, n ₂ + 1; t)(n ₂ + 1) { a ₂ (n ₁ − 1) + b ₂ }

− P (n 1 , n 2 ; t) { c 1 n 1 + f 1 + c 2 n 2 + d 1 n 1 + d 2 n 2

+n 1 b 1 + n 2 (a 2 n 1 + b 2 ) } (10)

となる。

図

1:

モデル

3.1

キュムラント法

 キュムラント母関数を用いて、マスター方程式を 解く。

確率母関数を

G(z 1 , z 2 ; t) = ∑

n ₁ ,n ₂

P(n 1 , n 2 ; t)z ⁿ ₁ ¹ z ₂ ^{n 2} (11)

と定義し、マスター方程式と確率母関数から確率母 関数

G

に対する偏微分方程式を導く。

また、キュムラント母関数を

K(θ 1 , θ 2 ; t) = lnG(e ^{θ 1} , e ^{θ 2} ) (12)

と定義する。

表

1:

キュムラント法でのパラメーター

  スタンダードパラメーター

γ = 0.1 r = 0, 0001 m = n = 0.01/r

b = 50/γ a = 0.0001 e = 0.001

次にキュムラント母関数をテイラー展開する。

K(θ 1 , θ 2 ; t) = − k 1 θ 1 − k 2 θ 2

+ 1

2 (k ₁₁ θ ² + 2k ₁₂ θ ₁ θ ₂ + k ₂₂ θ ₂ ² ) +・

・・

(13)

ここで

k i

はクラスター

i

のサイズ

n i

の平均、k

ii

は クラスターのサイズ

n i

の分散、k

12

は共分散である

(i = 1, 2)。

G

についての偏微分方程式を方程式

(13)

を使って、

常微分方程式

k ˙ 1 = 1 − (n + b)γ 1 k 1 + eγ 2 k 2 + aγ 2 A 0 , (14) k ˙ ₂ = − (m + e)γ ₂ k ₂ + bγ ₁ k ₁ − γ ₂ A ₀ , (15) k ˙ 11 =1 − 2(n + b)γ 1 k 11 + (2/γ + b − n)k 1

+eγ 2 (2k 12 + k 2 ) + 2aγ 2 (k 1 k 12 + k 2 k 11 ) + aγ 2 A 0 , (16) k ˙ 22 = − 2(m + e)γ 2 k 22 + γ 2 (2/γ + e − m)k 2

+bγ 1 (2k 12 + k 1 ) − 2aγ 2 (k 1 k 22 + k 2 k 12 ) + aγ 2 A 0 , (17)

k ˙ 12 = − ((n + b)γ 1 + (m + e)γ 2 )k 12 + bγ 1 (k 11 − k 1 ) + eγ 2 (k 22 − k 2 )

− aγ 2 (k 1 k 12 + k 2 k 11 − k 1 k 22 − k 2 k 12 ) − aγ 2 A 0 , (18)

が導かれ、この５つの方程式を解くことによってク ラスターの振る舞いを解析する。

式

(14)-(18)

の右辺をゼロとおいて解いた定常解は

k 1 = 252, k 2 = 499874, k 11 = 322, k 22 = 5

×

10 ⁷ , k 12 = 6494

となる。

表

2:

モンテカルロ法でのパラメーター パラメーター

c ₁ = 4.0 c ₂ = 7.2 d ₁ = 5.0 d ₂ = 7.1 f ₁ = 10.0 γ = 0.1 r = 0.0001 m = n = 0.01/γ b = 50/γ a = 0.0001 e = 0.001

3.2

シミュレーション

 初期値を

k 1 = 10000, k 2 = 50000, k 11 = 322, k ₂₂ = 5

×

10 ⁷ , k ₁₂ = 6494

にする。

結果：図

2〜図 6

初期値を定常解から大幅にずらしても、長期的に見 ると定常解近傍に解がいき、定常解の吸引域

(basin of attraction)

が大きいことがわかる。

4.1

モンテカルロ法

 モンテカルロ法を用いて遷移確率から直接クラス ターの振る舞いを見る。

4.2

シミュレーション

 初期値

n ₁ = 250, n ₂ = 500000

結果:図

7

モンテカルロ法において定常解の吸引域は確認で きなかった。キュムラント法との結果の違いとして、

キュムラント法はキュムラント母関数をテーラー展 開し、第

2

項までを扱って解いたのに対して、モン テカルロ法は遷移率から直接解いたためだと考えら れる。

5.

変動係数と共分散

 変動係数から各クラスターのゆらぎの大きさ、共 分散からセクター１,２の相関関係をみる。クラス ター

i(i = 1, 2)

の変動係数

r i

をキュムラント法にお いて

r i = √

k ii /k i , c o = k 12 /k 1 k 2 (19)

モンテカルロ法においては

r i =

√ < (n _i − < n _i >) ² >

< n i > , (20)

表

3:

変動係数・共分散

r ₁ r ₂ c _o

0.07 0.014 5.1 × 10 ^{− 5}

c o = < (n 1 − < n 1 >)(n 2 − < n 2 >) >

< n 1 >< n 2 > (21)

と定義する。

これより、キュムラント法における変動係数・共

分散は表

3、モンテカルロ法における変動係数・共

分散は図

8、図 9

となる。

6.

まとめ

 変動係数は平均値に対する相対的な散らばりの大

きさ

(ゆらぎ)

を表しており、一般の巨視的体系では、

非常に小さな量である。このモデルにおいて変動係 数が大きい場合、平均値だけではこの系の挙動をう まく説明することができない。これはクラスターサ イズの時間発展が過去の履歴に強く依存しているか らである。これを

non-self averaging phenomena

と 定義する。一方、変動係数が小さい場合、平均値に よってこの系の様子が分かる。これを

self averaging phenomena

と定義する。

^∗ [4]

本研究において、モンテカルロ法、キュムラント 法どちらにおいてもクラスター１は

non-self aver- aging phenomena、クラスター２は self averaging phenomena

という結果を得た。

経済におけるこの現象の例を挙げると、国民所得 の平均に基づいた経済政策を評価する場合、国民所 得の変動係数が大きいと平均値に基づいた経済政策 にそれほど信頼性がなく、経済政策の効果が単純に 予見できない。

参考文献

[1]

青木正直『異質的エージェントの確率動学入門』

共立出版

2003

[2]

青木正直『経済における確率的モデルへの招待』

サイエンス社

2004

[3] M.Aoki,T.Nakano,G,Yoshida,An Example of Schumperterian Dynamics,unpublished pa- per,2004

[4] M.Aoki,Thermodynamic Limits of Macroeco- nomic or Financial Models,EconPapers,2006 [5]

遠藤靖 『確率モデルの基礎』東京電機大学出版

局

2002

[6]

小倉久直『確率過程入門』森北出版株式会社

1998

[7]

阿部龍蔵『統計力学』東京大学出版会

1966

10000

8000

6000

4000

2000

0

k1

1000 800 600 400 200 0

tx10

³

k

₁

:252.424

図

2: k ₁

500

400

300

200

100

0

k2x103

1000 800 600 400 200 0

tx10

³

k

₂

:499873

図

3: k ₂

4000

3000

2000

1000

0

k11

1000 800 600 400 200 0

tx10

³

k

₁₁

:322.591

図

4: k ₁₁

55 50

45 40 35 30 25

k22x106

1000 800 600 400 200 0

tx10

³

k

₂₂

:4.999 x10

⁷

図

5: k 22

80

60

40

20

0

k12x103

1000 800 600 400 200 0

tx10

³

k

₁₂

:6494.63

図

6: k 12

503

502

501

500

499

x103

3000 2500 2000 1500 1000 500 blue : 10000 green : 50000 yellow:100000 pink :150000 purple:200000 black :250000

gray :300000 light blue:350000 brown :400000 navy :450000 red :500000

図

7: n ₁ , n ₂

100

80

60

40

20

0

r1,r2x10-3

500 400 300 200 100

tx10

³

blue :r

₁

green:r

₂

図

8:

変動係数

-20 -15 -10 -5 0

cox10-6

500 400 300 200 100

tx10

³

図

9:

共分散