$title	Dynamic model for double dividend analysis.
display "com: Dynamic model for double dividend analysis.";
$ontext
$Id: dd_model.gms,v 1.1 2005/08/04 08:21:38 st Exp $
Copyright (C)	2004 Shiro Takeda
Time-stamp:	<2005-07-12 21:21:15 Shiro Takeda>
Author:		Shiro Takeda
First-written:	<2003/04/24>

$offtext
*	------------------------------------------------------------------
*	Important parameters.

*	Scenario type
$if not setglobal scen $setglobal scen base

*	The number of periods.
$if not setglobal tyear $setglobal tyear 40

*	Emission reduction rate (%).
$if not setglobal rdrate $setglobal rdrate 100

*	Trade balance.  Non-zero means balanced trade at each period.
$if not setglobal bopcon $setglobal bopcon 0

*	Quadratic adjustment cost.
$if not setglobal phi $setglobal phi 0.5

*	Uncompensated elasticity of labor supply:
$if not setglobal eols $setglobal eols 0.19

*	Intertemporal elasticity:
$if not setglobal sigu $setglobal sigu 0.5

*	Coefficient for Armington elasticity:
$if not setglobal c_siga $setglobal c_siga 1

*	Coefficient for energy eos:
$if not setglobal c_sige $setglobal c_sige 1

*	Government budget balance.  Non-zero means balanced budget at each
*	period.
$if not setglobal gbb $setglobal gbb 0

*	Labor growth rate
$if not setglobal gr $setglobal gr 0

*	Solver choice for calibration.
$if not setglobal mcp $setglobal mcp 1

*	Emission permit or carbon tax
$if not setglobal ca_tax $setglobal ca_tax 0

*	Solver choice
option mcp = path;
* option mcp = miles;

*	Digit.
option decimals = 8;

option sysout = on;

*	Define parameters.
parameter
    sol_mcp	Flag for quadratic programming
    rdrate	Reduction rate;
sol_mcp = %mcp%;
rdrate = %rdrate%;
display sol_mcp, rdrate;

parameter
    bopcon		Balance of payments restriction		/ %bopcon% /
    capflo		Free capital flow
    phi			Adjustment speed			/ %phi% /
    c_siga		Coefficient for Armington elasticity	/ %c_siga% /
    c_sige		Coefficient for energy eos 		/ %c_sige% /
    gbb			Government balanced budget		/ %gbb% /
    dcalib		/ 0 /
    ca_tax		Flag for carbon tax			/ %ca_tax% /
;
*	Capital flow constraint: If balance of payment restriction is not
*	imposed, capital flow is free.
bopcon = 1$%bopcon%;
capflo = yes$(not bopcon);
display bopcon, capflo, phi, gbb, c_siga, c_sige;

display "com: Defining sets and parameters for dynamic model.";
set
    t		Period index		/ 1*%tyear% /
    tfirst(t)	First period
    tlast(t)	Last period
;
tfirst(t) = yes$(ord(t) eq 1);
tlast(t) =  yes$(ord(t) eq card(t));
display t, tfirst, tlast;
alias (t,tt);

*	------------------------------------------------------------------
*	Import set definitions.
display "com: Import set definitions.";
$include set_definition

*	------------------------------------------------------------------
*	Import elasticity parameters.
display "com: Import elasticity parameters";
$include elasticity_parameters.gms

*	------------------------------------------------------------------
*	Import benchmark data.
display "com: Import benchmark data.";
$include dataset.dat
*	Endogenous sector.
display "Endogenous sector";
display iov, iovc, iovnc, ioqc, ioqnc;

*	Value added sector.
display "Value added sector";
display nlabinc, ncapinc, labtax, captax, nitax;

*	Final demand sector.
display "Final demand sector";
display cons, gexp, inve, expo, impo, impt, ctax, csub;

*	Household data.
display "Household data";
display lincome, cincome, linctax, cinctax, hctax, hcsub, hsave, hntrans;

*	------------------------------------------------------------------
*	Check imported data.
display "com: Check imported data.";

parameter
    ioq		Quantity data.
    chk_fdva	Check FD and VA data
    chk_out	Check output data
    chk_out_p	-1 if no output;

*	Quantity data.
ioq(e,i) = ioqc(e,i) + ioqnc(e,i);
display ioq;

*	Final demand and value added.
chk_fdva(i,"va") = nlabinc(i) + labtax(i) + ncapinc(i) + captax(i)
    + nitax(i);
chk_fdva("sum","va") = round(sum(i, chk_fdva(i,"va")), 6);
chk_fdva(i,"fd")
    = cons(i) + csub(i) + gexp(i) + inve(i)
    + expo(i) + impo(i) + impt(i);
chk_fdva("sum","fd") = round(sum(i, chk_fdva(i,"fd")), 6);

chk_fdva("diff","fd-va")
    = round(chk_fdva("sum","fd")
    - chk_fdva("sum","va"), 6);
display chk_fdva;
abort$chk_fdva("diff","fd-va") "FD and VA are inconsistent.";

*	Ouput data.
chk_out(i,"sec")
    = sum(ii, iov(ii,i))
    + nlabinc(i) + labtax(i) + ncapinc(i) + captax(i)
    + nitax(i);
chk_out(i,"com")
    = sum(ii, iov(i,ii))
    + cons(i) + csub(i) + gexp(i) + inve(i)
    + expo(i) + impo(i) + impt(i);
chk_out("sum","sec") = sum(i, chk_out(i,"sec"));
chk_out("sum","com") = sum(i, chk_out(i,"com"));
chk_out(i,"sec_com") = round(chk_out(i,"sec") - chk_out(i,"com"), 6);
chk_out("sum","sec_com")
    = round(chk_out("sum","sec")
    - chk_out("sum","com"), 6);
display chk_out;
abort$chk_out("sum","sec_com") "Output data are inconsistent.";

chk_out(i,"sec") = round(chk_out(i,"sec"), 6);
chk_out(i,"com") = round(chk_out(i,"com"), 6);
chk_out_p(i,"sec") = -1;
chk_out_p(i,"com") = -1;
chk_out_p(i,"sec")$chk_out(i,"sec") = 1;
chk_out_p(i,"com")$chk_out(i,"com") = 1;
option chk_out_p:0;
display chk_out_p;

*	----------------------------------------------------------------------
*	Parameters for technology improvements:
display "com: Parameters for technology improvements:";
parameter
    tpf	Technology coefficient for primary factors
    tene		Technology coefficient for energy goods

    t_int	Intermediate
;
tpf(t) = 1;
tene(t) = 1;
display tpf, tene;

t_int(t,i) = 1;

*	------------------------------------------------------------------
*	Production side
display "com: Production side";
parameter
    int0(i,ii)		Intermediate inputs
    vld0(i)		Value of labor inputs (net of tax)
    vkd0(i)		Value of capital inputs (net of tax)
    vpf0(i)		Value of primary factor (value added)
    ve0(i)		Value of energy inputs
    vpfe0(i)		Value of primary factor and energy inputs

    outs0(i)		Sector output (gross of tax)
    outc0(i)		Commodity output (gross of tax)

    d0(i)		Value of domestic supply
    x0(i)		Value of export supply
    out0(i)		Value of output (gross of tax)

    tq0(i)		Net indirect tax rate on production
    tl0(i)		Labor tax rate on production
    tk0(i)		Capital tax rate on production
;
*	Intermediate inputs.
int0(i,ii) = iov(i,ii);

*	Primary factor input (value).
vld0(i) = nlabinc(i);
vkd0(i) = ncapinc(i);

*	Sector output.
outs0(i) = sum(ii, iov(ii,i))
    + nlabinc(i) + ncapinc(i)
    + labtax(i) + captax(i) + nitax(i);
outs0(i) = outs0(i);

*	Commodity output.
outc0(i) = sum(ii, iov(i,ii))
    + cons(i) + csub(i) + gexp(i)
    + inve(i) + expo(i) + impo(i) + impt(i);
outc0(i) = outc0(i);

*	Output to export market.
x0(i) = expo(i);

*	Output to domestic market.
d0(i) = round(outc0(i) - x0(i), 8);

*	Total ouput.
out0(i) = d0(i) + x0(i);

display int0, vld0, vkd0, outs0, outc0, x0, d0, out0;

chk_out_p(i,"sec") = -1;
chk_out_p(i,"com") = -1;
chk_out_p(i,"sec")$outs0(i) = 1;
chk_out_p(i,"com")$outc0(i) = 1;
display chk_out_p;

*	Gross value added.
vpf0(i)$outs0(i) = ncapinc(i) + captax(i) + nlabinc(i) + labtax(i);
display vpf0;

*	Labor tax rate.
tl0(i)$vld0(i) = labtax(i) / vld0(i);

*	Capital tax rate.
tk0(i)$vkd0(i) = captax(i) / vkd0(i);

*	Indirect tax rate.
tq0(i)$outs0(i) = nitax(i) / outs0(i);
display tl0, tk0, tq0;

*	------------------------------------------------------------------
*	Investment
display "com: Investment";
parameter
    i0		Total value of gross investment
    ishr(i)	Share of each good in investment
    id0(i)	Gross investment demand;

id0(i) = inve(i);
i0 = sum(i, id0(i));
ishr(i) = id0(i)/i0;
display id0, i0, ishr;

*	------------------------------------------------------------------
*	Government demand
display "com: Government demand";
parameter
    gd0(i)	Government demand (including gov investment)
    g0		Government total expenditure;

gd0(i) = gexp(i);
g0 = sum(i, gd0(i));
display gd0, g0;

*	------------------------------------------------------------------
*	Trade and Armington.
display "com: Trade and Armington";
parameter
    a0(i)	Armington supply
    m0(i)	Demand for imported good (net of import tax)
    tm0(i)	Import tax rate
    bop0	Balance of payments of each good (ex - im)
    bopt0	Balance of payments (ex - im);

m0(i) = - impo(i);
tm0(i)$m0(i) = - impt(i)/m0(i);
a0(i) = d0(i) + (1+tm0(i)) * m0(i);
bop0(i) = x0(i) - m0(i);
bopt0 = round(sum(i, bop0(i)), 8);
display m0, tm0, a0, bop0, bopt0;

