{
 "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 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": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "sig = 0.25\n",
    "x1 = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x1, np.sin(x1), (x1-5)**2))\n",
    "X = sm.add_constant(X)\n",
    "beta = [5., 0.5, 0.5, -0.02]\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + sig * np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.980\n",
      "Model:                            OLS   Adj. R-squared:                  0.978\n",
      "Method:                 Least Squares   F-statistic:                     741.1\n",
      "Date:                Sat, 07 Nov 2020   Prob (F-statistic):           6.23e-39\n",
      "Time:                        11:34:25   Log-Likelihood:                -3.0227\n",
      "No. Observations:                  50   AIC:                             14.05\n",
      "Df Residuals:                      46   BIC:                             21.69\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.2072      0.091     57.006      0.000       5.023       5.391\n",
      "x1             0.4696      0.014     33.337      0.000       0.441       0.498\n",
      "x2             0.4563      0.055      8.239      0.000       0.345       0.568\n",
      "x3            -0.0176      0.001    -14.223      0.000      -0.020      -0.015\n",
      "==============================================================================\n",
      "Omnibus:                       10.776   Durbin-Watson:                   2.145\n",
      "Prob(Omnibus):                  0.005   Jarque-Bera (JB):               13.024\n",
      "Skew:                          -0.755   Prob(JB):                      0.00149\n",
      "Kurtosis:                       4.992   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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.76736342  5.20902136  5.61506369  5.9606248   6.22981294  6.4183212\n",
      "  6.53413514  6.59622074  6.6314084   6.66998472  6.74071656  6.86612468\n",
      "  7.05878367  7.31925595  7.6359994   7.98726368  8.34466375  8.67784345\n",
      "  8.95946295  9.16969058  9.29946112  9.35196523  9.34212516  9.29414314\n",
      "  9.23752499  9.20223238  9.21376024  9.28894802  9.43321299  9.6396602\n",
      "  9.89021532 10.1585942  10.41462126 10.62918728 10.77903194 10.85056542\n",
      " 10.8421012  10.76413314 10.63761134 10.49050013 10.35318359 10.25347355\n",
      " 10.21204017 10.23901591 10.33233075 10.47805236 10.65267507 10.82698183\n",
      " 10.97084449 11.05817373]\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.06298811 10.94765308 10.73006132 10.45098799 10.16410754  9.92285236\n",
      "  9.76733063  9.71450645  9.75404645  9.85084976]\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": "iVBORw0KGgoAAAANSUhEUgAAA6MAAAHWCAYAAACCFl9HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAB3QUlEQVR4nO3de3yO9R/H8de1mRkbc2ZzPi0ipwmJHKrpQEhHHYR0oKgs9tPBoUJDKSlUEhWFVk7NIYccwrScWzkVQ4hhDLNdvz++jM0w3LuvHd7Px+N+zK7rvq/7s9zd9r6/3+/na9m2jYiIiIiIiIg7eThdgIiIiIiIiOQ+CqMiIiIiIiLidgqjIiIiIiIi4nYKoyIiIiIiIuJ2CqMiIiIiIiLidgqjIiIiIiIi4nZ5nC6gWLFidoUKFZwuQ0RERERERDLB2rVrD9q2XTztccfDaIUKFYiKinK6DBEREREREckElmX9nd5xTdMVERERERERt1MYFREREREREbdTGBURERERERG3UxgVERERERERt1MYFREREREREbdzvJvulRw9epT9+/eTmJjodCmSg3h5eVGiRAkKFizodCkiIiIiIrlSlg6jR48e5d9//yUwMBAfHx8sy3K6JMkBbNsmISGB2NhYAAVSEREREREHZOlpuvv37ycwMJD8+fMriIrLWJZF/vz5CQwMZP/+/U6XIyIiIiKSK2XpMJqYmIiPj4/TZUgO5ePjo+nfIiIiIiIOydJhFNCIqGQavbZERERERJyT5cOoiIiIiIiI5DwKoyIiIiIiIuJ2CqOZoHPnzliWhWVZKVuItGjRgo8++uiq1iguXrwYy7I4ePBgJlYrIiIiIiLifgqjmeT2229n79697Ny5k3nz5tGmTRvefPNNmjZtyvHjx50uT0RERERExFG5IoxGRMfSZOjPVOw3myZDfyYiOjbTn9Pb25tSpUoRGBhInTp1ePnll1m8eDG//fYb7777LgCTJ0+mQYMG+Pn5UaJECR544IGUvS937txJixYtAChevDiWZdG5c2cAfvrpJ5o2bUrhwoUpUqQIISEhbNmyJdN/JhEREREREVfJ8WE0IjqWsBkbiI1LwAZi4xIIm7HBLYE0rZo1a9K6dWumT58OwOnTpxk4cCDr1q1j1qxZHDx4kEceeQSAsmXLptxv06ZN7N27l1GjRgFw/PhxevfuzerVq1m8eDGFChWiTZs2nD592u0/k4iIyIWc+ABYRESypzxOF5DZwiNjSEhMSnUsITGJ8MgY2tUNdHs9NWrUYMGCBQB06dIl5XilSpX4+OOPqV69Ort376ZMmTIUKVIEgBIlSlCsWLGU+95///2prjlhwgQKFizI6tWrufXWW93wU4iIiFzs3AfA5/7dPfcBMODIv7kiIpK15fiR0T1xCVd1PLPZtp2yv+Vvv/3GfffdR/ny5fHz8yM4OBiAf/7557LX2LZtG48++iiVK1emYMGClCxZkuTk5Cs+TkREJDNd7gNgERGRtHJ8GA3w97mq45lt8+bNVKpUiePHjxMSEkL+/PmZNGkSa9as4aeffgK44nTbNm3acODAAcaOHcuqVauIjo4mT548mqYrIiKOymofAIuISNaW48NoaEgQPl6eqY75eHkSGhLk9lo2btzITz/9RMeOHfnjjz84ePAg77zzDs2aNeOGG25g//79qe6fN29eAJKSzn/K/N9//7Flyxb+97//cfvtt1O9enWOHTvGmTNn3PqziIiIpJXVPgAWySxaGy3iGjk+jLarG8iQDrUI9PfBAgL9fRjSoVamr105deoU+/btY8+ePaxbt46RI0fSvHlz6tevT58+fShXrhze3t6MHj2a7du3M3v2bF5//fVU1yhfvjyWZTF79mwOHDhAfHw8hQsXplixYowfP56tW7eyZMkSnn32WfLkyfHLf0VEJIvLSh8Ai2SWrNQcUyS7y/FhFEwgXd6vJTuG3sPyfi3d0kRhwYIFlC5dmnLlytGqVSt+/PFH3nzzTZYuXUqBAgUoXrw4EydOJCIigho1ajBw4EBGjhyZ6hqBgYEMHDiQ/v37U7JkSXr27ImHhwdTp05l/fr11KxZkx49ejB48GC8vb0z/WcSERG5HKc+ABZxJ62NFnEdy7ZtRwsIDg62o6Ki0j23ZcsWqlev7uaKJDfRa0xERESuRsV+s0nvt2cL2DH0HvONbUNyMiQmwpkzl/5asiT4+7uxehFnWJa11rbt4LTHMzS307KsZkAfoD4QADxl2/YXF5zvADwD1AOKAS1s2158/WWLiIiIiGQdVfMl0WBlJPdvXEjpowfJYyeRJykJLzsZPkg+HzYzwtsbHn4YevaE4It+TxfJ8TK60NAX2Ah8efaWVgFgBTD5EudFRERERLIn24bly+HTT5kzdSp5Tp5kS/EKLK1YjyQPT2wvLxpVLYFvYGHIkwe8vK781dMTVqyAiRPNrWFDE0ofeMCEVJFcIENh1LbtOcAcAMuyvkjn/KSz54q5sjgREREREcfs3w9ffgmffgoxMeDnR54nn2Rx07b0/8ebPUdOEuDvQ2hIEJWvZW10p07wzjvmOUaPhscfh5dfhu7d4ZlnoGxZ1/9MIlmIWrCKiIiIiJyTlAQLFpgA+sMPZsrtLbfA55/Dgw9CgQI0B5a76vkKFYIXXoAePWDhQhNK33kHhg6Fdu3MaOltt4FlueoZRbIMR7rpWpbV3bKsKMuyog4cOOBECSIiIiIi5/3zDwwcCJUqQevWsGiRCYmbNpkpuk89BQUKZN7ze3jAHXeYALxtG7zyiqmhRQuoVQs+/hji4zPv+UUc4EgYtW17nG3bwbZtBxcvXtyJEkRERCQLiYiOpcnQn6nYbzZNhv6sPRvFfbZsgXvugQoVYMAACAqCb7+F2FgYMQJq1LjkQ5OTM6mmihVh2DDYvduMyHp7w/PPQ2Ag9OplpgyL5AC5Yp9RERERyboiomMJm7GB2LgEbCA2LoGwGRsUSCVznT4NgwZBnTrw66/w2muwfTvMm3fZJkKHDsEnn0CTJmY2baby8TEjslFRsHIltGljRkhvuAFCQuCnnzK5AJHMpTAqIiIijgqPjCEhMSnVsYTEJMIjNfojmWTVKqhfH958Ezp0MKOjgwaZEcl0nDoF339v7lq6NDz3HMTFQUCAm+q1LGjUCCZPhl27YPBg2LgR7rrLTCXO6FYyIllMhsKoZVm+lmXVsSyrztnHlDv7fbmz54ucPVfz7EOqnD1fKjOKFhERkazryBEziPPXX3DsmNkV43L2xCVc1XGRa3b8OLz0EjRubNLkzJnwzTdQosRFdz23m8uzz5oA2qGD2YmlRw9Yu9ZkwS5d3P8jULKkGcXduRP69DENj1q1gn37HChG5PpktJtuMLDogu8Hnr1NBDoDbYEJF5wff8H9BlxXhSIiIpLlbd0Ks2aZ3+2XLoUzZ86fy5/f/P5cqtTFt5Ilwe9oSQ7ZR/EscAorz/lFeAH+Pg78JOJqEdGxhEfGsCcuIWUblHbXsg3K9Zo3z2yXsnOnGdocOhQKFrzobn/9ZQYgJ082s3Z9fKB9e7Pryu23m21CswQvLwgPNyO8XbqYr9OnmxFUkWwio/uMLgYu2U/atu0vgC9cUlE2Z12h7faTTz7JF1984Z5iREREMsmZM2b0c+ZMc/vjD3P8xhvNYM0tt5gR0n37Ut/+/NOE1f/+u/BqweaLZeNT9V8K1t+Bf6UjhIYEufvHEhc7tx743DTsc+uBAfcF0kOHzN6dEydCtWrmBdi0aaq7HDxoehZNmmSWj1oWtGwJb7xhRkT9/NxT6jV5+GGoXt0k5ttuMyOlTz/tdFUiGZJVPtvJMfbu3Zvy51mzZvH000+nOubjk/pT3sTERLy8vNxWn4iIyLU6cgQiI034nDPH/I7v5WV+/33uOdNb5RJL7i5y+jQcOHA+pM5dfZjvFx1jX1Qp/v2zFPmDTnOkal5O1bhkHxnJBi63HjjTw6htw3ffmTWVhw5B//5memu+fCl3SUw0g4uDBpl1oTVrmia2jz4KZcpkbnkuVbu2aXL0yCPQvbv58wcf6H8eyfIURl2sVKnzy2T9/f1THdu5cyelS5fm66+/Zvz48axcuZLw8HB8fX3p2bMn8RfsHbV48WJatGjBgQMHKFasGAArVqwgLCyMNWvWULhwYdq2bcuwYcMomM4UExEREVfYtw+mTDFTcJcsMSOiRYuanTDatDENPVP9MxQfD3Pnwpo15hdhHx8zTzd//lR/zps/P4E+PgTmzw/V83NP/fyMDi1CgmdevvoK3n8/L507w6uvmqD77LNmWq9kL46tB96922yFMnMmBAebKbq1a6e6y+rV0K0bbNgAHTuarFq7thkVzZaKFDGfEr3+OgwZAuvXw7RpZjsYkSxKYdQBYWFhDB8+nM8++wwvLy8WLFhwxcds2LCBO++8k4EDB/Lpp59y6NAhevfuTZcuXZg2bZobqhYRkdzEtuHrr02zliNHzFaLr7xiAmijRuDpecGdDx0yv/TPmGGGTk+dgrx5TXK9yo0YfRo1otsjj9B1/oMs3FSKUaNg4ECzhcYjj5gtFuvVc+3PKpknwN+H2HSCZ6atB05OhvHjzacYiYkwfLh50Vyw0PP4cZPXRo0yH3BERMB992VOOW7n6Wn+Z6lXDzp3NutIp02DW291ujKRdGW7MNq7N/z+u3ufs04deP99113vhRdeoGPHjlf1mPDwcB566CFeeeWVlGMff/wxdevWZf/+/ZRIpwuciIjItThwwIxGTp9u1n6OG2fWgqayd6/5LX7GDFi0CJKSoGxZM4TZoYPZhNHDw8zHTUiAEyfO3y78/sI///uv2T+jVy+sl17i9pYtuf2RR9g2oAOjJvozYQJ8+aX5vbp3b6BcLCMXZoHGOHJJoSFBqdaMAvh4eWbOeuC4ONNlaNYss+Bz3DioXDnVXSIjzUt0507zdehQKFTI9aU4rmPH8+tIW7Qwv8g+/3w2HvaVnCrbhdGcIDg4+Kofs3btWrZu3crUqVNTjtlne+Vv27ZNYVRERFzixx9N75PDh80v6n36XDAKumOHCZ8zZpjuRbZtGsK8+qoJoPXrX/zLrre3uZ1dunJFb74Jmzeb7Ta+/hq6dqVy3uf44O67GTL6UT7bdy/vj/WhY0fIU7AIvvVK41fnH2ca48gVnfu7yPRuuuvWwf33w99/m7WSPXumei3+95/Z0WXSJAgKSreHUc5z441mLvJjj5n/HmvWwMcfm+nyIllEtgujrhyhdEqBAgVSfe/h4ZESLM9JTLN5cXJyMt26deOll1666HqBWgsgIiLX6ehR88v655/DTTfB/PnmK//8Y4Yjp08/PzWpbl3T8aVDBzP64urRlho1YPBg8xxr1phgOmUKBSIieNHXl57t2vNMhUZ8s7UDcYurc2xtBYretR4qHnRPYxy5Ku3qBmbu38mkSWbLlsKFzcLmW25JOWXb5uXTu7f5gOW118za0At6GOVs/v7mE6aBA83/Txs3mg+TypVzujIRIBuG0ZyoePHinDhxgqNHj6Y0I/o9zVzkevXqsWnTJqpUqeJAhSIikpMtXmyWl+3aBWFhZnDS2ysZPh4LoaFmCu0tt8CIEWbaX0Zb5l4vy4Kbbza34cNNod98g8f06YyPm8RQnzBmVmvJ+/+Gse7bhvjW+Ru7+Rb31CbOO33afIIyZoxp6Tx1qtm49qx//jHTzefMMS+hhQuhVi0H63WKh4cJo/Xrm1HS+vXNPjYtWjhdmQgeThcg0LBhQwoUKEBYWBhbt25l+vTpjBkzJtV9+vbty+rVq3n22WeJjo5m69atzJo1i2eeecahqkVEJDuJiI6lydCfqdhvNk2G/kxEdCwJCeZ3+RYtzBYty5aZ3ife+/6GO+80a8yaNIHt283Jl192XxBNy9MTWrWCTz+Fffvo+/ggVpSvzUPb57L2aGM+K/4oPr/n49+Jt7FokTMlihvt3m0C6JgxZi75ggUpQTQpCT780AywL14M770HK1bk0iB6obZtzbTdYsXgjjtg5EgzdCziIIXRLKBIkSJ89dVXzJ8/n1q1ajFu3DgGDx6c6j433XQTS5cuZefOndx2223Url2bsLAwSl7wCaCIiEh6IqJjCZuxgdi4BGwgNi6B3qP/pmqNxJS+Jr//Do0b2abpS82asGqV+fNPP0GFCs7+AGl5e9P4pS682vF/NH5+Al/Wu4cnDk5lq1dlXkseReuWp3jxRdM1VXKgRYtMt9iNG80+ouHhKd1y//rLNLh68UXzddMmM0U3Vffn3OyGG8z/223amPbYnTqZJmIiDrHSrlV0t+DgYDsqKirdc1u2bKF69epurkhyE73GRCQ3aDL055TtNewkiyMrq3BkRRXy+p1m5nf5uPNOzBzdrl3NYtFWreCzz6B8eWcLv4KI6NiUxjiNTh/gvTWTKbVsIQcLVuKZo++yvnIHvpho0aSJ655HXXsdZNsmeIaFmS5EM2aYcHXWrFkmW+XJY7Zt6dQpY8uZc+Xfb3Ky2Yv09dfNZsHff5+LFtKKEyzLWmvb9kVdXLVmVEREJIfbczaIJh705eDs2pze50+BG3dT5PZN3HnHnfDZ52YKblKS6bb5zDPZYguIixvjdIb58yn28stM39iRVbub0ePWkTR/pT6DB19bE9Fzo8rntiZR116HHD0KTz1lAugDD5gPS/z8AJOrBg0yyyLr1jV3yehgfq79+/XwMJ2cSpY07bMfeMA0Kcub1+nKJJfRNF0REZEcLsDfh5O7C7N3UhPOHMlPsXZrKXbvOur5xMHdd0O3bmba44YNZvPFs0E0vXWmWd4dd0B0NHzyCQ38trCaBtQc0ZmQWntYvfrqLxceGZNqj0yAhMQkwiNjXFSwpOfC195jL33OsZvqwg8/mCZaU6emBNG4OLjvPhNEn3gCli+/ulnluf7vt1s3s+521ix4+GFIs5uDSGZTGBUREcnh2pSqxYHvbsazwElKP7WUAtX28vDmn/ludHez4eKHH5pWoxc0J0pvnWnYjA3ZI5DmyQPPPIPH1r/weDWUJ7y+4aftVZndaBADXj3BqVMZv9S5UeWMHpfrd+Fr754tSxn7UQ9OHTzEL59MMSP4Zz8s2bgRGjQwy5pHj4Yvvrj60W/9/WJaDn/wgZmq26kTnDnjdEWSiyiMioiI5GArV8I7LxanZCmo8+zvBLKHyT++zdCZI/GqWxvWr4eePc20vQvkiBGjQoVg2DA8/tiCV9u7GWi/SdfwIAZUmUz02uQMXSLAP/10c6njcv3CI2NIPHmK1xeOZ/SP77KlREXuefJ9+h0sknKfb7+FRo0gPt50zO3R49pmluvv96wXXjCjzt99Z4aYk5Ku/BgRF1AYFRERyaFWrjS9SUqWhDUr8rDmpt2s+qY3t/69Dt5/3/wWX7lyuo/NUSNGlSrhFfEdLF2KX5WSDNn9OKcaNOG7gZuvuLNFaEgQPl6pW7H6eHkSGhKUiQXnbqd37+GrKf3pGvUDE+q34ZFH3uFfv2LsiUvgzBmz9e1DD0Ht2rB2LdfVoEp/vxd4+WUYOhS++casz1UgFTdQAyMREZEc6NdfTRAtUcLshBH47Xvml83Gjc18xmrVLvv4AH+flA68aY9nW02b4h+zmmNjJlH95VeoO6Au02cO4J5Fofj4pf8r0bkmNrmu26pTli1j7pe9yX8ynhfb9OHHGs1TThX3Kkjr1mZG+XPPmc9TLuy3cy1dcfX3m0bfvmbd6Ouvm+nun3560awJEVdSGBUREclhVq06H0QXL4YyE9+G114zHTMnT85Qx8zQkKBUXUYhh4wYeXjg1/NJkjq0ZsvtPem49n9sLjmdgtM+p8zdN6X7kIu79orL2bZZu/zKK+QrXYaH7xrM+sLlUk5bBwuzc/bNHD0Mn39uBu4udD1dcfX3m8Zrr5lAOmgQeHmZDtsKpJJJ9MoSERHJQVatgjvvhGLFYNHPNmU+7m9+uXz8cfj66wxv3dCubiBDOtQi0N8HCwj092FIh1o55pd2z4CS1Nz8HWv7fUfxk/9Q4p5gYjoNhNOnnS4t9zl+HB57DHr1grvuwnd9NF2ebZvy2vPaVpHYSY3J55WHZcsuDqKQQ9Y4ZyUDBpj9XMeNM+tJrzSfXeQaaWRUREQkh1i9+nwQXbzIpux7L5u5jN27X9PoRk4cMbpoKueDjSnWcTNr7ujF3V8PYN+CGRSfNQHPBvWcLjV3+Osv6NABNm2Ct9+Gfv3Aw4N2df25+8ZAXnoJxkyDli1hyhQoXjz9y+SoNc5ZgWWZv4/ERBg+3IyQvvdetth/WLIXjYzKdevZsyfNmzdP+b5z587ce++913XNAQMGULNmzeusTEQk91izxgTRokVh0cJkyr5zdlFdr17wySeaZselt6uJ9jhFi9iv+OD2H0jefwAa3syJl/pzVXvAyNX74QcIDoa9e83+LP/7X8rrdO9eaNHCbIHZpw9ERl46iIK64mYKy4J33zXvIaNGwauvaoRUXE7/MmWS2NhYunfvTpkyZcibNy+BgYE8/fTT7N69O9X9rhTc1q1bx3333UepUqXIly8f5cqV4/777+fvv//O7B/hmo0aNYrJkydn6L47d+7EsiyioqJSHe/Tpw9LlizJjPJERHKcqCi44w4oUgQWLzhDuTc6w9ixZpqdRjNSXG4qp48PvDi/LQve38RX1uPkf/8dEqrXM/OexbWSkqB/f2jXzjTSWrvWfJJy1ooVUL8+/P67GQ0NDze9dC5HXXEziWWZ95AePcwIaf/+CqTiUgqjmWDHjh0EBwezceNGJk6cyNatW5k8eTKbNm2iQYMG7Ny5M0PXOXDgAK1atcLX15fZs2fzxx9/MGnSJCpXrszRo0ddWvNpF66RKVSoEP7+/td1DV9fX4oWLeqagkREcrC1a00QLVwYFkWeply/R2HSJHjrLXjnHQXRC2RkKucTvQpz4+oJPFliLgd3HCW58S3YfUIhQdM9wYwuNxn6MxX7zabJ0J+JiI69ugscPAh33WVem926wS+/QPnygMk4n3wCzZtD/vymI/RDD2Xssjl9jbOjLAs++MBM9x8yxKwnFXERhdFM0KNHDzw8PFiwYAGtWrWiXLlytGjRggULFuDh4UGPHj0ydJ3ly5dz+PBhJkyYQP369alQoQK33XYb7777LrVq1brk486Ntr711luULFkSX19fnnrqKRIu+Ie0efPmPPfcc/Tp04fixYvT5OwmXZs3b+aee+7Bz8+PEiVK8Mgjj7Bv376UxyUlJdGnTx8KFy5M4cKF6d27N0lp9qFKO9pr2zYjRoygatWqeHt7U6ZMGcLCwgCoWLEiAA0aNMCyrJTpvmmn6SYnJzN48GDKli2Lt7c3tWrV4ocffkg5f26Edfr06dxxxx3kz5+fGjVqMH/+/Az9txYRyY7WroXbbwd/f1j800nKv9LRbFo/YoQZwZBUMjqVs359GLm5Nb1v38R4uxvWiOEk164Dy5a5ocqs61LTnDMcSNesMf9xly6F8ePNLV8+AE6eNNn0uefMa3rNGrjMrzrpalc3kOX9WrJj6D0s79dSQdSVPDzMuvMuXUyX3bfecroiySEURl3s0KFD/PTTT/To0YP8+fOnOpc/f36ef/555s6dy+HDh694rVKlSpGcnMy0adOwr3JKxJIlS1i3bh0LFy5k+vTpzJs3j759+6a6z+TJk7Ftm19++YUvv/ySvXv30qxZM2rWrMnq1atZsGAB8fHxtG3bluTkZABGjBjB+PHjGTt2LCtXriQpKYmvvvrqsrX873//Y/DgwYSFhbFp0ya+++47ypYtC8Dq1asB+Omnn9i7dy8zZsxI9xqjRo0iPDycYcOGsWHDBtq3b0+HDh34/fffU92vf//+vPjii6xbt44GDRrw8MMPEx8ff1X/7UREsoPffjMjooUKweI5Jyj/QluYOdMssnv5ZafLy5KuZipn0aLwXWRB9g0Yy+0sYO/fp6FpU3jmGTh0yF0lZynX3LHWts208VtvNd8vW2aS51m7dkGzZmbLltdeMy/jwoVdXb1cNw8P01338cfNPqTDhjldkeQA2a+bbu/eZhGBO9WpY5pAZMBff/2FbdtUr1493fM1atTAtm3++usvbr755steq1GjRvzvf//jySefpEePHjRo0IDmzZvTqVMnyp+d0nIpnp6eTJgwAV9fX2rWrMmwYcPo2rUrQ4YMoUCBAoAZlRwxYkTKY9544w1q167NsAveXL788kuKFClCVFQUN998M++//z6vvvoqDz74IGBCYmRk5CXriI+P57333uP999+nS5cuAFSpUoXGjRsDUPxsN4KiRYtSqlSpS15n+PDh9OnTh0cffRSAQYMGsXTpUoYPH55qfepLL71EmzZtAHjnnXf48ssv+f3337n13D+AIiI5wJYtZvSoYEFYMusY5Z+5B5Yvhy++gCefdLq8LOvcSFmqbrohQZccQfPwgDffhJ8atqLRIxt4lQH0+PR9PL7/3qyfe/zxXDUN+po61u7caaZ3zp9v1oV+9ZVp93zWkiVm+9uTJ+H7780yUsnCPD1hwgQ4c8Z0Pvby0odfcl00MppJrEv843RuhPNS59N6++232bdvH+PGjaNWrVp89tln1KhRg4ULF172cTfddBO+vr4p3zdu3JjTp0+zbdu2lGP169dP9Zi1a9eydOlSfH19U27nRjC3bdvGkSNH2Lt3b0qQBPDw8KBhw4aXrGPz5s2cOnWKVq1aZejnTc/Ro0fZs2dPylTic2699VY2b96c6thNN53fsDwgIACA/fv3X/Nzi4hkNYcOQZs25nfAJRGHKd/tDli50uwhqiB6RdcylbN1a/gl2pev6g6nTvJvbPOoYv5bt2gBaf4dysmuqmNtcjKMHg01a5rX50cfwdy5KUHUtk2D1latTOOt1asVRLMNT0/48kvzKcIrr5i/Z5FrlP1GRjM4QumUqlWrYlkWmzZtol0676pbtmzBsiwqV66c4WsWLVqUBx54gAceeIAhQ4ZQt25dBg8efF0BD0gZIT0nOTmZe+65h+HDh19035IlS6ZM1b0aVzu9+HLSC/Bpj3l5eV107lrqFhHJihIT4cEHzbTGX2YcoPxTd5owNG0a3Hef0+XlaBUqmF47b7xxE9WGLuO10p/xenRf8tSuDaGhZn5pmuU5OU1oSBBhMzakmqqb7jTnmBgzDXfZMggJMVN0L5jRdeKEGSz96isTQCdONKP8ko3kyWP+AhMT4YUXoEABeOopp6uSbEgjoy5WpEgRQkJCGDNmDCdOnEh17sSJE3z00UfcddddFClS5JqunzdvXipXrnzFdZAbNmzg+PHjKd//+uuvKY+9lHr16rFp0ybKly9PlSpVUt38/PwoVKgQpUuX5tdff015jG3bKes+01OjRg28vb0vOZKbN29egIuaIF2oYMGCBAQEsCxN44hly5ZRo0aNSz5ORCSneeklWLgQJg4/wM2vNoc//oAff1QQdRMvL9NM9Kd5Hoyzn6biqRj+CO5kDt54I8ye7XSJmeqKHWvPnDHrCGvXhk2bzLTxuXNTBdGdO6FJEzOQP3gwTJ+uIJpteXmZvXfuvNN8+DB1qtMVSTaU/UZGs4HRo0dzyy23cPvtt/PWW29RtWpVtm3bRv/+/bFtm9FppjMcPXr0okY8/v7+bNy4kSlTpvDwww9TrVo1bNtm5syZzJkzh4EDB162hjNnztClSxfeeOMN9uzZQ79+/Xj66acvGg29UI8ePRg/fjwPPfQQffv2pXjx4mzfvp1vv/2WESNG4OfnR69evRgyZAjVqlWjVq1ajBkzhr1791K6dOl0r3nuMWFhYXh7e9OsWTP+++8/1q5dy3PPPUeJEiXw8fEhMjKSChUqkC9fPgoVKnTRdUJDQ3njjTeoWrUq9evXZ/Lkyfzyyy+sXbv2sv8dRERyio8/NjMd+/VK4OGv28L27eYX/bNdyMV97rgD1q2DJ58sTvWfvuC1Zl1489/nyHPvvdC+vZl/WrYsEdGxGV6fml20qxuY/s+wfr3ptLp2LXToYF6saXpBzJ8PDz9sZvDOnm12eJFsztvbLPZt3Roee8zMDjjbu0MkIxRGM0HlypWJiopi0KBBPP744+zfv5/ixYtz9913M3XqVMqUKZPq/r/88gt169ZNdez+++/n3XffxdfXlz59+rBr1y7y5MlDxYoVGT58OL169bpsDbfddhs33ngjLVq04MSJEynXu5yAgACWL19OWFgYrVu35uTJk5QrV44777wTb29vAF555RX27dtHt7Nd8B5//HE6derEli1bLnndIUOGULhwYQYPHszu3bspWbIkTzzxBAB58uThgw8+YNCgQQwcOJCmTZuyePHii67x4osvcuzYMV599VX+/fdfgoKCmD59OnXq1LnszyQikhP8/LOZCXfPXcm8s+txWLXKTM1VEHVMiRImUL33HoSFNePrktEsfGYkFb4cBNWrs7H7y7ye72aOJZslI+e2QQGyfSBN5dQpePttMzpcpIjZWqhjx1R3SUqC8HCz21CNGia7VKniUL3ievnzw6xZpqtax47mf4zbb3e6KskmLFeu6bsWwcHBdlRUVLrntmzZcsmutHJpnTt35uDBg8yaNcvpUrI8vcZEJKvbtg1uvtkMMv3W8hW8R4+EkSPNnF3JEqKizIjfjh0w8sWdvPjXC1izZ7GleAX639mD38qc/3cm0N+H5f1aOlitC61aZUZDN282nYXfe8/siXOBTZuga1dz14cegk8/hQv6K0pOcuiQ+YBs2zaIjDy/lY8IYFnWWtu2g9Me15pRERGRLOro0fMz3hZ3/NAE0RdfNNucSZYRHGz2fX3kEej9fgVaHPuRbq0HUvDkcWZ8Fcqn0wbS+O/1YNuX3wYluzhxwnRRveUW8yKdPdt0V70giCYmwltvQb16Jpt88425KYjmYEWKmLnYZcrAPfeYT2lErkDTdEVERLKgpCQTbv76C3578weKv9HLtB4dOTJX7W2ZXRQsCJMmmfWkPXpYLE8OY+Edt/HCofE88dssvpnyPzaWrMz02x6A03fA2QZ+ruC2tanHjpkfcsQIs2b52WdNw6I0HYh++80MmK5bZ0aMP/gAzm4rLjldyZKmy1rTpqaT8pIlZnsfkUvQNF3J1fQaE5GsKjQUhg+H6X1X0+GD5lCrFixalOO3D8kJ/vwT7rrvNNv/yItf/R2UbLKe9n/+TPeoCKoc/AcCAswi4O7dzWjSdYiIjk13u5VUXW6v1+bNpiHRl19CfDzUr29enGnWLJ88CQMHmvWhJUrAJ59A27auKUGymW3boFkz86naL79A1apOVyQO0zRdERGRbGLiRPO7/huPbafDhDZmwejMmQqi2US1arD597y0eTSeY2srsuPzEGbkf5Z1c5abDsg33ghhYVC2rAmlW7de83OFR8akCqIACYlJhEfGXN8PkZhommS1aGHq/fRT0yn4119hzZqLgujy5VCnDgwdCk8+afKrgmguVrkyLFhgwmirVvD3305XJFlUlg+jTo/cSs6l15aIZEUrVpgBs3bNDjFg9d1m78a5c81Qk2Qb3t7w41e+/PILNKiZj20zqvPmU2WYndQaO3KemcP64IMwdqxJr+3bw7JlcJX/Nl1qDeo1r03dtw8GDYIKFeCBB0xXpqFDYfduMzLasGGqaeLx8dCrl5mVefKk6Vvz2Wfg739tTy85SPXqZg3psWMmkO7d63RFkgVl6TDq5eVFQkIOWOgvWVJCQgJeXl5OlyEiuVhEdCxNhv5MxX6zaTL0Z8bP2Uf79lClzEm+TWyHtXMHRERAUJDTpco1uvVW8wHDtGlmF5R774WWLWHNqZtgwgQzYtS/PyxdahJdw4YwZYpJdhkQ4O9zVcfTZdtmKuXDD5vR2jffNNPCf/zRTLfs2zfdRZ8LF5q7ffAB9OgBGzfCnXdm/GklF6hTx3yY9u+/ZruXAwecrkiymCwdRkuUKEFsbCwnTpzQKJa4jG3bnDhxgtjYWEpopEFEHHJurV9sXAI2sGv/aV7okp/jx8+w4oan8Fr5ixmJatrU6VLlOlkW3H+/mbo6erTZ7uTmm02Dqu0JpWHwYNi1Cz7+GI4cMSd8fU3jl8ceM4sw58+H/fsvunZoSBA+Xp6pjvl4eRIacoUPMI4dM+lx7FgTGJo1M8OaL7xgFr3+9JNp5ezpedFDjxyBp5822cLLy+ToDz9Up1y5hEaNzDKD7dtNU6O4OKcrkiwkSzcwAjh69Cj79+8nMTHRjVVJTufl5UWJEiUomKYDoIiIuzQZ+jOxZ6dS2jYcjKjHib9KMe7GLjy98QvTpfTVV50tUjLF0aMmX44caZZm9ugBr712dmeU5GSYN88swvz9dzOdd9eu8w8uXRpq1zYB8uzXiPj8hC/Yer6b7p3VaFcmrxl1TXv75x/z9fDh89esU8cU8eijl12XHBsLn39uMvO//0KfPjBgAPhcxSCs5GJz58J995m9kObN06cXucylGhhl+TAqIiKSE1XsN5tz/wLH/VKVIyuq0TtoAO/FDDRbZowZoy1ccrg9e0yY++wz8PODfv3M+suLwt1//8H69efD6e+/m2HWcx/U+/iYUdTChc8HzrTLnPz8oHz5i29BQSbUXuK1lpRkBknHjYNZs0xWbtUKhgyBBg1c/B9Ecr4ZM8xa6WbNzP60+iQj11AYFRERyULOjYwe/6MUB3+oT/vyE/nuny78Wq0BTTYugzzaCjy32LzZBNGZM6FMGTNr9/HH050he97p07BlS+qAGh+ffuAsV850FLqKDzd27TIh+bPPTO+iEiXM3qHduplGqSLXbPJkeOIJuPtuE05duOeuZF0KoyIiIllIRHQsL4/dwc7PG3FzoWUsPHIXO4uWYfu02bRpUs3p8sQBS5aY/WXXrDGDnM2bm1HIVq3MAGZmD5SfOQNz5phR0LlzzfTxO+8060PbtFFmEBcaO9bMAHngAfj6a334lgtcKozqb15ERMQBLSoFcuqnElT02sEPJ+7nWIGC/DNxqoJoLnbbbbBqlWli++OPplvt99+bcwEB54Npq1ZmBNVVdu40I6Cff26mDpcubbZB7doVKlZ03fOIpHjmGTh+HF55xXTBmjBBn3bkUgqjIiIibpacDE8+Ccd3nyKmQkd8DybBsqXcdeONTpcmDrMs0+PlvvvMyOT27SaULlxoRisnTTL3q1btfDBt0QKKFEn/eqdPw6FDZtnpf/9d/Offfze9ZADuugs++gjuucfkA5FM9fLL5gUaFmY6Yk2fDoUKOV2VuJmm6YqIiLhARHQs4ZEx5zuahgTRrm5guvd95x14rX8yf9XsQOXNM03K0AaNcgXJybBhw/lwunSpWSZqWVC3LlSpcj5snvsaH3/p63l5mW1FH3vMrActX959P4tIii++MHPBa9Qw88QD03/flOxNa0ZFREQyybk9QxMSk1KO+Xh5MqRDrYsCaWSkGYGKqB5G281DYdQoePFFd5csOUBiIqxefT6c7ttntocpUsR8vfDP6R0rUEANmyWLmDfPbMbr728+nKtZ0+mKxMUURkVERDLJhXuGXijQ34fl/VqmfL9jh9li75n8k3hn9xPQvTt88okSgaS4mhF2kRzl999Nh90TJ8xi6RYtnK5IXOhSYdTDiWJERERykj3pBNG0xxMSoEMHqH96JW//2820Sh09+pqCaER0LE2G/kzFfrNpMvRnIqJjr7V0yULOjbDHxiVgA7FxCYTN2KC/X8kd6tSBX38103RDQkyXXcnxMhRGLctqZlnWj5ZlxVqWZVuW1TnNecuyrAGWZe2xLCvBsqzFlmWpC4OIiOQKAf7pb9x+7rhtm10M/vv9H2bmaYdVrixMm3ZNXWIUWHKu8MiYVFO9ARISkwiPjHGoIhE3K1cOli2Dxo2hUyd4913zBio5VkZHRn2BjUAvIL2Pf18FXgFeABoA+4H5lmX5uaJIERGRrCw0JAgfL89Ux3y8PAkNCQLg449h+pfxrC7ZFu/kkzBzplm0dw0UWHKujIywi+R4hQubNaQPPQR9+8ILL0BS0pUfJ9lShrZ2sW17DjAHwLKsLy48Z1mWBfQGhtq2Pf3ssScxgfRRYKzryhUREcl6zq3pS2+t34oV8FKvZBaXfJySBzbA7NlQvfo1P5cCS84V4O+T7trjS428i+RY3t5mmm65chAeDrt3m+/z53e6MnExV+wzWhEoBcw7d8C27QTLspYCt6AwKiIiuUC7uoEXNZrZtw86doT3fF+n8b8R8N570Lr1dT2PAkvOFRoSlG5X5nMj7K6kRkmS5Xl4mGm6ZctCr17QsqWZVVK8uNOViQu5ooFRqbNf/01z/N8LzomIiOQqiYnw4IMQcvArno97B7p1M79QXacrTQmW7Ktd3UCGdKhFoL8PFqYbc3rbA10vrTuWbOWFF2D6dFi3Dm65BbZtc7oicSFXjIyek3Z1sZXOMXPCsroD3QHKlSvnwhJERESyhtBQOPXLKj716gq33QYffeSSLVwuNyVYsr/0Rthd7XLrjvU6kiypfXuzmW7btqa50axZcPPNTlclLuCKMLrv7NdSwK4Ljpfg4tFSAGzbHgeMA7PPqAtqEBERSZcT0xG/+Qamj9rFpvz34Vkq0HTOzZvXZdd3R2CRnEvrjiVbuuUWWLHCLHVo3hymToU2bZyuSq6TK6bp7sAE0jvOHbAsKx/QFFjhguuLiIhcEyemI65fDy90Oc7CAm3x8zwBP/4IxYpl2vOJXK0rbUUkkmVVqwYrV8KNN0K7dvDJJ05XJNcpo/uM+lqWVceyrDpnH1Pu7PflbNu2gfeBfpZldbAsqybwBRAPaLdaERFxjLu3QYmLg/vbJzORJ6iasB5ryhTzS5NIFqJ1x5KtlSwJixfDXXfBc8+ZNREnTzpdlVyjjI6MBgPRZ28+wMCzfx509vy7wEjgIyAKKA3cadv2MZdWKyIichXcOR3xzBl49FF4aueb3HNyBlZ4ONx9t8ufR+R6uatRkkimKVAAIiLg2Wdh+HCzXdZ334Gt1X/ZjWU7/JcWHBxsR0VFOVqDiIjkTE2G/pzuNiiB/j4s79fSpc/Vqxcc+OBrvqYTdOkCn37qkoZFIiJyGQsXwssvmzUSt95qttAKDna6KknDsqy1tm1f9BfjijWjIiIiWZK7piOOGQPrP1jEl55PQbNm8PHHCqIiIu7QqhX89huMGwd//gkNGsCTT0KstirKDhRGRUQkx3LHdMR582DcCxuYlacdnkFVzNQxF3bOFRGRK/D0hKefhr/+gr59YcoU0+xo8GA4ccLp6jLXoUOwdKnTVVwzTdMVERG5Rps3wwMN/2HRycYUK2Hh8etKKFvW6bJERHK37dtNKJ02DcqUgWHD4JFHcs6MlYQEmDkTvv4a5swxa2j37QNvb6cruyRN0xUREXGhAweg012HmJFwF8V84vH4aa6CqIhIVlCpkmlotGQJlCgBnTpB48bw669OV3btzpwxU3GefNL8TA89BKtXwwsvmHWz2XRGjsKoiIjIVTp1Ch5ud5LRu+6jqsdWPH6IgFq1nC5LREQu1KwZrFkDEybAP/+YQProo+bP2YFtw6pVpkNemTIQEgI//GCC6MKFsGsXjBgB9epl21FfhVEREZGrYNvQvWsSz614jCb2MjwmfQktWjhdloiIpMfDAzp3Ns2NXnsNvv8egoLgjTcgPt7p6tIXEwNvvmnWvTZqBGPHmk7B06eb6biffgotW5q1stmcwqiIiMhVGPKOTfBXvenIdBg50nxCLSIiWZuvr2loFBMD7dubP1erBp99ZgKe0/bsMf+mBAfDDTeY+sqXh88/N/VNmwYdOkC+fE5X6lJqYCQiIpJB06bB6gfe5V36Yr/0MtbIEU6XJCIi12LlSnjpJTMNFqB0aahb10x5rVvX3CpUyJzpr4mJsGOH6f77558wezb8/LOZehMcbKYSP/QQBAS4/rkdcqkGRgqjIiIiGbBmDXzcZDKfJz5O0gMP4znlKzP9S0REsqfkZFixAqKiIDra7Fe6ZQskJZnz/v4XB9SgoIxNjz1zBv7+2wTOtLedO88/B0CVKiaAPvqouX4OpDAqIiJyjXbtgj615zP58N3YTZqSd+HcLN1CX0RErlFCAmzcaILpuYC6fr3pXAeQPz/cdNP5gHrTTXDkyMWBc8cOMwJ6jq8vVK2a/q1YsWzbgCijFEZFREQyKCI6lvDIGPbEJVDSx5diE3z45u9WeFWtSL41v0ChQk6XKCIi7pKYCH/8YcLpuYD6++9w9Gjq++XPb0Y50wucJUvm+MB5OZcKo3mcKEZERDLXhWEqwN+H0JAg2tUNdLqsbCEiOpawGRtISEzCToYjE3z4/u/7SCpcCL9FcxVERURyGy8vs31XrVrwxBPmWHIybN8OGzZA4cImcAYE5OrAeS0URkVEcpgLwxRAbFwCYTM2ACiQZkB4ZEzKfzuPBSWI+PsxfPIcp8eT7/NVoP77iYgIpmdAlSrmJtdMnRdERHKYC8PUOQmJSYRHxjhUUfayJy4BgDO/FeOr6F6Ut3bS7aHXWOFd0uHKREREchaFURGRHOZcmMrocUktwN+H0+uLM2b+2zTiV3q17UNUmRsJ8PdxujQREZEcRWFURCSHuVRoUpjKmIZHqjJxbn/aMpPXWj5P5A234OPlSWhIzmy3LyIi4hSFURGRHCY0JAgfr9R7oClMZcxn7/7HQ0M70oJFDOzQj28a3EWgvw9DOtTSelsREREXUwMjEZEc5lxoUjfdqzP+zd00GXQnVTy2Y387gzfvb8ubThclIiKSgymMiojkQO3qBip8XoVPXv6Tu967g+J54vD4KZI8rW5zuiQREZEcT2FURERyLduGT7r/xv2ftiafN+Rduog8N9dzuiwREZFcQWFURERyJduGsY8uodOUNiQWKEyBNfPxrF7N6bJERERyDTUwEhGRXMe24dO2P9J5SgjHC5eh8OblCqIiIiJuppFRERG5LhHRsdmqWVJyMkxs9SVPLe7CnpL1KLthDlbxYk6XJSIikusojIqIyDWLiI4lbMYGEhKTAIiNSyBsxgaALBlIk5JgauP3eWrNS2yt0IrK677HKujndFkiIiK5kqbpiojINQuPjEkJouckJCYRHhnjUEWXlnja5oebXufRNS+xpcb9VPljtoKoiIiIgxRGRUTkmu2JS7iq4045nZDEwht60GHzW6y7uRvV108Fb2+nyxIREcnVFEZFROSaBfj7XNVxJxzde5xfK3ei9Y6PWXt7X2r/Og48PZ0uS0REJNdTGBURkWsWGhKEj1fqYOfj5UloSJBDFaW2ZuhC4srWotneqfza4V3qzx8KluV0WSIiIoIaGImIZLrERFi3DlauNLfVq83WIiVKnL8VL576+3PHihUDLy+nf4JLO9ekKKt10z3y92E23PEKt/41gZ15q7J59BIaPdvM0ZpEREQkNcu2bUcLCA4OtqOiohytQUTElfbtOx88V66EqCg4edKcCwyERo3McsUDB2D/fnM7cADOnEn/ekWKmHBatSrccQfceSdUq6YBvkuJfm06AUN6UjT5AMsahdJozhvkK5x1pg2LiIjkNpZlrbVtOzjtcY2Miohch8RE+P3388Hz119h505zzssL6tWDZ5+Fxo3NrWzZ9K+TnAxxcakDatpbdDTMnGnuX66cCaV33gmtWpnAmtsd+WMvW1v3pP7fM9icry4Hxs2h+eN1nS5LRERELkEjoyIi1+D0aRg5Et5+G+LjzbHAwPOhs3FjqFsX8uW74EG7d8OyZbBihRkqLVDA3Hx9r/znAgXA359tOz2ZPx/mzYOFC+HoUTNC2qDB+XDaqFHWntrrcrbNhpc+p+wHffC2T7Kk+QCaz3yFfL76vFVERCQruNTIqMKoiMhVWrAAevaEmBho2xYee8yEzzJlLrhTcjJs2WLC5y+/mK9//23OFSgAfn5w/LhJshl9H86fH4KDoWFDaNiQM/UbsmZvGebNM+H011/N0/r5QYsWJpiGhECVKi7/T5BlHI3exu57ulNj78+syd8M74njualjNafLEhERkQsojIqIXKfdu+GVV+Dbb6FyZfjgA7j77rMnT52CtWvPB8/ly+HwYXOuVCm49VZza9oUbroJ8pwdtbNtM0p6/Pj5cJrmz9FbdrPktx34791Fg3//4oZ/t+OZeNo8PiAAbr4ZGjYk/saGLDoWzJxf/IiMhB07zF1q1IAOHaB9ezNamyPWmp45w5ZnR1Hh89dJtPPwc+tw7prxNN4+ahIvIiKS1SiMiohco9OnYdQoGDgQkpLgf/+D0FDId+wAfPghLF5sWuSeOmUeEBR0PnjeeitUqnTNCTAiOpawGRtISExKOVbQSuLD6ha3xe2AVavMc2/dak56eJj02bAh+ys2ZN6xxnz+aw2W/OJBcjKULw/t2plgeuut2XO7zaPL1nOwXVcq/RfFYr97KfzNx9S+p8yVHygiIiKOUBgVEbkGixZBjx5mxm3btvD++1CxTCJ89BEMGGBGL+vXPx88mzQxe7K4SJOhPxMbl3DR8UB/H5b3a3n+wH//mVC6atX527mR2WLFONXoNn4r1ILJsS34bEV1Tp22KFbM/Ezt28Ptt6dZ35rVnDjBv1/M5fDYb6m8fgaHKczi9h9w39cP4Z0vJwz1ioiI5FzqpisichX27IE+feCbb6BiRdPF9t57gchIuLs3/PGHWZD53ntQvXrm1ZFOEE33eNGicNdd5gZm+u/WrWa68KJFeC9aRONd02kMfFiiJLurNGdeYgvGTG3B559XxdfX4u67TTC9+24oWDDTfqSMO3mS2M9+Im7ct1Tc+CMlk49jUZxxBbow9eGOvNyjhoKoiIhINqYwKiJygcREM/P2zTfNn998E/r2BZ/YrdD2ZZNKq1QxX++5J9MXYAb4+6Q7Mhrgf4V9My3LbExatSp07mzC6fbtsGgRHosWUW7RIrrtnUo34GTRANYVbsGUn1rwv29b8ESeitSpa9GwYcpyVKpUMTOAM5t98hS7Po0kbty3VNr0I4HJx/CmKDOLPMy3FZoQFVwGq/Bp4DRhMzYA0K5uYOYXJiIiIi6naboiImctXWqm5G7caEYHP/gAKpc4Bm+9ZUZAvb3h9dehVy/zZzdIb82oj5cnQzrUur4QZtvw559mveuiRea2fz8AcQXLsiXPTaw6Vp0NiTewhersLXgD1RoVSRVQXTUb2T51mh3j5nNk/LdU2RSBX/JRDlGYFSU7kNjhIRqENuehqb9kbLqyiIiIZDmapisicgm2DYMGmSWg5ctDRAS0vTcZa/Ik6NcP9u0zo4tDhpjOuG50LnCGR8awJy6BAH8fQkOCrn800LJMo6WgIHjmGfMfYcsWWLQI/19+ofGWLTSKWYDF2aZMR+G/n0uwcZ4Jpwu4gcMlq+PX4AYqNy/LzY08KFoUzpy54JZokxx/Ao4dwz56DCv+GNaxo+Zr/DE8jh/DWvUrVTdHUCk5jsP4s6z0/SR1eJDgvq24t+z5zVIzPF1ZREREsg2NjIpIrmbbJoQOGgRPPgljxkD+DavgxRdNQ6CGDc0Q6c03O12q+yUlmb1Rt2wxa2S3bCFp8x8kb9qC19FDKXc7Tn7+pBqn8KYgR/HjWMrNk+TLPsURCrI6oB3J9z9Ivb53UDwwb7r3y3AjJxEREclyNDIqIpKGbZs1oYMHQ5cuMH7QXjye6wdffmlGQCdOhMcec89iyazI09NsS1OpklkfC3gCnrYNBw+mBFR77R8ErP2D5MRkkvKXJSm/H4cL+PFfAT+SfQti+/ph+/qBn5/pjOTnh1XQD49CfhS7oRh3lEg/gF4oNCQo3enKoSFBmfXTi4iISCZTGBWRXOnCINq1K4xrOQWPG542m4r27Qv9+5vwJBezLLNgtHhxaNoUX8A3k58y06Yri4iIiGMURkUk17FteOMN05eoa1cYV38sHo89Z/YInTDBtI7NgIjoWIUjN2pXN1D/fUVERHIQhVERyVUuDKLdusHYKu/i8XxfMw31u+/A5wpbppyVtsttbFyCthoRERERuQq5dCGUiORGtm12ZnnrLXi6m8244v3x6NcXHnoIZszIcBAFM130wvWLAAmJSYRHxri6bBEREZEcSSOjIpKjXGrq7Lkg+vbb0L1bMp9498IaMtoMj37yiWnWcxW01YiIiIjI9VEYFZEc41JTZ20b1s4I5O234dluZxhzqivWp1/CK69AeLhpyHOVAvx90t1qJMA/46OrIiIiIrmZpumKSI6R3tTZE6eTeOGVRN5+G57veooxBx/EmvSl2Vj0GoMomK1GfLxSj6ZqqxERERGRjHPZyKhlWX7AYKA9UAKIBnrZtr3GVc8hInI5aafI2jbE/RLE0ZUVeKHLcUb90wFr/jx4/33o1eu6niunbjWiDsEiIiLiLq6cpvspcBPwJLAbeAxYYFlWDdu2Y134PCIi6bpw6qxtQ9zSII7+WoWq9TYzKqY71sqV8Pnn8NRTqR53rQEsp201og7BIiIi4k4umaZrWZYPcD/Qz7btxbZtb7VtewCwFXjOFc8hInIl56bOXhhEq9T6jdXxD2KtXg1Tp6YbRMNmbCA2LgGb8wEsIjr3fYamDsEiIiLiTq4aGc0DeAIn0xxPAG510XOIiFzWudG7F0NPcfTXStSps4alhx/Bb/8e+PFHaN36osdcLoDlttFAd3YI1nRgERERcUkYtW37mGVZK4HXLMvaCOwDHgEaY0ZHU7EsqzvQHaBcuXKuKEFEBACf/YHsWgih7bcy7LcHsA4fhshIaNo03fu7e4uWrBzC3NUhWNOBRUREBFzbTfdxIBmzXvQU8CLwDZCU9o62bY+zbTvYtu3g4sWLu7AEEcnN9u6Fxx+HdlU2MmxFU6zjx2HRoksGUbh00MqMLVqy+pRgd3UI1nRgERERAReGUdu2t9m2fRvgC5S1bftmwAvY4arnEBG5lKQk6NQJCsT/y7fHWmN5esDSpVCv3mUf584tWrJ6CGtXN5AhHWoR6O+DBQT6+zCkQy2Xj1a6ezRaREREsiZXdtMFwLbt48Bxy7IKAyHAq65+DhGRtN55B35ZlMiuag/gtesQrFgB1atf8XHu3KIlO4Qwd3QIdtd0YBEREcnaXLnPaAhmpPUPoAoQDsQAE1z1HCIi6VmyBAYMgDlBfSgV8wtMngx16mT48e7aokUhzAgNCUq1ZhQybzRaREREsi5XrhktBIzGhNEvgWXAnbZtJ7rwOUREUjlwAB59FF4uMZmQmA+gVy8zXzcLcueU4KzMXdOBRUREJGuzbNt2tIDg4GA7KirK0RpEJHtKToZ774X/FkSz0uMWPBo1hPnzwcvL6dIuKSt30xURERHJDJZlrbVtOzjtcZevGRURcZeRI2HV3P/YUbQDHvmKwtSpWTqIgvumBIuIiIhkdQqjIpIt/for9O+XxOoSj+AXtwfm/AIlSzpdloiIiIhkkMKoiGQ7hw/Dww/D+wX6U3v/fBg/Hm6+2emyREREROQqKIyKSLZi29CtGzTcNY3nkodB9+7mgIiIiIhkKwqjIpJlpdfsZ8/KQP6YsYnovJ2hXiP44AOnyxQRERGRa6AwKiJZUkR0bKq9KGPjEnhpzE6OfZGfdfnb4+XnC9Omgbe3w5WKiIiIyLVQGBWRLCk8MiYliAIkn/Jk37RafM8jBJzegfXdzxCorrQiIiIi2ZXCqIhkSXviElL+bNtwaF4twuJGcg9zzNTcpk0drE5ERERErpeH0wWIiKQnwN8n5c/HN5Sh+eZoBjKQuXXvgJ49HaxMRERERFxBYVREsqTQkCB8vDw5fdCXIvPzMtnqxKaSlUn8aAxYltPliYiIiMh10jRdEcmS2tUNJDERerTx5PuklpDPIvazybRtXMXp0kRERETEBRRGRSTL+mt+AB/uf5gaHjFYP0Zy5+2NnC5JRERERFxE03RFJEvatg3+Gvg1D/Et1uDBcPvtTpckIiIiIi6kkVERyXJsG17rsocxiS9wun5j8vbt63RJIiIiIuJiCqMikuVM+cbm0aXP4OeVQJ6vvwBPT6dLEhEREREXUxgVEbeIiI4lPDKGPXEJBPj7EBoSRLu6gRfd7/BhWP7sJEYzi+QhI6FaNQeqFREREZHMpjAqIpkuIjqWsBkbSEhMAiA2LoGwGRsALgqkw16M5a1jLxJfpwm+vV90e60iIiIi4h5qYCQimS48MiYliJ6TkJhEeGRMqmMrlts0ndyd/HlO4/vtBE3PFREREcnBFEZFJNPtiUu44vHERJj70Bfcwxzsd4ZC1aruKk9EREREHKAwKiKZLsDf54rHx7+5mz6xvfnvxmZ4v9LTXaWJiIiIiEMURkUk04WGBOHjlXrKrY+XJ6EhQQDs2G5TdVg3vD3PUPSHz8FDb00iIiIiOZ1+4xORTNeubiBDOtQi0N8HCwj092FIh1q0qxuIbcMPbT/jjuRIEga8C5UrO12uiIiIiLiBZdu2owUEBwfbUVFRjtYgIs6ZNeYfmvWoydEq9SkTs1CjoiIiIiI5jGVZa23bDk57XL/1iYhj4g7b+L3UjTweyZSeq+m5IiIiIrmJfvMTEcfMaT+e207P58Cr4XhWqeh0OSIiIiLiRgqjIuKI6O930mbJK8SUbUX5t59xuhwRERERcTOFURFxu8RTyZx+oiuWBYFzP9P0XBEREZFcSL8BiojbLX5kLA3jf+avZ0fge2N5p8sREREREQcojIqIW+3+ZQeNvw8luvgd1P3oaafLERERERGHKIyKiNvYSckcateFJDwp8eNnYFlOlyQiIiIiDlEYFRG3WffMGG46tJhVD44ksFFZp8sREREREQcpjIqIWxxbv4Nqn/dluV9rWk7u4nQ5IiIiIuIwhVERyXy2zZ77niPJ9qDAV+PI46XpuSIiIiK5ncKoiGS6PaO+JWhnJDMbvk2dNpqeKyIiIiIKoyKS2eLi8OnXi9886tNyeg+nqxERERGRLEJhVEQy1T+dwih46gAbeo6jVKCn0+WIiIiISBahMCoimebMLyspN+cTvvTvxSPh9ZwuR0RERESykDxOFyAiOVRiInEPducEZSk1bhB58zpdkIiIiIhkJRoZFZFMET9oJMX2beTzeqNp3dHX6XJEREREJItRGBUR19u+nbxDBvK91Z6Hv2qLpZ1cRERERCQNhVERcS3b5uhjz3MqyZP1XT/ghhucLkhEREREsiKtGRURl7KnfkvBlZH8r8AoXg0v43Q5IiIiIpJFKYyKiOvExXHy2V5soj6VRvTA39/pgkREREQkq1IYFRGXSXz1f+Q9coCRQXOY1E17ioqIiIjIpWnNqIi4xsqVeI7/hFH04rnx9fBUFhURERGRy1AYFZHrl5jI6S7PEEsZ1t8/iKZNnS5IRERERLI6TdMVkev33nvk/WMDL+f9gZHvaU9REREREbkyhVERuT47dpD0xgB+oD21+relbFmnCxIRERGR7MAl03Qty/K0LGuwZVk7LMs6efbrW5ZlKeyK5GS2jf3c85xM9GRYwAeEhjpdkIiIiIhkF64Ki32BHsCTwAbgJmAicAoY7KLnEJGs5rvvsCJ/IoxR9Hm/DD4+ThckIiIiItmFq8LoLcBM27Znnv1+p2VZPwINXXR9Eclq4uJIfrEX6z3rs7FJD0Z1dLogEREREclOXNVNdxnQwrKsGwAsy6oBtATmuOj6IpLV/O9/sH8/3ZLH8d4HnliW0wWJiIiISHbiqpHRYYAfsNmyrKSz133btu0x6d3ZsqzuQHeAcuXKuagEEXGblSuxP/mED+hNg2fqUbu20wWJiIiISHZj2bZ9/RexrIeBcCAU2ATUAUYBobZtf3a5xwYHB9tRUVHXXYOIuMnp09j16nFw6xHq5dtM9FY/ihVzuigRERERyaosy1pr23Zw2uOuGhkNB4bbtj3l7PcbLMsqD4QBlw2jIpLNhIdjbdrEU8wk9F0FURERERG5Nq4Ko/mBpDTHknDdmlQRyQr+/BN78GDmFHiQHeXv5bnnnC5IRERERLIrV4XRmUA/y7J2YKbp1gVeBr500fVFxGnJydC9OyctH7oeH8XXH4KXl9NFiYiIiEh25aow+gJmP9ExQAlgLzAeGOSi64uI0z7/HJYs4RWv8TTtWIqWLZ0uSERERESyM5c0MLoeamAkkg3s2wfVq7PJqzYNji3ijxgLNcIWERERkYzI7AZGIpKT9epF0vEEOiSOpd1zx3jk6yj2xCUQ4O9DaEgQ7eoGOl2hiIiIiGQzCqMicnmzZsG33/JR8cEcyVuZtYV+5lRcIgCxcQmEzdgAoEAqIiIiIldF3W5F5NKOHYPnnuNg6Zr0OfAqhVtu5pSdmOouCYlJhEfGOFSgiIiIiGRXCqMicmmvvYYdG8vDR8bRMiQvCaV3pXu3PXEJbi5MRERERLI7hVERSd+qVfDhh/wc9DxLExszahQEFvZJ964B/ukfFxERERG5FIVREblYYiI8/TSniwXQ/o936N0bgoIgNCQIHy/PVHf18fIkNCTImTpFREREJNtSAyMRudjw4bBhA2FVfsA3T0Fef90cPtekKDwyRt10RUREROS6KIyKSGp//QUDB7KjfkdGrm3LpEng53f+dLu6gQqfIiIiInLdNE1XRM6zbXjmGWzvfLTZ8QFNmkCnTk4XJSIiIiI5kUZGReS8CRNg0SKmthzL5kWlWfshWJbTRYmIiIhITqQwKpKLRUTHpqz/vNEzgRmjXyaxXlMeX9yNZ56BunWdrlBEREREciqFUZFcKiI6lrAZG0hITAKg+/cfwvHjPHxqNAX9PXjrLYcLFBEREZEcTWFUJJcKj4xJCaLNt62h7ZalDA16hlmbbmLMGCha1OECRURERCRHUwMjkVxqT1wCAPlPJ/DWvDH8WaQcA2OHkbfEEbp3d7g4EREREcnxFEZFcqkAfx8AQpd+SZmjB+gREM7J+EJUa78VT0+HixMRERGRHE9hVCSXCg0JotWudTy1diaf39iOBVvup2DNPQx+ppTTpYmIiIhILqA1oyK5VLty+QiZ/wE7SpSjz4n38PC0eW+ERbu6AU6XJiIiIiK5gEZGRXIj24bu3fGJO0TsaxEc3lGBIYPz0OXO0k5XJiIiIiK5hEZGRXKjL76AGTM4NWgYT7xXl2rVoHdvp4sSERERkdxEYVQkt9m2DV58EZo3p+/+V9ixA5Ysgbx5nS5MRERERHIThVGR3OTMGXj8cfD0ZE3PiXzwgCc9ekCzZk4XJiIiIiK5jcKoSG7yzjuwciWnJ37DY/8rR7lyMGSI00WJiIiISG6kMCqSW/z6KwwaBJ068cbmh/nzT5g3D/z8nC5MRERERHIjhVGR3CA+Hh57DAIDie72EcNvhy5d4I47nC5MRERERHIrhVGR3KB3b9i+ncT5i3nyxUKUKAEjRjhdlIiIiIjkZgqjIjnd99/DZ59BWBjvLGvGhg3www/g7+90YSIiIiKSmymMiuRke/bA009DvXps7DiAtxvBo49C27ZOFyYiIiIiuZ3CqEhOlZwMTz0FJ05wZuJXPPVUXvz9YdQopwsTEREREVEYFcm5Ro827XLHjGHknBuIioKpU6FYMacLExERERFRGBXJmTZtgldfhXvuIabFs7xRB9q3hwcecLowERERERFDYVQkpzl1yiwMLViQ5PGf0fUBi/z54aOPwLKcLk5ERERExFAYFclpXnsN1q+HmTP5aFpJli+HL76A0qWdLkxERERE5DyFUZGc5OefzQaizz7Ljhvvpd9D0Lo1PPGE04WJiIiIiKSmMCqSU+zZA48/DtWqYQ8fwdP3gacnjB2r6bkiIiIikvV4OF2AiLhAfDxxrUI4cfAQd9/Sg6pP/M3ChRAeDuXKOV2ciIiIiMjFFEZFsrukJPbe0wG/mM30aNuX9flqsH1WFfKX/4/iwbFOVyciIiIiki6FUZHsrk8fSi+dz8BWT/NzpQYcmlcTkj3wD1nPiPkxTlcnIiIiIpIuhVGR7Oyjj+D99/m8flu+rN+G45sDSNhWEv9mf+BV+AR74hKcrlBEREREJF1qYCSSXc2ZAy++CG3aMKFRTxL/8eTwghvxDjiMX/2dAAT4+zhbo4iIiIjIJSiMimRH69bBQw9B7drw9df0WH+Uzu19ASh67+9YHuDj5UloSJDDhYqIiIiIpE9hVCS7iY2Fe+6BQoVg5kzsAr788KEvpw/Y3NB5HQmFTxDg70NoSBDt6gY6Xa2IiIiISLoURkWyk/h4aNMGjhyBZcsgMJARw2HKFBgyxKJfvzpAHYeLFBERERG5MoVRkewiKQkefdRM0Z05E2rXZv586NsXOnY0X0VEREREsguFUZHs4pVXTAgdPRruvpsdO+Dhh6FGDZgwASzL6QJFRERERDJOYVQkC4qIjiU8MoY9cQkE+Psw5vAKao8aBb16QY8enDgB7dtDcjJERICvr9MVi4iIiIhcHYVRkUyUNlRmpKlQRHQsYTM2kJCYBEC1tUupOX0we5vdQekRI7Bt6NoV1q83u7tUruyOn0RERERExLU8nC5AJKc6Fypj4xKwgdi4BMJmbCAiOvayjwuPjEkJojX+3c7oH4axpURFHm/ZCzw9GTnSNCx6+21o3doNP4iIiIiISCZQGBXJJBeGynMSEpMIj4y57OP2xCUAUPLYQT6bNpAj+Xzpcv8bbEuABQvg1VdNw6J+/TKtdBERERGRTOeSMGpZ1k7Lsux0brNdcX2R7OhcqMzo8XMC/H0oeDKez6YPxu/0Cbp2fIP9fkUpklyYhx6C6tXVsEhEREREsj9XrRltAHhe8H1pYC3wrYuuL5LtBPj7EJtO8Azw97ns496s6UOlEaGUO7SX7h36s6VEJbzx4vCPwWpYJCIiIiI5hktGRm3bPmDb9r5zN+Bu4CjwnSuuL5IdhYYE4ePlmeqYj5cnoSFBl37Q6tXc2bUd5U4d5eWuw1hSuQEBhXwoue5Wdv6Zl6+/hipVMrlwERERERE3cHk3XcuyLKArMNm27ROuvr5IdnGua26Gu+l+/z106gSlSpF3yRJG33ADo4GRI+GVn0zDorvucl/9IiIiIiKZybJt27UXtKw7gUigrm3bv1/p/sHBwXZUVJRLaxDJVmzbJM7QULj5ZvjxRyhRAoCFC+HOO82eot99p3WiIiIiIpL9WJa11rbt4LTHM6Ob7tPAmssFUcuyuluWFWVZVtSBAwcyoQSRbOLMGejRA/r0gQ4dYNGilCC6cyc89BDccIMaFomIiIhIzuPSabqWZZUA7gN6XO5+tm2PA8aBGRl1ZQ2SfUVEx2Z8SmtOcOwYPPwwzJljRkWHDgUPj5RT7dubrBoRAX5+zpYqIiIiIuJqrl4z2hk4BUxx8XUlh4uIjiVsxoaUfTlj4xIIm7EBIGcG0thYuPde2LABPvkEnnkm5dSBA2Zt6IYNZsZu1aoO1ikiIiIikklcNk33bOOibsAU27aPueq6kjuER8akBNFzEhKTCI+McaiiTLRuHTRsCFu3wqxZqYLoP//ArbfCpk3www9w990O1ikiIiIikolcOTLaHKgKPObCa0ousSed/TgvdzzbmjsXHnwQChWCZcugdu2UU5s3m2ZFx4/D/PkmlIqIiIiI5FQuGxm1bXuRbduWbdurXXVNyT0C/H2u6ni29Mkn0KaNmXe7alWqILpqFTRtCklJsGSJgqiIiIiI5HyZ0U1X5KqFhgTh4+WZ6piPlyehIUEOVeRCycmmQdFzz0Hr1rB0KQSeXwc7bx60agX+/rB8Odx0k3OlioiIiIi4i8KoZAnt6gYypEMtAv19sIBAfx+GdKiV/ZsXLV8OjRvD8OFmC5eICPD1TTk9darpY1SlirlrpUrOlSoiIiIi4k6u7qYrcs3a1Q3M/uHznO3boV8/+O47CAiAL7+Exx5LtVnomDHQs6eZkvvjj2ZkVEREREQkt9DIqIgrxcWZKbnVq8Ps2TBgAPz5Jzz+eEoQtW0YNMgMlN57L0RGKoiKiIiISO6jkVERV0hMhHHj4M034dAh6NwZ3nrLjIpeIDkZevWC0aPhySfh008hj/4vFBEREZFcSCOjItfDts1eobVqmTm3N90Ev/0Gn39+URA9fdrM1B09Gl5+2dxFQVREREREciuFUZFrtW4d3HGH2a7Fts3Cz4ULoU6di+56/Djcdx988w0MGWL6GXno/z4RERERycX067DI1dq7F7p2hbp1IToaPvgANm40ofSCBkXnHDhgMuu8eWYmb79+6d5NRERERCRX0SRByfYiomMJj4xhT1wCAf4+hIYEub4rb3IyrFwJ06ebRHn6tJlr278/FC6c7kPOnIGxY+H1183I6Lffwv33u7YsEREREZHsSmFUsrWI6FjCZmwgITEJgNi4BMJmbAC4/kB6+jT8/DN8/z388AP8+y94eUH79vDOO1C58iUf+ssvZgnp+vXQsqUZPL3xxusrR0REREQkJ1EYlWwtPDImJYiek5CYRHhkzLWF0fh4mDvXBNDZs+HoUfD1hbvvNiH07ruhYMFLPjw21uzs8s03ULas2Wb0/vs1LVdEREREJC2FUcnW9sQlXNXxdB08CDNnmgA6bx6cOgXFikHHjiaA3n475Mt32UucOgXvvw+DB5vpua+9ZtaGFihwFT+MiIiIiEguojAq2VqAvw+x6QTPAH+f9B9g27Bvn2k4tG6dGf1cutSsCS1XDp591gTQJk0yvO/K3Llm79C//oK2beG996BSpev5qUREREREcj6FUcnWQkOCUq0ZBfDx8iQ0JAgOH4ZNm2DDBhM+z90OHTp/gRo1ICzMBNB69a5qPu22bfDSS2ZQtWpVmDMH7rrLlT+diIiIiEjOpTAq2Vq7uoF4njhOxDcL8d/xF/WP7abVmf2UnLDVLOA8p2BBqFnTTL2tWRNq1TIdhYoXv+rnPHHC7BUaHm4GT4cOhd69wdvbdT+XiIiIiEhOpzAq2cO56bV//AFbtpivZ29tdu2izbn7eXub0c6WLU3oPBc8y5S57i5CBw7AjBnw9tuwaxc8+ii8+y4EungXGRERERGR3EBhVLKW5GSz+DKd0MmRI+fv5+sLN9wAt91mvt5wgwmelStneK1nRsTGmr5G06efX1papw5MngzNmrnsaUREREREch2FUbmiiOhYwiNj2BOXQIC/D6EhQde/h+eF/v3XdLGNjIT582H//vPnSpeG6tWhUycTOKtXN18DAzNtv5QdO0z4nD4dfv3VHKtRA/r3hw4doHZtbdUiIiIiInK9FEblsiKiY1M1CIqNSyBsxgaAaw+kp0/DihUmfEZGQnS0OV6sGISEnJ9iGxQEhQq54se4oi1bTPicMeN8OXXrwltvmX1Cb7jBLWWIiIiIiOQaCqNyWeGRMak61QIkJCYRHhlzdWF02zb46ScTPhctgvh4M532llvMIsyQEJP+PDxc/BOkLz7eNNadPduE0C1bzPHGjWH4cNNcV9uziIiIiIhkHoVRuaw96ezhebnjqSxfDl9/bQLotm3mWMWK8Nhj0Lo1tGhhutxmovh4EzQ3bYLNm83XTZvg77/NeQ8Ps/bz+edNAFUzIhERERER91AYlcsK8PchNp3gGeDvc+kHrVkDr79uQmj+/Gbabe/eZvSzSpVMWXB5LnReGDg3b4adO8/fJ29eM922cWPo1s2sA23a9Jp2dxERERERkeukMCqXFRoSlGrNKICPlyehIUEX33njRhNCIyKgaFGzEefzz5tA6gJJSWZEMyYG/vwz9dfdu8/f78LQ2bWr2U70xhvNtNvrabSb6Y2cRERERERyEYVRuaxzYeuyIeyvv+DNN2HKFPDzg0GDoFeva56C+99/JmBeGDZjYmDrVtP76Bx/f9PjqEUL87VGDdeEzvRkSiMnEREREZFczLJt29ECgoOD7aioKEdrkGv0998weDB88QV4e8OLL0JoKBQpclWXOXnS7OE5d67pcfTHH+fPeXmZmb3VqpnAee5rUJBpvuuuLVaaDP053enKgf4+LO/X0j1FiIiIiIhkQ5ZlrbVtOzjtcY2MytXbt890wB03znzfsyf06welSmX4Etu2mfA5d65prpuQYPJs8+bQpYsZ4axWDSpUcP0o57W4rkZOIiIiIiJykSzwa76403Wte/zvP3j3XfjwQzNftksXs0a0bNkrPvTECViy5HwA3brVHK9SxTQTuusuuO02ly0vdblrauQkIiIiIiKXpDCai1zzusekJBg2DIYONW1rO3Uya0SrVLns88XGwrRpJnwuWWKm4/r4mDWevXqZ3V2ucIks46oaOYmIiIiIyBUpjOYi4ZExqcIUQEJiEuGRMZcOo/HxZl/QH36Adu3grbfMHNrLiIuDIUNg1Cg4dcqs73z2WRM+mzUzgTS7yVAjJxERERERyTCF0Vzkqtc97t4NbdrA+vVmam7Pnpe9/unT8PHHpqfRoUPw+OPw2mtQter1Vp41tKsbqPApIiIiIuIiCqO5yFWte4yKgrZtzcjo7NlmWPMSbBu++w7CwmD7drj9drO0tG5dV1YvIiIiIiI5iYfTBYj7hIYE4ePlmepYuusep00z82m9vWHFissG0V9+gUaN4KGHoEABszXL/PkKoiIiIiIicnkKo7lIu7qBDOlQi0B/HyzMHplDOtQ6P/XUtuGdd+CBB0yaXLUKatZM91p//AH33Wcya2wsTJgA0dEQEuK+n0dERERERLIvTdPNZS657vHUKXj6aZg0yXTL/fRTyJfvorvt2wcDBpjT+fOb7Nq7d/ZsSiQiIiIiIs5RGBU4eBDat4dly2DQINN1yLJS3SU+HkaMgPBwk1t79DB3K17coZpFRERERCRbUxjN7bZsgXvvhT17YMoUs/gzjU2b4O674Z9/oGNHs21LdtkfVEREREREsiaF0dxs3jx48EEzHXfxYmjY8KK7rFhhsqq3txk4bdLE/WWKiIiIiEjOowZGudXHH5vhzvLlYfXqdIPozJlmm5ZixUwoVRAVERERERFXURjNbZKSoFcveP55uOsuM9xZrtxFd5swwSwjvfFGWL4cKlZ0oFYREREREcmxFEZzm/794YMP4KWXICIC/PxSnbZtGDoUunSBVq1g0SI1KRIREREREdfTmtHcZNIkGDYMnnsORo686HRyMrz8MowaBY8+akZH8+a9tqeKiI4lPDKGPXEJBPj7EBoSlP6WMiIiIiIikispjOYWK1dCt27QsqVJm2mcPg2dO8M335h9Q0eMAI9rHDePiI4lbMYGEhKTAIiNSyBsxgYABVIREREREQE0TTd3+OcfaNcOypaF774DL69Up48dg3vuMUF06FAzaHqtQRQgPDImJYiek5CYRHhkzLVfVEREREREchSNjGZTGZ4Ge/w43HcfnDxptm8pUiTV6f37TRCNjjbTcjt3vv7a9sQlXNVxERERERHJfRRGs6EMT4NNToYnnoD162H2bKhePdV1duyAO++E2FjTy+jee11TX4C/D7HpBM8Afx/XPIGIiIiIiGR7mqabDWV4Guybb8KMGWYBaOvWqU6tWwe33AL//QcLFrguiAKEhgTh4+WZ6piPlyehIUGuexIREREREcnWNDKaDWVoGuw338Bbb0HXrmZf0QssXmxm7hYsCAsXQo0arq3v3OisuumKiIiIiMilKIxmQ1ecBrtmjdkotGlTGDMGLCvlPitWmEHSSpUgMtL0NMoM7eoGKnyKiIiIiMglaZpuNnTZabCxsWbYs1QpmD491Uah+/fDgw9CmTLwyy+ZF0RFRERERESuxGVh1LKs0pZlTbQs64BlWScty9psWdZtrrq+nNeubiBDOtQi0N8HCwj092FIh1q0CypstnA5dgxmzoTixVMek5QEnTqZNaLTpkHRoo6VLyIiIiIi4pppupZl+QPLgWXAPcABoBKw3xXXl4tdNA3WtuGRR2DtWvjhB6hZM9X9Bw40jYo++wzq1HFvrSIiIiIiImm5as3oq8Be27afuODYDhddWzLirbdg6lQYNgzatEl1au5cGDwYnnrKLCUVERERERFxmqum6bYDVlmWNdWyrP2WZf1uWVZPy7qgc45knunT4Y03zJ6ioaGpTv3zDzz2GNx0E4we7VB9IiIiIiIiabgqjFYCnge2AyHAKGAo0CO9O1uW1d2yrCjLsqIOHDjgohJyqehoePxxaNwYxo5N1Tn31Cl44AE4c8asE82f38E6RURERERELuCqaboeQJRt22Fnv4+2LKsqJoxeNB5n2/Y4YBxAcHCw7aIacp99+6BtWyhWDL7/HvLlS3W6Tx9YvdoMnFat6lCNIiIiIiIi6XDVyOheYHOaY1uAci66vqRl29Ctm2mP++OPULJkqtNTpphpua+8Ah06OFSjiIiIiIjIJbhqZHQ5EJTmWDXgbxddX9KaMgVmz4aRIy9qj7tli8mpTZrAkCHOlCciIiIiInI5rhoZfQ9oZFlWf8uyqliW9QDwIvCRi64vFzp4EF58EW6+2Xy9QHw83H8/FChgmut6eTlUo4iIiIiIyGW4ZGTUtu01lmW1A94BXgf+Oft1jCuuL2n07g1HjphNQz09Uw7bNjzzDMTEwLx5EBh46UuIiIiIiIg4yVXTdLFtezYw21XXk0uYOxe++sps5VKzZqpTn3wCX39tthxt1cqh+kRERERERDLAsm1nm9kGBwfbUVFRjtaQbRw7BjfeCH5+8Ntv4O2dcmrNGrj1Vrj9dpg5EzxcNQFbRERERETkOliWtda27eC0x102MipuEBYGu3fD8uWpguh//5n9REuXhkmTFERFRERERCTrUxjNLpYtgzFjTMOixo1TDicnwxNPwJ49JqMWKeJgjSIiIiIiIhmkMJodnDxp9mopV84sCL3A0KEwZw589BE0aOBQfSIiIiIiIldJYTQ7eOst0yI3MhJ8fVMOL1kCr78Ojz4Kzz3nYH0iIiIiIiJXSasLs7p162DYMHjySbjzzpTDp0+bbVwqVoSxY8GyHKxRRERERETkKmlkNCs7cwa6djULQUeOTHVq1CgzWDp7dqrBUsmgiOhYwiNj2BOXQIC/D6EhQbSrq41ZRURERETcRWE0K3vvPVi7Fr79NlVnothYGDgQ2rSBu+92sL5sKiI6lrAZG0hITAIgNi6BsBkbABRIRURERETcRGHUxVw24rZ1K7zxBtx3H3TsmOpUaKgZNH3/fdfUnNuER8akBNFzEhKTCI+MURgVEREREXEThVEXctmIm23D00+bvUTHjEm1IHTxYvjmG5NTK1VyZfW5x564hKs6LiIiIiIirqcGRi50uRG3q/LppyZ1hodDQEDK4cREeOEFqFAB+vW7/npzqwB/n6s6LiIiIiIirqcw6kIuGXGLjYU+faB5c7O36AXGjIGNG81SUh/lpmsWGhKEj5dnqmM+Xp6EhgQ5VJGIiIiISO6jMOpC1z3iZtvQo4fZt2X8+FTTc/ftM1NzW7c2y0jl2rWrG8iQDrUI9PfBAgL9fRjSoZbWi4qIiIiIuJHWjLpQaEhQqjWjcJUjbtOmwQ8/wLvvQpUqqU716wcJCWZLF+0pev3a1Q1U+BQRERERcZDCqAudCzfX1E33v/+gZ0+oXx9eeinVqRUrYOJEE0irVcuMykVERERERNxLYdTFrnnE7dVX4dAhmDcP8pz/a0lKMjN3y5SB115zYaEiIiIiIiIOUhjNCtatgwkT4OWXoXbtVKfGjoXff4epU6FAAWfKExERERERcTU1MMoC/n22F0fz+VL7TAOaDP2ZiOhYAA4cgP79oWVLeOABh4sUERERERFxIY2MOmz5mK9p8usSBrfsxpF8vhyJSyBsxgYA5nwcSHw8fPihmhaJiIiIiEjOojDqpKQkSgx+nX8KlWRS3XtSDickJvHGp3vY+GkgL78MNWo4WKOIiIiIiEgm0DRdJ02aRNV923n3tic5nccr5bBtQ8z3VSlZ0uwtKiIiIiIiktNoZNQpCQnw2mtsLhPErBuapjoVv74sp/f6Ez4JChZ0qD4REREREZFMpDDqlPffh9hY/hs/DZ+/85CQmARAUoIXR5bcQPU6p+jUydvZGkVERERERDKJpuk64cABGDIE2ralabf7GdKhFoH+PljAmdU1sE958c0X3mpaJCIiIiIiOZZGRp0waBCcOAHDhgHQrm4g7eoGEh0NweHQo8dF242KiIiIiIjkKBoZdbe//oJPPoGnn4Ybbkg5nJwMPXtC0aImq4qIiIiIiORkGhl1t7Aw8PaGAQNSHZ48GVasgM8+A39/RyoTERERERFxG42MutOKFTB9Orz6KpQsmXL41Cno3x9uvhk6d3auPBEREREREXfRyKi72Db06QOlS8Mrr6Q69emnsHs3TJgAHvp4QEREREREcgGFUXeZMQNWroTx46FAgZTDJ0/CO+/ArbdCq1YO1iciIiIiIuJGCqPukJgI/frBjTfCU0+lOjVuHOzZY9aMaisXERERERHJLRRG3WHsWNi6FWbPBk/PlMMJCWa70dtugxYtHKxPRERERETEzRRGM9uRIzBwoEmbd92V6tQnn8C+fTBlikO1iYiIiIiIOETtcjLbsGFw8CCEh6eah3v8OAwdCi1bmpFRERERERGR3ERhNDPt2gXvvQedOkH9+qlOffwx7N9vBk1FRERERERyG4XRzPTGG2ZLl7ffTnU4Ph7efRfuuMN00RUREREREcltFEYzy7p1MHEivPgilC+f6tRHH8GBAxoVFRERERGR3EthNLO8+ioULgz/+1+qw8eOmeWjrVtD48YO1SYiIiIiIuIwddPNDPPmmdvIkeDvn+rUhx/Cf/9pVFRERERERHI3y7ZtRwsIDg62o6KiHK3BpZKToV49OHoUtmwBb++UU0ePQoUKcMstMGuWcyWKiIiIiIi4i2VZa23bDk57XCOjrjZ9ulkv+tVXqYIowAcfwOHDGhUVERERERHRyKgrJSXBTTeZDrobNoCnZ8qpuDioWBGaNYMffnCuRBEREREREXfSyKg7fPcdbN4MU6emCqIAo0aZQDpggCOViYiIiIiIZCkaGXWVpCSoWRPy5DHTdD3ONyo+fNisFW3VCmbMcK5EERERERERd9PIaGabMgX++AOmTUsVRAHee880L9KoqIiIiIiIiKF9Rl3hzBkYNMisF23fPtWpQ4fg/fehY0dzWkRERERERDQy6hpffw1//mnm4KYZFR0xAuLj4c03HapNREREREQkC9LI6PU6Nypaty60a5fq1MGDZjuXBx80y0lFRERERETEcEkYtSxrgGVZdprbPldcO8ubNAm2bTMLQi0r1anhw+H4cXjjDWdKExERERERyapcOU03Bmh+wfdJLrx21pSYCIMHQ/360KZNqlP798Po0fDII1CjhkP1iYiIiIiIZFGuDKNnbNvOHaOh50ycCDt2wIcfXjQqGh4OCQkaFRUREREREUmPK9eMVrIsK9ayrB2WZU2xLKuSC6+d9Zw+DW+9BTffDHffnerUv//CRx9Bp04QFORQfSIiIiIiIlmYq8LoKqAzcBfwNFAKWGFZVtH07mxZVnfLsqIsy4o6cOCAi0pwsy++gL//hoEDLxoVHTbMZNXXX3emNBERERERkazOsm3b9Re1LF9gOzDUtu2Rl7tvcHCwHRUV5fIaMtWpU1C1KgQGwooVqcLo/v1Qvjw89JDJqyIiIiIiIrmZZVlrbdsOTns8U/YZtW073rKsTUDVzLi+4z7/HHbtgs8+u2hU9MMPTVYNC3OoNhERERERkWwgU/YZtSwrH3ADsDczru+okyfh7behSRO4/fZUp+LjzVrRdu20VlRERERERORyXDIyalnWcGAm8A9QAngdKABMdMX1s5RPP4XYWPjyy4tGRcePh8OHoW9fh2oTERERERHJJlw1TbcM8A1QDDgA/Ao0sm37bxddP2tISIB33oFmzaBFi1SnTp+GkSPhttugYUOH6hMREREREckmXBJGbdt+2BXXyfLGjYO9e+Hrry8aFZ0yBXbvNncRERERERGRy8uUbrpXI9t00z1xAipVgho14OefU51KToabbgIPD1i37qKcKiIiIiIikmu5tZtujvTJJ/Dvv/DddxedmjMHNm2CSZMUREVERERERDIiU7rp5jjHj8PQoaZ7btOmF50eNgzKlTN7i4qIiIiIiMiVaWQ0I8aMgQMHYODAi06tWAHLlsGoUeDl5UBtIiIiIiIi2ZBGRq8kPh7efRdCQuCWWy46PWwYFC0KXbs6UJuIiIiIiEg2pTB6JaNHw8GD6Y6Kbt4MP/4IPXtCgQIO1CYiIiIiIpJNKYxeztGjEB4Od92V7uah4eHg42PCqIiIiIiIiGScwujlrF9v9m1JZ1R092746ivo1g2KFXOgNhERERERkWxMDYwu59ZbITYW8ue/6NR775mc+vLLDtQlIiIiIiKSzWlk9ErSCaKHD8O4cfDww1ChgvtLEhERERERye4URq/Bxx+bJruvvup0JSIiIiIiItmTwuhVSkgwe4q2bg033eR0NSIiIiIiItmTwuhVmjgR9u+Hvn2drkRERERERCT7Uhi9CklJMHw43Hwz3Hab09WIiIiIiIhkX+qmexkR0bGER8awJy6BAH8fbrHqsG1bEd59FyzL6epERERERESyL4XRS4iIjiVsxgYSEpMA2H04gQ+/9CSgfCL33eflcHUiIiIiIiLZm6bpXkJ4ZExKEAU4+XdRTu0rhE/9rXh6OliYiIiIiIhIDqAwegl74hJSfX/018p4+p7kTMWdzhQkIiIiIiKSgyiMXkKAv0/Kn0/tK8jJv4vjF7yDwGLeDlYlIiIiIiKSMyiMXkJoSBA+XmY+7tFVlbHyJlI8OJbQkCCHKxMREREREcn+1MDoEtrVDQRg8Nf/8HdMaQKb/cO7j1ZPOS4iIiIiIiLXTiOjl9GubiANjzcmr5fFmm/KK4iKiIiIiIi4iMLoZezfDxMmwBNPQOnSTlcjIiIiIiKScyiMXsauXVCpEoSGOl2JiIiIiIhIzqI1o5dRvz5s3AiW5XQlIiIiIiIiOYtGRq9AQVRERERERMT1FEZFRERERETE7RRGRURERERExO0URkVERERERMTtFEZFRERERETE7RRGRURERERExO0URkVERERERMTtFEZFRERERETE7RRGRURERERExO0URkVERERERMTtFEZFRERERETE7RRGRURERERExO0URkVERERERMTtFEZFRERERETE7RRGRURERERExO0URkVERERERMTtFEZFRERERETE7RRGRURERERExO0s27adLcCyDgB/O1rElRUDDjpdhOR6eh1KVqDXoWQVei1KVqDXoWQF2eF1WN627eJpDzoeRrMDy7KibNsOdroOyd30OpSsQK9DySr0WpSsQK9DyQqy8+tQ03RFRERERETE7RRGRURERERExO0URjNmnNMFiKDXoWQNeh1KVqHXomQFeh1KVpBtX4daMyoiIiIiIiJup5FRERERERERcTuFUREREREREXE7hdHLsCzrecuydliWddKyrLWWZTV1uibJPSzLGmBZlp3mts/puiTnsyyrmWVZP1qWFXv2ddc5zXnr7Otzj2VZCZZlLbYs60aHypUcKgOvwy/SeY/81aFyJYeyLCvMsqw1lmUdtSzrgGVZMy3LqpnmPnpPlEyVwddhtnxPVBi9BMuyHgJGAe8AdYEVwFzLsso5WpjkNjFA6QtutZwtR3IJX2Aj0AtISOf8q8ArwAtAA2A/MN+yLD+3VSi5wZVehwALSP0eebd7SpNcpDkwBrgFaAmcARZYllXkgvvoPVEyW3Ou/DqEbPieqAZGl2BZ1ipgvW3bT19w7C9gmm3bYc5VJrmFZVkDgI62bde80n1FMotlWfFAT9u2vzj7vQXsAUbbtv322WM+mF+++ti2PdapWiXnSvs6PHvsC6CYbdv3OlWX5D6WZfkCR4B2tm3P1HuiOCHt6/DssS/Ihu+JGhlNh2VZeYH6wLw0p+ZhPpEQcZdKZ6eo7bAsa4plWZWcLkhyvYpAKS54f7RtOwFYit4fxf1utSxrv2VZf1qWNd6yrBJOFyQ5nh/m9+fDZ7/Xe6I4Ie3r8Jxs956oMJq+YoAn8G+a4/9i3nBE3GEV0Bm4C3ga89pbYVlWUSeLklzv3Hug3h/FaT8BTwCtMFMkbwZ+tizL29GqJKcbBfwOrDz7vd4TxQlpX4eQTd8T8zhdQBaXdg6zlc4xkUxh2/bcC78/uwh9O/AkMNKRokTO0/ujOMq27SkXfLvBsqy1wN/APcAMZ6qSnMyyrJHArcCttm0npTmt90Rxi0u9DrPre6JGRtN3EEji4k+0SnDxJ18ibmHbdjywCajqdC2Sq53r6Kz3R8lSbNveA+xG75GSCSzLeg94BGhp2/b2C07pPVHc5jKvw4tkl/dEhdF02LZ9GlgL3JHm1B2YrroibmdZVj7gBmCv07VIrrYD88tXyvvj2ddmU/T+KA6yLKsYEIjeI8XFLMsaBTyKCQB/pDmt90Rxiyu8DtO7f7Z4T9Q03UsbCUyyLGs1sBx4FggAPnG0Ksk1LMsaDswE/sF8wvo6UACY6GRdkvOd7dJX5ey3HkA5y7LqAIds2/7Hsqz3gf6WZf0B/Am8BsQDXztQruRQl3sdnr0NAKZjftGqAAzBdDD93s2lSg5mWdZHwONAO+CwZVnnRkDjbduOt23b1nuiZLYrvQ7Pvl8OIBu+J2prl8uwLOt5zN5RpTF7nb1k2/ZSZ6uS3MKyrClAM0xDrQPAr8Drtm1vdrQwyfEsy2oOLErn1ETbtjuf3crgTeAZoDCm2VYP27Y3uq1IyfEu9zoEngMiMPuA+2N++VqEeY/c5ZYCJVewLOtSvygPtG17wNn76D1RMtWVXodntxOKIBu+JyqMioiIiIiIiNtpzaiIiIiIiIi4ncKoiIiIiIiIuJ3CqIiIiIiIiLidwqiIiIiIiIi4ncKoiIiIiIiIuJ3CqIiIiIiIiLidwqiIiIiIiIi4ncKoiIiIiIiIuJ3CqIiIiIiIiLjd/wFwH0ZAOKZ+UAAAAABJRU5ErkJggg==\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": {},
   "outputs": [],
   "source": [
    "from statsmodels.formula.api import ols\n",
    "\n",
    "data = {\"x1\" : x1, \"y\" : y}\n",
    "\n",
    "res = ols(\"y ~ x1 + np.sin(x1) + I((x1-5)**2)\", data=data).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the `I` to indicate use of the Identity transform. Ie., we do not want any expansion magic from using `**2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept           5.207157\n",
       "x1                  0.469631\n",
       "np.sin(x1)          0.456257\n",
       "I((x1 - 5) ** 2)   -0.017592\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.062988\n",
       "1    10.947653\n",
       "2    10.730061\n",
       "3    10.450988\n",
       "4    10.164108\n",
       "5     9.922852\n",
       "6     9.767331\n",
       "7     9.714506\n",
       "8     9.754046\n",
       "9     9.850850\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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
