{
 "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:                Fri, 10 Jul 2020                                         \n",
      "Time:                        05:46:36                                         \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": [
      "[ 4.81461806  0.55467676 -0.01470894]\n",
      "[0.47387364 0.0731597  0.0064735 ]\n",
      "[ 4.44689456  4.73087926  5.00996302  5.28414584  5.55342772  5.81780867\n",
      "  6.07728867  6.33186774  6.58154588  6.82632307  7.06619933  7.30117464\n",
      "  7.53124902  7.75642247  7.97669497  8.19206654  8.40253716  8.60810685\n",
      "  8.80877561  9.00454342  9.1954103   9.38137624  9.56244124  9.7386053\n",
      "  9.90986843 10.07623061 10.23769186 10.39425217 10.54591155 10.69266998\n",
      " 10.83452748 10.97148404 11.10353966 11.23069434 11.35294809 11.4703009\n",
      " 11.58275277 11.6903037  11.79295369 11.89070275 11.98355087 12.07149805\n",
      " 12.15454429 12.23268959 12.30593396 12.37427739 12.43771988 12.49626143\n",
      " 12.54990205 12.59864172]\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.74647706e+00  5.38741684e-01 -3.05990705e-03]\n",
      "[0.13056536 0.02015753 0.00178363]\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 0x7fa456bbc970>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAHSCAYAAADlm6P3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzde3zO9f/H8cdn581pzmHOcsphYwjRClEkoZKETkIHKiq/TjroW6kcQqQQFSqikzOTwxw2cwqrnM1ZNsNm23V9fn+8myHawbZrh+f9drtu43Nd1673Vcpz7+v1fr0s27YRERERESkI3Fy9ABERERGRnKLwKyIiIiIFhsKviIiIiBQYCr8iIiIiUmAo/IqIiIhIgaHwKyIiIiIFhkdOvlipUqXsKlWq5ORLioiIiEgBFBERcdK27dJXXs/R8FulShXCw8Nz8iVFREREpACyLGv/1a6r7EFERERECgyFXxEREREpMBR+RURERKTAyNGa36tJSkri0KFDJCQkuHop2cbHx4eAgAA8PT1dvRQRERGRAs3l4ffQoUMUKVKEKlWqYFmWq5eT5Wzb5tSpUxw6dIiqVau6ejkiIiIiBVqaZQ+WZU2xLOu4ZVnbr7j+jGVZUZZl/W5Z1geZXUBCQgIlS5bMl8EXwLIsSpYsma93tkVERETyivTU/E4DOlx6wbKs24B7gAa2bd8EfHg9i8ivwTdFfn9/IiIiInlFmuHXtu3fgL+vuDwAeM+27Qv/POZ4NqzNJYYPH86HH147y8+bN48dO3bk4IpEREREJKtktttDTaCVZVnrLctaaVlWk6xc1H+ZFxlNy/eWU/XlX2j53nLmRUbn1Eub11f4FREREcmzMht+PYDiwM3AUOBb6xqf7VuW1c+yrHDLssJPnDiRyZcz5kVGM2zuNqJj4rGB6Jh4hs3ddt0BeMSIEdSqVYu2bdsSFRUFwOTJk2nSpAkNGzakW7dunD9/nrVr1/Ljjz8ydOhQAgMD2b1791UfJyIiIiK5U2bD7yFgrm1sAJxAqas90Lbtz2zbDrZtO7h06X+NV86QkYuiiE9yXHYtPsnByEVRmf6eERERzJo1i8jISObOncvGjRsB6Nq1Kxs3bmTLli3UqVOHL774ghYtWtC5c2dGjhzJ5s2bqV69+lUfJyIiIiK5U2Zbnc0DbgdCLcuqCXgBJ7NsVddwOCY+Q9fTY9WqVdx77734+fkB0LlzZwC2b9/Oq6++SkxMDGfPnqV9+/ZXfX56HyciIiIirpeeVmczgTCglmVZhyzLegyYAlT7p/3ZLKCPbdt29i4Vyvv7Zuh6el2tYqNv376MGzeObdu28cYbb1yzVVl6HyciIiIirpeebg8P2rZdzrZtT9u2A2zb/sK27UTbtnvZtl3Ptu1Gtm0vz4nFDm1fC19P98uu+Xq6M7R9rUx/z9atW/PDDz8QHx9PXFwcP/30EwBxcXGUK1eOpKQkvv7664uPL1KkCHFxcRd/f63HiYiIiEju4/IJbxnRJagCYGp/D8fEU97fl6Hta128nhmNGjXigQceIDAwkMqVK9OqVSsA3n77bZo1a0blypWpX7/+xcDbo0cPnnjiCcaOHcv3339/zceJiIiISO5j5UC1wkXBwcF2eHj4Zdd27txJnTp1cmwNrlJQ3qeIiIhIbmBZVoRt28FXXs9stwcRERERkWtzOl29gqtS+BURERGR6xMbC6tXw8l/mn/98AMUKQLXOeMhOyj8ioiIiEjGHD0Kr78O99wDVauCvz+0agWLF5v7a9eGfv0gOdm167yKPHXgTURERERyyJkzsHWruW3ZYm4PPwxPPWVC7YgRUKsW3HyzCboNG0KzZua5derAqFGuXf81KPyKiIiIFGS2DXv3mnDr5wft20NiIpQqBUlJ5jHFi5tw6+9vfl+hApw9C77XN2vBFRR+RURERAqK5GTw+Cf+vfYarFhhdnZTWrXefrsJv15eMHYsBASY0BsQAJcOBbOsPBl8QeGXU6dO0aZNGwCOHj2Ku7s7pUuXBmDDhg14eXm5cnkiIiIimXP0KGzalFqysHkzeHrCtm3m/l27TIh9+GETcBs2hHr1Up/fv79r1p3NCnz4LVmyJJs3bwZg+PDhFC5cmCFDhly8Pzk5GQ+PAv+PSURERHKrpCQTZLdsgd9/h3ffNaH25Zfhyy/NY6pUMeG2cWNT5mBZ8N13Ll22qyjVXUXfvn0pUaIEkZGRNGrUiCJFilwWiuvVq8fPP/9MlSpV+Oqrrxg7diyJiYk0a9aMCRMm4O7unsYriIiIiGRCbKypy/X0hDlzTNDdvt3U6AJ4e8Ozz0K5cjB4MDz6KDRokFqrK7kr/A4ebHbks1JgIIwenfHn/fHHHyxduhR3d3eGDx9+1cfs3LmT2bNns2bNGjw9PRk4cCBff/01vXv3vr5Fi4iIiMTEwKpVJhyl3PbsgbAw02HB2xtKloRBg8yublAQ1KyZWtMbGOja9edSuSr85ib33Xdfmju4y5YtIyIigiZNmgAQHx9PmTJlcmJ5IiIikl84HBAVBZGR5talC9xyi6nN7dzZlCjUqGFKFh57zOzqAnTqZG6SIbkq/GZmhza7FCpU6OKvPTw8cF4yoi8hIQEA27bp06cP//vf/3J8fSIiIpIHxcebFmGlS8OpU3DXXSbkxseb+729oXp1E34bN4Y1a6B+fTMtTbKEJrylQ5UqVdi0aRMAmzZtYu/evQC0adOG77//nuPHjwPw999/s3//fpetU0RERHKZlSvh44+hVy+46SYoXBj+7//MfcWLQ4kSpqvC9OmpLccGDDD3+/lBixYKvlksV+385lbdunVj+vTpBAYG0qRJE2rWrAlA3bp1eeedd7jjjjtwOp14enoyfvx4Kleu7OIVi4iISI46ftyULGzalNppAeDxx+Gvv8xQiKAg6NbN9NIFcHODBQtct+ZsNC8ympGLojgcE095f1+Gtq9Fl6AKrl4WAJZt2zn2YsHBwXZ4ePhl13bu3EmdOnVybA2uUlDep4iISL5m23DiBKSc8XnlFbNre+hQ6mNatDDlCmDaj5Url/r4AmBeZDTD5m4jPslx8Zqvpzv/61o/RwOwZVkRtm0HX3ldO78iIiIi13L0qOmuEBGReouNNeUJXl5QrBjceis0amRugYGXtxVr2NB1a3eRkYuiLgu+APFJDkYuisoVu78KvyIiIiK2DQcOpAbcwYPNobSpU02Nrru7qdnt1MkcREtONuH3xRddvfJc53BMfIau5zSFXxERESlYbBucThNoN22CYcNM4D11ytzv7g7t25vw+9BDpka3QQPw9XXtuvOI8v6+RF8l6Jb3zx3//BR+RUREJH87dgzCw81t40bz9b33oG9fs3t77Bjcc4/Z0W3c+PKgW6mSueUj2X0YbWj7Wlet+R3avlaWvcb1UPgVERGR/OP0aRNuCxeG5s3N4bQbbjD3WRbUqQMdOkDVquZavXpZP142F7vyMFp0TDzD5m4DyLIAnPJ9/jdvNwf3ulO17oVc1e1B4VdERETytk8/hdWrza7un3+aa927w3ffmdKFCRNMvW6jRiYUF2A5cRht/374bUYFdkyuQNGisHJK6sTl3EBDLoBDhw5xzz33cOONN1K9enUGDRpEYmIioaGhdLrK2MCff/6ZoKAgGjZsSN26dZk0aZILVi0iIlKAOBywfbs5gDZgAPTrl3rflClmmES9evDuu7B0KXz2Wer9AwZA69YFPvhC9h5G27QJevY0A+rGjjWVJD//nLuCL2jnF9u26dq1KwMGDGD+/Pk4HA769evHK6+8QseOHf/1+KSkJPr168eGDRsICAjgwoUL7Nu3L+cXLiIikp8dOwZly5pfDxsG48aZscCQ2l4sxYoVCrbplNWH0WwbFi6EDz+E5cvNMLrBg2HQIKhY8XpXmz0KfPhdvnw5Pj4+PPLIIwC4u7szatQoqlatym233favx8fFxZGcnEzJkiUB8Pb2plat3FHALSIikiedP2+6LaxbB+vXm6/R0aZ+198fatY0h9OaNjW3G28009FSKPimW1YdRrtwAb75xoTeHTvMALuRI+GJJ8zPJgBhB8MI3RdKSJUQmldsnpVv47rkrvA7eHDWF50HBsLo0de8+/fff6dx48aXXStatCiVKlXir7/++tfjS5QoQefOnalcuTJt2rShU6dOPPjgg7i5qYJEREQkTU4n/PGHCbl33GGmn335JQwcaO6vVs3s6jZrZg6oATzyiLnJdUup681st4fTp2HiRFPWcPSomeExYwbcf79pnJFi6Z6l3PX1XSQ7k/Hx8GFZ72W5JgDnrvDrArZtY6X8x5WO6wCff/4527ZtY+nSpXz44YcsWbKEadOmZfNKRURE8qhjx2DSJFi71oTemBhz/euvTZHo3Xebz8ibNTMH1CRbdQmqkOHDbXv3mr3EL76Ac+dMG+QZM6BNm9SfUQAij0QyKWISUzdPJcmZBECiI5HQfaEKv1f1Hzu02eWmm25izpw5l107c+YMBw8epHr16td8Xv369alfvz4PP/wwVatWVfgVERGxbdNtISzMBN2QEHjwQUhKguHDzYG0++6Dm282t9q1zfMCAsxNcp2NG01pw/ffm0qTnj3hhRdMK+QU55PO8+3v3zIxfCLro9fj6+FLu2rtWLpnKcnOZLzcvQipEuKy93Cl3BV+XaBNmza8/PLLTJ8+nd69e+NwOHjhhRfo27cvfn5+/3r82bNnCQ8PJyQkBIDNmzdTuXLlHF61iIhILuB0mkTkdELXrqbdWMqUtGLFIOXvx4AAs9tbtKjr1irp5nTCr7+aGt7ffjP/2oYMgWefNbW9KXae2MmkiEl8ueVLYhJiqFOqDmM6jOHhBg9T3Le4an5zK8uy+OGHHxg4cCBvv/02TqeTu+66i3fffZewsDCWLVtGwCU/jc6cOZMPPviAJ598El9fXwoVKqRdXxERKRgOH4Y1a8xt9WooVcoc9XdzM/2sOneGFi3McIk6dS4/lKbgm+slJMBXX8FHH8GuXaYS5eOP4bHHzL++sINhTFm5FBub5XuXs3L/SjzdPOletzv9g/vTqlKry0pGm1dsnqtCb4oCH34BKlasyE8//fSv6yEhIcTH/7sdSKtWrXJiWSIiIq7jdJpCz5QSwF69TI0umNG/TZua3rkpvv8+59cowPWPKz51yswJ+eQTOH4cgoJMJ4fu3cHT0zxmzo459JjTg2RnMgDli5Tn/bbv0zewL2UKlcmOt5VtFH5FRETEbPtt2JC6s7t2LcTGmnKFIkXMrm6jRnDLLaaT0qVH+8Vlrmdc8e7dMGqUmRESHw933glDh5pSbcsCh9PBj1G/MDF8Igv+WnDxeW6WG081eYoXW76Ybe8rOyn8ioiIFESxsSbgNm0KJUuabgyDB5v76tSBbt2gZcvU0oX773fdWuWaMjOueP16U887d66pVunVC55/3pxHBDgSd4QvIr/gs4jPOHjmIOWLlOfRwEf5Zvs3JDmS8HL34rYq/56FkFco/IqIiBQEZ8/C4sWwapU5xbR5syltmDULHngA7r0XqlY1gfefQU6S+6V3XLHTCT/9ZDo3rF5tZoe89BI88wyUL29avC7bs5yJEROZt2seyc5k2lVrx5gOY7i71t14uHnweKPHc+UBtoxS+BUREcmPDh40IbdiRVObe+KE2c318TEH0l591Vy/+Wbz+EqVzE3ylLTGFcfHw/Tp5hDbn3+aBhyjR8NNHcLYeDyUbecbMTtsBxMjJvLHqT8o6VuSwc0G82Twk9QoUeOy75lbD7BllMKviIhIfmDb5qj+ihUQGmoOq4GZjNa6NVSpYsYGBwWpXjcfuda44n5N6/DmmzBuHJw8CY0bm03+bt1gw+G13D79di44Llx8TouKLXit9Wt0r9sdHw8fV7yVHKPwKyIikhcdOGBC7t9/m1pdy4IPPjDtyG69FQYNMl/r1zePtywzQU3ylSvHFRdPKkGZ3fXp/2FhEhKgUyfTo7d1azifdI6pW77hjdA3LgZfC4tnmj7DmDvHuPJt5Ci3tB+S/7m7uxMYGEi9evW4++67ifln7OK+ffuol1L9fYmUARhxcXEXrw0aNAjLsjh58mSOrVtERAqYRYvg0UehWjXz+XWfPqY/lW2b+5csMeUNc+ea8BsYCO7url2zZLsuQRX4oNXtBEZ1ZPOo5qz4sTC9esGOHabOt3TdHQxa+CzlPy5Pv5/74efph6ebJ+6WOz4ePvSo18PVbyFHKfwCvr6+bN68me3bt1OiRAnGjx+f5nNq1KjB/PnzAXA6naxYsYIKFTI2J1tEROSajh41zVafeMIUboI5rPbjj6Z0YexY2LrVFHKmDBa44YbLB0tIvuZwmJ9zWrQwHeh++w1eeQX274fxExPZ6phNyLQQbppwE5MiJnF3zbtZ/chq/nzmT1b2Xcnbt73Nst7L8kUdb0bkybKH7ByX17x5c7Zu3Zrm4x588EFmz55Nr169CA0NpWXLlixYsCDN54mIiFxTVBSMHw/LlpltOzDH8p9+Gho2hP/7P3jrLQXcAu78eZg2zUxf273bfBDw/MdhFG0YSoNytRm3I4LPN33OsXPHqOpflffbvs8jgY9QulDpi98jvxxey4xcFX4HLxzM5qOb//MxsRdi2XpsK07biZvlRoOyDSjmXeyajw+8IZDRHUan6/UdDgfLli3jscceS/OxN954I/Pnz+f06dPMnDmTXr16KfyKiEj6nTtnek4tXw4dO5qizJgY+OILaNXKlDS0aXN56YKfn2vXLC517Jj52WjCBDOVrWlTeO89KBu8hnZfteHCqtQ63k41OzEgeADta7THzdIPS5fKVeE3PWITYnHaTgCctpPYhNj/DL/pER8fT2BgIPv27aNx48a0a9cuXc/r2rUrs2bNYv369UyaNOm61iAiIgVAQoI5lLZ0qem8kJRk5seWK2fCb3AwnD6tbgz5yPWOHgbzgcBHH5mWZYmJZtjekCFQK+gk0zZPpd/s/112gG1oi6G83+797Hg7+UKuCr/p2aENOxhGm+ltSHQk4uXuxdddv77ubfuUmt/Y2Fg6derE+PHjefbZZ9N8Xo8ePWjUqBF9+vTBTR9BiYjIpWwb/vrLHEIDGDgQvL3NJLUbboDnnjM7uy1bQqFC5jHu7jqglo9cz+hh2zYfDHz4oSnz9vGBvn1h8GCbmMLr+TT8U2aPms0FxwUCbwjk3IlzOJwOvNy96FK7S3a/tTwtV4Xf9GhesTnLei/LlprfYsWKMXbsWO655x4GDBiQ5uMrVarEiBEjaNu2bZatQURE8rgFC2DOHLO7u3+/uXbLLSb8WpYJxL6+rl2j5IjMjB52OOCHH8z44Q0boFQpeOMN6PPEOZYdm0nPlROIPBpJEa8iPN7ocfoH96demXrZeh4qv8lz4Reyt0g7KCiIhg0bMmvWLFq1akVUVBQBAQEX7x81atRlj3/yySezZR0iIpIHXLgAa9bAypUmobi5meTy/fdw++1mfmy7dlC9eupzFHwLjPSOHgZTAj51qjnEtncv1KgBn34KN3eKYtr2TwmaMY3YC7HUL1OfTzt+ykP1H6KId5GLzy/IB9gyyrJTegNe6wGWNQXoBBy3bbveFfcNAUYCpW3bTrPBbXBwsB0eHn7ZtZ07d1KnTp2MrjvPKSjvU0Qk3ztyxPSXWrjQTFM7d87U7e7YYRJLTAwULgweeXJ/SbJQy/eWX3X0cAV/X9a8fDtgOtqNG2cOsZ0+bdqWdRywit99J/HHqSjCj4Tj6eZJ97rdGdhkIC0rtsRKaW0n/8myrAjbtoOvvJ6e/zKnAeOA6Vd8w4pAO+BAVixQREQkV4qLMyH3ppvMDm5EhGk9Vq2a6cjQoQOEhECRf3bh/P1dulzJPa41enho+1rs2GF2eWfMMOce770X+jx9hJ/iXueVyM8Bc3jtycZP8mbIm5QtXNZVbyPfSTP82rb9m2VZVa5y1yjgRWB+Fq9JRETEdWwbtmwx09QWLjRlDUlJpr/ua6+ZQ2p//ml2eUX+w5Wjh8sV8+WuUvX5/LXS/PKLqYB57HGblj1/48cjE+i2ei7JzuSLz3ez3KhcrLKCbxbL1GcylmV1BqJt296irXcREcnzYmMhOhrq1jV1vC1amKlqDRuargwdOphrYBKLgq+kU5egCnSqX4HvvzedG16LgNKl4ZU34yjaagbTd07g06W/U9ynOM82fZabA26mz7w+F7tahVQJcfVbyHcyHH4ty/IDXgHuSOfj+wH9wHRHuBrbtvN1/UpaddUiIpLDbNvU6P76q7mtXm2Cbni46Sk1bx7Uqwfly7t6pZKHxcWZmSWjR5vGHzVrwhvjfudIxQmM+X06Z387S+NyjZnSeQoP1HsAP08zxCSgaIA6N2SjNA+8AfxT9vCzbdv1LMuqDywDzv9zdwBwGGhq2/bR//o+VzvwtnfvXooUKULJkiXzZQC2bZtTp04RFxdH1apVXb0cEZGCKyHBBFuAxx6DKVPMrxs0gLvuMlPWbrnFdeuTfOPwYRg7FiZONB8q1O+0Cv+2E4jz2cXmo5vxdvfmgXoP8FSTp2hSvkm+zD+5wfUceLuMbdvbgDKXfON9QHB6uj1cTUBAAIcOHeLEiROZeXqe4OPjc1m7NBERySF798LPP8Mvv5h2ZH/9BRUqwAMPwM03w513gv7/LFlk2zYzie2bb0y/3rseiIbbX+Pn6KkQYw6wPRX8FMNvG04pv1KuXm6BlWb4tSxrJhAClLIs6xDwhm3bX2TVAjw9PbUjKiIiWWvDBnjkEVPaAFCrFgwYYModAO5IV+We5HNZMXrYtmH5cjOUYtEi8PWz6fj0Si7UH8+Cgz/giE7t9OBmuVGhaAUFXxdLT7eHB9O4v0qWrUZERCSjzpyBxYvhp59MqH3oIbObe8MN8MQT0KmTDqjJv1zP6GEwDUC+/dYcYtu8GUoHxHH3WzP4038C8/7+nRInSvDczc/RtEJTHWDLZdSBW0RE8h7bhvHj4ccfITTUJJHixc2hNTAH1ZYtc+kSJXfLzOhhMD9rTZ5sDrEdOgTVmu3gtpETCE+azk+JcTT2NgfYetTrga+nmeanA2y5i8KviIjkfk6n6cSwaxf07g2WBdOmwdmzMGgQ3H23aUWmqWqSThkZPQwm6I4ZA599Bmf8V1Gm00SqV9vF7vObOJTgxQM3mQNsTSs0/dcBNo0ezl30fwkREcmdEhPNZLX5883t8GEzNviBB8Db2xRaFi3q6lVKHlXe3/eqo4fL+/te9vstW8whtpkzwel3lCqPvMaZ4l9wHJsT5y36B/fnrZC3KF2odE4tXa6Tm6sXICIictGZMyb0Arz/vhku8eWX0Lw5TJ9umqV6e5v7FXzlOgxtXwtfT/fLrqWMHrZtU0Z+xx0QGGjz3bq1VH/xIdyHVGJP8c8Bc3DSzXKjUtFKCr55jHZ+RUTEtY4cMbW78+aZOt05c0wZQ69eEBRkxgn7+qb9fUQy4MrRw+X9fRl8ey3ObKtAwz6wbdd5irWcSfk3x3PYjuSId1EGBg6keUBzHpn/iA6w5WEKvyIi4honTpiQu369+X316vDss3Djjeb3Vauam0g26RJUgS5BFYiNNbW8z3WF6PO7KdnhUwrdN4VY52kqlq7H600m8lCDhyjsVRiASsUq6QBbHqbwKyIi2S9lnPCcOaZs4aWXoFQpc3vnHejSBerWNQfZRHLIgQPw0rgw5m5aQeI5H0p2W4ZVYgExlhtda3fl6aZP06pSKx1gy2cUfkVEJPts2QKzZ8PcuRAVZcLtvfea+yzLTF8TyWGRkaY/78yNi7B73A23JIEFtm8JXmvyGv0a96NC0YwNu5C8QwfeREQk6zidsG5d6iS1CRPggw/M0IkJEyA62uz+iuQw24YFC0wJeaM7t/Jd/JNYPe8GdxN83XDjuZuf483b3lTwzee08ysiItfH4YA1a+C770ywPXLEjBdu0gReew3efRdKlnT1KqWAunABvvkGPvw4iR3OeXi1Ggetf8Pd3Yf21duzZPcSkp3JeLl70aZqG1cvV3KAwq+IiGTetm3Qvr0JvL6+cOed0L071K5t7g8IcO36pMA6fRomTYJRk49xPGAynh0ngm805f2r8lSTkTwa9CglfEsQdjBMh9cKGIVfERFJH6cT1q6Fb7+FmjXh6adNZ4aQELjnHujY0QyhEHGhuRvC+HDuCiJWliax/EqsXt+CWxK3V2/P000ncmeNO3F3S+3vq8NrBY/Cr4iI/Ld168x4q++/N1PWfHzgqafMfT4+5jNlERcLD4cXx4eyIqA9+CRCB/B196Nf8EAGNhlIzZI1Xb1EySUUfkVE5HK2Ddu3Q/365vfvvQcLF5qShvvvh06doEgR165RBPNhxK+/wjufHGC941O4eSx4mAmBbrjx4i1DGR4y3LWLlFxH4VdEREzg3bYNZs0yt717Yd8+qFwZRo+GEiU0TlhyjYQEmDHD5u2vQjl4wzi4eR6WG7SocAvhR9dfPMDWvnp7Vy9VciGFXxGRgi4iAnr3NkMo3N2hbVt4/fXUDg1Vqrh0eSIpTp2CMRPOMXr5V8TVGQe3b6ewWwn6NxvK000HUNm/sg6wSZoUfkVECproaDN4onp1c1CtUiUoXRrGjzedGsqUcfUKRS6zeze8OXY3M3ePJ7n+FAiJpUahIIa1mcKD9Xrg6+l78bE6wCZpUfgVESkITp40B9ZmzoRVq0yZQ//+JvyWLg2hoa5eoci/fLZwDe8u/Jz9Z6MgYB1WsDsdKnbj1XbP0KJii3+NHRZJD4VfEZH8KjkZPP7533znzhAWBnXqwJtvwgMPmHZlIrmM0wmz551hyM9vcLjSGChugz88UOsRPu70DuWLlHf1EiWPU/gVEclPHA5YsQK++socg//zTyhWDN5/3xxYa9AAtFsmuVB8PIycGsWo1eOIqTINKp+9eJ+7mzsNA25U8JUsofArIpIfHDwIY8aYnrtHjpig2707nDtnwm+rVq5eochVHT/h5PkJC/j2wFiSKi3GutGTVv496NXqVgYveoZERyJe7l6EVAlx9VIln1D4FRHJq/btg6QkM2UtLg7GjoW77oKHHjK9eH190/wWIq4S8XsMz06dSljyeNsWoNkAACAASURBVOziu/G+oRyPVnmLEV37cUORssyLjKaG9T4HE8Op6BfMsZOVoKKrVy35gcKviEhecuYMfPcdTJ8Ov/0GDz5odnvr1oXjx8Hf39UrlAJgXmQ0IxdFcTgmnvL+vgxtX4suQRXSfF7YwTBGL/2WVVv3c8RvMRQ5R5kLLXjp5hE807Yrnu6eF7//sLnbiE+qRjGqceYMDJu7DSBdryPyXxR+RUTyihdegAkTTIf/mjXh7bfNLm8KBV/JAanB1AFAdEx8msE0McnB41M+ZsaRl8FyQjGobt/FhK5vc0f9Rv96/MhFURe/f4r4JAcjF0Up/Mp1U/gVEcmttm+Hb781Ayc8PEz/3UcfNQMpmjbVwTVxiYwE08N/m9KG+UfGkVxkz8Xr7m7uPHbbLVcNvgCHY+IzdF0kIxR+RURykxMnTC/eL7+ETZtM6O3cGYKD4aWXXL06kXQF09VRO3lu5ieEJ00Hr3MUTrqFTqUfZf7pEek6wFbe35foq7xOeX/Vscv1U/gVEckttmwxITc5GYKCYPRoU9OriWuSi1wrmJYr5s3E5T/zzpKxRPssAYc3leJ68sadz/BIhyAsC8IO3p6u0cND29e6rLQCwNfTnaHta2XLe5KCxbJtO8deLDg42A4PD8+x1xMRydW2bYMpU6BUKXjlFdOj9+23oVs3qF/f1asTuaqUmt8Yx3YS3Lbh6ahB0rlTxPv+xAW/PXCmAk0YyJjeT9C8Yenrep3MHKoTSWFZVoRt28FXXtfOr4hITjp92pQ1TJkCERHg6QmPP27uc3eH4cNdujyRtHQJqsCuvyP4v9WvYNuJJkn4gMfhlnTz+R9jhtxLhXKeWfI6CruSHRR+RUSym9MJbm7m188/D9OmQcOGZihFz55m51ckD3DaTuZtX8jHa4dgkwgWYFvc6jOQX0eNw8/P1SsUSZvCr4hIdtmzxwTdL7+E+fMhMNAcWnv2WVPTK5JHxF2IY+yqaXy8+hP+tv6EcyWxfD2w3Gy8Pb34X8+HFHwlz1D4FRHJShcumCEUX3wBoaGmHdkdd5h6XoDatV26PJGruVZ97e6/dzN84SfMjppCklscRDcjOPkbPnysG16VI9J1eE0kt9GBNxGRrBAbC8WKmTHD5cpB2bKpPXkraiar5F5XDq2wsXF6bqNQoVCiEpaA0x23XffT5YZnee+ZZtx4o4sXLJJOOvAmIpLV4uJg1iyYPNns+G7eDEWKmINsN96YWucrkoulDK2Id9vCGfefSHQcwOlxGP4ujd+OVxnQpD8vjS9P6cw3bpD8yuEwY9UPHTK36OjUX6f83t8fNm509Uovo/ArIpJR27aZw2qzZsG5c1C3LjzxhPmLwMMDaqkXqeQdB2L387fbDM57rTAH2NwsrPUvUtzuzqG5TfDVXImCybbh1CnYvx8OHEj9emm4PXw4taQrhacnVKgAAQHQuHGu/P+hwq+ISHr8/bcJtkWLmh3emTPhgQdM6L35Zo0allwlrR65tm2z+sBq3ls5hkPeP5igk8Jyo2jQUWoXiVPwzc+Sk02AvTTY7t+f+usDB+D8+cuf4+tryrgCAuC228zXlKCbcitVKtd/6qXwKyJyLbYNq1fDpEnw/fdmAMXQoXD//XDPPSYIi+QyV9bwRsfEM2zuNgDurF+KWdtn8f7KMeyMiYT44rDpBXyctUloORBIwsIdf7dATVPL62zbjEvfswf27jW3S3994MC/d23LlIFKleCmm+DOO6FyZfP7lK8lS+aLH/QVfkVErmTbMH48fPop7NhhQu7jj5u/DAC8vc1NJBdKqeG91Nmkkzzz8zDOLlhETNJxOF4Xz8iJPNq4F0M/LcS2M9G8vmAkB8+HU9EvmLfu7KoBE3lBYqIJsn/9ZW4p4Tbl65U7t2XLQtWq5tOqnj2hShUTbCtXNju6BWSrX+FXRARM4P3zT6hZ0+xs/PQTFC5sWpY98AAUKuTqFUoulBtH8B6OiQfggttOzrqFkmQd4YLbVnAmw66O+O8axPNd2jLwJ4uSJc1zqlOBLkHPuHDVck0JCSbM/vlnashNuR04YIbopChSxITbGjWgXTuoVs38vmpVE3T1/zFA4VdECrq4OPjmG5g4EbZuNfVuAQEwd67+opD/9F/lBa4MwOWKebHj7AxiPb8CywYb2NEd78hXGfdSQ3p9AT4+LlueXE1yMuzbB1FR8McfqV//+svU5V5ak12ihAm3LVqYVoo1aphb9epQunS+KEvIbgq/IlIwHToE77wDX38NZ8+accPjx5u2PKDgK2m6WnlBfJKDkYuiXBJ+T8efZvKmyUQxhlivwyb0AtjuFC1dlqmzStG1cY4vSy516pQJtlfedu82JQwpSpQwn0KFhJi2iZcG3BIlXLb8/ELhV0QKjvh4cwCkUiVzGvmbb6B7d3jySWjWTDsmkiEp5QXpvZ5ddp3cxdj1Y5ka+SUJjvNY+26DP5+F24aDeyJubh4M63oHXRurhjdH2Lb54XrHDti503zdsQN27TLhN4Wnpwm0tWrB3Xebrym3UqVct/4CQOFXRPK/PXvM4bUpUyAoCJYuhfLl4dixAnPAQ7JeeX9foq8SdMv7Z/+fKdu2Wbx7MaPXjWbh7oW4Ob1xbn4I362DeKJzAwbPgKMerTV+ODs5HKZU4cqQu3On+TQpRalSphd4t26XB9wqVUz7RMlx+qcuIvlXaCiMHAkLFpid3i5dYODA1PsVfOU6DG1f67KaXwBfT/csbxE2LzKa1xfM5eD5cMr71qdV3WRCo78k6tQuPBJugLVvUfLAkzzXrwxPTk39VLwqzRV6s4Jtw8GDZrjN9u2pt127zGG0FOXLm5D76KNQp475dZ06aDRe7qPwKyL5y8mTpkuDj48ZM7xpE7z2mhlGERDg6tVJPpJS15ud3R7mRUYzeO63HHB/GdsjkZgk2LEV3I83htUzuNF5P0Of96JnT3XfyxInT5pgmxJ0U77GxaU+pmJFqFcP2ra9POQWK+a6dUuGWPalJwizWXBwsB0eHp5jryciBciGDebA2uzZZihFnz6mxtfdHby8XL06kUxp8O4n7Er4mCS3fWb0sA2EP0nRXSOY9UlJOnRQqXqmJCaa8oQtW8xt61YTdI8dS31MiRJQv74Juilf69VTyM1DLMuKsG07+Mrr2vkVkbzL6YTp003oDQ83O76PPmoauIPKGiRPSnIkMWfnHEavG822pPXg9AM8TNsypyclKjSkaL113HlnR1cvNW84dcoE3M2bU8Pujh2QlGTu9/FJnWh2adi94Qb9ZJFPpRl+LcuaAnQCjtu2Xe+fayOBu4FEYDfwiG3bMdm5UBGRi86eNUHXsmDMGLOLM24cPPywRg5LnvV3/N98FvEZ4zaMIzouGt/zN8KKcbDjIfxaLMWjwUL8vGrjXbxSjhyqy3Ns2xxu3bQpNehu3gzR0amPKVfOtDXs0MF8DQw0rcR08KxASc+/7WnAOGD6JdeWAMNs2062LOt9YBjwUtYvT0TkH7YNq1fDJ5/A4sVmGEWxYrBokRq7S56288ROxq4fy5dbviQ+OZ5Cx9rA0kn4n7uT+7rHEXbrRhLdfYF7wZk9h+ryHIfDDIHYtCn1FhkJsbHmfnd3U4cbEmICbsOG5lamjEuXLblDmuHXtu3fLMuqcsW1xZf8dh3QPWuXJSLyj/h4mDnThN7Nm6F4cejXz0xEAv1lJnnS2gNrmbxpMjtP7mR99Ho88MZrVy9YPohqZesz5FXo0QO8vIoxL7JOrhuhnKOSkkyZwqVBd/NmOH/e3O/jY4Jtz57QqJFpZ1ivnk4AyjVlxT7/o8Dsa91pWVY/oB9ApUqVsuDlRKRAcDpNe7K//oLHHjM1eJ99Bg89BH5+rl6dSKbEJ8UzfOVwRq4ZiY0ZPeyx9XGSF73LbS1LM2QGtGt3+QcZXYIqFJyw63CYFmIbN5o6/o0bTfnChQvm/sKFTbh94gkTdBs1gtq1VbYgGXJdf1osy3oFSAa+vtZjbNv+DPgMTLeH63k9EcnnbBtWrTK7vIUKwbRpJvSGh5u/5FTaIHnU4bjDTNg4gYnhEzkVf8p0bbAApzv1AqoxbW1pGjZ09SpzWEqN7saNqWF306bUARFFikDjxvDMM+Zro0ZmIpqbm2vXLXlepsOvZVl9MAfh2tg52S9NRPKfCxdMi7LRo03dXokSlw+jaNzYdWsTuQ6bjmxi1LpRzN4+m2RnMiVP3APr20H7IVgeiXh7eTFhaAgNK7p6pTng6FFYv960JUwJu6dPm/u8vc2Obt++0KSJudWqpaAr2SJT4deyrA6YA2632rZ9PmuXJCIFzogR8Pbbplm8Shskj3M4HfwY9SOj14/mt/2/4W0VpkjUAP5e8Cw+harz4WBocFcQ4Sfz8ejhhATzg+z69bBunbnt32/uc3c3NbnduqUG3Xr1wNPTtWuWAiPNIReWZc0EQoBSwDHgDUx3B2/g1D8PW2fbdv+0XkxDLkQEMM3kR4+GBx80U5IOHTIHWq4sdhTJA8IOhhG6L5QmFZrw+/HfGbthLHtO78GfyiSveZazvz1GYJ1iDBkC99+fDzOebcPevakhd/16E3xT+uhWqgTNmpn+282amR1e/XArOSDTQy5s237wKpe/yJJViUjB4XTCr7+a0LtsmRlA0bixCb8BARo9LHlS2MEwbp9+OxeSL5gDbEC5pJZ4/fI+MVu70OEOD4b8Arffno9+rjt/3pQtrF1rbuvWmbHAYGr1g4Ph+edTw265cq5dr8gVdDxSRHJGu3awfDlUqADvvWdOa5co4epViWSKbduEHQrjyZ+fJCE54Z+LFoT35+TiCfTsCS9MN+c187zoaBNy16wxXyMjU1sN1q4Nd9+dGnRvukmdFyTX059QEckehw7B55/DK6+Yz3kffdQE3m7d8uHnvlJQJDuTmbNjDh+v+5gN0RvwcSuM5fQwu74OL3o1eJj3Jpmf8fKk5GRTlpQSdNeuTa3V9fU1AffFF6FFC2jeXD/ASp6k8CsiWSsyEj76yHRvsG1o0wZatTKH2ETyqJiEGCZHTOaTDZ9w8MxByrjXoEz4OI4v7kPZ+tto3D2U5+4NoW2tPHZ47fx5U6O7apW5rVuX2mqsQgVo2RKee858bdhQP7hKvqDwKyJZ49QpuO8+WLHCNKJ/+mkYNAiqVHH1yqSAmhcZfd2T0Xb/vZsx68cwJXIK55LOUYUQiv0ynuPhHWkU5MboaXDffc3x8Mgjoff0aTMmPCXsRkSYg2mWZcJt374m6LZoYQ6qieRDCr8iknkJCbB1KzRtaj7+9PaGDz4w5Q3+/q5enRRg8yKjGTZ3G/FJDgCiY+IZNncbwH8G4HmR0bz261z2JfyC5XmYs/Z2PNw8qHbuQfbNGsy+/UHcdRcMWQYhIXngEFt0dGrQXbUKtm83n8h4eZkWYy+8YD6ZadFC/81KgaHwKyIZd/IkTJgA48ebAHzokJnGtGCBq1cm+VhGdnJHLoq6GHxTxCc5GLko6prP+T5iH/3nv8opt2/AwwYnuO+9n6QfRrH3Qnl69TJNDG66KcvfWtY5cABCQ81t5UozQQ3MpzEtWphea61bm+Dr6+vKlYq4jMKviKTf/v2mU8O0aSb0duxodo4KF3b1yiSfy+hO7uGY+Kt+n6tdPx1/msmbJvPa0pEkup+ElPb3TnecB+pToXEiG7/JpR27Lg27oaGm3y6YT2JatzblR61bm5IGdWEQARR+RSQ9LlwwJQ0nT8LUqfDww2YLrE4dV69MCoiM7uSW9/cl+ipBt7x/6m7nX3//xZh1Y5i6eSrnks7hEdMMa+tz2M3fAfdEwJPSwYXwdP+dcuWqZPVbypz/Cru33gqDB5t6jHr1NBpY5BoUfkXk6pxO+OUXU8Nbo4YJvY0bw5EjULy4q1cn+UBGyhgyspMLMLR9rct2igF8Pd0ZckdNVu1fxcfrPmb+rvl4uHlQ19GT/bOeIyaqIV7lYvAtWw6qrsTXroe3e83LAnOOO3rU9MdetswcJlXYFbluCr8icrkLF+Drr2HkSNi1CypXhp49U+9X8JUskNEyhvTs5F4q5Xu8vmAuB8+HU8E3iPYNvXg34jXCD4fj71WCwLj/Y8e0p9hyqpyZ0zD0BNP3RpCQXAYc9wEmMA9tXytL3nO6xMSYWt1ly8xtxw5z3d/fhFyFXZHrpvArIpcbPtzU9QYGwjffmPZlqhWULJbRMoZr7eT+VzAtW+oAf9kvkeCZQGzyDH6PsKnoV5OGBz9ly/Te7HDzo3dvU8FTuzZAaepG1r/u9mhpuXTHu0ohN94pHUvL/VtM2I2IMJ+6+PqaLgx9+phe2YGB4O6epesQKaj0N5pIQRcdDWPGmBGlrVrBgAFw++3Qtm0e6OMkeVVGyxhSAmh6g+ne03sZsngI8cmp369UdG8Ofj6Vc8XdeO1leOopKFv236+T1WH3UvPCD/DNhB/o8lcEt+zfTKPonXg7knF6eODWrBm8+qoJu82amTp7EclyCr8iBdXOnaa04auvwOGAMmVM+K1USc3tJdtltIwB0hdM1x1ax0dhHzF351wsLCzcsZ2Awwvvrf0ZP86NPn2gUKHrfQcZsG8fLFkCixdz+y+L6BIfB8DvZarxZaO7WVu5IYfqBbP0jY45uCiRgkvhV6Qg6tcPJk82H63272/Gl1at6upVSQGSmTKGa3E4HczbNY+Pwj4i7FAYxbz8udkxlB1TnyHGPkDF1qEM6BDCi282z5nKgTNnzOG0xYtN6P3zT3M9IICFNZqxqkoQa6oE8rdfsYtPsa6+4S0i2UDhV6QgsG1zYvzWW039bv368Prr8MwzUKqUq1cnBVBGyxiuJu5CHFM3T2X0utHsjdlLxULVaB4zlojPHyHsXGG6dIEhQyrQokU2jx52OGDjRhN2Fy+GdevMtUKFzOG0p56CO+6A2rUZ8/6KDO94i0jWsmzbTvtRWSQ4ONgODw/PsdcTKfAcDpg71xxg27QJvv3WHGATyaPCDoYxP2o+h84c4uc/fib2Qiz1i7XEe9PzhH91Dz7e7vTtaz7MqFkzGxdy5AgsWmSmGi5ZAqdPmxr5xo1N0L3jDmje3IwRvsSVXS7A7Hj/r2v9bK01FimILMuKsG07+Mrr2vkVyY+Sk80Utg8+MB+53ngjfP45dO7s6pWJZNq0zdN4/MfHcdgmONbzu53ya95l25JmlCoFw9+AgQOhdOlsePGkJFi7FhYuNIF3yxZzvVw56NIF2rc3h0RLlvzPb5MVO94icn0UfkXyE6fT9P50c4MPP4QiReC77+Dee9UmSfIkp+1kwZ8L+CjsI1bsW3HJHe5s/7EtNx5rxqefmo5gvlldOXDggAm7CxfC0qUQF2fKhlq2NJ+mdOgADRpkuCtKdneUEJH/pvArkh+cPAljx8LMmRAZCYULm0b5ZcqoXZnkSQnJCczYMoOP133MrpO7KFcogMbJTxFhTwG3RNzwYkT/EIb2yMKf6xITYfVqM9lwwQLTEQVM95MHH4Q77zRtAIsWzaIXFBFXUPgVycsOH4aPPoKJE+H8ebPDGxtrwu+VDUxF8oAT507wafinjNswjhPnT1DHP4jWJ78i7L37OXrBk5BeD1G7fSi9W4fQvGIWHGQ7dswE3V9+MTW8cXGmTrd1a3j8cRN4a9fWD5Ei+YjCr0hetXev+UvZ4TC7UsOGQd26rl6VSKZEnYxi1LpRfLnlSxKSE2hWoiNVI15gw+wQ9vla9HvMTPatUaM5cB2h1+k0hz9/+cXcNm4018uXhx49oGNHM2SicOEseV8ikvso/IrkJbt2wYYN0Lu36cs7YgR07QrVqrl6ZSIZZts2n278lE82fsKuk7vwdvemReHeHP3hOdb/VofSpeGtt8zQwevqyBcXZzoy/PIL/PorHD1qdnKbNYO33zaBNzBQu7siBYTCr0heEBkJ774Lc+aAvz907w5+fjBkiKtXJpJhyc5k5uyYw/CVw9l1chcAbnhQZMEcVqztSK1a8Nln8PDD4OOTyRfZvx9++sncVqww3RqKFTOH1Dp2NF+zpS2EiOR2Cr8iudmuXfD886YmsWhRU9oweLAJviJ5TNyFOL6I/ILR60azP3Y/xb1Lgm2BZeN02BSuvpUpL3ekY0fTsCRDUsoZfvzR3FJakdWqBYMGQadO0KIFeHpm+fsSkbxF4Vckt7FtOHfO1By6u5u/0EeMMFOiihVL+/kiuUz0mWjGrh/LpIhJxF6IpVHJVlQ7OpaVC0pCr3ZYHol4eXrxzYgQmlfMwDeOjzeTC1N2eA8fNqm5ZUsYORLuvtuEXxGRSyj8iuQWtm1Om7/1lvk4dv58M5zi4EHtVkmetPXYVj4K+4hvtn2D03bSqkR3EsJfYP2cpvj5wdOPQ+vOy/jjQighVdLZveHECfj5Z7O7u3ix6XJSuLApY+jcGe66K81BEyJSsCn8iriabZtdq7ffhvBwqFjRHGizbXMAR8FX8oiwg2Gs2LcCP08/fv3zV5bsWUIhz0LcXmQgB74dzMqwqtxwgylff/JJKFECTOeGNELv3r0wbx788AOsWWNKHAICoG9fE3hDQsDbO9vfn4jkDwq/Iq72ySemJrFqVZg82QRfLy9Xr0rykXmR0dk+Tve3/b/Rbno7Ep2JAJTwKcldXv9j8+dPsnhPcerWhSlToGfPdORU24Zt20zY/eGH1Prd+vXh1VfhnnsgKEjdGUQkUxR+RXKaw2FGDpctC7fdBg89ZGp5e/bULq9kuXmR0Qybu434JAcA0THxDJu7DSBLAnBMQgyfRXzGO7+9czH4Yrtxdtkgfl3yMrfdBpM/MVUJ/3mIzeGAsDATdufNgz17TLht0cKM6u7SBapXv+71iogo/IrklORkM354xAiIijKDKW67zdQn9unj6tVJPjVyUdTF4JsiPsnByEVR1xV+98fsZ8z6MUzeNJmziWepUyyYqIRtOO1kcHpxa0Bb/hcOjRv/xze5cAGWLTOB98cf4fhx86lHmzbw8svmwNoNN2R6jSIiV6PwK5IT5syBl16C3buhQQOz89u1q6tXJXlURsoYDsfEZ+h6WjYd2cRHYR8xe/tsAFqX7MG5JS+wYX4QvjeGEdQ1lBe6hdC1yTXqeOPjzcHOOXNM4D1zBooUMQfVunQxX4sWzdTaRETSQ+FXJLskJZmvnp5w8qQZTjFvntnNynATUxEjo2UM5f19ib5K0C3v75vu17Rtm4V/LeTDsA9Zvnc5RbyK0K7IYPbOepYVGypRvjy8/z7069ccf/+rhN5z58xkte+/N1PWzp2D4sWhWzdza9tWB9ZEJMco/IpktaQkmD4d3nkHhg6FgQPh8cehXz8d0JHrltEyhqHta10WlgF8Pd0Z2j7t/rcr963kkw2fsOnIJvbG7KV84Qrc5fEBmyb1Y+H+YtSvD19+CT16XOWM5pkzpiXZ99/DwoVmx7d0aVPj3r276dCgGncRcQGFX5GskpRkksCIEbBvHzRpktpg393dpUuT/COjZQwpgTgj3R5iEmIYtmwYE8MnAmBh0SzhNbZ//Cq/nvGibVuYNgnuuCP157l5kdFMnBdOvfBQuuxZR/Pdm3BPSoRy5eCxx8wOb6tW+m9BRFxO4Vckq/ToAXPnmtA7fjzcead2eiXLZaaMoUtQhXQdbjsQe4DR60ZfPMSWwna4sXGtLz07e/HCCxAYeMmTYmOJGDuN4tO/5sc9kXg5kzlUtDQzGnWkWv8+tO6tMh8RyV0UfkUyKzHR7PTeey+UKgWDB5sdLoVeyUbXU8ZwLVuObmHk2pHM2j4LgFbFe3B0TTt21RgAbol4uHnx3cchdEnp3HDmjBnMMns2LFpE48REDhUtzdTgzvxaqyVbytUEy6LCUV/WKPiKSC6j8CuSUSmhd8QI2L/f1DI++6z5SFckm2WmjOFqbNtm6Z6ljFw7kiV7llDYqzC3F3qWfbMGExpeiQoV4KmgmpRsFEqH2iE0L14fZs0ygXfBAtOmLCAAnn6ae09UIPKfwHupzHaUEBHJTgq/IhkxZQq89ZYJvc2awcSJ0L69q1clBUx6yxguFXYwjNB9odxS6Rb2x+7nw7UfsuXYFsoWKkd7t/eInPgkSw7407AhzJgB998PXkkN4NdD8NxHpktDQgKULw/9+5sH3HwzuLlx/L3lcJ0dJUREcorCr0habDt1R+unn6BMGfj0UzOySuUNkgeEHQyjzfQ2JCQnAGBjc2OxurQ9O4W1H/dk0Rlv2reHIV9Am1aJWIsXQd+ZMH8+nD9vBk08/rgJvC1b/quGNztKMUREsovCr8i1OBypE9nmzTOdG7780jTkV+iVPOJI3BGGLRtGfHLqzmzlmN789dZU9rq50bMnvDDYQYPTK82f9/vnwOnTZvLgww+bg5xpdGnIqlIMEZGcoPArciWHw0xgGz7cjCFu0ABiY819mjwlecSuk7v4cO2HzNg6gyRHEm6447SBZC9OLOjP0BcsXmi1njLLZkLHb+HIEShc2ExZe/BB5peqywfL93B44TnKr1uZZpjNTCmGiIgrKPyKXMrhMLW8ERFw002mQf+996pVk+QJtm2z5uAaRq4dyY9RP+Lr4UsL78fZP/t59h47TrGGoTxbryKvdPoZ7zm9YOQeM52iY0d48EHz1c8vw1PkRETyEoVfEduGlSvNxCl3d/Mx79ChcN99Cr2SJzicDn6M+pEP1n7AukPrKOFTkhDeYPtnTxF6oDR33nSAubespOG2b7B+2W7+XLdtC6+9Zn64K1bssu+X0SlyIiJ5icKvFFy2Db/+Cq+/Dps2wW+/mdrGIUNcvTKRNIUdDGPJniWcSzzHD7t+4M+//6Ri4Wq0jhvHho8eITIukTfrf0/vhl9RfMtK+B1o3hw++cT8YFe27DW/d0anyImI5CUKv1Lw2DYsW2Z2vdatg2rVYNo0EwxE8oBFfy3i7pl3k+RMAqCyX22aHZhN5LSONHNbwqqA3jRKBVlxqgAAIABJREFU+Am3bYlQsya8/Tb07Gn+rKdDZqbIiYjkFQq/UvDEx5sg4O0Nn30GffuCp6erVyWSpgOxBxgVNorxG8dfDL7YbvBjK/pvWMb9Xv3xSzgN58vAwAHQqxc0bpzh7iRqXSYi+Vma4deyrClAJ+C4bdv1/rlWApgNVAH2Affbtn06+5Ypcp3Cw03QnTAB/Pxg0SKoUwd8fFy9MpE0bT22lZFrRzJz20wsy6KWRzt2nl+GbSXh7bSZuW8yN3v4YXXtagJvmzbgkfm9DbUuE5H8zLJt+78fYFmtgbPA9EvC7wfA37Ztv2dZ1stAcdu2X0rrxYKDg+3w8PAsWLZIOm3damp65883fUtXrjRdHERyOdu2Cd0XygdrP2DhXwsp5FmYRucfpt4X5eh79CccARtZUdXi1jJNaXnvM3DPPaZVmYiIAGBZVoRt28FXXk9za8C27d8sy6pyxeV7gJB/fv0lEAqkGX5FckxMjBnBOnu26c371lswaJD69EqGzYuMztAOaEYffyWH08HcnXP5YO0HhB8Op5R3Wboc6EX3Wae57/zneJHE2RqBFBrwMc17Pmimr4mISLpl9nOxsrZtHwGwbfuIZVllsnBNIpkXHw++vmYK29698H//By+8ACVKuHplkgdltN/t9fTHXbF3BaPXjSb8SDiH4w5TyaMSz4TdxovLthCQ/BWxvmU52/cZSjzXh8INGmTl2xQRKVCy/cCbZVn9gH4AlSpVyu6Xk4Lq8GF45x0zhnjXLrPDGxamPr1yXTLa7zYz/XFPnT/Fy0tf5vPIzwGwbIshK0vz3soDJNnH+LPuPRT6vz4Uf+CO66rjFRG5Htf7qVZuktn/kx6zLKvcP7u+5YDj13qgbdufAZ+BqfnN5OuJXN2pU/D++6Z3aXIyPPEEJP1zCl7BV65TRvvdZuT6/pj9jFo3ismbJnM+6TzYgAVuTpt4fPntgUk0HXkf9QOKZ3r9Iv/P3n3HN1luARz/vUkHBQplQ1mljDJlVaCAWFm1gCy3IIog18EQ2VtkShmCyhIEkSGCgAwBZZRl2WUP2aNlQ5ktHXnvH0eoaEGgI216vp9PLm2SJifXjJPnPc85SiUFR5v6+LTJ72LgHWD4X//+kmQRKfW4wsKgdGm4eVN2uH/22WP3MVXqcTxpv9vHuf7eC3sZsWkEc/bNwbCZvLHPSsCf8H5jg2grOFldeWP6j9QsrH2nlVKpg6NNffzPpTHDMOYAIYCPYRhnDcNogyS99QzDOALU++t3pZJfZCSsWSM/588vY4j37oUZMzTxVUmuW4APbs7WB857VL/bh12/a/0SrDu5jgbT6vLMxGdYtHM2nTbFcWCMMy/+/Cp/uq9idJUNDKo7hLXvrdHEVymVqjja1MfH6fbw5kMuqpPEsSj1cNHR8N13Mqnq0iU4fVp2uffta+/IlAN70n63/7x+vqwu1Cp7nOGrO7H1zp/kug2DN0OVbRX4JbYd41u+zgc9PWhxP5eukQKPSimlnoyjTX3U3RMqdYuLgzlzYMAAOH4catSA2bO1vZNKMU0r5n+iw3pNK+bHI9sRxv0+hu1nt/LFzpsUvQpBm92J2dmapZnaEdOjDP0/gtzaJ0cplQY42tRHTX5V6nbypIwfLlcOli2DwMAnHtWqVEq5fiWMXtNaMuFWMBhgAO+sK8nV4CFM9X6JTl8583srGTKolFJphaNNfdTkV6U+wcGwapW0LitaFLZsgYoVtXuDSp1Mk3MbV/DlLz2Z6LKHG65/v8jKWs9WjF3QnJdeAqv14TejlFKp2ZMeBUvNNJtQqceOHRAQAC+8AN9/D1evyvmVK2viq1Kfq1f588t+vN86J16/NWBkpj34nvOixOJeEOsGphVXZxd+HOZP06aa+CqlVGqhK7/K/sLD4ZNPYN48mcQ2ciR89JFMalMqNTFNWLeOrT8M44s7q1joY8O1oEGVU7XYtXQsIVEVaN0aBjd8iaOxwfh7+eNXUDs3KKVUaqLJr7Ifm01WdN3cZBpbv34yijhrVntHptQDQnYtZe3Sr3DbvIMlOa+wtghkiXWl9IGW7P91KIcz5qZ7Z/jwQ8iZE8Dvr5NSSqnURpNflfIuX4bhw+GPP2DjRsiWDY4dAxcXe0emVDybDVatYsMPQ6jrtZ5oK/AsZIl1J9eW3lxa/TExXu5MGg1vv60HKpRSKq3Q5FelnFu3YMwYCAqC27clY7h9G9zdNfFVqceFCzBtGpHfTWJ61pP0qWsQfe+d0mbhxvruPGf2ZMo8aNQoacrRF4WGOcwuaqWUA7DZ5DP7xg35uVAhOX/5crh4UQZORUZCVBQULChTVkEGT4WHx18eGQnz5987JJZqaPKrUsbevVC3rrxomjaFIUNkNLFSqYHNBmvXwqRJRCxfyPiKsYx91ZmLLuBxtyTEHgMjDqvhwsRedWj7YtLd9aLQsAf6Z4ZFRNJrwV4ATYCVUk8vLg6uXYObN6FIETlv+XLYvx+uXJHTtWuQIwdMnCiXN24sHZdu3oy/napVYfNm+blnT9iz58H7qV07Pvn94w/5nHdzk1OGDBAbm6wP82lo8quSj80GJ05IuzIfH6hfHz7+GKpVs3dkSolLl2D6dJg8mfALR/nSPwMTu1i4aUDOiDowqyd3L9Xi5Xab8fIP5uVKSb+BLWjl4QcaxwNExsQRtPKwJr9KqQdFRcH583DunPx7/rwksL17y+X9+8OPP8YntqYp03QuXJDLJ0+GRYvAyUmS3uzZoWTJ+Nv394dixSBLlvhTgQLxly9cKP/eS2zd3MD1b/0dN21K1oefVDT5VUnPNGHlSujVS15wR49KV/8ffrB3ZErd79jApEmwYAF/ukcT1DwPM/I7EUM0Wc+8Dsu6YzUrMLgDfPAB5MiRfBvYwhMYGfqo85VSDuz4cdi2Dc6cgbCw+AR32TL5HO3dW8oH/85iga5dpXwwVy6oVEkS2xw5pNzg76Mkp0yRVqLu7gkPjPr000fH5+2d+MeYCmjyq5LWli1yWCQ4WA6zBAXJt0Ol7C0iAmbMIOSn0QQbp8hjc2N510L87HIMKxG4Hnif6N+7kjePNyMHQ4sWKfPU9fRwIyyBRNfTQ3fQKeUQTFNOFgscOQK//irJ7d9PwcFylHThQklkATJlAk9PyJtX9sdkzAivvSYTT/PmjT/lyiUruQAdOsjpYXLkSPaHmxZo8quSzvbtUtKQKxeMGwf/+59uZFP2t2sXjB8Ps2bxR/Y71G5tcNcKEImzeQ7nLb2IXt+Rms/moesPMkE7JWeqdAvweaDmF8DN2Uq3AJ+UC0IplTimKSupYWHw88+ygnviRPy/S5dKSUFoqPS1z5BBNooVKCA1s/fedN56S0oECxaUtp//XJ2tVk1LB5OAJr8qcc6elUM0zZrJJLYpU+Sbqbu7vSNT6VlUlOwwHj8eQkKIy5iBRe9Vp1P+fdy9e1GuYzOI3dCF1/MMpEsw+PraJ9R7db3a7UGpNODGDamZPXJEEtt7ye3IkbLp68wZ6NRJVm29veVUt2586UHDhrLXIEeOhMsO8uWTk0pWhmmaKXZnvr6+5vbt21Ps/lQyunZNevWOGyffYM+elRe7UvZ0/LjU8k6dCleuEF2yODNbV2aE6w4ORxzB5W5+op0ugmHDyXBhboPVNK+iwyiUUn+x2aQ70ZEj8Oef8ae33oL27WVlt0ABmVdeqJCU9xUpAu+8A889B9HRcP261NomlNyqFGUYxg7TNP+1tKErv+rJREXB11/D0KFSQ9miBXz+uSa+yn7i4mDFClnlXb4cLBZuNWvI5Jc8GX15CWE3fyRDeEVYNZfs116mWfut5KwcTGApHT2sVLp18SIcOBB/8vGRWtm4ODkMdK89V/78UKKEdD0AqcE9dEhWdJ2d/3279zadqVRNk1/1ZI4cge7dISBAVn7Ll7d3RCq9unpVVnjHj4eTJyFfPi7378pXlWIYd+B7Ik5cwyXcH1Z/h3eGenTravDmm+DqqqOHlUoXTFM6JRw4IAs3DRvK+RUqwO7d8dfLkgVatZKfnZ2lrCF/fmn5lTnzg7dpGJIoqzRNyx7Uo5mm7EzdskVWeAEOHoRSpewbl0q/du+Gr76CWbMgKoqQlyqwqH4hTuZ2YcmRX4mMvYPTkabErutBHZ9qdO0q39X0CKRSDuz27fgjkMOHw5IlkvRGRMh5Pj6yYgvw5ZfyhlC6tJw8PVPVG4ROfEw6WvagntyWLdCjh/RELVFCfs6USRNflfJiYqQF0Ndfw4YN0li9VSv6+OZgaPgXcHkXXAaOBmJZNZLX65amy2KoWNHegSulktzp0xASIpPG9uyRL8QREXKyWGRDmbMzvPlmfIL794min3xiv9j/g058TBma/Kp/O3sWOneW3fJ58shh5bZtE65vUio5XbggE4kmTpR58d7eMGoU2xqUp8O6kWwJmyzXMwCbFTezDGO/y877DewatVIqKVy/Lq0KQ0MlyR0zRtp/TZkCgwZJb9uSJWWjWfny8iXZ1RVGjbJ35E9NJz6mDE1+Vbx7fQqdnGDjRhg4UKa9/LPmSanktmWLlDb89JN8oAUEYE6cyOqSLgzdFMTauV3grgccaA1l54A1BsNwImvhPEzfc4D3G+S19yNQSj2JK1ekc1CmTDIhtH17mQ56T5488nmUNSu0aQMvvyyJ799H6zoAnfiYMjT5VXDzpvQo3LxZds3nzSsbiBzsTUWlctHRkuyOGye9o93d4cMPsX34AYvMgwxZN4idO7dhuZ0PNgXhdPJNsla8hFNMGe6yhwy2crjaSumHhFKp3e3bMtFs58740+nTMHeu9InPnVtWclu3llG9FStK8ntP4cJyckA68TFlaPKbnsXEwLffygrvxYvypnP7tqz0auKrUsrFi9Kbd/x42Znt4wNff01MizeZdXIxg1c059j1QxjXisLGSbyQsxU9PsvAgO1rCL8eCZQkQ2zJ+zenHxJKpSIRETL9c+tW6bLQoIHU5DZqJEcaixeH6tVlpfdekX7FilJ2lw7pxMeUoclvenXoEDRuLK3Lnn9edsZWqWLvqFR6smsXjB0Lc+bA3bvw4ouEtKnPbzmvExF1lDmTKnAh6gycL4/ljzm8Wf4Vuk11ut9d73ZO/ZBQKlW5Vzpns8mq7ebNMiDini5dJPktXFg2rpYvr9NA/0EnPqYMTX7Tm2vXIFs2efPx9obRo6X3YSpq86IcWFycfNH68kvpIpIxI7z3HnTsyEqnUzSa04hY21/N5c89g9vmiXxcP5BOywwKFHjwpvRDQik7stkksd28WVZ1t26VsbxLlkjHhbNnpTPQO+/As8/K4Ihs2eRvDQNq1rRv/KlY04r59X0smWnym14cPAg9e8qO2UOHpFXUihX2jkqlF9evy0CKr7+GEydkLGhQELRpw3nnuwRtHMO4LV8Sy1+Jr81Cg8JvMGdEg/uDlRKiHxJKpZCbN6VvbtWq8nvjxrBsmfycJYskt9Wrx19/9eqUj1Gpx6TJr6MLD4fPPpPEI3NmSYBTcLCJSueOHJENbNOmST15zZqS9DZpwombZxi0pg8/7PuOWFsMnKiN4bUBwxqLq4sLfVv6PzLxVUolo7NnYe1a+OMP6am7d6+s6F6/Hn/Epnlz8POTOn2Lxd4RK/XYNPl1ZPv3y+Gm2FiZWd63L+TMae+olKMzTVi/XkpqliyR1nlvvgmdOkGlSuy/uJ++s1vzy/E5mDYL7HqX5yzdGdChGG7FQ1h3Khh/L3/8CuoIYqVSRHS0bErbuBHefVe6LcydC127Sk1utWry+eHnJ69nkMRXqTRKxxs7mpgY+YZeqZIkIX37yjf0okXtHZlydDExMG+eJL07dkCOHPDRR3LKm5ctZ7fQfckw1l/8BaIzYtn5AS/n/5T+n+SnbFl7B69UOnP2rAyP2bBB6nWjouT8X3+FwEA4dw4uX5bJaFarfWNV6ik9bLyxJr+OwjRl/GvPntIu6uRJyJ7d3lGp9ODaNWmZN24chIXJIdBPP4WWLfnj8i6m7JzKuj93czxyB0Rmw3VXRz6o2IHuHXLg6Wnv4JVKBy5ckCR3wwaoUye+00+pUtJW7LnnpCSpZk1Z9VXKQTws+dWyB0fwxx/QrZv8W6oUzJ4dv6tWqeRy7Ji0KvvuO6nnrV1b+vUGBmIzYEjwCAas742JCSa4H2xPH7+hfDTPXbsbKZXcYmPhf/+ThPfIETnPzU2GGDVuDMWKyRdXfTGqdEiT37Tuzz+hRg15Q5s8WXorOul/VpVMTFO+ZI0aBYsWxdfzdu4MFSoQa4vl25DZDFg1jEscABMwwGJY6fGRJz2e1w9apZLcmTMyMS04WEYEf/ONvDYPH5YFkXbtZFW3UiVwcZG/MQxNfFW6pVlSWnTpkrSReeMNKFFCRsIGBko3B6WSQ1wcLFggY7C3bpUjC716wccfg6cnUbFRjPh1IiNDRnDT6QRcLEO+m304X2QkJjFgOGGNLWPvR8Gi0DDtC6wcx/DhUnJ0/Lj8ni0bNGkSf/nGjfaJy470Na4ehya/acmdOzIcYPhw2Z1bu7bUZ736qr0jU47qzh1pUzZ6tHzAFismq0rvvAOZMnHz7k16zxnJt/tHcdf5PMaFKtRz/ZKGLzzLxN27MWM8iLLsJYOtHN8Hu1Aye5jdPogWhYY9MBEuLCKSXgv2AuiHo0rdzp2DNWvkFBICO3fKCq/NBs88Ax07gr8/lCuXrluO6WtcPS5NftOCuDiYMUM6N4SHyzf74cN1Y4JKPpcuyUCKb76BK1ek1dFf/XmxWrl48wodJgSxIGwcsc7XcDpXhxaeswga+AL58hnUGL6GyJg4XCmFq60UAJG2OIJWHn7oh1Byr9gErTz8wChkgMiYR8eklF2tXCmbRw8ckN+zZZNx9Fevgqcn9O5t3/hSGX2Nq8elyW9aEB4OH34oc9DnzIFatewdkXJUR47IKu/06dL6qHFj2UxZowYYBvP2LKbXkiCORW0HpygynmtCp7K9GDiuKpkyxd9MeERkgjf/sPNTYsXmSWNSKsVERUkt/erVchowQErZsmWDAgWk926dOlChQrpe2f0v+hpXj0uT39Rqzx6p5R08GAoWlDrLcuVkk4JSSW3zZlnZXbgQnJ2hVSvo0gVKlgRg29ETvD3zUw6zSK5vsdIh/w+M6dsywRagnh5uhCXwgePp4Zbg3afEis2TxqRUsrtyRfZubNwoCbDVClWqxF9epYqs/qrHoq9x9bj0K2Rqc/asdGyoUAEmTJDfQeq6NPFVSWRRaBg1h67i/Zf7s7tIOZnctHatHEY9dUo20ZQsybKtByjRoxVVfijOYdsSwABDPqPz+Zx5aO/7bgE+uDk/eKGbs5VuAT4JXj8lVmyeNCalktTZs1I//+ab0o8dZGU3JgY++ECmIV69KivAgYH2jTWZLAoNo8bwNRTpuYwaw9ewKDQsSW9fX+PqcenKb2px+zYMGyaHnOPiZNWtd2/t16uS3C9bT7Bl8NdM3/QTxa6e5WyW3Ayp/wHP9O/MSzVKYJrw3a876Pf7UM5lXQhObpS705EPXvSn6+Y3iI6LxsXqgr+X/0Pv495q7ePW8KbEis2TxqRUkhg2DGbOjK/bzZdPNoyClDAEB9sttJSUEqVN+hpXj0snvKUWN29C8eJS1zVkCHh52Tsi5Whu3YJvv+XCwGHkuX6JA7mLMLHqKywrWZM4ixXPLG4E5ndm1LYhXM+1EuNuVmq6dGBi646ULpwLgJAzIQSfDMbfyx+/gn5JFto/PxhBVmyGNS+nH1wqbTBNKVf77TfYtg3mzpWjdR06SD/2+vXlVLZsujyKV2P4mgS/4Ob3cGNTz9p2iEilBzrhLbUxTVi8WA6DzZ8vzcYPHtSVXpX0Ll2Cr76S7g3XrnG8UDm61W/P+iKVwDCIi7YQceYUp7LNICQmBIt7LhpnGsrEjh+RL1vW+zcj3RgiCY94hqUekXQLSLq2Zbpio9KsbdvktbVypYwRBklwL16EPHnktad0M5pKVTT5tYdt26BrV1i/Hnx8ZDpPkSKa+KqkdeqUTGKbMgUiI6FZM+jRg65rbxMWEcntqCNERK4h1uUwlP0T41Z+Wucby5et2uKeIeMDN5VShyw12VWpms0GoaGwfDk0bw6lS0sP3qVLISBATvXqSRsy9QDdjKZSE01+U9L169KybM4c6dE7fjy0bSu765X6S6L73e7dCyNGyPPMYoGWLaVdWSnpt/vqmTB6zA8iuuQ4yGaCaeB+7W2mvj2QVysXSfAmtX+mSreiouCXXyThXbEifnU3d25Jfhs0kFXeh+3+VIBsRkuotEk3oyl70OQ3JdhskoRkzixTsvr2he7dda66+pdErbBu2iTDT5YuhUyZoFMn6NwZChTANGHV2rt0mfE9e7N8AaWOw71yf8OgSfX8D018QQ9ZqnTENGHXLtmHUauW/P7uu5Axo9TsNmggK7z3hgw56cfo49DSJpWa6Ks2OUVHw8SJssK7eTN4eEgbG21Srh7iiVdYTVM22AwZAhs2QM6c8Pnn8PHHkD07sbEwe84d+i6czJmCI8ErjPw8y3vV2jFy+8D7nRs+8mv8yLj0kKVyaHfuyHCJpUvlFB4Ovr5SoubmJuOES5TQ1d1E0tImlVpo8pscTBMWLJBejkePSgeH69cl+dXEVz3CY6+w2mywaBEMHQo7dsgUqLFjpYwmY0Zu3oTxX95g+OpviCg1Bspcorjz84xuNo2GJetiGAaBpWs9ducGPWSpHM7585A3r/zcooW8njJnllXdRo0e7LX7V8mQUsoxJCr5NQyjM9AWOYC6F2htmmZUUgSWZt26BS++KIegS5eGX3+V39Nhaxv15P5zhTU2Vmp5hw2T7iDFismGtrffBhcXwsNhxKArTNw1lrsVxoHvdSq5v8iY5n2o5VXzgdv0K+j32O3K9JClSvPi4mQl997q7u7dEBYmm9O6dpWjJbVqgYuLvSNVSiWzp05+DcPID3QESpumGWkYxk/AG8D0JIotbbl+HbJmlZUDHx9pYt66tdaDqSfysBXWHv6FpYRmxAg4cUJGXc+ZA6+8Ak5O7NsHncYvYe3dEZj5tkO1KJ7P3YxRTfpQ2bNyksSmhyxVmhUcDK+9Jm3/rFaoWVPGebu6yuU1atg1PKVUykpsZuYEuBmGEQNkBMITH1Iac+2a1FtOmiSbJIoWhalT7R2VSqP+ucJa1A2+vP4HZZu/Jy2VqlaV8oaGDTENC2vXwudjT7Euw6dQagEAFouVGU1n0uKZFvZ8KErZx/nzMir4l1/g1VdlIaJ4cWlB9tJLUtagbSWVSteeOvk1TTPMMIyRwGkgEvjNNM3fkiyy1C46GiZMkM1F167JbuBMmewdlXIATSvmp6lXRmmOP3YsXL0KtWvLiNQXXiAm1mDejzB4whEOZh8O5WdgWG2YGPDX/56+ftreD0OplGOa8MUXUre7ZYuc5+UFTZvKz/nzw6xZdgtPOZZEt6NUdpeYsodsQBOgCBABzDMMo6VpmjP/cb12QDuAQoUKJSLUVCQ6GsqXh0OHZDUhKEh+VyqxLl+GMWNkYtSNG7JS1bs3VKvGzZvw7RgY8f0+LhQfCrXn4mxxoW3lD6lftDZvLXjrfvcGfy9/ez8SpZJPXJx0zjlyBN57T/ZUzJ8v/w4aBE2apNsxwip5pcTAH5X8ElP2UBc4YZrmJQDDMBYA1YEHkl/TNCcDkwF8fX3Nf95ImnL4sNTzurhAmzby5hoQoG+wKvEuXJBpbOPHS9ulV16BPn2gfHnCwmBcDxi/cAe3Kg2B5gtxs2SmfdWudKn+KXky52FRaBjFjC84E72dghl9uXC5EBS094NSKglFRcHvv8PChVLWcPmy7LNo2VLekzdskLZkSiUjHfjjGBKT/J4GqhmGkREpe6gDbE+SqFKbEyegVy+YO1dGEj/3nOwOViqxwsPlyMGkSXD3LrzxhiS9pUuzZw+Megdmrd9EXM3B0GIF7k4edK7en45VO5IjYw7g7ysR3mTFmxs30JUI5Rhu3IAMGSS5/eIL+OwzSXgbNpSShoCA+O4MmviqFKADfxxDYmp+txiGMR/YCcQCofy1wuswrl2TPqrjxskO4X79oEIFe0elHMHp0/JhPnWqtC97+23o1QuzeAlWr4Y+/f9ga/RUjPyhmO+Gkt01J91qDuOjZz8ii2uWB25KVyKUQ7l0CRYvll7pq1bBzz9L391WraBaNXjhBW1HpuxGB/44hkR1ezBNcwAwIIliSV3i4qBKFTh2TDazDRokmyaUekJ/3xxR2RbByKO/4rV0nlz47rvQsycxBb2ZOxeCXjPZk3EU1O8Bhg0T6FS1E0NqDyGTS8IbKnUlQjmEy5el3GfDBhni4uUF7dtLL2uAIkXkpJQd6cAfx6BNaP/u3qjYevVkpfeLL6R1mW5mU0/pXklCngunGREyj2b712CzWDjevAXeIz/netZCfPstjPnSRniWRWSoPxiyhd7/e6thJU+mPA9NfEFXIlQadeyYbFKzWKBbN8ieHZydpeyneXN539X9FCqV0YE/jsEwzZTbg+br62tu355Ky4K3bpU63g0bYN48WYFQKpFe7zqD11Z+T9P9wcRYnZhVIZBJVZpjcy/K8zE1mDwljluFfiJT4BBuZ9pP8ezFea3Ma4wOGX2/c8PqVqsfOYntn7uPQVYihjUvp2/IKnU5ehR++kmS3tC/vuS99JKUOSilVBIzDGOHaZq+/zxfV35PnpRWUnPmQO7cMkXrXm9IpZ7WkSMweDCzf5hJtNWZqb5NmFy1OeG3CnMj2Jvbh3Oy85npZG4/FFyO4JWrDH2em81rZV7DarHSsHhDgk8G4+/l/58jiHUlQqVqhw9DiRKyijtiBHz7rdTujhoFL78MhQvbO0KlVDqTvld+TROeeUZWI7p0gR49wN3d3lGptOyvpJeZM8HVldm+jRj1TGPOXizBja3eRJ11B98DUs5JAAAgAElEQVQpWP2HE+d2lop5K9K3Vl+almyKxbDYO3qlEs80Yf9+Wd2dP19+3rYNfH1lscFqhYLah08plfx05feemBiYNg1atJCJbFOnQr58+masEucfSS+dO3O3Yzc2T3Nh/zdO3M20HWq1goKbweUaRbNVYkzgJAKLBWJoXaNyFIcOQbNm8q9hQK1aMqnQy0suv/evUkrZUfpJfk1TGqN37y6H4ZycZDJQlSr2jkylZQkkvRHvd2PiwjyMqwbnrtwk88vduVtiIhgABm3Lfs7k5n016VVp3+HD0v/c0xPatpUSBm9v6NhRkuC8ee0doVJK/Uv6SH537JDNbMHBMqFt8WLpG6nU00og6T39ejdGz8rDlMpwO/Y6Rd/6CveiY7gZe/X+n1kNC965nTTxVWnX8eOS8M6dC7t3ywpv69aS/Lq5wbJl9o5QKaUeKX0kv926wb598M038P770k5Hqadx/Lj0fJ4xQ5LeTz5hV71uDJ+el3nVwMh4ldJtvuRE7nEci71OI+9GNCnRhI4rOt7v3uDv5W/vR6HUkzl/Pn4V95NP5ChatWowZgy8+qr2QFdKpSnpY8PbyZOQLZuMxVTqaZw+DYMHY5s2jRgszCwfyPSiHYg4XpN9213JlPsiZdqOZn/Gb7gde4vmpZrT97m+VMxXEYCQMyGP3b1BqVTh3DlpS/bjj7Bli7yPFiokG9gyZ9YuDUqpVC99b3jTTRbqaYWHy4jrb78lzoQ5zwQyLEtHjuytSswOd6x5j1O61yhOZJrGttgoXi/5On2e60PZ3GUfuBm/gn6a9Kq04cABqdldu1YmrZUvD0OGQMaMcnmZMvaNTymlEil9JL9KPamLF2H4cJgwAWJjiXqrNS+ceJutO5/Fli0US41+OHvtJSbbBg4YcbQq05LeNXvjk1NHXKo05s4dWLpUJqzVrQs5c0JYmExae/NNKFXK3hEqpVSS0uRXqb+7cgWCgqQ9U1QUN5u1YnSmfoyY782dO+D0/FRs/u2wYcMGuMVVIXvs+3zftK29I1fq8cXEwO+/y3CfRYvg1i2Zalm3rgz7OXjQ3hEqpVSy0eRXKYCICBg9Wjbw3L7N5fpvMsgYwNcLS2C1wkvvHGFtrs5cc/kVMKVtmWnB1SxJ4aze9o5eqSfTsKEkvx4e8MYb8NZb0pNXKaXSAU1+Vfp28yaMHSujViMiCK/+Cj2jPuOHlWXImhXe63GQq2WGsOjYHJwMZ9xja3LL2IJpxmLgRFajAt0CtNRBpWIHDkhLvoULZeNalizSsaF9ewgIkK4lSimVjmjyq9KnqCiYOFE2s126xMnyjfkk60B++aMChQtD91F7OZpvMFP/nIfbKTc+rfYpXap3YfOROPovX8CZO9spmNGXzwOb07SitnlSqcyVK9KOb+ZM2LkTLBaoXx8uX5bkt0EDe0eolFJ2kz5anSl1T0wMTJ8On38OZ89y3LsuH18bzIprValcGV7pEMoW10EsOryQzC6Z6VClA52rdSZXplz2jlypR7t1C65fl567Bw5IVwZfX2jZEl5/XaetKaXSnfTd6kylG4tCwwhaeZjwiEg8PdzoFuAjK7M2m0yk6t8fjh7leB4/PnKZwcrjL1D9tRAa1/mIK9Y99Dq5iayuWelfqz+dqnUiu1t2ez8kpR4uNlZqd2fOlI1rTZvCrFlQujQcPQpFi9o7QqWUSnU0+VUOY1FoGL0W7CUyJg6AsIhIev28h7zrfqPad2Ng715OZn2GDizht6sNebuVwWtVR/JTeHc4J0dA6hV+i5/e+AaPDB72fChK/bchQ2DcOGnLly0bvP02vPNO/OWa+CqlVIIs9g5AqaQStPLw/cQXwO/UbmZ915lqnd/jzJEoXudHKhuhPNO7EXO3bGRbxef4Kbwb3Cv9MS3sOu5C8MHbdnoESj3CuXOS7Mb99RyPjoaaNWUj27lzUsPup4NUlFLqv+jKr3IY4RGRAFQMO0TXdT9Q48xuThsFaMMU1uV5h46fWnmnbjBB2z5n6OJgnPHAPbYRt6y/3e/eYI0tQ9DKw7qJTaUOkZHwyy/w/ffw229SvuPrC9Wrw8CB9o5OKaXSJE1+lcOofvcCrX6ZQsCJEC6Qm46MZXqe18hb+yJfD1jLkE2fs3HeRvJmzsuYgDGMWVQQgwxkinueKMteMtjK4WordT+JVsquDh6EatXgxg0oWBB69ZLSBh9traeUUomhya9K+06f5nrnAcxYMINbZKYvg5jk/RaWaufJ6j0FS65FBM7ZQX73/HwV+BVtKrbBzdmNecFrCIuIxNVWCldb/AhXTw83Oz4YlW6dOCHtydzd4dNPoUQJaNUKmjUDf39pV6aUUirRNPlVaZZ56TLh7YeSe/43uNoMxlk7szGgI3tKryAyw/vYnI4SyWkK2QoxoeEEWldojatTfEP/bgE+D2yQA3BzturQCpVybt+G+fNh2jRYtw4MQ6atAVitMmZbKaVUktLkV6U5cddvceh/Y/CaF0Re223muLzLhQ8G0KJXQaKPjmDh6l6YmBgY9KzRk4EvDMTF6vKv27lX15tgazSlkotpSpIL0KGDJL7Fikn3hrfflhIHpZRSyUaTX5Vm3L4WzY7/Tab0gkGUibvIykzNuNp5ME16lOT3M4t4cVFjdl/Yff/6FsNCFtcsCSa+9zStmF+TXZUywsKkrGH6dFntLVdOyhtat5auDfcSYqWUUslKi8hUqnfhnI15TWdxKWdJas3rwOmMpVg7LITaEfNxemM/NX6owMs/vcydmDv0q9UPNyc3rIYVF6sL/l7+9g5fpWfR0TBvHgQGQqFC0Ls35MkDd+7I5WXLwnPPaeKrlFIpSFd+Var11fzz/NF7FT2PBPEqezic+RkODFxO+U51mXdwPhUmt+HApQP45PBhZrOZvF72dZwsTgQWCyT4ZDD+Xv74FdS+p8oOrl+HrFnh7l14913Inl26Nbz7rpQ4KKWUshtNflWqYpqwfj1M/WQj7+7qTwfWctylEB/V6Mva6lUILH6UFRM/5eDlg5TKWYrZzWfzWpnXsFqs92/Dr6CfJr0q5V2/DrNnw9SpMohi507p3LBlC5QqJRvYlFJK2Z0mvypViI2FBQtg9uDjvL63DzP4kctO2Wj9YnMWl8tArOUod6zTGLMzjLK5yzL3lbm8UvoVLIZW7ig7Cw2FL7+U8obISHjmGWjbVgZSWK1S2qCUUirV0ORX2dWtW/Ddd/D9yEu0OjOYn5iA4ezEuEpvMKZmKU5mHgTEggFWW15y3e3F7g8Ga9Kr7OvCBciQQUob9uyREcOtWknSW7my1vAqpVQqphmEsotz52TvT4kCdzjfaSjrw4vS0fI1zm3fheMHGfdSQU5nGgmGJL6YBpnj6lE8S11NfJV9xMbCsmXQvDkUKCDlDQBvvCFP6IkTZfSwJr5KKZWqaRahUtT+/fDee1C0cCwXh01l793iDKUPmRrVJnZPKNM+qkbJBbU5FhuE1cgCphOYFgyc8TAq6gAKlfJMEwYMAC8vaNQINm6ETz6RnwFcXSFTJruGqJRS6vFp2YNKdqYJwcEwciT8+qvJyy5LOZa5J/muHYAK1Yj5YhYzMh9jyKqmnIg4QeV8lVn8xmJib1dkwIqFnLmznYIZffk8sLn25FUpIzYWtm+HatVkJXfzZqnlHTdOkl6Xh/eOVkoplbpp8quSTUyM9PIfOVI2vtfPto2TXl0pfHI9eJUgevKPzChygyEbW3My4iS+nr58FfgVDYo3wPjr0HGzSh3s/ChUunLypJQzTJ0KFy/C6dPg6SnlDk76dqmUUo5A381Vkrt5E6ZMkQ3wp09Dbe+THKnSm2Jb5xBSIBs/BDXiZrmSzD3Yg1P7T1ElfxW+afANgcUC7ye9SqWoAwegSxdYuVJ+b9AA2rWD3Lnld018lVLKYeg7ukoyYWFyVHjSJGl52qB6BGt8h+K9dCzGOSvr+75NXZcfibm9FDYvpXTO0ixvsZyAogGa9KqUd/KktBspWxayZJEEuF8/aNNGprEppZRySJr8qkTbuxdGjZL+/nFx8HqzaIYVnkjh7z+HkKtEt2rBtJZl6LJlGDGxMX/9lUHFnC/xYrEX7Rq7SmfudWyYOFFWeevVk38LFIATJ8Cie4CVUsrR6Tu9eiqmCatWwYsvyj6g+fPhww9Mwr9ewOzdZSg8uhPRFZ9h0oLeFCu/jg829SImJlt89wbTmY378rAoNMzeD0WlF5MnQ5Ei0LSpfGPr1w++/Tb+ck18lVIqXdCVX/VEYmJg7lzoNyiak3+6YMkURaH6YUxscpDA2cPhq03cLVuSadM/ZujVxZzZvRa/An64R7bn1o3SRFsOEWXZSwZbObCVIGjlYe3goJKHacLatVC9ugykuHkTSpeGr76Sjg1ax6uUUumSYZpmit2Zr6+vuX379hS7P5V0rl+XRbKxY+HsWXDJeQv3Z49RMv9OemycxkuHNnA9V07m9KrPUDZw5sYZ/Ar4MdB/IHW96+Ld61cSeqYZwInhDVP64ShHdu0aTJ8upQ1//gk//AAtW0oyrLXlSimVbhiGscM0Td9/nq9LH+qRzpyRhHfyZFk4e+EFcK+zG4vHYdqHzKX1b4uJdLLw6muVWVL6FHdvzKZ6wepMbTyVut51729k8/RwIywi8l+37+nhltIPSTmqO3egfXv48UeIjAQ/P5gxA155RS7XxFcppRSa/KqHCA2VTWxz58qC2WuvSSeoyuVj6dtwMp1/msneXDeo+V4uduWNIta6A9e4Uvz+7hzqFKnzr+4N3QJ86LVgL5ExcffPc3O26sQ2lTh37sCuXVLa4OYGhw9Dq1bw4YdQvry9o1NKKZUKafKr7jNN2fg+ciSsXg2ZM0OHDtCpExQujFxYoQv9Du3nvSY5mF0OMC6BaZAtuh0lM79GXe86Cd72vbreoJWHCY+IxNPDjW4BPlrvq57O8eMwYYIMo4iJgfBwcHeX0cO6wquUUuoRNPlVREdD1+HXmDremTsXMuOS5S6tOkYxdmBWPDyQ/qeBXbj7+wqm1cvJZ69l54J5hfgiXgMnazTdXyz5yPtpWjG/Jrsqcfbsgd694ddfpTvDyy/Dxx/LNzXQxFcppdR/0uQ3HYuIkIEUI0bFcfVSNpxz3iBHg11kKh3O1gwW/tiZnwY/TyT624lMq+rCkP4enDEv45ffj4ZZA5i+fxg2MwaL4cynzzfTxFYlj2vX4PZt6cVrmrB9u7Qpa9cO8utzTiml1JPR5DcdOnVKRg9PmSIDrrIWu05u/yNkKHIZwwCX2BhabFpC1dE/Mrl0JEN6ZuK09RbV8ldgiv9A6nnXwzAM2latT/DJYPy9/PEr6Gfvh6Ucze7d8M03MGsWNGsGM2dKHe+ZM+DsbO/olFJKpVGJSn4Nw/AApgBlkYPg75mmGZIUgamkt3On1PP+9JMcHX7jDdnE1uzHEKlgME0CDofQZf13rPM6T5mPXLiQ2aRq/jJM9h9I/aL1H9jI5lfQT5NelfSWLIGgINiwQTaxtWghpQ33aOKrlFIqERK78jsWWGGa5iuGYbgAGZMgJpWEbDZYsUKS3rVrZU9Q587QsSMULCjX8VzhRtbD+2m8dwy/FT1BjTYWLmeCzIY3y98aQ0DRgH91b1AqSV2+DNmzSx3v+vUQFiZP2tat5XyllFIqiTz1kAvDMLIAuwFv8zFvRIdcpJy7d+Vo8ahRsl+tQAH45BNo2xayZv3bFS9d4s8POvFdxBxGPAemAZgGuWxtmNR0AM0qFbDXQ1Dpwa5dMG4czJ4Nv/wCAQFS35shA1it9o5OKaVUGpYcQy68gUvANMMwygM7gE6mad5OxG2qRLp2TQZbjRsH589LieQPP8Drr//jaHF0NLFfjeWHnwcwuEokx58hvnuDYRBQNqcmvip5xMXBokUyPWXDBsiYEd57D4oVk8szZbJvfEoppRyaJRF/6wRUAiaYplkRuA30/OeVDMNoZxjGdsMwtl+6dCkRd6ce5cQJ6cdbsKB0gqpQAVatkmEVLVs+mPjGLlvC940LUfJ0d94LiMTDuzRB9YJwc3bDalhxc3LlI7/G9nswyjHFxsq/Nps8Wc+ckdKGsDAYPx6KFrVvfEoppdKFxKz8ngXOmqa55a/f55NA8mua5mRgMkjZQyLuTyVg2zbJH+bPl6PEb70lm9jKlfv3deP272P28LcYlG0vR/ygYsai/NJ4NC+VeAnDMKhRsIZ2b1BJb98+ORSxdq3U4Dg7Q3AwFCmipQ1KKaVS3FMnv6ZpnjcM44xhGD6maR4G6gAHki409TA2GyxbJknv+vVSw9utm0xjS6jtadyVy/z4xdt8HrWCP4tBeYsnC5t9SZMyr2j3BpU8bDZYvhzGjJFxgRkywNtvS2+9bNniSxyUUkqpFJbYbg8dgFl/dXo4DrROfEjqYaKipH531Cg4fBgKFYLRo2UTm7v7v6+/8cR6vp77KZsvhXIqi41yLtn5OWAkTau+g8VITMWLUv9h3Tpo1Ei+jQ0bBu+/Dzly2DsqpZRSKnHJr2mau4B/7aJTSevKFZgwAb76Ci5ehEqVZHP8K68k3PLUZtoY+n1b+p+chmmA4Q6DS31Mr1fHadKrkkdYGHz9tWxe69cP/P1h4UJo2FD78iqllEpVdMJbKnbsmBw1/u47iIyEwEApb/D3lyEV/2QzbSzc+C2f/dqDfRmu3z/fYrFiyZdfE1+V9HbskCfp3LlS6vDOO3K+YUDTpvaNTSmllEqAZkOp0ObNsqpbvDhMniyT2Pbtg19/hRde+Hfia5omi/b8RKXP8/PKmg+IuX2Dgbxwv3uDi9UFfy9/uzwW5cAGDwZfX1i8GNq3h6NH5ZuaUkoplYrpym8qYbNJDjFyJGzaBB4e0LOnbGLLly/hvzFNk6WHl/DZ4s7sjDxO8Ssw804V3ugzhyURrhRbvoAz0dspmNGXC5cLQcGUfUzKwdy6Jclt7dpQtiy89JKUObRp84/JKUoppVTqpcmvnUVGwowZsontyBHw8pLe/++9B5kzJ/w3pmmy/OhyBizvwfZr+/C+Ct8f8eStT6fhVLc+i0LD6LVgL5Ex3mTFmxs3oNeCvQA0rZhAOwilHuXcOSk4nzABIiJg0CBJfsuXl5NSSimVhmjyayeXLklf/6+/hsuX5ejx3LnQvDk4/eO/yqLQMPovX8CZO9vxyJAZ10wbOXxjD14RMHWrG2+/OhjnkR3ubywKWnmYyJi4B24jMiaOoJWHNflVT6ZjRxkZGBcnT84uXaBaNXtHpZRSSj01TX5T2JEj0p5s+nRpXdaokWxie+65hDexLQoNo9OCuZy29gSnGCLiIFe4waS18G7F1rgsHA65cz/wN+ERkQne98POV+o+04SQEPDzkyekuzv873/wySc6gU0ppZRD0OQ3hfzxh9TzLloELi7QqhV07gylSj367/ou/4kwywgwYgCw2OCVgznYUnUg7b75KMG/8fRwIyyBRNfTwy3Rj0M5qJgYmDdPnqShoTKYonZtGDLE3pEppZRSSUq7PSSjuDj4+WdZRKtRQ/r+9+kDE5eeY7/3Ghp8v4waw9ewKDTsX3/7x5k/qDujLvtjPsUt5irOcWC1gcV0YmG57qxxL/zQ++0W4IOb84NjY92crXQL8Enyx6jSuMhIORRRtCi0aCG/T5kC1avbOzKllFIqWejKbzK4fVvKGsaMkV693t5S2/vuu/D7n/c2o0lNblhE5AOb0bac3cKA4AGsPLaS3BZ3hq52pf3maIY+/xzfVS6AxVIJV7PUI1dx79X1Bq08THhEJJ4ebnQL8NF6XxUvJkZqxC0WCAoCHx8pQm/QQM5TSimlHJRhmmaK3Zmvr6+5ffv2FLu/lHbhgiS548fD1atQtarU8zZtCta/FmJrDF+TYElC1ixnKOi1lGVHlpHDOSs9drvz0cKzRD5Tjfcqv8OebPF9ytycrQxrXk6TWfXkDh+WZHfDBmke7ewsuy9z5bJ3ZEoppVSSMgxjh2ma/5pErCu/SeDQITlyPGMGREdD48bQtauUOvxzE9s/N51FG8eJcJ7NqZjNnDntwdDrvrQfvx337Blh+iwyvfkm7+0K11VclThbtsAXX8QXnbduDXfuSH9eTXyVUkqlI5r8PiXTlMWzkSNhyRLIkEEmu376qRxBfph7m9FuWdZyw3k+MZZTWMxM1D/nx7xFB8lyORQ+6QL9+0OWLICUMWiyq57ahg1Qq5ZMTundWyan5Mlj76iUUkopu9Dk9wnFxsKCBZL0btsG7h5xFKhzGqPUUfblt3Lwjg8+PDxRfauGE91X9eOOEQqAYVqYsDQL7XaEgL+/1E2UKZNCj0Y5pNhY+OknuHlT2pTVqBE/J9vd3d7RKaWUUnalye9jujfZdcwYOHkSihWD//W+xtq4bdxF2pCFRTx8ktqxq8f4fP3nzNwzE4vVCjbAAEucjfM5b8KcOfD66wk3+1Xqcdy5I0/SUaPkSVq9OrRrJxvY3n/f3tEppZRSqYJu6/4P585Je7JChaBTJ/D0hIULpc53r3vo/cT3nnuT1O45FXGK9xe/j8/XPvy0/yc6V+3Eoqzv4xYL1jhwsThRb+xCWZXTxFc9rQULoHBhKWnw9IRffpFyB31OKaWUUg/Qld+HOHBAFtBmzpSuUM2ayWTXv7c/fdQktbAbYQzZMIQpO6dgGAYfPfsRvbI0JF/nfrBtG6ubVCL47Vr4V3kNv4J+KfSolEO5cEGaSXt6QsGCUKUK9Owp4wKVUkoplSBNfv/GNGUQRVAQ/PoruLlB27Yyia1YsX9fP6FJanFcIy7zQoqOW4bNtNGmYhv6lG9PgS8mwPhAGUU8axZ+b76Jn67Kqadx6pQUnU+ZAm+9BVOnwrPPwrJl9o5MKaWUSvU0+UX2B82fL/nEjh3S+WngQPjoI8iZ8+F/1y3Ah14L9hIRt487lm3EGheItG7GsMXSusK79H2uD17LNoFvbbh8Gdq3h0GDpL2UUk/q8GEYPlwOR4DMyO7Rw74xKaWUUmlMuk5+b96URbMvv5TFNB8f2RTfsqWs+v6XphXzE3ohmEFbemKacWBAxVz+/PTGtxQ7Hw3NWstScpUqsHw5VKqU/A9KOa6xY2HuXPlW1qWLFKIrpZRS6omkyw1v4eFSGlmwoJQ0FCok+4MOHJBN8Y+T+N64e4PP133OF9vbYCKJr9Ww8mpJf4qNmALly8OePTBpEoSEaOKrntzGjTJueMMG+X3AAOniMHasJr5KKaXUU0pXK7/79skmtlmzZJ/Qyy/LAlrVqo9/G7ejb/PV1q8I+iOIq5FXqVW4FlvPbiXGFoMLVvy7j4cdF2WC1hdf6PQs9WRME1avlvKY9eul7ub8eblMB1MopZRSiebwya9pwpo1Us+7YgVkzAgffACffALe3o9/O5ExkUzcPpFhG4dx6c4lGhRvwOf+n1PZszIhW+YT/G0f/Ff+iV+2PLDhZ6hZM/kelHJcTZrIyEBPT6nHef99edIqpZRSKkk4fPI7YQJ8/LEsmg0ZIolv9uyP//d3Y+8yZecUhm4cSvjNcOoUqcOgFwZJe7LoaBg6FL9Bg/BzcoKBo6TPqrNz8j0g5VhME1auhHr1wGqFhg0hMFCOHGTIYO/olFJKKYfj8Mnva6+Bqyu0aPH4uUTImRBWn1jN7ejbzN43m9PXT/NcoeeY3Xw2z3s9L1fasEEy6QMHoHlzqcMsUCD5HohyLDYb/PwzDB4steE//yzPo//9z96RKaWUUg7N4ZPfnDmhTZvHv/7G0xup830dom3RAJTKWYqVLVdSz7sehmHAlSvQvbuMkS1USA5RN2qUTNErhxMXJx0bBg+GgwelxciMGdC4sb0jU0oppdIFh09+H5fNtDH/wHzeX/zh/cQXDCrlbEz9ovXl8PT330PXrnDtmiTA/ftDpkx2jVulQZ99Bi4u8OOP8MorUu6glFJKqRSR7pNf0zRZ8ucS+q3tx54Le3Ay8yD/t9gwcGLjvjyscl5H3XEDpGevn5+0LytXzt6hq7QgJgZ++EEaSK9eLV+WVq+G/PnBki47DSqllFJ2lW6TX9M0+f347/Rd05dt4dsomq0oxZ16cfdmNaItfxJl2Yt7dEk+3biF54d1B/fMksC0aaNJi/pvsbGS9A4eDMePS5/nc+dkTnbBgvaOTimllEq30mXyu/7Uevqu6cuG0xsolLUQU16aQqvyrSjR5zcMwNVWijrHohj821d4RZxjYZkXaLbmR8id296hq7Tg8mWoVg2OHYOKFWWCyksvgWHYOzKllFIq3UtXye+Ws1vot7Yfvx//nbyZ8/J14Ne0rdQWVydXADw93LgbFk6/1VNocnAdx7N58tbrgzlVwY9mmviqR4mNhdBQePZZ2WVZvz4EBMhGNk16lVJKqVTD4ZPfkDMhzN43m9BzoWw6s4mcGXMyst5IPnz2QzI6/214gGnydeROvKcMJENMFF/WeJMJ1V7F4ubGsAAf+z0AlbrFxsKcOTKR7fRpOHEC8uWD8ePtHZlSSimlEuDQyW/ImRCen/48MbYYANpVasfI+iNxd3V/8IqHDsH//kfF9eu5XKkqHzz/AZtdcuHp4Ua3AB+aVsxvh+hVqhYXF5/0/vknlC8v3Rvy5rV3ZEoppZR6BIdOfoNPBmMzbQBYDSteHl4PJr5378Lw4TB0qOzCnzKFnK1bM0c3tKn/cvw4vPMOlC0LCxbIWGJ93iillFKpnkMnv0Z0GUzTCUwTEyeM6DLxF27YAO3ayarvm2/CmDEyA1mphJimbFzbvFm+MBUvDiEh4OurSa9SSimVhjjsp/ai0DC+D3Yh993BeMS2JPfdwXwf7MKydfsl6a1VC6KiYPlymD1bE1+VMNOEFStkI1uzZrBoEdy6JZdVqaKJr1JKKZXGOOzKb9DKw0TGxOFKKVxtpcA0qbN/LdW+/BYib8ikts8+0wlt6uEOH5a+zps2gZcXTJ8OLVqAk8O+bGqQ+UkAABE8SURBVJRSSimH57Cf4uERkfd/zn/9IoN+G0/t49vZk7cYOdavlv6rSiXk1i3InBmyZ4crV2DiRGjdWkYSK6WUUipNc9jk19PDjbCISF48vInRy0ZjYjCwzvusqv0qGzTxVQnZsQP694dLl2DLFsiVCw4c0D69SimllANx2ILFbgE+uDlbOZzLi+AilanXdjw/+jWnS2Bpe4emUpv9++Hll2XzWkiI/BwbK5dp4quUUko5FIdd+b3XmzdopQsfZ++Np4cMq9CeveoBK1dCYKCUOXz2GXzyCWTNau+olFJKKZVMDNM0U+zOfH19ze3bt6fY/SmVoAsXpE+vn5/0ev7iC/j4Y8iRw96RKaWUUiqJGIaxwzRN33+e77BlD0r9y/Xr0K8fFC0Kb70lU9pcXaXOVxNfpZRSKl3Q5Fc5vshICAoCb28YPBgaNpRyB6vV3pEppZRSKoVp8qsc36pV0L27DKrYsQPmzoUSJewdlVJKKaXswGE3vKl0zGaD+fOlR++HH0KjRtK6rEoVe0emlFJKKTvTlV/lWFatkhXe11+HGTMkETYMTXyVUkopBWjyqxzFoUMQEAD16smK74wZsHEjWPQprpRSSql4iS57MAzDCmwHwkzTbJT4kJR6AqYpK7t37kg976hR0rbM1dXekSmllFIqFUqKmt9OwEEgSxLcllKP5+pVGDIEbt2CSZOgUiU4cwbc3OwdmVJKKaVSsUQdEzYMowDQEJiSNOEo9R8iI2HECOnVO2aMjCG22eQyTXyVUkop9R8SWxD5JdAdsD3sCoZhtDMMY7thGNsvXbqUyLtT6dqmTeDjAz16QPXqsHs3TJ2qdb1KKaWUemxPnTUYhtEIuGia5o5HXc80zcmmafqapumbK1eup707lV6ZJty4IT8XKSKDKtauhWXLoFw5+8amlFJKqTQnMTW/NYDGhmE0ADIAWQzDmGmaZsukCU2le7t2QdeuUtqwdi14ekJwsL2jUkoppVQa9tQrv6Zp9jJNs4Bpml7AG8AaTXxVkggLg9atZRNbaCg0bx5f16uUUkoplQg64U2lLuvXQ2AgxMTAp59Cnz6QLZu9o1JKKaWUg0iS5Nc0zWAgOCluS6VDcXFw9iwULgy+vtCyJXTvLh0dlFJKKaWSkG6TV/a1apWUN9SpA9HRkDGj9O3VxFcppZRSyUCTX2UfBw5Aw4YyjvjGDRlY4exs76iUUkop5eC05lelvK1bpU9vpkwysKJDB8iQwd5RKaWUUiod0JVflTKio2HHXy2hfX1h8GA4ehS6ddPEVymllFIpRpNflbxMExYvhjJlpK43IkImsvXsCTr0RCmllFIpTJNflXz27YP69aFJE3Bygh9/BA8Pe0ellFJKqXRMa35V8jh+HCpUAHd3GDsWPvxQN7QppZRSyu505VclnZiY+PHD3t7SsuzoUejYURNfpZRSSqUKmvyqpLF8OTzzDNStCydOyHlt2kCOHPaNSymllFLqbzT5VYlz6BA0aCCn2FhYuBC8vOwdlVJKKaVUgrTmVz29a9ekbZnVCiNHSr9eFxd7R6WUUkop9VCa/KonY7PB2rXStixbNvjhB6hRA3LntndkSimllFL/Scse1OPbuRNq1pS63vXr5bxmzTTxVUoppVSaocmv+m9XrsAHH0iJw7Fj8N13kgQrpZRSSqUxWvagHs1mg+rVJent2BE++0wHVSillFIqzdLkVyVs2zaoVEk2s40eDYULQ9my9o5KKaWUUipRtOxBPejcOXj7bahSBaZPl/MaNvx/e3cfK1V953H8/V3EUpWKenVFClbsKhSbBZYa1vrUuvWpG3Rts4FsXaoYbaSJpJiI0pqatE2U0JilrcatStn40G7UXWto6tb0aW1pQXvFS/BZzLJlZddnahsEf/vH77B7O8xc7sBlzpk571cymTPnnPF+/eY3h0/O/OYcg68kSeoJhl9l77wDy5fDCSfA974HS5fC3LllVyVJkjSinPagbO5cuP/+fJb35pvhgx8suyJJkqQRZ/its61b4aCD4JBDYNEimD8f5swpuypJkqT9xmkPdbRzJ9xyC5x4ItxwQ1532mkGX0mS1PMMv3Wzbh3Mng1XXpmv5rBgQdkVSZIkdYzht06+9a18FYfNm+Huu+FHP4IpU8quSpIkqWMMv70uJdi2LS9//OP5RhVPPQXz5kFEubVJkiR1mOG3lw0MwBlnwGc/m19PmZKv5HDooaWWJUmSVBbDby96+2245hqYPh02bIBzz81ngCVJkmrOS531mt/8Bj71KXjxRbjkErjpJujrK7sqSZKkSjD89oqU8hzeiRNhwgS488485UGSJEn/x2kP3e7dd+G22+C88/L1e/v64Oc/N/hKkiQ1YfjtZgMD+eYUV1wB27fD66+XXZEkSVKlGX670R/+AEuXwowZ8PTTsHIlPPIIHHFE2ZVJkiRVmuG3W913H3zmM/mavfPne81eSZKkYTD8doutW+Gqq/INK8aMgbVr84/avJKDJEnSsBl+qy4lWLUKpk6FW2+FRx/N68eOLbcuSZKkLmT4rbJNm/INKubPz3dn6++Hc84puypJkqSuZfitsoUL4Re/gBUr8uXLpk4tuyJJkqSu5k0uqmbDBjj8cBg/Hr7xDRg1CiZNKrsqSZKknuCZ36rYvh1uuCFfvuzaa/O6444z+EqSJI0gz/xWwZo1cNll+azvvHmwbFnZFUmSJPUkz/yW7e674ZRT4I034KGH8usjjyy7KkmSpJ5k+C3L73+fn88+GxYvzmd9P/nJcmuSJEnqcYbfTnv9dViwAM48E3buzDepWLYM3ve+siuTJEnqeYbfTlq9Gk46CVauhLPOyuFXkiRJHWP47YS33oJLLsnTGsaNyz9w+9rX4MADy65MkiSpVgy/nTB6NKxbB0uXwmOPwUc+UnZFkiRJtWT43V9eew2uvhq2bYMxY3Lo/cpX4D3vKbsySZKk2trr8BsREyPixxGxMSI2RMRVI1lYV/v+92HaNLj5ZvjpT/M6pzhIkiSVbl/O/O4AFqeUpgKzgYUR8aGRKatLvfoqXHwxzJmTr9X76197+TJJkqQK2evwm1LaklJ6vFh+C9gITBipwrrSFVfAvffC9dfD2rUwc2bZFUmSJGmQSCnt+38k4gPAz4CTUkpvNmy7HLgcYNKkSX/x0ksv7fPfq5TXXoMdO/KZ3hdeyHdqmzGj7KokSZJqLSIeSynNaly/zz94i4hDgPuARY3BFyCldFtKaVZKadaRvXbb3ocfhg9/OJ/xBZg82eArSZJUYfsUfiNiNDn43pVSun9kSuoCv/sdXHklnHNOvjPbddeVXZEkSZKG4YC9fWNEBHA7sDGl9PWRK6niBgbgwgvzFIcvfCFfvuy97y27KkmSJA3DXodf4KPAxcCTEdFfrLsupbR638uqsKOPhqOOgttvhzPOKLsaSZIktWGvw29K6d+BGMFaqqu/P1+z99vfhr4+ePRRiHr8r0uSJPUS7/A2lB074KtfhZNPhh/+ME91AIOvJElSlzL8tvLMM3DqqfDFL8JFF+W5viecUHZVkiRJ2gf7Mue3d6UE8+bBiy/CPffA3LllVyRJkqQRYPgdbMuWfOmygw+GVavgsMPgmGPKrkqSJEkjxGkPuzzwQL5hxZIl+fW0aQZfSZKkHmP43bYNLrssz+s99lhYuLDsiiRJkrSf1Dv89vfD9Olwxx1w7bXwy1/ClCllVyVJkqT9pN5zfseOhTFj4Cc/gdNPL7saSZIk7Wf1O/P7/PPwpS/lKzocfzysX2/wlSRJqon6hN+U4M478zSHFSvyZcwA/qQ+LZAkSaq7eiS/V16BT38aLr0UZs3KZ3snTy67KkmSJHVY78/5TQk+8Yl8h7Ybb4TFi2HUqLKrkiRJUgl6P/xGwPLlMG4czJhRdjWSJEkqUe+HX4CPfazsCiRJklQB9ZjzK0mSJGH4lSRJUo0YfiVJklQbhl9JkiTVhuFXkiRJtWH4lSRJUm0YfiVJklQbhl9JkiTVhuFXkiRJtWH4lSRJUm0YfiVJklQbhl9JkiTVhuFXkiRJtWH4lSRJUm0YfiVJklQbhl9JkiTVhuFXkiRJtWH4lSRJUm1ESqlzfyziv4GXOvYH/18f8D8l/N1uZb/aZ8/aY7/aY7/aY7/aY7/aY7/aU2a/jk0pHdm4sqPhtywRsS6lNKvsOrqF/WqfPWuP/WqP/WqP/WqP/WqP/WpPFfvltAdJkiTVhuFXkiRJtVGX8Htb2QV0GfvVPnvWHvvVHvvVHvvVHvvVHvvVnsr1qxZzfiVJkiSoz5lfSZIkqbfCb0ScGxFPR8RzEbGkyfaIiH8otq+PiJll1FkFETExIn4cERsjYkNEXNVknzMj4o2I6C8e15dRa1VExKaIeLLoxbom2x1fhYg4cdC46Y+INyNiUcM+tR9fEXFHRGyNiIFB6w6PiH+LiGeL58NavHfI410vatGvZRHxVPGZeyAixrV475Cf317Uol9fjoj/HPS5O7/Fex1fed13B/VqU0T0t3hvHcdX0xzRFcewlFJPPIBRwPPAZOBA4AngQw37nA/8AAhgNvCrsususV/jgZnF8ljgmSb9OhN4qOxaq/IANgF9Q2x3fDXvyyjgv8jXWxy8vvbjCzgdmAkMDFp3E7CkWF4C3Niip0Me73rx0aJfZwMHFMs3NutXsW3Iz28vPlr068vA1Xt4n+Or+fblwPUtttVxfDXNEd1wDOulM78nA8+llF5IKW0H7gUuaNjnAmBVytYA4yJifKcLrYKU0paU0uPF8lvARmBCuVV1PcdXc2cBz6eUyrjBTaWllH4GvNqw+gLgO8Xyd4ALm7x1OMe7ntOsXymlh1NKO4qXa4D3d7ywimoxvobD8dUgIgL4W+CejhZVYUPkiMofw3op/E4A/mPQ683sHuaGs0/tRMQHgBnAr5ps/suIeCIifhAR0zpaWPUk4OGIeCwiLm+y3fHV3Fxa/4Ph+Nrdn6aUtkD+xwU4qsk+jrXmLiV/+9LMnj6/dfL5YprIHS2+knZ87e404OWU0rMtttd6fDXkiMofw3op/EaTdY2XshjOPrUSEYcA9wGLUkpvNmx+nPxV9Z8DK4B/6XR9FfPRlNJM4DxgYUSc3rDd8dUgIg4E5gD/3GSz42vvOdYaRMRSYAdwV4td9vT5rYtbgOOB6cAW8lf5jRxfu5vH0Gd9azu+9pAjWr6tybqOjbFeCr+bgYmDXr8f+O1e7FMbETGaPGDvSind37g9pfRmSmlbsbwaGB0RfR0uszJSSr8tnrcCD5C/thnM8bW784DHU0ovN25wfLX08q7pMsXz1ib7ONYGiYj5wF8Df5eKCYWNhvH5rYWU0ssppZ0ppXeBf6R5Hxxfg0TEAcBFwHdb7VPX8dUiR1T+GNZL4Xct8GcRcVxxtmku8GDDPg8Cf1/8Kn828MauU/N1U8xfuh3YmFL6eot9ji72IyJOJo+XVzpXZXVExMERMXbXMvlHNgMNuzm+dtfybInjq6UHgfnF8nzgX5vsM5zjXS1ExLnANcCclNLbLfYZzue3Fhp+h/A3NO+D4+uP/RXwVEppc7ONdR1fQ+SI6h/DOvXLuk48yL+2f4b8C8KlxbrPAZ8rlgP4ZrH9SWBW2TWX2KtTyV8xrAf6i8f5Df36PLCB/CvMNcApZdddYr8mF314ouiJ42vPPTuIHGYPHbTO8fXHPbqH/NXzO+QzIQuAI4BHgGeL58OLfY8BVg96727Hu15/tOjXc+S5g7uOY7c29qvV57fXHy369U/F8Wk9OWyMd3y17lexfuWu49agfR1frXNE5Y9h3uFNkiRJtdFL0x4kSZKkIRl+JUmSVBuGX0mSJNWG4VeSJEm1YfiVJElSbRh+JUmSVBuGX0mSJNWG4VeSJEm18b9XcejfYa2ThgAAAABJRU5ErkJggg==\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.40747839 0.40758736]\n",
      "[0.41236678 0.03553119]\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": [
      "[4.83396446 0.51155795]\n",
      "[0.10404716 0.00896512]\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3hU1dbA4d+ZSaUk9N57laqCiFIFAQFRBJQO0vwUr4pdUJCiYAEhFAEB6dJBEWlRxEgNvfckVCEBQnrmfH8sAoKB9CnJep/HBzIzyezcezhr9t5rr2WYpolSSiml7Mvi6AEopZRSWZEGYKWUUsoBNAArpZRSDqABWCmllHIADcBKKaWUA2gAVkoppRzAzZ5vli9fPrNUqVL2fEullFLKYXbt2vWPaZr5E3vOrgG4VKlS7Ny5055vqZRSSjmMYRhnH/ScLkErpZRSDqABWCmllHIADcBKKaWUA9h1DzgxsbGxBAcHExUV5eihZBgvLy+KFSuGu7u7o4eilFLKSTg8AAcHB5MzZ05KlSqFYRiOHk66M02Tq1evEhwcTOnSpR09HKWUUk7C4UvQUVFR5M2bN1MGXwDDMMibN2+mnuErpZRKOYcHYCDTBt8Emf33U0oplXJOEYCdyaeffsq4ceMe+PyKFSs4dOiQHUeklFIqM3K5ALwiMIQGYzZR+v2faTBmEysCQ+z7/hqAlVJKpQOXCsArAkP4YNl+QsIiMYGQsEg+WLY/zUF45MiRVKxYkWbNmnH06FEAvv/+ex599FFq1KjBCy+8QEREBH/99RerVq1iyJAh1KxZk5MnTyb6OqWUUiopLhWAx647SmRs/D2PRcbGM3bd0VT/zF27drFw4UICAwNZtmwZO3bsAKBDhw7s2LGDvXv3UrlyZWbMmMETTzxB27ZtGTt2LHv27KFs2bKJvk4ppZRKisOPIaXE+bDIFD2eHFu2bOH5558nW7ZsALRt2xaAAwcO8PHHHxMWFkZ4eDgtWrRI9PuT+zqllFJO7to1yJPHbm/nUjPgIrm8U/R4ciWWpdyzZ08mTpzI/v37GTZs2AOPESX3dUoppZzUhQvQoQMULSp/txOXCsBDWlTE2916z2Pe7laGtKiY6p/51FNPsXz5ciIjI7l58yarV68G4ObNmxQuXJjY2FjmzZt35/U5c+bk5s2bd75+0OuUUko5sagoOHhQ/p47Nxw9Cm+9BXasWOhSS9DtaxUFZC/4fFgkRXJ5M6RFxTuPp0bt2rXp1KkTNWvWpGTJkjRs2BCAESNG8Pjjj1OyZEmqV69+J+h27tyZV199lQkTJrBkyZIHvk4ppZQTCg6GyZNh2jTIkQNOnAAvLzhwAOxcs8EwTdNub1a3bl3z/n7Ahw8fpnLlynYbg6Nkld9TKaWc0r59MGoULFkCNhu0aweDB8PTT2do4DUMY5dpmnUTe86lZsBKKaVUskVHQ0wM5MwJZ87Ar79K0P2//wMnqM3vUnvASimlVJIuXoRPP4WSJeGLL+Sx1q1l+fmrr5wi+ILOgJVSSmUWO3fC+PGwaBHExkKrVtC0qTxntcqerxPRAKyUUsp1xcdLcAWZ7f76KwwYIMvMFSo4dmxJ0CVopZRSrufKFRg5UpaZDx+Wx77+WpaZJ0xw+uALOgNWSinlSvbskQA7f74kWT3zjCRaARQv7tixpVCSAdgwjJlAG+CyaZrV/vX468D/AXHAz6Zpvptho8xAV69epentPYKLFy9itVrJnz8/ANu3b8fDw8ORw1NKKZXgxg144gk5NtS7tywzV6ni6FGlWnJmwLOAicCchAcMw2gMtAMeMU0z2jCMAhkzvIyXN29e9uzZA0gv4Bw5cvDOO+/ceT4uLg43N10oUEopu7t6FaZPh+3bYelS8PGBFSvgsccgVy5Hjy7Nkowspmn+YRhGqfseHgiMMU0z+vZrLqf/0BynZ8+e5MmTh8DAQGrXrk3OnDnvCczVqlVjzZo1lCpVirlz5zJhwgRiYmJ4/PHH8fPzw2q1JvEOSimlHmj/fllmnjtXSkY2aQI3b8p53meecfTo0k1qp3YVgIaGYYwEooB3TNPckdbBvPmmLO+np5o14dtvU/59x44dY8OGDVitVj799NNEX3P48GEWLVrE1q1bcXd3Z9CgQcybN4/u3bunbdBKKZVVrVolVaq8vaF7d3j9dahWLenvc0GpDcBuQG6gHvAosNgwjDJmInUtDcPoB/QDKFGiRGrHaXcdO3ZMcia7ceNGdu3axaOPPgpAZGQkBQq47Gq8UkrZX2gozJgBhQvDK69As2Ywbhz06mXX1oCOkNoAHAwsux1wtxuGYQPyAVfuf6FpmtOAaSC1oB/2Q1MzU80o2bNnv/N3Nzc3bDbbna8TWg6apkmPHj0YPXq03cenlFIu7cAB+O47+PFHiIyUpKpXXoFs2eDttx09OrtI7TngFUATAMMwKgAewD/pNShnU6pUKXbv3g3A7t27OX36NABNmzZlyZIlXL4sW+DXrl3j7NmzDhunUkq5hHfegerVYc4cCbp798osOItJzjGkBUAjIJ9hGMHAMGAmMNMwjANADNAjseXnzOKFF15gzpw51KxZk0cffZQKtw94V6lShc8//5xnnnkGm82Gu7s7kyZNomTJkg4esVJKOZHQUJg5E7p2hYIFpTxk/vzQty/kzevo0TmMtiO0k6zyeyql1B2HDsky85w5EBEhR4r69HH0qOxK2xEqpZSyn7g46T7022/g6SnLzK+/LsdS1B1aC1oppVTahYXBsmXydzc3KFMGRo2S2swzZmjwTYTOgJVSSqXev5eZIyPh3DkoVgwmT3b0yJyezoCVUkql3JEj0Lw5VK0KP/wAnTrBrl0SfFWy6AxYKaVU8oSGwuXLULEi5M4Np0/LMvOrr0K+fI4encvRAKyUUurhDh68WzSjdm3YskWOEx0/Lp2JVKroEjQQHBxMu3btKF++PGXLlmXw4MHExMTg7+9PmzZt/vP6NWvWUKtWLWrUqEGVKlWYOnWqA0atlFIZbNMmObNbrRrMng2dO0sgTqDBN02yfAA2TZMOHTrQvn17jh8/zrFjxwgPD+ejjz5K9PWxsbH069eP1atXs3fvXgIDA2nUqJF9B62UUhklNFSSqUASrI4fh9GjIShIs5nTWZYPwJs2bcLLy4tevXoBYLVa+eabb5g5cyYRERH/ef3NmzeJi4sj7+3qLZ6enlSsWNGuY1ZKqXR34AD07w9Fi8pSM8je7qlT8P77usebAZxrD9gB/QgPHjxInTp17nnMx8eHEiVKcOLEif+8Pk+ePLRt25aSJUvStGlT2rRpQ5cuXbBYsvxnGaWUqzFNaXD/3XeweTN4eUnRjAYN5HlPT8eOz45WBIYwdt1RzodFUiSXN0NaVKR9raIZ+p5ZPmqYpomRyD7Ggx4HmD59Ohs3buSxxx5j3Lhx9O7dO6OHqZRS6ed2RzcMA8aMgZMn5c/gYCkXWbWqY8dnZysCQ/hg2X5CwiIxgZCwSD5Ytp8VgSEZ+r7ONQN2QD/CqlWrsnTp0nseu3HjBkFBQZQtW/aB31e9enWqV69Ot27dKF26NLNmzcrgkSqlVBrt3Suz3eXL4dgxaYSwdCkUKiTVq7KoseuOEhkbf89jkbHxjF13NENnwVl+Bty0aVMiIiKYM2cOAPHx8bz99tv07NmTbNmy/ef14eHh+Pv73/l6z5492v1IKeW84uJgyRJ4+mnZkps/H158EWJi5PlixbJ08AU4HyZJZ9GWw1x3W0y05fA9j2eUrP2/OmAYBsuXL2fQoEGMGDECm81Gq1atGDVqFAEBAWzcuJFi/6rssmDBAr788kv69++Pt7c32bNn19mvUsr5mKYsMZ88CR07QqlSMG6cNL7PndvRo0tXad2/LZLLm5M3dnHJ40MgHgN3CsaMpIxP7YwbNBqAAShevDirV6/+z+ONGjUiMvK/n4AaNmxoj2EppVTK7doly8xxcTB3rlSt+vNPqFcPrFZHjy7dJezfJiwhJ+zfAskKwqZp0rhmMNv+Hg1G3O3H4oh3O8iQFl0ybuDoErRSSrm+2FhYuFCyl+vWlSXnvHllFgzyeCYMvvDw/duk+J/xp8HMBozc1gOPeC+I9wCbFYvhzltPP5/hWdA6A1ZKKVf35Zfw8cdQrpwks/bsCb6+jh6VXTxon/Zh+7fbQ7bz0aaP2HBqAz4UJfumqdz6sxdPdtpJzfb+vFy/EfWL18+oId+hAVgppVyJacK2bbLM3KULtGkj+7q1akHLlpDFahIUyeVNSCLBtkgu7/88NnfvXEb9OYrD/xwmu5GPHH9+zQ3/gTzb3Ivh26Bu3fpAxgfeBBqAlVLKFURFweLFEnh37oScOSEhH6VwYfkvCxrSouI9e8AA3u5WhrS4W6HwxLUTvPbza/x26jd5IN6NW3MX0KxsM4b/AfXtF3PvoQFYKaVcQePG8PffUKkSTJwI3btLEM7iEvZpE8uCDroexIg/RjAzcCamaYBpgGGCYfLqZzuY1q2ZQ8euAVgppZyNaUrm8g8/gJ+flIh8/33Inl26E2kXonu0r1X0noSpy7cu89a6t/Db4UdcvI1shwZyM/AZLJ07gTUGTw8PejlBEx0NwEgDhurVqxMXF0fp0qX58ccfyZUrF2fOnKFNmzYcOHDgntf37NmTxYsXc+nSJXLe/gQ6ePBgJkyYwJUrV8inRcuVUqkRESGFMiZOlKpVuXPDa69BnTrQrp2jR+cwyTnnGxAUwNoTazl3/RxLDi0hIjaSHCd6EL1mKJXKl2L4BPCpspHfz/rTqJR9kqySogEY8Pb2Zs/tJhA9evRg0qRJD2xHmKBcuXKsXLmSrl27YrPZ2Lx5M0WLZmzKulIqEzt3TipVhYbCI4/A99/Dyy9DIhX5spLknPPddGoTLee1JNYWC0C2C00wl06ibLFKDJ8jeWqyaFCfJ0o4PvAmcMl0uYCgAEZvGU1AUEC6/+z69esTEpJ0Ae4uXbqwaNEiAPz9/WnQoAFuWbycm1IqBUwTNmyAadPk6+LFJZv599+lK1zfvlk++MLDz/lGx0Xz3bbvaLew3Z3gi81Kzn+asXRqJXbtgueec94Ve6eKGG/++iZ7Lj68HeH16Ovsu7QPm2nDYlh4pOAj+Ho++LxbzUI1+bZl8po8xMfHs3HjRvr06ZPka8uXL8/KlSsJDQ1lwYIFdO3albVr1ybrfZRSWdjNmzBnjiwzHzkCpUtDnz5SKGPcOEePzukkdp7XJJ6jN1dSYWIvzl0/h1doTch5GCxxeFg9WPp1IxqUcMBgU8ipAnByXI+6js20AWAzbVyPuv7QAJwckZGR1KxZkzNnzlCnTh2aN2+erO/r0KEDCxcuZNu2bUydOjVNY1BKZQGLFkmT+5s34dFHJRB37Jhpq1Slh4RzvtGWw0Ra9gE2bln9ibOEEHm2LqyaTlFLM1557288KvjTpLRz7O8mh1MF4OTMVAOCAmg6pykx8TF4WD2Y12Femv/HTtgDvn79Om3atGHSpEm88cYbSX5f586dqV27Nj169MCSxQ6/K6WSIT4e1q6FkiWhenWoUkWSqf7v/+Dxxx09OpcwpEVFBi9bxEXrB0AcGMCNkvDLcvJHtWPoJwbdu4Obm32LaKQHpwrAyVG/eH02dt+I/5n0z2Tz9fVlwoQJtGvXjoEDByb5+hIlSjBy5EiaNXPsWTKllJO5du3uEaJTp2DgQPl79erw44+OHp1L8c11jIgcX0GUNErAZsH7cB++ebM9vXqBh4djx5cWLheAQYJwRi0x1KpVixo1arBw4UIaNmzI0aNH72lH+M0339zz+v79+2fIOJRSLurdd2V/NzISnnoKRo+G55939Khczr/rNbvH5AOrOxg23C0erJ3YjKfLOHqEaeeSATi9hYeH3/P1v1sTxsbG/uf1HTt2TPTnnDlzJl3HpZRyAbGx8Msvkm5rsUjmcteucn63Rg1Hj87l7L+0n082f8LKoytxj80HG78m19mBdPlfIHnr+tO8nOvs8SZFA7BSSqXGxYtyVnfKFDh/Htavh2bN4NNPHT0yp5OcQho/HfyJkVtGsvfSXtzifeD3EeQ8Npj3/peT116D7Nldb483KRqAlVIqJa5ehcGDpTFCbCw88wxMnSq1mtV/JFVII+h6EP/3y/+x6tgqMAGbGx4rF/PpSy14Y03mLnetAVgppZISGQnHjsmSsq+vFMoYNEj+q1DB0aNzag8qpDHy17/5/dKf+O2YTGx8LCCNEgyLydvjdvNR8xaOGbAdJRmADcOYCbQBLpumWe2+594BxgL5TdP8J7WDME0Tw1lLlaQD0zQdPQSlVGqcPQuTJ8P06ZJue/YsuLvDvn1Zru9uat1fSMNGONfdlnMuaiW7/o7B3NMTj7OtsbXrimnE4OHmwbOVGjlmsHaWnBnwLGAiMOffDxqGURxoDpxLywC8vLy4evUqefPmzZRB2DRNrl69ipeXl6OHopRKrsBAGD4cVq2Sr9u3l7O7CeVmNfgmW0IhjUjLXm64LSHKOAKWSDjYEc+/hvPGy5UY8j2ciMqY46XOLMkAbJrmH4ZhlErkqW+Ad4GVaRlAsWLFCA4O5sqVK2n5MU7Ny8vrnqNMSikndPMmxMRA3ryyz/vnn/DeezBgAJRwgbqGTurNZqXov+ptQi0rpIiGzQIrZvJclY5M25mDQoXkdfnJuOOlzipVe8CGYbQFQkzT3JvUrNUwjH5AP5DCFfdzd3endOnSqRmGUkql3ZEjMGkSzJ4tZSK/+kp67gYFSR9elSpxtjhm75nNUP/PCLUGSYIVAAYdXj/O0t45HDk8p5DidRTDMLIBHwFDk/N60zSnmaZZ1zTNuvnz50/p2ymlVMb4+Wdo3hwqV5aORO3bQ5cu8pxhaPBNJZtpY+GBhVSaUJW+q/ty4XghjHXjsZreWA0r3h4evNP8OUcP0ymkZgZcFigNJMx+iwG7DcN4zDTNi+k5OKWUSldhYZArl/x98WKZ/Y4cKa3/ChRw7NicVHLO8AYEBbD5zGbcLe7MDpzLwav7MK5Uw7J5Od0fb8cncw0uuT+a5fZ4k5LiAGya5n7gzpVqGMYZoG5asqCVUipD7dwp5SEXLoS//oLateGbb8DH525ilfqPpM7wggTfxrMbEx0fLd90oyisn0eXGp34dIWV8uXl4TJZcI83KUkuQRuGsQAIACoahhFsGEbSzXKVUsrRYmJg7lyoV09a/y1ZIg3v8+SR5/Pk0eCbhAed4R277igA24K38crSbneDr81CtcgBHFr0MvN+vBt8VeKSkwXdJYnnS6XbaJRSKq1iYuTMbnS0FMooUgQmTIDu3aWIhkq2+8/wJjhz/RCt547nl5OrIMpXzkZbbHi5ezDtg6ZULm7ngboo/finlHJ9pgmbN8sy86lTco43Z05Zei5fXpKqVIolnOGNthwmyrIfN1shIizbiHD7gwuHfGDrCFrmGkyn1w9wwVP3d1NKA7BSynXdvAlz5sgxosOH5Qzvq6/K7NfLS8tEptGQFhV5c9liLlk/wiRGHozzhD/fo5HXEMaMzcNjj4E0SdDAm1IagJVSrsc0ZVa7apVUqKpbF2bNgk6d9PhQOnqigjue+eZjhsZIEQ3TIP/ZN1j+2WgaNHD06DLAzp0QHCxH0uxA66kppVxDXBwsXy5FMr78Uh578UXYtg127IAePTT4ppOwqDA+WP8xJb4qw7Fru8C0gs2Kp9WLlaOez1zBNzJSirA89pgk6w0ZIh/w7EBnwEop53b5sjRDmDJFqlMVLy4zXQBPT26vgap0cCvmFt8GTGDU718SYYbBgU7UCP2M3oOuEZ7fn8aZaY/31ClptDFzJly7JgVZvvsOunWzW86ABmCllHPr0wfWrJFm9xMmQJs2enwonf1+5nfG/jUO/5NbuWULhaNtqHJpBOPeqUnLlgnxKBME3vh4+PVXyRn49VdpqtG+Pbz2GjRqZPdkPb2KlVLOIzJSimVMmQI//SRNEEaOlCXnypUdPbpMJ84WxyebhjJm6xjABJuFonun4te7H889l4mSx//5B2bMkOvqzBkoVAg++QT69YOiRZP89oyiAVgp5Xj3LwdWqQLnz0sAfuQRR48u07GZNhYf/Im3Vw/lfMwxaZRggMViMOitq7R9ytEjTAemCdu3y2x38WLJjH/6afjiC3j+eTm77GAagJVSjhUWJrPb+Hi5Mb72mtwoM830y3mYpsmaYz8zeOXHnI7cC5eqke/8GG7U/ox4YvCwetC4dCNHDzNtIiJkFWXSJNi9G3LkkFrfAwdC1aqOHt09NAArpewrNBR++AEOHpRlwVy55CxvgwagfbMzREBQADN2z2TD0QDORh6Ea2UpcHAeX3TrRNdXrOy48JTrN0o4flyWmH/4Qa6xqlXBzw+6dpWiLE5IA7BSyj4CA2VWMn++7PU2aCB/envfzWpW6W76rhn0W9MPExuYBjn3v8vYZz+n9zj3O6uw9Yu7aKOE+HhpK+nnB+vWSXLeCy9ICdKGDZ1+FUUDsFIq482bJzMRb2/5c9AgqFnT0aPK1PZf2s+gJZ/w5z8r7+7xGhbefT0X/Rs5fv8zTS5fvptUde6cJFINHy5V0AoVcvTokk0DsFIq/QUFwdSpUKMGdOwIrVrB119Dz56QO7ejR5epHb96nNeWDGP9hYUQ7YP3qX7EVfkRmyF7vE3LNnL0EFPHNCEgQGa7P/0kTTeaNpW2km3buuTRNNcbsVLKOZkmbNoky8wrV8rX77wjATh3bvjf/xw9wkwpICgA/zP+VM5XmVkBP7Pq3A+YsZ5473+fD59+h/JtIxm5uRpBETspnq0ul/4pAa7UrejWLdm28PODPXukh/OAAZJUVamSo0eXJhqAlVLpo3NnOe6RN6+U8xswAEqVcvSoXNqKwBDGrjvK+bBIiuTyZkiLirSvdffcakBQAE3mNCEqLlo+8MS74XlgEO88/iEfLCjE+mMhfLBsP5GxZfClDDduwAfL9gPc83Oc0tGjcjRt1iy4fl2Oo02dCq+8AtmzO3p06UIDsFIqdQ4ckD24zz+XTOZu3aB1a3jpJa3JnA5WBCYEz3gAQsIi7wmeoZGhDF71IVGxUdIoAYMn3f/Hz7O/xMdHfsbYdUfvfH+CyNh4xq476pwBOC4OVq+WVZSNG+WsbseOcjStfn2nT6pKKQ3ASqnki4mBFSvkBvnHHxJo27aFZ56REpEq3TwoeI75dQ+/n57NpD1jibWGgWnFALzcPfiy+/N3gi/A+bDIRH/2gx53mIsXpd731KnSjah4camA1qcPFCzo6NFlGA3ASqnkSThbeeEClC4t5SF795YlZ5Xu7g+SJrGExWzk7I2FbNv/D9bTbehdYgQvvRzJ7muJn+EtksubkESCbZFc3hk69mQxTfjzT/kwt3SpzH6bN4eJE2UlxQWTqlIq8/+GSqnUMU3w94d9+2DwYEmk6tFDzle2bCmF7FWGSQieUZaDXDdXE8UR8P0HzjzNywVG8+3X9cmfX17b4gGNEoa0qHjPMjaAt7uVIS0q2uNXSNzNmzB3riRVHTgg2xevvy45AxUqOG5cDmCYdup7CFC3bl1z586ddns/pVQq3LgBP/4oN8hDh+Rc5enTuq9rZ8t2B9F34TBCs80CwwSbgffhT5jYtz+9nymS7J+TVCKX3Rw6JNfUnDkShGvVkr3dLl0gWzb7j8dODMPYZZpm3cSe0xmwUuqulSulUEZ4ONStK2X9OnXS4JtBEguO7WoWYc62Nbz988eEZt8nRTQADIMXXoxJUfAFSdhyWMJVbKzkDPj5yWqKh4dcT4MGweOPZ7qkqpTSAKxUVhYbC8uXSw3mJ56Q6lQdOsjMRBvdZ6jEspxfm7+Q3osWEeq9A26Uo7bbpxzK/QWxthg83DwYVL+tYwedXOfPw7Rp8t+FC1CyJIweLUlVCevmSgOwUllSSIjcHL//Xm6QPXtKAC5ZEmbPdvTosoSELOdoy2Fumr8TGXUWW679cL0YT9z6npmDe1CxnDsBQc+4RqME04Tff5ekquXLwWaTXIFp0+DZZ8FqdfQInY4GYKWymjfflExTm01ujNOnQ4sWjh5VlnM+LJIbNn9Cvb4Gwwae4Ha0DwVyvMTWac/ceZ3TN0q4P2cgTx6pejZgAJQt6+jROTUNwEpldtevSym/Pn1kD658eblBDhwIZco4enSZ1sOSn3afOc7lm1OJyL/6X99hIUeZWErlcJHb8v79EnR//FHKRf47Z8DbCY45uQAX+X9aKZVi+/fLcuDcuXKDLFlSmiK89pqjR5bpPaiKVfC188zeOo2d8T9ALk+sp18kvtQqIA4DN3yNmo49IpSUmBhYtkwC75Yt4OkpWcyDBsGjjzp6dC5HA7BSmc3Vq9C+vRQ58PK6e4Osm+hJiCwrI4/n/Ht/N8qyH7eYMly4eJLX/RcBUOrqa0zr9iG3vOIYuvapO40Shj/bwTlLRAYHS5Wq77+HS5dk5WTsWOjVSwuxpIEGYKUyg3PnZMbburXsweXKBePGyQ0yTx5Hj87pJFVnOa3Oh0USbTnMJY+PMIkFNxO8LVhPdGbpq6Np16jEnde2r/V6mt8vQ5im1GP284NVqyRnoHVr+TDXooUWYkkHGoCVclU2m9wgJ02SAva+vlJT18NDvlYPlNFNCvJnNwi8uhzTI0YaJZiQ7XpHalbvc0/wdUphYZIJP3mydCTKm1faSvbvLyVIVbrRAKyUK1q/XvZyjx+Xc5XvvQf9+knwVUnKqCYF4ZHR9J0ylZ3hIzFzXwabzBIN3Mif83Hn3t/du1c+zM2bBxERUK+eVK3q2FELsWQQDcBKuYrduyFnTslizpMH8uWDYcPgxRclGUYlW3o0Kfj3HnKhHJ5kizvE77fGEZ8jCJ+IJnQsOoRtN48THOnE+7vR0dIIYdIk+OsvyV5++WVZZq5d29Gjy/Q0ACvlzKKipMm9nx9s2wZ9+0oiTJ06csNUqZLWJgUJe8jXYvcTFrGRM7EHwCcEz+jHGF7nBz4Y2vR2lcWWGfMLpNXZs5JUNX06XLkC5crB119LQZbcuR09uiwjyQBsGMZMoA1w2TTNarcfGws8B8QAJ4FepmmGZeRAlcpyPv8cvv1WsporVYLx46F7dziFV/wAACAASURBVEePKlNImImmNgv6i1+OcO7ir0QUnQDeJpgG2c6/ziMl2/Fh56YZOfTUs9lk68LPD9askcfatJGtjGbNNKnKAZIzA54FTATm/Oux9cAHpmnGGYbxBfAB8F76D0+pLCQ+HjZtkpuhYUgyTOPGshzYqFGWL1yf3lLTpMA0YcTcTWwL/gCz2I67jRIw8Mh/nYs3otJ9nGkWGioFMiZPhhMnJGfg/fclqaqEkyeEZXJJBmDTNP8wDKPUfY/99q8v/wZeTN9hKZWFXL4MM2fClCmyNLh5swTcsWM16DoJ04SvF2/j0y0fEZ5/I2QvhveFXkQVmodpShENL1t152h0n2DXLpntLlgAkZHQoAEMHy7NNjRnwCmkxx5wb2BROvwcpbKWf/6Rusw//SQVhpo0kbO7DRrI8xp8He6vcwF88ctCNh3eRXierVhzFOCVPN/SukUbhv98nLCYKkRZ9uNlq04uazXHZznfnzOQLRt06yarKDVqOHZs6j/SFIANw/gIiAPmPeQ1/YB+ACV0uUNldeHhcOyYZJj6+MDOnVK0fsAAqFzZ0aNT/zJ00WJGHHoZjHjIDU9l68+Kd8aRO3sOALw9vRi7zoPzYZUd2+ge4PRpWUGZMeNuzsCECZIz4OvrmDGpJKU6ABuG0QNJzmpqmqb5oNeZpjkNmAZQt27dB75OqUzt4EHZg5szR6pUnTolZ3YPHdLkFyezZP05Xl86nIsFZ4JhggFWw0rLeiXvBF9wcKN7kKSqdevkCNEvv8h11K6dJFU1bqwrKC4gVQHYMIyWSNLV06ZpRqTvkJTKRLZsgU8+kT6pHh7SKWbgwLu9UTX4Oo1ffr/EwPmjOVdgMhSAmp4dOWKuIjY+Fg+rB41KNXL0EMXVq3dzBk6dgkKF4OOPpRBLsWKOHp1KgeQcQ1oANALyGYYRDAxDsp49gfWGfMr62zTNARk4TqVcR1CQBNuCBeHmTanT/OWXUpc5Xz5Hj079S0BQADP+WMuvAWcJ8V0ChaJ53KMXs/t8QsVCJQgICsD/jD+NSjVyfE/e7dtlb3fhQimg8dRTMHq0NN7QCmguyXjI6nG6q1u3rrlz5067vZ9SdmOzwYYNcoNcvRqGDIExY+Rx0JmuE5q2YSMDtjyLacSCAWWNZvzUcxK1SlRw9NDuioyUgOvnJ/kCOXLcTaqqVs3Ro1PJYBjGLtM0E21FppWwlEqr776TIhknT8oZy3fflTOWoIE3g6SlleDufVH09pvK3jwfgWcsABbDSp/GTZwn+J44IUvMM2fKOd7KlWHiRAm+Pj6OHp1KJxqAlUop05SkqoQZSEAAFC4sZyxfeEHPWGaw5LQSTCxAl/cqSO8Js9juPRwKB1HIrEOo9QBxtjjn2OONj5dkqkmTJLnKzQ2ef15mu08/rUlVmZAuQSuVXLduwfz5ks0cGCj9d6tVk/04Dbrp6mEz3AZjNiXaSKFoLm+2vt/kPwE6JtSLsKAjRJYbD3mPU8T2OBOfH8nzNZs6xx7vlStyfCihEEuRIpJQ9eqr8nfl0nQJWqm0uHJF6jLPng3Xr0vQ9fODkiXleQ2+6SqpGW5SrQQTev3eijhFWLg/cTn2wuMncbtZiZnPrqLro224nTxK/eL1HRN4TRP+/ltmuwmFWBo3lkIs7dqBu7v9x6TsTgOwUomJiYGQEGlA7uEBP/4IrVvLEaIGDXQ5MA2S2r9NCKD/Fhkbz9h1R2lfq2iSrQTPnTO5fHEnMdWHQ14bmAY5b3Ujj1tHuj32XMb+ckm5dUtKQ/r5ySpKzpwy0x00CKpUcezYlN1pAFbq34KCYNo0aflXuLD04PX1lWDs7UR1fl1UcvZvk5rhPqiVYJ/alXnprb8JjhgFNTbd0yjB6mFQNEf2FI81tYle/3HsmGxd/PCDrKJUry5fd+0qmc0qS9IArBTIcuCYMXKEyDShVSuZ7SbQ4JsukprdAknOcO9vJZjPzYfs5914dU5/bOVX454tH94RL3DTe/WdRgm+Rs0U1WlOzgeFJMXFwc8/yzLz+vWyrPziizLb1VUUhQZglZX984/s3+bMCUePSoP7hCNEpUo5enSZUlKzW3jwDPffAbR9raJ4uZ9j6MKf2RWyD1upDbjH+/Jmjc/5rNVgNhy8ztC1TxMUsZPi2eoy/NkOKZq9JueDwgNduiSN7qdOlRWVYsUkh6BPH6lapdRtGoBV1mKacmxo8mRJfhk1Ct56C7p0gc6dNaEqgyU1u4X/znDvX/4NDYV+Xy9niaUjuMVDSWhXshs/dBpPbu/ct39GDtrXej3V40zOB4V7mCZs3Sp7u0uWQGwsNG8uDRHatJEjRUrdR68KlTWYpuzrTpoE+/bJrLdvX3j2WXleS/nZRXJmt5B4o4MbN+Dzby4xPnAUMY9Mki5FtxslPF628p3gmx6S80EBkO5W8+ZJ4N23T/IFXntNultVdHBrQuX0tEyPytyCguRPw5A+qYYhS4Pnz0tlIW0BaFftaxVldIfqFM3ljYGc3R3dofpDl3XDw2HY6FAKvvwhY2PKEFtzEk8VeRYvdy+shjVDimgMaVERb3frPY/d80Hh8GF44w0oWlSCrcUiH/BCQuCbbzT4qmTRQhwq84mOhqVLZZk5IECaIRQpItmnPj6a/OIiNh8PYMxP6/g98CLR5RaC13WeKdKZ7zp8RoW8FTK8iMb9WdDvNilDu6BdMtvdtElWTV56SZKq6tXT60ol6mGFODQAq8zj0iWZfcycKcUzypWT2UnfvtqU3IVERcGQSf5MvP4MWKRRwiO5GjCn0yRqFKph/wFduCCz24SVk5Il5brq3RsKFLD/eJRL0UpYKvOKj5ds5oIFpXjGN9/cLZjRtKk2Q3AhMTEwbXosnyydTdhjb4NXQqMEC51rt7Zv8DVN+OMPme0uWyZHilq0kHKRrVrd7eesVBpoAFau6eJFOeoxbZpUEPr1VyheXB7PnX7JOCrjxcbCD7NsfDh/EVerD4WnTlAyWxUuRp+0f6OEGzdg7lwJvAcPyrX0xhvyga5cOfuMQWUZGoCVawkIgK+/hhUrZFbSvPnd1n+gwdeFxMXB3Lkm7/+wmktVP4ZG+ymVrTrjn1vFcxXb8Hfw3/ZrlHDggATdH3+UrK86dWQro1MnyJYtY99bZVkagJXzu3YNsmeXM7pbtkgCzODBEnjLl3f06FQKxcfD57MD+HrzdG74BECTwxT2LMe41vPpXK0TFkO2DTK8UUJMDCxfLoH3jz/k+urUSY4RPfqoJlWpDKdJWMo5mSZs3y6ZzIsWyXLzK69IMXuLRUtDuiCbTbZTX581jYt1BoLFBhi81+A9RjQejrvVTh2AgoPv1vu+eFEabgwcCL16Qb589hmDyjI0CUu5jvh4WfpL6LmbIwf07Am1a8vz2VNWUF+lXLo2IUA+S61aBe9+tY9jxT6GR1ffec5qWPD19Mn44GuasHmzFGJZuVI+DTz7rBwhatlSk6qUQ2gAVs7h0iXJZLZYpHyfxSJB+JVXpGqVsot0aUJwm2nC2rXw7hfHOJh/GDRdSDaLL6/U6M/c/XOIiY/J+ASr69dhzhxZZj5yBPLmhbfflu2LMmUy7n2VSgYNwMpxoqOlbu7kyVLGLyREgu3mzXKj1D24dJfWXrzJYZowYXkAX6xcwYXYQ9BoLR5WT96q/yHvPvkOub1z06tWj4xNsNq3T4Lu3LmybfH44zB7thTO8PJK//dTKhU0ACv7CwmB8eOlN+o//8jxjqFD7wZc3YfLEOnRizcp/v7w2vifOVS9PZSOAwNerPwSE1tNoGCOgndelyEJVjExUgFt0iRpjODlJU02XntNspqVcjIagJV9xMXJGcs8eSTofv01tGsnyS9NmmjBDDtIj168D7J1K3wwPJQtceOgwViwxgHSKKF24Zr3BN90d+7c3aSqy5flA91XX0nuQJ48Gfe+SqWRBmCVsUJCJIP5++8l0M6ZAzVqSHm//PkdPbosJb168f7btm3w0WfhbAyfgPHkWPAMo3qehhy4tg3TjMPEDSOmavr+IiBJVBs3ymx39e2krjZtJKmqeXP9QKdcggZglTH8/WWZefVquVm2bClnLBNo8LW79OjFm2DWhgBGLdrA8aCrGDUWQLbLtCzbhubF3mDaxngKxh8gyrIfL1t1Zvt7UClPSJoyqe8IDZW9XD8/OH5crqP33pOkqpIl0/7zlbIjDcAq/Vy5IslTFgusWSPrkkOGwKuvasapE0hLL94E+/bBa19s4c+yTaFoLBSDGgXqMPm5ldQrVo8GYzYRGRuDJ5XxtEmrx0hbypK4EhUYKLPd+fMhMhKeeAKGDYMXX5QCGkqlVXy8LOk88YTd3lLXaVTaJBStf/ll6Y26caM8/sknUvBg9GgNvk4iNb14Exw6BB1fslHjlQVsLd4erNKlyGJYeKnaC9QrVg9IexLXPaKipDRk/fpyDnzBAujaVYLx1q1yRE2Dr0oPp09LQZYGDWRlxU50BqxSJypKEl+mTJHm5L6+klBVtqw8r+3/nNLDZreJOXYMPhtuMn/naizNPoaq+ynpU4bzt8KJt8X/5xxvapO47nHmjFxXM2ZIwl6FCvDtt9CjB+TKlfyfo9SDxMfLIfVr16B7d9m+aNRIEkNLlbLbMDQAq+QzTSmYUagQuLnBF19IByItWp/pnDoFI0bA7C0bMZp+CF22U9K3HKOaLcAjpgGf/rqCoKidFM9Wl0v/lIDi8n0pTeK6w2aDdetkb/fnn+VIWtu2coRIs+RVegkOlg9206fL36tXh27d5PqaM8fuw9EArJIWHi7Lf5MnSwA+e1YC8J49mkyVyZw9CyNHwozAaZgNvsDsdooiOYrzaePv6VGjBz/vu8wHy/cTGVsGX8pw4wb3nCVObhLXHVevynnwyZMl6hcsCB99BP36yYc7pdLLV1/Bu+/KROKZZyRJ9LnnHFrwR5sxqAc7cUIa3P/4I9y8KZ8WBw6E3r117y2TCQmRwPv9qr3EtXwNim8FwN3izm/dfruzzNxgzKZEl5iL5vJm6/tNkv+GO3dKUtXChbKd0bChHCHq0AE8PNLjV1JZXVCQzHY7doSqVeHvvyU5tE8f2e+1E23GoJIvMlJuiLlzS2LCjBmyvNy/vyTDaHnITOXiRRgzBvwWHyPuyWGYry7E0+pFTLyBiYnNtBEQFHAnAKcpySoyUjpb+fnBjh3SWKNnT/lQ98gj6fdLqawrLk72dqdOlT9NUyrrVa0K9erJf05EA7ASR45IUtWsWTLDHTcOmjaVqVHevI4enUpnV67Al1/Cd3POEVN/OPSbhbe7F/+r/xFPlXyK9gvbJ9osIVVJVqdOyRLzzJmS9FK5Mnz3nSS/+PhkwG+nsiSbTVbpjhyRPJX334e+fe06202pJAOwYRgzgTbAZdM0q91+LA+wCCgFnAFeMk0zNOOGqTLMsmVyM/T3l33dDh2gfXt5zmLR4JvJXLsGb34VwMJda4j1OY5l4ErcrDDo0f/jgyc/uFMycmP3jYk2S0h2klV8PPz6qywz//qrXEvPPy/LzI0a6UqKSruE2e769bKfa7HI9VWsmFRFc7dTf+k0SHIP2DCMp4BwYM6/AvCXwDXTNMcYhvE+kNs0zfeSejPdA3YSISFyZhdkFrJliyS99O4tSTDKIZLThze1vXrDwmQ7f+yydUQ+3wYs0ijhuQrPMbHVREr4lkifcf7zj2xbTJkix4kKF5Zr69VX715zSqVFwt7ujBmSyVyoEOzeLdeaE0rTHrBpmn8YhlHqvofbAY1u/3024A8kGYCVA8XFSQLClCnw229ywdasKb13fXz0mIeDJadTUWp69d68KZODcRPCuV5pPJYOIyT4Io0S6hern6Lgm/Be97yfaUqCi58fLF4sbSYbNZI17vbtXWImolzExo2SwXx/JrOLXmOp3QMuaJrmBQDTNC8YhlHgQS80DKMf0A+gRImU/UNX6SBh6jN9Opw/L8szw4bd/bSohQ2cQnI6FaWkV++tWzBxInz5dRTXykzBo98ocL/CEyWeZEfIDuJscf/Z302xiAg5nubnJx/ocuaUPbdBg6BKldT/XKUSBAdL7kCxYrJC98QTckytV68U7+2mdvUoI2V4EpZpmtOAaSBL0Bn9fgrZfzt/Xs5RWq1SReiJJyQRplUr2etVTiU52cXJeU1kJLw/MYAZmzdyKyYcrz7zwDOYhqWb8nmTz6lXrB4BQQGJ7u8m2/Hjci398IN8wKtWTb5+5RUJwkqlRXy8FGWZOlVW7Ww2OYXRuzd4e8Pw4Sn+kalZPbKH1N6JLxmGUfj27LcwcDk9B6VS6cIF+bQ4bZqUgty7V26IZ8/qTNfJJSe7+GGviYqSjo/DZmwl9Lkm8FgMGFA2f1UmPDubJqXvntGtX7x+ygNvfLzcDP38ZAvDzQ1eeEEqVT35pCZVqfTTowfMmyf5KO+9J6sqaawnn5LVI3tK7cbfKqDH7b/3AFamz3BUqmzfLl1hSpSAjz+G8uVh6FDZJwENvi5gSIuKeLtb73ns/uzixF7jZXGj+s1alCtv8sakVYS3egHcYu40Sni5+sv3BN8Uu3QJRo2SG2D79nDwoMxAgoKkiEbDhhp8Veol1GTu0EGSQ0GS9n76Cc6du3vtpVG6NglJR8k5hrQASbjKZxhGMDAMGAMsNgyjD3AO6JiRg1SJuHJFKgb5+krFfH9/ePNNuXjLl3f06FQKJaeE479fE3I1CrdTpbkWUJ6pXr+To+OH4LudgjmLcTnC406jhMalGqd8MKYJf/0ls92ffoLYWDkT/s03Up9ZtzBUWiWs1n3/vazQ5c8v53eLFoWnnkr3t0uXJiEZQEtRuhLThN9/l72RpUulduCQIRATI/skXl6OHqHKYPHxsjo3fDicjA4gZ/uPuJlvM8V9ijPs6WH0qNmDHSE7UrfHGx4u/Xb9/GT7wsfnbqWqSpUy7HdSWUxYmBwdio6WRhv9+8vqSgaWIL1/DxhkhSm57TjTQktRujrTlONCkyfD0aOypDxokMxGQGvnZgHx8XLC57PP4KjnHDzbjQafI3hnL8DnDcfTv05/PN2kPneK93iPHJFra9YsuHFDykJOnSpJVdmzZ8wvpLKOS5ckYe/kSZnx5sol19uTT9pttS7FTULsRAOwszJNuTFWrix7bGvXQp48cpPs2FFb/2URNpsUK/v0Uzh48RjZXxoEBTcSjTRKmN9hPk3LNE35D46Lg1WrpFLVpk1yjrJjR0mq0prfKq1sNrmupk6FFSvkemvcWFbrPDzkGJGdpbQXtj1oAHY2YWHSfWjaNDh0SKoJFS8ud2ENulmGaUp8HDYM9p45i2+74VjKzCbGYsGw3W2UsD1k+z0BOMmzjhcuyJnwqVMl6aVECUl06dMHCjzwOL9SKTN1qqzS5ckDb7whuSkVk+gJnQVpAHYW587JNGfhQjnMWbeuXMR58sjzGnyzBNOUxY6hQ2HXkUv4PjcKt+enEGmFN+q+TrOyzei4uGOijRIeeNbRNGkffkr2dpculdlI8+Yy+23TRs6KK5VapilJoFOnyl5u585yKsPHR46qaW7KA2kAdqQbNySbuWxZySxdvhy6dZOkhNq1HT06ZUemCRs2wP++DuBg1FqyVz6DR9ulhBNN71q9+eSpTyjuKw3qH9Qo4f6zjtmjI3g+0J/q3w+CS6dl7+2NN2DAAM2UV2l39SrMni2B99gxub4aNpTn8ueXHAL1UBqAHWHnTllinj9f+lNu2ABFikhzVm10n+X8/jt88glsCd4IXZ8Fayy3gOalmzOp1STK5703WD4oySrhTGO5f87RLfBnOhzYRM6YSPYXLCuF6zt31pUUlX5atZIaBAnlITt2lEpVKtk0ANvTsmVydGj3brkRdu4ss90EGnyzlL/+ksC76Y8ofJpMwaPHx8QQC0ijhMalGv8n+D5QbCyvBG2n9Z/LqX9uP9FWd9ZUbsiPtVpzpXINtvZORaKWUgnCwmDOHDkDt369LC+PGyez3urVHT06l6UBOKPt2QMVKkjAPXNGihpMnAhdu0oRDZXlbN8ue7zr1seSs+Escn0ynDAzmLqF67L/8v6UNUo4f15WU6ZN4/MLFwj2LciYp3uy+JHmXMvmK2cdW+oZXpUKpgnbtskS86JFkpvy6KNyzfn43F1uVqmmhTgywq1bkkw1bZrcbX/4QQoaxMVJwose8ciSAgPh9S8C2Hp+M97ZY8lWfy5XzRPUK1aPkU1G0qR0k+Q1SkhIevHzk7wBmw1atoRBg1hR6BHGbjjhVGcdlYs6fFi6WuXIIfu5/ftDrVqOHpXLeVghDg3A6SkqCt5+G+bOlQSrKlXkou3WDXLndvTolIPs3y/HiZbv+At6NQZrDABlc5fl25bf0rp8a4zkfCi7fl2OqPn5yc0xTx7pEDNggCTyKZUWu3bJbNdikb7hAEuWQIsW2uUqDbQSVkaKiJCpTYMGsocbGCgVqvr3l8d0tptlHT4sJ8sWL4ZsVTfi07cnN5Dga8FC71q9aVOhTdI/aP9+OTI0d66srtStK6sqnTpp0otKm/BwWa2bOlWSQ729pUiGacq968UX7TIMZ+zVaw8agFPrwAFZYp4zR/Z1L1yQfZE//5RPkCrLSbiJnD1lEL2jMlf3FsSr3N+UHPoRZy2bKZCtAJFR7thMW9KNEmJiJGlv0iS5pjw9oUsXKW7w6KP2+6VU5pQQYEeOhDFjpKfzd99Jboqdu6c5a69ee9AAnFLbtsFbb0kKq4eHfELs3//uEo0G3yxpRWAIb884waXfy3PrYFEovAf3ga8QWWATkdkLMP52vebdF3Y/fI83KEg+2H3/vdTQLVMGxo6VWUnevPb/xVTmERkpyzFTp8qeSIsW8oGuTRs5SuSg1Tpn7dVrDxqAk+PgQUmeqlRJitNfvQpffQXdu0O+fI4enXKwc+dgQH+4tKshVFmE5c3PsPkeJd7MTglrHw69MZ7sHtLUINEzvAl1cydNkvqTpgmtW8vNsUUL/VCn0ubwYQm6s2fLcaKKFWXVDqTMbfHiDh2es/bqtQcNwA8SGSkJCFOnwtatkgU4d64s1Rw+rHu7ipAQKaP8/fcQmzMa66vNiS+0CRuAaSVfzHtYomrfCb7/ERYmzTUmT5ZKQvnySXvJ/v2hdGk7/iYq00lYYrbZ4JlnZDXlhRfk2nr6aae6fzlrr1570ACcmBEjpPl4aKiU7Bs7Fnr0uPu8E128yv4uXpRtsylTIM7rIpUGj+JQ9inEc3sZzQBMkxjLCcr5NPjvD9izR2a78+bJB7169SS7+cUXtW6uSpuTJ2ULY9062LFDulwtXCj3MSdttjGkRcVEe/UOaZH5mzdoAAa5Ca5aJTdAq1XO67ZoIR08GjXSgKsAKds9dqzUUYm2hFJt0FiO5x3PkfhompXoxKGTZTlv/QLTjMPADV+j5t2bSHS0rKhMmgQBAZJt+vLLssysdb9VWsTGwurV8olw/Xq5h7VtC9euQcGCchrDiTlrr157yNrngBP2RubMkdnu+vXQrJmjR6WczLVrsuU/fjzcKrCRQi9+wQ2fv4iMj6BL9S581ugzyuUpx4rAEIauXUZQxE6KZ6vL8Gc70D5PnFxj06dLBC9fHgYOlMIsejZcpUXCMvPGjXLfKlYMXn1VWksWzfzBy1VoIY77XbgAL70kxzvc3WVvRGe76j5hYfDtt7IbcSMiivID3+d43vEAWAwLs9vNpmuNrv/9RptNPsz5+cGaNfJY27Yy223aVJOqVOrFx8vy8pQpkhT65Zdyvf32mwRhN13UdDZaiAOkuf3p05JdWqCAHCFK2NvNn9/Ro1NO5OZNmDBBas2H3YilVu9ZnC8/nOORwXdeY2AQdCPo3m8MDZUCGZMnw4kTcp198IEkvjg401S5uIsXYeZM2d89e/bepWWLRUqRKpeTuQNwQibztGky2y1RQoKw1SrLNkr9y61bskX75Zdw9ZqNWt0WcqXaUAIjTlI/b30+rv4e765/l5j4mHubJezaJbPdBQvkmmvQAIYPl5UVDw+H/k7KhSUsMQO8/74cI2rSRCYO7drptZUJZN4A/OOPMHjwfzOZdflP3ScyUlb0RvwQQKjvZiq0tZKr6lwCww9QI2cNprRbQ6vyrTAMgzqF60ghjSL1qb/5OEx6UxpuZM8u58IHDoQaNRJ9n6xabk+l0NWrcjwtoQtRrVrSt/LDD6WzWjrT69JxMk8AjoyEn36C+vUl4BYvLuff+vfXvV2VqOhoOcM7ahRcsN5ulGCJ4ZgB+eKKsPCFhXSs2hGLcfdDW/24QtRfEwYzXpQbZaVKsl7dvftD20tm5XJ7KhlMU6rrTZki97HoaFlJiY6W5zOo2YZel47l+tPBgwdlplukiMxwFyyQxxs1kvNvjRtr8FX3iImRyUW5cvD66+BZ2R9rj47SpcgATANuNcMz9kkJvvHx8MsvXGzYFFvZssSNHYd/oSr8OWWR5Ba8/nqSvZ0fVm5PZWEJSbAREfDss7ByJfTtC/v2ybZZvXoZ+vZ6XTqW686ATVNmuBs2yF5IQibz0087emTKScXFyYmzESPgzBmo3nwPRYd8zLbQnzHMHGBaARMDN9ziqjNt2XbabzghSVWnT2PNnpvv6ndiQY2WXPTJh3ewldF7zidrppCVy+2pRAQGynW1b5+cC8+eHdaule2LHDnsNgy9Lh3LdQOwYcgSzbPPak1m9dB9rPh4mD8fPvALIMTNnzLVS9PwreVsubaYXJG5yBXbg5xxzxFrOU2UsZ/ql30Z/Pda2hzeAvGx8NRTDK33CguK1iHW6n7nPVNSMD4rl9tTt0VEyJ7ulCmSN+DtDZ07y+PZs6eqYEZa92/1unQs1w3AIM1WVZb3oH0smw1ijhfls8/gSHgARq8mYI3mFCbnr3vxccOPefuJt2n97W6uXgnlucNBdAv8i0cuniDcw5uf6z5Lh+kjoVo1fnz/ZxI7MZ/cR5q6uQAAIABJREFUmUJWLreX5SVkM69YAb17Q+XKUtWlW7c0FWNJj/1bvS4dy7UDsFL8dx/LNOHqgfx0+96XiEtQofZFirw8mPO2qDuvaV3mVUY0GQHHjzP74CLyLZlPrqhwjuUtwSfNB7C2ZjM+7lIPqsmNLK0zhaxcbi9Lio6G5ctlttuqFbz7rmyTFS0KTz2VLnkp6dHGT69Lx9IArFxewizUNCHyZAHCtlQg9rIv1sLnaDt+JOtujCc6LoqEnEMDN2y/3eLSxEYUDPidcm5uhDRuyUelmvBL7vIUyZ2Nj++7CaXHTKF9raJ6Y8vsTp+WugMzZkjp0dKl7/Zx9vRM1xyV9Nq/1evScTQAK4dL6z5WYV9vTgbmIOzPCsRcyIU1/0W8+35ITLHZrA6NIK+lCXmjO5EjKoQyV1cxYPtp2h2dyRWfvPDZZ9C3L0WLFGHSQ95DZwrqgf5dMGPQICkL+dxzMGCAJIpmUO0B3b91fRqAlUOlZR/LNKWPfeiiBlze64lReR3uL40kLtcuIo0IHivUguntvuTTd36ja+B8Wh39E8/4OLaWfIQB7QeysdzjHB/aLtlj1ZmCuseFC9JkY+ZM+P13qbT31VeQM6ddSo/q/q3r0wCsHCq1+1h//CHFgf74A4oWt1Dz48HscZtALAAGr1X+nInXC8CzPVi6Zw83PLIxv+azzK3ZipP55OZYVGcKKqVMEzZvliNEK1bI2bZmzeDGDXm+ShW7DUVXZVxfmgKwYRj/A/oCJrAf6GWaZtTDv0upu1K6jxUQAEOHyvHvQkXi6T5uIX+6DWNP2Mk7r7GaUPT7EbAhGqpXZ8+HY+gTV4Grxt3auTpTUCmSsMx86ZL0CvfxgTfflNoD5cs7bFi6KuPaUr05YRhGUeANoK5pmtUAK9A5vQamsoYH7Vfd//iOHZJM+sQTsHefSa8vVpD3w5rMCe9KDo9sjC3cHe94C9Z48IgzaVT8KdiyBfbupebI9/ik82MUzeWNgcx8R3eorjcu9XCmCdu2Se/mNm3ksUKFpNVkSIjUl3dg8FWuL61L0G6At2EYsUA24Hzah6SykqT2sfbsgWHDYNXuALyqbKbNiOxczD+fHy5up7xHGRZaOtFx9FYsQXNoUKcA/u1q0KjtYOrXaH3P++hMQSVbeLhUbpkyRSpW5cgBXbtKRRerVcrcqju0mUPqpToAm6YZYhjGOOAcEAn8Zprmb/e/zjCMfkA/gBIlSqT27VQm9aB9rHLuRXnxRVi6FLwqbIXeTYgyYlgTD7kv+zLjwmN0H7Ubt5hT0Lw5fDeR+q1bU18bkqvUSlhmnjFDlpcfeUT2el95RRKr1H9oM4e0MUwzsfo+yfhGw8gNLAU6AWHAT8AS0zTnPuh76tata+7cuTNV76eyhiNH5GTQokUy8ajX7U/883Yh1hoMgGGDz/zh/d0+uPftLe3/MqBFm8oioqPlU97kyVKlqlcvaWF6+LB0VtNGLg/VYMymRI9CFc3lzdb3mzhgRM7HMIxdpmnWTey5tBxQawacNk3zimmascAy4Ik0/DyVhZ04ISW9q1aF1avh1Q+O0nRyJ9YXaIjFdgG3eLDGg2e8wdECL9JyyHz45hsNvip1Tp+WJvfFi8sM98IF8PKS53LnlmQDDb5J0mYOaZOW9bpzQD3DMLIhS9BNAZ3eqv942B7R6dPw+ecwa0MAlrL+tBlSkew1VjH92By8j1r4+E/4v20ws05t5tYoyLUcTfizbGWMCAf/Usr1/Ltgxksvwe7d0LatrKI0a5ZhBTMyMy0GkjZp2QPeZhjGEmA3EAcEAtPSa2Aqc3jQHtE/F63sXFmIGTPAKCGNEuKNaFZh4n4I3tgBH5wpzNJSLWnVuxFXs+cCwPP2jon+A1fJdvGiFMyYP1/Osfn6SoJVwYJQrJijR+fStBhI2qQpY8U0zWHAsHQai8qE7i+0EXfTk5C/y9FvdAHcLNC931UOlRrA37fk+Lhhwv8uleaL18ZDq1YU3neRiGX7Qf+Bq5QwTalO5ecnTRHi4qBpU6nP7OsLdeo4eoSZghYDSRtNGVUZKmEvKP6WB9f/Lkv4npKYNoOCj+yk1/PD8Ytey41bNqwmgIGHmyftP5kHxesD+g9cpdLRo9C4seznvvEG9O+v+QIZRI/4pZ4GYJWh8rv7cHR9EW7uLokZZ6VOhV8oX/k9NlQ4yBgT2l3zZUTNtwlv8iT+F/+mUalG1L8dfBPoP3CVpF27JJMZZLm5UiXJ5mvaVBrfK+WENACrDBEaKnXp93/TgPiIWJrXfo+oatM5VPAGu7LDkzcKMK7+aB5/pved76lfvrEDR6xcTmSknFfz85NSadmySdWqhGSrhOpVSjkpDcAqzf6d5VzAKwdlLtVi7QIfclwPZlbVKewt9g1j6kVgGrLHO6DMB0zuNsrRw1aubvRoGDECKleGCRPkHJuvr6NHpVSyaQBWaZKQ5XzrFtzcVZZz20qTPXory/KPJ6zKGoY+beNggbuvt1islCilVYVUCsXFwapVssz8v/9JYfB+/aBJk/9v777DoyrTPo5/nyQkdAKKIt1FpClKsSBIUVdRkSZgF8GVtaxr2VfFVayroKxYMDSliW1FQBFUsLEuTaVXqYIQOkoVSDLzvH/coZpAJCFnZvL7XFcuJnPOJM/hzJw7Tzn3bUXudc+uRCEFYMmVXp8sY8PkKjD9JG7Z+z53FerLmmrLeeiyOGadFubMkqfzTP0u9Jzck7RQGonxiTSv2jzoZku0SE21Od1Bg2DdOkucsWuXbatYUbcRSVRTAJbjsmcPDBwI+3oWpdfeV6heaRgf1E2jXdVElpSFhHAZBrd+gVvPuZWEuAQu+9NlTFo1KctFViJZ8t56tytXWgnA/v2t56t83xIj9E6WP2TfPhjSfx8Lnh7FDdv6cT9TGFQvnsuuCeHjAJ9GifQ21Cl2J13rtTzwukaVGinwytH9+isMHw6jR1vB58RE6/lWqQLVqgXdOpE8pwAsOZKWBiP7/Mwvzw+k0843uYtNzK5dmSva1WZioUVwoKZHHIXjSvNIy7ODbK5Ek5kzbSXze+/Z0EqjRpabuUoVm+MViVEKwHJUGWlhvnr0S1z/fly/5xMA5rS4hO4d4nhry5cUSSjCJeU6MWnNR4R9BnGuEA82a6f7diVnvvsOLrzQbiG65RbLy3zuuUG3SiRfKAAXcNkVSght+ZXZ9w2j9Af9KVluGeMaFmVs0iX83D6Zz7d8TNwvcdx3wX10b9KdU4qdwrQ10zTHK8e2bJnlYS5e3OpOnn8+DBkC7dvrFiIpcI67HvDxUD3gyHJkoQSA+pt/4pF506k7czRF/B4GVa3N3bcsJRSXkbmH4/KqNzK4fS8qltQKVMmBjAwYN86Gmb/4whZRde1qq/hEYtzR6gGrB1yA7S+UkJSRxpU/TubmqRNp+OsCdlOUcck3M6xxE76q18eCrwM8lMhoy64NXRV8JeceegheecVuGXr2WfjLX6BcuaBbJRI4BeACLG7VTzwy+zM6zfmSk9K2s4QzebDoC4xpVo97HpvP51/cS9jtAB8HHhwJFAtfpGLbkj3v4dtvrbf70EPQsCHccQc0bw5XX61biEQOoU9DQRMOw4QJ+JQU/jv+UzyOj2nDwGJdmdH8FHzdEexO6s1DX26hVFwDiuy5EVyYvXHzKRw+m6RwLdXild/bsQNGjLDAu2iRVSG69loLwLVr25eIHEYBuKDYuhWGDMEPGIBbuZJfCp1KPx5jSIlb2XrFFEI1XmBf/CJCcb9Qs3RD+rceybZfqx+YI04K1wJUi1eyEApBnTqwdi2cdx4MHQrXXacqRCLHoAAc6374wXol778Pe/cyt1RTevI835dtx8OPF+Ki5H/y7tIXsBt5HZ1r92Boh6dxzkFV+xGqxSuH2bcPPvwQPvvMer3x8dC7N5xxhvV4RSRHFIBj0f4ybSkpMGMGvyUU4e3C19GXf5BaqDZPvRrHrX/+kqcnP8YPS3848LJ4F0eNckUs+GZSLV45YPVqu4Vo8GDYvNkC7vr1UL48XH990K0TiTpxQTdA8tCKFbbwpWJF6NKFret28ODJPSmXsYG7wgNJbRFP8Tte4424xrT64HI27t7IP5v8kyIJRYh38SqUINn79ls4/XR48UW46CKYOBGWLLHgKyLHRT3gaBcK2VBgv37w+ecQF8e2Fu3o/ds9PD+1GXGF0yl6zSiKnvUeuxOW82v8YjZsLU3fK/tyR/07SEpIotWZrZREQw63davN5ZYqZauYL7zQEmd07gyVKwfdOpGYoEQc0WrzZssgNGAArFoFp53GprbdeHz1HbzxaQWSk4GzlxB34TB+KfoiuDB4KB66ijLpXVnd69qgj0AijfeHrxnYtw9uvBHeeSfololELSXiiBXeW+7cfv1sjjctDZo3J/X+F3lkalveHVCI4sXhiSfg2ttX8+fh97MpNIGDlRLiSPAnUzG5dJBHIZHqH/+Al1+2NJFdusDdd8PZKqohcqJoDjga/PabLXxp0MAqxXz0EXTrxupPF9K58jdUfrAjH40vRPfuMH3BBn654F4aDq/Or3xDCX8xjkTwcTgSKOXO1W1EYpYts6C7apV9366dLdxLTbXauwq+IieUesCRbOlSuxAOGwbbtsFZZ0H//qxuchPPvlKCYddAfJVpXPTIJO5pX585276h4YjXSA+n0/XcrvRo1oMZKxxPfDaaNb/NoFLRhjxzZXutai7IQiEYP95GUSZMsMxU9etD1apw8cX2JSL5QnPAkSYjwy6QKSkHE9d36AB3382aKk14vqdj8GCIi4MGN37FtMpX49lnuZqBm86+iaeaP8UZZc4I9jgk8qSlWcKM5ctt9fJf/2oLrE47LeiWicQszQFHg40bbZh5wABYs+awxPXrfTl69rTiMd5bLvuzWq3k4e8ePRh8PZQOt6HD6S9wRhn1cAV7s0ybBl99BT16QGIi3H47nHkmtG6tvMwiAdMccJC8hylTbKVppUrw2GN2cRw9Gn76iWGNbqd8xzQqVA7RNyVM06t3s2hJOud0HcgDP1zA7rgfAAfe4UgkKaMJvScsCfqoJGi7d8Mbb9jQcuPG8NJLsGmTbeve3WrvKviKBE6fwiDs2gXvvmvzcHPn2r2Wd91lXzVrsnUr3HHHTj56pyw+I55iddZS8qIlLDx1Ik0+/ICNv60mKVSLUzMexFHosEIJqlRUwH37rfVut2+HunVt2OTGG21ls4hEFAXg/PTjjwcXVe3YAeecA4MG2QWyWDF+/RX69LDSqbt2F6dorXUUaTGGvaVHsSluMaG4jRTbcwbjbhjHc6MKs277XoADhRIAVSoqaDIy4JNPrPBBy5b2nmrTBrp1s4xVh6QVFZHIogB8omVkwNixtqjq66+hUCHo2BHuucduKXKOHTvglWegTx/ruHTsCJNL/pe9p77D1kKDwXnwjlLpN5Oc0Ymrz7ya9JapByoV7adKRQXIxo02zDxwoFUhuuoqC8ClSsHw4UG3TkRyQAH4RNmw4eAFMjXV5nife85WUJ1yCgDvTVlH92f3sPa/lQjvTeSCFnsY9EoRdiZPYcKw+9nh5x7MoYHDEUeF5GIAB24lUqWiAujpp+29lJ4Ol10GfftCq1ZBt0pE/iAF4LzkPfzvfza3O2qU9X4vv9y+v/pqK9uG5dW4+7HtvD3oZEK/JVKk2kZKNVnKxkpz6TplDDM3fU1yUllK/9aObXHj8T4jyyQaqlRUQOzebekgO3SAMmWgVi3LUnXXXVBDIx4i0SpXAdg5lwy8CZyF9dW6eu+n5UXDosrOnfD22xZoFyyA5GT4+9/hzjuhevUDu+3dax3inj1h48ZSFK66mVKXjSb91M/41f3EvoR5rN1Ugl6X9uJv5/+NLxZu44nPWiiJRkG1ZMnBNQPbt9v0RZcu0KmTfYlIVMttD/hV4HPvfQfnXCJQNA/aFD0WLrQL5FtvWRCuX9/u5b3+eih68L9i3z57+rnnYN06aNEC3GVTyag8ni2JvYAwAMUyLqNM+l94pMl1ALStV4y29e4N4sgkSPv2wTXXWCKWI9cMiEjMOO4A7JwrCTQFbgPw3qcBaXnTrAiWnm65mPv1g0mTICkJrrvOhgTPP/+wVafp6dZ5+de/4Oef7ZbMt9+GWudtoN7LvdkSGguEMxNpxFHIl6dSctmADkwCtWkTTJ5s9+gmJVl2qn/9y9YMnHpq0K0TkRMgNz3gPwGbgaHOuXOAmcB93vvdedKySJOaaouqBg2C9estd+4LL0DXrnDyyYftmpFhgfaZZ+Cnn+CCC+ylDZr8Qu+pL3L1q6+R5tMp7i9gt5uZ7RyvxLj9mapSUmDkSHtu/Xo46SStZBYpAHITgBOA+sC93vvvnHOvAt2BHofu5JzrBnQDqBxthby9t15uv34wZgyEw3DllRZNW7Y8sKhqv1DIyqg+/bQVmqlfH+7uOY2dp0xg5M5UOr72ATv37eSmujfxVLOnmL+6sAolFFQzZlge5jlzoGRJWy9w990WfEWkQMhNAF4LrPXef5f5/YdYAD6M934QMAisGEMufl/+2b4dRoywwLt4sa08ffBBS15frdrvdg+H4cMP4amnbPezz7Z4XaruJK5453LSf0wHoGmVpqRclcJZp5wFQLUyaI63IFm+HPbssTdI5q1o9O8PN9+sTFUiBdBx54L23m8A1jjn9o+ZXgosypNWBWXePOuJVKgA995rF8WhQy3RwYsv/i74em+B9txzbRoY4IMP4IeZ6WysOJC2/2lDetiCb7yLp2W1lgeCrxQQoRCMG2cjJtWrwyOP2POVK8Ps2fZ+U/AVKZByW4zhXuAd59w84Fzg+dw3KZ+lpcF771kd1HPOsbm3Tp3g++/t67bbLM3fIby3a2qDBrZmZu9em/OdMzfEvhpvU7t/Te4cfydVkquQFJ9EvIsnMT6R5lWbB3KIEpChQ+GMM2xF87x5NkTy5ptBt0pEIkSubkPy3s8BsqxzGPF+/tkWVL3xhq1ArVYN/v1vC7jZzMN5DxMnwhNPWGw+/XQYMgRuvtkzbvlH1H+jBws3L6ReuXqMv3E8+3bW5cnPx7Bmr83xbtxSGSrl72FKPpsxw4aYk5JgyxaoUsVGT9q2tVuKREQyOe/zb1q2YcOGfsaMGfn2+34nHLbaqP36WX5m7y1D1T33WMaquOwHBL7+2gLvlCmWVbJHD6hx6VSGzX+T6Wuns3jLYmqcVINnWzzLtbWvZeyc9Vnmau7Z/mwttIo1e/fa3ENKiv1lNmKEzeuGw0d9T4lI7HPOzfTeZ9lRLRhXh23brMRQrVoWaCdPhocfhpUrrZJMy5bZXignT7bEGZdearcUpaTYCudZJ/Wi2YgmDJ0zlMVbFtOh+t9ZcPcCOtbpSJyLo/eEJYcFX4A96SHV640l+/ZZfd1KlaBzZ6tw9dprNuQMCr4iclSxfYWYM8du9ShfHh54wIaWR4yANWssH2TVqtm+dPp0i9UXX2wrm195BVasgEbtZtN4yKUMmP+o9aABfBz//XEH4+ZuPPD67Oryql5vlAuHYelSe5yYCBMm2Jvkyy9h0SJbvFeqVLBtFJGoEHvFGPbts3uCUlIsyUGRInDTTZa4vn79Y7585kwbav70U0golkbp5iuo2mIj8fUct36SwshFI0mgBMUzrmJ3/JcHkmjEZ9Sh94QlB4aXyycXITWLYKt6vVFq+3ZLa9avnyXLSE2FEiXgu+8sEIuI/EGxE4BXr4YBAyzp8ubNdstHnz62qKp06WO+fN48ePJJyzJZvGSYsi2W4xqMZ2/SZOam/8x3X8+icEIRejTtwfAJZ+EoRvFQC/bGzadw+GySwrUO690+dEUN1euNBatWWcazESOsKtEFF9gCgKQk267gKyLHKboDcDhsCetTUmD8eHuudWvLKHTppTmag1u0yO4OGTnSRg6ffhrGh//HyrSv2Jz4HPsLJRQNNaVO0n0806I9X333Nanb9pAUrkVSuNaBn3Vo71b1eqNYerrN5550kq0fGDYMbrjBFus1aBB060QkRkRvAH7/feuJLF9uWYUefdQyVVXK2X0+S5daruZ334VixeDxxy3ZVThpKy8935cdiR9xaKGERF+Vzdut15PT3q3q9UaZDRvs1rSBA63Q/fDhlmVlwwbN64pInoveALxrF5QrZ1H02mtzPBS4ciU8+yy89ZaH+DAlzl9FjT+vo3brU+g77z+8NO0ldhTaSeGMeuyNnw8+hCOBwuGzD/Rw1buNMdOnw6uvwqhR1vu94gorKbmfgq+InADRex+w94eV/juW1autHu/QoRAX7yl27ioSGo0lrdgPhN12fov/lpDbQbua7WhW7h4GfJnOttCCA3O8yfFn6R7eWPLbb1C4sE1T/OMftnagSxdbrHfmmUG3TkRixNHuA47eHnAOg29qKjz/vCW8cs6ur9OKT2ENU9iY+DiQAQ4SQ2dSu3BPRl93NwBVSqbSe0Ii67bVUg83lqxYYQUQBg+25Bl//jM89piNpBQrFnTrRKQAid4AfAwbNkCvXrYwOhyG22+Hf/4TylcIUa7HGH5JGAQuw3b2jqLhC9m5o8qB12v+NoaEw3a/7uuvw2efWRnJ9u0PFrovUybY9olIgRRzAXjzZku9m5JidRY6d7a1WlWqeMb8OIYeA3qwJXERCeEKhP1eIPy7OV6JEaGQBdtQCP7yFwvETzwB3bpZchYRkQDFTADeuhVeeskyAe7ZY6l4e/SAatU8E1dMpOObjzNj3QxqnlyT/2s4gE+mV2J7aNFhc7y6RzdGzJtnvd1Jk2DhQiuC8MUXVplI9+2KSISI+gC8bRu8/LJ97dpldXmffBJq1oSBMwby59d6sWrbKqqUqsLQNkO5ue7NJMQl0Lh8Kr0nFNYcb6xIT7fizK+/Dv/7n2VAu/FGe1OULg21awfdQhGRw0RtAN6xw+4ceeklyxJ47bWWUOOss2DW+lk0evMepqdOB6BQXCGGtx1Os6rNDrxec7wxZtIk++vr9NOtrGSXLprbFZGIFrXFGHr0sOm8Zs1g9mxL/5xQ7kc6jexEg0ENmLNxDg5bKR32YaaumRpwiyXPeA9Tp1oP97HH7LlLL7WFVsuW2W1FCr4iEuGiNgA//DD88AN8/DEkV11Fl4+7UKdfHT5b/hlPNH2CRxoOwZEIPg7vE3BpdYJusuTWnj0wZIilg2zc2Cpm7C9yHxdn5avi44Nto4hIDkXtEPTP4Wl8sv0Ten2wlLFLxhIfF88DFz7AI40fYcrSNB4dPZ9TQv86sMhq+KREapZJ1bBzNPv73+HNN22eYcAAq3JVvHjQrRIROS5RmQlr2pppNBvWjPRwOgBta7Tl9atep0JJC66Ne32dZSnACslFmNL9klz/fskH4TB89RX07Wu5Q885xypnbN4MTZv+oSxoIiJBiblMWJNWTSIUtkII8S6e8yucfyD4QvZF77N7XiLIjh1WBCElBZYsgbJl4aefLABrJbOIxJCoDMDNqzYnKSGJtFAaifGJNK/a/LDt5ZOLZNkDVqKNCJeRYfePrV9vdXdHjICOHQ/W3hURiSFRGYAbVWrEV7d+xaRVk2hetTmNKjU6bHtOywVKwEIhGDfO0kP27w8JCZbGrEYNOO+8oFsnInJCRWUA/mh2Kr0n7GHdtrqMS97DQ1ccvrhK5QIj3Nattpiqf38rU1WxIqxbBxUqWAozEZECIOoC8EezUw/r3aZu28Ojo+cD/C4IK+BGoKlT7Z7dvXuhRQvo0wdat7ber4hIARJ19wH3nrDksKFlgD3pIXpPWBJQi+So0tLgvffgP/+x7xs0gDvvhPnz4euvrSqRgq+IFEBRd+XTCucosX49DBpk9+tu2GB1d6+7zhZUvfxy0K0TEQlc1PWAs1vJrBXOEeTFF6FyZUvOXa+eZaz6/POgWyUiElGiLgA/dEUNihQ6PN2gVjgHbO9eGDYMUlPt+7p14W9/s7zMn34KV15pqSJFROSAqLsqtq1XgZ7tz6ZCchEclt2qZ/uzteAqCD//DI8+aquYu3Q5OM/bsqUNM59xRrDtExGJYFE3Bwxa4Ry4cBiuvx5GjbLv27SxHm+LFsG2S0QkikRdD1gCsns3jB1rj+PirNzfQw/BypUwejRcconyM4uI/AFR2QOWfLRiheVlHjIEtm+3ed0zzrDVzSIictzUA5asrVwJrVpB9epWkahlS5g8GapVC7plIiIxIdc9YOdcPDADSPXet8p9kyQwO3ZYSsiaNSE52cr/9egBf/0rlC8fdOtERGJKXgxB3wcsBkrmwc+SIPz4I7z+upUBrFULvv/e5niXL9ftQyIiJ0iurq7OuYrA1cCbedMcyVfffguXX25B9403oF07m+/dT8FXROSEyW0P+BXgYaBEdjs457oB3QAqV66cy18nubZtGxQqBMWKWcH7hQvh2Wfhjjvg1FODbp2ISIFx3F0c51wrYJP3fubR9vPeD/LeN/TeNyxbtuzx/jrJrQULrAhChQq2ohmgc2dYtQoef1zBV0Qkn+WmB9wYaO2cuwooDJR0zr3tvVdB10jhPXz8Mbz2GnzzDRQuDDfeCM2b2/bExECbJyJSkB13APbePwo8CuCcaw78n4JvhNizB4oUscQYffpYL7dnTxtmPumkoFsnIiIoEUdsmTvX7tn98ENb2VyuHLz/PpxyimruiohEmDy5KnvvJwGT8uJnyR+UkXFwmPnbb63ne8stEArZdt2/KyISkdQtilbe2xDzmjXQsSNUqQK9e0PXrnYPr4iIRDQF4Ggze7YNM+/aBR98AKefDlOmwPnnQ3z8sV8vIiIRQZkWokF6OowcCRdfDPXrW93dsmWtLCBAo0YKviIiUUYBOBq8+ip06mR5ml96CdautYxVylQlIhK1NAQdiWbNsmHm1q0tPWTnzlCjBlx1lXq6IiIxQgE4UqSnw5gxtpp5yhRLFXnuubatbFm45poD85DCAAAI6UlEQVRg2yciInlKAThStGwJX38Nf/oTvPwy3HablQQUEZGYpEnEoMyaZbmZd++27++/Hz75BJYutccKviIiMU094PyU1TDzzTdDkyYaYhYRKWAUgPPLpk12C1Fq6sFh5i5doFSpoFsmIiIBUAA+kWbPtvzMt91m+Zg7dIDLLoMrr9RqZhGRAk4BOK+lp8NHH9kw8+TJcPLJcMMNkJQEr7wSdOtERCRCaBFWXpo40YaX9yfNePllWLbMgq+IiMgh1APOrTlzrLB97dpWEKFmTejXT0kzRETkqNQDPh4ZGTBqFDRrBvXqwbPP2vM1asAXX9iKZgVfERE5CgXgP2rgQBtm7tABfv4Z/v1v6/GKiIj8ARqCzomFC21oOT7ebiOqXt1yNbdqpZ6uiIgcFwXg7IRCMHasrWaeNMlWNrdpA089pSpEIiKSa4okR9q3z4aVq1WD9u1h5Up48UWrxQsKviIikifUA97v11+hdGlISLA53apVoU8fKwmYoP8mERHJWwU7soTD8OmnNsw8dy6sXg2FC1uhBBVDEBGRE6hgBuDt22HoUHj9dVixAipUsApEoZBtV/AVEZETrGAF4HDY5nDnzYMHHoCLLoLnnrO53kKFgm6diIgUILEfgMNhS47x6qu2sKpvXyv/N3cu1K0bdOtERKSAit0lvbt22WKqOnWgZUurTFSlim1zTsFXREQCFbs94O7dISUFzjsP3n4bOna0nM0iIiIRIDZ6wN7DN99Au3Ywfbo998ADMHUqfP893HSTgq+IiESU6A7Ae/bAm2/COefAJZdY/d01a2xbtWrQqFGw7RMREclG9A5Bh8MWeJcts/ncwYOt8H2RIkG3TERE5JiiNwDHxcGTT0LFitC0qS2sEhERiRLRG4DB5nZFRESiUHTPAYuIiESp4w7AzrlKzrlvnHOLnXMLnXP35WXDREREYlluhqAzgH9472c550oAM51zX3jvF+VR20RERGLWcfeAvffrvfezMh/vBBYDFfKqYSIiIrEsT+aAnXNVgXrAd3nx80RERGJdrgOwc644MAq433u/I4vt3ZxzM5xzMzZv3pzbXyciIhITchWAnXOFsOD7jvd+dFb7eO8Hee8beu8bli1bNje/TkREJGbkZhW0AwYDi733ffKuSSIiIrEvNz3gxsAtwCXOuTmZX1flUbtERERi2nHfhuS9nwwo/6OIiMhxUCYsERGRACgAi4iIBEABWEREJADOe59/v8y5zcDqPPyRJwNb8vDnBUnHEnli5ThAxxKpYuVYYuU4IO+PpYr3Pst7cPM1AOc159wM733DoNuRF3QskSdWjgN0LJEqVo4lVo4D8vdYNAQtIiISAAVgERGRAER7AB4UdAPykI4l8sTKcYCOJVLFyrHEynFAPh5LVM8Bi4iIRKto7wGLiIhEpagIwM65ls65Jc655c657llsd8651zK3z3PO1Q+incfinKvknPvGObfYObfQOXdfFvs0d85tPyS/9hNBtDUnnHOrnHPzM9s5I4vtEX9enHM1Dvm/nuOc2+Gcu/+IfSL2nDjnhjjnNjnnFhzyXBnn3BfOuWWZ/5bO5rVH/Vzlt2yOpbdz7sfM988Y51xyNq896nsxv2VzLE8551KPlTs/ks5LNsfxn0OOYZVzbk42r420c5Ll9TfQz4v3PqK/gHhgBfAnIBGYC9Q+Yp+rgM+w3NQXAt8F3e5sjuU0oH7m4xLA0iyOpTkwLui25vB4VgEnH2V7VJyXQ9obD2zA7tuLinMCNAXqAwsOee5FoHvm4+7AC9kc61E/VxFyLJcDCZmPX8jqWDK3HfW9GCHH8hTwf8d4XUSdl6yO44jtLwFPRMk5yfL6G+TnJRp6wOcDy733K733acD7QJsj9mkDvOXNdCDZOXdafjf0WLz36733szIf7wQWAxWCbdUJFRXn5RCXAiu893mZLOaE8t5/C/xyxNNtgOGZj4cDbbN4aU4+V/kqq2Px3k/03mdkfjsdqJjvDTsO2ZyXnIio83K048gsSdsJeC9fG3WcjnL9DezzEg0BuAKw5pDv1/L7oJWTfSKKc64qUA/4LovNjZxzc51znznn6uRrw/4YD0x0zs10znXLYnu0nZfryf5iEi3nBOBU7/16sIsOcEoW+0TbuQHoio2oZOVY78VI8bfM4fQh2Qx1RtN5uRjY6L1fls32iD0nR1x/A/u8REMAzqrk4ZFLt3OyT8RwzhUHRgH3e+93HLF5FjYEeg7QF/gov9v3BzT23tcHrgTucc41PWJ71JwX51wi0BoYmcXmaDonORU15wbAOfcYkAG8k80ux3ovRoL+QDXgXGA9Nnx7pGg6Lzdw9N5vRJ6TY1x/s31ZFs/l+rxEQwBeC1Q65PuKwLrj2CciOOcKYSf/He/96CO3e+93eO93ZT7+FCjknDs5n5uZI977dZn/bgLGYMM0h4qa84JdJGZ57zceuSGazkmmjfuH+jP/3ZTFPlFzbpxznYFWwE0+c0LuSDl4LwbOe7/Rex/y3oeBN8i6jVFxXpxzCUB74D/Z7ROJ5ySb629gn5doCMA/ANWdc6dn9lKuB8Yesc9Y4NbMVbcXAtv3DylEksw5k8HAYu99n2z2KZe5H86587FztDX/WpkzzrlizrkS+x9ji2UWHLFbVJyXTNn+NR8t5+QQY4HOmY87Ax9nsU9OPleBc861BB4BWnvvf8tmn5y8FwN3xPqHdmTdxqg4L8BlwI/e+7VZbYzEc3KU629wn5egV6bl5AtbTbsUW4X2WOZzdwJ3Zj52QErm9vlAw6DbnM1xNMGGLeYBczK/rjriWP4GLMRW2U0HLgq63dkcy58y2zg3s73RfF6KYgG11CHPRcU5wf5oWA+kY3+l3w6cBHwFLMv8t0zmvuWBTw957e8+VxF4LMuxubf9n5cBRx5Ldu/FCDyWEZmfg3nYxfu0SD8vWR1H5vPD9n8+Dtk30s9JdtffwD4vyoQlIiISgGgYghYREYk5CsAiIiIBUAAWEREJgAKwiIhIABSARUREAqAALCIiEgAFYBERkQAoAIuIiATg/wETa7hGdeCWmgAAAABJRU5ErkJggg==\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.4rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
