$title  Putty-SemiPutty Technology

*       Technology choice is based on an ex-ante isoquant which forms an envelope
*       for ex-post substitution possibilities.  That is, once a point on the ex-ante 
*       isoquant has been established, input-choices can be reoptimized in subsequent periods.

*       In this GAMS program, I evaluate the magnitude of approximation errors
*       which are introduced when installed technology is selected on the
*       basis of weighted or unweighted average factor prices over the life of the capital stock.

set     t       Time horizon /0*100/;


variables       
        obj     Variable profit -- present value return to a unit of capital
        Y(t)    Output
        L(t)    Employment
        K(t)    Capital stock 
        E(t)    Energy use
        Cbar    Reference cost
        Ebar    Reference energy use
        Lbar    Reference labor use
        Ybar    Reference output
        pKbar   Reference price of capital
        pLbar   Reference price of labor
        ALPHA   Reference energy value share;
                
*       Benchmark statistics:

parameter
        l0      Baseline labor demand /0.45/
        k0      Baseline capital stock
        e0      Baseline energy price /0.05/
        lamda   Capital survival share /0.93/
        r0      Interest rate /0.05/
        sigmaL  Long-run elasticity (dual) /0.5/
        sigmaS  Short-run energy elasticity (dual) /0.2/
        sigmaKL Short-run KL elasticity (dual) /0.5/

        Kbar    Reference capital use
        theta   Baseline energy value share
        rk0     Base year capital price
        rhoL    Long-run elasticity (primal)
        rhoS    Short-run energy elasticity (primal)
        rhoKL   Short-run elasticity (primal)
        srvshr  Survival share
        pref    Reference price
        kvs     Baseline capital value share
        py(t)   Projected output price
        pl(t)   Projected wage rate
        pe(t)   Projected energy price
        pinv    Price of new vintage capital;

*       Assign benchmark parameter values:

rhoL = 1 - 1/sigmaL;
rhoS = 1 - 1/sigmaS;
rhoKL = 1 - 1/sigmaKL;

rk0 = r0+1-lamda;
k0 = (1-e0-l0)/rk0;
kvs = rk0*k0/(rk0*k0+l0);

theta = e0;
pref(t) = 1/(1+r0)**ord(t);
py(t) = pref(t);
pl(t) = pref(t);
pe(t) = pref(t);
pinv = 1;
srvshr(t) = lamda**(ord(t)-1);

*       Fix the capital investment:

Kbar    = k0;


equations       objdef,output,capital,benchmark,alphadef,cbardef,lbardef,kbardef,ebardef;

*       Objective function: maximize present value of returns:

objdef..        obj =e= sum(t, py(t)*Y(t) - pl(t)*L(t) - pe(t)*E(t))/Kbar;

*       Period by period output is on the ex-post isoquant which is anchored at a reference point
*       on the ex-ante isoquant:

output(t)..     Y(t) =e= Ybar * ( ALPHA * (E(t)/Ebar)**rhos + 
                        (1-ALPHA) * (kvs * (K(t)/kbar)**rhokl + (1-kvs) * (L(t)/Lbar)**rhokl)**(rhos/rhokl))**(1/rhos);

*       Capital depreciates over time:

capital(t)..    K(t) =e= Kbar * srvshr(t);

*       Reference output level is on the ex-ante isoquant:

benchmark..     (theta * (Ebar/e0)**rhoL + (1-theta) * ((Kbar/k0)**kvs * (Lbar/L0)**(1-kvs))**rhoL)**(1/rhoL) =e= Ybar;

*       Benchmark value share, scaling the price of energy to unity:

alphadef..      ALPHA * (Ebar + pLbar*Lbar + rk0*pKbar*Kbar) =e= Ebar;

*       Benchmark cost (normalizing the price of energy to be unity):

cbardef..       Cbar =e= (theta + (1-theta)*(pKbar**kvs * pLbar**(1-kvs))**(1-sigmaL))**(1/(1-sigmaL));

Lbardef..       Lbar =e= Ybar * l0 * Cbar**sigmaL * (pKbar**kvs * pLbar**(1-kvs))**(1-sigmaL) / pLbar;

kbardef..       Kbar =e= Ybar * k0 * Cbar**sigmaL * (pKbar**kvs * pLbar**(1-kvs))**(1-sigmaL) / pKbar;

ebardef..       Ebar =e= Ybar * e0 * Cbar**sigmaL;

model puttyputty/all/;

*       Assign some initial values:

Y.L(t)  = srvshr(t);
L.L(t)  = l0 * srvshr(t);
K.L(t)  = k0 * srvshr(t);
E.L(t)  = e0 * srvshr(t);
Cbar.L  = 1;
Ybar.L  = 1;
Ebar.L  = e0;
Lbar.L  = L0;
ALPHA.L = theta;
pLbar.L = 1;
pKbar.L = 1;
obj.L = sum(t, py(t)*Y.L(t) - pl(t)*L.L(t) - pe(t)*E.L(t))/Kbar;

*       Install some bounds to avoid bad function calls:

Y.LO(t) = 0.01 * srvshr(t);
L.LO(t) = 0.01 * l0 * srvshr(t);
K.LO(t) = 0.01 * k0 * srvshr(t);
E.LO(t) = 0.01 * e0 * srvshr(t);
Cbar.LO = 0.01 * 1;
Ebar.LO = 0.01 * e0;
Lbar.LO = 0.01 * l0;
ALPHA.LO = 0.01 * theta;
pLbar.LO = 0.01 * 1;
pKbar.LO = 0.01 * 1;

