{
 "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, 11 Sep 2020                                         \n",
      "Time:                        06:24:40                                         \n",
      "No. Iterations:                    19                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "var_0        -41.0265      9.792     -4.190      0.000     -60.218     -21.835\n",
      "var_1          0.8294      0.111      7.472      0.000       0.612       1.047\n",
      "var_2          0.9261      0.303      3.057      0.002       0.332       1.520\n",
      "var_3         -0.1278      0.129     -0.994      0.320      -0.380       0.124\n",
      "==============================================================================\n",
      "\n",
      "If the model instance has been used for another fit with different fit parameters, then the fit options might not be the correct ones anymore .\n"
     ]
    }
   ],
   "source": [
    "huber_t = sm.RLM(data.endog, data.exog, M=sm.robust.norms.HuberT())\n",
    "hub_results = huber_t.fit()\n",
    "print(hub_results.params)\n",
    "print(hub_results.bse)\n",
    "print(hub_results.summary(yname='y',\n",
    "            xname=['var_%d' % i for i in range(len(hub_results.params))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huber's T norm with 'H2' covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-41.02649835   0.82938433   0.92606597  -0.12784672]\n",
      "[9.08950419 0.11945975 0.32235497 0.11796313]\n"
     ]
    }
   ],
   "source": [
    "hub_results2 = huber_t.fit(cov=\"H2\")\n",
    "print(hub_results2.params)\n",
    "print(hub_results2.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Andrew's Wave norm with Huber's Proposal 2 scaling and 'H3' covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [-40.8817957    0.79276138   1.04857556  -0.13360865]\n"
     ]
    }
   ],
   "source": [
    "andrew_mod = sm.RLM(data.endog, data.exog, M=sm.robust.norms.AndrewWave())\n",
    "andrew_results = andrew_mod.fit(scale_est=sm.robust.scale.HuberScale(), cov=\"H3\")\n",
    "print('Parameters: ', andrew_results.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See ``help(sm.RLM.fit)`` for more options and ``module sm.robust.scale`` for scale options\n",
    "\n",
    "## Comparing OLS and RLM\n",
    "\n",
    "Artificial data with outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "x1 = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x1, (x1-5)**2))\n",
    "X = sm.add_constant(X)\n",
    "sig = 0.3   # smaller error variance makes OLS<->RLM contrast bigger\n",
    "beta = [5, 0.5, -0.0]\n",
    "y_true2 = np.dot(X, beta)\n",
    "y2 = y_true2 + sig*1. * np.random.normal(size=nsample)\n",
    "y2[[39,41,43,45,48]] -= 5   # add some outliers (10% of nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1: quadratic function with linear truth\n",
    "\n",
    "Note that the quadratic term in OLS regression will capture outlier effects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 5.05647156  0.5175479  -0.01320007]\n",
      "[0.48610764 0.07504846 0.00664063]\n",
      "[ 4.72646991  4.98939267  5.24791725  5.50204364  5.75177184  5.99710184\n",
      "  6.23803366  6.47456729  6.70670273  6.93443999  7.15777905  7.37671992\n",
      "  7.5912626   7.8014071   8.0071534   8.20850151  8.40545144  8.59800318\n",
      "  8.78615672  8.96991208  9.14926925  9.32422822  9.49478901  9.66095161\n",
      "  9.82271602  9.98008224 10.13305027 10.28162011 10.42579176 10.56556523\n",
      " 10.7009405  10.83191758 10.95849648 11.08067718 11.1984597  11.31184402\n",
      " 11.42083016 11.52541811 11.62560786 11.72139943 11.81279281 11.899788\n",
      " 11.982385   12.06058381 12.13438443 12.20378686 12.2687911  12.32939716\n",
      " 12.38560502 12.43741469]\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.98020450e+00  5.04277495e-01 -2.32470185e-03]\n",
      "[0.14193966 0.02191357 0.00193901]\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 0x7f02e9d9d490>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAHSCAYAAADlm6P3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAACRmElEQVR4nOzdd7yO9R/H8dd1pmMee++VvY6dEhUhREOJlLJKiij9Ku1C2VtKQ0iJUKlEQ/besjn2OObZ9/X74xMnWQfnnPuM9/PxuB/nnOte37vOcX/u7/UZjuu6iIiIiIikBT7eXoCIiIiISFJR8CsiIiIiaYaCXxERERFJMxT8ioiIiEiaoeBXRERERNIMBb8iIiIikmb4JeWT5ciRwy1SpEhSPqWIiIiIpEErVqw46rpuzv8eT9Lgt0iRIixfvjwpn1JERERE0iDHcXZf7rjSHkREREQkzVDwKyIiIiJphoJfEREREUkzrpnz6zjOx0Az4LDruuX/dbw78AwQA8xxXbfPjSwgOjqaffv2ERERcSN3TxHSpUtHgQIF8Pf39/ZSRERERNK0+BS8TQRGAJ+dP+A4zh1AC6Ci67qRjuPkutEF7Nu3j0yZMlGkSBEcx7nRh0m2XNfl2LFj7Nu3j6JFi3p7OSIiIiJp2jXTHlzX/R04/p/DXYH3XdeN/Oc2h290AREREWTPnj1VBr4AjuOQPXv2VL2zLSIiIpJS3GjObymgnuM4SxzH+c1xnOpXuqHjOJ0cx1nuOM7yI0eOXOk2N7iMlCG1vz4RERGRlOJGg18/ICtQC+gNfOVcIcJzXXec67ohruuG5Mx5SZ/hZOf111/ngw8+uOL1M2bMYOPGjUm4IhERERFJKDca/O4DprtmKeABciTcsq5sxqpQ6r7/K0VfmkPd939lxqrQpHjauOdX8CsiIiKSYt1o8DsDaADgOE4pIAA4mkBruvKTrgql7/R1hIaF4wKhYeH0nb7upgPgd955h9KlS3PnnXeyZcsWAMaPH0/16tWpVKkSrVu35ty5c/z1119899139O7dm8qVK7N9+/bL3k5EREREkqdrBr+O40wGFgGlHcfZ5zhOR+BjoJjjOOuBKcBjruu6ibtUGDh3C+HRsRcdC4+OZeDcLTf8mCtWrGDKlCmsWrWK6dOns2zZMgBatWrFsmXLWLNmDWXKlGHChAnUqVOH5s2bM3DgQFavXk3x4sUvezsRERERSZ6u2erMdd2Hr3DVowm8lmvaHxZ+Xcfj448//uC+++4jffr0ADRv3hyA9evX88orrxAWFsaZM2do1KjRZe8f39uJiIiIiPelqAlv+YKDrut4fF2uVq9Dhw6MGDGCdevW0a9fvyu2Kovv7URERETE+1JU8Nu7UWmC/H0vOhbk70vvRqVv+DFvu+02vv32W8LDwzl9+jSzZs0C4PTp0+TNm5fo6GgmTZp04faZMmXi9OnTF36+0u1EREREJPmJz4S3ZKNllfyA5f7uDwsnX3AQvRuVvnD8RlStWpWHHnqIypUrU7hwYerVqwfAW2+9Rc2aNSlcuDAVKlS4EPC2adOGp556imHDhvH1119f8XYiIiIikvw4SVCndkFISIi7fPnyi45t2rSJMmXKJNkavCWtvE4RERGR5MBxnBWu64b893iKSnsQEREREbkZKSrtQURERESSoeho2LQJVq6Mu2TNCv/UUiUnCn5FREREJP4iImDdOvj7b3jkETt2330wZ459nyEDVKkClSp5b41XoeBXRERERK7ul19g0iTb0d2wAWJjwXGgeXPImBGefhratoWqVaFECfD1vfZjeomCXxERERGB48dh1aqLUxd+/BGKFrWUhu+/h2rV4N57LcitWtV2eQHuuce7a78OCn5FRERE0poDByzQLV8eChWy3Nx/ptwCdqxqVUtxAOjWDZ55xnZ7U7g0H/weO3aMhg0bAnDw4EF8fX3JmTMnAEuXLiUgIMCbyxMRERG5ecePw+DBcTu6Bw/a8eHDLaitVg3697eAt0oVyJ794vsn4zSG65Xmg9/s2bOzevVqAF5//XUyZszICy+8cOH6mJgY/PzS/H8mERERSe48Hti2zYLb8+kLd90FffqAv78Ft6VLw913xwW5VarYffPls9ulAYrqLqNDhw5ky5aNVatWUbVqVTJlynRRUFy+fHlmz55NkSJF+OKLLxg2bBhRUVHUrFmTUaNG4ZuKPh2JiIhIMhQTA1u2wKlTULs2uC4ULgz79tn1AQFQoQIEBdnPmTLZbdOl896ak4lkFfw+9xz8swmbYCpXhiFDrv9+W7du5ZdffsHX15fXX3/9srfZtGkTU6dOZeHChfj7+9OtWzcmTZpE+/btb2bJIiIiIpf69luYNw9WrIA1ayA83HZuV660XNxevSBzZtvVLVvWAuB/U+ALJLPgNzl54IEHrrmDO2/ePFasWEH16tUBCA8PJ1euXEmxPBEREUmNzvfQXbHCgtrdu2HuXLvuq6+sl26VKtCliwW51arF3fe557yy5JQmWQW/N7JDm1gynG/dAfj5+eHxeC78HPFP5aPrujz22GO89957Sb4+ERERSeHCw20Ht2pV26UdMAD+9z9LaQCbkFa1Kpw7B+nTw9ix1lPXx8e7607h9F8vHooUKcLKlSsBWLlyJTt37gSgYcOGfP311xw+fBiA48ePs3v3bq+tU0RERJKxvXutu0KHDpaPmymT5etu2GDXV6sGvXvD11/Djh1w7JgNl0if3q7PnFmBbwJIVju/yVXr1q357LPPqFy5MtWrV6dUqVIAlC1blrfffpu7774bj8eDv78/I0eOpHDhwl5esYiIiHjN2bNWxLRihV2eegpuvRU2b4Znn4VcuSzQbdnSvhYrZvdr2NAukqgc13WT7MlCQkLc5cuXX3Rs06ZNlClTJsnW4C1p5XWKiIikKefO2SVHDggNhUaNbBra+XTJPHlg0CB4+GFLczh+3NqKpYJhEcmd4zgrXNcN+e9x7fyKiIiIxIfrwpIlsHy5XVasgI0boWtXGDECcueGEiWgdWsICbFd3Xz54u4fFAT583tv/QIo+BURERG5VFSUdV1YvtyC3i5dbLf2gQesl27u3Bbc3nefDZIA8PODGTO8umy5NgW/IiIikrZ5PHGFZAMHWkuxtWstAIa41mJgxWj589tFqQtXFRlpdXvJLetTwa+IiIikHedHAC9bFnf5+284cAB8feHIEevC8NxzlroQEgJFisTdv2ZNb608xTh50rqyDRkCgYH2n9cvGUWcyWgpIiIiIgksNBSWLrXUhIwZ4d134dVX7bqgIOuj27atFa1lymS9duWG7N9vAe+YMXD6tP0n79PHPlMkJwp+RUREJPXYswc+/9wC3mXLbEcXbCxwgwbQvDnkzQvVq9sI4OS0JZlEZqwKZeDcLewPCydfcBC9G5WmZZUbL8TbtAk++MD+s8fGwoMPWrviqlUTcNEJSJ2SgX379tGiRQtKlixJ8eLF6dGjB1FRUSxYsIBmzZpdcvvZs2dTpUoVKlWqRNmyZRk7dqwXVi0iIpKGRUVZcDtiBLRvDz/8YMcPHYJXXoEtW6xn7tChsGgR1Klj11esCB072tc0Gvj2nb6O0LBwXCA0LJy+09cxY1XodT/WX39Zq+KyZWHyZOjUyVIcJk9OvoEvaOcX13Vp1aoVXbt2ZebMmcTGxtKpUyf+97//0bRp00tuHx0dTadOnVi6dCkFChQgMjKSXbt2Jf3CRURE0grXtcERGTNaQmmjRrBqVVxBWp48cPvt9n2VKnDiBAQHe225ydnAuVsIj4696Fh4dCwD526J1+6vxwNz5lh2yJ9/QrZs8Npr8MwzkDNnYq06YaX54PfXX38lXbp0PP744wD4+voyePBgihYtyh133HHJ7U+fPk1MTAzZs2cHIDAwkNKlSyfpmkVERFK148ctbWHxYrssXQpNmsAXX9iI3zx5oEcPKz6rUQMKFIjrvODnp8D3KvaHhV/X8fOiouDLL60ZxsaNUKiQbap37AgZMiTGShNP8gp+n3vOxgEmpMqVLfv6CjZs2EC1atUuOpY5c2YKFSrEtm3bLrl9tmzZaN68OYULF6Zhw4Y0a9aMhx9+GB/N2hYREbl+MTHWT3ffPrj3Xjt2++2wfr21HytfHu6/33Z7wYJc9dK9YfmCgwi9TKCbLzjosrc/dQrGj4fBg612sGJF+wzy4IPg75/Yq00cySv49QLXdXEu06fvSscBPvroI9atW8cvv/zCBx98wM8//8zEiRMTeaUiIiKpxJ9/wqxZtqu7fLl1WggOhmPHLODt3986MYSEWAcGuarrKWDr3ag0faevuyj1Icjfl96NLj6LffCg7eyOHm2ZJvXrw0cf2WeQlN7eOHkFv1fZoU0s5cqV45tvvrno2KlTp9i7dy/Fixe/4v0qVKhAhQoVaNeuHUWLFlXwKyIi8l/R0XZGd9EiC3Q/+gjSp4fvv7etxKpV4amnoFYtu5yPqpo08eqyU5LzBWzng9nzBWzAZQPg88euFCxv3WqdGz791P73tW5tnRtq1EiiF5QEklfw6wUNGzbkpZde4rPPPqN9+/bExsbSq1cvOnToQPr06S+5/ZkzZ1i+fDn169cHYPXq1RQuXDiJVy0iIpIMua4FsH/8Af/7n+3qhv9zir1AAdi928Z99eljVVLp0nl3vanAjRSwtayS/5Lrli61Dfdvv4WAAHj8cejVC0qWTLSle02aT1R1HIdvv/2WadOmUbJkSUqVKkW6dOl49913AZg3bx4FChS4cFm1ahUDBgygdOnSVK5cmX79+mnXV0RE0p7YWFizxiYatG9vUdL5XNyAAJtt27kzTJ1qvXf37o2bcxscrMA3gdxoARvYZ5UffrCUhpo14ddfoW9f+4wyZkzqDHxBO78AFCxYkFmzZl1yvH79+oSHX/rLU69evaRYloiISPJx9qyN7cqTxyqfypa1aiiA3Lmhdm3ImtV+rlkTlizx3lrTkOstYANLZ5g61dqVrVtnm/KDBsGTTyZsivWivYtYsGsB9YvUp3bB2gn3wDdJwa+IiIhcav9+WLjQitMWLrTc3Xbt4JNPIF8+63FVtaoNjyhaNOVXQaVQ8S1gAzhzxtKuBw2yjfhy5Sy3t00b26xPSIv2LqLBZw2Iioki0C+Qee3nJZsAWMGviIhIWufxWPPWvXvhnnvs2B13WPVTUJDt5L700sXtxgYN8t565YJrFbABHD4Mw4fDyJE2/+O222DUKKsrTOhOrR7Xw/yd83lu7nNExEQAEBUbxYJdCxT8ioiIiBetWwc//mjFaX/+aVFR9uxw5IgFt8OHW25ulSopt6FrGnG5AjaA7dvhww9tsz4y0kYR9+ljjTUS2pGzR5i4eiLjVo5j2/FtZArIhJ+PH67rEuAbQP0i9RP+SW+Qgl8REZHU7uxZazf2559W0RQYaOe7P/zQqppatYJ69eDWW+Puc/fd3luv3JQVKyyf9+uvbeBd+/bwwguQ0ANpXdfl992/M3bFWL7Z9A1RsVHcWuhW+t3ej/vL3s+qA6uU8ysiIiJJZMsWG831xx+wcqVNUnMcaNHCdnN79rSIKE8eb69UEoDrwi+/WLuyefNsCnTv3jYFOm/ehH2u4+HH+XT1p4xbOY7NRzeTJTALXap1oVO1TpTLVQ44P3gjnP1hFZkdHE7vRqFXbL2W1BT8ioiIpHRHj8Lvv8OCBTZ39tZb4dAhS12oUcOioHr1rDgtSxa7T758Xl2yJIyYGJg2zXZ6V6+2/60DBliXucyZE+55XNdl3IpxjFkxhg2HNxDtiaZWgVp80uITHiz3IOn942YjXO/gjaSm4Bfw9fWlQoUKxMTEULRoUT7//HOCg4PZtWsXzZo1Y/369RfdvkOHDnz11VccOnSITP/0BOnRowfDhg3jyJEj5MiRwxsvQ0RE0pIzZ6wIbcEC2LDBjgUFWQn/rbdaoHvypPrpplLnzsHHH1vmyq5dcMst9vMjj1hWS0I5GXGSz9d+zqBFg9gZthMAP8ePz1p+RrtK7S57nxsZvJGU0vyQC4CgoCBWr17N+vXryZYtGyNHjrzmfUqUKMHMmTMB8Hg8zJ8/n/z5vf8/VEREUqEDB2DyZOjSBd56y46lTw+zZkH+/PDuu9aOLCzMtvzAkj0V+KY6R4/CG29AoULQvbulNMyYYZ9/Hn88YQJf13VZFrqMjjM7km9QPrr/0J2o2CgcrJ2di8u+U/uueP+bGbyRFFLkzm9iNk2uXbs2a9euvebtHn74YaZOncqjjz7KggULqFu3Lj/88EOCrkVERNK4N96AKVNg82b7OVMmePRR+97HB3buTPheVZIs7dplu7wTJtjE6HvvhRdfhLp1E+45TkeeZvL6yYxZPoZVB1eRwT8DbSu0pXO1zkTFRtHws4ZExUZds3vDjQzeSErJKvh97sfnWH1w9VVvczLyJGsPrcXjevBxfKiYuyJZArNc8faV81RmSOMh8Xr+2NhY5s2bR8eOHa9525IlSzJz5kxOnDjB5MmTefTRRxX8iojIjTlxAn77DebPty28n3+24rQjR2yAxBNPWN/dypVtR/c8Bb7JnhV+XbkH77WsXm05vF99Zf+727a1FO6yZRNujasPrmbs8rF8se4LzkSdoUKuCoxqMoq2FduSOTAucXhe+3nx2ny8nsEb3pCsgt/4OBlxEo/rAayR8smIk1cNfuMjPDycypUrs2vXLqpVq8Zdd90Vr/u1atWKKVOmsGTJEsaOHXtTaxARkTRoxgx4+23rxuC6lqZw6602NjhLFhgxwtsrlJtwo4Vfrgu//mpB708/QcaM8NxzdilQIGHWNn/nfEYuG8mmI5vYeHQj6fzS8VC5h+hcrTO1CtTCuczEvtoFa8frjHt8Bm94U7IKfuOzQ7to76KLtt0ntZp006kP53N+T548SbNmzRg5ciTPPvvsNe/Xpk0bqlatymOPPYaPPn2LiMiVREXB4sXWg+qXX6wLQ9Wq4OsLGTJAv362s1uzZsJWK4lXXW/hV2wsfPONBb0rVkDu3JbO3bWrzRtJCJuObOL1Ba/z1cavAHBw6FGzB6/d/hrZgrIlzJNw5cEbyUGyCn7jo3bB2vHedr9eWbJkYdiwYbRo0YKuXbte8/aFChXinXfe4c4770ywNYiISCqyZ48Vqf32m5Xn+/hA9erWqQEscfPee727Rkk08S38Cg+HiRPhgw9gxw6bOzJ2rA2nSIiaxciYSL7d/C1jlo/ht92/4ePEbdj5OD7kzpA7QQPf5O6awa/jOB8DzYDDruuW/891LwADgZyu6x5NnCVeKr7b7jeiSpUqVKpUiSlTplCvXj22bNlCgX+dYxg8ePBFt+98vqpWRETStp0743Z2a9SwIRLZs8O+fZaz27Ah1K+fcFt4ctNuNh/3Wq5V+HX8OIwaBcOGWXp3jRowcKDNIfH1vfnn33FiB+NWjOPjVR9z5NwRigYX5f2G71M+d3ke+OqBeBWvpUbx2fmdCIwAPvv3QcdxCgJ3AXsSfllJ68z5T+D/mDVr1oXvo6OjL7n9Aw88cNnH2bVrV4KuS0REUoBevSx3d8cO+zlvXqhQwb7PkAHi0UFIkl5SDGK4UuFXh4plef55G8B39izcc491brjtNqtzvBkxnhhmb53NmOVjmLt9Lr6OL/eWvpcu1bpwV/G7Luz6JtZZ9JTgmsGv67q/O45T5DJXDQb6ADMTelEiIiLJTkwMLF1qFUi7d8Mnn9jxAwegfHmrRrrzTps2cLMRjCS6pBjE8N/CrywROci6vgLd3rdpaA8/bBOmK1a8+ef6bvN3DFs6jDWH1nD03FHyZ8rP67e/TseqHSmQ+dIqucQ8i57c3VDOr+M4zYFQ13XXXK4aUEREJNX46ScYM8ZSGk6dsrzdGjUgMtKK07780tsrlBuQVIMYWlTOT7bT+enfH77/3k4GPPMMPP+8Daq4GR7Xw8/bf+bdP9/l992/A5bD+37D9+lVpxd+PimutCtJXPd/Fcdx0gP/A+6O5+07AZ3ACsRERESSrdOnrcfU3LnQty8ULGi5vCtWwEMPwd13Q4MGkC3tFAelVok9iMHjgZkzoX9/WLIEcua04Xzdut38r8+Rs0f4ZPUnjF0xlh0ndpDBPwMODi4uDg4e16PA9ypu5L9McaAocH7XtwCw0nGcGq7rHvzvjV3XHQeMAwgJCXEv94Cu6162n1xq4bqXfdkiIpIcHD1qY7N+/BH+/NPSGzJmtKqjggXhySehUyelMqQyiTWIISICvvjCCte2boVixayorUMHCLqJuNp1XRbuXcjo5aP5euPXRMVGcVvh23inwTvkzZiXeybdk2YL2K7XdQe/ruuuA3Kd/9lxnF1AyI12e0iXLh3Hjh0je/bsqTIAdl2XY8eOkU7z1UVEkocTJ6wjQ44c1ls3OhpeeskSL3v1gsaNoU4dCAiw2ydE2b0kOwk9iCEszLJjhg6FgwetjfPUqdC69c39Cp2MOMkXa79gzIoxrD+8niyBWehcrTNdQrpQNmfcmLe0XMB2vZxr7Uo6jjMZqA/kAA4B/VzXnfCv63cRz+A3JCTEXb58+UXHoqOj2bdvHxEREde9+JQiXbp0FChQAH9/f28vRUQkbVq50hIuf/wRFi2yc9IPPmjRCcDhw5Ar19UfQ+QyQkNhyBDry3v6tGXG9Olj2TE3uqe3aO8iJq2bxN6Te5m3cx5no88Ski+ELtW60KZ8GzIEZEjQ15BaOY6zwnXdkEuOJ+Up+csFvyIiIgnu5ElrMVavnv1cowYsWwbVqllfqcaNbZqan/Ii5cZs2mSpDV98YZPZHnoIeveGKlVu/DHPRZ/j3T/e5b0/38PjegBoVqoZ/W7vR0i+S2I4uYYrBb/6qxcRkZTPdS0amTPHdnj//NMC22PHIH16a6iaN692d+WmLVxo44e/+85yeDt3tnkmRYve+GNuObqFMcvHMHHNRMIiwi4c93V8qVOgjgLfBKbgV0REUqaICGs7FhAAgwdbvi5Y7m7v3tC0qbUiA6hUyXvrlBTP44HZs61zw19/WbeG116zlmU5c97YY0bHRvPdlu8YvXw083bOw8/Hj1ZlWnFbodvo/XNvFa8lIgW/IiKScuzbZ1HInDnWd3fKFGje3ALdDBmgSRPr0CCSACIjrY3zwIF2YqFwYRtF/MQT9ut2I/ad2se4FeP4aOVHHDhzgEJZCvFOg3d4osoT5MmYB4CqeauqeC0RKfgVEZHk7/Bhy9Ndtcp+LloUOnaMO9dcurRdRBLAqVMwbpydUNi/304cTJoEDzwAN1K7vnDPQj5a+RHbT2xn4d6FuK7LPSXvYVzIOO4pcQ++Phe3g0jL09eSgoJfERFJXiIibNDErFmQPTu8/badWy5cGNq0gXvv1QhhSRQHDlirstGjLQBu0AA+/tg6ONzIr9uxc8fot6Afo5aNwsUaDLSr0I437niDollvIklYboqCXxERSR6++cZK53/6Cc6ds0ET7drZdY4D337r3fVJqrVlC3zwAXz2mc04ad3a0sarV7/+x3JdlyWhSxi9fDRT108lMjbywnW+ji9lcpZR4OtlCn5FRMQ7Nm+GH36AHj2scO3nn22McIcOlsdbv35cwZpIIli82Do3zJhhdZNPPGF1kyVKXP9jnY06y5frvmTU8lGsPriajAEZ6VilI7UL1KbT7E4qYEtG1OdXRESShsdjvXZnzLBd3C1b7Pi6dVC+PJw9a23JlM4gich1rRvegAHw++8QHAzdusGzz0Lu3Nf/eBuPbGT0stF8tvYzTkWeomLuinQN6UrbCm3JFJgJsKEVKmBLeurzKyIiSS8qynJ4M2eGuXOtG4Ofn+3qdu8OLVpAgQJ22xstnxeJh+homDzZOjesX29NQQYPtrrJTJmu77F+3/07o5ePZvPRzaw+uJoA3wAeKPsAXUO6UqdgHZz/fIBTAVvyouBXREQS1unTls4wY4a1JHv+eXj9dQt4P//c2pJlzerlRUpaceaMzTgZPBj27rWTDJ99ZrWT19u5Ye/Jvbw2/zUmrpkIgIND15CuvFH/DXJmuMGGv5LkFPyKiEjCcF3rBTVrlu345sxpPzdoYNcHBcGjj3p3jZJmHD5sPXlHjYITJ+D222HMGJtufT2ZNR7Xwy87fmHUslHM2jrrwthhAB/Hh4KZCyrwTWEU/IqIyI3Zv99ydzdssAjDcWyE8DPPQMuWUKcO+Ppe82FELmfGqlAGzt3C/rBw8gUH0btRaVpWyX/N+23fbp0bJk60IRX33Qd9+kDNmtf3/MfDjzNx9URGLx/NtuPbyJk+Jy/WfZGQfCE8Ov1RFbClYAp+RUQk/vbtg6++srZkixbZbm+ZMlasliEDDB/u7RVKKjBjVSh9p68jPDoWgNCwcPpOXwdwxQB4xQobP/zNN5ZW/thj1rnhemefLN+/nFHLRjF5/WQiYiKoW7Aub9R/g9ZlWhPoZ91H5rWfpwK2FEzdHkRE5Oq2brUUhqxZLXmyUyeoXNmaobZubcGvSAKq+/6vhIaFX3I8f3AQC19qcOFn17UOef3721yULFmga1fr3JA3b/yea9HeRfy842ciYyL5ecfPLNu/jAz+GXi04qN0DelKpTyVEuplSRJTtwcREYm/jRth2jT4+msrjR81yqKKBx+0HN7ixb29QknF9l8m8P338ZgYOwExYACsWQP58lkXh06drLFIfH21/ivaftuWGE8MAEWCizD8nuG0q9iOLOmy3PTrkORJwa+IiMQJD4caNSzgdRy49VYYMsSGToBtrWVRUCCJK19w0GV3fnOnz8Dw4TBoEOzaZScdPv4Y2ra1IRXxEeuJZc7fcxi1bBRzt8+9cNzH8eGpqk/xTI1nEuhVSHKl4FdEJC3bvNl2eI8ft15QQUFw223QubOlNMT33LFIAurdqPRFOb+x5/wJX12MjeuK8mwY1K1rnRyaNrXhgPFx+OxhJqycwJgVY9hzcg/5MuXjySpP8sW6L4iOjSbAN4A7ityReC9Kkg0FvyIiac327dbt/6uvbLqa41gqg8djkcTIkd5eoaRx54va3pqym79/ycfZdQXxRPvSvLl1bqhbN36P47oui/YtYtSyUUzbOI2o2CgaFG3AoLsH0bx0c/x9/XmiyhMqXktjFPyKiKQF27bZJLV06eDLL+G11yyCGDrUdnjzX7uFlEhSWbUKvhqYnzVf5cfHBzq0gxdeiF9t5aK9i/hp+0+ciznH3G1zWXNoDZkDM9O5Wme6hnSlTM6LH0TT19IeBb8iIqnV+bZkkyfD8uXWA6pVK0tpePzxuLHCIsmA61rHhgED4KefbORwz57Qo0f8P5tNXj+Z9t+2v1DAVjxrccY2G8sjFR4hY0DGRFy9pCQKfkVEUptjxyzI/eMPiyiqVbNS+Nr/7G7lyuXd9Yn8S0wMTJ9uQe+KFZAnD7z/vn1GCw6Ox/09MczaMouRy0Yyb+e8C8d9HB+eqPIEnap1SrzFS4qk4FdEJKU7dQpmzIDTp+HppyFbNkifHl5/Hdq0gVKlvL1CkUuEh9sUtg8+gB077Nd0/HibgJ0u3bXvf+jMIcavHM/YFWPZd2ofBTIXoFPVTny29jMVsMlVKfgVEUmJzp2DOXNgyhT7GhlpO7zdulkB2w8/eHuFIpd1/Li1jR42DI4csbHDH3xg3fSuNQ3bdV0W7l3IqGWj+Hrj10R7ormz2J0Mv2c4zUo1w8/Hjw6VO6iATa5Kwa+ISEoRG2vdGBwHXnrJRgnnyWPnh9u0gVq17DqRZGjPHuvP+9FHNg27SRN48UWoV+/Kv7aL9i5iwa4F1Mxfk20ntjFy2UjWHlpL5sDMdA3pSrfq3Sid4+L5xSpgk2tR8Csikpy5LqxcCV98Ybu806db7u7TT0OLFlC//rW3y0S8aN06y+edPNmC3Icfht69oUKFq99v0d5FNPisAZExkbi4AFTMXZGxzcbStkJbMgRkSILVS2qk4FdEJDk6fdrakH3xBWzZYuOrmjaFwEC7vnRpu4gkQ64Lv/8O/ftbBk6GDPDss/Dcc1Co0NXvG+uJZfbW2bzw0wtExEQA4ODQqVonRjcdjaOzG3KTFPyKiCQXR47YzNbq1S3YHTzYtsd69YL774esWb29QpGrio212ssBA2DpUsiZE95+G7p2tTrMqzly9ggfrfzowgS2HOlz4Ofjh+u6BPgG8FilxxT4SoJQ8Csi4k3h4TBrFnz2GcydC0WL2k5vYCDs3AmZM3t7hSLXFBFhv8IffAB//w3Fi8Po0fDYYzYx+0pc12Vp6FJGLhvJ1A1TiYqN4o4idzC40WCal27OstBlKl6TBKfgV0TEW0aPhr594eRJ6+Lfsye0bRtX/aPAV5K5sDD7NR46FA4dgpAQmDYN7rvv6qno4dHhTFk/hZHLRrLiwAoyBmTkqapP0a16N8rmLHvhdipek8Sg4FdEJKls3w6ffw4dOkCRIjZhrUULaN9ehWviVTNWhTJw7hb2h4WTLziI3o1K07LKlceq7dsHQ4bA2LFw5gw0amSdG+rXv3rDkW82fsPQJUNZfXA1p6NOUzZnWUY2GUm7iu3IFJgpwV+XyOUo+BURSUxhYbYV9umnsHChRQbFilnwe++9dhHxohmrQuk7fR3h0bEAhIaF03f6OoBLAuCNG21Y4KRJ4PHAQw9Z54bKla/8+B7Xw0/bf+Lt399m4d6FAPg6voy4ZwTdqndTHq8kOQW/IiKJ5exZK20/fRrKlIH33rO0hoIFvb0ykQsGzt1yIfA9Lzw6loFzt1wIfv/804rYZs2yHN4uXSxLp0iRKz9uWEQYn6z6hFHLR7Ht+DYy+mfEwbnQtuxU5CkFvuIVCn5FRBLKhg02r3XvXuvJmyGDVQBVrWrT1/RGL8nQ/rDwyx4PPRHOzJkW9P71F2TPbhOzn34acuS48uOtPbSWkUtH8sW6LzgXfY46BevwZv03yZ85P42/aExUbBQBvgHUL1I/UV6PyLUo+BURuRlhYda9/5NPYNky8POzVIboaPD3h06dvL1CkavKFxxE6L8CYDfGh7Mb83FueQlaDrDd3WHD4Ikn7PPc5UTHRvPt5m8ZsXQEf+z5g3R+6WhboS1PV3+aKnmrXLjdvPbz1L1BvE7Br4jI9YqNtUtAgPV36tHD+vEOGmRpDblyeXuFIvHWu1Fp+k5fx9kzDqdXF+L08qLEnklH0dJRvDMUHnjAPtP916K9i5i1dRaHzxzm+23fc+DMAYplLcYHd33A41UeJ1vQpY191b1BkgMFvyIi8fX335bW8Nlndv63Y0do1w7q1rXUBqU1SApUM09+yuzNzHdTgoiN9CNL8eM81/8M/brmuOyvtOu6jFsxjqe/f5pY13KFa+Wvxfh7x9O4RGN8fdS1RJI3Bb8iIlfj8Vinho8/tqofHx/r61SsmF2fNavl84qkMFu2WOeGzz+HmJhM3H8/9OkD1apdfhTbuehzTF43mRHLRrD64Op/XeNDicz1aVqqaZKsW+RmKfgVEfkv17Uxw0WLWrA7cqR1bHjvPdvpzX/l/qciyd2iRVbENnOmDRJ88knr3FC8+OVvv+PEDkYvG82EVRM4EXGCwpluIVvs/Zzw+Q6XGBz8+GNdLmaUDL1qb2BJo2Jjk10PcwW/IiLnHT0KX3wBEybYQIr9+yE4GH74wcrbldYgKZTHA99/b0HvH3/YCYtXXoFnnrl8irrH9fDz9p8ZsWwEc7bOwcfxoVWZVjxT4xle+jKa/VERBPjUJMJnHek8FcBT6qLWaJJGuC4cOwZ79lz5kj07rFvn7ZVeRMGviMjmzZbD++23EBUF1avb+KqAALs+Z05vrk7khkVFWTOSgQOtE1+hQvar3bEjZMx48W0X7V3Ej9t+5FTkKb7f9j1bj20lV4ZcvHLbK3Su1pn8mS2wPXByDgCBnjIEespcuP+VWqZJCubx2CbArl2wc6d93b374uA2/D//34OC7BetUCFo0gRKl/bGyq9Kwa+IpE1790JkJJQoYbsXP/9snfs7doSKFb29OpGbcvo0jB8PgwfbKOIKFSy396GHrAPff01aO4kOMzoQ48YAUC5nOSa1mkTrMq0J9Au86Lb/bY327+OSwrguHD4cF9ju3Hnx93v22Ceof8uTxwLbChWgadO4QLdwYfuaPXuyP0um4FdE0o7oaJg926KCuXPh/vth6lSbvnbgQNxOr0gKdeiQ9eQdOtzD2dM+pCt0jDId9vJG95zcV/XilIQYTwyztsxixLIR/Lrz1wvHfR1f2lZoyyMVHrnsc5xvjfbvqXBB/r70bpT8dvgE+3dv925L5dq2zb6ev+zYcenObc6c1ty5alVo3dq+L1rUvhYuDOnSeeFFJCwFvyKSNgwebOd+DxywgrWXX4bHH4+7XoGvpGB//w0ffmid+KKiXDKWPkye6tsIzHeSc8DL3x7EcaBllfwcPXeUCSsnMGr5KPac3EPBzAXpEtKFiasnEh0bfc3pa+fzegfO3cL+sHDyBQfRu1Fp5ft6U3i4Bbb/DW63bbPd29h/ja8OCrJuNSVKwN13W2B7PrgtUuTSfJhUyHFdN8meLCQkxF2+fHmSPZ+IpGHR0Vao1qSJdeh/7TVYtcomrt1zz+W79oukMMuWWRHbN9/Y57fHHoNlmRZx3O/4JbfNknkv1cot5ct1XxIZG8kdRe6ge43u3Fv6Xvx8/Fi0d5GmryVn0dGWjrB1q33a2bo17vs9ey6+bbZs1r6jeHELcs9/X7w45M2b7NMSEorjOCtc1w255LiCXxFJVXbuhI8+sr68Bw/CnDkWALtumvkHX1I317Wsnf79YcECyJIFunWDZ5+1dMyiL83BBSJ9NhHuswbwEOGzikjfTaT3T0/7iu15psYzlMtVzsuvRC7huvbv1ubNdvl3oLtzJ8TExN02OBhKlYq7lCxpl+LF7Tq5YvCrrQ8RSR2OHLEevD/9ZEFukya2y3v33Xa9Al9J4aKj4auvbKd37VrL3vngA3jqKcicOe52+YKD2HpqAUcD3gViwQFfT3aK+HZjVc93CE4X7K2XIOdFR1u+7ebNsGlTXLC7eTOcPBl3u6AgC2grVbI50+eD3FKlUkRhWXKl4FdEUq7du21M1d1322m+s2ctvaFjRyhY0NurE7msGatCrytf9uxZaz394Yd2drtMGfjkE3jkkYtT1V3XZfG+xfjmGMbRiFngeP65wiGr24zBTV9W4JvUzp2zgHbDhouD3G3bLAA+L18+uOUWaNvWvpYpYy3C8ue3QTuSoK4Z/DqO8zHQDDjsum75f44NBO4FooDtwOOu64Yl4jpFRIzHY+d8R42yrv25clkvJ19f694vkozNWBV6UaeE0LBw+k63AQD/DYCPHIERI+xy/DjUrWvfN216cTwUERPB1PVTGb50OCsOrCBzYGbq5m/KX/vn4rox+Dj+9Lq9lQrSElN4eFyQe/6ycaPt7p5PL/Xzs/zbW26Bli3t6y23WJCbJYtXl5/WxGfndyIwAvjsX8d+Bvq6rhvjOE5/oC/wYsIvT0TkX2bNgh49LPctVy546SVLbUhmozNFrmTg3C0XtQgDCI+OvWg62o4dMGiQpa2Hh0OLFtC7twW//7bv1D5GLxvN+JXjOXLuCGVylGFUk1G0q9SOjAEZVcCWGKKiLMhdt+7iQPe/QW6pUtYqrF07KFfOLiVKXL7JsiS5awa/ruv+7jhOkf8c++lfPy4G7k/gdYmI2JvJwoV2SrBYMdsdKVQI3nsP7rtP7ckkxbnSFLT9YeGsWmX5vF99ZZ/n2rWDF16wM+Dnua7LH3v+YPjS4Xy76VtcXO4tdS/da3SnQdEGOP/KAa1dsLaC3hvlunZGad06S7A+/3Xz5riiMz8/y7+tUgUeffTiIFf/NiVrCZHz+wQwNQEeR0TEnDoFX3wBo0fD+vXw3HPWp/e226y8XSSF+u90NNeFiN3ZiVxRiqr9IVMm6NXLfuXz5Yu734KdCxi2dBjrDq1j24ltZE2XlZ61e9KtejeKBBdJ8teRqpw+bcHtvwPddesgLCzuNucnmjVrZhMgK1Sw3V0FuSnSTQW/juP8D4gBJl3lNp2ATgCFChW6macTkbSgVy8YO9aqfKpWtbZlbdp4e1UiCeL8dLRzkR7ObcnDqSXFiTqUhaw5YunfHzp3vjj9c3fYbv736/+YtM7eZh0cXqr7Eq/e/irp/dN76VWkUK4L+/fD6tUXX7Zti7tN5swW2LZpExfkli+v1mGpzA0Hv47jPIYVwjV0r9Is2HXdccA4sD6/N/p8IpJKRUfDL79A48bWticmxlr6dO0K1aurlY8kqevtxHC97i6dnx9iMvDphEAijwcRlPMsT792gg9fzkpgoN3GdV1+2/0bw5YMY+aWmfz7LdbH8SFzYGYFvtcSHW2dYP4b6B47FnebEiWgcmXo0MEC3YoVbYdX/+akejcU/DqO0xgrcLvddd1zCbskEUkTDhyAceNsl/fAAViyBGrUgKFDvb0ySaOupxPD9Tp2zBqUDBsGR48GU6sW9OkDLVpkwMcnAwDnos/xxdovGL50OOsPryd7UHZerPsiNfLV4JHpjxAVG3XN0cNpUkSEpSmsXAkrVtjX9eshMtKuT5fOdm/vu8+C3cqVLdDNlMmbqxYvik+rs8lAfSCH4zj7gH5Yd4dA4Od/kusXu67bJRHXKSKpxYED0LMnfP217fI2bgzjx0O1at5emaRx8enEcL1277Z09fHjreVrs2YW9N56a9wG466wXYxcOpIJqyZwIuIElfNU5uPmH9OmfBuC/IMAmNd+njo3gP1HXLs2LshdscK6LZwvQsua1dKlnn02LtAtVUrjzOUi8en28PBlDk9IhLWISGp19qzNpC9XzhIaly6F7t0ttaFkSW+vTgS4eieG67V2rXVumDLFgty2ba1zQ/nydv1fe/5iwqoJ/H38bxbuXYiDQ6syrXi25rPULVj3oq4NkEY7N5w7Z6kKy5fHBbubNkHsPx9QcuSwD81Nm1rAW60aFC6stAW5Jn0UEpHEs327nev9+GPImdPaBKVPb7PqNbVIkpn/dmL49/Er+XeOcN4sQTTJUYGF03Py44+QMaO1pX7uubiBg2ejzvLmb28y8K+BuFgub/uK7Xmn4TsUyFwgMV5WyhAdbakKy5bFXdavjwt08+Sx4Pa+++xr1apQoIACXbkhCn5FJOEtXQpvvQVz5ljD0tat4emn496oFPhKMnS+E8O/Ux+C/H3p3aj0ZW9/Pkf4XGQs5/7Ow4olxVl0IJgs2WJ55x1funa1s/AQl9rw0aqPCIsIu/AYvo4vt+S4JW0Fvh4PbN16caC7erXl7oL9R6te3XJEqle3y7/7voncJAW/IpIwzp61XZrMmWHPHitge+UV6NJFb1ySIpzP641vt4f+s//m8LL8nFpajJgTGfDLepZsjdZRqu4xXn65Pq7rMv+f/rzfbfkOB4fWZVvToEgDnp/7fNopYDt0yP49WLzYvi5fbr28wc4EVasG3brFBbrFimlHVxKVc5UuZQkuJCTEXb58eZI9n4gkgV27YORI68f7/PPw2mtWfBIby4XeTSKpyIkTNn/l1Xcj8JxNR0CeMDLX3E76UgdxfMAlgv89cOyirg2dq3Wma/WuF3Z4U+3o4chIWLUqLtBdvNj+jQArOqtY0bq61KhhgW6ZMhpPLonGcZwVruuG/Pe4dn5F5Mb89pu1JZs503Zp7r/fOjeAvcmpulpSmb17YcgQ69B35gwElzpLYNXVUORPIn3Xcc6TnyifLZzz/5nOs09ftmvDeamigM11YefOiwPd1ashKsquL1gQataEZ56xr1Wr2k6viJfp3UlE4i86Gvz97ftBg2DhQnjpJevaUCAN5SxKmrJhAwwcCJMmWbzXpg307g07PRH0mP4Xe3xfBqL/ubVDnXxN6d/oxct2bUjRIiKs48Jff8VdDh2y69Knh5AQq+6rWdMu+RNuOIhIQlLwKyLXFhpqqQ3jx9vuTvHi1sUhWzYIunIlvEhK5brw55/Wrmz2bIvtunWzzJ4iRSA8Opzl637kVIb3ITIu8L2vRFemtx3pzaUnnIMHYdEiC3IXLrR2Y+d3dYsXh7vvhtq17VK+vM72SIqh31QRubLly61D/1dfWYV2y5ZxzeS1qyOpkMcD331nQe+iRdZK9o03rFlJ9uyw9+Re+v4yinErx3E8/DjFsxbnbHQYHtdDgG8AvW97NMHXlNgjlwF74Rs2WJC7cKEFvDt22HUBAbar26MH1KljwW7u3An7/CJJSMGviFze8eNQt64VrT3zjE1MKlrU26sSSRSRkfDFF5besGWL7e6OGAGPPw5BQS4L9y5k2K/DmL5pOi4uLW9pybM1nuW2wrexeN/iRCteS7SRy5GR9uH2jz9si3vhQggLs+ty57a//W7dLNitWlXFq5KqKPgVEXPypA2jWLHCooBs2WDWLKhVy9qXiaRCJ0/C2LFWyHbgAFSpApMnW/1mDBFMXT+VYUuHsfLASoLTBdOzdk+6Ve9GkeAiFx4jMYvXEmzk8smTtpX9xx92WbrUAmCAW26BBx6wmcu33mofclNTrrLIfyj4FUnrtm+H4cMt8D19GurVs1L2jBktp08kFdq/3wLeMWPs1/7OO+HTT+3rrK3fcfekIaw6sIqwyDDK5izLmKZjeLTio2QIyJC067zRkcuHDsHvv8cFu2vXWmqDr6/t5D79tP2t161r0xdF0hAFvyJp2ezZ0Ly5vSG2aWOV2tWqeXtVIolm0yb44AP4/HNrRf3AA9Cnj8WDS/YtodEXr/Lzjp8B8HF8GNp4KN1rdPda14Z4j1zev9/aD56/bN5sx9Ontxzd116zXd1atSBD0gbwIsmNgl+RtCQmBr75xjo0NG8O9etrCpukWv8uFMt4MjfpN5dj6YIg0qWDTp2gZ08oWDiarzd+TbePhrIkdAkBvgE4OLi4ODicjTrr1XZlVxq5/FqlTJae9NtvsGABbNtmV2bObEHu44/D7bdbVH++PaGIAAp+RdKG06dtAtvQobB7NzRrZsFvxozw5pveXp1IgpuxKpSXvlnH8c3ZObWkMpH7suETFMWDT51ixDuZIf0Rxq4Yy+gZo9l/ej8ls5VkWONhlMlZhuaTmyeb0cPn83o/m/I7Rdcv445Dm7jj4EYyvL3HbhAcbOkLXbpYsFu5slqOiVyDxhuLpHajRkHfvnDqlL1J9uoF994LPj7eXplIooiKgnJtN7J7fkGij2XCN/M5MlffScaKe8mWYy9Vyy3iy3VfEhkbyd3F76ZHzR40LtEYH8f+JpLF6OGDB2H+fPj1V7ucbzuWNasFuecvFStqPLDIFWi8sUhasnKlVWxnzWqNSu+5x4Le6tW9vTKRRHPqlM1hGTwYQkPL4p/zFJkfmohT7Fdcx+Gw70r2Rq/n7w3pebzy43Sv2Z2yOcte8jheGT18/LilL5wPeDdutOOZM1t60rPPwh132DAJfXAVuSna+RVJLTwe+P57+PBDexN97z0bPSySyh08CMOG2UmOkyctRjxYeDVH88ziSOAbQCw44OPJSkH/B1nV6z2yBmX17qJPn7YuDOd3dlevtrFy6dPbGZoGDeyFVKmiNAaRG6SdX5HUynWtTdkHH1iFd4EC1qn/qae8vTKRRPX33/Zr/+mnlurQurV1bshcdAvPzRnCpp2TwPmnUMx1yObey5Cmr3on8I2Ott66P/8Mv/wCS5ZYAWpAgA2SeOMNC3irV7djIpJoFPyKpFQREZAunTWj//pr+37SJOvdpOpuScWWLrXxw9OnW5zYoQP07Omy0+cn+i0Zyg/f/0CAbwDVct/GysN/4Lqx+Dj+9Lq9dcKPBb4S17W+ar/8YgHvggXWP9txbFTwCy9YU+E6daz7iogkGQW/IinNnj2W1DhxIqxaZXNYp0yx3EBNZZJUynVh7lzo39/iyOBgq+Ps2PUsPx38nBY/D2Xz0c3kyZiHN+q/QedqncmdMXfSFq8dOGDB7vnL/v12vEQJePRRC3bvuMOmJ4qI1yj4FUkp1qyxdIYpU+znhx+Ouy5LFu+sSSSRRUfD1Km207tuHeSsuogm7y7gkUalWXNsMdUmjScsIoyQfCF8ft/nPFjuQQJ849IGErV4LTzcpqjNnQs//QQbNtjxHDmgYUMLdu+80z6gikiyoeBXJCU4etRyAQMDrer7ueegUCFvr0ok0Zw9a62pBw2ykx3lysErY/5i4JEGfB8VyfezwAcfWpdtTY+aPahTsE7iD6NwXQtw5861y++/Q2Sk/V3WqwePPWbBbqVK6sggkowp+BVJjmJiYNo0+OsvGD7cdpK++cYmN2X1cpW6SCI6fBhGjICRI637V716MHREFCcLTqXnj32IjI3855YOLYp35qsHRiXugo4dsxSG87u7oaF2vEwZ6NoVGjWC226zLg0ikiIo+BVJTs6dg08+sRL2Xbvgllusd1OWLDaYQiSV2rHDuvR9/LFtprZsCR2fPcxKZyxdl4/i4MqD+Lm5sLctDw5+rNxSkhmrQhO2iC0mxjoxnN/dXbbMdnyDg21Xt1EjuxQsmHDPKSJJSsGvSHKxaBG0aAFHjkDt2jBkiCaxSaq3cqXl806bZu1s27eHZh3XMPPQUFr/aVPY7ilxD7t31+fMqbJE+Wwmwmcd6TwVwFOKgXO33Hzwe+gQ/Pgj/PCD7e6eOGF/dzVrQr9+FuyGhKjfrkgqob9kEW/at8869IeEWFLj7bdbTu+tt6pzg6RarmuZBAMG2NfMmaFX71jKNJ/NZ38PoeXcBaT3T0/HKh3pXrM7t+S4haIvzcEBAj1lCPSUufBY+8PCr38BsbHWL+2HH2wwzIoVdjx3bvsAes89cNddSjESSaUU/Ip4w+bN9s7/xRc2rnTlSosApk3z9spEEk1MjLWkHjDAuvTlzQtPvvsLe/INYdKxVez/eT+FshRiwJ0DeLLqkxcNo8gXHEToZQLdfMHx7JF75IilMXz/vX09ftx2d2vVgrfftoC3cmWdaRFJAxT8iiSlVavgrbdgxgwbStG5M/Tq5e1ViSSqc+csl/fDD+NS2d8bu53lWV7mo81fwS7wcXx4u8HbvFj3Rfx8Ln1r6t2oNH2nryM8OvbCsSB/X3o3Kn35J3Vd+3ubPRvmzInL3c2ZE5o1s2D37rvVc1ckDVLwK5LYXNe2vPz9rVHpggXwyivQvbu9EYukUkePWteG4cOtaUKt2i4d3/qNZb5DeHnrdzgH41J7HBx88Lls4AtcyOsdOHcL+8PCyRccRO9GpS/O9z13DubNs4B39mwbMuE4UKMGvP46NGkCVatqd1ckjXNc102yJwsJCXGXL1+eZM8n4lWxsTZ/9b334JFHbJxpdLSVsmfM6O3ViSSaXbusP++ECRaPNmkeSeVHp/D98SGsPria7EHZ6RLShRr5a9Dm6zZExUYR4BvAvPbzrn8gxZ49trM7ezb8+quN/c6UyYrUzu/w5sqVKK9TRJI3x3FWuK4b8t/j2vkVSWhRUZbL278/bN0KpUpB4cJ2nb+/XURSoTVrLJ936lTbcG3d/hDZG4/h692j+H7jYcrlLMf4e8fTtkJbgvyDmLEqlBJOf/ZGLadg+hAOHS0E1+ogdr5Y7fzu7tq1drxYMUsjatbM+u4GBFz9cUQkzVLwK5LQ2re3d/8qVayA7b77wNfX26sSSRSuC/Pn22e9nzYuIqDUAu7qlZ/MlRbw7Y5JRG2MomnJpjxX6zkaFm14YQrbjFWh/+TwFiMLxTh1CvpOXwdwaeuys2fh55/hu+8s4D1yxP6mbr3VRn43awalS6tDiojEi9IeRG7WiRM2kuqJJyB/fmuQHxZmxTR6M5ZU6nxWz4ABsHw5BFdYyOnWDYglCoBA30CerPok3Wt0p3SOS4vS6r7/62W7N+QPDmLhSw2sBeDs2TBzpvVDi4iwYS9Nmlj/68aN1YpMRK5KaQ8iCe3AARg8GEaPhjNnIE8eeOopa4wvkkqFh8Onn9oQwu3bodgtZ3ho0Cf8HPU6sREW+Do49KnbhzfvePOKj3NJf17XpcSxvdy9aDHM/J99iHRdSxnq1AmaN7d0BqUNichNUvArcr1c1wZRjB9vBWwPPQQvvQQVK3p7ZSKJ5sQJGDUKhg2Dw4eh4m27af78CH47PZ6pp05SPmd5zkSfIdYTS4BvAPeUuOeqj5cvOIgDx88QErqJu/5ezJ3bllD0xAG7slo1eOMNGzhRoYLOoIhIglLwKxJfoaGW1uA4VsLevj306QMlSnh7ZZLGzFgVevWWXwlo7147wTFuHJw961LrgUWUuX0wfxybzoZjDveXvZ/naz1PzQI1WbR3EQt2LaB+kfpX7toQEQHz5jFpyRdk/ul7sp07RaSvH4sKVeKzWq2o+exjNG5cPVFei4gIKOdX5NrWrIF33oFvvrGm+RUr2u6vdqPEC+IKxS4e9vBeqwoJGgCvX2+1ZF9+CR4nmtodv+ZkmSGsP7GU4HTBdKraiWdqPEPBLNdqzwCcOmWT1b791r6eOQOZM7O3zh2MzVaJGbnKkyV39kQN4kUk7VHOr8j1WrzYgt7Zs61v6IsvQr58dp0CX/GSgXO3XBT4AoRHxzJw7pabDhxdF/74w4rY5syBgLJzyddnMGcyrWBh5FFKOiUZ2WQk7Su1J2PANXpVHz5s3Rm+/dYK1qKirN/uww9bB5QGDSgYGMjbwNs3tWoRkeuj4Ffkck6dgrvusl6hb74JzzyjynJJFi4pFLvG8fjweKypwoAB9pkva8ktlHn9JTYxgz2AT5QPA+8aSM/aPfFxrjIdbfduC3anT4eFC+2Bixa1v5/77oPatdX2T0S8TsGvCNiW1w8/WAQwZgxkzmw7vtWqaRqbJCv5goMu2yIsX3DQdT9WZCR8/rmlN2zd6pK3zq+Uf3cw66PmcMrxxXEdXFwcHKJjoy8f+G7bZilB33wDy5bZsYoV4dVXLeCtWFFnSkQkWVHwK2mbxwMzZsDbb1s+b6FCsH+/Fbbdfru3Vydyid6NSl8257d3o0t76V7JyZP2GW/IEDh4NILCTSdT6PEh7IlcS4xfTvrV7kf1fNV5YNoDF0YP1y9SP+4BNm2yYPfrry0nHqB6dZt00aqVikBFJFlT8Ctp17Zt0LIlbNgAJUvCxx9D27YaiyrJ2vm83uvp9nC+O8SevR4860tyYkVBwjlK8YdGE1x8FLujD1M+S3km1JrAIxUeIZ1fuotHDweFEP7XGdjYzwLejRvtgevWhUGDLOA9P8JbRCSZU/AraUtMDOzYAaVKQcGCVsD2v//Bgw8qF1FSjJZV8se7uG3GqlB6jtvBgU2niPBshLCD+N8/B79i37LdjaRJkSY8X+v5S0cPf7OW4vs8dN4Szj1b3qfYif24Pj449erB8OGW0pA//0XPk1Tt10REboaCX0kboqLgs8/gvfdsMMW2bRAYCD/95O2ViSSav/6CJ54M4MQZH3jsIfCNAgeiXX9yO41Z0G0At+S4Je4Orgtr1nCq5wB+WL2AImEHiHF8WFSoIh/VuI+11eoz+537L3me/7ZfCw0Lp+/0dQAKgEUk2VHwK6lbRISlM/TvD3v2QEiIFeJoRKqkUh6PtSnr3x8WLgnHqTYNp+WruH42ehjXIXNMK4Ii2lng67rW1HfqVPjqK/j7bx5xfPircCVG1XqAn0vW5ET6LAA4sZd/zsRsvyYiktAU/ErqNm8ePP001KkDY8dCo0aqPJdUKSoKJk2yzg2b9h4ky50jyfDqGM66R/Hz5CfG9QM8OPiR3hNCnchD0K+fBbybN4OPD9xxB7zwAi325GBDTOAlz3GljhKJ0X5NRCSxKPiV1OXMGRg92t7Ie/WCJk3gzz8t+FXQK6nQqVM2enjIEAiNWUPWJoPxfehLThFD81LNqZ6jHZ//lp6TsRtIH/EnTbeG033pCEoe2W1/E7ffDj16WNFarlwAPHWFKXJX6iiRkO3XREQSm4JfSR1On4aRI+HDD+HoUWjd2o47jlWki6QyBw/C0KEwarSHU7m/J/j+wZD1V6L8M9C1cmd61OpBiWwlYNcuWn43loCvp1HywHY8jsPxKjWgX2/7O8mT55LHvt6OEgnRfk1EJKk4rute/QaO8zHQDDjsum75f45lA6YCRYBdwIOu65641pOFhIS4y5cvv8kli/zHN99Ap05w/Djccw+89hrUquXtVYkkiq1b4YMP4OP584kNGUJQsVWEB+ylQOYCdK/RnaeqPkXWk5EwbRpMngyLFtkda9WCNm3g/vsv6tKQUNTtQUSSG8dxVriuG3LJ8XgEv7cBZ4DP/hX8DgCOu677vuM4LwFZXdd98VqLUPArCebkSStmy50bVq60gPe116BGDW+vTCRRLFli44en/7wfn7v64qnwGTjg4NDv9n68XL4r/jNnWcA7f75VvlWsaAFvmzY2ZlhEJA254eD3nzsXAWb/K/jdAtR3XfeA4zh5gQWu617z/JaCX7lpYWF2rnfIEGjRAiZO9PKCRBLP+anbAwbAb1tW43/bYGLLTMbjRF+4jS8+vBVair4Tt1sbv+LF4eGH7VK2rBdXLyLiXVcKfm805ze367oHAP4JgHNd5Yk7AZ0AChUqdINPJ2ne8eMW8A4dahU+LVtakY5IIvPG6fzoaJgyBfoP8LAh8gcC6w+CO34lwD8DHSt15tbjGXls20Ci3FgCPB7qLzsKzzxjAW9IiIo7RUSuItEL3lzXHQeMA9v5Teznk1TqzTct8G3d2vr0Vqrk7RVJGpBUwxvOB9j7Dkfhs7UYp9fk4VieLwm4ewhk3kLOTAV4Ns/TPLUokuDBU+DoUQqUzcSCxmWof/tj1H6jsyYUiojEk9IeJHk6dgwGDbJWZXXrwoED1sWhQgVvr0zSkLrv/3rZFl75g4NY+FKDBHmOGatC6f35FvZvOso5z3IIOgxlvoN0J6iWtRy9jpXi/kmr8d++E9Klg+bNoW1baNwYAgISZA0iIqlRQqc9fAc8Brz/z9eZN7E2kTjng95hw+DsWciY0YLfvHntIpKEEnt4w/bt0KWry6HIPdD8yQsj1AqfzMNbPxfk0RUbcJyN0LAhvPKa9eLNnDlBnltEJK26ZvDrOM5koD6Qw3GcfUA/LOj9ynGcjsAe4IHEXKSkER98YOkNZ87Agw9aekO5ct5elaRhiTW8YcUKeL+/h2/WzMWt9SEUmwcu4ICvBzovP0jlfcVxPvzQOjXky3dTzyciInGuGfy6rvvwFa5qmMBrkbTo+HHIksXyFT0e69P76qtQvry3VyaSoMMbXBd+/hneHRDBbye+wKfuINxym8gWHkir5T5MquQhyhd8XF++qdCT2U0bs7BnwqRWiIhIHE14E+84fjwuvWHsWKtS791bVeqSrFzvpDO4tDtEz4alidqWn3eGHGZD0Gj8aw6D9MepeNSPXtPhvl0BzC7egG25i7I69xnSuRU4m708r2g6mohIolDwK0nr30HvmTPwwANQubJdp8BXkqGWVfLHu7PDv7tDeKJ92DQvF23GHcKv/EvE3D0F/GJotAV6LvWhftnGOK8+Bs2aEbjpGFFztxCs6WgiIolOwa8krcaNYdmyuJxepTdIKjJw7hbOnPTh+PYwIgO/J0vhxcTc+jd+0dBxNTx3uhy3tOoEQ9pArrj26NcTYIuIyM1R8CuJKywMRo60gRQZM1pRW9asalkmqc6uXbB2Wn6yBL9HeL1PcR047sIjawKodLQRfca8r4lrIiLJgIJfSRwnT9pQikGD7PtbbrEBFbfd5u2ViSSo1ath+Fsb2HG4F+nrzyM0c4x1bgAcHH4o8xC7MnWgjwJfEZFkQcGvJKzYWHj/ffjwQzhxAlq0gNdfj8vrFUkFXBd+mxvBrDfH8nfwIOZV2cO5inDbocxUCKvMrwUX47oxgB+ZfarcUHcIEUmbvDFSPa1R8CsJIyYG/PysZdm8eTaY4vXXoVo1b69MJMHExrjMH7iMxZPfZ0mZ2cy5Oxo/DzzkqUjvJu9QsXozZqwK5fAP09l7bjkF04fw5j2t9MYlIvGSVCPV07p4jTdOKBpvnAqdPQujRln3hsWLIX9+CA+HoJsbAiCSnCxYOpsvPx5Bhs0rWVziCIsLQpaoADoXbkWPhwaSL0sBby9RRFKBpBipnpYk9HhjSevCw2HMGEtxOHwYGjWCc+fsOgW+khpER3N66vd8M+59OtVfTHQeIC/kjMnC0Nv+R8e63cgQkMHbqxSRVCSxR6rfjNSUjqHgV67f2bNWwLZvHzRsaCOJ69Tx9qpEEsbWrZwaPIH90z9hYsUjDLkNon3tKh986HFXb569rbd31ygiqVJijVS/WaktHcPH2wuQFCIqCn74wb7PkAGefRZ++w1++UWBr6QYM1aFUvf9Xyn60hzqvv8rM1aF2hVnz8LEiZytWo81t5Wm26GBlO98lAF1faiY53YCfQPxdXwJ9AukQVGdehSRxNG7UWmC/H0vOnajI9UT0sC5Wy4a8w4QHh3LwLlbvLSim6OdX7m66Gj47DN46y3YvRs2bLBepb218yUpyyU7FyfOMXn4NCqFLaPA3O+Yl+8Mb9UO4q8W4O8Jon25J3nlrh4Uy1qMRXsXsWDXAuoXqU/tgrW9/EpEJLW6kZHqSSE5p2PcCAW/cnkxMfDll5bSsH07VK8OY8dCmTLeXpnIDTm/c5H13ElabZjPg2t/osiJPXxSMYC3O6QnNBdkIphXa7/G8/U6kzUo64X71i5YW0GviCSJ5DjxMbmmY9woBb9yeSdPwtNPQ8mSMGsWNG0KjuPtVYncGNel0JolvLjmR7KcXcgvxWPpWysHP5bOTHimU+TzLc34RkNpX7UNAb4B3l6tiEiy0rtR6YvOnEHySMe4UQp+xXg8MH06zJxpaQ7Zs8OyZVC6tIJeSbkOHYJPP4Xx45m8bRsTywdx32MePD6AcxTfM+Up49uODf/rjaPfcxGRODExsH8/hIXR8sxJ8mQ/yNxFW/g9cxEiipfk1cqZaDxtJHx6zuomzp2zS8+eUK+etT997DH4/XfIndvbr+YiCn7TOte1nd3XXoM1ayyt4dAhyJPHOjqIpDQejxVijh8PM2ZATAzfVKlErzYh7C79T59xB3B9yBYUwrtN2yrwFZHUKTYWjh+HI0fg6FH7WqwYVKliU1iffRbCwuxs78mT9n2fPnbmd+dOKFXqwkPV+ufCqFHQtYHNdh8wwIrgM2SA9Ontcvq03SFLFnueZPjvq4LftGzXLnjwQdvhLVECvvgC2rSxKW0iSeyme0ju3w+ffAITJsDOnURmzcbbd9zDoEL7OFdwFQGxWQnJ/AArTs/EdWPwcfzpeft9yS63TkTkqlw3LqCcMwdCQ+3fvwMHrO/+7bfDc89ZwXq6dLYh8G89e1pQGhAAf/4JwcF2KVrUvhYpYrfLlw8++siOZckSd8mb166vVMl2h6+kTBmYMiUhX3mCUfCbFh0+DLly2S9wunQWLLRvb+OJRbzghntIxsbC3LkwbhzMng2xseyvfjtPV67PjAK/Q/ZZZI4tSr/qw3jhzsfJGJBRnRtEJHlyXdt9PXUKChWyYx98AJs3xwW3Bw5YSsG0aXb944/bbi5AzpyWXlC1qv3s729F61myQI4cdsmZEwr8M5EyQwbb3b2SDBmgY8crX58Md3TjS+ON05KFC+HVV+Hvv+2SLp23VyQC3MBIz4MH4eOPLejdvRs3Vy4W13+Qjr4umwpOgfTHyE8N3mj0Ao/VuA8/H32wExEvi4qyXdoTJ+IC1HffhT/+gD177HLmjHVXWrrUrq9RwwZK5c1rO7F580LNmnFB6fr1tjObO7cFu3IRjTdOy5Yvt6D3xx9tx7dv3xT9iU1Sn3j1kHRdWLDAxmpPnw4xMXjuaMhr97Xkw9g/icgyFnxjqBDQnA9bv8Cdpeoql1dEko7r2i7s3r1QrZodGzLE/r3ascN2b13XAtj9++36v/+2+5QuDXfdZTu+pf/VQWHJkqu/X5cvn2gvJzVT8JvaLVkCtWpZ94b+/S2JPUMGb69K5CJX7SF5/Lh1bBgzBrZuhaxZiejcndcyVGXIyTFEBw8FF3wcXz5vOZlHKj3khVcgImlCbKzVyxQtCj4+8M03Vi+zfbsFuGfPWrAaEWE5tSdP2v3uvNNyaQsXtst5n3xy9efTB/hEoeA3Ndq0yU6FPPCAnTIZOxYefhgyZfL2ykQu65Iekq5L7cN/03/DX/D6LIiMhDp1OPp0XzqF+jHz1DA86QfjE5gecMBxcRzYfWqHV1+HiKQC59NBHcc6GkyZAlu22GX7dktf2LnTgtkDB2z3tlgxaNjQvhYrFvdY/frZRZIVBb+pyY4d8MYb9ik0d25o0cI+eXbq5O2ViVzV+aK2kd+tosZfP9Bh3VxKHthuH9g6dmTDne148velLN75BgTvIqNvSbpXGMPdVUrT5MsmRMVGEeAbQP0i9b37QkQkZTl+3NKpNm2yM0tbttjXmTOtsGzLFhg0yDoilSoF995raQlZstj9n3nGLpKiKPhNDfbvt6D344+tY0PPntanL0CTqsR7rqt12YYNtPxoFC0/+8wKPipXhtfH8n3BO3h+1idsXXIPBIeRJ6ou/W4fTKfbm+Pj+AAwr/08dW8QkSuLibEd2w0bYONG+9qhAzRqZJ0UWre22+XPbwHuQw9B1n/Gm993nw1uUDekVEX/N1OD81OsunSBl1+O68En4iXxal0WHW1DKEaNsp2XwEAWtW/A/PpFOBFdl8//+pVDe7tDrmhuce/jg/teoGnFS4Pb2gVrK+gVEetnu2MHrF1rnRFq1bK0hCJFLFXhvKJFoUkT+75KFeuscMstl08N1CZSqqTgNyU6dsymqpw6BaNH2x/vvn3Ww08kGRg4d8tFM+ABwqNjGTh3Cy1zOzZ9bdw4O2tRuDC8/z6/312Bht+1ImZrJDijIVcAt6Z/khGPPk+lgiW89EpEJFmKibHdWNe1tIOVK2HdOis4A+t/W6uWTSvt1ctSFcqVs8EL/y76Dgqy1mKSpij4TUlOnrTco8GD7dRw+/b2SdfHR4GvJCuXtC5zXWruXU/7GXPg1cX2xtW4MYwdy4k6d9HjoxlM+vIpPBkjbfQwDn1v6827d73tjeWLSHKyYwesWgVr1tiu7po1Fsz++KMVpa1cCYGB8MQTNnWsYkUoW9bu6zjWS1fkXxT8phRz51rHhhMn4P774fXX7VOsSDJ0vnVZhshz3LdhPu1WzaH00T2cCsoEPXpA167sCMxN13Ef88uvt+DJsosAv4J4HH9cPAT4BnDvLU29/TJEJCnFxFiB2cqVdlboxRft+JNPwvz5ttFTqpQNeahXL+5+ixZ5Z72SYin4Tc4iIqz5dcGC1si6fn0bVlGlirdXJnJVr9/iz8F3h9Jyzc9kigpnbZ4SvNzseWr37UaOTBl5esxwVvmOhqATZPerS58qg3nh3uYs2bdExWsiaUFUlE0kcxyYONFacq5ZA+H/nDXKnNmKt/394f33LfAtV87SFERuksYbJ0dRUdb4+q23rL3K/PlqdC3Jn8djpyGHD4cff8Tj589P5W9nbIV7OFymEnXyBPLNmonszPQ5+EZTJOI+3m/+Ag/VUZArkqrFxFiXhaVLYdkyWLHC8nN37rTCtFGjYOpUG/l7/lK6tDosyE3TeOOUIDYWJk2ylIadO6FOHWuOrcBXkrOTJ+3D2siRsG2bdRt54w18OnUiU9QOIudMZP3W6fwVtRAypaMKTzC8TU/q3lLS2ysXkYTmujYBbelSew8rWNCGRLRrZ9cHB9vo3+eei7tPt252EUkiCn6Tk3Hj7B+AqlUtkGjcWIGvJF8bN8KIEfDZZ1ZhXaeOna1o1YpIfGk39n2mHXsNHA9kgeq+HZnS+T2K5cnp7ZWLSEI6etT+LTi/s3v0qB0fO9aGLDVsaMOXatSws5kp/H3tunqYS7Kk4NebXBdmz7acpsaNrXtDnjzQsmWK/8dBUqnYWJgzx1IbfvnFKqzbtIHu3aFaNY6EnaPbkI+YcWgQMZm3X7ibr48v991ePMEDX70JiSShmBjrtrBokV3q1YPOnS0f9+23rY1Y8+YW5NaoYbUqYGeD2rb17toTSLx6mEuyp+DXG1wX5s2DV16BJUtsykzjxtZ78L77vL06kUudPAkTJtjuzs6dUKAAvPMOPPUU5MzJpj1H6PL66/wRMRI36CgZPTVpmfcJvjnydqKNHtabkEgii4iAdOnsPatJE/j9d5t2BhbQnu84lC2b9Z1Pn957a00iV+1hrn93UgwFv0lt2TIbPbxggeVCjR8Pjz3m7VWJXN7ff9su7yefWG/pW2+1ASstW4KfH/PXbOfZD19nvd8n4B9O7sh7ee3W3nRteiuO49Bt7x2J1r1Bb0IiCSg2Ftavh4UL43Z2c+a0r44DhQpZy7Hate1SqNDFZyjTQOALl+lhfo3jkjwp+E0qrmv/UPz9N2zaBMOGWS5UYKC3VyZyMdeFX3+FIUMsxcHPD9q0YX6jh+mxext7V88iw5JVRERu4ljWmeDvR8mIR/mwWS/urVX2oodKzNHDehMSuQnh4dZarFYt+7ltW+u4AJA7twW4t98ed/uxY5N+jcnQ+R7mlzsuKYeC38S2cSO89pqNT3zxRcuPbNHi4vGKIslBeDh8+aUFvevX267Pq69Cly7MOOjh2elT2OvTF/yiCfMH/DNQ8uRzTHn6BaqWzJvky9WbkMh1CAuDP/6AP/+0r8uXQ3Q0HDoEuXJBx47QrBnUrQtFiqju5Ap6Nyp9UboVQJC/L70blfbiquR6+Xh7AanW9u1WwFa+PPz0kxW1gRUGKPCV5GT/fss/L1jQTmv6+lqaw5498MYbnM2ajS5jP2ZvVH/wibbxw65DZk9LcuZp6pXAF+xNKMjf96JjehMS+ceePdZhITTUfv7mGytGGzzYAtvnn4fvvoNMmez6u+6CRx+FokUV+F5Fyyr5ea9VBfIHB+EA+YODeK9VBaVapTDa+U0Mw4ZBr14W8L7wguX45sjh7VWJXGz5ctvlnTrV8v1atLDRw7ffDo7DgRMn6TZiOLOODCE27344UQzShYETi4Mf6X0qezXF4Pybjbo9iACnT8P06VZP8ttvVpgKVqj6xBNw7712vHp1TUm7SS2r5Ne/Mymcgt+EcviwfVrOmdMaeHfuDP/7n1XEiiQXsbEwezZH33yPHCuXcDogiO9rtiDriz25+946AKzfHUqXiUP5K2oMbsBpMp9rSOYTr+JkLUB09GYifNaRzlOBQE8Zr6cY6E1I0iTXteB2wQKbkNa4saUtdehgnRduuy3ug2yFCnafXLnsIiIKfm/aiRMwcCAMHWpdG0aNspypunW9vTKROGfPwsSJttO7bRuRWXLxVoMnmVrxbs4Epido6SkejFzAtJUT2eT3JTix5Dv9IG807k3He6oyc/X5tmJlCPSUAZRiIJLkJk60/tq//Qb79tmxhx+24DdXLiumLlXK0utE5IoU/N6o06ct4P3gA+tv2KbNxeMaRZKD/futN++YMfZBrWZNXqnxCJPzV+Oc31bCfWbjOZWNc1FreXPDr+Ck55YzXRj80PM0rln0wsMoxUAkiYWGWteV0FB46SU7NnYs7NhhO7r169vXsv/qsHLLLV5ZqkhK47ium2RPFhIS4i5fvjzJni9RdesGo0dbnuRbb8WdWhJJDlavhkGDYMoUS3W47z7o2RPq1KHoS3M452zgcMAr4PxTwBaehcA9j7P0jVeoWDK7t1cvkjb98QdMnmxDkLZutWP58sHu3dZy8PhxyJpVBWki8eQ4zgrXdUP+e1znRuIrKsqC3dWr7ecXX7TpbDNmKPCV5MHjsb68DRtClSpW/NK1q/WW/vprqFOHMxHhhJ9ezGHev7hzA82pVvFeBb4iSeX0aft77dXLzsqADZj4/HMoUQI+/BBWrYK9ey3wBcvnVeArctOU9nAtsbHWLuaNN6zAoE8fqFwZChe2i4i3RUTY7+iHH8LmzTZ6eMAAGz0cHAzA/hPH6frxKOYcHUZsriNwrDQEnLnQuSGrfxXl74okttBQ677w88+weDHExNigo5YtoV496N49rlOQiCQaBb9XM3Om5Vpt3gxVq8LIkVZYIJIcHD9uZyOGD7dG9VWqWBD84IMX3jzX7dlDl08HsSjyI1z/s2Q52YTnb+lNuSYlePPnb9l7bjkF04fw5j2tlL8rkpBcF7Zts0C3fHnrwBAWBq+/bu8nL7wAd94JderEtR5TD3iRJKHg97/O50A7DqxYYVWz33xjOZM63STJwa5d1qh+wgTr4tC4MX+26MCLJ3Kwf10E+fb+wZ0lPUxbNZFNvlPAdcgX9ghvNn6BJ5pWuPBrfH/17l59GSKpjutaitHPP9tl1y47/txzFvyWLQtHjkB2pReJeNNNFbw5jvM88CTgAuuAx13XjbjS7ZN9wduvv9qkq5deskk4ERG2g+bre+37iiS2lSutrd60afZB7JFH4IUXmBGTjb7T13Eidj2noucTGbkfT/bVEJmRW851YshDz9GodkFvr14k9YmJsfSF0FB46CE7VqKEBbgNGtjUtLvusmPaPBFJclcqeLvhnV/HcfIDzwJlXdcNdxznK6ANMPGGV+ktixfbQIpff7V8ydh/ZnanS+fddYm4Lsyda0Hvr7/aKNLnn7cG9gUKAND/3Z/ZGzads9k/gkAXMoDfrkcom6sda4YoTUckQe3bBz/+aJdffoGTJyF3bks3chzb8S1YMK5ITUSSnZvt9uAHBDmO4wekB/bf/JKSWJcuULs2rF9vAwD+/ttSHES8KToaPvsMKlWCe+6BLVssAN67174WKMDZyHC6ThjD4qOdOZtjPHYCBsCHjPn9OOXEevMViKQOEREW0MbE2M+DBlkx6ZIl8MADluaweXPczm7Rogp8RZK5G/4LdV031HGcD4A9QDjwk+u6PyXYyhLT5s32D1RgoFXYFi5sVbYZM3p7ZZLWnTkDH31knRv27bNCmU8/tSEqAQEAHAg7TrdPRjHr8DBi0x3BCQ8h3YnmRBQejevG4OBHOk8Fr48eFkmx/v47bnd3/nwbHfznnza5s3t36NjR8neVyiCSIt1M2kNWoAVQFAgDpjmO86jrul/853adgE4AhQoVuvGVJoQdO+DNN62P4pAh9o9Y27beXZMIWI7g8OE2je3ECZvcNG6cdRf55w12/d49dP10MAsjxuP6nyVzWBN6VO9DxXuK8+qs9YRFFSDCZx3pPBUI9i2v1mUi8RUZCefO2QCJP/+0TRGAkiXhySft77BKFTtWtOiVH0dEUoSbOTdzJ7DTdd0jAI7jTAfqABcFv67rjgPGgRW83cTz3bjQUHj7bdtR8/OznMk2bbyyFJGL7Nplu7wTJtju0n332QCVmjUv3GTe+nX0mDqQDc5kAPKceJg37urNU6/FdW7w83cYODeA/WFlNHpYJD727YPvv7dBE/Pm2dTOAQOgRg37ENq4MRQv7u1VikgiuJngdw9Qy3Gc9FjaQ0MgebZyePhhK2p76ikrbMuvoEC8bO1ae6OdMsXa6bVrB717wy23APDXnr8Y/NtEftu0niOBi8CTgVKnn2HQg8/T9NZLz6C0rJJfwa5IfNWvD7/9Zt8XLgyPPQb33ms/BwTA0097bWkikvhuJud3ieM4XwMrgRhgFf/s8CY7w4dD5sw6XSXe5brwxx/Qv7/tOGXMaP0/n3vuQucGj+vh6S/7M+bvV8DxQACUOfcUnz/2PtXKZvPq8kVSnOPH4YcfYNYs2LMH/vrLjt95JzRrBk2aQJkyyt0VSWNuqiTVdd1+QL8EWkviqVTJ2yuQtMzjsTff/v1h0SLImdPScLp1sxxDIDwqkj5ffs6EjQMJz7D1wl19fXxp16yoAl+R6zF7tnVFWbjQWlfmzm2BbmSkFTq/8oq3VygiXnSzrc5E5EpiYmDSJKhYEVq2hAMHbET27t2WfpM1K0dOneTBoQPI/FpRRux+ipjwDNwT+CZB/kH4Or4E+AZQv0h9b78SkeQrOhoWLIBevayoGWyMcFiYDSxasgT274ePP7bAV0TSPDUjFEloERHWnqx/f9i5k1PFSzP4wb58XrgWuU9lpPfm45TL40PnT4ay4Oxo3IBTZDjdkJ7FPuWtvncSEOCwaO+dLNi1gPpF6lO7YG1vvyKR5OXcOZg5086o/PCDBboBAVCnDhQrZl18Hn3U26sUkWTqpsYbX69kP95Y5GacOQNjx1r3hgMHoEYNFj/chSeO5uaEZyMRPutwwgty7tQ2IvN8Cz4x5Dx6P6/c0Yfurasp7VDkanbvtiC3UiU4dgxy5YLs2aFpUytWu+sum4AoydaMVaEMnLuF/WHh6kojSSLBxxuLyD+OHbOiymHDrEdvw4bwxRdwxx306j+fE56VHAx4GYiGTEB6fwJ3P8Kkx1+h9R0lvL16keTJdWH1atvhnTnTvm/QwNqSZc8Oa9ZYsZqvr7dXKvEwY1UofaevIzzaJk+GhoXTd/o6AAXAkuQU/IrcqP37bdTpmDFw9iy0aAF9+17o0eu6Lpt3ruJEphEQGG33cR0ynnuAHPkeUOAr8l8ej7X+A3jwQRsd7Dg2WW3gQPsbO698ee+sUW7IwLlbLgS+54VHxzJw7hYFv5LkFPyKXK+dOy2f95NPrKjt4YetsOafN+Po2Bhenfw1I1f350y+1XA2B3j8wPHg4EfGwCrky3zl0cM6NShpypkzMHeu7e7+9JONnw8Otpzde+6xlmS5cnl7lXKT9oeFX9dxkcSk4FckvrZsgffes5QGX194/HHo08cKbIBT58J5duInTN71AVEZduIXW5q64SM4RFHORG+M1+hhnRqUNGPtWnj1VQt8IyMhWzbL3T192oLff+/ySoqXLziI0MsEuvmCr7wRIJJYFPyKXMv69fDOOzB1KqRLB927wwsvXJgUuO/YCbp8NIofTgzFE3SEoIiadC8yiAEvNSddoM8/O7np4jV6WKcGJdXavx9mzLAzJLfdZm3HVq2Czp2hVStLbfDTW1Jq1btR6Ys+2AME+ftecSNAJDHpXxqRK1m50oZRfPutTWPr0wd69mRR5HYW7PiMAqHlGPX97yyOHgsBZ8h69h56V3qRPg/ehq9vXOuG6xk9rFODkqps325/P9On24AXsImGt90GpUtbBwe1OUkTzv8bqJQuSQ4U/Ir81+LF8NZbNoI4SxZ47TV49lnInp1Fexdxx8QGRMZGAi64PhQ404Z3mvah/d03P0lQpwYlRXNdOHgQ8ua17xs1sgC4alX7INmqlXVoOE+Bb5pyPRsBIolJwa8I2Bv1779b0Hu+ldI778DTT1sADHw2byld53YkMn0EOIDr0CBXZ+a9MSrBlqFTg5LiuK61HZs2zbozHD1qAbC/v01VK1QIihTx9ipFRC5Q8Ctpm+vCL7/Am2/Cn39C7tzwwQeWh5gxIx6PS/9pPzHwr/c5ETwf/DKC6wdY54btB0syY1Vogu1m6NSgpCizZsHzz9vurq8v1K9vP8fGWvB7223eXqGIyCUU/Era5LpWZf7GG5bmUKCADaro2BGCgoiMiqXPuKmM39Sf8OBV+PjnI+uuvqQPrk4s2y90bsBTKsGL0XRqUJIl14WlS21396GHICQEsmaF4sWt1V+LFpAzp7dXKWmcWkVKfCj4lbTFdeHHHy3oXbLETsmOGQMdOkBgIMdPRfD0yLF8fWAgMZm3E+BTmg5ZJzD0hbZUeusXXMDPU4ZAT1zeoorRJNXyeOzv5HxKw969tqNbvLgFv7feah8iRZIBtYqU+FLwK6nKFT/1u64VsL3xBixbBoULw9ix0KEDiw6tYPq8d1mw7DArzn6Lm+EQGTzV6VVkIG8/2gI/X5s4pWI0SRP+XbQWEwNNmsC5c3D33Va0du+9tuMrksyoVaTEl4JfSTUu+6n/m7Xk+e1nak0aCcuXW+HN+PHQvj0EBDDhj9k8Ne8+XGLAgWC3Bq9XmUz3ZvXx8bm4El3FaJJqua619ps6Fb76ynrwbt4MAQEwZw6UK3eh8FMkuVKrSIkvBb+Salz0qd91uWvbEp5dOJkKh7bbFLYJE6BdO/D3Z/Zf23l+2gdsyzQefGLBAR986dO8JT3q3XHZx1cxmqRKkyfbpLXt223IxN13w4MPWsqDry/UqePtFYrEi87OSXwp+JVUY39YOLgud/+9mB4LJ1Pu8A52BeflhSbP8cGMAbh+/oyftZp+P/fnYLavIKMfxTzNCA2YS4wnmgDfAOoXqX/V51AxmqR4GzbYDu/jj0PRorbLW7w4vPwytGxpY4ZFUiCdnZP4UvArqYPr8uD+VbSf+wnlDu9gR9Z89Gz6PDPL1idXlgy8MekvBi17n1O5fsTJnIl6fi8w9vHnKFMgL4v2LmLBrgXUL1Kf2gVre/uViCS8bdtsh3fKFNi4EXx8oFQpC35btbKLSAqns3MSX47rukn2ZCEhIe7y5cuT7PkkDThfyNavH6xYwe6seRlapw0zy9YnKtaX0wd2ci54EtG5F+MbkZNmOZ5jTMdu5AkO9vbKRRJXdLR1Zjh50lqQRUdDvXrQpg20bm09rUVEUjHHcVa4rhvy3+Pa+ZWU6Xyf3n79rPdo0aLwySesLncHP8z+jr1HRhEduBlK7ybgXBGezDuSwe0fJ2M65X5JKnbsmLUkmzzZRgfPn2+Fal9+CTVrQsGC3l6hiIjXKfiVlOX8RLZ+/WDRImtZ9tFH0L49W/dF88b4V9gcMARyu+A6PFb4dca3+x/+vvpVl1Tsp59g6FD7GhMDpUtD27b29+I4cP/93l6hiEiy4ePtBYjE2/z5Ni717rth3z7r07t1KwurtKLKCwMoPboIqwIHA5bK4+vjQ+niAQp8JfWJjIQZMyAszH7evBnWrbPRwitXwqZN1sHBca72KCIiaZKCX0n+fv8d7rgDGjSAnTth5EjcrX/zbbF7Kf7sK9z6dWFWZ3uFIgHVebvWaIL8g/B1fOPVvUEkxfB47G+hUyfIkwfuuw9mzrTrunSBXbtgwACoUkVBr4jIVWhLTJKvxYtt9+qXX+zNftgwPB2fYuysfbz+zLMczjsRcsdQwechhrd5kdtLVwKgQdlK6t4gqcuJE1C5MuzZA+nTW+Dbti3ceaddHxDg1eWJiKQkCn4l+Vm1Cl57DWbPtir1QYOIfLwL73y5hSHPdOB0wWmQz49bMzzO6Ed7Uz5/8YvuXrtgbQW9krLt3WtFaydPwjvv2Djhli2hRg1o0QIyZvT2CkVEUiwFv5J8bNxohWxff21v9u+9x093hfDKD1+xqm8DYvIsxidfJu7N8QIjH32Oglnz3vRTzlgVqp6QkjyEhdnv/qRJ8NtvVqzWoEFc0drQod5eoYhIqqDgV5LMFQPN7dvh9dftTT9DBnjtNfY/9ByPfDGB375rBD4eyA1Nc3fm88feJ2tQcIKt59/TgELDwuk7fR2AAmBJGtHRNnDC19fydd97z4ZPvP46PPIIlCjh7RWKiKQ6KniTJHE+0AwNC8fFAs1hn8xjV6u21pbpm2+gd282/fg39Y/fQv4Rt/NbYG9wPAD4+vhSt1zhBAt8waYA/XsMJkB4dCwD525JsOcQuYTrwrJl8OyzkC+ftScD6NrVjm/ebGk/CnxFRBKFdn4lSfw70Mx55jjdFk/jkdU/4OBAt24sbPg8T3/3E2u+qgs5dpA1pgxtK7zKhE0fEBUblSidG/aHhV/XcZGbEh4OgwfD559bgBsYCM2bW1472AAKDaEQEUl0Cn4lSewPCyc4/BRdFn/NYyvn4B8bzVcV7qJ/wcc5dXwPu/6oA4UOkt+twTuNP6Rdjeb4OD48sveeROvckC84iNDLBLr5gjUFThLIqVOwZQtUr24dGcaMgSJFoGdPeOAB0JhtEZEkp+BXbsh1FYqdPs3LK77mod+/ImNUODPL3M4b+Z9kXcAaPOW7QFAYpf3u5MNWk2hyyx04/+pRmpidG3o3Kn1Rzi9AkL8vvRuVTpTnkzQiNtba8332GXz7LWTObENZ/PxgwwbIlMnbKxQRSdMU/Mp1i3ehWESE7XS9+y5PHTnCzyXq0L30bWzJtwBydQT/SEo4TfnssX7ULlw9yV/H+bWq24MkmGnTbMpaaKh1LHnsMWjf3graQIGviEgyoOBXrtvVCsVaVskPMTG26/X667B3L5G3NuDd5l14O+pTPMXe/+cePjxfZQSDmj+d5Ov/t5ZV8ivYlRsXFgZffWVjt2+5BXLksGEUQ4bAvfdaXq+IiCQr6vYg1+1KBWEHTpy1na/y5aFjR8Kz5qXn/aNJnz87bxZ4CLfoXOtX6oCvj0POrKeSeOUiCSA2FubOhYcftsmDnTtbegPYGO7Zs+H++xX4iogkU9r5let2SaGY63L7zpX0XfgFDPibc0XL0rPx24zN/ieU7EqAJzMdy/elaaU6PPDVA4nWvUEk0Xk8UK6cFbFlywZPPQUdOkDVqt5emYiIxJOCX7lu/y4Uq7pvEy/+/ik1967neLaCPFa/J58VWwyFXiG9m5Pna75H7/pdyZIuCwDz2s9LtO4NIgnu1CmYOhUWLoSJE20gRbdu1p9XaQ0iIimS47pukj1ZSEiIu3z58iR7Pkk886bNI7Dfq9y6aREH02WjW+V7+bbyKsizlqxOIV6u35tutZ8gvX/6G34OjR4Wr3Bd+OMPmDDB0njCw6FsWTuWLZu3VyciIvHkOM4K13VD/ntcO79yffbsgX79aPDpp8wrEUT5tmXZkOskZPmUvP638NbdE2lf5RH8ff1v6mk0eli85quvoE0b68zQrh088QTUqGH56iIikuIp+JX4OXYM3n0Xd+RITjouLZpU4feQleBsxMHhrTvepW+9F/FxEqaG8podJUQSQlSUFah9/DHcdRf06AHNmsGnn0Lr1pAhg7dXKCIiCUzBr1zd2bMwZAix7w/gZMxpnqpVmRk1d+FJv/LCTXwcH3wcEizwBY0elkS2YYMFvJ9/DkeOWA5vkyZ2XYYM1ptXRERSJQW/cnnR0TBhAtGvvMHxiIM8XasE39aIxRO4ivr5mvFw1Xt57sfnEq1zg0YPS4KLirIRwwDPPQe//QbNm1taQ6NGcYMoREQkVVPwKxdzXdyvpnGu1yscPfM3PevkZUbVAFy/HTQv8SBv3PkSlfJUAqBCrgqJ1rlBo4clQbguLFsG48fD11/Dxo2QNy8MHw7Zs0POnN5eoYiIJDEFvylAUnU9+HP6cOZ+2p+yO0KZGhLMdxV98fE9yqPlH+OV+n0omb3kRbevXbB2orUr0+hhuSmnTtmUwfHjYe1aCAqChx6y3V+waWwiIpImqdVZMvffrgdgO6DvtaqQYIHguSVr+fGFTjxy+xIi//k45OcE0DWkG33q9aJA5gIJ8jwiicp1bdxw1qywbx8ULmyjhp96yqaxZcni7RWKiEgSUquzFOpGux7EZ7f4+Jq97O7wGqeOT6RPcx8i/+lO5uDwYr3evN3g7QR/PSIJ7vBh687w0UdQtCj8+CMUKABbt0Lx4t5enYiIJDMKfpO5G+l6cK0eubvXnmTz4+8RcWoQA+vFsLAwZPLNhL97Do/rIcA3gKYlmyb8ixFJSIsXw5AhMH26FWjWrQuPPBJ3vQJfERG5jJsKfh3HCQY+AsoDLvCE67qLEmBd8o8b6Xpwpd3idz7aSdjir3DDX2NEvTOszAd5A/IxrMFLdKzakTUH12j0sCRvR4/a8InAQPj9d/jpJ3j6aUttKFvW26sTEZEU4GZ3focCP7que7/jOAHAjc+ylcu6ka4H/94Vdl2I3J2VO+evoXLe1+h/61E254Ri6Qox4e5+PFrxUQJ8rf1TYhawidyw8+OGx461jg0ffwxt20K3btC9uxWziYiIxNMNB7+O42QGbgM6ALiuGwVEJcyy5Lwb6XqQLziIfcfDObnvMNmPz6Oe5xf+fOQQk4KhfLqiTGn6HveXvR9fn7i+pknVUUIk3mJiYMQIGDcONm2ygrVOnSDkn9qFjBm9uz4REUmRbmbntxhwBPjEcZxKwAqgh+u6ZxNkZXJByyr54x2IhodDhdNVOPj9l5xp1Z2TPi47HMh+Ljsvhwzk7SYdcBznovtcK0dYJMm4LuzeDUWK2NCJCRMgc2bb7X3oIUivk0siInJzbmYerR9QFRjtum4V4Czw0n9v5DhOJ8dxljuOs/zIkSM38XRyNcePwzvvQJWSmwldeBf7Wz5DrK8LDoBD46pP8k7Txy8JfOHqHSVEksSZM5bWULkyVKhgPzsO/PmnFbY9/rgCXxERSRA3E/zuA/a5rrvkn5+/xoLhi7iuO8513RDXdUNyappSgtuzB55/HkqW3smcBfUIbV+G2XesoBp5CfQNxNfxJcgvHU/XaXHFx7iRjhIiCWL3bnjmGciXD7p0AR8f+PDDuFHD6s0rIiIJ7IbTHlzXPeg4zl7HcUq7rrsFaAhsTLilydWsWwcDB8Lk2bsoUftJznWexxJfePBUAfq2HEXFGveyaO+ieHVvuJGOEiI3LCoKTp600cInTlh/3gcftAK2mjVtx1dERCSR3Gy3h+7ApH86PewAHr/5JcmVuK51dxowAL5f/DcFbn8WnvmRbUD7/Tl46YGhlGwU1+c0vt0bbqSjhMh127PHUhs++ggaNbLxw5Urw8GDEBzs7dWJiEgacVPBr+u6q4FLxsZJwoqNhQFfLmLUDwvYt74gOcp/ic8zP3A0FrpsyUTvZu9S6PVudsr4BtxIRwmReJs/34ZRzJ5tPzdrBu3axV2vwFdERJKQJrwlYxER8Pnn8NbERextcAeUioTScDYKXlgRQM/bXiT3py9DunQ3/VzX01FC5JrCwqxLg48PzJljRWsvvWStygoX9vbqREQkDVPwmwyFhcGYMTBkqMuhwN8JfKAj+EUC4HigV0x13pr4PeTI4d2FivzX2rUwciR88QXMmAF33QWvvgrvvgsBAd5enYiIiILf5CQ01M4Ojxnrcib3jwQ/8g5kXkj6CAdPLHh8IMA/HU2eHKrAV5KPmBgLdIcPt6T0dOlsAluhQna9OjaIiEgyouA3Gdi40To3fDHJQ2zJGQQ//Q6kW0mm8ADe+h46xlRg9WtPsSDb6Wt2bhBJMlFRtpvrutCjh30/YAA88QRkz+7t1YmIiFyWgl8vWrgQ+veHWXNi8K8ylcwvvstxv41kj8rIBzPh0UNZCXj7PWjfntq+vijkFa9zXVi61HZ5//oLtmwBf3/47TcoWjSuP6+IiEgypeA3iXk8VvT+yphFrDv7C+kynyXrq9M44ewgn5uTETP9eWBzNH69XoEXX4SMGb29ZBGIjISvvoJhw2D5csiUyaaunTtnaQ0lSnh7hSIiIvGi4DeJREbCpEmW3rA5Yj60awS+0UQARXzz8MmcYO5dcgSfhx+BGe/F5UuKJAe//Qbt28Mtt1hBW7t2FgCLiIikMAp+E9mpUzBuHAweDPuPniZv89EEVniTSDcaAB8PtPvlIC2oBX8Nhlq1vLxiEWDFChg61D6Evf22dW2YPx9uv10T2EREJEW7sakIck0HDkDfvhY79H71BIF3v0nGVwtzoPyLVMpWnECPD76xEOhxuOOJNyx/UoGveFNMjKU23HorhITAt9/GBbqOA/XrK/AVEZEUTzu/CWzrVvjgA/j0U4gOOEzpxwYTk3ckO2NO06JIU/63IRvVe01lUeEAFrS7jfoP96V2ifreXrYIPPecpTQUK2anKh5/XG3KREQk1VHwm0CWLLEuT99+CwHZQ7mlxwdszTyWLbERPFjqAV4+VpaKL4yGQ4fgsceo/e671M6Xz9vLlrRs/XpLbejeHSpWhK5d4e67oWlTdW0QEZFUS8HvTXBd+OEH+N+YRaw+sYD04aWp8urPrPf7mA1uLI+Wf5S+QXdT+uUPYeVXULs2zJoF1at7e+mSVrkuzJ1rO7s//WQDKerVs+C3XDm7iIiIpGIKfv9jxqpQBs7dwv6wcPIFB9G7UWlaVsl/0W2io2HKFNvpXR+2CDo0AL8IzgFrffx4ssqT9CnalqJvj4QpbaFAAWv18PDDypkU73Fd+wC2ZAnkzQvvvAOdO2sghYiIpCkKfv9lxqpQ+k5fR3h0LAChYeH0nb4OgJZV8nPmDHz0EQwaBHv3Qok66yj4cBf2RkcA4ODwQsizvLc4Izx4twUbr70GffpAhgxee12Shh08CNOmwTPP2AevNm3g6afhoYdsIpuIiEga47ium2RPFhIS4i5fvjzJni8+u7j/Vvf9XwkNC7/keE6/zDR26jFyJJw4AVWbriDwrrdZFDaDAJ90RHuicV0XP9eHn2Zk4Y61xyy46N8fChdOzJcocnlr11pqw5df2qmK1asttUFERCSNcBxnheu6If89nmp3fq+1i3s5+/8T+EafSM+ppcXYs74AK2LhtrZ/EVXrbRYd/YHgiGDalO7Jig21yHl4IZX3fcmja4+QOSYrf3w0jnodWyXuCxS5nH37rEvDL79A+vTw1FPQoweULOntlYmIiCQLqTb4HTh3y4XA97zw6FgGzt1yxeA3X3AQoWHhRB7MzKnFxTm3NS/4xJL1rm8o1XwMvx2eT45zOXiv4Xt0q96Nh9/6mTdmjeOhtT9zLH0WBtz+LN+Ub0jeIxlZeJW1Xe+OtMhVRUTA9u1WrJYzp52eeO896NQJsmXz9upERESSlVQb/P53F/dax10X7sxUkQ/GOpyL2QolPiag8jl8Sn7HCb8N7DmXl0F3D6JTtU5kcAJg1Cgmf/A/gqIj+Kh6S4bXbcPpwAxXfQ64sR1pkcs6dgzGjIHhw61rw7ZtEBgIy5apsFJEROQKUm3we34X93LH/+38UKsBA2DNmhxkKv8HtGoITjRRDgQH5GDQnaN4vMrjpPNLB/PmwbPPwsaNbCxZnf/d/gTbsxe86nP8243sSItcZPdum6Ty8cdw7hw0bgy9e8f15lXgKyIickWpdrxx70alCfK/uFF/kL8vvRuVBuDsWdswK1kS2raFyOhYuoyYStBjD4BPNDjggw+96vaga/WupNt3EFq3hjvvtNPMM2dycMp09ucpcsXnuJzr3ZEWuSD2nw9Na9fC2LHw4IOwbp01m27QQEGviIhIPKTand/zu6j/za29tWB+Xn8dRoyws8Z1bo2hZb8v+eHMu4w5uoXCWQoT4BtArCeWAN8AGuatC/362dawj4/1Ru3ZE9KloyWA41xX/m58d6RFAPB44PvvYeBAuPVW+/1r2hR27QJNCBQREbluqTb4BQuAzweiO3daf95HJkB4ODRrEUW5Rz5j2sH3GLJ7B5VyV2LaA9NoVaYVS/YtYcGu+dTf4VK74WPW1Pfhhy0ALlDgis8RH70blb4o5xeuvVssaVBkJHzxBXz4IWzaBAULwqOP2nU+Pgp8RUREblCqDn4BNm6Et9+2vF4fH3i4XQSFW37Mp9v7M3vTHkLyhTC40WDuLXUvzj+njWufzETtfr/A/PlQqZJNZ6tXL0HWc6UdaeX7ykW6dbOc3sqV7ffvgQfA39/bqxIREUnxUn3wO3PFIr45vIDWL9SiTP3VjNswkAMrD1CnYB3GNhtLo+KNLgS9hIVxISciSxYYPdr6pPr6Xu0prtv17hZLGnDgAAwdaj16S5eG55+3sw0NGyqXV0REJAGl6uB30d5FvLWnIdH1IvgKF5bAHUXuYFKrSdQvUj8u6PV44PPPbQzxkSPQpQu89RZkz+7dFyCp399/Wz7vp59a65GiRS34LV/eLiIiIpKgUnXwu2DXAiJiInCxEc6dq3VmTLMxF99o1Sp45hn46y+oVcsq56tW9cJqJU1xXejQwT50BQTAE09Ar15QooS3VyYiIpKqpdpWZwD1i9QnnV86fB1fgvyCeKzSY3FXHj8OTz8NISG2+/bJJ7BwoQJfSTyuC4sW2VfHgfz54aWXrHPD6NEKfEVERJKA47pukj1ZSEiIu3z58iR7PrDUhwW7FlC/SH1qF6xtKQ4ffwx9+1oA/Mwz8MYbEBycpOuSNCQ2FqZPh/ffh5Ur4fffE6yAUkRERC7PcZwVruuG/Pd4qk57AKhdsLYFvWBjX59+2r7Wq2eFbRUreneBknpFR1u7svffh61boVQpGD8eatTw9spERETSrFQf/AJw9Kjt9E6YALlzW0DyyCOqopfEFRVlRZQFC8K0aXDffQneOURERESuT+oPfmfNgsceg9OnbTLba69B5szeXpWkRqdOWe7u7NmwYAFkyABLl0KRIvqgJSIikkyk6oI3wFpHVa8Oa9bABx8o8JWEd/QovPoqFCpkBWwZMtjsbLDfPwW+IiIiyUbq3/ktXx7mzvX2KiS1WrMG6tSBc+egVSt4+WWoVs3bqxIREZErSP07vyIJbccO+P57+758eejaFTZsgG++UeArIiKSzKX+nV+RhLJ1K7z7rhVM5slj/Xn9/CydRkRERFIE7fyKXMu2bfDoo1CmDEydCt27WyGbnz47ioiIpDR69xa5Eo8HfHzg4EH49lvrFtKrl+36ioiISIqk4Ffkv1atgrfftiB35Ei49VYIDdUUQBERkVRAaQ8i5y1bBs2bQ9WqMG8e5MsXd50CXxERkVRBO78iAIMHW1pD1qzw1lvwzDMKeEVERFIhBb+Sdi1ZApkyQdmycO+9EBkJ3bppEIqIiEgqprQHSXuWLoUmTaBWLXjnHTtWooRNZ1PgKyIikqop+JW0Y8UKaNYMata0APj992HsWG+vSkRERJKQ0h4k7ZgxA/76y3Z7u3e3lAcRERFJU7TzK6nX6tXQsiV895393KePTWV7+WUFviIiImmUgl9JfdauhVatoEoVWLAAjh6145kyKadXREQkjVPwK6nL889DpUrWp/f1122n94knvL0qERERSSaU8ysp386dNpAiMNAGVPzvfzaGOGtWb69MREREkpmb3vl1HMfXcZxVjuPMTogFicTbvn3QpQuUKgXjxtmxdu1sNLECXxEREbmMhNj57QFsApRMKUnj8GF47z0YPRo8HujUCVq39vaqREREJAW4qZ1fx3EKAE2BjxJmOSLx0KYNDBsGbdvC1q0wcqSlPYiIiIhcw83u/A4B+gBX7BvlOE4noBNAoUKFbvLpJE06fRqGD7cd3hw54MMPIUMGS3cQERERuQ43vPPrOE4z4LDruiuudjvXdce5rhvium5Izpw5b/TpJC2KiIBBg6BYMStimzXLjleposBXREREbsjNpD3UBZo7jrMLmAI0cBzniwRZlcjEiRbg9uplwe6SJfD4495elYiIiKRwNxz8uq7b13XdAq7rFgHaAL+6rvtogq1M0rbZsyFvXuvX+9NPUKOGt1ckIiIiqYCGXEjyMH8+1K0LmzbZz598AosXQ4MG3l2XiIiIpCoJEvy6rrvAdd1mCfFYksasXAmNGlmQu3cvHDhgxzNlAsfx7tpEREQk1dHOr3iH60KHDlCtGixfbh0ctm7VTq+IiIgkKo03lqR19Ki1K3McKFQIXnkFXngBsmTx9spEREQkDdDOrySNU6cs0C1UCH75xY69+Sa89ZYCXxEREUky2vmVxBUVBePGwRtv2K7vww9D8eLeXpWIiIikUQp+JfG4LtxxB/z1F9SvDwMHQkiIt1clIiIiaZjSHiThLV4MMTGW19ujB8yZA7/+qsBXREREvE7BrySczZuhZUuoXRs+/9yOPfggNGmitmUiIiKSLCj4lZt38CB07Qrly9sO79tvW9ArIiIikswo51duXsuWsGKFBcCvvgq5cnl7RSIiIiKXpeBXrl9sLHzxBdx3H2TODMOHQ3AwlCzp7ZWJiIiIXJXSHuT6nC9c69ABPv3UjlWvrsBXREREUgQFvxI/W7ZAixbQsCEcPw6TJ8Mzz3h7VSIiIiLXRWkPEj+9esHvv8N771n7sqAgb69IRERE5Lop+JXLi4qCkSMtr7dIERgxwgLe3Lm9vTIRERGRG6a0B7mY68K330K5ctCzJ0yZYseLFFHgKyIiIimegl+Js3KljSNu1QoCAuCHH+Cll7y9KhEREZEEo7QHiTN+PGzcCKNHw5NPgp9+PURERCR10c5vWhYVBR9+CEuW2M/vvgt//w1duijwFRERkVRJwW9aNWeOjSN+4QWYPt2OZc0KWbJ4d10iIiIiiUjBb1qzeTM0aQLNmoGPD3z/PfTv7+1ViYiIiCQJndtOa777DhYuhEGD4OmnrbBNREREJI1Q8JvaxcbCxx9DzpzQsqUNqOjQAXLl8vbKRERERJKc0h5Ssz/+gOrVoVMnmDrVjgUGKvAVERGRNEvBb2q0fz+0bQu33QZHj9qgii+/9PaqRERERLxOwW9q9Oef8M038MorVuD20EPgON5elYiIiIjX/b+9e4+Vor4COP49Iib4qlrQomItagr4oBpUfBEjWoWotA0SjfFV4ysaJbahxndMfLXWV2M0viKtxlcURQVjA8SmRowUUECsokG0IkINYq2xor/+MUO8LruXu9zrztyd7yeZ7NyZ37InJ2eGk9nfzDrnt13MmAHLlsEZZ8AJJ8BBB8GgQUVHJUmSVCpe+e3t3n8fJkyAI4+EW2/NbnCLsPGVJEmqw+a3t/ryS7j+ehgyBJ55Bq65Jvultj59io5MkiSptJz20FstWACXXpo9vuyWW2DXXYuOSJIkqfS88tubvPce3HNPtj5iRNYAT5li4ytJktRFNr+9wZdfwrXXwtChcPHF2ePLAPbaq9i4JEmSehmb37KbNQuGD88eWzZ2LCxaBP37Fx2VJElSr+Sc3zJbtQqOPRZ22AGmTYMxY4qOSJIkqVfzym/ZfP01TJ0KKWVXeKdPz6722vhKkiR1m81vmcybBwcfDOPGwcyZ2bZRo6Bfv2LjkiRJahM2v2WwZg1cdFH2BIelS+HBB+GII4qOSpIkqe0457doKWWN7ty5cO652VMdtt226KgkSZLaks1vUd59N/sJ4r59v214Dzig6KgkSZLamtMeWu2rr+C662DYMLj99mzb0Ufb+EqSJLWAV35bafZsOOssWLgQxo+Hk04qOiJJkqRK8cpvq9x0U/Ykh9Wr4emn4fHHYccdi45KkiSpUmx+v29r12avhx4KF1wAb7wBxx9fbEySJEkV5bSH78uHH8KFF2a/znbHHTByZLZIkiSpMF757WnffAN33QVDh8Jzz8EuuxQdkSRJknJe+e1JS5bA6afDSy/B6NFZE7z77kVHJUmSpJzNb0/aZBNYtgweeABOPRUiio5IkiRJHdj8dtfLL8Mjj8Ctt8LgwfDOO9kPV0iSJKl0NnrOb0QMiohZEbE4IhZFxEU9GVjpff45TJwIhxwCTz0FH32UbbfxlSRJKq3u3PC2FvhNSmkoMBI4PyKG9UxYJTdjBuy9N9x2G5x3XvajFQMHFh2VJEmSNmCjpz2klJYDy/P1zyJiMbAT8EYPxVZOX3wBp5wCW2wBL74Io0YVHZEkSZK6qEfm/EbErsC+wCs98e+V0syZcNhh0K8fTJ8Oe+wBm29edFSSJElqQref8xsRWwJPABNTSmvq7D87IuZExJyVK1d29+Nab9UqOPnk7NFl996bbRs+3MZXkiSpF+pW8xsRfcka34dSSk/WG5NSujulNCKlNGLAgAHd+bjWSgkeewyGDcter7oKzjyz6KgkSZLUDRs97SEiArgPWJxSurnnQiqJSZPgpptgxIhvb3CTJElSr9adOb+HAKcACyJifr7t0pTStG5HVZSUYO3a7HFl48fDgAFw8cWwqY9DliRJagfdedrD34H2+QmzFSvgnHOyR5bdeScceGC2SJIkqW10+4a3tvDoo7DnnvD889lTHFIqOiJJkiR9D6rd/K5cCRMmwIknwm67wfz52TSHaJ8L2pIkSfpWtZvfzz7Lbma7/np46SUYMqToiCRJkvQ9qt6dXJ98ApMnw8SJMHgwLF0KW21VdFSSJElqgWpd+X322Wxu76RJsGBBts3GV5IkqTKq0fyuXg2nnw7HHQfbbw+vvgr77FN0VJIkSWqx9p/2kBIcdRTMmweXXw5XXAGbbVZ0VJIkSSpA+ze/EXDDDbD11rD//kVHI0mSpAK1f/MLMHp00RFIkiSpBKox51eSJEnC5leSJEkVYvMrSZKkyrD5lSRJUmXY/EqSJKkybH4lSZJUGTa/kiRJqgybX0mSJFWGza8kSZIqw+ZXkiRJlWHzK0mSpMqw+ZUkSVJl2PxKkiSpMmx+JUmSVBk2v5IkSaoMm19JkiRVhs2vJEmSKsPmV5IkSZURKaXWfVjESuC9ln3gt/oDqwr43N7KfDXPnDXHfDXHfDXHfDXHfDXHfDWnyHz9OKU0oHZjS5vfokTEnJTSiKLj6C3MV/PMWXPMV3PMV3PMV3PMV3PMV3PKmC+nPUiSJKkybH4lSZJUGVVpfu8uOoBexnw1z5w1x3w1x3w1x3w1x3w1x3w1p3T5qsScX0mSJAmqc+VXkiRJaq/mNyKOiYh/RsSSiLikzv6IiNvz/a9HxH5FxFkGETEoImZFxOKIWBQRF9UZc3hEfBoR8/PlyiJiLYuIWBoRC/JczKmz3/rKRcRPO9TN/IhYExETa8ZUvr4i4v6I+DgiFnbYtl1E/DUi3s5ft23w3k7Pd+2oQb7+EBFv5sfclIjYpsF7Oz1+21GDfF0dEf/qcNyNbfBe6yvb9miHXC2NiPkN3lvF+qrbR/SKc1hKqS0WoA/wDjAY2Ax4DRhWM2YsMB0IYCTwStFxF5ivgcB++fpWwFt18nU48GzRsZZlAZYC/TvZb33Vz0sf4COy5y123F75+gJGAfsBCzts+z1wSb5+CXBjg5x2er5rx6VBvn4ObJqv31gvX/m+To/fdlwa5Otq4LcbeJ/1VX//H4ErG+yrYn3V7SN6wzmsna78HgAsSSm9m1L6H/AIMK5mzDjgzykzG9gmIga2OtAySCktTynNzdc/AxYDOxUbVa9nfdU3GngnpVTED9yUWkrpb8AnNZvHAZPz9cnAL+q8tSvnu7ZTL18ppRdSSmvzP2cDO7c8sJJqUF9dYX3ViIgAJgAPtzSoEuukjyj9Oaydmt+dgPc7/P0B6zdzXRlTORGxK7Av8Eqd3QdFxGsRMT0i9mxtZKWTgBci4h8RcXad/dZXfSfS+D8M62t9O6SUlkP2nwuwfZ0x1lp9vyb79qWeDR2/VXJBPk3k/gZfSVtf6zsMWJFServB/krXV00fUfpzWDs1v1FnW+2jLLoyplIiYkvgCWBiSmlNze65ZF9VDwf+BDzV4vDK5pCU0n7AGOD8iBhVs9/6qhERmwHHA4/X2W19bTxrrUZEXAasBR5qMGRDx29V3AnsBvwMWE72VX4t62t9J9H5Vd/K1tcG+oiGb6uzrWU11k7N7wfAoA5/7wx8uBFjKiMi+pIV7EMppSdr96eU1qSU/pOvTwP6RkT/FodZGimlD/PXj4EpZF/bdGR9rW8MMDeltKJ2h/XV0Ip102Xy14/rjLHWOoiI04BjgZNTPqGwVheO30pIKa1IKX2dUvoGuIf6ebC+OoiITYFfAY82GlPV+mrQR5T+HNZOze+rwB4R8ZP8atOJwNSaMVOBU/O78kcCn667NF81+fyl+4DFKaWbG4z5UT6OiDiArF7+3booyyMitoiIrdatk91ks7BmmPW1voZXS6yvhqYCp+XrpwFP1xnTlfNdJUTEMcDvgONTSv9tMKYrx28l1NyH8Evq58H6+q4jgTdTSh/U21nV+uqkjyj/OaxVd9a1YiG72/4tsjsIL8u3nQucm68HcEe+fwEwouiYC8zVoWRfMbwOzM+XsTX5ugBYRHYX5mzg4KLjLjBfg/M8vJbnxPracM42J2tmf9Bhm/X13Rw9TPbV81dkV0LOBH4IzADezl+3y8fuCEzr8N71znftvjTI1xKyuYPrzmN31ear0fHb7kuDfP0lPz+9TtZsDLS+Gucr3/7AuvNWh7HWV+M+ovTnMH/hTZIkSZXRTtMeJEmSpE7Z/EqSJKkybH4lSZJUGTa/kiRJqgybX0mSJFWGza8kSZIqw+ZXkiRJlWHzK0mSpMr4P0hCuGT7Vl/7AAAAAElFTkSuQmCC\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.58851503 0.38554724]\n",
      "[0.4180715  0.03602273]\n"
     ]
    }
   ],
   "source": [
    "X2 = X[:,[0,1]]\n",
    "res2 = sm.OLS(y2, X2).fit()\n",
    "print(res2.params)\n",
    "print(res2.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate RLM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.05362915 0.4841294 ]\n",
      "[0.1107057  0.00953885]\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+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABstUlEQVR4nO3dd3SUVRPA4d+bTaX3FnoH6QQUEOlSpKMICtIEBCsqghVQEBQ/sYGELtKbKDYUEFGIdJAaegu9hADp2fv9MYQAJiSk7W4yzzkcyO6b5C6b7Oy9d+6MZYxBKaWUUunLzdEDUEoppTIjDcBKKaWUA2gAVkoppRxAA7BSSinlABqAlVJKKQfQAKyUUko5gHt6frN8+fKZkiVLpue3VEoppRxm69atF40x+eO7L10DcMmSJdmyZUt6fkullFLKYSzLOp7QfboErZRSSjmABmCllFLKATQAK6WUUg6QrnvA8YmKiuLUqVOEh4c7eihpxtvbm6JFi+Lh4eHooSillHISDg/Ap06dInv27JQsWRLLshw9nFRnjOHSpUucOnWKUqVKOXo4SimlnITDl6DDw8PJmzdvhgy+AJZlkTdv3gw9w1dKKXX/HB6AgQwbfGNl9MenlFLq/jlFAHYmI0eO5JNPPknw/uXLl7N37950HJFSSqmMyOUC8PLtQTQYt4ZSw3+iwbg1LN8elL7fXwOwUkqpVOBSAXj59iDeXLaLoOAwDBAUHMaby3alOAiPGTOGChUq0Lx5cwIDAwGYOnUqderUoXr16nTp0oXQ0FA2bNjADz/8wNChQ6lRowaHDx+O9zqllFIqMS4VgMevDCQsKuaO28KiYhi/MjDZX3Pr1q0sWLCA7du3s2zZMjZv3gxA586d2bx5Mzt37qRSpUpMnz6d+vXr0759e8aPH8+OHTsoU6ZMvNcppZRSiXH4MaT7cTo47L5uT4q//vqLTp06kSVLFgDat28PwO7du3nnnXcIDg7m+vXrtGzZMt7PT+p1SimlnNzVq5AzZ7p9O5eaARfJ5XNftydVfFnKvXv35quvvmLXrl2MGDEiwWNESb1OKaWUE7LbIXbrcMsWOHMm3b61SwXgoS0r4ONhu+M2Hw8bQ1tWSPbXfOSRR/juu+8ICwvj2rVrrFixAoBr165RuHBhoqKimDt37q3rs2fPzrVr1259nNB1SimlnFhwMHz+OVSsCKNGyW1NmoCnZ7oNwaWWoDvW9AVkL/h0cBhFcvkwtGWFW7cnR61atXjyySepUaMGJUqUoGHDhgB88MEHPPjgg5QoUYKqVaveCrrdunWjf//+fPHFFyxZsiTB65RSSjmhXbtg4kT49luZ+darBw89JPe5uUHevOk2FMsYk27fzM/Pz9zdD3jfvn1UqlQp3cbgKJnlcSqllNOJiQHbzdXTp5+GZcvgqafg+eehVq00/daWZW01xvjFd59LLUErpZRSSXbmjCwvFy8uM1+AcePg1CmYPj3Ng29iXGoJWimllLonY2D9evjqK1i6FKKjoVUrmQUDFCvm2PHdJtEZsGVZMyzLOm9Z1u67bn/RsqxAy7L2WJb1cdoNUSmllEpE7Hbq9esScFeuhBdfhAMH4JdfoEYNhw4vPkmZAc8CvgJmx95gWVYToANQzRgTYVlWgbQZnlJKKXUPBw/CpEmwYwesWQPZs8Ovv0LNmpA1q6NHd0+JzoCNMeuAy3fdPAgYZ4yJuHnN+TQYm1JKKfVfMTHwww/QsiWULy/LzYUKxZ3nffhhpw++kPwkrPJAQ8uyNlqW9adlWXVSc1BKKaVUghYtgg4dYPduSbI6cQLmz3eJoHu75CZhuQO5gYeAOsAiy7JKm3jONFmWNQAYAFC8ePHkjjPNXLp0iWbNmgFw9uxZbDYb+fPnB2DTpk14puOhbKWUUvHYtEnO7tapAy+8AJ06wZIl0L49eHg4enTJltwAfApYdjPgbrIsyw7kAy7cfaExZgowBeQccHIHmlby5s3Ljh07AOkFnC1bNl5//fVb90dHR+PursniSimVrsLDYeFCWV7eskVmt+XKyX3e3tCli2PHlwqSG1mWA02BtZZllQc8gYupNShH6927N3ny5GH79u3UqlWL7Nmz3xGYq1Spwo8//kjJkiWZM2cOX3zxBZGRkTz44INMmjQJm82WyHdQSil1T7EFMypVgi+/hGeegRw5HD2qVJVoALYsaz7QGMhnWdYpYAQwA5hx82hSJNArvuXn+/XKK5LIlppq1IDPPrv/zztw4ACrVq3CZrMxcuTIeK/Zt28fCxcuZP369Xh4eDB48GDmzp3LM888k5IhK6VU5mK3w2+/wddfS0azry8MHSqVqpo0gXga5mQEiQZgY0z3BO7qkcpjcSpPPPFEojPZ1atXs3XrVurUkRy0sLAwChTQE1lKKZUkV67AzJkSeA8dggIFYP9+CcCx9ZkzMKfa3EzOTDWtZL0tm87d3R273X7r49iWg8YYevXqxdixY9N9fEop5dJCQqBECbh2DRo0kGzmLl3Ay8vRI0s3Wgs6CUqWLMm2bdsA2LZtG0ePHgWgWbNmLFmyhPPn5Rj05cuXOX78uMPGqZRSTisiAubOlaVlkP3cjz6C7dvh77+lOUImCr6gAThJunTpwuXLl6lRowZff/015cuXB6By5cqMHj2aRx99lGrVqtGiRQvOpGMzZ6WUcnonTsDbb0tDhB49YMUKuHFD7hs0yClLRKYXbUeYTjLL41RKqVuWLYMnnpB/t20rSVXNm0vf3UziXu0InWoPWCmllAu7cgW++QbKlIF27aBRIxg+HAYMkP1edYfM8zZEKaVU2tixA/r3l+zlIUPgp5/k9rx5YcwYDb4J0BmwUkqp5Hv2WWlu7+MjiVSDBzu80b2r0ACslFIq6U6cgKlTJZs5Rw7pSPTAA9C7N+TO7ejRuRQNwEoppe7NbodVq6RK1YoVctuDD0piVWySlbpvGoCVUkolLDgY6taVxvf580tS1cCBcqxIpYgmYQGnTp2iQ4cOlCtXjjJlyvDyyy8TGRnJ2rVradu27X+u//HHH6lZsybVq1encuXK+Pv7O2DUSimVRrZtk31dgFy54NFHYc4cOHlSkqo0+KaKTB+AjTF07tyZjh07cvDgQQ4cOMD169d5++23470+KiqKAQMGsGLFCnbu3Mn27dtp3Lhx+g5aKaVSW3g4fPut1GCuXVv2eMPC5L6vvpLuRJmsUlVay/QBeM2aNXh7e9OnTx8AbDYbEyZMYMaMGYSGhv7n+mvXrhEdHU3evHkB8PLyokKFCuk6ZqWUSlU//wzFiknLv+Bg+PxzOHJEMptVmnGuPWAH9CPcs2cPtWvXvuO2HDlyULx4cQ4dOvSf6/PkyUP79u0pUaIEzZo1o23btnTv3h23TFTZRSnl4ux2WLlS9nT9/KBCBWjYUCpVNW2aYdv/OZtMHzWMMVjx/LAldDvAtGnTWL16NXXr1uWTTz6hb9++aT1MpZRKuYsXYfx4KFsW2rSBL76Q28uUkbKRzZpl+uAbW6Y6PTjXDNgB/QgfeOABli5desdtISEhnDx5kjJlyiT4eVWrVqVq1ar07NmTUqVKMWvWrDQeqVJKpcCwYbK0HBEBjzwC48ZBx46OHpXT2LRJekYA/P57+nzPTD8DbtasGaGhocyePRuAmJgYXnvtNXr37k2WLFn+c/3169dZu3btrY937NhBCS2zppRyNqGhUpc5Kko+LlAA+vWDXbvgzz+ha1fw9HTsGFPJ8u1BNBi3hlLDf6LBuDUs3x6U5M/dvRs6dZJjzTt2QOvWskKfHpxrBuwAlmXx3XffMXjwYD744APsdjtt2rThww8/JCAggNWrV1O0aNFb18+fP5+PP/6YgQMH4uPjQ9asWXX2q5RyHgcOwNdfw6xZklCVLx889hi89pqjR5Ymlm8P4s1luwiLigEgKDiMN5ftAqBjTd8EP+/wYRg5UloUZ88O778vaUjZs6fDoG/K9AEYoFixYqyIre5ym8aNGxMWm4Z/m4YNG6bHsJRSKukuXYJu3aRilbs7PP641GV++GFHjyxNjV8ZeCv4xgqLimH8ysB4A3BQEIweDdOmgYeHnLZ64w3pG5HeNAArpZSrOntWlpRbtJA6zHa7RJd+/aBQIUePLl2cDv7vJCm+2y9ehI8+kiPNMTHSIfHtt6FIkbhrAk4GsPbYWhqXbEy9YvXSctiABmCllHItxsC6dVKXedkyaYhw5ozs565e7ejRpYnl24MYvzKQ08FhFMnlw9CWFW7Nbovk8iEoniBcJJecYQ4JgU8/lT83bkCPHrL0XKrUndcHnAyg6eymREZH4uXuxepnVqd5EM70SVhKKeUyVq2CKlWgcWNJ1X3pJfjnnwyTTBWf2D3eoOAwDHF7vLGJVkNbVsDHw3bH5/h42HipUQU++QRKl4ZRo6Sa5q5dkpd2d/A9dPkQz//8POHR4dixExkTydpja9P8sekMWCmlnNnOnZIZVLq0/J01K8ycCU8+mSkqVSW2xxs7E46dIRfKloXqoTV59fFcnD4t3RJHj5Z6I3c7cuUIo9eNZvbO2djcbLi7uWOMwdPmSeOSjdP8sWkAVkopZxMRAUuWyDLzhg3SfWjyZDkrs2mTo0eXrpKyx9uxpi/tqvkybx6MGAH+R6FBA5g3Dxo1+u/nHg8+zuh1o5m1cxY2y8YLdV9gWINhHAs+pnvASimVaY0ZI0WJLl6EcuVk87JXL0ePymES2+M1BpYvh3fegb17pfrwzz9Dq1b/Leq1fN9yxv49lq1ntmJzs/Fc7ecY/vBwfHPILLpw9sLpEnhj6R4w0oChRo0aVKlShXbt2hEcHAzAsWPHqFKlyn+ujy3Sce3atVu3vfzyy1iWxcWLF9Nr2EqpjCAmRpKnjJGPL1yQusy//w7798OQIZAnj2PH6EAJ7fG+/mgFfv9dWhV37iz/jQsXwtatUkzj9uAbFBLE44sep9OiTmw6vQkLi4VdFvJlmy9vBV9H0AAM+Pj4sGPHDnbv3k2ePHmYOHFiop9TtmxZvv/+ewDsdjt//PEHvr6OeyKVUi7m/HkpB1m2LDRvDn//LbdPmCDZzc2bgzZ5oWNNX8Z2ropvLh8swDeXD73K1OaL13x59FH5b5wxQypade1653/ZmWtnePmXlynzRRm+2/8dFhKVDYZ9F/fd8X1SUk0ruVzy2Q04GcDYv8YScDIg1b92vXr1CApK/D++e/fuLFy4EIC1a9fSoEED3N11RV8plYjz56W3brFi8OabULIkLF4sfXgh0zdDiE/Hmr6sH96U77o9Ru6/m/Jmn/zs2wdffimFv/r0kdojsc5dP8erK1+l9Belmbh5Ik9XfZpFjy/C290bm2X7T5JVYpnWacWpIsYrv77CjrM77nnN1Yir/HvuX+zGjpvlRrWC1cjplTPB62sUqsFnrT5L0vePiYlh9erV9OvXL9Fry5Urx/fff8+VK1eYP38+PXr04JdffknS91FKZTLXrkntwxo1IGdOSaR67jn5U6mSo0fn9A4cgPfekyXmXLlg7Fh48UVJCI8VcDKAnw7+xNHgo3y37zsiYiLoWa0n7z7yLmXySGOdItmLxJtkdb/VtFKLUwXgpLgafhW7kUrZdmPnavjVewbgpAgLC6NGjRocO3aM2rVr06JFiyR9XufOnVmwYAEbN27E398/RWNQSmVAu3dLXeZvv5WazIcOgZcXBAbq8nISnDghNZpnzQJvb3jrLSkdmSvXndf9evBX2i1oR7Q9GoCWZVryResvKJ+3/B3X1StWL94kq6RW00ptThWAkzJTDTgZQLPZzYiMicTT5sncznNTnLUWuwd89epV2rZty8SJE3nppZcS/bxu3bpRq1YtevXqhZv+MimlYq1dC+++K/u6Xl5yZnfw4LjlZX29uKdz5+DDD+XkFchs9803paHT7S6HXebTgE8Zv2H8reBrs2w0KtHoP8H3XhLLtE4rThWAk6JesXqsfmZ1mpzVypkzJ1988QUdOnRg0KBBiV5fvHhxxowZQ/PmzVNtDEopF3X0qKyJFigg9Q/PnIHx42WD0hGV/p3IvUpJ3u7KFfjkEzmFFREh/3XvvgvFi995XXB4MBMCJvDZxs8IiQihaammbDi5gaiYqGQV0RjassIdHZVAMq2HtqyQjEebdC4XgCHhZYTUULNmTapXr86CBQto2LAhgYGBd7QjnDBhwh3XDxw4ME3GoZRyATEx8Msvssz8yy+yRjp6NLRtK390ppukdoHXr8MXX8j7leBgaeo0ahSUv20SG3AygF8P/crpa6dZvHcxVyOu0qVSF0Y0GkHVglVT1Ejh7mpa93qTkJosE3v2LB34+fmZLVu23HHbvn37qJQJkhAyy+NUKtP45BNprXP8uHQe6t9f/hQr5uiROZUG49bEu7zrm8uHNUOa4u8vtUfOn5e2xWPGQPXqd1676sgq2sxtQ5Q9CoCGxRvyResvqFGoRjo8gpSxLGurMSaeQphJmAFbljUDaAucN8ZUueu+14HxQH5jjFagUEplXMZIXeYaNeTjzZuhTBmZtnXsKM1lM6CkLh8nJL5EJmO3CFyXl/JfS6JV48ZSzareXRPX65HXmbhpIu//+f6t4GuzbLQu29olgm9ikrIEPQv4Cph9+42WZRUDWgAnUn9YSinlJK5elSzmyZNhzx7Ytw8qVoQ5czJs0I2VlOXjxNye4GQMhO4vTPDf5Ym+nI3SdWD6dGjW7M7jz6FRoUzaPImP1n/ExdCLPOT7ENvPbifaHp1ujRLSQ6IB2BizzrKskvHcNQF4A/g+tQellFIOd+aMNI6dO1cayfr5SbSIzQjK4MEXUud87NCWFRi+dBdXAvNwZV0Fos7nxDP/NYb/7xIfDsl7K/AGnAzg9yO/cyn0Egv2LOD8jfM8WuZRRjYaSb1i9VK0x+uskpWEZVlWeyDIGLPTSqRqi2VZA4ABIFnD8THGkNjXcWXpuc+ulEqBsDA4fVqWlr29YelSOUI0aFD8/ewyuNQ4H5s7xBfrx3yc3+mFe64blOu6h7FDc9HFLy6Arz22lke/ffTWMrNfYT+WdV1Gg+INbl2Tlsm3jnLfKXqWZWUB3gbeS8r1xpgpxhg/Y4xf/vz5/3O/t7c3ly5dyrBByhjDpUuX8Pb2dvRQlFIJOXgQXnsNihaF7t3ltty5IShIZr2ZMPhCwudgk3I+dssW6cXbuDGEXPDi668h9HxWDix84FbwjYiOYNLmSXSY3+FW8HWz3OhcqfMdwTejSs4MuAxQCoid/RYFtlmWVdcYc/Z+v1jRokU5deoUFy5cSMZQXIO3t/cdR5mUUk7ijz+k4sOqVVJMuHNnKQ9pjGxKenk5eoQOlZzzsXv3ytndZcvk+PMnn0gNEp/bYnZkTCQzt89kzF9jOBlykmoFqxF4MTDD7fEm5r4DsDFmF3CrHollWccAv+RmQXt4eFCqVKnkfKpSSt2/U6ekHnP27FKfOTBQzu726yfHidQt93M+9uhR2TKfM0fqkYwaBa+8AjlyyP0BJwNYfXQ1NyJvMH/3fI5fPU69ovWY3n46zUs3559T/2S4Pd7EJHoO2LKs+UBjIB9wDhhhjJl+2/3HSGIAju8csFJKpTm7Xfrrfv01rFghpZZefBEiI8Fmkz8qWU6flrO7U6fKf+OLL8KwYXcW//rrxF80+6bZrWXmSvkq8WnLT2lZpmWGzv+BFJ4DNsZ0T+T+kskcl1JKpS1jZA108mQ4cgTy54c33oB27eR+T0/Hjs+FXboEH30kLQGjo+HZZ2XpuUiRuGui7dHM2zWPISuH3LHH26NaD1qVbeWgkTsPlyxFqZRSCTJG+tdVqCD7uD/9JMlVY8ZAp06Zfl83pa5dgwkT4H//k3/36CFLz6VLx10TY49h4Z6FjPpzFAcuHaBcnnJcj7xOjD0GT5snTUo2cdj4nYkGYKVUxhASElcwIzAQTp6EggWlRrNP2na1yQzCwmDSJBg3Di5elHy199+HBx6Iu8Zu7Czes5hRf45i38V9VC1QlWVdl9GhYgc2ntqY6fZ4E6MBWCnl2o4fl9ntvHlSMKNWLYkU2bPL/Rp8UyQqCmbMkGB7+jQ8+qjkrNWpE3fN+hPrmbR5EgGnAjgafJTK+Suz+InFdK7UGTdLTrtmxHO8KaUBWCnlesLCZBpWrJgsOc+bJwUznnvuzsigki0mBhYsgPfek+3z+vWlKFjjxnHXGGP46O+PeGvNWxgMFhajGo/i7YZvY3PTxLbEaABWSrmOwEBZYv7mG4kIP/4IJUtKK50sWRw9ugzBGPjhB3jnHdi9WzoT/fQTtG4dV6/ZGMOKAysYuXYk289uv/W5bpYbHm4eGnyTSJtVKqWc3y+/QNOm0gRh4kQpsTR0aNz9GnxTxerV8NBD0twpMhIWLoRt26BNGwm+xhh+PvgzdafVpcOCDoREhPDuI+/i4+6DzbJlqiIaqUFnwEop53TsmJxp8fSE7dul0sPYsdCnjyRXqSRJSjvBf/6Bt9+GNWtkVX/6dHjmGSkOBrDhxAambpvKpqBN7L24l5K5SjKj/Qx6VOuBh82D1mVba4JVMiRaiCM1aSEOpdQ9xcTAzz/LMvMvv8D8+bK3Gx4ugdhNF+3ux93tBEFKSY7tXJWONX35919Zal6xAgoUkCA8cGDcSS1jDF9u+pIhK4dgN3YsLIY1GMb7Td7Hw5ZBu0HZ7an6c5aiQhxKKZXmIiLg44+lnNLJk1C4sESDBjcL8mszk2RJqJ3g+3NPsGi8LwsWSKnIMWPgpZcgW7a469YeW8t7f7zHXyf+unWbm+VGDq8cGS/42u1SD9zfXzbBly1Ll2+rAVgp5Rh2u9RiLldOZrcLFsge72efSaWqTNBvN63d3TYwOsSbq+vLcXxXUQJ94M034fXXpfFTrHXH1zFi7QjWHltLkexFeK3ea0zaPInImMiMt8d7/jzMnAlTpkiqd968UtIrthlHGtMArJRKXxcvShazvz+cOyeHS7Nmhc2bNZkqlRXJ5UNQcBgxNzy5GlCWazukJ3vhekFsX1bs1lZ6wMkAvtn5DVtPb2XLmS0UylaIz1t9zoDaA/B296ZLpS4ZZ4/XGOmC5e8P330nB50bNZLDzZ07p2ulNA3ASqn0sX+/rHUuXixLzg0awIgRcTNdDb6pbnD9irz8VhiXN5XARNvIWuUUhRod5pO+5W4F32lbpzHwp4HYjR2Al+q+xLjm4/DxiCtgkiGKaNz+xu/gQZn2P/88DBgAlSo5ZEgagJVSaefqVQgNlT3d69flgGn//pLpU6WKo0eXYd24IU0SPvqoCMHBkLfqObwf3EeJMvZbWdBbTm9hxNoR/Hzw51ufZ7NsFMpW6I7g69KMgb/+kqC7ZImcrapfX7pGPP64w6ukaQBWSqW+LVskk3n+fOjWTc61+PnB2bMOf9HLyCIiJI9t9GhZ3X/sMfl3jRoFAZnybjuzjfbzB7HiwAry+ORhkN8gZu2YlbH2eK9cgdmzJfDu2yf9nwcMcLo3fhqAlVKpZ948+PRT2LpVlpSfegoGDYq7X4NvmoiOhjlzpCvR8ePwyCOwdGlcEnnAyQDm757PznM7WXd8Hbm9czO6yWhefPBFcnjloGe1nq6/x2uMHGiePBkWLZKjaw8+KIWsn3zSKbc4NAArpVJm717ZQ7Ms2LhRpmFffSV96nLmdPToMjS7XQLte+/JFrufnyT0tmgRl8Q759859F7emxgjx5Gerfksnzz6CTm9454bl97jvXpV3n34+8OuXXKWqndvme3WqOHo0d2TBmCl1P0LC5NkKn9/2LBBSig1aSId2r280uUIR2ZmDPz6qxyV3r4dKleWo6sdO8b91++9sJdRf45i0Z5Ftz7PZtkonbv0HcHXJRkTt82xYIHkGdSqJT+P3bvHdcJyclpWRimVdFevwquvgq8v9OolmaX/+59U7AcpmKHBN0399ZcsMbdpA8HBstX577/QqZP81wdeDOTpZU9TZVIVfj74M89Ueybj1Gq+dk2CbO3aULeuFKt+6ik5wrZ1q+zzukjwBZ0BK6USExEhBTMqV5Y93CVLZI1z4ECZ9WrATRdbt0rZyF9/laTySZOgXz/Yei6AjzespUzuMvx48Efm7pqLt7s3bzR4g9frv06+LPl4zu85197j3b5dAu/cuZJNX726/Ac8/bSU8nJRGoCVUvE7fFg2FGfOlJnt0aNSsergwXQtVpCRJaVRwr59cmpm6VLIkwfGj4fBgyWnKOBkAE1nNyU8OhwAT5snrz70KkMbDKVA1gK3voZL7vHeuCEz3MmTZYbr4yPJVAMHSnJVBnjjpwFYKXWnDRtg1Cj47Tew2aQs5MCBcS94GnxTxd2NEoKCw3hz2S4AOtb05ehReRq+/VaC7YgRMGRIXF7bseBjvPTLS7eCr4XFa/Ve48NmHzrk8aSaXbtktvvttxASIisvX3wBPXtCrlyOHl2q0gCslJKzK97e0uYvOFgym0eNkjVOX99EP13dv4QaJYxZfJRV032ZMkWa8gwZAsOGQf78cs2JqycYs24MM3bMwMLC3c0dYwyeNk/alW/ngEeSCu5O6vPygieekDd+DRpkiNlufDQAK5VZxbb+8/eXv994A8aNg1atpBevzeboEWZodzdKiAnzIGRjGU5sLckO5L3PO+9A0aJy/6mQU4z9ayxTt00FYECtAbzZ8E1OXj3puvu7+/bJz9/s2VI8o3x5Serr1UsaI2RwGoCVyow++ggmTryz9d+zz8p92nM3XcQ2SrBH2AjZUoqQTaUxke7kr3GOgCWFKFNGrvsh8AfG/T2Ozac3A9CvZj/eavgWxXNKY4WiOYq6VuCNiJAN7cmTJaXbwwO6dJHZbqNGGXa2Gx8NwEplBjExUiUotjTSzp1SPENb/znMy40r8OJ7IVz8uzT2MC98yp2lUJNDfPpcKcqUgbPXz/L4/AGsP70CDGC58YbfZD56rL+jh548Bw5IUt+sWXDpEpQpI28Ee/eGAgUS++wMSQOwUhnZ2bNSim/qVFlW3rVLauHOng3u+uvvCFFREoPef9+X86d8yVn2MlnqbaZU5UiGtqxA/fIevP7b63y5aSKR0RHySRZgYObGDdQr0uY/mdJOKzJSWv75+0sLQHd36NABnnsOmjbN9Kst+huoVEZ04gS89hosXy6Fgps2ldlG+fJyvwbfdGe3S9Gm996TE1716sn7oCZN8gAPczH0IuPXj+fpX74iPDqcvFYzTGRNLnl+hjHRWLhji36A8SsDnT8AHz4sb/pmzIALF6BkSWlF2bcvFCrk6NE5Df0tVCqjOH9emtvXqCFnVTZtgpdflupAsYFXpTtjYMUKSajatQuqVZOPH3tMtjt/Pfgr49aPY1PQJsKjw+letTvvPfIerT45hAHcI/MT7rYLb3tVvOyV/pO85TSioqTdpL8//P77nUfYWrTQpL54aABWypUZI0t7/v6y1FelCmzbJgH46NFMv8TnSMu3B/HOpPMc/KkkkWdyU7hYNPPnu9O1qzwtV8Ku8NpvrzFzx0xA6jTP6TyHp6o+BUCRXEEEBYfhZa+Elz2uYXyRXE7WUerYsbjZ7tmzUKyYHmFLIg3ASrmqhQulRNLBg1KgYPBgme3G0uDrMB99c573R3oTeqwmtuxh5Gn1LzlqnsG7QhWuRWbjs38+Y8I/E7gacfWOzzsefPzWv4e2rHBHoQ4AHw8bQ1tWSLfHkaDoaPjpJ3nj9+uvMpVv00Zmu61b62w3iTQAK+UqjJFjG1WqSE3C0FCpzvDOO1K0QHvtOty//8p7oh9+KIBblghyN9tD9honsNzthNpDeemnd7n223cEhwfTqWInOlboyHM/PUdkTOR/GiXE7vMmVqoyXZ06BdOmyZ+gIChSRB5wv35QvLjjxuWiLGNMun0zPz8/s2XLlnT7fkplCJcvS7bOlClSuOCzz2Rv15hMdWbSmR06JKUi58+X3gBWtUCy+x0lyns3YW7bibGuEGr7G7t1jfYV2jOy0UhqFq4JSD1npy6kERMjs1x/f5n1GgMtW8pst21bTehLhGVZW40xfvHdp/9zSjmr6GjJGl20SIoXPPSQNEbo2lXu1+CbppLSKOHUKfjgA5g+XaonDh8Or78O7aYEcShkG+c93wNiwALPmEpU9H6R77sNuuNrOG2jhNOn446wnTghZUqHDYP+/aFUKUePLkNINABbljUDaAucN8ZUuXnbeKAdEAkcBvoYY4LTcJxKZQ5XrsDff0v2qLu7tF7r10/2dmN77qo0l1ijhAsXYOxY6Yhnt8v2+1tvyQmb0KhQypf9k4DdH4N1c//WWOTkQUa1bu+oh5Q0djusWiVVqn74QWa/zZtLecj27aUblko1iS5BW5b1CHAdmH1bAH4UWGOMibYs6yMAY8ywxL6ZLkErFQ9jpAD9lCky242KgjNn4qrvq3TXYNwaguI57lPQOxst7Y2YMEG24Hv1kqXnEiUgLCqMKVunMG79OM5eP0vZnNU5fHUvxsTgZnkw5uFFDG/mpAH43DlZXZk6FY4cgXz5oE8feeNXtqyjR+fSUrQEbYxZZ1lWybtu++22D/8BHk/RCJXKrDZulBrMu3dD9uxSlm/gQA2+Dnb3WVt7lBvXtpXk5D9l2BQuuwCjRsGVrAHMPrKK4P3BLNizgNPXTtOkZBMWPb6IhiUaOvf+7t1H2KKioHFjKZjRqZO2nUwHqbEH3BdYmNCdlmUNAAYAFNcsOZXZGSM1mT08wM9PskizZJGZR7dukC2bo0eoiGuUYGIsru8sztUNZYm54U2uChdZMz8fNWvCuuPraD6rOVH2KACqF6zOnE5zaFKqya2v45T7uxcvSi3MKVPkCFuePPDCCzLbrVjR0aPLVFIUgC3LehuIBuYmdI0xZgowBWQJOiXfTymXFRwMc+bIi96uXVIPd/lyKVqwcaOjR5fhJCWB6l5ebV6BF9+/xPl1ZYm5mgWvopfw7bKDz18pRpVqUUzdOouhvw+9FXzdLDe6PtD1juDrVGKPsPn7w5IlUqP54YelLubjj0svaJXukh2ALcvqhSRnNTPpeZZJKVfz3nvwySfSdLx2bXkR7N7d0aPKsBJLoLoXY2DZMnj3XV/O7PMla5EQsj66iTI1r/Pqo6W5zEoqfDWao8FHqZyvMuHR4UTbo/G0edKkpBMG3ytX5Aibv78cYcuZU7Y4Bg6EBx5w9OgyvWQFYMuyWgHDgEbGmNDUHZJSLi44WA6E9ukjM4sCBaBnT1niq13b0aPL8MavDLyjehRAWFTMPZsYGAMrV0pNk61bpVPjkiVQuM4e/ji2hoiYCIavH8DhK4fxK+LHV22+onXZ1vxz6h/n2+M1BgICJOguWgTh4fDgg3Kk6MknZctDOYWkHEOaDzQG8lmWdQoYAbwJeAG/W3IW8R9jzHNpOE6lnJsxspQ8ZYq0vAkLg6JF5TjRCy84enSZSkLNCu6+PXaZ+shub0LXVybkWC5KloRvvoGnn4aAoL9pMrsZkTGRAJTPU54fuv1A2/Jtufm651x7vFevwrffxm1zZM8ubwIHDtQjbE4qKVnQ8a2VTU+DsSjlmi5cgGbN5EUvWzbo0UNe9HS26xCxCVTx3R5r+fYghkw8ztk/qhB+pAC2bOEUbLWHj0flokudwizas4iXfnnpVvB1s9zoVaMX7Sq0S7fHkSTGwObNMttdsEDORtWuLUG4e3dN6nNyWglLqfsVO9s9eFCWlvPlk/rML74omczZszt6hJlaYk0M9u+HZ3u5c2lXfdy8I8nVeB/Zax3D8ojmnd+3MGrLUvZc2EOpXKW4GnGVGHuM8+3xXrsGc+dK4N2xA7JmlWm7vvFzKRqAlUqquzOZixSRWYa7O8yb5+jRqZsSamJQM48vffvKErNxz0vO+gfIUfcollckYW7/EOwxj6joY1QylVj4+EIer/w4G09tdK493q1bJejOmwc3bsjS8qRJEnxz5HD06NR90mYMSiXF7NkyuwgPl/O7AwboEp+LOHtWakv4+0uHxuefh798/uJc9FVCbEsJcf8eu9sV3O1FKeXZi31vjsLm5kTt9K5fl+Vlf3/YskW6Xj35pPw8Pvig1gR3ctqMQan7deWKJLQ0aCBLetWrS93BAQOgVi1Hj04lwZUr8PHH8MUX0suiXz/pnOfra3jn1yOM3/g2UVYQGMC4U8T+Mh+36Z7k4JvSs8aJ2rlTgu6cObLk/MAD8mB69pT+z8rlaQBWKpYxsH69LDEvXiyz3ffeiwvAkyc7eoQqCa5fh88/h/HjISQEnnoKRo6EMmUMvx76lc7TRrD59GayeeYiKsoCywB2Hql6KckBNCVnje8pNFSODvn7S8U0Ly/p9TxwoLwZ1NluhuLm6AEo5TQaN4aGDaVCVZ8+sG2bFPxVLiE8XAJv6dJynrdRI8lP+vZbwxHrN+rPqE+beW04f+M809pN4+enf8DH3RubZcPH3YvB9ZLeKOFeZ42TZe9e6fHs6ys/e1euwKefStP7b7+VqlUafDMcnQGrzCm2NN+yZfJC5+YmFfafeUYymbNmdfQI1T3cvvxbOHsWakbU4IeZuTl1Sk6EjR4N9iIb+Hz7dDav38yu87solqMY/m396V2jN542aau3+pnVyUqySupZ43sKD5dqH/7+0oLSw0PKQg4cCI88ogE3E9AArDKXixcloWrKFAgMlMzRwYOhfHnJzlFOL3b5NzQyhtB9hdn8d3kCrmSjfJVIVs3ypFkzmLRpEi/OehG7sWNh8Vq91xjTdAxe7nd2+EluIY2knDVOUGCg/PzNmgWXL0u7v48/lk5Y2gUrU9EArDKPLVtkHy0yEurXlxfAJ57Q0nwu5uNfA7m0Ly/B6yoQdSEHHvlCyN95M3n9ruFVzpNms0ew5uiaW9e7WW7k9cn7n+CbEomdNf6PiAhp+efvD2vXytG1Tp1kttukiazAqExHA7DKuM6flyCbLZvMcqtXh1dflTOTVao4enQqGf74A7Z9VZOI07lxz3WDvG23k7XSaSLd97E9bB4NZ26nYNaCvPzgy0zZOoXImEg8bZ40Ltk4VceR0Fnj/yRgHT4ss92ZM6ViWqlS8OGHss9bqFCqjkm5Hj0HrDIWux1Wr5b+usuXS5PxJ5+Uc5TKZW3aBG+/DatWgWeOcLLVO4BH9d+54bGGSLfDRNoCcScX41q8w6A6g8jikYWAkwGOKaIRFQU//CBZ86tWgc0G7dvLbLdFC53tZjJ6DlhlHs8/Ly98sU3G+/eX1jbKJe3eLWd3ly+Xip+ffgq+9S7x8k+zOGv7TI4QGcgV05av231Ktzrlbn1uujdKOHZM3vjNmCHVP4oVg/fflwPIRYqk3ziUy9AArFxXTAz8/rss8Y0bJ4lUffpIBmmnTtpk3IUdPgwjRkjFxezZ4YMP5JTO4Rs7GLl2JGfdv4+72HKjbZUqdwTfdBMdDT/9JHu7v/4qmcuPPSaz3VatZParVAI0ACvXc+qUzDKmT4cTJ2RqdPCgBOC6deWPcklBQRJsp0+XUzlDh8KwYRAUtYvev4xk2b5l5PLORf9a/Znz75xbe7z3c4Y3VZw8CdOmyUCDgmSG++678OyzMvNVKgl0D1i5lhs3oGBB+bt5cykN2aEDeHo6emQqBS5ckEWMiRNlG79F5+scKTOHIOsXLM8ThJgd5PDKwZCHhvDKQ6+QyztX+u/xxsTILNffX2a9xkDLljLbbdtWMpuVuovuASvXdfy4zDL27IGlS6VAxowZ0hChdGlHj06l0NWrsq/76adShfGZZ6De42cZ+c9nnLGNB+xghxymBV81+ZSeD8Zlr6fbHu/p0/IzOG2arLgULAjDh0t+QcmSaf/9VYalAVg5n6goWLFCElpWrpTbWrWSV+gsWaRilXJpoaEy2x03TmpRPP645CvZChygwaTBXLStAQxYgHHDLaYwk/84T88H02mAdrvkF/j7S0ZzTIysuPzvf7Li4uGRTgNRGZkGYOU8jJEkltmzZS+taFHZV+vbF0qUcPToVCqIjJSJ5OjRcOaMvK8aPRpyljrEB+s+YM7iORi7B1liHiHMFoAx0Vi4422ven9lHpPr3DlZYZk6FY4elfyCV1+VrY6yZdP++6tMRQOwcqzYCkFTp0KXLlIw44knZJmvVSvdV8sgYmJg7lzJbD52DKo9FkDbJ9fSsnZ5Jh38mW9++gYPmwevPPgK67bV48JVHyJi9hHutgtve1W87JWSVuYxOex2qfAxebKcd4qOlsYcH34o2fReqVdBS6nb6aubcoz9+yXofvMNXLoke2k+N19gc+SQpBbl8oyR91fvvisNf2rWhJc+DuCtA03ZdSSCqUcMHm4evFD3BYY1GEbh7IVZXkBqPRNVCS+7nOG+Z5nH5LpwQSqlTZkChw7J2fEXX5TZbsWKqfu9lIqHBmCVfqKj42a0AwZAQAB07Cj/btZMKwRlIMbIFurbb0sJ7ooVpcVynWYnaTF7MOHR4beubVe6P5+1+uzWx0ku85jcga1bJ7PdZctkTfzhh2Vq/vjjenZcpSsNwCrt7dwps93Fi6W0Uf78MGkSFCggf1SGsn69BN4//5St+5kzoUn7IMb/M5buX00lOiaG2FbkFu5s3V+W5duD7giwHWv6pk7AjXX5sqy2TJkiqy85c8rxoQEDtC64chgNwCpt3LghZYymToXNm2UfrUsXSX8FfdHLgHbsgHfekSOyuasF0O7jtfRuXYV1p1bz3KTJxJgY8tIS94guxFiXbu3vYi/P+JWBqRtwQWa7GzZIJvOiRZJv8NBD8o6ga1ftgqUcTgOwSj3GSODNlk2ySQcMgAcegM8+g549ZY9NZTiBgfDeexLjcueGQR8GMNPelB9DI1ix1OCGG71r9OadR96h6Ud7MYC7KXBrfxfus5F9YoKD4dtvJfDu2SO1LPv1k5/H6tVT7/solUIagFXKXbkiKa5Tp0oy1fffS5GMPXukEYJlOXqE6qbl24NSbW/1xAkYNUrymHx8ZPbba9AF+q4cRvgJ2eO1sBhSbwifPPoJAEVyHU1+I/t7MUZaJvn7S+ersDAp1jJ1KnTrJm8KlXIymvWikm/TJildVKSIZI+6u9+ZvVy5sgZfJ7J8u2QXBwWHYYCg4DDeXLaL5duD7uvrnDsnjRHKlYM5c+Cll2DL3ktEN3qTGrNK8feJv7FZNmyWDW93b7pU6nLrc4e2rICPx50NClKU4RwSAl9/LenVDz0k0/CePWHrVtn6ePZZDb7KaekMWN2fCxckgcXTU6pULV8OvXtLWb5atRw9OnUP41cGEhYVc8dtYVExSd5/vXIFPvlEdhQiIuRpf3n4ZRae+JQ6cz/nRuQNulXpxnuN3uNK2JV46zSnWobz1q0y2503T7Y9atSQQPzUU3KMTSkXoM0YVOLubnI/b54c2QgJkXZrWbM6eoQqCUoN/4n4ftst4Oi4xxL8vOvX4YsvYPx42V7t3h1aPPcbc459xD+n/iE0KpQnKj/BiEYjeKDAA2k1fBnIggVyhGjrVln37t5dspnr1NHVFuWUtBmDSp6ICHnVnT5dyhflzStN7mvUkPt1puFSiuTySdL+a+w+cdDFCNwCy3B5QxmuXrbRrh0Mey+EWadfo+8f0wCwWTZmd5xNz+o9027gO3fKbHfOHLh2TTLov/wSevSAXLnS7vsqlcY0AKs7RUfLOckqVWSZec4cKFMGxo7VsnwubmjLCry5bNcdy9B3778u3x7E8CW7ubi9EMHryxETkoUsJS4x6pOruJddQLvVn3Al/Mqt62OM4bfAvfS8Lbk4VRK9QkNlP9ffH/75RwpkdO0qs9169XS2qzIEDcBKHD4sRehnzpQM0tOnZYlv61ZdYs4gEtt/tdth2PgrHFnZgOjL2fAsfIXcbf4husy3fHBiGdEnQvAr2JxT16tyzjYRgzRK+GtXAZaXk0IasYlesUE+NtHr9u9/T3v2SND99ltZ765YESZMkGQ/PcamMhgNwJndpk3w1luyx+vmBq1bS0JVbLs1Db4ZSnwVpoyBn3+W6lUHdlbBI18IOXtMJqLkNC7Z9mGsUHxiarN+4CSGfHsdr6gwCsYUibeQRrISvcLDYckSCbx//y0rL126yGz3kUd0tqsyLA3AmdGePTK7LV1aXn0PHZIDnX37SgtAlWn8+ae8/9qwQXYayjy5lXNlRnHVY8XNXrwWeSKfp2L2ztT1rcvp4J8A8LJXireQRkIFNeK9PTBQgu4330ipyHLlJOegd29pA6hUejt9Wo5VppNEzwFbljXDsqzzlmXtvu22PJZl/W5Z1sGbf+dO22GqFLtxQ5aY69eX/d1x4+T2unXhyBEpZaTBN9PYsgVatpSue8ePw1eTw3lh9pdcrvEY1z1X3Halhc1249Y+cUIFM2JvT+x+IiIkk7lJE1le/vJLacSxerXkHrz+ugZflb7OnpXAC/JaePx4un3rpBTimAW0uuu24cBqY0w5YPXNj5WzGjYMCheWcnyxhznHjJH7LEu7EGUie/bI6m6dOrK9/9EnEby+aBJjb5RlyO8vUbVQBfpXHY2b5QXGDTfLg1cbdbq1fJxYIY2E7h9ZyRPeeEPe5HXvLi9yH34Ip05JslXTpvpzqNJPdDT88AN06CA/kx9+KLc3aCDNYtJJokvQxph1lmWVvOvmDkDjm//+BlgLDEvNgakUuHIlrkBG7P5Z586yt1u/vu6pZUJHjsDIkZLU7lMugCYjV1PN7xpfHZjPyd9P8nDxh5ndaTZNSjbBsiz61GmarEIat99/4VIIXU9v56VDaygw+m85M96+vezttmihAVc5xujRMHGizHwLFoTXXpPtN5DXxnRs0pGkQhw3A/CPxpgqNz8ONsbkuu3+K8aYRJehtRBHGortczptmiS0hIdLglWdOo4emUqi1KzTHCsoSF5vpk2TSqEdX/6LpVmbEWWPAqBK/ip82vJTmpdujpVab8yOHpWiLTNmSN3KYsWkEULfvum6v6YUINtvv/wiSz+WBYMHy8pLv37Qpk1cwmkacWghDsuyBgADAIoXL57W3y5zOnBAajAfPChlIvv2lRq4NWs6emQqiVJ8fOe2rzN+ZSAnT0cTs6MCFzcVxx5j8eyAaMo/MYcPNr5KVLgEXzfLje5Vu9OiTIuUP4DoaFixQpKqfvtNXujatIHnnoNWrWT2q1R6MUZqgU+bJjkH167Jx35+Mvt1klXA5Abgc5ZlFTbGnLEsqzBwPqELjTFTgCkgM+Bkfj91u5gYqcMcGiolIUuWlK5D774r7/K0z6nLSWmdZpDg+8b8vZzfUIKQzaUwke5kr3qSLq9+z6prXzD5z0NUyFuBG1E3iLHH4GnzpEnJJikb+IkT8iI3bRqcOSMz3HfflTeAxYql7GsrlRyBgfK6uHu3nPZ44gmZ7dauLfc7SfCF5AfgH4BewLibf3+faiNSCTt2LK5YxqlT8OCD8oPm6SktAJXLuq/jO/EIC4NX3g7l5NpG2PNuxf2xmXiXOEtYjlXMOn6K6gWr832372lXvh3/nPon3v3dJIuJkSW9yZPlb2Nklvv11/DYY7LWrVR6iYmBVaskw759eyheXPZ2X3xRWlE6ccncRH9TLMuajyRc5bMs6xQwAgm8iyzL6gecAJ5Iy0EqJIPm/ffl361aweef39n6T7m0+63THLtPPKRpBS5u8+WDD+D06XJ41PkJe5uORBPNdQts9kLkj3iLbQM/wM2SpKd6xeolL/AGBUld8GnT4ORJKFQI3nxTZrslSybnYSuVfMePy2Rk5kxZiXnoIQnAPj4SkF1AUrKguydwV7NUHou63Z498mL34otQqpRUBBo5Evr00aW9DCipdZpjrzF2CFyfh6c+yk1UMNRvYMfW5X8E5f4QrGj5BGORLaYF5XI0uxV875vdLnu6/v6yxxsTIxnMEybIi10aJ7AoFa933ok7OtSihRRw6dDBsWNKBl0rcibXr8uZyGnTICBAXtzq1JEA3LSp/FEZUlL65I5fGUhoZAxhBwsS/FcFoi5mx71AMEUH+XO92ixOnvsXd5Mfu3EH7Fi4k8uqmbxm92fPynbH1Kmy9ZE/vxTJ6N9fSmYplZ7+/VcmJEOHyrndhx+GESPkqGWJEo4eXbJpP2BnERYmP1iXL0tC1bPPQs+e6XooXDkvY6Bwt41c+asCkbZ9WJVX4lM8mKgiPxLldphyecoxotEIvKMbMurX7zkZuoViWfx4v3XnpGdR2+2wZo3Mdpcvl8zmJk3k3G6nTpJroFR6CQmB+fMl8G7eLD9/c+dK3osL0X7AzujiRamKsHu3zHh9fOTAZvXq2m5N3WHDBmmUcG7tg1gV1sCTrTBWFKEW2Ox5KOs+jL3Pj8bdTX6du9R68f6+wYULMGsWTJkidcHz5oWXX5azu+XLp/4DUuoud+c2DG9YlHZt6kgQrlJFtjx69MhwZUo1AKcnu11q3k6fDt99B5GRUos5NFSODg0a5OgRKieyc6dsdf34I+QvYGg1fBmrPAcR7SbneDEWuU1bxrd58VbwTTJjpBODvz8sWyY/iw0bSp5Bly7Sf1epdLB8exCfzF5H2x2/UeLKGd5q9SJv/HqE4oPeoHrnFrINl0EnJBqA09PMmbK0nCePFCjo1w+qVXP0qJSTOXBAtrcWLICcuQx9R69mb4ER/Hp6A7k883E10h1j7LhZHrzWqMv9Vcu6fFm6D/n7y3nJXLnkjd+AAVC5cpo9JqX+IzoaVq4k9xtj+WP/P3jYY9hY9AE8o6MIAwbnqsf6unVT7dulRaW5lNIAnFYiI2XqMm2azCj69ZN6zFmzQseOOsNQ/3HihJw0mzULvLzg6bfXcqTEe8w4/RdFrxfl68e+pm/Nvmw9vfX+zvEaA+vXS9BdvFjOS9arJ9+oa1fZ/lAqvfn7wwsv8ECWXEz368Diai04nDfuhEdSz8AnRWpVmkttGoBT2/79ssT8zTeyt+brK4EXIHduORiuMqWE3oGfPy8nKr7+Guy+AVQcOgu3oluYe3EbhUMK82XrL3m21rN4u8ubtiSf4w0Ohm+/lRe6PXukIMGzz0pSVdWqaftglbpdeLhsu02bBr16wTPPwJNPQpEidN6blRPXov7zKQm1tkyO1Kg0lxY0AKeG6Oi46j+9e0uft3bt5MWuZUutg6vifQf+xrx9zJ2YnV8W5CAsDOq9MIWAPIPYgx0uwssPvszYZmPx8biPFyJjYONGCboLF0p2fZ068sLXrZuswCiVXv79V3725syRLm2lSsV1wcqXDzp14tWSQYmegU+plFaaSysagJPr9mLfy5fLflru3JJJWrCg/FHqptvfgdsjbVzbVpKQf8pwMMKD5r02El5vBH+fXXnreptlo2DWgkkPviEhckRj8mR50cuWTY6xDRwItWqlxUNSKn5RUXEFWnr1gr1747bhmjT5TxvKpJyBT6mkVppLbxqA71dwsCwvT5smR4iyZJFi3zduSADWpCoVj9PBYZhoN67tLM7VgDLYb3jjVWclNBnJqiz/kPdqXgb7DWbmjplExkTiafOkccnGiX/hLVtktjtvnmTT16wpQfippyB79jR/XEoBcXkG06bBzz/LcbYcOeS10tdXjrbdQ8eavmm6FJyUSnOOoAE4Kex2Cbx58sD58/DKK3J8yN/f6Yt9q9ST3CzK6GiwHSrNyd9LEJNjJ25tXsS91EYisuzCnex82PRDXqj7Atm9stOjWo/EE6yuXZMCBf7+sG2bvAns1k1muxn4yIZyQpcvS8W0adNkFTB7dujeXd4M5sjhNBOS9JhlJ4dWwrqXEyfiin37+Umje5B3d2XLOnZsKl3dvYcL8g56bOeqCf4S2+2wdKl05wsMBPeGk4lu+jxYdjCQ096aSW0/5am6FZM2iB07JOjOmSNlS6tUkeNsPXpIH2il0kNMDFy9KhOS3bsloa9BA1li7tpV8wzuopWw7tcvv0i3od9+k6WVFi3kXV0sDb6Zzv1kURojP0Jvvy0xs0y93dSbMIqAq0viLrLcaFeleuLB98YNSaby94dNm+T4WteuMtvVimkqPR07FjchadhQcg6qVIHDh6F0aUePziUls0VKBrRnjyQPgDRC2LtXpi5Hj0og7tLFseNTDpXULMp16+S16bHH4KK1jwc/6caRltXYHbaSPjX64OPug82y4ePuxeB67RP+hrt3SycsX1+ZWYSEwGefSUvAb76B+vU1+Kr08fPPcpqjdGn44AMp2PLEbR1oNfgmW+aeAYeESLmh6dNldvH999JibfhwKUWkx4fUTYllUW7dKjPelXsCyFpnMZU+3M3+yFUER2blzYff5LX6r5HHJw/9a/VPeI83LEwKZfj7SwFoT08pPD9woER1DbgqvezbJ3XAbTYpWbp/P7z3nrRDdeHuQ84mc+4Bh4RIsflFiyRZ4IEHZJbRs2eGK/atUkdCe8CDa9Rk7byCLFsG2eou5kab7hjkmh5VezCh1QTyZUnkZ2r/fgm633wjZyXLlZOg26uX/jyq9HPjRlw71A0bZB+lVSu53dtbJyTJpHvAIP1N9+6VnrrZssH27fD00xJ469bV2UUml1iG891ZlHnsucmxowavfJgFnyJHqDFyNDuZhUHe0NosG5XzV044+EZESBMEf3+ZYXh4SMW0gQOhcWP9eVTp5+pVGDZMjrJduwYVKkiD+9q15X5NqkozGTsAR0fL/sX06fDTT3JO9/RpebHbvl1f5BSQ9DqxHWv68mAhX0aPlj71Vp5jPPDmGPZ5zWK/mztdK3Tl+8DviYqJSvgc78GDUqxl1ixpSVm6NIwbJ0t7BQqkw6NVCllpCQyEhx6SCckff0jP5/79JaNZXxvTRcYNwEuXShLLmTNSleq11+RFLrZCi/6AqZuSkuF86RJ8/DF8tjSAqDLLKfFaIEFZfibQshhUexDDHx5OkexFCDgZ8N893shIyS/w95d2lDYbdOggs93mzf9TGUipNGGMZAlOmyZHKnPlgpMnpYzu3r26xOwAGScA37ghP1R160KlSlCokJzd7dcP2rSJC7xK3eVeGc7Xrkkv8P/9D0IKf4/V43GMFc0xoHP5znze+nOK5ih663PuaJRw5IhMlWfMkAIuxYvD6NHQty8ULpz2D0ypWCtWyCTk4EE5M963r9Sqj61hr8HXIVw7AMfWY54+XSoDXbsmR4fef1+WUX74wdEjVC4gvgxne5Qbtv1lKV0aLkacpvQz47hRYBIxRmbKNsuGXxG/O4IvINseK1bIbPe332SlpW1bme1qYw6VXmJiYOVKqFhRtjmyZpU3fe++K0cqs2Rx9AgVrhyAjZGjGevXx9Vj7tcPHn7Y0SNTLub2OrEmxuL6rqKEbChPtAmm+FND8Co6meMmitblWrPqyKr493hPnJClvenTJc+gaFE5ytavn/xbqfRw/LisuMyYAadOyZHKsWMl+bRp0zT5ls7Y6N5VuG4AtizZR3vmGa3HrFKkY01f7HZ4+au1nA46hP2YIXe7mYRW9OeUieCZas/wTsN3KJOnzJ17vEXqxs12f/lF3hS2bi2Nfdu0iVveUyqtGSNnxr/7Tj5+9FEp3NKuXZp+W2dtdO8qMuc5YKVuMkZi6JBPAzjycFNwjwDLYGHRo1oP3n3kXcrlLXfnJwUFyUx32jRJYilUCJ59lpUPteX9XTd0JqDSR2CgnO549VX5+M03wcsrXYtlNBi3Jt4CNb65fFg/PG1m3K5GzwErFY/Vq+Gtt2DTrst493oLPMIBsLB45aFX+LTlp3EXx8TInq6/P/z4o3x82yxj+e7zOhNQaS8sTJJNp02TjGZ3dzk+VKqULDWnM2dtdO8qNACrTOH2farsIQVhSxV27QojR8sJeD/2GeHmGjZLEqQ8bZ48UflmrdszZ2Q/bepU2V/Lnx9ef13OS5Ypc+vr30+zBqWS5Z9/pDLV1avyszd2LPTuLSswDuKsje5dhQZgleHF7lNdDcpC8F+1OXbSB+uRMXgO+4IQK4THKz3OiEYjuBZxTfZ3SzxCvf3XYcjjcn43OloSWD7+GDp2lBrNd9GZgEp1165JrfqcOaUDVtWq8vPXqxc0auQU58edtdG9q9AArDKEe2Vivj/3BCd+rkpo+D5o2AOKB2DcQ8nm1oA1/SdSvVB1+SLnz1Nvwzro/Yyc4c2bV2qGDxgghenvQWcCKlUYI41hpk6V4Hvjhhwbiu2zO2uWo0d4B2dtdO8qNAArl5dQJubFszY2fVeI7d88AI++AbX9wQKMG3kiXiaHvQXVC1aTMnyTJ0sGaVQUPPKItF3r3FmK0CeBzgRUqhgwQPZ3s2SR0x39+8ODDzp6VPfUsaavBtxk0gCsXN7d+68xoZ4EBZSh//+yYPMbj9vr47B7XobbEv69os/w6p4VUPFVOHBAyvINHiwvgJUr3/cYdCag7tvtpSHHj5e93CeekAp+3bvr0cpMQAOwcnmx+6z2CHdCNpXm6o6CUH0q1svjiPG+SM38jThzpgbnbF9jTBSeMbBg4VIeOREjje3fflte+HxStlysMwGVJOfPS+vJadPkzV/OnNIKtVAhyaxXmYYGYOXyCmbJSuCaggSfDgK/IdBoA3hdIYdVi5/7fE/9bJX498MvOLgyDwdynqXuaS+K1usGP74qiS1KpZcrV+SMbni4VO17+20poKGlITMlDcDKZUVGSq5K4MS6BFd/F5769OYer0XB6MEsLNOS+u9OgYULqRYeTrW6daHPGHjySe1xqtJHUFBcWUh/f2mJ+vnnUka3UiVHj045mAZg5fTuznB+tXkFru/xZcT7kRzPMwPPvmPA+9St6y1g4K7FNBo9SXqd9u4tzRBq1HDUQ1CZSXS0lCadMkX6kdvtsrQcFSVd2QYMcPQIlZNIUQC2LGsI8CyS3rIL6GOMCU+NgSkFd2Y4GwMH/slJ9wlZiCo6Dc/HP4AsJ/ArWp8nc3Rl+O7PiTQxeNoNra7nBv/RksySPbujH4bKTL78UspDFiwIb7whDTnKlnX0qJQTSnYAtizLF3gJqGyMCbMsaxHQDZiVSmNTivErAwmNjOHa2bNcvXAMe7gdnp4HOY9Rs5AfoyI78uikv7C2fUqdsl6sbVuTxm0GUW90X0cPXWUGUVHS9nTqVKnB/OST8NRTULKktKHUPuTqHlK6BO0O+FiWFQVkAU6nfEhKxTmyy5uLgWFEtx4MpSLBgqwRBRgUUJuP1+/Dur4FqlWDiROp9/TT1MuZ09FDVpnBoUOSxTxzpmQ1Fy0qSQkgM99OnRw7PuUSkh2AjTFBlmV9ApwAwoDfjDG/3X2dZVkDgAEAxYsXT+63U5nMtm3w1tsxnA06Cq1fBnd5cXOzw1t/nee1f4Kxnu4Ozz0nhQosy8EjVhmeMXE/Z506wb59Msvt319qNNtsjh2fcjnJLiZqWVZuoANQCigCZLUsq8fd1xljphhj/Iwxfvnz50/+SFWmsG8fPP6Endq9FvB7+SrQpQcF7GF4RoMtBjxjLI4X6MiqVdukLN9DD2nwVWlr/3547TXZx71xQ26bPh1OnIDly+GxxzT4qmRJyRJ0c+CoMeYCgGVZy4D6wJzUGJjKXJZtCmDMvD/Yts1gqzYPquyl0jUfRi2E9oej+aqBHzMq58Uq/Bit3+nMY1rwQqWl2LZ/U6fCX39J27+OHeUcb9asULeuo0eoMoCUBOATwEOWZWVBlqCbAVtSZVQq0zh7Fl74aD1LszWFXJHQDIqGWHy0GJ6IKorbwOegVy+G5M3LEEcPVmV8kZHS7WrfPnjmGZn1fvSRdCAqWNDRo1MZTEr2gDdalrUE2AZEA9uBKak1MJVxxNep6JESvnz0seGrn5dib9H/jj3e/mEP8OSXX0Djxrq8rNLejRuwcKHMditWlMSqWrWk/27duvozqNJMirKgjTEjgBGpNBaVAd3dqejkuUieHXKdbJcnE1NvJKFdzuEbAhdiIMZm4enpRdM3p0Cxeg4eucrwdu6ULlhz50rv3YoV71xadvIuRHDvNpzK+WklLJWmYjsVhdsDuX7uKOUOXyZL5RUE+l6j1BWYcao2PZ8YzeaK2Vl7fB2NSzamngZflVauXZM9XDc3SeKbNUt67fbvDw0auNRsN6E2nIAGYRdhGWMSvyqV+Pn5mS1bdJs4Myk59Gc89+3lSM3hxNhiwIL81y0eOfkQ8/+3AA/f4vouXqUtY2DrVikNOW8efP89NGsm53c9PKQ+swtqMG4NQTc7gd3ON5cP64c3dcCIVHwsy9pqjPGL7z6dAas0YY+I4u/hPzD4xzFMarudmNifNANRnt05U6PfreCr7+JVmoiIkEYIU6fC9u1xTe6LFJH7CxRw7PhS6HQ8wfdetyvnk+xzwErFxxw7zoEn32FF5UK8F/w4w3ps54KPF5axgXHDwpOcbn4MbVkBiFuivl1YVAzjVwY6YvjK1RkDZ87Iv93c4P33pRnCpElw+rSc380gXYiK5Iq/f3VCtyvnozNglXLR0fDzz1z60J+9p39mZBNY8wzksefm8xYjKOTZltG//czJ0C0Uy+LH+60735rd6rt4lSqCg2HOHFlmvnIFjh2T5eVt26TRvQvt7SbV0JYV7lg9AvDxsN16c6ucnwZglXynTsG0aVz7YgoLSp1hbAMbR/NCFnsexjd/h+frPoePh7wb7+r3Yrxfokgun3j3sfRdvEqSPXtg/HhYtEiKZ9SuDe++CzExUp2qcGFHjzDNxL6J1fwJ16UBWN2fmBhYuRL8/TE//simQnZ6d87J/mKAiQFjI699KGWzdL0VfO9F38Wr+3blivwc5ssHJ0/C0qVSNKN/fwnAmUjHmr4pDriaBOk4ugeskubMGRgzBsqUgcce48+df1H5ydI8NAD2FwkHY4EFYLhm9iV5D7djTV/Gdq6Kby4fLCSDc2znqvoCoO5kDKxfL4G2SBGZ9QK0aCE/m5MnZ7rgmxpikyCDgsMwxCVBLt8e5OihZQo6A87k7vnu126HVavA3196nkZHs6xqXV7qlp+gilvwNpAztDdetrJc8HwfY6KxcMfbXvW+9nBT4128ysAmT5Ym93v3Qvbs0nf36aflPpsNsmVz7Phc2L2SIPV3Mu1pAM7EEjoC5Hn5Im22/CoJLUeOsP6BnHz6VDV+dXMntOQmPO05Ger3Pm83e4k2n20lKDiMgpFjCHfbhbe9Kl72SrqHq5LPGDk2VKuWfPzXXxJ4p0+XhvdZszp2fBmIJkE6lgbgTOyOd7/GUO/ELp7e8QstxgZATDTRDRrxaqPWfFlsMrhtAwNPlO3DlC6fkss7FxC3h0tUJbzscrxD93BVsly6BLNnyxu//fvh33+halU5y+vl5ejRZUiaBOlYGoAzsdPBYeQOvUqX3avpvnMlZS4HEeydjW+qt+Vss2f54NQcIotPvHW9zc1GzeLlbgVf0ExMlQqCgmDYMGn/FxEhPZ5nzIDSpeV+Db5pRpMgHUsDcGZkDPz9N/6/TqDR7j/xiolmi28lXmndjflWfS5kW4LxaY+trA9tivZkzbnFRMVE4WnzpHHJxv/5crqHq+7bpUtyjK16dVleXrsWnn0WBgyAatUcPbpMQ99AO5bWgs5MrlyRJT5/f9i3j6hsOVhYsRGf1y7Hbi4QHnYOSq7CsnvRrfQLfN51KPmz5ifgZABrj63VRgkqZYyR/Vx/f5ntVq4se70Qd25XqQxGa0FnZsZIX1N/f+l5Gh4ubdZmzMC965Ms+WIRW8KfBUuWoB7w6sSqF76mUPa45uP1itXTwKtSZskSKZCxfz/kzCkz3f794+7X4KsyIQ3AGdXVq1Kaz98fdu2Soxq9e8PAgZjqNfh2xXFee/EVLhabDpYdLLBZNp5+uM4dwVepZLm5zUHlypA3L4SGQq5c0uy+a1dpjKBUJqeFODISY2DzZtlLK1IEXnhB6uH6+0sh+q+/5rtL+fAdMJheW8pxseg3+GXrhLeHNzbLluAer1JJdvkyfPYZPPAAPPKIbHkA9OwJAQHyJlCDr1KAzoAzhmvXpM+pv39c27WnnoKBA8HPj4CTAUz94UN+2XSQszlWQGHDw1n7MbPvW5TNX+yee7xapk4lSXQ09O0rNZlvz2R+8km5PwM2Q1AqpTQAu7A/5q/k6mdf0XzHarJFhnG1XCVyTpwoVYJy5gRg8pofGbyuE4ZoyA2VbO1Z9uwXVCxc4tbXSWiPV3v1qnu6fFmWmdu3B3d3uHED+vWTN36ayaxUojQAu5obN2DBAq5M+JIme3YS5u7FTxUfZm6N1uwvUZmx9arRMWdOtu4/T68pH7Mn2+fgFn1rj7dn44fuCL73omXq1H8YAxs2yGrL4sUQFSXbGwUKSFMEpVSSaQB2Ff/+Ky96c+ZASAjBBUrwebMBLKvSlBDvm7Vwo+2MXLKdDxZ+wTbbV5A9nJIxLTnj+QfR9oTP8SZEy9SpO2zcKDPcPXviajIPHCjBVyl13zQAO7EfAg7x76dTab3hB2qf3k+Mpxe2rk/AwIE0XRGMsSwi3PYR7rYLW3h5rp0/xfFi34J3KGXDnmb6U+/yyAPlk32OV8vUZXKxR9jc3aFOHfD1lcA7bZrs7WoTBKVSRAtxOKO9ezk8+n/k+24ROcOvczhPUebVaMXPNVowrEcDOtb0pcG4NRwJ2cY5z7cxRAHyPHqe6sj3gz6kVe1KKR7G3XvAIGXqtF1gBnf3Ebb27eH77x09KqVckhbicAXh4bKH5u8Pf/1FcZs7v5Svz7warfinWNVbWaSx+6/PNijEwGVLMJ6R0ofXQPbrTzD7+Qm0SqXgqGXqMqGRI+HjjyEsTPrrTpkC3bs7elRKZUg6A3a0AwfkRW7WLKmPW7YsDBhA7ZO+XMqS8z+X2+1h1Cq8jx8ufIrxCga7G1jgZnkw5uFFDG/WPt0fgnJh167B/PlyTtfHByZOhJ07ZW9XG9wrlWI6A3Y2kZHw3Xcy2/3jD9lj69hRXvSaNgU3N7zHrYHb9l9jTCiXr/xJaJ75nAy5TK7g9oxtOZLqtcK1TrO6f1u3ys/fvHmSWV+ggPwMPv+8o0emVKahATg9HT4MU6dKgYILF6BkSRgzRgoYFCp0x6WxbcIuRW8jOGoFUZ57wTcE71NtGN1sFK++53ertsG9Aq8W0lB3uHwZHn1UAnCWLNCtm7zxq1PH0SNTKtPRAJzWoqLghx9ktvH771J0vl07edF79FFwi78a6KNVcjNu1RL235gO3oDdjRZun/HL5JeTXLdeC2koQJaU9+6VvdzcuWWbo08f6NHjVsEWpVT60wCcVo4fl9nu9Olw9iwUKwbvvy+zXd+Eg19YVBhvf+fPxJ3jiPQ8d+t2m82iSePQ+2oao4U0MrHQUOl+5e8v53fz5YMuXcDTExYscPTolFJoM4bUFR0ts902baBUKRg7Fvz8YMUKOHpU2rElEHzDo8N5c/mX5BlVhgn7hmDOV6aL5yR8PHyS3ShBC2lkUosXSzOOvn3lSNGECRAYKMFXKeU0dAacTLfvrVa3rjPu8kYq/rQIgoLkxe/dd6VqUPHi9/w6fx77kw9Wf8rfxzYQ4X4R27lH6FNsHp9/2pjs2SHgZI1kJ1lpIY1MIjxc+u0+8ADUrAkVKsBjj8k2R8OG2ghBKSelATgZlm8P4u0lO6hzYDMjd/xK08ObsYzhXP1GFPzqK2jbVjKb7yEyJpJXfniXr3eOB8uAmxvNIr5gwdgXyJcv7gUzoUYJSRGbyHV3IY2hLSsk6+spJxMYGHeE7fJlGDJEAnC1ajB3rqNHp5RKhAbg+3XmDOeHvcvKjT9SNOQCF7LmYvKDXZhfvSWmZCnWd2x669L4MpAfq1aAietnM2L1B4S4Hb91rc1m0az59TuCb0ppIY0MrGtXWWp2d4dOneC556BJE0ePSil1H1IUgC3LygVMA6ogtRD7GmMCUmFczsVuh1WrJKHl++8ZEBPDXyVqMKZJP1aVe5AomwcA1m3LvXdnIJ8Kvs5ziz/j+rIl3HA/Bmf88PN6id0F3yHKHpmsPd6k6FjTVwNuRnDkiCRPDR8umfMPPiiz3T59/nOETSnlGlI6A/4c+NUY87hlWZ5AllQYk/M4f17O7E6dKi+A+fLBq6/S1V6VTe55/nP57XursRnI4W57CHFbQXjMQYzXOThTk/oRXzLl9cc4GHma937x4GT4Fopl8ePcxeJQLD0foHJqUVHw448weTL89pscYWvdWgLva685enRKqRRKdgC2LCsH8AjQG8AYEwlEps6wHMhul+pU/v6wfLm8CDZqBKNHQ+fO4OXFU9uD2JXI3mpQ8HWuuC0kxHO+7PG6W9i2vEP+PI+y/puGt82QS5OT0oSEoGd0VZzAQKmKdvo0FC0Ko0ZJUt89jrAppVxLSmbApYELwEzLsqoDW4GXjTE3UmVk6e3iRUlmmTIFDh6UggUvvAADBkDFindceq+9Vbuxs2j3Ek5bbxLldSS2SRFgkb36CUpni7r1uXpGV90SEwMrV8qxoe7doUwZCcBdu8qsN5GkPqWU60nJb7U7UAt40Riz0bKsz4HhwLu3X2RZ1gBgAEDxRI7kpDtjYN06me0uXSo1mhs0kCNEjz8uxekTcPfeqt3YWbJ3Ka+tGMmJ8N0QUhnbzneIqTMeiMLCnZxWjVuzZD2jqwA4c0a2OaZMgRMnoFYtKQ/p7g7ffuvo0Sml0lBKAvAp4JQxZuPNj5cgAfgOxpgpwBSQbkgp+H6p5/JlmD1bAu/+/VKOb+BAme1WqZLkLxNwMoA/jv2Bm2VjasB8joTuhIsVKBw4nwnPPoFH37OMXFmAk6Gyx/t+6863grae0VV8+ikMGyYFXJo3h//9Dzp00HO7SmUSyQ7AxpizlmWdtCyrgjEmEGgG7E29oaUyY2DDBgm6ixZBRIRkks6YAU8+KYXp78OGExtoMrsJkTE3t72vFiXv9jl81LMbvSbYbq4Y+tK59ovxfr6e0c2EYrc5OnSAcuWkStorr8gbv3LlHD06pVQ6S+nG0ovA3JsZ0EeAPikfUioLDoY5cyTw7t4N2bPL0Y2BA6FGjfv+csYYfj30K72X9Y8LvnY3Hiv0HEt/fRovr6R9HT2jm0kYA+vXSybz4sWyzeHhAS+/DI88In+UUplSigKwMWYHEG+jYYcyBjZtkqC7YAGEhUlz8SlTJMElW7ZkfEnD70d+541f3mPnpY0QUgiyemDZ7Hh7evJ296ZJDr6x9IxuBme3yyrLli2QI4e86Rs4UEpGKqUyvYyVWnntmpTg8/eHHTsga1ZpuTZwoATg+xS7x5vdMzuzty1ky/n1EFwcz41TeLVpL5p038rWi8mr06wyqC1b4Ndf4Z13pGBGp04waJBsc2TN6ujRKaWciGVM+uVF+fn5mS1btqT+F962TYLuvHlw/TpUry6l+Z56SmYeyRBwMoAm3zQhIiZCbriRD9tf7zPoob68+6YXBQqk4viVa7txA+bPl2Xm2Eb3Bw7omV2lFJZlbTXGxLtS7Loz4NgXPX9/mXX4+MjxjYEDoW7dFGWS/nX8L55Z1jsu+NrdqBX9Et8tHJRYcyOV2axfL+0nQ0Ikg/6rr7TRvVIqSVw3AL/9Nnz+ueynffmlvOjlypWiL7nh5AbeWT2CP46vgrDc4OkBbna8PTz5akhzimuZSBUeLmfGs2aFjh1ltaVLF6lSVb++HiFSSiWZ6y5BHz4MZ8+m6EUv4GQAa4+tJX+W/Czas4Tfj67ECsuPWTeM1vkH0e2VnQS56x6vAg4dktWWmTPh0iVo3x6+/97Ro1JKObmMuQRdpoz8Saa793itiByw7iMe9nqej/6XlXr1AOrd/KMytaFD4ZNPpDpVhw6SVKWt/5RSKeS6ATgFtp/ZTt/v+962x2tR+PgQvhn5Bs2a6SpipnfqFEybBoMHQ4ECclY3Z05ZZi5c2NGjU0plEBk2AC/fHvSfIheli1xi5NqRfLf/O9yisoGbO1gGT5sni8e2pL4mWGVedru0/Js8GVaskLPklSrJ8aF27eSPUkqlIpcMwPEF19sLWsS1+pMyj0ev7uOZ5SO45vY3tqgc8PdIip19mV6v78Oz/FqaltI93kwtNBSqVZO8ggIFpD7zgAFQsqSjR6aUysBcLgDfHVyDgsP+00c3ttXfDbc/ueq+mCi3YxCVDQLeJf/RIYwcnpu8NYOYsCaM04eq8XOuMIa2DNKqVJmFMfDXX1It7fXX5dxu166S0dypE3h6OnqESqlMwOUCcFL66B6/eohLHpMJt22XC+zu8N08chcqw5G9uVm5P/EgntgsW7mgq1elC9bkybB3L+TJI+fGs2eHDz909OiUUpmMm6MHcL/u1Uf30OVD9FreiyCvQYRbuwALLMAy5Gy/jEqPnsXH595BHOJm2UHBYRjiAvTy7UFp++BU2vnpJyhSBF56Sc7wTp8OJ09K8FVKKQdwuQAcX7/cKOssoVm/ouJXFZm3YzHWP6/AomVg9wLjhmXZyOVR7Varv3sFcUg8QCsXEBoqZ3ZXr5aPa9eW0qSbN8vSc9++992CUimlUpPLLUHH9tENjtlNqNs/RFknCbNtxS3aHWvri5i/h/HcU4Xwe/UMn235mJOhWyiWxY/3W3e+tYRcJJcPQfEE4djgnliAVk5s/35ZYv7mG2lF2acPNGsGhQrB1KmOHp1SSt3icgG4Y01f9l/eylt/D8cYmaXaDnQh5scv6NmxCCO3QenSAIXp++iL8X6N2CB++yzXx8N2a4acWIBWTqp/fzm/6+Eh5SEHDYKGDR09KqWUipfLBWAA47kHMLK/a7dRMUdtFgYUSXKb1duzpeNLskosQCsnceIEzJgBb7why8kNG8q7r759oWBBR49OKaXuySUDcOOSjfGwvIiMicTL3ZOpbzfmgftslNCxpm+CWc2JBWjlQHY7rFwJX38tiVXGSD3wRx+FZ55x9OiUUirJXLYZw4YTAfx5XBslZCoXLsCDD8LRo1Iw49lnpWBGiRKOHplSSsUrQzZjqF+8HvWLa+DN0IyRfrv790uwzZcPWraURggdO2rBDKWUS3PZAKwysJAQmDNHlpl375YGCL16SXLV1187enRKKZUqXO4csMrgFi4EX194/nmZ4U6dCgcPSvBVSqkMRGfAyrHCw2HxYqhcWYplVK0Kjz8uR4jq1NHekEqpDEsDsHKMw4fB31+OEV26JCUia9eWQDxzpqNHp5RSaU4DsEp/vXrBt9+Cmxt06CCz3aZNHT0qpZRKV7oHrNLe2bMwYQLE3CxsUq0avPceHD8OS5dC8+YSjJVSKhPRGbBKG8bAunWStbx0KURHyxne+vXhtdccPTqllHI4DcAq9Z04Aa1bS8/dXLngxRfhueegfHlHj0wppZyGBmCVOnbskCXlDh3kGFG5cvD66/Dkk9r2Tyml4qEBWCVfeDgsWQKTJkFAAJQsCe3agc0Gy5c7enRKKeXUNPNFJc/8+VC0KPTsKceIJkyAbds0mUoppZJIZ8AqaWJi4JdfoFIlKFNGlpkbNYLBg+UIkRbMUEqp+6LTFXVv58/D2LESdNu1k9KQAI88ItnNzZpp8FVKqWTQGbBKWP/+8M03EBUlHYg++USSrJRSSqWYzoBVnOvXpRlCrBw5pErV3r2wZo3UaNamCEoplSpSPAO2LMsGbAGCjDFtUz4kle727pVM5tmz4do1qFgRqleH//3P0SNTSqkMKzVmwC8D+1Lh66j0duQING4MDzwge7sdO8pxomrVHD0ypZTK8FIUgC3LKgo8BkxLneGoNHfqlARZgEKFZMb70Udy++zZ8NBDmlSllFLpIKVL0J8BbwDZE7rAsqwBwACA4sWLp/DbqWQxBlavlmXmH36QkpB79kiFqq1bHT06pZTKlJI9A7Ysqy1w3hhzz1dwY8wUY4yfMcYvf/78yf12KrmWLZM93RYtpDnC66/DTz/pLFcppRwsJTPgBkB7y7LaAN5ADsuy5hhjeqTO0FSybd8uVary55cuRHnyyPLyE0+At7ejR6eUUooUzICNMW8aY4oaY0oC3YA1GnwdKDxcmtzXqwe1aoG/v9z+xBOy59uzpwZfpZRyInoO2NUZA2+9JTPeZ56By5fhs8/ghRfkfl1qVkopp5QqlbCMMWuBtanxtVQS2O2SPFWnjgTYffukNOTzz2tdZqWUchFaitKVXLoEM2fC11/D0aNw6BCULi0tAW02R49OKaXUfdAlaFdw8iT06SPLzEOHyt+x7QBBg69SSrkgnQE7q7AwOHsWSpUCLy85v9unj7T/q1LF0aNTSimVQhqAnc3hwzB5MsyYISUi162DAgXgzBnw9HT06JRSSqUSDcDOYt06GDcOfv1VlpQ7dZLZbiwNvkoplaFoAHakixelHGSWLPDvv7BjB4wYIX14ixRx9OiUUkqlIU3CSm/GwMaN0KuXJFHNmSO3P/ssHD8uAViDr1JKZXg6A04vxsCsWTBxopzhzZZNgm6jRnK/VqlSSqlMRQNwWrt8WWoxW5YkV4WFSRDu2ROyJ9hESimlVAanATgtxMRIMtXEifDnn3DiBOTNK12I8ubVSlVKKaU0AKeq4GCYOjWuUlXhwlI4w+3mVnu+fA4dnlJKKeehATg1hIZKJvO5c/DGG1KXedw4OUrk4eHo0SmllHJCGoCTKzwcFi6UZebixaUec4UKcOSIVK9SSiml7kGPId2vY8dg2DA5QtS7N1y/Ds2axd2vwVcppVQS6Aw4Kex2+dvNTUpEfvIJdOwo7f+aNNGkKqWUUvdNZ8D3cuUKTJggS8srVshtL70ks+ClS7X3rlJKqWTTGXB8duyQvd25c+Xcbv36cWd2NZNZKaVUKtAAHMsYmc3a7dCli3QfevppWWauUcPRo1NKKZXBaAA+dQqmTIHvvoPNm6Uk5KJFULo05M7t6NEppZTKoDLnHrAx8Mcf8PjjULIkjB4NJUrApUtyf+3aGnyVUkqlqcw5A960SRKo8uSBIUNg0CCZ8SqllFLpJHME4P37Jakqa1apUFW3LixeDI89Bj4+jh6dUkqpTCjjLkFHR8Py5dC8OVSqJPu8ISFyn2XJ8rMGX6WUUg6ScQPwm29KLeYDB+DDD+HkSZg0ydGjUkoppYCMFIA3bYJnnoGNG+XjZ5+FZcukNvObb0KBAo4dn1JKKXUb1w7AYWEwaxbUqQMPPihHiQID5b4KFWQG7J45trmVUkq5FteNTsZIgYwDB2SPd+JE6NkzrmKVUkop5cRcNwBbFrz7Lvj6QuPGWpNZKaWUS3HdAAzQo4ejR6CUUkoli2vvASullFIuSgOwUkop5QAagJVSSikH0ACslFJKOUCyA7BlWcUsy/rDsqx9lmXtsSzr5dQcmFJKKZWRpSQLOhp4zRizzbKs7MBWy7J+N8bsTaWxKaWUUhlWsmfAxpgzxphtN/99DdgH+KbWwJRSSqmMLFX2gC3LKgnUBDamxtdTSimlMroUB2DLsrIBS4FXjDEh8dw/wLKsLZZlbblw4UJKv51SSimVIaQoAFuW5YEE37nGmGXxXWOMmWKM8TPG+OXPnz8l304ppZTKMFKSBW0B04F9xphPU29ISimlVMaXkhlwA6An0NSyrB03/7RJpXEppZRSGVqyjyEZY/4GtAWRUkoplQyWMSb9vpllXQCOp+KXzAdcTMWv50j6WJxPRnkcoI/FWWWUx5JRHgek/mMpYYyJNwEqXQNwarMsa4sxxs/R40gN+licT0Z5HKCPxVlllMeSUR4HpO9j0VrQSimllANoAFZKKaUcwNUD8BRHDyAV6WNxPhnlcYA+FmeVUR5LRnkckI6PxaX3gJVSSilX5eozYKWUUsoluUQAtiyrlWVZgZZlHbIsa3g891uWZX1x8/5/Lcuq5YhxJiYpPZQty2psWdbV24qbvOeIsSaFZVnHLMvadXOcW+K53+mfF8uyKtz2f73DsqwQy7Jeuesap31OLMuaYVnWecuydt92Wx7Lsn63LOvgzb9zJ/C59/y9Sm8JPJbxlmXtv/nz851lWbkS+Nx7/iymtwQey0jLsoISK1zkTM9LAo9j4W2P4ZhlWTsS+Fxne07iff116O+LMcap/wA24DBQGvAEdgKV77qmDfALUhjkIWCjo8edwGMpDNS6+e/swIF4Hktj4EdHjzWJj+cYkO8e97vE83LbeG3AWeTcnks8J8AjQC1g9223fQwMv/nv4cBHCTzWe/5eOcljeRRwv/nvj+J7LDfvu+fPopM8lpHA64l8nlM9L/E9jrvu/x/wnos8J/G+/jry98UVZsB1gUPGmCPGmEhgAdDhrms6ALON+AfIZVlW4fQeaGJM5uuh7BLPy22aAYeNMalZLCZNGWPWAZfvurkD8M3Nf38DdIznU5Pye5Wu4nssxpjfjDHRNz/8Byia7gNLhgSel6RwquflXo/jZj+ArsD8dB1UMt3j9ddhvy+uEIB9gZO3fXyK/watpFzjVKx791CuZ1nWTsuyfrEs64H0Hdl9McBvlmVttSxrQDz3u9rz0o2EX0xc5TkBKGiMOQPyogMUiOcaV3tuAPoiKyrxSexn0Vm8cHM5fUYCS52u9Lw0BM4ZYw4mcL/TPid3vf467PfFFQJwfPWm707dTso1TsO6dw/lbcgSaHXgS2B5Og/vfjQwxtQCWgPPW5b1yF33u8zzYlmWJ9AeWBzP3a70nCSVyzw3AJZlvQ1EA3MTuCSxn0Vn8DVQBqgBnEGWb+/mSs9Ld+49+3XK5ySR198EPy2e21L8vLhCAD4FFLvt46LA6WRc4xSsRHooG2NCjDHXb/77Z8DDsqx86TzMJDHGnL7593ngO2SZ5nYu87wgLxLbjDHn7r7DlZ6Tm87FLvXf/Pt8PNe4zHNjWVYvoC3wtLm5IXe3JPwsOpwx5pwxJsYYYwemEv8YXeJ5sSzLHegMLEzoGmd8ThJ4/XXY74srBODNQDnLskrdnKV0A36465ofgGduZt0+BFyNXVJwJjf3TO7ZQ9myrEI3r8OyrLrIc3Qp/UaZNJZlZbUsK3vsv5Fkmd13XeYSz8tNCb6bd5Xn5DY/AL1u/rsX8H081yTl98rhLMtqBQwD2htjQhO4Jik/iw53V/5DJ+Ifo0s8L0BzYL8x5lR8dzrjc3KP11/H/b44OjMtKX+QbNoDSBba2zdvew547ua/LWDizft3AX6OHnMCj+NhZNniX2DHzT9t7nosLwB7kCy7f4D6jh53Ao+l9M0x7rw5Xld+XrIgATXnbbe5xHOCvGk4A0Qh79L7AXmB1cDBm3/nuXltEeDn2z73P79XTvhYDiF7b7G/L5PvfiwJ/Sw64WP59ubvwb/Ii3dhZ39e4nscN2+fFfv7cdu1zv6cJPT667DfF62EpZRSSjmAKyxBK6WUUhmOBmCllFLKATQAK6WUUg6gAVgppZRyAA3ASimllANoAFZKKaUcQAOwUkop5QAagJVSSikH+D9McIxaRAL+LwAAAABJRU5ErkJggg==\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.6rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
