diff --git a/CHANGELOG.md b/CHANGELOG.md index 08a89dab..d1b162b8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- Add an option for producing SCC ranges that account for statistical uncertainty. ([PR #143](https://github.com/ClimateImpactLab/dscim/pull/143), [@davidrzhdu](https://github.com/davidrzhdu)) + ### Fixed - Fix concatenate_energy_damages netcdf saving functionality which was not clearing data encoding causing some coordinates to be truncated. ([PR #229](https://github.com/ClimateImpactLab/dscim/pull/229), [@JMGilbert](https://github.com/JMGilbert)) @@ -14,6 +18,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.5.0] - 2023-11-17 ### Added + - Add naive list of package dependencies to pyproject.toml.([PR #123](https://github.com/ClimateImpactLab/dscim/pull/123), [@brews](https://github.com/brews)) - CI, coverage, DOI badges on README. ([PR #134](https://github.com/ClimateImpactLab/dscim/pull/134), [@brews](https://github.com/brews)) diff --git a/src/dscim/menu/main_recipe.py b/src/dscim/menu/main_recipe.py index b7adafb2..0ab531d7 100644 --- a/src/dscim/menu/main_recipe.py +++ b/src/dscim/menu/main_recipe.py @@ -318,6 +318,8 @@ def scc(): self.logger.info("Processing SCC calculation ...") if self.fit_type == "quantreg": self.full_uncertainty_iqr + self.calculate_scc + self.stat_uncertainty_iqr else: if len(self.fair_aggregation) > 0: self.stream_discount_factors @@ -1121,7 +1123,9 @@ def discounted_damages(self, damages, discrate): def uncollapsed_sccs(self): """Calculate full distribution of SCCs without FAIR aggregation""" - md = self.global_consumption_no_pulse - self.global_consumption_pulse + md = ( + self.global_consumption_no_pulse - self.global_consumption_pulse + ) # this is for full uncertainty # convert to the marginal damages from a single pulse md = md * self.climate.conversion @@ -1135,9 +1139,17 @@ def uncollapsed_sccs(self): return sccs - # @cachedproperty - # def quantiles_sccs(self): - # return self.uncollapsed_sccs.quantile(self.scc_quantiles, dim=self.fair_dims) + @cachedproperty + @save(name="stat_uncertainty_iqr") + def stat_uncertainty_iqr(self): + """Calculate the distribution of quantile-weighted SCCs produced from + quantile regressions that have already been collapsed across other dimensions to give statistical-only uncertainty. + """ + return quantile_weight_quantilereg( + self.calculate_scc, + fair_dims=self.fair_dims, + quantiles=self.full_uncertainty_quantiles, + ) @cachedproperty @save(name="full_uncertainty_iqr") @@ -1146,7 +1158,9 @@ def full_uncertainty_iqr(self): quantile regressions. """ return quantile_weight_quantilereg( - self.uncollapsed_sccs, quantiles=self.full_uncertainty_quantiles + self.uncollapsed_sccs, + fair_dims=self.fair_dims, + quantiles=self.full_uncertainty_quantiles, ) def calculate_discount_factors(self, cons_pc): diff --git a/src/dscim/utils/utils.py b/src/dscim/utils/utils.py index 816d4173..c027a652 100644 --- a/src/dscim/utils/utils.py +++ b/src/dscim/utils/utils.py @@ -99,7 +99,7 @@ def get_weights(quantiles): return weights -def quantile_weight_quantilereg(array, quantiles=None): +def quantile_weight_quantilereg(array, fair_dims, quantiles=None): """Produce quantile weights of the quantile regression damages. Parameters @@ -113,10 +113,17 @@ def quantile_weight_quantilereg(array, quantiles=None): quantiles = [0.01, 0.05, 0.167, 0.25, 0.5, 0.75, 0.833, 0.95, 0.99] qr_quantiles = array.q.values + + to_stack = ["q"] + list(set(fair_dims).intersection(set(array.coords))) weights = xr.DataArray(get_weights(qr_quantiles), dims=["q"], coords=[qr_quantiles]) - ds_stacked = array.stack(obs=("simulation", "q")) - weights_by_obs = weights.sel(q=ds_stacked.obs.q) + if len(to_stack) > 1: + ds_stacked = array.stack(obs=to_stack) + weights_by_obs = weights.sel(q=ds_stacked.obs.q) + else: + ds_stacked = array.rename({to_stack[0]: "obs"}) + weights_by_obs = weights.sel(q=ds_stacked.obs) + dim = "obs" # these are quantiles of the statistical or full uncertainty, weighted by the quantile regression quantiles diff --git a/tests/test_fair_utils.py b/tests/test_fair_utils.py index 43b20dca..09b84ea0 100644 --- a/tests/test_fair_utils.py +++ b/tests/test_fair_utils.py @@ -1,9 +1,11 @@ import numpy as np import pandas as pd +import xarray as xr import pytest import statsmodels.formula.api as smf from dscim.utils.utils import modeler from dscim.utils.utils import get_weights +from dscim.utils.utils import quantile_weight_quantilereg def test_modeler(): @@ -57,3 +59,143 @@ def test_get_weights(): ) np.testing.assert_allclose(get_weights([1]), np.array([1])) np.testing.assert_allclose(get_weights([0]), np.array([1])) + + +def test_quantile_weight_quantilereg_dim(): + """ + input cases covered : + fair dim (simulation in this case) is included in uncollapsed_sccs + """ + ds_in = xr.Dataset( + { + "uncollapsed_sccs": ( + [ + "discount_type", + "model", + "ssp", + "rcp", + "simulation", + "gas", + "q", + "weitzman_parameter", + "fair_aggregation", + ], + np.ones((1, 2, 2, 2, 5, 1, 3, 1, 1)), + ), + }, + coords={ + "discount_type": (["discount_type"], ["Euler_Ramsey"]), + "model": (["model"], ["IIASA GDP", "OECD Env-Growth"]), + "ssp": (["ssp"], ["SSP3", "SSP4"]), + "rcp": (["rcp"], ["rcp245", "rcp585"]), + "simulation": (["simulation"], [0, 1, 2, 3, 4]), + "gas": (["gas"], ["CO2_Fossil"]), + "q": (["q"], [0.01, 0.5, 0.99]), + "weitzman_parameter": (["weitzman_parameter"], ["0.1"]), + "fair_aggregation": (["fair_aggregation"], ["uncollapsed"]), + }, + ) + + ds_out = quantile_weight_quantilereg( + ds_in.uncollapsed_sccs, ["simulation"], quantiles=[0.01, 0.5, 0.99] + ) + + ds_out_expected = xr.Dataset( + { + "uncollapsed_sccs": ( + [ + "discount_type", + "model", + "ssp", + "rcp", + "simulation", + "gas", + "q", + "weitzman_parameter", + "fair_aggregation", + ], + np.ones((1, 2, 2, 2, 5, 1, 3, 1, 1)), + ), + }, + coords={ + "discount_type": (["discount_type"], ["Euler_Ramsey"]), + "model": (["model"], ["IIASA GDP", "OECD Env-Growth"]), + "ssp": (["ssp"], ["SSP3", "SSP4"]), + "rcp": (["rcp"], ["rcp245", "rcp585"]), + "gas": (["gas"], ["CO2_Fossil"]), + "weitzman_parameter": (["weitzman_parameter"], ["0.1"]), + "fair_aggregation": (["fair_aggregation"], ["uncollapsed"]), + "quantile": (["quantile"], [0.01, 0.5, 0.99]), + }, + ) + + np.testing.assert_equal(ds_out, ds_out_expected) + + +def test_quantile_weight_quantilereg_nodim(): + """ + input cases covered : + fair dim (simulation in this case) is not included in uncollapsed_sccs + """ + ds_in = xr.Dataset( + { + "uncollapsed_sccs": ( + [ + "discount_type", + "model", + "ssp", + "rcp", + "gas", + "q", + "weitzman_parameter", + "fair_aggregation", + ], + np.ones((1, 2, 2, 2, 1, 3, 1, 1)), + ), + }, + coords={ + "discount_type": (["discount_type"], ["Euler_Ramsey"]), + "model": (["model"], ["IIASA GDP", "OECD Env-Growth"]), + "ssp": (["ssp"], ["SSP3", "SSP4"]), + "rcp": (["rcp"], ["rcp245", "rcp585"]), + "gas": (["gas"], ["CO2_Fossil"]), + "q": (["q"], [0.01, 0.5, 0.99]), + "weitzman_parameter": (["weitzman_parameter"], ["0.1"]), + "fair_aggregation": (["fair_aggregation"], ["uncollapsed"]), + }, + ) + + ds_out = quantile_weight_quantilereg( + ds_in.uncollapsed_sccs, ["simulation"], quantiles=[0.01, 0.5, 0.99] + ) + + ds_out_expected = xr.Dataset( + { + "uncollapsed_sccs": ( + [ + "discount_type", + "model", + "ssp", + "rcp", + "simulation", + "gas", + "q", + "weitzman_parameter", + "fair_aggregation", + ], + np.ones((1, 2, 2, 2, 5, 1, 3, 1, 1)), + ), + }, + coords={ + "discount_type": (["discount_type"], ["Euler_Ramsey"]), + "model": (["model"], ["IIASA GDP", "OECD Env-Growth"]), + "ssp": (["ssp"], ["SSP3", "SSP4"]), + "rcp": (["rcp"], ["rcp245", "rcp585"]), + "gas": (["gas"], ["CO2_Fossil"]), + "weitzman_parameter": (["weitzman_parameter"], ["0.1"]), + "fair_aggregation": (["fair_aggregation"], ["uncollapsed"]), + "quantile": (["quantile"], [0.01, 0.5, 0.99]), + }, + ) + + np.testing.assert_equal(ds_out, ds_out_expected)