Help from Patrons: Any GAMS users, or mathematicians familiar with GAMS optimization routines?

I’m slogging through Nordhaus’s DICE model in preparation for a review paper that I’ve been invited to submit to a leading science journal, as a follow-up to the Globalizations paper “The appallingly bad neoclassical economics of climate change” (click here for a Patreon post with the paper in PDF if you haven’t yet seen it: the publishers have moved my paper behind a paywall). This is why I haven’t been posting all that much in the last couple of weeks, because this is hard going, not because DICE is complicated, but because…

Well blow me down with a feather, if there isn’t an error in DICE!

DICE is a typical Neoclassical Ramsey growth model [1], where a utility-maximizing representative agent with perfect foresight trades off consumption (via work, which is a source of disutility) against leisure to determine an optimal growth path, and in which what is not consumed (in the process of maximizing utility) is saved. This savings determines investment, investment determines the capital stock, which determines output (along with assumed growth in “total factor productivity” and population, both of which are assumed to occur and are exogenous to the model). The dilemma (for Neoclassical economists) is that the equilibrium of the Ramsey model is unstable: it’s a “saddle point equilibrium” where one eigenvector has a negative eigenvalue, and another has a positive eigenvalue. So Ramsey’s solution was to impose a “benevolent dictator” who knew the future “bliss point” where the future equilibrium lay, could work out the stable eigenvector of this model, intrapolate (I can’t call it extrapolation, because that goes forwards in time, not backwards as here) it back to today’s conditions, and if the consumption and savings levels aren’t on the stable eigenvector, move them instantly to the consumption and savings combo that is on the stable eigenvector. Then the economy can move forwards from there to the future (unstable) bliss point.

Ah, for a moment there, it felt like I was back in Amsterdam (in a coffee shop, rather than an anti-lockdown riot). Yes, it’s insane. But that model has become the basis of Neoclassical modelling, with the idea of a “benevolent social planner” replaced by a rational representative agent, where rational means “someone who can predict the future”.

There’s just one problem with Nordhaus’s implementation of this Ramsay growth model (this Wikipedia explanation of the model is pretty good): the model depends on savings, but he forgot to include an equation for the savings rate!

And yet the model still runs in GAMS, . I’d like to work out why, and I would appreciate some help from Patrons, if there are any that are either (a) experienced users of GAMS and/or (b) mathematicians familiar with the solution to optimization routines that GAMS uses. My intuitions are (a) that the model works because the optimization routine works backwards from the endpoint of the simulation to its initial conditions, and (b) that endpoint is given by the parameter optlrsav, which is an imposed rate of savings needed to force the model to converge to its unstable equilibrium.

In the verbatim listing of the latest published version of the model (which he still calls a “beta version” even though there’s been no published update to it for five years) shown below, I’ve put in bold the equations relevant to determining output. They’re extracted here for easy reference:

Nordhaus’s economic growth equations

Output Equations

YGROSS(t) Gross world product GROSS of abatement and damages (trillions 2005 USD per year)

ygrosseq(t).. YGROSS(t) =E= (al(t)*(L(t)/1000)**(1-GAMA))*(K(t)**GAMA);

Total Factor Productivity equations (al(t))

al(t) Level of total factor productivity

a0 Initial level of total factor productivity /5.115 /

al(“1”) = a0; loop(t, al(t+1)=al(t)/((1-ga(t))););

ga(t) Growth rate of productivity from

ga(t)=ga0*exp(-dela*5*((t.val-1)));

ga0 Initial growth rate for TFP per 5 years /0.076 /

dela Decline rate of TFP per 5 years /0.005 /

Labor equations (L(t))

l(t) Level of population and labor

l(“1”) = pop0;

loop(t, l(t+1)=l(t););

loop(t, l(t+1)=l(t)*(popasym/L(t))**popadj ;);

pop0 Initial world population 2015 (millions) /7403 /

popadj Growth rate to calibrate to 2050 pop projection /0.134 /

popasym Asymptotic population (millions) /11500 /

gl(t) Growth rate of labor

gfacpop(t) Growth factor population

 

Capital Equations

K(t) Capital stock (trillions 2005 US dollars)

gama Capital elasticity in production function /.300 /

