CJL1

$ontext
annual version of CJL

Modifications:
(1) annual version when time interval of one period becomes a parameter
(2) get rid of CO2 time travel in Forcing
(3) use explicit finite difference scheme in temperature in atmosphere
(4) cancel fixed saving rate (s.fx(t) = .22)
(5) cancel extra variables and equations and simplify
(6) cancel terminal condition of capital
(7) use nonzero terminal value function
(8) cancel the fixed emission control rate for the first period
(9) change sum(t$(ord(t)<card(t)), E(t)) <= 6000 into sum(t, E(t)) <= 6000 (10) change participate rate at the first stage from 0.25372 to 1 (11) use the solution of CJL2 as initial guess (12) change scale1 from 194 to 8600 $offtext SETS T Time periods /1*600/ ; parameter deltat interval length of each period (number of years); deltat = 1; SCALARS ** Preferences B_ELASMU Elasticity of marginal utility of consumption / 2.0 / B_PRSTP Initial rate of social time preference per year / .015 / ** Population and technology POP0 2005 world population millions /6514 / GPOP0 Growth rate of population per year /.035 / POPASYM Asymptotic population / 8600 / A0 Initial level of total factor productivity /.02722 / GA0 Initial growth rate for technology per year /.0092 / DELA Decline rate of technol change per year /.001 / DK Depreciation rate on capital per year /.100 / GAMA Capital elasticity in production function /.300 / K0 2005 value capital trill 2005 US dollars /137. / ** Emissions SIG0 CO2-equivalent emissions-GNP ratio 2005 /.13418 / GSIGMA Initial growth of sigma per year /-.00730 / DSIG Decline rate of decarbonization per year /.003 / ELAND0 Carbon emissions from land 2005(GtC per year) / 1.1000 / ** Carbon cycle MAT2000 Concentration in atmosphere 2005 (GtC) /808.9 / MU2000 Concentration in upper strata 2005 (GtC) /1255 / ML2000 Concentration in lower strata 2005 (GtC) /18365 / b12 Carbon cycle transition matrix per year /0.0189288 / b23 Carbon cycle transition matrix per year /0.005 / ** Climate model T2XCO2 Equilibrium temp impact of CO2 doubling oC / 3 / FEX0 Estimate of 2000 forcings of non-CO2 GHG / -.06 / FEX1 Estimate of 2100 forcings of non-CO2 GHG / 0.30 / TOCEAN0 2000 lower strat. temp change (C) from 1900 /.0068 / TATM0 2000 atmospheric temp change (C)from 1900 /.7307 / C1 Climate-equation coefficient for upper level per year /.0220 / C3 Transfer coeffic upper to lower stratum per year /.300 / C4 Transfer coeffic for lower level per year /.0050 / FCO22X Estimated forcings of equilibrium co2 doubling /3.8 / ** Climate damage parameters calibrated for quadratic at 2.5 C for 2105 A1 Damage intercept / 0.00000 / A2 Damage quadratic term / 0.0028388 / A3 Damage exponent / 2.00 / ** Abatement cost EXPCOST2 Exponent of control cost function /2.8 / PBACK Cost of backstop 2005 000$ per tC 2005 /1.17 / BACKRAT Ratio initial to final backstop cost / 2 / GBACK Initial cost decline backstop pc per year /.005 / LIMMIU Upper limit on control rate / 1 / ** Availability of fossil fuels FOSSLIM Maximum cumulative extraction fossil fuels / 6000 / ; SET TLAST(T); TLAST(T) = YES$(ORD(T) EQ CARD(T)); parameters bb11 Carbon cycle transition matrix bb12 Carbon cycle transition matrix bb21 Carbon cycle transition matrix bb22 Carbon cycle transition matrix bb23 Carbon cycle transition matrix bb32 Carbon cycle transition matrix bb33 Carbon cycle transition matrix CC1 Climate-equation coefficient for upper level CC3 Transfer coeffic upper to lower stratum CC4 Transfer coeffic for lower level GGPOP0 Growth rate of population GGA0 Initial growth rate for technology DDELA Decline rate of technol change GGSIGMA Initial growth of sigma DDSIG Decline rate of decarbonization EELAND0 Carbon emissions from land 2005(GtC) GGBACK Initial cost decline backstop pc scale1 Scaling coefficient in the objective function; bb12 = b12 * deltat; bb23 = b23 * deltat; bb11 = 1 – bb12; bb21 = 587.473*bb12/1143.894; bb22 = 1 – bb21 – bb23; bb32 = 1143.894*bb23/18340; bb33 = 1 – bb32 ; CC1 = C1*deltat; CC4 = C4*deltat; CC3 = C3; GGPOP0 = GPOP0*deltat; GGA0 = GA0*deltat; DDELA = DELA*deltat; GGSIGMA = GSIGMA*deltat; DDSIG = DSIG*deltat; EELAND0 = ELAND0*deltat; GGBACK = GBACK*deltat; scale1 = POPASYM; *scale1 = 194; PARAMETERS L(T) Level of population and labor AL(T) Level of total factor productivity SIGMA(T) CO2-equivalent-emissions output ratio RR(T) Average utility social discount factor beta one-period discount factor GA(T) Growth rate of productivity from 0 to T FORCOTH(T) Exogenous forcing for other greenhouse gases GCOST1 Growth of cost factor GSIG(T) Cumulative improvement of energy efficiency ETREE(T) Emissions from deforestation COST1(t) Adjusted cost for backstop LAM Climate model parameter Gfacpop(T) Growth factor population ; * Important parameters for the model LAM = FCO22X/ T2XCO2; Gfacpop(T) = (exp(GGPOP0*(ORD(T)-1))-1)/exp(GGPOP0*(ORD(T)-1)); L(T)=POP0* (1- Gfacpop(T))+Gfacpop(T)*popasym; ga(T)=GGA0*EXP(-DDELA*(ORD(T)-1)); al(“1”) = a0; LOOP(T, al(T+1)=al(T)/((1-ga(T)));); gsig(T)=GGSIGMA*EXP(-DDSIG*(ORD(T)-1));sigma(“1”)=sig0;LOOP(T,sigma(T+1)=(sigma(T)/((1-gsig(T+1))));); cost1(T) = (PBACK*SIGMA(T)/EXPCOST2)* ( (BACKRAT-1+ EXP (-GGBACK* (ORD(T)-1) ) )/BACKRAT); ETREE(T) = EELAND0*(1-0.01*deltat)**(ord(T)-1); RR(t)=1/((1+B_PRSTP)**(deltat*(ord(T)-1))); beta = 1 / ((1+B_PRSTP)**deltat); FORCOTH(T)= FEX0+ 0.01*deltat*(FEX1-FEX0)*(ORD(T)-1)$((ORD(T)-1)*deltat<=100)+ 0.36$((ORD(T)-1)*deltat>100);

