{
 "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.988\n",
      "Model:                            OLS   Adj. R-squared:                  0.987\n",
      "Method:                 Least Squares   F-statistic:                     1287.\n",
      "Date:                Mon, 24 Feb 2020   Prob (F-statistic):           2.32e-44\n",
      "Time:                        22:49:06   Log-Likelihood:                 10.116\n",
      "No. Observations:                  50   AIC:                            -12.23\n",
      "Df Residuals:                      46   BIC:                            -4.584\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.0074      0.070     71.295      0.000       4.866       5.149\n",
      "x1             0.4960      0.011     45.791      0.000       0.474       0.518\n",
      "x2             0.5148      0.043     12.089      0.000       0.429       0.600\n",
      "x3            -0.0202      0.001    -21.285      0.000      -0.022      -0.018\n",
      "==============================================================================\n",
      "Omnibus:                        0.999   Durbin-Watson:                   2.043\n",
      "Prob(Omnibus):                  0.607   Jarque-Bera (JB):                0.631\n",
      "Skew:                           0.274   Prob(JB):                        0.730\n",
      "Kurtosis:                       3.041   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.50134616  4.98737075  5.43308135  5.81042435  6.10147053  6.30136081\n",
      "  6.41910454  6.47609935  6.50261552  6.53282271  6.59917592  6.72708331\n",
      "  6.93073163  7.21075541  7.55413274  7.93632467  8.32530693  8.68683143\n",
      "  8.99005311  9.21259769  9.34423772  9.38857291  9.36243871  9.29314028\n",
      "  9.21396636  9.15872     9.15616482  9.22529923  9.37223503  9.58919331\n",
      "  9.85578243 10.14234835 10.41484672 10.64043624 10.79287438 10.85682901\n",
      " 10.83039738 10.72541873 10.56552919 10.38227835 10.20994584 10.07990921\n",
      " 10.01548874 10.02811594 10.11545543 10.26178859 10.44059613 10.6189151\n",
      " 10.76275462 10.84268049]\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.82514732 10.6718622  10.40301182 10.06459909  9.71718005  9.42103759\n",
      "  9.22142214  9.13747276  9.15753126  9.2419965 ]\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXhM1xvA8e/JIostxJ7Y21K7yM9e+661VdVSpdXaqqiWRlVRlKLVTSlKUfuWqlbtamnR2JXag4QSIYjsyfn9cZM0iQmRTDKZyft5Hk/izp25781N3jlz7jnvUVprhBBCWB87SwcghBAifSSBCyGElZIELoQQVkoSuBBCWClJ4EIIYaUcsvJghQoV0mXKlMnKQwohhNU7dOjQLa114ZTbszSBlylTBj8/v6w8pBBCWD2l1GVT26ULRQghrJQkcCGEsFKSwIUQwkplaR+4KdHR0QQEBBAREWHpUEQaODs74+npiaOjo6VDESLHe2wCV0otAJ4Hbmqtq8RvewkYDzwL1NZap/vOZEBAAHnz5qVMmTIopdL7MiILaK0JDg4mICCAsmXLWjocIXK8tHSh/AC0SbHtJNAF2J3RACIiInB3d5fkbQWUUri7u8unJSGyice2wLXWu5VSZVJsOw2YLelK8rYecq2EyD4y/SamUqq/UspPKeUXFBSU2YcTQmRjwcEwZgycP2/pSGxDpidwrfVcrbW31tq7cOGHJhJlC/b29tSoUYPKlStTvXp1Pv/8c+Li4h75HH9/f5YtW5ZFEQpr43skkLof76LkoB3Un7wT3yOBlg7JohISd5kyMGUKbNli6Yhsg8VHoTwp3yOBTN98hmsh4ZRwc2Fk6wp0qumRodd0cXHh6NGjANy8eZOePXty9+5dJkyYkOpzEhJ4z549M3RsYXu+XH2D8Z/GcPdEA3SUIwEqjm7TI6jwdCS1qznRuDG88grY2fggXt8jgXyy7iJntxXn/uEy6Gh7XnpJMXYsVKli6ehsg1X9CvkeCWT0uhMEhoSjgcCQcEavO2HW1k2RIkWYO3cu33zzDVpr/P39ee655/Dy8sLLy4s//vgDAB8fH/bs2UONGjWYOXNmqvuJnCE6GlauhMaNYXi3otw7UpwOxVYwvNJHtKm0gmJFL3D5ZgQbN0KfPlC3Lvz1l+nX8j0SSIOpOyjr8wsNpu6wyta775FAhn4eyKHpdbm7vzzO5W5Spv8+evgEUqXSoz/dirRLyzDC5UAToJBSKgAYB9wGvgYKA78opY5qrVtnZqAA0zefITw6Ntm28OhYpm8+k+FWeFLlypUjLi6OmzdvUqRIEbZu3YqzszPnzp2jR48e+Pn5MXXqVGbMmMHGjRsBCAsLM7mfsC2mPgE2LuNBu3awfz80KHmFqaU/5ZXg1XhcSX7PJ9glHwX/V4NjjVvTcfe71KnjxBtvwCefQKFC/73+6HUnEn/PExopgFl/xzPC1M8ASLYtwK8QV9d741jwAUV7HSZX4VAKPrjD3SFfEua3EbQmJE8BnD2KU7BcSShSBIoVg9694ZlnLHyG1iMto1B6pPLQejPH8ljXQsKfaHtGJKwVGh0dzZAhQzh69Cj29vacPXvW5P5p3U9YL1PJdeTiM8T8UoQSl44QUG08JU78hgb2lq7B5Ob9OF7sKcrfDuSpW1eoHnqdDpHB1Fg9hosVfmRWk/mMWFCftWvhiy+M3JVVjZT0MvkzWH0MFETHarSGU7+VIOT3ijiVCqZIZz+KRd9gwPZ19Dq6CcfYGH6p2JAbeQpS6EEIRR7cpcrpc+Q/eBBu3jR+EN9/D926WfhMrYNV9YGXcHMh0ESyLuHmYtbjXLx4EXt7e4oUKcKECRMoWrQox44dIy4uDmdnZ5PPmTlzZpr2E9YrZXKNuedM4MratLq7mQ32L2F/0w3GjGFrnXYMO3A3cd8At2IcrFCbql2qQk0P+O037AcMYOiqhvTqMYQelz7h1VfzEBgIgXfCwcRIzcxopKSHqTeY6DijsaPjFLe3VSb0SGlcnw3Eo5Uf7+9byCtHfsU+Lpb1lZsxq95L+BdM/kbk4ebCPp9mEBAAL79s/Nu3D6ZPh1y5suzcrJFVJfCRrSske/cHcHG0T/wIZw5BQUEMHDiQIUOGoJTi7t27eHp6Ymdnx6JFi4iNNY6dN29e7t+/n/i81PYT1itlV0HSxkP0bVdurKxDl7Cf+FG/in3lqrB5MxQuTCtgiscjbra3aQMnT8IHH1Bw1jcsyreSN0p/zejR3SjSoBLODU6Rcri9uRsp6ZXaG4mOsSPop5qEny9GvjrnKVXnEPPXTaLO1ZOsqNaKWfW6cdWt2KNf09MTdu2CUaOMlvhff8GqVcZ2YZJVJfCEPwBzj0IJDw+nRo0aREdH4+DgQO/evRkxYgQAgwcP5sUXX2T16tU0bdqU3LlzA1CtWjUcHByoXr06ffv2TXU/YZ1MdRUoQANRQXm4sbIOr0f9wJzYtzhZuhLVd+6E/PkTn9+ppsejfy/z5sX3dR9W3i/LxJ9nsv5qL/qUC2HFvv7kC3XArdVxVPwQA3M3UjLC1KdgHQe3fqlO+PliFGhxkopP/8WipeMoc+c6414ew+Iy9Sjh5kKBqBjuhEWbfM1Ejo4wcybUqwf9+kHNmrB8ObRokdmnZpVUQl9vVvD29tYpb+ydPn2aZ599NstiEBmXE65Zg6k7THbX6UgHrv3QkGFh3zAjyofd5b25u3g5L9R/6rGvmbJFHxaf0PJGPmDRqo+o+u95+j01kyVn36ZglZvkbX0Ij0JOZmmkmEvKNzatIWRrFe4dKY1bk9PULbuNH1aPxyU6kuNfLaBh/26pPheMN6cpXaqaPr9//oGuXeHUKZgzB/r3z/Tzy66UUoe01t4pt1tVC1yIrJJaV0Hw1soMDZnNDHzYXrUxofMX0rH24wt7mWrRJ7jvlJtXu03kh9XjWHBuOI6VI1lw8j1qe7bF1xecnMxzTuaQ8lNwrF8l7h0pTadX71PKcSOTlowj3Dk3fot9adatxSOf+9hP0BUrwoED8OKLMGQIeHmB90M5LEeTFrh4YjnhmplqgYee9OCpXyL4QzXErksnY+C3vX2aJpel1qJPKndkGAvWTMA78DS7+y2h2fyedO0KK1aAvb3ZTzHDvv4ahg6F11+H+T13oNq1NYYAbtpk3n7r4GCoUcN4Jzt8GPLlM99rW4nUWuBWNZFHiKwysnUFXBz/y5rRd1yJ2VKK1bl6oEp5wvz5ick7LZPL0jKK5IGTK4N6TuJ2rTo0XdCbX3v+yJo1Rs9BFraz0mTZMiN5d+oE3/U/hOrU0Ujev/9u/puO7u5GP7i/f/b8YViQJHAhTOhU04MpXari4eYCsYp7v9Zijn4bz5grqGXLwM0NePS47aRSG0Xi5uKIh5sLCmM43Ufd/0fh37dB48a0XdmXRT1+Y8ECeO+97JO3fvvNmE3auDGsmHgOhxfaGkn2t9+gYMHMOWjDhvDxx8annvnzM+cYVkj6wIVIRcJIkvffh8CAH+nOMpgwAerXT9wnrZPLUhsCO75DZdN9wL6+0LgxvTd05cbLvzPq81oUKAAffmiec0uv/fuNLukqVWDDd9dxatPKeGfZsgU8MvlGq4+PMcxw6FBjlIoUVJEWuBCPsn07rJl2gXmOg41W4AcfJHs8tZZ1yu1JW/QJre1UR1+A0c/766+oQoV4b1d7RnS+xNix8NVX/+2S1TVTTp2C9u2heHHYvDKEfN3aQFAQ/Ppr1kx/t7ODJUuM4ZovvwwPHmT+MbO5HN8CDw4Opnnz5gD8+++/2Nvbk1D29uDBg+SSmWA5VkQEDH4zmrXOvXB2soOlS8Eh+Z/Mk0wue+zY8JSKF4dNm1ANGjDj7zbcarePYcMKYW8PHvWztmbKlSvQqpUxMXLrhnCKvNkRTp+GX36B//3P7MdLVdGi8OOPRjBDhxrT7nOwHJ/A3d3dE0vJjh8/njx58vDee+8l20drjdYaO1uv/ymSmTYNul/6hBocgMWroFSph/bJrMlliZ59Fn7+GdW8OQsLdiCs/XaGDHGhbIdQ4p7Nmpopt24Z+TI0FHbviKHsmJ6wZ49xJ7NlS7MeK01atDA+CU2eDM2bQw4u6SwZKRXnz5+nSpUqDBw4EC8vL65evYpb/I0rgBUrVvDGG28AcOPGDbp06YK3tze1a9dm//79lgpbmMnFi7Bi8gXG2E2B7t3hpZdS3bdTTQ/2+TTj0tT27PNpZv4WcIMGsGwZdgf2s1J1p/Pz0VzaUIH7Rx5+QzF3zZRbt6BtW7h8GX7eoKk2Z7DRP//ll8bPxVLGjzf6wYcNgzt3LBeHhWWrFvjw4RDfGDabGjWMsgrpcerUKRYuXMicOXOIiYlJdb+hQ4cyatQo6tati7+/P88//zwnT55MZ8QiOxg+HGbEDsfBxRE++8zS4eBbtg7nX3ib9zZ8Rf9q7djxzBJub6kKQN6aVxL3M2fNFH9/aN3aSN5r18JzWz+CefOMpXXefttsx0kXBwf49luoVctI5l9+adl4LCRbJfDspnz58vwvDf1727Zt48yZ/4aN3blzh/DwcFxcskcBIvF4SSfjOF/3pPzPV2nHRhg/HUqUsHhso9edIPzZVsSFhDBq92JmVe9H/6e+5/aWqugYe/J6X8I1l/lqphw7ZtTcioiAbdug4dFvYNIkeOMNmDjRLMfIsBo1YMAAmDUL3nwzR45KyVYJPL0t5cyStCCVnZ0dSWetRkREJH6vtZYbnlYs6TT3uGg7rvqW5Be7Ttwr9RT5hg61dHjJxpp/W68beSPDGHRgDRH1hzPK6Stu76iEXZA7n3wTS6eaGX+z2bnTmKCTLx/s3QuV/17136yd2bN5qFSiJU2caIwNHzrUGDKUnWLLAtIHnkZ2dnYUKFCAc+fOERcXx/r1/61n0aJFC2bNmpX4/6Pm7gcSZpVy+N2En/9OTJD3DpRn+L1vKRfnz4fNB2SLetQp+7U/bdyHJTXb0e+PlQT1/J5p0yDkVFHGvV6Cw4fTfxytjfuSbdpAyZLwxx9Q+fAS6NXLGEK5bNlDo3Aszt3d+GSwcyesWWPpaLKcJPAn8Omnn9KmTRuaN2+OZ5LpwrNmzWLfvn1Uq1aNSpUqMW/ePAtGKR7F1NT3hBKn0Xdcyf+nAx+oyWys0JCfC2WPei8P9WsrxUctB7K5RnPsxnzASDWD33dpIiON+3rffvvkszYPHzYGlPTqBbVrw57dmpLLp8Grr0KjRrBxI2TXLsH+/aF6dXj3XQgLs3Q0WUqKWYknZs3X7FFFpW6urcWPFwbRxv43mr85G7tSpYyVYiwstTKsU1+oSMfpI42WZ58+3Jo0h1f7O7NpkzFwpU8foxprgQKpv7a/vzG7c+lSY13OsWNhYP84co1+1+jTfPllWLQoe5VENGXPHuONZuxYY8q9jZFyskKQ+jC7cH93Gp4/Thd8mVb/VULcizElmyyikNpY8441PYz+34kTYfx4Cv39NxvXrGPWhpLMmmU0TIcMgeefh1deMYaxX7sG168bXy9eNJ5uZ2cMqx41CvI7R0LfvkYJxGHD4PPPjR2yu+eegx49jMH7r70GZR9f4tcWSAtcPDFrvmamWuA6ThH8Qz0OBjfANV8Ifd9dwPD2j5jmnh1t2GBkaWdnWL0a3agxhw4ZkxaXLzfWC05KKShcGF54wRiF5+mJMfTk7beN1uynn8LIkdZ1UzAwECpUMPqC1mf5muuZKt3lZJVSC5RSN5VSJ5NsK6iU2qqUOhf/9REf0oTIPlKWiQWIPFmaV4OWUynuNGUWfsuuD9tYV/IG6NABDh40qgG2aIH6YDTehfz54gsjr23dasy/OXjQWDt41f5AnnpnB9sL/cKwr5dxtf2LxvJlJ08a9UZGjbKu5A1GMa0PPzROdMsWS0eTJR7bAldKNQJCgcVa6yrx26YBt7XWU5VSPkABrfX7jzuYtMBtg7Vfs6Rjvos45eX6F5U4EVaRfA2ronbssL7EldTduzBwoNE3AsYc+P79jaa2oyNgnP+YNUfJdyeINw+up9fRX9HKjiu93+SZzycllsq1SpGRxnhwBwc4fjzxnK1duvvAtda7lVJlUmzuCDSJ/34RsAt4bAIXIjtIWlRqxAjwvPcu+dRt1MyZ1p28wajUt3y50QWyYIFR7OnFF40iUOXKwa1bNA34lxPhodihiVV2rKragi8a9sShZEn2pUjeaVltKFtxcjJmznbsaKyjaekZo5ksvTcxi2qtrwNora8rpYqYMSYhssTZs7Dpq3OcsPsa9drrxsw+W1GqlNG5PXassdDCwoUQEgJeXvyUN5TbLvm445KPvWVqcMG9JAAqxb0BU+t4ZmbFQ7N54QWj4NW4cUahK3d3S0eUaTL99rJSqr9Syk8p5RcUFJTZh0uXgIAAOnbsyNNPP0358uUZNmwYUVFRAOzatYvnn3/+oeds3LiRmjVrUr16dSpVqsR3332X6XH+8MMPDBkyBIA5c+awePHiVPf19/dn2bJlif/38/NjaDaYVZidvPsuTFOjsHNxMiaD2CJ7e6OI95o1xpz4FSv47qURfNGwF4tqvZCYvAHyuzimOsEpganVhrIdpYzRM3fvGm9iNiy9LfAbSqni8a3v4sDN1HbUWs8F5oLRB57O42UarTVdunRh0KBB/PTTT8TGxtK/f3/GjBnD9OnTTT4nOjqa/v37c/DgQTw9PYmMjMTf3z9dx4+NjcU+HSvWDhw48JGPJyTwnvGlNr29vfGWFb0TbdkCoRt38gK+MGEyFCtm6ZDMwlSXByQfgti0YmHWHgpMlpwd7RQPomIICTcmNT1qAWZzVzzMFFWrGnVSZs+GQYOgUiVLR5Qp0tsC3wD0if++D/CTecLJejt27MDZ2ZnXXnsNAHt7e2bOnMmCBQsIS2VW1/3794mJicE9/qOZk5MTFSo8PGZ4/Pjx9O7dm2bNmvH0008nztDctWsXTZs2pWfPnlStalSU+/HHH6lduzY1atRgwIABxMYaf1wLFy7kmWeeoXHjxuzbty/Za8+YMQMwSt+2aNGC6tWr4+XlxYULF/Dx8WHPnj3UqFGDmTNnJvskcfv2bTp16kS1atWoW7cux48fT3zN119/nSZNmlCuXDm+Srr8iw2JiIChb8UyK9cIdMlS8M47lg7JLEzNMh25+hgj1xxLtm3toUBerOWRbHWgPM4ORMemrX1lzoqHmerjjyFPHuNGR3ZZUNTMHtsCV0otx7hhWUgpFQCMA6YCq5RS/YArQOrFkp+EBerJ/v3339SqVSvZtnz58lGqVCnOnz9v8jkFCxakQ4cOlC5dmubNm/P888/To0cPkws+HD9+nP379/PgwQNq1qxJ+/btAWO1n5MnT1K2bFlOnz7NypUr2bdvH46OjgwePJilS5fSsmVLxo0bx6FDh8ifPz9NmzalZs2aDx2jV69e+Pj40LlzZyIiIoiLi2Pq1KnMmDGDjRs3AsabRoJx48ZRs2ZNfH192bFjB6+++mpi/ZZ//vmHnTt3cv/+fSpUqMCgQYNwtJE7+QmmTYMG53+gEkdh2vLsO0X8CZlaYDk67uHEFR4dy85/gpLNMi3r80uajpHaakPZUqFCRj/4iBHGsm/xf3u2JC2jUHqk8lBzM8diEVprlImRB6ltTzB//nxOnDjBtm3bmDFjBlu3buWHH354aL+OHTvi4uKCi4sLTZs25eDBg7i5uVG7dm3Kxs8W2759O4cOHUosXRseHk6RIkU4cOAATZo0SVzi7eWXX+bs2bPJXv/+/fsEBgbSuXNnAJydnR97znv37mXt2rUANGvWjODgYO7evQtA+/btcXJywsnJiSJFinDjxo1kdV+s3fnzMGdyMP/keh9qNzSmituIJ+naSLlvCTcXk90mbi6O5HZysJ5RKCm99ZYxGmXECGNIpY01RrLXVHoL1JOtXLlyYjJLcO/ePa5evUr58uUJDg5O9blVq1alatWq9O7dm7Jly5pM4CnfBBL+n7RUrdaaPn36MGXKlGT7+vr6PvJNJOG5T8rUcxKO45Sk5oW9vf0jF7KwNlobU8sn69HkjQsxqj5Z+7DBJFJLwqntm1Rqa3uO71DZuhJ2SrlyGcMKX3jBqBs+fLilIzIrKyhykLmaN29OWFhY4oiO2NhY3n33Xfr27Yurq6vJ54SGhibrkjh69CilS5c2ue9PP/1EREQEwcHB7Nq1y+QCEc2bN2fNmjXcjJ/vfPv2bS5fvkydOnXYtWsXwcHBREdHs3r16oeemy9fPjw9PfH19QUgMjKSsLAw8ubNy/37903G1KhRI5YuXQoYXSuFChUiX758qfyEHi+rV0dPrzVrIGTzfl6LnocaPty40WVDTM0ydbRTONonf5My1Q3SqaYHU7pUTdYvPqWLlZUTSE379kbre8IEY404G5K9WuAWoJRi/fr1DB48mIkTJxIXF0e7du345JNPEvfZvn17sm6E5cuXM23aNAYMGICLiwu5c+c22foGqF27Nu3bt+fKlSuMHTuWEiVKPNQNUqlSJSZNmkSrVq2Ii4vD0dGRWbNmUbduXcaPH0+9evUoXrw4Xl5eiTc3k1qyZAkDBgzgo48+wtHRkdWrV1OtWjUcHByoXr06ffv2TdZ3Pn78eF577TWqVauGq6srixYtSvfPL7uPFU4YlRFwI5qb8xvi5zQQXcgDNW6cpUMzu9SKXpnaZuraJJ3gZFMShhVWr26Mi58929IRmY0Us8pEqa1yb+2SXrPUyrN6uLlYvBRr0jeX29sq0ffQBr5iOAenzaH2yAEWjU1YwLBh8PXX8PvvRvVCK5LuYlZCPEpqN86yw1jhhFEZkf/mI/chJybZfcjvZb14J+YpS4cmLGHyZChTxig3++CBpaMxixzfhZKZxtv4LDBI/cZZdhgrfC0knLgoe279XJPF9v3IRSQftRzItbsRj3+ysD158hglBZo0gdGjwQbmOWSLFnhWduOIjEl5rVK7cRYWFWPxm5rF87twe0sVmtz+gx6xq5hTtyuXC5TIFm8uwkIaNzYWQP76a0gyEMFaWbwF7uzsTHBwMO7u7o8dMicsS2tNcHBwsrHmKW+c5Xdx5EFUTOI6k1l5UzPlNPKCgc+Q+29Hlji+woW8nsyu09W6JqKIzPHJJ8bEntdegxMnjJa5lbL4Tczo6GgCAgKIiJCPtdbA2dkZT0/PVGdnPslNTXOWKk05GiYqKC83FtVjl0trakfso1Pvz7j3TCXrm4giMsfevcYamgMHGvMBsrlsuyamo6Nj4oxEYf3SelPT3MMPk04jj4uyJ8jXiw/sPuG50J0wdy6b3nzziV9T2LCGDY0aOJ9/Dl26GOVnrVC26AMXtiO1/uWU203V7chIqdKENwit4faWKtS/7cf4mI/xrdQY3ngjXa8pbNykSfDMM9CvH9y7Z+lo0kUSuDCL4GCIjjZ9U9NUv7O5hx8mvEHcO1Ae179zsTzXy/gXKM6sbu/Z1HR5YUYuLrBokbFI6DvvWGXFQot3oQjrFBsLBw7ApjUPCFrzO0Wv/sV1ShBU6FnyF6lGkIM9Ma73KeV1hwn9SgBG/3hCf7ebq2Pijc6k0jtC5N2WFXhzSBR395dmS+5mFIy4Q/cXv+CtDg9XbxQiUd264ONj3NgsVsxolVvRG74kcPFEgoJg9sBj2G3+lfoPtjKWveQiSSK+Zfy7ZV+E43FVWLL/FUZtqsjdcoE4VYzDPo/R351QoyNpDer0jhCJjoYNX3sQsj+W74u+Rosbv/Np53foO6ij3LAUjzdxovGL/cknRvKeONFqkrjFR6EI63HiSAx7nxvJoAdG1cjTbhXQbZpRqV8XqF8fbtyA06eNf6dOEbv3D+zP/sMtO3fmxvVnNgMILu9MXq/LOJcNooBr+kuVJtY4CYrk/i//48GZvByu3JvKf6821kqbPt1q/ghFNhAXZ4xImTcPxozJdkk8tVEoaK2z7F+tWrW0sE6//hist9q10Br0ggqdtfdbS3Tp9zfqih9u0usPB5h+Ulyc1tu369+erqtjlJ2Oxl6vtu+svfDTDgXv64KtjuvQ0CePZf3hAF3xw03aY9A27eQRrPMSone719EatJ4+PWMnKnKu2Fit33jD+D368EPj9zebAPy0iZwqNzHFI2kN84ef5JlX/kejuN95r8kIJnTqR1CeAsBjRo4oBc2a8fHrk2k8YB7z63SipcMWDuHNz6GdeWZLKJ6e8N578OefRr96WoxfGMjV9VUI/K4pbtfD2Ze/NnXv+PHxSz7GiwmRHnZ28N13xqilSZOM1XzM1UNx9Wqm3CSVBC5StfavQF4rt4DuX9Ylj0Mo3bpPZU2dhysMPm7kyMjWFQguVIKpTV6jwaCFTGv0KrUdDvIn9fk9VwuOztxJ/fqaYsWgd29YscLojbl9G/791/jdv3ABNm6E5s3h2Ne1CTtTlPbll3Mg9/8oF3aFN178iIXlGmbWj0LkFAlJvF8/oxvFxyf9ha+0hh07oHNno4jWjh1mDRXkJqZIhe+RQFb13cxK/34cylOVQb3fJShfIZP7Pm7kSPLp9vBT2z48O+VDXvhzA9WmT2dbXDNCC5dhV+FufP7zy/T4sSZguv+xfLEHjK70GX1uLKDCucvczF2Ant0nc6xEBTykxokwBzs7mDvX+AQ5bZrx/euvw+DBUL78458fGgpLlsA338CpU8banO+/DxUrmj1UuYkpTGrSawurlr3CTWd3ug6eRISjUf9EAUl/Y1wc7TO2ckt4OKxaBStXwtatEBNDuOdTnCrfgSiX/NjlcoBcjtg5OuAeepmyu39A3b3L38WeYqHX8/xc8TkiHZ0yHocQKWkN+/YZiXjtWqOPr21bePttY4UfpYwJENeuwfXrxtcjR2DxYrh7F7y8jH27d4c0rFX7KKndxJQELh5y62Ycx0q0oF7sfl54dSYXiidf1NjDzSVzFrkNDob1641kvmsXpFiPM9rOnl1VGuH0zjBuV6vF9C1nrXexXWFdrl0zWuLffWf06xUoYHStREUl38/REbp2NRJ33bpmG8kiCVykidawuNJU+vwzmvfqv8ea55okezxLV9qJi+Mnv8uMX3eM6MhoYuzsiHB0lta2sJyoKFi3DrZtM7pGSpQw/hUv/t/XDM0CGdUAABaySURBVLa2TcmUYlZKqWHAmxifrOdprbN+WXlhVqtG7KfXPx/iV6ETvzRtBjFxiY9leSlWOzum7bjEHRzB6b/qhwkjXySBiyyXK5fRJdK9u6UjATIwCkUpVQUjedcGqgPPK6WeNldgIusd3x1CnS96EOziSa0/FzLlxWoWX6U8Oy/ZJoSlZaQF/iywX2sdBqCU+h3oDEwzR2Aia4WHaQLaD6AVVwlduwdVwI1OBdws3srNzku2CWFpGRkHfhJopJRyV0q5Au2Akil3Ukr1V0r5KaX8goKCMnA4kZk2vrWJdqGruPTax7i1rWfpcBKltbqhEDlRhm5iKqX6AW8BocApIFxr/U5q+8tNzOzpdlAs14vVpIBzGCVCTht30i3E1Co9gNlW7hHCGmXKTUyt9ffA9/EH+AQIyMjrCcvY0ncZ3eNOcGX8Cosnb1Or9EzpUjXrRr4IYUUyNJVeKVUk/mspoAuw3BxBiaxz9Xwk9X4dy9n8NegZ7W7RleTNvUqPELYuo1Pp1yql3IFo4C2t9R0zxCSy0J5ec+jJZfo0GULAvUgga1eST0pGnAjxZDLUAtdaP6e1rqS1rq613m6uoETWOLX/Hi0PTmKPWwN+r/hssscs0fJN63qaQgiDVCPMwY71+YzC3OLT1qYnJWR1y1dGnAjxZKQaYQ61/6cbvHD2M05X7cadGpUgG4y1Tl61UEacCPE4ksBzIK3h6sBJ1CKSMj9OYmSsa7LRH2C5lm+nmh6SsIVII0ngOdDe5Vfo+O93nGv0BpWqPU2n+O3S8hXCukgCz4FujPkKO+IoP3904jZp+QphfeQmZg5z8o97tPSfx5lq3XB6upSlwxFCZIAk8Bzm+LD55OcenjPftXQoQogMkgSegwT4x9DA70vOeTQmf7Nalg5HCJFBksBzkN1D11CaK+QdL61vIWyBJPAc4t5dTcVfPiMwbwWKvd7e0uEIIcxAEngO8duYPXjF+RE95B2wk8suhC2QYYQ5QHQ0uC34jBDHQpQZ+6rJmtsyhFAI6yMJPAfY/PVZ2oX/zPnuY9n1z22TNbchaysPCiEyTj5L2zitIWLKTKJVLp6a+ZbU3BbChkgCt3H7fwmm3a1FXGzQG7tiRaTmthA2RBK4jbswZgGuhFPms6GA1NwWwpZIArdhgVdiqX98Nhc8G+FSuyogNbeFsCWSwG3Yzvd/oxyXyDPqrcRtnWp6MKVLVTzcXFCAh5sLU7pUlRuYQlghpbXOsoN5e3trPz+/LDteThYZCbvztcdLHcH9/mWLrjYvhMgYpdQhrbV3yu3SArdRv82+SPOoTdzp2l+StxA2ShK4jQqdNps47Cg3tb+lQxFCZBJJ4Dbo8L5w2lxfwKXqnbHzLGHpcIQQmSRDCVwp9Y5S6m+l1Eml1HKllLO5AhPpd9RnBe7cpviktx6/sxDCaqU7gSulPIChgLfWugpgD3Q3V2AifYJvaarvm8W1ApXJ076xpcMRQmSijNZCcQBclFLRgCtwLeMhifTyPRLID4P/xFcf4pPq71HpqHE5pHCVELYp3Qlcax2olJoBXAHCgS1a6y0p91NK9Qf6A5QqJWswZhbfI4H4rDnBmEPruKfysLRGbaJWHwMF0bHGUFEpXCWEbclIF0oBoCNQFigB5FZKvZJyP631XK21t9bau3DhwumPVDzS9M1nUMfs6Ra9ltVl2vLAyZXoOJ2YvBNI4SohbEdGbmK2AC5prYO01tHAOqC+ecIST+paSDgv7duFE1Esa9rssfsKIaxfRhL4FaCuUspVKaWA5sBp84QlnlTBMDfeDFnMTrd6XChc8pH7SuEqIWxDuhO41voAsAY4DJyIf625ZopLPKFOBy9RkgB+bNAqcZujncLRXiXbTwpXCWE7MjQKRWs9DhhnplhEOt2+DU3+msd119Kcrd8AdS8qccQJyCgUIWyVLKlmAzZMOk7fuN1cHzCNPR+0fOhxSdhC2CZJ4FYqYWHiwNvhjJ21lnA7F4p/2M/SYQkhspAkcCvkeyQwcWHiXKec6RG1gvVl2+B6OZxOBS0dnRAiq0gxKyuUdGHiLvv24ko4ixs1k/HdQuQw0gK3QgnjuGNuuDIg5Af+yF+Lf4qVRcn4biFyFGmBW6GEcdwNdl2kLP4srt8m2XYhRM4gCdwKjWxdAccoZ964vJxAx2Jsq/w/Gd8tRA4kCdwKdarpQQd/F1rq7fxYrT1F3fPKwsRC5EDSB26F7t6Fuhu/IdwhD+9vnc77BQpYOiQhhAVIC9wKLZ90gS7RK7nXYyBI8hYix5IEbmXCwsDlm+nEKgeKTn3H0uEIISxIEriVWTHzOt0jFhL8fF8oIQsWC5GTSQK3IlFREPXp5zgQQ7HPR1k6HCGEhUkCtyJr5t2h1/053GjUDfVUeUuHI4SwMBmFYiViYyFo3DfkJZQ8X422dDhCiGxAWuBWYsPyB/QK/pJrXu1R1atZOhwhRDYgCdwKxMXBeZ/5FCKYojOl9S2EMEgCtwKjxgfSPXAGfxasSaM/IvE9EmjpkIQQ2YAk8Gxu+d5ruH46j5IE8GXLlwgMCWf0uhOSxIUQksCzu89HXuH9qOlsLNWE/WWMvu/w6Fip/S2EkASenZ0/D8MPfIW9imVq21eSPXZNan8LkeOlO4ErpSoopY4m+XdPKTXcnMHldPP67qOXXs6cWl0JcCuW7DGp/S2ESPc4cK31GaAGgFLKHggE1psprhxv86Y4uu0byu3cJVjcrGuyx6T2txACzNeF0hy4oLW+bKbXy9Gio2HP6wupxWHyzprOR93/h4ebCwrwcHOR2t9CCACU1jrjL6LUAuCw1vobE4/1B/oDlCpVqtbly5LjE/geCWT65jNcCwmnhJsLI1tXoFNND76dcpcXP3gGx4pPUfDUXlDK0qEKISxIKXVIa+2dcnuGp9IrpXIBHQCTM0y01nOBuQDe3t4Zf7ewEb5HAhm97kTi6vIJwwMD/R2IHTedwgShlvwqyVsIkSpzdKG0xWh93zDDa+UY0zefSUzeCR6Ews+v+TE4+gtWVWxNw213Zby3ECJV5ihm1QNYbobXyVFSDgPUGlw3uLPkbjsu5CvJxDav8yC+VQ5In7cQ4iEZaoErpVyBlsA684STc6QcBhj1ZwmWXhyIk30Eg7r58MDJFZBJO0KI1GUogWutw7TW7lrru+YKKKcY2boCLo72AIRfdGfqni+oxWGGdRzBRXfPZPvKpB0hhClSD9xCErpEJq64TEPfP+jDYmY914edT9d+aF+ZtCOEMEWm0luQd2EPqq2O5rPoUYQ274jHzImJrfIEMmlHCJEaaYFbyL59MOv5TcwK6Umk51PkWbeYTvnygZ2dybHhQgiRkiTwLJJ00o7DmbL02LCWH+PGEVWhGs6/+UK+fIDRtSIJWwiRFpLAM0HKGZZNKxZm7aFAwiLiiN1aki+P+dCBn7nYsgvlfJeAq6ulQxZCWCFJ4GZmaoblj39cIexCUTx3R7LiVlfKqkt81GwA25u+xD5J3kKIdJIE/gRSq12SVNIZlnFR9oQdL06D/ed488F42vErQU4F6dF1Mn6elVF3IyxxGkIIGyEJPI1Sq10CyWdJXrmiiQgoRpGL4XT+Zwdvxn5Paa7wr1MhvvZ6mSW12nMrdwFAhgcKITJGEngamapdEhaumTz7EuEFggjZfRzHv4/yw/1j1OAo7twG4Pdi3kyp24vtT9Uh2v6/H7cMDxRCZJQk8DS6FhJOrogY6v1xjhqXz/NMiD8Vos5TjovYYRRZjLR3IbBYJTbnf45TniXYW64G/gU9cHG05+VaHuz8J0iGBwohzMamE3ha+qwfKyaGSwt2MnXBItoEbSY/94gkF2cdy3PMvSI/F2vJv0+XwM+9MFHlyvNu20q4Ajvjj+shyVoIkUlsNoGntc86Yd+HEn2VIpzqPYWi676lbPQNCpKPn/O34acaDTjsXZpYBwcc7RQoiI6NL3N+L4rR604wpUtV9vk0y9LzFULkPDY7ld5Un7Wpyn4JiT4wJByNkei//m47p4s2oNLKcRxzqs3PfdcSffUGeXZ+zrU2VYhzcMDDzYU8zg7/Je9HHEMIITKDzbbAU6vgl3J7ykTf7uCfTNn1NVorlnZYyctru+EQ/1Pq5Jl8lmRZn1+e6NhCCGFONtsCT22IXsrtCcnWNSqcKcu+49udk/mbSjRvu4hv6xXi6Q9/ocHUHSZXxknrMYQQIjPYbAJPWm87gamheyXcXHCKjmTxvIm8fHUjn+R5jx79J3C7mkOybpXR6048lMTTegwhhMgMNtuFktDV8bhRKCNbPUNMh0F4hx6nV+k57H2pJHb20aRcfTmhbzvp89N6DCGEyAw2m8AhbZX9vBavolTgRiYVHM2+biXxLOhMYBr7z9N6DCGEyAw2ncAf587KLXh8OZJfXLsy8J/JfFhYAdBg6g6TSVz6toUQ2YnN9oE/TtSp8zi88jJ/qyqU3LqQQvHJG6RvWwhhHXJmC/zePYIbdsAxxh7/L3zpUD9Psoelb1sIYQ1yZAK/2HoQpe6cZX63rQwcVtbkPtK3LYTI7jLUhaKUclNKrVFK/aOUOq2UqmeuwDLLvc1/Um7/MpaW+oA3lja1dDhCCJFuGW2Bfwn8prXuqpTKBWTv5WW05nbfEYRSHO/V7yfOsBRCCGuU7ha4Uiof0Aj4HkBrHaW1DjFXYJnh369XU+bf/Wx5bhKVa+e2dDhCCJEhGelCKQcEAQuVUkeUUvOVUg9lRaVUf6WUn1LKLygoKAOHy6DISNRoH46rasxuUJ6yPqlPkRdCCGuQkQTuAHgBs7XWNYEHgE/KnbTWc7XW3lpr78KFC2fgcBlzacTXFA27xMdV3uGGDn3kFHkhhLAGGUngAUCA1vpA/P/XYCT0bCfu5i0KfTeJ3xxac7Bl0WSPSflXIYS1SncC11r/C1xVSiXMbmkOnDJLVGZ27tWJuMbeZ1K9N7FzjHvocSn/KoSwRhkdh/E2sDR+BMpF4LWMh2RekSfOUm7zt6x3f5O4NgXgnkyRF0LYhgyNA9daH43v366mte6ktb5jrsDM5WLvj4jAmaJzJjCqrUyRF0LYDpseCR15/ipPH1vD+tLv8FLX//q+ZYq8EMIW2HQCPzdiNs+iKTLhrcRtMkVeCGErbLYaoQ4Lx3PTXHbk6UijV8tYOhwhhDA7m03gFyYuwy0mmLA3hqLU4/cXQghrY5sJXGscZ3/FSbtqNP+4saWjEUKITGGTCfzmmt2UvnucUy2GkievNL+FELbJJm9i3vzwK+wpSJ0ve1o6FCGEyDQ21wIPP+3Ps2d92fVMf0pXlAk6QgjbZXMJ/Nw736JReEwabOlQhBAiU9lUAtehDyizdR478nehTteSlg5HCCEylU0l8DMfLSVfXAjRg2TooBDC9tlUAnf4YR4n7avR/KMGlg5FCCEync0k8JA/T/PUHT/ON+iLs4s0v4UQts9mEvjFCUuIxY6nxvawdChCCJElbCOBx8VRfOdS/sjTmsrNi1k6GiGEyBI2kcCv/Lib4lFXWF2+JeVGy2LFQoicwSYS+PlJ33OPvPxa+2lZrFgIkWNYfQKPCw3D+9xPrM/TnpiC/928lMWKhRC2zuoT+D/TNpCP+6yv2vShx2SxYiGELbP6BB6zYDFXKMmRug+vsiOLFQshbJlVViP0PRLI9M1niPD/lwOBW1hU+m1cc9sTHh2buI8sViyEsHUZaoErpfyVUieUUkeVUn7mCupRfI8EMnrdCQJDwmmz5zAOxLKseh1erOWBh5sLCvBwc2FKl6qy9qUQwqaZowXeVGt9ywyvkybTN59JbGm/eH4rf9l5cb5iXiL+CWKfT7OsCkMIISzO6vrAE25Mlrt8Ha+o46wp3Q6l5IalECLnyWgC18AWpdQhpVR/UzsopforpfyUUn5BQUEZPNx/NyZf+OMAMdizqUHNZNuFECKnyGgCb6C19gLaAm8ppRql3EFrPVdr7a219i5cuHAGDwcjW1fAxcGOToFb2e7UhHseTnLDUgiRI2UogWutr8V/vQmsB2qbI6hH6VTTg48LKMrGXubnsk3lhqUQIsdK901MpVRuwE5rfT/++1bAx2aL7BE8ftpJFI74rBqMZ9UCWXFIIYTIdjIyCqUosF4ZS984AMu01r+ZJapH0ZrS+1dxqEAL6knyFkLkYOlO4Frri0B1M8aSJpdW+1E22p9Tbcdl9aGFECJbsbphhNe/Wk0UjlQZ09HSoQghhEVZVwLXmtIHV3GoYEuKV5LuEyFEzmZVCfzCSj88oi8T1r6bpUMRQgiLs6oEfv2rVdJ9IoQQ8awmges4Tdm/VnPIvRVFK7hZOhwhhLA4q0ng55b9hUfMZSKef8nSoQghRLZgNQn8xjdG90nVD6X7RAghwAoSuO+RQOp/sp1SB1exw7Upe+8/sHRIQgiRLWTrBJ6weIPbvvOU1lfZUL6RrDYvhBDxsnUCT1i8oe3hv4jCkd/rPyurzQshRLxsvSZmwiINu57x4l+XQjxwc0q2XQghcrJsncBLuLkQGBLOoXqlOVSvdLLtQgiR02XrLpSRrSvg4mifbJss3iCEEIZs3QJPWKRh+uYzXAsJp4SbCyNbV5DFG4QQgmyewMFI4pKwhRDiYdm6C0UIIUTqJIELIYSVkgQuhBBWShK4EEJYKUngQghhpZTWOusOplQQcDmdTy8E3DJjONZAzjlnkHPOGTJyzqW11oVTbszSBJ4RSik/rbW3pePISnLOOYOcc86QGecsXShCCGGlJIELIYSVsqYEPtfSAViAnHPOIOecM5j9nK2mD1wIIURy1tQCF0IIkYQkcCGEsFJWkcCVUm2UUmeUUueVUj6WjicrKKX8lVInlFJHlVJ+lo4nMyilFiilbiqlTibZVlAptVUpdS7+awFLxmhuqZzzeKVUYPy1PqqUamfJGM1JKVVSKbVTKXVaKfW3UmpY/Habvc6POGezX+ds3weulLIHzgItgQDgL6CH1vqURQPLZEopf8Bba22zkx2UUo2AUGCx1rpK/LZpwG2t9dT4N+sCWuv3LRmnOaVyzuOBUK31DEvGlhmUUsWB4lrrw0qpvMAhoBPQFxu9zo84526Y+TpbQwu8NnBea31Rax0FrAA6WjgmYQZa693A7RSbOwKL4r9fhPGLbzNSOWebpbW+rrU+HP/9feA04IENX+dHnLPZWUMC9wCuJvl/AJn0w8hmNLBFKXVIKdXf0sFkoaJa6+tg/CEARSwcT1YZopQ6Ht/FYjPdCUkppcoANYED5JDrnOKcwczX2RoSuDKxLXv3+5hHA621F9AWeCv+o7ewTbOB8kAN4DrwmWXDMT+lVB5gLTBca33P0vFkBRPnbPbrbA0JPAAomeT/nsA1C8WSZbTW1+K/3gTWY3Ql5QQ34vsQE/oSb1o4nkyntb6htY7VWscB87Cxa62UcsRIZEu11uviN9v0dTZ1zplxna0hgf8FPK2UKquUygV0BzZYOKZMpZTKHX/zA6VUbqAVcPLRz7IZG4A+8d/3AX6yYCxZIiGRxeuMDV1rpZQCvgdOa60/T/KQzV7n1M45M65zth+FAhA/3OYLwB5YoLWebOGQMpVSqhxGqxuMhaeX2eI5K6WWA00wymzeAMYBvsAqoBRwBXhJa20zN/1SOecmGB+rNeAPDEjoH7Z2SqmGwB7gBBAXv/kDjD5hm7zOjzjnHpj5OltFAhdCCPEwa+hCEUIIYYIkcCGEsFKSwIUQwkpJAhdCCCslCVwIIayUJHAhhLBSksCFEMJK/R8CWrmE44QjbgAAAABJRU5ErkJggg==\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.007438\n",
       "x1                  0.496007\n",
       "np.sin(x1)          0.514754\n",
       "I((x1 - 5) ** 2)   -0.020244\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.825147\n",
       "1    10.671862\n",
       "2    10.403012\n",
       "3    10.064599\n",
       "4     9.717180\n",
       "5     9.421038\n",
       "6     9.221422\n",
       "7     9.137473\n",
       "8     9.157531\n",
       "9     9.241996\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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
