{
 "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.984\n",
      "Model:                            OLS   Adj. R-squared:                  0.983\n",
      "Method:                 Least Squares   F-statistic:                     935.7\n",
      "Date:                Fri, 10 Jul 2020   Prob (F-statistic):           3.22e-41\n",
      "Time:                        05:46:36   Log-Likelihood:                 1.9020\n",
      "No. Observations:                  50   AIC:                             4.196\n",
      "Df Residuals:                      46   BIC:                             11.84\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.0645      0.083     61.184      0.000       4.898       5.231\n",
      "x1             0.4762      0.013     37.303      0.000       0.451       0.502\n",
      "x2             0.4994      0.050      9.952      0.000       0.398       0.600\n",
      "x3            -0.0177      0.001    -15.775      0.000      -0.020      -0.015\n",
      "==============================================================================\n",
      "Omnibus:                        0.124   Durbin-Watson:                   2.203\n",
      "Prob(Omnibus):                  0.940   Jarque-Bera (JB):                0.321\n",
      "Skew:                           0.060   Prob(JB):                        0.852\n",
      "Kurtosis:                       2.627   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.6224677   5.0843021   5.50767542  5.86536907  6.13998753  6.32681633\n",
      "  6.43459668  6.48408926  6.50466334  6.52947147  6.59000261  6.71090859\n",
      "  6.9059539   7.17575419  7.50767513  7.87790822  8.25538243  8.6068691\n",
      "  8.90244128  9.12039072  9.25079438  9.29714482  9.27577636  9.21318142\n",
      "  9.14165785  9.09400251  9.09812284  9.17245174  9.32291911  9.54197756\n",
      "  9.80984244 10.0977424  10.37264644 10.60269082 10.76241417 10.83694071\n",
      " 10.82442433 10.73635188 10.59565592 10.43294695 10.2814842  10.17171127\n",
      " 10.12625425 10.15620441 10.25929587 10.42027789 10.61342056 10.80674211\n",
      " 10.96726342 11.06642646]\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.07678486 10.95641188 10.72489332 10.42686278 10.12107377  9.86601481\n",
      "  9.70558939  9.65836688  9.71303611  9.8311747 ]\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3zN9R/A8ddnZ7PNjA2bmMwld2OYS/wIC0VFSlHkGpKKENKFJKKbdEEUQiWFSCT3RMwl93vD5jLmOtvs9vn98d2Wbecw29nOZe/n47HH7Oz7Ped9dux9Pnt/P5/3R2mtEUII4XhcbB2AEEKInJEELoQQDkoSuBBCOChJ4EII4aAkgQshhINyzc8HK1mypC5fvnx+PqQQQji8HTt2XNRa+2W+PV8TePny5QkLC8vPhxRCCIenlDpp7nYpoQghhIOSBC6EEA5KErgQQjiofK2Bm5OYmEhERATx8fG2DkVkg4eHB2XLlsXNzc3WoQhR4Nk8gUdERODt7U358uVRStk6HHEbWmuio6OJiIigQoUKtg5HiALP5gk8Pj5ekreDUEpRokQJLly4YOtQhANYsiuSyasOc+ZKHGV8PBnetiod6wbYOiynYvMEDkjydiDyWonsWLIrklE/7yUuMRmAyCtxjPp5L4AkcSuSi5hCCKubvOpwevJOE5eYzORVh20UkXOSBA6YTCaCg4OpWbMmderU4aOPPiIlJeW254SHh7NgwYJ8ilAIx3LmStxd3S5yxi5KKHcjL+pqnp6e7N69G4CoqCieeeYZrl69ytixYy2ek5bAn3nmmVw9thDOqIyPJ5FmknUZH08bROO8HGoEnlZXi7wSh+a/utqSXZFWewx/f39mzJjBZ599htaa8PBwmjVrRr169ahXrx5//fUXACNHjmTTpk0EBwfz8ccfWzxOiIJoeNuqeLqZMtzm6WZieNuqNorIOTnUCPx2dTVrXhipWLEiKSkpREVF4e/vz+rVq/Hw8ODo0aN07dqVsLAwJk6cyAcffMDy5csBiI2NNXucEAVR2u9j5r+WAZpOXCszU6zEoRJ4ftbV0vYKTUxMZNCgQezevRuTycSRI0fMHp/d44QoKDrWDciQnGVmivU5VALPr7raiRMnMJlM+Pv7M3bsWEqVKsU///xDSkoKHh4eZs/5+OOPs3WcEAVVfv0FXZA4VA08P+pqFy5cYMCAAQwaNAilFFevXqV06dK4uLjw7bffkpxs/Af09vbm+vXr6edZOk4IYZCZKdbnUCNwS3W13L57x8XFERwcTGJiIq6urnTv3p1XX30VgIEDB/LEE0/w448/0rJlS7y8vACoXbs2rq6u1KlTh549e1o8TghhkJkp1qfSar0WD1Dqa+ARIEprXSv1ts7AGKA60FBrna2rdSEhITrzhb2DBw9SvXr1u49c2Iy8ZhnJkvHsSauB33vmOKWvR1M4IQ7f5Jt0qeFL7WImiIuD9u2hcWNbh2p3lFI7tNYhmW/Pzgh8NvAZMPeW2/YBnYDpVolOCAd1NxfmCnqi71gklkabPqb0ht8zfmP5Lf9+913o0gUmToTAwHyNzxHdMYFrrTcqpcpnuu0gSF8MISxdmBvzy/4MybplNT9+2hFZMGdgXL1qJOYpUyjt7g7vvQctWoC3NxQp8t/nxESYNAkmT4YlS2DoUBg50vieMCvPL2IqpfoppcKUUmHSxU44kiW7Imk6cS0VRv5K04lrzS4Ys3QB7kpcYoYFZ/O3nip4vUG0hlmzoHJl+PBD6N4djh7ldLdRxAXfD7VqQfnyUKIEuLsbifqdd+DIEXjiCRg/3jj3669BJgWYlecJXGs9Q2sdorUO8fPLsqmyEHYpu6t+s3sBztKVJqedgZGSAoMHQ9++ULUqhIVxcNgseo26h4oVYc6c25x7770wbx5s3Wok+D59oHNnSErKr+gdhkNNIxQiv2S3m565qa1aQ9y/JbnwSzBn5zbhwpJ6XF5bnWth5Yk9Woqk6+7pxzrlDIzEROjRAz79FAYPZuv7G+j4Tj1q1IAffoAXXoCHH87G/TRqBH/9ZZRUFi+GAQOMH65I51DTCIXIL9mds3zr1NZTpzTqSAXO7ShN4hVPXDwSKFTqGgkXvIk77o9OSk30pmSKhoRTqtm/ztcbJC4Onn4ali3j2vBxdNoxmjWfKHx94a23YNAguKs/xJWCYcP+q6P7+xs1dAFkI4Erpb4DWgAllVIRwNvAJWAq4Af8qpTarbVum5eB5pXo6GhCQ0MBOHfuHCaTibRSz7Zt2yhUqJAtwxM2cjdzltvXCmDVjAC2pM7Jqt0onqTi63nm+lyqXD7FFb8iXPcqRtlygfx9zpM/wpvwz98PwOHynK1oIikIXJ1hKHX1Kjz2GGzaRPjwz2kybyDXrxvl7379cnkt8p134MIFmDDBeAcYMsRqYTuy7MxC6WrhW4utHItNlChRIr2V7JgxYyhSpAjDhg1L/35SUhKuTvHbJe7G8LZVM0wPBPOrfmNijPLsypXw0ksw4uE9BPw8leRvv8V08ybnihSnaEIshROMTbs7AyjFpTZPM/jy2wwcWI2pU+Gzz6BVq3x8gtZ29arxBPbsYfPA+YR+2pXSpWHVKggKssL9KwWffw7R0fDqq1CypHFRtICTGrgZPXv25NVXX6Vly5aMGDGCMWPG8MEHH6R/v1atWoSHhwMwb948GjZsSHBwMP3795cl9E6iY90AJnQKIsDHEwUE+HgyoVNQhil/588bs+FWr4YlI7bw6d6WBLSrA/PnY+rRA/bu5Z7r0RS+GQfx8XDmDOzZAyNHUnzzMubsqMmO+k+jzu4l9EFNpccPWbU1cr65eRM6dkTv2cPcxxfzv8+70qgRbNtmpeSdxmQyLm6GhkKvXvDrr1a8c8dkV0PLwYMhdTBsNcHB8Mknd3/ekSNH+OOPPzCZTIwZM8bsMQcPHuSHH35g8+bNuLm5MXDgQObPn89zzz2Xu6CFXcjcTe9Whw8bF+LOn4ctr/9Cg0lPGX/aT5pkzJooXjzjCe7uULq08REUBEOGcGzYW1Rf8A3/JC9ihncvBi/5gj7nI0iZGkmn+g4yNzwlBZ57DtavZ0r9bxny4yP07WsMlvOk+ujublzQbNnS+NNn9Wpo2jQPHsgx2FUCtyedO3fGZDLd9pg1a9awY8cOGjRoABg9Vfz9/fMjPGFD27YZydtkgr3D5lBxfB8ICTFGhCVKZO9O/PzoWb0zCf2b8tJf3zNw1yzqev1Dxy3L6dcL2m0Du29oqbVRzli4kK+rTWLIjm58/DG88opR8cgz3t7w229G4u7c2firpmTJPHxA+2VXCTwnI+W8cmszKldX1wx7ZMbHG/VMrTU9evRgwoQJ+R6fsI2zZ43rdMWKwbYuH1HynaHQujX8/PNdX6U7cyUOXcSXt9q8wLZ7azJ5xRR2ugfz2N7lPPhgKZYuzf77gU188AFMmcIvlQbT59AwvvzSmOmXL/z8YOFCY6ph377GqLwArgyXGng2lC9fnp07dwKwc+dO/v33XwBCQ0NZtGgRUVFRAFy6dImTJ0/aLE6RtxIT4cH2N7kQncQQ726UnDCUyNaPwLJlOZpiceuMluXVm/Nkt0lo92Q2m5pS6e8FNGtmlGjs0rx58NprbC77NB2Pf8inn6r8S95pgoONWSlLl8LMmfn84PZBEng2PPHEE1y6dIng4GC+/PJLqlSpAkCNGjV49913adOmDbVr16Z169acPXvWxtGKvNK5z3UO7HLno8DneWnPfOYHP0SbBgNYcuBiju4v8yKg/aUq8VSfT4mpE8ycpGfpdfR1WrbQ9pfEf/8d3asX+/1b0ipiDh986MJLL9kolsGD4cEHjc+HnbgtgSVa63z7qF+/vs7swIEDWW4T9q0gvmY//aQ1aP10helag54X/JAOfG2ZDhyxXDeZsCbH97t4Z4RuMmGNLp96P4t3RmidkKB1v35ag/7Ktb+uVT1JnztnxSeTG2FhOqVIEX3St7YuyhU9caKtA9JaR0ZqXaKE1vXqaX3zpq2jyRNAmDaTU+2qBi6EPTpyBHr2hMp+/zD9zFD2lqrE2ND+6TXX3PQzsTjTZdo0KFGCvhMmUPTINdq0nMPv69woVSrHD5V7x4+j27UjmhI0vvwbw8cVY8QI41s2bZVbpoxRQnn8cWO558SJ+fO4dkASuBC3ceOG0RjPyy2Bn127onQKL3YYSYKrW/oxedLPRCljyXixYjw1ciTeh6/xcIsfWbnBE5tMdIqKQj/0EDeuJtH05ipMrW4yK+ZXfptoJ61yO3Y0lntOmgRt2xrTDAsAqYELcRtDh8L+/bD1gRHUOnuQNx4dwinf0unft/aerJCpja1uwO7XJ/CQXsEnRx6m3f+uEZnfa31iYtDt25MYHknrm8u52MQNl5AD9tcq96OPjPaz3bvDpUv5+9g2IglcCAs2bIDp0+HrRxcTuPgTePllWo1+4barM3PLXBvbri7B7Hj3U5qpP/nqeEsev/8cJ05YPv9OPczvSlwc+skn0Tt28kTSD5xsGkCR/x3MMGPPblrlennBggXG1J3+/QtE50IpoQhhRlwcPP88NL/3X3ps6AUNGsCkSXR0d8/TsoClNraveNVg89Kl1H7yKRZFNqZ7oxV8sb4GNWv+d9zdbO+WLdeuoR97DL1hI32ZSeCLj7LX61fI5nRrm7TKrV8fxo2DUaNg/nzo1i3/Y8hHksCFMGPcODh2NIW/a/dEXcNoZO3ufsfzcuu2bWzbt8e0aQNlHnqEpdFN6dLoZ870KMJ17yjK+HgSm5BksZRx1/tzXrhAStuHSNm9h+7MZ3lwE4LKrMU3yY3LsYlZ4lNkHInnRWkp24YPh+XL4cUXoXlzKFfONnHkAymhABEREXTo0IHKlStTqVIlXnnlFRISEgBYv349jzzySJZzli9fTt26dalTpw41atRg+vS839+5Z8+eLFq0CIC+ffty4MABi8euX7+ev/76K/3radOmMXfuXIvHi//8849xLWxW02/w3bPRWHFYoUK+PLalUWv67SEhuIZtRQeUZOmNhwiafpDY8BJEXokzm1gh65vCHXcbOn2apCbNSNh9gMf0Un5rVo/ibfZx5locMfFJuJkyDsE93Uw827hcnpaW7orJBHPnGn1aevY0PjupAp/AtdZ06tSJjh07cvToUY4cOUJMTAyjR4+2eE5iYiL9+vVj2bJl/PPPP+zatYsWLVrk6PFz2r1w5syZ1KhRw+L3MyfwAQMGSJOtbEhKMnpRVSseRY/9w40RXO/e+fb45nb4yTKaLV+err0/YFeZqsxP7sFLPyzjxrry6GTztY3Mbwq33W3o8GESGjYl9vhZ2rr8xvYOpfBpciy95p2YovEq5JolWb/bMYjNI1vx78T2bB7ZyvYbNVesaPTmWLcOpkyxbSx5qMAn8LVr1+Lh4UGvXr0AMJlMfPzxx3z99dfExsaaPef69eskJSVRIrVRhbu7O1WrZv1zccyYMXTv3p1WrVpRuXJlvvrqK8BIri1btuSZZ54hKCiI5ORkhg8fToMGDahdu3b6aF5rzaBBg6hRowbt27dPX7IP0KJFC8LCwgBYuXIl9erVo06dOoSGhhIeHs60adP4+OOPCQ4OZtOmTRla4u7evZvGjRtTu3ZtHn/8cS5fvpx+nyNGjKBhw4ZUqVKFTZs2WeNH7FCmTIEdO+DXqkNwib1hXMV0yb9fE0ttbIEMFycP3XTlua7vsLBGa97gPTZs60zFWS4kXsy4pN9cKcNcmcaUksyDaxaSWL8RV8/H80Tx9Rx91hWvaueyHHs1LtG+krUlvXsbjWtGjTKmEjkh+6qB26Cf7P79+6lfv36G24oWLUq5cuU4duyY2XOKFy/OY489RmBgIKGhoTzyyCN07doVFzO/6Hv27GHr1q3cuHGDunXr0r59e8DY7Wffvn1UqFCBGTNmUKxYMbZv387Nmzdp2rQpbdq0YdeuXRw+fJi9e/dy/vx5atSoQe9Mo8ELFy7w/PPPs3HjRipUqMClS5coXrw4AwYMyLA5xZo1a9LPee6555g6dSoPPPAAb731FmPHjuWT1J9RUlIS27ZtY8WKFYwdO5Y//vgjGz9k53DiBLz5JrzdeCWBfy6At9+GatXyPY7Mi3vMXZxUQIKrG689+gorajTlnRXTWXf5YaZ/3Y8vHurL1aAoAnw905N304lr0+vdPoUz1rEbnN7HmN++oubl46zmQaZUm8bXv1eiy/y1RF7JGp/D7OOpFHz1FdSqZVzM/PvvPOpxazsFfgSutUaZ6WJm6fY0M2fOZM2aNTRs2JAPPvggS2JN06FDBzw9PSlZsiQtW7Zk27ZtADRs2JAKqXXV33//nblz5xIcHEyjRo2Ijo7m6NGjbNy4ka5du2IymShTpgytzGzZsnXrVpo3b55+X8Uz96HO5OrVq1y5coUHHngAgB49erBx48b073fq1AmA+vXrp29aURBoDQMHQlHTDd6IfMHYSX3UKFuHBZgveWj+mwyyvlIIbftPZVbdzvTVM1n522MMmneY3uf8OHfaJUu9O62OXfraBT748VN+XDCSIpcTebbwj+x+/3cW7qjEvfdms5xj7/z9jSS+ezdY6OvvyOxrBG6DfrI1a9bkp59+ynDbtWvXOH36NJUqVSI6OtriuUFBQQQFBdG9e3cqVKjA7NmzsxyT+U0g7etb29VqrZk6dSpt22bcVnTFihW3fRNJO/dOx9wN99SZFiaTiaSkJKvdr7376Sdj+6+wVmNxXRtuTALPh1kn2WFpZorGKLGcuRJHcX9fSsz6GJfkEbh2foFh4UNhCpyf4o+PRzM2Bwbj4pFI1UvhVL96nOqxR7knKYqbFGK8aRTHer/EtA9L4+393/3fumGzTZbIW0uHDkY55f33oX17p9oAosCPwENDQ4mNjU2foZGcnMzQoUPp2bMnhQsXNntOTEwM69evT/969+7dBAYGmj126dKlxMfHEx0dzfr169M3f7hV27Zt+fLLL0lMNP6sPXLkCDdu3KB58+Z8//33JCcnc/bsWdatW5fl3Pvvv58NGzakt7i9lLoCzdvbm+vXr2c5vlixYvj6+qbXt7/99tv00XhBdf26Ub17qspu6m34yLiK2by5rcNKZ6lkEeDjmaUWrULq4/fvNjh5kouTvmaDf0MeSPiTqYffZMo/79Dr9EK8r8WzSrXm03sn8+WgA7x48T2+mVGaNceyLgLqWDfAMerdd/LJJxAYCM8+a2yO7CSysyv918AjQJTWulbqbcWBH4DyQDjwlNb6ct6FmXeUUixevJiBAwcybtw4UlJSaNeuHe+99176MWvWrKFs2bLpX3/33XdMmjSJ/v374+npiZeXl9nRNxilkvbt23Pq1CnefPNNypQpw5EjRzIc07dvX8LDw6lXrx5aa/z8/FiyZAmPP/44a9euJSgoiCpVqphNtH5+fsyYMYNOnTqRkpKCv78/q1ev5tFHH+XJJ59k6dKlTJ06NcM5c+bMYcCAAcTGxlKxYkW++eabXPwEHd/YsXA2MoWZtV5AFS9uzCG0I9ndYDmDcuUoObwXU5IDee1yLIEnL5Kc7Mbp0r7gqSnrayT/NFZfBGRvvL2NDSCaNYNOneCPP+zmL6zcUPoOy02VUs2BGGDuLQl8EnBJaz1RKTUS8NVaj7jTg4WEhOi0mRNpDh48SPXq1XMav10zt8u9M3Cm12zPHqhXD2Y1/ZoeG/vA7NnQo4etw8oip93+MidmMJJ/5nnaTSeuJdJMqSZtlO80fvgBunQxNkWeNcthdvFRSu3QWodkvv2OI3Ct9UalVPlMN3cAWqT+ew6wHrhjAhfCnqSkwAsvQIVil+i+f4RRG+3e3dZhmXW7DZbvdB7cuY592xWgzuTpp40phePGQc2aRrcyB5bTi5iltNZnAbTWZ5VSFhtcKqX6Af0AyjnxklZzLO1mL+zDnDnw119woOWbqA2X6FG/Bxtf/81xL9ZZkJ3kX8bH0+wI3GGmDN6NMWPgwAFjyX21asaFTQeV5xcxtdYztNYhWusQPz8/S8fkdRjCShzttbLUnS862vj97V1nB9XWf8m39R9lg2cZ80vLCwCnmDKYXS4uxrt3cDB07erQi3xymsDPK6VKA6R+jrrD8RZ5eHgQHR3tcImhINJaEx0djYeHh61DyZbb9fwYPhyuXk7h05QXueTlw+Smz2Q41yb9rG3I0gpQZ/krJAsvL/jlF+Pzo4/CxZzta2prOS2h/AL0ACamfl6a0wDKli1LREQEF5xoao8z8/DwyDAjx55Z6vkxeuoFDnwTwKJ23+C14m/eeGQo1929spzvdPXfO8hpnd1hlS0LS5bAAw8YSXzJEvJsz7qUlDxpyZCdaYTfYVywLKmUigDexkjcC5VSfYBTQOecBuDm5pa+ilAIazKXgFNumjj8UxUaVb5Ep79HQLNmbGvyMFyNz3KsU9Z/RUaNGhmbQDz7rFFS+f57I6Fbi9awaJHRo2HpUmOFrxXd8S1Ba91Va11aa+2mtS6rtZ6ltY7WWodqrSunfi4Y+xcJh2IuAV/eUI3k654srTYCdeUKfP45wx+qVnDqvyKrTp2MPilFi0KrVsZepNZoQbtlizGz6amnjB4sZhbW5VaBX4kpnFfmC3Pxp4sTs6s8o1uuoNSymTBsGAQFFbz6r8iqdm0ICzOS7ejRxsyUnNbFT5wwpis2aQL//gszZ8KuXRCSZRp3rt1xIY81mVvII0ReSlsAE3HhJlFzHuAe9ySOu9fFxdVk7NzgKWUScQutjRbCr7xiNMKaO9coqdypfp2SYqwK+/Zb+OwzcHU1pjkNGwZFitz+3GzI8UIeIRxZ2oW54cPhg2jY0HkYLj+egPXrJXmLrJSCAQOgYUPo3Nkoqfj4QOPGcP/9xkejRsbS/MOHYe1a42PdOrh0yTi/Vy9joVCZMnkfrozAhbPbutUoRb7bYTujljbm38efoVtIL8fusCfy3rVrRpvKLVuMj/37jRG6UuDrayRsMPbcDA01kn2rVnmSuC2NwCWBC6cWFWVsVO7hksAh7xASoqJo3uNzLpj+m8turjeIEFlcvQrbthkjgvBwYzTeqpWxX2oe91SREooocBITjb+Co6PhWM/3MX25l7efG5cheYPlnduFyKBYMWjd2viwEzILRTitoUNh40ZYOOYAZWa9C08/zY+l65o9tqAt2hHOQRK4cEqzZ8PUqTBq0HUemdPZmOP76acWF+fIoh3hiCSBC6ezfbsxkSC0lebds73h0CGjD7S/f8Fq2iScntTAhVM5d85YWHfPPbD0f5NxeWcRTJ5sXGzCifZ5FAJJ4CKXzp2DlSuNWVa+vkZ/oIAA43NgoLEWAnK+o8zdOHTIWEAXHQ17PvoDrxdHGSvrMjXtL3BNm4TTkgQugLtLsH//DcuWwYoVxgphMNY63LhhzPy4VaNGUKv5FTYkHibR3bhQmBf7La5bZ4y8CxWCzQtOcl/fLlC9ukNtmyXE3ZIELrK9oW18vLHCeMYMMJmMabDvvQcPPwx16hhrHC5ehIgI42P/fqP0PGuyD6iWeARexKtGJIWrnSUO603dmz0bnn8eqlSBFT/FEfhsJ+OdZPFiqyxjFsJeyUIeka0NbU+ehCefNPr9jBgBI0cao24SE9k4fwU7f/qdk3gSH3AvHTo24aHW9YwsDwT02UDMwTLEHggg6WphlHsiRWpF4B18isiZOWvduWRXJJN+O8yB5fdydUtl6jSKZ8MP8RTr8ySsWWP8ifDIIzn+mQhhT2Qhj7DoThvarloFzzwDSUlGz/sONY7Ctyvhjz9IXLOW5jdiaH7ridMhxdUVlwoVoEMHQotUYWPzGHyaHeFmRHGu7yrH9V2BXN9RgRbHoF8/aNcu9Q0hG+ZvPMOr469xKawhSZeKUKTOKYrXX098s3cofOYkI9sPYeu+wgwPiJRat3BqksCFxQ1tSxfz5J13jD1ga9WCn75PpPI3r0PHD4wDKlbk15otWFW6FmFla+J98wZlr54n4NoFqt+M5jmva/DJJ8xNSuKwf3l+qtGSX6o/wLnHLuGWcJj7k+vz57JiPPusMVhv1AjatIG2bY3Om66p/ztTUoyKyNGj8MUXMH2WPykJZShU5jIlH91JQ99NfDNnLB5JCTzXeSxbAutAHtTZhbA3UkIRWWrgYMyNbnzjfmZ/XIxu3WDGm6fx7Pm0Md1kwACjVWbFilQY+SuW/gcpoLpbAhOSDnLvrz9RfO9OkpULW6o1wmXgQJoMfIYUXNi8GX7/3Rjph4UZtfRChYxrj4mJGXvru7uDW+XTFKl3EvfSV2l1bBuf/fI+lzyL0bPzGI6VLJchhlvLQEI4KimhCIvMzY1uVzKItwYUpXitKExxX5AQ/BGuLhq37783mtWnsjR6B9DAgcRCdHGrx4Q5PejoHYdp9mz+N3MmvNSdG++/weygtsy67wE8y9zD69Or0qxcAGvWGItxXFzAzQ2OX7zGX/9eIIZYKjW8QqJbHP4nj9Jj5TKe3rOa/aUq0ueJt7lQxDdLDLJEXjgzGYGLLM6cgZq1k4nTN/iw4vO8GLaQA/4VePXJ0Qzo+1CGkoS50bs5GUbCCQls/2gWfPkFDU7tI8HFlR1lq7OlUn3q9n2all3apDfQv/X+TSnJPHjsb3rtWEbjU3uJc3VnUVAo77XoTXwhD7N/CcgIXDgDGYGLbElMNNa+XLsG79d6hRfDFrKgTlvGhvbjppt7lql/mUfvloYDGUbChQoxmKpEdp1IlQvhdNq/jmb/7uLVdbNh3WwY7Gc08DaZKHngDNPib+KWnETg5bMEXL9ARFE/PnmwD6sat+dQYiHK+HjSspofP+2IzFIGkiXywpnlagSulHoFeB6j3PmV1vqT2x0vI3D7N2QIfPIJPNPgE+ZvH8KiWqEMazc4w2KYAB9Piwt+sjMlETBbO/eLucz/Tu7mY98oo4ZiMrHvYjyJLq4kmFy55lGERbVC+aNyI1JcTPw7sX2G8/NjtacQtmD1EbhSqhZG8m4IJAArlVK/aq2P5jxMYUsLFxrJe1z3Iwxe+AZ7S1VidJuBGZK3gvQEbW7Bz/C2Vc1eEM08EjZXO79QxJdtTdvBLYm+/23eEDKTJfKioMlNN8LqwOGBWD8AABa1SURBVFatdazWOgnYADxunbBEfjt5Evr0gQcbXWf09o64ebozuPOb3HRzTz9GQZZRc9pmCGmyu8O7pa6ALav50XTiWiqM/JWmE9fSspqfdA8UwoLc1MD3AeOVUiWAOKAdkKU+opTqB/QDKFeuXOZvCzvx2muQnKRZWrwXavth3Fev5iXfqhlKEpZmm2Se6ZGdkbC5mS+Z69iRV+L4aUckT9QPYN2hC1IaESKT3NbA+wAvAjHAASBOaz3E0vFSA7dPmzZB8+awOnQiD64ZBR9+CK++muW47Na3cyqv718IR2WpBp6rDR201rO01vW01s2BS4DUvx1McrLRoKq9/3ZC174OXboYVzLNyOvNEO60pF8IkVGuphEqpfy11lFKqXJAJ+B+64Ql8svs2bB7Vwq/3/cSSvnD9OkW26/m9WYIlso0st2ZEObldh74T6k18ETgRa31ZSvEJPLJtWvw+uvwzn3fUvLY30Y2L1r0tufk5UyP7M5gEUIYcpXAtdbNrBWIyH/jx0Nc1DVGJI2Axo2he3ebxiPbnQlxd2QlZgF1/Lgx53tJzXG4HYiCT5elL1+3JZnLLUT22f43VtjEsGFQ03SIhw5/Ar17Q4MGtg5JCHGXZAReAG3ZAkuWaI5VHow6X9jYF00I4XBkBF4AvfRaPB0KLabS0VVMadaNJZGJdz5JCGF3ZARewHz0XRR7/vRhoccQjpa4l6nV2+AmO9cI4ZBkBF7AjH8PernOomL8Kca37EOSyTVLPxMhhGOQBF6AHDgA1/b5MNo0nl2lq7K+Yv3078lqRyEcjyTwAmTCBOhhmkO5m5FMado1w4pLWe0ohOORGngBcfw4/DA/iVNe77HXu3KG0besdhTCMckIvICYOBG6mxZwT0w4N157nQDfwrft1y2EsH+yqXEBcPo0VK6YzCmv6vhX8IKdOy02rBJC2B/Z1LgAmzwZOqf8gP/Vo/DWz5K8hXASksCdXFQUzJqRzDHvcVAuCDp0sHVIQggrkQTu5D7/HB65uYjSNw/BmwvtomGVEMI6JIE7sbg4+PLzFLZ7j4N7a8ATT9g6JCGEFUkCd2Jz50Kj6F8JZD+8Pk9G30I4GUngTiolBT76CL4vMgldIhD19NO2DkkIYWWSwJ3U8uVQ/MgW6vInjJ8CrvJSC+Fs5LfaSX34IbztORntWRzVp4+twxFC5IFcFUWVUkOUUvuVUvuUUt8ppTysFZjIubAwOLvxCG3jl6AGDgQvL1uHJITIAzlO4EqpAOBlIERrXQswAV2sFZjIuQ8/hFFuH0KhQjBokK3DEULkkdxOS3AFPJVSrkBh4EzuQxK5cfIkbFh4nm4pc1A9e0KpUrYOSQiRR3KcwLXWkcAHwCngLHBVa/175uOUUv2UUmFKqbALFy7kPFKRLVOmwCA9FdeUBBg61NbhCCHyUG5KKL5AB6ACUAbwUkp1y3yc1nqG1jpEax3i5+eX80jFHV29Ct99FcPLbl+gHn8cKldmya5Imk5cS4WRv9J04lqW7Iq0dZhCCCvJzSyUB4F/tdYXAJRSPwNNgHnWCEzc3pJdkUxedZgzV+Io4+PJ8LZVObI6gKdiZlGEyzB8OEt2RTLq573EJSYDEHkljlGy/6UQTiM3CfwU0FgpVRiIA0IB6RWbD8wl5pE/7uP8NF/+cf2QbffUZMj6WGIT9qcfkyZt/0tJ4EI4vtzUwP8GFgE7gb2p9zXDSnGJ25i86nCWxHzxn3tod3kx5ZJOM71RJyKvxHE5NtHs+bL/pRDOIVcLebTWbwNvWykWkU2ZE7DWcG1rBUaZnuGQbyBrKzW47fmy/6UQzkG6GzmgzAk47oQ/rS9tolbyAaY3egKtLL+ssv+lEM5DErgDGt62Kp5upvSvr/1dkVGm94jw9mNZ9eYZjvXxdCPAx1P2vxTCCUkvFAeUloAnrzrMvwcLUf/0If7HX7zbZABJpv9eUk83E2MeqykJWwgnJQncQXWsG0DHugE89RT0cn2MlKIlqPPWqwRsPJ1haqEkbyGclyRwB3biBBxatI+H9TJ4ZSyPNqnMo00q2zosIUQ+kRq4A/voI3iNSaQU9oIXX7R1OEKIfCYJ3EGdPw+rZ56kC9/h0r8flChh65CEEPlMEriDev99eDHhI0wmYMgQW4cjhLABqYE7oDNn4OcvznHY5SvUs8/CvffaOiQhhA3ICNwBTZwIwxLeoxAJ8MYbtg5HCGEjksAdTEQE/DbtJP3VdFTv3nDffbYOSQhhI5LAHcz48fB60jvGJvNvvmnrcIQQNiQJ3IGEh8OmmYd5jjnGZsVS+xaiQJME7kDGj4e3Ut5GeXrAqFG2DkcIYWMyC8XOmNtpp2PdAI4fh51f7+arlB9gyGjw97d1qEIIG5MEbkdutwXa0k8DeIc3SSnmg8uwYbYMUwhhJySB2xFzO+3EJSbz9teReM85SXu9HEa8Bz4+NopQCGFPJIHbEXNbnekUOPRjZda6dSKlmD8uL79sg8iEEPZILmLaEXNbnV3fFciT51fQNGE9Lm+/BV5eNohMCGGPJIHbkcw77SRdd8dtQwk+c30F3agRDBhgw+iEEPYmxwlcKVVVKbX7lo9rSqnB1gyuoOlYN4AJnYLSt0CL3xTER8mv4aOuoWbNwuhcJYQQhhzXwLXWh4FgAKWUCYgEFlsprgIrbaedlSvhk/dX8izz4I23oGZNW4cmhLAz1rqIGQoc11qftNL9FWixsTBsQAy/u/Yn5b7quLz+uq1DEkLYIWsl8C7Ad+a+oZTqB/QDKFeunJUezrmNHw/PnxxNaXUaNetPcHe3dUhCCDuU6wSulCoEPAaYXduttZ4BzAAICQnRuX08R2BpNWV2jvOMCmD9xK1sYqrR76RJExs8AyGEI7DGCPxhYKfW+rwV7svh3W415a1J3Nxxr04/Qew8Lza79YaSATBhQv4/ASGEw7DGNMKuWCifFESWVlNOXnX4tsclxxbi3HdBzI/vRMXko7jMmQ3e3vkRshDCQeVqBK6UKgy0BvpbJxzHk7kMEmlmNSVkXWV569c6yYWLP9Xlq+svEqrXwezZEBqal2ELIZxArkbgWutYrXUJrfVVawXkSNLKIJFX4tAYZRBl4djMqyzTvtYaolcG8faZD+imFzC9TW/o0SNvAxdCOAVZiZkL5solGrIkcU83E8PbVs1w2/C2VfFwNXF1c2V67F/KSN7n+3rtKDVhbN4GLYRwGtLMKhfMNZ8CI4kH+HhmmF0C0HTi2vTbXmpejdK7m1N98yqm8hKbqjfBc9oXdKhXNsv9ZXdWixCiYJEEnguWat4BPp5sHtkq/evMM07Cj5no86EHIy9OYLR6D9WwEc3WrobChbPcV3ZntQghCh4poeRC5uZTYL5ccmup5caBMnjPKcn6S214k3dxea476vdVZpN35nPTmJvVIoQoeCSB50Lm5lMBPp5M6BSUZWR85kocCReKcHF5EJ2XbWFncgPuczvCCx1HGTNOiha1+BiWyjSWbhdCFBxSQsmltOZT5iQmwpIlkPhDLVqHr6cPb9CcP1kfWJ/h7V+hUNk7l0AslWnM9Q4XQhQsksCt7MwZ+Ptv2L4xjqg5v/HQ5QUcYzke3ORU0VK80Xgg84IfxrOQKxMylVrMGd62aoYaOJgv0wghCh5J4DmlNZeOX+bIH6c4s+Uk13Ydx+XEMUrfOEowx3iMU5hIId6nFIW69WNDyIO8fsaLM1fjCbiLmSRpx8gsFCFEZkrr/OsvFRISosPCwvLt8azq7Fmu/rCS83NW4nF8H8VjTlFEx2Q45IabD9dLV8ZUrTK+DSrj2uJ/0KIFuMr7pBAi55RSO7TWIZlvl8xyO/v2kThnATE//obvyd0UA25Qmn1ejUmq3AaPKuUoWa8cgc0DKRZcAa8SJZAdK4UQ+UUSOFkXyoxsVpaHf5iGy2dTUFqxh6ZsKTaBwk88TOuhtWlTw9KCeSGEyD8FPoFnXihTY/s66r87A9cbUUyjP7ufHM/TA0vw2gPgIpMuhRB2pMAn8LSFMvdcu8jY1dNoe2wrewji6cILeePnBxjQ1tYRCiGEeQU+gZ+5EkfA1SgWffsaxWJvMJxJTL+vCz4PH6StJG8hhB0r8Am8pimOT797E48bSTRx3cTpB93wrb2Hsr6yUEYIYd8KXAK/9YJlZY9kvlnwJr7XLtHatIrzXePwLnNOFsoIIRxCgUrgt16wdE+8ybj5b+EXcZzHWEp0Fx88ypyy2P5VFs8IIexNgUrgaRcsXZOT+GLpRBpEHKAr33Ho4cqEz7sPCAKkhasQwjEUqIlxaR38xv4xjdDj23mBL1lR+34Iuv2GwyAtXIUQ9qdAJfAyPp7cf3IPz+5eyWSGMjewE8Xb7CMg0wVLaeEqhHAEud2V3geYCdTC2Emst9Z6izUCywsjWgQSNOk5jqsKvOMzAr+OOyns4ZLlgqWlFq7FPN2kLi6EsBu5HYFPAVZqrasBdYCDuQ8p7zy2Yg4VLkcyQE+nSLuj3HuPm9kNGMzttOPmoriRkJRhB/pRP+9lya7IfHwGQgjxnxyPwJVSRYHmQE8ArXUCkGCdsPLAgQPoiRNZ4NKNMt1as3qO5UPNtXCNTUjicmxihuPS6uIyChdC2EKO28kqpYKBGcABjNH3DuAVrfWNTMf1A/oBlCtXrv7JkydzFXCOpKSgmzfn+raD1Cl0iC3H/Ljnnru7iwojf8XcT0oB/05sb40ohRDCLEvtZHNTQnEF6gFfaq3rAjeAkZkP0lrP0FqHaK1D/Pz8cvFwufDVV6jNm3k58UNeHnf3yRssb2EmW5sJIWwlNwk8AojQWv+d+vUijIRuX86eRY8YwV8eLQmr0YNBg3J2N9ndgV4IIfJLjmvgWutzSqnTSqmqWuvDQChGOcW+jBhBckw8PZOnMf0zhZtbzu5GtjYTQtib3K7EfAmYr5QqBJwAeuU+JCs6dgw9fz5T1RDqPV2Fli1zd3e324FeCCHyW66mEWqtd6fWt2trrTtqrS9bKzBrCB/+Jje1G5PUYI5V2CxT/oQQTsVpV2KuWrGNgF9+5Cvdl9j68VzUV2TethDCqThtAo8Z9x5auzDZZSjeIeGA9DMRQjgX50zgZ8/SfvsK5vAcl2uZcC1yM/1b0s9ECOEsnLOd7Icf4paczERGUrThiQzfknnbQghn4TQJPG2nnbgz59g87XMWu3YhqlIRSpT4b2GozNsWQjgTpyihpG3AEHkljl5hv+CeeJN3k0bzZK8YAnw8UUCAj6fZxlVCCOGonGIEnrYBQ9H4GHrsWMbPrh05cU8pfPU/bB7ZytbhCSFEnnCKBJ52YbL7zl8pmhDLu7xF0cYn5IKlEMKpOUUJpYyPJ67JSXTf9Sur3EI5ULISnhWj5IKlEMKpOUUCH962Ku3Cw7gn5hKfJ75M0YbHKVxILlgKIZybU5RQOtYNoEnERiJcA1jpHkr9+7cxor1csBRCODeHTOBpUwbTugKOqe5G660beYt3eH20F2PeyGXXKiGEcAAOl8DTpgzGJSYDxt6Upyd+TpJyZabuyxb76ocohBB5xuESeNqUwTTuiTd5fPdqlpoepc6DpQkMtGFwQgiRjxzuImbmqYHtDm/GN/46nycNom9fGwUlhBA24HAJPPPUwG67VnDErRJ/ejXl0UdtFJQQQtiAwyXwW/emrB51gvpnDvFl4kDaP5FAoUI2Dk4IIfKRw9XAb92bstuqFcS7uDMnpSdbR3vbODIhhMhfDpfAIXVvyvuKosf/yc9eXQiqW5wqVWwdlRBC5K9cJXClVDhwHUgGkrTWIdYIKlvmzUPFxDCJFxgkFy+FEAWQNUbgLbXWF61wP3dn1ixO+tThcEpDnnwy3x9dCCFszuEuYgKwdy/s2MGU673p1l3hKT2rhBAFUG4TuAZ+V0rtUEr1s0ZA2TJnDskmN+YmP8Pzz+fbowohhF3JbQmlqdb6jFLKH1itlDqktd546wGpib0fQLly5XL5cEBiIsybx4YijxBYqSR16uT+LoUQwhHlagSutT6T+jkKWAw0NHPMDK11iNY6xM/PLzcPZ1i1Cs6f55OrPenePfd3J4QQjirHCVwp5aWU8k77N9AG2GetwCyaPZuYwn787vIwXbvm+aMJIYTdyk0JpRSwWCmVdj8LtNYrrRKVJdHR6F9+4XuPQbRs40apUnn6aEIIYddynMC11ieA/K1Af/cdKjGRTxN7MqJbvj6yEELYHcdaiTl7NqdK1OVEfG06drR1MEIIYVuOMw88de735zd68vjj4OVl64CEEMK2HGcEPmcOKa5uzIp/hgUy+0QIIRwkgafO/d7m9whuuiStWtk6ICGEsD27L6Es2RXJ8L6T4Px5Jpx/jgahMbg6xtuOEELkKbtO4GkbGDfesYYLhXxZkdKevR57WLIr0tahCSGEzdn1WDZtA+NRD71M0QslUMnxpJS4zORV8ekbOwghREFl1yPwtA2Mb8QUZUdUM7xqRqJU1o2NhRCiILLrBJ62gfGNA8Zo26vGmQy3CyFEQWbXCTxtA2PXIvF4BZ3GtVgcnm4mhretauvQhBDC5uy6Bp6+gbHXYc5ciaCMjyfD21aV+rcQQmDnCRxSNzCWhC2EEFnYdQlFCCGEZZLAhRDCQUkCF0IIByUJXAghHJQkcCGEcFCSwIUQwkFJAhdCCAeltNb592BKXQBO5vD0ksBFK4bjCOQ5FwzynAuG3DznQK21X+Yb8zWB54ZSKkxrHWLrOPKTPOeCQZ5zwZAXz1lKKEII4aAkgQshhINypAQ+w9YB2IA854JBnnPBYPXn7DA1cCGEEBk50ghcCCHELSSBCyGEg3KIBK6UekgpdVgpdUwpNdLW8eQHpVS4UmqvUmq3UirM1vHkBaXU10qpKKXUvltuK66UWq2UOpr62deWMVqbhec8RikVmfpa71ZKtbNljNaklLpXKbVOKXVQKbVfKfVK6u1O+zrf5jlb/XW2+xq4UsoEHAFaAxHAdqCr1vqATQPLY0qpcCBEa+20ix2UUs2BGGCu1rpW6m2TgEta64mpb9a+WusRtozTmiw85zFAjNb6A1vGlheUUqWB0lrrnUopb2AH0BHoiZO+zrd5zk9h5dfZEUbgDYFjWusTWusE4Hugg41jElagtd4IXMp0cwdgTuq/52D8x3caFp6z09Jan9Va70z993XgIBCAE7/Ot3nOVucICTwAOH3L1xHk0Q/Dzmjgd6XUDqVUP1sHk49Kaa3PgvGLAPjbOJ78MkgptSe1xOI05YRbKaXKA3WBvykgr3Om5wxWfp0dIYErM7fZd93HOppqresBDwMvpv7pLZzTl0AlIBg4C3xo23CsTylVBPgJGKy1vmbrePKDmeds9dfZERJ4BHDvLV+XBc7YKJZ8o7U+k/o5CliMUUoqCM6n1hDTaolRNo4nz2mtz2utk7XWKcBXONlrrZRyw0hk87XWP6fe7NSvs7nnnBevsyMk8O1AZaVUBaVUIaAL8IuNY8pTSimv1IsfKKW8gDbAvtuf5TR+AXqk/rsHsNSGseSLtESW6nGc6LVWSilgFnBQa/3RLd9y2tfZ0nPOi9fZ7mehAKROt/kEMAFfa63H2zikPKWUqogx6gZwBRY443NWSn0HtMBos3keeBtYAiwEygGngM5aa6e56GfhObfA+LNaA+FA/7T6sKNTSv0P2ATsBVJSb34doybslK/zbZ5zV6z8OjtEAhdCCJGVI5RQhBBCmCEJXAghHJQkcCGEcFCSwIUQwkFJAhdCCAclCVwIIRyUJHAhhHBQ/wfyVeoFkACOMQAAAABJRU5ErkJggg==\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.064507\n",
       "x1                  0.476216\n",
       "np.sin(x1)          0.499431\n",
       "I((x1 - 5) ** 2)   -0.017682\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.076785\n",
       "1    10.956412\n",
       "2    10.724893\n",
       "3    10.426863\n",
       "4    10.121074\n",
       "5     9.866015\n",
       "6     9.705589\n",
       "7     9.658367\n",
       "8     9.713036\n",
       "9     9.831175\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.4rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
