{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prediction (out of sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:57.728336Z",
     "iopub.status.busy": "2021-09-17T11:14:57.727898Z",
     "iopub.status.idle": "2021-09-17T11:14:58.223589Z",
     "shell.execute_reply": "2021-09-17T11:14:58.223973Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:58.227332Z",
     "iopub.status.busy": "2021-09-17T11:14:58.226425Z",
     "iopub.status.idle": "2021-09-17T11:14:59.027188Z",
     "shell.execute_reply": "2021-09-17T11:14:59.026600Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "plt.rc(\"figure\", figsize=(16, 8))\n",
    "plt.rc(\"font\", size=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Artificial data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:59.033454Z",
     "iopub.status.busy": "2021-09-17T11:14:59.033002Z",
     "iopub.status.idle": "2021-09-17T11:14:59.036254Z",
     "shell.execute_reply": "2021-09-17T11:14:59.035450Z"
    }
   },
   "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, 0.5, 0.5, -0.02]\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + sig * np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:59.040344Z",
     "iopub.status.busy": "2021-09-17T11:14:59.039907Z",
     "iopub.status.idle": "2021-09-17T11:14:59.058162Z",
     "shell.execute_reply": "2021-09-17T11:14:59.057772Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.987\n",
      "Model:                            OLS   Adj. R-squared:                  0.987\n",
      "Method:                 Least Squares   F-statistic:                     1206.\n",
      "Date:                Fri, 17 Sep 2021   Prob (F-statistic):           1.03e-43\n",
      "Time:                        11:14:59   Log-Likelihood:                 7.6619\n",
      "No. Observations:                  50   AIC:                            -7.324\n",
      "Df Residuals:                      46   BIC:                            0.3242\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          4.9421      0.074     66.995      0.000       4.794       5.091\n",
      "x1             0.5012      0.011     44.054      0.000       0.478       0.524\n",
      "x2             0.4840      0.045     10.823      0.000       0.394       0.574\n",
      "x3            -0.0202      0.001    -20.231      0.000      -0.022      -0.018\n",
      "==============================================================================\n",
      "Omnibus:                        0.486   Durbin-Watson:                   2.338\n",
      "Prob(Omnibus):                  0.784   Jarque-Bera (JB):                0.395\n",
      "Skew:                           0.208   Prob(JB):                        0.821\n",
      "Kurtosis:                       2.874   Cond. No.                         221.\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "olsmod = sm.OLS(y, X)\n",
    "olsres = olsmod.fit()\n",
    "print(olsres.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In-sample prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:59.062190Z",
     "iopub.status.busy": "2021-09-17T11:14:59.061726Z",
     "iopub.status.idle": "2021-09-17T11:14:59.067465Z",
     "shell.execute_reply": "2021-09-17T11:14:59.066858Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.43689064  4.91270848  5.35022686  5.72306605  6.01436661  6.21955939\n",
      "  6.34711618  6.41715787  6.4581486   6.5022192   6.5798882   6.7150477\n",
      "  6.92103797  7.1984557   7.53505593  7.9077639   8.28646614  8.63895805\n",
      "  8.93623501  9.15725774  9.29240905  9.34507386  9.33108309  9.27611268\n",
      "  9.21146513  9.16892654  9.17554429  9.2491832   9.39559042  9.60745115\n",
      "  9.86559048 10.14212377 10.40503796 10.62345132 10.77268715 10.83832814\n",
      " 10.81858506 10.72459075 10.57857103 10.41019301 10.25169096 10.1325702\n",
      " 10.07475956 10.08900847 10.17312089 10.31231608 10.4816571  10.65014779\n",
      " 10.78582543 10.86101207]\n"
     ]
    }
   ],
   "source": [
    "ypred = olsres.predict(X)\n",
    "print(ypred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a new sample of explanatory variables Xnew, predict and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:59.071631Z",
     "iopub.status.busy": "2021-09-17T11:14:59.071193Z",
     "iopub.status.idle": "2021-09-17T11:14:59.074899Z",
     "shell.execute_reply": "2021-09-17T11:14:59.074277Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10.84402874 10.69880682 10.4443285  10.12385184  9.79431964  9.5124179\n",
      "  9.32069706  9.23715423  9.25082681  9.32447655]\n"
     ]
    }
   ],
   "source": [
    "x1n = np.linspace(20.5, 25, 10)\n",
    "Xnew = np.column_stack((x1n, np.sin(x1n), (x1n - 5) ** 2))\n",
    "Xnew = sm.add_constant(Xnew)\n",
    "ynewpred = olsres.predict(Xnew)  # predict out of sample\n",
    "print(ynewpred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:59.079253Z",
     "iopub.status.busy": "2021-09-17T11:14:59.078813Z",
     "iopub.status.idle": "2021-09-17T11:14:59.510597Z",
     "shell.execute_reply": "2021-09-17T11:14:59.511274Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f05cc40e820>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6MAAAHWCAYAAACCFl9HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAB3YklEQVR4nO3deZzN1R/H8dcxBmMd+zKI0ES2yaCiyNJUKkvhVylKikiUKZJEZRtFC4rKUpJky1LWlKg0QpYoa4x9GetgzHx/fxxjnbHeud97Z97Px+M+mO+99zuf0bc7933POZ9jHMdBRERERERExJsyuF2AiIiIiIiIpD8KoyIiIiIiIuJ1CqMiIiIiIiLidQqjIiIiIiIi4nUKoyIiIiIiIuJ1CqMiIiIiIiLidRndLiBfvnxOiRIl3C5DREREREREUsHSpUv3Oo6T/8LjrofREiVKEB0d7XYZIiIiIiIikgqMMVuSO65puiIiIiIiIuJ1CqMiIiIiIiLidQqjIiIiIiIi4nUKoyIiIiIiIuJ1CqMiIiIiIiLida53072cQ4cOsXv3buLj490uRdKQwMBAChQoQM6cOd0uRUREREQkXfLpMHro0CF27dpFSEgIQUFBGGPcLknSAMdxiIuLIyYmBkCBVERERETEBT49TXf37t2EhISQNWtWBVHxGGMMWbNmJSQkhN27d7tdjoiIiIhIuuTTYTQ+Pp6goCC3y5A0KigoSNO/RURERERc4tNhFNCIqKQaXVsiIiIiIu7x+TAqIiIiIiIiaY/CqIiIiIiIiHidwmgqaNWqFcYYjDFnthC5++67GTJkyFWtUVywYAHGGPbu3ZuK1YqIiIiIiHifwmgqqVevHjt27GDz5s3Mnj2bBx98kJ49e3LnnXdy9OhRt8sTERERERFxVboIo1OWxVCj33xKdp1BjX7zmbIsJtW/Z+bMmSlUqBAhISFUrlyZl156iQULFvDnn38yYMAAAL788kuqVq1Kjhw5KFCgAE2bNj2z9+XmzZu5++67AcifPz/GGFq1agXADz/8wJ133knu3LnJkycPERER/P3336n+M4mIiIiIiHjKFYVRY8xdxpjvjDExxhjHGNPqgvubGGNmGWP2nL6/dirUek2mLIuh26SVxMTG4QAxsXF0m7TSK4H0QuXLl+fee+9l4sSJAJw8eZJevXqxYsUKpk+fzt69e3n00UcBKFas2JnHrV69mh07dvD+++8DcPToUTp16sSSJUtYsGABuXLl4sEHH+TkyZNe/5lEREQ8wY0PjkVExF0Zr/Bx2YFVwJjTtwtlAxYDX6Zwv2uiZq0jLj7hvGNx8QlEzVpHo7AQr9dTrlw55s6dC8DTTz995viNN97IsGHDKFu2LNu2baNo0aLkyZMHgAIFCpAvX74zj3344YfPO+fIkSPJmTMnS5YsoWbNml74KURERDwn6YPjpN/XSR8cA678rhYREe+4opFRx3FmOo7zmuM43wKJydz/heM4vYDvPV3g9doeG3dVx1Ob4zhn9rf8888/adiwITfccAM5cuQgPDwcgP/++++S59iwYQOPPfYYpUqVImfOnBQsWJDExMTLPk9ERMQXXeqDYxERSbuudGTUbxUJDiImmeBZJDjIhWpgzZo13HjjjRw9epSIiAjq1avHF198QYECBdi7dy933nnnZafbPvDAAxQtWpRPPvmEkJAQMmbMSLly5TRNV0RE/JKvfXAs6cuUZTFEzVrH9tg4igQHERkRqhF5ES9xpYGRMeZZY0y0MSZ6z549qfq9IiNCCQoMOO9YUGAAkRGhqfp9k7Nq1Sp++OEHHnnkEdauXcvevXvp06cPd911FzfffDO7d+8+7/GZMmUCICHh7KfF+/btY+3atbz22mvUq1ePsmXLcvjwYU6dOuXVn0VERMRTUvqA2K0PjiX98KXeIiLpkSth1HGc4Y7jhDuOE54/f/5U/V6NwkLo26QCIcFBGCAkOIi+TSqk+ideJ06cYOfOnWzfvp0VK1bw3nvvUbt2bapUqUKXLl0oXrw4mTNn5qOPPmLjxo3MmDGDHj16nHeOG264AWMMM2bMYM+ePRw5coTcuXOTL18+RowYwfr16/npp59o27YtGTOm+UFuERFJo3zpg2NJX653injiRYvXRORqpIutXRqFhbCoax029WvAoq51vDL1Yu7cuRQuXJjixYtTt25dvvvuO958801+/vlnsmXLRv78+Rk9ejRTpkyhXLly9OrVi/fee++8c4SEhNCrVy+6d+9OwYIF6dChAxkyZGD8+PH89ddflC9fnvbt2/PWW2+ROXPmVP+ZREREUoNbHxyLXMsU8f374eOP4Y474J13UqsykfTBOI5zdU8w5gjQwXGcUcnclw/YA9ztOM6CKzlfeHi4Ex0dnex9f//9N2XLlr2q+kSuhq4xERGR9KtGv/nJ9hYJCQ5iUdc6Z74+eRJ++AHGjIFp0+zXt9wCL78MTz3lzYpF/JMxZqnjOOEXHr+iuZ3GmOxA6dNfZgCKG2MqA/sdx/nPGJMHKA4En35MaWNMLLDTcZyd11m7iIiI+JEjR2DtWsiXDwoVgixZ3K5IJHmREaHnbSsEZ6eIOw5ER9sAOm4c7NsHBQrA88/Dk09C5cpweoMEEblGV7rQMBz48Zyve52+jQZaAQ8BI8+5f8Q5j3vzuioUERERn7dtmx0x+u47mD/fjhwlyZULChe2wfTC27nH8+WDDOliAZH4iqSp4Od2021VsRx//1CIbo/ZD1UyZ4ZGjWwArV8fAgPdrVkkLbmiMHp6ym2Kn/2cnrI7yiMViYiIiM9zHFi2zIbP776zfwcoXRo6dIAaNeDgQdi5E3bssH/u3AlLl9o/Dx+++JyZMsHDD0PnzlC1qnd/Hkm/GoWFULd0CBMnwhdfwHOv2ev7zjvtNNxHHoHgYLerFEmb1IJVRERErsjx4/DjjzZ8TpsGMTF2muIdd0D//vDQQxAaemVTF48ePRtQk25r1tgwMG6cDbOdO9sRqYCAy55O5JokJMDQofDaa3Z6ealS8Oab0KIF3HjjZZ586JD9nyI42H6SIiJXTWFUREREUhQbC1Om2AA6e7YNkdmyQUSEDZ/33w/Xsktbtmz2jX+pUucf79sXPv8cPvjAjkiVKAEvvACtW9vpvpI+TFkWc97U2ciIUI93V/77b3td/for3Hsv9OgBt99+mQ9TDh60/0OMHw9z5kDSPu9Zs9pQerW3YsW0qFrStavuputp6qYrbtI1JiKSsu+/h6eftqOWISE2fD74INx9dwrvn+Pj4ZdfYPp0WLEC8uSBggVTvgUFpfi9ExJsAB40CBYuhOzZbS0dO14cYCVtmbIsJtmmQp7a7ic+3o7kv/UW5MgBgwfD449fIoQeOWIvxvHjbUvdkyfhhhugWTMoXtx+YnO5W0JC8ufOmxdefx3atbOLU0XSqOvqpisiIiLpx9Gj0KWL3UuxfHn49ls7FTfZN+t799rUOn26faN+6JCdslipEmzdCrt22dGk5GTPfn44LVLEDofWqkVAgKFxY2jc2K4zHTzYTqf88ENo2NBO4d2XPYaBs1N39Ey8L2rWuvOCKEBcfAJRs9Zd93/f6Gg7GvrXX9C8uR2BL1AgmQceOwYzZ9oAOn26nY5bpIhtpfu//0G1alfeStdx7P9UScH0wAH75/79MHasvZgHD4bevW0q1rx0SUc0Mirpmq4xEZHz/fYbPPEEbNhgm7e89dYFo6COA6tX2zfo06fbOY6JibYdboMG8MADUK+eDZpJjh+H3bttME36M7nbf//ZUagKFWwXpMcft/N5T9u+HYYMsSF5/37IXOgg2atsIlvZ7ZgAx6OjZ+Kekl1nkNy7UwNs6tcgxeddamrvsWN2Lei779pLdehQ+6HGeU6csB+ojB9vR0KPHrVJtWlTm1xr1Eidds9z50LXrvZTlwoV7Fz1++/XvjGSpqQ0MqowKumarjERESs+3g7M9Oljl7GNHg21ap2+8/hxWLDgbADdssUer1LFhs8HHoBbb73+N+pxcbZ70YcfwvLldk1d69Z2NOqcbjLHjkGFFmvZ+nMI8ftyEJj/EPkarCBTwUOEBAexqGud66tDXFWj33xiYuMuOn6p/7aXmtobfDCENm1g/Xpo0wYGDLigO+6WLfbinzjRjuLnyWPbOjdvbv8nyOiFiYSJiTBhAnTvbj8JuusuO5f4tttS/3uLeEFKYVS7eYmIiKRzf/9t3/O+/bbdS/Gvv04H0cOH7Qhlvnxw330wciRUrgwjRthWutHRdrgpPNwzI0ZBQXZh6J9/2oWi99xjpy+WLm0XrM6ZA45D1qyQeNMGCrf+mfyNo0k8lokdY2oQu7g0MfuPX38d4qrIiFCCAs+fqhoUGEBkRGiKz0luau/RI4a2be0a58REmDcPhg8/J4gmLR4tV85+CNK4sZ1yvnOnfWDdut4JomD//2ne3P7POGQIrFtnuyk1bmyPiaRRCqMeZoy55K1Vq1ZulygiIgLYN+gffGAHNf/7DyZNsnkzZ07sSGjFinY+Y/Pmdv3c3r22k+gzz9j1c6nFGKhZ006X3LLFjhb9/rsNp+XKwZAhlA5yMAay3rSLwq1/JmvoDg4uDGXfuBqsXZt6pUnqaxQWQt8mFQgJDsJgR0QvN/16+wUjqcfWF2D7p3ex648ivPwyrFwJdc4dVF24EMLC7PTY+vVt4Bs50rbVDQxMnR/sSgQG2pkA69fbOfLz5tmF2888A9u2uVeXSCrRNF0P27lz55m/T58+nTZt2rBjx44zx4KCgsh1Tm/6+Ph4At180Uvn/PEaExHxhG3boFUr+163QQP49FO7lo5jx6BbN5tSS5eGUaPsWjm3nTgB33xjp/D+8Qfx2bIzrlxdRla+n015bEiJ/zeEQ/MqcupEBvr1s1vCpMYSP/E9SVN7E+IC2T/nFo79HUJg/kPc3HQdfw2pevaBe/bAK6/Y6/qGG+z19OCDrtV9WXv22LnzQ4fai7ljRxugc+d2uzKRq6Jpul5SqFChM7fg0/NAkr4+fvw4wcHBjBs3jjp16hAUFMQnn3zCqFGjyH5uowdgwYIFGGPYu3fvmWOLFy+mVq1aZM2alZCQENq1a8ehQ4e8+eOJiEgaMG6c7ZPy22/wyScwbdrpILp4sZ2G+8EHNsktX+4bQRTsthdPPAFLlsBvvxHYqCEtls3kxxHPMWj6u5QNPMmQ1/Pzz98ZqFsXOnWysyw3b3a7cPGGyIhQzP5gdo6uybF1hclVcx03PvMrvZ85PYKfmGg/cbn5ZvjyS3j1VduIy5eDKNhNfAcNstN2mzaFqCi7fnrAALvGWsTPKYy6oFu3bjz//POsWbOGRo0aXdFzVq5cyT333MNDDz3EihUrmDRpEsuXL+fpp59O3WJFRCTN2L/f7krx2GNQtqzdCvTZZ8GcOG5Hi+680+6hOH++DaTndLL1KdWrw5dfkmHbVujWjcbrFvL9iHY02vQ7hQvbcP3ZZ7Y5acWK9u8uTwSTVHZ8XQjbv7idDE4GCj++mFsaxNCvaXk7tfevv+y07zZt4JZb7Ics/fr57vWdnBIlYMwYW/sdd9gwXaaMHeHVxS1+zO/2Ge3Uyf5/6E2VK9v+CZ7ywgsv8Mgjj1zVc6KiomjevDkvv/zymWPDhg0jLCyM3bt3UyDZTbJERESszZuhdm3bd+idd2z2zJgR24SoZUtYs8Ym04EDIUcOl6u9QgUL2imMzZvDU0/ZDqjNmmE++oinn85PnTr28DPPwOTJtu/S7ztT3v5D/M+pU/Daa3bAsGbNDEyYkIVChWraO48csRvmDh5sp7WOHGmvdX/eMqViRZgxA37+2QbSp56yC2IHDvTvn0vSLb8Lo2lBePhF06Uva+nSpaxfv57x48efOZa03nfDhg0KoyIikqItW2wQPXQIFi2CatWwI6A93rJ7GhYqZPdXjIg473mX2rcxJdfynOtWqZJtcDRgAPTqZUd2hwyhRNOmzJtn+Ogj+779ppsTyVl3D4E32emNMbFxdJu0EkCB1A/t22dH+ufOtT1/Bg2CTJmwI4WTJ8OLL9rF0W3a2Os8b163S/acu+6y/zO/+CK89579n/vjjyEg4PLPFfEhfhdGPTlC6ZZsF0wLyZAhAxc2koqPjz/v68TERJ555hk6d+580flCQvQLVEREkrd1q93a4uBB26zo1lux83Nbtjz75+DBF2y8ePG+jVcS3K7lOR4TGGi77jZsaEeLmjeH8ePJMHQoHTsWJCICwusfZvvkymS9uQB57llFQFA8cfEJRM1apzDqZ1assLuexMTYadhnVi1t2wZt29rRw4oVbdOr22/3+Pd35UOXC2XIYKfT58plpzscPmyn8mbK5N06RK6D1oz6gPz583Ps2LHzmhEtv2Au8q233srq1aspXbr0RbegoCAvVywiIv4gJsYG0X37YPZsuLXiKbuZaNWqdi/FqVPtmrMLgigkv29jUnBLybU8x+PKl4dff7VrAmfMsFvBjB1L6E0Oef+3iOC71nLsn0Ls+PxOTsQEAxdvCyK+bfx4my9PnLCzVc8E0Vmz7NqqBQvg3XftouFUCqLdJq0kJjYOh7MfukxZFuPx73VZxtj/pwcMsP8wjRursZH4FYVRH1C9enWyZctGt27dWL9+PRMnTmTo0KHnPebVV19lyZIltG3blmXLlrF+/XqmT5/Oc88951LVIiLiy3bssEF09277Hn330VWsLFUJevRgTtmazBg3Fx56KMXnpxTQLhXcruU5qSJjRjsvd/lyuOkmaNECGjWiYsaj5Lp9A4WfXIQJcNj51e0cXlacwrn0oa4/SEiwa53/9z87wr90qe1lRUIC9OwJ990HhQvbO1566fSiaM/ziQ9dLhQZaVtjf/+9/XfQbgviJxRGfUCePHkYO3Ysc+bMoUKFCgwfPpy33nrrvMdUrFiRn3/+mc2bN1OrVi0qVapEt27dKFiwoEtVi4iIr9q50wbRHTvsUtB9caso+Wgjbor5hxcejKTNfS/TZf62S47kFAlOPqCldPxan5Oqbr4ZfvnFjpLNns23H7Xhf2vmk6nAQQq1/IUsJfayf3YFsvxaXYNJPm7/fpuxoqKgXTu7LLhQIeynLffeC7172ynnv/8OoaGpWovPfOhyoWefhbFj7VrSunXhnO0BRXyVuXCtoreFh4c70dHRyd73999/U7ZsWS9XJOmJrjERSWt277bNiv77zwbRmjfv5d+Kt1F8938817g7C0qdbaIXEhzEoq51kj3Phes/AYICA+jbpMIVrxm9kud4zb//2vmcv/zC4tDqvFynHaZoMYptDuebETkJC4OJE6FkSXfLlIv99Rc0amSnnQ8dCq1bn77jl1/s2uD9+2HIkHPm66auGv3mE5NM8LzU/09eNX06PPIIlCoFc+ZAkSJuVySCMWap4zgXdXHVyKiIiEgasWePHRDZvNkul6wZugfq1KHYnq20afL6eUEULj2S0ygshL5NKhASHITBvtG+XKi8lud4TZky8NNP8P773LF1Jb9+3ZnFleIYPzwn06bBxo0QHm6nNIvvuHB9aOvW2G65775rP3XJmhV++81rQRQgMiKUoMDzu9YGBQYQGZG6I7JX7IEH7HTd//6z+6tu3Oh2RSIp0siopGu6xkQkrdi3D+rUgX/+sUG0TvndNpmuX8+Lj/dmar5yFz3HZ0ZyvG39emjWDJYtg27doHdv1m/OSJMmsGoVvPWWPZxBH9m7JiHB7h86YADUqAHffnt6Wm5srO2WPGWK3Vf2s89sN1kv84luupezZImdwhwUZEdIy138GiDiLSmNjPrd1i4iIiJyvv37oV49WLcOpk2DOrfsgrvrwKZNMGMGd+cOZXYy02d9ZiTH20qXhsWL7R6NffvC4sWUHjeOX38tTJs28Prr9n38mDE25/hF8EhD9u+HRx+1HaDbtbM7D2XKBPz5JzRtakf8Bg2y//2McaXGRmEhvn8NVKtmh5Pr17f7ks6aBVWquF2VyHn0mZ+IiIgfO3AA7rkH1qyxg0X1K5zuXrR5M8ycCXXq+Pb0WbdkyWK7j44ZA3/8AWFhZPtjAWPH2vAzc6bdAef9Cbt8ZxuPdGDlSvvvvmABjBhh14hmCnRg+HC44w44edIGrE6dXAuifqV8eVi4ELJnt68LP//sdkUi51EYFRER8VMHD0JEhG3wMnky3Ftx+9nuRd9/b/9+WqOwEBZ1rcOmfg1Y1LVO+g6i53riCTsMmjs31K2L6deXF19IZP58uzvGS4/nY+9fBc57iuvbeKRR335r14fGxdkw+swzwNGj8OST8Nxz9npetixV9g5N00qXts2eQkLsC8b337tdkcgZCqMiIiJ+6NAhuxxs+XLbBfb+SjH2zXpMjG2je9ddbpfoP265xY6ONm9uFyo++CB3ltvHn39CYIGD7P3uVvbPK4uTcHYkzvVtPPzIlGUx1Og3n5JdZ1Cj3/yLRpUTEqB7dzsDt2JFiI4+nTfXrrUbiY4da7dumTkT8uVz54fwd0WL2lHRsmXt/sLffON2RSKAwqiIiIjfOXYM7r/fvmkfPx4eDNtmg+jOnXZdWM2abpfof7Jnt6Fn6FCYOxduvZUi25YQ1nYFOW7dzOHoG9k1vjqnDmcGXNw71c8kbfeT0jTn2Fibjfr0gTZt4McfT+9E8vXXtr3x7t128WiPHuoodb3y57f/wLfdZhflfvaZ2xWJKIyKiIj4E8eBVq1s/51x46Bxlf+gVq2zb9rvuMPtEv2XMbZjzqJFNvjUrMnHhxcRct/f5H1gOSd35mLHyDtJ+K9g+m3+dJWiZq07r3EWnJ3mvGaN7bEzezZ8/LFdFpqZE9Chgw1LlSvbabn16rlTfFqUK5f9wKp+fTsPetAgtyuSdE5hVERExI/06QMTJkC/fvBI1S12RHTfPrt1w223uV1e2hAebju33nsvlQb0YM6SYYSFb6Pwk4vIkiuebePCWTwuhPj4S5/mctNT04OUpjP/uyQn1avb6eY//miXhLJli51ePmQIvPyyvSNEa5s9LmtWmDrVbo3z0kt2NoCISxRG5bp16NCB2uc0yWjVqhUPPPDAdZ3zzTffpHz58tdZmYhI2jJ1qt125PHHIbLpZhtEDxyw00qrVXO7vLQld27bnrh/f4rOn8mcr14m5sU87N2YnbZtISoK7rzT7p6TnMtNT00vLpzO7DgQu7AMeyaHU66cnWpesya2qc6tt9p1ohMnwsCBEBjoTtHpQebMdir0Aw9Ax452eFrEBQqjqSQmJoZnn32WokWLkilTJkJCQmjTpg3btm0773GXC24rVqygYcOGFCpUiCxZslC8eHEefvhhtmzZkto/wjV7//33+fLLL6/osZs3b8YYQ3R09HnHu3Tpwk8//ZQa5YmI+KVVq6BFCztoN6LnNkztWrad7rx59qB4XoYM8MorMH8+HD4M1aoRNPpjhg11+OYb+PtvCAuzXWAvdKnpqelJZEQoQYEBACSeyMieSVU4uPgm6jx4lJ9+gqKFE+x60Pvvh2LFYOlSaNLE5arTiYwZ4auvbAOvpk3t/lAiXqYwmgo2bdpEeHg4q1atYvTo0axfv54vv/yS1atXU7VqVTZv3nxF59mzZw9169Yle/bszJgxg7Vr1/LFF19QqlQpDh065NGaT5486bFz5cqVi+Dg4Os6R/bs2cmbN69nChIR8XP79kHDhrbHztSvjhLU/CE7Ijpvnh1NktR111127WKtWnZNaePGNL17L8uXQ2iofR///PNw/PjZp6Q0PTW9deFN2uM2T3xedo6pQdzGArR5JZa5U7OR5dBuu9XI22/D00/Dr7/abUjEe3LkgGnTICjIjpLu2eN2RZLOKIymgvbt25MhQwbmzp1L3bp1KV68OHfffTdz584lQ4YMtG/f/orOs2jRIg4cOMDIkSOpUqUKJUqUoFatWgwYMIAKFSqk+Lyk0da3336bggULkj17dp566ini4s7+Aqxduzbt2rWjS5cu5M+fnxo1agCwZs0aGjRoQI4cOShQoACPPvooO3fuPPO8hIQEunTpQu7cucmdOzedOnUiISEh2e+fxHEc3n33XcqUKUPmzJkpWrQo3bp1A6BkyZIAVK1aFWPMmem+F07TTUxM5K233qJYsWJkzpyZChUqMHXq1DP3J42wTpw4kfr165M1a1bKlSvHnDlzrujfWkTEV8XHQ7NmsG0bTJ6YSJGuT8KKFbaNbliY2+WlHwUL2q1FBg2yU0orVqTkxnksXAhdusCwYXYXkrVr7cNT6rab3rrwJiZCzOIQ/vnkNoIDsrNgfgaG9w/GLF5kr99Fi2xX188+s4FIvK94cbsGYMcOaNwYTpxwuyJJRxRGPWz//v388MMPtG/fnqxZs553X9asWXn++ef5/vvvOXDgwGXPVahQIRITE/n2229xHOeq6vjpp59YsWIF8+bNY+LEicyePZtXX331vMd8+eWXOI7DwoULGTNmDDt27OCuu+6ifPnyLFmyhLlz53LkyBEaNmxIYmIiAO+++y4jRozgk08+4ddffyUhIYGxY8despbXXnuNt956i27durF69WomTJhAsWLFAFiyZAkAP/zwAzt27GDSpEnJnuP9998nKiqK/v37s3LlSho3bkyTJk1Yvnz5eY/r3r07HTt2ZMWKFVStWpX//e9/HDly5Kr+7UREfMnLL9tZosOHw23TX4dJk+C99+C++9wuLf3JkAE6dYLff4fgYKhfn0w9XiXqnZPMmAHbt0OVKjB69PnTU5MEBQakqy68GzZAnTq2OW6NGnYGbq27HHv91qplG+n89psdFRV3Va9uL9xFi+weO1f5vlPkWmV0u4Cr1qmT3eHbmypXhsGDr+ih//77L47jULZs2WTvL1euHI7j8O+//1LtMs0mbrvtNl577TVatmxJ+/btqVq1KrVr1+bxxx/nhhtuuORzAwICGDlyJNmzZ6d8+fL079+f1q1b07dvX7JlywbYUcl33333zHPeeOMNKlWqRP/+/c8cGzNmDHny5CE6Oppq1aoxePBgXnnlFZo1awbYkDhr1qwU6zhy5AiDBg1i8ODBPH36l03p0qW5/fbbAcifPz8AefPmpVChQimeZ+DAgXTp0oXHHnsMgN69e/Pzzz8zcODA89andu7cmQcffBCAPn36MGbMGJYvX05N7bknIn7o00/hww9tw8uWjIa+fW3b0Y4d3S4tfatc2XbeefllGDAA5s3j/q++Yvnym2jRwm6988QTIbzxnOGjhWvZHhtHkeAgIiNCaRSW9rvDJiba67ZbN8iUyQ56PvUUmEMH4eGnYPJkuy7088/tViPiG5o1g3Xr4I034Oab4bXX3K5I0gGNjPq4d955h507dzJ8+HAqVKjAZ599Rrly5Zg3b94ln1exYkWyZ89+5uvbb7+dkydPsmHDhjPHqlSpct5zli5dys8//0z27NnP3JJGMDds2MDBgwfZsWPHmSAJkCFDBqpXr55iHWvWrOHEiRPUrVv3qn7ucx06dIjt27efmUqcpGbNmqy5YLF9xYoVz/y9SJEiAOzevfuav7eIiFt++cWuQ7znHuj/4C92tKJuXfsu3xi3y5OsWe3c3MmTbUvdsDBCZn3O3DkOb74JY8fCm08XYci9ddjUrwGLutZJF0H0n3/sEttOneyo6OrVduDT/LXCNtr67jt4913b9UlB1Pcktevu3j35zlwiHuZ/I6NXOELpltKlS2OMYc2aNTRu3Pii+9esWYMxhtJXsUA/b968NG3alKZNm9K3b1/CwsJ46623rivgAWdGSJMkJibSoEEDBg4ceNFjCxYseGaqri8xF7whCzynDXzSfb5Yt4jIpfz3n90CsEQJ+KbfRjLe0xhKlrQbjGq7C9/SqBFUrQpPPgmtWxPw/ff0HD6c2rVz89hjduvXLl3g1Vdtr5i0KiHBLqft0cMu/RwzxnZ/NieOQ8++dlQ/f35YsOD0Xi7ik4yxUzI2bbLX9A032OtbJJVoZNTD8ubNS0REBEOHDuXYsWPn3Xfs2DGGDBnCfffdR548ea7p/JkyZaJUqVKXXQe5cuVKjh49eubr33777cxzU3LrrbeyevVqbrjhBkqXLn3eLUeOHOTKlYvChQvz22+/nXmO4zhn1n0mp2zZsmTOnDnFkdxMmTIBXNQE6Vw5c+akSJEiLFq06Lzjv/zyC+XKlUvxeSIi/ujYMZtvjh+HaV8eJFeLB+07/enT7d6X4ntCQmDOHOjf3+5NWqkStczPrFhhP1R45x3bJHbYMDh1yu1iPW/NGrsmNDLSNsddvRqeeALMzz/ZKc29e9uWw8uWKYj6gyxZ7Ih/wYLw0EOwdavbFUkapjCaCj766CNOnTpFvXr1mD9/Plu3bmXBggXUr18fx3H46KOPznv8oUOHWL58+Xm3zZs3M336dFq0aMH06dP5559/WLduHQMHDmTmzJnJjrqe69SpUzz99NOsXr2aOXPm0LVrV9q0aXPRaOi52rdvz8GDB2nevDm///47GzduZO7cuTz77LMcPnwYgBdffJEBAwbw7bffsm7dOjp16sSOHTtSPGeOHDl48cUX6datGyNHjmTDhg0sWbKEYcOGAVCgQAGCgoKYNWsWu3bt4uDBg8meJzIykoEDBzJu3Dj++ecf3njjDRYuXEiXLl0u+e8gIuJPHMeurVu+HL7+8hShPf9n5z1OnAhlyrhdnlxK0p6kv/5q38zXrk2+93swdlQ8S5bYJXjPPw8VKtiZqmmhP8ypU3bAMywM1q+3W1ZOngyFM++H1q2hdm04eRJ++MHOWy5QwO2S5UoVKGA/ADt61AZSNYOU1OI4jqu3KlWqOClZs2ZNivf5uv/++8955plnnCJFijgZM2Z0Chcu7LRu3drZunXreY9r2bKlA1x0e/jhh50NGzY4zz33nBMaGupkzZrVyZkzp1OpUiVn0KBBTmJiYorfu2XLlk6DBg2cXr16Ofnz53eyZcvmPPnkk87Ro0fPPKZWrVpO+/btL3ruP//84zz88MNOcHCwkyVLFuemm25yOnTo4Jw4ccJxHMeJj493OnXq5OTKlcvJlSuX06FDB6dt27ZOrVq1Lvr+SRISEpy+ffs6JUuWdAIDA52iRYs6r7322pn7R4wY4RQrVszJkCHDmfP07NnTueWWW847R+/evZ2iRYs6gYGBTvny5Z3JkyefuX/Tpk0O4Pzxxx/n/TyAM2HChBT/rfz5GhORtOfttx0HHKd/f8dxOna0Xwwf7nZZcrUOH3acp5+2//2qV3ec5cudxETHmTrVcUJD7eFatRzngl9ZfuWvvxynShX7szzyiOPs3Ok4TmKi44wd6zgFCjhOQIDjvPKK45zz3kP80PffO06GDI7z0EOOc+qU29WIHwOinWSyoHFc/mguPDzciY6OTva+v//+O8WutJKyVq1asXfvXqZPn+52KT5P15iI+IqpU+303Mcfhy9qfIx5vh107my3wRD/NGECtG0LBw7YBZS9exMfUoJPP4WePWHPHnjsMTuNt0QJt4u9MgcPwvvvw9tv291thgyxM3DZuBHatYPZs6FaNbsXUaVKbpcrnvDRR/DCC3bxc1SU29WInzLGLHUcJ/zC45qmKyIi4rJVq2xWCQ+Hzx6di3mhAzRooDd+/q5pUzt/9ZVXbDANDSXw1Zdo13Qv69fbhqWTJ0NoqH1IbKzbBSfPcezWqk8/DYUL2yD98MN2bWjTRvF2e5vy5WHxYvjgA/ungmja0aEDtG8PAwfa5kYiHqQwKiIi4qJ9+6BhQ8ieHaYPXEvmxx+BsmXtAryAALfLk+uVOzf06wf//mu7+rz/PpQqRc6P+vB2t6P8848dHR04EEqVspsGnDzpdtHWoUO26VJYmO0K/M039kOT6GgYNw7yb1piP0F59VXbuejvv+0Imq7btGfwYPvfuF07+PFHt6uRNERhNA0aNWqUpuiKiPiBhAQ7LXfbNpg2ah8FWz8AmTLBtGmQM6fb5YknFS1qR5VWroS777bDomXKUHTmcEaOOMWyZVClip2ZXbasnZ3911/gxu5k0dF2W9siRWzTJWNsKN2+3c6+rXLTYejY0SbUvXth0iQ7xFu0qPeLFe/ImBHGj4ebbrLD4v/843ZFkkYojIqIiLikVy+YNQuGDDpJeN+H7RYKU6b4zwJCuXrlytn/xgsX2r1jn3sOypen0sbJzJ7l8MMPdi3myy/bma6FC9uR088+gy1bUq+sw4dPB80qdlvJr76C//0PliyBP/+0S19zHthiL9qbb7brCNu3t6Ohl+nwL2lErly2w25AADzwAOzf73ZFkgb4fBh1u8GSpF26tkTETTNmwFtvQf5bt5P94wfhp5+IfmMg3HHHZZ87ZVkMNfrNp2TXGdToN58py2K8ULF4VM2a8MsvNphmyABNmsAddxAR9DNLl8J//8HIkXDPPXZW5DPP2M8oypSxwfDbb+0U7+tx6hQsXWrPV6SIzcWnTtmmRNu324HcquXjMOO+gnr1bHju1csG6l9/hQ8/1Ah+elOypL1mt2yBRx7xnTnl4rd8upvu+vXrKVKkCFmzZvVyVZIeHDt2jO3bt1O6dGm3SxGRdGbjRqhYOZFTWY/QrWwkPRcM58PbmzO0biv6NqlAo7CQFJ87ZVkM3SatJC4+4cyxoMCAyz5PfNipUzB6tO0MFBNjm1d17GgDa9asOA6sWQPz5sHcubBggR3JNMau56xXD+rWtUF1/34bUs+9JXds3z7bGRcgKAiaN7dhtHp1MDjwxx82DY8bZx9YsiS0agUtW8INN7j5ryW+YOxYu4C4dWsYMcJejCKXkFI3XZ8Oo4cOHWLXrl2EhIQQFBSE0YUuHuA4DnFxccTExFCwYEFy6lNdEfGiuDg7+LlybTyN633AuBmvMK90Ndo2fg3HZCAkOIhFXeuk+Pwa/eYTExt30fHLPU/8wLFjdrSxb18bADNlgttvt0mzbl07fzYwkPh4mxWTwumvv0J8fMqnzZUL8uSBvHkvvhUrZgdlg4OBXbvgyy/h889t+g0KsqNfTz0FtWrZEVyRJG+8Yad3DBxo55WLXIJfhlGwgXT37t3EX+pVVuQqBQYGUqBAAQVREfEqx7HbY4waBbfe9x2zFzzB/qy5aPTEuxzNbGcBGWBTvwYpnqNk1xkk95v7cs8TP3L0qF1TOm+evS1fbi+eHDngrrvOhtPy5SFDhjMP37794rCZOzcEBl7ie8XHw/ff2wA6Y4Ydpb3tNnuhNmtmk6xIchIT7cLiiRPtWtL77nO7IvFhKYXRjG4UczVy5sypwCAiImnCiBE2iPbuFkfT0V0ITDjFs01ePxNEAYoEB13yHEWCg5IdGb3c88SPZMsG995rb2Dn1P7449lwOmOGPZ4/P9SpQ7a6dbm3bl2IKGmDbGwsHDgAaw+c/XtKf65dC3v2QMGCtpXvU0/Zdr4il5Mhg53K/e+/NpT+/rttbiVyFXx+ZFRERCQt+OMPuwTw7toOM/M/iflqLM8368n3Jc5+UHwlaz+1ZlTYuvVsMJ03D3bssMcDAux+QZeSM6cdLg0Otn8WKWKDxL33XmYIVSQF//1n95sNDraBNHdutysSH3Rd03SNMXcBXYAqQBHgKcdxRp1zvwF6As8CuYHfgfaO46y+3LkVRkVEJK3buxduvdUOJKx+9n2yde8EvXsz5YGniZq1ju2xcRQJDiIyIvSKAuWUZTHX9DxJgxzHjm7Om2fn6Z4bNJP+TPp7rlw2sIp42i+/QJ06ULs2zJxp9yUVOcf1htH7gZrAn8AY4PkLwuirwOtAK2Ad8Mbpx4c6jnP4UudWGBURkbQsIcEupfr5Z1jx/gJC29eDBx+066zUEEZE0orPPrN7EHXuDO+953Y14mOua82o4zgzgZmnTzTqghMboBPQz3GciaePtQR2A48Bn1xP4SIiIv7szTdhzhz4qt9/hL7e1O6/MXq0gqiIpC2tW8PKlTBoEFSoYNcfi1yGJ8bQSwKFgNlJBxzHiTPG/AzcgcKoiIikU9Omwdtvw3NPxvHohCZ2g/gpU+y6PREXaaq3pIqBA2H1amjbFkJD7T5WIpfgiTBa6PSfuy44vgtI9lXNGPMsdn0pxYsX90AJIiIivmXDBnjiCbg1zOGjhLawdCl89519gybiogubYMXExtFt0kqAywZShVi5pIwZYfx4qF7dbmD7xx92M1uRFLgyR8hxnOGO44Q7jhOeP39+N0oQERFJNceOwcMP25m4sx/6iIxjx9j5ug8+6HZpIkTNWndeN2aAuPgEomatu+TzkkJsTGwcDmdD7JRlMalYrfidPHnsB2/HjkHDhvZPkRR4IozuPP1nwQuOFzznPhERkXTBcaBdO/jrL/i+60/kfbszPPQQ9OjhdmkiAGxPZp/aSx1Pcq0hVtKhsmVh3DhYvtyuHXV5K0nxXZ4Io5uwobN+0gFjTBbgTmCxB84vIiLiNz75BMaMgXc7baX6wKZQujR88YUaFonPKBIcdFXHk1xriJV0qkED6NcPvvkG+vRxuxrxUVf0m9EYk90YU9kYU/n0c4qf/rq4Y/eGGQy8aoxpYowpD4wCjgBfpU7ZIiIiqWfKshhq9JtPya4zqNFv/hVPQ1yyBF58ER665zidFjaB48fVsEh8TmREKEGB5+83GhQYQGTEpdczX2uIlXQsMhJatIDXX7evhSIXuNKPacOBZadvQUCv03/vffr+AcAgYAgQDRQG7rncHqMiIiK+5lrXxcXEQOPGUKSwwzd522Gio+HLL+Hmm71TuKRbV/vhSaOwEPo2qUBIcBAGCAkOom+TCpdtRHStIVbSMWNgxAioVs2G0pUr3a5IfIxxXJ7DHR4e7kRHR7tag4iISJIa/eYTk8y0w5DgIBZ1rZPsc44dg0pVT7JxfQZevTWSPr8NZu2znbn5E238Lqnrws64YAPilYTLa/1+6qYrV237dggPh8yZbYfdfPncrki8zBiz1HGc8AuPe2JrFxERkTTjatfFJSZC/UbHWL8miAdqjaDXwg+ZU7oaL+avR59lMXqjLlflasPepZoKpca11ygsRNe0XL0iRew03bvugkcegTlzIDDQ7arEB6ibgoiIyDmudl1c796weE5WKlefzaglkfwXXIiXHniZY6ccdRmVq3ItU8TVVEj8RrVq8Omn8NNP0LGj29WIj1AYFREROcfVrIv7+mvo1QuKllvJhH+fIYPj0PrhNzicORugQCBX51q2TlFTIfErLVrAK6/Axx/DsGFuVyM+QGFURETkHFfa3GXJErt9Xu07TjIl7kmKHdzJs026sznP2ccpEMjVuJZRTjUVEr/Tp4/d9qVjR/jxR7erEZdpzaiIyDVSI4+063Lr4rZtg0aNoFBBh++LP0eWxct5pWEkfxQrf+YxCgRytYoEByXbPOtSH2okXad6LRK/ERAAX30Ft90GTZvaT/ZuvNHtqsQlCqMiItfgwg6WSWu7AL0JTOOOHoWGDeHIEVj2VF+yfDAKevbkjoZtWKRAINchMiI02c64l/tQQ02FxO/kzAnffWfXkd5/v11HWrCg21WJC7S1i4jINbiW7T/E/yUmQrNmMHkyREeOJ6z//+Dxx+GLL+x+eiLXSTMuJF1ZuBAiIqB0aViwAPLkcbsiSSXa2kVExIPUwTJ96tkTJk6EcR1/JWxwS6hZEz77TEFUPEajnJKu3HknTJ0KDzxgQ+ncuZArl9tViRepgZGIyDXw9Q6WU5bFUKPffEp2nUGNfvMvuTWEXJmvvoK334auzTbSfFxDKFrUDpFmzux2aSIi/qt+ffj2W1i+3DY2OnrU7YrEixRGRUSugS93sLyWvQrl0n7/HZ5+GhrUiKXPigaYU6dg5kzIl8/t0kRE/N+DD9pP/H791S7KP37c7YrESxRGRUSuwZVu/+GGa9mrUFK2dat9b3RDkXgmBTyC2bjBjojedJPbpYmIpB1Nm8LIkTBvHjzyCJw86XZF4gVaMyoico18dW2X1rN6zpEj9gP7uGMOf9dsR6aJ82D0aKhVy+3SRETSniefhLg4aNsWHnsMvv4aMiqupGUaGRURSWN8fT2rv0hMhCeegJUrYckjA8g98TN4/XX7ZklERFLHc8/BoEG2W1yrVpCQcNmniP9SGBURSWN8eT2rP+neHaZMgWktvyV0ZFf43/+gd2+3yxIRSfs6dYJ33oGxY6FdO3B5K0pJPRr3FhFJY5KmDmuvwmvjODZz9usH/Rr/zn3jnoA77rBrmbSFi4iId7z2mu2s26cPBAXB4MF6DU6DFEZFRNIgX13P6uscB3r0sB/Id3lkM6/8/BCmSBE7RJoli9vliYikL2+/DceO2SCaNasNpgqkaYrCqIiICDaIdu0KAwZAlxY7GLD0fszJkzBjBuTP73Z5IiLpjzHw3nu2qVG/fpAtm127L2mGwqiIiKR7jgMvvWQ/fO/+5Fbe+rUuZvt2u5fozTe7XZ6ISPplDAwdakdIe/SwI6QvveR2VeIhCqMiIpKuJSZCx44wZAi81WoD3X+sizlwAGbPtmtFRUTEXRkywOef2xHSl1+2a0jbtXO7KvEAhVEREUm3EhPt+5nhw2HA02vp8kNdzPHjMH8+VKnidnkiIpIkY0bbXff4cXj+eTtC2rKl21XJddLWLiIiki4lJECbNjaIftB6BV2m3YVJSICfflIQFRHxRZkywYQJUK8ePP00jB/vdkVynRRGRUQk3UlIgKeesrO+hj+zhA6T7sZkzgw//wzly7tdnoiIpCRLFtvh/I47oEULmDrV7YrkOiiMiojIdZmyLIYa/eZTsusMavSbz5RlMW6XdEmnTtn3L198AaOe+YU24+thgoNh4UK46Sa3yxMRkcvJls12Og8Lg4cftptDnzrldlVyDRRGRUTkmk1ZFkO3SSuJiY3DAWJi4+g2aaXPBtL4eHj0Ufj6a/j6mbm0/CoCihSxQbRECbfLExGRK5UzJ8yZA//7H/TsCbVqwaZNblclV0lhVERErlnUrHXExSecdywuPoGoWetcqihlJ09Cs2bw7bcw+ZkZNP/iAShVyq4RDQlxuzwREblauXLBl1/axkarVkGlSjBmjN2vS/yCwqiIiFyz7bFxV3XcLSdO2JlcU6bAzKe/pdGoRlChAvz4IxQs6HZ5IiJyPR57DFasgMqVbYfdRx+FAwfcrkqugMKoiIhcsyLBQVd13A0HD0KjRjB9Osxr9QX3jWoO1avD3LmQN6/b5YmIiCeUKGE/YHznHZg40Y6SLljgdlVyGQqjIiJyzSIjQgkKDDjvWFBgAJERoS5VdL7Zs21z3NmzYeETw6kzuiXUrg2zZtnpXSIiknYEBMBrr8Hixbbrbp060LWrXachPimj2wWIiKR1p07BypX2d+Ovv8KSJZCYCAUKnL0VLHj+10m3PHns71Zf1SjMrrWMmrWO7bFxFAkOIjIi9Mxxtxw+DF262D1Ey5aFTS8OpvigznD//XbRaJDvjNyKiIiHVa0Kf/4JL70E/fvbRkdjx8LNN7tdmVzAOC4v8A0PD3eio6NdrUFExJP27YPffjs/fB49au8rVAhuu81+YLt799nb3r02oF4oQwbIl88G0zJl4J57ICICSpb07s/kT+bNg9atYetWeOP5vby+6wUCJnxtF41+9ZXdNF1ERNKHKVPgmWfg2DF47z147jkwxu2q0h1jzFLHccIvOq4wKiJy7RITYc0aGzqTwue6041kAwJsL4Xbb7d7c99+O9xwQ/K/AxMSbK+FcwPqubddu+yHvP/9Zx9fpowNpffcA3ffDdmze+1H9llHjsCrr8LQoXa70O+e/JbQ95+H2Fjo0QO6dYOMmhAkIpLubN8OTz1l12w8+CB89hnkz+92VemKwqiIiAedOgXDhkGvXnYkFOwI5rnBMzzc7st9HseBbdvshpfBwXbd4hXOw3Uc+Ocfu9xx1izbl+HYMQgMhBo1zobTypXtiGp68tNP9n3G5s3Q47nd9NjdgYyTJkCVKjBypO2cKyIi6VdiInzwgV1DGhxsfzfcd5/bVaUbCqMiIh7y22/Qrh0sXw7160OLFjaAliqVzKhnYiKsXg0LF569xcSc/5icOe0vxnNvuXNf/HWJErYLbJYsgN2u5Jdf7Ae9s2bZrvZgp/TWr382nKblnUuOHrW9Kj74AErd6DD9yW+4+aMOcOgQvPkmREZqNFRERM5audJuBbNqFbzwgl1T6u99BE6cgMyZ3a7ikhRGRUSu0759dqbniBEQEgKDB9tliOcF0JMnYenSs8Fz0aKze50VLgx33gk1a9oAGht79nbgQPJfHzp0fhGZM9tFp7Vq2a6wt9125pfozp02mCbd9uyxtd1xBzRpAo0bp621pr/8YkdD16+H7s/s4s3dz5Pxu0m2ccXIkXDLLW6XKCIivuj4cTtC+v77tsvdc89Bw4b2Q19/cPKk/WR89mzbnOnAATt1yocpjIqIXKPERBg9Gl55xb7ed+oEPXtCjhzYhYq//XY2fP72G8TF2SeWKWPDZ9LtxhuZsnz71XWeTUiwgfTAAbs4dcECe1u2zBaWKZMdLa1d+2w4zZqVxEQ7cjt9OkyebP8Odgpvkyb2Vq6cf/ZwiIuD11+HQYOgxA0O01t8TbmhHewwae/etnviOaOhU5bF+Fy3XxER8QGzZsHLL9sZTAAVK9pQ2rAh3Hqr7/ySdBzbkGLOHBtAf/zR/s4LCLC/9+vXh+7dfXomkMKoiMg1+OsvOyV38WK7LnPYsNPLD+PioF8/O73nxAm7SLNSpbPBs2ZN2zr3HFOWxdBt0kri4hPOHAsKDKBvkwpXH44OHrRDg0nh9M8/bTgNDDwbTmvVssOiWbOycaMNpZMm2Z8FbJOfpGAaHu47v3NTcuAAfPcd9OljPwDu2nIHb+1tR8YZU+3PPHKk/YT7HB79NxcRkbRp/XqYOtXeFi2yv0+LFoWHHrLBtHZt73di37fPtodPmu60das9Xrq0XYNTv77tYOgne2YrjIqIXIXDh+3o5wcf2OWaUVHw5JOnGwP98AN06AAbNkDz5tCqlQ19OXNe8pw1+s0nJjbuouMhwUEs6lrn+go+eND+Al2wwHbzWbrUjqpmymR/iTZoYG+lSrF9u/19O2mS/XA1IQGKFbPTeJs0sTnaV/Y23b/fduX/9luYO9f2fSpdymHyI2MpP7yj/VDgrbegc+dki07Vf3MREUl79uyBGTPsL8pZs+zvmZw57T7VDRvapkepEQBPnrSfFieNfi5dakdEg4Ohbl0bPuvXhxtv9Pz39gKFURGRK+A4MGGCzTY7dthlJO+8A3nyYBsPdepkk9FNN8GQIVCv3hWfu2TXGST3imuATf0aeOgnOO3QIRtO5861v1ST9pu56aazwfTOO9l/JBPTptlgOnu2XUaTL58dBb7tNjvgGB5+ekqyl+zdezaAzptnOxeXKAFNH3F4suJybpnQEzNtmv0A4PPPITQ0xXN59d9cRETSlrg4+3t06lSYNs3utRYYaD/kbdjQjpwWK5b8cxMT7VTaw4ftkp7Dh5O/HTpkZzctWHB26u3tt58d/QwP9+npt1dKYVRE5DL++ccOeM6ZY5eKDBsG1aph09AHH9ih0lOn7LqMyMir7lzn6ijdhg02lM6YYX/hnTxpE2b9+jaY3n8/R7IX4ocf7FTYX3+1s5bAjgaXK2eDaVJALVfOs6One/bYacQTJpwdrb3xxrMBtOzKbzDfTrA/R1CQ/YSgY8fLFqGRURER8YiEBPj997PTeZM+5K1c2Y5eXhgyjx61n3BfiTJlzp96e5mZVv5IYVRE5BJGjoS2bc/mnLZtT+ecxYvtotG//rJTdD788JqnyPjM+sUjR+yQ44wZMHPm2a1mqlSxP+P998Ott7LvcCaWLLG/e5NuSY2Bs2e3H9aeG1Dz5rVZPaVbfPzFx5YtswF0wQL7IXLp0mcDaOiKcwJoQICdptS0qZ1PnDfvFf2oPvNvLiIiacvatWen8p46ZT/gvfCWPXvyx8+9ZcuWLjYHVxgVEUnByJHQurXNOl98cbrv0L598Oqr8NlntonB++/bEHSdXX58rrOr49igfXrU1PntN0xiIiczZGRToRJkrhZOiXo1ISwMp2Il/t2e7bxwuny5/R18PUJDbQB9osJyyixLIYA2asSUrSeu6d/O5/7NRURE0hmFURGRZIwebfeqrF/ffsCZJVMijBpl93GJjbWLR3v2tJ9upnFTlsXQ/8tfuHX9MirsWs8tOzdQfvdGcsed3uvUGJscw8LsPOawMI6XDWPZljwsWWJnJGXMeBW3AIfQuOXcGJ1yACVfvjO1aYRTRETEPymMiohc4MsvbYfcunXtOsmgf/+C55+3jX/O28clfUh2faXjEGaOMPmOrLbBwrJl9pbUYh7ghhtsQM2Z03ZAOn7cbndzJX9PTEwxgF62NrT2U0RExB+kFEb9vzWTiMg1GDsWWra0fQKmToWgkUNtQ5zgYNuhtWXLdLGG41zbkwl7GMNycsBDDWzXwCR7954NpsuW2fm6x49Dliz2ljmz/TNXrouPZcnCv7EnWbTtKGuy5GFVldo826TaJUc4k63tEsdFRETE9ymMiki6M26cHRGtVct2av+39ztU6v8680pVZeCj3XiucjUapbMgClAkOCjZ0cciwUEXPzhfvrN7nl2lM1Nui5yechsP3SatBEgxkF5VbSIiIuIX0t+7LRFJ18aPhxYt4M47zw+ic0pXp23j1/g7PhPdJq1kyrIYt0v1usiIUIICz98qJSgwgMiIlPfxvBZRs9adt/YTIC4+gahZ61yvTURERLxHYVRE0o0JE+Dxx+1y0BkzINuoIWeC6PONuhIfEAhcPhilVY3CQujbpAIhwUEY7HrM1GgQdC1Tbr1Vm4iIiHiPpumKSLowcSI8+ijcfrvdWjPbqCHQocNFQTRJel2L2CgsJNUD3rVOufVGbSIiIuI9GhkVkTRv8mT43/+genUbRLOPtkGUhx7i7Za9LgqioLWIqUlTbkVERAQ8GEaNMTmMMYONMVuMMXHGmMXGmKqeOr+IyLWYOhWaNYOqVeH77yHHmLNBlAkT6NygvIKRl2nKrYiIiIBnp+l+ClQEWgLbgBbAXGNMOcdx0l8nEBFx3Xff2a0rq1SxQTTnF+cHUTJlOhOAomatY3tsHEWCg4iMCFUwSmWacisiIiLGcZzrP4kxQcBh4GHHcaaec3wp8L3jOK+n9Nzw8HAnOjr6umsQEQG7bUjUrHWsj87BnilVKF32FH/8kolcX14cREVEREQk9RljljqOE37hcU+NjGYEAoDjFxyPA2p66HuIiFxS0v6V+9flYfeUW8mU/xDcG83mfkuo1K+7gqiIiIiID/FIGHUc57Ax5lfgdWPMKmAn8ChwO7D+wscbY54FngUoXry4J0oQESFq1joO783Enu/CCMx7hALNf6fZiqlUmjNMQVRERETEx3iym+4TQCJ2vegJoCMw7vSx8ziOM9xxnHDHccLz58/vwRJEJD2L2XecPd+FgQP5Gy+l5ZqpvD1nGHNKV1cQFREREfExHgujjuNscBynFpAdKOY4TjUgENjoqe8hInIpp/4ox8ntucl770qe2vjtmSD6dsteCqIiIiIiPsbj+4w6jnPUcZwdxpjcQAQw9XLPERG5XnPmwPafbyBX2FaePjnqTBB9uWl3Ojco73Z5IiIiInIBj23tYoyJwIbbtUBpIOr030d66nuIiCRn1y544gkoW9YwqP0q7n5uGAtKVuHtlr3o3aC8thARERER8UGe3Gc0F9AXKArsByYC3R3Hiffg9xAROU9iIjz5JBw8CPMmHeSWJzpASBFqR8/mpzx53C5PRERERFLgsTDqOM43wDeeOp+IyJV4912YPRuGDXW4ZXAb2LIFFi4EBVERERERn+bJkVEREa/6/Xd47TV4+GF4jk9sx9z+/eH2290uzaOmLIshatY6tsfGUSQ4iMiIUE09FhEREb+nMCoifungQXj0UQgJgc9fXIGp3wnuuw+6dHG7NI+asiyGbpNWEhefAEBMbBzdJq0EUCAVERERv+bxbroiIqnNceDZZ+G//2D8Z0fI+UwzyJsXRo+GDGnrZS1q1rozQTRJXHwCUbPWuVSRiIiIiGdoZFRE/M5nn8E330Cfdxyqj2oH69fD/PmQP7/bpXnc9ti4qzouIiIi4i/S1hCCiKR5a9ZAx45Qty68WnAUfPklvPkm1Krldmmpokhw0FUdFxEREfEXCqMi4jfi4qB5c8ieHcb1WEOGF9pDnTq2i1EaFRkRSlBgwHnHggIDiIwIdakiEREREc/QNF0R8RsvvQSrVsGsycfI374Z5MgBY8dCQMDln+ynkpoUqZuuiIiIpDUKoyLiFyZOhI8/ts1y75nxop2vO3s2FCrkdmmprlFYiMKniIiIpDkKoyLi87ZsgWeegapVoW+Fr6Dlp9C9O9Sr53ZpIiIiInKNFEZFxKfFx9v9RBMSYGLff8jY6DmoWdM2LRIRERERv6UGRiLis6Ysi6HkPZv59Ve44d6l5GjfBDJnhnHjIKM+SxMRERHxZwqjIuKTpiyLofOQLcT8dAPZKmyly4aeBK9bza9vvAdFi7pdnoiIiIhcJ4VREfFJA75fx/YZ5ciQ9QSPFxtOyz9n8Em1JnQ5piAqIiIikhZonpuI+KR1Cwpwcmcwt9adTtS8QSwrHMrAu57gVGyc26WJiIiIiAdoZFREfM6OHXDw55vJVnwHI9Z0AeCFhq8QHxBIkeAgl6sTEREREU9QGBURn/PSS5DBCeDVYn24dcc6ut/zPNtyFSQoMIDIiFC3yxMRERERD1AYFRGfMmcOfP019G+/jW5/fs7i0OpML3sXIcFB9G1SgUZhIW6XKCIiIiIeoDWjIuIzjh+H55+HMqUdOv7bgQADd/zwNZtKlHC7NBERERHxMI2MiojP6NcP1q+HCY9PIWD6d9CrFyiIioiIiKRJCqMi4hP++Qf69oXWTQ9R6dMXoFIl6NTJ7bJEREREJJVomq6IuM5x7PTcoCB4P+frsH07TJoEGfUSJSIiIpJWaWRURFw3bhzMmwefPvcH2T7/CNq3h2rV3C5LRERERFKRcRzH1QLCw8Od6OhoV2sQEfccOAA33ww3Fj/F4lNVMbt3w99/Q86cbpcmIiIiIh5gjFnqOE74hcc1MioirureHfbuhQk138csXw4ffqggKiIiIpIOKIyKiGt+/x0+/hh6ttxM0eFvwEMPQePGbpclIiIiIl6gMCoirjh1Ctq2hcKFHLrFtAdj7KioMW6XJiIiIiJeoDAqIq746CNYvhy+fXQigbNnwltvQfHibpclIiIiIl6iBkYi4nXbtkHZshBx20EmrC6LKVzYztnVVi4iIiIiaU5KDYz0zk9EvK5zZztN97NCr2Hm74Jp0xRERURERNIZTdMVEa+aORO+/RY+bvkrucYOg44doUoVt8sSERERES/TNF0R8Zpjx6B8ecieOZ7lGauQIfYArFkDOXK4XZqIiIiIpBJN0xUR173zDmzaBBuefY8Mw1fC1KkKoiIiIiLplKbpiohXrF0LUVHwcuON3PhFL7uf6EMPuV2WiIiIiLhEYVREvOLxZ46TaOJ5YGlzjiYafni2m9sliYiIiIiLFEZFJNX1+Ggvfy7KQqvS71P7v2gG3PkEnRftY8qyGLdLExERERGXKIyKSKqKj4d3eweRN1cMfTa/w/LCZfgi7H7i4hOImrXO7fJERERExCUKoyKSqoYNg7g92XgnXxfyxh2ie0QHEjMEALA9Ns7l6kRERETELQqjIpJq9u2DN9+EasUX03rjN4yrFMHqgqXO3F8kOMi94kRERETEVQqjIpJqevaEg7EOE3J142jmrLx7Z4sz9wUFBhAZEepidSIiIiLiJoVREUkVq1fDxx/D0IipFF/5M1s6vkLWIoUwQEhwEH2bVKBRWIjbZYqIiIiIS4zjOK4WEB4e7kRHR7tag4h4luNARAT8teQ4McHlCMieFZYvh4wZ3S5NRERERLzMGLPUcZzwC4/rnaGIeNz06TBnDiy6/10CZm6CuXMVREVERETkPJqmKyIedfIkvPwy1Cq1jdsX9IEmTaBuXbfLEhEREREfozAqIh710Ufw77/wVbFXMQkJMHCg2yWJiIiIiA9SGBURj9mzB3r3hper/0KRBV/BK69AyZJulyUiIiIiPkhhVEQ8pkcPOHY4gbcPdYSiReHVV90uSURERER8lDqKiIhH/PUXjBgBY2t/Tpb5y2DcOMiWze2yRERERMRHaWRURK6b40CnTlA8ZyzNVrwGd94JzZu7XZaIiIiI+DCFURG5blOmwI8/wpTKb5LhwH744AMwxu2yRERERMSHeSSMGmMCjDFvGWM2GWOOn/7zbWOMpgGLpHEnTkCXLvBgqTVUXPgRtGkDlSu7XZaIiIiI+DhPhcVXgfZAS2AlUBEYDZwA3vLQ9xARH/T++7Bxo8OSsBcxOXLA22+7XZKIiIiI+AFPhdE7gGmO40w7/fVmY8x3QHUPnV9EfNCuXTZ7vh0+lbzRc+303Hz53C5LRERERPyAp9aM/gLcbYy5GcAYUw6oA8z00PlFxAd17w5O3HFe2fkS3HILtGvndkkiIiIi4ic8NTLaH8gBrDHGJJw+7zuO4wxN7sHGmGeBZwGKFy/uoRJExJuWLYPPP4fpNd4j8JdNMGouZNQycRERERG5Mp4aGW0OPAk8Btx6+u/PG2NaJ/dgx3GGO44T7jhOeP78+T1Ugoh4i+PAiy9ChdzbuO/Pd6BJE6hb1+2yRERERMSPeGoYIwoY6DjO16e/XmmMuQHoBnzmoe8hIj7i229h4UL4p9qrmBUJMHCg2yWJiIiIiJ/x1MhoViDhgmMJHjy/iPiIY8fsVi4tSy+izJKvIDISSpZ0uywRERER8TOeGhmdBnQ1xmwCVgNhwEvAGA+dX0R8RFQUbPsvgQ/LvABFi0LXrm6XJCIiIiJ+yFNh9AXsfqJDgQLADmAE0NtD5xcRH7BlC/TrB8PCPydH9DIYNw6yZXO7LBERERHxQ8ZxHFcLCA8Pd6Kjo12tQUSuTLNmsGjafv4LuomA8uXgp5/AGLfLEhEREREfZoxZ6jhO+IXHtaZTRK7Ijz/ChAkwpeIbBBw8AB99pCAqIiIiItdMYVRELuvUKbuVy32FlxMePQzat4eKFd0uS0RERET8mHaoF5HL+uQTWLnSYX7ZFzDxeaBXL7dLEhERERE/pzAqIpe0bx/06AHvlPuKfGt+gU8/hdy53S5LRERERPycwqiIXFKPHuAcPERkQBeoWhWeesrtkkREREQkDVAYFZEUrVhhp+jOrvQWgct3wczvIIOWmouIiIjI9dO7ShFJluPYpkXVcvxNnZWD4emn7cioiIiIiIgHaGRURJI1YQL89JPD1ps7YnZmh7593S5JRERERNIQhVERucixY9ClC7xUYjJF186FDz+E/PndLktERERE0hBN0xWRi/TvD3u3HuOduM5QoQK0bet2SSIiIiKSxmhkVETOs3kzDBgAX97Snyyr/4NvfoKMeqkQEREREc/SyKiInKdLF7iRjTRZ3x8eewzuusvtkkREREQkDdJwh4icMX8+TJwIa2/ujNkWCFFRbpckIiIiImmUwqiIAHDqlN3KpVXB7wld+51dOFqkyEWPm7IshqhZ69geG0eR4CAiI0JpFBbiQsUiIiIi4s8URkUEgI8/hn9WneDXwi/CTTdBp04XPWbKshi6TVpJXHwCADGxcXSbtBJAgVRERERErorWjIoIe/fCG2/AkNKDyL7jX7uVS6ZMFz0uata6M0E0SVx8AlGz1nmrVBERERFJIzQyKiL06AE5Dm7j6eNvQaNGcM89yT5ue2zcVR0XEREREUmJRkZF0rnly2H4cJhUOpIMTiK8916Kjy0SHHRVx0VEREREUqIwKpKOJSZC+/bwQPYFVPnna+jaFUqWTPHxkRGhBAUGnHcsKDCAyIjQ1C5VRERERNIYTdMVScfadovl98XZ+StHG7bnLkT0PY/z0CUen9SkSN10RUREROR6KYyKpFPDZ+zk00H5eClPP8rtX8+zjbuzcOZ6ErMEXTJcNgoLUfgUERERkeumaboi6ZDjwKudA7mBzfQ69A4LSlZhdpnb1BlXRERERLxGI6Mi6dAXX0Dsv3mYlOdBnCOG7hHtwRhAnXFFRERExDs0MiqSzuzaBZ06Qbu8Q7h7/6/0q/0UMbkKnLlfnXFFRERExBs0MiqSzrzwAuQ6EsOgwO4sKV6BsZXvPXOfOuOKiIiIiLdoZFQkHZk8GSZMcJh9Y1syO/Ec/GAoRXJnwwAhwUH0bVJBzYlERERExCs0MiqSTsTG2j1FuxYfR5l10+G996jfsCb1G7pdmYiIiIikRwqjIulEZCSwaxdv5XgBbrsNOnZ0uyQRERERScc0TVckHZg3Dz79FH4o8wIZ447A559DQIDbZYmIiIhIOqaRUZE07uhRePZZaF9oIhXXTYA+faBsWbfLEhEREZF0TiOjImncG29A7MZ9vHeyPdx6K3Tp4nZJIiIiIiIaGRVJy37/HQYPhkU3dSbTxn3w+SwIDHS7LBERERERjYyKpFUnT0Lr1tAi9wxu++cLeO01qFTJ7bJERERERACNjIqkWX37wtbVB/kj73NQvjx07+52SSIiIiIiZyiMiqRBq1bBO+/AD6UiCdq0A76fDJkyuV2WiIiIiMgZmqYrksYkJNjpuQ8GzaXOhhHw8stQtarbZYmIiIiInEcjoyJpzAcfwOolR5iXvw2UKQO9erldkoiIiIjIRRRGRdKQjRvt0tDxJV4j25YtMOlnCApyuywRERERkYtomq5IGpGQAG3awJ0s5MHNH2I6dICaNd0uS0REREQkWRoZFUkj3nwTFs+PY2fB1lCwBPTp43ZJIiIiIiIpUhgVSQOmToW334bvy79BrlX/wty5kD2722WJiIiIiKRI03RF/Ny6dfDEE9Cx9EwiVr9r5+rWret2WSIiIiIil6QwKuLHDh+Gxo3h5ozrGbT7MUylSjB4sNtliYiIiIhclqbpivgpx4GnnoJta4+w/YZGZDgUAJMnQ9asbpcmIiIiInJZCqMifioqCiZOdFhb8Smyr/qbRUPG8srXG9keu5oiwUFERoTSKCzE7TJFRERERJKlMCrih+bMgW7d4MuKUYT+9S2rXuzOMzHBxMXHARATG0e3SSsBFEhFRERExCdpzaiIn9m8GR59FFoXm81jq7pBs2Y8V/Bu4uITzntcXHwCUbPWuVOkiIiIiMhlKIyK+JG4OGjSBEJObmRY7P8wt9wCn3/O9oPHk3389tg4L1coIiIiInJlFEZF/ITjQNu2sG7ZURbmbUyAcWzDomzZKBIclOxzUjouIiIiIuI2j4RRY8xmY4yTzG2GJ84vIjBsGIwZ4/Bb+Tbk3LISxo2DUqUAiIwIJSgw4LzHBwUGEBkR6kapIiIiIiKX5akGRlWBc98JFwaWAt946Pwi6cqUZTFEzVrH9tg4igQH8UDBCrzxYn5GlBtMhVXj4J134N57zzw+qUnRuc9RN10RERER8WXGcRzPn9SY7kAkUNhxnEsuWgsPD3eio6M9XoOIp10YEFMr7E1ZFkO3SSvPNCQ6dSQzO0fXpFGWhUw4cB+mYUP49lswxuPfW0RERETE04wxSx3HCb/wuMe3djHGGKA18OXlgqiIv7gwIKbm1ilRs9ad+T5OgmHvlFspejyG4QnNMDfdBKNGKYiKiIiIiN9LjQZG9YGSwIiUHmCMedYYE22Mid6zZ08qlCDiWecGxCSptXXKuR1wD8wvh4kJYmq2B8iYcBKmTIEcOTz+PUVEREREvC01wmgb4A/HcVak9ADHcYY7jhPuOE54/vz5U6EEEc9KaYuU1Ng6JakD7pGVIRz+8wY+y/s4lQ6u5a3/vQY33eTx7yciIiIi4gaPhlFjTAGgIZcYFRXxR97cOiUyIpTEbfnZP7sCL+buz2P7JvPhXS2o0amVx7+XiIiIiIhbPD0y2go4AYzz8HlFXOXNrVMSN4Ww45uq1Msxh6iDr/NL2dsp9l4fdcYVERERkTTFYw2MTjcuegb42nGcI546r4gv8NbWKZ9/Dm3aQJMK6xm/vSUZSpei5q/fQ65cHv0+IiIiIiJu82Q33dpAGaCFB88p4jMahYWk6ujke+/Byy/Dc7etYNiGezBOom1YpCAqIiIiImmQx6bpOo7zo+M4xnGcJZ46p0h64Djw+us2iHavvYhhf9fCZM4ECxdC2bJulyciIiIikio8vs+oiFy5xETo0AGGDYNBET/w4s9NMEWLwpw5cMMNbpcnIiIiIpJqFEZFXBIfDy1bwrhx8OVD3/DY9y0wt9wCs2ZBgQJulyciIiIikqoURkVccOwYNG0KM2fC901GcO/k56BGDZg2DYKD3S5PRERERCTVeXprFxG5jNhYiIiA77+HX5tEce+kZ+Hee+2IqIKoiIiIiKQTGhkV8aJdu2zuXL3KYdVDr1FuUj9o3hzGjIFMmdwuT0RERETEaxRGRbxkyxaoXx+2b01g/T0dKD71Y3juORgyBAIC3C5PRERERMSrNE1XxAv+/tsuCT2wO55NNVpQfObH0LWrbaOrICoiIiIi6ZDCqEgqmzMH7rwTAuOPsaFiI/LP+xr694e+fcEYt8sTEREREXGFwqhIKtm61XbMveceKJ7rIGuK30vOX76H4cPhlVfcLk9ERERExFUKoyIeduKEHfS8+WaYPh3e67Kd6Jx3E7TiN/j6a2jTxu0SRURERERcpwZGIh70ww/QsSP8+y80aZjA8MpDyTvodYiPh+++s610RUREREREI6MinrB5MzRuDPfdZ79ePHgJE/+rSt5eHaF6dVixQkFUREREROQcCqMi1+H4cXjrLShbFmbPhvd6HODvWm25vfNtsHMnjB8Ps2ZBmTJulyoiIiIi4lM0TVfkGs2YAS++CBs2wCMPO3xc4wvy9u0C+/bZO3r1gpw53S5TRERERMQnaWRU5Cpt3AgPPggPPACBgbB4xGom7KlN3pdaQqlSsHQpDBqkICoiIiIicgkaGRW5AgkJsHgxfPMNjBgBGTPCoLeP8kLsWwS0exdy5LBbtrRuDRn0GY+IiIiIyOUojIqk4NQp+PlnmDgRJk2yS0AzZ7Z7hw6+e6ptTvTff/DUU9C/P+TP73bJIiIiIiJ+Q2FU/N6UZTFEzVrH9tg4igQHERkRSqOwkGs6V3w8zJ9vA+jkybB3LwQFwf33wyOPwAM3rSP7m5HQehqULw8LF0LNmh7+iURERERE0j6FUfFrU5bF0G3SSuLiEwCIiY2j26SVAFccSE+cgDlzbACdOhUOHIDs2e2a0Ecegfsq7yDrtPHw3lfwxx+QLRtERdkmRYGBqfaziYiIiIikZQqj4teiZq07E0STxMUnEDVrXYphND4e/v0X/voLpk+HadPg0CHIlQsaNoSHH4Z7qsWSZeYkGPYV/PgjJCZCWJgNoY8/DoULe+PHExERERFJsxRGxa9tj41L8fipU3bbldWr7W3VKvvnP//YQAqQJ49dA/rww1D3jjgyzZ0Jo8ZC0xlw8qTtjtu9Ozz6qN1MVEREREREPEJhVPxakeAgtu0/zqmDQcTvzUH8nhyc3JsdDuQk2yCbJ5PceCPccovdluWWW+yt/M2nCPzlR/jqK/jfJDtEWrAgtGsHjz0GVauCMe79gCIiIiIiaZTCqKQKTzYVSrJvH6xbd/7tvxV3sm1LACSe3U4lY644KlUw1Ln9bOgsW9Yu9SQ+HtassXuBfrIEpkyBXbvsnqAPP2wDaO3adu8WERERERFJNXrHLR53PU2FEhPtNNq1ay8Onvv2nX1cYCCULg3hlQKpWuswfx3exrGs+yle6hRdHyptv098vJ2Xu3QpfL7U/rlihe1YBBzNnJWfb6jMwnrtuK1DCx66rVSq/Ft4OpSLiIiIiKQFCqPicVfbVGj/fpg9G2bOhB9+gD17zt5XsCCEhkKTJvbPm2+2f5Yoce7gZQ6IL302eI4Yf1HwJEcOuPVW6NCBP/KW5I2tmVmboyCOsSOqk2f8S2LmLB4Nip7o9CsiIiIiklYpjIrHXaqpENjRz+XLbfj8/nv47Td7LG9eiIiAevXs1NrQUNvh9iKxsbBohQ2by5fbP1etOrtANGfOM8GTKlXsrXRpyGCDZ6d+84nJeX6Nl+vAey2updOviIiIiEh6oTAqHlckOIiYCwJp4vGMBO0pwlNP2QC6a5c9Hh4Or78O999v/x4QcM6THAc2bjo/dC5fDlu2nH1M/vxQuTJ07Hg2eJYqdSZ4JudyYdlTvPV9RERERET8kcKoeFxkRCjdJq3k0J5MHPu7CHEb83MiJjc4GdiT245+3nef/bNgwXOeGBdn5+vOnWtD519/2e62YMPlTTfBbbdB27ZQqZINoYUKXXW32+TCctJxT/LW9xERERER8UcKo+JxdxYPofSGXMz4JitOQgayFTnEg08fofPTOalW7YJGtbGxMH06TJ5sF4weO2bb3laqBC1anA2d5ctD1qweqS8pLJ87hTYoMIDIiFCPnN/b30dERERExB8pjIrHHD8OH30Eb78Nhw9np80zdgpusWI5z3/gjh0wdaoNoPPnw6lTULgwGxs8wgc5bmFGnlAK5M2Rap1nk86Z2l1uvfV9RERERET8kXEcx9UCwsPDnejoaFdrkOvjODB+PHTrBps32/WfAwbYJkRnbNhgw+ekSbZjkeNAmTLQuDE0bsyUwBC6TVl90Shi3yYVFN5ERERERPyYMWap4zjhFx7XyKhcl19+gZdfhiVL7IzaOXNsN1wA1q+HL76wIXSl3dKEsDDo3duG0HLlzqz3jOo3X51nRURERETSEYVRuSb//gtdu9qBzpAQGDXKLvEMCMCuA33rLfjgA7tnS82aMGgQNGpkNwhNhjrPioiIiIikLwqjclX27rUDm8OGQZYsdn1o586newslJMCIz6F7d/vA1q2hVy8oUuSy51XnWRERERGR9CXlzRhFznH8OERFQenSMGSIzZn//mtzZ9aswM8/241Cn30WQkMhOhpGjLiiIAq282xQYMB5x9R5VkREREQk7dLIqFzW+vW2KdG//0KDBrY5Ublyp+/csgUiI2HCBChWDL7+Gpo1u+q9P9V5VkREREQkfVEYlUtavhzurpfAkbgECjZbxoGwo/xzIpRyR4Ohf387XGqMnY7bpct17QXaKCxE4VNEREREJJ1QGJUU/fQT3P9AIicznKTAo0sIzHeEmAMOv/T+gIjFYwjavRMefdSG0mLF3C5XRERERET8iMKoJGvKFPjf/yBDrjgKPvwbGXMep8KOf+k5bzjhMX+zNuQmbl44wXbKFRERERERuUoKo3KRzz+HNm1sP6Idty0iS+Ax3vzhEx5b8QN7sgUTed+LTKxQl40KoiIiIiIico3UTVfOM2CA7ZRbrx7Mmwdl8sbz6cTePLbiB4ZXbczdbYYzoWJ9CufO5napIiIiIiLixzQyKgA4DrzyCgwcaKfnjh4NmY7s55uJPQnesoLI+15kQsX6gLZcERERERGR66cwKpw6Bc88YwNo+/bwwQeQYdcOuOce8v7zD79HfcLi+Bsx2nJFREREREQ8RGE0nYuLs9uCTp9ud2fp0QPMpo1Qvz7s2gUzZ1K9bl0WuV2oiIiIiIikKQqj6VhsLDz4ICxaBEOHQrt2wMqVEBEBJ07A/PlQrZrbZYqIiIiISBqkMJpO7dgB994Lf/8NX39tR0f57Te4/34ICoKFC6FcObfLFBERERGRNErddNOhDRvs9qAbNsCMGaeD6Jw5ULcu5M1rh0oVREVEREREJBVpZDSd2bABatSwTYt+/BGqVgW+/RYee8wG0FmzoGBBt8sUEREREZE0TiOj6cjx49C0qV0OunDh6SD66afQvLldG7pggYKoiIiIiIh4hcfCqDGmsDFmtDFmjzHmuDFmjTGmlqfOL9evc2dYtgzGjIGyZYGoKGjTxjYsmj0bgoPdLlFERERERNIJj0zTNcYEA4uAX4AGwB7gRmC3J84v1++rr+Djj+GVV+DBBxzo9hr062dHRceMgUyZ3C5RRERERETSEU+tGX0F2OE4zpPnHNvkoXPLdfr7b3j2Wdu06O3eidC2HQwfDm3bwkcfQUCA2yWKiIiIiEg646lpuo2A340x440xu40xy40xHYwxxkPnl2t09KhdJxoUZLdwCRzwjg2i3brZzUUVREVERERExAWeCqM3As8DG4EI4H2gH9A+uQcbY541xkQbY6L37NnjoRLkQo4Dzz8Pa9bYaboha+fh9OzJD5XrUTLhDmr0/5Epy2LcLlNERERERNIh4zjO9Z/EmJNAtOM4d5xzrA/Q2HGcspd6bnh4uBMdHX3dNcjFPv8cWreGnj3hzed2cLx8RbZlCOKhJ97jWKYgAIICA+jbpAKNwkJcrlZERERERNIiY8xSx3HCLzzuqZHRHcCaC479DRT30PnlKq1YAe3bQ9260KPbKXj0UZzDR2jbsNuZIAoQF59A1Kx1LlYqIiIiIiLpkafC6CIg9IJjNwFbPHR+uQqHDtl1orlzw9ixENC7J/z0E93veZ71+S7+fGB7bJwLVYqIiIiISHrmqW66g4DFxpjuwHggDOgIvOah88sVchx45hnYuBHmz4eCf34PffrAM8/we6kGkEzwLBIclMyZREREREREUo9HRkYdx/kD21G3GbAKeAfoAQz1xPnlyg0dChMmwDvvwF0lt0KLFlCpEnzwAZERoQQFnt89NygwgMiICwe1RUREREREUpenRkZxHGcGMMNT55Or98cf0LkzNGgAkS+ehLubQXy8TadBQWeaFEXNWsf22DiKBAcRGRGq5kUiIiIiIuJ1Hguj4q4DB+w60cKFYfRoyPBaV/jtN/jmGyhT5szjGoWFKHyKiIiIiIjrFEbTAMeBli1h+3ZYuBDy/jwZBg2CF16wCVVERERERMTHKIymAe++C9OmwfvvQ/V8GyDiKahaFaKi3C5NREREREQkWQqjfm7RIujaFR5+GF5ocxxqNgNj7PTczJndLk9ERERERCRZCqN+bM8eaN4cSpSAzz4D8/JL8OefMHWqPSgiIiIiIuKjFEb9WOfONpD+9hvkmjkOhg2DyEh46CG3SxMREREREbkkj+wzKt73yy8wdiy88gqEBa2FNm2gRg27waiIiIiIiIiP08ioH0pIsI1yixWDbi8eg7ubQlAQfP01BAa6XZ6IiIiIiMhlKYz6oREjYPlyGD8esr7SAVavhh9+gKJF3S5NRERERETkimiarp/Ztw+6d4e774am8V/ByJHw+utwzz1ulyYiIiIiInLFNDLqZ3r0gIMH4aN3DmIavwTVq0PPnm6XJSIiIiIiclUURv3I8uXwySfQoQOUm9ALdu+G6dMhIMDt0kRERERERK6KwqifcBzbtChPHujdfDXc9YHtoBse7nZpIiIiIiIiV01h1E+MG2e3cxkx3CFXj46QM6e2cREREREREb+lMOoHjhyByEg7CPp0zm9h/nwYMgTy5XO7NBERERERkWuiMOoH3n4btm+HyV8eJUPLl6FyZXjuObfLEhERERERuWYKoz7un3/gvfegZUuoNq8vbN0KX32lpkUiIiIiIuLXFEZ9XOfOkCULRD23HmpHQYsWULOm22WJiIiIiIhcF4VRHzZ9OsycCe++C/n7dIZMmWDAALfLEhERERERuW4Z3C5Aknf8OHTqBGXLQscbp9tk+uabULiw26WJiIiIiIhcN42M+qj33oMNG2Du9ONkfLET3Hyz3WhUREREREQkDVAY9UFbt9otRBs3hrorTqfS2bPtNF0REREREZE0QGHUB73yCiQmwvsv/wf134aHH4b69d0uS0RERERExGMURn3MTz/B119Dz55Q7P0u9uC77170uCnLYoiatY7tsXEUCQ4iMiKURmEhXq5WRERERETk2iiM+oCkYBmz/zh7vriL/IWy0LXqT9BrAvTuDTfccNHju01aSVx8AgAxsXF0m7QSQIFURERERET8grrpuiwpWMbExnFoWXGO7cxO0O1/crLj83DjjRAZedFzomatOxNEk8TFJxA1a523yhYREREREbkuGhl1WVKwTDgWyMGFoWS5YS+tY0eQc+M/MHUqZMly0XO2x8Yle66UjouIiIiIiPgajYy6LClAxi4MJfFkADfd8TOdF43lxxurwIMPJvucIsFBV3VcRERERETE1yiMuqxIcBCnDgZx5K9iZK/8Hz1WDiVTQjwfN+4IxiT7nMiIUIICA847FhQYQGREqDdKFhERERERuW4Koy6LjAjlyO9lwDjUKTGVR1bNY3T1Jjz6eN0Un9MoLIS+TSoQEhyEAUKCg+jbpIKaF4mIiIiIiN8wjuO4WkB4eLgTHR3tag1u2rwZSpdxKBi+hVmb6lEg7iC/z1zEgzVucrs0ERERERGR62aMWeo4TviFxzUy6rI+fSAgg2FVo1mU37WBAsM/UhAVEREREZE0T2HURZs3w8iR0P6pY+Qe1BPuuguaNXO7LBERERERkVSnrV1c1KcPZMgAPQt+DLt2wYQJKTYtEhERERERSUs0MuqSTZvsqGiHp46Sa1g/qFcP7rzT7bJERERERES8QmHUJUmjom/kGwp79kCvXm6XJCIiIiIi4jUKoy7YtAlGjYIOrY6Q65MBcM89cMcdbpclIiIiIiLiNQqjLjgzKpp3COzdq1FRERERERFJdxRGvSxpVPSFVofJNTwK7rsPbrvN7bJERERERES8SmHUy5JGRXvk/hD27YM333S7JBEREREREa9TGPWipFHRjq0OkWv4QGjQAKpVc7ssERERERERr1MY9aIzo6K5PoADB7RWVERERERE0i2FUS85Myra8iA5R7wLDz0EVaq4XZaIiIiIiIgrFEa95MyoaI7BEBurtaIiIiIiIpKuKYx6QdKo6IstY8n52SBo1AjCwtwuS0RERERExDUKo17Qpw8EBMDr2QbBwYMaFRURERERkXRPYTSVnRkVffIAOT8fDA8/DJUquV2WiIiIiIiIqxRGU9k779hR0e5Z34NDh6BnT7dLEhERERERcZ3CaCrauBFGj4ZOT+wj52eDoWlTqFDB7bJERERERERc55Ewaox50xjjXHDb6Ylz+7OktaLds7wLR49qVFREREREROS0jB481zqg9jlfJ3jw3H4naVS0S6u95Bj1ITRvDrfc4nZZIiIiIiIiPsGTYfSU4zjpfjQ0SdKo6GuZBtpR0TfecLskERERERERn+HJNaM3GmO2G2M2GWO+Nsbc6MFz+5WkUdGXWuwmx+iP4NFHoWxZt8sSERERERHxGZ4Ko78DrYB7gTZAIWCxMSZvcg82xjxrjIk2xkTv2bPHQyX4jqRR0W4ZoyAuTqOiIiIiIiIiFzCO43j+pMZkBzYC/RzHee9Sjw0PD3eio6M9XoNbtm6FG2+EyCd30WdcSXjkERgzxu2yREREREREXGGMWeo4TviFx1NlaxfHcY4Aq4EyqXF+XzZ4MDgOdM0wAE6ehB493C5JRERERETE56RKGDXGZAFuBnakxvl9VWwsDB8Ozz20g5xfDoUWLaBMusvjIiIiIiIil+WpfUYHGmNqGWNKGmOqA98C2YDRnji/v/j4YzhyBN4M6g/x8RoVFRERERERSYGntnYpCowD8gF7gN+A2xzH2eKh8/u848fh/fehaa3d5J/0CTz5JJQq5XZZIiIiIiIiPskjYdRxnP954jz+7MsvYedO6F/nAzhxArp2dbskERERERERn5Uqa0bTm8REGDgQalQ8TImZQ6BJE7jpJrfLEhERERER8VkKox4wbRqsWwcfVBiBiY2FV15xuyQRERERERGfpjDqAVFRULr4ScIWvAe1a0O1am6XJCIiIiIi4tMURq/TokX2NvTOrzAxMfDqq26XJCIiIiIi4vMURq9TVBTkzZ1I3aVRUKkSRES4XZKIiIiIiIjP89TWLunS2rXw3XfwRbMZZBi/BsaOBWPcLktERERERMTnaWT0Orz7LmTODE039YcSJaBZM7dLEhERERER8QsaGb2EKctiiJq1ju2xcRQJDiIyIpRGYSGA3VN0zBh4+75FZJq6CD78EDLqn1NERERERORKKD2lYMqyGLpNWklcfAIAMbFxdJu0EoBGYSF88AHEx8Pzh/tD3rzw1FNulisiIiIiIuJXNE03BVGz1p0Jokni4hOImrWOw4dh2DDoWG812eZPgxdegGzZXKpURERERETE/yiMpmB7bFyKxz/9FGJjoXtgFGTNCh06eLc4ERERERERP6cwmoIiwUHJHi+UIyuDBsEjt20j/+yx8MwzdpquiIiIiIiIXDGF0RRERoQSFBhw3rGgwACqO5XYuhX6FxoEjgMvvXTeY6Ysi6FGv/mU7DqDGv3mM2VZjDfLFhERERER8QtqYJSCpK6553bT7XJPKD1b5eH2mw9Qcu5wePRRuOGGM8+5XNMjERERERERsRRGL6FRWMh5IXLWLFi5EpY2GYpZewQiI897/KWaHimMioiIiIiInKVpuldhwAC4sXAcYQvfh/vug4oVz7v/Uk2PRERERERE5CyF0Su0dCnMnw8f3zYKs2cPvPrqRY9JqelRSsdFRERERETSK4XRKxQVBcE5Eqi7bCBUrw533XXRY1JqehQZEeqtMkVERERERPyC1oxegY0bYcIEGPPARDJ8txHeGwjGXPS45JoeRUaEar2oiIiIiIjIBRRGr8CgQRCQwaHppv4QGgoNG6b42AubHomIiIiIiMjFFEYvY+9e+Owz6FN3Hplm/QmffgoZNLtZRERERETkeihVXcbQoRAXB20P9YfChaFFC7dLEhERERER8XsKo5dw7Bh8+CG8WHMp2X+dC506QebMbpclIiIiIiLi9xRGL+GPP+DwYegeOABy5oTnnnO7JBERERERkTRBYfQSatWC7Qs3kO+nb6FdO8iVy+2SRERERERE0gQ1MLqMPCPfhYwZ4cUX3S5FREREREQkzdDI6KXs3g0jR0LLlrZ5kYiIiIiIiHiEwuilbN0KN94IXbq4XYmIiIiIiEiaomm6l1KlCqxaBca4XYmIiIiIiEiaopHRy1EQFRERERER8TiFUREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOoVRERERERER8TqFUREREREREfE6hVERERERERHxOuM4jrsFGLMH2OJqEZeXD9jrdhGS7uk6FF+g61B8ha5F8QW6DsUX+MN1eIPjOPkvPOh6GPUHxphox3HC3a5D0jddh+ILdB2Kr9C1KL5A16H4An++DjVNV0RERERERLxOYVRERERERES8TmH0ygx3uwARdB2Kb9B1KL5C16L4Al2H4gv89jrUmlERERERERHxOo2MioiIiIiIiNcpjIqIiIiIiIjXKYxegjHmeWPMJmPMcWPMUmPMnW7XJOmHMeZNY4xzwW2n23VJ2meMucsY850xJub0ddfqgvvN6etzuzEmzhizwBhzi0vlShp1BdfhqGReI39zqVxJo4wx3YwxfxhjDhlj9hhjphljyl/wGL0mSqq6wuvQL18TFUZTYIxpDrwP9AHCgMXA98aY4q4WJunNOqDwObcK7pYj6UR2YBXwIhCXzP2vAC8DLwBVgd3AHGNMDq9VKOnB5a5DgLmc/xp5v3dKk3SkNjAUuAOoA5wC5hpj8pzzGL0mSmqrzeWvQ/DD10Q1MEqBMeZ34C/Hcdqcc+xf4FvHcbq5V5mkF8aYN4FHHMcpf7nHiqQWY8wRoIPjOKNOf22A7cBHjuO8c/pYEPbNVxfHcT5xq1ZJuy68Dk8fGwXkcxznAbfqkvTHGJMdOAg0chxnml4TxQ0XXoenj43CD18TNTKaDGNMJqAKMPuCu2ZjP5EQ8ZYbT0/72WSM+doYc6PbBUm6VxIoxDmvj47jxAE/o9dH8b6axpjdxph/jDEjjDEF3C5I0rwc2PfPB05/rddEccOF12ESv3tNVBhNXj4gANh1wfFd2BccEW/4HWgF3Au0wV57i40xed0sStK9pNdAvT6K234AngTqYqdIVgPmG2Myu1qVpHXvA8uBX09/rddEccOF1yH46WtiRrcLEJHkOY7z/blfn16EvhFoCbznSlEiIj7CcZyvz/lypTFmKbAFaABMcqcqScuMMe8BNYGajuMkuF2PpE8pXYf++pqokdHk7QUSgIIXHC8IqJupuMJxnCPAaqCM27VIupb0GqjXR/EpjuNsB7ah10hJBcaYQcCjQB3HcTaec5deE8VrLnEdXsRfXhMVRpPhOM5JYClQ/4K76mO76op4nTEmC3AzsMPtWiRd24R9g3Xm9fH0tXknen0UFxlj8gEh6DVSPMwY8z5nA8DaC+7Wa6J4xWWuw+Qe7xeviZqmm7L3gC+MMUuARUBboAjwsatVSbphjBkITAP+AwoAPYBswGg365K073SXvtKnv8wAFDfGVAb2O47znzFmMPCaMWYt8A/wOnAE+MqFciWNutR1ePr2JjAR+0arBNAX28F0spdLlTTMGDMEeAJoBBwwxiStAz3iOM4Rx3EcvSZKarvcdXj69fJN/PA1UVu7XIIx5nns3lGFsXuddXYc52d3q5L0whjzNXAXtqHWHuA3oIfjOGtcLUzSPGNMbeDHZO4a7ThOq9NbGfQEngNyY5tttXccZ5XXipQ071LXIdAOmILdBzwY++brR+xr5FavFCjpgjEmpTfKvRzHefP0Y/SaKKnqctfh6e2EpuCHr4kKoyIiIiIiIuJ1WjMqIiIiIiIiXqcwKiIiIiIiIl6nMCoiIiIiIiJepzAqIiIiIiIiXqcwKiIiIiIiIl6nMCoiIiIiIiJepzAqIiIiIiIiXqcwKiIiIiIiIl6nMCoiIiIiIiJe939RDg2irxvX1AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x1, y, \"o\", label=\"Data\")\n",
    "ax.plot(x1, y_true, \"b-\", label=\"True\")\n",
    "ax.plot(np.hstack((x1, x1n)), np.hstack((ypred, ynewpred)), \"r\", label=\"OLS prediction\")\n",
    "ax.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting with Formulas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using formulas can make both estimation and prediction a lot easier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:59.514502Z",
     "iopub.status.busy": "2021-09-17T11:14:59.513584Z",
     "iopub.status.idle": "2021-09-17T11:14:59.521662Z",
     "shell.execute_reply": "2021-09-17T11:14:59.522302Z"
    }
   },
   "outputs": [],
   "source": [
    "from statsmodels.formula.api import ols\n",
    "\n",
    "data = {\"x1\": x1, \"y\": y}\n",
    "\n",
    "res = ols(\"y ~ x1 + np.sin(x1) + I((x1-5)**2)\", data=data).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the `I` to indicate use of the Identity transform. Ie., we do not want any expansion magic from using `**2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:59.525205Z",
     "iopub.status.busy": "2021-09-17T11:14:59.524332Z",
     "iopub.status.idle": "2021-09-17T11:14:59.531320Z",
     "shell.execute_reply": "2021-09-17T11:14:59.531942Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept           4.942116\n",
       "x1                  0.501201\n",
       "np.sin(x1)          0.484040\n",
       "I((x1 - 5) ** 2)   -0.020209\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": {
    "execution": {
     "iopub.execute_input": "2021-09-17T11:14:59.534777Z",
     "iopub.status.busy": "2021-09-17T11:14:59.533884Z",
     "iopub.status.idle": "2021-09-17T11:14:59.541031Z",
     "shell.execute_reply": "2021-09-17T11:14:59.541643Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    10.844029\n",
       "1    10.698807\n",
       "2    10.444329\n",
       "3    10.123852\n",
       "4     9.794320\n",
       "5     9.512418\n",
       "6     9.320697\n",
       "7     9.237154\n",
       "8     9.250827\n",
       "9     9.324477\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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}