Solucion_Tarea6MC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analítica y Ciencia de Datos\n",
    "\n",
    "## CIDE - Otoño 2015\n",
    "\n",
    "\n",
    "### Solución Tarea 6: Modelos de Clasificación."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sys\n",
    "sys.path.append('/Users/danielvaughan/Documents/Python Scripts')\n",
    "import olsdan as ols\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Punto 1: Vamos a limpiar los datos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = pd.read_excel('/Users/danielvaughan/Documents/Clase CIDE/Slides/Analitica/datasets/data_tarea6.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 287,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id_cliente</th>\n",
       "      <th>MOB</th>\n",
       "      <th>EDAD_CLIENTE</th>\n",
       "      <th>SEXO</th>\n",
       "      <th>UN_MESES_SINPAGO</th>\n",
       "      <th>FACT_PROM_6M_ANTES</th>\n",
       "      <th>PAYRATIO_PROM_6M_ANTES</th>\n",
       "      <th>BALANCE_PROM_6M_ANTES</th>\n",
       "      <th>USOVECES_PROM_6M_ANTES</th>\n",
       "      <th>TKTPROM_PROM_6M_ANTES</th>\n",
       "      <th>SCORE_LAST</th>\n",
       "      <th>SCORE_FIRST</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>32</td>\n",
       "      <td>43</td>\n",
       "      <td>M</td>\n",
       "      <td>0</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>771</td>\n",
       "      <td>767</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>61</td>\n",
       "      <td>41</td>\n",
       "      <td>M</td>\n",
       "      <td>0</td>\n",
       "      <td>2353.033333</td>\n",
       "      <td>1.652819</td>\n",
       "      <td>2357.733333</td>\n",
       "      <td>3.633333</td>\n",
       "      <td>88.201330</td>\n",
       "      <td>616</td>\n",
       "      <td>611</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>35</td>\n",
       "      <td>37</td>\n",
       "      <td>M</td>\n",
       "      <td>0</td>\n",
       "      <td>21.542857</td>\n",
       "      <td>0.203771</td>\n",
       "      <td>6774.142857</td>\n",
       "      <td>0.028571</td>\n",
       "      <td>17.514286</td>\n",
       "      <td>561</td>\n",
       "      <td>587</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>27</td>\n",
       "      <td>35</td>\n",
       "      <td>F</td>\n",
       "      <td>1</td>\n",
       "      <td>1150.444444</td>\n",
       "      <td>3.484363</td>\n",
       "      <td>1506.703704</td>\n",
       "      <td>1.185185</td>\n",
       "      <td>229.744092</td>\n",
       "      <td>605</td>\n",
       "      <td>616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>42</td>\n",
       "      <td>28</td>\n",
       "      <td>F</td>\n",
       "      <td>0</td>\n",
       "      <td>215.523810</td>\n",
       "      <td>1.291662</td>\n",
       "      <td>775.309524</td>\n",
       "      <td>0.142857</td>\n",
       "      <td>80.047619</td>\n",
       "      <td>752</td>\n",
       "      <td>758</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   id_cliente  MOB  EDAD_CLIENTE SEXO  UN_MESES_SINPAGO  FACT_PROM_6M_ANTES  \\\n",
       "0           1   32            43    M                 0            0.000000   \n",
       "1           2   61            41    M                 0         2353.033333   \n",
       "2           3   35            37    M                 0           21.542857   \n",
       "3           4   27            35    F                 1         1150.444444   \n",
       "4           5   42            28    F                 0          215.523810   \n",
       "\n",
       "   PAYRATIO_PROM_6M_ANTES  BALANCE_PROM_6M_ANTES  USOVECES_PROM_6M_ANTES  \\\n",
       "0                0.000000               0.000000                0.000000   \n",
       "1                1.652819            2357.733333                3.633333   \n",
       "2                0.203771            6774.142857                0.028571   \n",
       "3                3.484363            1506.703704                1.185185   \n",
       "4                1.291662             775.309524                0.142857   \n",
       "\n",
       "   TKTPROM_PROM_6M_ANTES  SCORE_LAST  SCORE_FIRST  \n",
       "0               0.000000         771          767  \n",
       "1              88.201330         616          611  \n",
       "2              17.514286         561          587  \n",
       "3             229.744092         605          616  \n",
       "4              80.047619         752          758  "
      ]
     },
     "execution_count": 287,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.12 (500, 12)\n",
      "(440, 12)\n"
     ]
    }
   ],
   "source": [
    "# Eliminar aquellos con filas cero: son 12% de los clientes\n",
    "ind_cero = data[[u'FACT_PROM_6M_ANTES',u'PAYRATIO_PROM_6M_ANTES', u'BALANCE_PROM_6M_ANTES',\n",
    "                 u'USOVECES_PROM_6M_ANTES',  u'TKTPROM_PROM_6M_ANTES']].sum(axis=1)==0\n",
    "print ind_cero.mean(), data.shape\n",
    "\n",
    "data.drop(data[ind_cero].index, axis=0, inplace=True)\n",
    "print data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.598465473146\n"
     ]
    }
   ],
   "source": [
    "# Género a una dummy: 59% de la muestra son hombres\n",
    "hombre_dum = 1.0*(data.SEXO=='M')\n",
    "hombre_dum[data.SEXO.isnull()==True] =  np.nan\n",
    "print hombre_dum.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.206818181818\n"
     ]
    }
   ],
   "source": [
    "# Variable dependiente: 20% de la muestra se ha atrasado por lo menos una vez\n",
    "yy = 1.0*(data.UN_MESES_SINPAGO>0)\n",
    "yy[data.UN_MESES_SINPAGO.isnull()] = np.nan\n",
    "print yy.mean()\n",
    "data.UN_MESES_SINPAGO.groupby(data.UN_MESES_SINPAGO).count()\n",
    "nobs = yy.shape[0]\n",
    "yy = yy.reshape((nobs,1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Arreglemos los regresores que vamos a utilizar\n",
    "antiguedad  = np.asarray(data.MOB).reshape((nobs,1))\n",
    "edad_clie   = np.asarray(data.EDAD_CLIENTE).reshape((nobs,1))\n",
    "edad2_clie  = edad_clie**2\n",
    "factu_antes = np.asarray(data.FACT_PROM_6M_ANTES).reshape((nobs,1))\n",
    "payratio    = np.asarray(data.PAYRATIO_PROM_6M_ANTES).reshape((nobs,1))\n",
    "saldo_antes = np.asarray(data.BALANCE_PROM_6M_ANTES).reshape((nobs,1))\n",
    "uso_antes   = np.asarray(data.USOVECES_PROM_6M_ANTES).reshape((nobs,1))\n",
    "tktprom_ant = np.asarray(data.TKTPROM_PROM_6M_ANTES).reshape((nobs,1))\n",
    "# SCORE: miremos que no haya negativos\n",
    "data.SCORE_LAST.loc[data.SCORE_LAST<0] = np.nan\n",
    "data.SCORE_FIRST.loc[data.SCORE_FIRST<0] = np.nan\n",
    "score_last  = np.asarray(data.SCORE_LAST).reshape((nobs,1))\n",
    "score_delta = np.asarray(data.SCORE_LAST - data.SCORE_FIRST).reshape((nobs,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xmat   = np.concatenate((antiguedad, edad_clie, factu_antes,\n",
    "                      payratio, saldo_antes, uso_antes, tktprom_ant, score_last, score_delta),axis=1)\n",
    "labs_x = np.array(['antiguedad', 'edad_clie', 'factu_antes','payratio', 'saldo_antes', 'uso_antes', \n",
    "          'tktprom_ant', 'score_last', 'score_delta'])\n",
    "\n",
    "xmat_cons = np.concatenate((np.ones((nobs,1)),antiguedad, edad_clie, factu_antes,\n",
    "                      payratio, saldo_antes, uso_antes, tktprom_ant, score_last, score_delta),axis=1)\n",
    "labs_x_cons = np.array(['constante','antiguedad', 'edad_clie', 'factu_antes','payratio', 'saldo_antes', 'uso_antes', \n",
    "          'tktprom_ant', 'score_last', 'score_delta'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(384, 1) (384, 9)\n"
     ]
    }
   ],
   "source": [
    "# Vamos a eliminar los nans\n",
    "# necesitamos concatenar Y y X, eliminar, y volver a separar\n",
    "predata = np.asarray(pd.DataFrame(np.concatenate((yy,xmat),axis=1)).dropna())\n",
    "nobs = predata.shape[0]\n",
    "yy   = predata[:,0].reshape((nobs,1))\n",
    "xmat = predata[:,1:]\n",
    "print yy.shape, xmat.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Punto 2: Lasso y Modelo Lineal de Probabilidad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Para no penalizar la constante: medias cero\n",
    "yy_std = yy - np.mean(yy)\n",
    "prexmat_std = xmat.copy()\n",
    "nobs, nvar = prexmat_std.shape\n",
    "prexmat_std -= prexmat_std.mean(axis=0).reshape((1,nvar))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>384.000000</td>\n",
       "      <td>3.840000e+02</td>\n",
       "      <td>3.840000e+02</td>\n",
       "      <td>3.840000e+02</td>\n",
       "      <td>3.840000e+02</td>\n",
       "      <td>3.840000e+02</td>\n",
       "      <td>3.840000e+02</td>\n",
       "      <td>3.840000e+02</td>\n",
       "      <td>384.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>-6.772360e-15</td>\n",
       "      <td>-4.556355e-13</td>\n",
       "      <td>-1.821460e-15</td>\n",
       "      <td>-1.005418e-12</td>\n",
       "      <td>-1.393561e-16</td>\n",
       "      <td>6.054416e-14</td>\n",
       "      <td>-3.789561e-14</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>38.269320</td>\n",
       "      <td>1.091924e+01</td>\n",
       "      <td>1.606062e+03</td>\n",
       "      <td>2.552374e+00</td>\n",
       "      <td>4.047146e+03</td>\n",
       "      <td>1.239815e+00</td>\n",
       "      <td>2.200555e+02</td>\n",
       "      <td>5.362701e+01</td>\n",
       "      <td>27.235766</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>-38.953125</td>\n",
       "      <td>-2.099740e+01</td>\n",
       "      <td>-8.172067e+02</td>\n",
       "      <td>-1.476013e+00</td>\n",
       "      <td>-1.204876e+04</td>\n",
       "      <td>-8.344817e-01</td>\n",
       "      <td>-1.538020e+02</td>\n",
       "      <td>-1.528724e+02</td>\n",
       "      <td>-110.140625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>-25.953125</td>\n",
       "      <td>-8.997396e+00</td>\n",
       "      <td>-7.020734e+02</td>\n",
       "      <td>-1.267160e+00</td>\n",
       "      <td>-1.909616e+03</td>\n",
       "      <td>-6.806357e-01</td>\n",
       "      <td>-1.212896e+02</td>\n",
       "      <td>-3.337240e+01</td>\n",
       "      <td>-7.140625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>-11.953125</td>\n",
       "      <td>-9.973958e-01</td>\n",
       "      <td>-4.706118e+02</td>\n",
       "      <td>-9.377056e-01</td>\n",
       "      <td>-1.167282e+03</td>\n",
       "      <td>-4.076962e-01</td>\n",
       "      <td>-7.164976e+01</td>\n",
       "      <td>8.127604e+00</td>\n",
       "      <td>1.859375</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>10.296875</td>\n",
       "      <td>9.002604e+00</td>\n",
       "      <td>5.899390e+01</td>\n",
       "      <td>1.805201e-01</td>\n",
       "      <td>3.836276e+02</td>\n",
       "      <td>1.863518e-01</td>\n",
       "      <td>2.289922e+01</td>\n",
       "      <td>4.112760e+01</td>\n",
       "      <td>7.859375</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>227.046875</td>\n",
       "      <td>2.700260e+01</td>\n",
       "      <td>1.366311e+04</td>\n",
       "      <td>2.329729e+01</td>\n",
       "      <td>3.633554e+04</td>\n",
       "      <td>9.765518e+00</td>\n",
       "      <td>1.518531e+03</td>\n",
       "      <td>9.712760e+01</td>\n",
       "      <td>149.859375</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                0             1             2             3             4  \\\n",
       "count  384.000000  3.840000e+02  3.840000e+02  3.840000e+02  3.840000e+02   \n",
       "mean     0.000000 -6.772360e-15 -4.556355e-13 -1.821460e-15 -1.005418e-12   \n",
       "std     38.269320  1.091924e+01  1.606062e+03  2.552374e+00  4.047146e+03   \n",
       "min    -38.953125 -2.099740e+01 -8.172067e+02 -1.476013e+00 -1.204876e+04   \n",
       "25%    -25.953125 -8.997396e+00 -7.020734e+02 -1.267160e+00 -1.909616e+03   \n",
       "50%    -11.953125 -9.973958e-01 -4.706118e+02 -9.377056e-01 -1.167282e+03   \n",
       "75%     10.296875  9.002604e+00  5.899390e+01  1.805201e-01  3.836276e+02   \n",
       "max    227.046875  2.700260e+01  1.366311e+04  2.329729e+01  3.633554e+04   \n",
       "\n",
       "                  5             6             7           8  \n",
       "count  3.840000e+02  3.840000e+02  3.840000e+02  384.000000  \n",
       "mean  -1.393561e-16  6.054416e-14 -3.789561e-14    0.000000  \n",
       "std    1.239815e+00  2.200555e+02  5.362701e+01   27.235766  \n",
       "min   -8.344817e-01 -1.538020e+02 -1.528724e+02 -110.140625  \n",
       "25%   -6.806357e-01 -1.212896e+02 -3.337240e+01   -7.140625  \n",
       "50%   -4.076962e-01 -7.164976e+01  8.127604e+00    1.859375  \n",
       "75%    1.863518e-01  2.289922e+01  4.112760e+01    7.859375  \n",
       "max    9.765518e+00  1.518531e+03  9.712760e+01  149.859375  "
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(prexmat_std).describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10b6708d0>"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEACAYAAACznAEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8FdX5+PHPk7BGspKQsCRCAIEALfhFFAWJVjAoCmhF\nUBHUilagqF0AkYJVERFqRdSqRQEpi/y0FkE2lwi1CqUFFyCEoIQlgJCEECCS7fn9kZvrTXKzB5MJ\nz/v1mheznDP33LnknJlzZuYRVcUYY4zxqe0CGGOMqRusQTDGGANYg2CMMcbFGgRjjDGANQjGGGNc\nrEEwxhgD1ECDICJxIpIgIntFZFIpaea5tn8pIj3LyysiM0TkkIhsd01x1S2nMcaYslWrQRARX2A+\nEAfEACNFpEuxNDcAHVS1IzAWeKUCeRX4s6r2dE3rqlNOY4wx5avuFUJvIElV96tqDrAcGFIszc3A\nIgBV3QIEiUhEBfJKNctmjDGmEqrbILQGDnosH3Ktq0iaVuXkneDqYlogIkHVLKcxxphyVLdBqOh7\nLyp7tv8K0A7oARwB5lYyvzHGmEpqUM38h4FIj+VICs70y0rTxpWmYWl5VfX7wpUi8jfg/eIfLCL2\nEiZjjKkCVfV6kl7dK4RtQEcRaSsijYDbgVXF0qwC7gYQkSuAk6p6rKy8ItLSI/8w4GtvH66qNrmm\n6dOn13oZ6tJkx8OOhR0P71NZqnWFoKq5IjIeWA/4AgtUdbeIPODa/qqqfiAiN4hIEnAGuKesvK5d\nPysiPSjokvoOeKA65TTGGFO+6nYZoaprgbXF1r1abHl8RfO61t9d3XIZY4ypHHtSuZ6IjY2t7SLU\nKXY8fmTHoig7HqWT8vqU6ioRUaeW3RhjaouIoOdpUNkYY0w9YQ2CMcYYwBoEY4wxLtYgGGOMAaxB\nMMYY42INgjHGGMAaBGOMMS7WIBhjjAGsQTDGGONiDYIxxhjA4Q3CkOVD2JS8qdxXuhpjjCmfoxuE\nQR0G8atVv6L333qz/cj22i6OMcY4muNfbpev+bz15VtM/mgyH939ETFhMbVdNGOMqbPKerldteMh\n1DYf8WF0j9H4iA/XL7meTWM20S64XW0XyxhjHKfaXUYiEiciCSKyV0QmlZJmnmv7lyLSs7y8IhIi\nIhtFJFFENohIUHnlGPXzUUzpO4Xr3rqOlMyU6n4tY4y54FSrQRARX2A+EAfEACNFpEuxNDcAHVS1\nIzAWeKUCeScDG1X1EuAj13K5HrrsIe6/9H5iF8aSlJZUna9mjDEXnOpeIfQGklR1v6rmAMuBIcXS\n3AwsAlDVLUCQiESUk9edx/Xv0IoWaHLfyfy2z2/p+0ZfNidvrur3MsaYC051xxBaAwc9lg8Bl1cg\nTWugVRl5w1X1mGv+GBDu7cP/vuwvhDaJIF8bcjr3DPn5+aBKoDZiXJN7mTD7VwxoN4Do4PaICELB\ngAoArn891xfdphRscS0WpvfchzdlbHKEipZfITcXfvihYP6n++DzrPC3rUv3WtSRQ1NnlPX3d6Er\ndmxUISfXh7xcIS/fh/z8so9ddRuEiv7ZVOQXFG/7U1UVEa+fc/fTj5LfSKERSFuQqGIJLoJvvk9E\nv69gKb1Rr7PGGFMNUspUbJsW/Kulpq/APg7+AAezK1Sq6jYIh4FIj+VICs70y0rTxpWmoZf1h13z\nx0QkQlWPikhLwGuVfvjDFEJDQ3n1i+vp3uZurm47uhpfxZlU4cwZSE8vmNLSfpz3nLytP3kSGjeG\n4ODSp5AQ7+sDAwvyGuN0qpCdDVlZcPZs0X8rMu+5riLr8/KgadOCyc+v4vMVTeu53KQJ+BQbGCir\nh6O6DcI2oKOItAVSgNuBkcXSrALGA8tF5ArgpKoeE5HUMvKuAkYDz7r+fc/bh0dERACQnZ9DSMOm\n1fwqtad4pV5WxV5828mT0KhR2RV4ly7e1wcFFeQ1pq7Jz69cJestTUUr6qws8PUtvfL1VtkW/tu8\nObRpU3YFXnxbw4Z1t9erWg2CquaKyHhgPeALLFDV3SLygGv7q6r6gYjcICJJwBngnrLyunY9C3hb\nRO4D9gPDyygDOfnZhDRqUp2vUm2qBf/RKnp27rnt5MmC/yRBQaVX7J07l77NKnXzU8jJqfzZc2Uq\naM/t2dkFZ7cVOQsuvj44uPJn0r6+tX106wbHP6mclpPD4n9fyd3dnyYkZGC19qkKp08XPfsurSL3\nNjVoULlul8JtQUHW/WIqT7VgUL+ilWxluzaK/wsVPwsuq7ujol0ddfUs2unq9ZPKKefO4Sf5iBT9\nKtnZcPAgHDgAJ054r9w91508WbRP3fNs3XP+kktKr9yb1O5FiqkD8vKq1p1RlQr6hx8Krg4rWsl6\nrm/RovJ91w0b1vbRNeeb4xuEg1nZNMjLZ/HihmzfDvv2wXffwfffQ+vWEBUFoaFFK/bIyJJ96dan\nXj95DhjW1Jl0Wdtzc0uvYMurfJs3r9wgY5Mm1tVhapajG4Tp02Fe4jmevTOfH35oQL9+MGYMtGtX\n0Bg0cPS3q7/y8wvObqvTnVGZitrHp+yz5dK6NIKDC/4fVeZMulEj6+owzuXoKjMtDUY9nE14bh4j\nft2QgIDaLpGz5ObCuXMFlXPh5LlckXnPPuyKVtDnzhV0zVXmTLpwOTCw8rfj2YmBMRXj6D+VF1+E\n8YnnaHQiDx8f53RwqsdTvtWtkKuSv3BSLeh2aNKkoIJu3LigAi1rvjBtYZdFQACEh1eui8TbvdHG\nmNrn6AYBICU7m4aUHFQuS2FlXHzyVmmWt87b9opUzj4+RStYz4rZ27y3bYGB5efzrLw9K/cmTezM\n2RhTlKOrhF//Gr44E0pgo1+zISSCs2cL7hTKyIBTpyAzs+DfwjsyCifwXlmWVnF6O1sOCSl9W0Uq\nequMjTF1jaOrpe7dYdnOc2iW0LSpEBlZcKdQQEDB2XNAAPj7/1ixF1bGdvucMcaU5OgH0/Ly82my\naRPxDe/k0ks/p0mTyPIzGmPMBaysB9McPbR3IieHAF9fVHMRsdN+Y4ypDkc3CCnnztGqcWNUcyo1\nqGyMMaYkZzcI2dm0atQI1RxH3XZqjDF1kaMbhCPZ2a4rhFy7QjDGmGpydC2acu6c+wqhcAyhzPCW\nxrGcevODMU7i7AYhO5vufn4lrhCs8qhfrJE35qfh6C6jgisEX8AHEUd/FWOMqXVVrkVFJERENopI\noohsEJGgUtLFiUiCiOwVkUnl5ReRtiKSJSLbXdPLpZUhJTublo187JZTY4ypAdU5rZ4MbFTVS4CP\nXMtFiIgvMB+IA2KAkSLSpQL5k1S1p2t6qLQCpJw7R8uGPjagbIwxNaA6DcLNwCLX/CJgqJc0vSmo\n3Perag6wHBhSifxlOp6TQ2gDsVtOjTGmBlSnQQhX1WOu+WNAuJc0rYGDHsuHXOvKy9/O1V0ULyJ9\nSytASIMGNPASPvNC4e/vz/79+3/yz42NjWXBggVVyjtmzBimTZtWwyUyxtSEMmtSEdkIRHjZNNVz\nQVVVRLzd2lN8nXhZVzx/ChCpqukicinwnoh0VdXMEntftIg/rW3GkSM/kJMTT2xsbFlfx9FiY2MZ\nNWoU9913n3tdZmbJQ/JTEJEq3/lTnbzGmMqLj48nPj6+QmnLbBBUdUBp20TkmIhEqOpREWkJfO8l\n2WHA841zbVzrALzmV9VsINs1/z8R2Qd0BP5XfOf/N348UzsGsGPHO/TpE1vWV3G8+lSJ2m3Bxvx0\nYmNji5wsP/HEE6WmrU6X0SpgtGt+NPCelzTbgI6uO4caAbe78pWaX0RCXYPRiEg0BY3Bt94K4MSn\nlGfNmkWHDh0ICAiga9euvPdewWFbuHAhffv25fe//z0hISFER0ezbt06AKZOncrmzZsZP348/v7+\n/OY3vwHAx8eHb78tODSpqancdNNNBAYG0rt3bx5//HH69esHwP79+/Hx8SE/P99djuLdPm+88QYx\nMTGEhIQQFxfHgQMH3Ns2btxI586dCQoKYsKECaiqu1Lft28f1157LaGhoYSFhXHXXXeRkZHhzrt9\n+3YuvfRSAgICGDFiBD8UBqQwxtQ51WkQZgEDRCQRuNa1jIi0EpE1AKqaC4wH1gO7gBWqurus/MDV\nwJcish1YCTygqie9FaD4U8pO0KFDB/71r39x6tQppk+fzl133cXRo0cB2Lp1K507dyY1NZU//OEP\n7u6hp59+mn79+vHSSy+RmZnJvHnzSux33Lhx+Pv7c+zYMRYtWsTixYvLvKrw7Lr55z//yTPPPMM/\n/vEPTpw4Qb9+/Rg5ciQAJ06c4NZbb2XmzJmkpqbSvn17PvvssyL7njp1KkeOHGH37t0cPHiQGTNm\nAJCdnc3QoUMZPXo06enp3Hbbbbzzzjv16mrHmHql8GzPaROgfz18WDMzv9ItW7pqoYKvVLaCaMLV\nn2pCjx499J///Ke++eab2qFDB/f6M2fOqIjosWPHVFU1NjZW//a3vxXJKyK6b98+zc3N1YYNG2pi\nYqJ72+OPP659+/ZVVdXvvvtORUTz8vLc22NjY3XBggWqqhoXF+eeV1XNy8tTPz8/TU5O1kWLFmmf\nPn2KfG6bNm2KpPf0j3/8Q3v27Kmqqp9++qm2atWqyPYrr7xSp02bVrGD41KR39QYUzGuvyev9aqj\nH++t6ptOa6pJqIrFixfTs2dPgoODCQ4O5ptvvuHEiROICBERP47f+/n5AXD69Gn3utLOrI8fP05u\nbi6RkT8O17Rp06bCZUpOTmbixInuMjVv3hyAw4cPc+TIkRL78vycY8eOMWLECNq0aUNgYCCjRo0i\nNTUVgJSUFFq3bl0k78UXX2xjCMbUUc5uENxjCM7oMkpOTmbs2LG89NJLpKWlkZ6eTrdu3SpUQZbV\nzRIWFkaDBg04ePDHO3w95y+66CIAzp49615X2E0FEBUVxWuvvUZ6erp7OnPmDH369KFly5ZF9qWq\nRZYfe+wxfH19+eabb8jIyOCtt95yj1W0bNmSw4cP4yk5Odm6jIypo5zdILjHEJwxqHzmzBlEhNDQ\nUPLz83nzzTf55ptvgPLvvAkPD2ffvn1et/n6+nLLLbcwY8YMsrKySEhI4K233nJXvGFhYbRu3Zq3\n3nqLvLw83njjjSL7evDBB5k5cya7du0CICMjg5UrVwJwww03sHPnTv7xj3+Qm5vLvHnzijQmp0+f\n5qKLLiIgIIDDhw/z3HPPubf16dOHBg0aMG/ePHJycnj33Xf5z3/+U4UjZ4z5KTi6QWjRqBH5+c4Z\nVI6JieG3v/0tffr0ISIigm+++Ya+ffu6B3iLnzl7Lk+cOJH/9//+HyEhITz88MMl9j1//nwyMjKI\niIhg9OjRjBw5kkaNGrm3v/766zz33HOEhoaya9currrqKve2oUOHMmnSJEaMGEFgYCDdu3dn/fr1\nAISGhrJy5UomT55MaGgoSUlJ9O3747OC06dP53//+x+BgYHcdNNN3Hrrre5yN2rUiHfffZeFCxfS\nvHlz3n77bW699daaOZjGmBonTu3PFRFVVdLSPuTAgWfo0eOjwvXWRw1MmjSJ77//njfffLO2i1Jt\n9psaU3Ncf09e+20dfYUAOO620/Nlz549fPXVV6gqW7du5Y033mDYsGG1XSxjjIM4o/O9DE57MO18\nyczMZOTIkaSkpBAeHs7vfvc7br755touljHGQRxfk1blttP6qFevXuzdu7e2i2GMcTDrMjLGGAPU\niwbBuoyMMaYmOL5BcNJtp8YYU5c5vkGwKwRjjKkZ9aBBsCsEY4ypCfWgQci1u4yMMaYG1IMGwTnv\nMqqM+Pj4Im8VrQxvAXGq89ndunVj06ZNVdqXMcY5HF+TWpfR+Vf4Aj5jTP1W5SsEEQkRkY0ikigi\nG0QkqJR0cSKSICJ7RWSSx/rbRGSniOSJyKXF8kxxpU8QkYFllcMGlY0xpmZUp8toMrBRVS8BPnIt\nF+GKjTwfiANigJEi0sW1+WtgGLCpWJ4YCmIvx7jyvSwipZbTabedpqSkcOutt9KiRQuio6N58cUX\nAcjKymLMmDGEhITQtWvXEq+JLi0WM0B+fj6/+93vCAsLo3379qxZs6ZCZUlLS+Oee+6hdevWhISE\nlPruo7Zt2/LRRwUvD1RVd1lCQ0O5/fbbSU9Pr8qhMMbUMdU5tb4Z6O+aXwTEU7JR6A0kqep+ABFZ\nDgwBdqtqgmtd8f0OAZapag6wX0SSXPv5wlshCgaVm1Sq4PJEzQRo0emVewNnfn4+N910E8OGDWPF\nihUcPHiQ6667jk6dOvHxxx/z3Xff8e2333L69Gni4uKKHJvCWMwRERG8/fbb3HXXXezbt4/w8HBe\ne+011qxZw44dO/Dz8+OWW26pUBCaUaNGERAQwK5du7jooov4/PPPvabzfDX3vHnzWLVqFZs2bSIs\nLIwJEyYwbtw4li5dWqljYYypg0qLrVneBKR7zIvnssf6XwKveyzfBbxYLM0nwKUeyy8Cd3os/w24\n1cu+VVU1KekPun//M8XjhdZJX3zxhUZFRRVZN3PmTB0zZoxGR0fr+vXr3etfe+01bdOmTan76tGj\nh65atUpVVa+55hp99dVX3ds2bNhQIoZycSkpKerj46MnT54sse2TTz4p8tlt27bVjz76SFVVu3Tp\n4p4v3E/Dhg3L/Kzqqsu/qTFOQxkxlcu8QhCRjUCEl01TizUqKiLeTpdr6iX2XvczY8YM0tI24evr\nz5AhVxAbG1tDH3d+JCcnk5KSQnBwsHtdXl4e/fr1IyUlpcidPVFRUUXyLl68mOeff579+/cDBZHK\nTpw4AcCRI0fKzOvNwYMHCQkJITAwsFLfYf/+/QwbNgwfnx978Ro0aMCxY8do2bJlpfZljDn/4uPj\niY+Pr1DaMhsEVR1Q2jYROSYiEap6VERaAt97SXYY8Lx3MhI4VE6Ziudp41pXwowZM9i7N42mTdvT\npk1sObutfVFRUbRr147ExMQS26Kjozlw4ABduhQMsRw4cMC9rTAW88cff0yfPn0QEXr27OkOGtOy\nZcsi6T3nSxMZGUlaWhoZGRmVahSioqJ488036dOnT4XzGGNqT2xsbJGT5SeeeKLUtNUZVF4FjHbN\njwbe85JmG9BRRNqKSCMKBotXeUnn2eG9ChghIo1EpB3QEdhaWiGcdNtp79698ff3Z/bs2WRlZZGX\nl8c333zDtm3bGD58OM888wwnT57k0KFD7sFmKDsWM8Dw4cOZN28ehw8fJj09nVmzZpVblpYtWzJo\n0CAeeughTp48SU5OToWeNXjwwQd57LHH3I3O8ePHWbXK209qjHGa6jQIs4ABIpIIXOtaRkRaicga\nAFXNBcYD64FdwApV3e1KN0xEDgJXAGtEZK0rzy7gbVf6tcBDWngq7IWTbjv18fFh9erV7Nixg+jo\naMLCwhg7diynTp1i+vTpXHzxxbRr1464uDjuvvtu90BuabGYC91///1cf/31/PznP6dXr15F4hqX\n5a233qJhw4Z07tyZ8PBw5s2b595WWv6JEydy8803M3DgQAICAujTpw9bt5baXhtjHMTxMZV37x5D\nUFB/Wra8p3C9xd+tZ+w3Nabm1POYyrmO6TIyxpi6rB40CPXzXUY1pVmzZvj7+5eYPvvss9oumjGm\njnF8TWoxlct2+vTp2i6CMcYh6sEVgnMGlY0xpi6rBw2Cc247NcaYuqweNAg2qGyMMTXB8Q1CwdtO\nrcvIGGOqy/ENgnUZGWNMzagHDYKzBpX37NlDjx49CAgIYP78+bVdHGOMcXNOTVoKp912Onv2bH7x\ni1+wY8eOKu9jzJgxREZG8uSTT9ZgySpn//79REdHk5ubW+TNp8YY53L8X7LTrhCSk5OJiYmp7WLU\nGHulhDH1Rz1oEJwzhnDttdcSHx/P+PHj8ff3Z968efTs2ZPAwECioqJKvJb2X//6F1deeSXBwcFE\nRUWxaNEiXn/9dZYuXcrs2bPx9/dnyJAhQMGL87799lt33jFjxjBt2rQyy3Py5EkGDx5MixYtCAkJ\n4aabbuLw4R/fNB4bG8sf//hH+vbtS0BAANdffz2pqakAXH311QAEBQXh7+/Pli1bAHjjjTeIiYkh\nJCSEuLi4Iq/ifuSRRwgPDycwMJCf/exn7Ny5sxpH0xhT40qLnFPXJ1xRtD7/vL2eOZNYPBpQeSGD\namaqgtjYWF2wYIGqqsbHx+s333yjqqpfffWVhoeH63vvvaeqqvv371d/f39dvny55ubmampqqu7Y\nsUNVVceMGaPTpk0rsl8R0X379rmXvaUpLjU1Vd99913NysrSzMxMve2223To0KHu7f3799cOHTro\n3r17NSsrS2NjY3Xy5Mnu8hWPyvbee+9phw4dNCEhQfPy8vSpp57SK6+8UlVV161bp//3f/+nGRkZ\nqqqakJCgR44cqdAxq9BvaoypEMqImFYPrhCq0GVUU01ClctckLd///507doVgO7duzNixAg+/fRT\nAJYuXcqAAQO4/fbb8fX1JSQkhJ///Ocl9lGRzylNSEgIw4YNo0mTJjRr1ozHHnvM/flQ8FbEe+65\nhw4dOtCkSROGDx/uHvvwtu+//vWvTJkyhU6dOuHj48OUKVPYsWMHBw4coFGjRmRmZrJ7927y8/Pp\n1KkTERHegvEZY2pLPWgQnNNlVKgw1sCWLVu45ppraNGiBUFBQbz66qvuLpmDBw8SHR19Xstx9uxZ\nHnjgAdq2bUtgYCD9+/cnIyOjSGXvWWk3bdq0zHcjJScnM3HiRIKDgwkODqZ58+YApKSkcM011zB+\n/HjGjRtHeHg4DzzwAJmZmefvyxljKq3KDYKIhIjIRhFJFJENIhJUSro4EUkQkb0iMslj/W0islNE\n8kTkUo/1bUUkS0S2u6aXyyqH0waVPd1xxx0MHTqUQ4cOcfLkSR588EF3ZRwVFcW+ffu85vMWvMbP\nz4+zZ8+6l48cOVJukJy5c+eSmJjI1q1bycjI4NNPP/XskiuTt31HRUXx2muvkZ6e7p7OnDnDFVdc\nAcCECRPYtm0bu3btIjExkeeee67czzHG/HSqc4UwGdioqpcAH7mWixARX2A+EAfEACNFpItr89fA\nMMBb3MYkVe3pmh4qqxBOu+3U0+nTpwkODqZRo0Zs3bqVpUuXurfdcccdfPjhh6xcuZLc3FxSU1P5\n8ssvAQgPDy8ygAzQo0cP/v73v5OXl8e6desqFA7z9OnTNG3alMDAQNLS0rzGWi2tcQgLC8PHx6dI\no/Xggw8yc+ZMdu3aBUBGRgYrV64EYNu2bWzZsoWcnBz8/Pxo0qQJvr6+5ZbRGPPTqU6DcDOwyDW/\nCBjqJU1vCir3/aqaAywHhgCoaoKqlow2X0lOfpfRyy+/zB//+EcCAgJ48sknuf32293boqKi+OCD\nD5g7dy7NmzenZ8+efPXVVwDcd9997Nq1i+DgYG655RYAXnjhBd5//32Cg4NZunQpw4YNK/fzH374\nYbKysggNDeXKK69k0KBBJc78PZdFxL3s5+fH1KlTueqqqwgODmbr1q0MHTqUSZMmMWLECAIDA+ne\nvTvr168H4NSpU4wdO5aQkBDatm1LaGgov//976t3AI0xNarKITRFJF1Vg13zAqQVLnuk+SVwvare\n71q+C7hcVSd4pPkE+K2q/s+13Bb4BtgLZACPq+q/vHy+qiqbNjXlqqtS8fX1K1xv98bXM/abGlNz\nygqhWWbnu4hsBLzdCjLVc0G14LZHL+mq8lecAkSqarprbOE9Eemqql5HIAtebufMKwRjjKlLymwQ\nVHVAadtE5JiIRKjqURFpCXzvJdlhINJjORI4VM5nZgPZrvn/icg+oCPwv+Jpp0+fTnJyHp988hTX\nXHMNsbGxZe36gjRz5kyeeeaZEuuvvvpq1qxZUwslMsb8lOLj44mPj69Q2up0Gc0GUlX1WRGZDASp\n6uRiaRoAe4BfUHDmvxUYqaq7PdJ8AvxOVf/rWg4F0lU1T0SiKRh07qaqJ4vtW/Pystm0qSmxsbme\n6617oZ6x39SYmlNWl1F1BpVnAQNEJBG41rWMiLQSkTUAqpoLjAfWA7uAFYWNgYgME5GDwBXAGhFZ\n69pvf+BLEdkOrAQeKN4YFFLNdewdRsYYU9dU+QqhtomI5uRk8PnnrenXL9NzvZ1N1jP2mxpTc87X\nFUKtc+JTysYYU1c5vEFw7lPKxhhT1zi6QbBbTmvG3//+d66//vraLoYxppY5egzh7Nnv2LGjP336\nJHuut/7mMjgx0pn9psbUHBtDuIDk5+dXKJ1VsMaY4hzeIDjrttO2bdsya9YsunbtSkhICPfeey/n\nzp0jPT291MhlK1eupFevXkX28+c//5mhQwteHTVmzBh+/etfc8MNN9CsWTPi4+NZs2ZNqZHYPCOd\nBQQE8MUXX7Bw4UL69evnTvPvf/+byy67jKCgIHr37s3nn39+vg+NMaYuKC1yTl2fAM3M/FK3bu3m\nLRpQnXTxxRdr9+7d9dChQ5qWlqZXXXWVPv7442VGLvvhhx80JCREd+/e7d5Pjx499N1331VV1dGj\nR2tgYKD++9//dqcvLxJb8Uhnb775pvbt21dVC6KoBQUF6ZIlSzQvL0+XLVumwcHBmpqaev4PUCnq\n8m9qjNNQRsQ0R9+iU9UuI6ngY9zlfn4lX5UhIowfP57WrVsDMHXqVCZMmMCTTz5Z5O2kjz32GNde\ney0AjRs3Zvjw4SxZsoSnnnqKnTt3kpyczODBg937HDp0KH369HGn79+/v3tfnpHYhgwZUm5X0Zo1\na+jUqRN33nknACNGjGDevHm8//77jB49ulLf1xjjLA5vEKp222llK/KaFBn546udoqKiSElJISsr\ni4cffpj169eTnp4OFMQqUFVEhNGjR3PHHXfw1FNP8dZbb3H77bfTsGFDr/uEgkhskydPZufOnWRn\nZ3Pu3DmGDx9eofKlpKQQFRVVZN3FF1/s7sIyxtRfjh5DcOJtpwcOHCgy36pVq3Ijl11xxRU0atSI\nTZs2sWzZMkaNGlXmZ3iLxFY42FxeFLXWrVuTnJxcZF1ycjJt2rSpytc1xjiIoxsEpwXHUVVefvll\nDh8+TFq95wDaAAAgAElEQVRaGk8//TQjRowgMzOz3Mhlo0aNYvz48TRq1Igrr7yyyD6L8xaJrbAh\n8BbpzNOgQYNITExk2bJl5ObmsmLFChISEtxdVMaY+svhDUKOo55UFhHuuOMOBg4cSPv27enYsSOP\nP/54hSKXjRo1ip07d3LXXXeV2GfxtGVFYvOMdBYSEsKWLVuK7KN58+asXr2auXPnEhoaypw5c1i9\nejUhISHn6agYY+oKRz+YduLEGg4ffpGf/Wyt5/o6e499u3btWLBggXvAuDKysrIIDw9n+/bttG/f\n/jyUru6qy7+pMU5Tjx9Mu3DeZfTKK6/Qu3fvC64xMMb8dBxdm14oTyq3bdsWEeG9996r7aIYY+ox\nhzcIzrpC+O6776qUb//+/TVbEGOM8aLKXUYiEiIiG0UkUUQ2iEhQKeniRCRBRPaKyCSP9c+JyG4R\n+VJE3hWRQI9tU1zpE0RkYGllcOJtp8YYU1dVZwxhMrBRVS8BPnItFyEivsB8IA6IAUaKSBfX5g1A\nV1X9OZAITHHliQFud6WPA14WEa/ldNptp8YYU5dVp0G4GVjkml8EDPWSpjeQpKr7VTUHWA4MAVDV\njapa+GrOLUDhk09DgGWqmqOq+4Ek135KcNptp8YYU5dVp0EIV9VjrvljQLiXNK2Bgx7Lh1zrirsX\n+MA138qVrrw8qOY46m2nxhhTl5V5ei0iG4EIL5umei6oqoqItxvFy715XESmAtmqurSMZF7347RB\nZWOMqcvKrE1VdUBp20TkmIhEqOpREWkJfO8l2WHA881rkXic/YvIGOAG4Bdl5GnjWlfC3LmrycvL\nJCRkBrGxscTW4kvrjDGmLoqPjye+gm94rvKTyiIyG0hV1WdFZDIQpKqTi6VpAOyhoMJPAbYCI1V1\nt4jEAXOB/qp6wiNPDLCUgnGD1sCHQActVlAR0eTkWeTkpNG+/bOe6+vdU60+Pj4kJSURHR1dYtvC\nhQtZsGABmzdvroWS/TTq429qTG05X08qzwIGiEgicK1rGRFpJSJrAFQ1FxgPrAd2AStUdbcr/4tA\nM2CjiGwXkZddeXYBb7vSrwUeKt4YFCq47dS6jGrTmDFjmDZtWm0XwxhTA6pcm6pqGnCdl/UpwI0e\ny2spqNiLp+tYxr5nAjPLL4M9h2CMMTXF3mX0E3v22Wdp06YNAQEBdO7cmY8//pitW7fSp08fgoOD\nadWqFRMmTCAnJ8dr/tTUVG6++WYCAwO5/PLLS7zGuirxkN98801iYmIICAigffv2vPbaa+5t8fHx\ntGnThj//+c+Eh4fTqlUrFi5cCMBrr73G0qVLmT17Nv7+/gwZMgQoCLJz66230qJFC6Kjo3nxxRfd\n+9u6dSu9evUiMDCQiIgIfvvb31b2EBpjzhNn1abFFNx2GlDpfPESXyOfH6uxlUq/Z88eXnrpJbZt\n20ZERAQHDhwgNzeXkydP8sILL9CrVy8OHjzIoEGDePnll5k4cWKJfYwbNw4/Pz+OHj3Kt99+y/XX\nX+8eW0hLS+PGG29k/vz5jBw5krfffpsbb7yRpKSkMl9fHR4ezpo1a2jXrh2bNm1i0KBBXHbZZfTs\n2ROAY8eOcerUKVJSUtiwYQO//OUvGTZsGGPHjuXzzz8nMjKSP/3pTwDk5+dz0003MWzYMFasWMHB\ngwe57rrr6NSpEwMHDmTixIk88sgj3HnnnZw9e5avv/66UsfQGHMelRZsua5PgO7d+4geODDXWwDp\nOmnv3r3aokUL/fDDDzU7O7vUdM8//7wOGzbMvSwium/fPs3NzdWGDRvqnj173Nsee+wx7du3r6qq\nLl68WC+//PIi++rTp48uXLiwUuUcOnSovvDCC6qq+sknn2jTpk01Ly/Pvb1Fixa6ZcsWVVUdM2aM\nPv744+5tX3zxhUZFRRXZ38yZM/Wee+5RVdWrr75ap0+frsePH69weeryb2qM07j+nrzWqw7vMnLW\noHKHDh34y1/+wowZMwgPD2fkyJEcOXKExMREBg8eTMuWLQkMDGTq1KmkpqaWyH/8+HFyc3NLxGUu\nVNV4yGvXruWKK66gefPmBAcH88EHHxT5/ObNm+Pj8+N/FT8/P06fPu11X8nJyaSkpBAcHOyennnm\nGb7/vuCu5AULFpCYmEiXLl3o3bs3a9asKbNsxpifjqMbBCe+3G7kyJFs3ryZ5ORkRIRJkybx0EMP\nERMTQ1JSEhkZGTz99NPuGMiewsLCaNCgQYm4zIWqEg/53Llz3HrrrfzhD3/g+++/Jz09nRtuuKHC\nt3kWj9YWFRVFu3btSE9Pd0+nTp1i9erVQEGjuHTpUo4fP86kSZP45S9/SVZWVoU+yxhzfjm6QXDa\noHJiYiIff/wx586do3HjxjRt2hQfHx8yMzPx9/fHz8+PhIQEXnnlFa/5fX19ueWWW5gxYwZZWVns\n2rWLRYsWuSvlqsRDzs7OJjs7m9DQUHx8fFi7di0bNmyo8HcKDw/n22+/dS/37t0bf39/Zs+eTVZW\nFnl5eXzzzTds27YNgCVLlnD8+HEAAgMDEZEiVx/GmNrj6L9Ep912eu7cOaZMmUJYWBgtW7bk+PHj\nzJo1izlz5rB06VICAgIYO3YsI0aMKHLm7Tk/f/58Tp8+TUREBPfeey/33nuve1tV4iH7+/szb948\nhg8fTkhICMuWLXPfLeTt84u777772LVrF8HBwdxyyy34+PiwevVqduzYQXR0NGFhYYwdO5ZTp04B\nsH79erp164a/vz+PPPIIy5cvp3HjxpU+lsaYmufomMo7d95BSMggIiLu8lxvT7XWM/abGlNz6nFM\nZXvbqTHG1BSHNwgWIKeimjVrhr+/f4nps88+q+2iGWPqCOeMyHrhtNtOa1Npt4kaY0whR18hOPG2\nU2OMqasc3SA47bZTY4ypyxzeINgVgjHG1BSHNwi5dpeRMcbUEIc3CDaobIwxNaXKDYKIhIjIRhFJ\nFJENIhJUSro4EUkQkb0iMslj/XMisltEvhSRd0Uk0LW+rYhkuaKouSOpeWNdRsYYU3Oqc4UwGdio\nqpcAH7mWixARX2A+EAfEACNFpItr8wagq6r+HEgEpnhkTVLVnq7podIKYIPK58eMGTMYNWpUbRfD\nGPMTq06DcDOwyDW/CBjqJU1vCir3/aqaAywHhgCo6kZVLXyl5xag9FdylsJuOzXGmJpTnQYhXFWP\nueaPAeFe0rQGDnosH3KtK+5e4AOP5Xau7qJ4EelbWgGcdoXg4+NT5M2gngHqT5w4weDBgwkODqZ5\n8+ZcffXV7vf37N69m9jYWIKDg+nWrRvvv/9+uZ+1Zs0aevbsSWBgIFFRUTzxxBPubfv378fHx4fF\nixdz8cUXExYWxsyZBSGs161bxzPPPMOKFSvw9/d3R03LyMjgvvvuo1WrVrRp04Zp06a5X9GdlJRE\n//79CQoKIiwsjBEjRtTMATPG/KTKrE1FZCMQ4WXTVM8F1YKoXl7SlftGMhGZCmSr6lLXqhQgUlXT\nReRS4D0R6aqqmcXzvv76cSIiXqJhwyBiY2OJjY0t7+MAiI8v/e2dlREbW70XromI+02ic+fOJTIy\nkhMnTgDwxRdfICLk5ORw00038atf/YoPP/yQzZs3M2TIELZt28Yll1xS6r6bNWvGkiVL6Nq1K19/\n/TUDBgygR48eRd5k+tlnn5GYmMiePXvo3bs3t956K3FxcTz22GPs27ePxYsXu9OOGTOGiIgI9u3b\nx+nTpxk8eDCRkZGMHTuWadOmERcXx6effkp2drb7VdfGmNoXHx9PfHx8xRKXFkqtvAlIACJc8y2B\nBC9prgDWeSxPASZ5LI8BPgOalPE5nwCXelmvn33WSn/44ZC38HB1UmEozEJjxozRadOmqarqH//4\nRx0yZIgmJSUVybNp0yaNiIgosm7kyJE6Y8aMSn32xIkT9ZFHHlFV1e+++05FRA8fPuze3rt3b12x\nYoWqqk6fPl3vuusu97ajR49q48aNNSsry71u6dKles0116iq6t13361jx47VQ4eK/hY1pS7/psY4\nDecphOYqYLRrfjTwnpc024COrjuHGgG3u/IhInHA74EhqvpDYQYRCXUNRiMi0UBH4NsSe6Z+3Haq\nrm6h3//+93To0IGBAwfSvn17nn32WaAgLKZnyEyoWFjMLVu2cM0119CiRQuCgoJ49dVXS4TljIj4\n8eKvvLCYOTk5tGzZ0h0W88EHH3QHupk9ezaqSu/evenWrRtvvvlm5Q6CMaZOqE5tOgt4W0TuA/YD\nwwFEpBXwuqreqKq5IjIeWA/4AgtUdbcr/4tAI2Cjq9vkcy24o6g/8ISI5AD5wAOqetJbAZx226mf\nnx9nz551Lx85csRd2Tdr1ow5c+YwZ84cdu7cybXXXstll11G69atOXjwIKrq7l5KTk6mc+fOZX7W\nHXfcwW9+8xvWr19Po0aNeOSRR9zdUeUpHhAnMjKSxo0bk5qa6jW6WXh4OK+99hpQ0A113XXX0b9/\nf6Kjoyv0ecaYuqHKVwiqmqaq16nqJao6sLDSVtUUVb3RI91aVe2kqh1U9RmP9R1V9WItdnupqr6j\nqt1c6/5PVUuNwu60QeUePXrw97//nby8PNatW8emTZvc21avXk1SUhKqSkBAAL6+vvj6+nL55Zfj\n5+fH7NmzycnJIT4+ntWrV5c7cHv69GmCg4Np1KgRW7duZenSpWVGPvMUERHB/v373VcvLVu2ZODA\ngTz66KNkZmaSn5/Pvn373OVfuXIlhw4dAiAoKMjCYhrjUI7+q3XabacvvPAC77//PsHBwSxdupRh\nw4a5tyUlJTFgwAD8/f258sorGTduHP3796dhw4a8//77rF27lrCwMMaPH89bb71V5oAywMsvv8wf\n//hHAgICePLJJ7n99tuLbC+rcbjtttuAgpCcvXr1AmDx4sVkZ2cTExNDSEgIt912G0ePHgVg27Zt\nXHHFFfj7+zNkyBDmzZtH27Ztq3KIjDG1yNEhND/5xJerr/4BH58Gnust3GI9Y7+pMTWn3obQhDxc\n48/GGGOqydENgkiDCveL10ddu3b1GhZz2bJltV00Y4wDOWdE1gsnDSifDzt37qztIhhj6hGHXyE4\nZ0DZGGPqOmsQjDHGAI5vEC7sLiNjjKlJDm8Q7ArBGGNqisMbBGdfIRS+hrrwNdLGGFObHN0g+Pg4\n7wqhbdu2fPzxx+Wmu1Aai4oeD2PM+efoBsGJVwiVfeq2qk/o6o+vCa/T7ClkY+oOhzcIzrpCGDVq\nFAcOHOCmm27C39+flStXFtn+zjvv0K5dO3bu3MnVV18NFLwsLiAggC+++IKFCxdy1VVXMWHCBIKC\ngujSpUuRs+vY2Fgef/xxrrrqKi666CK+++47/v3vf3PZZZcRFBRE7969+fzzz4uknzZtGldddRX+\n/v7cfPPNnDhxgjvvvJPAwEB69+5NcnJyud9r4sSJREVFERgYSK9evfjXv/7l3jZjxgyGDx/O6NGj\nCQgIoFu3bvz3v//1ejzmzJlTreNrjKmm0gIl1PUJ0P/859LSgj/UWW3bttWPPvpIVX8MVJObm6tv\nvPGGdujQwR1AZ//+/SoimpeX58775ptvaoMGDfQvf/mL5ubm6ooVKzQwMFDT09NVVbV///568cUX\n665duzQvL0+PHj2qQUFBumTJEs3Ly9Nly5ZpcHCwpqWludN37NhRv/32W83IyNCYmBjt0KGDfvTR\nR5qbm6t333233nPPPeV+pyVLlmhaWprm5eXp3LlzNSIiQs+dO6eqBcF2mjRpomvXrtX8/HydMmWK\nXnHFFV6PR2nq+m9qjJNwngLk1LqqdhkVhq6s7lRTnn/+eebMmcOnn37qjiGgpXSjtGjRgokTJ+Lr\n68vw4cPp1KkTq1evdn+vMWPG0KVLF3x8fNiwYQOdOnXizjvvxMfHhxEjRtC5c2dWrVrlTn/PPffQ\nrl07AgICGDRoEJdccgnXXnstvr6+3HbbbWzfvr3c8t95550EBwfj4+PDo48+yrlz59izZ497e79+\n/YiLi0NEuOuuu/jyyy+re8iMMeeBwxuEqnUZldY6VnaqKXPnzmXcuHG0atWq3LStW7cusnzxxRdz\n5MgR97JndLWUlBSioqJKpE9JSXEvh4eHu+ebNGlCixYtiiyXFkXN05w5c4iJiSEoKIjg4GAyMjKK\nBOPx/Aw/Pz9++OGHej9YbowTVblBEJEQEdkoIokiskFEgkpJFyciCSKyV0Qmeax/UkS+FJEdIvKR\niER6bJviSp8gIgNLL4MzB5WL27BhA0899RTvvvtumemAEqEzk5OTizQknvlat25dYgwgOTm5RKNS\n3meWZfPmzTz33HOsXLmSkydPkp6eTmBgYIUbzAv55YTG1DXVuUKYDGxU1UuAj1zLRbhiI88H4oAY\nYKSIdHFtnq2qP1fVHhTEY57uyhNDQezlGFe+l0XEazmdeNtpeHg4+/btK7Kua9eurFu3jnHjxvH+\n++8DEBYWho+PT4m033//PfPmzSMnJ4eVK1eSkJDADTfc4N7uWRHfcMMNJCYmsmzZMnJzc1mxYgUJ\nCQkMHjzYa/qqXPVkZmbSoEEDQkNDyc7O5k9/+hOnTp2qcH5vx8MYUzuq0yDcDCxyzS8ChnpJ0xtI\nUtX9qpoDLAeGAKhqpke6ZkBhH8MQYJmq5qjqfiDJtZ8SnHaXEcCUKVN46qmnCAkJ4Z133nGfIf/s\nZz9j9erV3H///axfvx4/Pz+mTp3KVVddRUhICFu2bEFEuPzyy9m7dy9hYWFMmzaNd955h+DgYPf+\nPc+4Q0JCWL16NXPnziU0NJQ5c+awevVqQkJCvKb3NjZS3hl8XFwccXFxXHLJJbRt25amTZsW6aYq\nb5+FxyM4OJg///nPFTmExpjzpMoR00QkXVWDXfMCpBUue6T5JXC9qt7vWr4LuFxVJ7iWnwZGAVlA\nb1XNEJEXgS9U9e+uNH8D1qrqO8X2rV99dRPdu68qXq56e1/7woULWbBgAZs3b67tovyk6vNvasxP\nrayIaWV2wovIRiDCy6apnguqqiLi7S+2zL9iVZ0KTBWRycBfgHtKS+pt5Suv7KVFixlAwT31sbGx\nZX2cMcZccOLj44mPj69Q2jIbBFUdUNo2ETkmIhGqelREWgLfe0l2GIj0WI4EDnlJtxT4oJQ8bVzr\nShg37md07Tqj9C9Qz9T07a4VtXnz5iLjFJ7lqcx4gTHmp1f8ZPmJJ54oNW11uoxmA6mq+qzrDD9I\nVScXS9MA2AP8AkgBtgIjVXW3iHRU1b2udBMo6DIa5RpUXkrBuEFr4EOggxYrqIjozp13EhOzpHi5\nrHuhnrHf1JiaU+Uuo3LMAt4WkfuA/cBw14e1Al5X1RtVNVdExgPrAV9ggaruduV/RkQ6AXnAPuDX\nAKq6S0TeBnYBucBDxRuDH7+Y8247NcaYuqrKVwi1TUQ0IeFXdOr0evH1djZZz9hvakzNKesK4YJ8\nUtkYY0xJDm8QrMvIGGNqisMbBLtCMMaYmuLwBsGuEIwxpqY4vEGwK4SaMGbMGKZNm1bbxTDG1DJH\nNwhOfLldXVQTD7wtXLiQfv361VCJjDG1wdENgnUZ/ai6MRrstk5jjMMbBOddITz77LO0adOGgIAA\nOnfuzMcff0x+fj4zZ86kQ4cOBAQE0KtXLw4dKnjDR3kxkYvHUE5ISGDAgAE0b96czp07l4jbXJ70\n9HQGDx5MixYtCAkJ4aabbioSg2HhwoW0b9+egIAAoqOjWbp0KQkJCTz44IN8/vnn+Pv7F3mbqjHG\nQWoqethPPQF64MBcLY4KxN+FmpkqKyEhQSMjI/XIkSOqqpqcnKz79u3T2bNna/fu3TUxMVFVVb/6\n6itNTU3V1NTUcmMie8ZQPnnypLZp00YXLlyoeXl5un37dg0NDdVdu3aVWa4xY8bo448/rqqqqamp\n+u6772pWVpZmZmbqbbfdpkOHDlVV1dOnT2tAQIC7nEePHtWdO3eqqurChQu1b9++lT8oFVCR39QY\nUzHU35jKVQ2hWTNTZfn6+nLu3Dl27txJTk4OUVFRREdHs2DBAp5++mk6duwIQPfu3QkJCWHNmjXl\nxkT2jKG8bt062rVrx+jRo/Hx8aFHjx7ccsstlbpKCAkJYdiwYTRp0oRmzZrx2GOP8emnn7q3+/j4\n8PXXX5OVlUV4eDgxMTGuY2pdTsY43QXZINSWDh068Je//IUZM2YQHh7OyJEjSUlJ4eDBg7Rv375E\n+orERPaMoZycnMyWLVsIDg52T0uXLuXYsWMVLuPZs2d54IEHaNu2LYGBgfTv35+MjAxUlYsuuogV\nK1bw17/+lVatWjF48GD27NlThSNhjKmLHN4gOG9QeeTIkWzevJnk5GREhEmTJhEZGUlSUlKJtBWJ\niex5d1BUVBT9+/cnPT3dPWVmZvLSSy+VW67C/cydO5fExES2bt1KRkYGn376aZEB64EDB7JhwwaO\nHj1K586duf/++0uUwxjjTI5uEJx222liYiIff/wx586do3HjxjRp0oQGDRrwq1/9imnTppGUlISq\n8tVXX5GWllbpmMiDBw8mMTGRJUuWkJOTQ05ODv/5z39ISEgos1yeFf7p06dp2rQpgYGBpKWlFXl3\n+vfff88///lPzpw5Q8OGDbnooovw9fUFCmIjHzp0iJycnJo8ZMaYn5CjGwSnXSGcO3eOKVOmEBYW\nRsuWLTlx4gTPPPMMjz76KMOHD2fgwIEEBgZy//3388MPP1Q6JnKzZs3YsGEDy5cvp3Xr1rRs2ZIp\nU6aQnZ1dZrk8n0N4+OGHycrKIjQ0lCuvvJJBgwa5t+Xn5/P888/TunVrmjdvzubNm3nllVcA+MUv\nfkHXrl2JiIigRYsWNX3ojDE/AUe//vrYsRW0aDG8+Hob4Kxn7Dc1pubU49dfO+sKwRhj6rIqNwgi\nEiIiG0UkUUQ2iEhQKeniRCRBRPaKyCSP9U+KyJciskNEPhKRSNf6tiKSJSLbXdPLpZfBWWMItalr\n1674+/uXmJYtW1bbRTPG1BHVjal8QlVnuyr6YC0ZU9mXgpjK1wGHgf/wY0xlf1XNdKWbAPxcVX8l\nIm2B91W1ezmfrydOrKV587ji6617oZ6x39SYmnO+uoxuBha55hcBQ72k6Q0kqep+Vc0BlgNDAAob\nA5dmwInKFsC6jIwxpuZUp0EIV9XCJ56OAeFe0rQGDnosH3KtA0BEnhaRA8BoYJZHunau7qJ4Eelb\nWgGcdtupMcbUZWWeYovIRiDCy6apnguqqiLi7Zq+zOt8VZ0KTBWRycDzwD1AChCpqukicinwnoh0\nLXZFAcCsWYto0uQToOBFb7GxsWV9nDHGXHDi4+OJj4+vUNrqjCEkALGqelREWgKfqGrnYmmuAGao\napxreQqQr6rPFksXBXygqt28fM4nwG9V9X/F1mtGxhYCAnoXT2/9zfWM/abG1JzzNYawioKuHlz/\nvuclzTago+vOoUbA7a58iEhHj3RDgO2u9aGuwWhEJBroCHzrrQB2l5ExxtSc6jQIs4ABIpIIXOta\nRkRaicgaAFXNBcYD64FdwApV3e3K/4yIfC0iO4BY4Leu9VcDX4rIdmAl8ICqnvRWABtUPj8qE/3M\nwm8aU39UuUZV1TQKbictvj4FuNFjeS2w1ku6X5ay33eBdytSBrtCqH2er72Ij49n1KhRHDx4sJxc\nxpi6yJ5Uric8X1BXG59tjHE+RzcITrzttC6G0ExNTeXmm28mMDCQyy+/nH379hXZXt4+RYSzZ88y\naNAgUlJS8Pf3JyAggKNHj7J161b69OlDcHAwrVq1YsKECfZGVGPqKEefYlf1CkGeqJl39+v0yp0Z\n79mzh5deeolt27YRERHBgQMHyM3NZe7cuSxfvpy1a9fSsWNHvv76a/z8/EhLS+PGG29k/vz5jBw5\nkrfffpsbb7yRffv2ERwcDMCSJUtYu3YtnTp1IjMzk27duvHUU0+xfv16vvrqKwYMGEC3bt3o0qVL\nqeUaN24cfn5+HD16lG+//Zbrr7+e6OhoAM6cOcOAAQPK3Keq4ufnx7p167jrrruKdBmlpKTwwgsv\n0KtXLw4ePMigQYN4+eWXmThxYmUPtzHmfCsttmZdnwA9d+5YafFC66S9e/dqixYt9MMPP9Ts7Gz3\n+k6dOumqVatKpF+8eLFefvnlRdb16dNHFy5cqKqqsbGxOn36dPe25cuXa79+/YqkHzt2rD7xxBOl\nlik3N1cbNmyoe/bsca977LHH3PGRy9unZzzmTz75RNu0aVPqZ6mqPv/88zps2LAy0xRXl39TY5yG\nMmIqO/wKwVldRp4hNHfu3Mn111/P3Llzz0sIzUK5ubncfffdpZbp+PHj5ObmFtmP52dWZZ+eEhMT\nefTRR/nvf//L2bNnyc3NpVevXhXKa4z5aTl6DMGJg8p1LYRmWFgYDRo04MCBA+51nvMV2WdhGbyF\n0fz1r39NTEwMSUlJZGRk8PTTT5Ofn1/WITLG1BKHNwjOukKoiyE0fX19ueWWW5gxYwZZWVns2rWL\nRYsWuSv3G2+8scx96o9deISHh5OamsqpU6fc+z99+jT+/v74+fmRkJDgjrBmjKl7HN4gOOsKoa6G\n0Jw/fz6nT58mIiKCe++9l3vvvde9zd/fv8x9ej6H0LlzZ0aOHEl0dDQhISEcPXqUOXPmsHTpUgIC\nAhg7diwjRozweiVhjKl9jg6hmZ+fX6Jysffe1D/2mxpTc+pxCE070zTGmJri6AbBVJyF0DTGlMfR\nXUbeym7dC/WP/abG1Jx622VkjDGm5liDYIwxBnD4u4xKY4PNxhhTeVW+QhCREBHZKCKJIrJBRIJK\nSRcnIgkisldEJnnZ/lsRyReREI91U1zpE0RkYGXKVdo7Omxy9mSMOf+q02U0GdioqpcAH7mWi3CF\nwpwPxAExwEgR6eKxPRIYACR7rIuhINRmjCvfyyJiXVvlqGgQ7QuFHY8f2bEoyo5H6apT0d4MLHLN\nL+7VLSEAAAQySURBVAKGeknTG0hS1f2qmgMspyB+cqE/A38olmcIsExVc1R1P5Dk2o8pg/0nL8qO\nx4/sWBRlx6N01WkQwlX1mGv+GBDuJU1rwDOe4iHXOkRkCHBIVb8qlqeVK12JPMYYY86fMgeVRWQj\nEOFl01TPBVVVEfHW0eu181dEmgKPUdBd5F5dRlGsE9kYY86zKj+YJiIJQKyqHhWRlsAnqtq5WJor\ngBmqGudangLkA2soGHc460raBjgMXA7cA6Cqs1x51gHTVXVLsX1bI2GMMVWgpTyYVp0GYTaQqqrP\nishkIEhVJxdL0wDYA/wCSAG2AiNVdXexdN8B/6eqaa5B5aUUjBu0Bj4EOqjdamKMMedVdcYQZgED\nRCQRuNa1jIi0EpE1AKqaC4wH1gO7gBXFGwMXd2WvqruAt13p1wIPWWNgjDHnn2PfZWSMMaZmOfL+\n/vIedqvPRCRSRD4RkZ0i8o2I/Ma1vkIPCtZXIuIrIttF5H3X8gV7PEQkSET+n4jsFpFdInL5hXo8\nXA+57hSRr0VkqYg0vlCPRUU4rkEo72G3C0AO8IiqdgWuAMa5vn+5DwrWcxMp6GYsvOS9kI/HC8AH\nqtoF+BmQwAV4PESkLXA/cKmqdgd8gRFcgMeiohzXIFD+w271mqoeVdUdrvnTwG4KBt8r8qBgvSQi\nbYAbgL/x4+3LF+TxEJFAoJ+qvgEF43iqmsGFeTxOUXAC5ee6wcWPgptbLsRjUSFObBD+f3v3zxpF\nFEZh/DmIATH2WqwkiLZisFARBEkjSEorJY29gijoF7ASU9lIsEghiEqIpWBhIYgSQVBsBPEPxIgS\nEFuPxb3iGghMobts7vlVOzO7cDkw++7O3PfOhs1uram/gA4AT+nWKLhZXQcuUqY0/9ZqHpPAF0m3\nJC1LuilpOw3mYfsbcA14TykEa7Yf0mAWXY1iQchdcEDSOHAPOGf7e/+xOiuriZwknQRWbb9gg+bG\nlvKgNJtOATdsTwE/WHdJpJU8JO0BzgMTlBUQxiWd7n9PK1l0NYoF4RPQ69vu8fdSF5uepK2UYrBg\ne7Hu/ixpZz2+C1gd1vgG7AgwU3tZbgPHJS3Qbh4fKUvCPKvbdykFYqXBPA4CT2x/rVPg7wOHaTOL\nTkaxIDwH9kqakDRGWRl1achjGhiVhz3MA69tz/UdWgJm6+tZYHH9Zzcj21ds92xPUm4YPrJ9hnbz\nWAE+SNpXd00Dr4AHtJfHG+CQpG31vJmmTDxoMYtORrIPQdIJYI4ya2De9tUhD2lgJB0FHgMv+fNX\n9zKlC/wOsBt4B5yyvTaMMQ6LpGPABdsz9fkaTeYhaT/lBvsY8JayHMwWGsxD0iXKl/5PYBk4C+yg\nwSy6GMmCEBER/94oXjKKiIj/IAUhIiKAFISIiKhSECIiAkhBiIiIKgUhIiKAFISIiKhSECIiAoBf\n2J0M2Xk2XQYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b31c6d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn import linear_model\n",
    "\n",
    "L = 100\n",
    "grid_lamb = np.linspace(0.001,7,L)\n",
    "BetaLasso = np.zeros((L,nvar))\n",
    "for l in range(L):\n",
    "    lambvec = grid_lamb[l]\n",
    "    clf = linear_model.Lasso(alpha=lambvec)\n",
    "    clf.fit(prexmat_std,yy_std)\n",
    "    #print pd.DataFrame(clf.coef_, index = labs_x, columns=['Scikit-Learn'])\n",
    "    BetaLasso[l,:] = clf.coef_\n",
    "pd.DataFrame(BetaLasso, columns= labs_x).plot()\n",
    "plt.legend(loc='best', fontsize=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+4VWWd9/H3h4OkAir+SBMxzDQtJUElFH8cf9WxdLCn\nKfPSscaci5nGtOlpUp+rjKbJcqZmbDJ5tLTLaSxtjHwMJQp1I6CCyBFQwBERBPyNSgqRB/k+f9x7\ny+Zwzv5x2Oesvff5vK7rXOy11r32+m7B777P977XvRQRmJlZ8xmQdQBmZtY7nODNzJqUE7yZWZNy\ngjcza1JO8GZmTcoJ3sysSVWU4CW1SVom6SlJl3dxfIKkhZLaJT0q6dT8/hGS7pf0hKTHJV1a6w9g\nZmZdU7l58JJagCeB04G1wCPAeRGxtKjN4IjYkH99JPCbiHi/pP2A/SLiMUlDgEeBc4rPNTOz3lFJ\nD34ssDwiVkZEB3AbMKG4QSG55w0BXsnvfyEiHsu/fhNYCuxfi8DNzKy0ShL8cGB10faa/L5tSDpH\n0lJgGrBdKUbSSGA0MLcngZqZWXUqSfAVrWUQEXdGxOHA2cDPi4/lyzN3AJfle/JmZtbLBlbQZi0w\nomh7BKkX36WImCVpoKS9ImKdpJ2AXwP/FRF3dm4vyYvhmJn1QESo1PFKevDzgUMkjZQ0CDgXuKu4\ngaSDJSn/ekz+wuvy+24ClkTEtSWCbNifb37zm5nH4Pizj6M/xt/IsTdD/JUo24OPiM2SLgGmAy3A\nTRGxVNLE/PEbgE8BF0rqAN4EPps/fTxwAbBIUnt+35UR8buKojMzs+28+GJl7Sop0RAR00iDp8X7\nbih6/S/Av3Rx3mx8M5WZWU3dfntl7Zx8d1Bra2vWIewQx5+tRo6/kWOHxo5/wYLK2pW90am3SYqs\nYzAzayQf/jAsWiSizCCrE7yZWQPZtAn23BP+9KfyCd4lGjOzBrJ4MRx6aGVtneDNzBpIezuMHl1Z\nWyd4M7MGsmABjBlTWVsneDOzBlJNgvcgq5lZg+jogD32SDc6DR3qQVYzs6axbBmMGAFDhlTW3gne\nzKxBtLdXXp4BJ3gzs4axYEHlM2jACd7MrGFUM8AKHmQ1M2sIW7bAsGHwzDPpTlbJg6xmZk3h6adT\ngt9zz8rPcYI3M2sA1ZZnwAnezKwhVDuDBip84IeZmVVu40Y46yx48820PWAA/PjHcPTRpc+7/Xb4\nwQ+6PrZ8Odx6a3VxOMGbmdXYQw/B+vUweXLa/tWv4Kc/LZ/gr70WLr4Yjjpq+2MtLV3vL8WzaMzM\nauyqq2DzZrj66rS9bBmccQY8+yyom3kvL72UlgF+6SUYNKj8NWo2i0ZSm6Rlkp6SdHkXxydIWiip\nXdKjkk4tOnazpBclLa7kWmZmjS6Xg+InAn7gA/Cud8GiRd2fM20anH56Zcm9UmUTvKQW4DqgDfgg\ncJ6kwzs1mxERH46I0cDngRuLjv0sf66ZWdP705/SjJfjj9+6T0o1+alTuz9v6tTUppYq6cGPBZZH\nxMqI6ABuAyYUN4iIDUWbQ4BXio7NAl6rQaxmZnXv4YfhyCO3XxCsVIJ/6y34wx/gzDNrG0slCX44\nsLpoe01+3zYknSNpKTANuLQ24ZmZNZbO5ZmCk06CpUtTjb2zWbPgsMNg331rG0slCb6iEdCIuDMi\nDgfOBn6+Q1GZmTWomTPh5JO33z9oUKqxT5u2/bHeKM9AZdMk1wIjirZHkHrxXYqIWZIGStorItZV\nEsSkSZPeed3a2kprV19/ZmZ1btMmmD8fxo/v+nihTPO5z23dFwG//S3ccUfp987lcuRyuariKTtN\nUtJA4EngNOA5YB5wXkQsLWpzMLAiIkLSGOC/I+LgouMjgd9GxJFdvL+nSZpZU8jl4IorUh2+Ky++\nmGbUFE+FfPJJOO00WL26+ymUXanJNMmI2AxcAkwHlgC3R8RSSRMlTcw3+xSwWFI78EPgs0VB/BJ4\nEDhU0mpJf135RzAzaxzdlWcK9t031dpnzdq6r1CeqSa5V8o3OpmZ1cgpp8DXvlZ6Nsw//3Pq4Z99\ndtqePDntq7YGX0kP3gnezKwGNm2CvfeG556D3Xbrvt3zz6eEvnlz2t59d/inf4Kdd67uek7wZmZ9\n5IEH4KtfhXnz+uZ6fuCHmVkfmTt327tX64ETvJlZDTz9NBxySNZRbMsJ3sysBp55Bt73vqyj2JYT\nvJlZDaxYAQcdlHUU2/Igq5nZDnr7bdh11/SQj2pnw/SUB1nNzPrA2rWwzz59l9wr5QRvZraDVqyo\nv/o7OMGbme2weqy/gxO8mdkOq8cZNOAEb2a2w1yiMTNrUi7RmJk1Kffgzcya0IYN8Mc/wn77ZR3J\n9pzgzcx2wDPPwMiRMKAOs2kdhmRm1jjqtTwDTvBmZjukXqdIghO8mdkOaegevKQ2ScskPSXp8i6O\nT5C0UFK7pEclnVrpuWZmja5ep0gCDCx1UFILcB1wOrAWeETSXRGxtKjZjIj4f/n2RwK/Ad5f4blm\nZg2tkUs0Y4HlEbEyIjqA24AJxQ0iYkPR5hDglUrPNTNrZBEpwddrD75cgh8OrC7aXpPftw1J50ha\nCkwDLq3mXDOzRvXii2kd+KFDs46kayVLNEBFT+KIiDuBOyWdCPxc0mHVBDFp0qR3Xre2ttLa2lrN\n6WZmmejL8kwulyOXy1V1TrkEvxYYUbQ9gtQT71JEzJI0ENgz366ic4sTvJlZo+jLGTSdO7/f+ta3\nyp5TrkQzHzhE0khJg4BzgbuKG0g6WJLyr8cARMS6Ss41M2tk9TxFEsr04CNis6RLgOlAC3BTRCyV\nNDF//AbgU8CFkjqAN4HPljq39z6KmVnfWrECxo/POoru+aHbZmY91NoKV10Fp55atmnN+aHbZma9\naNUqOPDArKPonnvwZmY9EAG77AKvvZb+7GvuwZuZ9ZJXXoHBg7NJ7pVygjcz64G1a2F4nd+66QRv\nZtYDa9bAAQdkHUVpTvBmZj3gHryZWZNygjcza1Jr17pEY2bWlNascQ/ezKwpuURjZtaknODNzJrQ\nhg2waRPsuWfWkZTmBG9mVqVC710lFwrInhO8mVmVGqE8A07wZmZVc4I3M2tSjbBMATjBm5lVzT14\nM7MmcOml8OMfb7vPCd7MrAn84Q8wY8a2+5omwUtqk7RM0lOSLu/i+PmSFkpaJGmOpFFFxy6TtFjS\n45Iuq3XwZma96eWXYeVKmD07PcGpoClq8JJagOuANuCDwHmSDu/UbAVwUkSMAr4N3Jg/9wjgYuBY\n4MPAWZIOrm34Zma9Z84cOPlkGDIEnnwy7du8OSX+/fbLNrZKlOvBjwWWR8TKiOgAbgMmFDeIiIci\nYn1+cy5Q+F47HJgbEZsi4m1gJvC/ahe6mVnvmjULTjwx/cyalfa9+CLstRfstFO2sVWiXIIfDqwu\n2l6T39edLwD35F8/DpwoaU9JuwKfYGvyNzOre7NnwwknpJ9Cgm+U+jvAwDLHo8zxd0g6BbgIGA8Q\nEUslXQP8HtgAtANbujp30qRJ77xubW2ltbW10suamfWKDRvg8cdh7Fh497vhe99L+7Oqv+dyOXK5\nXFXnKKL7HC5pHDApItry21cCWyLimk7tRgFTgLaIWN7Ne10NPBsR/7fT/igVg5lZFu67D77+dXjw\nwTTAus8+sHAhTJkCy5ZtP3Wyr0kiIkquhlOuRDMfOETSSEmDgHOBuzpd5EBScr+gc3KX9O6iNp8E\nflHdRzAzy0ah/g5pUbETTkglm6Yp0UTEZkmXANOBFuCmfOllYv74DcBVwDBgstLSah0RMTb/FndI\n2gvoAL4YEX/spc9hZlZTs2enm5wKCnX49evh8M5zCetUyRJNnwTgEo2Z1ZnNm9Na7ytXbl3zfe5c\nmDgRhg1LpZvTTss0xIpKNOUGWc3M+p3HHoMDD9z2gR6jR8PTT6c58Y1SovFSBWZmnRTX3wsGDYJj\nj4UXXnCCNzNrWIX5752dcAIMHZp+GoETvJlZkeefh5kzt+/BA5x0UirdNAoneDOzvKefTr30r3yl\n60R+2mlw9919H1dPeRaNmRmwaBGceSZ84xvwt3+bdTTleRaNmVkFli2DM86A//gPOPfcrKOpHZdo\nzKzfmzEDJkxoruQOTvBmZqxaBQcdlHUUtecEb2b93qpV8N73Zh1F7TnBm1m/5wRvZtakmjXBe5qk\nmfVrmzbB7rvDxo3Q0pJ1NJWrxXrwZmZN7dln09oyjZTcK+UEb2b9WrOWZ8AJ3sz6OSd4M7Mm5QRv\nZtaknODNzJpUv07wktokLZP0lKTLuzh+vqSFkhZJmiNpVNGxKyU9IWmxpF9IeletP4CZ2Y7otwle\nUgtwHdAGfBA4T1Ln54mvAE6KiFHAt4Eb8+eOBP4GGBMRRwItwGdrGbyZ2Y7YvBmeew5GjMg6kt5R\nrgc/FlgeESsjogO4DZhQ3CAiHoqI9fnNucAB+dd/BDqAXSUNBHYF1tYscjOzHfTcc7DPPvCuJq0t\nlEvww4HVRdtr8vu68wXgHoCIeBX4AfAs8BzwekTM6HmoZma11czlGSj/wI+K1xCQdApwETA+v30w\n8GVgJLAe+G9J50fErZ3PnTRp0juvW1tbaW1trfSyZmY91kgJPpfLkcvlqjqn5Fo0ksYBkyKiLb99\nJbAlIq7p1G4UMAVoi4jl+X3nAmdExMX57b8CxkXE33c612vRmFkmvvMdeOMN+N73so6kerVYi2Y+\ncIikkZIGAecCd3W6yIGk5H5BIbnnLQPGSdpFkoDTgSXVfggzs97SSD34niiZ4CNiM3AJMJ2UnG+P\niKWSJkqamG92FTAMmCypXdK8/LkLgf8kfUksyre9sRc+g5lZjzR7gvdywWbWbx12GPz61/ChD2Ud\nSfUqKdE4wZtZvxQBgwfDSy/BkCFZR1M9rwdvZtaNl1+GXXZpzOReKSd4M+uXmr3+Dk7wZtZPOcGb\nmTUpJ3gzsybVHxJ8uaUKzMwa1vTpcP31XR+bPx9+9KO+jaeveZqkmTWtCy+EvfaCrpa3GjAAzjgD\ndt65z8OqiUqmSboHb2ZNq70dbrkFxozJOpJsuAdvZk1p48bUe3/99eZc7903OplZv7V4cVqKoBmT\ne6Wc4M2sKbW399/STIETvJk1pQULYPTorKPIlhO8mTUl9+A9yGpmTaijA3bfPS0oNnhw1tH0Dg+y\nmlm/tGRJuku1WZN7pZzgzazpuDyTOMGbWdNZsMAJHpzgzawJtbd7Bg1UkOAltUlaJukpSZd3cfx8\nSQslLZI0R9Ko/P4P5B/CXfhZL+nS3vgQZmYFW7bAY485wUOZtWgktQDXAacDa4FHJN0VEUuLmq0A\nToqI9ZLagBuBcRHxJDA6/z4D8uf/phc+g5nZO556CvbZB4YNyzqS7JXrwY8FlkfEyojoAG4DJhQ3\niIiHImJ9fnMucEAX73M68HRErN7RgM3MSnF5ZqtyCX44UJyU1+T3decLwD1d7P8s8IvqQjMzq54H\nWLcqt1xwxXcgSToFuAgY32n/IOBsYLv6fcGkSZPeed3a2kprV4s3m5l1oaMDHnww1d4BHngAvvnN\nbGPqDblcjlwuV9U5Je9klTQOmBQRbfntK4EtEXFNp3ajgClAW0Qs73RsAvB3hffo4hq+k9XMemzK\nFPjiF+Hww9P2oEFw223NX4OvxQM/5gOHSBoJPAecC5zX6SIHkpL7BZ2Te955wC8rjNnMrCoPPABf\n/jJccUXWkdSfkjX4iNgMXAJMB5YAt0fEUkkTJU3MN7sKGAZMzk+HnFc4X9Jg0gDrlF6J3sz6vdmz\n4YQTso6iPnmxMTNrWG+8Ae95D6xb1/8e7OHFxsysqT38cJox09+Se6Wc4M2sYc2a5fJMKU7wZtaw\nZs+GE0/MOor65Rq8mTWkt96CvfaC1athjz2yjqbvuQZvZk2rvR0OPrh/JvdKOcGbWUNy/b08J3gz\na0iuv5fnGryZNZyItCTwwoUwvNTyh03MNXgza0rLlsFuu/Xf5F4pJ3gzazguz1Sm3GJjZmaZuP12\nuPjiro/9+c/wk5/0bTyNyDV4M6tLF14IY8fC5z7X9fEhQ0AlK9DNrZIavBO8mdWlgw6Ce+7Zus67\nbcuDrGbWkNasSStFHnZY1pE0Nid4M6s7hTXe+3MJphac4M2s7niWTG04wZtZ3fEyBLXhQVYzqyuv\nvw4jRsCrr8JOO2UdTf2qxUO3+9S0aTBzJjz7LKxalf6iC4YPTyPqA+sqYjOrtQcfhGOPdXKvhbIl\nGkltkpZJekrS5V0cP1/SQkmLJM2RNKro2B6S7pC0VNISSeNKXeurX01rPH/84/C976UbHX71q/Tz\n4ovw0EM9+5Bm1jhmzXL9vVZK9ocltQDXAacDa4FHJN0VEUuLmq0AToqI9ZLagBuBQiL/IXBPRPyl\npIHA4FLXW7cuJfn999/+2IQJMHWq/+LNmt3s2XDVVVlH0RzK9eDHAssjYmVEdAC3AROKG0TEQxGx\nPr85FzgAQNLuwIkRcXO+3eaidtuJSDW3vfbq+vhZZ6UEb2bNa9MmWLAAxpX8Xd8qVS7BDwdWF22v\nye/rzheAe/KvDwJelvQzSQsk/UTSrt2d+MYbMGhQ909HP+YYeOUVWLGiTMRm1rDmz093rg4dmnUk\nzaHckGXF01sknQJcBIwveu8xwCUR8Yika4ErgO1++Zo0aRKvvQYtLZDLtdLa2rrd+w8YAJ/4BNx9\nN3zpS5VGZWaNxPX37uVyOXK5XFXnlJwmmR8UnRQRbfntK4EtEXFNp3ajgClAW0Qsz+/bD3goIg7K\nb58AXBERZ3U6NyKCRx9NK8e1t3cf7JQpcMMNMH16VZ/RzOrU66/DP/wDdHSk7Tlz4Pvfh099Ktu4\nGkEt1qKZDxwiaaSkQcC5wF2dLnIgKblfUEjuABHxArBa0qH5XacDT3R3oXXruq+/F5xxRppC9cYb\nZaI2s4bwhz/A0qXQ1pZ+rr46/aZutVGyRBMRmyVdAkwHWoCbImKppIn54zeQSi7DgMlKC0d0RMTY\n/Ft8Cbg1/+XwNPDX3V2rkgQ/dGgafJkxAz75yUo+npnVs1wOPv1puOCCrCNpTnVzJ+t118GSJXD9\n9aXb//CHsGgR3HRT38RnZr3niCPgllvg6KOzjqTxNNRywZX04CFNl7z7btiypfdjMrPe89JLaVng\no47KOpLmVTc3/q9bB+97X/l2Bx8Mw4bB17+enqpejQMP9OCNWb144IG0oFhLS9aRNK+6SvDHHltZ\n2x/+MK1b8+yzlb9/BHzjG3D22Wm+vZlla+ZMOPnkrKNobnWV4Csp0QB89KPpp1r33Zfq98ccU/25\nZlZbuRzcfHPWUTS3hqvB74hx4+Dhh3v3GmZW3iuvpN/AR4/OOpLm5gRvZn3ugQdg/Hgv/93bnODN\nrM/lctDFiiRWY3WR4Ds6YMMG2H333r3OYYelXw1ffrl3r2NmpTnB9426SPCvvpqmPg7o5WgGDEgz\ndebO7d3rmFn31q1LT2wbMybrSJpfXVTASq0DX2vjxqUEf9ZZ5dta4/vzn7tfgtq69+c/p6nFveHe\ne+H4411/7wt18Z+4L+rvBePGwbXX9s21LFtvvAHvfW96hsAee2QdTeOYN6/3E/APftB7721b9bsE\n/5GPpH/Ab7/tO+ia3dy58Npr6RFw/o2tcvfem5658O//nnUktqPqogbflwl+773h3e+GZcv65nqW\nnTlzYPDgdMekVW72bD90o1n0uwQPni7ZX8yZA3/3d2nGhlXm7bfTMxfGjy/f1upf3ST4Pffsu+s5\nwTe/t99OJZovfzk9UGJ9t497t2JPPJF+w91336wjsVqomwTvHrzV0uLF8J73wPDhaWrsnDlZR9QY\nZs1KKzxac+iXCX7UqDSz4o9/7LtrWt+aM2drmaG11WWaSrn+3lz6ZYLfaae0yNEjj/TdNa1vOcFX\nL8I9+GZTNsFLapO0TNJTki7v4vj5khZKWiRpjqRRRcdW5ve3S5rX3TX6OsEDHHccPPRQ317T+k5x\ngv/IR9LjIP0bW2mrVqWxi4MPzjoSq5WSCV5SC3Ad0AZ8EDhP0uGdmq0AToqIUcC3gRuLjgXQGhGj\nix7EvZ0sEvzxx6fZAtZ81qyBjRvh0EPT9s47p2cAuA5fWqH3rpJP+bRGUq4HPxZYHhErI6IDuA2Y\nUNwgIh6KiMIchbnAAZ3eo+w/l6x68A8/vGPPdl21qvdu57aemzMnfYEXJyqXacpz/b35lEvww4HV\nRdtr8vu68wXgnqLtAGZImi/pb7o7aaedUi+rL+23X1q98n/+p2fnv/22Fy6rV13N4z75ZN/wVI7r\n782nXIKvuH8q6RTgIqC4Tj8+IkYDZwJ/L6nL/kFf994LdqRMM29eWnb40UdrG5PtuOL6e8G4cfD4\n42l9GtveK6/A2rVphpk1j3Jr0awFRhRtjyD14reRH1j9CdAWEa8V9kfE8/k/X5b0G1LJZ1bn8996\naxKTJqXXra2ttPbRQtGFgdaLLqr+3KlT0zzrBQtqH5f13Jtvphubjj562/277JL2feUrcEDnImIZ\nw4fDxRdXd868eenB8H3lnHPgwx/u+flz5qQvQa/wWL9yuRy5KuuMihJFZEkDgSeB04DngHnAeRGx\ntKjNgcB9wAUR8XDR/l2Bloh4Q9Jg4PfAtyLi952uEaedFsyYUVXcNbFgAVx4YerZVWvUqPTFcMst\n0N5euu2DD9bHF8GQIXDBBTv2P/Gdd6ZBzHq1alVKVl39ZjZzZnrwerX+9V9TD3fXXSs/5zOfSc8f\n+MAHqr9etRYvTuXGn/1s+2Pf+U56PF45zzyT/l/4+tdrH5/1DklERMkxzpIJPv8mZwLXAi3ATRHx\nXUkTASLiBkk/BT4JPJs/pSMixkp6HzAlv28gcGtEfLeL94/PfCa4/fZqPlptbN6cHjSyenV1y8mu\nWpVmZaxcCfvsk1YsLLXm+JgxcPjh6VpZmjYNrr8ePvaxnp2/Zg0ccUT6kqhnH/sYnH127d5vzBiY\nPDlNt6zU+98Pv/1t+nvvbQsXwrnnbr+A3pYtadmByZNht93Kv89xx1XWzupDJQm+bF8uIqYB0zrt\nu6Ho9cXAdr/ARsQK4KhKAs2qBj9wYErUc+dWl/TuvhvOPDOtVPj+96f1O7p7Os1rr8FTT6UZO4MG\n1SbunrrmmlRa6mmCv/vutOzuddfVNq56N3p0+i2t0gS/fj288MLWaZq97Ygj4Pnnt5+NtmxZStif\n/nTfxGH1py7uZM0qwUPqtVQ70Dp16tb1xUePLl1+mTUrXSPr5A7wiU+kXmVPp3YWf+7+ZMyY6kps\njz2WSnh99byBlpY0o6vz+kqe9mhO8FXe0bphQ0raH/1o2h4zpnQNfubMNEWvHnzoQ+nPJ56o/tyN\nG9Nn6Wnvv5FVm+AXLOj754129e/Y0x7NCf64rU94qsS996beUqFmX64HX09Pj5dSD3zq1OrPvf/+\nlLSyHkfIwqhRaWZOR0dl7dvb+z7BdzXld9Ys9+D7u36f4AtPeFqypLL2ncsURx2VZjF09QXx+uvp\nRqpjj61NrLXQ0wTfX8szkMZa3vveyv+NLFiQvvj70rhxMH9+mjgAaeLAhg19M4vH6le/T/CQej+V\nlGkitk90u+0G++8PTz65ffvZs9PAXD3U3wtaW2HRojTtr1Jdfe7+ptIyzcaNaSnqQjmsrwwblub3\nF6b8zp7tdWWsThJ8Xz7NqSsnnghf+1paRa/Uz0EHpVg7z47orkxTT+WZgp13hlNPhd/9rvJzFi1K\n00D7c2+w3FhLweLFaWpkFl/qxRMGCgne+re6uG8t6x785z9feSLeZ5/t9xV6d53nh8+cCf/2bzsa\nXe0VyjSVzmcv9N77c29w9GiYMqV8uyzKMwXHHZf+zX3xi6n+/rnPZROH1Y+6SPDV3GTUG1padmwN\n7DFj4Oqrt923fn2ahzy220WSs/Pxj8M//mMaNNxpp/Ltp06Fb3+79+OqZ6NHpxuK3n679PTHLGbQ\nFBx/PHz3u+nei2eeye6LxupHXZRoBtRFFD1XuBGmeOnh2bNTci91h2tW9t8/faHdfz/86U+lf1av\nTjNITjop66iztcceaTB++fLS7bKYQVNw2GHw6qtpOYmxYyv78rbmVhc9+Ea3995psPWZZ7b+JlBP\n89+78vnPw4QJZZu907aeBoqzUhhr6W4s4q230kybrFZkHDAgzab5/vfhL/8ymxisvjR437l+dJ5l\nUY8DrMUuuaR8773wM3ly1tHWh3IzaZYsSQPx1SxKVmvHHZfi8ACrgRN8zRxzDFx2WZoXf9RRqaxR\nj/V367lyM2myLM8UHHdcGiMYNy7bOKw+uERTI//7f287T3yPPfr+KVXWuwolmoiuZxRlOYOm4IQT\n4Ec/gqFDs43D6kPZ5YJ7PQApso7BrFLDh8P06TBixPbH2trSLJZ6Ls1Z86jJcsFmttXHPrb94wAL\nBg/OvgdvVsw9eDOzBlRJD96DrGZmTcoJ3sysSTnBm5k1qbIJXlKbpGWSnpJ0eRfHz5e0UNIiSXMk\njep0vEVSu6Tf1jJwMzMrrWSCl9QCXAe0AR8EzpPU+TnxK4CTImIU8G3gxk7HLwOWAE05kprL5bIO\nYYc4/mw1cvyNHDs0fvyVKNeDHwssj4iVEdEB3AZss4JJRDwUEevzm3OBAwrHJB0AfBz4KdCUi802\n+j8Sx5+tRo6/kWOHxo+/EuUS/HBgddH2mvy+7nwBuKdo+9+BfwS2dN3czMx6S7kEX3FZRdIpwEXA\n5fnts4CXIqKdJu29m5nVs5I3OkkaB0yKiLb89pXAloi4plO7UcAUoC0iluf3XQ38FbAZ2BnYDfh1\nRFzY6dymrM2bmfW2cjc6lUvwA4EngdOA54B5wHkRsbSozYHAfcAFEfFwN+9zMvDViDi76k9gZmY9\nUnItmojYLOkSYDrQAtwUEUslTcwfvwG4ChgGTFZaYq8jIrpaKNc9dTOzPpT5WjRmZtY7Mr2TtdxN\nVPVM0s2SXpS0OOtYekLSCEn3S3pC0uOSLs06pkpJ2lnSXEmPSVoi6btZx9QTjXwToKSV+Zsb2yXN\nyzqeakkhulL1AAACwklEQVTaQ9Idkpbm/w01zCNSJH0g/9+98LO+u/9/M+vB52+iehI4HVgLPEKn\n+n49k3Qi8CbwnxFxZNbxVEvSfsB+EfGYpCHAo8A5DfTff9eI2JgfJ5pNGuOZnXVc1ZD0FeBoYGhE\n/EXW8VRD0jPA0RHxatax9ISkW4CZEXFz/t/Q4KL7eRqGpAGk/Dk2IlZ3Pp5lD77sTVT1LCJmAa9l\nHUdPRcQLEfFY/vWbwFJg/2yjqlxEbMy/HEQaH2qoRNMkNwE2ZNySdgdOjIibIY01NmJyzzsdeLqr\n5A7ZJvhqb6KyXiJpJDCadCdyQ5A0QNJjwIvA/RGxJOuYqtToNwEGMEPSfEl/k3UwVToIeFnSzyQt\nkPQTSRk+Kn2HfBb4RXcHs0zwHt2tA/nyzB3AZfmefEOIiC0RcRRpaYyTJLVmHFLFmuQmwPERMRo4\nE/j7fMmyUQwExgDXR8QYYANwRbYhVU/SIOBs4L+7a5Nlgl8LFD/ZcgSpF299RNJOwK+B/4qIO7OO\npyfyv1rfDRyTdSxVOB74i3wd+5fAqZL+M+OYqhIRz+f/fBn4Dank2ijWAGsi4pH89h2khN9ozgQe\nzf8ddCnLBD8fOETSyPw30bnAXRnG068o3bRwE7AkIq7NOp5qSNpb0h7517sAZwDt2UZVuYj4PxEx\nIiIOIv2KfV/nO7zrmaRdJQ3Nvx4MfBRomNlkEfECsFrSofldpwNPZBhST51H6iB0K7OHbnd3E1VW\n8VRL0i+Bk4G9JK0GroqIn2UcVjXGAxcAiyQVkuOVEfG7DGOq1HuAW/IzCAYAP4+IezOOaUc0Wrly\nX+A3+RsbBwK3RsTvsw2pal8Cbs13Lp8G/jrjeKqS/2I9HSg5/uEbnczMmpQf2Wdm1qSc4M3MmpQT\nvJlZk3KCNzNrUk7wZmZNygnezKxJOcGbmTUpJ3gzsyb1/wG+9fXHD85HRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10df62250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Vamos a utilizar validación cruzada: 10-CV\n",
    "grp_cv = 10\n",
    "ind_grp = np.dot(np.random.multinomial(1,[1.0/grp_cv]*grp_cv, nobs),np.arange(1,grp_cv+1).reshape(grp_cv,1))\n",
    "CVRes = np.zeros((L,grp_cv))\n",
    "for g in range(grp_cv):\n",
    "    # asignemos en entrenamiento y validación\n",
    "    ind_train = (ind_grp!=g+1).flatten()\n",
    "    n_train   = ind_train.sum()\n",
    "    ind_valid = (ind_grp==g+1).flatten()\n",
    "    n_valid   = ind_valid.sum()\n",
    "    # y,X\n",
    "    yy_train = yy_std[ind_train]\n",
    "    yy_valid = yy_std[ind_valid]\n",
    "    xx_train = prexmat_std[ind_train,:]\n",
    "    xx_valid = prexmat_std[ind_valid,:]\n",
    "    # Para cada grupo tenemos que estimar todo el path de Lasso\n",
    "    for l in range(L):\n",
    "        # Lasso con la muestra de entrenamiento\n",
    "        lambvec = grid_lamb[l]\n",
    "        clf = linear_model.Lasso(alpha=lambvec)\n",
    "        clf.fit(xx_train,yy_train)\n",
    "        # Ahora obtengamos yy_est con la muestra de validación\n",
    "        pest_valid =  np.dot(xx_valid,(clf.coef_).reshape((nvar,1)))\n",
    "        # queremos que la clasificación entre clientes malos (yy=1) y cero coincida con el promedio muestral\n",
    "        # obtengamos el percentil 1-yy_mean\n",
    "        ptile = np.percentile(pest_valid,100*(1-yy.mean()))\n",
    "        #ptile = 0\n",
    "        # listos para clasificar\n",
    "        clas_est = yy_std.max()*(pest_valid>=ptile)  + yy_std.min()*(pest_valid<ptile)\n",
    "        # obtengamos el error de clasificación\n",
    "        misclas_err = np.not_equal(clas_est,yy_valid).mean()\n",
    "        # guardemos\n",
    "        CVRes[l,g] = misclas_err\n",
    "# veamos qué obtuvimos\n",
    "pd.DataFrame(CVRes, index=grid_lamb).mean(axis=1).plot()\n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Lasso_CV</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>antiguedad</th>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>edad_clie</th>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>factu_antes</th>\n",
       "      <td>-0.000012</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>payratio</th>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>saldo_antes</th>\n",
       "      <td>0.000009</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>uso_antes</th>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tktprom_ant</th>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_last</th>\n",
       "      <td>-0.001181</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_delta</th>\n",
       "      <td>-0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Lasso_CV\n",
       "antiguedad  -0.000000\n",
       "edad_clie   -0.000000\n",
       "factu_antes -0.000012\n",
       "payratio    -0.000000\n",
       "saldo_antes  0.000009\n",
       "uso_antes   -0.000000\n",
       "tktprom_ant  0.000000\n",
       "score_last  -0.001181\n",
       "score_delta -0.000000"
      ]
     },
     "execution_count": 218,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Veamos los lambdas asociados: escojamos el primero\n",
    "lambda_argmin = grid_lamb[pd.DataFrame(CVRes).mean(axis=1).values==pd.DataFrame(CVRes).mean(axis=1).min()][0]\n",
    "# Obtengamos los coeficientes estimados con este lambd\n",
    "clf = linear_model.Lasso(alpha=lambda_argmin)\n",
    "clf.fit(prexmat_std,yy_std)\n",
    "pd.DataFrame(clf.coef_, index = labs_x, columns=['Lasso_CV'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Encontramos que las variables relevantes son **facturación promedio anterior**, **Saldo Anterior**, **Último Score**.\n",
    "Clientes que facturaron más, tuvieron un menor saldo promedio en los últimos seis meses y un mayor score crediticio tienen una menor probabilidad de dejar de pagar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Punto 3: Función para estimar un logit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 280,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mylogit_fn(Ymat,Xmat):\n",
    "    # before defining attributes clean data\n",
    "    xmat    = np.asarray(Xmat)\n",
    "    K       = xmat.shape[1]\n",
    "    T       = xmat.shape[0]\n",
    "    ymat    = np.asarray(Ymat)\n",
    "    # vamos a limpiar los datos  de missing\n",
    "    predata = np.asarray(pd.DataFrame(np.concatenate((ymat,xmat), axis=1)).dropna())\n",
    "    xmat    = predata[:,1:]\n",
    "    nvar       = xmat.shape[1]\n",
    "    nobs       = xmat.shape[0]\n",
    "    ymat    = predata[:,0].reshape((nobs,1))\n",
    "\n",
    "    #----------------------------------------------\n",
    "    # NEWTON RAPHSON:\n",
    "    # Inicialicemos con OLS.  Están mal, pero tienen mismos signos\n",
    "    olsini  = ols.ols_dan(ymat,xmat)\n",
    "    beta_old = np.asarray(olsini.betahat()).reshape((olsini.nvar,1))\n",
    "    #beta_old = np.zeros((nvar,1))\n",
    "    # Listos para empezar el algoritmo\n",
    "    dist = 1\n",
    "    maxIter = 10000\n",
    "    counter = 0\n",
    "    while dist>0.0001 and counter<maxIter:\n",
    "        # Necesitamos la diagonal de W\n",
    "        xbeta = np.dot(xmat,beta_old)\n",
    "        Pmat  = np.divide(np.exp(xbeta), 1+np.exp(xbeta))\n",
    "        PreW  = np.multiply(Pmat,1-Pmat)\n",
    "        # Eliminemos ceros porque hacen que el problema sea inestable\n",
    "        PreW[PreW==0] = 0.0000001\n",
    "        # Y-p\n",
    "        Y_minusP = ymat - Pmat\n",
    "        # Usemos broadcasting para obtener el vector que queremos\n",
    "        w_inv  = 1.0/PreW\n",
    "        z      = xbeta + np.multiply(w_inv, Y_minusP)\n",
    "        # producto cruzado ponderado\n",
    "        # Listos: para invertir la matrix vamos a usar la descomposición de Cholesky \n",
    "        # Hace que el problema sea más estable\n",
    "        premat = np.dot(xmat.T, PreW*xmat)\n",
    "        newbeta = np.dot(np.linalg.inv(premat),np.dot(xmat.T,np.multiply(PreW,z)))\n",
    "        dist = np.max(np.abs(beta_old-newbeta))\n",
    "        # update beta_old\n",
    "        beta_old = newbeta\n",
    "        counter +=1\n",
    "    #return beta_old, counter\n",
    "    return beta_old"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4\n"
     ]
    }
   ],
   "source": [
    "# Simulemos un modelo\n",
    "np.random.seed(250679)\n",
    "N = 200\n",
    "\n",
    "# dos regresores independientes entre sí\n",
    "x0 = np.ones((N,1))\n",
    "x1 = 2*np.random.randn(N,1)\n",
    "x2 = np.random.randn(N,1)\n",
    "xmat = np.hstack((x0,x1,x2))\n",
    "beta_true = np.array([-2,3,5]).reshape((3,1))\n",
    "\n",
    "# vamos a simular residuos de un logit\n",
    "# Paso 1: simular N números aleatorios de una distribución U[0,1]\n",
    "# Paso 2: invertir la CDF de una v.a. logística\n",
    "u = np.random.rand(N,1)\n",
    "eps = np.log(u) - np.log(1-u)\n",
    "# variable dependiente latente\n",
    "y = np.dot(xmat,beta_true) + eps\n",
    "\n",
    "# necesitamos una variable dicótoma:\n",
    "d = 1*(y>0)\n",
    "# vamos la proporción\n",
    "\n",
    "print d.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.09552386],\n",
       "       [ 3.55413408],\n",
       "       [ 5.92810098]])"
      ]
     },
     "execution_count": 285,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# probemos la función\n",
    "mylogit_fn(d,xmat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Punto 4: Logit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 298,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Logit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>constante</th>\n",
       "      <td>9.078780</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>antiguedad</th>\n",
       "      <td>-0.018686</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>edad_clie</th>\n",
       "      <td>0.010848</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>factu_antes</th>\n",
       "      <td>0.000158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>payratio</th>\n",
       "      <td>-0.071825</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>saldo_antes</th>\n",
       "      <td>0.000083</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>uso_antes</th>\n",
       "      <td>-0.378073</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tktprom_ant</th>\n",
       "      <td>-0.000562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_last</th>\n",
       "      <td>-0.014456</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_delta</th>\n",
       "      <td>0.001197</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Logit\n",
       "constante    9.078780\n",
       "antiguedad  -0.018686\n",
       "edad_clie    0.010848\n",
       "factu_antes  0.000158\n",
       "payratio    -0.071825\n",
       "saldo_antes  0.000083\n",
       "uso_antes   -0.378073\n",
       "tktprom_ant -0.000562\n",
       "score_last  -0.014456\n",
       "score_delta  0.001197"
      ]
     },
     "execution_count": 298,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# estimemos un modelo logit para los datos del punto 1\n",
    "beta_logit = pd.DataFrame(mylogit_fn(yy,xmat_cons), index=labs_x_cons, columns=['Logit'])\n",
    "beta_logit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>OLS</th>\n",
       "      <th>Logit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>constante</th>\n",
       "      <td>1.839519</td>\n",
       "      <td>9.078780</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>antiguedad</th>\n",
       "      <td>-0.001630</td>\n",
       "      <td>-0.018686</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>edad_clie</th>\n",
       "      <td>0.001617</td>\n",
       "      <td>0.010848</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>factu_antes</th>\n",
       "      <td>0.000017</td>\n",
       "      <td>0.000158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>payratio</th>\n",
       "      <td>-0.006738</td>\n",
       "      <td>-0.071825</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>saldo_antes</th>\n",
       "      <td>0.000011</td>\n",
       "      <td>0.000083</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>uso_antes</th>\n",
       "      <td>-0.039579</td>\n",
       "      <td>-0.378073</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tktprom_ant</th>\n",
       "      <td>-0.000061</td>\n",
       "      <td>-0.000562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_last</th>\n",
       "      <td>-0.002310</td>\n",
       "      <td>-0.014456</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_delta</th>\n",
       "      <td>0.000067</td>\n",
       "      <td>0.001197</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  OLS     Logit\n",
       "constante    1.839519  9.078780\n",
       "antiguedad  -0.001630 -0.018686\n",
       "edad_clie    0.001617  0.010848\n",
       "factu_antes  0.000017  0.000158\n",
       "payratio    -0.006738 -0.071825\n",
       "saldo_antes  0.000011  0.000083\n",
       "uso_antes   -0.039579 -0.378073\n",
       "tktprom_ant -0.000061 -0.000562\n",
       "score_last  -0.002310 -0.014456\n",
       "score_delta  0.000067  0.001197"
      ]
     },
     "execution_count": 302,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# miremos los del modelo lineal de probabilidad\n",
    "linprob = ols.ols_dan(yy,xmat_cons)\n",
    "beta_linprob = pd.DataFrame(linprob.betahat(), index=labs_x_cons, columns=['OLS'])\n",
    "pd.concat((beta_linprob, beta_logit),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 325,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>OLS</th>\n",
       "      <th>Logit</th>\n",
       "      <th>EfectosMarg_Logit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>constante</th>\n",
       "      <td>1.839519</td>\n",
       "      <td>6.717859</td>\n",
       "      <td>0.946679</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>antiguedad</th>\n",
       "      <td>-0.001630</td>\n",
       "      <td>-0.025802</td>\n",
       "      <td>-0.003636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>edad_clie</th>\n",
       "      <td>0.001617</td>\n",
       "      <td>0.018905</td>\n",
       "      <td>0.002664</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>factu_antes</th>\n",
       "      <td>0.000017</td>\n",
       "      <td>0.000330</td>\n",
       "      <td>0.000046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>payratio</th>\n",
       "      <td>-0.006738</td>\n",
       "      <td>0.012575</td>\n",
       "      <td>0.001772</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>saldo_antes</th>\n",
       "      <td>0.000011</td>\n",
       "      <td>0.000092</td>\n",
       "      <td>0.000013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>uso_antes</th>\n",
       "      <td>-0.039579</td>\n",
       "      <td>-0.567066</td>\n",
       "      <td>-0.079911</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tktprom_ant</th>\n",
       "      <td>-0.000061</td>\n",
       "      <td>-0.001990</td>\n",
       "      <td>-0.000280</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_last</th>\n",
       "      <td>-0.002310</td>\n",
       "      <td>-0.010815</td>\n",
       "      <td>-0.001524</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_delta</th>\n",
       "      <td>0.000067</td>\n",
       "      <td>-0.004071</td>\n",
       "      <td>-0.000574</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  OLS     Logit  EfectosMarg_Logit\n",
       "constante    1.839519  6.717859           0.946679\n",
       "antiguedad  -0.001630 -0.025802          -0.003636\n",
       "edad_clie    0.001617  0.018905           0.002664\n",
       "factu_antes  0.000017  0.000330           0.000046\n",
       "payratio    -0.006738  0.012575           0.001772\n",
       "saldo_antes  0.000011  0.000092           0.000013\n",
       "uso_antes   -0.039579 -0.567066          -0.079911\n",
       "tktprom_ant -0.000061 -0.001990          -0.000280\n",
       "score_last  -0.002310 -0.010815          -0.001524\n",
       "score_delta  0.000067 -0.004071          -0.000574"
      ]
     },
     "execution_count": 325,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# vamos a obtener los efectos marginales: evaluados en las medias\n",
    "mean_x = np.nanmean(xmat_cons,axis=0)\n",
    "xmeanbeta = np.dot(mean_x.reshape((1,nvar+1)),\n",
    "                  beta_logit.values.reshape((nvar+1,1)))[0,0]\n",
    "Pmat  = np.divide(np.exp(xmeanbeta), 1+np.exp(xmeanbeta))\n",
    "margeff = Pmat*(1-Pmat)*beta_logit\n",
    "margeff.columns = ['EfectosMarg_Logit']\n",
    "LogitDF = pd.concat((beta_linprob, beta_logit,margeff),axis=1)\n",
    "LogitDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Una vez tenemos los efectos marginales y los comparamos con los del modelo lineal de probabilidad encontramos que son muy similares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 330,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>OLS</th>\n",
       "      <th>Logit</th>\n",
       "      <th>EfectosMarg_Logit</th>\n",
       "      <th>Error Estándar Coeficientes</th>\n",
       "      <th>Std Efecto Marginales</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>constante</th>\n",
       "      <td>1.839519</td>\n",
       "      <td>6.717859</td>\n",
       "      <td>0.946679</td>\n",
       "      <td>1.833378</td>\n",
       "      <td>0.234230</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>antiguedad</th>\n",
       "      <td>-0.001630</td>\n",
       "      <td>-0.025802</td>\n",
       "      <td>-0.003636</td>\n",
       "      <td>0.007051</td>\n",
       "      <td>0.000798</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>edad_clie</th>\n",
       "      <td>0.001617</td>\n",
       "      <td>0.018905</td>\n",
       "      <td>0.002664</td>\n",
       "      <td>0.014482</td>\n",
       "      <td>0.001858</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>factu_antes</th>\n",
       "      <td>0.000017</td>\n",
       "      <td>0.000330</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>0.000242</td>\n",
       "      <td>0.000029</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>payratio</th>\n",
       "      <td>-0.006738</td>\n",
       "      <td>0.012575</td>\n",
       "      <td>0.001772</td>\n",
       "      <td>0.088490</td>\n",
       "      <td>0.011569</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>saldo_antes</th>\n",
       "      <td>0.000011</td>\n",
       "      <td>0.000092</td>\n",
       "      <td>0.000013</td>\n",
       "      <td>0.000048</td>\n",
       "      <td>0.000006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>uso_antes</th>\n",
       "      <td>-0.039579</td>\n",
       "      <td>-0.567066</td>\n",
       "      <td>-0.079911</td>\n",
       "      <td>0.233824</td>\n",
       "      <td>0.030000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tktprom_ant</th>\n",
       "      <td>-0.000061</td>\n",
       "      <td>-0.001990</td>\n",
       "      <td>-0.000280</td>\n",
       "      <td>0.001040</td>\n",
       "      <td>0.000137</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_last</th>\n",
       "      <td>-0.002310</td>\n",
       "      <td>-0.010815</td>\n",
       "      <td>-0.001524</td>\n",
       "      <td>0.002745</td>\n",
       "      <td>0.000349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_delta</th>\n",
       "      <td>0.000067</td>\n",
       "      <td>-0.004071</td>\n",
       "      <td>-0.000574</td>\n",
       "      <td>0.004856</td>\n",
       "      <td>0.000628</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  OLS     Logit  EfectosMarg_Logit  \\\n",
       "constante    1.839519  6.717859           0.946679   \n",
       "antiguedad  -0.001630 -0.025802          -0.003636   \n",
       "edad_clie    0.001617  0.018905           0.002664   \n",
       "factu_antes  0.000017  0.000330           0.000046   \n",
       "payratio    -0.006738  0.012575           0.001772   \n",
       "saldo_antes  0.000011  0.000092           0.000013   \n",
       "uso_antes   -0.039579 -0.567066          -0.079911   \n",
       "tktprom_ant -0.000061 -0.001990          -0.000280   \n",
       "score_last  -0.002310 -0.010815          -0.001524   \n",
       "score_delta  0.000067 -0.004071          -0.000574   \n",
       "\n",
       "             Error Estándar Coeficientes  Std Efecto Marginales  \n",
       "constante                       1.833378               0.234230  \n",
       "antiguedad                      0.007051               0.000798  \n",
       "edad_clie                       0.014482               0.001858  \n",
       "factu_antes                     0.000242               0.000029  \n",
       "payratio                        0.088490               0.011569  \n",
       "saldo_antes                     0.000048               0.000006  \n",
       "uso_antes                       0.233824               0.030000  \n",
       "tktprom_ant                     0.001040               0.000137  \n",
       "score_last                      0.002745               0.000349  \n",
       "score_delta                     0.004856               0.000628  "
      ]
     },
     "execution_count": 330,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Finalmente, bootstrapping:\n",
    "B = 300\n",
    "BetaBoot   = np.zeros((B,nvar+1))\n",
    "BetaBootME = np.zeros((B,nvar+1))\n",
    "for b in range(B):\n",
    "    ind_b = np.random.choice(nobs,nobs,replace=True)\n",
    "    beta_logit = pd.DataFrame(mylogit_fn(yy[ind_b],xmat_cons[ind_b,:]), index=labs_x_cons, columns=['Logit'])\n",
    "    BetaBoot[b,:] = beta_logit.values.flatten()\n",
    "    # Efectos Marginales\n",
    "    xmeanbeta = np.dot(mean_x.reshape((1,nvar+1)),\n",
    "                      beta_logit.values.reshape((nvar+1,1)))[0,0]\n",
    "    Pmat  = np.divide(np.exp(xmeanbeta), 1+np.exp(xmeanbeta))\n",
    "    margeff = Pmat*(1-Pmat)*beta_logit\n",
    "    BetaBootME[b,:] = margeff.values.flatten()\n",
    "\n",
    "BetaBoot = pd.DataFrame(BetaBoot, columns = labs_x_cons)\n",
    "BetaBootStd = pd.DataFrame(BetaBoot.std())\n",
    "BetaBootStd.columns = [u'Error Estándar Coeficientes']\n",
    "# repitamos para efectos marginales\n",
    "BetaBootME = pd.DataFrame(BetaBootME, columns = labs_x_cons)\n",
    "BetaBootMEStd = pd.DataFrame(BetaBootME.std())\n",
    "BetaBootMEStd.columns = [u'Std Efecto Marginales']\n",
    "#pd.concat((beta_linprob, beta_logit,margeff),axis=1)\n",
    "pd.concat((LogitDF, BetaBootStd, BetaBootMEStd),axis=1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}