$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),
we instead
   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;