Repository: IndEcol/pymrio

Issue:
Ghosh implementation for downstream scope 3
The current pymrio has all relevant matrices related the Leontief attribution model implemented. This can be used to calculate the factors of production that are used as input to a sector or product, covering the upstream indirect inputs, also referred to as upstream scope 3. The scope 3 of a sector also includes a downstream component, which consists of the factors of production associated with the sector's input to other sectors. An attribution methodology that covers this so-called downstream scope 3 is the Ghosh model. Mathematically, there are many similarities between the Ghosh and Leontief models, and it thus makes sense to implement Ghosh in analogy with Leontief.

The proposal is to 
* implement the equivalent of A for Ghosh (often referred to A* in literature)
* implement the Ghosh inverse (often referred to G in literature)
* implement the downstream scope 3 multiplier, M_{down}, such the sum of the M+M_{down} is the full scope multiplier, with M the existing multiplier in pymrio that covers scope 1,2&3 upstream.
* a short addition to the pymrio background page that introduces the Ghosh model
* tests that test the functionality of the added functions

Hints from issue comments:
@konstantinstadler: That would be a fantastic contribution. I would love to add this to pymrio and happy to receive a pull request for this

Solution patch:
diff --git a/doc/source/api_references.rst b/doc/source/api_references.rst
index 3acf7e0..67bf8f3 100644
--- 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
index 73e9d3f..e9532e5 100644
--- 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
index 2b56957..c23a4b4 100644
--- 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
index 60844e9..8c397fe 100644
--- a/pymrio/core/mriosystem.py
+++ b/pymrio/core/mriosystem.py
@@ -32,11 +32,14 @@ from pymrio.core.constants import (
 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 @@ class Extension(BaseSystem):
         S=None,
         S_Y=None,
         M=None,
+        M_down=None,
         D_cba=None,
         D_pba=None,
         D_imp=None,
@@ -790,6 +796,7 @@ class Extension(BaseSystem):
         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 @@ class Extension(BaseSystem):
     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 @@ class Extension(BaseSystem):
         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 @@ class Extension(BaseSystem):
             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 @@ class Extension(BaseSystem):
                         "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 @@ class IOSystem(BaseSystem):
         Z=None,
         Y=None,
         A=None,
+        As=None,
         x=None,
         L=None,
+        G=None,
         unit=None,
         population=None,
         system=None,
@@ -1717,7 +1740,9 @@ class IOSystem(BaseSystem):
         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 @@ class IOSystem(BaseSystem):
         """
         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 @@ class IOSystem(BaseSystem):
             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 @@ class IOSystem(BaseSystem):
             )
             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
 
diff --git a/pymrio/tools/iomath.py b/pymrio/tools/iomath.py
index a92f635..323b0fd 100644
--- a/pymrio/tools/iomath.py
+++ b/pymrio/tools/iomath.py
@@ -4,7 +4,7 @@ All methods here should follow the functional programming paradigm
 
 Note
 ----
-To avoid namespace pollution everythin here starts with calc_
+To avoid namespace pollution everything here starts with calc_
 
 """
 
@@ -152,6 +152,50 @@ def calc_A(Z, x):
         return Z * recix
 
 
+def calc_As(Z, x):
+    """Calculate the As matrix (coefficients) from Z and x
+
+    As is a normalized version of the industrial flows of the transpose of Z, which quantifies the input to downstream
+    sectors.
+
+    Parameters
+    ----------
+    Z : pandas.DataFrame or numpy.array
+        Symmetric input output table (flows)
+    x : pandas.DataFrame or numpy.array
+        Industry output column vector
+
+    Returns
+    -------
+    pandas.DataFrame or numpy.array
+        Symmetric input output table (coefficients) As
+        The type is determined by the type of Z.
+        If DataFrame index/columns as Z
+
+    """
+    if (type(x) is pd.DataFrame) or (type(x) is pd.Series):
+        x = x.values
+    if (type(x) is not np.ndarray) and (x == 0):
+        recix = 0
+    else:
+        with warnings.catch_warnings():
+            # catch the divide by zero warning
+            # we deal wit that by setting to 0 afterwards
+            warnings.simplefilter("ignore")
+            recix = 1 / x
+        recix[recix == np.inf] = 0
+        recix = recix.reshape((1, -1))
+    # use numpy broadcasting - factor ten faster
+    # Mathematical form - slow
+    # return Z.dot(np.diagflat(recix))
+    if type(Z) is pd.DataFrame:
+        return pd.DataFrame(
+            np.transpose(Z.values) * recix, index=Z.index, columns=Z.columns
+        )
+    else:
+        return np.transpose(Z) * recix
+
+
 def calc_L(A):
     """Calculate the Leontief L from A
 
@@ -188,6 +232,65 @@ def calc_L(A):
         return np.linalg.inv(I - A)
 
 
+def calc_G(As, L=None, x=None):
+    """Calculate the Ghosh inverse matrix G either from As (high computation effor) or from Leontief matrix L and x
+    (low computation effort)
+
+    G = inverse matrix of (I - As) = hat(x)^{-1} *  L^T * hat(x)
+
+    where I is an identity matrix of same shape as As, and hat(x) is the diagonal matrix with values of x on the
+    diagonal.
+
+    Note that we define G as the transpose of the Ghosh inverse matrix, so that we can apply the factors of
+    production intensities from the left-hand-side for both Leontief and Ghosh attribution. In this way the
+    multipliers have the same (vector) dimensions and can be added.
+
+    Parameters
+    ----------
+    As : pandas.DataFrame or numpy.array
+        Symmetric input output table (coefficients)
+
+    Returns
+    -------
+    pandas.DataFrame or numpy.array
+        Ghosh input output table G
+        The type is determined by the type of As.
+        If DataFrame index/columns as As
+
+    """
+    # if L has already been calculated, then G can be derived from it with low computational cost.
+    if L is not None and x is not None:
+        if (type(x) is pd.DataFrame) or (type(x) is pd.Series):
+            x = x.values
+        if (type(x) is not np.ndarray) and (x == 0):
+            recix = 0
+        else:
+            with warnings.catch_warnings():
+                # catch the divide by zero warning
+                # we deal wit that by setting to 0 afterwards
+                warnings.simplefilter("ignore")
+                recix = 1 / x
+            recix[recix == np.inf] = 0
+            recix = recix.reshape((1, -1))
+
+        if type(L) is pd.DataFrame:
+            return pd.DataFrame(
+                recix * np.transpose(L.values * x), index=Z.index, columns=Z.columns
+            )
+        else:
+            # G = hat(x)^{-1} *  L^T * hat(x) in mathematical form np.linalg.inv(hatx).dot(L.transpose()).dot(hatx).
+            # it is computationally much faster to multiply element-wise because hatx is a diagonal matrix.
+            return recix * np.transpose(L * x)
+    else:  # calculation of the inverse of I-As has a high computational cost.
+        I = np.eye(As.shape[0])  # noqa
+        if type(As) is pd.DataFrame:
+            return pd.DataFrame(
+                np.linalg.inv(I - As), index=As.index, columns=As.columns
+            )
+        else:
+            return np.linalg.inv(I - As)  # G = inverse matrix of (I - As)
+
+
 def calc_S(F, x):
     """Calculate extensions/factor inputs coefficients
 
@@ -308,13 +411,38 @@ def calc_M(S, L):
     -------
     pandas.DataFrame or numpy.array
         Multipliers M
-        The type is determined by the type of D.
-        If DataFrame index/columns as D
+        The type is determined by the type of S.
+        If DataFrame index/columns as S
 
     """
     return S.dot(L)
 
 
+def calc_M_down(S, G):
+    """Calculate downstream multipliers of the extensions
+
+    M_down = S * ( G - I )
+
+    Where I is an identity matrix of same shape as G
+
+    Parameters
+    ----------
+    G : pandas.DataFrame or numpy.array
+        Ghosh input output table G
+    S : pandas.DataFrame or numpy.array
+        Direct impact coefficients
+
+    Returns
+    -------
+    pandas.DataFrame or numpy.array
+        Downstream multipliers M
+        The type is determined by the type of S.
+        If DataFrame index/columns as S
+
+    """
+    return S.dot(G - np.eye(G.shape[0]))
+
+
 def calc_e(M, Y):
     """Calculate total impacts (footprints of consumption Y)
 
@@ -333,7 +461,9 @@ def calc_e(M, Y):
         Multipliers m
         The type is determined by the type of M.
         If DataFrame index/columns as M
-    The calcubased on multipliers M and finald demand Y"""
+        The calculation is based on multipliers M and final demand Y
+
+    """
 
     return M.dot(Y)
 
@@ -357,8 +487,6 @@ def recalc_M(S, D_cba, Y):
         Multipliers M
         The type is determined by the type of D_cba.
         If DataFrame index/columns as D_cba
-
-
     """
 
     Y_diag = ioutil.diagonalize_columns_to_sectors(Y)
@@ -472,7 +600,6 @@ def calc_gross_trade(
             - totals: df with gross total imports and exports per sector
               and region
 
-
     """
 
     Z_trade_blocks = ioutil.set_dom_block(Z, value=0)
diff --git a/pymrio/tools/ioparser.py b/pymrio/tools/ioparser.py
index 2957025..e0443d8 100644
--- a/pymrio/tools/ioparser.py
+++ b/pymrio/tools/ioparser.py
@@ -1871,16 +1871,18 @@ def parse_eora26(path, year=None, price="bp", country_names="eora"):
                     )
 
         eora_data = {
-            key: pd.read_csv(
-                zip_file.open(filename),
-                sep=eora_sep,
-                header=None,
-            )
-            if filename in zip_file.namelist()
-            else pd.read_csv(
-                indices_file.open(filename),
-                sep=eora_sep,
-                header=None,
+            key: (
+                pd.read_csv(
+                    zip_file.open(filename),
+                    sep=eora_sep,
+                    header=None,
+                )
+                if filename in zip_file.namelist()
+                else pd.read_csv(
+                    indices_file.open(filename),
+                    sep=eora_sep,
+                    header=None,
+                )
             )
             for key, filename in eora_files.items()
         }
diff --git a/pymrio/tools/ioutil.py b/pymrio/tools/ioutil.py
index 9624cf9..6f300e7 100644
--- a/pymrio/tools/ioutil.py
+++ b/pymrio/tools/ioutil.py
@@ -3,6 +3,7 @@ Utility function for pymrio
 
 KST 20140502
 """
+
 import json
 import logging
 import os


Test patch:
diff --git a/tests/test_math.py b/tests/test_math.py
index e2ecb28..3b54cdb 100644
--- a/tests/test_math.py
+++ b/tests/test_math.py
@@ -16,12 +16,15 @@ sys.path.append(os.path.join(TESTPATH, ".."))
 # the function which should be tested here
 from pymrio.tools.iomath import calc_A  # noqa
 from pymrio.tools.iomath import calc_accounts  # noqa
+from pymrio.tools.iomath import calc_As  # noqa
 from pymrio.tools.iomath import calc_e  # noqa
 from pymrio.tools.iomath import calc_F  # noqa
 from pymrio.tools.iomath import calc_F_Y  # noqa
+from pymrio.tools.iomath import calc_G  # noqa
 from pymrio.tools.iomath import calc_gross_trade  # noqa
 from pymrio.tools.iomath import calc_L  # noqa
 from pymrio.tools.iomath import calc_M  # noqa
+from pymrio.tools.iomath import calc_M_down  # noqa
 from pymrio.tools.iomath import calc_S  # noqa
 from pymrio.tools.iomath import calc_S_Y  # noqa
 from pymrio.tools.iomath import calc_x  # noqa
@@ -232,6 +235,26 @@ def td_small_MRIO():
             columns=_Z_multiindex,
         )
 
+        As = pd.DataFrame(
+            data=[
+                [0.19607843, 0.0, 0.17241379, 0.1, 0.0, 0.12820513],  # noqa
+                [0.09803922, 0.13333333, 0.05172414, 0.0, 0.19230769, 0.0],  # noqa
+                [0.01960784, 0.0, 0.34482759, 0.0, 0.01923077, 0.0],  # noqa
+                [0.11764706, 0.0, 0.06896552, 0.02, 0.0, 0.02564103],  # noqa
+                [
+                    0.09803922,
+                    0.33333333,
+                    0.03448276,
+                    0.2,
+                    0.38461538,
+                    0.25641026,
+                ],  # noqa
+                [0.1372549, 0.2, 0.0, 0.18, 0.01923077, 0.25641026],  # noqa
+            ],
+            index=_Z_multiindex,
+            columns=_Z_multiindex,
+        )
+
         L = pd.DataFrame(
             data=[
                 [
@@ -287,6 +310,61 @@ def td_small_MRIO():
             columns=_Z_multiindex,
         )
 
+        G = pd.DataFrame(
+            data=[
+                [
+                    1.33871463,
+                    0.07482651,
+                    0.38065717,
+                    0.19175489,
+                    0.04316348,
+                    0.25230906,
+                ],  # noqa
+                [
+                    0.28499301,
+                    1.37164729,
+                    0.22311379,
+                    0.15732254,
+                    0.44208094,
+                    0.20700334,
+                ],  # noqa
+                [
+                    0.05727924,
+                    0.02969614,
+                    1.54929428,
+                    0.02349465,
+                    0.05866155,
+                    0.03091401,
+                ],  # noqa
+                [
+                    0.1747153,
+                    0.02185805,
+                    0.15941084,
+                    1.05420074,
+                    0.01404088,
+                    0.07131676,
+                ],  # noqa
+                [
+                    0.58648041,
+                    0.93542302,
+                    0.39473223,
+                    0.60492361,
+                    1.95452858,
+                    0.79595212,
+                ],  # noqa
+                [
+                    0.38121958,
+                    0.41222074,
+                    0.17907022,
+                    0.34854312,
+                    0.18081884,
+                    1.48492515,
+                ],  # noqa
+            ],
+            index=_Z_multiindex,
+            columns=_Z_multiindex,
+        )
+
         x = pd.DataFrame(
             data=[
                 [51],
@@ -346,6 +424,29 @@ def td_small_MRIO():
             columns=_Z_multiindex,
         )
 
+        M_down = pd.DataFrame(
+            data=[
+                [
+                    0.48172779,
+                    0.48999986,
+                    0.74944707,
+                    0.38438353,
+                    0.48030952,
+                    0.50914162,
+                ],  # noqa
+                [
+                    0.26832875,
+                    0.25724662,
+                    0.28759602,
+                    0.18651650,
+                    0.21856841,
+                    0.25216408,
+                ],  # noqa
+            ],
+            index=["ext_type_1", "ext_type_2"],
+            columns=_Z_multiindex,
+        )
+
         D_cba = pd.DataFrame(
             data=[
                 [
@@ -536,6 +637,17 @@ def test_calc_A_MRIO(td_small_MRIO):
     pdt.assert_frame_equal(td_small_MRIO.A, calc_A(td_small_MRIO.Z, x_Tvalues))
 
 
+def test_calc_As_MRIO(td_small_MRIO):
+    pdt.assert_frame_equal(td_small_MRIO.As, calc_As(td_small_MRIO.Z, td_small_MRIO.x))
+    # we also test the different methods to provide x:
+    x_values = td_small_MRIO.x.values
+    x_Tvalues = td_small_MRIO.x.T.values
+    x_series = pd.Series(td_small_MRIO.x.iloc[:, 0])
+    pdt.assert_frame_equal(td_small_MRIO.As, calc_As(td_small_MRIO.Z, x_series))
+    pdt.assert_frame_equal(td_small_MRIO.As, calc_As(td_small_MRIO.Z, x_values))
+    pdt.assert_frame_equal(td_small_MRIO.As, calc_As(td_small_MRIO.Z, x_Tvalues))
+
+
 def test_calc_Z_MRIO(td_small_MRIO):
     pdt.assert_frame_equal(td_small_MRIO.Z, calc_Z(td_small_MRIO.A, td_small_MRIO.x))
     # we also test the different methods to provide x:
@@ -551,6 +663,10 @@ def test_calc_L_MRIO(td_small_MRIO):
     pdt.assert_frame_equal(td_small_MRIO.L, calc_L(td_small_MRIO.A))
 
 
+def test_calc_G_MRIO(td_small_MRIO):
+    pdt.assert_frame_equal(td_small_MRIO.G, calc_G(td_small_MRIO.As))
+
+
 def test_calc_S_MRIO(td_small_MRIO):
     pdt.assert_frame_equal(td_small_MRIO.S, calc_S(td_small_MRIO.F, td_small_MRIO.x))
 
@@ -572,6 +688,12 @@ def test_calc_M_MRIO(td_small_MRIO):
     pdt.assert_frame_equal(td_small_MRIO.M, calc_M(td_small_MRIO.S, td_small_MRIO.L))
 
 
+def test_calc_M_down_MRIO(td_small_MRIO):
+    pdt.assert_frame_equal(
+        td_small_MRIO.M_down, calc_M_down(td_small_MRIO.S, td_small_MRIO.G)
+    )
+
+
 def test_calc_gross_trade_MRIO(td_small_MRIO):
     gt = calc_gross_trade(td_small_MRIO.Z, td_small_MRIO.Y)
 
diff --git a/tests/test_outputs.py b/tests/test_outputs.py
index c77c920..9b18702 100644
--- a/tests/test_outputs.py
+++ b/tests/test_outputs.py
@@ -14,7 +14,6 @@
 
 """
 
-
 import os
 import sys