Banking_Recovery.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Regression discontinuity: banking recovery\n",
    "<p>First, we'll load the banking dataset and look at the first few rows of data. This puts us in a good position to understand the dataset itself and begin thinking about how to analyze the data.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>id</th>\n",
       "      <th>expected_recovery_amount</th>\n",
       "      <th>actual_recovery_amount</th>\n",
       "      <th>recovery_strategy</th>\n",
       "      <th>age</th>\n",
       "      <th>sex</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2030</td>\n",
       "      <td>194</td>\n",
       "      <td>263.540</td>\n",
       "      <td>Level 0 Recovery</td>\n",
       "      <td>19</td>\n",
       "      <td>Male</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1150</td>\n",
       "      <td>486</td>\n",
       "      <td>416.090</td>\n",
       "      <td>Level 0 Recovery</td>\n",
       "      <td>25</td>\n",
       "      <td>Female</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>380</td>\n",
       "      <td>527</td>\n",
       "      <td>429.350</td>\n",
       "      <td>Level 0 Recovery</td>\n",
       "      <td>27</td>\n",
       "      <td>Male</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1838</td>\n",
       "      <td>536</td>\n",
       "      <td>296.990</td>\n",
       "      <td>Level 0 Recovery</td>\n",
       "      <td>25</td>\n",
       "      <td>Male</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1995</td>\n",
       "      <td>541</td>\n",
       "      <td>346.385</td>\n",
       "      <td>Level 0 Recovery</td>\n",
       "      <td>34</td>\n",
       "      <td>Male</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     id  expected_recovery_amount  actual_recovery_amount recovery_strategy  \\\n",
       "0  2030                       194                 263.540  Level 0 Recovery   \n",
       "1  1150                       486                 416.090  Level 0 Recovery   \n",
       "2   380                       527                 429.350  Level 0 Recovery   \n",
       "3  1838                       536                 296.990  Level 0 Recovery   \n",
       "4  1995                       541                 346.385  Level 0 Recovery   \n",
       "\n",
       "   age     sex  \n",
       "0   19    Male  \n",
       "1   25  Female  \n",
       "2   27    Male  \n",
       "3   25    Male  \n",
       "4   34    Male  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read in dataset\n",
    "df = pd.read_csv(\"bank_data.csv\")\n",
    "\n",
    "# Print the first few rows of the DataFrame\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Graphical exploratory data analysis\n",
    "<p>The bank has implemented different recovery strategies at different thresholds (\\$1000, \\$2000, \\$3000 and \\$5000) where the greater the Expected Recovery Amount, the more effort the bank puts into contacting the customer. Zeroing in on the first transition (between Level 0 and Level 1) means we are focused on the population with Expected Recovery Amounts between \\$0 and \\$2000 where the transition between Levels occurred at \\$1000. We know that the customers in Level 1 (expected recovery amounts between \\$1001 and \\$2000) received more attention from the bank and, by definition, they had higher Expected Recovery Amounts than the customers in Level 0 (between \\$1 and \\$1000).</p>\n",
    "<p>Here's a quick summary of the Levels and thresholds again:</p>\n",
    "<ul>\n",
    "<li>Level 0: Expected recovery amounts &gt;\\$0 and &lt;=\\$1000</li>\n",
    "<li>Level 1: Expected recovery amounts &gt;\\$1000 and &lt;=\\$2000</li>\n",
    "<li>The threshold of \\$1000 separates Level 0 from Level 1</li>\n",
    "</ul>\n",
    "<p>A key question is whether there are other factors besides Expected Recovery Amount that also varied systematically across the \\$1000 threshold. For example, does the customer age show a jump (discontinuity) at the \\$1000 threshold or does that age vary smoothly? We can examine this by first making a scatter plot of the age as a function of Expected Recovery Amount for a small window of Expected Recovery Amount, \\$0 to \\$2000. This range covers Levels 0 and 1.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXu8XVV1779DCO8kJIeAYEIT8sFIIhcw4aUolCih1WrAR6VqabWGezG9iKWExtvmRGt41Iq1qORIaCm1BRSsXu9HEGIBlQM1waBgeEtJECEJCMgzgXn/WGuFdXbWY67n3vvk9/181mfvvfacY4451zp7nDXmmGOacw4hhBAii9d0WwEhhBC9j4yFEEKIXGQshBBC5CJjIYQQIhcZCyGEELnIWAghhMilUWNhZnua2TfN7G4zW2tmR5vZRDO73szuC18nNKmDEEKI6jT9ZPEPwLXOuTcAhwBrgXOAlc65A4GV4WchhBA9jDW1KM/MxgF3AAe4WCNmdg9wnHPuUTPbF7jROTejESWEEELUwo4Nyj4A2AD8k5kdAqwGzgD2cc49ChAajL2TKpvZAmABwO677z77DW94Q4OqCiHE6GP16tUbnXOT6pDV5JPFHOBW4C3OudvM7B+Ap4E/d87tGSv3pHMuc95izpw5btWqVY3oKYQQoxUzW+2cm1OHrCbnLNYD651zt4Wfvwm8CXgsdD8Rvj7eoA5CCCFqoDFj4Zz7NbDOzKL5iLnAL4DvAKeG504Fvt2UDkIIIeqhyTkLgD8Hvm5mOwEPAn9KYKCuMrOPAQ8D729YByGEEBVp1Fg459YASf6yuVVlb968mfXr1/PCCy9UFdUz7LLLLkyePJkxY8Z0WxUhhBhB008WjbF+/XrGjh3L1KlTMbNuq1MZ5xybNm1i/fr1TJs2rdvqCCHECPo23ccLL7zAwMDAqDAUAGbGwMDAqHpSEkKMHvrWWACjxlBEjLb+CCFGD31tLIQQQrSDjIUQQohcZCyEEELkImNRgfnz5zN79mxmzZrF0NAQACtWrOD1r389xx13HB//+MdZuHAhABs2bOC9730vhx9+OIcffjg//vGPu6m6EEIUom9DZ3uBSy+9lIkTJ/L8889z+OGH8853vpPPfvaz3H777YwdO5bjjz+eQw45BIAzzjiDM888k2OOOYaHH36YefPmsXbt2i73QAgh/NiuniyG1w1z4r+eyPC64VrkfelLX+KQQw7hqKOOYt26dVx++eUce+yxTJw4kTFjxvD+97+6OP2GG25g4cKFHHroobz73e/m6aef5plnnqlFD9Esdd83ol10/ephu3qyWHrTUq574DoArv3wtZVk3Xjjjdxwww0MDw+z2267cdxxxzFjxozUp4VXXnmF4eFhdt1110rtivap874R7aPrVw/b1ZPFkmOXMG/6PJYcu6SyrKeeeooJEyaw2267cffdd3Prrbfy3HPPcdNNN/Hkk0+yZcsWrr766q3lTzjhBC666KKtn9esWVNZB9EOdd43on10/eqhsf0s6iRpP4u1a9dy0EEHdUkjePHFF5k/fz6PPPIIM2bMYMOGDQwODnLvvffy+c9/nv3224+DDjqIiRMn8rnPfY6NGzfyiU98grVr17Jlyxbe9ra3cfHFF28jt9v9EkKMHurcz2K7ckPVyc4778z3vve9bc7PmTOHBQsWsGXLFk466SROOOEEAPbaay+uvPLKttUUQoha2K7cUG0wODjIoYceyhvf+EamTZvG/Pnzu62SEEJURk8WNfP5z3++2yoIIUTt9PWTRT/MtxRhtPVHCDF66Ftjscsuu7Bp06ZR8wMb7Wexyy67dFsVIYTYhr51Q02ePJn169ezYcOGbqtSG9FOeUII0Wv0rbEYM2aMdpQTQoiW6Fs3lBBCiPaQsRBCCJGLjIUQQohcZCyEEH1NXVlllZ02m76d4BZCCKgvq6yy02YjYyGE6GuibLJVs8rWJWe00rdZZ4UQQmRTZ9ZZzVkIIYTIpVE3lJk9BDwDvAxscc7NMbOJwJXAVOAh4APOuSeb1EMIIUQ12niy+F3n3KGxR6FzgJXOuQOBleFnIYQQPUw33FDvAS4L318GaMMHIYTocZo2Fg74vpmtNrMF4bl9nHOPAoSvezesgxBCiIo0HTr7Fufcr8xsb+B6M7vbt2JoXBYA7L///k3pJ4QQwoNGnyycc78KXx8HvgUcATxmZvsChK+Pp9Qdcs7Ncc7NmTRpUpNqCiGEyKExY2Fmu5vZ2Og9cAJwJ/Ad4NSw2KnAt5vSQQghmmZ7SRPSpBtqH+BbZha182/OuWvN7CfAVWb2MeBh4P0N6iCEEI2yvaQJacxYOOceBA5JOL8JmNtUu0II0SbbS5oQ5YYSQogKHD3l6FH9RBGhdB9CiNI07a8fWj3EXhfsxdDqoUbkC39kLIQQpYn89UtvWtqI/MUrF7Pp+U0sXrm4EfnCHxkLIURplhy7hHnT5zXmr182dxkDuw6wbO6yRuQLf2QshBClifz1R085euu5Ol1TC2YvYOPZG1kwe0F+4T6mH8JvNcEthKiV7SWUtE76YcxkLIQQtbK9hJLWST+MmXbKE0KIUYp2yhNCCNEqMhZCCCFykbEQQgiRi4yFEEKIXGQsRN/TDzHq2zNVrk9b1zbeTq/cT72iR4RCZ0Xf0w8x6tszVa5PW9c23g7QE/dTr93XMhai7+mHGPXtmSrXp61rm9ROt++nnruvnXM9f8yePdsJUZZbHr7Fzbt8nrvl4Vu6KkM0Q5lr0+/X01d/YJWr6XdYTxZi1FPH43yvuQTEq5S5Nv1+Pbuhv4yFGPXU8Tjfcy4BsZUy16bfr2c39Fe6DyGEGKUo3YcQfUSvhUD6EOk8tHpohO5JfWm7f0Xaq6pbXX3rx3ugE7mhhGiYfvSPRzqv+tUqNj2/CQh0T+pL2/0r0l5V3erqWz/eA53IWAjRMP3oH490Pfmgk7lm7TXb9CEpxLSt/hVpr6pudfWtH++BTjRnIYQQoxTNWQjhyfC6YY665CiOuuSo1v3FPn7qusrUoVc3dKmCjy5Dq4fY64K9GFo91BP6lClbiZ3YvTZZdS3YaPLQojxRlnmXz3MM4hjEzbt8Xlfazmq3rjJ16NUNXargo8vA+QOOQdzA+QM9oU+ZslVggKecFuUJkc+SY5fwmxd+s/V9223ntVtXmTr06oYuVfDRZdncZSxeuZhlc5f1hD5lylbiGX5Vm6y6rE6Th54sRN20ne5h+arlbuD8Abd81fJW2qubsuPlW893fHzldZaLf47eL1+1vNZ7oO4xypPn0x41pvvouiHwOWQsRN207U5p0x3SBGXHy7ee7/j4yussF/8cvY/abNq9V7Zenjyf9uo0Fo27ocxsB2AV8Ihz7l1mNg24ApgI3A58xDn3UtN6CBGnbXdKm+6QJig7Xr71fMfHV15aqG+8XmdYcFXqHqM8eW3fw42HzprZp4A5wLjQWFwFXOOcu8LMLgbucM59NUuGQmeFEKI4fRM6a2aTgXcCl4SfDTge+GZY5DJgfpM6CNFJWyGi3djlLQmf8NF4iHFnio8mdW1qjNoM8e32Lnud7Y8IFa8xdLZpN9QXgbOBseHnAeA3zrkt4ef1wOuSKprZAmABwP7779+wmmJ7wif1Qj+lNc9rZ/HKxWx6fhOLVy5mwewFqTJue+Q2AO5/4v4RKT6a1LWpMWozvUa3d9nrbD+6jktvWgpj2a+udhozFmb2LuBx59xqMzsuOp1QNNEP5pwbAoYgcEM1oqTYLmkrRLSbu7zF8ZkPiIcYf/Swj9bqy8/StakxatOf3+1d9jrbj4eKX/fMdb0fOgucS/Dk8BDwa+A54OvARmDHsMzRwHV5shQNJZLo993O4vRLX3zCOY/82pHuyK8dOSJMtc1+ZYXNVpXlU66NPo+qnfKcc38F/BVA+GRxlnPuQ2b2DeB9BBFRpwLfbkoHMboZDZk8I/qlL3l6xt1ZS29aCnTXLVPV1eVbt21X1PayU94i4Aoz+1vgp8CKLuggRgG9tJq4Kv3SF59wzqQV8910y1QZ27KhumXbq1uvOlHWWSGEGKX0TeisECKdtndhayIcuJey0GbRVEbYtB0Fq8iuMqZNXg8lEhSiS7S9C1sT4cCjZa6latnOHQWryG5jjqUMMhZCdIm2d2FrIhx4tMy1VC3rmzqk6bDtJq+H5iyEEGKUojkLIXqUsjuz+c4FlPVJl90xsEkfeN19zKvb9vxKHbr01JxQXQs2mjy0KE/0C2VTkWel1M4qV1R+0bpNpnKvu495ddtOS1+HLlV1ph8W5QmxPVI2FbnvXECVNNhldgxs0gdedx/z6rY9v1KHLj01J1SX1Wny0JPF9kHTaRLqlN9L6Tl8Uk20pW+ZXeiq6tZU3zrHNZ7GpE796tA/bSdA7cEtRiVNh2HWKb+XQkZ9Uk20nQHXN5S0Dt3ayFwLI7O5FmnHJ0VKnSHNwNZrwE6MKyUwARkL0TM0/chdp/xecg/4pJpoOwNukV3oqurWVubaMm68JDlFvy/bxskHncxpF5z2dGmhHSh0VgghRikKnRWNUWeoXl2yonDURTcs2ib80zd9QlLYaFaqBp9Q06Rw1zJ1OndZS0sdUSYNxaIbFo0I5S2ic1J7UflZX5619SgajlukL/H20q6f71h39rdsio6i93X8/i2SGqQIZUOjC1HX5EeThya426PO8MK6ZEXhqGM+M2ab8E+fNtLCRqPzkfyk77Jkp4W7lqkz7/J5mfr49rWzbDRmUShvEZ2T2ouXjx9lQ3l9Q0ezrp/vWHfKSxvnPDlF7+vO+9e33SKkXUdqDJ3dYXBwsBkrVCNDQ0ODCxYkbwcp6mX6hOmsf3o9S45dwpTxU3pC1p677MnwumFOP+J0tryyhcnjJnPu3HOZMn6KVxvTJ0znzsfvHFEvrt/CIxay+eXNI2Sk1cnqX9k60efjph6Xqk/R8YzKzj9oPndvvJtlc5cxe7/ZhXROai8qP37n8UzabRKTdpvEARMOSO1vnn55fYm319lOnows/SePm8xZbz4rcZzz5BS9r+P37x477ZF6fauQdh2XLl366ODgYLEVoilozkIIIUYpmrMQtVMljURRv3FVispM8udmpeVI83V3zi3ktR9vI81XHp+nyJuP6ZRX1/xG2ng0mY68SiqMzuuZNYZ1kqZX2r1UdF6mLr0bSxFSlz+ryUNzFs1TNY1EEb9xVYrKTPLnZqXlyPJ1x+cW8tqPt5HmK++cp8jyjSfJq2N+I208qvrqs6iSCqPzemaNYZ2k6ZV2L5WZl6l7bNGchaibsvMLZfzGbeua5M+N/MiRLz9PftLcQl778TY+MOsDib7y+DzFhmc3ZM7HdMqra36jU9doPKr66rPIkuVzT8WvZ9YY1kmaXmn3UtF5mbr0jre74sIVtc1ZdP2pwefQk0V7+KQeKJPSoSjLVy13A+cPuOWrlo94H9chSr/gq0fRtApJKRTOvv5sN3D+gDv7+rO3GYNIz+i7pHay0kbEx3XmRTPd+HPHu+WrlqeOd5qsOq5P0pj7jFPS56L1i5zPKpM1zj7XIG2sl69anpn6o86/j6xrfOTXjnQzL5q59bukcWEn1rqafoc1wS1GcOK/nsh1D1zHvOnzUlMPRGUGdh1g0/ObMsuWZa8L9mLT85sY2HUAYOv7jWdvHKED4K2HT9/SykOQQmHMa8aw+ZXNW1/jbUcpLqLvktqJ6935fee4Rn2bs9+cxPFOk1XH9YmPfzTmPuMU16vMOCeNR9Z5ILMMpI9z1ndZYx2/PlnXuI6/j7xrHNE5FtF7/pGn3UY3vlTjHcgNJUbgG4qa5QKpg/ij/dwD5ia6SKqGQOaR5HqKQlGTwiAP3ufgESGSaS6WvFDVhUcs5MEnHuSll1/ivLefl+pyKhMSXGb8O910eeNaZZx9XGk+LsG8cfa5Bmljfdabz8p0d9X595F1jTtDiuNjEb1/4BsP/HLwrwe/XKrxDvRkIYQQoxSFzoquUmSHs6JhgU2F8CaVqRIa2hku6dP3Muk6onpFQyvT+uYbalsmdUkZfRbdsIjx545n1ldmJY6RTxs+1yJPVtXw16hsmTHL0zcp1Lro/VQLdU1+NHlogru3SAsJzEuvUCRNRRPpI/JCQYu03Rku6dP3Muk64vXKhAt39s031NanzbL9iX/uTOFS5pr4XIs8WVXDX+Nly6YgSdM3KdTa936ixtDZ/AKwD7AC+F74eSbwsboU8DlkLHqLItErRTeNKbsRTJEorrSonSJtd0YL+fS9iPyqG++k9S0tQiepfF6bZfsT/3z29We7ccvGuZlfnpk4Rj5t+FyLPFlp0V9F+1hmzPL0jV8zn0iwOG0bi+8BHwDuCD/vCPy8LgV8DhmL7lLnj2xUvsiP3/JVy91uf7ub2/1zu48IJfX9gywTDpz0o5r2x+vbl7xwYJ++pH2X1WefUM6sH9ykMM00/Yr+mNXR9zL1y7RdxAj5/jMRGcyk0PAy/es817ax+En4+tPYuTV1KeBzyFh0lzrdN/HyRV0+nauXfR/1i7io4u6INHdBklug6IrurBXkWX1J+y6rz3kuqE7d0ur7ZH4t6iapo+9l6pdpu4h7q4ibsjMzcJX+dZ6r01js6DGt8ayZDQAOwMyOAp7Kq2RmuwA3AzuHTyPfdM4tMbNpwBXAROB24CPOuZc89BBdonMXrrxXH3lFdh1bNncZZ157JmbGsrnLOHjvgzPrltGvc4e3pJ3e0naB8+3LsrnLWLxyMcvmLgMY8T5PryTZWTvkpe2clrZ7XaduSbvEPfPiM4zdeaxXe1m79vlStO9lyhSpmyYv6XzeeMbfH/LaQ1hx+4oR90LZ/jW5I2Ju6KyZvQn4R+CNwJ3AJOB9zrmf5dQzYHfn3G/NbAzwI+AM4FPANc65K8zsYgL31lezZCl0VgghitNq6Kxz7nbgWODNwGnArDxDEdZzzrnfhh/HhIcDjge+GZ6/DJhfQm+xnZEXctgZ+uiTRTXpfJEwxvh38Yyy8Z3kpn1xGnuet6eXHkntpIVN5snpJCljbV74ZtHQzLQxKiu7aqhzXrhwVubhImXy9C5Svko4ctL3DHCgt+J55PmpgJMTjrnA3h51dwDWAL8Fzgf2Au6PfT8FuDOl7gJgFbBq//33T/Xbie2DPP9sp4/YJ4tq0vm8Mnm76sXnV+KHrz86a+6kylxRUsbaPH982bmozjEqK7vqXFneXI3PvJHv3FKW3kXKF6nr1f99ca7FCe7/BzwBXB0em8Jz9xHMN/jI2BP4T+CtCcYiN7JKE9wiL/KjM/okLflgXjRLkTDG+HedSeZmXjTTzbxoppt64dStCQHz9EhqJy1sMk9OJ/Hx8I30KRPlljRGZWVXjcLLiwLziUgrklCxSN/SyhcN1c3rPwM85WoyFj5zFv8X+DPn3GPh532ArwJ/BtzsnHujzxOMmS0BngMWAa91zm0xs6OBQefcvKy6mrMQQojitJ3uY2pkKEIeB17vnHsC2JxWycwmmdme4ftdgbcDawmeMN4XFjsV+HYZxUV3qZpqoIr/Pc8H7pt2IZpbKJpGI2vHuzydkuYwOuXm+dZ951Wq+OrjspL6m9THtM9J4zy8bphZX56VOBa+FJkjKSKvzLxNXfMUPvMQZVPIVMUndPaHZvZd4Bvh5/cCN5vZ7sBvMurtC1xmZjsQGKWrnHPfNbNfAFeY2d8CPyVYHS76jKU3Ld2aIrlM+uXO+kXkJZWNnwO47ZHbtp5Pkrf0pqVby2SVy6p3/xP3b01VDXjrBEHY7ILZr2ZSjsvt1GXxysVsen7T1jppYxWdj1KlR6+d5TrlpfUzLquzv0l9TLqOaeO89Kal/GLjLxLHwpckHYveS0nysvpXpK5vW0Di+7z7dulNS3PL14mPsfgEwaT2MeHn/wL2dc49C/xuWiUXREwdlnD+QeCI4qqKXqJqPHeVtRo+8eV5ax+itR5Jawfy2o5kf/Swj6auxUjT6bHfPsaTLzy5zfqKrLUnPjH78c9Za0WS5KX1M5J16U8vze1v1mvSOC85dgnrnlrHI888krvWxEfHuF5l780q6xbKrDVKaytLTtp90sS6im3wmdgADgUuAB4icCMtrGvSxOfQBHfvUWWyNW1Cz1dG1qRj3uRxPG1F3k5mSX30yfWUlzokTa6vXr7c8vAtI3bcy+pj3u5vzvlP9ialB+lsJyl1SJKMrN3ofOT44Hs/pV3vMpPgaVRNkRKXsXzV8nYmuM3s9cAHgVMIIqCuBM5yzv1O8yZsJJrg7j2q7I6WttNZ/H2WjKxd3LL06NxdLG8ns7Q+wrY7k2Xttpa261qSXB+9fOnc3S1px7uk3fnS2vXdPS9rF7d4O1ltZe1ml9ZG2bHyvZ+iPnS2VWRXwTyK7jKYJWNg1wE2fXET7lfOKikVkuWGuhv4IfAHzrn7AczszDoaFf1PU24kHxlZbpS8FBFxl0iSWyVLlq8LIM8dlOeqydPLFx9XT5K7Ka1dH/dVVD/NxRe1k+f+y0sJU9aNWLRfPulLfMfFhzrSdcSv6WkXnPZ0ZaVCsp4sTiJ4sngzcC1BPqdLnHPT6mrcFz1ZCCFEcVoJnXXOfcs594fAG4AbgTOBfczsq2Z2Qh2Ni+5TJr1FU236tuuTtqOoPklhq3WGTSaFolYd8yopM7LaKBJaW8d9kTY2RcajrlDuptuqM6y3bt1yKTLBQZAp9jTgB3VNmvgcmuBujjLpLZpq07fdIqmii+hTNt1FkdQV8XQgVce8SsqMIunKy7btS9rYFBmPqvq01VbVe7ToPUGb+1n0wiFj0Rxl0ls01aZvu0U2ofHVJynyxldekdQVWVFHRce8SsqMIhshlW3bl7SxKTIeVfVpq60q92iZe6JOY5Gb7qMX0JyFEEIUp+10H0KMoE6f7dDqIcafO55ZX5k14lwRn7lvOoyIznTdnXMIvvUX3bDIO/VC9N2iGxZtreuTTqNTn6Q5iM60GVlpP8qmpWhi7iqrTZ90LVXmrbLusbw0InnpVNrAZ4zqRk8WojBVYsE760Yx6sA257Li1uNyonQPUfk8/eLy5+w3Z5t4ft/6Y14zhs2vbN5aF9LXiUQyozrRa3zdQdr6iqxY/6S1FFnrPLJ0zBrnutZ++JC3xiJeJmusqq7VSVsfk7Z+pk18xgjqfbLYYXBwsA45jTI0NDS4YEHx3DGiGaZPmM76p9ez5NglTBk/pVLdPXfZk5sfupnpE6fzhRO+sPXc8Lphls1dxuz9ZufKOXifg0eUz9MvLv8Dsz7AnY/fyeRxkzl37rlMGT/Fu/7pR5zOlle2bK173NTjUutFMucfNJ+7N97N6Ueczh477cHCIxay4dkNTB43mbPefBabX968Tf24Pp1tTJ8wneF1w7z08kuc9/bzRvR/4REL2fzy5q2vSfV9r1dcRtFrXpTpE6Zvc03S9Moaqyw9s+6xtD5njWvTY9KJzxgBLF269NHBwcFymRo7qWvyo8lDE9zVKDMR2qYOSROqRVMo+KQfWb5quRu3bJyb+sWphVNERKkzdv/c7lvr502Gdvah7MR+0nd1BSAkpSbpnID30bvKvdW0/DooO95F0uL4pBYpKq/OdB9dNwQ+h4xFNcqEWLapQ1KoZtEdytJCROPtde5iV6Sv8fDazl3x0uR09qFsyHDSd3WFNnfKSQrt9dG7yr3VtPw6KDvePvdlUtm0coXl1bhTntxQ2wE+j+ZVXEtVdUhyCfi4orLkJ7UXubz2HbsvMwZmZD6+J8kfXjfMc5uf43XjXseMgRmpbqO0fmWNQZa7J+m7yJ1U1RWS5FqJ3GI+brk67q2m5ddBWXecz32ZVDbNXVhU3oqvrHh68OzB8+oYA01wCyHEKEWhs6LnKbrrWFYoZBQaGoWb1rEjmk/Ia1Iairzd4sr0MU1G1dDMzhDhOnc2LHIuS2Z8XOM76lUJ9/XRP0+m7zWt0kadtNGWz+ZHQhSm6K5jWbu/RTuqXTh84dZQ1ao7ovnsNlZmt7g0WUX6H8lI2+nOl/iOePEQ4Tp2NixyLktm2q6DebsQVtU/T0/fa1qljTppoy0ZC9EIRXcdy9r9LUqzfdqc07jj13fUsiNa0VTjRXaLK9rHtO/SdrrzJZ46++C9D05t2wff61k0VX1818F46vK8XQir6p+np+81rdJGnbTSVl0z5U0eiobqPZJC9rJC/+rcAazuEMnO3FBpYZJ54bxp/Y9Cdmd+eWbpcNQy+ZrS+lXmmkShw1m77hXRNauvPt9VJe3+jcYrLUS2znxhRev5hhHHrzlKJCi6TVLIXlboXx3hjU2FSMbDYrPCUvPCedP6X0c4aplMsGn9KnNN4rLywpmLhD1XDZstS9b9mxUi66tTWd3rCiOOytVpLBQ6K0qRFLKXt9K4anhjUyGSnath08JS88J50/p/8D4Hb7NKvWi4aNFV7VEbSf0qc02SVoqnUSTsuWrYbFnS7t9ovHxW02fpVFb3OsKI49d8xYUralvBrdBZIYQYpSh0VjRO2RDQovhkmI23WyTTZpnwx856ad99+JoPs9Nnd2LRDYu26U9nFt1OPXwz1Xa2GWWs9Rkr3/6klasSaponK6lu3RlU6wh9rav9bjG8bhgGOLA2gXX5s5o8NGfRPm35ksv4t4v62bPSJ+TVS/vOBs0xiBvzmTGJ/UmbG8ibR8jSZ8xnxhRKgeLTn7RyRec3ssr71q1zfqLMta/zvm5yvqWQDjWm+1DorEikTJhrGeLhnT46JYW8ZpXvLFu2Xud3e+22F1fddRVnHn3miO+XzV3GX37/L5k8fnJqCKVP2G5Sm4e89hBW3L7Ca6x8+5NXrmyoaZGxL3pdfagj9LXu9ttmybFLuO6l656uS57mLIQQYpSiOQvRs5SdJygit866ZdKM1KFrnp++6DxFXalPyspPktOk3z7vmpbtU9061zFX1gvzH0BzcxbAFOA/gbXAXcAZ4fmJwPXAfeHrhDxZmrPoH8rOExSRW2fdtO/ifvRorsC37TLx8Gnf+85TVI3rryvNeN1ra8roU6VPdetcx1xZFZ2ocZ1Fk3MWW4C/cM7dbmZjgdVmdj3wJ8BK59x5ZnYOcA6wKEOO6CPKzhOUkVtH3TJpRuqHis7mAAAQCElEQVTQNc9PX3Seou7UJ0XlJ8lp0m+fd03L9qlunYumP0kq2wvzH0B70VDAt4F3APcA+4bn9gXuyaurJ4vRT1YqijKpE+rYna1Iyoe89oru/JfWXpHd1NJ23Tvya0e6mRfNLLxbYJI+ZSmaaqRz9z5fWVljlCQzr61OuWdff/Y217XqvVsn9Fu6D2Aq8DAwDvhNx3dPptRZAKwCVu2///41D6HoNaqEXvrIK/J99F2RlA957RXd+S+tPZ9xydI/7qJpwq3XhAxfl11ZN2Pa+JQNZ27KFVuGOo1F46GzZrYHcDXwSefc02bmVc85NwQMQRAN1ZyGoheoEnrpI6/I99G5tKyvZdwXRUKEfXTNGpcs/SMXTZTdtW63XhMyfF12Rd2MaSHMRVxYSW7Cplyx3abR0FkzGwN8F7jOOfeF8Nw9wHHOuUfNbF/gRufcjCw5Cp0VQoji9EXorAWPECuAtZGhCPkOcGr4/lSCuYyepWfC1lqm6X7XlX6jLb18QoKrpsjwORcPC412k+t8LZMKJb5LXV2hm1G9sjv+FU0hUmcobFMh4E3Qln5NuqHeAnwE+LmZrQnPLQbOA64ys48RzGO8v0EdKtPmble9RNP9LrqjWrd3HPPZOa3Kbmy+u88B2+wuF+2oV2RnvTSZ0XdZ/fIlaafBMvWjennXoOrOfWltJ7XVS7T199GYsXDO/QhIm6CY21S7ddMzYWst03S/6wgpbIIi4Yt5ZcrMm/ici+8ud83aa7bOTRTZWS9JZto8RpWwXCi/41/Reaw6Q2H7ad6hrb8PpfsQQohRSl/MWYjRRR2+3/i5vHQM8Xo+aczLpqEukpoj8umXTRVely5l01x01vVND99UqpWi+pdtr6j+bfW3Cl2ZR6krBrfJQ4vyuk/VdBJFY9nj9cpsKepbrmhqjiqpwuvUpUyai866vms/mkq1UlT/su0V1b+t/lbBtx36bVFe1UPGovvUuQG9zyrZeD2f1c+++mWtFM+qE1/1nLRqtwxldclarVxkpbPvqvIqq5DrXKldtr2i+rfV3yr4tlOnsdCchRBCjFI0ZyGEEKJVZCyEEELkImMhhBAiFxkLIYQQuchYCCGEyEXGQgghRC4yFkIIIXKRsRBCCJGLjIUQQohcZCyEEELkImMhhBAiFxkLIYQQuchYCCGEyEXGQgghRC4yFkIIIXKRsRBCCJGLjIUQQohcZCyEEELkImMhhBAiFxkLIYQQuchYCCGEyEXGQgghRC6NGQszu9TMHjezO2PnJprZ9WZ2X/g6oan2hRBC1EeTTxb/DJzYce4cYKVz7kBgZfhZCCFEj9OYsXDO3Qw80XH6PcBl4fvLgPlNtS+EEKI+2p6z2Mc59yhA+Lp3y+0LIYQoQc9OcJvZAjNbZWarNmzY0G11hBBiu6ZtY/GYme0LEL4+nlbQOTfknJvjnJszadKk1hQUQgixLW0bi+8Ap4bvTwW+3XL7QgghStBk6Oy/A8PADDNbb2YfA84D3mFm9wHvCD8LIYTocXZsSrBz7pSUr+Y21aYQQohm6NkJbiGEEL2DjIUQQohcZCyEEELkImMhhBAiFxkLIYQQuchYCCGEyEXGQgghRC4yFkIIIXKRsRBCCJGLjIUQQohcZCyEEELkImMhhBAiFxkLIYQQuchYCCGEyEXGQgghRC4yFkIIIXKRsRBCCJGLjIUQQohcZCyEEELkImMhhBAiFxkLIYQQuchYCCGEyEXGQgghRC4yFkIIIXKRsRBCCJGLjIUQQohcZCyEEELkImMhhBAil64YCzM70czuMbP7zeycbugghBDCn9aNhZntAHwZ+D1gJnCKmc1sWw8hhBD+dOPJ4gjgfufcg865l4ArgPd0QQ8hhBCe7NiFNl8HrIt9Xg8c2VnIzBYAC8KPL5rZnS3oVpW9gI3dVsKDftCzH3QE6Vk30rNeZtQlqBvGwhLOuW1OODcEDAGY2Srn3JymFauK9KyPftARpGfdSM96MbNVdcnqhhtqPTAl9nky8Ksu6CGEEMKTbhiLnwAHmtk0M9sJ+CDwnS7oIYQQwpPW3VDOuS1mthC4DtgBuNQ5d1dOtaHmNasF6Vkf/aAjSM+6kZ71Upue5tw20wVCCCHECLSCWwghRC4yFkIIIXLpaWPRS2lBzGyKmf2nma01s7vM7Izw/KCZPWJma8Lj92N1/irU/R4zm9eirg+Z2c9DfVaF5yaa2fVmdl/4OiE8b2b2pVDPn5nZm1rScUZszNaY2dNm9sleGE8zu9TMHo+v7SkzfmZ2alj+PjM7tSU9/87M7g51+ZaZ7Rmen2pmz8fG9eJYndnh/XJ/2Jek8PY6dSx8jZv+LUjR88qYjg+Z2ZrwfFfGMpSf9jvU/P3pnOvJg2Dy+wHgAGAn4A5gZhf12Rd4U/h+LHAvQbqSQeCshPIzQ513BqaFfdmhJV0fAvbqOHcBcE74/hzg/PD97wPfI1j/chRwW5eu9a+B3+mF8QTeBrwJuLPs+AETgQfD1wnh+wkt6HkCsGP4/vyYnlPj5Trk/BdwdNiH7wG/17COha5xG78FSXp2fP/3wN90cyxD+Wm/Q43fn738ZNFTaUGcc486524P3z8DrCVYjZ7Ge4ArnHMvOud+CdxP0Kdu8R7gsvD9ZcD82Pl/cQG3Anua2b4t6zYXeMA5998ZZVobT+fczcATCe0XGb95wPXOuSecc08C1wMnNq2nc+77zrkt4cdbCdYxpRLqOs45N+yCX5F/4dW+NaJjBmnXuPHfgiw9w6eDDwD/niWj6bEM9Uz7HWr8/uxlY5GUFiTrx7k1zGwqcBhwW3hqYfiId2n0+Ed39XfA981stQVpUwD2cc49CsENB+zdA3pGfJCRf4i9Np5QfPy6rS/ARwn+q4yYZmY/NbObzOyt4bnXhbpFtKVnkWvc7bF8K/CYc+6+2Lmuj2XH71Dj92cvGwuvtCBtY2Z7AFcDn3TOPQ18FZgOHAo8SvC4Ct3V/y3OuTcRZPb9hJm9LaNsV8fZgoWZ7wa+EZ7qxfHMIk2vbo/rp4EtwNfDU48C+zvnDgM+BfybmY2jO3oWvcbdvvanMPKfma6PZcLvUGrRFJ0K69rLxqLn0oKY2RiCC/R159w1AM65x5xzLzvnXgG+xquuka7p75z7Vfj6OPCtUKfHIvdS+Pp4t/UM+T3gdufcY9Cb4xlSdPy6pm84Wfku4EOhO4TQtbMpfL+aYA7g9aGecVdV43qWuMbdHMsdgZOBK6Nz3R7LpN8hWrg/e9lY9FRakNBvuQJY65z7Qux83L9/EhBFU3wH+KCZ7Wxm04ADCSa/mtZzdzMbG70nmPC8M9Qning4Ffh2TM8/DqMmjgKeih5nW2LEf229Np4xio7fdcAJZjYhdLOcEJ5rFDM7EVgEvNs591zs/CQL9pLBzA4gGL8HQ12fMbOjwnv8j2N9a0rHote4m78Fbwfuds5tdS91cyzTfodo4/6sc6a+7oNgJv9eAsv96S7rcgzBY9rPgDXh8fvA5cDPw/PfAfaN1fl0qPs91BwVkaHnAQTRIncAd0XjBgwAK4H7wteJ4Xkj2IzqgbAfc1oc092ATcD42LmujyeB8XoU2EzwH9jHyowfwZzB/eHxpy3peT+BLzq6Ry8Oy743vB/uAG4H/iAmZw7BD/YDwEWEmR0a1LHwNW76tyBJz/D8PwP/s6NsV8YylJ/2O9T4/al0H0IIIXLpZTeUEEKIHkHGQgghRC4yFkIIIXKRsRBCCJGLjIUQQohcZCxEIczsZRuZLbbxbMBmtqeZnV6i3qCZnZVyPsp6+gszO6UeTdvHzA4zM2ctZjVO0eNPzGy/buogmkXGQhTleefcobHjvBba3BMobCxyuNA5dyhBorXl4arYrhMuniryd3kK8KPwtZv8CSBjMYqRsRCVMbPxFuw1MCP8/O9m9vHw/W/N7O/N7HYzW2lmk8Lz083sWguSHf7QzN4Qnt/Hgn0Y7giPNwPnAdPDJ4G/C8v9pZn9xIJkdEtjunw61OUGYEae7i5IDvccQZrmonphZp8yszvD45PhufPjT0Lhk8xfpOltwf4Ia83sKwSLvP7azC6M1f+4mcVX60bnDXgfwQ/1CWa2S0ze3WZ2SajX183s7Wb2Ywv2LjgiLDfRzP4j1OVWM/sfMX3PirVzZygz0vNrFuyl8H0z29XM3kewGO3r4TXaNW/cRR9S90pIHaP7AF7m1ZWja4A/DM+/AxgmSMVwbay8I8hRBPA3wEXh+5XAgeH7I4EfhO+vJEiOBsE+BuPp2D+AIDXBEMHq1NcA3yXYj2A2wSrV3YBxBCtTk/ZNGIzOE+xh8MPYd0X0itrbHdiDYFXvYeFxU0zmL4D9M/SeCrwCHBWW351gxe2Y8PMtwMEJ/TgGWBm+/zfg5PD9VIIkggeH7awGLg3bfQ/wH2G5fwSWhO+PB9Z0jk/4+c5QZiT30PD8VcCHw/c30uLqfx3tHzsiRDGed4H7ZgTOuevN7P0EqQUOiX31Cq8mYftX4BoLMma+GfiGvbqR2M7h6/EEOXVwzr0MPGWvprCOOCE8fhp+3oMgP89Y4FsuzIlkZln5g84Mn34OIMzjX0KvY8L2ng3rXwO81Tn3JTPbO/ThTwKedM49bGb/O0Xvh4H/dsF+AzjnnjWzHwDvMrO1BEbj5wl9OIVgbwfC148AUWK5X0Z1zOwuAqPizOznBD/6EBib94Zt/sDMBsxsfMaYRXLXhO9Xx2SJUY6MhaiF0M9+EPA8we5b61OKOoL/dn+TZHR8mwPOdc4t79Dhk/inhL7QOfd5MzsZ+Bczm15Cr6wtM79J4CJ6La/+oKfpPRV4tqP+JcBi4G7gn7ZpOEhk917g3RakIzdgwMIkksCLseKvxD6/wqt/92lpqrcw0kW9S+x9XO7LgFxO2wmasxB1cSbBrl2nAJfGJoxfQ/CjCfBHwI9ckH//l+GTSDSpGz2NrAT+V3h+Bwv2CXiG4Kkh4jrgo+GTAGb2OjPbG7gZOCn0o48F/iBPaRekeF4FnFpCr5uB+Wa2mwUZfk8CfhiWv4LAJfc+AsORpXeSXrcRpJD+I5J3aHs7cIdzbopzbqpz7ncI0lYX2ZntZuBDoS7HARvDMXiIwD2HBXs2T/OQ1XmNxChDxkIUZVcbGTp7npm9Hvgz4C+ccz8k+BH6P2H5Z4FZZraawJXzmfD8h4CPmVmUHTfaJvMM4HdDd8lqYJYL9g74cTjR+nfOue8T+OiHw3LfBMa6YLvJKwnmUq7m1R/uPD4DfCp8Oiqi1+0EWUn/i2C3skuccz8FcM7dRfDj+Yh7dQezRL0z9LoK+LELtr3s5BSCvUriXE1gXHwZBOaY2c8IggiiFNdXAxPNbA2BgbzXQ9Y/Axdrgnv0oqyzolHM7LfOuT26rUc/YmbfJXCXrey2LkLoyUKIHsOCRYj3EgQTyFCInkBPFkIIIXLRk4UQQohcZCyEEELkImMhhBAiFxkLIYQQuchYCCGEyOX/A7gEMroEv105AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Scatter plot of Age vs. Expected Recovery Amount\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.scatter(x=df['expected_recovery_amount'], y=df['age'], c=\"g\", s=2)\n",
    "plt.xlim(0, 2000)\n",
    "plt.ylim(0, 60)\n",
    "plt.xlabel('Expected Recovery Amount')\n",
    "plt.ylabel('Age')\n",
    "plt.legend(loc=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Statistical test:  age vs. expected recovery amount\n",
    "<p>We want to convince ourselves that variables such as age and sex are similar above and below the \\$1000 Expected Recovery Amount threshold. This is important because we want to be able to conclude that differences in the actual recovery amount are due to the higher Recovery Strategy and not due to some other difference like age or sex.</p>\n",
    "<p>The scatter plot of age versus Expected Recovery Amount did not show an obvious jump around \\$1000.  We will be more confident in our conclusions if we do statistical analysis examining the average age of the customers just above and just below the threshold. We can start by exploring the range from \\$900 to \\$1100.</p>\n",
    "<p>For determining if there is a difference in the ages just above and just below the threshold, we will use the Kruskal-Wallis test which is a statistical test that makes no distributional assumptions.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KruskalResult(statistic=3.4572342749517513, pvalue=0.06297556896097407)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Import stats module\n",
    "from scipy import stats\n",
    "\n",
    "# Compute average age just below and above the threshold\n",
    "era_900_1100 = df.loc[(df['expected_recovery_amount']<1100) & \n",
    "                      (df['expected_recovery_amount']>=900)]\n",
    "by_recovery_strategy = era_900_1100.groupby(['recovery_strategy'])\n",
    "by_recovery_strategy['age'].describe().unstack()\n",
    "\n",
    "# Perform Kruskal-Wallis test \n",
    "Level_0_age = era_900_1100.loc[df['recovery_strategy']==\"Level 0 Recovery\"]['age']\n",
    "Level_1_age = era_900_1100.loc[df['recovery_strategy']==\"Level 1 Recovery\"]['age']\n",
    "stats.kruskal(Level_0_age,Level_1_age) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Statistical test:  sex vs. expected recovery amount\n",
    "<p>We were able to convince ourselves that there is no major jump in the average customer age just above and just \n",
    "below the \\$1000 threshold by doing a statistical test as well as exploring it graphically with a scatter plot.  </p>\n",
    "<p>We want to also test that the percentage of customers that are male does not jump as well across the \\$1000 threshold. We can start by exploring the range of \\$900 to \\$1100 and later adjust this range.</p>\n",
    "<p>We can examine this question statistically by developing cross-tabs as well as doing chi-square tests of the percentage of customers that are male vs. female.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sex                Female  Male\n",
      "recovery_strategy              \n",
      "Level 0 Recovery      108   139\n",
      "Level 1 Recovery      316   354\n",
      "0.3941650543686612\n"
     ]
    }
   ],
   "source": [
    "# Number of customers in each category\n",
    "crosstab = pd.crosstab(df.loc[(df['expected_recovery_amount']<2000) & \n",
    "                              (df['expected_recovery_amount']>=0)]['recovery_strategy'], \n",
    "                       df['sex'])\n",
    "print(crosstab)\n",
    "\n",
    "# Chi-square test\n",
    "chi2_stat, p_val, dof, ex = stats.chi2_contingency(crosstab)\n",
    "print(p_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Exploratory graphical analysis: recovery amount\n",
    "<p>We are now reasonably confident that customers just above and just below the \\$1000 threshold are, on average, similar in terms of their average age and the percentage that are male.  </p>\n",
    "<p>It is now time to focus on the key outcome of interest, the actual recovery amount.</p>\n",
    "<p>A first step in examining the relationship between the actual recovery amount and the expected recovery amount is to develop a scatter plot where we want to focus our attention at the range just below and just above the threshold. Specifically, we will develop a scatter plot of Expected Recovery Amount (Y) vs. Actual Recovery Amount (X) for Expected Recovery Amounts between \\$900 to \\$1100.  This range covers Levels 0 and 1.  A key question is whether or not we see a discontinuity (jump) around the \\$1000 threshold.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEKCAYAAAAvlUMdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYFOW59/HvDySyCAiIiiwBiYnLARFGAY2KEoF4JCrGKBpFY0KSE49b1IPxGNBEs5jENx49KhoElWiM+8mVSJAImgDqjCGoKBFXRgkgCm64APf7R1UPxdDT3TNT1dXL/bmuuqb7meqqu2t66u5nqadkZjjnnHNxa5N2AM455yqTJxjnnHOJ8ATjnHMuEZ5gnHPOJcITjHPOuUR4gnHOOZeIxBKMpL6SHpX0vKTnJJ0blneXNFfSi+HPbmG5JF0raYWkpZKGRrY1KVz/RUmTkorZOedcfJTUdTCSegG9zOxpSZ2BOuA44AzgbTP7qaQpQDcz+y9JRwP/CRwNDAd+bWbDJXUHaoEawMLtDDOzdxIJ3DnnXCwSq8GY2Sozezp8/B7wPNAbOBaYFa42iyDpEJbfZoHFwM5hkhoLzDWzt8OkMhcYl1Tczjnn4rFDMXYiqT9wAPAEsJuZrYIgCUnaNVytN7Ay8rL6sKyp8mz7mQxMBujUqdOwvffeO7434ZxzFa6uru4tM+sZ1/YSTzCSdgLuBc4zs3clNblqljLLUb59odl0YDpATU2N1dbWNj9g55yrUpJei3N7iY4ik9SOILnMNrP7wuLVYdNXpp9mTVheD/SNvLwP8GaOcueccyUsyVFkAn4DPG9mv4r86iEgMxJsEvBgpPz0cDTZCGBD2JQ2BxgjqVs44mxMWOacc66EJdlEdghwGvCMpCVh2Q+AnwJ3SzoLeB04MfzdHwlGkK0APgTOBDCztyX9CHgqXO8KM3s7wbidc87FILFhymnL1gfz6aefUl9fz0cffZRSVK4atW/fnj59+tCuXbu0Q3EuJ0l1ZlYT1/aKMoqsVNTX19O5c2f69+9PjsEGzsXGzFi3bh319fUMGDAg7XCcK6qqmirmo48+okePHp5cXNFIokePHl5rdlWpqhIM4MnFFZ1/5ly1qroE45xzrjg8wTjnnEuEJ5gSNX/+fBYuXNiqbey0004xRePyueqqq9IOwbmS4wmmRMWRYAq1efPmouwnn02bNqUdQot5gnFue55giuy4445j2LBh7LfffkyfPh2Ahx9+mKFDh7L//vszevRoXn31VW688UauueYahgwZwuOPP84ZZ5zBPffc07CdTO3k/fffZ/To0QwdOpRBgwbx4IMPZt1vY/Pnz+eII47glFNOYdCgQQDccccdHHTQQQwZMoRvf/vbDYmncXwAb7/9NscddxyDBw9mxIgRLF26lC1bttC/f3/Wr1/fsJ/Pfe5zrF69mrVr13LCCSdw4IEHcuCBB/K3v/0NgGnTpjF58mTGjBnD6aefzqGHHsqSJUsaXn/IIYewdOnSrO/hySef5OCDD+aAAw7g4IMPZvny5QDMnDmT4447jvHjxzNgwACuu+46fvWrX3HAAQcwYsQI3n47uE53yZIljBgxgsGDB3P88cfzzjvBHSBGjRpF5hqqt956i/79+zdsd8KECYwbN4699tqLiy++GIApU6awceNGhgwZwqmnnlrQ8XeuKphZRS7Dhg2zxpYtW7ZdWT4LX19oY28fawtfX9js12azbt06MzP78MMPbb/99rN//etf1qdPH3v55Ze3+f3UqVPt6quvbnjdpEmT7Pe//33D806dOpmZ2aeffmobNmwwM7O1a9fawIEDbcuWLdusk82jjz5qHTt2bNjvsmXL7JhjjrFPPvnEzMy++93v2qxZs2zNmjVZ4zv77LNt2rRpZmY2b94823///c3M7JxzzrEZM2aYmdnixYtt9OjRZmY2ceJEe/zxx83M7LXXXrO999674X0OHTrUPvzwQzMzmzlzpp177rlmZrZ8+XLL9nfM2LBhg3366admZjZ37lybMGGCmZndeuutNnDgQHv33XdtzZo11qVLF7vhhhvMzOy8886za665xszMBg0aZPPnzzczs8suu6xhv4cffrg99dRTDcf0s5/9bMN2BwwYYOvXr7eNGzdav3797PXXX897rDPH17lSB9RajOfhqrrQsiUuX3A5c14Kpj57+OsPt3p71157Lffffz8AK1euZPr06Rx22GENF+F17969WdszM37wgx/w2GOP0aZNG9544w1Wr17N7rvvnve1Bx10UMN+582bR11dHQceeCAAGzduZNddd2Xx4sVZ4/vrX//KvffeC8CRRx7JunXr2LBhAyeddBJXXHEFZ555JnfddRcnnXQSAI888gjLli1r2Pe7777Le++9B8BXvvIVOnToAMCJJ57Ij370I66++mpmzJjBGWec0WT8GzZsYNKkSbz44otI4tNPP2343RFHHEHnzp3p3LkzXbt2Zfz48QAMGjSIpUuXsmHDBtavX8/hhx8OwKRJkzjxxBOz7idq9OjRdO3aFYB9992X1157jb59++Z5lXPVyRNMHlMPn7rNz9aYP38+jzzyCIsWLaJjx46MGjWK/fffv6FpJ5cddtiBLVu2AEFS+eSTTwCYPXs2a9eupa6ujnbt2tG/f/+CL+rr1KlTw2MzY9KkSfzkJz/ZZp2HHnoo63UclmWKIUmMHDmSFStWsHbtWh544AH++7//G4AtW7awaNGihkTSVBwdO3bkqKOO4sEHH+Tuu+8m1y0XLrvsMo444gjuv/9+Xn31VUaNGtXwux133LHhcZs2bRqet2nTJm9fT/RYNz6W0e22bdu2rPuNnEua98HkMbLvSB7++sOM7Duy1dvasGED3bp1o2PHjrzwwgssXryYjz/+mAULFvDKK68ANPQPdO7cueEbPkD//v2pq6sD4MEHH2z4tr5hwwZ23XVX2rVrx6OPPsprr7Xsdg6jR4/mnnvuYc2aNQ1xvPbaa4wcOTJrfIcddhizZ88GgsS5yy670KVLFyRx/PHHc8EFF7DPPvvQo0cPAMaMGcN1113XsL9oP0tj3/zmNznnnHM48MADc9boNmzYQO/ewb3nZs6c2az327VrV7p168bjjz8OwO23395Qm4ke62i/Vy7t2rXbpgblXDlYtHIR4+4Yx6KVixLZvieYIho3bhybNm1i8ODBXHbZZYwYMYKePXsyffp0JkyYwP7779/QpDR+/Hjuv//+hk7+b33rWyxYsICDDjqIJ554ouFb/6mnnkptbS01NTXMnj2blt7Fc9999+XHP/4xY8aMYfDgwRx11FGsWrWqyfimTZtGbW0tgwcPZsqUKcyaNathWyeddBJ33HFHw7oQNA1m1t9333258cYbm4xl2LBhdOnShTPPPDNnzBdffDGXXHIJhxxySItGws2aNYuLLrqIwYMHs2TJEn74wx8CcOGFF3LDDTdw8MEH89ZbbxW0rcmTJzN48GDv5HdlJdMFcPmCyxPZflXNpvz888+zzz77pBSRK9Sbb77JqFGjeOGFF2jTpjK+A/lnz5WiRSsXcfmCy5l6+FRG9h0Z+2zKlfHf6yrGbbfdxvDhw7nyyisrJrk4V6ri7ALIxjv5K9wzzzzDaaedtk3ZjjvuyBNPPJFSRLmdfvrpnH766duU3Xrrrfz617/epuyQQw7h+uuvL2ZozrlmqroEY2ZVNbvtoEGDcnaol4Mzzzwzb39MKavUZmjn8qmqNoj27duzbt06/4d3RWPhDcfat2+fdijOFV1iNRhJM4BjgDVm9m9h2e+AL4Sr7AysN7MhkvoDzwOZC0IWm9l3wtcMA2YCHYA/AudaCzNEnz59qK+vZ+3atS16T861ROaWyc5VmySbyGYC1wG3ZQrMrGHcqqRfAhsi679kZkOybOcGYDKwmCDBjAP+1JKA2rVr57etdc65IkmsiczMHgPezvY7BZ0gXwPuzLUNSb2ALma2KKy13AYcF3eszjnn4pdWH8yhwGozezFSNkDS3yUtkHRoWNYbqI+sUx+WOeecK3FpjSKbyLa1l1VAPzNbF/a5PCBpPyDbcK8m+18kTSZoTqNfv34xhuucc665il6DkbQDMAH4XabMzD42s3Xh4zrgJeDzBDWWaO9oH+DNprZtZtPNrMbManr27JlE+M7FIuk5oJwrBWk0kX0JeMHMGpq+JPWU1DZ8vCewF/Cyma0C3pM0Iuy3OR0o7I5azpWwpOeAcq4UJDlM+U5gFLCLpHpgqpn9BjiZ7Tv3DwOukLQJ2Ax8x8wyAwS+y9Zhyn+ihSPInCslcd4GwrlSVVWTXTrnykfjiRhd8uKe7LLqpopxzpWHuO8m64rPE4xzriR5M2L58wTjnCtJmankXfmqqskunSslPlTZVTqvwTiXEu9jcJXOE4xzKfE+BlfpPME4lxLvY3CVzvtgnHPOJcITjHPOuUR4gnHOOZcITzDOOecS4QnGOedcIjzBOOecS4QnGOecc4nwBOOccy4RnmCcc84lIm+CkXRIIWXOOedcVCE1mP8psMw555xr0ORcZJJGAgcDPSVdEPlVF6Bt0oE555wrb7lqMJ8BdiJIQp0jy7vAV/NtWNIMSWskPRspmybpDUlLwuXoyO8ukbRC0nJJYyPl48KyFZKmNP8tOuecS0OTNRgzWwAskDTTzF5rwbZnAtcBtzUqv8bMfhEtkLQvcDKwH7AH8Iikz4e/vh44CqgHnpL0kJkta0E8zjnniqiQ6fp3lDQd6B9d38yOzPUiM3tMUv8C4zgWuMvMPgZekbQCOCj83QozexlA0l3hup5gnHOuxBWSYH4P3AjcAmyOYZ9nSzodqAW+b2bvAL2BxZF16sMygJWNyoc3tWFJk4HJAP369YshVOeccy1VyCiyTWZ2g5k9aWZ1maWF+7sBGAgMAVYBvwzLlWVdy1GelZlNN7MaM6vp2bNnC0N0zjkXh0ISzP9J+g9JvSR1zywt2ZmZrTazzWa2BbiZrc1g9UDfyKp9gDdzlDvnnCtxhTSRTQp/XhQpM2DP5u5MUi8zWxU+PR7IjDB7CPitpF8RdPLvBTxJUIPZS9IA4A2CgQCnNHe/zjnnii9vgjGzAS3ZsKQ7gVHALpLqganAKElDCBLUq8C3w308J+lugs77TcD3zGxzuJ2zgTkE197MMLPnWhKPc8654pJZk10awQpBh/x2zKzx8OOSUlNTY7W1tWmH4ZxzZUNSnZnVxLW9QprIDow8bg+MBp5m++tbnHPOuQaFNJH9Z/S5pK7A7YlF5JxzriK0ZLr+Dwk64Z1zzrkm5a3BSPo/tl570hbYB7g7yaCcc86Vv0L6YKLzhm0CXjOz+oTicc45VyHyNpGFk16+QDCTcjfgk6SDcs45V/4KuaPl1wguejwR+BrwhKS80/U755yrboU0kV0KHGhmawAk9QQeAe5JMjDnnHPlrZBRZG0yySW0rsDXOeecq2KF1GAeljQHuDN8fhLwp+RCcs45VwkKudDyIkknAIcQTD453czuTzwy55xzZa2QGgxmdq+kuZn1JXU3s7cTjcw551xZK+RCy28DVwAbgS0EtZgWTdfvnHOuehRSg7kQ2M/M3ko6GOecc5WjkNFgLxHMP+acc84VrJAazCXAQklPAB9nCs3snMSics45V/YKSTA3AX8BniHog3HOOefyKqSJbJOZXWBmt5rZrMySeGTOOee2sWjlIsbdMY5FKxclsl0+Q6c4t1tIgnlU0mRJvSR1zyz5XiRphqQ1kp6NlF0t6QVJSyXdL2nnsLy/pI2SloTLjZHXDJP0jKQVkq6VpBa9U+ecK3OXL7icOS/N4fIFlyeyXTqzR5zbLaSJ7JTw5yWRskKGKc8ErmPbWyvPBS4xs02SfhZu87/C371kZkOybOcGYDKwGPgjMA6fScA5V4WmHj51m59xb3fOe3PejHO7hUzXP6DxAuxdwOseA95uVPZnM9sUPl0M9Mm1DUm9gC5mtsjMjCBZHZdv384515SkmpmKYWTfkTz89YcZ2XdkItvlEz6Ic7sFT1qpwJGSbgFWxrDvb7BtTWSApL9LWiDp0LCsNxC9uVl9WNZUjJMl1UqqXbt2bQwhOucqTVLNTG57hVzJP5ygmex4oDvwPeCi1uxU0qUEd8ecHRatAvqZ2TpJw4AHJO1HMGtAY5alLPiF2XRgOkBNTU2T6znnqldSzUxue03WYCRdKelF4CqCIcoHAGvDUWTvtHSHkiYBxwCnhs1emNnHZrYufFxHcHHn5wlqLNFmtD5ArG2EzsWtnJtgyjn2QiXVzOS2l6uJbDKwmqCT/Y4wAbSqViBpHEGn/lfM7MNIeU9JbcPHewJ7AS+b2SrgPUkjwtFjpwMPtiYGV52KeeIs5yaYco7dlZ5cTWS7A2OAicD/k/Qo0EHSDpGO+iZJuhMYBewiqR6YSjBqbEdgbjjaeLGZfQc4DLhC0iZgM/CdyGzN3yUYkdaBoM/GR5C5ZmsYhglBZ2aCyrkJJo7YF61cxOULLmfq4VO9llDlFLZS5V5Jak/QrDUR+CIwz8xOyf2qdNXU1FhtbW3aYbgS4Se94hl3xzjmvDSHsQPHJp7MXbwk1ZlZTVzbK/R+MB8B9wD3SOpC0OHvXNloGIbpElfONTgXr4KHKWeY2bs+VYyrZNXQ0Z0k70R3GQXVYJyrJsXsr3GukjW7BuNcpZt6+FTGDhzrTTwucZVeWy7kQsta4Fbgt625/sW5cuH9Na4YFq1cxPg7x7Nu4zqgMmvLhdRgTgb2AJ6SdJeksT6jsXMuSdFv9pX2LT/zfs6fcz7rNq6jR4cesdWWS+1Y5a3BmNkK4FJJlxEMVZ4BbJE0A/h15HoV55yLRbQfDKioPrHMexvee3hDU2xcAyJKrf+woE5+SYOBM4GjgXsJ5hD7IsGdLrNNse+cy8KvxylMtqHOldInFn1vcX8GSm2IeN4LLSXVAeuB3wD3mtnHkd/dZ2YTkg2xZfxCS1eK/CJEV8rivtAyZx+MpDYESWW0mf02mlwASjW5OFeqfIRa8xWzX6HU+jDKXc4mMjPbEk5QeVWR4nGuovkIteYrZr9CqfVhlLtCRpHNlXShpL6SumeWxCNzzjmKW+trzr68tpNfIX0wr2QpNjPbM5mQ4uF9MM65JLW2Py2tAR+59lvUPhgAMxuQZSnp5OKcq0zZag1p1SRaW7NK6947xdxvIVfydwQuILil8WRJewFfMLM/JB6dc85FZOsjSavfJFd/WiG1k7SGFBdzv4VcB3MrUAccHD6vB34PeIJxzjUoRpNPrutjSmlkXiFJL60BH8XcbyEJZqCZnSRpIoCZbfSpYpxzjRWjJpHt5FiKI/NKMemloZBRZJ9I6gAYgKSBwMe5X+Jc6/gInfLj1/hslfY9cUrl/6eQBDMNeBjoK2k2MA+4uJCNS5ohaY2kZyNl3SXNlfRi+LNbWC5J10paIWmppKGR10wK139R0qTmvEFXntLqAHUtl/ZJ1W1VKv8/hYwi+zMwATgDuBOoMbP5BW5/JjCuUdkUYJ6Z7UWQrKaE5V8G9gqXycANECQkYCowHDgImJpJSq5y+bfh8lAq35TLSTGOWcn8/5hZzgV4CJgIdMq3bhOv7w88G3m+HOgVPu4FLA8f3wRMbLxeuO+bIuXbrNfUMmzYMHOumix8faGNvX2sLXx9YdH2Ofb2scY0bOztY4u2T7Ot7/Wm2puK/p5bK61jVgig1lpwnm9qKaSJ7JfAocAySb+X9FVJ7VuR03Yzs1VhclsF7BqW9wZWRtarD8uaKt+OpMmSaiXVrl27thUhOld+0mgWyfVNOclv6pn3+oN5PyiJpqDmiKt2UQ61x0KayBaY2X8AewLTga8BaxKIJdvINMtRvn2h2XQzqzGzmp49e8YanHOlrvGJqxgnoFz9LkkmvMx7vWr0VUVpCorjWGa2AcTSV1Uq/Sy5FHo/mA7AeOAkYCgwqxX7XC2pl5mtktSLrcmqHugbWa8P8GZYPqpR+fxW7N+5itR4uG7aEze2dqhurutqou918rDJrQu0AHEcy7j/HmUxFDpfGxrwO+BV4EbgSKBNc9rg2L4P5mpgSvh4CvDz8PG/A38iqLGMAJ4My7sDrwDdwuUVoHu+/VZrH0wa7fCuNJX7ZyHuvorWHI84jmXa+y8EMffBFJIgxgFtW7TxYNTZKuBTgprIWUAPgtFjL4Y/u4frCrgeeAl4hmC0WmY73wBWhMuZhey7WhNMKXcgOtccuU6qLTnhtuZ/o7Un+Na+vlj/12kkmHbAOcA94fKfQLs4g0hiqdYEU+7fWp0rREtOuIX8byx8faENv3m4Db95+DbrtfYE39rXV3IN5haCPpcjw+VW4JY4g0hiqdYE41rPk3TpS+pvlEkEjZNBHDWYbImr1MSdYAq5H8w/zGz/fGWlxu8H4wqRrSO5tff5cOVr0cpFnD/nfACuGXtNrLMSlMPnKu77wRQyimyzpIFm9lIYwJ7A5rgCcC5N2Ub2lMXoHJeIkX1HsvibixPZdjV+rgq50PIi4FFJ8yUtAP4CfD/ZsJwrjmwXvfmcWukrh4sIG8sXcxyfq3I7LnlrMGY2L3OTMYKRXi+Ymc+m7CpCKU717tK/hqclihFzuR2XQu5o+T1gtpktDZ93k3SWmf1v4tG1wgeffMC4O8YV/X7XzrnWK8fmpGLEXHbHJd8oAGBJlrK/xznSIImly4Aufj1IDHxEVfWqpr99Nb3XqMbvmxQmu2wTvYOlpLbAZxLJdjHaY6c98s5RVG7tmUnLdjzKYb4jl19LPuu5/vaV9r9TrZ/zpN93IaPI5gB3S7qRYJLJ7xDcgKykdfpMp7xtlOXWnpk0H1FVuVryWc/1t6+0/520P+e55l1LUuLvO18Vh2Ck2XcJruK/F/g2LZw6pphLIRdaVmu1uCmVfDwq+b0VIvr+k55Xq9jHuhL+tqUyxRPFvpI/2CcdgC/EueOkF7+S30WVyj9w3JKak6s1J+2kpnFpzf5KPQmVSnxFTzDAVwjuLvlK+HwI8FCcQSSxeIKpLvn+QUvlHzhuSZ3Miz0xZNL7q9QvGHFLI8HUAV2JjBwDlsYZRBKLJ5jqUq0nkKQSa6U1cxW6/WzrVeqXk2zSSDBPhD89wbiSVU0ngeao1sTbUtmOVzGPYdqf47gTTCHDlJ+VdArQVtJekv4HWNiqkQWuaCptOGlTfHqX7OK6/3u1yHa8inkMWztsuOT+3/NlIKAjcCXwVLhcCbSPM8slsVRKDaZcblTk0pX2N18Xj7T/3yl2DcbMPjSzS83swHC5FNgtwZxXdpL81tDabzRTD5/K8N7DWf/R+tL5VlNBSuUbY9oXCpbKcSh3ra2Jl1yNNVf2AUYCXwV2DZ8PBn4LrIwzyyWxFLMGk2QtIY5vpl6LSU7Sx7Y1ndNxrFso/4xVBorVyQ9cDTwP3EnQNDYVWA2cSyuayAhmZV4SWd4FzgOmAW9Eyo+OvOYSYAXBcOmxhewnzgRT7kNgG8dX6vGWk6SPZRIn7iS26Z+prcr5WBQzwSzLJBKgG7AR2CvWnUNb4F/AZ8MEc2GWdfYF/gHsCAwAXqKAmQTiTDCV9u2s0t5PPmn/w7dm/0nEnvbxiEMpv4e4/7+K+V7jTjC55iLbaGYfhc1o70habmYvNqf5rQCjgZfM7LXIfJqNHQvcZcE9aF6RtAI4CChaY+/Uw6ey/qP1Df0YpTxSqZA5jdKed6nY0p43qzX7T+J+NZVwD5y0/6a5xP3/VcrvNa+mMg+wHngosmzzPI7sBswAzg4fTwNeBZaG5d3C8uuAr0de8xvgq01sbzJQC9T269cvvrRu5fOtv1ziLKa0v+2mvf8kpfXeKvmYNlbONZhcJ//Dcy2t3nEw5f9bwG7h890ImszaEAyFnhGWX58lwZyQb/txd/KXywe6FOMsxZhKVbkdK/9CU1mKlmCSXgiavv7cxO/6A8+Gjy8BLon8bg4wMt/2K+U6mEpQ7JNQuZ2ko0rphF3IcSznY13JWvp3qaQEcxdwZuR5r8jj8wn6XQD2Y9tO/pcpcid/tUjqZJH09PCNldJJurlKae6wYh1HT1Lxa+nfriISDMHsAOuArpGy24Fnwj6YhxolnEsJRo8tB75cyD48wTRf2ifmXCea5sQW5wmrXE5+5Tz0OO3PXVxK6bNS9TWYpBdPMM2X9j9IrhNNWrGVy8kv7b9da5R67IXGVy6flVyK2cn/f41GkW2zxBlEEosnmPJTiieaUozJJSfb37vQxFEJn5W4E4yCbW5P0uFZfxEyswW5fp+2mpoaq62tTTsM58pCWveEL7V4xt0xjjkvzWHswLEN15yU2rFJkqQ6M6uJa3tNXmhZ6gmkuarpQ+Jcc5XaxXxpxZPtIslCL0z1c8z28s6mHN4D5h5JyyS9nFmKEVyc0p5t1iU742502z6zb/Nlm4U3zeOY1qzArZnNOMlzTDH+FotWLoIe7BXnNnNNFZNxK8FEl9cARwBnAk3O61Kqqm16lFKU5LfS6LaBkvo2Xg6yfUtPs1ZTjtPZJHmOKcbf4vIFl8Nn6BLrRvN10gB14c9nImWPx9kRlMSSRid/JXTyJSnJ4xPdtv8d4uHHsXQU42+x8PWFRg82WDE6+TMk/Q04FLgH+AvBlPo/NbMvxJrpYpZGJ3+2DkK3LW+ndmnzz2DT4u7kz9sHQ3Cvlo7AOcAw4DRgUlwBlIo42jhL7m5yJcj7wlzayu0zWNZ9inFWh0ppaW4TWVoXSVVbM0S1vd+m+HFIT7kd+2Kemyji/WAAkPQosF07mpkdmUTCS0tagwBKbXho0prbeVupzRnV9ncvJeU2gKCcBygVMorswsjj9sAJwKZkwklPWh+6cv7wFENSJ+K0E5f/3V2hyi0hbqMl1R5gQZzVqCSWXE1k5VZFLgWVdmOpSpg3yrm4kUITWffI0zYEHf27J5PuiqNSvxUnKXPM1n+0np3b71y095jUtzevQTiXvEKayOoI+mBE0DT2CnBWkkElLamTSyW3q2eO1fqP1lfEeyzrZgfnykQh18G0N7OPGpXtaGYfJxpZK3Xds6s9vKBlUz7kkquWUsk1mIxqeI/OVau4r4MpJME8bWZD85WVGu0hG/vz+C949IspnXPy1sAzAAASA0lEQVSVqmgXWkraXdIwoIOkAyQNDZdRBBdelrS2bdoyYZ8JsW8338WU0+ums8vPd2F63fTY952mSn1fzrnk5LofzCTgDKAGeIqtE1y+C8wys/uKEWBLJVWDyWeXn+/Cuo3r6NGhB29d/FZR952ETJPY4vrFbPh4Q1HelzfDOZeOotVgzGyWmR0BnGFmR5rZEeFybBzJRdKrkp6RtERSbVjWXdJcSS+GP7uF5ZJ0raQVkpZKyts816V9l1Z34rdkioarRl9Fjw49uGr0Va3ad6nIDFzo3bl30d5XuUzlsWjlIkbcMoIRt4woz2k8nEtavnHMwFXAzpHn3YAft3Z8NPAqsEujsp8DU8LHU4CfhY+PBv5EUIsaATyRb/txzKbs10qkc/1LKV2nlCuWzOej2j8jrnIQ83UwhSSCv2cpe7rVO86eYJYDvcLHvYDl4eObgInZ1mtq6TKgS6tPUK050ZXSSbIUlcuxzSSR4TcP3+52ADfV3mTDbx5uw28e7n9nVxHSSDBLgR0jzzsAz7V6x8H1NE8TXGczOSxb32idd8KffwC+GCmfB9Rk2eZkoBaopWu63yrTqP2UU1LLdXyaeh+Z8uE3Dy/asW28zx4/61HU/ZeCcvpcudaJO8EUcqHlHcA8SbcSXHD5DeC2Al6XzyFm9qakXYG5kl7IsW62O2hmm4BzOjAdgutg0rxKO40rxcvpQs9cx6ep95EpH957eNFui5C5IHPRykWMv3M86zau43N8rqpuy1BOnytXYgrJQsA44BfAL4GxcWa4cPvTCCbVjK2JLI07WqatUr5p5qvBZMqLdYfMfHFVump939WIYjeRbfcCOAS4vlU7hU5A58jjhWESu5ptO/l/Hj7+d7bt5H8y3z6qMcFUmySbIbP1vTRWaROAOhd3gimkiQxJQ4CJwEkEfSetHaa8G3C/JAjmQ/utmT0s6SngbklnAa8DJ4br/5FgJNkK4EPgzFbu31WAJJshC5l7La2mI2+ycuWiyQQj6fPAyQSJZR3wO4ILM49o7U7N7GVg/yzl64DRWcoN+F5r9+sqS5ITVkb7XjIXfTaW1ozMPhO0Kxe5ruTfAjwOnGVmK8Kyl81szyLG12I1NTVWW1ubdhgVya+0d64yFe1KfoI7V/4LeFTSzZJGk300l6sy5XKlfRxaMpuDq27+mdkq11Qx95vZScDewHzgfGA3STdIGlOk+FwJyjfhZznJdzKopmTq4uGfma1y1WAAMLMPzGy2mR0D9AGWEIzwckVUSt+KMv0TldA8Fj0ZZDvGlZRMXXH4Z2arvPeDKVeV1gfj96FJRrQ/KZNs/Bi7alXMPhhXQvxbUeGaU9uL1saae4xLqVbpXCnyGoyrOHHU9goZKee1Sldp4q7BFHShpXPlJNt1Is0dWl3IxYx+PYpzuXmCcSUjrutrsl2A2dyr3wtJHkle6OlcJfAE40pGklOgNLe24cnDudbzBONKRpJNTp4wnCs+H0XmSkYxrq/xkV/OFY8nmDKR1Imx2k64fpW1c8VT9QmmXE6wSZ0Yq+2E69cTOVc8Vd8HE2fHcpKzDCfVP1FtQ229L8a54qn6BBPnCTbJUVBJnRibu12fqt85V6iqTzBxnriroTbgd1N0zhWq6vtg4lQpswzn6pfyPgznXKGKnmAk9ZX0qKTnJT0n6dywfJqkNyQtCZejI6+5RNIKScsljS12zNUmV8d/pSRR51zy0mgi2wR838yeltQZqJM0N/zdNWb2i+jKkvYFTgb2A/YAHpH0eTPbXNSoq0g1NPU555JX9ARjZquAVeHj9yQ9D/TO8ZJjgbvM7GPgFUkrgIOA0h5XXMZ8pJVzLg6p9sFI6g8cADwRFp0taamkGZK6hWW9gZWRl9WTOyE555wrAaklGEk7AfcC55nZu8ANwEBgCEEN55eZVbO8POtNbCRNllQrqXbt2rUJRO2cc65QqSQYSe0IkstsM7sPwMxWm9lmM9sC3EzQDAZBjaVv5OV9gDezbdfMpptZjZnV9OzZM7k34JxzLq80RpEJ+A3wvJn9KlLeK7La8cCz4eOHgJMl7ShpALAX8GSx4nXOOdcyaYwiOwQ4DXhG0pKw7AfARElDCJq/XgW+DWBmz0m6G1hGMALtez6CzDnnSl8ao8j+SvZ+lT/meM2VwJWJBeWccy52fiW/c865RHiCcc45lwhPMM455xLhCcY551wiPME455xLhCcY55xzifAE45xzLhGeYJxzziXCE4xzzrlEeIJxzjmXCE8wzjnnEuEJxjnnXCI8wTjnnEuEJxjnnHOJ8ATjnHMuEZ5gnHPOJcITjHPOuUR4gnHOOZcITzDOOecSUTYJRtI4ScslrZA0Je14nHPO5VYWCUZSW+B64MvAvsBESfumG5VzzrlcyiLBAAcBK8zsZTP7BLgLODblmJxzzuWwQ9oBFKg3sDLyvB4Y3nglSZOByeHTjyU9W4TYWmMX4K20gyiAxxkvjzNeHmd8vhDnxsolwShLmW1XYDYdmA4gqdbMapIOrDXKIUbwOOPmccbL44yPpNo4t1cuTWT1QN/I8z7AmynF4pxzrgDlkmCeAvaSNEDSZ4CTgYdSjsk551wOZdFEZmabJJ0NzAHaAjPM7Lk8L5uefGStVg4xgscZN48zXh5nfGKNUWbbdWU455xzrVYuTWTOOefKjCcY55xziSjbBCPpXEnPSnpO0nlhWXdJcyW9GP7sFpZL0rXhNDNLJQ1NMcarJb0QxnG/pJ3D8v6SNkpaEi43FiPGHHFOk/RGJJ6jI+tfEh7L5ZLGphzn7yIxvippSVhe1OMpaYakNdFrr1ryeZQ0KVz/RUmTUozx1DC2pZIWSto/8ppXJT0THtdYh7W2IM5RkjZE/s4/jLwm0emlmhnnRZEYn5W0WVL38HdpHM8Tw/+jLZJqGq2f9f+7RcfTzMpuAf4NeBboSDBQ4RFgL+DnwJRwnSnAz8LHRwN/IrieZgTwRIoxjgF2CNf5WSTG/sCzJXQspwEXZll/X+AfwI7AAOAloG1acTZa55fAD9M4nsBhwNDoPpv7eQS6Ay+HP7uFj7ulFOPBmX0TTNH0ROQ1rwK7lMixHAX8Ics22oafzT2Bz4Sf2X3TirPR68YDf0n5eO5DcFHlfKAmUp71/7ulx7NcazD7AIvN7EMz2wQsAI4nmD5mVrjOLOC48PGxwG0WWAzsLKlXGjGa2Z/D5wCLCa7pSVNTx7IpxwJ3mdnHZvYKsIJgKp9U45Qk4GvAnUWIZTtm9hjwdqPi5n4exwJzzextM3sHmAuMSyNGM1sYxgBF/pw281g2JfHppVoR50SK+DnNFqeZPW9my7Os3tT/d4uOZ7kmmGeBwyT1kNSR4BthX2A3M1sFEP7cNVw/21QzvVOKMeobBN9kMwZI+rukBZIOTTi+QuI8O2wimZGp6pPOscwXJ8ChwGozezFSlsbxjGru5zGNY9tUjFFnse3n1IA/S6pTMD1TMeSKc6Skf0j6k6T9wrK0Pqc5j2f42R0H3BspTuN4NiXWz2ZZXAfTmJk9L+lnBN/w3ieorm3K8ZKCppqJU74YJV0aPp8dFq0C+pnZOknDgAck7Wdm76YU5w3AjwiO048Imp++QQrHMk+cGY2/FaZyPAvU1DFM5djmIukIggTzxUjxIWb2pqRdgbmSXgi/JafhaeCzZva+gn7CBwiaeEvuWIbGA38zs2iNopSOZ1PHLVtlJO/xLNcaDGb2GzMbamaHEVT/XgRWZ5q+wp9rwtVTmWqmiRgJO2+PAU61sOEzrJKuCx/XEbR3fj7pGJuK08xWm9lmM9sC3MzWZrDUpu3JcTx3ACYAv4usm9rxjGju5zGNY9tUjEgaDNwCHJs5lgBm9mb4cw1wP8VpIs0ap5m9a2bvh4//CLSTtAvpfU6bPJ6hk2nUPJbS8WxKrJ/Nsk0wYbZHUj+Ck8udBNPHZEbeTAIeDB8/BJyuwAhgQ6YaW+wYJY0D/gv4ipl9GFm3p4L73iBpT4JvYS8nHWOOOKN9VMcTNFFBcCxPlrSjpAFhnE+mFWf4qy8BL5hZfWTd1I5nRHM/j3OAMZK6hU2SY8KyoscYHuP7gNPM7J+ZlSV1ktQ58ziMsRizljcV5+5h/xuSDiI4p60jvemlmvqbI6krcHijsrSOZ1Oa+v9u2fFszeiENBfgcWAZQVPJ6LCsBzCP4JvtPKB7WC6CG5a9BDxDZNRECjGuIGjLXBIuN4blJwDPhes+DYxP+VjeHh6rpeEHqVdk/UvDY7kc+HKacYblM4HvNFq3qMeTINmtAj4l+LZ3Vks+jwTNkCvC5cwUY7wFeCfyOa0Ny/cMj+k/wuN7acrH8uzI33kxcHBkO0cD/wyPc6pxhuufQdCBHt1GWsfz+PDxx8BqYE5k/az/3y05nj5VjHPOuUSUbROZc8650uYJxjnnXCI8wTjnnEuEJxjnnHOJ8ATjnHMuEZ5gXOIUzBy7JLLEPrNtln3uLOk/WvC6aZIubKI8M7v0MkkT44m0+CQdIMlUxJmwm4jjDEl7pBmDS5YnGFcMG81sSGT5aRH2uTPQ7ASTxzVmNoRgkr+bJLWLefstEl6w2Zz/5YnAX8OfaToD8ARTwTzBuFRI6qrg3hJfCJ/fKelb4eP3Jf1S0tOS5knqGZYPlPRwOCng45L2Dst3U3BvnX+Ey8HAT4GBYY3j6nC9iyQ9pWACz8sjsVwaxvIIwRTmOVkwoeaHBFPqNzcuJF2g4J4gz2rrfW1+Fq1xhTWm7zcVt4L73Twv6X8JLiS9TNI1kdd/S9Kvshx3AV8lOLmPkdQ+sr0XJN0SxjVb0pck/U3BvU0OCtfrLumBMJbFCqaT2a7mF26jfyTOmxXcf+TPkjpI+ipQA8wO/0Yd8h13V4bivmrUF18aL8Bmtl4RvgQ4KSw/ClhEMO3Ew5H1jWCeNoAfAteFj+cR3gMGGE54Tw2COcjOCx+3BbrS6H4wBFNwTCe4ir4N8AeC+2QMI7iaviPQheAK+mz3wZmWKSe4t8bjkd81J67M/joBOxFcvX1AuCyIbHMZ0C9H3P2BLcCIcP1OBFdYtwufLwQGZXkfXwTmhY9/C0wIH/cnmDx0ULifOmBGuN9jgQfC9f4HmBo+PhJY0vj4hM+fDbeZ2e6QsPxu4Ovh4/kUaVYNX9JZynI2ZVd2NlrQtLQNM5sr6USCaVP2j/xqC1snrrwDuE/STgQ3wfp9OPUUBDdFguBEd3q4zc3ABm29vUDGmHD5e/h8J4J5ljoD91s4L5ykXPMrnR/WsvYkvE9LC+L6Yri/D8LX3wccambXSto17JPoCbxjZq9LOqeJuF8HXrPgfjKY2QeS/gIcI+l5gkTzTJb3MJHgXh6EP08jmHMM4JXMayQ9R5CITNIzBIkCggR1QrjPvyi4fULXHMcss90l4eO6yLZchfME41IT9hvsA2wkuItjfROrZqYLX58tURW6O+AnZnZToxjOo/Bp3K8xs19ImgDcJmlgC+LKNh16xj0EzVe7szUJNBV3f+CDRq+/BfgB8AJw63Y7Dib/PAH4ioLbRQjooXCyRYJ5qTK2RJ5vYeu5oqnp3DexbZN7+8jj6HY3A94cViW8D8al6XzgeYJv1TMineZtCE60AKcAf7XgPi6vhDWeTMd2ptYzD/huWN5WUhfgPYLaScYc4BthjQNJvRXMzvwYcHzYL9CZ4H4dOZnZfUAtMKkFcT0GHCepo4LZc48nmMQTgqRycvje78kTd7a4niCYUv0Ust8x8UvAP8ysr5n1N7PPEtz4Kt/dIaMeA04NYxkFvBUeg1cJmg6RNJTgdrv5NP4buQrjCcYVQwdtO0z5p5I+D3wT+L6ZPU5w4vrvcP0PgP0k1RE0M10Rlp8KnCUpM/Ns5pat5wJHhE05dcB+Fty/5G9hZ/PVZvZngj6HReF69wCdzexpgua4JQQn28zJPp8rgAvCWlhz4nqaYPbnJ4EngFvM7O8AZvYcwQn3Ddt6V8SsceeI626CG1q9k+V3EwnuNxJ1L0FCKtQ0oEbSUoKBFJmp6e8FuktaQpBU/5n95duYCdzonfyVy2dTdiVH0vtmtlPacZQjSX8gaMqbl3YsznkNxrkKoODC0n8SDKjw5OJKgtdgnHPOJcJrMM455xLhCcY551wiPME455xLhCcY55xzifAE45xzLhH/H5LWfO4fDpaWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Scatter plot of Actual Recovery Amount vs. Expected Recovery Amount \n",
    "plt.scatter(x=df['expected_recovery_amount'], y=df['actual_recovery_amount'], c=\"g\", s=2)\n",
    "plt.xlim(900, 1100)\n",
    "plt.ylim(0, 2000)\n",
    "plt.xlabel(\"Expected Recovery Amount\")\n",
    "plt.ylabel(\"Actual Recovery Amount\")\n",
    "plt.legend(loc=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Statistical analysis:  recovery amount\n",
    "<p>Just as we did with age, we can perform statistical tests to see if the actual recovery amount has a discontinuity above the \\$1000 threshold. We are going to do this for two different windows of the expected recovery amount \\$900 to \\$1100 and for a narrow range of \\$950 to \\$1050 to see if our results are consistent.</p>\n",
    "<p>Again, the statistical test we will use is the Kruskal-Wallis test, a test that makes no assumptions about the distribution of the actual recovery amount.</p>\n",
    "<p>We will first compute the average actual recovery amount for those customers just below and just above the threshold using a range from \\$900 to \\$1100.  Then we will perform a Kruskal-Wallis test to see if the actual recovery amounts are different just above and just below the threshold.  Once we do that, we will repeat these steps for a smaller window of \\$950 to \\$1050.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KruskalResult(statistic=30.246000000000038, pvalue=3.80575314300276e-08)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute average actual recovery amount just below and above the threshold\n",
    "by_recovery_strategy['actual_recovery_amount'].describe().unstack()\n",
    "\n",
    "# Perform Kruskal-Wallis test\n",
    "Level_0_actual = era_900_1100.loc[df['recovery_strategy']=='Level 0 Recovery']['actual_recovery_amount']\n",
    "Level_1_actual = era_900_1100.loc[df['recovery_strategy']=='Level 1 Recovery']['actual_recovery_amount']\n",
    "stats.kruskal(Level_0_actual,Level_1_actual) \n",
    "\n",
    "# Repeat for a smaller range of $950 to $1050\n",
    "era_950_1050 = df.loc[(df['expected_recovery_amount']<1050) & \n",
    "                      (df['expected_recovery_amount']>=950)]\n",
    "Level_0_actual = era_950_1050.loc[df['recovery_strategy']=='Level 0 Recovery']['actual_recovery_amount']\n",
    "Level_1_actual = era_950_1050.loc[df['recovery_strategy']=='Level 1 Recovery']['actual_recovery_amount']\n",
    "\n",
    "stats.kruskal(Level_0_actual,Level_1_actual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Regression modeling: no threshold\n",
    "<p>We now want to take a regression-based approach to estimate the impact of the program at the \\$1000 threshold using the data that is just above and just below the threshold. In order to do that, we will build two models.  The first model does not have a threshold while the second model will include a threshold.</p>\n",
    "<p>The first model predicts the actual recovery amount (outcome or dependent variable) as a function of the expected recovery amount (input or independent variable). We expect that there will be a strong positive relationship between these two variables.  </p>\n",
    "<p>We will examine the adjusted R-squared to see the percent of variance that is explained by the model.  In this model, we are not trying to represent the threshold but simply trying to see how the variable used for assigning the customers (expected recovery amount) relates to the outcome variable (actual recovery amount).</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/halilkacmaz/anaconda3/envs/tfdeeplearning/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.\n",
      "  return ptp(axis=axis, out=out, **kwargs)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>    <td>actual_recovery_amount</td> <th>  R-squared:         </th> <td>   0.261</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                      <td>OLS</td>          <th>  Adj. R-squared:    </th> <td>   0.256</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>Least Squares</td>     <th>  F-statistic:       </th> <td>   63.78</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>                <td>Sun, 05 May 2019</td>    <th>  Prob (F-statistic):</th> <td>1.56e-13</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                    <td>23:58:10</td>        <th>  Log-Likelihood:    </th> <td> -1278.9</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>         <td>   183</td>         <th>  AIC:               </th> <td>   2562.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>             <td>   181</td>         <th>  BIC:               </th> <td>   2568.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>                 <td>     1</td>         <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>nonrobust</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "              <td></td>                <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>                    <td>-1978.7597</td> <td>  347.741</td> <td>   -5.690</td> <td> 0.000</td> <td>-2664.907</td> <td>-1292.612</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expected_recovery_amount</th> <td>    2.7577</td> <td>    0.345</td> <td>    7.986</td> <td> 0.000</td> <td>    2.076</td> <td>    3.439</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>64.493</td> <th>  Durbin-Watson:     </th> <td>   1.777</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td> 185.818</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 1.463</td> <th>  Prob(JB):          </th> <td>4.47e-41</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 6.977</td> <th>  Cond. No.          </th> <td>1.80e+04</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.8e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                              OLS Regression Results                              \n",
       "==================================================================================\n",
       "Dep. Variable:     actual_recovery_amount   R-squared:                       0.261\n",
       "Model:                                OLS   Adj. R-squared:                  0.256\n",
       "Method:                     Least Squares   F-statistic:                     63.78\n",
       "Date:                    Sun, 05 May 2019   Prob (F-statistic):           1.56e-13\n",
       "Time:                            23:58:10   Log-Likelihood:                -1278.9\n",
       "No. Observations:                     183   AIC:                             2562.\n",
       "Df Residuals:                         181   BIC:                             2568.\n",
       "Df Model:                               1                                         \n",
       "Covariance Type:                nonrobust                                         \n",
       "============================================================================================\n",
       "                               coef    std err          t      P>|t|      [0.025      0.975]\n",
       "--------------------------------------------------------------------------------------------\n",
       "const                    -1978.7597    347.741     -5.690      0.000   -2664.907   -1292.612\n",
       "expected_recovery_amount     2.7577      0.345      7.986      0.000       2.076       3.439\n",
       "==============================================================================\n",
       "Omnibus:                       64.493   Durbin-Watson:                   1.777\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              185.818\n",
       "Skew:                           1.463   Prob(JB):                     4.47e-41\n",
       "Kurtosis:                       6.977   Cond. No.                     1.80e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.8e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Import statsmodels\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# Define X and y\n",
    "X = era_900_1100['expected_recovery_amount']\n",
    "y = era_900_1100['actual_recovery_amount']\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "# Build linear regression model\n",
    "model = sm.OLS(y, X).fit()\n",
    "predictions = model.predict(X)\n",
    "\n",
    "# Print out the model summary statistics\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Regression modeling: adding true threshold\n",
    "<p>From the first model, we see that the regression coefficient is statistically significant for the expected recovery amount and the adjusted R-squared value was about 0.26. As we saw from the graph, on average the actual recovery amount increases as the expected recovery amount increases. We could add polynomial terms of expected recovery amount (such as the squared value of expected recovery amount) to the model but, for the purposes of this practice, let's stick with using just the linear term.</p>\n",
    "<p>The second model adds an indicator of the true threshold to the model. If there was no impact of the higher recovery strategy on the actual recovery amount, then we would expect that the relationship between the expected recovery amount and the actual recovery amount would be continuous.     </p>\n",
    "<p>In this case, we know the true threshold is at \\$1000.  </p>\n",
    "<p>We will create an indicator variable (either a 0 or a 1) that represents whether or not the expected recovery amount was greater than \\$1000. When we add the true threshold to the model, the regression coefficient for the true threshold represents the additional amount recovered due to the higher recovery strategy.  That is to say, the regression coefficient for the true threshold measures the size of the discontinuity for customers just above and just below the threshold.</p>\n",
    "<p>If the higher recovery strategy did help recovery more money, then the regression coefficient of the true threshold will be greater than zero.  If the higher recovery strategy did not help recover more money than the regression coefficient will not be statistically significant.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>    <td>actual_recovery_amount</td> <th>  R-squared:         </th> <td>   0.314</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                      <td>OLS</td>          <th>  Adj. R-squared:    </th> <td>   0.307</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>Least Squares</td>     <th>  F-statistic:       </th> <td>   41.22</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>                <td>Sun, 05 May 2019</td>    <th>  Prob (F-statistic):</th> <td>1.83e-15</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                    <td>23:59:17</td>        <th>  Log-Likelihood:    </th> <td> -1272.0</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>         <td>   183</td>         <th>  AIC:               </th> <td>   2550.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>             <td>   180</td>         <th>  BIC:               </th> <td>   2560.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>                 <td>     2</td>         <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>nonrobust</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "              <td></td>                <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>                    <td>    3.3440</td> <td>  626.274</td> <td>    0.005</td> <td> 0.996</td> <td>-1232.440</td> <td> 1239.128</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expected_recovery_amount</th> <td>    0.6430</td> <td>    0.655</td> <td>    0.981</td> <td> 0.328</td> <td>   -0.650</td> <td>    1.936</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>indicator_1000</th>           <td>  277.6344</td> <td>   74.043</td> <td>    3.750</td> <td> 0.000</td> <td>  131.530</td> <td>  423.739</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>65.977</td> <th>  Durbin-Watson:     </th> <td>   1.906</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td> 186.537</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 1.510</td> <th>  Prob(JB):          </th> <td>3.12e-41</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 6.917</td> <th>  Cond. No.          </th> <td>3.37e+04</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 3.37e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                              OLS Regression Results                              \n",
       "==================================================================================\n",
       "Dep. Variable:     actual_recovery_amount   R-squared:                       0.314\n",
       "Model:                                OLS   Adj. R-squared:                  0.307\n",
       "Method:                     Least Squares   F-statistic:                     41.22\n",
       "Date:                    Sun, 05 May 2019   Prob (F-statistic):           1.83e-15\n",
       "Time:                            23:59:17   Log-Likelihood:                -1272.0\n",
       "No. Observations:                     183   AIC:                             2550.\n",
       "Df Residuals:                         180   BIC:                             2560.\n",
       "Df Model:                               2                                         \n",
       "Covariance Type:                nonrobust                                         \n",
       "============================================================================================\n",
       "                               coef    std err          t      P>|t|      [0.025      0.975]\n",
       "--------------------------------------------------------------------------------------------\n",
       "const                        3.3440    626.274      0.005      0.996   -1232.440    1239.128\n",
       "expected_recovery_amount     0.6430      0.655      0.981      0.328      -0.650       1.936\n",
       "indicator_1000             277.6344     74.043      3.750      0.000     131.530     423.739\n",
       "==============================================================================\n",
       "Omnibus:                       65.977   Durbin-Watson:                   1.906\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              186.537\n",
       "Skew:                           1.510   Prob(JB):                     3.12e-41\n",
       "Kurtosis:                       6.917   Cond. No.                     3.37e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 3.37e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create indicator (0 or 1) for expected recovery amount >= $1000\n",
    "df['indicator_1000'] = np.where(df['expected_recovery_amount']<1000, 0, 1)\n",
    "era_900_1100 = df.loc[(df['expected_recovery_amount']<1100) & \n",
    "                      (df['expected_recovery_amount']>=900)]\n",
    "\n",
    "# Define X and y\n",
    "X = era_900_1100[['expected_recovery_amount','indicator_1000']]\n",
    "y = era_900_1100['actual_recovery_amount']\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "# Build linear regression model\n",
    "model = sm.OLS(y,X).fit()\n",
    "\n",
    "# Print the model summary\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Regression modeling: adjusting the window\n",
    "<p>The regression coefficient for the true threshold was statistically significant with an estimated impact of around \\$278 and a 95 percent confidence interval of \\$132 to \\$424.  This is much larger than the incremental cost of running the higher recovery strategy which was \\$50 per customer. At this point, we are feeling reasonably confident that the higher recovery strategy is worth the additional costs of the program for customers just above and just below the threshold.  </p>\n",
    "<p>Before showing this to our managers, we want to convince ourselves that this result wasn't due just to us choosing a window of \\$900 to \\$1100 for the expected recovery amount. If the higher recovery strategy really had an impact of an extra few hundred dollars, then we should see a similar regression coefficient if we choose a slightly bigger or a slightly smaller window for the expected recovery amount.  Let's repeat this analysis for the window of expected recovery amount from \\$950 to \\$1050 to see if we get similar results.</p>\n",
    "<p>The answer? Whether we use a wide window (\\$900 to \\$1100) or a narrower window (\\$950 to \\$1050), the incremental recovery amount at the higher recovery strategy is much greater than the \\$50 per customer it costs for the higher recovery strategy.  So we can say that the higher recovery strategy is worth the extra \\$50 per customer that the bank is spending.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>    <td>actual_recovery_amount</td> <th>  R-squared:         </th> <td>   0.283</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                      <td>OLS</td>          <th>  Adj. R-squared:    </th> <td>   0.269</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>Least Squares</td>     <th>  F-statistic:       </th> <td>   18.99</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>                <td>Mon, 06 May 2019</td>    <th>  Prob (F-statistic):</th> <td>1.12e-07</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                    <td>23:34:39</td>        <th>  Log-Likelihood:    </th> <td> -692.92</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>         <td>    99</td>         <th>  AIC:               </th> <td>   1392.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>             <td>    96</td>         <th>  BIC:               </th> <td>   1400.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>                 <td>     2</td>         <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>nonrobust</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "              <td></td>                <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>                    <td> -279.5243</td> <td> 1840.707</td> <td>   -0.152</td> <td> 0.880</td> <td>-3933.298</td> <td> 3374.250</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expected_recovery_amount</th> <td>    0.9189</td> <td>    1.886</td> <td>    0.487</td> <td> 0.627</td> <td>   -2.825</td> <td>    4.663</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>indicator_1000</th>           <td>  286.5337</td> <td>  111.352</td> <td>    2.573</td> <td> 0.012</td> <td>   65.502</td> <td>  507.566</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>39.302</td> <th>  Durbin-Watson:     </th> <td>   1.955</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td>  82.258</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 1.564</td> <th>  Prob(JB):          </th> <td>1.37e-18</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 6.186</td> <th>  Cond. No.          </th> <td>6.81e+04</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 6.81e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                              OLS Regression Results                              \n",
       "==================================================================================\n",
       "Dep. Variable:     actual_recovery_amount   R-squared:                       0.283\n",
       "Model:                                OLS   Adj. R-squared:                  0.269\n",
       "Method:                     Least Squares   F-statistic:                     18.99\n",
       "Date:                    Mon, 06 May 2019   Prob (F-statistic):           1.12e-07\n",
       "Time:                            23:34:39   Log-Likelihood:                -692.92\n",
       "No. Observations:                      99   AIC:                             1392.\n",
       "Df Residuals:                          96   BIC:                             1400.\n",
       "Df Model:                               2                                         \n",
       "Covariance Type:                nonrobust                                         \n",
       "============================================================================================\n",
       "                               coef    std err          t      P>|t|      [0.025      0.975]\n",
       "--------------------------------------------------------------------------------------------\n",
       "const                     -279.5243   1840.707     -0.152      0.880   -3933.298    3374.250\n",
       "expected_recovery_amount     0.9189      1.886      0.487      0.627      -2.825       4.663\n",
       "indicator_1000             286.5337    111.352      2.573      0.012      65.502     507.566\n",
       "==============================================================================\n",
       "Omnibus:                       39.302   Durbin-Watson:                   1.955\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               82.258\n",
       "Skew:                           1.564   Prob(JB):                     1.37e-18\n",
       "Kurtosis:                       6.186   Cond. No.                     6.81e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 6.81e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Redefine era_950_1050 so the indicator variable is included\n",
    "era_950_1050 = df.loc[(df['expected_recovery_amount']<1050) & \n",
    "                      (df['expected_recovery_amount']>=950)]\n",
    "\n",
    "# Define X and y \n",
    "X = era_950_1050[['expected_recovery_amount','indicator_1000']]\n",
    "y = era_950_1050['actual_recovery_amount']\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "# Build linear regression model\n",
    "model = sm.OLS(y,X).fit()\n",
    "\n",
    "# Print the model summary\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}