add wellbeing loss calculation

unbreakable/data/writer.py

@@ -3,9 +3,10 @@
 
 import pandas as pd
 import numpy as np
+import yaml
 
 
-def get_outcomes(households, tot_exposed_asset, expected_loss_frac, years_to_recover) -> dict:
+def get_outcomes(households: pd.DataFrame, tot_exposed_asset: float, expected_loss_frac: float, years_to_recover: int, consump_util: float) -> dict:
     '''Calculate outcomes of interest from the simulation model.
 
     Args:
@@ -50,13 +51,14 @@ def get_outcomes(households, tot_exposed_asset, expected_loss_frac, years_to_rec
     annual_average_consumption_loss, annual_average_consumption_loss_pct = calculate_average_annual_consumption_loss(
         affected_households, years_to_recover)
 
-    r = calculate_resilience(affected_households)
+    tot_wellbeing_loss = calc_tot_wellbeing_loss(households, consump_util)
+    r = calculate_resilience(affected_households, tot_wellbeing_loss)
 
     # Get the weighted average vulnerability by consumption quintile and decile
     weighted_vuln_quint = get_weighted_vuln(affected_households, quintile=True)
     weighted_vuln_dec = get_weighted_vuln(affected_households, quintile=False)
 
-    return {
+    outcomes = {
         'total_population': total_population,
         'total_asset_loss': total_asset_loss,
         'total_consumption_loss': total_consumption_loss,
@@ -77,11 +79,18 @@ def get_outcomes(households, tot_exposed_asset, expected_loss_frac, years_to_rec
         'annual_average_consumption_loss_pct': annual_average_consumption_loss_pct,
         'r': r,
         'mean_recovery_rate': mean_recovery_rate,
+        'tot_wellbeing_loss': tot_wellbeing_loss,
         'weighted_vuln_quint': weighted_vuln_quint,
         'weighted_vuln_dec': weighted_vuln_dec,
         'years_in_poverty': years_in_poverty
     }
 
+    # Save outcome names in a yaml file to pick up in preprocessing
+    with open('analysis/outcome_names.yaml', 'w') as f:
+        yaml.dump(list(outcomes.keys()), f)
+
+    return outcomes
+
 
 def find_poor(households: pd.DataFrame, poverty_line: float, years_to_recover: int) -> tuple:
     '''Get the poor at the beginning of the simulation and the poor at the end of the simulation
