Simplified DICE model

A simplified Integrated Assesment Model calibration and projection.

Simplified DICE with Exogenous GDP projections

Material based on “Teaching Options for IAMs” by Barrage, Greenstone, Metcalf in AEA Continuing Education: Climate Change - January 2020

Check the interactive shiny app here

Key equations

1. We assume a particular GDP path: \(Y_{2000}, Y_{2010}, Y_{2020},...\)

2. Production causes carbon emissions \(E_t\) based on carbon intensity \(\sigma_t\)

\[ E_t=\sigma_t Y_t \]

3. Society can reduce carbon emissions by fraction \(\mu_t\)

\[E_{t}^{Net}=(1-\mu_t)E_t\]

4. Fraction \(\xi\) of emissions add to atmospheric carbon stock \(S_t\)

\[S_t=S_{2000}+\xi \sum_{j=2000}^t E_j^{Net}\]

4. Carbon concentrations lead to temperature change \(T_t\) (change since 1900)

\[T_t \approx T_{2000}+\eta \log_2\left( \frac{S_t}{S_{2000}}\right)\]

6. Damage Function \(D(T_t)\) due to climate change \(T_t\)

\[D(T_t)=\left(\phi_1 T_t^{\phi_2}\right) Y_t\]

7. Total mitigation costs \(MC_t\)

\[MC_t=\left(\theta_{1t}\mu_t^{\theta_2}\right) Y_t\]

8. Economy’s interest rate \(r_t\) follows Ramsey rule

\[r_t\approx \rho + \alpha g_t^*\]

9. Social Planner’s problem: choose Mitigation to minimize present value sum of climate damages and mitigation costs

\[\min_{\{\mu_t\}_{t=0}^{2250}} \sum_{j=2000}^{2500} \frac{D(T_t)+MC_t(\mu_t)}{(1+r_t)^{j-2000}} \]

Summary of the model

Exogenous series

• \(Y_t\): GDP. The series is constructed by assuming a growth rate for the GDP per Capita (\(g^{*}_t\)) and Population (\(g^{N}_t\)).
• \(r_t\): Market discount rate. Follows Ramsey Rule, equation 8.
• \(\sigma_t\): Carbon intensity. Series is constructed assuming a constant growth rate, \(g^{\sigma}\).
• \(\theta_{1t}\): Scale parameter of the total Mitigation costs function. Series is constructed assuming a constant growth rate, \(g^{\theta}\).

Endogenous series

• \(E_t\): Gross Carbon Emissions. Follows equation 2.
• \(E_t^{Net}\): Net of Mitigation Carbon Emissions. Follows equation 3.
• \(S_t\): Atmospheric CO2 Stock. Follows equation 4.
• \(T_t\): Equilibrium Temperature change. Follows equation 5.
• \(D(T_t)\): Damages costs due to climate change. Follows equation 6.
• \(MC_t\): Total Mitigation costs. Follows equation 7.

Choice Variable

• \(\mu_t\): Mitigation policy. Taken as given for initial projections then optimized in the Social Planner’s problem according to equation 9.

Calibrated Parameters

