make quantreg run with statistical uncertainty

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))
 

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):

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

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)