python-dml-401k-IV.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "narrative-sailing",
   "metadata": {
    "papermill": {
     "duration": 0.040767,
     "end_time": "2021-03-24T14:13:51.082088",
     "exception": false,
     "start_time": "2021-03-24T14:13:51.041321",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Inference on Predictive and Causal Effects in High-Dimensional Nonlinear Models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ready-appearance",
   "metadata": {
    "papermill": {
     "duration": 0.036878,
     "end_time": "2021-03-24T14:13:51.155414",
     "exception": false,
     "start_time": "2021-03-24T14:13:51.118536",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Impact of 401(k) on  Financial Wealth\n",
    "\n",
    "As a practical illustration of the methods developed in this lecture, we consider estimation of the effect of 401(k) eligibility and participation \n",
    "on accumulated assets. 401(k) plans are pension accounts sponsored by employers. The key problem in determining the effect of participation in 401(k) plans on accumulated assets is saver heterogeneity coupled with the fact that the decision to enroll in a 401(k) is non-random. It is generally recognized that some people have a higher preference for saving than others. It also seems likely that those individuals with high unobserved preference for saving would be most likely to choose to participate in tax-advantaged retirement savings plans and would tend to have otherwise high amounts of accumulated assets. The presence of unobserved savings preferences with these properties then implies that conventional estimates that do not account for saver heterogeneity and endogeneity of participation will be biased upward, tending to overstate the savings effects of 401(k) participation.\n",
    "\n",
    "One can argue that eligibility for enrolling in a 401(k) plan in this data can be taken as exogenous after conditioning on a few observables of which the most important for their argument is income. The basic idea is that, at least around the time 401(k)’s initially became available, people were unlikely to be basing their employment decisions on whether an employer offered a 401(k) but would instead focus on income and other aspects of the job. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "divine-phoenix",
   "metadata": {
    "papermill": {
     "duration": 0.036497,
     "end_time": "2021-03-24T14:13:51.229370",
     "exception": false,
     "start_time": "2021-03-24T14:13:51.192873",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Data\n",
    "\n",
    "The data set can be downloaded from the github repo\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "75662e18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import relevant packages\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.model_selection import cross_val_score, cross_val_predict, KFold\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.linear_model import LassoCV, RidgeCV, ElasticNetCV, LinearRegression, Ridge, Lasso, LogisticRegressionCV\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier\n",
    "from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier\n",
    "from sklearn.tree import DecisionTreeRegressor, DecisionTreeClassifier\n",
    "import patsy\n",
    "import warnings\n",
    "from sklearn.base import BaseEstimator, clone\n",
    "import statsmodels.api as sm\n",
    "from IPython.display import Markdown\n",
    "import wget\n",
    "import os\n",
    "import seaborn as sns\n",
    "warnings.simplefilter('ignore')\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6e66e77f",
   "metadata": {},
   "outputs": [],
   "source": [
    "file = \"https://raw.githubusercontent.com/CausalAIBook/MetricsMLNotebooks/main/data/401k.csv\"\n",
    "data = pd.read_csv(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "65a2d086",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ira</th>\n",
       "      <th>a401</th>\n",
       "      <th>hval</th>\n",
       "      <th>hmort</th>\n",
       "      <th>hequity</th>\n",
       "      <th>nifa</th>\n",
       "      <th>net_nifa</th>\n",
       "      <th>tfa</th>\n",
       "      <th>net_tfa</th>\n",
       "      <th>tfa_he</th>\n",
       "      <th>...</th>\n",
       "      <th>i3</th>\n",
       "      <th>i4</th>\n",
       "      <th>i5</th>\n",
       "      <th>i6</th>\n",
       "      <th>i7</th>\n",
       "      <th>a1</th>\n",
       "      <th>a2</th>\n",
       "      <th>a3</th>\n",
       "      <th>a4</th>\n",
       "      <th>a5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>...</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>3462.871609</td>\n",
       "      <td>3990.459304</td>\n",
       "      <td>63595.865759</td>\n",
       "      <td>30022.505093</td>\n",
       "      <td>33573.360666</td>\n",
       "      <td>1.392864e+04</td>\n",
       "      <td>1.041415e+04</td>\n",
       "      <td>2.156603e+04</td>\n",
       "      <td>1.805153e+04</td>\n",
       "      <td>5.162490e+04</td>\n",
       "      <td>...</td>\n",
       "      <td>0.209178</td>\n",
       "      <td>0.172668</td>\n",
       "      <td>0.121432</td>\n",
       "      <td>0.158548</td>\n",
       "      <td>0.077358</td>\n",
       "      <td>0.144024</td>\n",
       "      <td>0.209682</td>\n",
       "      <td>0.297428</td>\n",
       "      <td>0.215431</td>\n",
       "      <td>0.133434</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>9648.027853</td>\n",
       "      <td>12825.840676</td>\n",
       "      <td>73679.374267</td>\n",
       "      <td>40056.879558</td>\n",
       "      <td>51830.770524</td>\n",
       "      <td>5.490488e+04</td>\n",
       "      <td>5.602886e+04</td>\n",
       "      <td>6.256504e+04</td>\n",
       "      <td>6.352250e+04</td>\n",
       "      <td>9.325339e+04</td>\n",
       "      <td>...</td>\n",
       "      <td>0.406742</td>\n",
       "      <td>0.377979</td>\n",
       "      <td>0.326645</td>\n",
       "      <td>0.365272</td>\n",
       "      <td>0.267171</td>\n",
       "      <td>0.351132</td>\n",
       "      <td>0.407102</td>\n",
       "      <td>0.457150</td>\n",
       "      <td>0.411142</td>\n",
       "      <td>0.340061</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>-40000.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-5.023020e+05</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-5.023020e+05</td>\n",
       "      <td>-5.023020e+05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.000000e+02</td>\n",
       "      <td>-1.350000e+03</td>\n",
       "      <td>3.750000e+02</td>\n",
       "      <td>-5.000000e+02</td>\n",
       "      <td>2.000000e+02</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>48000.000000</td>\n",
       "      <td>7000.000000</td>\n",
       "      <td>10000.000000</td>\n",
       "      <td>1.635000e+03</td>\n",
       "      <td>2.700000e+02</td>\n",
       "      <td>3.325000e+03</td>\n",
       "      <td>1.499000e+03</td>\n",
       "      <td>1.850000e+04</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>200.000000</td>\n",
       "      <td>95000.000000</td>\n",
       "      <td>50000.000000</td>\n",
       "      <td>46650.000000</td>\n",
       "      <td>8.765500e+03</td>\n",
       "      <td>6.543500e+03</td>\n",
       "      <td>1.908450e+04</td>\n",
       "      <td>1.652450e+04</td>\n",
       "      <td>6.920000e+04</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>100000.000000</td>\n",
       "      <td>153000.000000</td>\n",
       "      <td>300000.000000</td>\n",
       "      <td>150000.000000</td>\n",
       "      <td>300000.000000</td>\n",
       "      <td>1.430298e+06</td>\n",
       "      <td>1.430298e+06</td>\n",
       "      <td>1.536798e+06</td>\n",
       "      <td>1.536798e+06</td>\n",
       "      <td>1.687115e+06</td>\n",
       "      <td>...</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>8 rows × 44 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                 ira           a401           hval          hmort  \\\n",
       "count    9915.000000    9915.000000    9915.000000    9915.000000   \n",
       "mean     3462.871609    3990.459304   63595.865759   30022.505093   \n",
       "std      9648.027853   12825.840676   73679.374267   40056.879558   \n",
       "min         0.000000       0.000000       0.000000       0.000000   \n",
       "25%         0.000000       0.000000       0.000000       0.000000   \n",
       "50%         0.000000       0.000000   48000.000000    7000.000000   \n",
       "75%         0.000000     200.000000   95000.000000   50000.000000   \n",
       "max    100000.000000  153000.000000  300000.000000  150000.000000   \n",
       "\n",
       "             hequity          nifa      net_nifa           tfa       net_tfa  \\\n",
       "count    9915.000000  9.915000e+03  9.915000e+03  9.915000e+03  9.915000e+03   \n",
       "mean    33573.360666  1.392864e+04  1.041415e+04  2.156603e+04  1.805153e+04   \n",
       "std     51830.770524  5.490488e+04  5.602886e+04  6.256504e+04  6.352250e+04   \n",
       "min    -40000.000000  0.000000e+00 -5.023020e+05  0.000000e+00 -5.023020e+05   \n",
       "25%         0.000000  2.000000e+02 -1.350000e+03  3.750000e+02 -5.000000e+02   \n",
       "50%     10000.000000  1.635000e+03  2.700000e+02  3.325000e+03  1.499000e+03   \n",
       "75%     46650.000000  8.765500e+03  6.543500e+03  1.908450e+04  1.652450e+04   \n",
       "max    300000.000000  1.430298e+06  1.430298e+06  1.536798e+06  1.536798e+06   \n",
       "\n",
       "             tfa_he  ...           i3           i4           i5           i6  \\\n",
       "count  9.915000e+03  ...  9915.000000  9915.000000  9915.000000  9915.000000   \n",
       "mean   5.162490e+04  ...     0.209178     0.172668     0.121432     0.158548   \n",
       "std    9.325339e+04  ...     0.406742     0.377979     0.326645     0.365272   \n",
       "min   -5.023020e+05  ...     0.000000     0.000000     0.000000     0.000000   \n",
       "25%    2.000000e+02  ...     0.000000     0.000000     0.000000     0.000000   \n",
       "50%    1.850000e+04  ...     0.000000     0.000000     0.000000     0.000000   \n",
       "75%    6.920000e+04  ...     0.000000     0.000000     0.000000     0.000000   \n",
       "max    1.687115e+06  ...     1.000000     1.000000     1.000000     1.000000   \n",
       "\n",
       "                i7           a1           a2           a3           a4  \\\n",
       "count  9915.000000  9915.000000  9915.000000  9915.000000  9915.000000   \n",
       "mean      0.077358     0.144024     0.209682     0.297428     0.215431   \n",
       "std       0.267171     0.351132     0.407102     0.457150     0.411142   \n",
       "min       0.000000     0.000000     0.000000     0.000000     0.000000   \n",
       "25%       0.000000     0.000000     0.000000     0.000000     0.000000   \n",
       "50%       0.000000     0.000000     0.000000     0.000000     0.000000   \n",
       "75%       0.000000     0.000000     0.000000     1.000000     0.000000   \n",
       "max       1.000000     1.000000     1.000000     1.000000     1.000000   \n",
       "\n",
       "                a5  \n",
       "count  9915.000000  \n",
       "mean      0.133434  \n",
       "std       0.340061  \n",
       "min       0.000000  \n",
       "25%       0.000000  \n",
       "50%       0.000000  \n",
       "75%       0.000000  \n",
       "max       1.000000  \n",
       "\n",
       "[8 rows x 44 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "00884061",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ira</th>\n",
       "      <th>a401</th>\n",
       "      <th>hval</th>\n",
       "      <th>hmort</th>\n",
       "      <th>hequity</th>\n",
       "      <th>nifa</th>\n",
       "      <th>net_nifa</th>\n",
       "      <th>tfa</th>\n",
       "      <th>net_tfa</th>\n",
       "      <th>tfa_he</th>\n",
       "      <th>...</th>\n",
       "      <th>i3</th>\n",
       "      <th>i4</th>\n",
       "      <th>i5</th>\n",
       "      <th>i6</th>\n",
       "      <th>i7</th>\n",
       "      <th>a1</th>\n",
       "      <th>a2</th>\n",
       "      <th>a3</th>\n",
       "      <th>a4</th>\n",
       "      <th>a5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>69000</td>\n",
       "      <td>60150</td>\n",
       "      <td>8850</td>\n",
       "      <td>100</td>\n",
       "      <td>-3300</td>\n",
       "      <td>100</td>\n",
       "      <td>-3300</td>\n",
       "      <td>5550</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>78000</td>\n",
       "      <td>20000</td>\n",
       "      <td>58000</td>\n",
       "      <td>61010</td>\n",
       "      <td>61010</td>\n",
       "      <td>61010</td>\n",
       "      <td>61010</td>\n",
       "      <td>119010</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1800</td>\n",
       "      <td>0</td>\n",
       "      <td>200000</td>\n",
       "      <td>15900</td>\n",
       "      <td>184100</td>\n",
       "      <td>7549</td>\n",
       "      <td>7049</td>\n",
       "      <td>9349</td>\n",
       "      <td>8849</td>\n",
       "      <td>192949</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2487</td>\n",
       "      <td>-6013</td>\n",
       "      <td>2487</td>\n",
       "      <td>-6013</td>\n",
       "      <td>-6013</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>300000</td>\n",
       "      <td>90000</td>\n",
       "      <td>210000</td>\n",
       "      <td>10625</td>\n",
       "      <td>-2375</td>\n",
       "      <td>10625</td>\n",
       "      <td>-2375</td>\n",
       "      <td>207625</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 44 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    ira  a401    hval  hmort  hequity   nifa  net_nifa    tfa  net_tfa  \\\n",
       "0     0     0   69000  60150     8850    100     -3300    100    -3300   \n",
       "1     0     0   78000  20000    58000  61010     61010  61010    61010   \n",
       "2  1800     0  200000  15900   184100   7549      7049   9349     8849   \n",
       "3     0     0       0      0        0   2487     -6013   2487    -6013   \n",
       "4     0     0  300000  90000   210000  10625     -2375  10625    -2375   \n",
       "\n",
       "   tfa_he  ...  i3  i4  i5  i6  i7  a1  a2  a3  a4  a5  \n",
       "0    5550  ...   1   0   0   0   0   0   1   0   0   0  \n",
       "1  119010  ...   0   1   0   0   0   0   0   0   1   0  \n",
       "2  192949  ...   0   0   0   1   0   0   0   0   1   0  \n",
       "3   -6013  ...   0   0   1   0   0   1   0   0   0   0  \n",
       "4  207625  ...   0   1   0   0   0   0   0   1   0   0  \n",
       "\n",
       "[5 rows x 44 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1d3dc85c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r",
      "  0% [                                                                                ]    0 / 2314\r",
      "100% [................................................................................] 2314 / 2314"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "# Pension 401(k) data set\n",
       "\n",
       "## Description\n",
       "Data set on financial wealth and 401(k) plan participation. The source of this data set is the `hdm` package in `R`. It was created on 2/20/2023, by calling the `data(pension)` command in `R` after loading the `hdm` package and then writing the created dataframe to a csv file.\n",
       "\n",
       "## Format\n",
       "Dataframe with the following variables (amongst others):\n",
       "\n",
       "**p401**\n",
       "\n",
       "participation in 401(k)\n",
       "\n",
       "**e401**\n",
       "\n",
       "eligibility for 401(k)\n",
       "\n",
       "**a401**\n",
       "\n",
       "401(k) assets\n",
       "\n",
       "**tw**\n",
       "\n",
       "total wealth (in US $)\n",
       "\n",
       "**tfa**\n",
       "\n",
       "financial assets (in US $)\n",
       "\n",
       "**net_tfa**\n",
       "\n",
       "net financial assets (in US $)\n",
       "\n",
       "**nifa**\n",
       "\n",
       "non-401k financial assets (in US $)\n",
       "\n",
       "**net_nifa**\n",
       "\n",
       "net non-401k financial assets\n",
       "\n",
       "**net_n401**\n",
       "\n",
       "net non-401(k) assets (in US $)\n",
       "\n",
       "**ira**\n",
       "\n",
       "individual retirement account (IRA)\n",
       "\n",
       "**inc**\n",
       "\n",
       "income (in US $)\n",
       "\n",
       "**age**\n",
       "\n",
       "age\n",
       "\n",
       "**fsize**\n",
       "\n",
       "family size\n",
       "\n",
       "**marr**\n",
       "\n",
       "married\n",
       "\n",
       "**pira**\n",
       "\n",
       "participation in IRA\n",
       "\n",
       "**db**\n",
       "\n",
       "defined benefit pension\n",
       "\n",
       "**hown**\n",
       "\n",
       "home owner\n",
       "\n",
       "**educ**\n",
       "\n",
       "education (in years)\n",
       "\n",
       "**male**\n",
       "\n",
       "male\n",
       "\n",
       "**twoearn**\n",
       "\n",
       "two earners\n",
       "\n",
       "**nohs, hs, smcol, col**\n",
       "\n",
       "dummies for education: no high-school, high-school, some college, college\n",
       "\n",
       "**hmort**\n",
       "\n",
       "home mortage (in US $)\n",
       "\n",
       "**hequity**\n",
       "\n",
       "home equity (in US $)\n",
       "\n",
       "**hval**\n",
       "\n",
       "home value (in US $)\n",
       "\n",
       "## Details\n",
       "The sample is drawn from the 1991 Survey of Income and Program Participation (SIPP) and consists of 9,915 observations. The observational units are household reference persons aged 25-64 and spouse if present. Households are included in the sample if at least one person is employed and no one is self-employed. The data set was analysed in Chernozhukov and Hansen (2004) and Belloni et al. (2014) where further details can be found. They examine the effects of 401(k) plans on wealth using data from the Survey of Income and Program Participation using 401(k) eligibility as an instrument for 401(k) participation.\n",
       "\n",
       "## References\n",
       "V. Chernohukov, C. Hansen (2004). The impact of 401(k) participation on the wealth distribution: An instrumental quantile regression analysis. The Review of Economic and Statistics 86 (3), 735â€“751.\n",
       "\n",
       "A. Belloni, V. Chernozhukov, I. Fernandez-Val, and C. Hansen (2014). Program evaluation with high-dimensional data. Working Paper."
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "readme = \"https://raw.githubusercontent.com/CausalAIBook/MetricsMLNotebooks/main/data/401k.md\"\n",
    "filename = wget.download(readme)\n",
    "display(Markdown(open(filename, 'r').read()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "looking-invention",
   "metadata": {
    "papermill": {
     "duration": 0.037509,
     "end_time": "2021-03-24T14:13:52.250583",
     "exception": false,
     "start_time": "2021-03-24T14:13:52.213074",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The data consist of 9,915 observations at the household level drawn from the 1991 Survey of Income and Program Participation (SIPP).  All the variables are referred to 1990. We use net financial assets (*net\\_tfa*) as the outcome variable, $Y$,  in our analysis. The net financial assets are computed as the sum of IRA balances, 401(k) balances, checking accounts, saving bonds, other interest-earning accounts, other interest-earning assets, stocks, and mutual funds less non mortgage debts. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "received-nutrition",
   "metadata": {
    "papermill": {
     "duration": 0.038843,
     "end_time": "2021-03-24T14:13:52.327055",
     "exception": false,
     "start_time": "2021-03-24T14:13:52.288212",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Among the $9915$ individuals, $3682$ are eligible to participate in the program. The variable *e401* indicates eligibility and *p401* indicates participation, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0d26d552",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGwCAYAAABIC3rIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAo7ElEQVR4nO3df3RU9Z3/8ddIyBgwuRJCZhwJNLRZNpjgj+AJk1ZBgQC7MevxHHEbHekx5cdiYUeg2MDWIqdNKq4B1xxZYG1RxKXn1E3b3dqUuC0pGAIxx5wCItU1u8CSIdGdTBKMExrm+4flfjsEEWOSSfg8H+fMOc6977nzuZwT8zx3fsQRiUQiAgAAMNg1sV4AAABArBFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADBeXKwXMFycP39ep0+fVmJiohwOR6yXAwAArkAkElFHR4c8Ho+uuebTrwMRRFfo9OnTSktLi/UyAABAH5w8eVLjx4//1P0E0RVKTEyU9Mk/aFJSUoxXAwAArkR7e7vS0tLs3+OfhiC6QhdeJktKSiKIAAAYZj7r7S68qRoAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPHiYr0ARMv59kuxXgIw5DQ8/XCslwDgKscVIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxYh5E//u//6uHHnpIY8eO1ahRo3TLLbeooaHB3h+JRLR+/Xp5PB4lJCRo5syZOnr0aNQxwuGwli9frpSUFI0ePVqFhYU6depU1EwwGJTP55NlWbIsSz6fT21tbYNxigAAYIiLaRAFg0F99atf1ciRI/WrX/1Kb7/9tp555hldf/319szGjRtVXl6uiooK1dfXy+12a86cOero6LBn/H6/KisrtXv3bu3fv1+dnZ0qKChQT0+PPVNUVKTGxkZVVVWpqqpKjY2N8vl8g3m6AABgiHJEIpFIrJ78O9/5jt544w3t27fvkvsjkYg8Ho/8fr8ef/xxSZ9cDXK5XHrqqae0ZMkShUIhjRs3Tjt37tQDDzwgSTp9+rTS0tL02muvae7cuTp27JimTJmiuro65ebmSpLq6urk9Xr1zjvvaPLkyZ+51vb2dlmWpVAopKSkpH76F+iNP90B9Maf7gDQV1f6+zumV4h+8YtfaNq0abr//vuVmpqqW2+9Vdu3b7f3NzU1KRAIKD8/397mdDo1Y8YM1dbWSpIaGhp07ty5qBmPx6OsrCx75sCBA7Isy44hSZo+fbosy7JnLhYOh9Xe3h51AwAAV6eYBtH777+vLVu2KCMjQ7/+9a+1dOlSrVixQi+99MlVkkAgIElyuVxRj3O5XPa+QCCg+Ph4jRkz5rIzqampvZ4/NTXVnrlYWVmZ/X4jy7KUlpb2xU4WAAAMWTENovPnz+u2225TaWmpbr31Vi1ZskSLFi3Sli1bouYcDkfU/Ugk0mvbxS6eudT85Y5TUlKiUChk306ePHmlpwUAAIaZmAbRDTfcoClTpkRty8zM1IkTJyRJbrdbknpdxWlpabGvGrndbnV3dysYDF525syZM72ev7W1tdfVpwucTqeSkpKibgAA4OoU0yD66le/quPHj0dt+8Mf/qCJEydKktLT0+V2u1VdXW3v7+7uVk1NjfLy8iRJOTk5GjlyZNRMc3Ozjhw5Ys94vV6FQiEdOnTInjl48KBCoZA9AwAAzBUXyyd/7LHHlJeXp9LSUi1YsECHDh3Stm3btG3bNkmfvMzl9/tVWlqqjIwMZWRkqLS0VKNGjVJRUZEkybIsFRcXa9WqVRo7dqySk5O1evVqZWdna/bs2ZI+ueo0b948LVq0SFu3bpUkLV68WAUFBVf0CTMAAHB1i2kQ3X777aqsrFRJSYk2bNig9PR0bd68WQ8++KA9s2bNGnV1dWnZsmUKBoPKzc3Vnj17lJiYaM9s2rRJcXFxWrBggbq6ujRr1izt2LFDI0aMsGd27dqlFStW2J9GKywsVEVFxeCdLAAAGLJi+j1EwwnfQwTEDt9DBKCvhsX3EAEAAAwFBBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwXkyDaP369XI4HFE3t9tt749EIlq/fr08Ho8SEhI0c+ZMHT16NOoY4XBYy5cvV0pKikaPHq3CwkKdOnUqaiYYDMrn88myLFmWJZ/Pp7a2tsE4RQAAMAzE/ArRTTfdpObmZvt2+PBhe9/GjRtVXl6uiooK1dfXy+12a86cOero6LBn/H6/KisrtXv3bu3fv1+dnZ0qKChQT0+PPVNUVKTGxkZVVVWpqqpKjY2N8vl8g3qeAABg6IqL+QLi4qKuCl0QiUS0efNmrVu3Tvfdd58k6cUXX5TL5dIrr7yiJUuWKBQK6YUXXtDOnTs1e/ZsSdLLL7+stLQ0vf7665o7d66OHTumqqoq1dXVKTc3V5K0fft2eb1eHT9+XJMnT77kusLhsMLhsH2/vb29v08dAAAMETG/QvTuu+/K4/EoPT1df/u3f6v3339fktTU1KRAIKD8/Hx71ul0asaMGaqtrZUkNTQ06Ny5c1EzHo9HWVlZ9syBAwdkWZYdQ5I0ffp0WZZlz1xKWVmZ/RKbZVlKS0vr1/MGAABDR0yDKDc3Vy+99JJ+/etfa/v27QoEAsrLy9OHH36oQCAgSXK5XFGPcblc9r5AIKD4+HiNGTPmsjOpqam9njs1NdWeuZSSkhKFQiH7dvLkyS90rgAAYOiK6Utm8+fPt/87OztbXq9XX/7yl/Xiiy9q+vTpkiSHwxH1mEgk0mvbxS6eudT8Zx3H6XTK6XRe0XkAAIDhLeYvmf250aNHKzs7W++++679vqKLr+K0tLTYV43cbre6u7sVDAYvO3PmzJlez9Xa2trr6hMAADDTkAqicDisY8eO6YYbblB6errcbreqq6vt/d3d3aqpqVFeXp4kKScnRyNHjoyaaW5u1pEjR+wZr9erUCikQ4cO2TMHDx5UKBSyZwAAgNli+pLZ6tWrdc8992jChAlqaWnR97//fbW3t2vhwoVyOBzy+/0qLS1VRkaGMjIyVFpaqlGjRqmoqEiSZFmWiouLtWrVKo0dO1bJyclavXq1srOz7U+dZWZmat68eVq0aJG2bt0qSVq8eLEKCgo+9RNmAADALDENolOnTunrX/+6PvjgA40bN07Tp09XXV2dJk6cKElas2aNurq6tGzZMgWDQeXm5mrPnj1KTEy0j7Fp0ybFxcVpwYIF6urq0qxZs7Rjxw6NGDHCntm1a5dWrFhhfxqtsLBQFRUVg3uyAABgyHJEIpFIrBcxHLS3t8uyLIVCISUlJQ3Y8+R8+6UBOzYwXDU8/XCslwBgmLrS399D6j1EAAAAsUAQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIw3ZIKorKxMDodDfr/f3haJRLR+/Xp5PB4lJCRo5syZOnr0aNTjwuGwli9frpSUFI0ePVqFhYU6depU1EwwGJTP55NlWbIsSz6fT21tbYNwVgAAYDgYEkFUX1+vbdu2aerUqVHbN27cqPLyclVUVKi+vl5ut1tz5sxRR0eHPeP3+1VZWandu3dr//796uzsVEFBgXp6euyZoqIiNTY2qqqqSlVVVWpsbJTP5xu08wMAAENbzIOos7NTDz74oLZv364xY8bY2yORiDZv3qx169bpvvvuU1ZWll588UV99NFHeuWVVyRJoVBIL7zwgp555hnNnj1bt956q15++WUdPnxYr7/+uiTp2LFjqqqq0r/8y7/I6/XK6/Vq+/bt+o//+A8dP348JucMAACGlpgH0aOPPqq//uu/1uzZs6O2NzU1KRAIKD8/397mdDo1Y8YM1dbWSpIaGhp07ty5qBmPx6OsrCx75sCBA7IsS7m5ufbM9OnTZVmWPXMp4XBY7e3tUTcAAHB1iovlk+/evVsNDQ168803e+0LBAKSJJfLFbXd5XLpf/7nf+yZ+Pj4qCtLF2YuPD4QCCg1NbXX8VNTU+2ZSykrK9OTTz75+U4IAAAMSzG7QnTy5En9/d//vXbt2qVrr732U+ccDkfU/Ugk0mvbxS6eudT8Zx2npKREoVDIvp08efKyzwkAAIavmAVRQ0ODWlpalJOTo7i4OMXFxammpkb/9E//pLi4OPvK0MVXcVpaWux9brdb3d3dCgaDl505c+ZMr+dvbW3tdfXpzzmdTiUlJUXdAADA1SlmQTRr1iwdPnxYjY2N9m3atGl68MEH1djYqEmTJsntdqu6utp+THd3t2pqapSXlydJysnJ0ciRI6NmmpubdeTIEXvG6/UqFArp0KFD9szBgwcVCoXsGQAAYLaYvYcoMTFRWVlZUdtGjx6tsWPH2tv9fr9KS0uVkZGhjIwMlZaWatSoUSoqKpIkWZal4uJirVq1SmPHjlVycrJWr16t7Oxs+03amZmZmjdvnhYtWqStW7dKkhYvXqyCggJNnjx5EM8YAAAMVTF9U/VnWbNmjbq6urRs2TIFg0Hl5uZqz549SkxMtGc2bdqkuLg4LViwQF1dXZo1a5Z27NihESNG2DO7du3SihUr7E+jFRYWqqKiYtDPBwAADE2OSCQSifUihoP29nZZlqVQKDSg7yfK+fZLA3ZsYLhqePrhWC8BwDB1pb+/Y/49RAAAALFGEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4Q/qPuwLA1eTEhuxYLwEYciY8cTjWS5DEFSIAAACCCAAAgCACAADGI4gAAIDxCCIAAGC8PgXR3Xffrba2tl7b29vbdffdd3/RNQEAAAyqPgXR3r171d3d3Wv7xx9/rH379n3hRQEAAAymz/U9RL///e/t/3777bcVCATs+z09PaqqqtKNN97Yf6sDAAAYBJ8riG655RY5HA45HI5LvjSWkJCg5557rt8WBwAAMBg+VxA1NTUpEolo0qRJOnTokMaNG2fvi4+PV2pqqkaMGNHviwQAABhInyuIJk6cKEk6f/78gCwGAAAgFvr8t8z+8Ic/aO/evWppaekVSE888cQXXhgAAMBg6VMQbd++XX/3d3+nlJQUud1uORwOe5/D4SCIAADAsNKnIPr+97+vH/zgB3r88cf7ez0AAACDrk/fQxQMBnX//ff391oAAABiok9BdP/992vPnj39vRYAAICY6NNLZl/5ylf03e9+V3V1dcrOztbIkSOj9q9YsaJfFgcAADAY+hRE27Zt03XXXaeamhrV1NRE7XM4HAQRAAAYVvoURE1NTf29DgAAgJjp03uIAAAAriZ9ukL0yCOPXHb/j370oz4tBgAAIBb6FETBYDDq/rlz53TkyBG1tbVd8o++AgAADGV9CqLKyspe286fP69ly5Zp0qRJX3hRAAAAg6nf3kN0zTXX6LHHHtOmTZv665AAAACDol/fVP1f//Vf+uMf/9ifhwQAABhwfXrJbOXKlVH3I5GImpub9ctf/lILFy7sl4UBAAAMlj4F0VtvvRV1/5prrtG4ceP0zDPPfOYn0AAAAIaaPgXRb3/72/5eBwAAQMz0KYguaG1t1fHjx+VwOPQXf/EXGjduXH+tCwAAYND06U3VZ8+e1SOPPKIbbrhBd955p+644w55PB4VFxfro48+6u81AgAADKg+BdHKlStVU1Ojf//3f1dbW5va2tr085//XDU1NVq1alV/rxEAAGBA9ekls1dffVU//elPNXPmTHvbX/3VXykhIUELFizQli1b+mt9AAAAA65PV4g++ugjuVyuXttTU1N5yQwAAAw7fQoir9er733ve/r444/tbV1dXXryySfl9Xr7bXEAAACDoU9BtHnzZtXW1mr8+PGaNWuWZs+erbS0NL3xxht69tlnr/g4W7Zs0dSpU5WUlKSkpCR5vV796le/svdHIhGtX79eHo9HCQkJmjlzpo4ePRp1jHA4rOXLlyslJUWjR49WYWGhTp06FTUTDAbl8/lkWZYsy5LP51NbW1tfTh0AAFyF+hRE2dnZevfdd1VWVqZbbrlFU6dO1Q9/+EO99957uummm674OOPHj9cPf/hDvfnmm3rzzTd1991362/+5m/s6Nm4caPKy8tVUVGh+vp6ud1uzZkzRx0dHfYx/H6/KisrtXv3bu3fv1+dnZ0qKChQT0+PPVNUVKTGxkZVVVWpqqpKjY2N8vl8fTl1AABwFXJEIpHI531QWVmZXC5Xr2+l/tGPfqTW1lY9/vjjfV5QcnKynn76aT3yyCPyeDzy+/328cLhsFwul5566iktWbJEoVBI48aN086dO/XAAw9Ikk6fPq20tDS99tprmjt3ro4dO6YpU6aorq5Oubm5kqS6ujp5vV698847mjx58iXXEQ6HFQ6H7fvt7e1KS0tTKBRSUlJSn8/vs+R8+6UBOzYwXDU8/XCsl9AvTmzIjvUSgCFnwhOHB/T47e3tsizrM39/9+kK0datW/WXf/mXvbbfdNNN+ud//ue+HFI9PT3avXu3zp49K6/Xq6amJgUCAeXn59szTqdTM2bMUG1trSSpoaFB586di5rxeDzKysqyZw4cOCDLsuwYkqTp06fLsix75lLKysrsl9gsy1JaWlqfzgsAAAx9fQqiQCCgG264odf2cePGqbm5+XMd6/Dhw7ruuuvkdDq1dOlSVVZWasqUKQoEApLU69NsLpfL3hcIBBQfH68xY8ZcdiY1NbXX86amptozl1JSUqJQKGTfTp48+bnOCwAADB99+h6iC2+gTk9Pj9r+xhtvyOPxfK5jTZ48WY2NjWpra9Orr76qhQsXqqamxt7vcDii5iORSK9tF7t45lLzn3Ucp9Mpp9N5pacBAACGsT4F0Te/+U35/X6dO3dOd999tyTpP//zP7VmzZrP/U3V8fHx+spXviJJmjZtmurr6/Xss8/a7xu6+GpUS0uLfdXI7Xaru7tbwWAw6ipRS0uL8vLy7JkzZ870et7W1tZLfpcSAAAwT59eMluzZo2Ki4u1bNkyTZo0SZMmTdLy5cu1YsUKlZSUfKEFRSIRhcNhpaeny+12q7q62t7X3d2tmpoaO3ZycnI0cuTIqJnm5mYdOXLEnvF6vQqFQjp06JA9c/DgQYVCIXsGAACYrU9XiBwOh5566il997vf1bFjx5SQkKCMjIzP/RLT2rVrNX/+fKWlpamjo0O7d+/W3r17VVVVJYfDIb/fr9LSUmVkZCgjI0OlpaUaNWqUioqKJEmWZam4uFirVq3S2LFjlZycrNWrVys7O1uzZ8+WJGVmZmrevHlatGiRtm7dKklavHixCgoKPvUTZgAAwCx9CqILrrvuOt1+++19fvyZM2fk8/nU3Nwsy7I0depUVVVVac6cOZI+uRLV1dWlZcuWKRgMKjc3V3v27FFiYqJ9jE2bNikuLk4LFixQV1eXZs2apR07dmjEiBH2zK5du7RixQr702iFhYWqqKjo87oBAMDVpU/fQ2SiK/0egy+K7yECeuN7iICr17D+HiIAAICrCUEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMF5Mg6isrEy33367EhMTlZqaqnvvvVfHjx+PmolEIlq/fr08Ho8SEhI0c+ZMHT16NGomHA5r+fLlSklJ0ejRo1VYWKhTp05FzQSDQfl8PlmWJcuy5PP51NbWNtCnCAAAhoGYBlFNTY0effRR1dXVqbq6Wn/84x+Vn5+vs2fP2jMbN25UeXm5KioqVF9fL7fbrTlz5qijo8Oe8fv9qqys1O7du7V//351dnaqoKBAPT099kxRUZEaGxtVVVWlqqoqNTY2yufzDer5AgCAockRiUQisV7EBa2trUpNTVVNTY3uvPNORSIReTwe+f1+Pf7445I+uRrkcrn01FNPacmSJQqFQho3bpx27typBx54QJJ0+vRppaWl6bXXXtPcuXN17NgxTZkyRXV1dcrNzZUk1dXVyev16p133tHkyZM/c23t7e2yLEuhUEhJSUkD9m+Q8+2XBuzYwHDV8PTDsV5CvzixITvWSwCGnAlPHB7Q41/p7+8h9R6iUCgkSUpOTpYkNTU1KRAIKD8/355xOp2aMWOGamtrJUkNDQ06d+5c1IzH41FWVpY9c+DAAVmWZceQJE2fPl2WZdkzFwuHw2pvb4+6AQCAq9OQCaJIJKKVK1fqa1/7mrKysiRJgUBAkuRyuaJmXS6XvS8QCCg+Pl5jxoy57Exqamqv50xNTbVnLlZWVma/38iyLKWlpX2xEwQAAEPWkAmib33rW/r973+vf/3Xf+21z+FwRN2PRCK9tl3s4plLzV/uOCUlJQqFQvbt5MmTV3IaAABgGBoSQbR8+XL94he/0G9/+1uNHz/e3u52uyWp11WclpYW+6qR2+1Wd3e3gsHgZWfOnDnT63lbW1t7XX26wOl0KikpKeoGAACuTjENokgkom9961v6t3/7N/3mN79Renp61P709HS53W5VV1fb27q7u1VTU6O8vDxJUk5OjkaOHBk109zcrCNHjtgzXq9XoVBIhw4dsmcOHjyoUChkzwAAAHPFxfLJH330Ub3yyiv6+c9/rsTERPtKkGVZSkhIkMPhkN/vV2lpqTIyMpSRkaHS0lKNGjVKRUVF9mxxcbFWrVqlsWPHKjk5WatXr1Z2drZmz54tScrMzNS8efO0aNEibd26VZK0ePFiFRQUXNEnzAAAwNUtpkG0ZcsWSdLMmTOjtv/4xz/WN77xDUnSmjVr1NXVpWXLlikYDCo3N1d79uxRYmKiPb9p0ybFxcVpwYIF6urq0qxZs7Rjxw6NGDHCntm1a5dWrFhhfxqtsLBQFRUVA3uCAABgWBhS30M0lPE9REDs8D1EwNWL7yECAAAYIggiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYLyYBtHvfvc73XPPPfJ4PHI4HPrZz34WtT8SiWj9+vXyeDxKSEjQzJkzdfTo0aiZcDis5cuXKyUlRaNHj1ZhYaFOnToVNRMMBuXz+WRZlizLks/nU1tb2wCfHQAAGC5iGkRnz57VzTffrIqKikvu37hxo8rLy1VRUaH6+nq53W7NmTNHHR0d9ozf71dlZaV2796t/fv3q7OzUwUFBerp6bFnioqK1NjYqKqqKlVVVamxsVE+n2/Azw8AAAwPcbF88vnz52v+/PmX3BeJRLR582atW7dO9913nyTpxRdflMvl0iuvvKIlS5YoFArphRde0M6dOzV79mxJ0ssvv6y0tDS9/vrrmjt3ro4dO6aqqirV1dUpNzdXkrR9+3Z5vV4dP35ckydPvuTzh8NhhcNh+357e3t/njoAABhChux7iJqamhQIBJSfn29vczqdmjFjhmprayVJDQ0NOnfuXNSMx+NRVlaWPXPgwAFZlmXHkCRNnz5dlmXZM5dSVlZmv8RmWZbS0tL6+xQBAMAQMWSDKBAISJJcLlfUdpfLZe8LBAKKj4/XmDFjLjuTmpra6/ipqan2zKWUlJQoFArZt5MnT36h8wEAAENXTF8yuxIOhyPqfiQS6bXtYhfPXGr+s47jdDrldDo/52oBAMBwNGSvELndbknqdRWnpaXFvmrkdrvV3d2tYDB42ZkzZ870On5ra2uvq08AAMBMQzaI0tPT5Xa7VV1dbW/r7u5WTU2N8vLyJEk5OTkaOXJk1Exzc7OOHDliz3i9XoVCIR06dMieOXjwoEKhkD0DAADMFtOXzDo7O/Xee+/Z95uamtTY2Kjk5GRNmDBBfr9fpaWlysjIUEZGhkpLSzVq1CgVFRVJkizLUnFxsVatWqWxY8cqOTlZq1evVnZ2tv2ps8zMTM2bN0+LFi3S1q1bJUmLFy9WQUHBp37CDAAAmCWmQfTmm2/qrrvusu+vXLlSkrRw4ULt2LFDa9asUVdXl5YtW6ZgMKjc3Fzt2bNHiYmJ9mM2bdqkuLg4LViwQF1dXZo1a5Z27NihESNG2DO7du3SihUr7E+jFRYWfup3HwEAAPM4IpFIJNaLGA7a29tlWZZCoZCSkpIG7Hlyvv3SgB0bGK4ann441kvoFyc2ZMd6CcCQM+GJwwN6/Cv9/T1k30MEAAAwWAgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPGMCqLnn39e6enpuvbaa5WTk6N9+/bFekkAAGAIMCaIfvKTn8jv92vdunV66623dMcdd2j+/Pk6ceJErJcGAABizJggKi8vV3Fxsb75zW8qMzNTmzdvVlpamrZs2RLrpQEAgBiLi/UCBkN3d7caGhr0ne98J2p7fn6+amtrL/mYcDiscDhs3w+FQpKk9vb2gVuopJ5w14AeHxiOBvrnbrB0fNwT6yUAQ85A/3xfOH4kErnsnBFB9MEHH6inp0culytqu8vlUiAQuORjysrK9OSTT/banpaWNiBrBPDprOeWxnoJAAZKmTUoT9PR0SHL+vTnMiKILnA4HFH3I5FIr20XlJSUaOXKlfb98+fP6//+7/80duzYT30Mrh7t7e1KS0vTyZMnlZSUFOvlAOhH/HybJRKJqKOjQx6P57JzRgRRSkqKRowY0etqUEtLS6+rRhc4nU45nc6obddff/1ALRFDVFJSEv/DBK5S/Hyb43JXhi4w4k3V8fHxysnJUXV1ddT26upq5eXlxWhVAABgqDDiCpEkrVy5Uj6fT9OmTZPX69W2bdt04sQJLV3KexMAADCdMUH0wAMP6MMPP9SGDRvU3NysrKwsvfbaa5o4cWKsl4YhyOl06nvf+16vl00BDH/8fONSHJHP+hwaAADAVc6I9xABAABcDkEEAACMRxABAADjEUQAAMB4BBFwkeeff17p6em69tprlZOTo3379sV6SQD6we9+9zvdc8898ng8cjgc+tnPfhbrJWEIIYiAP/OTn/xEfr9f69at01tvvaU77rhD8+fP14kTJ2K9NABf0NmzZ3XzzTeroqIi1kvBEMTH7oE/k5ubq9tuu01btmyxt2VmZuree+9VWVlZDFcGoD85HA5VVlbq3nvvjfVSMERwhQj4k+7ubjU0NCg/Pz9qe35+vmpra2O0KgDAYCCIgD/54IMP1NPT0+sP/rpcrl5/GBgAcHUhiICLOByOqPuRSKTXNgDA1YUgAv4kJSVFI0aM6HU1qKWlpddVIwDA1YUgAv4kPj5eOTk5qq6ujtpeXV2tvLy8GK0KADAYjPlr98CVWLlypXw+n6ZNmyav16tt27bpxIkTWrp0aayXBuAL6uzs1HvvvWffb2pqUmNjo5KTkzVhwoQYrgxDAR+7By7y/PPPa+PGjWpublZWVpY2bdqkO++8M9bLAvAF7d27V3fddVev7QsXLtSOHTsGf0EYUggiAABgPN5DBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEwyocffqjx48fL4XCora0tat/hw4c1Y8YMJSQk6MYbb9SGDRv051/m39zcrKKiIk2ePFnXXHON/H7/4C4ewIAhiAAYpbi4WFOnTu21vb29XXPmzJHH41F9fb2ee+45/eM//qPKy8vtmXA4rHHjxmndunW6+eabB3PZAAYYQQTgqhCJRLRx40ZNmjRJCQkJuvnmm/XTn/40ambLli1qa2vT6tWrez1+165d+vjjj7Vjxw5lZWXpvvvu09q1a1VeXm5fJfrSl76kZ599Vg8//LAsyxqU8wIwOAgiAFeFf/iHf9CPf/xjbdmyRUePHtVjjz2mhx56SDU1NZKkt99+Wxs2bNBLL72ka67p/b++AwcOaMaMGXI6nfa2uXPn6vTp0/rv//7vwToNADESF+sFAMAXdfbsWZWXl+s3v/mNvF6vJGnSpEnav3+/tm7dqunTp+vrX/+6nn76aU2YMEHvv/9+r2MEAgF96Utfitrmcrnsfenp6QN+HgBihyACMOy9/fbb+vjjjzVnzpyo7d3d3br11ltVUlKizMxMPfTQQ5c9jsPhiLp/4aWyi7cDuPoQRACGvfPnz0uSfvnLX+rGG2+M2ud0OnXPPffo8OHD9nuKLoROSkqK1q1bpyeffFJut1uBQCDqsS0tLZL+/5UiAFcvggjAsDdlyhQ5nU6dOHFCM2bM6LX/1VdfVVdXl32/vr5ejzzyiPbt26cvf/nLkiSv16u1a9equ7tb8fHxkqQ9e/bI4/H0eikNwNWHIAIw7CUmJmr16tV67LHHdP78eX3ta19Te3u7amtrdd1112nhwoVR8x988IEkKTMzU9dff70kqaioSE8++aS+8Y1vaO3atXr33XdVWlqqJ554Iuols8bGRklSZ2enWltb1djYqPj4eE2ZMmVQzhXAwHBE/vxbxwBgmIpEInruuef0/PPP6/3339f111+v2267TWvXrtWdd94ZNbt3717dddddCgaDdhBJn3wx46OPPqpDhw5pzJgxWrp0aa8gutT7iSZOnMgn0YBhjiACAADG43uIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGO//AVKl47No7XF5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.countplot(data['e401'],)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "material-sending",
   "metadata": {
    "papermill": {
     "duration": 0.039406,
     "end_time": "2021-03-24T14:13:53.435306",
     "exception": false,
     "start_time": "2021-03-24T14:13:53.395900",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Eligibility is highly associated with financial wealth:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a7741bb6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABCAAAAHpCAYAAABTM02+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABeYElEQVR4nO3deXyU5b3///c9WSZhC3sSBMIiO0ojVAgWkbJoqBxxqbQqgmutnlOVcvo17qCV4xEQ3FBPKSk/LWKLoFZcUIFoSVUocSmKoGAoJmKQJCSQbe7798cwA5N9mXsmk/v1fDzmMc4999xzzZzhXPTN5/pchmVZlgAAAAAAAGzkCvcAAAAAAABA20cAAQAAAAAAbEcAAQAAAAAAbEcAAQAAAAAAbEcAAQAAAAAAbEcAAQAAAAAAbEcAAQAAAAAAbEcAAQAAAAAAbEcAAQAAAAAAbEcAAQAAAAAAbEcAAaDJDh8+rN69e8swDBUWFgY89+mnn2rixImKj4/XaaedpoULF8qyLP/zeXl5uuKKKzRkyBC5XC7ddtttIRnzkSNHNHv2bCUkJCghIUGzZ8+uMXYAAIIlEufK3//+9xo/frzatWunzp07h+Q9ATgLAQSAJrvuuut05pln1jheXFysqVOnqlevXvroo4/0+OOPa/HixVq6dKn/nPLycvXo0UN33XWXRo0aFbIxX3HFFcrJydEbb7yhN954Qzk5OZo9e3bI3h8A4CyROFdWVFTo5z//uX7961+H7D0BOAsBBABJkmVZ+t///V8NGDBA8fHxGjVqlP7617/WOG/FihUqLCzU/Pnzazz3/PPPq6ysTJmZmRo5cqQuueQS3XnnnVq6dKn/X3b69eun5cuX6+qrr1ZCQoLtn0uSPv/8c73xxhv6wx/+oLS0NKWlpen//u//9Le//U27d+8OyRgAAJGvLc+VkrRgwQLdfvvtOuOMM0L2ngCchQACgCTp7rvv1qpVq7RixQr961//0u23366rrrpKW7du9Z+za9cuLVy4UKtXr5bLVfP/fWRnZ2vixIlyu93+Y+eff76+/fZb7d+/v0Xj69ChQ7239PT0Ol+bnZ2thIQEjR071n9s3LhxSkhI0LZt21o0LgCAc7TluRIAQiE63ANo7bKysvTII49ox44dysvL0/r16zVz5kzb3u/+++/XggULAo4lJiYqPz/ftvcESktLtXTpUr377rtKS0uTJA0YMEDvv/++nnnmGU2cOFHl5eX65S9/qUceeUR9+/bV119/XeM6+fn56tevX8CxxMRE/3P9+/dv9hhzcnLqfT4+Pr7O5/Lz89WzZ88ax3v27MmfLQBAo7T1uRIAQoEAogGlpaUaNWqUrrnmGl166aUhec8RI0bo7bff9j+OiooKyfvCuXbt2qWysjJNnTo14HhFRYVSU1MlSRkZGRo2bJiuuuqqeq9lGEbAY185afXjTXX66ae36PW1vb9lWS0eFwDAGZwwVwKA3QggGpCenl5vuVpFRYXuvvtuPf/88yosLNTIkSP18MMP67zzzmv2e0ZHRyspKanZrweayjRNSdJrr72m0047LeA5X4nou+++q08//dS/1tX3l6Xu3bvrrrvu0oIFC5SUlFSjouDQoUOSTv7rTnN16NCh3ucnTJig119/vdbnkpKS9N1339U4/v3337d4XAAAZ2jrcyUAhAIBRAtdc8012r9/v1544QX16tVL69ev1wUXXKBPP/1UgwYNatY19+zZo169esntdmvs2LF66KGHNGDAgCCPHDhp+PDhcrvdys3N1cSJE2s9Z926dTp+/Lj/8UcffaRrr71W7733ngYOHChJSktL05133qmKigrFxsZKkt566y316tWrRrlpU7WkrDQtLU1FRUX68MMPdfbZZ0uSPvjgAxUVFWn8+PEtGhcAwBna+lwJAKFAANECX331ldasWaN///vf6tWrlyRp/vz5euONN7Rq1So99NBDTb7m2LFjtXr1ag0ePFjfffedHnzwQY0fP17/+te/1K1bt2B/BECS1LFjR82fP1+33367TNPUT37yExUXF2vbtm3q0KGD5syZ4/+Lk09BQYEkadiwYf69wq+44gotWLBAc+fO1Z133qk9e/booYce0r333htQVur7C1JJSYm+//575eTkKDY2VsOHD69zjC0pKx02bJguuOAC3XDDDXrmmWckSTfeeKMuvPBCDRkypNnXBQA4R1ufKyUpNzdXP/zwg3Jzc+XxePxjOP300xusrgCARrHQaJKs9evX+x+/+OKLliSrffv2Abfo6Gjr8ssvtyzLsvbt22dJqvd2yy231PmeJSUlVmJiorVkyRK7Px4czjRNa/ny5daQIUOsmJgYq0ePHtb5559vbd26tdbzN2/ebEmyjhw5EnD8k08+sSZMmGC53W4rKSnJuv/++y3TNAPOqe3PQUpKik2fzOvw4cPWlVdeaXXs2NHq2LGjdeWVV9YYOwAA9Wnrc+WcOXNqfd/Nmzfb+r4AnMOwrBOL09AgwzACdsFYu3atrrzySv3rX/+q0SiyQ4cOSkpKUmVlpb766qt6r9ulS5d61/xNnTpVp59+ulasWNHizwAAAAAAQDiwBKMFUlNT5fF4dOjQIU2YMKHWc2JiYjR06NBmv0d5ebk+//zzOq8PAAAAAEAkIIBoQElJifbu3et/vG/fPuXk5Khr164aPHiwrrzySl199dVasmSJUlNTVVBQoHfffVdnnHGGpk+f3uT3mz9/vmbMmKG+ffvq0KFDevDBB1VcXKw5c+YE82MBAAAAABBSLMFowJYtWzRp0qQax+fMmaPMzExVVlbqwQcf1OrVq3Xw4EF169ZNaWlpWrBggc4444wmv98vfvELZWVlqaCgQD169NC4ceP0wAMP1NtwCAAAAACA1o4AAgAAAAAA2M4V7gEAAAAAAIC2jwCiFpZlqbi4WBSHAABQN+ZLAADQFAQQtTh69KgSEhJ09OjRcA8FAIBWi/kSAAA0BQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEAAAAAACwHQEEHO/zvGLtzj8a7mEAAAAAQJtGAAHH+983vtCjm74M9zAAAAAAoE2LDvcAgHArKa9SlMsI9zAAAAAAoE2jAgKOd6zCI49phXsYAAAAANCmEUDA8Y4TQAAAAACA7Qgg4HjHKz0ifwAAAAAAexFAwPGogAAAAAAA+xFAwPHKKj0yLQIIAAAAALATAQQczTQtlVWZBBAAAAAAYDMCCDhaWZVHkliCAQAAAAA2I4CAox2vOBFAUAEBAAAAALYigICjHa/0BhCmGeaBAAAAAEAbRwABR/NVQNADAgAAAADsRQABR/NVQNADAgAAAADsRQABR/P3gCCAAAAAAABbEUDA0Y6dqIBgBQYAAAAA2IsAAo5Wxi4YAAAAABASBBBwNHpAAAAAAEBoEEDA0Y5V+JZgEEAAAAAAgJ0IIOBoZZUswQAAAACAUCCAgKP5dsEwzTAPBAAAAADauLAGEFlZWZoxY4Z69eolwzC0YcOGes+fO3euDMOocRsxYoT/nMzMzFrPKSsrs/nTIBL5dsEwqYAAAAAAAFuFNYAoLS3VqFGj9MQTTzTq/OXLlysvL89/O3DggLp27aqf//znAed16tQp4Ly8vDzFxcXZ8REQ4Y6zCwYAAAAAhER0ON88PT1d6enpjT4/ISFBCQkJ/scbNmzQkSNHdM011wScZxiGkpKSgjZOtF2+HhAmu2AAAAAAgK0iugfEypUrNWXKFKWkpAQcLykpUUpKinr37q0LL7xQO3furPc65eXlKi4uDrjBGXy7YJA/AEDDmC8BAEBLRGwAkZeXp9dff13XX399wPGhQ4cqMzNTr7zyitasWaO4uDidc8452rNnT53XWrRokb+6IiEhQX369LF7+GgljlMBAQCNxnwJAABawrCs1rH43TAMrV+/XjNnzmzU+YsWLdKSJUv07bffKjY2ts7zTNPUWWedpXPPPVePPfZYreeUl5ervLzc/7i4uFh9+vRRUVGROnXq1KTPgchy5R/+ob/vPax2sVHatfCCcA8HAFo15ksAANASYe0B0VyWZemPf/yjZs+eXW/4IEkul0s//vGP662AcLvdcrvdwR4mIoB/G87WkcMBQKvGfAkAAFoiIpdgbN26VXv37tV1113X4LmWZSknJ0fJyckhGBkizcklGGEeCAAAAAC0cWGtgCgpKdHevXv9j/ft26ecnBx17dpVffv2VUZGhg4ePKjVq1cHvG7lypUaO3asRo4cWeOaCxYs0Lhx4zRo0CAVFxfrscceU05Ojp588knbPw8iz7EKj9zRLnnoAQEAAAAAtgprALF9+3ZNmjTJ/3jevHmSpDlz5igzM1N5eXnKzc0NeE1RUZHWrVun5cuX13rNwsJC3XjjjcrPz1dCQoJSU1OVlZWls88+274PgohVVulRfEyUissqwz0UAAAAAGjTWk0TytakuLhYCQkJNNVygDPvf1PtYqOVX1ymfYumyzCMcA8JACIG8yUAAGiKiOwBAQTL8UqP4mOjJEmswgAAAAAA+xBAwLGqPKYqPZbiYrx/DOgDAQAAAAD2IYCAY5VVebe+iI/xVUAQQAAAAACAXQgg4FjHKqokSXEEEAAAAABgOwIIOFZZRWAFBEswAAAAAMA+BBBwrOOVHkmnVkCEczQAAAAA0LYRQMCxaizBIIEAAAAAANsQQMCxfBUQ/iUY9IAAAAAAANsQQMCxynxLMGK9fwxoQgkAAAAA9iGAgGMdqwisgDDNcI4GAAAAANo2Agg41vGKwCaULMEAAAAAAPsQQMCxKjymDEkxUSeWYNCEEgAAAABsQwABxzJNSy6XIZdx4jEVEAAAAABgGwIIOJZpSS5DchneBMJDBQQAAAAA2IYAAo5lWpZchiHDXwER3vEAAAAAQFtGAAHHMi3JOKUCgiUYAAAAAGAfAgg4lmlaMmSwBAMAAAAAQoAAAo5lWpZcLp2yBIMAAgAAAADsQgABx/I2oTxZAWGaYR4QAAAAALRhBBBwLNOyZEj+bTg9VEAAAAAAgG0IIOBYpmkFVEDQAwIAAAAA7EMAAceqvgTDogICAAAAAGxDAAHHMi3rxDac3sdUQAAAAACAfQgg4FjWiQDC8C3BoAICAAAAAGxDAAHH8li+HhDex+QPAAAAAGAfAgg4lr8HhIsmlAAAAABgNwIIOJavB8SJAgiWYAAAAACAjQgg4FiWJW8TShe7YAAAAACA3Qgg4Fge05Khk9tweswwDwgAAAAA2jACCDiWaVlyubxVEBI9IAAAAADATgQQcCzLklynVECwBAMAAAAA7EMAAcfyNaH0bcNJE0oAAAAAsA8BBBzLY1oyjFN7QBBAAAAAAIBdCCDgWKblrX44uQQjzAMCAAAAgDaMAAKOZVmWXIZxcgkGFRAAAAAAYBsCCDiWdwnGyQoIekAAAAAAgH0IIOBYpiUZMvzbcJpUQAAAAACAbQgg4FiWZcnlkgzDkCFvIAEAAAAAsAcBBBzLtCx5owfJ5TJYggEAAAAANiKAgGN5LPmXXxgGSzAAAAAAwE4EEHAs07L8AUSUYcikAgIAAAAAbBPWACIrK0szZsxQr169ZBiGNmzYUO/5W7Zs8a7Xr3b74osvAs5bt26dhg8fLrfbreHDh2v9+vU2fgpEKt82nJJ3Jwy24QQAAAAA+4Q1gCgtLdWoUaP0xBNPNOl1u3fvVl5env82aNAg/3PZ2dmaNWuWZs+erY8//lizZ8/W5Zdfrg8++CDYw0eEM81qSzCogAAAAAAA20SH883T09OVnp7e5Nf17NlTnTt3rvW5ZcuWaerUqcrIyJAkZWRkaOvWrVq2bJnWrFnTkuGijfFYllw6WQFBAQQAAAAA2Ccie0CkpqYqOTlZkydP1ubNmwOey87O1rRp0wKOnX/++dq2bVud1ysvL1dxcXHADW2fdUoPCJchlmAAQAOYLwEAQEtEVACRnJysZ599VuvWrdNLL72kIUOGaPLkycrKyvKfk5+fr8TExIDXJSYmKj8/v87rLlq0SAkJCf5bnz59bPsMaD1MM7AHBLtgAED9mC8BAEBLhHUJRlMNGTJEQ4YM8T9OS0vTgQMHtHjxYp177rn+44bvn7VP8P5Ld+CxU2VkZGjevHn+x8XFxfylygHMU7bhZAkGADSM+RIAALRERAUQtRk3bpyee+45/+OkpKQa1Q6HDh2qURVxKrfbLbfbbdsY0Tp5TgmmXC7vYwBA3ZgvAQBAS0TUEoza7Ny5U8nJyf7HaWlp2rRpU8A5b731lsaPHx/qoaGVMy1LvroYgyUYAAAAAGCrsFZAlJSUaO/evf7H+/btU05Ojrp27aq+ffsqIyNDBw8e1OrVqyV5d7jo16+fRowYoYqKCj333HNat26d1q1b57/GrbfeqnPPPVcPP/ywLrroIr388st6++239f7774f886F1M03J5fL1gGAbTgAAAACwU1gDiO3bt2vSpEn+x751pXPmzFFmZqby8vKUm5vrf76iokLz58/XwYMHFR8frxEjRui1117T9OnT/eeMHz9eL7zwgu6++27dc889GjhwoNauXauxY8eG7oMhIpiWpehTmlCyBAMAAAAA7GNYFv+rq7ri4mIlJCSoqKhInTp1CvdwYJOfP71N7ugo3TLpdP32Lzn62RnJuutnw8M9LACIGMyXAACgKSK+BwTQXOyCAQAAAAChQwABxzItS65Tl2CQQAAAAACAbQgg4FimaelED0oZNKEEAAAAAFsRQMCxPJZ3+02JCggAAAAAsBsBBBzLsiydKIA4sQ1nWIcDAAAAAG0aAQQcy7QsuVwnKyBMEggAAAAAsA0BBBzLNOWvgDAMyUMPCAAAAACwDQEEHKv6Lhg0oQQAAAAA+xBAwLFMy9KJ/IElGAAAAABgMwIIOJZpyV8B4V2CEeYBAQAAAEAbRgABx6pRAcESDAAAAACwDQEEHMs0T+0BIZZgAAAAAICNCCDgWKYlfwWEYRjyEEAAAAAAgG0IIOBYgbtgeAMJAAAAAIA9CCDgWNYpFRD0gAAAAAAAexFAwLE8liVDp+yCQQkEAAAAANiGAAKO5W1C6f1vKiAAAAAAwF4EEHAsK6AHBE0oAQAAAMBOBBBwrMBdMFiCAQAAAAB2IoCAY3ksS8YpFRAswQAAAAAA+xBAwLEsS6f0gGAbTgAAAACwEwEEHMs8pQeEQQ8IAAAAALAVAQQcy7Qsfw8ImlACAAAAgL0IIOBYpikZ8vWAED0gAAAAAMBGBBBwLNOy5DrxJ8CgCSUAAAAA2IoAAo5lWYEVECzBAAAAAAD7EEDAsbxNKL3/7XIZogACAAAAAOxDAAFHsixLlrxLLyTvHwQqIAAAAADAPgQQcCRf1nBqBYSHEggAAAAAsA0BBBzJ13DSXwFBE0oAAAAAsBUBBBzJFzb4KyAM77acAAAAAAB7EEDAkXxhg68CwjBYggEAAAAAdiKAgCPVVgFhEUAAAAAAgG0IIOBI/h4QOtkDgl0wAAAAAMA+BBBwpOq7YBiGIZMAAgAAAABsQwABR/KFDSd3wTgZSgAAAAAAgo8AAo5UswcETSgBAAAAwE4EEHAkX7XDyV0wxBIMAAAAALARAQQcybfjhUEFBAAAAACEBAEEHMlTyxIM8gcAAAAAsA8BBBzp5C4YJ5tQsg0nAAAAANgnrAFEVlaWZsyYoV69eskwDG3YsKHe81966SVNnTpVPXr0UKdOnZSWlqY333wz4JzMzEwZhlHjVlZWZuMnQaSpvguGYRj+xpQAAAAAgOALawBRWlqqUaNG6YknnmjU+VlZWZo6dao2btyoHTt2aNKkSZoxY4Z27twZcF6nTp2Ul5cXcIuLi7PjIyBC+bKGEysw5HKxDScAAAAA2Ck6nG+enp6u9PT0Rp+/bNmygMcPPfSQXn75Zb366qtKTU31HzcMQ0lJScEaJtqgGj0gTkQRpmnJ5TsIAAAAAAiaiO4BYZqmjh49qq5duwYcLykpUUpKinr37q0LL7ywRoVEdeXl5SouLg64oW0z/QHEiR4QJ/4ksBMGANSN+RIAALRERAcQS5YsUWlpqS6//HL/saFDhyozM1OvvPKK1qxZo7i4OJ1zzjnas2dPnddZtGiREhIS/Lc+ffqEYvgIo5PbcPqaUJ6ogCCAAIA6MV8CAICWiNgAYs2aNbr//vu1du1a9ezZ03983LhxuuqqqzRq1ChNmDBBL774ogYPHqzHH3+8zmtlZGSoqKjIfztw4EAoPgLCyNfv4UTu4A8iTDNMAwKACMB8CQAAWiKsPSCaa+3atbruuuv0l7/8RVOmTKn3XJfLpR//+Mf1VkC43W653e5gDxOtmG/LTX8PiBP3LMEAgLoxXwIAgJaIuAqINWvWaO7cufrzn/+sn/3sZw2eb1mWcnJylJycHILRIVKYLMEAAAAAgJAKawVESUmJ9u7d63+8b98+5eTkqGvXrurbt68yMjJ08OBBrV69WpI3fLj66qu1fPlyjRs3Tvn5+ZKk+Ph4JSQkSJIWLFigcePGadCgQSouLtZjjz2mnJwcPfnkk6H/gGi1fDmDL3jwLcUw2YsTAAAAAGwR1gqI7du3KzU11b+F5rx585Samqp7771XkpSXl6fc3Fz/+c8884yqqqp0yy23KDk52X+79dZb/ecUFhbqxhtv1LBhwzRt2jQdPHhQWVlZOvvss0P74dCqnayA8D72BREeAggAAAAAsEVYKyDOO+88/24EtcnMzAx4vGXLlgav+eijj+rRRx9t4cjQ1vmChhP5wylLMMI0IAAAAABo4yKuBwQQDGa1JRi+JpT0gAAAAAAAexBAwJF8lTcne0CwBAMAAAAA7EQAAUfy5QxGtW04qYAAAAAAAHsQQMCR/D0gqjWhNM1wjQgAAAAA2jYCCDhS9SUYvgoIDxUQAAAAAGALAgg40skmlL57ekAAAAAAgJ0IIOBIvl4PvuaTrhNJRH3bwgIAAAAAmo8AAo7kW2pxogDCf88SDAAAAACwBwEEHMnfA8IVWAHBEgwAAAAAsAcBBBzJt9uFr/LB1wOCAggAAAAAsAcBBBzJrLYLhm87TiogAAAAAMAeBBBwpJNNKL2P/btgUAIBAAAAALYggIAjndyG80QPiBNBBLtgAAAAAIA9CCDgSNUrIHzbcXrMcI0IAAAAANo2Agg4Ul0VEPSAAAAAAAB7EEDAkUyz9h4QLMEAAAAAAHsQQMCRqu+C4a+AIIAAAAAAAFsQQMCRfCstauyCwRIMAAAAALAFAQQcyV8BIW/wYPiXYIRtSAAAAADQphFAwJFq9oDw3lMBAQAAAAD2IICAI5mWZOhk5YN/CQYlEAAAAABgCwIIOJJpWf7qB+lkAGFSAQEAAAAAtiCAgCNZluUPHaSTSzHIHwAAAADAHgQQcCSPWXsFBEswAAAAAMAeBBBwJNNSQAWE68SfBJZgAAAAAIA9CCDgSGa1JRj+HhBUQAAAAACALZoVQOzbty/Y4wBCyrIUsATDYBtOAAAAALBVswKI008/XZMmTdJzzz2nsrKyYI8JsJ2nrl0wqIAAAAAAAFs0K4D4+OOPlZqaqt/+9rdKSkrSr371K3344YfBHhtgG9OyFFVtCYYhyWOGb0wAAAAA0JY1K4AYOXKkli5dqoMHD2rVqlXKz8/XT37yE40YMUJLly7V999/H+xxAkHlXYJhBBxzuQx5TBIIAAAAALBDi5pQRkdH6+KLL9aLL76ohx9+WF999ZXmz5+v3r176+qrr1ZeXl6wxgkElVltG05Jchn0gAAAAAAAu7QogNi+fbtuvvlmJScna+nSpZo/f76++uorvfvuuzp48KAuuuiiYI0TCCqPZcmlahUQhiEP+QMAAAAA2CK6OS9aunSpVq1apd27d2v69OlavXq1pk+fLpfLm2f0799fzzzzjIYOHRrUwQLBYlqSq1r8FuUyZFIBAQAAAAC2aFYAsWLFCl177bW65pprlJSUVOs5ffv21cqVK1s0OMAulmXV7AFhGKoigAAAAAAAWzQrgNi0aZP69u3rr3jwsSxLBw4cUN++fRUbG6s5c+YEZZBAsJmWVW0BhrcHBNtwAgAAAIA9mtUDYuDAgSooKKhx/IcfflD//v1bPCjAbh7TW/FwKu8uGAQQAAAAACLP4cOH1bt3bxmGocLCwoDnPv30U02cOFHx8fE67bTTtHDhQlmn/ONrXl6errjiCg0ZMkQul0u33XabLWNsVgBh1fGvxCUlJYqLi2vRgIBQsCxLrmolEFEGAQQAAACAyHTdddfpzDPPrHG8uLhYU6dOVa9evfTRRx/p8ccf1+LFi7V06VL/OeXl5erRo4fuuusujRo1yrYxNmkJxrx58yRJhmHo3nvvVbt27fzPeTweffDBB/rRj34U1AECdjBr6wFBBQQAAACAVsayLD3yyCN6+umnlZeXp8GDB+uee+7RZZdd5j9nxYoVKiws1L333qvXX3894PXPP/+8ysrKlJmZKbfbrZEjR+rLL7/U0qVLNW/ePBmGoX79+mn58uWSpD/+8Y+2fZYmBRA7d+6U5P0CPv30U8XGxvqfi42N1ahRozR//vzgjhCwgWlJ1fIHuQzv9pwAAAAA0Frcfffdeumll7RixQoNGjRIWVlZuuqqq9SjRw9NnDhRu3bt0sKFC/XBBx/o66+/rvH67OxsTZw4UW6323/s/PPPV0ZGhvbv3x/SNgpNCiA2b94sSbrmmmu0fPlyderUyZZBAXbzmFYtAQTbcAIAAABoPUpLS7V06VK9++67SktLkyQNGDBA77//vp555hmNGzdOv/zlL/XII4+ob9++tQYQ+fn56tevX8CxxMRE/3OtNoDwWbVqVbDHAYSUtwdEzW04WYIBAAAAoLXYtWuXysrKNHXq1IDjFRUVSk1NVUZGhoYNG6arrrqq3utUX37u6+tY/bjdGt2E8pJLLlFxcbH/v+u7NVZWVpZmzJihXr16yTAMbdiwocHXbN26VaNHj1ZcXJwGDBigp59+usY569at0/Dhw+V2uzV8+HCtX7++0WOCM5hWbbtgSFUEEAAAAABaCdM0JUmvvfaacnJy/Lddu3bpr3/9q95991395S9/UXR0tKKjozV58mRJUvfu3XXfffdJkpKSkpSfnx9w3UOHDkk6WQkRKo2ugEhISPCnIwkJCUF589LSUo0aNUrXXHONLr300gbP37dvn6ZPn64bbrhBzz33nP7+97/r5ptvVo8ePfyvz87O1qxZs/TAAw/o4osv1vr163X55Zfr/fff19ixY4MybkQ+07JUPetzGYZMekAAAAAAaCV8/7Cem5uriRMn1nh+3bp1On78uP/xRx99pGuvvVbvvfeeBg4cKElKS0vTnXfeqYqKCn8fx7feeku9evWqsTTDbo0OIE5ddhGsJRjp6elKT09v9PlPP/20+vbtq2XLlkmShg0bpu3bt2vx4sX+AGLZsmWaOnWqMjIyJEkZGRnaunWrli1bpjVr1gRl3Ih83l0wAo+xBAMAAABAa9KxY0fNnz9ft99+u0zT1E9+8hMVFxdr27Zt6tChg+bMmRNwfkFBgSTv/1bu3LmzJOmKK67QggULNHfuXN15553as2ePHnroId17770BSzBycnIkSSUlJfr++++Vk5Oj2NhYDR8+PGifp1k9II4fPy7LsvzbcH7zzTdav369hg8frmnTpgVtcNVlZ2fXuP7555+vlStXqrKyUjExMcrOztbtt99e4xxfaFGb8vJylZeX+x/7lpqg7TLNmkswolyiAgIA6sF8CQBA6D3wwAPq2bOnFi1apK+//lqdO3fWWWedpTvvvLNRr09ISNCmTZt0yy23aMyYMerSpYvmzZunefPmBZyXmprq/+8dO3boz3/+s1JSUrR///6gfZZmBRAXXXSRLrnkEt10000qLCzU2WefrdjYWBUUFGjp0qX69a9/HbQBnio/P7/GGpXExERVVVWpoKBAycnJdZ5Tfc3LqRYtWqQFCxbYMma0TnVVQFR5CCAAoC7MlwAAhJ5hGPrNb36j3/zmNw2ee9555/kbTJ7qjDPOUFZWVr2vre11wdboJpSn+uc//6kJEyZIkv76178qKSlJ33zzjVavXq3HHnssqAOsrjHdO2s7p77unhkZGSoqKvLfDhw4EMQRozUyrZq/E8Mw5KECAgDqxHwJAABaolkVEMeOHVPHjh0leZtXXHLJJXK5XBo3bpy++eaboA7wVHV174yOjla3bt3qPae+7p5ut1tutzv4A0arZVqWXDUqICSTHhAAUCfmSwAA0BLNqoA4/fTTtWHDBh04cEBvvvmmvy/DoUOH1KlTp6AO8FRpaWnatGlTwLG33npLY8aMUUxMTL3njB8/3rZxIfJ4A4hqPSAMQ6zAAAAAAAB7NCuAuPfeezV//nz169dPY8eOVVpamiTv/9A/tXFFQ0pKSvz7mErebTZzcnKUm5sryVvqefXVV/vPv+mmm/TNN99o3rx5+vzzz/XHP/5RK1eu1Pz58/3n3HrrrXrrrbf08MMP64svvtDDDz+st99+W7fddltzPiraKNNSrdtwek7sswsAAAAACK5mLcG47LLL9JOf/ER5eXkaNWqU//jkyZN18cUXN/o627dv16RJk/yPfV0458yZo8zMTOXl5fnDCEnq37+/Nm7cqNtvv11PPvmkevXqpccee8y/BackjR8/Xi+88ILuvvtu3XPPPRo4cKDWrl2rsWPHNuejoo0ya+kLYrjENpwAAAAAYBPDCkWrywhTXFyshIQEFRUV2bqkBOHz6+d26GDhcWWkD/Mf+983vlC3Dm79Yc6YMI4MACIH8yUAAGiKZlVAlJaW6n/+53/0zjvv6NChQzKrla1//fXXQRkcYJdae0C4jBq/ZQAAAABAcDQrgLj++uu1detWzZ49W8nJyfVucQm0RnX1gKhiCQYAAAAA2KJZAcTrr7+u1157Teecc06wxwOEhFVbDwhD8rAiCQAAAEAIHSw8riOlFSF5ry7tY3Va5/iQvFdtmhVAdOnSRV27dg32WICQ8ZiWXNVKIFwuQx5WYAAAAAAIkYOFxzV5yRaVVYbmf4jExbj0zm/PC1sI0awA4oEHHtC9996rP/3pT2rXrl2wxwTYzrRUsweEYajc4wnTiAAAAAA4zZHSCpVVmrpl0um2hwIHC4/ryc17daS0osnv9dRTT+mRRx5RXl6eRowYoWXLlmnChAlNHkOzAoglS5boq6++UmJiovr166eYmJiA5//5z38257JAyHi34Qw85jLYhhMAAABA6J3WOV79u7cP9zBqtXbtWt1222166qmndM455+iZZ55Renq6du3apb59+zbpWs0KIGbOnNmclwGthsesLYAw6AEBAAAAAKdYunSprrvuOl1//fWSpGXLlunNN9/UihUrtGjRoiZdq1kBxH333declwGthmmpRhNKekAAAAAAwEkVFRXasWOH7rjjjoDj06ZN07Zt25p8PVdzB1JYWKg//OEPysjI0A8//CDJu/Ti4MGDzb0kEDKmZdXoAeEyDJkswQAAAAAASVJBQYE8Ho8SExMDjicmJio/P7/J12tWBcQnn3yiKVOmKCEhQfv379cNN9ygrl27av369frmm2+0evXq5lwWCBnLsmqkby5DqjIpgQAAAACAU1WvHrcsq8axxmhWBcS8efM0d+5c7dmzR3Fxcf7j6enpysrKas4lgZCqtQeEyxD5AwAAAAB4de/eXVFRUTWqHQ4dOlSjKqIxmhVAfPTRR/rVr35V4/hpp53WrDIMINRq7QFBE0oAAAAA8IuNjdXo0aO1adOmgOObNm3S+PHjm3y9Zi3BiIuLU3FxcY3ju3fvVo8ePZpzSSCkausBEcU2nAAAAADC4GDh8Vb7HvPmzdPs2bM1ZswYpaWl6dlnn1Vubq5uuummJl+rWQHERRddpIULF+rFF1+U5P2X5NzcXN1xxx269NJLm3NJIKRMy9vz4VTeXTAIIAAAAACERpf2sYqLcenJzXtD8n5xMS51aR/bpNfMmjVLhw8f1sKFC5WXl6eRI0dq48aNSklJafL7NyuAWLx4saZPn66ePXvq+PHjmjhxovLz85WWlqbf//73zbkkEFJmbT0gWIIBAAAAIIRO6xyvd357no6UVoTk/bq0j9VpneOb/Lqbb75ZN998c4vfv1kBRKdOnfT+++9r8+bN2rFjh0zT1FlnnaUpU6a0eEBAKJi1dG11GWIbTgAAAAAhdVrn+GaFApGoyQGEaZrKzMzUSy+9pP3798swDPXv319JSUnN3ooDCLXaekC4DJZgAAAAAIBdmrQLhmVZ+o//+A9df/31OnjwoM444wyNGDFC33zzjebOnauLL77YrnECQWVZqnUbTgIIAAAAALBHkyogMjMzlZWVpXfeeUeTJk0KeO7dd9/VzJkztXr1al199dVBHSQQbB7TqpG+0QMCAAAAAOzTpAqINWvW6M4776wRPkjST3/6U91xxx16/vnngzY4wC70gAAAAACA0GpSAPHJJ5/oggsuqPP59PR0ffzxxy0eFGC32rbhjHJRAQEAAAAAdmlSAPHDDz8oMTGxzucTExN15MiRFg8KsFttDVNpQgkAAAAA9mlSAOHxeBQdXXfbiKioKFVVVbV4UIDdPJZVowLCuwQjPOMBAAAAgLauSU0oLcvS3Llz5Xa7a32+vLw8KIMC7GZaqr0CgiUYAAAAAEKp8IB07HBo3qtdN6lzn9C8Vy2aFEDMmTOnwXPYAQORwKqtAuLEAdO0/P8NAAAAALYpPCA9+WOp8nho3i8mXrrlo7CFEE0KIFatWmXXOICQqqsCQpKqTEuxBBAAAAAA7HbssDd8mPBbKcHmUKDogPTeEu97NiGAyMrK0iOPPKIdO3YoLy9P69ev18yZM5s1hCYFEEBb4TEtVY8YfJmDyTIMAAAAAKGU0Efqdnq4R1Gr0tJSjRo1Stdcc40uvfTSFl2LAAKOZFk1l1n4KiDYCQMAAAAAvNLT05Wenh6UazVpFwygrTCtmj/+qBOBBI0oAQAAACD4CCDgSKZl1dkDwuMhgAAAAACAYCOAgCOZpiWj+i4YJx5TAQEAAAAAwUcAAUcyrZMVDz6+xyY9IAAAAAAg6Agg4EiWZan6TpsuekAAAAAAgG3YBQOOZFqqpQeE976KHhAAAAAAQqnoQKt9j5KSEu3du9f/eN++fcrJyVHXrl3Vt2/fJl2LAAKO5G1CGXjMvwSDCggAAAAAodCumxQTL723JDTvFxPvfc8m2L59uyZNmuR/PG/ePEnSnDlzlJmZ2aRrEUDAkUzLqtkDwrcEgx4QAAAAAEKhcx/plo+kY4dD837tunnfswnOO+88WUH6R1oCCDiSdwlG4LEoKiAAAAAAhFrnPk0OBSIVTSjhOL70zqU6ekBQAQEAAAAAQUcAAcfxLbGo0QOCJRgAAAAAYBsCCDiOL1+o0QPCtwTDDPWIAAAAAKDtI4CA4/h6PNTcBcN776EHBAAAAAAEXdgDiKeeekr9+/dXXFycRo8erffee6/Oc+fOnSvDMGrcRowY4T8nMzOz1nPKyspC8XEQAawGKiA8lEAAAAAAQNCFNYBYu3atbrvtNt11113auXOnJkyYoPT0dOXm5tZ6/vLly5WXl+e/HThwQF27dtXPf/7zgPM6deoUcF5eXp7i4uJC8ZEQATx1VUD4e0CEekQAAAAA0PaFNYBYunSprrvuOl1//fUaNmyYli1bpj59+mjFihW1np+QkKCkpCT/bfv27Tpy5IiuueaagPMMwwg4LykpKRQfBxHCtwSjZgWE954mlAAAAAAQfGELICoqKrRjxw5NmzYt4Pi0adO0bdu2Rl1j5cqVmjJlilJSUgKOl5SUKCUlRb1799aFF16onTt31nud8vJyFRcXB9zQdlknKhyqV0BE+ZpQ0gMCAGrFfAkAAFoibAFEQUGBPB6PEhMTA44nJiYqPz+/wdfn5eXp9ddf1/XXXx9wfOjQocrMzNQrr7yiNWvWKC4uTuecc4727NlT57UWLVqkhIQE/61Pnz7N+1CICP4KCFWrgDhRAlFFBQQA1Ir5EgAAtETYm1Aa1f4Z2rKsGsdqk5mZqc6dO2vmzJkBx8eNG6errrpKo0aN0oQJE/Tiiy9q8ODBevzxx+u8VkZGhoqKivy3AwcONOuzIDLU2QPCvw0nAQQA1Ib5EgAAtER0uN64e/fuioqKqlHtcOjQoRpVEdVZlqU//vGPmj17tmJjY+s91+Vy6cc//nG9FRBut1tut7vxg0dEO7kNJz0gAKApmC8BAEBLhK0CIjY2VqNHj9amTZsCjm/atEnjx4+v97Vbt27V3r17dd111zX4PpZlKScnR8nJyS0aL9oOX8AQVe3X798Fgx4QAAAAABB0YauAkKR58+Zp9uzZGjNmjNLS0vTss88qNzdXN910kyRvqefBgwe1evXqgNetXLlSY8eO1ciRI2tcc8GCBRo3bpwGDRqk4uJiPfbYY8rJydGTTz4Zks+E1q/KU9cuGL5tOAkgAAAAACDYwhpAzJo1S4cPH9bChQuVl5enkSNHauPGjf5dLfLy8pSbmxvwmqKiIq1bt07Lly+v9ZqFhYW68cYblZ+fr4SEBKWmpiorK0tnn3227Z8HkcG3BCPKxRIMAAAAAAgVw7KoN6+uuLhYCQkJKioqUqdOncI9HATZV9+XaPKSrbrvwuEamnzy/75llR5dk/mRlv/iR7roR6eFcYQAEBmYLwEAQFOEfRcMINR8FQ6uahUQvooIKiAAAAAAIPgIIOA4DfWAqCKAAAAAAICgI4CA4zTUA8IkgAAAAACAoCOAgOP4Khyq5Q8yDEOG2IYTAAAAAOxAAAHH8ZimpJoVEL5jVEAAAAAAQPARQMBxfD0gooyaAYTLMOgBAQAAAAA2IICA49S1C4b3GLtgAAAAAIAdCCDgOL4eD9V3wfAdM+kBAQAAAABBRwABx/EtsairB4THDPWIAAAAAKDtI4CA43g8te+C4T1m+JtUAgAAAACChwACjlNVXw8IQ1RAAAAAAIANCCDgOL4eD3XtguGhBwQAAAAABB0BBBynoR4QJrtgAAAAAEDQEUDAcXw9HuraBaOKAAIAAAAAgo4AAo5T5am7AsJliG04AQAAAMAGBBBwHF/AUNsuGIZhyEMFBAAAAAAEHQEEHKfKtOQyvGFDdVEuAggAAAAAsAMBBBzHY1q1bsEp+bbhJIAAAAAAgGAjgIDjVHmsWrfglNiGEwAAAADsQgABx/GYVq07YEiSyyW24QQAAAAAGxBAwHE8llXrDhjSiQoIAggAAAAACDoCCDiOxySAAAAAAIBQI4CA49ADAgAAAABCjwACjuMxTbnq+OWzCwYAAAAA2IMAAo7jseprQmnIpAICAAAAAIKOAAKOU9VAD4gqDwEEAAAAAAQbAQQcx+OppwLCEBUQAADUxjSl7KekooPhHgkAIEIRQMBxGqqAoAcEAAC1+PAZ6c0MaffGcI8EABChCCDgOPVuw+kigAAAoIb8z6RN93r/+3hhWIcCAIhc0eEeABBq3iaUtT/nMrwVEgAA4ATLktbfKHU6Taoqk44fCfeIAAARigoIOE79PSCogAAAIEDlMem7f0kjZkrxXQggAADNRgABx6kyLbnq6QFBE0oAAE5RXuK9j+kguTtIx38I73gAABGLAAKO4zFNRdVRARFFDwgAAAJVnAggYuOl2A5UQAAAmo0AAo7jsUQPCAAAGqv8qPc+pp0U25EKCABAsxFAwHE8pln/EgwCCAAATvJVQETHn1iCURjW4QAAIhcBBBynymPVuQTDMAx56AEBAMBJ/h4Q8ZK7o1RW5N0ZAwCAJiKAgON46mlCSQ8IAACq8S/BiPcuwfBUeHfGAACgiQgg4DhVZt0VEC5DBBAAAJyq4qgkQ4qO8y7BkGhECQBoFgIIOI7HtFRH/iCXQQUEAAAByku81Q+G4V2CIRFAAACahQACjlNlWoqqswmlRP4AAMApKkqk2Hbe/44lgAAANB8BBBzHY5p1LsGgBwQAANWUl0jRJwIIKiAAAC0Q9gDiqaeeUv/+/RUXF6fRo0frvffeq/PcLVu2yDCMGrcvvvgi4Lx169Zp+PDhcrvdGj58uNavX2/3x0AEqaqnCSVLMAAAqKbiqHcJhnSiEsIggAAANEtYA4i1a9fqtttu01133aWdO3dqwoQJSk9PV25ubr2v2717t/Ly8vy3QYMG+Z/Lzs7WrFmzNHv2bH388ceaPXu2Lr/8cn3wwQd2fxxECE89SzDYhhMAgGrKS6SYOO9/Gy5vFQQBBACgGcIaQCxdulTXXXedrr/+eg0bNkzLli1Tnz59tGLFinpf17NnTyUlJflvUVFR/ueWLVumqVOnKiMjQ0OHDlVGRoYmT56sZcuW2fxpECk8pqU68gd2wQAAoLqKEik6/uRjAggAQDOFLYCoqKjQjh07NG3atIDj06ZN07Zt2+p9bWpqqpKTkzV58mRt3rw54Lns7Owa1zz//PPrvWZ5ebmKi4sDbmi7qjyWXPX0gDAJIACgVsyXDlVWfHIJhiTFdpCOF4ZtOACAyBW2AKKgoEAej0eJiYkBxxMTE5Wfn1/ra5KTk/Xss89q3bp1eumllzRkyBBNnjxZWVlZ/nPy8/ObdE1JWrRokRISEvy3Pn36tOCTobWrbwkGPSAAoG7Mlw51ag8ISYptTwUEAKBZosM9AKPav0RbllXjmM+QIUM0ZMgQ/+O0tDQdOHBAixcv1rnnntusa0pSRkaG5s2b539cXFzMX6rasCrTrD+AoAcEANSK+dKhykukmHYnH7MEAwDQTGELILp3766oqKgalQmHDh2qUcFQn3Hjxum5557zP05KSmryNd1ut9xud6PfE5HNY9W9BMPlEkswAKAOzJcOVVFScwnGkf1hGw4AIHKFbQlGbGysRo8erU2bNgUc37Rpk8aPH9/o6+zcuVPJycn+x2lpaTWu+dZbbzXpmmjbvE0o6+gBQQUEAACBqgcQ7o5SGRUQAICmC+sSjHnz5mn27NkaM2aM0tLS9Oyzzyo3N1c33XSTJG+p58GDB7V69WpJ3h0u+vXrpxEjRqiiokLPPfec1q1bp3Xr1vmveeutt+rcc8/Vww8/rIsuukgvv/yy3n77bb3//vth+Yxofbw9IGp/zqAHBAAAJ3mqpKryWnbBKAzbkAAAkSusAcSsWbN0+PBhLVy4UHl5eRo5cqQ2btyolJQUSVJeXp5yc3P951dUVGj+/Pk6ePCg4uPjNWLECL322muaPn26/5zx48frhRde0N1336177rlHAwcO1Nq1azV27NiQfz60TvXtguFySabVcN8QAAAcoeKo9z6gAqKDVHnsRDDBkhwAQOOFvQnlzTffrJtvvrnW5zIzMwMe/+53v9Pvfve7Bq952WWX6bLLLgvG8NAGmVb9u2B4z5GiyB8AAE5XXuK9D+gB0dF7f7xQ6tj4vl0AAIStBwQQLlUN9ICQxDIMAAAkqdxXAVFtFwyJnTAAAE1GAAHH8faAqGsJBgEEAAB+FbVUQLg7eO8JIAAATUQAAcepbxcMXy7BThgAAKj2CgjfEoyywpAPBwAQ2Qgg4Dj17YLhX4LhIYAAAKDeCohjP4R+PACAiEYAAceprwIi+kQyUe7xhHJIAAC0TrU1oXRFe7flLCsKz5gAABGLAAKOYpqWLJ3s9VBdbPSJAKLSDOGoAABopSpKpKhYb+hwKncHlmAAAJqMAAKOUnWiuWRUHRUQsScqIMoqqYAAAEDlRwOrH3xi21MBAQBoMgIIOIpvd4u6dsHwVUCUUQEBAIC3AuLUBpQ+sR2k44UhHw4AILIRQMBRfLtb1NUDwl8BUUUFBAAA3gqIWgKImHYswQAANBkBBBzFt7uFq45ffmy0N5hgCQYAAPI2oYyJq3mcJRgAgGYggICjVJnepRV19YCIiWIJBgAAfhUl3h0vqnN3kI4fCf14AAARjQACjuLrAdHQLhhUQAAAoHqaUHagAgIA0GQEEHAUdsEAAKAJCCAAAEFEAAFHaWgXjCiXIZchlVWxBAMAAFUcrX0JRmx7qfKY5KkM/ZgAABGLAAKO4l+CUXv+IMMwFBvtUjkVEAAAnGhCWUcFhEQVBACgSQgg4ChVDfSAkLzLMFiCAQCApIrS2rfhjG3vvSeAAAA0AQEEHMXTQA8IyduIkl0wAACOZ1neXTBqrYDwBRCFIR0SACCyEUDAURraBUPyBRBUQAAAHK7yuGSZUkxczed8SzCOF4Z0SACAyEYAAUc52QOi/iUY5TShBAA4XeUx732tTSjpAQEAaDoCCDhKlekNFuraBUOiAgIAAEne5RdS7RUQMfGS4WIJBgCgSQgg4CiN6QERE+ViG04AACpKvfe1VUAYhrcKggoIAEATEEDAUU72gKj7HHbBAABApwQQtVRASN5GlPSAAAA0AQEEHKVRFRAswQAAoP4lGBIVEACAJiOAgKNUNWYXDCogAACQKnxNKOupgCCAAAA0AQEEHKVRu2BEu3S8kh4QAACH8y3BiKmlB4TEEgwAQJMRQMBRfBUQ9e6CEeVSORUQAACnqyiRjCjJFVP787Ht2QUDANAkBBBwlMb0gGAbTgAA5K2AiInz7nhRm9gOBBAAgCYhgICjNGYXDLbhBABAUuWx2rfg9KEJJQCgiQgg4ChVpjdYaKgHBBUQAADHqyipewcM6WQTSssK3ZgAABGNAAKO4ml0DwgqIAAADldR2nAFhFnlrZQAAKARCCDgKI3tAVHhMWWa/IsOAMDBKkqlaHfdz8e2996zDAMA0EgEEHAUXwBRT/6g2GjvH4ty+kAAAJysolSKrmcJhruD956tOAEAjUQAAUepMi25DMmorwIiyvvHgj4QAABHayiAoAICANBEBBBwFI9p1dv/QZJio73Pl1URQAAAHKzBJpQnKiDYihMA0EgEEHCUqsYEEP4KCJZgAAAcrKKkgSaUVEAAAJqGAAKOYppWvQ0opZM9IFiCAQBwtIpj9TehjIr1BhSlBaEbEwAgohFAwFEaUwERQw8IAAC8PSBi6qmAkKR2XaWS/NCMBwAQ8Qgg4Cge05Sr0RUQLMEAADhYZQNNKCUpvqt09LvQjAcAEPEIIOAoVaYlV2N7QNCEEgDgVKYpVR5vuAIivrN0NC8kQwIARD4CCDiK2ahdMLx/LMpZggEAcKrKY977+npASCcqIFiCAQBoHAIIOEpVk5pQsgQDAOBQFaXe+/p2wZBO9IBgCQYAoHHCHkA89dRT6t+/v+Li4jR69Gi99957dZ770ksvaerUqerRo4c6deqktLQ0vfnmmwHnZGZmyjCMGreysjK7PwoigMe01EABBE0oAQCoKPHexzSiB0R5sXfHDAAAGhDWAGLt2rW67bbbdNddd2nnzp2aMGGC0tPTlZubW+v5WVlZmjp1qjZu3KgdO3Zo0qRJmjFjhnbu3BlwXqdOnZSXlxdwi4trYAKFIzSmB4TLMBQTZRBAAACcq9EVEF289+yEAQBohOhwvvnSpUt13XXX6frrr5ckLVu2TG+++aZWrFihRYsW1Th/2bJlAY8feughvfzyy3r11VeVmprqP24YhpKSkmwdOyKTpxFLMCTvMoyyKpZgAAAcyt8DoqEKiG7e+6PfSV0H2DsmAEDEC1sFREVFhXbs2KFp06YFHJ82bZq2bdvWqGuYpqmjR4+qa9euAcdLSkqUkpKi3r1768ILL6xRIVFdeXm5iouLA25omzyNqICQJHdUFBUQAFAN86WDNHYJhq8Cgp0wAACNELYAoqCgQB6PR4mJiQHHExMTlZ/fuDK+JUuWqLS0VJdffrn/2NChQ5WZmalXXnlFa9asUVxcnM455xzt2bOnzussWrRICQkJ/lufPn2a96HQ6lWZlqIa8auPjXbRhBIAqmG+dJDGLsGIaS9FxdKIEgDQKGFvQmlUK4e3LKvGsdqsWbNG999/v9auXauePXv6j48bN05XXXWVRo0apQkTJujFF1/U4MGD9fjjj9d5rYyMDBUVFflvBw4caP4HQqvmMU25GvH7iok2VF5FBQQAnIr50kH8AUQD23AahtSuGxUQAIBGCVsPiO7duysqKqpGtcOhQ4dqVEVUt3btWl133XX6y1/+oilTptR7rsvl0o9//ON6KyDcbrfc7gYmWLQJVabVqAAiNooKCACojvnSQSpKJVe0FBXT8LnxXbw9IAAAaEDYKiBiY2M1evRobdq0KeD4pk2bNH78+Dpft2bNGs2dO1d//vOf9bOf/azB97EsSzk5OUpOTm7xmBH5TNNSVCN6QMRGu1RODwgAgFNVlEoxDSy/8InvIh1lFwwAQMPCugvGvHnzNHv2bI0ZM0ZpaWl69tlnlZubq5tuukmSt9Tz4MGDWr16tSRv+HD11Vdr+fLlGjdunL96Ij4+XgkJCZKkBQsWaNy4cRo0aJCKi4v12GOPKScnR08++WR4PiRalcZWQMREuVTGEgwAgFNVlDbc/8GnXVfp+y/tHQ8AoE0IawAxa9YsHT58WAsXLlReXp5GjhypjRs3KiUlRZKUl5en3Nxc//nPPPOMqqqqdMstt+iWW27xH58zZ44yMzMlSYWFhbrxxhuVn5+vhIQEpaamKisrS2effXZIPxtaJ+8uGA2fFxvlUlkFAQQAwKEqShveAcMnvqtUQgUEAKBhYQ0gJOnmm2/WzTffXOtzvlDBZ8uWLQ1e79FHH9Wjjz4ahJGhLaoyLbnUuAqI4/SAAAA4VUWJFN3IAKJdV6msSKo83vhlGwAARwr7LhhAKHk8ZqN6QLR3R6nweEUIRgQAQCtUeazxAUR8V+89W3ECABpAAAFH8ViSqxEBRJd2sfquuDwEIwIAoBWqKG1CANHFe89OGACABhBAwFGqTFNRjWhC2aVdrH4orVClh2UYAAAHKm/KEoxu3vujefaNBwDQJhBAwFGqPJYaUQChLu29+55/f5QqCACAA1WUNL4JZWwHKdotFf3b3jEBACIeAQQcxbsLRuMqICQpv7jM7iEBAND6VDZhCYZhSAl9pO+/sHdMAICIRwABR6kyLbkaswSjvTeAOEQAAQBwoopSKboJO1ok9JEO7bJvPACANoEAAo5yvMIjd3TDP/uO7mhFRxnKLyKAAAA4UEVp45dgSFKXFG8FhGXZNyYAQMQjgICjlJRXqV1sdIPnGYahru1i9R09IAAATtSUXTAkqXOK9zVFB+wbEwAg4hFAwFGOllUqPjaqUed6t+KkAgIA4DCeKslT4W0s2VidU7z3h+gDAQCoGwEEHMOyLJWWexQf07gAonO7GJZgAACc5/gR7727Y+Nf076HFBNPHwgAQL0IIOAYZZWmPJaldo2tgGhPBQQAwIGOFXjv4xIa/xrD8FZBsBMGAKAeBBBwjKPllZLU6AqILu1idaiYHhAAAIcpPRFAuJsQQEhS577Sd/8K/ngAAG0GAQQco6SsSpIa3QOia/tYHS2v0rGKKjuHBQBA69KcCghJSugrFeyWTE/wxwQAaBMIIOAYR08EEI1egtEuRpL0HVUQAAAnKS2QXFFSbPumva5LilRVLh3Zb8uwAACRjwACjlFSfqICoglLMCTRiBIA4CzHDnuXXxhG017n3wnj8+CPCQDQJhBAwDGONnEJhi+AOHSUAAIA4CClBU1ffiFJ8V28t39/FPwxAQDaBAIIOEZTKyDiY6PU3h2lAz8cs3NYAAC0LscKpLhOTX+dYUjJP5L2vh30IQEA2gYCCDhGSVmlYqNcio5q/M9+QPcO2nmg0L5BAQDQ2pQWSO5mBBCS1Oss6bvPpJJDwR0TAKBNIICAY5SUVzW6AaXPoMQO+uc3R2RZlk2jAgCglTnWzCUYktQr1Xv/1bvBGw8AoM0ggIBjHC2ranT/B59BPTvqyLFK7T/MMgwAgEOUHm5+ABHfWep2urT3naAOCQDQNhBAwDGONqMC4vSeHSRJ//zmiB1DAgCgdbEs6fgP3l0wmsvXB8I0gzYsAEDbQAABxygpq1JcIxtQ+nRwR6t3l3j9M5cAAgDgAGWFklnVvCaUPqed5Q0x8j8O2rAAAG0DAQQco6SsstE7YJxqUM8O2k4FBADACUoPe++buwRDknoMk2LaS7teCc6YAABtBgEEHONoeVUzA4iO2vPdUR0tq7RhVAAAtCLHCrz3LQkgomKkgT+VdmRKlWVBGRYAoG0ggIBjNKcJpSSN6NVJpiVt3v29DaMCAKAVKQ1CACFJQy/0LsP410stHxMAoM0ggIBjNDeA6NkpTqf37KBXcg7aMCoAAFqRYwWSDCm2Q8uuk3CadNpo6R8rvI0tAQAQAQQcpKS8Su2asQRDktIGdNOW3d+r6BjLMAAAbVhpgeTuJLmaN18GGHqhlP+JdOCDll8LANAmEEDAESzLUml58yogJGncgG7ymJbe/Fd+kEcGAEArcuxwy5df+Jw2WurcV3r391RBAAAkEUDAIcqrTFWZluJjo5v1+q7tYzW8VydtYBkGAKAtKy1o2RacpzJc0llzpP1Z0p5NwbkmACCiEUDAEY6WVUlSs3bB8JkwqLu2fXVY3xwuDdawAABoXXxLMIKl99lS0pnSW3dLnqrgXRcAEJEIIOAIvi00m7sEQ5LSBnRXB3e0nv8gN1jDAgCgdTn2ffAqICTJMKQx10oFu6Udq4J3XQBARCKAgCOUlLe8AiI22qWJg3to7UcHVFbpCdbQAABoPUoLgtcDwqfb6dLgC6S375MKDwT32gCAiEIAAUcoObEEo10LKiAkafKwnio6XqmX6QUBAGhrLEs6/kPwAwhJGn2NFB0vvXorDSkBwMEIIOAIR30VEC0MIJIT4pU2oJt+v/FzfVt4PBhDAwCgdSgtkKrKpfiuwb92bHsp7Rbpq3dYigEADkYAAUcoCUITSp9rz+mv2CiXbn1hJ0sxAABtx4EPvPc9hthz/d4/loZMlzb+t7Qvy573AAC0agQQcISS8irFRBmKiWr5T75DXLRumXS6duYW6qdLtuj1T/OCMEIAAMIsN1vq0FNq38O+9zj7RinpDGntVdL3u+17HwBAq0QAAUf4tvC4OsXFBO16Q5M66X8vPVPJneL16+f/qdvX5vgbXQIAEJFys6Uew+x9D1e0NPH/SfFdpJXTpNx/2Pt+AIBWhQACjrB59yEN7xXEbcUkJXeO1/zzh+jm8wbqjc/ydfnT2So8VhHU9wAAICQqjkl5H0s9R9j/XrEdpPP/R+rcR/rTDOmf/x+NKQHAIQgg0OYdLDyuL78rUWqfLrZcf8KgHrpvxnD9+8gx/fL//qEfSgkhAAAR5uAOyaySEoeH5v3cHaQpC6UB50mv/Kf051lSETtMAUBbRwCBNm/L7kNyGdKZvW3YVuyElG7tdffPhuvbwjJd9vQ2dsgAAESW3H94KxM69w3de0bFSON/I/30HunfH0mPp0pvL5COHwndGAAAIUUAgTbv3c8PaUhSR7V3R9v6Pn26ttN9M4arpKxKM5/8u97e9Z2t7wcAQNDkZks9h0lGGP5q2GesNPNpadhM6R9PSkuHSX+7XTrwoWSy2xQAtCVhDyCeeuop9e/fX3FxcRo9erTee++9es/funWrRo8erbi4OA0YMEBPP/10jXPWrVun4cOHy+12a/jw4Vq/fr1dw0crV1Jepb9/VaAf9e4ckvdLTojXfTNGqFfneF2/ersuW7FNK7Z8pZ25R1RRZYZkDAAANElZsXcLzp42N6CsT2w76azZ0iV/kEZcIv1rg7RyqvTIQOmFK6X3H5X2/12qpMIQACKZvf8k3IC1a9fqtttu01NPPaVzzjlHzzzzjNLT07Vr1y717VuzBHDfvn2aPn26brjhBj333HP6+9//rptvvlk9evTQpZdeKknKzs7WrFmz9MADD+jiiy/W+vXrdfnll+v999/X2LFjQ/0REUb7C0p103M7ZFnS2AHdQva+XdvH6nfnD9EH+37Qe3u+17K3v1R5lal2sVG68MxkzUw9TT/u1zUoW4ICANAiVeXe/4EvSSkTwjsWybs7xqhfSmdcLhV8Kf17u/T9F9JX73jDB1eMdNpZUv+JUsp4qfcYyd0x3KMGADSSYVnhazs8duxYnXXWWVqxYoX/2LBhwzRz5kwtWrSoxvn/7//9P73yyiv6/PPP/cduuukmffzxx8rOzpYkzZo1S8XFxXr99df951xwwQXq0qWL1qxZ06hxFRcXKyEhQUVFRerUKbg7JyB4TNPS0bIqlVd55HIZOlpWpbyi48o5UKisLwv04b7D6tkpTvOmDFafru3CNs4qj6n9h48p50Ch3tvzvQ4dLVf72CglJ8QrOspQ4fFKFR+vVFmlRy7DkDvapfbuaLV3Rys2yqXoKENRLkOd4mLUs6NbQ5M7anhygoYld1TX9rEyDCNsnw2AszFfRjBPpXTwn9LfH5X2viNNXSgljgz3qOpmeqTCb6Tv/iV995mU/6lUXuxdMtJ1gNR9sNQxyRtGWJZUVeYNLCqPe7f+dHeQ2veUEk6TOiZLHRKluE5SdLy3F4Ur2nuLipVcjfgHAss6sXOHJVnmyaUihuENSRpzDQBwoLBVQFRUVGjHjh264447Ao5PmzZN27Ztq/U12dnZmjZtWsCx888/XytXrlRlZaViYmKUnZ2t22+/vcY5y5Ytq3Ms5eXlKi8v9z8uKiqS5P2LVWu1M/cHXZu5Q5Ueyvrr0yvRrQ0f7gn3MPwGdnYpxrR04IejOnr0aK3nVEiq4ykAQdK/ezutuHK0eocxnAynjh07Niu8tH2+/OFr6YUrpKJ/B+d6aJzYjtKnb3hvESFG6pYqFR+UDu+Rvj1xA+BsV/xFSkkL2uWaO1eifmELIAoKCuTxeJSYmBhwPDExUfn5+bW+Jj8/v9bzq6qqVFBQoOTk5DrPqeuakrRo0SItWLCgxvE+ffo09uOglToQ7gEAaJUOSBpxd7hHET7NrVhgvmyrjkp6JdyDAICW+Z8Lgno5qvvsEdYeEJJqpEqWZdWbNNV2fvXjTb1mRkaG5s2b539smqZ++OEHdevWLWSpV3Fxsfr06aMDBw7wQw8hvvfw4HsPD7738GiN33vHjs1bMx/u+bI1fpdOwPceHnzv4cH3Hh6t8Xtv7lyJ+oUtgOjevbuioqJqVCYcOnSoRgWDT1JSUq3nR0dHq1u3bvWeU9c1Jcntdsvtdgcc69y5c2M/SlB16tSp1fyhcxK+9/Dgew8PvvfwaAvfe2uZL9vCdxmJ+N7Dg+89PPjew4Pvve0LW4ec2NhYjR49Wps2bQo4vmnTJo0fP77W16SlpdU4/6233tKYMWMUExNT7zl1XRMAAAAAANgvrEsw5s2bp9mzZ2vMmDFKS0vTs88+q9zcXN10002SvKWeBw8e1OrVqyV5d7x44oknNG/ePN1www3Kzs7WypUrA3a3uPXWW3Xuuefq4Ycf1kUXXaSXX35Zb7/9tt5///2wfEYAAAAAABDmAGLWrFk6fPiwFi5cqLy8PI0cOVIbN25USkqKJCkvL0+5ubn+8/v376+NGzfq9ttv15NPPqlevXrpscce06WXXuo/Z/z48XrhhRd0991365577tHAgQO1du1ajR07NuSfryncbrfuu+++GqWtsBffe3jwvYcH33t48L0HD99lePC9hwffe3jwvYcH37tzGJaviyMAAAAAAIBNwtYDAgAAAAAAOAcBBAAAAAAAsB0BBAAAAAAAsB0BBAAAAAAAsB0BRJgcOXJEs2fPVkJCghISEjR79mwVFhbW+5q5c+fKMIyA27hx40Iz4Aj21FNPqX///oqLi9Po0aP13nvv1Xv+1q1bNXr0aMXFxWnAgAF6+umnQzTStqUp3/uWLVtq/LYNw9AXX3wRwhFHtqysLM2YMUO9evWSYRjasGFDg6/ht95yTf3e+a03HfNlaDBXhgdzZegxX4YH8yV8CCDC5IorrlBOTo7eeOMNvfHGG8rJydHs2bMbfN0FF1ygvLw8/23jxo0hGG3kWrt2rW677Tbddddd2rlzpyZMmKD09PSA7V1PtW/fPk2fPl0TJkzQzp07deedd+o3v/mN1q1bF+KRR7amfu8+u3fvDvh9Dxo0KEQjjnylpaUaNWqUnnjiiUadz289OJr6vfvwW2885kv7MVeGB3NleDBfhgfzJfwshNyuXbssSdY//vEP/7Hs7GxLkvXFF1/U+bo5c+ZYF110UQhG2HacffbZ1k033RRwbOjQodYdd9xR6/m/+93vrKFDhwYc+9WvfmWNGzfOtjG2RU393jdv3mxJso4cORKC0bV9kqz169fXew6/9eBrzPfOb71pmC9Dg7kyPJgrw4/5MjyYL52NCogwyM7OVkJCgsaOHes/Nm7cOCUkJGjbtm31vnbLli3q2bOnBg8erBtuuEGHDh2ye7gRq6KiQjt27NC0adMCjk+bNq3O7zk7O7vG+eeff762b9+uyspK28baljTne/dJTU1VcnKyJk+erM2bN9s5TMfjtx5e/NYbh/nSfsyV4cFcGTn4vYcXv/e2hwAiDPLz89WzZ88ax3v27Kn8/Pw6X5eenq7nn39e7777rpYsWaKPPvpIP/3pT1VeXm7ncCNWQUGBPB6PEhMTA44nJibW+T3n5+fXen5VVZUKCgpsG2tb0pzvPTk5Wc8++6zWrVunl156SUOGDNHkyZOVlZUViiE7Er/18OC33jTMl/ZjrgwP5srIwe89PPi9t13R4R5AW3L//fdrwYIF9Z7z0UcfSZIMw6jxnGVZtR73mTVrlv+/R44cqTFjxiglJUWvvfaaLrnkkmaOuu2r/p029D3Xdn5tx1G/pnzvQ4YM0ZAhQ/yP09LSdODAAS1evFjnnnuureN0Mn7rocdv3Yv5svVhrgwP5srIwO899Pi9t10EEEH0n//5n/rFL35R7zn9+vXTJ598ou+++67Gc99//32NhLU+ycnJSklJ0Z49e5o8Vifo3r27oqKiavxLwqFDh+r8npOSkmo9Pzo6Wt26dbNtrG1Jc7732owbN07PPfdcsIeHE/ittx5O/K0zX7YezJXhwVwZOfi9tx783tsGAogg6t69u7p3797geWlpaSoqKtKHH36os88+W5L0wQcfqKioSOPHj2/0+x0+fFgHDhxQcnJys8fclsXGxmr06NHatGmTLr74Yv/xTZs26aKLLqr1NWlpaXr11VcDjr311lsaM2aMYmJibB1vW9Gc7702O3fu5LdtI37rrYcTf+vMl60Hc2V4MFdGDn7vrQe/9zYiXN0vne6CCy6wzjzzTCs7O9vKzs62zjjjDOvCCy8MOGfIkCHWSy+9ZFmWZR09etT67W9/a23bts3at2+ftXnzZistLc067bTTrOLi4nB8hIjwwgsvWDExMdbKlSutXbt2WbfddpvVvn17a//+/ZZlWdYdd9xhzZ4923/+119/bbVr1866/fbbrV27dlkrV660YmJirL/+9a/h+ggRqanf+6OPPmqtX7/e+vLLL63PPvvMuuOOOyxJ1rp168L1ESLO0aNHrZ07d1o7d+60JFlLly61du7caX3zzTeWZfFbt0tTv3d+603HfGk/5srwYK4MD+bL8GC+hA8BRJgcPnzYuvLKK62OHTtaHTt2tK688soa28xIslatWmVZlmUdO3bMmjZtmtWjRw8rJibG6tu3rzVnzhwrNzc39IOPME8++aSVkpJixcbGWmeddZa1detW/3Nz5syxJk6cGHD+li1brNTUVCs2Ntbq16+ftWLFihCPuG1oyvf+8MMPWwMHDrTi4uKsLl26WD/5yU+s1157LQyjjly+7aqq3+bMmWNZFr91uzT1e+e33nTMl6HBXBkezJWhx3wZHsyX8DEs60QXFQAAAAAAAJuwDScAAAAAALAdAQQAAAAAALAdAQQAAAAAALAdAQQAAAAAALAdAQQAAAAAALAdAQQAAAAAALAdAQQAAAAAALAdAQQAAK1EVlaWZsyYoV69eskwDG3YsKHJ17AsS4sXL9bgwYPldrvVp08fPfTQQ8EfLAAAYcJ8GbkIIABEnGeffVZ9+vSRy+XSsmXLwj0cIGhKS0s1atQoPfHEE82+xq233qo//OEPWrx4sb744gu9+uqrOvvss4M4SgCRgLkSbRnzZeQyLMuywj0IAG3f/fffrw0bNignJ6fRrzEMQ+vXr9fMmTP9x4qLi9W9e3ctXbpUl156qRISEtSuXbvgDxgIs9p+/xUVFbr77rv1/PPPq7CwUCNHjtTDDz+s8847T5L0+eef68wzz9Rnn32mIUOGhGfgAJqNuRJoOubLyEIFBICIkpubq8rKSv3sZz9TcnIyf6GCo1xzzTX6+9//rhdeeEGffPKJfv7zn+uCCy7Qnj17JEmvvvqqBgwYoL/97W/q37+/+vXrp+uvv14//PBDmEcOIJSYK+F0zJetmAUAjTBx4kTrv/7rv6z//u//trp06WIlJiZa9913n//5wsJC64YbbrB69OhhdezY0Zo0aZKVk5NjWZZlrVq1ypIUcFu1alW975eSkhJwfkpKSq3X2bdvn7V3717rP/7jP6yePXta7du3t8aMGWNt2rTJxm8DsJ8ka/369f7He/futQzDsA4ePBhw3uTJk62MjAzLsizrV7/6leV2u62xY8daWVlZ1ubNm60f/ehH1qRJk0I5dMCxmCuB0GO+jCzRoYs6AES6P/3pT5o3b54++OADZWdna+7cuTrnnHM0ZcoU/exnP1PXrl21ceNGJSQk6JlnntHkyZP15ZdfatasWfrss8/0xhtv6O2335YkJSQk1PteH330kXr27KlVq1bpggsuUFRUlDp06KA+ffpoypQp+vDDD9WnTx/16NFDn332maZPn64HH3xQcXFx+tOf/qQZM2Zo9+7d6tu3byi+GsB2//znP2VZlgYPHhxwvLy8XN26dZMkmaap8vJyrV692n/eypUrNXr0aO3evZsyUyAEmCuB8GK+bN0IIAA02plnnqn77rtPkjRo0CA98cQTeueddxQVFaVPP/1Uhw4dktvtliQtXrxYGzZs0F//+lfdeOON6tChg6Kjo5WUlNSo9+rRo4ckqXPnzgGv8U0cPXr08B8fNWqURo0a5T/nwQcf1Pr16/XKK6/oP//zP1v+wYFWwDRNRUVFaceOHYqKigp4rkOHDpKk5ORkRUdHB/yla9iwYZK8Jdn8hQqwH3MlEF7Ml60bAQSARjvzzDMDHicnJ+vQoUPasWOHSkpK/H/h8Tl+/Li++uor28dVWlqqBQsW6G9/+5u+/fZbVVVV6fjx48rNzbX9vYFQSU1Nlcfj0aFDhzRhwoRazznnnHNUVVWlr776SgMHDpQkffnll5KklJSUkI0VcDLmSiC8mC9bNwIIAI0WExMT8NgwDJmmKdM0lZycrC1bttR4TefOnW0f13//93/rzTff1OLFi3X66acrPj5el112mSoqKmx/byCYSkpKtHfvXv/jffv2KScnR127dtXgwYN15ZVX6uqrr9aSJUuUmpqqgoICvfvuuzrjjDM0ffp0TZkyRWeddZauvfZaLVu2TKZp6pZbbtHUqVNrlKICsAdzJWA/5svIRQABoMXOOuss5efnKzo6Wv369av1nNjYWHk8niZdNyYmplGvee+99zR37lxdfPHFkryT0v79+5v0XkBrsH37dk2aNMn/eN68eZKkOXPmKDMzU6tWrdKDDz6o3/72tzp48KC6deumtLQ0TZ8+XZLkcrn06quv6r/+67907rnnqn379kpPT9eSJUvC8nkAnMRcCQQP82XkIoAA0GJTpkxRWlqaZs6cqYcfflhDhgzRt99+q40bN2rmzJkaM2aM+vXr50+ne/furY4dO/rXwNalX79+euedd3TOOefI7XarS5cutZ53+umn66WXXtKMGTNkGIbuuecemaZpx0cFbHXeeefJsqw6n4+JidGCBQu0YMGCOs/p1auX1q1bZ8fwALQAcyUQPMyXkcsV7gEAiHyGYWjjxo0699xzde2112rw4MH6xS9+of379ysxMVGSdOmll+qCCy7QpEmT1KNHD61Zs6bB6y5ZskSbNm1Snz59lJqaWud5jz76qLp06aLx48drxowZOv/883XWWWcF7fMBANBSzJUAIBlWfdERAAAAAABAEFABAQAAAAAAbEcAASAsnn/+eXXo0KHW24gRI8I9PAAAwo65EkBbwxIMAGFx9OhRfffdd7U+FxMTwx7MAADHY64E0NYQQAAAAAAAANuxBAMAAAAAANiOAAIAAAAAANiOAAIAAAAAANiOAAIAAAAAANiOAAIAAAAAANiOAAIAAAAAANiOAAIAAAAAANju/wd2eFA2BGOG2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1058.88x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.displot(data=data, x='net_tfa', kind='kde', col='e401', hue='e401', fill=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "awful-antigua",
   "metadata": {
    "papermill": {
     "duration": 0.040645,
     "end_time": "2021-03-24T14:13:54.245822",
     "exception": false,
     "start_time": "2021-03-24T14:13:54.205177",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The unconditional APE of e401 is about $19559$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ef070d79",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19559\n"
     ]
    }
   ],
   "source": [
    "e1 = data[data['e401'] == 1]['net_tfa']\n",
    "e0 = data[data['e401'] == 0]['net_tfa']\n",
    "print(f'{np.mean(e1) - np.mean(e0):.0f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cross-priority",
   "metadata": {
    "papermill": {
     "duration": 0.041187,
     "end_time": "2021-03-24T14:13:54.437695",
     "exception": false,
     "start_time": "2021-03-24T14:13:54.396508",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Among the $3682$ individuals that  are eligible, $2594$ decided to participate in the program. The unconditional APE of p401 is about $27372$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "33cd014e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "27372\n"
     ]
    }
   ],
   "source": [
    "e1 = data[data['p401'] == 1]['net_tfa']\n",
    "e0 = data[data['p401'] == 0]['net_tfa']\n",
    "print(f'{np.mean(e1) - np.mean(e0):.0f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "suitable-vulnerability",
   "metadata": {
    "papermill": {
     "duration": 0.041678,
     "end_time": "2021-03-24T14:13:54.640532",
     "exception": false,
     "start_time": "2021-03-24T14:13:54.598854",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "As discussed, these estimates are biased since they do not account for saver heterogeneity and endogeneity of participation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a9354414",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['age', 'inc', 'fsize', 'educ', 'db', 'marr', 'male', 'twoearn', 'pira',\n",
       "       'nohs', 'hs', 'smcol', 'col', 'hown'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = data['net_tfa'].values\n",
    "D = data['e401'].values\n",
    "D2 = data['p401'].values\n",
    "D3 = data['a401'].values\n",
    "X = data.drop(['e401', 'p401', 'a401', 'tw', 'tfa', 'net_tfa', 'tfa_he',\n",
    "               'hval', 'hmort', 'hequity',\n",
    "               'nifa', 'net_nifa', 'net_n401', 'ira',\n",
    "               'dum91', 'icat', 'ecat', 'zhat',\n",
    "               'i1', 'i2', 'i3', 'i4', 'i5', 'i6', 'i7',\n",
    "               'a1', 'a2', 'a3', 'a4', 'a5'], axis=1)\n",
    "X.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4c3e489",
   "metadata": {},
   "source": [
    "### We define a transformer that constructs the engineered features for controls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a2cfdf4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defaulting to user installation because normal site-packages is not writeable\n",
      "Requirement already satisfied: formulaic in c:\\programdata\\anaconda3\\lib\\site-packages (0.5.2)\n",
      "Requirement already satisfied: wrapt>=1.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from formulaic) (1.14.1)\n",
      "Requirement already satisfied: scipy>=1.6 in c:\\programdata\\anaconda3\\lib\\site-packages (from formulaic) (1.9.1)\n",
      "Requirement already satisfied: pandas>=1.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from formulaic) (1.4.4)\n",
      "Requirement already satisfied: interface-meta>=1.2.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from formulaic) (1.3.0)\n",
      "Requirement already satisfied: typing-extensions>=4.2.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from formulaic) (4.3.0)\n",
      "Requirement already satisfied: numpy>=1.16.5 in c:\\programdata\\anaconda3\\lib\\site-packages (from formulaic) (1.21.5)\n",
      "Requirement already satisfied: astor>=0.8 in c:\\programdata\\anaconda3\\lib\\site-packages (from formulaic) (0.8.1)\n",
      "Requirement already satisfied: pytz>=2020.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from pandas>=1.0->formulaic) (2022.1)\n",
      "Requirement already satisfied: python-dateutil>=2.8.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from pandas>=1.0->formulaic) (2.8.2)\n",
      "Requirement already satisfied: six>=1.5 in c:\\programdata\\anaconda3\\lib\\site-packages (from python-dateutil>=2.8.1->pandas>=1.0->formulaic) (1.16.0)\n"
     ]
    }
   ],
   "source": [
    "!pip install formulaic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "53b1283d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.base import TransformerMixin, BaseEstimator\n",
    "from formulaic import Formula\n",
    "\n",
    "class FormulaTransformer(TransformerMixin, BaseEstimator):\n",
    "    \n",
    "    def __init__(self, formula, array=False):\n",
    "        self.formula = formula\n",
    "        self.array = array\n",
    "    \n",
    "    def fit(self, X, y=None):\n",
    "        return self\n",
    "    \n",
    "    def transform(self, X, y=None):\n",
    "        df = Formula(self.formula).get_model_matrix(X)\n",
    "        if self.array:\n",
    "            return df.values\n",
    "        return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "824ee320",
   "metadata": {},
   "outputs": [],
   "source": [
    "transformer = FormulaTransformer(\"0 + poly(age, degree=6, raw=True) + poly(inc, degree=8, raw=True) \"\n",
    "                                 \"+ poly(educ, degree=4, raw=True) + poly(fsize, degree=2, raw=True) \"\n",
    "                                 \"+ male + marr + twoearn + db + pira + hown\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b35bbd34",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>db</th>\n",
       "      <th>hown</th>\n",
       "      <th>male</th>\n",
       "      <th>marr</th>\n",
       "      <th>pira</th>\n",
       "      <th>poly(age, degree=6, raw=True)[0]</th>\n",
       "      <th>poly(age, degree=6, raw=True)[1]</th>\n",
       "      <th>poly(age, degree=6, raw=True)[2]</th>\n",
       "      <th>poly(age, degree=6, raw=True)[3]</th>\n",
       "      <th>poly(age, degree=6, raw=True)[4]</th>\n",
       "      <th>...</th>\n",
       "      <th>poly(fsize, degree=2, raw=True)[1]</th>\n",
       "      <th>poly(inc, degree=8, raw=True)[0]</th>\n",
       "      <th>poly(inc, degree=8, raw=True)[1]</th>\n",
       "      <th>poly(inc, degree=8, raw=True)[2]</th>\n",
       "      <th>poly(inc, degree=8, raw=True)[3]</th>\n",
       "      <th>poly(inc, degree=8, raw=True)[4]</th>\n",
       "      <th>poly(inc, degree=8, raw=True)[5]</th>\n",
       "      <th>poly(inc, degree=8, raw=True)[6]</th>\n",
       "      <th>poly(inc, degree=8, raw=True)[7]</th>\n",
       "      <th>twoearn</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>...</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9915.000000</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9.915000e+03</td>\n",
       "      <td>9915.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>0.271004</td>\n",
       "      <td>0.635199</td>\n",
       "      <td>0.205951</td>\n",
       "      <td>0.604841</td>\n",
       "      <td>0.242158</td>\n",
       "      <td>41.060212</td>\n",
       "      <td>1792.938981</td>\n",
       "      <td>82850.255371</td>\n",
       "      <td>4.022991e+06</td>\n",
       "      <td>2.036867e+08</td>\n",
       "      <td>...</td>\n",
       "      <td>10.581241</td>\n",
       "      <td>37200.623197</td>\n",
       "      <td>1.997590e+09</td>\n",
       "      <td>1.437534e+14</td>\n",
       "      <td>9.104568e+17</td>\n",
       "      <td>-1.945721e+16</td>\n",
       "      <td>-9.132838e+16</td>\n",
       "      <td>1.934747e+16</td>\n",
       "      <td>2.234841e+16</td>\n",
       "      <td>0.380837</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.444500</td>\n",
       "      <td>0.481399</td>\n",
       "      <td>0.404415</td>\n",
       "      <td>0.488909</td>\n",
       "      <td>0.428411</td>\n",
       "      <td>10.344505</td>\n",
       "      <td>899.134499</td>\n",
       "      <td>61725.701911</td>\n",
       "      <td>3.939869e+06</td>\n",
       "      <td>2.448502e+08</td>\n",
       "      <td>...</td>\n",
       "      <td>11.391608</td>\n",
       "      <td>24774.288006</td>\n",
       "      <td>3.026128e+09</td>\n",
       "      <td>4.129207e+14</td>\n",
       "      <td>3.249753e+18</td>\n",
       "      <td>5.290488e+18</td>\n",
       "      <td>5.317962e+18</td>\n",
       "      <td>5.317873e+18</td>\n",
       "      <td>5.262071e+18</td>\n",
       "      <td>0.485617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>25.000000</td>\n",
       "      <td>625.000000</td>\n",
       "      <td>15625.000000</td>\n",
       "      <td>3.906250e+05</td>\n",
       "      <td>9.765625e+06</td>\n",
       "      <td>...</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-2652.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-1.865179e+10</td>\n",
       "      <td>-9.222694e+18</td>\n",
       "      <td>-9.221958e+18</td>\n",
       "      <td>-9.221276e+18</td>\n",
       "      <td>-9.223372e+18</td>\n",
       "      <td>-9.221563e+18</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>32.000000</td>\n",
       "      <td>1024.000000</td>\n",
       "      <td>32768.000000</td>\n",
       "      <td>1.048576e+06</td>\n",
       "      <td>3.355443e+07</td>\n",
       "      <td>...</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>19413.000000</td>\n",
       "      <td>3.768646e+08</td>\n",
       "      <td>7.316072e+12</td>\n",
       "      <td>3.677463e+16</td>\n",
       "      <td>-4.537058e+18</td>\n",
       "      <td>-4.603302e+18</td>\n",
       "      <td>-4.539584e+18</td>\n",
       "      <td>-4.508422e+18</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>40.000000</td>\n",
       "      <td>1600.000000</td>\n",
       "      <td>64000.000000</td>\n",
       "      <td>2.560000e+06</td>\n",
       "      <td>1.024000e+08</td>\n",
       "      <td>...</td>\n",
       "      <td>9.000000</td>\n",
       "      <td>31476.000000</td>\n",
       "      <td>9.907386e+08</td>\n",
       "      <td>3.118449e+13</td>\n",
       "      <td>4.226261e+17</td>\n",
       "      <td>3.660147e+16</td>\n",
       "      <td>-2.199258e+17</td>\n",
       "      <td>3.665403e+16</td>\n",
       "      <td>4.447478e+16</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>48.000000</td>\n",
       "      <td>2304.000000</td>\n",
       "      <td>110592.000000</td>\n",
       "      <td>5.308416e+06</td>\n",
       "      <td>2.548040e+08</td>\n",
       "      <td>...</td>\n",
       "      <td>16.000000</td>\n",
       "      <td>48583.500000</td>\n",
       "      <td>2.360356e+09</td>\n",
       "      <td>1.146744e+14</td>\n",
       "      <td>2.117477e+18</td>\n",
       "      <td>4.485390e+18</td>\n",
       "      <td>4.471642e+18</td>\n",
       "      <td>4.613388e+18</td>\n",
       "      <td>4.509509e+18</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>64.000000</td>\n",
       "      <td>4096.000000</td>\n",
       "      <td>262144.000000</td>\n",
       "      <td>1.677722e+07</td>\n",
       "      <td>1.073742e+09</td>\n",
       "      <td>...</td>\n",
       "      <td>169.000000</td>\n",
       "      <td>242124.000000</td>\n",
       "      <td>5.862403e+10</td>\n",
       "      <td>1.419428e+16</td>\n",
       "      <td>9.222042e+18</td>\n",
       "      <td>9.218801e+18</td>\n",
       "      <td>9.219035e+18</td>\n",
       "      <td>9.223289e+18</td>\n",
       "      <td>9.220028e+18</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>8 rows × 26 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                db         hown         male         marr         pira  \\\n",
       "count  9915.000000  9915.000000  9915.000000  9915.000000  9915.000000   \n",
       "mean      0.271004     0.635199     0.205951     0.604841     0.242158   \n",
       "std       0.444500     0.481399     0.404415     0.488909     0.428411   \n",
       "min       0.000000     0.000000     0.000000     0.000000     0.000000   \n",
       "25%       0.000000     0.000000     0.000000     0.000000     0.000000   \n",
       "50%       0.000000     1.000000     0.000000     1.000000     0.000000   \n",
       "75%       1.000000     1.000000     0.000000     1.000000     0.000000   \n",
       "max       1.000000     1.000000     1.000000     1.000000     1.000000   \n",
       "\n",
       "       poly(age, degree=6, raw=True)[0]  poly(age, degree=6, raw=True)[1]  \\\n",
       "count                       9915.000000                       9915.000000   \n",
       "mean                          41.060212                       1792.938981   \n",
       "std                           10.344505                        899.134499   \n",
       "min                           25.000000                        625.000000   \n",
       "25%                           32.000000                       1024.000000   \n",
       "50%                           40.000000                       1600.000000   \n",
       "75%                           48.000000                       2304.000000   \n",
       "max                           64.000000                       4096.000000   \n",
       "\n",
       "       poly(age, degree=6, raw=True)[2]  poly(age, degree=6, raw=True)[3]  \\\n",
       "count                       9915.000000                      9.915000e+03   \n",
       "mean                       82850.255371                      4.022991e+06   \n",
       "std                        61725.701911                      3.939869e+06   \n",
       "min                        15625.000000                      3.906250e+05   \n",
       "25%                        32768.000000                      1.048576e+06   \n",
       "50%                        64000.000000                      2.560000e+06   \n",
       "75%                       110592.000000                      5.308416e+06   \n",
       "max                       262144.000000                      1.677722e+07   \n",
       "\n",
       "       poly(age, degree=6, raw=True)[4]  ...  \\\n",
       "count                      9.915000e+03  ...   \n",
       "mean                       2.036867e+08  ...   \n",
       "std                        2.448502e+08  ...   \n",
       "min                        9.765625e+06  ...   \n",
       "25%                        3.355443e+07  ...   \n",
       "50%                        1.024000e+08  ...   \n",
       "75%                        2.548040e+08  ...   \n",
       "max                        1.073742e+09  ...   \n",
       "\n",
       "       poly(fsize, degree=2, raw=True)[1]  poly(inc, degree=8, raw=True)[0]  \\\n",
       "count                         9915.000000                       9915.000000   \n",
       "mean                            10.581241                      37200.623197   \n",
       "std                             11.391608                      24774.288006   \n",
       "min                              1.000000                      -2652.000000   \n",
       "25%                              4.000000                      19413.000000   \n",
       "50%                              9.000000                      31476.000000   \n",
       "75%                             16.000000                      48583.500000   \n",
       "max                            169.000000                     242124.000000   \n",
       "\n",
       "       poly(inc, degree=8, raw=True)[1]  poly(inc, degree=8, raw=True)[2]  \\\n",
       "count                      9.915000e+03                      9.915000e+03   \n",
       "mean                       1.997590e+09                      1.437534e+14   \n",
       "std                        3.026128e+09                      4.129207e+14   \n",
       "min                        0.000000e+00                     -1.865179e+10   \n",
       "25%                        3.768646e+08                      7.316072e+12   \n",
       "50%                        9.907386e+08                      3.118449e+13   \n",
       "75%                        2.360356e+09                      1.146744e+14   \n",
       "max                        5.862403e+10                      1.419428e+16   \n",
       "\n",
       "       poly(inc, degree=8, raw=True)[3]  poly(inc, degree=8, raw=True)[4]  \\\n",
       "count                      9.915000e+03                      9.915000e+03   \n",
       "mean                       9.104568e+17                     -1.945721e+16   \n",
       "std                        3.249753e+18                      5.290488e+18   \n",
       "min                       -9.222694e+18                     -9.221958e+18   \n",
       "25%                        3.677463e+16                     -4.537058e+18   \n",
       "50%                        4.226261e+17                      3.660147e+16   \n",
       "75%                        2.117477e+18                      4.485390e+18   \n",
       "max                        9.222042e+18                      9.218801e+18   \n",
       "\n",
       "       poly(inc, degree=8, raw=True)[5]  poly(inc, degree=8, raw=True)[6]  \\\n",
       "count                      9.915000e+03                      9.915000e+03   \n",
       "mean                      -9.132838e+16                      1.934747e+16   \n",
       "std                        5.317962e+18                      5.317873e+18   \n",
       "min                       -9.221276e+18                     -9.223372e+18   \n",
       "25%                       -4.603302e+18                     -4.539584e+18   \n",
       "50%                       -2.199258e+17                      3.665403e+16   \n",
       "75%                        4.471642e+18                      4.613388e+18   \n",
       "max                        9.219035e+18                      9.223289e+18   \n",
       "\n",
       "       poly(inc, degree=8, raw=True)[7]      twoearn  \n",
       "count                      9.915000e+03  9915.000000  \n",
       "mean                       2.234841e+16     0.380837  \n",
       "std                        5.262071e+18     0.485617  \n",
       "min                       -9.221563e+18     0.000000  \n",
       "25%                       -4.508422e+18     0.000000  \n",
       "50%                        4.447478e+16     0.000000  \n",
       "75%                        4.509509e+18     1.000000  \n",
       "max                        9.220028e+18     1.000000  \n",
       "\n",
       "[8 rows x 26 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transformer.fit_transform(X).describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8c5e20e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "transformer = FormulaTransformer(\"0 + poly(age, degree=6, raw=True) + poly(inc, degree=8, raw=True) \"\n",
    "                                 \"+ poly(educ, degree=4, raw=True) + poly(fsize, degree=2, raw=True) \"\n",
    "                                 \"+ male + marr + twoearn + db + pira + hown\", array=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f751665e",
   "metadata": {},
   "source": [
    "# Effect of Eligibility on Financial Assets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "051af0ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "modely = make_pipeline(transformer, StandardScaler(), LassoCV())\n",
    "modeld = make_pipeline(transformer, StandardScaler(), LassoCV())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "53bc1c83",
   "metadata": {},
   "outputs": [],
   "source": [
    "resy = y - modely.fit(X, y).predict(X)\n",
    "resD = D - modeld.fit(X, D).predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "51244ab4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9018.119246620929"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(resy * resD) / np.mean(resD**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de8e13e7",
   "metadata": {},
   "source": [
    "# Instrumental Variables: Effect of 401k on Financial Assets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "93b53b7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "modely = make_pipeline(transformer, StandardScaler(), LassoCV())\n",
    "modeld = make_pipeline(transformer, StandardScaler(), LassoCV())\n",
    "modelz = make_pipeline(transformer, StandardScaler(), LassoCV())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "629efff8",
   "metadata": {},
   "outputs": [],
   "source": [
    "resy = y - modely.fit(X, y).predict(X)\n",
    "resD = D - modeld.fit(X, D).predict(X) # instrument is e401k (eligibility)\n",
    "resD2 = D2 - modelz.fit(X, D2).predict(X) # treatment is p401k (participation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "08b2d78e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12959.476983995874"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(resy * resD) / np.mean(resD2 * resD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "825db7e2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.13"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 537.547458,
   "end_time": "2021-03-24T14:22:44.931595",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2021-03-24T14:13:47.384137",
   "version": "2.3.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}