*	------------------------------------------------------------------
*	Household data
display "com: Household data";

parameter
    epi0	Expanded disposal income
    wg0		Household gross current expenditure
    wn0		Household net current expenditure
    c0		Household consumption expenditure
    ce0		Energy good consumption
    cne0	Non-energy good consumption
    cg0		Household gross consumption expenditure
    cd0		Household consumption demand

    li0		Value of labor income (gross of tax)
    ki0		Value of capital income (gross of tax)

    tli0	Labor income tax rate
    tci0	Capital income tax rate
    ts0		Consumption subsidy rate
    tc0		Consumption tax rate

    vos		Value of savings
    vot		Value of net transfer
;
li0 = lincome;
ki0 = cincome;
vos = hsave;
vot = hntrans;
cd0(i) = cons(i) + csub(i);
display li0, ki0, vos, cd0, vot;

ts0(i)$cd0(i) = - csub(i) / cd0(i);
tc0(i)$cons(i) = ctax(i) / cons(i);
c0 = sum(i, cons(i));
cg0 = sum(i, (1+tc0(i)) * (1+ts0(i)) * cd0(i));
tli0 = linctax / li0;
tci0 = cinctax / ki0;
display ts0, tc0, cg0, tli0, tci0;

cne0 = sum(nene, cons(nene));
ce0 = sum(ene, cons(ene));
display cne0, ce0;

*	------------------------------------------------------------------
*	Total labor endowment.
display "com: Total labor endowment.";
$ontext
Total labor endowment = a * labor supply
Leisure = total labor endowment - labor supply

Assumptions:

* Labor hours per month are 190 hours.
* Time available for leisure and work is 16 hours per day.

Then, total time available for leisure and work per month is given by

	16 * 30 = 480

Thus, parameter a is

	a = 480 / 190 = 2.53

$offtext
parameter
    lratio	Labor supply ratio		/ 2.53 /
    endl	Benchmark total labor endowment
    lei0	Benchmark leisure
;
endl = lratio * li0;
lei0 = endl - li0;
display lratio, endl, lei0;

*	------------------------------------------------------------------
*	Defining sets and parameters for dynamic model.

*	-------------------------------------------------------------------
*	Calibrate base year variables.
$ontext
Calibrate benchmark and baseline variables to make the baseline equilibrium
the steady state.
$offtext
parameter
    vks0	Benchmark total value of capital stock
    vai0	Benchamrk value of capital income (net of tax)
    vat0	Benchmark capital income tax;
vks0 = sum(i, vkd0(i));
vat0 = cinctax;
vai0 = cincome - vat0;
display vks0, vat0, vai0;

*	Exogenous parameter.
parameter
    gr		Baseline labor growth rate	/ %gr% /;
display gr;

*	Calibrated parameters
parameter
    r		Baseline interest rate
    delta	Depreciation rate
    j0		Benchmark net investment
    jd0		Benchmark net investment demand

    rk0		Benchmark rental price for producer
    rke0	Benchmark rental price for household
    ks0		Benchmark capital stock
    pk0		Benchmark capital price
    pka0	Reduced cost of investment (adjustment premium)

    ld0		Benchmark labor input (quantity)
    kd0		Benchmark capital input (quantity);

*	Target values for calibrated variables.
parameter
    rt0		Interest rate (target level)	    / 0.03 /
    delt0	Depreciation rate (target level)    / 0.07 /;
display rt0, delt0;

*	Variable definitions.
variables
    rt		Interest rate
    delt	Depreciation rate
    rket	Rental rate for household
    kt		Capital stock

    jt		Net investment
    loss	Loss function

    lam_1	Multiplier
    lam_2	Multiplier
    lam_3	Multiplier;

*	Equation definitions.
equations
    eq_obj	Objetive function
    eq_con1	Constraint 1
    eq_con2	Constraint 2
    eq_con3	Constraint 3
    eq_con4	Constraint 4

    eq_rt	FOC
    eq_delt	FOC
    eq_rket	FOC
    eq_kt	FOC;
*	------------------------------------------------------------------
*	NLP model.
eq_obj .. loss =e=
	  power((rt - rt0)/rt0, 2) + power((delt - delt0)/delt0, 2);

eq_con1 .. rket =e= (1+phi*(gr+delt))*(rt+delt)
	   - phi*(delt+gr)**2 / 2;

eq_con2 .. rket*kt =e= vai0;

eq_con3 .. i0 =e= (delt+gr) * kt * (1+phi*(delt+gr)/2);

model calib_nlp NLP model / eq_obj, eq_con1, eq_con2, eq_con3 /;

*	------------------------------------------------------------------
*	MCP model.
eq_rt .. 2*(rt-rt0)/(rt0**2) - lam_1*(1+phi*(gr+delt)) =e= 0;

eq_delt .. 2*(delt-delt0)/(delt0**2)
	   + lam_1*( phi*(gr-rt) - (1+phi*(gr+delt)) )
	   - lam_3*kt*(1+phi*(delt+gr)) =e= 0;

eq_rket .. lam_1 + lam_2*kt =e= 0;

eq_kt .. lam_2*rket - lam_3*(delt+gr)*(1 + phi*(delt+gr)/2) =e= 0;

model calib_mcp MCP model / eq_con1, eq_con2, eq_con3, eq_rt, eq_delt,
			    eq_rket, eq_kt /;

*	Initial values.
rt.l = rt0;
delt.l = delt0;
kt.l = vks0 / (delt.l + rt.l);
rket.l = rt.l+delt.l;
lam_1.l = 0;
lam_2.l = 0;
lam_3.l = 0;

parameter	calibres	Calibrted results;
if(sol_mcp,
    display "com: Calibration by MCP";
    solve calib_mcp using mcp;
    calibres("r","cal") = rt.l;
    calibres("rke","cal") = rket.l;
    calibres("delta","cal") = delt.l;
    calibres("k0","cal") = kt.l;
    calibres("pk0","cal") = 1+phi*(delt.l+gr);
    calibres("vks0","cal")
	= calibres("k0","cal") * calibres("rke","cal");
    calibres("j0","cal") = (delt.l+gr)*kt.l;
    calibres("pka0","cal") = phi*(delt.l+gr)**2 / 2;
    calibres("adjc","cal") = phi*(delt.l+gr) / 2;
    display calibres;
    else
    display "com: Calibration by NLP";
    solve calib_nlp using nlp minimizing loss;
    calibres("r","cal") = rt.l;
    calibres("rke","cal") = rket.l;
    calibres("delta","cal") = delt.l;
    calibres("k0","cal") = kt.l;
    calibres("pk0","cal") = 1+phi*(delt.l+gr);
    calibres("vks0","cal")
	= calibres("k0","cal") * calibres("rke","cal");
    calibres("j0","cal") = (delt.l+gr)*kt.l;
    calibres("pka0","cal") = phi*(delt.l+gr)**2 / 2;
    calibres("adjc","cal") = phi*(delt.l+gr) / 2;
    display calibres;
);
r = rt.l;
delta = delt.l;
rke0 = rket.l;
rk0 = rke0/(1-tci0);
ks0 = kt.l;
pk0 = calibres("pk0","cal");
pka0 = calibres("pka0","cal");
j0 = calibres("j0","cal");
display "Display calibrated parameters";
display r, delta, rke0, rk0, ks0, pk0, pka0, j0;

*	Benchmark capital inputs.
kd0(i) = ks0*vkd0(i)/vks0;
display kd0;

*	Benchmark net investment demand.
jd0(i) = ishr(i)*j0;
display jd0;

*	wrt labor, quantity = value.
ld0(i) = vld0(i);
display ld0;

*	------------------------------------------------------------------
*	Reference prices.
display "com: Reference prices.";
parameter
    pl0		Reference wage rate for industry
    rk00	Reference rental price for industry
    rk000	Reference rental price for industry

    ple0	Reference labor price for household (price of leisure)
    rke0	Reference rental price for household
    pm0		Reference price of imported goods
    pc0		Reference good price for household
    pc00	Reference good price for household

    pd0		Reference price for domestic supply
    px0		Reference price for export supply
    pq0		Reference price for output;

pl0(i) = 1+tl0(i);
rk00(i) = (1+tk0(i))*rk0;
rk000(i) = (1+tk0(i));

ple0 = 1-tli0;
rke0 = (1-tci0)*rk0;

pc0(i) = (1+ts0(i));
pc00 = 1.03;
pm0(i) = 1+tm0(i);

px0(i) = 1-tq0(i);
pd0(i) = 1-tq0(i);
pq0(i) = 1-tq0(i);
display pl0, rk00, ple0, rke0, pc0, pm0, px0, pd0;

*	Benchmark current expenditure.
wg0 = cg0 + ple0*lei0;
display wg0;

*	Benchmark expanded disposal income = labor income + capital income -
*	saving.
epi0 = ple0*endl + rke0*ks0 - vos + vot;
display epi0;
parameter	chk_revex	Check revenue and expenditure consistency;
chk_revex("rev-exp") = round(wg0 - epi0, 3);
display chk_revex;
abort$chk_revex("rev-exp") "Household revenue and expenditure are inconsistent";

*	------------------------------------------------------------------
*	Baseline path.
display "com: Baseline path.";
parameter
    qref(t)	Reference quantity path
    pref(t)	Reference price path
    wref(t)	Reference (period) utility path
    lref(t)	Labor endowment path
;
*	Reference quantity path.
qref(t) = (1+gr)**(ord(t)-1);

*	Reference price path.
pref(t) = (1/(1+r))**(ord(t)-1);

*	Reference utility path
wref(t) = wg0 * qref(t);

*	Labor endowment path
lref(t) = endl * qref(t);

display qref, pref, wref, lref;


parameter
    gqref(t)	Reference government period expenditure
    gbar	Benchmark government lifetime expenditure
;gqref(t) = g0*qref(t);
gbar = sum(t, pref(t) * gqref(t));
display gqref, gbar;

parameter    ubar	Benchmark lifetime utility;
*	Prsigent value of utility over the model horizon:
ubar = sum(t, wref(t) * pref(t));
display ubar;

parameter	adj0;
adj0 = phi*(delta+gr)/2 + sum(i, ishr(i)*phi*(delta+gr)/2);
display adj0;

