{
 "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.982\n",
      "Model:                            OLS   Adj. R-squared:                  0.981\n",
      "Method:                 Least Squares   F-statistic:                     847.1\n",
      "Date:                Tue, 08 Sep 2020   Prob (F-statistic):           3.05e-40\n",
      "Time:                        17:35:08   Log-Likelihood:                -1.5314\n",
      "No. Observations:                  50   AIC:                             11.06\n",
      "Df Residuals:                      46   BIC:                             18.71\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.9557      0.089     55.897      0.000       4.777       5.134\n",
      "x1             0.5119      0.014     37.438      0.000       0.484       0.539\n",
      "x2             0.5263      0.054      9.790      0.000       0.418       0.634\n",
      "x3            -0.0212      0.001    -17.678      0.000      -0.024      -0.019\n",
      "==============================================================================\n",
      "Omnibus:                        0.100   Durbin-Watson:                   2.074\n",
      "Prob(Omnibus):                  0.951   Jarque-Bera (JB):                0.207\n",
      "Skew:                           0.098   Prob(JB):                        0.902\n",
      "Kurtosis:                       2.754   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.42518256  4.92609653  5.38562042  5.77507388  6.07612714  6.28381257\n",
      "  6.40734081  6.46858748  6.49849914  6.53200877  6.60229647  6.73533798\n",
      "  6.94563701  7.23384226  7.5866409   7.97894578  8.37801721  8.74884198\n",
      "  9.05988595  9.28827498  9.42355316  9.46940059  9.44302885  9.3723531\n",
      "  9.29140576  9.23474503  9.23177713  9.30192489  9.45143667  9.67235974\n",
      "  9.94384694 10.23558174 10.51275923 10.74180443 10.89588855 10.95933693\n",
      " 10.93020432 10.82059445 10.65467147 10.46468982 10.28569493 10.14976533\n",
      " 10.08074229 10.09031269 10.17608898 10.32200139 10.50093818 10.67919999\n",
      " 10.82203659 10.89935622]\n"
     ]
    }
   ],
   "source": [
    "ypred = olsres.predict(X)\n",
    "print(ypred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a new sample of explanatory variables Xnew, predict and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10.87581224 10.71322024 10.43221789 10.07983582  9.71798289  9.40828874\n",
      "  9.19701469  9.10372738  9.11650817  9.19487121]\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": "iVBORw0KGgoAAAANSUhEUgAAA6MAAAHWCAYAAACCFl9HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABvy0lEQVR4nO3deZyN5f/H8dc9YzDWyZZMUVEiKkWbRKkmSkn7vmvRXoo2RaFo+7Wrvu17SUUoSdImhYRE0TLWFmsTY+b+/XEjZHfm3LO8no/HeZi5zzn3+cx0OnPe57quzxWEYYgkSZIkScmUEncBkiRJkqSSxzAqSZIkSUo6w6gkSZIkKekMo5IkSZKkpDOMSpIkSZKSzjAqSZIkSUq6UnEXUK1atXDHHXeMuwxJkiRJUgH4+uuvfw/DsPrax2MPozvuuCNjxoyJuwxJkiRJUgEIguDndR13mq4kSZIkKekMo5IkSZKkpDOMSpIkSZKSzjAqSZIkSUo6w6gkSZIkKeli76a7MQsXLmTu3Lnk5ubGXYqKkbS0NGrUqEGlSpXiLkWSJEkqkQp1GF24cCFz5swhMzOT9PR0giCIuyQVA2EYkpOTQ3Z2NoCBVJIkSYpBoZ6mO3fuXDIzMylXrpxBVAkTBAHlypUjMzOTuXPnxl2OJEmSVCIV6jCam5tLenp63GWomEpPT3f6tyRJkhSTQh1GAUdEVWB8bkmSJEnxKfRhVJIkSZJU/BhGJUmSJElJZxgtAOeccw5BEBAEwaotRA455BAefvjhzVqjOGLECIIg4Pfffy/AaiVJkiQp+QyjBeSwww5j1qxZzJgxg/fff5927drRrVs3WrRowZIlS+IuT5IkSZJiVSLC6ICx2TTvPZydugyiee/hDBibXeCPWaZMGWrWrElmZiZ77bUX11xzDSNGjOCbb77h7rvvBuCFF16gWbNmVKxYkRo1anDiiSeu2vtyxowZHHLIIQBUr16dIAg455xzABgyZAgtWrRgm222oUqVKmRlZTF58uQC/5kkSZIkKVGKfRgdMDabrv0nkD0/hxDInp9D1/4TkhJI19aoUSOOPPJI3nzzTQCWLVvG7bffzvjx4xk4cCC///47p556KgA77LDDqttNnDiRWbNm8cADDwCwZMkSrrrqKkaPHs2IESOoXLky7dq1Y9myZUn/mSRJSoQ4PjiWJMWrVNwFFLQ+Q6eQk5u3xrGc3Dz6DJ1C+yaZSa+nYcOGDBs2DIDzzjtv1fGdd96ZRx99lAYNGvDbb7+x/fbbU6VKFQBq1KhBtWrVVt32+OOPX+OcTz/9NJUqVWL06NEcdNBBSfgpJElKnJUfHK/8e73yg2Mglr/VkqTkKPYjozPn52zW8YIWhuGq/S2/+eYbjj32WOrUqUPFihVp2rQpAL/88ssGz/Hjjz9y2mmnUbduXSpVqsS2225Lfn7+Ru8nSVJhtKEPjiVJxVexD6O1MtI363hBmzRpEjvvvDNLliwhKyuLcuXK8fzzz/PVV18xZMgQgI1Ot23Xrh3z5s3j8ccf58svv2Ts2LGUKlXKabqSpCKpsH1wLElKjmIfRjtn1Sc9LXWNY+lpqXTOqp/0Wr777juGDBnCCSecwPfff8/vv/9Oz549Ofjgg9ltt92YO3fuGrcvXbo0AHl5/35a/McffzB58mRuvPFGDjvsMBo0aMCiRYtYvnx5Un8WSZISpbB9cCxJSo5iH0bbN8mkV4fGZGakEwCZGen06tC4wNegLF26lNmzZzNz5kzGjx/PvffeS6tWrdhnn3247rrrqF27NmXKlOGhhx7ip59+YtCgQdxyyy1rnKNOnToEQcCgQYOYN28eixcvZptttqFatWo88cQTTJs2jY8//piLL76YUqWK/fJfSVIxVZg+OJYkJU+JSDDtm2QmvQHCsGHD2G677UhNTSUjI4NGjRrRrVs3LrroIkqXLk358uV59tlnufHGG3n44YfZY489uPfeeznyyCNXnSMzM5Pbb7+dm266iQsuuICzzjqLZ555hldffZUrrriCRo0aUa9ePe65557/NDWSJKmoWPk3us/QKcycn0OtjHQ6Z9W3eZEkFXNBGIaxFtC0adNwzJgx67xu8uTJNGjQIMkVqSTxOSZJkiQVrCAIvg7DsOnax4v9NF1JkiRJUuFjGJUkSZIkJZ1hVJIkSZKUdIZRSZIkSVLSGUYlSZIkSUlnGJUkSZIkJZ1hVJIkSZKUdIZRSZIkSVLSGUYlSZIkSUlnGJUkSZIkJZ1hNMGCINjg5Zxzzom7REmSJEmKXam4CyhuZs2aterrgQMHcuGFF65xLD09fY3b5+bmkpaWlrT6JEmSJKkwcGQ0wWrWrLnqkpGRscaxf/75h4yMDF5++WUOPfRQ0tPTefzxx3nmmWeoUKHCGucZMWIEQRDw+++/rzr22Wef0bJlS8qVK0dmZiaXXHIJCxcuTOaPJ0mSJEkJYRiNQdeuXbn00kuZNGkS7du336T7TJgwgSOOOIJjjjmG8ePH079/f8aNG8d5551XsMVKkiRJUgHYpGm6QRAcDFwH7APUAs4Nw/CZ1a7vAFwE7A1UAw4Jw3BEoosFuOoqGDeuIM68fnvtBfffn7jzXX755ZxwwgmbdZ8+ffpw8sknc+2116469uijj9KkSRPmzp1LjRo1ElegJElbaMDYbPoMncLM+TnUykinc1Z92jfJjLssSVIhtKlrRisA3wHPrbisrTzwGfDCeq7Xapo2bbrZ9/n666+ZNm0ar7766qpjYRgC8OOPPxpGJUmxGzA2m679J5CTmwdA9vwcuvafAGAglST9xyaF0TAM3wPeAwiC4Jl1XP/8iuuqJbK4dUnkCGVcypcvv8b3KSkpq4LlSrm5uWt8n5+fzwUXXMDVV1/9n/NlZvoHXpIUvz5Dp6wKoivl5ObRZ+gUw6gk6T/splsIVK9enb///puFCxdSqVIlAMatNRd57733ZuLEidSrVy+GCiVJJdHmTrmdOT9ns45Lkko2GxgVAvvttx/ly5ena9euTJs2jTfffJNHHnlkjdvccMMNjB49mosvvpixY8cybdo0Bg4cyEUXXRRT1ZKk4mzllNvs+TmE/DvldsDY7PXep1ZG+mYdlySVbLGE0SAIOgZBMCYIgjHz5s2Lo4RCpUqVKrz44ot88MEHNG7cmH79+tGjR481brPHHnswcuRIZsyYQcuWLdlzzz3p2rUr2267bUxVS5KKsw1NuV2fzln1SU9LXeNYeloqnbPqF0iNkqSiLVh7reJG7xAEi4HLVu+mu9p11YB5bEY33aZNm4ZjxoxZ53WTJ0+mQYMGm1WftDl8jknSuu3UZRDreocQANN7H/XfK8IQJk7k22ffZMGg99l+9gz+yqhOtd13oXaTBlC7NtSpE1122AHSHS1V0ZSTA4MHw+uvw0EHQadOcVckFX5BEHwdhuF/uri6ZlSSJP1HrYx0stex1nONKbfTp8OHH8Lw4dFlzhz2AKhbFw7Zn53mzIFxX8B7/SE/f80T1agRBdPVQ+pOO8FhhxW6oOp2Nfr7738D6MCBsGQJVKsGe+8dd2VS0bap+4xWAFZ2zkkBagdBsBfwZxiGvwRBUAWoDWSsuE29IAjmA7PDMJyd0IolSVKB65xVf41tWgC2X7qQ+8IZcOHLUQidPj26omZNaN0aWrcmPLQ1Pyytw8SJUKVKlDlrbJNLlZxsUn79GX75BX7+Obr88gt89x2891403ATRO/xOneDSS6M7x8ztakqu9QXQM86AE0+Eli2hlMM60lbZpGm6QRC0Aj5ax1XPhmF4ThAE5wBPr+P628MwvG1D53aaruLkc0yS1m/A2Gze+d+7HPT5YFr+NoG6c1aEz8qVoVWrVQH0r5oNGP5RwNCh8P77Uc5cW0pK9Ea+Ro11XKqHZJb5nZ0WjKP+Bw+SMvBdKFMGzjoLrrkGdtstqT/36pr3Hr7OEeLMjHQ+7XJoDBWpIK0MoK+9BoMGRQG0enXo0GG1AJqSD198Ae++G32IUq3ampfq1aN/q1SBtLS4fySpUNiqabor1n8GG7j+GeCZLaxNkiQVNosW0f7pu2j/8ENRMDzoILiqI7RuzfI99uarb1J5/30YegF8+WU0C7dSpSifdukCzZrBwoUwd+6/l3nz/v3666+jfxcsgOgtRnXgcDIzD+eGS7/nnPn3UfH55+CJJ+Doo+G66+DggyFY79uRAuF2NcVffj68/Ta88ko0Avr331GePOMMOOmk6GlXKn8ZjBgBl/WPbjx7djQsWq5c9ERfn4yMNQPq2qF1771hzz2T9aNKhY6TCyRJ0preeSeaKpudHf175538/FelKHzeFc3QnT8/yoX77gs33QRZWdHXmzsQtHTpvyF16lR45hm48tHduCp4nJMPvYNbqz9C/Q8eImjVCvbZB669Fk44IWkjTpu0dlZF1ujRcNll8NVXUTY888zVAujSJTB0KJzTP0qpCxZE4bNtWzjuODjqqGiWwLJl8Pvv/73Mm7fm97/8At98Ex1ftuzfIk46CXr2jNZaSyXMZnfTTTSn6SpOPsckaTWzZsHll8Obb0KjRuQ+8gT3frY/Tz8NU1bs6LL99lHwzMqKRkGrVEl8GTNmwFNPRZdZs6BurRzubfI8bSbfQ9pPP0RNj668Ei64IBqOLUBrrxmFaLuaXh0au2a0CJszB7p2haefhu22g7vuglNPhVIL/4yCZ//+0ZzznJzoSX7MMVEAPfzwrW+wFYbR/N+5c6NPX+65B3Jz4ZJL4Oabo1QsFTPrm6ZrGFWJ5nNMkojmKT7xBNxwA/zzD3TrxufNr+PCS9OYOBEOPRTatYsC6G67JW+mbG5utG7v8cejAaogzOfWfQbRaek9VPvu4yiIXnhhNIW3Zs0Cq8NuusVHbi489BDcdluUM6+6Cm65fikV33wG3ngDPvoI8vIgMzMKn8cdt2KYtAAnE86aFRX05JNQvnw0z/2qq6JRWKmYMIxK6+BzTFKJN2kSdOwIn34Khx7Kwrsf44Ynd+Gxx6LtQB9+OAqicVt7tPTIamO4a9t7aPz96wSVK8Mjj8DJJ8ddpgqxDz+EK66InvJHHgn33w/1pw5k8SWXUeG3n/mpSiafNj6YHc4/jVant426biXT5MnRcO3bb0OtWtCjB5x9NqSmJrcOqQCsL4wm+f8ySZK00oCx2TTvPZydugyiee/hDBibnbwHXzECyl57weTJhP97mtcvGkb9o3ehX79oYGbSpMIRRAF23DF6b/7zz/DWW0DTpuw16WUa5E1kWsoucMopURj944+4S1UhM2MGHH98tIXtP/9ES6Lfu/8H6l/dFtq1Y05OPmee1J1DL3iMW/Y/nUt+KMWA8bOSX2iDBjBgAIwcGX0SdP75UXOjQYOiqb1SMWQY1Va77LLLaNWq1arvzznnHI4++uitOudtt91Go0aNtrIySSq8Vq5FzJ6fQ8i/+1cmJZCOHBmF0O7d4aST+PX9ybR78xxOOjlgu+2ipi733QcVKhR8KZsrLQ3at4+23/jpJ2h3XX0azx9Fj7J3kvfmW4SNGkVr/lTi5eTA7bdHGW/IELjzTpj4xSLajbqBoHEjGDWK/2t7MVnnPsgnO+29av55Tm4efYZOia/wFi3g88+jacNLl0bdpA89NOqyJBUzhtECkp2dTceOHdl+++0pXbo0mZmZXHjhhfz2229r3G5jwW38+PEce+yx1KxZk7Jly1K7dm2OP/54fl7XJm6FxAMPPMALL7ywSbedMWMGQRCw9lTt6667jo8//rggypOkQqHP0ClrNMWBJLwJ/uuvaI1ly5awdCl5g4Zw794v0KBlDT76KOqjMnp01LS2KNhxR+jTB775thTD9r2RvfO+YtrCGtFw7vnnb3jLDRVbYRj14GrQIFqKeeyx8P2kfG7c4XnK7rEr3H13tG/LDz9wX+OjWZ763/WgG9u6p8BnNQRBNJw7aVK0yHXixKhd9SmnwI8/JvaxpBgZRgvA9OnTadq0Kd999x3PPvss06ZN44UXXmDixIk0a9aMGTNmbNJ55s2bR+vWralQoQKDBg3i+++/5/nnn6du3bosTPAf2GWrtxjfSpUrVyYjI2OrzlGhQgWqVq2amIIkqRBK+v6V338fTfl7+mno3Jmxz3/Hvrdkce21UTadOBGuuaZg+7QUlAYNoi0gr3lmT1qWHU3voCv5Tz9DfqPGMHx43OUpiX78MWp4e8IJUX+rESPglc5fs8OpB8FZZ0WdmL/4Av73P6hZc71b9Gxo656kzmpIS4u2V/rxR7jlFnj33egJf+WV0XYxUhFnGC0AnTp1IiUlhWHDhtG6dWtq167NIYccwrBhw0hJSaFTp06bdJ5PP/2Uv/76i6effpp99tmHHXfckZYtW3L33XfTuHHj9d5v5WjrHXfcwbbbbkuFChU499xzycn59w1Oq1atuOSSS7juuuuoXr06zZs3B2DSpEkcddRRVKxYkRo1anDqqacye/bsVffLy8vjuuuuY5tttmGbbbbhqquuIi8vb52Pv1IYhtxzzz3ssssulClThu23356uXbsCsNNOOwHQrFkzgiBYNd137Wm6+fn59OjRgx122IEyZcrQuHFj3n777VXXrxxhffPNNzn88MMpV64cDRs25IMPPtik37UkJduWvAneYt98w9IDmvPHX4tpd9I91PrsMpq2LEd2Nrz6ajSrdccdE/+wyRQEUa+XCT+UYeq5PTkw/JTps8pG+89ccQX8/XfcJaqAvfBCNPt8zBh48EH4Zug8Wr7YEZo1i8Lc//4XTX/db79V9+mcVZ/0tDUbBKWnpdI5q/56HyeWWQ0VK0bT6qdNg3PPjUZL69aNRnnz8wvucaUCZhhNsD///JMhQ4bQqVMnyq3VkrtcuXJceumlDB48mL/++muj56pZsyb5+fm88cYbbG7X448//pjx48fz4Ycf8uabb/L+++9zww03rHGbF154gTAM+eSTT3juueeYNWsWBx98MI0aNWL06NEMGzaMxYsXc8wxx5C/4oXunnvu4YknnuDxxx/n888/Jy8vjxdffHGDtdx444306NGDrl27MnHiRF5//XV22GEHAEaPHg3AkCFDmDVrFv3791/nOR544AH69OnDXXfdxYQJEzjuuOPo0KED48aNW+N2N910E1dccQXjx4+nWbNmnHLKKSxevHizfneSlAxb8iZ4i4waRW7LVvyRX4qjmvdj8KALmfVpbSrt9Sv3vDaTk05K3lYtyVC1atRxt8/I/Tmp3lge4Ap48EFyGzeJRsRU7CxaFA16nnkmNGkC336znMvy/49SDXaJZgJcfTX88EMU4tbqkNu+SSa9OjQmMyOdAMjMSN/oHrJJn9Wwuu22i/Y6+u67aErDDTdEP9fy5QX/2FIBKHqTca66CtYKIAVur72i/t+bYOrUqYRhuN7tQho2bEgYhkydOpV99913g+faf//9ufHGGzn77LPp1KkTzZo1o1WrVpx++unUqVNng/dNTU3l6aefpkKFCjRq1Ii77rqL888/n169elG+fHkgGpW85557Vt3n1ltvZc899+Suu+5adey5556jSpUqjBkzhn333Zf777+f66+/npNOOgmIQuLQoUPXW8fixYu57777uP/++znvvPMAqFevHgcccAAA1Vds7Fy1alVqbmCPuL59+3Lddddx2mmnAdC9e3dGjhxJ375911ifevXVV9NuRevHnj178txzzzFu3DgOOuigDf6+JCnZVr7ZLdD9K4cOheOOY1b5qhzV4Gm+G9qatKqLqXb6Z5Td/i8e+Syd0w+ulbjHS7IN7f/ZogV8Pr4c9977AG26Hcvj089l+wObw/U3kHJ7NyhTJubqlQhjxsCpp0bNrG67DW46fDSljjkvmnd++OHwwAPRtNYNaN8kc7P+v6uVkU72OoJngcxqWJ8GDaItYO64A269FRYvhpde8nmtIseR0QISrOdj5pUjnOu7fm133nkns2fPpl+/fjRu3JinnnqKhg0b8uGHH27wfnvssQcVVmuDeMABB7Bs2TJ+XG3R+z5rdaj4+uuvGTlyJBUqVFh1WTmC+eOPP7JgwQJmzZq1KkgCpKSksN9q013WNmnSJJYuXUrr1q036eddl4ULFzJz5sxVU4lXOuigg5g0adIax/bYY49VX9eqFb3Bmjt37hY/tiQVpPZNMvm0y6FM730Un3Y5NLFB9M03oV07wvr1ab3zq3z3yWGk15tLzbNHUXb7aHZOUkZyCsimrNsrXRq6dIFHvj+Ua1p/y9PhOaTc1Yu/G+0L48cnpcbYtu4p5vLzoW9fOPDAqOHsiOH5dCvXh1Itm0dDpW+9FX0YUwB7iSdtVsPGBEG0jvT++6F/fzjmGFiyJLk1SFup6I2MbuIIZVx22WUXgiBg4sSJtG/f/j/XT548mSAIqFu37iafs2rVqpx44omceOKJ9OrViyZNmtCjR4+tCnjAqhHSlfLz8znqqKPo27fvf2677bbbrpqquzk2d3rxhqwrwK99LC0t7T/XbUndklSkPfssnHce4f7706XRIKb1y6B8o9+oeuS3BKn/vi4ndSQnwTa0bm/tUL/TTvD6+5V5882nOOei9vSediFp++xH+PCjlL7o3AKpb2VYXlnjyrAMJPZDhxJo9uxoffD770OHDvDknXPY5qqzo/B5wgnwxBOwlY0UNyQpsxo2x5VXRmtKL7wQsrKifUkrV46nFmkzOTKaYFWqVCErK4tHHnmEv9dqlvD333/z8MMP06ZNG6pUqbJF5y9dujR169bd6DrICRMmsGS1T8e++OKLVfddn7333puJEydSp04d6tWrt8alYsWKVK5cme22244vVltzE4bhqnWf69KwYUPKlCmz3pHc0qVLA/ynCdLqKlWqRK1atRg1atQax0eNGkXDhg3Xez9JKpEefBDOOYf8Q1tz/vbvc3e/DNqdvpjtj/1ujSAay0hOAm3uur0giHLK/01vx/9dMIGP8w6i9MXn8cfxHeGffxJeXyxNbkqAIUOiptAjR8Jjj8EbF33ANq32hI8/jtZSvvZagQbRlQp0VsOWOO88eOWVaG+mQw6BefPirUfaRIbRAvDQQw+xfPlyDjvsMIYPH86vv/7KiBEjOPzwwwnDkIceemiN2y9cuJBx48atcZkxYwYDBw7kjDPOYODAgfzwww9MmTKFvn378t5773HcccdtsIbly5dz3nnnMXHiRD744AO6dOnChRde+J/R0NV16tSJBQsWcPLJJ/Pll1/y008/MWzYMDp27MiiRYsAuPLKK7n77rt54403mDJlCldddRWzZs1a7zkrVqzIlVdeSdeuXXn66af58ccfGT16NI8++igANWrUID09naFDhzJnzhwWLFiwzvN07tyZvn378vLLL/PDDz9w66238sknn3Dttddu8PcgSSVGGEbrx664guXHHEf7lHd5+rXy9OwJbz9fgd7Hb16TlsJuS7sRV6oEPZ+oTvjeEB4s34Wq/Z9g9q4tyJ/xS0Lri7XJTTG0bBlcdx20aQM1asCYz3O5aEZXgiOzoq5VX30FHTsWr25cm+vEE6N1pJMnw8EHw1p720uFUdGbplsE1K1blzFjxtC9e3fOPPNM5s6dS/Xq1Wnbti2vvvoq22+//Rq3/+STT2jSpMkax44//njuvvtuKlSowHXXXcevv/5KqVKl2Gmnnejbty9XXnnlBmto2bIlu+++O4cccgh///33qvNtSK1atfj000/p2rUrRx55JP/88w+1a9fmiCOOoMyKBfHXXnsts2fP5oILLgDgzDPP5PTTT2fy5MnrPW+vXr3YZptt6NGjB7/99hvbbrstZ511FgClSpXi//7v/+jevTu33347LVq0YMSIEf85xxVXXMGiRYu4/vrrmTNnDvXr1+fNN99kr7322uDPJEklQhjC9ddD374sPflMjvj1f3zyeSkefzx6fw6b36SlsOucVX+NabCweaO9h7cpRZMZvejZdj86fXU2i3bdm+XPv0LVkw9LSH2FoslNMTF1atSk6Ouv4dJL4Z7LplP2vNOi7sgdO8J998FaOxiUWG3aRNOVjz466uI1bFi0BYxUSAWJXNO3JZo2bRqOGTNmnddNnjx5vV1ptX7nnHMOv//+OwMHDoy7lELP55ikIi8vL3qH3q8fS87pxIFj/o/vf0jhxRejaanF2Ya66W6qMIRXe/xA49s6sFs4mR/O6EGDZ7v8ZwuQLaltXWG5qI9IJ8uAsdncPWQKUz+twp8fNKJc2RSefzaF9rmvR2sjIVobeuKJ8RZaWI0ZE60fLVMGPvgAdt897opUwgVB8HUYhk3XPu7IqCRJRVVubrTB4iuv8NclN7LP4DuYOy9g0CA4LDEDfIVaIkZ7gwBOuXVXvj/qS94/7ELavHAT40d9yS6fPku5WhlbVRsUoiY3RciAsdlc//IkZg5qyJJJmZTZ4Q/qHPUVez7/AvR/EfbbD15+OepMpXVr2jRaWHv44dGU3SFDoFmzuKuS/sMwKklSUZSTE40KDRrEzCvvYu9Xric3F4YPh41sY6112G2f8uw060UGHLU/Rw2/lpk7NiPnxf7sdmLjLT5ncZsanSy3PjmT6S8dyPIF5ah80BSa1RvGQ2/cTZ0/fo326uneHVbrnq/12H13+OST6JOp1q1h4MAomEqFiA2MiqFnnnnGKbqSVJz9/Te0bQvvvce0ax+l4TPXk5YGo0YZRLdGmbIB7T+8gnH3jaBM3hJqn7Qfg05/EXcIS46Ve4dOeGwfwrwUtj31My4t/yDvvnANVXIWcvaJ3aFXL4Po5qhbN3phyMyMpu2+917cFUlrMIxKklSU5OVF3VxGjmTsNc/T+OGLqVkTPv0UXAK/cQPGZtO893B26jKI5r2HM2Bs9n9u0+yq5pSe8A3TqzbjqJfO4O0dr2DmjGUxVFtyzJkTfb7SuTNUafA79U8bzJNfd6Hn0IcZvf3utDn3QX5qcmDcZRZNmZnRlN0GDeDYY+H11+OuSFql0E/TDcOQoCS36VaBibt5lyRttjCEK66Ad95h9FkPceD9p7PXXjB4MFSvHndxhd/aTYWy5+fQtf8EgP9Mp63SsCbbzBzGhKO7cNwH9/LlLl8z/rHXaXN+raTXXRAS0fwpUd5/P1r6vGABPPoo7Jb+BfWuvJRtF/1Or1bn0G/fDpQtncbNRXhf3NhVrw4ffQRHHQWnnAKLFkV7k0oxK9Qjo2lpaeTkuB+XCkZOTg5pTvWRVJT07QuPPMKkozqz33OdaNkyen9pEN00fYZOWaO7LUBObh59hk5Z5+2D0mk0fv8eZt73Ko3zx7P3BU2467AP+OuvZFRbcFaG8uz5OYT8G8rXNUqcqMdb12j0smXRjkRZK7cK/TKfixfeTasLjmeb8qXpdNH99NvvBGptU94uxIlQuXK07cthh8H558P998ddkVS4t3ZZuHAhc+bMITMzk/T0dEdIlRBhGJKTk0N2djbbbrstlSpVirskSdq4V16BU09lVsuTqTPqJZq3SGHwYChbNu7Cio6dugxiXe96AmB676M2eN9l4ycx//ATqTZvMg9XuIG6L3Sn7bFF8wPN5r2Hr3MP1MyMdD7tcmhCH2t9W9xc0Wwvnr6jJl99BRddBPfeMIdyF58VDZOecEK0bUtGRkJr0QpLl8Jpp0H//nD77XDLLVFbaakAFcmtXVaGhJkzZ5KbmxtzNSpO0tLSDKKSio6PP4azz2ZRk4Np9NUzNNg9hQEDDKKbq1ZG+jpDWK2M9I3et/SeDakx4yvmnXEVl7/Vm8/af8x1J7zMzU/UKXKZaeY6fgcbOr411jUa/fv4bbm8TzUqpkfLF0/IGAYHngnz58Njj0HHjoajglSmDLz6KlxwAXTrBunp0WJdKQaFOoxCFEgNDJKkEmvyZGjfnmXb78w+v7xFxeplGTw4mnGnzdM5q/46R+k6b+paxHLlqN6/H7kvtGbvCzrS4I296Dz8KY5/sQNHHllARReArQnlm2v1gJu/LJU/P9idJd/tQJnMPxn3cSXq/K9b1CF3t92iUdHGW76VTmFaB1volSoF//tftEXUDTfALrtA+/ZxV6USqFCvGZUkqUSbNQvatCEvrQyH5Q7mT6owdCjUKh49dJKufZNMenVoTGZGOgHRtNQtWYuYdsbJlJ00lrSGu/DEn8fzY5tOXHzOPyxYUDB1J1rnrPqkp6WucWyzQvlmWBlwl82pxKxnD2LJd9tT+cCptDn1feqc1RJ69owa6Xz11VYH0WSugy0WUlLgmWegWTM4/XT45pu4K1IJVKjXjEqSVGItXgwtWxJOmcIZ23/MgF/3Yfhw2G+/uAvTKsuWsfyGmyh1f1/GswdXbvsqNz63G0ccEXdhG5esUcSXPpnJZZ2X8dfo2qSWX0a1dmM5Mfc97hn6IKXJh8cfj7Yq2krJXAdb7MyeHW1QnJcHo0dHW8FICba+NaOGUUmSCpvly+GYYwjff5+bG7/DXRPa8u670KZN3IVpnd57j9zTzyZ3wd9cGj5M6QvOpu89ASV5lVEYwssvw7XXwpw5ITWazqRis7H0/uYJjv/iHWjaNGrKVbduQh5va5pTCfj2W2jeHHbdNdqTtHz5uCtSMbO+MOo0XUmSCpMwhEsugcGDeWrvR+g5ri1PPWUQLdTatiVt4njKHLwfz3AuLZ88k/13X8SwYXEXFo/vvoNDDolmfm6/PXzxRcDs5xYxddRNURC99lr49NOEBVFY/3rXglgHWyztsUf04cC4cXDGGZCfH3dFKiEMo5IkFSY9e8KTT/JBsxu58KuO9O4NZ58dd1HaqFq1SP3wA+jRg9NSXua9OXtz/eHfcMkl8OefcReXHAsXwjXXwF57wYQJ0QzcLz7KYd/Bt0cHZ86EQYOi/XJLl07oYydzHWyxddRRcM89MGAA3Hhj3NWohDCMSpJUWDz/PNx8MxP3PoMjvrqDK6+E66//780GjM2mee/h7NRlEM17D7dJS2GRmgo330wwYgS1a/zD6JT9KfPYA+y8Yz633gp//RV3gQUjDOHFF6F+fbj/fjj/fPjhB+hYayCpe+wOt90WdWodPx7ati2QGhLVnKrEu/JKuPhiuOsuePrpuKtRCeCaUUmSCoMPP4Qjj2T2Li2oPXkIHU4uzUsvRQ0vV7eya+ja25P4xruQ+eOPqEvsO+8wbZtmnPPXvUyodBBXXglXXw3bbBN3gYnx3XfQqVO0zLBpU3jkEWhW9aco1AwcCA0awEMPwaE2ESoycnOjUdKPPoIPPoBWreKuSMWAa0YlSSqsJkyADh1YlFmfRj/05+DWpXn22f8GUYA+Q6esEUQBcnLz6DN0SpKK1SapWjWa7vjss9QrN5NRtGBoheN5qcc0dtwRunUr2iOlq0/J/e67f6fkNht0GzRsCCNGRNNxx483iBY1aWnw2mtQrx4cfzxMnRp3RSrGDKOSJMUpOxvatmVZmQrsM2cwOzTOoH9/KFNm3TefuY7tKzZ0XDEKAjjrrGjOavfu7L9gKD+UasgLNa7hwe5/stNO0QzW+fPjLnTT5ebCCy+sY0rudu9GU3Jvvx06dIDvv48aFaWlxV2ytkRGRjSyHQRw9NElZ+Gzks4wKklSXJYsgWOOIf/P+RyR+x7Lt9uBwYPZ4JYgdg0tgsqVg1tugalTSTnnbNr99ABzK9Xj3jr30/P2Zey4Y5ThtjSUFvQa4jCEMWPgiiuiLSjPPDPqkvvll/D4DT9R9Zx2cMwxkJ4Ow4fDSy+5V2VxULduNLo/YwaccAIsWxZ3RSqGDKOSJMUhPx/OOYdw7FgurPAyk9L2ZOhQqFlzw3eza2gRtt128MQTMG4cpfZvxnnfXs2iHRpyY/03ue22kJ12gu7dYcGCTT/lyjXE2fNzCIHs+Tl07T8hIYH0l1+gV69o1m2zZtFU3JYt4e231zMld9y4aE8XFR8HHQRPPhmtH7300uiTCSmBDKOSJMXhttvgjTe4P7MPry45mkGDYJddNn43u4YWA40bw9ChMHgwZSqV5frRJ7CoycFcsOdXdOsGO+4IHTtG02F//nnDp0r0GuKFC6MmqoceGtVx441QrVoURGfPhtdfzeeY8O01p+ROmeKU3OLszDPhppvgqaeirV+kBLKbriRJyfbyy3DaaQyrfS5Zvz7F2+8EHH103EUpFsuXw//+F03jnTuXP9ucxs1BT176tM6qEdLataFFCzj44Ojf3XaLlvIB7NRlEOt6JxcA03sftcklfPBBtLPQgAGQkxP1rjnzTDjjDNh5pxC+/Taafvvyy/Drr9GI6MMP22m1pMjPh1NOgTfegP79o616pM2wvm66hlFJkpJp9Gg4+GB+rLYvDbM/oO//leHyy+MuSrFbtAjuuou8vn1ZnpfPoF0P4qsdDuGfXU4je05dRo6EOXOim1avHs2ebNECnvnpK+anzyNIWfP9XGZGOp92ibrYhmE04vnXX/+9TJwY5cs5c6BKFTj55Kjn0n77QTBjehRAX3oJJk2CUqUgKwtOOw1OPNGR0JLm77+jDx8mToRRo6BJk7grUhFiGJUkKW6//QbNmrEwtyx1/xjNqZdX5//+L+6iVFgMGJvNA08Pp+PHL5A19Quq5Cwkn4D5jfZimw7H8GvjNgz7qykjP01l5EiYPj26X0rp5ZTO/JNSlXPI/yeNYFlptkuvBEtL8+efUWOk/Px1P2ZaWtQs9ayzoG1bKD1/brStx0svweefRzdq0SIKoCecEM3ZVck1ezbsu2/0hPrySxtVaZMZRiVJitOSJdCiBcu/n8Y+Sz+jdttGDBgAqakbvadKiOa9h5O9YouelPw8Gs+eRqufvubwX8bS6LfvoyHOqlXhiCOgTRtmNs7i48k1eO6txXw8EpYuLE3pcsvZYbtS7JxZmm22gW22iUY8V3699rFq1SB9+aJofu5LL0XzdfPyYI89ogB66qnRPGFppfHjo6H5XXeFkSOhfPm4K1IRYBiVJCku+flw0kmEb73FiWXe5cf6bfnkE6hQIe7CVJhscP1n5/2joDh4MAwZAnPnRlfusw+0aQNHHhk1RsrJiaZTbuyyZEn077Rp0X6SOTlRx6KVAbRRoyT+5CpyBg6EY4+N1o6+/jqk2BNVG7a+MFoqjmIkSSpRunWDN9+kR+V7+KJCW74caBDVf9XKSF81Mrr2capWjRrInHJK9OHGuHH/BtNeveCOOzb/AdPSogWo550Hp58O++//b2ckaUOOPjrqrHv11XDzzdCzZ9wVqYjapJHRIAgOBq4D9gFqAeeGYfjMatcHQDegI7AN8CXQKQzDiRs7tyOjkqRi7aWX4PTTGVD1fM7IeYJRnwbstVfcRakwWrln6OpbtaSnpW58656//oJhw6J9YMqVi6ZNliu3/kv58pCebgMibZ0whIsvhn794MUXo1F1aT22dmS0AvAd8NyKy9quB64FzgGmALcCHwRBUD8Mw0VbVLEkSUXdl18SnnceE6sezCl/PsIb7xhEtX4rA2efoVOYOT+HWhnpdM6qv/E9ZLfZJupuKyVTEMCDD8L338P550cbJTdrFndVKmI2e81oEASLgctWjoyuGBWdCTwUhuGdK46lA3OB68IwfHxD53NkVJJULP36K+y7L3/8nU79haPp9n/V3MJFUvEzb17UYXfpUhgzBmrVirsiFULrGxlNxGrjnYCawPsrD4RhmAOMBA5MwPklSSpaliyBY45h2fwlHLzwXU6/wiCq4mHA2Gya9x7OTl0G0bz3cAaMzY67JMWtenV4++1oM9v27aNmWNImSkQYrbni3zlrHZ+z2nVrCIKgYxAEY4IgGDNv3rwElCBJUiGRnw9nnkn47bcct/RV6rbbnXvvjbsoaeutXNOaPT+HEMien0PX/hMMpIq2AnrhBfjqK7jggmg9qbQJEtmHee1nXbCOY9ENw7BfGIZNwzBsWr169QSWIElSzG69Fd56i66l+jJrrza89JJ7iap46DN0yhrNlQBycvPoM3RKTBWpUGnfPurq/NJLcNddcVejIiIRYXT2in/XHgWtwX9HSyVJKr5efBHuvJOXyl/AC9Wu4t133cJFxcfMdWw7s6HjKoFuvDHafujGG+Hdd+OuRkVAIsLodKJAevjKA0EQlAVaAJ8l4PySJBV+H35IeO65jKnQkk7hwwwcFJC5kSaoUlFSKyN9s46rBAoCeOop2HvvaKuXiRvd5VEl3CaF0SAIKgRBsFcQBHutuE/tFd/XDqN2vPcDXYIg6BAEQSPgGWAx8FLBlC1JUiEybhzhccfxc9n6HLFkAC+8VtotXFTsdM6qT3ramnPO09NS6ZxVP6aKVCiVKxc1NKpQAY45Bv74I+6KVIht6shoU2Dsiks6cPuKr7uvuP5u4F7gYWAMsB1whHuMSpKKvenTCdu04c+8DA5aNJiej2Rw1FFxFyUlXvsmmfTq0JjMjHQCIDMjnV4dGm98H1SVPJmZMGAAZGfDCSdAbm7cFamQ2ux9RhPNfUYlSUXW778TNm/O3z/Po9nSUZzWoyE33xx3UZJUSDz/PJx1FlxyCTzySNzVKEbr22e0VBzFSJJU5C1ZAkcfzfKffuGI5cOoe2pt3ksdzlNdcqiVkU7nrPqOGEkq2c48EyZMgD59oHHjKJRKqzGMSpK0uZYvh5NPJn/0V5wYvkla2yZM3Xkk/yyItr1Yuf8iYCCVVLL16gWTJsHll8Nuu8Ehh8RdkQqRRO4zKklS8ReGcNFFMGgQncKHWX5Ue5Ye+CX/LHf/RUn6j9TUaO/RXXeN1o/++GPcFakQMYxKkrQ5unWD//2PO1NuYULzi3ntNZi96O913tT9FyUJqFQJ3nkn+jDvmGNg4cK4K1IhYRiVJGlTPfYY9OjBM6XO59WGt/Puu9EuBu6/KEkbUa8evPEGTJkCp58OeXkbv4+KPcOoJEmbYsAAwk6dGJp2NHdu/xhD3w/YZpvoKvdflKRNcOih8MADMHAgth4X2MBIkqSNGzWK/FNOZWypfbk441WGDSvFdtv9e/XKJkV9hk5h5ny76UrSel16adRht3dvaNQoGiVVieU+o5IkbcikSeQf2Jyf/65B67Kf8tYn1dhzz7iLkqQibNkyOOII+OILGDkS9t037opUwNa3z6jTdCVJWp/ffiM/60j+/LssRwZDeXaQQVSStlrp0tH60e22g6OPjrZ+UYlkGJUkaV3mzyf/yDbkzJpPVt5g+r6xIy1axF2UJBUT1arB0KHR1i+HHQbTpsVdkWJgGJUkaW0LFxIe3Y68SVNolzeAK5/ei3bt4i5KkoqZXXeFYcOiabutW8PPP8ddkZLMMCpJ0upmzyZs2ZK8z77gtPBF2t17KGedFXdRklRM7b47fPABLFgQBdLs7LgrUhIZRiVJxdqAsdk07z2cnboMonnv4QwYu4E3Oj/+SP6BzVk6YSpHhQOpf9OJXH118mqVpBKpSRMYMgTmzImm7M6dG3dFShLDqCSp2BowNpuu/SeQPT+HEMien0PX/hPWHUjHjiX/gANZ+OsCDs4bTtY9WdxxR9JLlqSSaf/9YdCgaKru4YfDn3/GXZGSwDAqSSq2+gydQk5u3hrHcnLz6DN0ypo3HD6c/INbMnt+WQ7iU655eV+uuSaJhUqS4OCD4e234fvvISsrmrqrYs0wKklJtFlTRrXVZs7P2fjx118n/8g2TPmnDoeW+YwH36/PKackqUBJ0poOPzza9mXcOGjbFhYvjrsiFSDDqCQlyWZNGVVC1MpI3/DxRx4hPPlkvsjflxOqj+T1zzI55JAkFihJ+q927eDll+GLL+CYYyBn3R8squgzjEpSkmzylFElTOes+qSnpa5xLD0tlc5H7Aq33gqdOjGQdly+6/sM+XIbGjeOqVBJ0ppOOAGefRZGjIAOHWDp0rgrUgEoFXcBklRUDRibTZ+hU5g5P4daGel0zqpP+yaZ6739Jk0ZVUKt/O+xxn+nw+px7GM94Il+PMn5vHTQYwx7uxTbbBNzsZKkNZ1xRjQq2rEjnHIKvPYapKXFXZUSyDAqSVtg5ZTblSOdK6fcAusNpLUy0sleR/Bc31RSJUb7Jpn//jf55x/CU08jGPAWd3AT357Qg/eeDyhbNt4aJUnrceGFUSC98ko46yx44QVITd34/VQkOE1XkrbAlky5Xe+U0az6BVKj1jJ/PvmHZxEOGMDl/B9/XHUHr7xqEJWkQu+KK6B3b3jlFbjgAsjPj7siJYgjo5K0BTZrym0YwvTptJ8+jt3mjGTup1+xZDksqlaThvs2pNH3f8Hi7WH77aFWLShTpoCrL4FmzWL54UcSTprMmbzEvvec4tYtklSU3HAD/P03dO8O6enw8MMQBHFXpa1kGJWkLbC+Kbd1KqTC2LFRS/px46Kvx4+HhQsB2C0lhd123TWal/Ltt/BZ//+evHr1KJhuvz1kZv779e67w957Q4qTWjbLxx+z/MxzWPbbPI5PGcTZLxzu1i2SVBTddls0ZbdPnyiQ9u1rIC3iDKOStAU6Z9Wn65vf0nDGd+w5ayoN5/7E7vOmU/+PX+GW3OhG5crBnnvC6adDkyaw117QqFH0B3SlhQvht9/+vWRn//v1L7/AZ5/BH3/8e/vq1aFNGzjqKDjiCMjISOaPXbTMnQudO8NzzzEzdUfOKjeCbu82desWSSqqggDuuisaIb333ujvbI8ecVelrWAYlaQt0H7Zb7QY1I2q474C4PeKVchrvAcp5570b/CsW3eNJgthGGXNb76JlrtUqwbVqlWiWo2GbFO/4fr7MeTkRHccPRoGDYKBA+G556JzN28ebQretm0UdIvIJ8Sb24l4s+Tnw5NPknd9F/IXLqYPXXm97s0890Y5t26RpKIuCOD//g/++QfuuANKl4abby4yf/+0piAMw1gLaNq0aThmzJhYa5CkTTZ9OnTtCq++CttuC7ffDsceCzVr/uems2fDmDFrXubMWfdpgwCqVFkZUNd9qVsX9t0XypTKgy+/hPfei8LpuHHRSXbY4d9geuihUKFCwf0etsLanYghauTUq0PjrQ+k48eTf9HFpHz5BZ+ktOSKtEc5+bYGXHNN9H5FklRM5OXB2WfDiy9Gf4efeCKaPaRCKQiCr8MwbPqf44ZRSdoEf/0FPXtGn8ampsJ110VTQCtWBOD33+Hrr9cMnr/9Ft01JQUaNoSmTaPLPvtEPYp+/z26zJv379fruuTm/ltG2bJw4IHQqlV02XdfKPN7NgwZEgXTDz6AxYuj5NWqVTSdt0OHaM1pIdG89/B1rrfNzEjn0y6HbtlJFy2Cbt0I/+//+JMqXJ3Xl4XHnMkD/xdQp85WFixJKpzy8+GBB6BLl+gT3WeegaysuKvSOhhGJWlLLFsGjz4ade/76y8455xofUpmJt98E/VO+PxzmDHj37vUr/9v8GzaNJqxu6WDlGEY5ax58+C772DEiOgyfnx0XXr6WuF0r2WU/vKTaNT0vffg+++jEx14IJx4IpxwQuzBdKcug1jXX54AmN77qM07WRjCm2+Sd8VVBLNm0o+OPLZDT3o8XIV27RJRrSSp0Bs/PurPMHFitB9p7964b1fhYhiVpM0RhvDWW1Er+WnToHXrKHnutRe//AI33RTtu12lChx22L/Bc++9oXLlgi/vzz9h5Mg1wylE4bR5c6ix6wK+zZtGtdxvOOW3Lzhl+udUnjo5utEBB8BJJ8UWTBM2MvrTT4SdLiMYMphvU/fiUh6l5Q37c9NNUU8LSVIJkpMT/c1+8MGoh8JLL2GjgMLDMCpJm2r0aLj2Whg1Kppf27cvHHkkCxYG9OoF998frfG86qpoZlAywufG/PHHv+H0ncG5zJiaBkCQtpyydX6n8q6/88Bx8zlpzsfw+uv/ptcDDvh3xHSHHZJS61avGV26FPr0If+OO8nJLcVN+T2Y2OoyHny0FLvtVoCFS5IKv8GD4dxzo9lMd90FV1zhlmiFgGFUkjZmxoyoOdErr0TNibp3h/POIzcsxWOPRd/+/juceWbUwK927bgLXrfmvYfzy6zlLP21CjkzqpPzY3XyFkZDhXvsEfU3Or7xDzT56Q1S33z93wZI++//bzAt4B9ui7rp/vQTDB5M3gMPkjp1Cq9zAj2r38/1D2Ryyik2UpQkrTB3LlxwAbz7Lhx+eLSWtFatuKsq0QyjkrQh/frB5ZdHzYmuvRauv56wQkUGDIhm/UydGjWo7dMnmopbmK29JjMMIfePCvzzYw32pAGjRsHy5bDNNlGfh5P3nsrh81+n/HurBdNddoGDDorm/B50EOy6a/LT3tKl8Em0/jV87z2CKVMA+CG1AVfl30u9y46kR4/CMTItSSpkwhAefxyuuSZau/Hkk9C+fdxVlViGUUlan7vuiubbHnlk1Bp+++358ssok376aTRTt08faNOmaIy+bWxN5oIFUdPdlT2O5syJfq5994XT95tGewaQOf0TUj77NJr/C9HeMiuDafPmUUvggtgr5ZdfoilWgwcTDhtGsGQJuallGJXakgHL2jIkaEvt1rtw112F/0MBSVIh8P33cNppMHYsXHgh3HcflC8fd1UljmFUktYWhnDjjVHXvVNPhWef5adf0+jaFV57LZqp26NHtPSkVKm4i910m7MmMz8/+vs8aFAUTEePjn4tZcvCHo1D2uw8hUNLj6Lh/E+pMmkUKT9Oi+5YtmyUXlcG1AMOiIZaN1duLnz22b/J+LvvAJhTtg4DlrXl3fy2fFPpEA5uU5527aIPBKpU2eJfjSSpJFq2DG65JfpkeZddor1Jm/4nF6kAGUYlaXX5+XDZZdG2LRdfzMKeD3Fbj1QeegjS0v7dRnRLt2SJ2xatySTaQmbYsGjP1G++iS4LFkTXlSoFLevP5rhtP6N5OIp6s0dRfupYguXLoxuUKxdNc155KVVqze/XvpQqRfjTTwQLF7I8JY2vy7XgtcVteY+2LK+7G+2OCWjXLsq6aWkF+MuSJJUMH30EZ50Fs2dHjSCuvz76e6QCZxiVpJVyc6PhzhdfhOuvZ/ZVvTmyTcCECdHh7t3tc7BSGML06f8G02++iYLq779H11cIlnBC7dG02eZzqgZ/Ql4eQV4e5OVBfh4pecujY/nR90F+dH0QRl//tGRb+ue04aOgNXu2qES7dnD00dFerUVhSrQkqYj56y+46KKos/zBB8Nzz0GdOnFXVewZRiUJ4J9/4OST4Z13oFcvfjyxC0ccEX1I2r9/1NBHGxaGkJ29ZkAdOxaWLPl3QHTloOjq/67r6+23j8Kn028lSUkThlEIveyy6JPP44+PmhsdfnjR2aj6l1+iqUwffhh1D/7gg7gr2iDDqCQtWgTHHhtN03n4YcYdeClHHhkNlL73Huy3X9wFSpKkpPnpJ7jttmgLmPnzIT09+lS6ffvok9KqVWMucDV//BG9f1kZQKet6OFQowYcdhg8+2yhbnCxvjBaeCuWpET6889o+O3rr+H55xlZ+wzatYRKlaLX9gYN4i5QkiQl1c47RyOkubkwciS89RYMGBBdUlOhRQs47rjog+xkT+VdsgRGjfo3fI4bF43oVqwILVtCp07QujU0alSk17U4Miqp+Js1C444An74AV57jbc5lpNPhp12gqFDoXbtuAuUJEmFQhhGH1yvDKUTJ0bHmzSJRkzbt4fGjRMfAHNz4auvouA5bBh8/nl0LC0NDjwwCp6tW0OzZkWyq5/TdCWVTDNmRNNXZs+Gt9/m6V9ac8EFUUf3QYOi7TO3tPOsJEkq5qZOhbffjoLpZ59FYXWnnf4Nps2br9mRNwwhJyca2VzX5e+/1/x+8eJo1HPEiOjrIIiC78rwedBBxWJfVMOopBJh9WB5wLK5PPXijaTnLoXBg7n74/244YaoP0H//tG2LZuzJ6ckSSrBZs+O1pcOGBCNXi5bFnXfq1x5zbC5ufmqXr3og/PWreGQQwrXWtUEMYxKKvZWD5aNZk/judduJS8llW+fepUhYw6nb9+oke5zz0Hp0tF9mvceTvb8nP+cKzMjnU+7HJrknyB+jhJH/D1IkjZo0SIYMgQGD45CafnyW3ZJTy/Saz43VYE3MAqCoCLQAzgOqAGMBa4Mw/CrRD2GJG1In6FTyMnNo+lvE3n69dtYULYip510J9/f3Yh530Rr/R94YM3ZNDPXEUQ3dDwOyQpGa48SZ8/PoWv/CQCFIoj5e5AkFRoVK8KJJ0YXbbGUBJ7rSSALOBtoDLwPDAuCwL/ckpJi5vwctl8whyfevIO5FarQ4eS+jBnRjnnfbMdtt8GDD64ZRAFqZaSv81zrO55sK4NR9vwcQv4NRgPGZif8sVaG+dXl5ObRZ+iUhD/W5vL3IElS8ZOQMBoEQTpwPNAlDMMRYRhOC8PwNmAacEkiHkOSNmbH8ik8+lZPUsN8zmnXg28HH0XOtBrsdMwUunVb9yyYzln1SU9bM6Gmp6XSOat+kqresGQGo8I8SuzvQZKk4idRI6OlgFTgn7WO5wAHJegxJGn9wpBnxjxL4zk/ckXrG/jyveNZOjODWh3Gc+9tFdZ7t/ZNMunVoTGZGekERGtFC1PzomQGo8I8SuzvQZKk4icha0bDMFwUBMHnwM1BEHwHzAZOBQ4gGh2VpIL15JPUeedVvjv7Sl4fcjXL55el4TnjufOyGhsNlu2bZBaa8Lm2Whnp62ywVBDBqHNW/XV2Fi4Mo8T+HiRJKn4SuWb0TCAf+A1YClwBvAzkrX3DIAg6BkEwJgiCMfPmzUtgCZJKpNGj4bLLCI/I4vaF9/DPvAq8/14pJj7VpNCGzE2VzGnEhXmU2N+DJEnFT8K3dgmCoDxQKQzDWUEQvApUCMPwqPXd3q1dJG2VefNgn30gJYW+p35N595VueceuOaauAtLHLcZifh7kCSpaEr6PqNBEGwDTAeuD8Ow3/puZxiVtMXy8iArC0aN4uOen9Lq2n0444xoH9ESsGVXoWFIlCRJG5KMfUaziKb9fg/UA/oAU4CnE/UYkrSGW26BDz9kZo+naHfbPjRtCv36GUSTyT05JUnSlkrkmtHKwENEYfQ5YBRwRBiGuQl8DEmKDBgAvXqx9OwLafXceaSnQ//+kG7D06RyT05JkrSlEjYyGobha8BriTqfJK3XDz/AWWcRNm3GibMeZPp0GD4cdtgh7sJKHvfklCRJWyqRI6OSVPAWL4YOHaB0ae7e9w3efb8MDz4ILVrEXVjJ5J6ckiRpSxlGJSXFgLHZNO89nJ26DKJ57+EMGJu9+ScJQ7jgApg8mREXv0KXR2rTsSNcfHHi69WmSeaWK5IkqXhJ2DRdSVqfhDW5eeABePVVZl3ek6PuO4wDD4QHHyyIirWpVv73s5uuJEnaXAW2tcumcmsXqfhr3ns42etYQ5iZkc6nXQ7dtJN88gkccghLs9qx28T+5C4PGDMGatZMcLGSJElKqALf2kWS1merm9zMmgUnnUS4886cuPgZZs0O+OQTg6gkSVJR5ppRSQVuq5rcLFsGJ54ICxdy1779eXdkZR5/HJo1S3CRkiRJSirDqKQCt6VNbgaMzea1Q06BTz/lqsa30/XFRlx5JZx9dkFWK0mSpGRwmq6kArclTW4GjM3m3b7P8tRnb/HkbifwwJirKbfjH7Q4/R/A5jiSJElFnQ2MJBVKh90+iGfuPY+lKaXZM3c8/5Qqy3ZnjaJ2rVKb3vRIkiRJsbOBkaQi5cx3HqPWwnm0rjaIv5dUouZJn5FaLpeZ83PjLk2SJEkJ4JpRSYXPJ59w9jcD6VfrDEb83oYqWRMoXWMRsIlNjyRJklToGUYlFS5//w3nnceCGrW5dtbDlKs/k/INZwKb1vRIkiRJRYPTdCUVLt26wbRpdN55GGnLy7HrieOYl8smNT2SJElS0WEYlVR4fPkl3HsvXzXpyBNjW/P223DMMS3irkqSJEkFwGm6kgqHpUvhvPNYWq0WWePv5uyz4Zhj4i5KkiRJBcWRUUmFwx13wKRJXLH9IMrXqsz998ddkCRJkgqSYVRS/MaOhV69+KrhWfSb1JahQyEjI+6iJEmSVJCcpispXrm5cN55LKtcjSMn3ccll8ARR8RdlCRJkgqaI6OS4nX33TBuHFdu25+Mnatw991xFyRJkqRkMIxKis/EidC9O1/XPYnHfzqOESOgQoW4i5IkSVIyOE1XUjzy8uD881lWtiJtfnyQq6+Ggw+OuyhJkiQliyOjkuLxwAPw5ZdcU+Ulqu5WgzvuiLsgSZIkJZNhVFLyTZ0KN93ENzscw2PZp/DZYEhPj7soSZIkJZPTdCUlV34+XHABuSllOPrXR+l6Y8C++8ZdlCRJkpLNkVFJyfXYYzByJNdV/B819qzFLbfEXZAkSZLiYBiVlDwzZhBefz3jtz2CR/84hzHPQenScRclSZKkOBhGJSVHGELHjizPCzh2Tj9u7xmwxx5xFyVJkqS4GEYlJccLL8AHH9A1/WG2268OnTvHXZAkSZLiZBiVVPD++ovwuuuYkrEfj+RczNhnoZSvPpIkSSWabwclFbybbyac9zunhkPoeV8K9evHXZAkSZLiZhiVVLDGjCF89FGeLHM5ZfdqwuWXx12QJEmSCgPDqKSCk5cHl1zC/LI16bKsBx/3g9TUuIuSJElSYWAYlVRw+vWDMWO4hJe55MZKNG4cd0GSJEkqLAyjkgrGnDmEXbvyWXprvq51Mk/fHHdBkiRJKkwMo5IKRufO5C36m/PyH+axxwPS0+MuSJIkSYVJStwFSCqGPv4Ynn+eu8PrOeDs+rRuHXdBkiRJKmwcGZWUWMuWEV56KbPK7MhjFW7km75xFyRJkqTCyDAqKbHuv59g0iQ68i49nyxHtWpxFyRJkqTCyDAqKXF++YX8227nvVLHsuyQozn99LgLkiRJUmFlGJW0WQaMzabP0CnMnJ9DrYx0OmfVp32TzOjKq65i2TK4ttQDvPcoBEG8tUqSJKnwMoxK2mQDxmbTtf8EcnLzAMien0PX/hMAaD9zHLz1Ft3ozXl31qFu3RgLlSRJUqFnN11Jm6zP0CmrguhKObl5PDDwW/I7Xc4PpRowrNHVXHNNTAVKkiSpyHBkVNImmzk/Z53HjxvyHCk/T+ciPuLRp0qTlpbkwiRJklTkODIqaZPVykj/z7Gd/szm4i/e5HnOYI8rWrHvvjEUJkmSpCInIWE0CILUIAh6BEEwPQiCf1b8e0cQBI68SsVI56z6pKel/nsgDLnjg8fIoRz3bdeHO+6IrzZJkiQVLYkKizcAnYCzgQnAHsCzwFKgR4IeQ1LMVnbNXdlN98xfv6T5jLF04iFue6wmFSvGXKAkSZKKjESF0QOBd8MwfHfF9zOCIHgH2C9B55dUSLRvkhmF0oULWb7LhXwT7M284y7mmGPirkySJElFSaLWjI4CDgmCYDeAIAgaAocC7yXo/JIKmfDWbqTMnc215R7l/gdTN34HSZIkaTWJGhm9C6gITAqCIG/Fee8Mw/CRdd04CIKOQEeA2rVrJ6gESUkzfjzh//0fj3MRJ/fdl1q14i5IkiRJRU2iRkZPBs4CTgP2XvH1pUEQnL+uG4dh2C8Mw6ZhGDatXr16gkqQlBT5+eSefxF/UJW39+1Jx45xFyRJkqSiKFEjo32AvmEYvrLi+wlBENQBugJPJegxJBUG/fqR9vWXdE59nnv+tw0pbhAlSZKkLZCoMFoOyFvrWB7uYyoVL7Nnk3tdF0ZyKHVuPJ3dd4+7IEmSJBVViQqj7wJdgiCYDkwEmgDXAM8l6PySCoHcy68hf0kO9+z8CG/dFMRdjiRJkoqwRIXRy4n2E30EqAHMAp4Auifo/JLi9sEHpL3xMrfTjZufr0+ZMnEXJEmSpKIsIWE0DMNFwFUrLpKKm3/+Iee8S/mVXZh/cRcOPDDugiRJklTUJWpkVFIxtrx7T9J/m0a36sPod3fZuMuRJElSMWCDIUkb9v33BHf15gVO58xnWlOxYtwFSZIkqTgwjEpavzBkyVmXsDC/PJ8edw9t28ZdkCRJkooLp+lKWq/8Z5+n/FcjuKX843R/fNu4y5EkSVIxYhiVtG5//MHSy65lLAew9yMXUL163AVJkiSpOHGarqR1WtTpBtKW/MULzR/j9DN9qZAkSVJi+Q5T0n+En4yi4qtP8WDaNdzw4h4EQdwVSZIkqbhxmq6kNS1bxoJTL2Y+dShzZzfq1Im7IEmSJBVHjoxKWsPiHveSkT2Rh+o/xEXXlI+7HEmSJBVThlFJ/5o+nbRe3Xkr6MC5bx5NamrcBUmSJKm4MoxKioQhc0/sxLK8VKZf9QC77x53QZIkSSrODKOSAMh54U1qfD2Yh7ftQade28ddjiRJkoo5GxhJgoULWXrJlUymCS1fv4wyZeIuSJIkScWdI6OSmHX+zVRaMovhJz3OAS38jEqSJEkFzzAqlXDLPhtDjTce5vkKl3LxU83iLkeSJEklhEMgUkm2dCl/dLgAqEGtp++kQoW4C5IkSVJJ4cioVIL9duHtbDdnPK8c0o/DT6gcdzmSJEkqQQyjUgm1+P3P2O75u3i1wvmcP6Bd3OVIkiSphHGarlQSLV7MkuPPYi512PGt+6hUKe6CJEmSVNI4MiqVQFOPvY7qi3/i43OfZb/DKsZdjiRJkkogw6hUwsz+33vsMvxxXs68jjP7tYi7HEmSJJVQhlGpBFk+5w9KXXw+E1Ma0fzDHpRyor4kSZJi4ltRqZgYMDabPkOnMHN+DrUy0umcVZ/2TTL/vUEYMqX1JeyS+wejew6hbf0y8RUrSZKkEs8wKhUDA8Zm07X/BHJy8wDInp9D1/4TAFYF0qndX2b3ia/z6l69OLnrnqvut8EAK0mSJBUQp+lKxUCfoVNWBdGVcnLz6DN0CgCLJv9Gje6dGFPmQI78sDPwb4DNnp9DyL8BdsDY7GSXL0mSpBLIMCoVAzPn56z/eH4+Px96Lqn5uQTPPUflKqnAxgOsJEmSVJAMo1IxUCsjfb3Hx5z3CI1mD+Ojo+9ln5PqrrpugwFWkiRJKmCGUakY6JxVn/S01DWOpaelcsU2pWn47PV8ltGWNv0vXOP6DQVYSZIkqaAZRqVioH2TTHp1aExmRjoBkJmRzp1tG9D4umv4J0hn+yFPUiotWOM+6wuwnbPqJ7FySZIklVR205WKifZNMtfohDvi0O40XPwVn1z+Gi32226dtwfspitJkqRYBGEYxlpA06ZNwzFjxsRag1TcfPu/MTQ8f3++3PEUmk9/Ie5yJEmSVIIFQfB1GIZN1z7uNF2pmFk4J4dyF5/JvNSaNBrxYNzlSJIkSetkGJWKmc9bdaVe7vf80fcZKtfZJu5yJEmSpHUyjErFyPCbh5P1/QN8ue/lNLrqsLjLkSRJktbLMCoVE5M/ms2uPc/m5/T67PNB77jLkSRJkjbIMCoVA9nTcvgn61iq8Cdl3nyZUpXKxV2SJEmStEGGUamIW7Qgn++ansOeuV8x554XqdmmSdwlSZIkSRtlGJWKsNxceGef28la8BpTz+vNTle3j7skSZIkaZMYRqUiKgzh2awXOf3H7kw58FzqP9k57pIkSZKkTWYYlYqoFy79jDM/Oo+farek/kePQRDEXZIkSZK0yQyjUhH0zv/NIOux9vxVoTY7jnkTSpeOuyRJkiRps5SKuwBJm+fTwQupe9XRpJfKpcxnA0mpXjXukiRJkqTN5sioVIR8/91yco45mV3DKYSvvUHpxvXjLkmSJEnaIgkJo0EQzAiCIFzHZVAizi8J5syBL5pfw2HLh7Cg5yNUOq513CVJkiRJWyxR03SbAamrfb8d8DXwWoLOL5Vof/8Nz+33MJ0XPsjs066hZtcL4y5JkiRJ2ioJCaNhGM5b/fsgCM4HFgKvJ+L8UkmWlwd3HzaUm3++ktlNj6bmc3fHXZIkSZK01RK+ZjQIggA4H3ghDMO/E31+qaS56+xJXP35ScyvtTs1h78Eqakbv5MkSZJUyBVEA6PDgZ2AJ9d3gyAIOgZBMCYIgjHz5s1b382kEu/xO+ZxyotHE5RLp9rn70LFinGXJEmSJCVEQYTRC4GvwjAct74bhGHYLwzDpmEYNq1evXoBlCAVfQNeXcrutxxHZsosKnz4DtSuHXdJkiRJUsIkdJ/RIAhqAMcCnRJ5XqmoGjA2mz5DpzBzfg61MtLpnFWf9k0yN3q/d98JWXLahbTnU5Y9+yop+++bhGolSZKk5En0yOg5wFLglQSfVypyBozNpmv/CWTPzyEEsufn0LX/BAaMzV7vfcIQHrh7KfOPPYvT859nyQ3dKX3GSckrWpIkSUqShIXRFY2LLgBeCcNwUaLOKxVVfYZOISc3b41jObl59Bk6ZZ23z82Fq8/+k71uOIIzeYFlt91J+V43J6NUSZIkKekSOU23FbALcEYCzykVWTPn52zy8b/+gsvb/sjNXxxFvdTp5D/7EqVPP7WgS5QkSZJik7AwGobhR0CQqPNJRV2tjHSy1xE8a2Wkr/H91Klwc+vPeejXY6hUPp9Sg4dBixbJKlOSJEmKRUF005UEdM6qT3ramnuCpqel0jmr/qrvR4yAO5u8wTO/HkqF7TMo883nBlFJkiSVCIZRqYC0b5JJrw6NycxIJwAyM9Lp1aHxqm66/3sqZHDrPjyz5ESCvfcmfeznsOuu8RYtSZIkJUlCt3aRtKb2TTL/s5VLXh7cdMNydrznMu7icZZ1OJmyLz4DZcuuus2WbgkjSZIkFRWGUSmJFi+GC05exNnvnUQbhpB/fRdK97oTUv6dpLByS5iVnXhXbgkDGEglSZJUbDhNV0qS336D4/f7ja7vHURWygfQrx8pd/VaI4jC5m8JI0mSJBVFjoxKSTBmDHRtM45n/ziKGumLSHlrEGRlrfO2m7MljCRJklRUOTIqFaClS+GBB+DO5u/x1h8tqF4jhVJfjFpvEIX/bv2yseOSJElSUWQY1UYNGJtN897D2anLIJr3Hs6AsdlF+nGSIT8fXn4Z9q2/gL+v6soby9pRptEupH3zJeyxxwbvuylbwkiSJElFndN0tUHJaqazNY+zJZ1nC7Jb7bBhcFPnZew37jFGpHZnG/4gPPNMUh95BCpU2Oj9V9ZhN11JkiQVZ0EYhrEW0LRp03DMmDGx1qD1a957ONnrWKuYmZHOp10Ojf1x1g6xEI0irr6fZyLusynGjYMbrg+p/MHr9CnVlTrLfyI85BCCPn1gn322+LySJElSURYEwddhGDZd+7jTdLVByWqms6WPsyWdZxPdrXbGDDjjDLiiyUh6Dt+f1ziZHeqXg/feI/jwQ4OoJEmStA6GUW1QsprpbOnjbEmITVTA/v13uPpqOHbXSZzy8jGMpCVNts2G//2PlPHjoE0bCILNOqckSZJUUhhGtUHJaqazpY+zJSF2awP2339Dr15w4E6zaHh/R8Yub0zb8h9Dr16kTP0Bzj0XUlM3fiJJkiSpBDOMaoPaN8mkV4fGZGakExCt4dzatZWJfJwtCbFbcp8whClToE8f2KvuInJvvJVvc+pxfqlnSLniclJ++hG6dIFy5Tb+w0qSJEmygZGKvoLqprtsGYwcCQMHwqCBIRV+HMfRDOTqtIeokjsXTj4Z7rwT6tYtyB9PkiRJKtLW18DIMCqtZu5ceO+9KIB+MXQB+y0eRruU92iXNpiqS2dFNzr00Gie7r77xlusJEmSVASsL4y6z6hKtDCE8eOj8Dnw3ZDFoyfRhve4pvR77Jc7ilSWE1asTJCVBW3bwpFHwrbbxl22JEmSVOQZRlVi5ObCTz/B999H6z8nT4ZP319C/ZnDact7vFX6PbbjFwDC3fYgaHsdtG1LcMABUKpUNLX36YnMnD9mk6cDS5IkSVo3w6gKxJas40yEMIy2XJkyJbqsDJ7ZkxaQOuNH6uT9SF1+pB7TuDBtKv3yviCNZeSXr0DK4YdB25uhTRuC7bf/z8/Ttf+EVfuTZs/PoWv/CQAGUkmSJGkLGEaVcMkIbn/9BdOmwY8/Rv9Omxoy97u5hFOnUX3Rv4HzhOBHdkn5kSp5v69x//zqNUipVxcOuAzatiXloIOgTJn1Pl6foVNW/Twr5eTm0WfoFMOoJEmStAUMo0q4RAS3MITZs/8Nmz9OC5n33RyW/vAzKb/MoOqSn6nDz+zIDNrzMzsFMygfLvn3/ikpLN9uB0rVr0dQr0PU8bZuXahXD3bemZSKFTfrZ5o5P2ezjkuSJEnaMMOoEm5LgtvcuTD8w5Dv3ppK3hdfUWb2DGrlRoHzAGZwMr+Qzj9r3GdZhW3I22FH0urtQqmdD4uC5orQGey4I2mlSyfsZ6qVkU72OuqvlZGesMeQJEmSShLDqBJuU4LbwoXw8YiQ8f1/ZPmwEeyS/RGtGMEpzFx1m78r1mDZdnVI3XkPSjc4BnauAzvuCHXqQJ06lK5UKRk/DgCds+qvMfUYID0tlc5Z9ZNWgyRJklScGEaVcOsKbmWCNLIyGnFPp59YOuQjak8fQctwBO34DYDFFbZl6QGHkNe+FakHN4edd6ZcuXKUi+uHWMvK6cVxNGWSJEmSiiPDqBJuZUDr/sKvLB+1hBbTJ7H/H2NomX8BtfkVgMXla7CkWStyOxxC2uGtqFC/PhWCIM6yN6p9k0zDpyRJkpQghlEl3JTvQ6ZdPZZ+H99GU74GYHF6NRbu3Yqc47qQ3qYVFRo0KPThU5IkSVLBMYwqYX6YEjLg0qG0Gn4r1/EVf2TszIKr76Nyh8OosPvuhk9JkiRJqxhGtdWmTQ1589IPaTHsVq7nc/6sVIeFtz1J1cvOgrS0uMuTJEmSVAgZRrXFfvwRXu80guZDb+UGPmF+he1ZeOtjVLnyXEjgtiqSJEmSih/DqDbbTz/Bq5d9wv6Du9GFj5hfvhYLb3qIjGsugDJl4i5PkiRJUhFgGNUmmzEDXrr8c5oN6kbX8AMWltuWBV3uJ+O6jpCevtH7S5IkSdJKhlFt1D//wANnfc0eb9zCjeFgFqVXZ0Hne6h8w8VQrrDsBCpJkiSpKDGMaoNmz4anD+xH5+mX8k/ZDOZfcxcZN3aC8uXjLk2SJElSEWYY1XqNHZPH6ENuoOvie5jdpA01P3oZKleOuyxJkiRJxUBK3AWocHr7pSX8uv8JXLT4Huad1Imao98xiEqSJElKGMOo1hCGcP/1M8k8vSVH5b3Dwh4PUP3Vh6CUg+iSJEmSEseEUcIMGJtNn6FTmDk/h1oZ6XTOqk/7JpkA5OTA7R3G02nI0VQv9Rd5r75NpQ5Hx1yxJEmSpOLIMFqCDBibTdf+E8jJzQMge34OXftPAGC/mpn0aTWI2384hbBSZcqMGEXQZK8Yq5UkSZJUnDlNtwTpM3TKqiC6Uk5uHt3+l80jDR6kzw/HkLfzrlSaPNogKkmSJKlAOTJagsycn/OfY/9MrsZ17z7HZeHDLGh1DBkDX3LbFkmSJEkFzpHREqRWRvqqr8MQln9ck6fe6c5l4cMsuegaKg/rbxCVJEmSlBSG0RKkc1Z90tNSyc9NoUz/Ggz84nyyGMrXnXtS/rF7IDU17hIlSZIklRBO0y1B2jfJZNGCgH4nT+XVucdSMXUhX97/LM0vOz3u0iRJkiSVMIbREiQ3Fz6/9leGzG1LUL0a5YZ/QfNGjeIuS5IkSVIJlLBpukEQbBcEwbNBEMwLguCfIAgmBUHQMlHn19YJQ7jp3Jnc/M1x5FfflnLffgkGUUmSJEkxScjIaBAEGcCnwCjgKGAesDMwNxHn19Z75L6lHPfi8VRLW0DpD4dCzZpxlyRJkiSpBEvUNN3rgVlhGJ612rHpCTq3ttKQwSFlr+3EAXxB/guvQ+PGcZckSZIkqYRL1DTd9sCXQRC8GgTB3CAIxgVBcFkQBEGCzq8tNGkSfNDhEc7nKZZ1vomUk06IuyRJkiRJSlgY3Rm4FPgJyAIeAHoDndZ14yAIOgZBMCYIgjHz5s1LUAla2++/Q4/DPqb3P1eR0/poSvfuHndJkiRJkgRAEIbh1p8kCJYBY8IwPHC1Yz2B48IwbLCh+zZt2jQcM2bMVtegNS1bBme0+JmHRzelQu2qpH/7JVSuHHdZkiRJkkqYIAi+DsOw6drHEzUyOguYtNaxyUDtBJ1fmyEM4coL/6bL6OOonL6M9PffNohKkiRJKlQS1cDoU6D+Wsd2BX5O0Pm1Ge67N+Tg585nL8aR8sZAqL/2fxpJkiRJileiRkbvA/YPguCmIAjqBUFwInAF8HCCzq9NNHAgzLmuD6fyCtx5J7RtG3dJkiRJkvQfCRkZDcPwqyAI2gM9gVuAX1b8+0gizq9N89138NSJQ3iTLizvcBKlunaJuyRJkiRJWqdETdMlDMNBwKBEnU+bZ+5cuPzIqQxYegrLG+5B6ef+B+6sI0mSJKmQSlgYVXyWLoXTj1nEIzOPpVylUqQNGgDly8ddliRJkiStl2G0iAtDuLhjPp2+PJP6KT+Q8tYHsOOOcZclSZIkSRtkGC3i+vaFOs91pz1vw30PwCGHxF2SJEmSJG2UYbQIGzwYPr/+LfpzO+E55xBcfnncJUmSJEnSJjGMFlELFsBdZ09iUMpZ5O+9LymPPmrDIkmSJElFRqL2GVWS3XBdHr3nnUfpyumkDOgPZcvGXZIkSZIkbTJHRougjz6ClCcfZ3++hAdfgMzMuEuSJEmSpM1iGC1i/v4bbjp3JkNTupLX8jBSTzst7pIkSZIkabM5TbeI6dYNrvr5KsqnLiW1n+tEJUmSJBVNhtEi5Kuv4Pt7BnESr5PS7RaoVy/ukiRJkiRpizhNt4hYtgwuO3cJb6R0Iq9eA1I7d467JEmSJEnaYobRIqJ3bzh+4u3swM/wxEgoXTrukiRJkiRpixlGi4CJE+GdHuP5MrgXzr8AWrSIuyRJkiRJ2iquGS3k8vLgwvPyeIyLCKpWgbvuirskSZIkSdpqjowWcg8+CHuNfpymfAn3vwBVqsRdkiRJkiRtNcNoIfbTT/Bg15l8W6orYavDCNxTVJIkSVIx4TTdQioMoWNH6LP8KsqlLiV41D1FJUmSJBUfjowWUk8/DWU+HEQHXofb7nBPUUmSJEnFimG0EJo5E26+eglfl+lEuHMDAvcUlSRJklTMGEYLmTCETp3guiW3s13ez/C4e4pKkiRJKn4Mo4XMG2/ATwPG82bKvXCBe4pKkiRJKp4Mo4XIH3/AFZ3yeL/cRQTl3VNUkiRJUvFlGC1ErrkGTvjjcRrnfwn93FNUkiRJUvFlGC0khgyBD56byU9lukKLw2Aje4oOGJtNn6FTmDk/h1oZ6XTOqk/7JplJqlaSJEmSto5htBBYtgwuuwyeqngVZZYthY3sKTpgbDZd+08gJzcPgOz5OXTtPwHAQCpJkiSpSEiJuwDBY49B/R8H0WbR6wS33LLRPUX7DJ2yKoiulJObR5+hUwqyTEmSJElKGEdGYzZ/PvS8PZcx6VcT1tltk/YUnTk/Z7OOS5IkSVJh48hozHr3hmP/fJrtc6YS3H33Ju0pWisjfbOOS5IkSVJhYxiN0c8/w2P35dAr/XY44AA4+uhNul/nrPqkp6WucSw9LZXOWfULokxJkiRJSjin6cbo5pvhkryHqLJsJvR+eYNNi1a3skmR3XQlSZIkFVVBGIaxFtC0adNwzJgxsdYQh2++gUP3mc/MsjtTrtV+MHhw3CVJkiRJUsIFQfB1GIZN1z7uNN0YhCF07gy3pPel3D9/Qc+ecZckSZIkSUnlNN0YDB4ME4fPZkjp++Dkk6FJk7hLkiRJkqSkcmQ0yZYvh+uvh7sr30mpvKXQo0fcJUmSJElS0jkymmTPPANLJk7njNTHCc4/H3bZhQFjs21GJEmSJKlEMYwm0ZIlcOut8FS1bgSLU+HWWxkwNpuu/SeQk5sHQPb8HLr2nwBgIJUkSZJUbDlNN4nuuQeqzprAkX+8QHDFFZCZSZ+hU1YF0ZVycvPoM3RKTFVKkiRJUsFzZDRJZs+Gu++Gj7a7meDvSnDDDQDMnJ+zztuv77gkSZIkFQeOjCbJbbfB3v98RrNZ70QdjKpUAaBWRvo6b7++45IkSZJUHBhGk2DyZHjyiZD/bdsVtt0Wrrxy1XWds+qTnpa6xu3T01LpnFU/2WVKkiRJUtI4TTcJbrgBjik7lHozR8JDD0H58quuW9mkyG66kiRJkkqSIAzDWAto2rRpOGbMmFhrKEgffwyHtMpn1nb7sG3ZBfD991C6dNxlSZIkSVJSBEHwdRiGTdc+7shoAcrPh+uug0uqvs62s8bB888bRCVJkiSJBK0ZDYLgtiAIwrUusxNx7qLs1Vdh3JhcepW6GRo1glNPjbskSZIkSSoUEjkyOgVotdr3eeu5XYnwzz/QtSvcuv3TVPptGjzxDqSmbvyOkiRJklQCJDKMLg/DsMSPhq700EMw9+e/6Vz1djjgADj66LhLkiRJkqRCI5Fbu+wcBEF2EATTgyB4JQiCnRN47iLlzz/hzjvhwV0fouwfM6F3bwiCuMuSJEmSpEIjUWH0S+AcoA1wIVAT+CwIgqrrunEQBB2DIBgTBMGYefPmJaiEwmHA2GwaHfML4fy/OOHnO5nd/BA4+OC4y5IkSZKkQiUhYTQMw8FhGL4WhuG3YRgOA45ece6z13P7fmEYNg3DsGn16tUTUUKhMGBsNtc9PZVZX2RyU/Vbqbx0IZc26MCAsdlxlyZJkiRJhUoip+muEobhYmAisEtBnL+w6jN0CnM+2Ylt8+fQ6a8neafBwXxTdSf6DJ0Sd2mSJEmSVKgUSBgNgqAssBswqyDOX1j9/DMs/nYHbq/aldJ5y7j3oNMBmDk/J+bKJEmSJKlwSdQ+o32DIGgZBMFOQRDsB7wBlAeeTcT5i4rcr3dje37h3D9f4rU9jmBGlUwAamWkx1yZJEmSJBUuidraZXvgZaAaMA/4Atg/DMOfE3T+Qu+nn+CPsdvxSLVzCX4PefiAkwBIT0ulc1b9mKuTJEmSpMIlIWE0DMNTEnGeouyOOyAzZQ7n/fUKg/c+gpmVa5CZkU7nrPq0b5IZd3mSJEmSVKgkamS0RJs2DZ57Dt7fsy+lxuXS7uX/o129enGXJUmSJEmFVoE0MCppevSA7dJ+p9XkR+HUU8EgKkmSJEkbZBjdSj/8AC+8AE/veT8p/+TAjTfGXZIkSZIkFXqG0a3UvTtsW2Y+h056EI4/Hho2jLskSZIkSSr0DKNbYfJkePlleGafB0lZtBBuuinukiRJkiSpSDCMboXu3aF62UUcNvF+aNcO9tor7pIkSZIkqUgwjG6hiRPh1Vfh6f0eJeWvPx0VlSRJkqTNYBjdQrffDtXL/03Wd/fA4YfDfvvFXZIkSZIkFRmG0S3w7bfw+uvwvwOfJGXeXLjllrhLkiRJkqQixTC6BW6/HapVXEqbCXfDwQdDixZxlyRJkiRJRUqpuAsoasaNg/794Z2jniFlUDY890zcJUmSJElSkePI6Ga67TaoWimXthN6R+tEW7eOuyRJkiRJKnIMo5vh66/h7bfhqdYvkvrLDLj5ZgiCuMuSJEmSpCLHMLoZunWDqhl5HP1tz2hP0aOOirskSZIkSSqSXDO6iUaPhkGD4K2TXyf11anwxhuOikqSJEnSFnJkdBN16wbVquTT7ts7oGFDOO64uEuSJEmSpCLLkdFN8PnnMGQIvHnG26S+MBFefBFSzPGSJEmStKVMVJugWzeoVjXk2O/ugHr14KST4i5JkiRJkoo0w+hGjBoFH3wAjx07mNRx38CNN0IpB5QlSZIkaWsYRjeiWzeoUT2k/Xc9oHZtOOOMuEuSJEmSpCLPMLoBI0fC8OHw8AkfkTr6C+jSBdLS4i5LkiRJkoo8w+gG7LMP3HcfHDexB9SqBeeeG3dJkiRJklQsGEY3oHx5uKrpKFJHjoDOnaFs2bhLkiRJkqRiwTC6MXfeCdWrQ8eOcVciSZIkScWGYXRDvvoq2mD02muhXLm4q5EkSZKkYsMwuiE77QS33AKXXBJ3JZIkSZJUrLhh5oZUqwbdu8ddhSRJkiQVO46MSpIkSZKSzjAqSZIkSUo6w6gkSZIkKekMo5IkSZKkpDOMSpIkSZKSzjAqSZIkSUo6w6gkSZIkKekMo5IkSZKkpCsVdwGF2YCx2fQZOoWZ83OolZFO56z6tG+SGXdZkiRJklTkGUbXY8DYbLr2n0BObh4A2fNz6Np/AoCBVJIkSZK2ktN016PP0CmrguhKObl59Bk6JaaKJEmSJKn4MIyux8z5OZt1XJIkSZK06Qyj61ErI32zjkuSJEmSNp1hdD06Z9UnPS11jWPpaal0zqofU0WSJEmSVHzYwGg9VjYpspuuJEmSJCWeYXQD2jfJNHxKkiRJUgFwmq4kSZIkKekKJIwGQXBjEARhEAQPFcT5JUmSJElFW8LDaBAE+wMXAt8m+tySJEmSpOIhoWE0CILKwIvA+cBfiTy3JEmSJKn4SPTIaD/gjTAMhyf4vJIkSZKkYiRh3XSDILgQqAecuQm37Qh0BKhdu3aiSpAkSZIkFREJGRkNgqA+0BM4PQzDZRu7fRiG/cIwbBqGYdPq1asnogRJkiRJUhGSqJHRA4BqwHdBEKw8lgocHATBxUD5MAyXJuixJEmSJElFXKLC6ABgzFrHngamEo2YbnS0VJIkSZJUciQkjIZhOB+Yv/qxIAiWAH+GYfhdIh5DkiRJklR8JHyfUUmSJEmSNiZh3XTXFoZhq4I6tyRJkiSpaHNkVJIkSZKUdIZRSZIkSVLSGUYlSZIkSUkXhGEYbwFBMA/4OdYiNq4a8HvcRajE83mowsDnoQoLn4sqDHweqjAoCs/DOmEYVl/7YOxhtCgIgmBMGIZN465DJZvPQxUGPg9VWPhcVGHg81CFQVF+HjpNV5IkSZKUdIZRSZIkSVLSGUY3Tb+4C5DweajCweehCgufiyoMfB6qMCiyz0PXjEqSJEmSks6RUUmSJElS0hlGJUmSJElJZxjdgCAILg2CYHoQBP8EQfB1EAQt4q5JJUcQBLcFQRCudZkdd10q/oIgODgIgneCIMhe8bw7Z63rgxXPz5lBEOQEQTAiCILdYypXxdQmPA+fWcdr5BcxlatiKgiCrkEQfBUEwcIgCOYFQfBuEASN1rqNr4kqUJv4PCySr4mG0fUIguBk4AGgJ9AE+AwYHARB7VgLU0kzBdhutUvjeMtRCVEB+A64EshZx/XXA9cClwPNgLnAB0EQVExahSoJNvY8BBjGmq+RbZNTmkqQVsAjwIHAocByYFgQBFVWu42viSpordj48xCK4GuiDYzWIwiCL4FvwzC8cLVjU4E3wjDsGl9lKimCILgNOCEMw0Ybu61UUIIgWAxcFobhMyu+D4CZwENhGN654lg60Zuv68IwfDyuWlV8rf08XHHsGaBaGIZHx1WXSp4gCCoAC4D2YRi+62ui4rD283DFsWcogq+JjoyuQxAEpYF9gPfXuup9ok8kpGTZecUUtelBELwSBMHOcRekEm8noCarvT6GYZgDjMTXRyXfQUEQzA2C4IcgCJ4IgqBG3AWp2KtI9P75rxXf+5qoOKz9PFypyL0mGkbXrRqQCsxZ6/gcohccKRm+BM4B2gAXEj33PguCoGqcRanEW/ka6Ouj4jYEOAtoTTRFcl9geBAEZWKtSsXdA8A44PMV3/uaqDis/TyEIvqaWCruAgq5tecwB+s4JhWIMAwHr/79ikXoPwFnA/fGUpT0L18fFaswDF9Z7dsJQRB8DfwMHAX0j6cqFWdBENwLHAQcFIZh3lpX+5qopFjf87CoviY6MrpuvwN5/PcTrRr895MvKSnCMFwMTAR2ibsWlWgrOzr7+qhCJQzDmcBv+BqpAhAEwX3AqcChYRj+tNpVviYqaTbwPPyPovKaaBhdhzAMlwFfA4evddXhRF11paQLgqAssBswK+5aVKJNJ3rzter1ccVzswW+PipGQRBUAzLxNVIJFgTBA8BpRAHg+7Wu9jVRSbGR5+G6bl8kXhOdprt+9wLPB0EwGvgUuBioBTwWa1UqMYIg6Au8C/xC9AnrLUB54Nk461Lxt6JLX70V36YAtYMg2Av4MwzDX4IguB+4KQiC74EfgJuBxcBLMZSrYmpDz8MVl9uAN4neaO0I9CLqYPpWkktVMRYEwcPAmUB74K8gCFaOgC4Ow3BxGIahr4kqaBt7Hq54vbyNIvia6NYuGxAEwaVEe0dtR7TX2dVhGI6MtyqVFEEQvAIcTNRQax7wBXBLGIaTYi1MxV4QBK2Aj9Zx1bNhGJ6zYiuDbsBFwDZEzbY6hWH4XdKKVLG3oechcAkwgGgf8AyiN18fEb1G/pqUAlUiBEGwvjfKt4dheNuK2/iaqAK1sefhiu2EBlAEXxMNo5IkSZKkpHPNqCRJkiQp6QyjkiRJkqSkM4xKkiRJkpLOMCpJkiRJSjrDqCRJkiQp6QyjkiRJkqSkM4xKkiRJkpLOMCpJkiRJSjrDqCRJkiQp6f4fFnvqUtdS9O8AAAAASUVORK5CYII=\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           4.955748\n",
       "x1                  0.511913\n",
       "np.sin(x1)          0.526253\n",
       "I((x1 - 5) ** 2)   -0.021223\n",
       "dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we only have to pass the single variable and we get the transformed right-hand side variables automatically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    10.875812\n",
       "1    10.713220\n",
       "2    10.432218\n",
       "3    10.079836\n",
       "4     9.717983\n",
       "5     9.408289\n",
       "6     9.197015\n",
       "7     9.103727\n",
       "8     9.116508\n",
       "9     9.194871\n",
       "dtype: float64"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.predict(exog=dict(x1=x1n))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
