{
 "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.987\n",
      "Model:                            OLS   Adj. R-squared:                  0.987\n",
      "Method:                 Least Squares   F-statistic:                     1208.\n",
      "Date:                Mon, 24 Feb 2020   Prob (F-statistic):           9.91e-44\n",
      "Time:                        22:49:06   Log-Likelihood:                 7.1613\n",
      "No. Observations:                  50   AIC:                            -6.323\n",
      "Df Residuals:                      46   BIC:                             1.325\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.9615      0.075     66.588      0.000       4.812       5.112\n",
      "x1             0.4948      0.011     43.061      0.000       0.472       0.518\n",
      "x2             0.5677      0.045     12.567      0.000       0.477       0.659\n",
      "x3            -0.0190      0.001    -18.864      0.000      -0.021      -0.017\n",
      "==============================================================================\n",
      "Omnibus:                        0.123   Durbin-Watson:                   1.994\n",
      "Prob(Omnibus):                  0.940   Jarque-Bera (JB):                0.021\n",
      "Skew:                          -0.041   Prob(JB):                        0.990\n",
      "Kurtosis:                       2.943   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.48569704  4.98751642  5.44597256  5.83012653  6.12020511  6.31084944\n",
      "  6.4119955   6.44724165  6.44997147  6.45786894  6.50672707  6.62456723\n",
      "  6.82703528  7.11483101  7.4735932   7.87625909  8.28751072  8.66957758\n",
      "  8.98844222  9.21942917  9.35125908  9.38790189  9.34792457  9.26144063\n",
      "  9.16516252  9.09636994  9.08678483  9.15735959  9.31483463  9.55063095\n",
      "  9.84225965 10.15701679 10.45735656 10.70705993 10.87718535 10.95082378\n",
      " 10.92587676 10.81540122 10.64546423 10.45086027 10.2693946  10.13567196\n",
      " 10.0754111  10.10121927 10.21052089 10.38598069 10.59835186 10.81128124\n",
      " 10.98728199 11.09389262]\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.0986701  10.95537979 10.68628459 10.3421188   9.98966661  9.69541101\n",
      "  9.50925641  9.45231034  9.51171553  9.6437979 ]\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3zO5f/A8de1g52MOW2YM+V8GAulJOecQiqHSCmnkEK/6cCiItRUpK9yjJDTVOQQkeQ0mbPJYdjIeQwb2+7r98dnY2ab2e7tPr2fj4fH7LPPfX/et7veu+7357rel9JaI4QQwvY4WToAIYQQ2SMJXAghbJQkcCGEsFGSwIUQwkZJAhdCCBvlkpcXK1q0qC5XrlxeXlIIIWzerl27Lmqti6U9nqcJvFy5coSFheXlJYUQwuYppU6md1xKKEIIYaMkgQshhI2SBC6EEDYqT2vg6UlISCAqKor4+HhLhyKywN3dnVKlSuHq6mrpUIRweBZP4FFRUXh7e1OuXDmUUpYOR2RCa82lS5eIioqifPnylg5HCIdn8RJKfHw8RYoUkeRtA5RSFClSRD4tCWElLD4CByR52xB5r7IvdHc0E9dEcCYmjpI+HoxoVZmOAf6WDkvYMKtI4ELYu9Dd0Yxcto+4hCQAomPiGLlsH4BjJXGTCVatgv/+g9u3ISHB+Hr7NiQlwfPPQ/Xqlo7SZjwwgSulZgLtgPNa6xrJx14AgoGqQH2ttU2vznF2dqZmzZokJCTg4uLCK6+8wtChQ3FyyrjCFBkZyd9//0337t3zMFJhqyauibiTvFPEJSQxcU2E4yTwAwegb1/4+++Mz/nkE+PP22+Ds3PexWajslIDnw20TnNsP9AZ+NPcAT1I6O5oGo3fQPmglTQav4HQ3dE5fk4PDw/Cw8M5cOAA69atY9WqVXz00UeZPiYyMpIff/wxx9cWjuFMTNxDHbcr8fHw4YcQEAARETBrFpw6ZYzCL1+G69fh1i04dw7atoURI6BpU4iMtHTkVu+BCVxr/SdwOc2xQ1rriFyLKgMpH0OjY+LQ3P0Yao4knsLX15fp06czZcoUtNZERkby1FNPUbduXerWrcvfyaOHoKAgNm/eTJ06dQgJCcnwPCEASvp4PNRxu7FxI9SqBR9/DF27wqFD0Ls3lC4Nfn5QqBB4eUG+fODrC0uXwuzZsHu38bjZs0F2DctQrs9CUUr1VUqFKaXCLly4kKPnyuxjqDlVqFABk8nE+fPn8fX1Zd26dfzzzz8sWrSIIUOGADB+/HieeuopwsPDefvttzM8TziwuDg4exa0ZkSryni43lsS8HB1ZkSryhYKLg8sXAjNmhm17bVrYe5cTEWKsWQJnD+fwWOUgldegb17jRH7q68adfGLF/M0dFuR6zcxtdbTgekAgYGBOfpVmpcfQ1P2Ck1ISGDQoEGEh4fj7OzMkSNH0j0/q+cJB3D6NEdGjaf4oh8oEBfLzXzuNClfkXVlKvB7YkF+K1GDqJqP2dUslLQzbD53i6ThiL7w5JOwahUJ+bxYMBfeG51AdKQrhZ4+RNXWZzP+NyhXDjZsgJAQeP99o7Ty55/g5pbXL82q2dQslJI+HkSnk6zN/TH0+PHjODs74+vry0cffYSfnx979uzBZDLh7u6e7mNCQkKydJ6wU1rDtm0weTKmJUupaNKEurZnk3qGShyh8skj1Dh9gF5xkfTWJmNk2e9zS0dtFmln2FT85y/qLhvL5ao18VzyK7PmejFhglHSdvONo2iHfXhWPkt0DJnPxHF2huHDoXx56NLF+PvXX+fhK7N+NpXAR7SqfM9/KGD+j6EXLlygf//+DBo0CKUUV69epVSpUjg5OTFnzhySkoxre3t7Exsbe+dxGZ0n7Et6c7lVYiLub/Sh1Z71XHUqwHd6KFMYxH9+3rj5XyExtgtJVz3QsZ44a837fEzQ7M/Qv/6Gy/Rp0LGjpV9WjqQubTY4tY/pyz/h3yJl6NdgEqfqeXP6NDRoAO6N93CzeBSplxJkaSbO888bs1JCQqBRI6OWLoAs1MCVUguArUBlpVSUUqqPUqqTUioKeBxYqZRak9uBgvFbelznmvj7eKAAfx8PxnWumeOPoXFxcdSpU4fq1avTvHlzWrZsyejRowEYOHAgc+bMoWHDhhw5cgQvLy8AatWqhYuLC7Vr1yYkJCTD84T9SO8metCif3Dq2ZtWe9Yzhg8pne8EYwIHktDnJMV7bKNQkwiKtQ+n+MtbKTlgPZdueuA+6ROaeO7kwAU/6NSJ2527wgPuD+XG7CtzSSlhBkQfZsbSMZwu6EeXev/j7zlNcHKCdetg61aIL3Fv8k77+Ex99pmRvF9/3bgRKgBQOg/v8AYGBuq0GzocOnSIqlWr5lkMIucc9T1rNH7DPSU8J1MSn6+YTKcjfzCCz5jZog35a51GuZjSfby/jwdbgpoCxk280e8lUHTGZ3zIGBKKlMArbJNR+00jbYkCjE+e5hi8mEOj8RtwOXGMn+e8zRWPAjxbZQGHtz1N/tJXObqzIH5+d89LrwSa+t8lU9HRxo3NokVhxw7In9/Mr8R6KaV2aa0D0x63qRKKEJZ0Jk3ynhj6NZ3+/YMgp0+Z27kJ3hXT3TQFuL/U5+sL0753ZfebHzCgd2sm7W1BQp0mFAzfxIorLveUaW7eTrSqRUBpy0itynrR9YtPMSkn2vot5vC2RnhXO8t335vw8yt453E5LoH6+8OCBdCypbEgaP580h3SOxBJ4EJkIG2i8vF05crNBJxMSXy27BueP/Y777uMYU7Xxnj431sC8fFwxcvN5YF9T04SzcEXY+lwYx4/H3uZM1Ub82Wfj4nOXxgg3RFrCkssArqvJcCVmzz2w0dUuniKzsV/4mBEI0o1PcGXE/LRud69rzfl9eeoH0yzZjBmDHzwgTHDZeBAs702WyQJXIh0pNe7xNVJ4eqsGLZiAS8cW8OofKOY17sx+YteJfUA2cPVmeAO1R+YmFJfQz8PnX6bQei+15j5vw/o+upYzhUukunjLbEIKO1ajDd2LOfZg38yuvAofjn3PHPmQK9eGbca7hjgn/NPDSNHGkX1oUMhMBDq18/Z89kwi7eTFcIapbdoLMGkafTvEfofXsg0tzdY9U5LvuxfiYldamfrxnrqaygFkW3c6NJgGoUTY1gwczTFzsdk+FhLLQJKPep/IjKcoE2zWebVjjGXg5k9G3r1yoMgnJxg7lwoWdJY9HP7dh5c1DrJCFyIdKRXnvCMjyN4WQhHVSWa/h3CgLp3ZxplZ1SZ3jWONinAi+5f8dOmIcyb8xFd+owjtrB7lksyuS1lLYb/1fNM+XkCEa6VeOXGj1ToGEHPnlXuOTdX2+cWLgxTp0K7djB5Mrz7rnme18bICFyIdKRXnhi+cCllkk5z7IPZVK6b82miGZVAjjQsSq9m4ylnOsmU2SF4xTsR3KE6W4KacmJ8W7YENbXY7JMRrSrjQyLTln+K8y1Nx9u/4NE8ms9Hed9zXl70LaJtW3juOaMmHhVlvue1IQ6fwC9dukSdOnWoU6cOxYsXx9/f/873tx34o5mjS9u7pMHfx3jt3EKWVBpCqzGNcuUaYJRGejQsw8XmlXmr4fs8nbCFT2f9QJNyJc1yzZzqWKckKw7Mp9a5o/Q0zSO+mQvTJ3jf9wslr/oWMXmy0Wtl2DDzPq+NcPgSSpEiRQgPDwcgODiY/PnzM3z48HvO0Vqjtc60P7iwL6lnTMT8G8vnf03gX9cqtN8xPleukW6ZIagpEb1deW7Oe8yoWYkXDo6hQAGzXT5b9NRvKPvrYoIZTeXhHVgxIf2ZfHnWt6hcOXjvPRg1Ct54A5o3N+/zWznJSBk4evQoNWrUoH///tStW5fTp0/j4+Nz5+cLFy7k9ddfB+DcuXN07tyZwMBA6tevz7Zt2ywVtjCjjgH+bBrelEm/Lqe4/o98P87Fo5B5e9x0DPDPtDRSeVYQkS1ep0/0WKYGzuL6dbNe/h4PXO3511+Y3hrKL7TjQv9RTMggeUMet88dMQIqVoRBgxzuhqZVjcCHDoXkwbDZ1KljfMrKjoMHDzJr1iy+/fZbEhMTMzxvyJAhvPvuuzRs2JDIyEjatWvH/v37sxmxsCaLev5Kj/Nz2N/xfWp0eSzvA1CKciu/4b/HTjF8T1+CGpXmoy3NyZ/fvDcJH7jl25kz3GjzAtGm8qx86Qe+meqU6RqavOhbdIe7u9Hkqk0bh7uhaVUJ3NpUrFiRxx578P+0v//+OxERd2t7V65cIS4uDg8PO2/Wb+cOhd+i4cK3iCpYnRqLRlkuEFdXiv+5mCs1nyJ4byeGPPEHzaaU5ON12d9jM23yz3S1Z7WinGvcBa/YWL5p8TtT5vnwoGqiWRbtPIxnn717Q7NbN2PDCAdgVQk8uyPl3JK6IZWTkxOp+8bEx8ff+bvWmh07dpAvX748jU/kHpMJ1neawiCOEzN9jbFjjCUVKEChrb9xI+BJJuxrTasXfuN6VyecPe4m3awur09vtJ2RM1ducrTtECod28ondRYz4dfquGQxa5hl0c7DmDwZqlY1bmj+9FPeXdeCpAaeRU5OThQqVIh///0Xk8nE8uXL7/ysefPmTJ069c734eauA4k89+NXF3k5ciynaz6Lz4stLR2OoWRJvP5eR/5Crqw43wn3eSVIunHvL5as3CRMb4ZIRt4JXUWl9dOZXzqId/7uYvHfY5kqV87Y/GHxYvj9d0tHkyckgT+Ezz77jNatW9OsWTNKlSp15/jUqVPZsmULtWrVolq1anz33XcWjFJkR+obeI99sJnr7waTn+uUWjDJ0qHdq2JF3DetxTvfNVZd7oyeX4HE2Lu71GTlJmFWkrw2QZ/5Gxl8ZBpryrxBx4OfYhMVweHDoVIlh7mhKe1kxUOzt/csbUmh8FIvdhxtzsFWL1Nr9WzLBpeBTbN/5rHXX+KI6VGaea7BpX0khSvFZGkZf0ZtXVNWe0ZduMXL83fyybnRbCjXlVF9XiP6+m2LrgB9KCtXGis0P/8c3nnH0tGYRUbtZGUELhxe6pJCXGQRgo9+wU1nd0YGtrdwZBl7uncHdn81g+rqIOvjW+Gx0I9q0Y1ITOCBGz9ktIAouEN1lr3alAHLTvLJudGEPdKBId1fJur67dxbTZkb2rY1ZqQEB8N//1k6mlwlCVw4vJSSginBidorL9KBX5j6xIscTLTufU0bDeyO668/Uzv/Sfa7BXJt1jZ6dvIi8oTKNOGmt7PV2PY1uRDmz/jKswiKHMB/ge0Y1msg19OUynNlNWVuCAmB+HhjkY8dkwQuHF5K3fj61vJMuP4BpzyLM7t+e4u0a31ozz6LCgsj/yMlWa1a89a5bzg7qxFXt1cg8Zp7hgk3ZQHRsU/bMrhcUyZ2zYdn3x6EXH2Naw2aU3zzYk5dT3/tgyX6kD+0Rx819tGcNQu2b7d0NLnGqqYRCmEJI1pVZtj3R3lu2xbqsIdBzd/FycPDIu1as6VSJdi6ld8ea8P4w+/TwG0rvTb+SPTGqriVukRstTPUv76Z8wnX8PPIT686VahRyI/ISJgcoqmyZyG/OQ/Bx+Uq+v1gCrw3EvLlu9N5MC2b+MUGxqYPP/wAgwfDtm08cPK6DZIELhxexwB/Pt/jwVjdg92+VdjdsCXjWlex/pt1qeXPz7heo9mz9keCNs3mjGtxfi7cih+u9mbd2rbsXKdQLiZOJriwAwBNeU7wvedbNOVXdL0GqJkzoHr1O0+Zp6spc4O3N0yYAD17wpw58Oqrlo7I7B6YwJVSM4F2wHmtdY3kY4WBRUA5IBJ4UWt9JffCFCL3rFwJDbfMoBTRlFr0B1uaNLF0SNkyonUVRt58gR2lq9M9fDXtI9bQ4/Zy/vMowhKfDlx0KkLlxH955NZxHrkZScHbsWg8ISQENXgwON97YzPPV1Pmhh494JtvICgIOneGggUf/BhbktJpL6M/QGOgLrA/1bEJQFDy34OAzx70PFpr6tWrp9M6ePDgfcfy2unTp3WHDh10pUqVdIUKFfSQIUP0rVu3tNZa//HHH7pt27b3PeaXX37RderU0bVq1dJVq1bV3377ba7HOWvWLP3mm29qrbWeNm2anjNnTobnnjhxQs+fP//O9zt37tSDBw82SxzW8J6ZS1yc1gHlLusrToV0UqtnLR1Oji3/J0o/MW69Lvd/v+omY1bpAc8F6TWPNNS3nFy0Bn0mfxH9V9laem5AG60nT9b6xAlLh5z7wsK0Vkrrd96xdCTZBoTpdHLqA0fgWus/lVLl0hx+DmiS/Pc5wEbg/3L6y8QStNZ07tyZAQMGsGLFCpKSkujbty/vv/8+EydOTPcxCQkJ9O3blx07dlCqVClu3bpFZGRktq6flJSEc5qRT1b0798/059HRkby448/0r17dwACAwMJDLxvGqnD+/xzeDHyMwqqGNRn4ywdTo6lXb7eaLwbq6o8Sf5bN9HADTdPwJh50vOtphaKMo/Vqwevvw5ffWV8taM1DNmt6vtprc8CJH/1NV9IeWvDhg24u7vzanJ9zNnZmZCQEGbOnMnNmzfTfUxsbCyJiYkUKWJsOuvm5kblyvfXBYODg+nZsydNmzblkUceubNCc+PGjTzzzDN0796dmjVrAjBv3jzq169PnTp16NevH0lJRt1x1qxZPProozz99NNs2bLlnueeNMlYJXj06FGaN29O7dq1qVu3LseOHSMoKIjNmzdTp04dQkJC2LhxI+3atQPg8uXLdOzYkVq1atGwYUP27t175zlfe+01mjRpQoUKFfjqq69y/O9rzU6dglkfR/O205eo7t2hdm1Lh2R2KXO+r7t53kneNlXHNpdPPoH8+Y0VmiaTpaMxm1y/iamU6gv0BShTpkzmJ1ugn+yBAweoV6/ePccKFChAmTJlOHr0aLqPKVy4MB06dKBs2bI0a9aMdu3a0a1bt3Q3fNi7dy/btm3jxo0bBAQE0LZtWwB27NjB/v37KV++PIcOHWLRokVs2bIFV1dXBg4cyPz582nRogWjR49m165dFCxYkGeeeYaAgID7rtGjRw+CgoLo1KkT8fHxmEwmxo8fz6RJk/j1118B45dGitGjRxMQEEBoaCgbNmygV69ed/q3HD58mD/++IPY2FgqV67MgAEDcHV1zfzf2EYNGwbvJXxEPuckGDvW0uHkCruoY5tDsWIwfjz072/spTl4sKUjMovsJvBzSqkSWuuzSqkSwPmMTtRaTwemg7GUPpvXyzVaa1Q6jY0zOp7i+++/Z9++ffz+++9MmjSJdevWMXv27PvOe+655/Dw8MDDw4NnnnmGHTt24OPjQ/369SlfvjwA69evZ9euXXda18bFxeHr68v27dtp0qQJxYoVA+Cll17iyJEj9zx/bGws0dHRdOrUCQB39wcvPvnrr79YunQpAE2bNuXSpUtcvXoVgLZt2+Lm5oabmxu+vr6cO3funr4v9uL332HvkggWOc1EDXgTkt8Le5TnXQGtVd++8MsvRr/wpk3vmXFjq7KbwH8GXgHGJ39dYZZoLNBPtnr16neSWYpr165x+vRpKlasyKVLlzJ8bM2aNalZsyY9e/akfPny6SbwtL8EUr5P3apWa80rr7zCuHH31mBDQ0Mz/SWS8tiHld5jUq7j5na3MZKzs3OmG1nYqrg4YyA2xet9lPIwOtgJ+6cUzJgBNWvCyy8bc8NT/fduix5YA1dKLQC2ApWVUlFKqT4YibuFUupfoEXy9zapWbNm3Lx5k7lz5wLGTcVhw4bRu3dvPD09033M9evX7ylJhIeHU7Zs2XTPXbFiBfHx8Vy6dImNGzemu0FEs2bNWLJkCefPGx9kLl++zMmTJ2nQoAEbN27k0qVLJCQksHjx4vseW6BAAUqVKkVoaCgAt27d4ubNm3h7exMbG5tuTI0bN2b+/PmAUVopWrQoBSy92WIe+ugjKHxsB61vLEUNHw6+NnsLRzwsPz8jiYeHG/to2riszELplsGPmpk5FotQSrF8+XIGDhzI2LFjMZlMtGnThk8//fTOOevXr7+njLBgwQImTJhAv3798PDwwMvLK93RN0D9+vVp27Ytp06d4sMPP6RkyZL3lUGqVavGxx9/TMuWLTGZTLi6ujJ16lQaNmxIcHAwjz/+OCVKlKBu3bp3bm6m9sMPP9CvXz9GjRqFq6srixcvplatWri4uFC7dm169+59T+08ODiYV199lVq1auHp6cmcOXNy+K9oO3bvhkkTNYd93wZ87aZbnXgI7dsb5ZSJE42mV08/bemIsk3ayeaijHa5t3W2+p4lJkLDhhD47wK+vdbdGIm99pqlwxKWcOMGBATArVuwZw+k2rDcGkk7WWFXHriDejq++goO7brBF67vGnODe/fO/UCFdfLygnnzIDramFpoo6QXSi4KDg62dAh26YE7qKfjxAn48EP4/pEJeP4bBSsW2mVzI/EQ6teH0aONWni7dtC1q6UjemhW8V9wXpZxRM5Yw3uV3p6OmfWp1tqYdVKWk7x0aoKxa3mjRnkRqrB2I0fC448bKzTXrrV0NA/N4gnc3d2dS5cuWUViEJnTWnPp0qUszTXPTRn1o87o+OzZxv+bKyqPwMlJwWef5WJ0wqa4uMDSpVCxorGTT/LsLLM6eBCefdYo15iZxUsopUqVIioqigsXLlg6FJEF7u7uubqwJ3R39ANXDT5Mn+o9e2DgQHirziYe2b3YmENYunSuxS9sUIkS8Oef0LGjMT/83DnzzE66dQvGjYNPPzVa2x45Av7mXVBl8VkowjFkJTGnrW2D0bcj7Ua9WT3vyhUIDITbcUkcL1wP19grcOgQoRFXZGm5uN+tW0bv8MWLjT4LEyZk/z7Jli3wxhtw6JDR0jYkxFjOn00ZzUKx+Ahc2L4HJees3nTMrLad+rys9PcwmaBXL6NhVcTgabiG7IFFiwiNuPLQN0CFg3BzgwULjMU+n39ubIg8cybky5f157h2zeg9Pm0alC0Lq1YZ5ZNcIglc5EhWknNWE/PD1LYf1N9j3Dj49VeY995BKnwxAlq1ghdeYOJnf2QpFuGgnJ2N+aYlSxobIkdGGmsFWraEjEqHiYmwdauRrOfOhbNnjcZ8Y8caHRBzkSRwkSNZSc5ZTczm2oNx7VpjymCvl27RfWV3o/44ezYo9dA3QIUDUsqYneLvb4ym+/QxjlerZgwEWrY0eor/8YeRtNeuhatXjRuijRvDsmXQoEGehGrxWSjCtmUlIWaUgNMeT+ldndrD9q4+cQK6dzcazX1fbCRqzx7jY3Dx4g8VixD06mXMHNm7FyZNMhL6N98YJZFy5Yw9Nv/6C55/HpYsgYsXYf36PEveICNwkUOFTYU48U8B8p+EK175cSkYj7N3PMVLmIiMNP47z+rmuDntXR0ebrS2SEqC1UNX4/p6iLHKLnkjC7CDjXpF3lLK6F5Ys6ZxYzMuDjZvhsOHjdF27drGOZYKT2ahiIcRujuaz746iu/Wk9SNOkrdG7tpyDaKcZH/8OMvnuQvnmQzT7GH2lSv5ULPnlCo5llm7j6UazM/1qyBLl2gUCFYO+88VV6sZdz137EDPO4dXWdlRowQ1iSjWSiSwEWWLd4WxY4XZvFx1Me4cRuACK8KXK5Zm398ilLy9FEaRB+iZMx/ANxyL8CCwoN458wwYlRhmjY1Zmm98AJk0Kk3W2bNMmZs1agBq1ZqSvbvAOvWwc6dxshJCBsn0whFjpz59wZJzYYx8eZPrPZpwg9NW7CvdEWuuefH38eDLUF3N8hdvXonW2evoH74JnpHfEo3z6/5+7FhDD32Nr17F+Dtt40b+wMGGAvgsktrY13ORx9BixawZLGmwKRRxvSTL7+U5C3snozAxQOFLTyKZ8/OVEncz/iq/fmu/bNodff+twJOjDf2+kw7rbDK+RMM//tHmkdsRRcuzInnRzDq4hAW/uyJyWTcDxo4EJo3z/rmKPHxxnTdr74y6t69e8P0b024vvu2cfDVV41WsRasTQphTtJOVjy00N3RvN5gGpW6BeKXFE3f9pOY3qHtPckb7p3BkXZa4WHf8rze8X1eGzQN1bAhFb4bybywKpz/aiEffqD55x/jHmPhwkaf/alT4fjx9OM5c8aYHlimjDGCT0w0JpjMnJ6Ia7/XjOQ9dCh8/70kb+EQpIQi0hW6O5rF/Vfww45BhLvVZFD3d7hQvBiuChKS7n5qSzuDI6NphX94lYaVK2HTJhg6lMJvduOjx7/igyWTWRtTn99+g99+M6ofYMzYcnGBhAS4fdv4GhtrlE3at4e33oJnngF1Kx5e7AahoUYt5cMPJXkLhyEJXKRryjfh/LgjmAOuVek2YBS33NzBpPHxcMXLzSXDGRwPXIzz9NMQFgZz5hD/bhDuTzYgoVoTYpt35/Mlbajh7c/q1cbkEaWMVcyursafwoWNXkN36uYXLxo9nNevN2reQ4bcd12ZcSLsmdTAxX1ir5rYUbQZjydup223rzlRpvidn6Wud6fnYRpSjV2wndc2L6RP2ArcE28T7l+FpH79qTe8731T/+6xY4dRa1m06G4dpVevbMcihLWTWSgiS7SGVU9/xkuJG3mnftA9yRsevGIxq4txJq6J4JKTGxOffoXp9Tvz/P4N9Aj/jYqjhsLkMfDSS0bb1wIFoGBB4+vFi/Dtt8b0wPz5jSXOb75pLHFOR1Z7sAhhq3KUwJVSbwFvYAzMvtNaTzZLVMJiVgRt5fk9H7KtYid+a9EYEk13fpbVFYsPajQF99bKr3p4M/Ox55gZ2IHHT+1jQdJuY3J3fPz9D6xaFaZMMSaUFyiQ5Wtk5bgQtibbCVwpVQMjedcHbgOrlVIrtdb/mis4kbf2/hlDwIRuXPQoTf2dsxgXeT3X6sfp1sqV4lTtBhA00vg+Pt5oz3n1qvHHycnYSTyLNynN1RxLCGuVkxF4VWCb1vomgFJqE9AJmGCOwETeunFdE9WmL9WI4vqyv3AqVJCOhQrmWqkhSz1J3N3B3Z3Q6AQm/n7a+EWyNibLv0ik74mwdzmZB74faKyUKqKU8gTaAPftVaWU6quUClNKhcm2adYrtO8q2txYzOk3xuLTumGuX69jgD/jOtfE38cDBfj7eKR7c9W/croAABa8SURBVDHlRmR0TByau/3GQ3ffv79g6O5oGo3fQPmglTQavwEgS9cQwlblaBaKUqoP8CZwHTgIxGmt387ofJmFYp2iT5u4VDaA/K7XaDVkMr5FvK1mul2j8RvSLYOkXb4vM06EPcuVlZha6xla67pa68bAZUDq3zZo9nNzqKX3MumpHiQ4u2Q6ys1rWb0RmdmMEyHsVY4SuFLKN/lrGaAzsMAcQYm8E77jNl13f8xez6r8Vu9uI3prSX5Z3YBBZpwIR5TTXihLlVIHgV+AN7XWV8wQk8gjWsPGnjOoyHE+b97tvh4n1pD8srpLj+y0IxxRjuaBa62fMlcgIu+kLC+P3ubK1iNj2F6oLpuqBNx3njUkv6wuDJIZJ8IRyUpMB5Nys+9mvIkeazdTgv8Y+uwwXF2cMm1SZUlZWRiU0+3YhLBFksAdTMrNPqedPoyIm8za4o3YXroqPvlcMm1SZQuykuiFsCeSwB3MmZg4TLdcGLhlBYWI4YvWLwFwNS6B8NEtLRydEOJhyIYODqakjwduO7wYkjiFZeVbcNivwp3jQgjbIgncwbzVpDI9wtbiSRxTm3YGrKveLYTIOimhOJhrYUUYcPt//FmqEceLlsbfRuvdQghJ4A4lMRGOjFmAH+fxnTGPEy1bWDokIUQOSAnFgSxZrHn5QghXy9RAtWhu6XCEEDkkCdxBaA0bPthALfbh/cFQ2fhXCDsgCdxBrF4NHY6HEOddDKeePSwdjhDCDCSBO4h5H0bQjpXkGzLA2ChBCGHzJIE7gC1boNGuL0l0zofz4IGWDkcIYSaSwB3AlDGXeYU56G7dwc/P0uEIIcxEphHauYMHoeza6XhxE0ZkuFmSEMIGyQjczn0dksBgpnD7qWZQq5alwxFCmJEkcDt26RJcmR2KP9EM9H2CRuM3WMU2aUII85ASih0bGnyVNxK/JcrLjw0VAjEl73UJyNJ5IeyAjMDt1O3bED7rFM3YwIK6LTE5GduSWctel0KInJMEbqeWLIHXbszktnLhp1r39vm2hr0uhRA5JyUUO6Q1fDPpJr+o2ayu/DgX8he65+fS+1sI+5CjBK6Ueht4HdDAPuBVrXW8OQITDy9ls+Lj+91pvTuCQsSw5LF295wjvb+FsB/ZLqEopfyBIUCg1roG4Ax0NVdg4uGkbFYcHRPHtbDyDFDf8G+RspTp1Bp/Hw8U4O/jwbjONeUGphB2IqclFBfAQymVAHgCZ3IeksiOlM2KE696UCUimvqEMTqgH39EXGRLUFNLhyeEyAXZHoFrraOBScAp4CxwVWu91lyBiYeTcmMy9p+yDGAaN13cWVajqdywFMKO5aSEUgh4DigPlAS8lFIvp3NeX6VUmFIq7MKFC9mPVGSqpI8HptvOOIUXoLv6kRXVGhPr5iU3LIWwYzmZRtgcOKG1vqC1TgCWAU+kPUlrPV1rHai1DixWrFgOLicyM6JVZW4dLEPP2wvw1PHMC2gjNyyFsHM5SeCngIZKKU+llAKaAYfME5Z4WO1r+aP2P8JA128IL/EoMVVqyg1LIexctm9iaq23K6WWAP8AicBuYLq5AhMP5+efoVr0ZqoQAZ/OYktvuXEphL3L0SwUrfVoYLSZYhE58Pnn8IHnVLRHEdRLL1k6HCFEHpCl9HZg+3aI3BJFy/gVqD59wENuXArhCCSB24GQEHjL7X8obYL+/S0djhAij0gCt3EnT8LPi2/R33k6qm1bKF/e0iEJIfKIJHAb9/XX0FkvxfvmeRg0yNLhCCHykHQjtFGhu6MZt+IoO79+gm0ek7lepBz5W7SwdFhCiDwkI3AblNK46t/NRal9ez/1b+5kSpWWhO45a+nQhBB5SBK4DZq4JoKbt0xc21WeIV6fE+fixo/VmspOO0I4GCmh2KAzMXHcOFAK76u3eMl5CctrNOGae35ipXGVEA5FRuA2qIS3J1f/rsQb+b/BM+kWP9RtC8hOO0I4GkngNqi+qTamGHcG6G/Z6V+NQ74VpHGVEA5IEriNSUqC3+YVpm+JJZS/EcWceu1lpx0hHJTUwG3MwoVw5IhmW8VJUKkSU5Z9As7Olg5LCGEBksBtSFISjB0Lr5ffQKFjYTB9uiRvIRyYJHAbsngxRETApprjoEQJ6NXL0iEJISxIauA2wmQyRt8vVdiJ37718M474OZm6bCEEBYkI3AbsWQJHDwIq+t/Bpd9oF8/S4ckhLAwGYHbgKQkGDMGni1/mFI7lxlNq7y9LR2WEMLCJIHbgG+/hQMHYGq5iSh3dxgyxNIhCSGsgCRwK3fuHLz/PnR9Mopyf/0AffpAsWKWDksIYQUkgVu5d9+FmzdhasUvUCYTDB9u6ZCEEFZCErgV27wZ5s6FMf2iKbz4f9C9O5Qta+mwhBBWItsJXClVWSkVnurPNaXUUHMG58gSE+HNN6FMGRh+4f+MO5nBwZYOSwhhRbI9jVBrHQHUAVBKOQPRwHIzxeXwvv4a9u2DjZ9sweX9+fDBB1ChgqXDEkJYEXOVUJoBx7TWJ830fA7tzBkYPRratk6i8ZLBULo0BAVZOiwhhJUxVwLvCixI7wdKqb5KqTClVNiFCxfMdDn7NmwY3L4NMxvNQO3eDZMmgZeXpcMSQliZHCdwpVQ+oAOwOL2fa62na60DtdaBxWT62wNNm2Z0HBzz9hV8J78HTz8NL7xg6bCEEFbIHCPwZ4F/tNbnzPBcDm3jRmONTps2MDx2FFy5Al99BUpZOjQhhBUyRy+UbmRQPhFZd/w4dOkClSrBwvf34fTUNzBgANSqZenQhBBWKkcjcKWUJ9ACWGaecBxTbCw895zRcfDn5Ul4jxwEPj5GAxQhhMhAjkbgWuubQBEzxeKQTCbo2RMOHYLVq+GR796FP/+EmTOhcGFLhyeEsGLSTtaCtIb33oMVK+DLL6H5yRnwxRdGt8FXX71zXujuaCauieBMTBwlfTwY0aqy7H8phJAEnhuyknBv3IA33oAFC4yvg2v/CS0GQMuWEBJyz3ONXLaPuIQkAKJj4hi5bB+AJHEhHJz0QjGzlIQbHROH5m7CDd0dfeecI0egYUNYuFBTpuUxDqkZXG3Tnlj/MrBoEbjc/b06cU3EneSdIi4hiYlrIvLqJQkhrJSMwM0ss4TbMcCfFSuMrSy1UxKluv6Dt18kM34YAyYTL7QZSeDG0/xx+J87o/fomLh0r3Mmg+NCCMchCdzMMkqsp06bGDrUqHUHBoKp6XaS4k/y9fIJlL8STa8Xx3DYuzgR206hkx8THROHgjvfp1bSxyO3XoIQwkZIAjez1KNmreFWdCFid5Uj7khxvtLQv79R4u7X63cm/PYlhW5eY2TrwWwtW9t4TJrn03BfEvdwdWZEq8p58XKEEFZMEriZvdOsMsOnR3I1siA39pXm9rmCOLkn0KHHDb4I9qZC8Zsw4v+Ys3gKR4qU4bUuozngVzHT59SAv4+HzEIRQtxDEngOaA1RUbBrF2zfDlu3ws6d/ty8aSRX16KxVOh4mE9GFKDr4yVgy1/Q9g04fJij3V/nxdLtiEn1FmRULvH38WBLUNO8eVFCCJshCZz0p/0B9x1rX8ufjRuNniVhYUbiTmmw6OICAQHGlpWPP27MMilXWKPWH4SZq6DLKjh7Fvz9Yd06KjVvTnCa6z5TpRhLd0XfcxNUyiVCiIwordMb8+WOwMBAHRYWlmfXy4q086wBXJ0UKEhIMv5tEi57ovb5UuXQJapdDceP//DwiqO0XxKPllb4FUmkqFccLvHX4Xryn2vXjK3kExOhQAFo1croUtWpExQsmGk8smhHCJGaUmqX1jow7XGHH4GnN+0vwaRBax7ffIxme3ZT9+ZearIPVxIBSFJOJMY7kxjlguuFfOTL5wIeHpA//90/pUsbi3LatoUnngBX1yzF0zHAXxK2ECJLHD6Bpzftr/T5c7y/ZD6tYzdwRfkQ7lOVaRVeYH/ZCuwp8QjnvIveOTer9WkZWQshzM3hE3jqaX/uCfH0W/8LA/cs4BZuvFd2JIueb0iSq3OGj8/KghpZDi+EyA0Ov5R+RKvKeLg688SJcNZ+M4S398xhqUsnmnSaz49dG91J3hltqZCVBTWyHF4IkRscPoF3DPDn23JxzFw8luvxBWlfNpTEtSGM+rAu/j4eKIwySY+GZfBIMxL3cHXmmSrFaDR+A+WDVtJo/IZ7ep6kyGiULsvhhRA54fAlFA4d4rFBfTipy/LbyM38/EmROzuYpS1vBJYtnOm0v4xKIxn1NJHl8EKInHDsBH72LNcbP8v1+HzMenE141Il7/SknSHSaPyGTBtXpRjRqvJ9UxVlfrcQIqfsOoFnOvMjNpYbTdrCxYt8GLCJqT+Uf+i9g7NaGkm5psxCEUKYk90m8ExnftTwJa5dF9yO7OV1v1+ZtLYe+fI9/DUepjQi87uFEOZmtzcxM5r5EfzzARY/+RIef66ln8u31J0SQNGiGTzJA6TMYElNSiNCiLyS013pfZRSS5RSh5VSh5RSj5srsJzKqLxR6sQhnt8RypcM5teOtZi2Z3e6M0eyomOAP+M617xntsq4zjVlpC2EyBM5LaF8CazWWndRSuUDPM0Qk1mkW97Qmg9XzeQCxZj0WG88Kp4lLoH7bjo+DCmNCCEsJdsjcKVUAaAxMANAa31bax1jrsByKr3yRocDm2hwYS8fuAWjnjx/57jMxxZC2KKclFAqABeAWUqp3Uqp75VSXmlPUkr1VUqFKaXCLqT0Xs0DacsbFT1h5PofCKMeK1rUxynf3fq4zMcWQtiinCRwF6AuME1rHQDcAILSnqS1nq61DtRaBxYrViwHl3t4HQP82RLUlBPj27Iydisl4s8xotjHeFT77845ctNRCGGrcpLAo4AorfX25O+XYCR063P8OM4hk5hHD9qF1KVUIbnpKISwfdm+iam1/k8pdVopVVlrHQE0Aw6aLzTzie0/HKckZ/Z0Hc/EHr4MQ7YnE0LYvpzOQhkMzE+egXIceDXnIZnZ+vV4r1vOWLePGT65lKWjEUIIs8lRAtdahwP3bfNjTa68+T4xlMM7eBh+fpaORgghzMduV2IC6K3bKBSxnblFhzHwHXdLhyOEEGZlt71QAM6/F4IbBSkzqne2ep0IIYQ1s98R+KlTFN20lPmeb9DtjfyWjkYIIczObhP4pY+moDUk9BuMu1RPhBB2yG5KKKl7f1f0gOU/fMcKp870GFnG0qEJIUSusIsReErv7+iYODTQcMNqvBNi2Ny4P3m8+FMIIfKMXSTw1L2/lTbRe8evbKc+fwZIjxMhhP2yiwSeuptgkyO7qBQfyVS/17iS77IFoxJCiNxlFwk8dTfB3n+u5jSl2NCkmnQZFELYNbtI4Cm9vyufO8HTl7czLX8f8le6IV0GhRB2zS5moaR0E7z90nfcwJN1zVsy/vmy0mVQCGHX7CKBA3SsUpgbx9exyvMF/v7pCVxdLR2REELkLrsooQBcnLECr6RYbnR5RZK3EMIh2M0IPObLOdygDE1GP23pUIQQIk/YxQjcFH2W8kfX8le5npSrYBcvSQghHsgust3xj3/EGRMF3uxp6VCEECLP2EUCz7dwLjudG9D8TZk2KIRwHDZZA0/duKpRzH/Mi9nL309N5TFZtyOEcCA2l8BTGlel9D55YvVGbuPKhW6yUbEQwrHYXAkldeMqZ1MSz59ew8p8rVkQc8bCkQkhRN6yuQSeunHV4+GH8DNdYEnllpy9GpfJo4QQwv7kqISilIoEYoEkIFFrnes71Jf08SA6OYl32rWZixRh61MVKOkj2+4IIRyLOUbgz2it6+RF8oa7jau8b9yg7eXfWezTAbfCTtK4SgjhcGzuJmZKg6p/Bk7HnVtsaNSMcZ1rSuMqIYTDUVrr7D9YqRPAFUAD/9NaT0/nnL5AX4AyZcrUO3nyZLavl9q+Ys/gceUs5eIO4eKqzPKcQghhjZRSu9KrcuS0hNJIa10XeBZ4UynVOO0JWuvpWutArXVgMTNtUHk14j+qX9zEsce6SvIWQjisHCVwrfWZ5K/ngeVAfXME9SAHxy7FCY3/0Bfy4nJCCGGVsp3AlVJeSinvlL8DLYH95gosM54rf+JIvupUf7F6XlxOCCGsUk5G4H7AX0qpPcAOYKXWerV5wsrYud1nqBmzmajHX0RJ9UQI4cCyPQtFa30cqG3GWLLk4NilPIOm3AgpnwghHJvNrcQstO4n/nWvSYW2VS0dihBCWJRNJfDjm6Opc/0vzjV50dKhCCGExdlUAj/y6RIAKgVJ+UQIIWwmgWsNvht/4qhXbYo/LcvmhRDCZhL4gdWnqRv/N5ebS/lECCHAhhL48QlG+eSR96R8IoQQYCMJ3GQC/y0/caxAAIXqP2LpcIQQwirYRALfsfgk9RK2EdtGyidCCJHCJhL46clG+eRRKZ8IIcQdVp/AQ3dHc8rtFL/7PkGLlScJ3R1t6ZCEEMIqWPWGDnd2oG/YEhq2hJg4Ri7bByAbOAghHJ5Vj8BT70CfIi4hiYlrIiwUkRBCWA+rTuCpd6DPynEhhHAkVp3AS/p4PNRxIYRwJFadwFN2oE/Nw9VZdqAXQgis/CZmyo3KiWsiOBMTR0kfD0a0qiw3MIUQAitP4GAkcUnYQghxP6suoQghhMiYJHAhhLBRksCFEMJGSQIXQggbJQlcCCFslNJa593FlLoAnMzmw4sCF80Yji2Q1+wY5DU7hpy85rJa62JpD+ZpAs8JpVSY1jrQ0nHkJXnNjkFes2PIjdcsJRQhhLBRksCFEMJG2VICn27pACxAXrNjkNfsGMz+mm2mBi6EEOJetjQCF0IIkYokcCGEsFE2kcCVUq2VUhFKqaNKqSBLx5MXlFKRSql9SqlwpVSYpePJDUqpmUqp80qp/amOFVZKrVNK/Zv8tZAlYzS3DF5zsFIqOvm9DldKtbFkjOaklCqtlPpDKXVIKXVAKfVW8nG7fZ8zec1mf5+tvgaulHIGjgAtgChgJ9BNa33QooHlMqVUJBCotbbbxQ5KqcbAdWCu1rpG8rEJwGWt9fjkX9aFtNb/Z8k4zSmD1xwMXNdaT7JkbLlBKVUCKKG1/kcp5Q3sAjoCvbHT9zmT1/wiZn6fbWEEXh84qrU+rrW+DSwEnrNwTMIMtNZ/ApfTHH4OmJP89zkY/+HbjQxes93SWp/VWv+T/PdY4BDgjx2/z5m8ZrOzhQTuD5xO9X0UufSPYWU0sFYptUsp1dfSweQhP631WTD+RwB8LRxPXhmklNqbXGKxm3JCakqpckAAsB0HeZ/TvGYw8/tsCwlcpXPMuus+5tFIa10XeBZ4M/mjt7BP04CKQB3gLPC5ZcMxP6VUfmApMFRrfc3S8eSFdF6z2d9nW0jgUUDpVN+XAs5YKJY8o7U+k/z1PLAco5TkCM4l1xBTaonnLRxPrtNan9NaJ2mtTcB32Nl7rZRyxUhk87XWy5IP2/X7nN5rzo332RYS+E7gEaVUeaVUPqAr8LOFY8pVSimv5JsfKKW8gJbA/swfZTd+Bl5J/vsrwAoLxpInUhJZsk7Y0XutlFLADOCQ1vqLVD+y2/c5o9ecG++z1c9CAUiebjMZcAZmaq0/sXBIuUopVQFj1A3GxtM/2uNrVkotAJpgtNk8B4wGQoGfgDLAKeAFrbXd3PTL4DU3wfhYrYFIoF9KfdjWKaWeBDYD+wBT8uH3MGrCdvk+Z/Kau2Hm99kmErgQQoj72UIJRQghRDokgQshhI2SBC6EEDZKErgQQtgoSeBCCGGjJIELIYSNkgQuhBA26v8Bb/N70jw/5M4AAAAASUVORK5CYII=\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           4.961530\n",
       "x1                  0.494829\n",
       "np.sin(x1)          0.567696\n",
       "I((x1 - 5) ** 2)   -0.019033\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.098670\n",
       "1    10.955380\n",
       "2    10.686285\n",
       "3    10.342119\n",
       "4     9.989667\n",
       "5     9.695411\n",
       "6     9.509256\n",
       "7     9.452310\n",
       "8     9.511716\n",
       "9     9.643798\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
}