solve puttyputty using nlp maximizing obj;

parameter       intensity       Energy intensity of output
                profit          Operating margin;

intensity(t,"Ref") = E.L(t)/(e0*Y.L(t));
profit(t,"Ref") = (py(t)*Y.L(t) - pl(t)*L.L(t) - pe(t)*E.L(t))/(srvshr(t)*pref(t));

*       Assume that energy prices on average increase at 3% per year up to a 
*       factor of 5 above the steady-state level:

loop(t,
        pe(t+1) = min(5*pe(t+1),pe(t)*1.03/(1+r0));
);

set tlbl(t) /0,20,40,60,80,100/;

$setglobal batch yes
$setglobal gp_opt3 "set key outside"
$setglobal gp_opt0 "set xlabel 'year'"
$if "%batch%"=="yes" $setglobal gp_opt1 "set term png"
$setglobal domain t
$setglobal labels tlbl
pe(t) = pe(t)/pref(t);
$if "%batch%"=="yes" $setglobal gp_opt2 "set output 'fig1.png'"
$libinclude plot pe

pe(t) = pe(t)*pref(t);

solve puttyputty using nlp maximizing obj;

profit(t,"Dynamic") = (py(t)*Y.L(t) - pl(t)*L.L(t) - pe(t)*E.L(t))/(srvshr(t)*pref(t));
intensity(t,"Dynamic") = E.L(t)/(e0*Y.L(t));
intensity(t,"Dyn_Ref") = Ebar.L/(e0*Ybar.L);

*       Now do a calculation in which we use a weighted average of energy prices
*       to select a point on the ex-ante isoquant:

pKbar.FX = sum(t,srvshr(t))/sum(t, srvshr(t)*pe(t)/pref(t));
pLbar.FX = sum(t,srvshr(t))/sum(t, srvshr(t)*pe(t)/pref(t));

*       After having fixed the point, we then compute the time path of response
*       to market prices:

solve puttyputty using nlp maximizing obj;
intensity(t,"Weighted_Ref") = Ebar.L/(e0*Ybar.L);
intensity(t,"Weighted") = E.L(t)/(e0*Y.L(t));
profit(t,"Weighted") = (py(t)*Y.L(t) - pl(t)*L.L(t) - pe(t)*E.L(t))/(srvshr(t)*pref(t));

pKbar.FX = card(t)/sum(t, pe(t)/pref(t));
pLbar.FX = card(t)/sum(t, pe(t)/pref(t));
solve puttyputty using nlp maximizing obj;
intensity(t,"Unweighted_Ref") = Ebar.L/(e0*Ybar.L);
intensity(t,"Unweighted") = E.L(t)/(e0*Y.L(t));
profit(t,"Unweighted") = (py(t)*Y.L(t) - pl(t)*L.L(t) - pe(t)*E.L(t))/(srvshr(t)*pref(t));


$if "%batch%"=="yes" $setglobal gp_opt2 "set output 'fig2.png'"
$libinclude plot intensity

$if "%batch%"=="yes" $setglobal gp_opt2 "set output 'fig3.png'"
$libinclude plot profit

*       Repeat the experiment with highly variable energy prices:

pe(t) = uniform(0.5,1.5)*pe(t);

pe(t) = pe(t)/pref(t);
$if "%batch%"=="yes" $setglobal gp_opt2 "set output 'fig4.png'"
$libinclude plot pe
pe(t) = pe(t)*pref(t);

pLbar.LO = 0.01 * 1;
pKbar.LO = 0.01 * 1;
pKbar.UP = +INF;
pLbar.UP = +INF;

solve puttyputty using nlp maximizing obj;

profit(t,"Dynamic") = (py(t)*Y.L(t) - pl(t)*L.L(t) - pe(t)*E.L(t))/(srvshr(t)*pref(t));
intensity(t,"Dynamic") = E.L(t)/(e0*Y.L(t));
intensity(t,"Dyn_Ref") = Ebar.L/(e0*Ybar.L);

*       Now do a calculation in which we use a weighted average of energy prices
*       to select a point on the ex-ante isoquant:

pKbar.FX = sum(t,srvshr(t))/sum(t, srvshr(t)*pe(t)/pref(t));
pLbar.FX = sum(t,srvshr(t))/sum(t, srvshr(t)*pe(t)/pref(t));

*       After having fixed the point, we then compute the time path of response
*       to market prices:

solve puttyputty using nlp maximizing obj;
intensity(t,"Weighted_Ref") = Ebar.L/(e0*Ybar.L);
intensity(t,"Weighted") = E.L(t)/(e0*Y.L(t));
profit(t,"Weighted") = (py(t)*Y.L(t) - pl(t)*L.L(t) - pe(t)*E.L(t))/(srvshr(t)*pref(t));

pKbar.FX = card(t)/sum(t, pe(t)/pref(t));
pLbar.FX = card(t)/sum(t, pe(t)/pref(t));
solve puttyputty using nlp maximizing obj;
intensity(t,"Unweighted_Ref") = Ebar.L/(e0*Ybar.L);
intensity(t,"Unweighted") = E.L(t)/(e0*Y.L(t));
profit(t,"Unweighted") = (py(t)*Y.L(t) - pl(t)*L.L(t) - pe(t)*E.L(t))/(srvshr(t)*pref(t));

$if "%batch%"=="yes" $setglobal gp_opt2 "set output 'fig5.png'"
$libinclude plot intensity

$if "%batch%"=="yes" $setglobal gp_opt2 "set output 'fig6.png'"
$libinclude plot profit