{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:00.212754Z",
     "iopub.status.busy": "2021-02-02T06:52:00.204791Z",
     "iopub.status.idle": "2021-02-02T06:52:00.676155Z",
     "shell.execute_reply": "2021-02-02T06:52:00.677341Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:00.683108Z",
     "iopub.status.busy": "2021-02-02T06:52:00.681497Z",
     "iopub.status.idle": "2021-02-02T06:52:01.766383Z",
     "shell.execute_reply": "2021-02-02T06:52:01.767615Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
    "\n",
    "np.random.seed(9876789)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS estimation\n",
    "\n",
    "Artificial data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:01.773600Z",
     "iopub.status.busy": "2021-02-02T06:52:01.771997Z",
     "iopub.status.idle": "2021-02-02T06:52:01.780177Z",
     "shell.execute_reply": "2021-02-02T06:52:01.781219Z"
    }
   },
   "outputs": [],
   "source": [
    "nsample = 100\n",
    "x = np.linspace(0, 10, 100)\n",
    "X = np.column_stack((x, x**2))\n",
    "beta = np.array([1, 0.1, 10])\n",
    "e = np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model needs an intercept so we add a column of 1s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:01.786230Z",
     "iopub.status.busy": "2021-02-02T06:52:01.784746Z",
     "iopub.status.idle": "2021-02-02T06:52:01.791358Z",
     "shell.execute_reply": "2021-02-02T06:52:01.792414Z"
    }
   },
   "outputs": [],
   "source": [
    "X = sm.add_constant(X)\n",
    "y = np.dot(X, beta) + e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:01.797272Z",
     "iopub.status.busy": "2021-02-02T06:52:01.795610Z",
     "iopub.status.idle": "2021-02-02T06:52:01.832720Z",
     "shell.execute_reply": "2021-02-02T06:52:01.833669Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       1.000\n",
      "Model:                            OLS   Adj. R-squared:                  1.000\n",
      "Method:                 Least Squares   F-statistic:                 4.020e+06\n",
      "Date:                Tue, 02 Feb 2021   Prob (F-statistic):          2.83e-239\n",
      "Time:                        06:52:01   Log-Likelihood:                -146.51\n",
      "No. Observations:                 100   AIC:                             299.0\n",
      "Df Residuals:                      97   BIC:                             306.8\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.3423      0.313      4.292      0.000       0.722       1.963\n",
      "x1            -0.0402      0.145     -0.278      0.781      -0.327       0.247\n",
      "x2            10.0103      0.014    715.745      0.000       9.982      10.038\n",
      "==============================================================================\n",
      "Omnibus:                        2.042   Durbin-Watson:                   2.274\n",
      "Prob(Omnibus):                  0.360   Jarque-Bera (JB):                1.875\n",
      "Skew:                           0.234   Prob(JB):                        0.392\n",
      "Kurtosis:                       2.519   Cond. No.                         144.\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "model = sm.OLS(y, X)\n",
    "results = model.fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quantities of interest can be extracted directly from the fitted model. Type ``dir(results)`` for a full list. Here are some examples:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:01.838974Z",
     "iopub.status.busy": "2021-02-02T06:52:01.837389Z",
     "iopub.status.idle": "2021-02-02T06:52:01.846526Z",
     "shell.execute_reply": "2021-02-02T06:52:01.847692Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [ 1.34233516 -0.04024948 10.01025357]\n",
      "R2:  0.9999879365025871\n"
     ]
    }
   ],
   "source": [
    "print('Parameters: ', results.params)\n",
    "print('R2: ', results.rsquared)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS non-linear curve but linear in parameters\n",
    "\n",
    "We simulate artificial data with a non-linear relationship between x and y:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:01.852389Z",
     "iopub.status.busy": "2021-02-02T06:52:01.851005Z",
     "iopub.status.idle": "2021-02-02T06:52:01.859463Z",
     "shell.execute_reply": "2021-02-02T06:52:01.860657Z"
    }
   },
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "sig = 0.5\n",
    "x = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x, np.sin(x), (x-5)**2, np.ones(nsample)))\n",
    "beta = [0.5, 0.5, -0.02, 5.]\n",
    "\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + sig * np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:01.865293Z",
     "iopub.status.busy": "2021-02-02T06:52:01.863952Z",
     "iopub.status.idle": "2021-02-02T06:52:01.884278Z",
     "shell.execute_reply": "2021-02-02T06:52:01.885187Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.933\n",
      "Model:                            OLS   Adj. R-squared:                  0.928\n",
      "Method:                 Least Squares   F-statistic:                     211.8\n",
      "Date:                Tue, 02 Feb 2021   Prob (F-statistic):           6.30e-27\n",
      "Time:                        06:52:01   Log-Likelihood:                -34.438\n",
      "No. Observations:                  50   AIC:                             76.88\n",
      "Df Residuals:                      46   BIC:                             84.52\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             0.4687      0.026     17.751      0.000       0.416       0.522\n",
      "x2             0.4836      0.104      4.659      0.000       0.275       0.693\n",
      "x3            -0.0174      0.002     -7.507      0.000      -0.022      -0.013\n",
      "const          5.2058      0.171     30.405      0.000       4.861       5.550\n",
      "==============================================================================\n",
      "Omnibus:                        0.655   Durbin-Watson:                   2.896\n",
      "Prob(Omnibus):                  0.721   Jarque-Bera (JB):                0.360\n",
      "Skew:                           0.207   Prob(JB):                        0.835\n",
      "Kurtosis:                       3.026   Cond. No.                         221.\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "res = sm.OLS(y, X).fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract other quantities of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:01.890264Z",
     "iopub.status.busy": "2021-02-02T06:52:01.888718Z",
     "iopub.status.idle": "2021-02-02T06:52:01.900083Z",
     "shell.execute_reply": "2021-02-02T06:52:01.901109Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [ 0.46872448  0.48360119 -0.01740479  5.20584496]\n",
      "Standard errors:  [0.02640602 0.10380518 0.00231847 0.17121765]\n",
      "Predicted values:  [ 4.77072516  5.22213464  5.63620761  5.98658823  6.25643234  6.44117491\n",
      "  6.54928009  6.60085051  6.62432454  6.6518039   6.71377946  6.83412169\n",
      "  7.02615877  7.29048685  7.61487206  7.97626054  8.34456611  8.68761335\n",
      "  8.97642389  9.18997755  9.31866582  9.36587056  9.34740836  9.28893189\n",
      "  9.22171529  9.17751587  9.1833565   9.25708583  9.40444579  9.61812821\n",
      "  9.87897556 10.15912843 10.42660281 10.65054491 10.8063004  10.87946503\n",
      " 10.86825119 10.78378163 10.64826203 10.49133265 10.34519853 10.23933827\n",
      " 10.19566084 10.22490593 10.32487947 10.48081414 10.66779556 10.85485568\n",
      " 11.01006072 11.10575781]\n"
     ]
    }
   ],
   "source": [
    "print('Parameters: ', res.params)\n",
    "print('Standard errors: ', res.bse)\n",
    "print('Predicted values: ', res.predict())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare the true relationship to OLS predictions. Confidence intervals around the predictions are built using the ``wls_prediction_std`` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:01.906299Z",
     "iopub.status.busy": "2021-02-02T06:52:01.904630Z",
     "iopub.status.idle": "2021-02-02T06:52:02.290660Z",
     "shell.execute_reply": "2021-02-02T06:52:02.291839Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5830a726d0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFlCAYAAAAzqTv+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABksklEQVR4nO3ddXxV9RvA8c/ZXdIpwgApJQQBAWUiMkFCSTHBxE7EQECQEP0BgoUoFqIIIoI0Ij1CRoc0EqKMjlHbWH1/fzwbDNjG4tzc83699mK7uzvne3bZfc63nscyxqCUUkop1/JzdwOUUkqpvEgDsFJKKeUGGoCVUkopN9AArJRSSrmBBmCllFLKDTQAK6WUUm7g78qTlShRwlSoUMGVp1RKKaXcZu3atceMMSXT+55LA3CFChVYs2aNK0+plFJKuY1lWfsy+p4OQSullFJuoAFYKaWUcgMNwEoppZQbuHQOOD0JCQns37+fuLg4dzfFaYKDgylbtiwBAQHubopSSikP4fYAvH//fgoWLEiFChWwLMvdzbGdMYbjx4+zf/9+Klas6O7mKKWU8hBuH4KOi4ujePHiPhl8ASzLonjx4j7dw1dKKZV9bg/AgM8G31S+fn1KKaWyzyMCsCfp378/w4YNy/D7U6dOZevWrS5skVJKKV/kdQF46vooGg1eSMWes2g0eCFT10e59vwagJVSStnAqwLw1PVR9Jq8iajoWAwQFR1Lr8mbch2EP/jgA6pWrcpdd93Fjh07APj2229p0KABtWvX5r777iMmJobly5czffp0unfvTp06ddi9e3e6z1NKKaWuxqsC8NA5O4hNSLrksdiEJIbO2ZHjY65du5ZffvmF9evXM3nyZFavXg1Ax44dWb16NRs3bqR69eqMGjWK2267jXbt2jF06FA2bNhA5cqV032eUkopdTVu34aUHQeiY7P1eFYsXbqUe++9l3z58gHQrl07ADZv3kyfPn2Ijo7m7NmztGzZMt2fz+rzlFJKebhTp6BwYZedzqt6wGWKhGTr8axKb5Xyk08+yYgRI9i0aRP9+vXLcBtRVp+nlFLKw/XpA1GuW1fkVQG4e8uqhAQ4LnksJMBB95ZVc3zMO+64gylTphAbG8uZM2eYMWMGAGfOnKF06dIkJCQwbty4C88vWLAgZ86cufB1Rs9TSinl4dasgfvug5Ur5euePcGFGQu9KgB3qBvKoI61CC0SggWEFglhUMdadKgbmuNj3nzzzTz00EPUqVOH++67j8aNGwMwcOBAbr31Vpo3b061atUuPP/hhx9m6NCh1K1bl927d2f4PKWUUh7IGIiIgBYtoEEDWLAA9uyR74WGwjXXuKwpljHGZSerX7++ubwe8LZt26hevbrL2uAueeU6lVLKo7VpA7NmQalS8MYb8MILUKiQ005nWdZaY0z99L531UVYlmV9D7QBjhhjaqY8NhRoC8QDu4Euxpho21qslFJK2eX8eQgKks9btYJ77oEuXSAkd+uHcisrQ9A/AK0ue2weUNMYcxOwE+hlc7uUUkqp3Fu3Dm66CcaOla9feQVeesntwReyEICNMUuAE5c9NtcYk5jy5QqgrBPappRSSuVMcjIMGwYNG8K5c1DW88KUHfuAnwIm2HAcpZRSKvcOHIAnnoD586FjR/jmGyhe3N2tukKuVkFbltUbSAQy3H9jWdZzlmWtsSxrzdGjR3NzOqWUUurqVq6E5csl8E6a5JHBF3LRA7Ys6wlkcVYzk8lSamPMN8A3IKugc3o+pZSy09T1UQyds4MD0bGUKRJC95ZVs76l0RjYvBmmTYN//oHHHoMmTSST0oIF8oZfrJj8W7z4xQVAynni4uDPP6FZM7j3Xti9G6691t2tylSOArBlWa2AHkATY4xXVx84fvw4zZo1A+DQoUM4HA5KliwJwKpVqwgMDHRn85RSTpBa2CU1t3xqYRcg8yB85gz06wdTp8LevfJYqVISfAG2b5fEDmkFBsK4cXD//TZfhbrg4EEJuhs2yOtSurTHB1/I2jak8UA4UMKyrP1AP2TVcxAwLyWN4wpjzAtObKfTFC9enA0bNgBSC7hAgQK89dZbF76fmJiIv79XpcxWSl1FZoVdLgnA587BnDkQEwOPPgr58sFvv0HNmpI1qW1bebNPVasWrF8PJ07A8ePy7+bNshAIYO5cyb705JNQpozzLzQvWLsW2reHkyfh558vfT083FUjizGmUzoP+3TJnyeffJJixYqxfv16br75ZgoWLHhJYK5ZsyYzZ86kQoUKjB07luHDhxMfH8+tt97Kl19+icPhuMoZlFLulKXCLhMnwosvSiCtU0cCsMMhQ5sZ3ZTnyyfPzcjChTBkCPTtC61bwzPPwN13Z3w8lblff5WbmRIlZPg5s9+9B/KoV71bNxlBsFOdOvDpp9n/uZ07dzJ//nwcDgf9+/dP9znbtm1jwoQJ/PnnnwQEBPDSSy8xbtw4Hn/88dw0WSnlZGWKhBCVThAuUyREeq2vvALjx0uqwokTISVFLZC7YDl4MDz9NIwaBT/8ANOnQ3g4LFqU82PmZX/9BXXrwuTJMhXgZTwqAHuSBx544Ko92QULFrB27VoaNGgAQGxsLNe4MI+oUipnureseskcMKQp7LJzJ0yZAgMHyjCz3b3T66+XQDxwIMycKQu6AOLjpQdyyy32ns/XnDsn87w1a8J770FCgtcucvOoAJyTnqqz5M+f/8Ln/v7+JCcnX/g6teSgMYYnnniCQYMGubx9SqmcS53nTV0FXSXE8L+gvTSo2woIlZXNzu5RBQTIwqFUX38NXbvK/tUhQ7yyR+d0+/bJfO+RI7Brlwz5e2nwBS+rhuQuFSpUYN26dQCsW7eOvSmrH5s1a8akSZM4cuQIACdOnGDfvn1ua6dSKus61A3lz55N2dsqP/N+eJUGfV69dGWzq3XpIj3un3+GG26Azz6DxMSr/1xesWyZjA7s3Qvffy/B18tpAM6C++67jxMnTlCnTh1GjhzJDTfcAECNGjV4//33adGiBTfddBPNmzfn4MGDbm6tUipL4uPh9dfhzjulN7p0KVSs6L72FCgAgwbBpk2yarpbN+kN53XGwJdfyutUuLAk2Wh1eXkC76TlCF0kr1ynUl7BGNmXO3myLLgaPBjSTDu5nTGy17hMGbj1Vjh9WuY6PTSjk1MZI0P1iYlSUKFIEXe3KFtyVY5QKaV8jmXJm3qzZlIZx9Okti9Vjx6yMOzLLyW3cV7w338SfMuXl2H54GDw861BW9+6GqWUykxMjMwlguzr9cTgm54XXpDe8H33wYMPyiIkX7Z4MdSrd3EIPl8+nwu+oAFYKZVXnD0ryS9atIBDh9zdmuypXVvmPj/4QPJP16ghOad9jTGyHaZZMxluHznStkNPXR9Fo8ELqdhzFo0GL2Tq+ijbjp1TGoCVUr7v9GlZuLNkCXz3nVfkCb5CQAC8846kuqxTBypVcneL7BUTI0UtXn9dUnyuXAnVqtly6NTc31HRsRgu5v52dxDWAKyU8m3R0dCypbyh//ILdO7s7hblTo0aUue2YkXpMXbuDKNHSwF6b5acLJmtBg6UfNuFCtl26Mxyf7uTBmCllG8bNUoS9k+cCA884O7W2Ov0adi/H556Cho1gtWr3d2i7ImJgQEDpMpUgQLS/j59bJ/vzVLubzfQAAzs37+f9u3bc/3111O5cmVee+014uPjiYiIoE2bNlc8f+bMmdStW5fatWtTo0YNvv76aze0WimVJW+8IW/sHTq4uyX2K1wYIiIkr/Q//0iiii5dJJ+1p1u4UKpH9e8Pv/8ujzkpq1WZIiHZetxV8nwANsbQsWNHOnTowN9//83OnTs5e/YsvXv3Tvf5CQkJPPfcc8yYMYONGzeyfv16wsPDXdtopVTmYmKgUyfYsUO29NSu7e4WOY+fn6wW3rED3n5bAltAgLtblbHoaKkC1ayZtD0iAh56yKmn7N6yKiEBl+b2v5D72428MwBHRkrGmMjIXB9q4cKFBAcH06VLFwAcDgeffPIJ33//PTExMVc8/8yZMyQmJlI8ZUN8UFAQVau690VUSqWRlASPPAITJsD27e5ujesUKiQ5pHfsgIIFJXFH+/YXe5ee4sUXpcfeo4fM+TZp4vRTdqgbyqCOtQgtEoIFhBYJYVDHWpfWfnYDz0vEkV5v8sEHZb9eTIzMc/z1l0zY+/nBTTfBa69JTchjxyS7TVoREZmebsuWLdSrV++SxwoVKkT58uXZtWvXFc8vVqwY7dq147rrrqNZs2a0adOGTp064eeDe9SU8jrGyCraqVMll3L79u5ukesFB8u/+/fLDUjr1nDPPRLwGjeWEQFXW71abhCqVpWtVG+9Jft8XahD3VC3B9zLeV/UOHXq4mq/5GT5OheMMVjp/IfM6HGA7777jgULFnDLLbcwbNgwnnrqqVy1QSllk08/hc8/lyDctau7W+NeFStKXulhw2D5culpVq8OBw645vyxsbI6u0EDmZv+8EN5vFIllwdfj2WMcdlHvXr1zOW2bt16xWOZWr7cmJAQYxwO+Xf58uz9/GXmzZtnGjdufMljp06dMsWKFTOzZs0yrVu3zvTnjx49agoUKHDV82T7OpVS2ZOQYMzttxtz333GJCW5uzWe5dw5Y374wZhOnYxJTpbHfvrJmMWLL35tp/feM6ZoUWPAmBo1jBkxwphTp+w/jxcA1pgMYqL39YDDwiQDzMCB8m9YWK4O16xZM2JiYhgzZgwASUlJvPnmmzz55JPkS6fc1dmzZ4lIM6y9YcMGrrvuuly1QSllA39/mDsXfvrJJ9MW5kq+fLJQ6+efZQg6ORn69r3YK/74Y9i8WXqt2WWMTP/Nnn1xdPLcObjrLli0SI778su27uv1FVoNCfjvv/946aWX2L59O8nJydxzzz0MGzaMyMhI7r777gsLrgDGjx/PoEGD2L17NyEhIeTPn5/PPvuM+vXTLXZxgSdcp1I+adcu6N0bvvlGtuWorImJkb3RX399cUFr376yLzc6WrJu3XADXH+9/BscLJWIChaUfdVDhsDu3fKROhU4f76sbjbGPXPNHkirIV1FuXLlmDFjxhWPh4eHE5vOHWHjxo1d0Syl1NUcOwZ33w0nT8LRoxqAsyO1V/zEE7JYa8MGybIFsoBr/HgJxGlNmCCLYuPi5PmVK8soZOXK8rOpK5q9MPhOXR/F0Dk7OBAdS5kiIXRvWdXpi7Y0ACulvFNsLLRrJ2XrFi6EKlXc3SLvVa3apXmXa9aUZB7Hj8POnfJx/rwsqALZjbJzp3va6gSpuaJT01Wm5ooGnBqENQArpbxPcjI8/jisWCHDqLfd5u4W+R7LghIl5MPHf7+Z5YrWAKyUUmkdOCDBd+hQqZHrBu4YslTO4a5c0R4RgE0me259gSsXuimVJ5QtK3tc3TTn664hS+UcZYqEEBUdy81R22j47yZWlK/FutDqTs8V7fa1+sHBwRw/ftxng5QxhuPHjxOcmp1GKZVzCxdKcYWkJFmR66Ybd08tb6dypnvLqtwRtZmJ43rwxtKxjPulN2GHdzo9V7Tbe8Bly5Zl//79HD161N1NcZrg4GDKli3r7mYo5d22b5fh5tBQqaDjxn2lnlreTmXT7t2wciUdOnfm+iLRWCYZB0BSIn0KHOFGX18FHRAQQMWKFd3dDKWUJzt6VPIZBwbCzJluT+qQOmSZ3uPKg0VGSgIny5I6AfPnQ0gItG3LjZ3bwajPID4e/8BA+drJ3B6AlVIqU3FxUsv34EF506xQwc0NkiHLtHPA4Bnl7VQmIiOl2E98vHxdqhS89x489ZQkF0nNshgRIc/LZZbFrNAArJTybOvXS9KHn36CW291d2uAiwutdBW0l0hIkMCalHLD5OcHr74qGdTSCgtzSeBNpQFYKeXZwsJgzx7psXgQTyxvpy5z+jS8/Tb8/bfUDwgMlB5wYCA0bXrhacnJ0kGeMEE+HzHCNc3TAKyU8kxjx0r2paef9rjgm+rECRkZv+YaKFYMHA53t0hd8Pvv8Pzzsme8WzeoX/+SIWbTMIzVqyToTpwoCdWCgqSkvKtSWWsAVkp5nrlzZW7u9tvhySc9JrIlJUlt+T/+gDlzwG/FcpqwiEU0ZZVfGCVLSjC+5hq5Z2iQGEn9sxFUfT6cku1cN7SZp0VHy/Dy2LFw440wadKFqQvTMIwNwWFMmAC/PgJ790JAALRqBYMGSWbTggVd11QNwEopz7J8Odx7ryT3/+03twffgwcl2P7xB8ybB2dPnOdOIvhf0a8JZypYFon+HzDnziHc9udQzsUU4MyeAiQlJFEjfgMGiP89iF5NFtC0dxjNmmm1RKeIjJTebYMGkiWtXz+p6BQYCMiC5zfekPwtDodUS3z3XVnfV7Soe5qsAVgp5Tk2bpTtRqGhEvXc9c6Y0pSRj0dS9K8IIgjndMkq/F7oFW6OnU1A7Bk4GwAYMIaA5Hja3LATyjSj2NmzcPYs5zZuxu9gMhZgiKfO8hHsbTGax0IfoW7Xxjz5lB8lSrjt8nxLRAS0aCETuIGBUps4pTLT/v3w5pvw669QqZJUX+zYEY/43WsAVkp5jsWLZY/v/Plum/eNi5P1OksHL2N+clP8SYSgYJg0B78XN8OjD8tYZYECcrOQuqinc+cLK2inro9iwvBf+X5sTwKSEkl0OIivdYYnNk3luahv+bdHOX54pxMnWj3Cw23PUet4BNad4VeswM0r+aZzdZ379sGjj8pKZ5DXY/ly4sOa8OmnstMoKUn+7d5dyhp7CsuVKSDr169v1qxZ47LzKaW8RNpVL6dOuS3H85Il8OyzYO3czrL8LSlx7l/5hsMhUblXr0t/IHXY87J9o40GL0w3t3DlfLCgyinOfDWOfMvm8A/XUdocJIh4rOBA/BYuuCSIp7fXeFDHWrYGYXcH+Vxd57x50KmTlKZMTJRIGxjI6iELePzLMLZvl3ulTz8Fd+V7sixrrTGmfnrf05kIpZR7HTki5e4iI+VrNwTfU6fghRegWZMEnjk6iC0BdSjhd1J6tg6H/BsefuUPhoVJUL6s55qaknJdaHW+DHuQdaHVAdgTA3TuTMEls3AcPki55+4h2IrHQRLExbLxjR9JTJRjuCLfdGrwi4qOxXCxqMTU9VG2neNqcnydo0bJ6qnSpWWfeEQEp94aSJ+wBdzSNYz4eJgxA6ZNc1/wvRodglZKuU90NLRsCTt2XEySkAV29tqmToWXX4ZDh+DTh1fz6i/vwP33M/v53kz5bSlVtq5hV4363BNcng5ZPGaWUlWWLEngE51hzCjM+Xis5CRqr/iaOddEU2bMENvyTWf2u3JXHdy0cnydjRrJCvnhwzH58vPjn9fzyogwkpJgwADZ/utJw83p0QCslHKPc+egTRvYsgWmT5ctR1lgVynA06dli/H0Sed5qmIET61oSYMGt8Hb65jKNXKOwpWYG1YJgKXZOEeWU1WmpD+0IiIwt9zK1q8Wc8ekoVhtp9Kt7DN8et/dmMuCSHbyTV/td+UJRSWynFc7NVPGyZPwww9QrRqMGsXZs/DSE5IoLTxcOsaVKrmk6bmmQ9BKKdeLj5fKRpGR8PPPMpSYRXYMzR47JomQik3+luP5y/PlvntoUGKvfLNu3Vyfo0PdUAZ1rEVokRAsILRISMZzminD2FazptSYOID4TTvYeP39XL9/D/t/aMqNqw7xUuSv3By1Ldv5pq92HRkFc1cWlejesiohAZduNbviOlPzOH/2GYwZI2PLyEr1evVg3Djp9c6f7z3BF7QHrJRyB4cD8uWDb7+V1EPZkNteW1SU7Fh5fMc7vJ08COscMsd76NCFyUI7eoY5TVVZuGY5bt05luWLE2jdeQXjF71CMHHE+wey4tuJNMnGMa92HZ5QVOKqebUTE2U8ObWIgsOB2byFr6La8frrkoFswYL0p+g9nQZgpZTr7Nkjwa5sWclQlIOMFLkpBbhnjyRgeCjqY3okDbr4jaQkWc2cspjKE8oN3tYkgEkvLsPv3XgsIDAxnjojP4dH7pb0TVlwtevwlKISmd6sPPAALFsG/v5gDCYwkD7zw/nfIhk4GTMGSpZ0aXNto0PQSinXWLkSGjaExx+Xr3OYDipLQ5bp2LJFpplPnYIXXnbIGHRISLqrnHN6Drv5NwvHLyQI4+cgGQeFV83ncIVbSF63IUs/n5Xr6FA3lD97NmXv4Nb82bOp5+0zfvZZmeBdsoT9zw/k/iILGLIkjCFDYNYs7w2+oPuAlVKu8NtvkiyhTBlJkl81d4Esu6ug16yBu1smU83xNyMXVqVmTWTv8YoVGdZ/dff+2AtS9hqfvjmcrwcc4rHIF1ldtiMNVn/Jtdde/cc95jqyyhgYOVKGnLt1u/DQp59Cjx6y6+iXX1xaNTBXMtsHrAFYKeU8xsDHH0sKoltvldXOLu6yLFkC7Vsn8q15ho78ht+WzXDddS5tg12MgZ8+O8FbvQKgYEGm9F1Po5B1spfaRUXknSomBl58UcaV27eHKVM4fsKi1b2xrFkaQr7rD1Gz0056dajs2TcRaWQWgHUOWCkP4XU9layIi5MtIx07yjBiiOvmUEE62506nmdSYCean50i+QjLl3dpG+xkWfB4t2LUbwEPPwyOV1/EsBIsCys4WFYjeWsQnjRJqhgdOiRLmvv04c/lFh3uT+T40SCK3rWFgjf/w+Hz5GjbmSfSOWClPIAnZCSy1dmzkh4wJESGeH/91eXBd8mQSFa3GcCfjsY0PzNFtrC8+65rCr06WY0aMqV+on5LDBaWMZjYWBm69Ubz5sliq0OHIDCQ5GbNGfyhH02awJn4BEo9+ieF6v1z4aWzOyOYu2gAVsoDuCLtoMvMmAE33SRZigCKF3d5/b2/vomkfs9mvGsGcGPMaujdG7p2dWkbnC0kBO4Z3orkwGAS8cNgwU8/kTTxN3c3LeuOHZN/16y5cGNkkpL4sUsEvXrJwEmpx5cQdO3pK3407RarqeujaDR4IRV7zqLR4IVec+OqAVgpZ4qJgbVrYetW+frYMRkivP9+eP11+Ogj+PVX/PfuSffHXZmRKNf27IG2bSX7fUiI5Hd0gx07YFq3CAKJl7Dk5wf587ulLU4XFoZ/xAJi33mf3o0W8zxf0Xx4O/bsQf7PpVYI8jRJSfDJJzIXP2+ezF8HB5Ps5yAuKZAf/gln5EhJfFW2VPpbrlK3Unnz6JHOAStlt/79JTn85s0SlIyR3uDo0dIbLFpU9sT88YekYwQ6t3iKQUU7cs2Z40we251V5W7kz+vqsPumW915JVk3dapUpfH3h2HDpLeZxb2qdjp8GFq3SuJ/SVvxCwqARDIupOArwsIoGBbG/wyMGdOYX7pCw1rn+Me/GSFli2GNHCmvRQarvV1u61bJAbpihaQirVGDuOKh/PTwAvb+EMHucuEMnxFG7dry9KslC/GEfNY5pauglbJDVJQUkQeoVUsKg994I9SsKR91615ZksUYKUawfz9/HDjP68uOUejEYd5dOIqG//5FiZhT8ryqVSVjVOPGLr2kLDl3TnqXBw7AO+/ABx9c/D242NmzEN7E8Mxfr/JC4heyZ6VwYc8IOi7077/w1FMQsmAG34W8SqnYfbLX2RgICnLfQq3ISHj/fZg7V16X4cOhUyfmzbd4+WX4+2947DH44gsoWPDSH81sgWLFnrNIL4pZwN7BrZ1+WVejq6CVcpatWyXwzJsn7yBlykiC2qzMeVqW9IaLFqVVLYi7JoqhcwJ5tX0PQgsF8V5lQ9OoTfKGmbrhc+xY+PxzaN5cqgg1bOiWniZ79sAbb8CJE7B4sVz3Dz+4vh0pEhPhwQeh+fohvGC+gLfegsGD3dYedypfXmLcl1+25cbuTZni35rbExdjgeytTZPxyyXOn5fge889sirezw9+/JEDdVvzRicZZr7+emlz8+bpHyKzTFmekLUsx4wxLvuoV6+eUcon/PuvMV26GOPnZ0yhQsa8/74xZ886/7yTJxtz223GOBzGgDEFCxrTvr0xsbHOP/e+fcZ8+KExDRvKufPlM2bwYGMSEpx/7kwkJxvzzDPGPMaP0q7OnY1JSnJrmzzFjh3GPHPjchNDsEnCzyQGhRizfLkxH39szC+/OPe127fPmF69jClZ0ph27S78n012OMyy1v8zBQsaExRkzIABufvvO2XdflOtz2xzXY+ZFz6q9Zltpqzbb9+15AKwxmQQE3UIWqnsOnoUKlSQbtcrr0hB9hIlXNuG6GhYtAjmzIF//pH5ZIDXXpNezu23Q/XqMnyd0wVIxkhvvlw5mbv+/nuZu6tXT5anPv645HS+jKv3Mw8cCAP6JvJvqQaUqVlcNv8GBjrtfJ4kK7/rxESY3iuSHd9EMP10OPnCb2Xa/pspsCvlte3aVVatr12b8+H6lGxdhIfLwsMRIyTpCsg8b+vW0K0bJj6e88mB3GkWULhlGCNGQJUquf0tePYees2EpZQdjh69mMXpu+9kvMzTMio9/jhMmSIToqk6dZKSfyDjfWXKSOCMiYEzZ+SaKleW4cFvvpHHDh2SRLt798qb6csvSwHdkyczvebL68+CLJjJsBRfLv3wA3TpIpf9wycnsRx+Mr+YB2T3dx0XB19/DYMGwZHDybxbdxZvWh9TaF2EPMHPT+aIZ8yQ17hSpYtTKWkDbNoAbQzMni03ZImJcuNTuTIcPCg5nF94Aa67jr/+gll9IjkzI4LNJcJ5fGQY993nE1uyr0oDsFK5NW0aPPKIZOvJRu1at4iPh507Yft22LZNejlPPilbUvLnv3JrSrdusiUkJuZibzkkRN5sO3aUlIBZTB/ZaPDCdOfjQouE8GfPprm6rMvNnQsv37OXT64dQoutnxJYKPjqP+RDcvq7PncOvvwShgyB48dhRoVXaf3PF1gYWazVubNkLStQQBYUliolN2PJyRJg+/aVef99++QjJubiwR0O2V43cCBHzwTz889yk7RhgyxVeOUVSXJ1+SIrX6aLsJTKjdGj4ZlnoEED+fB0gYEXV1+n5e8vC8W2bZMeSv78UKiQ9FhAgu6xY/LumMMhXDvq6GbF5s3wVfvZrOJxipw6j3X4TSh0va3ncIXcDJ3m9HedP7+k5n7hBVmIPHxwZ5oyikDiSbYC2Vj+Pm787A7y7foLNm4k6ffZOFJu2hLjznN4/lJCTx6RKY5WraQX/OWXkJSECQxkcbGOfPJQML//Lp3ievUuLHh2+UyNp9MArFRmhg6VYuAtWkhFnwIF3N2inLMsGVrMaAjZsmSuNxdcsSL1yBEYeFcEE+Pa4EcyVkKQ3Dhc710B+PIh5NQEEpC1HMe5/V0XLCgJwl5+OYyJPRcQ90cE46LCWfpBGJYFtWtD2epnsW6byK9LX8A/OZEEhz+vl2zFna8+QL3S13LypCyED3Q8SMCfEXy5NZw574Rx7bUysPLEE1feB6qLdAhaqYwsWJBSvf0hqc6SRxb25Iaz54DPn4dmTQ1fRtbhJvOXPOhwyEqsXr1yfXxXyu1wvTN+17Gxkh9jyRIZZV68NInkRAcNieROxwIWJTVjBekv0goKktmKJ56Q+1X/NN07T14k5Ww6BK1UTjRtChMnwr33ypu8uqrUN1VnvNkaI+t6DizfS7WQvZAYcHFe0gszXeV2uN4Zv+uQELjzTvkAqPDWXOIOFWb7v8XZEvMIfkEJFA3ZgiM4gdEv1KFoUShWjAv/pnePmtuevi+7agC2LOt7oA1wxBhTM+WxYsAEoALwD/CgMeak85qplIvExcFLL8Gbb0omq/vvd3eLvE5mSRNyY/BgWRs0YEAlAp/bKclAFi/22kxXdgzXO+t3nSq0RBBR/icJLnvp23tokRDatMnaMbw5VaSzZaUYww/A5cs+ewILjDHXAwtSvlbKu505I9l6Ro+WcTjlMSZPht/fWcovN/2Pd/sYyQx2220y7OyFwRckx3FIwKUjK2lzHHsCO9roqoV53uiqAdgYswQ4cdnD7YEfUz7/Eehgb7OUcrHkZFmmuWSJpHt8+ml3t0ilWLcO+j2yixn+9/JA3Bisc2ev/kNeoEPdUAZ1rEVokRAspFfprP3SOWVHGzPq0XtcqsgTJ2ShpQtlaRGWZVkVgJlphqCjjTFF0nz/pDGmaAY/+xzwHED58uXr7du3z4ZmK2Wzvn1lIc8XX8gQtPIIBw7AXfVOMuNYQyoWOo7fyhX2pE5SLuPq5CzZcvSoTGD7+0v+8I8+kj30Nq6oz2wRltPrARtjvjHG1DfG1C+Zxc38SrlUUhKsWiUlZF580d2tUSliYqBjm3hGHr2PitY/+E2dkuXg660F2n2Rx/X0Dx6UfcvNmslUxpIl8vjLL8Pq1S69wcvpKujDlmWVNsYctCyrNHDEzkYp5VIOh2T6SUzMG7nxvEByMnzQJpKn139PY78l+H3/Q5bLMeqqW8/j7MViWXLggOyRWrBAltRXrSqVzCpVku9XrHhlyVAny2kPeDrwRMrnTwDT7GmOUi4UHS1JhA8elCAcFOTuFinkvXF4p0h6L2rG09Zo/AIDLmbryoLMVt2qPCg5Wf4tXlw2OvftK6nUtm2TaacKFdzWtKsGYMuyxgORQFXLsvZblvU0MBhoblnW30DzlK+V8h7JyfDoozB+POze7e7WqDSGDoVrf/2UYM7jZ5Ikd3VERJZ/XlfdKkBGtD7/HOrUkQTYQUGwdCn07y9bDD1gtOuqQ9DGmE4ZfKuZzW1RynX695dh5xEjpHSfyjJnZjUaMwYW9ZjNDH7D8gMsR7YTbbisQPv58xdHTdasgfXrpeJAQIC0OSBAyvAFBNh7XnV1CxZIac4tW6RqWXS0JMH2gKCblmbCUnnPlCky9NSli654ziZnzq/+8Qd8+1Qk8/zuw++mm7CGfiiLYrKZaKN7y6rprrq1ZX/tkSNy4zZ9upRj2rpVcmvPmAHvvXfl80+flgA8bJiUBOrQQQoYeHNOcU927pzM8/72m8znTp0K7dp5XOBNpbmgVd6SnCwVjfz9JYtSsH0l7PJCvltnlRtctQpearKFBQmNKXBdcRzLl0kZvByy/bX46y9ZIR8ZKZPU5crJG3v37hKAz56FU6ekFGRCgnzEx0tFAz8/Cc7Dh0v9v+Bg6ZV16iQfyj7GQNu2kqTljTds/fvOKa0HrFRa0dGyGKN0adsO6dF7HdPIbWCq2HMW6b1jWMDewa1z1KadO6FRI3gn8T26Bo7EsWK5y1ejpuvPP+UN/fbb4fBhGU5u21YqDtSunf1eVWIiLFsmIzBTp8LNN8vnAGvXytce2lPzaNHRkjq2f3+5MTLGo36Pbt0HrJRHMEZSTMbFQZEitgZf8I6Vt6k3CVHRsRguDh9nZ4+s3VmNDh6Eli3l/bLtqndx/LXB/cH35El47jkJvP36yWOlSsk8b79+sqgnJ2/w/v4ynP7ZZ/DPP1KpHmDTJhmVadxYgr6Xceue69Wr5cZlzJiLvzsPCr5XowFY5Q1jx0qijTFjnHJ4b1h5a8dNgp35i0+dgo4tzvLZfx1Z+PkWqlxv5WrYOdeMkVXx1arB999LZqTp051zLsuCwoXl8+rV4auvpLjE7bfLPPHWrc45r83suKnLEWPgk09k6CQpSZJpPPywc8/pBBqAle/bvx9efVXe3JyU49kb8t3acZNgV1ajuDi4v108723pSBsznZr592br551ixgzo3Fn2ha5ZI/uh8ud3/nn9/aXH/fff8MEHsGgRNGwoxUE8nNtGfj75ROZ477lHVp97aUEOXQWtfJsxEnQTE2XIz0l1fZ268tYmdm3PyW1Wo1OnoHf4nwzb8CK12QTfjybLte3sFh8vvc06daQN48fDAw+4p/5z/vySmem552RotWBB+f/71VdyY5DaY/YgLh/5SUyUG5ZnnoFCheRv24uGnC+nPWDl277+WraLDBuWrWxK2eVx+W7T4Qnl7w4dgjfqLebTDU0k+AYESEpAd9i5U+YPmzaVuwI/PxnGdEfwTatECbj7bvl83TrZKle9uvOGw3PBZSM/xsiIRMOGsoCyUCEJwl4cfEEDsPJ1TZrIUNXzzzv9VB3qhvJnz6bsHdyaP3s29ajgC+6/Sdi1S6bsrtu3BAcpIwXJydnKcmWb1aulMYcPy7oAD+xdAlCvnrS1ZElZfd25Mxw75u5WXeCSm7pTp6BjR3j7bVmgl5ho37HdTLchKd/kYVsR8rp16+CJFgeJTQ5i2oc7uLFrMxn+DQyUrEWunMObMwfuuw+uuUY+t7H0nNPEx8OgQfD++5JGcf16j/n/7dT973/9Ja/VP/9ID/i11zzmurNK9wGrvOejj6Tn8MMPHrEZPy9bsAC6t9vBtPMtKRZWjfxL/5CEFhER2c5yZYunnpIANnu2lKPzJn/9Jck87rxTkn2cPCk3Er7IGEmosW8f/Pqr16aM1QCs8patW2Vu7+67YfLkLN0x54UsVu7w66/w+SMrmG7aUKioA8cfv8uwqjucOiVDzQkJF+cRvdn//ic3msOHy9C0l/UMMxQXJ8PMBQrA3r2QL597t6flkibiUHlHQoLkgi1YUFaPZjH4umUvo4/74gv46aGZzEtuSqHyRSTDlTuCb3KyZEqqXx9OnJCFX94efAHuvRduuEGqerVvL9vtvN0//0hCkmeeka8rVvTq4Hs1GoCVbxk8WPZwjhyZ5T9cb8hi5U3OnIGP7o/k4Cvv833wCwTWuVGCrxNXoWcoPl5qPn/8sYyIFCni+jY4S/Xqktpy2DCYP1/mhn/7zd2tyrnfU0ZH/v47z+TI1n3AynecOwdffinDcfffn+Ufc+ZexuRk2LFDpqNXr4aYBZFUPRjB1mvCOVQxjBIlZIFriRJc+Pzaa2UE3dumro2RFMc/P7uQH4+3IciKx48ArA9/dM885dmzsoBn7lwZru3Z03eGaVM5HNK779BBikWUL+/uFmXfqVOyU+H77+Gmm+QmokoVd7fKJTQAK9+RP78srgkMzNaP2Vk/9sgRyYq3erVU+Fm79mJCo/ZBfzAxvj0Ok0DSqQBGxA4mIqkxi07fwIGzMiTakEjCiaBXQDiO28No2lTW2zRokLXLctdc9r59kmyMGdMZ6/c4IcRhGQMJyC+imRvKh7/xhqwA+/57KT3pyypXlhuNVN26yQjQW295fj3i2FiYORN69YK+fb3vzjMXdBGW8g2rV8vwlV/2Z1XsqGR05gwMGQLLhkZyW/xC/nVUolL5JJoUXMf51vdy3SONqf7p8/h9982VPzxjBuebtyH2f59Q+L03AEi2/PmqVF/eP/Qsh7iW/PllEWhqQL755ivzRbijIlNiotQW+KbPvwxL6ErbpGmYipWwDkTJN92xzSjV8eOwcqWkK8xLkpIkocikSVK16bvvZP7bk5w4IWs0evaUv9kzZ2Tdhg/SVdDKt23ZAnXrQp8+cgedAzntOSYlyU6nPn2g7qFZzLTaYZlkLgx0BgfL/OOLL8q2m1atJDAFBEjkKl0abr1VhmhfeUWG0C/7m5z/0Uam7b2JbXP+ZcffFmXZT+t8ESQ3Cad6lzBatJDFvc6q1ZuRlSslv0mFjVMZ73iUoECD34D+0vtasyZX24xy3JNPrXr1yCMQFJTt8/qUKVPg5Zcl2cgbb8jfhicEuWnT4IUX4OhRmcNu2NDdLXIqDcDKdyUlSUaj3btl+1HJki479cKF8r62fWMcN4cFM+GmDyj3dR/5pp+fBNSPPpLctaky2/8aGSlDtakJKj7/XObHunaVY7zyCnzxBQYLAyQQQGtmsti/OXfcAevZSkjlIwQUO3fJYXNTq/dycXFy3WtGRHJ+9kI2lWjKy++VosXCHljDhklx+lzKVU++f38YMAC++QaefTbXbfF60dGSQWrsWKm2dO217ktSc+yY/F8eP1565qNHy42zj9MArHzXZ59Jj2vcOFl85QI7dkD37rB8xjGGFPyAhwN/I9/erVibN10aQHMy9JpZgN62TU48a9aFhxLzF+bdV04yc5bFyc37OUAZGhVcRLMCc4gsX4tN1StwXZVEInvfmePrPXZMTjl9Oiz+I5YnY0YwmF6STjIkBCub13m13m2Oe/IjR0re5C5dYNQo31twlRsHD8poizEyJF+njtw9uvCGlRYtpNLTu+/K0HM212p4Kw3Ayjft3Qs1a0qwmjnT6W+4cXHw1RORnJo4l3KOAzzi+IXAhLNYTz0l25+KF3d+hqe0vWR/f3kz690bjOHsNWWJjz5LocRzWCSTQCAtmMOKwDuoW8eiXj0ufFSrJodLTJRBhMTESz8/dUrW9EybJnXOmyXP5fXgr2iaOIegxBgM0rPG4YCBA2UBTRZkpXdbsecs0ntXyrQnP3EiPPSQVDSaPPnSUQd1UUyMZAL79VcICZGh4LfekuBsJ2NkePmrr2DECChaVG5IS5aUlc55iAZg5ZvWrpVhxmnToFw5p57q7FnodWckH65pSjBxEnzuuEPeYKpXd+q5r5BekE9MhJ9/5tT7gyn097YLc9B7y9/CFw+sZOPqeBqs/JyV5+sQQDx12UAE4axAft4imXzEEEg8QZwnnEU8yjhmVOxKqcda8szpjyk78WOs9u0ld/I77+Sop5+V3m22e8Dnzskq4CpV5K4hX74stSVP275dtmb9/LPcrMyeLav7cismRo45YgRs3Cj7rqdOlaIoeZQGYOW7XDCfdeKEjNo1WzWIgda7+CUnyRzv++9nuefnMml7yA6H3CB06QKbN0OtWheelvpX/2fz/qy+ux/l9y3lvs/uuPJ4gYES7G++WT5P/V3nsKefld5tjuaAt26VXlzRollui0LWTnz2mYzg5MsnQXnWLFkYFRYm/5Ytm/7PJiTI0HZ8vNz8nDghN2cnTkgv99VXZVooj98QZRaAdZxGeZ9Dh2SBUu/eTv/jPngQ2jc7yxs7nqNytzb4fRV4secXHu7Uc+dIWJj0SC8PjjVryqrTHj1g9GjZo2tZ3H79EW5/HYiqBGWHynUtXAgzZkgWkaQkOdblQTYsLEdD7FnZc50aZK+6CnrXLpl6eO01qFEj221RyMjB8OEXvy5ZUm6yvvhCVu+DzFds3SqPd+0KK1ZI2stDh+QGuGlT+T9XrJisxwgPlz1zOgd/VdoDVt7ngQckQGzc6NRi7nv3wkNNj/Llv62px1qsMWOgUiX3VfGxw+UrrdMbPs7Kc3LItr3Kx45J7+zUKend+3C+YLeIj5e/rxUr4PRpudkFydG8fz+EhkrPuGxZ+Ru8I53REwXoELTyJVOnShL6Dz6QeUgn2boVnr5zD2OPtaRCQBSOiROgbVunnc+lsjJ87MTFZLnO1nX+PNx1lyRfiYjw+X2kyrtpAFa+ITpahhqvuUbefJ2UYm/1anilxU5mnm5M0YKJ+M+e6Z29XV9kjFS7+ukn+OUXWfmslAfTcoTKN/TqJVl9Ro1yWvCNiJAprVOFyxPctgX+K/7U4OtJ1qyRPd/vvafBV3k9XYSlvEfXrpI5x0k1ZVd/HsmBbl/wcNnH6P9nSwqG/uSU86hcaNBAhijyQAYl5ft0CFp5vsREpydWODg5kuL3NSGQBIy/P9aSJdrz9SQrV0pxhbxWWEF5PR2CVt7t5ZelQHdyslMOHxsLy58dTQAJALJFJyLCKedSObBvH7RrJ1tc4uPd3RqlbKMBWHm2uXMlsX65cjkqNXg1xkDvzntpdmKC7Ft0ODx3j29edPq0pJc8f14ynuWR/MEqb9A5YOW5Tp2SfYfVqsmiGyf44gtoMfVFgoL9sH78RTIDeeseX1+TmCh1bbdtgz/+cH3KT6WcTAOw8lxvvglRUbB8udTVtdnSpfD669C5+WhavLsHGjey/RwqFyZOlBzFX38t+36V8jEagJVnOnpUkm68/bYUrLfZgQMwst1sKldowfCJpfErbHM1GJV7Dz8s9WvtKBKglAfSAKw8U8mSsGWLVFOxWXw8fN3kZ36OfoSDL31O4cKv2H4OZ8l1FilvMGWKTDtUr67BV/k0XYSlPM8ff0gRgFKlICjI9sN/0mkVvXY9xdHqd1C633O2H99ZUvMoR0XHYoCo6Fh6Td7E1PVR7m6afRYulAQbTkwzqpSn0ACsPMvMmXD33fDtt045/ISPo3hscgdiCpem5JLfvGpV7dA5Oy4pYgAQm5DE0Dk73NQim23cCB06wA03wOjR7m6NUk6nQ9DKc5w4Ac89J3Vru3Sx/fDbRi3nljcfoZhfNP4Rq6BECdvP4UwH0injl9njXmXfPrnxKlxYRkCcMPWglKfRAKw8x2uvyeKrmTNtH3qOXRhJhWfvIpDz+Pn7Y8WesfX4rpCVWrpea+BAyYiybFnGBeCV8jE6BK08w5QpMHas1B29+WbbD7+kx0wCTDwOkrFSi8x7me4tqxIS4LjksZAAB91bOq8mssuMGAGLF8ONN7q7JUq5jPaAlWcoXRrat3fK4psVkw/QcM3nkknLwmszXaWudnb2KmiXrbROSpKe72uvQdGicNNN9p9DKQ+mxRiUexkjKSCdJOacYU3Ju2kQtwRr9GiCD+zRTFeZSF1pnXaxV0iAg0Eda9kbhI2RHN8jR8KPP8Ljj9t3bKU8SGbFGLQHrNzr9dcl//KwYU4JxLPajOSB2Dns7PYlNzyh9WOvJrOV1rYFYGPgrbck+L79tgZflWfpHLByn99+g88+k6FIJwTfNeN20DriLbaWb8UNH79g+/F9kdNXWicmwtNPw8cfwyuvwKBB9hxXKS+kAVi5x65d8NRTcMst8OGHth8+Jga+7b6T4/6luG7B904d5vYlGa2otm2l9YkTkmyjXz8YPtwpFa6U8hY6BK1cLy4OHnhAhp5//dUpyTB694ZvDral87xWlKsSYMsx80IayO4tq6Y7B5zrldZnz0JICFxzjSTcKFw4ly1VyvtpAFaut3Yt7NwJEybAddfZcsi0wbHhzmMUnHKel196liZ32Rd80wam1DSQgE8FYaestD52DO65B+rXhy+/1OCrVAoNwMr1GjWCvXulN2SDtMExKCaeAdPeJ9iKY2WHcOAGW87hksVJHqJD3VD7rmn/fmjRAvbsgXffteeYSvkInYBRrrNtG/z0k3xuU/CFS4Nj919+o0rybrqF9+LTtfttO4dPp4F0lp075WZr/36YMwfatnV3i5TyKNoDVq5x7pzM+x45Im/ENub6TQ2CT8z9g2eOjuPXom1Zd8t1WDYGR59OA+kM8fHQqpWkl4yIcEp2M6W8nfaAlfMZAy+8AFu3wrhxtifaL1MkhFv3bKLf+hEYoO2ZOdwctc3W4GhbGsiEBDh0CM6cke1XviYmRq4rMBC++QaWLtXgq1QGNAAr5zJGki2MHQvvvQfNm9t+iu4tq1Jz1TEMflhAQFIit0dtuSQ4Tl0fRaPBC6nYcxaNBi/Mdg3dDnVDGdSxFqFFQrCA0CIh2csONWkSNGkiC5BKl4ZChcDfXwIWSErGGjWgQQP5HfXsCdOne1eQnjMHataUBBsAd90FVX0gT7VSTqJD0Mq5Vq+WLFcvvyx7g5yguqMEH+3vyJvWFwQST6J/APUev5cmKcHRrhXMWVqctHcvzJ0LkZGwfDnMmgXXXy/bcOLipNzi9dfL0Oy5c7I1B6QC0I03ymNHjkiiiu++k+pQAKNGSRKL226T53nS/tkjR+CNN2R0o1o1qFPH3S1SyitoLmjlfAsXSv5lJwQNc/Yc/11zM8N5jXcm1qXYXxFX5HpuNHhhuvO3oUVC+LNnU3sacuyY3GCMGiW91pIlJVi+/770CrMrNlaCeY0a8nWjRhLQQXrPd94JDz8sH+40ZQo884wMqb/zDvTqZXspSaW8meaCVq43aZKsdL7jDmhqU5BLx/b7+1A9did3vFWLYq3DoPWVRRZcsoLZGEmt+dJL8OqrUKXKFdm3spXIIyTkYvAFqZO7e7f0rJculeHeEiUkABsjgb5JE7nxCLBn73OmEhLkPEWKSDu//vrS9iqlrkp7wMp+qVtOwsPlcyelgTw1ZwUFW93GlFIvcu+BLzLsYDulBxwfL0FnzhyYMUOu8exZKFAg3afbXmXIGOkl58sne2yrVpUh6sKFZQ75rrugdWv7itsnJcGqVXKtM2dKoP/664tt0VSfSqUrsx6wB00kKZ+wfDl07CjzlL/+6rw35vPnOffw0+ynLFUnD8p0dNvWQvbJyZLBq0YN6NpVFlGdPCnfyyD4QuaJPHLEsiT4AlSqBMePSw/8/vvlNXjhBVi3Tr6/fr3M0U6YIMPa2b3p7tlTFo7ddpvk7S5eHG699dK2KKWyTYeglX3++kt6XaGh8Mcftm83Smvz139SNXonYzpO4+nbCmX6XNvSK/73nwS4VaukePzs2dCyZZYCkNOHwQsVkhufjh0lwKbNNLZli6xM/uQT+bpECZmX/u03KFZM6vH+/PPF67AsmcedMkU+N0Z61G3byt7eokXtabNSeZwGYGWfUaMgf36YNw9Klcryj2W3yEFCAnT6tikFyuxm/pjyWTqHLekVixSR6/vxR3jkESkmkUUuTeRhWdIrTvXoo/DQQ7Bpk9w8rFolWapSnT8Pp07J58bIh8MBBw9CmTIwZIj9bVRK5W4O2LKs14FnAANsAroYY+Iyer7OAfuopCR5w05Oljft0KwHumzPjSYl8fMry3nkq8ZMmwbt2tlxAVexYAE0bCjBN4fznbbPASulvIJT5oAtywoFugL1jTE1AQfg5j0RyuVmzIDateHwYdlmlI3gC9mfGz3Rfzidv7qDHndEuib4jhghxQQGDpSvczjfmetEHkopn5PbIWh/IMSyrAQgH3Ag901SXiE5Wba+9OsnqQYTEnJ0mOzMjZrde8j/v97McrTl5Z8a5uh8WZacDN27S0KM9u1tqeRja5UhpZTXy3EP2BgTBQwD/gUOAqeMMXPtapjyYKdPy2Kffv3gscdkj2oOt7tkNAd6xePLl3OuUXMSky2i3vmScuWduPI2NhYefFCCb9euslgpf37nnU8plSflZgi6KNAeqAiUAfJblvVoOs97zrKsNZZlrTmamlZPebcePWQv6KefyoKkkJwvJMrSFqHISEx4OAUO7yHQSuCp5v/l+HxZcviw3FR8+il89lm2FlsppVRW5WYf8F3AXmPMUWNMAjAZuO3yJxljvjHG1DfG1C9ZsmQuTqfcLjFR/v3gA0kv+dprud4DmqW50YgITMo8sb+VjP+yiFydM0OpK4ErVJBVwq+95pzzKKUUuZsD/hdoaFlWPiAWaAboEmdflDrfO2+erAguVkxSTNrkanOj4+Ju4l6CCCCeJIc/K0vfSBPbzp7i5ElZ6Xz//XKDUSjzvcVKKZVbuZkDXglMAtYhW5D8gG9sapfyBMZIson69WW+t2JFCcYutPKDL/nzk920CP6dj297jM4Pvc8LuwOzXU4wU/HxMqf9zz+SaEIppVxAc0Gr9B04IIn+ly6VwDtwIHTu7Nq0gwcPcqp8NTYl1uL+9h8TXO3iGgLbKhkZA08/DaNHw08/SdIKpZSyiVZDUll3+rQMv5YsKcHpiy+k3FxgoGvbYQznnniRwMR4XrpuCEFVL13AZ1sKx8GDJfj27avBVynlUhqAldi5U4aZlyyBv/+WRP9Ll7qtOWb8L+SfN43ujiGcuCcO/8s63ralcKxWTXrA/fvbczyllMoiDcB52cmTMHWqbCmaNk0S8L/+usvnea9w5gzxz7/Kem4lutvTFAxZTWyaPB85rmSUVkyM3GTce698KKWUi2k5wrwkOVkS8W/fLl/v2gVPPQUrVkjCiT17ZLVzJmX1XOFEQkGecIzl45qj+WpIcftTOP7zD9xwA4wfb1eTlVIq27QH7OtWrYIdO2D+fFnRfPSo1IodORLq1ZNasbVre05N15gYevTIx6SzrVg7VnJg2JrC8dQpaNMGzp2DunXtOaZSSuWABmBv9/ffUvt1zx7YvVv+DQ2F4cPl+w89JD2+YsVki80990gNW5DiCXXquKvlVzp2jPM16sDRfrzR/Vlq17b5+AkJkmJyxw6YM0fmf5VSyk00AHu6mTMhMhIOHbr4UagQLFok33/2WVi8WD4PDJQ6sGkzjo0fD4ULy5CrJ6dUjIwk+fkX8Tt6iP/KNOTTfk44x5tvwty58N130NSGLUxKKZULGoA9wbFjsHHjxY+9eyWoWhb8+qsE0VKl5OPaa6FKlYs/O2iQ9OwqVpSer99l0/oNnVw1yA6RkRAejl98PIn40/+ts/bXPjAGiheHt96SVc9KKeVmGoDd7Z13JIimKlNG5mTPnZPFUJ9/LvtUM+q9hoW5pp3ONG0aJj4eC3BYhoZxEYDN12VZss1KKaU8hK6CdrV162TYeP16+frhh+HDDyXP8pEjEBUFv/9+cSVy4cKePXRsg6TgfAAk4sAvKBDCw+07uDHQpYsMPSullAfRHrArxMbKUPLIkbBypZTva9xYVuHedJN85GEf+PVlKQ0Z9tBaar8Wbm+vfsQI+OEHaNAAWrSw77hKKZVLGoCdLSkJbrxR5nWrVZP6so8/DkWKuLtl7rdsGTtXRfPee214qHMLao+zOUBu3gzdu0Pr1vDii/YeWymlckkDsLP89x+ULSvDx+++KzVmw8M9Z7+tux07RvJDD+N/ND/lS7fkiy8C7D1+XJwUjyhcGEaN0t+7Usrj6Byw3YyBb76RbT9jx8pjXbrAnXdqEEiVMi+bdOgo9yf8wqgxAfYPCIwdC5s2yQK2UqVsPrhSSuWe9oDtdPo0PPccTJgAzZvrnGNGPvsMZs7kDYbT9M263HmnE87x9NNQtarMtSullAfSAGyXtWsvZp363/+gR48r9+Qq2L0b8/bbzAlqx+LrX2H1BzYf/9gxOHNG9kVr8FVKeTANwHb57z+Ij5cEGo0aubs1HstUrMRnNb9jyObWzP3ZIijIzoMbqV28fLmk5HRzUQmllMqMBuDciI6WgNu+PXToIDmWQ2yqU+trli+HGTOYabXj9fWPM2wY1Kpl8zm++07KKn70kQZfpZTHs4wxLjtZ/fr1zZo1a1x2PqeKjZWAu3Yt/PuvpDlU6Vr8w3Rue+Y+/JMSOU8Qz1Sbw5gtTewdod+zRyL6bbdJoQUd/ldKeQDLstYaY+qn9z19l8qJpCR45BFYtkxW2WrwzdDU9VEc/mwk/kmJkmqSROoUHcv0jVH2ncQYWfzmcMD332vwVUp5BR2Czi5j4OWXYcoUWc374IPubpHbTV0fxdA5OzgQHUuZIiF0b1n1Qv3e336czVdbIjBYJOJHgp8/yytXZcqcHfbV+D1/XhZd3X8/lCtnzzGVUsrJNABn15w58PXXssq5a1d3t8btpq6PotfkTcQmJAEQFR1Lr8mbAOhwrR9DR/XgtKMgjySMo/Y1q9nUojjrQqtjRcfa14jgYPj2W/uOp5RSLqBjddnVsiXMmHFpBaM8bOicHReCb6rYhCSGztkBJUsyr3JLWibOYXbppox+9A7WhVYHoEwRGxarGQM9e4KvrCtQSuUpGoCzavZsyS1sWdCmjWa1SnEgnZ5sSHwc8fsPsP+QP90O/MS2fNW4puMa/AKS5fsBDrq3rJr7k//2GwwZAgsX5v5YSinlYhqAs2L5cujYEd58090t8TiX92T9kxIZMX0Ik3/uwX2t40iIczDs29OUL+uHBYQWCWFQx1q5n/89flzm4uvVgzfeyN2xlFLKDXQO+Gq2bZMeb7lyF3M7qwu6t6x6cQ7YGP43ZwTNdq9mUJXhrNkczKxZ0KpVKV7D5nzMb7wBJ05InV9//W+slPI++s6Vmf37Zc43KEgWX5Us6e4WeZzUnuzQOTvoPO0rHtw0n/E1e/HO5lf54gto1coJJ120CMaMgT59oHZtJ5xAKaWcTwNwZgYMkGxXS5bINheVrg51Q+kw5WtYMZE9tdrTedMHdOsGL73kpBM2agTDh8veX6WU8lKaCSszCQmwfbsTcib6mMhIaNYMExtHLMH0bbSAIYvDcDiccK74eAgMdMKBlVLKfpoJK7t27pT5xYAADb6ZOXMGevWCefMw5+OxMAQSzwd3RTgn+C5bBpUrw8aNTji4Ukq5lgbgy8XFwb33wt13yz5Tlb5du6BhQxg6lKNxBYkzgSTgwC84kKCW4fafLy5OKh05HBKElVLKy+kc8OX69IGtW2XRle71Td8ff0CnTuBwsGnYHMIHNuOWQg0Z9VgEZTqHQ1iY/ed87z3YsUNeF610pJTyARqA01q8GD7+WFYPtWjh7tZ4pu+/l57oTTcxvcsUHni7IhUqwIjfwyhT2QmBF2DDBvjwQ3jySX1dlFI+QwNwqtOn5Q2+cmV5s1cXRUZCRASEh8Ptt2OeepqhoZ/So1t+7rhD6lIUK+bE848ZAyVKSJ1fpZTyERqAU50/DzVrwjvvQP787m6NbTKrVJQlkZFw552y+jg4mIQ/FvBc0rf88B48+ih8951sk3aqjz6Cbt2cHOWVUsq1NACnKllSiiz4kEwrFWUlCG/aBM8+KzcngImP58cuEfywJ4z+/aFvXydPk+/ZI4uurrsOypd34omUUsr1dBX00aPwwAOwb5+7W2K7TCsVZebQIalzfNNNsHcv+PtjHA7ikgP5cV84Y8ZAv35ODr7JyTIl0Lix7MdWSikfk7d7wMbA88/DrFnSnUsj10O3HiC9SkWXPz51fRS/fzOZKlvXsO+G2jR/6SE63FAE1q+HPn0w3V5n2agdLHkvgsVWOANnhBEe7oLGf/01LF0qi74CAlxwQqWUcq28HYDHjpUVRB9+eEnCjVwP3XqIMkVCiEonCKdWMJq6PoqJn47n+3HvEJCUAEssHotPhm6d6LB9O2vWO3jrPli8OIzq1cOYPBmqVXNBw//9F95+G+66S3rBSinlg/LuEPSRI9C1K9x++xXl7HI8dHuZqeujaDR4IRV7zqLR4IVMXR+V62ZnR/eWVQkJuDQl1YVavEuWUPyRBxk17h2CkhJS/iMY6u1Zz8Bf9vHI4w4aNJAt0V9+KcmnXBJ8jYEXXpAh6G++0b3YSimflXd7wIMHw9mz8O23XJ43MStDt1fjCb3oDnVDccTGsOaz0YStX0Sx5HjODXifO+uGwvS1lDu4lyUVbyZ8z1r8TDIJjgB+P96BdR/dytYAWRDeowcUKuSS5or4eCn9OGiQFsBQSvm0vBuA+/eHpk3T7dZdbeg2KzLrRacG4CzNM6fdg3t5hqnYWJm/njcPypSBkBBJD3nHHbB7N9xyC21PnKBt2p/5byNwN7RtyyM9ChAVHUvd/7ZTd+V/zPmvI8u3NabkzYdYO7U05cpl+XLtExQk879KKeXj8l4ATk6GpCTp1rVpk+5TLikyn+LC0G0WXa0XnaUe8ty50K6d9AodDrlhuO8+KcN39iwULHjlCd59VwJwqVKyknnPHpg/X67b4biQ3/rUaYsw6vDtzDim727O1LhAgq87RoXmy/nkpQoXgq9LF6O9+y60bw/10y0copRSPiXvzQGPGydF3A8cyPApHeqGMqhjLUKLhGABoUVCGNSxVrYCT0a95dTHU3vIN0dt46XIX6n/3xYq7d/J3K8myhOTkyX4nj8vQTMxEVatgoMH5fsFCkDLlhfnSB0O2Rv03nsXvz9ypPT0g4LA4SA5IJAJh8O5666UxFK9ikFUKYrXOE6pB1dS9/m/+OSlCpf00HtN3kRUdCyGizcJTpnLnj4d3n8fZs+2/9hKKeWB8lY94DNn4IYbJKlDZCT4Oe/+4/IeLkgvOjWQV+w5iya71/Dt5IH4J8tzLGBTqcrUOrRLfqBvX1mhnZgoNXAXLLh0GDqlDu+FGrlpvn/ihJQy3r4dzs6LJHlhBBOOhLOCMKpXl9jetq2MWGdUOrDR4IXpDsWHFgnhz55Nbfk9AXD4sKxCL1NGbjK03q9SykdkVg84bw1Bv/++JJmYNs2pwRcuDiNnNHxbpkgIvSK+JyAl+CYDM6o1Zkzb55mUepD33pOyiBnMASfUD+Pg6AXE/B7B+sLhLBodxvYeEnSPHr34vMDAMG67LYwHe8JPbaFKlaxdgx2L0a7KGOjSRW6Ofv5Zg69SKs/IOwF4xw745BN5s7/lFpecskPd0EuHrZcvh46vwtChdG9ZlW9XPMz7Mz/BPzmJBIc/4xvey6Odwi85RtItYfxdNIxdu2DXp/D331KKd9cuSd6VlBQGSGAuUULWlLVvL/+mflSokHEvNzN2LEa7qvHjZdj5iy+gRg37jquUUh4u7wTgL76QVcKDBrn2vMuWyVzsxo2wZQsULQrbttGhTRvo15WupctSZesadtWoz0PPdaR9nVB27pR1U/Pnw6JFEB198XCFC8P118s9ROfO0putUgWqVpUAbCc7FqNd1QMPyKK4Rx+175hKKeUF8s4ccGKiBMDatV13zoULJZuTMbJY6vXXZVj5smpLhw/L9G1q0P3vP3n8uuugeXNJh1y1qgTaYsVcm5vCaaugY2Ph3Dn77xqUUsqD5O054PPnISZGep6uDL4AK1de/NzPT4JNSvA1RoLugAHSSQZpYtOmkgCjeXOoVCnzYOuKLUJXDKPbpUcPmDxZbooKF7b/+Eop5eF8fxvSp5/KyudDh1xzvpgYePppGTsOD4fgYJmADQwktYrB0qVSYrd5c/jnH/jgA1i9WhZOTZokmRgrV7568HXZFiG7/f47fP65DD9r8FVK5VG+3QM+cAAGDpRh4Guvdf75du+WRBkbN8qCojfflG5uyirmVY4w3m0p+TWuvRaGD5dyu8HB2T9VVjJteaTDh2Uh3E03uX4+XimlPIhvB+AePWTu9+OPnX+u6dPh8cdlqHnWLLjnHnk8LIwNIWH07QszZsgo9NCh8NJLkC9fzk/nki1CdkvdcnT6tMyP5+TOQymlfITvBuDISCk32Lu3TKY66xwREVCkiETUevVkDLlCBQBOnoQXX4QJE+Qp778vBZjSyyCZXS7ZImS3c+dkXH3YMLjxRne3Riml3Mp3A/CMGTLO27Onc45/eRaq7t1lhXNKr27XLkk1vWcP9Okjo9FFith3epdsEbJbgQIwc6a7W6GUUh7Bdxdh/e9/sGGDvOk7w7x5spUmKUmCcNGiF4Lv0qWS4vHoUdlWNHCgvcEX7MlX7TInT8KTT8qKM8vSGr9KKYUv9oATEmD/fqklW6qUc86RlCSLq0DmfNOscP7pJ3jmGRmFnjUr62kfc8JpW4TsdP483HuvZAF78skLw/NKKZXX+V4PeNQoyVqxZYtzjm8MdOsGS5bIv++/DwsWkHxrGO++K+uwGjWCFSucG3y9QnKyLLpavBh++OHCTYpSSilf6wGfOSMl+Ro2dF5e4WHDYMQImdQdNgyQkegunWWx1dNPw5df2lNTwKW1eJ2hd2/J9TxokOTNVEopdYFvBeCPPoIjR2RLkDPmGRMTZVz5wQelTCCyrbVDB0l6NWSIrMWy49SXlzNMTbQBeEcQPntWXofnn5ftYEoppS7hOwH44EHpkT7wANx6q3PO4e8Pf/whn/v5sWuX5Pg4cgR++02mOu3itYk2UhUoIPO++fProiullEpHruaALcsqYlnWJMuytluWtc2yrLCr/5STLF4sc47/+5/9x968WSrYnzghK52Dgzl1SrYZnT0r08F2Bl/w0kQbIDk1n3wS4uIkzaS/79zjKaWUnXL77vgZ8Icx5n7LsgKBXOR2yqWHH5bkysWL23vc/fvh7rsluJ89C8WKkZQEnTpJ5skFC6B+unUucscrE23s2SN3JfnywalTmulKKaUykeMesGVZhYA7gFEAxph4Y0y0Te3Knq1b5V+7g++8eRJdT5yQAgLlywNSrWj2bCjfZjtP/D6LRoMX2l4EoXvLqoQEOC55zKMTbRw/LjcqCQnyy3HWFjCllPIRuRmCrgQcBUZblrXesqzvLMvKf/mTLMt6zrKsNZZlrTl69GguTpeBZcskreGECfYed/lyaNVKVlklJUmVIyS75YcfQpF6+0iquttplYi8KtHG6dPQvj3s2wfTpkG1au5ukVJKebzcDEH7AzcDrxpjVlqW9RnQE3g37ZOMMd8A3wDUr1/f5OJ8VzJGlh2XKQNt29p6aGbOlGFnkNXPERGscoTxzDNQqNJJCt156T5jZyyQ8opEGyAL4LZtkywkjRu7uzVKKeUVctMD3g/sN8akVp2fhARk15k8WTJevPde7koLpadtWwgJuVDL9+iN4dx7L5QuDYVbr8ZyXHkv4fELpOz2119yE1S1qsz/PvCAu1uklFJeI8cB2BhzCPjPsqzUSclmwFZbWpUVCQlSaOHGG+GJJ+w7bnIyfPUV1K0rK6wGDuT87wto80EYp07J1tZyZdIfOPDoBVJ2MkZqKtatK71ekBXPSimlsiy3q6BfBcalrIDeA3TJfZOyaPNmOHZMJmVzsNUlwyxTI0fCK69IQOnUCdMwjGefgFWrpMNdqxZ0T/TCSkR2iYmRZNfjx0tCkvvuc3eLlFLKK+UqABtjNgBO2ISTBXXrSnWdQoWy/aMZZZnK/+8emr/9tqzmffhhQJJr/fSTjHKn7vVNnZf16jSROfHvv5L2a8MG2W/ds6cm2VBKqRzy7iwJORz2TC/L1Pnz8ZR67S0ICoLvvgPL4o8/JIviAw9ITd+0vGaBlJ02b4a9e6XWcuvW7m6NUkp5Ne8OwDmU3mKpZ1ZP5aZ9W2DcOChThpMnpZDPjTfC6NF5uKN38CAsXAiPPAL33CMB2O7ixkoplQflyQCcXpapRZXqU8aK58lOnQDp+R45IrUX8l+xuzkP2LdPNjyPGiWLrlq0gJIlNfgqpZRNfK8ecBakzTJlGdnru79MJYp8NAQsi8WL4dtv4fXX4WbXbqxyvwMHpKZilSryS3jiCdnjW7Kku1umlFI+JU/2gNMuonpwxndUP3eY2G9G0b5uKHFxUkGvQgUYMMC97XSp+PiLRYwnT4aXXoK33oJy5dzbLqWU8lF5MgBDyiKqdX9A5C+ScvKWCoAs7t2xA+bM8fGh5+PHJd3m8uWSzjMgQOZ6y5SBqCj7E5sopZS6RJ4cggZg6VJ49lmZ34yIgMhItmyBwYPh0UdlytNnJCfD339f/Prll6FECSmx+NFHktSkSRNJuQkafJVSygXybA+YIUMk+ALEx5O8KIJnZ4ZRqBB8/LF7m5Zlyclw9KjM2x48CHfcAQUKSPf9u+/g0CH5OHBAEmgcPgzXXCNlG8uXh9tuk2pPIXkkg5dSSnmQvBmAExNh7Vrw85P9RYGBTDkRTmQk/PijB603Sk6GLVtg40YZFn7oIZmcnjlT5mgPHrzYawVYswbq1ZOgvHkzXHstNGgg/9aqdbE+b4cO7rgapZRSaeTNAOzvL0Ft2TLYsYMjNcLp8lgYd90Fjz3m7sYhw8U9esCSJTJXm6p6dQnA114Ld94JoaHyUaaMVIlILQP46KPyoZRSymPlvQB88qRk0LrmGujYEWPguXulI/nVV25IuLFzp2w2joiQOdmnn5Zh5A0bpCJTeDjcequsRk5dFVa/vnTVlVJKea28FYCNgfvvlznPmTMB2XEzbZpMCVeu7MK2bN4MAwfCxInSripVZDU2SG92zx4XNkYppZSr5a0APG2abLUZMQKA6Gh49VWoUwfeeMPFbXn2WQnCPXvKfG7Zsi5ugFJKKXfKOwH4/Hl4801J7vz884DEvsOHpcZvDioaZs/69dLN/vxzWeX1/fdQqhQUK+bkEyullPJEPhuAL6/3+83B+dy4Zw/MnQv+/mzcCN98A127ypSq06xbJym1pk+XuecNG2QbUPXqTjypUkopT+eTiThS6/1GRcdigAMnz+E3YQIHm7SQ4Af06iXxsG9fJzXCGOnx1q8vq5kHDJD6xSnnV0oplbf5ZA/48nq/xvKjw6PDqBKczCxg0SKYPVuK/Th1BHj9enjwQfj66xzXLlZKKeWbfDIAp633W+b0EU6EFCIuIJitidIx7dFD1jy98ooTTr5vHyQlQaVKslUoMDAPFxNWSimVEZ8cgi5TJCW1ojEMnz6UX8b3AmMoUySESZNg9Wp47z0nZGBcvFiGnB97TCJ9UJAGX6WUUunyyQCcWu+365/jqR+1jWXl6xAS6M/rTavyzjuyEPrxx208oTHwxRdw111S5GD0aA28SimlMuWTQ9Ad6oZSYvVyGv35MwZ4Zt10ar3wKFvXhrJrF8yYAQ6HTSc7f16qC40aJZmrxo6FQoVsOrhSSilf5ZM9YIDbJ4/CAiwgODmRhv9sYcAAaNwYWre28URJSbK1qE8fmDpVg69SSqks8ckeMABnzlxS7WhcVDiHD8OUKTaNDsfGSjc6Xz7480+Z71VKKaWyyHcD8LJlsGABrF7NiZvCee3hMDp2hLAwG46dlASPPAKnT0vtXQ2+Simlssn3AvCxYxIgS5WSRVF33UX/rtJh/d//bDpH9+7Slf7sMxsnk5VSSuUlvjcH3K+f1MU9fRqA3bulzODTT0PVqjYc//PP4ZNP4LXXJI+lUkoplQO+FYD37JEEzw8/fGEx1LvvSqGFfv1sOP6MGdCtG7RvDx99ZMMBlVJK5VW+FYD79YOAAIm6SB2E8eOl1GCZMjYcv2JFuPdeGDdOh56VUkrliu8E4E2bJDC++uqFaNujBxQvLlO2uXLqlCTbqFkTJk2C/Plz316llFJ5mu8E4PnzoUgRibrAwoXyUO/euayDEB0tS6d79bKjlUoppRTgSwH49ddh1y4oVgxjoH9/6Qi/+GIujhkfDx07ynFbtrSrpUoppZQPBGBjJEDChdqCERGwdCn07AnBwbk49muvSe3CUaPgzjtz3VSllFIqlfcH4Llz4YYbJCFGigEDoHRpePbZXBz3jz9k/9Jbb0l1I6WUUspG3p2IIzlZ5mavu+5CD3XxYvn49NNc9n7j4uD222HgQFuaqpRSSqXl3QF40iRYvx7GjJHC90id32uvheeey+WxO3SQ/b5aVlAppZQTeO8Q9NKlUgawYkXo3BmQ9M8LF8Lbb0NISA6PO2uWpJhMTtbgq5RSymm8MwBHRkLz5pL3OSoKVq0CZO63VCl4/vkcHvfECXjmGVl0lZhoX3uVUkqpy3jnEHRExMUAmZQEEREsN2HMnw/DhkmFwBx59VUJ6rNnXxjSVkoppZzBO3vA4eESIB0O+Tc8nAEDoGRJeOGFHB5z8mT4+WdJY1mnjo2NVUoppa7knT3gsDCp9RsRAeHhrLDCmDsXhgzJYZbI2Fh46SWoW1czXimllHIJ7wzAIEE4LAyAAXdDiRISQ3MkJAR++UUOEhBgXxuVUkqpDHhvAE6xapXkzBg0CAoUyMEBTp+W0oXh4XY3TSmllMqQd84BpzFggGSgfPnlHPzw4cNw/fUwcqTt7VJKKaUy49UBePVq+P13ePNNKFgwBwd46SUpNai9X6WUUi7mlQF46vooGg1eSJNOh/EPSaBC4wPZP8jcubLyuX9/qF7d9jYqpZRSmfG6ADx1fRS9Jm9iz/YAYneXIn/9PQyc+xdT10dl/SCJidJtrlRJyhgqpZRSLuZ1AXjonB3EJiQRs60MfkEJFKr3D7EJSQydsyPrB1m3Dv7+Gz78EIKCnNdYpZRSKgNetwr6QHQsAEXCt1Og7j78ghIveTxLbrkFdu+GMmWc0USllFLqqryuB1ymiFRZsCwIKBJ7xeNXtXMnGAOhoVpsQSmllNt4XQDu3rIqIQGOSx4LCXDQvWXVq//wnj1QqxZ8/LGTWqeUUkpljdcNQXeoGwrIXPCB6FjKFAmhe8uqFx7PVI8e4O8PnTo5uZVKKaVU5rwuAIME4SwF3LSWLoVJkyRzh879KqWUcjOvG4LOkeRk2W4UGgpvveXu1iillFLe2QPOtj174L//clksWCmllLJP3gjAVarArl05rFWolFJK2c/3h6DXrJHMVwULgp/vX65SSinv4NsRKSoKmjSBt992d0uUUkqpS/h2AH7nHen9vvqqu1uilFJKXcJ3A/D69TBmjKx+rljR3a1RSimlLuG7AbhvXyhaFHr1cndLlFJKqSvkOgBbluWwLGu9ZVkz7WiQLc6cgX37ZM9v4cLubo1SSil1BTu2Ib0GbAMK2XAsexQsCBs2yPyvUkop5YFy1QO2LKss0Br4zp7m2GDPHoiOli1HgYHubo1SSimVrtwOQX8KvA0k574pNnnmGbjtNik5qJRSSnmoHAdgy7LaAEeMMWuv8rznLMtaY1nWmqNHj+b0dFmzaJF8PP+81vpVSinl0SyTw56iZVmDgMeARCAYmQOebIx5NKOfqV+/vlmzZk2OzndVxsAdd8gQ9O7dEBzsnPMopZRSWWRZ1lpjTP30vpfjHrAxppcxpqwxpgLwMLAws+DrdPPmwbJl0Lu3Bl+llFIez3f2Ac+fD+XLw9NPu7slSiml1FXZEoCNMRHGmDZ2HCvHPvxQsl8FBbm1GUoppVRWeH8P2Bj491/5vFgx97ZFKaWUyiLvD8BTpkDlyhAZ6e6WKKWUUlnm3QE4ORn69YNKlaBBA3e3RimllMoyO1JRus/EibB5M/z8M/h796UopZTKW7y3B5yUBP37w403woMPurs1SimlVLZ4b7dxzRrYtQvGjweHw92tUUoppbLFewPwrbdKAC5Xzt0tUUoppbLNewMwwHXXubsFSimlVI547xywUkop5cU0ACullFJuoAFYKaWUcgMNwEoppZQbaABWSiml3EADsFJKKeUGGoCVUkopN9AArJRSSrmBBmCllFLKDTQAK6WUUm6gAVgppZRyAw3ASimllBtoAFZKKaXcwDLGuO5klnUU2GfjIUsAx2w8njvptXgeX7kO0GvxVL5yLb5yHWD/tVxnjCmZ3jdcGoDtZlnWGmNMfXe3ww56LZ7HV64D9Fo8la9ci69cB7j2WnQIWimllHIDDcBKKaWUG3h7AP7G3Q2wkV6L5/GV6wC9Fk/lK9fiK9cBLrwWr54DVkoppbyVt/eAlVJKKa/kFQHYsqxWlmXtsCxrl2VZPdP5vmVZ1vCU7/9lWdbN7mjn1ViWVc6yrEWWZW2zLGuLZVmvpfOccMuyTlmWtSHlo6872poVlmX9Y1nWppR2rknn+x7/uliWVTXN73qDZVmnLcvqdtlzPPY1sSzre8uyjliWtTnNY8Usy5pnWdbfKf8WzeBnM/27crUMrmWoZVnbU/7/TLEsq0gGP5vp/0VXy+Ba+luWFZXm/9E9Gfysx7wuGVzHhDTX8I9lWRsy+FlPe03Sff9169+LMcajPwAHsBuoBAQCG4Ealz3nHmA2YAENgZXubncG11IauDnl84LAznSuJRyY6e62ZvF6/gFKZPJ9r3hd0rTXARxC9u15xWsC3AHcDGxO89iHQM+Uz3sCQzK41kz/rjzkWloA/imfD0nvWlK+l+n/RQ+5lv7AW1f5OY96XdK7jsu+/xHQ10tek3Tff9359+INPeBbgF3GmD3GmHjgF6D9Zc9pD4wxYgVQxLKs0q5u6NUYYw4aY9alfH4G2AaEurdVTuUVr0sazYDdxhg7k8U4lTFmCXDisofbAz+mfP4j0CGdH83K35VLpXctxpi5xpjElC9XAGVd3rAcyOB1yQqPel0yuw7LsizgQWC8SxuVQ5m8/7rt78UbAnAo8F+ar/dzZdDKynM8imVZFYC6wMp0vh1mWdZGy7JmW5Z1o2tbli0GmGtZ1lrLsp5L5/ve9ro8TMZvJt7ymgCUMsYcBHnTAa5J5zne9toAPIWMqKTnav8XPcUrKcPp32cw1OlNr0tj4LAx5u8Mvu+xr8ll779u+3vxhgBspfPY5Uu3s/Icj2FZVgHgN6CbMeb0Zd9ehwyB1gY+B6a6uHnZ0cgYczNwN/CyZVl3XPZ9r3ldLMsKBNoBE9P5tje9JlnlNa8NgGVZvYFEYFwGT7na/0VPMBKoDNQBDiLDt5fzptelE5n3fj3yNbnK+2+GP5bOY7l+XbwhAO8HyqX5uixwIAfP8QiWZQUgL/44Y8zky79vjDltjDmb8vnvQIBlWSVc3MwsMcYcSPn3CDAFGaZJy2teF+RNYp0x5vDl3/Cm1yTF4dSh/pR/j6TzHK95bSzLegJoAzxiUibkLpeF/4tuZ4w5bIxJMsYkA9+Sfhu94nWxLMsf6AhMyOg5nviaZPD+67a/F28IwKuB6y3LqpjSS3kYmH7Zc6YDj6esum0InEodUvAkKXMmo4BtxpiPM3jOtSnPw7KsW5DX6LjrWpk1lmXltyyrYOrnyGKZzZc9zStelxQZ3s17y2uSxnTgiZTPnwCmpfOcrPxduZ1lWa2AHkA7Y0xMBs/Jyv9Ft7ts/cO9pN9Gr3hdgLuA7caY/el90xNfk0zef9339+LulWlZ+UBW0+5EVqH1TnnsBeCFlM8t4IuU728C6ru7zRlcx+3IsMVfwIaUj3suu5ZXgC3IKrsVwG3ubncG11IppY0bU9rrza9LPiSgFk7zmFe8JshNw0EgAblLfxooDiwA/k75t1jKc8sAv6f52Sv+rjzwWnYhc2+pfy9fXX4tGf1f9MBr+Snl7+Av5M27tKe/LuldR8rjP6T+faR5rqe/Jhm9/7rt70UzYSmllFJu4A1D0EoppZTP0QCslFJKuYEGYKWUUsoNNAArpZRSbqABWCmllHIDDcBKKaWUG2gAVkoppdxAA7BSSinlBv8Hr6hH8BWglosAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(x, y, 'o', label=\"data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "ax.plot(x, res.fittedvalues, 'r--.', label=\"OLS\")\n",
    "ax.plot(x, iv_u, 'r--')\n",
    "ax.plot(x, iv_l, 'r--')\n",
    "ax.legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS with dummy variables\n",
    "\n",
    "We generate some artificial data. There are 3 groups which will be modelled using dummy variables. Group 0 is the omitted/benchmark category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.296798Z",
     "iopub.status.busy": "2021-02-02T06:52:02.295402Z",
     "iopub.status.idle": "2021-02-02T06:52:02.306952Z",
     "shell.execute_reply": "2021-02-02T06:52:02.308110Z"
    }
   },
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "groups = np.zeros(nsample, int)\n",
    "groups[20:40] = 1\n",
    "groups[40:] = 2\n",
    "#dummy = (groups[:,None] == np.unique(groups)).astype(float)\n",
    "\n",
    "dummy = pd.get_dummies(groups).values\n",
    "x = np.linspace(0, 20, nsample)\n",
    "# drop reference category\n",
    "X = np.column_stack((x, dummy[:,1:]))\n",
    "X = sm.add_constant(X, prepend=False)\n",
    "\n",
    "beta = [1., 3, -3, 10]\n",
    "y_true = np.dot(X, beta)\n",
    "e = np.random.normal(size=nsample)\n",
    "y = y_true + e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.312963Z",
     "iopub.status.busy": "2021-02-02T06:52:02.311477Z",
     "iopub.status.idle": "2021-02-02T06:52:02.322276Z",
     "shell.execute_reply": "2021-02-02T06:52:02.323400Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.         0.         0.         1.        ]\n",
      " [0.40816327 0.         0.         1.        ]\n",
      " [0.81632653 0.         0.         1.        ]\n",
      " [1.2244898  0.         0.         1.        ]\n",
      " [1.63265306 0.         0.         1.        ]]\n",
      "[ 9.28223335 10.50481865 11.84389206 10.38508408 12.37941998]\n",
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 2 2 2 2 2 2 2 2 2 2]\n",
      "[[1 0 0]\n",
      " [1 0 0]\n",
      " [1 0 0]\n",
      " [1 0 0]\n",
      " [1 0 0]]\n"
     ]
    }
   ],
   "source": [
    "print(X[:5,:])\n",
    "print(y[:5])\n",
    "print(groups)\n",
    "print(dummy[:5,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.328606Z",
     "iopub.status.busy": "2021-02-02T06:52:02.327273Z",
     "iopub.status.idle": "2021-02-02T06:52:02.348012Z",
     "shell.execute_reply": "2021-02-02T06:52:02.349040Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.978\n",
      "Model:                            OLS   Adj. R-squared:                  0.976\n",
      "Method:                 Least Squares   F-statistic:                     671.7\n",
      "Date:                Tue, 02 Feb 2021   Prob (F-statistic):           5.69e-38\n",
      "Time:                        06:52:02   Log-Likelihood:                -64.643\n",
      "No. Observations:                  50   AIC:                             137.3\n",
      "Df Residuals:                      46   BIC:                             144.9\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             0.9999      0.060     16.689      0.000       0.879       1.121\n",
      "x2             2.8909      0.569      5.081      0.000       1.746       4.036\n",
      "x3            -3.2232      0.927     -3.477      0.001      -5.089      -1.357\n",
      "const         10.1031      0.310     32.573      0.000       9.479      10.727\n",
      "==============================================================================\n",
      "Omnibus:                        2.831   Durbin-Watson:                   1.998\n",
      "Prob(Omnibus):                  0.243   Jarque-Bera (JB):                1.927\n",
      "Skew:                          -0.279   Prob(JB):                        0.382\n",
      "Kurtosis:                       2.217   Cond. No.                         96.3\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "res2 = sm.OLS(y, X).fit()\n",
    "print(res2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare the true relationship to OLS predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.353894Z",
     "iopub.status.busy": "2021-02-02T06:52:02.352453Z",
     "iopub.status.idle": "2021-02-02T06:52:02.632129Z",
     "shell.execute_reply": "2021-02-02T06:52:02.633324Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFlCAYAAAAzqTv+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABrfUlEQVR4nO3dd3iT1RfA8e/btKVllj0KCIjssqcIFBCKgIigKCroTxEVFyIoOBBEBQVFBWSI4lZEoSAomzLLLrss2WXv1ZEm9/fHbWnBjrRN8iblfJ6nT9vkTd6bpM3JXecYSimEEEII4V4+ZjdACCGEuB1JABZCCCFMIAFYCCGEMIEEYCGEEMIEEoCFEEIIE0gAFkIIIUzg686TFStWTFWoUMGdpxRCCCFMs2nTprNKqeJpXefWAFyhQgU2btzozlMKIYQQpjEM43B618kQtBBCCGECCcBCCCGECSQACyGEECZw6xxwWqxWK8eOHSMuLs7spni1gIAAypYti5+fn9lNEUII4QDTA/CxY8coUKAAFSpUwDAMs5vjlZRSnDt3jmPHjlGxYkWzmyOEEMIBpg9Bx8XFUbRoUQm+OWAYBkWLFpVRBCGE8CKmB2BAgq8TyHMohBDexSMCsNksFgt169alZs2a1KlTh88++wy73Z7hbQ4dOsQvv/ziphYKIYTIbUyfA86q8KgYRi/Yw/GLsZQJCmRQWFW61gvO0X0GBgayZcsWAE6fPs1jjz3GpUuXGD58eLq3SQ7Ajz32WI7OLYQQ4vbkVT3g8KgYhszcTszFWBQQczGWITO3Ex4V47RzlChRgilTpjB+/HiUUhw6dIgWLVpQv3596tevz5o1awAYPHgwK1eupG7duowdOzbd44QQQoi0eFUPePSCPcRabTddFmu1MXrBnhz3glOrVKkSdrud06dPU6JECRYtWkRAQAD79u2jZ8+ebNy4kVGjRjFmzBjmzp0LwPXr19M8TgghhEiLVwXg4xdjs3R5TiilAL1P+aWXXmLLli1YLBb27t2b5vGOHieEEE515QqcPQuyBdHreFUALhMUSEwawbZMUKBTz3PgwAEsFgslSpRg+PDhlCxZkq1bt2K32wkICEjzNmPHjnXoOCGEcJp586BzZzAMSEgAX696S7/tedUc8KCwqgT6WW66LNDPwqCwqk47x5kzZ3j++ed56aWXMAyDS5cuUbp0aXx8fPjxxx+x2fQQeIECBbhy5cqN26V3nBBCON2FC/Dkkzr4AgwdCkmjdsJ7eNXHpeR5Xmevgo6NjaVu3bpYrVZ8fX3p1asXAwYMAKBfv350796dGTNm0Lp1a/LlywdA7dq18fX1pU6dOjz11FPpHieEEE5lt8M998CePfDOO/orTx6zWyWywVBu/NTUsGFDdevCpOjoaKpXr+62NuRm8lwKkYtduACFCoGPD/z1FwQHQ/36EBcHERFQtarMA3sgwzA2KaUapnWdVw1BCyHEbenPP6FaNZg0Sf9+//06+AJcvw733Qfh4aY1T2SPBGAhhPBUp09Djx7w0ENQtqweer5V4cKQNy8cO+b+9okckQAshBCe6K+/oEYNmD0bPvgA1q6F2rX/e5xh6OAsAdjreNUiLCGEuG3kzQuVK8M330DNmhkfW64cHD3qnnYJp8m0B2wYRoBhGOsNw9hqGMZOwzCGJ11exDCMRYZh7Ev6Xtj1zRVCiFxKKfjhB/joI/1727YQGZl58AXpAXspR4ag44E2Sqk6QF2gg2EYTYHBwBKl1F3AkqTfhRBCZNWxY3pP75NPwsKFkJioL3e0zOjgwTBnjuvaJ1wi0yFopfcpXU361S/pSwEPAKFJl38PRABvOr2FLnbu3Dnatm0LwMmTJ7FYLBQvXhyA9evX4+/vb2bzhBC5mVLw7bcwYIAOul98AS++CBZL5rdNrVo117RPuJRDc8CGYViATUBlYIJSap1hGCWVUicAlFInDMMo4cJ2ukzRokVvlCIcNmwY+fPnZ+DAgTeuT0xMxFfSuwkhXOHgQXjhBWjeHKZOhTvvzN79nDmjF2u1bw/lyzu3jcJlHIosSikbUNcwjCBglmEYtRw9gWEYfYG+AOW95A/jqaeeokiRIkRFRVG/fn0KFChwU2CuVasWc+fOpUKFCvz00098+eWXJCQk0KRJE7766issWf30KoS4fdjtsGgRhIVBpUp6dXPdujrBRnadOAHPPgszZkgA9iJZ6toppS4ahhEBdABOGYZROqn3Wxo4nc5tpgBTQGfCyuj++/eHpM6o09StC59/nvXb7d27l8WLF2OxWBg2bFiax0RHRzN9+nRWr16Nn58f/fr14+eff6Z37945abIQIrfavx+eeQZWrICVK/W+3uSEGmkIj4pxLPVuuXL6u6yE9iqZBmDDMIoD1qTgGwjcC3wMzAGeBEYlfZ/tyoa628MPP5xpT3bJkiVs2rSJRo0aATqndIkSXjkSL4RwJZtNz+++8w74++t53+bNM7xJeFQMQ2Zuv1EDPeZiLENmbgf4bxAOCpJkHF7IkR5waeD7pHlgH+B3pdRcwzAigd8Nw3gGOAI8nNPGZKen6iqpiyn4+vpit9tv/B4XFwfomsFPPvkkI0eOdHv7hBBepEsX+PtvvdJ50iSdxzkToxfsuRF8k8VabYxesOe/AdgwZC+wF8p00kEptU0pVU8pVVspVUsp9X7S5eeUUm2VUnclfT/v+uaao0KFCmzevBmAzZs3c/DgQQDatm3LH3/8wenTevT9/PnzHD582LR2CiE8SGKinu8F6N0bfvpJbxVyIPgCHE+j9nlGl8teYO8jqSgd0L17d86fP0/dunWZOHEiVapUAaBGjRp88MEHtG/fntq1a9OuXTtOnDhhcmuFEKbbsgUaN4aJE/XvjzwCjz/u+L5eoExQYJYu59tvZS+wl5H9Namkt9gqMDCQhQsXpnndI488wiOPPOLCVgkhvEZ8PHz4IYwcCUWL6l5pNg0Kq3rTHDBAoJ+FQWFV076BrH72OtIDFkIIZ9i4Ua9oHjECevaEnTvhgQeyfXdd6wUzslsIwUGBGEBwUCAju4WkvQoaYPdueO89OHs22+cU7iU9YCGEcIbLl/XX3LnQqZNT7rJrveD0A+6tDhyA99/XtYGLFXPK+YVrSQ9YCCGya+XKlO0bbdrofb5OCr5ZljzcLSuhvYYEYCGEyKorV+Cll6BlS/jqK4hNWpmcJ495bUpOxiErob2GBGAhhMiKhQuhVi0deF99FaKiIDCdlcnuJMk4vI7MAQshhKNOnID779c5nFetgrvvNrtFumZwRASEhupe8PHjZrdIOEh6wIDFYqFu3brUqlWLhx9+mOvXr2f7vp566in++OMPAPr06cOuXbvSPTYiIoI1a9Zk+RwVKlTgrKx0FMJ91q/X30uXhvnzda/XA4Lv9QUrsd3TCvtbb6PatoVx4+CXX8xulnCQBGD0Pt8tW7awY8cO/P39mTRp0k3X22y2dG6ZsalTp1KjRo10r89uABZCuMnp0zqJRpMmeugZoHVrCAgwt13AmvGbudzpESx2Kz4oSEjQW6GykOxDpBIbC5Mn6xrNbuKdATgyUm90j4x0+l23aNGC/fv3ExERQevWrXnssccICQnBZrMxaNAgGjVqRO3atZk8eTKg80G/9NJL1KhRg06dOt1ISwkQGhrKxo0bAZg/fz7169enTp06tG3blkOHDjFp0iTGjh1L3bp1WblyJWfOnKF79+40atSIRo0asXr1agDOnTtH+/btqVevHs899xzKjX8gQtyWlIKff4YaNSA8HD74QAdeD3D2lI15Nd+g8cuNyaPiSfTxx4oF5esP+fPrakuJiWY30/vs2aMX1q1a5b5zKqXc9tWgQQN1q127dt18QatW//2aMEFfd+2aUnXrKuXjoxTo73XrKjVtmr7+zJn/3tYB+fLlU0opZbVaVZcuXdRXX32lli1bpvLmzasOHDiglFJq8uTJasSIEUoppeLi4lSDBg3UgQMH1J9//qnuvfdelZiYqGJiYlShQoXUjBkzkh5KK7VhwwZ1+vRpVbZs2Rv3de7cOaWUUu+9954aPXr0jXb07NlTrVy5Uiml1OHDh1W1atWUUkq9/PLLavjw4UoppebOnasAdebMmcyfSyFE9vTqpd9jmjRRaudOs1ujlFLKblfql1+UKl7MrmYaD6oN9fqo2OPn1fYpa9RgPlKrRq9RavJk3e6jR81urne4dEmpH39M+T3pPdqZgI0qnZjofYuwLl1KSXBut+vfcyg2Npa6desCugf8zDPPsGbNGho3bkzFihUBWLhwIdu2bbsxv3vp0iX27dvHihUr6NmzJxaLhTJlytCmTZv/3P/atWtp2bLljfsqUqRImu1YvHjxTXPGly9f5sqVK6xYsYKZM2cC0KlTJwoXLpzjxyyEuEXy+4qPj05mUb8+vPwyZFKW1OUiI7k0/R92/radoac+oWLju6i88HdC6um376JN76ILA7i6uxJ0S9oLfOxYjtJg3hbmz4e+fSEmRk8x3HUXJL1Hu4vnBeCIiPSvy5tXDwu1bavnO/z99e/NmunrixXL+PbpSJ4DvlXqkoRKKcaNG0dYWNhNx/z9998Ymcy5KKUyPQbAbrcTGRlJYBpbGhy5vRAim/79F/r0ge7d9TBkz55mtwgA++pIVKtQCtoSaAb83LY6DRZ8hMWS8tZdonJBSrKOZXuioVx3feHRo9C0qTmN9nTnz8Nrr8EPP0D16rB6tQ6+JvC+OeBmzWDJEp1vdcmSlODrYmFhYUycOBGr1QrA3r17uXbtGi1btuS3337DZrNx4sQJli1blkaTm7F8+fIbZQzPn9eVGwsUKMCVK1duHNe+fXvGjx9/4/fkDwUtW7bk559/BuCff/7hwoULLnmMQtx2bDb47DMICYHNm6FAAbNbdMO+FSc41O5ZLLYEDAAfC43bFvhPh9wS6M9Zn5JYTqbq9cpe4LQlJuoPJr/8Au++q1ezm/hBxfN6wI5o1sxtgTdZnz59OHToEPXr10cpRfHixQkPD+fBBx9k6dKlhISEUKVKFVq1avWf2xYvXpwpU6bQrVs37HY7JUqUYNGiRdx///089NBDzJ49m3HjxvHll1/y4osvUrt2bRITE2nZsiWTJk3ivffeo2fPntSvX59WrVpRXqqeCJFzu3bB00/DunXQubMuHegBw7YJCfDxxxA47HNetu/DbvHDwI7h76/3+qbhXN6y5Dt3VCfjKFoUrl1za5s93pkzeoTU1xdGjYI774Q6dcxuFYZy44rahg0bquRVwcmio6OpXr2629qQm8lzKUQWLFumtxh98QU8+qhHbN/Z+ud+Php8md/316d3t6t8+vpxihnnUhJtpNPxWF+uG4VP7eGuhJ16BbcHPBaPoBR89x0MGABjxugV4m5mGMYmpVTDtK7zzh6wEEJkx/r1sHYtvPKK3lZ08CCkWuthlmuXEll6/1juXTmUwX61eTx8LV0eyA9U0QdkMuJ3+s67OXg8gMpK1ovccPCgXmS1eDG0aKG/PIz3zQELIURWXb8Or7+uA9lnn6UM0XpA8N3+xvdcKXoH9698g70Vwrhz60y6PJC1IHqw+0Aetf/CmTPAt996zCIy03z/vc7XvXatztkdEQFVqpjdqv+QACyEyN2WLdOLrD77DJ59FrZu9YjAe+4cfNd0ErVGP0VJ23Hsvv7U+fkNClZ3sP5vKsmFkI4eRff8fv/99k7GUawYtGoFO3fCCy/orWUeyCNa5c556NxKnkMh0nDqlN7T6+OjA/GkSVCokKlNUgpmfXOeGjVg//pzKAwMwG5LZOevf2XrPqvEbuUoZYmds0gvJLPb4eRJ5zbckyUk6J0xH36of+/UCebNAw9fsGp6AA4ICODcuXMSQHJAKcW5c+cI8ID8tEJ4hA0b9PeSJfUb8dat6a4gdqeYnRdZUOE5mvepxh1Bx9l1fzHiff1INHywWnz54GoJwqNisny/Je4qRFliiN97OGUl99GjTm69h9qwARo2hKFDYffulFzOXjAXbvoirLJly3Ls2DHOnDljdlO8WkBAAGU9YAuFEKY6fVovsJo+XWc6CgvTiXtMZrfDwn7h1JnSj3bqFFGhA8jXciub48vyeMEPaXpkO2vLh7C5ZBWOLNhD13pZG4YuGlIGOwa2w8egXGN9YW7fC3z9ug66Y8dCqVIwezZ06WJ2q7LE9ADs5+d3I0WjEEJki1Lw00/Qvz9cvQrvv+8ZxRMiIznz22KOfr+MDpeW8W/+Opz46S8aPtCAQ4PnAbA5uDqbg1O2Dx6/GJvl0xh5/DljKYnvyWN6QvjOO91a1ccU+/bBl1/qDGaffJLp1EJ4VAyjF+zh+MVYygQFMiisapY/6Dib6QFYCCFyrHdvHYCbNYOpU3UVo0y4+g05cWUkqm1bClsTKAxEhz5PtQVfYvj7AVAmKJCYNIJtmaD/pqJ1xLm85ch34SgULgz79+ek6Z7r4kXd033ySZ1IY98+uOOOTG8WHhXDkJnbibXq0rIxF2MZMnM7gKlB2PQ5YCGEyBa7PaWAQqdOuje0cqXDwXfIzO3EXIxFkfKGnJ3517RsD/+Xo/f+D8Majy82LBao3r78jeALMCisKoF+N+eVDPSzMCisavbOeeeDrLB73l5Xp5k9W7+2zzyT8gHDgeALMHrBnhvBN1ms1cboBXuc3coskQAshPA+u3dDy5YwYYL+/dFHs1S5yFVvyNcvJ/JXqzHc+WAIJRKOYvj6gsWSZhrJrvWCGdkthOCgQAwgOCiQkd1Cst0j29pxCEOuvYPNBrz5Zu7ZC3zqFPToAV27QvHiem9v5cpZuov0hvWzM9zvTDIELYTwHlarnu97/329l7do0WzdjSvekNdN2UrgK324P34jW+/oQsV/vsJy8UiGaSS71gt22hBouXKALZETMT6UPXVKjwZ4O6tVF0s4fhw++ADeeAP8/DK/3S2cPdzvLBKAhRDeISoKnnoKtm3TPaIvv9TbjLLBGW/I4VEx/D1lJndu28yKcw/Sdk8EvS1H2Dnsd+oMfShpG0ywSwvHpJ7H7rxuL/EMZPvKbZQtVw5OnNDJOHy98G0+JgbKlNHB9osvdBaratWyfXeDwqreNAcMORvudxYZghZCeIdLl3Qt19mz9TajbAZfyPn8a3hUDL998TvjprzO62t+YMaeR1leoz4r5y2h5nsPu2UP6q3z2Efz+mHBzoal0d6bjMNm09uKqlSBadP0ZV265Cj4gvOH+53FCz8aCSFuG0uX6jq9AwfqYdz9+yFPnhzfbfIbb3ZXQY/9YTtjfvuBPHZdH1wZ8TQvtogxUWV4MCzHzXPIrfPYp0oWBGD/hj3QvZ6+8OhRjyix6JAdO/SWonXroGNHuPdep969M4f7nUUCsBDC81y4AIMGwTff6N7PSy9BQIBTgm+y7Lwh2+2w6KXZ/DKxH6U4gdWwYKBItPiytnyIWxf13Hqu00UKYccg6OxpvQ84NNRjcyD/x7hxulhGoULw8896AZkXZLLKKQnAQgjP8uefOuCeOaMX3QwbpoOvyfbuhWnd/mLkzq7s8K/Bsx3fISD/tZQsVsHVCXbjop5b57Ftvr6cNEpSLvYMVK2qc197uuTaxZUr63n9sWP1SmcTJM+nx5yLI7hogFsSdUgAFkJ4jmPH4LHH9H7PefOgfn2zW4Q1QTH1nUO89mVF8ubpSKfeUzjbrx0H5+4m1mq7kcXK3Yt6khcW+V25RJHrlzhUJJiJhZ/lct6qPOG2VmTTlSvw1lt6FfuwYbpgxn33mdac8KgY3vhpN2VnXaX3pXVsv68IQ64lAK5N1CEBWAhhLrtdF01v317PV0ZE6OT62dhu4mz/fvwHfkMH0zPhDOs6H2DklKKULv2svtLfz9TUhl3rlqHs3zNp+MErALT8YCFrW77JjrX5+AJ0cpLixeG779zWJof8/Tc8/7z+sPXaa2a3BqXgzU8uUHyWL3Pje+JPPNbpvjz+6IeMXuDv2tdUKeW2rwYNGighhLhh3z6lQkOVAqWWLTO7NTdcu5yoVtZ7WdlB2UHZLH5KrVpldrNS7Nun1L336uetQAGlPv5YKatVDX/XqspxRMXHK6XatlWqaVOzW5ri9GmlHntMt7l6daXWrDG7RerIEaU6d9ZNmpC3j7LreKysho/6uGVvVeHNuTk+B7BRpRMTvWSGXgiRqyQm6oQaISF6f+/XX+sC6knCo2JoPmopFQfPo/mopU5LEemI5XOvsK94M+6JGgeAAfhghxUr3NaGDFmtutDE+vU6E9iFC3qu3NeXDrs/5wjlOb77ss7M4UklCWNiIDwc3ntPv+Yu3B+dGbsdvvoKataElUsSqNBpL8u63InN8LlRGnJt+RCXJ+qQIWghhPt16gQLF8KDD8L48TrpQhKzEudfOK8Y9IbBN98U4OeCDSjQszOVpo/Sxd7TSCXpdhs3Qr16emj+hx/0QqsyZSA2Vg/h16hBQGW95ejslmNUKFvW/GQchw/roPvqq1C3Lhw5ku3sZc6ye3fSbqfVVibd9Sm9E6ay4I2/GLSoLj0e//jGorroCrUY6eI5fekBCyHc4/p1nWgBoG9f+OMPmDnzpuAL5iTOjxixkpiS9VgzbQ9vvAEPnpxIpWlDYckSGDFCfzerx3bhgn6+GjVKSU7RunXK8xYXp5NVzJlDoZo6AF/eeVT3gO12HYTdzWbTGaxq1oR33tGpJMHU4JuQoLNZ1qkD/ts3cap8Y57ZNwS/hnXpXKM4I7uFcKpmfSY268GpmvXdkqhDesBCCNdbskQHkVdfhVdege7d0z3UnYnzz/62mPP93iH0wjpi/Csw6+vzVH0q1QHNmpkXeJWC337TNY7PndPJSNIqsBAUBHnzwtGjFO/QDYC4/cegWx3o1cutTQZuTqhx330wceJ/PmS52/r1ukmFtq8kqvjbVD+7CiNfKf0B8MEHAehazP2lCSUACyFc5/x5HTimTdN7PWvXzvQm7kicrxSse3gMjf98g6IobD6+lJwzleAw8+Yl/+PFF3XwatQIFizQQ7hpMQzd2z12jLyVy2DHQB05Ck2egSZN3Npkrl/XvXPwiIQa167Bu+/qzninIpHM8g/DciZWV8367ju98t5EMgQthHCNf/7R+3l/+EGXx9u2zaF5VGfXyb3V/v3Qpg1Y/wzHQGEAFkPhu3m9U+4/R6xWPacLKQUnIiPTD77JypbV23r8/BhdZiwReZLyYSql79PVoqL0ufLm1Xm6o6P1fm4Tg+/ChVCrFkwbe4Fl1fsx4/FZWGwJKQds2mRa25JJABZCuEbevDowbNgAo0ZBoGM9WFclzk+0Kub1+J7na64kKgouDvhAt8licckiqyyv5F67Fho00F020O1xtMZxqhXPK+u9yuJrST350qX1CmlXuXQJXnhBJ0z56Sd9WZs2UKyY686ZiXPndNGssDDoHP8np4rUoOXuKeQpkEe/zi56vbNDhqCFEM5ht8OkSbqA+vDhelvRhg3Z6gU5O3H+zrkHufzYc3S6soi85XpTdW0LypQJhYeWZFivN7uytJL70iUYMkQ/d8HB0LJl1k/41ls6lzJQo9hprq86BtTX88PHjuXgkWQgPFwPk588qRNqJM2lmiE8KoZP5u9h/9ogLi6tRcnYM+yo+jI198zSK8enJmVV69jRJa93dkkAFkLkXHQ0PPssrF6tux42m+5pmJlQPzIS64KlrJt7mnqbpqIMH6L6TKD15OdTxv5ctMgqo5XcNwXgiAg9VHvqlF6cNmIEFCiQ9RPeddeNHx/e+yHvXJrG9euXyVu2rGv2Ar/6qh4er11bB+JGjZx/DgeFR8Uw8Lu9lP0zlmdi1rG6yFmevuNrqu3/Gz7+GAYMSNmGZeaiujRIABZCZF9Cgh5e/vBDyJ9fL2zp3dv8SjaRkdhat8UnPp57sLOvxN0UX/wb9ULKueX0Dq/kLlUK7rgD5szR6Tez6/RpXcSiY0d87ihHwcgr7N99mcrlysGiRdm/39SU0nuK/fz06uZSpfQCOzekDE0ulHBr2k+7HV5//wol5/ryT+JD+JGA9bIvz977Dv90uJ/pb/R2edtyQgKwECL7Dh+Gjz7S24o+/xxKlDC7RVw8Fc+yJ37g/vgEfLFjw4eEh1sQlMXgm96bviPSW8ldtqC/fp62btUrw6tVgzVrcv6B5dQp6NcPfv+dwLv0XuAzUceo7KxkHHv36m1kLVroXnqHDvrLQTl5LtMbzo855Mv0sSU5svJOJgb0IyAxDgPAlkjtk/uZWNH8Qh6ZkUVYQoisuXwZpk7VP991lx5+/uUXjwi+yz9azengurQ+8AtWw5dEw4cEX18+uFYyS+ksk9/0Yy7Gokh503f0PtJayd3g7EHm/DRQz5eePq0TaIBzRgvK6qDLsWMUqqU/aFzedVQviBo8WI9UZEdCgh7dqF0btmyBSpWyfBc5fS5vHc5XNoOTKyrwco9i+GzZzJagOrSPW4bdzWkknUF6wEIIx82Zo3tax4/rPaYhIVCxotmt4tS+y0R1GEKHA19x2FKefve+S1wpUmr1lqzCkVvnXzPg8BxuOpKPGb1gDxdOX+DdDb/x6JqZGMWL6206Dz/s3GH6oCDIlw+OHqVYJ52MI37/MWj9TMq+3KzaskVPJ2zfDg89pOd8S5fO8t3k9LlMPWwff6IQ5/6pjfVMQTqUm8Hfx3sSX7gorzz0DscCCrk1jaQzSAAWQmTu5Em9SGjGDL258s8/dfA1mVLwyxdnCB1Qn/Yqhg13v8qTDVtwPTAA4EatXshaJi1nZOO6sZL7zBmo8bxepDZqlA6WzmYYN/YC56lYhhcL/kjBwOZ0UQouXtQL4goWzPp9XrkCs2frVJfZ5MhzmdEQdZmgQI6eTqDK/Ks0it7K5oCLRHUrxbX6hTAS3iagf3/aHLrO6AV7mBhcnTJBgYx0c2nI7JIALITImM2mt8YcOaKT6Q4apPdRmuzYryuY/85qvjkQCuUep8WYrjTq0ZTCo5ZyPYeZtHKcjevkSV2paNgwXZN3zx4oUsTh82dL8l5gPz/WV3mCoheBixf0eceMubFNKUN//62rPo0apZMm79uX40IOmT2XmW3Zah9Um8UfbmfO1W4EEIsRB30CPuD+jk9BPZ3Jqmvhwl4RcG8lc8BCiLQdOJCynWjCBJ3J6u23TQ++iVbF6nuHUuaxVjx14B1W+LWl568PUL5HU8A5mbSyfR92O0yZAtWrw+jRsHmzvtzVwRf0CvR//gGgRcGtlNodkTI0ndle4FOndNrITp3gr790zxecUkUps+cyvSHqj2b9y9NPw7AXitFfTSGQWHwAhcHAAhe8MuDeSgKwEOJmCQl6ZXONGrpoKkC7dlClirntAnb9fYj1xe6j+ZIRGIAvdvzsCfisiLhxjDMyaWXrPqKjdfKR557TqSO3bXPv/tjg4BvD208efp83j75409B0mpTSq7GrV9eFCYYN0x8asrMXOR2ZPZe3DlErBdd2l2LTZ42Z8f11tlftTvdrv2EYBlgs+AQGUP3J9It5eBMZghZCpFi3TpeN2bFDL7x56CGzWwTo9MiLHhhPm0WDucOAPZ1eo8rSSenW6nVGJq0s3YdSOnfz8ePw7bc6F6K790Lv2gU//givv05i6XJU/HcRly5BoVRpKv/jzBm9KjskJKXn7gIZPZeph6gTr+Th/KJaxO4rRb4yl1m5sgC13rfDUyOheXNYtcpjslg5gwRgIYQ2erQumhAcnOOFN860fLlev/TcvoNUCG5J+b8nEW2z8Gm5ilTetZH9NRrSMaA8Xc1o3KpVOtVhvny6+k+pUuZtxzp8WM/ddumC5Y6yFFx1hV3RlylUtuzNyTisVr1trFcv3dbISKhaFXzMGRAdFFaVwX9u58zGYC4sq8adtgOML9YT+1cjqVuvqe6ZJ3+YadHClDa6igRgIW53yfO8TZro3L4ffpj1FbNZ4GhShqvzlnO07/tMP96NxIovUvefUdQO8yV8y3G9aKdQJRY20/tSV6aXZ9lVLlzQRQ6mToX339cFFBwotehSqfYCB1bRe4HPbTkKjz+e0mNct05/mtm+XW8pat/eZb1eR9XIF4zln2JU37yedwM6EGpbgU9cIH7+F/QBZmdVcyEJwELcrk6c0Dl9g4Nh7Fi90tmBQgCuyGoENwfPqH5TqDPxeaqjGOezgoSp9Qlso4NITveV5ohS8Pvv+nk7e1avCHdkdbE7JAfgo0cJqqXrAF+JPgbPh+lFVa++CuPGQZkyeoTD5Fq4VqtenD18OPzP51e+Mp7GiFO6J/7zTzrdZS6X6ZiDYRjlDMNYZhhGtGEYOw3DeDXp8mGGYcQYhrEl6auj65srhMgxux2+/lr3fObMgZIlHb6ps7MaQUrwBDi1/wr/3PUydSY+h4ECdK3ewHURN453xh7dbBs6FB59NKXM4ief6LKLniDViueirWvTxlhGlH8TiI+Hd97RiTT69dNzxSZPL2zcqNenvfUWdO4MX1SbiKH0641hwM6dprbPXRwZ9E8EXldKVQeaAi8ahlEj6bqxSqm6SV9/u6yVQgjn2L9fZ0bq21fPXW7bplMVppJRHdvMAmhm0guSMRdimTYNRtX5lbD9E9hd66F0a/WmtxfXZakHbTZdMhDgiSfgs8907d569VxzvuxKXvF86hS+hQuwLziUfWeCdB5om03nnB4/3qXTC5m5fl0PGjRpApWOLmfJ6M388Qf4f/KhS2sze6pMh6CVUieAE0k/XzEMIxrw/g1YQtyO7HadXGHqVHj66f/Mr2U2RJzT3mdaSRkKnIin7OI8PH0cWt3zDEcGNKLGg/X04qA0arcOCqt6Uxsh6/t8HbZli14VXqEC/PGHXqxU1YNTHG7erAMZ8Gi+v/DZ4g/5wnTgNdmSJfpz39kDl1hW7Q1a7p4CK7vAwNlw7736AA+q1esOhkru9jtysGFUAFYAtYABwFPAZWAjupd8IY3b9AX6ApQvX77B4cOHc9xoIUQWrFqla7aOGaN/j4+HPHnSPLT5qKVpZi0KDgpk9eA2mV6fmeQAX/3QDpoc3o46kp9nDv+OzfBlzheH6PNiHocW4+ZkHtoh16/ryclPP4WiRfXwbY8eXrUgaH/RxpyIK0yLawtMbceFC7pq4a5vI3kz/1d09JmP/9Xzuk7v8OGeM4R/7ZpOZtKvn1NfZ8MwNiml0qw16XAANgwjP7Ac+FApNdMwjJLAWUABI4DSSqmnM7qPhg0bqo0bN2ap8UKIbLp4UQ8vT56sa85u2KDTImag4uB5pPWOYAAHR3X6Tw8ZdO8zK4kuln83hyZ9HsbfloAP8G/e6uSb/Tul7q3l8ENzqS1bdHnFAwd07/eTT6BwYbNb5ZhFi+CHH2DaNLZX74Hvv7upZttlyucGpXTK8JdegspnIokgFF97gg5uySMwnmTLFj0xHRGh9xw7SUYB2KGNX4Zh+AF/Aj8rpWYCKKVOKaVsSik78DXQ2FkNFkLkgFJ6uLR6db3YasAAvaglk+ALmc+v5jTLVFwcGJNXkScp+CrDh0pvPe45wRf0qvBSpWDZMv38eUvwBf2h4aef4MQJEsuUI1gd49w59zfj+HHo1g0eflhRv9gRpr8Qga+R9KHNx0envvQEp0+nlNasWxf+/depwTczjqyCNoBvgGil1GepLk9dl+pBYIfzmyeEyLLLl+GFF/Q+z/Xr9TBqvnwO3dSRHMhd6wWzenAbDo7qxOrBbRwOvqvnXaROHXhz7YPYDD+UxYIRkAejTeZD1y6llA5anTrpxUrFi8Pq1d65EKic3v/LsWP43lGWglwhJvqy206fOhX2nr//5WDldsw735TgBxrqxVWesshKKfj+e93QF1/UhUYAypd3azMc6QE3B3oBbW7ZcvSJYRjbDcPYBrQGXnNlQ4UQGbDZdBCx2aBQIZ0+av16aNAgS3fjjDzKt7p8NoG5jYYT0rk8Za/u5v2FzfBdvRxjxAi98MbMBTcHDkCHDjor1PnzmNJddKZUyTjyVtXB+PzWdNJQOtm+fdC2LfR7LpFRxcawwyeECqfWYwwdqq9YsgQ84TXfv1/nNn/qKR2Ao6LcHnhvUEq57atBgwZKCOFkUVFKNWqkFCj1xx9mt+YmKz6JVNG+NZUCtbHaY+rqwdNmN0mzWpUaM0apwEClChRQavx4pRITzW5Vzp0/r/8OPv1UndxzUVXkXzXhC6tLTzlj3THVtcnvaojPByrMb646VLqubkOXLkodPerSc2dZbKxSJUooVbCgUhMnKmWzufyUwEaVTkyUTFhCeKtbV+r++queePMA5+eu4dT/3qT52VWc8ivHnk/n0mBAJ7OblSIxUY+VtmunSy0m9xy9XVCQTqwSH0/xyoU45leII8ddd7pPfznNnJd38c/5J/EnAauPL+vzhnBq1EQav/Gc56wa37VL93YDAvRK59q19Vy/yaQcoRDe6uGH9Qrdp57SpfAefdT0Nzyl4LNnFxBw/71UPbsau2Fh72efUNUTgu/163oI9No1/Ua8Zo3enpVbgi/o1//kSRgyBB8fGFpgLEXW/eP001y/rut2DOxVjCcu/0Ygsfhiw89mZV1wDV5TVUz/WwTg6lW9CDEkRBfLAJ3i0gOCL0gAFsK7nD6t31RAFwCIiNCrON1R8D0ThzedZW6Jxykw9Q/8ScAHBSg2zVzgcKpKl1m0CGrV0qkk583TlxUt6hlBwoWevzqamtF/OPU+ly3THchJn1zi66AneTbxWxSQaBhYLb6sLR/inrSgmZk/X7/mY8fqDCD33292i/5DArAQ3kAp+OYbqFZNBxGApk11AXiT2RIVfz/xC3kbVqfD2d+JvyMBq6+FRMMHq8WXVcE1HU5V6XRnz0Lv3rrwgJ+f/sDSo4c5bXGXKVNu1HG+VKAsBS4dc8rdXryoCym1aQNtrs7hVLGa/O/CL3zdqCuPP/Ihn7XoxeOPfsjm4OquSwvqqEGDdE83MBBWroSJE/XiRA8jc8BCeLrdu+G552DFCl0P9dlnzW7RDdELjnDukRfoeOlv9gQ15rH2fdhXsQxLY0JoemQ7a8uHsDm4OoZZPaLnntMFJ95+WxckCAgwpx3udPQozJoFiYnEFitH8XPR2O05K/c7a5berXPqlI5tH174B7/1RYkYO5nP9liItdqIrFAHcGFa0MwopXcB+PrqbU558+pqD+lkffMEEoCF8GQ//6wzBuXNq5NCPP20aYXTU4uPiGTF+xFsjzjLcyqCzb3GUu/bl7k+ZjlcjGVzcHU2B6fUmXVrj+jQIR1oS5WCjz+GYcP0HODtolw5vSH3xAnspctSfs8iTp3S28Kz6sQJncnq+Mw1fF/oM+566z4qjHgGro0Bf39C/fwY6eq0oI44eFB/2GreHN57T+/p7uQB6w4yIQFYCE+UmKg/yTduDI88AqNHZ6lsoCvtGfYrFd//H61VIi0t/sRP+4X6vboCbi6UcCubTde7fftt6NpVf3ipXNn15/U0qfYC+1YsR8GIK+yNvkzp0o5XQVIKvv1W53DucPUPVhuP4HPJDiPDoWONm/bxdq0X7P6AmywxEb74Qk/LWCw3ht69hfkfpYUQKc6f10PMyfOUd92lc/t6QPC9fDaBuY3f587hvfBT8fhiIw8JFDwWfeMYVyTycMi2bToovPaaHn4cOdK15/Nkydmwjh7F1vcFAojl0HnHg+/+/TpvRt8+Nj4o8hk/8xg+yp5yQESEc9ubXTt26HUQAwfqBu/apRdbeRHpAQvhCZTS+3hfe01nYxowQPfoLJbMb+sGqz9bR9E3+9A5cQd7yt1LlTOrwGpNM62g23tEf/6pt2AVLgy//OIR27FMVbas3vPq40NwlXzEo6eFM5OYqBcMDx2qX9Zlvb+j5Q+v62HdTZvSfb1Nk5Cgt1z9/rvu+Xrja55ehg5XfEkmLCHScOyYUu3b6+xBjRsrtWWL2S264dQppQbet0PZMNQJ37IqevRf+oo1a5T66CP93Sxxcfr7mTNKvfiiUmfPmtcWD2W/fEVN8H1FfdVtUYbHRUUpVb++UnmIVa+02a5iYpRSCQlKzZmjlN3uGa+3UkotWaLU0KEpv8fHm9cWB5FBJiwJwEKY7exZpSpVUmrcOI9Jh2hfvUZt7zhIhRVco/z8lJrdZaqKP3PJ7GZpFy4o1bevUk2aeMzz5bESEpQNQ/1efWiaV1+/rtTgwUo191mjvg/oq64WK6/spUsrde2amxuaifPnlXr6aR2yKldW6vJls1vksIwCsAxBC2GGtWvhq6/0SpeiRWHPHr3oygOc+n4+xZ7qTE1szDbGcfynpVR87Bmzm6Wl3g8zYIAeN/WQYXqPMnAgHD4MM2Zwwb8Uec7+dy/w8uV6uUH5fYtZbtyHJS4R4g347DO96t4TJJfWfPllvaf7zTf1KudAk/cZO4kswhLCnS5dgn794O67dUqhQ4f05R4QfG2Jin+e/I18Tz2EDzYMwN/HSsXDEWY3TS9O695d57ouUUJXeho92qP3eJrq3DmIjATgUsGyFLqcMgl88aLesRMaCkXjj/NPvu5YVKK+0scHYj0gi1Wyc+fgmWd06siNG2HUqFwTfEECsBDuoRTMmKEzWU2eDK+8oldtesg2mZ3b7aws0Y37fujJuXzlUXnygMWC4SmLbvLm1aUDR46EDRuyXGbxtlO2rN7Em5jImfylKB5/jAqD/qZar21UusvGd19bGTgQluwqjV+nDvqDjKfU6rXb9UiHUlCsmE5As24d1K1rbrtcQAKwEO6QmKiXl5Yurd9MPv8cChQwu1XExyneew/qNfBhQ1xtNj3xGeUvbsdn2TLza7fu2wdPPKFzXwcE6B7Q4ME6paTIWFIyjgWLothGIfxJ4NQfjdjzUwgPxk/lTJGKjO53kLz5DJg+XY/GmP16g56KCQ3VIx3JObvr1vWIESKXSG9y2BVfsghL3FYSEvTCquQFI4cP6zq0HmLzL9FqY+A9qjVL1BNP6MXEHsFqVWrUKKUCApQqVEip1avNbpH3mTdPKVDPPv+lKvFwpGrKGvWx8bpaXyBEKVBb7qil1L59ZrcyRUKCUh98oJS/v1KFCys1bZpefZ0LIIuwhHCzdev0RNvWrbr31qcPlC9vdqsAuDZ/BUeeHUHNYxFc9ynAmHcuUX+E2a1Ksnmzfq6iouDBB2H8eChTxuxWeZ/KlaFdO05dT+Ruv9VMN4bgp6xwBaY06sqo1k9zwEOmPwD9Ws+bpxPQfPmlRySecQcZghbCmS5d0slzmzXTqzZnztSLSDzEtv7fEHhfKNWPLcbXUPj/9iP1RzxodrNSvPGGnrv880/93EnwzZ4qVWDhQs5Wr0PooY34KSsGYDN8uBhYkNKF85ndQl2XOT5e//zqq7o28/Tpt03wBQnAQjhXp04wYYLeNrFrl/5k7wEZes6cgccfh0NfzMZAAXrBa97920xuGXr+8fhx/fO0afp569bN3DblEoPCqrK2SmPiff1vlIeMqlTXnGpFqS1apAtkJKcMbdcOHnjA3DaZQIaghXCmefMgJgZq1DC7JYBeSLp08EJ++OoqM+K70eZ/b8Kvi8GaYP6K14sX9X7Vb77RW7MmTEjJYyxyLiyMrsWKwcBPeCXAl8q7NrK/RkMe6dvNvOIJ58/r/dvff6976W3amNMOD2HoOWL3aNiwodq4caPbzifE7ezolnPs7vQ67Y5/z/b8zfCJXE3NWobeHxoRoYOvWSteZ87UQ/WnT8Prr+uSgblof6dHaNcOrlzRSV88wcKF0KuX3tv75pvw7ru3RX1mwzA2KaUapnWd9ICFcJZ58/Sb3dChpm6Vsa1cw/7+4yi+eQGhXGFD+7epP/MdLPmShsKbNTN3q8m4cXofdN26MHcu1K9vXltys7Jl9VCvpyhRAipV0oG4Th2zW+MRJAAL4Sx//qkDyvvvm9aEAz9HUrZXG6qqeOwYnPv0OxoN6H3j+nCziqcrpYcfixaFnj11ZZ2XX5Y9va5UrtyNZBwZ7aN12d+E3Q5ff63n9L/4Qn/gWrPGI9ZEeApZhCWEs0RF6d6cCW8w8bF2JvfdxLe9I/BJSitoWHwoHh9z45jwqBiGzNxOzMVYFBBzMZYhM7cTHhWTzr06yYEDeji0QwcdDIoV0/OAEnxdq2xZHQRPnEj3EJf9Tezdq+d3n39e1+1NXu0swfcmEoCFcIb4eP1GU6+e208d9dsethcL5X9fN6PY3XdhCfBPM43k6AV7iLXabrptrNXG6AV7XNMwm00n9q9VS6eP7NNHL70W7lGvHjz9dIZBz+l/E1arztdcuzZs2QJTp8LixZKzOx0yBC2EM+zcCYmJvHvYj58Gz3PL8O7VC1YiOo3m3sj3iTcC2d1/Mv0/6w5rg9NcZHX8YtpJ9tO7PEeOHdNbiTZsgPvv15WfypZ1/nlE+ho10l8ZcPrfxJkz8NFHejve+PE69apIl3wcFcIJ1q7YxjX/QFYUKOfS4d3wqBj6vjCOkbWeJ6ZYLTpHvs2uyl2w7I2m9tj/Eb7lOM2Xx1LxUm2aL4+96fxlgtJeZZze5TlSrJhe1fzbbzB7tgRfsygFcXHpXu2Uv4m4OJgyRZ+rTBnYvl2vh5DgmykJwEI4wetx5anVfzqHg1LedJw9vBseFcP0Ub/z+eQ3GLRzKhXtB/iofj+O/D6W/JVLZTqfNyisKoF+N9fODfSzOC8pw5o1cN99KcUTIiLgkUdk3s9MxYvDkCHpXp3jv4mVK/WK5uee0z8D3HFHdlt725EALIQTHL8YizJ8/hNsnDW8qxSEv/YP438fQR6VgC82fAw7vvmv3gjymc3nda0XzMhuIQQHBWIAwUGBjOwWkvNh8qtXdSrBe+7RQ/EHD+rLJfCar3hxOHo03asd+ZsIj4qh+ailVBw8j+ajluoPdJcvw4svQsuWet530SL9s8gSmQMWIqdsNsKnD2ZqyH38VaPVTVc5Y3g3Zvt5dnYcyHfHpnHUJ5j8XEUpsFp8WVs+5EaQd2Q+r2u9YOfOSy9cCH37wpEj+g35o488osyiSFK2rJ6Pz0BGfxPJoyrJH+ySR1VazB1K0W2b4LXXdBnDfB6QW9oLSQAWIqf27aPOoR3kD2l/08U5Hd6122H+8+E0mPo8bdRZJtz5HF/c355aZw/Q9Mh21pYPYXNwdYKTgnyZoEBi0gjCLpnjBd0tHzVKDzevXAnNm7vmPCL7ypWDBQuyffPUoyqFr1/imn9eYoH3Gj3K+EnjoEkTJzX09iRD0ELkVFQUAO0e7+C04d3oaGjRAuZ9HcPl/MGcnreR4BnvYsmfl83B1fmqWQ82B1e/Kci7fI432Z9/6nzXhgG//KK3m0jw9Uxly+p9wFZrtm5+/GIsKEWXXctZPPUFXoz8HYB5RapK8HUC6QELkVObN0OePLTpFkqbHCaXsEasZn+/Mfy+py67g97juWkvUPnx5zD8fOmadEx6WYuSv7ss09WJEzp/88yZOn/zmDFQqpRz7lu4Rtu2eu+11ZqtxCd1jKu89Men3PvvBraUrsK8avqDlstGVW4zUoxBiJy6915d2SeHf9v7PpxOpXcew4Idm2Hh0l8rKdLJxJzNyZSC777T2avi4mD4cP1zBukNRS4waxbWXr1JjLcypmUvpjW4H7uPhUA/i3MW790mMirGIEPQQuTUnXdCx47Zvvm1i1bmNh9JhXcexwc7ABYfKLItwkkNzKExY3RGpZAQ2LoV3nhDgq+3UEpXnLp4Meu3rVQJv3uas3LmUua364nysThv5bwApAcshKkWLoTvey/h51P3sqdMKFXOr8WwWnWt3iVLzKtaZLPpsnElSujvM2fCM89IKklvc/EiFC4Mn36qRy0ykpgIn38O//4LEye6o3W3BekBC+EqNlvmx6Th3NHrfNJuEWFhsCmoLVsmr6NqzDKMpUv1tg4zg2/yCrCOHfWbctGi8OyzEny9UaFCeotQBnuBAZ29qlkzGDQIjh/P9qItkTXyHyVETrz/PpQv7/AblloTyb+t+xB/RxVeXdyZka+cYMsWqNu3sU54kE4aSbewWuHDD3XZuD17oH9/sFgyu5XwZIahtyKltxc4Ph7ee09X8Tp8GKZPh/BwqVTlJjKRI0ROREVB/vwOvWGd+WkBRXp34k5lw47B8cFfMnikTl2ZXsIDwD3zbUeOQJcueo63Rw8YN04PPwvvl1EyjrNn9bDzo4/C2LE6h7dwG+kBC5ETUVGZliC02+Hrz6/h36sHPkoHWMPiQ9mCV24c4/ZSgbcqWVK/+c6apXtBEnxzj3Llbh6Cvn4dJkzQC7SCg/WUw48/SvA1gQRgIbLr7Fnds8ggAO/deJlWraDva/mYf0dflH+eNGv1urVUYLJVq6B9e7hyRddrXbwYunZ13fmEOXr10usKAJYt06vZX3pJv/6gKxgJU0gAFiK7kjJgpRWArfF25nX9muKN7qDoliVMmwY9Do7GJ2JZmous3Foq8OpVeOUVnTx/3z499ydyr9atdW3m556DNm30vPCyZXqhnTCVzAELkV2lSulAljoAR0ZyYvwfXJi5jE5xUewsHsrUOXdQrGnS9c2apbm6eVBY1ZvmgMFFaSQXLdLFEw4fhpdf1ouu8ud37jmEZ4mL0yUJp0yBgQN1IpW8ec1ulUACsBDZFxICX3yR8ntkJLYWrShls1IK2Nt9MDVnfORQWT6Xp5GElOIJefJI8YTbic2mFwmuXw+NGpndGpGKBGAhsis6GipXvrEC2r4sAsOWiAEoi4UqDQreFHzDo2IyDLBOLxWYbM4caNBAL7j55Re9NzQgwPnnEZ4pX76bPygKjyFzwEJkx9WrULOm7lEmiakcShwB2H3+u8gqeZtRzMVYFCnbjFy61/fMGejZEx54QKeTBL3aWYKvEB5BArAQ2bF1qx7STTX/u+laNToyjxMv/HeRlVu3GSkFv/4KNWroFJIjRsAnnzj/PEKIHJEhaCGyY/Nm/T1VAC48dTTz+Qw+ugwF/W863K3bjL78UmexatIEvv1WB2IhhMeRACxEdkRF6WQVqfZQ5tsXxeGAqlS9JfiC3k4Uk0awddo2I6Xg/Hmdt7lXL523uV8/SSUphAeTIWghsiM5A1aqRVblz0Vxqkz9NA8fFFaVQL+bg6HTthkdPgxhYTqpRmIiFCmitxhJ8BXCo0kPWIjsGD36ppq453acoIT9FDtC0s6K5ZJtRna7Lhs3eLD+/ZNPpGKREF5EArAQ2XHvvTf9enT2ZooCBVuln5bSqduMTp7URRNWrtQ93ylT4I47nHPfQgi3kI/LQmTV1q2wYMFNtYA3XK9Jf8ZSoWtd97ShcGF9/m+/hfnzJfgK4YUkAAuRVVOmwMMP3zT/u/JoBf4I7k+xigVcd96dO+Ghh1KKJ6xaBf/7n0OZtoQQnkcCsBBZFRWli9anmm/Nt3I+odVOuuZ8VqvO2Vy/PkRE6AxcIIFXCC8nAViIrLDZ9BB0qv2/8ScvMPHQffSyfef8823ZAo0bwzvv6FKBu3bp34UQXi/TAGwYRjnDMJYZhhFtGMZOwzBeTbq8iGEYiwzD2Jf0vbDrmyuEyfbt0wXNUwXgI3O2ABB4d/oLsLJt8GA4cQL+/BOmT9d7j4UQuYIjPeBE4HWlVHWgKfCiYRg1gMHAEqXUXcCSpN+FyN3SqAF8Yam+LLizkwLw+vVw7Jj+eepU3evt1s059y2E8BiZBmCl1Aml1Oakn68A0UAw8ADwfdJh3wNdXdRGITxH9+56WDhVekdjaxQxRjAVGuewdxobC2++qXNIv/uuvqxsWZ1YQwiR62RpH7BhGBWAesA6oKRS6gToIG0YRprvPoZh9AX6ApQvXz5HjRXCdP7+UKfOTRcVO7KZg0H1CM5J4qk1a+Dpp2HPHujTJ6V6kRAi13J4EZZhGPmBP4H+SqnLjt5OKTVFKdVQKdWwePHi2WmjEJ5BKXjjDVi9+qaLuvuEs6zdyOzf76+/wj33QFwcLFoEX3+ta/YKIXI1hwKwYRh+6OD7s1JqZtLFpwzDKJ10fWngtGuaKIQHsNth/HidgnLr1hsXHzkCUVfvokSbWlm/z7g4/T0sDAYNgu3b/5NhSwiRezmyCtoAvgGilVKfpbpqDvBk0s9PArOd3zwhPMC//0LbtvDKK9CuHTz++I2rjv60nOeZSL1aVsfv78oVePFF3eu1WvUc78cfQwEXJvEQQngcR3rAzYFeQBvDMLYkfXUERgHtDMPYB7RL+l2I3GXPHqhdW9f//fprnYIy1fBw4Myf+ZC3qVnHweUUixdDSIguotCixU3pLIUQt5dM3zWUUquA9FLutHVuc4TwENevQ968UKWK3ov71FNQrtx/Dit0IIo9eevRLH8mWamuXoUBA3QQr1JFp5G8+27XtF0I4RUkE5YQqdls8OmnUKECHDyo0z2++26awRerlXKXtnO2rAP7f/38YN06vYhryxYJvkIICcBC3LB7t56XHThQ78UNCMjw8CsbdpNHxWOvWz/tAy5e1Pd1+bIunrB+vZ7rDQx0ftuFEF5HArAQSuli9nXrwt698PPPEB4OpUtneLNjS/YAULhNGj3guXOhZk34/HNdQAF0EBZCiCQSgIUwDDhwADp31mkfH3vMoUpDi4MeohAXqdypasqF589Dr15w//1QtKgedu7SxYWNF0J4qyxlwhIi17BaYdQovQe3cWMYN07P02bBli2Qp3ghSgenurBzZx10hw6Ft9/WmbOEECINEoDF7ScqShey37oV4uN1AM5i8MVu59FZj1Lmjt4YRueUyxct0tWLKld2bpuFELmODEGL20d8vF7R3LgxnDql53k/+CBbd2Xde5B2F2ZQu8TJm6/Il0+CrxDCIRKAxe3j2291wH3sMdi5Ex54INt3dWLeZgDyt0i1AGvOHL1n2JqFrFhCiNuWBGCRu8XGwrZt+udnn4WlS+H773Nc4u/y8iis+HJHp1Q5oGfNgmnTwFdmdoQQmZMALHKvVat06cAOHXQg9vWF1q2dcte+26OINmpQJSTV1qKoKKhf36EV1EIIIQFY5D7XrunCCS1b6uHgH35wevKLs9cC2V2sRUpnNz5eD2vXcyArlhBCIKugRW5z6pTOYnXwILz8Mnz0EeTP79RTKAXdmEmXLtAj+cIdOyAxUQKwEMJhEoBF7mC3g48PlCgBnTpBjx662pALnDgBZ87o0e2bLgwK0kPQQgjhABmCFt7vn3902scDB/T867hxLgu+AJeHjmEdjW+uAdy5s86CVamSy84rhMhdJAAL73X+PDz5JHTsCBaLLvnnDmsjKcwFQurfkrzDMGQBlhDCYRKAhXeaNQtq1IBfftHJNTZtgtq10z08PCqG5qOWUnHwPJqPWkp4VEy2T134UBR78tWnUKGkC2w2aNJEt0UIIRwkc8DCOy1aBGXKwPz5uopRBsKjYhgyczuxVhsAMRdjGTJzOwBd6wVndNP/unCBktcOcqFG35TL9u7VpQYlAYcQIgukByy8g1K6TOD69fr3MWN00YNMgi/A6AV7bgTfZLFWG6MX7MlyM2IjtwBgNEi12CoqSn+XBVhCiCyQHrDwfMeOwfPPw7x5es63cWPIm9fhmx+/GJvp5eFRMYxesIfjF2MpExTIoLCqafaO95/IRzQPU/TeVNuNNm/WtX6rVXP8MQkhbnvSAxaeSyn4+mu9wnnpUvj0U/jmmyzfTZmgtJNwJF+ePEQdczEWRcoQdVrzxGsSG/MIv1OjVfGUC6OiICQk6xWVhBC3NQnAwnP9/DP07auHdrdtgwED9GrnLBoUVpVAv5tvF+hnYVBYVSBrQ9Txfy1geMBIysdEplxYrRrcf3+W2yWEuL3JELTwLDabzmJVuTI88ojuVT78sE6ykU3JQ8npDTE7MkQNwIIFvDyvAwoD494AWLJEZ92aMCHbbRNC3L4kAAvPER0NzzyjE2rs2QOFCukg7ARd6wWnu+K5TFAgMWkE4ZuGrufNQz32GAA+KEhIgIgIaNhQF3mQ/b9CiCySIWhhPqtV52yuW1cH3tGjoWBBt50+wyHqs2d1/eDOnTmrihJPHuw+FvD3h9BQGDECgoNlC5IQIsskAAtzXbigVzW//TZ06QK7dkGvXm7tUXatF8zIbiEEBwViAMFBgYzsFqJ7zDYbcfOXMbrAcMpd3c3Ux5ahho9IGX6OioLChWUBlhAiy2QIWphDKR1kg4KgQQOdzapbtzQPdXSLUE7cNER9+DB89SUnS47k5VdLMu/Cv9xVOy+rvoGGDZsBzVJuGBWle8JCCJFF0gMW7rdmDTRqlFI8YerUDIOvo1uEcmzVKujcGVWtGtYvJtC9+k7++gve/SgvGzfq6d6bnDkDMTFSglAIkS0SgIX7XL0Kr74K99yjg9fp05nexJlZrDL000/QqhXMm4c9LoGH43/EUieErVthyJB0RpiTM2BJABZCZIMMQQv3WLwYnn0WDh2Cl17Si64KFMj0Zg5vEcoJux01YADY7RiAHYMhD+ym0cxMdj8FB8Prr0sAFkJkiwRg4R6zZumVwytX6h6wgxzaIkQ254nXrYPatdm6N5CZ+cbw5pnn8MeKJcCfJm+GZj4+VLOmzkkthBDZIAFYuE54OJQurUv1ffKJ7k4Gpp0WMj2DwqreVMkIbs5iBdmodnTlCrz1FmrCBBa3eJ+Oa96hSJHetBxxF218IjBah+oVzrc+nFuC/IiKNto80AICArL0mIQQAmQOWLjCqVPQowc8+CCMHasvy5cvy8EXMtkilMTheeLISHjqKahcGTVhAj8WeoluK17liSf07qe27zTDeGtIusE39WKwS6fO0ebR9kQPeDfLj0kIIUB6wMKZlIIff4T+/eH6dT3PO3Bgju82oyxW4OA8cWSkXmRltWLH4DkmsTioLzN/h3btMm/DrUG++ukDAHxzvTAyCC2EyA4JwMJ5pk/X5QLvvltXLXJTeb4M54mVgrg4iIjAnmjDB7DhQ9d7zvH5fN0xd8StQb7mKR2AV+Yrm9PmCyFuUzIELXLGbof9+/XPDz2ke8ArV7q1Nm56qSSHhuSDTp2I6/k/3lkcSpzKQyIWfPL40+mTUIeDL/x30VfNUwc4ky8I37LOTQgihLh9SAAW2bd7N7RsCc2bw6VLuijBE0/kqHJRdtw6T1yuoD/TEzbQ/pG2WJeuYPjCZoxe2ZTf+y6BESOwLFuS5jxvRm4N8jVP/0t0qcoM6uC+DxpCiNxFhqBF1lmtumDC8OF6DHfsWLcWT0hL13rBdI07AjP/gV/nw44dbCgaxsPnJlG2eQW2ToVq1W5JI5nF+4eUkoZfPfASjzSp4PSUmEKI24cEYJE1Fy5A69awdauu0/vll1CqlNmt0ous2rZFJSSAzc5I36GMjBvGxxMMnn/eOZ3ymxeDdcr5HQohbmsSgIVjUhdPaNwY3ntPbzPyBGvXwssvo+ITMOw2rFgIvjOAXYsMypVzwfk2b9Y5oDt2BIsl8+OFECINMgcsMrd8uU63+O+/OghPmeIZwTcpt7S6+26uRh8hzu6HFQuGvz+9vw11TfAFXTziiSfcWjJRCJH7SAAW6bt0CZ57Tpfbu3IFzp83u0Up5s+HmjVR48bxS6F+lL7+L6PvW0rC2yPwjViCcXf25nodEhUFdeu6fbGZECJ3kSFokbY5c+CFF+DkSV1w4P33IW9es1ulJSRg7/cSpy/npbtaSUyh5vzxG4SFZX+RlcNsNti2TReWEEKIHJAALNK2cCEUK6bzOTdq5NJTOVRIYc0amDAB+vRh7rXWfHz9HzZcLM8L/fMwYgTkz+/SJqbYu1dn+ZIKSEKIHJIALLTkNJJVqkDTprp4gp9fOoVwncehQgqzZukkH3Y7ib/O4EO1nKu1mrF8tq7z4FZSA1gI4SQyiSV0jd4OHXQaySlT9GV587o8+EImhRTsdpgwAfXooyi7HQCl7Hx4bwSbNpkQfAEeeQR27oQaNUw4uRAiN5EAfDuz2eDzz3Vd2zVrYPx4vcLXjTIspNC/P7z0ErsD6hJHAIlYsAT40+b9UPz93drMFBaLDr6+MngkhMgZCcC3sx9/hNde06ucd+6EF190+8reW3Ms+9msFIy7SumCgfyQ93me9f+exra1zBuwFJ8PR+CzNOtpJJ1GKf18rVxpzvmFELmKfIy/3cTH64VEISF6L2vRotC5s2l7WgeFVWXIzO1UP7SDB3cso8XBzewoWoMBef/hyR3+dOxYg50ToXx5N6xwzojVCp9+qkcMqlaFFi3Ma4sQIleQAHw7WbMG+vSBc+fgwAGdx/n++01tUtd6wRRbv4qmHw3GYrehgK8uD+BCUX9++QUefdQD8l2sXav3Q2/bBl26wOOPm9wgIURuIEPQt4MrV+Dll+Gee+DaNfjuO8cL4bra5s3cM/h5fO02DMCGhYZ1EomOhp49PSD4RkXp+sbnzunV2LNnQ4ECJjdKCJEbSA84tzt5UuduPnZMB+EPP3TjptnMXS0UzEWfchQjFl8S8cnjT4+vQqGYiY1SStc4vusunfHqq6/gscdMr/gkhMhdpAecW1mt+nvJkrpq0erV8MUX5gdfpeC33+DBB5n3l53qoSUpf2EbEx9ehu297NXqdapDh/SwfJ06cPiw7oI//7wEXyGE00kPOLdJTqjx9tuwbBlUrqwXD5ktMlJn1Vq9Glav5t8ijfhf+FlK1CzBjBnQtKnJi6wSE/UCq/fe00H3ww8hWGr9CiFcRwJwbnLwoO6tLVyo5y2VMrtF2urV0Lo1KqlXPiXPK7x6+TOGDLMwZAjm7elNFh+vn6/Nm3Xvd/x4KF/e5EYJIXI7CcC5xRdfwFtv6X28EybgtCr0zrBkCcpqxQASsWAvUYpN/1ioWdPkdiUk6OifJ48OvG+/rcssmr7ySwhxO/CQd2iRY/v3Q5s2sGsX9OtnfvBNSIAvvsB25Tq/X2hHHAFYsaD8/On7S6j5wTc8XA/PR0bq34cNg27dJPgKIdwm0x6wYRjfAp2B00qpWkmXDQOeBc4kHfaWUupvVzVSpCE2FkaMgE6doHlz+OwznR7REwLIhg3wzDOwfTsfjCvCsH970a7cLNqU/oMD9evQMV95uprVtmPH4KWX9Hai2rU9YPxbCHG7cmQI+jtgPPDDLZePVUqNcXqLROaWL9f1aPft08OnzZu7pXBCpq5dg3ffRX3xBVfyleIpy2wWn+tE6Qe2sqeqjb3GgwCsvLXakbtMngwDB+oc2B9/rNNKesLzJoS4LWU6TqmUWgGcd0NbRGYuXNCBNzRUB5FFi/SqXRcKj4qh+ailVBw8j+ajlhIeFZP2gZGReiHT2LFML9SXcld2kffRLlR7cQ3+1Y7d1DG/Ue3I3S5f1slIdu6EN96Q4CuEMFVOJgpfMgxjm2EY3xqGUdhpLRLp++knmDYNBg2C7dvh3ntderrkWr0xF2NRpNTqvSkInzsHf/+NatsW+7YdxJGH3/x68+u8Qvz0E5xJvJzmfadXBcmprl3Tz9WMGfr311+Hv/+GihVdf24hhMhEdgPwROBOoC5wAkh3o6lhGH0Nw9hoGMbGM2fOpHeYSE9MDKxYoX9+4QWdGvGTT3S9XhfLsFavUjB9OtSoweWXhmCLTcAHO35GItNfiKBjR338rdWOkqV3udP8848uszhmjH7OQC9M84Q5ciGEIJsBWCl1SillU0rZga+BxhkcO0Up1VAp1bB48eLZbeftx26HiRN17dnevXWiCF9fXcXITdLrpdqPHIUHHoBHH+VAYjlePPg6VsMfu4+u15snLPTGsYPCqhLoZ7np9oF+FgaFVXVNo0+e1BUcOnbUH1JWrICPPnLNuYQQIgeyFYANwyid6tcHgR3OaY4A9Faili31dqJGjWDJElMKwKfVS21+aAuLv+1H4oLFDM07hpqX13Lne72xLFuCzwcjdFtTpZLsWi+Ykd1CCA4KxACCgwIZ2S3EdQuwVq7URROGD9c9XykbKITwUIbKJFuSYRi/AqHo9PingPeSfq8LKOAQ8JxS6kRmJ2vYsKHauHFjTtqb+0VH6zzEBQrorUW9e5s2bJo8B1z90A6aHtnG2vK1OZbnTj4K/4Hnzo2heJM7mToVatUypXkpoqN1qcBHHtFD48eOQblyJjdKCCHAMIxNSqmGaV6XWQB2JgnAGTh1ShdOUArGjoUnnoASJcxuFSu+mUmzvj3wtduI8wngPr/FbLA056OP9HZaiyXz+3CZuDgYOVJ/lSgB//6rt2UJIYSHyCgASyYss126pIeaK1XS2awMAwYM8Ijgy8aNtHznRfySavX62q30LreCnTvh1VdNDr7Ll+tSge+/Dz166DzOEnyFEF5EArCZwsP1IqvJk6FvXyhVyuwWabGxMHAgqkkTrl1OJB5/rFgw/P353/ehVKhgcvsOHtRpNxMSYP58vT3LEz6wCCFEFkgxBjPY7Xq+8o8/dDrE8HC92OoW4VExjF6wh+MXYykTFMigsKpOX7yU5jmqFSF2xl/8VehZnr3wMf3b7WJQowjydw41r1avUrBpEzRsqPfx/vkntG/vlu1YQgjhChKA3UkpPcTs46OHnEeN0sPNaWRkSl4AlbwPNzkJBjgvhWPqRVY9Dmyi/MWTvHd6AN+cbE3EkU0Elc3Pzz9A584m1+o9dEgP08+fr/NMN2gAXbua1x4hhHACCcDusmuXLhH40Uc6HeLHH2d4eEZJMJwVgEcv2EP1gzv47dch+NkTAVi9vz2fJtxPv356bVPBgk45VfYkJsKXX8K77+oPLmPH6nlfIYTIBSQAu1p8vI5kH32ktxadO+fQzdJLguHMFI62o0f5eP6X+CcF30Qs5LHEU+rxNUyYcLfTzpMtSul53pUrdcWnr76C8uXNbZMQQjiRLMJypVWrdI9t+HC9Ujc6WmeQcoA7Ujh+tmQi5c6fJAE/rFhI8PFjV9eCVKoV57RzZFlsbMpQ/RNPwO+/w19/SfAVQuQ60gN2pfXr9V7Vf/6BDh2ydNNBYVVvmgMGJ6Vw3LsXChXimLUknxT4hu0qgOCi/9Khwh9srF6F3RVqMNJVaSIzM3++Hqb/+GO9SK1vX3PaIYQQbiA9YGdSSlfeCQ/Xv7/yCuzYkeXgCy5I4Wi1wqhRqNq1ie46mBo1YNGumrQaUJTEAb58fe8DnKpZ37VpItNz+jQ8/jjcdx8EBEDZsu49vxBCmEAyYTnLkSPw4oswd64OJH//bXaLUkybBm+/DSdOEFG0Oz3PjSOkXWkmT/aAynwzZuhe75Ur8NZbMGSIJNQQQuQaGWXCkiHonLLZYPx4HeCUgk8/1T1fT/H++6j33gMgAX9GJrzOqO9Km5li+r+qV4evv9bfhRDiNiEBOKcWL4b+/fUw88SJmJ8mKklcHAQEcPxgPKUw8EHha9j44+UICjzp3D29WUoYYrXqDykBAfp5e+gh6N5d740WQojbiLzrZcf167Bsmf65fXv9899/e0bwvXAB+vTB1qo1r/e38fD3nYkn4Eat3gKdQ516uuRkHjEXY1GkJAwJj4r578Hr1+tMVkOGwMaNNycmEUKI24y882XVggW6/l6nTnD2rA4goaGeMZ77559Qowb2ad8xdU9Lxn+RSO3nmmFbmHatXmfIKGHIDVeu6N5u06Z6H/SsWTp/syc8Z0IIYRIZgnbU6dM6beTPP0PVqnprUbFiZrcKIiN173vFClixgsNF6vKgfR7XStZn0Rxo2RKgGbRzTSpJhxKG7Nql58lfeMED0msJIYRnkADsiEuXdK/34kV47z3PWakbGQlt26ISEsBuZ1rAC7x46QsGvOXHu+/qaVZXKxMUSEwaQbimJRa++w6eegqaNNGlFj1hiF4IITyEBOCMnDkDxYtDoUIwbBi0bu05K3X374d+/VDxCRh2G1YsXC9SjrV/+1Gnjvua8Z+EIUrx2K6lDFvxLSTE6znyMmUk+AohxC0kAKclIUFnY/roI1i0SBdP6NfP6afJVrnBxET49FPUsGFY7Rbsdl8sAH7+PP9bKL5uDL6QUplp9II9+B38lzFLJtLw3yj9nE2ZooOvEEKI/5AAfKtVq3QKxOhonb/5zjtdcppslRuMioJnnoGoKFYWfZBHz42ne8PDDGsVQdHuoabV6u1aL5iuVQtD+R56m9HEifo5lNXNQgiRLgnAqQ0YoEve3XEHzJsHHTu67FQOlxuMjISICAgNxf7W21zfd4I+vn+wwNadz76Fp54qg2GYWKt3zx6oUgXy5oVvv9W1eoPdnMpSCCG8kATg5FSchqED7+uv6+pF+fK59LQOrR6OjNTzzomJ2Hz9eb3Uz3x/NZR7HypM9DgoVcqlTczYtWt6QdrYsfDbb/Dww9Cli4kNEkII73J7jxEeOqT38/76q/791VdhzBiXB19woNzgxYs6pWV8PNhs2OMTKHF+N9NmFWbGDJOD78KFEBKiM1o9+yy0a2diY4QQwjvdngHYatWBtmZNvX821nlF7h01KKwqgX6Wmy67UW5w1iyoUQO1aROJ+GLFgrL48/IfoXTt6vam3mzgQAgLAz8/WL4cJk2CoCCTGyWEEN7n9huC3rgR+vSBrVv1kOn48VCunNubkXr18E2roFfNhFde4UiROjyo/uLOsgmM6hBBpadD8TdpkRVK6S8fH7j7br3B+J133LPRWAghcqnbLwAfOaL3986cCV27mpoOsWu9YB2IlYLz51FFijJn86NE5Y9n5MVXeW2wH0OHQmCgiYusDh/WGaxatYI334Ru3fSXEEKIHLk9hqBnzdJbYwAefBD27tXfPSEX8YwZULky1rtb0r2rjQf6FGdOlYFEbvRj5EgITHuq2PVsNvjyy5RhekkfKYQQTpW7e8BHj8LLL8Ps2Tod4nPP6WFUJyyyylYSjdQSE6F/f9SECQAo/Dh/cC2ffNKc114DXzNfmehoePppWLvW88osCiFELpE7A7DNBhMmwNtv658/+URX43FSYghHk2ikG6SPHNE98M2bATAAH+z8+coKig5q7pQ25silS3DggC480bOnZ4wUCCFELpM7h6B37NAB9557YOdOGDRIr9p1EkdK8GVUJ9caVJyj5/Mx1PIBsQTqWr2B/jqblZOFR8XQfNRSKg6eR/NRS9Ou0wuwZg2MGqV/btpUb9F67DEJvkII4SK5JwBfuQK//65/rlMHNm3SZfoqVnT6qRxJonFrkG5yZDuTf3qLUZP/pXGrQMofWs7OB97m+hxdq9dwQa3ejD4E3HD5Mrz0kv6wMnmyfh7BxMlnIYS4PeSOIejZs3UQOX4cGjaESpWgXj2XnS69Enypk2skB+N7DkYxcMUP1D25j4N5ynFpclEuloI//zSSFhM3g/tds8o503SX8+bB889DTIxO+vHBB5A/v0vaIoQQ4mbeHYCPHdOLrMLDdb3e33/XwdfF/lOCj1RJNJKUCQqk8/wfGbz8OwCs+PK/+O8536gwexY6J3dFZgvBMuypnz6ti01UrKhXYjdtmvMGCSGEcJj3DkHHxUGjRrBggZ673LzZbdWAutYLZmS3EIKDAjGA4KBARnYLuSn4DWpfhSei/gH0IitQ3Fd3BpMnK6cF38yGl/+T7lIp7jkYRZlCAVCiBCxZop83Cb5CCOF23tsDDgiAr77S873Z6PXmdBvRjSQaqSkF338PYWHYDpbldfU13/EIfiRg8/Wl6asdaJWVrUoZcKSaUuqeevCl03y0YDytDm4m8p7v9Q0k8AohhGm8NwCD3sqTDdmqxZuZAwd0DdwlS5hebSiP7h5O3br3c/TVJVQ9EYFfaCitnNhDd2QhWNd6wWCzcWj4xzy74FsMw2DrGyNo9uITTmuHEEKI7PHuAJxNDtfizUxkJCxdCidPor75Bqvy442ASUw5+CyjRunywn5+zQDnD407shAMoOtH/WHuLLjvPpg0iTrlyzu9LUIIIbLutgzADtXizUxkJLRtq+eilWJ7wXu47/JvVAkNZusUuOsuJzU2HRkuBEtI0Pt3/fzgf/+D7t1lT68QQngY712ElQOZ1uLNTFwczJ6NSkgApbDhw6z4+xg2JZglS1wffCGDhWDxR6F+fV2rF+D+++HxxyX4CiGEh7kte8CObCNK14oV8OyzxCZaMOz+WEjAbvGn3++tKd7FhY1Ow00Lwa5e1SUCv/wSgoOhdm33NkYIIUSW3JYBON1avBnN/166BIMHw6RJnCtUkceuTMCvcD5GdYig1kuhFDerVi/AypXQq5cuHdivH4wcKdWLhBDCw92WARjS2UaUnp07ISwMdeIE3wYN4JWL79PzmXyMHg2FC5sYeJP5++sKTytX6pSSQgghPN5tG4AdsmYNLF/OlZBm7PVrygv2NzhXpDF//Qlt2pjYLqXgjz9g2zYYMUKXWty+3WnVnoQQQrieBOC0KAVDh8KHH2I3fPC1+/OKsYTQQY0ZNgzy5jWxbcePw4sv6vSbDRvqkosBARJ8hRDCy8i79q0OHoSwMPjgA5RS+Nht+JHA7/0i+OQTE4OvUjB1KtSoAfPn6xrHkZE6+AohhPA6EoCT2WwwdiyqVi0SVq7lS/+BxBKIzdC1eoMfDzW3fcePw6uv6tSb27bpGse+MoAhhBDeSt7BkxkGsb/NZktgG3qc+4pKLcvRtV83yh+IgNBQtxV6uInNpoeau3XTW4vWrdM9YBluFkIIr3d7B+CICPjgAxJffIWx+7vwyda/SPDPz+jJBn36gI+Pa9JIOmTHDnjmGVi/HhYv1lm3atUypy1CCCGc7vYNwJMm6T2zSqGWrmCmWk7zB5oxYYLubJomIUHv4/3wQyhUCH75xeQl10IIIVzh9gvAly/rhBoTJ6JIqtWr7EzpGUGtn5uZn7Gxc2dYtEjnbv78cyhe3OQGCSGEcIXbbzLx3XdRkyYxN3+PG4usfAP9CXk51Lzge/06WK3659deg7/+gp9/luArhBC5WK7tAYdHxdxINVnDN47+jUrSOPQehl98l9XqMc4Ub8L0DyNpdC3CvEVWAMuWQZ8++mvIEF02UAghRK6XKwNweFQMQ2Zup/rBHby1MZwWB6OILliFSsZGLp4rxmuvF2P4cMiXz8RFVpcuwRtvwJQpcOed5n0AEEIIYYpcGYBHL9hDyx0r+Gr2KCxKYcPg87ODiCt1nbVrC9CoUeb3kboH7VCxhqxYtkwXTzhxAgYOhOHDTU6vJYQQwt1yZQAutWMz4+Z8go9SANjxoXqFjUQ+FESjRh0zvX1yDzq5XGHMxViGzNwO4JwgnC8fFCsGM2dC48Y5vz8hhBBeJ3ctwoqPB+B4iQb8nSeMOAKwYiHR4suWe8oQXNSxtI2jF+y5qVYwQKzVxugFe7LXLqXgt9/grbf0740bw+bNEnyFEOI2ljt6wPHx8OGHqN9+Y9xTm9k8sRXd7S24t/FvtAlYxLryIURXqMXIsKoO3d3xi7FZujzjOzsOL7wAc+bogBsXJ8UThBBCZB6ADcP4FugMnFZK1Uq6rAgwHagAHAJ6KKUuuK6ZaYiM1JmsihbV+2Wjo/mnyBMMfTuRsPt96NLvFNO2lGHSxR6UCQpkZBbmcMsEBRKTRrAtExToePuUgm+/hddf1x8QRo+G/v0lf7MQQgjAsR7wd8B44IdUlw0GliilRhmGMTjp9zed37x0REbq1IxxcSiluJK3JD19/mGDpQOTf4MePcAwStOnQ+ls3f2gsKo3zQEDBPpZGORgDxpIKZ7QoIGuYnTXXdlqixBCiNwp03FQpdQK4PwtFz8AfJ/08/dAV+c2KxMRETplo1IoDD67/jzFe3UgOhoeeYQcJ9ToWi+Ykd1CCA4KxACCgwIZ2S0k8x603Q6zZuneb3AwrF2rVzxL8BVCCHGL7I6HllRKnQBQSp0wDKOEE9uUudBQrD7+YEsg0fDnvs/CaNLfuafoWi84ayue9+zRxRNWr4YFC6B9eymeIIQQIl0uXwlkGEZfwzA2Goax8cyZM86502bNiB63hIXNR8DiJTTpb2ISi8RE+PhjXad31y74/nto18689gghhPAKhkraK5vhQYZRAZibahHWHiA0qfdbGohQSmU6QdqwYUO1cePGHDbZw3TponM3d+8O48dDqVJmt0gIIYSHMAxjk1KqYVrXZbcHPAd4MunnJ4HZ2bwf75SQkFI84fnnYcYM+OMPCb5CCCEclmkANgzjVyASqGoYxjHDMJ4BRgHtDMPYB7RL+v32sGGDXtk8erT+vWNHeOghc9skhBDC62S6CEsp1TOdq9o6uS2eLTYW3nsPPv0USpfWc75CCCFENklWCEesXw9PPAH79sGzz+reb6FCZrdKCCGEF/PKAOzSSkVpMQy9t3fxYp0ARAghhMghrwvALq9UlGzxYr2n9733oFEjiI6WNJJCCCGcxusqAji9UtGtLl3Sw8zt2sGvv8LVq/pyCb5CCCGcyOsCsFMrFd1q3jyoWVMXUXjjDYiKgvz5c36/QgghxC28rlvnlEpFaTl7VieSrlBB53Nu1Chn9yeEEEJkwOt6wIPCqhLoZ7npsixXKkpt1Sq9wKpYMViyBDZtkuArhBDC5bwuAGe7UtGtzpzRPd4WLSA8XF/WpAnkyePsJgshhBD/4XVD0JCNSkWpKQXTp8PLL8Ply/DBB9C5s3MbKIQQQmTCKwNwjvTrB5Mm6WHmadP0oishhBDCzW6PAKwU2O1gsUCnTlCpErz2mmwtEkIIYRqvmwPOspgYXTJwVFK9iM6dYdAgCb5CCCFMlXsDsFLwzTdQo4Ze3RwUZHaLhBBCiBtyZzfwyBGdzWrhQmjVSgfiO+80u1VCCCHEDbmzB3zqFKxbBxMmwNKlEnyFEEJ4nNzTAz54EObO1duLGjWCo0ehQAGzWyWEEEKkyft7wHY7jB8PISHwzju69wsSfIUQQng07w7A+/dD69a613vPPbB9O5QsaXarhBBCiEx57xB0bCw0bw7x8bp60VNPgWGY3SohhBDCId4bgAMDdSarOnUgOJtpKYUQQgiTeG8ABujY0ewWCCGEENni3XPAQgghhJeSACyEEEKYQAKwEEIIYQIJwEIIIYQJJAALIYQQJpAALIQQQphAArAQQghhAgnAQgghhAkkAAshhBAmkAAshBBCmEACsBBCCGECCcBCCCGECSQACyGEECYwlFLuO5lhnAEOO/EuiwFnnXh/ZpLH4nlyy+MAeSyeKrc8ltzyOMD5j+UOpVTxtK5wawB2NsMwNiqlGprdDmeQx+J5csvjAHksniq3PJbc8jjAvY9FhqCFEEIIE0gAFkIIIUzg7QF4itkNcCJ5LJ4ntzwOkMfiqXLLY8ktjwPc+Fi8eg5YCCGE8Fbe3gMWQgghvJJXBGDDMDoYhrHHMIz9hmEMTuN6wzCML5Ou32YYRn0z2pkZwzDKGYaxzDCMaMMwdhqG8Woax4QahnHJMIwtSV9DzWirIwzDOGQYxvakdm5M43qPf10Mw6ia6rneYhjGZcMw+t9yjMe+JoZhfGsYxmnDMHakuqyIYRiLDMPYl/S9cDq3zfD/yt3SeSyjDcPYnfT3M8swjKB0bpvh36K7pfNYhhmGEZPq76hjOrf1mNclnccxPdVjOGQYxpZ0butpr0ma77+m/r8opTz6C7AA/wKVAH9gK1DjlmM6Av8ABtAUWGd2u9N5LKWB+kk/FwD2pvFYQoG5ZrfVwcdzCCiWwfVe8bqkaq8FOInet+cVrwnQEqgP7Eh12SfA4KSfBwMfp/NYM/y/8pDH0h7wTfr547QeS9J1Gf4teshjGQYMzOR2HvW6pPU4brn+U2Col7wmab7/mvn/4g094MbAfqXUAaVUAvAb8MAtxzwA/KC0tUCQYRil3d3QzCilTiilNif9fAWIBoLNbZVLecXrkkpb4F+llDOTxbiUUmoFcP6Wix8Avk/6+Xugaxo3deT/yq3SeixKqYVKqcSkX9cCZd3esGxI53VxhEe9Lhk9DsMwDKAH8KtbG5VNGbz/mvb/4g0BOBg4mur3Y/w3aDlyjEcxDKMCUA9Yl8bVzQzD2GoYxj+GYdR0b8uyRAELDcPYZBhG3zSu97bX5VHSfzPxltcEoKRS6gToNx2gRBrHeNtrA/A0ekQlLZn9LXqKl5KG079NZ6jTm16XFsAppdS+dK732Nfklvdf0/5fvCEAG2lcduvSbUeO8RiGYeQH/gT6K6Uu33L1ZvQQaB1gHBDu5uZlRXOlVH3gPuBFwzBa3nK917wuhmH4A12AGWlc7U2viaO85rUBMAzjbSAR+DmdQzL7W/QEE4E7gbrACfTw7a286XXpSca9X498TTJ5/033ZmlcluPXxRsC8DGgXKrfywLHs3GMRzAMww/94v+slJp56/VKqctKqatJP/8N+BmGUczNzXSIUup40vfTwCz0ME1qXvO6oN8kNiulTt16hTe9JklOJQ/1J30/ncYxXvPaGIbxJNAZeFwlTcjdyoG/RdMppU4ppWxKKTvwNWm30SteF8MwfIFuwPT0jvHE1ySd91/T/l+8IQBvAO4yDKNiUi/lUWDOLcfMAXonrbptClxKHlLwJElzJt8A0Uqpz9I5plTScRiG0Rj9Gp1zXysdYxhGPsMwCiT/jF4ss+OWw7zidUmS7qd5b3lNUpkDPJn085PA7DSOceT/ynSGYXQA3gS6KKWup3OMI3+Lprtl/cODpN1Gr3hdgHuB3UqpY2ld6YmvSQbvv+b9v5i9Ms2RL/Rq2r3oVWhvJ132PPB80s8GMCHp+u1AQ7PbnM7juAc9bLEN2JL01fGWx/ISsBO9ym4tcLfZ7U7nsVRKauPWpPZ68+uSFx1QC6W6zCteE/SHhhOAFf0p/RmgKLAE2Jf0vUjSsWWAv1Pd9j//Vx74WPaj596S/18m3fpY0vtb9MDH8mPS/8E29Jt3aU9/XdJ6HEmXf5f8/5HqWE9/TdJ7/zXt/0UyYQkhhBAm8IYhaCGEECLXkQAshBBCmEACsBBCCGECCcBCCCGECSQACyGEECaQACyEEEKYQAKwEEIIYQIJwEIIIYQJ/g+Un5D1BjmJ5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(x, y, 'o', label=\"Data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "ax.plot(x, res2.fittedvalues, 'r--.', label=\"Predicted\")\n",
    "ax.plot(x, iv_u, 'r--')\n",
    "ax.plot(x, iv_l, 'r--')\n",
    "legend = ax.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joint hypothesis test\n",
    "\n",
    "### F test\n",
    "\n",
    "We want to test the hypothesis that both coefficients on the dummy variables are equal to zero, that is, $R \\times \\beta = 0$. An F test leads us to strongly reject the null hypothesis of identical constant in the 3 groups:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.638159Z",
     "iopub.status.busy": "2021-02-02T06:52:02.636664Z",
     "iopub.status.idle": "2021-02-02T06:52:02.647455Z",
     "shell.execute_reply": "2021-02-02T06:52:02.648549Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 0 0]\n",
      " [0 0 1 0]]\n",
      "<F test: F=array([[145.49268198]]), p=1.2834419617289997e-20, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "R = [[0, 1, 0, 0], [0, 0, 1, 0]]\n",
    "print(np.array(R))\n",
    "print(res2.f_test(R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also use formula-like syntax to test hypotheses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.653610Z",
     "iopub.status.busy": "2021-02-02T06:52:02.652089Z",
     "iopub.status.idle": "2021-02-02T06:52:02.662492Z",
     "shell.execute_reply": "2021-02-02T06:52:02.663667Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[145.49268198]]), p=1.2834419617290716e-20, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res2.f_test(\"x2 = x3 = 0\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Small group effects\n",
    "\n",
    "If we generate artificial data with smaller group effects, the T test can no longer reject the Null hypothesis: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.668675Z",
     "iopub.status.busy": "2021-02-02T06:52:02.667193Z",
     "iopub.status.idle": "2021-02-02T06:52:02.675061Z",
     "shell.execute_reply": "2021-02-02T06:52:02.676150Z"
    }
   },
   "outputs": [],
   "source": [
    "beta = [1., 0.3, -0.0, 10]\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + np.random.normal(size=nsample)\n",
    "\n",
    "res3 = sm.OLS(y, X).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.681190Z",
     "iopub.status.busy": "2021-02-02T06:52:02.679639Z",
     "iopub.status.idle": "2021-02-02T06:52:02.688586Z",
     "shell.execute_reply": "2021-02-02T06:52:02.689742Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[1.22491119]]), p=0.30318644106321035, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res3.f_test(R))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.694363Z",
     "iopub.status.busy": "2021-02-02T06:52:02.692922Z",
     "iopub.status.idle": "2021-02-02T06:52:02.703048Z",
     "shell.execute_reply": "2021-02-02T06:52:02.704135Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[1.22491119]]), p=0.30318644106321035, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res3.f_test(\"x2 = x3 = 0\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multicollinearity\n",
    "\n",
    "The Longley dataset is well known to have high multicollinearity. That is, the exogenous predictors are highly correlated. This is problematic because it can affect the stability of our coefficient estimates as we make minor changes to model specification. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.708835Z",
     "iopub.status.busy": "2021-02-02T06:52:02.707250Z",
     "iopub.status.idle": "2021-02-02T06:52:02.730856Z",
     "shell.execute_reply": "2021-02-02T06:52:02.732164Z"
    }
   },
   "outputs": [],
   "source": [
    "from statsmodels.datasets.longley import load_pandas\n",
    "y = load_pandas().endog\n",
    "X = load_pandas().exog\n",
    "X = sm.add_constant(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit and summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.737688Z",
     "iopub.status.busy": "2021-02-02T06:52:02.736136Z",
     "iopub.status.idle": "2021-02-02T06:52:02.763170Z",
     "shell.execute_reply": "2021-02-02T06:52:02.764205Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                 TOTEMP   R-squared:                       0.995\n",
      "Model:                            OLS   Adj. R-squared:                  0.992\n",
      "Method:                 Least Squares   F-statistic:                     330.3\n",
      "Date:                Tue, 02 Feb 2021   Prob (F-statistic):           4.98e-10\n",
      "Time:                        06:52:02   Log-Likelihood:                -109.62\n",
      "No. Observations:                  16   AIC:                             233.2\n",
      "Df Residuals:                       9   BIC:                             238.6\n",
      "Df Model:                           6                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const      -3.482e+06    8.9e+05     -3.911      0.004    -5.5e+06   -1.47e+06\n",
      "GNPDEFL       15.0619     84.915      0.177      0.863    -177.029     207.153\n",
      "GNP           -0.0358      0.033     -1.070      0.313      -0.112       0.040\n",
      "UNEMP         -2.0202      0.488     -4.136      0.003      -3.125      -0.915\n",
      "ARMED         -1.0332      0.214     -4.822      0.001      -1.518      -0.549\n",
      "POP           -0.0511      0.226     -0.226      0.826      -0.563       0.460\n",
      "YEAR        1829.1515    455.478      4.016      0.003     798.788    2859.515\n",
      "==============================================================================\n",
      "Omnibus:                        0.749   Durbin-Watson:                   2.559\n",
      "Prob(Omnibus):                  0.688   Jarque-Bera (JB):                0.684\n",
      "Skew:                           0.420   Prob(JB):                        0.710\n",
      "Kurtosis:                       2.434   Cond. No.                     4.86e+09\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The condition number is large, 4.86e+09. This might indicate that there are\n",
      "strong multicollinearity or other numerical problems.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/travis/miniconda/envs/statsmodels-test/lib/python3.7/site-packages/scipy/stats/stats.py:1604: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    }
   ],
   "source": [
    "ols_model = sm.OLS(y, X)\n",
    "ols_results = ols_model.fit()\n",
    "print(ols_results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Condition number\n",
    "\n",
    "One way to assess multicollinearity is to compute the condition number. Values over 20 are worrisome (see Greene 4.9). The first step is to normalize the independent variables to have unit length: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.769047Z",
     "iopub.status.busy": "2021-02-02T06:52:02.767587Z",
     "iopub.status.idle": "2021-02-02T06:52:02.778455Z",
     "shell.execute_reply": "2021-02-02T06:52:02.779604Z"
    }
   },
   "outputs": [],
   "source": [
    "norm_x = X.values\n",
    "for i, name in enumerate(X):\n",
    "    if name == \"const\":\n",
    "        continue\n",
    "    norm_x[:,i] = X[name]/np.linalg.norm(X[name])\n",
    "norm_xtx = np.dot(norm_x.T,norm_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we take the square root of the ratio of the biggest to the smallest eigen values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.784211Z",
     "iopub.status.busy": "2021-02-02T06:52:02.782919Z",
     "iopub.status.idle": "2021-02-02T06:52:02.791448Z",
     "shell.execute_reply": "2021-02-02T06:52:02.792634Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56240.86912789531\n"
     ]
    }
   ],
   "source": [
    "eigs = np.linalg.eigvals(norm_xtx)\n",
    "condition_number = np.sqrt(eigs.max() / eigs.min())\n",
    "print(condition_number)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dropping an observation\n",
    "\n",
    "Greene also points out that dropping a single observation can have a dramatic effect on the coefficient estimates: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.797161Z",
     "iopub.status.busy": "2021-02-02T06:52:02.795912Z",
     "iopub.status.idle": "2021-02-02T06:52:02.807700Z",
     "shell.execute_reply": "2021-02-02T06:52:02.808728Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage change 4.55%\n",
      "Percentage change -2228.01%\n",
      "Percentage change 154304695.31%\n",
      "Percentage change 1366329.02%\n",
      "Percentage change 1112549.36%\n",
      "Percentage change 92708715.91%\n",
      "Percentage change 817944.26%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ols_results2 = sm.OLS(y.iloc[:14], X.iloc[:14]).fit()\n",
    "print(\"Percentage change %4.2f%%\\n\"*7 % tuple([i for i in (ols_results2.params - ols_results.params)/ols_results.params*100]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at formal statistics for this such as the DFBETAS -- a standardized measure of how much each coefficient changes when that observation is left out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.813353Z",
     "iopub.status.busy": "2021-02-02T06:52:02.811909Z",
     "iopub.status.idle": "2021-02-02T06:52:02.818826Z",
     "shell.execute_reply": "2021-02-02T06:52:02.819830Z"
    }
   },
   "outputs": [],
   "source": [
    "infl = ols_results.get_influence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general we may consider DBETAS in absolute value greater than $2/\\sqrt{N}$ to be influential observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.824596Z",
     "iopub.status.busy": "2021-02-02T06:52:02.822947Z",
     "iopub.status.idle": "2021-02-02T06:52:02.831665Z",
     "shell.execute_reply": "2021-02-02T06:52:02.832725Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2./len(X)**.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:52:02.837266Z",
     "iopub.status.busy": "2021-02-02T06:52:02.835657Z",
     "iopub.status.idle": "2021-02-02T06:52:02.871484Z",
     "shell.execute_reply": "2021-02-02T06:52:02.872631Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    dfb_const  dfb_GNPDEFL       dfb_GNP     dfb_UNEMP     dfb_ARMED  \\\n",
      "0   -0.016406  -169.822675  1.673981e+06  54490.318088  51447.824036   \n",
      "1   -0.020608  -187.251727  1.829990e+06  54495.312977  52659.808664   \n",
      "2   -0.008382   -65.417834  1.587601e+06  52002.330476  49078.352378   \n",
      "3    0.018093   288.503914  1.155359e+06  56211.331922  60350.723082   \n",
      "4    1.871260  -171.109595  4.498197e+06  82532.785818  71034.429294   \n",
      "5   -0.321373  -104.123822  1.398891e+06  52559.760056  47486.527649   \n",
      "6    0.315945  -169.413317  2.364827e+06  59754.651394  50371.817827   \n",
      "7    0.015816   -69.343793  1.641243e+06  51849.056936  48628.749338   \n",
      "8   -0.004019   -86.903523  1.649443e+06  52023.265116  49114.178265   \n",
      "9   -1.018242  -201.315802  1.371257e+06  56432.027292  53997.742487   \n",
      "10   0.030947   -78.359439  1.658753e+06  52254.848135  49341.055289   \n",
      "11   0.005987  -100.926843  1.662425e+06  51744.606934  48968.560299   \n",
      "12  -0.135883   -32.093127  1.245487e+06  50203.467593  51148.376274   \n",
      "13   0.032736   -78.513866  1.648417e+06  52509.194459  50212.844641   \n",
      "14   0.305868   -16.833121  1.829996e+06  60975.868083  58263.878679   \n",
      "15  -0.538323   102.027105  1.344844e+06  54721.897640  49660.474568   \n",
      "\n",
      "          dfb_POP      dfb_YEAR  \n",
      "0   207954.113588 -31969.158503  \n",
      "1    25343.938290 -29760.155888  \n",
      "2   107465.770565 -29593.195253  \n",
      "3   456190.215133 -36213.129569  \n",
      "4  -389122.401700 -49905.782854  \n",
      "5   144354.586054 -28985.057609  \n",
      "6  -107413.074918 -32984.462465  \n",
      "7    92843.959345 -29724.975873  \n",
      "8    83931.635336 -29563.619222  \n",
      "9    18392.575057 -29203.217108  \n",
      "10   93617.648517 -29846.022426  \n",
      "11   95414.217290 -29690.904188  \n",
      "12  258559.048569 -29296.334617  \n",
      "13  104434.061226 -30025.564763  \n",
      "14  275103.677859 -36060.612522  \n",
      "15 -110176.960671 -28053.834556  \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/stats/outliers_influence.py:693: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return self.resid / sigma / np.sqrt(1 - hii)\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/stats/outliers_influence.py:733: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_internal * np.sqrt(hii / (1 - hii))\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/stats/outliers_influence.py:762: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_external * np.sqrt(hii / (1 - hii))\n"
     ]
    }
   ],
   "source": [
    "print(infl.summary_frame().filter(regex=\"dfb\"))"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}