*	------------------------------------------------------------------
*	Baseline tax rate.
display "com: Baseline tax rate..";
parameter
    tq(i)	Indirect tax rate
    tm(i)	Import tax rate
    tl(i)	Tax rate on labor input in production
    tk(i)	Tax rate on capital input in production

    tli		Labor income tax rate
    tci		Capital income tax rate
    ts(i)	Consumption subsidy rate
    tc		Consumption tax rate

    tca(t,e)	Carbon tax rate;

*	Initial values.
tq(i) = tq0(i);
tm(i) = tm0(i);
tl(i) = tl0(i);
tk(i) = tk0(i);
tc = 0.03;
tli = tli0;
tci = tci0;
ts(i) = ts0(i);
tca(t,e) = 0;
display tq, tm, tl, tk, tc, tli, tci, ts, tca;

*	------------------------------------------------------------------
*	Calibrate leisure-consumtpion elasticity of substitution.
display "com: Calibrate leisure-consumtpion elasticity of substitution.";
$ontext
Elasticity of substitution between consumption and leisure is calibrated from
labor supply elasticity.
$offtext

sigcl = ( eols * li0 / lei0 + ple0 * li0 / epi0 )
    / (1 - ple0 * lei0 / epi0 );
display eols, sigcl;

*	Check calibration.
parameter sigcl_, sh_l, sh_le, sh_led;

sh_l = ple0*li0 / epi0;
sh_led = ple0*endl / epi0;
sh_le = ple0*lei0 / epi0;
sigcl_("sigcl") = ( eols*sh_l/(sh_led - sh_l) + sh_l ) / (1 - sh_le);
sigcl_("diff") = round(sigcl_("sigcl") - sigcl, 6);
display sh_l, sh_le, sh_led, sigcl_;
abort$sigcl_("diff") "Inconsistent in calibration";

*	------------------------------------------------------------------
*	Check data.
display "com: Check data.";

*	Check final demand and value added.
chk_fdva(i,"va") = rk00(i)*kd0(i) + pl0(i)*ld0(i)
    + tq0(i)*outs0(i);
chk_fdva("sum","va") = round(sum(i, chk_fdva(i,"va")), 6);

chk_fdva(i,"fd") = cd0(i) + id0(i) + gd0(i) + bop0(i) - tm0(i)*m0(i);
chk_fdva("sum","fd") = round(sum(i, chk_fdva(i,"fd")), 6);

chk_fdva("diff","fd-va")
    = round(chk_fdva("sum","fd") - chk_fdva("sum","va"), 6);
display chk_fdva;
abort$chk_fdva("diff","fd-va") "FD and VA are inconsistent.";

*	Check output data.
chk_out(i,"sec") = sum(ii, int0(ii,i)) + chk_fdva(i,"va");
chk_out(i,"com") = sum(ii, int0(i,ii)) + chk_fdva(i,"fd");
chk_out("sum","sec") = sum(i, chk_out(i,"sec"));
chk_out("sum","com") = sum(i, chk_out(i,"com"));
chk_out(i,"sec_com") = round(chk_out(i,"sec") - chk_out(i,"com"), 6);
chk_out("sum","sec_com")
    = round(chk_out("sum","sec") - chk_out("sum","com"), 6);
display chk_out;
abort$chk_out("sum","sec_com") "Output data are inconsistent.";

*	------------------------------------------------------------------
*	Check investment and savings consistency.
parameter chk_is	check investment and saving;
chk_is("save") = vos;
chk_is("ginv") = i0;
*	Saving = investment because benchmark balance of payment and
*	government fiscal balance are zero.
chk_is("save-ginv")
    = round(chk_is("save") - chk_is("ginv") - bopt0, 6);
display chk_is;
abort$chk_is("save-ginv") "Inconsistent";

*	------------------------------------------------------------------
*	Check household income and expenditure.
parameter chk_ie	Check household income and expenditure;
*	Expanded income.
chk_ie("inc") = ple0*endl + rke0*ks0 - vos + vot;
*	Expanded expenditure.
chk_ie("exp") = cg0 + ple0*lei0;
chk_ie("inc-exp")
    = round(chk_ie("inc") - chk_ie("exp"), 6);
display chk_ie;
abort$chk_ie("inc-exp") "Income and expenditure are inconsistent!";

*	------------------------------------------------------------------
*	Check government revenue and expenditure
parameter chk_gov	Check government revenue and expenditure;

*	Government revenue.
chk_gov("rev")
    = sum(i, tk0(i) * rk0 * kd0(i)
    + tl0(i) * ld0(i)
    + tq0(i) * (x0(i)+d0(i)))
    + tli0 * li0
    + tci0 * rk0 * ks0
    + sum(i, tm0(i) * m0(i))
    + sum(i, ts0(i) * cd0(i))
    + sum(i, tc0(i) * (1 + ts0(i)) * cd0(i))
    - vot;
chk_gov("exp")= g0;
chk_gov("rev-exp") = round(chk_gov("rev") - chk_gov("exp"), 6);
display chk_gov;
abort$chk_gov("rev-exp") "Government expenditure and revenue are inconsistent";

*	------------------------------------------------------------------
*	Adjust energy (quantity) data.
display "com: Adjust energy (quantity) data.";
*	------------------------------------------------------------------
*	Check imported energy data.
display "com: Check imported energy data.";
display ioqc, ioqnc, fdq;

parameter cacoef "Carbon coefficient (t-C/10^7kcal)" /
	  lim	1.2000000000000E-01
	  coc	1.0451149310000E+00
	  sla	1.0154523000000E+00
	  cru	7.9152314900000E-01
	  nat	5.8507665700000E-01
	  pet	7.8289055874027E-01
	  cok	1.2310959200000E+00
	  gas	5.9652317200000E-01
	  /
	  ioca	Carbon emissions (MtC);

*	Calculate carbon emissions.
ioca(e,i) = cacoef(e) * ioqc(e,i);
ioca(e,"cons") = cacoef(e) * fdq(e,"cons");

ioca(e,"end") = sum(i, ioca(e,i));
ioca(e,"fd") = ioca(e,"cons");
ioca(e,"sum") = ioca(e,"end") + ioca(e,"fd");

ioca("sum",i) = sum(e, ioca(e,i));
ioca("sum","cons") = sum(e, ioca(e,"cons"));
ioca("sum","sum") = sum(e, ioca(e,"sum"));
display ioca;

*	------------------------------------------------------------------
*	Check consistency between value and quantity data.
display "com: Check consistency between value and quantity data.";
*	------------------------------------------------------------------
*	Check the value and quantity consistency
display "com: Check the value and quantity consistency";
$ontext
There are three cases:

(1) Both value data and quantity data are positive.
(2) Value data are positive, but quantity data are zero.
(3) Value data are zero, but quantity data are positive.

Set quantity data to zero for case (2) and (3).

By this adjustment, carbon emissions become less than the real value.
$offtext
parameter
    ch_iov
    ch_iovc
    ch_iovnc

    ch_ioq
    ch_ioqc	Int input for combustion
    ch_ioqnc	Int input for non-combustion
    ch_fdq	Final demand (quantity)

    ch_flag	Flag for check: (0 or 1)
    ch_flag_c	Flag for check: (0 or 1)
    ch_flag_nc	Flag for check: (0 or 1)
    ch_flag_fd	Flag for check: (0 or 1);

*	Check sign consistency.
ch_iov(e,i) = iov(e,i);
ch_iovc(e,i) = iovc(e,i);
ch_iovnc(e,i) = iovnc(e,i);

ch_ioq(e,i) = ioq(e,i);
ch_ioqc(e,i) = ioqc(e,i);
ch_ioqnc(e,i) = ioqnc(e,i);

ch_flag(e,i) = 0;
ch_flag_c(e,i) = 0;
ch_flag_nc(e,i) = 0;

*	Intermediate inputs.
loop((e,i)$(ch_iov(e,i) or ch_ioq(e,i)),
    if((ch_iov(e,i) * ch_ioq(e,i) > 0),
	ch_flag(e,i) = 0;
	elseif (ch_iov(e,i) * ch_ioq(e,i) = 0),
	ch_flag(e,i)$ch_iov(e,i) = 1;
	ch_flag(e,i)$ch_ioq(e,i) = -1;
    ));
option ch_flag:0;
display ch_flag;

*	Combustion purpose.
loop((e,i)$(ch_iovc(e,i) or ch_ioqc(e,i)),
    if((ch_iovc(e,i) * ch_ioqc(e,i) > 0),
	ch_flag_c(e,i) = 0;
	elseif (ch_iovc(e,i) * ch_ioqc(e,i) = 0),
*		Value data only.
	ch_flag_c(e,i)$ch_iovc(e,i) = 1;
	ioqc(e,i)$ch_iovc(e,i) = 0;
*		Quantity data only.
	ch_flag_c(e,i)$ch_ioqc(e,i) = -1;
	ioqc(e,i)$ch_ioqc(e,i) = 0;
    ));

*	Non-combustion purpose.
loop((e,i)$(ch_iovnc(e,i) or ch_ioqnc(e,i)),
    if((ch_iovnc(e,i) * ch_ioqnc(e,i) > 0),
	ch_flag_nc(e,i) = 0;
	elseif (ch_iovnc(e,i) * ch_ioqnc(e,i) = 0),
*		Value data only.
	ch_flag_nc(e,i)$ch_iovnc(e,i) = 1;
	ioqnc(e,i)$ch_iovnc(e,i) = 0;
*		Quantity data only.
	ch_flag_nc(e,i)$ch_ioqnc(e,i) = -1;
	ioqnc(e,i)$ch_ioqnc(e,i) = 0;
    ));

option ch_flag_c:0, ch_flag_nc:0;
display ch_flag_c, ch_flag_nc;
display "After adjustment", ioqc, ioqnc;

*	Calculate carbon emissions.
ioca(e,i) = cacoef(e) * ioqc(e,i);
ioca(e,"cons") = cacoef(e) * fdq(e,"cons");

ioca(e,"end") = sum(i, ioca(e,i));
ioca(e,"fd") = ioca(e,"cons");
ioca(e,"sum") = ioca(e,"end") + ioca(e,"fd");

