{
 "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 statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\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:                Mon, 24 Feb 2020   Prob (F-statistic):          2.83e-239\n",
      "Time:                        22:49:06   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",
      "Warnings:\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:                Mon, 24 Feb 2020   Prob (F-statistic):           6.30e-27\n",
      "Time:                        22:49:06   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",
      "Warnings:\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+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzddXiW5RfA8e+zd0mnCKNRCUkBZSAyGiUFC7DA+KEoIoqAoiCogGCAhKCECAYojFJJBwijQxpJZXRsMLaxun9/HEaOsXhzO5/r2rW99Tz3u8F7nrvOsYwxKKWUUsq5vFzdAKWUUio70gCslFJKuYAGYKWUUsoFNAArpZRSLqABWCmllHIBDcBKKaWUC3g782SFChUypUuXduYplVJKKZfZuHHjaWNM4ZQec2oALl26NBs2bHDmKZVSSimXsSzr8K0e0yFopZRSygU0ACullFIuoAFYKaWUcgGnzgGnJD4+niNHjhAbG+vqpjiUv78/xYsXx8fHx9VNUUop5QZcHoCPHDlC7ty5KV26NJZlubo5DmGM4cyZMxw5coQyZcq4ujlKKaXcgMuHoGNjYylYsGCWDb4AlmVRsGDBLN/LV0oplXYuD8BAlg6+ybLDe1RKKZV2bhGA3cnAgQMZMWLELR8PCQlh586dTmyRUkqprMjjAnDI5nDqDV1Gmb4LqDd0GSGbw517fg3ASiml7MCjAnDI5nD6zdpGeEQMBgiPiKHfrG2ZDsIff/wx5cuXp0mTJuzZsweAb775htq1a1OtWjU6dOhAdHQ0q1evZu7cufTu3Zvq1auzf//+FJ+nlFJK3Y5HBeDhC/cQE5943X0x8YkMX7gnw8fcuHEjP/30E5s3b2bWrFmsX78egPbt27N+/Xq2bt1KxYoVmThxInXr1qVNmzYMHz6cLVu2UK5cuRSfp5RSSt2Oy7chpcfRiJh03Z8WK1eu5NFHHyVHjhwAtGnTBoDt27fTv39/IiIiiIqKonnz5im+Pq3PU0op5eYiIyFvXqedzqN6wMXyBaTr/rRKaYXy888/z+jRo9m2bRsDBgy45RaitD5PKaWUm+vfH8Kdt67IowJw7+blCfCxXXdfgI+N3s3LZ/iYDz30ELNnzyYmJoYLFy4wb948AC5cuEDRokWJj49n+vTpV56fO3duLly4cOX2rZ6nlFLKzW3YAB06wNq1crtvX3BitkKPCsDtagQypH0VAvMFYAGB+QIY0r4K7WoEZviY9913H08++STVq1enQ4cO1K9fH4DBgwfzwAMP0LRpUypUqHDl+U899RTDhw+nRo0a7N+//5bPU0op5YaMgdBQaNYMateGpUvhwAF5LDAQ7rjDaU2xjDFOO1mtWrXMjfWAd+3aRcWKFZ3WBlfKTu9VKaXcUqtWsGABFCkCvXpBt26QJ4/DTmdZ1kZjTK2UHrvtIizLsiYBrYCTxpjKl+8bDrQG4oD9QBdjTIT9mqyUUkrZyaVL4OcnP7doAY88Al26QEDm1g9lVlqGoKcALW64bzFQ2RhTFdgL9LNzu5RSSqnM27QJqlaFadPk9muvwauvujz4QhoCsDFmBXD2hvsWGWMSLt9cAxR3QNuUUkqpjElKghEjoE4duHgRirtfmLLHPuCuwM92OI5SSimVeUePwnPPwZIl0L49TJgABQu6ulU3ydQqaMuy3gMSgFvuv7Es62XLsjZYlrXh1KlTmTmdUkopdXtr18Lq1RJ4f/nFLYMvZKIHbFnWc8jirMYmlaXUxpgJwASQVdAZPZ9SStlTyOZwhi/cw9GIGIrlC6B38/Jp39JoDGzfDnPmwKFD8Mwz0KCBZFJaulQ+8AsUkO8FC15dAKQcJzYWVq2Cxo3h0Udh/364805XtypVGQrAlmW1APoADYwxHl194MyZMzRu3BiA48ePY7PZKFy4MADr1q3D19fXlc1TSjlAcmGX5NzyyYVdgNSD8IULMGAAhITAwYNyX5EiEnwBdu+WxA7X8vWF6dPhscfs/TZUsmPHJOhu2SJ/l6JF3T74Qtq2If0IBAOFLMs6AgxAVj37AYsvp3FcY4zp5sB2OkzBggXZsmULILWAc+XKxdtvv33dc4wxGGPw8vKovCVKqVtIrbDLdQH44kVYuBCio+HppyFHDvj1V6hcWbImtW4tH/bJqlSBzZvh7Fk4c0a+b98uC4EAFi2S7EvPPw/Fijn+jWYHGzdC27Zw7hz88MP1fw83d9sAbIzpmMLdWb7kz759+2jXrh0PPvgga9euJSQkhGrVqhERIdudf/rpJ5YsWcK3337LiRMneOWVV/j333/x8vJi1KhR1En+D6eUcjtpKuwycya88ooE0urVJQDbbDK06X2Lj84cOeS5t7JsGQwbBh98AC1bwosvwsMP3/p4KnUzZsjFTKFCMvyc2u/eDbnVX71nTxlBsKfq1eHLLzP22p07dzJ58mS+/vprEhISbvm8Hj168M4771CnTh0OHTpEq1at2L59ewZbrJRytGL5AghPIQgXyxcgvdbXXoMff5RUhTNnwuUUtUDmguXQofDCCzBxIkyZAnPnQnAw/Plnxo+Znf39N9SoAbNmyVSAh3GrAOxuypUrR+3atW/7vCVLlrBnz9WaxOfOnSMmJoYAN9jorZS6We/m5a+bA4ZrCrvs3QuzZ8PgwTLMbO/e6d13SyAePBjmz5cFXQBxcdIDuf9++54vq7l4UeZ5K1eGQYMgPt5jF7m5VQDOaE/VUXLmzHnlZy8vL65d7H1t2UFjjC7YUsqDJM/zJq+CvivA8InfQWrXaAEEyspmR/eofHxk4VCy8eOhRw/ZvzpsmEf26Bzu8GGZ7z15EvbtkyF/Dw2+4GHVkFzJy8uL/Pnz888//5CUlMTs2bOvPNakSRPGjBlz5fYWe4+jK6Xsrl2NQFb1bcTBFjlZPOV1avd//fqVzc7WpYv0uH/4Ae65B0aOhFSmvrKdv/6S0YGDB2HSJAm+Hk4DcDoMGzaMFi1a0LhxY4pfk9ZszJgxrFq1iqpVq1KpUiW++eYbF7ZSKZUmcXHw5pvQsKH0RleuhDJlXNeeXLlgyBDYtk1WTffsKb3h7M4YGDtW/k5580qSjRY3lifwTFqO0Imy03tVyq0ZI/tyZ82SBVdDh8I1U04uZ4zsNS5WDB54AM6fl7lON83o5FDGyFB9QoIUVMiXz9UtSpdMlSNUSqksx7LkQ71xY6mM426S25esTx9ZGDZ2rOQ2zg7++0+Cb8mSMizv7w9ZLBdD1no3SimVmuhomUsE2dfrjsE3Jd26SW+4Qwd44glZhJSVLV8ONWteHYLPkSPLBV/QAKyUyi6ioiT5RbNmcPy4q1uTPtWqydznxx9L/ulKlSTndFZjjGyHadxYhtvHjbPboUM2h1Nv6DLK9F1AvaHLCNkcbrdjZ5QGYKVU1nf+vCzcWbECvv3WI/IE38THB959V1JdVq8OZcu6ukX2FR0tRS3efFNSfK5dCxUq2OXQybm/wyNiMFzN/e3qIKwBWCmVtUVEQPPm8oH+00/QqZOrW5Q5lSpJndsyZaTH2KkTTJ4sBeg9WVKSZLYaPFjybefJY7dDp5b725U0ACulsraJEyVh/8yZ8Pjjrm6NfZ0/D0eOQNeuUK8erF/v6halT3Q0fPihVJnKlUva37+/3ed705T72wU0AANHjhyhbdu23H333ZQrV4433niDuLg4QkNDadWq1U3Pnz9/PjVq1KBatWpUqlSJ8ePHu6DVSqk06dVLPtjbtXN1S+wvb14IDZW80ocOSaKKLl0kn7W7W7ZMqkcNHAi//Sb3OSirVbF8KacFvtX9zpLtA7Axhvbt29OuXTv++ecf9u7dS1RUFO+9916Kz4+Pj+fll19m3rx5bN26lc2bNxMcHOzcRiulUhcdDR07wp49sqWnWjVXt8hxvLxktfCePfDOOxLYfHxc3apbi4iQKlCNG0vbQ0PhyScdesrezcsT4GO77r4rub9dyDMDcFiYZIwJC8v0oZYtW4a/vz9dunQBwGaz8cUXXzBp0iSio6Nvev6FCxdISEig4OUN8X5+fpQv79o/olLqGomJ0Lkz/Pwz7N7t6tY4T548kkN6zx7InVsSd7Rte7V36S5eeUV67H36yJxvgwYOP2W7GoEMaV+FwHwBWEBgvgCGtK9yfe1nF3C/RBwp9SafeEL260VHyzzH33/LhL2XF1StCm+8ITUhT5+W7DbXCg1N9XQ7duygZs2a192XJ08eSpYsyb59+256foECBWjTpg2lSpWicePGtGrVio4dO+KVBfeoKeVxjJFVtCEhkku5bVtXt8j5/P3l+5EjcgHSsiU88ogEvPr1ZUTA2davlwuE8uVlK9Xbb8s+XydqVyPQ5QH3Rp4XNSIjr672S0qS25lgjMFK4R/kre4H+Pbbb1m6dCn3338/I0aMoGvXrplqg1LKTr78Er76SoJwjx6ubo1rlSkjeaVHjIDVq6WnWbEiHD3qnPPHxMjq7Nq1ZW7600/l/rJlnR583ZYxxmlfNWvWNDfauXPnTfelavVqYwICjLHZ5Pvq1el7/Q0WL15s6tevf919kZGRpkCBAmbBggWmZcuWqb7+1KlTJleuXGk6V7rfq1Iq7eLjjXnwQWM6dDAmMdHVrXEvFy8aM2WKMR07GpOUJPd9/70xy5dfvW1PgwYZkz+/MWBMpUrGjB5tTGSk/c/jAYAN5hYx0fN6wEFBkgFm8GD5HhSUqcM1btyY6Ohopk6dCkBiYiJvvfUWzz//PDlSKHcVFRVF6DXD2lu2bKFUqVKZaoNSyg68vWHRIvj++yyZtjBTcuSQhVo//CBD0ElJ8MEHV3vFn38O27dLrzW9jJHpv99/vzo6efEiNGkCf/4px+3e3a77erMKrYYE/Pfff7z66qvs3r2bpKQkHnnkEUaMGEFYWBgPP/zwlQVXAD/++CNDhgxh//79BAQEkDNnTkaOHEmtWikWu7iOO7xXpbKcffvgvfdgwgTZlqPSJjpa9kaPH391QesHH8i+3IgIybp1zz1w993y3d9fKhHlzi37qocNg/375St5KnDJElndbIxr5prdkFZDuo0SJUowb968m+4PDg4mJoUrwvr16zujWUqp2zl9Gh5+GM6dg1OnNACnR3Kv+LnnZLHWli2SZQtkAdePP0ogvtbPP8ui2NhYeX65cjIKWa6cvDZ5RbMHBt+QzeEMX7iHoxExFMsXQO/m5R2+aEsDsFLKM8XEQJs2UrZu2TK46y5Xt8hzVahwfd7lypUlmceZM7B3r3xduiQLqkB2o+zd65q2OkByrujkdJXJuaIBhwZhDcBKKc+TlATPPgtr1sgwat26rm5R1mNZUKiQfGXx329quaI1ACul1LWOHpXgO3y41Mh1AVcMWSrHcFWuaLcIwCaVPbdZhTMXuymV5RUvLntcXTTn66ohS+UYxfIFEB4Rw33hu6jz7zbWlKzCpsCKDs8V7fK1+v7+/pw5cyZLByhjDGfOnME/OUONUipjli2T4gqJibIi10UX7u5a3k5lTO/m5XkofDszp/eh18ppTP/pPYJO7HV4rmiX94CLFy/OkSNHOHXqlKub4lD+/v4UL17c1c1QynPt3i3DzYGBUkHHhftK3bW8nUqn/fth7VraderE3fkisEwSNoDEBPrnOsm9WX0VtI+PD2XKlHF1M5RS7uzUKcln7OsL8+e7PKlD8pBlSvcrNxYWJgmcLEvqBCxZAgEB0Lo193ZqAxNHQlwc3r6+ctvBXB6AlVIqVbGxUsv32DH50Cxd2tUtonfz8tfNAYN7lLdTqQgLk2I/cXFyu0gRGDQIunaV5CLJWRZDQ+V5mcyymBYagJVS7m3zZkn68P338MADrm4NcHWhla6C9hDx8RJYEy9fMHl5weuvSwa1awUFOSXwJtMArJRyb0FBcOCA9FjciDuWt1M3OH8e3nkH/vlH6gf4+koP2NcXGjW68rSkJOkg//yz/Dx6tHOapwFYKeWepk2T7EsvvOB2wTfZ2bMyMn7HHVCgANhsrm6RuuK33+B//5M94z17Qq1a1w0xmzpBrF8nQXfmTEmo5ucnJeWdlcpaA7BSyv0sWiRzcw8+CM8/7zaRLTFRasv/8QcsXAhea1bTgD/5k0as8wqicGEJxnfcIdcMtRPCqBUVSvn/BVO4jfOGNrO1iAgZXp42De69F3755crUhakTxBb/IH7+GWZ0hoMHwccHWrSAIUMks2nu3M5rqgZgpZR7Wb0aHn1Ukvv/+qvLg++xYxJs//gDFi+GqLOXaEgon+QfTzAhYFkkeH/MwobDqLtqOBejc3HhQC4S4xOpFLcFA8T95ke/Bktp9F4QjRtrtUSHCAuT3m3t2pIlbcAAqejk6wvIgudevSR/i80m1RLff1/W9+XP75omawBWSrmPrVtlu1FgoEQ9V30yXm7KuGfDyP93KKEEc77wXfyW5zXui/kdn5gLEOUDGDAGn6Q4Wt2zF4o1pkBUFERFcXHrdryOJWEBhjiqrx7NwWaTeSawMzV61Of5rl4UKuSyt5e1hIZCs2YygevrK7WJL1dmOnIE3noLZsyAsmWl+mL79rjF714DsFLKfSxfLnt8lyxx2bxvbKys11k59C+WJDXCmwTw84dfFuL1ynZ4+ikZq8yVSy4Wkhf1dOp0ZQVtyOZwfh41g0nT+uKTmECCzUZclQs8ty2El8O/4d8+JZjybkfOtujMU60vUuVMKFbD4JtW4GaXfNOZep+HD8PTT8tKZ5C/x+rVxAU14MsvZadRYqJ8791byhq7C8uZKSBr1aplNmzY4LTzKaU8xLWrXiIjXZbjecUKeOklsPbu5q+czSl08V95wGaTqNyv3/UvSB72vGHfaL2hy1LMLVwuByy9K5ILX08nx18LOUQpippj+BGH5e+L17Kl1wXxlPYaD2lfxa5B2NVBPlPvc/Fi6NhRSlMmJEik9fVl/bClPDs2iN275Vrpyy/BVfmeLMvaaIypldJjOhOhlHKtkyel3F1YmNx2QfCNjIRu3aBxg3hePDWEHT7VKeR1Tnq2Npt8Dw6++YVBQRKUb+i5Jqek3BRYkbFBT7ApsCIAB6KBTp3IvWIBthPHKPHyI/hbcdhIhNgYtvb6joQEOYYz8k0nB7/wiBgMV4tKhGwOt9s5bifD73PiRFk9VbSo7BMPDSXy7cH0D1rK/T2CiIuDefNgzhzXBd/b0SFopZTrRERA8+awZ8/VJAlpYM9eW0gIdO8Ox4/Dl0+t5/Wf3oXHHuP3/73H7F9XctfODeyrVItH/EvSLo3HTFOqysKF8X2uE0ydiLkUh5WUSLU141l4RwTFpg6zW77p1H5XrqqDe60Mv8969WSF/KhRmBw5+W7V3bw2OojERPjwQ9n+607DzSnRAKyUco2LF6FVK9ixA+bOlS1HaWCvUoDnz8sW47m/XKJrmVC6rmlO7dp14Z1NhHCHnCNvWRYFlQVgZTrOkeZUlZfTH1qhoZj7H2Dn18t56JfhWK1D6Fn8Rb7s8DDmhiCSnnzTt/tduUNRiTTn1U7OlHHuHEyZAhUqwMSJREXBq89JorTgYOkYly3rlKZnmg5BK6WcLy5OKhuFhcEPP8hQYhrZY2j29GlJhFRg1jecyVmSsYcfoXahg/JgjRqZPke7GoEMaV+FwHwBWEBgvoBbz2leHsa2Gjei0swPidu2h613P8bdRw5wZEoj7l13nFfDZnBf+K5055u+3fu4VTB3ZlGJ3s3LE+Bz/Vazm95nch7nkSNh6lQZW0ZWqtesCdOnS693yRLPCb6gPWCllCvYbJAjB3zzjaQeSofM9trCw2XHyrN73uWdpCFYF5E53uPHr0wW2qNnmNFUlXkrl+CBvdNYvTyelp3W8OOfr+FPLHHevqz5ZiYN0nHM270Pdygqcdu82gkJMp6cXETBZsNs38HX4W14803JQLZ0acpT9O5OA7BSynkOHJBgV7y4ZCjKQEaKzJQCPHBAEjA8Gf45fRKHXH0gMVFWM19eTOUO5QbrNvDhl1f+wuv9OCzANyGO6uO+gs4PS/qmNLjd+3CXohKpXqw8/jj89Rd4e4MxGF9f+i8J5pM/ZeBk6lQoXNipzbUbHYJWSjnH2rVQpw48+6zczmA6qDQNWaZgxw6ZZo6MhG7dbTIGHRCQ4irnjJ7D3rwbB+MV4IfxspGEjbzrlnCi9P0kbdqSpten5X20qxHIqr6NODi0Jav6NnK/fcYvvSQTvCtWcOR/g3ks31KGrQhi2DBYsMBzgy/oPmCllDP8+qskSyhWTJLkl89cIEvvKugNG+Dh5klUsP3DuGXlqVwZ2Xu8Zs0t67+6en/sFZf3Gp+/L5jxHx7nmbBXWF+8PbXXj+XOO2//crd5H2llDIwbJ0POPXteuevLL6FPH9l19NNPTq0amCmp7QPWAKyUchxj4PPPJQXRAw/Iamcnd1lWrIC2LRP4xrxIe37Fa8d2KFXKqW2wF2Pg+5FnebufD+TOzewPNlMvYJPspXZSEXmHio6GV16RceW2bWH2bM6ctWjxaAwbVgaQ4+7jVO64l37tyrn3RcQ1UgvAOgeslJvwuJ5KWsTGypaR9u1lGDHAeXOoIJ3tju0v8YtvR5pGzZZ8hCVLOrUN9mRZ8GzPAtRqBk89BbbXX8GwFiwLy99fViN5ahD+5RepYnT8uCxp7t+fVast2j2WwJlTfuRvsoPc9x3ixCUytO3MHekcsFJuwB0yEtlVVJSkBwwIkCHeGTOcHnxXDAtjfasPWWWrT9MLs2ULy/vvO6fQq4NVqiRT6mdrNcdgYRmDiYmRoVtPtHixLLY6fhx8fUlq3JShn3rRoAFciIunyNOryFPz0JU/nb0zgrmKBmCl3IAz0g46zbx5ULWqZCkCKFjQ6fX3/p4QRq2+jXnffMi90evhvfegRw+ntsHRAgLgkVEtSPL1JwEvDBZ8/z2JM391ddPS7vRp+b5hw5ULI5OYyHddQunXTwZOijy7Ar87z9/00mu3WIVsDqfe0GWU6buAekOXecyFqwZgpRwpOho2boSdO+X26dMyRPjYY/Dmm/DZZzBjBt4HD6T4cmdmJMq0AwegdWvJfh8QIPkdXWDPHpjTMxRf4iQseXlBzpwuaYvDBQXhHbqUmHc/4r16y/kfX9N0VBsOHED+zSVXCHI3iYnwxRcyF794scxf+/uT5GUjNtGXKYeCGTdOEl8VL5LylqvkrVSePHqkc8BK2dvAgZIcfvt2CUrGSG9w8mTpDebPL3ti/vhD0jECnZp1ZUj+9txx4QyzpvVmXYl7WVWqOvurPuDSt5JmISFSlcbbG0aMkN5mGveq2tOJE9CyRSKfJO7Ey88HErh1IYWsIiiI3EFBfGJg6tT6/NQD6lS5yCHvxgQUL4A1bpz8LW6x2tvpdu6UHKBr1kgq0kqViC0YyPdPLeXglFD2lwhm1LwgqlWTp98uWYg75LPOKF0FrZQ9hIdLEXmAKlWkMPi990LlyvJVo8bNJVmMkWIER47wx9FLvPnXafKcPcH7yyZS59+/KRQdKc8rX14yRtWv79z3lBYXL0rv8uhRePdd+Pjjq78HJ4uKguAGhhf/fp1uCWNkz0revO4RdJzo33+ha1cIWDqPbwNep0jMYdnrbAz4+bluoVZYGHz0ESxaJH+XUaOgY0cWL7Ho3h3++QeeeQbGjIHcua9/aWoLFMv0XUBKUcwCDg5t6fC3dTu6ClopR9m5UwLP4sXyCVKsmCSoTcucp2VJbzh/flpUgdg7whm+0JfX2/YhMI8fg8oZGoVvkw/M5A2f06bBV19B06ZSRahOHZf0NDlwAHr1grNnYflyed9Tpji/HZclJMATT0DTzcPoZsbA22/D0KEua48rlSwpMW7s2Nbc27sRs71b8mDCciyQvbXXZPxyikuXJPg+8oisivfygu++42iNlvTqKMPMd98tbW7aNOVDpJYpyx2ylmWYMcZpXzVr1jRKZQn//mtMly7GeHkZkyePMR99ZExUlOPPO2uWMXXrGmOzGQPG5M5tTNu2xsTEOP7chw8b8+mnxtSpI+fOkcOYoUONiY93/LlTkZRkzIsvGvMM30m7OnUyJjHRpW1yF3v2GPPivatNNP4mES+T4BdgzOrVxnz+uTE//eTYv93hw8b062dM4cLGtGlz5d9sks1m/mr5icmd2xg/P2M+/DBz/3xnbzpiKvT/3ZTqM//KV4X+v5vZm47Y771kArDB3CIm6hC0Uul16hSULi3drtdek4LshQo5tw0REfDnn7BwIRw6JPPJAG+8Ib2cBx+EihVl+DqjC5CMkd58iRIydz1pkszd1awpy1OffVZyOt/A2fuZBw+GDz9I4N8itSlWuaBs/vX1ddj53ElaftcJCTC3Xxh7JoQy93wwOYIfYM6R+8i17/LftkcPWbW+cWPGh+svZ+siOFgWHo4eLUlXQOZ5W7aEnj0xcXFcSvKloVlK3uZBjB4Nd92V2d+Ce++h10xYStnDqVNXszh9+62Ml7lbRqVnn4XZs2VCNFnHjlLyD2S8r1gxCZzR0XDhgryncuVkeHDCBLnv+HFJtHvwoHyYdu8uBXTPnUv1Pd9YfxZkwcwtS/Fl0pQp0KWLvO0pX5zDsnnJ/GI2kN7fdWwsjB8PQ4bAyRNJvF9jAW9Zn5NnU6g8wctL5ojnzZO/cdmyV6dSrg2w1wZoY+D33+WCLCFBLnzKlYNjxySHc7duUKoUf/8NC/qHcWFeKNsLBfPsuCA6dMgSW7JvSwOwUpk1Zw507izZetJRu9Yl4uJg717YvRt27ZJezvPPy5aUnDlv3prSs6dsCYmOvtpbDgiQD9v27SUlYBrTR9YbuizF+bjAfAGs6tsoc+/rBosWQfdHDvLFncNotvNLfPP43/5FWUhGf9cXL8LYsTBsGJw5A/NKv07LQ2OwMLJYq1MnyVqWK5csKCxSRC7GkpIkwH7wgcz7Hz4sX9HRVw9us8n2usGDOXXBnx9+kIukLVtkqcJrr0mSqxsXWWVlughLqcyYPBlefBFq15Yvd+fre3X19bW8vWWh2K5d0kPJmRPy5JEeC0jQPX1aPh0zOIRrjzq6abF9O3zd9nfW8Sz5Ii9hnXgL8txt13M4Q2aGTjP6u86ZU1Jzd+smC5FHDe1EIybiSxxJli9bS3bg3k+5CawAACAASURBVJEPkWPf37B1K4m//Y7t8kVbQuwlTixZSeC5kzLF0aKF9ILHjoXERIyvL8sLtOeLJ/357TfpFNeseWXBs9NnatydBmClUjN8uBQDb9ZMKvrkyuXqFmWcZcnQ4q2GkC1L5nozwRkrUk+ehMFNQpkZ2wovkrDi/eTC4W7PCsA3DiEnJ5CAtOU4zuzvOnduSRDWvXsQM/suJfaPUKaHB7Py4yAsC6pVg+IVo7DqzmTGym54JyUQb/PmzcItaPj649QseifnzslCeF/bE/isCmXszmAWvhvEnXfKwMpzz918Haiu0iFopW5l6dLL1duflOos2WRhT2Y4eg740iVo3MgwNqw6Vc3fcqfNJiux+vXL9PGdKbPD9Y74XcfESH6MFStklHn5ykSSEmzUIYyGtqX8mdiYNaS8SMvPT2YrnntOrle9r+neufMiKUfTIWilMqJRI5g5Ex59VD7k1W0lf6g64sPWGFnXc3T1QSoEHIQEn6vzkh6Y6Sqzw/WO+F0HBEDDhvIFUPrtRcQez8vufwuyI7ozXn7x5A/Ygc0/nsndqpM/PxQowJXvKV2jZrann5XdNgBbljUJaAWcNMZUvnxfAeBnoDRwCHjCGHPOcc1UykliY+HVV+GttyST1WOPubpFHie1pAmZMXSorA368MOy+L68V5KBLF/usZmu7DFc76jfdbLAQn6Ee5/Dv/j1H++B+QJo1Sptx/DkVJGOlpZiDFOAG5d99gWWGmPuBpZevq2UZ7twQbL1TJ4s43DKbcyaBb+9u5Kfqn7C+/2NZAarW1eGnT0w+ILkOA7wuX5k5docx+7AHm101sI8T3TbAGyMWQGcveHutsB3l3/+Dmhn53Yp5VxJSbJMc8UKSff4wguubpG6bNMmGNB5H/O8H+Xx2KlYF6Nu/yIP0K5GIEPaVyEwXwAW0qt01H7pjLJHG2/Vo3e7VJFnz8pCSydK0yIsy7JKA/OvGYKOMMbku+bxc8aY/Ld47cvAywAlS5asefjwYTs0Wyk7++ADWcgzZowMQSu3cPQoNKl5jnmn61Amzxm81q6xT+ok5TTOTs6SLqdOyQS2t7fkD//sM9lDb8cV9aktwnJ4PWBjzARjTC1jTK3CadzMr5RTJSbCunVSQuaVV1zdGnVZdDS0bxXHuFMdKGMdwitkdpqDr6cWaM+K3K6nf+yY7Ftu3FimMlaskPu7d4f16516gZfRVdAnLMsqaow5ZllWUeCkPRullFPZbJLpJyEhe+TG8wBJSfBxqzBe2DyJ+l4r8Jo0Jc3lGHXVrftx9GKxNDl6VPZILV0qS+rLl5dKZmXLyuNlytxcMtTBMtoDngs8d/nn54A59mmOUk4UESFJhI8dkyDs5+fqFinks3FUxzDe+7MxL1iT8fL1uZqtKw1SW3WrsqGkJPlesKBsdP7gA0mltmuXTDuVLu2ypt02AFuW9SMQBpS3LOuIZVkvAEOBppZl/QM0vXxbKc+RlARPPw0//gj797u6Neoaw4fDnTO+xJ9LeJlEyV0dGprm1+uqWwXIiNZXX0H16pIA288PVq6EgQNli6EbjHbddgjaGNPxFg81tnNblHKegQNl2Hn0aCndp9LMkVmNpk6FP/v8zjx+xfICLFu6E204rUD7pUtXR002bIDNm6XigI+PtNnHR8rw+fjY97zq9pYuldKcO3ZI1bKICEmC7QZB91qaCUtlP7Nny9BTly664jmdHDm/+scf8E3XMBZ7dcCralWs4Z/Koph0Jtro3bx8iqtu7bK/9uRJuXCbO1fKMe3cKbm1582DQYNufv758xKAR4yQkkDt2kkBA0/OKe7OLl6Ued5ff5X53JAQaNPG7QJvMs0FrbKXpCSpaOTtLVmU/O1Xwi475Lt1VLnBdevg1QY7WBpfn1ylCmJb/ZeUwcsgu/8t/v5bVsiHhckkdYkS8sHeu7cE4KgoiIyUUpDx8fIVFycVDby8JDiPGiX1//z9pVfWsaN8KfsxBlq3liQtvXrZ9f93Rmk9YKWuFREhizGKFrXbId16r+M1MhuYyvRdQEqfGBZwcGjLDLVp716oVw/eTRhED99x2Nasdvpq1BStWiUf6A8+CCdOyHBy69ZScaBatfT3qhIS4K+/ZAQmJATuu09+Bti4UW67aU/NrUVESOrYgQPlwsgYt/o9unQfsFJuwRhJMRkbC/ny2TX4gmesvE2+SAiPiMFwdfg4PXtk7Z3V6NgxaN5cPi9br3sf299bXB98z52Dl1+WwDtggNxXpIjM8w4YIIt6MvIB7+0tw+kjR8KhQ1KpHmDbNhmVqV9fgr6Hceme6/Xr5cJl6tSrvzs3Cr63owFYZQ/TpkmijalTHXJ4T1h5a4+LBHvmL46MhPbNohj5X3uWfbWDu+62MjXsnGnGyKr4ChVg0iTJjDR3rmPOZVmQN6/8XLEifP21FJd48EGZJ9650zHntTN7XNRliDHwxRcydJKYKMk0nnrKsed0AA3AKus7cgRef10+3ByU49kT8t3a4yLBXlmNYmPhsTZxDNrRnlZmLpVzHkzX6x1i3jzo1En2hW7YIPuhcuZ0/Hm9vaXH/c8/8PHH8OefUKeOFAdxcy4b+fniC5njfeQRWX3uoQU5dBW0ytqMkaCbkCBDfg6q6+vQlbd2Yq/tOZnNahQZCe8Fr2LElleoxjaYNJk017azt7g46W1Wry5t+PFHePxx19R/zplTMjO9/LIMrebOLf9+v/5aLgySe8xuxOkjPwkJcsHy4ouQJ4/83/agIecbaQ9YZW3jx8t2kREj0pVNKb3cLt9tCtyh/N3x49Cr5nK+3NJAgq+Pj6QEdIW9e2X+sFEjuSrw8pJhTFcE32sVKgQPPyw/b9okW+UqVnTccHgmOG3kxxgZkahTRxZQ5skjQdiDgy9oAFZZXYMGMlT1v/85/FTtagSyqm8jDg5tyaq+jdwq+ILrLxL27ZMpu1KHV2Dj8khBUlK6slzZzfr10pgTJ2RdgBv2LgGoWVPaWriwrL7u1AlOn3Z1q65wykVdZCS0bw/vvCML9BIS7HdsF9NtSCprcrOtCNndpk3wXLNjxCT5MefTPdzbo7EM//r6StYiZ87hLVwIHTrAHXfIz3YsPecwcXEwZAh89JGkUdy82W3+fTt0//vff8vf6tAh6QG/8YbbvO+00n3AKvv57DPpOUyZ4hab8bOzpUuhd5s9zLnUnAJBFci58g9JaBEamu4sV3bRtasEsN9/l3J0nuTvvyWZR8OGkuzj3Dm5kMiKjJGEGocPw4wZHpsyVgOwyl527pS5vYcfhlmz0nTFnB2yWLnCjBnwVec1zDWtyJPfhu2P32RY1RUiI2WoOT7+6jyiJ/vkE7nQHDVKhqY9rGd4S7GxMsycKxccPAg5crh2e1omaSIOlX3Ex0su2Ny5ZfVoGoOvS/YyZnFjxsD3T85ncVIj8pTMJxmuXBF8k5IkU1KtWnD2rCz88vTgC/Doo3DPPVLVq21b2W7n6Q4dkoQkL74ot8uU8ejgezsagFXWMnSo7OEcNy7N/3E9IYuVJ7lwAT57LIxjr33EJP9u+Fa/V4KvA1eh31JcnNR8/vxzGRHJl8/5bXCUihUlteWIEbBkicwN//qrq1uVcb9dHh35559skyNb9wGrrOPiRRg7VobjHnsszS9z5F7GpCTYs0emo9evh+ilYZQ/FsrOO4I5XiaIQoVkgWuhQlz5+c47ZQTd06aujZEUxz+8tIzvzrTCz4rDCx+sT79zzTxlVJQs4Fm0SIZr+/bNOsO0yWw26d23ayfFIkqWdHWL0i8yUnYqTJoEVavKRcRdd7m6VU6hAVhlHTlzyuIaX990vcye9WNPnpSseOvXS4WfjRuvJjRq6/cHM+PaYjPxJEb6MDpmKKGJ9fnz/D0cjZIh0TqEEUwo/XyCsT0YRKNGst6mdu20vS1XzWUfPizJxpg3l2lezxJALJYxEI/8Ihq7oHx4r16yAmzSJCk9mZWVKycXGsl69pQRoLffdv96xDExMH8+9OsHH3zgeVeemaCLsFTWsH69DF95pX9WxR6VjC5cgGHD4K/hYdSNW8a/trKULZlIg9ybuNTyUUp1rk/FL/+H17cTbn7xvHlcatqKmE++IO+gXgAkWd58XeQDPjr+Ese5k5w5ZRFockC+776b80W4oiJTQoLUFpjQ/19GxPegdeIcTJmyWEfD5UFXbDNKduYMrF0r6Qqzk8RESSjyyy9Stenbb2X+252cPStrNPr2lf+zFy7Iuo0sSFdBq6xtxw6oUQP695cr6AzIaM8xMVF2OvXvDzWOL2C+1QbLJHFloNPfX+YfX3lFtt20aCGBycdHIlfRovDAAzJE+9prMoR+w//JJZ9tZc7Bquxa+C97/rEozhFa5gglqUEwFbsE0ayZLO51VK3eW1m7VvKblN4awo+2p/HzNXh9OFB6Xxs2ZGqbUYZ78slVrzp3Bj+/dJ83S5k9G7p3l2QjvXrJ/w13CHJz5kC3bnDqlMxh16nj6hY5lAZglXUlJkpGo/37ZftR4cJOO/WyZfK5tntrLPcF+fNz1Y8pMb6/POjlJQH1s88kd22y1Pa/hoXJUG1ygoqvvpL5sR495BivvQZjxmCwMEA8PrRkPsu9m/LQQ7CZnQSUO4lPgYvXHTYztXpvFBsr73vD6DAu/b6MbYUa0X1QEZot64M1YoQUp8+kTPXkBw6EDz+ECRPgpZcy3RaPFxEhGaSmTZNqS3fe6bokNadPy7/lH3+UnvnkyXLhnMVpAFZZ18iR0uOaPl0WXznBnj3QuzesnneaYbk/5infX8lxcCfW9m3XB9CMDL2mFqB37ZITL1hw5a6EnHl5/7VzzF9gcW77EY5SjHq5/6RxroWElazCtoqlKXVXAmHvNczw+z19Wk45dy4s/yOG56NHM5R+kk4yIAArne/zdr3bDPfkx42TvMldusDEiVlvwVVmHDsmoy3GyJB89epy9ejEC1aaNZNKT++/L0PP6Vyr4ak0AKus6eBBqFxZgtX8+Q7/wI2Nha+fCyNy5iJK2I7S2fYTvvFRWF27yvanggUdn+Hp2l6yt7d8mL33HhhD1B3FiYuIIk/CRSySiMeXZixkje9D1KhuUbMmV74qVJDDJSTIIEJCwvU/R0bKmp45c6TOeeOkRbzp/zWNEhbilxCNQXrW2GwweLAsoEmDtPRuy/RdQEqfSqn25GfOhCeflIpGs2ZdP+qgroqOlkxgM2ZAQIAMBb/9tgRnezJGhpe//hpGj4b8+eWCtHBhWemcjWgAVlnTxo0yzDhnDpQo4dBTRUVBv4ZhfLqhEf7ESvB56CH5gKlY0aHnvklKQT4hAX74gciPhpLnn11X5qAPlryfMY+vZev6OGqv/Yq1l6rjQxw12EIowaxBXm+RRA6i8SUOPy4RzJ88zXTmlelBkWea8+L5zyk+83Ostm0ld/K772aop5+W3m26e8AXL8oq4LvukquGHDnS1JZsbfdu2Zr1ww9ysfL777K6L7Oio+WYo0fD1q2y7zokRIqiZFMagFXW5YT5rLNnZdSu8bohDLbexyspUeZ4P/oozT0/p7m2h2yzyQVCly6wfTtUqXLlacn/61c1Hcj6hwdQ8vBKOox86Obj+fpKsL/vPvk5+XedwZ5+Wnq3GZoD3rlTenH586e5LQpZOzFypIzg5MghQXnBAlkYFRQk34sXT/m18fEytB0XJxc/Z8/KxdnZs9LLff11mRbK5hdEqQVgHadRnuf4cVmg9N57Dv/PfewYtG0cRa89L1OuZyu8vva92vMLDnbouTMkKEh6pDcGx8qVZdVpnz4webLs0bUsHrz7JA++CYSXheLD5X0tWwbz5kkWkcREOdaNQTYoKEND7GnZc50cZG+7CnrfPpl6eOMNqFQp3W1RyMjBqFFXbxcuLBdZY8bI6n2Q+YqdO+X+Hj1gzRpJe3n8uFwAN2ok/+YKFJD1GMHBsmdO5+BvS3vAyvM8/rgEiK1bHVrM/eBBeLLRKcb+25KabMSaOhXKlnVdFR97uHGldUrDx2l5TgbZba/y6dPSO4uMlN59Fs4X7BJxcfL/a80aOH9eLnZBcjQfOQKBgdIzLl5c/g8+lMLoiQJ0CFplJSEhkoT+449lHtJBdu6EFxoeYNrp5pT2Ccc282do3dph53OqtAwfO3AxWaazdV26BE2aSPKV0NAsv49UeTYNwCpriIiQocY77pAPXwel2Fu/Hl5rtpf55+uTP3cC3r/P98zeblZkjFS7+v57+OknWfmslBvTcoQqa+jXT7L6TJzosOAbGipTWpF5S+Lfuhnea1Zp8HUnGzbInu9BgzT4Ko+ni7CU5+jRQzLnOKim7PqvwjjacwxPFX+Ggauakzvwe4ecR2VC7doyRJENMiiprE+HoJX7S0hweGKFY7PCKNihAb7EY7y9sVas0J6vO1m7VoorZLfCCsrj6RC08mzdu0uB7qQkhxw+JgZWvzQZH+IBZItOaKhDzqUy4PBhaNNGtrjExbm6NUrZjQZg5d4WLZLE+iVKZKjU4O0YA+91Okjjsz/LvkWbzX33+GZH589LeslLlyTjWTbJH6yyB50DVu4rMlL2HVaoIItuHGDMGGgW8gp+/l5Y3/0kmYE8dY9vVpOQIHVtd+2CP/5wfspPpRxMA7ByX2+9BeHhsHq11NW1s5Ur4c03oVPTyTR7/wDUr2f3c6hMmDlTchSPHy/7fpXKYjQAK/d06pQk3XjnHSlYb2dHj8K4Nr9TrnQzRs0sildeO1eDUZn31FNSv9YeRQKUckMagJV7KlwYduyQaip2FhcH4xv8wA8RnTn26lfkzfua3c/hKJnOIuUJZs+WaYeKFTX4qixNF2Ep9/PHH1IEoEgR8POz++G/6LiOfvu6cqriQxQd8LLdj+8oyXmUwyNiMEB4RAz9Zm0jZHO4q5tmP8uWSYINB6YZVcpdaABW7mX+fHj4YfjmG4cc/ufPw3lmVjui8xal8IpfPWpV7fCFe64rYgAQE5/I8IV7XNQiO9u6Fdq1g3vugcmTXd0apRxOh6CV+zh7Fl5+WerWduli98Pvmria+9/qTAGvCLxD10GhQnY/hyMdTaGMX2r3e5TDh+XCK29eGQFxwNSDUu5GA7ByH2+8IYuv5s+3+9BzzLIwSr/UBF8u4eXtjRVzwa7Hd4a01NL1WIMHS0aUv/66dQF4pbIYHYJW7mH2bJg2TeqO3nef3Q+/os98fEwcNpKwkovMe5jezcsT4GO77r4AHxu9mzuuJrLTjB4Ny5fDvfe6uiVKOY32gJV7KFoU2rZ1yOKbNbOOUmfDV5JJy8JjM10lr3Z29Cpop620TkyUnu8bb0D+/FC1qv3PoZQb02IMyrWMkRSQDhJ90bCh8MPUjl2BNXky/kcPaKarVCSvtL52sVeAj40h7avYNwgbIzm+x42D776DZ5+137GVciOpFWPQHrByrTfflPzLI0Y4JBAvaDWOx2MWsrfnWO55TuvH3k5qK63tFoCNgbffluD7zjsafFW2pXPAynV+/RVGjpShSAcE3w3T99Ay9G12lmzBPZ93s/vxsyKHr7ROSIAXXoDPP4fXXoMhQ+xzXKU8kAZg5Rr79kHXrnD//fDpp3Y/fHQ0fNN7L2e8i1Bq6SSHDnNnJbdaUW23ldZnz0qyjQEDYNQoh1S4UspT6BC0cr7YWHj8cRl6njHDIckw3nsPJhxrTafFLShxl49djpkd0kD2bl4+xTngTK+0joqCgAC44w5JuJE3byZbqpTn0wCsnG/jRti7F37+GUqVssshrw2OdfaeJvfsS3R/9SUaNLFf8L02MCWngQSyVBB2yErr06fhkUegVi0YO1aDr1KXaQBWzlevHhw8KL0hO7g2OPpFx/HhnI/wt2JZ2y4YuMcu53DK4iQ30a5GoP3e05Ej0KwZHDgA779vn2MqlUXoBIxynl274Pvv5Wc7BV+4Pjj2/ulX7kraT8/gfny58YjdzpGl00A6yt69crF15AgsXAitW7u6RUq5Fe0BK+e4eFHmfU+elA9iO+b6TQ6Czy36gxdPTWdG/tZsur8Ulh2DY5ZOA+kIcXHQooWklwwNdUh2M6U8nfaAleMZA926wc6dMH263RPtF8sXwAMHtjFg82gM0PrCQu4L32XX4Gi3NJDx8XD8OFy4INuvsproaHlfvr4wYQKsXKnBV6lb0ACsHMsYSbYwbRoMGgRNm9r9FL2bl6fyutMYvLAAn8QEHgzfcV1wDNkcTr2hyyjTdwH1hi5Ldw3ddjUCGdK+CoH5ArCAwHwB6csO9csv0KCBLEAqWhTy5AFvbwlYICkZK1WC2rXld9S3L8yd61lBeuFCqFxZEmwANGkC5bNAnmqlHESHoJVjrV8vWa66d5e9QQ5Q0VaIz4605y1rDL7EkeDtQ81nH6XB5eBorxXMaVqcdPAgLFoEYWGwejUsWAB33y3bcGJjpdzi3XfL0OzFi7I1B6QC0L33yn0nT0qiim+/lepQABMnShKLunXlee60f/bkSejVS0Y3KlSA6tVd3SKlPILmglaOt2yZ5F92QNAwURf57477GMUbvDuzBgX+Dr0p13O9octSnL8NzBfAqr6N7NOQ06flAmPiROm1Fi4swfKjj6RXmF4xMRLMK1WS2/XqSUAH6T03bAhPPSVfrjR7Nrz4ogypv/su9Otn91KSSnkyzQWtnO+XX2Sl80MPQSM7BbkU7H6sPxVj9vLQ21Uo0DIIWt5cZMEpK5iNkdSar74Kr78Od911U/atdCXyCAi4GnxB6uTu3y8965UrZbi3UCEJwMZIoG/QQC48fOyz9zlV8fFynnz5pJ3jx1/fXqXUbWkPWNlf8paT4GD52UFpICMXriF3i7rMLvIKjx4dc8sOtkN6wHFxEnQWLoR58+Q9RkVBrlwpPt3uVYaMkV5yjhyyx7Z8eRmizptX5pCbNIGWLe1X3D4xEdatk/c6f74E+vHjr7ZFU30qlaLUesBuNJGksoTVq6F9e5mnnDHDcR/Mly5x8akXOEJxys8akurotl0L2SclSQavSpWgRw9ZRHXunDx2i+ALqSfyyBDLkuALULYsnDkjPfDHHpO/QbdusGmTPL55s8zR/vyzDGun96K7b19ZOFa3ruTtLlgQHnjg+rYopdJNh6CV/fz9t/S6AgPhjz/svt3oWtvHr6J8xF6mtp/DC3XzpPpcu6VX/O8/CXDr1knx+N9/h+bN0xSAHD4MniePXPi0by8B9tpMYzt2yMrkL76Q24UKybz0r79CgQJSj/eHH66+D8uSedzZs+VnY6RH3bq17O3Nn98+bVYqm9MArOxn4kTImRMWL4YiRdL8svQWOYiPh47fNCJXsf0smVoyTeewS3rFfPnk/X33HXTuLMUk0sipiTwsS3rFyZ5+Gp58ErZtk4uHdeskS1WyS5cgMlJ+Nka+bDY4dgyKFYNhw+zfRqVU5uaALct6E3gRMMA2oIsxJvZWz9c54CwqMVE+sJOS5EM7MO2BLt1zo4mJ/PDaajp/XZ85c6BNG3u8gdtYuhTq1JHgm8H5TrvPASulPIJD5oAtywoEegC1jDGVARvg4j0RyunmzYNq1eDECdlmlI7gC+mfGz07cBSdvn6IPg+FOSf4jh4txQQGD5bbGZzvzHQiD6VUlpPZIWhvIMCyrHggB3A0801SHiEpSba+DBggqQbj4zN0mPTMjZr9B8j5yXsssLWm+/d1MnS+NEtKgt69JSFG27Z2qeRj1ypDSimPl+EesDEmHBgB/AscAyKNMYvs1TDlxs6fl8U+AwbAM8/IHtUMbne51RzoTfevXs3Fek1JSLIIf3csJUo6cOVtTAw88YQE3x49ZLFSzpyOO59SKlvKzBB0fqAtUAYoBuS0LOvpFJ73smVZGyzL2nAqOa2e8mx9+she0C+/lAVJARlfSJSmLUJhYZjgYHKdOICvFU/Xpv9l+HxpcuKEXFR8+SWMHJmuxVZKKZVWmdkH3AQ4aIw5ZYyJB2YBdW98kjFmgjGmljGmVuHChTNxOuVyCQny/eOPJb3kG29keg9omuZGQ0Mxl+eJva0kvP8KzdQ5byl5JXDp0rJK+I03HHMepZQic3PA/wJ1LMvKAcQAjQFd4pwVJc/3Ll4sK4ILFJAUk3Zyu7nR6bFVeRQ/fIgj0ebN2qL30sBuZ7/s3DlZ6fzYY3KBkSf1vcVKKZVZmZkDXgv8AmxCtiB5ARPs1C7lDoyRZBO1asl8b5kyEoydaO3HY1n1xX6a+f/G53WfodOTH9Ftv2+6ywmmKi5O5rQPHZJEE0op5QSaC1ql7OhRSfS/cqUE3sGDoVMn56YdPHaMyJIV2JZQhcfafo5/hatrCOxWycgYeOEFmDwZvv9eklYopZSdaDUklXbnz8vwa+HCEpzGjJFyc76+zm2HMVx87hV8E+J4tdQw/Mpfv4DPbikchw6V4PvBBxp8lVJOpQFYib17ZZh5xQr45x9J9L9ypcuaY378iZyL59DbNoyzj8TifUPH224pHCtUkB7wwIH2OZ5SSqWRBuDs7Nw5CAmRLUVz5kgC/jffdPo8700uXCDuf6+zmQeI6PkCuQPWE3NNno8MVzK6VnS0XGQ8+qh8KaWUk2k5wuwkKUkS8e/eLbf37YOuXWHNGkk4ceCArHZOpayeM5yNz81ztml8XnkyXw8raP8UjocOwT33wI8/2qvJSimVbtoDzurWrYM9e2DJElnRfOqU1IodNw5q1pRasdWquU9N1+ho+vTJwS9RLdg4TXJg2DWFY2QktGoFFy9CjRr2OaZSSmWABmBP988/Uvv1wAHYv1++BwbCqFHy+JNPSo+vQAHZYvPII1LDFqR4QvXqLmv6TU6f5lKl6nBqAL16v0S1anY+fny8pJjcswcWLpT5X6WUchENwO5u/nwIC4Pjx69+5ckDf/4pj7/0EixfLj/7+kod2Gszjv34I+TNK0Ou7pxSMSyMpP+9gtep4/xXrA5fDnDAfI0xqgAAIABJREFUOd56CxYtgm+/hUZ22MKklFKZoAHYHZw+DVu3Xv06eFCCqmXBjBkSRIsUka8774S77rr62iFDpGdXpoz0fL1umNav4+CqQfYQFgbBwXjFxZGANwPfjrJ/7QNjoGBBePttWfWslFIupgHY1d59V4JosmLFZE724kVZDPXVV7JP9Va916Ag57TTkebMwcTFYQE2y1AnNhSw8/uyLNlmpZRSbkJXQTvbpk0ybLx5s9x+6in49FPJs3zyJISHw2+/XV2JnDevew8d20Gifw4AErDh5ecLwcH2O7gx0KWLDD0rpZQb0R6wM8TEyFDyuHGwdq2U76tfX1bhVq0qX9nYx14fsJI6jHhyI9XeCLZvr370aJgyBWrXhmbN7HdcpZTKJA3AjpaYCPfeK/O6FSpIfdlnn4V8+VzdMtf76y/2rotg0KBWPNmpGdWm2zlAbt8OvXtDy5bwyiv2PbZSSmWSBmBH+e8/KF5cho/ff19qzAYHu89+W1c7fZqkJ5/C+1ROShZtzpgxPvY9fmysFI/ImxcmTtTfu1LK7egcsL0ZAxMmyLafadPkvi5doGFDDQLJLs/LJh4/xWPxPzFxqo/9BwSmTYNt22QBW5Eidj64UkplnvaA7en8eXj5Zfj5Z2jaVOccb2XkSJg/n16MotFbNWjY0AHneOEFKF9e5tqVUsoNaQC2l40br2ad+uQT6NPn5j25Cvbvx7zzDgv92rD87tdY/7Gdj3/6NFy4IPuiNfgqpdyYBmB7+e8/iIuTBBr16rm6NW7LlCnLyMrfMmx7Sxb9YOHnZ8+DG6ldvHq1pOR0cVEJpZRKjQbgzIiIkIDbti20ayc5lgPsVKc2q1m9GubNY77Vhjc3P8uIEVClip3P8e23Ulbxs880+Cql3J5ljHHayWrVqmU2bNjgtPM5VEyMBNyNG+HffyXNoUrR8ilzqftiB7wTE7iEHy9WWMjUHQ3sO0J/4IBE9Lp1pdCCDv8rpdyAZVkbjTG1UnpMP6UyIjEROneGv/6SVbYafG8pZHM4J0aOwzsxQVJNkkD1/NOYuzXcficxRha/2WwwaZIGX6WUR9Ah6PQyBrp3h9mzZTXvE0+4ukUuF7I5nOEL93A0IoZi+QLo3bz8lfq9v373O1/vCMVgkYAX8V7erC5XntkL99ivxu+lS7Lo6rHHoEQJ+xxTKaUcTANwei1cCOPHyyrnHj1c3RqXC9kcTr9Z24iJTwQgPCKGfrO2AdDuTi+GT+zDeVtuOsdPp9od69nWrCCbAitiRcTYrxH+/vDNN/Y7nlJKOYGO1aVX8+Ywb971FYyyseEL91wJvsli4hMZvnAPFC7M4nLNaZ6wkN+LNmLy0w+xKbAiAMXy2WGxmjHQty9klXUFSqlsRQNwWv3+u+QWtixo1UqzWl12NIWebEBcLHFHjnLkuDc9j37PrhwVuKP9Brx8kuRxHxu9m5fP/Ml//RWGDYNlyzJ/LKWUcjINwGmxejW0bw9vveXqlridG3uy3okJjJ47jFk/9KFDy1jiY22M+OY8JYt7YQGB+QIY0r5K5ud/z5yRufiaNaFXr8wdSymlXEDngG9n1y7p8ZYocTW3s7qid/PyV+eAjeGThaNpvH89Q+4axYbt/ixYAC1aFOEN7JyPuVcvOHtW6vx66z9jpZTn0U+u1Bw5InO+fn6y+KpwYVe3yO0k92SHL9xDpzlf88S2JfxYuR/vbn+dMWOgRQsHnPTPP2HqVOjfH6pVc8AJlFLK8TQAp+bDDyXb1YoVss1FpahdjUDazR4Pa2ZyoEpbOm37mJ494dVXHXTCevVg1CjZ+6uUUh5KM2GlJj4edu92QM7ELCYsDBo3xsTEEoM/H9RbyrDlQdhsDjhXXBz4+jrgwEopZX+aCSu99u6V+UUfHw2+qblwAfr1g8WLMZfisDD4EsfHTUIdE3z/+gvKlYOtWx1wcKWUci4NwDeKjYVHH4WHH5Z9pipl+/ZBnTowfDinYnMTa3yJx4aXvy9+zYPtf77YWKl0ZLNJEFZKKQ+nc8A36t8fdu6URVe61zdlf/wBHTuCzca2EQsJHtyY+/PUYeIzoRTrFAxBQfY/56BBsGeP/F200pFSKgvQAHyt5cvh889l9VCzZq5ujXuaNEl6olWrMrfLbB5/pwylS8Po34IoVs4BgRdgyxb49FN4/nn9uyilsgwNwMnOn5cP+HLl5MNeXRUWBqGhEBwMDz6I6foCwwO/pE/PnDz0kNSlKFDAgeefOhUKFZI6v0oplUVoAE526RJUrgzvvgs5c7q6NXaTWqWiNAkLg4YNZfWxvz/xfyzl5cRvmDIInn4avv1Wtkk71GefQc+eDo7ySinlXBqAkxUuLEUWspBUKxWlJQhv2wYvvSQXJ4CJi+O7LqFMORDEwIHwwQcOniY/cEAWXZUqBSVLOvBESinlfLoK+tQpePxxOHzY1S2xu1QrFaXm+HGpc1y1Khw8CN7eGJuN2CRfvjsczNSpMGCAg4NvUpJMCdSvL/uxlVIqi8nePWBj4H//gwULpDt3jUwP3bqBlCoV3Xh/yOZwfpswi7t2buDwPdVo+uqTtLsnH2zeDP37Y3q+yV8T97BiUCjLrWAGzwsiONgJjR8/HlaulEVfPj5OOKFSSjlX9g7A06bJCqJPP70u4Uamh27dRLF8AYSnEISTKxiFbA5n5pc/Mmn6u/gkxsMKi2fikqBnR9rt3s2GzTbe7gDLlwdRsWIQs2ZBhQpOaPi//8I770CTJtILVkqpLCj7DkGfPAk9esCDD95Uzi7DQ7c3CNkcTr2hyyjTdwH1hi4jZHN4ppudHr2blyfA5/qUVFdq8a5YQcHOTzBx+rv4JcZf/odgqHlgM4N/OkznZ23Uri1boseOleRTTgm+xkC3bjIEPWGC7sVWSmVZ2bcHPHQoREXBN99wY97EtAzd3o479KLb1QjEFhPNhpGTCdr8JwWS4rj44Uc0rBEIczdS4thBVpS5j+ADG/EyScTbfPjtTDs2ffYAO31kQXifPpAnj1OaK+LipPTjkCFaAEMplaVl3wA8cCA0apRit+52Q7dpkVovOjkAp2me+do9uDdmmIqJkfnrxYuhWDEICJD0kA89BPv3w/330/rsWVpf+5r/tgIPQ+vWdO6Ti/CIGGr8t5saa/9j4X/tWb2rPoXvO87GkKKUKJHmt2s/fn4y/6uUUllc9gvASUmQmCjdulatUnzKdUXmL7sydJtGt+tFp6mHvGgRtGkjvUKbTS4YOnSQMnxRUZA7980neP99CcBFishK5gMHYMkSed8225X81pHnLYKozjfzY5m7vykhsb74lzpN6aar+eLV0leCr1MXo73/PrRtC7VSLByilFJZSvabA54+XYq4Hz16y6e0qxHIkPZVCMwXgAUE5gtgSPsq6Qo8t+otJ9+f3EO+L3wXr4bNoNZ/Oyh7ZC+Lvp4pT0xKkuB76ZIEzYQEWLcOjh2Tx3PlgubNr86R2myyN2jQoKuPjxsnPX0/P7DZSPLx5ecTwTRpcjmxVL8CEF6EgpXOUOSJtdT439988Wrp63ro/WZtIzwi5v/t3Xl4VeW59/Hvw4wgRETB4ICIpcixFU3VaNXYI7V1OAynUlAUfUGxKg5HESiVQaiK4tBajx4VbfUocklRZosMQVBUkEFBQEKOKIMBZVQQSPK8f9w7EEISQrL2Xnvt/D7XlSubPaz1LFb2vvczrPvGc+BLQlzmsidOhBEjYNq04LctIpKEqlc94J074Sc/saQO8+dDjfh9/yjZwwXrRRcF8lMHTOGSNQt5YfxwahXacxzwWbPTOPObHHvB4MG2Qjs/32rgzpx58DB0rA7v/hq5xR7fssVKGa9cCd+/O5/CWdmM3ZTFh2TStq3F9quvthHrskoHXvjIrFKH4luk1ef9Ab8K5P8JgLw8W4Wenm5fMlTvV0RSRHn1gKvXEPSIEZZkYsKEuAZfODCMXNbwbXpafQZmv0TtWPAtBCb99CJeuboP44o28uCDVhaxjDngfRmZbHx5JrumZrO4cRazX85kZX8Lups3H3henTqZXHBBJl0HwKtXQ+vWFTuGIBajHZb3cNNN9uXo9dcVfEWk2qg+AXjVKnjySfuwP/fchOyyU/sWBw9bf/ABdOkLjz1Gv8vb8MKH3Rgx+UlqFRawr2YtxpzfmR7dsw7aRsG5maw+JpOcHMh5ClavtlK8OTmWvKugIBOwwNy0qa0p69jRfhf9tGxZdi+3PEEsRjusMWNs2PmZZ+CMM4LbrohIkqs+AfiZZ2yV8MMPJ3a/8+bZXOzSpbB8ORxzDKxYQaerroIhd3LnCSfS+vOF5JyRwe9v6ULHs1rwxRe2bmrGDJg9G7ZtO7C5xo3h9NPtO8S111pvtnVraNPGAnCQgliMdljXXGOL4nr0CG6bIiIRUH3mgPPzLQD+/OeJ2+esWZbNyXtbLHXPPTasXKLaUl6eTd8WBd2vv7b7TzkFOnSwdMht2ligbdIksbkp4rYKevdu+OGH4L81iIgkkeo9B7xnD+zaZT3PRAZfgI8+OnC7Rg0LNrHg670F3WHDrJMM1sRf/coSYHToAK1alR9sE3GJ0CHD6EHp3x/Gj7cvRY0bB799EZEkl/qXIT31lK18/uabxOxv1y7o1cvGjrOyoF49m4CtU4eiKgZz51qJ3Q4d4Msv4c9/hgULbOHUuHGWifG00w4ffBN2iVDQpk6Fp5+24WcFXxGpplK7B7xhAwwfbsPAzZvHf39r1liijKVLbUHRvfdaNze2ivnjmpk8cLnl12jeHP76Vyu3W6/eke+qIpm2klJeni2E+9nPEj8fLyKSRFI7APfvb3O/TzwR/31NnAg33GBDzVOmwBVX2P2ZmSypn8ngwTBpko1CP/YY3HYbHHVU5XeXkEuEglZ0ydGOHTY/XplvHiIiKSJ1A/D8+VZucNAgm0yN1z6ysyEtzSLqOefYGHLLlgBs3Qp/+AOMHWtPGTHCCjCVlkHySCXkEqGg/fCDjauPGgXt2oXdGhGRUKVuAJ40ycZ5BwyIz/ZLZqHq189WOMd6dTk5lmo6Nxf+9CcbjU5LC273CblEKGgNG8LkyWG3QkQkKaTuIqyHHoIlS+xDPx7efdcupSkosCB8zDH7g+/cuZbicfNmu6xo+PBggy8Ek686YbZuhRtvtBVnzqnGr4gIqdgD3rcP1q2zWrLNmsVnHwUFtrgKbM632ArnV1+F3r1tFHrKlIqnfayMuF0iFKQ9e6BzZ8sCduON+4fnRUSqu9TrAY8ebVkrli+Pz/a9h7vvhvfes98jRsDMmRSel8kDD9g6rAsvhA8/jG/wjYTCQlt0NWcO/P3v+7+kiIhIqvWAd+60knznnx+/vMKjRsHf/maTuqNGATYSfdO1ttiqVy/47/8OpqZAQmvxxsOgQZbr+eGHLW+miIjsl1oB+PHHYdMmuyQoHvOM+fk2rty1q5UJxC5r7dTJkl6NHGlrsYLYdclyhkWJNoBoBOHvv7fz0KePXQ4mIiIHSZ0AvHGj9UivuQbOOy8++6hVC955x27XqEFOjuX42LQJ/vlPm+oMSmQTbRRp2NDmfRs00KIrEZFSVGkO2DmX5pwb55xb6Zxb4ZzLPPyr4mTOHJtzfOih4Le9bJlVsN+yxVY616vH9u12mdH339t0cJDBFyKaaAMsp+aNN8KPP1qayVqp8x1PRCRIVf10/Avwjvf+d865OkAVcjtVUbdullz52GOD3e66dfDb31pw//57aNKEggLo3t0yT86cCRml1rmomkgm2sjNtW8lRx0F27cr05WISDkq3QN2zjUCLgZGA3jv93rvt5X/qjj5/HP7HXTwffddi65btlgBgZNPBqxa0bRpcPJVK+k5dQoXPjIr8CII/S5vQ/3aNQ+6L6kTbXz3nX1R2bfP/nPidQmYiEiKqMoQdCtgM/Cyc26xc+5F51yDkk9yzt3inFvonFu4efPmKuyuDPPmWVrDsWOD3e4HH8BvfmOrrAoKrMoRlt3y0Uch7Zy1FLRZE7dKRJFKtLFjB3TsCGvXwoQJ8NOfht0iEZGkV5Uh6FrA2UBf7/1Hzrm/AAOAB4o/yXv/PPA8QEZGhq/C/g7lvS07Tk+Hq68OdNNMnmzDzmCrn7Oz+bhmJr17Q6NWW2l06cHXGcdjgVQkEm2ALYBbscKykFx0UditERGJhKr0gNcB67z3RVXnx2EBOXHGj7eMFw8+WLXSQqW5+mqoX39/Ld/N7bLo3BlOOAEaX7kAV/PQ7xJJv0AqaJ9+al+C2rSx+d9rrgm7RSIikVHpAOy9/wb42jlXNCn578DngbSqIvbts0IL7dpBz57BbbewEJ57Dtq3txVWw4ezZ+pMrvpzJtu326WtJ6WXPnCQ1AukguS91VRs3956vWArnkVEpMKqugq6L/BabAV0LnBT1ZtUQcuWwbff2qRsJS51KTPL1LPPwh13WEDp3h1/fiY394SPP7YO95lnQr/8CFYiCsquXZbseswYS0jyn/8ZdotERCKpSgHYe78EiMNFOBXQvr1V12nU6IhfWlaWqQZf5dLh/vttNW+3boAl13r1VRvlLrrWt2heNtJpIivjq68s7deSJXa99YABSrIhIlJJ0c6SUMlhz9KyTO3Zs5dmd90HdevCiy+Cc7zzjmVRvOYaq+lbXGQWSAVp2TL4v/+zWstXXhl2a0REIi3aAbiSSlss1XvB2/xs7XJ47TVIT2frVivk064dvPxyNe7obdwIs2bBddfBFVdYAA66uLGISDVULQNwaVmmZrfKIN3t5cbu3QHr+W7aZLUXGhxydXM1sHatXfA8erQtuvr1r+G44xR8RUQCknr1gCugeJYp5+1a33XprUh7fCQ4x5w58MILcM89cHZiL6wK34YNVlOxdWv7T+jZ067xPe64sFsmIpJSqmUPuPgiqq6TXqTtD3nsfn40Hdu34McfrYJey5YwbFi47UyovXsPFDEePx5uuw3uuw9OOincdomIpKhqGYAhtohq0Tsw/w1LOXluS8AW965aBf/6V4oPPX/3naXb/OADS+dZu7bN9aanw/r1wSc2ERGRg1TLIWgA5s6Fm2+2+c3sbJg/n+XL4ZFHoEcPm/JMGYWFsHr1gX/ffjs0bWolFh9/3JKaXHKJpdwEBV8RkQSotj1gRo604Auwdy+Fs7O5eXImjRrBE0+E27QKKyyEzZtt3nbjRrj4YmjY0LrvL74I33xjPxs2WAKNvDw4/ngr23jyyXDBBVbtqX41yeAlIpJEqmcAzs+HTz6BGjXs+qI6dXhrSxbz58M//pFE640KC2H5cli61IaFf/97m5yePNnmaDduPNBrBVi4EM45x4LysmXQvDn84hf2+8wzD9Tn7dQplMMREZEDqmcArlXLgtq8ebBqFZvOyOKm6zO57DK4/vqwG4cNF/fvD++9Z3O1Rdq2tQDcvDlceim0aGE/6elWJaKoDGCPHvYjIiJJq/oF4K1bLYPW8cdDly54D7d0to7kc8+FkHDjiy/sYuPsbJuT7dXLhpGXLLGKTFlZcN55thq5aFVYRoZ11UVEJLKqVwD2Hn73O5vznDwZsCtuJkywKeHTTktgW5Ytg+HD4c03rV2tW9tqbLDebG5uAhsjIiKJVr0C8IQJdqnN3/4GwLZt0LcvnHUW/Nd/JbgtN99sQXjAAJvPPfHEBDdARETCVH0C8J49cO+9lty5Tx/AYl9entX4rURFwyOzeLF1s59+2lZ5vfQSNGsGTZrEecciIpKMUjYAl6z3+/zGGbTLzYXp06FWLZYuheefhzvvtCnVuFm0yFJqTZxoc89LlthlQG3bxnGnIiKS7FIyEUdRvd/123bjgQ1bf6DG2LFsvOTXFvyAgQMtHg4eHKdGeG893owMW808bJjVL47tX0REqreU7AGXrPfrXQ069RhF63qFTAFmz4Zp06zYT1xHgBcvhq5d4X/+p9K1i0VEJDWlZAAuXu83fccmttRvxI+16/F5vnVM+/e3NU933BGHna9dCwUF0KqVXSpUp041LiYsIiJlSckh6PS0WGpF7/nrxMd4Y8xA8J70tPqMGwcLFsCDD8YhA+OcOTbkfP31Funr1lXwFRGRUqVkAC6q93vn+2PIWL+CeSefRf06tbjnV2344x9tIfQNNwS4Q+/hmWfgssusyMHLLyvwiohIuVJyCLpT+xY0XfABF77/Oh7ovWgiZ97ag88/aUFODkyaBDVrBrSzPXusutDo0Za56n//Fxo1CmjjIiKSqlKyBwzwy/GjcYAD6hXmc/6Xyxk2DC66CK68MsAdFRTYpUV/+hO8/baCr4iIVEhK9oAB2LnzoGpHr63PIi8P3noroNHh3butG33UUfD++zbfKyIiUkGpG4DnzYOZM2HBArb8LIu7umXSpQtkZgaw7YICuO462LHDau8q+IqIyBFKvQD87bcWIJs1s0VRl13G0Dutw/rQQwHto18/60r/5S8BTiaLiEh1knpzwEOGWF3cHTsAWLPGygz26gVt2gSw/aefhiefhLvusjyWIiIilZBaATg31xI8d+u2fzHUAw9YoYUhQwLY/qRJcPfd0LEjPP54ABsUEZHqKrUC8JAhULu2RV2sDsKYMVZqMD09gO2feip07gyvvaahZxERqZLUCcCffWaBsW/f/dG2f3849libsq2S7dst2ca//RuMGwcNGlS9vSIiUq2lTgCeMQPS0izqArNm2V2DBlWxDsK2bbZ0euDAYNopIiJCKgXge+6BnBxo0gTvYehQ6wj/4Q9V2ObevdCli2338suDaqmIiEgKBGDvLUDC/tqC2dkwdy4MGAD16lVh23fdZbULR4+GSy+tclNFRESKRD8AT58OP/mJJcSIGTYMTjgBbr65Ctt95x27fum++6y6kYiISICinYijsNDmZk85ZX8Pdc4c+3nqqSr2fn/8EX75Sxg+PJi2ioiIFBPtADxuHCxeDK+8YoXvsTq/zZvDLbdUcdudOtn1viorKCIicRDdIei5c60M4KmnwrXXApb+edYsuP9+qF+/ktudMsVSTBYWKviKiEjcRDMAz58PHTpY3uf16+HjjwGb+23WDPr0qeR2t2yB3r1t0VV+fnDtFRERKSGaQ9DZ2QcCZEEBZGfzgc9kxgwYNcoqBFZK374W1KdN2z+kLSIiEg/R7AFnZVmArFnTfmdlMWwYHHcc3HprJbc5fjy8/rqlsTzrrCBbKyIicoho9oAzM63Wb3Y2ZGXxoctk+nQYObKSWSJ374bbboP27ZXxSkREEiKaARgsCGdmAjDst9C0qcXQSqlfH954wzZSu3ZwbRQRESlDdANwzMcfW86Mhx+Ghg0rsYEdO6x0YVZW0E0TEREpUzTngIsZNswyUN5+eyVenJcHp58Ozz4beLtERETKE+kAvGABTJ0K994LRx9diQ3cdpuVGlTvV0REEiySAfjtxeu58JFZXNI9j1r199Hyog1HvpHp023l89Ch0LZt4G0UEREpT+QC8NuL1zNw/GfkrqzN7jXNaJCRy/Dpn/L24vUV30h+vnWbW7WyMoYiIiIJFrkA/Ni/VrF7XwG7VqRTo+4+Gp3zJbv3FfDYv1ZVfCOLFsHq1fDoo1C3bvwaKyIiUobIrYLesG03AGlZK2nYfi016uYfdH+FnHsurFkD6enxaKKIiMhhRa4HnJ5mVRacg9ppuw+5/7C++AK8hxYtVGxBRERCE7kA3O/yNtSvXfOg++rXrkm/y9sc/sW5uXDmmfDEE3FqnYiISMVEbgi6U/sWgM0Fb9i2m/S0+vS7vM3++8vVvz/UqgXdu8e5lSIiIuWLXAAGC8IVCrjFzZ0L48ZZ5g7N/YqISMgiNwRdKYWFdrlRixZw331ht0ZERCSaPeAjlpsLX39dxWLBIiIiwakeAbh1a8jJqWStQhERkeCl/hD0woWW+eroo6FG6h+uiIhEQ2pHpPXr4ZJL4P77w26JiIjIQVI7AP/xj9b77ds37JaIiIgcJHUD8OLF8Mortvr51FPDbo2IiMhBUjcADx4MxxwDAweG3RIREZFDVDkAO+dqOucWO+cmB9GgQOzcCWvX2jW/jRuH3RoREZFDBHEZ0l3ACqBRANsKxtFHw5IlNv8rIiKShKrUA3bOnQhcCbwYTHMCkJsL27bZJUd16oTdGhERkVJVdQj6KeB+oDCAtgSjd2+44AIrOSgiIpKkKh2AnXNXAZu8958c5nm3OOcWOucWbt68ubK7q5jZs+2nTx/V+hURkaTmfCV7is65h4HrgXygHjYHPN5736Os12RkZPiFCxdWan+H5T1cfLENQa9ZA/XqxWc/IiIiFeSc+8R7n1HaY5XuAXvvB3rvT/TetwS6AbPKC75x9+67MG8eDBqk4CsiIkkvda4DnjEDTj4ZevUKuyUiIiKHFUgA9t5ne++vCmJblfboo5b9qm7dUJshIiJSEdHvAXsPX31lt5s0CbctIiIiFRT9APzWW3DaaTB/ftgtERERqbBoB+DCQhgyBFq1gl/8IuzWiIiIVFgQqSjD8+absGwZvP461Ir2oYiISPUS3R5wQQEMHQrt2kHXrmG3RkRE5IhEt9u4cCHk5MCYMVCzZtitEREROSLRDcDnnWcB+KSTwm6JiIjIEYtuAAY45ZSwWyAiIlIp0Z0DFhERiTAFYBERkRAoAIuIiIRAAVhERCQECsAiIiIhUAAWEREJgQKwiIhICBSARUREQqAALCIiEgIFYBERkRAoAIuIiIRAAVhERCQECsAiIiIhcN77xO3Muc3A2gA32RT4NsDthUnHknxS5ThAx5KsUuVYUuU4IPhjOcV7f1xpDyQ0AAfNObfQe58RdjuCoGNJPqlyHKBjSVapciypchyQ2GPRELSIiEgIFIBFRERCEPUA/HzYDQiQjiX5pMpxgI4lWaXKsaTKcUACjyXSc8AiIiJRFfUesIiISCRFIgA7537jnFvlnMtxzg0o5XHnnPtr7PFPnXNnh9GgABKcAAAEqElEQVTOw3HOneScm+2cW+GcW+6cu6uU52Q557Y755bEfgaH0daKcM596Zz7LNbOhaU8nvTnxTnXptj/9RLn3A7n3N0lnpO058Q595JzbpNzblmx+5o45951zq2O/T6mjNeW+75KtDKO5THn3MrY389bzrm0Ml5b7t9iopVxLEOdc+uL/R1dUcZrk+a8lHEcY4sdw5fOuSVlvDbZzkmpn7+hvl+890n9A9QE1gCtgDrAUuCMEs+5ApgGOOB84KOw213GsZwAnB27fTTwRSnHkgVMDrutFTyeL4Gm5TweifNSrL01gW+w6/YicU6Ai4GzgWXF7nsUGBC7PQAYWcaxlvu+SpJj+TVQK3Z7ZGnHEnus3L/FJDmWocB9h3ldUp2X0o6jxOOPA4Mjck5K/fwN8/0ShR7wuUCO9z7Xe78XeAPoWOI5HYFXvPkQSHPOnZDohh6O936j935R7PZOYAXQItxWxVUkzksx/w6s8d4HmSwmrrz37wFbStzdEfhH7PY/gE6lvLQi76uEKu1YvPfTvff5sX9+CJyY8IZVQhnnpSKS6ryUdxzOOQd0BcYktFGVVM7nb2jvlygE4BbA18X+vY5Dg1ZFnpNUnHMtgfbAR6U8nOmcW+qcm+aca5fQhh0ZD0x3zn3inLullMejdl66UfaHSVTOCUAz7/1GsA8d4PhSnhO1cwPw/7ARldIc7m8xWdwRG05/qYyhziidl4uAPO/96jIeT9pzUuLzN7T3SxQCsCvlvpJLtyvynKThnGsI/BO423u/o8TDi7Ah0J8DTwNvJ7p9R+BC7/3ZwG+B251zF5d4PDLnxTlXB/gP4M1SHo7SOamoyJwbAOfcICAfeK2MpxzubzEZPAucBpwFbMSGb0uK0nnpTvm936Q8J4f5/C3zZaXcV+XzEoUAvA44qdi/TwQ2VOI5ScE5Vxs7+a9578eXfNx7v8N7/33s9lSgtnOuaYKbWSHe+w2x35uAt7BhmuIic16wD4lF3vu8kg9E6ZzE5BUN9cd+byrlOZE5N865nsBVwHU+NiFXUgX+FkPnvc/z3hd47wuBFyi9jZE4L865WkAXYGxZz0nGc1LG529o75coBOAFwOnOuVNjvZRuwMQSz5kI3BBbdXs+sL1oSCGZxOZMRgMrvPdPlPGc5rHn4Zw7FztH3yWulRXjnGvgnDu66Da2WGZZiadF4rzElPltPirnpJiJQM/Y7Z7AhFKeU5H3Veicc78B+gP/4b3fVcZzKvK3GLoS6x86U3obI3FegMuAld77daU9mIznpJzP3/DeL2GvTKvID7aa9gtsFdqg2H23ArfGbjvgmdjjnwEZYbe5jOP4JTZs8SmwJPZzRYljuQNYjq2y+xC4IOx2l3EsrWJtXBprb5TPy1FYQG1c7L5InBPsS8NGYB/2Lb0XcCwwE1gd+90k9tx0YGqx1x7yvkrCY8nB5t6K3i/PlTyWsv4Wk/BYXo29Dz7FPrxPSPbzUtpxxO7/e9H7o9hzk/2clPX5G9r7RZmwREREQhCFIWgREZGUowAsIiISAgVgERGRECgAi4iIhEABWEREJAQKwCIiIiFQABYREQmBArCIiEgI/j+pYk3Y5LCFVgAAAABJRU5ErkJggg==\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 = sm.categorical(groups, drop=True)\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:                Mon, 24 Feb 2020   Prob (F-statistic):           5.69e-38\n",
      "Time:                        22:49:06   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",
      "Warnings:\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+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeVxUZRfA8d9lAAE3zF1w3xcU1zSzzCVNzcyyxdJ6K7VsNzV9K9s0LUwzM0stq7esLJVMy13KfQV3EXdARVRwA2GYed4/HhZRkAFmQ8/38+ED3Llz79Ngc+bZzjGUUgghhBDCuTxc3QAhhBDiViQBWAghhHABCcBCCCGEC0gAFkIIIVxAArAQQgjhAhKAhRBCCBfwdObNypUrp2rUqOHMWwohhBAus23btjNKqfI5PebUAFyjRg22bt3qzFsKIYQQLmMYxrHcHpMhaCGEEMIFJAALIYQQLiABWAghhHABp84B58RsNhMTE8OVK1dc3ZQizcfHh8DAQLy8vFzdFCGEEDZweQCOiYmhZMmS1KhRA8MwXN2cIkkpxdmzZ4mJiaFmzZqubo4QQggbuHwI+sqVK5QtW1aCbyEYhkHZsmVlFEEIIYoQlwdgQIKvHchrKIQQRYtbBGBXM5lMBAcH07hxY5o1a8akSZOwWq03fM7Ro0eZM2eOk1oohBDiZuPyOeD8Cg2PJWRpJCcSk6ni78uIbvXp0zygUNf09fUlIiICgNOnT9O/f3/Onz/P+++/n+tzMgJw//79C3VvIYQQt6Yi1QMODY9l9PxdxCYmo4DYxGRGz99FaHis3e5RoUIFZsyYwRdffIFSiqNHj9KhQwdatGhBixYtWL9+PQCjRo1izZo1BAcHM3ny5FzPE0IIIXJSpHrAIUsjSTZbsh1LNlsIWRpZ6F7w1WrVqoXVauX06dNUqFCB5cuX4+PjQ1RUFI8//jhbt25lwoQJTJw4kUWLFgGQlJSU43lCCCFETopUAD6RmJyv44WhlAL0PuWXXnqJiIgITCYTBw4cyPF8W88TQgi7ungRzpwB2YJY5BSpAFzF35fYHIJtFX9fu97n8OHDmEwmKlSowPvvv0/FihXZsWMHVqsVHx+fHJ8zefJkm84TQgi7WbwYevUCw4DUVPAsUm/pt7wiNQc8olt9fL1M2Y75epkY0a2+3e4RHx/P888/z0svvYRhGJw/f57KlSvj4eHB//73PywWPQResmRJLl68mPm83M4TQgi7S0iAp57SwRdgzBhIH7UTRUeR+riUMc9r71XQycnJBAcHYzab8fT0ZMCAAQwbNgyAoUOH8tBDD/Hbb79xzz33ULx4cQCaNm2Kp6cnzZo14+mnn871PCGEsCurFe68EyIj4e239VexYq5ulSgAQznxU1OrVq3UtQuT9u3bR8OGDZ3WhpuZvJZC3MQSEqB0afDwgD//hIAAaNECrlyBsDCoX1/mgd2QYRjblFKtcnqsSA1BCyHELWnePGjQAL76Sv9+//06+AIkJcF990FoqOvaJwpEArAQQrir06fhkUfg4YchMFAPPV+rTBnw84OYGOe3TxSKBGAhhHBHf/4JjRrBH3/A2LGwcSM0bXr9eYahg7ME4CKnSC3CEkKIW4afH9SpA998A40b3/jcqlUhOto57RJ2k2cP2DAMH8MwNhuGscMwjD2GYbyffvw2wzCWG4YRlf69jOObK4QQNyml4Icf4KOP9O+dO8OGDXkHX5AecBFlyxB0CtBJKdUMCAa6G4bRFhgFrFRK1QVWpv8uhBAiv2Ji9J7ep56CZcsgLU0ft7XM6KhRsHCh49onHCLPIWil9yldSv/VK/1LAQ8AHdOPfw+EAW/avYUOdvbsWTp37gzAqVOnMJlMlC9fHoDNmzfj7e3tyuYJIW5mSsG338KwYTroTpkCL74IJlPez71agwaOaZ9wKJvmgA3DMAHbgDrANKXUJsMwKiqlTgIopU4ahlHBge10mLJly2aWInzvvfcoUaIEw4cPz3aOUgqlFB4esmZNCGFHR47ACy9A+/YwaxbUrl2w68TH68Va994L1arZt43CYWyKKEopi1IqGAgE2hiG0cTWGxiGMdgwjK2GYWyNj48vaDud7uDBgzRp0oTnn3+eFi1aEB0djb+/f+bjv/zyC8899xwAcXFx9O3bl1atWtGmTRs2btzoqmYLIdyd1QpLl+qfa9XSq5tXrix48AU4eRIGDYLNm+3TRuEU+VoFrZRKNAwjDOgOxBmGUTm991sZOJ3Lc2YAM0BnwrrR9V97DdI7o3YTHAyffVaw5+7du5fZs2fz1VdfkZYxJ5ODV155hZEjR9K2bVuOHj1Kr1692L17dwFbLIS4aR08CM8+C//+C2vW6H29GQk1chAaHmtb6t2qVfV3WQldpOQZgA3DKA+Y04OvL9AF+BhYCDwFTEj//ocjG+oKtWvXpnXr1nmet2LFCiIjIzN/T0hIIDk5GV9f+1ZpEkIUURaLnt99+23w9tbzvu3b3/ApoeGxjJ6/K7MGemxiMqPn7wK4Pgj7+0syjiLIlh5wZeD79HlgD2CuUmqRYRgbgLmGYTwLHAf6FbYxBe2pOsrVBRU8PDy4Om/2lStXMn9WSsmCLSFE7nr3hr/+0iudv/pK53HOQ8jSyMzgmyHZbCFkaeT1AdgwZC9wEZTnHLBSaqdSqrlSqqlSqolS6oP042eVUp2VUnXTv59zfHNdx8PDgzJlyhAVFYXVamXBggWZj3Xp0oVp06Zl/h5h73F0IUTRk5am53sBBg6EH3/UW4VsCL4AJ3KofX6j47IXuOiRZb358PHHH9O9e3c6d+5MYGBg5vFp06axbt06mjZtSqNGjZg5c6YLWymEcLmICGjTBqZP178/+ig88YTt+3qBKv45T2Hldpxvv5W9wEWMlCO8ichrKYSLpaTAuHEwfjyULQtffw0PPFCgS107Bwzg62VifN+gQtdAF84j5QiFEMLRtm7VK5o//BAefxz27Clw8AW90Gp83yAC/H0xgAB/3xsH3/374d134cyZAt9TOJcUYxBCCHu4cEF/LVoEPXva5ZJ9mgfY3ts9fBg++EDXBi5Xzi73F44lPWAhhCioNWuytm906qT3+dop+OZbxroUWQldZEgAFkKI/Lp4EV56Ce66C778EpLTVyYXK+a6NmUk45CV0EWGBGAhhMiPZcugSRMdeF99FcLDwR2S7kgyjiJH5oCFEMJWJ0/C/ffrHM5r18Idd7i6RbpmcFgYdOyoe8EnTri6RcJG0gMGTCYTwcHBNGnShH79+pGUlFTga4WFhdGrVy8AFi5cyIQJE3I9NzExkS+//DLf93jvvfeYOHFigdsohMinjCIHlSvDkiW61+sGwTdp6Rosd96N9b9voTp3hqlTYc4cVzdL2EgCMODr60tERAS7d+/G29ubr776KtvjSimsGRlt8qF3796MGjUq18cLGoCFEE5y+rROonH77XroGeCee8DHx7XtAtZ/sZ0LPR/FZDXjgYLUVL0VKh/JPsRVkpP1vm0n5sYomgF4wwa90X3DBrtfukOHDhw8eJCjR4/SsGFDhg4dmlmOcNmyZbRr144WLVrQr18/Ll26BMCSJUto0KABd955J/Pnz8+81nfffcdLL70E6JKFDz74IM2aNaNZs2asX7+eUaNGcejQIYKDgxkxYgQAISEhtG7dmqZNm/Luu+9mXmvcuHHUr1+fLl26ZCv8IIRwAKXgp5+gUSMIDYWxY3XgdQNn4iwsbjySNi+3oZhKIc3DGzMmlKc3lCihqy3doHqbyEVkpF5Yt3at027pfgG4Y8frvzJ6iUlJ0Ly5LuH13//q782bw3ff6cfPnLn+ufmQlpbG33//TVBQEACRkZEMHDiQ8PBwihcvztixY1mxYgXbt2+nVatWTJo0iStXrjBo0CD+/PNP1qxZw6lTp3K89iuvvMLdd9/Njh072L59O40bN2bChAnUrl2biIgIQkJCWLZsGVFRUWzevJmIiAi2bdvGv//+y7Zt2/jll18IDw9n/vz5bNmyJZ8vqhAiX556Cp58EurU0cPNb70FXl4ubZJS8PPP0KiJB6n7DhLR/D/4xhxk/1dhjOFDNoxdqVdhf/st5PI+JK5x4YLO0Q26du2BA9Chg9NuX/QWYZ0/n5Xg3GrVvxdScnIywcHBgO4BP/vss5w4cYLq1avTtm1bADZu3MjevXtpn15CLDU1lXbt2rF//35q1qxJ3bp1AXjyySeZMWPGdfdYtWoVP/zwA6DnnEuXLk1CQkK2c5YtW8ayZcto3rw5AJcuXSIqKoqLFy/y4IMP4ufnB+ihbSGEnWW8r3h46GQWLVrAyy+DyeTadm3YwPlf/2bPL7sYE/cJNdvUpc6yuQQ112/fZdvWpTfDuLS/FvRN3wscE5O1L1jkbMkSGDwYYmP1FEPdulCzplOb4H4BOCws98f8/PSwUOfOer7D21v/3q6dfrxcuRs/PxcZc8DXurocoVKKrl278vPPP2c7JyIiAsNOcy5KKUaPHs2QIUOyHf/ss8/sdg8hRA4OHYLnnoOHHtLDkI8/7uoWAWBdtwF1d0dKWVJpB/zUuSEtl36EyZT11l2hTikqsonVkfug6kP6YHQ0pHcexDXOnYPXX4cffoCGDWHdOh18XcD9hqDz0q4drFyp862uXJkVfB2sbdu2rFu3joMHDwKQlJTEgQMHaNCgAUeOHOHQoUMA1wXoDJ07d2Z6emUUi8XChQsXKFmyJBcvXsw8p1u3bnz77beZc8uxsbGcPn2au+66iwULFpCcnMzFixf5888/HfmfKsStw2KBSZMgKAi2b4eSJV3dokxR/57kaNdBmCypGAAeJtp0Lnldh9zk680Zj4qYTl3V65W9wDlLS9MfTObMgXfe0dMLLvyg4n49YFu0a+e0wJuhfPnyfPfddzz++OOkpKQAMHbsWOrVq8eMGTPo2bMn5cqV484772T37t3XPX/KlCkMHjyYb775BpPJxPTp02nXrh3t27enSZMm3HfffYSEhLBv3z7apf+3lShRgh9//JEWLVrw6KOPEhwcTPXq1engxDkKIW5ae/fCM8/Apk3Qq5cuHegGw7apqfDxx+D73me8bI3CavLCwIrh7Z3rupazfoEUPxutk3GULQuXLzu30e4uPl6PkHp6woQJULs2NGvm6lZJOcKbibyWQuTD6tV6i9GUKfDYY26xfWfHvIN8NOoCcw+2YGDfS3z6xgnKGWezEm3k0vHYXLUvZeIiqZu6R6/WcoP/FreglF6kO2wYTJyoV4g72Y3KERbNHrAQQhTE5s2wcSO88oreVnTkCFy11sNVLp9PY9X9k+myZgyjvJryROhGej9QAqinT8hjxO907Ts4csKHOgpZL5LhyBG9yGrFCr2y2Q1HDoveHLAQQuRXUhK88YYOZJMmZQ3RukHw3TXyey6Wrc79a0ZyoEY3au+YT+8H8hdEjzw0nMesc4iPR29DcpNFZC7z/fc6X/fGjXoba1gY1Kvn6lZdRwKwEOLmtnq1XmQ1aRIMGgQ7drhF4D17Fr5r+xVNQp6mouUEVk9vmv00klINbaz/e5WMQkjR0eie39y5t3YyjnLl4O67Yc8eeOEFvbXMDblFq5w5D32zktdQiBzExek9vR4eOhB/9RWULu3SJikFC745R6NGcHDzWRQGBmC1pLHn54LtcKiXvINoAkleuFwvJLNab61kHKmpemfMuHH69549YfFiqFbNte3Kg8sDsI+PD2fPnpUAUghKKc6ePYuPG+SnFcItZGSLq1hRvxHv2JHvzHiOELsnkaU1htD+uQZU9z/B3vvLkeLpRZrhgdnkydhLFQgNj833dSvULU0gsaQcOJa1kjs62s6td1NbtkCrVjBmDOzfn5XLuQjMhbt8EVZgYCAxMTHEx8e7uilFmo+PD4FusIVCCJc6fVovsPr1V53pqFs3nbjHxaxWWDY0lGYzhtJVxRHecRjF79rB9pRAnig1jrbHd7GxWhDbK9bj+NJI+jTP3zB02aAqWDGwHIuBqm30wZt9L3BSkg66kydDpUrwxx9QxLIEujwAe3l5UdPJ6b+EEDcZpXRO39deg0uX4IMP3KN4woYNxP+ygujvV9P9/GoOlWjGyR//pNUDLTk6ajEA2wMasj0ga/vgicTkfN/GKOZNvKkinqdi9IRw7dpOrerjElFR8PnnOoPZJ5/kObUQGh5LyNJITiQmU8XflxHd6uf7g469uTwACyFEoQ0cqANwu3Ywa5auYpQHR78hp63ZgOrcmTLmVMoA+zo+T4Oln2N466IOVfx9ic0h2Fbx9y3Q/c76VaV4QjSUKQPpGftuOomJuqf71FM6kUZUFFSvnufTQsNjGT1/F8lmCwCxicmMnr8LwKVB2OVzwEIIUSBWa1YBhZ49dW9ozRqbg+/o+buITUxGkfWGXJD515zsCj1EdJf/YJhT8MSCyQQN762WGXwBRnSrj69X9rySvl4mRnSrX7B71n6Qf63ut9fVbv74Q/9tn3026wOGDcEXIGRpZGbwzZBsthCy1LWlXSUACyGKnv374a67YNo0/ftjj+WrcpGj3pCTLqTx590Tqf1gEBVSozE8PcFkyjGNZJ/mAYzvG0SAvy8GEODvy/i+QQXuke3oMZrRl9/GYgHefPPm2QscFwePPAJ9+kD58npvb506+bpEbsP6BRnutycZghZCFB1ms57v++ADvZe3bNkCXcYRb8ibZuzA95XnuD9lKzuq96bm319iSjx+wzSSfZoH2G0ItGpVwJLGyVgPAuPi9GhAUWc262IJJ07A2LEwcmSB6jLbe7jfXiQACyGKhvBwePpp2LlT94g+/1xvMyoAe7whh4bH8teM+dTeuZ1/zz5I58gwBpqOs+e9uTQb83D6NpgAhxaOuXoeu9emA6QwnF1rdhJYtSqcPKmTcXgWwbf52FioUkUH2ylTdBarBg0KfLkR3epnmwOGwg3324sMQQshiobz53Ut1z/+0NuMChh8ofDzr6HhsfwyZS5TZ7zBG+t/4LfIx/inUQvWLF5J43f7OWUP6rXz2NF+XpiwsmXVvqKbjMNi0duK6tWD2bP1sd69CxV8wf7D/fZSBD8aCSFuGatW6Tq9w4frYdyDB6FYsUJfNuONt6CroCf/sIuQX/5HMasZAGWk0L7cciaGV+HBboVunk2unceOq1gKgINbIuGh5vpgdLRblFi0ye7dekvRpk3Qowd06WLXy9tzuN9eJAALIdxPQgKMGAHffKN7Py+9BD4+dgm+GQryhmy1wvKXF/LTl0OpzAnMhgkDRZrJk43Vgpy6qOfae52+rTRWDPzPnNb7gDt2dNscyNeZOlUXyyhdGn76SS8gKwKZrApLArAQwr3Mm6cDbny8XnTz3ns6+LrYgQMw+6FFjN/9ALu9GzG4x1v4lLiclcUqoCEBTlzUc+08tsXTk1NGRaomx0P9+jr3tbvLqF1cp46e1588Wa90doGM+fTYs1cIKOvjlEQdEoCFEO4jJgb699f7PRcvhhYtXN0izKmKWW8f5fXPa+JX7D56DpzBmaFdOLJIDwFnZLFy9qKejIVFXhfPc1vSeY7eFsD0MoO44FefJ53WigK6eBH++1+9iv2993TBjPvuc1lzQsNjGfnjfgIXXGLg+U3suu82Rl9OBRybqEMCsBDCtaxWXTT93nv1fGVYmE6uX4DtJvZ26OPf8RozisdT49nU6zDjZ5SlcuVB+kFvb5emNuwTXIXAv+bTauwrANw1dhkb73qT3RuLMwV0cpLy5eG775zWJpv89Rc8/7z+sPX6665uDUrByI8TKR/qyaKUx/EmBfOvnjzx2DhClno79m+qlHLaV8uWLZUQQmSKilKqY0elQKnVq13dmkyXL6SpNc1fVlZQVlAWk5dSa9e6ullZoqKU6tJFv24lSyr18cdKmc3q/XfMqirHVUqKUqpzZ6XatnV1S7OcPq1U//66zQ0bKrV+vatbpI4fV6pnT92kaX7PKauOx8pseKiP7xqoary5qND3ALaqXGJiEZmhF0LcVNLSdEKNoCC9v3fmTF1APV1oeCztJ6yi5qjFtJ+wym4pIm3xz6KLRJVvx53hUwEwAA+s8O+/TmvDDZnNutDE5s06E1hCgp4r9/Sk+/7POE41Tuy/oDNzuFNJwthYCA2Fd9/Vf3MH7o/Oi9WqX7pGjWDtqlRq9DzA6t61sRgemaUhN1YLcniiDhmCFkI4X8+esGwZPPggfPGFTrqQzlWJ8xPOKUaMNPjmm5L8VKolJR/vRa1fJ+hi7zmkknS6rVuheXM9NP/DD3qhVZUqkJysh/AbNcKnjt5ydCYihhqBga5PxnHsmA66r74KwcFw/HiBs5fZy759MGgQbFpn5qs6nzLQPIulI/9kxPJgHnni48xFdftqNGG8g+f0pQcshHCOpCSdaAFg8GD4/XeYPz9b8AXXJM4P+3ANsRWbs352JCNHwoOnplNr9hhYuRI+/FB/d1WPLSFBv16tW2clp7jnnqzX7coVnaxi4UJKN9YB+MKeaN0Dtlp1EHY2i0VnsGrcGN5+W6eSBJcG39RU/acMDgbvXduIq9aGZw+OxqtVML0alWd83yDiGrdgertHiGvcwimJOqQHLIRwvJUrdRB59VV45RV46KFcT3Vm4vwzv6zg3NC36ZiwiVjvGiyYeY76T191Qrt2rgu8SsEvv+gax2fP6mQkORVY8PcHPz+IjqZ8974AXDkYA32bwYABTm402RNq3HcfTJ9+3YcsZ9u0STfJf/cawsu/RcMzazGKV9IfAB98EIA+5ZxfmlACsBDCcc6d04Fj9my917Np0zyf4ozE+UrBpn4TaTNvJGVRWDw8qbhwFgHdXDcveZ0XX9TBq3VrWLpUd91yYhi6txsTg1+dKlgxUMej4fZn4fbbndvmpCTdOwe3SKhx+bLugE+ZAr3KbmCBdzdM8cm6atZ33+mV9y4kQ9BCCMf4+2+9yuWHH3R5vJ07bZpHtXed3GsdPAidOoF5XigGCgMwGQrP7Zvtcv1CMZv1nC5kFZzYsCH34JshMFBv6/HyIqTKZMKKpefDVEpf09HCw/W9/Px0nu59+/R+bhcG32XLoEkT+O6zBMIaDWVu/wWYLKlZJ2zb5rK2ZZAALIRwDD8/HRi2bIEJE8DXth6soxLnp5kVix/5nucbryE8HBKHjdVtMpkcssgq3yu5N26Eli3hnXf07x072l7j+KoVz2uav8qKy+k9+cqV9QppRzl/Hl54QSdM+fFHfaxTJyhXznH3zMPZs7poVrdu0CtlHnG3NeKu/TPwKVVM/50d9PcuCBmCFkLYh9UKX32lC6i//77eVrRlS4F6QfZOnL9n0REu9B9Cz4vL8as6kPobO1ClSkd4eOUN6/UWVL5Wcp8/D6NH69cuIADuuiv/N/zvf3UuZaBRudMkrY0BWuj54ZiYwvyn5C40VA+TnzqlE2qkz6W6Qmh4LJ8sieTgRn8SVzWhYnI8u+u/TOPIBXrl+Kz0rGo9ejjk711QEoCFEIWXsbdj3Trd9bBYdE/DlQn1N2zAvHQVm/48TfPts1CGB+HPTeOer5/PGvtz0CKrG63kzhaAw8L0UG1cnF6c9uGHULJk/m9Yt27mj/0OjOPt87NJSrqAX2CgY/YCv/qqHh5v2lQH4tat7X8PG4WGxzL8uwMEzkvm2dhNrLvtDM9Un0mDg3/Bxx/DsGFZ27BcuaguBxKAhRAFl5qqh5fHjYMSJfTCloEDXV/JZsMGLPd0xiMlhTuxElXhDsqv+IXmQVWdcnubV3JXqgTVq8PChTr9ZkGdPq2LWPTogUf1qpTacJGD+y9Qp2pVWL684Ne9mlJ6T7GXl17dXKmSXmDnhJShGYUSrk37abXCGx9cpOIiT/5OexgvUjFf8GRQl7f5u/v9/DpyoMPbVhgSgIUQBXfsGHz0kd5W9NlnUKGCq1tEYlwKq5/8gftTUvHEigUPUvt1wD+fwTe3N31b5LaSO7CUt36dduzQK8MbNID16wv/gSUuDoYOhblz8a2r9wLHh8dQx17JOA4c0NvIOnTQvfTu3fWXjQrzWuY2nB971JNfJ1fk+JraTPcZik/aFQwASxpNTx1kek3XF/LIiyzCEkLkz4ULMGuW/rluXT38PGeOWwTffz5ax+mAYO45PAez4Uma4UGqpydjL1fMVzrLjDf92MRkFFlv+rZeI6eV3C3PHGHhj8P1fOnp0zqBBthntCBQB11iYijdRH/QuLA3Wi+IGjVKj1QURGqqHt1o2hQiIqBWrXxforCv5bXD+cpicOrfGrz8SDk8IrYT4d+Me6+sxurkNJL2ID1gIYTtFi7UPa0TJ/Qe06AgqFnT1a0iLuoC4d1H0/3wlxwzVWNol3e4UomsWr0V63H82vnXG7B5DjcXGeeELI0k4XQC72z5hcfWz8coX15v0+nXz77D9P7+ULw4REdTrqdOxpFyMAbueTZrX25+RUTo6YRdu+Dhh/Wcb+XK+b5MYV/Lq4ftU06W5uzfTTHHl6J71d/468TjpJQpyysPv02MT2mnppG0BwnAQoi8nTqlFwn99pveXDlvng6+LqYU/Px5PHe/3oJ7VSxb7niVp1p1IMnXByCzVi/kL5OWPbJxZa7kjo+HRs/rRWoTJuhgaW+GkbkXuFjNKrxY6n+U8m1Pb6UgMVEviCtVKv/XvHgR/vhDp7osIFteyxsNUVfx9yX6dCr1llyi9b4dbPdJJLxvJS63KI2R+hY+r71Gp6NJhCyNZHpAQ6r4+zLeyaUhC0oCsBDixiwWvTXm+HEYOxZGjND7KF0s5ud/WfL2Or453BFV9Qk6TOxD60faUmbCKpIKmUmr0Nm4Tp3S5Xbee0/X5I2MhNtus/n+BZKxF9jLi831nqRsIpCYoO87cWLmNqUb+usvXfVpwgRo1gyiogpdyCGv1zKvLVv3+jdlxbhdLLzUFx+SMa7Acz5jub/H09BcZ7LqU6ZMkQi415I5YCFEzg4fztpONG2azmT11lsuD75pZsW6LmOo0v9unj78Nmu8OvP4zw9Q7ZG2gH0yaRX4GlYrzJgBDRtCSAhs366POzr4gl6B/vffAHQotYNK+8Oyhqbz2gscF6fTRvbsCX/+qXu+YJcqSnm9lrkNUX+04KHJrt0AACAASURBVBDPPAPvvVCO19RMfEnGA1AYDC+ZUCQD7rUkAAshsktN1SubGzWCL7/Ux7p2hXr1XNsuYO9fR9lc7j7ar/wQA/DEiqc1FY9/wzLPsUcmrQJdY98+nXxkyBCdOnLnTufujw0IyBzefurYB7wZ/WK2oekcKaVXYzdsqAsTvPee/tBQkL3Iucjrtbx2iFopuLy/EtsmteG375PYVf8hHrr8M4ZhgMmEh68PDZ/KvZhHUSJD0EKILBllY3bv1gtvHn7Y1S0CdHrk5Q98Qaflo6hmGET2Gka9ldNzrdVrj0xa+bqGUjp384kT8O23Oheis/dC790L//sfvPEGaZWrUvPQcs6fh9JXpam8Tny8XpUdFJTVc3eAG72WVw9Rp10sxrnlTUiOqkTxKhdYs6YkTT6wwtPjoX17WLvWbbJY2YMEYCGEFhKiiyYEBBR64Y09/fOPXr80JOoINQLvotrir9hnMfFpYA3q7N3KwUat6OFTjT6uaNzatTrVYfHiuvpPpUqu24517Jieu+3dG1P1QEqtvcjefRcoHRiYPRmH2ay3jQ0YoNu6YQPUrw8erhkQHdGtPqPm7SJ+awAJqxtQ23KYL8o9jvXL8QQ3b6t75hkfZjp0cEkbHUUCsBC3uox53ttv17l9x43L/4rZfLA1KcOlxf8QPfgDfj3RF0utF2m+ZAJN7/UkNOKEXrRTuhbL2ul9qWtyy7PsKAkJusjBrFnwwQe6gIINpRYd6qq9wL719F7gsxHR8MQTWT3GTZv0p5ldu/SWonvvdViv11aNigdg+rscDbdv5h2f7nS0/IvHFV+8vBP0Ca7OquZAEoCFuFWdPKlz+gYEwOTJeqWzDYUAHJHVCLIHz/ChM2g2/Xkaopjq8S/mWS3wuUcHkcLuKy0UpWDuXP26nTmjV4TbsrrYGTICcHQ0/k10HeCL+2Lg+W56UdWrr8LUqVClih7hcHEtXLNZL85+/3142uMXphv/wbiidE/8px91usubXJ5jDoZhVDUMY7VhGPsMw9hjGMar6cffMwwj1jCMiPSvHo5vrhCi0KxWmDlT93wWLoSKFW1+qr2zGkFW8ASIO3iRv+u8TLPpQzBQgK7V67MxLPN8e+zRLbAxY+Cxx7LKLH7yiS676A6uWvFc9p6mdDJWE+59O6Sk6Kr0n3+uk6js3evy6YWtW/X6tP/+F3r1gs8bfImh9N8bw4A9e1zaPmexZdA/DXhDKdUQaAu8aBhGo/THJiulgtO//nJYK4UQ9nHwoM6MNHiwnrvcuVOnKrzKjerY5hVA85JbkIxNSGb2bJjQ7Ge6HZrG/iYP51qrN7e9uA5LPWix6JKBAE8+CZMm6dq9zZs75n4FlbHiOS4OzzIliQroSFS8v84DbbHonNNffOHQ6YW8JCXpQYPbb4da0f+wMmQ7v/8O3p+Mc2htZneV5xC0UuokcDL954uGYewDiv4GLCFuRVarTq4waxY888x182t5DREXtveZU1KGkidTCFxRjGdOwN13PsvxYa1p9GBzvTgoh9qtI7rVz9ZGyP8+X5tFROhV4TVqwO+/68VK9d04xeH27TqQAY8V/xOPCG8o3k0HXhdbuVJ/7jtz+DyrG4zkrv0zYE1vGP4HdOmiT3CjWr3OYKiMbr8tJxtGDeBfoAkwDHgauABsRfeSE3J4zmBgMEC1atVaHjt2rLBtFkLkx9q1umbrxIn695QUKFYsx1PbT1iVY9aiAH9f1o3qlOfjeckI8A2P7ub2Y7tQx0vw7LG5WAxPFk45ynMvFrNpMW5h5qFtkpSkJyc//RTKltXDt488UqQWBB0s24aTV8rQ4fJSl7YjIUFXLdz77QbeLPElPTyW4H3pnK7T+/777jOEf/myTmYydKhd/86GYWxTSuVYa9LmAGwYRgngH2CcUmq+YRgVgTOAAj4EKiulnrnRNVq1aqW2bt2ar8YLIQooMVEPL3/9ta45u2WLTot4AzVHLSandwQDODKh53U9ZNC9z/wkuvjnu4Xc/lw/vC2peACH/BpS/I+5VOrSxPb/NkeKiNDlFQ8f1r3fTz6BMmVc3SrbLF8OP/wAs2ezq+EjeB7aTwPLXpd8blBKpwx/6SWoE7+BMDriaU3VwS1jBMadREToiemwML3n2E5uFIBt2vhlGIYXMA/4SSk1H0ApFaeUsiilrMBMoI29GiyEKASl9HBpw4Z6sdWwYXpRSx7BF/KeXy1slqkrV8D4ei3F0oOvMjyo9d8n3Cf4gl4VXqkSrF6tX7+iEnxBf2j48Uc4eZK0KlUJUDGcPev8Zpw4AX37Qr9+ihbljvPrC2F4Gukf2jw8dOpLd3D6dFZpzeBgOHTIrsE3L7asgjaAb4B9SqlJVx2/ui7Vg8Bu+zdPCJFvFy7ACy/ofZ6bN+th1OLFbXqqLTmQ+zQPYN2oThyZ0JN1ozrZHHzXLU6kWTN4c+ODWAwvlMmE4VMMo1PeQ9cOpZQOWj176sVK5cvDunVFcyFQVb3/l5gYPKsHUoqLxO674LTbX50KO/KvQxyp05XF59oS8EArvbjKXRZZKQXff68b+uKLutAIQLVqTm2GLT3g9sAAoNM1W44+MQxjl2EYO4F7gNcd2VAhxA1YLDqIWCxQurROH7V5M7Rsma/L2COP8rUunEllUev3CepVjcBL+/lgWTs81/2D8eGHeuGNKxfcHD4M3bvrrFDnzuGS7qI9XZWMw6++DsbnduSShtLOoqKgc2cYOiSNCeUmstsjiBpxmzHGjNEPrFwJ7vA3P3hQ5zZ/+mkdgMPDnR54MymlnPbVsmVLJYSws/BwpVq3VgqU+v13V7cmm38/2aD2eTZWCtTWBv3VpSOnXd0kzWxWauJEpXx9lSpZUqkvvlAqLc3VrSq8c+f0v4NPP1WnIhNVTQ6paVPMDr3lb5tiVJ/b56rRHmNVN69F6mjlYN2G3r2Vio526L3zLTlZqQoVlCpVSqnp05WyWBx+S2CryiUmSiYsIYqqa1fq/vyznnhzA+cWrSfuP2/S/sxa4ryqEvnpIloO6+nqZmVJS9NjpV276lKLGT3Hos7fXydWSUmhfJ3SxHiV5vgJx93u0zmnWfjyXv4+9xTepGL28GSzXxBxE6bTZuQQ91k1vnev7u36+OiVzk2b6rl+F5NyhEIUVf366RW6Tz+tS+E99pjL3/CUgkmDluJzfxfqn1mH1TBxYNIn1HeH4JuUpIdAL1/Wb8Tr1+vtWTdL8AX99z91CkaPxsMDxpSczG2b/rb7bZKSdN2O4QPK8eSFX/AlGU8seFnMbApoxOuqnsv/LQJw6ZJehBgUpItlgE5x6QbBFyQAC1G0nD6t31RAFwAIC9OrOJ1R8D0Px7adYVGFJyg563e8ScUDBSi2zV9qc6pKh1m+HJo00akkFy/Wx8qWdY8g4UDPXwqh8b7f7XrN1at1B/KrT84z0/8pBqV9iwLSDAOzyZON1YKckxY0L0uW6L/55Mk6A8j997u6RdeRACxEUaAUfPMNNGiggwhA27a6ALyLWdIUfz05B79WDel+Zi4p1VMxe5pIMzwwmzxZG9DY5lSVdnfmDAwcqAsPeHnpDyyPPOKatjjLjBmZdZzPlwyk5PkYu1w2MVEXUurUCTpdWkhcucb8J2EOM1v34YlHxzGpwwCeeGwc2wMaOi4tqK1GjNA9XV9fWLMGpk/XixPdjMwBC+Hu9u+HIUPg3391PdRBg1zdokz7lh7n7KMv0OP8X0T6t6H/vc8RVbMKq2KDaHt8FxurBbE9oCGGq3pEQ4boghNvvaULEvj4uKYdzhQdDQsWQFoayeWqUv7sPqzWwpX7XbBA79aJi9OxbVzC33htLkvY5K+ZFGki2WxhQ41mgAPTguZFKb0LwNNTb3Py89PVHnLJ+uYOJAAL4c5++klnDPLz00khnnnGZYXTr5YStoF/PwhjV9gZhqgwtg+YTPNvXyZp4j+QmMz2gIZsD8iqM+vUHtHRozrQVqoEH38M772n5wBvFVWr6g25J09irRxItcjlxMXpbeH5dfKkzmR1Yv56vi89ibr/vY8aHz4LlyeCtzcdvbwY7+i0oLY4ckR/2GrfHt59V+/p7ukG6w7yIAFYCHeUlqY/ybdpA48+CiEh+Sob6EiR7/1MzQ/+wz0qjbtM3qTMnkOLAX0AJxdKuJbFouvdvvUW9OmjP7zUqeP4+7qbq/YCe9asSqmwixzYd4HKlW2vgqQUfPutzuHc/dLvrDMexeO8FcaHQo9G2fbx9mke4PyAmyEtDaZM0dMyJlPm0HtR4fqP0kKILOfO6SHmjHnKunV1bl83CL4XzqSyqM0H1H5/AF4qBU8sFCOVUjH7Ms9xRCIPm+zcqYPC66/r4cfx4x17P3eWkQ0rOhrL4BfwIZmj52wPvgcP6rwZg5+zMPa2SfxEfzyUNeuEsDD7tregdu/W6yCGD9cN3rtXL7YqQqQHLIQ7UErv4339dZ2Nadgw3aMzmfJ+rhOsm7SJsm8+R6+03URW7UK9+LVgNueYVtDpPaJ58/QWrDJlYM4ct9iO5VKBgXrPq4cHAfWKk4KeFs5LWppeMDxmjP6zrh74HXf98IYe1t22Lde/t8ukpuotV3Pn6p5vUfyb55ahwxFfkglLiBzExCh17706e1CbNkpFRLi6RZni4pQaft9uZcFQJz0D1b6QP/UD69cr9dFH+rurXLmiv8fHK/Xii0qdOeO6trgp64WLaprnK+rLvstveF54uFItWihVjGT1SqddKjZWKZWaqtTChUpZre7x91ZKqZUrlRozJuv3lBTXtcVG3CATlgRgIVztzBmlatVSaupUt0mHaF23Xu3qMUJ1K7VeeXkp9UfvWSol/ryrm6UlJCg1eLBSt9/uNq+X20pNVRYMNbfhmBwfTkpSatQopdp7rFff+wxWl8pVU9bKlZW6fNnJDc3DuXNKPfOMDll16ih14YKrW2SzGwVgGYIWwhU2boQvv9QrXcqWhchIvejKDcR9v4RyT/eiMRb+MKZy4sdV1Oz/rKubpV29H2bYMD1u6ibD9G5l+HA4dgx++40E70oUO3P9XuB//tHLDapFreAf4z5MV9IgxdCpTf38XNDoHGSU1nz5Zb2n+8039SpnXxfvM7YTWYQlhDOdPw9Dh8Idd+iUQkeP6uNuEHwtaYq/n/qF4k8/jAcWDMDbw0zNY2GubppenPbQQzrXdYUKutJTSIhb7/F0qbNnYcMGAM6XCqT0haxJ4MREvWOnY0com3KCv4s/hEml6Qc9PHTRZndx9iw8+6xOHbl1K0yYcNMEX5AALIRzKAW//aYzWX39Nbzyil616SbbZPbssrKmQl/u++FxzhavhipWDEwmDHdZdOPnp0sHjh8PW7bku8ziLScwUG/iTUsjvkQlyqfEUGPEXzQYsJNadS18N9PM8OGwcm9lvHp21x9k3KVWr9WqRzqUgnLldAKaTZsgONi17XIACcBCOENaml5eWrmyfjP57DMoWdLVrSLliuLdd6F5Sw+2XGnKticnUS1xFx6rV7u+dmtUFDz5pM597eOje0CjRumUkuLG0pNxLF0ezk5K400qcb+3JvLHIB5MmUX8bTUJGXoEv+IG/PqrHo1x9d8b9FRMx456pCMjZ3dwsFuMEDlEbpPDjviSRVjilpKaqhdWZSwYOXZM16F1E9vn7FNbfe9U97BSPfmkXkzsFsxmpSZMUMrHR6nSpZVat87VLSp6Fi9WCtSg5z9XFfptUG1Zrz423lCbSwYpBSqiehOloqJc3cosqalKjR2rlLe3UmXKKDV7tl59fRNAFmEJ4WSbNumJth07dO/tueegWjVXtwqAy0v+5figD2kcE0aSR0kmvn2eFh+6ulXptm/Xr1V4ODz4IHzxBVSp4upWFT116kDXrsQlpXGH1zp+NUbjpcxwEWa07sOEe57hsJtMfwD6b714sU5A8/nnbpF4xhlkCFoIezp/XifPbddOr9qcP18vInETO1/7Bt/7OtIwZgWehsL7l//R4sMHXd2sLCNH6rnLefP0ayfBt2Dq1YNlyzjTsBkdj27FS5kxAIvhQaJvKSqXKe7qFuq6zCkp+udXX9W1mX/99ZYJviABWAj76tkTpk3T2yb27tWf7N0gQ098PDzxBByd8gcGCtALXv0O7nRxy9DzjydO6J9nz9avW9++rm3TTWJEt/psrNeGFE/vzPKQ4bWCXVOt6GrLl+sCGRkpQ7t2hQcecG2bXECGoIWwp8WLITYWGjVydUsAvZB01ahl/PDlJX5L6Uun/7wJP68Ac6rrV7wmJur9qt98o7dmTZuWlcdYFF63bvQpVw6Gf8IrPp7U2buVg41a8ejgvq4rnnDunN6//f33upfeqZNr2uEmDD1H7BytWrVSW7duddr9hLiVRUecZX/PN+h64nt2lWiHx4Z1NG5i6P2hYWE6+Lpqxev8+Xqo/vRpeOMNXTLwJtrf6Ra6doWLF3XSF3ewbBkMGKD39r75Jrzzzi1Rn9kwjG1KqVY5PSY9YCHsZfFi/WY3ZoxLt8pY1qzn4GtTKb99KR25yJZ736LF/LcxFU8fCm/XzrVbTaZO1fugg4Nh0SJo0cJ1bbmZBQbqoV53UaEC1KqlA3GzZq5ujVuQACyEvcybpwPKBx+4rAmHf9pA4IBO1FcpWDE4++l3tB42MPPxUFcVT1dKDz+WLQuPP64r67z8suzpdaSqVTOTcdxoH63D/k1YrTBzpp7TnzJFf+Bav94t1kS4C1mEJYS9hIfr3pwL3mBSkq18PXgb3w4MwyM9raBh8qB8SmzmOaHhsYyev4vYxGQUEJuYzOj5uwgNj83lqnZy+LAeDu3eXQeDcuX0PKAEX8cKDNRB8OTJXE9x2L+JAwf0/O7zz+u6vRmrnSX4ZiMBWAh7SEnRbzTNmzv91uG/RLKrXEf+M7Md5e6oi8nHO8c0kiFLI0k2W7I9N9lsIWRppGMaZrHApEnQpIlOH/ncc3rptXCO5s3hmWduGPTs/m/CbNb5mps2hYgImDULVqyQnN25kCFoIexhzx5IS+OdY178OGqxU4Z3LyWYCesZQpcNH5Bi+LL/ta95bdJDsDEgx0VWJxKTc7xObscLJSZGbyXasgXuv19XfgoMtP99RO5at9ZfN2D3fxPx8fDRR3o73hdf6NSrIlfycVQIO9j4704ue/vyb8mqDh3eDQ2PZfALUxnf5HliyzWh14a32FunN6YD+2g6+T+ERpyg/T/J1DzflPb/JGe7fxX/nFcZ53a8UMqV06uaf/kF/vhDgq+rKHXD6kZ2+Tdx5QrMmKHvVaUK7Nql10NI8M2TBGAh7OCNK9Vo8tqvHPPPetOx9/BuaHgsv06Yy2dfj2TEnlnUtB7moxZDOT53MiXqVMpzPm9Et/r4emWvnevrZbJfUob16+G++7KKJ4SFwaOPyryfK5UvD6NH5/pwof9NrFmjVzQPGaJ/BqhevaCtveVIABbCDk4kJqMMj+uCjb2Gd5WC0Nf/5ou5H1JMpeKJBQ/DimeJS5lBPq/5vD7NAxjfN4gAf18MIMDfl/F9gwo/TH7pkk4leOedeij+yBF9XAKv65UvD9HRuT5sy7+J0PBY2k9YRc1Ri2k/YZX+QHfhArz4Itx1l573Xb5c/yzyReaAhSgsi4XQX0cxK+g+/mx0d7aH7DG8G7vrHHt6DOe7mNlEewRQgksoBWaTJxurBWUGeVvm8/o0D7DvvPSyZTB4MBw/rt+QP/rILcosinSBgXo+/gZu9G8iY1Ql44NdxqhKh0VjKLtzG7z+ui5jWNwNcksXQRKAhSisqCiaHd1NiaB7sx0u7PCu1QpLng+l5azn6aTOMK32EKbcfy9Nzhym7fFdbKwWxPaAhgSkB/kq/r7E5hCEHTLHC7pbPmGCHm5eswbat3fMfUTBVa0KS5cW+OlXj6qUSTrPZW8/koF3Wz/GF19Nhdtvt1NDb00yBC1EYYWHA9D1ie52G97dtw86dIDFM2O5UCKA04u3EvDbO5hK+LE9oCFftnuE7QENswV5h8/xZpg3T+e7NgyYM0dvN5Hg654CA/U+YLO5QE8/kZgMStF77z+smPUCL26YC8Di2+pL8LUD6QELUVjbt0OxYnTq25FOhUwuYQ5bx8GhE5kbGcx+/3cZMvsF6jwxBMPLkz7p5+SWtSjju8MyXZ08qfM3z5+v8zdPnAiVKtnn2sIxOnfWe6/N5gIlPmlmXOKl3z+ly6EtRFSux+IG+oOWw0ZVbjFSjEGIwurSRVf2KeS/7ahxv1Lr7f6YsGIxTJz/cw239XRhzuYMSsF33+nsVVeuwPvv659vkN5Q3AQWLMA8YCBpKWYm3jWA2S3vx+phwtfLZJ/Fe7eIGxVjkCFoIQqrdm3o0aPAT7+caGbRnROo8fYTeGAFwOQBt+0Ms1MDC2niRJ1RKSgIduyAkSMl+BYVSumKU4mJ+X9urVp43dmeNfNXsaTr4ygPk/1WzgtAesBCuNSyZfD9wJX8FNeFyCodqXduI4bZrGv1rlzpuqpFFosuG1ehgv4+fz48+6ykkixqEhOhTBn49FM9anEjaWnw2Wdw6BBMn+6c9t0CpAcshKNYLHmfk4Oz0Ul80nU53brBNv/ORHy9ifqxqzFWrdLbOlwZfDNWgPXood+Uy5aFQYMk+BZFpUvrLUI32AsM6OxV7drBiBFw4kSBF22J/JH/o4QojA8+gGrVbH7DUus3cOie50ipXo9XV/Ri/CsniYiA4MFtdMKDXNJIOoXZDOPG6bJxkZHw2mtgMuX9POG+DENvRcptL3BKCrz7rq7idewY/PorhIZKpSonkYkcIQojPBxKlLDpDSv+x6XcNrAntZUFKwYnRk9l1Ec6dWVuCQ8A58y3HT8OvXvrOd5HHoGpU/Xwsyj6bpSM48wZPez82GMwebLO4S2cRnrAQhRGeHieJQitVpj52WW8BzyCh9IB1jB5EFjyQuY5Ti8VeK2KFfWb74IFuhckwffmUbVq9iHopCSYNk0v0AoI0FMO//ufBF8XkAAsREGdOaN7FjcIwAe2XuDuu2Hw68VZUn0wyrtYjrV6nVoqMMPatXDvvXDxoq7XumIF9OmT9/NE0TJggF5XALB6tV7N/tJL+u8PuoKRcAkJwEIUVHoGrJwCsDnFyuI+MynfujplI1YyezY8ciQEj7DVOS6ycmqpwEuX4JVXdPL8qCg99yduXvfco2szDxkCnTrpeeHVq/VCO+FSMgcsREFVqqQD2dUBeMMGTn7xOwnzV9PzSjh7yndk1sLqlGub/ni7djmubh7RrX62OWBwUBrJ5ct18YRjx+Dll/WiqxIl7HsP4V6uXNElCWfMgOHDdSIVPz9Xt0ogAViIggsKgilTsn7fsAFLh7upZDFTCTjw0Cga//aRTWX5HJ5GErKKJxQrJsUTbiUWi14kuHkztG7t6taIq0gAFqKg9u2DOnUyV0BbV4dhWNIwAGUyUa9lqWzBNzQ89oYB1u6lAjMsXAgtW+oFN3Pm6L2hPj72v49wT8WLZ/+gKNyGzAELURCXLkHjxrpHmS62Tkeu4IPV4/pFVhnbjGITk1FkbTNy6F7f+Hh4/HF44AGdThL0amcJvkK4BQnAQhTEjh16SPeq+d9tlxvQg8WcfOH6RVZO3WakFPz8MzRqpFNIfvghfPKJ/e8jhCgUGYIWoiC2b9ffrwrAZWaFsIRJ8NEFKOWd7XSnbjP6/HOdxer22+Hbb3UgFkK4HQnAQhREeLhOVnHVHsriUeEc86lP/WuCL+jtRLE5BFu7bTNSCs6d03mbBwzQeZuHDpVUkkK4MRmCFqIgMjJgXbXIqtrZcOKqtMjx9BHd6uPrlT0Y2m2b0bFj0K2bTqqRlga33aa3GEnwFcKtSQ9YiIIICclWE/fs7pNUsMaxOyjnrFgO2WZkteqycaNG6d8/+UQqFglRhEgAFqIgunTJ9mv0H9spC5S6O/e0lHbdZnTqlC6asGaN7vnOmAHVq9vn2kIIp5CPy0Lk144dsHRptlrAW5Ia8xqTqdEn2DltKFNG3//bb2HJEgm+QhRBEoCFyK8ZM6Bfv2zzv2uia/B7wGuUq1nScffdswcefjireMLatfCf/9iUaUsI4X4kAAuRX+Hhumj9VfOtxdcsoWODU465n9mscza3aAFhYToDF0jgFaKIkwAsRH5YLHoI+qr9vymnEph+9D4GWL6z//0iIqBNG3j7bV0qcO9e/bsQosjLMwAbhlHVMIzVhmHsMwxjj2EYr6Yfv80wjOWGYUSlfy/j+OYK4WJRUbqg+VUB+PjCCAB878h9AVaBjRoFJ0/CvHnw669677EQ4qZgSw84DXhDKdUQaAu8aBhGI2AUsFIpVRdYmf67EDe3HGoAJ6zSxwJ62SkAb94MMTH651mzdK+3b1/7XFsI4TbyDMBKqZNKqe3pP18E9gEBwAPA9+mnfQ/0cVQjhXAbDz2kh4WvSu9o7Agn1gigRptC9k6Tk+HNN3UO6Xfe0ccCA3ViDSHETSdf+4ANw6gBNAc2ARWVUidBB2nDMHJ89zEMYzAwGKBatWqFaasQruftDc2aZTtU7vh2jvg3J6AwiafWr4dnnoHISHjuuazqRUKIm5bNi7AMwygBzANeU0pdsPV5SqkZSqlWSqlW5cuXL0gbhXAPSsHIkbBuXbZDD3mEsrrr+IJf9+ef4c474coVWL4cZs7UNXuFEDc1mwKwYRhe6OD7k1JqfvrhOMMwKqc/Xhk47ZgmCuEGrFb44gudgnLHjszDx49D+KW6VOjUJP/XvHJFf+/WDUaMgF27rsuwJYS4edmyCtoAvgH2KaUmXfXQQuCp9J+fAv6wf/OEcAOHDkHnzvDKK9C1KzzxROZD0T/+w/NMp3kTs+3Xu3gRXnxR93rNZj3H+/HHUNKBSTyEEG7Hlh5we2AA0MkwjIj0rx7ABKCrYRhRQNf034W4uURGQtOmuv7vzJk6BeVVw8O+839iHG/RuJmNyylWQo27CQAAIABJREFUrICgIF1EoUOHbOkshRC3ljzfNZRSa4HcUu50tm9zhHATSUng5wf16um9uE8/DVWrXnda6cPhRPo1p12JPLJSXboEw4bpIF6vnk4jeccdjmm7EKJIkExYQlzNYoFPP4UaNeDIEZ3u8Z13cgy+mM1UPb+LM4E27P/18oJNm/QirogICb5CCAnAQmTav1/Pyw4frvfi+vjc8PSLW/ZTTKVgDW6R8wmJifpaFy7o4gmbN+u5Xl9fBzReCFHUSAAWQildzD44GA4cgJ9+gtBQqFz5hk+LWRkJQJlOOfSAFy2Cxo3hs890AQXQQVgIIdJJABbCMODwYejVS6d97N/fpkpDK/wfpjSJ1OlZP+vguXMwYADcfz+ULauHnXv3dmDjhRBFVb4yYQlx0zCbYcIEvQe3TRuYOlXP0+ZDRAQUK1+aygFXHezVSwfdMWPgrbd05iwhhMiBBGBx6wkP14Xsd+yAlBQdgPMZfLFaeWzBY1SpPhDD6JV1fPlyXb2oTh37tlkIcdORIWhx60hJ0Sua27SBuDg9zzt2bIEuZT5whK4Jv9G0wqnsDxQvLsFXCGETCcDi1vHttzrg9u8Pe/bAAw8U+FInF28HoESHqxZgLVyo9wyb85EVSwhxy5IALG5uycmwc6f+edAgWLUKvv++0CX+LvwTjhlPqve8Kgf0ggUwezZ4ysyOECJvEoDFzWvtWl06sHt3HYg9PeGee+xyac9d4ewzGlEv6KqtReHh0KKFTSuohRBCArC4+Vy+rAsn3HWXHg7+4Qe7J784c9mX/eU6ZHV2U1L0sHZzG7JiCSEEsgpa3Gzi4nQWqyNH4OWX4aOPoEQJu95CKejLfHr3hkcyDu7eDWlpEoCFEDaTACxuDlYreHhAhQrQsyc88oiuNuQAJ09CfLwe3c520N9fD0ELIYQNZAhaFH1//63TPh4+rOdfp051WPAFuDBmIptok70GcK9eOgtWrVoOu68Q4uYiAVgUXefOwVNPQY8eYDLpkn/OsHEDZUggqMU1yTsMQxZgCSFsJgFYFE0LFkCjRjBnjk6usW0bNG2a6+mh4bG0n7CKmqMW037CKkLDYwt86zJHw4ks3oLSpdMPWCxw++26LUIIYSOZAxZF0/LlUKUKLFmiqxjdQGh4LKPn7yLZbAEgNjGZ0fN3AdCnecCNnnq9hAQqXj5CQqPBWccOHNClBiUBhxAiH6QHLIoGpXSZwM2b9e8TJ+qiB3kEX4CQpZGZwTdDstlCyNLIfDcjeUMEAEbLqxZbhYfr77IASwiRD9IDFu4vJgaefx4WL9Zzvm3agJ+fzU8/kZic5/HQ8FhClkZyIjGZKv6+jOhWP8fe8cGTxdlHP8p2uWq70fbtutZvgwa2/zcJIW550gMW7kspmDlTr3BetQo+/RS++Sbfl6nin3MSjozjGUPUsYnJKLKGqHOaJ16f1oZHmUuju8tnHQwPh6Cg/FdUEkLc0iQAC/f1008weLAe2t25E4YN06ud82lEt/r4emV/nq+XiRHd6gP5G6JO+XMp7/uMp1rshqyDDRrA/ffnu11CiFubDEEL92Kx6CxWderAo4/qXmW/fjrJRgFlDCXnNsRsyxA1AEuX8vLi7igMjC4+sHKlzro1bVqB2yaEuHVJABbuY98+ePZZnVAjMhJKl9ZB2A76NA/IdcVzFX9fYnMIwtmGrhcvRvXvD4AHClJTISwMWrXSRR5k/68QIp9kCFq4ntmsczYHB+vAGxICpUo57fY3HKI+c0bXD+7VizOqLCkUw+phAm9v6NgRPvwQAgJkC5IQIt8kAAvXSkjQq5rfegt694a9e2HAAKf2KPs0D2B83yAC/H0xgAB/X8b3DdI9ZouFK0tWE1Lyfape2s+s/qtR73+YNfwcHg5lysgCLCFEvskQtHANpXSQ9feHli11Nqu+fXM81dYtQoWRbYj62DH48nNOVRzPy69WZHHCIeo29WPtN9CqVTugXdYTw8N1T1gIIfJJesDC+davh9ats4onzJp1w+Br6xahQlu7Fnr1QjVogHnKNB5quIc//4R3PvJj61Y93ZtNfDzExkoJQiFEgUgAFs5z6RK8+irceacOXqdP5/kUe2axuqEff4S774bFi7FeSaVfyv8wNQtixw4YPTqXEeaMDFgSgIUQBSBD0MI5VqyAQYPg6FF46SW96KpkyTyfZvMWocKwWlHDhoHVigFYMRj9wH5az89j91NAALzxhgRgIUSBSAAWzrFggV45vGaN7gHbyKYtQhRwnnjTJmjalB0HfJlffCJvxg/BGzMmH29uf7Nj3uNDjRvrnNRCCFEAEoCF44SGQuXKulTfJ5/o7qRvzmkhczOiW/1slYwgexYr+H97dx5n9dj/cfx1zVZTKC1Uo372bi22u1JKJqQk2kjoxn0LqZCUuwWVtAndKC2IuN1KtCGVSgtNaKFtRJu0qEQLNc3MOdfvj+ukacxymjkz33Nm3s/HYx5zlu855/rOd875nOv6Xtfnk4dqR4cOQd++2NGjmXfV07RY+gTlyt1F40EXcE3UQkyTRDfDOfPuZAryg87xcU2rq6BkyZPaJxER0DlgKQi7d0P79tCmDYwc6W4rXfqkgy/kskQoIOjzxElJcM89cP752NGjebtMN9oufoSOHd3qp2ufaIDp2yfb4JtxMtiB3fu4psP1JPd48qT3SUQE1AOWULIW3n4buneHw4fded6ePfP9tDllsYIgzxMnJblJVmlp+DE8wFjmlb2fqe9B06a5tyFzkL9oz2YAXj98OhqEFpG8UACW0Jk82ZULvPJKV7WokMrz5Xie2FpISYGFC/Gn+4gCfETRutE+/jPbdcyDkTnI19ztAvCS0mflt/kiUkxpCFryx++HjRvd5VtucT3gJUsKtTZudqkkn6pdGm68kZTb/8kT8xJJsSVIJ5qoEnHc+Gxi0MEX/jrpq+buzewtXZaYs0KbEEREig8FYMm7776Dxo2hYUM4cMAVJejYMV+Vi/Ii83niqqfFMTn1a66/7VrSFixm4NwGjFhSn/funw+DBhH92fwsz/PmJHOQr7lnE8mVzqdX88L7oiEiRYuGoOXkpaW5ggkDB7ox3JEjC7V4QlZaX5ZA65RtMPUTeHc2rF3L1+Wbceu+sZzV8Gy+fQ3+9rdMaSRP8vnheEnDV1p147Yrzg55SkwRKT4UgOXk/PYbNGkC337r6vS+9BJUquR1q9wkq2uvxaamgs/P0JinGJoygOGjDZ07h6ZTfuJksBvz/4QiUqwpAEtwMhZPqFcP+vd3y4zCwbJl8NBD2KOpGL+PNKJJOK8k6z81VK1aAK+3cqXLAd2iBURH5769iEgWdA5YcrdokUu3uGmTC8Ljx4dH8A3klrZXXsnvydtI8ceSRjQmLo67JiQWTPAFVzyiY8dCLZkoIkWPArBk78ABeOABV27v0CH49VevW3Tc7NlQsyb25Zf5X5kuVD68iRE3LCC13yBiFs7HXJm3c71BWbUKLr200CebiUjRoiFoydrMmfDgg/Dzz67gwNNPQ6lSXrfKSU3F36Ubew6Wop1dwo4yDXl/EjRrlvdJVkHz+WD1aldYQkQkHxSAJWtz50KFCi6fc926BfpSQRVSWLoURo+GTp346I8mDD/8CV/vr8aD3UswaBCcckqBNvG47793Wb5UAUlE8kkBWJxjaSQvvBDq13fFE2JjsymEGzpBFVKYNs0l+fD7SX93CoPtIn6v1YBFM1ydh0KlGsAiEiI6iSWuRm/z5i6N5Pjx7rZSpQo8+EIuhRT8fhg9GtuhA9bvB8BaP4OvW8iKFR4EX4DbboN166BGDQ9eXESKEgXg4szng//8x9W1XboURo1yM3wLUY6FFLp3h27d+K7kpaRQknSiiS4ZxzVPJxIXV6jNPC462gXfGA0eiUj+KAAXZ2+/DY8+6mY5r1sHXbsW+szezDmWY31pnJbyO5VPi+etUp25L24i9XzL+LjHAqIGDyJqwcmnkQwZa93fa8kSb15fRIoUfY0vbo4edROJatd2a1nLl4eWLT1b09qrWXX6TF3DRVvX0mbtZ1y1ZSVry9egR6lPuHttHC1a1GDdGKhWrRBmOOckLQ2ef96NGFSvDldd5V1bRKRIUAAuTpYuhU6dYN8+2LzZ5XG+6SZPm9T6sgQqfPU59Yf0JtrvwwKvHOzBb+Xj+N//oEOHMMh3sWyZWw+9ejXcfDPceafHDRKRokBD0MXBoUPw0EPQqBH88Qe8+WbwhXAL2sqVNOrdmRi/DwP4iKbOJekkJ8Ptt4dB8F21ytU33rfPzcaeMQNOPdXjRolIUaAecFH3888ud/P27S4IDx5ciItmc/d7mQT2R1WlAkeIIZ2oEnG0fyURKnjYKGtdjeMLLnAZr155Be64w/OKTyJStKgHXFSlpbnfZ57pqhZ98QW8+KL3wddamDQJ2rTh4w/9XJR4JtV+W82YWz/D1z9vtXpDautWNyx/ySXw44+uC965s4KviIScesBFzbGEGv36wWefwfnnu8lDXktKclm1vvgCvviCTeXq8s/pv3BGzTOYMgXq1/d4klV6uptg1b+/C7qDB0OCav2KSMFRAC5KtmxxvbW5c915S2u9bpHzxRfQpAk20Ct/tcRDPHxwJH0GRNOnD96t6T3m6FH391q50vV+R42CatU8bpSIFHUKwEXFiy9C375uHe/o0YSsCn0ozJ+PTUvDAOlE4zujMis+iaZmTY/blZrqon+JEi7w9uvnyix6PvNLRIqDMPmElnzbuBGuuQbWr4cuXbwPvqmp8OKL+A4d5r3fmpJCSdKIxsbGcf//Er0PvtOnu+H5pCR3fcAAaNtWwVdECk2uPWBjzASgJbDHWlsrcNsA4D5gb2CzvtbaWQXVSMnCkSMwaBDceCM0bAgvvODSI4ZDAPn6a7j3XlizhmdeLseATf+gadVpXFP5fTZffgktSlejtVdt274dunVzy4kuvjgMxr9FpLgKZgj6TWAU8Fam20daa58LeYskd4sWuXq0P/zghk8bNiyUwgm5+uMPePJJ7Isvcqh0Je6JnsG8fTdSudW3bKju43vTBoAlmasdFZZx46BnT5cDe/hwl1YyHP5uIlIs5TpOaa1dDPxaCG2R3Pz2mwu8iYkuiHz6qZu1W4Cmr9pBw2ELOKf3xzQctoDpq3ZkvWFSkpvINHIkk8vcT9VD6ynV4Wb+1nUpcX/bfkLH/M9qR4Xt4EGXjGTdOnj8cQVfEfFUfk4UdjPGrDbGTDDGnB6yFkn2/vtfeOMN6NUL1qyB664r0Jc7Vqt3x/4jWI7X6j0hCO/bB7NmYa+9Fv/qtaRQgkmxd/Hux2X4739hb/rBLJ87uypIIfXHH+5vNWWKu/7YYzBrFpxzTsG/tohILvIagMcA5wGXAruAbBeaGmPuN8YsN8Ys37t3b3abSXZ27IDFi93lBx90qRGffdbV6y1gOdbqtRYmT4YaNTjYrQ++I6lE4SfWpDP5wYW0aOG2z1zt6Jjsbg+ZTz5xZRafe879zcBNTAuHc+QiIuQxAFtrd1trfdZaP/AqUC+Hbcdba+tYa+tUrFgxr+0sfvx+GDPG1Z696y6XKCImxlUxKiTZ9VL9236CVq2gQwc2p1el65bHSDNx+KNcvd4SzRL/3LZXs+rEx0af8Pj42Gh6NateMI3++WdXwaFFC/clZfFiGDKkYF5LRCQf8hSAjTGVM1xtA6wNTXMEcEuJGjd2y4nq1oX58z0pAJ9VL7Xh1m+YN6EL6XPm8VSp56h5cBnn9b+L6M/mE/XMINfWDKkkW1+WwNC2tUkoG48BEsrGM7Rt7YKbgLVkiSuaMHCg6/mqbKCIhCljc8mWZIx5F0jEpcffDfQPXL8UsMBW4AFr7a7cXqxOnTp2+fLl+WpwkZec7PIQn3qqW1p0112eDZseOwd80da11N+2mmXVLmZ7ifMYMv0tHtj3HBWvOI/XXoNatTxp3nHJya5U4G23uaHx7duhalWPGyUiAsaYFdbaOlnel1sADiUF4Bzs3u0KJ1gLI0dCx45wxhlet4rFr0+lwf3tifH7SIkqyQ2x8/g6uiFDhrjltNHRuT9HgUlJgaFD3c8ZZ8CmTW5ZlohImMgpACsTltcOHHBDzeee67JZGQM9eoRF8GX5cho/0ZXYQK3eGH8ad1VdzLp18MgjHgffRYtcqcCnn4b27V0eZwVfEYkgCsBemj7dTbIaNw7uvx8qVfK6Rc6RI9CzJ/aKK/jjYDpHiSONaExcHP+cmMjZZ3vcvi1bXNrN1FSYPdstzwqHLywiIidBxRi84Pe785Xvv+/SIU6f7iZbZTJ91Q5GzNnAzv1HqFI2nl7Nqod88lKWr/G3chyZ8iEflrmP+34bTvem6+lVdyGntEz0rlavtbBiBdSp49bxfvABXH99oSzHEhEpCArAhclaN8QcFeWGnIcNc8PNWWRkOjYB6tg63GNJMCB0KRwzTrJqv3kF1fb/TP89PXj95yYs3LaCsmedwjtvQcuWHtfq3brVDdPPnu3yTP/979Das2zSIiIhoQBcWNavdyUChwxx6RCHD89x85ySYIQqAI+Ys4GLtqxl0rt9iPWnA/DFxmY8n3oTXbq4uU2nnRaSl8qb9HR46SV48kn3xWXkSHfeV0SkCFAALmhHj7pINmSIW1q0b19QD8suCUYoUzj6fvqJ4bNfIi4QfNOJpkR0CpXuXMro0VeG7HXyxFp3nnfJElfx6ZVXoFo1b9skIhJCmoRVkD7/3PXYBg50M3WTk10GqSAURgrHF+aPodqvP5NKLGlEkxoVy/rWp3FurZSQvcZJO3Lk+FB9x47w3nvw4YcKviJS5KgHXJC++sqtVf3kE2je/KQe2qtZ9RPOAUOIUjh+/z2UKcP2tDMZceprrLElqVJ+M83Pfp/lF13Id2fXYGhBpYnMzezZbph++HA3Se3++71ph4hIIVAPOJSsdZV3pk931x9+GNauPengCwWQwjEtDYYNw158Mcmte1OjBsxdX4vGPSqQ3iOGV69rxe6alxdsmsjs7NkDd94JN9wAJUvCWWcV7uuLiHhAmbBCZds26NoVPvrIBZJZs7xu0XFvvAH9+sGuXSws347b971M7aaVGTcuDCrzTZnier2HDkHfvtCnjxJqiEiRkVMmLA1B55fPB6NGuQBnLTz/vOv5hounn8b27w9AKnEMTX2MYW9W9jLF9F9ddBG8+qr7LSJSTCgA59e8edC9uxtmHjMG79NEBaSkQMmS7NxylEoYorDEGB/vP7SQU+8O7Zrek0oYkpbmvqSULOn+brfcAu3aubXRIiLFiD718uLwYfjsM3f5+uvd5VmzwiP4/vYbdOqE7+omPNbdx60TW3KUkn/W6j21ZWJIX+5YMo8d+49gOZ4wZPqqHX/d+KuvXCarPn1g+fITE5OIiBQz+uQ7WXPmuPp7N94Iv/ziAkhiYniM537wAdSogf+NN3ltQ2NGvZjOxQ80wDc361q9oZBTwpA/HTrkerv167t10NOmufzN4fA3ExHxiIagg7Vnj0sb+c47UL26W1pUoYLXrYKkJNf7XrwYFi/mx3KX0sb/MX+ceTmfzoTGjQEaQNOCSSUZVMKQ9evdefIHHwyD9FoiIuFBATgYBw64Xu/+/dC/f/jM1E1Kgmuvxaamgt/PGyUfpOuBF+nRN5Ynn3SnWQtalbLx7MgiCNeMPgJvvgn33ANXXOFKLYbDEL2ISJhQAM7J3r1QsSKUKQMDBkCTJuEzU3fjRujSBXs0FeP3kUY0h8tVZdmsWC65pPCa8ZeEIdZyx/oFDFg8AVKPunPkVaoo+IqIZKIAnJXUVJeNacgQ+PRTVzyhS5eQv0yeyg2mp8Pzz2MHDCDNH43fH0M0QGwcnSclElOIwReOV2YaMWcDsVs28dz8MdTZtMr9zcaPd8FXRET+QgE4s88/dykQk5Nd/ubzziuQl8lTucFVq+Dee2HVKpaUb0OHfaNoV+dHBly9kPLtEj2r1dv6sgRaVz8dqrV3y4zGjHF/Q81uFhHJlgJwRj16uJJ3//d/8PHH0KJFgb1U0OUGk5Jg4UJITMTftx+Hf9hFp5j3meNrxwsT4J57qmCMh7V6N2yACy+EUqVgwgRXqzehkFNZiohEIAXgY6k4jXGB97HHXPWi0qUL9GWDmj2clOTOO6en44uJ47FK7zDx90Suu+V0kl+GSpUKtIk5++MPNyFt5EiYNAluvRVuvtnDBomIRJbiPUa4datbz/vuu+76I4/Ac88VePCFIMoN7t/vUloePQo+H/6jqZzx63e8Me10pkzxOPjOnQu1a7uMVvfdB02betgYEZHIVDwDcFqaC7Q1a7r1s0dCV+Q+WL2aVSc+NvqE2/4sNzhtGtSogV2xgnRiSCMaGx3HQ+8n0rp1oTf1RD17QrNmEBsLixbB2LFQtqzHjRIRiTzFbwh6+XLo1Am+/dYNmY4aBVWrFnozMs4ePmEW9OdT4eGH2VbuEtrYDznvrFSGNV/Iuf9KJM6jSVZY636iouDKK90C4yeeKJyFxiIiRVTxC8Dbtrn1vVOnQuvWnqZDbH1ZggvE1sKvv2LLlWfmyg6sOuUoQ/c/wqO9Y3nqKYiP93CS1Y8/ugxWV18N//43tG3rfkREJF+KxxD0tGluaQxAmzbw/ffudzjkIp4yBc4/n7QrG9OutY9WnSoy88KeJC2PZehQiM/6VHHB8/ngpZeOD9MrfaSISEgV7R7wTz/BQw/BjBkuHeIDD7hh1BBMsspTEo2M0tOhe3fs6NEAWGL5dcsynn22IY8+CjFeHpnkZPjXv2DZsvArsygiUkQUzQDs88Ho0dCvn7v87LOuGk+IEkMEm0Qj2yC9bZvrga9cCYABovDzwcOLKd+rYUjamC8HDsDmza7wxO23h8dIgYhIEVM0h6DXrnUBt1EjWLcOevVys3ZDJJgSfDnVyU0rW5Gffi1N/+hnOEK8q9UbH+eyWYXY9FU7aDhsAef0/piGwxZkXacXYOlSGDbMXa5f3y3RuuMOBV8RkQJSdALwoUPw3nvu8iWXwIoVrkzfOeeE/KWCSaKROUhfsW0N4/7bl2HjNlHv6niqbV3E2lb9ODzT1eo1BVCrN6cvAX86eBC6dXNfVsaNc39H8PDks4hI8VA0hqBnzHBBZOdOqFMHzj0XLruswF4uuxJ8GZNrHAvGjbasoufit7j05x/YUqIqB8aVZ38l+OADE5hM3ABuKphZzrmmu/z4Y+jcGXbscEk/nnkGTjmlQNoiIiIniuwAvH27m2Q1fbqr1/veey74FrC/lOAjQxKNgCpl42k5+216L3oTgDRi+OfRifxa93Q2zA1N7orcJoLl2FPfs8cVmzjnHDcTu379/DdIRESCFrlD0CkpULcuzJnjzl2uXFlo1YBaX5bA0La1SSgbjwESysYztG3tE4Jfr+svpOOqTwA3yQosN1w6hXHjbMiCb27Dy39Jd2ktjbasokqZknDGGTB/vvu7KfiKiBS6yO0BlywJr7zizvfmodeb32VEfybRyMhamDgRmjXDt+UsHrOv8ia3EUsqvpgY6j/SnKtPZqlSDoKpppSxp55wYA9D5ozi6i0rSWo00T1AgVdExDORG4DBLeXJgzzV4s3N5s2uBu78+Uz+21N0+G4gl156Ez89Mp/quxYSm5jI1SHsoQczEaz1ZQng87F14HDumzMBYwzfPj6IBl07hqwdIiKSN5EdgPMo6Fq8uUlKggUL4Oefsa+/TpqN5fGSYxm/5T6GDXPlhWNjGwChHxoPZiIYQOsh3eGjaXDDDTB2LJdUqxbytoiIyMkrlgE4qFq8uUlKgmuvdeeirWXNaY244eAkLkxM4NvxcMEFIWpsNnKcCJaa6tbvxsbCP/8J7dppTa+ISJiJ3ElY+ZBrLd7cpKTAjBnY1FSwFh9RTDt6AwPGJzB/fsEHX8hhItjRn+Dyy12tXoCbboI771TwFREJM8WyBxzMMqJsLV4M993HkfRojD+OaFLxR8fR5b0mVLy5ABudhRMmgv3+uysR+NJLkJAAF19cuI0REZGTUiwDcLa1eHM6/3vgAPTuDWPHsq/MOdxxaDSxp5dmWPOF1OqWSEWvavUCLFkC//iHKx3YpQsMHarqRSIiYa5YBmDIZhlRdtatg2bNsLt2MaFsDx7e/zS331uaESPg9NM9DLzHxMW5Ck9LlriUkiIiEvaKbQAOytKlsGgRh2o34PvY+jzof5x95erx4QdwzTUetstaeP99WL0aBg1ypRbXrAlZtScRESl4CsBZsRaeegoGD8Zvoojxx/GwmU9ir3oMGAClSnnYtp07oWtXl36zTh1XcrFkSQVfEZEIo0/tzLZsgWbN4JlnsNYS5fcRSyrvdVnIs896GHythddegxo1YPZsV+M4KckFXxERiTgKwMf4fDByJLZWLVKXLOOluJ4cIR6fcbV6E+5M9LZ9O3fCI4+41JurV7saxzEawBARiVT6BD/GGI5MmsE38dfQft8rnNu4Kq27tKXa5oWQmFhohR5O4PO5oea2bd3Soi+/dD1gDTeLiES84h2AFy6EZ54hvevDjNx4M89++yGpcacwYpyhUyeIiiqYNJJBWbsW7r0XvvoK5s1zWbdq1fKmLSIiEnLFNwCPHevWzFqLXbCYqXYRDVs1YPRo19n0TGqqW8c7eDCUKQP/+5/HU65FRKQgFL8AfPCgS6gxZgyWQK1e62f87Qup9U4D7zM2tmwJn37qcjf/5z9QsaLHDRIRkYJQ/E4mPvkkduxYPjql/Z+TrGLi46j9UKJ3wffwYUhLc5cffRQ+/BDeeUfBV0SkCCuyPeDpq3b8mWqyRkwK3eueSb3ERgzc/yRf2DvYW/EKJg9Oou4fC72bZAXw2WfQqZP76dPHlQ0UEZEir0gG4OmrdtBn6hou2rKWvsunc9WWVSSfdiHnmuXs31eBRx+rwMCBULq0h5OsDhyAxx+H8ePhvPO8+wIgIiJghAeaAAALJklEQVSeKJIBeMScDTReu5hXZgwj2lp8GP7zSy9SKh1m2bJTqVs39+fI2IMOqljDyfjsM1c8Ydcu6NkTBg70OL2WiIgUtiIZgCutXcnLM58lyloA/ERx0dnLSbqlLHXrtsj18cd60MfKFe7Yf4Q+U9cAhCYIly4NFSrA1KlQr17+n09ERCJO0ZqEdfQoADvP+DuzSjQjhZKkEU16dAzfNKpCQvng0jaOmLPhhFrBAEfSfIyYsyFv7bIWJk2Cvn3d9Xr1YOVKBV8RkWKsaPSAjx6FwYOxkybx8j0rWTnmatr5r+K6epO4puSnfFmtNsln12Jos+pBPd3O/UdO6vacn2wnPPggzJzpAm5KiooniIhI7gHYGDMBaAnssdbWCtxWDpgMnA1sBdpba38ruGZmISnJZbIqX96tl01O5pNyHXmqXzrNbori5i67eeObKozd354qZeMZehLncKuUjWdHFsG2Stn44NtnLUyYAI895r4gjBgB3bsrf7OIiADB9YDfBEYBb2W4rTcw31o7zBjTO3D936FvXjaSklxqxpQUrLUcKnUmt0d9wvKY5oyfDLfeCsZUplPzynl6+l7Nqp9wDhggPjaaXkH2oIHjxRP+/ndXxeiCC/LUFhERKZpyHQe11i4Gfs10cytgYuDyRKB1iNuVs4ULXcpGa7EYXjjcmYr/aM769dC+PflOqNH6sgSGtq1NQtl4DJBQNp6hbWvn3oP2+2HaNNf7TUiAZcvcjGcFXxERySSv46FnWmt3AVhrdxljzghhm3KXmEhaVBz4Ukk3cdzwQjOu6B7al2h9WcLJzXjesMEVT/jiC5gzB66/XsUTREQkWwU+E8gYc78xZrkxZvnevXtD86QNGpD88nzmNhwE8+ZzRXcPk1ikp8Pw4a5O7/r1MHEiNG3qXXtERCQiGBtYK5vjRsacDXyUYRLWBiAx0PutDCy01uZ6grROnTp2+fLl+WtxuLn5Zpe7uV07GDUKKlXyukUiIhImjDErrLV1srovrz3gmcDdgct3AzPy+DyRKTX1ePGEzp1hyhR4/30FXxERCVquAdgY8y6QBFQ3xmw3xtwLDAOaGmN+AJoGrhcPX3/tZjaPGOGut2gBt9zibZtERCTi5DoJy1p7ezZ3XRvitoS3I0egf394/nmoXNmd8xUREckjZYUIxldfQceO8MMPcN99rvdbpozXrRIRkQgWkQG4QCsVZcUYt7Z33jyXAERERCSfIi4AF3ilomPmzXNrevv3h7p1ITlZaSRFRCRkIq4iQMgrFWV24IAbZm7aFN59F37/3d2u4CsiIiEUcQE4pJWKMvv4Y6hZ0xVRePxxWLUKTjkl/88rIiKSScR160JSqSgrv/wCt90GZ5/t8jnXrZu/5xMREclBxPWAezWrTnxs9Am3nXSloow+/9xNsKpQAebPhxUrFHxFRKTARVwAznOlosz27nU93quugunT3W1XXAElSoS8zSIiIplF3BA05KFSUUbWwuTJ8NBDcPAgPPMMtGwZ2gaKiIjkIiIDcL506QJjx7ph5jfecJOuREREClnxCMDWgt8P0dFw441w7rnw6KNaWiQiIp6JuHPAJ23HDlcycFigXkTLltCrl4KviIh4qugGYGvh9dehRg03u7lsWa9bJCIi8qei2Q3cts1ls5o7F66+2gXi887zulUiIiJ/Kpo94N274csvYfRoWLBAwVdERMJO0ekBb9kCH33klhfVrQs//QSnnup1q0RERLIU+T1gvx9GjYLateGJJ1zvFxR8RUQkrEV2AN64EZo0cb3eRo1gzRo480yvWyUiIpKryB2CPnIEGjaEo0dd9aJ77gFjvG6ViIhIUCI3AMfHu0xWl1wCCXlMSykiIuKRyA3AAC1aeN0CERGRPInsc8AiIiIRSgFYRETEAwrAIiIiHlAAFhER8YACsIiIiAcUgEVERDygACwiIuIBBWAREREPKACLiIh4QAFYRETEAwrAIiIiHlAAFhER8YACsIiIiAeMtbbwXsyYvcCPIXzKCsAvIXw+L2lfwk9R2Q/QvoSrorIvRWU/IPT78n/W2opZ3VGoATjUjDHLrbV1vG5HKGhfwk9R2Q/QvoSrorIvRWU/oHD3RUPQIiIiHlAAFhER8UCkB+DxXjcghLQv4aeo7AdoX8JVUdmXorIfUIj7EtHngEVERCJVpPeARUREIlJEBGBjTHNjzAZjzEZjTO8s7jfGmJcC9682xlzuRTtzY4ypaoz5zBiTbIxZZ4x5JIttEo0xB4wx3wR+nvKircEwxmw1xqwJtHN5FveH/XExxlTP8Lf+xhhz0BjTPdM2YXtMjDETjDF7jDFrM9xWzhjzqTHmh8Dv07N5bI7vq8KWzb6MMMZ8F/j/mWaMKZvNY3P8Xyxs2ezLAGPMjgz/Ry2yeWzYHJds9mNyhn3Yaoz5JpvHhtsxyfLz19P3i7U2rH+AaGATcC4QB3wL1Mi0TQvgE8AA9YEvvW53NvtSGbg8cPlU4Pss9iUR+Mjrtga5P1uBCjncHxHHJUN7o4Gfcev2IuKYAI2By4G1GW57FugduNwbGJ7Nvub4vgqTfbkeiAlcHp7VvgTuy/F/MUz2ZQDQM5fHhdVxyWo/Mt3/PPBUhByTLD9/vXy/REIPuB6w0Vq72VqbCkwCWmXaphXwlnWWAWWNMZULu6G5sdbustauDFw+BCQDCd62qkBFxHHJ4Fpgk7U2lMliCpS1djHwa6abWwETA5cnAq2zeGgw76tCldW+WGvnWmvTA1eXAWcVesPyIJvjEoywOi457YcxxgDtgXcLtVF5lMPnr2fvl0gIwAnATxmub+evQSuYbcKKMeZs4DLgyyzubmCM+dYY84kxpmahNuzkWGCuMWaFMeb+LO6PtOPSgew/TCLlmACcaa3dBe5DBzgji20i7dgA/As3opKV3P4Xw0W3wHD6hGyGOiPpuFwF7LbW/pDN/WF7TDJ9/nr2fomEAGyyuC3z1O1gtgkbxphTgA+A7tbag5nuXokbAr0EeBmYXtjtOwkNrbWXAzcAXY0xjTPdHzHHxRgTB9wMTMni7kg6JsGKmGMDYIzpB6QD72SzSW7/i+FgDHAecCmwCzd8m1kkHZfbybn3G5bHJJfP32wflsVt+T4ukRCAtwNVM1w/C9iZh23CgjEmFnfw37HWTs18v7X2oLX298DlWUCsMaZCITczKNbanYHfe4BpuGGajCLmuOA+JFZaa3dnviOSjknA7mND/YHfe7LYJmKOjTHmbqAlcKcNnJDLLIj/Rc9Za3dba33WWj/wKlm3MSKOizEmBmgLTM5um3A8Jtl8/nr2fomEAPw1cIEx5pxAL6UDMDPTNjOBuwKzbusDB44NKYSTwDmT14Fka+0L2WxTKbAdxph6uGO0r/BaGRxjTGljzKnHLuMmy6zNtFlEHJeAbL/NR8oxyWAmcHfg8t3AjCy2CeZ95TljTHPg38DN1trD2WwTzP+i5zLNf2hD1m2MiOMCXAd8Z63dntWd4XhMcvj89e794vXMtGB+cLNpv8fNQusXuK0z0Dlw2QCjA/evAep43eZs9qMRbthiNfBN4KdFpn3pBqzDzbJbBlzpdbuz2ZdzA238NtDeSD4upXABtUyG2yLimOC+NOwC0nDf0u8FygPzgR8Cv8sFtq0CzMrw2L+8r8JwXzbizr0de7+Mzbwv2f0vhuG+vB14H6zGfXhXDvfjktV+BG5/89j7I8O24X5Msvv89ez9okxYIiIiHoiEIWgREZEiRwFYRETEAwrAIiIiHlAAFhER8YACsIiIiAcUgEVERDygACwiIuIBBWAREREP/D9SIhgJs+ctbwAAAABJRU5ErkJggg==\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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/numpy/core/fromnumeric.py:2495: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.\n",
      "  return ptp(axis=axis, out=out, **kwargs)\n"
     ]
    }
   ],
   "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:                Mon, 24 Feb 2020   Prob (F-statistic):           4.98e-10\n",
      "Time:                        22:49:06   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",
      "Warnings:\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:1449: 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\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": "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/scipy/stats/_distn_infrastructure.py:901: RuntimeWarning: invalid value encountered in greater\n",
      "  return (a < x) & (x < b)\n",
      "/usr/lib/python3/dist-packages/scipy/stats/_distn_infrastructure.py:901: RuntimeWarning: invalid value encountered in less\n",
      "  return (a < x) & (x < b)\n",
      "/usr/lib/python3/dist-packages/scipy/stats/_distn_infrastructure.py:1892: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= _a)\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:761: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_external * np.sqrt(hii / (1 - hii))\n"
     ]
    },
    {
     "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"
     ]
    }
   ],
   "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.2rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
