This is the source code to generate terminal value function for CJL.
=======================================================
$ontext
Get the terminal value function of Cimate project
$offtext
SET T Time periods /1*800/ ;
scalar deltat interval length of each period (number of years) / 1 /;
set n index /1*5/;
alias(n,n2,n3,n4,n5,n6,d,d2,d3,d4,d5,d6);
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
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 /
** 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 /
** Carbon cycle
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 /
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 /
;
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
GGBACK Initial cost decline backstop pc
;
***********************************
set dim indices of dimensions /1*6/;
parameters
SVmin(dim) minimal state variables
SVmax(dim) maximal state variables
LBrate lower bound rate
UBrate upper bound rate
;
LBrate = 0.8;
UBrate = 1.5;
SVmin(‘1’) = 80000*0.5;
SVmax(‘1’) = 80000*2;
SVmin(‘2’) = 900*LBrate;
SVmax(‘2’) = 900*UBrate;
SVmin(‘3’) = 1600*LBrate;
SVmax(‘3’) = 1600*UBrate;
SVmin(‘4’) = 20000*LBrate;
SVmax(‘4’) = 20000*UBrate;
SVmin(‘5’) = 1;
SVmax(‘5’) = 3;
SVmin(‘6’) = 1.5;
SVmax(‘6’) = 3.5;
parameter z(n) standard Chebyshev nodes;
z(n) = – cos((2*ord(n)-1)/(2*card(n))*pi);
parameter ext(dim) expanded Chebyshev parameter;
*ext(dim) = -(z(‘1’)+1)*(SVmax(dim)-SVmin(dim))/2/z(‘1’);
ext(dim) = 0;
parameter extSVmin(dim) expanded minimal state value;
extSVmin(dim) = SVmin(dim)-ext(dim);
parameter extSVmax(dim) expanded maximal state value;
extSVmax(dim) = SVmax(dim)+ext(dim);
parameter dz(dim) derivative of z over x;
dz(dim) = 2/(extSVmax(dim)-extSVmin(dim));
parameter ChebyT(n,d) Chebyshev polynomial values at nodes;
ChebyT(n,’1′) = 1;
ChebyT(n,’2′) = z(n);
loop(d$(ord(d)>=2 and ord(d)<card(d)), ChebyT(n,d+1) = 2*z(n)* ChebyT(n,d) – ChebyT(n,d-1); ); parameter coefs(d,d2,d3,d4,d5,d6) complete Chebyshev coefficients; ************************************ ** compute terminal value function ************************************ parameters beta one-period discount factor LAM Climate model parameter sigmaEnd CO2-equivalent-emissions output ratio cost1End Adjusted cost for backstop abatecost1End Cost of emissions reductions InvestRatio investment ratio of capital ALend Level of total factor productivity FORCEs(T) Radiative forcing in watts per m2 CAs(T) Carbon concentration in atmosphere GtC CSs(T) Carbon concentration in shallow oceans Gtc CDs(T) Carbon concentration in lower oceans GtC TAs(T) Temperature of atmosphere in degrees C TLs(T) Temperature of lower oceans degrees C Ks(T) Capital stock trillions US dollars Qs(T) Gross world product net of abatement and damages Cs(T) Consumption trillions US dollars DamageFuns(T) damage function utils(T) utility AAL(T) Level of total factor productivity SSIGMA(T) CO2-equivalent-emissions output ratio GGA(T) Growth rate of productivity from 0 to T GGSIG(T) Cumulative improvement of energy efficiency ; parameter totalUtil(n,n2,n3,n4,n5,n6) total utility; 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; GGA0 = GA0*deltat; DDELA = DELA*deltat; GGSIGMA = GSIGMA*deltat; DDSIG = DSIG*deltat; * Important parameters for the model LAM = FCO22X/ T2XCO2; gga(T)=GGA0*EXP(-DDELA*(ORD(T)-1)); aal(“1”) = a0; LOOP(T, aal(T+1)=aal(T)/((1-gga(T)));); ggsig(T)=GGSIGMA*EXP(-DDSIG*(ORD(T)-1)); ssigma(“1”)=sig0;LOOP(T,ssigma(T+1)=(ssigma(T)/((1-ggsig(T+1))));); beta = 1 / ((1+B_PRSTP)**deltat); ******************************** sigmaEnd = sum(T$(ord(T)=card(T)), SSIGMA(T)); cost1End = (PBACK*sigmaEnd/EXPCOST2) * ((BACKRAT-1)/BACKRAT); abatecost1End = cost1End; InvestRatio = (1-(1-DK)**deltat)/deltat; ALend = sum(T$(ord(T)=card(T)), AAL(T)); loop(n, Ks(“1”) = extSVmin(‘1’) + (1+z(n))/dz(‘1’); loop(n2, CAs(“1”) = extSVmin(‘2’) + (1+z(n2))/dz(‘2’); loop(n3, CSs(“1”) = extSVmin(‘3’) + (1+z(n3))/dz(‘3’); loop(n4, CDs(“1”) = extSVmin(‘4’) + (1+z(n4))/dz(‘4’); loop(n5, TAs(“1”) = extSVmin(‘5’) + (1+z(n5))/dz(‘5’); loop(n6, TLs(“1”) = extSVmin(‘6’) + (1+z(n6))/dz(‘6’); loop(T$(ord(T)<card(T)), Ks(T+1) = Ks(T); FORCEs(T) = FCO22X*((log((CAs(T))/596.4)/log(2)))+0.3; CAs(T+1) = CAs(T)*bb11+CSs(T)*bb21; CSs(T+1) = CAs(T)*bb12+CSs(T)*bb22+CDs(T)*bb32; CDs(T+1) = CDs(T)*bb33+bb23*CSs(T); TAs(T+1) = TAs(T)+CC1*(FORCEs(T)-LAM*TAs(T)-CC3*(TAs(T)-TLs(T))); TLs(T+1) = TLs(T)+CC4*(TAs(T)-TLs(T)); ); DamageFuns(T) = (1 / (1+A1*TAs(T)+A2*TAs(T)**A3))$(TAs(T)>0) + 1$(TAs(T)<=0); Qs(T) = (1-abatecost1End)*ALend*POPASYM**(1-GAMA)*Ks(T)**GAMA * DamageFuns(T); Cs(T) = Qs(T)-InvestRatio*Ks(T); utils(T) = ((Cs(T)/POPASYM)**(1-B_ELASMU)-1)/(1-B_ELASMU); totalUtil(n,n2,n3,n4,n5,n6) = SUM(T,beta**(ord(T)-1)*utils(T)); ); ); ); ); ); ); coefs(d,d2,d3,d4,d5,d6) = 0; coefs(d,d2,d3,d4,d5,d6)$(ord(d)+ord(d2)+ord(d3)+ord(d4)+ord(d5)+ord(d6)<=card(d)+card(dim)-1) = 64*sum((n,n2,n3,n4,n5,n6), ChebyT(n,d)* ChebyT(n2,d2)* ChebyT(n3,d3)* ChebyT(n4,d4)* ChebyT(n5,d5)* ChebyT(n6,d6)* totalUtil(n,n2,n3,n4,n5,n6)) / (card(n)*card(n2)*card(n3)*card(n4)*card(n5)*card(n6)); coefs(‘1’,d2,d3,d4,d5,d6) = coefs(‘1′,d2,d3,d4,d5,d6)/2; coefs(d,’1′,d3,d4,d5,d6) = coefs(d,’1′,d3,d4,d5,d6)/2; coefs(d,d2,’1′,d4,d5,d6) = coefs(d,d2,’1′,d4,d5,d6)/2; coefs(d,d2,d3,’1′,d5,d6) = coefs(d,d2,d3,’1′,d5,d6)/2; coefs(d,d2,d3,d4,’1′,d6) = coefs(d,d2,d3,d4,’1′,d6)/2; coefs(d,d2,d3,d4,d5,’1′) = coefs(d,d2,d3,d4,d5,’1’)/2; ******************************** execute_unload ‘coefs_Climate_termValue’,SVmin, SVmax, dz, coefs; file coefs_Climate_termValue; put coefs_Climate_termValue; coefs_Climate_termValue.nw = 12; coefs_Climate_termValue.nr = 2; coefs_Climate_termValue.nz = 1e-15; loop(dim, put SVmin(dim):14:6 /; put SVmax(dim):14:6 /; ); loop((d,d2,d3,d4,d5,d6), put coefs(d,d2,d3,d4,d5,d6):14:6 /; );