notebooks/Fake-data-Senegal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Microsimulation with artificial survey data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt  # For graphics\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np  # linear algebra and math\n",
    "import pandas as pd  # data frames\n",
    "\n",
    "from openfisca_core.model_api import *\n",
    "from openfisca_senegal import CountryTaxBenefitSystem as SenegalTaxBenefitSystem  # The Senegalese tax-benefits system\n",
    "\n",
    "from openfisca_senegal.survey_scenarios import SenegalSurveyScenario"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%config Completer.use_jedi = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Building the artificial data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Sénégal is composed by almost 15 millions people dispatched in around 1.6 million household"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "household_weight = 100\n",
    "size = int(1.6e6 / household_weight)\n",
    "print(f\"Size of the sample: {size}\")\n",
    "np.random.seed(seed = 42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We assume that 2/3 of the household heads are married and that only married households do have children. The mean number of children per household is 5 and is normally distributed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "est_marie = np.random.binomial(1, .66, size = size)\n",
    "est_celibataire = np.logical_not(est_marie)\n",
    "nombre_enfants = np.maximum(\n",
    "    est_marie * np.round(np.random.normal(5, scale = 3, size = size)),\n",
    "    0,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We assume that 80% of the population are wage earners.\n",
    "We choose a mean wage of 5,000,000 CFA with a log normal ditribution.\n",
    "Since \n",
    "$$ \\text{mean wage}  = e^{\\mu + \\frac{\\sigma ^ 2}{2}} $$ \n",
    "and\n",
    "$$ \\text{median wage} = e^\\mu $$\n",
    "we can compute the distribution according to the following expressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mean_wage = 5e6\n",
    "median_wage = .75 * mean_wage\n",
    "est_salarie = np.random.binomial(1, .8, size = size)\n",
    "mu = np.log(median_wage)\n",
    "sigma = np.sqrt(2 * np.log(mean_wage / median_wage))\n",
    "salaire = (\n",
    "    est_salarie * \n",
    "    np.random.lognormal(mean = mu, sigma = sigma, size = int(size))\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We choose a mean pension of 2 500 000 CFA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mean_pension = 2.5e6\n",
    "median_pension = .9 * mean_pension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mu = np.log(median_pension)\n",
    "sigma = np.sqrt(2 * np.log(mean_pension / median_pension))\n",
    "pension_retraite = (\n",
    "    np.logical_not(est_salarie) *\n",
    "    np.random.lognormal(mean = mu, sigma = sigma, size = int(size))\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "input_data_frame = pd.DataFrame({\n",
    "    'est_marie': est_marie,\n",
    "    'est_celibataire': est_celibataire,\n",
    "    'nombre_enfants': nombre_enfants,\n",
    "    'pension_retraite': pension_retraite,\n",
    "    'salaire': salaire,\n",
    "    'household_id': range(size),\n",
    "    'household_role_index': 0,\n",
    "    'household_weight': household_weight,\n",
    "    'person_weight': household_weight,\n",
    "    })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "input_data_frame.salaire.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "input_data_frame.groupby(['est_marie', 'nombre_enfants'])['salaire'].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "input_data_frame.salaire.hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "input_data_frame.salaire.hist(bins = 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "input_data_frame.pension_retraite.hist(bins=100, range = [.0001, 1e7])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Microsimulation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Scenario"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "As with test case, we can build a scenario with survey data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "data = dict(input_data_frame = input_data_frame)\n",
    "scenario = SenegalSurveyScenario(\n",
    "    data = data,\n",
    "    year = 2017,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "scenario.tax_benefit_system.variables.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can compute the value of any variable for the whole population an draw distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.simulation.calculate('impot_revenus', period = 2017)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame(\n",
    "    {'impot_revenus': scenario.simulation.calculate('impot_revenus', period = 2017)}\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "df.hist(bins = 100, range = [.00001, 2.5e7])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Special scenario methods allow access to aggregates and pivot tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "scenario.weight_variable_by_entity = {'person': 'person_weight', 'household': 'household_weight'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "scenario.compute_aggregate('impot_revenus', period = 2017) / 1e9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "df = scenario.compute_pivot_table(\n",
    "    aggfunc = 'sum', \n",
    "    values = ['impot_revenus'], \n",
    "    columns = ['nombre_enfants'],\n",
    "    period = 2017,\n",
    "    )\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "df.transpose().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "scenario.compute_pivot_table(\n",
    "    aggfunc = 'sum', \n",
    "    values = ['impot_revenus'], \n",
    "    columns = ['nombre_enfants'],\n",
    "    period = 2017,\n",
    "    ).transpose().reset_index().plot(x = 'nombre_enfants', kind = 'bar')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Evaluate the financial impact of a reform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Parametric reform\n",
    "\n",
    "Write a parametric reform tthat increseases the top marginal tax rates and evaluate how much revenue can be collected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "year = 2017\n",
    "def modify_parameters(parameters):\n",
    "    parameters.prelevements_obligatoires.impots_directs.bareme_impot_progressif[5].rate.update(period = period(year), value = .5)\n",
    "    return parameters\n",
    "    \n",
    "class tax_the_rich(Reform):\n",
    "    name = \"Tax last bracket at 50%\"\n",
    "\n",
    "    def apply(self):\n",
    "        self.modify_parameters(modifier_function = modify_parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "senegal_tax_benefit_system = SenegalTaxBenefitSystem()\n",
    "tax_the_rich_tax_benefit_system = tax_the_rich(senegal_tax_benefit_system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario = SenegalSurveyScenario(\n",
    "    data = data,\n",
    "    tax_benefit_system = tax_the_rich_tax_benefit_system,\n",
    "    baseline_tax_benefit_system = senegal_tax_benefit_system,\n",
    "    year = 2017,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('reform tax the rich: ', scenario.compute_aggregate('impot_revenus', period = 2017) / 1e9)\n",
    "print('baseline: ', scenario.compute_aggregate('impot_revenus', period = 2017, use_baseline = True) / 1e9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.compute_aggregate('impot_revenus', period = 2017, difference = True) / 1e9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Structural reform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from openfisca_senegal.entities import Person\n",
    "\n",
    "def build_ultimate_reform_tax_benefit_system(threshold = 0, marginal_tax_rate = .4):\n",
    "    year = 2017\n",
    "\n",
    "    senegal_tax_benefit_system = SenegalTaxBenefitSystem()\n",
    "    class impot_revenus(Variable):\n",
    "        def formula(person, period):\n",
    "            impot_avant_reduction_famille = person('impot_avant_reduction_famille', period)\n",
    "            reduction_impots_pour_charge_famille = person('reduction_impots_pour_charge_famille', period)\n",
    "            impot_apres_reduction_famille = impot_avant_reduction_famille - reduction_impots_pour_charge_famille\n",
    "            impot_revenus = max_(0, impot_apres_reduction_famille)\n",
    "            return impot_revenus * (impot_revenus > threshold)            \n",
    "\n",
    "    def modify_parameters(parameters):\n",
    "        parameters.prelevements_obligatoires.impots_directs.bareme_impot_progressif[5].rate.update(\n",
    "            period = period(year), value = marginal_tax_rate)\n",
    "        return parameters\n",
    "\n",
    "\n",
    "    class ultimate_reform(Reform):\n",
    "        name = f\"Tax the rich and save the poor taxpayers (tax < {threshold})\"\n",
    "\n",
    "        def apply(self):\n",
    "            self.update_variable(impot_revenus)\n",
    "            self.modify_parameters(modifier_function = modify_parameters)\n",
    "\n",
    "    return ultimate_reform(senegal_tax_benefit_system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "reformed_tax_benefit_system = build_ultimate_reform_tax_benefit_system(\n",
    "    threshold = 100000, \n",
    "    marginal_tax_rate = .5\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario = SenegalSurveyScenario(\n",
    "    data = data,\n",
    "    tax_benefit_system = reformed_tax_benefit_system,\n",
    "    baseline_tax_benefit_system = SenegalTaxBenefitSystem(),\n",
    "    year = 2017\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print('reform: ', scenario.compute_aggregate('impot_revenus', period = 2017) / 1e9)\n",
    "print('baseline: ', scenario.compute_aggregate('impot_revenus', period = 2017, use_baseline = True) / 1e9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario = SenegalSurveyScenario(\n",
    "    data = data,\n",
    "    tax_benefit_system = reformed_tax_benefit_system,\n",
    "    baseline_tax_benefit_system =  SenegalTaxBenefitSystem(),\n",
    "    year = 2017\n",
    "    )\n",
    "print(scenario.compute_aggregate('impot_revenus', period = 2017) / 1e9)\n",
    "print(scenario.compute_aggregate('impot_revenus', period = 2017, use_baseline = True) / 1e9)\n",
    "cost = - (\n",
    "    scenario.compute_aggregate('impot_revenus', period = 2017) -\n",
    "    scenario.compute_aggregate('impot_revenus', period = 2017, use_baseline = True)\n",
    "    ) / 1e9\n",
    "print(cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Budget neutral reforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def compute_reform_cost(threshold = 0, marginal_tax_rate = .4):\n",
    "    reformed_tax_benefit_system = build_ultimate_reform_tax_benefit_system(\n",
    "        threshold = threshold, \n",
    "        marginal_tax_rate = float(marginal_tax_rate)  \n",
    "        # We need to convert to float here since fsolve use numpy array which are not accepted as \n",
    "        # legislation parameters\n",
    "        )\n",
    "    scenario = SenegalSurveyScenario(\n",
    "        data = data,\n",
    "        tax_benefit_system = reformed_tax_benefit_system,\n",
    "        baseline_tax_benefit_system = SenegalTaxBenefitSystem(),\n",
    "        year = 2017\n",
    "        )\n",
    "    cost = - (\n",
    "        scenario.compute_aggregate('impot_revenus', period = 2017) - \n",
    "        scenario.compute_aggregate('impot_revenus', period = 2017, use_baseline = True)\n",
    "        ) / 1e9\n",
    "    return cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def compute_reform_cost_from_marginal_tax_rate(marginal_tax_rate = .4):\n",
    "    return compute_reform_cost(threshold = 100000, marginal_tax_rate = marginal_tax_rate)\n",
    "\n",
    "def compute_reform_cost_from_threshold(threshold = 0):\n",
    "    return compute_reform_cost(threshold = threshold, marginal_tax_rate = .41)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import fsolve\n",
    "balancing_marginal_tax_rate = fsolve(compute_reform_cost_from_marginal_tax_rate, .40)\n",
    "balancing_threshold = fsolve(compute_reform_cost_from_threshold, 1000)\n",
    "print(f\"balancing_marginal_tax_rate to finance a 100000 threshold: {float(balancing_marginal_tax_rate)}\")\n",
    "print(f\"balancing_threshold financed by an increase from 40% to 41%: {float(balancing_threshold)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(compute_reform_cost_from_marginal_tax_rate(marginal_tax_rate = balancing_marginal_tax_rate))\n",
    "print(compute_reform_cost_from_threshold(threshold = balancing_threshold))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "reformed_tax_benefit_system = build_ultimate_reform_tax_benefit_system(\n",
    "    threshold = balancing_threshold, \n",
    "    marginal_tax_rate = .41\n",
    "    )\n",
    "scenario = SenegalSurveyScenario(\n",
    "    data = data,\n",
    "    tax_benefit_system = reformed_tax_benefit_system,\n",
    "    baseline_tax_benefit_system = SenegalTaxBenefitSystem(),\n",
    "    year = 2017\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.compute_pivot_table(\n",
    "    aggfunc = 'sum', \n",
    "    values = ['impot_revenus'], \n",
    "    columns = ['nombre_enfants'],\n",
    "    period = 2017,\n",
    "    difference = True,\n",
    "    ).stack().reset_index().plot(x = 'nombre_enfants', kind = 'bar')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "scenario.compute_pivot_table(\n",
    "    aggfunc = 'mean', \n",
    "    values = ['impot_revenus'], \n",
    "    columns = ['nombre_enfants'],\n",
    "    period = 2017,\n",
    "    difference = True,\n",
    "    ).stack().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.compute_pivot_table(\n",
    "    aggfunc = 'sum', \n",
    "    values = ['impot_revenus'], \n",
    "    columns = ['nombre_enfants'],\n",
    "    period = 2017,\n",
    "    difference = True,\n",
    "    ).stack().sum() / 1e9\n",
    "# The result is close to zéro and almost equal to the one computed above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Diaporama",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}