***********************************
* get terminal value function
***********************************

set dim indices of dimensions /1*6/;

parameters
SVmin(dim) minimal state variables
SVmax(dim) maximal state variables
dz(dim) derivative of z over x
;

set d index /1*5/;
alias(d,d2,d3,d4,d5,d6);

parameter coefs(d,d2,d3,d4,d5,d6) complete Chebyshev coefficients;

execute_load ‘coefs_Climate_termValue’,SVmin, SVmax, dz, coefs;

*****************************************************************

VARIABLES
K(T) Capital stock trillions US dollars
MAT(T) Carbon concentration in atmosphere GtC
MU(T) Carbon concentration in shallow oceans Gtc
ML(T) Carbon concentration in lower oceans GtC
TATM(T) Temperature of atmosphere in degrees C
TOCEAN(T) Temperature of lower oceans degrees C
E(T) CO2-equivalent emissions GtC
C(T) Consumption trillions US dollars
MIU(T) Emission control rate GHGs
UTILITY;

POSITIVE VARIABLES K, MAT, MU, ML, TATM, TOCEAN, C, MIU, E;

free variable zPlus(dim) normalized tomorrow state variables;
zPlus.lo(dim) = -1;
zPlus.up(dim) = 1;
zPlus.l(dim) = 0;

free variable ChebyBasis(dim,d) Chebyshev polynomial basis;

variable SVend(dim) terminal state variables;
variable termUtil terminal utility function

*****************************************************************

EQUATIONS

KK0 Initial condition for capital
MMAT0 Starting atmospheric concentration
MMU0 Initial shallow ocean concentration
MML0 Initial lower ocean concentration
TATM0EQ Initial condition for atmospheric temperature
TOCEAN0EQ Initial condition for lower ocean temperature

KK(T) Capital balance equation
MMAT(T) Atmospheric concentration equation
MMU(T) Shallow ocean concentration
MML(T) Lower ocean concentration
TATMEQ(T) Temperature-climate equation for atmosphere
TOCEANEQ(T) Temperature-climate equation for lower oceans
EE(T) Emissions equation

