Sectoral Certainty Equivalent Calculations Undefined (Appear as 0.0 …

…or -0.0) (#68)

* Fix indexing keys for norm_cpc_values_ce and update README

* Remove sectoral cpc norm ce values

* Fill with undef instead of 100 which was used for testing

* Simplify the NamedTuple for norm_cpc_values_ce and clean up commenting, add skipmissing warning

* Delete test file

* Add testing

* Make clarifying notes about ag equity weighitng

* Update how-to guide file, final edits to scc.jl

* Throw error for negative cpc in consumption per capita ew setting

README.md

@@ -338,8 +338,8 @@ You can also pass `compute_scc` a vector of `NamedTuple`s to the `discount_rates
 - label - a `String` label for the discount rate
 - prtp - a `Float64` for the pure rate of time preference Ramsey parameter
 - eta - a `Float64` for the risk aversion Ramsey parameter
-- ew - a member of `[nothing, :gdp, :consumption_region, :consumption_country]` indication whether to equity weight, and if so, whether to use gdp or consumption to do so
-- ew_norm_region - a `String` dictating the normalization region for equity weighting (a country if using `:gdp` or `:consumption_country` or a FUND region if using `:consumption_region`)
+- ew - a member of `[nothing, :gdp_region, :gdp_country, :consumption_region, :consumption_country]` indication whether to equity weight, and if so, whether to (1) use gdp or consumption (2) use regional or country-level spatial aggregation
+- ew_norm_region - a `String` dictating the normalization region for equity weighting (a country if using `:gdp_country` or `:consumption_country` or a FUND region if using `:gdp_region` or `:consumption_region`)
 
 For example:
 

test/runtests.jl

@@ -12,7 +12,7 @@ ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
     @info("test_regression_deterministic.jl")
     @time include("test_regression_deterministic.jl")
 
-    if VERSION == v"1.10" # random number generator not alwasy stable between versions
+    if VERSION == v"1.10" # random number generator not always stable between versions
         @info("test_regression_mcs.jl")
         @time include("test_regression_mcs.jl")
     end 

test/save_validation_data.jl

@@ -114,19 +114,13 @@ for gas in [:CO2, :CH4, :N2O]
     outdir = joinpath(@__DIR__, "validation_data", "validation_data_$validation_label", "default_model_MCS_SCC_2020", "$gas")
     isdir(outdir) || mkpath(outdir)
     m = MimiGIVE.get_model()
-    save_scc_mcs_data(seed, outdir, n; 
-                    m = m, year = 2020, discount_rates = discount_rates, gas = gas,
-                    save_list = save_list, save_md = true, save_cpc = true, save_slr_damages = true,
-                    compute_sectoral_values = true, compute_domestic_values = true)
+    save_scc_mcs_data(seed, outdir, n; m = m, year = 2020, discount_rates = discount_rates, gas = gas, save_list = save_list)
 end
 
 # SSP245, SC-CO2 and SC-CH4 and SC-N2O in year 2020
 for gas in [:CO2, :CH4, :N2O]
     outdir = joinpath(@__DIR__, "validation_data", "validation_data_$validation_label", "SSP245_model_MCS_SCC_2020", "$gas")
     isdir(outdir) || mkpath(outdir)
     m = MimiGIVE.get_model(; socioeconomics_source = :SSP, SSP_scenario = "SSP245")
-    save_scc_mcs_data(seed, outdir, n; 
-                    m = m, year = 2020, discount_rates = discount_rates, gas = gas,
-                    save_list = save_list, save_md = true, save_cpc = true, save_slr_damages = true,
-                    compute_sectoral_values = true, compute_domestic_values = true)
+    save_scc_mcs_data(seed, outdir, n; m = m, year = 2020, discount_rates = discount_rates, gas = gas, save_list = save_list)
 end

test/test_regression_mcs.jl

@@ -6,35 +6,6 @@ include("utils.jl")
 # label of validation data to compare AGAINST
 validation_label = "current"
 
-##------------------------------------------------------------------------------
-## Validate Model Data
-##------------------------------------------------------------------------------
-
-savevars = [
-    (compname = :DamageAggregator, varname = :total_damage),
-    (compname = :DamageAggregator, varname = :total_damage_share),
-    (compname = :DamageAggregator, varname = :total_damage_domestic),
-    (compname = :DamageAggregator, varname = :cromar_mortality_damage),
-    (compname = :DamageAggregator, varname = :agriculture_damage),
-    (compname = :DamageAggregator, varname = :energy_damage),
-    (compname = :DamageAggregator, varname = :cromar_mortality_damage_domestic),
-    (compname = :DamageAggregator, varname = :agriculture_damage_domestic),
-    (compname = :DamageAggregator, varname = :energy_damage_domestic),
-    (compname = :global_netconsumption, varname = :net_consumption),
-    (compname = :global_netconsumption, varname = :net_cpc),
-    (compname = :global_netconsumption, varname = :global_gdp),
-    (compname = :global_netconsumption, varname = :global_population),
-    (compname = :temperature, varname = :T),
-    (compname = :glaciers_small_icecaps, varname = :gsic_sea_level) ,
-    (compname = :antarctic_icesheet, varname = :ais_sea_level),
-    (compname = :greenland_icesheet, varname = :greenland_sea_level),
-    (compname = :thermal_expansion, varname = :te_sea_level),
-    (compname = :landwater_storage, varname = :lws_sea_level)
-]
-
-##------------------------------------------------------------------------------
-## Validate SCC MCS Data
-##------------------------------------------------------------------------------
 discount_rates = [
                     # Constant discount rates
                     (label = "CR 1%", prtp = 0.01, eta = 0.0), (label = "CR 2%", prtp = 0.02, eta = 0.0), (label = "CR 2.5%", prtp = 0.025, eta = 0.0), (label = "CR 3%", prtp = 0.03, eta = 0.0), (label = "CR 5%", prtp = 0.05, eta = 0.0),

test/utils.jl

@@ -91,7 +91,23 @@ function save_scc_data(outdir::String;
 end
 
 """
-    save_scc_mcs_data()
+    save_scc_mcs_data(seed::Int, outdir::String, n::Int;
+                            m::Model=MimiGIVE.get_model(), 
+                            year::Union{Int, Nothing} = nothing, 
+                            last_year::Int = MimiGIVE._model_years[end], 
+                            discount_rates=nothing,
+                            certainty_equivalent=true,
+                            gas::Symbol = :CO2,
+                            save_list::Vector = [],
+                            save_md::Bool = true,
+                            save_cpc::Bool = true,
+                            save_slr_damages::Bool = true,
+                            compute_sectoral_values::Bool = true,
+                            compute_domestic_values::Bool = true,
+                            CIAM_foresight::Symbol = :perfect,
+                            CIAM_GDPcap::Bool = false,
+                            pulse_size::Float64 = 1.,
+                        )
 
 Save the data from the user-specifified  SCC Monte Carlo Simulation computation 
 using model `m` into output folder `outdir`.
@@ -101,17 +117,17 @@ function save_scc_mcs_data(seed::Int, outdir::String, n::Int;
                             year::Union{Int, Nothing} = nothing, 
                             last_year::Int = MimiGIVE._model_years[end], 
                             discount_rates=nothing,
-                            certainty_equivalent=false,
+                            certainty_equivalent=true,
                             gas::Symbol = :CO2,
                             save_list::Vector = [],
-                            save_md::Bool = false,
-                            save_cpc::Bool = false,
-                            save_slr_damages::Bool = false,
-                            compute_sectoral_values::Bool = false,
-                            compute_domestic_values::Bool = false,
+                            save_md::Bool = true,
+                            save_cpc::Bool = true,
+                            save_slr_damages::Bool = true,
+                            compute_sectoral_values::Bool = true,
+                            compute_domestic_values::Bool = true,
                             CIAM_foresight::Symbol = :perfect,
                             CIAM_GDPcap::Bool = false,
-                            pulse_size::Float64 = 1.
+                            pulse_size::Float64 = 1.,
                         )
 
     Random.seed!(seed)
@@ -148,16 +164,22 @@ function save_scc_mcs_data(seed::Int, outdir::String, n::Int;
     df_expected_scc =       DataFrame(:region => [], :sector => [], :dr_label => [], :prtp => [], :eta => [], :scc => [])
     df_se_expected_scc =    DataFrame(:region => [], :sector => [], :dr_label => [], :prtp => [], :eta => [], :se => [])
     df_sccs =               DataFrame(:region => [], :sector => [], :dr_label => [], :prtp => [], :eta => [], :scc => [],  :trial => [])
+    df_ce_scc =             DataFrame(:region => [], :sector => [], :dr_label => [], :prtp => [], :eta => [], :scc => [])
+    df_ce_sccs =            DataFrame(:region => [], :sector => [], :dr_label => [], :prtp => [], :eta => [], :scc => [],  :trial => [])
 
     for (k,v) in results[:scc]
         append!(df_expected_scc, DataFrame(:region => k.region, :sector => k.sector, :dr_label => k.dr_label, :prtp => k.prtp, :eta => k.eta, :scc => v.expected_scc))
         append!(df_se_expected_scc, DataFrame(:region => k.region, :sector => k.sector, :dr_label => k.dr_label, :prtp => k.prtp, :eta => k.eta, :se => v.se_expected_scc))
         append!(df_sccs, DataFrame(:region => k.region, :sector => k.sector, :dr_label => k.dr_label, :prtp => k.prtp, :eta => k.eta, :scc => v.sccs, :trial => collect(1:length(v.sccs))))
+        append!(df_ce_scc, DataFrame(:region => k.region, :sector => k.sector, :dr_label => k.dr_label, :prtp => k.prtp, :eta => k.eta, :scc => v.ce_scc))
+        append!(df_ce_sccs, DataFrame(:region => k.region, :sector => k.sector, :dr_label => k.dr_label, :prtp => k.prtp, :eta => k.eta, :scc => v.ce_sccs, :trial => collect(1:length(v.ce_sccs))))
     end
 
     df_expected_scc|> save(joinpath(scc_outdir, "expected_scc.csv"))
     df_se_expected_scc|> save(joinpath(scc_outdir, "se_expected_scc.csv"))
     df_sccs|> save(joinpath(scc_outdir, "sccs.csv"))
+    df_ce_scc|> save(joinpath(scc_outdir, "ce_scc.csv"))
+    df_ce_sccs|> save(joinpath(scc_outdir, "ce_sccs.csv"))
 
     # marginal damages
     mds_outdir = joinpath(outdir, "mds")
@@ -307,14 +329,14 @@ end
                             year::Union{Int, Nothing} = nothing, 
                             last_year::Int = MimiGIVE._model_years[end], 
                             discount_rates=nothing,
-                            certainty_equivalent=false,
+                            certainty_equivalent::Bool=true,
                             gas::Symbol = :CO2,
                             save_list::Vector = [],
-                            save_md::Bool = false,
-                            save_cpc::Bool = false,
-                            save_slr_damages::Bool = false,
-                            compute_sectoral_values::Bool = false,
-                            compute_domestic_values::Bool = false,
+                            save_md::Bool = true,
+                            save_cpc::Bool = true,
+                            save_slr_damages::Bool = true,
+                            compute_sectoral_values::Bool = true,
+                            compute_domestic_values::Bool = true,
                             CIAM_foresight::Symbol = :perfect,
                             CIAM_GDPcap::Bool = false,
                             pulse_size::Float64 = 1.
@@ -329,14 +351,14 @@ function validate_scc_mcs_data(seed::Int, validationdir::String, n::Int;
                                 year::Union{Int, Nothing} = nothing, 
                                 last_year::Int = MimiGIVE._model_years[end], 
                                 discount_rates=nothing,
-                                certainty_equivalent=false,
+                                certainty_equivalent::Bool=true,
                                 gas::Symbol = :CO2,
                                 save_list::Vector = [],
-                                save_md::Bool = false,
-                                save_cpc::Bool = false,
-                                save_slr_damages::Bool = false,
-                                compute_sectoral_values::Bool = false,
-                                compute_domestic_values::Bool = false,
+                                save_md::Bool = true,
+                                save_cpc::Bool = true,
+                                save_slr_damages::Bool = true,
+                                compute_sectoral_values::Bool = true,
+                                compute_domestic_values::Bool = true,
                                 CIAM_foresight::Symbol = :perfect,
                                 CIAM_GDPcap::Bool = false,
                                 pulse_size::Float64 = 1.
@@ -373,7 +395,7 @@ function validate_scc_mcs_data(seed::Int, validationdir::String, n::Int;
                         pulse_size = pulse_size
                     )
 
-    # save list - just compare model_1 for now, model1 is sufficiently tested
+    # save list - just compare model_1 for now, model_2 is sufficiently tested
     # by testing the scc values
     for el in save_list
         validation_df = load(joinpath(validationdir, "results", "model_1", "$(el[1])_$(el[2]).csv")) |> DataFrame
@@ -388,6 +410,8 @@ function validate_scc_mcs_data(seed::Int, validationdir::String, n::Int;
     validation_df_expected_scc = load(joinpath(validationdir, "scc", "expected_scc.csv")) |> DataFrame
     validation_df_se_expected_scc = load(joinpath(validationdir, "scc", "se_expected_scc.csv")) |> DataFrame
     validation_df_sccs = load(joinpath(validationdir, "scc", "sccs.csv")) |> DataFrame
+    validation_df_ce_scc = load(joinpath(validationdir, "scc", "ce_scc.csv")) |> DataFrame
+    validation_df_ce_sccs = load(joinpath(validationdir, "scc", "ce_sccs.csv")) |> DataFrame
 
     for (k,v) in results[:scc]
         validation_vals = validation_df_expected_scc |> 
@@ -407,6 +431,18 @@ function validate_scc_mcs_data(seed::Int, validationdir::String, n::Int;
             DataFrame
         validation_vals = validation_vals.scc
         @test validation_vals ≈ v.sccs atol = atol
+
+        validation_vals = validation_df_ce_scc |> 
+            @filter(_.dr_label == k.dr_label && _.region == String.(k.region) && _.sector == String.(k.sector)) |> 
+            DataFrame
+        validation_vals = (validation_vals.scc)[1]
+        @test validation_vals ≈ v.ce_scc atol = atol
+
+        validation_vals = validation_df_ce_sccs |> 
+            @filter(_.dr_label == k.dr_label && _.region == String.(k.region) && _.sector == String.(k.sector)) |> 
+            DataFrame
+        validation_vals = validation_vals.scc
+        @test validation_vals ≈ v.ce_sccs atol = atol
     end
 
     # marginal damages

test/validation_data/validation_data_current/SSP245_model_MCS_SCC_2020/CH4/scc/ce_scc.csv

@@ -0,0 +1,71 @@
+"region","sector","dr_label","prtp","eta","scc"
+globe,agriculture,"CR 5%",0.05,0.0,312.9824182922838
+domestic,energy,"OtherRamsey",0.01,1.0,0.0
+domestic,total,"CR 2.5%",0.025,0.0,113.36780655189646
+domestic,total,"OtherRamsey",0.01,1.0,0.0
+globe,cromar_mortality,"CR 3%",0.03,0.0,552.4607931259561
+globe,energy,"CR 5%",0.05,0.0,34.578654996919035
+domestic,slr,"OtherRamsey",0.01,1.0,0.0
+domestic,energy,"CR 1%",0.01,0.0,6.936701164082456
+domestic,total,"DICE2016",0.015,1.45,0.0
+domestic,slr,"CR 1%",0.01,0.0,4.846629009972227
+globe,total,"OtherRamsey",0.01,1.0,1007.4106039835181
+globe,slr,"CR 2.5%",0.025,0.0,13.447182080986886
+globe,energy,"CR 3%",0.03,0.0,46.633925729644155
+domestic,cromar_mortality,"CR 2.5%",0.025,0.0,97.88085709423693
+domestic,total,"CR 5%",0.05,0.0,82.3800879929708
+domestic,agriculture,"CR 2%",0.02,0.0,8.505259389060251
+domestic,agriculture,"CR 2.5%",0.025,0.0,7.92871908030003
+domestic,cromar_mortality,"CR 3%",0.03,0.0,91.04184198818115
+domestic,slr,"CR 2.5%",0.025,0.0,2.324768054789354
+globe,total,"CR 3%",0.03,0.0,1039.6858588372527
+globe,energy,"CR 2.5%",0.025,0.0,51.86758550411205
+domestic,cromar_mortality,"CR 5%",0.05,0.0,71.42405700333973
+globe,cromar_mortality,"CR 1%",0.01,0.0,984.2774220830879
+globe,total,"DICE2016",0.015,1.45,779.1792391639924
+globe,cromar_mortality,"CR 5%",0.05,0.0,411.9309034025832
+globe,cromar_mortality,"OtherRamsey",0.01,1.0,533.5845962652041
+domestic,cromar_mortality,"DICE2016",0.015,1.45,0.0
+globe,cromar_mortality,"CR 2.5%",0.025,0.0,607.8336172067003
+domestic,agriculture,"OtherRamsey",0.01,1.0,0.0
+domestic,total,"CR 2%",0.02,0.0,123.1329395328541
+globe,cromar_mortality,"DICE2016",0.015,1.45,416.0266924328222
+globe,energy,"CR 2%",0.02,0.0,59.549300758886396
+globe,agriculture,"CR 1%",0.01,0.0,767.3609970667426
+globe,agriculture,"CR 2%",0.02,0.0,535.8724972801336
+globe,total,"CR 1%",0.01,0.0,1881.5806217327954
+globe,total,"CR 2%",0.02,0.0,1294.3247270837326
+domestic,energy,"DICE2016",0.015,1.45,0.0
+domestic,slr,"DICE2016",0.015,1.45,0.0
+domestic,agriculture,"CR 3%",0.03,0.0,7.422840854134208
+domestic,energy,"CR 2.5%",0.025,0.0,5.233462322487242
+domestic,slr,"CR 2%",0.02,0.0,2.821624061610267
+domestic,agriculture,"CR 5%",0.05,0.0,5.87861230324269
+globe,agriculture,"OtherRamsey",0.01,1.0,417.1875669347992
+domestic,agriculture,"DICE2016",0.015,1.45,0.0
+globe,agriculture,"CR 3%",0.03,0.0,429.3722190474898
+domestic,energy,"CR 5%",0.05,0.0,3.861711393985601
+domestic,cromar_mortality,"CR 1%",0.01,0.0,131.47890035135606
+domestic,cromar_mortality,"CR 2%",0.02,0.0,106.14798247608756
+domestic,total,"CR 3%",0.03,0.0,105.3187550606997
+globe,energy,"CR 1%",0.01,0.0,97.8711750533289
+domestic,slr,"CR 3%",0.03,0.0,1.9736020817986102
+globe,slr,"CR 1%",0.01,0.0,32.07102752931768
+domestic,cromar_mortality,"OtherRamsey",0.01,1.0,0.0
+globe,total,"CR 5%",0.05,0.0,766.216509611285
+domestic,energy,"CR 3%",0.03,0.0,4.880470136506528
+globe,slr,"CR 5%",0.05,0.0,6.7245329191658945
+globe,slr,"OtherRamsey",0.01,1.0,11.232975585182519
+globe,total,"CR 2.5%",0.025,0.0,1148.3770065412075
+domestic,slr,"CR 5%",0.05,0.0,1.2157072923527905
+domestic,agriculture,"CR 1%",0.01,0.0,9.97959457754908
+domestic,energy,"CR 2%",0.02,0.0,5.658073606013818
+globe,energy,"OtherRamsey",0.01,1.0,45.40546519774148
+globe,agriculture,"CR 2.5%",0.025,0.0,475.2286217486864
+globe,slr,"DICE2016",0.015,1.45,7.298594836236386
+globe,slr,"CR 2%",0.02,0.0,16.772923806389514
+domestic,total,"CR 1%",0.01,0.0,153.241825103034
+globe,slr,"CR 3%",0.03,0.0,11.218920933538953
+globe,agriculture,"DICE2016",0.015,1.45,320.92596642914845
+globe,energy,"DICE2016",0.015,1.45,34.927985465391934
+globe,cromar_mortality,"CR 2%",0.02,0.0,682.1300052375107