diff --git a/src/components/EquityWeighting.jl b/src/components/EquityWeighting.jl index 1a323873..0a81265a 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/src/compute_scc.jl b/src/compute_scc.jl index d6a29401..c5e8975e 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 @@ -100,7 +118,7 @@ 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 @@ -108,7 +126,7 @@ function compute_scc( 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 \ No newline at end of file +end diff --git a/test/test_EquityWeighting.jl b/test/test_EquityWeighting.jl index 7ff6a33d..64949424 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/test_standard_api.jl b/test/test_standard_api.jl index 294f018c..e6615cd8 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 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