dk Depreciation rate on capital (per year) /.100 /

I(t) Investment (trillions 2005 USD per year)

S(t) Gross savings rate as fraction of gross world product

kk(t+1).. K(t+1) =L= (1-dk)**tstep * K(t) + tstep * I(t);

seq(t).. I(t) =E= S(t) * Y(t);

cc(t).. C(t) =E= Y(t) – I(t);

 

You will (if you’re a modeler, or an astute reader) notice that there is no equation for S(t). So how the hell does the program still run? I think it’s because of what Neoclassicals call the “transversality condition” that is used to force equilibrium on the model…

If all else fails…

set lag10(t) ;

lag10(t) = yes$(t.val gt card(t)-10);

S.FX(lag10(t)) = optlrsav;

Diagnosing why an incomplete model still runs

My guess is that GAMS uses the definition for S.FX (where FX is a GAMS operator that provides a fixed value for the variable preceding it, if no other value exists), and backward iterates from this to the start of the simulation, thus enabling the model to work even though it is incomplete.

But I need to be sure of this. I’m still waiting on word from some of my academic colleagues who might have the relevant background, but I’d be mighty pleased if there was someone here who could help suss this out. The optimization routine that Nordhaus uses is “nlp”:

solve co2 maximizing utility using nlp;

 

So, can anyone assist? The basic DICE models are here, and need GAMS to run.S

Nordhaus’s DICE model (verbatim)

$ontext

This is the beta version of DICE-2016R. The major changes are outlined in Nordhaus,

“Revisiting the social cost of carbon: Estimates from the DICE-2016R model,”

September 30, 2016,” available from the author.

 

Version is DICE-2016R-091916ap.gms

$offtext

 

$title DICE-2016R September 2016 (DICE-2016R-091216a.gms)

 

set t Time periods (5 years per period) /1*100/

 

PARAMETERS

** Availability of fossil fuels

fosslim Maximum cumulative extraction fossil fuels (GtC) /6000/

**Time Step

tstep Years per Period /5/

** If optimal control

ifopt Indicator where optimized is 1 and base is 0 /0/

** Preferences

elasmu Elasticity of marginal utility of consumption /1.45 /

prstp Initial rate of social time preference per year /.015 /

** Population and technology

gama Capital elasticity in production function /.300 /

pop0 Initial world population 2015 (millions) /7403 /

popadj Growth rate to calibrate to 2050 pop projection /0.134 /

popasym Asymptotic population (millions) /11500 /

dk Depreciation rate on capital (per year) /.100 /

q0 Initial world gross output 2015 (trill 2010 USD) /105.5 /

k0 Initial capital value 2015 (trill 2010 USD) /223 /

a0 Initial level of total factor productivity /5.115 /

ga0 Initial growth rate for TFP per 5 years /0.076 /

dela Decline rate of TFP per 5 years /0.005 /

** Emissions parameters

gsigma1 Initial growth of sigma (per year) /-0.0152 /

dsig Decline rate of decarbonization (per period) /-0.001 /

eland0 Carbon emissions from land 2015 (GtCO2 per year) / 2.6 /

deland Decline rate of land emissions (per period) / .115 /

e0 Industrial emissions 2015 (GtCO2 per year) /35.85 /

miu0 Initial emissions control rate for base case 2015 /.03 /

** Carbon cycle

* Initial Conditions

mat0 Initial Concentration in atmosphere 2015 (GtC) /851 /

mu0 Initial Concentration in upper strata 2015 (GtC) /460 /

ml0 Initial Concentration in lower strata 2015 (GtC) /1740 /

mateq Equilibrium concentration atmosphere (GtC) /588 /

mueq Equilibrium concentration in upper strata (GtC) /360 /

mleq Equilibrium concentration in lower strata (GtC) /1720 /

* Flow paramaters

b12 Carbon cycle transition matrix /.12 /

b23 Carbon cycle transition matrix /0.007 /

* These are for declaration and are defined later

b11 Carbon cycle transition matrix

b21 Carbon cycle transition matrix

b22 Carbon cycle transition matrix

b32 Carbon cycle transition matrix

b33 Carbon cycle transition matrix

