diff --git a/config.default.yaml b/config.default.yaml
index a19516d1..50c07f46 100644
--- a/config.default.yaml
+++ b/config.default.yaml
@@ -79,10 +79,15 @@ sector:
   'time_dep_hp_cop' : True
   'space_heating_fraction' : 1.0 #fraction of space heating active
   'retrofitting' :
-      'retro_exogen': True # space heat demand savings exogenously
+      'retro_exogen': False # space heat demand savings exogenously
       'dE': 0.4            # reduction of space heat demand (applied before losses in DH)
-      'retro_endogen': False  # co-optimise space heat savings
-      'retrofitting-cost_factor' : 1.0
+      'retro_endogen': True  # co-optimise space heat savings
+      'cost_factor' : 1.0
+      'interest_rate': 0.04  # for investment in building components
+      'annualise_cost': True  # annualise the investment costs
+      'tax_weighting': False   # weight costs depending on taxes in countries
+      'construction_index': True   # weight costs depending on labour/material costs per ct
+      'l_strength': ["0.076", "0.197"] # additional insulation thickness[m], determines number of retro steps(=generators per bus) and maximum possible savings
   'tes' : True
   'tes_tau' : 3.
   'boilers' : True
diff --git a/config.myopic.yaml b/config.myopic.yaml
index a0692ff5..fb7bf50a 100644
--- a/config.myopic.yaml
+++ b/config.myopic.yaml
@@ -82,7 +82,12 @@ sector:
       'retro_exogen': True # space heat demand savings exogenously
       'dE': [0.0, 0.1, 0.2, 0.5]  # reduction of space heat demand per year (applied before losses in DH)
       'retro_endogen': False  # co-optimise space heat savings
-      'retrofitting-cost_factor' : 1.0
+      'cost_factor' : 1.0
+      'interest_rate': 0.04  # for investment in building components
+      'annualise_cost': True  # annualise the investment costs
+      'tax_weighting': False   # weight costs depending on taxes in countries
+      'construction_index': True   # weight costs depending on labour/material costs per ct
+      'l_strength': ["0.076", "0.197"] # additional insulation thickness[m], determines number of retro steps(=generators per bus) and maximum possible savings
   'tes' : True
   'tes_tau' : 3.
   'boilers' : True