ioca("sum",i) = sum(e, ioca(e,i));
ioca("sum","cons") = sum(e, ioca(e,"cons"));
ioca("sum","sum") = sum(e, ioca(e,"sum"));
display ioca;

*	------------------------------------------------------------------
*	Contribution ratio.
display "com: Contribution ratio";
parameter	rat_ene		Contribution ratio;
rat_ene(e,i)$(ioqc(e,i)+ioqnc(e,i))
    = ioqc(e,i)/(ioqc(e,i)+ioqnc(e,i));
display rat_ene;

*	------------------------------------------------------------------
*	Define energy inputs and consumption.
display "com: Define energy inputs and consumption.";
parameter
    veic	Value of energy input (combustion purpose)
    veict	Sumal value of energy input (combustion purpose)
    veinc	Value of energy input (non-combustion purpose)
    vec		Value of energy consumption
    vect	Value of energy for combustion purpose;

veic(e,i) = rat_ene(e,i) * int0(e,i);
veic("ele",i) = int0("ele",i);
veic(e,"sum") = sum(i, veic(e,i));
veinc(e,i) = int0(e,i) - veic(e,i);
display veic, veinc;

vec(e) = cd0(e);
display vec;

veict(i) = sum(ec, veic(ec,i)) + veic("ele",i);
display veict;

vect(e) = veic(e,"sum") + vec(e);
display vect;

*	Value of energy inputs
ve0(i)$outs0(i) = sum(ec, veic(ec,i)) + sum(ele, int0(ele,i));
display ve0;

*	Primary factor and energy.
vpfe0(i)$outs0(i) = vpf0(i) + ve0(i);
display vpfe0;

*	------------------------------------------------------------------
*	Deriving carbon emission data.
display "com: Deriving carbon emission data.";

parameter
    carbi0	Carbon emissions from intermediate inputs (MtC)
    carbc0	Carbon emissions from final consumption (MtC)
    carbt0	Carbon emissions (MtC)

    carbl	Carbon emission restriction (MtC)
    carbl0	Benchmark carbon emissions (MtC)

    scale_ca	Scale factor 	/ 1000 /
;
*	Carbon emissions from intermediate inputs.
carbi0(e,i) = cacoef(e)*ioqc(e,i) / scale_ca;
carbi0(e,"sum") = sum(i, carbi0(e,i));
carbi0("sum",i) = sum(e, carbi0(e,i));
carbi0("sum","sum") = sum((e,i), carbi0(e,i));

*	Carbon emissions from final demand (consumption).
carbc0(e) = cacoef(e)*fdq(e,"cons") / scale_ca;
carbc0("sum") = sum(e, carbc0(e));
display carbi0, carbc0;

carbt0(e) = carbi0(e,"sum") + carbc0(e);
carbt0("sum") = sum(e, carbt0(e));
display carbt0;

*	Benchmark carbon emissions.
carbl0 = sum((e,i), carbi0(e,i)) + sum(e, carbc0(e));
*	Carbon emission restriction is set to zero tentatively.
* carbl(t) = 0;
carbl(t) = carbl0 * 1e+6;
option carbl0:6, carbl:6;
display carbl0, carbl;

parameter	chk_carbon	Check carbon emissions (MtC);

*	Carbon emissions from intermediate inputs.
chk_carbon(e,i) = carbi0(e,i);
chk_carbon(e,"int") = sum(i, chk_carbon(e,i));
chk_carbon("sum","int") = sum(e, chk_carbon(e,"int"));
chk_carbon("sum",i) = sum(e, chk_carbon(e,i));
*	Carbon emissions from consumption.
chk_carbon(e,"con") = carbc0(e);
chk_carbon("sum","con") = sum(e, chk_carbon(e,"con"));
chk_carbon(e,"sum") = chk_carbon(e,"int") + chk_carbon(e,"con");
*	Total carbon emissions.
chk_carbon("sum","sum") = sum(e, chk_carbon(e,"sum"));
option chk_carbon:3;
display chk_carbon;

*	------------------------------------------------------------------
*	Scale parameter
display "com: Scale parameter";
parameter
    scale_k	Scale for benchmark capital stock		/ 1 /
    scale_l 	Scale for labor endowment			/ 1 /
    scale_t 	Scale for transfer to household			/ 1 /
    scale_g 	Scale for government lifetime expenditure	/ 1 /;
display scale_k, scale_l, scale_t, scale_g;

*	------------------------------------------------------------------
*	Tax instruments.
display "com: Tax instruments.";
set	taxinst		Tax instruments
	/
 	lum	Lump-sum tax
	lin	Labor income tax
	lab	Labor tax on production
	cin	Capital income tax
	cap	Capital tax on production
	ctx	Consumption tax
	imp	Import tax
* 	bau	No carbon limit
* 	labu	Labor tax on production (uniform)
* 	capu	Capital tax on production (uniform)
	  /
	tins(taxinst)	Tax instrument;
tins(taxinst) = no;
tins("lum") = yes;
display taxinst, tins;

*	------------------------------------------------------------------
*	MCP code
display "com: MCP code";

*	------------------------------------------------------------------
*	Share parameters
display "com: Share parameters";

parameter
    cost_q(i)		Unit cost for production

    sh_x(i)		Share of export in output
    sh_fe(i)		Share of primary factor and energy in production cost
    sh_f(i)		Share of primary factor in PF-E aggregation
    sh_fl(i)		Share of labor in primary factor
    sh_e		Share of energy goods in production
    sh_nen(ii,i) 	Share of non-energy intermediate inputs
    sh_nc(e,i)		Share of energy for non-combustion purpose
    sh_cl(cl,i)		Share of cok and lim for combustion purpose

    sh_nene		Share of non-energy goods in consumption.
    sh_c(i)		Share of each good in non-energy consumption.
    sh_ec(i)		Share of each good in energy consumption.
    sh_lei		Share of leisure
    sh_u(t)		Share of period utility
    sh_g(t)		Share of each good in government expenditure

    sh_ad(i)		Share of domestic good in Armington aggregation
    sh_inv(i)		Share of input in investment activity

    sh_g_gdp(t)		Share of government expenditure in GDP
;
*	Production cost (net of tax).
cost_q(i)$out0(i) = (1-tq0(i))*out0(i);

*	Share of export.
sh_x(i)$out0(i) = x0(i) / out0(i);

*	Share of primary factor-energy composite
sh_fe(i)$out0(i) = vpfe0(i) / cost_q(i);

*	Share of PF in PFE composite
sh_f(i)$vpfe0(i) = vpf0(i)/ vpfe0(i);

*	Share of labor in PF composite
sh_fl(i)$vpf0(i) = pl0(i)*ld0(i) / vpf0(i);

*	Share of energy in energy composite.
sh_e(ec,i)$veict(i) = veic(ec,i) / veict(i);

*	Share of electricity in energy composite
sh_e(ele,i)$veict(i) = veic(ele,i) / veict(i);

*	Share of non-energy inputs.
sh_nen(nen,i)$out0(i) = int0(nen,i) / cost_q(i);

*	Share of energy input for non-combustion purpose.
sh_nc(e,i)$out0(i) = veinc(e,i) / cost_q(i);

*	Share of cok and lim for combustion purpose.
sh_cl(cl,i)$out0(i) = veic(cl,i) / cost_q(i);

*	Government expenditure
sh_g(t) = pref(t) * qref(t) / sum(tt, pref(tt) * qref(tt));

display cost_q, sh_x, sh_fe, sh_f, sh_fl, sh_e, sh_nen, sh_nc, sh_cl, sh_g;

*	Share of domestic goods in Armington aggregation.
sh_ad(i)$a0(i) = d0(i) / a0(i);
display sh_ad;

*	Share of each good in investment.
sh_inv(i)$j0 = jd0(i) / j0;
display sh_inv;

*	Share of leisure in total expenditure.
sh_lei = ple0*lei0 / wg0;

*	Share of non-energy goods in consumption.
sh_nene = cne0 / c0;

*	Share of each good in non-energy consumption.
sh_c(nene) = pc0(nene) * cd0(nene) / cne0;

*	Share of each good in energy consumption.
sh_ec(ene) = pc0(ene) * cd0(ene) / ce0;
display sh_lei, sh_nene, sh_c, sh_ec;

*	Share of each period utility.
sh_u(t) = pref(t)*wref(t) / ubar;
display sh_u;

*	Share of government expenditure in GDP:
sh_g_gdp(t) = g0 / chk_fdva("sum","fd");
display sh_g_gdp;

parameter
    pgle0	Reference price of lifetime government expenditure;
pgle0 = sum(t, sh_g(t) * pref(t));
display pgle0;

parameter
    roc_gdp	Target rate of change in GDP / 0 /
    roc_ca	Target rate of chagne in carbon emissions / 0 /
;
*	------------------------------------------------------------------
*	MCP formulation
display "com: MCP formulation";
variables
*	Unit demand and supply
    a_x(t,i)		Export supply
    a_d(t,i)		Domestic supply
    a_l(t,i)		Demand for labor
    a_k(t,i)		Demand for capital
    a_e(t,ii,i)		Demand for energy (exc. cok and lim) in production

    a_ec(t,i)		Consumption demand for energy goods
    a_c(t,i)		Consumption demand for non-energy goods
    a_cc(t)		Demand for aggregate consumption
    a_lei(t)		Demand for leisure
    a_u(t)		Demand for period utility

    a_ad(t,i)		Demand for domestic good from Armington activity
    a_am(t,i)		Demand for imported good from Armington activity

*	Activity variables
    q_q(t,i)		Activity level for production
    q_a(t,i)		Armington aggregation
    q_x(t,i)		Export activity
    q_m(t,i)		Import activity
    q_cc(t)		Consumption aggregation activity
    q_w(t)		Period utility
    q_j(t)		Net investment
    q_u			Lifetime utility
    q_k(t)		Capital accumulation (investment) activity
    q_g(t)		Government expenditure at each period
    q_cnc		Carbon emissions
    q_sl(t)		Labor supply

*	Unit cost.
    c_cc(t)		Unit cost of consumption aggregation
    c_a(t,i)		Unit cost of Armington aggregation
    c_w(t)		Unit cost of period utility (expenditure function)
    c_u			Unit cost of lifetime utility (expenditure function)
    c_q(t,i)		Unit cost of production
    c_e(t,i)		Price index of emission sources for combustion purpose

