apis/examples_python/nordhausdice_8py_source.html

if __name__ == "__main__":


    # [1]

    from gams.magic import GamsInteractive

    gams = GamsInteractive()


    # [2]

    gams.gams('''

# 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


# DICE-2016R September 2016 (adapted from 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 capital 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

model CO2 / all /;''')


    # [3]

    gams.gams('''

MIU.FX('1') = miu0;

solve CO2 maximizing UTILITY using nlp;''')


    # [4]

    import pandas as pd

    m = gams.exchange_container

    results = {}

    for d in ['TATM','RI','EIND','CPRICE']:

        results[d] = pd.DataFrame(columns=['opt','bas'])

        results[d]['opt'] = m[d].records['level'].copy()


    # [5]

    gams.gams('''

# 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.

option bRatio = 1; MIU.LO('1') = 0; MIU.UP('1') = inf; # Reset fixing from opt run and discard basis

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));

option bRatio=0.25;

solve CO2 maximizing UTILITY using nlp;''')


    # [6]

    tstep = int(m['tstep'].toValue())

    l = [ 2010 + tstep*n for n in range(len(m['t'].records))]

    for d in ['TATM','RI','EIND','CPRICE']:

        results[d]['bas'] = m[d].records['level'].copy()

        # remap index from 1*100 to 2010,2015,2020,...,2505

        results[d].set_index(pd.Index(l[:len(m[d].records['level'])]), inplace=True)


    # [7] - [10] - skipped


    # [11]

    gams.gams('''

## 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;


Parameter rep(*,t);

rep("Industrial Emissions GTCO2 per year", t) = EIND.L(t);

rep("Atmospheric concentration C (ppm)", t) = MAT.L(t)/2.13;

rep("Atmospheric Temperature", t) = TATM.L(t);

rep("Output Net Net)", t) = Y.L(t);

rep("Climate Damages fraction output", t) = DAMFRAC.L(t);

rep("Consumption Per Capita", t) = CPC.L(t);

rep("Carbon Price (per t CO2)", t) = CPRICE.L(t);

rep("Emissions Control Rate", t) = MIU.L(t);

rep("Social cost of carbon", t) = scc(t);

rep("Interest Rate", t) = RI.L(t);

rep("Population", t) = l(t);

rep("TFP", t) = al(t);

rep("Output gross,gross", t) = YGROSS.L(t);

rep("Change tfp", t) = ga(t);

rep("Capital", t) = K.L(t);

rep("s", t) = S.L(t);

rep("I", t) = I.L(t);

rep("Y gross net", t) = YNET.L(t);

rep("damages", t) = DAMAGES.L(t);

rep("damfrac", t) = DAMFRAC.L(t);

rep("abatement", t) = ABATECOST.L(t);

rep("sigma", t) = sigma(t);

rep("Forcings", t) = FORC.L(t);

rep("Other Forcings", t) = forcoth(t);

rep("Period utilty", t) = PERIODU.L(t);

rep("Consumption", t) = C.L(t);

rep("Land emissions", t) = etree(t);

rep("Cumulative ind emissions", t) = CCA.L(t);

rep("Cumulative total emissions", t) = CCATOT.L(t);

rep("Atmospheric concentrations Gt", t) = MAT.L(t);

rep("Atmospheric concentrations ppm", t) = ppm(t);

rep("Total Emissions GTCO2 per year", t) = E.L(t);

rep("Atmospheric concentrations upper", t) = MU.L(t);

rep("Atmospheric concentrations lower", t) = ML.L(t);

rep("Atmospheric fraction since 1850", t) = atfrac(t);

rep("Atmospheric fraction since 2010", t) = atfrac2010(t);''')


    # [12]

    import copy


    rep = copy.deepcopy(m['rep'])

    rep.records['t'] = rep.records['t'].map({s:2010+tstep*int(s) for s in m['t'].records['uni']})

    print(rep.pivot(index='t', columns='uni'))


    # [13]

    gams.gams_cleanup(closedown=True)