Correspondence_Matrix.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building Correspondence Matrix for Footprint Calculation\n",
    "\n",
    "This Jupyter notebook demonstrates how to calculate the household Carbon Footprint, using Germany as an example. Prior to the calculation, the datasets of national household consumption and consumer expenditure survey are collected from Eurostat. The matrices of product taxes and trade margins are collected from EXIOBASE Supply-Use Tables (SUTs).\n",
    "\n",
    "This work is adapted from Hardadi, et al. (2020)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Step: Import of required libraries\n",
    "\n",
    "In this step all the libraries needed within the script are imported"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xlrd, xlsxwriter\n",
    "import numpy as np\n",
    "import scipy.io\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Step: Import Exiobase 3.4\n",
    "\n",
    "Now, Exiobase 3.4 will be imported. This includes the following:\n",
    "\n",
    "L-Matrix containing the Leontief-Inverse\n",
    "\n",
    "S-Matrix containing the emissions\n",
    "\n",
    "Y-Matrix containing the final demands\n",
    "\n",
    "FDE-Matrix containing the direct emissions caused by the final demands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Import L-Matrix (Leontief-Inverse)\n",
      "Import S-Matrix (Emissions).\n",
      "Import Y-Matrix (Final Demands)\n",
      "Import FDE-Matrix (Direct Emissions from Final Demand).\n",
      "Import the Names of Industry Sectors.\n",
      "Import the Names of Extension Codes.\n",
      "Import the Names of Regions.\n"
     ]
    }
   ],
   "source": [
    "Filestring_Matlab_in = \"C://Users//Gilang Hardadi//Documents//EXIOBASE_2//EXIOBASE3_10_ITC_pxp.mat\"\n",
    "#This variable is the location of the previously parsed EXIOBASE dataset\n",
    "\n",
    "print('Import L-Matrix (Leontief-Inverse)')\n",
    "MRIO_L = scipy.io.loadmat(Filestring_Matlab_in)['EB3_L_ITC']\n",
    "\n",
    "print('Import S-Matrix (Emissions).')\n",
    "MRIO_S = scipy.io.loadmat(Filestring_Matlab_in)['EB3_S_ITC']\n",
    "\n",
    "print('Import Y-Matrix (Final Demands)')\n",
    "MRIO_Y = scipy.io.loadmat(Filestring_Matlab_in)['EB3_Y']\n",
    "\n",
    "print('Import FDE-Matrix (Direct Emissions from Final Demand).')\n",
    "MRIO_FDE = scipy.io.loadmat(Filestring_Matlab_in)['EB3_FinalDemand_Emissions']\n",
    "\n",
    "print('Import the Names of Industry Sectors.')\n",
    "MRIO_Prod = scipy.io.loadmat(Filestring_Matlab_in)['EB3_ProductNames200']\n",
    "\n",
    "print('Import the Names of Extension Codes.')\n",
    "MRIO_Ext = scipy.io.loadmat(Filestring_Matlab_in)['EB3_Extensions']\n",
    "\n",
    "print('Import the Names of Regions.')\n",
    "MRIO_Reg = scipy.io.loadmat(Filestring_Matlab_in)['EB3_RegionList']\n",
    "\n",
    "NoofCountries  = 49\n",
    "NoofProducts   = 200\n",
    "NoofIndustries = 163\n",
    "NoofFinDemands = 7\n",
    "\n",
    "Population = 80.645605 #million people\n",
    "Country_ID = 5 #Index position for Germany in EXIOBASE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Step: Import characterisation factors\n",
    "\n",
    "In order to calculate the environmental footprints, characterisation factors are needed\n",
    "to convert the embodied emissions for each type of Greenhouse Gases (GHGs) to midpoint indicators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['CO2 - GWP100                               ', 'CH4 - GWP100                               ', 'N2O - GWP100                               ', 'SF6 - GWP100                               ', 'Global Warming Potential 100               ', 'Land use - crops                           ', 'Land use - pasture                         ', 'Land use - forest                          ', 'Land use                                   ', 'Total Emission relevant energy use         ', 'Total Energy inputs from nature            ', 'Total Energy supply                        ', 'Total Energy Use                           ', 'Domestic extraction  - primary crops       ', 'Domestic extraction  - crops residues      ', 'Domestic extraction  - fodder crops        ', 'Domestic extraction - grazing              ', 'Domestic extraction - wood                 ', 'Domestic extraction - aquatic animals      ', 'Domestic extraction - metal ores           ', 'Domestic extraction - nonmetallic minerals ', 'Domestic extraction - fossil fuels         ', 'Domestic extraction used Total             ', 'Domestic extraction unused                 ', 'Material extraction Total - used and unused', 'Water Consumption Green - Agriculture      ', 'Water Consumption Blue - Agriculture       ', 'Water Consumption Blue - Livestock         ', 'Water Consumption Blue - Manufacturing     ', 'Water Consumption Blue - Electricity       ', 'Water Consumption Blue - Domestic          ', 'Water Consumption Blue - Total             ', 'Water Withdrawal Blue - Manufacturing      ', 'Water Withdrawal Blue - Electricity        ', 'Water Withdrawal Blue - Domestic           ', 'Water Withdrawal Blue - Total              ']\n"
     ]
    }
   ],
   "source": [
    "ImpactFile  = xlrd.open_workbook('Characterization_EB36.xlsx')\n",
    "ImpactSheet = ImpactFile.sheet_by_name('Emissions')\n",
    "ImpactCategory_Names = []\n",
    "for m in range(0,36):\n",
    "    ImpactCategory_Names.append(ImpactSheet.cell_value(0,m+1))\n",
    "    \n",
    "MRIO_Char = np.zeros((36,1707))\n",
    "for m in range(0,36):\n",
    "    for n in range(0,1707):\n",
    "        MRIO_Char[m,n] = ImpactSheet.cell_value(n+1,m+1)\n",
    "\n",
    "print(ImpactCategory_Names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Step: Import initial correspondence matrix\n",
    "\n",
    "Initial correspondence matrix is a binary matrix which maps the consumption categories from one classification (in this case COICOP) to another classification (EXIOBASE). In this matrix, 1 represents that the initial category (COICOP) is being assigned to the target category (EXIOBASE), while 0 means that the initial category is not mapped to the target category. This matrix is built according to the United Nation classification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0. ... 0. 0. 0.]\n",
      " [1. 0. 0. ... 0. 0. 0.]\n",
      " [1. 0. 0. ... 0. 0. 0.]\n",
      " ...\n",
      " [0. 0. 0. ... 0. 0. 1.]\n",
      " [0. 0. 0. ... 0. 0. 1.]\n",
      " [0. 0. 0. ... 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "CorrespondenceMatrix_file = xlrd.open_workbook('Correspondence_Matrix.xlsx')\n",
    "CorrespondenceMatrix_sheet = CorrespondenceMatrix_file.sheet_by_name('Matrix')\n",
    "\n",
    "# Use this correspondence matrix as iteration point for obtaining matrix (matrices) transforming \n",
    "# expenditures in CES sectors into EXIOBASE sectors.\n",
    "CorrespondenceMatrixInitial = np.zeros((200,104))\n",
    "for m in range(0,200):\n",
    "    for n in range(0,104):\n",
    "        CorrespondenceMatrixInitial[m,n] = CorrespondenceMatrix_sheet.cell_value(m+1,n+1)       \n",
    "\n",
    "print(CorrespondenceMatrixInitial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Step: Import data on taxes on products and trade and transport margins\n",
    "\n",
    "This data is taken from EXIOBASE SUTs, which is necessary to later convert the consumption data in purchaser's price to basic price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "TaxTradeMargin_file = xlrd.open_workbook('DE_2010.xlsx')\n",
    "TradeMargins_sheet = TaxTradeMargin_file.sheet_by_name('trade_margins_init')\n",
    "TransMargins_sheet = TaxTradeMargin_file.sheet_by_name('transport_margins_init')\n",
    "ProductTaxes_sheet = TaxTradeMargin_file.sheet_by_name('product_taxes_init')\n",
    "\n",
    "TradeMargins = np.zeros((200))\n",
    "TransMargins = np.zeros((200))\n",
    "ProductTaxes = np.zeros((200))\n",
    "for m in range(0,200):\n",
    "    TradeMargins[m] = TradeMargins_sheet.cell_value(m+14,168)/1e6 #to convert the unit from Euro to MEuro\n",
    "    TransMargins[m] = TransMargins_sheet.cell_value(m+14,168)/1e6\n",
    "    ProductTaxes[m] = ProductTaxes_sheet.cell_value(m+14,168)/1e6\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Step: Import national household consumption data in COICOP classification\n",
    "\n",
    "The national household consumption data in COICOP classification is collected from Eurostat database. The data is available for 42 consumption categories. In this datasheet, category IDs are also assigned for each consumption category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Food', 'Non-alcoholic beverages', 'Alcoholic beverages', 'Tobacco', 'Narcotics', 'Clothing', 'Footwear', 'Actual rentals for housing', 'Imputed rentals for housing', 'Maintenance and repair of the dwelling', 'Water supply and miscellaneous services relating to the dwelling', 'Electricity, gas and other fuels', 'Furniture and furnishings, carpets and other floor coverings', 'Household textiles', 'Household appliances', 'Glassware, tableware and household utensils', 'Tools and equipment for house and garden', 'Goods and services for routine household maintenance', 'Medical products, appliances and equipment', 'Out-patient services', 'Hospital services', 'Purchase of vehicles', 'Operation of personal transport equipment', 'Transport services', 'Postal services', 'Telephone and telefax equipment', 'Telephone and telefax services', 'Audio-visual, photographic and information processing equipment', 'Other major durables for recreation and culture', 'Other recreational items and equipment, gardens and pets', 'Recreational and cultural services', 'Newspapers, books and stationery', 'Package holidays', 'Education', 'Catering services', 'Accommodation services', 'Personal care; other services n.e.c.', 'Personal effects n.e.c.', 'Social protection', 'Insurance', 'Financial services n.e.c.', 'Other services n.e.c.']\n",
      "[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23\n",
      " 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41]\n",
      "[124616.  14897.  19865.  24906.      0.  56081.  10768. 102183. 126878.\n",
      "  10187.  36564.  64840.  33887.   6970.  12554.   6884.   7574.  17415.\n",
      "  31579.  23540.  14614.  55066.  89512.  28396.   2657.   4456.  28793.\n",
      "  29013.   3059.  28225.  31972.  24347.  23953.  10358.  52271.  13987.\n",
      "  23882.  10150.  22035.  43115.  44457.      0.]\n"
     ]
    }
   ],
   "source": [
    "HH_expenditure_file = xlrd.open_workbook('nama_10_co3_p3.xls')\n",
    "HH_expenditure_sheet = HH_expenditure_file.sheet_by_name('Data')\n",
    "\n",
    "HHE_CategoryNames = []\n",
    "HHE_CategoryAggID = np.zeros(42, dtype=int)\n",
    "HH_Expenditure_DE = np.zeros(42)\n",
    "for m in range(0,42):\n",
    "    HHE_CategoryNames.append(HH_expenditure_sheet.cell_value(10,m+1))\n",
    "    HHE_CategoryAggID[m] = HH_expenditure_sheet.cell_value(9,m+1)\n",
    "    try:\n",
    "        HH_Expenditure_DE[m] = HH_expenditure_sheet.cell_value(15,m+1)\n",
    "    except ValueError:\n",
    "        HH_Expenditure_DE[m] = 0\n",
    "\n",
    "print(HHE_CategoryNames)\n",
    "print(HHE_CategoryAggID)\n",
    "print(HH_Expenditure_DE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Step: Import household expenditure survey data\n",
    "\n",
    "The household expenditure survey data in COICOP classification is collected from Eurostat database. This data is available at a higher resolution level than the national household consumption data, for 104 consumption categories. In this datasheet, category IDs are assigned accordingly to the expenditure categories in national household consumption data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Bread and cereals', 'Meat', 'Fish and seafood', 'Milk, cheese and eggs', 'Oils and fats', 'Fruit', 'Vegetables', 'Sugar, jam, honey, chocolate and confectionery', 'Food products n.e.c.', 'Non-alcoholic beverages', 'Alcoholic beverages', 'Tobacco', 'Narcotics', 'Clothing materials', 'Garments', 'Other articles of clothing and clothing accessories', 'Cleaning, repair and hire of clothing', 'Shoes and other footwear', 'Repair and hire of footwear', 'Actual rentals paid by tenants', 'Other actual rentals', 'Imputed rentals of owner-occupiers', 'Other imputed rentals', 'Materials for the maintenance and repair of the dwelling', 'Services for the maintenance and repair of the dwelling', 'Water supply and miscellaneous services relating to the dwelling', 'Electricity', 'Gas', 'Liquid fuels', 'Solid fuels', 'Heat energy', 'Furniture and furnishings', 'Carpets and other floor coverings', 'Repair of furniture, furnishings and floor coverings', 'Household textiles', 'Major household appliances whether electric or not', 'Small electric household appliances', 'Repair of household appliances', 'Glassware, tableware and household utensils', 'Tools and equipment for house and garden', 'Non-durable household goods', 'Domestic services and household services', 'Pharmaceutical products', 'Other medical products', 'Therapeutic appliances and equipment', 'Medical services', 'Dental services', 'Paramedical services', 'Hospital services', 'Motor cars', 'Motor cycles', 'Bicycles', 'Animal drawn vehicles', 'Spare parts and accessories for personal transport equipment', 'Fuels and lubricants for personal transport equipment', 'Maintenance and repair of personal transport equipment', 'Other services in respect of personal transport equipment', 'Passenger transport by railway', 'Passenger transport by road', 'Passenger transport by air', 'Passenger transport by sea and inland waterway', 'Combined passenger transport', 'Other purchased transport services', 'Postal services', 'Telephone and telefax equipment', 'Telephone and telefax services', 'Equipment for the reception, recording and reproduction of sound and picture', 'Photographic and cinematographic equipment and optical instruments', 'Information processing equipment', 'Recording media', 'Repair of audio-visual, photographic and information processing equipment', 'Major durables for outdoor recreation', 'Musical instruments and major durables for indoor recreation', 'Maintenance and repair of other major durables for recreation and culture', 'Games, toys and hobbies', 'Equipment for sport, camping and open-air recreation', 'Gardens, plants and flowers', 'Pets and related products', 'Veterinary and other services for pets', 'Recreational and sporting services', 'Cultural services', 'Games of chance', 'Books', 'Newspapers and periodicals', 'Miscellaneous printed matter', 'Stationery and drawing materials', 'Package holidays', 'Education', 'Restaurants, cafés and the like', 'Canteens', 'Accommodation services', 'Hairdressing salons and personal grooming establishments', 'Electrical appliances for personal care', 'Other appliances, articles and products for personal care ', 'Prostitution', 'Jewellery, clocks and watches', 'Other personal effects', 'Social protection', 'Insurance connected with the dwelling', 'Insurance connected with health', 'Insurance connected with transport', 'Other insurance', 'Financial services n.e.c.', 'Other services n.e.c.']\n",
      "[ 0  0  0  0  0  0  0  0  0  1  2  3  4  5  5  5  5  6  6  7  7  8  8  9\n",
      "  9 10 11 11 11 11 11 12 12 12 13 14 14 14 15 16 17 17 18 18 18 19 19 19\n",
      " 20 21 21 21 21 22 22 22 22 23 23 23 23 23 23 24 25 26 27 27 27 27 27 28\n",
      " 28 28 29 29 29 29 29 30 30 30 31 31 31 31 32 33 34 34 35 36 36 36 37 37\n",
      " 37 38 39 39 39 39 40 41]\n",
      "[ 567.  679.  119.  469.   85.  286.  368.  247.  140.  338.  265.  201.\n",
      "    0.   21.  911.   30.   29.  236.    8. 2621.   53. 3100.  259.  116.\n",
      "  120.  679.  677.  426.  353.   35.  149.  461.   53.   12.  103.  153.\n",
      "   56.   13.   93.  152.  106.   78.  251.   69.  225.  279.  144.   55.\n",
      "   76. 1111.   33.   51.    0.  189. 1144.  408.  663.  124.   22.  108.\n",
      "   12.  134.    8.   65.   38.  662.  137.   44.  142.   82.   11.   55.\n",
      "    0.   10.  103.   66.  192.  155.    0.  219.  269.  126.  145.  195.\n",
      "  160.    0.  712.  236.  909.  113.  297.  275.   11.  333.    0.   94.\n",
      "   56.   59.    0. 1142.  405.  476.   62.  123.]\n"
     ]
    }
   ],
   "source": [
    "CES_expenditure_file = xlrd.open_workbook('hbs_exp_t121.xls')\n",
    "CES_expenditure_sheet = CES_expenditure_file.sheet_by_name('Data')\n",
    "\n",
    "CES_CategoryNames = []\n",
    "CES_CategoryAggID = np.zeros(104, dtype=int)\n",
    "CES_DEexpenditure = np.zeros(104)\n",
    "for m in range(0,104):\n",
    "    CES_CategoryNames.append(CES_expenditure_sheet.cell_value(10,m+1))\n",
    "    CES_CategoryAggID[m] = CES_expenditure_sheet.cell_value(9,m+1)\n",
    "    CES_DEexpenditure[m] = CES_expenditure_sheet.cell_value(15,m+1)\n",
    "\n",
    "print(CES_CategoryNames)\n",
    "print(CES_CategoryAggID)\n",
    "print(CES_DEexpenditure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Step: Build tuple matrix to aggregate expenditures with identical indices\n",
    "\n",
    "In this step, a function is built to later aggregate expenditures with identical indices. This function is taken from PySUT library developed by Konstantin Stadler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define ancillary functions to construct aggregation matrices\n",
    "\n",
    "def MI_Tuple(value, Is): \n",
    "    \"\"\"\n",
    "    Define function for obtaining multiindex tuple from index value\n",
    "    value: flattened index position, Is: Number of values for each index dimension\n",
    "    Example: MI_Tuple(10, [3,4,2,6]) returns [0,0,1,4]\n",
    "    MI_Tuple(138, [100,10,5]) returns [2,7,3]\n",
    "    MI_Tuple is the inverse of Tuple_MI.\n",
    "    \"\"\"\n",
    "    IsValuesRev = []\n",
    "    CurrentValue = value\n",
    "    for m in range(0,len(Is)):\n",
    "        IsValuesRev.append(CurrentValue % Is[len(Is)-m-1])\n",
    "        CurrentValue = CurrentValue // Is[len(Is)-m-1]\n",
    "    return IsValuesRev[::-1]    \n",
    "\n",
    "def Tuple_MI(Tuple, IdxLength): \n",
    "    \"\"\"\n",
    "    Function to return the absolution position of a multiindex when the index tuple\n",
    "    and the index hierarchy and size are given.\n",
    "    Example: Tuple_MI([2,7,3],[100,10,5]) returns 138\n",
    "    Tuple_MI([0,0,1,4],[3,4,2,6]) returns 10\n",
    "    Tuple_MI is the inverse of MI_Tuple.\n",
    "    \"\"\"\n",
    "    # First, generate the index position offset values\n",
    "    IdxShift =  IdxLength[1:] +  IdxLength[:1] # Shift 1 to left\n",
    "    IdxShift[-1] = 1 # Replace lowest index by 1\n",
    "    IdxShift.reverse()\n",
    "    IdxPosOffset = np.cumproduct(IdxShift).tolist()\n",
    "    IdxPosOffset.reverse()\n",
    "    Position = np.sum([a*b for a,b in zip(Tuple,IdxPosOffset)])\n",
    "    return Position\n",
    "\n",
    "def build_Aggregation_Matrix(Position_Vector): # from PySUT\n",
    "    \"\"\"Turn a vector of target positions into a matrix that aggregates \n",
    "    or re-arranges rows of the table it is multiplied to from the left \n",
    "    (or columns, if multiplied transposed from the right)\"\"\"\n",
    "    AM_length = Position_Vector.max() + 1 # Maximum row number of new matrix (+1 to get the right length, as 0 is the smallest target position entry.)\n",
    "    AM_width  = len(Position_Vector) # Number of rows of the to-be-aggregated matrix\n",
    "    Rearrange_Matrix = np.zeros((AM_length,AM_width))\n",
    "    for m in range(0,len(Position_Vector)):\n",
    "        Rearrange_Matrix[Position_Vector[m].item(0),m] = 1\n",
    "        # place 1 in aggregation matrix at [PositionVector[m],m], so that column m is aggregated with Positionvector[m] in the aggregated matrix\n",
    "    return Rearrange_Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using build_Aggregation_Matrix function, the household expenditure survey data (in 104 categories) is aggregated accordingly to the national household consumption data (in 42 categories)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2960.  338.  265.  201.    0.  991.  244. 2674. 3359.  236.  679. 1640.\n",
      "  526.  103.  222.   93.  152.  184.  545.  478.   76. 1195. 2404.  408.\n",
      "   65.   38.  662.  416.   65.  516.  614.  500.  712.  236. 1022.  297.\n",
      "  619.  150.   59. 2023.   62.  123.]\n"
     ]
    }
   ],
   "source": [
    "CES_DE_agg = build_Aggregation_Matrix(CES_CategoryAggID).dot(CES_DEexpenditure)\n",
    "\n",
    "print(CES_DE_agg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the national household consumption data is disaggregated accordingly (into 104 categories) based on the household expenditure survey data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 23870.7         28585.9          5009.9         19744.9\n",
      "   3578.5         12040.6         15492.8         10398.7\n",
      "   5894.          14897.          19865.          24906.\n",
      "      0.           1188.39656912  51553.77497477   1697.70938446\n",
      "   1641.11907164  10414.95081967    353.04918033 100157.68249813\n",
      "   2025.31750187 117094.90919917   9783.09080083   5007.16949153\n",
      "   5179.83050847  36564.          26766.26829268  16842.58536585\n",
      "  13956.41463415   1383.7804878    5890.95121951  29699.44296578\n",
      "   3414.46958175    773.08745247   6970.           8652.08108108\n",
      "   3166.77477477    735.14414414   6884.           7574.\n",
      "  10032.55434783   7382.44565217  14543.72293578   3998.07522936\n",
      "  13037.20183486  13739.87447699   7091.54811715   2708.57740586\n",
      "  14614.          51195.25188285   1520.65104603   2350.09707113\n",
      "      0.           7037.34109817  42596.39267887  15191.72046589\n",
      "  24686.54575707   8630.15686275   1531.15686275   7516.58823529\n",
      "    835.17647059   9326.1372549     556.78431373   2657.\n",
      "   4456.          28793.           9554.76201923   3068.68269231\n",
      "   9903.47596154   5718.90865385    767.17067308   2588.38461538\n",
      "      0.            470.61538462   5634.06007752   3610.1744186\n",
      "  10502.3255814    8478.43992248      0.          11403.69381107\n",
      "  14007.27687296   6561.02931596   7060.63         9495.33\n",
      "   7791.04            0.          23953.          10358.\n",
      "  46491.52544031   5779.47455969  13987.          10609.93537964\n",
      "    424.39741519  12847.66720517      0.           6360.66666667\n",
      "   3789.33333333  22035.              0.          24338.76915472\n",
      "   8631.52496293  10144.70588235  44457.              0.        ]\n"
     ]
    }
   ],
   "source": [
    "HHE_DE_dag = np.zeros((104))\n",
    "\n",
    "for i in range(104):\n",
    "    idx = CES_CategoryAggID[i]\n",
    "    if CES_DE_agg[idx] != 0:\n",
    "        HHE_DE_dag[i] = CES_DEexpenditure[i]/CES_DE_agg[idx]*HH_Expenditure_DE[idx]\n",
    "    else:\n",
    "        HHE_DE_dag[i] = 0\n",
    "\n",
    "print(HHE_DE_dag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9. Step: Calculate household final demand data for each EXIOBASE category in basic and purchaser's price\n",
    "\n",
    "The MRIO_Y matrix contains national final demand data (household, government, capital formation, etc.) in 9800 region-products (200 products from 49 regions). In this step, the household final demand for German household (in basic price) is extracted and aggregated into 200 products."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.83012439e+00 1.22059515e+03 3.07370145e+02 1.22834965e+04\n",
      " 7.61501082e+02 0.00000000e+00 2.44661457e-05 2.15723580e+03\n",
      " 0.00000000e+00 0.00000000e+00 2.49159199e+03 4.15914350e+02\n",
      " 3.27180124e+02 5.11700480e+01 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 1.15481870e+03 3.84112240e+02 6.16371465e+01\n",
      " 0.00000000e+00 1.76790782e+02 1.02845415e-04 2.78248336e+01\n",
      " 0.00000000e+00 3.24660292e+02 1.46274558e-04 0.00000000e+00\n",
      " 1.25500914e+04 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 1.48667191e+02 2.49268368e+03 5.55179302e+03\n",
      " 7.01694183e+03 1.17963142e+04 4.99298482e+02 1.55820270e+04\n",
      " 1.74224723e+02 1.46640037e+03 5.65244341e+04 9.65656597e+03\n",
      " 7.39217581e+03 3.65131686e+03 1.25984629e+04 1.88603482e+04\n",
      " 6.51960189e+03 3.41460397e+03 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 5.79152016e+03 2.46127399e+04 8.17294376e+00\n",
      " 0.00000000e+00 9.09953930e-05 2.30674216e+04 6.08536413e-01\n",
      " 0.00000000e+00 2.26742370e+02 1.31087450e+01 9.62830880e+03\n",
      " 0.00000000e+00 4.72232831e-05 8.10317558e+02 3.03087125e-05\n",
      " 7.59222499e-05 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 8.21173585e+01\n",
      " 1.02832234e+02 2.70567295e+04 9.14357362e-05 1.33297752e+02\n",
      " 3.06662056e+02 7.56690670e+01 1.24085403e+02 8.45325210e+03\n",
      " 0.00000000e+00 0.00000000e+00 1.76240592e+03 6.61948955e+00\n",
      " 4.70581708e+01 0.00000000e+00 3.08354549e+03 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 7.29275433e+03 2.92778471e+03 5.29261677e+03 3.38383912e+03\n",
      " 1.27271405e+04 1.02952483e+04 6.29993618e+04 3.87966037e+03\n",
      " 3.36476573e+04 0.00000000e+00 0.00000000e+00 3.12016084e+03\n",
      " 1.06453151e+03 1.69573511e+03 2.44360443e+02 4.56789330e+02\n",
      " 1.01533395e+02 4.75482896e+02 1.42107938e+02 0.00000000e+00\n",
      " 0.00000000e+00 1.51314595e+00 1.56810636e+02 3.08933232e+03\n",
      " 1.57684038e+04 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 1.55951313e-04 0.00000000e+00 7.45665329e+03 6.59012612e+02\n",
      " 4.61925714e+03 4.58302765e+03 0.00000000e+00 2.65640202e+04\n",
      " 3.04633339e+02 1.62428763e+03 6.56318221e+02 8.00877780e+04\n",
      " 1.16205046e+04 1.38834022e+04 0.00000000e+00 5.51048201e+02\n",
      " 4.46749728e+02 1.91466310e+04 8.79931877e+03 3.59123285e+04\n",
      " 4.16561901e+04 3.79072238e+04 6.37888452e+02 2.33767398e+05\n",
      " 4.90587614e+03 3.21082853e+02 1.79207427e+02 1.43060479e+04\n",
      " 1.10570408e+04 1.69167815e+04 3.55278670e+04 1.96974948e+03\n",
      " 2.78926484e+03 2.74938790e+03 3.30850977e+03 4.53014507e+02\n",
      " 1.68622714e+03 2.63829891e+03 5.03823917e+02 6.22909344e+01\n",
      " 2.74338851e+03 1.48990347e+03 3.20464017e+01 7.53579578e+02\n",
      " 1.39989249e+03 2.78819802e+01 6.31461224e+01 5.23349459e+01\n",
      " 6.59900291e+02 5.23674936e+00 3.37837563e+01 6.23288496e+02\n",
      " 2.62600333e+04 3.52357852e+04 2.90008870e+03 0.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "HH_FinalDemand_bp = MRIO_Y[:,NoofFinDemands*Country_ID]\n",
    "HH_FinalDemand_bp = HH_FinalDemand_bp.reshape(49,200).sum(axis=0)\n",
    "print(HH_FinalDemand_bp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, the trade and transport margins and product taxes are added to convert the national household demand from basic price to purchaser's price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.89142148e+00 1.26706845e+03 3.21613939e+02 1.25223196e+04\n",
      " 8.23570613e+02 0.00000000e+00 2.44661457e-05 2.20003671e+03\n",
      " 0.00000000e+00 0.00000000e+00 2.62235062e+03 4.59566290e+02\n",
      " 3.29645391e+02 5.21324128e+01 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 1.17521055e+03 3.95117877e+02 6.27123984e+01\n",
      " 0.00000000e+00 1.78294282e+02 1.02845415e-04 2.95115167e+01\n",
      " 0.00000000e+00 3.39638884e+02 1.46274558e-04 0.00000000e+00\n",
      " 1.35023557e+04 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 1.52016006e+02 2.68900188e+03 5.79972027e+03\n",
      " 7.33103620e+03 1.36730117e+04 5.14010038e+02 1.67138395e+04\n",
      " 1.84103188e+02 1.54853487e+03 6.06836365e+04 1.08044907e+04\n",
      " 9.77654138e+03 1.69284068e+04 1.50129826e+04 2.60286391e+04\n",
      " 8.69610333e+03 3.57830317e+03 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 5.88276498e+03 2.61814576e+04 1.58577919e+01\n",
      " 0.00000000e+00 9.09953930e-05 3.21207661e+04 7.80519132e-01\n",
      " 0.00000000e+00 2.84386826e+02 1.96870372e+01 1.71414554e+04\n",
      " 0.00000000e+00 4.72232831e-05 1.18264675e+03 3.03087125e-05\n",
      " 7.59222499e-05 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 8.54308946e+01\n",
      " 1.11414708e+02 2.83850257e+04 9.14357362e-05 1.36026626e+02\n",
      " 3.17359986e+02 1.01059944e+02 1.51057204e+02 8.87914192e+03\n",
      " 0.00000000e+00 0.00000000e+00 1.81059427e+03 6.75927769e+00\n",
      " 4.80082696e+01 0.00000000e+00 3.15267664e+03 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 7.73787405e+03 2.99916902e+03 6.25256716e+03 3.57745933e+03\n",
      " 1.38445518e+04 1.17151409e+04 6.46408491e+04 4.12252549e+03\n",
      " 4.08347991e+04 0.00000000e+00 0.00000000e+00 3.50487463e+03\n",
      " 1.22050690e+03 1.89044664e+03 2.69146011e+02 5.10466082e+02\n",
      " 1.14695498e+02 5.29265573e+02 1.58875609e+02 0.00000000e+00\n",
      " 0.00000000e+00 1.58774905e+00 1.75174215e+02 3.44688696e+03\n",
      " 1.75640338e+04 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      " 1.55951313e-04 0.00000000e+00 8.19876616e+03 7.47625957e+02\n",
      " 5.69093970e+03 5.05540028e+03 0.00000000e+00 2.81745023e+04\n",
      " 3.07226772e+02 8.71956274e+02 4.51523140e+02 9.26474744e+04\n",
      " 1.10461492e+04 1.31472925e+04 0.00000000e+00 5.51815234e+02\n",
      " 4.47397227e+02 1.96872504e+04 8.86484462e+03 4.10800969e+04\n",
      " 4.20936404e+04 4.27135298e+04 6.44779252e+02 2.35463872e+05\n",
      " 5.03901770e+03 3.27776998e+02 1.84841957e+02 1.50904644e+04\n",
      " 1.11871540e+04 1.72921784e+04 3.56586947e+04 1.98115794e+03\n",
      " 2.80554155e+03 2.76539989e+03 3.32762231e+03 4.55613249e+02\n",
      " 1.69601140e+03 2.65369930e+03 5.06652235e+02 6.26273664e+01\n",
      " 2.75899430e+03 1.49889626e+03 3.21748645e+01 7.59812580e+02\n",
      " 1.40889304e+03 2.87434338e+01 6.39798906e+01 5.26866630e+01\n",
      " 6.64202024e+02 5.35689237e+00 3.40858764e+01 6.48030703e+02\n",
      " 2.76546501e+04 3.58558386e+04 2.97297594e+03 0.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "HH_FinalDemand_pp = HH_FinalDemand_bp + TradeMargins + TransMargins + ProductTaxes\n",
    "print(HH_FinalDemand_pp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The national household final demand data in EXIOBASE is not much different compared to the national household consumption data from Eurostat, with less than 3% difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1280882.5257171914\n",
      "1316506.0\n"
     ]
    }
   ],
   "source": [
    "print(HH_FinalDemand_pp.sum(axis=0))\n",
    "print(HHE_DE_dag.sum(axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10. Step: Calculate Tax and Margin Rates\n",
    "\n",
    "The tax rates are calculated by diving the taxes by the household final demand in purchaser's price, while the margin rates are calculated by dividing the margins by the household final demand in producer's price (purchaser's price without the taxes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.02418355  0.03126361  0.0377614   0.0145143   0.05942865  0.\n",
      "  0.          0.01714576  0.          0.          0.03937407  0.07077579\n",
      "  0.00171403  0.01534796  0.          0.          0.          0.01699467\n",
      "  0.02473512  0.01628453  0.          0.0081102   0.          0.05689677\n",
      "  0.          0.04296694  0.          0.          0.06516156  0.\n",
      "  0.          0.          0.          0.          0.          0.\n",
      "  0.          0.          0.          0.          0.          0.0194266\n",
      "  0.06995704  0.0412926   0.02760297  0.13658788  0.02815349  0.06393749\n",
      "  0.05141695  0.05092256  0.06474677  0.10432469  0.23061846  0.78211018\n",
      "  0.15145742  0.26338579  0.23981205  0.04090444  0.          0.\n",
      "  0.          0.01317137  0.05440478  0.48111202  0.          0.\n",
      "  0.28068241  0.21975195  0.          0.20178177  0.33402865  0.43671875\n",
      "  0.          0.          0.31357627  0.          0.          0.\n",
      "  0.          0.          0.          0.          0.          0.\n",
      "  0.          0.          0.          0.03688768  0.07311211  0.0440122\n",
      "  0.          0.01983822  0.03156533  0.24662225  0.17393123  0.04580287\n",
      "  0.          0.          0.02026626  0.01623217  0.01558819  0.\n",
      "  0.01649887  0.          0.          0.          0.          0.\n",
      "  0.          0.          0.          0.          0.          0.\n",
      "  0.          0.          0.05453341  0.02123474  0.14849457  0.05162499\n",
      "  0.07709392  0.11576526  0.02249599  0.05737728  0.16833161  0.\n",
      "  0.          0.10976535  0.12779558  0.10299763  0.09208967  0.10515243\n",
      "  0.11475693  0.10161756  0.10553962  0.          0.          0.0469867\n",
      "  0.10483037  0.10373263  0.10223335  0.          0.          0.\n",
      "  0.          0.          0.15584151  0.11852631  0.18831381  0.09343921\n",
      "  0.          0.08889204  0.15045515  0.          0.          0.13556437\n",
      " -0.05199598 -0.05598945  0.          0.00139002  0.00144726  0.02746038\n",
      "  0.00739165  0.12579738  0.01039231  0.11252421  0.01068707  0.00720481\n",
      "  0.02642213  0.02042286  0.03048296  0.05198094  0.01163059  0.02170906\n",
      "  0.00366889  0.00575848  0.00580163  0.00579012  0.0057436   0.00570383\n",
      "  0.00576898  0.00580337  0.00558236  0.00537196  0.00565633  0.00599961\n",
      "  0.00399265  0.00820334  0.00638838  0.02997045  0.01303172  0.00667564\n",
      "  0.00647654  0.02242775  0.0088635   0.03818061  0.05042974  0.01729295\n",
      "  0.02451659  0.        ]\n",
      "[0.00842823 0.00558892 0.00678323 0.00462462 0.01694473 0.\n",
      " 0.         0.00234916 0.         0.         0.01091899 0.02605323\n",
      " 0.00577441 0.00316055 0.         0.         0.         0.00036316\n",
      " 0.00319805 0.00087548 0.         0.00032512 0.         0.0002721\n",
      " 0.         0.00118552 0.         0.         0.00573813 0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.00265433\n",
      " 0.00328027 0.00151823 0.01567416 0.00077335 0.00048121 0.00403776\n",
      " 0.0023617  0.00223119 0.00405488 0.00214414 0.01724496 0.01008821\n",
      " 0.01104407 0.01631024 0.01377636 0.00504984 0.         0.\n",
      " 0.         0.00237038 0.0058295  0.00674175 0.         0.\n",
      " 0.0016278  0.00075881 0.         0.00114704 0.00017221 0.00281206\n",
      " 0.         0.         0.00182217 0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.00197119 0.00422886 0.00291162\n",
      " 0.         0.00022762 0.00221368 0.00613696 0.00559557 0.00226613\n",
      " 0.         0.         0.00647972 0.00452216 0.00426866 0.\n",
      " 0.00551996 0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.00316394 0.0026223  0.00591239 0.00263322\n",
      " 0.00391952 0.00614796 0.00296466 0.00162785 0.00922688 0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.         0.         0.         0.         0.\n",
      " 0.         0.        ]\n"
     ]
    }
   ],
   "source": [
    "Tax_Rate = np.zeros((200))\n",
    "Mar_Rate = np.zeros((200))\n",
    "\n",
    "for m in range(0,200):\n",
    "    if HH_FinalDemand_pp[m] == 0:\n",
    "        Tax_Rate[m] = 0\n",
    "        Mar_Rate[m] = 0\n",
    "    else:\n",
    "        Tax_Rate[m] = ProductTaxes[m] / HH_FinalDemand_pp[m]\n",
    "        HH_FinalDemand_pdp = HH_FinalDemand_pp[m] - ProductTaxes[m]\n",
    "        if TradeMargins[m] < 0:\n",
    "            Mar_Rate[m] = 0\n",
    "        else:\n",
    "            Mar_Rate[m] = (TradeMargins[m] + TransMargins[m]) / HH_FinalDemand_pdp\n",
    "            \n",
    "print(Tax_Rate)\n",
    "print(Mar_Rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 11. Step: Calculate the correspondence matrix\n",
    "\n",
    "The correspondence matrix (or concordance matrix) is a matrix that assigns the household consumption from home categories (in this case COICOP) to the target categories (EXIOBASE) according to the proportion of the consumption in the target categories. While the initial correspondence matrices for different countries should be similar, the correspondence matrix for each country will be different since the consumption structure for each country is different.\n",
    "\n",
    "The correspondence matrix is calculated using RAS method (Miller and Blair, 2009), an iterative method to balance a matrix, where the desired result is that the axial sum of the matrix is identical to the target vector (in this case the household consumption in EXIOBASE categories, purchaser's price)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RAS(i,j,k):\n",
    "    IterationMatrix0 = i\n",
    "    HH_Expenditure = j\n",
    "    Final_Demand = k\n",
    "    IterationMatrix1 = IterationMatrix0.dot(np.diag(np.nan_to_num(HH_Expenditure/(IterationMatrix0.sum(axis=0)))))\n",
    "    counter = 0\n",
    "    while True:\n",
    "        IterationMatrix2 = np.diag(np.nan_to_num(Final_Demand/(IterationMatrix1.sum(axis=1)))).dot(IterationMatrix1)\n",
    "        IterationMatrix1 = IterationMatrix2.dot(np.diag(np.nan_to_num(HH_Expenditure/(IterationMatrix2.sum(axis=0)))))\n",
    "        counter = counter + 1\n",
    "        if counter == 500:\n",
    "            return IterationMatrix1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-16-8e7f7f2fefc2>:5: RuntimeWarning: invalid value encountered in true_divide\n",
      "  IterationMatrix1 = IterationMatrix0.dot(np.diag(np.nan_to_num(HH_Expenditure/(IterationMatrix0.sum(axis=0)))))\n",
      "<ipython-input-16-8e7f7f2fefc2>:8: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  IterationMatrix2 = np.diag(np.nan_to_num(Final_Demand/(IterationMatrix1.sum(axis=1)))).dot(IterationMatrix1)\n",
      "<ipython-input-16-8e7f7f2fefc2>:8: RuntimeWarning: invalid value encountered in true_divide\n",
      "  IterationMatrix2 = np.diag(np.nan_to_num(Final_Demand/(IterationMatrix1.sum(axis=1)))).dot(IterationMatrix1)\n",
      "<ipython-input-16-8e7f7f2fefc2>:9: RuntimeWarning: invalid value encountered in true_divide\n",
      "  IterationMatrix1 = IterationMatrix2.dot(np.diag(np.nan_to_num(HH_Expenditure/(IterationMatrix2.sum(axis=0)))))\n"
     ]
    }
   ],
   "source": [
    "Correspondence_Matrix = RAS(CorrespondenceMatrixInitial, HHE_DE_dag, HH_FinalDemand_pp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the result from RAS procedure is being normalized to obtain the correspondence matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[8.19420616e-05 5.48932122e-02 1.39332822e-02 ... 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00]\n",
      " ...\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-18-78a08803f9e1>:1: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Correspondence_Matrix = np.nan_to_num(Correspondence_Matrix/Correspondence_Matrix.sum(axis=0)).transpose()\n"
     ]
    }
   ],
   "source": [
    "Correspondence_Matrix = np.nan_to_num(Correspondence_Matrix/Correspondence_Matrix.sum(axis=0)).transpose()\n",
    "print(Correspondence_Matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12. Step: Bridge the household expenditure survey data into EXIOBASE\n",
    "\n",
    "Initially, the expenditure survey data (in purchaser's price) is bridged into EXIOBASE categories (200 products) using the correspondence matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.04646115 31.12445133  7.900171   ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " ...\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]]\n"
     ]
    }
   ],
   "source": [
    "# The expenditure vector is diagonalized in order to obtain the footprint value specifically for each product\n",
    "Expenditure_pp = np.diag(CES_DEexpenditure).dot(Correspondence_Matrix)\n",
    "print(Expenditure_pp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Afterwards, the expenditures are converted into basic price using the tax and margin rates information. The margins are later assigned to the \"Wholesale trade\" and \"Retail trade\" categories, to take into account the emissions happening due to trade activities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.04495544 30.11577167  7.64414266 ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " ...\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.        ]]\n"
     ]
    }
   ],
   "source": [
    "Expenditure_bp = np.zeros((104,200))\n",
    "\n",
    "for m in range(0,104):\n",
    "    for n in range(0,200):\n",
    "        Expenditure_bp[m,n] = (1-Mar_Rate[m]) * (1-Tax_Rate[m]) * Expenditure_pp[m,n]\n",
    "    Expenditure_bp[m,153] = 0.68 * (Mar_Rate[0:124] * (1-Tax_Rate[0:124])).dot(Expenditure_pp[m,0:124]) \n",
    "    Expenditure_bp[m,154] = 0.32 * (Mar_Rate[0:124] * (1-Tax_Rate[0:124])).dot(Expenditure_pp[m,0:124])\n",
    "    #0.68 and 0.32 are the proportions of wholesale and retail trade in EXIOBASE for German households\n",
    "    #The proportions for other countries might differ than these values.\n",
    "\n",
    "print(Expenditure_bp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the expenditure survey data in basic price is split into 9800 region-products. First of all, the regional shares of household consumption for each consumption category (200 categories) are calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.         0.00465635 0.00766211 ... 0.00021632 0.         0.        ]\n"
     ]
    }
   ],
   "source": [
    "MRIO_Y_RegionalShares = np.zeros(9800)\n",
    "counter = 0\n",
    "for m in range(0,49):\n",
    "    for n in range(0,200):\n",
    "        if HH_FinalDemand_bp[n] == 0:\n",
    "            MRIO_Y_RegionalShares[counter] == 0\n",
    "        else:\n",
    "            MRIO_Y_RegionalShares[counter] = MRIO_Y[counter,NoofFinDemands*Country_ID]/HH_FinalDemand_bp[n]\n",
    "        counter = counter + 1\n",
    "\n",
    "print(MRIO_Y_RegionalShares)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that the households purchase their consumption similar to the national household consumption basket, their expenditure is split into 9800 region-products based on the previously calculated regional shares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(9800, 104)\n"
     ]
    }
   ],
   "source": [
    "EXIO_Expenditure = np.zeros((9800,104))\n",
    "for m in range(0,104):\n",
    "    Expenditure_Product = Expenditure_bp[m,:]\n",
    "    counter = 0\n",
    "    for n in range(0,49):\n",
    "        for o in range(0,200):\n",
    "            EXIO_Expenditure[counter,m] = Expenditure_Product[o] * MRIO_Y_RegionalShares[counter]\n",
    "            counter = counter+1\n",
    "            \n",
    "print(EXIO_Expenditure.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12. Step: Calculating the household Carbon Footprint\n",
    "\n",
    "Household carbon footprint shows the total GHG emissions embodied in the household consumption, measured in CO<sub>2</sub> equivalent. The carbon footprint is calculated using this formula.\n",
    "\n",
    "$$ ICF = C \\times (S. L. Y_i)$$\n",
    "\n",
    "Where ICF is the Indirect Carbon Footprint, C is the Characterization Factor matrix for Carbon Footprint, L is the Leontief-Inverse matrix representing the production recipe, and Y is the final demand matrix, which value here is the household expenditure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.24224218e+02 4.92452796e+02 1.02931680e+02 4.18730083e+02\n",
      " 8.15729186e+01 1.66086642e+02 2.13705890e+02 1.75395616e+02\n",
      " 9.86938175e+01 1.71604387e+02 7.61684404e+01 9.42967417e+01\n",
      " 0.00000000e+00 2.10963048e+01 9.44700408e+02 1.67797594e+01\n",
      " 4.63319194e+00 1.49125593e+02 1.24252103e+00 3.45261348e+02\n",
      " 6.74232972e+00 4.11979610e+02 3.29483660e+01 6.09579158e+01\n",
      " 8.37175707e+00 6.39588882e+02 1.61574640e+03 5.28880059e+02\n",
      " 4.99208386e+02 1.09952597e+02 1.73350477e+02 2.08918462e+02\n",
      " 2.96430305e+01 1.91718287e+00 8.84738500e+01 3.13791254e+01\n",
      " 1.90790566e+01 2.07694811e+00 5.04878434e+01 3.60088531e+01\n",
      " 1.05695760e+02 1.21871656e+01 2.70702967e+02 5.19231812e+01\n",
      " 7.78041838e+01 5.06711303e+01 2.94459369e+01 6.90476139e+00\n",
      " 1.51403637e+01 2.26359067e+02 6.61352635e+00 2.35486286e+01\n",
      " 0.00000000e+00 1.16878630e+01 1.18519500e+03 7.81236065e+01\n",
      " 1.85481273e+02 2.40132259e+01 2.96108586e+00 1.46427615e+02\n",
      " 3.60166267e+01 1.03110327e+02 2.12819219e+00 2.72153256e+00\n",
      " 2.78174377e+01 5.37801949e+01 7.20299533e+01 2.51143684e+01\n",
      " 9.24691475e+01 8.09301571e+00 1.44845127e+00 1.26064260e+01\n",
      " 0.00000000e+00 1.35382876e+00 8.83581850e+00 3.06983981e+01\n",
      " 8.72164161e+01 1.46183682e+02 0.00000000e+00 4.05404844e+01\n",
      " 4.97963028e+01 2.54854562e+01 3.41371729e+01 2.41342640e+01\n",
      " 3.76686046e+01 0.00000000e+00 4.28520870e+02 2.84015063e+01\n",
      " 1.35957771e+02 1.74549199e+01 4.81293346e+01 4.32678423e+01\n",
      " 3.90612032e+00 2.19616587e+02 0.00000000e+00 4.12380606e+01\n",
      " 8.99128311e+00 1.02956461e+01 0.00000000e+00 1.75021082e+02\n",
      " 6.18594932e+01 7.41718775e+01 7.72045517e+00 0.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "#Expenditure has to be converted from Euro to MEuro\n",
    "Footprint = MRIO_Char[4].dot(MRIO_S.dot(MRIO_L).dot(EXIO_Expenditure))/1e6\n",
    "print(Footprint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the footprint calculated before only considers the embodied emissions from the production stage and excludes the direct emissions from heating and transportation, household direct emissions need to be calculated additionally by multiplying the emissions intensity (in kg CO<sub>2</sub>/Euro) with the expenditure.\n",
    "\n",
    "The emissions intensity data for German households are available from Hardadi, et al. (2020). To calculate the emissions intensity for other countries, divide the emissions intensity in physical unit by the fuel prices. The emission intensities in physical unit are:\n",
    "\n",
    "Heating oil:\t11.8 kg CO<sub>2</sub>/gallon\n",
    "\n",
    "Heating gas:\t55.9 ton CO<sub>2</sub>/TJ\n",
    "\n",
    "Gasoline fuel:\t3.18 ton CO<sub>2</sub>/ton\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1921.9199999999998\n",
      "48.65\n",
      "1496.72\n",
      "1427.1000000000001\n",
      "4894.39\n"
     ]
    }
   ],
   "source": [
    "DirEmis_Fuel = 1.68 * CES_DEexpenditure[54] #Expenditure of transportation fuel\n",
    "DirEmis_Wood = 1.39 * CES_DEexpenditure[29] #Expenditure of solid heating fuel (wood or coal)\n",
    "DirEmis_HOil = 4.24 * CES_DEexpenditure[28] #Expenditure of heating oil\n",
    "DirEmis_HGas = 3.35 * CES_DEexpenditure[27] #Expenditure of heating gas\n",
    "Dir_Emission = DirEmis_Fuel + DirEmis_Wood + DirEmis_HOil + DirEmis_HGas\n",
    "\n",
    "print(DirEmis_Fuel)\n",
    "print(DirEmis_Wood)\n",
    "print(DirEmis_HOil)\n",
    "print(DirEmis_HGas)\n",
    "print(Dir_Emission)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the total household carbon footprint is calculated as the sum of indirect carbon footprint and direct emissions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18287.63  kg CO2e/household\n"
     ]
    }
   ],
   "source": [
    "HH_Footprint = Footprint.sum(axis=0) + Dir_Emission\n",
    "print(round(HH_Footprint, 2), \" kg CO2e/household\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References:\n",
    "\n",
    "Hardadi, G., Buchholz, A., Pauliuk, S. (2020). Implications of the Distribution of German Household Environmental Footprints across Income Groups for Integrating Environmental and Social Policy Design. Journal of Industrial Ecology, 00(0), 1–11. https://doi.org/10.1111/jiec.12670\n",
    "\n",
    "Miller, R. E., & Blair, P. D. (2009). Input-Output Analysis. Cambridge University Press.\n",
    "\n",
    "Min, J., & Rao, N. D. (2017). Estimating Uncertainty in Household Energy Footprints. Journal of Industrial Ecology, 00(0), 1–11. https://doi.org/10.1111/jiec.12670\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}