*	Price variables
    p_q(t,i)		Price index (unit revenue) for output
    p_f(t,i)		Price index of primary factor composite
    p_ee(t,i)		Price index of energy composite in production
    p_fe(t,i)		Price index of primary factor-energy composite
    p_c(t)		Price index of non-energy composite in consumption
    p_ec(t)		Price index of energy composite in consumption
    p_i(t)		Price index of investment good
    p_e(t,i)		Price index of emission sources for combustion purpose

    p_cc(t)		Price index of aggregate consumption
    p_d(t,i)		Price of domestic good
    p_x(t,i)		Price of export good
    p_m(t,i)		Price of import good
    p_fxv		Exchange rate (intertemporal)
    p_fx		Exchange rate
    p_l(t)		Wage rate for producer
    p_le(t)		Wage rate for household
    p_w(t)		Price of period utility
    p_u			Price of lifetime utility
    p_a(t,i)		Price of Armington good
    p_g(t)		Price of government expenditure
    p_k(t)		Shadow price of capital
    p_rk(t)		Rental price for industry
    p_rke(t)		Rental price for household
    p_gle		Price index for government lifetime consumption
    p_ca(t)		Price of emission permits
    p_ka(t)		Adjustmnet premium
    p_disc		Discount factor
    p_intr		"Interest rate 1/(1+r)"

*	Income variables
    v_inc_h		Extended lifetime income of household
    v_inc_g(t)		Government income (period)
    v_inc_gl		Government income (lifetime)

    q_tcap		post-terminal capital stock
    m_lum		Multiplier for lump-sum tax
    m_tl		Multiplier for labor tax
    m_tk		Multiplier for capital tax
    m_tli		Multiplier for labor income tax
    m_tci		Multiplier for capital income tax
    m_tc		Multiplier for consumption tax
    m_tm		Multiplier for import tariff

    t_pf(t)		Technology parameter for primary factor
    t_ene(t)		Technology parameter for energy
    roc_t_pf		Rate of change in t_pf
    roc_t_ene		Rate of change in t_ene

    v_gdp(t)		GDP
    c_emi(t)		Carbon emissions

    m_lum_p(t)		Lump-sum tax (period)
    share_g(t)		Share of period government expenditure
;
equations
*	Unit demand and supply
    e_a_x(t,i)		Export supply
    e_a_d(t,i)		Domestic supply
    e_a_l(t,i)		Demand for labor
    e_a_k(t,i)		Demand for capital
    e_a_e(t,ii,i)	Demand for energy (exc. cok and lim) in production

    e_a_ec(t,i)		Consumption demand for energy goods
    e_a_c(t,i)		Consumption demand for non-energy goods
    e_a_cc(t)		Demand for aggregate consumption
    e_a_lei(t)		Demand for leisure
    e_a_u(t)		Demand for period utility

    e_a_ad(t,i)		Demand for domestic good from Armington activity
    e_a_am(t,i)		Demand for imported good from Armington activity

*	Zero profit conditions
    e_q_q(t,i)		Activity level for production
    e_q_a(t,i)		Armington aggregation
    e_q_x(t,i)		Export activity
    e_q_m(t,i)		Import activity
    e_q_cc(t)		Consumption aggregation activity
    e_q_w(t)		Period utility
    e_q_g(t)		government expenditure at each period
    e_q_j(t)		Net investment
    e_q_u		Lifetime utility
    e_q_k(t)		Capital accumulation (investment) activity
    e_q_cnc		Carbon emissions
    e_q_sl(t)		Labor supply

*	Unit cost.
    e_c_cc(t)		Unit cost of consumption aggregation
    e_c_a(t,i)		Unit cost of Armington aggregation
    e_c_w(t)		Unit cost of utility (expenditure function)
    e_c_u		Unit cost of lifetime utility (expenditure function)
    e_c_q(t,i)		Unit cost of production
    e_c_e(t,e)		Price index for energy

*	Price variables
    e_p_q(t,i)		Price index (unit revenue) for activity
    e_p_f(t,i)		Price index of primary factor composite
    e_p_ee(t,i)		Price index of energy composite in production
    e_p_fe(t,i)		Price index of primary factor-energy composite
    e_p_c(t)		Price index of non-energy composite in consumption
    e_p_ec(t)		Price index of energy composit in consumption
    e_p_i(t)		Price index for investment
    e_p_e(t,i)		Price index of emission sources for combustion purpose

    e_p_cc(t)		Price index for aggregate consumption
    e_p_d(t,i)		Price of domestic good
    e_p_x(t,i)		Price of export good
    e_p_m(t,i)		Price of import good
    e_p_fxv		Exchange rate (intertemporal)
    e_p_fx		Exchange rate
    e_p_l(t)		Wage rate for producer
    e_p_le(t)		Wage rate for household
    e_p_w(t)		Price of utility
    e_p_u		Price of lifetime utility
    e_p_a(t,i)		Price of Armington good
    e_p_g(t)		Price of government expenditure
    e_p_k(t)		Shadow price of capital
    e_p_rk(t)		Rental price for industry
    e_p_rke(t)		Rental price for household
    e_p_gle		Price index for government lifetime consumption
    e_p_ca(t)		Price of emission permits
    e_p_ka(t)		Adjustmnet premium
    e_p_disc		Discount factor
    e_p_intr		Interest rate

*	Income balance
    e_v_inc_h		Extended lifetime income of household
    e_v_inc_g		Government income (period)
    e_v_inc_gl		Government income (lifetime)

    e_q_tcap		Post-terminal capital stock

    e_m_lum		Multiplier for lump-sum tax
    e_m_tl		Multiplier for labor tax
    e_m_tk		Multiplier for capital tax
    e_m_tli		Multiplier for labor income tax
    e_m_tci		Multiplier for capital income tax
    e_m_tc		Multiplier for consumption tax
    e_m_tm		Multiplier for import tariff

    e_t_pf(t)		Technology parameter for primary factor
    e_t_ene(t)		Technology parameter for energy
    e_roc_t_pf		Rate of change in t_pf
    e_roc_t_ene		Rate of change in t_ene

    e_v_gdp(t)		GDP
    e_c_emi(t)		Carbon emissions

    e_m_lum_p(t)	Lump-sum tax (period)
    e_share_g(t)	Share period government expenditure
;
*	------------------------------------------------------------------
*	Unit cost and price index.

*	Unit cost of consumption aggregation.
e_c_cc(t) ..

    c_cc(t) =e=
    ((sh_nene * p_c(t)**(1-sigcc) + (1-sh_nene) * p_ec(t)**(1-sigcc)
    )**(1/(1-sigcc)))$(sigcc ne 1)
    +
    ( p_c(t)**sh_nene * p_ec(t)**(1-sh_nene) )$(sigcc = 1);

*	Unit cost of Armington aggregation
e_c_a(t,i)$a0(i) ..

    c_a(t,i) =e=
    ( (1$(1-d0(i)) + p_d(t,i)$d0(i))**sh_ad(i)
    * (1$(1-m0(i)) + (p_m(t,i))$m0(i))**(1-sh_ad(i)) )$(siga(i) = 1)
    +
    ( ( (sh_ad(i) * p_d(t,i)**(1-siga(i)) )$d0(i)
    + ( (1-sh_ad(i)) * (p_m(t,i))**(1-siga(i)) )$m0(i)
    )**(1/(1-siga(i))) )$(siga(i) ne 1);

*	Unit cost of period utility
e_c_w(t) ..

    c_w(t) =e=
    ( (p_le(t))**(sh_lei)
    * ((1 + tc*m_tc) * p_cc(t) / pc00)**(1-sh_lei)
    )$(sigcl = 1)
    +
    ((sh_lei * p_le(t)**(1-sigcl)
    + (1-sh_lei) * ((1 + tc*m_tc) * p_cc(t)
    / pc00)**(1-sigcl))**(1/(1-sigcl))
    )$(sigcl ne 1);

*	Unit cost of lifetime utility
e_c_u ..

    c_u =e=
    ( prod(t, (p_w(t) / pref(t))**sh_u(t) ) )$(sigu = 1)
    +
    ( ( sum(t, sh_u(t) * (p_w(t) / pref(t))**(1-sigu) ) )**(1/(1-sigu))
    )$(sigu ne 1);

*	Unit cost of production.
e_c_q(t,i)$out0(i) ..

    c_q(t,i) =e=
    sum(nen$sh_nen(nen,i), sh_nen(nen,i) * p_a(t,nen))
    + sum(e$sh_nc(e,i), sh_nc(e,i) * p_a(t,e))
    + sum(cl$sh_cl(cl,i), sh_cl(cl,i) * p_e(t,cl))
    + sh_fe(i) * p_fe(t,i);

*	Price index of output.
e_p_q(t,i)$out0(i) ..

    p_q(t,i) =e=
    ( sh_x(i) * p_x(t,i)**(1+eta)
    + (1-sh_x(i)) * p_d(t,i)**(1+eta) )**(1/(1+eta));

*	Price index of primary factor composite.
e_p_f(t,i)$vpf0(i) ..

    p_f(t,i) =e=
    (
    ((1 + tl(i)*m_tl) * p_l(t) / pl0(i))**sh_fl(i)

    * ((1 + tk(i)*m_tk) * p_rk(t) / rk000(i))**(1-sh_fl(i))
    )$(sigpf(i) = 1)
    +
    (
    ( sh_fl(i)
    * ((1 + tl(i)*m_tl) * p_l(t) / pl0(i))**(1-sigpf(i))
    + (1-sh_fl(i))
    * ((1 + tk(i)*m_tk) * p_rk(t) / rk000(i))**(1-sigpf(i))
    )**(1/(1-sigpf(i))) )$(sigpf(i) ne 1);

*	Price index for energy composite in production.
e_p_ee(t,i)$veict(i) ..

    p_ee(t,i) =e=
    (( sum(ec$sh_e(ec,i), sh_e(ec,i) * p_e(t,ec)**(1-sigee))
    + sh_e("ele",i) * p_a(t,"ele")**(1-sigee))**(1/(1-sigee)) )$(sigee ne 1)
    +
    ( prod(ec$sh_e(ec,i), p_e(t,ec)**(sh_e(ec,i)))
    * p_a(t,"ele")**(sh_e("ele",i)) )$(sigee = 1);

