{
 "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.989\n",
      "Model:                            OLS   Adj. R-squared:                  0.988\n",
      "Method:                 Least Squares   F-statistic:                     1326.\n",
      "Date:                Tue, 17 Dec 2019   Prob (F-statistic):           1.18e-44\n",
      "Time:                        23:39:41   Log-Likelihood:                 8.4408\n",
      "No. Observations:                  50   AIC:                            -8.882\n",
      "Df Residuals:                      46   BIC:                            -1.234\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.8426      0.073     66.676      0.000       4.696       4.989\n",
      "x1             0.5258      0.011     46.943      0.000       0.503       0.548\n",
      "x2             0.5333      0.044     12.111      0.000       0.445       0.622\n",
      "x3            -0.0219      0.001    -22.256      0.000      -0.024      -0.020\n",
      "==============================================================================\n",
      "Omnibus:                        3.262   Durbin-Watson:                   2.218\n",
      "Prob(Omnibus):                  0.196   Jarque-Bera (JB):                2.222\n",
      "Skew:                          -0.402   Prob(JB):                        0.329\n",
      "Kurtosis:                       3.648   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.29539098  4.80736761  5.27727471  5.67604936  5.98511733  6.19944471\n",
      "  6.32836499  6.3940456   6.42784586  6.46516473  6.53962487  6.6775485\n",
      "  6.8936328   7.18853511  7.54876503  7.94890079  8.35576599  8.7338804\n",
      "  9.05128925  9.2848133   9.42385717  9.47215005  9.44713309  9.37709381\n",
      "  9.29651856  9.24042651  9.23861612  9.31076928  9.46321769  9.68790257\n",
      "  9.96369866 10.25988507 10.54119263 10.77359851 10.92991608 10.9942616\n",
      " 10.96466386 10.85338805 10.68492063 10.49194626 10.30997784 10.17152178\n",
      " 10.10073735 10.10946733 10.19529231 10.34192817 10.52190155 10.70106359\n",
      " 10.84420053 10.92081924]\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.89466562 10.72741816 10.43998979 10.08003848  9.71029886  9.39322261\n",
      "  9.17568813  9.07752305  9.08664973  9.16204231]\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd1yV1R/A8c/hshHFAZm4sJRf7oGmWe6t5cpKzdEyLXOUlqalNtTSshy/cuQqtRxJqak/c1RqahjurTnAhSAukHE5vz8uEOMiCBfu4Pt+vXgpD89zn+/jg9977nnO+R6ltUYIIYT9cbJ2AEIIIXJHErgQQtgpSeBCCGGnJIELIYSdkgQuhBB2yrkgT1aqVCldsWLFgjylEELYvb17917TWvtm3F6gCbxixYqEhIQU5CmFEMLuKaXOmdsuXShCCGGnJIELIYSdkgQuhBB2qkD7wM1JSEggLCyMu3fvWjsUkQPu7u6ULVsWFxcXa4ciRKFn9QQeFhaGt7c3FStWRCll7XDEPWitiYyMJCwsjICAAGuHI0ShZ/UEfvfuXUnedkIpRcmSJYmIiLB2KDYvODScKRuPczE6ljI+HoxsG0iXOv7WDks4GKsncECStx2Re5W94NBwRv94kNgEIwDh0bGM/vEggCRxYVE2kcCFcCRTNh5PTd4pYhOMjP/5cKFqlcunkPyX7SgUpdR8pdRVpdShNNt6KKUOK6WSlFJB+Rti/jMYDNSuXZtq1apRq1YtPv/8c5KSku55zNmzZ1m6dGkBRSjsycXoWLPbo2MTCI+ORfNvqzw4NLxggysgKZ9CCsv1WktOhhEuBNpl2HYI6Ab8bumAshMcGk7jyVsIGLWOxpO3WOQXwsPDg3379nH48GE2bdrEL7/8woQJE+55jCRwkZUyPh452i82wciUjcfzORrryOpTiKNer7Vk24Witf5dKVUxw7ajUPD9oQXRt+jn58ecOXOoX78+48eP59y5c/Tp04c7d+4AMHPmTB577DFGjRrF0aNHqV27Nv369aNr165m9xOFzztNy/PblHk8ue9XSt+6xslS5TnuW4HjvhU5XqoCF3weQCtT2ymr1rq9y+q6LkbHkrIImDxOybt87wNXSg0ABgCUL18+T691r3d1S/atVapUiaSkJK5evYqfnx+bNm3C3d2dkydP0rNnT0JCQpg8eTJTp05l7dq1AMTExJjdTzi+4NBwpmw4ht/R/fQ9sY1Oh3/jqds3uezjx+GSFah75RRPHvsjdf/TJfwZ3e4N9pSrnuPWur0p4+NBeIYkrjW4h5enfn2YNAlat7ZScA4k3xO41noOMAcgKCgoTwtw3utd3dJS1gpNSEhg8ODB7Nu3D4PBwIkTJ8zun9P9hGMJDg1nwg9/8emqSbQ+tYdYZzfW/6cxRQcM4HLJHpw568Sl8nDp5j/s2raeKtFHGLRnJcuXjmJZvY54T5tq7UvIFyPbBqZ+WtZGxZ3D/tzZXZE2Udt43et9qn3kBpsrQsU0XxUqgIdjvqHlF7sahWLuXT1luyWdOXMGg8GAn58fEyZM4IEHHmD//v0kJSXh7u5u9php06blaD9hP8yNogDSbTPcuM7c796jbvgxJjftz/xyPbh89BFi3ylLYmzaR0wBwGsopyRmP/g2X1YYxfN/f4/qfQBmz4b27a1yjfmlSx1/kpJgxAe3ubalJM/fXs5Il+4EcBrt44+66AGf/wgJCekPfPVV+PJLcHOzTuB2xq4SeNp39RQeLobU/1iWEBERwcCBAxk8eDBKKW7cuEHZsmVxcnJi0aJFGI2mc3t7e3Pr1q3U47LaT9gnc89bRq7YDwoSjKZPZ4kXLvDN8nEEXA/n5cc+ZsnJgcT/VhwMRjwrX2HDzDI0bAgXLsC5c6avf/5xYuHC8vQ9sJRTbYcw9sxLuHToAAMHwsyZYDBY87ItJiYGlk/0o+/PcxnuPINiRKHrNIC3JqK6dQNnZ0hKgkuX4OxZ09eOHfDVV7B/P6xaBWXKWPsybJ5K6SrIcgellgHNgFLAFWAcEAXMAHyBaGCf1rptdicLCgrSGfuFjx49yiOPPJLjgPNjbKnBYKBGjRokJCTg7OxMnz59ePPNN3FycuLkyZN0794dT09PmjdvzowZM7h9+zYJCQm0a9eOa9eu0b9/fzp16mR2P0d0v/fMHjWevMXsp70UAVHhfPvDexS7e4vna03n55AXcfaJxbvOWbyqhVO+jDM7RrUwe+zt2zB5Mnz2GbgRx7o6Y2n851To3RsWLjQlNzsWHg59O0Yydv/TNGcbunNn1IgR0Lhx9k8uV66E/v3B29v098aNCyRmW6eU2qu1zjRkO9sEbkmWSODC+grDPQsYtY6s/mdUu3KaxT+8h1aK7lXm8/u+brgHXMW3y984uRrxcDEwqVuNbBsW587B86/GsH2jJ2M9xvFh7AeEte1M2bUr7TaJ790Lw9sfY/61JwkwnMcwfx706ZP68xw1wA4fhi5dTP9A06ebulUK+ZCVrBK4ff6WCJHPsnreUiz2FnNWfcRdZzc6ll3O/n3NKF7rIoE9jnP5lvG+PhWGRoVzrf5BHijuw8TgMSS4uDJ541gONu3I6x3e4sKtBLuawbhqFSzotYk1CT3w8nHFsHYrpBlKm+NhwNWqwZ49pk8kgwZBSAjMmiX94mZIAhfCDHPPW1wUTNkwHd8712nhv4b9R5tR4rEzzJvhRte6ze/7HCnDYt0rRFL6+Z1MXTGcxFuuTN35NqOv3WbIUyMJj8Yu6qgM+zCKuPe/J5gh/ONbgTMLltL2sUfT7XNfw4CLF4c1a2DcOPj4Y4iOhhUrCn1LPCNZ0EEIM7rU8WdStxr4+3igAH8fD1bqfbQ58ScTSr3PjgttCXjyBN/MdKNr3dwl1rTDX11K3uHB53cyw+8VhjKN9id28uXPU3BKMqbWUbH0DGRLGflJFPr9xXzF62yrVI8uvSczbNf1TDHe9zBggwE++gg++cTUvJ8929Kh2z1J4EJkoUsdf3aMasE/kzuyo10Jak37gNP/6cjEq2P54gvFmZ+r5KlVnHH4q6FIHA/03MWcSr0Yzud0OLGT9zfPBa1tto7KTz9B9OiVfMlwfnm4Ma92H8MdN0+z0+azGu6b7TDgESOgbVsYPhwOHbr3voWMJHAhsnP7Njz7LHFFS/HYiYX06KEYMiTvLzuybSAeLumHDRpcjfh1D2FujaeZylv0/3str+5ZlelYW6grsmEDrOq+lNl6IFsq1mdolxEYnf69nowta3PXm6NhwE5OsGgRFC0Kzz1nGqMoAEngQmTv9dfRp07RM2kpxR4qxbx5lumKNddN07theTzdnCjZ/gATqr/BMp5j9LaFPHVkW6bjrVlHZcsWWPDUauYb+xJaoRaDuo0iwZB+mb2MLWtz15uT0ToAPPAAfPutaYTKW29Z8ErsW6F/iBkZGUnLli0BuHz5MgaDAV9fXwD27NmDq6urNcMT1rZkCSxezOKK41h/uSm7t5oagpbSpY5/pgQWVKEEUzYeR7U/xEjDJErvv8zUtV8S4VWcPyvUSt3PXNdDQdTgnrQggj9e2UOw8VmO+v+HdZO+wulINORggp25682xNm1g5EiYMgVatYLu3XN7CQ5DxoGnMX78eIoUKcKIESPSbddao7XGyUk+sIBt3bN8FRUFgYGcd69MQNgfzJtv4IUXCjYEoxGeaneRyb+2pYLhLM/0ncQxvwCzY80zDtMDcjwmPSsZ3xACDeX5593zrEvsxJlSZejV+2MSvIvRvZ4/W49F5P/iDfHx8PjjcPIk7Ntnqp9SCGQ1DlwyUhZOnTpF9erVGThwIHXr1uXChQv4+Pik/vz777/n5ZdfBuDKlSt069aNoKAgGjRowK5du6wVtrCkd98lKeo6ncK+pt8LBZ+8wTQQ46f1ZRjfZBU3jcVYuGQCDYzXzSZlS9fgzrgowz8nnDky9go/J3YmrLgvfXt+wE33IsQmGNl6LOLfB76jWuTfkEdXV1i2zPTO1rs3JCbmz3nshE11oQwbZnpTtaTateGLL3J37JEjR1iwYAFff/01iff4RRkyZAhvv/02DRs25OzZs3Tq1IlD8rTcvu3Zg54zh2+8hpFUoSYzZ1ovFGdn+GFzFd558n+M3tCE2fPew2fo9kyt46ym/ue2rzztG0L8tSL4LfPkl/g2RHl706fnB0R5FsvzOXLloYfg669NCfzjj01jxQspm0rgtuahhx6ifv362e7366+/cvz4v62c69evExsbi4eUxrRPRiMMGsRNz9KMuD2ejfPA09O6ITk7wydrqvJRlw28ua4FYdVb8skL47nk5g2YhhYqMDv9P7fVOlOScsJ1T4otK8HGu62I8zLQu9eHXPEuZZFz5FqvXrB2ramweJ8+UKlSwZ7fRthUAs9tSzm/eHl5pf7dycmJtM8L7t69m/p3rbU88HQkX38Nf//Nq3xP38FFadjQ2gGZODvD+z8H8eXTaxm4ui3/nf0xfQa8x2130++phkxJPC/VOsv4ePDPCWfclz/AxthWGNzj6dVzEuE+pdPtZ+mKoDk2ZQr8/LNpVMrq1QV/fhsgfeA55OTkRPHixTl58iRJSUmsTvML06pVK2bNmpX6/T5L9wOJgnPlCnrMGP4s0ortZZ7h44+tHVB6Tk4wbFUT+gTOpHrsMWbP/Qy3uLjUn2vI3TA9M1p41aTod6XYFNOeoi7R9On5ARdLV6R3w/IWO0ee+PvDu+9CcDD8+mvBn98G2FQL3NZ98skntGvXjvLly1O1alXikv/jzJo1i0GDBrFgwQISExNp3rx5uoQu7MjIkSTdjqG/cSazvlMWHTJoKUrBxX4BDFj8Kd8ce4t5X09lcO/h3Cjlib+PR5ZlbLOTtk9dHX2Y8msu86fqhvZIon/3CdysUo1JtlZY68034ZtvYOhQ0wM0F5fsj3EkKUPkCuKrXr16OqMjR45k2iZsm8Pes23btAY92fCu7tbN2sHc2+q/w/R/xq7XL9eYrGNx0/+o8rpNu+/06r/D8vR65Ueu1d71zujnWazjcNHRFR7W+vRpC0dvYatXaw1aT59u7UjyDRCizeRU6UIRAiAxET1kCFfcyzPNcwzTp1s7oHtLmdV4pFd9uj41HWdlZPWGV7gxZjtpHs/k2JSNx7lx0YOry+vz1t6FfEtfQso/wrP9ptr+A8LOnU0Te95/H65ds3Y0BUq6UIQAmDcPdeAAr7OCcZ974m9DvQRZ+XdWYwti/3mKi42fpt/655hXfj+Byz+kcRMDOZl7duoU7P82kFJH4lmsBtGDH1leoxVj2r5OYpwdpAilTOto1qwJ771nWpatkJAWuBDXr5M0Ziw7XZoQ/mh3Xn3V2gHdP4+A0jx0dgvn2w3g5YhJFGkexDu+83l3WAwhIZBxwnVMDBw5AgMGQNvAs0w7Op6TqgpPqTV80rQfb7cfSoLBpeCHB+ZW1arw+uswZ45pTc1CQhK4EB98AFFRvJbwBece2c4Tn9pWve0cc3Wl/PrZxM75lgD/BKZEvcTIL/35vf6btK5wgqZN4ZEqRsp5RxPodYHe1UJp8M2rnKAyLzov5od67WkycC5fNewBSllveGBujR9vWghiyJDM71gOShK4KNyOHsU4YyZzeZkzdYvj+sBNm6q3nRserzyPz4WD8NtveHZtyzDDDH69EMiG7V4cPenMhdvFuUB5QqnLS84LMQx6FcOZ0xSZ8xXO5cpZf3hgbhUvbpqZ+fvvptV7CoGcrEo/H+gEXNVaV0/eVgL4AagInAWe0Vpfz+5ktl7MSuSMI90z3a49N/63k0D3w7gNOIyT+78lE/IyJM+mXL5sqqd97RoUKQLe3oReT2T1qZtsKl4Zp/Ll7GbdzWwZjRAUZCpEduwYOMhs6LwUs1oItMuwbRSwWWtdGdic/L3dCgsLo3PnzlSuXJmHHnqIoUOHEh8fD8C2bdvo1KlTpmPWrl1LnTp1qFWrFlWrVmV2ASz3tHDhQgYPHgzA119/zeLFi7Pc9+zZsyxdujT1+5CQEIZYYhUCR/LLL6iNG5igx2NscS1d8gbr1tu2qNKl4Z13TDMXx40juPmz9KImiys04lLRUnb/iSMdgwGmTYPz5+Hzz60dTb7LNoFrrX8HojJs7gwsSv77IqCLheMqMFprunXrRpcuXTh58iQnTpzg9u3bjBkzJstjEhISGDBgAGvWrGH//v2EhobSrFmzXJ3faDRmv5MZAwcOpG/fvln+PGMCDwoKYrqtj40rSPHxGIcM56QhkMUVeuNVLXPyspsHeNkIDg1Pt57mhDWHLVq10OY0awZdupjqpFy+bO1o8lVu+8Af0FpfAkj+0y+rHZVSA5RSIUqpkIiIiFyeLv9s2bIFd3d3XkiuFWowGJg2bRrz588nJoulm27dukViYiIlS5YEwM3NjcDAzA97xo8fT58+fWjRogWVK1dm7ty5gKlV37x5c3r16kWNGjUA+O6772jQoAG1a9fm1VdfTU3sCxYsoEqVKjRt2pQdO3ake+2pU6cCptK3rVq1olatWtStW5fTp08zatQo/vjjD2rXrs20adPSfZKIioqiS5cu1KxZk4YNG3LgwIHU13zxxRdp1qwZlSpVcuyEP306htMnGK6n8d4nGk/XXCz1ZQcyloQNj47lekyC2X0d5hMHmD5txMfD2LHWjiRf5fsgT631HGAOmPrA77mzFerJHj58mHr16qXbVrRoUcqXL8+pU6fMHlOiRAmeeuopKlSoQMuWLenUqRM9e/Y0u+DDgQMH2LVrF3fu3KFOnTp07NgRMK32c+jQIQICAjh69Cg//PADO3bswMXFhddee40lS5bQunVrxo0bx969eylWrBjNmzenTp06mc7Ru3dvRo0aRdeuXbl79y5JSUlMnjyZqVOnsnbtWsD0ppFi3Lhx1KlTh+DgYLZs2ULfvn1T67ccO3aMrVu3cuvWLQIDAxk0aBAujjY9+fx5jO+P5xc6UXlIe4Y9CxWr1Mj3lWyswVyN8Kw4yicOAB5+2DQa5fPPYfBgUx5wQLlN4FeUUg9qrS8ppR4ErloyqIKktUaZWeAwq+0p5s2bx8GDB/n111+ZOnUqmzZtYuHChZn269y5Mx4eHnh4eNC8eXP27NmDj48PDRo0ICAgAIDNmzezd+/e1NK1sbGx+Pn5sXv3bpo1a5a6xNuzzz7LiRMn0r3+rVu3CA8Pp2vXrgC4u7tne83bt29n1SrTQrktWrQgMjKSGzduANCxY0fc3Nxwc3PDz8+PK1euULZs2Wxf054kDR5CfJzmA9+Z/DretC1PS33ZsJy2qh3lE0c6Y8eaHt6++SZs3myZhUxtTG4T+M9AP2By8p8/WSQaK9STrVatWmoyS3Hz5k0uXLjAQw89RGRkZJbH1qhRgxo1atCnTx8CAgLMJvCMbwIp36ctVau1pl+/fkyaNCndvsHBwfd8E0k59n6ZOyblPG5ubqnbDAbDPReysEs//YTTmp94n095d3YFihXL/hB7ltVCDz4eLni5OTvcJ450fHxgwgTTBJ+ffzZNuXcw2faBK6WWAX8CgUqpMKXUS5gSd2ul1EmgdfL3dqlly5bExMSkjugwGo289dZb9O/fH88sqvjfvn07XZfEvn37qJDF2nw//fQTd+/eJTIykm3btpldIKJly5asXLmSq1dNH2SioqI4d+4cjz76KNu2bSMyMpKEhARWmBnbWrRoUcqWLUtwcDAAcXFxxMTE4O3tza1bt8zG1KRJE5YsWQKYulZKlSpFUVssu2dpt2+TMOgNDqoaXOg2jOQPLQ5tZNtAPFwy9++Pf6pawSyBZm0DBphmaY4YYeoTdzA5GYXSU2v9oNbaRWtdVmv9jdY6UmvdUmtdOfnPjKNU7IZSitWrV7NixQoqV65MlSpVcHd3Z+LEian7bN68mbJly6Z+hYaG8umnnxIYGEjt2rUZN26c2dY3QIMGDejYsSMNGzbkvffeo0yZMpn2qVq1Kh999BFt2rShZs2atG7dmkuXLvHggw8yfvx4GjVqRKtWrahbt67Zc3z77bdMnz6dmjVr8thjj3H58mVq1qyJs7MztWrVYtq0aen2Hz9+PCEhIdSsWZNRo0axaNEis6/raPT743C5dIGRRb7mi1kO1q+fhZSiVzZRv9sanJ3hs89MBV8csMSzrEqfj7Ja5d7e2eU927ePpHpBzE16CbcFs+nf39oBiQLVvj3s2mVK5Mmjx+yJrEovbE7G8cn5NpHEaCTuxYFE6hL8r9kkitUsoPMK2/HZZ3DrlsMtgGwHtSLt1/jx460dgs1KGZ+cMsQtZTYgYPGP93rGTNxCdzPK9VvavBXPu6sL5rzChlStCoMGmbpReveGRo2sHZFF2EQLvCC7cUTeWOpemRufnC+zAXfvJmnESNbQiRqTe7P40BHHnoUosjZxIpQrBy++SK5WvbBBVk/g7u7uREZGShK3A1prIiMjczTWPDtZjU8Oj461XPdGVBTxXZ4hLKkMsxst4o0hKsvzOtQsRGGetzfMm2cqcuUgn46t3oVStmxZwsLCsMVp9iIzd3f3TBN70i6Gm9MxxVmNT1aQuj1P3RtJScT37Iu+fJnXfLczf3UJDIasz+tQsxBF1lq3hpdfNk21794dzAzrtSdWH4Ui7FvGvmwwjTPObqiaueMUYO63MTdlXZMmTsZpzGiGGmbSa8frPPpo3uIVDuTGDaheHYoVg717Ic3kNVslo1BEvshtX7a58clZNSXuu3vjt99g7BiW8RzVZr2WmryzOq8k70KmWDHT0muHD8NHH1k7mjyxeheKsA256QYBUOfO0eufvdS+eJwIr+KcLV6GMyX8OXenjGlZq3uUAshYf6Tx5C157964dIm7XZ/jrK7Mjr5zaNXgIo0nZ74uSdiFXPv20K+fqeRs166QxSQ5WyddKOL+uhW0hq1bYe1a2LABjh4FINKjKEXj7uCSlKY1XqIEPPecafhW9eoWjcPsG07Ru8S1aEfi+Yu8XG0XXeeVYNw66S4RWbh+HapVA19f+OsvcHW1dkRZyqoLRRK4yLLlm6nvOS4OBg6EhQvBzQ3dpClbSj7O2BP12XetHi5ud6nkeopA5xM8os5QMzqUjmHbcEtM4FrtBswMbMWysvUpVapoli18c4kZSLet+X98WbU3PF1iDoo4zaIVH3L3tqZf8TV8ta8Rzy3J4XWJwiulyNW4cTY9MkUSuMhSwKh1ZvufFfDPZFP9ci5dgm7dYNcuZlV7lQnGt7h+riKJsS44OWk8HrxFfKwTxlvuJCX82zPn5xXGK8X/y8s3FlPxVjjXPIuxqG4nlj/amdG9H8u2JZyTh53NTocwc/UnRBgfoEuJH4l85joVKmA2eWe6LiGefx6WLoUlS6BnT2tHY5YkcJGlrFrgKSVHfY/sY+5PEykaF8OLxWey9PKLOHnG4VEpgmKVr/HJUD+eb2oq0vXYpC2EXU0g8aY7Cde8uXO0DLGn/VBJ0NZ7DcOdp9Hm+m/cdvXg50c70ev7L8BMga/sYkvR/cCvfLJ+BgeoSVf/b9FPh6WubWnJUS3CgcXEmPrEd+wwrWZvg2UqZRSKyJK5kqMuToo78Yk8un0tPywdRWySC41ct7L0an9KtD5E2cG/UqrjflyqhPPVn8dSj7t0IxYnt0RcfW/j9cgl/Lrtpewbmyje9hBbfR6n7fVt1CSUdW7teHb7KnRAgGlc7m+/mVYUz8DsCBSteezsPr5e9TGfrf+CLbSgwyNLoee5dAsTa0xJPC2HXLhA5I2np+mZTv368Oyz8Msv1o4oxySBC7ND64q4O9PyyHY+XzeNPT61qBe7jwO6FqV7/Yl33XPpBpekTbLmRowY3BPxrn2B0r12UWbAVs418KZXwvdU1if53vMlEhcvMS1E6+8Pr70GW7ZA8kISaV/PKy6GPn+v5X9zX2fpD2Opd+okH/AevRt9jtuTp1GGzO1tnXw9MmRQ3JO3N6xfDzVqmLoKN2+2dkQ5Il0owqzHB81n3YIhnDI8RKOYv3Aqdwvfzn9j8MpcFD9tl4S5PmsXJwUKEoz//q654UJbz/rsXlecw3tu04FfeKHISlrcXYdbYgzayws8PblrhBtxSSQpJ3xib+FpvMtfBDGDwWwKbIGuG8ZDNWOJiU80u1ivdJeI+xIZaWpMnDkDGzfC449bOyJA+sDF/UhI4ODDtSkffp46xv1cq6co3vwoyqAz9SubG5aXk5EkaUehHD9uGpG4cSPs3hpD07sbaMLvuBKPASMGjDiRxB28+N7tOXyfqcm0Cd4kLymaek6ZYSks4soVaNLE9OB+82abmG4vCVzk3JgxMHEiz7GUtdUfp2SHAyhlSojd6/mz9VhEvq2lePcubN9uGpbr5AQuLv9+lSgBnTqZuizNye1kJCEyCQszJfGrV01DZ4cPN3Xx3a+7d2HdOvjuO9Oav1ksvZgdSeAiZzZvRrduzRLXFxlb7itK99nBlZg7khBF4XPhAoweDd9/b2pN9O0LI0dCYDYPwZOS4I8/TEl7xQpT7ZXSpU3ft2yZq1AkgYvsXb2KrlWLc9E+NHQO4fe9XlSpYu2ghLCyf/4xrejzzTemyWzdupla5S4ucOdO+q9z5+CHH+D8efDyMlU8fP55aNECDIbsz5WFfEngSqmhwCuYRmvN1Vp/ca/9JYHbMK2hQwcSNm2lnnEP76+oydNPWzsoIWzI1aswfTrMnGlqVZtjMECbNqak3bmzKYlbQFYJPNfFrJRS1TEl7wZAPLBBKbVOa30y92EKq1mxAjZsYDgzaDlMkrcQmfj5maoXvv22adKPqysUKWJK0ilfRYuCBRY8yam8VCN8BNiltY4BUEr9BnQFPrVEYKIAJSSQ8Pa7HHeqwf5HB7FF7qAQWSta1DRz0wbkZSLPIaCJUqqkUsoT6ACUy7iTUmqAUipEKRUiq+7YqLlzcTl3mrGGyXy3zICLi7UDEkLkRK4TuNb6KPAJsAnYAOwHEs3sN0drHaS1DvL19c11oCKf3L5N/NgJbKMpNd5un9tRTkIIK8jTVHqt9Tda67pa6yZAFCD933Ym6bPPcb1+lSklP+GdUVkvviCEsD15WpFHKeWntb6qlCoPdAMaWSYsUSCuXsU4eQqr6UaPqdEzMvcAABcySURBVI9SpIi1AxJC3I+8Lqm2SilVEkgAXtdaX7dATKKAJIz/CHU3liVVJ7Kyr7WjEULcrzwlcK31E5YKRBSw06dxmvM183iJYV8F4iR1KYWwO/LftpCKGfEecUZn/mo/jiZNrB2NECI3JIEXRocP4xm8jBlOwxg9I+vVcIQQtk0SeCF0+I2PicGDr2p1o++KLQSHhls7JCFELkgCL2TWbwrl4W2rWOjUD+MTUYRHxzL6x4OSxIWwQ5LAC5nw97/ERScwt1oPDB6mFWxiE4xM2XjcypEJIe6XJPDC5M4duv0VzGq6EPFEUrofmV08WAhh0ySBFyLRXyyghPEGX1V+HmfvuHQ/M7cYsRDCtkkCLyyMRhKnTGMnjfinTcl0P/JwMaSuWymEsB+SwAuJW98GU+rGGf5qMoLPXnoYfx8PFKZV22XhXyHskyypVhhozYVyjYgLv0b8geNUrZH7pZ2EEAUvqxV5pAVeCMT8upNy4bvZUmu4JG8hHEhei1kJGxIcGs6Ujce5GB2bbhX58OFTKUEJ6k3vb+0QhRAWJAncQQSHhjP6x4PEJhgBUifouJw+S/vDP/F9wLv0amKZBVaFELZBulAcxJSNx1OTd4rYBCOX3p2HEQPlJ79mpciEEPlFEriDMDcRxy3+Lt1PBbOlWDca95CiVUI4GkngDsLcRJwOf/xNcR2Nfu11lKyWJoTDkQTuIEa2DcTDJc0IE63pf+BnDhtq0GKcrLshhCOSBO4gutTxZ1K3GqkTdJqdvUSt+EOc6/g6rm7S/BbCEckoFAfSpY5/6ozK3Q/3IppiRA9oSePJWzINLRRC2D9J4A7o2qHL1Dm9krUPD+DDXacyDS0EJIkL4QDy1IWilBqulDqslDqklFqmlHK3VGAi9w4OmYsrCSx6oqnZoYVS+1sIx5DrBK6U8geGAEFa6+qAAXjOUoGJ3Ll7K4H//PY1e33bcsDP0+w+UvtbCMeQ14eYzoCHUsoZ8AQu5j0kkRc73w7mwaSLOA95Pcsa31L7WwjHkOsErrUOB6YC54FLwA2t9f8y7qeUGqCUClFKhUREROQ+UpEtrcF78SzCXStSc1SHzEMLkdrfQjiSvHShFAc6AwFAGcBLKfV8xv201nO01kFa6yBfX9/cRyqytXP2QerH/MalLoNQzoZMQwul9rcQjiUvo1BaAf9orSMAlFI/Ao8B31kiMHH/bnw0g1jcqfnFS6nb0g4tFEI4lrz0gZ8HGiqlPJVSCmgJHLVMWOJ+Hfo9imbh33Gi/vO4Plgy+wOEEHYvL33gu4GVwN/AweTXmmOhuMR9OjT8GzyJJeDzN6wdihCigMiSag4g7JyRhIoPY/SvwMNh26wdjhDCwmRJNQf267C1BHCWomOk9S1EYSIJ3M7dvAkBa6ZzzbMcfq90tnY4QogCJAnczgV/fJimxi3EvvAaOEtpGyEKE0ngdiwhAZz+O4M4J3fKjX/Z2uEIIQqYJHA79tPC63S9/S1XWvSCUqWsHY4QooDJZ247pTWcn7AAL2Lw+EQeXgpRGEkL3E5tWGekS/hMLlV+Aqe6ta0djhDCCiSB2yGtYfuwlVTiH0p9NMza4QghrEQSuB1a/4umx+lJnCpSicC9LjSevIXg0HBrhyWEKGCSwO2M1vDjaz9Tm/389/HOJCmn1KXSJIkLUbjIQ0w7EBwazpSNx7kYHYvH5bLMPj+VC+4P8nP1pqn7pCyVJpUHhSg8pAVu44JDwxn940HCo2NJ0lB6TThPsJ25j3Um0ZD+/VeWShOicJEEbuOmbDyeujBx7Bk/RkbNIMK1OMtrt8q0ryyVJkThIgncxqW0qrWGiltjac8G5j/6FHdd3NPtJ0ulCVH4SAK3cSmt6tjTfoyInMlN5yJ8V68DPh4uslSaEIWcPMS0cSPbBjJq1UGK/KZ4mpX8t14PEosU5cOnqknCFqKQkwRu47rU8WfTak/qXhtBnJMr61s9I61tIQQgCdzmXbsGIdOv8qVajGHQQNZ91N3aIQkhbIT0gdu4d96BcTffRHl5ot5/z9rhCCFsiLTAbdiOHXBp/i904BeY8Bn4+Vk7JCGEDcl1C1wpFaiU2pfm66ZSSiorWUhCAgweEM8M5+EkVQ6EwYOtHZIQwsbkugWutT4O1AZQShmAcGC1heIq9KZPhxZHZvAQJ+DLX8DV1dohCSFsjKW6UFoCp7XW5yz0eoXahQsw6/0rHHT+AN2mA6p9e2uHJISwQZZ6iPkcsMzcD5RSA5RSIUqpkIiICAudznFpDUOHwvtxY/AkBvX559YOSQhho/KcwJVSrsBTwApzP9daz9FaB2mtg3x9ffN6Ooc3dSqcW72XfknzUUOHQqBMjxdCmGeJLpT2wN9a6ysWeK1Cbc0aePftRA6VfAMMvvCeDBsUQmTNEgm8J1l0n4icO3gQevWC+Q+MJvDKn/Ddd1CsmLXDEkLYsDwlcKWUJ9AaeNUy4TiutIsylPHxYGTbwNTp8FevwpNPQi/n5fS5MhVeew1697ZyxEIIW6e01gV2sqCgIB0SElJg57MVKYsypNT1BlP51+71/Nl86BqhX9eiyuULhDg/hmu92rB1qwwbFEKkUkrt1VoHZdwuU+kLQNpFGVLEJhj5bscFDiyrgnu4E8HuT3LDxY0N46ZL8hZC5IhMpS8A5pY6S7zpzrWf6hJ/sRg/lmhBueiLPNdzEqf/jOTDv7eY7WoRQoi0JIHn0b36tlOU8fEgPE0Sjz1bkmtr6qATDUypMoC2J35jbOtB7C1bFWITiI5NAEhdbR6QJC6EyES6UPIg7YLDmn8TbnBoeLr9RrYNxMPFgNZwY+fDXF3+KO7ut1lcqStvnfiG5TVa8V2dDmbPkbLavBBCZCQJPA+y6tvOmHCfquXPM6WDiF75GNF/BFL1kYMcLtWK54+t5b8Nn2ZUuzdAqSzPI6vNCyHMkS6UPMgqsaZsv3EDFi6EWbPg5MlS+PnB/Hf+ott3XVHXr7Nn8lcs0VXQ0bH4+3gQE5/I9ZiETK8nq80LIcyRBJ4HGfu2dZIiIcIbj2hf+vSB1avhzh1o1Ag+GBtP9xvf4DJyODz4IOzcSYNatdiR5vWyGm4oq80LIcyRBJ5LSUnQ+5GqfLzoMrfDvYm75EP85WLoBNM/6S1f6NEDhvWOoNau2TDqv3DpErRoAT/8AKVKZXrNlAeV2T0UFUIIkAR+X06dgm+/hW3bIDQUbt0qDZRGGZJw8btJ6fqXeK6jF4O7ehEQtRc1/xvotATi4qBdO1iwAFq3BqesHz10qeMvCVsIkSOSwLMRHQ3Ll8OiRbBzpyn31q8PffpAvXpQt46mquc5XPftgV27YM0umPA3xMeDpye8+CIMGQL/+Y+1L0UI4WAkgWN+LHeryv68/TbMn29qQFetCn2H3OCKcSsPX/iThofO0SYkHL+Th+H6ddMLeXhAUJCpoHfDhtC8ORQvbt2LE0I4rEKfwDM+OAyPjmXYzHPcWl+KqMuuFKl5nuqNwnizyC7Krv6exmf+xqCTiDM4c9IvgDvNOxDQ5glT4q5ZE1xcrHxFQojCotAn8LRjubWGm3sqEf17IAavOBp2XMHLl5fSddFWSsTe5KJ3KWY17MHGKo044VuBBIML/j4e7Hi1hZWvQghRGBX6BJ4yZjvxlhuR62pz91wpPCtfZFD5z/hwwywMSUb+V7khy2u2ZnvF2iQ5GcweL4QQBa3QJ/AyPh6cv2jkypLHMMa4UrbVLqZdGsPTm7ewo0JNhnUaQUSREvc8XgghrKHQJ/BhzQPp97QnxjtuPNF+GXN3juHhyAt80bgn0x97Ll2LWwFpq6fLJBshhDUV6looWsOaWf7EhhVn8KNfsn7Dy/jevcmicV8zu3nfdMnbw8VA74bl8ffxQAH+Ph5M6lZDxmwLIaymULfAJ040Tcz55uU/eXHRKAgKwnPlSl4oU4biOSgTK4QQ1lRol1RbsQKeeQZe73aJGX/WQ3l4wF9/QYms+7uFEMIaslpSrVC2wP/6C/r2hSYN45l+8WnUjRuwcaMkbyGEXclTH7hSykcptVIpdUwpdVQp1chSgeWXuDjo1QseeADWBw7FaddOU42SGjWsHZoQQtyXvLbAvwQ2aK2fVkq5Ap4WiClfffGFqSjVwaHz8Pzya3j7bVNfihBC2JlcJ3ClVFGgCdAfQGsdD8RbJqz8ER4OH34II57YTfWvXjdVBpw40dphCSFEruSlC6USEAEsUEqFKqXmKaW8Mu6klBqglApRSoVERETk4XR59847QEICH4f3gzJlYNkyMBiyPU4IIWxRXhK4M1AX+EprXQe4A4zKuJPWeo7WOkhrHeTr65uH0+XNjh2wZAlMC5qB65njvFK3D43n7s+0ALEQQtiLvCTwMCBMa707+fuVmBK6zTEaYfBgqOAbRbd9H7CrXHU2Pdwgy1XkhRDCHuQ6gWutLwMXlFIpc8lbAkcsEpWFzZsH+/bBO76jKBlzg4+bv5S6Cry5VeTzKjg0nMaTtxAwah2NJ2+RNwghRL7I6yiUN4AlySNQzgAv5D0ky4qKgjFj4OlHL/DC3kUEV23KwQcrp9vHkhUFzdUXH/3jQQCZySmEsKg8jQPXWu9L7t+uqbXuorW+bqnALOWDD0wL5sz2G4tCM7VJ30z7WLKiYNr64inyo5UvhBAOXczq2jWYMwfe6xRKibXfcq73y0SWKpNuH0tXFMyqNS91w4UQlubQCXzWLIiN1Yy8MgJKlKDKFxOZ1K1GvlYUzKo1L3XDhRCW5rC1UGJjYeZMGFd/PV67t8CXX4KPD13q+ORrX/TItoHp+sBB6oYLIfKHwybwRYtMXShv+k6ASpVg4MACOW/Km4OUohVC5DeHTOBGI3z2GfR/ZDdFj+4xNcVdXQvs/F3q+EvCFkLkO4dM4D/9ZCpY9evjM6BoUVPtWCGEcDAO9xBTa5gyBRqUu0T53cvhhRfA29vaYQkhhMU5XAt8507YtQt2t5+NCks0zaEXQggH5DAJPDh5DcvQ+VXxdC9Cjd3/hQ4d4OGHrR2aEELkC4dI4CnT129edif25AO8UvkzPE5GsDDoSeZO3iKjQYQQDskhEnjK9PWbIQFgSGLQzXmcKlGWCbH+6LumGZBSk0QI4Wgc4iHmxehYkuIN3DlShqYVfqHOleMsqtcJnVxxMIXUJBFCOBKHSOBlfDyIOf4gOt6FN4wzuOnqyY/VWpjdV2qSCCEchUMk8JFtA4k5WI6yPqd46sIWVtRsTYyb+fWVpSaJEMJROEQCf8TTn9gLJXiz5HQMSUlsbNqd3g3L4+GSfr1LqUkihHAkDvEQc/58cHZKYkjczzi1ac3yqX0ACKpQQmqSCCEcll0m8JQx3xejYynt7cmxb5owsuF2DDvPweSPU/eTmiRCCEdmdwk845Jlp//2JjrSQFfjXChSBLp0sXKEQghRMOyuDzzjkmW3D5TDy/M6j+wLhu7dwcvLitEJIUTBsbsEnnYYYOJtN2JP+/K0/7cUibsDffpYMTIhhChYeepCUUqdBW4BRiBRax1kiaDupYyPB+HJSfzOobKgnegdt4yrRUvh16xZfp9eCCFshiVa4M211rULInmDacy3h4sBrU3dJ2XLHKNZ2F9Ed+0BBkP2LyCEEA7C7h5ipowqGfvVFRKvezGgwmJcLhqpMuI1K0cmhBAFK68tcA38Tym1Vyk1wNwOSqkBSqkQpVRIREREHk9n0qWOP0HxdSlaFN512gR16kD16hZ5bSGEsBd5TeCNtdZ1gfbA60qpJhl30FrP0VoHaa2DfH1983g6k5gYWLUKBrc6huHvEHl4KYQolPKUwLXWF5P/vAqsBhpYIqjsrFkDt2/DK+7fgpMT9OxZEKcVQgibkusErpTyUkp5p/wdaAMcslRg97JkCfg/mESF7d9BmzZQunRBnFYIIWxKXlrgDwDblVL7gT3AOq31BsuElbXISFi/HsY88Tvq/HlZcV4IUWjlehSK1voMUMuCseTIihWQmAg9EpaYps537lzQIQghhE2wu5mYS5ZAjf8kUPK3H03J29N83W8hhHB0dpXAz52D7dth9KNbUFFR8Mwz1g5JCCGsxq4S+NKlpj87xa4Ab2/TA0whhCik7CaBa23qPnmiYQLev642dZ+4u1s7LCGEsBq7SeAHDsDhwzCy3haIioIePawdkhBCWJXdJPAlS8DZGVpfXw5Fi0r3iRCi0LOLBJ6UBMuWQYfWCbivXw1PPSXdJ0KIQs8uEvjvv0NYGAyvuRmuX5fRJ0IIgZ0k8CVLTCulNb60QrpPhBAimV3UAx89Gjq1icfl1eTRJ25u1g5JCCGszi4SeKVKUOl4cveJjD4RQgjATrpQAFMRFOk+EUKIVDafwINDw2n60UZuLF3O+oceJfjINWuHJIQQNsGmu1CCQ8MZ/eNBHj22m2Jxd1hRqRF//ngQ+HdtTCGEKKxsugU+ZeNxYhOMdDy2nZuunmyvWIfYBCNTNh63dmhCCGF1Nt0CvxgdC8CnTfvyU9WmxDu7pNsuhBCFmU0n8DI+HoRHxxJRpAQRRUqk2y6EEIWdTXehjGwbiIeLId02DxcDI9sGWikiIYSwHTbdAk95UDll43EuRsdSxseDkW0D5QGmEEJg4wkcTElcErYQQmSW5y4UpZRBKRWqlFpriYCEEELkjCX6wIcCRy3wOkIIIe5DnhK4Uqos0BGYZ5lwhBBC5FReW+BfAG8DSVntoJQaoJQKUUqFRERE5PF0QgghUuQ6gSulOgFXtdZ777Wf1nqO1jpIax3k6+ub29MJIYTIQGmtc3egUpOAPkAi4A4UBX7UWj9/j2MigHO5OiGUAgpbJSu55sJBrrlwyMs1V9BaZ2oB5zqBp3sRpZoBI7TWnfL8YlmfI0RrHZRfr2+L5JoLB7nmwiE/rtmmZ2IKIYTImkUm8mittwHbLPFaQgghcsaeWuBzrB2AFcg1Fw5yzYWDxa/ZIn3gQgghCp49tcCFEEKkIQlcCCHslF0kcKVUO6XUcaXUKaXUKGvHUxCUUmeVUgeVUvuUUiHWjic/KKXmK6WuKqUOpdlWQim1SSl1MvnP4taM0dKyuObxSqnw5Hu9TynVwZoxWpJSqpxSaqtS6qhS6rBSamjydoe9z/e4ZovfZ5vvA1dKGYATQGsgDPgL6Km1PmLVwPKZUuosEKS1dtjJDkqpJsBtYLHWunrytk+BKK315OQ36+Ja63esGaclZXHN44HbWuup1owtPyilHgQe1Fr/rZTyBvYCXYD+OOh9vsc1P4OF77M9tMAbAKe01me01vHA90BnK8ckLEBr/TsQlWFzZ2BR8t8XYfrFdxhZXLPD0lpf0lr/nfz3W5gql/rjwPf5HtdscfaQwP2BC2m+DyOf/jFsjAb+p5Taq5QaYO1gCtADWutLYPqPAPhZOZ6CMlgpdSC5i8VhuhPSUkpVBOoAuykk9znDNYOF77M9JHBlZptt9/tYRmOtdV2gPfB68kdv4Zi+Ah4CagOXgM+sG47lKaWKAKuAYVrrm9aOpyCYuWaL32d7SOBhQLk035cFLloplgKjtb6Y/OdVYDWmrqTC4EpyH2JKX+JVK8eT77TWV7TWRq11EjAXB7vXSikXTIlsidb6x+TNDn2fzV1zftxne0jgfwGVlVIBSilX4DngZyvHlK+UUl7JDz9QSnkBbYBD9z7KYfwM9Ev+ez/gJyvGUiBSElmyrjjQvVZKKeAb4KjW+vM0P3LY+5zVNefHfbb5USgAycNtvgAMwHyt9cdWDilfKaUqYWp1g6lezVJHvGal1DKgGaYym1eAcUAwsBwoD5wHemitHeahXxbX3AzTx2oNnAVeTekftndKqceBP4CD/Lvwy7uY+oQd8j7f45p7YuH7bBcJXAghRGb20IUihBDCDEngQghhpySBCyGEnZIELoQQdkoSuBBC2ClJ4EIIYackgQshhJ36PzziAAdF5eA1AAAAAElFTkSuQmCC\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.842590\n",
       "x1                  0.525808\n",
       "np.sin(x1)          0.533273\n",
       "I((x1 - 5) ** 2)   -0.021888\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.894666\n",
       "1    10.727418\n",
       "2    10.439990\n",
       "3    10.080038\n",
       "4     9.710299\n",
       "5     9.393223\n",
       "6     9.175688\n",
       "7     9.077523\n",
       "8     9.086650\n",
       "9     9.162042\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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}