diff --git a/scripts/build_retro_cost.py b/scripts/build_retro_cost.py
new file mode 100644
index 00000000..a70ab0cc
--- /dev/null
+++ b/scripts/build_retro_cost.py
@@ -0,0 +1,477 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Jan 20 14:57:21 2020
+
+@author: bw0928
+
+*****************************************************************************
+This script calculates cost-energy_saving-curves for retrofitting
+for the EU-37 countries, based on the building stock data from hotmaps and
+the EU building stock database
+*****************************************************************************
+
+Structure:
+
+    (1) set assumptions and parameters
+    (2) read and prepare data
+    (3) calculate (€-dE-curves)
+    (4) save in csv
+
+*****************************************************************************
+"""
+
+import pandas as pd
+import matplotlib.pyplot as plt
+pd.options.mode.chained_assignment = None
+
+#%% ************ FUCNTIONS ***************************************************
+
+# windows ---------------------------------------------------------------
+def window_limit(l, window_assumptions):
+    """
+    define limit u value from which on window is retrofitted
+    """
+    m = (window_assumptions.diff()["u_limit"] /
+         window_assumptions.diff()["strength"]).dropna().iloc[0]
+    a = window_assumptions["u_limit"][0] - m * window_assumptions["strength"][0]
+    return m*l + a
+
+def u_retro_window(l, window_assumptions):
+    """
+    define retrofitting value depending on renovation strength
+    """
+    m = (window_assumptions.diff()["u_value"] /
+         window_assumptions.diff()["strength"]).dropna().iloc[0]
+    a = window_assumptions["u_value"][0] - m * window_assumptions["strength"][0]
+    return max(m*l + a, 0.8)
+
+def window_cost(u, cost_retro, window_assumptions):
+    """
+    get costs for new windows depending on u value
+
+    """
+    m = (window_assumptions.diff()["cost"] /
+         window_assumptions.diff()["u_value"]).dropna().iloc[0]
+    a = window_assumptions["cost"][0] - m * window_assumptions["u_value"][0]
+    window_cost = m*u + a
+    if annualise_cost:
+        window_cost = window_cost * interest_rate / (1 - (1 + interest_rate)
+                      ** -cost_retro.loc["Windows", "life_time"])
+    return window_cost
+
+#  functions for intermediate steps (~l, ~area)  -----------------------------
+def calculate_new_u(u_values, l, l_weight, k=0.035):
+    """
+    calculate U-values after building retrofitting, depending on the old
+    U-values (u_values).
+    They depend for the components Roof, Wall, Floor on the additional
+    insulation thickness (l), and the weighting for the corresponding
+    component (l_weight).
+    Windows are renovated to new ones with U-value (function: u_retro_window(l))
+    only if the are worse insulated than a certain limit value
+    (function: window_limit).
+
+    Parameters
+    ----------
+    u_values: pd.DataFrame
+    l: string
+    l_weight: pd.DataFrame (component, weight)
+    k: thermal conductivity
+
+    """
+    return u_values.apply(lambda x:
+                                 k / ((k / x.value) +
+                                      (float(l) * l_weight.loc[x.type][0]))
+                                 if x.type!="Windows"
+                             else (min(x.value, u_retro_window(float(l), window_assumptions))
+                                   if x.value>window_limit(float(l), window_assumptions) else x.value),
+                             axis=1)
+
+def calculate_dE(u_values, l, average_surface_w):
+    """
+    returns energy demand after retrofit (per unit of unrefurbished energy
+    demand) depending on current and retrofitted U-values, this energy demand
+    is weighted depending on the average surface of each component for the
+    building type of the assumend subsector
+    """
+    return u_values.apply(lambda x: x[l] / x.value *
+                                     average_surface_w.loc[x.assumed_subsector,
+                                                           x.type],
+                                     axis=1)
+
+
+def calculate_costs(u_values, l, cost_retro, average_surface):
+    """
+    returns costs for a given retrofitting strength weighted by the average
+    surface/volume ratio of the component for each building type
+    """
+    return u_values.apply(lambda x: (cost_retro.loc[x.type, "cost_var"] *
+                                     100 * float(l) * l_weight.loc[x.type][0]
+                                     + cost_retro.loc[x.type, "cost_fix"]) *
+                              average_surface.loc[x.assumed_subsector, x.type] /
+                              average_surface.loc[x.assumed_subsector, "surface"]
+                              if x.type!="Windows"
+                              else (window_cost(x[l], cost_retro, window_assumptions) *
+                                               average_surface.loc[x.assumed_subsector, x.type] /
+                                               average_surface.loc[x.assumed_subsector, "surface"]
+                                    if x.value>window_limit(float(l), window_assumptions) else 0),
+                             axis=1)
+
+
+# ---------------------------------------------------------------------------
+def prepare_building_stock_data():
+    """
+    reads building stock data and cleans up the format, returns
+    --------
+    u_values:          pd.DataFrame current U-values
+    average_surface:   pd.DataFrame (index= building type,
+                                     columns = [surface [m],height [m],
+                                                components area [m^2]])
+    average_surface_w: pd.DataFrame weighted share of the components per
+                       building type
+    area_tot:          heated floor area per country and sector [Mm²]
+    area:              heated floor area [Mm²] for country, sector, building
+                       type and period
+
+    """
+
+    building_data = pd.read_csv(snakemake.input.building_stock,
+                                usecols=list(range(13)))
+
+    # standardize data
+    building_data["type"].replace(
+                {'Covered area: heated  [Mm²]': 'Heated area [Mm²]',
+                 'Windows ': 'Windows',
+                 'Walls ': 'Walls',
+                 'Roof ': 'Roof',
+                 'Floor ': 'Floor'}, inplace=True)
+
+    building_data.country_code = building_data.country_code.str.upper()
+    building_data["subsector"].replace({'Hotels and Restaurants':
+                                        'Hotels and restaurants'}, inplace=True)
+    building_data["sector"].replace({'Residential sector': 'residential',
+                                     'Service sector': 'services'},
+                                    inplace=True)
+    # extract u-values
+    u_values = building_data[(building_data.feature.str.contains("U-values"))
+                             & (building_data.subsector != "Total")]
+
+    components = list(u_values.type.unique())
+
+    country_iso_dic = building_data.set_index("country")["country_code"].to_dict()
+
+    # add missing /rename countries
+    country_iso_dic.update({'Norway': 'NO',
+                            'Iceland': 'IS',
+                            'Montenegro': 'ME',
+                            'Serbia': 'RS',
+                            'Albania': 'AL',
+                            'United Kingdom': 'GB',
+                            'Bosnia and Herzegovina': 'BA',
+                            'Switzerland': 'CH'})
+
+    # heated floor area ----------------------------------------------------------
+    area = building_data[(building_data.type == 'Heated area [Mm²]') &
+                         (building_data.subsector != "Total")]
+    area_tot = area.groupby(["country", "sector"]).sum()
+    area["weight"] = area.apply(lambda x: x.value /
+                                          area_tot.value.loc[(x.country, x.sector)],
+                                axis=1)
+    area = area.groupby(['country', 'sector', 'subsector', 'bage']).sum()
+    area_tot.rename(index=country_iso_dic, inplace=True)
+
+    # add for some missing countries floor area from other data sources
+    area_missing = pd.read_csv(snakemake.input.floor_area_missing,
+                               index_col=[0, 1], usecols=[0, 1, 2, 3])
+    area_tot = area_tot.append(area_missing.unstack(level=-1).dropna().stack())
+    area_tot = area_tot.loc[~area_tot.index.duplicated(keep='last')]
+
+    # for still missing countries calculate floor area by population size
+    pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
+    pop_layout["ct"] = pop_layout.index.str[:2]
+    ct_total = pop_layout.total.groupby(pop_layout["ct"]).sum()
+
+    area_per_pop = area_tot.unstack().apply(lambda x: x / ct_total[x.index])
+    missing_area_ct = ct_total.index.difference(area_tot.index.levels[0])
+    for ct in missing_area_ct:
+        averaged_data = pd.DataFrame(
+            area_per_pop.value.reindex(map_for_missings[ct]).mean()
+            * ct_total[ct],
+            columns=["value"])
+        index = pd.MultiIndex.from_product([[ct], averaged_data.index.to_list()])
+        averaged_data.index = index
+        averaged_data["estimated"] = 1
+        if ct not in area_tot.index.levels[0]:
+            area_tot = area_tot.append(averaged_data, sort=True)
+        else:
+            area_tot.loc[averaged_data.index] = averaged_data
+
+    # u_values for Poland are missing -> take them from eurostat -----------
+    u_values_PL = pd.read_csv(snakemake.input.u_values_PL)
+    area_PL = area.loc["Poland"].reset_index()
+    data_PL = pd.DataFrame(columns=u_values.columns, index=area_PL.index)
+    data_PL["country"] = "Poland"
+    data_PL["country_code"] = "PL"
+    # data from area
+    for col in ["sector", "subsector", "bage"]:
+        data_PL[col] = area_PL[col]
+    data_PL["btype"] = area_PL["subsector"]
+
+    data_PL_final = pd.DataFrame()
+    for component in components:
+        data_PL["type"] = component
+        data_PL["value"] = data_PL.apply(lambda x: u_values_PL[(u_values_PL.component==component)
+                                                               & (u_values_PL.sector==x["sector"])]
+                                         [x["bage"]].iloc[0], axis=1)
+        data_PL_final = data_PL_final.append(data_PL)
+
+    u_values = pd.concat([u_values,
+                          data_PL_final]).reset_index(drop=True)
+
+    # clean data ---------------------------------------------------------------
+    # smallest possible today u values for windows 0.8 (passive house standard)
+    # maybe the u values for the glass and not the whole window including frame
+    # for those types assumed in the dataset
+    u_values[(u_values.type=="Windows") & (u_values.value<0.8)]["value"] = 0.8
+    # drop unnecessary columns
+    u_values.drop(['topic', 'feature','detail', 'estimated','unit'],
+                  axis=1, inplace=True, errors="ignore")
+    #  only take in config.yaml specified countries into account
+    countries = ct_total.index
+    area_tot = area_tot.loc[countries]
+
+    # average component surface --------------------------------------------------
+    average_surface = (pd.read_csv(snakemake.input.average_surface,
+                                   nrows=3,
+                                   header=1,
+                                   index_col=0).rename(
+                                       {'Single/two family house': 'Single family- Terraced houses',
+                                        'Large apartment house': 'Multifamily houses',
+                                        'Apartment house': 'Appartment blocks'},
+                                    axis="index")).iloc[:, :6]
+    average_surface.columns = ["surface", "height", "Roof",
+                               "Walls", "Floor", "Windows"]
+    # get area share of component
+    average_surface_w = average_surface[components].apply(lambda x: x / x.sum(),
+                                                          axis=1)
+
+    return (u_values, average_surface,
+            average_surface_w, area_tot, area, country_iso_dic, countries)
+
+
+def prepare_cost_retro():
+    """
+    read and prepare retro costs, annualises them if annualise_cost=True
+    """
+    cost_retro = pd.read_csv(snakemake.input.cost_germany,
+                         nrows=4, index_col=0, usecols=[0, 1, 2, 3])
+    cost_retro.index = cost_retro.index.str.capitalize()
+    cost_retro.rename(index={"Window": "Windows", "Wall": "Walls"}, inplace=True)
+
+    window_assumptions = pd.read_csv(snakemake.input.window_assumptions,
+                                     skiprows=[1], usecols=[0,1,2,3], nrows=2)
+
+    if annualise_cost:
+        cost_retro[["cost_fix", "cost_var"]] = (cost_retro[["cost_fix", "cost_var"]]
+                                                .apply(lambda x: x * interest_rate /
+                                                       (1 - (1 + interest_rate)
+                                                        ** -cost_retro.loc[x.index,
+                                                                           "life_time"])))
+
+    return cost_retro, window_assumptions
+
+
+def calculate_cost_energy_curve(u_values, l_strength, l_weight, average_surface_w,
+                                average_surface, area, country_iso_dic,
+                                countries):
+    """
+    returns energy demand per unit of unrefurbished (dE) and cost for given
+    renovation strength (l_strength), data for missing countries is
+    approximated by countries with similar building stock (dict:map_for_missings)
+
+    parameter
+    -------- input -----------
+    u_values:          pd.DataFrame current U-values
+    l_strength:        list of strings (strength of retrofitting)
+    l_weight:          pd.DataFrame (component, weight)
+    average_surface:   pd.DataFrame (index= building type,
+                                     columns = [surface [m],height [m],
+                                                components area [m^2]])
+    average_surface_w: pd.DataFrame weighted share of the components per
+                       building type
+    area:              heated floor area [Mm²] for country, sector, building
+                       type and period
+    country_iso_dic:   dict (maps country name to 2-letter-iso-code)
+    countries:         pd.Index (specified countries in config.yaml)
+    -------- output ----------
+    res:               pd.DataFrame(index=pd.MultiIndex([country, sector]),
+                                    columns=pd.MultiIndex([(dE/cost), l_strength]))
+    """
+
+    energy_saved = u_values[['country', 'sector', 'subsector', 'bage', 'type']]
+    costs = u_values[['country', 'sector', 'subsector', 'bage', 'type']]
+
+    for l in l_strength:
+        u_values[l] = calculate_new_u(u_values, l, l_weight)
+        energy_saved[l] = calculate_dE(u_values, l, average_surface_w)
+        costs[l] = calculate_costs(u_values, l, cost_retro, average_surface)
+
+    # energy and costs per country, sector, subsector and year
+    e_tot = energy_saved.groupby(['country', 'sector', 'subsector', 'bage']).sum()
+    cost_tot = costs.groupby(['country', 'sector', 'subsector', 'bage']).sum()
+
+    # weighting by area -> energy and costs per country and sector
+    # in case of missing data first concat
+    energy_saved = pd.concat([e_tot, area.weight], axis=1)
+    cost_res = pd.concat([cost_tot, area.weight], axis=1)
+    energy_saved = (energy_saved.apply(lambda x: x * x.weight, axis=1)
+                    .groupby(level=[0, 1]).sum())
+    cost_res = (cost_res.apply(lambda x: x * x.weight, axis=1)
+                .groupby(level=[0, 1]).sum())
+
+    res = pd.concat([energy_saved[l_strength], cost_res[l_strength]],
+                    axis=1, keys=["dE", "cost"])
+    res.rename(index=country_iso_dic, inplace=True)
+
+    res = res.loc[countries]
+
+    # map missing countries
+    for ct in map_for_missings.keys():
+        averaged_data = pd.DataFrame(res.loc[map_for_missings[ct], :].mean(level=1))
+        index = pd.MultiIndex.from_product([[ct], averaged_data.index.to_list()])
+        averaged_data.index = index
+        if ct not in res.index.levels[0]:
+            res = res.append(averaged_data)
+        else:
+            res.loc[averaged_data.index] = averaged_data
+
+    return res
+
+
+# %% **************** MAIN ************************************************
+if __name__ == "__main__":
+    #  for testing
+    if 'snakemake' not in globals():
+        import yaml
+        import os
+        from vresutils.snakemake import MockSnakemake
+        snakemake = MockSnakemake(
+            wildcards=dict(
+                network='elec',
+                simpl='',
+                clusters='37',
+                lv='1',
+                opts='Co2L-3H',
+                sector_opts="[Co2L0p0-168H-T-H-B-I]"),
+            input=dict(
+                building_stock="data/retro/data_building_stock.csv",
+                u_values_PL="data/retro/u_values_poland.csv",
+                tax_w="data/retro/electricity_taxes_eu.csv",
+                construction_index="data/retro/comparative_level_investment.csv",
+                average_surface="data/retro/average_surface_components.csv",
+                floor_area_missing="data/retro/floor_area_missing.csv",
+                clustered_pop_layout="resources/pop_layout_{network}_s{simpl}_{clusters}.csv",
+                cost_germany="data/retro/retro_cost_germany.csv",
+                window_assumptions="data/retro/window_assumptions.csv"),
+            output=dict(
+                retro_cost="resources/retro_cost_{network}_s{simpl}_{clusters}.csv",
+                floor_area="resources/floor_area_{network}_s{simpl}_{clusters}.csv")
+        )
+        with open('config.yaml', encoding='utf8') as f:
+            snakemake.config = yaml.safe_load(f)
+
+#  ******** (1) ASSUMPTIONS - PARAMETERS **********************************
+    retro_opts =  snakemake.config["sector"]["retrofitting"]
+    interest_rate = retro_opts["interest_rate"]
+    annualise_cost = retro_opts["annualise_cost"]  # annualise the investment costs
+    tax_weighting = retro_opts["tax_weighting"]   # weight costs depending on taxes in countries
+    construction_index = retro_opts["construction_index"]   # weight costs depending on labour/material costs per ct
+    # additional insulation thickness, determines maximum possible savings
+    l_strength = retro_opts["l_strength"]
+
+    k = 0.035   # thermal conductivity standard value
+    # strenght of relative retrofitting depending on the component
+    # determined by historical data of insulation thickness for retrofitting
+    l_weight = pd.DataFrame({"weight": [1.95, 1.48, 1.]},
+                            index=["Roof", "Walls", "Floor"])
+
+
+    # mapping missing countries by neighbours
+    map_for_missings = {
+        "AL": ["BG", "RO", "GR"],
+        "BA": ["HR"],
+        "RS": ["BG", "RO", "HR", "HU"],
+        "MK": ["BG", "GR"],
+        "ME": ["BA", "AL", "RS", "HR"],
+        "CH": ["SE", "DE"],
+        "NO": ["SE"],
+        }
+
+# %% ************ (2) DATA ***************************************************
+
+    # building stock data -----------------------------------------------------
+    (u_values, average_surface, average_surface_w,
+     area_tot, area, country_iso_dic, countries) = prepare_building_stock_data()
+
+    # costs for retrofitting -------------------------------------------------
+    cost_retro, window_assumptions = prepare_cost_retro()
+
+    # weightings of costs
+    if construction_index:
+        cost_w = pd.read_csv(snakemake.input.construction_index,
+                             skiprows=3, nrows=32, index_col=0)
+        # since German retrofitting costs are assumed
+        cost_w = ((cost_w["2018"] / cost_w.loc["Germany", "2018"])
+                  .rename(index=country_iso_dic))
+
+    if tax_weighting:
+        tax_w = pd.read_csv(snakemake.input.tax_w,
+                            header=12, nrows=39, index_col=0, usecols=[0, 4])
+        tax_w.rename(index=country_iso_dic, inplace=True)
+        tax_w = tax_w.apply(pd.to_numeric, errors='coerce').iloc[:, 0]
+        tax_w.dropna(inplace=True)
+
+# %% ********** (3) CALCULATE COST-ENERGY-CURVES ****************************
+
+    # for missing weighting of surfaces of building types assume MultiFamily houses
+    u_values["assumed_subsector"] = u_values.subsector
+    u_values.assumed_subsector[
+        ~u_values.subsector.isin(average_surface.index)] = 'Multifamily houses'
+
+    dE_and_cost =  calculate_cost_energy_curve(u_values, l_strength, l_weight,
+                                          average_surface_w, average_surface, area,
+                                          country_iso_dic, countries)
+
+    # weights costs after construction index
+    if construction_index:
+        for ct in list(map_for_missings.keys() - cost_w.index):
+            cost_w.loc[ct] = cost_w.reindex(index=map_for_missings[ct]).mean()
+        dE_and_cost.cost = dE_and_cost.cost.apply(lambda x: x * cost_w[x.index.levels[0]])
+
+    # weights cost depending on country taxes
+    if tax_weighting:
+        for ct in list(map_for_missings.keys() - tax_w.index):
+            tax_w[ct] = tax_w.reindex(index=map_for_missings[ct]).mean()
+        dE_and_cost.cost = dE_and_cost.cost.apply(lambda x: x * tax_w[x.index.levels[0]])
+
+    # get share of residential and sevice floor area
+    sec_w = (area_tot / area_tot.groupby(["country"]).sum())["value"]
+    # get the total cost-energy-savings weight by sector area
+    tot = dE_and_cost.apply(lambda col: col * sec_w, axis=0).groupby(level=0).sum()
+    tot.set_index(pd.MultiIndex.from_product([list(tot.index), ["tot"]]),
+                  inplace=True)
+    dE_and_cost = dE_and_cost.append(tot).unstack().stack()
+
+    summed_area = pd.DataFrame(area_tot.groupby("country").sum())
+    summed_area.set_index(pd.MultiIndex.from_product(
+                          [list(summed_area.index), ["tot"]]), inplace=True)
+    area_tot = area_tot.append(summed_area).unstack().stack()
+
+# %% ******* (4) SAVE   ************************************************
+
+    dE_and_cost.to_csv(snakemake.output.retro_cost)
+    area_tot.to_csv(snakemake.output.floor_area)
+
+
+