<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Microsimulation with artificial survey data
Step1:  Building the artificial data
Sénégal is composed by almost 15 millions people dispatched in around 1.6 million household
Step2:  We assume that 2/3 of the household heads are married and that only married houshold do have children. The mean number of children per household is 5 and is normally distributed
Step3:  We assume that 80% of the population are wage earners.
We choose a mean wage of 5 000 0000 CFA with a log normal ditribution.
Since 
$$ \text{mean wage}  = e^{\mu + \frac{\sigma ^ 2}{2}} $$ 
and
$$ \text{median wage} = e^\mu $$
we can compute the distribution according to the following expressions.
Step4:  We choose a mean pension of 2 500 000 CFA
Step5:  Microsimulation
As with test case, we can build a scenario with survey data
Step6:  We can compute the value of any variable for the whole population an draw distributions
Step7:  Special methods allow access to aggregates and pivot tables
Step8:  Evaluate the financial impact of a reform
Write a parametric reform tthat increseases the top marginal tax rates and evaluate how much revenue can be collected
<ASSISTANT_TASK:>
Python Code:
import matplotlib.pyplot as plt  # For graphics
%matplotlib inline
import numpy as np  # linear algebra and math
import pandas as pd  # data frames
from openfisca_core.model_api import *
from openfisca_senegal import SenegalTaxBenefitSystem  # The Senegalese tax-benefits system
from openfisca_senegal.survey_scenarios import SenegalSurveyScenario
Explanation: Microsimulation with artificial survey data
End of explanation
year = 2017
household_weight = 100
size = int(1.6e6 / household_weight)
print("Size of the sample: {}".format(size))
np.random.seed(seed = 42)
Explanation: Building the artificial data
Sénégal is composed by almost 15 millions people dispatched in around 1.6 million household
End of explanation
est_marie = np.random.binomial(1, .66, size = size)
est_celibataire = np.logical_not(est_marie)
nombre_enfants = np.maximum(
    est_marie * np.round(np.random.normal(5, scale = 3, size = size)),
    0,
    )
Explanation: We assume that 2/3 of the household heads are married and that only married houshold do have children. The mean number of children per household is 5 and is normally distributed
End of explanation
mean_wage = 5e6
median_wage = .75 * mean_wage
est_salarie = np.random.binomial(1, .8, size = size)
mu = np.log(median_wage)
sigma = np.sqrt(2 * np.log(mean_wage / median_wage))
salaire = (
    est_salarie * 
    np.random.lognormal(mean = mu, sigma = sigma, size = int(size))
    )
Explanation: We assume that 80% of the population are wage earners.
We choose a mean wage of 5 000 0000 CFA with a log normal ditribution.
Since 
$$ \text{mean wage}  = e^{\mu + \frac{\sigma ^ 2}{2}} $$ 
and
$$ \text{median wage} = e^\mu $$
we can compute the distribution according to the following expressions.
End of explanation
mean_pension = 2.5e6
median_pension = .9 * mean_pension
mu = np.log(median_pension)
sigma = np.sqrt(2 * np.log(mean_pension / median_pension))
pension_retraite = (
    np.logical_not(est_salarie) *
    np.random.lognormal(mean = mu, sigma = sigma, size = int(size))
    )
input_data_frame = pd.DataFrame({
    'est_marie': est_marie,
    'est_celibataire': est_celibataire,
    'nombre_enfants': nombre_enfants,
    'pension_retraite': pension_retraite,
    'salaire': salaire,
    'id_menage': range(size),
    'role_menage': 0,
    })
input_data_frame.salaire.hist(bins=100)
input_data_frame.pension_retraite.hist(bins=100, range = [.0001, 1e8])
Explanation: We choose a mean pension of 2 500 000 CFA
End of explanation
data = dict(input_data_frame = input_data_frame)
scenario = SenegalSurveyScenario(
    year = year,
    data = data,
    )
Explanation: Microsimulation
As with test case, we can build a scenario with survey data
End of explanation
scenario.simulation.calculate('impot_revenus', period = year)
df = pd.DataFrame({'impot_revenus': scenario.simulation.calculate('impot_revenus', period = year)})
df.hist(bins = 100)
Explanation: We can compute the value of any variable for the whole population an draw distributions
End of explanation
scenario.compute_aggregate('impot_revenus', period = year) / 1e9
scenario.compute_pivot_table(
    aggfunc = 'sum', 
    values = ['impot_revenus'], 
    columns = ['nombre_enfants'],
    period = year,
    ).stack().reset_index().plot(x = 'nombre_enfants', kind = 'bar')
Explanation: Special methods allow access to aggregates and pivot tables
End of explanation
def modify_parameters(parameters):
    parameters.bareme_impot_progressif[5].rate.update(period = period(year), value = .5)
    return parameters
    
