{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prediction (out of sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:33.146379Z",
     "iopub.status.busy": "2021-02-02T06:51:33.145621Z",
     "iopub.status.idle": "2021-02-02T06:51:33.379390Z",
     "shell.execute_reply": "2021-02-02T06:51:33.378912Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:33.386356Z",
     "iopub.status.busy": "2021-02-02T06:51:33.383012Z",
     "iopub.status.idle": "2021-02-02T06:51:34.303579Z",
     "shell.execute_reply": "2021-02-02T06:51:34.304838Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "plt.rc(\"figure\", figsize=(16,8))\n",
    "plt.rc(\"font\", size=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Artificial data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:34.309751Z",
     "iopub.status.busy": "2021-02-02T06:51:34.308995Z",
     "iopub.status.idle": "2021-02-02T06:51:34.317285Z",
     "shell.execute_reply": "2021-02-02T06:51:34.318250Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:34.322852Z",
     "iopub.status.busy": "2021-02-02T06:51:34.321478Z",
     "iopub.status.idle": "2021-02-02T06:51:34.353434Z",
     "shell.execute_reply": "2021-02-02T06:51:34.354390Z"
    }
   },
   "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:                     1353.\n",
      "Date:                Tue, 02 Feb 2021   Prob (F-statistic):           7.53e-45\n",
      "Time:                        06:51:34   Log-Likelihood:                 9.2333\n",
      "No. Observations:                  50   AIC:                            -10.47\n",
      "Df Residuals:                      46   BIC:                            -2.818\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.9077      0.071     68.652      0.000       4.764       5.052\n",
      "x1             0.5179      0.011     46.976      0.000       0.496       0.540\n",
      "x2             0.4952      0.043     11.426      0.000       0.408       0.582\n",
      "x3            -0.0212      0.001    -21.864      0.000      -0.023      -0.019\n",
      "==============================================================================\n",
      "Omnibus:                        0.349   Durbin-Watson:                   2.124\n",
      "Prob(Omnibus):                  0.840   Jarque-Bera (JB):                0.521\n",
      "Skew:                           0.059   Prob(JB):                        0.771\n",
      "Kurtosis:                       2.514   Cond. No.                         221.\n",
      "==============================================================================\n",
      "\n",
      "Notes:\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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:34.358836Z",
     "iopub.status.busy": "2021-02-02T06:51:34.357459Z",
     "iopub.status.idle": "2021-02-02T06:51:34.368689Z",
     "shell.execute_reply": "2021-02-02T06:51:34.369677Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.37860807  4.86941945  5.3208846   5.70601499  6.00756216  6.22085149\n",
      "  6.35455029  6.42924376  6.47405292  6.52185006  6.60385799  6.74452017\n",
      "  6.95748471  7.24336207  7.58962472  7.97266541  8.36167575  8.72370795\n",
      "  9.02908805  9.25629127  9.39547846  9.45011268  9.43639015  9.38057936\n",
      "  9.31470514  9.27128714  9.2779969   9.3531115   9.50251071  9.71871094\n",
      "  9.98209478 10.26413409 10.53207707 10.75432932 10.90564469 10.97127316\n",
      " 10.94938421 10.85136759 10.69996201 10.52551924 10.36101737 10.23664264\n",
      " 10.17482992 10.1865767  10.26963625 10.40888677 10.57881582 10.74771197\n",
      " 10.88287489 10.95598792]\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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:34.377097Z",
     "iopub.status.busy": "2021-02-02T06:51:34.372678Z",
     "iopub.status.idle": "2021-02-02T06:51:34.382005Z",
     "shell.execute_reply": "2021-02-02T06:51:34.382955Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10.93372496 10.78001952 10.51429186 10.18079828  9.83779567  9.54327813\n",
      "  9.34077806  9.24870787  9.25585197  9.32411278]\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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:34.387208Z",
     "iopub.status.busy": "2021-02-02T06:51:34.385835Z",
     "iopub.status.idle": "2021-02-02T06:51:34.757276Z",
     "shell.execute_reply": "2021-02-02T06:51:34.758305Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd9b3540150>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6MAAAHWCAYAAACCFl9HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAB0TElEQVR4nO3dd3RU1RrG4d9OCBBq6CVUQQMoCBKwIE3AqFRRBLGAClhARYoSsICoVLkgRQUpKoiFEhUQBJFqwUBQekckdDDUACnn/rFpwYQ6mTNJ3metWSFnzpx84507yTt7728bx3EQERERERER8SY/twsQERERERGRjEdhVERERERERLxOYVRERERERES8TmFUREREREREvE5hVERERERERLxOYVRERERERES8LpPbBeTPn98pVaqU22WIiIiIiIhIKli+fPkBx3EKXHzc9TBaqlQpIiMj3S5DREREREREUoEx5u/kjmuaroiIiIiIiHidwqiIiIiIiIh4ncKoiIiIiIiIeJ3CqIiIiIiIiHidwqiIiIiIiIh4nevddC/nyJEj7Nu3j7i4OLdLkXQkICCAggULkitXLrdLERERERHJkHw6jB45coS9e/cSHBxMYGAgxhi3S5J0wHEcYmNjiY6OBlAgFRERERFxgU9P0923bx/BwcFky5ZNQVQ8xhhDtmzZCA4OZt++fW6XIyIiIiKSIfl0GI2LiyMwMNDtMiSdCgwM1PRvERERERGX+HQYBTQiKqlGry0REREREff4fBgVERERERGR9EdhVERERERERLxOYTQVtG3bFmMMxphzW4jUrVuXkSNHXtUaxQULFmCM4cCBA6lYrYiIiIiIiPcpjKaS+vXrs3v3brZv386PP/5I48aNeeutt6hZsybHjx93uzwRERERERFXZYgwGhEVTY3+8yndYyY1+s8nIio61X9mlixZKFy4MMHBwVSuXJkuXbqwYMECVqxYwcCBAwGYOHEi1apVI2fOnBQsWJAWLVqc2/ty+/bt1K1bF4ACBQpgjKFt27YAzJ49m5o1a5InTx7y5s1LWFgY69atS/XnJCIiIiIi4inpPoxGREUTPm0V0TGxOEB0TCzh01Z5JZBe7JZbbuG+++5j6tSpAJw+fZo+ffrw559/MmPGDA4cOMCjjz4KQPHixc+dt2bNGnbv3s2wYcMAOH78OJ07d2bZsmUsWLCA3Llz07hxY06fPu315yQiIuIJbnxwLCIi7srkdgGpbdCcDcTGJSQ5FhuXwKA5G2hWJdjr9VSoUIF58+YB8PTTT587fsMNN/Dhhx9Svnx5du7cSbFixcibNy8ABQsWJH/+/OfOfeihh5Jcc/z48eTKlYtly5Zx9913e+FZiIiIeM7ZD47P/r4++8Ex4MrvahER8Y50PzK6Kyb2qo6nNsdxzu1vuWLFCpo2bUrJkiXJmTMnoaGhAOzYseOS19iyZQutW7emTJky5MqVi0KFCpGYmHjZx4mIiPiiS31wLHKlNLoukvak+zBaNCjwqo6ntrVr13LDDTdw/PhxwsLCyJYtG59//jl//PEHs2fPBrjsdNvGjRuzf/9+Pv74Y37//XeioqLIlCmTpumKiEia5GsfHEva40vLskTkyqX7MNo9LITAAP8kxwID/OkeFuL1WlavXs3s2bN5+OGHWb9+PQcOHOC9996jVq1alCtXjn379iU5P3PmzAAkJJz/tPjgwYOsW7eOnj17Ur9+fcqXL8/Ro0eJj4/36nMRERHxFF/74FjSHo2ui6RN6T6MNqsSTL/mFQkOCsQAwUGB9GteMdXXoJw6dYo9e/awa9cu/vzzT4YMGUKdOnWoWrUq3bp1o0SJEmTJkoURI0awdetWZs6cyRtvvJHkGiVLlsQYw8yZM9m/fz/Hjh0jT5485M+fnzFjxrB582YWLlzIc889R6ZM6X75r4iIpFO+9MGxpE0aXRdJmzJEgmlWJdjrDRDmzZtHkSJF8Pf3JygoiFtuuYW33nqLZ599lsyZM5M9e3Y+/fRTevbsyciRI6lUqRJDhgzhvvvuO3eN4OBg+vTpQ69evWjXrh1PPvkkEyZM4KuvvuKll17illtuoWzZsrz//vv/aWokIiJyoYioaAbN2cCumFiKBgXSPSzEZ5oDna3DV+sT31c0KJDoZIKnRtdFfJtxHMfVAkJDQ53IyMhk71u3bh3ly5f3ckWSkeg1JiIZwcXdasGOPHpjppCIN+g1LuLbjDHLHccJvfh4up+mKyIiktFpPZ2kd24tyxKR65MhpumKiIhkZFpPJxmBG8uyROT6aGRUREQknVO3WhER8UUKoyIiIumcutWKiIgvuqIwaoypZYz5zhgTbYxxjDFtL7q/uTFmjjFm/5n766RCrSIiInINtJ5ORER80ZWuGc0BrAY+O3O7WHbgF2BiCveLiIiIi7SeTkREfM0VhVHHcWYBswCMMROSuf/zM/fl92RxIiIiIiIikj5pzaiIiIiIiIh4nSth1BjTwRgTaYyJ3L9/vxsliIiIiIiIiItc2WfUcZzRwGiA0NBQx40aRERExMO2bIEffrC3Vasgf34oUgQKF7a3s/++8Fj27G5XLXJNYo8ncvyEIX8B43YpImmWK2E0PTPm0m9Ibdq0YcKECd4pRkREJDWdPAkLF54PoBs32uNly0KtWvDvv7BnD6xcCXv3QkLCf6+RM+f5YFqsGDz8MDRpApn0J4r4nvjDx1k9eDYnJk2nwrYZBPmdgEL57QcvZ28FCiT9/uJjWbO6/TREfIbe6T1s9+7d5/49Y8YM2rdvn+RYYGDSDcbj4uIICAjwWn0iIiLX5cLRz59/hthY+8d1nTrQqRPcf78NoxdLSICDB2H3bhtQ9+w5/++zX3/+GSZPhhIl4IUXoF07yJfP609R5ELOwUNsGfo9JyZN58Ztc6jMSQ6afKy+sTkl7yhM8Sz74cABe1u1Cvbvh0OHwElh8l/27OcDatWq8Prr9oMYkQzoisKoMSYHcPY3ix9QwhhTGTjkOM4OY0xeoAQQdOacssaYGGCP4zh7PFqxjytcuPC5fwcFBSU5tn37dooUKcIXX3zBmDFj+PXXXxk0aBA5cuSgU6dOHDt27NxjFyxYQN26ddm/fz/589smxb/88gvh4eH88ccf5MmThyZNmjBgwABy5crlvScoIiIZy6lTsGBB8qOf7drZ8FmnDpz5sDUmBn78Gv76y+bIQoXsrXBhfwoVKkjeigXxu/XW5H9WfDx8/z0MHw49ekDv3tC6Nbz4IlSunPrPVeSsXbvY/WEEJyZOo+T2BZQlgZ0UY0HZ9uR88kGqdanJ3dkv8Wd0QoKdGbD/gqB64EDS7/ftg/Hj4dNP4ZVX4LXXIHdu7z1HER9wpSOjocDPF3zf58ztU6At0AQYf8H9Yy44r/d1VZgOhYeHM3jwYMaOHUtAQADz5s277GNWrVrFvffeS58+ffjkk084dOgQnTt35umnn2bKlCleqFpERDKaxZ9MofRrL1Ps0C5OZcpMTPW7KHTR6KfjwNq1MHOmvS1dmvxs3LMyZbIDQoULnw+qNqxC0aKZuO++B8n94IOwejWMGAGffw7jxsHdd8NLL0GzZqAZRZIaNm0iZvx0TkyaTtEdv1EEWE8IU0p3J8cTzanVJZT7c1/h+lB///PTci9l+3Y7MtqvH4weDW+8Ac8/D5kzX++zEUkTrnSf0QVAiv/vcxxnAjDBIxVdRufOdumJN1WuDEOHeu56L774Ig8//PBVPWbQoEG0bNmSrl27njv24YcfUqVKFfbt20fBggU9V6CIiGRsR4+y9ZkXqfnNp/wdVJgOD/ZiYenb8MuWjX53V+TeosH8fCZ8zpoFf/9tH3brrfDqq9CwIdx+Oxw9amff7t17/nbx96tX269xcfYa2bPDk09Cx463cPNHH9k/0seNg5Ej4ZFHIDjY/rHevj3od594QPy3MzncqRf5dv5JELCZqnxX7B2yP9Gce18uT6tCqfjDS5WCiROhSxf7f57OneGDD+C996BFC/DTLoySvmnNqAtCQ0Ov+jHLly9n8+bNfPXVV+eOOWfWImzZskVhVEQkA4mIimbQnA3siomlaFAg3cNCaFYl2DMXnzsX2ren1I4djA1tyuCaTxCbOSvxhwOJXVWQp77OwskdtndRtmxQvz707AkPPPDfZW958thb+fKX/pGOY6f3rl8PY8bY7Pnhh1C3LnTqlIcmL3clU+fONvkOH25Hkt5+G1q1slN4r+H3qgibNnGo7Svk/WUm+wnh4wL/I8cTD9LwhZI8V8bLtdx2m/3/3o8/2lDaqhW8/z4MHGinwYukU2kujHpyhNIt2S9qY+/n53cuWJ4Vd/Yj4jMSExNp164dr7zyyn+uFxzsoT9ARETE50VERRM+bRWxcXYubHRMLOHTVgFcXyA9fBi6dYNPPoGbbuKR1gP4o2gFjq4oxbGVJYg7mBOATEHHeaGDHf2sVcszjUGNsaH1zjvtbeBAGDsWRo2Chx6C4sXh+ef9adeuMQUaN7apdcQIu9bus8+gcWObXvX7UK7EsWOcfP0dMg0fQqbErPTNPZhbP3mR8Icyc5lNEVKXMRAWZj/hmTjRTtmtW9f+n23AALj5ZheLE0kdGvv3AQUKFODEiRMcOXLk3LGVF81Fvu2221izZg1ly5b9z+3iDr0iIpJ+DZqz4VwQPSs2LoFBczZc+0VnzbJ/6I4bZ0dlVq5ka6Ga7P3iTv796Wb8ssaR5561FG2/gGqv/s6wYXDvvam3Q0X+/LaXy9atEBEBN91kR1+LFYM2bSDyWDkbRnfuhP79Yd688/Wn1MFUxHFwJn3BiRIhZB02gEmJrRny7EZe2dmVJg+7HEQv5O9vX+gbNtgQumQJVKpkG4ZFR7tdnYhHKYz6gNtvv53s2bMTHh7O5s2bmTp1KqNGjUpyzmuvvcayZct47rnniIqKYvPmzcyYMYNnn33WpapFRMQNu2Jir+r4JR06ZP/obdgQgoLgt99IeG8A748KZM3I24k7mIN8jVZS6LFfyVVtG7kKneTV+0Ku7wlcBX9/aNrUZs21a+0y0WnToFo1uOMOmPh9bk51fs227r31VnjmGbjvPtixw2s1ShoRFUVstZqYxx9j7b9FeKb8L1RaMYHeHxUmRw63i0tBYKD9cGjLFnj5ZdvM68YboVcvO5NBJB1QGPUBefPmZdKkScydO5eKFSsyevRo+vbtm+ScSpUqsWjRIrZv307t2rW59dZbCQ8Pp1Ch1FxVLyIivqZoUPKzYVI6nqKICDuaOGmSXYO5fDkbclWjVi07WzfsXsOYbw8SUuMQfgaCgwLp17yi59amXqXy5e1gaHS07e/y77/wxBN2S9JhM8uSMO9n2+Ro6VL7vD76CBITXalVfMjBgyR0eJ7EqqEcXb6Rjlk+4bdhyxi96k6qVLn2y0ZERVOj/3xK95hJjf7ziYhKxRHLfPlgyBA7Pf3BB21zo7Jl7frp06dT7+eKeIG5eK2it4WGhjqRkZHJ3rdu3TrKX67rgch10GtMRNKai9eMAgQG+F95UNy/326R8uWXtl38uHEkVKrC0KE2kwYG2rD32GP4zrTFZCQmwk8/waBBtu/LXXfZdablsm63Q6jz5tn1dp98Ajfc4Ha54m0JCfDxx8T1eB1z9Agj6ciKJn3o92EQRYte36Wv+/+D12v5cjtiOn++ncM+axaU8XbHJZGrY4xZ7jjOf7rNaWRUREQkDWlWJZh+zSsSHBSI4SpHLKdNs6OGU6dC376wbBkbslWhZk07GtqgAaxZA48/fu1B1FsjRn5+tt45c2yvl/XrbbYe8FUp4mf9aNvyLl8OFSvadK1R0oxj0SLiK1eFjh1ZfLQyDxRZSdkZw/j02+SD6NW+ZlNl3fbVqFrVftgycyYcPAg1a9r/44qkQWmum66IiEhG16xK8NWPwIwcCZ062W1Q5s8nofwtSUZDJ06E1q2vbzQ01Tr9XuLnnd3i5obnchHweyg9egTyzTeGcePaUWl1GDz7rF1v9/XXtsHRTTd5vA7xEbGxOF26Yj76kN1+Jejq9w2luz7E9LcMF21kcM61vGY9um77Whlj91NatMh2361Vy34yo22OJI3RyKiIiEh69/77Nog2bQpLlrAh4JZzo6H33msHVTwxLdebI0ZnQ0R0TCwOsD/hCP/euZDuAw/yzz928Kj32OKcnj7TbgGzZo1tcjR4sJ3CKenLunUkVLsD89GHDKYrT1Rdx+tRDzNgYMpBFK7tNeuxddueUKGC7babOzfcc48NpyJpiMKoiIhIevbOOzZ1tmhBwpff8P6ILFSubKe1Tpxo+xgVKeKZH+XNEaOUQsTShD9ZuxZatYI+fSC0miGywpO2HW9YGHTvbheYrl3r8ZrEBY4D48eTeFsoMet309DvB8zgwcz/LRuVKl3+4dfymu0eFkJggH+SY4EB/nQP816n6SRuuAEWL7Z7H4WFwQ8/uFOHyDVQGBUREUmPHMfOwX3jDXjiCWJGfcE9YQEeHw29kDdHjC4VIvLls7tgfP+9XVJ3++3QY1gRYr+YDpMn260yQkNhwgSP1yVedOSIXeD89NMsPH0H9fOtJHzhfXTtatcUX4lrec1e17rt1BIcDAsX2rbTTZvCN9+4V4vIVVAYFRERSW8cx44AvvsutGvHoSETqH9fJn791eYvT46GXsibI0ZXEiIaNbKh++mnYcAAqHKbYWnxVrB6Ndx5Jzz1lL3zxAmP1yfX57JNhSIjSaxyG4mTv6QX7/BurR+Z/VdR7r776n7Otb5mm1UJZmmPe9jWvyFLe9zjbhA9q0AB+PlnqF7dTg0YP97tikQuS2FUREQkPUlMhBdfPLdO9MC7H3NPfT9WrYLp06FNm9TbssWbI0ZXGiKCgmxj3R9/hJMnbePRzv0Lc3z6j3bUeMIEuOMO2OClTqhyWRevBz7bVCgiKtp+0DJ0KM5dd7Hvn1PUchbihPdi9lx/rmXrdZ8c5bweuXPbRkb169sPWj74wO2KRC5J+4xKhqbXmIikKwkJtnvs2LHQrRt7uw6kfgPD5s12NDQszO0CPevCbrpFgwLpHhZyyRBx7BiEh8OIEVC2rJ3JWHnvHDvV8+RJuydpy5ZefAaSnBr95xOdzDTsCgGnmbVyHMyYwcxMTemUbRwfTMxL48YuFOnrTp2CRx+1n0C98w707OnbGwdLupfSPqPa2kVERCQ9iI+3004nToQ33mD3c324p65hxw67HeE997hdoOdd7RY3OXLA8OHw8MN2vewdd8CwYWF0WBGFadXSTm1ctAiGDIEsWVKxcrmU5NYD375jFcO+H0x87FE6M5xfK3bkpymGG25wocC0IEsWu53R00/bteOHD9u56gqk4mM0TVeuW6dOnahTp86579u2bUujRo2u65q9e/fmlltuuc7KREQyiLg4u0noxInw7rvs7PA2tesYdu6E2bPTZxC9HrVrQ1QU1KkDzz0Hj/coxrEZC6BrVxg1CmrUgG3b3C4zw7pw3a9fYgKdl0ziiy97cfRUHkITfud0+04s/UVB9LIyZbLT0F94AQYNguef17ZG4nMURlNJdHQ0HTp0oFixYmTOnJng4GDat2/Pzp07k5x3ueD2559/0rRpUwoXLkzWrFkpUaIEDz30EH///XdqP4VrNmzYMCZOnHhF527fvh1jDBdP1e7WrRsLFy5MjfJERNKXU6fsUN8338D777O9dU9q1YK9e+3SsZo13S7QNxUoALNm2RmMX34JoXcGsKrNYDufefNmuO02+O47t8vMkM6uBy509ACTv+xF56WT+cL/UW43K3hlQmVGj4asWd2uMo3w87Pz0nv0gI8/hieftB9eifgIhdFUsG3bNkJDQ1m9ejWffvopmzdvZuLEiaxZs4Zq1aqxffv2K7rO/v37qVevHjly5GDmzJmsX7+ezz//nDJlynDkyBGP1nz69GmPXSt37twEBQVd1zVy5MhBvnz5PFOQiEh6FRtrt3H47jsYOZKtzbpQuzb8+y/MnWu305SU+flBr17w0092FuPtt8P4Q01hxQooU8b+t+3eXX+8e1mzKsGMK7CPORNe4uZdW3jSjOe1IhNY8EdO2rRxu7o0yBjo1w/eew+++MJ+eHXypNtViQAKo6miY8eO+Pn5MW/ePOrVq0eJEiWoW7cu8+bNw8/Pj44dO17RdZYuXcq///7L+PHjqVq1KqVKlaJ27doMHDiQihUrpvi4s6Ot77zzDoUKFSJHjhw89dRTxMaeX4NRp04dnn/+ebp160aBAgWoUaMGAGvXrqVhw4bkzJmTggUL8uijj7Jnz55zj0tISKBbt27kyZOHPHny0LlzZxIumvJx8Wiv4zi8//773HjjjWTJkoVixYoRHh4OQOnSpQGoVq0axphz030vnqabmJhI3759KV68OFmyZKFixYp8++235+4/O8I6depUGjRoQLZs2ahQoQJz5869ov/WIiJpzrFj0LChbRM7diwb679ArVr28E8/2d0d5MrUqQMrV9rdXp5+Gtq+fQPH5yyx0xsHD7YnXDSzSVLJ6dPQtSt3vtyGfwNKUzUhiuMPtmXdqkxUquR2cWnc2e5d331n9z06dsztikQURj3t0KFDzJ49m44dO5ItW7Yk92XLlo0XXniBH374gX///fey1ypcuDCJiYlMmTKFq+16vHDhQv78809++uknpk6dyo8//shrr72W5JyJEyfiOA6LFy/ms88+Y/fu3dSqVYtbbrmFZcuWMW/ePI4dO0aTJk1ITEwE4P3332fMmDF8/PHH/PrrryQkJDBp0qRL1tKzZ0/69u1LeHg4a9as4ZtvvqF48eIALFu2DIDZs2eze/dupk2bluw1hg0bxqBBgxgwYACrVq3iwQcfpHnz5qxcuTLJeb169eKll17izz//pFq1arRq1YpjerMVkfTm2DG47z7bbOfzz1l359PUrm3/jv/5ZzvDVK5OoUI217/1Fnz2GVSvlZV1nUbakaQ//4QqVey8Z7lql90z9KytW+Huu2HIECbm7sgtR3+jw+AQpkyxO5aIB3TsCJ9+at8o7r3XTqMQcVHa66bbubP9+NKbKleGoUOv6NRNmzbhOE6K24VUqFABx3HYtGkT1S/zsfUdd9xBz549adOmDR07dqRatWrUqVOHxx57jJIlS17ysf7+/owfP54cOXJwyy23MGDAAJ555hn69etH9uzZATsq+f777597zJtvvsmtt97KgAEDzh377LPPyJs3L5GRkVSvXp2hQ4fy6quv8sgjjwA2JM65xC/nY8eO8b///Y+hQ4fy9NNPA1C2bFnuvPNOAAoUKABAvnz5KFy4cIrXGTx4MN26daN169YAvP322yxatIjBgwcnWZ/6yiuv0PhMj/f33nuPzz77jJUrV3L31e6CLSLiqxISbLOiX3+FL79kVbkW1Kttp5wuWAAVKrhdYNrl7w+9e9v+RY89BqGh8PHHj/J4ZBVo0QLuvx+6dbMLTTNndrvcNOHsnqGxcXYW1dk9Q4GknZC//hrat+d0vKFN5qksCGzO7O+15jlVPPmkbS3dqhXUrQvz5kH+/G5XJRmURkZTiUmhdfbZEc6U7r/Yu+++y549exg9ejQVK1Zk7NixVKhQgZ9++umSj6tUqRI5cuQ49/2dd97J6dOn2bJly7ljVatWTfKY5cuXs2jRInLkyHHudnYEc8uWLRw+fJjdu3efC5IAfn5+3H777SnWsXbtWk6dOkW9evWu6Pkm58iRI+zatevcVOKz7r77btauXZvkWKUL5vAULVoUgH379l3zzxYR8Tldu8L338Pw4ay8sQV160JAACxcqCDqKQ0a2M+9Q0PhiSegQc/i1G42iEm3hsGgQcTcGgobNrhdZpowaM6Gc0H0rNi4BAbNOfPfLzbWtjRu2ZJtWctz04mV7LqjOVFRCqKpqnlz+z6yfr3996lTblckGVTaGxm9whFKt9x4440YY1izZg3NmjX7z/3r1q3DGEOZMmWu+Jr58uWjRYsWtGjRgn79+lGlShX69u17XQEPODdCelZiYiINGzZk8ODB/zm3UKFC56bqXo2rnV58KckF+IuPBQQE/Oe+a6lbRMQnjRgBw4ZB585EVn+Be++BnDlh/nzbb0c8p2hRu/a2VYejTB2fk4Ald7G5WVYWlq7KwNkfEF+5CpmGfwDPPKO9Gy8huT1Dzx1fvx4eeQRWreLTwq/Sbs87dHk1gHfftbuSSCoLC4Px4+1Mi2eftf/Wa1m8TCOjHpY3b17CwsIYNWoUJ06cSHLfiRMnGDlyJPfffz958+a9putnzpyZMmXKXHYd5KpVqzh+/Pi573/77bdzj03Jbbfdxpo1ayhZsiRly5ZNcsuZMye5c+emSJEi/Pbbb+ce4zjOuXWfyalQoQJZsmRJcSQ385lpThc3QbpQrly5KFq0KEuWLElyfMmSJVTQMICIZBQzZ8LLL0OTJvz15GDq1bPr6BYuVBBNLZkywe6b/qDgw8tIOJaV3Z/ezbRTD3Nv2+GsDC4H7dvbzqSHDrldqtdc8frPMy7cM/Qcx6HdloVQtSqnduyhRY4feOnEAL6ZHsCAAQqiXvXoo3ah9KefwsCBblcjGZDCaCoYMWIE8fHx1K9fn/nz5/PPP/+wYMECGjRogOM4jBgxIsn5R44cYeXKlUlu27dvZ8aMGTz++OPMmDGDjRs3smHDBgYPHsysWbN48MEHL1lDfHw8Tz/9NGvWrGHu3Ln06NGD9u3b/2c09EIdO3bk8OHDtGzZkt9//52tW7cyb948OnTowNGjRwF4+eWXGThwIFOmTGHDhg107tyZ3bt3p3jNnDlz8vLLLxMeHs748ePZsmULy5Yt48MPPwSgYMGCBAYGMmfOHPbu3cvhw4eTvU737t0ZPHgwkydPZuPGjbz55pssXryYrl27XvK/g4hIurByJbRsCZUrs2fIFzRq6k/OnDaIlirldnHp266YWALL7KdI28VkLnSYg7Mqs2Z+Ax5q+K794/3776FSJdsQJp07u/4zOiYWh/PrPy8VSM/uGXpW9lMn+GDW/+g1ZRDbClTnhsMr2XjDfURGQjITysQb3nrLrh/t0QNSaCQpklr02VMqKFOmDJGRkbz99ts88cQT7Nu3jwIFCvDAAw/w1VdfUaxYsSTnL168mCpVqiQ59tBDDzFw4EBy5MhBt27d+Oeff8iUKROlS5dm8ODBvPzyy5esoXbt2tx8883UrVuXEydOnLvepRQtWpSlS5cSHh7Offfdx8mTJylRogT33nsvWbJkAaBr167s2bOHdu3aAfDEE0/w2GOPsW7duhSv269fP/LkyUPfvn3ZuXMnhQoV4sknnwQgU6ZMfPDBB7z99tv06dOHmjVrsmDBgv9c46WXXuLo0aO8+uqr7N27l5CQEKZOnUrlypUv+ZxERNK86Gi7DUOePMR+/T2NW2Xn0CFYvBhKlHC7uPSvaFAg0TGxZMp1kkKtfuPIsjLELL6JPXvy8vP0+6n76z12mmO9evDaa9Cnz7nmRhFR0Qyas4FdMbEUDQqke1hI0qY9acyl1n+m9LzOHh80ZwO5N6zhoxkDKXZoF5+V6cNTW3rxRBt/Ro2CizYgEG8yBsaNg23b4PHH7ZvLRX1FRFKL8eSavmsRGhrqREZGJnvfunXrUuxKKylr27YtBw4cYMaMGW6X4vP0GhMRn3bsGNSqBZs2kbh4KS36VmL6dPj2WzjTOFxS2cXdYAHYn4e4n6qxe0cA3btD3x7HyfxqZ/jkE9v16IsviDiW7T+PCwzwp1/zimk2kJbuMZPk/mo0wLb+DVN+YEICjBoF3btzOlc+HuMLvjtcmxEjoF07LVP0GXv2wO23Q3w8LFsGwWnzdSq+yRiz3HGc0IuPa5quiIiIL0pIsOu5/vwTvv6anl9WYto0GDJEQdSbmlUJpl/zigQHBWKA4KBAhnYqycY1AbRvb2fq3lk/O+u7joEpU2DLFqhShbXvDiX2dHySayXpIpsGJbv+8xLHAYiMhDvugJde4u+y9Sj170ois9fml1/sklsFUR9SuLCddn7kCDRpAhf0HhFJLQqjIiIivqhLF5gxA4YPZ+yu+xkwAJ5/3vYwEu9qViWYpT3uYVv/hiztcQ/NqgSTPTt8/DFERMDff8Ntt8FH+x/C+fMvqFaNnlMHM+LbAeQ6mbThYErdZdOCi9d/gh3t7R4W8t+T//0XXngBqlcn8e9/GH7HJEqtmUGVewuwfLlmgfqsSpXgyy/tOvXHHwftSCCpTGE0HZowYYKm6IqIpGXDh8MHH8ArrzC/3As89xzce689pJEk39K0KaxaZffEfP55aNqxGPsnz+PDsHaEbfqVH8a9yH0blsKZZVGXHEX0ccmNEv9n2nFiIkyYACEhOB9/zPoGL1I2YQNdIlvzzjuG77+Ha9xQQLylYUM7BSMiAnr2dLsaSefUwEhERMSXzJgBnTtDkyasf2YQD90NN90EX3+tLS88ITWaChUpAj/8YD8seO01qFjZnw6v96N18C28890QPoroR1SREIbUf4aHWrb20DM5z5uNkppVCU752n/9ZUdDly7lZJU7eaXMj3z0Y2Vuvx2+HQMVK6ZKSZIaXnrJ7gM7YACEhMBTT7ldkaRTPt/AqFy5chh9DCypwHEc1q9frwZGIuI7oqLsEFtICAemLeKOetk5ehR+/11buHhCcs2IPN1U6K+/4LHHYPVqaPjoMWJu+p27ln1Ht6VfUPDIAdsZuV8/uOUWj/y8a31OHg2wR45A797wwQc4efIwp+4AHvq+LX6Z/OjXz44Y+/tf9iria+Li4IEH7B5Sc+dC7dpuVyRpWJpsYBQQEEBsbNpdWyG+LTY2loCAALfLEBGxdu60QSVvXk5N+Z4HH8/Ozp22c66CqGdcamsST6lUCf74w67tnTk5B4en1qNRj2EU3P039O9vt8249VZ4+mn455/r/nnX8pyuZb/QZDmOXV9YrhwMHcr+pu2oW2QD93/zNPXv9WPtWujUSUE0zQoIgG++gTJloHlz2LzZ7YokHfLpMFqwYEGio6M5ceIEbo/gSvrhOA4nTpwgOjqaggULul2OiIjdwqVxYzh6FOf7GbR7syhLlsCnn9pGpOIZKTUP8nRToaxZYehQmD3b9vGpXRseeDgbK8Nes912O3eGSZPs/OvXXrMnXaNreU4eCeXr10ODBvDooyQULsrQVr9ROOIjNuzPy5Qpdrlh8eJXfjnxUUFBdumAMfbDsut4rYokx6dXn+TKlQuAXbt2ERcX53I1kp4EBARQqFChc68xERHXJCRAq1a2C86MGbzzXSUmToR33oGWLd0uLn0pGhRIdDIhLbWaCoWFwaZNMGKEnZlbpQo8+mg+3n77fcq++CK8+SYMGgRjxthGMZ062SR7Fa7lOV1XKN+50z6hIUMge3bWdhpF4+87sHWyP88+awd/g4KutHpJE8qUgWnToH59aNHCLpDWzDLxEJ9eMyoiIpLudekC//sfjBrF5KDnad0annzSNiRVywTP8saa0ZTExNjcOXQonD5t99h84w0osu9P6NHDDqOWKAF9+9pFp1c4t/VanlON/vOTDbDBQYEs7XHPfx9w5AhMnQqffw4LFoDjENuyDV3iBvLRtIKUKwejR9vlzpKOTZhgGxl16AAffaQ3KLkqaXLNqIiISLr26ac2iL74Ir/c+jxPPQW1atk/7PV3nudd0dYkqSQoCN591y6769DBDoaWKQM9v7qVmMk/wE8/QYEC0KaNHUIdMwa2b7/sda/lOV3RfqFxcXZ6ZqtWUKiQXeO6cycJb/RmSv/NFJs7gbHfF+Stt+yWlAqiGUDbtvaDk9GjYdgwt6uRdEIjoyIiIm74/Xe7mLBGDbaOms0dNQPInRt++w3y5XO7OEltW7bYWbpffAF58ti/8Tu9kEi2md/YIdNNm+yJZcvatZkNGkDduh6bA5tsN93KRW33pc8/t42JDhyA/PmhVSv+rvk4I/+ozmefG/buhRo1bCapUMEj5UhakZgIDz9sFwV/951dRypyBVIaGVUYFRER8bZduyA0FLJm5fC8P7izUT727LFB9Kab3C5OvGnlSujVC2bNgqJF4a234Km2DgGb19ntNObOtVNjjx8HPz+oVu18OL3jDsic+fqL2LrVNlSaOBE2brTrVps0IfbhJ/gqJoxPPg1g6VI7c7hhQ3jmGZtB/DS/LmM6ftxO4di4EZYutS2kRS5DYVRERMQXnDxpR0TXrCFhya80Cq/IvHk2c9Sp43Zx4pZFiyA8HH75BW68EZ591g6E3nor+CectiPpZ8PpsmV2hCp7dvuiadDANpcpV84GhaNH7TrPo0cvf9uwwf5QY6BOHZzHHiey5EOM/io3X35pGz3fdJMNoE88AUWKuP1fSnzCrl1Qvbr9RGLZMihc2O2KxMcpjIqIiLjNcWwDkE8/hWnT6PH7gwwYYHuBPPus28WJ2xzHLtN86y2IirLH8uSxg1B168I998DNN4PfkRj4+efz4fRq93/084OcOe2tcGF46CEOhD3GZz8XZ+xYWLsWsmWDRx6xIbRGDa1hlmSsWGEXC1eqZEfvs2RxuyLxYQqjIiIibvvf/2z33N69mXzTW7RuDc89Bx9+6HZh4muio+3f9z//bG9bt9rj+fPbwdC6de2tXDkwf2+HefPstitnQ2ZKt1y5cLIGEnvScOgQ/PknjBtnl//Fx8Ptt9sA2rIlaPczuaypU+0a0rZt7QtJn1pIChRGRURE3DR3Ltx3HzRrRlTPb6hR04+qVW0TVU8s+5P0bccOG0rnz7df//nHHi9c+Hw4LVMG/v0XDh2Cgwft1wtvFx47der8tfPnt1Nwn3nGjryKXJXevaFPH7v37CuvuF2N+CiFUREREbds3mzXVxUrxv5vfyG0Tg4SEyEy0u6aIXI1HMeOlJ4dNf35Z9i9+7/nBQZC3rznb/nyJf0+b14oVgzq1dMHInIdEhOhRQvbYXfWLAgLc7si8UEKoyIiIm44csR2Pd27l7hfI6nfvjTLlsGSJVC1qtvFSXrgOLax6Z49SYNmYKDblUmGceyYXVz899+2oZHagstFUgqjmdwoRkREJENITITHH7dJ4ccf6TysNIsW2V00FETFU4yBkBB7E3FFjhx24XG1atCkid2nykN74kr6ph2iREREUstbb8H338PQoYzZcg+jRkG3btC6tduFiYh4WMmStqHRli3w6KOQkOB2RZIGXFEYNcbUMsZ8Z4yJNsY4xpi2F91vjDG9jTG7jDGxxpgFxhgtgRcRkYzrm2/gnXfgmWdYWrkjHTvCvfdC//7nT4mIiqZG//mU7jGTGv3nExEV7V69IiLXq2ZNGDUKZs+GHj3crkbSgCsdGc0BrAZeBmKTuf9VoCvwIlAN2AfMNcbk9ESRIiIiacrKlXarg7vuYmf4SB562FCyJHz5Jfj721MioqIJn7aK6JhYHCA6JpbwaasUSEUkbWvfHjp1gsGD7Z7KIpdwRWHUcZxZjuP0dBxnCpB44X3GGAN0Bvo7jjPVcZzVQBsgJ6CJSCIikrHs3w9Nm0LevMROnMqDrbJw/LhtNJknz/nTBs3ZQGxc0mlssXEJDJqzwbv1ioh42pAhcM890KGDXT8qkgJPrBktDRQGfjx7wHGcWGARcJcHri8iIpI2nD5tN4Dftw9negTPvlWYyEiYOPG/+zfuikluolHKx0VE0oyAAPj6a7t30IMPws6dblckPsoTYbTwma97Lzq+94L7kjDGdDDGRBpjIvfv3++BEkRERFzmOPDii7BoEYwdy9DFVfn8c7sXfNOm/z29aFDy+26kdFxEJE3Jl8922D12DJo1g1h90Cb/5cluuhdvWGqSOWZPdJzRjuOEOo4TWqBAAQ+WICIi4pKBA2H0aAgPZ26B1nTrZgcEXn89+dO7h4UQGOCf5FhggD/dw7Q/h4ikEzffDF98AStWwDPP2A/tRC7giTC658zXi0dBC/Lf0VIREZH058svbefIVq3Y8tQ7tGwJ5cvb3h1+KfymbVYlmH7NKxIcFIgBgoMC6de8Is2qBHu1dBGRVNW4Mbz7LkyeDAMGuF2N+JhMHrjGNmwgbQD8AWCMyQrUBLp74PoiIiK+a9EiaNMGatXi6PAJNK1j0+e330LOy/SUb1YlWOFTRNK/Hj3gr7+gZ087Wtq4sdsViY+4ojBqjMkBlD3zrR9QwhhTGTjkOM4OY8xQoJcxZj2wEXgdOAZ84fGKRUREfMX69XYtVOnSJEyZzhPtsrBuHcyZA2XKuF2ciIiPMAbGjoVNm6B1a9th9+KubpIhXek03VAg6swtEOhz5t9vn7l/IDAEGAlEAkWAex3HOerRakVERHzF3r3wwAO2a+QPP9D13bx8+y38739Qv77bxYmI+Jhs2eweVzlyQJMmcPCg2xWJD7jSfUYXOI5jkrm1PXO/4zhOb8dxijiOk9VxnNpn9hsVERFJf44ft9PM9uyBGTMY9l1phg2Dl1+Gl15yuzgRER9VrBhMn263ennkEYiLc7sicZknu+mKiIikfwkJdprZ8uXw5ZdERFfjlVfsbN3333e7OBERH3fHHbbz+Pz50LWr29WIyzzRwEhERCRjcBzo3NnunTd8OMsKN6F1HahWDSZNAn//y11ARERo0wZWrbKf4N1yC3To4HZF4hKFURERkSv1v//BiBHQtSvbGnai8R1QuLDNptmyuV2cpAcRUdEMmrOBXTGxFA0KpHtYiDouS/o0YACsXQvPPw9BQXbarmQ4CqMiIiJXYupU6NYNHn6YQz0Gcv/ddrnTrFlQqJDbxUl6EBEVTfi0VcTGJQAQHRNL+LRVAAqkkv74+8M338B998Fjj0FgoLZ8yYC0ZlRERORyfvkFHn8c7ryTU2M+o/nDfmzbZhtDlivndnGSXgyas+FcED0rNi6BQXM2uFSRSCrLnh1mzoTKlaFFC5g3z+2KxMsURkVERC5l0ya7DUHx4jgR3/J0x0AWLoTx46FWrfOnRURFU6P/fEr3mEmN/vOJiIp2r2ZJk3bFxF7VcZF0IVcuuznzTTdB06awZInbFYkXKYyKiIik5MABu5eoMTBrFm9+kJ8vvoB337UNdc86O70yOiYWh/PTKxVI5WoUDQq8quMi6UbevDB3LhQvDg0bQmSk2xWJlyiMioiIJCc21o6I7twJ333HuEVleecdaNcOwsOTnqrpleIJ3cNCCAxI2pI5MMCf7mEhLlUk4kWFCtlpunnzQliY7bYr6Z7CqIiIyMXi4mxDjd9+g0mT+PHonXToAPfeC6NG2YHSC2l6pXhCsyrB9GtekeCgQAwQHBRIv+YV1bxIMo5ixeCnn2wzowYNYONGtyuSVKZuuiIiIheKjbWNNGbOhGHD+Ktscx6+GypUsI0fAwL++5CiQYFEJxM8Nb1SrlazKsEKn5Kx3XCDHSGtVQvq1YPFi6FUKberklSikVEREZGzjh61a0RnzYIPPyT6oZdo2BBy5rSHcuVK/mGaXiki4kHlytlAevy4DaTRWn+fXimMioiIABw8eP5T+IkTOfrYczRqBDExdpC0WLGUH6rplSIiHlapEsyeDfv3Q/36sG+f2xVJKtA0XRHJ8CKiohk0ZwO7YmIpGhRI97AQhYiMZvduuyB00yaYPp2TDRrTopntnzFjht0C73I0vVJExMOqV7efBoaF2ffon3+GPHncrko8SCOjIpKhaUsOYft2qFkTtm2DWbM4WqcxDzxgt737+GO47z63CxQRycBq1oSICFi3Du6/3y6nkHRDYVREMjRtyZHBrVsHd98Nhw7BTz9xqPI9NGgACxc53PjIGt7ZNJMa/efrwwkRETfdey98/bXdf7RRIzhxwu2KxEMURkUkQ9OWHBnYihW2W2N8PCxYwJ6St1OnDixf4VCkeRSnS2/XaLmIiK9o2hQmTrTr+ps3h1On3K5IPEBhVEQytJS23tCWHKkrIiqaGv3nU7qHSyOPS5ZA3bqQLRssXszfuStRsyZs2QI3PbmSTDfsTnK6RstFRHxAq1bwySd2HUWrVnZPaEnTFEZFJEPTlhze5/o63Tlz7JSvIkVgyRI2JN7I3XfDgQN2J4ET+Xcl+zCNlouI+ICnn4bhw+060jZtICHhsg8R36UwKiIZmrbk8D5X1+lOnQqNG0NICCxaxMqDxalZE06fhgUL4M47NVouIuLzOnWC/v1h8mRo186+iUuapK1dRCTD05Yc3uXaOt0JE+CZZ+COO2DmTH5ZG8QDD0CuXHZE9Kab7Gndw0IIn7YqSWDWaLmIiI957TU4eRJ694a//rLB9OwbuaQZGhkVERGvcmXk8YMP4KmnoF49+PFH5v4RRIMGULCgXT564d8vGi0XEUkj3nrLTtfdvh2qVIFx48Bx3K5KroJxXP4fLDQ01ImMjHS1BhERXxYRFc2gORvYFRNL0aBAuoeFXDYYXctjvOXsmtGLRx5TJfDt3QsvvgjffAMPPgiTJxPxQxZatoRy5eDHH6FQIc/+SBER8bLoaHjySZg/H1q0sJtE58njdlVyAWPMcsdxQi8+rmm6IiJedLUh8eLgdrbZD5Di467lMdfjap/T2ftSNSw7jt0CoHNnOHYM3nkHXnuNzydn4qmnoFo1mDVLf6uIiKQLwcEwdy4MHgy9esFvv8GkSVCzptuVyWVoZFRExEsioqIJn/oX+Q7uJiAhnqNZshGXIxd9HqmaYhCr0X8+0cmspQwOCmRpj3s89phr5dVRziu1Ywc89xz88IPtSDR2LJQvz6hR0LEj3HMPfPst5MjhTnkiIpKK/vgDWreGrVttMH3zTcik8Te3aWRURMQNCQmwahUsWULO0VNZsO0vCh07lOSUU4MDIG8eyJ3b3oKCzv37qb8OciRrDo5mycbO3IX4o1gFYgJzXbLZjzcbBF2qM66nw+hlR2ATE+Gjj2xTC8ex60RfeAHHz5/+/aBnT2jSBL76CrJm9WhpIiLiK6pVgxUr4KWXoG9f26Fu0iQoXdrtyiQZCqMiIp506pT9VHbxYnv75Rc4fBiAcjkL8GuJikQGV+BolmzkOnWcnKdOkOvUcZ67NZ897/BhiImx618OH6b1gUNkO30yyY9Yn78ka8reCl8ft1OQihRJcn/RoMBkR0ZTo0GQt4LvZaceb9xo2/svXgz168OYMVCqFDt2wIOtT7JiaVayV4hmX/WNzF53k8+snxURkVSQMyeMHw9hYfDss1C5Mnz4oR0xFZ+iMCoicj0OH7aB82z4/OMPG0gBypeHli2hZk1OVK1J2NCDRG8M5PSO3IDBL9tp/ANPkbegQ557QihYEAoUsLd8+eysoh+jonnjmyj8jx+l7MF/qP7PGu7auZqmf86Dlt/Zn3PjjVCr1rlb93tvInz6aq9sTeKt4JvSCOyQWWtpNneS7aiYNavtpNi2LYmO4cOR0O3VRE7FZSJPvTXkrLqd3UdJ1fWzIiLiQ1q1stt5PfaYvc2eDSNG2D29xCdozaiIyLU4ftw2xXn/fYiLs8nxttugZk0Sa9RkU8Ea/LIxP7/9Br//DqtX2xm7AJlynwD/RBJjM5MYGwCY/1zeGNtcp2BB8As8xZ7Th4nPcYRitxzlzWcK0rJaIYiKgkWL7G3xYvj3X/vg4sX555ZQvshaijn5y3GqdBm631cu1ZoXeWPNaOkeM7n4t1X5fVsZ8MMHVNqz2XbKHTkSihRh/Xpo395u2RJ040Gy1/uTTLmTBubUWD8rIiI+Kj7e/s7u2xdKlbJ7klav7nZVGUpKa0YVRkVErobj2O43L79sG+U88QT/NmnDr84d/PJndn77zQ6OHjliT8+d2/6+u+MOuP12OBC4m9HL1p1b99ilXgg1Swazfz/s2wf79//3dvb4pk0292bJAjVq2Nmo9evbDOxvEm3iPRtOFy2y25qAXSdz//3wwANQty5ky+bR/yTe2EbmwqZMmePj6PTLlzz/+xSOZstJ3nGj4eGHiYs3DBwIb78N2bPD//4Hb62dmVzWxwDb+jf0aI0iIuLjliyxI6S7dtlfFq++Cv7+bleVISiMiohcry1bbEOEWbNwKlZkTpNRvPT13WzaZO/294dKlWzoPBs+b7oJ/Pw88+OPHbMDoPPm2dtff9njQUG2Q+zZcFq2LBgcm15/+sl2lf3pJzhxwibZOnXOh9Mbb/RMcaksIiqafpN+ocb6X3nut6ncdHAHERXrETBsKA3r3kJkJDzzjP1v0qIFDB9u9w/1ZmdhERFJA2Ji7DrSr7+2H9B+9hkUK+Z2VemewqiIyLU6eRIGDIB+/SAggH/av82jSzuxdFkAoaF2Wejtt0PVqh4fdLykvXvt/t7z5tnt1f75xx4vUcKG0nr17K1QoTPPYfFiG0xnzYING+zJZcqcD6Z16kCg55scXZcDB+xI9JQpJM77Cb/4OP4OKswHzV6i5ktPcm9IML1729nShQrBqFHQrNn5h/vk1jMiIuIux4EJE+DFFyEgAB59FB5+2PZeSIvbwMTH+3zdCqMiItfihx/sL6stWzjZrBW9sr7P/74qSsGC0L8/PPmk50Y+r4fjwObN50dN58+3H/4aY4Ny06b2Vq6cPca2beeD6fz5EBtrGwDVrQv33WcfVLGid9P1WXv3wvTpMGUKLFhgF9uWLm3/UHj4Ydu23xh+/tmuDd2yxX4dONCOEl/MG9OIRUQkDdq4EV5/HWbOtLOH8ue3n2g+/LCdchQQ4HaFyfv3X1i69PyynMOHYd06t6u6JIVREZGrsWMHdO4M06fj3BTCtHojeXpSPU6csMtF33jDrgf1VQkJdpu1H36A776D5cvt8bJlzwfTu+46s1Tm5ElYuNCe/MMP9pcz2JQdEgJVqti2+GdvBQp4vuDoaJg2zQbQxYttur7ppvMBtHLlMynahuzu3eGTT+zA7pgxNkOLiIhckxMnbKfdKVPg++/tupg8eewvy4ceggYN7DIXt+zbl7QnxF9/2d+TmTPbxhS1akGfPj49OqowKiJyJU6fhiFDbMc9YFOrN3j4ly78tT4zYWEwdKgdXUxrdu60ofTbb+Hnn20jpPz5oVEj+7u2QQPb9AewQXzFCli50nbsXbnSHjsrONiGwwtDaunS/x0ijo+329ycPJn064X/Xr7c/vL/9Vf7mJtvPh9Ab775XAAFWL/eDpgOH24HT7t0sb973Ri8FRGRdOrkSbv2ZcoU+0vz8GG7FUzjxvZ3U1hY6i9p+eefpOFz/Xp7PFs2+0ny2e3cqlf3veU1KVAYFRG5nPnzoWNHWL+e4/c+SGeG8smPJbjhBhtCGzVKko3SrCNH7AfA335rZ+nGxNgZuvXr22DauPGZdaYXOnTIhtILA+q6def3q8mZ094uDJ4JCVyRypXtL/iHHkqS9B3HZtXp0+3t7AykGjXs/x6h//mVJiIi4kGnT9sGgFOn2l9Ehw7ZT24bNrS/tx544IJPcq/R2XU2F4bP7dvtfblzw913nw+fVav67tThy1AYFRFJiePYebfvvkti6RuYWH04HSIeIFMm6NULXnnFhrX0KC7Ozor99lt7+/tve/zGG+3SzLO3KlWSGYE8edJuJ7NyJfz5p/0+a1Y7lSmlrxcfK1kSbrjh3CXj423n/enTISLCDsj6+0Pt2nYr0WbN1PRQRERcEBdnl7RMnWqXlezbZ0clw8LsJ7inTtnwei1f4+LszyhQ4HzwrFXL9m5IJ1vPKIyKiCTHcaBHDxg4kC11n+G+jcPZHB3IY4/ZBrrBGajPjePYZSizZsGyZXa/1Ohoe5+/v501ezacVq8Ot9zimQ9oT560TZemT7eB+OBBm1PvvdcG0MaNIV++6/85IiIiHpGQYD85nTIFZsw4v3Va5szJf73UfZkzQ6lS9lPXkJD0MQUrGQqjIiIXcxy74fXgwcwv9zz114+gchU/PvjAzooR2L3bhtILb4cO2fuyZrUzbM8G1Lx57cjmhbeEhP8eu/D25582/B47ZpfkNGpkA+h990GOHK4+dREREfEQhVERkQs5DnTtCv/7Hz/e2JGwTcN5/XVD797pZkZMqnAcuyvMheF0+XI4fvzarlewoJ16++CDtot+5sweLVdERER8QEph1Hf7/4qIpBbHsQtBhw1jZpmXaLRpKO+9ZwgPd7sw32eMXeJ5ww3QsqU9lpAAGzbYQJop0/mbv3/S75O7L3PmdDsjSURERC5DYVREMhbHgZdeghEjiLjhFR7c8j5DhhheecXtwnxDRFQ0g+ZsYFdMLEWDAukeFkKzKpdeOOvvDxUqeKlAERERSTcURkUk40hMhE6d4MMPmVKqGy22DmTkSMMLL7hdmG+IiIomfNoqYuPslizRMbGET1sFcNlAKiIiInK1/C5/iohIOpCYCM8/Dx9+yOQSr/LI9oF88omC6IUGzdlwLoieFRuXwKA5G1yqSERERNIzjYyKSPqXmAjPPguffMLnxcJp+8+7fPa54fHHr++y1zKl1Zftiom9quMiIiIi10MjoyKSviUmQvv28MknjCv6Ok/tepfJX3omiIZPW0V0TCwO56e0RkRFe6RsNxQNCryq4yIiIiLXQ2FURNKvhAR4+mkYN46PC7/Fc/veZspUwyOPXP+l0+OU1u5hIQQGJN3XJjDAn+5hIS5VJCIiIumZx8KoMSanMWaoMeZvY0ysMeYXY0w1T11fROSqJCTAU0/Bp58ysmAfXv63NxHfGpo188zl0+OU1mZVgunXvCLBQYEYIDgokH7NK6bpqcciIiLiuzy5ZvQToBLQBtgJPA7MM8ZUcBwn7c5bE5E0JSIqmiGz1tLli3dptnYhA4L60Ofom8yYAfXre+7nFA0KJDqZ4JnWp7Q2qxKs8CkiIiJe4ZGRUWNMIPAQ0MNxnAWO42x2HKc3sBl43hM/Q0TkciKiogmf+hcdvxpEs7UL6RnYh54netLrg/0eDaKgKa0iIiIi18tT03QzAf7AyYuOxwJ3e+hniIhc0qA5G2i0YjYtV83lvSyv0j8hnIKP/M6sA6s8/rM0pVVERETk+nhkmq7jOEeNMb8CrxtjVgN7gEeBO7GjoyIiqS7npnX0/fEj5meqzVv0plCr38lS5DC7YlLn52lKq4iIiMi182Q33SeAROx60VPAS8BkIOHiE40xHYwxkcaYyP3793uwBBHJsI4d46PvB3KY3LSOn0y+5ivJUuQwkPbXcYqIiIikRx4Lo47jbHEcpzaQAyjuOE51IADYlsy5ox3HCXUcJ7RAgQKeKkFEMirHgeeeo8SBnbRK+IpTtQ6TtcQhQOs4RURERHyVx/cZdRznuOM4u40xeYAw4FtP/wwRkSTGjoVJk3jb9Obo3XdSPmyX1nGKiIiI+DiPbe1ijAnDhtv1QFlgELABGO+pnyEi8h9//YXz4ossyVqfTwv0JDLCn3z57nG7KhERERG5DE/uM5ob6AcUAw4BU4FejuPEefBniIicd/QoTosWHHLy0CphEtOn+JMvn9tFiYiIiMiV8FgYdRzna+BrT11PROSSHAeefRZn02aaO/MJH16Q6tXdLkpERERErpQnR0ZFRLxn9GiYPJk3zTsUfqQ2HTu6XZCIiIiIXA2FURFJe6KicF5+mQWZw5hSMpw/PgFj3C5KRERERK6GwqiIpC1HjuA88ggHnHy08fucmVP9yJnT7aJERERE5GopjIpI2uE40L49iVu20dz5mb4TClCxottFiYiIiMi1UBgVkbTjww/h66/pRT9CnqlJmzZuFyQiIiIi10phVETShhUrcF55hXkB9zOn/Kv8MtztgkRERETkeiiMiojvO3wY5+EW7KMg7TN/xtwpfgQGul2UiIiIiFwPhVER8W2OA888Q+L2v2nuLOT9L/Jz441uFyUiIiIi10thVER828iRMHUqPRhI9c41eOghtwsSEREREU9QGBUR37VhA4ldu/Gjf0OWhnZlwQC3CxIRERERT/FzuwARkWQlJrL/0TYcjgvk6YBRnK7zK7PWRLtdlYiIiIh4iEZGRcQnRfUZQpWo32nHGOIb7uUAMYRPWwVAsyrBLlcnIiIiItdLI6Mi4nv27KHMgLdZSE2+LBdGtrL7AIiNS2DQnA0uFyciIiIinqAwKiI+x3m5M1lOneK5zCPJU399kvt2xcS6VJWIiIiIeJLCqIj4lpkzMV9/xTu8zv4GCfhnP53k7qJB2mBUREREJD1QGBUR33HsGPEdnmedXwV+qP4y+SvvSXJ3YIA/3cNCXCpORERERDxJYVREfIbT63X8du2kY8AYvvkyF/0fqkhwUCAGCA4KpF/zimpeJCIiIpJOqJuuiPiGZctg+Ad8yPM0eu8uSpeG0gQrfIqIiIikUwqjIuK+uDjin27PflOEKZXfY97LbhckIiIiIqlNYVREvCIiKppBczawKyaWokGBdA8LOT/qOWQImdb8RSe/6Qwdnxt/f3drFREREZHUpzAqIqkuIiqa8GmriI1LACA6JpbwaasAaJYzloQ3e/MtD1KuRzNuvdXNSkVERETEW9TASERS3aA5G84F0bNi4xIYNHs9Ce2f40RcZv5XejhvvOFSgSIiIiLidRoZFZFUtysmNtnjdyyZif+Cn3iVUbw7IZisWb1bl4iIiIi4RyOjIpLqigYF/udY3hOHeeOnsSzlLujwLLVquVCYiIiIiLhGYVREUl33sBACA5J2Jeo9/xOynYrl9QKj6T9Qb0UiIiIiGY2m6YpIqjvbNfdsN91m+1bTZM3PvM0bdB5zM7lzu1ygiIiIiHidwqiIeEWzKsE2lB4/TlxIJzaYENY368mbTd2uTERERETcoDAqIl7lvNWbgOjtdMmxkLGj1LFIREREJKNSGBUR71mxAmfIEMbQnoeG1aJwYbcLEhERERG3KIyKiHfEx3O6bXv+pQAzaw4g4im3CxIRERERNymMioh3fPQRmVetoFvAl7w/Lg/GuF2QiIiIiLhJYVREUt/+/Zzu8QaLqEfFvo9QtqzbBYmIiIiI2xRGRSTVxb/aE3P8GB+U+YCpXTQkKiIiIiKgneZFJHVFRuI/YSzDeZGuYysQEOB2QSIiIiLiCxRGRST1JCZysn0n9lKQNQ+9Re3abhckIiIiIr5C03RFJPV89hlZV/7OW1km8Paw3G5XIyIiIiI+RGFURFLH4cOceuU1VnAHZfs8QXCw2wWJiIiIiC9RGBWRVBH/em8CYvYzuNQsJr+iFQEiIiIikpT+QhQRz1uzBr9RwxlDe174pCqZM7tdkIiIiIj4GoVREfEsxyG2w0vEJOZiWZN3qVfP7YJERERExBdpmq6IeNbUqQT+Mp+eASPoPSK/29WIiIiIiI9SGBURzzlxgtgXurCRShR681mKF3e7IBERERHxVQqjIuIx8X37Ebj/HwYUm8j47np7EREREZGUac2oiHjG1q0weBCTaE2bT2qRJYvbBYmIiIiIL1MYFRGPOPHcK5yMz8TP9w0kLMztakRERETE12kenYhcv9mzyTb3O17P1J83Pgp2uxoRERERSQMURkXk+pw6xfF2LxHNjWTv1ZmSJd0uSERERETSAo9M0zXG+Btj+hpjthljTp75+o4xRmFXJJ2Lf38o2aM30b/wMF7poYWiIiIiInJlPBUWXwM6Am2AVUAl4FPgFNDXQz9DRHxNdDSJffryLU146JP7yZrV7YJEREREJK3wVBi9C/jecZzvz3y/3RjzHXC7h64vIj4iIiqaQXM2sCsmllHfD6Pu6Xhm3PM/xjR0uzIRERERSUs81U13CVDXGFMOwBhTAbgHmOWh64uID4iIiiZ82iqiY2IJ/Wc196+dy2DTlaqds7ldmoiIiIikMZ4KowOAz4G1xpg4YA3wqeM4o5I72RjTwRgTaYyJ3L9/v4dKEJHUNmjOBmLjEvBPTKD3rE/YQXFG3v4In69Z63ZpIiIiIpLGeCqMtgSeBFoDt5359wvGmGeSO9lxnNGO44Q6jhNaoEABD5UgIqltV0wsAK1XzObmmM10z/YeAXftOndcRERERORKeWrN6CBgsOM4X575fpUxpiQQDoz10M8QEZcVDQrkZPRuui6cxDzq8XNYCNkC9lE0KNDt0kREREQkjfFUGM0GJFx0LAHPjbyKiA/oHhZC3OODyB5/gu7FehN44z4CA/zpHhbidmkiIiIiksZ4Kox+D/QwxmzDrhetAnQBPvPQ9UXEBzQ7uQPWzmag6cb+hqcpmSeQ7mEhNKsS7HZpIiIiIpLGeCqMvojdT3QUUBDYDYwB3vbQ9UXEbQkJHH3iBQ4TTEL4W+x8N4fbFYmIiIhIGuaRMOo4zlGg85mbiKRDcSM+IueWlfQq9BUD31AQFREREZHr46mRURFJz/btI77H6yykHg0ntCBrVrcLEhEREZG0Tg2GROSyjnbqQaaTx/i+wXDC7jNulyMiIiIi6YDCqIhc2q+/kvOb8QwP6EK3seXdrkZERERE0glN0xWRlCUkcPjxjhwlGL8336B4cbcLEhEREZH0QmFURFJ0evjH5N4aRZ/grxjwmpoWiYiIiIjnaJquiCRv3z4Swnsxj3o0ndiCgAC3CxIRERGR9ERhVESSdfh527RoXpPh1K6jpkUiIiIi4lkKoyLyH84vv5J72nhGZu7CK6PVtEhEREREPE9rRkUkqYQEYh7ryHGCCXzvDQoVcrsgEREREUmPNDIqIkmcHPYxebZHMbzUENp1VtMiEREREUkdCqMict7+/ST27MVP3EPzyS3w93e7IBERERFJrxRGReScQ8/2IODUMRa1GMHtd6hpkYiIiIikHoVREQFs06K808fxYdYuvPShmhaJiIiISOpSAyMRgYQEDrXuSCzB5B78BvnyuV2QiIiIiKR3GhkVEU7872Py/R3FRzcO4Ynn1bRIRERERFKfwqhIRrd/P/SyTYse/qoFfnpXEBEREREv0J+dIhncgae6E3D6GMueGEHlKmpaJCIiIiLeoTAqkoHFf/8D+Wd+yqjsr/LCcDUtEhERERHvUQMjkYzq8GFOPNaef6hA6fFvkju32wWJiIiISEaikVGRDOpAm65kP7qbr+6bQJMWWdwuR0REREQyGI2MiqQTEVHRDJqzgV0xsRQNCqR7WAjNqgQne+7pGT+S/9uxjMj+Gq98Uc3LlYqIiIiIKIyKpAsRUdGET1tFbFwCANExsYRPWwXw30B65AgnHmvHFspR5vPe5Mnj7WpFRERERDRNVyRdGDRnw7kgelZsXAKD5mz4z7n72nQn55FopjYcz/0PZvVWiSIiIiIiSSiMiqQDu2Jir+j46VnzKBgxmjE5u/LipDu8UZqIiIiISLIURkXSgaJBgZc/fvQox1s9w3pCuPGLPuqeKyIiIiKuUhgVSQe6h4UQGOCf5FhggD/dw0LOfb/nyVfJffQfvm06nnqNkg+vIiIiIiLeogZGIunA2SZFKXXTPTnzJwpHfMSYXF154fM73SxVRERERAQA4ziOqwWEhoY6kZGRrtYgkq4dO8bB4IocOJKZPT+spPZ9GhUVEREREe8xxix3HCf04uMaGRVJ56Kf6EGRI3/z2UOLeeWiIHo1e5OKiIiIiHiSwqhIOhb7wwKCI0YyLvcrdPi0RpL7rmpvUhERERERD1MDI5H06vhxjrd8mk2UpfzUd8iePendV7M3qYiIiIiIpymMiqRTOx4LJ+/R7fzYchx31sv2n/uvdG9SEREREZHUoDAqkg4dm7mQEt8OZ2KeF3l6fM1kz7mivUlFRERERFKJwqhIenP8OCcefYYt3ED5ae8RmEK2vJK9SUVEREREUosaGImkM9se60Xpo1v49rEFtK+TPcXzLrc3qYiIiIhIalIYFUlHjvywhJLffsAXeTvx5Njalz2/WZVghU8RERERcYXCqEg64cQcJvaRthygFBW+7UeWLG5XJCIiIiKSMq0ZFUkPEhP5++7HyHvsbxY89RmV787hdkUiIiIiIpekMCqSDmx78i1KrZnJJxWH0faTu90uR0RERETkshRGRdK4PSOnUnrSO0wNeobHFj+Pn/5fLSIiIiJpgP5sFUnDjv22mlwvtiHS/3aq/DKSXLmN2yWJiIiIiFwRNTASSaMSDvzLkXrNME5OTn0xlRvKq2ORiIiIiKQdGhkVSYsSEthc/VHyn9jBr12nUuMRbc8iIiIiImmLwqhIGrSmWS9Cts3h61ojaT74LrfLERERERG5agqjImnM5ne/4uYZA/i2yHO0nNfe7XJERERERK6JwqhIGrJv7p8UfeNp/shSgxqRwwgIcLsiEREREZFrozAqkkbE7jxIfKNmxBBEjh+mkL9oZrdLEhERERG5Zh4Jo8aY7cYYJ5nbTE9cXySjc+Li2RLaknynd7FpwHTK1y3sdkkiIiIiItfFU1u7VAP8L/i+CLAc+NpD1xfJ0P645zWq7/2J7x8cR+Pu1d0uR0RERETkunkkjDqOs//C740xzwBHgG88cX2RjGxF10lUXzKE2WU70WjqU26XIyIiIiLiER5fM2qMMcAzwETHcU54+voiGcmWb1ZQfkg7lueoRa3IIRjjdkUiIiIiIp6RGg2MGgClgU9SOsEY08EYE2mMidy/f39Kp4lkaIfW7yPw0WYc9CtAkcXfkC23WueKiIiISPqRGmG0PfCH4zgrUzrBcZzRjuOEOo4TWqBAgVQoQSRtO33sNDvufIQ8Cfv5d+x0ilYu6HZJIiIiIiIe5dEwaowpCDQFxnjyuiIZyZHdx/mrZGMqxywkssMYKrat6nZJIiIiIiIe56luume1BU4BX3r4uiIZwp61h9hbrRFVTvzO0qfHUvPjx90uSUREREQkVXhsZPRM46J2wJeO4xz11HVFMorNi3cTU7k25U4s5+t243j1xlKU7jGTGv3nExEV7XZ5IiIiIiIe5cmR0TrAjYCGckSu0opvtpC3VQMKO/v5tsck3s6Ug9iYWACiY2IJn7YKgGZVgt0sU0RERETEYzw2Muo4zs+O4xjHcZZ56poiGcGCD/6i6CN3k4sjxEydz7DceYmNS0hyTmxcAoPmbHCpQhERERERz0uNbroicoW+7fELt75cGxPgj7NgEcUerMauMyOiF0vpuIiIiIhIWqQwKuICx4GJj8+mwYD6nMhWgBwrl5KvZgUAigYFJvuYlI6LiIiIiKRFCqMiXhYfD6PrfcUjk5qwP28IBTcsJnuFkufu7x4WQmCAf5LHBAb40z0sxNulioiIiIikGk9v7SIil3D8OEy48yOeX/UCO4rfTcm/vscE5U5yztkmRYPmbGBXTCxFgwLpHhai5kUiIiIikq4ojIp4yYH9Dt/c1o+OO3uxvWIjSv3+NQQmP/W2WZVghU8RERERSdc0TVfEC7Ztdfjupm48v7MX/9R6jFLLp6UYREVEREREMgKNjEqaFxEV7dNTWn+eG8/eJu15+uQEdj38IsW/Ggp++hxIRERERDI2hVFJ0yKiogmfturcvpzRMbGET1sFcMlA6o0Au3UrfPhMJI8seJ66RLK/Y2+KDn8TjPHozxERERERSYs0PCNp2qA5G84F0bNi4xIYNGdDio85G2CjY2JxOB9gI6KiPVLTkSPQu3MMc27sxIAF1SmfYyenP/uSAiPeUhAVERERETlDYVTStF0xsVd1HK4twF6JhAQY+4lDj+KTeG5YOZ51PmRd06do3vVjQtbkoEb/+R4LvCIiIiIiaZ3CqKRpRYOSbwKU0nG4tgB7OYsXwyOV1lO6fT1GHXmcHBVKsPDzGTxcsQUbT/qnygisiIiIiEhapjAqaVr3sBACA/yTHAsM8Kd7WEiKj7mWAJuSv/+GJx46waJavZi8thJ3Z4vCGfUhOf76ldf/yZIqI7AiIiIiIumBGhhJmna26dDVNCPqHhaSpOkRXD7AXuzYMRgwANYMmMGQ+BcpxXbiH32CTP8bBIUKAakzAisiIiIikl4ojEqa16xK8FV1wr2WAHtWYiJMmgQfdNtBz30v05cI4sqWh08WkKl27STnFg0KJDqZ4HktI7AiIiIiIumNwqhkSFcTYBMT4fffYfp0+H5aHI22DGWRX28yZ3GgT38CXnkFMmf+z+M8MQIrIiIiIpJeKYyKJOPUKfj5ZxtAv41wKLzvT570m8jigMnkZxdOo6aYD4ZByZIpXuN6RmBFRERERNI7hVGRM44cgR9+gIgImDkT8hz9m6cyf8EfWSdSnLU4fpkw994Pz3+Cuf/+K7rm1U4hFhERERHJKBRGJVVEREWniRHBPXvgu+9sAP3pJ8h++hBP5/iG5dkncuPRJXAaqH43PPYhpkULyJfP7ZJFRERERNIFhVHxuIio6CRrJc/urwm4FkgTEmDbNli3DtautV9XrYKoKMjixPJMwRkMLD6JCn/Pwu9YHBQvDy+9C61bQ6lSrtQsIiIiIpKeKYyKxw2asyHF/TVTO4yePg2bNp0PnGe/bthg14GeVaHwIZoU+p1PKn9DpU1T8d93BPyLwMsvwWOPQeXKYEyq1ioiIiIikpEpjIrHpfb+mo4De/fCxo02eG7caG/r1sHmzXYUFGyWLFUKqoScoEOFFYTyB2UO/kGezcvw37YF9gA5c8LDD8Hjj0OdOuDv75EaRURERETk0hRGxeM8tb/mv/8mDZwXfj169Px5mTNDmTJw883Qsnkcd+RYzc3Hl1E0+g8yRf0Bc9ecT6jFi0O1atChnf16110QqH0/RURERES8TWFUPO5a9tc8eRIWL4Yff4RffrGh88CB8/f7+dlRzptugho17NeQMvGUZx1F96zAL2o5/PEHzFppLwaQNy9Urw5Nm9rgWa0aFC6cOk9aRERERESuisKoeNyV7K/pOHZa7Zw5NoAuXAixsXaUs3p1aN7cBs4bb7RfSxc5SZZNq2HFCtt16PMV8Ndf54Nn9uxw223wwgv2AtWqQenSWvcpIiIiIuKjFEYlVSS3v+ahQzBv3vkAunOnPR4SAu3bQ66yB5j772oOxhyE2F3U33eYW9dstQF0zRqIj7cPyJ3bBs+OHe3X226zqVXrPUVERERE0gyFUUk18fHw++82fM6ZY2fROo7NkvXrw5tvwr33QsmSMOeHP1g+dCzD1i4hNHod/k4iAKfy5CNL9VB44IHzwdOlEc+0sneqiIiIiEhaoDAqHpeYCJMnwxtv2L09/fzszNk334SwMDuDNlMmYMsW+GoqTJ1K2LJlhAHrCpRi1B0t+LPITawuVAb/4sVYGl7P7afkk3unioiIiIikZQqj4jGOAz/8AOHhdjln5co2lIaFQZ48Z05atw76T4UpU+DPP+2xqlUZULsNs2+6i215kwY7c/ikV59DStzcO1VEREREJD1SGBWP+PVX6NEDFi2CG26AL76Ali3Bzzg2dA6xI6CsW2cfcNdd8P77tlNRqVJ813++R7aDSS2pvXeqiIiIiEhG4+d2AZK2rVkDzZrZbLlhA4wcafPmow+exG/QANtYqEoVeO89KFQIRoyA6GhYuhS6dLH7tWC3gwkMSNqA6HLbwXhTSqHYV8KyiIiIiEhaozAq1+Tvv6FtW6hYEX7+Gd55xy4BfeF5h8wzp0OFCnaotFQpGD0adu+2J3bsCEWL/ud6zaoE0695RYKDAjFAcFAg/ZpX9JkpsL4elkVERERE0hpN05WrcuAAvPsujBplG9p26WLXiObLB6xdCy+/bPdvuflm+7XelTcfSm47GF9xJXunioiIiIjIlVMYlSty7BgMGQKDB8Px43ZUtHdvKF4c+PdfeLm3naObMyd88AE8//yZlrnphy+HZRERERGRtCZ9pQVJFX/+CY0awc6d8OCDdmS0fHkgIQE+/gR69bKBtEMH6NsX8ud3u2QREREREfFxCqNySfPm2Ya3uXLZnkN33XXmjkWL7JTclSuhVi07GnrrrW6WKiIiIiIiaYgaGEmKPvsM7r/f9iD67bczQXTHDrtnS+3acPAgfPUVLFigICoiIiIiIldFYVT+w3Fsd9w2bWzmXLwYiuWLhT59oFw5+O47eOstWL8eHnnEdjISERERERG5CpqmK0nEx8MLL8CYMfD44zB2LGT+Zwvcdx9s3gwtWsCgQVCypNulioiIiIhIGqYwKuccO2Zn4M6aBT172tFRs+ovCAuDuDj46Se45x63yxQRERERkXRAYVQA2LsXGjaEqCj46CN49lngl1/swezZbRCtUMHtMkVEREREJJ1QGBU2bLCNivbuhW+/tdu4MHu2baNbrBjMnatpuSIiIiIi4lFqYJTBnd2u5dgx2xS3USNsh9wmTSAkxHYvUhAVEREREREPUxjNwKZOhXr1IF8++PVXqFYNGD0aHn0U7rjDptNChdwuU0RERERE0iGF0Qxq6FDbGPe22+zS0DI3ONC/v10s+sADdppu7txulykiIiIiIumUwmgGk5gIXbrAK69As2a2L1H+fA68+iqEh0Pr1jB9OmTL5napIiIiIiKSjqmBUQbTrx/873/w4ov2q78TD+2ehXHjoGNH+OAD8Ev6GUVEVDSD5mxgV0wsRYMC6R4WQrMqwS49AxERERERSQ88NjJqjClijPnUGLPfGHPSGLPWGFPbU9eX6zd7NrzxBjz+OAwbBv7xp+zGouPG2TuGD082iIZPW0V0TCwOEB0TS/i0VURERbvzJEREREREJF3wSBg1xgQBSwEDNATKAy8C+zxxfbl+W7faGbiVKsHHH4M5fsy2zp02zQ6Rvv02GPOfxw2as4HYuIQkx2LjEhg0Z4O3ShcRERERkXTIU9N0XwV2O47z5AXHtnno2nKdTpyABx+0/542DbLFHoSGDSEyEj79FJ58MsXH7oqJvarjIiIiIiIiV8JT03SbAb8bY74yxuwzxqw0xnQyJpmhNvEqx4H27WHVKvjiC7gh90GoXRtWrrTJ9BJBFKBoUOBVHRcREREREbkSngqjNwAvAFuBMGAY0B/o6KHryzUaPtyG0L594b57E+GJJ2DTJvjhB2jS5LKP7x4WQmCAf5JjgQH+dA8LSa2SRUREREQkA/DUNF0/INJxnPAz30cZY27EhtERF59sjOkAdAAoUaKEh0qQiy1aBF27QtOmdtcW+vWzIfTDD6Fu3Su6xtmuueqmKyIiIiIinmQcx7n+ixjzNzDXcZx2Fxx7AvjIcZzsl3psaGioExkZed01SFLR0VC1KuTODcuWQe7l86FBA/4Ja0Krmp3YdfikgqWIiIiIiKQ6Y8xyx3FCLz7uqZHRpcDF8zZvAv720PXlKpw6BQ8/DMeOwfz5kPv4Lnj0UY6ULEOzik9y8PBJ4Pw2LYACqYiIiIiIeJWn1oz+D7jDGNPLGFPWGNMCeAkY6aHry1Xo3Bl++w0mTIAKN8XDo4/CsWM837QHB03mJOdqmxYREREREXGDR8Ko4zh/YDvqPgKsBt4F3gBGeeL6cuXGjYOPPoLXXrOjo7z+ul08+vHH/JKlULKP0TYtIiIiIiLibZ4aGcVxnJmO49zqOE5Wx3FuchznA8cTC1LlikVGwgsvQP368M47wIwZMGAAPPssPP64tmkRERERERGf4bEwKu7avx+aN4dChWDyZMi0c7vdQ/S222DoUEDbtIiIiIiIiO/wVAMjcVH8mWWh+/bB0qWQP+cpuL8FJCbCN99A1qyAtmkRERERERHfoTCaDvTqBT/9BOPH2+1c6NTVztmdPh1uuCHJuc2qBCt8ioiIiIiI6zRNN4375hsYOBCefx7atgW+/BJGjoSuXaFZM5erExERERERSZ7CaBoVERVNlS6/0fKxeHIUP0zdttGwfj20awc1akC/fm6XKCIiIiIikiKF0TQoIiqaHlNXse6rEEymRHI3/oN+U//gSKNmEBhoR0cDAtwuU0REREREJEVaM5oGDZqzgf3Li3JqVx7yPfAnmXKcpNes4eTYuhHmzIFixdwuUURERERE5JI0MpoG/bMrnpgF5chS7CDZb9lJy79+5KHV8/ngrlbQoIHb5YmIiIiIiFyWwmgadOrXW0g8nYm8967m5n1beXvuRywuWZmpDzzldmkiIiIiIiJXRGE0jVm6FPZFFiXv7dvJl2svoyL68W9gTno0f42u91dwuzwREREREZErojWjaUhcHDz3HJQoAQMGZyWw82cUO7yXFzsMofsTNbV/qIiIiIiIpBkKo2nIsGGwejVEREDTgGiInAWvvMKH77/sdmkiIiIiIiJXRWE0jfjnH+jdGxo1giaNEuHOjlCoELz1ltuliYiIiIiIXDWF0TSic2dITIThw8GMHwd//AETJ0KuXG6XJiIiIiIictUURtOAWbNg2jR47z0olesQ9OgBNWtC69ZulyYiIiIiInJNFEZ93IkT0KkTlC8PXbsCnV+HmBgYMQKMcbs8ERERERGRa6Iw6uPeew+2bYOff4bMq1fARx/Biy9CpUpulyYiIiIiInLNFEZ92Pr1MHAgPPEE1KmVCHd3ggIFoE8ft0sTERERERG5LgqjPspxoGNHyJ4dBg0CPvsMfv0Vxo+HoCC3yxMREREREbkuCqM+avJkmD8fRo2CQlli4NVX4c474ckn3S5NRERERETkuimM+qCYGOjSBapVgw4dgC5vwcGDMGcO+Pm5XZ6IiIiIiMh1Uxj1Qa+/Dvv32y1d/Nf8ZTvnPvccVKnidmkiIiIiIiIeoTDqYyIj7dTcTp3gtioO1OoIefJA375ulyYiIiIiIuIxCqM+JCEBnn8eChU6kz0nTYIlS2DMGMib1+3yREREREREPEZh1Id8/LEdGZ08GXKbI9C9u104+vTTbpcmIiIiIiLiUQqjPiAiKpp3p2xj+fu3k7vMUbLcFAt9hsDevfDdd2paJCIiIiIi6Y5SjssioqIJn7aKdd+VxIn3I/s9fzFu9AwShw2Ddu3syKiIiIiIiEg6ozDqskFzNnB4VyDHVxcjV9XtBOQ5Rs8fRnI0czZ47z23yxMREREREUkVCqMu2xUTy+ElN2Eyx5Prji00Wr+YO3esYmDNJyB/frfLExERERERSRUKoy7LdawgJzYWIVf1reT0P0Kv+WNZVagMi2o3c7s0ERERERGRVKMw6jL/FRXxDzxNrtDtvPTLlxQ5dpD37u9I1/sruF2aiIiIiIhIqlE3XRctXgxRv2alTefDnDK7efqPb5lR9T5avvQIzaoEu12eiIiIiIhIqlEYdYnjQM+eUKQIfPhebgIf/gZyZqfRrE+hYEG3yxMREREREUlVCqMumTMHliyBUaMgcOWvMGsW9O+vICoiIiIiIhmCwqgLHAdefx1KlYJnngEavgkFCkCnTm6XJiIiIiIi4hUKoy6YPh2WL4cJEyDzb4tg3jx4/33Int3t0kRERERERLzCOI7jagGhoaFOZGSkqzV4U0ICVKoEiYmwejX4168L69fDli2QLZvb5YmIiIiIiHiUMWa54zihFx/XyKiXTZ4Ma9fC11+D/8L5sGABfPCBgqiIiIiIiGQoGhn1org4KFcOcueGyD8c/GrXhO3bYfNmyJrV7fJEREREREQ8TiOjPmDcONi6FWbOBL+f5sLSpbadroKoiIiIiIhkMBoZ9ZKTJ6FsWShZEpYsdjB33Qm7d8PGjZAli9vliYiIiIiIpAqNjLrsww8hOhomTgTzwyz4/XcYM0ZBVEREREREMiSFUS84ehTeew8aNIA6tR0IfRNKl4Y2bdwuTURERERExBUKo14wbBgcOADvvAN89x2sWAHjx0NAgNuliYiIiIiIuEJhNJX9+y8MHgxNm0L10ESo8ibceCM8/rjbpYmIiIiIiLhGYTSVDRoER45A377AtGnw11924Wgm/acXEREREZGMS4koFe3da6foPvooVKyQAK3eshuNtmrldmkiIiIiIiKu8vPERYwxvY0xzkW3PZ64dlr23ntw6hT06QN8/TWsXQu9e4O/v9uliYiIiIiIuMqTI6MbgDoXfJ/gwWunOTt2wEcfwVNPQdlS8dCwN9xyC7Ro4XZpIiIiIiIirvNkGI13HCfDj4ae1bev/frGG8DkybBxI0ydCn4eGYwWERERERFJ0zyZjG4wxkQbY7YZY740xtzgwWunKZs22Z1bnnsOShSNt/N0K1eGZs3cLk1ERERERMQneGpk9HegLbAeKAi8DvxijLnZcZyDF59sjOkAdAAoUaKEh0rwHb17Q5Ys0LMn8NlnsGWL3V9Uo6IiIiIiIiIAGMdxPH9RY3IAW4H+juMMudS5oaGhTmRkpMdrcMv69VChArz2GvTrcxpCQiB/fli2DIxxuzwRERERERGvMsYsdxwn9OLjqbK1i+M4x4wxa4AbU+P6vmzAAMiaFbp0ASZMgO3bYdQoBVEREREREZELpMq8UWNMVqAcsDs1ru+r/v4bJk6EDh2gQK5T8M47cMcdcN99bpcmIiIiIiLiUzwyMmqMGQx8D+zArhl9A8gOfOqJ66cVgwbZAdBu3YBPPoF//oFx4zQqKiIiIiIichFPTdMtBkwG8gP7gd+AOxzH+dtD1/d5e/bY/NmmDRTLFwvvvQc1a0K9em6XJiIiIiIi4nM8EkYdx2nlieukZUOHQlwcvPoq8PHHsGsXTJqkUVEREREREZFkaK8RD/j3X9uj6JFH4MZisdC/P9xzD9Sp43ZpIiIiIiIiPilVuulmNCNGwNGjEB6O7aC7dy9Mnux2WSIiIiIiIj5LYfQ6HTtmp+g2bgyVKsTDg4Ph9ts1KioiIiIiInIJCqPXafRoOHQIevYEpkyBrVvh/fe1VlREREREROQSjOM4rhYQGhrqREZGulrDtTp1CkqXhvLl4ad5DlSpYg+uWQN+Wo4rIiIiIiJijFnuOE7oxcc1MnodJkyA3bvh88+BOXPgzz9h/HgFURERERERkctQGL1G8fEwYABUr24b51K3PxQrBq1bu12aiIiIiIiIz1MYvUZffQXbttnmReb332DhQhgyBDJndrs0ERERERERn6cweg0SE+G99+CWW6BRI+ChAZAnD7Rv73ZpIiIiIiIiaYLC6DX47jtYuxYmTQK/DesgIgLefBNy5CAiKppBczawKyaWokGBdA8LoVmVYLdLFhERERER8SkKo1fJceyoaJky8MgjQPuBEBgIL75IRFQ04dNWERuXAEB0TCzh01YBKJCKiIiIiIhcQG1fr9K8efDHH/Daa5Bp9z8wcaKdnps/P4PmbDgXRM+KjUtg0JwNLlUrIiIiIiLimxRGr9J770FwMDz5JLZhEUCXLgDsiolN9jEpHRcREREREcmoFEavwi+/wIIF0K0bZDl2EEaPhkcfhZIlASgaFJjs41I6LiIiIiIiklEpjF6F996DfPnONM0dMQJOnIBXXz13f/ewEAID/JM8JjDAn+5hIV6uVERERERExLepgdEVWrkSZs6Evn0hO8dh+HBo3Nju73LG2SZF6qYrIiIiIiJyaQqjV6h/f8iZEzp2BMaOhYMHoUeP/5zXrEqwwqeIiIiIiMhlaJruFdi4Eb7+2gbRPDniYPBgqFkT7rrL7dJERERERETSJI2MXoEBAyBLFujcGZg8Gf75Bz76yO2yRERERERE0iyNjF7Gjh3w2WfQrh0UKpBok2nFinD//W6XJiIiIiIikmZpZPQyBg+2X7t3B2bMgLVrYeJEMMbVukRERERERNIyhdFL2LcPxoyBJ56AEsUdaNUfSpWCli3dLk1ERERERCRNUxi9hH/+gRtugNdeA5YsgV9/tfuLZtJ/NhERERERkeuhVHUJVavC6tVnZuR26Q8FCsBTT7ldloiIiIiISJqnBkaXYQzw118waxa8/DJky+Z2SSIiIiIiImmewuiVGDAAcuSAF15wuxIREREREZF0QWH0EiKiomnx6kQSJn/JF1XuJ2L7CbdLEhERERERSRe0ZjQFEVHRhE9bRfi8L0kwfgyr2JAj01YB0KxKsMvViYiIiIiIpG0aGU3BoDkbyBZzkEdWzWXaLfewN2d+YuMSGDRng9uliYiIiIiIpHkKoynYFRNLkaMH2JG7MKOrN09yXERERERERK6PpummoGhQIKspy73PjDzTUvf8cREREREREbk+GhlNQfewEAID/JME0cAAf7qHhbhYlYiIiIiISPqgkdEUnG1SNGjOBnbFxFI0KJDuYSFqXiQiIiIiIuIBCqOX0KxKsMKniIiIiIhIKtA0XREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOoVRERERERER8TrjOI67BRizH/jb1SIuLz9wwO0iJMPT61B8gV6H4iv0WhRfoNeh+IK08Dos6ThOgYsPuh5G0wJjTKTjOKFu1yEZm16H4gv0OhRfodei+AK9DsUXpOXXoabpioiIiIiIiNcpjIqIiIiIiIjXKYxemdFuFyCCXofiG/Q6FF+h16L4Ar0OxRek2deh1oyKiIiIiIiI12lkVERERERERLxOYVRERERERES8TmH0EowxLxhjthljThpjlhtjarpdk2Qcxpjexhjnotset+uS9M8YU8sY850xJvrM667tRfebM6/PXcaYWGPMAmPMzS6VK+nUFbwOJyTzHvmbS+VKOmWMCTfG/GGMOWKM2W+M+d4Yc8tF5+g9UVLVFb4O0+R7osJoCowxLYFhwHtAFeAX4AdjTAlXC5OMZgNQ5IJbRXfLkQwiB7AaeBmITeb+V4GuwItANWAfMNcYk9NrFUpGcLnXIcA8kr5HPuCd0iQDqQOMAu4C7gHigXnGmLwXnKP3REltdbj86xDS4HuiGhilwBjzO/CX4zjtLzi2CZjiOE64e5VJRmGM6Q087DjOLZc7VyS1GGOOAZ0cx5lw5nsD7AJGOI7z7pljgdg/vro5jvOxW7VK+nXx6/DMsQlAfsdxGrlVl2Q8xpgcwGGgmeM43+s9Udxw8evwzLEJpMH3RI2MJsMYkxmoCvx40V0/Yj+REPGWG85MUdtmjPnSGHOD2wVJhlcaKMwF74+O48QCi9D7o3jf3caYfcaYjcaYMcaYgm4XJOleTuzfz/+e+V7vieKGi1+HZ6W590SF0eTlB/yBvRcd34t9wxHxht+BtsD9QHvsa+8XY0w+N4uSDO/se6DeH8Vts4EngXrYKZLVgfnGmCyuViXp3TBgJfDrme/1nihuuPh1CGn0PTGT2wX4uIvnMJtkjomkCsdxfrjw+zOL0LcCbYAhrhQlcp7eH8VVjuN8ecG3q4wxy4G/gYbANHeqkvTMGDMEuBu423GchIvu1nuieEVKr8O0+p6okdHkHQAS+O8nWgX57ydfIl7hOM4xYA1wo9u1SIZ2tqOz3h/FpziOswvYid4jJRUYY/4HPArc4zjO1gvu0nuieM0lXof/kVbeExVGk+E4zmlgOdDgorsaYLvqinidMSYrUA7Y7XYtkqFtw/7xde798cxrsyZ6fxQXGWPyA8HoPVI8zBgzDGiNDQDrL7pb74niFZd5HSZ3fpp4T9Q03ZQNAT43xiwDlgLPAUWBj1ytSjIMY8xg4HtgB/YT1jeA7MCnbtYl6d+ZLn1lz3zrB5QwxlQGDjmOs8MYMxToZYxZD2wEXgeOAV+4UK6kU5d6HZ659QamYv/QKgX0w3Ywne7lUiUdM8aMBJ4AmgH/GmPOjoAecxznmOM4jt4TJbVd7nV45v2yN2nwPVFbu1yCMeYF7N5RRbB7nb3iOM4id6uSjMIY8yVQC9tQaz/wG/CG4zhrXS1M0j1jTB3g52Tu+tRxnLZntjJ4C3gWyINtttXRcZzVXitS0r1LvQ6B54EI7D7gQdg/vn7Gvkf+45UCJUMwxqT0h3Ifx3F6nzlH74mSqi73OjyznVAEafA9UWFUREREREREvE5rRkVERERERMTrFEZFRERERETE6xRGRURERERExOsURkVERERERMTrFEZFRERERETE6xRGRURERERExOsURkVERERERMTrFEZFRERERETE6xRGRURERERExOv+D8c86GemK2bWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x576 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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:34.763089Z",
     "iopub.status.busy": "2021-02-02T06:51:34.761717Z",
     "iopub.status.idle": "2021-02-02T06:51:34.774030Z",
     "shell.execute_reply": "2021-02-02T06:51:34.775034Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:34.779526Z",
     "iopub.status.busy": "2021-02-02T06:51:34.778179Z",
     "iopub.status.idle": "2021-02-02T06:51:34.788582Z",
     "shell.execute_reply": "2021-02-02T06:51:34.789551Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept           4.907719\n",
       "x1                  0.517908\n",
       "np.sin(x1)          0.495210\n",
       "I((x1 - 5) ** 2)   -0.021164\n",
       "dtype: float64"
      ]
     },
     "execution_count": 1,
     "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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T06:51:34.793911Z",
     "iopub.status.busy": "2021-02-02T06:51:34.792574Z",
     "iopub.status.idle": "2021-02-02T06:51:34.803099Z",
     "shell.execute_reply": "2021-02-02T06:51:34.804074Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    10.933725\n",
       "1    10.780020\n",
       "2    10.514292\n",
       "3    10.180798\n",
       "4     9.837796\n",
       "5     9.543278\n",
       "6     9.340778\n",
       "7     9.248708\n",
       "8     9.255852\n",
       "9     9.324113\n",
       "dtype: float64"
      ]
     },
     "execution_count": 1,
     "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}