sig0 Carbon intensity 2010 (kgCO2 per output 2005 USD 2010)

** Climate model parameters

t2xco2 Equilibrium temp impact (oC per doubling CO2) / 3.1 /

fex0 2015 forcings of non-CO2 GHG (Wm-2) / 0.5 /

fex1 2100 forcings of non-CO2 GHG (Wm-2) / 1.0 /

tocean0 Initial lower stratum temp change (C from 1900) /.0068 /

tatm0 Initial atmospheric temp change (C from 1900) /0.85 /

c1 Climate equation coefficient for upper level /0.1005 /

c3 Transfer coefficient upper to lower stratum /0.088 /

c4 Transfer coefficient for lower level /0.025 /

fco22x Forcings of equilibrium CO2 doubling (Wm-2) /3.6813 /

** Climate damage parameters

a10 Initial damage intercept /0 /

a20 Initial damage quadratic term

a1 Damage intercept /0 /

a2 Damage quadratic term /0.00236 /

a3 Damage exponent /2.00 /

** Abatement cost

expcost2 Exponent of control cost function / 2.6 /

pback Cost of backstop 2010$ per tCO2 2015 / 550 /

gback Initial cost decline backstop cost per period / .025 /

limmiu Upper limit on control rate after 2150 / 1.2 /

tnopol Period before which no emissions controls base / 45 /

cprice0 Initial base carbon price (2010$ per tCO2) / 2 /

gcprice Growth rate of base carbon price per year /.02 /

 

** Scaling and inessential parameters

* Note that these are unnecessary for the calculations

* They ensure that MU of first period’s consumption =1 and PV cons = PV utilty

scale1 Multiplicative scaling coefficient /0.0302455265681763 /

scale2 Additive scaling coefficient /-10993.704/ ;

 

* Program control variables

sets tfirst(t), tlast(t), tearly(t), tlate(t);

 

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 rate

ga(t) Growth rate of productivity from

forcoth(t) Exogenous forcing for other greenhouse gases

gl(t) Growth rate of labor

gcost1 Growth of cost factor

gsig(t) Change in sigma (cumulative improvement of energy efficiency)

etree(t) Emissions from deforestation

cumetree(t) Cumulative from land

cost1(t) Adjusted cost for backstop

lam Climate model parameter

gfacpop(t) Growth factor population

pbacktime(t) Backstop price

optlrsav Optimal long-run savings rate used for transversality

scc(t) Social cost of carbon

cpricebase(t) Carbon price in base case

photel(t) Carbon Price under no damages (Hotelling rent condition)

ppm(t) Atmospheric concentrations parts per million

atfrac(t) Atmospheric share since 1850

atfrac2010(t) Atmospheric share since 2010 ;

* Program control definitions

tfirst(t) = yes$(t.val eq 1);

tlast(t) = yes$(t.val eq card(t));

* Parameters for long-run consistency of carbon cycle

b11 = 1 – b12;

b21 = b12*MATEQ/MUEQ;

b22 = 1 – b21 – b23;

b32 = b23*mueq/mleq;

b33 = 1 – b32 ;

* Further definitions of parameters

a20 = a2;

sig0 = e0/(q0*(1-miu0));

lam = fco22x/ t2xco2;

l(“1”) = pop0;

loop(t, l(t+1)=l(t););

loop(t, l(t+1)=l(t)*(popasym/L(t))**popadj ;);

 

ga(t)=ga0*exp(-dela*5*((t.val-1)));

al(“1”) = a0; loop(t, al(t+1)=al(t)/((1-ga(t))););

gsig(“1”)=gsigma1; loop(t,gsig(t+1)=gsig(t)*((1+dsig)**tstep) ;);

sigma(“1”)=sig0; loop(t,sigma(t+1)=(sigma(t)*exp(gsig(t)*tstep)););

 

pbacktime(t)=pback*(1-gback)**(t.val-1);

cost1(t) = pbacktime(t)*sigma(t)/expcost2/1000;

 

etree(t) = eland0*(1-deland)**(t.val-1);

cumetree(“1”)= 100; loop(t,cumetree(t+1)=cumetree(t)+etree(t)*(5/3.666););

 

