{
 "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 statsmodels.api as sm"
   ]
  },
  {
   "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.986\n",
      "Model:                            OLS   Adj. R-squared:                  0.985\n",
      "Method:                 Least Squares   F-statistic:                     1065.\n",
      "Date:                Fri, 24 Apr 2020   Prob (F-statistic):           1.73e-42\n",
      "Time:                        14:17:04   Log-Likelihood:                 4.2343\n",
      "No. Observations:                  50   AIC:                           -0.4686\n",
      "Df Residuals:                      46   BIC:                             7.180\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.0458      0.079     63.868      0.000       4.887       5.205\n",
      "x1             0.5012      0.012     41.134      0.000       0.477       0.526\n",
      "x2             0.5964      0.048     12.452      0.000       0.500       0.693\n",
      "x3            -0.0200      0.001    -18.693      0.000      -0.022      -0.018\n",
      "==============================================================================\n",
      "Omnibus:                        1.977   Durbin-Watson:                   1.969\n",
      "Prob(Omnibus):                  0.372   Jarque-Bera (JB):                1.174\n",
      "Skew:                          -0.332   Prob(JB):                        0.556\n",
      "Kurtosis:                       3.351   Cond. No.                         221.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\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.54585161  5.06544943  5.53948885  5.93546501  6.23260388  6.42527541\n",
      "  6.52391845  6.55332576  6.54857056  6.54924415  6.59295115  6.70913117\n",
      "  6.91422195  7.20895873  7.5782534   7.99367346  8.4181134   8.81189104\n",
      "  9.13926716  9.37431739  9.50519151  9.53606038  9.48643071  9.38794015\n",
      "  9.2791593   9.1992547   9.18155384  9.24806959  9.40588352  9.64598241\n",
      "  9.94473913 10.26779449 10.57570235 10.83041052 11.00151263 11.07124381\n",
      " 11.0373994  10.91369716 10.72752352 10.51543403 10.31714739 10.16901979\n",
      " 10.0980717  10.11754847 10.22474415 10.40144609 10.61692719 10.83299434\n",
      " 11.01026337 11.1146298 ]\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.11028937 10.95038634 10.6583104  10.28736367  9.90771053  9.58919888\n",
      "  9.38425896  9.31506661  9.36811384  9.49751605]\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": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deVyU1RrA8d9hQEDBcMENUqyruCMKmmnuaakpbrnknpmZpVak1nUru+7ZbTUzzUzL3L1qamkuuYPgknsuKW6IGyIoy7l/vGCAoAjDDDM838/Hj/LOO+/7DCPPHM7yHKW1RgghhO1xsHYAQgghskcSuBBC2ChJ4EIIYaMkgQshhI2SBC6EEDbK0ZI3K168uPbx8bHkLYUQwuaFhoZe0Vp7pj9u0QTu4+NDSEiIJW8phBA2Tyl1JqPj0oUihBA2ShK4EELYKEngQghhoyzaB56R+Ph4zp07R1xcnLVDEVng4uKCt7c3Tk5O1g5FiHzP6gn83LlzuLu74+Pjg1LK2uGIB9BaExUVxblz5yhfvry1wxEi37N6Ao+Li5PkbSOUUhQrVozIyEhrh2KTlodFMGXdUc5fj6WMhyvBLX0J8veydljChlk9gQOSvG2IvFfZszwsgpFLDxAbnwhAxPVYRi49ACBJXGSbDGIKYQFT1h29l7xTxMYnMmXdUStFJOyBJHDAZDJRs2ZNqlatip+fHx9//DFJSUkPfM7p06dZsGCBhSIUtu789dhHOi5EVuSJLpRHkRv9iK6uroSHhwNw+fJlunfvzo0bNxg3blymz0lJ4N27d8/RvUX+UMbDlYgMknUZD1crRGM9Mg5gXjbVAk/pR4y4Hovmn37E5WERZrtHiRIlmDlzJp9//jlaa06fPs0zzzxDrVq1qFWrFtu3bwdgxIgRbN26lZo1azJ9+vRMzxP50/KwCOpP3Ej5EaupP3EjTSp54upkSnOOq5OJ4Ja+VorQ8izx85vf2FQL/EH9iOb8FH/iiSdISkri8uXLlChRgl9//RUXFxeOHz9Ot27dCAkJYeLEiUydOpVVq1YBcPv27QzPE/lPRgOWS0Ij6Fjbi9+PRD609WmvrVRL/fzmJzaVwC3Zj5iyV2h8fDyDBw8mPDwck8nEsWPHMjw/q+cJ+5dZovr9SCTbRjR94HPtebaKjAOYn00lcEv1I548eRKTyUSJEiUYN24cJUuWZN++fSQlJeHi4pLhc6ZPn56l84T9y0misudWahkPV65EXqfDn7/jfieGeJMj8Q6OuLm5wneR4OICLVpA0aLWDtVmPDSBK6VmA22Ay1rrasnHOgNjgcpAHa21RfoKglv6pmmdgPn7ESMjIxk4cCCDBw9GKcWNGzfw9vbGwcGBuXPnkpho3Nvd3Z3o6Oh7z8vsPJH/pDQ0CsfdoujtG5wpUhqtHDJsaKTvLsmogQL20UqdWugcXpOGU/b6xfsfXJb8t7c3LFwITz9t0dhsVVZa4N8BnwPfpzp2EOgAfJ0LMWUqpQVi7v7B2NhYatasSXx8PI6OjvTs2ZO33noLgEGDBtGxY0cWLVpEkyZNKFSoEAA1atTA0dERPz8/+vTpk+l5Ip/Rmgmlorm24DOeO7QF58R4bhVw5UjJJynSoC7Mi4CmTcHLK8PuEgXoDC5r07NVzpyBoUOpt3w50T7/4s2OU9hQ2IfH3RwZ0tCH5ysVh7t34e+/4eWXoVEjmDQJhg0DWTj2QCqlr/eBJynlA6xKaYGnOr4JeCerLfCAgACdfmDv8OHDVK5cOYvhirxA3rMM3LoF338PM2bAgQPEF3JjZfWm7PbwIfDG3zSJOUux44cgNhYKFYLPP6f+hceJuHF/Ebf0SdzVycSEDtVtrwtFa5g6FcaONb4ePdpIygUKZP6c69ehXz9YtgyCgmDOHPDwsEi4eZlSKlRrHZD+uE31gQuRJ125As2awf79UKsWfPMNTl270tHNjY6pz0tIgD//hKFDoW9fRlR6hvdbvs5NF7c0l9OAl4erbc9C0Rreegs++cRIxP/9L5QtCxiNbUdHcMhoErOHByxZYpwfHGx8Pxctgtq1LRu/jcj1BK6UGgAMACib/AYKYQsyms4Habvw3q9TnFZDe8Dx47BqFbRunfkFHR3Bzw9++w0mTeL5UaPwP3+EoS+8Q4h31XuneXm4PnS2Sl6T5nv1mAtzji+j4pwvYMgQmD4dlCImBr75BqZNM/Jzhw6ZXEwp40Oubl148UWjP/yHH6BzZ4u+JluQ6wlcaz0TmAlGF0pu308Ic8iofzp40T5QEJ9o/DeOi7jAv7q+TEL0RRxXrTJa4ckuXjR6VE6eNBqe5cr988fLy4TDe++xzbs65YcMYOGCkUxq1JuZdTva5OKe9N+rjmvmUPGP+Zzq1JPy06dz9Zris8/g00/h6lVo2BBKlszChevVg7AwaNcO+vY1PvwqVszdF2NjpAtFiAxkNJ0vPumf9kexmOss+Ok9Hr9+iWF9PuKzZs1ISIC1a2HWLKMxnphozIi7ejXttcuXh4kToXPPF1hd/neOvz6Q9zbNwdWtIOU/GGlz3SWpv1cDdy7mrT/ms6hac8b7vsSVxn9zYYcXOt6RwIax/G+C66NNMCle3JiV4ucH3bvD9u0P7kPPZx66lF4p9SOwA/BVSp1TSr2slGqvlDoH1ANWK6XW5XagQljSg6btFY+5xo8/vkfZ65fo12k0q4pV4ttvjdb1Cy/Ajh3w9ttw5AhERRnjm4cOwS+/wJdfgrs7dOkCzzwDJZwr0Wzvb9C+PcNWfUHQn79b8FWaR8r3qt+eFYzY/B3LqzTinYZvc/S7QM5vKUfBCpco3W8L0Q03c9k1G8vmvb3h228hNBT+/W8zR2/bHtoC11p3y+ShZZkcF8LmZTYnW+kkvlgxicdvXKJfpzFsL1uDuO2V6T8JnL2v4tvjHOPfKEanOv+0ogsVgsqVjT8AAwYYkyv+/W+jm/ellxyZOHkB3jdaQ58+xkBemzYWeqU5V8bDlaq7NzJ64zesqfg0bwSO5sIPT5F025ni7UIpVMmY9x0bT/YXJAUFwcCBMGUKPPus8UdIF0pUVBTNkvsuL168iMlkwtPTE4Ddu3dTQH5dy5cyWjTm5KDovWcVdc8eJPj5IWz39uPaLzWIPuCNW42/KdryIHEOmlH/O4+jU+ZL300m6N/faIVPmAAffwzr1rmw6sfl1L3Z1BisW7/eaKI/gLVqpqS/74uFonl59ceEl67Aq5UnE7GgLg7O8ZR8aTvOpW6meW6OFiRNmwZbtkCvXsaMn+Sf0/zMpqoR5oZixYoRHh5OeHg4AwcOZNiwYfe+LlCgAAkJCdYOUVhBkL8XEzpUx8vDFYUxM+TLOu6M+ON7/qj0FD/7tuTmirpEH/DmsfrHKPrcAZSD0Uee1Y0a3N3hP/8xcpGHBzR+wZ3Vr68x+mLatIHkEscZsVZlv/T3vXEpihfGDQYXF16r+C3nVjxNoZK3qfTq7vuSN+RwQVLBgvDjj3DtmjGomYU1LPYu37fAM9KnTx+KFi1KWFgYtWrVwt3dHTc3N9555x0AqlWrxqpVq/Dx8eGHH37g008/5e7du9StW5cvv/wSk8n0kDsIWxDk7/VPizYx0WgRu7riu3Apnv1KE3ESij23Hze/s/c991FamhUrGmNzbdvCC/08mTXmV/p9W9+YkhgSAqVL3/cca9VMSX1fpZP4ePXHlL16nm61fmDn5mfo0gXmzCnMuiP/yp2yFzVqGN0ob74Jn38Ob7yRs+vZuDyVwIcOfWCjI1tq1jTWEjyqY8eO8dtvv2EymRibspIsncOHD7Nw4UK2bduGk5MTgwYNYv78+fTq1StnQYu8Z9o02LGD+O/m07p/aQ4fhuXLYcLBK0Rcv//0R21penrChg3Qowe8PPZxonqs4p2l9VAdOsCmTeDsnOZ8a1X2S339wdsX0uL4ToY/OZLFod3o08cYa3RwyL2yF8aNB8O6dcZCn0aNjKSeT+WpBJ6XdO7c+aEt6Q0bNhAaGkpgYCBg1FQpUaKEJcITlvTnnzBqFHTowNBd3dizx1jp3aYNJHiZr8BawYLGosO33oJ3P62B49NzGba9MwwaZMxNTFUXxFo7/KTct+mJ3Qz7YwE/lWrL5L8+orjfRWbNKpVmdWWa32DMSSljFLhGDWPZ/e7dmSzrtH95KoFnp6WcW1IXo3J0dEyzR2ZcnFG/QmtN7969mTBhgsXjExYSHw+9e0Phwixp/hVfDlK8844xKQLM39I0mYxVit7e8Na7nShX5d90mD2efcXLM6jI0/fu0aSSJ0tCI3K1MmdGglv6MmPWWj5ZNY0D7pXoe/En3CtfYsY3iVi059DT0+hK6dkT5s0z3qN8KH9+bD0iHx8f9u7dC8DevXs5deoUAM2aNWPx4sVcvnwZgKtXr3LmzBmrxSlywbRpEBrK2fdn0OudEjRoYAw8phbk78W2EU05NbE120Y0NUurMzjYyE+dDo3jD89WVJ0yhrL7dt0bsEzZ4Sf1IKslCl4F/aswP6+dTIJ2ol30Glwr32L2d4l0DLTC4qPu3aFOHXjvPYiJsfz984A81QLPqzp27Mj3339PzZo1CQwMpGLyct4qVaowfvx4WrRoQVJSEk5OTnzxxReUK1fOyhELs7h4EcaPJ6FNEC2+7oibm7Eo0MnJMrd/5x1ITHSg1Ygf2VUggC9WTKRt7+mce6xklnf4MSutoU8f3E6doIVeT5XnfVi27L7uectxcDDqrNSvD5MnwwM2IbdbWmuL/aldu7ZO79ChQ/cdE3lbvnnPBgzQSY6OemjrY9rBQeuNG60ThkfDw7oCR/U1h8L6cPFyutrQhbrc8FXaZ/gqywbyn/9oDfotNU03a6b17duWvX2munTR2tVV67//tnYkuQYI0RnkVOlCESIjBw/CrFkcbPg6n6yuwIcfQpMm1gmlyvPniXxG0ylpCU9eieCLZRNxTEyw7CYPv/yCfv99flTd2FN/GCtWgGte2WNi0iRISoKRI60dicVJAhciI8HBJLkVpvWu0TRvDiNGWDGUlr6UanSK0GceZwAzaXgmjAnrvyK4hYUq8504QXzn7uzXNfg6YBar1yjy1IZT5coZxWfmz4ddu6wdjUVJAhcCY4Vh/YkbKT9iNcP6TYa1a5njPYooXZSZM607Sy1lVWiVVhGsaFqHDxhF5/3raPW/7x/+5JyKiuJW8yBuxjgwqtoylq8viLt77t/2kY0YAaVKGTv+5KMVmpLARb6Xenm4SkpkwKqvOF3Im0GHXueDD4zyr9aWMtPl6oaqeH0zjnn0oMDY97n97YLcu+n161wNbIHjmROMeGIRczaVz7u7m7m7w0cfGaUgFy60djQWIwlc5Hupl4d3PvAblSNP827CJJy87jBkiJWDy8DL/RUu82axmUY4vtKXGys2mf0e+sZNIqq1xO3UAUZVWcqEXU0pVszstzGv3r2NpdfvvmvsPZoPSAIX+V7K8vCCd2N5e+sP7CxYm0V3ulD42X045tGJtp17OBO7YBl/8STO7Z/n6EeLc3zNlG6kKm8uZbdXc0pE7OWLRosYv7cVxYubIejcZjIZ0wrPnjXm7+cDksCBc+fO0a5dOypUqMCTTz7JkCFDuHv3LgCbNm2iTQa1mVetWoW/vz9+fn5UqVKFr7/+Otfj7NOnD4sXGz+o/fv359ChQ5meu2nTJrZv337v6xkzZvD99xboM7VBKbM5Xt21lBIx1xh6+3MK1zlNed94K0f2YM91K8KtVZv5s4A/vv/uzK8tppAQn73+35RupMtnbvHVzOkExIQysMpUfD4OsN487+xo3Bjatze2PLp40drR5Lo82r6wHK01HTp04LXXXmPFihUkJiYyYMAA3n//faZMmZLhc+Lj4xkwYAC7d+/G29ubO3fucPr06WzdPzExMVvVC2fNmvXAxzdt2oSbmxtPJ+9fNXDgwGzFlx8Et/Rlyryt9N+zjMVOQYQWqs4TTbYR3LLqw59sZYGtPLl5ZgO76vbm2V/f5adSx/i4b08iHaMfuKw/fU3vmNhEnP/QfLV1Gs8kbuO1Wh+y4dl/cWT9UdrXsq0t3pg0Cf73Pxg7FmbMsHY0uSrft8A3btyIi4sLffv2BcBkMjF9+nRmz57N7du3M3xOdHQ0CQkJFEvuFHR2dsbX9/4aFGPHjqVnz540bdqUChUq8M033wBGcm3SpAndu3enevXqJCYmEhwcTGBgIDVq1LjXmtdaM3jwYKpUqULr1q3vLdkHaNy4MSEhIQCsXbuWWrVq4efnR7NmzTh9+jQzZsxg+vTp1KxZk61btzJ27FimTp0KQHh4OE899RQ1atSgffv2XLt27d41hw8fTp06dahYsSJbt241x7c4zwvy9+L7i7/inBDP+/ETqdjhGJO6VLWZvSkLl3Sl7smf2NBoCF2vzuKDjz8icWtJzpwhwxrhqQdtk5LgxFZ3np0SwrZNnQhMCuHN+u+z7tmaQO5XN8wVFSoYBcC++cbYy86O5a0WuBXqyf7555/Url07zbHChQtTtmxZTpw4keFzihYtStu2bSlXrhzNmjWjTZs2dOvWDYcM5prt37+fnTt3EhMTg7+/P61btwaM3X4OHjxI+fLlmTlzJo899hh79uzhzp071K9fnxYtWhAWFsbRo0c5cOAAly5dokqVKvTr1y/N9SMjI3nllVfYsmUL5cuX5+rVqxQtWpSBAwemqWG+YcOGe8/p1asXn332GY0aNWL06NGMGzeOT5K/RwkJCezevZs1a9Ywbtw4fvvttyx8k23cqVM8sXges1U/Arv78sO31g4oGxwcGP1cW9bhwH82f0r49j3M3t6P2Z7dGBZ2m4AZRlnxK1dg7JwIrl4sSsINVyrtjOSLm69Qk32sK9WAce36cd7jn4qaFl0sZE6jRsHcucaA5qpV1o4m1+StBG4FWmtUqjKdDzueYtasWRw4cIDffvuNqVOn8uuvv/Ldd9/dd167du1wdXXF1dWVJk2asHv3bjw8PKhTpw7lk+enrV+/nv3799/r375x4wbHjx9ny5YtdOvWDZPJRJkyZWja9P66Fzt37qRhw4b3rlW0aNEHvt4bN25w/fp1GjVqBEDv3r3p3Lnzvcc7dOgAQO3atbPdLWRr9Jgx3E00MaXgGDZNtXY02Xf+eiw/P9WMY2W9eXnbSoac+i/vRE5j3eoWvP74IE47PEmZpLPU4SyPc5ZKHKEjS7jo4smrz73Huor10pSstUR1w1xTvDi8/76RwDdsgORtE+3NQxO4Umo20Aa4rLWulnysKLAQ8AFOAy9qra/lOBor1JOtWrUqS5YsSXPs5s2bnD17lieffJKoqKhMn1u9enWqV69Oz549KV++fIYJPP2HQMrXqcvVaq357LPPaNmyZZpz16xZ88APkZTnPuycR+GcPGJlMpnyx3ZyBw7ADz/wiQ5m0EdelCpl7YCyL6VWd3gZX97oHIznrZfpum8dPcLXs+JWEPxTEZlE5cClgsWYU+kFpj/Tg1vOBfFwdaKQs6PF99jMNW+8AV98YVQFCw21y5rhWXlF3wHPpTs2Atigta4AbEj+2iY1a9aM27dv35uhkZiYyNtvv02fPn0oWLBghs+5desWmzZtuvd1eHh4phUIV6xYQVxcHFFRUWzatOne5g+ptWzZkq+++or4eGPWw7Fjx4iJiaFhw4b89NNPJCYmcuHCBX7//ff7nluvXj02b958r8Tt1atXAXB3dyc6Ovq+8x977DGKFClyr3973rx591rj+VHC8Pe5SWHWVBvOoEHWjiZnglv64ur0z4B4pFtRZjXuwc6Ne4xuhJ9+gm3bWLdmN34jV/L04Dl82HwAt5wL4upkYmzbqmYvi2tVLi7GrtHh4fDDD9aOJlc8NIFrrbcAV9MdbgfMTf73XCDIzHFZjFKKZcuWsWjRIipUqEDFihVxcXHhP6mKPm/YsAFvb+97f8LCwpg8eTK+vr7UrFmTMWPGZNj6BqhTpw6tW7fmqaeeYtSoUZQpU+a+c/r370+VKlWoVasW1apV49VXXyUhIYH27dtToUIFqlevzmuvvZZhovX09GTmzJl06NABPz8/unTpAsALL7zAsmXL7g1ipjZ37lyCg4OpUaMG4eHhjB49OgffQRu2bRuOv/yPiXo4E2cWzbNzvrMqo42YJ3SojnZ0pP4BV8qHuVF/SxyxpcowvlNNi9cSt4ouXSAw0OhOyWRSgi1TOgt1A5RSPsCqVF0o17XWHqkev6a1LpLJcwcAAwDKli1bO/2GB4cPH6Zy5crZjT9PGzt2bJqBRHvxsPcs/RS1PPmruNbEBDTk5t4TfNjrBF/OzUvVmcwnZcZJ+p177DZhZ2TLFmPvzI8+MjZ/sEFKqVCtdUD647neKaS1nqm1DtBaB3h6eub27YSVpZ6ilrJ7TEZT2axNr1pNob1/MLXgGD6YZp/JGx68e31WpC7yVX/ixjz3PmZJw4bGHngTJsClS9aOxqyym8AvKaVKAyT/ffkh5+dLY8eOtbvW98PkNGFYwordpzn50hCOUYFVz9bnj7M2mJSyKCe719vKh3GWTJoEcXEwZoy1IzGr7CbwlUDKLqK9gRU5CSIr3Tgib3jYe5WThGEJy8MiOPDuFJ6MPsnwoh8Q5/u37SalLMhsHndW5nfbwodxllWsaCzumTkTtm2zdjRm89AErpT6EdgB+CqlzimlXgYmAs8qpY4DzyZ/nS0uLi5ERUVJErcBWmuioqJwcXHJ9JycJAxLmLl0N69vnc86nmVXUCmUsuGklAXpZ6ZA1ud35/UP40f20UfG5g/9+tlNtcKHjrtrrbtl8pBZZsZ7e3tz7tw5IiMjzXE5kctcXFzw9vbO9PHglr4ZDprllQUhnefPxz3pFqNrDaGA5z87mdtsUnqIlIHK7Awqp8wrz+i4TXJzg1mzoHlzGD0aMql1ZEusPnHKycnp3ipCYftykjBy281tB+h1ahlfub7CxaYKxT+/9dlsUsqCIH+vbH3/8/qHcbY0awYDBsDHH0PHjvDUU9aOKEesnsCF7cto2uC2Efcv+7cqrYnoPJQSeDCrYxDKZEdJKZfk5Q/jHJkyBX75Bfr2hbAwY8GPjZIELnIk/TzjlJkKQJ76Qd/34Qr8Lmxk+bOfM/at6vaXlHJJdlvveVrhwkalwueeg3HjjOmFNipLC3nMJSAgQKeUQBX2of7EjRn2k3p5uOaZVnjM1TtcKVmFOw6uPH4lHFd3abcI4OWX4bvvYOdOY7VmHma1hTzCvuX1mQpaw2/NJ1Iu4SR3Jn4iyVv8Y9o0o8Zu375w5461o8kWSeAiR3IybdASq/xWjNhB67APCa/6EtWHNTf79YUN8/Aw5oX/+Sd8+KG1o8kWSeAiR7I7z9gSq/x2/3oDv8ndOefkTc9mbW13KbjIPa1aQZ8+xhxxG9x+TX6fFI8uMRFOnoT9+wnav58623YRc+QEoZ5PEl6jPg1e7UKbhwx8PWiVnzkGzS5cgNMvDKIWZ+nYfgrRroWIzqMDrMLKvvrK2KrotdeMrpQhQ8x/D63TbJZhLtICF49m/nxjt5OKFaFTJxg/njLnT1Ohanm6ntrBxB/G0KZpdWOxxPTpcPNmhpfJzb7zu3dhZqP5vHhnAdP9e7GvfMV7j9nzqkuRTS4usGQJdOhgbOs4ebJ5r795szFIevKkea+LJHCRVXFx8Oqr0KMHVK0Ks2dDSAjcugVHjsCvvxqtmN9/N34ILl2Ct94yEv233xqt9lRyc8n9+H4nGXb8NXYWrcmM5veXqs8rA6wiDylQwNjwomtXGD7cPH3ikZFG90zjxhAVlSuVECWBi4f76y94+mljwGf4cNi0yRi5r10bXF3R2uiyiE0sYPxnnTzZ2Kps1y544gno399ogaTaWCInNToyozVMmZDAc/N74OTswOT+75LoYLrvPHtedSlywMnJ2LmnVy9jqf2oUcZ/qkeVlGT8rPj6woIFMHKkMVBar57ZQ5Y+cPFgy5cbrQgHB1i5El54AYCYGKOx/csvsGYNpOx/XKwYPP44eHuDj08dXhizjWaRP2F6b7hRl7lzZ/j4Y4L8jXoq5lpQk5AAQwcnEPB1f55mB4nf/kivKs9w2N6WgovcZTLBnDlGi3z8eLh+3Wi0PKD+Txr79hl96Tt2GJtIfPklVKmSe/FqrS32p3bt2lrYkJUrtQatAwK0PnVKa631sWNat2qltbOz8VChQlq3bav19Olajx+v9cCBWrdpo3XNmlq7uRnnlCyp9TuDYvS5AeN0kqur1u7uWn/2mdYJCWYJ8+ZNrdu1uK2X01Zr0Iljx917bNnec/rpCRu0z/BV+ukJG/SyvefMck9h5xITtX7jDeM/MGhdr57W06Zpffp02vNu3tR6zRqt331X68BArR0ctPb01HruXK2TkswWDhCiM8ipksBFxs6e1bpoUa39/bWOjdVaa718udaFCxuHhw3T+tdftY6Ly/wSsbFaL1midYcOWhcoYPxva+rzlz5V4Vnji7p1td6/P0dhnjundf1q1/VmGuokpbT+/PMcXU+INI4c0fqjj4yfg5RkHhhoJPe6dbU2mYxjTk5aN2ig9ZgxWkdFmT2MzBK4LKUX90tMhGbNSNi9hx6vfclOU3ES9lQlYpMP/6pyl8JtQriqrj1St8f167B0qbFyeetWTW/HBXxqGopbwnXUO++ghr8LRTLcVjVDCQmwYgWMH3yRuZefo5rDIRzmfW8MQgmRG/76y5itsnixsdN9YCA0aWKM+9SrB4Vyb2u+zJbSSwIX9/vgAxgzhhFt32Z+uRZcWelP3BlPCvv9TfGWh0hUOdsg9+BBY83EqrlRjLn1Dn35jjsF3Ijs8CqlJg7DsVzm14qMNEo6z/gyiQrnNvKt40Aed7yAw4plLPesKkWqhGXk0rzuzEgCF1mzdSs0bszaGk3pX28EF+fXIzHGmaLP/om739kMn5LdwlXR0cYg/Y6vwmm5bzIvspBETGwt24PzLw7lpldl7iQ5ER8P8fFw7BhsX3iW7nfn8JrLHErHnUaXKIFasYLlzo/L7uvCbkkCFw8XFQU1a4KLC1VbTeDEombER7lRsusunEvfyPRpCjg1sXWObn39OuxYcIq70ybS4uT3uBJHEopIPBhtd94AABjJSURBVInAi/OUoaDpDo0TN+CANgrz9+9v7Dbu4mITVRGFyK7MErhMIxQGrY2EeOkS7NjB1TdKcPdCEYq3C31g8gbzzKv28IA79Qowsnt7xl5vRLNDIZS+dQWvuCgaut3B//Z5uBOH6vRvYw56ul2c8npVRCFygyRwYVi50pjzPXUqC47W5uIOKFL3FIUqXbx3ipODAgXxif/81pbVedUZ7dqTvmsjpT5KbKHHWBz4z5arGbWi01/Po6AT127H33dfWbQj7FmOVmIqpYYopQ4qpf5USg01V1DCwhIT4d//hooV+bP5EF55BRo0gJmfFsDLwxWFkUSndPZjSie/NMey0sec1cqDWW1FZ3S9W3EJOJnSDirJoh1h77LdAldKVQNeAeoAd4G1SqnVWuvj5gpOWMiPP8LBg8R+t5COXRxxc4OFC6FMGS861bk/OT/qoGBWKw9mdRf0jK4Xn6TxcHWikLOjzEIR+UZOulAqAzu11rcBlFKbgfaAmUt5iVx19y6MGYOuWZO+qztx/Dhs2ABlypjvFlltWWd1F/TMrncjNp7wMS1yGK0QtiMnXSgHgYZKqWJKqYJAK+Dx9CcppQYopUKUUiGRkZE5uJ0wl9Q74Ux58V04eZJtrf7DwkUOfPSRsS7BnLJaeTDI34sJHao/tIsmNysZCmFLcjSNUCn1MvA6cAs4BMRqrYdldr5MI7S+1LvIu8THsXnmAM56lKZNUhglijmxdy84mnloO/3O9fBoc7TTD1g2qeTJktAImfMt8o1c2dRYa/2t1rqW1rohcBWQ/u88LnX/ca+9qyl56yqjiw3n0nknpk83f/KGrLesM5LRgOWS0Ag61vbK1vWEsCc5+nFVSpXQWl9WSpUFOgDmL3grzCql/9j9Tgyv7VzMxrKBrD3cFdd/XaRZs1K5dt8gf69sJdjMBkB/PxIpC3REvpfT9tYSpVQxIB54XWt9zQwxiVyUMtOj/+5lFImLZpTzWHSCiUrtTgG5l8CzSxboCJG5nHahPKO1rqK19tNabzBXUCL3BLf0pXR8DP33LGelTxO2H3+eIoFnGP1SWWuHliEZsBQic7KlWj4T5O/FrPgwCsXHMTb2AxwLxvPZZOc823+cG1uvCWEvZCl9fhMfT9Wl87hcozlh+xvw5ZfwUkMzTvo2s5QPFikTK8T9JIHnN0uWQEQEIxO/pmpVeOUVawf0cNkdABXC3kkCz28++YQbJSow5+LzrJ6dO9MGhRCWIX3g+cnOnbBrF9OThlCrtgPPPWftgIQQOSHtr/zkv//lbsHHmHqlN/Nmpt0RKivlXoUQeYsk8Pzi3Dn0okUs8BiKT3k32rX756H0S91Tyr3Co1ceFEJYjnSh5BdffAFaMzZqMCNHgkOqd/5B5V6FEHmXJHA7tzwsgqYfrOH6J5+zwrkNcd7edOmS9hxZ7SiEbZIEbsdSukae2r4Gj7hbTIt9B1XjCKsOpN0JR1Y7CmGbJIHbsSnrjhJ7N4G+ISvZ61SDHYVqU6Dy2fu6RmS1oxC2SQYx7dj567HU+/sAFaLO0psPKfzMKZRj0n1dI7LaUQjbJAncjpXxcOWlFb8Q5VCExQXaUcxv573j6clqRyFsj3Sh2LF/BxSl5dHtzE3qQ4E6F3EokChdI0LYEWmB27Hn9/wCOpFvHPvj7n8GL+kaEcKuSAK3V0lJJHw1k62qCc0HVuGzT6pYOyIhhJlJF4q9Wr8ex7OnmaEHMnSotYMRQuQGaYHbkdT1TL5bMYHKqgS6bRBPPmntyIQQuUFa4HYi9e7tJW9eocGRHXyr+1Gz7Q1rhyaEyCWSwO1E6nomL+5fj0Iz1/NFVl06YOXIhBC5JUcJXCk1TCn1p1LqoFLqR6WUi7kCE48mZXGOKSmRrqG/sY6WRNXTXLgh9UyEsFfZTuBKKS/gTSBAa10NMAFdzRWYeDQpi3Oa/rWHMnGX+ca1DwV9L0o9EyHsWE67UBwBV6WUI1AQOJ/zkER2pNQz6bbrVyIow9Y6lSno7CCLdoSwY9lO4FrrCGAq8DdwAbihtV6f/jyl1AClVIhSKiQyMjL7kYoHCvL34r91CtM4YjffOvTlicbXmNChuizaEcKO5aQLpQjQDigPlAEKKaV6pD9Paz1Tax2gtQ7w9PTMfqTioZ5as4IkHNAvv8ausY0keQth53LShdIcOKW1jtRaxwNLgafNE5Z4ZHFxmL77lpW0o8/7kriFyA9yksD/Bp5SShVUSimgGXDYPGGJR3V77iIKxUVxuMkgypWzdjRCCEvI9kpMrfUupdRiYC+QAIQBM80VmHg01//zJX/jy/NTmlo7FCGEheRoForWeozWupLWuprWuqfW+o65AhNZF79rL2X+3smGCq9Rq7aydjhCCAuRlZh24PTwr7iNKxU/6m3tUIQQFiTFrGxUSuGq6AuR7Nw8nxWFu9G1k4e1wxJCWJC0wG1Q6sJVbbfsoiCxzPFrwYrwiIc/WQhhNySB26B7hau0pufB/7HToQ5H6xS+b7d5IYR9kwRug1IKVwXuO4Fvwgm+q9ghw93mhRD2TRK4DUopUPXS9l+5QjE2NKuS5rgQIn+QBG6Dglv6UvZ8NG2i1zGvVEfi3Rxkt3kh8iGZhWKDgvy9SPpjM0k4sKBFY9ltXoh8ShK4DTq5+wrPnZpLWOWX2DO3m7XDEUJYiXSh2KD9Az6jILE8MeNda4cihLAiSeA25sS+GJ7Z9zkHnmyHZ8PK1g5HCGFFksBtzM7+syjGVbz+O9zaoQghrEwSuA05ejCehiHTOOn9DEVb17N2OEIIK5MEbkN+H/AjZTlL0YnS+hZCyCyUPCelSNX567GUSTU98MihJBrsmMyF4tUo3b2VtcMUQuQBksDzkJQiVbHxiQBEXI9l5NIDABwfGkYwfxI97ntQUvNbCCFdKHnKvSJVqcTGJxL8+VnqbZnI2QJlaH2tJMvDpOqgEEISeJ6SUTEqnaTwW/IXDdjGjAbt+Ts6npFLD0gSF0JIAs9LMipGdXd3aSbfGMVh9yf4qXYLwGiVS+lYIYQk8DwkuKUvrk6me18n3HSh39bVPMEpPmzVj0SHfx6T0rFCiGwncKWUr1IqPNWfm0qpoeYMLr8J8vdiQofqeHm4ooAiG4rxXtIEVpd/hu0+NdOcK6VjhRDZnoWitT4K1ARQSpmACGCZmeLKt4L8vQjy92L1arg+6SWcTIlMb9UvzTlSOlYIAebrQmkG/KW1PmOm6+Vrt2/DnP7beIkFmN4NZvDLz95rlXt5uDKhQ3UpHSuEQGmtc34RpWYDe7XWn2fw2ABgAEDZsmVrnzkjOf5h3h+RSMdJgVTxjMTl1BEoVMjaIQkhrEgpFaq1Dkh/PMcLeZRSBYC2wMiMHtdazwRmAgQEBOT808IGZLaaMivWr4crk2dTizD49EdJ3kKITJljJebzGK3vS2a4ls170GrKhyXxI0cguNNJNjuMIKHuMzh26ZLr8QohbJc5+sC7AT+a4Tp2IbPVlA+btx0VBV1b3WTh7Rdwd9M4zpklS+aFEA+UowSulCoIPAssNU84ti+z+dkPmrd99y682CGBiae7UpFjmJYuhooVcytEIYSdyFEXitb6NlDMTLHYhTIerkRkkKwfc3Wi/sSN9/WLaw2DB0ObLcE8xy/w5dfQtKkVIhdC2BqpRpiJjAYigYcOTga39E3TBw7g5KCIuZvA9dh4IG2/+F8bveCbmQzjExg6FAYMsNArFELYOrNMI8yqgIAAHRISYrH7ZVf6gUgwkjAK4hP/+X65OpkynJOdPvnfvpvAtdvxac7RCQ7EbKlKjT0nWEdL9lSszeX5PxMUUDZ3X5wQwubk2jRCe5TRQGR80v0fdCmDk+kTeMpqyhTlR6xO83jCDVeuLPfj5Ys/MdkhmJNFvBjQ6m0SVxwCk0kW6QghskSKWWXgUQpFZeXc1HVLYv/ypMBsb1Zf6swXDCakXGV6dP2QW84FpcqgEOKRSALPwKMUisrKucEtfXGMLcjVDZVou3gf++JrU9dxJyNaDqZ353FEuhW9d65UGRRCZJV0oWQgs4HIjPrAH1RUSmvYuuEuIR+dptvmP2il11CPnYT4+DO63VAOOd8/gUeqDAohskoSeAZS+qBTBiKrON5hZNkE3P86xqlNu+BWNI6urlQt70n5uUXhJ2dwcACtSYjXXLqkOX82ibuhB6gd8wcNuQ3Ancp+MPgLAgYOZMC+C/d9SEiVQSHEo5AEnomg4kkEHZ4La9fC5cv3jvsVKQJFi8LNO+i/75AYewfu3EEnJpGkFUlaURRFERTnXZ7gTNN+PNm/Cc4tGuFc7J8Wd/oPiUetmSKEEJLA04uJgcmTYcoUow/kxRehZk2oVo3YJ6qyKrQ0m7co9u6F8HCITe6yVk4JuHtH06yRIz3bulO3LvyrzINvlX62ihBCPApJ4Bjztqf+cpiAHWt5b8tcSty8Al26wKRJ6LLl2LYNvv8efv4ZbtwAd3fw94fm7W8RcuskFL+GU7EYlIPmsJMJ5VOdMmUkMQshcle+T+DLwyIYtSiM6T9/QPO/9rC/1L8YFjScjm92I2KlF598AidPGlVdO3aEXr2gcWMwmaD+xN0USDdrJLO54UIIYW75PoFPWXuEsf+bTvO/9vBhk5eZHdiO+FuurO7qyo0T0KABjB5tJG83t7TPzU7hKiGEMJd8n8C7/28mHQ9u5OMGL/FtnfbcPlaSqF9qoBMdmDkT+vfPvKprZoWrZCqgEMIS8vdCni+/5PUdP7PAryWfBLxE1NpqRC4LwNHjNjXfCOGVVx5ckju4pS+uTqY0x2QqoBDCUvJvC3zZMhg8mAsNn2V83SFc+uEp7l7woHDdE5Ru+hejO1d76CVkKqAQwpryZwLftg26dYO6dSm5ejkeQYq7F13wDNpLhbo3CG5ZLctJWKYCCiGsJd8l8JU7ThDYpiNxBYvyevN3cX0zkR0b3Jk2Dd56q7a1wxNCiCzLVwl8eVgEZ0Z+SNvrl+jW9T/sDn2CqF/cadExhmHDZPd3IYRtyVeDmN/9/AevbF/ImopP87tqTNS66rj4RHLTf7fsHyyEsDn5qgXeZ+UMTElJfOA/mMhltXEqEoNnu71cjE6wdmhCCPHIcrorvYdSarFS6ohS6rBSqp65AjO77dsJOrSJmYEdCNvUCpTGs9MeHFwSZN62EMIm5bQF/l9grda6k1KqAFDQDDGZX1ISvPkmsSVKMd39be5eeozibffi5BEr87aFEDYr2wlcKVUYaAj0AdBa3wXumicsM5s7F0JDSfjqByJH1KbwE9coVOkCXjJvWwhhw3LSAn8CiATmKKX8gFBgiNY6JvVJSqkBwACAsmWtsOP6zZswciTUq8e74d2JvaUI/6MI1aq1tnwsQghhRjnpA3cEagFfaa39gRhgRPqTtNYztdYBWusAT0/PHNwumyZOhEuXODrov3w9U/H661Dt4YsshRAiz8tJAj8HnNNa70r+ejFGQs87oqPhiy/QL77IyzMCKV4cxo2zdlBCCGEe2U7gWuuLwFmlVMoIYDPgkFmiMpe5c+HmTdZVHsa2bUZj3MPD2kEJIYR5KK31w8/K7MlK1QRmAQWAk0BfrfW1zM4PCAjQISEh2b7fI0lKgkqVSCxcBO+IXZQtCzt2GHsPCyGELVFKhWqtA9Ifz1E601qHJ/dv19BaBz0oeVvc2rVw/DhDYl7i4kWIDdjDyn0R1o5KCCHMxm5XYl76cDIUKso3f/WnYOUIbrpfZuTSKACZNiiEsAt2k8CXh0Xcq8td7+5lFuzczHifIdw9XZBiTxndNrJfpRDCnthFAl8eFsHIpQeIjU8EoNXvi4gzOfH55bdx8YmkQInoe+fKfpVCCHthF0N6U9YdvZe8H4uNpsOfG1lcqjWXbj9O4bp/pTlX6p4IIeyFXSTw1K3qrvvXUTD+Dh9HD8epxA1cykXde0zqnggh7IldJPCUVrUpKZFeoav5o3htwm4+RelnzuBdxBUFeHm4MqFDden/FkLYDbvoAw9u6cvIpQdofOAPvKIjGVpgGo6PxfLx8GJ0DKxh7fCEECJX2EUCT2lVe/34Pn8XKs3SqO70DY6mY6C0toUQ9ssuulAAgjw1gSfD2PH4KzxWxMQno2XNvBDCvtlNAmfBAtCaUUd6MGgQuLlZOyAhhMhd9pPAf/iBkyWf4m/nCrzxhrWDEUKI3GcfCXzfPjhwgM+v96RLFyhZ0toBCSFE7rOLQUzmzSPR5MTcO11Y3t/awQghhGXYfgJPTIQFC9j2WCuKFy9GgwbWDkgIISzD9rtQNmyACxf49GoP+vcHpawdkBBCWIbtJ/B584h1foy1pjb06mXtYIQQwnJsO4HfuoVeupRF6kVatHWRwUshRL5i2wl8+XLU7dt8E9eT/jJ4KYTIZ2x7EHPePC66+nC6SH1atrR2MEIIYVm22wK/cAH922/Miu1Bn34OmEzWDkgIISwrRy1wpdRpIBpIBBIy2jU51yxYgEpKYh49WNvPYncVQog8wxxdKE201lfMcJ0sWx4Wge/0b4hx8OdSeU/2XY+gPFJ5UAiRv9hcF8rysAg+m/0blSOO8nNSV5yqnGbk0gMsD4uwdmhCCGFROU3gGlivlApVSg3I6ASl1AClVIhSKiQyMjKHtzP2v2x6cAsASwu0pWDFS/d2mxdCiPwkpwm8vta6FvA88LpSqmH6E7TWM7XWAVrrAE9Pzxzeztj/stXhbYRSi0tVXVGOSfeOCyFEfpKjBK61Pp/892VgGVDHHEE9iL++if+loyyiM4Wq/NNtIrvNCyHym2wncKVUIaWUe8q/gRbAQXMFlplxCUcAWOrWBmeva4DsNi+EyJ9yMgulJLBMGdWjHIEFWuu1ZonqAXw3r2cv/sQGuuGojJZ3cEtf2W1eCJHvZDuBa61PAn5mjOXhzp6lQOhOFvMR67/woXJlH4veXggh8hLbmka4eDEAh6t2pnJlK8cihBBWZlO1UG7PW8wx/Gj4cgVrhyKEEFZnOy3wc+coGLadxaozXbtaOxghhLA+m0ngevESAM7X60Tp0lYORggh8gCb6UKJnrOI01Sn8asyXVAIIcBWWuARERTev43ljp1p397awQghRN6Q5xP48rAIpr06EYBfnmjChhNStEoIISCPJ/DlYRGMXHqApJNx7CaQv2o4S+VBIYRIlqcT+JR1R4mNT2ScwxiedtmMq0+kVB4UQohkeXoQM6XCYNGWB4iPckOZdJrjQgiRn+XpFnhKhUEHpyScS92877gQQuRneTqBB7f0xdUp7W7FUnlQCCEMeboLJaXC4JR1Rzl/PVYqDwohRCp5OoGDkcQlYQshxP3ydBeKEEKIzEkCF0IIGyUJXAghbJQkcCGEsFGSwIUQwkYprbXlbqZUJHAmm08vDlwxYzi2QF5z/iCvOX/IyWsup7X2TH/Qogk8J5RSIVrrAGvHYUnymvMHec35Q268ZulCEUIIGyUJXAghbJQtJfCZ1g7ACuQ15w/ymvMHs79mm+kDF0IIkZYttcCFEEKkIglcCCFslE0kcKXUc0qpo0qpE0qpEdaOxxKUUqeVUgeUUuFKqRBrx5MblFKzlVKXlVIHUx0rqpT6VSl1PPnvItaM0dwyec1jlVIRye91uFKqlTVjNCel1ONKqd+VUoeVUn8qpYYkH7fb9/kBr9ns73Oe7wNXSpmAY8CzwDlgD9BNa33IqoHlMqXUaSBAa223ix2UUg2BW8D3WutqyccmA1e11hOTP6yLaK2HWzNOc8rkNY8Fbmmtp1ozttyglCoNlNZa71VKuQOhQBDQBzt9nx/wml/EzO+zLbTA6wAntNYntdZ3gZ+AdlaOSZiB1noLcDXd4XbA3OR/z8X4j283MnnNdktrfUFrvTf539HAYcALO36fH/Cazc4WErgXcDbV1+fIpW9GHqOB9UqpUKXUAGsHY0EltdYXwPhBAEpYOR5LGayU2p/cxWI33QmpKaV8AH9gF/nkfU73msHM77MtJHCVwbG83e9jHvW11rWA54HXk3/1FvbpK+BJoCZwAZhm3XDMTynlBiwBhmqtbz7sfHuQwWs2+/tsCwn8HPB4qq+9gfNWisVitNbnk/++DCzD6ErKDy4l9yGm9CVetnI8uU5rfUlrnai1TgK+wc7ea6WUE0Yim6+1Xpp82K7f54xec268z7aQwPcAFZRS5ZVSBYCuwEorx5SrlFKFkgc/UEoVAloABx/8LLuxEuid/O/ewAorxmIRKYksWXvs6L1WSingW+Cw1vrjVA/Z7fuc2WvOjfc5z89CAUiebvMJYAJma60/snJIuUop9QRGqxuMjacX2ONrVkr9CDTGKLN5CRgDLAd+BsoCfwOdtdZ2M+iXyWtujPFrtQZOA6+m9A/bOqVUA2ArcABISj78HkafsF2+zw94zd0w8/tsEwlcCCHE/WyhC0UIIUQGJIELIYSNkgQuhBA2ShK4EELYKEngQghhoySBCyGEjZIELoQQNur/bKXMPsIwMnIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 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.045799\n",
       "x1                  0.501192\n",
       "np.sin(x1)          0.596429\n",
       "I((x1 - 5) ** 2)   -0.019998\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.110289\n",
       "1    10.950386\n",
       "2    10.658310\n",
       "3    10.287364\n",
       "4     9.907711\n",
       "5     9.589199\n",
       "6     9.384259\n",
       "7     9.315067\n",
       "8     9.368114\n",
       "9     9.497516\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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
