$TITLE SAM-Drud-Kendrick
* Developed by Ruben Mercado

options limrow = 4;
option NLP=MINOS;
* Note to PC users-- From the proceeding line, remove the asterisk.
* i.e., on a PC, this code must be run with the option NLP=MINOS.

SETS
i general index /labor, capital, rural, urban, food, clothing/
s(i) sectors /food, clothing/
f(i) factors /labor, capital/
h(i) households /rural, urban/;
ALIAS (i,ip);
ALIAS (i,iq);

PARAMETERS
b(s) technical coefficients
a(i,ip) share coefficients;
b('food') = 1.2;  b('clothing') = 1;

TABLE sam(i,ip)
           labor capital rural urban food clothing
labor                                  75   85
capital                                50   60
rural        90     30
urban        70     80
food                      60    65
clothing                  60    85             ;

a(i,ip)= sam(i,ip) / sum(iq, sam(iq,ip));
DISPLAY a;

POSITIVE VARIABLES
p(i) price
q(i) quantiy
y(i) income
t(i,ip) payment
c(i,ip) commodity ;

VARIABLES
j  performance index;

EQUATIONS
eph(h)
eqs(s)
eys(s)
eyf(f)
eyh(h)
etfs(f,s)
ethf(h,f)
etsh(s,h)
eetsh(s,h)
ecfs(f,s)
eeys(s)
eeyf(f)
eeyh(h)
jd performance index definition;
* performance index equation
jd..    j =E= 0;

*sectors
eqs(s)..      q(s) =E= b(s)* prod(f, c(f,s)**a(f,s));
ecfs(f,s).. c(f,s) =E= a(f,s) * q(s) * p(s) / p(f);
eys(s)..      y(s) =E= p(s) * q(s);
etfs(f,s).. t(f,s) =E= p(f) * c(f,s);

*factors
eyf(f)..      y(f) =E= p(f) * q(f);
ethf(h,f)..  t(h,f)=E= a(h,f) * y(f);

*households
etsh(s,h)..       t(s,h) =E= a(s,h) * y(h);
eph(h)..             p(h)=E= prod(s, p(s)**a(s,h));
eetsh(s,h)..       t(s,h)=E= p(s) * c(s,h);
eyh(h)..            y(h) =E= p(h) * q(h);

*linkage
eeys(s)..             y(s) =E= sum(h,t(s,h));
eeyf(f)..             y(f) =E= sum(s,t(f,s));
eeyh('rural').. y('rural') =E= sum(f,t('rural',f));
*notice that we eliminate one of the linkage equations (Walras law)

*initial values to facilitate solver convergence
p.l(i) = 1;  q.l(i) = 1;  y.l(i) = 1;

*lower bound to avoid division by zero
p.lo(f) = 0.001;

*lower bounds to avoid undefined derivatives in exponential functions
p.lo(s) = 0.001;  c.lo(f,s) = 0.001;

*exogenous variables
q.fx('labor') = 2;  q.fx('capital') = 1;

*numeraire
p.fx('urban') = 1;

MODEL SAMDK /all/;
option iterlim = 10000;
SOLVE SAMDK MAXIMIZING J USING NLP;

PARAMETER REPORT;
REPORT(i, "price") = p.l(i);
REPORT(i, "quantity") = q.l(i);
REPORT(i, "income") = y.l(i);

DISPLAY REPORT;  DISPLAY t.l, c.l;