Merge pull request #215 from anthofflab/equity_proportion

use `equity_proportion` variable for ramsey discounting

src/components/EquityWeighting.jl

@@ -56,8 +56,8 @@
     # Equity weighted impact totals
     rcons_percap_dis = Parameter(index=[time, region], unit="\$/person")
 
-    wit_equityweightedimpact = Variable(index=[time, region], unit="\$million")
-    widt_equityweightedimpact_discounted = Variable(index=[time, region], unit="\$million")
+    wit_partiallyweighted = Variable(index=[time, region], unit="\$million")
+    widt_partiallyweighted_discounted  = Variable(index=[time, region], unit="\$million")
 
     yagg_periodspan = Parameter(index=[time], unit="year")
 
@@ -122,19 +122,20 @@
             if p.equity_proportion == 0
                 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]
+
+                v.wit_partiallyweighted[tt, rr] = (p.cons_percap_aftercosts[tt, rr] - p.rcons_percap_dis[tt, rr]) * p.pop_population[tt, rr] # equivalent to emuc = 0
+                v.widt_partiallyweighted_discounted[tt, rr] = v.wit_partiallyweighted[tt, rr] * v.dfc_consumptiondiscountrate[tt] # apply Ramsey discounting
             else
                 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]
+
+                ## Equity weighted impacts (end of page 28, Hope 2009)
+                v.wit_partiallyweighted[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_partiallyweighted_discounted[tt, rr] = v.wit_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_utilitydiscountfactor[tt]
-
-            v.addt_equityweightedimpact_discountedaggregated[tt, rr] = v.widt_equityweightedimpact_discounted[tt, rr] * p.yagg_periodspan[tt]
+            v.addt_equityweightedimpact_discountedaggregated[tt, rr] = v.widt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]
             v.aact_equityweightedadaptation_discountedaggregated[tt, rr] = v.wacdt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]
         end
 

src/compute_scc.jl

@@ -42,7 +42,7 @@ end
 
 @defcomp ExtraEmissions begin
     e_globalCO2emissions = Parameter(index=[time],unit="Mtonne/year")
-    pulse_size = Parameter()
+    pulse_size = Parameter(unit="Mtonne CO2")
     pulse_year = Parameter()
     e_globalCO2emissions_adjusted = Variable(index=[time],unit="Mtonne/year")
 
@@ -74,7 +74,7 @@ and ptp_timepreference in the model. If no values are provided, the discount fac
 PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.
 
 The size of the marginal emission pulse can be modified with the `pulse_size` keyword argument, in metric 
-tonnes (this does not change the units of the returned value, which is always normalized by the `pulse_size` used).
+tonnes of CO2 (this does not change the units of the returned value, which is always normalized by the `pulse_size` used).
 
 By default, `n = nothing`, and a single value for the "best guess" social cost of CO2 is returned. If a positive 
 value for keyword `n` is specified, then a Monte Carlo simulation with sample size `n` will run, sampling from 
@@ -88,7 +88,8 @@ function compute_scc(
         year::Union{Int, Nothing} = nothing,
         eta::Union{Float64, Nothing} = nothing,
         prtp::Union{Float64, Nothing} = nothing,
-        pulse_size = 100000.,
+        equity_weighting::Bool = true,
+        pulse_size = 100_000.,
         n::Union{Int,Nothing}=nothing,
         trials_output_filename::Union{String, Nothing} = nothing,
         seed::Union{Int, Nothing} = nothing
@@ -97,21 +98,29 @@ function compute_scc(
     year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
     !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing 
 
-    if eta != nothing
+    if eta !== nothing
         try
             set_param!(m, :emuc_utilityconvexity, eta)      # since eta is a default parameter in PAGE, we need to use `set_param!` if it hasn't been set yet
         catch e
             update_param!(m, :emuc_utilityconvexity, eta)   # or update_param! if it has been set
         end
     end
 
-    if prtp != nothing
+    if prtp !== nothing
         try
             set_param!(m, :ptp_timepreference, prtp * 100)      # since prtp is a default parameter in PAGE, we need to use `set_param!` if it hasn't been set yet
         catch e
             update_param!(m, :ptp_timepreference, prtp * 100)   # or update_param! if it has been set
         end
     end
+
+    if !equity_weighting
+        try
+            set_param!(m, :equity_proportion, 0)
+        catch e
+            update_param!(m, :equity_proportion, 0)
+        end
+    end
 
     mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
 

test/test_EquityWeighting.jl

@@ -34,7 +34,7 @@ pcdt = m[:EquityWeighting, :pcdt_partiallyweighted_discounted]
 wacdt = m[:EquityWeighting, :wacdt_partiallyweighted_discounted]
 aact = m[:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated]
 
-wit = m[:EquityWeighting, :wit_equityweightedimpact]
+wit = m[:EquityWeighting, :wit_partiallyweighted]
 addt = m[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated]
 addt_gt = m[:EquityWeighting, :addt_gt_equityweightedimpact_discountedglobal]
 te = m[:EquityWeighting, :te_totaleffect]
@@ -50,7 +50,7 @@ pcdt_compare = readpagedata(m, "test/validationdata/pcdt_partiallyweighted_disco
 wacdt_compare = readpagedata(m, "test/validationdata/wacdt_partiallyweighted_discounted.csv")
 aact_compare = readpagedata(m, "test/validationdata/aact_equityweightedadaptation_discountedaggregated.csv")
 
-wit_compare = readpagedata(m, "test/validationdata/wit_equityweightedimpact.csv")
+wit_compare = readpagedata(m, "test/validationdata/wit_partiallyweighted.csv")
 addt_compare = readpagedata(m, "test/validationdata/addt_equityweightedimpact_discountedaggregated.csv")
 addt_gt_compare = 204132238.85242900
 te_compare = 213208136.69903600

test/test_standard_api.jl

@@ -34,4 +34,19 @@ sccs2 = MimiPAGE2009.compute_scc(year=2020, n=10, seed=350)
 sccs3 = MimiPAGE2009.compute_scc(year=2020, n=10, seed=351)
 @test sccs3 != sccs1
 
+#Test equity weighting options
+
+# Test that regional equity weighting inreases the SCC
+@test MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp = 0.01, equity_weighting = true) >
+        MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp = 0.01, equity_weighting = false)
+
+# Test that when eta==0, equity_weighting does not change the SCC value
+@test MimiPAGE2009.compute_scc(year=2020, eta=0., prtp = 0.03, equity_weighting = true) ≈
+        MimiPAGE2009.compute_scc(year=2020, eta=0., prtp = 0.03, equity_weighting = false) atol=1e3
+
+# Test Monte Carlo w/ and w/o equity weighting, with the same seed
+scc10a = MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp = 0.01, equity_weighting = true, n = 10, seed=350)
+scc10b = MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp = 0.01, equity_weighting = false, n = 10, seed=350)
+@test all(scc10a .> scc10b)
+
 end