Parameter	Notation	Units	Source	Value
Baseline global GDP	\(Y_{2000}\)	(trillions of 2000-US$/year)	World Bank	32.21
Baseline global population	\(N_{2000}\)	(billions)	World Bank	6.08
Population growth rate until 2100	\(g^{N}_t\)	(% per year)	Assumed	1%
Population growth rate after 2100	\(g^{N}_t\)	(% per year)	Assumed	0%
Consumption (GDP per capita) growth rate, until 2050 (inclusive),	\(g^*_t\)	(% per year)	Assumed	2.75%
Consumption (GDP per capita) growth rate after 2050 (i.e., starting in 2060)	\(g^*_t\)	(% per year)	Assumed	2.00%
Base year (2000) CO2 emissions	\(E_{2000}\)	(billion tons of carbon)	World Bank	6.73
Base year (2000) CO2 emissions intensity	\(\sigma_{2000}\)		Calculated	0.21
CO2 emissions intensity growth rate	\(g^{\sigma}\)	(% per year)	Assumed	-2%
Base year (2000) CO2 concentrations	\(S_{2000}\)	(billion tons of carbon)	Nordhaus (2011) (approx)	800
Fraction of CO2 emissions remaining in atmosphere	\(\xi\)		Assumed	0.50
Temperature increase since 1900	\(T_{2000}\)	(degrees Celsius)	Nordhaus (2011) (approx)	1
Climate sensitivity parameter	\(\eta\)		Nordhaus (2011) (approx)	3.00
Damage function scale parameter	\(\phi_1\)		Assumed	0.0033
Damage function curvature parameter	\(\phi_2\)		Assumed	2.00
Mitigation cost function scale parameter	\(\theta_1\)		Assumed	0.1
Mitigation cost function curvature parameter	\(\theta_2\)		Assumed	3.3
Mitigation cost scale parameter growth rate	\(g^{\theta}\)	(% per year)	Assumed	-2.00%
Utility discount rate (= pure rate of social time preference) per year	\(\rho\)		Assumed	1.50%
Elasticity of marginal utility of consumption	\(\alpha\)		Assumed	1.50

Projecting Variables

The first step is to define calibrated parameters in the code.

(Click ‘Show Code’ to the see the source code)

Show Code

#Parameters
TimeSqn=seq(2000,2500,by=10)  # Projection Time range
Y2000=32.209312564935         # Baseline global GDP (trillions of constant 2000 US$/year)
E2000=6.73                    # Base year (2000) CO2 emissions (billion tons of carbon)
T2000=1                       # Temperature increase since 1900 (degrees Celsius)
S2000=800                     # Base year (2000) CO2 concentrations (billion tons of carbon)
Pop2000=6.084959036           # Base year (2000) global population (billions)
sigma2000=E2000/Y2000         # Base year (2000) CO2 emis intensity (bill tons of carbon/tril. US$)
sigma_gr=-0.02                # CO2 emissions intensity growth rate (per year)
xi=0.5                        # Fraction of CO2 emissions remaining in atmosphere
eta=3                         # Climate sensitivity parameter
phi1=0.0033                   # Damage function scale parameter
phi2=2                        # Damage function curvature parameter
pop_gr2100=0.01               # Population growth rate (per year) until 2100
pop_gr2500=0.00               # Population growth rate (per year) after 2100
consumption_gr2050=0.0275     # Consumption (GDP per capita) growth rate (per year), until 2050
consumption_gr2500=0.02       # Consumption (GDP per capita) growth rate (per year) after 2050
rho=0.015                     # Utility discount rate (= pure rate of social time preference) per year
alpha=1.5                     # Elasticity of marginal utility of consumption
theta1=0.1                    # Mitigation cost function scale parameter
theta2=3.3                    # Mitigation cost function curvature parameter 
Mitigation_gr=-0.02            # Mitigation cost scale parameter growth rate

Economic Variables

Assuming a GDP per capita growth rate and a population growth rate we can project a series for the GDP, \(Y_t\).

Variable description and units:

• GDP: (trillions of constant 2000 USD/year), \(Y_t\)
• Population: (billions, beginning of year), \(N_t\)
• GDP_Pop: GDP/Population (thousands USD/year), \(Y_t/N_t\)
• DiscRate: Market Discount Rate (percent per year), \(r_t\)

Show Code

