一般均衡モデルと不動点アルゴリズム２００３．１２．３０．関東学園大学犬童健良（M消費者、N産業）の場合 • 消費者の効用最大化 – 効用関数U（・）は厳密擬凸かつ微分可能とすると、解（需要関数）は価格に対し0次同次となる。 • 生産者の費用最小化 – 一次同時の生産関数Q から導かれる要素需要は要素価格に対し0次同次となる。   U m X 1m ,  , X Nm  max. s.t. m m m p X  wL  rK i 1 i i N where X im , Lm , K m , etc. are initial endowm ents.   Q j L j , K j  min. j  1,, N .    L j  L j r , w, Q j ,   K j  K j r , w, Q j . 均衡の条件 • • • 均衡の定義：全財・要素の超過需要関数が0以下となる価格ベクトル。利潤０条件（相補スラックネス条件）：産出量が正のとき、利潤は0である。 Walras法則：均衡条件のラグランジアンが非正。 Equilibrium : m j   X p ,  , p , r , w  Q  0, m1 j 1 N M   r, w, Q    j j m L r , w , Q  L  j 1 m1  0, N  j 1 K j N M j M m K 0 m 1 for j  1,, N and zero profit conditionsbelow. Zero profit conditions: Q j  0  p j Q j  wL j  rK j for j  1,, N . Walras law :        j 1 p j Q j  m1 X mj  w  j 1 Lj  m1 Lm  r  j 1 K j  m1 K m  0. N M N M N M 均衡の計算手順 1. 2. 3. 4. 5. 計算に先立って、各財（企業）ｊについてその生産量Qjで生産要素LjとKjを除した新しい変数を作る。（解の次元を下げる工夫）利潤0条件（相補性条件）から、財価格を要素価格の関数として表す。各個人の需要を求め、市場需要を集計する。各財の市場需要に見合う生産量を求める。各要素の派生需要を求め、その集計的超過需要を0と置いて連立方程式を解く。ただしWalras 法則によりその合計は0である。 1. Lj Kj j  l r , w,1, j  k j r , w,1 j Q Q 2. j  1,, N . p j r , w  wl j r , w,1  rk j r , w,1 j  1, , N . 3. X mj r , w  X mj  p1 r , w,, p N r , w, r , w j  1, , N , m  1, , M . 4. Q j r , w  X mj  p1 r , w, , p N r , w, r , w L j r , w  l j r , w,1Q j r , w j  1, , N , K j r , w  k j r , w,1Q j r , w 5. j  1, , N .  k r , w   j 1 K j r , w  m 1 K m , N M  l r , w   j 1 L j r , w  m 1 Lm . N M j  1, , N , 不動点アルゴリズム • 一般均衡やゲームの解、つまりナッシュ均衡の存在は、連続写像の不動点として証明されている（いわゆる存在定理）。 • しかし複雑なモデルでは手計算が事実上ムリであるので、具体的に均衡を求めるには、その不動点を近似する計算方法（アルゴリズム）が必要になる。 • 一般均衡モデルの文脈では、Scarf(1967,1973) の古典的アルゴリズムとそれを改良発展させたMerrill(1972)、Van der Laan & Talman(1979)などが知られる。 • これらの方法では隣接する分割単体を探索するため、三角形分割(Hansen-Kuhn-McKinnon)を利用する。 • 例えば、２生産要素ＫとＬの価格ベクトル（ｒ、ｗ）を正規化した基本単体（ｒ＋ｗ＝１）を３角形分割し、超過需要がすべて０に十分近くなる頂点が見つかるまで分割単体間を移動する。定義（単体） • 定義（ｍ次元単体; simplex）：一次独立な整数ベクトルｂj j=1,…,m+1をの頂点とする凸包。すなわちＳ＝｛X⊆Zm｜X＝∑ajbj, ∑aj=1,aj≧0}. • 定義（基本単体）：各頂点ｂjがいずれも単位ベクトルのとき。 • 定義（辺単体）：ｍ－１次元単体S’の頂点がｍ次元単体Sの部分集合のとき。 3角形分割とピボッティング 1. 基点 – 頂点ｂ1をｎ＋１次元基本単体の中から選ぶ。 2. 単体分割 – – – 各分割に属するｎ＋１個の頂点を定める。ただしベクトルφは１からｎまでの整数の順列とする。また各ｅjを次のように定義する。 3. ピボッティング – – – 隣接する単体へ移動する。取り除く頂点をｂjとする。追加する頂点は図のようである。ただしｂｎ＋１はｂ１と隣接する。ＬｅｍｋｅとＨｏｗｓｏｎの非巡回の議論が適用され、解への収束が保証される。 1.   b1  b11 / D,, bn11 / D  I n 1 where b1j  0, b1j  Z   0 j  1,, n  1,  j 1 b1j  D. n 1 2. j b j 1  b j  e where D j  1,, n φ  1 ,, n  is a perm utation of 1,, n, and e1  1,1,0,,0, e 2  0,1,1,0,,0,, e n  0,,0,1,1. ｂｊ＋１ｂｊ’ 3.  b  b j 1  b j 1  b j . j ｂｊｂｊ－１Ｓｃａｒｆアルゴリズム 1. 2. 3. 4. 5. 格子数Dを選び3角形分割を施す。初期単体を選ぶ。初期単体の全頂点をその規則(*)にしたがってラベル付けする。複数の頂点に同じラベルがある場合、以下の手順４を行う。全ての頂点に異なるラベルが付いている場合は手順５に飛ぶ。ピボッティングしてからラベルを計算しなおし、手順3へ戻る。ただし初期単体のときは、基本単体の頂点を除去する。初期単体以外のときは、同じラベルをもつ頂点のうち、直前の手順で導入された頂点を残す。完全にラベル付けされた単体から任意の点を選ぶ。超過需要ベクトルをすべて計算し、もし十分な近似であったならば停止する。さもなくば格子数Dを増やし、新たな初期単体を選んでから、手順2に戻る。 (*) ラベル付け規則： • 基本単体の境界上の頂点に、最初に０となる座標の番号 i （又は識別記号）をつける。 • それ以外の頂点に、超過需要が正となる最初の財や要素の番号ｊ（同上）をつける。アルゴリズムの実行過程の例 • • • • D=5、1次元基本単体１＝ｒ＋ｗとする。超過需要関数のグラフが下左の図のようであったとしよう。このとき、各頂点のラベルは下右の図のようである。また頂点(1/5,4/5) を基点とする単体からのピボッティングを図中に示した。（(0,1)側から出発したとする。） • 移動先の単体では完全ラベル付け済みなので、均衡点の近似が得られた。精度がこれで十分なら、停止する。 w ρl 1 D  5, φ  1, b1  1, 4 / D, and e1  1,1. Then b 2  b1  e1 D  1 5 , 4 5  1 5 ,  1 5 w ρｋ let  2 5 , 3 5. k Pivoting:  b1  b 2  b 2  b1  4 5 , 6 5  1 5 , 4 5  3 5 , 2 5. k k l l 0 l 1r 0 1 r Prologを用いた実装と実験 • まずPrologを用いてShoven and Whalleyのテキストに出ている２財２消費者モデルをモデリングし、テキストに示された均衡の数値例を検証する。 • 次に、同上モデルに対するScarfアルゴリズムを実装し、均衡を実際に近似してみよう。 Shoven＆Whalleyの２×２市場経済 • Prologプログラムgeq00.plとして作成。筆者のホームページで公開している。 (http://www.us.kanto-gakuen.ac.jp/indo) • 図に実行結果とパラメータのモデル例を示す。テキストの数値と若干異なったが、均衡条件は十分満たされている。 %均衡の数値例（Shoven & Whalley より） price_of_product(firm(m),1.399). price_of_product(firm(s),1.093). price_of_capital(1.373). price_of_labor(1.000). consumer(r). %rich consumer(p). %poor ?- production(A,[K,L|B],C). firm(m). %maker A = firm(m) K = 6.21451 L = 26.3591 B = [1.5, 0.6, 2.0, 1.373, 1.0] C = [24.9405, 24.9405, 5.04871e-029, true] ; firm(s). %sales A = firm(s) K = 18.7894 L = 33.6307 B = [2.0, 0.7, 0.5, 1.373, 1.0] C = [54.3762, 54.3762, 5.04871e-029, true] ; market_share(firm(s),consumer(p),0.7). V = value(34.8918=r:16.1048+p:18.787) ; No ?- utility(A,B,C). endowment(consumer(p),capital,0). F = firm(s) A = consumer(r) B = [11.5116, 16.6699] C = 27.8641 ; endowment(consumer(p),labor,60). ?- [geq00]. % geq00 compiled 0.01 sec, 0 bytes Yes ?- total_expendenture_for(F,Q,V). F = firm(m) Q = quantity(24.9405=r:11.5116+p:13.4289) Q = quantity(54.3762=r:16.6699+p:37.7063) V = value(59.4332=r:18.2202+p:41.213) ; No market_share(firm(m),consumer(r),0.5). market_share(firm(s),consumer(r),0.5). market_share(firm(m),consumer(p),0.3). demand_substitution_elasticity(consumer(r),1.5). demand_substitution_elasticity(consumer(p),0.75). endowment(consumer(r),capital,25). endowment(consumer(r),labor,0). scale_of_business(firm(m),1.5). scale_of_business(firm(s),2.0). A = consumer(p) B = [13.4289, 37.7063] C = 50.8946 ; No ?- inputs_weight(firm(m),0.6). inputs_weight(firm(s),0.7). supply_substitution_elasticity(firm(m),2.0). supply_substitution_elasticity(firm(s),0.5). プログラムについての説明 • • • • 述語utility/3は消費者の効用最大化と商品需要、述語production/3は生産者の費用最小化する要素需要をそれぞれ計算する（ともにCES型関数）。実行に先立って、前頁右上下図のようなモデルパラメータと均衡価格のデータを与える必要がある。 total_expendenture_for/3は消費者の効用最大化の結果に基づき、企業別の商品需要に集計する。均衡条件のチェックは、production/3の中でも行っている。すなわち、上記のパラメータの下で、まず各企業にとって費用最小化する生産量と要素の派生需要が計算される。 production/3の第３引数（図では変数Cでバインドされたリスト）によって、消費者の効用最大化の結果である商品需要が、最適生産量と一致するか否かをチェックしている。図のように要素の超過需要を検証するプログラムを作成し実行すると、わずかな不均衡（資本側に正の超過需要）が検出されたものの、ほぼ均衡が達成されていることが分かる。 ?- excess_demand(A,B,C). A = capital(0.00388917=25.0039-25, price(1.373)) B = labor(-0.0102005=59.9898-60, price(1.0)) C = disequilibrium No ?- excess_demand(capital(EDK=K1-K,price(R)),labor(EDL=L1L,price(W)),T):endowment(consumer(r),capital,Kr), endowment(consumer(p),capital,Kp), endowment(consumer(r),labor,Lr), endowment(consumer(p),labor,Lp), K is Kr + Kp, L is Lr + Lp, production(firm(m),[Km,Lm,_,_,_,R,W],_), production(firm(s),[Ks,Ls,_,_,_,R,W],_), K1 is Km + Ks, L1 is Lm + Ls, EDK is K1 - K, EDL is L1 - L, equilibrium_check([EDK,EDL],T). equilibrium_check([EDK,EDL],disequilibrium):member(X,[EDK,EDL]), X>0, !. equilibrium_check(_,equilibrium). 表計算によるチェック • • 数値例が若干文献と異なるので念のため表計算でShoven＆Whalleyの例題とそこで主張された近似均衡を検算した。下図のように、先述のPrologモデルの出力とほぼ同じ結果が得られている。練習問題：図はエクセルのシーと埋め込みになっている。図をダブルクリックするか、エクセルにコピーしてから、資本要素価格ｒの値を1.3735に変えて実験してみよ。またSolverを用いて均衡を見出せ。 Shaven-Whalley Economy 2004/1/4 需要サイド Demand by consumers 消費者1 utility of 1 生産者m 生産者s consumer1 11.12618447 1.42492 1.714703 消費者2 utility of 2 consumer2 -2.28614561 0.139881 0.334425 prices subst. elasticity μ１ μ２ market share α1m α1s 1.5 0.5 0.5 α2m α2s 0.75 0.3 0.7 endowments / 供給サイド firm m Supply 派生要素需要Lm demands for inputs 26.35908 Km 6.214506 scale param φm 1.5 input weight δm 0.6 elasticity σm 2 firm s Ls 33.63072 Ks 18.78938 φs 2 δs 0.7 σs 0.5 Demands for products budget deficit α1[m,s]*p[m,s]^(1-μ1) x1m x1s px1 b1=px0_1 0.422728 0.478255 11.5116375 16.66991687 34.325 34.325 0 α2[m,s]*p[m,s]^(1-μ2) x2m x2s px2 b2=px0_2 0.326269 0.715736 13.42886495 37.70632931 60 60 0 xm xs 価格ベクトル 24.94050245 54.37624617 r w pm ps 1.373 1 1.399 1.093 budget deficit K0_1 L0_1 x1m x1s px1 b1=px0_1 25 0 11.514 16.674 34.33277 34.325 0.007768 K0_2 L0_2 x2m x2s px2 b2=px0_2 0 60 13.428 37.704 59.99624 60 -0.00376 K0 L0 xm xs <-- data for verification 25 60 24.942 54.378 intermediate terms for production functions K L profits (to be zeros) 25.00389 3.080465596 0.020814303 pQj-rKj-wLj=0 L K firm m firm s 59.9898 0.997156461 0.015966464 0.0001672 0.004693 超過需要 production functions excess demands 生産関数m 生産関数s Warlas law (=<0) K 24.94050245 54.37624617 -0.004861 0.003889 excess demands L 超過需要m 超過需要s Target -0.0102 0.000 0.000 0.0001192 Scarfアルゴリズムの実装 • 基本モデルでは予め価格のデータセットを与えて均衡条件を検証するのみだった。 • そこでScarfの不動点アルゴリズムを組み込むことによって、その均衡価格を近似できるようにしたい。 • 例題は基本単体が1次元のため、実装は容易である（ただし財価格は固定される）。 • またラベル貼りとそれに必要な超過需要関数の計算には先の例題モデルで作成した述語を再利用する。均衡近似の実験（１） • • • このPrologによるプログラムでは文献の数値との比較のため、正規化した要素価格に対応する 1次元単体の三角形分割を擬似的に用いたが、直接述語allocation/3を用いて出発頂点を数え上げており、明示的なピボッティングはしていない。以下の図は均衡近似の実験結果である。格子数10および5000のときの近似値は Shaven&Whallyの本のそれと近いことが分かる。ただし、財の価格は、要素価格の正規化スケールによって、以下のように調整する必要がある。 price_of_product(firm(m),A):- A is 1.399/(1.3735+1). price_of_product(firm(s),B):- B is 1.093/(1.3735+1). price_of_capital(A):-A is 1.3735/(1.3735+1). % 初期値 price_of_labor(B):-B is 1/(1.3735+1). % 初期値 ?- listing(log_simple). log_simple(id(1), prices([10/10, 0/10]), excess_demands([-20.6774, 660272.0]), label(2)). % 2. grid(10). ?- test_simple. approximated solution is : id(6), prices([5/10, 5/10]), excess_demands([7.63971, -2.76429]), label(1) log_simple(id(2), prices([9/10, 1/10]), excess_demands([-16.9082, 11.6681]), label(2)). log_simple(id(3), prices([8/10, 2/10]), excess_demands([-13.2364, 6.54272]), label(2)). Yes log_simple(id(4), prices([7/10, 3/10]), excess_demands([-8.38792, 3.61095]), label(2)). log_simple(id(5), prices([6/10, 4/10]), excess_demands([-1.73013, 0.687809]), label(2)). log_simple(id(6), prices([5/10, 5/10]), excess_demands([7.63971, -2.76429]), label(1)). % 2. grid(5000). ?- test_simple. log_simple(id(7), prices([4/10, 6/10]), excess_demands([21.1433, -6.91299]), label(1)). log_simple(id(8), prices([3/10, 7/10]), excess_demands([41.187, -11.7324]), label(1)). log_simple(id(9), prices([2/10, 8/10]), excess_demands([72.3653, -16.9863]), label(1)). log_simple(id(10), prices([1/10, 9/10]), excess_demands([126.664, -22.1559]), label(1)). log_simple(id(11), prices([0/10, 10/10]), excess_demands([2.05073e+006, -28.6762]), label(1)). Yes approximated solution is : id(2108), prices([2893/5000, 2107/5000]), excess_demands([0.00852542, 0.00268448]), label(1) Yes プログラミングの改良：練習問題 • • • 前頁に示したように、三角形分割とピボッティングを使わなくとも、Prologで格子点を数え上げるのは易しい。練習として、(前出のプログラムを見ずに！）N人にD個の同じ品物を配分する方法を枚挙する述語をalloc/3として作成しなさい。三角形分割とピボッティングをプログラムし、上記のallocを用いた結果と比較しなさい。次元数１，２，３の各場合について、実験してみること。上記の２種類のプログラムを用いて、例題の均衡を近似するプログラムを作りなさい。両者の性能を比較しなさい。述語allocationを用いた縦型探索：配分（価格）を数え上げながら要素超過需要が十分0に近いところを見つけ、近似が悪ければ格子数Dを増やしてやり直す。（実験（１）で作成済みのプログラムを手直しする。） Scarfアルゴリズム：三角形分割を用い、要素価格を更新しながら、各頂点の超過需要関数とそのラベルを計算する。ただし、ピボットしたときは新しい頂点についてのみ計算すればよい。すべてのラベルが異なる分割単体（つまり完全ラベル付け済み）で、その超過需要が十分0に近い頂点が見つかるまでピボッティングを続ける。近似が悪ければ格子数Dを増やしてやり直す。ピボッティングの実装 % test run for pivot/3. ?- B=out(1),D=out(2),F=out(3), H=out(2),J=out(1), A=[[4,0,0],[3,0,1],[3,1,0]], pivot(A,B,C),pivot(C,D,E), pivot(E,F,G),pivot(G,H,I), pivot(I,J,K). % 1-Dimension simplicies % N-Dimensional simplicies pivot([V1,V2],out(1),[V3,V2]):- pivot(V,out(K),T):- vector_diff(V2,V1,D), length(V,N), N > 2, vector_sum(D,V2,V3), nth1(K,V,VJ), \+ (member(X,V3),X<0). pivot(ring,K/N,[K0,K1]), nth1(K1,V,V1), pivot([V1,V2],out(2),[V1,V3]):- B = out(1) D = out(2) F = out(3) H = out(2) J = out(1) A = [[4, 0, 0], [3, 0, 1], [3, 1, 0]] C = [[2, 1, 1], [3, 0, 1], [3, 1, 0]] E = [[2, 1, 1], [2, 2, 0], [3, 1, 0]] G = [[2, 1, 1], [2, 2, 0], [1, 2, 1]] I = [[2, 1, 1], [1, 1, 2], [1, 2, 1]] K = [[0, 2, 2], [1, 1, 2], [1, 2, 1]] ; No ?- nth1(K0,V,V0), vector_diff(V1,V2,D), vector_diff(V1,VJ,D), vector_sum(D,V1,V3), vector_sum(D,V0,V2), \+ (member(X,V3),X<0). \+ (member(X,V2),X<0), substitute(K/N,V,_->V2,T). pivot(ring,1/N,[N,2]):-!. pivot(ring,K/N,[K0,K1]):- K > 0, K < N, !, K0 is K -1, K1 is K + 1. pivot(ring,N/N,[K0,1]):- K0 is N - 1, !. substitute(K/N,V,VK->V2,T):- length(V,N), nth1(K,V,VK), findall(W, (nth1(J,V,X), (J \= K -> W=X; W=V2)), T). 三角形分割の実装 %------------------------------------------------% subdivision of unit simplex %------------------------------------------------simplical_subdivision(dim(N), pivot(K),perm(Pi),(1/D)*(S->T) ):N1 is N + 1, vertex_of_simplical_subdivision( dim(N), vertex(N1),perm(Pi),(1/D)*S ), pivot(S,out(K),T). initial_simplical_subdivision( dim(N),(1/D)*S ):grid(D), N1 is N + 1, allocation(N1,D,S). vertex_of_simplical_subdivision(dim(N), vertex(1),perm(Pi),(1/D)*[S] ):permutation(dim(N),Pi), grid(D), initial_simplical_subdivision( dim(N),(1/D)*S ). ?- simplical_subdivision( vertex_of_simplical_subdivision( dim(N), vertex(K), perm(Pi), (1/D)*[S|[S0|O]] ):grid(D), length(Q,N), N1 is N + 1, nth1(K0,Q,_), (K0 > N-> !, fail;true), K is K0 + 1, vertex_of_simplical_subdivision( dim(N), vertex(K0), perm(Pi),(1/D)*[S0|O] ), member(K0->KP,Pi), findall(W, ( nth1(J,S0,X), exchange( dim(N),row(KP),col(J),DK ), W is X + DK, W >= 0 ), S), length(S,N1). dim(2),pivot(K),perm(Pi),(1/D)*(S->T) ). K=1 Pi = [ (1->2), (2->1)] D = 10 S = [[10, 0, 0], [9, 1, 0], [9, 0, 1]] T = [[8, 1, 1], [9, 1, 0], [9, 0, 1]] ; K=1 Pi = [ (1->2), (2->1)] D = 10 S = [[9, 1, 0], [8, 2, 0], [8, 1, 1]] T = [[7, 2, 1], [8, 2, 0], [8, 1, 1]] ; K=2 Pi = [ (1->2), (2->1)] D = 10 S = [[9, 1, 0], [8, 2, 0], [8, 1, 1]] T = [[9, 1, 0], [9, 0, 1], [8, 1, 1]] Yes ?- 参考文献 • J.B. Shoven and J. Whalley (1992). Applying General Equilibrium. Cambridge University Press.[応用一般均衡分析―理論と実際. 小平裕 (訳). 東洋経済新報社.1993.] • H. E. Scarf (1982). The computation of equilibrium prices: An exposition. In K.J. Arrow and M.D. Intriligator (eds.), Handbook of Mathematical Economics, vol. II, North-Holand, Chapter 21.