CCACCA total cumulative carbon emissions GTX
UTIL Objective function

KendC Terminal capital
CAendC Terminal atmospheric carbon concentration
CSendC Terminal shallow ocean concentration
CDendC Terminal lower ocean concentration
TAendC Terminal atmospheric temperature
TLendC Terminal lower ocean temperature
normalizeSVend(dim) normalized state variables
ChebyBasisFunDeg0(dim) Chebyshev polynomial basis with degree 0
ChebyBasisFunDeg1(dim) Chebyshev polynomial basis with degree 1
recursiveCheby(dim,d) recursive formula for Chebyshev polynomial basis
TerminalUtilEQ terminal utility equation
;

** Equations of the model

KK0.. K(‘1’) =E= K0;
MMAT0.. MAT(‘1’) =E= MAT2000;
MMU0.. MU(‘1’) =E= MU2000;
MML0.. ML(‘1’) =E= ML2000;
TATM0EQ.. TATM(‘1’) =E= TATM0;
TOCEAN0EQ.. TOCEAN(‘1’) =E= TOCEAN0;

KK(T).. K(T+1) =L= (1-DK)**deltat *K(T)+deltat*( ( 1 – cost1(T)*(MIU(T)**EXPcost2) )*
(AL(T)*L(T)**(1-GAMA)*K(T)**GAMA) / (1+A1*TATM(T)+ A2*TATM(T)**A3) – C(T) );
MMAT(T+1).. MAT(T+1) =E= MAT(T)*bb11+MU(T)*bb21 + E(T);
MML(T+1).. ML(T+1) =E= ML(T)*bb33+bb23*MU(T);
MMU(T+1).. MU(T+1) =E= MAT(T)*bb12+MU(T)*bb22+ML(T)*bb32;
TATMEQ(T+1).. TATM(T+1) =E= TATM(t)+CC1*(FCO22X*(log(MAT(T)/596.4)/log(2))+FORCOTH(T) – LAM*TATM(t)-CC3*(TATM(t)-TOCEAN(t)));
TOCEANEQ(T+1).. TOCEAN(T+1) =E= TOCEAN(T)+CC4*(TATM(T)-TOCEAN(T));
EE(T).. E(T) =E= deltat*SIGMA(T)*(1-MIU(T))* (AL(T)*L(T)**(1-GAMA)*K(T)**GAMA) + ETREE(T);

CCACCA.. sum(T, E(T)) =L= FOSSLIM;
UTIL.. UTILITY =E= sum( T, (((C(T)/L(T))**(1-B_ELASMU)-1)/(1-B_ELASMU)) * (RR(T)*L(T)*deltat/scale1) ) + termUtil* (beta*sum(TLAST,RR(TLAST))*POPASYM/scale1) ;