@@ -226,13 +235,14 @@ def calculate_average_annual_consumption_loss(affected_households: pd.DataFrame,
     return annual_average_consumption_loss, annual_average_consumption_loss_pct
 
 
-def calculate_resilience(affected_households: pd.DataFrame) -> float:
+def calculate_resilience(affected_households: pd.DataFrame, tot_wellbeing_loss: float) -> float:
     '''Calculate socio-economic resilience of affected households.
 
     Socio-economic resilience is a ratio of asset loss to consumption loss.
 
     Args:
         affected_households (pd.DataFrame): Affected households.
+        tot_wellbeing_loss (float): Total wellbeing loss.
 
     Returns:
         float: Socio-economic resilience
@@ -247,7 +257,8 @@ def calculate_resilience(affected_households: pd.DataFrame) -> float:
         r = np.nan
 
     else:
-        r = total_asset_damage / total_consumption_loss
+        # r = total_asset_damage / total_consumption_loss
+        r = total_asset_damage / tot_wellbeing_loss
 
     return r
 
@@ -277,3 +288,21 @@ def get_weighted_vuln(affected_households: pd.DataFrame, quintile: bool) -> dict
         pop_by_d = df.groupby('decile').sum(numeric_only=True)[['popwgt']]
         average_v_by_d = v_weighted_by_d['v_weighted'].div(pop_by_d['popwgt'])
         return average_v_by_d.to_dict()
+
+
+def calc_tot_wellbeing_loss(households: pd.DataFrame, consump_util: float) -> float:
+    '''Calculate wellbeing loss.
+
+    Args:
+        households (pd.DataFrame): Households.
+        consump_util (float): Consumption utility.
+
+    Returns:
+        float: Total wellbeing loss.
+    '''
+    x = households['aeexp'].multiply(
+        households['popwgt']) / households['popwgt'].sum()
+    W = x**(-consump_util)
+    wellbeing = households['wellbeing']
+    wellbeing_loss = - wellbeing / W
+    return wellbeing_loss.sum()

unbreakable/model.py

@@ -50,7 +50,7 @@ def model(**params) -> dict:
 
     # Uncertainties
     poverty_bias = params['poverty_bias']
-    consumption_utility = params['consumption_utility']
+    consump_util = params['consump_util']
     discount_rate = params['discount_rate']
     lambda_increment = params['lambda_increment']
     income_and_expenditure_growth = params['income_and_expenditure_growth']
@@ -99,16 +99,16 @@ def model(**params) -> dict:
                                 .pipe(calculate_exposure, poverty_bias, calc_exposure_params)
                                 .pipe(identify_affected, ident_affected_params)
                                 .pipe(apply_policy, my_policy)
-                                .pipe(calculate_recovery_rate, consumption_utility, discount_rate, lambda_increment, years_to_recover)
-                                .pipe(calculate_wellbeing, consumption_utility, discount_rate, income_and_expenditure_growth, years_to_recover, add_income_loss, cash_transfer))
+                                .pipe(calculate_recovery_rate, consump_util, discount_rate, lambda_increment, years_to_recover)
+                                .pipe(calculate_wellbeing, consump_util, discount_rate, income_and_expenditure_growth, years_to_recover, add_income_loss, cash_transfer))
 
         if save_households:
             Path(f'../experiments/households/').mkdir(parents=True, exist_ok=True)
             households.to_csv(
                 f'../experiments/households/{district}_{random_seed}.csv')
 
         array_outcomes = np.array(list(get_outcomes(
-            households, tot_exposed_asset, expected_loss_frac, years_to_recover).values()))
+            households, tot_exposed_asset, expected_loss_frac, years_to_recover, consump_util).values()))
 
         # * To check whether we have different households affected in different runs
         # if district == 'Castries':

unbreakable/modules/integrator.py

@@ -6,15 +6,15 @@
 import pickle
 
 
-def calculate_recovery_rate(households: pd.DataFrame, consumption_utility: float, discount_rate: float, lambda_increment: float, years_to_recover: int) -> pd.DataFrame:
+def calculate_recovery_rate(households: pd.DataFrame, consump_util: float, discount_rate: float, lambda_increment: float, years_to_recover: int) -> pd.DataFrame:
     '''Calculates the recovery rate for each affected household.
 
     Recovery rate is a function of the household vulnerability (v), 
     consumption utility, discount rate, median productivity and years to recover.    
 
     Args:
         households (pd.DataFrame): Households data frame.
-        consumption_utility (float): Consumption utility.
+        consump_util (float): Consumption utility.
         discount_rate (float): Discount rate.
         lambda_increment (float): Lambda increment for the integration.
         years_to_recover (int): Number of years to recover.
@@ -42,7 +42,7 @@ def calculate_recovery_rate(households: pd.DataFrame, consumption_utility: float
 
     # Calculate the recovery rate for each affected household
     affected_households['recovery_rate'] = affected_households['v'].apply(lambda x: integrate_and_find_recovery_rate(
-        x, results, consumption_utility, discount_rate, average_productivity, lambda_increment, years_to_recover))
+        x, results, consump_util, discount_rate, average_productivity, lambda_increment, years_to_recover))
 
     # Set recovery rate to zero for unaffected households
     households['recovery_rate'] = 0
