{
 "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.985\n",
      "Model:                            OLS   Adj. R-squared:                  0.984\n",
      "Method:                 Least Squares   F-statistic:                     1027.\n",
      "Date:                Sun, 09 Aug 2020   Prob (F-statistic):           3.93e-42\n",
      "Time:                        21:12:11   Log-Likelihood:                 2.9050\n",
      "No. Observations:                  50   AIC:                             2.190\n",
      "Df Residuals:                      46   BIC:                             9.838\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.9967      0.081     61.588      0.000       4.833       5.160\n",
      "x1             0.4902      0.013     39.174      0.000       0.465       0.515\n",
      "x2             0.5709      0.049     11.606      0.000       0.472       0.670\n",
      "x3            -0.0183      0.001    -16.678      0.000      -0.021      -0.016\n",
      "==============================================================================\n",
      "Omnibus:                        2.560   Durbin-Watson:                   2.102\n",
      "Prob(Omnibus):                  0.278   Jarque-Bera (JB):                1.454\n",
      "Skew:                          -0.055   Prob(JB):                        0.483\n",
      "Kurtosis:                       2.172   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.5386542   5.03704758  5.49210698  5.87271985  6.15900199  6.3455645\n",
      "  6.44239911  6.47323659  6.47164785  6.47552834  6.52087193  6.6358572\n",
      "  6.83621783  7.12265755  7.48073452  7.88323413  8.29464025  8.67697041\n",
      "  8.99601627  9.22696389  9.35847064  9.39452873  9.35380932  9.2665949\n",
      "  9.16980393  9.1009253   9.09185897  9.16367489  9.32315131  9.56166109\n",
      "  9.85658913 10.17504798 10.47928113 10.73286632 10.90669959 10.98377701\n",
      " 10.96198836 10.85446383 10.68741675 10.49583687 10.31774159 10.18792983\n",
      " 10.13226465 10.16342406 10.27881801 10.46101396 10.680601   10.90102196\n",
      " 11.08457965 11.19863041]\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.21218748 11.07725494 10.81622059 10.48010344 10.13606243  9.85095361\n",
      "  9.67496145  9.62931189  9.70107555  9.84633324]\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+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3WUlEQVR4nO3deZyN1R/A8c+ZOztjxmiIGSFbyb5labHEJJUtFSXKElJCihbRqpBESFJKsu/JTvwoGmYsYVChGdtYZmyzz/n98cxMY9xrhnvv3GW+79fLy8xzn3uf7+OO75x7zveco7TWCCGEcD0ejg5ACCHErZEELoQQLkoSuBBCuChJ4EII4aIkgQshhIvyLMiL3Xbbbbp8+fIFeUkhhHB5O3fuPKu1Dsl9vEATePny5YmIiCjISwohhMtTSh0zd1y6UIQQwkVJAhdCCBclCVwIIVxUgfaBm5OamkpMTAxJSUmODkXkg6+vL2FhYXh5eTk6FCEKPYcn8JiYGAICAihfvjxKKUeHI25Aa825c+eIiYmhQoUKjg5HiELP4V0oSUlJlChRQpK3C1BKUaJECfm0JISTcHgCByR5uxB5r4RwHg7vQhFCFDIXLkBsLMTHG19n/X3pEjz2GNSs6egIXYYkcMBkMlGjRg1SU1Px9PSke/fuvPrqq3h4WP6AcvToUbZt20bXrl0LMFIhXFhqKnz8MXzwgfG1OSNHwuuvwzvvgK9vgYbnilwugS+JjGXM6mhOxCdSJsiPoeFVaV8n1KrX9PPzIyoqCoAzZ87QtWtXEhISGDVqlMXnHD16lNmzZ0sCFyI/9u6F7t0hMpK1NZuztHxDPG8Lpn3ze2jWqCoULw7p6Uby/ugjWLgQvvkGmjZ1dOROzSn6wPNrSWQswxftJTY+EQ3ExicyfNFelkTG2uwaJUuWZNq0aUyaNAmtNUePHuX++++nbt261K1bl23btgEwbNgwtmzZQu3atRk/frzF84Qo1NLS4MMPoV49ko79y8tPvE3vNkNYcff9LAm5h37RJpZc8oMSJaBkSfjuO1i1CpKS4P774eWXja4VYZbKa0s1pdQM4FHgjNa6euaxzsBI4G6godY6Xwuc1K9fX+deC+XAgQPcfffd+Qq26egNxMYnXnc8NMiPrcNa5Os1zClatCiXL1++5ljx4sU5ePAgAQEBeHh44Ovry+HDh+nSpQsRERFs2rSJsWPHsmLFCgCuXr1q9jx3dDPvWWFl7pMiYPNPj87C3P16xV+g/HOduSfmIOtqNuPDh/vzj/K/7rlm//9evgxvvw1ffAFhYTBnDjRpUkB343yUUju11vVzH89PF8p3wCTg+xzH9gEdga9sEl0+nTCTvG903BpZv9hSU1MZMGAAUVFRmEwmDh06ZPb8/J4n3F/WJ8XE1HTA+KQ4dP5uUJCarrOPDV+0F8Dlk7i5+x056zdmznmLSqf/oX+7Yay86z6Lzzf7/7doUfj8c3jqKejWzRjcjIgAmX9wjTwTuNZ6s1KqfK5jB6DgS8rKBPmZbYGXCfKz6XX+/vtvTCYTJUuWZNSoUZQqVYrdu3eTkZGBr4WBlfHjx+frPOH+xqyOzk5mWVIzrv+km5iazpjV0dclcHuM89hT7vv1Tktl0sIPuOfkEfp2eIt1le+94fNv+P+3cWOjS6VBA+jYEbZuBf/rW/GFld0HMZVSfYA+AHfccYdVrzU0vOo1v+kB/LxM2R9PbSEuLo6+ffsyYMAAlFIkJCQQFhaGh4cHM2fOJD3duHZAQACXcvTNWTpPFD4n4hOpcfIwPSOWUOH8CWKLhfBv0O3EBJbk38BS7L29MueKBGWfm5O51qyzt9Rz3oMpI50Jy8dw37HdDGo72GzyTr/ijYdPGsozI3//fytVgp9+gkcegV694McfQeYjAAWQwLXW04BpYPSBW/NaWT/Atm6dJCYmUrt27ewywm7dujF48GAA+vfvT6dOnZg/fz7NmzenSJEiANSsWRNPT09q1apFjx49LJ4nCpGMDPjlFxbOf4u6f+/morc/UWWqUuXscVr8HYFvWgoAV7x8GXd/N76r9yilg4te8xLmWu+WWurOIvuTsdZ8tGoSbQ5tY1TL3iyufm2/dlqCH0k7K3EhMpTiLf+karOz+f//+/DDxmDom29CvXowZIid7sa15DmICZDZhbIiaxAzx/FNwGsFNYgpnIO8Z2asXg2DBsGBA1wtVYZx5Tox5WofLlwqhWdgIj7FrlLaL4Y79WFeOfItLY/+wb7SlYkbN5HmXcKzX6bCsJ8x9z9SAf+Mbltgt3Mzsj41vLrma17csYgJTbow6cFns/v8U84W5eL2ilz5swyeJkX37oo33oDKlW/yQlpD586weLHx7/3QQ3a5H2dkzSCmEOJGVqyAjh1Jv7MS67v/yJuRndm5wwsPr3Q8b7tE8j/FuHzJh3PcxT4eYrlXT/rVnsq44x9QvVtb2PkqjBoFRYoU2DiPLbWvE0rZ5fOpt2MR39dty7xHezLm4bu4lKAY9IqJc3tL4eGVzmNdrjL5k6KEhd3ihZQyygwPHjQGN2VQM+86cKXUT8BvQFWlVIxSqqdSqoNSKgZoDPyslFpt70CFcAZLImNpOnoDFYb9TNPRG9g2aRa6UydiQ2oRenQb4TO7gpcXU6bA+TgTySeCSLnoQ2IiREcb43GPPuzN5KhXaBwQzT8tXoBx4+Dee+H0aYaGV8XPy3TNNW09zmNze/ZQ7+Ph0Lw5z21fwtbhLalWJJT3+5ThUnQp3noLTsWaWPajFck7S9GisGSJ0V3VsSNcvWqLO3BZeSZwrXUXrXVprbWX1jpMa/2N1npx5tc+WutSWuvwvF5HCFeXeyJZ+ajfqPfqCxzxvZvqJ9bQrH0Qu3YZDcO+fSEw8L/n+vpClSoQHg7LlsGaNZBatDh3rp1Gh9sXcDX6L47WaIhP3Gk+7liD0CA/FEaN9Mcdazht/zcJCdCpEwQHGwONnp5s2ACNGhnLm6xfb8ycD7luO14rVKoEs2fD7t3GRJ9CzKVmYgrhSDkHGBsd38P0he8TbapMo4vrCXjoAr+X+5kBqzfka2Zwq1Yw6ttYSrX5k2UJj/Gw/oXbLpzl7i6P43v6BFuHteCf0W3ZOqyF8yZvraFHDzh6FObNg1KlmDbN+CVVujRs3w73WS7/tk6bNsa0+xkz4Ndf7XQR5ycJXAgLcneXZPVN1405wIz5o/ibCrRM2QjhsXjUOwDq5pZ3+Gx9NL41j1Km569sD6lN67Q1lEiIp8Yz7eD4cXvfnvXGjjW6M8aOJaNxU4YMgRdfNMYWt22DO++08/VHjDD6wPv2heRkO1/MOUkCF8IMc+vuKCAg+QoTl44hNiOMlno9quNxAmr/e81zs8r+8pJVP20qkkKpLr8TdWdlHkpfT9GEy+gHHzRatgUo9y+sJZGxZo8BsGkTDBsGTz6JfvkVBg6Ezz6DAQNg+fJru4/sxt8fvvzSGNQcM6YALuh8Cn0Vyrlz52jZsiUAp06dwmQyEZLZYbdjxw68vb0dGZ5wEHP12Bp4a90Mbr98lqZei+HJo/iHXTD7/Pws75Cz4sTDO52QThHsX12dlns2sOlEa4q0aIHatg1uv93q+8nLzUz/9407xcPPPQ1VqqC/ns4bwxSTJsHgwUajvEDn2LRpY5QWfvABPP200T9eiBT6FniJEiWIiooiKiqKvn37MmjQoOzvvb29SUtLc3SIwgHMJeAH/4rg6X2r+YTXOdrRg4rVkyjub35z5/yU/eWuOFEemtDH9lO5bxVapKwi+fhpMh5pWyCr8Vma/p+VvLOPJSVTuk8PY7GphQsZ9VkAY8ZAv34OSN5ZPv8cfHygf3+jX74QKfQtcHN69OhBcHAwkZGR1K1bl4CAAIoWLcprr70GQPXq1VmxYgXly5dn1qxZfPHFF6SkpHDvvfcyefJkTCZTHlcQzi53PXaxpMuMXv4l+7iHjLdHcup9H+D6livkv+zP8sziYkyr05BOL85jWVQ7Mjp1xuPn5eBl/peFLeR3Qbjhm76l1rF9MGcOnyyvxqhR8PzzMGmSA2e3lyljzNJ8+WVj1cIuXRwUSMFzqgT+6quQua+CzdSubfyCvlmHDh1i3bp1mEwmRo4cafacAwcOMHfuXLZu3YqXlxf9+/fnxx9/5LnnnrMmZOEEcq+789aK7whJPs/g6kuYNdIn+zxrl3doXyfU7Ll9+kBGRlte7DeV6Wt7k9GzNx4zv7VblrQ0gSintge20DNiKfOadOTU6acYNszIlV9/DTfYvKpg9OsHM2cas2HbtIGgIAcHVDCcKoE7k86dO+fZkl6/fj07d+6kQYMGgLGmSsmSJQsiPGFnORNz1d+28NRfqxjj/xYT1jch94+FpSRsrb59YVJaL0a+HMPIH0aRERaGx0cf2Pw6YH6hOC8Pld0HXvHsv3z6ywR2hVVjy4OTmTQQOnQwcqZTfOA0meCrr4xVC4cPhylTHB1RgXCqBH4rLWV7ybkYlaenJxkZGdnfJyUlAcaa4d27d+fjjz8u8PiE/bWvE0q7O/yIH/UMu6lJkxUjKOjfzwMGwPiUd5k+JIZeH39IeplQTAP62fw6lj5JAHy5LJIpiz8k2ceXpV3m8uXoUrRpY/RW2LFX5+bVrWt0o3zxhbF9W6NGjo7I7hz9wccllC9fnl27dgGwa9cu/vnnHwBatmzJggULOHPmDADnz5/n2LFjDotT2N7Bx1+naNJZol6dSdPmjqlIGjRYEf/xVJbzKOrll0ifOcsu12lfJ/S6CUTta5Vm7Z/fUyn+BHuHLOKDcdVp1szYstIpC7Tef9/oE+/b19jOzc1JAs+HTp06cf78eWrXrs2UKVOoUqUKANWqVeODDz6gdevW1KxZk1atWnHy5EkHRyts5dTq3VTdNoNl5V6h27jaDo3ltWGeHHp/HhtpDj26kzpnoVWvZ7G+OyetjYGp+fM50OMTHvqwOY0aGUsB+Dnr2loBAcZH+d27jRpxN5ev5WRtRZaTdQ+F4T3TGZrdpVpR9mwklyOPUK52cUeHBMCUMZep+Xo4DdUfZCxcgk+HR655PD+7+ViqnLluzZV334X33uNYp8FUXT6W6jUU69fDxr+dfMcgrY31w3/7zVhBrHRpR0dkNUvLyUoLXAgztr29ktpn1xPV7l2nSd4A/YYW5a8vVrJH14AnOnH1543Zj5mbPWpuWv+NNo3INn48vPce/zR/gbtWjKVyFcXq1Ubyzs81HEopo64xORkyS3/dlSRwIXJJOJdGyJihHPWpwoM/2X7A0FrPvRzI0amrOZxRER5/jEu//A/IZ2ImH5uDz5gBgwdzsMYTVNo4jTp1FRs2QIkS+b+Gw1WuDG+8YaxauHFj3ue7KEngQuSyssPXVEk7QNqHn+Lp5xxlFrn7rE0Nkzn+zTpiMkLxeaQFx4d9yYkL5tfGzp2wLc0SLRPkBwsWoHv3Zl9YOLX2zqLzUyY2bPhvOdg8k78zGT7cWOzqpZcgJcXR0diFJHAhcti+JoGHtozgSOiDVBr8uKPDASx3jWwLjqNHt6msNbXkjk8GMH7mVPyTrk/iuRO2uU0jgkhj5p9z0U8+yf6ARtwbs5Chb/kwe7axlrml18rruEP5+cHEiXDggNEl5IYkgQuRKSUF9j7zMSU4R5k5nznNzueWui1+/P04J8pcpV//gYy4bTiPnf6FxV++ScWjJ7LPMzetv32d0Gs2jXjo0jE2z32NSrOmMTugL00vr+bLb4vwwQfXz7B0uR2D2raFdu3gvfdcY4nemyQJXIhM098+Srez44lt3g3/++o6OpxslronsurHPPzTmflCUzo0mEZg2iWWzx1InwW/UP9igsXdfNrXCWXroKb87bGNr6e+zNXYJFqxho/CJrN0fVF69DAfS+7k7/Q7BgFMmGBUpgwa5OhIbE4SOBATE0O7du2oXLkyFStWZODAgaRk9plt2rSJRx999LrnrFixgjp16lCrVi2qVavGV199Zfc4e/TowYIFCwDo1asX+/fvt3jupk2b2LZtW/b3U6dO5fvvv7d7jK7q+HEI/uwtMJko+/2Hjg7nGvnpnlAK9rS4nTZdJxMV2Iw3//qSBVOeodYDbYjqMIpLv+2Dv/4idd4iEl59lwvNOpBY+k7Uxx/xbUZ3wkP30f2HVuzZAw8+eONrmZvw49TKlYN33oFFi+CXXxwdjU0V+gSutaZjx460b9+ew4cPc+jQIS5fvsxbb71l8Tmpqan06dOH5cuXs3v3biIjI2nWrNktXT89PT3vk8yYPn061apVs/h47gTet29fWWTrBiY8H8XT6bNJenEg1u+8a1vmui0sde741wii8YWV/LnqX5a1nMCZ1CBqLhlFQJMaUKkSXk91ouiEDzj96wGWX2jKcyV+JnXqDCIOBfLss06yrok9DBkCd91lLDmbkODoaGwmP7vSz1BKnVFK7ctxLFgptVYpdTjzb+cplL1JGzZswNfXl+effx4Ak8nE+PHjmTFjBlct7Hh96dIl0tLSKFGiBAA+Pj5UrXp9H+DIkSPp1q0bLVq0oHLlynz99deAkVybN29O165dqVGjBunp6QwdOpQGDRpQs2bN7Na81poBAwZQrVo12rZtmz1lH6BZs2ZkTYpatWoVdevWpVatWrRs2ZKjR48ydepUxo8fT+3atdmyZQsjR45k7NixAERFRdGoUSNq1qxJhw4duHDhQvZrvvHGGzRs2JAqVaqwZcsWW/wTO71Vq6DlhjdJ9CtO4IdvODqc65jrtnim0R0W+6KVgnvCw3h83Ss0TNzMlAl7GFppNAPKjubRBosZ/PpxDi87SIUd85gW8wgvvuhka5rYg7c3fPst/Puvse+bm6wbnp/FrL4DJgE5P38PA9ZrrUcrpYZlfm/9T74D1pP9888/qVev3jXHihUrxh133MGRI0fMPic4OJjHH3+ccuXK0bJlSx599FG6dOmCh5k1Nffs2cPvv//OlStXqFOnDm3btgWM3X727dtHhQoVmDZtGoGBgfzxxx8kJyfTtGlTWrduTWRkJNHR0ezdu5fTp09TrVo1XnjhhWtePy4ujt69e7N582YqVKjA+fPnCQ4Opm/fvtesYb5+/frs5zz33HNMnDiRBx98kBEjRjBq1Cg+z/w3SktLY8eOHaxcuZJRo0axbt26vP6FXVpyMvzQ61d+5BfS3vnUaZchNbfiYf1ywXnOiFwaFcukuBgSO1XPPvaX116ah0EDZ+/6sLVGjYyde4YPN3aV7tnT0RFZLc8ErrXerJQqn+twO6BZ5tczgU3YIoE7gNYaZabawNLxLNOnT2fv3r2sW7eOsWPHsnbtWr777rvrzmvXrh1+fn74+fnRvHlzduzYQVBQEA0bNqRChQoArFmzhj179mT3byckJHD48GE2b95Mly5dMJlMlClThhYtWlz3+r///jsPPPBA9msFBwff8H4TEhKIj4/nwcyOzu7du9O5c+fsxzt27AhAvXr1OFrAezI6wrixmpdj3yDxtjD8Xh3g6HBuSn6Wsb3RxBun77u2h9dfh/XrjVULGzeGG3RDuoJbXU62lNb6JIDW+qRSyjaLbDpgPdl77rmHhQuvXRjo4sWL/Pvvv1SsWJFz585ZfG6NGjWoUaMG3bp1o0KFCmYTeO5fAlnf51yuVmvNxIkTCQ8Pv+bclStX3vCXSNZz8zrnZvj4GJsVmEwmt99O7tgx2PPeEt5kO4ye7sQrNN06l5p4UxA8POD776FWLWMPze3bXfp9t/sgplKqj1IqQikVERcXZ+/L3bSWLVty9erV7AqN9PR0hgwZQo8ePfD39zf7nMuXL7Np06bs76OioihXrpzZc5cuXUpSUhLnzp1j06ZN2Zs/5BQeHs6UKVNITU0FjN2Arly5wgMPPMCcOXNIT0/n5MmTbDQzJbhx48b8+uuv2Uvcnj9/HoCAgAAumdlLMTAwkOLFi2f3b//www/ZrfHCImtWY41WMbyb8iYXwioZ60e7IZeaeFNQSpc2kvjevcbgpgu71QR+WilVGiDz7zOWTtRaT9Na19da18/a7d2ZKKVYvHgx8+fPp3LlylSpUgVfX18++uij7HPWr19PWFhY9p/IyEg+/fRTqlatSu3atXn33XfNtr4BGjZsSNu2bWnUqBHvvPMOZcqUue6cXr16Ua1aNerWrUv16tV58cUXSUtLo0OHDlSuXJkaNWrQr18/s4k2JCSEadOm0bFjR2rVqsVTTz0FwGOPPcbixYuzBzFzmjlzJkOHDqVmzZpERUUxYsQIK/4FXUvWrMYjO4vS+fBq7uYgIxo9zZK9px0dml243MSbgvLww8ZCV1OmGOWFLipfy8lm9oGv0FpXz/x+DHAuxyBmsNb69bxep7AtJzty5MhrBhLdhSu/Z01Hb+DfMymcn96Ag5drcub2YnTsNobQ4v5sHXb9GIM7yM8Ss4VSSgrcdx8cPmwUT1j4FO0MLC0nm2cfuFLqJ4wBy9uUUjHAu8BoYJ5SqidwHOhs+RWEcJzcySs2PpH4zdXof+kbwjjBkGYvg1LExifSdPQGt0xy9tqz0+V5e8NPP0GdOtC5s1FPmkcRwC05eBBeecVY5dHGcwzyU4XSxcJDLW0aiRuytJu9KBi5Ny6IjU8k5UQQ/jt9edvjfdZVaMD2O2oAxsSYrF3ZsxaLAiTxubuKFeGHH+DJJ40ywxUrIHPHLatpbewKNHQo+PvDX3/ZPIE7xUzMgtwVSFjHld6r3CV0Ol1x9pcafOQ1DF+S+KBFL8BI3rnvyinXuBb20a4dbNgAFy4YSXzDButf88QJaNPGKFds3hz27ct7jYJb4PAE7uvry7lz51wqMRRWWmvOnTuHb871Rcnn/ooOkLtULuH3itQ8e4geqT+woGlHjgWHEhrkd13ytvR84caaNoUdO4wKlfBwmD791l9rwQKoUQM2b4bJk+Hnn+22rdut1oHbTFhYGDExMThjiaG4nq+vL2E5Pgaa66Zwlu6HrD5vgJSzRUnYWomJfo2INwXSdfk0ugYGAsbAZqyZZF2oS+0KowoVYNs2eOop6N3b6Lv+5JP8LxBz6JAx0/OHH6BBA5g1y3bdMRY4PIF7eXllzyIUrsfSTL8h83YzaG6UQwcEh4ZXZfiivVxNTuf8LzXp6jmbxokRRI4YQ3Bm8s55Xu5Nfgt9qV1hFBho9IMPGgTjxsHOnfDEE0bXSs2a1y8ac+wYzJ0Lc+ZAZKSR7EeOhDffLJAFZhy+K71wLeaqOvJidsfzArIkMpZX37pK3C9lOOJXCb9yIQTti7yuVSWlduI6U6YYG0GcOmV87+sL9esbybxkSVi8GH77zXjs3nuNmZ2dO0Oo7X9uLJURSgIX+Za7uwTMDwCaExrk55A6602b4KGHYHbld3ny4HuwZYtR+ytEfmhtLBb/++///dm1y6ghz5qO/+STcOeddg3jluvAhchirrtEk78k7ogBwZgYozvzgfLH6Xz0U+MbSd7iZihlTPApV874+QFjCcuzZ+3S0r5ZDq9CEa7jRlt7Za1VbbKwsFZBDwgmJxtdl4lXMlhWqpex4NennxZoDMJN+fg4RfIGaYGLm2Cpzztn94i5bhZrBwRvpX/61VeNheZ2d/+cojPXwtSpcMcdtxyDEM5IWuAi3/KzMJKtN73N+oUQG5+I5r8yxRvVmn/7rZGvP+8RRc2fhhsTNfr0uaXrC+HMZBBT3JSCrtawVKNtaVB02zZo0QIeanKV5SfroRISYM8euO02u8UohL3JIKawiZwLI12+DFsWxTHxzZUkRR3k1G33cOGuJoTcGUDZslC+vDGLOMfeFTftZjYkWLwYunaFsmVhfrkhqI0HYc0aSd7CbUkCFzcl5WISv/eaTvqWrZQ7tYM2/P3fg6cgfZ8HUaoOm/X9fE0znvV/hMc6evHss9CyJXje5E+cpX733IOin38Ogwcb5bir+i/D77mpxmL9rVrdwl0K4RqkD1zk2+ntR/mrzH08MP9lqp7dyqXKdTnS51NS126C06dh7VpM77xF3QcDeNV3KktpzyHPuymy6AceeTidsDBjgtv+/fm/Zl797unpMHCg8bodO8KGWScIHNzT2Mz6ww9td/NCOCOtdYH9qVevnhauae8nP+vzqri+QKDe8tqSvJ+QnKz1kiVa166tNeiLoXfpsQ3naB+vdA1aN2+u9fz5Wqek5P1Si3fF6CYfr9fl31ihm3y8Xi/eFaO11vrCBa3bt9catB48WOv0f45pXaWK1v7+Wu/fb90NC+FEgAhtJqdKAhc3lJGapreHv6M16P3etXT0yiM39wLp6VovWKD1PfdoDTr17hp67gurdLlyxk9fmTJajxih9c6dxqn5sX+/1v36aV2kiNZKaf3FF1rr6Gity5bVulgxrbdsudnbFMKpWUrgUoUiLEpNuMrBuztQ4+Qa1pV9ngY7viTw9lubkLMk4jiRY6fx/OoZlI8/yakmzZnzwIeMWlSF+EMlAAgMTuexR0yEh8MDDxj95amp//05dMhYH3/tWmMuRZcuRvdJbaKgdWvjQqtXGzusCOFGZC0UcdO21+hFg30z+Pmxr2i7pDcetzhiknNyj3daKs/tWs7AbXPwS0liVp1HGFenBydPVST1aEmILcXFeMvLd4aGQv/+xmqfISHA1q3Qti0UK2Zk9qqygqBwP5LAxU35Y/BPNBjfldX13iQ8wrrBQHO13MFXExj0vx/pGrWKK95+LKjekh9rtyHxzip80boFf/wBHh7GipxZf4KDjbJELy9YuWEvxz79gu4bfuBsYAgHZi4k/JGGVsUphLOSBC7y7dj6IxR/qC5/F63J3ac24VPEumrTCsN+trjYVZW4o7yybS6tD/2Gd0Yav5etTqNP3jRKSnx8rn9CRATH3hvD7SsX45OeyubydRj86GCuBN3msCVrhbA3SeAi241mUybGJ/NPmaaUTvqbK/+LIqyJ9euHWJpNmVOJK/F03ruObntWE3rhJBQtavSRFCtm/AkIgDNnICKCq96+LLynBd/XacvhkHLZrxHk50URH09Z01u4HbvMxFRKDQR6Y6wo+rXW+nNrXk/YX15boG29fxgPJe4k4q3F1LdB8gbzO954eShQkJpuNCDOFQli5v1Pcde49wg9dwCWL4f4eLh4ES5dMurMTSaYMIFG/9zORZ/rp3fGJ6YSn5hq9r6EcEe3nMCVUtUxkndDIAVYpZT6WWt92FbBCduztAXamNXRBMzcxUP7Pmdb/Zdp8kF7m10zK4HmbvWbO2acW/a/qhIzAkZv4GI+1hfPui9J4MJdWdMCvxv4XWt9FUAp9SvQAZBFl52YpbVFLh2Ip/b3vTlUtA4NN42x+XVzrqGS+3hO+Vksy1yL3hLZWV64M2um0u8DHlBKlVBK+QOPAGVzn6SU6qOUilBKRcjO845naWOF/itWEkgCfdsO4MGJW2+4XKu95HfpWHNL1hb3N7+BrOwsL9yZVYOYSqmewEvAZWA/kKi1HmTpfBnEdDxzGy6Ujk7i1yVdmF2qPSN7dAeM9UY61Qtl48E4p106NidLG0lIZYpwB3YZxNRafwN8k3mBj4AYa15P2F/u/ujSxfx4afUPpGPiy8cfyz4vMTWdH38/nl3+VxCDgjezdGxulvrZJXkLd2ZtFUpJrfUZpdQdQEegsW3CEvaUsz96+Yd7aJs4j0l3dicuOPia83J/NrP3oGB+l461xFI/uxDuytrlZBcqpfYDy4GXtNYXbBCTKCBXroDf+2+SoIrx9aMP5+s59hwUzM+WbUKI/1jbhXK/rQIRBW/+K1vokfwzGzu9TVqxQMjRf6y4vgUO9h0UlG4QIW6OzMQspGJjNP+Wa0oV72MEnzvMkugL1yTO5neFsHBnrAwKCuEEZE9McY35zy3n1YzfOPvOV+DvT/s6/tcl5vrlgqU1LIQTkxZ4IfTnnnR0rVqUKp5CyOk/jeX9hBBOy1ILXPbELITWvryM6vyJ36fvSfIWwoVJAi9k/v4b6m/+jPPFylO0xxOODkcIYQVJ4IXMvNd2cB//wzR4oLFnmRDCZUkCL0ROnoQKS8eT6F2MwEEvODocIYSVJIEXIjNGHqdTxnySu/U2NkkQQrg0SeCFxPnzUOTbiSgFQe+87OhwhBA2IAm8kPhq7CWeT53GpfAnoFy5vJ8ghHB6MopVCFy+DJcmfEMgF2HUYEeHI4SwEWmBFwJfT0mj99UJXKx1HzRs6OhwhBA2IgnczSUnw/6PllCBoxR7V1rfQrgTSeBubs4ceD7+M66WvhMef9zR4QghbEgSuBvTGlZ9EEETfsPvjYFgMuX9JCGEy5AE7sY2boRWRyaT6lMElbnXpRDCfUgCd2Nff3KeLvyE6vYsBAY6OhwhhI1JGaGbio6G0mu+w48kGNDP0eEIIexAWuBu6ovPM+jPFFIaNoVatRwdjhDCDiSBu6Hz5+H4jHVU4gjeA/s7OhwhhJ1YlcCVUoOUUn8qpfYppX5SSvnaKjBx6776CnqmTCateAh06uTocIQQdnLLCVwpFQq8AtTXWlcHTMDTtgpM3JqUFFj0+XEeYzmefXuBj4+jQxJC2Im1XSiegJ9SyhPwB05YH5Kwxvz50O7MNDyUhhdfdHQ4Qgg7uuUErrWOBcYCx4GTQILWek3u85RSfZRSEUqpiLi4uFuPVORJa5g4LoUXTdOh7aOy6qAQbs6aLpTiQDugAlAGKKKUejb3eVrraVrr+lrr+iEhIbceqcjT5s1QLnIxIemnUS/J4KUQ7s6aOvCHgH+01nEASqlFQBNgli0CEzdv3DgY5jmZjLJ34tG6NQBLImMZszqaE/GJlAnyY2h4VdrXCXVwpEIIW7AmgR8HGiml/IFEoCUQYZOoxE2LjoZ/lu+lCZuh36fg4cGSyFiGL9pLYmo6ALHxiQxftBdAkrgQbsCaPvDtwAJgF7A387Wm2SgucZM++wwGekxC+/pCz54AjFkdnZ28sySmpjNmdbQjQhRC2JhVU+m11u8C79ooFnGL4uJg2cwLfOExC/XMMxAcDMCJ+ESz51s6LoRwLTIT0w1Mngxdkr/FJ+0qDBiQfbxMkJ/Z8y0dF0K4FkngLi4xESZPymCo/5dw331Qu3b2Y0PDq+Lnde0a4H5eJoaGVy3gKIUQ9iCrEbq4WbOg/tlfKM3fMOCjax7LGqiUKhQh3JMkcBeWkWGUDn4TMAldtDTLKtzLp6M3XJesJWEL4Z4kgbuwlSshI/oQTVnFgWeGMGz5QSkZFKIQkT5wFzZuHLxedDLay4vXgxtJyaAQhYy0wF3Utm3wx6bLrPL9FtW5M/vSzVeWSMmgEO5LWuAu6qUhSfTw/hafpIv0KdaIIH8vs+dJyaAQ7kta4C5o9HdxRP1+G7P9J7C3eEXWFKuAV1IaXiZFarrOPk9KBoVwb9ICd0GffOTBo95LufvqX3xb/3FQitQMTRFvT0KD/FBAaJAfH3esIQOYQrgxaYG7mG3bIP5wMO8UG0WMbwjL7n4w+7GExFSi3m3twOiEEAVJWuAuZtQoaO67noYXo/jq3k6kmf77HSz93UIULtICdyHbtsGaNbC73GjOng1iXo1W2Y9Jf7cQhY+0wF3IqFHQImgnNY+t5/QLfbktJEj6u4UoxKQF7gKWRMbyzrST7FtTn59LvUVq0WLc8/4wtgYGOjo0IYQDSQvcyWXtqnNkVTnu9tnDw6fX8E2tR1jy92VHhyaEcDBJ4E5uzOpozh0oQdLREN4u/i7Jnt5Mq/uYTJEXQkgXirOLiUvm/LpGVAz6k85nlvNDnbac9w9EyRR5IQo9SeAOlJ8d4zMi7yI9wZ/hlUeiLyq+btgBkJJBIYQkcIfJz47x0dFw6n/lqV7ld575exlL7mnGyWIhUjIohACs6ANXSlVVSkXl+HNRKfWqDWNza3ntGK+1sb1lEX/FittGoZUHn9/3jJQMCiGy3XILXGsdDdQGUEqZgFhgsW3Ccn957Rg/bx6sWwfL+q+i3ORV8NFHbBveowAjFEI4O1tVobQE/tJaH7PR67m9G+0Yf/EiDBoEjeok8+jaV6ByZRg8uIAjFEI4O1sl8KeBn8w9oJTqo5SKUEpFxMXF2ehyru9GO8aPGAGnTsG8Jp+jDh+GL74AHx8HRSqEcFZKa533WTd6AaW8gRPAPVrr0zc6t379+joiIsKq67kTc1Uo2zf5MnpIMFWrbWfX4eYkNHmQ0htXOTpUIYQDKaV2aq3r5z5uiyqUNsCuvJK3q8lPiZ+1cu8Y//nc03zyRiDetyfwadobeGRk0O2epxkQGSuDlkKI69iiC6ULFrpPXFVWiV9sfCKa/0r8lkTG2u2aZ87AsL7F8PBJo13D6TwevZkp9z7BkaIhMutSCGGWVQlcKeUPtAIW2SYc55BXiZ+tJSdDx46Qctmbco9t5sP/fcHxwFJMvbcTIBsTCyHMs6oLRWt9FShho1icRl4lfrakNfTrB1u3wl1P7WPijrepcD6W7k++R7KXMXApsy6FEObIYlZm3KjEz9bGjoVvv4URI2C5aRIP/fUHIx96ka3lawOyUYMQwjJJ4GbcqMTPVpKTjZmWr78OTzwB75acTKXZ0/mrS082tnhCNmoQQuRJ1kLBfMXJxx1r2K0K5d9/oXNn2L4dXnsNPm62Go92r0DbtlT84Su2mkx5v4gQotCzug78ZjhjHXjuRaXAaG3bq+W7di106QIpKUbXSaeq+6BpUyhfHv73PwgIsPk1hRCuzVIdeKHvQimoipO4OHjzTQgPh9tvhz/+gE5FVkGrVuDvDytWSPIWQtyUQt+FcjMVJ+a6WoAbdrVERMDEiTBnjtHq7tYNpnxykSIjhsD06VCtmvFg2bL2uUEhhNsq9Am8TJAfsWaSde6KE3Prdw+dvxsUpKbr7GPDFu7l+F+e+J4vxYwZRj93kSLQqxe89BJUO7UBGj0PMTHwxhswciT4+tr9PoUQ7qfQJ/Ch4VXN9oHnrjgx19WSkq5Jv+xDyulAUk4EkXwiiOSTQQxM8QKMRQQnTIAej52j2M6N8PFSmDULqlQx+rsbN7b/DQoh3FahT+BZ3R15VZyciE/EJzmJ8vsS8Iv1ICW+KMkXAkhN8iaDZDyIpUjQfvzCLlCk+AU+ebYMpU9Fob5fC6/uMmbsBAQYy8K+/77R7y2EEFYo9FUoNxQRAatWcXHrHs78upMKiUcxkZH/53t5Ga3shx4y/jRoAJ6F/nemEOIm2XM1Qvdz/Dh62DDUT8YaXWeoyF5qMa94Rw5VKsnFOxUmUzreaDzQ6LR00j08SDV5YvL2pm+ru2leI9QoDSxa1LH3IoRwW5LAc7p8GT79FD1mDKmp8Alv83PVIXTuHcQzz4A6GcvP+ahCaS4zJ4UQBUASeJalS6F/fzhxguVFujAwZTQ937+DrcMha2Jk+9tDzU7ukanuQghHkAQOsGIFulMnTpWqRWePBcSGNGb2bNsWiRTEBhFCiMKl0CXw3In041IXuf+lrhwJqEPdExt4vGsAP0+GwEDbXaP5XSEs3Bl7TQ358EV7AWm9CyFuXaGaSp97p50ihw9Qq383Yj1DaRK/kg8nBPDjj9Yn79y7+fz4+/EC3SBCCFE4FKoEnnMyTljCab6fN4Kr2p/7Lq2l26AQXnnFttfIYqlQU3baEUJYo1B1oWQlzBJX4vl+7jv4JqdyX+pGzlT0Y8wY214jP2SnHSGENQpVC7xMkB9ozRfLP6X0xbM8ynIOhZSnZreD2GoJbktJWeX6XnbaEUJYq1Al8KHhVel4eCtNj+1hqPcn/ObViDue2sWwxyvZ9BrmdvN5ptEdhAb5yU47QgibsaoLRSkVBEwHqmN09b6gtf7NBnHZRfvKgYRv/Y7dPjWYktyfe/rs4b0XKts0keZ3bRUhhLCWtX3gE4BVWusnlFLegHOv0PTBB/idOUVfFjFmnDeDB1+3tIBNtK9jfsKPEELY0i0ncKVUMeABoAeA1joFSLFNWHYQHY3+7DMWBvQgIawxL7/s6ICEEMI61vSB3wnEAd8qpSKVUtOVUkVyn6SU6qOUilBKRcTFxVlxOStoDa+8QoqHHy9dGs2ECcZCgUII4cqsSeCeQF1gita6DnAFGJb7JK31NK11fa11/ZCQECsuZ4UlS2DNGt7R79G4XSlatXJMGEIIYUvWJPAYIEZrvT3z+wUYCd1pLImMpcV7K4nt8SL7fe9iou7PZ585OiohhLCNW07gWutTwL9Kqaxi5pbAfptEZQNZU9ofXzOL0Itx9E36iiINj7MnIdbRoQkhhE1YWwf+MvCjUmoPUBv4yOqIbGTM6mjUlcs8H7GMJT6Psq1oA/wbHpb1R4QQbsOqMkKtdRRgn1o8K52IT6Tb3nUEJl/hE96keOuDeHiny/ojQgi34bYzMcOKedPzj6X8ZrqXyDJ34X/3CUDWHxFCuA+3XcxqrM8xyiWcYggTKdb4CErJ+iNCCPfiti3whou/47hXeVbd1hr/imdk/REhhNtxzxb49u2orVsZx+eM/6gYvXu3dXREQghhc26TwHNuYzZ95VgamgJZXvwFPunm6MiEEMI+3CKBZ9V8J6amE5Zwmmb7NjNOD6HxkxpfX0dHJ4QQ9uEWfeA5tzHrEbGMDO3BRM/+HA6OcmxgQghhR27RAs+q7Q5IvsJTu9cylydJqAmeqZccHJkQQtiPW7TAs2q7n9y9hoDUq3zGYALq/yM130IIt+YWCXxoeFX8PD3oErWaLaop0VXLUKxkstR8CyHcmlsk8PZ1QplaOZVKF2KYoXtSsUWs1HwLIdyeW/SBAzywbSVXlT//NniCPZMDHB2OEELYnVu0wElMJH32HBboTjzTV5K3EKJwcI8EvnQpnpcTmOvbg86dHR2MEEIUDLfoQkmb/h0n1B2U6dqMokUdHY0QQhQM12+Bx8bisWEt3+nuvNDL9W9HCCHyy/Uz3g8/4KEz2FqxO40aOToYIYQoOK7dhaI1ydO+Yzv307pfRZRydEBCCFFwXLsFvn07Pv9EM8ujO91k1UEhRCHj0i3w9BnfkYwfVx7pTMmSjo5GCCEKllUJXCl1FLgEpANpWuuC2+A4KYn0H+ewkE507VuswC4rhBDOwhYt8OZa67M2eJ2bs3Qp3lcTWFGiBz+GF/jVhRDC4Vy2CyXp6x84Q1kq9W6Op8vehRBC3DprBzE1sEYptVMp1cfcCUqpPkqpCKVURFxcnJWXy3T+PF4bVzOXp3i+p2uPwwohxK2yNvs11VrXBdoALymlHsh9gtZ6mta6vta6fkhIiJWXy3zNRYsxZaRxsPbTVKpkk5cUQgiXY1UC11qfyPz7DLAYaGiLoPJyafocjlCRhi/WLYjLCSGEU7rlBK6UKqKUCsj6GmgN7LNVYBadOUPRHRuYp57mic4yc0cIUXhZM/xXClisjOmPnsBsrfUqm0R1AxnzFuChMzj54FOUKGHvqwkhhPO65QSutf4bqGXDWPLl4tdziaUa9/WtXtCXFkIIp+JaJRyxsRTbs4XF3k/x2OPSfSKEKNxcqoI67af5eKI51PhxWn2xgRPxiZQJ8mNoeFXZ/1IIUei4VAK/OG0Ox6jNxjJpmOITAYiNT2T4or0AksSFEIWK63Sh/PMPwYe3s8D7CTzCzlzzUGJqOmNWRzsoMCGEcAyXSeCJ388DYFmV+1Ee+rrHT2S2yIUQorBwmQR+dcZcttMQdb+/2cfLBPkVcERCCOFYrpHADx2ixPFI1pV4mlG9SuPnZbrmYT8vE0PDqzooOCGEcAyXGMS8NH0uAYDPs53pUDcUpWDM6mipQhFCFGpOn8CXRMaydbUXVejFbFMMlSIV7euESsIWQhR6Tp3Al0TGMnzRXk4FtyOp/G2U8trB8EXxgJQMCiGEU/eBj1kdTWJqOoH3/k2pp3YAUjIohBBZnDqBWyoNlJJBIYRw8gRuqTRQSgaFEMLJE/jQ8KpSMiiEEBY49SBm1kCllAwKIcT1nDqBA1IyKIQQFjh1F4oQQgjLJIELIYSLkgQuhBAuShK4EEK4KEngQgjhopTW12+OYLeLKRUHHLvFp98GnLVhOK5A7rlwkHsuHKy553Ja65DcBws0gVtDKRWhta7v6DgKktxz4SD3XDjY456lC0UIIVyUJHAhhHBRrpTApzk6AAeQey4c5J4LB5vfs8v0gQshhLiWK7XAhRBC5CAJXAghXJRLJHCl1MNKqWil1BGl1DBHx1MQlFJHlVJ7lVJRSqkIR8djD0qpGUqpM0qpfTmOBSul1iqlDmf+XdyRMdqahXseqZSKzXyvo5RSjzgyRltSSpVVSm1USh1QSv2plBqYedxt3+cb3LPN32en7wNXSpmAQ0ArIAb4A+iitd7v0MDsTCl1FKivtXbbyQ5KqQeAy8D3Wuvqmcc+Bc5rrUdn/rIurrV+w5Fx2pKFex4JXNZaj3VkbPaglCoNlNZa71JKBQA7gfZAD9z0fb7BPT+Jjd9nV2iBNwSOaK3/1lqnAHOAdg6OSdiA1nozcD7X4XbAzMyvZ2L84LsNC/fstrTWJ7XWuzK/vgQcAEJx4/f5Bvdsc66QwEOBf3N8H4Od/jGcjAbWKKV2KqX6ODqYAlRKa30SjP8IQEkHx1NQBiil9mR2sbhNd0JOSqnyQB1gO4Xkfc51z2Dj99kVErgyc8y5+31so6nWui7QBngp86O3cE9TgIpAbeAkMM6h0diBUqoosBB4VWt90dHxFAQz92zz99kVEngMUDbH92HACQfFUmC01icy/z4DLMboSioMTmf2IWb1JZ5xcDx2p7U+rbVO11pnAF/jZu+1UsoLI5H9qLVelHnYrd9nc/dsj/fZFRL4H0BlpVQFpZQ38DSwzMEx2ZVSqkjm4AdKqSJAa2DfjZ/lNpYB3TO/7g4sdWAsBSIrkWXqgBu910opBXwDHNBaf5bjIbd9ny3dsz3eZ6evQgHILLf5HDABM7TWHzo2IvtSSt2J0eoGY+Pp2e54z0qpn4BmGMtsngbeBZYA84A7gONAZ6212wz6WbjnZhgfqzVwFHgxq3/Y1Sml7gO2AHuBjMzDb2L0Cbvl+3yDe+6Cjd9nl0jgQgghrucKXShCCCHMkAQuhBAuShK4EEK4KEngQgjhoiSBCyGEi5IELoQQLkoSuBBCuKj/A7Z+nr4pZhBCAAAAAElFTkSuQmCC\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.996705\n",
       "x1                  0.490160\n",
       "np.sin(x1)          0.570882\n",
       "I((x1 - 5) ** 2)   -0.018322\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.212187\n",
       "1    11.077255\n",
       "2    10.816221\n",
       "3    10.480103\n",
       "4    10.136062\n",
       "5     9.850954\n",
       "6     9.674961\n",
       "7     9.629312\n",
       "8     9.701076\n",
       "9     9.846333\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
}
