From 5a6d2ab6c5517a472c4b6bc553d4e1cfa99397c8 Mon Sep 17 00:00:00 2001 From: Cora Kingdon Date: Tue, 7 Jul 2020 11:48:54 -0400 Subject: [PATCH 1/5] change "df_utilitydiscountrate" name to "df_utilitydiscountfactor" --- src/components/EquityWeighting.jl | 10 +++++----- test/test_EquityWeighting.jl | 4 ++-- ...tydiscountrate.csv => df_utilitydiscountfactor.csv} | 0 3 files changed, 7 insertions(+), 7 deletions(-) rename test/validationdata/{df_utilitydiscountrate.csv => df_utilitydiscountfactor.csv} (100%) diff --git a/src/components/EquityWeighting.jl b/src/components/EquityWeighting.jl index 2b2e6745..64274965 100644 --- a/src/components/EquityWeighting.jl +++ b/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,8 +123,8 @@ 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] @@ -132,7 +132,7 @@ ## 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] diff --git a/test/test_EquityWeighting.jl b/test/test_EquityWeighting.jl index 148a99f0..69715b80 100644 --- a/test/test_EquityWeighting.jl +++ b/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") diff --git a/test/validationdata/df_utilitydiscountrate.csv b/test/validationdata/df_utilitydiscountfactor.csv similarity index 100% rename from test/validationdata/df_utilitydiscountrate.csv rename to test/validationdata/df_utilitydiscountfactor.csv From c2a0567b22f678f698e8777454f9ead5e2d2d794 Mon Sep 17 00:00:00 2001 From: Cora Kingdon Date: Tue, 7 Jul 2020 12:03:40 -0400 Subject: [PATCH 2/5] add undiscounting function for the SCC --- src/compute_scc.jl | 70 ++++++++++++++++++++++++---------------------- 1 file changed, 37 insertions(+), 33 deletions(-) diff --git a/src/compute_scc.jl b/src/compute_scc.jl index bd5874ee..1c8c9049 100644 --- a/src/compute_scc.jl +++ b/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 @@ -79,15 +97,15 @@ 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 - eta == nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta) - prtp == nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.) + eta === nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta) + prtp === nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.) mm = get_marginal_model(m, year=year, pulse_size=pulse_size) # Returns a marginal model that has already been run 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 @@ -95,7 +113,7 @@ function compute_scc( simdef = getsim() seed !== nothing ? Random.seed!(seed) : nothing si = run(simdef, mm, n, trials_output_filename = trials_output_filename) - scc = si[:EquityWeighting, :td_totaldiscountedimpacts].td_totaldiscountedimpacts + scc = si[:EquityWeighting, :td_totaldiscountedimpacts].td_totaldiscountedimpacts ./ undiscount_scc(mm.base, year) end return scc @@ -114,17 +132,17 @@ function compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = noth 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 - eta == nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta) - prtp == nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.) + eta === nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta) + prtp === nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.) 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. @@ -136,27 +154,13 @@ function get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = mm = create_marginal_model(m, pulse_size) - add_comp!(mm.marginal, PAGE_marginal_emissions, :marginal_emissions; before = :co2emissions) - connect_param!(mm.marginal, :co2emissions=>:er_CO2emissionsgrowth, :marginal_emissions=>:er_CO2emissionsgrowth) - connect_param!(mm.marginal, :AbatementCostsCO2=>:er_emissionsgrowth, :marginal_emissions=>:er_CO2emissionsgrowth) - - 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 + add_comp!(mm.marginal, ExtraEmissions, :extra_emissions; after=:co2emissions) + connect_param!(mm.marginal, :extra_emissions => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions) + set_param!(mm.marginal, :extra_emissions, :pulse_size, pulse_size) + set_param!(mm.marginal, :extra_emissions, :pulse_year, year) - # Marginal emissions model - update_param!(mm.marginal, :marginal_emissions_growth, marginal_emissions_growth) - run(mm.marginal) + connect_param!(mm.marginal, :CO2Cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted) + run(mm) return mm end \ No newline at end of file From 88ba7f61b57593713327780d0d45c003f6efa3b6 Mon Sep 17 00:00:00 2001 From: Cora Kingdon Date: Tue, 7 Jul 2020 12:05:59 -0400 Subject: [PATCH 3/5] update test file --- test/test_standard_api.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_standard_api.jl b/test/test_standard_api.jl index adc21cd3..bcdf5f1a 100644 --- a/test/test_standard_api.jl +++ b/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 From fa317fd0b2b76947ba4a6093e7c5ffa66447b195 Mon Sep 17 00:00:00 2001 From: Cora Kingdon Date: Tue, 7 Jul 2020 13:29:29 -0400 Subject: [PATCH 4/5] fix co2cycle typo --- src/compute_scc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compute_scc.jl b/src/compute_scc.jl index 1c8c9049..ab503307 100644 --- a/src/compute_scc.jl +++ b/src/compute_scc.jl @@ -159,7 +159,7 @@ function get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = set_param!(mm.marginal, :extra_emissions, :pulse_size, pulse_size) set_param!(mm.marginal, :extra_emissions, :pulse_year, year) - connect_param!(mm.marginal, :CO2Cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted) + connect_param!(mm.marginal, :co2cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted) run(mm) return mm From b17f7ae4d0c289a31fa1128cf899f0bb32f15bbd Mon Sep 17 00:00:00 2001 From: Cora Kingdon Date: Tue, 22 Sep 2020 16:44:25 -0400 Subject: [PATCH 5/5] use mm.modified instead of mm.marginal --- src/compute_scc.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/compute_scc.jl b/src/compute_scc.jl index ce34c833..c5e8975e 100644 --- a/src/compute_scc.jl +++ b/src/compute_scc.jl @@ -180,13 +180,13 @@ function get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = mm = create_marginal_model(m, pulse_size) - add_comp!(mm.marginal, ExtraEmissions, :extra_emissions; after=:co2emissions) - connect_param!(mm.marginal, :extra_emissions => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions) - set_param!(mm.marginal, :extra_emissions, :pulse_size, pulse_size) - set_param!(mm.marginal, :extra_emissions, :pulse_year, year) + 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) - connect_param!(mm.marginal, :co2cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted) + connect_param!(mm.modified, :co2cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted) run(mm) return mm -end \ No newline at end of file +end