The DICEModel Module

DICEModel

Implementation of the DICE/RICE models

Notes:

• Optimization based on DICE2023-b-4-3-10.gms and included files (Nonco2-b-4-3-1.gms and FAIR-beta-4-3-1.gms)
• Variable casing has been harmonized that all parameters and post-optimization computation have lower cases, and all optimization variables have upper case.
• DICE has been generalized to work with multiple regions, and it is present in particular a version that uses the structure of DICE 2023 and the regional distribution of initial values (production, emissions, capital, population...) of RICE2020.

DICEParameters

Row and computed parameters for the optimization function.

This structure contains the "default of the default" parameters, which can eventually be modified using either keyword arguments in specific functions (each defining its own "defaults") or directly in the run_dice(pars) function (e.g. run_dice(a2base = [0.01])). This second method overrides the specific defaults of the DICE2013 function.

The structure first defines some "raw" parameters, and then some "computed" parameters (mostly matrices of ntsteps x nregions length). Both can be overridden with keyword arguments. In particular, "computed" parameters can be overridden in two ways: either by overriding the raw parameters from which they are computed, or by computing the parameter in a different way (outside the model) and overriding the computed parameter.

Available parameters:

• tstep: Years per period

• ntsteps: Number of time periods

• regions: Name of the regions

• weights: Utility weights to assign to each region. Note these are exogenous (default to equal weights), are NOT the Negishi weights.

• gamma: Capital elasticity in production function

• pop1: Initial world population 2020 (millions)

• popadj: Growth rate to calibrate to 2050 population projection

• popasym: Asymptotic population (millions)

• dk: Depreciation rate on capital (per year)

• q1: Initial world output 2020 (trill 2019 USD)

• al1: Initial level of total factor productivity

• ga1: Initial growth rate for TFP per 5 years

• dela: Decline rate of TFP per 5 years

• gsigma1: Initial growth of sigma (per year)

• delgsig: Decline rate of gsigma per period

• asymgsig: Asymptotic sigma

• e1: Industrial emissions 2020 (GtCO2 per year)

• miu1: Emissions control rate historical 2020

• cumemiss0: Cumulative emissions 2020 (GtC)

• a1: Damage intercept

• a2base: Damage quadratic term

• a3: Damage exponent

• expcost2: Exponent of control cost function

• pback2050: Cost of backstop in 2019$ per tCO2 (2050)

• limmiu2070: Emission control limit from 2070

• limmiu2120: Emission control limit from 2120

• limmiu2200: Emission control limit from 2220

• limmiu2300: Emission control limit from 2300

• delmiumax: Emission control delta limit per period

• betaclim: Climate beta

• elasmu: Elasticity of marginal utility of consumption

• prstp: Pure rate of social time preference

