$title  Market Equilibrium with a Convex Piecewise Linear Supply Schedule

set     it              Individual points on the supply schedule /i1*i5/,
        i(it)           Points corresponding to constraints /i1*i4/;

table data(it,*)        Supply schedule (convex)

                q             p
        i1      0             0
        i2     10             0.3
        i3     200            0.5
        i4     500            0.8
        i5    1000            1;


parameter       a(i)    Intercept,
                b(i)    Slope;

a(i) = data(i,"p");
b(i(it)) = (data(it+1,"p")-data(it,"p"))/(data(it+1,"q")-data(it,"q"));
display a,b;

parameter       epsilon         Price elasticity of demand
                pbar            Reference price
                qbar            Reference quantity
                tau             Tax rate /0/;

nonnegative
variables       Q(it)   Quantity,
                P       Price;

equations       profit, market;

profit(i)..     a(i) + b(i)*Q(i) =g= P;

market..        sum(i, Q(i)) =g= qbar * ((P+tau)/pbar)**(-epsilon);

Q.UP(i(it)) = data(it+1,"q") - data(it,"q");

*       Fix demand at 350 to compute the benchmark price:

pbar = 1;
qbar = 350;
epsilon = 0;
model equilibrium /profit.Q, market.P/;
solve equilibrium using mcp;

*       Calibrate the demand function at this point and then compute
*       the results:

pbar = P.L;
epsilon = 0.5;
equilibrium.iterlim = 0; 
solve equilibrium using mcp;
abort$(equilibrium.objval > 1e-5) "Benchmark replication fails";
equilibrium.iterlim = 10000; 

set     tauval /0*100/;

parameter       pivotdata       Pivot data;

equilibrium.solvelink = 2;
loop(tauval,
        tau = tauval.val/100;
        solve equilibrium using mcp;
        pivotdata(tauval,"qbar") = qbar;
        pivotdata(tauval,"pbar") = pbar;
        pivotdata(tauval,"epsilon") = epsilon;
        pivotdata(tauval,"q") = sum(i,Q.L(i));
        pivotdata(tauval,"p") = P.L;
        pivotdata(tauval,"tau") = tau;
);
execute_unload 'pivotdata.gdx',pivotdata;
execute 'gdxxrw i=pivotdata.gdx o=piecewise.xlsm par=pivotdata rng=PivotData!a2 cdim=0';