{
 "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.988\n",
      "Model:                            OLS   Adj. R-squared:                  0.987\n",
      "Method:                 Least Squares   F-statistic:                     1248.\n",
      "Date:                Mon, 07 Dec 2020   Prob (F-statistic):           4.71e-44\n",
      "Time:                        17:22:22   Log-Likelihood:                 8.8151\n",
      "No. Observations:                  50   AIC:                            -9.630\n",
      "Df Residuals:                      46   BIC:                            -1.982\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.0615      0.072     70.213      0.000       4.916       5.207\n",
      "x1             0.4840      0.011     43.534      0.000       0.462       0.506\n",
      "x2             0.4200      0.044      9.610      0.000       0.332       0.508\n",
      "x3            -0.0184      0.001    -18.830      0.000      -0.020      -0.016\n",
      "==============================================================================\n",
      "Omnibus:                       13.108   Durbin-Watson:                   1.819\n",
      "Prob(Omnibus):                  0.001   Jarque-Bera (JB):               17.300\n",
      "Skew:                           0.880   Prob(JB):                     0.000175\n",
      "Kurtosis:                       5.282   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.60194513  5.03816284  5.44086647  5.7871662   6.06243306  6.26270236\n",
      "  6.39532512  6.4777603   6.53470635  6.59404332  6.68225225  6.82006441\n",
      "  7.0190552   7.27974236  7.59150091  7.93430879  8.28203642  8.60673974\n",
      "  8.88325133  9.09331534  9.22858689  9.29200308  9.29730047  9.26675817\n",
      "  9.22753747  9.20721931  9.22927278  9.30919928  9.45198562  9.65128472\n",
      "  9.89045845 10.14531117 10.38806506 10.5919239  10.73547564 10.80621031\n",
      " 10.80257539 10.73423089 10.62046224 10.48701172 10.36184911 10.2705763\n",
      " 10.23222101 10.25611057 10.34033946 10.4720824  10.62970147 10.78630122\n",
      " 10.91414745 10.98922392]\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": [
      "[10.98615378 10.87138808 10.66139777 10.39371801 10.11775831  9.8827053\n",
      "  9.72548019  9.66169913  9.68184982  9.7536205 ]\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": "iVBORw0KGgoAAAANSUhEUgAAA6MAAAHWCAYAAACCFl9HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAB0dklEQVR4nO3de3zO5R/H8de1GcYwZzbn00qImg4OOSSHJEPnFBVKFMqxE1GhlSj0SwdCiUqTVHNOKJpTSKuEsjmznIbZvr8/LsOY8737u3t7Px+P+zH73vf93efm293e93Vdn8s4joOIiIiIiIiIN/m5XYCIiIiIiIhkPwqjIiIiIiIi4nUKoyIiIiIiIuJ1CqMiIiIiIiLidQqjIiIiIiIi4nUKoyIiIiIiIuJ1OdwuoEiRIk65cuXcLkNEREREREQywIoVK3Y7jlP0zOOuh9Fy5coRExPjdhkiIiIiIiKSAYwxW9I7rmm6IiIiIiIi4nUKoyIiIiIiIuJ1CqMiIiIiIiLidQqjIiIiIiIi4nUKoyIiIiIiIuJ1rnfTvZD9+/ezc+dOkpKS3C5FspCAgACKFStG/vz53S5FRERERCRbytRhdP/+/ezYsYPQ0FACAwMxxrhdkmQBjuOQmJhIXFwcgAKpiIiIiIgLMvU03Z07dxIaGkqePHkURMVjjDHkyZOH0NBQdu7c6XY5IiIiIiLZUqYOo0lJSQQGBrpdhmRRgYGBmv4tIiIiIuKSTB1GAY2ISobRtSUiIiIi4p5MH0ZFREREREQk61EYFREREREREa9TGM0AHTt2xBiDMebkFiKNGjVizJgxl7RGceHChRhj2L17dwZWKyIiIiIi4n0XFUaNMbcYY742xsQZYxxjTMcz7m9rjIk2xuw6cX/DDKjVpzRp0oRt27axefNmZs+eTatWrRg4cCD169fn0KFDbpcnIiIiIiLiqosdGQ0C1gE9gMR07s8LLAWe8VBdHhW1Ko66w+ZTvv8s6g6bT9SquAz/mbly5aJEiRKEhoZSs2ZNnnnmGRYuXMjKlSt5/fXXAZg8eTK1a9cmX758FCtWjLvvvvvk3pebN2+mUaNGABQtWhRjDB07dgTg+++/p379+hQsWJBChQrRrFkzNmzYkOGvSURERERExFMuKow6jvOt4zjPOY7zBZCSzv2THMd5GfjO0wVeqahVcQyYvpa4hEQcIC4hkQHT13olkJ6pWrVqNG/enC+//BKAY8eO8fLLL7NmzRq++eYbdu/ezf333w9A6dKlTz5u/fr1bNu2jVGjRgFw6NAhevbsyfLly1m4cCEFChSgVatWHDt2zOuvSURE5HRufAAsIiK+KYfbBWS0yOhYEpOS0xxLTEomMjqWiFqhXq+natWqzJ07F4BHH3305PEKFSrw7rvvcvXVV7N161ZKlSpFoUKFAChWrBhFihQ5+dh27dqlOef48ePJnz8/y5cvp169el54FSIiImdL/QA49f+7qR8AA678P1cko0StiiMyOpb4hERCggPp0yxM17jIZXAljBpjugBdAMqUKZOhPys+Ib1Zxec+ntEcxzm5v+XKlSt5+eWXWb16NXv37sVxHAD++ecfSpUqdc5zbNy4kRdffJFly5axa9cuUlJSSElJ4Z9//vHKaxAREUlPZvsAOLtQMPKuqJh/iBr5KV02LKbw4f0kGz/MJ/78U64QZYrmA3//tDc/v7OPpR6vVAnatYPcud1+WSKucCWMOo4zDhgHEB4e7mTkzwoJDiQuneAZEhyYkT/2nH777TcqVKjAoUOHaNasGU2aNGHSpEkUK1aM3bt3U79+/QtOt23VqhWhoaG89957hIaGkiNHDqpWrappuiIi4qrM9gFwdqDRaC9JTobFi2HaNG6Z+CkRBxM4HJCL+HxF8XNS8HdSCNjqQFCAfezpt5SU9I+l6tEDOnWCrl2hbFn3XqOIC7L8NN0+zcLSvEkDBAb406dZmNdrWbduHd9//z0vvPACv//+O7t37+a1116jfPnyAEyfPj3N43PmzAlAcvKp2vfs2cOGDRsYM2bMyQZHK1eu5Pjx4156FSIiIunLbB8AZwcajc5AKSmwZAlMmwZffAHbt0OePCwtfT3TC97O1/vuwilzhHzXbQHAAJuGtby4czuOPf8PP8Do0RAZaW+tWkG3btCkCZyYSSeSlWX5MJr6Ruzt6StHjx5l+/btpKSksGvXLubNm8drr73G9ddfT+/evTl8+DC5cuVi9OjRdOvWjQ0bNvDiiy+mOUfZsmUxxjBr1ixatWpFYGAgBQsWpEiRIrz//vuULl2auLg4+vTpQ44cWf6fUkREMrnM9AFwdqHRaA9LSYGff7YB9PPPIT4ecufGadmSzbXv4b2tLRkxIQdJsbnwy32M/CU3nnzqJX3oYoydqtu4sb398w+89x6MGwczZkBYGDz5JHToAAUKZMALFckcLirBGGOCgEonvvUDyhhjagJ7Hcf5xxhTCCgDBJ94TCVjTAKw3XGc7R6t+DJE1Ar1+qeDc+fOpWTJkvj7+xMcHEy1atUYOHAgjz/+ODlz5iRv3rx8/PHHPPfcc4wZM4YaNWowYsQImjdvfvIcoaGhvPzyyzz//PN06tSJhx9+mAkTJjB16lSefvppqlWrRqVKlXjzzTfPamokIiLiban/r339u1j+jU+hVAl/+rWsohG6DKTRaA/55ReYMsUG0K1bIVcuaNGCPbfew8d7W/HRtCDWfwkBAXBdvUS2Fl6Jf9kdmBx2uu0Vf+hSpgy8+iq8+KKtYcwYO333uefg4YftaOk113joxYpkHia1ac55H2RMQ2BBOnd97DhOR2NMR2B8Ove/7DjOoPOdOzw83ImJiUn3vg0bNnD11VdfsD6Ry6VrTETEMw4dgjlzYOZMmDULduywxwsVghIlzr4VL572+yJFbD8XuTRnrhkFG4yGtq2uDwEuxp9/wrPP2gs3Z05o3pzEVvcwPakVH32RnwUL7IzaOnXgoYfgnnvsNe2VplG//GJD6WefwdGj0LChDaWtW9tULOJDjDErHMcJP+v4xYTRjKQwKm7SNSYicvn+/Re++cb+Hj9/vv19uUABaN7c/vL+3392md3pt23bIDGdGaT+/lCyJNx9Nzz1FJxopyAXQd10L0NCArzyCrz9NuTOTfKA55lf6XHGfxVMVJS9RitWtAG0fXv7Z9fs3g0ffghjx9rpvKGh8MQT0Lmz/VRHxAcojIqkQ9eYiMjFS0mBFSts+Jw5E1avtscrVrR9V1q1gvr1zz9o4zhw8KANpjt2pA2qGzbA11/bn9O6NfTsac+nPi7iMcnJ8MEHdjrs7t04jzzKjNqv0G1ICeLj7ajnvffaEHrTTZns2ktOttMORo+20xACAuCBB+Ctt6BgQberEzmvc4VRdb0RERGRczpyBGbPPjX9dts2O522Th0YPtwG0Kuuuvhf2o2BfPnsrXLlU8ejVsWxIDqWEsUdnA0VmTO/DF995UetWnbp3H332WV8IudywRHiBQvsJxy//gr167O9/0g6jb2OWV2hdm07I/b22+1s3UzJ3x/uvNPeYmNtwf/7n+3I+8UXcP31blcocsm0OkNERETStXix7ZnSujVMnQr16sHEiXZE88cfoW9fuPrqKx89Sl33GJeQiH/+I+S4cT0hT8yl6wv7OHoUOna02y++/PKptagipzv9GnI4td9q1Ko4+PtvaNvWdq3dv5+UqZ8z5p4fqHzvdSxYYAcWf/oJIiIycRA9U1iYnWL84492xLROHduJ1+UZjyKXSmFURERE0jh6FPr1g1tusb/bfv21XbY2bZqdvlikiGd/Xnp7ZR4liTV5V7FunR2ZDQ+HQYNs09GOHU9NERaB9K8h/4MHONCrt/3EZPZsePVVfv9qA/VH3UX3pwx16sD69Xaw1N/fnbqv2I03wsqV0KgRPP643Qrm8GG3qxK5aJqmKyIiIietXm13kli71vZHefNNO6X2oiQn2+HUr76CZcvsfN6AAMiRI+3tjGPdV23juJ8/R/0DiClVlUXlr+NwzkDiExIxBm67zd7++APeeQfGj4ePP4YGDWyQaNXKh8OEeMTp+6r6pSRz19p59PlxIkUPJUCHDhwb9BrDJobw6o0QFGRH+Nu3t6P6Pt8AqkgRO4f+lVfs9IFVq+DLL6FKFbcrE7kghVERERHh+HGIjISBA6FwYdslt2XLi3jikSMwb54NoDNm2CHUXLnsiE1AgD3xkSP2a1KS/XrGrUnCIczx4+Q9doROMTM46h/AkrLXElOjHsTXgpAQwP5u/c47MGQIfPSRnaXYpo1d7/fxx3YATLKn1P1Wq2/7k6HRo6m2YyMxoVfTv+NrdG7flU532FHQ+++HkSOhWDH7vDO3xkmd3gv4ViD197f/8d50Ezz4oJ1K8NFHcNddblcmcl7qpivZmq4xERG71WKHDnbd3N13w7vv2kB6Tvv3w3ffwfTp8O23tj1u/vw2vbZpY/d2uejh1FOB4NjRY9Teup4mfy6j6V/LKJOw3T6gdm27cPXOO6FatZOLVI8ft2tZe/aEAwfswFCvXholzY6iYv5hc58X6bboE3bnCea1Ro8y+6omVNlah1mfBREaaq/rO+5I+7y6w+YTl3D2XkOhwYEs6d/YS9V72D//2A1Rly2z/0EMH659ScV16qYrIiIiaTiO/QW9Tx/buOXTT23X2nQbEu3caRePfvUVzJ0Lx47Z4aX777cBtHHjy253mzoCFRkdyzK/Gvxb40aKNK1CmVz/2dHWr7+GF16wt/LlbSht3Zoc9erx4IMBNGlit13s0weiomDCBKhU6bL/WsTXbNpERM+HYMkS5tZoyLONniDpvzD2TarGrB05ePJJGDo0/c9H4tMJouc77hPKlIFFi+DZZ213pmXL7ILvUB8a6ZVsQyOjkq3pGhOR7OLMdXGP1rqaaSNKMns2NG1qZ/Sd9buq49ihx3fftWtBU1KgXDnbmbRNG7j5Zu8NQ27bZucOz5hhw/DRoxAcbPfi6NQJp2EjPvkEnnrK5uThw+HJJ+2yVcmiHMcu/nzqKfsJytix7G3+AD16GiZPttO2338f6tY99ymy5Mjo6T77DDp1gjx5YMoUuPVWtyuSbOpcI6N6i/YwY8x5bx07dnS7RBERyWZO3/YixYE/lhakS0QRfliUwtix8P336QTRLVvstNv777ejoi+8YBuj/P237WpUr55358OWLGk7Kn3zDezZY6cIR0RAdDQ0boy5tTHtyy1m3Trb2Oipp6BJE9i82Xslihft2WOnonbsCLVqwa+/sqbag4TXNkydCi+9ZC/X8wVRgD7NwggMSHsdBwb406dZWMbV7k333Qe//GKbHDVtCq+9Zj9UEskkFEY9bNu2bSdv77///lnHRo0alebxSUlJbpQpIiLZSOq2F8mJAeyeUYvdM2sRUOggVbsto2vXM6blJifbLkHXXMPxhT/w1h3dqNhqOHUDGxDlFL3yTUU9IW9eOzI7fjxs3QqjRsFvv0H9+oQ+2oxZLy3jgw8gJgaqV7ejY9p+MQuZMwdq1LCj5MOHw/z5fLqkLDffbEfFFy2yTWUvZtZ4RK1QhratTmhwIAY7Ijq0bXXfal50IVdfDcuXw733wvPP22nue/e6XZUIoDDqcSVKlDh5Cw4OTnPsyJEjBAcHM2XKFBo3bkxgYCDvvfceEyZMICgoKM15Fi5ciDGG3bt3nzy2dOlSGjRoQJ48eQgNDaVr167s37/fmy9PRER8UHxCIkn78rBtfH0O/1mC4Aa/U/zBpezLccYvpCcCHU8/zY4a4TR7bCyjrmlBsp//yS6jUavi3HkR55I7Nzz9tB2xjYyElSsxN9/EY1/dQeyUldxwA3TpYmfzbt3qdrFyRY4csQ15mjaFAgVg2TKSevWlV2//kw1kV6ywDWUvRUStUJb0b8ymYS1Z0r9x1gqiqYKC4JNPYPRou+fq9dfbT2tEXKYw6oIBAwbw5JNP8ttvvxEREXFRz1m7di1NmzblzjvvZM2aNUyfPp3Vq1fz6KOPZmyxIiLi8wolF2LHlJtwjvtT4qElFLhpI8bPbocB2OGkwYPtdMfYWJg4kbatXmRj3iJpzpOYlExkdKwLr+Ai5MkDvXvDpk12KuLSpZS843rm5mvDlAG/smiRbcQ7caJGSX3SmjU2bY4caedgr1jBztBa3HabPfT003aHoeLF3S40EzMGunWDH3+0MyDq1oX33tN/EOIqn+um27On3ZDbm2rWtG90nvLUU09x1yXu+xQZGcm9997Ls88+e/LYu+++S61atdi5cyfFUjfMEhEROc3ff8M/k2vD8RSK3/czOYsdAE5bF7dsGTz22FmbMMb3n5Xu+TJ9l9GgIBgwwP7SPXIk5s03uW9GFC1b3M2TOwbRoUNVeg7bRd7GayhTyo8+zcKy5khYVpGSAiNG2OmlhQrZLYWaN+eXX2wfrd277QcMDz3kdqE+5MYbYeVKaN/etqHetMm2G84MU/Al29HIqAvCw89qJHVBK1asYPLkyQQFBZ281T2xKn/jxo2eLlFERLKATZugUSNIOZaDyI8SKF/l+Ml1ca83r0DEx5G2I+5//8HMmXZvlxMfbp4cNT3DuY5nOvnz2y42mzfDCy+Q78fvmLiqGp+XuJPCsXuIH1+fjWvyZM6px2L9+6/tQtWnj51nvXYtNG/ORx/Z2eT+/rBkiYLoZSlSxO4R/MQTdt1t794aIRVX+NzIqCdHKN2SN2/eNN/7+flx5hY7ZzY2SklJoVOnTvTq1eus84Vq3ygRETnD5s02iB44YKcv1qpVjF6c2Kpi9mx4sKntmJu6CWP+/Gme36dZGAOmryUxKfnkMZ/sMlqwIAwZAj168Em77rT9KYoIZxYTk9vT/7PX2d5wL6/nidXoaGYzdaoNSklJ8OGH8MgjHEsy9OgK//ufzahTpthMJZfJzw/GjoWAADv6nJRkm4FphFS8yOfCaFZUtGhRDh8+zP79+8l/4peB1WfMRb7uuutYv349lbSLt4iIXMCWLdCwIezfb7fkrFXrxB179sAzz9h5jWFhdu1YvXrpniM1nJ2+N6lPT2ktUoQXb36IkTXuoOvPX/Dwyk9o4zeDFxa+ysS41uzvdlYeFzfs3w/du8OkSbYT0aRJUKkS8fFw113w00/Qty+8+irk0G+xV84YG0BPD6RjxmiDXvEa/WecCdx4443kzZuXAQMG0KtXL9asWcPYsWPTPKZfv37cdNNNPPHEEzz++OPky5eP33//nZkzZ/Lee++5VLmIiGQ2//xjg+h//9kget11J+44McWRnTvt+rsXXrCdaM8jolao74bPdIQEBxIHDLm1M59e25yX5/6PMVu60+nPD+hSbQwvfleHa65xu8psLHXO7ZYtMHCgvUZz5GDxYrj7bjvKP3Wq3V5UPMgYeOMNG0iHD4fjx21jIwVS8QJdZZlAoUKF+OSTT5gzZw7Vq1dn3LhxDBkyJM1jatSowaJFi9i8eTMNGjTg2muvZcCAARRX2zgRETkhNYgmJNitGK+//sQdS5bALbfYPy9fDq+8csEgmhX1aRZGYIA/ABuLlKb9va/Qs+1zVAzewWf/1mXltY8w/X87Xa4yG0pKsut7b7nFBqPFi2HQIBz/HIwda6ebBwXBzz8riGYYY+x0/RdegA8+gEcftR13RTKYOXOtoreFh4c7MefY52jDhg1cffXVXq5IshNdYyKSVfz7rw2ie/bYIFq79ok7Zs2y8xvLlLFrRcuWdbNM10Wtijt76nHlAhzo/wq5x47gkJOHObe8QuvvniBnHk0gy3B//QUPPmg/JOnYEd5+G/Ll48gR6NoVJkyAli1h8mQ4sX27ZLTBg+3I9IMP2n8AzYcWDzDGrHAc56wurrq6REREfNzWrXb0aPfuM4LoxIl2hKNWLds5s2hRV+v0tHSD5QWmFZ9r6nG+0cNIerwjO+54irsXPcWfRT+g4KdjKNK6bkaVn705Dnz0EfToATlzwuef2w9NsCP87dpBTIwdMB04UDNGveqll2wAff55O2V30iQ7hVckA+g/bRERER8WF2eD6K5dduDzhhtO3DFiBHToYIdL58/PkkF0wPS1xCUk4gBxCYlXvE1LQPWrCNs8myU9PycwcQ9FIuqxvXlH2LHDY3ULdvi+XTvo1MnuefnrryeD6IIFdnr5+g0pXPXQr0w8Oov6r8/X9jve9txz8PrrdpHu/ffDsWNuVyRZlMKoiIiIj0oNojt2QHS0/b0ex4H+/eHZZ+0v+LNmQb58bpfqcZHRsWm2nQFITEomMjr2yk5sDHXfuosDy39nXKH+FIr+lCPlwnBGvW1HiTJY1Ko46g6bT/n+s6g7LAuGsDlzoHp1e12+8Yb9vlQpHMdu33fbbZArKIkSDy0hMeRfj33QIJehTx946y348ku7WPfoUbcrkixIYVRERMQHxcdD48awfbsNojfdhA1LnTrZjphPPAGffQa5crldaoaIT0i8pOOX6urwvNy3aSi9m67lhyM3Ynr2ILnW9XbT1gySEaO9mcaRI9CrFzRtavd+Xb7cfmDi58fhw9C+vb27VSsI7bCUlAL70zzdIx80yKXr2RNGj4YZM+xo9pEjblckWYzCqIiIiI/ZtcsG0fh4+P57uPlmIDHR/rL40Ud2kd3YseDv73apGSYkOPCSjl+O/Plh1PdhrIv8nrv9vmDb7/uhSRObmH7/3WM/J1WGjfa6LXXu7ciR8NRTdjHotdcCsGkT1K0LU6bYJs9ffgk7jxxM9zSe+qBBLlG3bnarl1mzICLCvteIeIjCqIiIiA85dszOvt2yBb77DurUwe7l0qwZzJwJ77wDgwbZrRqysNO3aUkVGOBPn2ZhHv05xsCzvQ3d57ejfpENPOc/nKNzF+FUq2aD1e7dHvtZGT3a63X//gv33ms/OUlMtBfs229DoP3AYM4cCA+3gfSbb2y/HD8/73zQIJeoSxf48EO7MP3OO+HwYbcrkixCYVRERMRHOA507w6LFtkB0Hr1gG3boEEDuwnjp5/aB2QDEbVCGdq2OqHBgRggNDiQoW2rX7Cb7uVq0ACW/5qbNc36UvrIn3xfugvO2LFQqZJd++iB9XRZJoQdOQKvvgpXXQVff223Clm/Hpo3B+x1/Prr9tuQEDtQevvtp57urQ8a5BI9+qjd6mX+fLvfzsH0R7BFLoXCqIiIiI8YPRref982urz/fmDjRptIN260Q0v33ed2iV4VUSuUJf0bs2lYS5b0b5xhQTRV0aJ28HnAiGK0jhtLk6K/sveqOrbRy9VXwxdf2KR1mbJECPvmG7jmGnjhBWjRwk5nfvHFk6OhBw/ay7RfPzur/KefbJ4/nbc/aJBL8PDDdquXRYvsv++BA25XJD5O+4yKiIj4gDlzbC+R1q1hyBBg9Wo7tHT8uG2qc+ONLleYPfj52UY7t9wC9913DUV/+ZYJ7WfTfvWzmLvvtgsgR4w4bY+di5cati5179TLcTl7tJ7Xn3/aC/Tbb20wnzPHrq89zV9/QZs28NtvtsdWnz7nnk1+rv1gJRN44AG77+j999uGVN9/DwUKuF2V+CjjXMEneJ4QHh7uxMTEpHvfhg0buPrqq71ckVyq7t27s27dOhYuXAhAx44d2b17N998881ln3PQoEF88cUXrFu3zkNVpk/XmIj4gj/+sFmzdGlYuhSCfo+BW2+1vwBGR9tf/sXrDhyAJ5+EyZOh0S3JfNnyIwqOeNHutfPAAzB0KJQp43aZZ0nt2nt6s6TAAP+LGn08M8T2r1+KVt+MtwE8Vy67Xvmpp2xYOc2338KDD9rw+dlnNsOIj5s+3a4Jvv56u5Y0f363K5JMzBizwnGc8DOPa5puBomLi6NLly6UKlWKnDlzEhoaSufOndm6dWuax3Xs2JE77rjjnOdZs2YNrVu3pkSJEuTOnZsyZcrQrl07tmzZktEv4bKNGjWKyZMnX9RjN2/ejDGGMz+Q6N27Nz/88ENGlCci4lMSEmy/kBw57PK7oK2/2+lxhQrBkiUKoi7Kl8/OWPz4Y1i+wp9Kwzsza+SfthPP9OlQpYqdU71rl9ulpnG5XXvTbD3jOFy39Htqt6gDw4bZUbI//oBnnkkTRP/7Dzp3tksMy5a160MVRLOItm3t1PQVK+wsDU3ZlcugMJoBNm3aRHh4OOvWrePjjz/mr7/+YvLkyaxfv57atWuzefPmizrPrl27uPXWWwkKCmLWrFn8/vvvTJo0iYoVK7J///4Ln+ASHDt2zGPnKlCgAMHBwVd0jqCgIAoXLuyZgkREfNTx43Z93d9/2y0vyuXYCk2bciQF7m33MuXH/ErdYfOzxj6UPuzhh2HlShu27rg/H0/vf4Uja2Jt2+OhQ6FUKTssuHjxFa0p9ZTL7dqbGmLDdm3msykDeGdmJLvyBNPlibdtY5sSJdI8PnX56EcfQd++dn1ohQqeehWSKbRuDVOn2n1jtYZULoPCaAbo1q0bfn5+zJ07l1tvvZUyZcrQqFEj5s6di5+fH926dbuo8yxZsoR9+/Yxfvx4rr/+esqVK0eDBg14/fXXqV69+jmflzra+sorr1C8eHGCgoJ45JFHSDxtX6iGDRvStWtXevfuTdGiRalbty4Av/32Gy1btiRfvnwUK1aM+++/n+3bt598XnJyMr1796ZgwYIULFiQnj17kpycnO7PT+U4Dm+++SaVK1cmV65clCpVigEDBgBQvnx5AGrXro0xhoYNGwJ2mm61atVOniMlJYUhQ4ZQunRpcuXKRfXq1ZkxY8bJ+1NHWL/88ktuu+028uTJQ9WqVZkzZ85F/V2LiGRGffvaWbhjx8It1+yBpk1J2ruPB9oOZFmOwjhAXEIiA6avVSDNQFGr4qg7bD7l+886Z/ivUsWGrR497O46N99bhtgXJ9suso8/bpNZ/fpQo4b9B/Xwh8qX4rK69h4+TO0l3zL5s+f57qOnqLL7HwY0607rh0cwp0DahLlnD7Rvb7djLVjQNnoePvxkDyPJatq2tXOvf/5ZXXblkimMetjevXv5/vvv6datG3ny5ElzX548eXjyySf57rvv2Ldv3wXPVaJECVJSUvjiiy+41LW9P/zwA2vWrGHevHl8+eWXzJ49m379+qV5zOTJk3Echx9//JGJEyeybds2brnlFqpVq8by5cuZO3cuBw8e5M477yQlJQWAN998k/fff5/33nuPn376ieTkZD755JPz1vLcc88xZMgQBgwYwPr16/n8888pXbo0AMuXLwfg+++/Z9u2bUyfPj3dc4waNYrIyEiGDx/O2rVradOmDW3btmX16tVpHvf888/z9NNPs2bNGmrXrs19993HQb0piogP+ugjeOstePpp6HTfQftL3t9/0+uBl1lZuHyax17MFEu5PGmmpnL+8J8rF4wcaTvu/vuvXUo3YXlVnFFvQ3y8bYWcMyd062b3NHniCduIyssuumuv49iA8fjjULIkI795k3L7tvF23fto1Pk9ptRsToqff5oQ+8UXULWqHSwbNMjO4IzLceEwLz7urrvs1lJLl9r3qkOH3K5IfITvddPt2dP7b9w1a9r/u1yEP//8E8dxztkUp2rVqjiOw59//skNF+i0d9NNN/Hcc8/RoUMHunXrRu3atWnYsCEPPvggZcuWPe9z/f39GT9+PEFBQVSrVo3hw4fz2GOPMXToUPLmzQvYUck333zz5HNeeuklrr32WoYPH37y2MSJEylUqBAxMTHccMMNjBw5kr59+3LPPfcANiRGR0efs46DBw/y1ltvMXLkSB599FEAKlWqxM033wxA0aJFAShcuDAlzpjec7o33niD3r1788ADDwAwePBgFi1axBtvvJFmfWqvXr1o1aoVAK+99hoTJ05k9erV1KtX77x/XyIimcnixTan3HYbvDn0GLRpB7/8Al9+yayfA9J9zoWmWMrlOd/6ynM1+7njDlizxs7MfeQRu6709dfzcn2nTvDYY/bf8t137WLT996Dm26Crl3hnnsgd+4Mf00X7NobH2+LnjDBbs0SGAh33cXienfQ5Z8gDh8/9QF5aojdvt1m7OnTbQifM8cOAp/ZLCk1zJ9eh2QR99wDKSn2wm/Vys4GOGNgRuRMGhnNIOYcvcpTRzjPdf+ZXn31VbZv3864ceOoXr06H374IVWrVmXevHnnfV6NGjUICgo6+f3NN9/MsWPH2Lhx48lj119/fZrnrFixgkWLFhEUFHTyljqCuXHjRv777z+2bdt2MkgC+Pn5ceN5thP47bffOHr0KLfeeutFvd707N+/n/j4+JNTiVPVq1eP3377Lc2xGjVqnPxzSEgIADt37rzsny0i4m1btthZb+XLw9QpKeR4rIPtVDluHEREXN4US7lsl7u+MjQUnnojjvKt/mDhT8cID4cGtx9m8xZjt30ZPx7i4mwX2r17oUMH+6Teve02KRnsrD1aqxaxw5otW9q2zf37Q+HC8MEHsH07TJxIvS738Fq7a9Ps//lam+rsXxtK1aowa5btZfTzzzaIwuU3SxIfdd99MHEi/PCDDaSHD7tdkWRyvjcyepEjlG6pXLkyxhjWr19PRETEWfdv2LABYwwVK1a86HMWLlyYu+++m7vvvpuhQ4dSq1YthgwZckUBDzg5QpoqJSWFli1b8sYbb5z12OLFi5+cqnspPLl1UHoB/sxjAad18Eu973LqFhFxw8GDth/IsWPw9QyHggOftmuxhg2zI2rYKZbpbctx1hRL8YiQ4EDi0gmeFwr/UavieGHGWlKqJhNacRP//VyRH+eUp3IVhx5PG557DgoVKmQ3Le3ZExYssKOlI0fCm2/aLsnVqqW9VawI/v7n/bmXJDnZzjabMMFOsdy71wbi/v1tOK5S5aynnL7/57//2hm8330HderYqeVhZ1yGlxvmxYc9+KAdIe3Qwb6hff21FgzLOfleGM3kChUqRLNmzRg7diy9evVKs2708OHDjBkzhhYtWtj/AV2GnDlzUrFiReLj48/7uLVr13Lo0KGTgfPnn38++dxzue6665g2bRply5ZNE+pOV7JkSX7++WcaN24M2LC5fPlySpYsme7jq1atSq5cuZg3bx6VK1dO9/UAZzVBOl3+/PkJCQlh8eLFJ38uwOLFi6lateo5nyci4ktSUmxX1rVr7S/3YVMHw5gx8OyztpPRCRecYikedbnh//QRQb9cxynYIJZ8tbaQtPxqRowI4cMP7Q4w3btD7twGGje2t/h4O3132TLboveLL0514M2d+1RIrV79VEgtVcpu4Hm6I0fsyOvWrWm/nv7nbdtsIM2VCyIi7JziJk0uGHhTUuzy1z597NNHjbJTdNN72uWGefFxDz1kL5RHHrHX1owZXpmCLr5HYTQDjB49mjp16tCkSRNeeeUVKleuzMaNG3n++edxHIfRo0enefz+/fvPasQTHBzMunXr+Oyzz7jvvvuoUqUKjuMwc+ZMvv32W15++eXz1nD8+HEeffRRXnrpJeLj4+nfvz+dO3c+azT0dN26deP999/n3nvvpV+/fhQtWpS///6badOm8eabb5IvXz569OjB0KFDqVKlCtWrV2fs2LFs27btnGE09TkDBgwgV65c3HLLLezZs4cVK1bQtWtXihUrRmBgINHR0ZQrV47cuXNToECBs87Tp08fXnrpJSpXrsz111/P5MmT+fHHH1mxYsV5/x5ERHzFoEHw1Ve2aVHTv8baAx06QGTkWUHj9NEpyViXG/7TG/nLkf8IAU1WsfqjEPr3t2HunXfg1VfhgQfAzw/b2OhEx3nANoLZsAHWrTt1mz/frulMlT+/DaUFCpwKm3v3nl1Uvnw2uIaG2tAZGgqVKtnRq4IFL/h3kZJi14IOHWpnYTZubEPp+bZr0Uh+Ntahg/204rHHoE0b+wanQCpnUBjNABUrViQmJobBgwfz0EMPsXPnTooWLcrtt9/O1KlTKVWqVJrH//jjj9SqVSvNsXbt2vH6668TFBRE7969+ffff8mRIwfly5fnjTfeoEePHuetoUGDBlxzzTU0atSIw4cPnzzf+YSEhLBkyRIGDBhA8+bNOXLkCGXKlKFp06bkypULgGeffZbt27fTqVMnAB566CEefPBBNmzYcM7zDh06lIIFCzJkyBC2bt1K8eLFefjhhwHIkSMHb7/9NoMHD+bll1+mfv36LFy48KxzPP300xw4cIC+ffuyY8cOwsLC+PLLL6lZs+Z5X5OIiC+YOhWGDLG/s/Uo/hk82N2ut/rgg7NHvMTrLif8n29EsEYN+PZbmDfPBtKHHrJLRyMj4awVOHnzQni4vZ1u3z67bUxqQF271q7tLFvWzplNDZ2pX0NDbWi9DPHxdgruBx/YNc1Fi9olzJ06Xfjy1Eh+Nvfoo/ZTjM6doV072+HqxO+UIgDGk2v6Lkd4eLgTExOT7n0bNmw4Z1daObeOHTuye/duvvnmG7dLyfR0jYmI21auhHr1bAfS+f1nE9DmDttdNTpa66x82JldZMGOCA5tWz1NEEtJgSlT7JTdLVugeXO7J+dp/fhckZxsL8Fx42xT1ORkO5japYsdSD2xykbk4owbZxcYt2plp57rAsp2jDErHMcJP/O4RkZFRERcsnu37ZxbpAjMeG4ZAXe1sZs0quGHz7vYEUE/P9vvpV07u0T4lVfsjnKNG9vwd+utcN11nu1bdD5bt54aBf33XyhWzDb47dzZ9k8SuSxduthPNJ580m4BM22aAqkACqMiIiKuSE6G+++3MytjJv5Gofa3Q4kS8P33EBzsdnniAZcyvTd3btur6pFHbDPdr78+tXQ0OBgaNrTB9NZb4aqrPDt7OznZNs0aN85uz5KSAk2b2qnDd96pzCAe0rWrvbi6d7dbwEydCudomCnZh8JoFjRhwgS3SxARkQt48UWYOxfqNZ9HoU73scdxWPHWRJqWKOF2aeKiQoVsQ6NXX4UdO2yvonnz7C0qyj4mJORUML31Vrss9GIcOwZ79py67d1rd3b56CM7Ilq8OPTrZ9eCnq8pkchl69bNBtKnn7afxk2ZokCazSmMioiIeFlUlO1IWqb6Oj785VGCjh7m3geHsfmXAwwtHafmLgLYcHj//fYG8Pffp4Lpd9+daqhbpYoNpZUr24B5etg8PXweOnT2zzAGmjWz27O0aqVcIF7w1FN2OL5XLztH/dNPIYciSXaV6f/lHcfBqJOgZAC3m3eJSPYUG2v3Ey0YupspCQ9Sbl88He55mQ3FKkBSMpHRsQqjkq4KFeytc2c7uLR27alwOmkSHDxo16AWLAiFC9tbaKjdkjT1+0KFTv25cGE7ylqsmNuvTLKdnj3tRfzss/brpElaJ59NZeowGhAQQGJiInny5HG7FMmCEhMTCdBHwCLiRQcP2oZFuXI6/C/v/dT541d63vEsP5W99uRj0tufUuRMfn5w7bX29swzkJRkr68CBU7sVyqS2T3zjL1Ye/Wy+wfNmGH3DZJsJVO/XRUrVoy4uDgOHz6sUSzxGMdxOHz4MHFxcRTTx8Ei4iWOY7fc+/13WN7iJdr9MZfI+g8RdU2jNI8LCdbogFy6gAA7IqogKj6lZ0+71cuqVXZLq99/d7si8bJMPTKa/8TmzPHx8SQlJblcjWQlAQEBFC9e/OQ1JiKS0UaMgM8/h+/afUD5ya+wuc39fFT1fjiecvIxgQH+9GkW5mKVIiJe1q6d7cJ1551Qpw589RU0aOB2VeIl5mJGHI0xtwC9geuBEOARx3EmnHa/AQYCXYCCwDKgm+M46y907vDwcCcmJuayihcREfEFCxbYPSNfvvl7nv/5DkyTJjBzJlHrdl5wH0oRkWxh0yZo2RL++gs+/BAeesjtisSDjDErHMcJP/P4xY6MBgHrgIknbmfqCzwLdARigZeAOcaYMMdxDlxWxSIiIlnAv//CvfdC6zKreH7N3Zjq1e0QaUDAJe1DKSKSpZUvD0uWwF132S5vGzfCwIGe3VRXMp2LWlngOM63juM85zjOF0DK6fedGBXtCQxzHOdLx3HWAR2AfMADHq5XRETEZxw9an+vKnL4H6YeaokpWBBmzYJ8+dwuTUQk8ylY0O5b1LEjvPyyDaVHj7pdlWQgTyxzLw+UAGanHnAcJxFYBNTxwPlFRER8Uo8eELs8gZ8K3U7A0UPw7bd2Lw0REUlfzpzw0UfwyisweTI0bWo3zZUsyRNhtMSJrzvOOL7jtPvSMMZ0McbEGGNidu3a5YESREREMpfx4+Gj944RU7YdBbb/YZtyVKvmdlkiIpmfMfD88/Dpp/Dzz3DzzXYtqWQ5nmwAfmYnJJPOMftAxxnnOE644zjhRbWfkIiIZDErVkDXJxxmlehEpS3zbTOOxo3dLktExLfcfz/Mmwd79thAunSp2xWJh3kijG4/8fXMUdBinD1aKiIikqXt3m13KhiWeyC3bZ8EQ4aoK6SIyOWqVw9++smuJ23cGKZOdbsi8SBPhNFN2EB6W+oBY0xuoD6gjy9ERCTbSE62H+Q32/ohPfcPgcces1PNRETk8lWubANp7dpw333w2mtwEdtTSuZ3UWHUGBNkjKlpjKl54jllTnxfxrEblY4E+htj2hpjqgETgIPApxlTtoiISObzwgvgNzead53HoVkzePddbUsgIuIJhQvD3LnwwAP2Q75OnSApye2q5Apd7D6j4cCC075/+cTtY+zeoq8DgcAYoCCwDGiqPUZFRCS7mDgRvhu2mp8D7sLvmlN7iYqIiIfkymU77FasaJdAbNkCX3wBwcFuVyaXyTguD3GHh4c7MTExrtYgIiJyJX78Ebo0/oslfvUoWCwAs2yZtnAREclIEyZA585QoQKMGQNNmrhdkZyHMWaF4zjhZx73ZDddERGRbOfvv6HrnXHM4TaCg45jZs9WEBURyWgdO8KcOXD8ONx2G7Rure1ffJDCqIiIyGX67z94sPkevtjflJK59uAX/T1cfbXbZYmIZA8NG8L69TB0KMyfD1WrQt++sH+/25XJRVIYFRERuQzHj0OHtgd4+88WVPbfiP83X0P4WTOQREQkI+XODf37wx9/QPv28MYbtvvuhx/aFueSqSmMioiIXIbe3Y/w9PzWXO+3Ev8vP7ef0IuIiDtKloSPPoLly6FSJdttt3Ztu6hfMi2FURER8bqoVXHUHTaf8v1nUXfYfKJWxbld0iUZM+o4Dd+7j8YswO/jCdCqldsliYgI2BkqixfDlCmwaxfccgvce6/tvCuZjsKoiIh4VdSqOAZMX0tcQiIOEJeQyIDpa30mkH7/bQr5ej5GBDNIGfWOnRYmIiKZhzFw330QGwuDBsHMmXDVVfDSS3DokNvVyWkURkVExKsio2NJTEq7jicxKZnI6FiXKrp469c5bGrTi4eZyNEXBuP3dHe3SxIRkXPJkwcGDrShtE0buzdpWJjdqzQlxe3qPGPfPjs12UcpjIqIiFfFJyRe0vHMYtcumFN/MF2Pvc2BTr3INfgFt0sSEZGLUbo0fPqpnb5bogQ89BDUreubIS4xEebOhQED7JrYIkWgRQufDdcKoyIi4lUhwYGXdDwzOHoUPrnxbXomDGJXy47ke+8NOw1MRER8R2oAHT8eNm+GG2+0Sy2iouwnjpnR8eOwbBm89hrceisULGj3VX3jDdtJ+KWXYMYMt6u8bMZxHFcLCA8Pd2JiYlytQUREvCd1zejpU3UDA/wZ2rY6EbVCXawsfY4D4+pN5PGlHdhauw2llk6DHDncLktERK7EgQM24I0cCUeO2GNhYVCv3qlbxYre/+DRcWDDBpg3z94WLrSbWgNcey00aWJDaf36EBTk3dqugDFmheM4Z+1/pjAqIiJeF7UqjsjoWOITEgkJDqRPs7BMGUQBPm8/gzaftOOfCg2psP4b+0m0iIhkDUeOwIoVdgrv4sWwZIldhwlQrFjacFqzJgQEeL6Gf/89FT7nzYNt2+zxChVs8GzSBBo1gqJFPf+zvURhVERE5BItenkBNwxqwdbCNam4aS4m37k/hfalgC0iIueQkmJHJpcsORVQN22y9+XJAzfdZKf71qtn/5w/f9rnHzsGe/fCnj3nvp1+/+7dsHOnfW6xYtC48anRz3LlvPrSM5LCqIiIyCX47eNfKN2xMbsCyxLy1yJyhxQ652N9beqxiIhcgvj4tOF09WobWv38oFo1u3QjNVwePHju8+TKBYULn3276iobQKtVy7L9CBRGRURELtKmSYsJ7nAnB/yDybNiMUVqhJz38XWHzScunW7AocGBLOnfOKPKFBERNxw4YJsKLV5sv/r5pR8yCxeGQoVO/TlPniwbNi/kXGFUHRhEREROs/mtryjxzAPE+ZeB76Ipc4EgCr67XY2IiFyGfPnsSGaTJm5X4vO0tYuIiMgJWwa8S+ln7uK3nDUxS5ZQqUm5i3qeL25XIyIi4jaFUREREcchruPzlB32JAsDb6fwqnlUvLHIRT+9T7MwAgP80xwLDPCnT7MwT1cqIiKSZWiaroiIZG9JSWyPeJzQb8fzWf7O1Fk1ljIVLu1/j6lNitRNV0RE5OIpjIqISPZ18CC7b72HEsu/450ig2i3+iVCQs1lbdMSUStU4VNEROQSKIyKiEj2tHMn/9VrScE/VzIoZBxdV3amePGzt2mJS0hkwPS1AAqbIiIiHqQ1oyIikv1s3MjBmnUJ+HM9z1aIovsaG0TBTrU9fb9QgMSkZCKjY10oVEREJOvSyKiIiGQvMTEcufV2ju5Ppu/V83hjyc0ULHjqbm3TIiIi4h0aGRURkewjOpqkeg3Zvj8P3Wst5a2f0wZR0DYtIiIi3qIwKiIi2cPEiaS0vIP1Ryvx7M0/8f6iMPLnP/th2qZFRETEOxRGRUQkS4tauZX/Ne8EHTowP7kB3cJnM2luSYKC0n98RK1QhratTmhwIAYIDQ5kaNvqal4kIiLiYVozKiIiWdbsmUsp1qM7EZtW8Sn307XCcEo0Wc3s2GvOGy61TYuIiEjGUxgVEfGiy9m/Ui7D8eMwahS39H+eYykBPMkYPq4cQZHWqznqOERGx+rvXURExGUKoyIiXuA48P6323n5g90c2FoOHENC3qN0XbyP5U1y0+rGwhQrBsWKQVAQGON2xT5s9Wro1AlWrGCOf0u6OmM5WPcIRW5ejfF3AHXGFRERyQwURkVEMkBiIqxYAT/9dOq2fXsJoAQmRzL4peAcCwBg6CwYetpzc+e2obRoUU4G1GLFoEoVuO02KFvWlZeU+SUmwuDBOJGR7A8oTCem8W3RZhRovobg4vvTPFSdcUVERNynMCoicplSp9zG7UukMAVpWPAqjm8rxE8/2cG5pCT7uIoVoUkT+HbHOnKF7iOg6AGMn4Nz3I/kwzlJOZyTD++tz86dnHXbsQPWrbNfjx2z5wsLg6ZN7a1hQ87ZiCdbWbAAunSBv/5iSuCj9DwWSbdBhbinRRwvzTxEYtKph6ozroiISOagMCoichmiVsXR63+b2LnoKo5uLcSWg7lZCeTKncJNN/rx7LNw881w0012VBOg7rCdxJ02PdTkSCFH/iOEljG0aHH+n+c4sGEDzJ5tbx98AO+8AwEBUKfOqXBaqxb4+5//XFnKvn3Qty988AE7girwAHPZU+VWZk+AmjUBQgkIQOt0RUREMiHjOI6rBYSHhzsxMTGu1iAicin274erbv+HbUtL4ZfrOIEVdpErdB85Q/ZRrvJxfnq+UbrPi1oVx4Dpa0lMSj55LDDA/7K2DTlyBJYuPRVOV62yxwsXtqOwTZvaKb2lS1/2y8zcHAemT4fu3UnZuYuxuZ/luaMDefbFPAwYADlzul2giIiIpDLGrHAcJ/zM4xoZFRG5SI4DU6bAs8/C9u2lCar5D8G3xOIfeGoO6PYD535+auD0xChd7tzQuLG9DRtmp/TOnQvR0TacTp1qH1e1KkREQJs2cP31WaQxUlwcdO8OUVFsLliLtimzcKpcx48T4Npr3S5ORERELpZGRkVELsJvv9n8s2ABhIfD0Rt+YX++nWc9LjQ4kCX9G3vuBycn8/XyTYyO/o0tB1MoUiT/BQOs49h1prNnw6xZsGgRJCfbUdKICGjbFurVgxw+9nHk/Klz2Db2A1r9/A05k5N5LdfLvH70WQa8mIPnnrNTlkVERCTzOdfIqMKoiMh5HDwIgwfDW29Bvnzw2mvQuTPM/PUyptz+/Td89x3MmXOqI9GFbikpJ5+egmFbviJsLVSS4rWqUu6GGrY7UuotODjdH7tnD8ycCV99ZUdOjx6103nvvNOOmN52mx1pzZS2boVPP+W/98dT4K/fSTL+fBfUhGcOjGZr8WIMHZVIr3uLu12liIiInIfCqIjIJXAc+OIL6NXLzgp97DEYOtRut5IqtZvuOafcJibaYcnvvrO3P/6wxytUsOExZ860t1y5zjr2/s9b2ZMESf45CDp6mDIJ2yibsJ3y/22n8MF9aYsuVOhUMK1U6dTXatVOBtWDB+H7720w/eYbu/41KAhatLDBtGVLyJ8/Y/9uL+i//+DLL2HyZFi4EByHmAK1+Dj5QaYc7MAev0IUuPkvCtz8F6UK5/bsSLSIiIh4nMKoiMhFio2Fp56yA5g1a8LYsbYz7kX5669T4XPhQhtIc+e2e7C0aGFvlStfdC3l+88ivXdpA2x6oYEdbd248dTtr7/s13/+sXNzU5UqBTVqQPXqJ2/HKlzF/MU5+eormDHDDtYGBMCtt8Itt8CNN9opyV4Jp8eO2b+zyZNxZs7EHD3K7kKV+cy/PW/tepC/qUiu0L3kCdtGnrDt5Mh/5NTfw7CWXihQRERELpcaGImInEPqCOfWncdIXnUVOxaXJW9ew+jR8MQTF9gq5fBhGzq/+84OOf71lz1eubKdz9uiBTRoAIGBl1VbSHBgmu1gTj9OUJANmDVqnP3EpCTYssWOxq5de+o2Z87JDVBz5shB87Awmlevzrvdq/N7QHW++qs6E38oy/ff205HxsDVV9tgesMN9mv16h5ab+o4tiXw5Mk406Zh9u7lYJ6iTA98nHeOtmfF3nBuaWDodRdMiFvMbue/9P8eRERExCdpZFREsrXU7Vb2bcrP7pk1Sd6fh/w14hj1ph8dm5Q89xOPHYM334RXX4VDh2zYbNTo1OhnxYoere9St4M55xTipKSzA+qvv9rgmiooiOQSofyXuzjbnOJsPFiCtTuL83diCXZQnP9yFado9RJUqlOM8Lq5uOEGKFv2tE69KSmwezfEx8O2badup33vxMfD9u2YY8c4liOQ73K14X+H2rPArwn1Gwdw11222VLx4lf29yAiIiLu0zRdEZF01B02n43rcrPz8xvwz3uUwi3WkLv0vvN3xZ03D7p1s/N527aFxx+381ozqAvQBdempvP4Sw5u+/fbFrxr19rWwdu22Xm727fbr/+dPSoJsI9gtlOCfQHFCPJPpNjxeAof30EAx8967B4Ksd2UZJtTkjhC2EZJ1lKdb/3v5Kbb8nHXXdC6NRQp4pm/BxEREckcFEZFRNJRsv1Sdnx+A/5BRyh+/8/kCDoKnGMtYnw8PPOM3cSzYkUYPRqaN/d+0RdQd9j8dKf2XtG2M0eO2FCaetu+nePxO9i7fjv7/9pBctwOEk0eEvKUZH+ekuwPCuFAUEkO5AvhUP6SHM5fAnLnJkcO0tzKlIE77oCCBa/wRYuIiEimpTWjIiJnWLwYOyJ6RhCFM9YiHj8O77wDAwfa6bmDBkG/fpl2P5T4dILo+Y5flNy57VzcsmVPHsoBFDtxExEREblUCqMiki0tXmyXdpYo6ZAn4heScp8KooEB/vRpFnbqgU8+aaev3n47vP22x9aDZpTzNj0SERERyST83C5ARMTbUoNoSAjE/BRAZIcqhAYHYrBTWYe2rU5EaAA88gjUrw8JCTB9ut2YM5MHUYA+zcIIDEjbAjhNwBYRERHJBDQyKiLZyulBdOFCKFkSIkqGnmqEk5wM48bBc8/ZLrn9+8MLL0DevK7WfSlSX4ua/YiIiEhmpjAqItlGekE0jZgY6NrVfm3UCMaMsZts+qCIWqEKnyIiIpKpaZquiGQLFwyiw4fDDTfA1q3wySd2+xYfDaIiIiIivkAjoyKS5Z03iDoODBhgw+i998J770GBAm6VKiIiIpJteGxk1BiTzxgz0hizxRiTaIxZaoyp7anzi4hcjvMG0ZQU6NbNBtHHH4dPP1UQFREREfEST07T/QBoBnQAqgOzgbnGGC1aEhFXnDeIHj8OHTrAu+9C3772q59WLoiIiIh4i0d+8zLGBALtgP6O4yx0HOcvx3EGAX8BXT3xM0RELkbUqjjqDptPyfZLaXDrcfIXTjo7iB49CnffDZMnw6uvwrBhYIxbJYuIiIhkS54aBsgB+ANHzjieCNTz0M8QETmvqFVxDJi+lo3rcrPj8xvwCzpCnojFLNsed+pBhw7BHXdAVBS8/bbdwkVBVERERMTrPBJGHcc5APwEvGCMCTXG+Btj2gM3A2f2rMQY08UYE2OMidm1a5cnShARITI6lgP7crDrq+vxz3uE4vf/TFLuw0RGx9oHJCTAbbfB/PkwYQI89ZSb5YqIiIhka55cIPUQkAJsBY4CTwNTgOQzH+g4zjjHccIdxwkvWrSoB0sQkewsbm8ie76piXMsB0XbriBH0FEA4hMSYedOaNjQ7iE6bZpdLyoiIiIirvHY1i6O42wEGhhj8gL5HcfZZoyZCmzy1M8QETmflFVXc2RLEQq3WEPOIgdPHq/JAbjlFvjnH5g5E5o1c7HKrCdqVRyR0bHEJyQSEhxIn2ZhRNRS7zoRERE5P4/vM+o4ziHgkDGmILa7bl9P/wwRkTMtWgRx88uTv1o8eatvPXn8qgPb+WT6QDh8AGbPhnpaxu5Jqet0E5PsJJi4hEQGTF8LoEAqIiIi5+XJfUabGWNaGGPKG2NuAxYAscB4T/0MEZH07N4NDzwAFSsa/vculCoYiAHqH9lG1NQB5Dl+FBYsUBDNAJHRsSeDaKrEpORT63RFREREzsGTI6MFgKFAKWAv8CXwvOM4SR78GSIiaTgOdOwIu3bBzz9DrVoh3F8vBJYvh+YPQWAgzJ0LV1/tdqlZUnxC4iUdFxEREUnlyTWj04BpnjqfiMjFeOstmDUL3nkHatU6cfCHH+z2LcWK2SBavryrNWZlIcGBxKUTPEOCA12oRkRERHyJJ7vpioh41fLl0L8/tGkD3bqdOPjnn9C6NZQuDT/+qCCawfo0CyMwwD/NscAAf/o0C3OpIhEREfEVHm9gJCLiDf/9B/fdByVLwocfgjHAgQM2iObIAd9+CyEhbpeZ5aU2KVI3XREREblUCqMi4nMcBzp3tju1/PgjFCwIpKTAQw/BH3/Yrrnlyrldpk+6nG1aImqFKnyKiIjIJVMYFRGfM24cfP45DBsGN9984uCQITBjBowcCY0bu1mez9I2LSIiIuJNWjMqIj7l11+hRw9o1gz69DlxcMYMGDQIHn4Ynn7azfJ8mrZpEREREW9SGBURn3HoENx7r52WO3Ei+PkBv/0G7dtDeDj8738nFo/K5dA2LSIiIuJNCqMi4jO6d4fYWPjkE7trCwkJtmFRnjzw1Vd2T1G5bOfajkXbtIiIiEhGUBgVEZ8waRJMmAAvvHBiSWhyMjzwAGzeDF9+CaVKuVyh79M2LSIiIuJNamAkIplebCx07Qq33AIvvXTi4AsvwHffwbvvQr16rtaXVWibFhEREfEm4ziOqwWEh4c7MTExrtYgIpnXkSNw002wdSusWQOhocC0aXbxaJcu8N57bpcoIiIiIudhjFnhOE74mcc1MioimVrv3jaEfvPNiSC6Zg088gjUqQPvvON2eSIiIiJymbRmVEQyrR9+gDFjoFcvaNkS2L0bIiIgONiuE82Z0+UKRURERORyaWRURDKlo0fhiSegfHl45RXg+HE7NTc+Hn78EUqUcLtEEREREbkCCqMikim9/jr8/jt8+63duYVn+sL8+TB+PNxwg9vliYiIiMgV0jRdEcl0/vwTXn0V7rkHWrQAJk6Et96Cp5+Gjh3dLk9EREREPEBhVEQyFceBJ5+EXLlg5EggJsZ2zW3YEN54w+XqRERERMRTNE1XRDKVKVNg7lwYPRpKBiZA27Z2fei0aRAQ4HZ5IiIiIuIhCqMikmns22c759aubZsX8cjTtmHRTz9B0aJulyciIiIiHqQwKiKZRv/+sGcPREeDf9SXMGkSDBxo06mIiIiIZCkKoyKSKSxdCuPGwTPPQM0S26HJ43D99fD8826XJiIiIiIZQGFURFyXlASPPw6lS8PLgxx4sAscPGi76GqdqIiIiEiWpDAqIq4bMQLWrYMZMyDoiwkwc6Y9WLWq26WJiIiISAZRGBURV23aBC+/DBERcGeNzVCjBzRoAD16uF2aiIiIiGQghVERcY3jQLdu4O8Pb49MgQ4d7R0TJoCftkEWERERycoURkXEK6JWxREZHUt8QiIhwYH0aRZG0l+hfPednZFb+qu34Ycf4MMPoVy5cz4nolaouy9ERERERDzCOI7jagHh4eFOTEyMqzWISMaKWhXHgOlrSUxKPnksZ3Iudk9oRLnS/iyf8Bs5brgOmja1C0eNSfc5gQH+DG1bXYFURERExIcYY1Y4jhN+5nHNgxORDBcZHZsmVAJsn1+JfXv8GDcmiRyPPgz58sH774Mx53xOYlIykdGxXqtbRERERDKOpumKSIaLT0hM8/3R+AIcWFmWfNdtIXz2x7BiBXzxBRQvfs7nXOi4iIiIiPgWjYyKSIYLCQ48+WcnxbAnujr+QUdpc300DBkC7dtDu3bnfM7FHBcRERER36IwKiIZrk+zMAID/AE4EFOOpJ0FKNtkFaPnvAElSsA775z3OakCA/zp0yzMKzWLiIiISMbSNF0RyXCpDYdembqZfxZXIThsN9F5PiLf5r9g9mwIDj7nc9RNV0RERCRrUjddEfGaiAibPf/+cAElHmhsNxkdPdrtskREREQkA6mbroi46ptv7K4tr/bbT4n+HaFyZRg+3O2yRERERMQlmqYrIhnu2DHo1QvCwuDpTT1h61ZYsgTy5nW7NBERERFxiUZGRSTDvf02/PUXfHLv1/h/PB7694ebbnK7LBERERFxkdaMikiG2rHDzshteeNupvx6DZQsCcuXQ86cbpcmIiIiIl5wrjWjmqYrIhnq+echMRHey9ML9u2DuXMVREVERERE03RFJOOsWAEffQRj7/ye/F9PhgEDoHp1t8sSERERkUxA03RFJEM4DtSrB/F/HGRjYDX88gbC6tWQK5fbpYmIiIiIF2maroh41WefwdKl8GuTF/GbuwUWL1YQFREREZGTNE1XRDzu0CHo2xceDltGtXmj4MknoW5dt8sSERERkUxEYVREPG74cNix9RhjkzpjQkJg6FC3SxIRERGRTEbTdEXEozZvhshImFQ9krxr18KMGZA/v9tliYiIiEgmo5FREfGovn0hzPmde2IHwz33wJ13ul2SiIiIiGRCGhkVEY/54Qf44vMUNpfpjDmQF95+2+2SRERERCSTUhgVEY9IToYePaB/wXGU+Wex3WC0eHG3yxIRERGRTErTdEXEIz74AHatiePlI32hcWPo2NHtkkREREQkE9PIqIhckqhVcURGxxKfkEhIcCB9moXRoFwozz/nEFX4SXIcSoJx48AYt0sVERERkUxMYVRELlrUqjgGTF9LYlIyAHEJiQyYvpbKfxeg0d7vqcfX8PrrULGiy5WKiIiISGbnkWm6xhh/Y8wQY8wmY8yRE19fMcYo7IpkIZHRsSeDaKr92wJZPOUoHwR2h1q1oFcvl6oTEREREV/iqbDYD+gGdADWAjWAj4GjwBAP/QwRcVl8QmKa7x0H9s6/mtGmD/mP7YYPv4Mc+gxKRERERC7MU7811gFmOo4z88T3m40xXwM3euj8IpIJhAQHEndaIE3cWIybN63jMcZDv352ZFRERERE5CJ4qpvuYqCRMeYqAGNMVaAx8K2Hzi8imUCfZmEEBvgD4CQbEueV532/zhwsVRYGDnS5OhERERHxJZ4aGR0O5AN+M8Yknzjvq47jjE3vwcaYLkAXgDJlynioBBHJaBG1QgG7dnRDdEleSIikIn/DxPkQGOhydSIiIiLiS4zjOFd+EmPuAyKBPsB6oCYwCujjOM6H53tueHi4ExMTc8U1iIj3bN8Od1dcycLEG/B/tKPdZFREREREJB3GmBWO44SfedxTI6ORwBuO43x24vu1xpiywADgvGFURHzPC/2P8/bhTjhFikBkpNvliIiIiIgP8lQYzQMkn3EsGc+tSRWRTGLZMij08QhqsQre/RwKFnS7JBERERHxQZ4KozOB/saYTdhpurWAZ4CJHjq/iGQCKSkQ2fkPJjGQpJatCWjXzu2SRERERMRHeSqMPoXdT3QsUAzYBrwPDPbQ+UUkExj/YQpPre2Mf55cBIwbC8a4XZKIiIiI+CiPhFHHcQ4APU/cRCQLSkiA2Gfe4zEW4bz9IYSEuF2SiIiIiPgwrekUkYsy8pl/ePFgX/bf2ATz6CNulyMiIiIiPs5T03RFJAtbt9bhxvFPkDNHCrmmjNP0XBERERG5YhoZFZHzchz4+t5PaMF3JA16DcqXd7skEREREckCFEZF5Ly++Wgnj2/owbbyNxPUv7vb5YiIiIhIFqFpuiJyTocPA08/RT5zEL+oD8Df3+2SRERERCSL0MioiJxTVMcoWh2eRvyjL5KjRlW3yxERERGRLERhVETStXl1Ag0/f5ItwTUo924/t8sRERERkSxG03RFsrGoVXFERscSn5BISHAgfZqFEVErFIA/I3rTiJ3s/XQmBAS4XKmIiIiIZDUKoyLZVNSqOAZMX0tiUjIAcQmJDJi+FoBSs3/nti0f8tMt/bi5xfVulikiIiIiWZSm6YpkU5HRsSeDaKrEpGRGfbWG4i92ZlNAZa6bMdCl6kREREQkq1MYFcmm4hMS0z1+z/iplE7axLYhH5IrONDLVYmIiIhIdqEwKpJNhaQTNGv8tZHHt07im7JPUqdffReqEhEREZHsQmFUJJvq0yyMwIBT+4bmPJ7E69+MZiuluOqrYS5WJiIiIiLZgcKoSDYVUSuUoW2rExociAGeWfgNVx39k3n3jKNSrXxulyciIiIiWZxxHMfVAsLDw52YmBhXaxDJ7pJXrsG5PpyvAh+gxc6PCQpyuyIRERERySqMMSscxwk/87i2dhHJ7o4fZ2+bx0ihEH4jRyiIioiIiIhXaJquSDZ3+NURFP1nBaOvGkPbzoXdLkdEREREsgmFUZHsbNUqcg5+gem05a4p7TDG7YJEREREJLtQGBXJrg4f5nCbB9iRUpRfOo/j2ppKoiIiIiLiPVozKpJFRK2KIzI6lviEREKCA+nTLIyIWqHnfPzxp58h95ZY+hafw3sjND1XRERERLxLYVQkC4haFceA6WtJTEoGIC4hkQHT1wKkH0ijosjx4Xu8Th8emXyrmhaJiIiIiNdpmq5IFhAZHXsyiKZKTEomMjr27AfHx5PUsRMruI6NHV+hSRMvFSkiIiIichqNjIpkAfEJiRd3PCWFlPYPc3x/Ir2KfcrXb+X0QnUiIiIiImfTyKhIFhASHHhxx0eMwG/BPJ52RtLngzCCgzO+NhERERGR9CiMimQBfZqFERjgn+ZYYIA/fZqFnTqwciUpA54jyrTh8P2daNXKy0WKiIiIiJxG03RFsoDUJkXn7KZ76BDOAw+wyxSjX6H3WfK2tnEREREREXcpjIpkERG1Qs+9lcszz8Aff3C/M5fBYwpTpIh3axMREREROZPCqEhW99VXMG4cI/z7UqBVY+65x+2CREREREQURkWytrg4nE6diA26nuH+Q1gzFoxm6IqIiIhIJqAwKpJVpaRAhw4cP3iE1sc+4fXxOSlZ0u2iREREREQshVGRrOrNN2HePHrmfJ/yzcLo0MHtgkRERERETlEYFcmKVq7Eef55Fhdry8RDj7HuPU3PFREREZHMRWFUJKs5dAjuv5/DQcVovfN9ho8xlC3rdlEiIiIiImkpjIpkNb164fz5J/flmUf1WwrxxBNuFyQiIiIicjaFUZGs5Kuv4P33mV6pH3O3NuLXD8DPz+2iRERERETOpjAqklXExsJjj7G3/PXc/9dgXouEypXdLkpEREREJH0aMxHJCnbvhpYtSfHPQdOEadSsnZOePd0uSkRERETk3BRGRXzd0aPQpg3O1q0MvHYGvx6swEcfQQ7NexARERGRTExhVMSXOQ489hgsXsycBz/mlXk3M3AgVKvmdmEiIiIiIuenMCriywYPhk8+YXPnV2g58V7uuAMGDHC7KBERERGRC1MYFfFVn3wCgwZx6O4O3DTjOcqXh0mT1D1XRERERHyDVpWJ+KIff4RHHyWlfgOabR7HocOG+QsgONjtwkRERERELo7GUER8zV9/QZs2UK4cz5SbzpJfcvLxx1C1qtuFiYiIiIhcPIVREV+ydy+0bAnAlPazGDWpEM89B23bulyXiIiIiMgl0jRdEV9x7JhNnZs3s3bkPDr0qESLFraHkYiIiIiIr1EYFfEFjgOdO8MPP7Dvnck0HVyPMmVsDyN/f7eLExERERG5dAqjIr7gtddg4kSOvziIlp8+yIEDMGcOFCzodmEiIiIiIpdHYVQks5s6FV54Adq356mdL/HTTzBtGlSr5nZhIiIiIiKXzyMNjIwxm40xTjq3WZ44v0i2tXQpdOgA9evzUZ0P+N97hn794O673S5MREREROTKeGpktDZw+sq1ksAKYJqHzi+S/fz9N7RuDaVLE/P8V3S9MxdNm8Krr7pdmIiIiIjIlfNIGHUcZ9fp3xtjHgP2A5974vwi2c7OnXYLl5QUdk2YRet7ChMaClOmqGGRiIiIiGQNHl8zaowxwGPAZMdxDnv6/CJZ3saN0Lw5xMWR9PV3tO1fhYQE+OknKFTI7eJERERERDzDI2tGz3AbUB744FwPMMZ0McbEGGNidu3ada6HiWQ/MTFw882wbx/Mn0+vqAYsXgwffgg1arhdnIiIiIiI52REGO0M/OI4zupzPcBxnHGO44Q7jhNetGjRDChBxAd99x00bAh588KSJYzfcBNjxkDv3nDffW4XJyIiIiLiWR4No8aYYkBr4H1Pnlcky5swAVq1gipV4Kef+H5TGF27wq23wtChbhcnIiIiIuJ5nh4Z7QgcBT7z8HlFsibHgddeg0cegUaNYOFCpiwoQatWULWq3WI0h3YDFhEREZEsyGO/5p5oXNQJ+MxxnAOeOq+IL4taFUdkdCzxCYmEBAfSp1kYEbVC7Z3JyfD00zB2LLRvDx9+yJj3c/LUU1C/Pnz9NRQo4G79IiIiIiIZxZNjLg2BykB7D55TxGdFrYpjwPS1JCYlAxCXkMiA6WsBiLiqEDz4IHz1FfTrh/Pqawx+xY9Bg+DOO+GzzyAw0MXiRUREREQymMfCqOM4CwDjqfOJ+LrI6NiTQTRVYlIy//vqFyLmvwFLl8Lbb5PS7Sl69IDRo6FDB/jgA/hm7XlGVEVEREREsgCtRhPJIPEJiWcdC9m/k3emDYQDO2DqVJIi7qZDe5gyBZ59Fl5/Hb5ec54RVQVSEREREckiMmJrFxEBQoLTzrO9aucmpk/qTYlDe2H2bA63vJvWrW0QHTYMIiPBz+/cI6qR0bHeLF9EREREJEMpjIpkkD7NwggM8Afg5i2/Mu2TfjjGj+UTZ7CvRgNuuw2io2HcOOjXD8yJSe7pjaie77iIiIiIiC/SNF2RDBJRKxSSk9n08us8+d17xBUpxZ/jp1H92trccgv88QdMmwbt2qV9XkhwIHHpTfENVkcjEREREck6NDIqklFiYoh4oi29vhlDrqa3UWHDSipWrk29erB5M3z77dlBFNKOqKYKDPCnT7Mw79QtIiIiIuIFCqMinvbff/DUU3DDDRAXB1OnwqxZrN5SkHr1YP9+mD8fbr01/adH1AplaNvqhAYHYoDQ4ECGtq2u5kUiIiIikqVomq6IpzgOfP459OwJ27dDt27wyitQoACLFkGrVpA/PyxcCFdddf5TRdQKVfgUERERkSxNI6MinrBxI7RoAffeCyVLwvLl8M47JAcVYNw4aNbMHl6y5MJBVEREREQkO1AYFbkSR4/a0c9q1WDpUnj7bRtEw8NZtgxuugkefxxuvBF+/BHKlHG7YBERERGRzEFhVORyLVwINWvCiy/aObi//w5PPcWO3f488ogNovHx8MknsGABFC3qdsEiIiIiIpmHwqjIpdq1Czp0gEaN7Mjot9/CtGkkFQ1h1CioUsUG0H79bD594IFTe4iKiIiIiIilBkYiF2vPHvj4Yzst9+BBeO45eP55yJOHBQtsA9316+360FGjIEw7sYiIiIiInJPCqMj5OA788AO8/z58+aUdCW3UCEaPhqpV+fdf6P0ITJsG5ctDVBTceadGQkVERERELkRhVLKlqFVxREbHEp+QSEhwIH2ahaXdSmXnTpgwAT74AP78EwoUgM6d7a1GDY4ehTdfg1dfhZQUePll6NMHAgNde0kiIiIiIj5FYVSynahVcQyYvpbEpGQA4hISGTB9LaSkELFngx0FnTEDkpKgXj144QW46y7IkweAWbOgRw+7m0vbtvDmm1CunIsvSERERETEBymMSrYTGR17MogCFDuwh7vXzqX2mLmwbxsULgzdu0OnTlC1KgBbtsD06fD55/DTT3av0Nmz4bbb3HoVIiIiIiK+TWFUfN4Fp9yeIT4hkbxHD1N3yxruWjePxn8tJ4eTwtIyNQgd8ya0aQO5c/Pnn/DlMLtUNCbGPvfaa2HkSOjaFXLm9M7rExERERHJihRGxaedc8otpA2kSUmwbBnMnUvUZ19yzT+/kcNJYVeeYD64oQ2f1WjKsfKVeK9a45MBdN06+9QbboDhw+2U3EqVvP0KRURERESyJoVR8WlnTrkFSExKJvL734nImQBz59rbwoV2OxZjKFO1Bh/efDcLy9Tgl9CqHNpVhGPrQ8gVXYYaz9lOuPXq2RHQtm2hdGk3XpmIiIiISNamMCo+LT4h8eSfix/YTd0ta6i3eTX1Nq+G5/bZOypXhocegiZN2FujIeviC7Eqeh8L5h1i97eFSN6fBz9/h1sbGwY+B61bQ4kS7rweEREREZHsQmFUfNeePTyw9ReqbIihzpZfqbznXwB25ynAivLXU77lAyzP34Rl28vy22+wfjrs2JH65IIEBRWkRUNo1w7uvNNQqJBbL0REREREJPtRGBXf8d9/sGgRzJ8PCxbAmjW8ChzKEchPwdcxvtRDzE2+jVX7b+D4+jyw3j4tKMg2xb39dvv1mmvsrXRpOyVXRERERES8T2FUMq+DB2HxYo7NXkBS9HwCN6zEz0nhmH9uVuepw/c5X+H7Y4345Xhtju8OwOQ8Tt5ih7ilsUOLW2zgrFoVypRR6BQRERERyWwURiVDXOp2KwApyQ6/f7aKAxO/osiv8ym7Yzk5nONAACu4iQW8wAIasy30RspdlZuwMHigCrwcBmFhULp0Dowp4J0XKCIiIiIiV0RhVDzuordbAfbsgUVTt3Hkw0+oueZjqiavIxk/VvrXZnGJ3myv2hhTtw4VquelbRj0rQSBgV5/SSIiIiIi4mEKo+Jx59xuJTqWO68NZcUKmPN1IomfzaDuXx9zJ7PxJ4U/C9/IsjvGUqbPvYRXLURtTa0VEREREcmyFEbF407fbgUg+XAAiZuKsmZjUe4ctpjW/31MN6ZRgP38l78029v0p0Tfh6lcNcylikVERERExNsURsXjQoID2brnCAdWl+XQ+hBCtu2jG5N42HxMRWcTSbnykty6HTzegQING1LAz8/tkkVERERExMsURrOZy2ksdCkcB+qYmox6Pxd3/vcdPQI6U5+fSMGwO7wOdB9EQNu2BAQFeexnioiIiIiI71EYzUYupbHQ5Vi8GHr3hgPLtrEwsDt1WMjfQSG8V/tRyvV4nGa333DFP0NERERERLIGzY/MRs7XWOhKxMZCmzZwe/39PLLuWdb61+Tm3Gtg7Fgq7PqHx6M/VBAVEREREZE0NDKajZzZWOhCxy9k504YNAjGvefwaM7JxAf1Je+hHZjOneHVV6FIkSuoNvPJ6CnOIiIiIiLZicJoNhISHEhcOsEzJPjSNu48fBhGjIDhw6HK4dXEFutOxe1L4IYbYPTXULu2p0rONDJ6irOIiIiISHajabrZSJ9mYQQG+Kc5FhjgT59mF7elSnIyfPQRVK4Mb764j6nFuhPD9VQ8HgsffAA//XRFQTRqVRx1h82nfP9Z1B02n6hVcZd9Lk/LqCnOIiIiIiLZlcJoNhJRK5ShbasTGhyIAUKDAxnatvpFjezNnQu1akGnx1LolvtDdhaowu2b38U8+ST88Qc89hhcwRYtqSOPcQmJOJwaecwsgdTTU5xFRERERLI7TdPNZiJqhV7ytNK334YePSAi9BcWVuxGoY2/QN26MHo01KzpkbrON/KYGabBemqKs4iIiIiIWBoZlXNyHHj+eejd4xjflu/G9PgbKXToX5g0CX780WNBFDL/yOOVTnEWEREREZG0NDIq6Tp+HJ54AqZ+eIDVoW2pummuHR4dPBjy5/f4z8vsI4+po7PqpisiIiIi4hkKo3KWxES47z74+esdxJa4nZLb18D48dCxY4b9zD7NwtJ0q4XMN/J4OVOcRUREREQkfQqjksa+fXDnnbBj8Z/8UaQ5BfZvh5kzoUWLDP25GnkUEREREcleFEblpLg4aN4c8v3+C7/mb0luHFiwwO4f6gUaeRQRERERyT7UwEgA+P13qFMHKm/8nkU5GpG7UF5YssRrQVRERERERLIXhVFh+XKoVw/u2DeJL4+1IsdVleGnn6BKFbdLExERERGRLEphNJuLjobGjRz6OK8z5sDDmIYN4IcfoEQJt0sTEREREZEsTGtGs7FPPoFHOqQwvuAzPLh7lG2hO2EC5MqV5nFRq+LUWEhERERERDxKI6PZ1MiR8Gj7o0QXut8G0Z49bTpNJ4gOmL6WuIREHCAuIZEB09cStSrOjbJFRERERCSLUBjNZhwH+veHgb3+I6ZoCxrtmgaRkTBiBPidfTlERsem2fsTIDEpmcjoWG+VLCIiIiIiWZCm6WYzI0fCx8O3sa5wC0rtWw+TJkH79ud8fHxC4iUdFxERERERuRgKo9nIkiUwqs9WVuepT7EjuzCzZkHTpud9TkhwIHHpBM+Q4MCMKlNERERERLIBj03TNcaUNMZ8bIzZZYw5Yoz5zRjTwFPnlyuzcyd0uCeRmTkiKOa/B7Nw4QWDKECfZmEEBvinORYY4E+fZmEZVKmIiIiIiGQHHhkZNcYEA0uAxUBLYBdQAdjpifPLlUlOhgcfcBi8vQvVnJWYz2dAePhFPTe1a6666YqIiIiIiCd5appuX2Cb4zgPn3Zsk4fOLVdo8GCoMW8EDzAZhgyBVq0u6fkRtUIVPkVERERExKM8NU03AlhmjJlqjNlpjFltjOlujDEeOr9cpu+/h58HzybS9IV27eD5590uSURERERExGMjoxWAJ4G3gGFATeCdE/eNPvPBxpguQBeAMmXKeKiE7CVqVdwFp87+8w+8cP9G5vnfB1ddAxMmgD4fEBERERGRTMA4jnPlJzHmGBDjOE6d0469BrRxHOfq8z03PDzciYmJueIaspOoVXEMmL42zf6fgQH+DG1b/WQgPXYMmtc9wOiVNxOWfxv+K36BChXcKllERERERLIpY8wKx3HOalrjqWm624Dfzji2AdCwZwaIjI5NE0QBEpOSiYyOPfl9394pdI/pwFX8jv8X0xRERUREREQkU/HUNN0lwJl7fVQBtnjo/HKa+HT2/Tz9+LRpUOCdIbTlK3jzLbj1Vm+WJyIiIiIickGeGhl9C7jJGPO8MaaSMeZu4GlgjIfOL6cJCQ485/HYWPiqQxQvM4iUhzpAjx5erk5EREREROTCPBJGHcf5BdtR9x5gHfAq8CIw1hPnl7T6NAsjMMA/zbHAAH+euuUq+rZcz/tHH+JYzRvwG/c/NSwSEREREZFMyVPTdHEcZxYwy1Pnk3NLbVJ0ejfd3k3DmD08N29ubE1AoSByfjMdcud2uVIREREREZH0eSyMindF1ApNs5XLB/87TsS0lpT3/wf/mQshNPSczxUREREREXGbwmgWsHIl7O8+gGbMJuXdD6BOnQs/SURERERExEWeamAkLtm3DyY1/4Rnkt8g8dFu+HV+zO2SRERERERELkhh1Ic5DgyJWMFruzqxv1YDAv/3ltsliYiIiIiIXBSFUR/20dAd9FoUQVLBYuSP/hwCAtwuSURERERE5KJozaiP2rHdocTALhTz203OuUuhaFG3SxIREREREbloGhn1UV/c/wUtj3/N/t5DMNfVcrscERERERGRS6Iw6oN+id7LXQu7s7XE9RR9tafb5YiIiIiIiFwyTdP1McnJsO2BZ6nFHo59GQ059E8oIiIiIiK+RyOjPubbZ+Zy594J/NG6L3nq1HS7HBERERERkcuiMOpD9vxziOqju/BPYBWunvKS2+WIiIiIiIhcNoVRH7Ky1UuUS9nE8bHvYwJzu12OiIiIiIjIZVMY9REbJv5C419HsrT641ToeIvb5YiIiIiIiFwRhVEfkHI0iYCuj7HTrwTXfDPc7XJERERERESumMKoD1j9wOtUOryWDU+9S4EyBdwuR0RERERE5IopjGZy+5f/zjXTBzOv8N00HHGn2+WIiIiIiIh4hMJoZpaSwq42nTlEXop99g5++tcSEREREZEsQvEmE4t76T0qxi/mm0YjqN6kuNvliIiIiIiIeIzCaCbl/LuV4GH9WBjQhDs+7+B2OSIiIiIiIh6lMJoZOQ7xrbtCcjLbBr5HocLG7YpEREREREQ8SmE0E0qcMJXQVd8wrtQQ7ulfwe1yREREREREPE5hNLPZs4fk7k+znNrUmdoDf3+3CxIREREREfE8hdFM5r/HniHX4X3MiviAG+soiYqIiIiISNaUw+0CBKJWxREZHUvFlYuZOGMiwwOeo9t7NdwuS0REREREJMNoZNRlUaviGDB9Lft27uWVWf/jd8IY0+AulsbFuV2aiIiIiIhIhlEYdVlkdCyJScn0XPQpZQ5v44ngUZha24iMjnW7NBERERERkQyjMOqy+IREyu2N45GVM3mfTvx+e36Mn0N8QqLbpYmIiIiIiGQYhVGXhQQH0mf+JI46uRhaqRu5S+89eVxERERERCSrUhh12dDi+2m5cTHD6cvRBv8BEBjgT59mYS5XJiIiIiIiknEURt3kONw8ehjxlGRi9fbkLHKQ0OBAhratTkStULerExERERERyTDa2sVNn39OwMplvGQ+ZM4XFahSpYLbFYmIiIiIiHiFwqhbjh7leN8BbDDVSXqgA1WquF2QiIiIiIiI9yiMumXsWHJs+Zu+ftG8PdDf7WpERERERES8SmHUDfv2kTJ4CPP8mlL8oaZUrux2QSIiIiIiIt6lMOqGV1+FhAT6+kXyxYtuFyMiIiIiIuJ96qbrbX//jfPOO0z0e4TrOtagYkW3CxIREREREfE+jYx623PPkZTiz0AzmIUvuF2MiIiIiIiIOzQy6k3LlsHUqUSm9KbpI6GUL+92QSIiIiIiIu7QyKi3OA707s3+wOK8cawPq553uyARERERERH3KIx6S1QULF7MAP//cW+nfJQr53ZBIiIiIiIi7lEY9bCoVXFERscSn5BISHAgfZqFEVGtGPTrx7aCV/PRgceIfc7tKkVERERERNylMOpBUaviGDB9LYlJyQDEJSQyYPpayn8Ww7V//skT/t/wSJcclCnjcqEiIiIiIiIuUxj1oMjo2JNBNFWOg/spOy6SDSGN+X7X7WzUqKiIiIiIiIjCqCfFJySedezJnz4n+PB+mh6NpPMThlKlXChMREREREQkk9HWLh4UEhyY9vv9O3k0ZgbTCrfj1xzXMWCAS4WJiIiIiIhkMgqjHtSnWRiBAf4nv++9aBIOhr773uTxxyE01MXiREREREREMhFN0/WgiFo2bUZGx1Lo97W0Xb+ATys9y46tZenf3+XiREREREREMhGFUQ+LqBVKRM0QaDyE5IJF6P73izzxNJQs6XZlIiIiIiIimYdHpukaYwYZY5wzbts9cW6fNGsWLFzIJ5UHcSRXAfr1c7sgERERERGRzMWTI6OxQMPTvk8+x+OytuPHoW9fjpWrQudfuvDUM1CihNtFiYiIiIiIZC6eDKPHHcfJvqOhqT79FDZs4H/1PyfHzgD69nW7IBERERERkczHk910Kxhj4owxm4wxnxljKnjw3L4hKQkGD+bI1TV5ZnFbunWDYsXcLkpERERERCTz8VQYXQZ0BFoAnYESwFJjTOH0HmyM6WKMiTHGxOzatctDJWQCEyfCxo2MLjqY3Hn86NPH7YJEREREREQyJ+M4judPakwQ8DcwzHGcEed7bHh4uBMTE+PxGrzu2DGoUoUj+YuRZ+0y+vQ1DB/udlEiIiIiIiLuMsascBwn/Mzjnpyme5LjOAeB9UDljDh/pvThh7BlC+NKDSZnLkOvXm4XJCIiIiIiknllSBg1xuQGrgK2ZcT5M50jR+DVVzlWuw595zWjQwd10BURERERETkfT+0z+oYxpoExprwx5kbgCyAv8LEnzp/pjRsHcXFMrDSEY0mG3r3dLkhERERERCRz89TWLqWAKUARYBfwM3CT4zhbPHT+zOvwYXjtNY7Xb0jvbxvTrh1Uzj6Tk0VERERERC6LR8Ko4zj3eeI8PmnsWNixgy/afc5/P6J9RUVERERERC5ChqwZzTYOHIDhw0m59TZ6z6hPo0ZQu7bbRYmIiIiIiGR+CqNXYvRo2L2bb28eQlycRkVFREREREQuVobsM3opfHaf0f/+g/LlcW6uwzWbviFnTli1CoxxuzAREREREZHMw6v7jGYLI0fCvn0sajKYDRvsqKiCqIiIiIiIyMVRGL0ce/fCiBEQEcHzX15H2bJwzz1uFyUiIiIiIuI7FEYvx4gRsH8/qyJeZskSePZZyOGpTXJERERERESyAYXRS7V7N4waBffcw8Ava1C4MDz6qNtFiYiIiIiI+BaF0Uv1+utw+DB/tR/EzJnQvTvkzet2USIiIiIiIr5FYfRSbN9ut3N54AFenX41gYE2jIqIiIiIiMil0UrHSzFsGBw7xvYuL/HJrfDEE1CkiNtFiYiIiIiI+B6NjF6suDj43//g4Yd5Y0ZlUlLgmWfcLkpERERERMQ3aWT0Yr32GiQn89/TL/Jefbj3XihXzu2iREREREREfJNGRi/Gli3w/vvw2GOM+bY8Bw9Cnz5uFyUiIiIiIuK7NDJ6MV55BYzhyLPPM6oeNGsGNWu6XZSIiIiIiIjv0sjohWzcCOPHQ5cufDy/NDt3Qr9+bhclIiIiIiLi2zQyeiFDhkBAAMn9nuONRhAeDg0bul2UiIiIiIiIb1MYPZ/YWJg0CXr25KufS/LXX/D552CM24WJiIiIiIj4NoXR8wkIgHvvxenbj+F3QKVK0KaN20WJiIiIiIj4PoXR86lQAT79lAXzISbGbjPq7+92USIiIiIiIr5PDYwuwuuvQ/Hi0KGD25WIiIiIiIhkDQqjF7B6NURHQ48ekDu329WIiIiIiIhkDQqjF/D66xAUBE884XYlIiIiIiIiWYfC6Hls2gRTp8Ljj0PBgm5XIyIiIiIiknUojJ7HkSPQvDn07Ol2JSIiIiIiIlmLuumex9VXw6xZblchIiIiIiKS9WhkVERERERERLxOYVRERERERES8TmFUREREREREvE5hVERERERERLxODYzOI2pVHJHRscQnJBISHEifZmFE1Ap1uywRERERERGfpzB6DlGr4hgwfS2JSckAxCUkMmD6WgAFUhERERERkSukabrnEBkdezKIpkpMSiYyOtalikRERERERLIOhdFziE9IvKTjIiIiIiIicvEURs8hJDjwko6LiIiIiIjIxVMYPYc+zcIIDPBPcywwwJ8+zcJcqkhERERERCTrUAOjc0htUqRuuiIiIiIiIp6nMHoeEbVCFT5FREREREQygKbpioiIiIiIiNcpjIqIiIiIiIjXKYyKiIiIiIiI1ymMioiIiIiIiNcpjIqIiIiIiIjXKYyKiIiIiIiI1ymMioiIiIiIiNcpjIqIiIiIiIjXKYyKiIiIiIiI1ymMioiIiIiIiNcpjIqIiIiIiIjXKYyKiIiIiIiI1ymMioiIiIiIiNcpjIqIiIiIiIjXKYyKiIiIiIiI1xnHcdwtwJhdwBZXi7iwIsBut4uQbE/XoWQGug4ls9C1KJmBrkPJDHzhOizrOE7RMw+6HkZ9gTEmxnGccLfrkOxN16FkBroOJbPQtSiZga5DyQx8+TrUNF0RERERERHxOoVRERERERER8TqF0Yszzu0CRNB1KJmDrkPJLHQtSmag61AyA5+9DrVmVERERERERLxOI6MiIiIiIiLidQqjIiIiIiIi4nUKo+dhjHnSGLPJGHPEGLPCGFPf7Zok+zDGDDLGOGfctrtdl2R9xphbjDFfG2PiTlx3Hc+435y4PuONMYnGmIXGmGtcKleyqIu4Diek8x75s0vlShZljBlgjPnFGLPfGLPLGDPTGFPtjMfoPVEy1EVehz75nqgweg7GmHuBUcBrQC1gKfCdMaaMq4VJdhMLlDztVt3dciSbCALWAT2AxHTu7ws8CzwF1AZ2AnOMMfm8VqFkBxe6DgHmkvY98nbvlCbZSENgLFAHaAwcB+YaYwqd9hi9J0pGa8iFr0PwwfdENTA6B2PMMuBXx3E6n3bsT+ALx3EGuFeZZBfGmEHAXY7jVLvQY0UyijHmINDdcZwJJ743QDww2nGcV08cC8T+8tXbcZz33KpVsq4zr8MTxyYARRzHucOtuiT7McYEAf8BEY7jzNR7orjhzOvwxLEJ+OB7okZG02GMyQlcD8w+467Z2E8kRLylwokpapuMMZ8ZYyq4XZBke+WBEpz2/ug4TiKwCL0/ivfVM8bsNMb8YYx53xhTzO2CJMvLh/39ed+J7/WeKG448zpM5XPviQqj6SsC+AM7zji+A/uGI+INy4COQAugM/baW2qMKexmUZLtpb4H6v1R3PY98DBwK3aK5A3AfGNMLlerkqxuFLAa+OnE93pPFDeceR2Cj74n5nC7gEzuzDnMJp1jIhnCcZzvTv/+xCL0v4EOwAhXihI5Re+P4irHcT477du1xpgVwBagJTDdnaokKzPGjADqAfUcx0k+4269J4pXnOs69NX3RI2Mpm83kMzZn2gV4+xPvkS8wnGcg8B6oLLbtUi2ltrRWe+Pkqk4jhMPbEXvkZIBjDFvAfcDjR3H+fu0u/SeKF5znuvwLL7ynqgwmg7HcY4BK4DbzrjrNmxXXRGvM8bkBq4Ctrldi2Rrm7C/fJ18fzxxbdZH74/iImNMESAUvUeKhxljRgEPYAPA72fcrfdE8YoLXIfpPd4n3hM1TffcRgCTjDHLgSXAE0AI8D9Xq5JswxjzBjAT+Af7CeuLQF7gYzfrkqzvRJe+Sie+9QPKGGNqAnsdx/nHGDMSeN4Y8zvwB/ACcBD41IVyJYs633V44jYI+BL7i1Y5YCi2g+lXXi5VsjBjzBjgISAC2GeMSR0BPeg4zkHHcRy9J0pGu9B1eOL9chA++J6orV3OwxjzJHbvqJLYvc56OY6zyN2qJLswxnwG3IJtqLUL+Bl40XGc31wtTLI8Y0xDYEE6d33sOE7HE1sZDAQeBwpim211cxxnndeKlCzvfNch0BWIwu4DHoz95WsB9j3yX68UKNmCMeZcvyi/7DjOoBOP0XuiZKgLXYcnthOKwgffExVGRURERERExOu0ZlRERERERES8TmFUREREREREvE5hVERERERERLxOYVRERERERES8TmFUREREREREvE5hVERERERERLxOYVRERERERES8TmFUREREREREvE5hVERERERERLzu/9viFi+/S0wFAAAAAElFTkSuQmCC\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           5.061452\n",
       "x1                  0.483995\n",
       "np.sin(x1)          0.420003\n",
       "I((x1 - 5) ** 2)   -0.018380\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    10.986154\n",
       "1    10.871388\n",
       "2    10.661398\n",
       "3    10.393718\n",
       "4    10.117758\n",
       "5     9.882705\n",
       "6     9.725480\n",
       "7     9.661699\n",
       "8     9.681850\n",
       "9     9.753621\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.9.1rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
