diff --git a/unbreakable/data/writer.py b/unbreakable/data/writer.py index 1b219f1..9e26ecc 100644 --- a/unbreakable/data/writer.py +++ b/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() diff --git a/unbreakable/model.py b/unbreakable/model.py index 95918fb..8426e52 100644 --- a/unbreakable/model.py +++ b/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,8 +99,8 @@ 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) @@ -108,7 +108,7 @@ def model(**params) -> dict: 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': diff --git a/unbreakable/modules/integrator.py b/unbreakable/modules/integrator.py index f81e61b..4eaf8a8 100644 --- a/unbreakable/modules/integrator.py +++ b/unbreakable/modules/integrator.py @@ -6,7 +6,7 @@ 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), @@ -14,7 +14,7 @@ def calculate_recovery_rate(households: pd.DataFrame, consumption_utility: float 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