rr(t) = 1/((1+prstp)**(tstep*(t.val-1)));

forcoth(t) = fex0+ (1/17)*(fex1-fex0)*(t.val-1)$(t.val lt 18)+ (fex1-fex0)$(t.val ge 18);

optlrsav = (dk + .004)/(dk + .004*elasmu + prstp)*gama;

 

*Base Case Carbon Price

cpricebase(t)= cprice0*(1+gcprice)**(5*(t.val-1));

 

VARIABLES

MIU(t) Emission control rate GHGs

FORC(t) Increase in radiative forcing (watts per m2 from 1900)

TATM(t) Increase temperature of atmosphere (degrees C from 1900)

TOCEAN(t) Increase temperatureof lower oceans (degrees C from 1900)

MAT(t) Carbon concentration increase in atmosphere (GtC from 1750)

MU(t) Carbon concentration increase in shallow oceans (GtC from 1750)

ML(t) Carbon concentration increase in lower oceans (GtC from 1750)

E(t) Total CO2 emissions (GtCO2 per year)

EIND(t) Industrial emissions (GtCO2 per year)

C(t) Consumption (trillions 2005 US dollars per year)

K(t) Capital stock (trillions 2005 US dollars)

CPC(t) Per capita consumption (thousands 2005 USD per year)

I(t) Investment (trillions 2005 USD per year)

S(t) Gross savings rate as fraction of gross world product

RI(t) Real interest rate (per annum)

Y(t) Gross world product net of abatement and damages (trillions 2005 USD per year)

YGROSS(t) Gross world product GROSS of abatement and damages (trillions 2005 USD per year)

YNET(t) Output net of damages equation (trillions 2005 USD per year)

DAMAGES(t) Damages (trillions 2005 USD per year)

DAMFRAC(t) Damages as fraction of gross output

ABATECOST(t) Cost of emissions reductions (trillions 2005 USD per year)

MCABATE(t) Marginal cost of abatement (2005$ per ton CO2)

CCA(t) Cumulative industrial carbon emissions (GTC)

CCATOT(t) Total carbon emissions (GtC)

PERIODU(t) One period utility function

CPRICE(t) Carbon price (2005$ per ton of CO2)

CEMUTOTPER(t) Period utility

UTILITY Welfare function;

 

NONNEGATIVE VARIABLES MIU, TATM, MAT, MU, ML, Y, YGROSS, C, K, I;

 

EQUATIONS

*Emissions and Damages

EEQ(t) Emissions equation

EINDEQ(t) Industrial emissions

CCACCA(t) Cumulative industrial carbon emissions

CCATOTEQ(t) Cumulative total carbon emissions

FORCE(t) Radiative forcing equation

DAMFRACEQ(t) Equation for damage fraction

DAMEQ(t) Damage equation

ABATEEQ(t) Cost of emissions reductions equation

MCABATEEQ(t) Equation for MC abatement

CARBPRICEEQ(t) Carbon price equation from abatement

 

*Climate and carbon cycle

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

 

*Economic variables

YGROSSEQ(t) Output gross equation

YNETEQ(t) Output net of damages equation

YY(t) Output net equation

CC(t) Consumption equation

CPCE(t) Per capita consumption definition

SEQ(t) Savings rate equation

KK(t) Capital balance equation

RIEQ(t) Interest rate equation

 

* Utility

CEMUTOTPEREQ(t) Period utility

PERIODUEQ(t) Instantaneous utility function equation

UTIL Objective function ;

 

** Equations of the model

*Emissions and Damages

eeq(t).. E(t) =E= EIND(t) + etree(t);

eindeq(t).. EIND(t) =E= sigma(t) * YGROSS(t) * (1-(MIU(t)));

ccacca(t+1).. CCA(t+1) =E= CCA(t)+ EIND(t)*5/3.666;

ccatoteq(t).. CCATOT(t) =E= CCA(t)+cumetree(t);

force(t).. FORC(t) =E= fco22x * ((log((MAT(t)/588.000))/log(2))) + forcoth(t);

damfraceq(t) .. DAMFRAC(t) =E= (a1*TATM(t))+(a2*TATM(t)**a3) ;

