{
 "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.983\n",
      "Model:                            OLS   Adj. R-squared:                  0.982\n",
      "Method:                 Least Squares   F-statistic:                     893.9\n",
      "Date:                Sun, 16 Aug 2020   Prob (F-statistic):           9.05e-41\n",
      "Time:                        18:00:44   Log-Likelihood:               -0.79447\n",
      "No. Observations:                  50   AIC:                             9.589\n",
      "Df Residuals:                      46   BIC:                             17.24\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.8511      0.087     55.528      0.000       4.675       5.027\n",
      "x1             0.5194      0.013     38.548      0.000       0.492       0.546\n",
      "x2             0.4956      0.053      9.357      0.000       0.389       0.602\n",
      "x3            -0.0216      0.001    -18.278      0.000      -0.024      -0.019\n",
      "==============================================================================\n",
      "Omnibus:                        0.646   Durbin-Watson:                   2.055\n",
      "Prob(Omnibus):                  0.724   Jarque-Bera (JB):                0.543\n",
      "Skew:                          -0.246   Prob(JB):                        0.762\n",
      "Kurtosis:                       2.862   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.31050146  4.80386045  5.25769481  5.64499445  5.94849709  6.16352437\n",
      "  6.29875055  6.37477728  6.42074887  6.46956393  6.55247024  6.69393079\n",
      "  6.90760463  7.1941027   7.54088745  7.92433275  8.31360556  8.67573173\n",
      "  8.98101357  9.20790909  9.34657126  9.4004659   9.38580221  9.32886966\n",
      "  9.26171874  9.21689507  9.22209235  9.29560227  9.44330932  9.65772391\n",
      "  9.91921294 10.19922545 10.46498344 10.68486723 10.83361049 10.89645138\n",
      " 10.87155794 10.77032902 10.61552152 10.43751148 10.26930339 10.14110762\n",
      " 10.07537699 10.08311779 10.1620815  10.2971342  10.46274307 10.62717129\n",
      " 10.75769239 10.82596705]\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.79748741 10.63723969 10.36465968 10.02403908  9.67368138  9.37162707\n",
      "  9.16144338  9.06155743  9.06074468  9.12087707]\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+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2b0lEQVR4nO3dd3zN1//A8dfJjkjECLWF1h5BqkrtGq1NUfpTo61V/VaVogtdFEWrqrWqWqNqFq1NtdSI0cYKRUpiRZCQJeP8/vgkacINSe5Nbu7N+/l4eJDP/dz7eV+Xd07en3PeR2mtEUIIYXscrB2AEEKI7JEELoQQNkoSuBBC2ChJ4EIIYaMkgQshhI1yys2LFStWTFeoUCE3LymEEDbv0KFD17XWPvcez9UEXqFCBQICAnLzkkIIYfOUUv+aOi4lFCGEsFGSwIUQwkZJAhdCCBuVqzVwU+Lj4wkJCSE2NtbaoYhMcHNzo0yZMjg7O1s7FCHyPasn8JCQEDw9PalQoQJKKWuHIx5Aa014eDghISH4+vpaOxwh8j2rJ/DY2FhJ3jZCKUXRokUJCwuzdih5ytojoUzdHMSlWzGU8nZndNsqdKlb2tphiXzA6gkckORtQ+SzSm/tkVDGrQ4kJj4RgNBbMYxbHQggSVzkOLmJKYQZpm4OSk3eKWLiE5m6OchKEYn85KEJXCm1UCl1TSl1LM2xHkqp40qpJKWUf86GmPMcHR3x8/OjRo0a1KlTh+nTp5OUlPTA5wQHB7N06dJcilDkVZduxWTpuBCWlJkR+CKg3T3HjgHdgN2WDuhh1h4JpfHkHfiO3UjjyTtYeyTU7Nd0d3fn6NGjHD9+nK1bt/LLL78wceLEBz5HErgAKOXtnqXjQljSQxO41no3cOOeYye11rn+M2JKvTH0Vgya/+qNlkjiKYoXL87cuXP58ssv0VoTHBxMkyZNqFevHvXq1WPv3r0AjB07lt9//x0/Pz9mzJiR4XnCvo1uWwV3Z8d0x9ydHWlR1SdTA42cGJCI/CPHb2IqpQYBgwDKlStn1ms9qN5oyRtGFStWJCkpiWvXrlG8eHG2bt2Km5sbZ86coXfv3gQEBDB58mSmTZvGhg0bAIiOjjZ5nrBvKf/u0s5CaVHVh1WHQh96YzO/3gA9dQrKlIGCBa0die3L8QSutZ4LzAXw9/c3awPO3Kw3puwVGh8fz/Dhwzl69CiOjo6cPn3a5PmZPU/Yny51S6dLuI0n78jUQCO3BiR5xZ9/wpQpsG4dTJ8OI0ZYOyLblyemEWZWKW93Qk0ka0vXG8+dO4ejoyPFixdn4sSJlChRgr/++oukpCTc3NxMPmfGjBmZOk/YDlPzu4GHzvkOuR5HbFAJ/A+ep+TtcII9SxHsXZLjRb1YWByefBKqVbP/G6Brj4QyZVMQZw97EhvwGJHB3hQuDO++Cy+8YO3o7INNJfDRbauk+5ETjHpjyn8sSwgLC2PIkCEMHz4cpRQRERGUKVMGBwcHvvvuOxITjWt7enpy+/bt1OdldJ6wTabKG6N/+gsUxCfq1GNpSx6HD2m2TDnKiFVbeD5xOSW5YrxYNHAVoijAP3sfZSvNmdP7PYpVKERYUsR917aHG6Brj4Ty5ryzXFxVn/hrhXD0jKF465PMHF+I3o1LWTs8u/HQBK6UWgY0B4oppUKA8Rg3NWcBPsBGpdRRrXXbnAwUTNcbLbHqLSYmBj8/P+Lj43FycqJv376MHDkSgGHDhtG9e3d++uknWrRogYeHBwC1a9fGycmJOnXq0L9//wzPE7bJVHkjPun+CmBMfCKT15/h7JA9tDswkbGcIF45s71MIybWH8DJEr6Uu3WFypFX6OV9l0cvh1Jz92wily3mQ6cJLH2iDS4N/sXBxbiWpQck1vL2rDDOL30S5ZBE0fZH8ah2CeWo+W5LAr3PRUKBAlCxovGrUCFrh2uzVEqtNzf4+/vre2/snTx5kmrVquVaDMJ8+eEz8x27kcz8z0i848zg77bx7p2pXHnED8/Rg/Ho14O1F2IzHmgcP0704BEU2LON41RnlPsUAp4uReVGN212GX5KuSn0ZgzqWGWCf30U52K3Kd49ACevaBpeDKTXX1t4NmgPronx6Z9cpAhUqmT8GjIEmjWzzpvIw5RSh7TW9625sakSihC5JaP7LWk5XHNm0g/z6R3/E/+2GkD5X74GFxcAuhR9wEySGjUo8PsWWL+eSsNG8mtoB5au782Z6gvo7Gd75ZOUclN0jCZ8c22ijpWlQOXLlGh3iIGBa3jh6K/43rxMpKsHa+q1Y+vj7bgcnYhf/A36FIunZlw4nDsHO3fCihUwfjy88w44Oj784vmcJHAhTDB1v8XZQaXWwD3+gXmrx9NI7+O3Hu/Q7McPISt9YpSCTp1wa9uWxE8+pc8H4/nj0wsM/udnZi0pgqtrDrypHDJ1cxBRtxXXVjYgLrQIhRqfpmqtvXy5agr+oSfZX7Ymnzfuw7aqjYlzcTXuIXjCCSqyxtmRSf1qGd/s7tyBoUONBL57NyxZAiVKWPvt5WnSC0UIE7rULc2kbrUo7e2OAkp7uzO1Rx2mPleHaqcSWL/qNerqo2x+ax7NVnyUteSdlqsrjhPfRy//kYaOB3ljVWP6NgkmPNyibydHhVyN5+qKJ4i77E2xzofo9sgSfl30P6qGBfP+8+/wfJ/JHGz8LM6eHqk3gFOk6xtTsCAsXgwLFsCePeDnBzt25P4bsiFSAxdZlp8/s6A/b+DS2B8vxygcN67Hu00Dy734b79xt30XwqPcGFLuVz7b7sejj1ru5XNCZCSUrhnBnVBPSnbaz/hL0xh8YDXHi1fkg77j+XHai6nnZnRfQQHnJ7dPfzAwEHr2hKAgY0T+7rv5uqSSUQ1cRuBCZNKdiETCWvemlA4lcfXPlk3eAM2a4bL/D4oUd2bJhSaMe3wbJ09a9hKWdPs2PPMMxFz2okaHXfx88CUGH1jN4rrt6TNwBr1faJXu/Cz1jalVCw4eNCaMT5gAAwdCLg42bYUkcCEyQWvY1vBdnorawrmRsyne8YmcuVCNGrge/hOnyr58f6sD4xtt5cSJvNcz5c4dePZZ2L8f1n11mb3HhlLr6lmGd3qLb3qOZGLP+vfdxM2ob0yG0yZTSirvv2/8/tVXOfV2bFa+v4kZHh5Oq1bGSOHKlSs4Ojri4+MDwIEDB3BJnlUg8reNA1fR5dRkjjYYhN9nL+fsxUqXxm3vTmIbt+S7oE70avgzx/q4k+RtzIqxds+UqCjo0MFYGr9u5nnaT34arl2DrZv5snnzDJ+XrXUcShkllMOHjbX3detCo0aWfUM2TGrgaUyYMIGCBQsyatSo1GMJCQk4OeX773Pp5KXPLDccWnycKv2eILRwLSpf2oVyy6UpImFhxDVuQcKZ87R3XcfpF5xx8bmT+nBpb3f2jG350Jex5JZv0dFG8v7tN/h58gnaz2wNMTGwaRM0sHBJKa1bt8Df3wjg8GF45JGcu1YeJDXwLOjfvz8jR46kRYsWjBkzhgkTJjBt2rTUx2vWrElwcDAAP/zwAw0aNMDPz4/BgwfLEno7c+XULbwHdiXGsSAl967KveQN4OOD6+/bueRdjA1xXaj0QxJ3wzxTH85MzxRLtmCOiYHOnWHXLtgw8RDtP20KSUnGlL+cTN4A3t6werWRyHv2hPj4hz0jX8hTQ8sRI+DoUcu+pp8fzJyZ9eedPn2abdu24ejoyIQJE0yec/LkSX788Uf27NmDs7Mzw4YNY8mSJbz44osmzxe2RWsIajGERonnubBoJ5WqmtfDI1sj4RIleHPoNGbOHsnGyM60+2Ej53sVwrVURKZ6pliq42FsLHTtCtu3w8Zxf/DM1PZQuDBs20auTZWpXRvmzzdubL71FsyYkTvXzcNkBJ6BHj164PiQaUvbt2/n0KFDPP744/j5+bF9+3bOnTuXSxGKnPbbiDU0u/IjAe3HU6nfU2a9ljkj4YE9GjGw/yeEe3my6W57Hl2WQFKIT6Z6plii42FcHHTvDps3wy9vbOWZGW2gZEn444/cS94p+vSB//3PGJUtW5a7186D8tQIPDsj5ZySthmVk5NTuj0yY2NjAaNneL9+/Zg0aVKuxydyVljQDarOGsZpDz8arBxj9uuZMxI2Hm/JGwWm89nXo9h0qz1dl69jf71KTN2844EjenNbMC/fe4lhL7lw81QxBj/xJa1nvQnVq8GWLVC8eKZew+KmTTPq4C+/DDVrGlMO8ykZgWdChQoVOHz4MACHDx/m/PnzALRq1YqVK1dy7do1AG7cuMG///5rtTiF5ZxoN5Ki+jpOi7/F0c3Z7NfLykjY1JTBLnVLs3ZSLyqdCsC1WkXWJ3Xgrzf/JOiPIg8c0Wd56l4a8zddpn93D26eKsqgWpP58sDrBBavyMbPl1gveQM4Oxs9U7y8jHJKPq6HSwLPhO7du3Pjxg38/PyYM2cOlStXBqB69ep89NFHtGnThtq1a9O6dWsuX75s5WiFufaN/5Vmwd+xr9lYKnbzs8hrZnYRy0NLLSVK4LR7J+cfKc9autF04zki9lVC63uWpScz1RJgUjdjxPqgeeUHDsCwnkW4e6MAo+qOY07g2+wvV5M+PT/kkz+vWuTvxCwlS8LXXxsrNqdPt3Y0ViPTCEWW2fNnFnEhgijfmkQ7F6Jc2CFcPC0z6+TeDSLAGAlP6lbrvu3YTJU87p0yWHvEjyxaMYHal8/wCvNYVq4jxdr/hbNX7P3L0rMYy48/Qv/+4OByk69K9affqZ/ZXulxhnUZR5yTi+ml79bSrZsxhfHYMaO3uJ2SaYRCZMJfbd+iRNIl4r9ZaLHkDRmPhO+tWWe21OL5iA99e33A/nI1+JaBfH1xJBEL6uEcXOGhsWRUj5/88z+8/TY8/zy0r/kvBzwa0e/Uz3z9RHcGdXuXOCdjUVshd+e8syp01ixwcjK6GObDpfZ56iamENZ0dPoOmp6ay67HR9O8n+XnNd+7+bEpmb3pmNLu9sVeH/K/Pct5be8PPJm0j14/rqC3MvJasWKmr3HvN4OkeAfuHCnPxf2V2B8N09psZmRAHxLuxvPac++yvlLD1HOdHRRRdxO4FWPUna29KpTSpeGTT+C112D5cujdO/djsCIZgQsBxN2KwXvsYIKdHqXBrxOtFkdmbzqmjOgfKVKQmU1e4PWXp1LOO5IAp4b4rJhNyUc0zz4LixbBzZvpr5HyzUAnOBAZUIHQb1pwc2d1SvlcZF7dAbyx5RnOuhRi1w8bafX2kHQ/NRR0c3pwS1hrGDrUWEg0YgTcuGG9OKxAauAiy+zxM1vj/xZdD03luaZfcvmZalbd2izbS9/Dwozi9S+/cKHkE8y5+xKzw3sR6+xF69bGYPXKFTh57i7BFxNJuOMKSQ5UKXWI0Z6f8vy59XjEx7KyZivebTMUVcDjvjJPllrC5qa//oL69WHAAJg3z3px5BDZUk2IDCybvZvuh2ay1Ks7AU9WACuXBTJTajHJxwfWr4e5cyn3xRdMOjmIj1xf56+K3ZhxZAAHDlajfLEo2hSOwqlQGBERF+h6exkdLu4k6Yri5+rNmP94F04WT74ZaGKeurnzynNMnTowciRMnQovvghNmlg3nlzy0BG4Umoh0AG4prWumXysCPAjUAEIBnpqrW9m9BopZARuH+zpM9OJSfzp2ZiqMUG06P8NESUKpD6W2WZReZLWRj/tb781VixGRJg+r1AhGDyYJ29X57LX/UXze0fWmZ1NYxVRUcbCHjc3oyeHLe1L9xDmzEJZBLS759hYYLvW+jFge/LXNiskJITOnTvz2GOPUalSJV5//XXu3r0LwK5du+jQocN9z9mwYQN169alTp06VK9enW+++SbH4+zfvz8rV64E4OWXX+bEiRMZnrtr1y727t2b+vXXX3/N4sWLczxGW3PotUU0itnHhGoj0iVvyNpy8zxHKaMuPGeOUTdZudKYN/3990ZTqC1bjG3LLl6ETz/FoVxZky9z78g6s7NprMLDw+gZfuoUfPqptaPJFQ8toWitdyulKtxzuDPQPPnP3wG7APPXG1uB1ppu3boxdOhQ1q1bR2JiIoMGDeKdd95h6tSpJp8THx/PoEGDOHDgAGXKlCEuLi61O2FWJSYmPrTniinz589/4OO7du2iYMGCNErunTxkyJBsxWdLslo7vn0ujEpfj2avy5P83L7efY9bvSxgKW5uRjOTBzC1iXNGKzazXeLJDc88A716GTNT+vTJ/V4tuSy7s1BKaK0vAyT/nuG6WqXUIKVUgFIqICwsLJuXyzk7duzAzc2NAQMGAODo6MiMGTNYuHAh0dHRJp9z+/ZtEhISKFq0KACurq5UqXL/P/QJEybQt29fWrZsyWOPPca85Jsru3btokWLFvTp04datWqRmJjI6NGjefzxx6ldu3bqaF5rzfDhw6levTrt27dPXbIP0Lx5c1LKUZs2baJevXrUqVOHVq1aERwczNdff82MGTPw8/Pj999/T9cS9+jRozRs2JDatWvTtWtXbiZPU2jevDljxoyhQYMGVK5cmd9//90Sf8W5IivNolKWqq9vMBgPfZuNw9/D3S39WCazy83zmuzu3JOnR9ZZNX06uLjA8OF2Pzc8x29iaq3nAnPBqIE/8GQr9JM9fvw49evXT3fMy8uLcuXK8c8//5h8TpEiRejUqRPly5enVatWdOjQgd69e+PgcP/3w7///pt9+/YRFRVF3bp1ad/eqCceOHCAY8eO4evry9y5cylUqBAHDx4kLi6Oxo0b06ZNG44cOUJQUBCBgYFcvXqV6tWrM3DgwHSvHxYWxiuvvMLu3bvx9fXlxo0bFClShCFDhqTbnGL79u2pz3nxxReZNWsWzZo14/3332fixInMTP47SkhI4MCBA/zyyy9MnDiRbdu2PexvOE/IbLOolERf8+Ap+oSvYUqJ4awuoOhevzQ7T4VZZNMDa7m3Pp3VOdp5emSdFaVKwYcfGvlk9eqH/vRhy7KbwK8qpUpqrS8rpUoC1x76jDxKa41SKtPHU8yfP5/AwEC2bdvGtGnT2Lp1K4sWLbrvvM6dO+Pu7o67uzstWrTgwIEDeHt706BBA3x9fQHYsmULf//9d2p9OyIigjNnzrB792569+6No6MjpUqVomXL+2+o7du3j6ZNm6a+VpEiRR74fiMiIrh16xbNmjUDoF+/fvTo0SP18W7dugFQv379bJeFrCGzKxinbg6CyCimbJvJGVWJec+1Jj4+kZ2nwmz3hmUyS/X+tguvvmpMgn/9dWjTBjw9H/oUW5TdBP4z0A+YnPz7OotEY4V+sjVq1GDVqlXpjkVGRnLx4kUqVapEeHh4hs+tVasWtWrVom/fvvj6+ppM4Pd+E0j5Om27Wq01s2bNom3btunO/eWXXx74TSTluQ87Jytck+/cOzo6kpCQYLHXzWmZnd526VYMb/+0mvKJF+jcZC7xBR1Tj9s6S/T+thtOTsYN3CefhIkTjRa0duihNXCl1DLgT6CKUipEKfUSRuJurZQ6A7RO/tomtWrViujo6NQZGomJibz55pv079+fAgUKmHzOnTt32LVrV+rXR48epXz58ibPXbduHbGxsYSHh7Nr1y4ef/zx+85p27Ytc+bMIT65Lebp06eJioqiadOmLF++nMTERC5fvszOnTvve+6TTz7Jb7/9ltri9kbySjRPT09u37593/mFChWicOHCqfXt77//PnU0bssyu4Lx6X/+5ZXLS/mq8EsENiqZetweblhmtuNhvtGwIbzyijEwDAy0djQ5IjOzUDJqLtDKwrFYhVKKNWvWMGzYMD788EOSkpJ49tln+eSTT1LP2b59O2XKlEn9etmyZUyZMoXBgwfj7u6Oh4eHydE3QIMGDWjfvj0XLlzgvffeo1SpUpw+fTrdOS+//DLBwcHUq1cPrTU+Pj6sXbuWrl27smPHDmrVqkXlypVNJlofHx/mzp1Lt27dSEpKonjx4mzdupWOHTvy3HPPsW7dOmbNmpXuOd999x1DhgwhOjqaihUr8u2335rxN5g3ZGbH84Rbd3j/5ymcUZX4vFen1OO2esPyXlmZSZJvTJoEa9YYy+137wYT96lsmSylz0Gmdrm3Bw/7zCy5C7q50sYycdky/u/CUr4avIFlFdzyRHyWlpf+7vOMb7+FgQNh4UJjqb0NkqX0IleYOxMip2JpEBjEixeW8GXhVyg9qA576tlnUrObmSSW1K8fLFgAo0dDp06QPP3XHtjXzxN5zIQJE+xu9P0wD5oJYa1Y3GNjmLL5c06rR/m8RyembbFi5zyR+xwcjBWat27BWJteNH6fPJHAc7OMI8xzIyqOyxGxGS4UyUszIS7digGtmbhsMeUSL/LGk28TX1jlz1kZ+V3t2vDmmzB/PmzcaO1oLMbqCdzNzY3w8HBJ4jbgRlQcwSFXOHsjLsPVjhnNeLDGLi6lvN3p/+sWel5bzwfFRnPsqeIPjFHYuQ8+MBL5wIFwNQ/s62kBVq+BlylThpCQEPLiMnuR3uWIWM7eiGPW/v8aT967UMTUTAhr7eIyNvoqzwbOZqVrFxb+X3McVJLMysjPXF1h6VKjb/jAgbBhg9H0y4ZZPYE7OzunriIUeduzGTTzT1uSMDWdL/puAjej49M9J6dXCN45fJpmHw/luGNtZgwdjqNrrMzKEFCjhtEz/H//Mxb6DBtm7YjMYvVphMJ2ZHbH9HtltItLynMtPeUt6WYEl8o9geudcM4sDaBRb9OLrEQ+pTU8+yzs2gWHD4MNTGOWXemF2TK72vFeGdWcFWSqe2CWJCZy5okXKHHnLL+9tkqSt7ifUsbc8IIFjZazcXHWjijbJIGLTHtQy9G4uIw7d5pK/AruG5WbPd0wLo5/mvSnypmNLG04i+6fN83+awn79sgjxtzwo0fhvfesHU22Wb0GLmxLykKRhATY//tdTk/ZyYptH1P9+m5ilAd33H2I9SpOYtHiqNIlKTOgDZ17GT+ipq2LmyrFQPanGyZdDeOif1ceDdnDwoof0WfnEFu/PyVyWqdOMHiw0eiqXTsw0e0zr5MauMiSC8FJrHnhJ8oeXM3T8b/ixW2iHTy4WLEZOj4R14hrFIgOw/vuNVwxtqU77VydkCeeo/zI7lTqUguUynY93ZTYwyeIaNoBr6jLLG75HQM39cTZ2SJvV9i7qCioVw8iI42auImNWfICqYELs21fdYuTlTvz+t7naeX0G9dbPc+d5RsoEHWdKmc2UjV4E743D1Mi7iKuSbFEHLvIvhdmEV2wOM3++IhK3eoQ7F6VPzt+zJuPFcp0Pf1Bu8zcXLGV+AaN0FHRrHptF4O2SfIWWeDhYWz6kJgIzZvDyZPWjihLZAQuHioxEb5+7Tit53TFl/PceH8mJcYPTe3slpkGStcCr3Ls47UU2riM+nd+IxEHDpV4muV1W7KhZhV8inqZfF5Gu6B/UcOJRxetpOzm+ZxU1Qmds552g+WGpcimEyeMEorWsH27sbt9HpLRCFwSuHig69dhTquVvPF3fxLcCuK2fiVuTz+V+nhGCfZB+ykGrjvLvxMX4Xd0EWV0CLccixBSuRVuz7akfP8WONeonLrAIm2pxSkxgVbHA+i7dzNPRRwkBjfWFOxLtY2fUbepfe64InJRUBC0aAHx8UYSr13b2hGlkgQusuxCcBLr67zDq5GTuerbkOK7V6LKpE/K5tSyI28msvu9ragVy6kdtp2yhAAQ7laK8EcbohwU50Ku45p0F5eku/hGh1Ay4SrnqcA87340mvka7V4oipPciheWcuaMMRKPjoZt26BuXWtHBEgCF1kUGwvzH53M8NBxhHUdhM+yL4ylyPfIaJGOAs5Pbp/p6928oTmw7CxhK3bidWgHlaOOEI8zcbgSixuxuBHuUJhVvu3Y29SXYuUT8HB1kr7XwvLOnTNG4pGRsHUr+N+XN3Od9AMXWfJl9528EfoOIU2ep8yqrzPsGZHZvSgfpnARRdtXH4VXH0XrV7h507jk5pOXmLAhkJjEhNQQ3B1iibpLrvdWEflExYrw229GEm/a1NiWbeRIyGDbxExJTDT+QVt4RyCZhSLus3TaJfr+8jzhRStT5pd5D2z4k93VmQ+iFBQpAoULw/ONSjG5R03KFP5v8VBBNyfiE9OP+63Vc1zYqQoV4I8/oFcvo5f4o4/Ciy/CsWOZf40rV2DxYmO1Z4kSkAPVB7NKKEqp14FXMH5inqe1nvmg86WEkvcd2hdPbKOW1HM4wu9L1zPxnH5omSK3t/GyVNlGiEy5eBGmT4d584x54x06GE2wChc2Rhtpf0VGGrXzTZvgyBHj+cWLGwuFRo2CWrWyFYLFa+BKqZrAcqABcBfYBAzVWp/J6DmSwPO269dhdcVRDLr9Gbvfmc1gVTFLs0tyiyUXAQmRaeHhMHs2fPGF8eeMODpC48ZG0m7XDurUMbt0khM18GrAPq11dPIFfgO6AlPMeE1hJYmJMLvVasbf/oxrPV5lXMGqxNyTJHO6BWxmye7rwiqKFoX33zd29tm7FxISjHnjaX+5uEDDhlCoUK6EZE4CPwZ8rJQqCsQAzwL3Da+VUoOAQQDlypUz43IiJ/00/SIj/h7AVd8nKPH9Z1wav83keXlhOzJTPcdlForINR4e0Lq1taMAzEjgWuuTSqlPga3AHeAvIMHEeXOBuWCUULJ7PZFzIiLA6b1xuKk4vLYtA1dXi80uyarM1tNl93UhzJyForVeoLWup7VuCtwAMqx/i7zru6H7eC5uCTcGjEJVNHZHyonZJQ+TsqrT4j3ChbBTZiVwpVTx5N/LAd2AZZYISuSeoJNJPLFsBLcKlKTk52NTjz+o93dOmbo5KF1dG2R6oBAPYu5CnlXJNfB44FWt9c2HPUHkLet7L2UU+4mYtMjYoSSN3C5TZFRfzwt1dyHyIrMSuNa6iaUCEblv8+oonv9rLFfK+LOvcUumTt5h1ZuC1qq7C2GrZCVmPnX3Lpx+ZSplCOXE++8wbu1xq9eerVF3F8KWSQLPpxZ9eJGXbkzhUpNevBfulSdqz9aouwthy6SZVT4UHg7ek8fi6KAp9f2nXJpjur+DNWrPMj1QiMyTEXg+NPmVP+iZsJSvaz1H42Vn8S5geg8yqT0LkbfJCDyfWbL7Ek+sm0GEgycLmncg8lYMzg4KZ0eVrsOf1J6FyPtkBJ7PLBh7jG5Ja/i2Zhci3Yxpg/FJGg8XJ6k9C2FjZASej9y+DQP2LyZGufFds3bpHouIiefo+DZWikwIkR0yAs9Hlnx4jt5Jy/m+aiduFkjfLU3q3ULYHkng+UR0NBSY9SlJypGlbTqme0zq3ULYJkng+cTSKSH0il1EeKeBjOzXXOrdQtgBqYHnA7GxoKdOw4lESs4cQ5cKMtdaCHsgI/B8YNnn13ghei5X2/yfsVmrEMIuSAK3c/HxEPXJTNyIpeTn46wdjhDCgiSB27nVC2/xYuSXXGnSA1VVblQKYU+kBm7HtIZLE+fhxW0KTh9j7XCEEBYmI3A7tn1TPM9d/oJLVVvg4F/P2uEIISxMErgdCxizgrKE4DPpTWuHIoTIAVJCsSNpd3QvFF2UBYHTue5TlWKdnrF2aEKIHCAjcDtx747uJddfpD6H+afvAHCQj1kIeyT/s+1E2h3dEyLdGBb8HWFORRnnXcvKkQkhcopZCVwp9YZS6rhS6phSaplSys1SgYmsSbt7ziN7EunIRhbXaU9wVJIVoxJC5KRsJ3ClVGngf4C/1rom4Ag8b6nARNakdBNMinNi0LGfiFWuLH2yjXQZFMKOmVtCcQLclVJOQAHgkvkhiexI2dHd9aAHfZN+YOWjbYj2LkqLqj40nrwD37EbaTx5R67vNC+EyDnZnoWitQ5VSk0DLgAxwBat9ZZ7z1NKDQIGAZQrVy67lxMP0aVuaRLi4cTUr3Enll86dKV7/dKsOhSaWhsPvRXDuNWBqecLIWybOSWUwkBnwBcoBXgopf7v3vO01nO11v5aa38fH5/sRyoeKv54UYbEfcNV/2dZOn0AO0+FpSbvFDHxiUzdHGSlCIUQlmROCeVp4LzWOkxrHQ+sBhpZJiyRVVrDmfcWU5yw1IU7aW9sppXRcSGEbTEngV8AGiqlCiilFNAKOGmZsERWbf41iV6hnxFevh4OrVoAGW+TJjc2hbAP2U7gWuv9wErgMBCY/FpzLRSXyKI/xqynCqcp9NFoUAr478ZmWrJ9mhD2Q2mtc+1i/v7+OiAgINeul18cPAhxDZ6iVuEQCl37B5z+uzeddnl9KW93RretIjcwhbAxSqlDWmv/e49LLxQ7sHbsn3zMHmLHzEyXvMGYbSIJWwj7JEvpbdw//0D9HdOIcfPG7dWXrB2OECIXSQK3cd+/f4YurCHxlaFQsKC1wxFC5CJJ4Dbs2jUotWIGiQ7OFBz3mrXDEULkMkngNmzB5DBeTPyWqG59oWRJa4cjhMhlchPTRkVGgprzFe7E4v6B7LgjRH4kI3Ab9dXUKF6K/ZJbT3WAatWsHY4QwgpkBG6Dbt6EqGlf4cN1+HSctcMRQliJjMBt0KxJd/hf7BRuP9kGGkn7GSHyK0ngNub6dUj4fDY+XOdNv47S51uIfEwSuI35/KPb/O/uVHaVfYItXr5o/uvzLUlciPxFErgNuXIFmD2bYoQzs3nPdI9Jn28h8h9J4DZkxge3GZEwlR1lG3C01P0dBaXPtxD5iyRwGxESAq5zZ1GUGyzt2M/kOdLnW4j8RRK4jfhsfCQjEqcR3bI9HV7uLH2+hRAyD9wWHD8Onou+oAg34dMJqe1hpc+3EPmbJPA8LikJxgwM4wf9GXFtO+Lqb/R0lz7fQghJ4HncggXQ7cAYvBzv4DB9srXDEULkIZLA87CrV2H1yD/4lW/Rb46B6tWtHZIQIg/J9k1MpVQVpdTRNL8ilVIjLBhbvjd6RDxTo4YSX7Is6v33rB2OECKPyfYIXGsdBPgBKKUcgVBgjWXCElu2QLHls6jJMfhqDXh4WDskIUQeY6kSSivgrNb6Xwu9Xr4WEwMTXwlhixpPYrv2OHbubO2QhBB5kKXmgT8PLDP1gFJqkFIqQCkVEBYWZqHL2bePPoLXL4zEzTkBx9mzQClrhySEyIPMHoErpVyAToDJxtRa67nAXAB/f39t7vWsae2R0Byfe71pExyetJmP+Qne/wh8fS36+kII+2GJEsozwGGt9VULvFaetfZIKONWBxITnwj81wEQyHYSv/cbQs9KNZg20JW9zsNIKlcZh1GjLBa/EML+WKKE0psMyif2ZOrmoNTkncKcDoAp3xBCb8WggQshSbz1kgcLYp+nZMK/9HziZRrP2CMtYoUQGTJrBK6UKgC0BgZbJpy8K6NOf9ntAJj2G0JSvAPXVtfn46jxtEnazLi2wwkoUx0sMMoXQtgvs0bgWutorXVRrXWEpQLKqzLq9JfdDoApiV9rCN9Yh56X1zMqaTqL67ZnmV+71POkz7cQIiPSjTCTRretYtEOgKW83dEabv1WlRpBIcxTr/BnuVp80OqV+86VPt9CCFNkKX0mZaUDYGZmqwx9sipDBzviGZTIOqeOhHkU5tXOY0lwvP8jkT7fQghTJIFnQWY6AGZmtkpAALzXvxSFLl7iV682eMVEMqT/LJ5tXotVh0LT3SyVPt9CiIworXNvara/v78OCAjIteuZI7tzvhtP3kGoiZKHt7szBVycCNpRgps7q9G60D7WO/TAJfoWrFgB7dubdV0hhP1SSh3SWvvfe1xG4CaYM+fbVL1aa7gS5Enk/krEnCvOS8VnMTv8TeJLPILLtr1Qp07qudLnWwiRWZLATXjQnO+HJddS3u6pI/CkOCfuHCvNnSPliQ/3xNk1mtll+jEsZDF7ytfm474T+CVN8s6IjMqFEKZIAjchu3O+796FrmVrMu23MO5cLET06UfQ8U64lLhJlwbzeDtkJo+HnGBh/U583PIlkuIdH/h6kDMrQIUQ9kESuAlpR9H3Hk8RHw8nTsCCtTdZueUON4ILcveaFzqhOFAc54J3KVgllP5lf+LlwAXUPnCKMA9vRrZ/g9U1WwFQOhOzS8z5aUAIYd8kgZswum2VdKNeADcnR7qVrcnIkbB7NwQGGiNuKIxy9sSlRCSedf+lYNkIPu0EL8QFoL74HH49TlTpckx4ZjjLqrckzskFyPzsEkuvABVC2A9J4CaknfN94YLG4VwFok6XY9R5Z1xc4Kmn4PXXYe2F40R5hlHA8xaPXzpOs/OHafr7YaqtDjZeqEYN+OEHPHr1wi/wKluzUcfOzE8DQoj8SRJ4BhqVLs0jB0uzd7XxdePG8O4Y6NkTCnslwtGjFBn1JY32/8XjISdwT4gjztGJg2VqMKn5AMZ98QbUrJnayzu7s0tM/TQgc8OFECAJ3KQNG+CllyAiAt59FwYMgIoVgUOHYMQXsH493LzJWCCoWDmW1WnLHxX8+LNcbWJc3Cjt7c64WrUsEktWVoAKIfIXSeBpREXBm2/CN98YU7O3b4eaVRNg7Vp4cSbs2QMFC0L37tC6Nb/6VGPk79dyfHQsc8OFEKZIAk926BD07g3//AOjR8OHH2hcF8+D9h/DhQvGEHzGDGM4XqgQYOxkEecjc7SFENYhCRw4eRJatQIvL9ixA5r734EBr8Dy5cYdyy++gA4dwPH+edsyOhZCWEu+T+BXr0LzpxOITkzEs/0evvn5AvUGfoLXv2dh0iR46y1wkK67Qoi8J18n8OhoaPL0XcKuOVKiz0E6Xt7KlF8/J87JhT9mL+WpIb2sHaIQQmQo3ybwpCTo2xfOHHfGp8shxpz9ihF7lnG4VBWGdR6H4y0f9lg7SCGEeIB8m8DHjIHVq6Fwy5O8cHcZI/YsY2XNVoxrN5x4R2eUrHQUQuRx5m5q7A3MB2oCGhiotf7TAnHlqPnzYdo0ePVViNG7mPzNLP4sV4ux7V5L3RFHVjoKIfI6c0fgnwObtNbPKaVcgAIWiClHhYbCG2/A00/DzDGXia//AdcLFmZYmu3MZKWjEMIWZDuBK6W8gKZAfwCt9V3grmXCyjmjRhmdBL/5Ig6nXt1xirpN0MK1FDjvyC2Zyy2EsCHmjMArAmHAt0qpOsAh4HWtdZRFIssBO3caU7snjNdU/OxV+PNP+OknWj7XmpbWDk4IIbLInAnOTkA9YI7Wui4QBYy99ySl1CClVIBSKiAsLMyMy5knPt6oefv6wrjCX8OCBfDee/Dcc1aLSQghzGFOAg8BQrTW+5O/XomR0NPRWs/VWvtrrf19fHzMuJx5Pv/cWHH5zYTLuLz7FrRtCxMmWC0eIYQwV7YTuNb6CnBRKZVyt68VcMIiUVlYaChMnAgdO0LrnW9DXBx8+aWssBRC2DRzM9hrwBKl1N+AH/CJ2RHlgJQbl18NDIBFi1jyZFd85wfRePIO1h4JtXZ4QgiRLWZNI9RaHwX8LRNKzkh749J9wlCue3gzqd5zaGSDYCGEbbPrGoLWRi+qChVgnO9yiv4VwNQmfbnj+t909ZQNgoUQwtbY9VL6HTsgIAAWzIrG5d23OFaiEj/Vevq+82SDYCGELbLrEfinn8Ijj0Dfq1MhJISvOg8nyeH+nt6ybF4IYYvsNoEfPgxbt8J7/S/i/Nmn0KMHbYb0wN05fQKXZfNCCFtlNyWUtUfSb23msLMhXl4FePnf94zesVOm0KWCbBAshLAfdpHA1x4JZdzqwNTNhYPPKy5tdefl7qdw+WkJDBtm3MlEtkATQtgPuyihTN0clG5n+MgDFcEhiSYXPzBG3yNGWC84IYTIIXaRwNPOIkmMcuFOYBlKVj1FlyProHt3owGKEELYGbtI4GlnkUQe8oVEBwZ5foXX3Wh4800rRiaEEDnHLhL46LZVcHd2JCnOkTuHy+P5WAivnFzJdb8G8MQT1g5PCCFyhF3cxEy5KfnGu9EkxTkzuML3lD5zjX0vTqLz5B0y40QIYZfsYgQO0KlOaRyDHqPJU5qpN9dwp5wvA8MfIfRWTLq+J9K8SghhL+wmge/cCWfPwnvNf4eAAObU60x0gk53jvQ9EULYE7tJ4PPmQeHC0PLoZ1C0KAsrPmXyPOl7IoSwF3aRwK9fhzVrYFTHIBw3/AzDhlHEp7DJc6XviRDCXthFAl+8GO7ehSFxn4OrK7z6aurMlLSk74kQwp7Y/CwUrWHuXGjeIJoivy6Bnj2hRAm6lDAel74nQgh7ZfMJ/I8/ICgI5r68Cg5EwksvpT4mfU+EEPbM5kso8+aBlxc0DloIlSpB06bWDkkIIXKFTSfwmzfhp5/g9Q5ncfx9FwwcCEpZOSohhMgdZpVQlFLBwG0gEUjQWufqBsc//ACxsTDMbSE4OEC/frl5eSGEsCpL1MBbaK2vW+B1skRro3zyeL1EHtm0CNq1g9JS7xZC5B82W0I5cAACA2Fio81w6ZJRPhFCiHzE3ASugS1KqUNKqUGmTlBKDVJKBSilAsLCwsy83H/mzwcPD3j64kLw8YGOHS322kIIYQvMTeCNtdb1gGeAV5VS900B0VrP1Vr7a639fXx8zLycISYGVqyA/u3DcP7lZ+jbF1xcLPLaQghhK8xK4FrrS8m/XwPWAA0sEdTDvP/FDSIjocjpjyA+nu1Pts+NywohRJ6S7QSulPJQSnmm/BloAxyzVGAZWXsklDkLEnD0iOGFkLUcKVmF4X/dlTaxQoh8x5wReAngD6XUX8ABYKPWepNlwsrYJ6vOEXW2GE3K/0qV6xdYUbu1tIkVQuRL2Z5GqLU+B9SxYCyZcnZ/EUhyoH/Sd8Q4ubK+mlF2lzaxQoj8xuamEd49VZYCPmF0Ct7J5soNueNaAJA2sUKI/MemEnhQENwJ8aJrqZUUjr3N2urNAWkTK4TIn2yqG+H33xsr5j/w3sRNj0LsqVCX0tImVgiRT9lMAk9KgiVLoGPz21T8cysMGMCZqZ2tHZYQQliNzZRQ9uyB4GAY/dhaYyVPnz7WDkkIIazKZhL4998bS+efOLsEKlSARo2sHZIQQliVTSTw2NjkpfPPXMVp5zZj9C19v4UQ+ZxNJPANGyAiAoYXXwGJiVI+EUIIbCSBb9wIpUpBlYAlUKcO1Khh7ZCEEMLqbCKBL1gAf/5wFnVgv4y+hRAimU0kcAcHKPfHUuOL3r2tG4wQQuQRNpHA0dqYBN60KZQta+1ohBAiT7CNBH7kiLGO/oUXrB2JEELkGbaRwJcsAWdneO45a0cihBB5hm0spe/QwdhxvkgRa0cihBB5Rp4fga89Ekrj/Rrfa1VoPHmH7LwjhBDJ8vQIfO2RUMatDiQmPhGA0FsxjFsdCCDdB4UQ+V6eHoFP3RyUmrxTyPZpQghhyNMJPKNt0mT7NCGEsEACV0o5KqWOKKU2WCKgtDLaJk22TxNCCMuMwF8HTlrgde4zum0V3J0d0x2T7dOEEMJgVgJXSpUB2gPzLRNOel3qlmZSt1qU9nZHAaW93ZnUrZbcwBRCCMyfhTITeAvwND8U07rULS0JWwghTMj2CFwp1QG4prU+9JDzBimlApRSAWFhYdm9nBBCiHuYU0JpDHRSSgUDy4GWSqkf7j1Jaz1Xa+2vtfb38fEx43JCCCHSynYC11qP01qX0VpXAJ4Hdmit/89ikQkhhHigPD0PXAghRMYsspRea70L2GWJ1xJCCJE5SmudexdTKgz4N5tPLwZct2A4tkDec/4g7zl/MOc9l9da33cTMVcTuDmUUgFaa39rx5Gb5D3nD/Ke84eceM9SAxdCCBslCVwIIWyULSXwudYOwArkPecP8p7zB4u/Z5upgQshhEjPlkbgQggh0pAELoQQNsomErhSqp1SKkgp9Y9Saqy148kNSqlgpVSgUuqoUirA2vHkBKXUQqXUNaXUsTTHiiiltiqlziT/XtiaMVpaBu95glIqNPmzPqqUetaaMVqSUqqsUmqnUuqkUuq4Uur15ON2+zk/4D1b/HPO8zVwpZQjcBpoDYQAB4HeWusTVg0shyU3CfPXWtvtYgelVFPgDrBYa10z+dgU4IbWenLyN+vCWusx1ozTkjJ4zxOAO1rradaMLScopUoCJbXWh5VSnsAhoAvQHzv9nB/wnnti4c/ZFkbgDYB/tNbntNZ3MTofdrZyTMICtNa7gRv3HO4MfJf85+8w/uHbjQzes93SWl/WWh9O/vNtjN27SmPHn/MD3rPF2UICLw1cTPN1CDn0l5HHaGCLUuqQUmqQtYPJRSW01pfB+I8AFLdyPLlluFLq7+QSi92UE9JSSlUA6gL7ySef8z3vGSz8OdtCAlcmjuXtuo9lNNZa1wOeAV5N/tFb2Kc5QCXAD7gMfGbVaHKAUqogsAoYobWOtHY8ucHEe7b452wLCTwEKJvm6zLAJSvFkmu01peSf78GrMEoJeUHV5NriCm1xGtWjifHaa2vaq0TtdZJwDzs7LNWSjljJLIlWuvVyYft+nM29Z5z4nO2hQR+EHhMKeWrlHLB2DziZyvHlKOUUh7JNz9QSnkAbYBjD36W3fgZ6Jf8537AOivGkitSElmyrtjRZ62UUsAC4KTWenqah+z2c87oPefE55znZ6EAJE+3mQk4Agu11h9bN6KcpZSqiDHqBqNn+1J7fM9KqWVAc4w2m1eB8cBaYAVQDrgA9NBa281Nvwzec3OMH6s1EAwMTqkP2zql1FPA70AgkJR8+G2MmrBdfs4PeM+9sfDnbBMJXAghxP1soYQihBDCBEngQghhoySBCyGEjZIELoQQNkoSuBBC2ChJ4EIIYaMkgQshhI36f7rpWe2BgCEZAAAAAElFTkSuQmCC\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.851059\n",
       "x1                  0.519373\n",
       "np.sin(x1)          0.495606\n",
       "I((x1 - 5) ** 2)   -0.021622\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.797487\n",
       "1    10.637240\n",
       "2    10.364660\n",
       "3    10.024039\n",
       "4     9.673681\n",
       "5     9.371627\n",
       "6     9.161443\n",
       "7     9.061557\n",
       "8     9.060745\n",
       "9     9.120877\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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
