emp11.gms : Test EMP formulations of scarfmcp

Description

```Scarf's Activity Analysis Example

Scarf, H, and Hansen, T, The Computation of Economic Equilibria.
Yale University Press, 1973.

Rather than form the MCP explicitly (as in the GAMSLIB model scarfmcp),
max_p  sum (h, i(h) * log(expend_h(p))) - p'*sum(h, endow(.,h))
s.t    A'p <= 0, p >= 0

Here expend_h is the expenditure function defined by:
expend_h(p)  = min_c  p'*c   s.t.  u_h(c) >= 1

This is detailed in Rutherford's 1992 paper entitled "Sequential Joint Maximization"
```

Small Model of Type : GAMS

Category : GAMS Test library

Main file : emp11.gms

``````\$Title Test EMP formulations of scarfmcp (EMP11,SEQ=500)

\$ontext
Scarf's Activity Analysis Example

Scarf, H, and Hansen, T, The Computation of Economic Equilibria.
Yale University Press, 1973.

Rather than form the MCP explicitly (as in the GAMSLIB model scarfmcp),
max_p  sum (h, i(h) * log(expend_h(p))) - p'*sum(h, endow(.,h))
s.t    A'p <= 0, p >= 0

Here expend_h is the expenditure function defined by:
expend_h(p)  = min_c  p'*c   s.t.  u_h(c) >= 1

This is detailed in Rutherford's 1992 paper entitled "Sequential Joint Maximization"

\$Offtext

sets
c       commodities

/capeop, nondurbl, durable, capbop, skillab, unsklab/

h       consumers

/agent1,agent2,agent3,agent4,agent5/

s       sectors

/d1, d2, n1, n2, n3, cd, c1, c2/;

alias (c,cc);

table e(c,h)  commodity endowments

agent1     agent2     agent3     agent4     agent5
capbop         3         0.1         2          1           6
skillab        5         0.1         6          0.1         0.1
unsklab        0.1       7           0.1        8           0.5
durable        1         2           1.5        1           2

table d(c,h) reference demands

agent1     agent2     agent3     agent4     agent5
capeop         4         0.4          2         5          3
skillab        0.2                    0.5
unsklab                  0.6                    0.2        0.2
nondurbl       2         4           2          5          4
durable        3.2       1           1.5        4.5        2

table data(*,c,s)  activity analysis matrix

d1          d2          n1          n2          n3

output.nondurbl                                 6.0         8.0         7.0
output.durable          4.0         3.5
output.capeop           4.0         4.0         1.6         1.6         1.6
input .capbop           5.3         5.0         2.0         2.0         2.0
input .skillab          2.0         1.0         2.0         4.0         1.0
input .unsklab          1.0         6.0         3.0         1.0         8.0

+          cd          c1          c2

output.capeop           0.9         7.0         8.0
input .capbop           1.0         4.0         5.0
input .skillab                      3.0         2.0
input .unsklab                      1.0         8.0;

parameter       alpha(c,h)      demand function share parameter;
alpha(c,h) = d(c,h) / sum(cc, d(cc,h));

parameter  a(c,s)  activity analysis matrix;
a(c,s) = data("output",c,s) - data("input",c,s);

parameters
pLev(c)  'optimal price level' /
capeop    1.11786500051246
nondurbl  0.546345872736342
durable   1.02036919464262
capbop    1.23055906867928
skillab   0.871845012358443
unsklab   0.287283691841206
/
;

positive variables
p(c)    commodity price,
y(s)    production,
i(h)    income;

equations
mkt(c)          commodity market,
profit(s)       zero profit,
income(h)       income index;

mkt(c)..        sum(s, a(c,s) * y(s)) + sum(h, e(c,h)) =g=

sum(h,
i(h) * alpha(c,h) / p(c));

profit(s)..     -sum(c, a(c,s) * p(c)) =g= 0;

income(h)..     i(h) =g= sum(c, p(c) * e(c,h));

model scarf / mkt.p, profit.y, income.i/;

* Now set up the equivalent model using emp
variables z;
equations objdef;

objdef..  z =e=   sum(h, i(h)*sum(c, alpha(c,h)*log(p(c))))
- sum(c, sum(h, e(c,h))*p(c));

model scarfemp /objdef, profit, income/;

p.lo(c)  = 0.00001\$(smax(h, alpha(c,h)) gt eps);

* fix the price of numeraire commodity:
i.fx(h)\$(ord(h) eq 1) = sum{c, e(c,h)};

* solve and check MCP model
p.l(c) = 1;
y.l(s) = 1;
i.l(h) = sum(c, p.l(c) * e(c,h));
solve scarf using mcp;
abort\$[scarf.solveStat <> 1] 'bad solveStat for MCP', scarf.solveStat;
abort\$[scarf.modelStat <> 1] 'bad modelStat for MCP', scarf.modelStat;
abort\$[smax{c, abs(pLev(c)-p.l(c))} > 1e-5] 'bad p level',
p.l, pLev;

* Set initial values for production (y) and other variables
profit.m(s) = 1;
p.l(c) = 1;
i.l(h) = sum(c, p.l(c) * e(c,h));

* Now reproduce the MCP solution with EMP, using the "equilibrium"  keyword
file myinfo / '%emp.info%' /;

put myinfo / 'equilibrium';
put / 'max z p objdef profit';
put / 'vi';
loop(h\$(i.lo(h) le i.up(h)),
put / income(h) i(h) );
putclose;

solve scarfemp using emp;
abort\$[scarfemp.solveStat <> 1] 'bad solveStat for EMP-equilibrium', scarfemp.solveStat;
abort\$[scarfemp.modelStat  > 2] 'bad modelStat for EMP-equilibrium', scarfemp.modelStat;
abort\$[smax{c, abs(pLev(c)-p.l(c))} > 1e-5] 'bad p level',
p.l, pLev;

y.l(s) = profit.m(s);
display y.l;

* Now reproduce the MCP solution with EMP,
* using "max z" in the solve statement and the "dualEqu" keyword
* Set initial values for production (y) and other variables
profit.m(s) = 1;
p.l(c) = 1;
i.l(h) = sum(c, p.l(c) * e(c,h));

put myinfo / 'dualEqu';
loop(h\$(i.lo(h) le i.up(h)),
put / income(h) i(h) );
putclose;

solve scarfemp using emp max z;
abort\$[scarfemp.solveStat <> 1] 'bad solveStat for EMP-dualEqu', scarfemp.solveStat;
abort\$[scarfemp.modelStat  > 2] 'bad modelStat for EMP-dualEqu', scarfemp.modelStat;
abort\$[smax{c, abs(pLev(c)-p.l(c))} > 1e-5] 'bad p level',
p.l, pLev;

y.l(s) = profit.m(s);
display y.l;
``````