dameq(t).. DAMAGES(t) =E= YGROSS(t) * DAMFRAC(t);

abateeq(t).. ABATECOST(t) =E= YGROSS(t) * cost1(t) * (MIU(t)**expcost2);

mcabateeq(t).. MCABATE(t) =E= pbacktime(t) * MIU(t)**(expcost2-1);

carbpriceeq(t).. CPRICE(t) =E= pbacktime(t) * (MIU(t))**(expcost2-1);

 

*Climate and carbon cycle

mmat(t+1).. MAT(t+1) =E= MAT(t)*b11 + MU(t)*b21 + (E(t)*(5/3.666));

mml(t+1).. ML(t+1) =E= ML(t)*b33 + MU(t)*b23;

mmu(t+1).. MU(t+1) =E= MAT(t)*b12 + MU(t)*b22 + ML(t)*b32;

tatmeq(t+1).. TATM(t+1) =E= TATM(t) + c1 * ((FORC(t+1)-(fco22x/t2xco2)*TATM(t))-(c3*(TATM(t)-TOCEAN(t))));

toceaneq(t+1).. TOCEAN(t+1) =E= TOCEAN(t) + c4*(TATM(t)-TOCEAN(t));

 

*Economic variables

ygrosseq(t).. YGROSS(t) =E= (al(t)*(L(t)/1000)**(1-GAMA))*(K(t)**GAMA);

yneteq(t).. YNET(t) =E= YGROSS(t)*(1-damfrac(t));

yy(t).. Y(t) =E= YNET(t) – ABATECOST(t);

cc(t).. C(t) =E= Y(t) – I(t);

cpce(t).. CPC(t) =E= 1000 * C(t) / L(t);

seq(t).. I(t) =E= S(t) * Y(t);

kk(t+1).. K(t+1) =L= (1-dk)**tstep * K(t) + tstep * I(t);

rieq(t+1).. RI(t) =E= (1+prstp) * (CPC(t+1)/CPC(t))**(elasmu/tstep) – 1;

 

*Utility

cemutotpereq(t).. CEMUTOTPER(t) =E= PERIODU(t) * L(t) * rr(t);

periodueq(t).. PERIODU(t) =E= ((C(T)*1000/L(T))**(1-elasmu)-1)/(1-elasmu)-1;

util.. UTILITY =E= tstep * scale1 * sum(t, CEMUTOTPER(t)) + scale2 ;

 

*Resource limit

CCA.up(t) = fosslim;

 

* Control rate limits

MIU.up(t) = limmiu;

MIU.up(t)$(t.val<30) = 1;

 

** Upper and lower bounds for stability

K.LO(t) = 1;

MAT.LO(t) = 10;

MU.LO(t) = 100;

ML.LO(t) = 1000;

C.LO(t) = 2;

TOCEAN.UP(t) = 20;

TOCEAN.LO(t) = -1;

TATM.UP(t) = 20;

CPC.LO(t) = .01;

TATM.UP(t) = 12;

 

* Control variables

set lag10(t) ;

lag10(t) = yes$(t.val gt card(t)-10);

S.FX(lag10(t)) = optlrsav;

 

* Initial conditions

CCA.FX(tfirst) = 400;

K.FX(tfirst) = k0;

MAT.FX(tfirst) = mat0;

MU.FX(tfirst) = mu0;

ML.FX(tfirst) = ml0;

TATM.FX(tfirst) = tatm0;

TOCEAN.FX(tfirst) = tocean0;

 

** Solution options

option iterlim = 99900;

option reslim = 99999;

option solprint = on;

option limrow = 0;

option limcol = 0;

model CO2 /all/;

 

* For base run, this subroutine calculates Hotelling rents

* Carbon price is maximum of Hotelling rent or baseline price

* The cprice equation is different from 2013R. Not sure what went wrong.

If (ifopt eq 0,

a2 = 0;

solve CO2 maximizing UTILITY using nlp;

photel(t)=cprice.l(t);

a2 = a20;

cprice.up(t)$(t.val<tnopol+1) = max(photel(t),cpricebase(t));

);

 

miu.fx(‘1’)$(ifopt=1) = miu0;

solve co2 maximizing utility using nlp;

solve co2 maximizing utility using nlp;

