{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prediction (out of sample)"
   ]
  },
  {
   "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 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": {},
   "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.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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.984\n",
      "Model:                            OLS   Adj. R-squared:                  0.983\n",
      "Method:                 Least Squares   F-statistic:                     970.0\n",
      "Date:                Thu, 05 Nov 2020   Prob (F-statistic):           1.43e-41\n",
      "Time:                        07:28:38   Log-Likelihood:                 1.5386\n",
      "No. Observations:                  50   AIC:                             4.923\n",
      "Df Residuals:                      46   BIC:                             12.57\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.9649      0.083     59.546      0.000       4.797       5.133\n",
      "x1             0.4954      0.013     38.528      0.000       0.470       0.521\n",
      "x2             0.4615      0.051      9.129      0.000       0.360       0.563\n",
      "x3            -0.0190      0.001    -16.796      0.000      -0.021      -0.017\n",
      "==============================================================================\n",
      "Omnibus:                        4.447   Durbin-Watson:                   2.049\n",
      "Prob(Omnibus):                  0.108   Jarque-Bera (JB):                3.703\n",
      "Skew:                           0.662   Prob(JB):                        0.157\n",
      "Kurtosis:                       3.163   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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.49083841  4.95047481  5.3736972   5.73535453  6.01937263  6.22139516\n",
      "  6.34949933  6.42286885  6.46864209  6.51745333  6.59839963  6.73426028\n",
      "  6.93775419  7.20945026  7.53767386  7.900425    8.26899278  8.61267261\n",
      "  8.9038108   9.12234804  9.25911506  9.31733912  9.31211379  9.26791911\n",
      "  9.21459959  9.18246084  9.19729045  9.27612116  9.42443242  9.63525015\n",
      "  9.89029254 10.16297372 10.42277173 10.64024318 10.79186077 10.86387868\n",
      " 10.85459102 10.77461212 10.64513255 10.49443765 10.35326029 10.24973155\n",
      " 10.20475886 10.22859099 10.31913416 10.46229611 10.63430134 10.80559744\n",
      " 10.94571039 11.0282514 ]\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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[11.02549655 10.90062849 10.67174529 10.38009018 10.07995374  9.8253817\n",
      "  9.65694264  9.59179541  9.61948802  9.70451667]\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": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6MAAAHWCAYAAACCFl9HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABz0UlEQVR4nO3de3zO5R/H8de1GTbDnNlQTq2EHCaFQqlVKHTSkRSRU2GiUkk551DSUengJxUmlAkh5DBRchg5dJhzmeMw2/f3x2UYw3Dv/t7b3s/H437Mvvd9f+/P9vt2/+73ruv6XMZxHERERERERES8yc/tAkRERERERCTnURgVERERERERr1MYFREREREREa9TGBURERERERGvUxgVERERERERr1MYFREREREREa/L5XYBRYsWda688kq3yxAREREREZFMsGLFij2O4xQ787jrYfTKK68kNjbW7TJEREREREQkExhj/kzvuKbpioiIiIiIiNcpjIqIiIiIiIjXKYyKiIiIiIiI1ymMioiIiIiIiNcpjIqIiIiIiIjXud5N90L279/Prl27SEpKcrsUyUYCAgIoXrw4BQoUcLsUEREREZEcyafD6P79+9m5cydhYWEEBgZijHG7JMkGHMchMTGR+Ph4AAVSEREREREX+PQ03V27dhEWFkZQUJCCqHiMMYagoCDCwsLYtWuX2+WIiIiIiORIPh1Gk5KSCAwMdLsMyaYCAwM1/VtERERExCU+HUYBjYhKptG1JSIiIiLiHp8PoyIiIiIiIpL9KIyKiIiIiIiI1ymMZoI2bdpgjMEYc3ILkUaNGvHOO+9c1BrFefPmYYxhz549mVitiIiIiIiI9ymMZpLGjRuzfft2tm7dyqxZs2jWrBmvvPIKN910E4cOHXK7PBEREREREVfliDAavTKeeoPmUq73DOoNmkv0yvhMf808efJQsmRJwsLCqF69Ot27d2fevHn88ssvDBkyBIAvvviC2rVrkz9/fooXL879999/cu/LrVu30qhRIwCKFSuGMYY2bdoAMHPmTG666SYKFSpE4cKFiYyMZN26dZn+M4mIiIiIiHhKtg+j0Svj6TN5NfEJiThAfEIifSav9kogPVOVKlW44447mDRpEgDHjh2jX79+/Prrr0yfPp09e/bw0EMPAVCmTJmTj1uzZg3bt29n1KhRABw6dIhnn32WZcuWMW/ePAoWLEizZs04duyY138mERERkZzGjYEOkewol9sFZLahMXEkJiWnOZaYlMzQmDia1wjzej2VK1dm9uzZALRt2/bk8fLly/Puu+9yzTXX8M8//1C6dGkKFy4MQPHixSlatOjJx957771pzvnJJ59QoEABli1bRv369b3wU4iIiIjkTKkDHamfL1MHOgBXPluKZGXZfmR0W0LiRR3PbI7jnNzf8pdffuGee+7hiiuuIH/+/ERERADw119/nfccmzZt4uGHH6ZChQoUKFCAEiVKkJKScsHniYiIiMjlOd9Ah4hcnGwfRkNDAi/qeGZbu3Yt5cuX59ChQ0RGRhIUFMTnn3/O8uXLmTlzJsAFp9s2a9aM3bt38/7777N06VJWrlxJrly5NE1XREREJJP52kCHSFaW7cNoVGQ4gQH+aY4FBvgTFRnu9Vp+//13Zs6cyX333cf69evZs2cPAwYM4Oabb+bqq69m165daR6fO3duAJKTT/317d9//2XdunW88MILNG7cmGuuuYYDBw5w/Phxr/4sIiIiIjmRrw10iGRl2T6MNq8RxsCWVQkLCcQAYSGBDGxZNdPn9B89epQdO3awbds2fv31V4YPH07Dhg2pVasWPXv2pGzZsuTJk4fRo0ezefNmZsyYQd++fdOc44orrsAYw4wZM9i9ezcHDx6kUKFCFC1alA8//JA//viD+fPn06FDB3LlyvbLf0VERERc50sDHSJZXY5IMM1rhHl9Qfns2bMpVaoU/v7+hISEUKVKFV555RWefvppcufOTb58+fj000954YUXeOedd6hWrRrDhw/njjvuOHmOsLAw+vXrx4svvshTTz3F448/zrhx45g4cSJdu3alSpUqVKxYkTfffPOspkYiIiIi4nmpnymHxsSxLSGR0JBAoiLD1bxI5BIYx3FcLSAiIsKJjY1N975169ZxzTXXeLkiyUl0jYmIiIiIZC5jzArHcSLOPJ7tp+mKiIiIiIiI71EYFREREREREa9TGBURERERERGvUxgVERERERERr1MYFREREREREa/LEVu7iIiIiEj2Fr0yXtutiGQxCqMiIiIikqVFr4ynz+TVJCYlAxCfkEifyasBFEhFfFiGpukaY242xnxrjIk3xjjGmDZn3N/SGBNjjNl94v6GmVCriIiIiMhZhsbEnQyiqRKTkhkaE+dSRSKSERldMxoM/A50AxLTuT8fsBjo7qG6REREREQyZFtCeh9Pz31cRHxDhqbpOo7zHfAdgDFmXDr3f37ivqKeLE5ERERE5EJCQwKJTyd4hoYEevaFHAd+/RUmTICvv4bdu8HPD/z97dfT/52RYxUrwjPPQMOGYIxnaxXJArRm1MPMBd5IWrduzbhx47xTjIiIiEgWdLHNiKIiw9OsGQUIDPAnKjLcMwVt3GgD6IQJsH495MoFt90GLVpASgokJ5//a3rHjh+HH3+Eb76BKlWga1d45BEICvJMzSJZgCth1BjTHmgPULZsWTdKyDTbt28/+e/p06fTrl27NMcCA9P+hS4pKYmAgACv1SciIiLiyy6lGVHqcY92042Ph4kTbQCNjbXHbr4ZunWD++6Doh6YEJiYaM//1lvQvj08/zy0awedOkE2+4wskh5X9hl1HOcDx3EiHMeJKFasmBslZJqSJUuevIWEhKQ5duTIEUJCQpgwYQK33HILgYGBvP/++4wbN47g4OA055k3bx7GGPbs2XPy2OLFi2nQoAFBQUGEhYXRsWNH9u/f780fT0RERCRTXWozouY1wljU+xa2DGrCot63XFoQ/fdfeP99O222TBno0cNOzR02DP7+G+bPhw4doGhRHAeOHLn4l0gjMBDatoWVK+25b7nFvla5cjbwLlhgX18km3IljOZ0ffr04ZlnnmHt2rU0b948Q89ZvXo1t99+O3fffTe//vorkydPZtWqVbRt2zZzixURERHxIq83I0pOhi+/hKZNoWRJGzZ37IBXXoG4ODsq2qMHlC59cslonz5QoQL06+ehGoyxo67ffAObN0NUFMydCw0aQM2a8MknHki+Ir4ny60ZffZZWLXKu69ZvTqMHOm583Xp0oX77rvvop4zdOhQHnzwQXr06HHy2LvvvkuNGjXYtWsXxYsX91yBIiIiIi7xWjMigIULoUsX++GydGn7QfPhh+2Hv9P6gKxbZ2fsfvmlzaf+/tC4MVx/vedL4oorYNAgePllGD/eTuFt2xZ69bJTeTt2tLWKZAMZCqPGmGCg4olv/YCyxpjqwH+O4/xljCkMlAVCTjymojEmAdjhOM4Oj1acDURERFz0c1asWMEff/zBxIkTTx5zTkzb2LRpk8KoiIj4hIMHbU+W1auhSBEoUQKKFz/1NTjYc01DL7bJjWQNmd6MCGD7dhvuvvjCBrsvv4T777ddbk/YvPlUAP3tN3vdNmwIzz0H997rmSWj5xUUZNePPvUUzJtnQ+nAgTB4sC2ga1eoW1ddeCVLy+jIaATw42nf9ztx+xRoA9wNfHLa/R+e9rhXL6vCM3hyhNIt+fLlS/O9n5/fyWCZKikpKc33KSkpPPXUUzz33HNnnS8sTP/HKyIi7nAcWLsWZs6E77+Hn36CY8fO/fjAwLTh9PSvpUrB7bfDiZYL53UpTW4ka8iUZkSpjh2zoa5fP/vvF16wtxOfzf7+G776yobQ5cvtU+rWtU+57z57jXqdMdCokb1t2QJjxsBHH9lC69SBzz+HSpVcKEzk8mV0n9F5wDn/7OI4zjhgnEcqyoGKFSvG4cOH2b9/PwUKFABg1RlzkWvWrMmaNWuoWLFiOmcQERHxnv37Yc4cGz5nzrQf4AGuvdYO1txxh/2MvG8f7NoFO3em//Wff2DFCvvv5BODYEFBdpZkx452qdy5nK/JjcJo1te8Rpjn/3ecNcteoHFx0KSJHeGoWJGjR2HsGPjf/2DRIvvQiAgYOhQeeMDHmtqWK2cLe/VVG0JffBFq17YjvE2bul2dyEXLcmtGs6M6deqQL18++vTpw3PPPcevv/7KmDFj0jzm+eef54YbbqBDhw48/fTT5M+fn/Xr1zNt2jTef/99lyoXEZGcwHHsNMXU8Llokd0isUABu26ub18bQMuUOg6LF8O0afDGCoLz5CEsKMgmzNNv4UFQ49T3KXmDOJgSxF978/PO0gg++19ePvrIrsfr2BEefNCOqJ7O601uJOvauhW6d4cpU2zXoenTbRgFYmKgc2f44w+oWhXeeMMGUJ//23++fLbR0h132Cm7zZrZNaavvJJmqrGIr1MY9QGFCxdm/PjxREVF8fHHH3PzzTfTv39/HnvssZOPqVatGgsWLOCll16iQYMGJCcnU758eVq0aOFi5SIikp0lJsKAATB2rF1iB7avS8+ecOedcOONEJC43444vfgtzJgB//0HAQF2WPPQITv8efjwqduhQ6eGQU/wAwoAVYB38+dnVNN7mFnwQfouuI0nnshD9+7wxBP2s3fqbESvNrmRrCkxEYYMsc2A/Pxs0uzeHfLm5e+/7drPSZPsNRUTY6eIZzlXXmmbMHXsCK+9Zjv/fvEFFCrkdmUiGWLOXKvobREREU5s6kbCZ1i3bh3XXHONlyuSnETXmIhI+n76yfZN2bAB7rkHmjeHyMgTa+b++suOfn77re1WlJQEhQvb0aa777af6k8sO0lXUlLagJp627kTpk61I1h79+IULMiOG1sw7tAD9F98K4nJuWnc2H7udsrE89K3Zze5Gdiyqqbp5nSOY6+j556zo6IPPGD37ixThmPH7Ozc116DlBR46SW7a0uePG4XfZkcB957D7p1s/ujTpkC1aq5XZXIScaYFY7jnNXFVWFUcjRdYyIiae3fD717w7vv2uVpH3wAjW9JgV9+ORVAU/saVKpkk+rdd9th0lwemnB17BjMnm0btERHw759pIQU4rcKLRn21wN8ufsWSoTm4qam+/mj8G/86+xTN12xNm+GZ56xQ53XXgtvv20b/2D/btKpk92m5Z57bCi98kpXq/W8xYttp6WEBNvk6OGH3a5IBDh3GNU0XREREQHsLNsOHWDbNjuo1P/FI+QbPRhaf2AP+vlBvXp26uPdd0O4B7faOF3u3HDXXfZ29CjMmoXfV19RfepXfHFgLGMLFGWOf0ve/OBBfvNrQOu2/gx5XjMTc7zx4+2wuTE2aT7zDAQEsH27Hf2cMMH+geW0JaPZT9269g9HDzwAjzxiWwIPGWKnzov4IK1wFhERyeF277YDKE2bQsGCdnBleIufyFevuu3aWbMmfPqpnUa7YAFERZ0ziEavjKfeoLmU6z2DeoPmEr0y/vKKy5PHNmf5/HPbdnfKFPLc1Zi7/hvPHG7l3zyhlB37ChHhB/jqKztbUXKYAwegdWt49FHbhejXX6FbN46bAEaOtJfq5Mm2t8+aNdk4iKYqWdK2u+7a1Ybyxo3tf7siPkhhVEREJIdyHDuYdM018M03duvFX37cR51xHeHmm+HIEds+d9o0ePxxKFr0vOdL3fszPiERh1N7f2YkkGYoxObNaxevTphgg+k335D/9hvp67zGsr0VmffgGFo0TTq51YzkAMuX2z+WfPGF7SY7fz5ceSULF9rDzz0H9evD77/bv6uc2ZU52woIgFGj7O8l9Xe0ZInbVYmcRWFUREQkB/r7bzsS+uijdunnypXw8nVTyV29sl0o+txz9hN8ZGSGz3m+vT/P55JCbFCQ3dIiOhqWLqXwjVczhk4M+b4Kva+azNtvOWc27ZXsJCXFTj+tW9dO5Z43D/r14999uWjTBm66ye5zO2WKnX7u81u1ZJZHHoGff7YzDG6+2TY50vQB8SEKoyIiIjlISgqMGQOVK9vP76NGwcKvt3Ptq/fbUceiRe0IyvDhEBx8Uee+1L0/LzXEnnT99Zj582DaNK6s4M/4I/dSq1s9OlZbxOrVGTuFZCHbt9s/kjz/vO1E9OuvcNNNLFxotx763/+gTx9Yu9Ze0sa4XbDLrrvObvly6612Te2TT9pZDyI+QGFUREQkh4iLgwYNbEfRunVhze8OXYM+wr/KNXYq7oAB9kNr7drAxa//PNcenxfa+/NSQ2waxkDTpuRe9xvOBx9SPWQrH6ytzx/XtWRkx7iTn709vqZVvGv6dLtlyaJFdgT/669JKViIQYOgYUM7ALhkib2U8+XL2ClzxDVRuLD93fXtC598Yucu//mn21WJKIyKiIjkBHPnQq1atoHLuHEw8+2NXNn2FmjXzg4n/fabHU460XXzUqbORkWGExjgn+ZYYIA/UZHn77p7qSE2XblyYdo9RdA/Gzncpz93+P9A5/euZXLJZxgxfM0lr2kVlx05YhvyNGsGYWGwYgW0a8fuPYYmTeyle999tpFszZoZP+3lrHPOcvz97Qar334LGzfaN4TZs92uSnI4hVEREZFsbvp0u0tKuXLw+8okWm8fhKlW1S4U/eADm1SvuirNcy5l6mzzGmEMbFmVsJBADBAWEsjAllUvuPfnpYbY88qXj6ABLxEYv4ltd3fg/n0f0q5HHR79cBZ59idl+GcSH7BuHdSpY/cM7drVDn1ecw0LFti/o/z4o10KOWECFChwcae+7CniWVGzZrapUcmSdrrzp5+6XZHkYAqjctk6d+5Mw4YNT37fpk0bmjZtelnnfPXVV6lSpcplViYiIhMnQosWdseLn0auIPSe2nYYqUkTu6iuXTu7f+gZLnXqbPMaYSzqfQtbBjVhUe9bLhhEU59zKSE2Q4oXp+zU0RxftYY5herz4s4RzH+vAy1ilmCclAz9TOISx7F/LKlVy+5zO306jBpFckBeXn8dGjWyy5qXLoWnn760taEemSKeFV11lQ31t94KTzwBn33mdkWSQymMZpL4+Hjat29P6dKlyZ07N2FhYbRr145//vknzeMuFNx+/fVX7rnnHkqWLEnevHkpW7Ys9957L3/68Dz/UaNG8cUXX2TosVu3bsUYQ2xsbJrjPXv2ZP78+ZlRnohIjjF2LDz0ENx4I8zv8g0hd9W1W6JMmmRvoaHnfK5Hp85mwKWE2IsReN1VDOnVi2Z3jWFTrisZsep1vnjndcrtjM+0n0kuw969cP/9NmXWr2+nkTdpws6ddjCvb197bcfG2v48l8rb17lPCQ6GqVNtIG3Txu7lK+JlCqOZYMuWLURERPD777/z6aef8scff/DFF1+wZs0aateuzdatWzN0nt27d3PrrbcSHBzMjBkzWL9+PZ9//jkVKlRg//79Hq352LFjHjtXwYIFCQkJuaxzBAcHU6RIEc8UJCKSzWSk4crIkfDUU/aD+w8PfEBQmwfsCNPq1dCy5QVfI1OmzrosKjKcP2qW4+GufelUfgDXHNrIzHFdGLZyJnjw/wflMv30k02YU6fa7VtmzoRSpZg71x5etAg++shmp/z5L++lsuN1flECA+3vuVEjaN3a7ksq4kUKo5mgU6dO+Pn5MXv2bG699VbKli1Lo0aNmD17Nn5+fnTq1ClD51m0aBF79+7lk08+oVatWlx55ZU0aNCAIUOGULVq1XM+L3W09fXXX6dEiRIEBwfzxBNPkJh4aspJw4YN6dixIz179qRYsWLUq1cPgLVr19KkSRPy589P8eLFeeihh9ixY8fJ5yUnJ9OzZ08KFSpEoUKFePbZZ0k+YyO3M0d7HcfhzTffpFKlSuTJk4fSpUvTp08fAMqVKwdA7dq1McacnO575jTdlJQU+vfvT5kyZciTJw9Vq1Zl6tSpJ+9PHWGdNGkSt912G0FBQVSuXJkffvghQ79rEZGs4kINVxwH+ve324Te29Jh2o0DyNPlabjjDvjhB8jgH/oydeqsS1J/ptJFA/nu/mo0f/grZuRuzo1fDeW/K2vgLFrsdok52/Hj0K+fbYubOzcsXgxRUSQ7frzyCjRuDIUK2eWOTz7pmS1bsuN1ftGCgmw37dRAOn682xVJTuI4jqu3WrVqOeeydu3ac97nq/7991/HGOO88cYb6d7/+uuvO8YY57///nMcx3Fat27tNGnSJN3H/vzzzw7gjB8/3klJSclwDa1bt3aCg4Od++67z1m9erUzc+ZMJzQ01OnSpcvJxzRo0MAJDg52unfv7qxbt85Zu3ats23bNqdIkSJOr169nLVr1zq//vqr07RpU6d27dpOcnKy4ziOM3jwYKdAgQLOxIkTnXXr1jmdO3d28ufP7zRo0CDN65/+M/Xu3dspWLCgM3bsWGfjxo3O4sWLnXfeecdxHMdZtmyZAzgzZ850tm/f7vz777+O4zjOK6+84lx77bUnzzF8+HAnf/78zvjx4524uDinb9++jp+fn7Ny5UrHcRxny5YtDuCEh4c73377rbNhwwbn8ccfdwoXLuwcOHDgnL+rrHiNiUjOVnfgHOeK56efdas7cI6TkuI4PXs6DjhO68eSneSuz9pvHn7YcY4dc7t0n7Rjh+P0rTnd+ZMyTjLGSXzyGcfZt8/tsnKev/5ynJtustfrY485zv79juM4zrZtjtOw4YlrurXjHDzobpnZ2qFDjtOokeP4+TnO+PFuVyPZDBDrpJMFc7kbhS/Bs8/CqlXefc3q1e18pwzYuHEjjuNwzTXXpHt/5cqVcRyHjRs3cv3115/3XDfccAMvvPACrVu3plOnTtSuXZuGDRvyyCOPcMUVV5z3uf7+/nzyyScEBwdTpUoVBg8ezJNPPsnAgQPJd2LjrXLlyvHmm2+efM7LL7/Mddddx+DBg08e++yzzyhcuDCxsbFcf/31jBw5kl69evHAAw8Adn1oTEzMOes4ePAgI0aMYOTIkbRt2xaAihUrcuONNwJQrFgxAIoUKULJkiXPeZ5hw4bRs2dPHn74YQBee+01FixYwLBhw9KsT33uuedo1qwZAAMGDOCzzz5j1apV1K9f/7y/LxGRrOJcjVXi9ybSsSO8/z506ZDEqANtMW99YbuPjhiRbpMigRIl4NXlTXh7wBr8Xn6JTmPf5ui3U8nz4Ttwzz1ul5czTJ5s55QnJdm5t48+CtiB/EcfhYMH7daYbdq4W2a2lzpC2rQpPPaYHXp+6CG3q5JsTv/PlEnMOeaO2D8MnPv+M73xxhvs2LGDDz74gKpVqzJ27FgqV67MnDlzzvu8atWqERwcfPL7G2+8kWPHjrFp06aTx2rVqpXmOStWrGDBggUEBwefvJUpUwaATZs2sW/fPrZv334ySAL4+flRp06dc9axdu1ajh49yq233pqhnzc9+/fvZ9u2bSenEqeqX78+a9euTXOsWrVqJ/8deqIxx65duy75tUVEfE16jVWcFMOhWbV4/33o2+Mwo/5ugRn/hZ2vO3KkgugF+PlBt5fyU2fJKO4L/Zm43YWheXOce++D7dvdLi/LO+ca58RE6NgR7r0XKlSwWw09+ijHj9sGRZGRUKyYnZarIOol+fLZrsU33WT/EjBxotsVSTaX9UZGMzhC6ZZKlSphjGHNmjU0b978rPvXrVuHMYYKFSpk+JxFihTh/vvv5/7772fgwIHUqFGD/v37X1bAA06OkKZKSUmhSZMmDBs27KzHlihRgpSUlIt+jdTw7QnpBfgzjwWc2Kz99PsupW4REV8VFRlOn8mrT+6N6Bz347/pNTkYV4Lhfffy3Jym8PPP8O670KGDy9VmLddfD5+srUOn9iso89Uw+kX3I9fs2fgNHWJH7hTqL1rqGufU6zV1jXP+P9Zza79usGYN9Opl/3CSOzfbtsHDD8P8+dC2rd1aNCjI5R8ip8mXD2bMsJsTP/KIHSE9MSNOxNP0ruphhQsXJjIykjFjxnD48OE09x0+fJh33nmHO++8k8KFC1/S+XPnzk2FChU4ePDgeR+3evVqDh06dPL7JUuWnHzuudSsWZM1a9ZwxRVXULFixTS3/PnzU7BgQUqVKsWSJUtOPsdxHJYtW3bOc1auXJk8efKccyQ3d+7cAGc1QTpdgQIFCA0NZeHChWmOL1y4kMqVK5/zeSIi2dHpDVecY/7sm1qHg3El+Pj1bTw35WY7jDRxooLoJSpYED7/MoBKY/tQK2A1ixJr2u1FGjaE9evdLi/LGRoTdzKIAuA43LtsGvUfvgv27IGYGBg8GHLnZtYsuzJq+XK77eXYsQqirkkNpHXr2r8OfP212xVJNqUwmglGjx7N8ePHady4MXPnzuXvv/9m3rx53HbbbTiOw+jRo9M8fv/+/axatSrNbevWrUyfPp1HH32U6dOns2HDBuLi4hg2bBjfffcdLVq0OG8Nx48fp23btqxZs4YffviB3r17065du7NGQ0/XqVMn9u3bx4MPPsjSpUvZvHkzs2fPpn379hw4cACAbt26MWTIEL755hvi4uJ49tln2X6eKUz58+enW7du9OnTh08++YRNmzaxbNky3n33XQCKFy9OYGAgMTEx7Ny5k3379qV7nqioKIYNG8aECRPYsGEDL7/8Mj/99BM9evQ47+9BRCQ7al4jjO863kKZZXdwYHNhvhn0B098VA+2bIHvvrP7M8olM8aOyn21shKdr55DW8ZyePnvONWqwQsvwGl/7JXzO32Nc0jift6f8gavzxrDz2Wq2r1Db7+d48fhpZdsw+cSJezeoY895mLRYgUH2/eTG2+0a0cVSCUTZL1pullAhQoViI2N5bXXXuOxxx5j165dFCtWjLvuuouJEydSunTpNI//6aefqFGjRppj9957L0OGDCE4OJiePXvy999/kytXLsqVK8ewYcPo1q3beWto0KAB1157LY0aNeLw4cMnz3c+oaGhLFq0iD59+nDHHXdw5MgRypYty+23306ePHkA6NGjBzt27OCpp54C4LHHHuORRx5h3bp15zzvwIEDKVSoEP379+eff/6hRIkSPP744wDkypWLt956i9dee41+/fpx0003MW/evLPO0bVrVw4cOECvXr3YuXMn4eHhTJo0ierVq5/3ZxIRyY4SEuw+9atXQ8yglTQedgckJ8OPP0Lt2m6Xl21ccw0sWWro2bMt5cY04eMiUTQZONA22Rk+HO67zzP7i2RjoSGBxCckUuev1YyY/iZFDyXQ/5aniLn1QRYWL058vM05P/1kZ0KPGqXRUJ+SGkjvvNP+D+XnZ9f4iniI8eSavksRERHhxMbGpnvfunXrztmVVs6tTZs27Nmzh+nTp7tdis/TNSYiWc3Ro3YEadEiWNB/PjcMuNvOLZ01C66+2u3ysq0pU2xYqn5wIV8W60yx+F/tXwTeftumVknXd3N+JbFHFPf++gObC4XS5e5ebC4TzsCWVQncFcajj9o+Ru+9d7KJrviiAwfsG8+yZXYZQMuWblckWYwxZoXjOBFnHtc0XRERkSwiJcVOH503D2Z3juaGVyIhLMwmUwXRTNWiBfz+O+RtXJ9S8bG8ddVoUpavgGrVICrKfliXU5KT4d13ueu+hrRY8yOfN2hF0zZvkXB1VV6/uyrLvwnjjjugZEk7LVdB1Mflzw/ff29nXjz4oP3rjIgHKIyKiIhkES++CP/7H3zzaDQ3j7oXrrvOzm88sQ2XZK5SpeyuF+9+kIsX4jtxlRPHxnqtYdgwCA+3/+O4POPMJyxbBnXqwDPPQM2a+P32G4/Nm8Da4ffy1WO3MKpHGAMG2JHmpUv1d5Qso0ABmDnTBtIHHoDoaLcrkmxAYTQbGjdunKboiohkM++9B4MGwbBm82n5dSv7gXDOHChSxO3SchRjoF0723un1HXFuWr+R/RuuISk4qF2G4yGDe1i3pzo339t5+EbboBt22DCBJg9++Q05pkzbbfcX36BL76ADz/U+tAsJzWQRkTYRmlTp7pdkWRxCqMiIiI+bvp06NQJOt/0K93n340pX95uuxAc7HZpOVb58na69JAhMGJxHa7YvpRfn3nfzuWtUQOefZbpC9ZRb9BcyvWeQb1Bc4leGe922ZkjJQU++siODo8dC889Z7fBadUKjCExEZ5/3vbAKVUKVqywuV2yqNRAWquWDaTTprldkWRhPh9G3W6wJNmXri0RyQqWL7dLtJpW3syouEhMgQJ2b0aNiLrO398uF42NheKl/Kk+pj3d7tzAsdbtcN56ixvuqssNP03HpCQTn5BIn8mrs18gXbkS6tWzw8WVK9vv33zTBhZsZqlSxYb29u3ttNzwcJdrlstXsKB9H6pe3XaVjolxuyLJonw6jAYEBJCYmHjhB4pcgsTERAICAtwuQ0TknLZsgaZNoXKRnUw6eDt+x5Psh75LWCMavTI+Z4zSuaBqVRuy+vSB0ROKED73XR5o8Tl/FSzBm9+NYO6HHWgT+y1+Bw8wNCbO7XI9IyEBunSx0zU3b4bPPoP58+0vA4iPt8sK77wTAgJg7lx4/30IDHS3bPGg1EBauTI0b263lhK5SD69tcv+/fvZuXMnYWFhBAYGYrSXl3iA4zgkJiYSHx9PiRIlKHDir7ciIr7k33/tgFPizv3ElWpI3j/j7BrRG2646HNFr4ynz+TVJCYlnzwWGODPwJZVaV4jzJNlZ0vRK+MZGhPHtoREQkMCiYoMP+fvbfFiePxx2LTJoWDEJh4q+TFProwmIn4d+3MH8VW123jqqxFQrpyXfwoPSUqC8ePtvNs9e2yTov79ISQEgOPH4Z13oG9f+9CXXoKePeHEduWSHe3eDY0a2b+ezZpl37hEznCurV18OoyCDaS7du0iKSnJi1VJdhcQEEDx4sUVREXEJx05Ao0bw+rlR/jz2rsIWf0TfPutHWa6BPUGzSU+4eyZRmEhgSzqfcvllputXUqQP3gQKjb+h51LSxNQ9ACFbl3DDbkX88SKb2m6fiG5cOCee6BbN7j5ZtsVydetWQOffAKffw67dtk/iowZY9fHnrBsGXToYGfq3nEHjB4NFSq4WLN4z44d0KABbN9um1Zdf73bFYmPOVcYzeVGMRejQIECCgwiIpJjpKTYkbWfFyXz5/WPELLsR9t69BKDKMC2dILo+Y7LKUNj4tIEUYDEpGSGxsSdM4wGB8N77xq6DI1l+/Rr2TXxBmLKVmJVo7rkHZ6LOxZMsXNWp0yxa+6efdY2+7nE4cOLGbm9KPv2wZdf2hC6dCnkygV33203u73zTvCzq70SEuCFF2zH51Kl4Ouv4d57s0bGFg8pWdLO3Lj5ZoiMtFN2q1d3uyrJAnx6zaiIiEhO06sXfP21w4obOlF62WQYPvyyW4+GhqS/UO9cx+WUSw3yzWuE8XZUKWpHLaVw4zWk/BfMX5/ewKiREcS2eAP++gs++ACOHYM2baBsWejXD3buvKj6Ukdu4xMSceDyGyWlpNgg8dhjNll26ACHDsGIEXa7lkmToEkT8PPDcezWqldfbbN1t26wbp3tZ6MgmgOVLm0XB+fPD7fdZkfTRS7A56fpioiIZFdnjmhV3luTj4aE8F3tl7lzeX/o3RsGDvTI62jN6KXx1BTnw4ftWsrBg+164ObN4bXXoGoVx05rHDXKbteTOzc89JBNdDVqQGjoeZPdpdZ35rX3crVgImNnwrhxdu1fwYLw8MN2FLRWrbNq2LDBLhedM8fOyHzvvTQzdiUn++MPO0KakgILFsBVV7ldkfiALLtmVEREJDs6MyAe3lCC3VNqMaTiCKL+6GFDwEcfeWyIKdOmcmZzng7y+/fb3DlsGBw4YLftefXVE9udxMXB22/babGHD9snFC1qU1716vZrjRpQqZLdVwYo13sG6X2SM8CWQU3O+TP1/3IZBfbupuqOP7hv9Wzqb12FHw7cequ99lq0SLf17a5ddi3o4MH27kGD7K4uJ8oRsdats2tIc+e2gbR8ebcrEpcpjIqIiPiQ00e0jsaHsPPLG3goeByf72uHX7NmdjpkLp9v7ZAjZEaQ/+8/ux3nqFGQmGhnxb788onP7IcOwapVthNQ6tfff7dTegGCgqBaNahencE7A1mUvwxxRa/gaEAeAo8dofih/7iWQ4xpVNI2lNm27dRt+3YObf2bfEcPn6zlnwLF+abqrfxUrymTBj98Vq3HjsH06XbQ9PvvbcfcRx6x9ZcocVm/BsnOfvvNdtktUMAG0kvYkkqyD4VRERERH5I6opW0N4gdn9flNr9ZfJvYklWh4Vy/Ybk2ZMwhdu2yo4xjxtiQ9+STdjuU0qXPeGBSkh1tOj2grlplmwwBx40fhwPyUuDY4TNfAvLmtes/Q0MhNJRPNh9hZ3ARdgYX5q+QkvwSdjWO8Uszmuo49iXGjbPrQv/9157iscegdWu7taTIBa1YAbfcAsWL20BaqpTbFYlLFEZFRER8SL1Bc/l7RxLbP69HzYO/MTflVv4qVJJnO44ipt/dbpcnXrZtGwwYYHsa+flB06Z2lmODBlClysnGtWk5DmzdytJJs1n33QLYv4/EoiW4vl4Vat1Q+WT4JCQkzXTv860znfzELYwfb0Po6tW2wW/z5jaA3nabBuvlEvz8s714ypaFefNsMJUcR2FURETEh0xaHk/rVnkI3bKHn3PfyOG8eXikzXC6t26gtZw52J9/2pHS77+HrVvtscKF4aabToXT6667vDWaZ66DdY77cXxrSUJ3Xc3KxYEkJ0OdOrbJ74MPQqFCl/1jSU43f77dDqhSJdutuXBhtysSL8uy+4yKiIhkR/M/DyNo8y5m57sdUlLo+dQQuj+iIJrTXXGFnbILNpjOn3/qNnWqPV6wYNpwWqPG+UcsU1Js46S9e+1a1eC9YTQNDmTykl3s/isvievCOJ4YgFMKeva0o6DXXJP5P6vkIA0awLff2iH/22+3bZgLFnS7KvEBGhkVERHxsnffhe7PJLK+ZCOuSPjVjhTccIPbZYmP++eftOF0wwZ7PH9+qFfPNj/au/dU6Pzvv1Pfp6Skf87Uabht2kDjxpqGK5lsxgzbqTkiAmJi7MUrOYKm6YqIiPiAOXPgjttTmFfsfurumoL55hto2dLtsiQL2r7d9oSZP98uxdu5006pLVzY3jLy7yJF7O4bIl4zeTI88ADUrw/ffWe7Q0u2pzAqIiLisg0b7ADocL8etPl3OAwfDs8953ZZIiLeNWGC3R/ottvs/PO8ed2uSDKZ1oyKiEiOlBl7RF6KvXuhWTNof2w0bQ4Nhy5d4NlnvV6H+DZfuV5FMtVDD8HRo/DEE3ae+DffQHCw21WJC9JrFH4WY8zNxphvjTHxxhjHGNPmjPuNMeZVY8w2Y0yiMWaeMebaTKlYREQkg1K7hsYnJOIA8QmJ9Jm8muiV8V6tIykJ7r8fKm+axsDEbjaVjhiRZrsNEV+5XkW8ok0bGDsWZs+2DY527HC7InFBhsIoEAz8DnQDzt6YCnoBPYAuQG1gF/CDMUarkkVExDVDY+JObl+RKjEpmaExcV6t49lnIWFOLF/5t8LUrGmnqF3O3hySLfnK9SriNW3b2i67cXF2DcO6dW5XJF6WoTDqOM53juO84DjON0CafmzGGAM8CwxyHGeS4zi/A62B/MDDHq5XREQkw7YlpPf303MfzwzvvAPTx/zJj0FNCShVDKZNg3z5vPb6knX4wvUq4nV33WW7cB05AnXr2q5ckmNkdGT0fMoBJYFZqQccx0kEFgB1PXB+ERGRSxIaEnhRxz3thx/g5a4J/BR8F8EBR2znyJIlvfLakvW4fb2KuKZWLViyxL4/3nYbfPml2xWJl3gijKb+v+rOM47vPO2+NIwx7Y0xscaY2N27d3ugBBERkbNFRYYTGJB2OmxggD9RkeGZ/tpxcfDwfcf4LrAlZY5uxEyZApUrZ/rrStbl5vUq4rorr4TFi+103YcegiFDwOVdPyTzeSKMpjrzajHpHLMPdJwPHMeJcBwnolixYh4sQURE5JTmNcIY2LIqYSGBGCAsJJCBLatmenfS//6DZk0dRh9tR51DP2LGjoVGjTL1NSXrc+t6FfEZhQrBrFnw4IPw/PPQqRMkJ1/4eZJleWJrl9TWVyWBv087XpyzR0tFRES8qnmNMK9+mE/tnPvY5n48mPIZvPYaPPaY115fsjZvX68iPidPHvjf/+CKK+zo6D//2KZvWmufLXliZHQLNpDelnrAGJMXuAlY7IHzi4iIZAmOY7cPLT33U/qm9LNbF7z0kttliYhkLX5+MHgwjB4NM2bYmSU7NcaVHWV0n9FgY0x1Y0z1E88pe+L7so7jOMBIoLcxpqUxpgowDjgI/C9zyhYREfE9b78NG9+fw8d+T8Gtt8L772svURGRS9WpE0yZAr//DjfeaBfjS7ZinAwsDDbGNAR+TOeuTx3HaXNie5dXgKeBQsBSoNOJbV7OKyIiwomNjb2YmkVEJAeKXhnP0Jg4tiUkEhoSSFRkuE9NZ/zmG3j1/jUszVWPoPDSmEWLoGBBt8sS8Shf/+9Qsqlly6BpU7t+9NtvoV49tyuSi2SMWeE4TsRZxzMSRjOTwqiIiFxI9Mp4+kxeTWLSqUYWgQH+PtPcZcECeOq2P1lIPYoWTsFv6RIoW9btskQ8ytf/O5RsbtMmuPNO+OsvGD8e7r3X7YrkIpwrjHqym66IiEimGBoTl+YDMEBiUjJDY9yfsrV6NbRtuotZzm0UDTyEX8xMBVHJlnz5v0PJASpUsFu/1Kplu8SNGOF2ReIBCqMiIuLztiUkXtRxb/n7b3ggch/RiZGU9f8Hv+9mQLVqrtYkkll89b9DyUGKFoXZs6FlS+jeHZ59Vlu/ZHEKoyIi4vNCQwIv6rg3/PcfNL/9MB/takZl1uAXPQXq1nWtHpHM5ov/HUoOFBgIX30Fzz0Ho0ZBixawbZvbVcklUhgVERGfFxUZTmCAf5pjgQH+REWGu1JPYiK0bJZE/7gHqJuyEL8vPofISFdqEfEWX/vvUHIwPz8YPtyG0Vmz4KqrYMAAOHLE7crkIimMioiIz2teI4yBLasSFhKIAcJCAl1rmpKcDI8+nEK7xW24y5mBefddePBBr9ch4m2+9N+hCABdu8LatXD77fDii3DttRAdbTd9lixB3XRFREQyyHGg0zMOld/rQmfegYEDoXdvt8sSEZE5c6BbN1izxu7zPHIkVKnidlVygrrpioiIXKYBA6Dke6/YIBoVBc8/73ZJIiICNoCuWgVvvw2//ALVq0OXLnaBv/gshVEREZEM+OQT2P3SSF6mP07bJ2HwYDDG7bJERCRVrlzQuTNs3AgdOsCYMVCpkv16/Ljb1Uk6FEZFREQuYMYMWPDkp4zkOVJa3Iv54H0FURERX1WkCIwebUdKr7sOOnWCmjXhxx/drixzZOGgrTAqIiJyHsuWwWcto/nQeZLjjRrjN2E8+Ptf+IkiIuKuqlXtWtJJk+DAAbjlFrjvPti61e3KLo/j2MZNQ4dCw4Z2SnIWpTAqIiJyDhs2wIDbfuSzYw+SUjOCXN9OgTx53C5LREQyyhho2dKGt9dfh++/h6uvhr594dAht6vLuMREW3vnzlC+vO0c3KsX7N0Ld98NSUluV3hJ1E1XREQkHTt2wNO1Yhm/vREBFa8gz8/z7dQvERHJuv75xzaf+9//ICzMdkW/5x4oUMDtys729992nciMGXaENzERgoJss6YmTeCuu6BMGberzJBzddNVGBURETnDuDnbeatVPLP23MHR/LlZ/c007ri91nmfE70ynqExcWxLSCQ0JJCoyHDtvygi4qsWLbJbwaxYAX5+dhuYunXhxhvt1woVvN8b4PhxWLLkVABdvdoeL1fOhs8mTey03Lx5vVuXByiMioiIZMD703fwwcN/Mu1AC/zzHuP+xwexu3hpBrases5wGb0ynj6TV5OYlHzyWGCA/3mfIyIiLktJsSOOixbB4sWwdCns32/vK1r0VDi98UaoXduOSnrS0aOwfbt9/RkzYOZMO+02Vy6oX/9UAL366izfNE9hVERE5AJ++w163jCTrxMf5FBQHlq3epWNxa4AICwkkEW9b0n3efUGzSU+IfGs4+d7joiI+JjkZLu29OefbTj9+WfbPABsQLzuurSjp2XLnh0SHccGyh07bNA8/euZx/buPfW8YsXstNsmTeD226FgQe/93F5wrjCay41iREREfM1PP8EnkV8yPfFxNoeE8cRDr7C9QLGT929LJ2xe6L7zPUdERHyMv7/twFu1KrRvb4/t2WOnzqaG07Fj4e237X2lSkGdOjaAnh40jx07+9x589rHlyxpRzobNjz1fbVqEBFhpwvnMAqjIiLZkNYvXpxvpzosvm84Hx/vyYoy1XmiZW/25w1O85jQkMBzPj80JDDdkdHzPUdERLKAokWhaVN7A7uu87ffToXT2FjbZb1kSQgPt19TQ+bp/y5QIMtPtc0MCqMiItnMmesX4xMS6TPZNkFQID3bJ2NT2N+uB4OckRy95wHi+wwkafoGOGP9Z1Rk+DnPERUZnu6a0fM9R0REsqBcuaBmTXvr3NntarK8nDcWLCKSzQ2NiUsTigASk5IZGhPnUkW+yXHgzTeOEPTUQ3RzRnKs07PkmTyBu+uUZ2DLqoSFBGKw6z4v1IioeY2wi36OiIhITqeRURGRbEbrFy8sJQVe6bqXxu80pwELOD74TXL36n7y/uY1wi46SF7Kc0RERHIyhVERkWxG6xfPLykJej30N09OupOr/TaQ8vkEcj3cyu2yREREchxN0xUR8aLolfHUGzSXcr1nUG/QXKJXxnv8NaIiwwkM8E9zTOsXrcOH4dlbV9Nj0o1UzPM3/j/E4KcgKiIi4gqNjIqIeIm3GgulnkvddNP67z94+eZ5vLGmOQEF85H3p4W2fb+IiIi4QmFURMRLztdYyNNBUesX04qPhxE3TOTNfx7naOmK5Fv0vd2sXERERFyjaboiIl6ixkLumD/P4eMqwxn2TysSq9ahwG8LFURFRER8gMKoiIiXhIYEknLMnyP/FOLI34VI+jcfyYcDKFVAjYUyw7//Qq/7t3C40V30TejB3lvvI2TZLChUyO3SREREBE3TFRHJNMeOwerVsHw5LFsG//x0E39vygWOSfO4bX4OxUZAsWJQtGjar6f/u0IFKF8ejDnHCwpg9w8dPy6JTZ1H8OrhV/HP7c+x10dSqHtn8Pe/8AlERETEKxRGRSRbiV4Z70rjnpQU2LjxVPBctgxWrYKjR+39RYtC7doB3HjLfn49soW9R45QwARzc9kwQvOGsHs37NkDu3dDXBwsXGi/T0lJ+zplysAtt5y6lS6d6T9alrJxI7z18BLaxbbnUVaz79YWBI17S78oERERH6QwKiLZhre61abauRPeeQcWL4bYWNi3zx4PCoJataBzZ6hdG66/Hq68MnVEswBwXYbOn5ICCQmnQupvv8HcuTB9Onz6qX1MpUqngmnDhlC8uMd/zCzh2DEY9do+Cgzsw6iU9zgcEkbKx9EUbHGP26WJiIjIORjHcVwtICIiwomNjXW1BhHJHuoNmkt8Os2AwkICWdT7Fo+9zpEjMGoUvPGG3bfyuutOhc7ateGaayBXJv6pLyXFTv+dO9fe5s+HAwfsfVWrngqnN98MISGZV4evWPiTwzcPTaJXfFdKmp0cfqorwW++Bvnzu12aiIiIAMaYFY7jRJx5XCOjIpJtZHa3WseBSZOgVy/YsgWaNYOhQyE8/CJOcvy4TZIHDkDu3BAQYL+e79+5cqVZKOrnZwPwddfBc8/ZU65YcSqcfvCBDct+fhARAU2aQNOmUKNG9lpvuncvDO64lfoTOzOSGeyrUBO/idMIrlXL7dJEREQkAxRGRSTbCA0JTHdkNDTk8rvVrlhhg99PP9nRxx9+gMaNM/DEpCT75Pnz7W3hwlPDmBcjIMB2Mape3d6uu85+rViRXLn8qFMH6tSBPn3sOtWlS20wjYmBV1+FV16BsDAbSps2tSOnQUEXX4YvcByYOP446zuOpO/BV8gVYDj6+ggKdu+cuUPSIiIi4lGapisi2caZa0YBAgP8Gdiy6iWvGd22DV54wa7RLFYMXn8dnnzyPE1Zjx613Yvmz4cFC+yC0kOH7H3XXAMNGtj5s8WL24WOSUn264X+ffSoLWbVKli7FpJP/Iz58kG1amkDatWqaZLmrl3w3Xd2rWlMDBw8CIGBcOutp8Jp2GUsqfVW06ikJJvlp/RZRtul7anOr+xrcDcFPx9tOzuJiIiITzrXNF2FURHJVjwVjA4fhmHDYPBgOw322WdtKC1Y8IwHJibCkiWnRj6XLLGLSsGGwgYN0gZQTzhyxAbSX3+14TT1tn+/vd/PD666yobTmjXhhhvsfN2gII4etRl52jR727rVPqVGDTvtOLjSLqb89Tvb92fs95cZfwA43YEDMHMmzPnqX3J/F81dh7/mdmZxOCSUwLGj8W/Z/LJfQ0RERDKXwqiISAakpMCECdC7N/zzD9x7LwwZYvf3TOPQIXvH0KE2kBpjRyVTw+dNN0GRIt4r3HHgzz9PBdPUoJqaNnPlsvXVrWtvN96IU7oM69Ybpk2zo6aLFzukpBj88h0hT2gCeUruI7j0AV5/KozWt5RK92Uzo2lUfDx8+y3M/fpfiiyYQovkr7mFuQRwnEMlyhPQ+mFyvxgFBQpc0vlFRETEuxRGRUQu4Oef7QjosmV2QHHECDugmYbj2LTaq5dNTQ88AI8/DvXq+Wbr2j177Gjt4sX2B1y2zA77gp2be+ONJ8NpvamHiPu9CEc2F+PojoIc/y/45GnKlrXb1dSqZQdZa9Wye6eW6z2D9P5fxABbBjXJUImOA7//DlOnwrxv9nDlr9Hcz9fcyhxykUxiaHnyPHI/fq0eyH5dmERERHIAhVERkXNwHOjb127VUqoUDBhg86Wf3xkPjI2Fbt1ssKtZ07asrV/flZovWVKS3bD055/tz7F4sR1RBY76B7C6ZEV+Cb2aTUVKszVfaf5IqsSf+ytyW8mqrFgBGzeeOlXZsnAw/y5SivxH7hL78AtMghSD4xiKBuVl5AM1OX7cLm89fvzU7fTvf/sNfpqyh5p/TeF+7AhoLpI5VqY8AQ8/gHngfgVQERGRLE5hVEQkHSkp0LUrvPMOtG1r82Vw8BkP2r7dLhgdN86u+xw4ENq0SSetZlHbtsHPPzNh5Jdctfl3quz8gzzJx0/efcw/gNwVykGFChwtXYG/cpVn9eEKLN5ZgUmry7I1Pv39PP1IpgD7KcReQkhI91bXLKERc/F3kjl+ZQVytbof7lcAFRERyU4URkVEznD8uA2gn38OPXvaJaBp8s/RozBypG2he/SoncP70kvZdq1iajOiY0ePUWr/bq5I2EHF/Tt5pPhxrjq4CzZtsrcztqY5XLQk6/xLcjTZUDj5ICVNIsFH95Pr0P7zvp5jDE6Fivjdd6+d7ly9ugKoiIhINnSuMKoN2UQkRzp6FFq1guho6N8fXnzxtBzkOHYBY8+eNnzdfbdtrVupkpslZ7rU7rdDY+KI9/PHubIc90eGc9XpXXEdB/791/5eNm+GTZsI2rSJWps325HikCvs2tkM3Ez+/JjsMrosIiIiF01hVERynEOHoHlzmD3bTsvt2vW0O3//3Y6AzpkDlSvbjTlvvz3d83hrf01val4j7Pw/gzG2c1HRolCnjvcKExERkWxHYVREcpSEBLjrLli6FD75xC79BGDfPjs8+u67djPRt96CDh0gICDd85y5v2Z8QiJ9Jq8GyPKBVERERMQbFEZFJMfYtcsOcq5dC199ZfcQBWyDojvusKOiHTtCv34X3CN0aEzcySCaKjEpmaExcR4Po9lxBFZEREREYVREcoS//4bGje3XadMgMvLEHZs3w223wc6d8P3355ySe6ZtCYkXdfxSaQRWREREsit1jhCRbG/jRrsd6I4dMGvWaUH0t9+gXj07d3fOnAwHUYDQkMCLOn6pzjcCKyIiIpKVeSyMGmPyG2NGGmP+NMYkGmMWG2Nqe+r8IiKX4rff4Kab4PBh+PFHG0oBWLQIGjQAf3/46aeLbsYTFRlOYIB/mmOBAf5ERYZ7qHLLWyOwIiIiIt7myZHRj4BIoDVQFZgFzDbGaB6ZiLhiyRKbN3PlsnmzZs0Td3z3nZ2aW6yYDaWVK1/0uZvXCGNgy6qEhQRigLCQQAa2rOrxqbPeGoEVERER8TbjOM7ln8SYQOAAcK/jOFNPO74C+N5xnJfO9dyIiAgnNjb2smsQETndnDlwzz1QsqTdwuXKK0/c8b//QevWULUqzJwJxYu7WeYFnblmFOwIbGYEXxEREZHMYIxZ4ThOxJnHPTUymgvwB46ccTwRqH/2w0VEMs+MGXb7lnLl7IjoySA6ejQ8+qhdJzpvns8HUfDeCKyIiIiIt3lkZBTAGLMYSAZaATuAh4BPgT8cxwk/47HtgfYAZcuWrfXnn396pAYRkTVr7PLPq6+2zYoKFwYcB157DV591Q6Xfvkl5M178jnaOsX79DsXERHJOc41MurJMFoB+Bi4GRtKfwE2ADUdxznngixN0xURT9m3D2rXhv374ZdfIDQUSEmBbt3sqGibNvDhh3YR6QmaBut9+p2LiIjkLJk9TRfHcTY5jtMACAbKOI5zPRAAbPHUa4iInEtKil0KunkzfPXViSCalGSn5Y4eDT16wNixaYIoaOsUN+h3LiIiImDXenqU4ziHgEPGmELY7rq9PP0aIiJnGjQIpk6FkSPh5puxe7ncdx98/z0MHAjPPw/GnPU8bZ3iffqdi4iICHgwjBpjIrEjreuBisBQIA74xFOvISKSnlmz4KWX4OGHoWtXYO9eaNYMfv4ZPvgA2rU753NDQwKJTycEaeuUzKPfuYiIiIBn9xktCIzGhtHPgIXA7Y7jJHnwNURE0ti6FR56CKpUsbnTHD0Cd94Jy5fDxInnDaIAUZHhBAb4pzkWGOBPVGT4OZ4hl0u/cxEREQEPjow6jvMV8JWnziciciGJidCyJSQnw+TJkC/IgSc6wNKlMGmSvfMCUhvmqLOr9+h3LiIiIpAJa0ZFRLzBceCZZ2DlSpg2DSpWBEaOgk8/tVu4ZCCIpmpeI0xByMv0OxcRERFPTtMVEfGa99+HcePg5ZehaVNg9mzbMbdFC+jb1+3yREREROQCFEZFJMtZssQ2KrrzTnjlFWDTJnjgAahcGT77DPz01iYiIiLi6/SJTUSylJ074d57oUwZ+OIL8Dt0AO65x27bMnUqBAe7XaKIiIiIZIDWjIpIlnH8ODz4oN255eefoXBICrR8DNavh5gYKF/e7RJFREREJIMURkUky3j+eZg/Hz7/HK67Dnilnx0NHTUKbr3V7fJERERE5CJomq6IZAkTJ8Lw4dC5Mzz6KHbrltdeg7ZtoUsXt8sTERERkYukMCoiPu/33+HJJ6FePXjzTeC33+Dxx+GGG2DMGLteVERERESyFIVREfFp+/bZLUPz54evv4bc+/fYhkWFCsHkyZAnj9slioiIiMgl0JpREfFZjgNt2sCWLfDjj1CqaBJEPgDbt8NPP0GpUm6XKCIiIiKXSGFURHzWp59CdLSdmlu/PtClu02ln30GtWu7XZ6IiIiIXAaFURHxSdu2QeeuKRS4cj9v7VzEsXvn0nvyaOjRAx57zO3yREREROQyac2oiPgcx4F7WiVyONEh+PaV1IxfR/fot1hYvibRrbq6XZ6IiIiIeIDCqIj4nPHjIfanQEJujqOs/5+8H/0G8QWL8UyzXgyd/Yfb5YmIiIiIB2iaroj4lB07oGtXyBP2H4Wrb+DdCQPIm3SUVq0Gsj9vMAcSEt0uUUREREQ8QGFURHyG40DHjpCYCJWf2ECbn7+k+vaNPN3iBTYVLQNAaEigy1WKiIiIiCcojIqIz5g40XbPHTIEri++i/ojvuarqo2JuaouAIEB/kRFhrtbpIiIiIh4hMKoiPiEXbugc2eoUwe6P30I/4ieHC4Zyoctu2KO2BHRqMhwmtcIc7tUEREREfEAhVER8QmdOsGBA/Dxx+Dfpxf88QdBP/7IDw0auF2aiIiIiGQChVERcd3XX8M338CAAVD57xgYMwa6dwcFUREREZFsyziO42oBERERTmxsrKs1iIh7du+Ga6+FsmVhyXf/katGVQgJgRUrIG9et8sTERERkctkjFnhOE7Emcc1MioiruraFRISYO5cyPVsZ7t4dNo0BVERERGRbE5hVERcM2UKfPklvPYaVFkzESZMgP79oWZNt0sTERERkUymaboi4op//7XTc0uVgmXR2wioUQWuugoWLoRc+juZiIiISHahaboi4lOefdYG0piZDgFPt4UjR+CzzxRERURERHIIfeoTEa+bNg2++AJefhmuW/I+xMTAO+/YkVERERERyREURkXEq/buhaefhqpV4cUHNsL1PeD226FjR7dLExEREREvUhgVEa/q3t02zJ0efZzc7VpD7tzw8cdgjNuliYiIiIgX+bldgIjkDNEr46n8xCrGjYOSN20l72cvwc8/w5gxEBbmdnkiIiIi4mUKoyKS6aJXxtPrf+vYMCmcgCIHqFZ6OpXeHcY/tzWDVq3cLk9EREREXKAwKiKZbmhMHDt+rEDygbyERi5jxMyh/BdUkHZ1ntD0XBEREZEcSmFURDLdlvUBHPjlSvLX/JM+f4whfM9fPH9HV9Yn5Xa7NBERERFxicKoiGSq5GTYP/s6/IKOElnua9otm8IX1e9kXoUIQkMC3S5PRERERFyiMCoimeq99+BQfAEqNPyF4bPf5K+Qkgxo1JbAAH+iIsPdLk9EREREXKKtXUQk0+zYAS+8AI0bw3v5PiBs/y7uf3gIhYoXJioynOY11EVXREREJKdSGBWRTNO9Oxw5AmPbLqLsI59Bly5MGtXT7bJERERExAdomq6IZIrZs2HCBHgp6ihl+7eDMmXgjTfcLktEREREfIRGRkXE444cgWeegYoVoY8zANatg+++g+Bgt0sTERERER+hMCoiHjdkCGzcCIve/51cnQfCI4/AnXe6XZaIiIiI+BBN0xURj/rjDxgwAFrdn0zdj5+CggVh5Ei3yxIRERERH6ORURHxGMeBTp0gTx54r8po+HopjB8PRYu6XZqIiIiI+BiNjIqIx3z1FcyaBW9130rBIS/aqbkPPeR2WSIiIiLigxRGRcQj9u2D556DWjUdHl/cwR587z0wxt3CRERERMQnKYyKiEf07Qs7dsBX94zHzIqBgQOhbFm3yxIRERERH6UwKiKXbcUKeOcd6PXEbsq/9SzceKPd20VERERE5Bw8EkaNMf7GmP7GmC3GmCMnvr5ujFGDJJFsLjkZOnSA4sXhtQPPwv798OGH4O/vdmkiIiIi4sM8FRafBzoBrYHVQDXgU+Ao0N9DryEiPuj99yE2FuZFzSD30P/Bq6/Ctde6XZaIiIiI+DjjOM7ln8SY6cC/juO0Pu3Yp0ARx3Ganu+5ERERTmxs7GXXICLet2MHhIdDg5oHmLrpWkyBAvDLL5A7t9uliYiIiIiPMMascBwn4szjnlozuhBoZIy5+sSLVQZuAb7z0PlFxAd17w5Hj8JnpV/A/POPnZ6rICoiIiIiGeCpabqDgfzAWmNM8onzvuE4zpj0HmyMaQ+0ByirbpsiWdLs2TBhAnzUdjEhn7wDXbrYxkUiIiIiIhngqWm6rYChQBSwBqgOjAKiHMcZe77napquSNYSvTKewdM3EjsigiCTxF8F6lAw5SisWQPBwW6XJyIiIiI+5lzTdD01MjoUGOY4zpcnvl9tjLkC6AOcN4yKSNYRvTKePpNXs/3H8hzfG8xLlZ+i4NqNLH77M+oqiIqIiIjIRfDUmtEgIPmMY8kePL+I+IChMXHs35GXfUsqUKv8XLqt/5QplRsSdTDM7dJEREREJIvx1MjoNKC3MWYLdppuDaA78JmHzi8iPiB+byL/zrqBXP5JfHCoEwfyBPHare1ISEh0uzQRERERyWI8FUa7YPcTHQMUB7YDHwKveej8IuIDAjaX4+hfReh99fPUXL+ers16sjeoIGEhgW6XJiIiIiJZjEfCqOM4B4BnT9xEJBvaswd2zg7nmpIr6btpFHPLR/DtNQ0IDPAnKjLc7fJEREREJIvRmk4RyZCePSHxgB8zi3XD8fPnpchOhBUKYmDLqjSvoTWjIiIiInJxPDVNV0SysR9/hE8/hcl3jaXsdz/Be++x+Ok2bpclIiIiIlmYRkZF5LyOHIGnn4a6V8TTfGEPaNgQ2rVzuywRERERyeIURkXkvAYOhI0bHaJLdsAkJcFHH4Gf3jpERERE5PJomq6InNP69TaMjq47gWKLp8Pw4VChgttliYiIiEg2oOENEUlXSoqdnntl0C46ru8KN9wAXbu6XZaIiIiIZBMKoyKSrnHjYMECmHlVV/wOHoCxY8Hf3+2yRERERCSbUBgVkbPs2mW3culzTTTll0+Evn2hcmW3yxIRERGRbERrRkXkLD16QK4De+nn3xGuuw6ef97tkkREREQkm1EYFZE0Zs+GL76A2Ot6EPD7boiZAQEBbpclIiIiItmMpumKyEmJidCxI7QJnUWtXz+BXr2gZk23yxIRERGRbEhhVEROGjAAtv9xkDHJ7eHqq+Hll90uSURERESyKU3TFREA1q6FwYNhRngfAjf8BZMXQt68bpclIiIiItmURkZF5OSeorfl/Ynb4kZDly5Qt67bZYmIiIhINqaRURHh448hdmEiO0o8CUWuhDfecLskEREREcnmFEZFcridOyEqCsaWfpWC/2y07XSDg90uS0RERESyOYVRkRyuWzcIPxDLQ/uHwVNPwa23ul2SiIiIiOQACqMiOdiECTB54jH+LtEW418Shg1zuyQRERERySEURkVyqH/+gWeegTFhAygRvxq+/RYKFnS7LBERERHJIRRGRXKglBRo0wZqJC7myQOvw6OPQrNmbpclIiIiIjmIwqhIDjR6NCyfs4+/izyCKVDWHhARERER8SKFUZEcZu1aeL6Xw8zQDuTf+TdMX6jpuSIiIiLidQqjIjnIsWPw2GPwVMCnNNj2Jbz+Otxwg9tliYiIiEgO5Od2ASLiPa+9Bvt/2ciI452hYUPo3dvtkkREREQkh9LIqEgOsXgxDBtwjPVFHiJXSm74/HPw93e7LBERERHJoRRGRXKAgwft9NxR+V/iyn9XwOTJULq022WJiIiISA6mMCqSA3TvDhU2/8DTDIWnn4YWLQCIXhnP0Jg4tiUkEhoSSFRkOM1rhLlcrYiIiIjkBAqjItnctGkw5cPdbMr3OFxRGYYPB2wQ7TN5NYlJyQDEJyTSZ/JqAAVSEREREcl0amAkko3t2gVPPekwKf8T5D++FyZMgKAgAIbGxJ0MoqkSk5IZGhPnRqkiIiIiksNoZFQkm3IcaN8eHv5vNDcnz4C33oJq1U7evy0hMd3nneu4iIiIiIgnaWRUJJv65BPYPPU3hpooaNIEOndOc39oSGC6zzvXcRERERERT1IYFcmGtmyB3l0P823QQ/gXLWSTqTFpHhMVGU5gQNqtXQID/ImKDPdmqSIiIiKSQ2markg2k5wMjz8OA4/14MqktRA9C4oVO+txqU2K1E1XRERERNygMCqSzQwbBkUXTuFJ3oOoKLjttnM+tnmNMIVPEREREXGFwqhINrJqFbz30j/8lvspnKq1MK+/7nZJIiIiIiLp0ppRkWziwAF47OFkvvB7jOCAo5gJEyB3brfLEhERERFJl0ZGRbKB48ehVSu4Z/1g6jnz4P1PoFIlt8sSERERETknjYyKZAM9eoD5bjqv0dem0tat3S5JREREROS8NDIqksWNGQOL3oplUa4H8buuBnz00VnbuIiIiIiI+BqFUZEsLCYG3uyyldg8TcldshhMnw758rldloiIiIjIBSmMimRRa9ZA+/v3MjfgLkICj2K+/xFKlnS7LBERERGRDFEYFcmCdu2Clk2OMuFIC8rzByb6B7jmGrfLEhERERHJMIVRkSzmyBFofo/Dq/88Sd3k+fDFF9CggdtliYiIiIhcFIVRkSzEcaBtW7hrSV8eYjy88QY88ojbZYmIiIiIXDRt7SKShfTrB0ETPuIl3oCnnoI+fdwuSURERETkkngkjBpjthpjnHRuMzxxfhGB//0PlvSbyfumA05kpN3TRVu4iIiIiEgW5alpurUB/9O+LwWsAL7y0PlFcrTFi2FE61XM978fc20VzNdfQ0CA22WJiIiIiFwyj4RRx3F2n/69MeZJYD/wtSfOL5KTbdkCzzT7mxinCXlKhOD33QzIn/+sx0WvjGdoTBzbEhIJDQkkKjKc5jXCXKhYREREROTCPN7AyBhjgCeBLxzHOezp84vkJPv2Qas79/G/hLsoGnQQ/5kLIezsgBm9Mp4+k1eTmJQMQHxCIn0mrwZQIBURERERn5QZDYxuA8oBH53rAcaY9saYWGNM7O7du8/1MJEcLSkJWt2bxBsb7uNqsx7/KZOgatV0Hzs0Ju5kEE2VmJTM0Jg4b5QqIiIiInLRMiOMtgOWO46z6lwPcBznA8dxIhzHiShWrFgmlCCStU35JZ4y9f/mgTntaezMZlXfIdC48Tkfvy0h8aKOi4iIiIi4zaNh1BhTHLgH+NCT5xXJSSYtj6dt+2TaL/uEJxjHqLoP8cjxa4heGX/O54SGBF7UcRERERERt3l6ZLQNcBT40sPnFckRDh6Epx7NTdcVn/Aar/DNtbcwov7DF5xyGxUZTmCAf5pjgQH+REWGZ3bJIiIiIiKXxGMNjE40LnoK+NJxnAOeOq9ITrFtG7S46yhvbejOY3zBpCq30OeOLif3Ej3flNvUJkXqpisiIiIiWYUnu+k2BCoBj3rwnCI5wurV8PAd//HujhbUZwFv1n+Et+u2OhlE4cJTbpvXCFP4FBEREZEsw2Nh1HGcHwFzwQeKSBo//ABRLf4g+kgTyvlvJfa1t/joSCU4rTuuptyKiIiISHaTGd10RSSDPv4YXr9zEXOP3MiV+ffgN2c2ES92YWDLqoSFBGKAsJBABrasqlFPEREREclWPDlNV0QyyHGgb1/4440v+cGvDf5XlsF/5ndQqRKgKbciIiIikv1pZFTEy44ehUcfcXDeeIMveYhcN9bGf9mSk0FURERERCQn0MioSAZFr4y/7G61//0H9919jMcWPc0TjMN55BH8xo6FPHkyqWoREREREd+kMCqSAdEr4+kzeTWJJ5oKxSck0mfyaoAMB9LNm6FV5F4Gb7qXRvwIr7yCeeWVNB1zRURERERyCk3TFcmAoTFxJ4NoqsSkZIbGxGXo+UuWwIO1N/PF5ro08F8In34Kr76qICoiIiIiOZZGRkUyYFtC4kUdT5WUBO++C5OiljDz+N2E5D+O39QfoEGDzChTRERERCTL0MioSAaEhgRe1HHHgehouK5yEn90e4tZxxsRUiY//kt/VhAVEREREUFhVCRDoiLDCQzwT3MsMMCfqMjwsx67fDk0bAjvtZjJtL+u4y26kbvxzfgvXwLhZz9eRERERCQn0jRdkQxIbVJ0vm66f/4JL7wAsf+LY3Tu7tzGdzhlK8KbUzHNmml9qIiIiIjIaRRGRTKoeY2wdDvn7tsHAwfCuBF7eTH5NT7zG41f3iAYMAzTubO2bRERERERSYem6YpcoqQkeOcdCK9wnP2Dx7DRVKJzyij8n2qL2bgRevRQEBUREREROQeNjIrPiF4Zf95psL7CcWDaNOjVC0rHzWZxvucoz+9wY0MYORKuu87tEkVEREREfJ7CqPiE6JXx9Jm8+uRenvEJifSZvBrAZwJpUhIsWACvvw7/zNvIe8E9uZVvcUqUh2GToXlzrQsVEREREckgTdMVnzA0Ju5kEE2VmJTM0Jg4lyqy9u2DiRPhkUegeHG4r/Fe7l0axXr/a7mFuTBoEGbNGmjRQkFUREREROQiaGRUfMK2hMSLOp6Ztm6103C//RbmzYPcxw/xcP5pzC40keqHvsfvyDHME0/AG29AyZJer09EREREJDtQGBWfEBoSSHw6wTM0JDBTXu/09amlCgTSomwV9scV59tv4bffIA9HaBf2PcMqfEmVP6fjf+Aw5A+FTh3hiSegWrVMqUtEREREJKdQGBWfEBUZnmbNKEBggD9RkeEef61vlsXTa9xmDm4vSOKWCvy9qQQ/H8xLHnOUbpVn879aE7l6fTT+8QegWDF4ojW0agX164OfZraLiIiIiHiCwqj4hNQmRZ7sprtvH6xff/YtbkMpnBR73lwBR7izRDQPFppIs91zKLjmAISEQKsH4MEHoVEjyKX/TEREREREPE2fssVnNK8RdtHhMzkZ/voLNm48O3Ru337qcblyQaVKUCU8iaK5v6eWE0vtxF+IjF9I0X/2cSB3ID9UuoGWg3vAbbdB7twe/ulEREREROR0CqPi844ft4Hzjz9s6Dz96+bNdsuVVAULwjXXQGQkXH1VCrULbuCag8sp/ncs/iuWQ8xKOHIEgP158rHgyhpMq3wz88rVomixEFo2ucWjtWeVvVNFRERERLxNYVR8RkoKrFoFP/+cNnRu2ZI2cAYFQcWKcO21dmvPihWhYgWHa4P/pOiW5ZjY5bB8OUxaAQcOnHpSzZrQsSOxRcvzUnwgccHFcYxdA5oZ61Ozwt6pIiIiIiJuURgVV23dCrNnww8/wJw58O+/9nhq4Kxa1W7hWamS/b5SJShVIgWzeROsXAm//AJfn/i6Z499ckAAXHcdPPoo1K4NERF2uPTE2s8IoIMXRizPt3eqwqiIiIiI5HQKo+JVe/fCjz/a8Dl7th39BAgNhSZN7HLNBg2gdGkwBjskun69DZtTToTOVatOjXgGBNgh0mbNbPCsXdsm2Dx5zlvHpaxPvVi+tHeqiIiIiIivURiVTHX0qJ12mxo+Y2PtdNzgYGjYELp0gcaN7cClMdh0+t0cGzpXrrSbfh49ak8WFGRHPB9/HGrUsNNuK1e+YPB0i7f3ThURERERyUoURiVTbNgAr70GU6bA4cPg7w916sBLL9nRzzp17KAmKSk2ob40FaKjYe1ae4KQEBs2O3e2X2vUgKuusifKIry5d6qIiIiISFajMCoetXUr9O8Pn35qByxbt4Y77rBTbwsWPPGgo0dh9lyYOhW+/dbuweLvDzffDO3b2/m6FSqcGCrNujJj71QRERERkexCYVQ8Yts2eOMN+PBD8POz029794YSJU48YO9eGP+dDaDffw8HD0K+fDap3nOPDaCFC7v6M2QGb6xNFRERERHJihRG5bLs3g2DB8M779j9QJ98El58EcqUAXbsgLe/tgF0/nz7gBIl4KGHbAC99VbIm9ftH0FERERERFygMCqXJCEB3nwTRo60a0IffRReftnOrmXfPnhpKAwfDomJcPXV0KOHDaB16tihUxERERERydEURuWiHDwIb70FQ4faQHr//dCvn+2Gy9GjMOo9u2j033+hVSvo29d2vBURERERETmNhqgkQxITYcQIKF/eTsOtX9/uvPLVV3BNeAqMH29HQJ99FqpXtx1yJ0xQEBURERERkXQpjMoFjfp6J4XCEuneHY4V+I9B43YxbZrNnMyaBbVq2Xm6ISEQEwOzZxPtV5J6g+ZSrvcM6g2aS/TKeJd/ChERERER8SWapivn9co7e3i9R2FMQDLFH1xC4JX/Mu4Pf64f79Bo3HCYPRuuvNKOjLZqBX5+RK+MT7O/ZnxCIn0mrwZQZ1kREREREQEURuU83nsPXutamIAiByl+33JyFThCmYQdRC34jEbrFkCRIraDUYcOdlPRE4bGxJ0MoqkSk5IZGhOnMCoiIiIiIoDCqKQjORl69bLNcAPL76Ho3b9Q7Pi/dJ49kUdWfk+ynz+jb3yQzt+/DwULnvX8bQmJ6Z73XMcvV/TKeIbGxLEtIZHQkECiIsMVekVEREREfJzCqKRx6BA88ojdGrRLF4gtsYYya1bx3pQB5D96iInX3c6oug8RUKY0ndMJogChIYHEpxM8Q0MCPV6vpgSLiIiIiGRNamAkJ23bBg0awLRpdvuWt96Ct1N+4/OJffk3qCB3tB3Ni5GdOVC4OFGR4ec8T1RkOIEB/mmOBQb4n/c5l+p8U4JFRERERMR3aWRUAPj1V2jaFPbuhW+/hSZ3pkDfV6j1+uvsur4+ne6IYtNRf8IyMA029T5vTJ319pRgERERERHxDIVR4bvv4MEH7fLPhQuhengiPPwETJwITz5J8XffJSYg4KLO2bxGmFemyXpzSrCIiIiIiHiOpunmcKNHQ7NmUKkSLF0K1UN3wS232CA6ZAh8+CFcZBD1Jm9OCRYREREREc/RyGgOlZwM3bvbdaF33223CQ3+ay00aQI7d8KkSdCypdtlXpA3pwSLiIiIiIjnKIzmQAcPwkMPwfTpNpAOGQL+c3+A++6DoCCYPx9q13a7zAzz1pRgERERERHxHE3TzWH27IGbboLvv4cxY+DNN8H/4w/hzjvhiivsXN0sFERFRERERCRr0shoDpKUBA88AOvW2VHRO25PgajnYdgwG0a//BIKFHC7TBERERERyQE8NjJqjClljPnUGLPbGHPEGLPWGNPAU+eXy9ejB/z4I3zwAdxx0yG4914bRDt1svu5KIiKiIiIiIiXeGRk1BgTAiwCFgJNgN1AeWCXJ84vl2/sWHj7bbtG9PHG26DB3bByJYwaBV27ul2eiIiIiIjkMJ6aptsL2O44zuOnHdvioXPLZVq8GDp2hNtug8EdtkCdm2HvXpg6FZo2dbs8ERERERHJgTw1Tbc5sNQYM9EYs8sYs8oY09kYYzx0frlEf/9td2i54gqY+OF+cjVvCocOwcKFCqIiIiIiIuIaT4XR8sAzwGYgEhgFDAI6pfdgY0x7Y0ysMSZ29+7dHipBzpSYCC1a2Ow5dXIyhZ55CDZsgG++gerV3S5PRERERERyME9N0/UDYh3H6XPi+5XGmErYMDr6zAc7jvMB8AFARESE46Ea5DSOA+3awS+/QHQ0VB7XC777Dt59F265xe3yREREREQkh/NUGN0OrD3j2Dqgm4fOLxdp2DAYPx5efx3u3j0Whg+HLl2IrtOMoYPmsi0hkdCQQKIiw2leI8ztckVEREREJIfxVBhdBISfcewq4E8PnV8uwsyZ8PzzcP/98EL9BXBbR7j9dqY+1oM+k1eTmJQMQHxCIn0mrwZQIBUREREREa/y1JrREcANxpgXjTEVjTH3A12Bdzx0fsmgDRugVSuoVg3GvbwZc29LKF8eJk5kyJxNJ4NoqsSkZIbGxLlUrYiIiIiI5FQeGRl1HGe5MaY5MADoC/x14usYT5w/u4teGc/QmLjLnjq7bx/cfTcEBMDUz/YR9GAzSEmBadMgJIRtCYnpPu9cx0VERERERDKLp6bp4jjODGCGp86XU0SvjPfI1NnkZHjkEdi0CWbHJHNF7xOdc2fNgkqVAAgNCSQ+neAZGhLogZ9EREREREQk4zw1TVcu0dCYOI9Mne3bF2bMgFGjoMH0KPj+exg9Gho1OvmYqMhwAgP80zwvMMCfqMgzl/uKiIiIiIhkLo+NjMql8cTU2S+/hIEDoX176BjwEYwYAV27wtNPp3lc6kirJ6YEi4iIiIiIXA6FUZdd7tTZX36Btm2hXj0Yfd88zF0dITIS3nwz3cc3rxGm8CkiIiIiIq7TNF2XXc7U2V27oHlzKFoUot/cRECre6FiRTtUmkt/ZxAREREREd+lxOKyy5k627Ej7N4NP8/cR9EnmtmDJzrnioiIiIiI+DKFUR9wKVNnp06FyZNh8BvHqT6oFWzcaDvnVqyYSVWKiIiIiIh4jsJoFrR/P3TqBNWqQY+dUTBzJrz/fprOuSIiIiIiIr5Ma0azoJdegm3bYHKLz/F/ayR062Zb6YqIiIiIiGQRCqNZzNKldvvQF9pso8LILlC/Pgwb5nZZIiIiIiIiF0VhNAtJSoJ27SC0lMMrOzvC0aMwdqw654qIiIiISJajMJqFDB8Oq1dD9EMTCfjuW+jfH666yu2yRERERERELprCaBaxaRO8+iq0abKbiE+7wPXXw3PPuV2WiIiIiIjIJVEYzQIcBzp0gIAAeMe/C+zbBx9/DP7+bpcmIiIiIiJySRRGs4Dx42H2bJj4UDRB306Evn3h2mvdLktEREREROSSGcdxXC0gIiLCiY2NdbUGX7ZnD1xzDdQst5eZf1fGlCgBy5fbYVIREREREREfZ4xZ4ThOxJnH1YbVx/XsCQkJ8GVod8wvu2HGDAgIIHplPENj4tiWkEhoSCBRkeE0rxHmdrkiIiIiIiIZojDqw+bMgU8/hY9bxVDoy3HQpw/UrEn0ynj6TF5NYlIyAPEJifSZvBpAgVRERERERLIETdP1UYmJUK0aBKfsZ8WxKvgF54OVKyFvXuoNmkt8QuJZzwkLCWRR71tcqFZERERERCR9mqabxbz+OvzxB/zdrDd+0/+BRYsgb14AtqUTRM93XERERERExNeom64PWr0ahgyBgZHzKD3tXXj2WbjxxpP3h4YEpvu8cx0XERERERHxNQqjPiYlBdq3h1IFDxMV9xSUL2+HSU8TFRlOYEDaPUYDA/yJigz3ZqkiIiIiIiKXTNN0fcx778GSJbDmjr74z9wEc+dCUFCax6Q2KVI3XRERERERyaoURn1IfDz07g1dai/hmpgR0KEDNGqU7mOb1whT+BQRERERkSxL03R9SJcuYI4dZdh/bTGlS8PgwW6XJCIiIiIikik0Muojpk6FKVNgUaP+5P5xHXz/PRQo4HZZIiIiIiIimUIjoz7g6FHbMPf+iiu5ccEgaN0a7rjD7bJEREREREQyjUZGfcD778M/W5P4tUJbTLFiMHy42yWJiIiIiIhkKoVRl+3fD/37w3vlh1Jg0yo7V7dwYbfLEhERERERyVQKoy57803IvSeeNgdeh5YtoXlzt0sSERERERHJdFoz6qKdO20Y/eKKl/B3kmHoULdLEhERERER8QqFURe9/jqEJ66i4V+fQteuUL682yWJiIiIiIh4habpumTTJnjvXYfVJXtgjhaGF190uyQRERERERGvURh1Sd++cI//dK7eNhfefhtCQtwuSURERERExGsURl2wciV8PSGJbYWjoFg4PP202yWJiIiIiIh4lcKoC3r3hueCPqDYf3HwyVQICHC7JBEREREREa9SGPWyuXNh6awEooNegYYNoVkzt0sSERERERHxOnXT9SLHsaOig/IPIG/if3ZfF2PcLktERERERMTrNDLqRZMmwe7lW2iXaxTm8cehZk23SxIREREREXGFRka9JCkJXngB3i3QG78Af3jjDbdLEhERERERcY1GRj0semU8Q2Pi2JaQSGhIIFGR4TSvEcbHH0ORjT9zB1/ZfV3CwtwuVURERERExDUKox4UvTKePpNXk5iUDEB8QiJ9Jq/mSKKh36ulmJW/O06+kphevVyuVERERERExF0Kox40NCbuZBBNlZiUzPP9Eqm/42uqsARGfATBwS5VKCIiIiIi4hsURj1oW0LiWceSEwPY82Nx3g56HipWgzZtvF+YiIiIiIiIj1EDIw8KDQk869j+nyvSOeldShzeardy8ff3fmEiIiIiIiI+RmHUg6IiwwkMOBU2j+/PS54V+Xgl1xtw113QuLGL1YmIiIiIiPgOj4RRY8yrxhjnjNsOT5w7K2leI4yBLasSFhKIAY4tq8yrvEY+5xAMHep2eSIiIiIiIj7Dk2tG44CGp32ffI7HZWvNa4TRvEYYv/8ODwxZz9O8j2nfDipXdrs0ERERERERn+HJMHrccZwcNxp6Li+8AG/698IvbxD06+d2OSIiIiIiIj7Fk2tGyxtj4o0xW4wxXxpjynvw3FnKwoVwcNpc7jw+DfPiC1C8uNsliYiIiIiI+BTjOM7ln8SYO4H8wHqgOPAScDVwreM4/6bz+PZAe4CyZcvW+vPPPy+7Bl/hONCgfjJjlkVQudR/+G2Ig7x53S5LRERERETEFcaYFY7jRJx53CPTdB3H+f6MF1sCbAZaA8PTefwHwAcAERERl5+Gfcjs2VB+8edUYRUMGq8gKiIiIiIikg5Prhk9yXGcg8aYNUClzDi/r3IcGPhyIv/ze5GUWtfj16qV2yWJiIiIiIj4pEzZZ9QYkxc7TXd7ZpzfV82bB9WXvEvJlG34DRkMftrGVUREREREJD0eGRk1xgwDpgF/YdeM9gXyAZ964vxZxZCXD/KF30CSGzbGv2FDt8sRERERERHxWZ6aplsamAAUBXYDS4AbHMfJPp2JLmDBAqix8C2KsAfe6O92OSIiIiIiIj7NUw2McvziyOEvJzDODCX5jqb433CD2+WIiIiIiIj4tExpYJTTLF4MNeaPIIQEeOM1t8sRERERERHxeQqjHjDypT2MNSM4fs995KpRw+1yREREREREfJ7C6GVauhRq/TiUYHMQ80Y/t8sRERERERHJEhRGL9PbL+7gA97m+IOPEFC5stvliIiIiIiIZAkKo5dhxQqoPWcgef2O4df/FbfLERERERERyTL83C4gK3v3hb/pwHscf6QNVKzodjkiIiIiIiJZhkZGL9GqVVB71uvk8nfw79/X7XJERERERESyFI2MXqIPnt9EWz7m+BPt4Yor3C5HREREREQkS1EYvQSrV0OdWa9Brlzk6feC2+WIiIiIiIhkOQqjl+DjXut5lC84/nQnCA11uxwREREREZEsR2H0Iq1dC3VmvsrxgEACX3ne7XJERERERESyJIXRi/RZz99oxUSSOz8LxYq5XY6IiIiIiEiWpDB6EeLi4MbvXyYxT0GC+vZwuxwREREREZEsS2H0IvzvueXcw1SSn+0JhQq5XY6IiIiIiEiWpTCaQX/8AXW/78uhvEUIfrGb2+WIiIiIiIhkaQqjGfR115+IJIbkqN6QP7/b5YiIiIiIiGRpCqMZsGWzQ73vX2JfUEkK9H7G7XJERERERESyPIXRDIjuMoebWYDT50UICnK7HBERERERkSxPYfQC/tzqUPf7l/gvuAwhUe3cLkdERERERCRbUBi9gBnPzKCOsxT6vgx58rhdjoiIiIiISLagMHoe//yVwo0zX2ZXgQoUfq612+WIiIiIiIhkG7ncLsCXJX33AzWclex+9XMICHC7HBERERERkWxDYfQ8yj19O5SbSbHGjd0uRUREREREJFtRGD0fYyAy0u0qREREREREsh2tGRURERERERGvUxgVERERERERr1MYFREREREREa9TGBURERERERGvUxgVERERERERr1MYFREREREREa9TGBURERERERGvUxgVERERERERr1MYFREREREREa/L5XYBvix6ZTxDY+LYlpBIaEggUZHhNK8R5nZZIiIiIiIiWZ7C6DlEr4ynz+TVJCYlAxCfkEifyasBFEhFREREREQuk6bpnsPQmLiTQTRVYlIyQ2PiXKpIREREREQk+1AYPYdtCYkXdVxEREREREQyTmH0HEJDAi/quIiIiIiIiGScwug5REWGExjgn+ZYYIA/UZHhLlUkIiIiIiKSfaiB0TmkNilSN10RERERERHPUxg9j+Y1whQ+RUREREREMoGm6YqIiIiIiIjXKYyKiIiIiIiI1ymMioiIiIiIiNcpjIqIiIiIiIjXZUoYNca8YIxxjDGjM+P8IiIiIiIikrV5PIwaY24A2gG/efrcIiIiIiIikj14NIwaYwoC44Engb2ePLeIiIiIiIhkH54eGf0A+MZxnLkePq+IiIiIiIhkI7k8dSJjTDugIvBYBh7bHmgPULZsWU+VICIiIiIiIlmER0ZGjTHhwADgEcdxjl3o8Y7jfOA4ToTjOBHFihXzRAkiIiIiIiKShXhqZPRGoCjwuzEm9Zg/cLMxpgOQz3Gcox56LREREREREcniPBVGo4HYM459AmzEjphecLRUREREREREcg6PhFHHcRKAhNOPGWMOAf85jvO7J15DREREREREsg+P7zMqIiIiIiIiciEe66Z7JsdxGmbWuUVERERERCRrM47juFuAMbuBP10t4sKKAnvcLkJyPF2H4gt0HYqv0LUovkDXofiCrHAdXuE4zlnbqLgeRrMCY0ys4zgRbtchOZuuQ/EFug7FV+haFF+g61B8QVa+DrVmVERERERERLxOYVRERERERES8TmE0Yz5wuwARdB2Kb9B1KL5C16L4Al2H4guy7HWoNaMiIiIiIiLidRoZFREREREREa9TGBURERERERGvUxg9D2PMM8aYLcaYI8aYFcaYm9yuSXIOY8yrxhjnjNsOt+uS7M8Yc7Mx5ltjTPyJ667NGfebE9fnNmNMojFmnjHmWpfKlWwqA9fhuHTeI5e4VK5kU8aYPsaY5caY/caY3caYacaYKmc8Ru+JkqkyeB1myfdEhdFzMMY8CIwCBgA1gMXA98aYsq4WJjlNHFDqtFtVd8uRHCIY+B3oBiSmc38voAfQBagN7AJ+MMbk91qFkhNc6DoEmE3a98i7vFOa5CANgTFAXeAW4Dgw2xhT+LTH6D1RMltDLnwdQhZ8T1QDo3MwxiwFfnMcp91pxzYC3ziO08e9yiSnMMa8CtznOE6VCz1WJLMYYw4CnR3HGXfiewNsA0Y7jvPGiWOB2A9fPR3Hed+tWiX7OvM6PHFsHFDUcZymbtUlOY8xJhjYBzR3HGea3hPFDWdehyeOjSMLvidqZDQdxpjcQC1g1hl3zcL+RULEW8qfmKK2xRjzpTGmvNsFSY5XDijJae+PjuMkAgvQ+6N4X31jzC5jzAZjzIfGmOJuFyTZXn7s5+e9J77Xe6K44czrMFWWe09UGE1fUcAf2HnG8Z3YNxwRb1gKtAHuBNphr73FxpgibhYlOV7qe6DeH8VtM4HHgVuxUySvB+YaY/K4WpVkd6OAVcDPJ77Xe6K44czrELLoe2IutwvwcWfOYTbpHBPJFI7jfH/69ycWoW8GWgPDXSlK5BS9P4qrHMf58rRvVxtjVgB/Ak2Aye5UJdmZMWY4UB+o7zhO8hl36z1RvOJc12FWfU/UyGj69gDJnP0XreKc/ZcvEa9wHOcgsAao5HYtkqOldnTW+6P4FMdxtgH/oPdIyQTGmBHAQ8AtjuNsPu0uvSeK15znOjxLVnlPVBhNh+M4x4AVwG1n3HUbtquuiNcZY/ICVwPb3a5FcrQt2A9fJ98fT1ybN6H3R3GRMaYoEIbeI8XDjDGjgIexAWD9GXfrPVG84gLXYXqPzxLviZqme27Dgc+NMcuARUAHIBR4z9WqJMcwxgwDpgF/Yf/C2hfIB3zqZl2S/Z3o0lfxxLd+QFljTHXgP8dx/jLGjAReNMasBzYALwEHgf+5UK5kU+e7Dk/cXgUmYT9oXQkMxHYwneLlUiUbM8a8AzwGNAf2GmNSR0APOo5z0HEcR++JktkudB2eeL98lSz4nqitXc7DGPMMdu+oUti9zp5zHGeBu1VJTmGM+RK4GdtQazewBOjrOM5aVwuTbM8Y0xD4MZ27PnUcp82JrQxeAZ4GCmGbbXVyHOd3rxUp2d75rkOgIxCN3Qc8BPvh60fse+TfXilQcgRjzLk+KPdzHOfVE4/Re6Jkqgtdhye2E4omC74nKoyKiIiIiIiI12nNqIiIiIiIiHidwqiIiIiIiIh4ncKoiIiIiIiIeJ3CqIiIiIiIiHidwqiIiIiIiIh4ncKoiIiIiIiIeJ3CqIiIiIiIiHidwqiIiIiIiIh4ncKoiIiIiIiIeN3/AQ/7Rz8x5efDAAAAAElFTkSuQmCC\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": {},
   "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": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept           4.964910\n",
       "x1                  0.495433\n",
       "np.sin(x1)          0.461495\n",
       "I((x1 - 5) ** 2)   -0.018963\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": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    11.025497\n",
       "1    10.900628\n",
       "2    10.671745\n",
       "3    10.380090\n",
       "4    10.079954\n",
       "5     9.825382\n",
       "6     9.656943\n",
       "7     9.591795\n",
       "8     9.619488\n",
       "9     9.704517\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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