class tax_the_rich(Reform):
    name = u"Tax last bracket at 50%"
    def apply(self):
        self.modify_parameters(modifier_function = modify_parameters)
senegal_tax_benefit_system = SenegalTaxBenefitSystem()
tax_the_rich_tax_benefit_system = tax_the_rich(senegal_tax_benefit_system)
scenario = SenegalSurveyScenario(
    data = data,
    tax_benefit_system = tax_the_rich_tax_benefit_system,
    baseline_tax_benefit_system = senegal_tax_benefit_system,
    year = year,
    )
print('reform tax the rich: ', scenario.compute_aggregate('impot_revenus', period = year) / 1e9)
print('baseline: ', scenario.compute_aggregate('impot_revenus', use_baseline = True, period = year) / 1e9)
from openfisca_senegal.entities import Person
def build_ultimate_reform_tax_benefit_system(seuil = 0, marginal_tax_rate = .4):
    senegal_tax_benefit_system = SenegalTaxBenefitSystem()
    class impot_revenus(Variable):
        def formula(individu, period):
            impot_avant_reduction_famille = individu('impot_avant_reduction_famille', period)
            reduction_impots_pour_charge_famille = individu('reduction_impots_pour_charge_famille', period)
            impot_apres_reduction_famille = impot_avant_reduction_famille - reduction_impots_pour_charge_famille
            impot_revenus = max_(0, impot_apres_reduction_famille)
            return impot_revenus * (impot_revenus > seuil)            
    def modify_parameters(parameters):
        parameters.bareme_impot_progressif[5].rate.update(period = period(year), value = marginal_tax_rate)
        return parameters
    class ultimate_reform(Reform):
        name = u"Tax the rich and save the poor taxpayers (tax < {})".format(seuil)
        def apply(self):
            self.update_variable(impot_revenus)
            self.modify_parameters(modifier_function = modify_parameters)
    return ultimate_reform(senegal_tax_benefit_system)
reformed_tax_benefit_system = build_ultimate_reform_tax_benefit_system(seuil = 100000, marginal_tax_rate = .45)
scenario = SenegalSurveyScenario(
    data = data,
    tax_benefit_system = reformed_tax_benefit_system,
    baseline_tax_benefit_system = SenegalTaxBenefitSystem(),
    year = year,
    )
print('reform: ', scenario.compute_aggregate('impot_revenus', period = year) / 1e9)
print('baseline: ', scenario.compute_aggregate('impot_revenus', use_baseline = True, period = year) / 1e9)
scenario = SenegalSurveyScenario(
    data = data,
    tax_benefit_system = reformed_tax_benefit_system,
    baseline_tax_benefit_system =  SenegalTaxBenefitSystem(),
    year = year,
    )
print(scenario.compute_aggregate('impot_revenus', period = year) / 1e9)
print(scenario.compute_aggregate('impot_revenus', use_baseline = True, period = year) / 1e9)
cost = - (
    scenario.compute_aggregate('impot_revenus', period = year) -
    scenario.compute_aggregate('impot_revenus', use_baseline = True, period = year)
    ) / 1e9
print(cost)
def compute_reform_cost(seuil):
    reformed_tax_benefit_system = build_ultimate_reform_tax_benefit_system(
        seuil = seuil, 
        marginal_tax_rate = .41
        )
    scenario = SenegalSurveyScenario(
        data = data,
        tax_benefit_system = reformed_tax_benefit_system,
        baseline_tax_benefit_system = SenegalTaxBenefitSystem(),
        year = year,
        )
    cost = - (
        scenario.compute_aggregate('impot_revenus', period = year) - 
        scenario.compute_aggregate('impot_revenus', use_baseline = True, period = year)
        ) / 1e9
    return cost
compute_reform_cost(seuil = 100000)
from scipy.optimize import fsolve
x = fsolve(compute_reform_cost, 100000)
x
compute_reform_cost(seuil = x)
reformed_tax_benefit_system = build_ultimate_reform_tax_benefit_system(
        seuil = x, 
        marginal_tax_rate = .41
        )
scenario = SenegalSurveyScenario(
    data = data,
    tax_benefit_system = reformed_tax_benefit_system,
    baseline_tax_benefit_system = SenegalTaxBenefitSystem(),
    year = year,
    )
scenario.compute_pivot_table(
    aggfunc = 'sum', 
    values = ['impot_revenus'], 
    columns = ['nombre_enfants'],
    period = year,
    difference = True,
    ).stack().reset_index().plot(x = 'nombre_enfants', kind = 'bar')
Explanation: Evaluate the financial impact of a reform
Write a parametric reform tthat increseases the top marginal tax rates and evaluate how much revenue can be collected
End of explanation
<END_TASK>