solve co2 maximizing utility using nlp;

 

** POST-SOLVE

* Calculate social cost of carbon and other variables

scc(t) = -1000*eeq.m(t)/(.00001+cc.m(t));

atfrac(t) = ((mat.l(t)-588)/(ccatot.l(t)+.000001 ));

atfrac2010(t) = ((mat.l(t)-mat0)/(.00001+ccatot.l(t)-ccatot.l(‘1’) ));

ppm(t) = mat.l(t)/2.13;

 

* Produces a file “Dice2016R-091916ap.csv” in the base directory

* For ALL relevant model outputs, see ‘PutOutputAllT.gms’ in the Include folder.

* The statement at the end of the *.lst file “Output…” will tell you where to find the file.

 

file results /Dice2016R-091916ap.csv/; results.nd = 10 ; results.nw = 0 ; results.pw=20000; results.pc=5;

put results;

put /”Results of DICE-2016R model run using model Dice2016R-091916ap.csv”;

put /”This is optimal if ifopt = 1 and baseline if ifopt = 0″;

put /”ifopt =” ifopt;

put // “Period”;

Loop (T, put T.val);

put / “Year” ;

Loop (T, put (2010+(TSTEP*T.val) ));

put / “Industrial Emissions GTCO2 per year” ;

Loop (T, put EIND.l(T));

put / “Atmospheric concentration C (ppm)” ;

Loop (T, put (MAT.l(T)/2.13));

put / “Atmospheric Temperature ” ;

Loop (T, put TATM.l(T));

put / “Output Net Net) ” ;

Loop (T, put Y.l(T));

put / “Climate Damages fraction output” ;

Loop (T, put DAMFRAC.l(T));

put / “Consumption Per Capita ” ;

Loop (T, put CPC.l(T));

put / “Carbon Price (per t CO2)” ;

Loop (T, put cprice.l(T));

put / “Emissions Control Rate” ;

Loop (T, put MIU.l(T));

put / “Social cost of carbon” ;

Loop (T, put scc(T));

put / “Interest Rate ” ;

Loop (T, put RI.l(T));

put / “Population” ;

Loop (T, put L(T));

put / “TFP” ;

Loop (T, put AL(T));

put / “Output gross,gross” ;

Loop (T, put YGROSS.L(t));

put / “Change tfp” ;

Loop (T, put ga(t));

put / “Capital” ;

Loop (T, put k.l(t));

put / “s” ;

Loop (T, put s.l(t));

put / “I” ;

Loop (T, put I.l(t));

put / “Y gross net” ;

Loop (T, put ynet.l(t));

put / “damages” ;

Loop (T, put damages.l(t));

put / “damfrac” ;

Loop (T, put damfrac.l(t));

put / “abatement” ;

Loop (T, put abatecost.l(t));

put / “sigma” ;

Loop (T, put sigma(t));

put / “Forcings” ;

Loop (T, put forc.l(t));

put / “Other Forcings” ;

Loop (T, put forcoth(t));

put / “Period utilty” ;

Loop (T, put periodu.l(t));

put / “Consumption” ;

Loop (T, put C.l(t));

put / “Objective” ;

put utility.l;

put / “Land emissions” ;

Loop (T, put etree(t));

put / “Cumulative ind emissions” ;

Loop (T, put cca.l(t));

put / “Cumulative total emissions” ;

Loop (T, put ccatot.l(t));

put / “Atmospheric concentrations Gt” ;

Loop (T, put mat.l(t));

put / “Atmospheric concentrations ppm” ;

Loop (T, put ppm(t));

put / “Total Emissions GTCO2 per year” ;

Loop (T, put E.l(T));

put / “Atmospheric concentrations upper” ;

Loop (T, put mu.l(t));

put / “Atmospheric concentrations lower” ;

Loop (T, put ml.l(t));

put / “Atmospheric fraction since 1850” ;

Loop (T, put atfrac(t));

put / “Atmospheric fraction since 2010” ;

Loop (T, put atfrac2010(t));

putclose;

 

 

[1] Ramsey, F.P. 1928 A Mathematical Theory of Saving. The Economic Journal
38, 543-559. (doi:10.2307/2224098).