title(A):- A=[ '% A Prolog Modeling of General Equilibrium :', '% Scarfs Fixed Point Algorithm in Shaven-Whalley Economy', '% created: 28--29 Dec 2003.(mere verification of equilibrium)', '% modified: 30 Dec 2003--7 Jan 2004. (fixed point algorithm)', '% modified: 15 Jul 2004. (corrected the positions of dynamic predicates for SWI-prolog 5.2.13.)', '% main predicates:', '?- test_simple.', '?- excess_demand(A,B,C).', '?- total_expendenture_for(F,Q,V).', '?- production(A,B,C).', '?- utility(A,B,C).' ]. % references: % J.B. Shoven and J. Whalley (1992). Applying General Equilibrium. % Cambridge University Press. % 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. %===================================== % Economy with 2 consumers and 2 firms % by Shoven and Whalley(1992). %===================================== % descriptive modeling. consumer(r). %rich consumer(p). %poor firm(m). %maker firm(s). %sales :- dynamic price_of_product/2. :- dynamic price_of_capital/1. :- dynamic price_of_labor/1. /* price_of_product(firm(m),1.399). price_of_product(firm(s),1.093). price_of_capital(1.373). price_of_labor(1.000). */ 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). %price_of_capital(0.5786). %price_of_labor(0.4214). :-dynamic market_share_at_time/4. market_share_at_time(0,A,B,C):- market_share(A,B,C). %--------------------------------------- % model parameters: demand side %--------------------------------------- market_share(firm(m),consumer(r),0.5). market_share(firm(s),consumer(r),0.5). market_share(firm(m),consumer(p),0.3). %market_share(firm(m),consumer(p),P):- % integer_between(0,10,K), % P is K/10. market_share(firm(s),consumer(p),Q):- market_share(firm(m),consumer(p),P), Q is 1 -P. %substitution elasticity in consumption. demand_substitution_elasticity(consumer(r),1.5). demand_substitution_elasticity(consumer(p),0.75). endowment(consumer(r),capital,25). endowment(consumer(p),capital,0). endowment(consumer(r),labor,0). endowment(consumer(p),labor,60). %--------------------- % utility function %--------------------- utility(Consumer,[Xm,Xs]/[R,W],U):- model_parameters_of(Consumer,[K,L,Am,As,E,Pm,Ps,CD]), price_of_capital(R), price_of_labor(W), demand_function_of(Consumer,[R,W],[K,L,Am,E,Pm,CD],Xm), demand_function_of(Consumer,[R,W],[K,L,As,E,Ps,CD],Xs), ces_utility_function([Xm,Xs,Am,As,E],U). model_parameters_of(Consumer,[K,L,Am,As,E,Pm,Ps,CD]):- endowment(Consumer,capital,K), endowment(Consumer,labor,L), market_share(firm(m),Consumer,Am), market_share(firm(s),Consumer,As), demand_substitution_elasticity(Consumer,E), price_of_product(firm(m),Pm), price_of_product(firm(s),Ps), common_denominator_for_demand_functions([Am,As,E,Pm,Ps],CD). common_denominator_for_demand_functions([Am,As,E,Pm,Ps],CD):- avoid_zero_for(Pm,Pm1), avoid_zero_for(Ps,Ps1), C1 is Am * Pm1^(1-E), C2 is As * Ps1^(1-E), CD is C1 + C2. demand_function_of(consumer(_),[R,W],[K,L,A,E,P,CD],X):- avoid_zero_for(CD,CD1), avoid_zero_for(P,P1), B is W*L + R*K, X is A * B/(P1^E * CD1). ces_utility_function([Xm,Xs,Am,As,E],U):- avoid_zero_for(Xm,Xm1), avoid_zero_for(Xs,Xs1), Ym is Am^(1/E) * Xm1^((E-1)/E), Ys is As^(1/E) * Xs1^((E-1)/E), U is (Ym+Ys)^(E/(E-1)). %----------------------------- % The model of demand %----------------------------- budget_of(Consumer,capital(endow(K),value(A)),labor(endow(L),value(W)),C):- endowment(Consumer,capital,K), endowment(Consumer,labor,L), price_of_capital(R), price_of_labor(W), A is R * K, B is W * L, C is A + B. expendenture_of(Consumer,[Xm,Xs],[Pm,Ps],Y=Ym+Ys):- utility(Consumer,[Xm,Xs],_U), price_of_product(firm(m),Pm), price_of_product(firm(s),Ps), Ym is Pm * Xm, Ys is Ps * Xs, Y is Ym + Ys. feasible_expendenture(Consumer,[Xm,Xs],[Pm,Ps],Y== Y. expendenture_by(firms(prices:[Pm,Ps]),quantity(X=m:Xm+s:Xs),value(Y=m:Ym+s:Ys)):- expendenture_of(consumer(r),[Xmr,Xsr],[Pm,Ps],_Yr=Ymr+Ysr), expendenture_of(consumer(p),[Xmp,Xsp],[Pm,Ps],_Yp=Ymp+Ysp), Xm is Xmr + Xmp, Xs is Xsr + Xsp, X is Xm + Xs, Ym is Ymr + Ymp, Ys is Ysr + Ysp, Y is Ym + Ys. expendenture_by(consumers,quantity(X=r:Xr+p:Xp),value(Y=r:Yr+p:Yp)):- expendenture_of(consumer(r),[Xmr,Xsr],[Pm,Ps],Yr=Ymr+Ysr), expendenture_of(consumer(p),[Xmp,Xsp],[Pm,Ps],Yp=Ymp+Ysp), Xr is Xmr + Xsr, Xp is Xmp + Xsp, X is Xr + Xp, Yr is Ymr + Ysr, Yp is Ymp + Ysp, Y is Yr + Yp. total_expendenture_for(firm(m), quantity(X=r:Xr+p:Xp), value(Y=r:Yr+p:Yp), prices(capital:R,labor:W) ):- utility(consumer(r),[Xr,_]/[R,W],_), utility(consumer(p),[Xp,_]/[R,W],_), price_of_product(firm(m),Pm), X is Xr + Xp, Yr is Pm * Xr, Yp is Pm * Xp, Y is Yr + Yp. total_expendenture_for(firm(s), quantity(X=r:Xr+p:Xp), value(Y=r:Yr+p:Yp), prices(capital:R,labor:W) ):- utility(consumer(r),[_,Xr]/[R,W],_), utility(consumer(p),[_,Xp]/[R,W],_), price_of_product(firm(s),Ps), X is Xr + Xp, Yr is Ps * Xr, Yp is Ps * Xp, Y is Yr + Yp. %--------------------------------------- % model parameters: supply side %--------------------------------------- scale_of_business(firm(m),1.5). scale_of_business(firm(s),2.0). inputs_weight(firm(m),0.6). inputs_weight(firm(s),0.7). %substitution elasticity in resources. supply_substitution_elasticity(firm(m),2.0). supply_substitution_elasticity(firm(s),0.5). %--------------------- % production function %--------------------- production(Firm,[Kz,Lz,S,D,E,R,W],[Q,Q1,SQE,T]):- total_expendenture_for(Firm,quantity(Q=_),value(_),prices(capital:R,labor:W)), firm_model_parameters(Firm,[S,D,E,R,W]), demand_function_of_capital([Q,S,D,E,R,W],Kz), demand_function_of_labor([Q,S,D,E,R,W],Lz), ces_production_function([S,D,E,Kz,Lz],Q1), SQE is (Q - Q1)^2, %nl,write((Firm, [S,D,E,R,W,Kz,Lz],SQE is (Q - Q1)^2)), (SQE < 0.1-> T = true; T= false). firm_model_parameters(Firm,[S,D,E,R,W]):- scale_of_business(Firm,S), inputs_weight(Firm,D), supply_substitution_elasticity(Firm,E), price_of_capital(R), price_of_labor(W). demand_function_of_capital([Q,S,D,E,R,W],K):- avoid_zero_divides([R,S,D,W],[R1,S1,D1,W1]), A is ((1-D1)*W1/(D1*R1))^(1-E), K is Q / S1 * (D1 * A + (1-D1))^(E/(1-E)). demand_function_of_labor([Q,S,D,E,R,W],L):- avoid_zero_divides([R,S,D,W],[R1,S1,D1,W1]), B is (D1*R1/((1-D1)*W1))^(1-E), L is Q / S1 * (D1 + (1-D1)*B)^(E/(1-E)). avoid_zero_divides([R,S,D,W],[R1,S1,D1,W1]):- avoid_zero_for(R,R1), avoid_zero_for(S,S1), avoid_zero_for(D,D1), avoid_zero_for(W,W1), (D=1 -> D1 is 1-10^(-10);D1=D). ces_production_function([S,D,E,Kz,Lz],Q):- avoid_zero_for(Kz,Kz1), avoid_zero_for(Lz,Lz1), LH is D * Lz1^((E-1)/E),%nl,write((z1,LH is D * Lz1^((E-1)/E))), KH is (1-D) * Kz1^((E-1)/E),%nl,write((z2,KH is (1-D) * Kz1^((E-1)/E))), Q is S * (LH+KH)^(E/(E-1)). %--------------------- % equilibrium condition %--------------------- excess_demand(capital(EDK=K1-K,price(R)),labor(EDL=L1-L,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). %-------------------- % test run %-------------------- /* ?- [geq00]. % geq00 compiled 0.01 sec, 0 bytes Yes ?- total_expendenture_for(F,Q,V,C). F = firm(m) Q = quantity(24.9447=r:11.5158+p:13.4289) V = value(34.8976=r:16.1106+p:18.787) C = prices(capital:1.3735, labor:1.0) ; F = firm(s) Q = quantity(54.3823=r:16.676+p:37.7063) V = value(59.4399=r:18.2269+p:41.213) C = prices(capital:1.3735, labor:1.0) ; No ?- production(A,[K,L|B],C). 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] ; 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] ; No ?- utility(A,B,C). A = consumer(r) B = [11.5158, 16.676]/[1.3735, 1.0] C = 27.8742 ; A = consumer(p) B = [13.4289, 37.7063]/[1.3735, 1.0] C = 50.8946 ; No ?- excess_demand(A,B,C). A = capital(0.00167711=25.0017-25, price(1.3735)) B = labor(0.00533747=60.0053-60, price(1.0)) C = disequilibrium ; No ?- */ %------------------------------------ % naiive fixed point approximation: % 1-dimesional simplecies case %------------------------------------ label_of_vertics(EDL,L):- nth1(L,EDL,X), X > 0,!. label_of_vertics(_,0). update_price(capital,S/D):- retractall(price_of_capital(_)), assert(price_of_capital(S/D)), !. update_price(labor,S/D):- retractall(price_of_labor(_)), assert(price_of_labor(S/D)), !. update_prices([Sk,Sl],D):- update_price(capital,Sk/D), update_price(labor,Sl/D), !. %-------------------------------- % a test routine for 1-D simplecies %-------------------------------- :- dynamic log_simple/4. :- dynamic log_simple_last_id/1. test_simple:- init_scarf_log, test_simple(_,_,_), fail. test_simple:- nl, %write(end), %listing(log_simple), log_simple(id(1),prices(_),excess_demands(_),label(L)), log_simple(id(N),prices([P,Q]),ED,label(L1)),L1\=L,!, grid(D), P = Sk/D, Q =Sl/D, update_prices([Sk,Sl],D), write('approximated solution is :'), nl,write((id(N),prices([P,Q]),ED,label(L1))). test_simple(prices([Sk/D,Sl/D]),excess_demands([EDk,EDl]),label(L)):- grid(D), allocation(2,D,[Sk,Sl]), update_prices([Sk,Sl],D), excess_demand(capital(EDk=_,_Pk),labor(EDl=_,_Pl),_T), label_of_vertics([EDk,EDl],L), update_scarf_log_id(_->N), assert(log_simple(id(N),prices([Sk/D,Sl/D]),excess_demands([EDk,EDl]),label(L))). log_simple_last_id(0). update_scarf_log_id(N0->N):- log_simple_last_id(N0), N is N0 + 1, retractall(log_simple_last_id(_)), assert(log_simple_last_id(N)). init_scarf_log:- abolish(log_simple/4), retractall(log_simple_last_id(_)), assert(log_simple_last_id(0)). %---------- % test run %---------- /* % 1. grid(10). ?- test_simple. endapproximated solution is : id(6), prices([5/10, 5/10]), excess_demands([7.63971, -2.76429]), label(1) Yes ?- listing(log_simple). :- dynamic log_simple/4. log_simple(id(1), prices([10/10, 0/10]), excess_demands([-20.6774, 660272.0]), label(2)). 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)). 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)). 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 ?- excess_demand(A,B,C). A = capital(7.63971=32.6397-25, price(5/10)) B = labor(-2.76429=57.2357-60, price(5/10)) C = disequilibrium ; No ?- % 2. grid(5000). ?- test_simple. endapproximated solution is : id(2108), prices([2893/5000, 2107/5000]), excess_demands([0.00852542, 0.00268448]), label(1) Yes ?- excess_demand(A,B,C). A = capital(0.00852542=25.0085-25, price(2893/5000)) B = labor(0.00268448=60.0027-60, price(2107/5000)) C = disequilibrium Yes */ %------------------------------------ % Scarf's fixed point approximation: % 1 or more-dimesional simplecies %------------------------------------ % stop criteria tolerance(0.01). :-dynamic dimension_of_simplex/1. dimension_of_simplex(4). :-dynamic grid/1. grid(10). exchange(dim(1),row(1),col(1),1). exchange(dim(1),row(1),col(2),-1). exchange(dim(2),row(1),col(1),1). exchange(dim(2),row(1),col(2),-1). exchange(dim(2),row(1),col(3),0). exchange(dim(2),row(2),col(1),0). exchange(dim(2),row(2),col(2),1). exchange(dim(2),row(2),col(3),-1). exchange(dim(N),row(J),col(K),D):- %dimension_of_simplex(N), N >2, integer(N), length(L,N),nth1(J,[_|L],_), J1 is J + 1, nth1(K,[_|L],_), (K = J -> (D = 1); true), (K = J1 -> (D = -1); true), (\+ member(K,[J,J1])-> (D = 0);true). permutation(dim(1),[1->1]). %permutation(dim(2),[1->1,2->2]). permutation(dim(2),[1->2,2->1]). permutation(dim(3),[1->1,2->2,3->3]). permutation(dim(3),[1->2,2->1,3->3]). permutation(dim(3),[1->1,2->3,3->2]). permutation(dim(3),[1->2,2->3,3->1]). permutation(dim(3),[1->3,2->2,3->1]). permutation(dim(3),[1->3,2->1,3->2]). permutation(dim(4),[1->X,2->Y,3->Z,4->W]):- member(X,[1,2,3,4]), member(Y,[1,2,3,4]),Y \=X, member(Z,[1,2,3,4]),Z \=Y, Z \=X, member(W,[1,2,3,4]),\+ member(W,[X,Y,Z]). %------------------------------------------------- % triangulation : subdivision of unit simplex %------------------------------------------------- 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]):- %dimension_of_simplex(N), permutation(dim(N),Pi), grid(D), initial_simplical_subdivision(dim(N),(1/D)*S). vertex_of_simplical_subdivision(dim(N),vertex(K),perm(Pi),(1/D)*[S|[S0|O]]):- %dimension_of_simplex(N), 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 %(W > 0 -> true;(nl,write(minus(W is X + DK)),fail)) ), S), length(S,N1). simplical_subdivision(dim(N),pivot(K),perm(Pi),(1/D)*(S->T)):- %dimension_of_simplex(N), N1 is N + 1, vertex_of_simplical_subdivision(dim(N),vertex(N1),perm(Pi),(1/D)*S), pivot(S,out(K),T). /* % test run. ?- simplical_subdivision(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 ?- simplical_subdivision(dim(3),pivot(K),perm(Pi),(1/D)*(S->T)),S=[[10,0,0,0],S2,S3,S4]. K = 1 Pi = [ (1->3), (2->2), (3->1)] D = 10 S = [[10, 0, 0, 0], [9, 1, 0, 0], [9, 0, 1, 0], [9, 0, 0, 1]] T = [[8, 1, 0, 1], [9, 1, 0, 0], [9, 0, 1, 0], [9, 0, 0, 1]] S2 = [9, 1, 0, 0] S3 = [9, 0, 1, 0] S4 = [9, 0, 0, 1] ; No ?- */ %------------------------------------------------- % pivotting vertex in the triangulation method %------------------------------------------------- pivot([V1,V2],out(1),[V3,V2]):- vector_diff(V2,V1,D), vector_sum(D,V2,V3), \+ (member(X,V3),X<0). pivot([V1,V2],out(2),[V1,V3]):- vector_diff(V1,V2,D), vector_sum(D,V1,V3), \+ (member(X,V3),X<0). pivot(V,out(K),T):- length(V,N), N > 2, nth1(K,V,VJ), pivot(ring,K/N,[K0,K1]), %nl,write(pivot(ring,K/N,[K0,K1])), nth1(K1,V,V1), nth1(K0,V,V0), vector_diff(V1,VJ,D), vector_sum(D,V0,V2), \+ (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, !. /* % 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). 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 ?- */ %================== % utility programs %================== avoid_zero_for(R,10^(-10)):- 0 is R, !. avoid_zero_for(R,R). integer_between(L,U,K):- L1 is integer(L), U1 is integer(U), N1 is U1 - L1 + 1, length(W,N1), nth1(K,W,_). vector_sum_2([S0k,S0l],[Dk,Dl],[Sk,Sl]):- Sk is S0k + Dk, Sl is S0l + Dl. vector_sum([],[],[]). vector_sum([S0|P],[D|Q],[S|R]):- vector_sum(P,Q,R), S is S0 + D. vector_diff([],[],[]). vector_diff([S0|P],[D|Q],[S|R]):- vector_diff(P,Q,R), S is S0 - D. % substitution % ----------------------------------------------------------- % 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). % sum % ----------------------------------------------------------- % sum([],0). sum([X|Members],Sum):- sum(Members,Sum1), Sum is Sum1 + X. % maximal solution for given goal clause : a naive solver %--------------------------------------------------------- max(X,Goal):- % X: the objective variable, % Goal: the objective function and constraints, setof((X,Goal),Goal,Z), member((X,Goal),Z), \+ ( member((Y,_),Z), Y > X ). % probability % ----------------------------------------------------------- % probability(W,P):- N1 is integer(1/W) + 1, length(L,N1),nth1(K,L,_), K1 is K - 1, P is K1/(N1 - 1). % allocation % ----------------------------------------------------------- % allocation(N,A,[X|Y]):- allocation(N,A,A,[X|Y]). allocation(0,_,0,[]). allocation(N,A,B,[X|Y]):- integer(A), length([X|Y],N), allocation(_N1,A,B1,Y), K is A - B1 + 1, length(L,K), nth0(X,L,X), B is B1 + X. % % probability (percentile) by using allocation % ----------------------------------------------------------- % probabilities(0,[]). probabilities(N,[X|Y]):- integer(N), length([X|Y],N), allocation(N,100,[X|Y]). :- title(A), forall(member(B,A),(nl,write(B))).