{
 "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:                Fri, 24 Apr 2020   Prob (F-statistic):          2.83e-239\n",
      "Time:                        14:17:04   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:                Fri, 24 Apr 2020   Prob (F-statistic):           6.30e-27\n",
      "Time:                        14:17:04   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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzde3zP9RfA8ddn313dr4m5q1xyv5SRDLmUa3RDKbr8lJKUUIqoEF3IJcolUaGYW+XaEHNfcldIGXLdmG12e//+OOY6s8v3up3n47HH9r19Pu/vxvd83rdzLGMMSimllHIuL1c3QCmllMqJNAArpZRSLqABWCmllHIBDcBKKaWUC2gAVkoppVxAA7BSSinlAt7OPFmRIkVM2bJlnXlKpZRSymW2bt16yhhTNLXHnBqAy5Yty5YtW5x5SqWUUsplLMs6fLPHdAhaKaWUcgENwEoppZQLaABWSimlXMCpc8CpSUhI4MiRI8TFxbm6KQ7j7+9PyZIl8fHxcXVTlFJKuQmXB+AjR46QN29eypYti2VZrm6O3RljOH36NEeOHKFcuXKubo5SSik34fIh6Li4OAoXLpwtgy+AZVkULlw4W/fwlVJKZZzLAzCQbYNviuz+/pRSSmWcWwRgdzJkyBBGjx5908dDQkLYvXu3E1uklFIqO/K4ABwSHkHDEasoN2AJDUesIiQ8wrnn1wCslFLKDjwqAIeERzBw3g4iImMxQERkLAPn7chyEP7ggw+oWLEiDzzwAPv27QPgyy+/pF69etSoUYNOnToRExPD+vXrWbhwIf369aNmzZocOHAg1ecppZRSt+JRAXjU0n3EJiRdc19sQhKjlu7L9DG3bt3K999/T3h4OPPmzWPz5s0AdOzYkc2bN7N9+3YqV67MlClTaNCgAe3atWPUqFH8/vvvVKhQIdXnKaWUUrfi8m1IGXE0MjZD96fH2rVrefjhh8mVKxcA7dq1A2Dnzp0MGjSIyMhIoqOjadmyZaqvT+/zlFJKubmoKMif32mn86gecIkCARm6P71SW6X8zDPPMG7cOHbs2MHgwYNvuo0ovc9TSinl5gYNggjnrSvyqADcr2VFAnxs19wX4GOjX8uKmT7m/fffz/z584mNjeX8+fMsWrQIgPPnz1O8eHESEhKYNWvW5efnzZuX8+fPX759s+cppZRyc1u2QKdOsHGj3B4wAJyYsdCjAnCHWoEM71iNwAIBWEBggQCGd6xGh1qBmT5m7dq1efzxx6lZsyadOnWiUaNGAAwbNox7772X5s2bU6lSpcvPf+KJJxg1ahS1atXiwIEDN32eUkopN2QMhIZCixZQrx6sXAkHD8pjgYFw221Oa4pljHHayerWrWuurwe8Z88eKleu7LQ2uEpOeZ9KKeXW2rSBJUugWDHo2xd69oR8+Rx2Osuythpj6qb22C0XYVmWNRVoA5wwxlS9dN8ooC0QDxwAuhtjIu3XZKWUUspOLl4EPz/5uVUreOgh6N4dArK2fiir0jMEPR1odd19y4GqxpjqwH5goJ3bpZRSSmXdtm1QvTrMnCm3X34ZXnrJ5cEX0hGAjTFrgDPX3bfMGJN46eYGoKQD2qaUUkplTnIyjB4N9evDhQtQ0v3ClD32AfcAZtvhOEoppVTWHT0KTz8NK1ZAx44weTIULuzqVt0gS6ugLct6G0gEbrr/xrKsFyzL2mJZ1paTJ09m5XRKKaXUrW3cCOvXS+D94Qe3DL6QhR6wZVlPI4uzmpk0llIbYyYDk0FWQWf2fEopZU8h4RGMWrqPo5GxlCgQQL+WFdO/pdEY2LkTFiyAv/+Gp56Cxo0lk9LKlfKBX6iQfC9c+MoCIOU4cXGwbh00awYPPwwHDsDtt7u6VWnKVAC2LKsV0B9obIzx6OoDp0+fplmzZgAcP34cm81G0aJFAdi0aRO+vr6ubJ5SygFSCruk5JZPKewCpB2Ez5+HwYMhJAQOHZL7ihWT4Auwd68kdriary/MmgWPPGLvt6FSHDsmQff33+XvUry42wdfSN82pO+AYKCIZVlHgMHIqmc/YPmlNI4bjDE9HdhOhylcuDC///47ILWA8+TJwxtvvHH58cTERLy9PSpltlLqFtIq7HJNAL5wAZYuhZgYePJJyJULfvwRqlaVrElt28qHfYpq1SA8HM6cgdOn5fvOnbIQCGDZMsm+9MwzUKKE499oTrB1K7RvD2fPwrffXvv3cHO3jCzGmM6p3J2tS/4888wzFCpUiPDwcGrXrk3evHmvCcxVq1Zl8eLFlC1blpkzZzJ27Fji4+O59957mTBhAjab7RZnUEq5UroKu8ydCy++KIG0Zk0JwDabDG3e7KI8Vy557s2sWgUjR8K770Lr1vDcc/Dggzc/nkrbnDlyMVOkiAw/p/W7d0Nu9Vfv00dGEOypZk347LOMv27//v2sWLECm83GkCFDUn3Onj17mD17NuvWrcPHx4eXXnqJWbNm0a1bt6w1WinlUCUKBBCRShAuUSBAeq0vvwzffSepCufOhUspaoGsBcsRI+DZZ2HKFJg+HRYuhOBg+PXXzB8zJ/vjD6hVC+bNk6kAD+NWAdidPProo7fsya5cuZKtW7dSr149AGJjY7nNiXlElVKZ069lxWvmgOGqwi7798P8+TBsmAwz27t3euedEoiHDYPFi2VBF0B8vPRA7rnHvufLbi5ckHneqlVh6FBISPDYRW5uFYAz01N1lNy5c1/+2dvbm+Tk5Mu3U0oOGmN4+umnGT58uNPbp5TKvJR53pRV0HcEGD70O0S9Wq2AQFnZ7OgelY+PLBxKMWkS9O4t+1dHjvTIHp3DHT4s870nTsBff8mQv4cGX/CwakiuUrZsWbZt2wbAtm3bOHRp9WOzZs344YcfOHHiBABnzpzh8OHDLmunUir9OtQKZN2AphxqlZvl01+h3qBXrl3Z7Gzdu0uP+9tv4a67YMwYSEy89etyit9+k9GBQ4dg6lQJvh5OA3A6dOrUiTNnzlCzZk0mTpzIXXfdBUCVKlV4//33adGiBdWrV6d58+YcO3bMxa1VSqVLfDy89ho0aSK90bVroVw517UnTx4YPhx27JBV0336SG84pzMGJkyQv1P+/JJko9X15Qk8k5YjdJKc8j6V8gjGyL7cefNkwdWIEXDVtJPLGSN7jUuUgHvvhXPnZK7TTTM6OZQxMlSfmCgFFQoUcHWLMiRL5QiVUirbsSz5UG/WTCrjuJuU9qXo318Whk2YILmNc4J//5XgW7q0DMv7+4NX9hq0zV7vRiml0hITI3OJIPt63TH4pqZnT+kNd+oEjz0mi5Cys9WroU6dK0PwuXJlu+ALGoCVUjlFdLQkv2jRAo4fd3VrMqZGDZn7/OADyT9dpYrknM5ujJHtMM2ayXD7xIl2O3RIeAQNR6yi3IAlNByxipDwCLsdO7M0ACulsr9z52Thzpo18NVXHpEn+AY+PvDWW5LqsmZNKF/e1S2yr5gYKWrx2muS4nPjRqhUyS6HTsn9HREZi+FK7m9XB2ENwEqp7C0yElq2lA/077+HLl1c3aKsqVJF6tyWKyc9xi5dYNo0KUDvyZKTJbPVsGGSbztfPrsdOq3c366kAVgplb1NmSIJ++fOhUcfdXVr7OvcOThyBHr0gIYNYfNmV7coY2Ji4L33pMpUnjzS/kGD7D7fm67c3y6gARg4cuQI7du3584776RChQq8+uqrxMfHExoaSps2bW54/uLFi6lVqxY1atSgSpUqTJo0yQWtVkqlS9++8sHeoYOrW2J/+fNDaKjklf77b0lU0b275LN2d6tWSfWoIUPgp5/kPgdltSpRICBD9ztLjg/Axhg6duxIhw4d+PPPP9m/fz/R0dG8/fbbqT4/ISGBF154gUWLFrF9+3bCw8MJDg52bqOVUmmLiYHOnWHfPtnSU6OGq1vkOF5eslp43z54800JbD4+rm7VzUVGShWoZs2k7aGh8PjjDj1lv5YVCfC5Nrf/5dzfLuSZATgsTDLGhIVl+VCrVq3C39+f7t27A2Cz2fj000+ZOnUqMTExNzz//PnzJCYmUvjShng/Pz8qVnTtH1EpdZWkJOjaFWbPhr17Xd0a58mXT3JI79sHefNK4o727a/0Lt3Fiy9Kj71/f5nzbdzY4afsUCuQ4R2rEVggAAsILBDA8I7Vrq397ALul4gjtd7kY4/Jfr2YGJnn+OMPmbD38oLq1eHVV6Um5KlTkt3maqGhaZ5u165d1KlT55r78uXLR+nSpfnrr79ueH6hQoVo164dZcqUoVmzZrRp04bOnTvjlQ33qCnlcYyRVbQhIZJLuX17V7fI+fz95fuRI3IB0ro1PPSQBLxGjWREwNk2b5YLhIoVZSvVG2/IPl8n6lAr0OUB93qeFzWioq6s9ktOlttZYIzBSuUf5M3uB/jqq69YuXIl99xzD6NHj6ZHjx5ZaoNSyk4++ww+/1yCcO/erm6Na5UrJ3mlR4+G9eulp1m5Mhw96pzzx8bK6ux69WRu+qOP5P7y5Z0efN2WMcZpX3Xq1DHX27179w33pWn9emMCAoyx2eT7+vUZe/11li9fbho1anTNfVFRUaZQoUJmyZIlpnXr1mm+/uTJkyZPnjy3PE+G36dSKmMSEoy57z5jOnUyJinJ1a1xLxcuGDN9ujGdOxuTnCz3ffONMatXX7ltT0OHGlOwoDFgTJUqxowbZ0xUlP3P4wGALeYmMdHzesBBQZIBZtgw+R4UlKXDNWvWjJiYGGbMmAFAUlISr7/+Os888wy5Uil3FR0dTehVw9q///47ZcqUyVIblFJ24O0Ny5bBN99ky7SFWZIrlyzU+vZbGYJOToZ3373SK/7kE9i5U3qtGWWMTP/9/POV0ckLF+CBB+DXX+W4vXrZdV9vdqHVkIB///2Xl156ib1795KcnMxDDz3E6NGjCQsL48EHH7y84Argu+++Y/jw4Rw4cICAgABy587NmDFjqFs31WIXl7nD+1QqW/rrL3j7bZg8WbblqPSJiZG90ZMmXVnQ+u67si83MlKybt11F9x5p3z395dKRHnzyr7qkSPhwAH5SpkKXLFCVjcb45q5Zjek1ZBuoVSpUixatOiG+4ODg4lN5YqwUaNGzmiWUupWTp2CBx+Es2fh5EkNwBmR0it++mlZrPX775JlC2QB13ffSSC+2uzZsig2Lk6eX6GCjEJWqCCvTVnR7IHBNyQ8glFL93E0MpYSBQLo17KiwxdtaQBWSnmm2Fho107K1q1aBXfc4eoWea5Kla7Nu1y1qiTzOH0a9u+Xr4sXZUEVyG6U/ftd01YHSMkVnZKuMiVXNODQIKwBWCnleZKToVs32LBBhlEbNHB1i7Ify4IiReQrm/9+08oVrQFYKaWudvSoBN9Ro6RGrgu4YshSOYarckW7RQA2aey5zQ6cudBNqRyhZEnZ4+qiOV9XDVkqxyhRIICIyFhqR+yh/j872FC6GtsCKzs8V7TL1+r7+/tz+vTpbBukjDGcPn0a/5TsNEqpzFu1SoorJCXJilwXXbi7a3k7lTn9Wlbk/oidzJ3Vn75rZzLr+7cJ+m+/w3NFu7wHXLJkSY4cOcLJkydd3RSH8ff3p2TJkq5uhlKebe9eGW4ODJQKOi7cV+qu5e1UBh04ABs30qFLF+4sEIllkrEBJCUyKM8J7s7uq6B9fHwoV66cq5uhlHJnJ09KPmNfX1i82OVJHVKGLFO7X7mxsDBJ4GRZUidgxQoICIC2bbm7SzuYMgbi4/H29ZXbDubyAKyUUmmKi5NavseOyYdm2bKubhH9Wla8Zg4Y3KO8nUpDWJgU+4mPl9vFisHQodCjhyQXScmyGBoqz8tilsX00ACslHJv4eGS9OGbb+Dee13dGuDKQitdBe0hEhIksCZdumDy8oJXXpEMalcLCnJK4E2hAVgp5d6CguDgQemxuBF3LG+nrnPuHLz5Jvz5p9QP8PWVHrCvLzRtevlpycnSQZ49W34eN845zdMArJRyTzNnSvalZ591u+Cb4swZGRm/7TYoVAhsNle3SF3200/wv//JnvE+faBu3WuGmE39IDZvkqA7d64kVPPzk5LyzkplrQFYKeV+li2Tubn77oNnnnGbyJaUJLXlf/kFli4Frw3racyv/EpTNnkFUbSoBOPbbpNrhnqJYdSNDqXi/4Ip2s55Q5s5WmSkDC/PnAl33w0//HB56sLUD+J3/yBmz4Y5XeHQIfDxgVatYPhwyWyaN6/zmqoBWCnlXtavh4cfluT+P/7o8uB77JgE219+geXLIfrMRZoQyocFJxFMCFgWid4fsLTJSBqsG8WFmDycP5iHpIQkqsT/jgHif/JjYOOVNH07iGbNtFqiQ4SFSe+2Xj3JkjZ4sFR08vUFZMFz376Sv8Vmk2qJ77wj6/sKFnRNkzUAK6Xcx/btst0oMFCinqs+GS81ZWK3MAr+EUoowZwregc/5XuZ2rE/4xN7HqJ9AAPG4JMcT5u79kOJZhSKjoboaC5s34nXsWQswBBPzfXjONRiGk8FdqVW70Y808OLIkVc9vayl9BQaNFCJnB9faU28aXKTEeOwOuvw5w5UL68VF/s2BG3+N1rAFZKuY/Vq2WP74oVLpv3jYuT9TprR/zGiuSmeJMIfv7ww1K8XtwJTz4hY5V58sjFQsqini5dLq+gDQmPYPbYOUydOQCfpEQSbTbiq53n6R0hvBDxJf/0L8X0tzpzplVXnmh7gWqnQ7GaBN+wAjen5JvO0vs8fBiefFJWOoP8PdavJz6oMZ99JjuNkpLke79+UtbYXVjOTAFZt25ds2XLFqedTynlIa5e9RIV5bIcz2vWwPPPg7V/L7/lbkmRC//IAzabROWBA699Qcqw53X7RhuOWJVqbuEKuWDlHVGc/2IWuX5byt+Uobg5hh/xWP6+eK1aeU0QT22v8fCO1ewahF0d5LP0Ppcvh86dpTRlYqJEWl9fNo9cSbcJQezdK9dKn30Grsr3ZFnWVmNM3dQe05kIpZRrnTgh5e7CwuS2C4JvVBT07AnNGifw3Mnh7PKpSRGvs9Kztdnke3DwjS8MCpKgfF3PNSUl5bbAykwIeoxtgZUBOBgDdOlC3jVLsP13jFIvPIS/FY+NJIiLZXvfr0lMlGM4I990SvCLiIzFcKWoREh4hN3OcSuZfp9TpsjqqeLFZZ94aChRbwxjUNBK7ukdRHw8LFoECxa4Lvjeig5BK6VcJzISWraEffuuJElIB3v22kJCoFcvOH4cPntiM698/xY88gg//+9t5v+4ljt2b+GvKnV5yL80HdJ5zHSlqixaFN+nu8CMKZiL8VjJSdTYMImlt0VSYsZIu+WbTut35ao6uFfL9Pts2FBWyI8di8mVm6/X3cnL44JISoL33pPtv+403JwaDcBKKde4cAHatIFdu2DhQtlylA72KgV47pxsMV74w0V6lAulx4aW1KvXAN7cRgi3yTnyl2dZUHkA1mbgHOlOVXkp/aEVGoq55152f7Ga+38YhdU2hD4ln+OzTg9irgsiGck3favflTsUlUh3Xu2UTBlnz8L06VCpEkyZQnQ0vPS0JEoLDpaOcfnyTml6lukQtFLK+eLjpbJRWBh8+60MJaaTPYZmT52SREiF5n3J6dylmXD4IeoVOSQP1qqV5XN0qBXI8I7VCCwQgAUEFgi4+ZzmpWFsq1lTqsx9j/gd+9h+5yPceeQgR6Y35e5Nx3kpbA61I/ZkON/0rd7HzYK5M4tK9GtZkQCfa7ea3fA+U/I4jxkDM2bI2DKyUr1OHZg1S3q9K1Z4TvAF7QErpVzBZoNcueDLLyX1UAZktdcWESE7Vrrte4s3k4djXUDmeI8fvzxZaI+eYWZTVeavWop7989k/eoEWnfZwHe/vow/ccR7+7Lhy7k0zsAxb/U+3KGoxC3zaicmynhyShEFmw2zcxdfRLTjtdckA9nKlalP0bs7DcBKKec5eFCCXcmSkqEoExkpslIK8OBBScDweMQn9E8afuWBpCRZzXxpMZU7lBts0NiHH178Da934rEA38R4ak78HLo+KOmb0uFW78NdikqkebHy6KPw22/g7Q3GYHx9GbQimA9/lYGTGTOgaFGnNtdudAhaKeUcGzdC/frQrZvczmQ6qHQNWaZi1y6ZZo6Kgp69bDIGHRCQ6irnzJ7D3rybBeMV4IfxspGMjfybVvBf2XtI3vZ7ul6fnvfRoVYg6wY05dCI1qwb0NT99hk//7xM8K5Zw5H/DeORAisZuSaIkSNhyRLPDb6g+4CVUs7w44+SLKFECUmSXzFrgSyjq6C3bIEHWyZTyfYnE1dVpGpVZO/xhg03rf/q6v2xl13aa3yudjCT3jvOU2EvsrlkR+ptnsDtt9/65W7zPtLLGJg4UYac+/S5fNdnn0H//rLr6PvvnVo1MEvS2gesAVgp5TjGwCefSAqie++V1c5O7rKsWQPtWyfypXmOjvyI166dUKaMU9tgL8bAN2PO8MZAH8ibl/nvhtMwYJvspXZSEXmHiomBF1+UceX27WH+fE6fsWj1cCxb1gaQ687jVO28n4EdKrj3RcRV0grAOgeslJvwuJ5KesTFyZaRjh1lGDHAeXOoIJ3tzh0v8oNvZ5pHz5d8hKVLO7UN9mRZ0K1PIeq2gCeeANsrL2LYCJaF5e8vq5E8NQj/8INUMTp+XJY0DxrEuvUWHR5J5PRJPwo+sIu8tf/mv4tkatuZO9I5YKXcgDtkJLKr6GhJDxgQIEO8c+Y4PfiuGRnG5jbvsc7WiObn58sWlnfecU6hVwerUkWm1M/UbYnBwjIGExsrQ7eeaPlyWWx1/Dj4+pLcrDkjPvKicWM4H59AsSfXka/O35f/dPbOCOYqGoCVcgPOSDvoNIsWQfXqkqUIoHBhp9ff+2NyGHUHNOMd8x53x2yGt9+G3r2d2gZHCwiAh8a2ItnXn0S8MFjwzTckzf3R1U1Lv1On5PuWLZcvjExSEl93D2XgQBk4KdZtDX63n7vhpVdvsQoJj6DhiFWUG7CEhiNWecyFqwZgpRwpJga2boXdu+X2qVMyRPjII/Daa/DxxzBnDt6HDqb6cmdmJMqygwehbVvJfh8QIPkdXWDfPljQJxRf4iUseXlB7twuaYvDBQXhHbqS2Lfe5+2Gq/kfX9B8bDsOHkT+zaVUCHI3SUnw6acyF798ucxf+/uT7GUjLsmX6X8HM3GiJL4qWSz1LVcpW6k8efRI54CVsrchQyQ5/M6dEpSMkd7gtGnSGyxYUPbE/PKLpGMEurTowfCCHbnt/GnmzezHplJ3s65MTQ5Uv9elbyXdQkKkKo23N4weLb3NdO5Vtaf//oPWrZL4MGk3Xn4+kMjNCylkF0FB5A0K4kMDM2Y04vveUL/aBf72bkZAyUJYEyfK3+Imq72dbvduyQG6YYOkIq1ShbjCgXzzxEoOTQ/lQKlgxi4KokYNefqtkoW4Qz7rzNJV0ErZQ0SEFJEHqFZNCoPffTdUrSpftWrdWJLFGClGcOQIvxy9yGu/nSLfmf94Z9UU6v/zB0ViouR5FStKxqhGjZz7ntLjwgXpXR49Cm+9BR98cOX34GTR0RDc2PDcH6/QM3G87FnJn989go4T/fMP9OgBASsX8VXAKxSLPSx7nY0BPz/XLdQKC4P334dly+TvMnYsdO7M8hUWvXrBn3/CU0/B+PGQN++1L01rgWK5AUtILYpZwKERrR3+tm5FV0Er5Si7d0vgWb5cPkFKlJAEtemZ87Qs6Q0XLEirahB3WwSjlvrySvv+BObzY2gFQ9OIHfKBmbLhc+ZM+PxzaN5cqgjVr++SniYHD0LfvnDmDKxeLe97+nTnt+OSxER47DFoHj6SnmY8vPEGjBjhsva4UunSEuMmTGjL3f2aMt+7NfclrsYC2Vt7VcYvp7h4UYLvQw/JqngvL/j6a47Wak3fzjLMfOed0ubmzVM/RFqZstwha1mmGWOc9lWnTh2jVLbwzz/GdO9ujJeXMfnyGfP++8ZERzv+vPPmGdOggTE2mzFgTN68xrRvb0xsrOPPffiwMR99ZEz9+nLuXLmMGTHCmIQEx587DcnJxjz3nDFP8bW0q0sXY5KSXNomd7FvnzHP3b3exOBvkvAyiX4Bxqxfb8wnnxjz/feO/dsdPmzMwIHGFC1qTLt2l//NJtts5rfWH5q8eY3x8zPmvfey9s93/rYjptKgn02Z/osvf1Ua9LOZv+2I/d5LFgBbzE1iog5BK5VRJ09C2bLS7Xr5ZSnIXqSIc9sQGQm//gpLl8Lff8t8MsCrr0ov5777oHJlGb7O7AIkY6Q3X6qUzF1PnSpzd3XqyPLUbt0kp/N1nL2fedgweO/dRP4pVo8SVQvL5l9fX4edz52k53edmAgLB4axb3IoC88Fkyv4XhYcqU2evy79bXv3llXrW7dmfrj+UrYugoNl4eG4cZJ0BWSet3Vr6NMHEx/PxWRfmpiV5G8ZxLhxcMcdWf0tuPcees2EpZQ9nDx5JYvTV1/JeJm7ZVTq1g3mz5cJ0RSdO0vJP5DxvhIlJHDGxMD58/KeKlSQ4cHJk+W+48cl0e6hQ/Jh2quXFNA9ezbN93x9/VmQBTM3LcWXRdOnQ/fu8ranf3oWy+Yl84s5QEZ/13FxMGkSDB8OJ/5L5p1aS3jd+oR820LlCV5eMke8aJH8jcuXvzKVcnWAvTpAGwM//ywXZImJcuFToQIcOyY5nHv2hDJl+OMPWDIojPOLQtlZJJhuE4Po1ClbbMm+JQ3ASmXVggXQtatk68lA7VqXiI+H/fth717Ys0d6Oc88I1tScue+cWtKnz6yJSQm5kpvOSBAPmw7dpSUgOlMH9lwxKpU5+MCCwSwbkDTrL2v6yxbBr0eOsSnt4+kxe7P8M3nf+sXZSOZ/V1fuAATJsDIkXD6NCwq+wqt/x6PhZHFWl26SNayPHlkQWGxYnIxlpwsAfbdd2Xe//Bh+YqJuXJwm0221w0bxsnz/nz7rVwk/f67LFV4+WVJcnX9IqvsTBdhKZUV06bBc89BvXry5e58fa+svr6at7csFNuzR3oouXNDvnzSYwEJuqdOyadjJodw7VFHNz127oQv2v/MJrpRIOoi1n+vQ7477XoOZ8jK0Glmf9e5c0tq7p49ZSHy2B7o3+MAACAASURBVBFdaMoUfIkn2fJle+lO3D3mfnL99Qds307STz9ju3TRlhh3kf9WrCXw7AmZ4mjVSnrBEyZAUhLG15fVhTry6eP+/PSTdIrr1Lm84NnpMzXuTgOwUmkZNUqKgbdoIRV98uRxdYsyz7JkaPFmQ8iWJXO9WeCMFaknTsCwB0KZG9cGL5KxEvzkwuFOzwrA1w8hpySQgPTlOM7q7zpvXkkQ1qtXEHMHrCTul1BmRQSz9oMgLAtq1ICSlaOxGsxlztqeeCcnkmDz5rWirWjyyqPUKX47Z8/KQnhf22P4rAtlwu5glr4VxO23y8DK00/feB2ortAhaKVuZuXKS9XbH5fqLDlkYU9WOHoO+OJFaNbUMCGsJtXNH3KnzSYrsQYOzPLxnSmrw/WO+F3Hxkp+jDVrZJR59dokkhNt1CeMJraV/JrUjA2kvkjLz09mK55+Wq5Xva/q3rnzIilH0yFopTKjaVOYOxceflg+5NUtpXyoOuLD1hhZ13N0/SEqBRyCRJ8r85IemOkqq8P1jvhdBwRAkybyBVD2jWXEHc/P3n8KsyumK15+CRQM2IXNP4FpPWtSsCAUKsTl76ldo2a1p5+d3TIAW5Y1FWgDnDDGVL10XyFgNlAW+Bt4zBhz1nHNVMpJ4uLgpZfg9dclk9Ujj7i6RR4nraQJWTFihKwNeu+98vi+sF+Sgaxe7bGZruwxXO+o33WKwCJ+RHifxb/ktR/vgQUCaNMmfcfw5FSRjpaeYgzTgeuXfQ4AVhpj7gRWXrqtlGc7f16y9UybJuNwym3Mmwc/vbWW76t/yDuDjGQGa9BAhp09MPiC5DgO8Ll2ZOXqHMfuwB5tdNbCPE90ywBsjFkDnLnu7vbA15d+/hroYOd2KeVcycmyTHPNGkn3+Oyzrm6RumTbNhjc9S8WeT/Mo3EzsC5E3/pFHqBDrUCGd6xGYIEALKRX6aj90plljzberEfvdqkiz5yRhZZOlK5FWJZllQUWXzUEHWmMKXDV42eNMQVv8toXgBcASpcuXefw4cN2aLZSdvbuu7KQZ/x4GYJWbuHoUXigzlkWnapPuXyn8dq4wT6pk5TTODs5S4acPCkT2N7ekj/8449lD70dV9SntQjL4fWAjTGTjTF1jTF1i6ZzM79STpWUBJs2SQmZF190dWvUJTEx0LFNPBNPdqKc9TdeIfPTHXw9tUB7duR2Pf1jx2TfcrNmMpWxZo3c36sXbN7s1Au8zK6C/s+yrOLGmGOWZRUHTtizUUo5lc0mmX4SE3NGbjwPkJwMH7QJ49nwqTTyWoPX1OnpLseoq27dj6MXi6XL0aOyR2rlSllSX7GiVDIrX14eL1fuxpKhDpbZHvBC4OlLPz8NLLBPc5RyoshISSJ87JgEYT8/V7dIIZ+NYzuH8favzXjWmoaXr8+VbF3pkNaqW5UDJSfL98KFZaPzu+9KKrU9e2TaqWxZlzXtlgHYsqzvgDCgomVZRyzLehYYATS3LOtPoPml20p5juRkePJJ+O47OHDA1a1RVxk1Cm6f8xn+XMTLJEnu6tDQdL9eV90qQEa0Pv8cataUBNh+frB2LQwZIlsM3WC065ZD0MaYzjd5qJmd26KU8wwZIsPO48ZJ6T6Vbo7MajRjBvza/2cW8SOWF2DZMpxow2kF2i9evDJqsmULhIdLxQEfH2mzj4+U4fPxse951a2tXCmlOXftkqplkZGSBNsNgu7VNBOWynnmz5ehp+7ddcVzBjlyfvWXX+DLHmEs9+qEV/XqWKM+kkUxGUy00a9lxVRX3dplf+2JE3LhtnChlGPavVtyay9aBEOH3vj8c+ckAI8eLSWBOnSQAgaenFPcnV24IPO8P/4o87khIdCundsF3hSaC1rlLMnJUtHI21uyKPnbr4RdTsh366hyg5s2wUuNd7EyoRF5yhTGtv43KYOXSXb/W/zxh6yQDwuTSepSpeSDvV8/CcDR0RAVJaUgExLkKz5eKhp4eUlwHjtW6v/5+0uvrHNn+VL2Ywy0bStJWvr2tev/78zSesBKXS0yUhZjFC9ut0O69V7Hq2Q1MJUbsITUPjEs4NCI1plq0/790LAhvJU4lN6+E7FtWO/01aipWrdOPtDvuw/++0+Gk9u2lYoDNWpkvFeVmAi//SYjMCEhULu2/AywdavcdtOemluLjJTUsUOGyIWRMW71e3TpPmCl3IIxkmIyLg4KFLBr8AXPWHmbcpEQERmL4crwcUb2yNo7q9GxY9CypXxett30DrY/fnd98D17Fl54QQLv4MFyX7FiMs87eLAs6snMB7y3twynjxkDf/8tleoBduyQUZlGjSToexiX7rnevFkuXGbMuPK7c6PgeysagFXOMHOmJNqYMcMhh/eElbf2uEiwZ/7iqCjo2CKaMf92ZNXnu7jjTitLw85ZZoysiq9UCaZOlcxICxc65lyWBfnzy8+VK8MXX0hxifvuk3ni3bsdc147s8dFXaYYA59+KkMnSUmSTOOJJxx7TgfQAKyyvyNH4JVX5MPNQTmePSHfrT0uEuyV1SguDh5pF8/QXR1pYxZSNfehDL3eIRYtgi5dZF/oli2yHyp3bsef19tbetx//gkffAC//gr160txEDfnspGfTz+VOd6HHpLV5x5akENXQavszRgJuomJMuTnoLq+Dl15ayf22p6T1axGUVHwdvA6Rv/+IjXYAVOnke7advYWHy+9zZo1pQ3ffQePPuqa+s+5c0tmphdekKHVvHnl3+8XX8iFQUqP2Y04feQnMVEuWJ57DvLlk//bHjTkfD3tAavsbdIk2S4yenSGsilllNvlu02FO5S/O34c+tZZzWe/N5bg6+MjKQFdYf9+mT9s2lSuCry8ZBjTFcH3akWKwIMPys/btslWucqVHTccngVOG/kxRkYk6teXBZT58kkQ9uDgCxqAVXbXuLEMVf3vfw4/VYdagawb0JRDI1qzbkBTtwq+4PqLhL/+kim7MofXYOPSSEFycoayXNnN5s3SmP/+k3UBbti7BKBOHWlr0aKy+rpLFzh1ytWtuswpF3VRUdCxI7z5pizQS0y037FdTLchqezJzbYi5HTbtsHTLY4Rm+zHgo/2cXfvZjL86+srWYucOYe3dCl06gS33SY/27H0nMPEx8Pw4fD++5JGMTzcbf59O3T/+x9/yN/q77+lB/zqq27zvtNL9wGrnOfjj6XnMH26W2zGz8lWroR+7fax4GJLCgVVIvfaXyShRWhohrNc2UWPHhLAfv5ZytF5kj/+kGQeTZpIso+zZ+VCIjsyRhJqHD4Mc+Z4bMpYDcAqZ9m9W+b2HnwQ5s1L1xVzTshi5Qpz5sDnXTew0LQhX0Ebtl9+kmFVV4iKkqHmhIQr84ie7MMP5UJz7FgZmvawnuFNxcXJMHOePHDoEOTK5drtaVmkiThUzpGQILlg8+aV1aPpDL4u2cuYzY0fD988vpjlyU3JV7qAZLhyRfBNTpZMSXXrwpkzsvDL04MvwMMPw113SVWv9u1lu52n+/tvSUjy3HNyu1w5jw6+t6IBWGUvI0bIHs6JE9P9H9cTslh5kvPn4eNHwjj28vtM9e+Jb827Jfg6cBX6TcXHS83nTz6REZECBZzfBkepXFlSW44eDStWyNzwjz+6ulWZ99Ol0ZE//8wxObJ1H7DKPi5cgAkTZDjukUfS/TJH7mVMToZ9+2Q6evNmiFkZRsVjoey+LZjj5YIoUkQWuBYpwuWfb79dRtA9beraGElx/O3zq/j6dBv8rHi88MH66GvXzFNGR8sCnmXLZLh2wIDsM0ybwmaT3n2HDlIsonRpV7co46KiZKfC1KlQvbpcRNxxh6tb5RQagFX2kTu3LK7x9c3Qy+xZP/bECcmKt3mzVPjZuvVKQqP2fr8wN749NpNAUpQP42JHEJrUiF/P3cXRaBkSrU8YwYQy0CcY231BNG0q623q1Uvf23LVXPbhw5JsjEULmenVjQDisIyBBOQX0cwF5cP79pUVYFOnSunJ7KxCBbnQSNGnj4wAvfGG+9cjjo2FxYth4EB4913Pu/LMAl2EpbKHzZtl+Mor47Mq9qhkdP48jBwJv40Ko0H8Kv6xlad86SQa593GxdYPU6ZrIyp/9j+8vpp844sXLeJi8zbEfvgp+Yf2BSDZ8uaLYu/y/vHnOc7t5M4ti0BTAnLt2jfmi3BFRabERKktMHnQP4xO6E3bpAWYcuWxjkbIg67YZpTi9GnYuFHSFeYkSUmSUOSHH6Rq01dfyfy3OzlzRtZoDBgg/2fPn5d1G9mQroJW2duuXVCrFgwaJFfQmZDZnmNSkux0GjQIah1fwmKrHZZJ5vJAp7+/zD+++KJsu2nVSgKTj49EruLF4d57ZYj25ZdlCP26/5MrPt7OgkPV2bP0H/b9aVGSI7TOFUpy42Aqdw+iRQtZ3OuoWr03s3Gj5Dcpuz2E72xP4udr8HpviPS+tmzJ0jajTPfkU6pede0Kfn4ZPm+2Mn8+9OolyUb69pX/G+4Q5BYsgJ494eRJmcOuX9/VLXIoDcAq+0pKkoxGBw7I9qOiRZ126lWr5HNt7/Y4agf5M7v6B5SaNEge9PKSgPrxx5K7NkVa+1/DwmSoNiVBxeefy/xY795yjJdfhvHjMVgYIAEfWrOY1d7Nuf9+CGc3ARVO4FPowjWHzUqt3uvFxcn73jIujIs/r2JHkab0GlqMFqv6Y40eLcXpsyhLPfkhQ+C992DyZHj++Sy3xeNFRkoGqZkzpdrS7be7LknNqVPyb/m776RnPm2aXDhncxqAVfY1Zoz0uGbNksVXTrBvH/TrB+sXnWJk3g94wvdHch3ajbVzx7UBNDNDr2kF6D175MRLlly+KzF3ft55+SyLl1ic3XmEo5SgYd5faZZnKWGlq7GjclnK3JFI2NtNMv1+T52SUy5cCKt/ieWZmHGMYKCkkwwIwMrg+7xV7zbTPfmJEyVvcvfuMGVK9ltwlRXHjsloizEyJF+zplw9OvGClRYtpNLTO+/I0HMG12p4Kg3AKns6dAiqVpVgtXixwz9w4+Lgi6fDiJq7jFK2o3S1fY9vQjRWjx6y/alwYcdneLq6l+ztLR9mb78NxhB9W0niI6PJl3gBi2QS8KUFS9ngez+1alrUqcPlr0qV5HCJiTKIkJh47c9RUbKmZ8ECqXPeLHkZr/l/QdPEpfglxmCQnjU2GwwbJgto0iE9vdtyA5aQ2qdSmj35uXPh8celotG8edeOOqgrYmIkE9icORAQIEPBb7whwdmejJHh5S++gHHjoGBBuSAtWlRWOucgGoBV9rR1qwwzLlgApUo59FTR0TCwSRgfbWmKP3ESfO6/Xz5gKld26LlvkFqQT0yEb78l6v0R5Ptzz+U56EOl72H8oxvZvjmeehs/Z+PFmvgQTy1+J5RgNiCvt0gmFzH4Eo8fFwnmV55kFovK9abYUy157twnlJz7CVb79pI7+a23MtXTT0/vNsM94AsXZBXwHXfIVUOuXOlqS462d69szfr2W7lY+flnWd2XVTExcsxx42D7dtl3HRIiRVFyKA3AKvtywnzWmTMyatds03CGWe/glZwkc7zvv5/unp/TXN1DttnkAqF7d9i5E6pVu/y0lP/165oPYfODgyl9eC2dxtx/4/F8fSXY164tP6f8rjPZ009P7zZTc8C7d0svrmDBdLdFIWsnxoyREZxcuSQoL1kiC6OCguR7yZKpvzYhQYa24+Pl4ufMGbk4O3NGermvvCLTQjn8giitAKzjNMrzHD8uC5Teftvh/7mPHYP2zaLpu+8FKvRpg9cXvld6fsHBDj13pgQFSY/0+uBYtaqsOu3fH6ZNkz26lsV9d57gvteAiPJQcpS8r1WrYNEiySKSlCTHuj7IBgVlaog9PXuuU4LsLVdB//WXTD28+ipUqZLhtihk5GDs2Cu3ixaVi6zx42X1Psh8xe7dcn/v3rBhg6S9PH5cLoCbNpV/c4UKyXqM4GDZM6dz8LekPWDleR59VALE9u0OLeZ+6BA83vQkE/5pTR22Ys2YAeXLu66Kjz1cv9I6teHj9Dwnk+y2V/nUKemdRUVJ7z4b5wt2ifh4+f+1YQOcOycXuyA5mo8cgcBA6RmXLCn/B+9PZfREAToErbKTkBBJQv/BBzIP6SC7d8OzTQ4y81RLyvpEYJs7G9q2ddj5nCo9w8cOXEyW5WxdFy/CAw9I8pXQ0Gy/j1R5Ng3AKnuIjJShxttukw9fB6XY27wZXm6xn8XnGlEwbyLePy/2zN5udmSMVLv65hv4/ntZ+ayUG9NyhCp7GDhQsvpMmeKw4BsaKlNaUflL49+2Bd4b1mnwdSdbtsie76FDNfgqj6eLsJTn6N1bMuc4qKbs5s/DONpnPE+UfIoh61qSN/Abh5xHZUG9ejJEkQMyKKnsT4eglftLTHR4YoVj88Io3KkxviRgvL2x1qzRnq872bhRiivktMIKyuPpELTybL16SYHu5GSHHD42FtY/Pw0fEgBki05oqEPOpTLh8GFo1062uMTHu7o1StmNBmDl3pYtk8T6pUplqtTgrRgDb3c5RLMzs2Xfos3mvnt8c6Jz5yS95MWLkvEsh+QPVjmDzgEr9xUVJfsOK1WSRTcOMH48tAh5ET9/L6yvv5fMQJ66xze7SUyUurZ79sAvvzg/5adSDqYBWLmv11+HiAhYv17q6trZ2rXw2mvQpfk0WrxzEBo1tPs5VBbMnSs5iidNkn2/SmUzGoCVezp5UpJuvPmmFKy3s6NHYWK7n6lQtgVj5xbHK7+dq8GorHviCalfa48iAUq5IQ3Ayj0VLQq7dkk1FTuLj4dJjb/l28iuHHvpc/Lnf9nu53CULGeR8gTz58u0Q+XKGnxVtqaLsJT7+eUXKQJQrBj4+dn98J923sTAv3pwsvL9FB/8gt2P7ygpeZQjImMxQERkLAPn7SAkPMLVTbOfVaskwYYD04wq5S40ACv3sngxPPggfPmlQw4/+5MInprXgZj8xSm65kePWlU7aum+a4oYAMQmJDFq6T4XtcjOtm+HDh3grrtg2jRXt0Yph9MhaOU+zpyBF16QurXdu9v98HumrOee17tSyCsS79BNUKSI3c/hSEdTKeOX1v0e5fBhufDKn19GQBww9aCUu9EArNzHq6/K4qvFi+0+9By7Koyyzz+ALxfx8vbGij1v1+M7Q3pq6XqsYcMkI8pvv928ALxS2YwOQSv3MH8+zJwpdUdr17b74df0X4yPicdGMlZKkXkP069lRQJ8bNfcF+Bjo19Lx9VEdppx42D1arj7ble3RCmn0R6wcg/Fi0P79g5ZfLNh3lHqb/lcMmlZeGymq5TVzo5eBe20ldZJSdLzffVVKFgQqle3/zmUcmNajEG5ljGSAtJBYi4YthR9kHpxa7CmTcP/6EHNdJWGlJXWVy/2CvCxMbxjNfsGYWMkx/fEifD119Ctm/2OrZQbSasYg/aAlWu99prkXx492iGBeEmbiTwau5T9fSZw19NaP/ZW0lppbbcAbAy88YYE3zff1OCrciydA1au8+OPMGaMDEU6IPhumbWP1qFvsLt0K+76pKfdj58dOXyldWIiPPssfPIJvPwyDB9un+Mq5YE0ACvX+Osv6NED7rkHPvrI7oePiYEv++3ntHcxyqyc6tBh7uzkZiuq7bbS+swZSbYxeDCMHeuQCldKeQodglbOFxcHjz4qQ89z5jgkGcbbb8PkY23psrwVpe7wscsxc0IayH4tK6Y6B5zlldbR0RAQALfdJgk38ufPYkuV8nwagJXzbd0K+/fD7NlQpoxdDnl1cKy//xR551+k10vP0/gB+wXfqwNTShpIIFsFYYestD51Ch56COrWhQkTNPgqdYkGYOV8DRvCoUPSG7KDq4OjX0w87y14H38rjo0dgoG77HIOpyxOchMdagXa7z0dOQItWsDBg/DOO/Y5plLZhE7AKOfZswe++UZ+tlPwhWuDY7/vf+SO5AP0CR7IZ1uP2O0c2ToNpKPs3y8XW0eOwNKl0Latq1uklFvRHrByjgsXZN73xAn5ILZjrt+UIPj0sl947uQs5hRsy7Z7ymDZMThm6zSQjhAfD61aSXrJ0FCHZDdTytNpD1g5njHQsyfs3g2zZtk90X6JAgHce3AHg8PHYYC255dSO2KPXYOj3dJAJiTA8eNw/rxsv8puYmLkffn6wuTJsHatBl+lbkIDsHIsYyTZwsyZMHQoNG9u91P0a1mRqptOYfDCAnySErkvYtc1wTEkPIKGI1ZRbsASGo5YleEauh1qBTK8YzUCCwRgAYEFAjKWHeqHH6BxY1mAVLw45MsH3t4SsEBSMlapAvXqye9owABYuNCzgvTSpVC1qiTYAHjgAaiYDfJUK+UgOgStHGvzZsly1auX7A1ygMq2Inx8pCOvW+PxJZ5Ebx/qdHuYxpeCo71WMKdrcdKhQ7BsGYSFwfr1sGQJ3HmnbMOJi5Nyi3feKUOzFy7I1hyQCkB33y33nTghiSq++kqqQwFMmSJJLBo0kOe50/7ZEyegb18Z3ahUCWrWdHWLlPIImgtaOd6qVZJ/2QFBw0Rf4N/bajOWV3lrbi0K/RF6Q67nhiNWpTp/G1gggHUDmtqnIadOyQXGlCnSay1aVILl++9LrzCjYmMlmFepIrcbNpSADtJ7btIEnnhCvlxp/nx47jkZUn/rLRg40O6lJJXyZJoLWjnfDz/ISuf774emdgpyqdj7yCAqx+7n/jeqUah1ELS+sciCU1YwGyOpNV96CV55Be6444bsWxlK5BEQcCX4gtTJPXBAetZr18pwb5EiEoCNkUDfuLFcePjYZ+9zmhIS5DwFCkg7J026tr1KqVvSHrCyv5QtJ8HB8rOD0kBGLd1A3lYNmF/sRR4+Ov6mHWyH9IDj4yXoLF0KixbJe4yOhjx5Un263asMGSO95Fy5ZI9txYoyRJ0/v8whP/AAtG5tv+L2SUmwaZO818WLJdBPmnSlLZrqU6lUpdUDdqOJJJUtrF8PHTvKPOWcOY77YL54kQtPPMsRSlJx3vA0R7ftWsg+OVkyeFWpAr17yyKqs2flsZsEX0g7kUemWJYEX4Dy5eH0aemBP/KI/A169oRt2+Tx8HCZo509W4a1M3rRPWCALBxr0EDydhcuDPfee21blFIZpkPQyn7++EN6XYGB8Msvdt9udLWdk9ZRMXI/Mzou4NkG+dJ8rt3SK/77rwS4TZukePzPP0PLlukKQA4fBs+XTy58OnaUAHt1prFdu2Rl8qefyu0iRWRe+scfoVAhqcf77bdX3odlyTzu/PnyszHSo27bVvb2FixonzYrlcNpAFb2M2UK5M4Ny5dDsWLpfllGixwkJEDnL5uSp8QBVswona5z2CW9YoEC8v6+/hq6dpViEunk1EQeliW94hRPPgmPPw47dsjFw6ZNkqUqxcWLEBUlPxsjXzYbHDsGJUrAyJH2b6NSKmtzwJZlvQY8BxhgB9DdGBN3s+frHHA2lZQkH9jJyfKhHZj+QJfhudGkJL59eT1dv2jEggXQrp093sAtrFwJ9etL8M3kfKfd54CVUh7BIXPAlmUFAr2BusaYqoANcPGeCOV0ixZBjRrw33+yzSgDwRcyPjd6ZshYunxxP/3vD3NO8B03TooJDBsmtzM535nlRB5KqWwnq0PQ3kCAZVkJQC7gaNabpDxCcrJsfRk8WFINJiRk6jAZmRs1Bw6S+8O3WWJrS69v6mfqfOmWnAz9+klCjPbt7VLJx65VhpRSHi/TPWBjTAQwGvgHOAZEGWOW2athyo2dOyeLfQYPhqeekj2qmdzucrM50BvuX7+eCw2bk5hsEfHWBEqVduDK29hYeOwxCb69e8tipdy5HXc+pVSOlJUh6IJAe6AcUALIbVnWk6k87wXLsrZYlrXlZEpaPeXZ+veXvaCffSYLkgIyv5AoXVuEwsIwwcHk+e8gvlYCPZr/m+nzpct//8lFxWefwZgxGVpspZRS6ZWVfcAPAIeMMSeNMQnAPKDB9U8yxkw2xtQ1xtQtWrRoFk6nXC4xUb5/8IGkl3z11SzvAU3X3GhoKObSPLG3lYz3b6FZOudNpawELltWVgm/+qpjzqOUUmRtDvgfoL5lWbmAWKAZoEucs6OU+d7ly2VFcKFCkmLSTm41NzorrjoP44cP8STZvNlY/G4a2+3sl5w9KyudH3lELjDypb23WCmlsiorc8AbgR+AbcgWJC9gsp3apdyBMZJsom5dme8tV06CsRNt/GAC6z49QAv/n/ikwVN0efx9eh7wzXA5wTTFx8uc9t9/S6IJpZRyAs0FrVJ39Kgk+l+7VgLvsGHQpYtz0w4eO0ZU6UrsSKzGI+0/wb/SlTUEdqtkZAw8+yxMmwbffCNJK5RSyk60GpJKv3PnZPi1aFEJTuPHS7k5X1/ntsMYLjz9Ir6J8bxUZiR+Fa9dwGe3FI4jRkjwffddDb5KKafSAKzE/v0yzLxmDfz5pyT6X7vWZc0x331P7uUL6GcbyZmH4vC+ruNttxSOlSpJD3jIEPscTyml0kkDcE529iyEhMiWogULJAH/a685fZ73BufPE/+/VwjnXiL7PEvegM3EXpXnI9OVjK4WEyMXGQ8/LF9KKeVkWo4wJ0lOlkT8e/fK7b/+gh49YMMGSThx8KCsdk6jrJ4znEnIy9O2mXxSdRpfjCxs/xSOf/8Nd90F331nryYrpVSGaQ84u9u0CfbtgxUrZEXzyZNSK3biRKhTR2rF1qjhPjVdY2Lo3z8XP0S3YutMyYFh1xSOUVHQpg1cuAC1atnnmEoplQkagD3dn39K7deDB+HAAfkeGAhjx8rjjz8uPb5ChWSLzUMPSQ1bkOIJNWu6rOk3OHWKi1VqwsnB9O33PDVq2Pn4CQmSYnLfPli6VOZ/lVLKRTQAu7vFiyEsDI4fv/KVLx/8+qs8/vzzsHq1/OzrK3Vgr8449t13kD+/DLm6c0rFrfGGNgAAIABJREFUsDCS//ciXieP82+J+nw22AHneP11WLYMvvoKmtphC5NSSmWBBmB3cOoUbN9+5evQIQmqlgVz5kgQLVZMvm6/He6448prhw+Xnl25ctLz9bpuWr++g6sG2UNYGAQH4xUfTyLeDHkj2v61D4yBwoXhjTdk1bNSSrmYBmBXe+stCaIpSpSQOdkLF2Qx1Oefyz7Vm/Veg4Kc005HWrAAEx+PBdgsQ/24UMDO78uyZJuVUkq5CV0F7WzbtsmwcXi43H7iCfjoI8mzfOIERETATz9dWYmcP797Dx3bQZJ/LgASseHl5wvBwfY7uDHQvbsMPSullBvRHrAzxMbKUPLEibBxo5Tva9RIVuFWry5fOdgHXu+ylvqMfnwrNV4Ntm+vftw4mD4d6tWDFi3sd1yllMoiDcCOlpQEd98t87qVKkl92W7doEABV7fM9X77jf2bIhk6tA2Pd2lBjVl2DpA7d0K/ftC6Nbz4on2PrZRSWaQB2FH+/RdKlpTh43fekRqzwcHus9/W1U6dIvnxJ/A+mZvSxVsyfryPfY8fFyfFI/LnhylT9PeulHI7Ogdsb8bA5Mmy7WfmTLmve3do0kSDQIpL87JJx0/ySML3TJnhY/8BgZkzYccOWcBWrJidD66UUlmnPWB7OncOXngBZs+G5s11zvFmxoyBxYvpy1iavl6LJk0ccI5nn4WKFWWuXSml3JAGYHvZuvVK1qkPP4T+/W/ck6vgwAHMm2+y1K8dq+98mc0f2Pn4p07B+fOyL1qDr1LKjWkAtpd//4X4eEmg0bChq1vjtky58oyp+hUjd7Zm2bcWfn72PLiR2sXr10tKThcXlVBKqbRoAM6KyEgJuO3bQ4cOkmM5wE51arOb9eth0SIWW+14Lbwbo0dDtWp2PsdXX0lZxY8/1uCrlHJ7ljHGaSerW7eu2bJli9PO51CxsRJwt26Ff/6RNIcqVaunL6TBc53wTkrkIn48V2kpM3Y1tu8I/cGDEtEbNJBCCzr8r5RyA5ZlbTXG1E3tMf2UyoykJOjaFX77TVbZavC9qZDwCP4bMxHvpERJNUkiNQvOZOH2CPudxBhZ/GazwdSpGnyVUh5Bh6Azyhjo1Qvmz5fVvI895uoWuVxIeASjlu7jaGQsJQoE0K9lxcv1e3/8+me+2BWKwSIRLxK8vFlfoSLzl+6zX43fixdl0dUjj0CpUvY5plJKOZgG4IxauhQmTZJVzr17u7o1LhcSHsHAeTuITUgCICIyloHzdgDQ4XYvRk3pzzlbXromzKLGbZvZ0aIw2wIrY0XG2q8R/v7w5Zf2O55SSjmBjtVlVMuWsGjRtRWMcrBRS/ddDr4pYhOSGLV0HxQtyvIKLWmZuJSfizdl2pP3sy2wMgAlCthhsZoxMGAAZJd1BUqpHEUDcHr9/LPkFrYsaNNGs1pdcjSVnmxAfBzxR45y5Lg3fY5+w55clbit4xa8fJLlcR8b/VpWzPrJf/wRRo6EVauyfiyllHIyDcDpsX49dOwIr7/u6pa4net7st5JiYxbOJJ53/anU+s4EuJsjP7yHKVLemEBgQUCGN6xWtbnf0+flrn4OnWgb9+sHUsppVxA54BvZc8e6fGWKnUlt7O6rF/LilfmgI3hw6XjaHZgM8PvGMuWnf4sWQKtWhXjVeycj7lvXzhzRur8eus/Y6WU59FPrrQcOSJzvn5+sviqaFFXt8jtpPRkRy3dR5cFX/DYjhV8V3Ugb+18hfHjoVUrB5z0119hxgwYNAhq1HDACZRSyvE0AKflvfck29WaNbLNRaWqQ61AOsyfBBvmcrBae7rs+IA+feCllxx0woYNYexY2furlFIeSjNhpSUhAfbudUDOxGwmLAyaNcPExhGLP+82XMnI1UHYbA44V3w8+Po64MBKKWV/mgkro/bvl/lFHx8Nvmk5fx4GDoTlyzEX47Ew+BLPBw+EOib4/vYbVKgA27c74OBKKeVcGoCvFxcHDz8MDz4o+0xV6v76C+rXh1GjOBmXlzjjSwI2vPx98WsZbP/zxcVJpSObTYKwUkp5OJ0Dvt6gQbB7tyy60r2+qfvlF+jcGWw2doxeSvCwZtyTrz5TngqlRJdgCAqy/zmHDoV9++TvopWOlFLZgAbgq61eDZ98IquHWrRwdWvc09Sp0hOtXp2F3efz6JvlKFsWxv0URIkKDgi8AL//Dh99BM88o38XpVS2oQE4xblz8gFfoYJ82KsrwsIgNBSCg+G++zA9nmVU4Gf075Ob+++XuhSFCjnw/DNmQJEiUudXKaWyCQ3AKS5ehKpV4a23IHduV7fGbtKqVJQuYWHQpImsPvb3J+GXlbyQ9CXTh8KTT8JXX8k2aYf6+GPo08fBUV4ppZxLA3CKokWlyEI2kmalovQE4R074Pnn5eIEMPHxfN09lOkHgxgyBN5918HT5AcPyqKrMmWgdGkHnkgppZxPV0GfPAmPPgqHD7u6JXaXZqWitBw/LnWOq1eHQ4fA2xtjsxGX7MvXh4OZMQMGD3Zw8E1OlimBRo1kP7ZSSmUzObsHbAz873+wZIl0566S5aFbN5BapaLr7w8Jj+CnyfO4Y/cWDt9Vg+YvPU6HuwpAeDgMGoTp8xq/TdnHmqGhrLaCGbYoiOBgJzR+0iRYu1YWffn4OOGESinlXDk7AM+cKSuIPvromoQbWR66dRMlCgQQkUoQTqlgFBIewdzPvmPqrLfwSUqANRZPxSdDn8502LuXLeE23ugEq1cHUblyEPPmQaVKTmj4P//Am2/CAw9IL1gppbKhnDsEfeIE9O4N9913Qzm7TA/dXickPIKGI1ZRbsASGo5YRUh4RJabnRH9WlYkwOfalFSXa/GuWUPhro8xZdZb+CUlXPqHYKhzMJxh3x+mazcb9erJlugJEyT5lFOCrzHQs6cMQU+erHuxlVLZVs7tAY8YAdHR8OWXXJ83MT1Dt7fiDr3oDrUCscXGsGXMNILCf6VQcjwX3nufJrUCYeFWSh07xJpytQk+uBUvk0yCzYefTndg28f3sttHFoT37w/58jmluSI+Xko/Dh+uBTCUUtlazg3AQ4ZA06apdutuNXSbHmn1olMCcLrmma/eg3t9hqnYWJm/Xr4cSpSAgABJD3n//XDgANxzD23PnKHt1a/5dzvwILRtS9f+eYiIjKXWv3uptfFflv7bkfV7GlG09nG2hhSnVKl0v1378fOT+V+llMrmcl4ATk6GpCTp1rVpk+pTrikyf8nlodt0ulUvOl095GXLoF076RXabHLB0KmTlOGLjoa8eW88wTvvSAAuVkxWMh88CCtWyPu22S7nt446ZxFETb5cHMfCA80JifPFv8wpyjZfz6cvlb0cfJ26GO2dd6B9e6ibauEQpZTKVnLeHPCsWVLE/ejRmz6lQ61AhnesRmCBACwgsEAAwztWy1DguVlvOeX+lB5y7Yg9vBQ2h7r/7qL8kf0s+2KuPDE5WYLvxYsSNBMTYdMmOHZMHs+TB1q2vDJHarPJ3qChQ688PnGi9PT9/MBmI9nHl9n//b+9Ow+vqjz3Pv59mCMIEVEwOCBiqXJsRVM1WjX2QG0dDsOpFBRFX1CsisNRBEplEKqiOLTWo0dFWz2KXlKU2SJDEBQVFFAQkJAjymBAGRUEkjzvH/cOhJCEkKy91147v8915cpmD2s9i5W97/0M676z6dAhllhqUFNY15yjT/+O5t0+pH3fT3n8llYH9NAHjf+MdVt34dn/JSEuc9kTJ8LIkTBtWvDbFhFJQjWrHvCOHfCTn1hSh/nzoVb8vn+U7uGC9aKLA/nJA6dw8eqFPDd+BHWK7DkO+Kz5KZzxTa69YMgQW6FdUGA1cGfOPHAYOlaHd1+N3BKPb95spYxXrIDv35lP0awcXt+YzQdkcdppFtuvvNJGrMsrHXjBQ7PKHIpvmZ7GewN/Fcj/EwD5+bYKPSPDvmSo3q+IpIiK6gHXrCHokSMtycSECXENvrB/GLm84duM9DQG5bxA3VjwLQIm/fRCXrqyL+OKN3L//VYWsZw54L2ZWWx4cSY7p+awqEk2s1/MYsUAC7qbNu1/Xr16WZx/fhbdBsLLV0KbNpU7hiAWox2S93DDDfbl6NVXFXxFpMaoOQF45Up4/HH7sD/nnITssnP7lgcOW7//PnTtB488Qv9L2/LcB90ZOflx6hQVsrd2Hcae14WePbIP2EbhOVmsOiqL3FzIfQJWrbJSvLm5lryrsDALsMDcrJmtKevUyX4X/7RqVX4vtyJBLEY7pLFjbdj5qafg9NOD266ISJKrOQH4qadslfCDDyZ2v/Pm2VzskiWwbBkcdRQsX07nK66Aobdz+3HH0+bzheSensnvb+pKpzNb8sUXtm5qxgyYPRu2bt2/uSZN4NRT7TvE1Vdbb7ZNG2jb1gJwkIJYjHZIV11li+J69gxumyIiEVBz5oALCiwA/vznidvnrFmWzcl7Wyx11102rFyq2lJ+vk3fFgfdr7+2+086CTp2tHTIbdtaoG3aNLG5KeK2CnrXLvjhh+C/NYiIJJGaPQe8ezfs3Gk9z0QGX4APP9x/u1YtCzax4Ou9Bd3hw62TDNbEX/3KEmB07AitW1ccbBNxidBBw+hBGTAAxo+3L0VNmgS/fRGRJJf6lyE98YStfP7mm8Tsb+dO6N3bxo6zs6FBA5uArVeP4ioGc+daid2OHeHLL+HPf4YFC2zh1LhxlonxlFMOHXwTdolQ0KZOhSeftOFnBV8RqaFSuwe8fj2MGGHDwC1axH9/q1dboowlS2xB0d13Wzc3tor5o9pZ3Hep5ddo0QL++lcrt9ugweHvqjKZtpJSfr4thPvZzxI/Hy8ikkRSOwAPGGBzv489Fv99TZwI111nQ81TpsBll9n9WVksTstiyBCYNMlGoR95BG65BY44ouq7S8glQkErvuRo+3abH6/KNw8RkRSRugF4/nwrNzh4sE2mxmsfOTmQnm4R9eyzbQy5VSsAtmyBP/wBXn/dnjJypBVgKiuD5OFKyCVCQfvhBxtXHz0a2rULuzUiIqFK3QA8aZKN8w4cGJ/tl85C1b+/rXCO9epycy3VdF4e/OlPNhqdnh7c7hNyiVDQGjWCyZPDboWISFJI3UVYDzwAixfbh348vPOOXUpTWGhB+Kij9gXfuXMtxeOmTXZZ0YgRwQZfCCZfdcJs2QLXX28rzpxTjV8REVKxB7x3L6xda7VkmzePzz4KC21xFdicb4kVzi+/DH362Cj0lCmVT/tYFXG7RChIu3dDly6WBez66/cNz4uI1HSp1wMeM8ayVixbFp/tew933gnvvmu/R46EmTMpOjeL++6zdVgXXAAffBDf4BsJRUW26GrOHPj73/d9SRERkVTrAe/YYSX5zjsvfnmFR4+Gv/3NJnVHjwZsJPqGq22xVe/e8N//HUxNgYTW4o2HwYMt1/ODD1reTBER2Se1AvCjj8LGjXZJUDzmGQsKbFy5WzcrE4hd1tq5syW9GjXK1mIFsevS5QyLE20A0QjC339v56FvX7scTEREDpA6AXjDBuuRXnUVnHtufPZRpw68/bbdrlWL3FzL8bFxI/zznzbVGZTIJtoo1qiRzfs2bKhFVyIiZajWHLBzLt05N845t8I5t9w5l3XoV8XJnDk25/jAA8Fve+lSq2C/ebOtdG7QgG3b7DKj77+36eAggy9ENNEGWE7N66+HH3+0NJN1Uuc7nohIkKr76fgX4G3v/e+cc/WAauR2qqbu3S258tFHB7vdtWvht7+14P7999C0KYWF0KOHZZ6cORMyy6xzUT2RTLSRl2ffSo44ArZtU6YrEZEKVLkH7JxrDFwEjAHw3u/x3m+t+FVx8vnn9jvo4PvOOxZdN2+2AgInnghYtaJp0+DEK1bQa+oULnhoVuBFEPpf2pa0urUPuC+pE2189519Udm71/5z4nUJmIhIiqjOEHRrYBPwonNukXPueedcw9JPcs7d5Jxb6JxbuGnTpmrsrhzz5llaw9dfD3a7778Pv/mNrbIqLLQqR1h2y4cfhvSz11DYdnXcKhFFKtHG9u3QqROsWQMTJsBPfxp2i0REkl51hqDrAGcB/bz3Hzrn/gIMBO4r+STv/bPAswCZmZm+Gvs7mPe27DgjA668MtBNM3myDTuDrX7OyeGj2ln06QONW2+h8SUHXmccjwVSkUi0AbYAbvlyy0Jy4YVht0ZEJBKq0wNeC6z13hdXnR+HBeTEGT/eMl7cf3/1SguV5corIS1tXy3fTe2y6dIFjjsOmly+AFf74O8SSb9AKmiffmpfgtq2tfnfq64Ku0UiIpFR5QDsvf8G+No5Vzwp+e/A54G0qjL27rVCC+3aQa9ewW23qAieeQbat7cVViNGsHvqTK74cxbbttmlrSdklD1wkNQLpILkvdVUbN/eer1gK55FRKTSqrsKuh/wSmwFdB5wQ/WbVElLl8K339qkbBUudSk3y9TTT8Ntt1lA6dEDf14WN/aCjz6yDvcZZ0D/gghWIgrKzp2W7HrsWEtI8p//GXaLREQiqVoB2Hu/GIjDRTiV0L69Vddp3PiwX1pelqmGX+XR8d57bTVv9+6AJdd6+WUb5S6+1rd4XjbSaSKr4quvLO3X4sV2vfXAgUqyISJSRdHOklDFYc+yskzt3r2H5nfcA/Xrw/PPg3O8/bZlUbzqKqvpW1JkFkgFaelS+L//s1rLl18edmtERCIt2gG4ispaLNVnwVv8bM0yeOUVyMhgyxYr5NOuHbz4Yg3u6G3YALNmwTXXwGWXWQAOurixiEgNVCMDcFlZpma3ziTD7eH6Hj0A6/lu3Gi1FxoedHVzDbBmjV3wPGaMLbr69a/hmGMUfEVEApJ69YAroWSWKeftWt+1Ga1Jf3QUOMecOfDcc3DXXXBWYi+sCt/69VZTsU0b+0/o1cuu8T3mmLBbJiKSUmpkD7jkIqpuk57ntB/y2fXsGDq1b8mPP1oFvVatYPjwcNuZUHv27C9iPH483HIL3HMPnHBCuO0SEUlRNTIAQ2wR1Sdvw/zXLOXkOa0AW9y7ciX8618pPvT83XeWbvP99y2dZ926NtebkQHr1gWf2ERERA5QI4egAZg7F2680eY3c3Jg/nyWLYOHHoKePW3KM2UUFcGqVfv/feut0KyZlVh89FFLanLxxZZyExR8RUQSoMb2gBk1yoIvwJ49FM3O4cbJWTRuDI89Fm7TKq2oCDZtsnnbDRvgoougUSPrvj//PHzzjf2sX28JNPLz4dhjrWzjiSfC+edbtae0GpLBS0QkidTMAFxQAB9/DLVq2fVF9erx5uZs5s+Hf/wjidYbFRXBsmWwZIkNC//+9zY5PXmyzdFu2LC/1wqwcCGcfbYF5aVLoUUL+MUv7PcZZ+yvz9u5cyiHIyIi+9XMAFynjgW1efNg5Uo2np7NDddm0aEDXHtt2I3DhosHDIB337W52mKnnWYBuEULuOQSaNnSfjIyrEpEcRnAnj3tR0REklbNC8BbtlgGrWOPha5d8R5u6mIdyWeeCSHhxhdf2MXGOTk2J9u7tw0jL15sFZmys+Hcc201cvGqsMxM66qLiEhk1awA7D387nc25zl5MmBX3EyYYFPCp5ySwLYsXQojRsAbb1i72rSx1dhgvdm8vAQ2RkREEq1mBeAJE+xSm7/9DYCtW6FfPzjzTPiv/0pwW2680YLwwIE2n3v88QlugIiIhKnmBODdu+Huuy25c9++gMW+/Hyr8VuFioaHZ9Ei62Y/+aSt8nrhBWjeHJo2jfOORUQkGaVsAC5d7/fZDTNol5cH06dDnTosWQLPPgu3325TqnHzySeWUmviRJt7XrzYLgM67bQ47lRERJJdSibiKK73u27rLjywfssP1Hr9dTZc/GsLfsCgQRYPhwyJUyO8tx5vZqatZh4+3OoXx/YvIiI1W0r2gEvX+/WuFp17jqZNgyKmALNnw7RpVuwnriPAixZBt27wP/9T5drFIiKSmlIyAJes95uxfSOb0xrzY90GfF5gHdMBA2zN0223xWHna9ZAYSG0bm2XCtWrV4OLCYuISHlScgg6Iz2WWtF7/jrxEV4bOwi8JyM9jXHjYMECuP/+OGRgnDPHhpyvvdYiff36Cr4iIlKmlAzAxfV+b39vLJnrljPvxDNJq1eHu37Vlj/+0RZCX3ddgDv0Hp56Cjp0sCIHL76owCsiIhVKySHozu1b0mzB+1zw3qt4oM8nEznj5p58/nFLcnNh0iSoXTugne3ebdWFxoyxzFX/+7/QuHFAGxcRkVSVkj1ggF+OH4MDHNCgqIDzvlzG8OFw4YVw+eUB7qiw0C4t+tOf4K23FHxFRKRSUrIHDMCOHQdUO3plXTb5+fDmmwGNDu/aZd3oI46A996z+V4REZFKSt0APG8ezJwJCxaw+WfZ3NE9i65dISsrgG0XFsI118D27VZ7V8FXREQOU+oF4G+/tQDZvLktiurQgWG3W4f1gQcC2kf//taV/stfApxMFhGRmiT15oCHDrW6uNu3A7B6tZUZ7N0b2rYNYPtPPgmPPw533GF5LEVERKogtQJwXp4leO7efd9iqPvus0ILQ4cGsP1Jk+DOO6FTJ3j00QA2KCIiNVVqBeChQ6FuXYu6WB2EsWOt1GBGRgDbP/lk6NIFXnlFQ88iIlItqROAP/vMAmO/fvui7YABcPTRNmVbLdu2WbKNf/s3GDcOGjasfntFRKRGS50APGMGpKdb1AVmzbK7Bg+uZh2ErVtt6fSgQcG0U0REhFQKwHfdBbm50LQp3sOwYdYR/sMfqrHNPXuga1fb7qWXBtVSERGRFAjA3luAhH21BXNyYO5cGDgQGjSoxrbvuMNqF44ZA5dcUu2mioiIFIt+AJ4+HX7yE0uIETN8OBx3HNx4YzW2+/bbdv3SPfdYdSMREZEARTsRR1GRzc2edNK+HuqcOfbzxBPV7P3++CP88pcwYkQwbRURESkh2gF43DhYtAheeskK32N1flu0gJtuqua2O3e2631VVlBEROIgukPQc+daGcCTT4arrwYs/fOsWXDvvZCWVsXtTpliKSaLihR8RUQkbqIZgOfPh44dLe/zunXw0UeAzf02bw59+1Zxu5s3Q58+tuiqoCC49oqIiJQSzSHonJz9AbKwEHJyeN9nMWMGjB5tFQKrpF8/C+rTpu0b0hYREYmHaPaAs7MtQNaubb+zsxk+HI45Bm6+uYrbHD8eXn3V0lieeWaQrRURETlINHvAWVlW6zcnB7Kz+cBlMX06jBpVxSyRu3bBLbdA+/bKeCUiIgkRzQAMFoSzsgAY/lto1sxiaJWkpcFrr9lG6tYNro0iIiLliG4AjvnoI8uZ8eCD0KhRFTawfbuVLszODrppIiIi5YrmHHAJw4dbBspbb63Ci/Pz4dRT4emnA2+XiIhIRSIdgBcsgKlT4e674cgjq7CBW26xUoPq/YqISIJFMgC/tWgdFzw0i4t75FMnbS+tLlx/+BuZPt1WPg8bBqedFngbRUREKhK5APzWonUMGv8ZeSvqsmt1cxpm5jFi+qe8tWhd5TdSUGDd5tatrYyhiIhIgkUuAD/yr5Xs2lvIzuUZ1Kq/l8Znf8muvYU88q+Vld/IJ5/AqlXw8MNQv378GisiIlKOyK2CXr91FwDp2Sto1H4NteoXHHB/pZxzDqxeDRkZ8WiiiIjIIUWuB5yRblUWnIO66bsOuv+QvvgCvIeWLVVsQUREQhO5ANz/0rak1a19wH1pdWvT/9K2h35xXh6ccQY89licWiciIlI5kRuC7ty+JWBzweu37iIjPY3+l7bdd3+FBgyAOnWgR484t1JERKRikQvAYEG4UgG3pLlzYdw4y9yhuV8REQlZ5Iagq6SoyC43atkS7rkn7NaIiIhEswd82PLy4Ouvq1ksWEREJDg1IwC3aQO5uVWsVSgiIhK81B+CXrjQMl8deSTUSv3DFRGRaEjtiLRuHVx8Mdx7b9gtEREROUBqB+A//tF6v/36hd0SERGRA6RuAF60CF56yVY/n3xy2K0RERE5QOoG4CFD4KijYNCgsFsiIiJykGoHYOdcbefcIufc5CAaFIgdO2DNGrvmt0mTsFsjIiJykCAuQ7oDWA40DmBbwTjySFi82OZ/RUREklC1esDOueOBy4Hng2lOAPLyYOtWu+SoXr2wWyMiIlKm6g5BPwHcCxQF0JZg9OkD559vJQdFRESSVJUDsHPuCmCj9/7jQzzvJufcQufcwk2bNlV1d5Uze7b99O2rWr8iIpLUnK9iT9E59yBwLVAANMDmgMd773uW95rMzEy/cOHCKu3vkLyHiy6yIejVq6FBg/jsR0REpJKccx977zPLeqzKPWDv/SDv/fHe+1ZAd2BWRcE37t55B+bNg8GDFXxFRCTppc51wDNmwIknQu/eYbdERETkkAIJwN77HO/9FUFsq8oeftiyX9WvH2ozREREKiP6PWDv4auv7HbTpuG2RUREpJKiH4DffBNOOQXmzw+7JSIiIpUW7QBcVARDh0Lr1vCLX4TdGhERkUoLIhVleN54A5YuhVdfhTrRPhQREalZotsDLiyEYcOgXTvo1i3s1oiIiByW6HYbFy6E3FwYOxZq1w67NSIiIoclugH43HMtAJ9wQtgtEREROWzRDcAAJ50UdgtERESqJLpzwCIiIhGmACwiIhICBWAREZEQKACLiIiEQAFYREQkBArAIiIiIVAAFhERCYECsIiISAgUgEVEREKgACwiIhICBWAREZEQKACLiIiEQAFYREQkBM57n7idObcJWBPgJpsB3wa4vTDpWJJPqhwH6FiSVaocS6ocBwR/LCd5748p64GEBuCgOecWeu8zw25HEHQsySdVjgN0LMkqVY4lVY4DEnssGoIWEREJgQKwiIhICKIegJ8NuwEB0rEkn1Q5DtCxJKtUOZZUOQ5I4LFEeg5YREQkqqLeAxYREYmkSARg59xvnHMrnXO5zrlCAS3kAAAEuElEQVSBZTzunHN/jT3+qXPurDDaeSjOuROcc7Odc8udc8ucc3eU8Zxs59w259zi2M+QMNpaGc65L51zn8XaubCMx5P+vDjn2pb4v17snNvunLuz1HOS9pw4515wzm10zi0tcV9T59w7zrlVsd9HlfPaCt9XiVbOsTzinFsR+/t50zmXXs5rK/xbTLRyjmWYc25dib+jy8p5bdKcl3KO4/USx/Clc25xOa9NtnNS5udvqO8X731S/wC1gdVAa6AesAQ4vdRzLgOmAQ44D/gw7HaXcyzHAWfFbh8JfFHGsWQDk8NuayWP50ugWQWPR+K8lGhvbeAb7Lq9SJwT4CLgLGBpifseBgbGbg8ERpVzrBW+r5LkWH4N1IndHlXWscQeq/BvMUmOZRhwzyFel1TnpazjKPX4o8CQiJyTMj9/w3y/RKEHfA6Q673P897vAV4DOpV6TifgJW8+ANKdc8cluqGH4r3f4L3/JHZ7B7AcaBluq+IqEuelhH8HVnvvg0wWE1fe+3eBzaXu7gT8I3b7H0DnMl5amfdVQpV1LN776d77gtg/PwCOT3jDqqCc81IZSXVeKjoO55wDugFjE9qoKqrg8ze090sUAnBL4OsS/17LwUGrMs9JKs65VkB74MMyHs5yzi1xzk1zzrVLaMMOjwemO+c+ds7dVMbjUTsv3Sn/wyQq5wSgufd+A9iHDnBsGc+J2rkB+H/YiEpZDvW3mCxuiw2nv1DOUGeUzsuFQL73flU5jyftOSn1+Rva+yUKAdiVcV/ppduVeU7ScM41Av4J3Om9317q4U+wIdCfA08CbyW6fYfhAu/9WcBvgVudcxeVejwy58U5Vw/4D+CNMh6O0jmprMicGwDn3GCgAHilnKcc6m8xGTwNnAKcCWzAhm9Li9J56UHFvd+kPCeH+Pwt92Vl3Fft8xKFALwWOKHEv48H1lfhOUnBOVcXO/mveO/Hl37ce7/de/997PZUoK5zrlmCm1kp3vv1sd8bgTexYZqSInNesA+JT7z3+aUfiNI5ickvHuqP/d5YxnMic26cc72AK4BrfGxCrrRK/C2Gznuf770v9N4XAc9RdhsjcV6cc3WArsDr5T0nGc9JOZ+/ob1fohCAFwCnOudOjvVSugMTSz1nInBdbNXtecC24iGFZBKbMxkDLPfeP1bOc1rEnodz7hzsHH2XuFZWjnOuoXPuyOLb2GKZpaWeFonzElPut/monJMSJgK9Yrd7ARPKeE5l3lehc879BhgA/If3fmc5z6nM32LoSq1/6ELZbYzEeQE6ACu892vLejAZz0kFn7/hvV/CXplWmR9sNe0X2Cq0wbH7bgZujt12wFOxxz8DMsNucznH8Uts2OJTYHHs57JSx3IbsAxbZfcBcH7Y7S7nWFrH2rgk1t4on5cjsIDapMR9kTgn2JeGDcBe7Ft6b+BoYCawKva7aey5GcDUEq896H2VhMeSi829Fb9fnil9LOX9LSbhsbwcex98in14H5fs56Ws44jd//fi90eJ5yb7OSnv8ze094syYYmIiIQgCkPQIiIiKUcBWEREJAQKwCIiIiFQABYREQmBArCIiEgIFIBFRERCoAAsIiISAgVgERGREPx/r6hH8AGBemwAAAAASUVORK5CYII=\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:                Fri, 24 Apr 2020   Prob (F-statistic):           5.69e-38\n",
      "Time:                        14:17:04   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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeViU5dfA8e/DAAJuuC/gmrmjiFtm7lu5ZZqVluavzMqyxbT0rSyz0tK00rLMsk1LSyPTcpfcd3AXNVdQEVFcWYaZ+/3jZhEFGWA29HyuiwvmmZnnuR1wztzbOYZSCiGEEEI4l4erGyCEEELciSQACyGEEC4gAVgIIYRwAQnAQgghhAtIABZCCCFcQAKwEEII4QKezrxY6dKlVdWqVZ15SSGEEMJltm/ffk4pVSar+5wagKtWrcq2bduceUkhhBDCZQzDOJ7dfTIELYQQQriABGAhhBDCBSQACyGEEC7g1DngrJjNZqKiokhMTHR1Uwo0Hx8fAgMD8fLycnVThBBC2MDlATgqKoqiRYtStWpVDMNwdXMKJKUUcXFxREVFUa1aNVc3RwghhA1cPgSdmJhIqVKlJPjmg2EYlCpVSkYRhBCiAHF5AAYk+NqBvIZCCFGwuEUAdjWTyURwcDD16tWjYcOGTJ48GavVesvnHDt2jDlz5jiphUIIIW43Lp8Dzq3Q8GgmLo3kVHwCFf19GdmlFr0aBeTrnL6+vkRERABw9uxZ+vfvz8WLFxk7dmy2z0kLwP3798/XtYUQQtyZClQPODQ8mtELdhMdn4ACouMTGL1gN6Hh0Xa7RtmyZZkxYwbTpk1DKcWxY8do1aoVISEhhISEsGHDBgBGjRrF2rVrCQ4OZsqUKdk+TgghhMhKgeoBT1waSYLZkulYgtnCxKWR+e4FX6969epYrVbOnj1L2bJlWb58OT4+Phw6dIh+/fqxbds2JkyYwKRJk1i0aBEA165dy/JxQgghRFYKVAA+FZ+Qq+P5oZQC9D7lF198kYiICEwmEwcPHszy8bY+Tggh7OryZTh3DmQLYoFToAJwRX9forMIthX9fe16nSNHjmAymShbtixjx46lXLly7Ny5E6vVio+PT5bPmTJlik2PE0IIu1m8GLp3B8OA5GTwLFBv6Xe8AjUHPLJLLXy9TJmO+XqZGNmllt2uERsby3PPPceLL76IYRhcvHiRChUq4OHhwU8//YTFoofAixYtyuXLl9Ofl93jhBDC7i5cgCef1MEXYMwYSB21EwVHgfq4lDbPa+9V0AkJCQQHB2M2m/H09GTAgAEMHz4cgKFDh9KnTx9+++032rVrR+HChQFo0KABnp6eNGzYkEGDBmX7OCGEsCurFe67DyIj4a239FehQq5ulcgDQznxU1OTJk3UjQuT9u/fT506dZzWhtuZvJZC3MYuXIDixcHDA/76CwICICQEEhMhLAxq1ZJ5YDdkGMZ2pVSTrO4rUEPQQghxR5o/H2rXhq++0rd79NDBF+DaNXjgAQgNdV37RJ5IABZCCHd19iw88gg8/DAEBuqh5xuVKAF+fhAV5fz2iXyRACyEEO7or7+gbl348094/33YtAkaNLj5cYahg7ME4AKnQC3CEkKIO4afH9SoAd9+C/Xq3fqxlSrByZPOaZewmxx7wIZh+BiGscUwjJ2GYew1DGNs6vGShmEsNwzjUOr3Eo5vrhBC3KaUgh9/hA8/1Lc7dICNG3MOviA94ALKliHoJKC9UqohEAzcbxjGPcAoYKVS6m5gZeptIYQQuRUVpff0PvkkLFsGKSn6uK1lRkeNgoULHdc+4RA5DkErvU/pSupNr9QvBTwItE09/gMQBrxh9xY6WFxcHB06dADgzJkzmEwmypQpA8CWLVvw9vZ2ZfOEELczpeC772D4cB10P/sMXngBTKacn3u92rUd0z7hUDbNARuGYQK2AzWAL5RSmw3DKKeUOg2glDptGEZZB7bTYUqVKpVeivDdd9+lSJEijBgxIv3+lJQUPCW9mxDCEY4eheefh5YtYeZMuOuuvJ0nNlYv1urcGSpXtm8bhcPYFFmUUhYg2DAMf+APwzDq23oBwzCGAEMAKheQP4xBgwZRsmRJwsPDCQkJoWjRopkCc/369Vm0aBFVq1bl559/5vPPPyc5OZnmzZvz5ZdfYsrtp1chxJ3DaoXly6FLF6heXa9uDg7WCTby6vRpeOYZ+O03CcAFSK66dkqpeMMwwoD7gRjDMCqk9n4rAGezec4MYAboTFi3Ov8rr0BqZ9RugoPh009z/7yDBw+yYsUKTCYT7777bpaP2b9/P3PnzmX9+vV4eXkxdOhQZs+ezcCBA/PXaCHE7enwYXj6aVizBtau1ft60xJqZCE0PNq21LuVKunvshK6QMkxABuGUQYwpwZfX6Aj8BGwEHgSmJD6/U9HNtTZ+vbtm2NPduXKlWzfvp2mTZsCOqd02bIFciReCOFIFoue333rLfD21vO+LVve8imh4dGMXrA7vQZ6dHwCoxfsBrg5CPv7SzKOAsiWHnAF4IfUeWAPYJ5SapFhGBuBeYZhPA2cAPrmtzF56ak6yvXFFDw9PbFarem3ExMTAV0z+Mknn2T8+PFOb58QogDp2RP+/luvdP7qK53HOQcTl0amB980CWYLE5dG3hyADUP2AhdAOU46KKV2KaUaKaUaKKXqK6XeSz0ep5TqoJS6O/X7ecc31zWqVq3Kjh07ANixYwdHjx4FoEOHDvz++++cPatH38+fP8/x48dd1k4hhBtJSdHzvQADB8LPP+utQjYEX4BTWdQ+v9Vx2Qtc8EgqShv06dOH8+fPExwczPTp06lZsyYAdevW5f3336dz5840aNCATp06cfr0aRe3VgjhchER0KwZTJ+ubz/6KDz+uO37eoGK/r65Os5338le4AJG9tdcJ7vFVr6+vixbtizL+x599FEeffRRB7ZKCFFgJCXBBx/A+PFQqpTulebRyC61Ms0BA/h6mRjZpVbWT5DVzwWO9ICFEMIetm3TK5rHjYN+/WDvXnjwwTyfrlejAMb3DiLA3xcDCPD3ZXzvoKxXQQMcOADvvAPnzuX5msK5pAcshBD2cOmS/lq0CLp1s8spezUKyD7g3ujIEXjvPV0buHRpu1xfOJb0gIUQIq/Wrs3YvtG+vd7na6fgm2tpw92yErrAkAAshBC5dfkyvPgitG4NX34JCakrkwsVcl2b0pJxyEroAkMCsBBC5MayZVC/vg68L78M4eHgm83KZGeSZBwFjswBCyGErU6fhh49dA7ndevg3ntd3SJdMzgsDNq21b3gU6dc3SJhI+kBAyaTieDgYOrXr0/fvn25du1ans81aNAgfv/9dwAGDx7Mvn37sn1sWFgYGzZsyPU1qlatyjlZ6SiE82zZor9XqABLluherxsE32tL12K5rw3W/3sT1aEDTJ0Kc+a4ulnCRhKA0ft8IyIi2LNnD97e3nz11VeZ7rdYLNk889ZmzpxJ3bp1s70/rwFYCOEkZ8/qJBrNm+uhZ4B27cDHx7XtAjZM28Glbo9isprxQEFyst4KlYtkH+I6CQnw9de6RrOTFMwAvHGj3ui+caPdT92qVSsOHz5MWFgY7dq1o3///gQFBWGxWBg5ciRNmzalQYMGfP3114DOB/3iiy9St25dunXrlp6WEqBt27Zs27YNgCVLlhASEkLDhg3p0KEDx44d46uvvmLKlCkEBwezdu1aYmNj6dOnD02bNqVp06asX78egLi4ODp37kyjRo149tlnUU78AxHijqQUzJ4NdetCaCi8/74OvG7gXIyFxfVep9mwZhRSSaR4eGPGhPL0hiJFdLWllBRXN7PgiYzUC+vWrXPeNZVSTvtq3LixutG+ffsyH2jT5uavL77Q9129qlRwsFIeHkqB/h4crNSsWfr+2Nibn2uDwoULK6WUMpvNqmfPnurLL79Uq1evVn5+furIkSNKKaW+/vprNW7cOKWUUomJiapx48bqyJEjav78+apjx44qJSVFRUdHq+LFi6vffvst9Z/SRm3dulWdPXtWBQYGpp8rLi5OKaXUO++8oyZOnJjejn79+qm1a9cqpZQ6fvy4ql27tlJKqWHDhqmxY8cqpZRatGiRAlRsbGzOr6UQIm8GDNDvMc2bK7V3r6tbo5RSympVas4cpcqUtqoFxkNqa6PBKuHUebV7xgY1ig/VuokblPr6a93ukydd3dyC4eJFpX76KeN26nu0PQHbVDYxseAtwrp4MSPBudWqb+dTQkICwcHBgO4BP/3002zYsIFmzZpRrVo1AJYtW8auXbvS53cvXrzIoUOHWLNmDf369cNkMlGxYkXat29/0/k3bdpE69at089VsmTJLNuxYsWKTHPGly5d4vLly6xZs4YFCxYA0K1bN0qUKJHvf7MQ4gZp7yseHjqZRUgIDBsGOZQldbiNG7k49x/2/rqbMTEfU63Z3dRYNo+gRvrtu9Q9d9OT4Vw5UB16p+4FjorKVxrMO8KSJTBkCERH6ymGu++G1PdoZ3G/ABwWlv19fn56WKhDBz3f4e2tb7dooe8vXfrWz89G2hzwja4vSaiUYurUqXTp0iXTY/7++2+MHOZclFI5PgbAarWyceNGfLPY0mDL84UQefTffzB4MPTpo4ch+/VzdYsAsK7fiGrTlmKWZFoAszvUofHSDzGZMt66y9YoRjk2szpyP1Tqow+ePAn33OOaRru78+fh1Vfhxx+hTh1Yv14HXxcoeHPALVrAypU63+rKlRnB18G6dOnC9OnTMZvNABw8eJCrV6/SunVrfv31VywWC6dPn2b16tVZNLkF//77b3oZw/PndeXGokWLcvny5fTHde7cmWnTpqXfTvtQ0Lp1a2bPng3AP//8w4ULFxzzjxTiTmOxwOTJEBQEO3ZA0aKublG6Q2tOc6zTM5gsyRgAHiaadSh6U4fc5OvNOY9ymM5c1+uVvcBZS0nRH0zmzIG339ar2V34QcX9esC2aNHCaYE3zeDBgzl27BghISEopShTpgyhoaE89NBDrFq1iqCgIGrWrEmbNm1uem6ZMmWYMWMGvXv3xmq1UrZsWZYvX06PHj14+OGH+fPPP5k6dSqff/45L7zwAg0aNCAlJYXWrVvz1Vdf8c4779CvXz9CQkJo06YNlaXqiRD5t28fPPUUbN4M3bvr0oFuMGybnAwffQS+737KMOshrCYvDKwY3t56r28W4vwCKRx3UifjKFUKrl51bqPdXWysHiH19IQJE+Cuu6BhQ1e3CkM5cUVtkyZNVNqq4DT79++nTp06TmvD7UxeSyFyYfVqvcXos8/gscfcYvvOzvmH+XDUJeYdDmFg7yt88topShtxGYk2sul4bKnUmxIxkdydvFev4HaDf4tbUAq+/x6GD4dJk/QKcSczDGO7UqpJVvcVzB6wEELkxZYtsGkTvPSS3lZ09Chct9bDVa5eTGFVjyl0XDuGUV4NeDx0Ez0fLALU1A/IYcTv7F33cvSUDzWUrBdJd/SoXmS1YgW0aqW/3EzBmwMWQojcunYNXntNB7LJkzOGaN0g+O5+/Qcul6pCj7Wvc7BqF+7auYCeD+YuiB7tM4LHrHOIjQW++85tFpG5zA8/6HzdmzbpnN1hYVCzpqtbdRMJwEKI29vq1XqR1eTJ8MwzsHOnWwTeuDj4/p6vqD9xEOUsp7B6etNw9usUq2Nj/d/rpBVCOnkS3fObN+/OTsZRujS0aQN798Lzz+utZW7ILVrlzHno25W8hkJkISZG7+n18NCB+KuvoHhxlzZJKfjj2/PUrQuHt8ShMDAAqyWFvb/8ladz1kzYyUkCSVi4XC8ks1rhzBn7NtydJSfrnTEffKBvd+sGixeDmy9YdXkA9vHxIS4uTgJIPiiliIuLw8cN8tMK4Ra2btXfy5XTb8Q7d2a7gtiZovfGs7Tqs7QcXJsq/qfY16M0SZ5epBgemE2evH+lLKHh0bk+b9m7ixNINEkHj2es5D550s6td1Nbt0KTJjBmDBw4kJHLuQDMhbt8EVZgYCBRUVHExsa6uikFmo+PD4FusIVCCJc6e1YvsJo7V2c66tJFJ+5xMasVlg0NpeGMoXRSMYS3HU7h1jvZkRTI48U+4J4Tu9lUOYgd5WpyYmkkvRrlbhi6VFBFrBhYjkdBpWb64O2+F/jaNR10p0yB8uXhzz+hZ09XtypXXB6Avby80lM0CiFEnigFP/8Mr7wCV67Ae++5R/GEjRuJ/XUFJ39Yzf0XV/NfkYac/vkvmjzYmGOjFgOwI6AOOwIytg+eik/I9WWMQt7EmsrheSZKTwjfdZdTq/q4xKFD8PnnOoPZxx/nOLUQGh7NxKWRnIpPoKK/LyO71Mr1Bx17c3kAFkKIfBs4UAfgFi1g5kxdxSgHjn5DTlm7EdWhAyXMyZQA9rd9jtpLP8fw9gKgor8v0VkE24r+N6eitUWcXyUKXzgJJUrA4cP5abr7io/XPd0nn9SJNA4dgipVcnxaaHg0oxfsJsGsS8tGxycwesFuAJcGYZfPAQshRJ5YrRkFFLp1072htWttDr6jF+wmOj4BRcYbcl7mX7OyO/Q/Tnb8H4Y5CU8smExQp3Pl9OALMLJLLXy9MueV9PUyMbJLrbxd866HWGN1v72udvPnn/p3+/TTGR8wbAi+ABOXRqYH3zQJZgsTl0bau5W5IgFYCFHwHDgArVvDF1/o2489lqvKRY56Q752KYW/2kziroeCKJt8EsPTE0ymLNNI9moUwPjeQQT4+2IAAf6+jO8dlOce2c6uoxl99S0sFuCNN26fvcAxMfDII9CrF5Qpo/f21qiRq1NkN6yfl+F+e5IhaCFEwWE26/m+997Te3lLlcrTaRzxhrx5xk58XxpMj6Rt7KzSk2r/fIkp/sQt00j2ahRgtyHQSpUASwqnoz0IjInRowEFndmsiyWcOgXvvw+vvw5eXjk/7wb2Hu63FwnAQoiCITwcBg2CXbt0j+jzz/U2ozywxxtyaHg0f89YwF27drAm7iE6RIYx0HSCve/Oo+GYh1O3wQQ4tHDM9fPY3TcfJIkR7F67i8BKleD0aZ2Mw7MAvs1HR0PFijrYfvaZzmJVu3aeTzeyS61Mc8CQv+F+e5EhaCFEwXDxoq7l+uefeptRHoMv5H/+NTQ8ml8/m8fUGa/x2oYf+S3yMf6tG8LaxSup905fp+xBvXEe+6SfFyasbF21v+Am47BY9LaimjVh1ix9rGfPfAVfsP9wv70UwI9GQog7xqpVuk7viBF6GPfwYShUKN+nTXvjzesq6Ck/7mbirz9RyKrrgysjiZallzMpvCIPdcl382xy4zx2TLliABzeGgl9GumDJ0+6RYlFm+zZo7cUbd4MXbtCx452Pb09h/vtRQKwEML9XLgAI0fCt9/q3s+LL4KPj12Cb5q8vCFbrbB82EJmfzmUCpzCbJgwUKSYPNlUOcipi3puvNbZksWxYuB/7qzeB9y2rdvmQL7J1Km6WEbx4jB7tl5AVgAyWeWXBGAhhHuZP18H3NhYvejm3Xd18HWxgwdhVp9FjN/zIHu86zKk65v4FLmakcUqoA4BTlzUc+M8tsXTkzNGOSolxEKtWjr3tbtLq11co4ae158yRa90doG0+fTouEQCSvk4JVGHBGAhhPuIioL+/fV+z8WLISTE1S3CnKyY+dYxXv28Gn6FHqDbwBmcG9qRo4v0EHBaFitnL+pJW1jkdfkiJa9d5FjJAKaXeIZLfrV4wmmtyKPLl+H//k+vYn/3XV0w44EHXNac0PBoXv/5AIF/XGHgxc3sfqAko68mA45N1CEBWAjhWlarLpreubOerwwL08n187DdxN7+++h3vMaMol9yLJu7H2H8jFJUqPCMvtPb26WpDXsFVyTw7wU0ef8lAFq/v4xNrd9gz6bCfAY6OUmZMvD9905rk03+/huee05/2Hr1VVe3BqXg9Y/iKRPqyaKkfniThHmuJ48/9gETl3o79neqlHLaV+PGjZUQQqQ7dEiptm2VAqVWr3Z1a9JdvZSi1jYapqygrKAsJi+l1q1zdbMyHDqkVMeO+nUrWlSpjz5SymxWY982q0qcUElJSqkOHZS65x5XtzTD2bNK9e+v21ynjlIbNri6RerECaW6ddNN+sJvsLLqeKzMhof6qPVAVfWNRfm+BrBNZRMTC8gMvRDitpKSohNqBAXp/b3ffKMLqKcKDY+m5YRVVBu1mJYTVtktRaQt/l10mUNlWnBf+FQADMADK6xZ47Q23JLZrAtNbNmiM4FduKDnyj09uf/Ap5ygMqcOXNKZOdypJGF0NISGwjvv6N+5A/dH58Rq1S9d3bqwblUyVbsdZHXPu7AYHumlITdVDnJ4og4ZghZCOF+3brBsGTz0EEybppMupHJV4vwL5xUjXzf49tuizC7WmKL9ulN97gRd7D2LVJJOt20bNGqkh+Z//FEvtKpYERIS9BB+3br41NBbjs5FRFE1MND1yTiOH9dB9+WXITgYTpzIc/Yye9m/H555BjavN/NVjU8YaJ7J0tf/YuTyYB55/KP0RXX7q9ZnvIPn9KUHLIRwjmvXdKIFgCFD4PffYcGCTMEXXJM4P2zcWqLLNWLDrEhefx0eOjOd6rPGwMqVMG6c/u6qHtuFC/r1ato0IzlFu3YZr1tiok5WsXAhxevpAHxp70ndA7ZadRB2NotFZ7CqVw/eekunkgSXBt/kZP2rDA4G793biancjKcPj8arSTDd65ZhfO8gYuqFML3FI8TUC3FKog7pAQshHG/lSh1EXn4ZXnoJ+vTJ9qHOTJx/7tcVnB/6Fm0vbCbauyp/fHOeWoOue0CLFq4LvErBr7/qGsdxcToZSVYFFvz9wc8PTp6kzP29AUg8HAW9G8KAAU5uNJkTajzwAEyfftOHLGfbvFk3yX/PWsLLvEmdc+swCpfXHwAfegiAXqWdX5pQArAQwnHOn9eBY9YsvdezQYMcn+KMxPlKwea+k2g2/3VKobB4eFJu4UwCurhuXvImL7ygg1fTprB0qe66ZcUwdG83Kgq/GhWxYqBOnITmT0Pz5s5t87VruncObpFQ4+pV3QH/7DPoXmojf3h3wRSboKtmff+9XnnvQjIELYRwjH/+0atcfvxRl8fbtcumeVR718m90eHD0L49mOeHYqAwAJOh8NyxxS7nzxezWc/pQkbBiY0bsw++aQID9bYeLy8mVpxCWKHUfJhK6XM6Wni4vpafn87TvX+/3s/twuC7bBnUrw/ff3qBsLpDmdf/D0yW5IwHbN/usralkQAshHAMPz8dGLZuhQkTwNe2HqyjEuenmBWLH/mB5+qtJTwc4oe/r9tkMjlkkVWuV3Jv2gSNG8Pbb+vbbdvaXuP4uhXPaxu9zIqrqT35ChX0CmlHuXgRnn9eJ0z5+Wd9rH17KF3acdfMQVycLprVpQt0T5pPTMm6tD4wA59ihfTv2UG/77yQIWghhH1YrfDVV7qA+tixelvR1q156gXZO3H+3kVHudT/WbpdXo5fpYHU2tSKihXbwsMrb1mvN69ytZL74kUYPVq/dgEB0Lp17i/4f/+ncykDdUuf5dq6KCBEzw9HReXnn5K90FA9TH7mjE6okTqX6gqh4dF8vCSSw5v8iV9Vn3IJseypNYx6kX/oleMzU7Oqde3qkN93XkkAFkLkX9rejvXrddfDYtE9DVcm1N+4EfPSVWz+6yyNdsxEGR6ED/6Cdl8/lzH256BFVrdayZ0pAIeF6aHamBi9OG3cOChaNPcXvPvu9B/7HvyAty7O4tq1S/gFBjpmL/DLL+vh8QYNdCBu2tT+17BRaHg0I74/SOD8BJ6O3sz6kud4qso31D78N3z0EQwfnrENy5WL6rIgAVgIkXfJyXp4+YMPoEgRvbBl4EDXV7LZuBFLuw54JCVxH1YOlb2XMit+pVFQJadc3uaV3OXLQ5UqsHChTr+ZV2fP6iIWXbviUaUSxTZe5vCBS9SoVAmWL8/7ea+nlN5T7OWlVzeXL68X2DkhZWhaoYQb035arfDae5cpt8iTf1IexotkzJc8eabjW/xzfw/mvj7Q4W3LDwnAQoi8O34cPvxQbyv69FMoW9bVLSI+JonVT/xIj6RkPLFiwYPkvq3wz2Xwze5N3xbZreQOLOatX6edO/XK8Nq1YcOG/H9giYmBoUNh3jx879Z7gWPDo6hhr2QcBw/qbWStWule+v336y8b5ee1zG44P/qYJ3OnlOPE2ruY7jMUn5REDABLCg3OHGZ6NdcX8siJLMISQuTOpUswc6b++e679fDznDluEXz//XA9ZwOCaXdkDmbDkxTDg2RPT96/Wi5X6SzT3vSj4xNQZLzp23qOrFZyNz53lIU/j9DzpWfP6gQaYJ/RgkAddImKonh9/UHj0r6TekHUqFF6pCIvkpP16EaDBhARAdWr5/oU+X0tbxzOVxaDM2uqMuyR0nhE7CDCvyGdE1djdXIaSXuQHrAQwnYLF+qe1qlTeo9pUBBUq+bqVhFz6BLh94/m/iNfctxUmaEd3yaxPBm1esvV5MSN86+3YPMcbjbSHjNxaSQXzl7g7a2/8tiGBRhlyuhtOn372neY3t8fCheGkycp3U0n40g6HAXtns7Yl5tbERF6OmH3bnj4YT3nW6FCrk+T39fy+mH7pNPFifunAebYYtxf6Tf+PtWPpBKleOnht4jyKe7UNJL2IAFYCJGzM2f0IqHfftObK+fP18HXxZSCXz6Ppc2rIXRW0Wy992WebNKKa74+AOm1eiF3mbTskY0rfSV3bCzUfU4vUpswQQdLezOM9L3AhapV5IViP1HMtyU9lYL4eL0grlix3J/z8mX480+d6jKPbHktbzVEXdHfl5Nnk6m55ApN9+9kh0884b3LczWkOEbym/i88grtj11j4tJIpgfUoaK/L+OdXBoyryQACyFuzWLRW2NOnID334eRI/U+SheL+mUNS95az7dH2qIqPU6rSb1o+sg9lJiwimv5zKSV72xcZ87ocjvvvqtr8kZGQsmSNl8/T9L2Ant5saXmE5SKB+Iv6OtOmpS+TemW/v5bV32aMAEaNoRDh/JdyCGn1zKnLVud/Ruw4oPdLLzSGx8SMBJhsM/79Og6CBrpTFa9SpQoEAH3RjIHLAChH8QAACAASURBVITI2pEjGduJvvhCZ7J6802XB98Us2J9xzFU7N+GQUfeYq1XB/r98iCVH7kHsE8mrTyfw2qFGTOgTh2YOBF27NDHHR18Qa9A/+cfAFoV20n5A2EZQ9M57QWOidFpI7t1g7/+0j1fsEsVpZxey+yGqD/84z+eegrefb40r6hv8CUBD0BhMKLohQIZcG8kAVgIkVlysl7ZXLcufPmlPtapE9Ss6dp2Afv+PsaW0g/QcuU4DMATK57WZDzWhKU/xh6ZtPJ0jv37dfKRZ5/VqSN37XLu/tiAgPTh7SePv8cbJ1/INDSdJaX0auw6dXRhgnff1R8a8rIXORs5vZY3DlErBVcPlGf75Gb89sM1dtfqQ5+rv2AYBphMePj6UOfJ7It5FCQyBC2EyJBWNmbPHr3w5uGHXd0iQKdHXv7gNNovH0VlwyCy+3Bqrpyeba1ee2TSytU5lNK5m0+dgu++07kQnb0Xet8++OkneO01UipUotp/y7l4EYpfl6byJrGxelV2UFBGz90BbvVaXj9EnXK5EOeX1yfhUHkKV7zE2rVFqf+eFQaNh5YtYd06t8liZQ8SgIUQ2sSJumhCQEC+F97Y07//6vVLzx46StXA1lRe/BX7LSY+CaxKjX3bOFy3CV19KtPLFY1bt06nOixcWFf/KV/edduxjh/Xc7c9e2KqEkixdZfZt/8SxQMDMyfjMJv1trEBA3RbN26EWrXAwzUDoiO71GLU/N3Ebgvgwura3GU5wrTS/bB+OZ7gRvfonnnah5lWrVzSRkeRACzEnS5tnrd5c53b94MPcr9iNhdsTcpwZfG/nBzyHnNP9cZS/QUaLZlAg86ehEac0ot2ildnWQu9L3VtdnmWHeXCBV3kYOZMeO89XUDBhlKLDnXdXmDfmnovcFzESXj88Ywe4+bN+tPM7t16S1Hnzg7r9dqqbuEATP+Ups6OLbztcz9tLWvwSPTFy/uCfoCrs6o5kARgIe5Up0/rnL4BATBlil7pbEMhAEdkNYLMwTN86AwaTn+OOiimeqzBPDMEn3Y6iOR3X2m+KAXz5unX7dw5vSLcltXFzpAWgE+exL++rgN8eX8UPNdFL6p6+WWYOhUqVtQjHC6uhWs268XZY8fCII9fmW78DyNR6Z747J91usvbXI5jDoZhVDIMY7VhGPsNw9hrGMbLqcffNQwj2jCMiNSvro5vrhAi36xW+OYb3fNZuBDKlbP5qfbOagQZwRMg5vBl/qkxjIbTn8VAAbpWr8+msPTH22OPbp6NGQOPPZZRZvHjj3XZRXdw3YrnUu0a0N5YTbh3c0hK0lXpP/9cJ1HZt8/l0wvbtun1af/3f9C9O3xe+0sMpX/fGAbs3evS9jmLLYP+KcBrSqk6wD3AC4Zh1E29b4pSKjj162+HtVIIYR+HD+vMSEOG6LnLXbt0qsLr3KqObU4BNCfZBcnoCwnMmgUTGv5Cl/++4ED9h7Ot1ZvdXlyHpR60WHTJQIAnnoDJk3Xt3kaNHHO9vEpb8RwTg2eJohwKaMuhWH+dB9pi0Tmnp01z6PRCTq5d04MGzZtD9ZP/snLiDn7/Hbw//sChtZndVY5D0Eqp08Dp1J8vG4axHyj4G7CEuBNZrTq5wsyZ8NRTN82v5TREnN/eZ1ZJGYqeTiJwRSGeOgVt7nuaE8ObUvehRnpxUBa1W0d2qZWpjZD7fb42i4jQq8KrVoXff9eLlWq5cYrDHTt0IAMeK/wXHhHeULiLDrwutnKl/tx37shFVtd+ndYHZsDanjDiT+jYUT/AjWr1OoOh0rr9tjzYMKoCa4D6wHBgEHAJ2IbuJV/I4jlDgCEAlStXbnz8+PH8tlkIkRvr1umarZMm6dtJSVCoUJYPbTlhVZZZiwL8fVk/qn2O9+ckLcDXObaH5sd3o04U4enj87AYniz87BiDXyhk02Lc/MxD2+TaNT05+cknUKqUHr595JECtSDocKlmnE4sQaurS13ajgsXdNXCfd9t5I0iX9LVYwneV87rOr1jx7rPEP7VqzqZydChdv09G4axXSmVZa1JmwOwYRhFgH+BD5RSCwzDKAecAxQwDqiglHrqVudo0qSJ2rZtW64aL4TIo/h4Pbz89de65uzWrTot4i1UG7WYrN4RDODohG439ZBB9z5zk+ji3+8X0nxwX7wtyXgA//nVofCf8yjfsb7t/zZHiojQ5RWPHNG9348/hhIlXN0q2yxfDj/+CLNmsbvOI3j+d4Daln0u+dyglE4Z/uKLUCN2I2G0xdOarINb2giMO4mI0BPTYWF6z7Gd3CoA27TxyzAML2A+MFsptQBAKRWjlLIopazAN0AzezVYCJEPSunh0jp19GKr4cP1opYcgi/kPL+a3yxTiYlgfL2OQqnBVxkeVP+/x90n+IJeFV6+PKxerV+/ghJ8QX9o+PlnOH2alIqVCFBRxMU5vxmnTkHv3tC3ryKk9AnmPh+Gp5H6oc3DQ6e+dAdnz2aU1gwOhv/+s2vwzYktq6AN4Ftgv1Jq8nXHr69L9RCwx/7NE0Lk2qVL8Pzzep/nli16GLVwYZueaksO5F6NAlg/qj1HJ3Rj/aj2Ngff9YvjadgQ3tj0EBbDC2UyYfgUwmif89C1Qymlg1a3bnqxUpkysH59wVwIVEnv/yUqCs8qgRTjMtH7Lznt8tenwo78+z+O1ujE4vP3EPBgE724yl0WWSkFP/ygG/rCC7rQCEDlyk5thi094JbAAKD9DVuOPjYMY7dhGLuAdsCrjmyoEOIWLBYdRCwWKF5cp4/asgUaN87VaeyRR/lGl84ls6jpWIK6VybwygHeW9YCz/X/YowbpxfeuHLBzZEjcP/9OivU+fO4pLtoT9cl4/CrpYPx+Z3ZpKG0s0OHoEMHGPpsChNKT2KPRxBVY7ZgjBmj71i5Etzhd374sM5tPmiQDsDh4U4PvOmUUk77aty4sRJC2Fl4uFJNmyoFSv3+u6tbk8majzeq/Z71lAK1rXZ/deXoWVc3STOblZo0SSlfX6WKFlVq2jSlUlJc3ar8O39e/x188ok6ExmvqvGf+uIzs0Mv+dvmKNWr+Tw12uN91cVrkTpWIVi3oWdPpU6edOi1cy0hQamyZZUqVkyp6dOVslgcfklgm8omJkomLCEKqhtX6v7yi554cwPnF20g5n9v0PLcOmK8KhH5ySIaD+/m6mZlSEnRY6WdOulSi2k9x4LO318nVklKokyN4kR5FefEKcdd7pM5Z1k4bB//nH8Sb5Ixe3iyxS+ImAnTafb6s+6zanzfPt3b9fHRK50bNNBz/S4m5QiFKKj69tUrdAcN0qXwHnvM5W94SsHkZ5bi06Mjtc6tx2qYODj5Y2q5Q/C9dk0PgV69qt+IN2zQ27Nul+AL+vd/5gyMHo2HB4wpOoWSm/+x+2WuXdN1O0YMKM0Tl37FlwQ8seBlMbM5oC6vqpou/1sE4MoVvQgxKEgXywCd4tINgi9IABaiYDl7Vr+pgC4AEBamV3E6o+B7Do5vP8eiso9TdObveJOMBwpQbF+w1OZUlQ6zfDnUr69TSS5erI+VKuUeQcKBnrsykXr7f7frOVev1h3Irz6+yDf+T/JMyncoIMUwMJs82VQ5yDlpQXOyZIn+nU+ZojOA9Ojh6hbdRAKwEAWBUvDtt1C7tg4iAPfcowvAu5glRfH3E3Pwa1KH+8/NI6lKMmZPEymGB2aTJ+sC6tmcqtLuzp2DgQN14QEvL/2B5ZFHXNMWZ5kxI72O88WigRS9GGWX08bH60JK7dtD+ysLiSldj/9dmMM3TXvx+KMfMLnVAB5/7AN2BNRxXFpQW40cqXu6vr6wdi1Mn64XJ7oZmQMWwt0dOADPPgtr1uh6qM884+oWpdu/9ARxjz5P14t/E+nfjP6dB3OoWkVWRQdxz4ndbKocxI6AOhiu6hE9+6wuOPHmm7oggY+Pa9rhTCdPwh9/QEoKCaUrUSZuP1Zr/sr9/vGH3q0TE6Nj2wcX/sFrSynCpnzN5EgTCWYLG6s2BByYFjQnSuldAJ6eepuTn5+u9pBN1jd3IAFYCHc2e7bOGOTnp5NCPPWUywqnXy8pbCNr3gtjd9g5nlVh7BgwhUbfDePapH8hPoEdAXXYEZBRZ9apPaJjx3SgLV8ePvoI3n1XzwHeKSpV0htyT5/GWiGQypHLiYnR28Jz6/Rpncnq1IIN/FB8Mnf/3wNUHfc0XJ0E3t609fJivKPTgtri6FH9YatlS3jnHb2nu5sbrDvIgQRgIdxRSor+JN+sGTz6KEycmKuygY4U+e4vVHvvf7RTKbQ2eZM0aw4hA3oBTi6UcCOLRde7ffNN6NVLf3ipUcPx13U31+0F9qxWiWJhlzm4/xIVKtheBUkp+O47ncP5/iu/s954FI+LVhgfCl3rZtrH26tRgPMDbpqUFPjsMz0tYzKlD70XFK7/KC2EyHD+vB5iTpunvPtundvXDYLvpXPJLGr2HneNHYCXSsITC4VIpljU/vTHOCKRh0127dJB4dVX9fDj+PGOvZ47S8uGdfIkliHP40MCx87bHnwPH9Z5M4YMtvB+ycnMpj8eyprxgLAw+7Y3r/bs0esgRozQDd63Ty+2KkCkByyEO1BK7+N99VWdjWn4cN2jM5lyfq4TrJ+8mVJvDKZ7yh4iK3WkZuw6MJuzTCvo9B7R/Pl6C1aJEjBnjltsx3KpwEC959XDg4CahUlCTwvnJCVFLxgeM0b/WlcP/J7WP76mh3W3b8/29+0yycl6y9W8ebrnWxB/59ll6HDEl2TCEiILUVFKde6sswc1a6ZURISrW5QuJkapEQ/sURYMddozUO2f+Je+Y8MGpT78UH93lcRE/T02VqkXXlDq3DnXtcVNWS9dVl94vqS+7L38lo8LD1cqJESpQiSol9rvVtHRSqnkZKUWLlTKanWP37dSSq1cqdSYMRm3k5Jc1xYbcYtMWBKAhXC1c+eUql5dqalT3SYdonX9BrW760jVpdgG5eWl1J89Z6qk2IuubpZ24YJSQ4Yo1by527xebis5WVkw1Lw6Y7K8+9o1pUaNUqqlxwb1g88QdaV0ZWWtUEGpq1ed3NAcnD+v1FNP6ZBVo4ZSly65ukU2u1UAliFoIVxh0yb48ku90qVUKYiM1Iuu3EDMD0soPag79bDwpzGVUz+volr/p13dLO36/TDDh+txUzcZpncrI0bA8ePw229c8C5PoXM37wX+91+93KDyoRX8azyAKTEFkgyd2tTPzwWNzkJaac1hw/Se7jfe0KucfV28z9hOZBGWEM508SIMHQr33qtTCh07po+7QfC1pCj+efJXCg96GA8sGIC3h5lqx8Nc3TS9OK1PH53rumxZXelp4kS33uPpUnFxsHEjABeLBVL8UsYkcHy83rHTti2USjrFP4X7YFIp+k4PD1202V3ExcHTT+vUkdu2wYQJt03wBQnAQjiHUvDbbzqT1ddfw0sv6VWbbrJNZu9uK2vL9uaBH/sRV7gyqlAhMJkw3GXRjZ+fLh04fjxs3ZrrMot3nMBAvYk3JYXYIuUpkxRF1ZF/U3vALqrfbeH7b8yMGAEr91XAq9v9+oOMu9TqtVr1SIdSULq0TkCzeTMEB7u2XQ4gAVgIZ0hJ0ctLK1TQbyaffgpFi7q6VSQlKt55Bxo19mBrYgO2PzGZyvG78Vi92vW1Ww8dgiee0LmvfXx0D2jUKJ1SUtxaajKOpcvD2UVxvEkm5vemRP4cxENJM4ktWY2JQ4/iV9iAuXP1aIyrf9+gp2LattUjHWk5u4OD3WKEyCGymxx2xJcswhJ3lORkvbAqbcHI8eO6Dq2b2DFnv9rme59qx0r1xBN6MbFbMJuVmjBBKR8fpYoXV2r9ele3qOBZvFgpUM8897kq23ejuocN6iPjNbWlaJBSoCKq1Ffq0CFXtzJDcrJS77+vlLe3UiVKKDVrll59fRtAFmEJ4WSbN+uJtp07de9t8GCoXNnVrQLg6pI1nHhmHPWiwrjmUZRJb10kZJyrW5Vqxw79WoWHw0MPwbRpULGiq1tV8NSoAZ06EXMthXu91jPXGI2XMsNlmNG0FxPaPcURN5n+APTvevFinYDm88/dIvGMM8gQtBD2dPGiTp7booVetblggV5E4iZ2vfItvg+0pU7UCjwNhfevPxEy7iFXNyvD66/rucv58/VrJ8E3b2rWhGXLOFenIW2PbcNLmTEAi+FBvG8xKpQo7OoW6rrMSUn655df1rWZ5869Y4IvSAAWwr66dYMvvtDbJvbt05/s3SBDT2wsPP44HPvsTwwUoBe8+h3e5eKWoecfT53SP8+apV+33r1d26bbxMgutdhUsxlJnt7p5SHDqwe7plrR9ZYv1wUy0lKGduoEDz7o2ja5gAxBC2FPixdDdDTUrevqlgB6IemqUcv48csr/JbUm/b/ewN+WQHmZNeveI2P1/tVv/1Wb8364ouMPMYi/7p0oVfp0jDiY17y8aTGvm0crtuER4f0dl3xhPPn9f7tH37QvfT27V3TDjdh6Dli52jSpInatm2b064nxJ3sZEQcB7q9RqdTP7C7SAs8Nq6nXn1D7w8NC9PB11UrXhcs0EP1Z8/Ca6/pkoG30f5Ot9CpE1y+rJO+uINly2DAAL2394034O2374j6zIZhbFdKNcnqPukBC2EvixfrN7sxY1y6VcaydgOHX5lKmR1LactltnZ+k5AFb2EqnDoU3qKFa7eaTJ2q90EHB8OiRRAS4rq23M4CA/VQr7soWxaqV9eBuGFDV7fGLUgAFsJe5s/XAeW991zWhCOzNxI4oD21VBJWDOI++Z6mwwem3x/qquLpSunhx1KloF8/XVln2DDZ0+tIlSqlJ+O41T5ah/1NWK3wzTd6Tv+zz/QHrg0b3GJNhLuQRVhC2Et4uO7NueANJinBytdDtvPdwDA8UtMKGiYPyiRFpz8mNDya0Qt2Ex2fgAKi4xMYvWA3oeHR2ZzVTo4c0cOh99+vg0Hp0noeUIKvYwUG6iB4+nS2D3HY38TBg3p+97nndN3etNXOEnwzkQAshD0kJek3mkaNnH7p8F8j2V26Lf/7pgWl770bk493lmkkJy6NJMFsyfTcBLOFiUsjHdMwiwUmT4b69XX6yMGD9dJr4RyNGsFTT90y6Nn9b8Js1vmaGzSAiAiYORNWrJCc3dmQIWgh7GHvXkhJ4e3jXvw8arFThnevXDAT1m0iHTe+R5Lhy4FXvuaVyX1gU0CWi6xOxSdkeZ7sjudLVJTeSrR1K/TooSs/BQba/zoie02b6q9bsPvfRGwsfPih3o43bZpOvSqyJR9HhbCDTWt2cdXblzVFKzl0eDc0PJohz09lfP3niC5dn+4b32RfjZ6YDu6nwZT/ERpxipb/JlDtYgNa/puQ6foV/bNeZZzd8XwpXVqvav71V/jzTwm+rqLULasb2eVvIjERZszQ16pYEXbv1ushJPjmSAKwEHbwWmJl6r8yl+P+GW869h7eDQ2PZu6EeXz69euM3DuTatYjfBgylBPzplCkRvkc5/NGdqmFr1fm2rm+Xib7JWXYsAEeeCCjeEJYGDz6qMz7uVKZMjB6dLZ35/tvYu1avaL52Wf1zwBVquS1tXccCcBC2MGp+ASU4XFTsLHX8K5SEPrqP0ybN45CKhlPLHgYVjyLXEkP8jnN5/VqFMD43kEE+PtiAAH+vozvHZT/YfIrV3Qqwfvu00PxR4/q4xJ4Xa9MGTh5Mtu7bfmbCA2PpuWEVVQbtZiWE1bpD3SXLsELL0Dr1nred/ly/bPIFZkDFiK/LBZC545iZtAD/FW3Taa77DG8G737PHu7juD7qFmc9AigCFdQCswmTzZVDkoP8rbM5/VqFGDfeelly2DIEDhxQr8hf/ihW5RZFKkCA/V8/C3c6m8ibVQl7YNd2qhKq0VjKLVrO7z6qi5jWNgNcksXQBKAhcivQ4doeGwPRYI6Zzqc3+FdqxWWPBdK45nP0V6d44u7nuWzHp2pf+4I95zYzabKQewIqENAapCv6O9LdBZB2CFzvKC75RMm6OHmtWuhZUvHXEfkXaVKsHRpnp9+/ahKiWsXuertRwLwTtPHmPbVVGje3E4NvTPJELQQ+RUeDkCnx++32/Du/v3QqhUs/iaaS0UCOLt4GwG/vY2piB87AurwZYtH2BFQJ1OQd/gcb5r583W+a8OAOXP0dhMJvu4pMFDvAzab8/T0U/EJoBQ99/3LipnP88LGeQAsLllLgq8dSA9YiPzasQMKFaJ977a0z2dyCXPYeg4PncS8yGAO+L/Ds7Oep8bjz2J4edIr9THZZS1K++6wTFenT+v8zQsW6PzNkyZB+fL2ObdwjA4d9N5rszlPiU8aGld48fdP6PjfViIq1GRxbf1By2GjKncYKcYgRH517Kgr++Tzb/vQB3Op/lZ/TFixGCYu/rWWkt1cmLM5jVLw/fc6e1ViIowdq3++RXpDcRv44w/MAwaSkmRmUusBzGrcA6uHCV8vk30W790hblWMQYaghcivu+6Crl3z/PSr8WYW3TeBqm89jgdWAEweUHJXmJ0amE+TJumMSkFBsHMnvP66BN+CQildcSo+PvfPrV4dr/tasnbBKpZ06ofyMNlv5bwApAcshEstWwY/DFzJ7JiORFZsS83zmzDMZl2rd+VK11Utslh02biyZfX3BQvg6acllWRBEx8PJUrAJ5/oUYtbSUmBTz+F//6D6dOd0747gPSAhXAUiyXnx2Qh7uQ1Pu60nC5dYLt/ByK+3kyt6NUYq1bpbR2uDL5pK8C6dtVvyqVKwTPPSPAtiIoX11uEbrEXGNDZq1q0gJEj4dSpPC/aErkj/6OEyI/33oPKlW1+w1IbNvJfu8EkVanJyyu6M/6l00REQPCQZjrhQTZpJJ3CbIYPPtBl4yIj4ZVXwGTK+XnCfRmG3oqU3V7gpCR45x1dxev4cZg7F0JDpVKVk8hEjhD5ER4ORYrY9IYV+/NSSg7sxl3KghWDU6OnMupDnboyu4QHgHPm206cgJ499RzvI4/A1Kl6+FkUfLdKxnHunB52fuwxmDJF5/AWTiM9YCHyIzw8xxKEVit88+lVvAc8gofSAdYweRBY9FL6Y5xeKvBG5crpN98//tC9IAm+t49KlTIPQV+7Bl98oRdoBQToKYeffpLg6wISgIXIq3PndM/iFgH44LZLtGkDQ14tzJIqQ1DehbKs1evUUoFp1q2Dzp3h8mVdr3XFCujVK+fniYJlwAC9rgBg9Wq9mv3FF/XvH3QFI+ESEoCFyKvUDFhZBWBzkpXFvb6hTNMqlIpYyaxZ8MjRiXiErc5ykZVTSwVeuQIvvaST5x86pOf+xO2rXTtdm/nZZ6F9ez0vvHq1XmgnXErmgIXIq/LldSC7PgBv3Mjpab9zYcFquiWGs7dMW2YurELpe1Lvb9Eiy9XNI7vUyjQHDA5KI7l8uS6ecPw4DBumF10VKWLfawj3kpioSxLOmAEjRuhEKn5+rm6VQAKwEHkXFASffZZxe+NGLK3aUN5ipjxwsM8o6v32oU1l+RyeRhIyiicUKiTFE+4kFoteJLhlCzRt6urWiOtIABYir/bvhxo10ldAW1eHYVhSMABlMlGzcbFMwTc0PPqWAdbupQLTLFwIjRvrBTdz5ui9oT4+9r+OcE+FC2f+oCjchswBC5EXV65AvXq6R5kqukZbEvHB6nHzIqu0bUbR8QkoMrYZOXSvb2ws9OsHDz6o00mCXu0swVcItyABWIi82LlTD+leN/+7/WpturKY08/fvMjKqduMlIJffoG6dXUKyXHj4OOP7X8dIUS+yBC0EHmxY4f+fl0ALjFzIkuYDB9egmLemR7u1G1Gn3+us1g1bw7ffacDsRDC7UgAFiIvwsN1sorr9lAWPhTOcZ9a1Loh+ILeThSdRbC12zYjpeD8eZ23ecAAnbd56FBJJSmEG5MhaCHyIi0D1nWLrCrHhRNTMSTLh4/sUgtfr8zB0G7bjI4fhy5ddFKNlBQoWVJvMZLgK4Rbkx6wEHkxcWKmmrhxe05T1hrDnqCss2I5ZJuR1arLxo0apW9//LFULBKiAJEALERedOyY6ebJP3dQCijWJvu0lHbdZnTmjC6asHat7vnOmAFVqtjn3EIIp5CPy0Lk1s6dsHRpplrAW6/V4xWmULVXsHPaUKKEvv5338GSJRJ8hSiAJAALkVszZkDfvpnmf9eerMrvAa9QulpRx1137154+OGM4gnr1sH//mdTpi0hhPuRACxEboWH66L11823Fl67hLa1zzjmemazztkcEgJhYToDF0jgFaKAkwAsRG5YLHoI+rr9v0lnLjD92AMMsHxv/+tFRECzZvDWW7pU4L59+rYQosDLMQAbhlHJMIzVhmHsNwxjr2EYL6ceL2kYxnLDMA6lfi/h+OYK4WKHDumC5tcF4BMLIwDwvTf7BVh5NmoUnD4N8+fD3Ll677EQ4rZgSw84BXhNKVUHuAd4wTCMusAoYKVS6m5gZeptIW5vWdQAvrBKHwvobqcAvGULREXpn2fO1L3e3r3tc24hhNvIMQArpU4rpXak/nwZ2A8EAA8CP6Q+7Aegl6MaKYTb6NNHDwtfl97R2BlOtBFA1Wb57J0mJMAbb+gc0m+/rY8FBurEGkKI206u9gEbhlEVaARsBsoppU6DDtKGYWT57mMYxhBgCEDlypXz01YhXM/bGxo2zHSo9IkdHPVvREB+Ek9t2ABPPQWRkTB4cEb1IiHEbcvmRViGYRQB5gOvKKUu2fo8pdQMpVQTpVSTMmXK5KWNQrgHpeD112H9+kyH+niEsrrT+Lyf95df4L77IDERli+Hb77RNXuFELc1mwKwYRhe6OA7Wym1IPVwjGEYFVLvrwCcdUwThXADVitMm6ZTUO7cmX74xAkIv3I3ZdvXz/05ExP19y5dYORI2L37pgxbQojbly2roA3gW2C/UmrydXctBJ5M/flJ4E/7eWaT5wAAIABJREFUN08IN/Dff9ChA7z0EnTqBI8/nn7XyZ//5Tmm06i+2fbzXb4ML7yge71ms57j/egjKOrAJB5CCLdjSw+4JTAAaG8YRkTqV1dgAtDJMIxDQKfU20LcXiIjoUEDXf/3m290Csrrhod9F8zmA96kXkMbl1OsWAFBQbqIQqtWmdJZCiHuLDm+ayil1gHZpdzpYN/mCOEmrl0DPz+oWVPvxR00CCpVuulhxY+EE+nXiBZFcshKdeUKDB+ug3jNmjqN5L33OqbtQogCQTJhCXE9iwU++QSqVoWjR3W6x7ffzjL4YjZT6eJuzgXasP/Xyws2b9aLuCIiJPgKISQAC5HuwAE9LztihN6L6+Nzy4df3nqAQioJa3BI1g+Ij9fnunRJF0/YskXP9fr6OqDxQoiCRgKwEErpYvbBwXDwIMyeDaGhUKHCLZ8WtTISgBLts+gBL1oE9erBp5/qAgqgg7AQQqSSACyEYcCRI9C9u0772L+/TZWGVvg/THHiqdGtVsbB8+dhwADo0QNKldLDzj17OrDxQoiCKleZsIS4bZjNMGGC3oPbrBlMnarnaXMhIgIKlSlOhYDrDnbvroPumDHw5ps6c5YQQmRBArC484SH60L2O3dCUpIOwLkMvlitPPbHY1SsMhDD6J5xfPlyXb2oRg37tlkIcduRIWhx50hK0iuamzWDmBg9z/v++3k6lfngUTpd+I0GZc9kvqNwYQm+QgibSAAWd47vvtMBt39/2LsXHnwwz6c6vXgHAEVaXbcAa+FCvWfYnIusWEKIO5YEYHF7S0iAXbv0z888A6tWwQ8/5LvE36V/wzHjSZVu1+WA/uMPmDULPGVmRwiRMwnA4va1bp0uHXj//ToQe3pCu3Z2ObXn7nD2G3WpGXTd1qLwcAgJsWkFtRBCSAAWt5+rV3XhhNat9XDwjz/aPfnFuau+HCjdKqOzm5Skh7Ub2ZAVSwghkFXQ4nYTE6OzWB09CsOGwYcfQpEidr2EUtCbBfTsCY+kHdyzB1JSJAALIWwmAVjcHqxW8PCAsmWhWzd45BFdbcgBTp+G2Fg9up3poL+/HoIWQggbyBC0KPj++UenfTxyRM+/Tp3qsOALcGnMJDbTLHMN4O7ddRas6tUddl0hxO1FArAouM6fhyefhK5dwWTSJf+cYdNGSnCBoJAbkncYhizAEkLYTAKwKJj++APq1oU5c3Ryje3boUGDbB8eGh5NywmrqDZqMS0nrCI0PDrPly5xLJzIwiEUL556wGKB5s11W4QQwkYyBywKpuXLoWJFWLJEVzG6hdDwaEYv2E2C2QJAdHwCoxfsBqBXo4BbPfVmFy5Q7upRLtQdknHs4EFdalAScAghckF6wKJgUEqXCdyyRd+eNEkXPcgh+AJMXBqZHnzTJJgtTFwametmJGyMAMBofN1iq/Bw/V0WYAkhckF6wML9RUXBc8/B4sV6zrdZM/Dzs/npp+ITcjweGh7NxKWRnIpPoKK/LyO71Mqyd3z4dGH205dSHa/bbrRjh671W7u27f8mIcQdT3rAwn0pBd98o1c4r1oFn3wC336b69NU9M86CUfa8bQh6uj4BBQZQ9RZzRNvSGnGo8yjbpsyGQfDwyEoKPcVlYQQdzQJwMJ9zZ4NQ4bood1du2D4cL3aOZdGdqmFr1fm5/l6mRjZpRaQuyHqpL+WMtZnPJWjN2YcrF0bevTIdbuEEHc2GYIW7sVi0VmsatSARx/Vvcq+fXWSjTxKG0rObojZliFqAJYuZdji+1EYGB19YOVKnXXriy/y3DYhxJ1LArBwH/v3w9NP64QakZFQvLgOwnbQq1FAtiueK/r7Ep1FEM40dL14Map/fwA8UJCcDGFh0KSJLvIg+3+FELkkQ9DC9cxmnbM5OFgH3okToVgxp13+lkPU587p+sHdu3NOlSKJQlg9TODtDW3bwrhxEBAgW5CEELkmAVi41oULelXzm29Cz56wbx8MGODUHmWvRgGM7x1EgL8vBhDg78v43kG6x2yxkLhkNROLjqXSlQPM7L8aNXZcxvBzeDiUKCELsIQQuSZD0MI1lNJB1t8fGjfW2ax6987yobZuEcqPTEPUx4/Dl59zptx4hr1cjsUX/uPuBn6s+xaaNGkBtMh4Yni47gkLIUQuSQ9YON+GDdC0aUbxhJkzbxl8bd0ilG/r1kH37qjatTF/9gV96uzlr7/g7Q/92LZNT/dmEhsL0dFSglAIkScSgIXzXLkCL78M992ng9fZszk+xZ5ZrG7p55+hTRtYvBhrYjJ9k37C1DCInTth9OhsRpjTMmBJABZC5IEMQQvnWLECnnnm/9u77/ioqvSP45+TRgIiSFEgws/OSrEtIAhiUDGIKE2xsequqIioSHEpKiDSVVYFKSqKriuIUqyAgBQ1qBSpEQsgggiIUhRCkpnz++MMEmLKkExyZ5Lv+/XKK1PuzJybm5lnzrnnPA9s2QLdu7tJV+XL5/uwoJcIFYbfj+3ZE/x+DODH0K/t1zSckc/qp8RE6NVLAVhECkQBWIrHzJlu5vDSpa4HHKSglghRwPPEn38O553H6m8SmFHuSf69+x7iyCA6Po6L/52U//hQ3bouJ7WISAEoAEvRmTULqld3pfpGjXLdyYSc00Lmpk9y7WMqGcGxWaygANWODhyA/v2x48Yx/9LHaf3ZI1SqdBvNh5zN5VGLMC2S3Azn7LuTLcgPOd3H5W0vhfj449onERHQOWApCjt3QqdO0L49jBnjbitX7riDL+SzRCgg6PPEKSlwxx1w1lnYceN4rUJ3Oix5kM6d3eqnKx5pgunfL9fgm3Uy2L6de7j8pqtI7fnoce+TiAioByyhZC289hr06AEHD7rzvL17F/pp88piBUGeJ05JcZOsMjLwY7iHCcyveDcz3oSWLfNvQ/Ygf+6uTQC8dPAkNAgtIgWhACyhM22aKxd4ySWualExlefL8zyxtZCWBosW4c/0EQX4iKJdsz38Z47rmAcje5Cvu9MF4KXlTi1s80WklNIQtBSO3w/ffecuX3+96wEvXVqstXFzSyX5WP1ycM01pN38Tx6Zn0SaLUMm0USVieOaUUlBB1/466Svujs3sbtcRWJODW1CEBEpPRSApeC+/hqaN4emTWHfPleUoHPnQlUuKojs54lrnhjHtPQvuerGK8hYuITB85owemlj3rx7AQwZQvTHC3I8z5uX7EG+7q7vSa12Fn1aFd8XDREpWTQELccvI8MVTBg82I3hjhlTrMUTctLuwkTapW2FGR/CG3Ng3Tq+rJzMDXsmcGrT01j9Ivztb9nSSB7n88PRkobPt+3OjRefFvKUmCJSeigAy/H57Tdo0QJWr3Z1ep99FqpV87pVbpLVFVdg09PB52d4zGMMTxvEyHGGrl1D0yk/djLYNYV/QhEp1RSAJThZiyc0agQDB7plRuFg2TK4/37s4XSM30cG0SSeGc+Gjww1axbB661c6XJAt24N0dH5by8ikgOdA5b8LV7s0i1+/70LwpMmhUfwDeSWtpdcwu+pW0nzx5JBNCYujtsmJxVN8AVXPKJz52ItmSgiJY8CsORu3z645x5Xbu/AAfj1V69bdNScOVC3Lva55/hfhW5UP/g9o69eSPqAIcQsWoC5pGDneoOyahVccEGxTzYTkZJFQ9CSs3fegXvvhZ9/dgUHHn8cypb1ulVOejr+bt3Ztb8sHe1StldoyltTITm54JOsgubzwZo1rrCEiEghKABLzubNgypVXD7nhg2L9KWCKqTw2Wcwbhx06cJ7f7Rg5MEP+XJvLe7tUYYhQ+CEE4q0iUd9843L8qUKSCJSSArA4hxJI3nOOdC4sSueEBubSyHc0AmqkMLMmS7Jh99P5hvTGWoX83u9Jiye7eo8FCvVABaRENFJLHE1elu1cmkkJ01yt5UtW+TBF/IppOD3w7hx2Jtuwvr9AFjrZ+iVi1ixwoPgC3DjjbB+PdSp48GLi0hJogBcmvl88J//uLq2n30GY8e6Gb7FKM9CCj16QPfufB1/AWnEk0k00fFxXP54EnFxxdrMo6KjXfCN0eCRiBSOAnBp9tpr8NBDbpbz+vVw333FPrM3e47lWF8GJ6b9TvUTE3i1bFfuiptCI98y3u+5kKihQ4haePxpJEPGWvf3WrrUm9cXkRJFX+NLm8OH3USi+vXdWtbKlaFNG8/WtPZJrk2/GWs5d8s62q/7mEs3r2Rd5Tr0LPsht6+Lo3XrOqwfD7VqFcMM57xkZMBTT7kRg9q14dJLvWuLiJQICsClyWefQZcusGcPbNrk8jhfe62nTWp3YSJVvviExsP6Eu33YYHn9/fkt8px/O9/cNNNYZDvYtkytx56zRq47jq49VaPGyQiJYGGoEuDAwfg/vuhWTP44w945ZXgC+EWtZUrada3KzF+HwbwEU2D8zNJTYWbbw6D4LtqlatvvGePm409ezaUL+9xo0SkJFAPuKT7+WeXu3nbNheEhw4txkWz+fu9QiJ7o2pShUPEkElUmTg6PZ8EVTxslLWuxvHZZ7uMV88/D7fc4nnFJxEpWdQDLqkyMtzvU05xVYs+/RSeecb74GstTJ0K7dvz/rt+zk06hVq/rWH8DR/jG1iwWr0htWWLG5Y//3z44QfXBe/aVcFXREJOPeCS5khCjQED4OOP4ayz3OQhr6WkuKxan34Kn37K95Ua8s9Zv3By3ZOZPh0aN/Z4klVmpptgNXCgC7pDh0Kiav2KSNFRAC5JNm92vbV589x5S2u9bpHz6afQogU20Ct/ocz9PLB/DP0GRdOvH96t6T3i8GH391q50vV+x46FWrU8bpSIlHQKwCXFM89A//5uHe+4cYSsCn0oLFiAzcjAAJlE4zu5Ois+jKZuXY/blZ7uon+ZMi7wDhjgyix6PvNLREqDMPmElkL77ju4/HLYsAG6dfM++KanwzPP4DtwkDd/a0ka8WQQjY2N4+7/JXkffGfNcsPzKSnu+qBB0KGDgq+IFJt8e8DGmMlAG2CXtbZe4LZBwF3A7sBm/a21HxRVIyUHhw7BkCFwzTXQtCk8/bRLjxgOAeTLL+HOO2HtWp54rhKDvv8HLWvO5PLqb7HpovNpXa4W7bxq27Zt0L27W0503nlhMP4tIqVVMEPQrwBjgVez3T7GWvtkyFsk+Vu82NWj/fZbN3zatGmxFE7I1x9/wKOPYp95hgPlqnFH9Gzm77mG6m1Xs7G2j29MewCWZq92VFwmToTevV0O7JEjXVrJcPi7iUiplO84pbV2CfBrMbRF8vPbby7wJiW5IPLRR27WbhGatWo7TUcs5PS+79N0xEJmrdqe84YpKW4i05gxTKtwNzUPbKDsTdfxt/s+I+5v247pmP9Z7ai47d/vkpGsXw8PP6zgKyKeKsyJwu7GmDXGmMnGmJNC1iLJ3X//Cy+/DH36wNq1cOWVRfpyR2r1bt97CMvRWr3HBOE9e+CDD7BXXIF/zTrSKMPU2Nt44/0K/Pe/sDtzf47PnVsVpJD64w/3t5o+3V3v1Qs++ABOP73oX1tEJB8FDcDjgTOBC4AdQK4LTY0xdxtjlhtjlu/evTu3zSQ327fDkiXu8r33utSIo0a5er1FLM9avdbCtGlQpw77u/fDdyidKPzEmkym3buI1q3d9tmrHR2R2+0h8+GHrszik0+6vxm4iWnhcI5cRIQCBmBr7U5rrc9a6wdeABrlse0ka20Da22DqlWrFrSdpY/fD+PHu9qzt93mEkXExLgqRsUkt16qf+uP0LYt3HQTmzJrct/mXmSYOPxRrl5vmeSkP7ftk1ybhNjoYx6fEBtNn+TaRdPon392FRxat3ZfUpYsgWHDiua1REQKoUAB2BhTPcvV9sC60DRHALeUqHlzt5yoYUNYsMCTAvA59VKbbvmK+ZO7kTl3Po+VfZK6+5dx5sDbiP54AVFPDHFtzZJKst2FiQzvUJ/EigkYILFiAsM71C+6CVhLl7qiCYMHu56vygaKSJgyNp9sScaYN4AkXHr8ncDAwPULAAtsAe6x1u7I78UaNGhgly9fXqgGl3ipqS4PcfnybmnRbbd5Nmx65BzwuVvW0XjrGpbVOo9tZc5k2KxXuWfPk1S9+ExefBHq1fOkeUelprpSgTfe6IbGt22DmjU9bpSICBhjVlhrG+R4X34BOJQUgPOwc6crnGAtjBkDnTvDySd73SqWvDSDJnd3IsbvIy0qnqtj5/NldFOGDXPLaaOj83+OIpOWBsOHu5+TT4bvv3fLskREwkReAViZsLy2b58baj7jDJfNyhjo2TMsgi/Ll9P8kfuIDdTqjfFncFvNJaxfDw8+6HHwXbzYlQp8/HHo1MnlcVbwFZEIogDspVmz3CSriRPh7ruhWjWvW+QcOgS9e2Mvvpg/9mdymDgyiMbExfHPKUmcdprH7du82aXdTE+HOXPc8qxw+MIiInIcVIzBC36/O1/51lsuHeKsWW6yVTazVm1n9NyN/LT3EDUqJtAnuXbIJy/l+Bp/q8Sh6e/yboW7uOu3kfRouYE+DRdxQpsk72r1WgsrVkCDBm4d79tvw1VXFctyLBGRoqAAXJysdUPMUVFuyHnECDfcnENGpiMToI6swz2SBANCl8Ix6ySrTptWUGvvzwzc1ZOXfm7Boq0rqHjqCbz+KrRp43Gt3i1b3DD9nDkuz/Tf/w7tPMsmLSISEgrAxWXDBlcicNgwlw5x5Mg8N88rCUaoAvDouRs5d/M6pr7Rj1h/JgCffpfMU+nX0q2bm9t04okheamCycyEZ5+FRx91X1zGjHHnfUVESgAF4KJ2+LCLZMOGuaVFe/YE9bDckmCEMoWj78cfGTnnWeICwTeTaMpEp1Ht1s8YN+6SkL1OgVjrzvMuXeoqPj3/PNSq5W2bRERCSJOwitInn7ge2+DBbqZuaqrLIBWE4kjh+PSC8dT69WfSiSWDaNKjYtnQ7kTOqJcWstc4bocOHR2q79wZ3nwT3n1XwVdEShz1gIvSF1+4taoffgitWh3XQ/sk1z7mHDCEKIXjN99AhQpsyziF0eVfZK2Np0blTbQ67S2Wn3sOX59Wh+FFlSYyP3PmuGH6kSPdJLW77/amHSIixUA94FCy1lXemTXLXX/gAVi37riDLxRBCseMDBgxAnveeaS260udOjBvQz2a96xCZs8YXriyLTvrXlS0aSJzs2sX3HorXH01xMfDqacW7+uLiHhAmbBCZetWuO8+eO89F0g++MDrFh318sswYADs2MGiyh25ec9z1G9ZnYkTw6Ay3/Tprtd74AD07w/9+imhhoiUGHllwtIQdGH5fDB2rAtw1sJTT7meb7h4/HHswIEApBPH8PRejHilupcppv/q3HPhhRfcbxGRUkIBuLDmz4cePdww8/jxeJ8mKiAtDeLj+WnzYaphiMISY3y8df8iyt8e2jW9x5UwJCPDfUmJj3d/t+uvh44d3dpoEZFSRJ96BXHwIHz8sbt81VXu8gcfhEfw/e036NIF32Ut6NXDxw1T2nCY+D9r9ZZvkxTSlzuSzGP73kNYjiYMmbVq+183/uILl8mqXz9YvvzYxCQiIqWMPvmO19y5rv7eNdfAL7+4AJKUFB7juW+/DXXq4H/5FV7c2Jyxz2Ry3j1N8M3LuVZvKOSVMORPBw643m7jxm4d9MyZLn9zOPzNREQ8oiHoYO3a5dJGvv461K7tlhZVqeJ1qyAlxfW+lyyBJUv4odIFtPe/zx+nXMRH70Dz5gBNoGXRpJIMKmHIhg3uPPm994ZBei0RkfCgAByMfftcr3fvXhg4MHxm6qakwBVXYNPTwe/n5fh7uW/fM/TsH8ujj7rTrEWtRsUEtucQhOtGH4JXXoE77oCLL3alFsNhiF5EJEwoAOdl926oWhUqVIBBg6BFi/CZqfvdd9CtG/ZwOsbvI4NoDlaqybIPYjn//OJrxl8ShljLLRsWMmjJZEg/7M6R16ih4Csiko0CcE7S0102pmHD4KOPXPGEbt1C/jIFKjeYmQlPPYUdNIgMfzR+fwzRALFxdJ2aREwxBl84Wplp9NyNxG7+nicXjKfB96vc32zSJBd8RUTkLxSAs/vkE5cCMTXV5W8+88wieZkClRtctQruvBNWrWJp5fbctGcsHRv8wKDLFlG5Y5JntXrbXZhIu9onQa1ObpnR+PHub6jZzSIiuVIAzqpnT1fy7v/+D95/H1q3LrKXCrrcYEoKLFoESUn4+w/g4Lc76BLzFnN9HXl6MtxxRw2M8bBW78aNcM45ULYsTJ7savUmFnMqSxGRCKQAfCQVpzEu8Pbq5aoXlStXpC8b1OzhlBR33jkzE19MHL2qvc6U35O48vqTSH0OqlUr0ibm7Y8/3IS0MWNg6lS44Qa47joPGyQiEllK9xjhli1uPe8bb7jrDz4ITz5Z5MEXgig3uHevS2l5+DD4fPgPp3Pyr1/z8syTmD7d4+A7bx7Ur+8yWt11F7Rs6WFjREQiU+kMwBkZLtDWrevWzx4KXZH7YPVJrk1CbPQxt/1ZbnDmTKhTB7tiBZnEkEE0NjqO+99Kol27Ym/qsXr3huRkiI2FxYthwgSoWNHjRomIRJ7SNwS9fDl06QKrV7sh07FjoWbNYm9G1tnDx8yC/mQGPPAAWyudT3v7Lmeems6IVos4419JxHk0yQpr3U9UFFxyiVtg/MgjxbPQWESkhCp9AXjrVre+d8YMaNfO03SI7S5MdIHYWvj1V2ylyryz8iZWnXCY4Xsf5KG+sTz2GCQkeDjJ6ocfXAaryy6Df/8bOnRwPyIiUiilYwh65ky3NAagfXv45hv3OxxyEU+fDmedRcYlzenYzkfbLlV555zepCyPZfhwSMj5VHHR8/ng2WePDtMrfaSISEiV7B7wjz/C/ffD7NkuHeI997hh1BBMsipQEo2sMjOhRw/suHEAWGL5dfMyRo1qykMPQYyXRyY1Ff71L1i2LPzKLIqIlBAlMwD7fDBuHAwY4C6PGuWq8YQoMUSwSTRyDdJbt7oe+MqVABggCj9vP7CEyn2ahqSNhbJvH2za5ApP3HxzeIwUiIiUMCVzCHrdOhdwmzWD9euhTx83azdEginBl1ed3IyKVfnx13IMjH6CQyS4Wr0JcS6bVYjNWrWdpiMWcnrf92k6YmHOdXoBPvsMRoxwlxs3dku0brlFwVdEpIiUnAB84AC8+aa7fP75sGKFK9N3+ukhf6lgkmhkD9IXb13LxP/2Z8TE72l0WQK1tixmXdsBHHzH1eo1RVCrN68vAX/avx+6d3dfViZOdH9H8PDks4hI6VAyhqBnz3ZB5KefoEEDOOMMuPDCInu53ErwZU2ucSQYN9u8it5LXuWCn79lc5ma7JtYmb3V4O23TWAycRO4tmhmOeeb7vL996FrV9i+3SX9eOIJOOGEImmLiIgcK7ID8LZtbpLVrFmuXu+bb7rgW8T+UoKPLEk0AmpUTKDNnNfou/gVADKI4Z+Hp/Brw5PYOC80uSvymwiWZ0991y5XbOL0091M7MaNC98gEREJWuQOQaelQcOGMHeuO3e5cmWxVQNqd2EiwzvUJ7FiAgZIrJjA8A71jwl+fa46h86rPgTcJCuwXH3BdCZOtCELvvkNL/8l3aW1NNu8ihoV4uHkk2HBAvd3U/AVESl2kdsDjo+H559353sL0Ost7DKiP5NoZGUtTJkCycn4Np9KL/sCr3AjsaTji4mh8YOtuOx4lirlIZhqSll76on7djFs7lgu27ySlGZT3AMUeEVEPBO5ARjcUp4CKFAt3vxs2uRq4C5YwLS/PcZNXw/mgguu5ccHF1B7xyJik5K4LIQ99GAmgrW7MBF8PrYMHsldcydjjGH1w0Nocl/nkLVDREQKJrIDcAEFXYs3PykpsHAh/Pwz9qWXyLCxPBw/gUmb72LECFdeODa2CRD6ofFgJoIBtBvWA96bCVdfDRMmcH6tWiFvi4iIHL9SGYCDqsWbn5QUuOIKdy7aWtae2Iyr90/lnKREVk+Cs88OUWNzkedEsPR0t343Nhb++U/o2FFrekVEwkzkTsIqhHxr8eYnLQ1mz8amp4O1+Ihi5uGrGTQpkQULij74Qh4TwQ7/CBdd5Gr1Alx7Ldx6q4KviEiYKZU94GCWEeVqyRK46y4OZUZj/HFEk44/Oo5ub7ag6nVF2OgcHDMR7PffXYnAZ5+FxEQ477zibYyIiByXUhmAc63Fm9f53337oG9fmDCBPRVO55YD44g9qRwjWi2iXvckqnpVqxdg6VL4xz9c6cBu3WD4cFUvEhEJc6UyAEMuy4hys349JCdjd+xgcsWePLD3cW6+sxyjR8NJJ3kYeI+Ii3MVnpYudSklRUQk7JXaAByUzz6DxYs5UL8J38Q25l7/w+yp1Ih334bLL/ewXdbCW2/BmjUwZIgrtbh2bciqPYmISNFTAM6JtfDYYzB0KH4TRYw/jgfMApL6NGLQIChb1sO2/fQT3HefS7/ZoIEruRgfr+ArIhJh9Kmd3ebNkJwMTzyBtZYov49Y0nmz2yJGjfIw+FoLL74IderAnDmuxnFKigu+IiIScRSAj/D5YMwYbL16pC9dxrNxvTlEAj7javUm3prkbft++gkefNCl3lyzxtU4jtEAhohIpNIn+BHGcGjqbL5KuJxOe57njOY1adetA7U2LYKkpGIr9HAMn88NNXfo4JYWff656wFruFlEJOKV7gC8aBE88QSZ9z3AmO+uY9Tqd0mPO4HREw1dukBUVNGkkQzKunVw553wxRcwf77LulWvnjdtERGRkCu9AXjCBLdm1lrswiXMsItp2rYJ48a5zqZn0tPdOt6hQ6FCBfjf/zyeci0iIkWh9AXg/ftdQo3x47EEavVaP5NuXkS915t4n7GxTRv46COXu/k//4GqVT1ukIiIFIXSdzLx0UexEybw3gmd/pxkFZMQR/37k7wLvgcPQkaGu/zQQ/Duu/D66wq+IiIlWIntAc9atf2oQGQSAAALmElEQVTPVJN1YtLo0fAUGiU1Y/DeR/nU3sLuqhczbWgKDf9Y5N0kK4CPP4YuXdxPv36ubKCIiJR4JTIAz1q1nX4z1nLu5nX0Xz6LSzevIvXEczjDLGfvnio81KsKgwdDuXIeTrLatw8efhgmTYIzz/TuC4CIiHiiRAbg0XM30nzdEp6fPYJoa/Fh+M8vfUirdpBly8rTsGH+z5G1Bx1UsYbj8fHHrnjCjh3QuzcMHuxxei0RESluJTIAV1u3kufeGUWUtQD4ieLc05aTcn1FGjZsne/jj/Sgj5Qr3L73EP1mrAUITRAuVw6qVIEZM6BRo8I/n4iIRJySNQnr8GEAfjr573xQJpk04skgmszoGL5qVoPEysGlbRw9d+MxtYIBDmX4GD13Y8HaZS1MnQr9+7vrjRrBypUKviIipVjJ6AEfPgxDh2KnTuW5O1aycvxldPRfypWNpnJ5/Ed8Xqs+qafVY3hy7aCe7qe9h47r9ryf7Ce491545x0XcNPSVDxBRETyD8DGmMlAG2CXtbZe4LZKwDTgNGAL0Mla+1vRNTMHKSkuk1Xlym69bGoqH1bqzGMDMkm+Norruu3k5a9qMGFvJ2pUTGD4cZzDrVExge05BNsaFROCb5+1MHky9OrlviCMHg09eih/s4iIAMH1gF8BxgKvZrmtL7DAWjvCGNM3cP3foW9eLlJSXGrGtDSstRwoewo3R33I8phWTJoGN9wAxlSnS6vqBXr6Psm1jzkHDJAQG02fIHvQwNHiCX//u6tidPbZBWqLiIiUTPmOg1prlwC/Zru5LTAlcHkK0C7E7crbokUuZaO1WAxPH+xK1X+0YsMG6NSJQifUaHdhIsM71CexYgIGSKyYwPAO9fPvQfv9MHOm6/0mJsKyZW7Gs4KviIhkU9Dx0FOstTsArLU7jDEnh7BN+UtKIiMqDnzpZJo4rn46mYt7hPYl2l2YeHwznjdudMUTPv0U5s6Fq65S8QQREclVkc8EMsbcbYxZboxZvnv37tA8aZMmpD63gHlNh8D8BVzcw8MkFpmZMHKkq9O7YQNMmQItW3rXHhERiQjGBtbK5rmRMacB72WZhLURSAr0fqsDi6y1+Z4gbdCggV2+fHnhWhxurrvO5W7u2BHGjoVq1bxukYiIhAljzAprbYOc7itoD/gd4PbA5duB2QV8nsiUnn60eELXrjB9Orz1loKviIgELd8AbIx5A0gBahtjthlj7gRGAC2NMd8CLQPXS4cvv3Qzm0ePdtdbt4brr/e2TSIiEnHynYRlrb05l7uuCHFbwtuhQzBwIDz1FFSv7s75ioiIFJCyQgTjiy+gc2f49lu46y7X+61QwetWiYhIBIvIAFyklYpyYoxb2zt/vksAIiIiUkgRF4CLvFLREfPnuzW9AwdCw4aQmqo0kiIiEjIRVxEg5JWKstu3zw0zt2wJb7wBv//ublfwFRGREIq4ABzSSkXZvf8+1K3riig8/DCsWgUnnFD45xUREckm4rp1IalUlJNffoEbb4TTTnP5nBs2LNzziYiI5CHiesB9kmuTEBt9zG3HXakoq08+cROsqlSBBQtgxQoFXxERKXIRF4ALXKkou927XY/30kth1ix328UXQ5kyIW+ziIhIdhE3BA0FqFSUlbUwbRrcfz/s3w9PPAFt2oS2gSIiIvmIyABcKN26wYQJbpj55ZfdpCsREZFiVjoCsLXg90N0NFxzDZxxBjz0kJYWiYiIZyLuHPBx277dlQwcEagX0aYN9Omj4CsiIp4quQHYWnjpJahTx81urljR6xaJiIj8qWR2A7duddms5s2Dyy5zgfjMM71ulYiIyJ9KZg945074/HMYNw4WLlTwFRGRsFNyesCbN8N777nlRQ0bwo8/QvnyXrdKREQkR5HfA/b7YexYqF8fHnnE9X5BwVdERMJaZAfg776DFi1cr7dZM1i7Fk45xetWiYiI5Ctyh6APHYKmTeHwYVe96I47wBivWyUiIhKUyA3ACQkuk9X550NiAdNSioiIeCRyAzBA69Zet0BERKRAIvscsIiISIRSABYREfGAArCIiIgHFIBFREQ8oAAsIiLiAQVgERERDygAi4iIeEABWERExAMKwCIiIh5QABYREfGAArCIiIgHFIBFREQ8oAAsIiLiAWOtLb4XM2Y38EMIn7IK8EsIn89L2pfwU1L2A7Qv4aqk7EtJ2Q8I/b78n7W2ak53FGsADjVjzHJrbQOv2xEK2pfwU1L2A7Qv4aqk7EtJ2Q8o3n3RELSIiIgHFIBFREQ8EOkBeJLXDQgh7Uv4KSn7AdqXcFVS9qWk7AcU475E9DlgERGRSBXpPWAREZGIFBEB2BjTyhiz0RjznTGmbw73G2PMs4H71xhjLvKinfkxxtQ0xnxsjEk1xqw3xjyYwzZJxph9xpivAj+PedHWYBhjthhj1gbauTyH+8P+uBhjamf5W39ljNlvjOmRbZuwPSbGmMnGmF3GmHVZbqtkjPnIGPNt4PdJuTw2z/dVcctlX0YbY74O/P/MNMZUzOWxef4vFrdc9mWQMWZ7lv+j1rk8NmyOSy77MS3LPmwxxnyVy2PD7Zjk+Pnr6fvFWhvWP0A08D1wBhAHrAbqZNumNfAhYIDGwOdetzuXfakOXBS4XB74Jod9SQLe87qtQe7PFqBKHvdHxHHJ0t5o4Gfcur2IOCZAc+AiYF2W20YBfQOX+wIjc9nXPN9XYbIvVwExgcsjc9qXwH15/i+Gyb4MAnrn87iwOi457Ue2+58CHouQY5Lj56+X75dI6AE3Ar6z1m6y1qYDU4G22bZpC7xqnWVARWNM9eJuaH6stTustSsDlw8AqUCit60qUhFxXLK4AvjeWhvKZDFFylq7BPg1281tgSmBy1OAdjk8NJj3VbHKaV+stfOstZmBq8uAU4u9YQWQy3EJRlgdl7z2wxhjgE7AG8XaqALK4/PXs/dLJATgRODHLNe38degFcw2YcUYcxpwIfB5Dnc3McasNsZ8aIypW6wNOz4WmGeMWWGMuTuH+yPtuNxE7h8mkXJMAE6x1u4A96EDnJzDNpF2bAD+hRtRyUl+/4vhontgOH1yLkOdkXRcLgV2Wmu/zeX+sD0m2T5/PXu/REIANjncln3qdjDbhA1jzAnA20APa+3+bHevxA2Bng88B8wq7vYdh6bW2ouAq4H7jDHNs90fMcfFGBMHXAdMz+HuSDomwYqYYwNgjBkAZAKv57JJfv+L4WA8cCZwAbADN3ybXSQdl5vJu/cblsckn8/fXB+Ww22FPi6REIC3ATWzXD8V+KkA24QFY0ws7uC/bq2dkf1+a+1+a+3vgcsfALHGmCrF3MygWGt/CvzeBczEDdNkFTHHBfchsdJauzP7HZF0TAJ2HhnqD/zelcM2EXNsjDG3A22AW23ghFx2Qfwves5au9Na67PW+oEXyLmNEXFcjDExQAdgWm7bhOMxyeXz17P3SyQE4C+Bs40xpwd6KTcB72Tb5h3gtsCs28bAviNDCuEkcM7kJSDVWvt0LttUC2yHMaYR7hjtKb5WBscYU84YU/7IZdxkmXXZNouI4xKQ67f5SDkmWbwD3B64fDswO4dtgnlfec4Y0wr4N3CdtfZgLtsE87/ouWzzH9qTcxsj4rgAVwJfW2u35XRnOB6TPD5/vXu/eD0zLZgf3Gzab3Cz0AYEbusKdA1cNsC4wP1rgQZetzmX/WiGG7ZYA3wV+GmdbV+6A+txs+yWAZd43e5c9uWMQBtXB9obycelLC6gVshyW0QcE9yXhh1ABu5b+p1AZWAB8G3gd6XAtjWAD7I89i/vqzDcl+9w596OvF8mZN+X3P4Xw3BfXgu8D9bgPryrh/txyWk/Are/cuT9kWXbcD8muX3+evZ+USYsERERD0TCELSIiEiJowAsIiLiAQVgERERDygAi4iIeEABWERExAMKwCIiIh5QABYREfGAArCIiIgH/h9Z65Hz3lzQhwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(x, y, 'o', label=\"Data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "ax.plot(x, res2.fittedvalues, 'r--.', label=\"Predicted\")\n",
    "ax.plot(x, iv_u, 'r--')\n",
    "ax.plot(x, iv_l, 'r--')\n",
    "legend = ax.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joint hypothesis test\n",
    "\n",
    "### F test\n",
    "\n",
    "We want to test the hypothesis that both coefficients on the dummy variables are equal to zero, that is, $R \\times \\beta = 0$. An F test leads us to strongly reject the null hypothesis of identical constant in the 3 groups:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 0 0]\n",
      " [0 0 1 0]]\n",
      "<F test: F=array([[145.49268198]]), p=1.2834419617286096e-20, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "R = [[0, 1, 0, 0], [0, 0, 1, 0]]\n",
    "print(np.array(R))\n",
    "print(res2.f_test(R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also use formula-like syntax to test hypotheses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[145.49268198]]), p=1.2834419617285855e-20, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res2.f_test(\"x2 = x3 = 0\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Small group effects\n",
    "\n",
    "If we generate artificial data with smaller group effects, the T test can no longer reject the Null hypothesis: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = [1., 0.3, -0.0, 10]\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + np.random.normal(size=nsample)\n",
    "\n",
    "res3 = sm.OLS(y, X).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[1.22491119]]), p=0.30318644106317366, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res3.f_test(R))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<F test: F=array([[1.22491119]]), p=0.3031864410631729, df_denom=46, df_num=2>\n"
     ]
    }
   ],
   "source": [
    "print(res3.f_test(\"x2 = x3 = 0\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multicollinearity\n",
    "\n",
    "The Longley dataset is well known to have high multicollinearity. That is, the exogenous predictors are highly correlated. This is problematic because it can affect the stability of our coefficient estimates as we make minor changes to model specification. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/numpy/core/fromnumeric.py:2542: 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:                Fri, 24 Apr 2020   Prob (F-statistic):           4.98e-10\n",
      "Time:                        14:17:04   Log-Likelihood:                -109.62\n",
      "No. Observations:                  16   AIC:                             233.2\n",
      "Df Residuals:                       9   BIC:                             238.6\n",
      "Df Model:                           6                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const      -3.482e+06    8.9e+05     -3.911      0.004    -5.5e+06   -1.47e+06\n",
      "GNPDEFL       15.0619     84.915      0.177      0.863    -177.029     207.153\n",
      "GNP           -0.0358      0.033     -1.070      0.313      -0.112       0.040\n",
      "UNEMP         -2.0202      0.488     -4.136      0.003      -3.125      -0.915\n",
      "ARMED         -1.0332      0.214     -4.822      0.001      -1.518      -0.549\n",
      "POP           -0.0511      0.226     -0.226      0.826      -0.563       0.460\n",
      "YEAR        1829.1515    455.478      4.016      0.003     798.788    2859.515\n",
      "==============================================================================\n",
      "Omnibus:                        0.749   Durbin-Watson:                   2.559\n",
      "Prob(Omnibus):                  0.688   Jarque-Bera (JB):                0.684\n",
      "Skew:                           0.420   Prob(JB):                        0.710\n",
      "Kurtosis:                       2.434   Cond. No.                     4.86e+09\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The condition number is large, 4.86e+09. This might indicate that there are\n",
      "strong multicollinearity or other numerical problems.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/scipy/stats/stats.py:1534: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    }
   ],
   "source": [
    "ols_model = sm.OLS(y, X)\n",
    "ols_results = ols_model.fit()\n",
    "print(ols_results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Condition number\n",
    "\n",
    "One way to assess multicollinearity is to compute the condition number. Values over 20 are worrisome (see Greene 4.9). The first step is to normalize the independent variables to have unit length: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_x = X.values\n",
    "for i, name in enumerate(X):\n",
    "    if name == \"const\":\n",
    "        continue\n",
    "    norm_x[:,i] = X[name]/np.linalg.norm(X[name])\n",
    "norm_xtx = np.dot(norm_x.T,norm_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we take the square root of the ratio of the biggest to the smallest eigen values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56240.868936151135\n"
     ]
    }
   ],
   "source": [
    "eigs = np.linalg.eigvals(norm_xtx)\n",
    "condition_number = np.sqrt(eigs.max() / eigs.min())\n",
    "print(condition_number)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dropping an observation\n",
    "\n",
    "Greene also points out that dropping a single observation can have a dramatic effect on the coefficient estimates: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage change 4.55%\n",
      "Percentage change -2228.01%\n",
      "Percentage change 154304695.31%\n",
      "Percentage change 1366329.02%\n",
      "Percentage change 1112549.36%\n",
      "Percentage change 92708715.91%\n",
      "Percentage change 817944.26%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ols_results2 = sm.OLS(y.iloc[:14], X.iloc[:14]).fit()\n",
    "print(\"Percentage change %4.2f%%\\n\"*7 % tuple([i for i in (ols_results2.params - ols_results.params)/ols_results.params*100]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at formal statistics for this such as the DFBETAS -- a standardized measure of how much each coefficient changes when that observation is left out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "infl = ols_results.get_influence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general we may consider DBETAS in absolute value greater than $2/\\sqrt{N}$ to be influential observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2./len(X)**.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    dfb_const  dfb_GNPDEFL       dfb_GNP     dfb_UNEMP     dfb_ARMED  \\\n",
      "0   -0.016406  -169.822675  1.673981e+06  54490.318088  51447.824036   \n",
      "1   -0.020608  -187.251727  1.829990e+06  54495.312977  52659.808664   \n",
      "2   -0.008382   -65.417834  1.587601e+06  52002.330476  49078.352378   \n",
      "3    0.018093   288.503914  1.155359e+06  56211.331922  60350.723082   \n",
      "4    1.871260  -171.109595  4.498197e+06  82532.785818  71034.429294   \n",
      "5   -0.321373  -104.123822  1.398891e+06  52559.760056  47486.527649   \n",
      "6    0.315945  -169.413317  2.364827e+06  59754.651394  50371.817827   \n",
      "7    0.015816   -69.343793  1.641243e+06  51849.056936  48628.749338   \n",
      "8   -0.004019   -86.903523  1.649443e+06  52023.265116  49114.178265   \n",
      "9   -1.018242  -201.315802  1.371257e+06  56432.027292  53997.742487   \n",
      "10   0.030947   -78.359439  1.658753e+06  52254.848135  49341.055289   \n",
      "11   0.005987  -100.926843  1.662425e+06  51744.606934  48968.560299   \n",
      "12  -0.135883   -32.093127  1.245487e+06  50203.467593  51148.376274   \n",
      "13   0.032736   -78.513866  1.648417e+06  52509.194459  50212.844641   \n",
      "14   0.305868   -16.833121  1.829996e+06  60975.868083  58263.878679   \n",
      "15  -0.538323   102.027105  1.344844e+06  54721.897640  49660.474568   \n",
      "\n",
      "          dfb_POP      dfb_YEAR  \n",
      "0   207954.113589 -31969.158503  \n",
      "1    25343.938289 -29760.155888  \n",
      "2   107465.770565 -29593.195253  \n",
      "3   456190.215132 -36213.129569  \n",
      "4  -389122.401700 -49905.782854  \n",
      "5   144354.586054 -28985.057609  \n",
      "6  -107413.074918 -32984.462465  \n",
      "7    92843.959345 -29724.975873  \n",
      "8    83931.635336 -29563.619222  \n",
      "9    18392.575057 -29203.217108  \n",
      "10   93617.648517 -29846.022426  \n",
      "11   95414.217290 -29690.904188  \n",
      "12  258559.048569 -29296.334617  \n",
      "13  104434.061226 -30025.564763  \n",
      "14  275103.677859 -36060.612522  \n",
      "15 -110176.960671 -28053.834556  \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/stats/outliers_influence.py:693: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return self.resid / sigma / np.sqrt(1 - hii)\n",
      "/usr/lib/python3/dist-packages/scipy/stats/_distn_infrastructure.py:903: RuntimeWarning: invalid value encountered in greater\n",
      "  return (a < x) & (x < b)\n",
      "/usr/lib/python3/dist-packages/scipy/stats/_distn_infrastructure.py:903: RuntimeWarning: invalid value encountered in less\n",
      "  return (a < x) & (x < b)\n",
      "/usr/lib/python3/dist-packages/scipy/stats/_distn_infrastructure.py:1912: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= _a)\n",
      "/usr/lib/python3/dist-packages/statsmodels/stats/outliers_influence.py:733: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_internal * np.sqrt(hii / (1 - hii))\n",
      "/usr/lib/python3/dist-packages/statsmodels/stats/outliers_influence.py:761: RuntimeWarning: invalid value encountered in sqrt\n",
      "  dffits_ = self.resid_studentized_external * np.sqrt(hii / (1 - hii))\n"
     ]
    }
   ],
   "source": [
    "print(infl.summary_frame().filter(regex=\"dfb\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