#Dataframe of Economic Variable
EconomicDF=data.frame(Year=TimeSqn, GDP=0, Population=0, GDP_Pop=0, DiscRate=0)
#Initial Values
EconomicDF$GDP[1]=Y2000; EconomicDF$Population[1]=Pop2000; 
EconomicDF$GDP_Pop[1]=Y2000/Pop2000; EconomicDF$DiscRate[1]=rho+alpha*consumption_gr2050;
#Evolution of Variables
for(t in 2:51){
    EconomicDF$GDP_Pop[t]=EconomicDF$GDP_Pop[t-1]*(1+ifelse(EconomicDF$Year[t]<=2050,consumption_gr2050,consumption_gr2500))^10
    EconomicDF$Population[t]=EconomicDF$Population[t-1]*(1+ifelse(EconomicDF$Year[t]<=2100,pop_gr2100,pop_gr2500))^10
    EconomicDF$GDP[t]=EconomicDF$Population[t]*EconomicDF$GDP_Pop[t]
    EconomicDF$DiscRate[t]=rho+alpha*ifelse(EconomicDF$Year[t]<=2050,consumption_gr2050,consumption_gr2500);
}

(Click on the series name to show/hide series plot) (Notice the units in each series)

Show Code

library(dplyr); library(ggplot2); library(highcharter);
#Plot
EconomicDF %>% mutate(across(where(is.numeric), round, 3)) %>% 
  reshape2::melt(id.vars="Year") %>% 
  hchart("line", hcaes(y=value, x=Year, group=variable)) %>%
  hc_add_theme(hc_theme_538()) %>% 
  hc_yAxis(title = list(text = "Level of Economic Variables")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Projected Economic Variables") %>% 
  hc_plotOptions(line=list("visible"=FALSE))

Emissions Series

*With constant mitigation series: \(\mu_t=0.1\) for all \(t\). This will be optimize below.

Variable description and units:

• Intensity: CO2 Emissions Intensity [(bil. MT carbon/year) / (tril. 2000 USD/year)], \(\sigma_t\)
• Gross_Emissions: CO2 Emissions (billion metric tons of carbon/year), \(E_t\)
• Net_Emissions: Net of Mitigation CO2 Emissions (bil. MT carbon/year), \(E_t^{Net}\)
• Atmospheric_CO2: Atmospheric CO2 Concentration (bil. Tons carbon, end of decade), \(S_t\)
• EqTemp_Change: Equilibrium Temperature Change (°C), \(T_t\)
• Mitigation_Policy: (percentage of emissions reduction), \(\mu_t\)

Show Code

#Dataframe of Emissions Series
EmissionsDF=data.frame(Year=TimeSqn, Intensity=0, Gross_Emissions=0, Net_Emissions=0, Atmospheric_CO2=0, EqTemp_Change=0, Mitigation_Policy=0.1)
#Initial Values
EmissionsDF$Intensity[1]=sigma2000
EmissionsDF$Gross_Emissions[1]=E2000
EmissionsDF$Net_Emissions[1]=(1-EmissionsDF$Mitigation_Policy[1])*EmissionsDF$Gross_Emissions[1]
EmissionsDF$Atmospheric_CO2[1]=S2000+xi*EmissionsDF$Net_Emissions[1]*10
EmissionsDF$EqTemp_Change[1]=T2000+eta*log(EmissionsDF$Atmospheric_CO2[1]/S2000)/log(2)
#Evolution of Variables
for(t in 2:51){
  EmissionsDF$Intensity[t]=EmissionsDF$Intensity[t-1]*(1+sigma_gr)^10
  EmissionsDF$Gross_Emissions[t]=EmissionsDF$Intensity[t]*EconomicDF$GDP[t]
  EmissionsDF$Net_Emissions[t]=(1-EmissionsDF$Mitigation_Policy[t])*EmissionsDF$Gross_Emissions[t]
  EmissionsDF$Atmospheric_CO2[t]=EmissionsDF$Atmospheric_CO2[t-1]+xi*EmissionsDF$Net_Emissions[t]*10
  EmissionsDF$EqTemp_Change[t]=T2000+eta*log(EmissionsDF$Atmospheric_CO2[t]/S2000)/log(2)
}

Show Code

#Plot
EmissionsDF %>% mutate(across(where(is.numeric), round, 3)) %>% 
  reshape2::melt(id.vars="Year") %>% 
  hchart("line", hcaes(y=value, x=Year, group=variable)) %>%
  hc_add_theme(hc_theme_538()) %>% 
  hc_yAxis(title = list(text = "Level of Emissions Variables")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Projected Emissions Variables") %>% 
  hc_plotOptions(line=list("visible"=FALSE))

Damage Costs Projections

Variable description and units:

• Damages_GDP: Damages/GDP (percent per year)
• Damages: (trillions of constant 2000 US$/year)
• PV_Damages: PV of damages (trillions of constant 2000 US$/decade)

Show Code

#Dataframe of Damages Projections
DamagesDF=data.frame(Year=TimeSqn, Damages_GDP=0, Damages=0, PV_Damages=0)

#Evolution of Variables
for(t in 1:51){
  DamagesDF$Damages_GDP[t]=phi1*(EmissionsDF$EqTemp_Change[t]^phi2)
  DamagesDF$Damages[t]=DamagesDF$Damages_GDP[t]*EconomicDF$GDP[t]
  DamagesDF$PV_Damages[t]=10*DamagesDF$Damages[t]/(1+EconomicDF$DiscRate[t])^((t-1)*10)
}

Show Code

#Plot
DamagesDF %>% mutate(across(where(is.numeric), round, 3)) %>% 
  reshape2::melt(id.vars="Year") %>% 
  hchart("line", hcaes(y=value, x=Year, group=variable)) %>%
  hc_add_theme(hc_theme_538()) %>% 
  hc_yAxis(title = list(text = "Level of Damages Variables")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Projected Damages") %>% 
  hc_plotOptions(line=list("visible"=FALSE))

Mitigation Cost Projections

Variable description and units:

• Mitigation_GDP: Total Mitigation costs / GDP (percent per period)
• Mitigation: Total Mitigation costs (tril. 2000 US$/decade)
• PV_Mitigation: PV of total Mitigation costs (tril. 2000 US$/decade)

Show Code

#Dataframe for Mitigation Series
MitigationDF=data.frame(Year=TimeSqn, Mitigation_GDP=0, Mitigation=0, PV_Mitigation=0)
for(t in 1:51){
  MitigationDF$Mitigation_GDP[t]=(theta1*(1+Mitigation_gr)^(10*(t-1)))*EmissionsDF$Mitigation_Policy[t]^theta2
  MitigationDF$Mitigation[t]=MitigationDF$Mitigation_GDP[t]*EconomicDF$GDP[t]*10
  MitigationDF$PV_Mitigation[t]=MitigationDF$Mitigation[t]/(1+EconomicDF$DiscRate[t])^((t-1)*10)
}

Show Code

#Plot
MitigationDF %>% 
  reshape2::melt(id.vars="Year") %>% 
  hchart("line", hcaes(y=value, x=Year, group=variable)) %>%
  hc_add_theme(hc_theme_538()) %>% 
  hc_yAxis(title = list(text = "Level of Mitigation Variables")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Projected Mitigation Series") %>% 
  hc_plotOptions(line=list("visible"=FALSE))

Solving for Mitigation policy

1. Wrap previous projections in a functions of the Mitigation series.
2. Optimize the objective function

Show Code

#Preliminary Function to retrieve series
projectSeries=function(mu, rho, alpha){
  #Parameters
  TimeSqn=seq(2000,2500,by=10)  # Projection Time range
  Y2000=32.209312564935         # Baseline global GDP (trillions of constant 2000 US$/year)
  E2000=6.73                    # Base year (2000) CO2 emissions (billion tons of carbon)
  T2000=1                       # Temperature increase since 1900 (degrees Celsius)
  S2000=800                     # Base year (2000) CO2 concentrations (billion tons of carbon)
  Pop2000=6.084959036           # Base year (2000) global population (billions)
  sigma2000=E2000/Y2000         # Base year (2000) CO2 emis intensity (bill tons of carbon/tril. US$)
  sigma_gr=-0.02                # CO2 emissions intensity growth rate (per year)
  xi=0.5                        # Fraction of CO2 emissions remaining in atmosphere
  eta=3                         # Climate sensitivity parameter
  phi1=0.0033                   # Damage function scale parameter
  phi2=2                        # Damage function curvature parameter
  pop_gr2100=0.01               # Population growth rate (per year) until 2100
  pop_gr2500=0.00               # Population growth rate (per year) after 2100
  consumption_gr2050=0.0275     # Consumption (GDP per capita) growth rate (per year), until 2050
  consumption_gr2500=0.02       # Consumption (GDP per capita) growth rate (per year) after 2050
  theta1=0.1                    # Mitigation cost function scale parameter
  theta2=3.3                    # Mitigation cost function curvature parameter 
  Mitigation_gr=-0.02
  #Dataframes
  EconomicDF=data.frame(Year=TimeSqn, GDP=0, Population=0, GDP_Pop=0, DiscRate=0)
  EmissionsDF=data.frame(Year=TimeSqn, Intensity=0, Gross_Emissions=0, Net_Emissions=0, Atmospheric_CO2=0, EqTemp_Change=0, Mitigation_Policy=mu)
  DamagesDF=data.frame(Year=TimeSqn, Damages_GDP=0, Damages=0, PV_Damages=0)
  MitigationDF=data.frame(Year=TimeSqn, Mitigation_GDP=0, Mitigation=0, PV_Mitigation=0)
  #Initial Values
  EconomicDF$GDP[1]=Y2000;
  EconomicDF$Population[1]=Pop2000; 
  EconomicDF$GDP_Pop[1]=Y2000/Pop2000;
  EconomicDF$DiscRate[1]=rho+alpha*consumption_gr2050;
  
  EmissionsDF$Intensity[1]=sigma2000
  EmissionsDF$Gross_Emissions[1]=E2000
  EmissionsDF$Net_Emissions[1]=(1-EmissionsDF$Mitigation_Policy[1])*EmissionsDF$Gross_Emissions[1]
  EmissionsDF$Atmospheric_CO2[1]=S2000+xi*EmissionsDF$Net_Emissions[1]*10
  EmissionsDF$EqTemp_Change[1]=T2000+eta*log(EmissionsDF$Atmospheric_CO2[1]/S2000)/log(2)
  #Projecting variables
  for(t in 1:length(TimeSqn)){
    if(t>1){
      EconomicDF$GDP_Pop[t]=EconomicDF$GDP_Pop[t-1]*(1+ifelse(EconomicDF$Year[t]<=2050,consumption_gr2050,consumption_gr2500))^10
      EconomicDF$Population[t]=EconomicDF$Population[t-1]*(1+ifelse(EconomicDF$Year[t]<=2100,pop_gr2100,pop_gr2500))^10
      EconomicDF$GDP[t]=EconomicDF$Population[t]*EconomicDF$GDP_Pop[t]
      EconomicDF$DiscRate[t]=rho+alpha*ifelse(EconomicDF$Year[t]<=2050,consumption_gr2050,consumption_gr2500);
      
      EmissionsDF$Intensity[t]=EmissionsDF$Intensity[t-1]*(1+sigma_gr)^10
      EmissionsDF$Gross_Emissions[t]=EmissionsDF$Intensity[t]*EconomicDF$GDP[t]
      EmissionsDF$Net_Emissions[t]=(1-EmissionsDF$Mitigation_Policy[t])*EmissionsDF$Gross_Emissions[t]
      EmissionsDF$Atmospheric_CO2[t]=EmissionsDF$Atmospheric_CO2[t-1]+xi*EmissionsDF$Net_Emissions[t]*10
      EmissionsDF$EqTemp_Change[t]=T2000+eta*log(EmissionsDF$Atmospheric_CO2[t]/S2000)/log(2)
    }
    DamagesDF$Damages_GDP[t]=phi1*(EmissionsDF$EqTemp_Change[t]^phi2)
    DamagesDF$Damages[t]=DamagesDF$Damages_GDP[t]*EconomicDF$GDP[t]
    DamagesDF$PV_Damages[t]=10*DamagesDF$Damages[t]/(1+EconomicDF$DiscRate[t])^((t-1)*10)
    
    MitigationDF$Mitigation_GDP[t]=(theta1*(1+Mitigation_gr)^(10*(t-1)))*EmissionsDF$Mitigation[t]^theta2
    MitigationDF$Mitigation[t]=MitigationDF$Mitigation_GDP[t]*EconomicDF$GDP[t]*10
    MitigationDF$PV_Mitigation[t]=MitigationDF$Mitigation[t]/(1+EconomicDF$DiscRate[t])^((t-1)*10)
  }
  #Objective function value
  ObjValue=sum(DamagesDF$PV_Damages)+sum(MitigationDF$PV_Mitigation)

  #Implied Optimal Carbon Price
  optCarbonPrice=(EmissionsDF$Intensity^(-1))*(theta1*(1+Mitigation_gr)^(0:50*10))*theta2*(EmissionsDF$Mitigation_Policy^(theta2-1))*1000
  
  #Collect Dataframes into one
  Projections=cbind(EconomicDF,EmissionsDF[,-1],DamagesDF[,-1],MitigationDF[,-1],optCarbonPrice)
  #Return objects
  return(list("Value"=ObjValue, "Damages"=sum(DamagesDF$PV_Damages),"MitigationCosts"=sum(MitigationDF$PV_Mitigation), "Projections"=Projections))
}

#Test=projectSeries(rep(0.1,51), 0.015, 1.5)

#Create Objective Function
SocialWelfare=function(mu1){
  Res=projectSeries(mu1, 0.015, 1.5)
  return(Res$Value)
}

#Perform Optimization
optimization=optim(par=rep(0.5,51),SocialWelfare,NULL,method = "L-BFGS-B", lower = rep(0, 51), upper = rep(1, 51))

OptimalSeries=projectSeries(optimization$par, rho, alpha)$Projections

Optimal Policy

Show Code

OptimalSeries %>% dplyr::select(Year, Mitigation_Policy) %>% 
  filter(Year<=2300) %>% 
  mutate(across(where(is.numeric), round, 3)) %>% 
  hchart("line", hcaes(y=Mitigation_Policy*100, x=Year)) %>%
  hc_add_theme(hc_theme_538()) %>% 
  hc_yAxis(title = list(text = "Mitigation Policy (%)")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Optimal Mitigation Policy")

Projections under optimal policy

Click on the name of one or more series to plot them

Show Code

OptimalSeries %>% 
  mutate(across(where(is.numeric), round, 3)) %>% 
  reshape2::melt(id.vars="Year") %>% 
  hchart("line", hcaes(y=value, x=Year, group=variable)) %>%
  hc_add_theme(hc_theme_538()) %>% 
  hc_yAxis(title = list(text = "Projected Level")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Projected Series under Optimal Policy") %>% 
  hc_plotOptions(line=list("visible"=FALSE))

Alternative parameters: Stern discounting, \(\rho=0.001\)

Click on the name of one or more series to plot them

Show Code

SocialWelfare2=function(mu1){
  Res=projectSeries(mu1, 0.001, 1.5)
  return(Res$Value)
}

#Perform Optimization
optimization2=optim(par=rep(0.5,51),SocialWelfare2,NULL,method = "L-BFGS-B", lower = rep(0, 51), upper = rep(1, 51))
OptimalSeries2=projectSeries(optimization2$par, rho, alpha)$Projections

Show Code

OptimalSeries2 %>% 
  mutate(across(where(is.numeric), round, 3)) %>% 
  reshape2::melt(id.vars="Year") %>% 
  hchart("line", hcaes(y=value, x=Year, group=variable)) %>%
  hc_add_theme(hc_theme_538()) %>% 
  hc_yAxis(title = list(text = "Projected Level")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Projected Series under Optimal Policy - Stern discounting") %>% 
  hc_plotOptions(line=list("visible"=FALSE))