{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Robust Linear Models"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation\n",
    "\n",
    "Load data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sm.datasets.stackloss.load(as_pandas=False)\n",
    "data.exog = sm.add_constant(data.exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huber's T norm with the (default) median absolute deviation scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-41.02649835   0.82938433   0.92606597  -0.12784672]\n",
      "[9.79189854 0.11100521 0.30293016 0.12864961]\n",
      "                    Robust linear Model Regression Results                    \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                   21\n",
      "Model:                            RLM   Df Residuals:                       17\n",
      "Method:                          IRLS   Df Model:                            3\n",
      "Norm:                          HuberT                                         \n",
      "Scale Est.:                       mad                                         \n",
      "Cov Type:                          H1                                         \n",
      "Date:                Mon, 24 Feb 2020                                         \n",
      "Time:                        22:49:06                                         \n",
      "No. Iterations:                    19                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "var_0        -41.0265      9.792     -4.190      0.000     -60.218     -21.835\n",
      "var_1          0.8294      0.111      7.472      0.000       0.612       1.047\n",
      "var_2          0.9261      0.303      3.057      0.002       0.332       1.520\n",
      "var_3         -0.1278      0.129     -0.994      0.320      -0.380       0.124\n",
      "==============================================================================\n",
      "\n",
      "If the model instance has been used for another fit with different fit parameters, then the fit options might not be the correct ones anymore .\n"
     ]
    }
   ],
   "source": [
    "huber_t = sm.RLM(data.endog, data.exog, M=sm.robust.norms.HuberT())\n",
    "hub_results = huber_t.fit()\n",
    "print(hub_results.params)\n",
    "print(hub_results.bse)\n",
    "print(hub_results.summary(yname='y',\n",
    "            xname=['var_%d' % i for i in range(len(hub_results.params))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huber's T norm with 'H2' covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-41.02649835   0.82938433   0.92606597  -0.12784672]\n",
      "[9.08950419 0.11945975 0.32235497 0.11796313]\n"
     ]
    }
   ],
   "source": [
    "hub_results2 = huber_t.fit(cov=\"H2\")\n",
    "print(hub_results2.params)\n",
    "print(hub_results2.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Andrew's Wave norm with Huber's Proposal 2 scaling and 'H3' covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [-40.8817957    0.79276138   1.04857556  -0.13360865]\n"
     ]
    }
   ],
   "source": [
    "andrew_mod = sm.RLM(data.endog, data.exog, M=sm.robust.norms.AndrewWave())\n",
    "andrew_results = andrew_mod.fit(scale_est=sm.robust.scale.HuberScale(), cov=\"H3\")\n",
    "print('Parameters: ', andrew_results.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See ``help(sm.RLM.fit)`` for more options and ``module sm.robust.scale`` for scale options\n",
    "\n",
    "## Comparing OLS and RLM\n",
    "\n",
    "Artificial data with outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "x1 = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x1, (x1-5)**2))\n",
    "X = sm.add_constant(X)\n",
    "sig = 0.3   # smaller error variance makes OLS<->RLM contrast bigger\n",
    "beta = [5, 0.5, -0.0]\n",
    "y_true2 = np.dot(X, beta)\n",
    "y2 = y_true2 + sig*1. * np.random.normal(size=nsample)\n",
    "y2[[39,41,43,45,48]] -= 5   # add some outliers (10% of nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1: quadratic function with linear truth\n",
    "\n",
    "Note that the quadratic term in OLS regression will capture outlier effects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 5.03549759  0.5237213  -0.01283698]\n",
      "[0.4438335  0.0685219  0.00606313]\n",
      "[ 4.71457319  4.97859421  5.23833801  5.4938046   5.74499398  5.99190615\n",
      "  6.23454112  6.47289887  6.70697941  6.93678275  7.16230887  7.38355779\n",
      "  7.60052949  7.81322398  8.02164127  8.22578134  8.42564421  8.62122986\n",
      "  8.81253831  8.99956955  9.18232357  9.36080039  9.53499999  9.70492239\n",
      "  9.87056758 10.03193556 10.18902632 10.34183988 10.49037623 10.63463537\n",
      " 10.77461729 10.91032201 11.04174952 11.16889982 11.29177291 11.41036879\n",
      " 11.52468746 11.63472892 11.74049317 11.84198021 11.93919004 12.03212266\n",
      " 12.12077807 12.20515627 12.28525726 12.36108104 12.43262762 12.49989698\n",
      " 12.56288913 12.62160407]\n"
     ]
    }
   ],
   "source": [
    "res = sm.OLS(y2, X).fit()\n",
    "print(res.params)\n",
    "print(res.bse)\n",
    "print(res.predict())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate RLM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.97950565e+00  5.06850953e-01 -2.40390537e-03]\n",
      "[0.14730291 0.02274158 0.00201228]\n"
     ]
    }
   ],
   "source": [
    "resrlm = sm.RLM(y2, X).fit()\n",
    "print(resrlm.params)\n",
    "print(resrlm.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare OLS estimates to the robust estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fa04669ca00>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAHSCAYAAADlm6P3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3gU1dfA8e9kSUiogdBrqKHX0HtHQUBQBKWoKEVQLCCivooNlSIQmqAoRYqgdPSHAiIIoRM6oYYSekmT1N15/zgmoSSQhCS7Sc7HZ55kZ2d378rCnrlz7jmGaZoopZRSSimVFTjZewBKKaWUUkqlFw1+lVJKKaVUlqHBr1JKKaWUyjI0+FVKKaWUUlmGBr9KKaWUUirL0OBXKaWUUkplGdnS88UKFChgenp6pudLKqWUUkqpLGjv3r03TNMseP/+dA1+PT092bNnT3q+pFJKKaWUyoIMwziX0P5Hpj0YhvGDYRjXDMM4fN/+1w3D8DcM44hhGONSa6BKKaWUUkqllaTk/M4FOt69wzCMVkBXoIZpmlWBCak/NKWUUkoppVLXI4Nf0zS3ALfu2z0E+Mo0zcj/jrmWBmNTSimllFIqVaU057ci0MwwjC+ACGCEaZq7U/JE0dHRXLx4kYiIiBQOJWNwdXWlRIkSODs723soSimllFJZVkqD32xAPqAhUA9YahhGWdM0zfsPNAxjIDAQoFSpUg880cWLF8mdOzeenp4YhpHC4Tg20zS5efMmFy9epEyZMvYejlJKKaVUlpXSOr8XgeWm2AXYgAIJHWia5mzTNL1N0/QuWPCBahNERETg4eGRaQNfAMMw8PDwyPSz20oppZRSji6lwe9KoDWAYRgVARfgRkoHkZkD31hZ4T0qpZRSSjm6pJQ6Wwz4Al6GYVw0DGMA8ANQ9r/yZ0uA/gmlPGREY8aMYcKExItXrFy5kqNHj6bjiJRSSimlVGp5ZM6vaZq9E7mrTyqPJUlW7g9k/Hp/LgWFU8zdjZEdvOhWu3j6vf7KlXTu3JkqVaqk22sqpZRSSqnUkdK0B7tYuT+Q0csPERgUjgkEBoUzevkhVu4PfKzn/eKLL/Dy8qJt27b4+/sD8N1331GvXj1q1qxJjx49uHPnDtu3b2f16tWMHDmSWrVqcfr06QSPU0oppZRSjilDBb/j1/sTHm29Z194tJXx6/1T/Jx79+5lyZIl7N+/n+XLl7N7t1Rs6969O7t37+bAgQNUrlyZOXPm0LhxY7p06cL48ePx8/OjXLlyCR6nlFJKKaUcU0pLndnFpaDwZO1Piq1bt/L000+TI0cOALp06QLA4cOH+fDDDwkKCiIsLIwOHTok+PikHqeUUkoppewvQ838FnN3S9b+pEqoEsOLL77ItGnTOHToEB9//HGiZcqSepxSSimllLK/DBX8juzghZuz5Z59bs4WRnbwSvFzNm/enBUrVhAeHk5oaChr1qwBIDQ0lKJFixIdHc3ChQvjjs+dOzehoaFxtxM7TimllFJKOZ4MlfYQW9UhNas91KlTh+eee45atWpRunRpmjVrBsBnn31GgwYNKF26NNWrV48LeHv16sWrr76Kj48Pv/zyS6LHKaWUUkopx2OkZ3leb29vc8+ePffsO3bsGJUrV063MdhTVnqvSimllFL2ZBjGXtM0ve/fn6HSHpRSSimlVAZhs9l7BAnS4FcppZRSSj0e04SrV+Nvv/EGuLvDjRv2G1MiMlTOr1JKKaWUcgC3bsGuXbLt3i0/g4MhJARcXKBcOejbF6Ki7D3SB2jwq5RSSimlEvfvv7BvnwS4L70E+fPDt9/CBx+AYUClSvDEE1CvHsTESPA7fLi9R50oDX6VUkoppdS9jh+HCRMk4D1yJD5/t2ZNaNsWeveGRo2gbl3Ik8e+Y00mDX6VUkoppbIi04RTp+5NXXj9dQlso6Jg5UqZzX36aflZrx4ULiyPLVNGtgwoywe/N2/epE2bNgBcuXIFi8VCwYIFAdi1axcuLi72HJ5SSimlVOq4ckVSGMqVg6AgKFsWbt+W+9zcZBY3Nu6pXh2uX5e0hkwmywe/Hh4e+Pn5ATBmzBhy5crFiBEj7jnGNE1M08TJSYtjKKWUUiqD2LIFfH3jZ3UvXIBnnoFly6QSw4AB4OUF9etDlSqQ7a6wMBMGvbGyfPCbmFOnTtGtWzeaNm3Kzp07WblyJTVr1iQoKAiAJUuWsGHDBr7//nuuXr3KkCFDOH/+PE5OTvj4+NCwYUM7vwOllFJKZQlRUXDokAS4QUEwerTsf/tt2LtXZnqbNJEg979OtgCMH2+f8dqZQwW/b74J/03CpppatWDy5JQ99ujRo/z44498++23xMTEJHrcG2+8wbvvvkvDhg0JCAigc+fOHD58OIUjVkoppZRKhGnGz8rOng0//gj790NkpOwrWxbee0+OmT9fcnQ9POw3XgfkUMGvoylXrhz16tV75HEbNmzA398/7vbt27cJDw/Hzc0tLYenlFJKqczu2rX4tIWdO2HPHjh9GvLmlXxdZ2cYNkxmdevXh9Kl44PjKlXsO3YH5VDBb0pnaNNKzpw54353cnLCNM242xEREXG/m6api+OUUkop9XjCw6WebqVKMls7d67U1QVwcoJq1aB7d7hzR4LfUaNkU8miK7iSyMnJiXz58nHy5ElsNhsrVqyIu69t27ZMnz497rZfauduKKWUUirzCQ6WAHfIkPh6uU2bwh9/yP1Nmkhe7pYt0jntwAFJdSha1K7DzugcaubX0X399dd07NiRUqVKUaVKFSL/y6+ZPn06Q4YM4ccffyQmJoZWrVrdEwwrpZRSKou7dk3SFnbulHq5XbtK8PvSS5A7t6QsvPvuvYvSKlSA+ypQqcdn3H0pP615e3ube/bsuWffsWPHqFy5crqNwZ6y0ntVSimlsqzYRWmmCX36wPbtEBAg91ksMHIkfPml3H/8uJQb03Kqqc4wjL2maXrfv19nfpVSSimlUso04exZmdHdsUO2ggVh7VoJgG/flpneYcOgQQOoUwdy5JDHGgZk4kmxM2fg2DHo1MneI7mXBr9KKaWUUkkVGgpHjkBsPf9nnoHly+X3HDkk0L271v9vv6X/GO3s8GH46itYskTOA86fl6IUjkKDX6WUUkqpxJw/Dxs3yoyur69EdoYhC9By5oQXXoB27STgrVbt3i5pWcyOHZLNsXo1uLrZKNToIpYaJ2g50YmRHbzoVru4vYcIaPCrlFJKKSWCg6Werq8vDBwIRYrIrO5bb0lpsYYNpdRYw4bxU5ndu9t3zHZmmnJuMHYs/PUX5MsHvQaFsDfnbqKcpSxsYBCMXn4IwCECYA1+lVJKKZV1nT4t05U7dsDRo/GL1erXh44doXdv6NBBF6Xdx2aDVask6N2zR6qvTZgg5wwdp+8hKijinuPDo62MX++vwa9SSimlVLoIC5NOadu3y8xur15SicEwZHa3YUN47jn5Wb++zPSCtAcuXNi+Y3cg0dGweLHk9B47BuXKwaxZ0L8/ZM8ux1wKCk/wsYntT296CgNcvHiRrl27UqFCBcqVK8fw4cOJiopi8+bNdO7c+YHj165dS+3atalZsyZVqlRh1qxZdhi1UkoppRJkmhLsAkRFSQOJvHmhdWv48EOpzhAVJfeXKQM3bsjCtP/7P8nfjQ18VZzwcJg+XUoP9+8vqc2LFkmltoED4wNfgGLubgk+R2L701uWD35N06R79+5069aNkydPcuLECcLCwvjggw8SPD46OpqBAweyZs0aDhw4wP79+2nZsmX6DloppZRS8aKipNTYN99Ajx5QrJjM4gK4uEDNmvDBBxLg3rol1RpeflnuNwxNZ3iI4GCZ5fX0lGptxYrBmjXSbK5374TX943s4IWbs+WefW7OFkZ28EqfQT9Clk972LRpE66urrz0X+9si8XCpEmTKFOmDK1atXrg+NDQUGJiYvDw8AAge/bseHk5xh+mUkoplSXcvi0BbNOmcvvJJ2XVFUiU1qYNtG8ff/wPP6T7EDO6a9dgyhSYNk0KW3ToAKNHQ/Pmcr7wMLF5vePX+3MpKJxi7m5a7SFRb74Jfn6p+5y1asHkyYnefeTIEerWrXvPvjx58lCqVClOnTr1wPH58+enS5culC5dmjZt2tC5c2d69+6Nk541KqWUUmkjMFCC223bZDtyRKotBAeDm5vED4MHQ+PGMjWpUuz8eVm49v33EBEhE+nvvSeZI8nRrXZxhwl27+dYwa8dmKaJkcApTGL7Ab7//nsOHTrEhg0bmDBhAn/++Sdz585N45EqpZRSWUBMDBw8CFu3yoI0Dw9JLn33XcnFbdRIrrc3aRJ/zT2B9TkqeY4fh6+/hp9+ktt9+8r/8kqV7DuutOBYwe9DZmjTStWqVfn111/v2RcSEsKFCxcoV65coo+rXr061atXp2/fvpQpU0aDX6WUUiqlLl2S1IStW6USQ2io7C9fXnrj9ukjZceqVtX83FS2d69Uelu+HFxd4bXX4J13oFQpe48s7WT5T1CbNm24c+cO8+fPB8BqtfLOO+/w4osvkiO29/ZdwsLC2Lx5c9xtPz8/SpcunV7DVUoppTK2mzelBdjIkbBunewLDYWPPoLLl2XKcdEiuHBBAl+QIrLVq2vgm0pMEzZvljxeb2/YsAHefx/OnZM839QKfLed38aoP0fhe8E3dZ4wlTjWzK8dGIbBihUreO211/jss8+w2Ww8+eSTjB07Fl9fXzZu3EiJEiXijl+8eDHjxo1j0KBBuLm5kTNnTp31VUoppR4mJgaGD4e//5Z8XZAqDO7uEuBWrChBcb589h1nJmeasHYtjPwwEv+D2XHKEUmpjhf48v0cPN8s9XKlb9y5wZi/xjBjzwxMTKbumsrGfhtpVLJRqr3G48jywS9AyZIlWbNmzQP7W7ZsSXj4gwWZmzVrlh7DUkoppTIW04SAANiyRbacOcHHR3Jzd+6EkiXh+eehWTOoV0+us4OUD9DAN83ExMDSpVKy7NAhyJbXRv52h8lZ/QKGs43P/rCQI5f5WAvUTNNkV+AuZuyZwc+HfybSGhl3X5Q1is0BmzX4VUoppVQm8sknMGeOpCsA5M8PXbvG379796NrZKlUFRkJ8+bJQrYzZ6ByZSj/7FGiSgdgWMy44x6n9fCd6DssPrSYGXtmsO/yPnK75GZA7QE0KtmIgWsGEmWNwsXiQkvPlqn4zh6PBr9KKaWUShrTlJ62mzfLtmsX+PtLey8XF6nEMGqUFIO9f3GaBr5pZuX+wHtq6g5rVolLO4oxcaKkUderBxMnQpcuUO79syT0J5Hc1sP+N/yZuWcmc/3mEhwZTPVC1ZnZaSYvVH+B3NlzA1AuXzk2B2ympWdLh5n1BQ1+lVJKKZUY0wSbDSwW+OUXKQVw/brcV6IEtGwpHRAKFpQOCCrdrdwfyOjlhwiPtmINd+bIPyXo80kBbBHSzXnBAvkZe+5RzN2NwAQC3aS0Ht56fiuz9szC/4Y/ey7vwdnJmWeqPMNr9V6jSckmD5SIbVSykUMFvbE0+FVKKaWUuH9md/NmmDtXOqiVLi3lxlq2lK1MGZ3NdQDj1/sTeisbIbsrEuZXCjM6G24VrlC+7UU2zvB+4PiRHbziguVYj2o9fDn0Mh/99RFz9s/BxMTAYHDdwYxpOYbCuQqnyftKSxr8KqWUUlmVaUobLzc3ae3l7X3vzG7HjlCggNyuVw/+KwuqHMOpU3BwSXnCDhcHm0HOypfJ0/AULgXDCE3kMUltPWyaJlvObWHGnhksP7acGFtM3H1OhhOl8pbKkIEvaPCrlFJKZS3nz8Nff8GmTbJ16gTffivB7tNPQ4MGOrPr4A4ckMoNS5eC6VScXDUukqf+aZzd49MZHpbG8LDWwyGRIfx08Cdm7J7BketHyOeaj+ENhtOgeAP6r+zvkAvYkkuDX8BisVC9enViYmIoU6YMCxYswN3dnYCAADp37szhw4fvOf7FF19k6dKlXL16ldy5Jal7+PDh+Pj4cP36dQrEniUrpZRS9hYWBrlyye+tW0vgCzKj27o1tGkjt52cYNYs+4xRJcm2bdKNbd06+SMdMQIqtb3GuC3HkpXGkJDD1w4zc/dM5h+cT1hUGHWL1uWHLj/Qq1ov3JwlkC6Rp4RDLmBLLg1+ATc3N/z8/ADo378/06dP54MPPnjoY8qXL8+qVavo06cPNpuNv/76i+LFU14fTymllEoVwcHSTGLDBpnZvXULAgNlFrdbNyk/1rq1tgrOIEwT1q+XoHfLFjln+ewzGDo0tjRyUfIVsD0yjeF+vhd82Xh2IwB/nvmTLee2kN2SnV7VejG03lDqFa/3wGMcdQFbcmXI4Nf3gm+anXk0atSIgwcPPvK43r178/PPP9OnTx82b95MkyZN+P3331N1LEoppdQjRUaCs7MEshMmwHvvgdUqebzNmkm74OhoKUX2xhv2Hm2Wc38ZsqQEpiB/hMuXS9C7f79kpUyeDK+8Ir1D7vawNIaErD6+mh7LesTl8RbLXYzx7cbzUq2X8Mjhkaz3lxE5VPD75v/exO+K30OPCY4M5uDVg9hMG06GEzUK1yBv9ryJHl+rSC0md5ycpNe3Wq1s3LiRAQMGPPLYChUqsGrVKm7fvs3ixYvp06ePBr9KKaXSns0mbbo2bJBtyxapylCvHtSvLyXH2raFhg2l/q6ym7vLkAEEBoUzevkhgESD1ago+OknaUxx4oR0fp4zB/r0kfOXlDJNk80Bm5m+ezrLjy3HRJpcOBlODK03lBGNR6T8yTMYhwp+kyI4IhibaQPAZtoIjgh+aPCbFOHh4dSqVYuAgADq1q1Lu3btkvS47t27s2TJEnbu3MkszZNSSimVVqxWqbV76JAEtteuyf5KleCllyBPHrndvLlsyiGMX+9/Ty4uJN5N7d9/4fvvZfL+4kWoXVsWtHXvLn/0KRUSGcL8A/OZsXsGx24cw8PNg+erP8+vx34l2hqNi8WFVp6tUv4CGZBDBb9JmaH1veBLm/lt4lYbLuy+8LFTH2JzfoODg+ncuTPTp0/njSRcGurVqxd16tShf//+OGnelFJKqdQSFiZ5u3/8Iduzz8Knn0K5ctChgyxSa9NGroUrh5VY17S799++DdOnS0rDzZty7vL999C+/eMV2zh87TAzds9gwcEFhEWFUa9YPeZ2nUvPqj1xc3Zj6IWhmWLxWko8Mvg1DOMHoDNwzTTNavfdNwIYDxQ0TfNG2gzxXo1KNmJjv41p8geWN29efHx86Nq1K0OGDHnk8aVKleKLL76gbdu2qTYGpZRSWVzXrvD775Kn6+YGLVrI4jSAHDm01m4G8rBualeuwKRJMHMmhIZKxbnRo6FJk5S9VmIL2HpX781r3q89sIAtsyxeS4mkzPzOBaYB9/xtMwyjJNAOOJ/6w3q4tPwDq127NjVr1mTJkiU0a9YMf39/Stx1Zj1p0qR7jh80aFCajEMppVQmFxgIf/4pM7vXrkn+Lsjs7ltvydRfkybg6mrfcaoUS6ibmiUsF3nO1MVzjJzf9OwpaxRr1kz566z2X02PpfEL2IrmKsq4tuN4ufbLWWIBW3I9Mvg1TXOLYRieCdw1CXgXWJXKY0p3YWFh99xes2ZN3O/R0dEPHP/ss88m+DwBAQGpOi6llFKZ0E8/SYeCI0fkduHCksoQEwPZssE339h3fCrV3N1NLeCkheh9Xtw8VJiLFoP+/eHdd6F8+ZQ9t2mabD2/lem7p/PL0V/i1kPFLmAb2WRkar2NTCdFOb+GYXQBAk3TPGA8IiHFMIyBwECQNAGllFIqyzh7Fv73P0ljmDoVSpeWag1FisCLL8rsbvXq2knNTlJahiw5ikYVp+CO4mxfJSXK3hwOb78NKW0NEBYVFteB7dC1Q+RzzcdzVZ9jxfEVcQvYWpdpnarvIbNJdvBrGEYO4AOgfVKON01zNjAbwNvb20zu6ymllFIZSmAgjB8vAe+JE7LP0xMuXJDgt18/2ZRdpaQMWVKZJmzcKDV6N22SZhQffwyvvw4eKcxC8L/hz4zdM5h7YC4hkSHULlKbOV3m0KtaL3I450jTHgiZTUpmfssBZYDYWd8SwD7DMOqbpnklNQenlFJKObxTp+C336BsWejcWepSzZ4tC9Veew2eeAIqVNDZXQeTnDJkSWWzwapVEvTu3g1Fi0rpsoEDIXfupD9PbCDbrHQzbt65yfTd0/nzzJ84OznzbNVnGVZvGA1LNOTuq+9ZeQFbciU7+DVN8xBQKPa2YRgBgPfjVHswTZNHpU9kdKapk95KKZVpbNoEa9ZI0Bs7u/vKKxL8Fiki9au0wYRDS0oZsqSKjobFiyWV+9gxOQ+aNQv690/+x8D3gi+t57cmMiYSABOTEnlK8Hmrz3mlzisUzlU42eNT90pKqbPFQEuggGEYF4GPTdOck1oDcHV15ebNm3h4eGTaANg0TW7evImrrthVSqmM6dIlOHBAZnFBalIdOACtWsGwYfDkk1KlIZYGvg7vYWXIkio8HH74QbJczp2T9O1Fi6Qsc7YUXFvfHbibob8NJSImIm7fC9VfYG63uWRzcqjWDBlaUqo99H7E/Z6PM4ASJUpw8eJFrl+//jhP4/BcXV3vKZmmlFLKgVmtct163TrZ9u+XgPbmTVm19NNPUKyY/K4ypITKkLk5WxjZweuRjw0Jkfq833wjVeoaNYJp06RWb3Ln8SJiIlh6ZCnTd09nV+Au3LK5kc3IhomJi8WFofWGauCbyuz+f9PZ2ZkyZcrYexhKKaWyupAQCXCzZ4eJE2HUKHBygsaNJYmzUydpMgGSw6sytLvLkCW12sP16zBligS6wcFSrOP996UrW3KD3vPB5/l2z7d8t+87bty5QaUClZj2xDT61uzLkWtHdPFaGjLSMxfV29vb3LNnT7q9nlJKKfVQAQGSu7tmDWzeDEuXQrducOYM7Nwp9Xfz57f3KJWdnT8v50PffQcREdC9u2S+1K2bvOcxTZNNZzcxffd0VvlLm4QuXl0YVm8Yrcu0zrTpn/ZiGMZe0zS9799v95lfpZRSKt1duwZt28IhKW2FlxcMHy4/QVYslS1rv/Eph+DvD19/DQsWyO0+feSCQKVKyXueDWc24LPTh4NXD3Iu+BwFchRgVJNRDPYeTKm82gMhvWnwq5RSKnP7919pI7x6NRQsKNFMwYISwbz4Ijz1lKYxOKD0aECRmH37JNPl11+lu/SQITBiBCS3V9fxG8f5cNOH/HrsVwAMDP6v+f/xfrP3cc2mi+DtRYNfpZRSmdOyZTB/vgS+kZGQNy/07Sv3GYakOCiHlJYNKBJjmrBlC4wdC3/8IR+X0aPlgkChQo9+fCyrzcraE2uZtnsaG85swGJY4u5zMpxwy+amga+dOdl7AEoppVSqOHkSfHyk0wBIDu+hQzBokLTbun5dWgwrh/ewBhSpzTRh7Vpo0gRatgQ/P6nXe+4cfPFF0gPfm3duMm7bOMr5lKPbz904fuM4X7T+glW9V+GWzQ2LYcHF4kJLz5ap/h5U8ujMr1JKqYzJNGHPHli5UrajR2V/8+ZQq5YUX502TTurZUCp2YAiMTExcnHgq6/g4EHpPD19Orz0ErglvdQv+y/vZ+quqSw+vJiImAhaerZkYvuJdK3UNa5E2cZ+G7V6gwPR4FcppVTGER0tnQXy5JFr0x07Sjvh5s2lh2zXruDpKcfGliVTGU5qNKBITGQkzJsH48bB6dNQubLc7t0bnJ2T9hxbzm1hxu4ZHL1+lEPXDpHDOQf9a/ZnaL2hVC9c/YHjtfWwY9HgVymllGMLD5dAd/lyWbQ2bBh89hm0aCFRS6dO4OFh71GqVPQ4DSgSExYmLYe/+UYa9tWrBxMmQJcuUs45KS6HXubDTR/yg98PgCxgG95gOGNajsHd1T3FY1PpS4NfpZRSjuvFF+GXX6RiQ758MrPbpo3c5+oK/frZdXgqbaSkAUVibt6UVG8fH7h9G1q3lnWQrVsnLSPGNE18L/oydddUfjn6CzG2mLj7nAwnCucsrIFvBqPBr1JKKcdw86Y0m/Dzg8mTZZ+bm1Ro6N5dViMl9bq0yvC61S7+WJUdAgNllnfWLDl36tZNqjfUr5+0x0fERLDk8BKm7prKvsv7yJM9D8PqDaNxycb0X9mfKGuULmDLoDT4VUopZT9Xr0o6w6+/SnUGq1VWHn3yidSamjnT3iNUGcypU5LPO2+efJyef14aU1StmrTHnw8+z8zdM/lu33fcDL9JlYJVmPHkDPrW7Esul1wAlMhTQhewZWAa/CqllEpfgYGQMye4u8O6dfDaa1CxIrz7LvToAXXqaIUGlWwHD0pjiqVL5QLBK69IY4oyZR7+ON8LvvwV8Bd5sufhr4C/WHl8JSBth1+v/zqtPFs90HZYF7BlbBr8KqWUSnuBgTK7u2wZbNsmCZjDhkmw26ABVKmiAa9KkW3bJOhdtw5y55aA9623oEiRRz9205lNdFzYkWhbNAC5XXIzsvFIhngPobR76TQeubIXDX6VUkqlnagoaNcOtm6VurzVqsGYMfDEE3J/3ryyKZUMpikFQMaOla5sBQpIAZChQ2Vd5KOcuX2GGbtnMH339LjA1wknRjQewUctPkrj0St70+BXKaVU6rl4UWZ4L12Cr78GFxcoXx7atoVnn4VKlew9QpWBWa2SIv7ll7B/P5QoIWsjX3lFMmkexjRNNpzZwNRdU1l7Yi1OhhMtPFuw7fw2YmwxuFhcaFe2Xfq8EWVXGvwqpZR6PFevSjmyJUvgn39kn7e3RCoWC8yZY9/xqQwvKgp++knOp06ckBTxOXOgTx85v3qYsKgw5h+Yz7Rd0zh24xgFcxTkg2YfMNh7MMXzFMf3gq8uXstiNPhVSimVfLduyVRb9uzw/ffw4YeS0vDZZ9Czp0QnSj2mO3fk4zV+vFxUqF1bFrR17y7nVQmJDWbL5y/Ptgvb+NHvR0IiQ/Au5s28bvPoWbUnrtlc447XxWtZjwa/SimlkiYkBFatkhneP/6AxYvhmWdgwAApoprUWlJKPcLt2zB9OkyZAjduSPfq776DDh0evi5y2/lttJnfhkhrJAAWw8Jz1Z7j9fqv06B4gweqNqisSYNfpZ+rWnUAACAASURBVJRSDxcaKp3W1q2DyEgoVQrefhtq1JD7ixRJ2tJ6pR7hyhWYNEnKO4eGwpNPSmOKpk0f/riQyBDm+c3j0y2fxgW+BgbvNnmXsW3GpsPIVUaiwa9SSql7RUfLzO6VKzKrmysXXL8OgwZBr17QsKGWJVOp6uxZSW344Qf5+D37LLz3HtSq9fDHnbx5kmm7pvGj34+ERoVSpUAVQiJDsNqsuFhceKriU+nzBlSGosGvUkopsNnA1xcWLpSkyps3oVw5ePllCXS3bLH3CFUmdOQIfPWVZNA4OUH//tLrpEKFxB9jM238cfoPfHb68Pup33F2co5LbahfvL4uYFOPpMGvUkplZaYpwe2HH0r9KDc36NpVesI+KsFSqRTauVM+bqtWQY4c8MYbkklTokTij4lNbZi2exonbp6gSK4ijGkxhkHegyiSKz7tRhewqUfR4FcppbKa8+dlqm3hQkmubNIEeveGypVl4Vru3PYeocqETBM2bZLGFJs2STOKjz+G118HD4+EH+N7wZdfjv7CheAL/O/0/wiNCqVB8QYs7L6QZ6o8g4vlEXXOlEqABr9KKZUVRETAvHkS8G7dKvsaNYKYGPm9enXZlEplNhusXi1B7+7dULQoTJgAAwcmfp5lM21M2TGFEX+OwGbaAOhQrgOftvqU+sXrp+PoVWakwa9SSmVWkZEQEABeXpK+8P77ULgwfP65zPSWLWvvEapMLDpaquJ99RUcPSoft2+/lbxeV9eEHxMaGcr8A/OZumsq/jf94/ZbDAstSrfQwFelCg1+lVIqMzFN2LEDFiyAn3+Wa8snT0ozikOHZNpN83hVGgoPhx9/hHHj4Nw56X2yaJFUcMiWSNRx+tZppu2axg9+PxASGUK9YvX4uPnHjNs+jihrFC4WF1p6tkzX96EyLw1+lVIqs1i2TGZ3T52SqbWnn4a+feMXtRUrZu8RqkwsJERSyCdNko7XjRrB1KnQqZNUcrifaZpsPLsRn50+rD2xFouThZ5Ve/JG/TdoUKIBAB3Kd9DKDSrVafCrlFIZ1e3bEvB26AClS4OzM5QsKQFwjx6QJ4+9R6iygOvXpRPbtGkQHAzt28tHsHnzBy8y+F7w5Y/Tf/Bv9L+sO7mOo9ePUjBHQT5s/iGDvQdTLPe9J2hauUGlBQ1+lVIqI4mJkQYUc+fKKqLISJg8GYYPl0oN3brZe4Qqi7hwASZOhNmzZT1l9+7Sja1u3YSPX35sOc/98hwxNllkWdGjIvO6zeO5qs+RPVv2dBy5yuo0+FVKqYwiKkqq/58/DwUKSMe1vn0TjzaUSgP+/vD11/DTT5JR06cPjBoFlSo9eKxpmvx97m98dvqw8vhKTEwAnAwn+tfsT7+a/dJ59Epp8KuUUo7r1i1ZLn/4MMyYAS4u8NprUr3hySfltlLpZN8+aUzx66+yfnLwYHjnHcm4uV9ETASLDi3CZ6cPB64ewMPNg741+rL06FKirdG4WFxo5dkq/d+EUmjwq5RSjuXutIZVq2S2t1YtuHNHWmGNGmXvEarHtHJ/IOPX+3MpKJxi7m6M7OBFt9rF7T2sBJmmlIUeOxbWr5c08tGjJcumUKEHjw8MCWTG7hnM3jebG3duUL1Qdb5/6nuer/48bs5uDPYerAvYlN1p8KuUUo4gtiLDnDkypVaggMzy9u8vwa/KFFbuD2T08kOER1sBCAwKZ/TyQwAOFQCbJqxbJzO927dLoPvllzBkCOTNK8f4XvBlc8BmWpRugWEY+Ozy4Zejv2C1WelaqSvDGwyPuy+WLmDL5KxWSQY/fTp+c3KSD48D0eBXKaXsJSREavHOmQOvvgoDBkgx1MKFNa0hkxq/3j8u8I0VHm1l/Hp/hwh+Y2KkgMhXX8HBg1CqlFRxePllcHOLP873gi9t5rchMiYS87//8mbPy/AGwxlabyhl8pWx35tQaSs8HM6cuTfAjd0CAqS7SSxnZ/D2tttQE6PBr1JKpbetWyXgXbZM0hmqVInv85o/v1ZsyMQuBYUna396iYyU7tfjxkkMU7my3O7dW+KXu1379xof/fUR4THxY+5SsQsLeywkl0uudB65ShORkRLgnjghTXJif548CYGB9x6bJw+UKwc1a0rJj3Ll4rcSJcBisc97eAgNfpVSKj2EhUGu/wKDd96BY8fghRdkSq1BA+26lkUUc3cjMIFAt5i7WwJHp72wMJg1C775Bi5dkkm65cuha9cHG1P4XfFjys4pLDq0iChrFE6GHJDdkp33mr6ngW9GY7XKTO3dwW3sz3PnwGaLP7ZAAahYEdq2hfLl44PbsmXBwyPD/fulwa9SSqWV6Gj4/XeZ5d20SXLh3N2l9XCJEpAzp71HqNLZyA5e9+T8Arg5WxjZwStdx3HzpnRf8/GRXimtW8tMb5s298YxVpuVVf6rmLJzClvObSGnc05eqf0Krzd4ndvht3XxWkYQFib16Y4fv3c7cUIW1MbKk0cC3IYNpYRixYpSWrFCBWmTnolo8KuUUqnt0iVpeTV/Ply5AkWKwNChklAJUqpMZUmxeb32qvYQGCizvLNmwb//ygzv6NFy8eFut8NvM2f/HKbtmsa54HN4unsyod0EBtQZgLure9xxGvQ6CNOEy5flitL9Qe7Fi/HHWSwyY1upkqwr8PKSILdiRShYMMPN4KaUYZpmur2Yt7e3uWfPnnR7PaWUSjcRETKFVrSofOFUry5fLgMGwBNPPJg4qVQ6OnVK8nnnzZOr3b17S9W8atXij/G94Muyo8sICApg/en13Im+Q4vSLRjeYDhdvLpgcYrP3cxI5doyFdOUk+ujR+HIEdmOHpUtKCj+uNy5JcCtVEkSuGN/L1cuSy2kNQxjr2maD6y405lfpZR6HEePwnffySxvmzawdKl8yVy5IrlwStlBbHAacCIbUXu9uHW4EM7OBgMGwMiRUOauYgw208aUHVMY8ecIbKbkeXaq0InPW39OrSIPltnLKOXaMrTYIPfuADf2Z3Bw/HEeHlC1qpzNVKkiW6VKchKeRWZxU0KDX6WUSokVK2DiRNi2TWZ1n35a2g3H0sBX2cnK/YEMn3KB6/9UJfx0YQyXGPI1OMv4MW683L5o3HH/Rv3LgoMLmLJzCsdvHI/bbzEsNCnZJMHAFxy/XFuGExIChw49uN09k1uwoAS2L7wgP6tWlZ8JdRpRj6TBr1JKJdXhwzKrki0b7N0L167B+PHQr59+CSm7M01pDth/SA5CzjbEyS2KvM38yV0nAItrDHP2SfB7Pvg803dNZ/a+2QRFBOFdzJuPm3/MuO3jiLJG4WJxoaVny0Rfx1HLtTm86GhZeHZ/kHvuXPwxefJIylSvXpKTUq2aBLkFC9pv3JmQBr9KKfUwd+5II4rZs2HHDli9Gp56Cj78ED77TC8tKruzWuVCxNixsH8/WHK7kq/NEXLVuICTi8zQmpicCdlHz2XfsvzYckxMelTuwfAGw2lcsjGr/C5R3vDgQtQeSubw5uqNUlAy4ddztHJtDunWLThwAPz8ZDtwQFIWYhtAZMsmJ9KNG8sVo+rVZStVSv9NSQca/CqlVEJCQyXAnTdPcuwqV5Zl8o0by/2urvYdn8ryoqJg4ULpxnbihFSkmjMHvru0g8thd4h0Oka40wFMrERY9hDldJI/z7jzdqO3GVpvKKXdSwN35/CWJS9lCQnhoTm8jlKuzSHYbHD27L1Brp+flDWMVbSoNIDo0AFq1JAgt1KlLLXwzNFo8KuUUrEiI+WyZI0akCMHrF8vFRsGD4ZmzXRGRjmEO3fg++9hwgSJsWrVknWW3btLJav8+ysybPkcrlg+BaxgQDazIIOqj2Vi5zfI6XJvfenk5vDau1yb3URHy6Kzfftk8/OTHtChoXK/k5MEtc2ayR9KzZqyFS5s33GrB2jwq5RSZ89K4dMffpDbFy5A9uySj6clypSDCAqC6dNh8mS4cUNirNmzZUIx9rzs0NVDrL0whcvOc8GMDWgNnqv8Et92H53g86Ykh7db7eKZO9iNiJC//7GB7r59EujGNoXIlUsC3H795GetWrIIzU1TPzICDX6VUlnXnj3wf/8nM7yGAV26wJAh8QGvBr7KAVy9CpMmwYwZMsn45JPSmKJpU7nfZtr47cTvTN45mQ1nNuCWzY2nKnZm/en1RFujcbG4MLRxt0SfP8vn8N65I7O4dwe6R47EN6Vxd4e6dWH4cKhTR7by5R/s/6wyjEcGv4Zh/AB0Bq6Zplntv33jgaeAKOA08JJpmkGJP4tSSjmIy5clT694cQgPl9mcjz6CV16RlsNKOYiAACkmMmeOXHF/5hkJemv9V4EsLCqMeX7zmLJzCidvnaR47uJ82eZLXq3zKh45PPC94Juk9sNZKoc3KkpmdHfvjt+OHJF/EwAKFJBA98kn5WedOuDpqSlPmcwjO7wZhtEcCAPm3xX8tgc2maYZYxjG1wCmaY561ItphzellF2YJmzdKteMly+X1dXTpsn+mBid4VUO5ehRWcS2aJFMLvbvD+++KwvaAM4FnWParml8t+87giODqV+8Pm81fIselXvgbEnZZzlTdmyzWqXb4t2B7oED8akL+fNDvXqyeXtLsFu8uAa6mUiKO7yZprnFMAzP+/b9cdfNHcAzjztApZRKE3PnSjOKw4fl8uUbb0hqA8iXnAa+6jGkZtC4axd8+SWsXCnrLd94A95+Wy5IbD+/nYlr53Pi5gn+Pvc3BgbPVHmGNxu+ScMSDR/7fWT4HF7TlFz9nTtl27VL0hf+/Vfuz51bgts33ogPeHVGN8tKjZzfl4GfU+F5lFIqdZw+DWXLyhfb9u1SUmjOHCkcnyOHvUenMonUaPNrmrBpk9To3bQJ8uWTLJzXX5cr8FHWKD7Z/CWf/P0JJnKltk/1PoxtM5aSeRMpxJsVhIZKzn5ssLtjh7QUB1msWrs2vPRSfKDr5aU5uirOYwW/hmF8AMQACx9yzEBgIECpUqUe5+WUUipxMTHSgGLGDNi4EXx9oWFD8PGRL0Od4VGp7HHa/Nps8nEdO1auxhctKvm9gwbJJOXNOzf5cutspu2exqXQS3GPsxgWqhSskrUCX6sVjh2LD3J37rw3T7dCBWjbFho0kL/zNWpoDV31UCkOfg3D6I8shGtjPiRx2DTN2cBskJzflL6eUkolKDRUAtxvv4WLF6VD0tixshobtBmFSjMpKREWHQ1LlkhO79GjcoHi228lr9fVFY5dP8bkzZNZcHAB4THhtCvbjrcbvs3//fV/SWo9nCkEBUmQu327nMTu3BlfSzdfPqhfX4oaN2ggv3t42He8KsNJUfBrGEZHYBTQwjTNO6k7JKWUSoJbt2TBisUindfq1JFFbJ06SetQpdJYckqEhYfDjz/CuHFw7hxUqybd2Xr2BIvF5I/TfzB552T+d+p/ZLdkp0+NPrzZ8E2qFaoGQOOSjZNUuSHDMU1pT+frK8Hu9u1yVmCakqZQowb06SMzug0byiyvXsVRjykp1R4WAy2BAsBV4GNgNJAduPnfYTtM0xz8qBfTag9KqccSFQW//iozvTdvykpuJye4fVtmhJRKR/fn/IKUCPuye/W4tIeQEJg5U+r0Xr0q8dv770P+mr5sPPsn/0b/yxr/NRy7cYwiuYowtN5QBtUdRMGcBe31ttLWnTuS5xEb6Pr6yt9lkAWpjRpJC/HGjSVXN3du+45XZWiPU+2hdwK756TKqJRSKimuXZMObDNnSp3eChVg2DDJ83Vx0cBX2cXD2vxevw5TpsjFiOBgaN9eavS2aAFrTqym5bwexNikiULF/BWZ320+Pav2JHu27PZ8S6nvxg345x/Ztm6VCgyxzSO8vKSxTGywW6mSLkpT6UKvDSqlHJfVKmkNW7bIEviOHaVqQ4cO+iWpHML9JcLOn5dGYN99Jx1yu3eXoLduXdh3eR/9Vk5i0aFF2ExZrOVkONG/Vn/61uxrr7eQekxTOnNs3Rof7B4/Lve5uEh+7ogR0KSJzPBqrq6yEw1+lVKO5e7UhieflPbD3bqBvz9UrGjv0SmVIH9/+PprWLBAbvfpA6NGQYWKVtacWMPbcyex5dwWcrnkokelHqw5uSau9XArz1b2HXxKWa3SLS020P3nH7j0X2WKvHklyO3XD5o1kyYSuvhUOQgNfpVSjuHGDUltmD49PrWhdGm5L1s2DXyVQ9q3TxpT/PqrVNQbMkQmN/MXCePH/T/SedoUTt8+Tam8pZjQbgID6gzA3dU9ya2HHUpMDOzfD3//Lds//0hlBpDOaM2bS6DbtKms6NOrM8pBPXLBW2rSBW9KqUT16CGthzt0kOvGmtqgHFRst+yxY2H9epnkHDpUPrYR2c8zdefUuNbDjUo04q2Gb/F05afJ5pTB5puio6WRRGywu21bfMmxChUkgbl5c9lKldIqDMrhJLbgTYNfpVT6s9ng999lRdDMmVCunJQ3AqhSxb5jUyoRpgnr1slM7/btUKgQvPUW1O3qy/Iz/7UeDvgbgB5VevBWw7dSpfVwuomMlLbAscHu9u1SnQGgcmUJdmMD3mLF7DtWpZIgxdUelFIq1YSFwfz5EvSeOCFfoGfPSvCrQa9yACv3Bz5QvaFz9eIsWyaNKQ4elEnOadOg34sx+Oz9mo7LPopbwPZ8tef5su2XlMqbATqaxsTIzO6mTfDXXzKzG/5f3eIaNWDAAAl2mzWTSF+pTEKDX6VU+oiMlEulV65I/c5Fi+CZZ8DZ2d4jUwp4sG7vxRuRDPnwNkMOFubKxWxUrgzz5sGTTwcz//Acqn/nw7ngc3GPtxgWqhWq5riBr9UKBw7EB7tbtsgJKUiwO3AgtGolwW7+/PYdq1JpSINfpVTa2bFD0hs++URWA330EdSsKWWOND9QOZjx6/0Jj7Zii7IQ5leKkN1lsYa5krN4CMuX56FG87NM2+3DsGlzCI0KpXnp5gz2Hsynf3+apq2HE5qNvru8WqJME44ciQ92N2+OX6Dm5QV9+0Lr1jK7WzCTNtVQKgEa/CqlUpfVCitXwsSJ0r0pb14YPBiKFpWl8Eo5qAuXYwjZW4HQvZ7YIlzIXuoG+Tvtx6nMNhZG7+CZ6StwMpx4rupzvNXwLeoWqwtAi9It0qxyw/2z0YFB4Yxefggg4QD43DnYsEG2jRvh+nXZX6aMLCpt1Uo2zdlVWZgueFNKpZ4DB6Sq/5kz8mX75pvw8suQK5e9R6ZUoi5dgm++gUlTY7BFZcOl4Vpc6q3Gkudfwi17iHI6QT7XfAyqO4hh9YdRPE8SZl1TSZOvNhEYFP7A/uLubmx7r7W09v7rLwl2//wTTp2SA4oUgTZtZGvVCjw9023MSjkKXfCmlEobly5BYKDk8ZYtC+XLS7X/p5+W7mxKOajTp2HcOJg7Vy5YNG0fhX85H656vEsUVjDAYhZgYPUv+KbzcHK65Ez3MV66L/DNHhNFncBjNP3bD5a/JwvWTFNOMFu2lLbfbdvKAlJNLVIqQRr8KqVS5sABmS5bvBgqVZLbuXNL4VOlHNjBg1K54eefZb3lyy/Dc4PPsPKyD9t3zQDT+t+RBr0qv8ys7u/bbazF8rqS++Qxmp/dR9MAP+pfPIJrTBQxThZo3Ag+/liC3fr1dfGoUkmkwa9SKnm2bZMv3I0bIWdOyed9802dZVIOb/t2qdG7dq1MlL4zwqRpL1/mn/yGNqskn7dt2bb8FfAXMbYYXCwuDG3cLf0Heu2apDGsX8/G3/6H641rAJzwKMWimh3ZXb4OnYb1onNTr/Qfm1KZgAa/SqlHi4yU68I5csiCmmPHZOps4EDIl8/eo1MqUaYJf/whQe/ff4OHB3z8SQylOixn9qFvGL96J/lc8zGqySiG1htK8TzF07/1cFSURObr18tg9+2T/fnz49quHfsq1eeTyOIcNHPF1x5OSrUHpVSCdMGbUipxwcEwa5Y0pXj9dXjvPSmMb7OBi4u9R6dUoqxWWLFCgt59+6B4cRg2Ihhqz2Hm/imcDz5P+fzleavhW/Sv2T/983lPnowPdv/6S+rtZssmZQDbt5f23nXqaN68Uo9BF7wppZLu4kWYPBlmz4bQUFkx3vC/Nq3Z9J8N5biiomDhQllz6e8PJRr50mHCcnKXPM/Y078TujmUFqVbMPWJqXSu2Bknwyl9BnbnjtTZ/f132U6flv1ly0q93fbtpeZunjzpMx6lsjD9FlNKPWjQIJmV6tkTRoyQGSil7OxhzR7u3IHvv4cJE+DCBahVCwZ++x1zrg7hYpgVjkH7su0Z22ZsXH3eNHfyZHywu3kzRESAm5uUHnvzTejYUaqjKKXSlQa/SmV1pilfzBMnwvTpULo0jB8vv2ttUOUgEmv28G+oQcDWYkyeDDduQNNmVvp9tZJN4ROZfdE37vEWw0JLz5ZpG/gmNrtbsaKcUD7xhHRTc3VNuzEopR5Jg1+lsqqYGFi+XAqd7t0LhQrBiRMS/FapYu/RKXWP2NbDsaz/unBpdxn6jy+ENRLadw6l8vM/svraZL44eZYy7mV4quwrrDszH5sZg0k2jKiqqT+ws2elfMRvv907u9u6Nbz1lgS8Zcum/usqpVJMg1+lsqKoKKhRQ5IiK1aURW39+umMlHJYsc0eYoLdCNlVlrCDJTGtTrjW2kvfd5ax4sIs/jgRTOOSjZnQfgJmuDcfrjhKIasXEU6HcLVVZ95mFyrlD0y4LXBSxcRIZYa1a2HdOjh6VPZXqKCzu0plEBr8KpVVBAXJl/ULL0ilhn79oHJl6NoVnNJp0Y9SKeQe6cHJP0vwb8g58FyES0cXLFX+R0SOzcw/bdKjcg/ebvQ2DUvIwswmX20iPNpKdiqT3VYZgHCblfHr/ZMf/N68Cf/7n/z9+f13+bvk7AzNm8Orr0KnThL8KqUyBA1+lcrsLl2Syg3ffiuVG+rXly/q9+3XtUqppNq1S8qV+a1sCJ5b4MW24BRNlAGG6ULnsi/h89QHeLp73vO4+9sCP2r/PUwTjhyRYHftWpnptdmgYEHo1g06d4Z27bQyg1IZlAa/SmVW167Bhx/CvHlyqbZnTxg1SmeolMMzTSl9O3asNBJ0LxBOp4/ms8XlfUJjov87yuBZr9f4ufekBJ+jmLsbgQkEusXc3RJ+0eho2LIFVq+WLSBA9teuDR98IAGvt7deJVEqE9DgV6nMJjQUcueG7Nlh1SoYMADeeQfKlbP3yJR6KJtN4s4vv5QZ34Ke12j7xXT2Z5vBuvAbeOX14mxQGFabFReLC2827Znoc43s4HVPdQgAN2cLIzvc1RI4KEjSGVavlgVrwcGSq9uuHYweLekMxbWTmlKZjQa/SmUGpgmbNknUcP06+PlB3rzSilgX3igHFxMDS5bIx/foUShe6yhNvp7EnqgFbIiO5KkyT/FOo3doXro5Oy7uSFLr4di83gfqAueLhqlT5cTw77/lxQsWhB49oEsXaNsWcqZztzelVLrS9sZKZWSxU2Vjx8Lu3VC0qJRXevNNWZCjlAOLiIAff5RqewEBJp6tN5G340QO3Pkd12yuvFjzRd5s+CZeBbwe/WSJMU0p5bd6tQS8Bw/K/sqVJdjt2lXy4LWNsFKZjrY3Vioz+uUXeO45qSM6axb07y/pDko5sJAQWX/5zTdwNfsW8j8xjaKl/QiIOkkhoxCftvyUIfWGUCBHgZS9QGz+7sqVsl28KLm6zZpJM5enntLcd6WyMA1+lcpIIiJg7lzIkUNKlXXrBj//DN27Qzb966wc2/Xr4OMD06ZBUEQQJfp/AEVncgsTI8pgdJPRfNTyI1yzpSBV584d+OMPWLEC1qyB27el2USHDvDFF5K/6+GR+m9KKZXh6LelUhlBaKjM7E6cCFeuSLDbr5/U6+2Z+KIfpRzBhQvy0Z09G8KzB1D+hSlEFfmei9awuGOcDCdyZ8+dvMD31i0pRbZiBaxfD+HhkC+fzOw+/TS0by8nikopdRcNfpVydPPnSw7v7duyGGfRImjZ0t6jUuqRTpyAr7+GBQvAVmQ3JYdO5HzuZQQYTvSq2ov8TnWZ5vceNjM66e2HAwMl2F25UtoJW61QooRUNXn6aUlt0Hx3pdRDaPCrlCO6dElmdQsUgMKFpV3q6NGyMEepdLByf+CDlRKS2Blt3z6p3PDLrzacq66l4LsTuOS8lVvZ8zCi7gheb/A6e04bjF5+iELWzx/dfjggAH79VXLcd+yQfZUrS93qbt2k/q5hpO7/AKVUpqXVHpRyJOfOydL3OXPg9ddh/Hh7j0hlQSv3ByZYI/fL7tUTDYBNE2au9WXyys2c3NgI15LHydF6ErecTlAqbyneavgWA2oPIHf23IC0H06oCUVxdze2vddapo1jA959++TOOnWkJFn37lCpUuq/caVUpqLVHpRyZCdPwldfSYqDYcCLL8LgwfYelcqixq/3vyfwBQiPtjJ+vf8Dwa9pSn+I0dN9OVSnNZSMhJdMIoBqxbyZ0WgJPar0IJvTvV83D7QZNk0q3jjHE/9sh4XD4fBh2d+woZwE9ugBZcqk9ltVSmVBGvwq5Qg++khyGIcMgZEjoWRJe49IZSLJTWF4IDBNYL/VCsuWSXrDwUB/XJ4bDs4Rcfe/WudVZnWehZFIOkIxdzcCb9+h6rUzPOG/jSf8t1Pu1kVshiF5u1OmyAxviRIpfNdKKZUwDX6Vsof9+6X80pgxUK2arAqaNAmKFLH3yFQmc38KQ2BQOKOXHwJINAAu5u6WYEpCMXc3IiPlAsVXX5ucif6HXO0nQInVmE7OWExpFOFiceGlWi8lHPiaJhw4wA8nVuC68ldK375MjOGEb6kaLGjQjQZvvcwT7Wqn0rtXSqkHafCrVHrasQM+/xzWrYM8eaBXLwl+S5Wy98hUJpWcFIZYIzt4PZDzm910oeL1upQpF8Nl9+Xk7DwR8u0iu5sHb9f7iKH1h3L61umEWw+bpnRWW7pUpotPnsTLrR9ZDgAAIABJREFUYuFqvaZ8Vaw3S4vXxa1YEUZ28OKJJC6qU0qplNLgV6n0YJqyKn31aim0//nnMHQouLvbe2Qqk0tKCsP9YoPi8ev9uXA5Bo6UJ/BAAeaWm4trn0ngFkDR/OV5u+EM+tfqTw5nqaVbKGeh+KDXNCVvd+lS2U6ckBbCrVvDu+9Ct24ULlCA94D3UvUdK6XUw2nwq1RaMU3YuRMaNJBFbPXrQ/PmMGgQ5Mpl79GpLOJhKQwPU79wcUrfOs+uE2uIyf0TzoPXQLYg6pZswojGk3iq4lNYnCz3Psg04ciR+Bne48elrXCrVvDOO1KHt2DB1Hx7SimVbBr8KpXaTFParH76KWzfDhs2QJs28MEH9h6ZyoISSmFwc7YwsoNXgsefPv1ftb1dC7F26Q8NrGBAk9ItGdtm7L3pDLFOnYIlS2DxYjh6VALeli1h+HBZtFaoUBq9O6WUSj4NfpVKLaYJv/8uQe/OnbJKffp0aNLE3iNTWdjdKQwPq/Zw8CB8+ZXJzzs3Q+PxmN1+j7vPYlhoX679vYHvhQsyw7t4MezdK/uaNZPPfI8e0pxFKaUckAa/SqWWyEhpsZo9O8yaBf37y+9K2Vm32sUTXdy2fTt88WUMvwX8glPT8Zj99lHArRBPV36VBQcXEG2NxsXiQkvPlnD9uqQzLFkCW7fKE3h7w8SJ0LOnliVTSmUIGvwqlVI2G6xaBfPmSUDg6iopDhUqSGtipRxUbGbOZ+NC2XZnDk6NJ4P3Ocq6e/Fu09n0rdkX12yuvFTrJTYf/52WZ2w0emUMbNwoBX6rVIHPPpNqJeXL2/vtKKVUsmjwq1Ry2WywfLl8+R88KF/+587Jz6pV7T06pRJltcKKFfDpN5c5lMMHo8G3kD2IhiWaMaqpD50rdsbJcILwcFixjEaLFtHot98gKgrKloVRoyTgrV7d3m9FKaVSTINfpZLj0iVo315WtHt5wU8/wXPPQTb9q6QcV1QUfPKDLzP+/pkg4xS0/QPDYuVpr+6823QEDUo0kMh44yZYuFBO7kJCoGhReO016N0b6tWTqiWJSG4XOaWUshf9xlbqUWw2qVFaqZJ0YKtcWSo39OwpdUuVclB37sB335mMWTeDoEZvQCUbAD0q9eDrdl9TLl9Z2LcPvnlb8ngvX5bmKz16wAsvSMWGJHzGU9JFTiml7EWDX6USE5ve8MknsrI9IECaUixbZu+RqSzuUbOsQUHgMy2GCet+JbT6BGiyJ+4+i2Ghrqsn5aYvkllef39wdoZOnSTg7dQJ3B5eA/h+Kekip5RS9qLBr1L3uzvoPXxYZnxnzoTcue09MqUeOsvaqFhxvp70LzN8fyCy9iToeJaSOSrSu/a7TN3hQ5Q1CpcYk5bDJsJFoEULaT7xzDOQL1+Kx5SSLnJKKWUvGvwqdT8/P3j2WQl6Fy3S9AaV5pKTL5vQLGvoDRcGvnmNW9lnYq0zA1rfpmb+xoxp8RVd/MHJ5ye6HYxmc0kbLbOVo9HrAyWPt2TJVBl/SrvIKaWUPWjwq5TNJkvgT5yA0aOhTh0pWZbEfEelHkdy82VjZ1MjnY7xb/gJIo/VIsplPzRfAJYo2pXsxpjC7Wi8xg8+GiQ5EEWL0qjP2zTq2zdNKjUkt4ucUkrZ0yODX8MwfgA6A9dM06z23778wM+AJxAA9DRN83baDVOpNGCasHIlfPwxHDokQcGIEZL/2KaNvUensojk5ssWc3fj+OkT3Cr6AeSNgkaANRueES1Z7+ZFxbG/w5kV/D97dx5mc/nGcfz9ndVYx5bsa5aQvX6T1CCNPUubEpIUpUUUEkKIEhUJpaSSSrRYsySyZA2Jshv7MgxmzHK+vz9uS8o+Y87Mmc/rus5l5syZc55zXWPmc55zP/dNxow2WrhVK6hV67q+kLvSKXIiIqnBlez8fgy8B4z/x3XdgDmu6w5yHKfb6c9fTv7liVwnK1dCu3awahWULGktyx56SDu9kuKutF7WdeGnuYlE/nqAw8UHgn/c6S9A29+z8eHUn8CZYy/ceve24Js58/Ve/lmXmiInIpKaXDb8uq67wHGcIv+6+l4g/PTHnwDzUfiV1M51rXdptmyQIwfExsL48Vb7qD694iWXq5f1eOCbqbG89Nl4tuV7E6r8RZa47MR4HFzXJcgDLfZmgjdegocf1ohhEbkufKmX97X+xc/juu4eANd19ziOc0MyrkkkebkuzJ4NvXrZifbp06FIERtUcYmm/SIp4WL1sp1rl2LUJ4fp/f377C/6DpTfz83Reeg1KxvNlxzht5tDmd/gZsLrPEFYv9b6WRaR68bXenlf9+0ux3HaA+0BChUqdL0fTuQc14V58yz0LloEhQrB44/b9Y6jsCCpwr/rZW/MlIliMTlo++ZgjpX4EMqf4O5dWenxNYTvOojTqDF824awevUICwz08upFJD3wtV7e1xp+9zmOk/f0rm9eYP/Fbui67mhgNEDVqlXda3w8kas3bpyF3fz5YeRIaNsWgoO9vSqR/2hSKT+ZQ3bQ6+tvWbHtT5bmn4tfKZdH10PXhVA+f3F4oY2V6OTO7e3likg642u9vK81/H4HtAYGnf53arKtSCQpfv0VEhLgzjttRGtMjAXgDBm8vTKRCzpwwOWRd99httMZ/DxQAFqshUGrclCoSWvo1RoqVPD2MkUkHfO1Xt5+l7uB4zhfAIuBUo7j7HIc53Es9NZxHOcvoM7pz0W8Z8UKqF8fqleHPn3sumzZ4OmnFXwlVdq6PYHGL4zl5l6Fme3/PDgeAPyB8hGPUmjjXhg6VMFXRLyua0QpQgLP74aUlnt5X0m3hxYX+ZIaoYr3rV8Pr75qQypy5IBBg+CZZ7y9KpGLWrUuml5v9WBJ6DgOhp6gzAFouzYX75Y/SpzjISggiPB6HazftIhIKuBrvbzV30nStmXLYM4ceO01eP55yJrV2ysSuaCff1jDsAnPM6/IAo4W8XD7Tn9GnIjgvnav4VftVprsWsL8bfMJLxJOWMEwby9XRLwktbYU86Ve3o7rptwZtKpVq7rLly9PsccTH7R5M/TtC1WrQqdOVt977Jjt+oqkMu6pOGYMHsPY3wczrdQOTgXA3Tty0aV8B+55sptNYRMROe3fLcXAygsGNiuftoKn60J0tI1Xz5/fawOkHMdZ4bpu1X9fr51fSRt27ID+/a2DQ0CATWUD+1jBV1KZX3/5gpFfDmLbiT/4tXACgaUhYn9F+rR8g8o17vH28kQklUp1LcUSEizEZs9un0+bBtu2WaiNioIjR6yN6Kuv2tdr1oTff7eveewcA7t2WQBORRR+JfV7913o0sVeST71FHTvDvnyeXtVIuc7eZLEiZPoP7UPfStux5MbyAU1Y2ox9vlPKJZXk9dE5NJSpKXYyZOwe7dd9uyBuDh49FH7Wq9e1h//wAG7HD4M5ctboAV4/XXrqgR2mDx7drjjjnP3HRYGN99s14eG2r9ZsiTf2pOJwq+kTocOgZ+f/ccpWxZat4aePe0Vpkhqsno1saNG8fmK8bxdLYZ1lYHT1WT+fv7UaXi3gq+IXJEktxQ7cAC2bIHIyHMB99Ah+OAD+/rjj8NHH53/PTlynAu/x49bicItt1hP8dy5oWjRc7f98ks7jJst24U7KQ0YcGXr9DKFX0ldjh2Dt9+Gt96Cdu2s1VOtWnYR8YILHj4pkRUmTuTo2FGM81vJW2GwqyFkjyrKfbla8mPUm8QlxhHkH0R4kXBvPwURSSMuNu78bEuxPXtg9WrYvv3cZccOmDEDMmeGIUPsckZAAOTNC7GxFlYbNYKbbrJ3T/Pls6/dcMO52w8deukFFvCNF/IKv5I6nDwJI0bAG2/Yq9RmzewVqogXnXf4xHXJ9cdq4r8cwq4dPzOqYgzv1vTnWAjkPFKd/mV60P2+eny3Zjcbp+dmZ9xyCmasyr6DhaCgt5+JiKQFTQpl4IY8B1gy6zey7dlBiZhDVEg8Smirb+wGX30Fzz1nHwcEQMGCULiw1eVmzmw7uHfdZcE2f37IlcveRT37AE1S/kmlQgq/kjq88AKMHg0REXawrep/DmeKpLghMzcScPwYj66fx8OrZxCYuI1Bt/vzcEOXBH+HvFH3MjT8JR6vexvwz7BcjGwU49gx6D55LUDaOqktItfHqVOwaRP8/bd1Lzrz72uvwe23w4IF3N6pFbeDdYMpXBiKFIb4ePv+Zs2gShW7Pm/e/3ZRKF/eLnJJCr/iHQkJ8Omn8L//QZky0LUrtGwJNWp4e2UiZvlynv3iDRpt+JmJZU/R6P5gNucGEv0J2H4f3z3Vm0a3lzzvW1LdSW0RSXmJidYR4c8/YeNG+/e+++Cee2DtWqhW7dxtc+aEEiUg5nSd7113wcKFULw45MkDjnP+fRco4DOlB96k8Cspy+Oxt21697ZfCi+9ZKUOJUrYReQ6umzz+OPH4YsvYNQo3JUrCSkZyM3tMrM9zylwT4EbQPZjvShzc9h/gi+k0EltEUkdjh2zv2MbN1r5wV13wcGDVm4QF3fudrlyQeXK9nGZMnZorHhxu4SGnn+fOXJA9eop9xzSKYVfSTkzZsDLL1vLlHLlbCTxvfd6e1WSTvy7eXxkVMy5kgT/QzBqFEyYQMKJaCaEF+DV5wqwK/suiE0E1wHHBTwEZPubrhFtLvgYST6pLSKpz/Hj1vKrUCFrudm4MaxcaZ0UzjhTa5szJ7z4om3mlC4NpUrZdWdkygQPPJDyz0HOo/ArKWfuXDhxAiZMgIce8trEF0mf/l2SEBx/ivprF1J8/Iuw4w9OZApiWNMKvJF3F9GZduF/uAxNnXHc/L+sDFz9MB43Hj8nkM53Nb1oCcNlT2qLSOo3dSosWQLr1tll2zYb3jB3rpUhZMoEd99tu7ilS9ulWDH7XsdJM+2+0jONN5brZ8UKeOUVO8wWEWHBNyjIegSKpLCi3X7EBYod2sUjq6fTfN0cQmOPszzvjXz4QAU+zPgb8cGHCdp7O48We5m3OzQkS2Y7Jb1452Lmb5tPeJFwwgqGXfJxLltaISLe5brnWoaduRw/btPLAOrWhTlzbNe2XDm7VKtmf8ckTdF4Y0k5GzbYqMNvvrG3ew4ftuszZfLuuiT9Skigxa7l1F/4Lf6Ja5hT1I/Bt5diYqG8bM2/BAJnknFXIzpXepnXulUnOPj8bw8rGHbZ0HtGk0r5FXZFUouEBKvJXbMGWrSwndlOnay15hnFikGlShaKHQfGj7da3KAg761briuFX0leL78Mb75pLVp694bOnSFrVm+vStKrPXtg7Fj44AMGREbyXalQmj/gR4KfB9gAnk1k3vYIfeq8xPM9y6oSR8QXLF4MH35oO7rr1ll7MbDRu0WLWq/bUqWgYkWbZJYt2/nf/8+hD+KTFH4l6fbtsxOqgYF2evX556F7dzvhKpLSXBd++cV2diZPhoQE3HvqMPO1jrTc9hEJ/lGnb+cQlulpFn0y/D/dhEQklYuNtd3cFStg+XL794MPrH1mZKTV7VasaLu8FStChQrWkQGsXvfuu727fvEqhV+5dlFRNkZx2DC7PPEEtG/v7VVJehUdbb2jR46E9eshNBRPp2eYcNdNdF/5Kbt3vQLxoTh+gTh+HoIDg3jr4YcUfEVSu9hY64+bM6eVKKxYYSE3IcG+nivX+YORmjWD5s3/2yNX5DSFX7l6J0/Cu+/CoEEWgFu0gPBwb69K0qv16y3wjh9vh1YqVyZu7AcMz+th4OLhHFk9DI4WpXLsCN5r9xjkWX3Fh9dExAvi42HiRFi61LourFljQbdHD3j9dStZeOklm3RWtart6P4z6P5znK/IBajbg1y9iAiYNQsaNLBfRBUqeHtFkt4kJMB338F778G8eRAcDA8+SPSTbei3bwUjVg7jZEAk7K1IeMDLjHj6Pm4urdf6IqnOoUMWcpcutZ3dZ5+1YUjZs9u/1arBrbfav2FhkC+ft1csacjFuj0o/MrleTw2kaZePTsBu3ChXX/HHd5dl6Q/Bw/aAbb334cdO6BQIRY/1YAfygezZu8RZm2fSnxAFH7ba9I458u88+w9FCyotz5FUoUz3RTAdm4nT4bNm+1zPz8rV/jqK/t82zbb0dUpVEkCtTqTq+e6MH26vdW0Zg0MHWo9exV6JaWtWmWlNp9/bie3a9WCYcP4qngcD33bEs9yq/3z3xPO4/kGM2hwNZ23FPG2I0es88Kvv8KiRbBrF2zaZAE4Lg7Kl4d27eC226x8IUuWc99bpIjXli2+T+FXLmzhQuvYsHChHTD47DObyiaSUuLjrVf0u+/aH8+MGeGxx+CZZ1iVPYHnvnqDX9ZMBFxwwMGfni3uoU/tat5euUj647rw99/WSiwgwErieva0r/n7Wx/dBg0s9AYH2yFpES9R+JULGzjQfpGNHAmPP65m35Jy9u6F0aNh1Cjr01u8OAwditumDfMO/84Lk7vy+8npcCoLmXa1IK74ZDxOPEH+QUSUDPf26kXSh/h467rwyy+2q/vrr3DggLUdq1IFatSA/v2henWr19WQI0lFFH7FbNkCffrYYIrixWHMGKvvzZjR2yuT9GL5chg+3OrL4+NtxOjYsXgi7uHr9T/QbWR9tiYsgeM3kGfrAPrd24HH+oby254rHz0sItcoJsYOpRUoACVKwPz5cM899rUSJaB+fQu6hQrZdXfeaReRVEjhN73btw/69bPm4IGB0LChhV+dqJWUkJAA335roXfRIsicGZ56Cp55hrjiRfhw2ef0fr08B9w/4XAxiux+n8GPtKb5GyH4+cGUVZEMmRnD7qhb+CE0hq4RkRotLJIcEhJg9mzb2V2wAJYtsxelr74KffvC7bfb4bQaNSBPHm+vVuSqqNtDetavH7zxhjUQb9cOevVS6JWUcfiwdW147z3YudPqyjt1YnFEWWbuW8Tuw0f48vfJHHN2wd4KlD/ajaHt7qN2zYCzh8WnrIqk++S1xMQnnr3bkEB/BjYrrwAscrWOHIGff7ba3aZNLfxmz25/H6pWPbeTW726vSsokgao24OY+Hjb4QXb9W3QwOqybrrJu+uS9OGPP+Cdd2wgRUwM1KxpB9oaNmT6llk0/qI+CZ4EcIA9lageO5a3n7mHatX+265syMyN5wVfgJj4RIbM3KjwK3Il5s6FH3+0XtmrV1vwvfVWC78BARaGS5VSva74HI1BSS8SE2HcOCtpONOn9513rL5SwVeuJ4/HWuZFREDZsvDxxzYVcM0amDuXHeGVaDWpMw0mNCbBPR18XT+ej7ifhZ9EXDD4AuyOirmq60XStRMnYOZMO8x8xgcf2LsvWbPaeY+ff7YShzMqV1bwFZ+knV9f57owdSq88ortulWtem7nVyMg5Xo6cQI++cTqeTdtgrx57V2G9u0hd27+OPAH3T5qww87PsP1gN/We/ArPhf84gkKCOKBauGXvPt8oSFEXiDo5gsNuU5PSCSNWb8evvjCdnaXLbNShoAAaNvW6nSHDbMXoyH6P3M17KzBRnZHxZAvNISuEaX0blMao/Tj6xo2tLewEhPh66/tF+Btt3l7VeLLdu+2wSgFC8LTT0O2bNYnets2eOUVlpzaTI0RTSg7sizfb/6KwJXP8HTCFj5sPZabgwaTJe4RSjhvsO9goUs+TNeIUoQEnj/9KSTQn64Rpa7jkxNJpTweWLkSBg+27j1gpQyDBtnv/y5dbOc3KurcAbW8eRV8r9KZswaRUTG4QGRUDN0nr2XKqkhvL02ugnZ+fdH69VCmjO3sNm1qlzZt7BW/yPWyahW8/TZMnAiJiewOj6DvTXWZma0YebeHUHvuJKb+MZa10T/DyRxkXNebztU70eXTnMzbEkn3yWuIiS9GNopx7Bh0n7wW4KI7Kmeu1w6MpFvR0fb/bfZsq989dMiuz5PHDpE2bQqNGllZgyQLnTXwDer24Eu2brWODZ99BhMmwMMPe3tF4us8HjswM3So9f3MnBkef5xZtR/guWXHOJL4O9F+PxKbuA1Phh1wtAChG16kZ/12dHg889k20tUHzb1gCUP+0BAWdauVss9JJLU6dMhCbpYs1gc7Kgpy5rQd3Dp14O67oXZtuPFGb6/Ua653SULRbj9yodTkAFsHNUi2x5HkoW4PvuzAAaulfP99GyP50ktQr563VyW+7ORJq+cdNszqeQsWhCFDrGVeaCi9B84gMnEs0UGTwXEhwMFZ1pVimR5mw3cVz5adn6HDayIXsWgRTJtmJQsrV9o5jgYNLPyGhsLmzVC4MGd7AKZj/25/eKYkAS7+DtLV0lkD36Ca37TOdS3ojhhhpQ1//201Xtmze3tl4ot277bDkwULQseO9nbqF1/YH+AuXTgW4kf/eUNYGv0o0cHfwNk9Ej+yVtxHYrHI/wRfuPgfDv1BkXRnyxb45ptzn/fsaf3YM2SwKZy//gpTppz7epEiCr6nXaokIbnorIFv0M5vWhQXZ23LWra0NjTDhkGuXFC6tLdXJr5q3Tp46y0rqUlIgCZNoHNna3jvOBw4cYBBs4YzcvkIYomCXXfjv782ibf1AeJx8CfELX/RMNs1otQFB1boD4r4vOPHrWRoxgzb3f37bzuvceiQ7eyOHg25c2uwxBVIiXeQdNbANyj8piUejx1uePVV2x0IDrbd3jvu8PbKxBe5rrVIGjLE/jBnzAhPPgnPP2/9ooHtUdt57ac3+XTdhyQQCxuaUcN5mfp1i/DJ5hVExfcn1m8tGTzlCfUvd9Ewqz8okm64rr2YLFzY3jkZM8ZeSGbMCOHh0KmT9cTOls1urz7sVyylShKaVMqv301pnMJvWuC6MGsWdOtmrWsqVDg3NEAkucXHW1u8N9+0GsMbbrCa8qeessM1wB8H/qDnjDeYuvlzPB7g90dpEPoSA3uWpnx5u5vSq8ozZGYQu6PKXFGY1R8U8VnHj8OcOVa7O20a7Npl76I8/DA88ADccottYgQHe3ulaZreQZIrpW4PaYHrQq1a1ie1f3+bjqUBFT7Law3Uo6Nh7Fgro9mxw8aaduli5TUZMgAwdsVYBv48lC3RGyAuI36r29OicGdee7Hgmc1gEXFdOxSaKRPs2WO7vPHx1qWhTh2oX996sJ/ptyvJRgMo5J8u1u1B4Te1+usveO01O7xWoIAdNMqZUzsDPu7fp5XBdi4GNit//X6B795to65HjYKjR6FGDeja1U6U+/nhui4/bfmJZ6Z2Y1P0SjvD5gmgedx3vPtcPfLmvT7LEklTTp60MqHp021397bb7DAowIABEBZmNfJBQd5dp0g6olZnacXevdCvnx1yCA62t8QKFIB8+by9MkkBKdpAff16K2347DObANWsme30np4A6HE9TP7jG3pMG8hfJ1ZAbBYIdsBx8Q9wqVJnNXnzqqWeCM88Ax9+CLGxVrtbu7bt7p7Ro4f31iYi/6Hwm1q4LvTta4eLTp2C9u3tYFs6blaeHqVIv9uFC6110g8/2GjTfx1ii0uM49M1n9Fr5hvsjtsIh0oQun4MYVXyMtNpjseNxyUAJ65s8q1JJC3weGD5cvj+e/j5Z6vjDQyEokXtd3aDBnDnnWfLhEQkdVL49bbERBtM4Tiwfbv98uzfXyd806nrdlrZ47Gw+8Yb1ic0Z07rGfr009YmDzgRd4IPfhvL6/Pe4nDiTthbgTx/TaTfQ/cR+sheev2wlhtOneve8Mn8IErniFQ9nfi+tWutFv7HH2HfPjtzUb26fVygALz4ordXKCJXQeHXW/7ZtmzSJKhSxVre+Ptf/nvFZyX7aeW4OPj8c3tH4Y8/7ODNO+9A27YsPvw78zeMoXLeyvy6/TfeWjScE+5B2F6DopEfMOCxutz/noO/P1QfZOUYwZQh2FMGgBiP5tmLj9qxw3Z3q1eHihXhyBEbPFG3LjRqZIOFcuTw9ipF5Bop/HrDrFnw8svWtqxiRRsaAAq+knz9bqOj7cXU229bW6VbboEJE6yGPDCQxTsXU2t8LU4lnMI9M4VtUwPKR3Vn0NPVqVfv/KFRGj8sPs11Yc0amDrVpqetXm3X9+1rv6OrV7cx8hcaTygiaY7Cb0pyXTtUNGWKjaScMEFty+Q/ktTvdv9+29kdOdJ2q8LDLQRHRJxNs1uObOHpH54nNj4WHMB1KLy3E58+NZwaNS58t5pnLz4nPh527oRixaz8rFYtiIqyrgxvvAH33mvt/sA2JrQ5IeIzFH5Two4dULCghY8777RA8tRTalsmyWbWj0uI7j+QBstnEJQYz96adck3oPfZzg0A6/av49WZg5i6eSJuogN+9sc82D+IL155iLCCF79/NY8XnxAdbdMKp061+t1cuWDTJggIgMmToUwZ9d4VSQcUfq+n/fvt8NqoUfDll9C0KbzwgrdXJb5k/Xp2vNybWtO+xeP4MblcLUbf2ow9NxZmYFABmgBLdy2lx/SBzN09FeIy4bfyeR4s2Jl7W29nq2c+4UXCCSsYdsmH0fhhSfMGD7YzFnFxduDz3nvt4rq2MREe7u0VikgK0ZCL6+H4cRg61A4ZxcRAu3bQq5d69UryWbYMBg6EKVM4GZSBzyrUZWy1JuzLYp0bXFwyZ/2TjFmns/zQXIjJTsCKZ2lbrhM9O+ek4CV2eUXSvK1b4dtvbTf3o4+gZEkbPvHTTxZ4b7/ddntFxKdpyEVKcV1rcL5sGTRvDq+/fq5uTCQpXBfmzrVpUXPnQvbs0Ls31Y+U4khIVk75bSDGbw64fpxIXEZC/AbYlpcMq96k0+3t6fp5FnLntrvSCFDxOYcPw4gRFnjPHFirUMEOqpUsaR0a6mkoi4gkMfw6jvMC0A4beLoWeMx13djkWFia4rpWQ1a3rjU3f+01CA2F//3P2ysTX+DxwHffWej97TfIm9cms7VvD1mykHHQXPYcW86+oFeABDvEdiwf/r+8R997H+eZbzOQNeu5u/v3COXIqBi6T14LoAAsaYfr2v+qkAwTAAAgAElEQVSHuDi44w4rXejfH6pWtf8fTZvaYTYRkX+55jYDjuPkB54FqrquWw7wBx5KroWlGT//bCG3aVPr3gAWghV8Jani4+HTT6F8efv5OnQIPvgAtmyxpvpZshCbEEvZEsvY798fnNPB1+NHpqPN+fytJvR46fzgC5ceoSySqiUmwvz58OyzUKiQHejs2dO+lj27jYdftMj+fyj4ishFJLXsIQAIcRwnHsgI7E76ktKI33+H7t1h2jSb8PPRR9CqlbdXJanQVZcYxMbCuHF2QGfbNgu/n38O999/tk4x+lQ07y35gEEL3uKYZy8cvhly/g1+ifj5BdCzxd088L8LP4Z69kqa4vGcawfZvLm9y5Yhg20yDBgADRueu2327N5Zo4ikKdccfl3XjXQc501gBxADzHJdd1ayrSy1e+op2LDBAsozz0CI+p3Kf11VicGJEzB6tB2U3LPH3j14910beX26R++hk4d485d3Gb7kHWI4ApvvpvTBzxnYPpzclZewYPvluzeoZ6+kevHxVtf+9dc2lnvdOuvQ0KEDPPII1K8PmTJ5e5UikkZdc/h1HCc7cC9QFIgCvnIcp6XruhP+dbv2QHuAQoUKJWGpXnb4sIWSzp0hd274+GPrEakRl3IJlyoxOBt+jx2zoRRDh9rhnJo1rYSmZs2zoXd39G5enzuUMatGEe+cgA1NqHaqO288eyvh4WduFkb1QpduWQbq2Sup2N9/227ulCk2pCVLFhsnHB1t4TciwtsrFBEfkJSyh7uBra7rHgBwHGcycDtwXvh1XXc0MBqs1VkSHs87YmJs923gQDh6FMqVs52HkiW9vTJJAy5ZYnD4sE1jGz7cJkvVqwevvALVq7N452LmLxzETTlvYsrvPzHxz3EkugmwrgW1g7sx8MVyVKt2bWtSz15JNU6dsnHvuXPbOx2OY90aGjeG++6De+6xEgcRkWSUlPC7A/if4zgZsbKH2oBvNfH9+GNrir5rl73NNnAg3HKLt1claciFSgxynoji+bU/QOGHrCd006YWeqtUAWDxzsXUGl+L2IRY66OSGICz5nGa5n6J/v2KcfPNSV9XkkYoiyRFXBzMmWODf6ZMsU2Fli0t/BYvbu9+BAZ6e5Ui4sOSUvO71HGcr4GVQAKwitM7vD5j+nRrK/Xpp5r+I9fknyUGeaIP8uTSybRYM5MMnnh48EHo0cPeTTjtt8jfaDmpnQVfABwqx3Vh8rsDKVzYO89BJMnOTFEDuOsuWLIEsmWzF34PPGC90c9Q8BWR60wT3v5p6VILI++8A2XLWp1Z5sznfmmLXIOZPy4lpm9/6q2Ygb/rIbLBfRR+s9/Z0hnXdVmwfQFdvnud5Udmw6ksEBCD4+cSHBDE3NZzLjt+WCTVSUy0VpCTJtnhtXXrICgIvvkGgoOhTh37V0TkOtGEt0vZtMlC7zffQJ48sGOHhd8sWby9MknLtm6FAQOI+PhjewHV7nHo1o3CRYoAFnqn/TWdrt+/zobjv8LxPGT6/Q263NWB6k3Xsfzg5Ts3iKQ6mzbZBsLXX8O+fdaVoVEjO8CWJ4+1K5N0SZMlJbVQ+O3c2X5RZ8gAffpYc/TMmb29KknLNm+2E+vjx1t/0iefhJdfhoIFAUj0JPL1+sl0+3EA206thqhChP7xHq82bMtTr4WQMSNAGHVQ6JU0wHVh5UorYyhRAg4etL7nDRpYaU/9+pz+oZZ0TJMlJTVJn+H35Enry+s49u+TT0KvXrYrIXKt/voLXn/d2pQFBkLHjvDSS5A/P4t3LmbOgk84HhvLJ799zd6EjXCwJDf+PY5+DzxCqzcCCQry9hMQuQobNsAXX8DEifaz37EjjBhhB9f279cmgpznito+iqSQ9BV+4+JgzBjo29d25SIiLKyIJMXGjdC/v01hCw620atdu9phSWD+1vncM+Ee4j3xdvtDJSi2dRIDWzej+TB//P29uHaRa1G7ttXxOo71o375ZWjWzL7m56fgK/+hyZKSmqSP8Ou68NVXVte7ebOdNtYuryTVhg0WeidOtLKZzp2hS5ezP1vRp6IZ+ssoBizqSzyng6/rR5uKj/HR8Pt1jlLShv377ffnggX2s+44Nlq4cWPr1HD6RZ7IpWiypKQm6SP8Nm5sIzLLl4cff7RhAkoecq3++MPePZg0yWoZu3a14HvDDQAcjjnMgDnvMmL5cGKdI7CrKn751oJfAsGBQbSvU1M/fpK6HT8O334Ln30GP/1knRvKlbMgnCeP/cyLXAVNlpTUJH2E3xYtbFpQy5boPWa5Zhs2WOj98ks7wd6tm4XeXLkA2Hd8H71mDGXc2pHE+x2HTY25J6QHbz5/G8dDFzN/m7o3SCoWH28T1zJnhpkzoVUrKFLE6tYffvi8ftQiV0uTJSU1UZ9fkcvZtMlC7+ef207vs89aV5CcOQHYcXQH3b4bwqS/x5JIHH4bHqDZDd0Z9MItFC/u5bWLXIrr2sCJzz6zF3XPPQc9e0JsLCxfDtWr610yEUmz1OdX5Gr99Rf062fBIEMGe6u3SxcWx/7N/D9GUzx7cT5dMpMfd43H9UDA+lY8XrwbfYfdRL583l68yGX07w/jxsGWLfbzfe+9cMcd9rUMGc59LCLiYxR+Rf5t82YLvRMm2ESqzp0t+N5wA4t3LqbmJ7U4lRgLLpAYRPD6p+hYoSs9xxQiRw5vL17kIvbtg3nz4KGH7PMVK6BYMWvz2LQpZM3q3fWJiKQQhV9J1/45caiq5whvbvyBwj9+bX16n3vO6h1Pd29Yuus37p/wuAVfABxqZ3yRKeMHqLOTpE4xMfDdd9baceZMO7hWowbkz28T2HQGQkTSIYVfSbfOTBzKeXA3A379kvvWzcHj+LH5gTYUf6vf2RZO8zb/Qqev+7M+dhacygwBAeDnwc8J5O5b/6fgK6nTTz/ZKOFjx6BAAXsh17KlBV9Q8BWRdEvhV9KtcZMW0mPmeB5cMwvXcZhQqT7v33YfgQULsPDGG/l+/Sw6T3mdzQkL4PgNZP/rdYID7sQtsYRTAWvJ4CnPJ/ODKJ0jUieWxfs2bYJPP4WKFS30VqhggycefRTCw234hIiIqNuDpEP79sGgQZx6dwSO6/JlhXsY8b8H2Js1Fy4eTiauIDjLNCLd3+BoAYru6cqbD7fjzXVL2H30v03a84eGsKhbLS88EUn3Dh+2Lg3jx1vXBj8/G7TyxhveXpmIiNep24PI4cMwZAi88w7ExjK70j0Mqnofm3McIcZvDp44OOFZhifLRjhUjHJRoxnWthW17woG4MVFGs8pqYDrnms/1qCBhd5y5exn++GHUasREZFLU/gV33f0KLz9tl2io23oSe/eJJzIxM7JX7DXvzuQAIHAkcLctOlDPuvcimpVzv/vofGc4lW//w4ff2yT11avhmzZYPBgyJLFShzUj1dE5Ioo/Eqq9c9ODNc0Dej4cXjvPQsIR45Y/eNrr0G5csQmxLJq+sdE0h+cBLu960fDMg/z/bC2F7w7jeeUFBcVZXW8H38MK1daF5JGjez6bNmsc4OIiFwVhV9Jlc50YjgTNCOjYug+eS3A5QNwTAyMGgUDB8KBA/bWcN++ULkyJ+JO8OqkYbz/+xBiA3fjHCmLX64YHL9EggKC6FG/0UXvVuM5JUXEx1u4zZ0bdu+2iYKVK1u5TosWZ8dpi4jItVH4lVRpyMyN5+2wAsTEJzJk5saLh834eFb3GUr+94aQ+9ghfitemZODx3JXm8YcjT1Gl08G8smmt4kPOkDA3pq0zPcpQ3vX5O/YJczfNp/wIuGEFQy75LqaVMqvsCvXx7p1NnFtwgSoVQu++AJuvhk2boSSJb29OhERn6HwK6nSxQ6RXfD6xESYOJET3V6h4q7tLM9fhk4NXmRJoVsI/PskeQZ1ZcnxsSQGRhG0px4dS77CwKHVzw60yk3YZUOvyHXz+edWj758ufWQbtTI2pOdoeArIpKsFH4lVbqiw2WuC99/D6+8AuvWsTtvcQbc15sZJTJz0u83Yo/9QlzWX/j71HEy7m7Ks5VeoffLVciQIQWfiMi/eTywYIHV6/r7w9q1cOoUDBtm3Rpy5/b2CkVEfJrCr6RKlz1cNncu9OgBS5fCTTfBxIlErMxINL9xMKgbOImQE9hZk9BjHTjwxf0E6KddvGnnTju4Nm4cbN0KP/4I9evbIcwBA9StQUQkhWjkj6RKTSrlZ2Cz8uQPDcHBBkkMbFaeJvGRUKcO1K4NkZEwdiz88Qdrq9/KgejPORg8wIKvA+BHthsLcnOVnAq+4j3790PdulC4MPTqBUWLwmefQc2a9vWgIAVfEZEUpEggqdZ5h8vWrYNXn4EpU+y0+9tvw1NPsXTXDp7s9wRr+BRy+uO/K4LEgj8BCTgEEOpXUa3IJOWtXWu7u40bQ86c1navZ0947DELvyIi4jUKv5K6bd0KvXvbCfgsWaxl2fPP89PW7TzTry0bA76ExCCKRz3DsAe6kJDFodf0yew8uZyCGavSt14zdWeQlBEdbaOGx4yBZcugUCE7vObvDwsXent1IiJymsKvpE7790P//tav198funaFl17i27+207l/G7ZlnAxkouyxLox8tDN3Vs5z9lubVOrkvXVL+vThh/D887bDe/PN9s7Eo4+qnEFEJBVS+JXU5dgxeOstu8TGsrhDQ+bXK0NMXDHGDmzNniw/gl82qp54ldFtn6NS6ZzeXrGkR0eOWN3u3XdD6dLWjuz++6FdOwgLU+gVEUnFFH4ldTh1ynZ5+/eHgwfh/vtZ/Hxz7vqpNfHLptoBtqCs3BnfnzEdnqZkoVBvr1jSG9e1FmVjx8LXX0NsLAwaZOG3Rg2NGhYRSSMUfsW7EhNtB61XL9i+HWrXxjNgIAM3R9Hvu/bEh5w63bnB4eU7X2BQ3Ve8vWJJjzweqFYNVq6ErFnt4NoTT0ClSt5emYiIXCWFX/EO17U+p927WyeHKlVIGDWGVyLjeO+zTpzMsRQ/cuFPIDgegvyDuLdshLdXLemF68LPP8OMGTBwIPj5wUMPwbPPWnlDxozeXqGIiFwjhV9JeYsWQbdudgK+RAniJnzBC4cC+XBaN07lXElAUGFa5Xif9x5/jHWHVjJ/23zCi4RrBLFcf4cOwSefwOjRsHEjhIZCp06QP78duhQRkTRP4VdSzh9/2E7vd9/BjTdy8u0RdIzNxmfLXichxzqCg0vQMf9HvNWqJRmCAgEIyxim0CspY8ECuOceqz8PC7NpbNrlFRHxOQq/cv3t3m29ej/6iMU3ZWDmq3fzC1WZv+0dPNk3kjG4DC+WmEC/hx4k0F8/kpJCjhyBTz+1oSkPPwxVq0KHDtC2LZQv7+3ViYjIdaKkIdfP0aMweLD1PE1IYGaHZtTPNRUPP4HzExlDitPjlkl0b9IcP0eTtiUFuC4sXWqdRb780jo2PPywXTJmtJ9VERHxaQq/kvzi4ixc9OsHBw9yqPEDPJq/HNOzvAVOPAB++PFK/bb0qHG/lxcr6UqHDvDBB5A5M7RpA08+CRUrentVIiKSghR+Jfl4PPDVV9CjB2zZwoGwu2jZsBqzcnwBWSeRPb4sx/1i8LiJBPkHUbNITW+vWHzdhg32QqxbN8ibF5o3t7D7yCM2LltERNIdhV9JHvPmwUsvwfLl7C9VjpYPtWN2ge8h888USLyTIfU/5sGqtVmya4m6N8j1FR8PU6bAyJEwfz4EBUF4ODRtCnXq2EVERNIthV9JmrVrbVdt2jT258tPq6b3MbPEfMg0lmLUZlizL2lU/q6zNw8rqO4Nch2dPGmjhiMjoUgR69Hbti3ccIO3VyYiIqmEwq+PmrIqkiEzN7I7KoZ8oSF0jShFk0r5k+8Bdu9mcb8nmbf+ByoczMh7jeow4+blEPI1pf3r8u79r3J3qduT7/FELsTjgTlz7BBbz552aK1DByttqFsX/P29vUIREUllHNd1U+zBqlat6i5fvjzFHi+9mrIqku6T1xITn3j2upBAfwY2K5/0AHz8OAwZwq+fv0HtB08Re+blkwMVMjRixIOvUr1ItaQ9hsjlHD1qwyhGjIBNmyBPHvs3a1Zvr0xERFIJx3FWuK5b9d/Xq7+UDxoyc+N5wRcgJj6RITM3XvudJiTAmDG4xUtw4M2+tL8zqwVfxy4dqz7D6pe/U/CV62/mTJu49txzkCMHTJgA27cr+IqIyBVR2YMP2h0Vc1XXX5LrwowZJL7YlYM71vPiHfn4okoGPIEH8HP8cRwI8g+i5S0PJ3HVIheRkGBTAUNDoVYtqFwZHngAOna0wRQiIiJXQeHXB+ULDSHyAkE3X2jI1d3R6tWceq4rh1b9xCvVs/FpsyASA/ZS58YWDG/+ClGxUercINfP/v0wdiy8/z7s2mVtymrVgty54aOPvL06ERFJoxR+fVDXiFIXrPntGlHqyu5g1y6Ov/Aqh2d9TL/qwYx7LgCP/3EaFmrJW/e+wk05bzp7U4VeuS569YI33rCBKXXqwHvvQcOG3l6ViIj4AIVfH3TmUNtVd3uIjubgS4M5PGkIb4bF89Gz/rj+iTQv3oZBDbtTLHux6794SZ/i4mxAyr332vS1YsVs+lrHjlC6tLdXJyIiPkTdHgQSE9nW6yPm/9CFEbcdY+WN4PgF0qLM4/SP6Ebh0MLeXqH4qn37bALbqFGwdy98/DG0bu3tVYmIiA+4WLcH7fymY64La976icPDn+atGpuY1vTMV/zpduv7DKz/eLI+3nXvPSxpR2wstG8PX35pu7716ln3Bk1fExGR60ytztIhjwfmjtjAxOLhDP21DnXabmL6Tad/FBwAlw+XLGTKqshke8wzvYcjo2JwgcioGLpPXpusjyGpXHw8/PabfZwhg+36PvkkbNwI06ZBRAT46VeSiIhcX9r5TUcSEmDqhwfZ/OZzrLrlCyY96hLkBpKPxpyKL8vBoEG4bgIOAfgnlGXIzI3JtjN7qd7D2v31cQcPwujRMHKkfbxrF+TKBTNmgON4e3UiIpLOJCn8Oo4TCowFygEu0NZ13cXJsTBJPrGxMH5sHOuGv8LOisOZ+kg8IZ4AOld6iq51XuW2fr/hD+SJe51Yv7Vk8JQn2FPmkn2Br7aEIVl7D0vasH079OtnQyhOnYK777bShuzZ7esKviIi4gVJ3fkdDsxwXfc+x3GCgIzJsCZJJtHRMOp9l0Vj3ia2ci9mtjxB5nh/upV+gs6NB5ArYy7gXF/gYE8Zgj1lzn7/xfoC/3t88pkSBuCiATjZeg9L6ubx2Ojh7Nkt8E6cCG3awLPPws03e3t1IiIi1x5+HcfJCtwJtAFwXTcOiEueZfm2633w69AheOcd+HJmH2KrDGF7y5NkjfOjZ74WvNDyPXKE5Djv9lfbF/haShiS3HtYUreTJ2H8eBg+3FqTffstlCxpHRwyZ/b26kRERM5Kys5vMeAAMM5xnArACuA513VP/PNGjuO0B9oDFCpUKAkP5xuuZdf0Su3aBW+9BV9NnUFIeCv+rncAXAjAj69a/8A9Jetd8Puuti/wtZQwXHPvYUndIiNhxAj44AM4fBiqVLHRw2co+IqISCqTlPAbAFQGOrmuu9RxnOFAN+DVf97Idd3RwGiwPr9JeDyfcD0Ofv31lw3D+nz2L+S74wkiW28kQzw4LrgOuI7Din2rLxp+wcLplT7+tZYwXM1jSCrnulazO3o0DBoETZrACy/AHXeolldERFK1pPQV2gXscl136enPv8bCsFxCch78Wr0aHnwQStVazI+nKhPT9k6OFNjI67vL8H3EODIEhuDv+BPkH0R4kfAkrvycrhGlCAn0P+86lTCkA4mJVs5QowZMmWLXdeoEf/8Nkyfb9Qq+IiKSyl3zzq/runsdx9npOE4p13U3ArWBP5Jvab4pOQ5+LVwIAwbA9HW/ElqrC267xSScgEFr89DxiTFkiWgEwJwipZi/bT7hRcIJKxiWbM9BJQzpzMmTNnlt6FDYvBkKF7YgDNayLFcury5PRETkaiRpvLHjOBWxVmdBwBbgMdd1j1zs9hpv/N+aX7Bd04HNyl8yPLqutUUdMAAW7viFjHe/wslCv5D7BHRdlZEOTfqT+clOEKDWzZLMbr3VhlPcdhu8+CI0baqfMxERSfWuy3hj13VXA/+5U7m4q901TUyEb76BnqMW81emcQRVWA53ryLzCYe+s/14qloHMk3qd653qkhSrV8Po0bB4MEQEgJ9+kC2bHD77SprEBGRNC9JO79XSzu/Vy4uDj791A6y/RU6Auo/C44HXHh2KQwMiCDjkGHWVkokqVwX5s2DN9+E6dMt9M6aZQfYRERE0qDrsvMrye/ECTtA/+ZbLruD5pG5wWuQa4HNzwP8XThV634ydp3k3YWK7zh0COrUgVWr4IYbbCrbU0+plldERHxSUro9SDI6fBj69oVChV06vzuHE/ffBa1rkyl0Oc8thgwJ4HgcPE4QM46FMWVVpLeXLGlZdLTt9ALkyAFlysCYMTaSuGdPBV8REfFZ2vn1sj177BD9+6NcTtwwh+xt+kCWRWT2C6X//My0XXCcKWXrMv3k/9ifcQsZPOXBUzJJfYElHdu718b/jRwJ8fGwe7fV8372mbdXJiIikiIUfr1kyxY7T/TROJeEQrPJ+fRrnAj5lUzBN9D/t4K0/W4nGf53Bw1a3s/6PMUByJZwrmzlWvoCSzq2Ywf07w+ffGKht3lz6NrVgq+IiEg6orKHFLZ2LTzyCBSv+Stj9j9G5q7lcVtGEJJrGyP338rfPffTcamHDOM/hwULiCpV7oL3czV9gSUdi421f6OjYcIEaNsWNm2Cr76yFmYiIiLpjMJvClmyBBo3hltucflq59vwWA08FT4mKvAPuvrfyV/9ougwdjXB3V6BP/+EFi3AcTRNTa6e68K0aRAeDq1b23Vly1qNzfvvQ4kSXl2eiIiINyn8XkeuC7NnQ82aEBbmMm/XDPL3CSO+dmdrWwb4eVyyz1xAcK17YMMGe2s6c+az99GkUn4GNitP/tAQHCB/aMhlB2JIOhUfb/3xbrkFGjSwaWzVq5/7ukocREREVPN7PXg8MGWKTWNbscIlx60zKdS3Dzs8S8mZrTAv52vHO+s+Ig4PQa5D+AtvQ7PnLnp/TSrlV9iVyxs4EHr3hnLlYPx4eOghCAz09qpERERSFYXfZBQfD59/boMpNmxwubHGDIq83odt8cvIkqUwo299h9bf7SDopeHcWywD89uEE/5gN8KK1vD20iUtOnIERoyAsDCoXRvat4cqVaB+fU1iExERuQiF32QQEwMffWTdG3bscClSZwbFW/dhc+wyCmcszJgao2m1KYSgB7pBZCS0aUPYoEGE5cnj7aVLWrRnD7z9ttXvHj8OPXpY+L3xRit3EBERkYtS+E2Co0etXeqQLxZzJNs8ClfPQMlOX7LpxDKKZCjCmDpjaBVQhaBnX4Cff4bKle2UfViYt5cuaVXv3vbWQnw8PPggdOtmNb4iIiJyRRR+r8H+/TBsmL3jfCzrrzhta4JfHNuBG50bGdtoLK0KNyaw3+vw3lN20GjUKGjXDvz9L3v/IudZuxZKlYKgIMiTxzo4vPQSFC/u7ZWJiIikOQq//zJlVSRDZm5kd1QM+UJD6BpR6uxhs+3b4c03YexYiD3lUr31DP4q3Z79sXEA+OHH01U78vjv/lCvHBw4YHWYr78OOXN682lJWrRokR1i+/FH+PBD69HbsaO3VyUiIpKmKfz+w5RVkXSfvJaY+EQAIqNi6D55LTu3BLDi+zx89hm4uNR6Ygb7yvRh0eFl3BhwI4F+gXhcD0FOALUHTYIZ6+B//7Neq1WqePlZSZriujBnjrW8+/lne9HUrx80bertlYmIiPgEhd9/GDJz49ngC3BqTzYOLCnOs6/fQIYMLvWfm8nOYn2YfWAphRMLM6bRGFpVaMWKDXOZ/3FvwictI+zUfhg3Dlq1Aj+1UZZr0KWLvWvw9tvwxBOQKZO3VyQiIuIzFH7/YXdUDK4Lp3bm4OjiEsRuyw3BcWRs+BFl6o/hu31LKRxXmNENR9O6YmuCnAD46CPCunUjLCoKnnkO+vSB0FBvPxVJKzwemDwZ3n0Xpk61n51vvoECBSA42NurExER8TkKv6d5PJBhTwG2/FSIOL8/oeT7ZKh+Ek/RqZwM+JMDsYXOhV7/IFi50uovly6FO+6w0286dS9XKiEBJk60SSgbNkDJkrBtG1SsqINsIiIi11G6D78JCTBpkp0r+nNdBfxLz4UH6oITT6wD/m52OtwyiGGNX7DQGxUFPTtbj9VcueCTT+DRRzVUQK7c0aNWC755M5QvbyH4vvvUCURERCQFpNui1NhY+OAD6yD1yCOQ6HF56YPZZG/zKPjFgwPg8ECZxxnZ9GWC/AJtZGypUhZ8O3aEjRuttlfBVy4nJgZmz7aPs2WDJk1sBvbq1davV8FXREQkRaS7nd/oaAu9b70Fe/dCtVtdWr82h1lxfRi8cxE3ZLzhXPcG/yA63d7M+qx27AgLF1oXh+nTbWCFyOWcPGk/cIMH2yG27dshf37rmSciIiIpzqfD7z979uYOzELRfRWZ/XVWjhyBWrVdXnxvHlOP9qb35oUUyFqAkfVH0rZSW1buWcn8bfMJz12NsKGT7DBSaKg1+H3sMXVxkMs7edLeIRgyBPbtg5o14YsvLPiKiIiI1/hs+D3Tszf6cADHlpVh+5pC/BYfwG01Y2jVdQlf7utD13ULyJclH+/Ve492ldsRHGCn68MK/I+wX7bCA49acGnf3g4m5cjh5WclacaBA9C9O9x5p420rlHD2ysSERERfDj8nunZG7v9RqJXFCHTzbsJvvMbNtw4jqeXrSFv5ry8U/cdnqjyBBkCMpz7xj//tBKHefPsUNJ330G1at57IpI2HD9uHT9+/x0++wwKF7aa8KJFvb0yERER+QefDb+7o2IACCg7hz9543EAABKXSURBVIw39yQuw5+c8N+Mvyc7w+sNp32V9ueH3pgYG0M8eLANFRg50nZ8dRBJLiU6Gt57z4rIDx2CiAj7WQoJUfAVERFJhXw2/OYLDWHLsZXsC+oOTgK4kCW+MTdneopnb6t3/o2nTYNnnoGtW61t2ZAhkCfPFT3OP+uK84WG0DWiFE0qqa4zXVi0CBo3hsOHoV496N0bbrvN26sSERGRS/DZk1tdI0qRGLAe8Jy+xo8MfjnoVvcfgyh27YLmzaFBA5umNW+etTO7iuDbffJaIqNicIHIqBi6T17LlFWRyf10JLWIjoY//rCPy5eHOnVs0Mm0aQq+IiIiaYDPht8mlfLT+a6m+DmB4Prh5wTS+a6mtiubkABDh0Lp0hZaBgyANWsgPPyqHuNMXfE/xcQnMmTmxmR8JpIqnDhhJTFFi8IDD4DrQtasNqDi1lu9vToRERG5Qj5b9gDQrXZj7io5z9qWFQknrGAY/PordOhgB5MaNLA2ZtdYm3mmrvhKr5c0KCYGRo2CQYNg/36oWxdee02DTURERNIonw6/AGEFwyz0HjoETzxhvXoLFIDJk23KVhJCTL7QECIvEHTzhYYkZcmSmnzzDXTuDHffbaH39tu9vSIRERFJAp8tezjLdeHjj63EYdw46NIFNmyApk2TvHvXNaIUIYHnd4MICfSna0SpJN2veNGpUzacYtw4+/yhh2yy3+zZCr4iIiI+wPfD74gRNpWtZElYudI6OWTOnCx33aRSfgY2K0/+0BAcIH9oCAOblVe3h7QoPt7eFShZ0vo8f/edXR8QANWre3dtIiIikmx8vuyBNm3sYFLLltdlLHGTSvkVdtO6GTPg6adhyxY7vDZ6NNxzj7dXJSIiIteB74ffzJmhVStvr0JSG4/HShxCQiAoCEJD4YcfoH59HWYTERHxYb5f9iDyT64L06fbyOoePey6WrVg+XLr/qHgKyIi4tMUfiX9WLAA7rzTdnePHDm/P69Cr4iISLqg8Cvpw8CBcNddVtf7/vvw55/QooW3VyUiIiIpzPdrfiX9Wr/exlaXKAHNmlltb8eOVucrIiIi6ZJ2fsX3bN4Mjz4K5ctDr152XalS8OKLCr4iIiLpnMKv+I69e61lWenSNpmta1cbXy0iIiJymsoexHcMGwYffGBjrF99FfLl8/aKREREJJXRzq+kXbGx8PbbMHeuff7yyza6+v33FXxFRETkghR+Je1JTISPP7Y63s6dYepUuz57drjpJq8uTURERFI3hV9JW2bNggoV4LHHIHdumD0bhg/39qpEREQkjVDNr6QNrmuDKP76C+LiYNIkaN4c/PT6TURERK6ckoOkbuvXQ+PGMHasfd6+vV13//0KviIiInLVlB4kddq927o23HKLjSVOTLTrAwPtIiIiInINVPYgqc8HH8ALL0BCAjz7LPTsCTlzentVIiIi4gMUfiV1iI+3sBsSAgULWqnDgAFQrJi3VyYiIiI+JMllD47j+DuOs8pxnB+SY0GSzrguTJkC5cpB//52Xf36MHGigq+IiIgku+So+X0O2JAM9yPpzZIlcOed0LSpHV67/XZvr0hERER8XJLCr+M4BYAGwNjkWY6kG4MHQ1iYtS4bNQrWroUGDby9KhEREfFxSa35HQa8BGS52A0cx2kPtAcoVKhQEh9O0rTDh61H7403Qr16cPIkvPgiZLnoj4+IiIhIsrrmnV/HcRoC+13XXXGp27muO9p13aqu61bNnTv3tT6cpGVxcTBsGJQoYeOIAcqXhz59FHxFREQkRSWl7KE60NhxnG3ARKCW4zgTkmVV4hvOHGYrW9Zal1WrBj16eHtVIiIiko5dc/h1Xbe767oFXNctAjwEzHVdt2WyrUzSvqFD7TBbUBBMnw4zZ1pXBxEREREvUZ9fSV67d8OxY1C6NDz6KGTMaJPaAvSjJiIiIt6XLOONXded77puw+S4L0mjTpyAvn3hppvgqafsuhtugA4dFHxFREQk1UiW8CvpmMcD48dDqVLQu7cNqPjwQ2+vSkREROSCtCUnSTNuHLRrB1Wr2lS2O+7w9opERERELkrhV67ezp12uf12eOQRyJwZ7r/fprSJiIiIpGJKK3LlTp60ut5SpaBtWyt5yJABHnxQwVdERETSBCUWuTzXhUmToEwZq+tt2NDalinwioiISBqj9CKXN3u27e5mzw7z51sQLlzY26sSERERuWoKv3JhBw7Y7i5AnTrw9dewYgXcdZd31yUiIiKSBAq/cr74eBg2zPr1tmhh/XsdB5o3B39/b69OREREJEkUfuWcmTPhllvghRfgtttg0SLIlMnbqxIRERFJNmp1JubPP6FuXShRAr7/Hho0sB1fERERER+ind/0LDYWpk2zj0uXttC7bp11c1DwFRERER+k8Jte/fADlC0LjRrB1q12XcOGEBzs3XWJiIiIXEcKv+nN5s0Wchs1sqA7axYULertVYmIiIikCNX8pifHj0PVqvy/vXuPlaK+Ajj+PRE0UvDFw1IftZhasSZVJAaqGEN9IBix9REJsRhNjFGjREzEGPEVjVjR+GjaaPBV0RoftNqosanGWqM8VFAQ6oNApFKk1YgFAS/8+scM8brsXu5yLzt7d76fZHJnZ37LnpycHU5mfzNDWxvcfjtcdhn07l10VJIkSQ1j89vqUoLXXoNRo6BvX3jwwexODoMHFx2ZJElSwzntoZUtXQonn5w9mGLrAytOP93GV5IklZbNbytavx6uvjq7Z+/cuXD33XDCCUVHJUmSVDinPbSalGD0aJgzByZNgunTYd99i45KkiSpKdj8topPPsmmM/TqBddeC/36wXHHFR2VJElSU3HaQ0+3aVN2dvfQQ+Gee7Jt48bZ+EqSJFXhmd+e7NVX4eKL4f33swvZzjij6IgkSZKammd+e6obb4Tjj88ubnvuOZg9Gw48sOioJEmSmppnfnuSLVtg40bYfffsoraNG+Gaa6BPn6IjkyRJ6hE889tTvPMOjBwJV12VvT72WLj5ZhtfSZKkOtj8Nrt16+DKK7PHEi9fDiNGFB2RJElSj+W0h2b2xhswYQKsWAEXXgi33gp77110VJIkST2WzW8zGzQI+veHRx/NpjlIkiSpS2x+m8mWLXD//fD66/Dww3DwwTB/PkQUHZkkSVJLcM5vs1i8OHswxUUXZU9rW7cu227jK0mS1G1sfou2YUP2OOIjj4QlS+DBB+Hll6Fv36IjkyRJajlOeyja119nUx3OOQdmzICBA4uOSJIkqWV55rcIn38O110HbW3Z3RsWLYJHHrHxlSRJ2slsfhspJXjqKRg6FG65JbuVGcCAAcXGJUmSVBI2v42yahWccQacdRbsvz/MmwejRhUdlSRJUqk457dRzjoL3noLpk+HK66AXqZekiSp0ezAdqZly7J5vP36wb33Qp8+cMghRUclSZJUWk572Bk2b4Y774TDD4dp07JtRxxh4ytJklQwz/x2t0WL4IILYO5cGDcOpkwpOiJJkiTlPPPbnWbNgmHDsukOjz0Gzz2XXdwmSZKkpmDz2x22bMn+jhwJEyfC++/DhAk+mliSJKnJ2Px2xYYNcPXVcOaZ2T18hwzJHk/swyokSZKaks3vjpo7F446Cm69NXtK2zffFB2RJEmStsPmt15bz/aOHAlr18ILL8DMmbDrrkVHJkmSpO2w+a3X+vXw0ENw3nnZnR3GjCk6IkmSJHWSzW9nbNgAd90FbW2wzz5Z0ztzJuy5Z9GRSZIkqQ42v9szb142t3fyZHjxxWxb//7FxiRJkqQdYvNby8aN2dzeESO+ndt76qlFRyVJkqQu8AlvtUycCE8/DeefD3fc4RQHSZKkFrDDZ34j4oCIeCUilkTE4oi4vDsDK8SmTdkFbQBTp357JwcbX0mSpJbQlWkPbcCUlNJQYARwSUQc1j1hFWDx4myKw5Qp2evhw72TgyRJUovZ4eY3pbQqpfR2vv4VsATYr7sCa5jNm2HGjOyitpUr4aSTio5IkiRJO0m3zPmNiIOAI4E53fHvNcyKFXDuufDaazB+PNx3HwwaVHRUkiRJ2km6fLeHiOgLPA1MTimtrbL/woiYHxHz16xZ09WP616bN8PHH2cPrZg928ZXkiSpxXWp+Y2I3mSN76yU0jPVxqSU7kspDU8pDR84cGBXPq57rFoFt9wCKcGQIbBsGUyaBBFFRyZJkqSdrCt3ewhgJrAkpXRH94W0Ez35JBx+ONx0Eyxdmm3bbbdiY5IkSVLDdOXM7zHAucDoiFiQL2O7Ka7u9cUX2X17zz47O9v7zjswdGjRUUmSJKnBdviCt5TSP4DmnyuQEpx4IixcCDfckD21rXfvoqOSJElSAVr/CW8RcNttsMce2b17JUmSVFqt3/wCjB5ddASSJElqAl2+1ZkkSZLUU9j8SpIkqTRsfiVJklQaNr+SJEkqDZtfSZIklYbNryRJkkrD5leSJEmlYfMrSZKk0rD5lSRJUmnY/EqSJKk0bH4lSZJUGja/kiRJKg2bX0mSJJWGza8kSZJKw+ZXkiRJpWHzK0mSpNKw+ZUkSVJp2PxKkiSpNCKl1LgPi1gDrGjYB35rAPCfAj63pzJf9TNn9TFf9TFf9TFf9TFf9TFf9SkyXz9MKQ2s3NjQ5rcoETE/pTS86Dh6CvNVP3NWH/NVH/NVH/NVH/NVH/NVn2bMl9MeJEmSVBo2v5IkSSqNsjS/9xUdQA9jvupnzupjvupjvupjvupjvupjvurTdPkqxZxfSZIkCcpz5leSJElqreY3IsZExD8j4qOImFplf0TE3fn+dyNiWBFxNoOIOCAiXomIJRGxOCIurzLm+Ij4MiIW5Mu0ImJtFhGxPCLey3Mxv8p+6ysXET9pVzcLImJtREyuGFP6+oqIByLis4hY1G7bPhHx14j4MP+7d433dni8a0U18vWbiFiaf+dmR8ReNd7b4fe3FdXI1/UR8a9237uxNd5rfWXbnmiXq+URsaDGe8tYX1X7iB5xDEsptcQC7AJ8DAwBdgUWAodVjBkLvAAEMAKYU3TcBeZrMDAsX+8HfFAlX8cDfyk61mZZgOXAgA72W1/V87IL8G+y+y223176+gKOA4YBi9ptuw2Ymq9PBabXyGmHx7tWXGrk6ySgV74+vVq+8n0dfn9bcamRr+uBK7fzPuur+v4ZwLQa+8pYX1X7iJ5wDGulM79HAx+llJallDYBfwTGV4wZDzySMm8Ce0XE4EYH2gxSSqtSSm/n618BS4D9io2qx7O+qvsF8HFKqYgH3DS1lNLfgc8rNo8HHs7XHwZOr/LWzhzvWk61fKWUXkopteUv3wT2b3hgTapGfXWG9VUhIgI4G3i8oUE1sQ76iKY/hrVS87sf8Em71yvZtpnrzJjSiYiDgCOBOVV2j4yIhRHxQkT8tKGBNZ8EvBQRb0XEhVX2W1/VnUPt/zCsr23tm1JaBdl/LsCgKmOsterOJ/v1pZrtfX/L5NJ8msgDNX6Str62NQpYnVL6sMb+UtdXRR/R9MewVmp+o8q2yltZdGZMqUREX+BpYHJKaW3F7rfJfqr+GXAP8KdGx9dkjkkpDQNOAS6JiOMq9ltfFSJiV+A04Mkqu62vHWetVYiIa4A2YFaNIdv7/pbF74CDgSOAVWQ/5VeyvrY1gY7P+pa2vrbTR9R8W5VtDauxVmp+VwIHtHu9P/DpDowpjYjoTVaws1JKz1TuTymtTSn9L19/HugdEQMaHGbTSCl9mv/9DJhN9rNNe9bXtk4B3k4pra7cYX3VtHrrdJn872dVxlhr7UTEJOBUYGLKJxRW6sT3txRSSqtTSptTSluA+6meB+urnYjoBfwKeKLWmLLWV40+oumPYa3U/M4DfhwRP8rPNp0DPFsx5lng1/lV+SOAL7eemi+bfP7STGBJSumOGmO+n48jIo4mq5f/Ni7K5hER34uIflvXyS6yWVQxzPraVs2zJdZXTc8Ck/L1ScCfq4zpzPGuFCJiDHAVcFpKaX2NMZ35/pZCxXUIv6R6Hqyv7zoBWJpSWlltZ1nrq4M+ovmPYY26sq4RC9nV9h+QXUF4Tb7tIuCifD2A3+b73wOGFx1zgbk6luwnhneBBfkytiJflwKLya7CfBP4edFxF5ivIXkeFuY5sb62n7M+ZM3snu22WV/fzdHjZD89f0N2JuQCoD/wN+DD/O8++dgfAM+3e+82x7tWX2rk6yOyuYNbj2O/r8xXre9vqy818vWH/Pj0LlmzMdj6qp2vfPtDW49b7cZaX7X7iKY/hvmEN0mSJJVGK017kCRJkjpk8ytJkqTSsPmVJElSadj8SpIkqTRsfiVJklQaNr+SJEkqDZtfSZIklYbNryRJkkrj/3aICN9S7LHhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(x1, y2, 'o',label=\"data\")\n",
    "ax.plot(x1, y_true2, 'b-', label=\"True\")\n",
    "prstd, iv_l, iv_u = wls_prediction_std(res)\n",
    "ax.plot(x1, res.fittedvalues, 'r-', label=\"OLS\")\n",
    "ax.plot(x1, iv_u, 'r--')\n",
    "ax.plot(x1, iv_l, 'r--')\n",
    "ax.plot(x1, resrlm.fittedvalues, 'g.-', label=\"RLM\")\n",
    "ax.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: linear function with linear truth\n",
    "\n",
    "Fit a new OLS model using only the linear term and the constant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.55290631 0.39535154]\n",
      "[0.38369939 0.03306109]\n"
     ]
    }
   ],
   "source": [
    "X2 = X[:,[0,1]]\n",
    "res2 = sm.OLS(y2, X2).fit()\n",
    "print(res2.params)\n",
    "print(res2.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate RLM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.07093614 0.48553947]\n",
      "[0.10553732 0.00909352]\n"
     ]
    }
   ],
   "source": [
    "resrlm2 = sm.RLM(y2, X2).fit()\n",
    "print(resrlm2.params)\n",
    "print(resrlm2.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare OLS estimates to the robust estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFlCAYAAAAzqTv+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3yTVRfA8d+TtKUtewqUUWZlDwEZIlsQWSIqKIIMBQSEV0VAkCF7uJiCCooiQ2QIiqggyCjIhiIgQ0bLHgVKd/K8fxzaAha6k7Q9Xz9+oEna3NT4nNx7zz3HME0TpZRSSjmWxdkDUEoppTIjDcBKKaWUE2gAVkoppZxAA7BSSinlBBqAlVJKKSfQAKyUUko5gZsjnyxfvnymr6+vI59SKaWUcprdu3dfMU0zf3z3OTQA+/r6smvXLkc+pVJKKeU0hmGcftB9ugStlFJKOYEGYKWUUsoJNAArpZRSTuDQPeD4REVFERgYSHh4uLOHkqY8PT0pUqQI7u7uzh6KUkopF+D0ABwYGEj27Nnx9fXFMAxnDydNmKbJ1atXCQwMpESJEs4ejlJKKRfg9CXo8PBw8ubNm2GDL4BhGOTNmzfDz/KVUkolntMDMJChg2+MzPAalVJKJZ5LBGBXMmrUKKZOnfrA+1euXMnff//twBEppZTKiNJdAF65N4h6EzdQYshP1Ju4gZV7gxz7/BqAlVJKpYJ0FYBX7g1i6PKDBAWHYQJBwWEMXX4wxUF43Lhx+Pn50bRpU44ePQrA559/Ts2aNalSpQrPPfccoaGhbNu2jR9//JFBgwZRtWpVTpw4Ee/jlFJKqYSkqwA8Zd1RwqJs99wWFmVjyrqjyf6Zu3fvZvHixezdu5fly5ezc+dOANq3b8/OnTvZv38/5cqV48svv6Ru3bq0adOGKVOmsG/fPkqVKhXv45RSSqmEOP0YUlKcCw5L0u2JsXnzZp599lm8vb0BaNOmDQABAQEMHz6c4OBgQkJCaN68ebzfn9jHKaWUcnHXrkGePA57unQ1Ay6cyytJtydWfBnKr776KjNmzODgwYOMHDnygUeIEvs4pZRSLsg0YcsW6NABfHzg/HmHPXW6CsCDmvvh5W695zYvdyuDmvsl+2c++eSTrFixgrCwMG7dusXq1asBuHXrFoUKFSIqKoqFCxfGPj579uzcunUr9usHPU4ppVQ6cOAA1K8PGzbAwIHgwGqF6WoJul01H0D2gs8Fh1E4lxeDmvvF3p4c1atX58UXX6Rq1aoUL16c+vXrAzBmzBgef/xxihcvTqVKlWKDbseOHXnttdeYNm0ay5Yte+DjlFJKuaArV2DuXLh9G8aNgypVYOlSaNkSsmZ16FAM0zQd9mQ1atQw7+8HfPjwYcqVK+ewMThTZnqtSinlUv7+Gz79FBYsgPBwaNMGVq6ENC6SZBjGbtM0a8R3X4JL0IZhzDMM45JhGAH33d7fMIyjhmEcMgxjcmoNVimllEpVM2dChQoSfDt3hoAAWLUqzYNvQhKzBP0VMANYEHODYRiNgLZAZdM0IwzDKJA2w1NKKaWSKDRUgm316lCrFjz1FIwdC6+/DvnzO3t0sRKcAZum+Sdw7b6b+wATTdOMuPOYS2kwNqWUUirxAgNh6FAoUgT69JG9XYAyZWDYMJcKvpD8LOiyQH3DMHYYhrHJMIyaqTkopZRSKkkGDIASJWDyZGjcGDZvhilTnD2qh0puAHYDcgO1gUHAUuMB7X4Mw3jdMIxdhmHsunz5cjKfTimllLpLdDSsWAG2O9URixaF/v3h+HFYtgyeeMLpe7wJSe4xpEBguSkp1H8ZhmEH8gH/ibCmac4F5oJkQSd3oEoppRTXr8Pnn8OMGXD2LKxeDa1awTvvOHtkSZbcALwSaAxsNAyjLOABXEm1UTnQ1atXadKkCQAXLlzAarWS/84+wV9//YWHh4czh6eUUgrg1i0YPBi+/lqSrBo3luzmp5929siSLcEAbBjGIqAhkM8wjEBgJDAPmHfnaFIk0NV05IHiVJQ3b1727dsHSC/gbNmy8c59n6RM08Q0TSyWdFU4TCml0jfThH//hZIlwdsbNm2CF1+U/d4qVZw9uhRLTBZ0J9M0C5mm6W6aZhHTNL80TTPSNM3OpmlWNE2zummaGxwxWEc6fvw4FStWpHfv3lSvXp2zZ8+SK1eu2PsXL15Mz549Abh48SLt27enRo0a1KpVi+3btztr2Eoplf6Fhkq1qooV5RhRWBhYrbB/P8yblyGCL7hYKcqBA+HOZDTVVK0Kn3ySvO/9+++/mT9/Pp999hnR0dEPfNybb77Ju+++S+3atTl16hStWrUiICDggY9XSikVj/PnYfp0mDNHOhNVqwYffSTBF8DNpUJWimWsV5PKSpUqRc2aCZ+w+v333zl6NK4n8fXr1wkLC8PLK2VdmpRSKlOIjAQPDzh5EiZNgnbtZEaWDjKZU8KlAnByZ6ppJetdhbktFgt3b3Pf3XbQNE1N2FJKqaSIioLly+XCX60azJoFdevCqVNypCgT0KyiRLJYLOTOnZtjx45ht9tZsWJF7H1NmzZl5syZsV/vS+11dKWUyiiuXZNZbsmS0LEjXL4sARhktptJgi9oAE6SSZMm0aJFC5o0aUKRIkVib585cyZbt26lcuXKlC9fns8//9yJo1RKKRc2fDgMGQJ+fnKG959/4LXXnD0qp9B2hA6UmV6rUkpht8O6ddIGcORIqFNHlphv3YJKlZw9Ood4WDtCl9oDVkoplQGEhEg3omnT4OhRKFQILl6U+3x9nTo0V6IBWCmlVOqx22VP9/hxqFkTFi6EDh0ky1ndQwOwUkqp5DNN2LoVliyRpWaLBcaMgeLFoXbtDH2MKKU0CUsppVTSRUTAN99AjRpQvz58+y0cOyb3dewo+70afB9KA7BSSqmkOXJE9nK7dJEykbNnQ2CgZDarRNMlaKWUUgnbtw/OnIE2baB0aXjqKXj5ZWjWTGe6yaQzYCAwMJC2bdtSpkwZSpUqxYABA4iMjGTjxo20atXqP49fs2YN1apVo0qVKpQvX545c+Y4YdRKKZXGbDZpet+woSRWvfWW7Pm6uUlbwKee0uCbApk+AJumSfv27WnXrh3Hjh3jn3/+ISQkhGHDhsX7+KioKF5//XVWr17N/v372bt3Lw0bNnTsoJVSKq2tXi0z3fbt5ezulCmwc6cG3FSU6ZegN2zYgKenJ926dQPAarXy8ccfU6JECRo1avSfx9+6dYvo6Gjy5s0LQJYsWfDTfQ+lVEbwzz+QNSv4+EC2bFIW8sMPZdk5g3UicgWu9Rt1Qj/CQ4cO8dhjj91zW44cOShWrBjHjx//z+Pz5MlDmzZtKF68OE2aNKFVq1Z06tQJiyXTLyYopdIj04TffpPr5Nq1ch3++GNo1Ej+zUSOH4d//5VtbUfI9FHDNE2MeJZUHnQ7wBdffMH69eupVasWU6dOpXv37mk9TKWUSn3z5kGFCtC8OezZA6NGSZ3mTOb8eWjxfAhl/ey0fOE2dcdvYOXeoDR/XteaATuhH2GFChX44Ycf7rnt5s2bnD17llKlSj3w+ypVqkSlSpV45ZVXKFGiBF999VUaj1QppVLB+fNSGhJg82bw8pKEqhdfhCxZnDs2B7txAyZPhg8/shMR6U32qmfIWfcY525GMnT5QQDaVfNJs+fP9DPgJk2aEBoayoIFCwCw2Wy8/fbbvPrqq3h7e//n8SEhIWzcuDH263379lG8eHFHDVcppZLONCXYdugARYrA7t1y+6xZsGuXnOfNRME3PFy2tkuWhPHjIXvZyxTuuYk8zQ5hzRoJQFiUjSnrjqbpODJ9ADYMgxUrVvD9999TpkwZypYti6enJ+PHjwdg/fr1FClSJPbfvXv3MnnyZPz8/KhatSojR47U2a9SyjVFRkpThBo14MknYcMGeOcdKFxY7vfyylRZzTYbzJ8PZcvKr6FmTVl5z/b0LtxzhxJhOcwNt6VEWA4DcC44LE3H41pL0E5StGhRVq9e/Z/bGzZsSFjYf/8D1K9f3xHDUkqp5ImOlqzliAjo31+ymj/7DDp3liznTMY04ccf4b334O+/JfB+/XVcjlnhdV6cuLmLix7DgGgMPHgkchwlc1RP03Fl+hmwUkplGLt2wSuvQK1aEnWyZ5fl5kOHoFevTBl8N2+GJ56Adu3kc8n338OOHXHB127aqV3hMJc9xoERDQaYRGNzO8Sg5ml7xFQDsFJKpWfR0bB0KdSrJ1O7lSulOULM6l3p0plqmTnGwYPQqpWsvJ86BXPnyueQDh3k12GaJquOrKLqZ1X5aE8/CmTNjYEbmBYshjtvNXg2TROwQJeglVIqfVuxQjKYS5aUkyTdukGOHM4eldOcOgUjRkhzppw5YeJEWYX39gb/s/78ceoPvN29+e7gd+w8t5Myecqw6LlFvFDhBXYE7mDjqY009G1InaJ10nysGoCVUio9CQiAadOgfHkpmtG2rZSNfPppsFqdPTqnWLk3iPE/nOTIL0UI2VccdzeDQYMMBg+GPHnkMf5n/Wn0dSMibBEAPJL1Eb5o/QVdq3bFzSKhsE7ROg4JvDE0ACullKuz2eCnn6Th/YYNkr38zjtyn4eHrLVmUou2nOPNYaFc3V4HM8pKtkpneaThCep0LEuePLKEvOf8Hl5d+Wps8LVgoW/NvvSo3sOZQ9cArJRSLu+11+T8TNGisqbasyfcqUefWUVGwpw58NbQfETfLox32fPkevIo7nlvEw1MWXeUskVuMOKPEfxw+Aeye2TH3eKO3bTjYfWgacmmzn4JGoBBGjBUqlSJ6OhoSpQowTfffEOuXLk4deoUrVq1IiAg4J7Hv/rqqyxdupSLFy+SPXt2AAYMGMC0adO4fPky+fLlc8bLUEplFEePwowZMsstXhxefx1atpRU3kzYFGHl3iCmrDvKueAwCuXw4nGzKqu+zMO//4JnsVvka3+ELIWDYx8fZVxg/+2FVJy1kaweWRnZYCT/q/0//r78t0P3eBOS+f5LxsPLy4t9d5pAdO3alZkzZz6wHWGM0qVLs2rVKjp37ozdbuePP/7AxydtM+aUUhmY3Q6//irLzL/8IkvLTzwhAbh2bWePzmlW7g1i6PKDhEbaCP83P7s2+eF/KSe+ZSNZu9aDD/bu59yNMCIshwm17CDKOEuYdRcWw8rbdd5m8BODyectkyJH7/EmJF0eQ/I/68+EzRPwP+uf6j+7Tp06BAUlXIS7U6dOLFmyBICNGzdSr1493DLhJ1OlVCqIipLObU8/LR3hRo+GM2ckuzmTm7LuKMGns3NxUW0ufV8LM8KdfK33UqjrFlq0gHdb+BHtvpsLHoO56baMMOsOspm1+fypbUx5akps8HVFLhUxBv4ykH0XHt6O8EbEDQ5cPIDdtGMxLFR+pDI5s+R84OOrFqzKJy0S1+TBZrOxfv16evRIeGO+TJkyrFq1iuvXr7No0SI6d+7M2rVrE/U8SinFyZPS/q9vX3B3lwOqgwfD88/L7Fdx5Ajsm1+e0H8KYvGOIHfTALJXPYNhNTl/E4LDg9l1fTYX3SeDab/zXRbala9P9zqPPfRnu4J0NwO+EX4D+51ftN20cyP8Rop/ZlhYGFWrViVv3rxcu3aNZolsBtm+fXsWL17Mjh07tDylUiphpilZzG3bSoGMgQMhMFDuGzECXn5Zgy/yK+nZUzolhp/KR84njuLT6w9yPHYaw2piJwx7tuWU+LQE4zaP40nf+ni6eWI1rHi5ZeGNOm2c/RISxaVmwImZqfqf9afJgiZE2iLxsHqwsP3CFK/px+wB37hxg1atWjFz5kzefPPNBL+vY8eOVK9ena5du2KxpLvPMkopRzp4EDp1knJM+fPDsGHQu7fUac7g7k6iKpzLi0HN/eKtMnXtmiR5T58uJ6/694dqbS4z8Y9/CbYFEGbZh51bhLr9ic0WTKtSrRjTaAxVC1bF/6y/SyVYJYZLBeDEqFO0Duu7rE+TX3TOnDmZNm0abdu2pU+fPgk+vlixYowbN46mTZ2fzq6UckFnzsDFi1IismhRyJ1bjhN17Aiens4enUPEJFGFRdkACAoO+0+v3dBQCboTJ0qP3s6d4YMPwNcXoBCBdn/e3zoU05RazaVyVubbDmupXSQuOc3VEqwSI90FYEjbX3S1atWoUqUKixcvpn79+hw9epQiRYrE3v/xxx/f8/hevXqlyTiUUumUacKWLZLNvGIFVKsmTRJy5ZLOAJnMlHVHY4NvjJheu60q+TB/PowaBefOyUmrCROgcmV5nM1uY3HAYj7a/SYmEnwthoUej3W8J/imVwkGYMMw5gGtgEumaVa87753gClAftM0r6TNENNeSEjIPV/f3ZowKirqP49//vnn4/05p06dStVxKaXSmTVr4P33JZM5Tx4YNAjeeMPZo3Kq+HrqmiYc25GTChXgn3+gTh1YtEgaJ8j9JiuPrOT9P97n0OVDlM5dmpCoEGx2Gx5WDxr6NnTsi0gjiZkBfwXMABbcfaNhGEWBZsCZ1B+WUkqlE0FB0vYvRw64ckW6E33+Obz0knQAyOQK5/Ii6K4gHH46L9c3+RF5Pjf5y0vzpjZtYHugP+M3/4GXmxcLDy5k9/nd+OX1Y0mHJXQo38HhjRIcIcEAbJrmn4Zh+MZz18fAu8CqVB6TUkq5NtMEf39pivDDDzBpErz1lvTi7do1Q7X/S2wC1YMMau7H0OUHuRGYleub/Aj/twBuOcLoN/I6n7yfG6v1v40SCmYryPy28+lcubPTGiU4QrL2gA3DaAMEmaa530jgjWYYxuvA6yBJS0oplW6ZJnzzjQTe3btlX3fAAGjfXu7PYN2IEpNAlZBKOXzIvys3R37xxuIZSfGWx5j8vjcv1Jbv33VuF11WdvlPo4RXq76a+i/IxST57IxhGN7AMGBEYh5vmuZc0zRrmKZZI3/+/A96TFKHke5khteoVIZ186b8aRjw9deStjt7thxYnTo1Jl03w3lYAlVCLlyAfv3g0Udh1yZv3nsPrp734NRPZXihtg8BlwJov6Q9NT+vyaWQS7hb3LEaVrK4ZaFJiSZp9ZJcSnJmwKWAEkDM7LcIsMcwjFqmaV5I6g/z9PTk6tWr5M2bl4Rm0+mVaZpcvXoVz0xy7ECpDME0YccOme2uWiXZQj4+sHSpJFhl0OvV3eJLoHrY7SCfVaZMgY8/hvBwaeQ0YgQUKiT3H792nFEbR/Hdwe/I5pGN0Q1HM7D2QA5dOpTh9ngTkuQAbJrmQaBAzNeGYZwCaiQ3C7pIkSIEBgZy+fLl5Hx7uuHp6XnPcSallIuKiIDvv5fAu3OnJFf16hW3vJyJ2gDen0B19+33i4iAWbNg3Di4ehVeeAHGjoUyZeT+lYdXMn7LeHaf200WtywMqjuId+u9S15v+X1mxD3ehCTmGNIioCGQzzCMQGCkaZpfptYA3N3dKVGiRGr9OKWUSh67HSwWKZzRtSuULQszZ0pi1Z22o5lNTALV3cvQXu5WBjX3i/3aZoNvv5VZ7pkz0LSpnOWtUUPuvxhykTfXvsnSv5cC4Ga4saTDElr7tXboa3FFicmC7pTA/b6pNhqllHK0mGXmkBBZai5WTBKsKleWgJyJxSRaxZcFbZpy7Pm99yAgAKpXhy+/lAAMcD3sOlO2TeHTHZ8SFhWGgYF555+ASwEagEmnlbCUUipFIiPjlpn/+ktmuD16xM2Cq1Z19ghTTUqPEbWr5vOfx2/dKo2btm6VnhJLlkgzJ4sFbkXc4pPtn/Ch/4fcjLhJx4odaevXlm6rusXW8M8ohTRSSgOwUirzmTlTzu2WLStFiLt2zZDLzKlxjOhuAQEy4129GgoWlETwHj1g1wV/xm7+jcu3L7P40GKuhF6hrV9bxjQaQ6VHKgFQLGexTJdklRANwEqpjG/HDgm0bdpIdlCXLlC+PDRrlqGXmR92jCgpAfjMGRg5Uk5gZc8uiVYDBkDWrPDn6T9puqApUXYp21uzcE1+euknavnUuudnZMYkq4Rk3HeeUipzi4yEhQvh8cehdm348UdJsALJZG7ePEMHX0jeMaK7XbkiCwVlykit5rffhpMnZRbs6WXj631f025xu9jgazEsPPvos/8Jvip+OgNWSmVMTz8NGzaAn1+GXmZ+mKQcI7rb7dtyjnfKFMlNe/VV6VhUtCjYTTvfH/qBERtHcOTKEcrmLcvtqNsZp1FCVBS4uzvkqTQAK6Uyhh07YM4caQOYPbt0Inr33Qy/zPwwiTlGdLeoKOkj8cEHsljQrp0sN9/I7s83p/7A46wH3x38jr0X9lIuXzmWPb+M9uXasz1we/re37XbYd06mDFDvv7pJ4c8rQZgpVT6FV/RjO7d4YknoEULZ4/O6R52jOhudrsU+Bo+HE6ckLaAK1ZIm8D7GyUUzl6YBe0W8FKll7BapDhJut3fDQ6Gr76SpLzjxyWzrHdvqYLmgEpnGoCVUunT5ctQqZJM1R59VGYvXbpkumXmhMR3jCiGacJvv8HQobBnj/w6f/pJVu8NA3YE7uCVFa/c0yihT40+vFLlFUe+hNQXECDvl2++kbredevKtP+558DDw2HD0ACslEofTBO2b4f9+2WWkj+/bE42bizVHzLpMnNy7dwJQ4bINnnx4rBggbQwtlrhwMUDvP/H+/x49EdyZsmJu8Udu2nHw+qRfhslREdLIt706bBxI3h6ygvu21eqiDiBBmCllGuLiJBKDzEtAAsUkMDr6QkTJzp7dOnO0aOy1LxsGeTLB598Ip9nsmSBf67+w8iNI1kSsIQcWXIwptEYBjw+gIBLAel3j/fSJfjii7juVcWLS//mHj2cXtdbA7BSynX98otkL1+6BOXKSbX/V16R4KuSJCgIRo+GefPAy0vO9b79tqzYLz+8nPGbx7P3/F483T0Z8sQQBtUdRG6v3EA63ePduVOWmRcvliNpTZvKXu8zz7hM32YNwEop12GaUt8wWzYpB1mmjJzh7d8fmjTJFC0AU9v16zB5siSHR0fDG2/IDLhAAbgQcoFuS/vxw+EfAHCzuLG0w1KeKfuMk0edTDFJedOnS4nRbNng9dflRZcr5+zR/YdumiilnC8sTKZm1atD/fqyRAhQqpQ0SGjaVINvEoWFSeAtVUp+ne3bw5EjspJvzXaVwb8NpuSnJVl+eDkG8rs1TZMDFw84eeTJEBgonyqKFZMVkhs3JAgHBcmfLhh8QQOwUsrZPvxQKjz06CFTtDlzZM9OJUt0tHQlKlNGGibUrg1790rLwHw+Nxm9cTQlp5VkyrYptC/XniUdluDp5onVsKavQhqmCZs2wfPPg68vjB8vVc9+/RUOH4Z+/eRYmgvTJWillGOZJmzeLEc/3Nxk2bBBA1lmbtBAZ7rJZJqwcqWUiTxyRGLRwoXyK/3j3z94+tvJbDu7jZuRN3n20Wf5oNEHVCxQEYAiOYqknySr27flhc2YAQcPQp48spndp48E4nREA7BSyjHuv3AuWybnLocO1aCbQps2yZGi7dvlSPTy5VLFKtIWwdvr3uOj7R8BUqv5yzZf0r1a93u+P10kWZ04IUl48+ZJAY2qVWWq36mTZJWlQxqAlVJpKzQURoyQi+XdF86WLeV+Db7Jtn+/fH5ZuxZ8fGTlvmtXwBLNV/u+YfSm0Zy+cTr28QYGF0MuOm/ASXV3ici1ayV7+fnnZXm5Tp10/97RPWClVOozTZmxgBwZWrsWnnpKlp737JFykel01uIKTp6Ezp2hWjWZ9U6eDMeOQbfudn44soSKsyrS/cfu5M+an4+e+ggvN6/0tccbHCwHlP385IPanj1yburMGfjuO9m+SOfBF3QGrJRKTbduSUmlGTOkRGRgIHh7SxaQA0v8ZVQXL8LYsZKn5uYmSVaNuviz68ofzNjjzsKDC9l/cT8V8ldg+QvLafdoOwzDoHaR2uljjzcgQM7qfvONbFnUqwdjxkgKdwZ8/2gAVkql3OnTks381VcShGvWlLVRtzuXmAx48XSkmzfl1/vhhxAeLgnjI0fCadu9jRJ8svvw7bPf0rFix9hGCeDie7wPKhHZr59M8TMwDcBKqeSx26VZbI4ccP48fPYZvPCCZDM//rizR5chRETIr3XsWLhyBTp0kL/7+UmXos7LO9/TKKF3jd68XPllJ486kVy4RKSjaABWSiVNcDDMny9LhY0bw9y5EnCDgqRBgkoxm022OkeMgFOn5Nc8caIsLOy7sI9W3w3np2M/kdszd/prlHB/ichmzVyuRKSjaABWSiXOoUOyTHh3C7eYnruGkeGC78q9QQn20U1tpgk//yyr9wcPygrsnDkSo45ePcKLy0ay9NBScnnmYnzj8fR/vD8HLx50/f3dB5WI7NtXzk1lUhqAlVIPFh0tsxLDkDOYX30Vtz/npBZujrBybxBDlx8kLMoGQFBwGEOXHwRIsyDs7y9JVZs3S/nIRYtkRX/l0eXU+Hwce8/vJatHVobXH87bdd8ml2cuwMX3dwMDZQ197lzp3+znJ0G4SxeXr1LlCIZpmg57sho1api7du1y2PMppZLpyhX4/HPZn1u4UOozX7ggSVX58jl7dGmu3sQNBAWH/ed2n1xebB3SOFWf6++/pXrVqlXwyCOSXNWzJ1wOP0f/n/uz/MhyQBolrOq4ipZlWqbq86c604Q//5Rl5hUrJFegdWv50JYJa3obhrHbNM0a8d2nM2ClVJzdu+XCuWiRLBs2bhyXwVywoHPH5kDn4gm+D7s9Oc6elWD79deyIjt2LAwcCGHGFd77YxIzds4g0haJgYGJiWma7L+w33UDcAYqEekoGoCVUiIyUvZ0w8KgWzeZsVSo4OxROUXhXF7xzoAL50p58ZCrV2HCBIlTpilBd+hQcM92g8n+H/Hx9o+5HXWbzpU707psa7qs6EKkLdJ1i2hkwBKRjqIBWKnM6tw52Zv77TdZMvTwkHXQ8uUhVy5nj86pBjX3u2cPGMDL3cqg5n7J/pm3b0tP3kmT5PRWly4wejQcj95A51+kUcKtyFt0KN+B0Q1HUz5/eVbuDaK0MYmzkbso6l2Di1eKQdHUeIUpFF+JyOeekyNoGaRKlSNoAFYqMzFN2LZNLpzLlsl5l1Sbh0AAACAASURBVJYt4do1yWKuW9fZI3QJMYlWKcmCjsmiDroajvV4Sa5vKcP1K1batJHOeaX9Ihj6+1A+3vExII0S5redz6tVX439fvkQUJKclOTmTdI8ESxBwcGSiDdzJhw/LtsSI0dKRnOhQs4ZUzqmAVipzGT9ejnTkjMnvPkmvPGGpNyq/2hXzSfZgW7l3iCG/HCQqwcLELy5LNHXs+FV9Brj59kY1DU3X+/7mpYzPuDMjTOx32NgcP7W+divp6w7es8MHCAsysaUdUcdH4APHZIPbZmkRKSjaABWKiM7dUr25/Lnh0GDoFEjmcF06ABZszp7dBnWsJmX+Pen2kReyIV7vlvkf24nnqUuMP/cDubPXMqxa8eo5VOLpkW689WhCdjNKEzcMCLj9twdkQj2UA8qEdm3b4Y+guZIGoCVymhME37/XWYsq1eDxSJLhCB7dV27Ond8Gdju3dKX9+/fq2HNEUqOF77GKL2BKCwEu20kKvoUldwqsarjKmy3q/PeigAK2MYSbjmIp70SX2/04NE8QbSr5pOmiWAPdfly3BG0TFoi0lE0ACuV0QwYILOWAgVg2DDo1QuKFHH2qDK0Y8dg+HBYulRilO8z/xBafjWXvIcA0WCA1Z6Psm7D2Nf7AyyGhXoTNxAWZSML5chiLwdAmD1uiTktEsEeSktEOpz2A1YqvTt8WI4MHT0qX3fpIi0Bz5yRfToNvmnm/Hk55lq+PPz0E7z/vvTqfX3g31zzngKGBF9Mg9xmSya17IPFkMtuQkvM7ar5MKF9JXxyeWEgRUAmtK+Uuvu/ERHw7bdSy7tWLVi+XFZLDh+GX3+FNm00+KahBGfAhmHMA1oBl0zTrHjntilAayASOAF0M00zOC0HqpS6i80my8szZkhilYcH1K4tpf5q1JB/VZq5cQMmT5ae8ZGRssjw/vsQZN9Dx9XDWXt8Ldk8cnE7yg3TtGMx3Hm7wXP3BM/ELDGnJBHsobREpEtIzBL0V8AMYMFdt/0GDDVNM9owjEnAUGBw6g9PKfUfNhtUrAhHjkDRonKmpWfPDNcMIa0lp9lCeLisyo4fLye3mnXzp3LbjVQvWZy+m5bzw+EfyOOVh0lNJ9GvVj/2X9j/wEYJDl9iflCJyP79oUkTPbvrBAkGYNM0/zQMw/e+236968vtQIfUHZZS6h67dsGaNTBqlCwJ9uolyTGtW8c1vVeJltRmCzabrOqPHCklJJs3h46D/OmzvTG/7QuHfeDt5s3IBiP5X+3/kdMzJ/DwRgmpcdY4UbREpMtKjf9zuwNLUuHnKKXuFtPCbcYM2LFDjg316CGz3oEDnT26dC2xZ2xNU07ivPeeNE2oWVNqN5d5LJDnlgwkPDockDO8/6vzP0Y1HJWkcaTZEjNIiciZM6VE5I0bWiLSBaUoCcswjGFANLDwIY953TCMXYZh7Lp8+XJKnk6pzGPXLgm0r7wC16/DtGnS8L6oK9QhTP8Sc8Z282apN9GunRyJXbYMVm+4xOrItyg9rTS7z+/GzXDDaljxdPPkmTLPOGr4D2a3S2nIZ56BMmVkX/fpp2HLFtizB7p31+DrQpI9AzYMoyuSnNXEfEhPQ9M05wJzQdoRJvf5lMrQTFOKHYSHywWzXDnZl+veXf606IGF1PSwBKiDB6U5wk8/QeHCMGcOtH8pmE/+mkrXaZ8QFh1G1ypdGdFgBOdvnX/gHq9DxVcicsQI2arQEpEuK1H9gO/sAa+5Kwu6BfAR0MA0zURPa7UfsFL3CQmR8n4zZsgaZ926sHWrs0eV4d2/BwxgvZ2NYv/WZNPP3uTMKQU1KrZZz7Rdk9kWuI2QyBBeqPACoxuO5tF8jzpx9HcJCJCgu2ABhIbK+6dfP2mMoCUiXUKK+gEbhrEIaAjkMwwjEBiJZD1nAX4zJHNuu2mavVNtxEplBrNmyVTr5k0p7Td/Prz4orNHlSncnQB19pyN6N2PcvmvIpxzMxg0CAa+E86k3UMYsvRTQBolfNX2K7pWdYEqYtHR0rVqxgxZNcmSRfZ1+/fXEpHpTGKyoDvFc/OXaTAWpTI2m03WNWvXlipVhQpJFnO/flIIQY+BOFTTMj4cWO3D1K8kUbh7d3jv/Sh+uzyfWt+OIfBmYOxjDQzO3TrnvMECXLoEX3wRVyKyWDGYOFES8/Llc+7YVLLoxpJSae3KFamlW6oUtG0ry4UAzz4rVYhq19bg60CRkZKbVLKkHCtq1gwOHLTRoP+3NFtZjl5relE0R1GmPz0dLzcvrIYVD6sHDX0bOmfAO3dK/e6iRaW0qJ8frFwpJbcGD9bgm47pAUKl0ordDq+9JmcwIyKkE9GHH0oQVg5nt8OiRVKx6t9/oWFDeGnwNjaFzqL1b9v4N/hfqjxShTWd1tCyTEsMw+CxQo85J8kq5gja9Onw11+QLZu8l/r2lQQ9lSFoAFYqNUVEyJGPmMzl0FBZInzjDahQIeHvV6ni7ipXhXJ60SR7ZVZ/kY/9+6FKFfj5Z5PDOT+h129vY2JiYDCm0Rjeq/9ebK1meHghjTRxf4nIsmXlCFrXrloiMgPSAKxUajhzRi6cX3whS84nT0qVoUWLnD2yTOfuDOeIoFzs+e5R/M/m5ZEi0Sxc6EbhOn8yYuNwNv+1OfZ7LIYFq2G9J/gm5nlSpYpVfCUiW7WKKxGpR9AyLA3ASqXEiRPwzjtSLgkkqapvXykTqZxiyrqj3DzvxfXNfoT9UxCLdwS5mwaQ/8nNfG1fya8LfqVQtkK8U+cdZu6cSaQtMsl7vEktZRmv+0tE5s4Nb70lJSJLlEjqy1bpkAZgpZLq5k24eFEqDWXLJkkygwfH1WdWThMYCAcWlSbkYFHw3UKWzpPxLmQl3HMLAfbtnD+Xl6nNpvJGzTfwcveifbn2ydrjTWwpy3jdXyKyShVZOenUCby9k/OyVTqlAVipxDp0KK7oQfXqsmz4yCNw+rT2THWya9fkRM706RARVQTvJssJrfcSEUQSAUAWilpf5dCAaWTPkj32+5K7x5uYUpb3sNth3TqZ7a5dK++X556TZea6dTULPpPSAKxUQn77DcaNg02bpOhBx46yzBxDg6/ThIZKjtLEibIw8corUK7tHkbtHwpEggGYkNfenmnPjL0n+KZEYnr5Ag8uEfn661LnUmVquruvVHwuXJCrO0jf3VOn5CofGCgX1Jo1nTm6TC8qShKFy5SRYmL168OGHRfJ1XEAI/+uR7T1XzCsYFqwGFl4p0HHVO06NKi5H17u937wuqeXb0CA7OUWKQL/+58UXvnuO1ktGTVKg68CdAasVBzTlDrMM2fCDz/In6+9Jnu7b7yhM10XYJryn2bYMPjnH6hTB+Z+c40t9ik889s0IqIj6Fa1G+83eJ+gm0FpdoY3vl6+7zYpRdt/d8Bb0+NKRL70klQ60xKRKh6JasaQWrQZg3JJdrskwcycCQcOQM6c0K2bXDhLlXL26DKl+I745Ljuw5AhkvNWvjy8OPR3/rRPwj9wG2FRYXSs2JHRDUdTJm8Zxw728mX4/PO4EpHFi8sHth49IG9ex45FuZwUNWNQKsO6ckXK+FksEoANQ9Y1X3oJsmZ19ugyrfuP+Px71J3On2fh9kmpxjjnyzBWWvoy8sT8O99hYUC1T/ikzZuOHejOnZJUtXix1Lds2lS+btVKV0tUomgAVpmLzQZr1shsd8sWmbHkySOZqXnyaDaqC4g54hN13ZvgzWUJPeyDxTOS4s8c4q1hfzL0z9FcC78IJrFJVt/t2kvDokGpus8br/hKRPbsKaslWiJSJZEGYJU5XLsmndU/+0yqVvn4SPZOTJUhXSp0GWcC7QRvq0DI/mIYVjvZ6xzB/YkvCPL6lgG/XiS7UZHckc8S7D4P04zGwA1rdIXEncFNrvtLRPr5SRDu0kVLRKpk0wCsMi7TlEzmrFnlovnee9IQ4aOPpCGCm779XcmNGzB1Kpyb2wi7zcCz4QqMWjMJ9QjAZrlMVqMMP3aaT58vbYBBFrMU4ZaDeNorkcVe7sFncJMrvhKRMe0jtUSkSgV6BVIZT1iY1GCeOVN6zn3/vcxYTp3SSlUuKDxc8pfGjYOrV6HeU+EcrvQu17LNBsME0yCfrStz243l6TJF8Mm1gaDgMLLYy5HFHrfs+58zuMmlJSKVg+hHOJVxnDwpdZl9fCQDNSICnnoq7n4Nvi7FZpMj1X5+Et+qV4fPftmI+eJTXMs+S4IvgGHQomIhnq1eBEjEGdzkOn5cBuLjI0fPYpLzAgNh8mQNvirVaQBW6ZvNJkuDILV1P/1UslE3bpTZy2uvOXV46r9ME1avlhLI3bpJjYrpK3ZgvtKM3tsbcTr4NO/WfRcvNy+shhUvtyy8UadN7Pe3q+bDhPaV8MnlhQH45PJiQvtKydv/tdslAe+ZZ6T13/Tp0KKFJOjt3Ssf5LQ+s0ojeg5YpU9XrkjA/ewz+Phj2dO9ckWOg2iVIZe1ZQsMGQJbT/uTt8ZGOrYswpls37P6n9Xk987P0CeG0qdmHzzdPPE/659mhTTiLRHZq5eWiFSpTs8Bq4zjr7/korlkiSwxN2gge3QgZ3qVSwoIkKTzNWsgb1V/3Ho25irhzDwHWd2zMrbRWAbUHkA2j2yx35PcRgkJDmTGDPjmG0nQq1MHPvhAGiN4eKTucymVAA3AyvWZppzPtdulEcLly7I0+MYbUKGCs0eX6SSlEf3p09J74Jtv5LTOoHGnWebZi6u3wu88wqClb0+GPTks7QYcHQ2rVkng1RKRyoVoAFau68QJWWJevRr27QNPTykEXKqUnr10ksQ2or9yBcaPl8UKw4Be75wnqvY4Pg6YS/RNOzHpJwZu/HW4OCv3pkERjftLRBYrJg01evTQ1RLlEjQAK9dis0lSzKxZ8MsvUtLv2Wdlz65gQahWzdkjzPAeNsNNqBF9SIhsyU+ZIqd5Ona/SvanJ/H14RlEHYoir9Ect4gO2IwrsWd4sZdN3SIa95eIbNJES0Qql6QBWLmWHTuk2EGhQjBypGQxa1KMwyQ0w31QsYugq+HMnAljxsDFi1C38+/Y605k5TV/wg6G8XLllxnVYBRNJx/BBNzM/Pec4U1xEY34SkS+9ppsU5Qvn7KfrVQa0QCsnMc045Kq8ueHDz+UpJg1a+T8rru7s0eY6SQ0w72/Eb1pQujhwtza+ij9rkG9hqE8Me5dfgicCZfAalj55tlveLnyywAUznU6cY3sE+v+EpFly8K0adC1q25TqKSJjobt2+GJJxz2lHoOWDleaCh8+SXUqAG1a8PKlXHB1jDkTKYGX6d40Ew05vaYIhimCWEn83Hh6ye4sroaBfPb6DNvJidalZLge5czN87E/j1VimiYJmzaBM8/D76+stlcpw78+iscPgz9+2vwVYl35Yos3ZQoAfXrw9GjDntqDcDK8QYPlg4yERGy1xsUJMkxyukeNBONub1dNR+6la1O8LK6XPr+cYiy0HT0dOhdkdln+lE2b1lmtZwVW0TDw+pBQ9+GsT8nRUU0bt+WmW6VKtCwIaxfD//7nyTrrVoFzZppfWaVOHa7NGgBKUI+apRsVSxfLuVrHUQLcai0FR0NP/0kgXb0aJnxHjsG58/Lp01t/+dS7t8DBpmhTmhfiXLePgwbJonoOStupVSnGVz23MbZW2eoUbgG4xqPo1nJZhiGkbpFNE6ckG2KefPkYlmlisxyO3XSKlUqaa5ckQIsc+bIEcaVK+X2c+fSLNdEC3Eox7t4UerozpkDZ89Kfd2LF+W+MmXkX+Vy7s52jsmC7l6tHD/PLkSHeeDpZdJs1FR+ZzB7okyMKIMJTSYwuN5gjLs+TKW4iIbdLkvK06dLVrzVKsUy+veHunX1g5tKmp07pUzt999LZnz9+vIBLoaTEj01AKvUZ7NB1apw4YIcAfn0U8ls1vZ/6UK7aj60q+bD9eswaRL0ayf/SdsMWM/p0sP47dKO2MdaDAumad4TfFMkvhKRI0ZoiUiVdDduyAqJu7scaVy9Wt5HvXpBxYrOHh2ge8AqNdy6JZmozz4rMxerVb4+fBh+/11u1+CbboSFSfOfUqXkzwadt1FjemNW5GjK5fBzDKk35IF7vCDL2PUmbqDEkJ+oN3EDK/cGJfykAQHS7q9IEdnXLVBAWgKePi37cxp8VWLt2RN3fHHVKrltwABZZp4+3WWCL+gMWKXEoUNSZWjBAgnCVavKMnOhQtIcQaUr0dEy+Rw1CoIMf4q/+B1Fy+1m3XV/CtwqwKctPuX1x17H082TNn5t4t3jTWylrNgn/PFHuShqiUiVEtHRUu909mxZbvbykvdSuTtnzV00K16TsFTy/PEHNG4sBexffFFmL7Vr696ci4uvylXbqj6sWAHDhsGRI1C63SJOVn0FOxJEe9fozdRmU8nqkTXBn19v4oZ4z/n65PJi65DG8sX9JSKLF5eCGT16QN68qfp6VQZ39aq8Z0wTHn1UVt/69IFXXoFcuZw9OkCTsFRqCAyUIyD580sizBNPSM3Bl1+W25RLeFgZyfhmp29+FMTgffn5J8CDUjX+pdGno9l4fQEm8sHcalgplqNYooIvJHCO+P4SkU2baolIlXSRkbBihWxz7d8vxxi9vORs+COPpKtJgAZg9WCmKWctZ8+WvRS7XWYpIIkNAwc6d3zqHgkt/95d5SryYg6ub/Ij/N8CuBc+ScMpU9ka9gWBNy10rNiRFUdWEGWLineP92Hur5TlER1Fy6NbeG3/zzDpcFyJyL5945YHlUqMc+fkA9uXX8KlS1KE5d13JUMQJGEvnUkwABuGMQ9oBVwyTbPindvyAEsAX+AU8IJpmtfTbpjKKfr2leCbNy+8/Tb07i3VYpRLSqiM5LngMKKCvbj6dwgR5hbwvE6WHiuILPIdW8Ls9KzWk+FPDscnh0+yz/EOau7H0OUHyXn1Ii/vW0un/b+QL/QGt4qX1BKRKulsNggJgZw55TjjpEmyYtKnj5SrTeeFVxIzA/4KmAEsuOu2IcB60zQnGoYx5M7Xg1N/eMqhdu+WgDtkCJQuDV26SIm/55+XVoDKpT1s+ffiRQj/szIXgs7AK63ALRwwicAgv7Up2/t9RsnccRWAknWO1zRpd/M4NbfP4JENa7HYTbaWq4Olfz/q9Xox3V8slQNduCCFV+bOlQpnn38OtWpJEE5mRnxS+lg7SoIB2DTNPw3D8L3v5rZAwzt//xrYiAbg9CksDJYulUpVf/0l5+ZatJAAXLu2/KvShfuXfwHsEW6YB8pSqhSE2nLj1q0z0W5hYAAm5La3ZW7bGZTMnYIL0e3bcmRoxgw4eBCf3LllxaRPH+rriolKii1bJCt++XLJbG7SRGoIgOztpiD4Jjo734GSuwf8iGma5wFM0zxvGEaBBz3QMIzXgdcBihUrlsynU2kiMlIC7blzkkE4bZrMenPmdPbIVDLELP+GRdkwoy3c2leMm/5lsEWYVO05jcAS47kSfhGwgAkWw513G3RL/gXoxAn54DZvnhTQqFJFqp9piUiVFMHBcs0xDFiyBH77TRI9e/eW7lapIKHtGWdJ8yQs0zTnAnNBjiGl9fOph7DZpC7zn3/C1KlyhGjoUKmJ2rBhusoeVP/VrpoPNhsMmXyDU7/6En3LA592nxNdawL7ws/SsGBDxjVejoGR/DrNdjusWyezXS0RqVJi507Z8lq8WCpVPfmkHEKfPFmymlNRQl2+nCW5AfiiYRiF7sx+CwGXUnNQKpVduCCZg3fXZR46VJKr+vVz9uhUKjBN+PlnGD3Uh+PBp8j1wge4ldhCUGQQj+d9nIWN59O4ROPYkpFJDrxaIlKlhogI+PZbCby7d8tKSefOUrwH0uwceHzbMzG3O1NysyJ+BLre+XtXYFXqDEelul9/hWLFYPhw8POTvZVTp7TgQQbi7w8NGkCrViZBJSZh9KhPsM8SrkaeY3LTyfj38KdJySbJq9ccECBLgT4+cSUiv/tOS0SqpLl5M+7vQ4dCeLisopw7B3PnsjLEO+nlS5MgVfpQp4HEHENahCRc5TMMIxAYCUwElhqG0QM4AzyfloNUSXDzppRkK1hQlgZr14Y335Szl37OfbOp1PX33/Dee7BqlUmux36j+NjhnI7eGXu/xbAQbY9OeuCNjpZz3zNmSIlIT08p69e3r5aIVIkXGSnt/mbPlg/9x49LudFdu6Bo0djtCkckSMXX5csVsqC1FGVGsX+/vNG//VayUjt3lkCsMpyzZ2HkSPj6a/Asu4UCHYdxij8plrMYL1d6mU+2f0KkLRIPqwfru6xP/HLzpUuSRKUlIlVKBAXJe+iLL6Q2vK+vdCAaODDe44yJKl+ajmkpyozujTfkDe/pKRmoffpAzZrOHpVKZVevwoQJMG2FP7byCynw3i4uuO0gPFtBptefzmvVXyOLWxZal22dtCQrLRGpUspul2Vlb2/pRjRhAjzzjFyLmjd/6BlwV02QcgQNwOnRyZNyQP2ddyBfPmjZUhrcd+0KefI4e3QqmR5UKOD2bfjkE0kOven7LXTuChY7F4Cni3dh2cuz8XaPO/aTqCIaERHSnHz6dDn/rSUiVXJcvizH0ObMkW2KsWPlenTypKygJIKrJkg5ggbg9MJmkzTXWbPkGIjFAo8/Lr12W7Vy9uhUCsW3Dzbk+wB++T4rq+bn4kLECfK+9D48sijum0wLu05a+DXgOu2qJfLcbWCgXCznzpUl57JltUSkSrpt2yQjftkyWTVp0CCuaI/VmujgC/eeX4/hCglSjqABOD0ICZEm0qdPS7r+iBEyW/FxbgKBSj13FwowTQg9UoigzX4ctV2n4AuDsRadx3XTind0A8Ks/phmNAZuWKMrJFxMwDTl7PeMGdJFxm6XD239+0ulIS0RqRIjLCzufO4nn8hEoFcvyZIvXz7ZP9ZVE6QcQQOwK4q5YO7aJSX9smWTvd0aNaBNG+lEpNKdh9WijdnvunnhPDcvnsV2ohaWBl9C5S+46mbS+7FerPrzcazkIcJ2mHDLQTztlchiL/fgvbL7SkSSOze89Zbsy2mJSJVYBw/GJXju2CFbFB9/LP12syauTWVC2lXzyRQB934agF3JjRuwYIH0ufz7b9nf7dVLAvCECc4enUqBhI5a5AgpwJEdV4hs/ib4hsPjJnYMClibs6P/bHxz+bL3gGSLZrGXI4s9bp/2P3tlJ07I8uC8efKe0hKRKqkiI2V5edYs2LpVEjxffDHuw7+uvqUKXXtyFWvWSFGDN9+UT5Xz5smSc7Zszh6ZSgUPqkX7wbdneOEFOPB5OSIrfAhuYWDI0cDc9nbMafUFvrl8gQSKCdjtUhrymWckIW/6dGmqsWUL7N0rR4k0+KqERETIn7dvQ8+ecoxo6lTJHfjqK6kdr1KNzoCdJaYLUfHiUoe5WrW4I0SPPebs0alUdv8ycXRIFm5sLcPpgPwcqvsx3kMnEGpc5t5GCa/esywX317Ze3UL8sym7+EFLRGpkslmkw9vs2fDlSuyzJw7t2yBPfqo5gikIQ3AjnbsmCwxf/UVXLsG3bpJAPbxkWVClSHFHLWwh7txY0cpbu4pApXnY3nnAyI9L1C/RBPGNh6LaZoPPcMbu1cWECB7u2O+gdBQaYTwwQdS/czDwwmvUKU7ly7F1YiPSfDs2VMqobm5pSixSiWOBmBH6tVLjn+4ucnxoT59JPiqDG9AQz/6Db/F5VMXoNp78ORm8L5AmdyPMbv1dzQq0Sj2sQ88w3t/icgsWeTsZb9+WiJSJY5pyozXzQ1+/FFqmTZuLMvMbdtqgqeDaQBOS0FBMtN96y1J369dW2qg9ugR1/1DZWjR0ZJXN2JkYS77ToDOw+/s8Rp0Lf8+8zuMTrhW8/0lIosVg4kT5X2UL59DXodK527dkizmWbPkCOObb8qHt3r1tPCKE2kATm12O6xfLxfLH3+UT5s1akg5tm7dnD065SCmKZPVoe+ZHIn6Be/nh0POPbH3Ww0LfgW9Hh587y8R2aSJlohUSRNzhOibb6SeQLVqMgkAScrT4OtUGoBT06VL8ony+HGZmbz9tiw7lyzp7JEpB/rzTxgyBPzPbcLrmeFQYAsFcvnSudJwPvT/MLZRQkPfhv/9Zi0RqVLKbo9LnOrfH7ZvlyNEb7wBtWrFdiFSzqcBOCVMUzIGjxyBV1+F/PnhySelT2qHDrJHpzKNAweg9zh//EO+xb3aTmi+k1zZCvHhk7PoUb0HHlYPWpZpGX+SVWCgJOfNnSv1dbVEpEqq06fl/bNggayeFCwo76n8+bWblYvSdoTJceuWNCWfPVvaABYqJG9+TWDIlP79V07+fLv/W3hWGiUA9KvZj8nNJuPl/oCi8vGViGzdWma7TZvq8Q+VMLsdfv1VrkVr1shtrVvDlClyHlw53cPaEer/4Um1ZIkcGerdW77+7DM4elSDbyZ06RIMGABlax/nu4iXof0rYNjv3GshOCRb/MH39m2ZqVSpIlnw69dLot6JE7Jx/NRTGnzVw8VMnM6ele5D27fD0KHyaXDlSg2+6YQuQSckPFz25CpWlASG8uWhfXsJwI8/rvspmdCtW/DRRzB5zhnCao6BPvNxt3qQJboxIcaW2EYJmw8WYGWZoLhiGidOSBbqvHkQHKwlIlXSmKbkBcyaJR/ili2TQj4bNsg5cD3/ne5oAH6Q48flgPr8+dIJfeBACcCVKsnRIpXpREbCkBn+zPljDaGex7C8tgo3N+hd4w389z3J5dteZLU8HdsoAXtZpq49TLsLB2SZee1ayV5u316SY+rV0w9wKmGhobBokQTePXsge3bJDYhJttJaAumWBuD4vPyy7PFarXI4vU8fOayuMiW7Xa5/b89ax8UmreCxaDCgZdlWzGw5k2I5i1Fi008AsY0ScoSH0OHgKl7Zuwaun9cSkSr5pk+XtPpKlWSv9+WXJQi7iId1+VIPpwEYJAN10SI5NmSxkhBDfgAAIABJREFUyEz30Uel0IFeLDOF+C4ibav68Msv8O77twjI+gk0HgcWCb5gIZe1IsVyFgPiSk2WvXyKLnt+4tlDf5A1KpwDxSvAzA+1RKRKnOhoWL1aZrs9e8rxoe7dZYn5iSdcbsUkoS5f6uEybwC22+H33+UT5erV8nWDBnJO7p13nD065UDxXUQGTjvDe4eyctjtSyzNJoLnFbzsFQnjKJi2e/d4Kz3CR24nYPEMHj99gAirO6vKN2BJrba80vdZKuuFSCXk/HnJB5gzRyroFS0qwRjkGFH+/Gn21CmZwT6oy9eUdUc1ACdC5gzAJ09Cs2byZ/78MGiQFDvQghmZ0t0Xkagr2bi2tQTh2X6CBh0g+zka+TbjUlA7bt4sToTlcOweb56QR7g4ZAT8/SuPBwYSWqgIM1u8xpdlGuJVuKAuxanEe/ppOdLYvLnMflu2lHrNaSylM9j7u3wldLu6V+YIwKYpfVEvXpQCGcWKScu/ceOkKYIWzMjUzgWHEX3Tk6uHgwnPPx+e2QrZL5Alujzrun5HA98GlBgSt8dbM8hK1z1raHX4T7LYouXM7owZeLdqRV+rlb5Ofj3Kxd28KXWZFy6UM7xZs0qS3iOPOPz4UEpnsDFbL/HdrhKWsQPwjRtSA/Wzz+DQIdnXfe45+WS5dKmzR6dcwLVrELGtPOeNb6HRSGmUYBrkjuzBo9k60cC3AQDFs1mpuuN3Xt29mqrn/yHEw4vFVZqzrsFzfPeR1vhWiXB/XeYaNWS5uWxZ2d91gpTOYAc197tnBg3g5W5lUHO/VBlfRpdxA/CcOVLcIDQUataUvpcvvuhySQzKOUJD4ZNPTMZ9/zOhjw+HQvsgtiicgZvVxrstHo0tEblu1mdkuX6VE3l8GNm0Fz9UbIItW3YmtK/kzJeh0ovDh6FyZVlt69RJ6jLXrOnsUaV4BhszS9Ys6OTJOAH49m3pGtOgAZQuLbPdTp3kCNFjjzl7dOohHHmMISpK6mC89/kfXKs6DNr5UyRrSer6DGTZP7Oxm1FYDHdGP1KaduMHxpaIzNK6NVtbdGTw9fwE3YzQC416uLNnpdpZZCRMmiSNNObPlzKRLlSXOTVmsO2q+ej/B8mU/mtBHzokS8wLFsjeyrhx0mRapQv3J4GAXAAmtK+Uqv9TmyaM/8afKesWcCPrLvDZRf4sPoxtNoJuVbvhbnXH/9gGNq6eTsMfD1Bn00nInVuOgvTpAyVKpNpYVAZlmlJWdNYsKSlqmpJzsmSJS6+86TnetPWwWtDpNwCbJrRoIUkMHh7w/PNyoaxb16Xf7Ope9SZuiHcJzCeXF1uHpE7xk/Xr4aWp87hU6zVplGDCMyW6s+zlmXi6eUrVs5gSkTduSInI/v21RKRKmg8+gJEjZYbbs6eUq/X1dfaolJM9LACn3yVow5Azu82aSSvAfPmcPSKVDGl5jGHPHug/+ijbPEZC7SV37fFa+OuEjT2ffU/dtYvjSkQ+9xz066clIlXi7NsnH9w6d5Y2pC+9JCslzz8Pnp7OHp1KB9JvAAYYM8bZI1AplBbHGI4fh/+NOs2akNFQ9WsM0xOvyCaEuW/GNKNxsxvMXraGuoe/1hKRKmkiIqQJwqxZsG0beHlJjsmTT0ruSenSzh6hSkfSdwBW6V5qHmNYs8+fEQtXs/fMP+D3I1aLhV7VB7B622M8ejmYuv+GYNh38dQJG+5mAd5s3Z1py8ZqiUiVOKYpq24HDkig/egjWX3LndvZI1PplAZg5VSJOcaQUJLIjRvQa+ovLDFaQ9ZoKActfNvyeZtPKPLHbnp8P4bqJ/cTYXXnx3KNGNukFYcKlsYnl5cGX/VgMeVqlyyRjGarVXru5skjxVe0Z7NKIQ3AyukedozhYaXyWpTz4aOZNxm7/iPCHpsQ2yjBalh58mwYRarUh8BAHi3ow4eNu/FthaZc984JaLEA9RDXr0vL0dmz4dgxKVd77JgcbezY0dmjUxmIfoRTLi2+UnmhETbennCBQh0mM+xSCcIeH03JbDWx4obFDh6RNhrO/RX8/GDlSrwDT1Nq6hi8CxfEQDKsU/uYk8ogAgLAx0eK+OTPLyUjz56V4KtUKkvRDNgwjP8BPZH80oNAN9M0w1NjYEpBXDZ0hOUwYcZBOFuXW4G3sFfvCdkvUCvnU/S7VZnK01YTaolmXSk3bG61GNuqK6+9/kxskNViASpeMUlVoaHSkKVcOcmE79RJ2pIqlYaSfQ7YMAwfYAtQ3jTNMMMwlgI/m6b51YO+J00KcSiXl5KD/vUmbuDkzT1ccB8ORqTcaJjkCKnEmizVqD93LVy+zIk8Piyo3oofKjYhJIuc3U3Ns8QqgzlzRsrVfvEFXLok9QO2bNHjZyrVpeU5YDfAyzCMKMAbOJfCn6cymJS2O+tQshyD1k6H4hFgACY0O1mAtQsPYbUHQOvWvOJZky2+VTCNe3dUtCWaitfHH8f1/G7VCvr2laQqDb7KwZK9B2yaZhAwFTgDnAdumKb56/2PMwzjdcMwdhmGsevy5cvJH6lKlx7W7uxhTp0yadznR976owU235WAFLHyiobhO25jfettOHECVq3iZLW6/wm+oC3R1B3BwTBtmjREAJntDh4s/cBXrfp/e/cdH1Wx/nH8MwlJiLRIjdIRASuCdCwoCKIgzV5AUUEvIKKA+FMRy7V7UQRRRFHAC6gURSkigqI3IEgLRToSOihFBUmb3x+zgRA3hSTk7G6+79crLzZ7zm6ew9ndZ2fOzDPQurVGNIsncv2qM8acCXQAqgNnA8WMMXdm3M9aO8pa28Ba26BcuXK5j1SC0qlWutq/H7oMmEuNF5oyL7YDFYtvY9yMaBZ8AM+vjWXu2Y9zxdK98Morx+szD2hTm+iI8JOeR6OchZUroWdPN6iqb1/44gt3f+PG8MILULWqt/FJoZeXLuhWwBZr7T4AY8wUoBkwPj8Ck9CQ00pXc9fHMWjCRyzdvZjU2KXEUpRnv4C7V/1BROeb4D+9uSyTOt9aEk1OYi20bQuzZ7uSkLff7rqZ69f3OjKRk+QlAW8DmhhjzgCOAi0BjbCSk2RX6SoxEe576yPGHe7u+pgrwMNx8OKqUhS97zGYnrMSkRrlXMjt2OGWjuzVy31Ja9bM1Ym/5x5XOEMkAOU6AVtrFxljPgOWAsnAMmBUfgUmoSGz1ukNdSvy2ke/8Pq3fdld42s3wAoIt1C+/S0U/WKsqlRJ1qyF776DESOOr9tMy5ZuKtHgwV5HVyC0lGBwC97lCCUoWQvjpm3glSk9WFtjPmckQZe1YXxysSExDCKLRDK361yaVm7qdagSyNatc6tXrV7tWrj33uuW/6tRw+vICkxBraUteROayxFKUIlLiGPUzImsiJvDqsprCa8KvVaU5Ikmfakwui89j6xn/tb5tKjWQslX/Fu7FnbtgquvhipV4Kyz4NFHXXnI6MI34j2rGQZKwMFBCVhOuw8/eJP7fu1HirFQFdpsL8c7LV+i2jPdXIF7oGmZpkq88k/JyW708ogR8O23cMEFEB/vEu6cOV5H56nTuZa2FAxNfpPT49gxtr45ij7tKtJzy8Mu+RoIN2Fc2a0f1W7rfjz5ivg1YYKbatali1vk+YUXXBJWwQwg83numv8ePJSAJX9t387+/gN5pm0Z6u/qyfCGO6l1uBZR4VGEm3Aii0TRoloLr6OUQGQtxMXBnj3u96JF3YCqadNc0YzHH4fy5b2NMYBo/nvwUxe05J21xE0fydczRnDg1zVMuBD2Xgl1fmvM2C4jaNfwUuIS4nSNV/w7csS1dkeMgGXL4Nln4amnoFMn9yN+af578NMoaMlSltMc/voLPv6Y7//7Mq2u2ExSGGCg9OELGNpuFF2vauZp7BLgrHWt2lGj3Bq8F17oViK64w4oXtzr6ETyhUZBS65ktpBCsW1buOa7KSSPeZ9JlQ/Tq204SWk9YTaMto3aKfmKf6mpsGQJNGrkruVu2uQWQujdGy6/XNd3pVBRApZMpZ/mYGwqV25eStelX3Llv5cw+bwwHr2rOL+WAQ5WhaIJYFIwFOGHVRWYVmuHusLkhAMH4MMP4e233YCq9evh3HNh0iQthCCFlhKwZGrnwaOU/PtPboyfy13LvqTagV1MOr84VXtWYEfsHoocjKX8tsFElalFUtI6/g6Lp2jqRZBaS3MRxdm5E4YMgfHj4ehRaN4cnnvuxEIISr5SiCkBi3+rVjF03ju0XjaHFRWO8WirGOaXr8iBcjvgUGlujnqf0c935eLnZmOBqNTziEo97/jDNRexEEtKcgUzqlSBiAj47DN3XbdXL7jkEq+jEwkYSsAh7pRqxSYnu/VRhw+H+fNpHxlFn2vOZWSD1RB2EOxBim7tzshuT3J3C7cUYE5XO5JCYOdON6Bq1CioVg3+9z8oV87dX7So19GJBBwl4BCW2SAqODGFYdqyHbw3eREtFnxO1xUzqXBoH1StytIh/eiatJ7VEV/B8YHyYdzYtvzx5AvZr3YkhcCSJfDaazB5svsS17atG1RlrRtUpeQr4pcScAjLrlbs/PFfwStDmbLmO6JSkllQ9RIeaXsri2rsZpN9A1KLU+1Ad3aVnUByaiKRRSL5V9MbTno+zUUspP76y12/jY6GxYvd2rsPPQQPPgg1a3odnUhQ0DzgEFZ90FdkPLuRyUlcv+4Hhu5dAD/9xJ+R0bx8eT3GXlycXSnHSIr5AVIiqbqnN6PveYxWzcqoiIacsHEjjBwJH3zgSkM++KAbXGUtnHGG19GJBBzNAy6k0l+fjT28nzuWz+S2FbMoe+QQ1K7N0616MqbumSQUfxXwTTfafANlUu9n6/h2x5+naWUtlFCoWQszZ7qxAbNmuRreXbpAw4ZueyFciUgkPygBh7ABrWsx5Y3/cstPX9BmfRxh1jKvVmOKPvwQ53dtzcghA9gX/TqQAgawYZSsHEOtEmrJCHDsGERFueu4zz8PW7a4he579ICzz/Y6OpGgpwQcinwlIjsOH07H+HgOR5fgg4Yd+frKznS5qSFxv37KyJfOIbn4Idh6NVT5AUjGUISYsEs0gKqwW7nS1WX+7DP45Rc3knniRIiNhchIr6MTCRlKwKFk40ZXaeiDD+DQIahbF0aPZvUVNdiz8wdK7vmJ7tP78Xf4Por82p7etZ+j+f3leGH+ZBKOLKHyGQ14tm1nDaAqjJKS3KpDw4fD99+7kct33AGJiW57lSrexicSgpSAg11qqhuBOny4u06Xdn2uTx9o1owF236g5diWJKUkuW7mPQ24rdQXvDW0CWXKuKe4tVkfTw9BPJQ2VWjrVrj5Zjd/99VXoXt3KF3a6+hEQpoScLA6eNDV1h0xwrV8Y2NPuj6XkprCmMXj6TOzL0kkHb/GO7B9Z15u18Tr6MVL1sLChe5LW3Kyq8d87rmucEajRu5LnIicdkrAwWbVKvfBOW6cW0e1WTN45hm48UaIjMRayyfxk+n3+WB2payB/ecSVvovTFgKkRGRdKzbwusjEK8cPeqS7VtvwdKlULIk3HvviVZwU410FylISsDBIEOJSKKi4PbbXbWh+vUBsNYyY/1M+kx5ki3HlsK+OtTZ/Qnv9etCeJVFmscrMHQoPPEEnH++Gytw111ad1fEQ0rAgWzvXhg92hU+2L7drSDz8suu1eK7gBuXEMeY5WOY+8siNh9ZCQeqUWnjh4x44A7aX1/Et7yq5vEWOtbCvHnuS1vXrtCxo3vdNG0KLVpo3V2RAKAEHIgWL3YfnBMnulGorVq5bsP27U+6Pvf+0vfpMb0HqaSChRLx/Xmzw7/p+nqkLuMVVn/84S5PDB8Oa9dC2bLudQNQoYL7EZGAoAQcKI4dg08/dYn2p59c1+D997sl3M4776Rd4/fE02/6U8zd8blbKMGAMeEM6FOae1ponmahdtVV8PPPcOmlbpDeLbdoMQSRAKUE7LXt2+Gdd9wSbvv2Qe3aLgl37eoGyaSz4bcNDJz5NNM2ToRjJQhfdx/m4o+xJpHI8EhandPCm2MQb6SkuKlnH30EY8e6kpDPPw8xMdC4sbqZRQKcErAXrHXFDoYPh6lT3Vze9u3doKqWLd0qMz5xCXFM+2UaK3etYfbmmdikKMIWP0b3OgN4fmRpNid21wCrwubAAVds5e23YfNmVxZy/XpXeOXaa72OTkRySAn4NJq2bMdJy/QNuqIy7Vf5BsbEx8OZZ8Ijj7gVZapX/8fjv1z/JZ0mdiLZJruu5tU30Sl6GK+9E0uNGm6fChpgVbhs2QIXXuimoF1+Obz4InTqBBERXkcmIqdICfg0mbZsx/GF6qsc2EXXuV9y5ZBv4NhfcMkl8P77cOutfpdw+/3o77y04BX+E/c6KTbZFdEgnN431eOtm2ML/FjEQ2lT0LZvh759XaWqgQOhQwf3OhKRoKUEfJq8NnMtjX9ZRNelX9Ji88+khIUxs3ZzZlzRhXdG9vF7fe7wscMMjXuDVxa8zpGUP2BTK0z1BZjwJKIiIrm9aYuCPxDxxr598N57J6agXXCBu0QRHg5PP+11dCKSD5SA85uvROS4116h+oFd7CsWw7Dmt/Lfuteyt0QZ15hNl3zjEuKYs3kOe//ay7hlEzmc/Bus7Ui1Lc/y+mMXEdsgju9+1TXeQmXsWDcCPm0K2ogRcP31KhEpEmKUgPNLhhKRh6teSJ/L7mRW7WYkhZ+4Pnd2zInFy7//9XtajW1FUmqSu2N7I8otnskLvRty991QpAhAU5pVUeINaYmJbgra+edDvXrQoEGmU9BEJHQoAedFxhKRRYu6EpG9erHFVOCbKfEkJaUc3z06IpwBbWqTnJrM+JXj6Tuz34nkmxrGtdU7MmVYQ6Kj/f+5jDIO8hrQpraWEgwmO3fCu++6nz174KGHXAI+/3z3mhKRkKYEnBsZS0RWqQIvvQT33Xe8RGRH367pE+Sjrc8lMfJHzn/raTYc/AX21YYzj0BYCmGmCFc2vuCUkm/aIC+AHQeP8viUePe3lYQDX+/eLvGmpEDbtm75yNatvY5KRApQnhKwMSYGGA1ciJso091aG5cfgQWkn35yLZNJk1y3YcuW7vd27fxen+tYryId61V0CyVsmMH/fXMDK/ctJ2z/+Zh5n1EysiERl08jucTPFE29iI/mR1Kn9I4cJdBXZ687nnzTHE1K4dXZ65SAA9Hff7s537fc4uZ5V6zoknCvXlCzptfRiYgH8toCfhOYZa290RgTCfxzTk2wO3YMPvnEJdpsSkT6M/yn4bzy46skHN5G2KEaMHccnWrdxuY2i/k9Ih44B5LPAeBo6skJNKsu5p0Hj/r9e5ndLx7Zts1VOnvvPdi/H8qVcwOrHn/c68hExGO5TsDGmJLAFcDdANbaRCAxf8Ly3qxZi9n9ypu0W/QlZY8c4o9q51AikxKR/izcvpA+M/qwZNcS1zeQGsHFmz/k3fcup1EjqD7oN7+PS0ug2XUxnx0TzQ4/yTb9IC/x0P790LMnTJvmfr/hBtfNfNVV3sYlIgEjLPtdMlUD2AeMMcYsM8aMNsYUy6e4vGEtfPcdO1pdT6vrmtB13n9ZdnZt7rz5ORrfMYxpzTtlm3xX7F5B+/+2p+n7Tfl521pINWAgLDyVmx79gUaN3H6ZJcq0+7PqYgYY0KY20REnd3unDfISjxw54hZCAFePecsW6N8fNm1y3c9XX636zCJyXF66oIsA9YE+1tpFxpg3gUHAU+l3Msb0AHoAVKlSJQ9/7jT66y8YP951M69aRYnoEoxu2JHx9a5je4yv8lSyzfL66rr963h6/tNMWj2JIkkx8P2/qZDSmN/btCcFt1jCVdVaHN9/QJvaJ7Vw4eQEml0Xc1ocGgUdALZudXN1338fIiNdt3NkpEvGSrgikom8JODtwHZr7SLf75/hEvBJrLWjgFEADRo0sHn4e/lv0yb3wfnBB3DokCtmP3o0jdfEcDTin0u4ZUyKcQlxTP1lKmv2rWHmxpmY5Gj48Qli1vdnyKAY7r8fft4z1+9iCdkl0Jx0MacN8hKPLF8OQ4bA9Oku0Xbu7KYSpdVlVvIVkSzkOgFba3cbYxKMMbWtteuAlsCa/AvtNElNhdmzXWt35kw3erlLF3d9rlkzMIbSL32bbfKbvm46nSd1TrdQwi0U+34Yg/qU5+HP3VgtgKaVM18sIasEml0LWTzy559uYF6ZMm5Voh9/hEGD3IIalSp5HZ2IBJG8joLuA3zsGwG9Gbgn7yGdJr4SkYwYARs3QmwsDB4MPXq45dzSySr5/XbkN17+8WWGLhxKcqpvoQQbTvOadZk2sjxly+ZPuOpiDjDpe0u6dYM334QWLSAhQQvei0iu5CkBW2uXAw3yKZbTI0OJSJo1g2efda3eyEi/D/GX/HpdfTbLD71H1zf/w5+Jf2I2t4bK30F4EsYUoV2rC/It+aaPQwnXY99+C0OHwldfud6SG2901c7AdTEr+YpILoVmJayMJSKjotyHZu/eUL9+jp6iY72KVCi7ja83L2Dvn3vpM38ivx/9naKbO2NnPkuxmHJEnfEZKWUWnXIRDQlwR46cWCZy/Hg3//upp9y0ogy9JSIiuWWsLbhxUQ0aNLBLliw5fX/AX4nIf/0L7r2XU22efrf1O64Zd83xWs1R+xpzbMoImte4lEMX/swfMbv/8ZiKMdH8OOjqfDkU8cCGDe5L25gxruXboIGbz1uihPsSJyJyiowxP1tr/fYUh0YLePFi98E5cWKOSkRmJTk1mbErxvLI7EdOWiih1O4OjH73Utq1gxqP/zP5gqpQBaXUVPj6axg2zA3Ki4iAm292SRdO+YubiEhOBW8CPnbMLeH21lu5KhGZUapN5ZPVn/D0/KdZ/9t6zvirDkT9DWHJRIZHMvn1FlxW1e2rKlQhwFp3DffoUbj1Vnct9+mnXTfzWWd5HZ2IFALBm4AHDnStllq13L/duuWoRGR6cQlxzNs6j4iwCMbHj2flnpWU+vtCmDaVovs6cN9jCynbYD6tap48h1dThILY+vUn6nrHxUGxYq67+cILMx2UJyJyOgRvAu7VC66/3hW2Dzv1ippxCXFc9dFVHEs5BsAZyRUJm/4xiRtv4cl+4fTvD6VKNQX+OYdXU4SCTNrc72HDYNYs1818yy1uTm+JEjkemCcikp+CNwHXquV+cuF/Cf/jril3HU++pIZxdMEDPNj8dp6a5aYIZ0dThILI55+7KlWxsa5yVc+eOTvJIiKnUfAm4FxYtmsZT857khkbZnAGZ0JKBJBKuIlkwostucl/wSoJNmndzOecA337usF4EydCp07qZhaRgBGyCTj9WroxJfcSXXYycbu+oljYmZRa/BKHvu5No/YraXjzfO5o1iLTcpESJPx1Mz/8sNuW1uUsIhJAQjIBp62luz/lew5FTmJr4hbM9mhKrnqCwzMG0PDiUrw0C66+2v81XglCDzzgFr0/6yxX6axHD6hQweuoREQyFZIJ+PlZP/Arb3A0cqG7wxbBTviUxINX8un4YnTpooVqgl5a0Yx+/aBaNbjnHrfYfRYlRkVEAklIJeB9f+3jpR9eYumx4djwZHenAVIt0ZfPpnxZuPHG6zyNUfIgNRXmzHHdzDNmuK7lJk1cAm7a1P2IiASJkEjAB/8+yOv/e503Fr3BkaQjRCbcyLG4m6BzVwhPxISFUyq2MhVLqnB+0EpOdtOF4uPdCOZnnnHdzBrNLCJBKmgTcFxCHF9v+prdf+5m4uqJHPz7IOcm3sTWMc9gf69D6fpbCD/6b5Kil1M09SJiwi9UoYxgs3Gja+k+9BAUKeJKRA4a5FYkUjeziAS5oEzAL839gv/7oQvWuvV4yx5rzBmfjGTTlnp07+4qCi7ZF8Wrs+uy82AtFcoIJtbCN9+49XZnzHCJt3Nnt9j9k096HZ2ISL4JugQ8bdkO/vPdVGxYiu/6bhj7F3SgSdU6jPkS6tRx+1WqpEIZQWflSleXee1aKF8eBg9WbWYRCVlBl4Bfnb2O8OQLoEgUkAQ2gtJ1ShN2Thx16mgpwKCzZQvs2weNGkHVqi7xPv64627WEoAiEsKCLgHvPHiUKM6j7IFXOJqyluKlalC0TKVTXgowfaEOdVEXMGth3jzXzTx9Olx6qVtSslQpmD/f6+hERApE0CXgtKUAixWrQTFqgD1xf06lFepIW81ox8GjPD4lHkBJ+HT7/HN44glYvRrKlXO3H3jA66hERArcqS8j5LEBbWoTHRF+0n2nuhTgq7PXnbSUIMDRpBRenb0uX2KUDLZtg8OH3e2DB9383TFj3P3PPQcV9aVHRAqfoEvAHetV5MXOF1ExJhoDVIyJ5sXOF51SyzWz7upT7caWLFgL33/vpgxVr+7KRALcdRcsXQp33w1FNS9bRAqvoOuChrwvBZjWje3vfskja+HDD9313RUroHRpGDgQbrrJbc/F2s0iIqGoUH4a5kc3tmRw6JD71xj4+GNISYFRoyAhAV58EapU8TY+EZEAE5Qt4LxKaz1rFHQeWQsLF7razJ9/7hZIqFgRPv0UYmK04oWISBYKZQKGvHdjF2qJifDJJy7xpk0fevBBV7UK4MwzvY1PRCQIFNoELLlgrWvV7t7tBlHVrAkjRkDXrlC8uNfRiYgEFSVgyd7SpW5Q1YED8MUX7nru4sVQt64GVYmI5JI+PcW/5GR3Lfeyy1ylqilT3Lq7Kb750/XqKfmKiOSBPkHFv7ffdvWYd+2CoUNh+3Z3zTc8PPvHiohIttQFLc6qVS7BXnONm7N7552uxXv99Uq6IiKngRJwYZaSAl995a7vfvutq0x17rluW+nScMMN3sYnIhLClIALs06d3GpElSq5Yhn33w9lyngdlYhIoaBrwIXJhg3w6KPwxx/u9x49YNIk2LwZBg1S8hURKUBqAYc6a+Gbb1w384wZrlhGmzbQujW0a+d1dCIihZYScCg7dAiaNYM1a6B6AmATAAAKNUlEQVR8eRg82K29GxvrdWQiIoVenrugjTHhxphlxpgv8yMgyaNt21y3MrgSkZdfDh995O4fMkTJV0QkQORHC7gvsBYomQ/PJblhLfz4o+tmnjoVIiOhbVsoWRLeecfr6ERExI88tYCNMZWA64HR+ROOnLLFi6FhQ9fS/eYbeOQRWLvWJV8REQlYeW0BvwEMBEpktoMxpgfQA6CK1oTNH3v2uJHMNWu6kct//w0jR8Jdd0GxYl5HJyIiOZDrFrAxph2w11r7c1b7WWtHWWsbWGsblCtXLrd/TgCWLXOrEFWpAv36uftq1ID4eDe4SslXRCRo5KUF3By4wRhzHVAUKGmMGW+tvTN/QpPjZs2CF16ABQtckr3/fujT58R2LXwvIhJ0cp2ArbWPA48DGGNaAP2VfPPRwYNujd0iRdx13m3b4LXX4N57ISbG6+hERCSPVAkr0KxfD717u/KQU6e6+/r3h40bXRUrJV8RkZCQL4U4rLXzgfn58VyFkrUwZ86JalWRkXDbbXDBBW57dLS38YmISL5TJSwvpaa6Re2tddd0Dx2CZ56Bnj2hQgWvoxMRkdNICdgLCQkwYgR89hmsWOEGVk2fDlWrQlSU19GJiEgBUAIuKNbCwoXwxhswebL7vVMnN9iqWDGoVcvrCEVEpAApAReUlSvdwggxMW4Ob+/ersUrIiKFkhLw6bJ/P4waBUeOwPPPQ926MGGCWwKweHGvoxMREY8pAee3VavcaObx412JyA4dXHezMXDrrV5HJyIiAULzgPPTsGFw0UXw8cfQrRusXg3TpqlSlYiI/INawHnx55/w4YfQqJH7ufZaePFFVyqyTBmvoxMRkQCmFnBubN3qqlNVquTm76ZVrKpVCwYNUvIVEZFsqQV8qh54AN57z3Ur33QT9O0LTZp4HZWIiAQZtYCzk5gIkyZBSor7vWZNGDjQtYInTFDyFRGRXFELODP79sE778Dbb8Pu3TBzprvG27+/15GJiEgIUAs4o0OH4L77oHJlGDzYzd+dNQvatPE6MhERCSFKwOAWRdiwwd0uXhwWLYJ77oE1a04kX00lEhGRfFS4u6DTphG9+aaryZyQAEWLwvLlEB7udXQiIhLCCmcLeMcOGDDgxDSismVh+HAo4vs+ouQrIiKnWeFpAVvrRjRHRcHmzTB0KNx4Izz8sEYyi4hIgQv9BJyYCJ9+6pYBbNzYtXQvuwy2bYOzz/Y6OhERKaRCtwv6t9/ghRegenW480744w+49FK3zRglXxER8VTotoCfeALefRdat4bRo91I5rDQ/b4hIiLBJTQykrVuutC117opROBqMsfHw+zZ0Latkq+IiASU4G4BHzkC48a5aURr10JsLOzc6bZVq+ZpaCIiIlkJ3gScmuqqVG3cCPXru0R8880QGel1ZCIiItkK3gQcFgZDhkCVKm5UsypViYhIEAneBAxwxx1eRyAiIpIrGpkkIiLiASVgERERDygBi4iIeEAJWERExANKwCIiIh5QAhYREfGAErCIiIgHlIBFREQ8oAQsIiLigVwnYGNMZWPMPGPMWmPMamNM3/wMTEREJJTlpRRlMvCotXapMaYE8LMxZo61dk0+xSYiIhKyct0CttbustYu9d3+A1gLVMyvwEREREJZvlwDNsZUA+oBi/Lj+UREREJdnldDMsYUByYDD1trD/vZ3gPo4fv1T2PMurz+zXTKAvvz8fm8pGMJPKFyHKBjCVShciyhchyQ/8dSNbMNxlqb62c1xkQAXwKzrbX/yfUT5f7vL7HWNijov3s66FgCT6gcB+hYAlWoHEuoHAcU7LHkZRS0Ad4H1nqRfEVERIJZXq4BNwfuAq42xiz3/VyXT3GJiIiEtFxfA7bW/gCYfIwlN0Z5/Pfzk44l8ITKcYCOJVCFyrGEynFAAR5Lnq4Bi4iISO6oFKWIiIgHgiIBG2OuNcasM8ZsNMYM8rPdGGOG+bavNMbU9yLO7OSkfKcxpoUx5lC66+qDvYg1J4wxW40x8b44l/jZHvDnxRhTO93/9XJjzGFjzMMZ9gnYc2KM+cAYs9cYsyrdfaWNMXOMMRt8/56ZyWOzfF8VtEyO5VVjzC++189UY0xMJo/N8rVY0DI5liHGmB3ZjZkJpPOSyXFMSncMW40xyzN5bKCdE7+fv56+X6y1Af0DhAObgBpAJLACOD/DPtcBM3HXpJsAi7yOO5NjOQuo77tdAljv51haAF96HWsOj2crUDaL7UFxXtLFGw7sBqoGyzkBrgDqA6vS3fcKMMh3exDwcibHmuX7KkCOpTVQxHf7ZX/H4tuW5WsxQI5lCNA/m8cF1HnxdxwZtr8ODA6Sc+L389fL90swtIAbARuttZuttYnARKBDhn06AGOtsxCIMcacVdCBZscWvvKdQXFe0mkJbLLW/up1IDllrf0e+D3D3R2Aj3y3PwI6+nloTt5XBcrfsVhrv7bWJvt+XQhUKvDAciGT85ITAXVesjoO31TUm4EJBRpULmXx+evZ+yUYEnBFICHd79v5Z9LKyT4BxWRdvrOpMWaFMWamMeaCAg3s1Fjga2PMz8ZVPMso2M7LrWT+YRIs5wSggrV2F7gPHaC8n32C7dwAdMf1qPiT3WsxUPT2dad/kElXZzCdl8uBPdbaDZlsD9hzkuHz17P3SzAkYH9TnTIO3c7JPgHDZF2+cymuC7Qu8BYwraDjOwXNrbX1gbZAL2PMFRm2B815McZEAjcAn/rZHEznJKeC5twAGGOewK3A9nEmu2T3WgwEI4FzgEuAXbju24yC6bzcRtat34A8J9l8/mb6MD/35fm8BEMC3g5UTvd7JWBnLvYJCMaV75wMfGytnZJxu7X2sLX2T9/tGUCEMaZsAYeZI9banb5/9wJTcd006QXNecF9SCy11u7JuCGYzonPnrSuft+/e/3sEzTnxhjTDWgH3GF9F+QyysFr0XPW2j3W2hRrbSrwHv5jDIrzYowpAnQGJmW2TyCek0w+fz17vwRDAl4MnGuMqe5rpdwKfJFhny+Arr5Rt02AQ2ldCoHEd80ky/KdxphY334YYxrhztFvBRdlzhhjihm3DjTGmGK4wTKrMuwWFOfFJ9Nv88FyTtL5Aujmu90N+NzPPjl5X3nOGHMt8Bhwg7X2SCb75OS16LkM4x864T/GoDgvQCvgF2vtdn8bA/GcZPH56937xeuRaTn5wY2mXY8bhfaE774HgAd8tw0wwrc9HmjgdcyZHMdluG6LlcBy3891GY6lN7AaN8puIdDM67gzOZYavhhX+OIN5vNyBi6hlkp3X1CcE9yXhl1AEu5b+r1AGWAusMH3b2nfvmcDM9I99h/vqwA8lo24a29p75d3Mh5LZq/FADyWcb73wUrch/dZgX5e/B2H7/4P094f6fYN9HOS2eevZ+8XVcISERHxQDB0QYuIiIQcJWAREREPKAGLiIh4QAlYRETEA0rAIiIiHlACFhER8YASsIiIiAeUgEVERDzw/8LnxBaq3UvqAAAAAElFTkSuQmCC\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",
    "ax.plot(x1, y2, 'o', label=\"data\")\n",
    "ax.plot(x1, y_true2, 'b-', label=\"True\")\n",
    "ax.plot(x1, res2.fittedvalues, 'r-', label=\"OLS\")\n",
    "ax.plot(x1, iv_u, 'r--')\n",
    "ax.plot(x1, iv_l, 'r--')\n",
    "ax.plot(x1, resrlm2.fittedvalues, 'g.-', label=\"RLM\")\n",
    "legend = ax.legend(loc=\"best\")"
   ]
  }
 ],
 "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
}