*	Price index for primary factor-energy composite
e_p_fe(t,i)$vpfe0(i) ..

    p_fe(t,i) =e=
    ((p_f(t,i)/t_pf(t))**sh_f(i) * (p_ee(t,i)/t_ene(t))**(1-sh_f(i)))$(sigpfe = 1)
    +
    ((sh_f(i) * (p_f(t,i)/t_pf(t))**(1-sigpfe)
    + (1-sh_f(i)) * (p_ee(t,i)/t_ene(t))**(1-sigpfe))**(1/(1-sigpfe)))$(sigpfe ne 1);

*	Price index for non-energy composite in consumption.
e_p_c(t) ..

    p_c(t) =e=
    ( prod(nene$sh_c(nene),
    ((1+ts(nene))*p_a(t,nene) / pc0(nene))**sh_c(nene)) )$(sigc = 1)
    +
    (( sum(nene$sh_c(nene),
    sh_c(nene) * ((1+ts(nene))*p_a(t,nene) / pc0(nene))**(1-sigc))
    )**(1/(1-sigc)) )$(sigc ne 1);

*	Price index for energy composite in consumption.
e_p_ec(t) ..

    p_ec(t) =e=
    ( prod(ec$sh_ec(ec),
    ((1+ts(ec))*p_e(t,ec) / pc0(ec))**sh_ec(ec))
    * (((1+ts("ele"))*p_a(t,"ele") / pc0("ele"))**sh_ec("ele"))$sh_ec("ele")
    )$(sigec = 1)
    +
    (( sum(ec$sh_ec(ec), sh_ec(ec)
    * ((1+ts(ec))*p_e(t,ec) / pc0(ec))**(1-sigec))
    + (sh_ec("ele") * ((1+ts("ele"))*p_a(t,"ele") / pc0("ele"))**(1-sigec)
    )$sh_ec("ele")
    )**(1/(1-sigec)) )$(sigec ne 1);

*	Price index of emission sources.
e_c_e(t,e) ..

    c_e(t,e) =e= p_a(t,e) + carbt0(e) * p_ca(t) / vect(e);

*	Price of investment good.
e_p_i(t) ..

    p_i(t) * j0 =e= sum(i, p_a(t,i) * jd0(i));

*	Wage rate for household.
e_p_le(t) ..

    ple0 * p_le(t) =e= (1-tli*m_tli) * p_l(t);

*	Rental price for household.
e_p_rke(t) ..

    rke0 * p_rke(t) =e=
    (1-tci*m_tci) * rk0 * p_rk(t);

*	Adjustment premium.
e_p_ka(t)$phi ..

    p_ka(t) * pka0 =e= phi * (q_j(t)*j0/(q_k(t)*ks0))**2 * p_i(t) / 2;

*	------------------------------------------------------------------
*	Zero profit conditions

*	Production activity.
e_q_q(t,i)$out0(i) ..

    c_q(t,i) =g= (1-tq(i)) * p_q(t,i) / pq0(i);

*	Armington aggregation activity.
e_q_a(t,i)$a0(i) ..

    c_a(t,i) =g= p_a(t,i);

*	Export activity.
e_q_x(t,i)$x0(i) ..

    p_x(t,i) =g= p_fxv$capflo + p_fx(t)$bopcon;

*	Import activity.
e_q_m(t,i)$m0(i) ..

    ( (1+tm(i)*m_tm) * p_fxv  )$capflo
    +
    ( (1+tm(i)*m_tm) * p_fx(t) )$bopcon
    =g= pm0(i) * p_m(t,i);

*	Consumption aggregation activity.
e_q_cc(t) ..

    c_cc(t) =g= p_cc(t);

*	Period utility production activity.
e_q_w(t) ..

    p_disc(t) * c_w(t) =g= p_w(t);

*	Lifetime utility production activity.
e_q_u ..

    c_u =g= p_u;

*	Government expenditure (period)
e_q_g(t) ..

    sum(i, p_a(t,i) * gd0(i)) =g= p_g(t) * g0;

*	Investment good production.
e_q_j(t) ..

    p_i(t) * ( 1 + (phi * (q_j(t)*j0/(q_k(t)*ks0)))$phi ) =g= pk0 * p_k(t);

*	Capital accumulation (investment) activity.
e_q_k(t) ..

    (q_k(t) - scale_k)$tfirst(t)
    +
    (pk0 * p_k(t-1) / p_intr(t) -
    rke0 * p_rke(t) - (1-delta) * pk0 * p_k(t) - (p_ka(t) * pka0)$phi)$(not tfirst(t))
    =g= 0;

*	Carbon emissions.
e_q_cnc(t,e)$vect(e) ..

    c_e(t,e) =g= p_e(t,e);

*	------------------------------------------------------------------
*	Unit demand and supply

*	Export supply.
e_a_x(t,i)$x0(i) ..

    a_x(t,i) =e= x0(i) * (p_x(t,i) / p_q(t,i))**eta;

*	Domestic supply.
e_a_d(t,i)$d0(i) ..

    a_d(t,i) =e= d0(i) * (p_d(t,i) / p_q(t,i))**eta;

*	Labor demand.
e_a_l(t,i)$ld0(i) ..

    a_l(t,i) =e=
    ld0(i) * (p_f(t,i) / ((1 + tl(i)*m_tl)*p_l(t) / pl0(i)))**sigpf(i)
    * (p_fe(t,i) / p_f(t,i))**sigpfe * t_pf(t)**(sigpfe - 1);

*	Capital demand.
e_a_k(t,i)$kd0(i) ..

    a_k(t,i) =e=
    kd0(i) * (p_f(t,i) / ((1 + tk(i)*m_tk)*p_rk(t) / rk000(i)))**sigpf(i)
    * (p_fe(t,i) / p_f(t,i))**sigpfe * t_pf(t)**(sigpfe - 1);

*	Demand for energy (exc. cok and lim) in production.
e_a_e(t,i,ii)$(ene(i) and veic(i,ii)) ..

    a_e(t,i,ii) =e=
    (veic(i,ii) * (p_ee(t,ii) / p_e(t,i))**sigee
    * (p_fe(t,ii) / p_ee(t,ii))**sigpfe * t_ene(t)**(sigpfe - 1))$ec(i)
    +
    (veic(i,ii) * (p_ee(t,ii) / p_a(t,i))**sigee
    * (p_fe(t,ii) / p_ee(t,ii))**sigpfe * t_ene(t)**(sigpfe - 1))$ele(i);

*	Consumption demand for energy goods
e_a_ec(t,i)$(ene(i) and cd0(i))..

    a_ec(t,i) =e=
    (cd0(i) * (p_ec(t) / ((1+ts(i))*p_e(t,i) / pc0(i)))**sigec
    * (p_cc(t) / p_ec(t))**sigcc)$ec(i)
    +
    (cd0(i) * (p_ec(t) / ((1+ts(i))*p_a(t,i) / pc0(i)))**sigec
    * (p_cc(t) / p_ec(t))**sigcc)$ele(i);

*	Consumption demand for non-energy goods
e_a_c(t,i)$(nene(i) and cd0(i)) ..

    a_c(t,i) =e=
    ( (cd0(i) * (p_c(t) / ((1+ts(i))*p_a(t,i) / pc0(i)))**sigc
    * (p_cc(t) / p_c(t))**sigcc) )$nen(i)
    +
    ( (cd0(i) * (p_c(t) / ((1+ts(i))*p_e(t,i) / pc0(i)))**sigc
    * (p_cc(t) / p_c(t))**sigcc) )$cl(i);

*	Demand for aggregate consumption.
e_a_cc(t) ..

    a_cc(t)
    =e=
    c0 * (c_w(t) / ((1+tc*m_tc) * p_cc(t) / pc00))**sigcl;

*	Demand for leisure.
e_a_lei(t) ..

    a_lei(t) =e= lei0 * (c_w(t) / p_le(t))**sigcl;

*	Demand for period utility.
e_a_u(t) ..

    a_u(t) =e= wref(t) * (c_u / (p_w(t) / pref(t)))**sigu;

*	Demand for domestic good from Armington activity.
e_a_ad(t,i)$d0(i) ..

    a_ad(t,i) =e= d0(i) * (c_a(t,i) / p_d(t,i))**siga(i);

*	Demand for imported good from Armington activity.
e_a_am(t,i)$m0(i) ..

    a_am(t,i) =e= m0(i) * (c_a(t,i) / p_m(t,i))**siga(i);

*	Labor supply.
e_q_sl(t) ..

    li0 * q_sl(t) =e= scale_l * lref(t) - a_lei(t) * q_w(t);

*	------------------------------------------------------------------
*	Market clearing conditions.

*	Market for aggregate consumption.
e_p_cc(t) ..

    c0 * q_cc(t) =g= a_cc(t) * q_w(t);

*	Market for domestic goods.
e_p_d(t,i)$d0(i) ..

    a_d(t,i) * q_q(t,i) =g= a_ad(t,i) * q_a(t,i);

*	Market for export goods
e_p_x(t,i)$x0(i) ..

    a_x(t,i) * q_q(t,i) =g= x0(i) * q_x(t,i);

*	Market for import goods.
e_p_m(t,i)$m0(i) ..

    m0(i) * q_m(t,i) =g= a_am(t,i) * q_a(t,i);

*	Exchange rate (intertemporal).
e_p_fxv$capflo ..

    sum(t, sum(i$x0(i), pref(t) * x0(i) * q_x(t,i)))
    =g= sum(t, sum(i$m0(i), pref(t) * m0(i) * q_m(t,i)));

*	Exchange rate (non-intertemporal).
e_p_fx(t)$bopcon ..

    sum(i$x0(i), x0(i) * q_x(t,i)) =g= sum(i$m0(i), m0(i) * q_m(t,i));

*	Market for labor.
e_p_l(t) ..

    li0 * q_sl(t) =g= sum(i$ld0(i), a_l(t,i) * q_q(t,i));

*	Market for period utility.
e_p_w(t) ..

    wg0 * q_w(t) =g= a_u(t) * q_u;

