diff --git a/py/openfisca-analysis.ipynb b/py/openfisca-analysis.ipynb
new file mode 100644
index 0000000..2b56530
--- /dev/null
+++ b/py/openfisca-analysis.ipynb
@@ -0,0 +1,202 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "source": [
+    "from openfisca_us import Microsimulation, reforms\n",
+    "from openfisca_us.api import *\n",
+    "import numpy as np\n",
+    "import pandas as pd"
+   ],
+   "outputs": [],
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "source": [
+    "# Calculate poverty rate\n",
+    "baseline = Microsimulation(year=2020)\n",
+    "baseline.calc(\"in_poverty\", map_to=\"person\").mean()"
+   ],
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": [
+       "0.11652875394831709"
+      ]
+     },
+     "metadata": {},
+     "execution_count": 2
+    }
+   ],
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "source": [
+    "from openfisca_us.entities import *\n",
+    "\n",
+    "\n",
+    "class basic_income(Variable):\n",
+    "    value_type = float\n",
+    "    entity = TaxUnit\n",
+    "    definition_period = YEAR\n",
+    "    label = \"Basic income\"\n",
+    "\n",
+    "    def formula(tax_unit, period, parameters):\n",
+    "        # Extract FPG parameter\n",
+    "        fpg_params = parameters(period).poverty.fpg\n",
+    "        # Includes first_person and additional_person\n",
+    "        # Extract tax unit number of people\n",
+    "        nb_people = tax_unit.nb_persons()\n",
+    "        # Calculate FPG\n",
+    "        fpg = (\n",
+    "            fpg_params.first_person.contiguous_US\n",
+    "            + (nb_people - 1) * fpg_params.additional_person.contiguous_US\n",
+    "        )\n",
+    "        # Extract taxable income\n",
+    "        taxable_income = tax_unit(\"taxable_income\", period)\n",
+    "        # Calculate basic income phased out at 50%\n",
+    "        return np.maximum(fpg - taxable_income * 0.5, 0)\n",
+    "        # TODO: Make 50% a parameter"
+   ],
+   "outputs": [],
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "source": [
+    "class taxable_income(Variable):\n",
+    "    value_type = float\n",
+    "    entity = TaxUnit\n",
+    "    definition_period = YEAR"
+   ],
+   "outputs": [],
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "source": [
+    "class SPM_unit_net_income(Variable):\n",
+    "    value_type = float\n",
+    "    entity = SPMUnit\n",
+    "    label = u\"SPM unit net income\"\n",
+    "    definition_period = YEAR\n",
+    "\n",
+    "    def formula(spm_unit, period):\n",
+    "        INCOME_COMPONENTS = [\n",
+    "            \"SPM_unit_total_income\",\n",
+    "            \"SPM_unit_SNAP\",\n",
+    "            \"SPM_unit_capped_housing_subsidy\",\n",
+    "            \"SPM_unit_school_lunch_subsidy\",\n",
+    "            \"SPM_unit_energy_subsidy\",\n",
+    "            \"SPM_unit_WIC\",\n",
+    "        ]\n",
+    "        EXPENSE_COMPONENTS = [\n",
+    "            \"SPM_unit_FICA\",\n",
+    "            \"SPM_unit_federal_tax\",\n",
+    "            \"SPM_unit_state_tax\",\n",
+    "            \"SPM_unit_capped_work_childcare_expenses\",\n",
+    "            \"SPM_unit_medical_expenses\",\n",
+    "        ]\n",
+    "        income = add(spm_unit, period, *INCOME_COMPONENTS)\n",
+    "        expense = add(spm_unit, period, *EXPENSE_COMPONENTS)\n",
+    "        basic_income = sum_contained_tax_units(\n",
+    "            \"basic_income\", spm_unit, period\n",
+    "        )\n",
+    "        return income - expense + basic_income"
+   ],
+   "outputs": [],
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "source": [
+    "reform = (\n",
+    "    reforms.restructure(SPM_unit_net_income),\n",
+    "    reforms.new_variable(basic_income),\n",
+    "    reforms.new_variable(taxable_income),\n",
+    ")\n",
+    "reform_sim = Microsimulation(reform, year=2020)\n",
+    "from openfisca_us_data import RawCPS\n",
+    "\n",
+    "taxable_income_values = RawCPS.load(2020, \"tax_unit\").TAX_INC.values\n",
+    "reform_sim.simulation.set_input(\"taxable_income\", 2020, taxable_income_values)"
+   ],
+   "outputs": [],
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "source": [
+    "df = reform_sim.df(\n",
+    "    [\"SPM_unit_net_income\", \"taxable_income\", \"in_poverty\", \"basic_income\"]\n",
+    ")"
+   ],
+   "outputs": [],
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "source": [
+    "df.basic_income.groupby(df.taxable_income.decile_rank()).mean()"
+   ],
+   "outputs": [
+    {
+     "output_type": "execute_result",
+     "data": {
+      "text/plain": [
+       "1.0     17615.118505\n",
+       "2.0     18213.645279\n",
+       "3.0     17982.893939\n",
+       "4.0     12921.723542\n",
+       "5.0      7550.313452\n",
+       "6.0      4565.691330\n",
+       "7.0      3264.406229\n",
+       "8.0      2298.077669\n",
+       "9.0      1967.494285\n",
+       "10.0     1499.863151\n",
+       "dtype: float64"
+      ]
+     },
+     "metadata": {},
+     "execution_count": 8
+    }
+   ],
+   "metadata": {}
+  }
+ ],
+ "metadata": {
+  "orig_nbformat": 4,
+  "language_info": {
+   "name": "python",
+   "version": "3.8.8",
+   "mimetype": "text/x-python",
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "pygments_lexer": "ipython3",
+   "nbconvert_exporter": "python",
+   "file_extension": ".py"
+  },
+  "kernelspec": {
+   "name": "python3",
+   "display_name": "Python 3.8.8 64-bit ('base': conda)"
+  },
+  "interpreter": {
+   "hash": "d8fe82497dc3af1dafdfcaf67c3f347e622d9ec55d37e96a4812404db83e4772"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
\ No newline at end of file