Ghosh implementation for downstream scope 3 (#136)

* Added Ghosh implementation for downstream calculation to pymrio, including tests and an addition to the documentation. * Small changes to theoretical description and related comments added to python code. * Small changes to theoretical description and related comments added to python code. * Update environment.yml Exclude pytest 5.8.0 because this version has an issue. * formatting correction in iomath.py * Changes related to new version of black (24.3.0). Changes are in files that are not related to the functionality added with this PR. --------- Co-authored-by: BeckebanzeF1 <felix.beckebanze@rabobank.com>
IndEcol · Apr 2, 2024 · a9896ff · a9896ff
1 parent ca36bd7
commit a9896ff
Show file tree

Hide file tree

Showing 9 changed files with 363 additions and 26 deletions.
diff --git a/doc/source/api_references.rst b/doc/source/api_references.rst
@@ -161,10 +161,13 @@ numpy array.
    calc_x
    calc_Z
    calc_A
+   calc_As
    calc_L
+   calc_G
    calc_S
    calc_F
    calc_M
+   calc_M_down
    calc_e
    calc_accounts
 

diff --git a/doc/source/math.rst b/doc/source/math.rst
@@ -176,6 +176,54 @@ Similarly, total requirements (footprints in case of environmental requirements)
 Internally, the summation are implemented with the `group-by <https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.groupby.html>`_ functionality provided by the pandas package.
 
 
+
+Upstream and Downstream scope 3
+------------
+
+In the context of impact analyses the factors of production are often categorized into scope 1, 2 and 3, with scope 3
+sub-divided into upstream and downstream.
+For a MRIO the scope 1 is the direct impact of the industries. The factors of production scope 1 associated
+with some product or service in sector 'i' of monetary value 'r' is given by :math:`S e_i r',
+where :math:`e_i' is the 'i^{th}' unit vector.
+Scope 2 is the indirect impact through directly consumed energy (mostly electricity). The precise defintion of scope 2 in an MRIO depends
+ on the list of MRIO sectors that are classified as scope 2 energy suppliers. Scope 2 is therefore included in the
+upstream scope 3, which we refer to as upstream indirect impacts. The upstream multipliers are
+.. math::
+
+    \begin{equation}
+        M_{up} = S ( L - I ) = M - S.
+    \end{equation}
+The downstream scope 3 consists of the factors of production associated with the sectors' output
+that is input to other sectors. The downstream impact can be attributed with the Ghosh methodology
+(`Lenzen, 2010 <https://www.sciencedirect.com/science/article/abs/pii/S092180091000128X>`_ ).
+The downstream attribution according to Ghosh is done by the input share matrix
+
+.. math::
+
+    \begin{equation}
+        A^{*} = Z^{T} \hat{x}^{-1}
+    \end{equation}
+
+Note that we have defined this matrix in analogy with :math:`A`, meaning that the factors of production coefficient
+are applied from the right-hand side. The full downstream multiplier (including scope 1) is given
+by :math:`S G` where
+
+.. math::
+
+    \begin{equation}
+        G = (\mathrm{I}- Z^{T}\hat{x}^{-1})^{-1}
+    \end{equation}
+is the transpose of the Ghosh inverse matrix.
+The pure downstream multiplier (excluding scope 1) is given by
+
+.. math::
+
+    \begin{equation}
+        M_{down} = S((\mathrm{I}- Z^{T}\hat{x}^{-1})^{-1} -I) = S(G-I) = S(\hat{x}^{-1} L^{T} \hat{x} - I)
+    \end{equation}
+
+The sector's total impact multiplier is simply the sum of :math:`M_{up}`, :math:`S` and :math:`M_{down}`.
+
 Aggregation
 ------------
 
@@ -254,6 +302,3 @@ and final demand impacts by
     \begin{equation}
         F_{Y, agg} = F_Y(B_k \otimes \mathrm{I})^\mathrm{T}
     \end{equation}
-
-    
-
diff --git a/environment.yml b/environment.yml
@@ -17,7 +17,7 @@ dependencies:
   # testing and developing
   - pdbpp
   - country_converter >= 0.8.0
-  - pytest >= 5.4.3
+  - pytest >= 5.4.3, < 5.8.0
   - isort >= 5.6.0
   - pytest-black
   - nbval

diff --git a/pymrio/core/mriosystem.py b/pymrio/core/mriosystem.py
@@ -32,11 +32,14 @@
 from pymrio.tools.iomath import (
     calc_A,
     calc_accounts,
+    calc_As,
     calc_F,
     calc_F_Y,
+    calc_G,
     calc_gross_trade,
     calc_L,
     calc_M,
+    calc_M_down,
     calc_S,
     calc_S_Y,
     calc_x,
@@ -738,6 +741,8 @@ class Extension(BaseSystem):
         Direct impact (extensions) coefficients of final demand. Index as F_Y
     M : pandas.DataFrame
         Multipliers with multiindex as F
+    M_down : pandas.DataFrame
+        Downstream multipliers with multiindex as F
     D_cba : pandas.DataFrame
         Footprint of consumption,  further specification with
         _reg (per region) or _cap (per capita) possible
@@ -776,6 +781,7 @@ def __init__(
         S=None,
         S_Y=None,
         M=None,
+        M_down=None,
         D_cba=None,
         D_pba=None,
         D_imp=None,
@@ -790,6 +796,7 @@ def __init__(
         self.S = S
         self.S_Y = S_Y
         self.M = M
+        self.M_down = M_down
         self.D_cba = D_cba
         self.D_pba = D_pba
         self.D_imp = D_imp
@@ -842,7 +849,7 @@ def __init__(
     def __str__(self):
         return super().__str__("Extension {} with parameters: ").format(self.name)
 
-    def calc_system(self, x, Y, Y_agg=None, L=None, population=None):
+    def calc_system(self, x, Y, Y_agg=None, L=None, G=None, population=None):
         """Calculates the missing part of the extension plus accounts
 
         This method allows to specify an aggregated Y_agg for the
@@ -852,7 +859,7 @@ def calc_system(self, x, Y, Y_agg=None, L=None, population=None):
         Calculates:
 
         - for each sector and country:
-            S, S_Y (if F_Y available), M, D_cba, D_pba_sector, D_imp_sector,
+            S, S_Y (if F_Y available), M, M_down, D_cba, D_pba_sector, D_imp_sector,
             D_exp_sector
         - for each region:
             D_cba_reg, D_pba_reg, D_imp_reg, D_exp_reg,
@@ -878,6 +885,9 @@ def calc_system(self, x, Y, Y_agg=None, L=None, population=None):
             Leontief input output table L. If this is not given,
             the method recalculates M based on D_cba (must be present in
             the extension).
+        G : pandas.DataFrame or numpy.array, optional
+            Ghosh input output table G. If this is not given,
+            M_down is not calculated.
         population : pandas.DataFrame or np.array, optional
             Row vector with population per region
         """
@@ -927,6 +937,17 @@ def calc_system(self, x, Y, Y_agg=None, L=None, population=None):
                         "Recalculation of M not possible - cause: {}".format(ex)
                     )
 
+        if self.M_down is None:
+            if G is not None:
+                self.M_down = calc_M_down(self.S, G)
+                logging.debug("{} - M_down calculated based on G".format(self.name))
+            else:
+                logging.debug(
+                    "Calculation of M_down not possible because G is not available.".format(
+                        self.name
+                    )
+                )
+
         F_Y_agg = 0
         if self.F_Y is not None:
             # F_Y_agg = ioutil.agg_columns(
@@ -1699,8 +1720,10 @@ def __init__(
         Z=None,
         Y=None,
         A=None,
+        As=None,
         x=None,
         L=None,
+        G=None,
         unit=None,
         population=None,
         system=None,
@@ -1717,7 +1740,9 @@ def __init__(
         self.Y = Y
         self.x = x
         self.A = A
+        self.As = As
         self.L = L
+        self.G = G
         self.unit = unit
         self.population = population
 
@@ -1795,7 +1820,7 @@ def calc_all(self):
         """
         Calculates missing parts of the IOSystem and all extensions.
 
-        This method call calc_system and calc_extensions
+        This method calls calc_system and calc_extensions
 
         """
         self.calc_system()
@@ -1842,10 +1867,18 @@ def calc_system(self):
             self.A = calc_A(self.Z, self.x)
             self.meta._add_modify("Coefficient matrix A calculated")
 
+        if self.As is None:
+            self.As = calc_As(self.Z, self.x)
+            self.meta._add_modify("Coefficient matrix As calculated")
+
         if self.L is None:
             self.L = calc_L(self.A)
             self.meta._add_modify("Leontief matrix L calculated")
 
+        if self.G is None:
+            self.G = calc_G(self.As)
+            self.meta._add_modify("Ghosh matrix G calculated")
+
         return self
 
     def calc_extensions(self, extensions=None, Y_agg=None):
@@ -1877,7 +1910,12 @@ def calc_extensions(self, extensions=None, Y_agg=None):
             )
             ext = getattr(self, ext_name)
             ext.calc_system(
-                x=self.x, Y=self.Y, L=self.L, Y_agg=Y_agg, population=self.population
+                x=self.x,
+                Y=self.Y,
+                L=self.L,
+                G=self.G,
+                Y_agg=Y_agg,
+                population=self.population,
             )
         return self