*	Market for lifetime utility.
e_p_u ..

    v_inc_h =g= ubar * q_u * p_u;

*	Market for Armington goods.
e_p_a(t,i) ..

    a0(i) * q_a(t,i) =g=
    ( (a_c(t,i) * q_cc(t))$cd0(i)
    + sum(ii$int0(i,ii), int0(i,ii) * q_q(t,ii)) )$nen(i)

    + ( a_ec(t,i) * q_cc(t)
    + sum(ii$veic(i,ii), a_e(t,i,ii) * q_q(t,ii)) )$ele(i)

    + ( vect(i) * q_cnc(t,i)
    + sum(ii$veinc(i,ii), veinc(i,ii) * q_q(t,ii)) )$cl(i)

    + ( vect(i) * q_cnc(t,i)
    + sum(ii$veinc(i,ii), veinc(i,ii) * q_q(t,ii)) )$ec(i)

    + sh_inv(i) * j0 * q_j(t)

    + (sh_inv(i) * phi * ((q_j(t)*j0)**2 /(2*q_k(t)*ks0)))$phi

    + gd0(i) * q_g(t);

*	Government expenditure (period).
e_p_g(t) ..

    (g0 * q_g(t) - scale_g * gqref(t))$(not dcalib)
    + (v_inc_g(t) - g0 * p_g(t) * q_g(t))$dcalib =g= 0;

*	Shadow price of capital.
e_p_k(t) ..

    j0 * q_j(t) + (1-delta) * ks0 * q_k(t)
    =g= ks0 * (q_k(t+1)$(not tlast(t)) + q_tcap$tlast(t));

*	Market for renting capital.
e_p_rk(t) ..

    ks0 * q_k(t) =g= sum(i$kd0(i), a_k(t,i) * q_q(t,i));

*	Price index for government lifetime consumption.
e_p_gle$(not dcalib) ..

    pgle0 * p_gle =e= sum(t, share_g(t) * p_disc(t) * p_g(t));

*     p_gle =e= sum(t, sh_g(t) * p_disc(t) * p_g(t) / pref(t));

*	Share of period government expenditure:
e_share_g(t)$(not dcalib) ..

    share_g(t) =e= p_disc(t) * g0 * p_g(t) * q_g(t)
    / sum(tt, p_disc(tt) * g0 * p_g(tt) * q_g(tt));

*	Emission demand:
e_c_emi(t) ..

    c_emi(t) =e= sum(e, carbt0(e) * q_cnc(t,e));

*	Market for emission permit.
e_p_ca(t)$(not ca_tax) ..

    carbl(t) =g= c_emi(t);

*	Price index of emission sources.
e_p_e(t,i)$e(i) ..

    vect(i) * q_cnc(t,i) =g=
    ((a_c(t,i) * q_cc(t))$cd0(i)
    + sum(ii$veic(i,ii), veic(i,ii) * q_q(t,ii)))$cl(i)
    +
    ((a_ec(t,i) * q_cc(t))$cd0(i)
    + sum(ii$veic(i,ii), a_e(t,i,ii) * q_q(t,ii)))$ec(i);

*	Interest rate.
e_p_intr(t)$(not tfirst(t)) ..

    p_intr(t) =e=
    (p_w(t) / p_w(t-1))$bopcon
*     (1/(1+r - p_w(t)/p_w(t-1)))$bopcon
    + (1/(1+r))$capflo;

*	Discount factor.
e_p_disc(t) ..

    p_disc(t) =e=
    (prod(tt$((not tfirst(tt)) and (ord(tt) le ord(t))),
    p_intr(tt)))$(not tfirst(t))
    + 1$tfirst(t);

*	------------------------------------------------------------------
*	Income constraints.

*	Household lifetime income.
e_v_inc_h ..

    v_inc_h =e=

    sum(t, p_disc(t) * p_le(t) * ple0 * scale_l * lref(t))

    + sum(t$tfirst(t), (rke0 * p_rke(t) + (1-delta) * pk0 * p_k(t)
    + (pka0 * p_ka(t))$phi)) * ks0 * scale_k

    - sum(tlast, p_disc(tlast) * pk0 * p_k(tlast) *  ks0 * q_tcap)

    + sum(t, p_disc(t) * p_g(t) * vot * scale_t * qref(t))

    - (pgle0 * p_gle * m_lum)$(not dcalib)

    - sum(t, p_disc(t) * p_g(t) * m_lum_p(t))
;
*	Government income (period).
e_v_inc_g(t) ..

    v_inc_g(t) =e=

    sum(i$ld0(i),
    tl(i) * m_tl * p_l(t) * a_l(t,i) * q_q(t,i))

    + sum(i$kd0(i),
    tk(i) * m_tk * rk0 * p_rk(t) * a_k(t,i) * q_q(t,i))

    + sum(i$out0(i), tq(i) * p_q(t,i) * out0(i) * q_q(t,i))

    + sum(i$m0(i), tm(i) * m_tm
    * (p_fxv$capflo + p_fx(t)$bopcon) * m0(i) * q_m(t,i))

    + tli * m_tli * p_l(t) * li0 * q_sl(t)

    + tci * m_tci * rk0 * p_rk(t) * ks0 * q_k(t)

    + tc * m_tc * p_cc(t) * a_cc(t) * q_w(t)

    + ( sum(nen$cd0(nen), ts(nen) * p_a(t,nen) * a_c(t,nen))
    + sum(ele$cd0(ele), ts(ele) * p_a(t,ele) * a_ec(t,ele))
    + sum(cl$cd0(cl), ts(cl) * p_a(t,cl) * a_c(t,cl))
    + sum(ec$cd0(ec), ts(ec) * p_a(t,ec) * a_ec(t,ec)) ) * q_cc(t)

    - p_g(t) * vot * scale_t * qref(t)

    + p_ca(t) * carbl(t)

    + p_g(t) * m_lum_p(t)
;
*	Government income (lifetime).
e_v_inc_gl ..

    v_inc_gl =e= sum(t, p_disc(t) * v_inc_g(t))

    + (pgle0 * p_gle * m_lum)$(not dcalib)
;
*	------------------------------------------------------------------
*	Terminal condition.
e_q_tcap ..

    sum(t$tlast(t), q_j(t) * q_w(t-1) - q_j(t-1) * q_w(t)) =e= 0;

*	------------------------------------------------------------------
*	Constraints

e_m_lum$((not dcalib) and tins("lum")) ..

    v_inc_gl =e= sum(t, p_disc(t) * p_g(t) * g0 * q_g(t));

e_m_tl$tins("lab") ..

    v_inc_gl =e= sum(t, p_disc(t) * p_g(t) * g0 * q_g(t));

e_m_tk$tins("cap") ..

     v_inc_gl =e= sum(t, p_disc(t) * p_g(t) * g0 * q_g(t));

e_m_tli$tins("lin") ..

    v_inc_gl =e= sum(t, p_disc(t) * p_g(t) * g0 * q_g(t));

e_m_tci$tins("cin") ..

    v_inc_gl =e= sum(t, p_disc(t) * p_g(t) * g0 * q_g(t));

e_m_tc$tins("ctx") ..

    v_inc_gl =e= sum(t, p_disc(t) * p_g(t) * g0 * q_g(t));

e_m_tm$tins("imp") ..

    v_inc_gl =e= sum(t, p_disc(t) * p_g(t) * g0 * q_g(t));

*	----------------------------------------------------------------------
*	Technology parameters:

e_v_gdp(t) ..

    v_gdp(t)
    =e= (p_cc(t) * a_cc(t) * q_w(t)
    - (sum(nen$cd0(nen), ts(nen) * p_a(t,nen) * a_c(t,nen))
    + sum(ele$cd0(ele), ts(ele) * p_a(t,ele) * a_ec(t,ele))
    + sum(cl$cd0(cl), ts(cl) * p_a(t,cl) * a_c(t,cl))
    + sum(ec$cd0(ec), ts(ec) * p_a(t,ec) * a_ec(t,ec))) * q_cc(t)
    + p_g(t) * g0 * q_g(t)
    + p_i(t) * j0 * q_j(t) * (1 + phi * j0 * q_j(t) / (2*q_k(t)*ks0))
    + sum(i, (p_fxv$capflo + p_fx(t)$bopcon) * x0(i) * q_x(t,i))
    - sum(i$m0(i),
    (1+tm(i))*(p_fxv$capflo + p_fx(t)$bopcon) * m0(i) * q_m(t,i))) / p_cc(t);
    ;

e_t_pf(t)$dcalib ..

    t_pf(t) =e= 1$tfirst(t)
    +  (((1 + roc_t_pf) * t_pf(t-1)))$(not tfirst(t));

e_t_ene(t)$dcalib ..

    t_ene(t) =e= 1$tfirst(t)
    +  (((1 + roc_t_ene) * t_ene(t-1)))$(not tfirst(t));

e_roc_t_pf$dcalib ..

    (v_gdp("38")/v_gdp("1"))**(1/37) - 1 =e= roc_gdp;

e_roc_t_ene$dcalib ..

    (c_emi("38")/c_emi("1"))**(1/37) - 1 =e= roc_ca;

*	----------------------------------------------------------------------

e_m_lum_p(t)$dcalib ..

    p_g(t) * g0 * q_g(t) / (p_cc(t) * v_gdp(t)) =e= sh_g_gdp(t);

*	------------------------------------------------------------------
*	Model declaration.
model dynamic_mcp Dynamic model in MCP format /

      e_a_x.a_x, e_a_d.a_d, e_a_l.a_l, e_a_k.a_k, e_a_e.a_e, e_a_ec.a_ec,
      e_a_c.a_c, e_a_cc.a_cc, e_a_lei.a_lei, e_a_u.a_u, e_a_ad.a_ad,
      e_a_am.a_am,

      e_q_q.q_q, e_q_a.q_a, e_q_x.q_x, e_q_m.q_m, e_q_cc.q_cc, e_q_w.q_w,
      e_q_u.q_u, e_q_g.q_g, e_q_j.q_j, e_q_sl.q_sl, e_q_k.q_k,
