Merge pull request #200 from anthofflab/NPV-SCC

Change SCC calculation to return the net present value in the emissions year

src/components/EquityWeighting.jl

@@ -42,7 +42,7 @@
     yp_yearsperiod = Variable(index=[time], unit="year") # defined differently from yagg
     dfc_consumptiondiscountrate = Variable(index=[time, region], unit="1/year")
 
-    df_utilitydiscountrate = Variable(index=[time], unit="fraction")
+    df_utilitydiscountfactor = Variable(index=[time], unit="fraction")
 
     # Discounted costs
     pcdt_partiallyweighted_discounted = Variable(index=[time, region], unit="\$million")
@@ -81,7 +81,7 @@
             v.te_totaleffect = 0
         end
 
-        v.df_utilitydiscountrate[tt] = (1 + p.ptp_timepreference / 100)^(-(p.y_year[tt] - p.y_year_0))
+        v.df_utilitydiscountfactor[tt] = (1 + p.ptp_timepreference / 100)^(-(p.y_year[tt] - p.y_year_0))
 
         for rr in d.region
 
@@ -123,16 +123,16 @@
                 v.pcdt_partiallyweighted_discounted[tt, rr] = v.pct_partiallyweighted[tt, rr] * v.dfc_consumptiondiscountrate[tt, rr]
                 v.wacdt_partiallyweighted_discounted[tt, rr] = p.act_adaptationcosts_total[tt, rr] * v.dfc_consumptiondiscountrate[tt, rr]
             else
-                v.pcdt_partiallyweighted_discounted[tt, rr] = v.pct_partiallyweighted[tt, rr] * v.df_utilitydiscountrate[tt]
-                v.wacdt_partiallyweighted_discounted[tt, rr] = v.wact_partiallyweighted[tt, rr] * v.df_utilitydiscountrate[tt]
+                v.pcdt_partiallyweighted_discounted[tt, rr] = v.pct_partiallyweighted[tt, rr] * v.df_utilitydiscountfactor[tt]
+                v.wacdt_partiallyweighted_discounted[tt, rr] = v.wact_partiallyweighted[tt, rr] * v.df_utilitydiscountfactor[tt]
             end
 
             v.pcdat_partiallyweighted_discountedaggregated[tt, rr] = v.pcdt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]
 
             ## Equity weighted impacts (end of page 28, Hope 2009)
             v.wit_equityweightedimpact[tt, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_aftercosts[tt, rr]^(1 - p.emuc_utilityconvexity) - p.rcons_percap_dis[tt, rr]^(1 - p.emuc_utilityconvexity)) * p.pop_population[tt, rr]
 
-            v.widt_equityweightedimpact_discounted[tt, rr] = v.wit_equityweightedimpact[tt, rr] * v.df_utilitydiscountrate[tt]
+            v.widt_equityweightedimpact_discounted[tt, rr] = v.wit_equityweightedimpact[tt, rr] * v.df_utilitydiscountfactor[tt]
 
             v.addt_equityweightedimpact_discountedaggregated[tt, rr] = v.widt_equityweightedimpact_discounted[tt, rr] * p.yagg_periodspan[tt]
             v.aact_equityweightedadaptation_discountedaggregated[tt, rr] = v.wacdt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]

src/compute_scc.jl

@@ -27,12 +27,30 @@ function getperiodlength(year)      # same calculations made for yagg_periodspan
     return (last_year - start_year) / 2
 end
 
-@defcomp PAGE_marginal_emissions begin 
-    er_CO2emissionsgrowth = Variable(index=[time,region], unit = "%")
-    marginal_emissions_growth = Parameter(index=[time,region], unit = "%", default = zeros(10,8))
+"""
+Applies undiscounting factor to get the SCC, discounted to the emissions year instead of the base year.
+"""
+function undiscount_scc(m::Model, year::Int)
+    df = m[:EquityWeighting, :df_utilitydiscountfactor]
+    consfocus0 = m[:GDP, :cons_percap_consumption_0][1]
+    consfocus = m[:GDP, :cons_percap_consumption][:, 1]
+    emuc = m[:EquityWeighting, :emuc_utilityconvexity]
+    sccii = getpageindexfromyear(year)
+
+    return df[sccii] * ((consfocus[sccii] / consfocus0)^-emuc)
+end
+
+@defcomp ExtraEmissions begin
+    e_globalCO2emissions = Parameter(index=[time],unit="Mtonne/year")
+    pulse_size = Parameter()
+    pulse_year = Parameter()
+    e_globalCO2emissions_adjusted = Variable(index=[time],unit="Mtonne/year")
+
     function run_timestep(p, v, d, t)
-        if is_first(t)
-            v.er_CO2emissionsgrowth[:, :] = p.marginal_emissions_growth[:, :]
+        if gettime(t) == p.pulse_year
+            v.e_globalCO2emissions_adjusted[t] = p.e_globalCO2emissions[t] + p.pulse_size / getperiodlength(p.pulse_year)
+        else
+            v.e_globalCO2emissions_adjusted[t] = p.e_globalCO2emissions[t]
         end
     end
 end
