{
 "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.980\n",
      "Model:                            OLS   Adj. R-squared:                  0.979\n",
      "Method:                 Least Squares   F-statistic:                     754.1\n",
      "Date:                Fri, 21 Feb 2020   Prob (F-statistic):           4.21e-39\n",
      "Time:                        13:57:24   Log-Likelihood:                -4.1308\n",
      "No. Observations:                  50   AIC:                             16.26\n",
      "Df Residuals:                      46   BIC:                             23.91\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.0831      0.093     54.429      0.000       4.895       5.271\n",
      "x1             0.4907      0.014     34.071      0.000       0.462       0.520\n",
      "x2             0.5216      0.057      9.212      0.000       0.408       0.636\n",
      "x3            -0.0189      0.001    -14.954      0.000      -0.021      -0.016\n",
      "==============================================================================\n",
      "Omnibus:                        2.536   Durbin-Watson:                   2.416\n",
      "Prob(Omnibus):                  0.281   Jarque-Bera (JB):                1.892\n",
      "Skew:                          -0.473   Prob(JB):                        0.388\n",
      "Kurtosis:                       3.122   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.61036341  5.09172228  5.53276677  5.90507171  6.19047043  6.38403949\n",
      "  6.49490757  6.54475568  6.56425502  6.58802778  6.64895881  6.77279276\n",
      "  6.97390428  7.25293632  7.59669448  7.98031483  8.37134897  8.73509521\n",
      "  9.04029999  9.2642927   9.39671042  9.44120056  9.41482152  9.3452401\n",
      "  9.26618593  9.21190973  9.21155609  9.28437501  9.43655927  9.6602269\n",
      "  9.93471621 10.22998043 10.5115243  10.74607157 10.90703217 10.97887089\n",
      " 10.95965979 10.86139491 10.70802541 10.53151874 10.36660852 10.24508788\n",
      " 10.19058594 10.21468562 10.31502061 10.47566411 10.6697456  10.86386558\n",
      " 11.02358297 11.11907377]\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.11980579 10.98378856 10.73147613 10.40948069 10.07916026  9.80159615\n",
      "  9.62263807  9.56167858  9.60690506  9.71819191]\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deViV1fbA8e9mEBBFHMAJB1IzNUGUHMvZLPWmaWampWVXzUyzsjTrqj9LLU3TspwyczYnHG/mXNpVQ8GccsgRnElEBAQO+/fHARQ8RxEOnGl9nodHeHmH9Xp0nX32u/faSmuNEEII++Ni7QCEEELkjiRwIYSwU5LAhRDCTkkCF0IIOyUJXAgh7JRbQV6sVKlSunLlygV5SSGEsHv79u27prX2y769QBN45cqVCQ8PL8hLCiGE3VNKnTW1XbpQhBDCTkkCF0IIOyUJXAgh7FSB9oGbkpKSQlRUFElJSdYOReSAp6cnAQEBuLu7WzsUIZye1RN4VFQURYsWpXLlyiilrB2OuA+tNTExMURFRREYGGjtcIRwelZP4ElJSZK87YRSipIlS3L16lVrh2JTwiKimbDxGBdiEynn68XQttXpFFLe2mEJJ2D1BA5I8rYj8lplFRYRzfCVB0lMMQAQHZvI8JUHASSJi3xnEwlcCHs1YeOxzOSdITHFwISNx5w+gcsnk/wno1AAV1dX6tSpQ61atQgODmbSpEmkpaXd95gzZ86waNGiAopQ2KoLsYkPtd1ZZHwyiY5NRHPnk0lYRLS1Q3ModpfAwyKiaTJ+K4HD1tNk/FaL/IPw8vIiMjKSw4cPs2nTJjZs2MDo0aPve4wkcAFQztfrobY7i/t9MhGWY1cJvCDe1f39/Zk5cybffPMNWmvOnDnDU089Rd26dalbty6///47AMOGDeO3336jTp06TJ482ex+wrENbVsdL3fXLNu83F0Z2rZ6jo7PjwaJLZBPJgXDrvrAC6q/8ZFHHiEtLY0rV67g7+/Ppk2b8PT05MSJE3Tv3p3w8HDGjx/PxIkTWbduHQAJCQkm9xOOLePfXW76eh35AWg5Xy+iTSTrYl7uNBm/VfrFLcSuEnhBvqtnrBWakpLCwIEDiYyMxNXVlePHj5vcP6f7CcfTKaR8rpKQIz8AHdq2epY3JwB3F8Wt5FRiE1MAx3rDsha7SuDm3tUt3d946tQpXF1d8ff3Z/To0ZQuXZoDBw6QlpaGp6enyWMmT56co/2EyODI3QymPpkkJKdyPSGFMnHXKHY7HgClNavmnqeTegJcXKBGDZBZvjlmVwnc1Lv6w/Q35sTVq1fp378/AwcORCnFjRs3CAgIwMXFhR9//BGDwXjtokWLcvPmzczjzO0nhDkF1SCxluyfTBoPmMvHv86jy+Ft9+48Nf3Pxx+H77+H+vULJkg798AErpSaA3QArmitH0/f1hUYBdQA6mutC6SzNy/9jfeTmJhInTp1SElJwc3NjVdeeYV3330XgAEDBtClSxeWLVtGixYt8Pb2BiAoKAg3NzeCg4Pp3bu32f2EyJB9XHSLx/xYsS86XxskNiEuDsaPZ9usL0FrpjfowoEy1dDpk8JKeHswtksQ/PMPjBwJjRrB4MEwZgzI/6P7Uhl9vWZ3UKopEA/MuyuB1wDSgBnA+zlN4KGhoTr7g72jR49So0aNXIQurEVes4eX/YElGJN1l3rl2fbXVcd8qKc1zJwJn3wCV69yvl1nej/yHH97l8rcxcvdlXGda9+557g4GDYMvvsOAgONx7dubaUbsB1KqX1a69Ds2x/YAtda/6qUqpxt29H0k1oqPiEcT0oK3LwJxYubfWC57a+r7BrW0koB5qO0NBg0CKZNg6ZNYcMGKoSG8vaDZmf6+MC338JLL8Ebb0CbNvDaazBpEvj6Wu9+bFS+94ErpfoCfQEqVqyY35cTwvrOn4cZM2DWLLhyBby8WOBZnItFS3HRpxTh5WuyrHZrUl3dHOKB5T1SUqB3b1i0CN5/H774AtIbe80Dy3NUl6dbX3jkkfuco2lTOHAA/u//YMIEOHwYtm8HL8d4PmAp+T6RR2s9U2sdqrUO9fO7Z01OIRyD1rB5M3TuDJUrw7hxpNRryOUPviTl3wM4FfAohQwpND5zgHEbv+G/P7xNs1P7HOaBZaaEBHj+eWPyHjcuM3lfumTsGalUCT76CNavz8G5vLyM51i+HPbuhddfN/49izu01g/8AioDh0xs3w6E5uQcWmvq1aunszty5Mg924Rtk9csm1u3tG7XTmvQaSVL6jPdh+nBHU9rDw+tjRlH66K+qdqjTKz2qnZBv1z7a33at6zWoC8+2VLro0etfQeWERur9VNPaa2U1tOna621/vtvrfv319rDQ2sXF627ddM6IiIX5x471vgXOWaMZWO2E0C4NpWbTW28ZydJ4OIu8prdJTZW6yef1GkuLnpTu0m6RmCiBq2LFdP6rbe0nj/fmHv69dM6pHGi9vKLT0/o1/Wadv+n03x8tHZz0/qjj7Q2GKx9N7kXG6t13brGe1myRGut9bx5WhcqZPzq21frEyfycP60NK179jSmrOXLLROzHcl1AgcWAxeBFCAK6AM8n/79beAysPFB59GSwB2GvGbpLl/WOiREp7m76zdL/aRB62bNjEl7ye9RuvG4Lbryh+t043Fb9Kr9UZmHbd6sdVCQ8X9fu9DL+kr73sYfunXTOinJeveTW7dva92ypTF5r1+vDQatP/zQeEstWmgdFfXgU+RIYqLWjRpp7eWl9b59FjqpfTCXwHMyCqW7mV+teuj+GhsUExNDq1atALh06RKurq5k9NXv3buXQoUKWTM8YavOn4fWrTGcPc9LHmv41eUZdu+GBg0eXOOkVSvYvx9++AE+/tgf//A5LG9Qky5LP4DLlyEsDIoVs+bd5VxamrFveutWmDeP+Kbt6NkZVq+G/v1h6lQLTqz09IRVq4yTfJ57Dv74A8qWtdDJ7ZNdVSPMDyVLliQyMpLIyEj69+/PkCFDMn/OSN5a6wfWBxdO5MQJePJJUqIu8bT+hQNln+F//zMmb8hZKVVXV+MouRMn4L33FC/sGcr4WvPRO3fCU09BdNaqhDZbtfCjj2DhQhg7lrNNX6FJE1i7Fr7+2jga0OKz4kuXhjVrIDYWOnaERAccxfMQnD6Bm3Py5Ekef/xx+vfvT926dTl//jy+d41DXbJkCW+88QYAly9fpnPnzoSGhlK/fn12795trbBFfrt1C557jsR/EmiYuJ1bIU+ya1fWIXEPU+OkaFGYOBHmzoVPjvWkf4UNpJ06bZyNeOQIYMOLI0ybBp9/Dm++ydGOw2jQAM6ehf/+FwYOzBw5aHnBwcY3jfBweOutfLqIfbCpWijvvAORkZY9Z5068NVXuTv2yJEj/PDDD0yfPp3U1FSz+w0aNIgPPviAhg0bcubMGTp06MChQ4dyGbGwaYMGoY8do4PeRLkOISxdCoULZ90lNzVOevWCcuWgS5c2nCn8K+sS2+HerBns2MGEjZdsqmphWEQ0eyZ9z2cLRrGrRmMOtR7BiGZpxN9Oxe/F3YzebyDBL59nlHbsaGz9f/YZdO0Kzz6bf9eyYTaVwG1NlSpVeOKJJx643+bNmzl27M7H4+vXr5OYmIiXTDpwLIsXw5w5jGUEJbu2YtEicDPxPyi3RdfatIGdO6FduxDqx/7Kbs9meLRsiee/RkPJgHv2t/QkoJysYRkWEc3SKUv5YclnRJZ7lN6NRnCudwm0IRX/7rtxLxVPdCwFUyb2k0+MfeL9+sGhQ8ZZnE7GphJ4blvK+eXuglQuLi6ZNcIBkpKSMr/XWssDT0f399+k9e3HH25NWFptFL/PMZ28IW9F14KCYPduePbZatQ/toW93s1Z8tMIurw0nnPFsz6ws+QkoJwuLrF0wWa+XTqKi0VL0vvpzzi7vDk6xZXS3XdTqFR85n4F8gnBwwPmzIHGjY2zhL79Nv+uZaNsKoHbMhcXF4oXL86JEyeoUqUKq1atyhyt0rp1a6ZNm8aQIUMAiIyMpE6dOtYMV1hScjJp3V7iVqIrvQstZOUKN4oUuf8huV3kASAgwDhrvFWrGjQ5uplf3ZuzZMkIur48nuhi/oCxRd/iMb97VreB3L1xmHvwOmrN4czzPe6ayLfff0iqiys9243naNgzpN12p/RLuynkf/OecxZImYAGDYx9r5MmwYsvQvPm+X9NGyIPMR/C559/zjPPPEOrVq0ICLjzkXbatGns2rWLoKAgatasyaxZs6wYpbC4ESNw2RdOL8P3jJxTiYIoxFiypHFmfnL12rRM2Uzx1Nv8tHQEZeOuUd7Xiy71yrNiX3SWB5tDlx1g6PIDuXrYaS7ZxiamEB2biFdyImO/H0bJhFheaz+Gvf/tjCHBg9Jd9+JRJs7ksQVWJmDMGKhSxTisJyGhYK5pIx5YTtaSpJysY3Cq12zjRnjmGaYxgKNvTeObbwr28levQsuWUPzEXra6tsatjB9s3UqTxX+bfFBqSnlfrwdWPGwyfqvZ87kZUpm9YgxPnong9Q6jWfK/QaTeKIx/1714VriOu4sCBSmGO7nknjKx+W37dmjRAt59F778smCuWYDMlZOVFrgQ5ty6RUqffhx1qcnS0IlWyQt+frBlC8RUqU8rwyZSrvwDTz1FoVMnc3yOnHRlDG1bHS9313u2K53GuJ+/ofnpfQxv+TZL9rxNaqw3NXsdxKvCdcr7ejGhazATXgimvK8XCuMbRoEmbzB2nfTvb3yQ5kTDeKUPXAgzjg38kOrRZ+lXaCu3Wkbw3yOVrDJsz9/fONGxRYsGND21lR1xT7NsyXC6dx3DCb9KDzw+J10Zph68Jicm8dHySTx/ZDuTGvZk6pFRpMQUocYrBzk0KwQIMXkOq/n8c2OZw9dfh4gI40NOBycJXAgTtqzYRtMfZ/ADvTnSpjhFVDTDV8YRfvaffF9Bx9xwvu3boW3bEJ44tIPfCrdk6ZLhvPLiGA6XrgJgtisjp0u0ZXnwmpjIxWc6UvbIdsY3eY1RZyaSfNmH8i9EMHZQGYver8X4+BhX8Hn2WWO/+KefWjuifCddKEJkpzXF3n2Pm7oo/6nwAd61jA8BE1MMLNx9Ll9nRN5v1qW/P2zbBt5P1OSJ+J3g7sPiJSMIjTpi2a6MuDh49lnK/raZvUPGMu7Cl9y+4MujLx3mm2FlrN/Svp9nnoFXXzW2xv/809rR5DtpgQuRjV6wkNBzEfR3+Za0dhdxvWtKePZH/pYe73y/OiqdQsrj6wu//AIdO1YlZOtODvi1ZvnSj6DGKAhqBq6ueYvl2jVjCzYyksuTFtJ7ZncSzsOihdC9e+283VxBmTTJOJ//3/+G3383Fp5xUNICF4I7xaKC3llKTN/B7KE+Pz3VFHffBz8AtOR455zUUSlSxNjVG/yvSgRe3cvBx16Ajz82Psg7cyb3F//5Z2jYEA4d4sjYMB7/rDuXLsGmTdDdXE1SW1SypPFh5t69xnotDkwSOBAVFUXHjh2pVq0aVapUYfDgwSQnJwOwfft2OnTocM8x69atIyQkhODgYGrWrMmMGTPyPc65c+cycOBAAKZPn868efPM7nvmzBkWLVqU+XN4eDiDBg3K9xjt0d3dFu9uW0jxpFgGl/iCEo3OZdnPXG0mS453Nneu7Ns9PWHFCnju1eIEHVrEmOrzSTvwp7HQ04IFD1x67O7qhi++P5+LzZ42trxdXPhl6CZCPm5P8eLGAR12OTeme3djd8pHH8G5cw/e3045fQLXWtO5c2c6derEiRMnOH78OPHx8YwYMcLsMSkpKfTt25e1a9dy4MABIiIiaJ7Lf+UGg+HBO5nQv39/Xn31VbO/z57AQ0NDmTp1aq6u5egyui0ev3SSV/9cxzQGcPZfLhTxdM3Sn9yjYcV7hto9zEPCnDA1nM/cNdzdjVUM589XTLjQkzoc4ErZIHjlFWNJ2q++MpYHzCbjDSv2cgwfbJ/L/K9ex2f3TiIHfMzHnQ7RdsyTNGoE//sfPPqoxW6tYCkF331nfCN7803HXUvT1CoP+fVliyvybN68WT/11FNZtt24cUOXKFFC37p1S2/btk23b98+y+9jYmK0n5+fTkhIuO+5R44cqXv27KlbtGihq1atqmfOnKm11nrbtm26efPmunv37rpGjRpaa63nz5+vn3jiCR0cHKz79u2rU1NTtdZaz5kzR1erVk03bdpUv/HGG/qtt97KPPeECRO01lqfOHFCt2rVSgcFBemQkBB98uRJ3aBBA+3j46ODg4P1pEmTstxHTEyM7tixo65du7Zu0KCBPnDgQOY5X3vtNd2sWTMdGBiop0yZYvK+rP2aWVrlD9fpyh+s0eEla+qLlNYBIZG60ofrdOUP192z76r95lfasZTcXOPUKa0bN9bahVQ9v+4knVrj8TsLctata1xLcupUrQcN0jurN9AnSwTo2y5uWoP+qVZr/diT27R7kSQNWr/+unGRHYcwebLx72DxYmtHkifkdkWeAmWFerKHDx+mXr16Wbb5+PhQsWJFTp40PVmiRIkSPPfcc1SqVIlWrVrRoUMHunfvjovLvR9o/vzzT3bv3s2tW7cICQmhffv2gHG1n0OHDhEYGMjRo0dZunQpu3btwt3dnQEDBrBw4ULatGnDyJEj2bdvH8WKFaNFixaEhITcc40ePXowbNgwnn/+eZKSkkhLS2P8+PFMnDiRdevWAcauoAwjR44kJCSEsLAwtm7dyquvvkpk+t/7X3/9xbZt27h58ybVq1fnzTffxN3iVfltSzlfL+r/up56MUd43Ws6qvnFzO3Z5aXGSU7l5hqBgbBjB4wd60rv/xvCG25DeK31Sfr5r6L2yVW4fvKJccciRfAt7M+xUhX5pUoDVrl0YuuRLhgOe+FZ8RrbN3rQuHE+3JS1vP02LFoEgwfD009DiRLWjsiibCuBW4HWGmWi8ry57Rlmz57NwYMH2bx5MxMnTmTTpk3MnTv3nv06duyIl5cXXl5etGjRgr179+Lr60v9+vUJDAwEYMuWLezbty+zdG1iYiL+/v7s2bOH5s2bZxbN6tatG8ePH89y/ps3bxIdHc3zzz8PgKen5wPveefOnaxYsQKAli1bEhMTw40bNwBo3749Hh4eeHh44O/vz+XLl7PUfXFEw58szxMj57GH+qx/tg5eha5ZvGukILi5wX/+YyyVPWcOLFtWlembh+LlNZSX/3WFsmU0J+P82fBHLAnX3DGc8UAnu+NR7jol2x+ganACjRvff8q93XF1hVmzoF49eP9941+MA7GtBG6FerK1atXKTGYZ4uLiOH/+PFWqVCEmJsbssbVr16Z27dq88sorBAYGmkzg2d8EMn6+u1St1ppevXoxbty4LPuGhYXd900k49iHZeqYjOt43DV7zdXV9b4LWTiKpot/wCf5Gj0rz6dwtWv5NkGnoAQHw5QpxtF0O3fC0qWwYoU/CQlQpgyU9ylMlMd1qHwNr0eu4PnIVQoXcmVoWzsZJviwgoNh6FAYPx569jQWl3EQD3yIqZSao5S6opQ6dNe2EkqpTUqpE+l/Fs/fMPNPq1atSEhIyBzRYTAYeO+99+jduzeFsy+1ki4+Pj5Ll0RkZCSVKpme0rx69WqSkpKIiYlh+/btJheIaNWqFcuXL+fKlSsA/PPPP5w9e5YGDRqwfft2YmJiSElJYdmyZfcc6+PjQ0BAAGFhYQDcvn2bhIQEihYtys2b95b4BGjatCkLFy4EjF0rpUqVwscJi+EDcPIkXtMns8DlVWZufYbT49uza1hLu03ed3N1hWbNjGWyL12CmzeNa3AeifBg3kIDtbueonCVqwQUt0LtkoL2n/9A1arQp49xPU0HkZNRKHOBZ7JtGwZs0VpXA7ak/2yXlFKsWrWKZcuWUa1aNR599FE8PT0ZO3Zs5j5btmwhICAg8ysiIoIvvviC6tWrU6dOHUaOHGmy9Q1Qv3592rdvT8OGDfnkk08oV67cPfvUrFmTTz/9lKeffpqgoCDatGnDxYsXKVu2LKNGjaJRo0a0bt2aunXrmrzG/PnzmTp1KkFBQTRu3JhLly4RFBSEm5sbwcHBTJ48Ocv+o0aNIjw8nKCgIIYNG8aPP/6Y+79AO3fllXdJSitEzNDxpPdoOaTsH+Q6hZRn17CWDvWGdV9eXjBvHkRFwWuvOc6oFFNPNrN/AZWBQ3f9fAwom/59WeBYTs5ji6NQ8tPdI0UciaO8ZrfX/qw16AmlxuukJGtHIwrEpEnGUSmTJlk7koeCmVEouR0HXlprfTH9DeAi4G9uR6VUX6VUuFIq/OrVq7m8nBAWlpJC3GvvcIKq1Jn7jjMUrhNgHOn2/PPwwQfGafZ2Lt8fYmqtZwIzwbigQ35fz5aMGjXK2iFY3PWEZC7dSKLdsPV2/bAvZtRUSl37i9lN1jKsvWRvp6GUcSRKvXrQrZux7GypUtaOKtdym8AvK6XKaq0vKqXKAlfyEoR+wJA9YRuuJyQT9U8CKWlpWSrlgQ3UgjbDVGnWjsVT8f78P/zXtT09F7e3doiioPn6wrJl0KiRcdbq+vVgYg6HPcht1GuAXunf9wJW5zYAT09PYmJicjUcThSsS7GJJN+6wdnYlMxtGZXybJHJ0qwr/uRw234YDHDlP9MIqCANB6dUt65xrOXPPxuHF9qpB7bAlVKLgeZAKaVUFDASGA/8pJTqA5wDuuY2gICAAKKiopD+cdt3/noCZ2NT+HrP9Szbc1KNz9wiBfnJVGnWp/bt5vHjG/mm8kQGfPzg1WyEA+vXD377DT75xNgab9HC2hE9tAcmcK21uUKSrSwRgLu7e+aMRGHb3jCz8O2DqvFltIQzkmlBdb1kf2MpevsWozfNIII6tFk32F4/NT80a7x52gWlYMYM2L8funaF1auhSRNrR/VQnOSfsLCEh6mUd7f7LVKQn7K/sQxevYLShisMf2IEr6/9lcBh62kyfqtFV9SxNfdb4UdgLK6+dq2xRkrLlsa6KXZEErjIsU4h5RnXufZDL9mVk0UK8sPdbzi1T//N66eX8Z3XvznZ0ttpEpq13jztStWqxtq5jRpBjx4wapRlJ/rs2gVt2sD585Y7ZzrbqoUibF5uKuWV8/XKVddLXmXEOWnDYT5bM4MLlOP77l1J5XaW/Sy9LJotsdabp90pWdK4Vl2/fjB6NBw/bhxumIPicGYdOAAjRhhHuZQpAydPQoUKlosZaYGLApDbrhdL6BRSntl/7SEo6Qg7unzNdb/bJvdz1ISW0xV+BFCokDFpjxsHixdDq1aQm8EVf/9tbMmHhBhb3+PGGZN3PjwklQQu8l1uu14s4eKCLVRZMIoNxbrTZX4np0to1nzztEtKwbBhxnHi+/dD7dowYACsWwe3bpk/LjbW2ILv3x8eewxWrYIPP4RTp4znu6v6qEXDLcjx16GhoTo8PLzAriesw1ZGPSSejCaxRghX0krhvn8vVYKL3DMiBowJzZGr8dnK62F3wsPh009h82Zj8vbwMC4Q2q6dsXV98CDs2WP8Opb+TMHNDf79b+PQxLJlLRaKUmqf1jr0nu2SwIUl2UqC1MkpHA9oQfmrkez77g+a9a+RJUZJaCLHbt82jhffsMH4deyuB8D+/tCgwZ2vJ56AYsUsHoIkcFEgmpgZK17e14tdw3JXSD83CXd/83epu2MyK7supvNPL+XqukKYdPIkHDkCQUFQqdK9tXrzgbkELqNQhEVZetRDbiYBHR69nLo7JrM+cCCdlkjyFhZWtarxywZIAhe5p7WxJbJlC+zbB+XK8cbfBsI9S/N3yQDiPItk7prTh4TZW9sJyalmxzGbSuAXNh+h4ujXifRswFN7v3Sa2ZbCOUkCFw8nNRXmzzc+cd+6FdKXgbvuWYaiyTGMSLtT6OpS4ZKsq/kkK0Pb07dbuwee2lRr2xxTLfrjc3/H7/V/kYgXRdb/hE+pQg97d0LYFUngIucuXkS/9BLq11+54V2W7S5PE0ZLttKSSx6l0TqNSrejeYy/qM4xGibs5tXwNbwRvpqkyOYwuJ+xmL6Z1RNMzRo0J3uLPnzEKmqNfZlLbgEkr/6Z6i0r5vVuhbB5ksBFzuzYgX6xG8kxN+nDfFak9qBFS0VA7Vjc4g9QtqhxzevbKS4cSSpG/eBeXLz9Pl3mX6bWHz/Qd/dMHtndncQifqi3B+L53lvG2W93yWk/efZxzL92m8aTP73N4cL18du9lsDafpa7byFsmPQQivvTGiZMQLdqxdnYYoQa9vDo6J7E/KPYsAEOF9+PoWh85u4u7mkYit7k538O8vbbsHZvafqdGsaSMSfpW/FntsTXx3PcSG6XqcjV7oPg9OnMY831k/t6ud8zCQigyaebmVy2P01/GshvJdoReHorZSR5CyciwwiFeUlJ8PLLsGoVYW4vMLjI98xY7MMzz9zZJXDYekz9C1LA6fFZV7vR2jg3YvXYwzy6diIvGRbiioFz9V+g5MCX2V61Fh9uOPnAMeSLNxznwJD59D75E4+lHWd2qR5MfqM7n71YR8ZzC4ck48DFQ9N9+6FmzWQIk9hR5x1WrFRkL92e23Hf16/D8inRqKlT6HJ9FsWJJdGlMH9WaEnYI0+wsUY1fEsX451Wj/Js7fKkpCrObj/NxU++JfTIj/hwk72F6vJ9UGd+bhGMclF5GmsuhC2TBC4ezrx50KsX4xjGsV7j+O478DLRw5HXmZdaw95dKRyfuR3vX1bR8HIY5bhodv/bFGKFbwcWPPk0R2oFZPmdqVa/EI5AJvKInDt4EEPf/vxKcyI7j2HJD+Ynm2Uk6dxMTc8y5rumF0OHjMCr4jdsm7mXhE27cMGAm4vG1UXj4gruvkWoOepFpq06ZJXytELYGmmBi6zi4kgNCSXmzE3altuLR59TXEmKt3jNkLy03G2l3ooQBcVcC1xGoYg7tCbttT5w6hQ93BYT1+4sl5Pi82XlmrysFGPN8rRC2JI8daEopQYD/8bY/ThLa/2VRaIS1jF1Ki4rl/MhX3Cukx9pxc9k+bUlV67Ja82U3KwMJISjyXUCV0o9jjF51weSgZ+VUuu11icsFZzIHyar+xWKJe2991lDRxLefJ8Unw0mj7XUyjXWWtUUuGcAABavSURBVGZNCEeSly6UGsBurXWC1joV2AE8b5mwRH4xt0r5kZf6c8NQhAGlJ7K/3HZ8C7ubPN5SCVZWihEi7/KSwA8BTZVSJZVShYF2wD0rdiql+iqlwpVS4Vdzs76csChTfc/1j+2l5qGdjFEjUJ3OcTE+gfikVNxdsw49sWSClX5sIfIuT6NQlFJ9gLeAeOAIkKi1HmJufxmFYn3ZZ066phlYP/Md3G640KjBarybn8n8na+XO94ebrJyjRBWli/jwLXW3wPfp19gLBCVl/OJ/Je977lb5EYeu3GaLh5L8Gp0Psu+NxJTiBz5dEGHKITIoTwNI1RK+af/WRHoDCy2RFAi/9zd91z09i2G7FjMDpqyo0V1XDyydq3IA0UhbFteZ2KuUEqVBFKAt7TW1y0Qk8hHd8+c7LHyR/ySr/NZ+bGUqneJpNQ7+8kDRSFsn8zEdFanT5Na7TEWGl4icMeP/FNUVmoXwlZJLRSRRcLgYWBw5fd2n9GrKYBMjBHC3kgCd0YRERRe+xOfuXzCsG8CHry/nTA5QUnelIQDkwTuhGKHjUPhQ/Kg9+6p722vTC2IPHzlQQBJ4sJhSTErZ3P8OD6/LGe2+wAGf1LM2tFYTF6KYwlhr6QF7mRujPiCQniQ0PcdSpSwdjSWk9fiWELYI2mBO5OoKLxXzuNH19fp+0lpa0djUebGrMtYduHIJIE7kbhRkyAtjYs9hlLasfK3FMcSTkm6UJxFTAwe82ayVHWnz5jK1o7G4vKytJsQ9koSuJOIH/c1RVJu8dfzw+hR0drR5A9Z5EE4G+lCcQbx8bhMm8oanuPVz2tZOxohhIVIAncCCV/NpHDSdfa2Gk61ataORghhKZLAHV1yMqlffMk2mvPipIbWjkYIYUGSwB3c7XlL8bl5gR0NPiQoyNrRCCEsSR5iOjKtuTF6MteowdNftrV2NEIIC5MWuAMzbP8N/6gI1gS+Q+Mm6sEHCCHsirTAHdjl4ZPxoASP/l9Pa4cihMgH0gJ3VKdOUWbPauYX6cvEs3sIHLaeJuO3EhYRbe3IhBAWIgncQV36aCoGXPkhpC0XbiaguVNiVZK4EI5BErgjiovDZ8UcfnLpyvV6t7P8SkqsCuE4JIE7oOsTv6dw6k1m1uiGi0fqPb+XEqtCOIY8JXCl1BCl1GGl1CGl1GKllKelAhO5ZDCgp0xlJ09yq0Mpk7tIiVUhHEOuE7hSqjwwCAjVWj8OuAIvWSowkTsJi1dTIu4Me5sM4eNulaXEqhAOLK9dKG6Al1LKDSgMXMh7SCK3wiKiOTxoLKepzPygcgCM61yb8r5eKKC8rxfjOteWin1COIhcjwPXWkcrpSYC54BE4Bet9S/Z91NK9QX6AlSs6KB1TG1AWEQ0C75eyfLr+xhabAzXfWIYvjKWcZ1rs2tYS2uHJ4TIB3npQikOdAQCgXKAt1LqnhkjWuuZWutQrXWon59f7iMV9zVh4zFe/mU9N/BhxVNPAjLiRAhHl5culNbAaa31Va11CrASaGyZsMTD0mfP8lz0Zr736IWh5q3M7TLiRAjHlZcEfg5oqJQqrJRSQCvgqGXCEg+r766tAMxt0AF1V9kTGXEihOPKdQLXWu8BlgP7gYPp55ppobjEw4iL44Uja1jm2oXYUEPmZhlxIoRjy1MxK631SGCkhWIRuXTp09mUMcRxqPObVPBLk0V9hXASSmtdYBcLDQ3V4eHhBXY9p5CaSoxvFY4kVqbG5R2UMj13Rwhhx5RS+7TWodm3y1R6O3d99gpK3jrH4bbvSfIWwslIPXB7pjXx//clV6lG6686WDsaIUQBkxa4HUvYtIsKF/9ge50hVH1UXkohnI38r7dTYRHR/PbqaK5RktmhtaXGtxBOSLpQ7EBYRDQTNh7LHF3S4jE/Ijb8j7WXt/B50Xe5UvIGw1ceBJBRJ0I4EWmB27iwiGiGrzxIdGxi5qo6C3efo8+6VdzCmx+bGVebl2nzQjgfSeA2bsLGYySmGLJsq3L5PJ0ubuRb73+TWDM5c7tMmxfCuUgCt3GmkvKbP6/mFt780PJZmTYvhBOTBG7jsiflqpfP8/yl9NZ3jTutb5k2L4TzkQRu44a2rZ5lVZ2M1vfJfj0IKC4LNQjhzGQUio3LSMoTNh7D69hxnr+0kZl+Q5k5sV6W7hMhhPORBG4HOoWUp1NIeU48MYtbeFN92lBJ3kIISeC2JvuY74yKgsmRR6gSvpQF5YfxygtS9EQIIQncpmSM+c4YNhgdm5g5Qafmv8dQFm8qTXlXWt9CCEAeYtoUU2O+E1MMrP8mjEf3LSEs4G2adpbWtxDCSBK4DTE15ts1zcCgxd9wngCqzRkurW8hRCZJ4DbE1EScV3ZspHbiX6xtPZWGbYpaISohhK2SBG5Dso/5LnPjKu//8SP/dWtPt8WdrBiZEMIW5TqBK6WqK6Ui7/qKU0q9Y8ngnE2nkPKM61yb8r7GCTqj1yzERadxa9w3lCwlfSdCiKxyPQpFa30MqAOglHIFooFVForLaWWM+b6xaB3FLmxmRuVx9H2vsrXDEkLYIEsNI2wF/K21Pmuh89kNc+O28yQhgZR+AzlMTZ5aJcMGhRCmWaoP/CVgsYXOZTdM1eoevvJgnlfHOd93DKXiz7Lz5e+oWaeQZYIVQjicPCdwpVQh4DlgmZnf91VKhSulwq9evZrXy9kUc+O287KwQuKKDZRb+AU/eb9Gz5lN8xqiEMKBWaIL5Vlgv9b6sqlfaq1nAjMBQkNDtQWuZzPMLaAQHZtIk/FbH7pbxbD/AHTrxgGC8Vs8FW9vS0cshHAkluhC6Y4Tdp+A+QUUFDx8t8qFC9xs3oEYQzEOjFlLi38VsXi8QgjHkqcErpQqDLQBVlomHPuSfdw2GJN39o8ZD+xWuXWLa02ew+3mdeZ1XUfvEVLXWwjxYHlK4FrrBK11Sa31DUsFZE+yj9su7+t1T/LOYHa9SoOBa8/2pPiZCMYFL2Howjoy6kQIkSNSjTCPMsZtZ2gyfivRJpK1ye4Wg4Ebfd6l1G9hfOo3hfe3dcDdPT+jFUI4EplKb2GmulVMrlcZFcX5Os0o9uNUprgOZOWrT7DjTN6GHwohnIu0wC3s7iXQzI5CCQvj1suvUzwxmT6FZrD+hVp4uv3D8JU3spxDCCHuRxJ4PsjerZIpIQHD4HdxnT2Do9Sjt/+3xL4Qh2fRWODOw05J4EKInJAEXhCOHIGFC0mdOx+3C+f5gqF8Ue91vFucws0162NPsw87hRAiG0ng+SE5GU6dgrVrYeFCOHAAAy5sUa35wm020d3KUuHRc1xPuHfMirmx5UIIkZ0k8NxISYELF+D8+Ttfp07B338bv86ehbQ0AA4VacAMprLMrTO3Hk/Fp8HfuPuew5CkcHdVpBjuJHGTDzuFEMIMSeA5kZwMv/0GGzYYv44dA5219ZzmW4JYv6qc9WzEnwE92XWpCluTm+BStiop1Y9TqMohPD1TM/dPSdP4ernj7eFm2UqGQgin4XQJPMflX7WGNWtg7lzYvBni48HDA928OTefeZFTKRU4GFuB3dEV2HysAscvFoVYcHeHkBBo2Am+bQ+tW0OVj06YHK95IzGFyJFP5/ctCyEclFMl8IzyrxkVBDPqlEC2oXvR0fDWW7B6NQQEQM+eXAxpx5wzLZmz1JtTG427KQWPPQb1W8Gb9aBhQ6hTBzw9s163nK9Xzif3CCFEDjlVAr9f+ddOIeWN/dazZ8PQoZCSQtKYCSz0e4e5C9zYOd2YsFu2hHfegXr1IDiYHFUMHNq2epY3DpD+biFE3jlVAjc3RO9CbKLx4WOfPrBjB7RowW+vzuLF4VW4dMnYyh43Dnr0gAoVHv66OZrcI4QQD8mpEri5rox6abHQpAkkJZHy7Sw+PN6Hya8pataE5cuhcWPyXGDK7OQeIYTIJaeqhWKqTknZlHh++Ok/kJzMyXm/Ezr9DSZ/pRg4EMLDjXldqgMKIWyRU7XAs3dlBBZWLFsxgcIXo3mpySx+6lIdV49kRkyJ49NBpawcrRBC3J9TJXC4qyvDYIAXXkAf3E+/xuNZur0nHhVi8PtXJCuuphIaUVu6PIQQNs2pulAyaQ2DBkFYGJ83HcLsP96lkH8c/l3CcS1yO88LEwshREFwuhY4AJMnw7ffEv/WB3w8bxQuhVLx6/IHLoXuDPMztTAxyEgSIYTtUFoX3ELxoaGhOjw8vMCuZ3LWpU8S1KpFatt2NLmwnPADGv+Xf8ejTFyWY7OvbenuokBxT+2ScZ2lq0UIkb+UUvu01qHZtztsF0rGrMssq8Ov+JNLvfui3d0ZYPiGP/a58MHnsfhWuJXlWFMLE6ek6SzJG3KwWLEQQuSjvK5K76uUWq6U+kspdVQp1chSgeWVqVmXTx79nTI7t/JL49HMWl+OCRNg3JCSOV6Y2BSp3y2EsJa89oFPAX7WWr+glCoEFLZATBaRPbF6piQxcvNM/ioeyPNb36ZHD3j3XePvcrowsSlSz0QIYS25boErpXyApsD3AFrrZK11rKUCy6vsiXXA/5YREHeF97yn4F7YnS+/ND9Bx9SEH3cXY/3uu0k9EyGENeWlC+UR4Crwg1IqQik1WymVg9JOBePuJFz5n2j67V3Biort2BD1L0aPhtKlzR/bKaT8Pd0qE7oGM+GF4Czb5AGmEMKacj0KRSkVCuwGmmit9yilpgBxWutPsu3XF+gLULFixXpnz57NY8g5FxYRzYSf/2LsrA+od+EYjYr/RVqJACIjjXW7hRDCHuTHKJQoIEprvSf95+VA3ew7aa1naq1Dtdahfn5+ebjcw+sUUp5d1W/Q7PR+drYay6FLAXz9tSRvIYRjyHUC11pfAs4rpTI6gVsBRywSlaWkpcGIESRXf5yuWwfQtauxnrcQQjiCvI5CeRtYmD4C5RTwWt5DsqB16+Cvv5hWfzEG5cbEidYOSAghLCdPCVxrHQnc0y9jM774gsQylRm69wVGjYGKFa0dkBBCWI7DzsRk1y7YtYuvC71HxUA33n/f2gEJIYRlOW4xqy++IMWnJKPPvcak6fcuNCyEEPbOMVvgR4/CmjX85D+QIv7evPqqtQMSQgjLc8wEPnEiaZ5eDD45kLffBi+Z7S6EcECOl8AvXID589lSqQ+JhUvx5pvWDkgIIfKH4yXwKVPQBgMDTr7LG29AyZLWDkgIIfKHYyXwGzdg+nQiH32R0wQyZIi1AxJCiPzjWKNQZsyAuDgGpQ6la1eoXNnaAQkhRP5xnASelgbTp3PukebsPFWXfUOtHZAQQuQvh+lC2TlrOZw+zfDo3hSr8g/nVLS1QxJCiHzlEC3wsIho3L7+juvuPqy43Y1idQ8wfOV1AKnXLYRwWA7RAp+xci9t/trJIveXSCuVjGfgNVlwWAjh8BwigTf83894GFKZkfAWReqcz1wqTRYcFkI4MvtP4FrT89Am9hUO4qBrLbxr3un7lgWHhRCOzP4T+B9/UOXyaWYn96Nw9Uu4eqUAsuCwEMLx2X8C//57Uj0KszC1J4GNr8iCw0IIp2Hfo1Di42HRIjYVfxG/wj78OSMEF5cQa0clhBAFwr5b4MuWQXw8n156gz59wMW+70YIIR6KXaa8sIhomozfSvgnEzhRuAq7VSN69bJ2VEIIUbDsLoGHRUQzfOVBvE4eIzT6KLNS++BV9Rp/XJGZl0II55KnBK6UOqOUOqiUilRKhVsqqPuZsPEYiSkGXvxzEynKlbnJfShc+5xM2hFCOB1LPMRsobW+ZoHz5MiF2ERc0gx0PLqDnwu3JgYfAh65woVYXVAhCCGETbC7LpRyvl40OneQ0vH/8OOt1ynyeBTKVcukHSGE08lrAtfAL0qpfUqpvqZ2UEr1VUqFK6XCr169msfLwdC21elydAdxrkVYx3MUCTovk3aEEE4pr10oTbTWF5RS/sAmpdRfWutf795Baz0TmAkQGhqa536OTo+VIOXv/7HIvTOUTaDyI5qhbWXSjhDC+eQpgWutL6T/eUUptQqoD/x6/6PyaO1a3G/F8yO9mPJxCfr1a5mvlxNCCFuV6y4UpZS3UqpoxvfA08AhSwVm1oIF3PAux+9uzejaNd+vJoQQNisvLfDSwCplrN3qBizSWv9skajMiYlB//e/LPJ4h2fau1KiRL5eTQghbFquE7jW+hQQbMFYHmzZMlRqKjNSe/BxzwK9shBC2Bz7Kma1YAFRvrU4bQimQwdrByOEENZlP+PAz5yBXbuYndCDF7oqPD2tHZAQQliX/STwRYsA+CH5ZXr0sHIsQghhA+yjC0VrWLCAwyWewuBViWbNrB2QEEJYn320wCMj4ehRvontSffu4Opq7YCEEML6bD6Bh0VEs3jQWJKVG0vTXqBc6GVrhySEEDbBprtQMmp/B1eoy36fOsS7uzPj4F4CH5Wp80IIYdMt8Iza37/5NGDijRF414omKdUgtb+FEAIbT+AXYhMBSDhibG1717iQZbsQQjgzm07gGTW+XYsk4V37PG7FErNsF0IIZ2bTCXxo2+p4ubtSJCiKUu3+BJDa30IIkc6mH2JmPKicsPEYF2ITKefrxdC21eUBphBCYOMJHIxJXBK2EELcy6a7UIQQQpgnCVwIIeyUJHAhhLBTksCFEMJOSQIXQgg7JQlcCCHslNJaF9zFlLoKnM3l4aWAaxYMxx7IPTsHuWfnkJd7rqS19su+sUATeF4opcK11qHWjqMgyT07B7ln55Af9yxdKEIIYackgQshhJ2ypwQ+09oBWIHcs3OQe3YOFr9nu+kDF0IIkZU9tcCFEELcRRK4EELYKbtI4EqpZ5RSx5RSJ5VSw6wdT0FQSp1RSh1USkUqpcKtHU9+UErNUUpdUUodumtbCaXUJqXUifQ/i1szRkszc8+jlFLR6a91pFKqnTVjtCSlVAWl1Dal1FGl1GGl1OD07Q77Ot/nni3+Ott8H7hSyhU4DrQBooA/gO5a6yNWDSyfKaXOAKFaa4ed7KCUagrEA/O01o+nb/sC+EdrPT79zbq41vpDa8ZpSWbueRQQr7WeaM3Y8oNSqixQVmu9XylVFNgHdAJ646Cv833u+UUs/DrbQwu8PnBSa31Ka50MLAE6WjkmYQFa61+Bf7Jt7gj8mP79jxj/4TsMM/fssLTWF7XW+9O/vwkcBcrjwK/zfe7Z4uwhgZcHzt/1cxT59JdhYzTwi1Jqn1Kqr7WDKUCltdYXwfgfAfC3cjwFZaBS6s/0LhaH6U64m1KqMhAC7MFJXuds9wwWfp3tIYErE9tsu9/HMpporesCzwJvpX/0Fo7pO6AKUAe4CHxp3XAsTylVBFgBvKO1jrN2PAXBxD1b/HW2hwQeBVS46+cA4IKVYikwWusL6X9eAVZh7EpyBpfT+xAz+hKvWDmefKe1vqy1Nmit04BZONhrrZRyx5jIFmqtV6ZvdujX2dQ958frbA8J/A+gmlIqUClVCHgJWGPlmPKVUso7/eEHSilv4Gng0P2PchhrgF7p3/cCVlsxlgKRkcjSPY8DvdZKKQV8DxzVWk+661cO+zqbu+f8eJ1tfhQKQPpwm68AV2CO1vozK4eUr5RSj2BsdQO4AYsc8Z6VUouB5hjLbF4GRgJhwE9AReAc0FVr7TAP/czcc3OMH6s1cAbol9E/bO+UUk8CvwEHgbT0zR9h7BN2yNf5PvfcHQu/znaRwIUQQtzLHrpQhBBCmCAJXAgh7JQkcCGEsFOSwIUQwk5JAhdCCDslCVwIIeyUJHAhhLBT/w+iFSHZm/LsmQAAAABJRU5ErkJggg==\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.083128\n",
       "x1                  0.490733\n",
       "np.sin(x1)          0.521571\n",
       "I((x1 - 5) ** 2)   -0.018911\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.119806\n",
       "1    10.983789\n",
       "2    10.731476\n",
       "3    10.409481\n",
       "4    10.079160\n",
       "5     9.801596\n",
       "6     9.622638\n",
       "7     9.561679\n",
       "8     9.606905\n",
       "9     9.718192\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
}