*       e_q_gle.q_gle,
      e_q_cnc.q_cnc,

      e_c_cc.c_cc, e_c_a.c_a, e_c_w.c_w, e_c_u.c_u, e_c_q.c_q, e_c_e.c_e,

      e_p_q.p_q, e_p_f.p_f, e_p_ee.p_ee, e_p_fe.p_fe, e_p_c.p_c, e_p_ec.p_ec,
      e_p_e.p_e, e_p_k.p_k, e_p_rke.p_rke, e_p_rk.p_rk, e_p_gle.p_gle,
      e_p_i.p_i, e_p_ka.p_ka, e_p_intr.p_intr, e_p_disc.p_disc,

      e_p_cc.p_cc, e_p_d.p_d, e_p_x.p_x, e_p_m.p_m, e_p_fxv.p_fxv,
      e_p_fx.p_fx, e_p_l.p_l, e_p_le.p_le, e_p_w.p_w, e_p_u.p_u, e_p_a.p_a,
      e_p_g.p_g, e_p_ca.p_ca,

      e_v_inc_h.v_inc_h, e_v_inc_g.v_inc_g, e_v_inc_gl.v_inc_gl,
      e_q_tcap.q_tcap,

      e_m_lum.m_lum,
      e_m_tl.m_tl, e_m_tk.m_tk, e_m_tli.m_tli, e_m_tci.m_tci,
      e_m_tc.m_tc, e_m_tm.m_tm,

      e_t_pf.t_pf, e_t_ene.t_ene, e_roc_t_pf.roc_t_pf, e_roc_t_ene.roc_t_ene,
      e_v_gdp.v_gdp, e_c_emi.c_emi,

      e_m_lum_p.m_lum_p,

      e_share_g.share_g

      /;
*	------------------------------------------------------------------
*	Lower bounds for variables.

*	Unit demand.
a_x.lo(t,i) = 0;
a_d.lo(t,i) = 0;
a_l.lo(t,i) = 0;
a_k.lo(t,i) = 0;
a_e.lo(t,ii,i) = 0;
a_ec.lo(t,i) = 0;
a_c.lo(t,i) = 0;
a_cc.lo(t) = 0;
a_lei.lo(t) = 0;
a_u.lo(t) = 0;
a_ad.lo(t,i) = 0;
a_am.lo(t,i) = 0;

*	Activity.
q_q.lo(t,i) = 0;
q_a.lo(t,i) = 0;
q_x.lo(t,i) = 0;
q_m.lo(t,i) = 0;
q_cc.lo(t) = 0;
q_w.lo(t) = 0;
q_g.lo(t) = 0;
q_j.lo(t) = 0;
q_sl.lo(t) = 0;
q_u.lo = 0;
q_k.lo(t) = 0;
q_cnc.lo(t,e) = 0;
q_tcap.lo = 0;

*	Unit cost.
c_cc.lo(t) = 1e-6;
c_a.lo(t,i) = 1e-6;
c_w.lo(t) = 1e-6;
c_u.lo = 1e-6;
c_q.lo(t,i) = 1e-6;

*	Price.
p_q.lo(t,i) = 1e-6;
p_f.lo(t,i) = 1e-6;
p_ee.lo(t,i) = 1e-6;
p_fe.lo(t,i) = 1e-6;
p_c.lo(t) = 1e-6;
p_ec.lo(t) = 1e-6;
p_e.lo(t,e) = 1e-6;
p_cc.lo(t) = 1e-6;
p_d.lo(t,i) = 1e-6;
p_x.lo(t,i) = 1e-6;
p_m.lo(t,i) = 1e-6;
p_fxv.lo = 1e-6;
p_fx.lo(t) = 1e-6;
p_l.lo(t) = 1e-6;
p_le.lo(t) = 1e-6;
p_w.lo(t) = 1e-6;
p_u.lo = 1e-6;
p_a.lo(t,i) = 1e-6;
p_g.lo(t) = 1e-6;
p_k.lo(t) = 1e-6;
p_rk.lo(t) = 1e-6;
p_rke.lo(t) = 1e-6;
p_gle.lo = 1e-6;
p_ca.lo(t) = 0;
p_ka.lo(t) = 1e-6;
p_i.lo(t) = 1e-6;

*	Income
v_inc_h.lo = 0;
v_inc_g.lo(t) = 0;
v_inc_gl.lo = 0;

v_gdp.lo(t) = 0;
c_emi.lo(t) = 0;
t_pf.lo(t) = 0;
t_ene.lo(t) = 0;
roc_t_pf.lo = 0;
roc_t_ene.lo = 0;

share_g.lo(t) = 0;

*	------------------------------------------------------------------
*	Initial values for variables.

*	Unit demand.
a_x.l(t,i) = x0(i);
a_d.l(t,i) = d0(i);
a_l.l(t,i) = ld0(i);
a_k.l(t,i) = kd0(i);
a_e.l(t,ene,i) = veic(ene,i);
a_ec.l(t,ene) = cd0(ene);
a_c.l(t,nene) = cd0(nene);
a_cc.l(t) = c0;
a_lei.l(t) = lei0;
a_u.l(t) = wref(t);
a_ad.l(t,i) = d0(i);
a_am.l(t,i) = m0(i);

*	Activity.
q_q.l(t,i) = qref(t);
q_a.l(t,i) = qref(t);
q_x.l(t,i) = qref(t);
q_m.l(t,i) = qref(t);
q_cc.l(t) = qref(t);
q_w.l(t) = qref(t);
q_g.l(t) = qref(t);
q_j.l(t) = qref(t);
q_sl.l(t) = qref(t);
q_u.l = 1;
q_k.l(t) = qref(t);
q_cnc.l(t,e) = qref(t);
q_tcap.l = sum(tlast, qref(tlast)) * (1+gr);

*	Unit cost.
c_cc.l(t) = 1;
c_a.l(t,i) = 1;
c_w.l(t) = 1;
c_u.l = 1;
c_q.l(t,i) = 1;
c_e.l(t,i) = 1;

*	Price.
p_q.l(t,i) = 1;
p_f.l(t,i) = 1;
p_ee.l(t,i) = 1;
p_fe.l(t,i) = 1;
p_c.l(t) = 1;
p_ec.l(t) = 1;
p_e.l(t,e) = 1;
p_cc.l(t) = 1;
p_d.l(t,i) = 1;
p_x.l(t,i) = 1;
p_m.l(t,i) = 1;
p_fxv.l = 1;
p_fx.l(t) = 1;
p_l.l(t) = 1;
p_le.l(t) = 1;
p_u.l = 1;
p_a.l(t,i) = 1;
p_g.l(t) = 1;
p_k.l(t) = 1;
p_rk.l(t) = 1;
p_rke.l(t) = 1;
p_gle.l = 1;
p_ca.l(t) = 0;
p_ka.l(t) = 1;
p_i.l(t) = 1;
p_w.l(t) = pref(t);

*	Income.
v_inc_h.l = ubar;
v_inc_g.l(t) = gqref(t);
v_inc_gl.l = gbar;

*	Tax multiplier.
m_lum.l = 0;
m_tl.fx = 1;
m_tk.fx = 1;
m_tli.fx = 1;
m_tci.fx = 1;
m_tc.fx = 1;
m_tm.fx = 1;

*	Technology parameters
t_pf.fx(t) = 1;
t_ene.fx(t) = 1;
roc_t_pf.fx = 0;
roc_t_ene.fx = 0;

v_gdp.l(t)
    = (p_cc.l(t) * a_cc.l(t) * q_w.l(t)
    - (sum(nen$cd0(nen), ts(nen) * p_a.l(t,nen) * a_c.l(t,nen))
    + sum(ele$cd0(ele), ts(ele) * p_a.l(t,ele) * a_ec.l(t,ele))
    + sum(cl$cd0(cl), ts(cl) * p_a.l(t,cl) * a_c.l(t,cl))
    + sum(ec$cd0(ec), ts(ec) * p_a.l(t,ec) * a_ec.l(t,ec))) * q_cc.l(t)
    + p_g.l(t) * g0 * q_g.l(t)
    + p_i.l(t) * j0 * q_j.l(t) * (1 + phi * j0 * q_j.l(t) / (2*q_k.l(t)*ks0))
    + sum(i, (p_fxv.l$capflo + p_fx.l(t)$bopcon) * x0(i) * q_x.l(t,i))
    - sum(i$m0(i),  (1+tm(i))*(p_fxv.l$capflo + p_fx.l(t)$bopcon) * m0(i) * q_m.l(t,i))
    )/p_cc.l(t);
    ;

c_emi.l(t) = sum(e, carbt0(e) * q_cnc.l(t,e));

m_lum_p.fx(t) = 0;

*	Numeraire (fix household lifetime income)
v_inc_h.fx = v_inc_h.l;
display v_inc_h.l;
parameter
    chk_numeraire	Numeraire check;

*	Interest rate.
p_intr.lo(t) = 1e-6;
p_intr.l(t) = 1 / (1 + r);
display p_intr.l;

*	Discount factor.
p_disc.lo(t) = 1e-6;
p_disc.l(t) =
    (prod(tt$((not tfirst(tt)) and (ord(tt) le ord(t))),
    p_intr.l(tt)))$(not tfirst(t))
    +
    1$tfirst(t);
display p_disc.l;

share_g.l(t)
    = p_disc.l(t) * g0 * p_g.l(t) * q_g.l(t)
    / sum(tt, p_disc.l(tt) * g0 * p_g.l(tt) * q_g.l(tt));
display share_g.l;

dynamic_mcp.workspace = 500;

$ontext
*	------------------------------------------------------------------
*	Benchmark replication with MCP
display "com: Benchmark replication with MCP";
dynamic_mcp.workspace = 500;
dynamic_mcp.iterlim = 0;
dynamic_mcp.optfile=1;
solve dynamic_mcp using mcp;
chk_numeraire = v_inc_h.m;
display chk_numeraire;

*	------------------------------------------------------------------
*	Cleanup calculation with MCP
display "com: Cleanup calculation with MCP";
dynamic_mcp.iterlim = 9000;
solve dynamic_mcp using mcp;
chk_numeraire = v_inc_h.m;
display chk_numeraire;
$offtext


* --------------------
* Local Variables:
* mode: gams
* fill-column: 78
* End: