chapter13.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 13"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup and Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING *** OLE2 inconsistency: SSCS size is 0 but SSAT size is non-zero\n"
     ]
    }
   ],
   "source": [
    "nhefs_all = pd.read_excel('NHEFS.xls')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1629, 64)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs_all.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 13.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subset the data as in the last chapter (though it's not necessary to restrict on all of these columns)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "restriction_cols = [\n",
    "    'sex', 'age', 'race', 'wt82', 'ht', 'school', 'alcoholpy', 'smokeintensity'\n",
    "]\n",
    "missing = nhefs_all[restriction_cols].isnull().any(axis=1)\n",
    "nhefs = nhefs_all.loc[~missing]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1566, 64)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 13.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>qsmk</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>1163</td>\n",
       "      <td>403</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "qsmk      0    1\n",
       "count  1163  403"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table = nhefs.groupby('qsmk').agg({'qsmk': 'count'}).T\n",
    "table.index = ['count']\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 13.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add squared features as in chapter 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in ['age', 'wt71', 'smokeintensity', 'smokeyrs']:\n",
    "    nhefs['{}^2'.format(col)] = nhefs[col] * nhefs[col]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add constant term. We'll call it `'one'` this time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs['one'] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A column of zeros will also be useful later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs['zero'] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add dummy features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "edu_dummies = pd.get_dummies(nhefs.education, prefix='edu')\n",
    "exercise_dummies = pd.get_dummies(nhefs.exercise, prefix='exercise')\n",
    "active_dummies = pd.get_dummies(nhefs.active, prefix='active')\n",
    "\n",
    "nhefs = pd.concat(\n",
    "    [nhefs, edu_dummies, exercise_dummies, active_dummies],\n",
    "    axis=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add interaction term"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs['qsmk_x_smokeintensity'] = nhefs.qsmk * nhefs.smokeintensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the model to estimate $E[Y|A=a, C=0, L=l]$ for each combination of values of $A$ and $L$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = nhefs.wt82_71\n",
    "X = nhefs[[\n",
    "    'one', 'qsmk', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2',\n",
    "    'qsmk_x_smokeintensity'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>one</th>                   <td>   -1.5882</td> <td>    4.313</td> <td>   -0.368</td> <td> 0.713</td> <td>  -10.048</td> <td>    6.872</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>                  <td>    2.5596</td> <td>    0.809</td> <td>    3.163</td> <td> 0.002</td> <td>    0.972</td> <td>    4.147</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sex</th>                   <td>   -1.4303</td> <td>    0.469</td> <td>   -3.050</td> <td> 0.002</td> <td>   -2.350</td> <td>   -0.510</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>race</th>                  <td>    0.5601</td> <td>    0.582</td> <td>    0.963</td> <td> 0.336</td> <td>   -0.581</td> <td>    1.701</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_2</th>                 <td>    0.7904</td> <td>    0.607</td> <td>    1.302</td> <td> 0.193</td> <td>   -0.400</td> <td>    1.981</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_3</th>                 <td>    0.5563</td> <td>    0.556</td> <td>    1.000</td> <td> 0.317</td> <td>   -0.534</td> <td>    1.647</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_4</th>                 <td>    1.4916</td> <td>    0.832</td> <td>    1.792</td> <td> 0.073</td> <td>   -0.141</td> <td>    3.124</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_5</th>                 <td>   -0.1950</td> <td>    0.741</td> <td>   -0.263</td> <td> 0.793</td> <td>   -1.649</td> <td>    1.259</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>exercise_1</th>            <td>    0.2960</td> <td>    0.535</td> <td>    0.553</td> <td> 0.580</td> <td>   -0.754</td> <td>    1.346</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>exercise_2</th>            <td>    0.3539</td> <td>    0.559</td> <td>    0.633</td> <td> 0.527</td> <td>   -0.742</td> <td>    1.450</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>active_1</th>              <td>   -0.9476</td> <td>    0.410</td> <td>   -2.312</td> <td> 0.021</td> <td>   -1.752</td> <td>   -0.143</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>active_2</th>              <td>   -0.2614</td> <td>    0.685</td> <td>   -0.382</td> <td> 0.703</td> <td>   -1.604</td> <td>    1.081</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>                   <td>    0.3596</td> <td>    0.163</td> <td>    2.202</td> <td> 0.028</td> <td>    0.039</td> <td>    0.680</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age^2</th>                 <td>   -0.0061</td> <td>    0.002</td> <td>   -3.534</td> <td> 0.000</td> <td>   -0.009</td> <td>   -0.003</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wt71</th>                  <td>    0.0455</td> <td>    0.083</td> <td>    0.546</td> <td> 0.585</td> <td>   -0.118</td> <td>    0.209</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wt71^2</th>                <td>   -0.0010</td> <td>    0.001</td> <td>   -1.840</td> <td> 0.066</td> <td>   -0.002</td> <td> 6.39e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeintensity</th>        <td>    0.0491</td> <td>    0.052</td> <td>    0.950</td> <td> 0.342</td> <td>   -0.052</td> <td>    0.151</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeintensity^2</th>      <td>   -0.0010</td> <td>    0.001</td> <td>   -1.056</td> <td> 0.291</td> <td>   -0.003</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeyrs</th>              <td>    0.1344</td> <td>    0.092</td> <td>    1.465</td> <td> 0.143</td> <td>   -0.046</td> <td>    0.314</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeyrs^2</th>            <td>   -0.0019</td> <td>    0.002</td> <td>   -1.209</td> <td> 0.227</td> <td>   -0.005</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk_x_smokeintensity</th> <td>    0.0467</td> <td>    0.035</td> <td>    1.328</td> <td> 0.184</td> <td>   -0.022</td> <td>    0.116</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = sm.OLS(y, X) \n",
    "res = ols.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at an example prediction, for the subject with ID 24770"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "example = X.loc[nhefs.seqn == 24770]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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>one</th>\n",
       "      <th>qsmk</th>\n",
       "      <th>sex</th>\n",
       "      <th>race</th>\n",
       "      <th>edu_2</th>\n",
       "      <th>edu_3</th>\n",
       "      <th>edu_4</th>\n",
       "      <th>edu_5</th>\n",
       "      <th>exercise_1</th>\n",
       "      <th>exercise_2</th>\n",
       "      <th>...</th>\n",
       "      <th>active_2</th>\n",
       "      <th>age</th>\n",
       "      <th>age^2</th>\n",
       "      <th>wt71</th>\n",
       "      <th>wt71^2</th>\n",
       "      <th>smokeintensity</th>\n",
       "      <th>smokeintensity^2</th>\n",
       "      <th>smokeyrs</th>\n",
       "      <th>smokeyrs^2</th>\n",
       "      <th>qsmk_x_smokeintensity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1581</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>26</td>\n",
       "      <td>676</td>\n",
       "      <td>111.58</td>\n",
       "      <td>12450.0964</td>\n",
       "      <td>15</td>\n",
       "      <td>225</td>\n",
       "      <td>12</td>\n",
       "      <td>144</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 21 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      one  qsmk  sex  race  edu_2  edu_3  edu_4  edu_5  exercise_1  \\\n",
       "1581    1     0    0     0      0      0      1      0           1   \n",
       "\n",
       "      exercise_2  ...  active_2  age  age^2    wt71      wt71^2  \\\n",
       "1581           0  ...         0   26    676  111.58  12450.0964   \n",
       "\n",
       "      smokeintensity  smokeintensity^2  smokeyrs  smokeyrs^2  \\\n",
       "1581              15               225        12         144   \n",
       "\n",
       "      qsmk_x_smokeintensity  \n",
       "1581                      0  \n",
       "\n",
       "[1 rows x 21 columns]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1581    0.342157\n",
       "dtype: float64"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.predict(example)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick look at mean and extremes of the predicted weight gain, across subjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimated_gain = res.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   min     mean      max\n",
      "------------------------\n",
      "-10.88     2.64     9.88\n"
     ]
    }
   ],
   "source": [
    "print('   min     mean      max')\n",
    "print('------------------------')\n",
    "print('{:>6.2f}   {:>6.2f}   {:>6.2f}'.format(\n",
    "    estimated_gain.min(),\n",
    "    estimated_gain.mean(),\n",
    "    estimated_gain.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick look at mean and extreme values for actual weight gain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   min     mean      max\n",
      "------------------------\n",
      "-41.28     2.64    48.54\n"
     ]
    }
   ],
   "source": [
    "print('   min     mean      max')\n",
    "print('------------------------')\n",
    "print('{:>6.2f}   {:>6.2f}   {:>6.2f}'.format(\n",
    "    nhefs.wt82_71.min(),\n",
    "    nhefs.wt82_71.mean(),\n",
    "    nhefs.wt82_71.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A scatterplot of predicted vs actual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9WYxkWXrf9zt3vzfW3LPWrullpqdnODMkR1xAwQBFyAYs2xIMyTBAGxRMQNCbbRGwZAM2vLzINuDl0YT4QC+y5UUGZcGGLZukbHGZmZ6Z7uFweq3uWnPP2OPu5xw/nIjorKzMquzuqszqqvMDClUZGXHjxFLfd863/D+htcZisVgslqM4F70Ai8VisTx7WOdgsVgsloewzsFisVgsD2Gdg8VisVgewjoHi8VisTyEd9ELeBKsrq7qGzduXPQyLBaL5QvF97///QOt9dpJv3sunMONGzd48803L3oZFovF8oVCCHH7tN/ZsJLFYrFYHsI6B4vFYrE8hHUOFovFYnkI6xwsFovF8hDWOVgsFovlIZ6LaiWLxfJsMilqtgcZaVmTBB6XujHN0JqdLwL25GCxWJ4Kk6Lmg90xlVQ0Q49KKj7YHTMp6otemuUMWOdgsVieCtuDjNBziHwXIQSR7xJ6DtuD7KKXZjkD9nxnsVieCmlZPxRCCj3nqZ8cbCjryWBPDhaL5amQBB5FrR64ragVSfD0DLUNZT05rHOwWCxPhUvdmKJW5JVEa01eSYpacakbP7XntKGsJ4d1DhaL5anQDD1e22jhuyaU5LsOr220nmqIJy1rQu9BsxZ6DmlpTw6fFhuIs1gsT425gzgv5qGsyHcXtz3tUNbzij05WCyW54aLCGU9r1jnYLFYnhsuIpT1vGLfMYvF8lxx3qGs5xV7crBYLBbLQ1jnYLFYLJaHsM7BYrFYLA9hnYPFYrFYHsI6B4vFYrE8hHUOFovFYnkI6xwsFovF8hDWOVgsFovlIS7cOQghXCHED4UQ/3D285eEEN8RQnwghPh7QojgotdosVgsLxoX7hyAfx1458jP/zHwn2utXwP6wK9fyKosFovlBeZCnYMQ4irwF4C/M/tZAH8O+J9nd/lt4C9dzOosFovlxeWiTw7/BfBvAfNxUSvAQGs9F1+/B1w56YFCiL8mhHhTCPHm/v7+01+pxWKxvEBcmHMQQvxzwJ7W+vtHbz7hrvqkx2utf1Nr/W2t9bfX1taeyhotFovlReUiVVl/CfgXhBD/LBABbcxJoiuE8Ganh6vA1gWu0WKxnDOTomZ7kJGWNUngcakbW8ntC+DCTg5a639ba31Va30D+JeB39Va/yrwe8Bfnt3t14DfuaAlWiyWc2ZS1HywO6aSimboUUnFB7tjJoUd83neXHTO4ST+JvA3hBAfYnIQv3XB67FYLOfE9iAj9ByU1tzupXx8MGV3mHFzb3zRS3vheCbOalrr3wd+f/bvj4Cfu8j1WCyWiyEtaxwBtw9TAs+hEbgUUvLuzphX1u1Et/PkWTw5WCyWF5Qk8Ljfzwg8h9BzEUIgcOjEHtuD7KKX90JhnYPFYnlmuNSNGeY1WoNGU9SKUiqudBPS0uYdzhPrHCwWyzNDM/T4ymYLjSYtJb4ruLHSwHUESWBDSueJfbctFsszxavrLbSG0HMIPYeiVhS14vpK46KX9kJhnYPFYnmmeguaocdrGy22BxmTwqzn+krDJqPPGftuWywvKHOHcDAp2B8XXO5GLCUBRW16C17buLjqoLmDsFwcNudgsbyAHG02y8oaIQQ7w5yskkS+S+g5tjroBceeHCyWF5B5s1nku+S1ohW5lLVib1xwY8Uj9JxTu5KfpRCU5elhTw4WywvE/MTw1t0+20Nj4GPfpaw1geuQlRKAolYnVgdZeYsXB+scLJYXhKOGfaURkJaKWwdTGqFHKRXjQhL5DnklKWrFpW780DWOnjiEEDYE9RxjnYPF8oJw1LCvt2OEAI1gnFdstiOUMqcF33VOTUanZU3oPWg2Qs+xDWrPITZQaLG8IKRlvTD4SeByY6XB3ijjcFpyfbnB65c2H5s7SAKPolZEvru47bQQlOWLjf1ELZYL5rwSvMcNexK4bHZiri03zlw2eqkb88GuUUi1DWrPNzasZLFcIOeZ4L3UjSlqRV5JtNaPzC2cxrz/wHdNNdOjQlCWLzb2E7VYLpCjeQBg8ff2IHviTWBPqvPYNqi9GFjnYLFcIEfzAHMe1WPweXnShv157nl4nl/bWbBhJYvlApnnAY7yRUnwPs89D8/zazsr1jlYLBfIk8gDXBTPc8/D8/zazop1DhbLBfJFTvA+zz0Pz/NrOyvP/jfQYnnO+aImeJ/nnofn+bWdFXtysFgsn4kvckjscTzPr+2svDhu0GKxPFEuaijPeVQR2YFD1jlYLJbPwXmHxOZVRKHn0Ay9pzqY6Isa7ntSWOdgsVwwL3o9/afhPJsGX3TsN9BiuUDOcyf8NDkvB3feTYMvsuO2CWmL5QJ5VD393HG8fbf/TDdgnWfD2Hk2Db7ojXDWOVgsF8hp9fQHk+ILY5jOs2HsPKuIXvRGOOscLJYL5LSd8KSovzCG6Twbxs6zafBFb4R7MYJnFsszymnzEZqhd6Jh+jQnhycZL3/Utc67Yey8qog+6+s6y/v+Rchl2JODxXKBnLYTXm2Gnyu2fpZ4+VlzGo+71llDPV+UHMqczxLC+jTv+7MeMhRa64tew+fm29/+tn7zzTcvehkWyxPjaBXT0RPF5W7MKKseu+OcG5/5rjcta+71UiqleX2zTTv22ZrF1I9e/6QQzfFrAeSVXDiy+XoftROevx6pNKOsZJTXCOAXXlllox09hXfwyfBpd/hnea/Ocp/zQgjxfa31t0/6nT05WCzPICedKC53Y7YG2Zl2nEfj5WlZc+tginDAEVBJxR/fPEAqdaacxqNi73Ojf3PPhMZeWW+d6GC2BxlSaXZGObWCbuLjOA7f+ejwmdsxH2X+OXzz2tKZchtnyVN8UXIZz1aQy2KxLDgeW5+fJM7SAHY0Xr43Lgg8BxAkgXEEGhhmFcuNcPGY03Iap8Xe9ZE1Pa5HIy1rRllJ4DoLw9gKXfpp+YVtYDvpVHGWPEUSePTTinFekVWS2HdpRT6d2L+Il3Eq9uRgsTzjzHfnb93tsz3MHthhnrbjnMfLe9OCm3sTPtqfcOdwSjM0BqgdeYyy6oHHnJbTOC32LmbPf5bTRxJ4jPKawBOL20qpaMf+M7djPgun5Q3asf/YPEU79vlwb8KkkMS+w6SQfLg3oW2dg8ViOStHjdBKIyAtFbcOpguDeppBb4Yel7sx9/sZtVJEnst6O2J3nJOWknYc4DrOmZOtjoCb+xN+sj2ilGqx0z9reORSN0YA49w8X1FLylotdsufJlH9LCS2T+uBGGXVY0ttR1nFq+sNGqFLWioaocur642HnPVFc2FhJSHENeC/BjYBBfym1vq/FEIsA38PuAHcAv4lrXX/otZpsVwkR43Qejvm1uEUrQW7o5xLHbOjX2mGfLA7fihpOjdUr6w3uXUwJfActIZ7/Skb7Ziff3mFUVY9UnX0aGL8jUvtReIaMGGlvTFKQ+y7rLUiHMGpzuoXXlnlOx8d0k9L2rHPZieirDVFbYzoWeRDzio3cpYE+ecpJT1JxkMqzc390eKar6yf/BrSsmYpCR4I6Wmtn7ncy0XmHGrgN7TWPxBCtIDvCyH+EfBXgf9Ha/23hRB/C/hbwN+8wHVaLMDF1KYfNUJJ4HJjpcHeKONwWnJ9ucFKM1xUHR03lvPHCiG4sdpgb1yQFjUaFsb0cZVCpwnd3dwbkxaSaSFpRi5lrXh/d8xGO+Sb15ZOvNZGO+LPfXXjgfdwWtYE7tmF9M4ivPc4B/Ik9KzmjjEtJWlR4wiHQVZyuRPRDD36ack72yPWWiGrzfBC+0I+KxcWVtJab2utfzD79xh4B7gC/EXgt2d3+23gL13MCi2WT7io2vTjHdRJ4LLZifnWrHpmlFWnxv2PPjYJPG6sNHh5rcnrm+0zG8HTKmvu9FK6ic9XNlv4rkutNIlvjO2jrn28+meetzh+/dPyEGep9Hmc7MXnlcWYFDVpITmclOyOcvJKcXN/TFbW1EpzOC3YGeYIIcjK+jP3hVw0z0TOQQhxA/hp4DvAhtZ6G4wDAdZPecxfE0K8KYR4c39//7yWanlBuSidnccZkkcZyydhhE6T95g/z9zpfPVS+zNVHH1aIb3j909L47Tv9NKFAX6cAzlrKelpuY0P98aM8pJhVjHJK2qlcRxBtxHQjn3e3xkTeA6tyCWfnRCOfle+KHPDL3w1Qogm8L8A/4bWeiSEeNxDANBa/ybwm2Ca4J7eCi2W85eKnvO4iWSPClE8iWlm7djnj28eoDEVTiaRLbi2nDyR0Mhx+ZB+WrI1yFlrhYvfH13v0ftLpXh/d4IQgtfWm4sduhA8cm1nCeucFnq63I15b2dMN/YJfYdmGFMpzaVORCU1geswzGquLjmUtSaePcfx78pZJUAuUmbjQk8OQggf4xj+O63135/dvCuEuDT7/SVg76LWZ7HMOU+p6OM8qhHrcaeDT9vEdZRJUXNzb4LWmv1xzvs7Yz7an3C5G/PqeuuJhEaO7qL3xgX3+xmXuzHrrfDE0N3R+986TGmELl/eaNEIvcUOXcAj13aWE9VpJ8Uf3RvQiTyEgNhzEQgC16FS4ArBwaQgL2t+vDV8oHT4s3xXLlpm48KcgzBHhN8C3tFa/2dHfvUPgF+b/fvXgN8577VZLMd5VuPETzNEcXNvzPYwIwk9Xltv8fJ6E9cR7AyzJ/q882slgTHAu6Oc270UpfWJobv5/a8vJ7y2bh43Zx4uOqm7fHuQ8fbdPtsD44AetfbTQk+DtOTKUkxZK5qRR1FLxnnJdj9Fas3bdwe0Yo/Ac+gmATujjN60/EzflYuWDL/IsNIvAf8q8CdCiLdmt/07wN8G/kchxK8Dd4C/ckHrs1gWXNTA+bOEFZ6WSumdXkozcgm9eWjEJXck3/n4EOCJhjkmRc27O2O6iUfDdSml6ed4aSWhkurEx5wlpDa/9vEQ0dasumlevXRSp3M/NRpQWSmJA5d25NFNTH/IvPprnNfc381YagQ0fIdvXFsi9h3W2yGTQjJMS/ppyS++snrq+3TaZ3xRocw5F+YctNb/BDgtwfAr57kWi+UsnPfA+YseIZpXkoOxKX2NPIdG6LE1yBCCJ76e7UFGJ/YQOAghFg7pfj/jlfWT3/PT5M6vrzQeuvZp5a/zaxx/jzuxz4d7UxqhRzM0Xcx7o4Kf+9Iyw1mF2EvLCWUlia91+fJGizu9KUlgynonheTGSgO9nDApHjbycx71GV90yeuFJ6Qtli8S55kgPEtN/+fhUa9lUtQoDb1piVKKUml645LVdsQ3r3YWYY4ntZ60rLnSTbjdSwEIPIHWMMzrU8MxjzvNzV/fW3f7LDcCNtrRwrDOd+Cn9nHsT3h1vbnQP2qGLpc6EVLpB56zUprX1pskgUvsu5S1Rin4uDchrySOgM3O6eGksziu+XpPc35PC+scLJYzct47+bOGFT6Lw3rca9keZFzpxNztTZEafEdQa8UoK1htBSeu5/M4ziQwCdcbKw32ZxIfjoCvbD76vZ07iPlz39wbkwTeA5LkK42A6Ux25MZqY7EjT4LTQzeDtOTl1QbLjU9e67yL+egJcr5ugLVWxLs7Q7aHBc3AxXVgkkvSQp56enjUZ3xRocw51jlYLGfkSe3kjxpRzSex1c8yYe3TOqz5c7+7M8J3BVeXGieeAtKyppSSL2+0mRY1eS3xPAe0ZlJIVpsPrudRpZ8nzZ847kiOGvPrszLZola8ekpI6Si7o5zvfHSIVEbIL3ArfvfdXTqRz0YnohUFpFX+kOzI9ZUG24PsxPe4mwSPVKI9vm6A2HfwHAcXUwLsuy5f2UxwhDj1O3L0M05Lyf44Z5hVNGbv0XmHMo9inYPFckaeRILwqBF1hOCDWanolzc+qdOfG/ZHhRUeMPKO4OpycmqoZ37fg0nB3ijnylKCAASCW4dTbqw0FpVCk8LMaNgbF7x1u4/nChzHQWuFVpq8UuwOc146YsDnRva44zycFPyvPzhgox3Rjn06sZmNPZ9LcTxBPHckn2aXPClq/vjmAY7jsJQEjPKKPz1IqZWkGbhUUrM7ztloRYzzciE7Mr/2ae/xN652F0Z/fns/LREIaqkYZhWjbILrOHz9SgepzKmiqhXXVxPEkXTqo74j8+dPS8n2MAMEroClxOdHdwckoYvgySb/z4p1DhbLGTkuNNcMfQ6nOZXUZ/7Pe9SI3j6czu6v2Z+U3JjFkueG/bSwArCYqrY7zJBac7ef8c1rXVab4UOhnrkzSssax3HYGeU4QiAEBK7D/jjnpZXGoo/jg90xsecyKSqUgnFR04l9hAPXujGjvGJvXLDaDBdG9rjjPJjk/L/v75GVknbk4zkuWZmz2Yn40b0BV7rxQyewuVDgUR4XqtoeZGigFZmeg0le0whdRpliWkquekZ59t3dEa3QVBsdvcajQjeNWXhtfvu812BnmBN4xhmNC8mf3BvwK29sAnDrwAgjNsNPKq42OxGdODj1tby20eKPbh5QS0UnCVhvJQDcPhzTCF1em/WUHD8RPu3812OvJIT4JeDfB16a3V8AWmv98hNbhcXyjDPX05kLzY2yih/eGbDZDvmpq52Hdv2ncdSIZpUkCVzyUi8SmEprslI+8B/+uMGcO4adUY7vOgQOVDW8fW/AL768+oAy6lFnlFeKVuhSSk2tFGWt8D1BWh6Z0SBmMXel+NqVLm/e6qOBspZEgcfdfs7rm00OJgUAt2YOblLUi9GXdw6nfPfjHofTkmtLMVLD3iRnvRkxzKpFPP8oxxVN27HPzjDj3Z0xndjjSjc58T1Oy5p25FHWmtAT5LUi8V2mpUQIwSir2Rtl5LWiG3ksJf5D1zgtdHP89j+6ecCHuyPSStGKfFYaAa3IZZCqRe/BlaWEnVFOKTWB51DUmq1BzrXlxiPDf+utkJdXTYhv8b5GLlJx4onwPPJfZ7nKbwH/JvB9QD6RZ7VYvmBsDzK6iU879tgbF9ybZDQjj04joBH6D9zvUTHiozHm2HcZZTVbw4zIdXAQ3DyY4gCvb7YeMoZHq2+yStKNfTY7Mff6Kb5rjPxcjnt+wjjqjOLA7GYDz6Eu4cZqgw93x+xNCpSGa8sJvUkJaNNzEPtcXYrJq5q7vZTlRoDvOeyMC7buDugmAWhN6LlcXor5eH+yMFQCczLZGuRkpaKoa7Z6KWutkCvHpDcOJyXfu3WIEGZuxCiruddPaYQeV5YiHOFwu5dyY6WBVIo/unnAeiskCTwj6xEH7IxyAEJPsD8pmRQVq42A9/dGCGC9FfKltSZJ4JFX8jPlifZGOYNMstrwqZXiXt+8nnbkLXSZlhKfyHcXSfVG4BAH3gMCifBwvup4fikrJZ4jiP1PGvGOngifdiUbnM05DLXW/8cTeTaL5QvKA/LXK8bAxL5DWn7SoHWW/MPRGPO0rPneRwcI1+Gnr3XZHeX4juByN3ogzPTh3hgBi110ErjcPUy5189YTgIC34ROXEdQSf3A7vGo0Vlvhdw6mFLUmkZgHpPXmp+5vkzoOdzcn/DHHx1ydSkmCVyySpGWklIqbqw2WG/FbA1S8lpRS8XuICMOXLKqYJhWFFISeg6u41AryTQ3jq+oJBstU0bajH1i32OQVnQTc2J489YhtYRrKxE39ybcPkwpq5rbPdOdfWMl4Y3LHRwBUiqkhuaqCYOlhUSj2WxH7I4ytvo57+2O+eqlFkuJTy8tUQq+vNl+qIz107A9yLiylHB3kHK7l+E6IJXpBfn5l1cX1zYJepeXZp9dXpn3792d0eJEt9aKHsjxPPi9qBlmFXcPpyjg2zeWF2s4WoxwHg1yZ3EOvyeE+E+Bvw8U8xvnctsWy4vA8Z1d7LtMCkkj/HQNSs3QTGibV9dsdmJcV3Cvn+O6mi+tNok9ExYBYzzf2xnTiTy6iWkSO0xzDiYlroD7RYXvuUil+aVX1/jS6oNJ3KMJ19h32exEbA1y4tkc41fXm0S+y63DKcO05PpyzOGkQmnwhCAJXO7sTPja5Q6FlKSVxHcdpNRMSkknCdBKc7c/JfBdEs8hDj3evjsgLRSBp3EdwSA1FUud2OdyN6KUCt91uLk/ms2XaLAzzBnnNf1pya3eBLRAo83tRcX1pQZvXOnSiT6Rk+gmUElTRTQpJCutgD+/skFWSm7uT4lCl04UMC1qVpvhmT+n4xxMCnrTgnFWsz8qaMcezdBlmteUR6Qx5o5/lJnu6qJWrDZCfEcgHKikXhQBHA3/zb8Xc5HDayuJcRK9lNg3Dvdoj8N5NMid5Uo/P/v720du08Cfe2KrsFiecY5XtbQin91RwaVOOBt7efYGJTMmsjkzbgGV1IBmb5TjCodSKuKZXtD9QUon8pAaGq4xitOiRgjN7V6GEIKvXWoRxh7vbI/45rXuA891POHquw7XVxKyUvLxwZRRUpLWik7kk1VmdGcpJUXpID3BaiPg6nKC40BZmyqd9VbEpKwJXVO62Z+dWqRU9MsaIRzSQoIAjYMQmlbsEwcOP7o35MZKQhx4fOva0izHYJzb3rjAm5V9gsNay6eoNf20ohX5vL875o0rHdZanwwoCj1n0WfQiT0OJwVQgxD4c20krRlm1af+nObMQ0p7o5wr3Yi1ZsjBpCDwTdI6Cd2FQz7q+NuxD3nFICu5thyzNyoIPI3viIfCf/PvxWsbrYXBT0vJvf6UW4cpr2+2H6jeOo8Gucc6B631Lz+xZ7NYzoknXclx3Mh2Yp9ffn19UXoJIASLJqxHPd/RkMBaK+LW4RTfEcSBy6SoF6WteSXZHRWsNk1M3RMOncRne5BT1orXNppklSSXitXQXcS2j093O9ooNk9m3++nHIxy7h5McV1BN/EZZRLXgdVGyFLT/PzVy11e3Whz+9DkAEJPcOsgJS0lUeKQVzXTQiGEplaaUVbjuRWR7yCVwgE6jRBPQOC6OGimpWJSFIsqoMCteHt7wDvbQ6paMalqfCGQGoQ2oZlRVqAF3O2l3DlM6cQeX95s0Qx9NCxktCPf4W4vw3cFlzoRpdSUUhHN3tvP0ki2CCn1M0Jf0JydXBSan7uxjNKfTAw46vgB3t0ZmWa4Qp46je+k7wWYwU6vrZvP7Xge4Twa5E69khDiX9Fa/7dCiL9x0u+PKalaLM8Mn6aS49M4kZOqWjba0QPPN9/FPapy5GhIYD76815/SuC5XFs2PQhKm/h67Jsh9EJrbh5MZobADJFxypp26OMIwZ1eSuA5HE4K0rJezC8++vr2xgWR5/L+3oh3tkaz0IfE1ZoodEjCAK01VzoJSegSeoL7/ZSXVhsLKYmqEZJVGq01h5MCqRSBB+NCs5L4LCUBtVKgwXVdIt9lMCkoFKwlkktLMQJNNwn43Xd3GaYVb93p4wozL+L2YUpVK6QjmGQVGsFyEqC0SbpKDe3QJa8l/98HB3xlo4UGDicFWSkppcb3XDwBh9OC1UZEHHhcW04eGl86KWo+3BtzdybZcX05YbPzcNOemfnsc6kb8dHehKJWNAKXjXa06H94+26fJPBm19LktVrMckALslLOBiOZfNV8ZvZp34s5jwoVPe0GuUe5mfn55GLa8yyWz8hZKznO4kTO4jw+beXI8ZCAI2CjHT/kTD7YHdNNfP7gwwNcR9AMHPZHNYNpMavsEShZsjsumBY1f/ZVo/z54d6EaSF5Zb35QLPZD2/3+fhgyt1+yiSvaEcBie9wOCn5qJey0Q755a9uMJyW/PDOgNfWm3RjRVZK1lshy42Al1YaXF0q+NE9wTiXLDUilpOQrVFGOwyIAkFvWpJEJsRTVhWjvKKsFVXlc20tYZRVvLc7oj+tuLYU0Yl9emlBHHpc6cZUM4Nba4h9wagoGU4rXtloMUxL0sKhE/uEvst7WyOS2PSK3O1NSUvJpXbE9jhjb1hwfUUSjnJ+NJvm9tpGa9F1/fbdPrujgkbggtD8+P6Qt+4MeGWtSSklH+dTfrI1ZL0dmbyGglbss+G7lFKyOy6o5JCvXTFjV/tpxY/vD813JHK5fTCllxU4ONxYSR4Kax13TqvNEKU1S0lwIVpKxznVOWit/6vZ3//B+S3HYvn8nLWS43FG/awnkE9bOXLWkMDBpODWwYRL7YhCmqauMHCoaonUgr1xRlZKskLiey4396d87apLOVEMs5LDacEra02U1ry3M+KHd/sgBPsj04k7ylIqpQgcwWY7otaag0mB7zh840qX0HdoxR774wIpNZVS9KclvWlJNwn481/fwBOCH28N+TM3ltkb5exPSiqpSDyHe1XNJK9xHYfVhnEC725PiTyPVujRCFx2xyVFpRhnFTf3JtS1QipYTgImVYWUDmUlQQjSskIpjXRAo3G1Zm9a8EY7mjX0ubx9b8A7u2MutWMudQVVrdkdpWx2Qg7GBZXUHIwL0kqy1U+JPJckiIk9j6rOmZSSd3dGXF9pEHkOd2fjRwPf5XIn5sZyg91xRl5pYt+hHXssN0yie5xXLDcC3rrTY5TXRJ6gGQZEARRSP9A4CA87p61BRjcOKENl3sNz1lI6zlma4CLg14GvAYtgptb6X3uK67JYPjNnPZ4/zqif9UTwWSpHzhISmBQ1RQ1rzQAhBEUlSfOa2wcpjVmfxM4gRSMIPcG9XorrCK50IxxhSmD3Rjm7o5y80hSlYm+Sc6efEroOgeeilGYqHDpAJ/Zp+B6eI0hCl1Fec3WpwdYg4//6yQ7r7ZCqqnFcj7SUfPNaB6lnr19qXlpt0ogKPt6bcl/muAiWGiGx7xJ4DmklyUvJD+70+fZLy1xfTnhne8SP7w2RaJLAxfMcZFkzzBXXVhr4QlBKzf6kQCkjN7LWDJkUpgLqpaUG3STgne0RCk3gCsrKJOwT32OU13jCnHJG6YRSmkTxWiukGbo4jgmdXVlKkFozySr8pgmv3R9kZhSo9inqmv60wBFwfbnBWivi9uGEI+kG+tPSzG+YlrRiD6UFw9z0N7yy1qATB4vP/IPdMeO8ovO/R/gAACAASURBVBl6i6FCIhJUUtIILk5P6ShncUn/DfAu8M8A/yHwq8A7T3NRFsvn4ayVHI8z6mc9ETytypFm6BH6gl5qduP3+hllXZPXmrwsQDi4rmcax3JJK3Fphx53exn9tEBrzXvbE64ux2RFzcG0JCslkecitWZUSJqBQ+ybmLgA3t0e0YpMh7LpYtbcOpzQTXyWYp8/HWRkpUn4frQ/5quX23QTj3e2xgzTilxKxnlNXWsqoOOZjt+0kOS1JHAhLyVKmc7htKyRKBwcikrhOYJSaSqlUUrTTHwGeUUn9pmWFY4jkKpiWrjsDwve2GzxzvaQslbUynSWe55DM/IYZzV1LRHCGO7NdsgoL5HK5Eyi1SaBEASuS29SIBXsjDMqpbjXm+K5pgppqRHgez7rrZhm+EkPg+sIU7c5Y1pJemlFI/BYb8aMi5pxXjEoykU3+5y0rJFKE4WfaDAFrsOkrhcNdRfNWZzDq1rrvyKE+Ita698WQvxd4P982guzPJ+cxzyEs4ZtHmfUz3oieBKVI7ujnB/dGzBITcjmG1eNTtKNlSZ/8OG+KRVVirLWVLUkCc1QHE869NKStDSVRpOy4k5/SjJLdmdlybSUbPWNBlFeK5LQoyglrtAoLfCEw7Qq6TQCsrqmThUf7I1ZbYW8vzNme5ATeg7jQjLKaw7HBWVVk4Qe72yNUBpC3+VOPyXNaxxXsNYKuddT9KVCSkXouQhHgBDEnkMhFVllXk/kuzgIxkWF1kaQMPIEh2PT9ZxWkmbgEfshroB7g4qlBK6vJPzJ/RHXlxt0E5+8cjkYlyzFPnmlcFzBKK8ZFjW+J+hNC0a5pBsLXNehqhVTDbEn6E9q0koicMgrszYXhdYaz3O4uhSDVgwztcgdtCIfgSCvTPOfL4T5XkceO+PMSJu4gsh1+WB3zLWlZPF5J4GH64iF7AdAKY1zPK9hPo/jLKuoZn8PhBBfB3aAG09tRZbnls+rB3OaY3mUoNmjeJxR/zQngs9TObI7yvm9d/fMrrmSvL875if3h/zKGxvkleSnrnQoasX3Pj6kqI38RVaaiqC8qnGFoB16ZLXkvZ0JtdSsN3wO0gohzGyCg0mB7zosNXzGeYXrOiSeYDIrhe02QlyMdMVSEhJ5DnujnI8OpoSuQyvy+O5HB9zpmfe51ppu4DHMKrJS0YlNk146GxLUn5S4rjG6ANOiJvIAfG6sxOyNM2LPY5gVKAkKTei6RK6DcM1JJq9rAt9hsx1xt5eB0Gx0IjqJoJME9KYFdw5SosChlCFpKckqyTAzBnac1/TSiryqiX2X/dJ0L/uui1KK5SSg0zBaS3d7GdeWY66vxAympsEwqySjQcHeOCcran7mpSWEEIvvyrz6af79WWmFfONql4/2JwykMfrNJAA0kecePWRwqRuzN87ZHRVobXIOk1xyqRNf+FzyOWdxDr8phFgC/l3gHwBN4N97qquyPJd8Hj2YR80LOC7//GkczmlGfe5w8kpyMCloht4DKqRPkh/dG+A5Jj4duGY4zSCt+IMPD8zatMsgrWhEPo3QYzgtkUqhlMZxHEDTijwC1yMOHXqTksO0pihrmnHA4ahAKo3jKJQGIRxi32GpGRJMS7QQtCOPWmp0bSS5398dEvseK42AG6sNDsYFh2lFKSUCjZSKXEqmUyMeeDCt8F0j3zFIM6QWdGMf4dTUClwgDIz89OGkQClNI5Q0Qp9JkXE4KQl9h1QbgSWljCOoSomKNGImnTHPM2g0Va1pxB63D1J2gpJm4NGJfYpK8cHelEvdmMvdiLt9k4tpRj7dhk9ztmuvtTJS3lnNRitksxPxwZ5xhtNS8oPbfaLA4421hEkleevOgL/87Wu8vNZ84Dsy35R842rMzb0JO6OMr7faRta7qFlvRfzsS8sPzERuhsa5HK1WenW9uShBfhY4SxPc35n98x8DVonV8pn5PHowpzmW7358iO+Yhqk4MPpBx8XZPm3oaneUL2QM2pHHUhLgOuKJhcCOG5WdYY6DJnAdgllcOvIc3tsdczApWGkE+K7DSysJB+OCpWaISM1cgY5rErYaeGk1IS9rRmlJP63whKCuFRXKdCsrgSME37jSZn9SstVLaUYunitIS8XeMENpo9jqOFDUNfsTzfYoJ/HNznec12ilWWmFOI6gNykXcx5MHF2RVZpaasq6IHAFjcCs1UyQEAyzkrJUDHOXV1ZdQs9dqNFKBVHgsBQHjLKKSirWOxG/8PIyd3optVJEgcNhWpEWFZHnsjXJ6QoH7SmGuUZruNKNSEIXB3CFxncFDpAVksT3aAVG4rufllzpRtw+nPL+7oRBWlFUkq1BynonpOF7VAi6jYBr3YSb+xNeXmsyKWp+dHfAKCuplcZzBO044JX1JofTgvv9lKVGwFevdHhp2Uhl+K5zYm/FL7++8cw4hKOcpVrppCa4IfB9rfVbT35JlueVz6MHc5JjkcoImn3jaoeGaxRH390eoTX4nrMQZ/u0oat//N4ew7TEcQTDtKQZVbw0G2hzdCTlwaRYTOtabYafeTxnWtZM85rLS/HiPXn7/oCtXkYn8Tgc+Sihub6cMC4qrizFXOpG7A8zppXkajMkcD1eXmvwnY975LUmLSp81yH0BFmhaAQencRnWkj2p7P5DEKwmvjsTyruDaeUUpJVRtjOwWjze64k8Bx6VUUtBa6jqRVGvVWA0IpMO7gODNIKrT6Rbq4keEKjEfiei+MKskpSF5L1VoTjuNzqZYDCd1wqKYl8B6EF/WmJ5wrWWxGh57LcDAkDl4/2Jnx8MMUVgiT02BqabvH+NOdwYrqXl5KA270poWcSys3IqMceTEqascvrm21WWgFvXO6QlZKDcU4nDhhlKZO8ZJAZ8bsrkUfgOSwnAZudiElecfvehBsrDQ4mBdvDjFbk0fTNd297mNEIXX759Y2HmiKLWrHSDHn7bp/bh1PK2kik3+2l7I8LfuGV1WfOQZxlNd+e/fnfZj//BeB7wF8XQvxPWuv/5GktzvJ88Xmqeo47lrSs+d6tHlkp2RnmbLYTM31sOuFwUrLRjrjTS1lrRYSec2Yp45t7Y+4PUlYaAYHrUinF/rgg8h18VxyRoFAcjHMQDnlVLHaFp3Vh39wbc6eXsj3MWWuGvLbRXIjHvXGpzT/6yS7JTJ7iB3cG/GRryEY3NqJtrsuf7owY5zUrjZB2w2eQVRRSMWtGppOYwTS9cU4jcKmkTzUTylNaMSnMIJnNTsSVpZit/hTfdVntxBykFXmpyGpFPrPsHqCAUmGSswIQGuGY300KidIQOOB5irTSyCNB9fn407SGxFFobbSCHDNplGFe0YmgrGrySqG0pqg1TmVE9CLPwXE8fM9hd5AbLatxjlCa/XFB6AmkNg13UmvyuqaqFIEbLUaZbmcFrdgFBVHostGNuLqU0I59kyCffTYIh3bkEYcunTigqmGrlqavxHMpasnhxPRIbDRDKqn43sc9XlqNCT139n120ZHmTi/lm9eWuNyNHyowGGVmSNIolzQCl07sMS1Mjmm1FfKtYx3cF43z+LuwAvyM1vo3tNa/gXEUa8A/BfzVp7g2yzPC3CC+Pesy/ayywPP4fiUVP9kecXN/ghCPfxwYx1LUirySTIuK93bGDLOar19uk5WKjw+n9Kcl93oZWSVnw2GMAqZU+szlgXd6qYmVY6qBAtclCU0zVBJ4i/DWKK8JfZd25BF6LuNZPft86Mucefjhw70JgSdQSrE7ynh3e7xY0+VuzM+/vEIlFe/ujNkfF6x3QtqhR1pp0rrm2rLpHC5rzXc+OuCD3RGtKOCb15ZwtODH94bc6k0BI70Reg6x55DXijw3ZaRbg5T9UUF/WlLWyugyzZroPAH6E/VxaoxzmBuIUhtHMSmMwa+0cQAaUFJTa3NikLPH6Nm/FTApzSlCSTOUSErIcpM8npaKUipqpZlFnow0t1IEnqBWir1JwcHEiPKFgcvBOOdeP6M3zcnKmmkuGWYS4QgCz2UwrTicFEyqikkuTTNf6DFOK+73p9w6nBL5DmkpF5P4itpUIXUSo8/0szdW2GhFxIHLx/sTpkXNpKxxHJc/unnA+7sj/uTOkLw6MuJGm8qlt+72+f339pBS8fpmiyuzvNjBpOBwnJP4JnwoEDQCF0eIRZjpWeIsJ4frQHnk5wp4SWudCSGKUx5jeU54GhOnlIZX1ppn0iGac7Sy6Ob+hEbo8pWNFp4raITzqWEj4sBlpRkuVE3zSvLdW71T8w/H4/95JVluhItJZ74rKCtNWWsudWNu7o1phh5ZaXZ/AIEnSEtTzrg3No+bX+9gUvDuzpBKaoraJ/QclDbKqnvjghsr5j398kaLK0sxO8OMf/j2Flpr2lHAtKzJa0h8MRO5q2bXAKXNBLerywmjvCKtJFIrPMdBKcXh1MTDtQO+AwjNtKrYHeWUlURqzV7oUGlNoTRHRlMsmBv5o9SzE4KafZbzhPNJk8Dc2TUcAdKBcHZDqWCc1bgupCZChYJFU5kCBmlJWkiWkoBpIekkPlqb70JvWpBVNUIbZ1MpmBaK+/0pruMQ+A5KarKypoxcJnnNpJRMiopLnZj7g5xpIfnyZotLnZhxXqG1CY1tdszM60bo8f7OmFKa2deb7ZhRUeFiTjZvbw8ppOZnX+qShD6H0xLXEewMM7PBEPDezoQkdI2ibVoyLSTBkQE+lfyklPVZ4yz/u/8u8MdCiN+Z/fzPA/+9EKIB/OSprczyTPCkJ06ddr0P98Y0Au+R/Q9zBzHPP2SV5NbBlMAzMtSDrOJyN0YgKGqJUrA1SMkqyc/dWDpxstpxx6c0TEvJWiviYFJwp28kKuZS2PPw1nyqWui5lLUm9l36acX+uKAT+wutnX/ywT6hJ1hrRtRKkVcmxFK7grSo6U0LtgY5a62Qu/2URuCy1orYH+dMipqDifm7rrUJjXRCpAxxHGjHPtuDjO98dMi0qGlGHqvtiKxQ7E9KNJpGaMaDCsDFYZrXuE5JJ/IYTSr+8MMergtF/UA/14K5wXcwBnuOuR6gZ0b9hMeAcQa1xgzHqSBX4HvQiV2yUlLU5lpyfj0FatZb5qHxPFOuuzfMGKQFw6LmYJqTl+a6oQPObHG1NgKAvitxhKbWAtcxA4GG04oa05sgpSYva+rArCHyjYjeajPk/35nFw0sNQNc4fDVKx3euNLmzqz73JuFspYbIaXUTIuKt+4O+Nb1Lq5jNj1745wkcBikFe/sjBAIbqxE5FXNzijnXi9lvR3SSQKU1rQin+vLCc8aZ6lW+o+EEP878Gcxn+Nf11q/Ofv1rz7NxVkunic9cerkxLIZaPNTVzpnEsDbGxdUUrHcCBcyyIO0pBv7fGm1SeSbHfzHvQmO4/DlzcaJozxPclSvrDW5uT8mLWruHExBwHo7YrMdPVA6245MiKmoTXBluRFzv5+y3AjZHeVkleRwWhKHHmlmyjkDz6UZmRGQ00ryRzcP6Kcl7chnvR1yt5eh0Kw1Q/JK0s8qJlmFKwSuA75jwlu+axrctgb5osa+FRkpjP1RiVKKsq7IK43vaJqxqfmfh7GysqaWpknMEZDN4gIO4AljZI8fIo47Do25z/FvwdyJzPfCuTT3nZ9KHMCbzYaQM6eBa24PBOS1uXgcCNAapQTS0UyrmmntsNNPzfPO1pgp84SCT045joRJbmS6kyAgLWuyWtGMfHzXpVaK24cpH+5N+M7Hh/zTX7uE1JqikiShSyfyZ0UNgsvdBnlVc3NvQi4ke6Oc+72UQmmWE58kMB3UjdCnHZlKrO1BRq0Uu8OS0DMVYuNcsj0qeHWtwd2ekTyfFBlfvdTi+oopYX3WOFNcQGv9fcwMacsLxuedOHU8bDN//NHrzQfanFUAr5KaD/cmvLoOS0nAZjtiKQkWhtsRgpeWjfqnK+Cl5YS0lLO5vqZJqx37J45uXEp81loR7++MWGlFLCcerdhnmFXEgbsYyLI9yKikXlQrdWKfUeYzyEpCzxjvu70aR2uUUEzLGo2L1prdUUElJRq41InppxXfv9UncB18z+F2keIIQT47DVzqJpSVZJBW7I8KKqWRSqOYVwNp6lqz3Ay410+ZljVKmcSu4wpi36UqNQOl8AUgzHtcSaMyWk4rPNO8bPII9YPOweVhZ8Ept+lj/z7uVOY5iMjVpudidsfYB4RA1qYsNys1sQ+Rb4oApmWFUsaZaPVJwlsfea45EnOfAE3iO+TCpVaQuILQc/nx1ojQc1iKfcaF5Hd+cJ+fubHEL726SlFLPtybst6KWEp8ilrhOg5/5ksrfPfjQ36yNaSbBFxumiFN/bTiW9e7NEOPvVGO4zhc6Sa8eavH3iTnSifCcwXDrOJyOyIKXL5xbYmlRsAwq2iEHt+41n3mKpXgjM7B8uLyeSqMTgrbTIoagaCbfHK9YVbzxqUHd04nCeApbRLGWWWSjNtDo1dztLO5EXqL3XQj9FhKfEBw63BK4Dp4rqA/KfkfvnuHspJoDXHkEnkTvnVtiUZoppJtdmK6iZFHADPE5Z2tIVHwSff08bDarcMpIBZCaq3IZ5RVrLUiLnVitocZVa0QDsSOy8e9lHFqBuPktanY8WbxZ60dLnUjQt/jy+tN3tsZm0ltUiHQTIuavJIzxxgyKSpuHk7wBLQij9j3GOUVAtjtZ+iZJfUCh7VmTFnXRoZaQ+CxqHqqFXgzbzDzI4vk82kcNdLzx5zkOOZoIJvFnpz5rr8CR+jFdVwHpIb9cUnkCjzHYVTOqqacTxLlx/McsWvyFo5jJM33JgVoTeA5CEcwyatFh3fkOSzFHn7g8pP7Q3xPsNmO0Vrxg9uHLDXCRaVRI/T43q0ea62IRuDhuAKtNN3EZXdUEMw6oG8fTgk8TFlwWfPu7phffHkFgE7kMyokl7tG+lzrTzYXc85DXuasWOdgeSSfRzfopLDNUhIs5gfPr/f6ZgvXebBwrj9r5Hr7bp87vZSNVsTepCBwHZLAxRWCQVY91FF6tON57pzu9c2kNdCMs5qDSYkAsqoGBGWqcBJ489Yhr220GaRGXXOYlWy2E0CzO8rIKs3XlhIqqfjR3QFJaMohj/4nzsucopYErqmQORjndJOIr2y0uLFi+i7+8MMD/nRrgKeNeJvSMEwr01sgTCdvXtZoAeO8ZFr0CV2XbuLTjDz6k5JWJHAch2ZsylM/2BkjK4XjO2iMXIVAMy5MSCv2Hbotj3FZM05LEJrYc3Ec0xOSFRW1BI0yFUWY8E/kCQKtmVSnO4ijtz/KKZzE/P7yyAMdTAJbY6qblNK0IheNQgvjyLyZ85g/+Sw6RRCY3IbUUFYSzxHEvsN6O6aXFriOQy1rAs8jrRXLTY/EE9waTvnR3QFbzZS9cUkSuvyL15ZwHcF3PjpkrRVSSckrqw0+PJhSpopW6BEHAfujnGbosdGOWGn4vHmrz+44ZyUxYSwQBK5gmNe4rliMOT1+Aj+psW5/XFzYycI6B8tj+ay6QaflKyqpThy6M//91iDjh3cGrLeM0FotNd+/3ePKUkIYGidS1oppUfF77+7y+mabdmx26Xd7KXd7UxzHoZuYipMPdidEvsPlpcRo63iOMb5a4HuC24dTbu1r1lohjdBjvRUSug5bw5z3dkZkVc0kl8ShSzPyUVovGp5eW2/RTyt+sjVcCKc5Shvp5cjjp68vkdeKvfEnDXMfHUwYTksqDftj0yNRVjWjvMTzTA7GdwWDaUVz9rqgMnLcSwmd2RruHqZM8xo6GoGm1hpfQ+K7+J7LVn+KVNCcrbtWCs9xGRf1LKHtcKkTc7c3RUq1yD3MjX2tYFRqXIzxPaka6Wmg4ZN6WMyo0GlZ0wyNUJ1U5jQTzBIbvgflLKEutEsUQFlLOo2I9UZApTTjvGatGTEuakrHJQpcmoGH1pqdUUlaKraHuSlKkOb08Qcf7tONQ3ZG5sQ3LSXN0OXV9Sa1NHMnpkXFIKvoTkvySrIzLNjohKy1InZHGUnkIoRmPGsg+ZnLy8S+Q17Jh07gN/fGpzbWHZ9idx48akzomJM3CwLQWuv2U1uV5cJ4ksfaT6NqOm8a2hnm3DqY8KXVJpvdyBgDrRllFbFf0I5McvXjgykvrcY4AoZZxfdv91luBLy7PTLCdEXBO1tDpDJzAZaSwMg81Ipm6HGvlzItJa7jsNYKSQtJIRXbg4zXL7W4c5hSScmHu2MGWc0rqw1eW2uZuLKAZuQiZw1iO7NYc+IKiloznJa0E59hWuI6Dl+/0uGj/Sl3eimjtOTj/SmjtMRzXRq+Qz8tmBQSFwfXMbOOpdKEroMLeMLMLthcikgL0707zGpCX9BOTNNWoWC54SOEQymhlhKtzY418FwTkqqVmSDnwEozouE77AyM9PdaM+ROvwABrjYlqvP//OflFI4iZr0aQitKCa6CUmtToopp1ss1LEWCTuIzSmsqrZBKstqMqRWstAKuLcXszMpWBZpr3YTdcY7WUFQ1pe+yOzK9KautELQ5tXiOKUP13ZRryzGu79JJAv7ww31eWWvxjasdhllFLy3Z7MS0Qpf7w5ztYcZS4lMrqKVitRHhCsHL601+7ksri5njJ53A78zkTE5rrDtvHjUJ7tlLn1ueKk+6p+G0fMVKM+SD3fHCAbVjn62B+U+1O0jJa8ndXkpzpmvUjnxWmz7ZrGt1mJXcWE0AQS+t2Bnm5LXkB7d7NEOfOHC410+JfJflRsQwqygqReQ7DDMzG+AwLQkch0YgkMrEpH3PwRGae72U7VHBB7sTskqiMclf3xMEnsOdXsq1pYTYd9gf5zNNJMG0lFxbivnerUOGecVXNlq044A/eH+fe8OUThyyO8oYZSWDoqIdCQLfpRGFuJ6iEzrUGtYaAQeTglpKehPF9dWYSmnSvGaiK1aaIVUtOZgYjaIvX27z6lqDW/tTBlmFN9e+0NoYOlfgaIxEhdYmXq8USRgShR6TzCR8TdOZfqDf4LyYHxTmfyppGuG0frDaSQtFI3TxSklaQ1r9/+y9WYwlWX7e9zvnxB53zT1r33q6umfp5rA5wyFtSz00YUK0SNuQYdovhmyDLza8vdiCXgwbAmwYfjC8wCAEw34QLOtFsgxJpERzTNHCcDjD4Wyc6a26q2vN/ebdYo84fjhxb2VVZWVlVmVWVXbfD2h05RZx7okTZ/n/v//3aZy0ohXYVNqcKLWG+cCi69us9dP6GZqdf8fXfPVCh/fXRoxS4+0Q2Ibum5YVSVbgWYpxrrnTi7m22ERogVd7Rl9ZCNkaJ+xEAWle8NULHfKyYpyWtHyHWzsRm6OMpZZLN/QJXUXoWFxfabHc8lhumZDSfhIs9/sxy02vVq9l8oE5OONzcjj0Gy+EWOJhJ7hbJ9KiGV4ajrumYb98xXzDfUxF9Tsfb9MNHPpxRpRXLDZcxmnJe/cGvH2xi2dLQs9hwVZcXjCJvLIyO6qL8z5345zeKON2L+Kr5zsMooxRUmArie9ItsY5TddinBb4jql8ng8deqOUUQaOElyYC+knBULAj+70WWq5aDRSGHrMJ9tjNkcJP3+hwygt2R6ndAKHm1tjOr5N03MIPcUozbm62KCoNJcWjHrnzZ0xldbsjFI2R2Zx6sc5wzhjpdWk45uFZaXlYFvKmOFoI1HRcBS3dkz4aJyVtH0bS0qEkpybC9G6ohflpHmBlAJUHZMXkm5oMc5LRolxRpMYHwNbCXaTnHFt/NNPcoqqQghNUT5OT30ReCxXUSeVSwy9Nq8T10kBgpK8NDkGWxi57zSv8CyFb5vNxvmFgE83I2xL4DkSXdkICZ3QwbEVb19oM4wLsyHaGDHKC6zYmARtjTI8S1KWFaWuWBumXFkIWR8kFJXJ1fyF1xf5dNvUpQySgigzhYUrLZc7uwmjpKDj2YzSkqLSD8lwP0mCpeXa3Ngac21R0PItskIzzkquLL5iHtITCCF+A/jvgDPABnAR4wT3xZNt2gwnif3CR8dd0wCP5ysmJ5O9C1BZVXy6NWSp7dP0bKQQpPVObnuY4NrGN2Cp6XF31xSlDdKCi/M+bd/h482RmfCDB1WqrgVJXtW6PdAJDHWwG9gErsUgLsg9UyvR9C3avsNi0+O9+31AMEyM5k9ZaVzbQghdy3QYn99hUtL0jGrrICnYGuf84pU5tkYpShqF2CgruLU95r37AyO3UedbxrkRmMvyysT/BVxdaBgNIym51x8xSDKKAsYYldE4Kyjr7XyUF+gKVhdd0rKi4Vms7Zqw0eVOSFoaMb9BklMUJVGq0TVVNbA1aW52o0le1L7MIKRJ/r4KHmQW4Dt1e2oWlUnrmsUiKky+wbbAdSwsIUiz0tS+BC4NT7E5yEiLEikVXc+BquTMnKE8DxPj7xDlFYMk59J8wKfbEbtRRpZX+I6prvZsi97IbCw+3hqy0vLxbQvPgffXhni2MUDKypL5hsNHG0Pu7sR0GjZn2z5xXjFKE7525WFRvckmbG2Q4dqq1m+qUIGgQrM1SlDSQ0nBcsvl2kuqgTjMyeG/An4R+H2t9c8JId4F/s2TbdYMJ4n9wkc/vN0zImaVpu3bU97/UWoaDoP9FqCWb3NrJ+LcXEjgWNzajowrWKX5ZNsk6L56ocOZWl8pdCyijSGuLdFao6QgLyuuLrb4aGuIZ1vYyrR9mBS0XYub2yPW+ilvnWsz33R5Y7XJT+5WnOv6zDcd0Gb33m14jOOUn9ztEziKSoBrCaoK5pqOqTr2FElW0I9SpDBy1BL40e0eQhpNo2vLDd67P+CT7TFxVrI7Tuu6BElWTuL/grKssCzJ3d0YW4KuLHpRQZaVCCGpkCR5gaUknhT4tsXWMMVzFd+9uYPGnIJ8W9IObLoNh1s7Y3pxTlmUJJmhh07oqSZ8YnbqrtGkQ0koagmKVwG2hVmtRIWFWbBk/f9JkMW2TFGhrUweCwm2FLiORGuJpqIT2Piuw0rb4Ypj9BU/LwAAIABJREFUoTVsj1OSomJUn5wcS+Laik5o49mKKDee05Zr8fVrC3yyOSbKKpIiM9dzLC4uBmwOE/JK0/JsHCWZC23mAhe0IPQsWoHDOd/GUYJK64fCqFujlKWm+5gES5HBl892uLk95lw3OBVU1lxrvS2EkEIIqbX+lhDivznphgkhfg347zGbhr+ptf6vT/qenxc8Gj6qNKwPUgSGPjdKS0bpiNW2j5Limb2Q9zud7JekdpRimBR8/9MeaVkR2IoKo0MUZSVfvdDhbG2x6NmKTgBzoUOUVtwZx4zTkvmGi2Mpri41SPOKjaHmTMPm2mKDn64NGCQ5b55pETg293oxi02XN8+0sZRhwIBmte2RFyWZb/PB+pi81AzTnMrTWChCR5JWJo4fuoqllscgzukENnlZEWUlZzrme3d7MdvjlHu9mOW2xzgt6MUZu2mGoxSObfIUu1FuJpaGS+hYplZCV4SuTT8piMZmAZK2wrYs4rRgmOYMkxwhhfFWTnI2hxUtz0bWcmeBrUioyEoTtC8rswDUBcUGlRHUK0sTnhEYOYr0EIvEZNJ+Eh6V2zgsBFAUYKkK15aMysrIakgTLqvLHCi1SRprBJYCx1a0XLOpsZWkH5emjiGwsKXEUpKtYcKHGyNeW2zUyq/QizNzQiwrmq5CSaN3dKHrsRC6dAKbD9dG+IXCUYpzcwGLDY8538iCbw5SQlfhWoKmp2h6Pq5t4duShYZLP87YGRsP8MlGbMJQ20+CRUnB9ZXWM7sKHicOszjsCiEawD8F/pYQYoMTPn0KIRTwPwG/CtwBviuE+Pta65mW0zHg0d37Zi3zXFSai/MBG8OUfpTRizK+8Yw6809zbgMTsupFOXd3Y752eY7vfrJNkpvQydXFBoFjmDtZPfFuDo0shWcZn+F+nFJpTcu3WB8YWe1fub5iagGijIZrcWsn4ucvzvHJ1piWb2FLQVYaQ59fvDLPKC1ZaLhTWY650OVna0POdDzWhwlFrtlIE750psNuVHJpKTQLmzZsklGagNac7YbYSnBxPmRnnPH9Wz3u7yZYUnBhPqAfG7vKKAUlNaEt6YSGJeVait4oZygLJNByjDx1VtQ8IaFJC5PwHGclEkGJxrclgW0xzjM0FaEv6UUZeVERFyVRWjG9RE37rDAv/YQuqmrhvFybXZhvC8j1UxeIp00Az7o4aIyyZ5mCb1XTSuiyFvcDcC1z2nEsU3G+1ApQShDaCkcana331h8UCUZZQZxV5KXmtcUGlxcb3Osn9KIET5nTgtaCKC9peQ5xXhDaFnd7MW+eaWGtGkab2fF7aA2f9iI8S5HbFYM4p6w0viNperZxw4tzVjs+lpQstx6u/j/T8bjbizjb9R+TYDlsgemLwGHe+t8EEuA/wWgptYH/8iQbBXwN+Ehr/TGAEOJv1+2YLQ7HgEd373FeTuPkgWNxad5CzwWPVW8eBU9Kbu+VnxilBb0o49pSyFzomjh/VjDKKrKi5M0zLTYGxmc3zqtpAdwwKbm5NaLhKmwpyasKWxihsz/+ZJtfuDT3ULuvLTVwbcnmICEqS1P1W1T85K45Ufz8xS6uJflwfcjGIGG+4XFuPqzpiJool5QaXFugK5MX8G011dGJ84Ki0ry23OTT7XGtYaQJHUVWlNzcjnCVZLEVoJQiSgvSQtOLMiwhSHIYJAVJUWIJkzCWFszbDrtJZpg8AsqiRGgjpldqTdO1ySqNLSSthounFENdMh863OvHU40JXVNTJ2WGk/BMVYsg2fX3pIASgW1psuxwHJnJhD3RWprgeSNUAiMGaMkHEhuWhIZnYStJy3PwXUmRmxNGqUFZkpWOR1FpVpo+bd8myiuklCw1XZqeouNbaCSukjRcxU6UklcVAmES3BZcWWwjJbiWRVJWvHW+y+1eVBdvwq1t07fd0OF+PyEvNasdizivuLUd49lm3Kd5yTDJubb0MA21GxjpjbbvPCbB8jLDSI/iMMJ74z1f/u8n2Ja9OAvc3vP1HeDre39BCPHbwG8DXLhw4QU167OBRymmUsAoKXl95YEy5GFzDU+qi3iSwN6NzcH0dydiY5Pf64YOhWdxtg4nGZqrwwfrQxqeid9mhZGPSLKCOCto+objvtLxOK981ocJ3/1kh8Wmw1LLoywrPtwYcb4bUGmIai+IsjKT88X5Bje3zBAPHMsouuYFgeNwcSGgE1hkpcZWirZnYVuCbuiw1PD4we0eWVnScGzavs2P7vRZbrmmFmGc8/76EFtJeqOMrKzYGCSgK7qhi60EO1FGpoWh6KZGoa6UBVkpqMoKbQsarpHwaPmK3TjHlqrepSpcWxJnJY5t0QkdPFtiWYreOMN1bMpKG2e38nGRPL3nH1KBXYIS1P7QD4renhQ+miw01YPLPITnIV+GNaW2qCB0jcuclGpKGdYCQtfiwryHZxvK6jDJsSxFUZiw5MpSgzdX27y31qcwxtlcnA9wlGJjkFBUFW+stvjxvT5K5DRci5WOh28r/rnXFtkYpNiWoCjNYrzcMs9sN8opdMXFhYDeOKcbOOzGOSAM9VbATpTxxdUWCHBtYxYU7nkX0qJioeHy2nLzlQgfPQmHYSvtLYZzMBuN8QkXwe0ncP7QeNNa/w7wOwDvvPPOyyECn1I8SjFdaftEaYkUAq2N8clulBO4ih/e7j0xMXZQXcTjzm0lH26MCOyHf1eKB0J8S02Xm1tj0kITOnIqZ21Jye44IStsuqFDJ7D57iclSoGbGwmInXFGWCuVXl9pUdXidEVlQjLb44Slhss/urlDVpScm/O5vR3x7Y83idIS31GstH02hwm7iaIdGI68Vwu0+a5iLnC404v5aG3IfTciLSrmGh5d3yLKSpSAH9/tGyMXrenUvHchzMLo2ZJxZgrRkrxCVybOvD6IsSxDsTUJbqNIV5UVzYZP6FkMYmP7eabjI4VmZ5STVyVpWdJPcjzbBItCS3EzyXCkwLYkaa1LoSoecmoTPDhNlKUJ0ygFvm1RqAqykryY5oUf0k6aLBwC4wSntcldTH4Oh18cbB6cOirMBIMwDC3fBiUEWkqWmh7bo5RhXrLQcFDK1KP8ypsrfPVil61hxv1+TFYa/aWvXpxjLnTIKxNOAk1eaXSt9rcxiPlwfYjvKP7i64tGTqN+hqO05NJCyJ1aPfXubkzDtTg/F0zHikCYBUkJ3lhtsRunDOKClZZHayHk7QvmtLAzzri3GxM41pG1yV42DnNyeGhpE0L8K5iwz0niDnB+z9fngHsnfM/PFR6lmE5OABPaqq4N781uf8j3b/W4vtJ8SMvooLqIR08nd3pjtNacmwsQwjA41vsxg9QUHZ3peFOq5+0dUxC3Ocq4utig4SozqaJZarpsDFMWmy67cWYkDRxFUlTc2DQ1B23fHPFdS9HyTSjro/URW6NtdkYpV5aapLkpXIpT8/9+nCOFSRImZcWNjSHDuOD11ZbxRx5n3O8nKGBzZJLEb55tstjwyGp5jPv9mNu9mCsLIeOsQAmTOA9d076tYcZyLS53q2dOWxqNbUt8pYhLo1QaOLUQH5rAkZSFmfAsKRjExkN6PnC5208YJXEdwiiIcxhbko7nkBQVji3pSEFWViRZCTm4gqk+UYGZjCvMAlDULmxKSCxZohxB4Ch642I6+e+d9B0Jnl0vQMWDn0sOXhwmpxEL8GyzSCW1r0MBiBKEAteRNFwbJc2S03Kt2hvC6FE1PYe80HywNqQbulxebLDYcHhvbVgLLsJi0+PmttHW0nWx44frQ87PN7D7MQ1X8We3dllu+SRlyZsrLaK0QLY8Wr5D07fpBs50Yh9T0K4XfVsJ5huOSSQrj6Zr5Foc64FOWLcmK+zVEnuSNtmrJLoHz6CtpLX+e0KI//wkGrMH3wVeE0JcBu4CvwX8Wyd8z881Gq4ZjPd3jaOarYxF5npdAdypPYorzbRi+qC6iMni89HGkBubIz7dGnN5MQQEUVZwc2uMbRmbxOWWz43NEWlestRyeff6End7EVFuThSe7XNze4zWgvVBwiA2bl2d0ObeTsIwK3CVYaScnw+Isgq/dtuqKrixOabpWqy0jaLmp9vjOtabszMyAnutwDEcd0vyxdU2O+MUWyk8x6pjx4VhHEUZZzo+/Tjn1nZCYNs0XBtXZfTGGcsth9BRbAwT7vUSs8jlJUuhg0AT5SWuZfGFxQa3tsdEuSlEQ0qalqjtMSVCQmhZZuKWMNdwSAo9le5WSnJuzmchdNgcpxSVJrAstqOUNK/MqUpJPNchL0ruZiUW5vRgq9oZrjAWjwITa7cU7MY5Uhj5ioYnaXi2uV5hKLBSmEUhKU3fVphjhUUtlc3jC8N+yWkLaPpG58pWksE4JavptKpedDzLQgnBmY5nHNSUWfiXW77xXygKPtzo88WzXc51DfPng/WRWSzr02jgKC7Nh2Zzgtm4vHNpnsWmy8/WBny0PjJChVnOW+fm2BwlSCmM2q9rQlmPCkfaSvLlcx1ubAx5b22Ia1VcnAu42zOSLBNmHTwcQjoIJ+G4+Lw4TFjpX9vzpcR4SJ9oGEdrXQgh/gPg9zA5r/9Va/3nJ3nPzzv2Dk4juyz44e0eSy0P15VoLRiX5dQneb/QETycqxinBZ9umZey41vEWcXN7bExdqk1HgJHMBc6bI0s8Cy+sGyilaWG0FFsDhMuzodcmg/ZGMRsj41hezdwanmMMR9vjnEto7HT9hzu9Y2zmkZzuzcmyUquLITcqXdl26OU3XFOqStUfToaRTmbQrDaaSOEoOnZvL7S5PxcwP/xnU/ZGqa0fAcwYnnDxAiurbRcxmlOb5xTVBrHstgeZfh14dz2OMWxJIlvTetGOoGFpyTd0KWKEnSlyUuN45qCKMeykEDbt7AsRZzmZBUoUeIqiRCmb0ZJaXwiKmi4NpY0RXljSqSCpKgYpGan2/IqMqskLkpTNFZosygAnmtorkpItKxQ0kiKSGkmRs9RhL5FlpekRYlGYJeaAkM7nWgRFfUKoDAT/FTyggc5CkuYn4Wehasky02XSgjank2U5YyTAtdRKASeY5lCsySn0gLXEni28evOSs32KMW1a60QgRm1QrDQcE3tAw9yasstI7H+d79/hzQv+HgrZ2uQUpQVjqWIUpPYXmn7XFkMeW25yQ9v96by6xPs3fy8db7L1aXmdLe/2vEZ1UWNk/DsYUNIx61OcBw4zJL0l/f8uwBuYphDJwqt9T8E/uFJ32cGg72DM3AtozqpNaPEyD1kpbHG3FsxfZDXw/og4e/92R2yQjMXWji2xb1+wlxgc3tnTDd0sKScCoqVlX5Iz8d3FHlpTOgBgjoncH4unN5XCjjfDciKkls7Md3AodTG1H2jH/HJ1oidUcZK26Ufm0rXLC9NXUU6ptLQ9h2Wmh5JUTLKSjZHGa8tN1luuWit2RikSCE52/VJC2My1AocPMssoP06F5CVFf/8awt8/9YuCGjWYoJ3dmLmAps8LwkcZZgxWiCEpOXbLDddNkcpSk2WZCOh3Q1seuOMcZSxMcwoq5JuXewW5SaBPbHHVJZgZ5SB1pQY5o1E8aUzDYZpweWFBn/03ibj1FBl00yTYyZs1wIlFIUuiWsLU1tqhARLKjq+TVoUbI9y0BrPsRAaRFWg87qCWT+cJCwAqR/4R08tQAHXESgh6PgWaaHZTQqWWi6uZZLrS22f850AISHKKuK8oCw186FFkhvJlIkpkpmILW7tjGl6Fqttj9eWGlRaTyftvaEcMDRuk89QbI0zQFNWpiBxY5DwlXOd6Wc5jHDkQeHZo8jbn4Q6wfPiMIvD39Ra/7O93xBC/DJGSmOGF4iTjEnuHZyTxLCjJMO0IC1MRe+kQnnycjzJ6wHgj29skZcVC6FDUcEwz+n4Fre2RwzSAqkEgWXxg9s93q518+Os4tPt8XRB2B3nzDXsx3Zhe0NW768NaXsW776+iJKS+7sJW+OUhu9wphswSAp+eGuXpZbLXOjiWJL7gwRLKZQwlam2lGRVhS8hLysuzodIKbjbi7GU0bnZHGVsDVI8WzCIUnY0OEJwYS5gsWVxYS7g8mLI+jBFa02lNYNE8qWzHa4smtqRzVEGUuJZiusrTcZZyc4ooaQiK4xirGtJ3j7fYRBn3N5JkEITOoK8tNgZpcSF5kzLpeEqLCm4ORoT2BZJkZPmGlcJysr4GgyTkoXQIasriC2l0VpSFBWOMAtDpWGQlPjWg5BQVsFi4NAOLEZZZfIBQqOUQNQucghBUeUkpQlRORbomv5q1V/nhfGMtjCnC9syVqdt30ZKhWcb6m5ZmhNc01MoKRmm+ZRMsNR2mW+4tByb79/dRWtt6l5y4/18balB07WI0oLFpocU4Cm1r8z8h+tD3lxt8We3dtmNMgJbEBeCqhK8e30R37HYiVKuLpq/exajqyfJ2z/t3X1ex8WTwGHu/D8AXz3E92Y4QRwlJvksi8jewRk4FpcWwjo5HFHVxXFSiMdejie9hBrD5y60xqkliLdGMZ5jpAUsaWoWsqLieze3OV/belpK0nAlo7RimGbYSvD9W72pI9fkczRci9Cx+PLZ9kMvVF6WtP0H4akoK/lgbcjmMCX0jELPQsPlXMdnkOZG4royzKLAkcw3TDL3G1cXkEIQZ0VdAJVhKxilZqEyRi8Wa4MIxxJ86Uybu70YKQVZZnbYCw2Pa0shncDh4mKDhmvxzz7aYmeUYkm4PB/Q8izQgpu9MU3HInAln2wbh7hLCx6jpGKcSG5uj8mqEoGgqM2BWoHNhW7IdpQQOg62yCi02ckvt1zmGjZxVtJLClqeg8D4HNiqQkqBxlzHs6irvm3SwmgUDeKMwLWwlQlvlUVBUmrIK5QWppreErjCeCm3PIftKKsL9zR5oadS344y1FMhBVKYcJUELi82GaQFu+McITReXdRWaXPa8G3BcsvnrXMdelHG0q5DkZvk/MX5kIZrM85zKq0pqoo7vTHLLf+Jk3eUFZzp+Li25A9+toFUAqsQrMyZmoi0LOnHxVQk73mMrvbiMO/u8zgunhQO8nP4BvBLwKIQ4j/d86MWD06JM7wgHDYm+ayJrcdrHwSXFhr80rVFBrGxVvSUfOrLMUoL3lsbsDVMKLVx8Wp4hmmzOcoJHcm1pS6OJdkZpeTiwY71i2dbtcKloRQ2PYemZ17QtKi4txsTutb0/vsdxYtKT6uBwYSjvni2xUcbI/K6EvXKYoPtkfFivjTvG9OWvOJ81+fNMy2Wmi4N12Kh4ZKXNkstnyTX3ChL4qyaekAstzzGtU/C5ijl6mKIUoIP14fYSvBLV+fZiXKGScEXlhukhWY+dKiqipvbEaO0ZKnt8caZFllVIaUmK2AU52wME3YjyaXFkDMdD2UJPlwzz9WSAiFNSGu15fLacos4L/h0y1S8LbY9BAJdJ3htAd3AIs4KHE+y2rGI0ortOCF0FUJrorwiqcupBSbnU2qNqCpkJepgVUUFdYjRInQEozQnr2CYGjMiTynSqqTUmo6nUEqaZHQFlgUSU+Q3HzpsjTNz1crQaKUSuEoReBYX5kLmQ5szHSOu2Isyzs/5xuNAa851A0P3HZvTiJLSVEAfMM4nG6CFhsfXr8yTl0Z8cBAbnSUl4PrKk50FnxWHeXePayE6Thx0Zwdo1L+zt3cGwF85yUbN8DgOG5N81sTW3sG517VsUMsAPG2QjtJiyt4Ypzktz2EnzslL43vcz0zM/Ww3oOWbYrOz3YC0KFFSsDPO6AYOc6ELGC9eJUU92T+gvt7pRVxfaT1Rp8mS4uEAOLDc8ijqicO1JB9uDPnKuTYfbowZZwWerVgIHfyaOTU5yk8WzLLStDxFqQVFVXFhLqDl2SRFhS0lc6GhOs6FJnR1ad5w5IdpwZXFEIHxQv7JvV3irGKp6XJtuUFRQpSW9JOM+YZNWQr6OmOh5bIzTujFBYOkQEll/IpDF1eKui6lAKlJS82XF0McpZAI+kk+NbKptAkz7UQVGs21pZBb2xFbwxRLSXRlfJYrAQ3PRlcVUgkSrXEthSWEqa5G0/Is+rFGKEnTVTQ8i51xToA5jRSlybMIpSgzWGgZbayq0ggNUV5S5MbhbrHp0nYt+mnB9jij4ZvkellpWg2L15ZbfP3KvHnmg5QkLxEaNJKO75AUJb3YyLInecncSovrKy1W2t6hN0ALDZcPN0ZobcZFL8rMqaEey8c5KR/23T2Oheg4cZDZzx8CfyiE+N+01p++wDbNsA8OG5N8nsTWhM46Sgvavj093j7t5DE5raz3YzqBhW8rbm5HrLY8xpkp+DrfDfhLX+7w4zu7Ru7aNdTDrKhYaXtUmoc+30TSI80L/uD9dT7eGNP2bM51PfJyf52mtKho+UY5NcnL6feUlPzi1YWpC1deat65NMcXz7b5we0+WhupC0tKlBQPhRXOdHy+8/E2oWtzfs4ndI18eEVOw7G5ttxgrR+zNUz4dHtcq9mal3yUFlxdanJjY8hP7vYZpQWLDRdLSraGmfElrkp644yVms5bahjFhtZpJmgTjhvEBY4USAlnuz5JVpIUJZ5tEboW9/sp/TgHJkwgxW6UEVHh2cY9b3MEji2RyiHJSmwpkaIiyTSeK8kro2HkWRZLLZckr3AdSVGYEE9ZGVG6KDMeBbYlaNS1Bmldhp2UJYWumHdM8nyQFKRFYRbt2jdjuekRehZOXFBUFVFa0vQsvnyuwy9cmqOoKgJHcXXJ5JVu70SErkXgWnQChz+/2+OPP9nBVZKvXuiy0jJS7l9YOXhi3bsBykuzcMeZkTdpexZvrhov8+OmkL6K+YTD4FAJaSHEv6613gUQQnSBv621/pdOtmkz7MVhY5LPOxCf5eQx+ZtSQ6gUniW4NB8wiDMWGh4aePf6sskTuBbf+XibXpTR8m1W2h5KSr5yrvPQRC8F3NuJ2a6F5OYDs1P/YGPEayut6alm71EcIHAV8Z6q1oWGOz2eT1y4AscyyfKGzS9dtaZCgw3PfmxSGMQ515YaU8G03/3JOrYyzKeVlscHG0NatVTFB+sjPtoc8dae3Mhk0XRtwSiDnShnpW0UPodxTsuz2bFz4rzCkhh6ZaXpBjZLLQ9LSe73E5q+Q1kZiemPN8d1bkbxpXNtqrLizs4YEFRVyUbfCBTaylAzz3cDPt4aGVlwJfEdQVVpKq0oK0FY10KUlaDQFVeXmiy1XDaHqdGhqu3YupZJZpda0wkclBQM44I3lkP6ScZOlFMUFaM0J86qetE1YUCpoR06vHN5jg/WhswrhfYq3EhRVYLrK00jj1KWtYaW+bu3z3d5u2a0rQ8SvvPxNq5t8e7rSwhhHONCR7HUchnE+fQZPwn7+Yu0ffuh9+Vp4/2oeBXzCYfBYWaMhcnCAKC17tWucDO8QBw2Jvm8A/FZTh6Tv9krQdzyLSwluDQfYis5veZyy+ObbyzvmzAPXeshSY8P10c0XJuBNvLUtmVCOB+sDfjG1YWHiu325lqaTSNilhbVviGxvX3k24qVljGB2bswTMJkf/TBJq6tONvx0ELwjavG0OfOTkxSlKw2PXRd9RzYRv77ezd3ePt8FyGgH2d8vDWu4/mmhqQfZSw1PfpJwXzD5ZdfW+D2TsQHG+YU0wksmq7NfMNjmOR8YaXBmXbAH7y3QakrQBO4kssLTZabHr//s/vMhw7LTVM1vbMd0fEthmnJIM7ZsBNA0/YsltseO5EhsoaOotTGb9pvuoxSo9201HSIsxJLSUZxgesYccOGa7Ew52ApE965stTgj29sGfJB6NDwbO73E5KiYphlaCRJUeBbNp3Q4ecudmn5FotNlxJYaHhsDrI6uS+p0Hy4MeLtc11avv2QB8Jqx58u1Eqa4klREyRsJegGzjPRPl8EhfRVzCccBodpXSWEuDCxBRVCXORlmZp+znGYmOTzDsRnOXlM/mZCgQXDmJnoJl2YDw/FoHr08/34Th+JZlDH0VdaPrbC2GI+0qb9TjxRVvDtG1ssNd2pV/UkuW52nRV5WT3WR6O04Ee3d7nfj7EtweYw5fb2mLLS/MLlOc51Qy4vNIhzQ0UdZxWX5kO2xxl5aYqghDCG8RM5DxHDziiui8UsAjenKDUt3+HL5zpcWWxweztiN85MMV+tKJrkBaFr5KJ9JUgriS2g4ShanuLDzSFVJegnBZaQNF0LKSWbwxTLUrR8m7I08fw4M7vytm+DX/GzuwOSXDMX2gihqUrN2Tlnqv/kSM0ozchKZdhY/ZhPd2przCgnLUsansV79wdUQNuzaXiKcSZoVTaakqoU+K7NxfmQi/MBa/2UuYaDLQUCwfXVFlIK1gfmVBM4psr7x3f6FJVpt5KCjWGCFIKlpvvQJsSxBFFWPnOY5qDxfpzU8Vctn3AYHOaT/nXg/xNC/GH99b9ArYY6w4vBUQfp8wzEZzl5TP7GtSQX5wPu9mL6ScHrK82pxeGzMKhW2h5JXrHQ9LjTi8iKgru9jKKq+GhjxNevzE9/99EdYJSZPis1NBZCelHOn37a49pSSDdwpp9r0obJyWPi6zBKchwlqSqT4A5dm61hwg9u7/Llcx1eX2mxOUzoRSVztYvY2boOJC9Lfnx3l6rUxGlBJ3RI8oqzHZ+dcVrH2Su+dLZN4CpubAy5tRMZuQbLo9CG49/xLBwV0A0cc4JxFSIv2RplrA9TkszE/TuB2bFv1f4EDUfSH1csBg5JVhB6Fuc6xrhpe5RhW4JRWhC6FhpT/GcpwWrXw3McOoGgEzjc2o440xHEecHWOGVnnHF2ziewFPcHCZtRxpX5cCoxMkgN48dRioU5d2poU1Sa0JXsjrO6xkGw2HD5aHPIuTkfzzEFbC3fxlaCH93uc2cn5tpSSMu3yQqTmLaVMTc6aBNyXON9vuE+dcy+alpIxw35tF/QWv8upqbh/wT+DvDzWuvfO+mGzWAwmbQmTlKTZOyzHHsn1/rh7d4TrzHPFc93AAAgAElEQVRZWCZCYbaST53E9/5NpeHqUpNf/8oZ3j7fpVGHiia7eiGMBMJEhuMgfOVch3FakJW6rqyO2Y1S3lxtTZPR01xDvQOcYGOYgjC7ZCGMgmboGr/nR9vwaB+P04I7OxEbo4SGpzjXCfAdieeaiU5KgW8bY5dKa1xLojH1D4Mkpx8V09qQpChZ68fM1ZLarq34jbfP8m987QLd0MFRJr9yY2NEVlSoWgLDlpJOaBM6JieitWYY5WyNUrJSkxUl76/32Y1ybEvRDhwGccHd3Yi80AglGCcF21HGxiDGdSyur5hcQl5UbA5SpJSstj2Wmy4rLZ9f++IqriVQSvLJ1hCkYLHlsdT02RpldHyL3bFJ6Lccm7ZrTJUWmyakFTpG1Xe17WFJozD7xpk2XzrbQlcwTEveONPi179yhm++scy7ry/zxmob3zYnHNdS5CWM85KWZzFKTV2Ha5kQWJybE4Ks5berypjqrLT9Z04gP2m8D+L8wDF7nO/lq4rD9maJqYj2gDeFkXb+pyfXrBkmOC7NlaPUPzzLyeOgv3mWuO4oNfzzpabL7Z0xW6OMCwshXznbYb5hNI6SvJz2w6M7wH6UYSnJYtMkKOO8pOEan4hH2/BoH7d9mw3H4v5uwrXFBkIIuqFL23dYaDpEWclP7w8A+PK5NpuDhN0op+UZoTjPUXTCECUllxca3O8nbI8zznV8vnpxjrfOd6fPwrMVN7cTUzBYaaTQuLbFMMnZHmX8/MU5fnSnT8O1+LA/Ag1N32akK6LMyD7sjFKKEmwlGGWaIslxpaBR14S4tuR816fhWUgxZjfOWG157MT5VJrdsRV/crPHXFAL9fUTbClZbQdozGK02nLZrE8eUWHc7CwlONMOGCcljq2M5HpestL2uLLY4MqiqVG5vFCx2jY04Rsbw2mY795uTD8y/sxpUZGVFR3PnCAGcTZ9zoJH2UYmcX4cu/X9xu7TxuyrqIV03DiM8N6/B/xHGNnsHwC/CHwb+ObJNm0GOL6E2csczEfNY+xdyC7OB6y0PX58t8+bq62HTFP29sOjuZaGZ6SWg9rA3bcVo7QkdB9vw94+NpaSJcM4r3f8RnU1yksWQpeWZ1OUmquLjWkYYrLgAoyzkkvzIUIIbm6Pa1+JkltbxhugXSdO994zzkqWWx53dyOyUnN50SPNjYnM1aUmG8OUOzsuljUmzytGqXHLsyVIUdcCZEYQz0hVC9K6ruTKUgPPFkR5hW0Zb4N+ZHIuncCwrBypGCcZrU5A6Fms7eYshS6Vhu1RgpRwphvguRZL0tBRN0cpUVISOorQVXy0MWSh4fHNN5aIc839fkTgWCR5yTibnASKKTlhUtR4puOzNUrpx4Y+faYTIgX89F6fuP4MlhDsxsX0Wb6IyfdpY/ZlaiG9qHDWU8NKmIXhF4BPtdbvAj8HbB57S2bYF4+GS+DZONJRVuyrMGksLU8Wq3UsPslLtDY1CBMm0QR7Q17fvrFVm+M8ONK3PYu7u9FD132SCNpb57t84+oCSorpPZuezTgtaHnWY22Y9PFEStxSgmtLDV5bbnJ7J6IXpZzv+lxaCOhFOWe7wUNtM4uQUem8vtJCSeODsNR0ubUTcWNzjG1LXEvwnY+3+Qc/vDtNooIRGZTSeA8ErjGd0Wher6t1ry426Cc5jlSEriQrSiTm96OsZHuUcHkx5GzXR0lJy7O4utxASsHOOCXOSvKiNDUGLSN1PnGea3oWQsEgNXLp28OEotJ4tgmvjbOStmszHzoM45KF0EPrCksIHEuy3PKMT0XLZ7ntstLyeGOlyTsXu/STnKysuLIYstB0p0q6e8M0gzjnG1cXuDhv2ubbxtgHLegGLuOsZK1vnNscy+RnDsJhQqfHMWaP6708Kl5kOOswnyTRWidCCIQQrtb6PSHE68fekhn2xXFxpI+zEOegncv6IOFHd4yw2UQPabnlHcigejTk9fHmiKQWVpvsFs92fX56f/hQcdtB/fDoSaLt27x7fWlaCLe3DZM+Xq/ZSWDi7r/yxgpJXhpdn5rxlJd6aiQzwZOUaodJjmsJw7KyDBPHtzX9pDCaVV2fbuCw2HD4YH2EEIK3z3VQ0lA0J8n8stJ88/oS/+BH97i5E9ENbOYaLoPYJJWLSlOUFT93ocvd3YhPNiMjZ60sXl9tsd5P+Nn9AZWGrKi4ON9gEBd1XyparoUS5gQjpeT15YDtcY5UgmGcEXo2Ld/mG1cXeO9en5/dHzIXOnzzjWVE7X739mqLS/PhdDxpbbyRJ6q7P7zdo6z0VFjRt1UtT1Ls+6y+erFLnBd8shWx3PZo+za7UcY//qnp272mU3vH5XF5IjyN9feyahdeZARAaH0wK1UI8XeBvwr8x5hQUg+wtdZ/6Vhb8hx455139Pe+972X3YwTw3EcI/e+OHsH81HZFwddZ5wWfOu9jVr/yIjnjdOCd68vPVSc9Oh9xlnxkKnKze0x47Sk4SoWa+e3SQ7hwpwxUjnu4/QoLfjWe+tIYa5tKp3VvpPchAXl16cDKcQ0kbn38/3gdo9+lNXmMcY4RmvNOCuYCx2uLDamYS0w/HCxz2f74e0eDddie5zyf/3ZPZKixJYmyf7FM+2aNhszHxppiQ/WhthKcX7OQ2vB9jih7TustDzu7MaUleZsx+fObkxaVHQD2+RLaqaT5ygEknFmKrqhthEVAoQkzQuUkmSFkUcHsPc8P4CdcUovyqeL6tYo5X4/oeFaOJapBxnV8iKTIre9+HB9yI2NoREI1II7u5HxuxaYJPo+SejJjnpvO5K6GPAkQlEvg600GQtij4DYo2P0KBBC/KnW+p39fnYYm9B/tf7nfyGE+BbQBn73yK2Y4ZlxHBzpp+2EDrvrOmjncnN7TOhaRmkUaHkmjPWjO7v86psrT7zP+2tDLs0HrA0S4qxEoEmygnGaM0pyEMblbbXtP+REd1gctsbi+krrscll7+lqlBYmX5BXhLXfxPtrQ1bbPl8533msrwFubA65tRMTOOZlzkuNEoKWb04fT3quj1JrTUW3xy9fm2djmJEWJaMk59xcQJKX3O5F5JUx8Gl6Nr6jKEvjPHemE3K243G7FxG6irXdhPVBwrmuT5yWOLbk4lzIhfmA7XHGx5tjmh5c7Bq21XLLZ2ucUlaatmeKBjdGCev9mD/6IGGp7ZFkJUtNt/a9qPhoY8y1JaNE24sy/uSTbZKsZLnts9Jy60lf72sWD2Zn/u0bW1S6Yr2fojFEgcsLIaXmIdOpCV50HuBl1C68SCmOI12x1lua4ZTioMF82OPqQS/gbpSxUDOJHtxTsjXKDryPq0Rdg9AgrAuc8lIzqkXx2p6a7uT3MpQOg6OEGp4WKri/awyF2r7D5tBIVISuInDVvovVasdnY5hQlGNGiTkdRXnBQsOl7duHSsibuLIxGbq2BOfnAuK8oiVsfu58lzu7MVvDjLfOd4jTip3YVF/nZckwLZkLHc62PfJSmzqFjs+Zjg9CYEnFVy40mQ/dOsRTsdj0TMJ5mNCPcxqezVfOG02sOCvojTP+/N6AqjIT1MYwptKCbmCRFCU/vT8kcIzPwlzoEGUFa/0ES0kWGgq0+SxXFkO+sNygOiBw4dmKfmSEG/16vOSlNtpR+0z6p1XD6Ch4keGsz06vzfBcOOyu66AXsBM4jNJqemIAGKUVneDBgrHffex6kO+1e/QcC8uSfPls+6Ej9FF3gkeJ0T7tdDVpuxCCi/XLODnS74eGa5LUUgi++8kOmaq4MB+w0jJ6UnsT8ge1eS50uLbUmOY+ri01piGoK4shZS1U1/QEtqWwZMIgFoyzgp1xSlZbYS42vWk4bLUTsNBwqCrNQsN9RFLEJJo7taQIGGtUIQRxbk4sFYJRZixIu4FNWRnRw5WWx43N0TQvszE0Nqld32aYFlxbbJAWRppDSYmn9ufE3N+NubIYstZXOJaaVrTfHyT80uL8vpP+adUwOgpepBTH53px+KxXOB4Fh911HfQCKin4v39wF8uSdDwL1zbJ0q9dXjrwPllR8YXlBrYyUgi+rXhtqcHN7fFDvxtlJXd6Y/JSH/p5HTXUcNDp6ll3pgsNl1+4PDfVgmr7zoFtf5RauzFMidICzf6J2LAWEvxke8y9XkSuNZaQzDc8RmnOMMlpe7DSCrjdi1louMwFNlpDPymmbXnSpGMUcD3W+gmjtEAASmjW+xlvnW9jK0FWGnmOCSNu0k9xZuiuDc+eMrRsJdiNjET7QcY8E3aTFIKPt8aEjqLtW/uaTu19dqdNw+ioeFHhrM9Wrx0Bx8ls+CzgsLuuJ72AAP0456sXu9zcGrE5SvGdkl99c+WhZPR+91FSstzymQsfnDCSvOTCXDClC5aVnurvf2G5MaXwPe15HWeo4ag7071jbKnpHigGuF+bK62NXaslsZUkzkv+n5+usdTyprv9CdvqR7d3+WBtQFaWeI4kySrmQ4cvnWlzZzdia5hRas3bFzr4tqIC5B66LDx50tk7Ue9GGZtD8GzBXOgga5MdJQS+o0iLivN7nptnS4ZpiZSCt853GaX5NFx1GGOewLG4vtriwnw43RjYB5hOnUYNo1cVT2UrnQY8C1vpRTMbTgOe5yR1lP589D6TStn9GFBgQgzvrQ2wpeDcXDCd2A/zvA7D0joKjtJHzzrG9vpjCAkCyTDJ0WhcyyJwTHJ+7+f4we0ef/ppz8iDW6YyXCAoK+OXsdLx0XUi96j9sPdzRFnBe/cH3O8nOJZhX2UFrLYdLi82UFI+9Ny2Rikbg4Sz3YBuXQV9mPse93ObYX88F1vps4qXWeH4vDjucNij19svdPG0v7u1E3FpPmCvg+yT+nO/3d1eue5HwwGvLTcfivcfdP39+uY4Qw2TnfrkHvd34yf2/7OOsUn/3OlFSA2BYyQ5bCkpq4pPtsYktZeyFPDW+S4C+NqlOd5fH7A1MqY/FSW9OOcLK61pzcSz9MOjEueXFkKTSBYgNHiOqpPsD4fLXltuTuXUj3rfz0uI6FXG7ORwyk4OJ7ETfpbrPfp3H24MGaclr680H9vZ751Mn3Uxe3T3uteg5xtXF2i41mNt6kUZ93YTFpvuQ2GY58FR+ut5x9jev39vbYAQ8OnWGM+2uDQfkJYlu1HBr3/lzFRvqNLw6c6Ytfrrs91garS09zM86XmM0mLqvgZwYS7g6p6FZZaf+2zhoJPD53ZxOK3H1uNe1J50vaysCOsCrf0mkG/f2GJcV7NOxO0+WB8S2KYdk/6cqKc+bz9PnldZVUYZU0hAs9o2UtR7BdkmC8jNrTEaQegYR7Sj3PdJE+hRw2eTNvfjnEGco6Tk61fmn+pYtvfvXUuy1o/5aHNMVVVcXjDOdCYvUXF18YHw4NP6+aBxD6bI6tPtiDQvScuSsoTrKy2+fnX+WN6L4zz1zgglz4+DFofDaCt9JjE5th5FmvpVwHFrJO13vbLSvL+2v37LZHIZJTlt33Dwb24bbf3XlhrklT5Q+rjSsD6I+dZ760fShJk8r16UG48GV3F5wXDpJwVRW6OUtb7JT/zgVg+tBU1PkdTJ3fX+4e57kH7NUfq/4RoP6ru9mM1ByjAtiNKCP76xxfogOfRnNr4IFklWstr2cG05VTA92wmm4avDjOeD5NPv78ZsDhMGcY4lJR3PyIy/vzZ4qqbRYXAS8vOfZcnsl41XeyY8YZxGZsNxF/rsd727u8Zwfb/aADCTYTtwyGs3LoDNYcJyy+P6SuuJRXNRVnJze1y7gOlDM44maLgWS02XxkL4WO5hY5hOufhNT3E7q0iKGI2HpQQ3t8bYlkBqnnrfg2ojjtr/g9gI9a0NEhwlcSzBMCn5zsfbfPON5YdOY/vtgidj9LXlJr6jWOvHU7rvRMHUU2raP08bz0/Lg2yNUgJH4dQLoKHJZtzaiZ5JnmEvjlMX6PMgmf2y8bk9OZxWHEbh9Hmv148LznYfvt5kdzzZOS81XbKiIi1KbGVorPu1Y6965ebQTJBCQOBahzb9edL1JkgLI2N9puMh0GRFRcORFKU2O3QNjiURGPG7p933oNPBUfs/ygoGcYajTBhHIGi66kF4jMPvgq8tNVlu+VyaD7kwF0wd0I7y7A9SEw0ci6J8OMycVxWedTzTxHGeel+myvDnBbPF4ZThuMNh+13v+koTJR8eGnsnkMm/Ly2EWEqyGxmXtSdJUkwm0ygr0FRkhfGbhqO/0E+anBuuRTdwpm3yHAshoelaVFqjNWS1PMTT7nvQBHrU/g8c4z7nWA9OOllZ0fLt6f0P65R3HM/+oMVttePTDR36SUFalKRlwTgtcW2L87Xg4fPgOGWuX5Zk9ucJs548hTjucNij15vsZGH/Yq8f3d5lEI8pKo0lBcstIzx3UFHS/d2YSoOu4NLCA2nno77QT6I43q8VRgPH4tK8xaX5cKoMmuTGH8FISqun3vdpxW5P6q/9EqOrHZ+f3uszTEqartGNygpTe/AyjGOeRhH9i68v8/++v0FvnGJJwUrbY7HpTqmwT8NBSeLjlLf4PEhlvGx8btlKMxyMJ73ko7Tgh7d7DJOcstIoKWh6Nm/VftFPu+ZJMcSexsI56n0Py4Q5zGdaHyR85+PtWn/Ipu3b02KxozCgXhTD7llZQMclC3/S7ZzhAWZU1hmODc9LpT3JF/pp/P2TuO9+/XG3F/Hx5ojAtaaGR5Mivye17TCT/qtem/Oqt2+GxzGrkJ7h2PC8IZCTZIgddO2Tuu+j/bE1SvizWz2UklyYDxilFd96b4N3ry89tW1PqwZ+1av6X/X2zXA0zBaHGY6Ez4Nm/lHwaH98sDbEsRRt30YKua/h0X44zOL1qvf9q96+GY6GGVtphiPhuKm0T8Ik1PK8RvEnjUf7Y2uUYkmYC93p7zRcyW6UHXCVZ7vXSfX9s+JVb98MR8NscTgGnJaJ7DjwIirLT1P166P90Q0cOoGLZz/Z8Oi47vWqVfW/6u2b4Wh4KU9NCPHfAn8ZyIAbwF/VWu/WP/trwL8LlMB/qLX+vZfRxsPi8+gLcdKV5aet+nVvf6y0fb713gaDpKDhSkZpxTgtHjI8Oq57vQp4kgruDKcfL+vk8E+AL2mtvwJ8APw1ACHEm8BvAV8Efg34n4UQ6olXeQVw2AKmGQ6P01z9utzyePf6Ep5tvLM9W/Lu9aVDCe2dNpymE94MR8dL2dpqrf/xni//GPgr9b9/E/jbWusU+EQI8RHwNeDbL7iJh8ZhGBozPvbR8KonNp/2PJdb3oHJ56Ne71XFaTvhzXA0vAo5h38H+Ef1v88Ct/f87E79vccghPhtIcT3hBDf29zcPOEmPhlPK+Of7a6OjuNMbB53Pui4n+dpHh+n+YQ3w9NxYouDEOL3hRA/2ee/39zzO38dKIC/NfnWPpfat0pPa/07Wut3tNbvLC4uHv8HOCSeNpHNwk5Hx3ElNk9i4j3u53max8dM3+izjRN7ilrrf/Ggnwsh/m3gXwZ+RT8o074DnN/za+eAeyfTwuPB0wqYPguFQS8j7HEcideTCHtsjVKirCDJK3xHsdR08W31zM/zNI+Pmb7RZxsvJawkhPg14D8DfkNrHe350d8HfksI4QohLgOvAX/yMtp4FEwmsrfOdx/b4Z723dUs7PEAo7RgY5AQZRWhoyjKiptbY3pR9tx+GntxWsbHjLr62cbLyjn8j0AT+CdCiB8IIf4XAK31nwN/B/gp8LvAv6+1Ll9SG48Fp70waBb2eID7uzFnuwFCQFZqHEuiEdzbTY7VT+M0jY+DNkYznG68LLbStQN+9jeAv/ECm3OiYZPD6ua8qpiFPR4gygq6gY1nKzaHCVFWEjrGwvN5/TRO6/iY4bOLz/0IfBFFbK9a4dJR8KrTSg/CfhPvfMN95o3AA6MjxcV6gZmojh5HO2eY4VXCq0Blfak4zWGTF4HPUthjteNzbzd+5vzJae+LGWY4Cj73i8OMq30wPktJx+fdCHyW+mKGGZ6Gz/2oPs1hkxeFz0rY4zjyJ5+Vvphhhqfhcz8Dzrjapw/PSiCYbQRmmOHw+NyHlWahgtOF56m7mOUMZpjh8JjNgMxCBacJz1P1PKONzjDD4TF7K2Y4VXiVPaxnmOGzhM99WGmG04XTLDcxwwynCbPFYYZThVneYIYZXgxmi8MMpwozAsEMM7wYzN6oU4DT6hR2UpjlDWaY4eQxOzm84jjNktkzzDDD6cXnd/t5SjDz6X02zE5bM8zwfJi9La84TpNk9qsyIb8Ipd0ZZvisYxZWesVxWqibr1L4a6a0O8MMz4/Z4vCK47RQN1+lCXmmtDvDDM+P2eLwiuO0UDdfpQn5tJy2ZpjhVcbsbTkFOA3UzVdJ8XSmtDvDDP9/e/cfq2VZx3H8/YEDB/lRRPwQBYQZxiiREJkFLdkIwa3IFk5yi8qiH2i5xsqi9WOsrbWay5YYOiZlgejmYmomMprBokACBMXBQgVBsCUIAYeOfPvjvk48nvs5PRzOj/s5z/m8trPnfq77x/U9165zvue67+dcV9t55GDtoppuf3WV0ZZZNfNPi7WLapvxtCuMtsyqmZODtRv/QjarHb6tZGZmOU4OZmaW4+RgZmY5Tg5mZpbj5GBmZjlODmZmluPkYGZmOU4OZmaW4+RgZmY5Tg5mZpbj5GBmZjlODmZmluPkYGZmOYUmB0mLJIWkwem9JN0taa+kHZImFRmfmVl3VVhykDQS+CjwSknxbGBs+loALC0gNDOzbq/IkcNdwDeBKCmbA/w6MpuAgZKGFxKdmVk3VshiP5I+DrwaEdslle66FNhf8v5AKjtU5hoLyEYXjBo1quOCNSvYiYZGDh09xckz2Qp7wwde5CVPrcN1WA+T9DRwcZldi4HvADPLnVamLMqUERHLgGUAkydPLnuMWVd3oqGRPYePU1/Xg/71dTQ0nmXP4eNeE9s6XIf1roiYUa5c0pXAGKBp1DAC2CppCtlIYWTJ4SOAgx0Vo1m1O3T0FPV1PejTqyfA/14PHT3lJVmtQ3X6M4eIeC4ihkbE6IgYTZYQJkXEa8Aa4DPpU0vXAsciIndLyay7OHmmkfq6t/+Y1tf14OSZxoIisu6i2salTwA3AHuBk8Dnig3HrFh9e2e3kppGDAANjWfp27vafnSt1hTew9LooWk7gIXFRWNWXYYPvIg9h48D2YihofEsDY1nGfXufgVHZrXO/yFtVsX619cxdtgAevXswYmGRnr17OGH0dYp3MPMqlxTgjDrTB45mJlZjpODmZnlODmYmVmOk4OZmeU4OZiZWY6yfy3o2iS9DrzciVUOBv7ZifV1VW6nytxGlbmNKrvQNrosIoaU21ETyaGzSdoSEZOLjqPauZ0qcxtV5jaqrCPayLeVzMwsx8nBzMxynBwuzLKiA+gi3E6VuY0qcxtV1u5t5GcOZmaW45GDmZnlODmYmVmOk0MrSJoraZeks5ImN9v3bUl7Jb0o6fqiYqwmkn4g6VVJ29LXDUXHVC0kzUp9Za+kO4uOp1pJeknSc6n/bCk6nmogabmkI5J2lpQNkrRW0p70+q621uPk0Do7gU8Cz5QWShoP3Ay8D5gF3COpZ/70bumuiJiYvp4oOphqkPrGL4HZwHhgXupDVt701H/8vw6ZB8h+z5S6E1gXEWOBdel9mzg5tEJEvBARL5bZNQdYFRENEbGPbJnTKZ0bnXUhU4C9EfGPiDgDrCLrQ2YVRcQzwL+aFc8BVqTtFcAn2lqPk0P7uBTYX/L+QCozuE3SjjQUbvNQt0a4v5y/AJ6S9KykBUUHU8WGRcQhgPQ6tK0X9EpwzUh6Gri4zK7FEfH7lk4rU9YtPiP8/9oLWAosIWuLJcDPgM93XnRVq9v2lwswNSIOShoKrJW0O/3lbB3MyaGZiJhxAacdAEaWvB8BHGyfiKrb+baXpPuAxzo4nK6i2/aX1oqIg+n1iKRHyW7JOTnkHZY0PCIOSRoOHGnrBX1bqX2sAW6WVC9pDDAW+FvBMRUuddImN5I90DfYDIyVNEZSb7IPM6wpOKaqI6mfpAFN28BM3IdasgaYn7bnAy3d5ThvHjm0gqQbgV8AQ4DHJW2LiOsjYpek1cDzQCOwMCLeKjLWKvETSRPJbpm8BHyp2HCqQ0Q0SroN+CPQE1geEbsKDqsaDQMelQTZ76rfRcSTxYZUPEkrgeuAwZIOAN8HfgyslnQr8Aowt831ePoMMzNrzreVzMwsx8nBzMxynBzMzCzHycHMzHKcHMzMLMfJwawCSSfS6yWSHqlw7B2S+rby+tdJapd/EJR0vyfxs/bg5GDd0oXMmhsRByPiUxUOuwNoVXJoTxHxhYh4vqj6rXY4OVhNkTRa0m5JK9KEf480/SWf1gb4nqQNwFxJl0t6Mk3q9mdJ49JxYyT9RdJmSUuaXXtn2u4p6adprYEdkm6X9DXgEmC9pPXpuJnpWlslPSypfyqfleLcQDYNfLnvpa+k1en6D0n6a9M6IpKWStqS1hf5Yck5fyo55oSkH0naLmmTpGHt3+JWq5wcrBa9F1gWEROAN4Gvluw7HRHTImIV2aLst0fE1cAi4J50zM+BpRFxDfBaC3UsAMYAH0j1/DYi7iabI2l6REyXNBj4LjAjIiYBW4BvSOoD3Ad8DPgw5ScuJMX9Rrr+EuDqkn2L0/oGE4CPSJpQ5vx+wKaIuIpsPqIvtlCPWY6Tg9Wi/RGxMW0/CEwr2fcQQPoL/kPAw5K2Ab8CmuaCmgqsTNu/aaGOGcC9EdEIEBHN59cHuJZsMZ+NqY75wGXAOGBfROyJbIqCB1uoYxrZWg9ExE5gR8m+myRtBf5OtshUuecMZzg32eGzwOgW6jHL8dxKVouazwlT+v7f6bUHcDQiJp7nNZrTeR6zNiLmva3w3HxTlZSb2ps0ueMi4JqIeEPSA0CfMof+J87Nj/MW/nm3VvDIwWrRKEkfTNvzgA3ND4iIN4F9kuYCKFxHH6AAAAELSURBVHNV2r2RbKZUgFtaqOMp4MuS6tL5g1L5cWBA2t4ETJX0nnRMX0lXALuBMZIuL4mxnA3ATenc8cCVqfwdZEnuWHqOMLuF880umJOD1aIXgPmSdgCDyBYdKucW4FZJ24FdnFuq8+vAQkmbgXe2cO79ZLNf7kjnfzqVLwP+IGl9RLwOfBZYmWLZBIyLiNNkzyweTw+kX26hjnuAIencb5HdVjoWEdvJbiftApaTJTOzduVZWa2mSBoNPBYR7y84lDZLH7ftFRGn0yhjHXBFWnfarEP5HqRZ9epL9rHYXmTPH77ixGCdxSMHMzPL8TMHMzPLcXIwM7McJwczM8txcjAzsxwnBzMzy/kv9F6QvaDqVsgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fix, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(estimated_gain, nhefs.wt82_71, alpha=0.15)\n",
    "ax.set_xlabel('predicted gain')\n",
    "ax.set_ylabel('actual gain');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 13.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 13.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 4 steps to the method are\n",
    "1. expansion of dataset\n",
    "2. outcome modeling\n",
    "3. prediction\n",
    "4. standardization by averaging\n",
    "\n",
    "I am going to use a few shortcuts. \"Expansion of dataset\" means creating the 3 blocks described in the text. The first block (the original data) is the only one that contributes to the model, so I'll just build the regression from the first block. The 2nd and 3rd blocks are only needed for predictions, so I'll only create them at prediction-time.\n",
    "\n",
    "Thus my steps will look more like\n",
    "1. outcome modeling, on the original data\n",
    "2. prediction on expanded dataset\n",
    "3. standardization by averaging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame({\n",
    "    'name': [\n",
    "        \"Rheia\", \"Kronos\", \"Demeter\", \"Hades\", \"Hestia\", \"Poseidon\", \n",
    "        \"Hera\", \"Zeus\", \"Artemis\", \"Apollo\", \"Leto\", \"Ares\", \"Athena\", \n",
    "        \"Hephaestus\", \"Aphrodite\", \"Cyclope\", \"Persephone\", \"Hermes\", \n",
    "        \"Hebe\", \"Dionysus\"\n",
    "    ],\n",
    "    'L': [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n",
    "    'A': [0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n",
    "    'Y': [0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0]\n",
    "})\n",
    "\n",
    "df['A_x_L'] = df.A * df.L\n",
    "\n",
    "df['zero'] = 0\n",
    "df['one'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>one</th>   <td>    0.2500</td> <td>    0.255</td> <td>    0.980</td> <td> 0.342</td> <td>   -0.291</td> <td>    0.791</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>     <td> 1.665e-16</td> <td>    0.361</td> <td> 4.62e-16</td> <td> 1.000</td> <td>   -0.765</td> <td>    0.765</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>L</th>     <td>    0.4167</td> <td>    0.390</td> <td>    1.069</td> <td> 0.301</td> <td>   -0.410</td> <td>    1.243</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A_x_L</th> <td> 9.714e-17</td> <td>    0.496</td> <td> 1.96e-16</td> <td> 1.000</td> <td>   -1.051</td> <td>    1.051</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = sm.OLS(df.Y, df[['one', 'A', 'L', 'A_x_L']])\n",
    "res = ols.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make predictions from the second block, with all A = 0, which also makes A x L = 0.\n",
    "\n",
    "\"the average of all predicted values in the second block is precisely the standardized mean in the untreated\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "standardized mean: 0.50\n"
     ]
    }
   ],
   "source": [
    "A0_pred = res.predict(df[['one', 'zero', 'L', 'zero']])\n",
    "print('standardized mean: {:>0.2f}'.format(A0_pred.mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make predictions from the third block, with all A = 1, which also makes A x L = L.\n",
    "\n",
    "\"To estimate the standardized mean outcome in the treated, we compute the average of all predicted values in the third block.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "standardized mean: 0.50\n"
     ]
    }
   ],
   "source": [
    "A1_pred = res.predict(df[['one', 'one', 'L', 'L']])\n",
    "print('standardized mean: {:>0.2f}'.format(A1_pred.mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 13.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We repeat the steps done above, but for the `nhefs` data:\n",
    "\n",
    "1. outcome modeling, on the original data\n",
    "2. prediction on expanded dataset\n",
    "3. standardization by averaging"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 1: Outcome modeling on the `nhefs` data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: For convenience, I've moved `qsmk` to the end of the variable list, next to `qsmk_x_smokeintensity`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = nhefs.wt82_71\n",
    "X = nhefs[[\n",
    "    'one', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2',\n",
    "    'qsmk', 'qsmk_x_smokeintensity'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols = sm.OLS(y, X) \n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>one</th>                   <td>   -1.5882</td> <td>    4.313</td> <td>   -0.368</td> <td> 0.713</td> <td>  -10.048</td> <td>    6.872</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sex</th>                   <td>   -1.4303</td> <td>    0.469</td> <td>   -3.050</td> <td> 0.002</td> <td>   -2.350</td> <td>   -0.510</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>race</th>                  <td>    0.5601</td> <td>    0.582</td> <td>    0.963</td> <td> 0.336</td> <td>   -0.581</td> <td>    1.701</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_2</th>                 <td>    0.7904</td> <td>    0.607</td> <td>    1.302</td> <td> 0.193</td> <td>   -0.400</td> <td>    1.981</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_3</th>                 <td>    0.5563</td> <td>    0.556</td> <td>    1.000</td> <td> 0.317</td> <td>   -0.534</td> <td>    1.647</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_4</th>                 <td>    1.4916</td> <td>    0.832</td> <td>    1.792</td> <td> 0.073</td> <td>   -0.141</td> <td>    3.124</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_5</th>                 <td>   -0.1950</td> <td>    0.741</td> <td>   -0.263</td> <td> 0.793</td> <td>   -1.649</td> <td>    1.259</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>exercise_1</th>            <td>    0.2960</td> <td>    0.535</td> <td>    0.553</td> <td> 0.580</td> <td>   -0.754</td> <td>    1.346</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>exercise_2</th>            <td>    0.3539</td> <td>    0.559</td> <td>    0.633</td> <td> 0.527</td> <td>   -0.742</td> <td>    1.450</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>active_1</th>              <td>   -0.9476</td> <td>    0.410</td> <td>   -2.312</td> <td> 0.021</td> <td>   -1.752</td> <td>   -0.143</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>active_2</th>              <td>   -0.2614</td> <td>    0.685</td> <td>   -0.382</td> <td> 0.703</td> <td>   -1.604</td> <td>    1.081</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>                   <td>    0.3596</td> <td>    0.163</td> <td>    2.202</td> <td> 0.028</td> <td>    0.039</td> <td>    0.680</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age^2</th>                 <td>   -0.0061</td> <td>    0.002</td> <td>   -3.534</td> <td> 0.000</td> <td>   -0.009</td> <td>   -0.003</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wt71</th>                  <td>    0.0455</td> <td>    0.083</td> <td>    0.546</td> <td> 0.585</td> <td>   -0.118</td> <td>    0.209</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wt71^2</th>                <td>   -0.0010</td> <td>    0.001</td> <td>   -1.840</td> <td> 0.066</td> <td>   -0.002</td> <td> 6.39e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeintensity</th>        <td>    0.0491</td> <td>    0.052</td> <td>    0.950</td> <td> 0.342</td> <td>   -0.052</td> <td>    0.151</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeintensity^2</th>      <td>   -0.0010</td> <td>    0.001</td> <td>   -1.056</td> <td> 0.291</td> <td>   -0.003</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeyrs</th>              <td>    0.1344</td> <td>    0.092</td> <td>    1.465</td> <td> 0.143</td> <td>   -0.046</td> <td>    0.314</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeyrs^2</th>            <td>   -0.0019</td> <td>    0.002</td> <td>   -1.209</td> <td> 0.227</td> <td>   -0.005</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>                  <td>    2.5596</td> <td>    0.809</td> <td>    3.163</td> <td> 0.002</td> <td>    0.972</td> <td>    4.147</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk_x_smokeintensity</th> <td>    0.0467</td> <td>    0.035</td> <td>    1.328</td> <td> 0.184</td> <td>   -0.022</td> <td>    0.116</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 2: Prediction on expanded dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "block2 = nhefs[[\n",
    "    'one', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2',\n",
    "    'zero', 'zero'  # qsmk = 0, and thus qsmk_x_smokeintensity = 0\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "block2_pred = res.predict(block2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "block3 = nhefs[[\n",
    "    'one', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2',\n",
    "    'one', 'smokeintensity'  # qsmk = 1, and thus qsmk_x_smokeintensity = smokeintensity\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "block3_pred = res.predict(block3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 3: Standardization by averaging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "original mean prediction: 2.64\n",
      "\n",
      " block 2 mean prediction: 1.76\n",
      " block 3 mean prediction: 5.27\n",
      "\n",
      "  causal effect estimate: 3.52\n"
     ]
    }
   ],
   "source": [
    "orig_mean = res.predict(X).mean()\n",
    "block2_mean = block2_pred.mean()\n",
    "block3_mean = block3_pred.mean()\n",
    "est_diff = block3_mean - block2_mean\n",
    "\n",
    "print('original mean prediction: {:>0.2f}'.format(orig_mean))\n",
    "print()\n",
    "print(' block 2 mean prediction: {:>0.2f}'.format(block2_mean))\n",
    "print(' block 3 mean prediction: {:>0.2f}'.format(block3_mean))\n",
    "print()\n",
    "print('  causal effect estimate: {:>0.2f}'.format(est_diff))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 13.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use bootstrap to calculate the confidence interval on the previous estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this runs too slowly, you can reduce the number of iterations in the `for` loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "boot_samples = []\n",
    "common_X = [\n",
    "    'one', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2'\n",
    "]\n",
    "\n",
    "for _ in range(2000):\n",
    "    sample = nhefs.sample(n=nhefs.shape[0], replace=True)\n",
    "    \n",
    "    y = sample.wt82_71\n",
    "    X = sample[common_X + ['qsmk', 'qsmk_x_smokeintensity']]\n",
    "    block2 = sample[common_X + ['zero', 'zero']]\n",
    "    block3 = sample[common_X + ['one', 'smokeintensity']]\n",
    "    \n",
    "    result = sm.OLS(y, X).fit()\n",
    "    \n",
    "    block2_pred = result.predict(block2)\n",
    "    block3_pred = result.predict(block3)\n",
    "    \n",
    "    boot_samples.append(block3_pred.mean() - block2_pred.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "std = np.std(boot_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "               estimate   95% C.I.\n",
      "causal effect      3.5   (2.6, 4.5)\n"
     ]
    }
   ],
   "source": [
    "lo = est_diff - 1.96 * std\n",
    "hi = est_diff + 1.96 * std\n",
    "\n",
    "print('               estimate   95% C.I.')\n",
    "print('causal effect   {:>6.1f}   ({:>0.1f}, {:>0.1f})'.format(est_diff, lo, hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This confidence interval can be slightly different from the book values, because of bootstrap randomness"
   ]
  }
 ],
 "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}