• pi_val: Capital risk premium (renamed to avoid conflict with Julia's pi)

• k0: Initial capital stock (10^12 2019 USD)

• siggc1: Annual standard deviation of consumption growth

• srf: Scaling factor for discounting

• scale1: Multiplicative scaling coefficient

• scale2: Additive scaling coefficient

• eland0: Carbon emissions from land 2015 (GtCO2 per year)

• deland: Decline rate of land emissions (per period)

• f_misc2020: Non-abatable forcings 2020

• f_misc2100: Non-abatable forcings 2100

• f_ghgabate2020: Forcings of abatable non-CO2 GHG in 2020

• eco2eghgb2020: Emissions of abatable non-CO2 GHG (GtCO2e) in 2020

• eco2eghgb2100: Emissions of abatable non-CO2 GHG (GtCO2e) in 2100

• emissrat2020: Ratio of CO2e to industrial CO2 in 2020

• emissrat2100: Ratio of CO2e to industrial CO2 in 2100

• fcoef1: Coefficient of non-CO2 abateable emissions

• fcoef2: Coefficient of non-CO2 abateable emissions

• yr0: Calendar year that corresponds to model year zero

• emshare0: Carbon emissions share into Reservoir 0

• emshare1: Carbon emissions share into Reservoir 1

• emshare2: Carbon emissions share into Reservoir 2

• emshare3: Carbon emissions share into Reservoir 3

• tau0: Decay time constant for Reservoir 0

• tau1: Decay time constant for Reservoir 1

• tau2: Decay time constant for Reservoir 2

• tau3: Decay time constant for Reservoir 3

• teq1: Thermal equilibration parameter for box 1

• teq2: Thermal equilibration parameter for box 2

• d1: Thermal response timescale for deep ocean

• d2: Thermal response timescale for upper ocean

• irf0: Pre-industrial IRF100

• irc: Increase in IRF100 with cumulative carbon uptake

• irt: Increase in IRF100 with warming

• fco22x: Forcings of equilibrium CO2 doubling

• mat0: Initial concentration in atmosphere in 2020 (GtC)

• res00: Initial concentration in Reservoir 0 in 2020 (GtC)

• res10: Initial concentration in Reservoir 1 in 2020 (GtC)

• res20: Initial concentration in Reservoir 2 in 2020 (GtC)

• res30: Initial concentration in Reservoir 3 in 2020 (GtC)

• mateq: Equilibrium concentration in atmosphere (GtC)

• tbox10: Initial temperature box 1 change in 2020 (°C)

• tbox20: Initial temperature box 2 change in 2020 (°C)

• tatm0: Initial atmospheric temperature change in 2020 (°C)

• times: Time periods sequence (0,5,10,...,400)

• tidx: Time periods index sequence (1,2,3,...,81)

• t0idx: Time periods index sequence (0,1,2,...,80)

• nreg: Number of regions

• ridx: Regions index sequence

• rartp: Risk-adjusted rate of time preference

• miuup: Upper bounds on miu

• varpcc: Variance of per capita consumption

• rprecaut: Precautionary rate of return

• rr1: STP factor without precautionary factor

• rr: STP with precautionary factor

• optlrsav: Optimal long-run savings rate used for transversality

• l_temp: Level of population and labor (temp, used only for its first value)

• l: Level of population and labor

• ga: Growth rate of Total Factor Productivity

• al_temp: Level of total factor productivity (temp, used only for its first value)

• al: Level of total factor productivity

• pbacktime: Backstop price 2019$ per ton CO2

• sig1: Carbon intensity 2020 kgCO2-output 2020

• gsig: Change in sigma (rate of decarbonization)

• sigma_temp: CO2-emissions output ratio (temp, used only for its first value)

• sigma: CO2-emissions output ratio

• eland

• co2e_ghgabateb

• f_misc

• emissrat

• sigmatot

• cost1tot

DICE2023()

Parameters constructor with defaults aligned to DICE2023.

Create a parameters struct with the defaults of the DICE2023 model. Different parameters, either the "raw" ones or the "computed ones", can be specified using keywork arguments.

See DICEParameters for a complete list of available parameters.

Example

alt_dam_scen = DICE2023(a2base = [0.01]) # single value array parameter for the "world region"

Note

• Even if DICE2023 treats the World as a single region, the model works with a "regional" dimension, so all parameters (except those related to the carbon cycle) should be entered as a single-value array for static data or a (ntime_periods by nregions) matrix for (computed) dynamic parameters.

DICE2023_NREG(n;kwargs...)

Build parameters for a DICE2023 world partitioned in n equal regions (unless parameters are overrided)

Create a parameters struct where the world is partitioned in n regions, and where each region is equal, with the same coefficients but with 1/n of initial emissions, capital, production, population, ...

Data is from DICE2013, so the output of DICE2023_NREG(1) is the same as DICE2023().

Using the keyword arguments one can specify individual parameters that differ from DICE2023, with eventually a regional specification.

See DICEParameters for a complete list of available parameters and run_dice to run the optimization with the parameter struct created with this function.

Example

four_regions_parameters = DICE2023_NREG(4) 
alt_dam_parameters      = DICE2023_NREG(2, a2base = [0.01, 0.02]) # two regions differing only for the a2base parameter
results = run_dice(four_regions_parameters)

RICE2023(;kwargs...)

Build parameters calibrated to have the DICE2023 world totals and the 12-regions RICE2020 regional distribution.

Create a parameters struct where the world is partitioned in 12 regions, with coefficients of DICE2023 but a regional variance derived in most cases from RICE2020 (initial emissions, capital, production, population) or assumed (e.g. damages)

Note that the utility weights to provide to each region are exogenous (default to equal weights), they are NOT the Negishi weights. The run_dice function can eventually be used iteractively to look for these weights.

See DICEParameters for a complete list of available parameters and run_dice to run the optimization with the parameter struct created with this function.

Example

w_rich  = [5,4,3,3,1,1,3,2,2,1,1.5,1] #
w_equal = fill(1,12)
res_cbopt_12r_poor = run_dice(RICE2023(;weights=w_poor))
res_cbopt_12r_rich = run_dice(RICE2023(;weights=w_rich)) 

Notes

• The default 12 regions are: ["USA", "EUS", "JPN", "OHI", "RUS", "EEC", "CHN", "IND", "MDE", "SSA", "LAA", "ROW"]

run_dice(pars;optimizer,bounds)
run_dice(;optimizer,bounds,kwargs...)

Run the D(R)ICE models (currently with the structure and, by default, the data of DICE2023), possibly with custom optimiser, bounds or parameters.

This function runs the DICE model and returns the results as a named tuple. Note that starting from DICEModel v0.2, a regional dimension is always present, and DICE is simply treated as RICE with a single region.

Function arguments

Positional:

• pars: An istance of the DICEParameters struct containing the needed parameters

Keyword arguments:

• optimizer: The optimiser to use and possibly its options. Defaults to: [optimizer = optimizer_with_attributes(Ipopt.Optimizer,"print_level" => 0, "max_wall_time"=>10.0^20, "max_cpu_time" => 10.0^20, "max_iter" => 3000, "acceptable_tol" =>10^-6, "acceptable_iter" => 15, "acceptable_dual_inf_tol" =>10.0^10, "acceptable_constr_viol_tol" => 0.01, "acceptable_compl_inf_tol" =>0.01, "acceptable_obj_change_tol" =>10.0^20)]. All, except the print levels, are the Ipopt defaults.
• bounds: A dictionary of equality or inequality constraints. Each constraint should be specified with the variable name as key and a two-element tuple as value. The first element is either "<=", ">=" or "==", and the second element is the right-hand side of the constraint (a single value, a vector of ntimesteps length or a matrix of ntimesteps x nregions). Default: (empty dictionary). See the source code for the names of the model variables.
• kwargs: Keyword arguments to override the default parameter values of DICE2023. See the documentation for the DICEParameters structure for the available model parameters.

WARNING: Sometimes changing a parameter doesn't lead to the expected behavior. This is because the model (in its original GAMS form that has been re-implemented in this package) performs some calibrations with the parameters, so several parameters have to be changed together. For example all the scenarios that test different discount rates don't change only the prstp parameter, but compute several other parameters, sometimes in a different matter than the default model, and have different calibration for initial conditions. Always check the source code to make sure that the parameter you want to change doesn't have other side effects in the model.

Outputs

• A named tuple containing the following fields: solved, status, times, tidx, the post_process computed values, the optimisation variables, the parameters structure (pars).

Examples:

res = run_dice()
ECO2_opt = res.ECO2
plot(res.times[1:11] .+ 2020,ECO2_opt[1:11],ylim=(0,80), title="CO₂ emissions",ylabel="GtCO₂/yr",label="C/B optimal", markershape=:circle, markercolor=:white)

res_crazy = run_dice(optimizer=optimizer_with_attributes(Ipopt.Optimizer,"print_level" => 0), bounds = Dict("MIU"=>("==",1.0), "TATM"=>("<=",15), "Y" =>(">=",[fill(floatmin(Float64),10);fill(0.1,71)]), "ECO2" =>("<=",10000)), a2base = 0.01)

w_dev_country_priority = [1,1,1,1.5,2,2,3,3,2,5,5,5] # Utility weights by region
res_12regions_devprior = run_dice(RICE2023(;weights=w_dev_country_priority))

Notes

• The bounds add constraints to the problem, but do not replace hard written bounds in the model. In particular, the miuup parameter should be used instead for the upper limit of emission controls.
• Bounds are always intended for the full time steps. If you need a bound for a subset of time steps (e.g. the first time step), you still need to assemble your full time array of the bound using floatmin(Float64) or floatmax(Float64) as appropriate.
• The version with keywords arguments is a tiny wrapper (calls) the version with the parameters struct, itself built using the DICE2023 function with the provided keyword arguments.
• The weights used in the DICEParameters struct are exogenous, are NOT the Negishi weights. The run_dice function can eventually be used iteractively to look for these weights.

run_dice_scenario(scenario::String)

Run one of the "official" 11 scenarios in Nordhous's DICE 2023 model:

• cbopt: The C/B optimal scenario
• t2c: The temperature constrained to max 2 °C scenario
• t15c: The temperature constrained to max 1.5 °C scenario
• altdam: The alternative damage scenario
• parisext: The Paris extended scenario
• base: The base (current policies) scenario
• r5: The scenario with discount rate of 5%
• r4: The scenario with discount rate of 4%
• r3: The scenario with discount rate of 3%
• r2: The scenario with discount rate of 2%
• r1: The scenario with discount rate of 1%

To run "your own" scenarios, use the function run_dice, where you can set the input parameters and constraints as you like.

Output

• A named tuple containing the following fields: solved, status, times, tidx, the post_process computed values and the optimization variables.

scaleweights(w)

Scale a vector of weights such that their sum is 1

@fields_to_vars(t,x)

Utility macro to convert struct fields to local variables (for readibility, so that we can write parameterx instead of using everywhere p.parameterx).