{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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": {},
   "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": {},
   "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": {},
   "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:                Thu, 05 Nov 2020   Prob (F-statistic):          2.83e-239\n",
      "Time:                        07:28:38   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": {},
   "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": {},
   "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": {},
   "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:                Thu, 05 Nov 2020   Prob (F-statistic):           6.30e-27\n",
      "Time:                        07:28:38   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": {},
   "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": {},
   "outputs": [
    {
     "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": {},
   "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": {},
   "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": {},
   "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:                Thu, 05 Nov 2020   Prob (F-statistic):           5.69e-38\n",
      "Time:                        07:28:38   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": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFlCAYAAAAzqTv+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABrfUlEQVR4nO3dd3iT1RfA8e/btKVllj0KCIjssmSKYAGhCIgIiqKC/BRRcSMoOBiigoKiAoKI4kIFBAqCsimzgEDZZckue5TZkSb398dtoWBH2iZ5k3I+z9OnbfIm703S5uSucwylFEIIIYRwLx+zGyCEEELcjiQACyGEECaQACyEEEKYQAKwEEIIYQIJwEIIIYQJJAALIYQQJvB158mKFSumKlSo4M5TCiGEEKbZtGnTWaVU8bSuc2sArlChAhs3bnTnKYUQQgjTGIZxOL3rZAhaCCGEMIEEYCGEEMIEEoCFEEIIE7h1DjgtVquVY8eOER8fb3ZTvFpAQABly5bFz8/P7KYIIYRwgOkB+NixYxQoUIAKFSpgGIbZzfFKSinOnTvHsWPHqFixotnNEUII4QDTh6Dj4+MpWrSoBN8cMAyDokWLyiiCEEJ4EdMDMCDB1wnkORRCCO/iEQHYbBaLhbp161KzZk3q1KnD559/jt1uz/A2hw4d4tdff3VTC4UQQuQ2ps8BZ1V4VAyjFu7heGwcZYICGRBWlc71gnN0n4GBgWzZsgWA06dP88QTT3Dx4kWGDRuW7m1SAvATTzyRo3MLIYS4PXlVDzg8KoZBs7YTExuHAmJi4xg0azvhUTFOO0eJEiWYNGkS48aNQynFoUOHaN68OfXr16d+/fqsXbsWgIEDB7Jq1Srq1q3LmDFj0j1OCCGESItX9YBHLdxDnNV202VxVhujFu7JcS84tUqVKmG32zl9+jQlSpRg8eLFBAQEsG/fPrp3787GjRsZOXIko0ePZt68eQBcu3YtzeOEEEKItHhVAD4eG5ely3NCKQXofcovv/wyW7ZswWKxsHfv3jSPd/Q4IYRwqsuX4exZkC2IXserAnCZoEBi0gi2ZYICnXqeAwcOYLFYKFGiBMOGDaNkyZJs3boVu91OQEBAmrcZM2aMQ8cJIYTTzJ8PHTuCYUBiIvh61Vv6bc+r5oAHhFUl0M9y02WBfhYGhFV12jnOnDnDCy+8wMsvv4xhGFy8eJHSpUvj4+PDzz//jM2mh8ALFCjA5cuXr98uveOEEMLpLlyAp5/WwRdg8GBIHrUT3sOrPi6lzPM6exV0XFwcdevWxWq14uvrS48ePejXrx8Affv2pWvXrsyYMYOWLVuSL18+AGrXro2vry916tShV69e6R4nhBBOZbfDvffCnj3w3nv6K08es1slssFQbvzU1KBBA3XrwqTo6GiqV6/utjbkZvJcCpGLXbgAhQqBjw/8+ScEB0P9+hAfDxERULWqzAN7IMMwNimlGqR1nVcNQQshxG1p5kyoVg0mTtS/P/igDr4A167BAw9AeLhpzRPZIwFYCCE81enT0K0bPPIIlC2rh55vVbgw5M0Lx465v30iRyQACyGEJ/rzT6hRA+bMgQ8/hHXroHbt/x5nGDo4SwD2Ol61CEsIIW4befNC5crw3XdQs2bGx5YrB0ePuqddwmky7QEbhhFgGMYGwzC2Goax0zCMYcmXFzEMY7FhGPuSvxd2fXOFECKXUgp++gk+/lj/3ro1REZmHnxBesBeypEh6ASglVKqDlAXaGcYRhNgILBUKXUXsDT5dyGEEFl17Jje0/v007BoESQl6csdLTM6cCDMneu69gmXyHQIWul9SleSf/VL/lLAQ0Bo8uU/AhHA205voYudO3eO1q1bA3Dy5EksFgvFixcHYMOGDfj7+5vZPCFEbqYUfP899Oung+6XX8JLL4HFkvltU6tWzTXtEy7l0BywYRgWYBNQGRivlFpvGEZJpdQJAKXUCcMwSriwnS5TtGjR66UIhw4dSv78+enfv//165OSkvCV9G5CCFc4eBBefBGaNYPJk+HOO7N3P2fO6MVabdtC+fLObaNwGYcii1LKBtQ1DCMImG0YRi1HT2AYRh+gD0B5L/nD6NWrF0WKFCEqKor69etToECBmwJzrVq1mDdvHhUqVOCXX37hq6++IjExkcaNG/P1119jyeqnVyHE7cNuh8WLISwMKlXSq5vr1tUJNrLrxAl47jmYMUMCsBfJUtdOKRVrGEYE0A44ZRhG6eTeb2ngdDq3mQRMAp0JK6P7f/11SO6MOk3duvDFF1m/3d69e1myZAkWi4WhQ4emeUx0dDTTpk1jzZo1+Pn50bdvX6ZOnUrPnj1z0mQhRG61fz88+yysXAmrVul9vSkJNdIQHhXjWOrdcuX0d1kJ7VUyDcCGYRQHrMnBNxC4H/gEmAs8DYxM/j7HlQ11t0cffTTTnuzSpUvZtGkTDRs2BHRO6RIlvHIkXgjhSjabnt997z3w99fzvs2aZXiT8KgYBs3afr0GekxsHINmbQf4bxAOCpJkHF7IkR5waeDH5HlgH2C6UmqeYRiRwHTDMJ4FjgCP5rQx2empukrqYgq+vr7Y7fbrv8fHxwO6ZvDTTz/NiBEj3N4+IYQX6dQJ/vpLr3SeOFHncc7EqIV7rgffFHFWG6MW7vlvADYM2QvshTKddFBKbVNK1VNK1VZK1VJKfZB8+TmlVGul1F3J38+7vrnmqFChAps3bwZg8+bNHDx4EIDWrVvzxx9/cPq0Hn0/f/48hw8fNq2dQggPkpSk53sBevaEX37RW4UcCL4Ax9OofZ7R5bIX2PtIKkoHdO3alfPnz1O3bl0mTJhAlSpVAKhRowYffvghbdu2pXbt2rRp04YTJ06Y3FohhOm2bIFGjWDCBP37Y4/Bk086vq8XKBMUmKXL+f572QvsZWR/TSrpLbYKDAxk0aJFaV732GOP8dhjj7mwVUIIr5GQAB99BCNGQNGiuleaTQPCqt40BwwQ6GdhQFjVtG8gq5+9jvSAhRDCGTZu1Cuahw+H7t1h50546KFs313nesGM6BJCcFAgBhAcFMiILiFpr4IG2L0bhgyBs2ezfU7hXtIDFkIIZ7h0SX/NmwcdOjjlLjvXC04/4N7qwAH44ANdG7hYMaecX7iW9ICFECK7Vq26sX2jVSu9z9dJwTfLUoa7ZSW015AALIQQWXX5Mrz8MrRoAV9/DXHJK5Pz5DGvTSnJOGQltNeQACyEEFmxaBHUqqUD72uvQVQUBKazMtmdJBmH15E5YCGEcNSJE/DggzqH8+rVcM89ZrdI1wyOiIDQUN0LPn7c7BYJB0kPGLBYLNStW5datWrx6KOPcu3atWzfV69evfjjjz8A6N27N7t27Ur32IiICNauXZvlc1SoUIGzstJRCPfZsEF/L10aFizQvV4PCL7XFq7Cdu992N95F9W6NYwdC7/+anazhIMkAKP3+W7ZsoUdO3bg7+/PxIkTb7reZrOlc8uMTZ48mRo1aqR7fXYDsBDCTU6f1kk0GjfWQ88ALVtCQIC57QLWjtvMpQ6PYbFb8UFBYqLeCpWFZB8ilbg4+OYbXaPZTbwzAEdG6o3ukZFOv+vmzZuzf/9+IiIiaNmyJU888QQhISHYbDYGDBhAw4YNqV27Nt988w2g80G//PLL1KhRgw4dOlxPSwkQGhrKxo0bAViwYAH169enTp06tG7dmkOHDjFx4kTGjBlD3bp1WbVqFWfOnKFr1640bNiQhg0bsmbNGgDOnTtH27ZtqVevHs8//zzKjX8gQtyWlIKpU6FGDQgPhw8/1IHXA5w9ZWN+zbdo9Eoj8qgEknz8sWJB+fpD/vy62lJSktnN9D579uiFdatXu++cSim3fd19993qVrt27br5gvvu++/X+PH6uqtXlapbVykfH6VAf69bV6kpU/T1Z87897YOyJcvn1JKKavVqjp16qS+/vprtXz5cpU3b1514MABpZRS33zzjRo+fLhSSqn4+Hh19913qwMHDqiZM2eq+++/XyUlJamYmBhVqFAhNWPGjOSHcp/6559/1OnTp1XZsmWv39e5c+eUUkoNGTJEjRo16no7unfvrlatWqWUUurw4cOqWrVqSimlXnnlFTVs2DCllFLz5s1TgDpz5kzmz6UQInt69NDvMY0bK7Vzp9mtUUopZbcr9euvShUvZlezjIfVP/V6q7jj59X2SWvVQD5Wq0etVeqbb3S7jx41u7ne4eJFpX7++cbvye/RzgRsVOnERO9bhHXx4o0E53a7/j2H4uLiqFu3LqB7wM8++yxr166lUaNGVKxYEYBFixaxbdu26/O7Fy9eZN++faxcuZLu3btjsVgoU6YMrVq1+s/9r1u3jhYtWly/ryJFiqTZjiVLltw0Z3zp0iUuX77MypUrmTVrFgAdOnSgcOHCOX7MQohbpLyv+PjoZBb168Mrr0AmZUldLjKSi9P+Zufv2xl86lMqNrqLyoumE1JPv30XbXIXnejHld2VoEvyXuBjx3KUBvO2sGAB9OkDMTF6iuGuuyD5PdpdPC8AR0Skf13evHpYqHVrPd/h769/b9pUX1+sWMa3T0fKHPCtUpckVEoxduxYwsLCbjrmr7/+wshkzkUplekxAHa7ncjISALT2NLgyO2FENn077/Quzd07aqHIbt3N7tFANjXRKLuC6WgLZGmwNTW1bl74cdYLDfeuktULkhJ1rN8TzSU66ovPHoUmjQxp9Ge7vx5eOMN+OknqF4d1qzRwdcE3jcH3LQpLF2q860uXXoj+LpYWFgYEyZMwGq1ArB3716uXr1KixYt+P3337HZbJw4cYLly5en0eSmrFix4noZw/PndeXGAgUKcPny5evHtW3blnHjxl3/PeVDQYsWLZg6dSoAf//9NxcuXHDJYxTitmOzweefQ0gIbN4MBQqY3aLr9q08waE2z2GxJWIA+Fho1LrAfzrklkB/zvqUxHIyVa9X9gKnLSlJfzD59Vd4/329mt3EDyqe1wN2RNOmbgu8KXr37s2hQ4eoX78+SimKFy9OeHg4Dz/8MMuWLSMkJIQqVapw3333/ee2xYsXZ9KkSXTp0gW73U6JEiVYvHgxDz74II888ghz5sxh7NixfPXVV7z00kvUrl2bpKQkWrRowcSJExkyZAjdu3enfv363HfffZSXqidC5NyuXfDMM7B+PXTsqEsHesCwbWIifPIJBA79glfs+7Bb/DCwY/j7672+aTiXtyz5zh3VyTiKFoWrV93aZo935oweIfX1hZEj4c47oU4ds1uFody4orZBgwYqZVVwiujoaKpXr+62NuRm8lwKkQXLl+stRl9+CY8/7hHbd7bO3M/HAy8xfX99ena5wmdvHqeYce5Goo10Oh4bynWh8Kk93JW4U6/g9oDH4hGUgh9+gH79YPRovULczQzD2KSUapDWdd7ZAxZCiOzYsAHWrYNXX9Xbig4ehFRrPcxy9WISyx4cw/2rBjPQrzZPhq+j00P5gSr6gExG/E7feQ8HjwdQWcl6kesOHtSLrJYsgebN9ZeH8b45YCGEyKpr1+DNN3Ug+/zzG0O0HhB8t7/1I5eL3sGDq95ib4Uw7tw6i04PZS2IHuzan8ftv3LmDPD99x6ziMw0P/6o83WvW6dzdkdEQJUqZrfqPyQACyFyt+XL9SKrzz+H556DrVs9IvCeOwc/NJlIrVG9KGk7jt3XnzpT36JgdQfr/6aSUgjp6FF0z2/69Ns7GUexYnDffbBzJ7z4ot5a5oE8olXunIfOreQ5FCINp07pPb0+PjoQT5wIhQqZ2iSlYPZ356lRA/ZvOIfCwADstiR2/vZntu6zStxWjlKWuLmL9UIyux1OnnRuwz1ZYqLeGfPRR/r3Dh1g/nzw8AWrpgfggIAAzp07JwEkB5RSnDt3jgAPyE8rhEf45x/9vWRJ/Ua8dWu6K4jdKWZnLAsrPE+z3tW4I+g4ux4sRoKvH0mGD1aLLx9eKUF4VEyW77fEXYUoSwwJew/fWMl99KiTW++h/vkHGjSAwYNh9+4buZy9YC7c9EVYZcuW5dixY5w5c8bspni1gIAAynrAFgohTHX6tF5gNW2aznQUFqYT95jMbodFfcOpM6kvbdQpokL7ka/FVjYnlOXJgh/R5Mh21pUPYXPJKhxZuIfO9bI2DF00pAx2DGyHj0G5RvrC3L4X+No1HXTHjIFSpWDOHOjUyexWZYnpAdjPz+96ikYhhMgWpeCXX+D11+HKFfjgA88onhAZyZnfl3D0x+W0u7icf/PX4cQvf9Lgobs5NHA+AJuDq7M5+Mb2weOxcVk+jZHHnzOWkviePKYnhO+8061VfUyxbx989ZXOYPbpp5lOLYRHxTBq4R6Ox8ZRJiiQAWFVs/xBx9lMD8BCCJFjPXvqANy0KUyerKsYZcLVb8hJqyJRrVtT2JpIYSA69AWqLfwKw98PgDJBgcSkEWzLBP03Fa0jzuUtR74LR6FwYdi/PydN91yxsbqn+/TTOpHGvn1wxx2Z3iw8KoZBs7YTZ9WlZWNi4xg0azuAqUHY9DlgIYTIFrv9RgGFDh10b2jVKoeD76BZ24mJjUNx4w05O/Ovadke/i9H7/8fhjUBX2xYLFC9bfnrwRdgQFhVAv1uzisZ6GdhQFjV7J3zzodZafe8va5OM2eOfm2fffbGBwwHgi/AqIV7rgffFHFWG6MW7nF2K7NEArAQwvvs3g0tWsD48fr3xx/PUuUiV70hX7uUxJ/3jebOh0MokXgUw9cXLJY000h2rhfMiC4hBAcFYgDBQYGM6BKS7R7Z1vaDGHT1PWw24O23c89e4FOnoFs36NwZihfXe3srV87SXaQ3rJ+d4X5nkiFoIYT3sFr1fN8HH+i9vEWLZutuXPGGvH7SVgJf7c2DCRvZekcnKv79NZbYIxmmkexcL9hpQ6DlygG2JE7E+FD21Ck9GuDtrFZdLOH4cfjwQ3jrLfDzy/x2t3D2cL+zSAAWQniHqCjo1Qu2bdM9oq++0tuMssEZb8jhUTH8NWkWd27bzMpzD9N6TwQ9LUfYOXQ6dQY/krwNJtilhWNSz2N3XL+XBPqzfdU2ypYrBydO6GQcvl74Nh8TA2XK6GD75Zc6i1W1atm+uwFhVW+aA4acDfc7iwxBCyG8w8WLupbrnDl6m1E2gy/kfP41PCqG37+czthJb/Lm2p+YsedxVtSoz6r5S6k55FG37EG9dR77aF4/LNj5Z1m09ybjsNn0tqIqVWDKFH1Zp045Cr7g/OF+Z/HCj0ZCiNvGsmW6Tm///noYd/9+yJMnx3eb8sab3VXQY37azqjffyaPXdcHV0YCzYotZnRUGR4Oy3HzHHLrPPapkgUB2P/PHuhaT1949KhHlFh0yI4dekvR+vXQvj3cf79T796Zw/3OIgFYCOF5LlyAAQPgu+907+fllyEgwCnBN0V23pDtdlj8ylymft2X0hzHalgwUCRZfFlXPsSti3puPdfpIoWwYxB09rTeBxwa6rE5kP9j7FhdLKNQIZg6VS8g84JMVjklAVgI4VlmztQB98wZvehm6FAdfE22dy9M6TqPETseYod/Dfq0f5eA/FdvZLEKrk6wGxf13DqPbfP15aRRknJxZ6BqVZ372tOl1C6uXFnP648Zo1c6myBlPj3mXDzBRQPckqhDArAQwnMcOwZPPKH3e86fD/Xrm90irImKye8d4o2vKpI3zwN06DmJs33v5+A8PQScksXK3Yt6UhYW+V2+SJFrFzlUJJgJhZ/jUt6qPOW2VmTT5cvwzjt6FfvQobpgxgMPmNac8KgY3vplN2VnX6HnxfVsf6AIg64mAq5N1CEBWAhhLrtdF01v21bPV0ZE6OT62dhu4mz/fvIHfoMH0j3xDOs7HmDEpKKULv2cvtLf39TUhp3rlqHsX7No8OGrALT4cBHrWrzNjnX5+BJ0cpLixeGHH9zWJof89Re88IL+sPXGG2a3BqXgrU9iKR7uy7yE7viTgHWaL08+/hGjFvq79jVVSrnt6+6771ZCCHHdvn1KhYYqBUotX252a667eilJrar3irKDsoOyWfyUWr3a7GbdsG+fUvffr5+3AgWU+uQTpaxWNex9qyrHEZWQoJRq3VqpJk3MbukNp08r9cQTus3Vqyu1dq3ZLVJHjijVoYNu0vi8vZVdx2NlNXzUJy16qgpvz8vxOYCNKp2Y6CUz9EKIXCUpSSfUCAnR+3u//VYXUE8WHhVDs5HLqDhwPs1GLnNaikhHrJh3mX3Fm3Jv1FgADMAHO6xc6bY2ZMhq1YUmNmzQmcAuXNBz5b6+tNv9BUcoz/Hdl3RmDk8qSRgTA+HhMGSIfs1duD86M3a7fupq1IDVyxKp0GEvyzvdic3wuV4acl35EJcn6pAhaCGE+3XoAIsWwcMPw7hxOulCMrMS5184rxjwlsF33xVgasG7KdC9I5WmjdTF3tNIJel2GzdCvXp6aP6nn/RCqzJlIC5OD+HXqEFAZb3l6OyWY1QoW9b8ZByHD+ug+9prULcuHDmS7exlzhIdDc89B+vXWJlY+TN6Wiez8K0/GbC4Lt2e/OT6orroCrUY4eI5fekBCyHc49o1nWgBoE8f+OMPmDXrpuAL5iTOjxi+ipiS9Vg7ZQ9vvQUPn5xApSmDYelSGD5cfzerx3bhgn6+Gja8kZyiZcsbz1t8vE5WMXcuhWrqAHxp51HdA7bbdRB2N5tNZ7CqWRPee0+nkgRTg29ion4p69YF/+2bOFW+Ec/uH4Rfg7p0rFGcEV1COFWzPhOaduNUzfpuSdQhPWAhhOstXaqDyGuvwauvQteu6R7qzsT5Z39fwvm+7xF6YT0x/hWY/e15qvZKdUDTpuYFXqXg9991jeNz53QykrQKLAQFQd68cPQoxdt1ASB+/zHoUgd69HBrk4GbE2o88ABMmPCfD1nutn69blLQjlVEFX+X6mdXY+QrpT8APvwwAJ2Lub80oQRgIYTrnD+vA8eUKXqvZ+3amd7EHYnzlYL1j46m0cy3KIrC5uNLybmTCQ4zb17yP156SQevhg1h4ULddUuLYeje7rFj5K1cBjsG6shRaPwsNG7s1iZz7ZrunYNHJNS4elV3wL/8EjoWjWS2fxiWM3G6atYPP+iV9yaSIWghhGv8/bde5fLTT7o83rZtDs2jOrtO7q3274dWrcA6MxwDhQFYDIXv5g1Ouf8csVr1nC7cKDgRGZl+8E1Rtqze1uPnx6gyY4jIk5wPUyl9n64WFaXPlTevztMdHa33c5sYfBctglq14IcvLhBRoy/Tn5iNxZZ444BNm0xrWwoJwEII18ibVweGf/6BkSMh0LEerKsS5ydZFfO7/cgLNVcRFQWx/T7UbbJYXLLIKssrudetg7vvhvff17+Hhjpe4zjViudV9V5jydXknnzp0nqFtKtcvAgvvqgTpvzyi76sVSsoVsx158zEuXO6aFZYGHRMmMmpIjVosXsSAQXz6NfZRa93dsgQtBDCOex2mDhRF1AfNkxvK/rnn2z1gpydOH/nvINceuJ5OlxeTN5yPam6rjllyoTCI0szrNebXVlayX3xIgwapJ+74GBo0SLrJ3znHZ1LGahR7DTXVh8D6uv54WPHcvBIMhAerofJT57UCTWS51LNEB4Vw6cL9rB/XRCxy2pRMu4MO6q+Qs09s/XK8cnJWdXat3fJ651dEoCFEDmXsrdjzRrd9bDZdE/DzIT6kZFYFy5j/Z+nqbd5MsrwIar3eFp+88KNsT8XLbLKaCX3TQE4IkIP1Z46pRenDR8OBQpk/YR33XX9x0f3fsR7F6dw7dol8pYt65q9wK+9pofHa9fWgbhhQ+efw0HhUTH0/2EvZWfG8WzMetYUOcszd3xLtf1/wSefQL9+N7ZhmbmoLg0SgIUQ2ZeYqIeXP/oI8ufXC1t69jS/kk1kJLaWrfFJSOBe7OwrcQ/Fl/xOvZBybjm9wyu5S5WCO+6AuXN1+s3sOn1aF7Fo3x6fO8pRMPIy+3dfonK5crB4cfbvNzWl9J5iPz+9urlUKb3Azg0pQ1MKJdya9tNuhzc/uEzJeb78nfQIfiRiveTLc/e/x9/tHmTaWz1d3rackAAshMi+w4fh44/1tqIvvoASJcxuEbGnElj+1E88mJCIL3Zs+JD4aHOCshh803vTd0R6K7nLFvTXz9PWrXpleLVqsHZtzj+wnDoFffvC9OkE3qX3Ap+JOkZlZyXj2LtXbyNr3lz30tu1018Oyslzmd5wfswhX6aNKcmRVXcyIaAvAUnxGAC2JGqf3M+EiuYX8siMLMISQmTNpUswebL++a679PDzr796RPBd8fEaTgfXpeWBX7EaviQZPiT6+vLh1ZJZSmeZ8qYfExuH4sabvqP3kdZK7rvPHmTuL/31fOnp0zqBBjhntKCsDrocO0ahWvqDxqVdR/WCqIED9UhFdiQm6tGN2rVhyxaoVCnLd5HT5/LW4XxlMzi5sgKvdCuGz5bNbAmqQ9v45djdnEbSGaQHLIRw3Ny5uqd1/LjeYxoSAhUrmt0qTu27RFS7QbQ78DWHLeXpe//7xJfiRq3eklU4cuv8awYcnsNNR8oxoxbu4cLpC7z/z+88vnYWRvHiepvOo486d5g+KAjy5YOjRynWQSfjSNh/DFo+e2NfblZt2aKnE7Zvh0ce0XO+pUtn+W5y+lymHrZPOFGIc3/XxnqmIO3KzeCv491JKFyUVx95j2MBhdyaRtIZJAALITJ38qReJDRjht5cOXOmDr4mUwp+++oM971Rn7Yqhn/ueY2nGzTnWmAAwPVavZC1TFrOyMZ1fSX3mTNQ4wW9SG3kSB0snc0wru8FzlOxDC8V/JmCgc3opBTExuoFcQULZv0+L1+GOXN0qstscuS5zGiIukxQIEdPJ1JlwRUaRm9lc0AsUV1KcbV+IYzEdwl4/XVaHbrGqIV7mBBcnTJBgYxwc2nI7JIALITImM2mt8YcOQIffggDBuh9lCY79ttKFry3hu8OhKLKPUnz0Z1p2K0JhUcu41oOM2nlOBvXyZO63M7Qobom7549UKSIw+fPlpS9wH5+bKjyFEVjgdgL+ryjR1/fppShv/7SVZ9GjoQ6dWDfvhwXcsjsucxsy1bboNos+Wg7c690IYA4jHjoHfAhD7bvBfV0JqvOhQt7RcC9lcwBCyHSduDAje1E48frTFbvvmt68E2yKtbcP5gyT9xHrwPvscqvNd1/e4jy3ZoAzsmkle37sNth0iSoXh1GjYLNm/Xlrg6+oFeg//03AM0LbqXU7ogbQ9OZ7QU+dUqnjezQAf78U/d8wSlVlDJ7LtMbov549r888wwMfbEYr6tvCSQOH0Bh0L/ABa8MuLeSACyEuFliol7ZXKMGfP21vqxNG6hSxdx2Abv+OsSGYg/QbOlwDMAXO772RHxWRlw/xhmZtLJ1H9HROvnI88/r1JHbtrl3f2xw8PXh7acPf8DbR1+6aWg6TUrp1djVq+vCBEOH6g8N2dmLnI7Mnstbh6iVgqu7S7Hp80bM+PEa26t2pevV3zAMAywWfAIDqP50+sU8vIkMQQshbkgpG7Njh15488gjZrcI0OmRFz80jlaLB1LeMNjTsR9Vlk5It1avMzJpZek+lNK5m48fh++/17kQ3b0Xetcu+PlnePNNkkqXo+K/i7l4EQqlSlP5H2fO6FXZISE3eu4ukNFzmXqIOulyHs4vrkXcvlLkK3OJVasKUOsDO/QaAc2awerVHpPFyhkkAAshtFGjdNGE4OAcL7xxphUr9Pql5/cdpELZFpSfP5Fom4XPylag8q6N7K/RgPYB5elsRuNWr9apDvPl09V/SpUybzvW4cN67rZTJyx3lKXg6svsir5EobJlb07GYbXqbWM9eui2RkZC1argY86A6ICwqgycuZ0zG4O5sLwad9oOMK5Yd+xfj6BuvSa6Z57yYaZ5c1Pa6CoSgIW43aXM8zZurHP7fvRR1lfMZoGjSRmuzF/B0T4fMO14F2yVXqLegpHUbutL+JbjetFOoUosaqr3pa5KL8+yq1y4oIscTJ4MH3ygCyg4UGrRpVLtBQ6sovcCn9tyFJ588kaPcf16/Wlm+3a9pahtW5f1eh1VI18wlr+LUX3zBt4PaEeobSU+8YH4+V/QB5idVc2FJAALcbs6cULn9A0OhjFj9EpnBwoBuCKrEdwcPKP6TqLOhBeojmKsz0qsk+sT0FIHkZzuK80RpWD6dP28nT2rV4Q7srrYHVIC8NGjBNXSdYAvRx+DF8L0oqrXXoOxY6FMGT3CYXItXKtVL84eNgx6+fzOBON/GPFK98Sn/qLTXeZymY45GIZRzjCM5YZhRBuGsdMwjNeSLx9qGEaMYRhbkr/au765Qogcs9vh2291z2fuXChZ0uGbOjurEdwIngCn9l/m78qvUGfC8xgoQNfqDVgXcf14Z+zRzbbBg+Hxx2+UWfz0U1120ROkWvFctGVtWhnLifJvDAkJuir9V1/pJCq7dpk+vbBxo16f9s470LEjfFXtawylX28MA3buNLV97uLIoH8S8KZSqjrQBHjJMIwaydeNUUrVTf76y2WtFEI4x/79OjNSnz567nLbNp2qMJWM6thmFkAzk16QjLkQx5QpMLLOb4T9O57dtR5Jt1ZventxXZZ60GbTJQMBnnoKPv9c1+6tV88158uulBXPp07hW7gA+4JD2XcmSOeBttl0zulx41w6vZCZa9f0oEHjxlDp6AqWjtrMH3+A/6cfubQ2s6fKdAhaKXUCOJH882XDMKIB79+AJcTtyG7XyRUmT4ZnnvnP/FpmQ8Q57X2mlZShwIkEyi7JwzPH4b57n+VIv4bUeLieXhyURu3WAWFVb2ojZH2fr8O2bNGrwitUgD/+0IuVqnpwisPNm3UgAx7P9yc+W/whX5gOvCZbulR/7jt74CLLq71Fi92TYFUn6D8H7r9fH+BBtXrdwVAp3X5HDjaMCsBKoBbQD+gFXAI2onvJF9K4TR+gD0D58uXvPnz4cI4bLYTIgtWrdc3W0aP17wkJkCdPmoc2G7kszaxFwUGBrBnYKtPrM5MS4Ksf2kHjw9tRR/Lz7OHp2Axf5n55iN4v5XFoMW5O5qEdcu2anpz87DMoWlQP33br5lULgvYXbcSJ+MI0v7rQ1HZcuKCrFu76PpK3839Ne58F+F85r+v0DhvmOUP4V6/qZCZ9+zr1dTYMY5NSKs1akw4HYMMw8gMrgI+UUrMMwygJnAUUMBworZR6JqP7aNCggdq4cWOWGi+EyKbYWD28/M03uubsP//otIgZqDhwPmm9IxjAwZEd/tNDBt37zEqiixU/zKVx70fxtyXiA/ybtzr55kyn1P21HH5oLrVliy6veOCA7v1++ikULmx2qxyzeDH89BNMmcL26t3w/Xc31Wy7TPncoJROGf7yy1D5TCQRhOJrT9TBLWUExpNs2aInpiMi9J5jJ8koADu08cswDD9gJjBVKTULQCl1SillU0rZgW+BRs5qsBAiB5TSw6XVq+vFVv366UUtmQRfyHx+NadZpuLjwfhmNXmSg68yfKj0zpOeE3xBrwovVQqWL9fPn7cEX9AfGn75BU6cIKlMOYLVMc6dc38zjh+HLl3g0UcV9YsdYdqLEfgayR/afHx06ktPcPr0jdKadevCv/86NfhmxpFV0AbwHRCtlPo81eWp61I9DOxwfvOEEFl26RK8+KLe57lhgx5GzZfPoZs6kgO5c71g1gxsxcGRHVgzsJXDwXfN/Fjq1IG31z2MzfBDWSwYAXkwWmU+dO1SSumg1aGDXqxUvDisWeOdC4HK6f2/HDuG7x1lKchlYqIvue30qVNh7/nrXw5WbsP8800IfqiBXlzlKYuslIIff9QNfeklXWgEoHx5tzbDkR5wM6AH0OqWLUefGoax3TCMbUBL4A1XNlQIkQGbTQcRmw0KFdLpozZsgLvvztLdOCOP8q0unU1kXsNhhHQsT9kru/lgUVN816zAGD5cL7wxc8HNgQPQrp3OCnX+PKZ0F50pVTKOvFV1MD6/NZ00lE62bx+0bg19n09iZLHR7PAJocKpDRiDB+srli4FT3jN9+/Xuc179dIBOCrK7YH3OqWU277uvvtuJYRwsqgopRo2VAqU+uMPs1tzk5WfRqpo35pKgdpY7Ql15eBps5ukWa1KjR6tVGCgUgUKKDVunFJJSWa3KufOn9d/B599pk7uiVUV+VeN/9Lq0lPOWH9MdW48XQ3y+VCF+c1Th0rX1W3o1Empo0ddeu4si4tTqkQJpQoWVGrCBKVsNpefEtio0omJkglLCG9160rd337TE28e4Py8tZz639s0O7uaU37l2PPZPO7u18HsZt2QlKTHStu00aUWU3qO3i4oSCdWSUigeOVCHPMrxJHjrjvdZ7+eZu4ru/j7/NP4k4jVx5cNeUM4NXICjd563nNWje/apXu7AQF6pXPt2nqu32RSjlAIb/Xoo3qFbq9euhTe44+b/oanFHz+3EICHryfqmfXYDcs7P38U6p6QvC9dk0PgV69qt+I167V27NyS/AF/fqfPAmDBuHjA4MLjKHI+r+dfppr13Tdjv49ivHUpd8JJA5fbPjZrKwPrsEbqorpf4sAXLmiFyGGhOhiGaBTXHpA8AUJwEJ4l9On9ZsK6AIAERF6Fac7Cr5n4vCms8wr8SQFJv+BP4n4oADFplkLHU5V6TKLF0OtWjqV5Pz5+rKiRT0jSLjQC1dGUTP6D6fe5/LlugM58dOLfBv0NM8lfY8CkgwDq8WXdeVD3JMWNDMLFujXfMwYnQHkwQfNbtF/SAAWwhsoBd99B9Wq6SAC0KSJLgBvMluS4q+nfiVvg+q0OzudhDsSsfpaSDJ8sFp8WR1c0+FUlU539iz07KkLD/j56Q8s3bqZ0xZ3mTTpeh3niwXKUuDiMafcbWysLqTUqhW0ujKXU8Vq8r8Lv/Jtw848+dhHfN68B08+/hGbg6u7Li2oowYM0D3dwEBYtQomTNCLEz2MzAEL4el274bnn4eVK3U91OeeM7tF10UvPMK5x16k/cW/2BPUiCfa9mZfxTIsiwmhyZHtrCsfwubg6hhm9Yief14XnHj3XV2QICDAnHa409GjMHs2JCURV6wcxc9FY7fnrNzv7Nl6t86pUzq2fXThb/w2FCVizDd8vsdCnNVGZIU6gAvTgmZGKb0LwNdXb3PKm1dXe0gn65snkAAshCebOlVnDMqbVyeFeOYZ0wqnp5YQEcnKDyLYHnGW51UEm3uMod73r3Bt9AqIjWNzcHU2B9+oM+vWHtGhQzrQlioFn3wCQ4fqOcDbRblyekPuiRPYS5el/J7FnDqlt4Vn1YkTOpPV8Vlr+bHQ59z1zgNUGP4sXB0N/v6E+vkxwtVpQR1x8KD+sNWsGQwZovd0d/CAdQeZkAAshCdKStKf5Bs1gsceg1GjslQ20JX2DP2Nih/8j5YqiRYWfxKm/Er9Hp0BNxdKuJXNpuvdvvsudO6sP7xUruz683qaVHuBfSuWo2DEZfZGX6J0acerICkF33+vczi3u/IHa4zH8LlohxHh0L7GTft4O9cLdn/ATZGUBF9+qadlLJbrQ+/ewvyP0kKIG86f10PMKfOUd92lc/t6QPC9dDaReY0+4M5hPfBTCfhiIw+JFDwWff0YVyTycMi2bToovPGGHn4cMcK15/NkKdmwjh7F1udFAojj0HnHg+/+/TpvRp/eNj4s8jlTeQIfZb9xQESEc9ubXTt26HUQ/fvrBu/apRdbeRHpAQvhCZTS+3jfeENnY+rXT/foLJbMb+sGaz5fT9G3e9MxaQd7yt1PlTOrwWpNM62g23tEM2fqLViFC8Ovv3rEdixTlS2r97z6+BBcJR8J6GnhzCQl6QXDgwfrl3V5zx9o8dObelh306Z0X2/TJCbqLVfTp+uerze+5ull6HDFl2TCEiINx44p1batzh7UqJFSW7aY3aLrTp1Sqv8DO5QNQ53wLauiR/2pr1i7VqmPP9bfzRIfr7+fOaPUSy8pdfaseW3xUPZLl9V431fV110WZ3hcVJRS9esrlYc49Wqr7SomRimVmKjU3LlK2e2e8XorpdTSpUoNHnzj94QE89riIDLIhCUBWAiznT2rVKVKSo0d6zHpEO1r1qrt7QeosIJrlZ+fUnM6TVYJZy6a3SztwgWl+vRRqnFjj3m+PFZiorJhqOnVB6d59bVrSg0cqFQzn7Xqx4A+6kqx8speurRSV6+6uaGZOH9eqWee0SGrcmWlLl0yu0UOyygAyxC0EGZYtw6+/lqvdClaFPbs0YuuPMCpHxdQrFdHamJjjjGW478so+ITz5rdLC31fph+/fS4qYcM03uU/v3h8GGYMYML/qXIc/a/e4FXrNDLDcrvW8IK4wEs8UmQYOjUpnnzmtDoNKSU1nzlFb2n++239SrnQJP3GTuJLMISwp0uXoS+feGee3RKoUOH9OUeEHxtSYq/n/6dfL0ewQcbBuDvY6Xi4Qizm6YXp3XtqnNdlyihKz2NGuXRezxNde4cREYCcLFgWQpdujEJHBurd+yEhkLRhOP8na8rFpWkr/Tx0UWbPcW5c/Dsszp15MaNMHJkrgm+IAFYCPdQCmbM0JmsvvkGXn1Vr9r0kG0yO7fbWVWiCw/81J1z+cqj8uQBiwXDUxbd5M2rSweOGAH//JPlMou3nbJl9SbepCTO5C9F8YRjVBjwF9V6bKPSXTZ++NZK//6wdFdp/Dq00x9kPKVWr92uRzqUgmLFdAKa9euhbl1z2+UCEoCFcIekJL28tHRp/WbyxRdQoIDZrSIhXjFkCNS724d/4muz6anPKR+7HZ/ly82v3bpvHzz1lM59HRCge0ADB+qUkiJjyck4Fi6OYhuF8CeRU380ZM8vITycMJkzRSoyqu9B8uYzYNo0PRpj9usNeiomNFSPdKTk7K5b1yNGiFwivclhV3zJIixxW0lM1AurUhaMHD6s69B6iM2/RquNgfeqlixVTz2lFxN7BKtVqZEjlQoIUKpQIaXWrDG7Rd5n/nylQD33wleqxKORqglr1SfGm2pDgRClQG25o5ZS+/aZ3cobEhOV+vBDpfz9lSpcWKkpU/Tq61wAWYQlhJutX68n2rZu1b233r2hfHmzWwXA1QUrOfLccGoei+CaTwFGv3eR+sPNblWyzZv1cxUVBQ8/DOPGQZkyZrfK+1SuDG3acOpaEvf4rWGaMQg/ZYXLMKlhZ0a2fIYDHjL9AejXev58nYDmq688IvGMO8gQtBDOdPGiTp7btKletTlrll5E4iG2vf4dgQ+EUv3YEnwNhf/vP1N/+MNmN+uGt97Sc5czZ+rnToJv9lSpAosWcbZ6HUIPbcRPWTEAm+FDbGBBShfOZ3YLdV3mhAT982uv6drM06bdNsEXJAAL4VwdOsD48XrbxK5d+pO9B2ToOXMGnnwSDn05BwMF6AWvefdvM7ll6PnH48f1z1Om6OetSxdz25RLDAiryroqjUjw9b9eHjKqUl1zqhWltnixLpCRkjK0TRt46CFz22QCGYIWwpnmz4eYGKhRw+yWAHoh6bKBi/jp6yvMSOhCq/+9Db8tAWui+SteY2P1ftXvvtNbs8aPv5HHWORcWBidixWD/p/yaoAvlXdtZH+NBjzWp4t5xRPOn9f7t3/8UffSW7Uypx0ewtBzxO7RoEEDtXHjRredT4jb2dEt59jd4U3aHP+R7fmb4hO5hpq1DL0/NCJCB1+zVrzOmqWH6k+fhjff1CUDc9H+To/Qpg1cvqyTvniCRYugRw+9t/ftt+H992+L+syGYWxSSjVI6zrpAQvhLPPn6ze7wYNN3SpjW7WW/a+PpfjmhYRymX/avkv9We9hyZc8FN60qblbTcaO1fug69aFefOgfn3z2pKblS2rh3o9RYkSUKmSDsR16pjdGo8gAVgIZ5k5UweUDz4wrQkHpkZStkcrqqoE7Bic++wHGvbref36cLOKpyulhx+LFoXu3XVlnVdekT29rlSu3PVkHBnto3XZ34TdDt9+q+f0v/xSf+Bau9Yj1kR4ClmEJYSzREXp3pwJbzAJcXa+6bOJ73tG4JOcVtCw+FA8Ieb6MeFRMQyatZ2Y2DgUEBMbx6BZ2wmPiknnXp3kwAE9HNqunQ4GxYrpeUAJvq5VtqwOgidOpHuIy/4m9u7V87svvKDr9qasdpbgexMJwEI4Q0KCfqOpV8/tp476fQ/bi4Xyv2+bUuyeu7AE+KeZRnLUwj3EWW033TbOamPUwj2uaZjNBp9/DrVq6fSRvXvrpdfCPerVg2eeyTDoOf1vwmrV+Zpr14YtW2DyZFiyRHJ2p0OGoIVwhp07ISmJ9w/78cvA+W4Z3r1ywUpEh1HcH/kBCUYgu1//htc/7wrrgtNcZHU8Ni7N+0nv8hw5dkxvJfrnH3jwQV35qWxZ559HpK9hQ/2VAaf/TZw5Ax9/rLfjjRunU6+KdMnHUSGcYN3KbVz1D2RlgXIuHd4Nj4qhz4tjGVHrBWKK1aJj5LvsqtwJy95oao/5H+FbjtNsRRwVL9am2Yq4m85fJijtVcbpXZ4jxYrpVc2//w5z5kjwNYtSGVY3csrfRHw8TJqkz1WmDGzfrtdDSPDNlARgIZzgzfjy1Hp9GoeDbrzpOHt4Nzwqhmkjp/PFN28xYOdkKtoP8HH9vhyZPob8lUtlOp83IKwqgX43184N9LM4LynD2rXwwAM3iidERMBjj8m8n5mKF4dBg9K9Osd/E6tW6RXNzz+vfwa4447stva2IwFYCCc4HhuHMnz+E2ycNbyrFIS/8Tfjpg8nj0rEFxs+hh3f/FeuB/nM5vM61wtmRJcQgoMCMYDgoEBGdAnJ+TD5lSs6leC99+qh+IMH9eUSeM1XvDgcPZru1Y78TYRHxdBs5DIqDpxPs5HL9Ae6S5fgpZegRQs977t4sf5ZZInMAQuRUzYb4dMGMjnkAf6scd9NVzljeDdm+3l2tu/PD8emcNQnmPxcQSmwWnxZVz7kepB3ZD6vc71g585LL1oEffrAkSP6Dfnjjz2izKJIVrasno/PQEZ/EymjKikf7FJGVZrPG0zRbZvgjTd0GcN8HpBb2gtJABYip/bto86hHeQPaXvTxTkd3rXbYcEL4dw9+QVaqbOMv/N5vnywLbXOHqDJke2sKx/C5uDqBCcH+TJBgcSkEYRdMscLuls+cqQebl61Cpo1c815RPaVKwcLF2b75qlHVQpfu8hV/7zEAUMaPs64iWOhcWMnNfT2JEPQQuRUVBQAbZ5s57Th3ehoaN4c5n8bw6X8wZyev5HgGe9jyZ+XzcHV+bppNzYHV78pyLt8jjfFzJk637VhwK+/6u0mEnw9U9myeh+w1Zqtmx+PjQOl6LRrBUsmv8hLkdMBmF+kqgRfJ5AesBA5tXkz5MlDqy6htMphcglrxBr29x3N9D112R00hOenvEjlJ5/H8POlc/Ix6WUtSvnuskxXJ07o/M2zZun8zaNHQ6lSzrlv4RqtW+u911ZrthKf1DGu8PIfn3H/v/+wpXQV5lfTH7RcNqpym5FiDELk1P3368o+Ofzb3vfRNCq99wQW7NgMCxf/XEWRDibmbE6hFPzwg85eFR8Pw4bpnzNIbyhygdmzsfboSVKCldEtejDl7gex+1gI9LM4Z/HebSKjYgwyBC1ETt15J7Rvn+2bX421Mu/ekVR470l8sANg8YEi2yKc1MAcGj1aZ1QKCYGtW+GttyT4eguldMWp2Nis37ZSJfzubcaqWctY0KY7ysfivJXzApAesBCmWrQIfuy5lKmn7mdPmVCqnF+HYbXqWr1Ll5pXtchm02XjSpTQ32fNgmeflVSS3iY2FgoXhs8+06MWGUlKgi++gH//hQkT3NG624L0gIVwFZst82PScO7oNT5ts5iwMNgU1Jot36ynasxyjGXL9LYOM4Nvygqw9u31m3LRovDccxJ8vVGhQnqLUAZ7gQGdvappUxgwAI4fz/aiLZE18h8lRE588AGUL+/wG5ZaG8m/LXuTcEcVXlvSkRGvnmDLFqjbp5FOeJBOGkm3sFrho4902bg9e+D118FiyexWwpMZht6KlN5e4IQEGDJEV/E6fBimTYPwcKlU5SYykSNETkRFQf78Dr1hnfllIUV6duBOZcOOwfFBYxn4sU5dmV7CA8A9821HjkCnTnqOt1s3GDtWDz8L75dRMo6zZ/Ww8+OPw5gxOoe3cBvpAQuRE1FRmZYgtNvh2y+u4t+jGz5KB1jD4kPZApeuH+P2UoG3KllSv/nOnq17QRJ8c49y5W4egr52DcaP1wu0goP1lMPPP0vwNYEEYCGy6+xZ3bPIIADv3XiJ++6DPm/kY8EdfVD+edKs1evWUoEpVq+Gtm3h8mVdr3XJEujc2XXnE+bo0UOvKwBYvlyvZn/5Zf36g65gJEwhAViI7ErOgJVWALYm2Jnf+VuKN7yDoluWMmUKdDs4Cp+I5WkusnJrqcArV+DVV3Xy/H379NyfyL1attS1mZ9/Hlq10vPCy5frhXbCVDIHLER2lSqlA1nqABwZyYlxf3Bh1nI6xEexs3gok+feQbEmydc3bZrm6uYBYVVvmgMGF6WRXLxYF084fBheeUUvusqf37nnEJ4lPl6XJJw0Cfr314lU8uY1u1UCCcBCZF9ICHz55Y3fIyOxNb+PUjYrpYC9XQdSc8bHDpXlc3kaSbhRPCFPHimecDux2fQiwQ0boGFDs1sjUpEALER2RUdD5crXV0Dbl0dg2JIwAGWxUOXugjcF3/ComAwDrNNLBaaYOxfuvlsvuPn1V703NCDA+ecRnilfvps/KAqPIXPAQmTHlStQs6buUSaLqRxKPAHYff67yCplm1FMbByKG9uMXLrX98wZ6N4dHnpIp5MEvdpZgq8QHkECsBDZsXWrHtJNNf+76Wo12jOfEy/+d5GVW7cZKQW//QY1augUksOHw6efOv88QogckSFoIbJj82b9PVUALjx5FAv4HD6+BAX9bzrcrduMvvpKZ7Fq3Bi+/14HYiGEx5EALER2REXpZBWp9lDm2xfF4YCqVL0l+ILeThSTRrB12jYjpeD8eZ23uUcPnbe5b19JJSmEB5MhaCGyIyUDVqpFVuXPRXGqTP00Dx8QVpVAv5uDodO2GR0+DGFhOqlGUhIUKaK3GEnwFcKjSQ9YiOwYNeqmmrjndpyghP0UO0LSzorlkm1GdrsuGzdwoP7900+lYpEQXkQCsBDZcf/9N/16dM5migIF70s/LaVTtxmdPKmLJqxapXu+kybBHXc4576FEG4hH5eFyKqtW2HhwptqAf9zrSavM4YKneu6pw2FC+vzf/89LFggwVcILyQBWIismjQJHn30pvnfVUcr8Efw6xSrWMB15925Ex555EbxhNWr4X//cyjTlhDC80gAFiKroqJ00fpU8635Vi0gtNpJ15zPatU5m+vXh4gInYELJPAK4eUkAAuRFTabHoJOtf834eQFJhx6gB62H5x/vi1boFEjeO89XSpw1y79uxDC62UagA3DKGcYxnLDMKINw9hpGMZryZcXMQxjsWEY+5K/F3Z9c4Uw2b59uqB5qgB8ZO4WAALvSX8BVrYNHAgnTsDMmTBtmt57LITIFRzpAScBbyqlqgNNgJcMw6gBDASWKqXuApYm/y5E7pZGDeALy/RlwR2dFIA3bIBjx/TPkyfrXm+XLs65byGEx8g0ACulTiilNif/fBmIBoKBh4Afkw/7EejsojYK4Tm6dtXDwqnSOxpbo4gxgqnQKIe907g4ePttnUP6/ff1ZWXL6sQaQohcJ0v7gA3DqADUA9YDJZVSJ0AHacMw0nz3MQyjD9AHoHz58jlqrBCm8/eHOnVuuqjYkc0cDKpHcE4ST61dC888A3v2QO/eN6oXCSFyLYcXYRmGkR+YCbyulLrk6O2UUpOUUg2UUg2KFy+enTYK4RmUgrfegjVrbrqoq084y9uMyP79/vYb3HsvxMfD4sXw7be6Zq8QIldzKAAbhuGHDr5TlVKzki8+ZRhG6eTrSwOnXdNEITyA3Q7jxukUlFu3Xr/4yBGIunIXJVrVyvp9xsfr72FhMGAAbN/+nwxbQojcy5FV0AbwHRCtlPo81VVzgaeTf34amOP85gnhAf79F1q3hldfhTZt4Mknr1919JcVvMAE6tWyOn5/ly/DSy/pXq/Vqud4P/kECrgwiYcQwuM40gNuBvQAWhmGsSX5qz0wEmhjGMY+oE3y70LkLnv2QO3auv7vt9/qFJSphocDZ03lI96lZh0Hl1MsWQIhIbqIQvPmN6WzFELcXjJ911BKrQbSS7nT2rnNEcJDXLsGefNClSp6L26vXlCu3H8OK3Qgij1569E0fyZZqa5cgX79dBCvUkWnkbznHte0XQjhFSQTlhCp2Wzw2WdQoQIcPKjTPb7/fprBF6uVche3c7asA/t//fxg/Xq9iGvLFgm+QggJwEJct3u3npft31/vxQ0IyPDwy//sJo9KwF63ftoHxMbq+7p0SRdP2LBBz/UGBjq/7UIIryMBWAildDH7unVh716YOhXCw6F06QxvdmzpHgAKt0qjBzxvHtSsCV98oQsogA7CQgiRTAKwEIYBBw5Ax4467eMTTzhUaWhJ0CMUIpbKHareuPD8eejRAx58EIoW1cPOnTq5sPFCCG+VpUxYQuQaViuMHKn34DZqBGPH6nnaLNiyBfIUL0Tp4FQXduyog+7gwfDuuzpzlhBCpEECsLj9REXpQvZbt0JCgg7AWQy+2O08PvtxytzRE8PoeOPyxYt19aLKlZ3bZiFEriND0OL2kZCgVzQ3agSnTul53g8/zNZdWfcepM2FGdQucfLmK/Llk+ArhHCIBGBx+/j+ex1wn3gCdu6Ehx7K9l2dmL8ZgPzNUy3AmjtX7xm2ZiErlhDitiUBWORucXGwbZv++bnnYNky+PHHHJf4u7QiCiu+3NEhVQ7o2bNhyhTwlZkdIUTmJACL3Gv1al06sF07HYh9faFlS6fcte/2KKKNGlQJSbW1KCoK6td3aAW1EEJIABa5z9WrunBCixZ6OPinn5ye/OLs1UB2F2t+o7ObkKCHtes5kBVLCCGQVdAitzl1SmexOngQXnkFPv4Y8ud36imUgi7MolMn6JZy4Y4dkJQkAVgI4TAJwCJ3sNvBxwdKlIAOHaBbN11tyAVOnIAzZ/To9k0XBgXpIWghhHCADEEL7/f33zrt44EDev517FiXBV+AS4NHs55GN9cA7thRZ8GqVMll5xVC5C4SgIX3On8enn4a2rcHi0WX/HOHdZEU5gIh9W9J3mEYsgBLCOEwCcDCO82eDTVqwK+/6uQamzZB7drpHh4eFUOzkcuoOHA+zUYuIzwqJtunLnwoij356lOoUPIFNhs0bqzbIoQQDpI5YOGdFi+GMmVgwQJdxSgD4VExDJq1nTirDYCY2DgGzdoOQOd6wRnd9L8uXKDk1YNcqNHnxmV79+pSg5KAQwiRBdIDFt5BKV0mcMMG/fvo0broQSbBF2DUwj3Xg2+KOKuNUQv3ZLkZcZFbADDuTrXYKipKf5cFWEKILJAesPB8x47BCy/A/Pl6zrdRI8ib1+GbH4+Ny/Ty8KgYRi3cw/HYOMoEBTIgrGqaveP9J/IRzaMUvT/VdqPNm3Wt32rVHH9MQojbnvSAhedSCr79Vq9wXrYMPvsMvvsuy3dTJijtJBwpl6cMUcfExqG4MUSd1jzx2qRGPMZ0atxX/MaFUVEQEpL1ikpCiNuaBGDhuaZOhT599NDutm3Qr59e7ZxFA8KqEuh38+0C/SwMCKsKZG2IOuHPhQwLGEH5mMgbF1arBg8+mOV2CSFubzIELTyLzaazWFWuDI89pnuVjz6qk2xkU8pQcnpDzI4MUQOwcCGvzG+HwsC4PwCWLtVZt8aPz3bbhBC3LwnAwnNER8Ozz+qEGnv2QKFCOgg7Qed6wemueC4TFEhMGkH4pqHr+fNRTzwBgA8KEhMhIgIaNNBFHmT/rxAii2QIWpjPatU5m+vW1YF31CgoWNBtp89wiPrsWV0/uGNHzqqiJJAHu48F/P0hNBSGD4fgYNmCJITIMgnAwlwXLuhVze++C506wa5d0KOHW3uUnesFM6JLCMFBgRhAcFAgI7qE6B6zzUb8guWMKjCMcld2M/mJ5ahhw28MP0dFQeHCsgBLCJFlMgQtzKGUDrJBQXD33TqbVZcuaR7q6BahnLhpiPrwYfj6K06WHMErr5Vk/oV/uat2XlZ/Bw0aNAWa3rhhVJTuCQshRBZJD1i439q10LDhjeIJkydnGHwd3SKUY6tXQ8eOqGrVsH45nq7Vd/Lnn/D+x3nZuFFP997kzBmIiZEShEKIbJEALNznyhV47TW4914dvE6fzvQmzsxilaFffoH77oP587HHJ/Jows9Y6oSwdSsMGpTOCHNKBiwJwEKIbJAhaOEeS5bAc8/BoUPw8st60VWBApnezOEtQjlht6P69QO7HQOwYzDood00nJXJ7qfgYHjzTQnAQohskQAs3GP2bL1yeNUq3QN2kENbhMjmPPH69VC7Nlv3BjIr32jePvM8/lixBPjT+O3QzMeHatbUOamFECIbJAAL1wkPh9Kldam+Tz/V3cnAtNNCpmdAWNWbKhnBzVmsIBvVji5fhnfeQY0fz5LmH9B+7XsUKdKTFsPvopVPBEbLUL3C+daHc0uQH17RRquHmkNAQJYekxBCgMwBC1c4dQq6dYOHH4YxY/Rl+fJlOfhCJluEkjk8TxwZCb16QeXKqPHj+bnQy3RZ+RpPPaV3P7V+rynGO4PSDb6pF4NdPHWOVo+3Jbrf+1l+TEIIAdIDFs6kFPz8M7z+Oly7pud5+/fP8d1mlMUKHJwnjozUi6ysVuwYPM9ElgT1YdZ0aNMm8zbcGuSrnz4AwHfXCiOD0EKI7JAALJxn2jRdLvCee3TVIjeV58twnlgpiI+HiAjsSTZ8ABs+dL73HF8s0B1zR9wa5Gue0gF4Vb6yOW2+EOI2JUPQImfsdti/X//8yCO6B7xqlVtr46aXSnJwSD7o0IH47v/jvSWhxKs8JGHBJ48/HT4NdTj4wn8XfdU8dYAz+YLwLevchCBCiNuHBGCRfbt3Q4sW0KwZXLyoixI89VSOKhdlx63zxOUK+jMt8R/aPtYa67KVDFvUlFGrmjC9z1IYPhzL8qVpzvNm5NYgX/P0v0SXqsyAdu77oCGEyF1kCFpkndWqCyYMG6bHcMeMcWvxhLR0rhdM5/gjMOtv+G0B7NjBP0XDePTcRMo2q8DWyVCt2i1pJLN4/3CjpOHXD73MY40rOD0lphDi9iEBWGTNhQvQsiVs3arr9H71FZQqZXar9CKr1q1RiYlgszPCdzAj4ofyyXiDF15wTqf85sVgHXJ+h0KI25oEYOGY1MUTGjWCIUP0NiNPsG4dvPIKKiERw27DioXgOwPYtdigXDkXnG/zZp0Dun17sFgyP14IIdIgc8AicytW6HSL//6rg/CkSZ4RfJNzS6t77uFK9BHi7X5YsWD4+9Pz+1DXBF/QxSOeesqtJROFELmPBGCRvosX4fnndbm9y5fh/HmzW3TDggVQsyZq7Fh+LdSX0tf+ZdQDy0h8dzi+EUsx7sneXK9DoqKgbl23LzYTQuQuMgQt0jZ3Lrz4Ipw8qQsOfPAB5M1rdqu0xETsfV/m9KW8dFWriCnUjD9+h7Cw7C+ycpjNBtu26cISQgiRAxKARdoWLYJixXQ+54YNXXoqhwoprF0L48dD797Mu9qST679zT+x5Xnx9TwMHw7587u0iTfs3auzfEkFJCFEDkkAFlpKGskqVaBJE108wc8vnUK4zuNQIYXZs3WSD7udpN9m8JFawZVaTVkxR9d5cCupASyEcBKZxBK6Rm+7djqN5KRJ+rK8eV0efCGTQgp2O4wfj3r8cZTdDoBSdj66P4JNm0wIvgCPPQY7d0KNGiacXAiRm0gAvp3ZbPDFF7qu7dq1MG6cXuHrRhkWUnj9dXj5ZXYH1CWeAJKwYAnwp9UHofj7u7WZN1gsOvj6yuCRECJnJADfzn7+Gd54Q69y3rkTXnrJ7St7b82x7GezUjD+CqULBvJT3hd4zv9HGtnWMb/fMnw+Go7PsqynkXQapfTztWqVOecXQuQq8jH+dpOQoBcShYTovaxFi0LHjqbtaR0QVpVBs7ZT/dAOHt6xnOYHN7OjaA365f2bp3f40759DXZOgPLl3bDCOSNWK3z2mR4xqFoVmjc3ry1CiFxBAvDtZO1a6N0bzp2DAwd0HucHHzS1SZ3rBVNsw2qafDwQi92GAr6+1I8LRf359Vd4/HEPyHexbp3eD71tG3TqBE8+aXKDhBC5gQxB3w4uX4ZXXoF774WrV+GHHxwvhOtqmzdz78AX8LXbMAAbFhrUSSI6Grp394DgGxWl6xufO6dXY8+ZAwUKmNwoIURuID3g3O7kSZ27+dgxHYQ/+siNm2Yzd6VQMLE+5ShGHL4k4ZPHn25fh0IxExullK5xfNddOuPV11/DE0+YXvFJCJG7SA84t7Ja9feSJXXVojVr4MsvzQ++SsHvv8PDDzP/TzvVQ0tS/sI2Jjy6HNuQ7NXqdapDh/SwfJ06cPiw7oK/8IIEXyGE00kPOLdJSajx7ruwfDlUrqwXD5ktMlJn1VqzBtas4d8iDflf+FlK1CzBjBnQpInJi6ySkvQCqyFDdND96CMIllq/QgjXkQCcmxw8qHtrixbpeUulzG6RtmYNtGyJSu6Vf5vnFV69NIZBQy0MGoR5e3pTJCTo52vzZt37HTcOypc3uVFCiNxOAnBu8eWX8M47eh/v+PE4rQq9MyxdirJaMYAkLNhKlGbT3xZq1jS5XYmJOvrnyaMD77vv6jKLpq/8EkLcDjzkHVrk2P790KoV7NoFffuaH3wTE+HLL7Fdvsb0C22IJwArFpSfP31+DTU/+IaH6+H5yEj9+9Ch0KWLBF8hhNtk2gM2DON7oCNwWilVK/myocBzwJnkw95RSv3lqkaKNMTFwfDh0KEDNGsGn3+u0yN6QgD55x949lnYvp0PxxZh6L89aFNuNq1K/8GB+nVon688nc1q27Fj8PLLejtR7doeMP4thLhdOTIE/QMwDvjplsvHKKVGO71FInMrVuh6tPv26eHTZs3cUjghU1evwvvvo778ksv5StHLMocl5zpQ+qGt7KlqY6/xMACrbq125C7ffAP9++sc2J98otNKesLzJoS4LWU6TqmUWgmcd0NbRGYuXNCBNzRUB5HFi/WqXRcKj4qh2chlVBw4n2YjlxEeFZP2gZGReiHTmDFMK9SHcpd3kffxTlR7aS3+1Y7d1DG/Xu3I3S5d0slIdu6Et96S4CuEMFVOJgpfNgxjm2EY3xuGUdhpLRLp++UXmDIFBgyA7dvh/vtderqUWr0xsXEobtTqvSkInzsHf/2Fat0a+7YdxJOH3/168tv8QvzyC5xJupTmfadXBcmprl7Vz9WMGfr3N9+Ev/6CihVdf24hhMhEdgPwBOBOoC5wAkh3o6lhGH0Mw9hoGMbGM2fOpHeYSE9MDKxcqX9+8UWdGvHTT3W9XhfLsFavUjBtGtSowaWXB2GLS8QHO35GEtNejKB9e338rdWOUqR3udP8/bcuszh6tH7OQC9M84Q5ciGEIJsBWCl1SillU0rZgW+BRhkcO0kp1UAp1aB48eLZbeftx26HCRN07dmePXWiCF9fXcXITdLrpdqPHIWHHoLHH+dAUjleOvgmVsMfu4+u15snLPT6sQPCqhLoZ7np9oF+FgaEVXVNo0+e1BUc2rfXH1JWroSPP3bNuYQQIgeyFYANwyid6teHgR3OaY4A9FaiFi30dqKGDWHpUlMKwKfVS212aAtLvu9L0sIlDM47mpqX1nHnkJ5Yli/F58Phuq2pUkl2rhfMiC4hBAcFYgDBQYGM6BLiugVYq1bpognDhumer5QNFEJ4KENlki3JMIzfgFB0evxTwJDk3+sCCjgEPK+UOpHZyRo0aKA2btyYk/bmftHROg9xgQJ6a1HPnqYNm6bMAVc/tIMmR7axrnxtjuW5k4/Df+L5c6Mp3vhOJk+GWrVMad4N0dG6VOBjj+mh8WPHoFw5kxslhBBgGMYmpVSDNK/LLAA7kwTgDJw6pQsnKAVjxsBTT0GJEma3ipXfzaJpn2742m3E+wTwgN8S/rE04+OP9XZaiyXz+3CZ+HgYMUJ/lSgB//6rt2UJIYSHyCgASyYss128qIeaK1XS2awMA/r184jgy8aNtHjvJfySa/X62q30LLeSnTvhtddMDr4rVuhSgR98AN266TzOEnyFEF5EArCZwsP1IqtvvoE+faBUKbNbpMXFQf/+qMaNuXopiQT8sWLB8Pfnfz+GUqGCye07eFCn3UxMhAUL9PYsT/jAIoQQWSDFGMxgt+v5yj/+0OkQw8P1YqtbhEfFMGrhHo7HxlEmKJABYVWdvngpzXNUK0LcjD/5s9BzPHfhE15vs4sBDSPI3zHUvFq9SsGmTdCggd7HO3MmtG3rlu1YQgjhChKA3UkpPcTs46OHnEeO1MPNaWRkSlkAlbIPNyUJBjgvhWPqRVbdDmyifOxJhpzux3cnWxJxZBNBZfMz9Sfo2NHkWr2HDulh+gULdJ7pu++Gzp3Na48QQjiBBGB32bVLlwj8+GOdDvGTTzI8PKMkGM4KwKMW7qH6wR38/tsg/OxJAKzZH8ZniQ/St69e21SwoFNOlT1JSfDVV/D++/qDy5gxet5XCCFyAQnArpaQoCPZxx/rrUXnzjl0s/SSYDgzhaPt6FE+WfAV/snBNwkLeSzxlHpyLePH3+O082SLUnqed9UqXfHp66+hfHlz2ySEEE4ki7BcafVq3WMbNkyv1I2O1hmkHOCOFI6fL51A+fMnScQPKxYSffzY1bkglWrFO+0cWRYXd2Oo/qmnYPp0+PNPCb5CiFxHesCutGGD3qv699/Qrl2WbjogrOpNc8DgpBSOe/dCoUIcs5ZkVIHJbFcBlCl6gHYV/mBj9SrsrlCDEa5KE5mZBQv0MP0nn+hFan36mNMOIYRwA+kBO5NSuvJOeLj+/dVXYceOLAdfcEEKR6sVRo5E1a5NdOeB1KgBi3bVokW/YiT18+Xb+x/iVM36rk0TmZ7Tp+HJJ+GBByAgAMqWde/5hRDCBJIJy1mOHIGXXoJ583Qg+esvs1t0w5Qp8O67cOIEEUW70v3cWELalOabbzygMt+MGbrXe/kyvPMODBokCTWEELlGRpmwZAg6p2w2GDdOBzil4LPPdM/XU3zwAWrIEAAS8WdE4puM/KG0mSmm/6t6dfj2W/1dCCFuExKAc2rJEnj9dT3MPGEC5qeJShYfDwEBHD+YQCkMfFD4Gjb+eCWCAk87d09vlhKGWK36Q0pAgH7eHnkEunbVe6OFEOI2Iu962XHtGixfrn9u21b//NdfnhF8L1yA3r2x3deSN1+38eiPHUkg4Hqt3gIdQ516upRkHjGxcShuJAwJj4r578EbNuhMVoMGwcaNNycmEUKI24y882XVwoW6/l6HDnD2rA4goaGeMZ47cybUqIF9yg9M3tOCcV8mUfv5ptgWpV2r1xkyShhy3eXLurfbpIneBz17ts7f7AnPmRBCmESGoB11+rROGzl1KlStqrcWFStmdqsgMlL3vleuhJUrOVykLg/b53O1ZH0Wz4UWLQCaQhvXpJJ0KGHIrl16nvzFFz0gvZYQQngGCcCOuHhR93pjY2HIEM9ZqRsZCa1boxITwW5nSsCLvHTxS/q948f77+tpVlcrExRITBpBuKYlDn74AXr1gsaNdalFTxiiF0IIDyEBOCNnzkDx4lCoEAwdCi1bes5K3f37oW9fVEIiht2GFQvXipRj3V9+1Knjvmb8J2GIUjyxaxlDV34PiQl6jrxMGQm+QghxCwnAaUlM1NmYPv4YFi/WxRP69nX6abJVbjApCT77DDV0KFa7BbvdFwuAnz8v/B6KrxuDL9yozDRq4R78Dv7L6KUTaPBvlH7OJk3SwVcIIcR/SAC+1erVOgVidLTO33znnS45TbbKDUZFwbPPQlQUq4o+zOPnxtG1wWGG3hdB0a6hptXq7VwvmM5VC0P5bnqb0YQJ+jmU1c1CCJEuCcCp9eunS97dcQfMnw/t27vsVA6XG4yMhIgICA3F/s67XNt3gt6+f7DQ1pXPv4devcpgGCbW6t2zB6pUgbx54fvvda3eYDenshRCCC8kATglFadh6MD75pu6elG+fC49rUOrhyMj9bxzUhI2X3/eLDWVH6+Ecv8jhYkeC6VKubSJGbt6VS9IGzMGfv8dHn0UOnUysUFCCOFdbu8xwkOH9H7e337Tv7/2Gowe7fLgCw6UG4yN1SktExLAZsOekEiJ87uZMrswM2aYHHwXLYKQEJ3R6rnnoE0bExsjhBDe6fYMwFarDrQ1a+r9s3HOK3LvqAFhVQn0s9x02fVyg7NnQ40aqE2bSMIXKxaUxZ9X/gilc2e3N/Vm/ftDWBj4+cGKFTBxIgQFmdwoIYTwPrffEPTGjdC7N2zdqodMx42DcuXc3ozUq4dvWgW9eha8+ipHitThYfUnd5ZNZGS7CCo9E4q/SYusUEp/+fjAPffoDcbvveeejcZCCJFL3X4B+MgRvb931izo3NnUdIid6wXrQKwUnD+PKlKUuZsfJyp/AiNiX+ONgX4MHgyBgSYusjp8WGewuu8+ePtt6NJFfwkhhMiR22MIevZsvTUG4OGHYe9e/d0TchHPmAGVK2O9pwVdO9t4qHdx5lbpT+RGP0aMgMC0p4pdz2aDr766MUwv6SOFEMKpcncP+OhReOUVmDNHp0N8/nk9jOqERVbZSqKRWlISvP46avx4ABR+nD+4jk8/bcYbb4Cvma9MdDQ88wysW+d5ZRaFECKXyJ0B2GaD8ePh3Xf1z59+qqvxOCkxhKNJNNIN0keO6B745s0AGIAPdma+upKiA5o5pY05cvEiHDigC0907+4ZIwVCCJHL5M4h6B07dMC9917YuRMGDNCrdp3EkRJ8GdXJtQYV5+j5fAyxfEgcgbpWb6C/zmblZOFRMTQbuYyKA+fTbOSytOv0AqxdCyNH6p+bNNFbtJ54QoKvEEK4SO4JwJcvw/Tp+uc6dWDTJl2mr2JFp5/KkSQatwbpxke2880v7zDym39pdF8g5Q+tYMdD73Jtrq7Va7igVm9GHwKuu3QJXn5Zf1j55hv9PIKJk89CCHF7yB1D0HPm6CBy/Dg0aACVKkG9ei47XXol+FIn10gJxvcejKL/yp+oe3IfB/OU4+I3RYktBTNnGsmLiZvCg65Z5Zxpusv58+GFFyAmRif9+PBDyJ/fJW0RQghxM+8OwMeO6UVW4eG6Xu/06Tr4uth/SvCRKolGsjJBgXRc8DMDV/wAgBVf/pfwI+cbFmbPIufkrshsIViGPfXTp3WxiYoV9UrsJk1y3iAhhBAO894h6Ph4aNgQFi7Uc5ebN7utGlDnesGM6BJCcFAgBhAcFMiILiE3Bb8BbavwVNTfgF5kBYoH6s7gm2+U04JvZsPL/0l3qRT3HoyiTKEAKFECli7Vz5sEXyGEcDvv7QEHBMDXX+v53mz0enO6jeh6Eo3UlIIff4SwMGwHy/Km+pYfeAw/ErH5+tLktXbcl5WtShlwpJpS6p568MXTfLxwHPcd3EzkvT/qG0jgFUII03hvAAa9lScbslWLNzMHDugauEuXMq3aYB7fPYy6dR/k6GtLqXoiAr/QUO5zYg/dkYVgnesFg83GoWGf8NzC7zEMg61vDafpS085rR1CCCGyx7sDcDY5XIs3M5GRsGwZnDyJ+u47rMqPtwImMungc4wcqcsL+/k1BZw/NO7IQjCAzh+/DvNmwwMPwMSJ1Clf3ultEUIIkXW3ZQB2qBZvZiIjoXVrPRetFNsL3ssDl36nSmgwWyfBXXc5qbHpyHAhWGKi3r/r5wf/+x907Sp7eoUQwsN47yKsHMi0Fm9m4uNhzhxUYiIohQ0fZic8wNBJwSxd6vrgCxksBEs4CvXr61q9AA8+CE8+KcFXCCE8zG3ZA3ZkG1G6Vq6E554jLsmCYffHQiJ2iz99p7ekeCcXNjoNNy0Eu3JFlwj86isIDobatd3bGCGEEFlyWwbgdGvxZjT/e/EiDBwIEydyrlBFnrg8Hr/C+RjZLoJaL4dS3KxavQCrVkGPHrp0YN++MGKEVC8SQggPd1sGYEhnG1F6du6EsDDUiRN8H9SPV2M/oPuz+Rg1CgoXNjHwpvD31xWeVq3SKSWFEEJ4vNs2ADtk7VpYsYLLIU3Z69eEF+1vca5II/6cCa1amdgupeCPP2DbNhg+XJda3L7dadWehBBCuJ4E4LQoBYMHw0cfYTd88LX786qxlNABjRg6FPLmNbFtx4/DSy/p9JsNGuiSiwEBEnyFEMLLyLv2rQ4ehLAw+PBDlFL42G34kcj0vhF8+qmJwVcpmDwZatSABQt0jePISB18hRBCeB0JwClsNhgzBlWrFomr1vGVf3/iCMRm6Fq9wU+Gmtu+48fhtdd06s1t23SNY18ZwBBCCG8l7+ApDIO43+ewJbAV3c59TaUW5ejctwvlD0RAaKjbCj3cxGbTQ81duuitRevX6x6wDDcLIYTXu70DcEQEfPghSS+9ypj9nfh0658k+udn1DcGvXuDj49r0kg6ZMcOePZZ2LABlizRWbdq1TKnLUIIIZzu9g3AEyfqPbNKoZatZJZaQbOHmjJ+vO5smiYxUe/j/egjKFQIfv3V5CXXQgghXOH2C8CXLumEGhMmoEiu1avsTOoeQa2pTc3P2NixIyxerHM3f/EFFC9ucoOEEEK4wu03mfj++6iJE5mXv9v1RVa+gf6EvBJqXvC9dg2sVv3zG2/An3/C1KkSfIUQIhfLtT3g8KiY66kma/jG83rDkjQKvZdhse+zRj3BmeKNmfZRJA2vRpi3yApg+XLo3Vt/DRqkywYKIYTI9XJlAA6PimHQrO1UP7iDdzaG0/xgFNEFq1DJ2EjsuWK88WYxhg2DfPlMXGR18SK89RZMmgR33mneBwAhhBCmyJUBeNTCPbTYsZKv54zEohQ2DL44O4D4UtdYt64ADRtmfh+pe9AOFWvIiuXLdfGEEyegf38YNszk9FpCCCHcLVcG4FI7NjN27qf4KAWAHR+qV9hI5CNBNGzYPtPbp/SgU8oVxsTGMWjWdgDnBOF8+aBYMZg1Cxo1yvn9CSGE8Dq5axFWQgIAx0vczV95wognACsWkiy+bLm3DMFFHUvbOGrhnptqBQPEWW2MWrgne+1SCn7/Hd55R//eqBFs3izBVwghbmO5oweckAAffYT6/XfG9trM5gn30dXenPsb/U6rgMWsLx9CdIVajAir6tDdHY+Ny9LlGd/ZcXjxRZg7Vwfc+HgpniCEECLzAGwYxvdAR+C0UqpW8mVFgGlABeAQ0E0pdcF1zUxDZKTOZFW0qN4vGx3N30WeYvC7SYQ96EOnvqeYsqUME2O7USYokBFZmMMtExRITBrBtkxQoOPtUwq+/x7efFN/QBg1Cl5/XfI3CyGEABzrAf8AjAN+SnXZQGCpUmqkYRgDk39/2/nNS0dkpE7NGB+PUorLeUvS3edvNvq2Y9I0ePRRMIzS9G5XOlt3PyCs6k1zwACBfhYGONiDBm4UT7j7bl3F6K67stUWIYQQuVOm46BKqZXA+Vsufgj4MfnnH4HOzm1WJiIidMpGpVAYfH7tBYr3aMeuXdCtGzlOqNG5XjAjuoQQHBSIAQQHBTKiS0jmPWi7HWbP1r3f4GBYt06veJbgK4QQ4hbZHQ8tqZQ6AaCUOmEYRgkntilzoaFYffzBlkiS4c8Dn4fR+HXnnqJzveCsrXjes0cXT1izBhYuhLZtpXiCEEKIdLl8JZBhGH0Mw9hoGMbGM2fOOOdOmzYleuxSFjUbDkuW0vh1E5NYJCXBJ5/oOr27dsGPP0KbNua1RwghhFcwVPJe2QwPMowKwLxUi7D2AKHJvd/SQIRSKtMJ0gYNGqiNGzfmsMkeplMnnbu5a1cYNw5KlTK7RUIIITyEYRiblFIN0rouuz3gucDTyT8/DczJ5v14p8TEG8UTXngBZsyAP/6Q4CuEEMJhmQZgwzB+AyKBqoZhHDMM41lgJNDGMIx9QJvk328P//yjVzaPGqV/b98eHnnE3DYJIYTwOpkuwlJKdU/nqtZObotni4uDIUPgs8+gdGk95yuEEEJkk2SFcMSGDfDUU7BvHzz3nO79FipkdquEEEJ4Ma8MwC6tVJQWw9B7e5cs0QlAhBBCiBzyugDs8kpFKZYs0Xt6hwyBhg0hOlrSSAohhHAar6sI4PRKRbe6eFEPM7dpA7/9Bleu6Msl+AohhHAirwvATq1UdKv586FmTV1E4a23ICoK8ufP+f0KIYQQt/C6bp1TKhWl5exZeOwxqFBB53Nu2DBn9yeEEEJkwOt6wAPCqhLoZ7npsixXKkpt9Wq9wKpYMVi6FDZtkuArhBDC5bwuAGe7UtGtzpzRPd7mzSE8XF/WuDHkyePsJgshhBD/4XVD0JCNSkWpKQXTpsErr8ClS/Dhh9Cxo3MbKIQQQmTCKwNwjvTtCxMn6mHmKVP0oishhBDCzW6PAKwU2O1gsUCHDlCpErzxhmwtEkIIYRqvmwPOspgYXTJwZHK9iI4dYcAACb5CCCFMlXsDsFLw3XdQo4Ze3RwUZHaLhBBCiOtyZzfwyBGdzWrRIrjvPh2I77zT7FYJIYQQ1+XOHvCpU7B+PYwfD8uWSfAVQgjhcXJPD/jgQZg3T28vatgQjh6FAgXMbpUQQgiRJu/vAdvtMG4chITAe+/p3i9I8BVCCOHRvDsA798PLVvqXu+998L27VCypNmtEkIIITLlvUPQcXHQrBkkJOjqRb16gWGY3SohhBDCId4bgAMDdSarOnUgOJtpKYUQQgiTeG8ABmjf3uwWCCGEENni3XPAQgghhJeSACyEEEKYQAKwEEIIYQIJwEIIIYQJJAALIYQQJpAALIQQQphAArAQQghhAgnAQgghhAkkAAshhBAmkAAshBBCmEACsBBCCGECCcBCCCGECSQACyGEECYwlFLuO5lhnAEOO/EuiwFnnXh/ZpLH4nlyy+MAeSyeKrc8ltzyOMD5j+UOpVTxtK5wawB2NsMwNiqlGpjdDmeQx+J5csvjAHksniq3PJbc8jjAvY9FhqCFEEIIE0gAFkIIIUzg7QF4ktkNcCJ5LJ4ntzwOkMfiqXLLY8ktjwPc+Fi8eg5YCCGE8Fbe3gMWQgghvJJXBGDDMNoZhrHHMIz9hmEMTON6wzCMr5Kv32YYRn0z2pkZwzDKGYax3DCMaMMwdhqG8Voax4QahnHRMIwtyV+DzWirIwzDOGQYxvbkdm5M43qPf10Mw6ia6rneYhjGJcMwXr/lGI99TQzD+N4wjNOGYexIdVkRwzAWG4axL/l74XRum+H/lbul81hGGYaxO/nvZ7ZhGEHp3DbDv0V3S+exDDUMIybV31H7dG7rMa9LOo9jWqrHcMgwjC3p3NbTXpM0339N/X9RSnn0F2AB/gUqAf7AVqDGLce0B/4GDKAJsN7sdqfzWEoD9ZN/LgDsTeOxhALzzG6rg4/nEFAsg+u94nVJ1V4LcBK9b88rXhOgBVAf2JHqsk+Bgck/DwQ+SeexZvh/5SGPpS3gm/zzJ2k9luTrMvxb9JDHMhTon8ntPOp1Setx3HL9Z8BgL3lN0nz/NfP/xRt6wI2A/UqpA0qpROB34KFbjnkI+Elp64AgwzBKu7uhmVFKnVBKbU7++TIQDQSb2yqX8orXJZXWwL9KKWcmi3EppdRK4PwtFz8E/Jj8849A5zRu6sj/lVul9ViUUouUUknJv64Dyrq9YdmQzuviCI96XTJ6HIZhGEA34De3NiqbMnj/Ne3/xRsCcDBwNNXvx/hv0HLkGI9iGEYFoB6wPo2rmxqGsdUwjL8Nw6jp3pZliQIWGYaxyTCMPmlc722vy+Ok/2biLa8JQEml1AnQbzpAiTSO8bbXBuAZ9IhKWjL7W/QULycPp3+fzlCnN70uzYFTSql96Vzvsa/JLe+/pv2/eEMANtK47Nal244c4zEMw8gPzAReV0pduuXqzegh0DrAWCDczc3LimZKqfrAA8BLhmG0uOV6r3ldDMPwBzoBM9K42pteE0d5zWsDYBjGu0ASMDWdQzL7W/QEE4A7gbrACfTw7a286XXpTsa9X498TTJ5/033ZmlcluPXxRsC8DGgXKrfywLHs3GMRzAMww/94k9VSs269Xql1CWl1JXkn/8C/AzDKObmZjpEKXU8+ftpYDZ6mCY1r3ld0G8Sm5VSp269wptek2SnUob6k7+fTuMYr3ltDMN4GugIPKmSJ+Ru5cDfoumUUqeUUjallB34lrTb6BWvi2EYvkAXYFp6x3jia5LO+69p/y/eEID/Ae4yDKNici/lcWDuLcfMBXomr7ptAlxMGVLwJMlzJt8B0Uqpz9M5plTycRiG0Qj9Gp1zXysdYxhGPsMwCqT8jF4ss+OWw7zidUmW7qd5b3lNUpkLPJ3889PAnDSOceT/ynSGYbQD3gY6KaWupXOMI3+Lprtl/cPDpN1Gr3hdgPuB3UqpY2ld6YmvSQbvv+b9v5i9Ms2RL/Rq2r3oVWjvJl/2AvBC8s8GMD75+u1AA7PbnM7juBc9bLEN2JL81f6Wx/IysBO9ym4dcI/Z7U7nsVRKbuPW5PZ68+uSFx1QC6W6zCteE/SHhhOAFf0p/VmgKLAU2Jf8vUjysWWAv1Ld9j//Vx74WPaj595S/l8m3vpY0vtb9MDH8nPy/8E29Jt3aU9/XdJ6HMmX/5Dy/5HqWE9/TdJ7/zXt/0UyYQkhhBAm8IYhaCGEECLXkQAshBBCmEACsBBCCGECCcBCCCGECSQACyGEECaQACyEEEKYQAKwEEIIYQIJwEIIIYQJ/g9Z65HzW5RAkAAAAABJRU5ErkJggg==\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": {},
   "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.2834419617286096e-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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[145.49268198]]), p=1.2834419617285855e-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": {},
   "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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[1.22491119]]), p=0.30318644106317366, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res3.f_test(R))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[1.22491119]]), p=0.3031864410631729, 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": {},
   "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": {},
   "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:                Thu, 05 Nov 2020   Prob (F-statistic):           4.98e-10\n",
      "Time:                        07:28:38   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": [
      "/usr/lib/python3/dist-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\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": {},
   "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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56240.868936151135\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": {},
   "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": {},
   "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": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2./len(X)**.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "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.113589 -31969.158503  \n",
      "1    25343.938289 -29760.155888  \n",
      "2   107465.770565 -29593.195253  \n",
      "3   456190.215132 -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": [
      "/usr/lib/python3/dist-packages/statsmodels/stats/outliers_influence.py:693: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return self.resid / sigma / np.sqrt(1 - hii)\n",
      "/usr/lib/python3/dist-packages/statsmodels/stats/outliers_influence.py:733: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_internal * np.sqrt(hii / (1 - hii))\n",
      "/usr/lib/python3/dist-packages/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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
