{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prediction (out of sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:20.077942Z",
     "iopub.status.busy": "2022-02-08T18:18:20.077485Z",
     "iopub.status.idle": "2022-02-08T18:18:21.113748Z",
     "shell.execute_reply": "2022-02-08T18:18:21.114527Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:21.117971Z",
     "iopub.status.busy": "2022-02-08T18:18:21.117005Z",
     "iopub.status.idle": "2022-02-08T18:18:22.512680Z",
     "shell.execute_reply": "2022-02-08T18:18:22.511380Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "plt.rc(\"figure\", figsize=(16, 8))\n",
    "plt.rc(\"font\", size=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Artificial data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:22.516701Z",
     "iopub.status.busy": "2022-02-08T18:18:22.515728Z",
     "iopub.status.idle": "2022-02-08T18:18:22.522247Z",
     "shell.execute_reply": "2022-02-08T18:18:22.522927Z"
    }
   },
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "sig = 0.25\n",
    "x1 = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x1, np.sin(x1), (x1 - 5) ** 2))\n",
    "X = sm.add_constant(X)\n",
    "beta = [5.0, 0.5, 0.5, -0.02]\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + sig * np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:22.526179Z",
     "iopub.status.busy": "2022-02-08T18:18:22.525179Z",
     "iopub.status.idle": "2022-02-08T18:18:22.548854Z",
     "shell.execute_reply": "2022-02-08T18:18:22.549583Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.986\n",
      "Model:                            OLS   Adj. R-squared:                  0.985\n",
      "Method:                 Least Squares   F-statistic:                     1056.\n",
      "Date:                Tue, 08 Feb 2022   Prob (F-statistic):           2.09e-42\n",
      "Time:                        18:18:22   Log-Likelihood:                 3.5963\n",
      "No. Observations:                  50   AIC:                            0.8074\n",
      "Df Residuals:                      46   BIC:                             8.455\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.9445      0.080     61.793      0.000       4.783       5.106\n",
      "x1             0.5051      0.012     40.926      0.000       0.480       0.530\n",
      "x2             0.5786      0.049     11.927      0.000       0.481       0.676\n",
      "x3            -0.0201      0.001    -18.548      0.000      -0.022      -0.018\n",
      "==============================================================================\n",
      "Omnibus:                        5.168   Durbin-Watson:                   2.286\n",
      "Prob(Omnibus):                  0.075   Jarque-Bera (JB):                4.064\n",
      "Skew:                          -0.646   Prob(JB):                        0.131\n",
      "Kurtosis:                       3.530   Cond. No.                         221.\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "olsmod = sm.OLS(y, X)\n",
    "olsres = olsmod.fit()\n",
    "print(olsres.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In-sample prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:22.553872Z",
     "iopub.status.busy": "2022-02-08T18:18:22.552344Z",
     "iopub.status.idle": "2022-02-08T18:18:22.560837Z",
     "shell.execute_reply": "2022-02-08T18:18:22.561528Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.44209751  4.95658678  5.4266468   5.82074414  6.1187256   6.31512937\n",
      "  6.42008236  6.45763636  6.46181645  6.47103071  6.52175996  6.64256397\n",
      "  6.84938926  7.14294909  7.50860625  7.91877796  8.33746763  8.72617897\n",
      "  9.05024085  9.28450366  9.41747131  9.45319003  9.41058356  9.32034405\n",
      "  9.2198894   9.14721566  9.13465432  9.20356039  9.36080379  9.59764068\n",
      "  9.89115002 10.20799947 10.50992173 10.76000182 10.92874206 10.99890851\n",
      " 10.96836237 10.85041146 10.67162376 10.46746233 10.27645888 10.1338831\n",
      " 10.06594792 10.08550286 10.18992285 10.3615395  10.57054403 10.77988471\n",
      " 10.9513544  11.05186784]\n"
     ]
    }
   ],
   "source": [
    "ypred = olsres.predict(X)\n",
    "print(ypred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a new sample of explanatory variables Xnew, predict and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:22.565511Z",
     "iopub.status.busy": "2022-02-08T18:18:22.564079Z",
     "iopub.status.idle": "2022-02-08T18:18:22.573539Z",
     "shell.execute_reply": "2022-02-08T18:18:22.574200Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[11.04644133 10.88975255 10.60449217 10.24236938  9.87145168  9.55949953\n",
      "  9.35737629  9.28659508  9.33405158  9.45523248]\n"
     ]
    }
   ],
   "source": [
    "x1n = np.linspace(20.5, 25, 10)\n",
    "Xnew = np.column_stack((x1n, np.sin(x1n), (x1n - 5) ** 2))\n",
    "Xnew = sm.add_constant(Xnew)\n",
    "ynewpred = olsres.predict(Xnew)  # predict out of sample\n",
    "print(ynewpred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:22.578250Z",
     "iopub.status.busy": "2022-02-08T18:18:22.576737Z",
     "iopub.status.idle": "2022-02-08T18:18:22.836517Z",
     "shell.execute_reply": "2022-02-08T18:18:22.837233Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd1f43624c0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6MAAAHWCAYAAACCFl9HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB1S0lEQVR4nO3dd3hURRvG4d9JCBBq6CUU6YJ0goqICoJBEEGK2EUpfgrSJEKwoYKAQRQUC4gUlaKAEQFBaYKFEppIL6KS0CHUACGZ748hSDChZbMn5bmva6+Ys7tn38Rls8/OzDuOMQYRERERERERb/JxuwARERERERHJfBRGRURERERExOsURkVERERERMTrFEZFRERERETE6xRGRURERERExOsURkVERERERMTrsrhdQMGCBc0NN9zgdhkiIiIiIiKSClatWnXQGFPo0uOuh9EbbriBiIgIt8sQERERERGRVOA4zl9JHdc0XREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOte76V7JsWPH2L9/P7GxsW6XIhmIn58fhQsXJk+ePG6XIiIiIiKSKaXpMHrs2DH27dtHYGAg/v7+OI7jdkmSARhjiImJITIyEkCBVERERETEBWl6mu7+/fsJDAwkR44cCqLiMY7jkCNHDgIDA9m/f7/b5YiIiIiIZEppOozGxsbi7+/vdhmSQfn7+2v6t4iIiIiIS9J0GAU0IiqpRs8tERERERH3pPkwKiIiIiIiIhmPwqiIiIiIiIh4ncJoKujQoQOO4+A4zoUtRBo2bMioUaOuaY3i4sWLcRyHgwcPpmK1IiIiIiIi3qcwmkoaN27Mnj172LVrFz/88AMtWrTgtddeo0GDBpw8edLt8kRERERERFyVKcJo+JpI6g9ZSJl+s6k/ZCHhayJT/TGzZctG0aJFCQwMpGbNmvTu3ZvFixezevVq3n77bQC++OIL6tatS+7cuSlcuDDt2rW7sPflrl27aNiwIQCFChXCcRw6dOgAwNy5c2nQoAH58uUjf/78BAcHs2nTplT/mURERERERDwlw4fR8DWRhM5YT2R0DAaIjI4hdMZ6rwTSS1WtWpWmTZsyffp0AM6ePcvrr7/OunXrmDVrFgcPHuThhx8GoGTJkhdut2HDBvbs2cOIESMAOHnyJD179mTFihUsXryYvHnz0qJFC86ePev1n0lEREQks3FjoEMkI8ridgGpLWzeFmJi4xIdi4mNI2zeFlrVCvR6PVWqVGH+/PkAPP300xeOly1blo8++ojKlSuze/duSpQoQf78+QEoXLgwBQsWvHDbNm3aJDrnuHHjyJMnDytWrOD222/3wk8hIiIikjklDHQkvL9MGOgAXHlvKZKeZfiR0ajomGs6ntqMMRf2t1y9ejUtW7akdOnS5M6dm6CgIAD+/vvvy55jx44dPPLII5QrV448efJQpEgR4uPjr3g/EREREUmZyw10iMi1yfBhtHiA/zUdT20bN26kbNmynDx5kuDgYHLkyMHnn3/OypUrmTt3LsAVp9ved999HDhwgE8++YTly5ezZs0asmTJomm6IiIiIqksrQ10iKRnGT6MhgRXwt/PN9Exfz9fQoIreb2WP/74g7lz59K2bVs2b97MwYMHeeutt7jjjju48cYb2b9/f6LbZ82aFYC4uH8/fTt06BCbN2+mf//+NG7cmMqVK3P8+HHOnTvn1Z9FREREJDNKawMdIulZhg+jrWoFMrh1NQID/HGAwAB/Breulupz+s+cOcPevXuJiopi3bp1DB8+nLvuuos6derQp08fSpUqRbZs2fjggw/YuXMns2fP5pVXXkl0jtKlS+M4DrNnz+bAgQOcOHGCfPnyUbBgQcaMGcP27dv56aef+N///keWLBl++a+IiIiI69LSQIdIepcpEkyrWoFeX1A+f/58ihUrhq+vLwEBAVStWpUBAwbQpUsXsmbNSs6cOZkwYQL9+/dn1KhRVK9eneHDh9O0adML5wgMDOT111/npZdeolOnTjzxxBOMHz+eqVOn0r17d6pWrUr58uV55513/tPUSEREREQ8L+E9Zdi8LURFx1A8wJ+Q4EpqXiRyHRxjjKsFBAUFmYiIiCSv27RpE5UrV/ZyRZKZ6DkmIpI2hK+J1Jt7EZEMynGcVcaYoEuPZ4qRUREREUm7tFWGiEjmlOHXjIqIiEjapq0yREQyJ4VRERERcZW2yhARyZw0TVdERESS5K11nMUD/IlMInhqqwwRkYxNI6MiIiLyHwnrOCOjYzD8u44zfE2kxx9LW2WIiGROCqMiIiLyH95cx+nWnuAiIuIuTdMVERGR//D2Ok439gQXERF3XVUYdRznDqAPUAcoDjxljBl/0fWtgWeA2kBBoKExZrGnixURERHv0DpOkcs4cQImT4Y9e8DP78qXLFmSPl6mDBQo4PZPI+Kaqx0ZzQX8AUw8f7lUTuBX4ItkrhcREZF0JCS4UqK9P+Hq1nF6q+mRiCt27YIPPoBPP4WjR1N+vhw5oE8fCAmBXLlSfj6RdOaqwqgxZg4wB8BxnPFJXP/5+esKerI4ERERcUdCgLyWYJnQ9CghwCY0Pbr4fCLpjjHwyy/w3nvwzTfgONC2LfToATffDLGxV3c5dy7x92fPwpdfwhtvwOjR8Oab8NRT4Ot7xZJEMgqtGfUwx3Eue/2TTz7J+PHjvVOMiIhIClzrOs6weVs4c+YstaO2UvnAn0TlKcTfeYsyYpZRGJX05+xZ+PprG0IjIiBfPjuC2bUrlCz57+18fSF79ut7jPvvh5494YUXoHNnGDkShg2De+7xxE8gkua5EkYdx+kCdAEoVaqUGyWkmj179lz471mzZtG5c+dEx/z9E6+1iY2Nxc/Pz2v1iYiIeNzJk/Djj/SYNIpGO1ZS8FQS0xc/Lg5lyyZ9KVrUjjaJpAUHD8Inn8CoUXZNaKVK8NFH8PjjkDOn5x+vXj078jptGvTtC8HB0LQphIVB1aqefzyRNMSVMGqMGQ2MBggKCjJu1JBaihYteuG/AwICEh3btWsXxYoVY9KkSYwZM4bffvuNsLAwcuXKRbdu3Thx4sSF+y5evJiGDRty4MABCha0s59//fVXQkNDWblyJfny5eP+++9n6NCh5MmTx3s/oIiIpEseX8u5dy989x3MnAnz58Pp09ybPScLywQxv/zNrA6sTJEThygZvZeqpw/RuQSwcycsXAiff26nPibw97eNXCpUgE6doHlzhVPxvj/+gBEj4Isv4PRpGwo/+8yOUvr8dzfEc+cgJgZy5/bAYzsOtGtnR0pHjbJTdmvUgI4d7TTei95fimQk2mfUBaGhoTz33HNs3LiRVq1aXdV91q9fzz333MP999/PunXrmDFjBmvXruXpp59O3WJFRCTdS1jLGRkdg+HftZzhayKv/iTGwIYNMHgw3HorFCsGXbrA+vX26/z5LPxpPf3a9OO7KncSmbcwqwMr80PNxhR6eyCMGwc//QT//AOnTsHmzTBnjm0G87//2SC6ahW0aGFHin78MXFgFUktO3fakchq1WwQfeIJ+1yfO9cevyiIGgPLl9vloiVK2JzoUdmyQe/esH07PP+8/XdToQIMHGj/3YhkMOluzWjPnrB2rXcfs2ZNu1zAU55//nnatm17TfcJCwujffv2vPDCCxeOffTRR9SqVYv9+/dTuHBhzxUoIiIZSti8LYm64gLExMYRNm/LlUdHN260nUO//da+aQeoW9eO3Nx/v30Df34UsyVg/LJeeQQ2e3Y79bHSJZ15Y2NhwgT7Dv+ee+COO+zj3HFHCn56kcuYNs2OPvr4wFtv2Q9WkthqZfNmmDTJXnbsgKxZ4b77oHHjVKqrQAH75rNrV+jXD155BT7+GAYNstOFkxipFUmPrnaf0VxA+fPf+gClHMepCRw2xvztOE5+oBQQcP425R3HiQb2GmP2erTiDCAoKOia77Nq1Sq2b9/O1KlTLxwz5z8x3rFjh8KoiIgkKyqJ/UIvdxywQ0AffmgbqwDcfTe8+KJ9Bx6YfIBtdlMguQ4F8vvv9v109n2wZo2dZViokN1uMVl+fnaa7uOPw5gx9o33nXdCkyY2lN5yy1X8tCJX4fRpu6XKqFG2I+7UqXDDDYluEhkJU6bYALp6tf3MpVEj6N8fWreG86uxUleFCjB9Oixdav8tduhgQ+o779hiRNK5qx0ZDQIWXfT96+cvE4AOwP3AuIuuH3PR7QakqMJLeHKE0i05L1n87uPjcyFYJoiNjU30fXx8PJ06daJXr17/OV/gZd4UiIiIFA/wJzKJ4Fk8wD+JWwNHjtjRom++gXvvhfHj4TIfekZF2Rm3c+bY2bUXtUBIxHFsQC1SxIbTIkX+vRQtCsWL20HQbNmyQbdu8PTTtnHMkCF2avB999lR01q1ruO3IHLe9u3w4IP2U5Leve3U86xZAfvUnz7dBtDFi+1nMkFB8O679i7Fi7tUc4MGsGyZDc2hofbDofvus02ObrzRpaJEUu5q9xldDCTbScAYMx4Y75GKMqFChQpx6tQpjh07dqEZ0dpL5iLXrl2bDRs2UL58+STOICIikryQ4EqJ9v8E8PfzJSS40n9v/Ouv8PDDtovoO+/Y9TGXTAmMi4OVK2H2bHtZs8YeL1ECHn3U9h+qVw+OHrV9jvbt+/dy8ffLltnvL14KV6iQ3eHi2WehRIkcdjSoSxd4/337xrt2bbvH4+uvQ5UqqfDbkgztq6/s6HuWLHbq+f33Ex8P4TNsX605c+yOLhUqwKuvwiOPQMWKbhd9no+P/bf5wAO20dJbb9m1ZOHhdm2rSHpkjHH1UqdOHZOcjRs3JntdevD1118b+yu2/vzzTwOYlStXJrrdoUOHTM6cOc1zzz1ntm3bZqZNm2ZuuOEGA5gDBw4YY4xZt26d8ff3N88884xZvXq12bZtm/nuu+9Mly5dvPozZTTp/TkmInK1vlm929w2eIG5oe8sc9vgBeab1bsT3yAuzpi33jLG19eYsmWNWbEi0dWHDxszebIxjz9uTMGCxoAxPj7G3H67MYMHG/P778bEx19fbcePG7N9uzGzZhnTsqUxjmPLaNPGmMWLLzrvkSPGvPqqMblz2xs9+qgxW7de34NK5hITY8yzz9on7q23GrNrlzHGmD/+MKZ+fXu4aFFjevY0ZuXK638ue9W+fcbUrGlMtmzG/Pij29WIXBYQYZLIggqjqehqw6gxxoSHh5sKFSqY7Nmzm3vuucd8/vnnicKoMcasXLnSBAcHm9y5c5scOXKYqlWrmldeecUrP0tGld6fYyIiHrFnjzGNG9u3Be3bGxMdbYyxb8g//9yYBg1sOARjChSwGXDSJGMOHUqdcv7805iQEGPy5bOPWa2aMZ98YsyJE+dvcPCgMX37GpMjhy0sJMSY2NjUKUbSv61bbWgD+1w5e9acOmVMaKgxWbLY5/Rnnxlz7pzbhV6HgweNqV7dmOzZjVmwwO1qRJKVXBh1jMtt04OCgkxERESS123atInKlSt7uSLJTPQcE5FM74cfbMOg48dh5Ei7VtRxiIyEZ56x03CrVLEzA5s3t71efH29U9qpUzB5sp2hu26dbRjz9NPw3HNQrhx2ru/LL9tuv3feadfTFSnineIkVXlsX9wpU+y876xZYeJEaN6cefPsc2jnTtsPKCwMzm/pnj4dOGCbGe3YAd9/b/8tiKQxjuOsMsb8p4ur+kKLiIhkRrGxdsuI4GC7UHPlSujUCYPDhAlQtSosXGgbt6xfb7c5rFfPe0EUIEcOm43XrLHNRIODbV6uUMH2bpm7pgjxn4yxIWPFCrue9NdfPV5H+JpI6g9ZSJl+s6k/ZOG17c8q18wj++LGxNhPUx5+GKpXh7Vr2VunOY88YpdXZslin9/jxqXzIAr23++CBVCmjP3EaOlStysSuWoKoyIiIpnNrl22be3QofYN+4oVcNNNREVBixZ2tOimm+xoZBL9i7zOceD22+0g165ddsvFiAjb6PfGG2GS7+OY35aBv78dFXr/fdsG1QM8EozkmlxuX9yrsmWL7b48ejT07Uv8wsWM/r4klSvbTrkDBsDvv0PDhp6v3TWFC9tAWqIENGuWKh/KiKQGhVEREZHMZPp024Fz40bbWfTjjzH+OZg40QbQhNHQn36yI5BpTWCgbaT711/w5ZeQO7ft4Hv/y9WJmnk+oXbvDo89BidPpvjxUhyM5Jpd1764CWbOtHuxREbC7Nn88dgQGjTy45ln7I5Av/8Or70G2bJ5uOi0oGhR+w+4WDE7/Lt8udsViVyRwqiIiEhmEB9vQ1rbtlCpEqxdC+3aERUF998PTz6ZeDTUm9Nxr0e2bHbbjRUrbHhesAAq1wtgbItwzMBBdrHprbfC1q0pepwUBSO5Lsntf5vsvrgJxo61i5srVybmt7WELm1GrVp2oHTCBPscqZTEbkYZSvHisGiRHSm95x47/V4kDVMYFRERyeiMsft1vv++TZo//4y5ocyF0dD582H48LQ7Gno5vr72R/r9dzvy1amLD8E/9WffxHl2r9S6de0+jNfpuoORXLeQ4Er4+yX+NCTZfXHBPr+HDLH7hzZpwo+hC7kpuARDhtjeXJs3wxNP2OnemUJgoA2kBQrYQLpqldsViSRLYVRERCSjGzwY3nsPevSA4cOJOuCXaDT099+hV6+0Pxp6OeXL2xmKH34Iv/0G5Z9twsSeqzEVK9rRsn794Ny5az7vNQcjSbFWtQIZ3LoagQH+OEBggD+DW1dLuptufDz07g2hocQ9+DBdS87knta5yJYNFi+Gzz7LAA2KrkfJkjaQ5s0LTZrYmRAiaZC2dpFMTc8xEUmPrmnbi9GjbZOixx7DjJ/AF5N86N4dTp+Gt96yM3fTcwhNyl9/QZcudteaxrefZlpgD/JOHW23v5g82U5hvAYe22ZEPCs2Fp56Cr78kmMdunP37+8SsdqHvn3tuuIMuS70Wv35p23qdeqU/bSmenW3K5JMKrmtXRRGJVPTc0xE0puE7q4XN9Xx9/NNeuRo+nR48EFo2pSYyeE89pQfM2bAbbfZLS0qVvRy8V5kjP0Ze/eGs2fhm/vHcc+3z+EULAjTpsEtt7hdoqTEyZN2/fPcuWx5YhD1ZoYSbxwmTrRroOUiO3bYQHrmjB0trVrV7YokE9I+oyIiIhnAVXd3XbjQdvi59VZOTfiaFq39+OYbu5vLkiUZI4hebv9Px4Gnn4YNG+Duu6Hp1KfoUOFXzuIHDRrYEWNJnw4dgrvvxvzwA+HNx3DjxP6UvsFh1SoF0SSVK2dDqJ+fnR2wcaPbFYlcoDAqKdatWzfuuuuuC9936NCB++67L0XnHDBgAFX1yZ2IyH9cVXfXVaugZUuoWJETU2bRrG0OFi2C8ePhxRczxrTcq93/MzDQ7vbx5ZcwK7IWpfZFsKNMYzt1OTTUY/uRipf88w80aIBZu5bXq07jgdmd6NjRbqtZrpzbxaVhFSrYQOrjYwPp5s1uVyQCKIymmsjISLp06UKJEiXImjUrgYGBdO7cmd27dye63ZWC27p162jZsiVFixYle/bslCpVijZt2vDXX3+l9o9w3UaMGMEXX3xxVbfdtWsXjuNw6VTtPn368NNPP6VGeSIi6doVu7tu3Wr32ixQgONfz6Xpw/n4+Wf4/HPbUTSjuJb9Px3HDhJv3Ajl7spOpa0zGZOjAwwZwt8t2tl5vJL2bdoEt93Gub8jaZd7HkO3PsDYsfDpp+Cv5sZXVqmSnTFhjA2k27a5XZGIwmhq+PPPPwkKCuKPP/5gwoQJbN++nS+++IINGzZQt25ddu3adVXnOXDgAHfffTe5cuVi9uzZbN68mc8//5xy5cpx7Ngxj9Z81oN/iPPmzUtAQECKzpErVy4KFCjgmYJERDKQy3Z3jYy0nTOBY9N/5J6nAlm+3PbseeQRN6pNPdez/+dvUZEcvmUJ+Vqu45lzn/CK7wBKzZ7O/jubwPHjqVWqeMLy5Zjbb+fk0VhuPf0Ta/Lcya+/2qnYcg2qVLEbrsbGQsOGsH272xVJJqcwmgq6du2Kj48P8+fP5+6776ZUqVI0bNiQ+fPn4+PjQ9euXa/qPL/88gtHjhxh3Lhx1KlThxtuuIE777yTt99+m2rVqiV7v4TR1oEDB1KkSBFy5crFU089RUzMv3+g77rrLp599ln69OlDoUKFqF+/PgAbN26kefPm5M6dm8KFC/Pwww+zd+/eC/eLi4ujT58+5MuXj3z58tGzZ0/i4uKSfPwExhjeeecdKlSoQLZs2ShRogShoaEAlClTBoC6deviOM6F6b6XTtONj4/nzTffpGTJkmTLlo1q1arx7bffXrg+YYR1+vTpNGnShBw5clClShV+/PHHq/pdi4ikF8lue1HaH4KD4cgRjk6dy93/q8CqVfDVV9CundtVe9717P+ZMJqa88a9FHvqZ94u1IOn+Iz8y37B3HEnXPT3zi2XWwebac2di2nUiH1n81H9+C8ENq/JqlV2X1m5DlWr2kB6+rQNpP/843ZFkokpjHrY4cOHmTt3Ll27diVHjhyJrsuRIwfPPfcc33//PUeOHLniuYoWLUp8fDzTpk3jWrse//TTT6xbt44FCxYwffp0fvjhB/r27ZvoNl988QXGGJYuXcrEiRPZs2cPd9xxB1WrVmXFihXMnz+fEydO0LJlS+Lj4wF45513GDNmDJ988gm//fYbcXFxfPnll5etpX///rz55puEhoayYcMGvv76a0qWLAnAihUrAJg7dy579uxhxowZSZ5jxIgRhIWFMXToUNavX88DDzxA69atWXvJvlkvvfQS3bt3Z926ddStW5eHHnqIEydOXNPvTkQkrWtVK5Bf+jXizyHN+aVfI1pVDID77oNt2zg68Vvu6l2b33+HGTPsFpsZ0fXs/3nxqKlfQAxFH/2N6bUb0oLviPl9K+durgdb/jvN11uudh1spjJpEqZFCzbFVaT2yZ/pMqQc33wDKZyAJdWrw/z5cPSo7Up85ozbFUkmlcXtAq5Zz57e37i3Zk27WfhV2LZtG8aYZLcLqVKlCsYYtm3bxs0333zZc916663079+fJ598kq5du1K3bl3uuusuHn30UUqXLn3Z+/r6+jJu3Dhy5cpF1apVGTp0KB07dmTw4MHkzJkTsKOS77zzzoX7vPrqq9SoUYOhQ4deODZx4kTy589PREQEN998M++99x4vvvgiDz74IGBD4rx585Kt48SJE7z77ru89957PH1+Lk358uWpV68eAIUKFQKgQIECFC1aNNnzDBs2jD59+vDI+Xlmb7zxBkuWLGHYsGGJ1qf26tWLFi1aAPDWW28xceJE1q5dy+23337Z35eISLoVG2uHPpcv5+inX3P7Kw3Zvt027QkOdru41JOwjc217P9ZPMCfyIsCqZMlnvxNNrCzUmnu/XYxX+9uTkDd28g6bxac/zvlTZdbB5sp9zUdMQJ69mSpz108HRDOl1/lpWFDzz9Mpt1HtmZNu/9R27bQqxd8+KHbFUkmpJHRNG7QoEHs3buX0aNHU61aNcaOHUuVKlVYsGDBZe9XvXp1cuXKdeH7evXqcfbsWXbs2HHhWJ06dRLdZ9WqVSxZsoRcuXJduCSMYO7YsYOjR4+yZ8+eC0ESwMfHh1sus1fbxo0bOXPmDHffffc1/dwXO3bsGFFRURemEie4/fbb2XhJe/LqF23mXLx4cQD2799/3Y8tIpKmxcdDhw7w/fdED/2YemGt2bEDZs3K2EE0wX9GiK8QIJIbTR0ako8xa4J4quKv/HU8P7F3NCJuxrfJnCX1XM862AzJGOL69oeePZnBAwy45XuWrEu9IJqpR6PbtIEXXoCPPrJdzkS8LP2NjF7lCKVbypcvj+M4bNy4kQeSmBu1ceNGHMehfPnyV33OAgUK0K5dO9q1a8fgwYOpVasWb775ZooCHnBhhDRBfHw8zZs3Z9iwYf+5bZEiRS5M1U1LHMdJ9L2fn99/rkuLdYuIpJgxdjRj0iSO9nuLW8Z0JjISvv/e7m8v/3Wl0dRpa8rRv/OvPPTlfQS1aU304FEE9Puf1+q7dOT24uOZxrlznHj8f+SaMpZP6MK2nh8y721fLvrz7lEajQaGDIGVK+12RzVq2Cm8Il6ikVEPK1CgAMHBwXz44YecOnUq0XWnTp1i1KhR3HvvveTPn/+6zp81a1bKlSt3xXWQ69ev5+TJkxe+X7Zs2YX7Jqd27dps2LCB0qVLU758+USX3LlzkzdvXooVK8ayZcsu3McYc2HdZ1IqV65MtmzZkh3JzZo1K8B/miBdLE+ePBQvXpxffvkl0fGff/6ZKlWqJHs/EZEM7a23YORIjnXsRZ2v+rFnD8ybpyB6JZcbTfX3h3e/KMTOMQuZ59OMgNBn+fPRl722F+n1rIPNUGJi2H9nW3JNGcvQrK9Q4KuPGfZu6gVR0Gg0AFmywNSpkDevHSk9etTtiiQTURhNBR988AHnzp2jcePGLFy4kH/++YfFixfTpEkTjDF88MEHiW5/7Ngx1q5dm+iya9cuZs2axWOPPcasWbPYunUrW7ZsYdiwYcyZMyfJUdeLnTt3jqeffpoNGzbw448/0q9fPzp37vyf0dCLde3alaNHj9K+fXuWL1/Ozp07mT9/Pl26dOH4+Zb3PXr04O2332batGls2bKFnj17smfPnmTPmTt3bnr06EFoaCjjxo1jx44drFixgo8++giAwoUL4+/vz7x589i3bx9Hk3kBDAkJYdiwYUyePJmtW7fy6quvsnTpUvr06XPZ34OISIb01Vfw8sucaPUYNecP4+Ahhx9+gEtWM8h1erhTTm5Y8w1fB3SmzKRBrKn9NHGnY1P9cZPtlJwJRujiD0fzd+VgCv46k0FF36fV72/Qtp1z5Tte5Ho6EV9PV+YMqWhR+7ry55926r+XPoARSX/TdNOBcuXKERERwRtvvMHjjz/O/v37KVSoEM2aNWPq1KmUKFEi0e2XLl1KrUv6k7dp04a3336bXLly0adPH/755x+yZMlCmTJlGDZsGD169LhsDXfeeSc33XQTDRs25NSpUxfOdzkJo4+hoaE0bdqU06dPU6pUKe655x6yZcsGwAsvvMDevXvp1KkTAI8//jiPPvoomzZtSva8gwcPJl++fLz55pvs3r2bIkWK8MT5ndezZMnCyJEjeeONN3j99ddp0KABixcv/s85unfvzvHjx3nxxRfZt28flSpVYvr06dSoUeOyP5OISIazdSt07EhMrXrUiBhL9Ekf5s+HoCC3C8tYqlTPQul/PmH67SVos/Y1VgbuofTKaRQum+vKd06BVrUCM0X4vFj0pj0cuaUpgcc3Mar+ZHrOa89lPjtPUsLaz4QptwlrP4HL/j5Dgisluh9kstHoizVoAGFh0Ls3vP02XLILg0hqcK51yxBPCwoKMhEREUlet2nTpmS70kryOnTowMGDB5k1a5bbpaR5eo6JSLpy6hTceitxu6O4Nesa/jxXkvnzbVNMSR3GwM9PjaXehGfYmKUG0Z9/xx0PFXe7rAxjQ/g28rS7h3znDrCw2ze0GNkE59oGRAGoP2RhkuttAwP8+aVfo8veN9N2002KMdC+PUyfbrd+SY2uUZIpOY6zyhjzn49NNTIqIiKSXjz/PKxfz3Olvmfb0ZIsXQrVqrldVMbmONBgfEf+rFuMct3bE/1wXd797lue+yyI85OG5DqFv7KK2wbei6+PYde4xdzf4fqH91Oy9jMzjkYny3Fg7FhYvx4eeghWr4ZA/W4k9WjNqIiISHowfjx89hmTy77EZ1FNmT5dQdSbynRths+vv5A9tx/PTLqDlyt+xYYNbleVPsXEQFjTBdw98C7is+fA+fUXqqYgiILWfnpU7twwYwacPGn3MD571u2KJANTGM2Axo8frym6IiIZyfr1mOeeY0vxhjy283U+/BBSuLuXXAf/W6pTYPsKTleuTdjf7ZlR83VGfWDU6+UKLm4sVCfkN/qUn0z3ec04VfgGCm39lfy3VEzxY2T6TsSeVrmyHSH97TcICXG7GsnAFEZFRETSsuPHoV07Tvnl5c6oSbwQ4kvnzm4XlYkVLkz+NQuIad+BV84NoODzD9Hm3lPs2+d2YWlTQmOhyOgYTm4vzK3vRfB+1KPsKVuLIpuX4FvSM+tvM3Mn4lTTvj306AEjR8LkyW5XIxlUml8zaozBuZ6V7CJX4HbzLhGRKzIGunTBbNtGC7OAeq2KMmSI20UJ2bLhP/kzTJ2beLDvi5T7YSf33BTO4ImBNGvmdnFpS9i8LZyKMUQvrUivFV8ygNf5oXQ9Bnd4hUX58nn0sbT2MxWEhUFEBHTqBNWrw003uV2RZDBpemTUz8+PmJhMtOmweFVMTAx+qbmTtohISn30EUyZwgDfgRytdRdffAE+afovdybiODghfXBmzqSW/2Z+PFqXV5tH0K2bXRMp1s71/hz87GbeXjGMAbzO1Kr38OyD/dh1Mt7t0uRq+PnZ/Udz5YLWreHYMbcrkgwmTf9JK1y4MJGRkZw6dUqjWOIxxhhOnTpFZGQkhQsXdrscEZGkRURgevViQfZmfFaoL999xzXvvShecN99+C77lYKB2fjNtwEHRk0lKAjWrXO7MHdFR0OXLuA/qQhLjjXhGUYz6tZ29G32PHE+vmoslJ4ULw5Tp8KOHfD002iRtHhSmp6mmydPHgCioqKIjY11uRrJSPz8/ChSpMiF55iISJpy5AjxbduxnyI85TORmbN8KK6tLdOuatXwWbkCn9atmfrzQwz7axO31H2Vt4b40LNn5hrNNsZuUfn883D7vums93uacz7xdLn/JX6oWA9QY6F06a67YPBgePFFGD4cXnjB7Yokg0jTYRRsIFVgEBGRTMMYzJMdiPs7kgfMEkZNK0CtWm4XJVdUqBDMnw/PPkufca9Tr/hGmrwwnu+/z8FHH0H58m4XmPr++Qe6dYO5M88woVAfHjIfQK2b+eHVEWxYfwonOobiAf6EBFfS2s70qE8fWLYM+vaFunXhjjvcrkgyAMft6a9BQUEmIiLC1RpERETSjGHDICSEHrzHDcN70KuX2wXJNTEGhg/HhITwd5Ea3H7wW3bHl6BI3SjeesOXp+8p5naFHhcXZ5c3h4ZCqdgdLCz4IEUiV0Pv3nY0LWtWt0sUTzl61AbR48dh9WoolvGez5I6HMdZZYz5z4bCmWjiiIiISBr3yy/E9+3HNNpw9pnu9OzpdkFyzRwHXniBZe+Np8DhLazMWov7y3zBvpXF6dS8MG2fPk50tNtFes4ff8Dtt9tpuX3Lfs16v9oUOfUnfPstvPOOgmhGkzcvzJhhGxk9+CBoGZ2kkMKoiIiIB4SviaT+kIWU6Teb+kMWEr4m8tpOcOAAp1u1Z0d8Gb64aywj33fQzmbpV59TJXjgsWGczJ6Db3c8yTdlGnND2d+ZPi43ZcvaHTPSc9fd06fhlVegdm34Z9tpNt/dlZd/fxCfmyrDmjVw//1ulyippWpVGD0afv4Z+vVzuxpJ5xRGRUREUih8TSShM9YTGR2DASKjYwidsf7qA2lcHCcfeAwOHqRf2a8Z/01etPNU+hYVHcO2QqW596n3ebf+I9y7aylr/2pA31tfoN4t8bz4IlSoAGPGwLlzbld79cLXRFK1y2pyB55g4EB4sP7v7Cp+G5UWfGib2ixZAqVLu12mpLZHH4WuXW0zo2nT3K5G0jGFURERkRQKm7eFmNi4RMdiYuMIm7flqu5/8qVB5PzlB/rnep9h82sSEJAKRYpXJWxdcjaLHyNuf4SmT49ifdFyDFk2nNlHbmPlp+soVcpuf3LTTfb9fFreMePoUXj2pSM8fF9ONoypDXE+dL51IB//eivxf+2EmTPtemdNy808hg+HW26x271s2+Z2NZJOKYyKiIikUFR00vMtkzt+sbNzF+I/dABf+jxO27mdKFPG09WJG0KCK+Hv53vh+z/zB9LpsSFEvDkCdu4k6Jk6/FKvD7OmnMDPD9q1g5tvtg150wpj7EDnk0/aPjUfv5WPuFgfit61ms9Kt2P0slfYUrA0Dz37EbRo4Xa54m1Zs8JXX3HhCZye552LaxRGRUREUihhFOxqjycw+w9wqvVjbKESfp9+xG31tUg0o2hVK5DBrasRGOCPAwQG+DO4TXWCXu4OmzfD00/jDH+H5iFVWDdwJhMmwIED0KQJNG4M4eFw5Ig7te/ZA0OGQKVKcOedtpYnnoDijy/h0QbvsHBDGx7/fQ4f39KG9o8MYQ3agi/TKlUKJk6EdeugRw+3q5F0SFu7iIiIpFDCmtGLp+r6+/kyuHW15PdTNIZd1VpQbMOPfNp5BV1H1/BStZJm/PIL/O9/tiVtq1acCRvJJ3NKMnCgDaaOAzVrQsOG0KgRNGgAqbX1emwszJkDY8far3Fx9vE6dYK2950mxzdf8mf/Nyiz/28icxfi5eDnWFSuLmCD9i/9GqVOYZI+hIbaTzA+/xwee8ztaiQNSm5rF4VRERERDwhfE0nYvC1ERcdQPMCfkOBKyQdRYFv396nwfndG3zSCTr93x0dzlTKn2Fi79u7118HXF954g7PPPM/yVVlYtAgWLYLffoMzZ+zVderYcNqwod1SJWfO639oY2DrVvjsM5gwAfbtg6JF7bTcp5+GigUPw8cfw8iRsG8f0ZVuYmDl5oSXv41zvlmAq/jQRTKHc+fsJyarVsHKlVClitsVSRqjMCoiIpJGRH2/jvzNbmFZzsbUjvyOPHk1PTfT+/NP6NbNDkvWqmX3TWnSBHLlIibGBtKEcLp8uX3vnyWL7R+TEE7LlbNTew8fhkOH7NeLL0kdSwi5zZtDx47QrBlk+edPeO89O0x68iQEB0NICDRqRPjaqGv60EUykagoO5RfqBCsWJGyT0okw1EYFRERSQNOHTzF3hJ18D97lFO/rqPcrYXcLknSCmNg+nS79i4qyjaIuesuuO8+eznf3erkSTvDd9EiWLgQIiIgPj750/r7Q/78UKCA/XrxpUQJaNvWNigiIsJ2xP36a/DxgUcesdu1VK/ulR9fMoD58+Gee+xU3QkT0GbJkkBhVERExGXGwILyz9Bo5xhWDvqRW/rf7XZJkhbFxsLPP8Ps2TBrFmw5v0VQlSo2lDZvDrfdZodGgWPHYOlS2Lv335CZEDzz5bNhNFnx8fD99zaELl5sF6U+8wx0726Tqsi1ev11GDAAPv3UDreLoDAqIiLium8em84DX7bltzv7Um/xELfLkfRi2zYbTGfPhp9+smE1Xz5o2tSG06ZNbfJMytmzthvS/v3/fr34snw5bNxog2fPntC5c+p1SZLMIS7OPid//hmWLYMaas4mCqMiIiKuWjThb2p2qMGh/BUot+cXnKx+19z0SIRjx+DHH+2I6Zw5NlD6+NiR0sqV4eDBxMEzOjrp8/j5QeHCdurvM89A+/b2mIgn7N9v14/mymWnf+sDjkwvRWHUcZw7gD5AHaA48JQxZvxF1zvAa0AXIB+wHOhqjNlwpXMrjIqISEa3dVMcB6s1pLpZi8/aNeSoVu76toMRuVh8vH2jP2uWvezebQNmwqVQocTfX3wsb16t55PUtWSJ7bDbujVMnarnWyaXXBjNcpX3zwX8AUw8f7nUi8ALQAdgC/Aq8KPjOJWMMcevq2IREZEM4NgxmHfHIJ6PW8qB4Z9TqFo5AMLmbUkURAFiYuMIm7dFYVSujo8P3HyzvbzxhtvViCR2xx0waBD06wd33gldu7pdkaRBV7WrmTFmjjGmvzFmGpCoX9v5UdGewBBjzHRjzB/Ak0Bu4BEP1ysiIpJuxMfDoHt/5rmDr7O3yWMU6vXvZvBR0TFJ3ie54yIi6U5IiG241auX3X9U5BKe2GK7DFAU+CHhgDEmBlgC3OaB84uIiKRLQ/oe4blfH+V4wTIUnTYq0XXFA5JucZrccRGRdMfHx27xUqwYPPig3QhX5CKeCKNFz3/dd8nxfRddl4jjOF0cx4lwHCfiwIEDHihBREQkbZkx3VB+2DMEOlHknT35Pw08QoIr4e/nm+iYv58vIcGVvFmmiEjqKlAAvvoKIiOhQwe7x5XIeZ4Io9fMGDPaGBNkjAkqVEibfYuISMayfj0seGQsD/I1vDkQ5+a6/7lNq1qBDG5djcAAfxwgMMBfzYvEa8LXRFJ/yELK9JtN/SELCV8T6XZJkpHdcguEhcHMmTB8uNvVSBpytQ2MLmfv+a9FgL8vOl7koutEREQyhcOHoXezzXwb24MzDRqTLTQk2du2qhWo8Cled2kn58joGEJnrAfQ81FST/futsNu375w661Qv77bFUka4ImR0T+xobNJwgHHcbIDDYBfPXB+ERGRdOHcOXis7WmGRT6EX94cZJs60a6ZEklDLtfJWSTVOA589hmULm33tdVSPeEqR0Ydx8kFlD//rQ9QynGcmsBhY8zfjuO8B/R3HGczsBV4GTgBTPJ4xSIiImlUSAjcs6gfNVgHX8yyTTtEUln4mkjC5m0hKjqG4gH+hARXuuwIpzo5i2vy5oWvv4bbboPHH4c5c/SBXSZ3tf/3g4A15y/+wOvn/zthU6u3gXeBUUAEUAy4R3uMiohIZvHBB7DlvTn0ZISdjta8udslSSaQMOU2MjoGw79Tbi+3BlSdnMVVtWvDiBEwbx4MHux2NeKyq91ndLExxkni0uH89cYYM8AYU8wYk90Yc+f5/UZFREQyvO++g+HddzHF73FMjRowdKjbJUkmcT1TbtXJWVzXpQs88gi8+iosWuR2NeIijYuLiIikwKpV8GT708z2b0PuHHE406ZB9uxulyWZxPVMuVUnZ3Gd48Ann0DFitCuHWzReuXMyhPddEVERDKlv/6C++6DD327UfnEarttQfnyV76jiIcUD/AnMongeaUpt+rkLK7LlQtmzbLrR4OD4ddfoXhxt6sSL9PIqIiIyHWIjoZmzaDt0U956MRYePllaNHC7bIkk9GUW0nXypWzTYwOHYKmTe0Lq2QqCqMiIiLX6OxZaNMG8mxZyYhzXe2n+gMGuF2WZEKacivpXp06MGMGbN4MrVrB6dNuVyRe5BhjXC0gKCjIREREuFqDiIjIxS63VYYx8NRTMGvCQXYVqEOuXI5dOFqggMtVi4ikY5Mn26ZGbdrA1Kng63vl+0i64TjOKmNM0KXHNTIqIiJykSttlfHmm/D5hDiWl32YXCf2wfTpCqIiIin18MMwfLh9Te3Rw37yJxmeGhiJiIhc5HJbZRxbH8hrr8HMaq9Sbv18GDvWTjETEZGU69ULoqJg2DAoVgxeesntiiSVKYyKiIhcJLktMXasy0GnV+Hlat/SYv1b0LkzPP20l6sTEcnghg6FvXttU7iiRaFjR7crklSkMCoiInKRpLbKOHswFwfD69Ck1Dbe2PUEBAXByJEuVSgikoH5+MBnn8GBA9ClCxQurE7lGZjWjIqIiFzk0q0y4k5k48C0mymaM4Zvs7TGyeoH06ZB9uwuVikikoH5nX+drV0b2reH335zuyJJJQqjIiIiF7l4qwxz1pcj4Tfjezor64KeJevWDbbjY+nSbpcpIpKx5coFs2dDYCDcdx9s2uR2RZIKFEZFREQu0apWIEtCGlFzR1Ni9uRh5ZMfUmDeJBg4EJo0cbs8EZHMoXBhmDfPjpQGB8Pu3W5XJB6mMCoiInIJY6BnT5g5E6b2+IVqn/WGli2hXz+3SxMRyVzKloXvv4foaLj3XjhyxO2KxIMURkVERC4SHw9du8IHH8Brz+yl7ZR2cMMNMGGCbawhIiLeVasWfPMNbNliPxiMSbrruaQ/+qsqIiJy3rlz8OST8NFH0O+FWF7b+CAcPQozZkDevG6XJyKSed19N0ycCEuXwqOPQlzcle8jaZ7CqIiICHDmDLRrB198AYMGwVuxIThLl8Knn0K1am6XJyIiDz0E771nR0m7drVrKiRd0z6jIiKS6Z08Ca1awfz5dvvQ548NgpEjoEcPePhht8sTEZEEPXpAVBS8/TYULw6vvup2RZICCqMiIpKpRUdD8+awbBmMHw9PHhoOL78Mjz0G77zjdnkiInKpIUNg71547TUoUgSeecbtiuQ6KYyKiHhR+JpIwuZtISo6huIB/oQEV6JVrUC3y8q09u+3uwVs2ABffQVt9o6CF16w83XHjQNfX7dLFBGRSzmOXUJx4AD873/2RfzttyF7drcrk2ukNaMiIl4SviaS0BnriYyOwQCR0TGEzlhP+JpIt0vLlP75B+64wzZn/O47aBM9Frp1s50av/wSsujzWhGRNMvPzzaX69kT3n8fbr7ZhlJJVxRGRUS8JGzeFmJiE3f/i4mNI2zeFpcqyry2b4cGDWDPHrufevCBL6BzZ2jaFKZOtW9yREQkbcueHd59F2bPttN2g4JsO3Q1Nko3FEZFRLwkKjrpfdGSOy6p448/bBA9cQIWLoQGe7+2+7ncdZf9lD1bNrdLFBGRa9GsGfz+O9x5Jzz3nO1Id/Cg21XJVVAYFRHxkuIB/td0XDxv5Ur7XsXHB5YsgTq7v4VHHoHbbrNzdf31/0JEJF0qWhTmzLEjpXPnQvXqsGCB21XJFSiMiohcp/A1kdQfspAy/WZTf8jCK679DAmuRI4sDsWP7SfXmVMA+Pv5EhJcyRvlZno//QSNGkHevHbP9Cp/fW8bFdWubad45czpdokiIpISPj52DemyZfbFvkkT6NcPzp51u7LUs2+fDd/plLoziIhch4RmRAlrQBOaEQG2O+7Jk7B1K2zebC9bttBq82ZabN6M75kz9j75i+HUrElxv1vhn5pQsyaULGm7BIpHzZkDbdpAmTLw448QuHkBtG4NVavaP+J58rhdooiIeEqtWhARAb17w9ChdoR00iSoUMHtyjxn1SoYMcL2OfD3t00Q0uHsHse4vMA3KCjIREREuFqDiMi1qj9kIZHRMRQ8eYQKB/+m3KHdlDu8myrH9nDLmf3w99//3thxbAq68UZ7qVABDh2CtWth3TrbTSfhtTggwIbSGjX+/VqlitYxXqczZ2DwYBg0yP4q586FgpuW2kZFZcvCokVQsKDbZYqISGqZMQM6dbKjox98YHsEpNcPfWNj4ZtvbAj99VfIlQs6dIDnn4eKFd2u7rIcx1lljAm69LhGRkVErsOZyD0MXjKR9r//iA82SJ7I6s/O/IFwd4N/g+eNN0L58pA9OydO2A9qV64E4weFmkGhJ6FIzhMUP7SeArvXkm3TOpx1a2HMGDhlp/KSJQtUrmz3IenY0X7iK1e0bJn9dW3cCI89Zt+D5N20zDa6KFkS5s9XEBURyehat7bbvjz+ODz1lP1U8uOP7Ye/6cXBgzB6NHz4IURG2g9T333X/jx587pdXYpoZFRE5Fqc/2T1RP9XyRZ7mom172NBubrsKFCCfbkKEJgvB7/0a0RcHGzaBMuX28uyZXb7s/j4y58+e3YoVAiKFIyjRs7t1HTWceOZtZSJXssNuxbhe/Y01KljtyF5+GGPTy8NXxNJ2LwtREXHUDzAn5DgSnbacTpy4gS8/DKMHAklSsAnn8C99wKrV9tFowUL2gWkgenr5xIRkRSIi7NTdl991f5x+PJLqF/f7aoub906+8fsyy/tVJ/GjaFHD/tHzdfX7equSXIjowqjIiJXa84c6NULtm5l7+2NeKraw2zKUwyAuBPZMPvzc2uuChz5KzcrV8Lx4/Zu+fLBLbf8e7n5Zjvr9sCBxJf9+/97LOH4qVMQwBEe40ue9x9DxZjficueA/Nge7L8rzPcemuKpx1dug4WbIOlwa2rpZtA+uOP0KUL7NoFXbvaKbq5cwPr19utW3Lntm10S5VyuVIREXHF8uW2i/quXTaYvvSSnYGUVsTFwcyZdiruTz/ZdaBPPAHdu9tlO+mUwqiIyPXassU2QZgzx67JePddaNaMEV/vY+BbhiM78xB3LAdg/57VqPFv8Lz1VrtENKXLU06etB+Qzp8P8380nPttJR3iPuVhJpObExwochOnH+1E8b6P41u4wHU9RsI62EsFBvjzS79GKfsBPOByo7aHD8MLL8D48VCpEnz6Kdx++/k7LlhgR5GzZrV/2MuVc+1nEBGRNODYMejWDT7/3P6h7tjRfmBZrpx760mPHLF/vEaNgr/+gtKlbY1PPw3587tTkwcpjIqIXKujR+GNN+wUmRw54LXXoFs3jpzMyoAB9u9F7ty2c/ytt9rwWbu2d5rZnThhB/iWzjlO1m+m0ixqDLewgjNkZWXJ1hxp05kqz91FuQpXv4NXmX6zSeovggP8OaS5x2q/HpcbtY3bGUjXrnZJTb9+dopu9uzY4eR+/eD9921C/fZb+1VERATs9NeQENuJFuz03bvugoYN7dcyZVI3nB47ZmfufP65vZw6ZR+3e3do0SJtjdimkMKoiMjViouDceOgf3+bcDp2hIEDiStYhM8+s4cPH4ZnnoE334QC1zcQ6VF798Kqcb+T9fNPqbv5cwJMNNspx4z8nYhu25kmDxWgQYPL/11LyyOjSdV27kQ2YhZX5/CGwtSuDWPH2gbEgJ2G9cQTdnudHj3sfN102PJeRERSmTF2BtSiRbB4sb3s32+vK1XKhsOEgHrDDdf3GCdP2m56GzbAH3/Yrxs2wD//2OuzZbOd9p5/3k6vyoAURkVEknHx9M+m0dsZ/NNYAjavt40NRo6E2rX59Vf7N2L1amjQwB6+EHzSGHMqhj0fzsB8MobA7T8RQ3Ym8gTj8/Sg/P1VaNkSgoPPr6W8SFpeM3rxqK0xcOL3khxZVBnifBj6li+9ep0P2mfP2tHswYPtJ9zjxtmmRSIiIlfDGNuB8OJwevCgva506X9HTRs2/G//gZgYe9+EsJkQPnft+vc22bLZTvs33WT3ur7pJvt+Iy18sp2KFEZFRJKQEMByRB/itQWjuX/TEvbkLsju/q9Tt++zRO1x6NsXvvjCNl8dNgzat09HW5T98Qex74zEZ9Ln+J49zcIs9xB2rieL/IJp1NiHli3h/vuhmO3DdF3ddL3Rgbf+kIXsPnyaM//k5+hv5Tn9V0GylTxElQe3snpYvQs/K48/bvdvfeopu7Y3nbe8FxERl8XH21HNi8Pp4cP2ujJlbIOCY8ds8Nyx4999w/387NKQm25KHDzLls1Q02+vlsKoiEgS6g9ZSEzUXqZMCqV09B4+vqUtH9/ShoD8+Wnl15A337R7TIeEQGgo5MzpdsXX6fweZeaDD3D27GFf/ht536cH7x58nFPk5OaboWVLe6lS5erDdmqPpp4+bZs2vTvmJD/96EdcTFacrLHka7iZQkGRDGlTjVbVi8Lw4XaxaECA3YutZcsUP7aIiMh/xMfbDz8Twumvv9oGQxcHzptust0L/fzcrjbNUBgVEUlCrR6TmTS5Pzcc2cNT7V5jWanqxOwoxOEFVTh3JBctW8I772SgBqxnz8LXX9tRw1WrOJcnHytqdmHw0a7MWlcSsB/a3n471K1rLzVqnG8IlITUWGd67JhtXDxjBnz/vW3WlDcv1LjtFPvy7yCmcCQlCme1I7B5TkOHDvDzz3Zj848/thu1ioiISJqhMCoicqkjR9hS9WZu2PcXT7d9jcV5buPIwirE7CiCf6GTfPN5ToKD3S4ylRhjP8197z2b+hyHmOZtmXtjT8ZvvpXly2HfPntTPz+oVu3fcFq3rh09zZLFcx149++326p9840dCT17FooUgVat4IEH7NKcrFkvqX/0aLufS5Ys8MEH8Oij6Wj+tIiISOahMCoicrGjR6FxY+LW/c7/2r7CjMP3c/jHqjhZ4ih0xw4+GJSbdje727THa/76y4a5MWPs7+XmmzGPPsaeW1qxLLIkK1dy4XLsmL1Ljhx2G5sd5m9i8x8ia7FofLLHQryDiXcomtufr5+pz7lzJLrExSX+ft06m4V//tnOfCpTxg5wPvCA3S7H1zeJeqOibIfjuXOhcWP47DMoWdKrvzIRERG5egqjIiIJjh2z7WRXrYIZM+i05HbGhgXgX2Y/Nz28lZfalnG9e6wrTpyACRPgo49sIwaAoCCbDB94gPhKldm2jUThdNVqw9kzKRuNrFbNPkTr1lC9ejKDm3Fx8NtvMGuWHRE9fRrCwuDZZ8Hn6vdSFREREe9TGBURARu4mja1+1B+/TVh21rx4os2DE2ZcslU0HQsxR1ut2yxc2a/+QZWrLDHKlW6EEypWxcch9hY+GD6fkZ9fYgjR+PJmzMLwVULE1Q2H76+dgZtUpeE60qVusx63GPHYN48+O47u4j00CF7pyZN7PTiihVT+msSERERL1AYFRE5dQqaNbNzQidPZuCWdrzyCjz0EEycmHGa3nm8w21kJISH22C6eLEdpQwMtAs6W7eGO+7wXJv6P/+04fO77+Cnn2wr4/z57f+3Fi3siLa2axEREUlXFEZFJHOLibFhZtEizOdf8Oqmhxk4EJ54wi45THJtYjqVGh1uLzh82E6V/eYbO2oZE2PDYosWtqV9njz/veTNa7/mzv3f0BoXB8uW2fA5a9a/04MrV4b77rPnrVcvU+7JJiIiklEkF0b1111EMr7Tp+3U0oULMePG03ftw4SFQadO8MknGW/JYVQSQfRyx69J/vw2wT/xhB1pnjfPBtNvv7XrTa8kR47EQXXXLrsHapYsdoS1UycbQsuXT3mtIiIikqYpjIpIxnb2LLRtC/PmYT4dS8/VTzByJHTtCiNHZrwgClA8wD/JkdHiAf6efaAcOf5dQxofDydP2nWe13KpXBmaN7fTbwMCPFufiIiIpGkKoyKSccXGQvv2MHs28R9+zHMrn+aTT6B3bxg2LONuSRkSXCnJNaMhwZVS70F9fOw03Ny57XpSERERkStQGBWRjOncOXjkEQgPJ37E+3Rc8Qzjx0NoKAwalHGDKHChSVGKuumKiIiIpDKFURHJeOLi4PHHYdo04sKG88TybkyaBK+/Dq+8krGDaIJWtQIVPkVERCRN89hqKcdxcjuO857jOH85jhPjOM6vjuPU9dT5RUSuWrduMGUKcW8N5aHlvZg0CQYPhldfzRxBVERERCQ98GTrjk+BYOBJoBrwAzDfcRx9NC8i3jNjBnz8Med69aH1sheZNg2GD4d+/dwuTEREREQu5pEw6jiOP9AG6GeMWWyM2W6MGQBsB571xGOIiFxRZCR07kx8rTq03jiImTNh1Cjo1cvtwkRERETkUp5aM5oF8AVOX3I8BrjdQ48hIpK8+Hjo0AFz+jTd8n3JrB+yMmaM3bZSRERERNIej4RRY8xxx3F+A152HOcPYC/wMFAPOzqaiOM4XYAuAKVKlfJECSKS2Y0YAfPnM6/Vx3wUXokRIxRE07LwNZHq9isiIpLJOcYYz5zIccoBnwF3AHHAamArUMcYUzm5+wUFBZmIiAiP1CAimdS6dXDzzRys25Qiv4bTtp3DlClqVpRWha+JTHIf1MGtqymQioiIZECO46wyxgRdetxjDYyMMTuMMXcCuYCSxpibAT9gp6ceQ0Qyl/A1kdQfspAy/WZTf8hCwtdE/vdGMTHw6KPEBeTnrm2fUr6Cw6efKoimZWHztiQKogAxsXGEzdviUkUiIiLiBo/vM2qMOQmcdBwnH7a77ouefgwRyfguHT2LjI4hdMZ6gMSjZ/36wYYN9K/+PTu3FWL5fMid242K5WpFRcdc03ERERHJmDy5z2iw4zj3Oo5TxnGcJsAiYDMwzlOPISKZx1WNns2dCyNH8lvQ87z9e1M+/BCqVfNyoXLNigf4X9NxERERyZg8uc9oXuADbACdCPwMBBtjYj34GCKSSVxx9OzAAejQgeOlb6JRxFCefho6dPBefd50VdOV05GQ4Er4+/kmOubv50tIcCWXKhIRERE3eGyarjHmK+ArT51PRDK34gH+RCYRSIsH+IMx0KkT5sgRmvv/QMXq/nzwgQtFesFVT1dORxLqVjddERGRzM3ja0ZFRDwhJLhSkh1XQ4IrwZgxMHMmI0u/w9rD1Yn4Gvwz6AzPy01XTs/hrVWtwHRdv4iIiKScwqiIpEnJjp7lPAm9erGlVGN6/dWTqV9BxYouF5uK1OxHREREMiqFURFJs/4zenb2LNzWkjM+2Wn093i6Pe9Du3bu1ecNl52uLCIiIpKOebKBkYhI6howAFatomPcGErcHMiwYW4XlPrU7EdEREQyKo2MiohXhK+JTFnDmiVLMEOGEJ7/aeaY1qz5CrJmTb160wo1+0m5FD/3REREJFU4xhhXCwgKCjIRERGu1iAiqevSjrBgR/cGt652daEgOhpq1GB/tB9lj61l6qxcNG+eevVKxpHi556IiIikmOM4q4wxQZce1zRdEUl1l+sIe1W6diV+dyQtjn1J91AFUbl6KX7uiYiISKrRNF0RSXUp6gj75ZcwaRIDs7yBf4NbeOMNDxcnGZq6EYuIiKRdGhkVkVSXXOfXK3aE3b8f07Urq7Pfxif5Q5k8GbLoIzS5Btf93BMREZFUpzAqIqnuujvChoYSd+wkj50ZyxdTslCsWCoWKRmSuhGLiIikXRpjEJFUd10dYVesgM8+41368NDrN9KwoZeKlQxF3YhFRETSLnXTFZG0Jz6ec3Vv5dDaf2hdZQuLV+fBz8/tokRERETkeiTXTVcjoyKS9owfT5bVK3nRmciIcQqiIiIiIhmRwqiIpC3R0Zx9oR8ruY1CvR4j6D+foYmIiIhIRqAwKiJpyrmXXiNL9EGGFJ/LlDcct8sRERERkVSibroiknb88QfOR6MYTRe6j69NzpxuFyQiIiIiqUUjoyKSNhjDiaee54zJy/r2g/hfk+s7TfiaSHVOFREREUkHFEZFJE2In/o1uSIW80auD3ljVIHrOkf4mkhCZ6wnJjYOgMjoGEJnrAdQIBURERFJYzRNV0Tcd/IkJ599gTXUpM4nXShwfVmUsHlbLgTRBDGxcYTN2+KBIkVERETEkxRGRcR10X0Hkzt6N5Prvc+DD/te93miomOu6biIiIiIuEdhVERcZbbvIMeHYUzJ8ijdptyOk4IGusUD/K/puIiIiIi4R2FURFwV+WAvzpisxAx4m1KlUnaukOBK+PslHln19/MlJLhSyk4sIiIiIh6nBkYi4pqjU76nxJrvGFVqKP/rVzzF50toUqRuuiIiIiJpn2OMcbWAoKAgExER4WoNIuKCM2fYU6gax487nFuznio1s7pdkYiIiIikAsdxVhljgi49rmm6IuKKLf97l2LHt7Hi0REKoiIiIiKZkMKoiHjdya2RlJgwkAW5W9JubFO3yxERERERFyiMiojXbW4Rgq85R77PhpMtm9vViIiIiIgbFEZFxKs2fryEOlsns6jui9RuW9btckRERETEJQqjIuI1Z0+dw7fX8+z2LcXt3/VzuxwRERERcZHCqIh4zYL2n1Dp9O/sC3mH3EVyuF2OiIiIiLhIYVREvGL7soPcMusV/ijciDpvtXG7HBERERFxmcKoiKQ6Y2B9+4Hk5ShFp44Ax3G7JBERERFxmcKoiKS6Hz/eQfO/P2RTvY4UvKuq2+WIiIiISBqgMCoiqerUKTjd5yXiHD9unDLA7XJEREREJI1QGBWRVPX58yu4/9RU9j/+AllKFXe7HBERERFJIxRGRSTV7NhuqDIuhKPZC1P6gxC3yxERERGRNERhVERSzZePzKKBWYJ5dQDkzu12OSIiIiKShiiMikiqmP3tOdqu7MvhQhUJ6NPJ7XJEREREJI3J4nYBIpLxnD4Nv3b6jOZs4tyoGeDn53ZJIiIiIpLGaGRURDzuvYEn6HbwNaKr1idL21ZulyMiIiIiaZBGRkXEo3btgrNDh1OMvTB6BjiO2yWJiIiISBqkkVER8agBz+6j97m3OXVvG6hXz+1yRERERCSN0sioiHjM3Llwy9wB+PucwXfEYLfLEREREZE0TGFURDzizBl473+bmcUYeOZZqFDB7ZJEREREJA3TNF0R8Yh334Vn/gqFHDnwHfCK2+WIiIiISBqnMCoiKfbPPzB/wM88QDhZ+veFwoXdLklERERE0jiPhFHHcXwdx3nTcZw/Hcc5ff7rQMdxNA1YJBN4obdh0NkQzhUpDr16uV2OiIiIiKQDngqLfYGuwJPAeqA6MAE4A7zpoccQkTRowQKInzadW1gGb42FHDncLklERERE0gFPhdHbgO+MMd+d/36X4zgzgVs8dH4RSYPOnoXe3c7ybZZQ4ivehM+TT7pdkoiIiIikE55aM/oz0NBxnBsBHMepAjQC5njo/CKSBo0cCQ02j+aGc9vxCXsbfH3dLklERERE0glPjYwOBXIDGx3HiTt/3kHGmA89dH4RSWOiouDdAUfZmPV1uL0R3Huv2yWJiIiISDriqZHR9sATwCNA7fP//ZzjOB2TurHjOF0cx4lwHCfiwIEDHipBRLwpJAS6n36bvGcPwttvg+O4XZKIiIiIpCOOMSblJ3Gcf4BhxpgRFx17GehgjCl/ufsGBQWZiIiIFNcgIt6zZAk8cudu/sxSAb8HW8OXX7pdkoiIiIikUY7jrDLGBF163FMjozmAuEuOxXnw/CKSRsTGwnPPwfCcr5LFJx4GDXK7JBERERFJhzy1ZvQ7oJ/jOH8CG4BaQG9goofOLyJpxMiR4GxYTztnPE7v3nDDDW6XJCIiIiLpkKfC6PPY/UQ/BAoDe4AxwBseOr+IpAG7d8OAAbC4UAhObF7o39/tkkREREQknfLINFpjzHFjTE9jTGljjL8xpqwxpr8x5rQnzi8iaUPv3nDXme+pc2Aeb9ZuS/3RawlfE+l2WSIiIiKSDnlqZFREMrgffoAZX59jfY4e7MpZjIm1mxMbHUPojPUAtKoV6HKFIiIiIpKeqMGQiFzRmTPQrRs8k/NjKp/axpC7OhDr6wdATGwcYfO2uFyhiIiIiKQ3CqMickXDhsGebccZcO5VVpSowtyKtyW6Pio6xqXKRERERCS9UhgVkcvatcvu3jL+xqEUOnOEQQ07guMkuk3xAH93ihMRERGRdEthVEQuq0cPKMk/tP7zHf5p2pKtpaskut7fz5eQ4EouVSciIiIi6ZXCqIgk67vvYOZMmFG5Pw6Gkh+PYHDragQG+OMAgQH+DG5dTc2LREREROSaqZuuiCTp1Cno3h3alYngptVfQL9+ULo0rUqrc66IiIiIpJzCqIgkafBg2LXLsKbGC1CoEISGul2SiIiIiGQgCqMi8h/btsHbb8O7d35LwE9L4MMPIU8et8sSERERkQxEa0ZFJBFj7J6iubOdpds/L0LlytC5s9tliYiIiEgGo5FREUlk+nT44Qf4qfVHZJmxDWbPhix6qRARERERz9LIqIhccPw49OwJDaoeocHiN6BxY7j3XrfLEhEREZEMSGFURC544w2IjIQp1QbiHDkCw4aB47hdloiIiIhkQAqjIgLAhg3w3nsQ+uAOik97H556CmrUcLssEREREcmgFEZFBGPguedsw9zXTvUFPz948023yxIRERGRDExhVET48ktYsgTGdfyZbLOmw4svQvHibpclIiIiIhmYY4xxtYCgoCATERHhag0imVl0NNx4I5QpHc+v1MPZvRu2boWcOd0uTUREREQyAMdxVhljgi49rpFRkUzu5ZfhwAGY1HIqzooVMGiQgqiIiIiIpDqFUZFMbPFiGDUKej8bQ5lP+kHNmvDEE26XJSIiIiKZgHayF8mkTpywDXMrVIBBhUfA33/DuHHgo8+oRERERCT1KYyKZFIhIfDXX7D8u/1kffgtuO8+aNTI7bJEREREJJPQEIhIJvTjj/Dxx/DCC1D3u1fh1CkIC3O7LBERERHJRDQyKpLJHDsGHTvaDroDm/8GjUZD9+72gIiIiIiIlyiMimQyvXtDZCT8tiSWbP/rAoGB8OabbpclIiIiIpmMwqhIJvL99zB2LISGws0/D4c//oDwcMid2+3SRERERCSTURgVySSOHIFOneCmm+C1x3dCndehVSto2dLt0kREREQkE1IYFckkevaEfftg5reGbL2egyxZ4P333S5LRERERDIphVGRTGDmTJg4EV59FepsmwLz5sHIkYQfcAj7YiFR0TEUD/AnJLgSrWoFul2uiIiIiGQCCqMiGdyhQ/DMM1CjBrz07GGo0RPq1iW8XktCZ6wnJjYOgMjoGEJnrAdQIBURERGRVKd9RkUyuOeft4F0wgTI+mo/+83o0YTN334hiCaIiY0jbN4WlyoVERERkcxEYVQkA5s+HSZPttNzaxxbCmPG2MWjNWsSFR2T5H2SOy4iIiIi4kkKoyIZ1IED8OyzUKcO9O111s7VLV0aXn8dgOIB/kneL7njIiIiIiKepDAqkgEZA889B0eP2um5fu++DZs2wYcfQs6cAIQEV8LfzzfR/fz9fAkJruRGySIiIiKSyaiBkUgGEb4mkrB5W4iKjiHr36XZOq0qQ4bATX5bYeBAaNcOmjW7cPuEJkUJ91E3XRERERHxJoVRkQwgfE3khc64cSey8fc3FckeGE35hifsXN3s2WHEiP/cr1WtQIVPEREREXGFpumKZABh87YQExuHMXDoh6rEx/qS/961bBz2ASxcCEOGQLFibpcpIiIiInKBwqhIBpDQAffkxkBithUl3x1bKOwfRddZH0G9etCli8sVioiIiIgkpjAqkgEUD/An9nBODv94E9kCD5M76E/6LxpHnjMn4ZNPwEf/1EVEREQkbdE7VJEM4LnbbuTgjCAcn3gKtljLbf/8Trs/5rPzyf9BtWpulyciIiIi8h8KoyLpXFwcfDG4OHFHc1Ll8T/ImfMoQ+d/yMkSpak4Kszt8kREREREkqQwKpLO9e0Lc+fCRx86rP+kDltyrqXUwd3kHDsa/P3dLk9EREREJEkKoyLp2IQJ8M470K0bdO4MbNoEgwfDI4/APfe4XZ6IiIiISLIURkXSqWXLbJPcRo1g+HDsfN0uXSBXLnj3XbfLExERERG5rCxuFyAi1273bnjgAShZEr76Cvz8gP6vwM8/2+HSwoXdLlFERERE5LIURkXSmVOnoFUrOHkS5s+HAgWA8HA7PbdLF3jiCZcrFBERERG5MoVRkXTEGOjYEVavhm+/hZtuArZuhSefhLp1YeRIt0sUEREREbkqCqMi6ciQITBlih0EbdECOzzaurWdpzttGmTL5naJIiIiIiJXxSMNjBzH2eU4jkniMtsT5xcRmDkTXnrJNsrt2xc7TNq5s+2gO2UKlCrldokiIiIiIlfNUyOjdQHfi74vBqwCvvLQ+UUytQ0b4NFHoU4d+PRTcBxg5PsweTK89RY0bux2iSIiIiIi18QjYdQYc+Di7x3H6QgcQ2FUJMUOHYL777c7toSHg78/tmvuCy9Ay5bnh0lFRERERNIXj68ZdRzHAToCXxhjYjx9fpHMJDYW2rWDyEj46ScIDAT27LEHb7jBbuPio+2CRURERCT9SY0GRk2AMsCY5G7gOE4XoAtAKa1zE0lWr16waJHNnLfcgk2n7dvDsWPwww+QN6/bJYqIiIiIXJfUGFLpDKw0xqxL7gbGmNHGmCBjTFChQoVSoQSR9O+TT2DUKOjT56KtQ198EZYutQtHq1VztT4RERERkZTwaBh1HKcw0JLLjIqKyJV9/TV06wb33mu3cwFsx9z33oPu3eHhh90sT0REREQkxTw9MtoBOANM9vB5RTKNkSPtTNxbbrHNcn19se10O3aE+vUhLMztEkVEREREUsxjYfR846JOwBRjzAlPnVcks4iPt41xe/SwTXJ//PH8ktCjR6F1a8iTB776CrJmdbtUEREREZEU82QDo7uACsBjHjynSKZw9qwd+PziC3j2WXj//fMjosZAhw6wY4ftZFS8uNulioiIiIh4hMfCqDFmEeB46nwimcXx49C2rW2OO3Ag9O8PTsK/pLfftpuLvvsuNGjgZpkiIiIiIh6VGlu7iMhV2rcPmjWDdetg7Fh4+umLrlywwCbT9u3t3F0RERERkQwkNbZ2EZGrsG0b3HYbbN4M336bOIjOn76YIy3bsDVfCRpXfpzwtVHuFSoiIiIikgo0MiqSisLXRBI2bwtR0TEUD/AnJLgSrWoFsmIFNG9ub7NoEdx887/3+Wn8TIKee4xYH1/+90B/dsZA6Iz1ALSqFejCTyEiIiIi4nkaGRVJJeFrIgmdsZ7I6BgMEBkdQ+iM9bzywUEaNoTcueGXXxIHUWbN4pYuD3Ikey5aPzaMnQVKABATG0fYvC2u/BwiIiIiIqlBYVQklYTN20JMbFyiYwdWF2Ngj/xUqgS//goVK1505WefQatWbClQiraPhvFPQNFE942KjvFC1SIiIiIi3qFpuiKp5OLwaAwc+6080Usrkf2GA/z0UyFy577oyrfegpdfhuBgXqj3PIeSyJ3FA/y9U7iIiIiIiBdoZFQklSSERxMPh3+8ieillchZJZKaHf/4N4jGxcHzz9sg+thjMHMm3e6vib+fb6Jz+fv5EhJcycs/gYiIiIhI6lEYFUklIcGVMHsLsPfz+pxYcwN5bt5BiQfW07f5+bm5p0/DQw/BqFHQpw9MmABZs9KqViCDW1cjMMAfBwgM8Gdw62pqXiQiIiIiGYqm6Yqkgr//hqlvB/L3lECy5jlNwRZrqHjbEUKCz4fKo0ehVStYvBjeeQd69050/1a1AhU+RURERCRDUxgV8aCTJ+HttyEszC4FfeUV6Ns3Ozlz1vr3RlFRcO+9sGkTfPklPPKIewWLiIiIiLhEYVTEA4yBSZOgb1+IjIT27WHoUChd+pIbbtkCwcFw6BDMng1NmrhSr4iIiIiI27RmVCSFVq6E+vVt/6EiRWDpUpgyJYkguny5veGpU3Z6roKoiIiIiGRiCqMi1ykqCp58Em6+GXbuhLFjbTC9/fYkbvz999CoEeTNazcYrVPH6/WKiIiIiKQlmqYr6V74mkjC5m0hKjqG4gH+hARXStXmPzExMHw4DB4MsbF2am7//pAnTxI3PnbMNigaNAhq1IA5c+zwqYiIiIhIJqcwKula+JpIQmesJyY2DoDI6BhCZ6wH8Hgg3bcPwsNhyBDYtQseeMA2KipXLokbx8TAhx/axHroEDz4IHz6Kf9uMCoiIiIikrlpmq6ka2HztlwIogliYuMIm7fFI+ffudMObDZoAMWKwf/+Z0dAFyyAGTOSCKLnzsGYMVChgt07tHZtO3d36lQFURERERGRi2hkVNK1qOiYazp+JcbA2rV2BPSbb2C9HWSlRg147TW7NWj16uA4l9wxPh6+/tru5bJtG9x6K3z+OTRseF11iIiIiIhkdAqjkq4VD/AnMongWTzA/6rPERcHP/9sA2h4uJ2C6zi2EdHw4dCyJZQtm8ydjbHNiV56yabYqlXh22+hRYskEquIiIiIiCRQGJVU4a2mQiHBlRKtGQXw9/MlJLhSsveJjYUdO+yo5/ffw3ffwcGDkC0bNG4ML79ss2Thwld48J9/htBQ+7VMGTsS+vDD4OvroZ9ORERERCTjUhgVj/NmU6GE8yUVfGNiYOtW2LgRNm369+u2bTaQgl3/ed99dvpt06aXX9aZELDzbtnAy799yW1blkPRorZRUceOkDWrR382EREREZGMzDHGuFpAUFCQiYiIcLUG8az6QxYmOXU2MMCfX/o18vjjRUfb0Hlx4Ny40TYfSnh6+/jYZkNVqkDlyv9+rV79KjLkgQMs+3I2a6fNpcY/G6n393qis+di7G3tqPhmf1rcVt7jP5OIiIiISEbhOM4qY0zQpcc1Mioe5+mmQgCnTsH27TZ0bt1qRzcTvh448O/tsmaFihWhTh147LF/Q2eFCpA9+1U80Llzdv7ub7/Zy7JlsH07twJBjg+bCpdhxG0PM7ZuS45lz0Xgkr8VRkVEREREroPCqHhcSpoKHTxo819C6EwInLt3X3Ku4jZ0tmplv1aoYENn2bKQ5Vqe1fv22QdMCJ4rV9rkC1CkCNSrB5060f53WFe0PKf9EifalARsEREREZHMTGFUPO5amgrFxtocOG+evWxYdZpAduNHLIXyxlK2ZCx3VjvHDU1jKVX8HCWLxRJYKBZ/v3P2zrGxdjTzcCz8eMoGyVOn4OTJK//38eOwf78tJEsWqFULOnWy27LUqwelS1/oiLt7yEJOp7Brr4iIiIiI/EthVDzuck2FwHayXTDzJFu+3czx5Rspc3ojt7CRLlk3UII/8SXenujo+csf11GEvz/kyGEvOXP++99589ph1YTvK1a0wbN2bXufZFxP114REREREUmeGhhJ6jp6lJMRm9j67UYOLdlIlm0bKX1qE2XYdeEm8Vn8OFayLL9lL8zWfCX4O6AoZ7L44ePnxyO3l+PWSkXtyKWfn70k998JAdPf33Ys8jBvbVcjIiIiIpKRJNfASGFUPO/cOc5Nnc7hl9+h8K6VFw6fJhtReW4ktkIVCtxehQINquDcVAXKlaP+O0u92oFXRERERES8Q910JfWdPEn82M84NXA4uQ7s4ggVmVx4EAENqlGpVWVqtS5D2Ry+Sd41NTrwioiIiIhI2qUwKim3fz9m5PvEjvyQrMcPs47b+KrUuzR5/366t/BJ6AF0WSnpwHs9NOVWRERERMRdnl9YJ5nH1q3wzDPElyyFGTSI2cfvoHWRX9g58ReG72zFffdfXRAF2yDI3y/xqGlqNQgKXxNJ6Iz1REbHYIDI6BhCZ6wnfE2kxx9LRERERESSpjAq1+6336B1a8yNN3L20wmMOfsktwVs4u/3vmHyX7fx+OPgm/Rs3GS1qhXI4NbVCAzwx8GuFR3culqqjFaGzduSqCsuQExsHGHztnj8sUREREREJGmapitXFL4mkmHfb6LKqp/ouiqcGn/9wcls+RhBfz7N/jxPhBThxxcgd+6UPU6rWoFemSqr9akiIiIiIu5TGJXLCl8TyQdjf2T0129QZf+f/JU1kB4+7zLhXEeefD43y16CwoXdrvLaeHt9qoiIiIiI/Jem6cplzRk9gynjelLs8CEeyzKBcmf/ZGzltpTrvZ4RI9JfEAXvrk8VEREREZGkaWRUkvfVV7w/5gV2+xan2bl5/F0mH4Xv+o2shY9z2O3aUiBhKrC66YqIiIiIuEdhVP7LGBg8GF56iVXZb+X+099x7vZDFL5t5YXuuOl9Squ31qeKiIiIiEjSNE1XEjt7Fjp2hJdeYk7AIzSOXYRfiz0E1N9+IYhqSquIiIiIiKSUwqj868gRaNoUxo3j3Tyv0T72C8LnZOej1/N7ZcsVERERERHJPDRNN5MJXxOZ9FrJnTuheXPit+/gWf+JfJfzcZbOgZo1ATSlVUREREREPEthNBMJXxNJ6Iz1xMTGARAZHUPojPXkXxfBHSGdOHM6nmbMZ88Nd/Db91C6tMsFi4iIiIhIhqVpuplI2LwtF4JogsbrFnJL53YcIYCqJ5Zx7rY7+OUXBVEREREREUldCqOZSFR0zL/fGEPXX6fy/ndhRGSvSfmDy6j9YAXmzYN8+dyrUUREREREMgeF0UwkYTsWv7hYhs15j5ClnzMpdxsanVjKUy8UYPJkyJ7d5SJFRERERCRTUBjNREKCK1HAnGXiV6/S9o8FvJm7L48e/4rH+sQwbBj46NkgIiIiIiJeogZGmUirmsUJCp1AsX828GSOMXwe8xQvhh1haJ8CbpcmIiIiIiKZjMJoZjJqFCXmfcugnG8xK1snlv4A9esriIqIiIiIiPcpjKZTye4XmpxlyzC9e/NT7vt4z68vv/wCN97ovXpFREREREQu5rFVgo7jFHMcZ4LjOAccxzntOM5Gx3Hu9NT55V8J+4VGRsdg+He/0PA1kUnf4eBBzIMPsj9rCdqcmMjkqT4KoiIiIiIi4iqPhFHHcQKAXwAHaA5UBp4H9nvi/JJYUvuFxsTGETZvy39vHBcHjz5K3J793HtyGi8Ozkfjxl4qVEREREREJBmemqb7IrDHGPPERcf+9NC55RKJ9gu90vE334QffqCrM5qybWrz4oupXJyIiIiIiMhV8NQ03VbAcsdxpjqOs99xnLWO43RzHMfx0PnlIgn7hV7x+Ny5mDfeYEq2J1lSsRPjxoH+j4iIiIiISFrgqTBaFngO2AkEAyOAIUDXpG7sOE4Xx3EiHMeJOHDggIdKyDxCgivh7+eb6Ji/ny8hwZX+PfDXX5hHH2V79mr08PuQb8Idcuf2cqEiIiIiIiLJ8NQ0XR8gwhgTev77NY7jVMCG0Q8uvbExZjQwGiAoKMh4qIZMI6FrbrLddM+cgXbtOH38HM1ip/HxjBxqWCQiIiIiImmKp8LoHmDjJcc2AT08dH65RKtagclv5dK7N6xcyaNMp11oBR54wLu1iYiIiIiIXImnwugvQKVLjlUE/vLQ+eVqTZoEH37Iuz4vcOLu1rz5ptsFiYiIiIiI/Jenwui7wK+O47wETAVqAd2B/h46v1yNDRswnTuzMuvtjCo6mGWTwNf3yncTERERERHxNo80MDLGrMR21H0Q+AMYBLwCfOiJ88tVOH4c06YNR87l5iFnKl9940fBgm4XJSIiIiIikjRPjYxijJkNzPbU+eQaGAOdOmG2bqO1WcCr44pTu7bbRYmIiIiIiCTPY2FUXPT++/DVV4QyhCrP3kWHDm4XJCIiIiIicnkKo+ndb79hXniB2b73syToRX56z+2CRERERERErkxhND07cID4dg+y2ylF73wTWDTdIWtWt4sSERERERG5MoXRdMz06Encnv20ZhmfTgsgMJltR0VERERERNIaj3TTFRf88gvO5EkMjQ/h8eG1uOMOtwsSERERERG5ehoZTY/i4jj7bHcOOIGsbhLK9O5uFyQiIiIiInJtFEbTo3HjyLp+NS/5fck7H+fEcdwuSERERERE5NoojKY30dGc6dOfldSnwqsPU6aM2wWJiIiIiIhcO4XRdObcq2/gd/Qgw0vPZXKIhkRFRERERCR9UgOj9GTTJpxR7zOWjnT7rDbZsrldkIiIiIiIyPVRGE0vjOFkl14cj8/J6taDaNTI7YJERERERESun6bpphPmu1nk/Hke/bK9y2ujCrtdjoiIiIiISIoojKYHZ85woksv/qEypYZ2pWhRtwsSERERERFJGYXRdOD0kPfIvW8Ho8rPY2Q3P7fLERERERERSTGtGU3r9uyBQQP5lvvpMOkefH3dLkhERERERCTlFEbTuENd+uHEnmXVI8OpW9ftakRERERERDxDYTQNi/9tOQVmTeSTHL3p9UE5t8sRERERERHxGK0ZTavi4zn4yPOcoxiF3u1PvnxuFyQiIiIiIuI5GhlNo459MJHCu1Yy7sahPNQ5t9vliIiIiIiIeJTCaFp07Bimbz+WcSsPfP0ojuN2QSIiIiIiIp6lMJoG/fPMQPKe3seqJ0dSpar+F4mIiIiISMajpJPGxG7YStGp7/FVzqd46kO1zxURERERkYxJYTSN+bttb2JMdvKOeoscOdyuRkREREREJHUojKYh+yd8T7nNs5le5VWCnyzqdjkiIiIiIiKpRmE0rTh7lthuPdnqVKRReHe3qxEREREREUlVCqNpxMZn3yfwxFbWPvEupStkdbscERERERGRVKUwmgaciz5BsfGDWZozmFajm7ldjoiIiIiISKpTGE0DVnf5iHzxh4h/ZQBZNSgqIiIiIiKZgMKoy84cOUWZ6cNYnqcxd7x4q9vliIiIiIiIeIXCqMuWdx5Dofj9+Lz6Co7jdjUiIiIiIiLeoTDqolOHT1Phm7dZm/cOgnrf4XY5IiIiIiIiXqMw6qKfO46jWHwUvgM0KioiIiIiIpmLwqhLjh44S+WZQ9iU91aq9bjb7XJERERERES8SmHUJYs7fk7J+L/J8voraFhUREREREQymyxuF5AZHdp3jmqzBrM9oA4Vut9L+JpIwuZtISo6huIB/oQEV6JVrUC3yxQREREREUk1Ghl1wQ9PT6Gs2UHW118mfG0UoTPWExkdgwEio2MInbGe8DWRbpcpIiIiIiKSahRGvWzP7jhqfT+IvwOqUarb/YTN20JMbFyi28TExhE2b4tLFYqIiIiIiKQ+hVEvm9NxOjeazWR9/WXw8SEqOibJ2yV3XEREREREJCNQGPWiXTvjufmHgewJuJGiXdsAUDzAP8nbJndcREREREQkI1AY9aLvOs+kGuvJ9vpL4OsLQEhwJfz9fBPdzt/Pl5DgSm6UKCIiIiIi4hXqpuslWzYbblv4JgfylqPQcw9dOJ7QNVfddEVEREREJDNRGPWSb7p8Tz9Wc2zAWMiS+NfeqlagwqeIiIiIiGQqmqbrBevWGu5c+iZH8pYmT9fH3S5HRERERETEdQqjXjDt2QXUYxlZX+0Hfn5ulyMiIiIiIuI6hdFUtmwZNF72JsfyBJKz61NulyMiIiIiIpImKIymsqldl3AnS8j28ouQLZvb5YiIiIiIiKQJCqOpaNEiaLb6TU7mLkK2bp3dLkdERERERCTN8EgYdRxngOM45pLLXk+cO70yBiZ1X0YT5pM1tA/4+7tdkoiIiIiISJrhya1dtgB3XfR9nAfPne7MmQOt/niTmJwF8H/+f26XIyIiIiIikqZ4MoyeM8Zk6tHQBPHx8GXvVUxiDnF9B0GuXG6XJCIiIiIikqZ4cs1oWcdxohzH+dNxnCmO45T14LnTlenTod3WgZzJEYBvj25ulyMiIiIiIpLmeCqMLgc6AE2BzkBR4FfHcQokdWPHcbo4jhPhOE7EgQMHPFRC2hAXB1/2/Z0HCMfvhR6QJ4/bJYmIiIiIiKQ5jjHG8yd1nFzATmCIMWb45W4bFBRkIiIiPF6DW6ZNg3PtHqJN9tn4Rf4F+fO7XZKIiIiIiIhrHMdZZYwJuvR4qmztYow5AWwAKqTG+dMqY2Diazt4kK/wfb6rgqiIiIiIiEgyUiWMOo6THbgR2JMa50+r5s2DxhtHYHyz4NOzu9vliIiIiIiIpFke6abrOM4w4Dvgb6Aw8AqQE5jgifOnFyNeP8I05zN45BEoXtztckRERERERNIsT42MlgAmY/canQGcAW41xvzlofOneUuXQrVlY8hpTvJ0vvqU6Teb+kMWEr4m0u3SRERERERE0hyPjIwaYx7yxHnSs6EDY/nEGcmvpWryk78dFY2MjiF0xnoAWtUKdLM8ERERERGRNCVV1oxmNqtXQ54fvibQRDK6bstE18XExhE2b4tLlYmIiIiIiKRNCqMeMPgtw4s+77A9X0l+KlvnP9dHRce4UJWIiIiIiEjapTCaQps3w/7pS6kZv5rpd7XFOP/9lRYP8HehMhERERERkbRLYTSFhgyBEJ93iC9QkMohz+Hv55voen8/X0KCK7lUnYiIiIiISNqkMJoCu3bB8s+30iz+O3y6Psf99cozuHU1AgP8cYDAAH8Gt66m5kUiIiIiIiKX8Eg33cxq2DDozgicrH7w3HOA7Zqr8CkiIiIiInJ5Ghm9Tnv3wvQxh3naGYfz2GNQpIjbJYmIiIiIiKQbGhm9Tu+9B0/FfkI2EwO9erldjoiIiIiISLqiMHodjhyBMaPOsj3b+3DHPVC1qtsliYiIiIiIpCsKo9dh1ChodmIq+dgDvce5XY6IiIiIiEi6ozB6jU6cgPfeNSzPMxxK3gT33ON2SSIiIiIiIumOwug1GjMGqh9eRDnWQu+x4DhulyQiIiIiIpLuKIxegzNn7HYu0/IPhyyF4ZFH3C5JREREREQkXdLWLtdg4kTIHbWZeodnQ9eukD272yWJiIiIiIikSwqjV+ncORg6FN4q9B4mWzZ49lm3SxIREREREUm3NE33Kn39NUTvOEjLrBNwnnwCChVyuyQREREREZF0SyOjVyE+Ht56C14t9DG+Z09Dr15ulyQiIiIiIpKuKYxehVmzYOsfZ+h89gNo1gwqV3a7JBERERERkXRNYfQKjIFBg+D5ApPxP7oPevd2uyQREREREZF0T2tGr2DRIlixwjAncDhUrw6NGrldkoiIiIiISLqnMHoFb70FD+abT4HI9TBoPDiO2yWJiIiIiIike5qmexlvT9zPggXQ2XmDg7nzM/PGBm6XJCIiIiIikiEojCYjfE0kn65bzy21ZtH48M+Mq9mcvrO2Er4m0u3SRERERERE0j2F0WSEzdvCWb/T9Ih7n5gs2fiy1r3ExMYRNm+L26WJiIiIiIikewqjyYiKjqHAyWhab1jEtGp3E+2f58JxERERERERSRmF0WQUD/Cn2PGD/BVQjM+CWiY6LiIiIiIiIimjbrrJCAmuROjJs9zTcdSFDrr+fr6EBFdyuTIREREREZH0T2E0Ga1qBQJ27WhUdAzFA/wJCa504biIiIiIiIhcP4XRy2hVK1DhU0REREREJBVozaiIiIiIiIh4ncKoiIiIiIiIeJ3CqIiIiIiIiHidwqiIiIiIiIh4ncKoiIiIiIiIeJ3CqIiIiIiIiHidwqiIiIiIiIh4ncKoiIiIiIiIeJ3CqIiIiIiIiHidwqiIiIiIiIh4ncKoiIiIiIiIeJ3CqIiIiIiIiHidwqiIiIiIiIh4ncKoiIiIiIiIeJ3CqIiIiIiIiHidwqiIiIiIiIh4nWOMcbcAxzkA/OVqEVdWEDjodhGS6el5KGmBnoeSVui5KGmBnoeSFqSH52FpY0yhSw+6HkbTA8dxIowxQW7XIZmbnoeSFuh5KGmFnouSFuh5KGlBen4eapquiIiIiIiIeJ3CqIiIiIiIiHidwujVGe12ASLoeShpg56HklbouShpgZ6Hkhak2+eh1oyKiIiIiIiI12lkVERERERERLxOYVRERERERES8TmH0MhzHec5xnD8dxzntOM4qx3EauF2TZB6O4wxwHMdcctnrdl2S8TmOc4fjODMdx4k8/7zrcMn1zvnnZ5TjODGO4yx2HOcml8qVDOoqnofjk3iNXOZSuZJBOY4T6jjOSsdxjjmOc8BxnO8cx6l6yW30miip6iqfh+nyNVFhNBmO47QHRgBvAbWAX4HvHccp5WphktlsAYpddKnmbjmSSeQC/gB6ADFJXP8i8ALwPFAX2A/86DhObq9VKJnBlZ6HAPNJ/BrZzDulSSZyF/AhcBvQCDgHzHccJ/9Ft9FroqS2u7jy8xDS4WuiGhglw3Gc5cDvxpjOFx3bBkwzxoS6V5lkFo7jDADaGmOqXum2IqnFcZwTQDdjzPjz3ztAFPCBMWbQ+WP+2DdffYwxn7hVq2Rclz4Pzx8bDxQ0xtznVl2S+TiOkws4CrQyxnyn10Rxw6XPw/PHxpMOXxM1MpoEx3GyAnWAHy656gfsJxIi3lL2/LSfPx3HmeI4Tlm3C5JMrwxQlIteH40xMcAS9Poo3ne74zj7HcfZ6jjOGMdxCrtdkGR4ubHvn4+c/16vieKGS5+HCdLda6LCaNIKAr7AvkuO78O+4Ih4w3KgA9AU6Ix97v3qOE4BN4uSTC/hNVCvj+K2ucATwN3YKZI3Awsdx8nmalWS0Y0A1gK/nf9er4nihkufh5BOXxOzuF2AiCTNGPP9xd+fX4S+E3gSGO5KUSIiaYQxZspF3653HGcV8BfQHJjhTlWSkTmOMxy4HbjdGBPndj2SOSX3PEyvr4kaGU3aQSAOKHLJ8SKAupmKK4wxJ4ANQAW3a5FMLeE1UK+PkqYYY6KA3eg1UlKB4zjvAg8DjYwxOy+6Sq+J4jWXeR7+R3p5TVQYTYIx5iywCmhyyVVNsF11RbzOcZzswI3AHrdrkUztT+wbrAuvj+efmw3Q66O4yHGcgkAgeo0UD3McZwT/BoDNl1yt10Txiis8D5O6fbp4TdQ03eQNBz53HGcF8AvwP6A48LGrVUmm4TjOMOA74G+gMPAKkBOY4GZdkvGd79JX/vy3PkApx3FqAoeNMX87jvMe0N9xnM3AVuBl4AQwyYVyJYO63PPw/GUAMB37RusGYDC2g+k3Xi5VMjDHcUYBjwOtgCOO4ySsAz1hjDlhjDF6TZTUdqXn4fnXywGkw9dEbe1yGY7jPIfdO6oYdq+zXsaYJe5WJZmF4zhTgDuwDbUOAMuAV4wxG10tTDI8x3HuAhYlcdUEY0yH81sZvAY8A+TDNtvqaoz5w2tFSoZ3uech8CwQjt0HPAD75msR9jXyH68UKJmC4zjJvVF+3Rgz4Pxt9JooqepKz8Pz2wmFkw5fExVGRURERERExOu0ZlRERERERES8TmFUREREREREvE5hVERERERERLxOYVRERERERES8TmFUREREREREvE5hVERERERERLxOYVRERERERES8TmFUREREREREvE5hVERERERERLzu/3TmTBg7djYEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x1, y, \"o\", label=\"Data\")\n",
    "ax.plot(x1, y_true, \"b-\", label=\"True\")\n",
    "ax.plot(np.hstack((x1, x1n)), np.hstack((ypred, ynewpred)), \"r\", label=\"OLS prediction\")\n",
    "ax.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting with Formulas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using formulas can make both estimation and prediction a lot easier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:22.841711Z",
     "iopub.status.busy": "2022-02-08T18:18:22.840138Z",
     "iopub.status.idle": "2022-02-08T18:18:22.856433Z",
     "shell.execute_reply": "2022-02-08T18:18:22.857122Z"
    }
   },
   "outputs": [],
   "source": [
    "from statsmodels.formula.api import ols\n",
    "\n",
    "data = {\"x1\": x1, \"y\": y}\n",
    "\n",
    "res = ols(\"y ~ x1 + np.sin(x1) + I((x1-5)**2)\", data=data).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the `I` to indicate use of the Identity transform. Ie., we do not want any expansion magic from using `**2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:22.861189Z",
     "iopub.status.busy": "2022-02-08T18:18:22.859753Z",
     "iopub.status.idle": "2022-02-08T18:18:22.871701Z",
     "shell.execute_reply": "2022-02-08T18:18:22.872368Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept           4.944539\n",
       "x1                  0.505053\n",
       "np.sin(x1)          0.578604\n",
       "I((x1 - 5) ** 2)   -0.020098\n",
       "dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we only have to pass the single variable and we get the transformed right-hand side variables automatically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T18:18:22.876494Z",
     "iopub.status.busy": "2022-02-08T18:18:22.874998Z",
     "iopub.status.idle": "2022-02-08T18:18:22.894041Z",
     "shell.execute_reply": "2022-02-08T18:18:22.894403Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    11.046441\n",
       "1    10.889753\n",
       "2    10.604492\n",
       "3    10.242369\n",
       "4     9.871452\n",
       "5     9.559500\n",
       "6     9.357376\n",
       "7     9.286595\n",
       "8     9.334052\n",
       "9     9.455232\n",
       "dtype: float64"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.predict(exog=dict(x1=x1n))"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}