@@ -54,13 +54,13 @@ def calculate_recovery_rate(households: pd.DataFrame, consumption_utility: float
     return households
 
 
-def integrate_and_find_recovery_rate(v: float, results: pd.DataFrame, consumption_utility: float, discount_rate: float, average_productivity: float, lambda_increment: float, years_to_recover: int) -> float:
+def integrate_and_find_recovery_rate(v: float, results: pd.DataFrame, consump_util: float, discount_rate: float, average_productivity: float, lambda_increment: float, years_to_recover: int) -> float:
     '''Find recovery rate (lambda) given the value of `v` (household vulnerability).
 
     Args:
         v (float): Household vulnerability.
         results (pd.DataFrame): Data frame to store integration results.
-        consumption_utility (float): Consumption utility.
+        consump_util (float): Consumption utility.
         discount_rate (float): Discount rate.
         average_productivity (float): Average productivity.
         lambda_increment (float): Lambda increment for the integration.
@@ -92,7 +92,7 @@ def integrate_and_find_recovery_rate(v: float, results: pd.DataFrame, consumptio
                 factor = average_productivity + _lambda
                 part1 = average_productivity - \
                     factor * v * np.e**(-_lambda * _t)
-                part1 = part1**(-consumption_utility)
+                part1 = part1**(-consump_util)
                 part2 = _t * factor - 1
                 part3 = np.e**(-_t * (discount_rate + _lambda))
                 dwdlambda += part1 * part2 * part3 * dt
@@ -107,12 +107,12 @@ def integrate_and_find_recovery_rate(v: float, results: pd.DataFrame, consumptio
             _lambda += lambda_increment
 
 
-def calculate_wellbeing(households: pd.DataFrame, consumption_utility: float, discount_rate: float, income_and_expenditure_growth: float, years_to_recover: int, add_income_loss: bool, cash_transfer: dict = {}) -> pd.DataFrame:
+def calculate_wellbeing(households: pd.DataFrame, consump_util: float, discount_rate: float, income_and_expenditure_growth: float, years_to_recover: int, add_income_loss: bool, cash_transfer: dict = {}) -> pd.DataFrame:
     '''Calculates consumption loss and well-being for each affected household.
 
     Args:
         households (pd.DataFrame): Households data frame.
-        consumption_utility (float): Consumption utility.
+        consump_util (float): Consumption utility.
         discount_rate (float): Discount rate.
         income_and_expenditure_growth (float): Income and expenditure growth.
         years_to_recover (int): Number of years to recover.
@@ -134,9 +134,8 @@ def calculate_wellbeing(households: pd.DataFrame, consumption_utility: float, di
                'net_consumption_loss_NPV',
                'c_t',
                'c_t_unaffected',
-               'w_final',
                'weeks_pov',
-               'w_final2']
+               'wellbeing']
 
     # Set values to zero
     households[columns] = 0
@@ -236,13 +235,12 @@ def calculate_wellbeing(households: pd.DataFrame, consumption_utility: float, di
         affected_households.loc[affected_households['c_t']
                                 < poverty_line_adjusted, 'weeks_pov'] += 1
 
-        # Integrate wellbeing
-        # !!!: Fix wellbeing integration, currently it is nan
-        affected_households['w_final'] += dt * (affected_households['c_t'])**(1 - consumption_utility) * \
-            np.e**(-discount_rate * _t) / (1 - consumption_utility)
-
-        affected_households['w_final2'] += affected_households['c_t_unaffected']**(1-consumption_utility)/(1-consumption_utility)*dt*(
-            (1-((affected_households['c_t_unaffected'] - affected_households['c_t'])/affected_households['c_t_unaffected'])*np.e**(-affected_households['recovery_rate']*_t))**(1-consumption_utility)-1)*np.e**(-discount_rate*_t)
+        # Integrate well-being
+        affected_households['wellbeing'] = affected_households['c_t_unaffected']**(1 - consump_util)\
+            / (1-consump_util) * dt\
+            * ((1 - ((affected_households['c_t_unaffected'] - affected_households['c_t']) / affected_households['c_t_unaffected'])
+                * np.e**(-affected_households['recovery_rate'] * _t))**(1 - consump_util) - 1)\
+            * np.e**(-discount_rate * _t)
 
         # Use to examine individual consumption recovery
         # Save consumption recovery value at the time _t