KendC.. SVend(‘1’) =e= sum(TLAST, (1-DK)**deltat*K(TLAST)+deltat*( ( 1 – cost1(TLAST)*(MIU(TLAST)**EXPcost2) )*
(AL(TLAST)*L(TLAST)**(1-GAMA)*K(TLAST)**GAMA) / (1+A1*TATM(TLAST)+ A2*TATM(TLAST)**A3) – C(TLAST) ));
CAendC.. SVend(‘2’) =e= sum(TLAST, MAT(TLAST)*bb11+MU(TLAST)*bb21 + E(TLAST));
CSendC.. SVend(‘3’) =e= sum(TLAST, MAT(TLAST)*bb12+MU(TLAST)*bb22+ML(TLAST)*bb32);
CDendC.. SVend(‘4’) =e= sum(TLAST, ML(TLAST)*bb33+bb23*MU(TLAST));
TAendC.. SVend(‘5’) =e= sum(TLAST, TATM(TLAST)+CC1*( FCO22X*(log(MAT(TLAST)/596.4)/log(2))+FORCOTH(TLAST) – LAM*TATM(TLAST)-CC3*(TATM(TLAST)-TOCEAN(TLAST))));
TLendC.. SVend(‘6’) =e= sum(TLAST, TOCEAN(TLAST)+CC4*(TATM(TLAST)-TOCEAN(TLAST)));
normalizeSVend(dim).. zPlus(dim) =e= (SVend(dim)-SVmin(dim))*dz(dim)-1;
ChebyBasisFunDeg0(dim).. ChebyBasis(dim,’1′) =e= 1;
ChebyBasisFunDeg1(dim).. ChebyBasis(dim,’2′) =e= zPlus(dim);
recursiveCheby(dim,d)$(ord(d)>=2 and ord(d)<card(d)).. ChebyBasis(dim,d+1) =e= 2*zPlus(dim)*ChebyBasis(dim,d) – ChebyBasis(dim,d-1); TerminalUtilEQ.. termUtil =e= sum((d,d2,d3,d4,d5,d6)$(ord(d)+ord(d2)+ord(d3)+ord(d4)+ord(d5)+ord(d6)<=card(d)+card(dim)-1), coefs(d,d2,d3,d4,d5,d6)*ChebyBasis(‘1’,d)*ChebyBasis(‘2’,d2)*ChebyBasis(‘3’,d3)*ChebyBasis(‘4’,d4)*ChebyBasis(‘5’,d5)*ChebyBasis(‘6’,d6)); ******************************* ** Upper and Lower Bounds K.lo(T) = 100; MAT.lo(T) = 10; MU.lo(t) = 100; ML.lo(t) = 1000; TOCEAN.up(T) = 20; TOCEAN.lo(T) = -1; TATM.lo(t) = 0; TATM.up(t) = 20; C.lo(T) = 20; miu.up(t) = LIMMIU; SVend.lo(dim) = SVmin(dim); SVend.up(dim) = SVmax(dim); ********************************** ** Solution options option iterlim = 10000; option reslim = 10000; option solprint = on; option limrow = 0; option limcol = 0; option nlp = conopt; model CO2 / all /; *** get the solution of CJL2 as initial guess *** set Thalf Time periods of CJL2 / 1*300 /; parameters C_CJL2(Thalf) consumption in CJL2 miu_CJL2(Thalf) emission control rates in CJL2 K_CJL2(Thalf) capital in CJL2; execute_load ‘sol_CJL2’, C_CJL2, miu_CJL2, K_CJL2; loop(Thalf, C.l(T)$(ord(T)=2*ord(Thalf)-1) = C_CJL2(Thalf); C.l(T)$(ord(T)=2*ord(Thalf)) = (C_CJL2(Thalf)+C_CJL2(Thalf+1))/2; miu.l(T)$(ord(T)=2*ord(Thalf)-1) = miu_CJL2(Thalf); miu.l(T)$(ord(T)=2*ord(Thalf)) = (miu_CJL2(Thalf)+miu_CJL2(Thalf+1))/2; K.l(T)$(ord(T)=2*ord(Thalf)-1) = K_CJL2(Thalf); K.l(T)$(ord(T)=2*ord(Thalf)) = (K_CJL2(Thalf)+K_CJL2(Thalf+1))/2; ); C.l(TLAST) = sum(Thalf$(ord(Thalf)=card(Thalf)),C_CJL2(Thalf)); miu.l(TLAST) = sum(Thalf$(ord(Thalf)=card(Thalf)), miu_CJL2(Thalf)); K.l(TLAST) = sum(Thalf$(ord(Thalf)=card(Thalf)),K_CJL2(Thalf)); E.l(T) = deltat*SIGMA(T)*(1-MIU.l(T))* (AL(T)*L(T)**(1-GAMA)*K.l(T)**GAMA) + ETREE(T); MAT.l(‘1’) = MAT2000; MU.l(‘1’) = MU2000; ML.l(‘1’) = ML2000; TATM.l(‘1’) = TATM0; TOCEAN.l(‘1’) = TOCEAN0; loop(T$(ord(T)<card(T)), MAT.l(T+1) = MAT.l(T)*bb11+MU.l(T)*bb21 + E.l(T); ML.l(T+1) = ML.l(T)*bb33+bb23*MU.l(T); MU.l(T+1) = MAT.l(T)*bb12+MU.l(T)*bb22+ML.l(T)*bb32; TATM.l(T+1) = TATM.l(t)+CC1*(FCO22X*(log(MAT.l(T)/596.4)/log(2))+FORCOTH(T) – LAM*TATM.l(t)-CC3*(TATM.l(t)-TOCEAN.l(t))); TOCEAN.l(T+1) = TOCEAN.l(T)+CC4*(TATM.l(T)-TOCEAN.l(T)); ); SVend.l(‘1’) = sum(TLAST, (1-DK)**deltat*K.l(TLAST)+deltat*( ( 1 – cost1(TLAST)*(MIU.l(TLAST)**EXPcost2) )* (AL(TLAST)*L(TLAST)**(1-GAMA)*K.l(TLAST)**GAMA) / (1+A1*TATM.l(TLAST)+ A2*TATM.l(TLAST)**A3) – C.l(TLAST) )); SVend.l(‘2’) = sum(TLAST, MAT.l(TLAST)*bb11+MU.l(TLAST)*bb21 + E.l(TLAST)); SVend.l(‘3’) = sum(TLAST, MAT.l(TLAST)*bb12+MU.l(TLAST)*bb22+ML.l(TLAST)*bb32); SVend.l(‘4’) = sum(TLAST, ML.l(TLAST)*bb33+bb23*MU.l(TLAST)); SVend.l(‘5’) = sum(TLAST, TATM.l(TLAST)+CC1*( FCO22X*(log(MAT.l(TLAST)/596.4)/log(2))+FORCOTH(TLAST) – LAM*TATM.l(TLAST)-CC3*(TATM.l(TLAST)-TOCEAN.l(TLAST)))); SVend.l(‘6′) = sum(TLAST, TOCEAN.l(TLAST)+CC4*(TATM.l(TLAST)-TOCEAN.l(TLAST))); zPlus.l(dim) = (SVend.l(dim)-SVmin(dim))*dz(dim)-1; ChebyBasis.l(dim,’1′) = 1; ChebyBasis.l(dim,’2’) = zPlus.l(dim); ChebyBasis.l(dim,d+1)$(ord(d)>=2 and ord(d)<card(d)) = 2*zPlus.l(dim)*ChebyBasis.l(dim,d) – ChebyBasis.l(dim,d-1); termUtil.l = sum((d,d2,d3,d4,d5,d6)$(ord(d)+ord(d2)+ord(d3)+ord(d4)+ord(d5)+ord(d6)<=card(d)+card(dim)-1), coefs(d,d2,d3,d4,d5,d6)*ChebyBasis.l(‘1’,d)*ChebyBasis.l(‘2’,d2)*ChebyBasis.l(‘3’,d3)*ChebyBasis.l(‘4’,d4)*ChebyBasis.l(‘5’,d5)*ChebyBasis.l(‘6’,d6)); UTILITY.l = sum( T, (((C.l(T)/L(T))**(1-B_ELASMU)-1)/(1-B_ELASMU)) * (RR(T)*L(T)*deltat/scale1) ) + termUtil.l * (beta*sum(TLAST,RR(TLAST))*POPASYM/scale1) ; ** solve the problem solve CO2 maximizing UTILITY using nlp ; ***************************** parameters C_CJL1(T) consumption in CJL1 miu_CJL1(T) emission control rates in CJL1 K_CJL1(T) capital in CJL1 ; C_CJL1(T) = C.l(T); miu_CJL1(T) = miu.l(T); K_CJL1(T) = K.l(T); execute_unload ‘sol_CJL1’, C_CJL1, miu_CJL1, K_CJL1; Parameters opt_tax(t) opt_mcemis(t); opt_tax(t)=-1*ee.m(t)*1000/(kk.m(t)+.00000000001); opt_mcemis(t)= expcost2*cost1(t)*miu.l(t)**(expcost2-1)/sigma(t)*1000; parameter CCA total carbon emission; CCA = sum(T, E.l(T)); display SVend.l, CCA, opt_tax, opt_mcemis; file CJL1_sol; put CJL1_sol; CJL1_sol.nw = 12; CJL1_sol.nr = 2; CJL1_sol.nz = 1e-15; loop(t, put t.tl:4:0; put K.l(t):14:6; put MAT.l(t):14:6; put MU.l(t):14:6; put ML.l(t):14:6; put TATM.l(t):14:6; put TOCEAN.l(t):14:6; put C.l(t):14:6; put miu.l(t):14:6; put E.l(t):14:6; put opt_tax(t):14:6; put opt_mcemis(t):14:6; put /; ); file params_CJL1; put params_CJL1; params_CJL1.nw = 12; params_CJL1.nr = 2; params_CJL1.nz = 1e-15; loop(t, put K.l(t):14:6 /; ); loop(t, put MAT.l(t):14:6 /; ); loop(t, put MU.l(t):14:6 /; ); loop(t, put ML.l(t):14:6 /; ); loop(t, put TATM.l(t):14:6 /; ); loop(t, put TOCEAN.l(t):14:6 /; ); loop(t, put C.l(t):14:6 /; ); loop(t, put miu.l(t):14:6 /; ); loop(t, put E.l(t):14:6 /; );

• CJL1.lst (530k)
• sol_CJL2.gdx (9k)