{
 "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:                Tue, 17 Dec 2019   Prob (F-statistic):          2.83e-239\n",
      "Time:                        23:42:53   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:                Tue, 17 Dec 2019   Prob (F-statistic):           6.30e-27\n",
      "Time:                        23:42:54   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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzddXiW5RfA8e+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/t3Xl8VOW9x/HPw44iBETB4IKIpci1FU3VaNXYui8FvZWC4nbBpSouVxEoFUGs+9Zar14VbbWKvKQomygIBEFRQBYFUQm5oiwGlFVBIMlz//hNIIRMGJIzc+ZMvu/XK68MM2fOeR5O5vzmWc7zw7PzS0JSxrLHjIF774UJE4Lft4hIGqpd+YA3bYKf/cwWdZg5E+ok7/tHxRYuWCu6LJAf3n88py2dw3OjhlKv1LZxwKetjuDobwvsDYMG2Qzt4mLLgTt58q7d0LE8vDty5JZ7fe1aS2X8+efww6SZlE7JZ8TqPD4kl44dLbZfeKH1WMdLHXjyA1Mq7Ypvk9WY9/v/JpD/JwCKimwWena2fclQvl8RyRBV5QOuXV3Q995ri0yMHp3U4As7u5Hjdd9mZzVmQP4L1I8F31Jg7M9P4aULr2Nk2U7uucfSIsYZA96ek8uqFyez+a185jXLY+qLuXzez4LumjU7t2vQIJeTTsqlW394+UJo3z6xOgQxGW2PvIerr7YvR6++quArIrVG7QnAX3wBjz9uF/vjj0/JIbt2brNrt/UHH8DFfeDhh+l7dgee+7A79457nHqlJWyvW4/hJ15Ezx55u+yj5PhcljTPpaAACp6AJUssFW9BgS3eVVKSC1hgbtnS5pR16WK/y37ato3fyq1KEJPR9mj4cOt2fuopOOqo4PYrIpLmak8AfuopmyV8//2pPe6MGTYWu2ABLFoEzZvD4sV0veACuPtmbj7oYNp/NoeCo3L4w7UX0+WYNnz5pc2bevddmDoV1q/fubtmzeDII+07xKWXWmu2fXvo0MECcJCCmIy2R5dcYpPievYMbp8iIhFQe8aAi4stAP7yl6k75pQptpqT9zZZ6rbbrFu5QraloiIbvi0Lut98Y88fdhiceaYth9yhgwXaFi1SuzZF0mZBb9kCP/4Y/LcGEZE0UrvHgLduhc2breWZyuAL8NFHOx/XqWPBJhZ8vbegO2SINZLBivib39gCGGeeCe3aVR1sU3GL0G7d6EHp1w9GjbIvRc2aBb9/EZE0l/m3IT3xhM18/vbb1Bxv82bo1cv6jvPyoFEjG4Bt0ICyLAbTp1uK3TPPhK++gr/8BWbPtolTI0faSoxHHLHn4JuyW4SC9tZb8OST1v2s4CsitVRmt4BXroShQ60buHXr5B9v6VJbKGPBAptQdPvt1syNzWKeVTeXu8629TVat4a//c3S7TZqtPeHSmSlrbRUVGQT4X7xi9SPx4uIpJHMDsD9+tnY72OPJf9YY8bAFVdYV/P48XDeefZ8bi7zG+cyaBCMHWu90A8/DDfcAPvsU/3DpeQWoaCV3XK0caONj1fnm4eISIbI3AA8c6alGxw40AZTk3WM/HzIyrKIetxx1ofcti0A69bBH/8II0bYJvfeawmYKltBcm+l5BahoP34o/WrP/IIdOoUdmlEREKVuQF47Fjr5+3fPzn7r7gKVd++NsM51qorKLClpgsL4c9/tt7orKzgDp+SW4SC1qQJjBsXdilERNJC5k7Cuu8+mD/fLvrJMGmS3UpTUmJBuHnzHcF3+nRb4nHNGrutaOjQYIMvBLNedcqsWwdXXWUzzpxTjl8RETKxBbx9OyxfbrlkW7VKzjFKSmxyFdiYb7kZzi+/DL17Wy/0+PGJL/tYHUm7RShIW7fCRRfZKmBXXbWje15EpLbLvBbwsGG2asWiRcnZv/dw663w3nv2+957YfJkSk/I5a67bB7WySfDhx8mN/hGQmmpTbqaNg3+8Y8dX1JERCTTWsCbNllKvhNPTN66wo88An//uw3qPvIIYD3RV19qk6169YL/+Z9gcgqkNBdvMgwcaGs933+/rZspIiI7ZFYAfvRRWL3abglKxjhjcbH1K3frZmkCsdtau3a1Ra8efNDmYgVx6IrpDMsW2gCiEYR/+MHOw3XX2e1gIiKyi8wJwKtWWYv0kkvghBOSc4x69eDtt+1xnToUFNgaH6tXw7//bUOdQYnsQhtlmjSxcd9999WkKxGRStRoDNg5l+WcG+mc+9w5t9g5l7vndyXJtGk25njffcHve+FCy2C/dq3NdG7UiA0b7DajH36w4eAggy9EdKENsDU1r7oKfvrJlpmslznf8UREglTTq+Nfgbe99793zjUAarC2Uw11726LK++/f7D7Xb4czj3XgvsPP0CLFpSUQI8etvLk5MmQU2mei5qJ5EIbhYX2rWSffWDDBq10JSJShWq3gJ1zTYFTgWEA3vtt3vv1Vb8rST77zH4HHXwnTbLounatJRA49FDAshVNmACHXvA5V741npMfmBJ4EoS+Z3egcf26uzyX1gttfP+9fVHZvt3+c5J1C5iISIaoSRd0O2AN8KJzbp5z7nnn3L4VN3LOXeucm+Ocm7NmzZoaHC6OGTNsWcMRI4Ld7wcfwDnn2CyrkhLLcoStbvnQQ5B13DJKOixNWiaiSC20sXEjdOkCy5bB6NHw85+HXSIRkbRXky7oesCxQB/v/UfOub8C/YG7ym/kvX8WeBYgJyfH1+B4u/Peph1nZ8OFFwa6a8aNs25nsNnP+fnMqptL797QtN06mp6+633GyZggFYmFNsAmwC1ebKuQnHJK2KUREYmEmrSAlwPLvfdlWedHYgE5dUaNshUv7rmnZqmFKnPhhdC48Y5cvms65XHRRXDQQdDs/Nm4urt/l0j7CVJB++QT+xLUoYON/15ySdglEhGJjGoHYO/9t8A3zrmyQcnfAp8FUqpEbN9uiRY6dYIrrwxuv6Wl8Mwz0LmzzbAaOpStb03mgr/ksmGD3dp6SHblHQdpPUEqSN5bTsXOna3VCzbjWUREElbTWdB9gFdiM6ALgatrXqQELVwI331ng7LVuNUl7ipTTz8NN91kAaVHD/yJuVxzJcyaZQ3uo4+GvsURzEQUlM2bbbHr4cNtQZL//M+wSyQiEkk1CsDe+/lAEm7CSUDnzpZdp2nTvX5rvFWm9v26kDPvvNNm83bvDtjiWi+/bL3cZff6lo3LRnqZyOr4+mtb9mv+fLvfun9/LbIhIlJN0V4loZrdnpWtMrV16zZa3XIHNGwIzz8PzvH227aK4iWXWE7f8iIzQSpICxfC//2f5Vo+//ywSyMiEmnRDsDVVNlkqd6z3+QXyxbBK69Adjbr1lkin06d4MUXa3FDb9UqmDIFLrsMzjvPAnDQyY1FRGqhWhmAK1tlamq7HLLdNq7q0QOwlu/q1ZZ7Yd/d7m6uBZYtsxuehw2zSVdnnQUHHKDgKyISkMzLB5yA8qtMOW/3+i7PbkfWow+Cc0ybBs89B7fdBsem9saq8K1caTkV27e3/4Qrr7R7fA84IOySiYhklFrZAi4/iarb2Ofp+GMRW54dRpfObfjpJ8ug17YtDBkSbjlTatu2nUmMR42CG26AO+6AQw4Jt1wiIhmqVgZgiE2imvs2zHzNlpw8vi1gk3u/+ALeeSfDu56//96W2/zgA1vOs359G+vNzoYVK4Jf2ERERHZRK7ugAZg+Ha65xsY38/Nh5kwWLYIHHoCePW3IM2OUlsKSJTv/feON0LKlpVh89FFb1OS002zJTVDwFRFJgVrbAubBBy34AmzbRunUfK4Zl0vTpvDYY+EWLWGlpbBmjY3brloFp54KTZpY8/355+Hbb+1n5UpbQKOoCA480NI2HnoonHSSZXtqXEtW8BIRSSO1MwAXF8PHH0OdOnZ/UYMGvLE2j5kz4Z//TKP5RqWlsGgRLFhg3cJ/+IMNTo8bZ2O0q1btbLUCzJkDxx1nQXnhQmjdGn71K/t99NE78/N27RpKdUREZKfaGYDr1bOgNmMGfPEFq4/K4+rLcznjDLj88rALh3UX9+sH771nY7VlOna0ANy6NZx+OrRpYz/Z2ZYloiwNYM+e9iMiImmr9gXgdetsBa0DD4SLL8Z7uPYia0g+80wIC258+aXdbJyfb2OyvXpZN/L8+ZaRKS8PTjjBZiOXzQrLybGmuoiIRFbtCsDew+9/b2Oe48YBdsfN6NE2JHzEESksy8KFMHQovP66lat9e5uNDdaaLSxMYWFERCTValcAHj3abrX5+98BWL8e+vSBY46B//7vFJflmmssCPfvb+O5Bx+c4gKIiEiYak8A3roVbr/dFne+7jrAYl9RkeX4rUZGw70zb541s5980mZ5vfACtGoFLVok+cAiIpKOMjYAV8z3++yqd+lUWAgTJ0K9eixYAM8+CzffbEOqSTN3ri2pNWaMjT3Pn2+3AXXsmMSDiohIusvIhTjK8v2uWL8FD6xc9yN1Roxg1WlnWfADBgyweDhoUJIK4b21eHNybDbzkCGWvzh2fBERqd0ysgVcMd+vd3Xo2vMR2jcqZTwwdSpMmGDJfpLaAzxvHnTrBv/7v9XOXSwiIpkpIwNw+Xy/2RtXs7ZxU36q34jPiq1h2q+fzXm66aYkHHzZMigpgXbt7FahBg1qcTJhERGJJyO7oLOzYksres/fxjzMa8MHgPdkZzVm5EiYPRvuuScJKzBOm2ZdzpdfbpG+YUMFXxERqVRGBuCyfL83vz+cnBWLmXHoMTRuUI/bftOBP/3JJkJfcUWAB/QennoKzjjDkhy8+KICr4iIVCkju6C7dm5Dy9kfcPL7r+KB3nPHcPT1Pfns4zYUFMDYsVC3bkAH27rVsgsNG2YrV/3rX9C0aUA7FxGRTJWRLWCAX48ahgMc0Ki0mBO/WsSQIXDKKXD++QEeqKTEbi3685/hzTcVfEVEJCEZ2QIGYNOmXbIdvbIij6IieOONgHqHt2yxZvQ++8D779t4r4iISIIyNwDPmAGTJ8Ps2az9RR63dM/l4oshNzeAfZeUwGWXwcaNlntXwVdERPZS5gXg776zANmqlU2KOuMMBt9sDdb77gvoGH37WlP6r38NcDBZRERqk8wbA777bsuLu3EjAEuXWprBXr2gQ4cA9v/kk/D443DLLbaOpYiISDVkVgAuLLQFnrt33zEZ6q67LNHC3XcHsP+xY+HWW6FLF3j00QB2KCIitVVmBeC774b69S3qYnkQhg+3VIPZ2QHs//DD4aKL4JVX1PUsIiI1kjkB+NNPLTD26bMj2vbrB/vvb0O2NbJhgy228R//ASNHwr771ry8IiJSq2VOAH73XcjKsqgLTJliTw0cWMM8COvX29TpAQOCKaeIiAiZFIBvuw0KCqBFC7yHwYOtIfzHP9Zgn9u2wcUX237PPjuokoqIiGRAAPbeAiTsyC2Ynw/Tp0P//tCoUQ32fcstlrtw2DA4/fQaF1VERKRM9APwxInws5/ZghgxQ4bAQQfBNdfUYL9vv233L91xh2U3EhERCVC0F+IoLbWx2cMO29FCnTbNfp54ooat359+gl//GoYODaasIiIi5UQ7AI8cCfPmwUsvWeJ7LM9v69Zw7bU13HfXrna/r9IKiohIEkS3C3r6dEsDePjhcOmlgC3/PGUK3HknNG5czf2OH29LTJaWKviKiEjSRDMAz5wJZ55p6z6vWAGzZgE29tuqFVx3XTX3u3Yt9O5tk66Ki4Mrr4iISAXR7ILOz98ZIEtKID+fD3wu774LjzxiGQKrpU8fC+oTJuzo0hYREUmGaLaA8/IsQNata7/z8hgyBA44AK6/vpr7HDUKXn3VlrE85pggSysiIrKbaLaAc3Mt129+PuTl8aHLZeJEePDBaq4SuWUL3HADdO6sFa9ERCQlohmAwYJwbi4AQ86Fli0thlZL48bw2mu2k/r1gyujiIhIHNENwDGzZtmaGfffD02aVGMHGzda6sK8vKCLJiIiElc0x4DLGTLEVqC88cZqvLmoCI48Ep5+OvByiYiIVCXSAXj2bHjrLbj9dthvv2rs4IYbLNWgWr8iIpJikQzAb85bwckPTOG0HkXUa7ydtqes3PudTJxoM58HD4aOHQMvo4iISFUiF4DfnLeCAaM+pfDz+mxZ2op9cwoZOvET3py3IvGdFBdbs7ldO0tjKCIikmKRC8APv/MFW7aXsHlxNnUabqfpcV+xZXsJD7/zReI7mTsXliyBhx6Chg2TV1gREZE4IjcLeuX6LQBk5X1Ok87LqNOweJfnE3L88bB0KWRnJ6OIIiIiexS5FnB2lmVZcA7qZ23Z7fk9+vJL8B7atFGyBRERCU3kAnDfszvQuH7dXZ5rXL8ufc/usOc3FxbC0UfDY48lqXQiIiKJiVwXdNfObQAbC165fgvZWY3pe3aHHc9XqV8/qFcPevRIcilFRESqFrkADBaEEwq45U2fDiNH2sodGvsVEZGQRa4LulpKS+12ozZt4I47wi6NiIhINFvAe62wEL75pobJgkVERIJTOwJw+/ZQUFDNXIUiIiLBy/wu6DlzbOWr/faDOplfXRERiYbMjkgrVsBpp8Gdd4ZdEhERkV1kdgD+05+s9dunT9glERER2UXmBuB58+Cll2z28+GHh10aERGRXWRuAB40CJo3hwEDwi6JiIjIbmocgJ1zdZ1z85xz44IoUCA2bYJly+ye32bNwi6NiIjIboK4DekWYDHQNIB9BWO//WD+fBv/FRERSUM1agE75w4GzgeeD6Y4ASgshPXr7ZajBg3CLo2IiEilatoF/QRwJ1AaQFmC0bs3nHSSpRwUERFJU9UOwM65C4DV3vuP97Ddtc65Oc65OWvWrKnu4RIzdar9XHedcv2KiEhac76aLUXn3P3A5UAx0AgbAx7lve8Z7z05OTl+zpw51TreHnkPp55qXdBLl0KjRsk5joiISIKccx9773Mqe63aLWDv/QDv/cHe+7ZAd2BKVcE36SZNghkzYOBABV8REUl7mXMf8LvvwqGHQq9eYZdERERkjwIJwN77fO/9BUHsq9oeeshWv2rYMNRiiIiIJCL6LWDv4euv7XGLFuGWRUREJEHRD8BvvAFHHAEzZ4ZdEhERkYRFOwCXlsLdd0O7dvCrX4VdGhERkYQFsRRleF5/HRYuhFdfhXrRroqIiNQu0W0Bl5TA4MHQqRN06xZ2aURERPZKdJuNc+ZAQQEMHw5164ZdGhERkb0S3QB8wgkWgA85JOySiIiI7LXoBmCAww4LuwQiIiLVEt0xYBERkQhTABYREQmBArCIiEgIFIBFRERCoAAsIiISAgVgERGRECgAi4iIhEABWEREJAQKwCIiIiFQABYREQmBArCIiEgIFIBFRERCoAAsIiISAue9T93BnFsDLAtwly2B7wLcX5hUl/STKfUA1SVdZUpdMqUeEHxdDvPeH1DZCykNwEFzzs3x3ueEXY4gqC7pJ1PqAapLusqUumRKPSC1dVEXtIiISAgUgEVEREIQ9QD8bNgFCJDqkn4ypR6guqSrTKlLptQDUliXSI8Bi4iIRFXUW8AiIiKRFIkA7Jw7xzn3hXOuwDnXv5LXGzpGKQs5AAAFIElEQVTnRsRe/8g51zb1pdwz59whzrmpzrnFzrlFzrlbKtkmzzm3wTk3P/YzKIyyJsI595Vz7tNYOedU8rpzzv0tdl4+cc4dG0Y5q+Kc61Du/3q+c26jc+7WCtuk7Tlxzr3gnFvtnFtY7rkWzrlJzrklsd/N47z3ytg2S5xzV6au1JWLU5eHnXOfx/5+3nDOZcV5b5V/i6kWpy6DnXMryv0dnRfnvVVe71IpTj1GlKvDV865+XHem27npNLrb6ifF+99Wv8AdYGlQDugAbAAOKrCNjcAz8QedwdGhF3uOHU5CDg29ng/4MtK6pIHjAu7rAnW5yugZRWvnwdMABxwIvBR2GXeQ33qAt9i9+1F4pwApwLHAgvLPfcQ0D/2uD/wYCXvawEUxn43jz1unoZ1OQuoF3v8YGV1ib1W5d9imtRlMHDHHt63x+td2PWo8PqjwKCInJNKr79hfl6i0AI+Hijw3hd677cBrwFdKmzTBfhn7PFI4LfOOZfCMibEe7/Kez839ngTsBhoE26pkqoL8JI3HwJZzrmDwi5UFX4LLPXeB7lYTFJ5798D1lZ4uvzn4Z9A10reejYwyXu/1nu/DpgEnJO0giagsrp47yd674tj//wQODjlBauGOOclEYlc71KmqnrErrHdgOEpLVQ1VXH9De3zEoUA3Ab4pty/l7N70NqxTezDugHYPyWlq6ZYN3ln4KNKXs51zi1wzk1wznVKacH2jgcmOuc+ds5dW8nriZy7dNKd+BeTqJwTgFbe+1VgFx3gwEq2idq5AfgvrEelMnv6W0wXN8W601+I09UZpfNyClDkvV8S5/W0PScVrr+hfV6iEIAra8lWnLqdyDZpwznXBPg3cKv3fmOFl+diXaC/BJ4E3kx1+fbCyd77Y4FzgRudc6dWeD0y58U51wD4HfB6JS9H6ZwkKjLnBsA5NxAoBl6Js8me/hbTwdPAEcAxwCqs+7aiKJ2XHlTd+k3Lc7KH62/ct1XyXI3PSxQC8HLgkHL/PhhYGW8b51w9oBnV6/5JOudcfezkv+K9H1Xxde/9Ru/9D7HHbwH1nXMtU1zMhHjvV8Z+rwbewLrPykvk3KWLc4G53vuiii9E6ZzEFJV19cd+r65km8icm9iElwuAy3xsQK6iBP4WQ+e9L/Lel3jvS4HnqLyMkTgvsevsxcCIeNuk4zmJc/0N7fMShQA8GzjSOXd4rJXSHRhTYZsxQNmstN8DU+J9UMMUGzMZBiz23j8WZ5vWZePXzrnjsXP0fepKmRjn3L7Ouf3KHmOTZRZW2GwMcIUzJwIbyrp60lDcb/NROSfllP88XAmMrmSbd4CznHPNY12hZ8WeSyvOuXOAfsDvvPeb42yTyN9i6CrMf7iIysuYyPUuHZwBfO69X17Zi+l4Tqq4/ob3eQl7ZloiP9hs2i+x2YEDY8/dg30oARphXYcFwCygXdhljlOPX2PdFp8A82M/5wHXA9fHtrkJWITNfvwQOCnscsepS7tYGRfEylt2XsrXxQFPxc7bp0BO2OWOU5d9sIDarNxzkTgn2JeGVcB27Ft6L2z+w2RgSex3i9i2OcDz5d77X7HPTAFwdZrWpQAbeyv7vJTd7ZANvFXV32Ia1uXl2OfgE+yif1DFusT+vdv1Lp3qEXv+H2Wfj3Lbpvs5iXf9De3zopWwREREQhCFLmgREZGMowAsIiISAgVgERGRECgAi4iIhEABWEREJAQKwCIiIiFQABYREQmBArCIiEgI/h9hu03SyqZCAgAAAABJRU5ErkJggg==\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:                Tue, 17 Dec 2019   Prob (F-statistic):           5.69e-38\n",
      "Time:                        23:42:54   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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd1zV1RvA8c+XCwhOzC2490CR1DSzzJGmZmZlZWlTM9umpg2z0rQwbZmlltWvzCyVLHMr5Z7gFnEDKqKCi3W59/z+OAwHyAXuQp/36+VLuOt7ArvPPec853kMpRRCCCGEcC4PVw9ACCGEuBlJABZCCCFcQAKwEEII4QISgIUQQggXkAAshBBCuIAEYCGEEMIFPJ15sfLly6uaNWs685JCCCGEy2zduvW0UqpCTvc5NQDXrFmTLVu2OPOSQgghhMsYhnE0t/tkCVoIIYRwAQnAQgghhAtIABZCCCFcwKl7wDkxm83ExMSQkpLi6qEUaT4+PgQEBODl5eXqoQghhLCBywNwTEwMpUqVombNmhiG4erhFElKKc6cOUNMTAy1atVy9XCEEELYwOVL0CkpKZQrV06CbyEYhkG5cuVkFUEIIYoQlwdgQIKvHcjPUAghiha3CMCuZjKZCAoKokmTJjRv3pxJkyZhtVqv+5wjR44wa9YsJ41QCCHEjcble8D5FRoeS8iSSI4nJlPVz5fhXRvQu4V/oV7T19eXiIgIAE6dOkW/fv04d+4c77//fq7PyQzA/fr1K9S1hRBC3JyK1Aw4NDyWUfN2EpuYjAJiE5MZNW8noeGxdrtGxYoVmTZtGl999RVKKY4cOUL79u0JDg4mODiYdevWATBy5EhWr15NUFAQkydPzvVxQgghRE6K1Aw4ZEkkyWbLFbclmy2ELIks9Cz4crVr18ZqtXLq1CkqVqzIsmXL8PHxISoqiscee4wtW7YwYcIEJk6cyN9//w1AUlJSjo8TQgghclKkAvDxxOR83V4YSilAn1N+6aWXiIiIwGQysX///hwfb+vjhBDCri5cgNOnQY4gFjlFKgBX9fMlNodgW9XP167XOXToECaTiYoVK/L+++9TqVIltm/fjtVqxcfHJ8fnTJ482abHCSGE3SxcCD17gmFAWhp4Fqm39JtekdoDHt61Ab5epitu8/UyMbxrA7tdIz4+nsGDB/PSSy9hGAbnzp2jSpUqeHh48L///Q+LRS+BlypVigsXLmQ9L7fHCSGE3SUkwJNP6uALMHo0ZKzaiaKjSH1cytzntXcWdHJyMkFBQZjNZjw9Penfvz9Dhw4FYMiQITz44IP8/vvv3H333ZQoUQKAZs2a4enpSfPmzXnqqadyfZwQQtiV1Qp33AGRkfDOO/pPsWKuHpUoAEM58VNTy5Yt1dWJSXv37qVRo0ZOG8ONTH6WQtzAEhKgTBnw8IC//gJ/fwgOhpQUCAuDBg1kH9gNGYaxVSnVMqf7itQStBBC3JTmzoWGDeGbb/T3992ngy9AUhLcey+EhrpufKJAJAALIYS7OnUK+vaFhx6CgAC99Hy1smWheHGIiXH++EShSAAWQgh39Ndf0Lgx/PknjB0LGzZAs2bXPs4wdHCWAFzkFKkkLCGEuGkULw5168J330GTJtd/bLVqEB3tnHEJu8lzBmwYho9hGJsMw9huGMZuwzDez7i9lmEYGw3DiDIM4zfDMLwdP1whhLhBKQU//QQffaS/79QJ1q/PO/iCzICLKFuWoFOBjkqp5kAQ0M0wjDbAx8BkpVQ9IAF41nHDFEKIG1hMjD7T++STsHQppKfr221tMzpyJCxY4LjxCYfIcwla6XNKFzO+9cr4o4COQGYroB+BMcBU+w/Rsc6cOUOnTp0AOHnyJCaTiQoVKgCwadMmvL1lYi+EcBCl4PvvYehQHXQ//xxefBFMpryfe7mGDR0zPuFQNu0BG4ZhArYCdYEpwEEgUSmV8TGNGMB+3RCcqFy5clmtCMeMGUPJkiUZNmzYFY9RSqGUwsNDctaEEHZ0+DC88AK0awczZkCdOgV7nfh4nax1zz1Qvbp9xygcxqaIopSyKKWCgACgNZBTtYccK3oYhjHIMIwthmFsiY+PL/hInezAgQM0bdqUwYMHExwcTHR0NH5+fln3z549m+eeew6AuLg4+vTpQ8uWLWndujUbNmxw1bCFEO7OaoUlS/TXtWvr7OYVKwoefAFOnICBA2HTJvuMUThFvrKglVKJhmGEAW0AP8MwPDNmwQHA8VyeMw2YBroS1vVe/7XXIGMyajdBQfDZZwV77p49e5g5cybffPMN6Zl7Mjl45ZVXGDFiBG3atOHIkSP07NmTXbt2FXDEQogb1oED8Oyz8N9/sHq1PtebWVAjB6HhsbaV3q1WTf8tmdBFSp4B2DCMCoA5I/j6Ap3RCVirgIeA2cCTwJ+OHKgr1KlTh1atWuX5uOXLlxMZGZn1fUJCAsnJyfj62rdLkxCiiLJY9P7uO++At7fe923X7rpPCQ2PZdS8nVk90GMTkxk1byfAtUHYz0+KcRRBtsyAqwA/ZuwDewBzlFJ/G4axB5htGMZYIBz4rrCDKehM1VEub6jg4eHB5XWzU1JSsr5WSknClhAid716wT//6Eznb77RdZzzELIkMiv4Zko2WwhZEnltADYMOQtcBOW5B6yU2qGUaqGUaqaUaqqU+iDj9kNKqdZKqbpKqYeVUqmOH67reHh4ULZsWaKiorBarcyfPz/rvs6dOzNlypSs7yPsvY4uhCh60tP1fi/AgAHw88/6qJANwRfgeA69z693u5wFLnokrTcfPv74Y7p160anTp0ICAjIun3KlCmsXbuWZs2a0bhxY6ZPn+7CUQohXC4iAlq3hqkZJzMfeQQef9z2c71AVb+ct7Byu53vv5ezwEWMtCO8gcjPUggXS02FceNg/HgoVw6+/Rbuv79AL3X1HjCAr5eJ8X0CC90DXTiPtCMUQghH27JFZzR/+CE89hjs3l3g4As60Wp8n0D8/XwxAH8/3+sH33374L334PTpAl9TOJc0YxBCCHs4f17/+ftv6NHDLi/Zu4W/7bPdQ4fggw90b+Dy5e1yfeFYMgMWQoiCWr06+/hGx476nK+dgm++ZealSCZ0kSEBWAgh8uvCBXjpJbjzTvj6a0jOyEwuVsx1Y8osxiGZ0EWGBGAhhMiPpUuhaVMdeF99FcLDwR2K7kgxjiJH9oCFEMJWJ07AfffpGs5r1sDtt7t6RLpncFgYdOigZ8HHc6wKLNyQzIABk8lEUFAQTZs25eGHHyYpKanArxUWFkbPnj0BWLBgARMmTMj1sYmJiXz99df5vsaYMWOYOHFigccohMinzCYHVarA4sV61usGwTdpyWosd9yF9a23UZ06wZdfwqxZrh6WsJEEYMDX15eIiAh27dqFt7c333zzzRX3K6WwZla0yYdevXoxcuTIXO8vaAAWQjjJqVO6iMZtt+mlZ4C77wYfH9eOC1j31TbO93gEk9WMBwrS0vRRqHwU+xCXSU7W57adWBujaAbg9ev1Qff16+3+0u3bt+fAgQMcOXKERo0aMWTIkKx2hEuXLqVt27YEBwfz8MMPc/HiRQAWL15Mw4YNueOOO5g3b17Wa/3www+89NJLgG5Z+MADD9C8eXOaN2/OunXrGDlyJAcPHiQoKIjhw4cDEBISQqtWrWjWrBnvvfde1muNGzeOBg0a0Llz5ysaPwghHEAp+OUXaNwYQkNh7FgdeN3A6TgLC5uMoPXLrSmmUkn38MaMCeXpDSVL6m5L1+neJnIRGakT69ascdol3S8Ad+hw7Z/MWWJSErRooVt4vfWW/rtFC/jhB33/6dPXPjcf0tPTWbRoEYGBgQBERkYyYMAAwsPDKVGiBGPHjmX58uVs27aNli1bMmnSJFJSUhg4cCB//fUXq1ev5uTJkzm+9iuvvMJdd93F9u3b2bZtG02aNGHChAnUqVOHiIgIQkJCWLp0KVFRUWzatImIiAi2bt3Kf//9x9atW5k9ezbh4eHMmzePzZs35/OHKoTIlyefhCeegLp19XLz22+Dl5dLh6QU/PorNG7qQdreA0S0eBrfmAPs+yaM0XzI+rErdBb2999DLu9D4irnz+sa3aB71+7fD+3bO+3yRS8J69y57ALnVqv+vpCSk5MJCgoC9Az42Wef5fjx49SoUYM2bdoAsGHDBvbs2UO7jBZiaWlptG3bln379lGrVi3q1asHwBNPPMG0adOuucbKlSv56aefAL3nXKZMGRISEq54zNKlS1m6dCktWrQA4OLFi0RFRXHhwgUeeOABihcvDuilbSGEnWW+r3h46GIWwcHw8stgMrl2XOvXc+63ReyevZPRcZ9Qq3U96i6dQ2AL/fZdrk09ejGUi/tqQ5+Ms8AxMdnngkXOFi+GQYMgNlZvMdSrB7VqOXUI7heAw8Jyv694cb0s1KmT3u/w9tbft22r7y9f/vrPz0XmHvDVLm9HqJSiS5cu/Prrr1c8JiIiAsNOey5KKUaNGsXzzz9/xe2fffaZ3a4hhMjBwYPw3HPw4IN6GfKxx1w9IgCsa9ej7upAaUsabYFfOjXi1iUfYTJlv3VXrFuaSmxkVeReqPagvjE6GjImD+IqZ8/C66/DTz9Bo0awdq0Ovi7gfkvQeWnbFlas0PVWV6zIDr4O1qZNG9auXcuBAwcASEpKYv/+/TRs2JDDhw9z8OBBgGsCdKZOnToxNaMzisVi4fz585QqVYoLFy5kPaZr1658//33WXvLsbGxnDp1ijvvvJP58+eTnJzMhQsX+Ouvvxz5nyrEzcNigUmTIDAQtm2DUqVcPaIsUf+d4EiXgZgsaRgAHiZadyp1zYTc5OvNaY9KmE5eNuuVs8A5S0/XH0xmzYJ339XbCy78oOJ+M2BbtG3rtMCbqUKFCvzwww889thjpKbq1sdjx46lfv36TJs2jR49elC+fHnuuOMOdu3adc3zP//8cwYNGsR3332HyWRi6tSptG3blnbt2tG0aVPuvfdeQkJC2Lt3L20z/ttKlizJzz//THBwMI888ghBQUHUqFGD9k7coxDihrVnDzzzDGzcCD176taBbrBsm5YGH38MvmM+42VrFFaTFwZWDG/vXPNazhQPoMSZaF2Mo1w5uHTJuYN2d/HxeoXU0xMmTIA6daB5c1ePStoR3kjkZylEPqxapY8Yff45PPqoWxzf2T73AB+NPM+cA8EM6HORT984TnnjTHahjVwmHpuq9aFsXCT10nbrbC03+G9xC0rpJN2hQ2HiRJ0h7mTXa0dYNGfAQghREJs2wYYN8Mor+ljR4cNwWa6Hq1w6l87K+ybTefVoRno14/HQDfS6vyRQXz8gjxW/U3Vu5/BxH+oqJF8k0+HDOslq+XKd2eyGK4dFbw9YCCHyKykJ3nhDB7JJk7KXaN0g+O4c8SMXytXgvtUj2F+zK3W2z6PX/fkLoocfHMaj1lnEx6OPIblJEpnL/Pijrte9YYM+xhoWBvXru3pU15AALIS4sa1apZOsJk2CgQNh+3a3CLxnzsAPbb6hachTVLIcx+rpTfNfRlC6kY39fy+T2QgpOho985sz5+YuxlG+PNx1F+zeDS+8oI+WuSG3GJUz96FvVPIzFCIHcXH6TK+Hhw7E33wDZcq4dEhKwfzvztK4MRzYdAaFgQFYLens/rVgJxzqJ28nmgCSFyzTiWRW681VjCMtTZ+MGTdOf9+jByxcCNWru3ZceXB5APbx8eHMmTMSQApBKcWZM2fwcYP6tEK4hcxqcZUq6Tfi7dvzXRnPEWJ3J7Kk5vO0e64hNfyOs+e+8qR6epFueGA2eTL2YkVCw2Pz/boV65UhgFhS9x/NzuSOjrbz6N3U5s3QsiWMHg379mXXci4Ce+EuT8IKCAggJiaG+Ph4Vw+lSPPx8SHADY5QCOFSp07pBKvfftOVjrp21YV7XMxqhaVDQmk+bQhdVBzhHYZS4s7tbEsN4PHS42hzbCcbqgeyrVJ9ji2JpHeL/C1DlwusihUDy9EYqNZa33ijnwVOStJBd/JkqFwZ/vwTiliVQJcHYC8vL2o5ufyXEOIGo5Su6fvaa3DxInzwgXs0T1i/nvjZy4n+cRXdzq3iYMnmnPj5L1refytHRi4EYJt/I7b5Zx8fPJ6YnO/LGMW8iTdVwvNkjN4QrlPHqV19XCIqCr74Qlcw++STPLcWQsNjCVkSyfHEZKr6+TK8a4N8f9CxN5cHYCGEKLQBA3QAbtsWZszQXYzy4Og35PTV61GdOlHWnEZZYG+HwTRc8gWGt27qUNXPl9gcgm1VP98CXe9M8WqUSIiGsmUho2LfDScxUc90n3xSF9KIioIaNfJ8Wmh4LKPm7STZbAEgNjGZUfN2Arg0CLt8D1gIIQrEas1uoNCjh54NrV5tc/AdNW8nsYnJKLLfkAuy/5qTnaEHie78NIY5FU8smEzQ6J7qWcEXYHjXBvh6XVlX0tfLxPCuDQp2zToP8J/V/c662s2ff+rf7bPPZn/AsCH4AoQsicwKvpmSzRZClri2tasEYCFE0bNvH9x5J0yZor9/9NF8dS5y1Bty0vl0/rprInUeCKRiWjSGpyeYTDmWkezdwp/xfQLx9/PFAPz9fBnfJ7DAM7Lt3Ucx6tI7WCzAm2/eOGeB4+Kgb1/o3RsqVNBne+vWzddL5LasX5DlfnuSJWghRNFhNuv9vg8+0Gd5y5Ur0Ms44g1547Tt+L78HPelbWF7jV7UWvQ1psRj1y0j2buFv92WQKtVAyzpnIj1ICAuTq8GFHVms26WcPw4jB0LI0YUqC+zvZf77UUCsBCiaAgPh6eegh079Izoiy/0MaMCsMcbcmh4LP9Mm0edHdv478wDdIoMY4DpGLvHzKH56IcyjsH4O7RxzOX72D037ieVYexcvYOAatXgxAldjMOzCL7Nx8ZC1ao62H7+ua5i1bBhgV9ueNcGV+wBQ+GW++1FlqCFEEXDuXO6l+uff+pjRgUMvlD4/dfQ8Fhmfz6HL6e9wRvrfuL3yEf5t3EwqxeuoMl7DzvlDOrV+9jRxb0wYWXzyr1FtxiHxaKPFdWvDzNn6tt69SpU8AX7L/fbSxH8aCSEuGmsXKn79A4bppdxDxyAYsUK/bKZb7wFzYKe/NNOQmb/j2JWMwDKSKVd+WVMDK/KA10LPTybXL2PHVepNAAHNkfCgy30jdHRbtFi0Sa7dukjRRs3Qvfu0LmzXV/ensv99iIBWAjhfhISYPhw+O47Pft56SXw8bFL8M1UkDdkqxWWvbyAX74eQhWOYzZMGCjSTZ5sqB7o1KSeq6916pYyWDHwO31KnwPu0MFtayBf48svdbOMMmXgl190AlkRqGRVWBKAhRDuZe5cHXDj43XSzZgxOvi62P79MPPBvxm/6352eTdmUPe38Sl5KbuKlX8j/J2Y1HP1PrbF05OTRiWqJcdDgwa69rW7y+xdXLeu3tefPFlnOrtA5n567JkU/Mv5OKVQhwRgIYT7iImBfv30ec+FCyE42NUjwpymmPHOEV7/ohbFi91LjwHTOD2kM4f/1kvAmVWsnJ3Uk5lY5HXhHLcknePILf5MLTuQ88Ub8ITTRlFAFy7AW2/pLPYxY3TDjHvvddlwQsNjGfHzPgLmX2TAuY3svPcWRl1KAxxbqEMCsBDCtaxW3TT9nnv0fmVYmC6uX4DjJvZ28OM/8Bo9ksfS4tnY8xDjp5WjSpWB+k5vb5eWNuwdVJWAf+bRcuwrANw5dikb7nyTXRtK8Dno4iQVKsAPPzhtTDb55x8YPFh/2Hr9dVePBqXgzU8SqDDfk79TH8ObVMy/efL4o+MIWeLt2N+pUsppf2699VYlhBBZoqKU6tBBKVBq1SpXjybLpfPpanWLl5UVlBWUxeSl1Jo1rh5WtqgopTp31j+3UqWU+vhjpcxm9f67ZlWNYyo1VSnVqZNSbdq4eqTZTp1Sql8/PeZGjZRat87VI1LHjinVs6ce0pTizymrjsfKbHioj+8coGq++XehrwFsUbnExCKyQy+EuKGkp+uCGoGB+nzv9Om6gXqG0PBY2k1YSa2RC2k3YaXdSkTa4t+/LxBVoS13hH8JgAF4YIX//nPaGK7LbNaNJjZt0pXAEhL0XrmnJ932fcYxqnN833ldmcOdWhLGxkJoKLz3nv6dO/B8dF6sVvj6a2jSBFavSKNmj/2s6lUHi+GR1RpyQ/VAhxfqkCVoIYTz9egBS5fCAw/AV1/pogsZXFU4P+GsYvgIg+++K8UvpW+l1GM9qf3bBN3sPYdSkk63ZQu0aKGX5n/6SSdaVa0Kycl6Cb9xY3zq6iNHpyNiqBkQ4PpiHEeP6qD76qsQFATHjhW4epm97NuXcdpprZlv6n3KgLQZLBnxF8OXBdH38Y+zkur21mzKeAfv6csMWAjhHElJutACwKBB8McfMG/eFcEXXFM4P+zD1cRWasG6mZGMGAEPnJxK7ZmjYcUK+PBD/berZmwJCfrn1apVdnGKu+/O/rmlpOhiFQsWUKaJDsDnd0frGbDVqoOws1ksuoJVkybwzju6lCS4NPimpelqls2bg/fOrcRVb82zUaPwahlEz8YVGN8nkLgmwUxt25e4JsFOKdQhM2AhhOOtWKGDyKuvwiuvwIMP5vpQZxbOPz17OWeHvEOHhI3Eetdk/vSzNHjqsge0beu6wKsUzJ6texyfOaOLkeTUYMHPD4oXh+hoKnTrA0DKgRjo0xz693fyoLmyoMa998LUqdd8yHK2TZv0kMrsXE14hbdpdHoNRonK+gPgAw8A0Lu881sTSgAWQjjO2bM6cMycqc96NmuW51OcUThfKdj48ERazx1BORQWD08qLZiBf1fX7Ute48UXdfBq1QqWLNFLuDkxDD3bjYmheN2qWDFQx6LhtmfhttucO+akJD07B7coqHHpErz7rp6M97hlPfO9u2KKT9Zds374QWfeu5AsQQshHGPRIn2e96efdHu8HTts2ke1d5/cqx04AB07gnluKAYKAzAZCs9tm+zy+oViNus9XchuOLF+fe7BN1NAgD7W4+VFSNXJhBXLqIeplH5NRwsP19cqXlzX6d67V5/ndmHwXboUmjaFmZMTWNVoCL8/Ph+TJS37AVu3umxsmSQACyEco3hxHRg2b4YJE8DXthmsowrnp5sVC/v+yOAmqwkPh8ShY/WYTCaHJFnlO5N7wwa49VY9ZQM9Hlt7HF+W8by6xassv5Qxk69SRWdIO8q5c/DCC7pgys8/69s6doTy5R13zTycOQNPPgldu8J9aXOJu6Uxd+6bRrFSxfTv2UG/74KQJWghhH1YrfDNN7qB+vvv62NFmzcXaBZk78L5u/8+zPl+z9PjwjKKVxtAgw3tqVq1Azy04rr9egsqX5nc587BqFH6Z+fvD3femf8LvvWWrqUMNC5/iqQ1MUCw3h+OiSnMf0ruQkP1MvnJk7qgRsZeqiuEhsfyyeJIotb7kbiyKZVT4tnV4GWaRM7XmeMzMqqqde/ukN93QUkAFkIU3t69MHAgrF2rpx4Wi55puLKg/vr1mJesZONfp2ixbQbK8CD8uSnc/e3g7LU/ByVZXS+T+4oAHBaml2rj4nRy2ocfQqlS+b9gvXpZXz68fxzvnJtJUtJ5igcEOOYs8Kuv6uXxZs10IG7Vyv7XsFFoeCzDfthPwNxknovdyNpbTvNMzek0PPAPfPwxDB2afQzLlUl1OZAALIQouLQ0vbw8bhyULKkTWwYMcH0nm/XrsdzdCY/UVO7ASlTF26mwfDYtAqs55fI2Z3JXrgw1asCCBbr8ZkGdOqWbWHTvjkeNapRef4ED+85Tt1o1WLas4K97OaX0mWIvL53dXLmyTrBzQsnQzEYJV5f9tFph6JiLVP7Hk0XpD+FFGubzngzs/A6Lut3HbyMGOHxshSEBWAhRcEePwkcf6WNFn30GFSu6ekQkxqWy6omfuC81DU+sWPAg7eH2+OUz+Ob2pm+L3DK5A0p765/T9u06M7xhQ1i3rvAfWOLiYMgQmDMH33r6LHB8eAx17VWMY/9+fYysfXs9S+/WTf+xUWF+lrkt58cc9uS3yZU4tqYO3/i8gE96CgaAJZ1mJw8wtZbrG3nkRZKwhBD5c/48zJihv65XTy8/z5rlFsH334/Wcso/iLsPzcJseJJueJDm6cnYS5XyVc4y800/NjEZRfabvq2vkVMm962nD7Pg52F6v/TUKV1AA+yzWhCggy4xMZRpqj9onN8TrROiRo7UKxUFkZamVzeaNYOICKhdO98vUdif5dXL+cpicOLfWrzStwIe28PZ7teMe1JWYXVyGUl7kBmwEMJ2Cxbomdbx4/qMaWAg1Krl6lERF3We8G6j6Hboa46aqjOk87ukVCa7V2+l+hy7ev/1Omzew81F5mNClkSScCqBdzfP5tF18zAqVNDHdB5+2L7L9H5+UKIEREdTvocuxpF6IAbufjb7XG5+RUTo7YSdO+Ghh/Seb5Uq+X6Zwv4sL1+2Tz1RhjOLmmGOL023ar/zz/HHSC1bjlceeocYnzJOLSNpDxKAhRB5O3lSJwn9/rs+XDl3rg6+LqYUzPo8ng5Dg7lHxbL59ld5smV7knx9ALJ69UL+KmnZoxpXViZ3fDw0HqyT1CZM0MHS3gwj6yxwsVpVebH0/yjt245eSkFiok6IK106/6954QL8+acudVlAtvwsr7dEXdXPl+hTadRfdJFW+7azzSeR8D6VuRRcBiPtbXxee42OR5IIWRLJVP9GVPXzZbyTW0MWlARgIcT1WSz6aMyxY7qY7vDh+hyli8X8+h+L31nLd4c6QLXHaT+xN636tqHshJUkFbKSVqGrcZ08qTsVjRmje/JGRsItt9h8/QLJPAvs5cWm+k9QLhFITNDXnTgx65jSdf3zj+76NGGCLpocFVXoRg55/SzzOrJ1j18zlo/byYKLffAhGSMFnvMZy33dn4IWupJV77Jli0TAvZrsAQshcnboUPZxoilTdCWrt992efBNNyvWdh5N1X4deOrQO/zn1YnHfr2f6n3bAPappFXg17BaYdo0aNQIQkJg2zZ9u6ODL+gM9EWLAGhfejuV94VlL03ndRY4Lk6XjezRA/76S898wS5dlPL6Wea2RP3R/NUNaYkAACAASURBVIM8/TSMeaE8r6np+JKMB6AwGFYqoUgG3KtJABZCXCktTWc2N26sm6YCdOkC9eu7dlzAnn+OsKn8vbRb8SEGCk+seFnT8PgvLOsx9qikVaDX2LtXFx95/nldOnLHDueej/X3z1refvLoB7wZ/eIVS9M5UkpnYzdqpBsTjBmjPzQU5CxyLvL6WV69RK0UXNpbha2f3sYfPyWxs8GDPHjpVwzDAJMJD18fGj2ZezOPokSWoIUQ2TZu1G1jdu3SiTcPPeTqEQG6PPKy+7+i47KR1DAgssfr1F/5Ta69eu1RSStfr6GUrt18/Dh8/z089ZTzz0Lv2QP/+x+88QbpVapR6+Ayzp2DMpeVqbxGfLzOyg4MzJ65O8D1fpaXL1GnXyjG2aVNST5QmRL+51nzdymavm+Fp8ZDu3awZo3bVLGyBwnAQggtJEQ3TfD3L3TijT39+6/OX3o+6jA1/e+k+j/fsNdi4tNqtai7ZwsHGreku091erticGvW6FKHJUro7j+VK7vuONbRo3rvtlcvTDUCKL3mAnv2nqdMQMCVxTjMZn1srH9/Pdb166FBA/BwzYLo8K4NGDl3J/Fb/ElY1ZA6lkN8Vf4xrFPG0zyojZ6ZZ36Yad/eJWN0FAnAQtzsMvd5b7tN1/YdNy7/GbP5YGtRhosL/yV60Af8drwP6bVeJGjRBJp19SQ04rhO2ilTm6Vt9bnU1bnVWXaUhATd5GDGDPjgA91AwYZWiw512Vlg3/r6LPCZiGh4/PHsGePGjfrTzM6d+kjRPfc4bNZrq8Yl/DEtKk+jbZt416cbHSz/4ZHii5d3gn6Aq6uqOZAEYCFuVidO6Jq+/v4webLOdLahEYAjqhrBlcEzfMg0mk8dTCMUX3r8R9qMYHw76iBS2HOlhaIUzJmjf26nT+uMcFuyi50hMwBHR+PXVPcBvrA3BgZ31UlVr74KX34JVavqFQ4X98I1m3Vy9vvvw1Mes5lqPI2RovRM/JefdbnLG1yeaw6GYVQzDGOVYRh7DcPYbRjGqxm3jzEMI9YwjIiMP90dP1whRKFZrTB9up75LFgAlSrZ/FR7VzWC7OAJEHfgAovqvUzzqc9joADdq9d3Y1jW4+1xRrfARo+GRx/NbrP4ySe67aI7uCzjudzdzehorCLc+zZITYV33tGFNIYM0XvFLt5e2LJF56e99Rb07AlfNPwaQ+nfN4YBu3e7dHzOYsuifzrwhlKqEdAGeNEwjMYZ901WSgVl/PnHYaMUQtjHgQO6MtKgQXrvcscOXarwMtfrY5tXAM1LbkEyNiGZmTNhQvNf6XpgCvuaPpRrr97czuI6rPSgxaJbBgI88QRMmqR797Zo4ZjrFVRmxnNcHJ5lSxHl34GoeD9dB9pi0TWnv/rKodsLeUlK0osGt90GtaP/ZUXINv74A7w/GefQ3szuKs8laKXUCeBExtcXDMPYCxT9A1hC3IysVl1cYcYMeOaZa/bX8loiLuzsM6eiDKVOpBKwvBjPHIe77niWY0Nb0fiBFjo5KIfercO7NrhijJD/c742i4jQWeE1a8Iff+hkpQZuXOJw2zYdyIBHS/yFR4Q3lOiqA6+LrVypt59PHzrHqoYjuHPfNFjdC4b9CZ07wwrH9GZ2Z4bKnPbb8mDDqAn8BzQFhgJPAeeBLehZckIOzxkEDAKoXr36rUePHi3smIUQ+bFmje7ZOnGi/j41FYoVy/Gh7SaszLFqkb+fL2tHdszz/rxkBvhGR3Zx29GdqGMlefboHCyGJws+P8JzLxazKRm3MPvQNklK0puTn34K5crp5du+fYtUQtCBcq05kVKW9peWuHQcCQm6a+Ge79fzZsmv6W5ajPeFs7pP7/vvu88S/qVLupjJkCF2/T0bhrFVKZVjr0mbA7BhGCWBf4FxSql5hmFUAk4DCvgQqKKUeuZ6r9GyZUu1ZcuWfA1eCFFAiYl6efnbb3XP2c2bdVnE66g1ciE5vSMYwOEJPa6ZIYOefean0MW/PyzgtucextuShgdwsHgjSvw5h8qdm9r+3+ZIERG6veKhQ3r2+8knULasq0dlm2XL4KefYOZMdjbqi+fBfTS07HHZ54a5c3Vifd349YTRAU9rmg5umSsw7iQiQm9Mh4XpM8d2cr0AbNPBL8MwvIC5wC9KqXkASqk4pZRFKWUFpgOt7TVgIUQhKKWXSxs10slWQ4fqpJY8gi/kvb9a2CpTKSlgfLuGYhnBVxke1H7rcfcJvqCzwitXhlWr9M+vqARf0B8afv4ZTpwgvWo1/FUMZ844fxjHj0OfPvDQQ4rg8sf47YUwPI2MD20eHrr0pTs4dSq7tWZQEBw8aNfgmxdbsqAN4Dtgr1Jq0mW3X96X6gFgl/2HJ4TIt/Pn4YUX9DnPTZv0MmqJEjY91ZYayL1b+LN2ZEcOT+jB2pEdbQ6+6/5JJCgI3tzwABbDC2UyYfgUw+iY99K1Qymlg1aPHjpZqUIFWLu2aCYCVdPnf4mJwbNGAKW5QOze8067fGaCfePGsG/hQQ7X7cLCs23wv7+lTq5ylyQrpeDHH/WH1Bdf1I1GAKpXd+owbJkBtwP6Ax2vOnL0iWEYOw3D2AHcDbzuyIEKIa7DYtFBxGKBMmV0+ahNm+DWW/P1Mvaoo3y186fT+LvV+zTtUR3/C/v4YGlbPNf+i/HhhzrxxpUJN4cOQbduuirU2bO4ZLpoT5cV4yjeQAfjs9tzKUNpZ1FR0KkTvDAonfHlJrLLI5CacZswRo/Wd6xYAe7wOz9wQNc2f+opHYDDw50eeLMopZz259Zbb1VCCDsLD1eqVSulQKk//nD1aK7w3yfr1V7PJkqB2tqwn7p4+JSrh6SZzUpNnKiUr69SpUop9dVXSqWnu3pUhXf2rP538Omn6mRkoqrFQTXlc7NDL/nHphjV+7Y5apTHWNXV6291pEqQHkOvXkpFRzv02vmWnKxUxYpKlS6t1NSpSlksDr8ksEXlEhOlEpYQRdXVmbq//qo33tzA2YXriHvqTdqdXkOcVzUiP/2b4KE9XD2sbOnpuvlAly661WLmzLGo8/PThVVSU6lQtwwxXmU4dtxxl/t01ikWvLyHRWefxJs0zB6ebCoeSNyEqbQe8bz7ZI3v2aNnuz4+OtO5WTO91+9i0o5QiKLq4Yd1hu5TT+lWeI8+6vI3PKVg0sAl+PTsTIPTa7EaJvZP+oQG7hB8k5L0EuilS/qNeN06fTzrRgm+oH//J0/CqFF4eMDoUpO5ZeMiu18mKUn37RjWvzxPnJ+NL8l4YsHLYmajf2NeV/Vd/m8RgIsXdRJiYKBulgG6xKUbBF+QACxE0XLqlH5TAd0AICxMZ3E6o+F7Ho5uPc3Cio9TasYfeJOGBwpQbJ23xOZSlQ6zbBk0bapLSS5cqG8rV849goQDDb4YQpO9f9j1NVet0hPIbz45x3S/JxmY/j0KSDcMzCZPNlQPdE5Z0LwsXqx/55Mn68pv993n6hFdQwKwEEWBUvDdd9CwoQ4iAG3a6AbwLmZJV/zzxCyKt2xE19NzSK2RhtnTRLrhgdnkyRr/JjaXqrS706dhwADdeMDLS39g6dvXNWNxlmnTsvo4nysVQKlzMXZ52cREXcmqY0fodGkBceWb8HTCLKa36s3jj4xjUvv+PP7oOLb5N3JcWVBbDR+uZ7q+vrB6NUydqpMT3YzsAQvh7vbtg+efh//+0/1QBw509Yiy7Ft6jDN9X6D7uX+I9GtNv3ueI6pWVVbGBtLm2E42VA9km38jDFfNiJ5/XjecePtt3ZDAx8c143Cm6GiYPx/S00kuX40KZ/ZitRau3e/8+fq0Tlycjm3jEhbhtakcYZO/ZVKkiWSzhfU1mwMOLAuaF6X0KQBPT33MqXhx3e0hl6pv7kACsBDu7JdfdMWg4sX1ActnnnFZ4/TLpf27nn/fD2Nn2GmeV2Fs6z+ZFt+/TNLEfyExmW3+jdjmn91n1qkzoiNHdKCtXBk+/hjGjNF7gDeLatX0gdwTJ7BWCaB65DLi4vSx8Pw6cQJeegmOz1vHj2UmUe+te6n54bNwaSJ4e9PBy4vxji4LaovDh/WHrXbt4L339JnuHm6Qd5AHCcBCuKP0dP1JvnVreOQRCAnJV9tAR4oc8yu1Pniau1U6d5q8SZ05i+D+vQEnN0q4msWi+92+/Tb07q0/vNSt6/jrupvLzgJ71qpG6bAL7N97nipVbO+CpBR8/72u4dzt4h+sNR7B45wVxodC98ZXnOPt3cLf+QE3U3o6fP653pYxmbKW3osK13+UFkJkO3tWLzFn7lPWq6dr+7pB8D1/Oo2/W39Anff746VS8cRCMdIoHbM36zGOKORhkx07dFB4/XW9/Dh+vGOv584yq2FFR2MZ9AI+JHPkrO3B98ABXTdj0HMWxt4yiV/oh4eyZj8gLMy+4y2oXbt0HsSwYXrAe/boZKsiRGbAQrgDpfQ53tdf19WYhg7VMzqTKe/nOsHaSRsp9+Zz9Ezfxb5qnWkQvwbM5hzLCjp9RjR3rj6CVbYszJrlFsexXCogQJ959fDAv34JUtHbwnlJT9cJw6NH61/rqgE/cOdPb+hl3a1bc/19u0xamj5yNWeOnvkWxd95bhU6HPFHKmEJkYOYGKXuuUdXD2rdWqmICFePKEtcnFLD7t2lLBjqhGeA2hvyl75j3TqlPvpI/+0qKSn67/h4pV58UanTp103FjdlPX9BTfF8RX3dZ9l1HxcerlRwsFLFSFavdNypYmOVUmlpSi1YoJTV6h6/b6WUWrFCqdGjs79PTXXdWGzEdSphSQAWwtVOn1aqdm2lvvzSbcohWteuUzu7D1f3lFqnvLyU+rPXDJUaf87Vw9ISEpQaNEip225zm5+X20pLUxYMNafR6BzvTkpSauRIpdp5rFM/+gxSF8tXV9YqVZS6dMnJA83D2bNKPfOMDll16yp1/ryrR2Sz6wVgWYIWwhU2bICvv9aZLuXKQWSkTrpyAyd/WEyFp3vSBAsLjC85/vNKavV71tXD0i4/DzN0qF43dZNlercybBgcPQq//06Cd2WKnb72LPC//+p0g+pRy/nXuBdTSjqkGjBpks66dweZrTVfflmf6X7zTZ3l7Ovic8Z2IklYQjjTuXMwZAjcfrsuKXTkiL7dDYKvJV2x6MnZlHz6ITywYADeHmZqHQ1z9dB0ctqDD+pa1xUr6k5PISFufcbTpc6cgfXrAThXOoAy57M3gRMT9YmdDh2gXOpxFpV4EJNK13d6eECyG1SxynTmDDz7rC4duWULTJhwwwRfkAAshHMoBb//ritZffstvPKKztp0k2Myu3ZYWV2xD/f+9BhnSlRHFSsGJhOGuyTdFC+uWweOHw+bN+e7zeJNJyBAH+JNTye+ZGUqpMZQc/g/NOy/g9r1LPww3cywYbBiTxW8enTTH2TcpVev1apXOpSC8uV1AZqNGyEoyLXjcgAJwEI4Q3q6Ti+tUkW/mXz2GZQq5epRkZqieO89CG7pweaUZmx9YhLVE3fisWqV63u3RkXBE0/o2tc+PnoGNHKkLikpri+jGMeSZeHsoAzepBH3Rysifw7kgdQZxN9Si5AhhylewoDfftOrMa7+fYPeiunQQa90ZNbsDgpyixUih8htc9gRfyQJS9xU0tJ0YlVmwsjRo7oPrZvYNmuv2uJ7h7qbFeqJJ3QysVswm5WaMEEpHx+lypRRau1aV4+o6Fm4UClQAwd/oSo+vF61YZ362HhDbSoVqBSoiBpNlYqKcvUos6WlKTV2rFLe3kqVLavUzJk6+/oGgCRhCeFkGzfqjbbt2/Xs7bnnoHp1V48KgEuL/+PYwA9pEhNGkkcpJr5zjuAPXT2qDNu26Z9VeDg88AB89RVUrerqURU9detCly7EJaVzu9dafjNG4aXMcAGmterNhLuf4ZCbbH8A+ne9cKEuQPPFF25ReMYZZAlaCHs6d04Xz23bVmdtzpunk0jcxI7XvsP33g40ilmOp6Hwnv0/gj98wNXDyjZihN67nDtX/+wk+BZM/fqwdCmnGzWnw5EteCkzBmAxPEj0LU2VsiVcPULdlzk1VX/96qu6N/Nvv900wRckAAthXz16wJQp+tjEnj36k70bVOiJj4fHH4cjn/+JgQJ0wmvxAztcPDL0/uPx4/rrmTP1z61PH9eO6QYxvGsDNtRvTaqnd1Z7yPDaQa7pVnS5Zct0g4zMkqFdusD997t2TC4gS9BC2NPChRAbC40bu3okgE4kXTlyKT99fZHfU/vQ8ek34dflYE5zfcZrYqI+r/rdd/po1pQp2XWMReF17Urv8uVh2Ce84uNJ3T1bONC4JY8M6uO65glnz+rz2z/+qGfpHTu6ZhxuwtB7xM7RsmVLtWXLFqddT4ibWXTEGfb1eIMux39kZ8m2eKxfS5Omhj4fGhamg6+rMl7nzdNL9adOwRtv6JaBN9D5TrfQpQtcuKCLvriDpUuhf399tvfNN+Hdd2+K/syGYWxVSrXM6T6ZAQthLwsX6je70aNdelTGsnodB177kgrbltCBC2y+5y2C572LqUTGUnjbtq49avLll/ocdFAQ/P03BAe7biw3soAAvdTrLipWhNq1dSBu3tzVo3ELEoCFsJe5c3VA+eADlw3h0C/rCejfkQYqFSsGZz79gVZDB2TdH+qq5ulK6eXHcuXgscd0Z52XX5YzvY5UrVpWMY7rnaN12L8JqxWmT9d7+p9/rj9wrVvnFjkR7kKSsISwl/BwPZtzwRtMarKVbwdt5fsBYXhklBU0TB5USI3NekxoeCyj5u0kNjEZBcQmJjNq3k5Cw2NzeVU7OXRIL4d266aDQfnyeh9Qgq9jBQToIHjiRK4Pcdi/if379f7u4MG6b29mtrME3ytIABbCHlJT9RtNixZOv3T47Eh2lu/A09PbUv72eph8vHMsIxmyJJJks+WK5yabLYQsiXTMwCwWXdi/aVNdPvK553TqtXCOFi3gmWeuG/Ts/m/CbNb1mps1g4gImDEDli+Xmt25kCVoIexh925IT+fdo178PHKhU5Z3LyaYCesRQuf1H5Bq+LLvtW95bdKDsME/xySr44k5F9nP7fZCiYnRR4k2b4b77tOdnwIC7H8dkbtWrfSf67D7v4n4ePjoI30c76uvdOlVkSv5OCqEHWz4bweXvH35r1Q1hy7vhobHMuiFLxnfdDCx5ZvSc/3b7KnbC9P+vTSb/DShEcdp928ytc41o92/yVdcv6pfzlnGud1eKOXL66zm2bPhzz8l+LqKUpCSkuvddvk3kZIC06bpa1WtCjt36nwICb55kgAshB28kVKdpq/9xlG/7Dcdey/vhobH8tuEOXz27QiG755BLeshPgoewrE5kylZt3Ke+3nDuzbA1+vK3rm+Xib7FWVYtw7uvTe7eUJYGDzyiOz7uVKFCjBqVK53F/rfxOrVOqP5+ef11wA1ahR0tDcdCcBC2MHxxGSU4XFNsLHX8q5SEPr6Ir6a8yHFVBqeWPAwrHiWvJgV5PPaz+vdwp/xfQLx9/PFAPz9fBnfJ7Dwy+QXL+pSgnfcoZfiDx/Wt0vgdb0KFSA6Ote7bfk3ERoeS7sJK6k1ciHtJqzUH+jOn4cXX4Q779T7vsuW6a9FvsgesBCFZbEQ+ttIZgTey1+N77riLnss78buPMvu7sP4IWYm0R7+lOQiSoHZ5MmG6oFZQd6W/bzeLfztuy+9dCkMGgTHjuk35I8+cos2iyJDQIDej7+O6/2byFxVyfxgl7mq0v7v0ZTbsRVef123MSzhBrWliyAJwEIUVlQUzY/somTgPVfcXNjlXasVFg8O5dYZg+moTjOlzvN8ft89ND19iDbHdrKheiDb/BvhnxHkq/r5EptDEHbIHi/oafmECXq5efVqaNfOMdcRBVetGixZUuCnX76qUjbpHJe8i5MMvNfqUb765ku47TY7DfTmJEvQQhRWeDgAXR7vZrfl3b17oX17WDg9lvMl/Tm1cAv+v7+LqWRxtvk34uu2fdnm3+iKIO/wPd5Mc+fqeteGAbNm6eMmEnzdU0CAPgdsNhfo6ccTk0Epeu35l+UzXuDF9XMAWHhLAwm+diAzYCEKa9s2KFaMjn060LGQxSXMYWs5MGQicyKD2Of3Hs/PfIG6jz+P4eVJ74zH5Fa1KPNvh1W6OnFC12+eN0/Xb544ESpXts9rC8fo1EmfvTabC1T4pLlxkZf++JTOBzcTUaU+CxvqD1oOW1W5yUgzBiEKq3Nn3dmnkP+2o8b9Ru13+mHCisUwce6v1dzSw4U1mzMpBT/8oKtXpaTA++/rr69T3lDcAObPx9x/AOmpZibe2Z+Zt96H1cOEr5fJPsl7N4nrNWOQJWghCqtOHejevcBPv5Ro5u9246n5zuN4YAXA5AG37Aiz0wALaeJEXVEpMBC2b4cRIyT4FhVK6Y5TiYn5f27t2njd0Y7V81ayuMtjKA+T/TLnBSAzYCFcaulS+HHACn6J60xk1Q7UP7sBw2zWvXpXrHBd1yKLRbeNq1hR/z1vHjz7rJSSLGoSE6FsWfj0U71qcT3p6fDZZ3DwIEyd6pzx3QRkBiyEo1gseT8mB2eik/ikyzK6doWtfp2I+HYjDWJXYaxcqY91uDL4ZmaAde+u35TLlYOBAyX4FkVlyugjQtc5Cwzo6lVt28Lw4XD8eIGTtkT+yP9RQhTGBx9A9eo2v2Gpdes5ePdzpNaoz6vLezL+lRNEREDQoNa64EEuZSSdwmyGceN027jISHjtNTCZ8n6ecF+GoY8i5XYWODUV3ntPd/E6ehR++w1CQ6VTlZPIRo4QhREeDiVL2vSGFf/zEm4Z0IM6yoIVg+Mjv2DkeF26MreCB4Bz9tuOHYNevfQeb9++8OWXevlZFH3XK8Zx+rRedn70UZg8WdfwFk4jM2AhCiM8PM8WhFYrTP/sEt79++KhdIA1TB4ElL6Q9Rintwq8WqVK+s13/nw9C5Lge+OoVu3KJeikJJgyRSdo+fvrLYf//U+CrwtIABaioE6f1jOL6wTg/VvOc9ddMOj1EiyuMQjlXSzHXr1ObRWYac0auOceuHBB92tdvhx69877eaJo6d9f5xUArFqls9lfekn//kF3MBIuIQFYiILKqICVUwA2p1pZ2Hs6FVrVoFzECmbOhL6HQ/AIW5VjkpVTWwVevAivvKKL50dF6b0/ceO6+27dm/n556FjR70vvGqVTrQTLiV7wEIUVOXKOpBdHoDXr+fEV3+QMG8VPVLC2V2hAzMW1KB8m4z727bNMbt5eNcGV+wBg4PKSC5bppsnHD0KL7+sk65KlrTvNYR7SUnRLQmnTYNhw3QhleLFXT0qgQRgIQouMBA+/zz7+/XrsbS/i8oWM5WB/Q+OpMnvH9nUls/hZSQhu3lCsWLSPOFmYrHoJMFNm6BVK1ePRlxGArAQBbV3L9Stm5UBbV0VhmFJxwCUyUT9W0tfEXxDw2OvG2Dt3iow04IFcOutOuFm1ix9NtTHx/7XEe6pRIkrPygKtyF7wEIUxMWL0KSJnlFmiK3bgRR8sHpcm2SVecwoNjEZRfYxI4ee9Y2Ph8ceg/vv1+UkQWc7S/AVwi1IABaiILZv10u6l+3/br3UkO4s5MQL1yZZOfWYkVLw66/QuLEuIfnhh/DJJ/a/jhCiUGQJWoiC2LZN/31ZAC47I4TFTIKPzkNp7yse7tRjRl98oatY3XYbfP+9DsRCCLcjAViIgggP18UqLjtDWSIqnKM+DWhwVfAFfZwoNodga7djRkrB2bO6bnP//rpu85AhUkpSCDcmS9BCFERmBazLkqyqnwknrmpwjg8f3rUBvl5XBkO7HTM6ehS6dtVFNdLT4ZZb9BEjCb5CuDWZAQtRECEhV/TEPbPrBBWtcewKzLkqlkOOGVmtum3cyJH6+08+kY5FQhQhEoCFKIjOna/4NvrPbZQDSt+Ve1lKux4zOnlSN01YvVrPfKdNgxo17PPaQginkI/LQuTX9u2wZMkVvYA3JzXhNSZTs3eQc8ZQtqy+/vffw+LFEnyFKIIkAAuRX9OmwcMPX7H/uzq6Jn/4v0b5WqUcd93du+Ghh7KbJ6xZA08/bVOlLSGE+5EALER+hYfrpvWX7beWWL2YDg1POuZ6ZrOu2RwcDGFhugIXSOAVooiTACxEflgsegn6svO/qScTmHrkXvpbfrD/9SIioHVreOcd3Spwzx79vRCiyMszABuGUc0wjFWGYew1DGO3YRivZtx+i2EYywzDiMr4u6zjhyuEi0VF6YbmlwXgYwsiAPC9PfcErAIbORJOnIC5c+G33/TZYyHEDcGWGXA68IZSqhHQBnjRMIzGwEhghVKqHrAi43shbmw59ABOWKlv8+9ppwC8aRPExOivZ8zQs94+fezz2kIIt5FnAFZKnVBKbcv4+gKwF/AH7gd+zHjYj0BvRw1SCLfx4IN6Wfiy8o7G9nBiDX9qti7k7DQ5Gd58U9eQfvddfVtAgC6sIYS44eTrHLBhGDWBFsBGoJJS6gToIG0YRo7vPoZhDAIGAVSvXr0wYxXC9by9oXnzK24qf2wbh/1a4F+YwlPr1sEzz0BkJDz3XHb3IiHEDcvmJCzDMEoCc4HXlFLnbX2eUmqaUqqlUqplhQoVCjJGIdyDUjBiBKxde8VND3qEsqrL+IK/7q+/wh13QEoKLFsG06frnr1CiBuaTQHYMAwvdPD9RSk1L+PmOMMwqmTcXwU45ZghCuEGrFb46itdgnL79qybjx2D8Iv1qNixaf5fMyVF/921KwwfDjt3XlNhSwhx47IlC9oAvgP2KqUmXXbXAuDJjK+fBP60//CEcAMHD0KnTvDKK9ClCzz+eNZd0T//y2Cm0qKp2fbXu3ABXnxRz3rNZr3H+/HHUMqBRTyEEG7HlhlwO6A/0NEwjIiMP92BCUAXwzCigC4Z3wtxY4mMhGbNdP/f6dN1CcrLlod95/3CON6mSXMbcPexKgAAIABJREFU0ymWL4fAQN1EoX37K8pZCiFuLnm+ayil1gC5ldzpZN/hCOEmkpKgeHGoX1+fxX3qKahW7ZqHlTkUTmTxFrQtmUdVqosXYehQHcTr19dlJG+/3TFjF0IUCVIJS4jLWSzw6adQsyYcPqzLPb77bo7BF7OZaud2cjrAhvO/Xl6wcaNO4oqIkOArhJAALESWffv0vuywYfosro/PdR9+YfM+iqlUrEHBOT8gMVG/1vnzunnCpk16r9fX1wGDF0IUNRKAhVBKN7MPCoL9++GXXyA0FKpUue7TYlZEAlC2Yw4z4L//hiZN4LPPdAMF0EFYCCEySAAWwjDg0CHo2VOXfezXz6ZOQ8v9HqIMidTt0SD7xrNnoX9/uO8+KFdOLzv36uXAwQshiqp8VcIS4oZhNsOECfoMbuvW8OWXep82HyIioFiFMlTxv+zGnj110B09Gt5+W1fOEkKIHEgAFjef8HDdyH77dkhN1QE4n8EXq5VH5z9K1RoDMIye2bcvW6a7F9Wta98xCyFuOLIELW4eqak6o7l1a4iL0/u8Y8cW6KXM+w/TJeF3mlU8eeUdJUpI8BVC2EQCsLh5fP+9Drj9+sHu3XD//QV+qRMLtwFQsv1lCVgLFugzw+Z8VMUSQty0JACLG1tyMuzYob8eOBBWroQffyx0i7/z/4ZjxpMaPS6rAT1/PsycCZ6ysyOEyJsEYHHjWrNGtw7s1k0HYk9PuPtuu7y0585w9hqNqR942dGi8HAIDrYpg1oIISQAixvPpUu6ccKdd+rl4J9+snvxi9OXfNlXvn32ZDc1VS9rt7ChKpYQQiBZ0OJGExenq1gdPgwvvwwffQQlS9r1EkpBH+bRqxf0zbxx1y5IT5cALISwmQRgcWOwWsHDAypWhB49oG9f3W3IAU6cgPh4vbp9xY1+fnoJWgghbCBL0KLoW7RIl308dEjvv375pcOCL8D50RPZSOsrewD37KmrYNWu7bDrCiFuLBKARdF19iw8+SR07w4mk2755wwb1lOWBAKDryreYRiSgCWEsJkEYFE0zZ8PjRvDrFm6uMbWrdCsWa4PDw2Ppd2EldQauZB2E1YSGh5b4EuXPRJOZIlgypTJuMFigdtu02MRQggbyR6wKJqWLYOqVWHxYt3F6DpCw2MZNW8nyWYLALGJyYyatxOA3i38r/fUayUkUOnSYRIaD8q+bf9+3WpQCnAIIfJBZsCiaFBKtwnctEl/P3GibnqQR/AFCFkSmRV8MyWbLYQsicz3MJLXRwBg3HpZslV4uP5bErCEEPkgM2Dh/mJiYPBgWLhQ7/m2bg3Fi9v89OOJyXneHhoeS8iSSI4nJlPVz5fhXRvkODs+cKIEe3mYcp0vO260bZvu9duwoe3/TUKIm57MgIX7UgqmT9cZzitXwqefwnff5ftlqvrlXIQj8/bMJerYxGQU2UvUOe0Tr0tvzSPMofFdFbJvDA+HwMD8d1QSQtzUJAAL9/XLLzBokF7a3bEDhg7V2c75NLxrA3y9rnyer5eJ4V0bAPlbok79awnv+4yneuz67BsbNoT77sv3uIQQNzdZghbuxWLRVazq1oVHHtGzyocf1kU2CihzKTm3JWZblqgBWLKElxd2Q2FgdPaBFSt01a0pUwo8NiHEzUsCsHAfe/fCs8/qghqRkVCmjA7CdtC7hX+uGc9V/XyJzSEIX7F0vXAhql8/ADxQkJYGYWHQsqVu8iDnf4UQ+SRL0ML1zGZdszkoSAfekBAoXdppl7/uEvXp07p/cM+enFblSKUYVg8TeHtDhw7w4Yfg7y9HkIQQ+SYBWLhWQoLOan77bejVC/bsgf79nTqj7N3Cn/F9AvH388UA/P18Gd8nUM+YLRZSFq8ipNT7VLu4jxn9VqHe/zB7+Tk8HMqWlQQsIUS+yRK0cA2ldJD184Nbb9XVrPr0yfGhth4RKowrlqiPHoWvv+BkpfG8/GolFiYcpF6z4qz5Dlq2bAu0zX5ieLieCQshRD7JDFg437p10KpVdvOEGTOuG3xtPSJUaGvWQM+eqIYNMX8+hQcb7eavv+Ddj4qzZYve7r1CfDzExkoLQiFEgUgAFs5z8SK8+irccYcOXqdO5fkUe1axuq6ff4a77oKFC7GmpPFw6v8wNQ9k+3YYNSqXFebMClgSgIUQBSBL0MI5li+HgQPhyBF46SWddFWqVJ5Ps/mIUGFYraihQ8FqxQCsGIy6fx+t5uVx+snfH954QwKwEKJAJAAL55g/X2cOr16tZ8A2sumIEAXcJ964EZo1Y/t+X+aVmMib8c/jjRmTjze3vdkh7/WhJk10TWohhCgACcDCcUJDoUoV3arvk/+3d+dxVo/9H8df12w1hfaUUTeFbi22u1JKJmmRkJDQz3LLliiUW4VK2tGN0nITcbuVUFkrlRZMaKFtRJuUVKKFmmY51++P6wzTmOU0c875njPzfj4e85izfM8513e+c87nXNf3uj6f0a47mZh3Wsj89Gtf76hKRnB0FisoQrWjgwdhwADs+PHMv/BxOn72CJUr30SroadzccwiTOtkN8M59+7kCvJDT83i4isvhLJlj2mfRERA54AlFHbtgq5d4aqrYOxYd1v58sccfKGQJUJ+AZ8nTkmBW26B007Djh/PqxV60WVJb7p3d6uf2jzSHDOgf77BN+dksP279nJxt3akPvDoMe+TiAioByzBZC28+ir06QOHDrnzvH37FvtpC8piBQGeJ05JcZOsMjLwYbiTicyveAdvvwFt2xbehtxB/szdmwF48VAlNAgtIkWhACzBM326Kxd4wQWualGYyvMVeJ7YWkhLg0WL8GVmEQNkEUPnlnv59xzXMQ9E7iDfYJcLwEvLn1zc5otIKaUhaCkenw82bnSXr7nG9YCXLg1rbdz8Ukk+1qg8XHYZadffyiPzk0mzZcgklpgyCVw2Ojng4At/nfTVYNdm9pSvSNzJwU0IIiKlhwKwFN0330CrVtCiBezf74oSdO9erMpFRZH7PHGtExKYnv4l7a5rQ8bCJQyZ15wxS5vxxh0LYOhQYj9ekOd53oLkDvINdm8itcZp9OsQvi8aIlKyaAhajl1GhiuYMGSIG8MdOzasxRPy0vncJDqnbYO3P4TX58DatXxZpT3X7p3IyS1O4esX4O9/z5VG8hifH/4safj8lb247vxTgp4SU0RKDwVgOTa//gqtW8PXX7s6vc8+CzVqeN0qN8mqTRtsejpk+RgR9xgj0gYzarzhrruC0yk/ejLYZcV/QhEp1RSAJTA5iyc0bQqDBrllRpFg2TK4917skXSML4sMYkmqW5b1Hxlq1QrB661c6XJAd+wIsbGFby8ikgedA5bCLV7s0i1u2uSC8OTJkRF8/bml7QUX8FvqNtJ88WQQi0lI4KYpyaEJvuCKR3TvHtaSiSJS8igAS/7274c773Tl9g4ehF9+8bpFf5ozBxo0wD73HP+r0JOahzYx5tKFpA8cStyiBZgLinauNyCrVsE554R9spmIlCwagpa8vfMO3H03/PSTKzjw+ONQrpzXrXLS0/H17MXuA+W42i5lR4UWvDkN2rcv+iSrgGVlwerVrrCEiEgxKABL3ubNg6pVXT7nJk1C+lIBFVL47DMYPx569OC931sz6tCHfLmvNnf3KcPQoXDccSFt4p++/dZl+VIFJBEpJgVgcbLTSJ5xBjRr5oonxMfnUwg3eAIqpDBzpkvy4fOR+foMhtnF/NawOYtnuzoPYaUawCISJDqJJa5Gb4cOLo3k5MnutnLlQh58oZBCCj4fjB+P7dYN6/MBYK2PYZcsYsUKD4IvwHXXwbp1UL++By8uIiWJAnBplpUF//63q2v72Wcwbpyb4RtGBRZS6NMHevXim7LnkEZZMokltmwCFz+eTEJCWJv5p9hYF3zjNHgkIsWjAFyavfoq3H+/m+W8bh3cc0/YZ/bmzrEcn5XBCWm/UfOERF4pdxe3J0yladYy3n9gITHDhhKz8NjTSAaNte7vtXSpN68vIiWKvsaXNkeOuIlEjRq5taxVqkCnTp6tae3Xvh79317DmVvXctXaj7lwy0rWVqnPA+U+5Oa1CXTsWJ91E6B27TDMcC5IRgY89ZQbMahXDy680Lu2iEiJoABcmnz2GfToAXv3wubNLo/z5Zd72qTO5yZR9ctPaTb8YWJ9WVjg+QMP8GuVBP73P+jWLQLyXSxb5tZDr14NV1wBN97ocYNEpCTQEHRpcPAg3HsvtGwJv/8OL78ceCHcUFu5kpb/upM4XxYGyCKWxmdnkpoK118fAcF31SpX33jvXjcbe/ZsOP54jxslIiWBesAl3U8/udzN27e7IDxsWBgXzRbutwpJ7IupRVUOE0cmMWUS6Pp8MlT1sFHWuhrHp5/uMl49/zzccIPnFZ9EpGRRD7ikyshwv0880VUt+vRTeOYZ74OvtTBtGlx1Fe+/6+PM5BOp/etqJlz7MVmDilarN6i2bnXD8mefDd9/77rgd92l4CsiQacecEmTnVBj4ED4+GM47TQ3echrKSkuq9ann8Knn7KpchNunfUz1RtUZ8YMaNbM40lWmZlugtWgQS7oDhsGSar1KyKhowBckmzZ4npr8+a585bWet0i59NPoXVrrL9XPrnMffQ+8DT9B8fSvz/erenNduSI+3utXOl6v+PGQe3aHjdKREo6BeCS4plnYMAAt453/HiCVoU+GBYswGZkYIBMYvFVr8GKD2Np0MDjdqWnu+hfpowLvAMHujKLns/8EpHSIEI+oaXYNm6Eiy+G9euhZ0/vg296OjzzDFkHD/HGr21JoywZxGLjE7jjf8neB99Zs9zwfEqKuz54MHTpouArImFTaA/YGDMF6ATsttY29N82GLgd2OPfbIC19oNQNVLycPgwDB0Kl10GLVrA00+79IiREEC+/BJuuw3WrOGJ5yozeNP/0bbWTC6u+SabzzubjuVr09mrtm3fDr16ueVEZ50VAePfIlJaBTIE/TIwDngl1+1jrbVPBr1FUrjFi1092u++c8OnLVqEpXBCoX7/HR59FPvMMxwsX4NbYmczf+9l1LzyazbUy+JbcxUAS3NXOwqXSZOgb1+XA3vUKJdWMhL+biJSKhU6TmmtXQL8Eoa2SGF+/dUF3uRkF0Q++sjN2g2hWat20GLkQk59+H1ajFzIrFU78t4wJcVNZBo7lukV7qDWwfWU63YFf7/nMxL+vv2ojvkf1Y7C7cABl4xk3Tp46CEFXxHxVHFOFPYyxqw2xkwxxlQKWoskf//9L7z0EvTrB2vWwCWXhPTlsmv17th3GMuftXqPCsJ798IHH2DbtMG3ei1plGFa/E28/n4F/vtf2JN5IM/nzq8KUlD9/rv7W82Y4a4/+CB88AGcemroX1tEpBBFDcATgLrAOcBOIN+FpsaYO4wxy40xy/fs2ZPfZpKfHTtgyRJ3+e67XWrE0aNdvd4QK7BWr7UwfTrUr8+BXv3JOpxODD7iTSbT715Ex45u+9zVjrLld3vQfPihK7P45JPubwZuYloknCMXEaGIAdhau8tam2Wt9QH/AZoWsO1ka21ja23jatWqFbWdpY/PBxMmuNqzN93kEkXExbkqRmGSXy/Vt+0HuPJK6NaNzZm1uGfLg2SYBHwxrl5vmfbJf2zbr309EuNjj3p8Ynws/drXC02jf/rJVXDo2NF9SVmyBIYPD81riYgUQ5ECsDGmZo6rVwFrg9McAdxSolat3HKiJk1gwQJPCsDn1UttsfUr5k/pSebc+TxW7kkaHFhG3UE3EfvxAmKeGOramiOVZOdzkxjRpRFJFRMxQFLFREZ0aRS6CVhLl7qiCUOGuJ6vygaKSIQytpBsScaY14FkXHr8XcAg//VzAAtsBe601u4s7MUaN25sly9fXqwGl3ipqS4P8fHHu6VFN93k2bBp9jngM7eupdm21SyrfRbby9Rl+KxXuHPvk1Q7vy4vvAANG3rSvD+lprpSgddd54bGt2+HWrU8bpSICBhjVlhrG+d5X2EBOJgUgAuwa5crnGAtjB0L3btD9epet4olL75N8zu6EufLIi2mLJfGz+fL2BYMH+6W08bGFv4cIZOWBiNGuJ/q1WHTJrcsS0QkQhQUgJUJy2v797uh5jp1XDYrY+CBByIi+LJ8Oa0euYd4f63eOF8GN9Vawrp10Lu3x8F38WJXKvDxx6FrV5fHWcFXRKKIArCXZs1yk6wmTYI77oAaNbxukXP4MPTtiz3/fH4/kMkREsggFpOQwK1TkznlFI/bt2WLS7uZng5z5rjlWZHwhUVE5BioGIMXfD53vvLNN106xFmz3GSrXGat2sGYuRv4cd9hTqqYSL/29YI+eSnP1/h7ZQ7PeJd3K9zO7b+Ook/b9fRrsojjOiV7V6vXWlixAho3dut433oL2rULy3IsEZFQUAAOJ2vdEHNMjBtyHjnSDTfnkZEpewJU9jrc7CQYELwUjjknWXXdvILa+35i0O4HePGn1izatoKKJx/Ha69Ap04e1+rdutUN08+Z4/JM/+Mf0NmzbNIiIkGhABwu69e7EoHDh7t0iKNGFbh5QUkwghWAx8zdwJlb1jLt9f7E+zIB+HRjO55Kv5yePd3cphNOCMpLFU1mJjz7LDz6qPviMnasO+8rIlICKACH2pEjLpINH+6WFu3dG9DD8kuCEcwUjlk//MCoOc+S4A++mcRSJvYINW78jPHjLwja6xSJte4879KlruLT889D7dretklEJIg0CSuUPvnE9diGDHEzdVNTXQapAIQjhePTCyZQ65efSCeeDGJJj4lnfecTqNMwLWivccwOH/5zqL57d3jjDXj3XQVfESlx1AMOpS++cGtVP/wQOnQ4pof2a1/vqHPAEKQUjt9+CxUqsD3jREYf/yJrbFmSqmyiwylvsvzMM/jmlPqMCFWayMLMmeOG6UeNcpPU7rjDm3aIiISBesDBZK2rvDNrlrt+332wdu0xB18IQQrHjAwYORJ71lmkdn6Y+vXho/UNuOiBKmQ+EMd/LrmSXQ3OC22ayPzs3g033giXXgply8LJJ4f39UVEPKBMWMGybRvccw+8954LJB984HWL/vTSSzBwIOzcyaIqV3P93udo1LYmkyZFQGW+GTNcr/fgQRgwAPr3V0INESkxCsqEpSHo4srKgnHjXICzFp56yvV8I8Xjj2MHDQIgnQRGpD/IyJdrepli+q/OPBP+8x/3W0SklFAALq7586FPHzfMPGEC3qeJ8ktLg7Jl+XHLEWpgiMESZ7J4895FHH9zcNf0HlPCkIwM9yWlbFn3d7vmGrj6arc2WkSkFNGnXlEcOgQff+wut2vnLn/wQWQE319/hR49yLqoNQ/2yeLaqZ04Qtk/avUe3yk5qC+Xncxjx77DWP5MGDJr1Y6/bvzFFy6TVf/+sHz50YlJRERKGX3yHau5c139vcsug59/dgEkOTkyxnPfegvq18f30su8sKEV457J5Kw7m5M1L+9avcFQUMKQPxw86Hq7zZq5ddAzZ7r8zZHwNxMR8YiGoAO1e7dLG/naa1CvnltaVLWq162ClBTX+16yBJYs4fvK53CV731+P/E8PnoHWrUCaA5tQ5NKMqCEIevXu/Pkd98dAem1REQigwJwIPbvd73efftg0KDImambkgJt2mDT08Hn46Wyd3PP/md4YEA8jz7qTrOG2kkVE9mRRxBuEHsYXn4ZbrkFzj/flVqMhCF6EZEIoQBckD17oFo1qFABBg+G1q0jZ6buxo3Qsyf2SDrGl0UGsRyqXItlH8Rz9tnha8ZfEoZYyw3rFzJ4yRRIP+LOkZ90koKviEguCsB5SU932ZiGD4ePPnLFE3r2DPrLFKncYGYmPPUUdvBgMnyx+HxxxALEJ3DXtGTiwhh84c/KTGPmbiB+yyaeXDCBxptWub/Z5Mku+IqIyF8oAOf2yScuBWJqqsvfXLduSF6mSOUGV62C226DVatYWuUquu0dx9WNv2fwRYuocnWyZ7V6O5+bROd6laB2V7fMaMIE9zfU7GYRkXwpAOf0wAOu5N3f/gbvvw8dO4bspQIuN5iSAosWQXIyvgEDOfTdTnrEvcncrKt5egrccstJGONhrd4NG+CMM6BcOZgyxdXqTQpzKksRkSikAJyditMYF3gffNBVLypfPqQvG9Ds4ZQUd945M5OsuAQerPEaU39L5pJrKpH6HNSoEdImFuz3392EtLFjYdo0uPZauOIKDxskIhJdSvcY4datbj3v66+76717w5NPhjz4QgDlBvftcyktjxyBrCx8R9Kp/ss3vDSzEjNmeBx8582DRo1cRqvbb4e2bT1sjIhIdCqdATgjwwXaBg3c+tnDwStyH6h+7euRGB971G1/lBucORPq18euWEEmcWQQi41N4N43k+ncOexNPVrfvtC+PcTHw+LFMHEiVKzocaNERKJP6RuCXr4cevSAr792Q6bjxkGtWmFvRs7Zw0fNgv7kbbjvPrZVPpur7LvUPTmdkR0WUeefySR4NMkKa91PTAxccIFbYPzII+FZaCwiUkKVvgC8bZtb3/v229C5s6fpEDufm+QCsbXwyy/YylV4Z2U3Vh13hBH7enP/w/E89hgkJno4yer7710Gq4sugn/9C7p0cT8iIlIspWMIeuZMtzQG4Kqr4Ntv3e9IyEU8YwacdhoZF7Ti6s5ZXNmjGu+c0ZeU5fGMGAGJeZ8qDr2sLHj22T+H6ZU+UkQkqEp2D/iHH+Dee2H2bJcO8c473TBqECZZFSmJRk6ZmdCnD3b8eAAs8fyyZRmjR7fg/vshzssjk5oK//wnLFsWeWUWRURKiJIZgLOyYPx4GDjQXR492lXjCVJiiECTaOQbpLdtcz3wlSsBMEAMPt66bwlV+rUIShuLZf9+2LzZFZ64/vrIGCkQESlhSuYQ9Nq1LuC2bAnr1kG/fm7WbpAEUoKvoDq5GRWr8cMv5RkU+wSHSXS1ehMTXDarIJu1agctRi7k1Iffp8XIhXnX6QX47DMYOdJdbtbMLdG64QYFXxGRECk5AfjgQXjjDXf57LNhxQpXpu/UU4P+UoEk0cgdpM/ftoZJ/x3AyEmbaHpRIrW3LmbtlQM59I6r1WtCUKu3oC8BfzhwAHr1cl9WJk1yf0fw8OSziEjpUDKGoGfPdkHkxx+hcWOoUwfOPTdkL5dfCb6cyTWyg3HLLavou+QVzvnpO7aUqcX+SVXYVwPeesv4JxM3h8tDM8u50HSX778Pd90FO3a4pB9PPAHHHReStoiIyNGiOwBv3+4mWc2a5er1vvGGC74h9pcSfORIouF3UsVEOs15lYcXvwxABnHcemQqvzSpxIZ5wcldUdhEsAJ76rt3u2ITp57qZmI3a1b8BomISMCidwg6LQ2aNIG5c925y5Urw1YNqPO5SYzo0oikiokYIKliIiO6NDoq+PVrdwbdV30IuElWYLn0nBlMmmSDFnwLG17+S7pLa2m5ZRUnVSgL1avDggXu76bgKyISdtHbAy5bFp5/3p3vLUKvt7jLiP5IopGTtTB1KrRvT9aWk3nQ/oeXuY540smKi6NZ7w5cdCxLlQoQSDWlnD31pP27GT53HBdtWUlKy6nuAQq8IiKeid4ADG4pTxEUqRZvYTZvdjVwFyxg+t8fo9s3QzjnnMv5ofcC6u1cRHxyMhcFsYceyESwzucmQVYWW4eM4va5UzDG8PVDQ2l+T/egtUNERIomugNwEQVci7cwKSmwcCH89BP2xRfJsPE8VHYik7fczsiRrrxwfHxzIPhD44FMBAPoPLwPvDcTLr0UJk7k7Nq1g94WERE5dqUyAAdUi7cwKSnQpo07F20ta05oyaUHpnFGchJfT4bTTw9SY/NR4ESw9HS3fjc+Hm69Fa6+Wmt6RUQiTPROwiqGQmvxFiYtDWbPxqang7VkEcPMI5cyeHISCxaEPvhCARPBjvwA553navUCXH453Hijgq+ISIQplT3gQJYR5WvJErj9dg5nxmJ8CcSSji82gZ5vtKbaFSFsdB6Omgj222+uROCzz0JSEpx1VngbIyIix6RUBuB8a/EWdP53/354+GGYOJG9FU7lhoPjia9UnpEdFtGwVzLVvKrVC7B0Kfzf/7nSgT17wogRql4kIhLhSmUAhnyWEeVn3Tpo3x67cydTKj7Affse5/rbyjNmDFSq5GHgzZaQ4Co8LV3qUkqKiEjEK7UBOCCffQaLF3OwUXO+jW/G3b6H2Fu5Ke++BRdf7GG7rIU334TVq2HoUFdqcc2aoFV7EhGR0FMAzou18NhjMGwYPhNDnC+B+8wCkvs1ZfBgKFfOw7b9+CPcc49Lv9m4sSu5WLasgq+ISJTRp3ZuW7ZA+/bwxBNYa4nxZRFPOm/0XMTo0R4GX2vhhRegfn2YM8fVOE5JccFXRESijgJwtqwsGDsW27Ah6UuX8WxCXw6TSJZxtXqTbkz2tn0//gi9e7vUm6tXuxrHcRrAEBGJVvoEz2YMh6fN5qvEi+m693nqtKpF555dqL15ESQnh63Qw1GystxQc5cubmnR55+7HrCGm0VEol7pDsCLFsETT5B5z32M3XgFo79+l/SE4xgzydCjB8TEhCaNZEDWroXbboMvvoD5813WrYYNvWmLiIgEXekNwBMnujWz1mIXLuFtu5gWVzZn/HjX2fRMerpbxztsGFSoAP/7n8dTrkVEJBRKXwA+cMAl1JgwAYu/Vq/1MfmGRTT8b3PvMzZ26gQffeRyN//731CtmscNEhGRUCh9JxMffRQ7cSLvHdf1j0lWcYkJNOqV7F3wPXQIMjLc5fvvh3ffhddeU/AVESnBSmwPeNaqHX+kmqwfl0afJifSNLklQ/Y9yqf2BvZUO5/pw1Jo8vsi7yZZAXz8MfTo4X7693dlA0VEpMQrkQF41qod9H97DWduWcuA5bO4cMsqUk84gzpmOfv2VuX+B6syZAiUL+/hJKv9++Ghh2DyZKiibgRiAAALpklEQVRb17svACIi4okSGYDHzN1Aq7VLeH72SGKtJQvDv3/uR1qNQyxbdjxNmhT+HDl70AEVazgWH3/siifs3Al9+8KQIR6n1xIRkXArkQG4xtqVPPfOaGKsBcBHDGeespyUayrSpEnHQh+f3YPOLle4Y99h+r+9BiA4Qbh8eahaFd5+G5o2Lf7ziYhI1ClZk7COHAHgx+r/4IMy7UijLBnEkhkbx1ctTyKpSmBpG8fM3XBUrWCAwxlZjJm7oWjtshamTYMBA9z1pk1h5UoFXxGRUqxk9ICPHIFhw7DTpvHcLStZOeEirvZdyCVNp3Fx2Y/4vHYjUk9pyIj29QJ6uh/3HT6m2wt+sh/h7rvhnXdcwE1LU/EEEREpPAAbY6YAnYDd1tqG/tsqA9OBU4CtQFdr7a+ha2YeUlJcJqsqVdx62dRUPqzcnccGZtL+8hiu6LmLl746iYn7unJSxURGHMM53JMqJrIjj2B7UsXEwNtnLUyZAg8+6L4gjBkDffoof7OIiACB9YBfBsYBr+S47WFggbV2pDHmYf/1fwW/eflISXGpGdPSsNZysNyJXB/zIcvjOjB5Olx7LRhTkx4dahbp6fu1r3fUOWCAxPhY+gXYgwb+LJ7wj3+4Kkann16ktoiISMlU6DiotXYJ8Euum68EpvovTwU6B7ldBVu0yKVstBaL4elDd1H9pg6sXw9du1LshBqdz01iRJdGJFVMxABJFRMZ0aVR4T1onw9mznS936QkWLbMzXhW8BURkVyKOh56orV2J4C1dqcxpnoQ21S45GQyYhIgK51Mk0DHse1p2ju4L9H53KRjm/G8YYMrnvDppzB3LrRrp+IJIiKSr5DPBDLG3GGMWW6MWb5nz57gPGnz5qQ+t4B5LYbC/AU07e1hEovMTBg1ytXpXb8epk6Ftm29a4+IiEQFY/1rZQvcyJhTgPdyTMLaACT7e781gUXW2kJPkDZu3NguX768eC2ONFdc4XI3X301jBsHNWp43SIREYkQxpgV1trGed1X1B7wO8DN/ss3A7OL+DzRKT39z+IJd90FM2bAm28q+IqISMAKDcDGmNeBFKCeMWa7MeY2YCTQ1hjzHdDWf710+PJLN7N5zBh3vWNHuOYab9skIiJRp9BJWNba6/O5q02Q2xLZDh+GQYPgqaegZk13zldERKSIlBUiEF98Ad27w3ffwe23u95vhQpet0pERKJYVAbgkFYqyosxbm3v/PkuAYiIiEgxRV0ADnmlomzz57s1vYMGQZMmkJqqNJIiIhI0UVcRIOiVinLbv98NM7dtC6+/Dr/95m5X8BURkSCKugAc1EpFub3/PjRo4IooPPQQrFoFxx1X/OcVERHJJeq6dUGpVJSXn3+G666DU05x+ZybNCne84mIiBQg6nrA/drXIzE+9qjbjrlSUU6ffOImWFWtCgsWwIoVCr4iIhJyUReAi1ypKLc9e1yP98ILYdYsd9v550OZMkFvs4iISG5RNwQNRahUlJO1MH063HsvHDgATzwBnToFt4EiIiKFiMoAXCw9e8LEiW6Y+aWX3KQrERGRMCsdAdha8PkgNhYuuwzq1IH779fSIhER8UzUnQM+Zjt2uJKBI/31Ijp1gn79FHxFRMRTJTcAWwsvvgj167vZzRUret0iERGRP5TMbuC2bS6b1bx5cNFFLhDXret1q0RERP5QMnvAu3bB55/D+PGwcKGCr4iIRJyS0wPesgXee88tL2rSBH74AY4/3utWiYiI5Cn6e8A+H4wbB40awSOPuN4vKPiKiEhEi+4AvHEjtG7ter0tW8KaNXDiiV63SkREpFDROwR9+DC0aAFHjrjqRbfcAsZ43SoREZGARG8ATkx0mazOPhuSipiWUkRExCPRG4ABOnb0ugUiIiJFEt3ngEVERKKUArCIiIgHFIBFREQ8oAAsIiLiAQVgERERDygAi4iIeEABWERExAMKwCIiIh5QABYREfGAArCIiIgHFIBFREQ8oAAsIiLiAQVgERERDxhrbfhezJg9wPdBfMqqwM9BfD4vaV8iT0nZD9C+RKqSsi8lZT8g+PvyN2tttbzuCGsADjZjzHJrbWOv2xEM2pfIU1L2A7Qvkaqk7EtJ2Q8I775oCFpERMQDCsAiIiIeiPYAPNnrBgSR9iXylJT9AO1LpCop+1JS9gPCuC9RfQ5YREQkWkV7D1hERCQqRUUANsZ0MMZsMMZsNMY8nMf9ZYwx0/33f26MOSX8rSycMaaWMeZjY0yqMWadMaZ3HtskG2P2G2O+8v885kVbA2GM2WqMWeNv5/I87jfGmGf9x2W1MeY8L9pZEGNMvRx/66+MMQeMMX1ybROxx8QYM8UYs9sYszbHbZWNMR8ZY77z/66Uz2Nv9m/znTHm5vC1Om/57MsYY8w3/v+fmcaYivk8tsD/xXDLZ18GG2N25Pg/6pjPYwv8vAunfPZjeo592GqM+Sqfx0baMcnz89fT94u1NqJ/gFhgE1AHSAC+Burn2qYnMNF/uRsw3et257MvNYHz/JePB77NY1+Sgfe8bmuA+7MVqFrA/R2BDwEDNAM+97rNhexPLPATbt1eVBwToBVwHrA2x22jgYf9lx8GRuXxuMrAZv/vSv7LlSJwX9oBcf7Lo/LaF/99Bf4vRsi+DAb6FvK4Qj/vvN6PXPc/BTwWJcckz89fL98v0dADbgpstNZuttamA9OAK3NtcyUw1X/5TaCNMcaEsY0BsdbutNau9F8+CKQCSd62KqSuBF6xzjKgojGmpteNKkAbYJO1NpjJYkLKWrsE+CXXzTnfD1OBznk8tD3wkbX2F2vtr8BHQIeQNTQAee2LtXaetTbTf3UZcHLYG1YE+RyXQATyeRc2Be2H/zO2K/B6WBtVRAV8/nr2fomGAJwE/JDj+nb+GrT+2Mb/Zt0PVAlL64rIP0x+LvB5Hnc3N8Z8bYz50BjTIKwNOzYWmGeMWWGMuSOP+wM5dpGkG/l/mETLMQE40Vq7E9yHDlA9j22i7dgA/BM3opKXwv4XI0Uv/3D6lHyGOqPpuFwI7LLWfpfP/RF7THJ9/nr2fomGAJxXTzb31O1AtokYxpjjgLeAPtbaA7nuXokbAj0beA6YFe72HYMW1trzgEuBe4wxrXLdHzXHxRiTAFwBzMjj7mg6JoGKmmMDYIwZCGQCr+WzSWH/i5FgAlAXOAfYiRu+zS2ajsv1FNz7jchjUsjnb74Py+O2Yh+XaAjA24FaOa6fDPyY3zbGmDigAkUb/gk5Y0w87uC/Zq19O/f91toD1trf/Jc/AOKNMVXD3MyAWGt/9P/eDczEDZ/lFMixixSXAiuttbty3xFNx8RvV/ZQv//37jy2iZpj45/w0gm40fpPyOUWwP+i56y1u6y1WdZaH/Af8m5jVBwX/+dsF2B6fttE4jHJ5/PXs/dLNATgL4HTjTGn+nsp3YB3cm3zDpA9K+0aYGF+b1Qv+c+ZvAikWmufzmebGtnnr40xTXHHaG/4WhkYY0x5Y8zx2Zdxk2XW5trsHeAm4zQD9mcP9USgfL/NR8sxySHn++FmYHYe28wF2hljKvmHQtv5b4soxpgOwL+AK6y1h/LZJpD/Rc/lmv9wFXm3MZDPu0hwCfCNtXZ7XndG4jEp4PPXu/eL1zPTAvnBzab9Fjc7cKD/tsdxb0qAsrihw43AF0Adr9ucz360xA1brAa+8v90BO4C7vJv0wtYh5v9uAy4wOt257Mvdfxt/Nrf3uzjknNfDDDef9zWAI29bnc++1IOF1Ar5LgtKo4J7kvDTiAD9y39Ntz8hwXAd/7flf3bNgZeyPHYf/rfMxuBWyN0Xzbizr1lv1+yVzucBHxQ0P9iBO7Lq/73wWrch37N3Pviv/6Xz7tI2g//7S9nvz9ybBvpxyS/z1/P3i/KhCUiIuKBaBiCFhERKXEUgEVERDygACwiIuIBBWAREREPKACLiIh4QAFYRETEAwrAIiIiHlAAFhER8cD/A53BEXgj07ZTAAAAAElFTkSuQmCC\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.2834419617291377e-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.2834419617291078e-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.3031864410632063, 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.3031864410632063, 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": [
      "/home/travis/miniconda/envs/statsmodels-test/lib/python3.7/site-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:                Tue, 17 Dec 2019   Prob (F-statistic):           4.98e-10\n",
      "Time:                        23:42:54   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": [
      "/home/travis/miniconda/envs/statsmodels-test/lib/python3.7/site-packages/scipy/stats/stats.py:1450: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    }
   ],
   "source": [
    "ols_model = sm.OLS(y, X)\n",
    "ols_results = ols_model.fit()\n",
    "print(ols_results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Condition number\n",
    "\n",
    "One way to assess multicollinearity is to compute the condition number. Values over 20 are worrisome (see Greene 4.9). The first step is to normalize the independent variables to have unit length: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "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.8709117813\n"
     ]
    }
   ],
   "source": [
    "eigs = np.linalg.eigvals(norm_xtx)\n",
    "condition_number = np.sqrt(eigs.max() / eigs.min())\n",
    "print(condition_number)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dropping an observation\n",
    "\n",
    "Greene also points out that dropping a single observation can have a dramatic effect on the coefficient estimates: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage change 4.55%\n",
      "Percentage change -2228.01%\n",
      "Percentage change 154304695.31%\n",
      "Percentage change 1366329.02%\n",
      "Percentage change 1112549.36%\n",
      "Percentage change 92708715.91%\n",
      "Percentage change 817944.26%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ols_results2 = sm.OLS(y.iloc[:14], X.iloc[:14]).fit()\n",
    "print(\"Percentage change %4.2f%%\\n\"*7 % tuple([i for i in (ols_results2.params - ols_results.params)/ols_results.params*100]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at formal statistics for this such as the DFBETAS -- a standardized measure of how much each coefficient changes when that observation is left out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "infl = ols_results.get_influence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general we may consider DBETAS in absolute value greater than $2/\\sqrt{N}$ to be influential observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2./len(X)**.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    dfb_const  dfb_GNPDEFL       dfb_GNP     dfb_UNEMP     dfb_ARMED  \\\n",
      "0   -0.016406  -169.822675  1.673981e+06  54490.318088  51447.824036   \n",
      "1   -0.020608  -187.251727  1.829990e+06  54495.312977  52659.808664   \n",
      "2   -0.008382   -65.417834  1.587601e+06  52002.330476  49078.352378   \n",
      "3    0.018093   288.503914  1.155359e+06  56211.331922  60350.723082   \n",
      "4    1.871260  -171.109595  4.498197e+06  82532.785818  71034.429294   \n",
      "5   -0.321373  -104.123822  1.398891e+06  52559.760056  47486.527649   \n",
      "6    0.315945  -169.413317  2.364827e+06  59754.651394  50371.817827   \n",
      "7    0.015816   -69.343793  1.641243e+06  51849.056936  48628.749338   \n",
      "8   -0.004019   -86.903523  1.649443e+06  52023.265116  49114.178265   \n",
      "9   -1.018242  -201.315802  1.371257e+06  56432.027292  53997.742487   \n",
      "10   0.030947   -78.359439  1.658753e+06  52254.848135  49341.055289   \n",
      "11   0.005987  -100.926843  1.662425e+06  51744.606934  48968.560299   \n",
      "12  -0.135883   -32.093127  1.245487e+06  50203.467593  51148.376274   \n",
      "13   0.032736   -78.513866  1.648417e+06  52509.194459  50212.844641   \n",
      "14   0.305868   -16.833121  1.829996e+06  60975.868083  58263.878679   \n",
      "15  -0.538323   102.027105  1.344844e+06  54721.897640  49660.474568   \n",
      "\n",
      "          dfb_POP      dfb_YEAR  \n",
      "0   207954.113590 -31969.158503  \n",
      "1    25343.938291 -29760.155888  \n",
      "2   107465.770565 -29593.195253  \n",
      "3   456190.215133 -36213.129569  \n",
      "4  -389122.401699 -49905.782854  \n",
      "5   144354.586054 -28985.057609  \n",
      "6  -107413.074918 -32984.462465  \n",
      "7    92843.959345 -29724.975873  \n",
      "8    83931.635336 -29563.619222  \n",
      "9    18392.575057 -29203.217108  \n",
      "10   93617.648517 -29846.022426  \n",
      "11   95414.217290 -29690.904188  \n",
      "12  258559.048569 -29296.334617  \n",
      "13  104434.061226 -30025.564763  \n",
      "14  275103.677859 -36060.612522  \n",
      "15 -110176.960671 -28053.834556  \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/stats/outliers_influence.py:693: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return self.resid / sigma / np.sqrt(1 - hii)\n",
      "/home/travis/miniconda/envs/statsmodels-test/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py:901: RuntimeWarning: invalid value encountered in greater\n",
      "  return (a < x) & (x < b)\n",
      "/home/travis/miniconda/envs/statsmodels-test/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py:901: RuntimeWarning: invalid value encountered in less\n",
      "  return (a < x) & (x < b)\n",
      "/home/travis/miniconda/envs/statsmodels-test/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py:1892: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= _a)\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/stats/outliers_influence.py:733: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_internal * np.sqrt(hii / (1 - hii))\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/stats/outliers_influence.py:761: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_external * np.sqrt(hii / (1 - hii))\n"
     ]
    }
   ],
   "source": [
    "print(infl.summary_frame().filter(regex=\"dfb\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}