@@ -100,15 +118,15 @@ function compute_scc(
     if n===nothing
         # Run the "best guess" social cost calculation
         run(mm)
-        scc = mm[:EquityWeighting, :td_totaldiscountedimpacts]
+        scc = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
     elseif n<1
         error("Invalid `n` value, only values >=1 allowed.")
     else
         # Run a Monte Carlo simulation
         simdef = getsim()
         seed !== nothing ? Random.seed!(seed) : nothing
         si = run(simdef, mm, n, trials_output_filename = trials_output_filename)
-        scc = getdataframe(si, :EquityWeighting, :td_totaldiscountedimpacts).td_totaldiscountedimpacts
+        scc = getdataframe(si, :EquityWeighting, :td_totaldiscountedimpacts).td_totaldiscountedimpacts ./ undiscount_scc(mm.base, year)
     end
 
     return scc
@@ -144,13 +162,13 @@ function compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = noth
     end
 
     mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
-    scc = mm[:EquityWeighting, :td_totaldiscountedimpacts]
+    scc = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
 
     return (scc = scc, mm = mm)
 end
 
 """
-get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = nothing)
+    get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = nothing)
 
 Returns a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of CO2 in year `year`.
 If no Model m is provided, the default model from MimiPAGE2009.get_model() is used as the base model.
@@ -162,27 +180,13 @@ function get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} =
 
     mm = create_marginal_model(m, pulse_size)
 
-    add_comp!(mm.modified, PAGE_marginal_emissions, :marginal_emissions; before = :co2emissions)
-    connect_param!(mm.modified, :co2emissions=>:er_CO2emissionsgrowth, :marginal_emissions=>:er_CO2emissionsgrowth)
-    connect_param!(mm.modified, :AbatementCostsCO2=>:er_emissionsgrowth, :marginal_emissions=>:er_CO2emissionsgrowth)
+    add_comp!(mm.modified, ExtraEmissions, :extra_emissions; after=:co2emissions)
+    connect_param!(mm.modified, :extra_emissions => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions)
+    set_param!(mm.modified, :extra_emissions, :pulse_size, pulse_size)
+    set_param!(mm.modified, :extra_emissions, :pulse_year, year)
 
-    i = getpageindexfromyear(year) 
-
-    # Base model
-    run(mm.base)
-    base_glob0_emissions = mm.base[:co2cycle, :e0_globalCO2emissions]
-    er_co2_a = mm.base[:co2emissions, :er_CO2emissionsgrowth][i, :]
-    e_co2_g = mm.base[:co2emissions, :e_globalCO2emissions]  
-
-    # Calculate pulse 
-    ER_SCC = 100 * -1 * pulse_size / (base_glob0_emissions * getperiodlength(year))
-    pulse = er_co2_a - ER_SCC * (er_co2_a/100) * (base_glob0_emissions / e_co2_g[i])
-    marginal_emissions_growth = copy(mm.base[:co2emissions, :er_CO2emissionsgrowth])
-    marginal_emissions_growth[i, :] = pulse
-
-    # Marginal emissions model
-    set_param!(mm.modified, :marginal_emissions_growth, marginal_emissions_growth)
-    run(mm.modified)
+    connect_param!(mm.modified, :co2cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted)
+    run(mm)
 
     return mm
-end
+end

test/test_EquityWeighting.jl

@@ -24,7 +24,7 @@ set_leftover_params!(m, p)
 run(m)
 
 # Generated data
-df = m[:EquityWeighting, :df_utilitydiscountrate]
+df = m[:EquityWeighting, :df_utilitydiscountfactor]
 wtct_percap = m[:EquityWeighting, :wtct_percap_weightedcosts]
 pct_percap = m[:EquityWeighting, :pct_percap_partiallyweighted]
 dr = m[:EquityWeighting, :dr_discountrate]
@@ -40,7 +40,7 @@ addt_gt = m[:EquityWeighting, :addt_gt_equityweightedimpact_discountedglobal]
 te = m[:EquityWeighting, :te_totaleffect]
 
 # Recorded data
-df_compare = readpagedata(m, "test/validationdata/df_utilitydiscountrate.csv")
+df_compare = readpagedata(m, "test/validationdata/df_utilitydiscountfactor.csv")
 wtct_percap_compare = readpagedata(m, "test/validationdata/wtct_percap_weightedcosts.csv")
 pct_percap_compare = readpagedata(m, "test/validationdata/pct_percap_partiallyweighted.csv")
 dr_compare = readpagedata(m, "test/validationdata/dr_discountrate.csv")

test/test_standard_api.jl

@@ -23,7 +23,7 @@ mm[:ClimateTemperature, :rt_realizedtemperature]
 
 # Test compute_scc_mm
 result = MimiPAGE2009.compute_scc_mm(year=2050)
-@test result.scc < scc1
+@test result.scc > scc1
 @test result.mm isa Mimi.MarginalModel
 
 # Test Monte Carlo SCC support