{
 "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.983\n",
      "Model:                            OLS   Adj. R-squared:                  0.982\n",
      "Method:                 Least Squares   F-statistic:                     884.9\n",
      "Date:                Fri, 11 Sep 2020   Prob (F-statistic):           1.14e-40\n",
      "Time:                        06:24:40   Log-Likelihood:              -0.080263\n",
      "No. Observations:                  50   AIC:                             8.161\n",
      "Df Residuals:                      46   BIC:                             15.81\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.0021      0.086     58.081      0.000       4.829       5.175\n",
      "x1             0.4995      0.013     37.607      0.000       0.473       0.526\n",
      "x2             0.4943      0.052      9.466      0.000       0.389       0.599\n",
      "x3            -0.0200      0.001    -17.150      0.000      -0.022      -0.018\n",
      "==============================================================================\n",
      "Omnibus:                        1.687   Durbin-Watson:                   1.755\n",
      "Prob(Omnibus):                  0.430   Jarque-Bera (JB):                1.341\n",
      "Skew:                          -0.206   Prob(JB):                        0.511\n",
      "Kurtosis:                       2.312   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.50212574  4.98049192  5.41996224  5.79360049  6.0841916   6.28707001\n",
      "  6.41088624  6.47618561  6.51203272  6.5512362   6.62495835  6.75759527\n",
      "  6.96276865  7.24108775  7.58004933  7.95609185  8.33846662  8.69428967\n",
      "  8.99394453  9.21594813  9.35048038  9.40099762  9.3836647   9.32469921\n",
      "  9.25606396  9.21021558  9.214772    9.28797491  9.43569268  9.65045622\n",
      "  9.91268612 10.19390951 10.46143804 10.68373852 10.8356137  10.90234214\n",
      " 10.88209687 10.78624558 10.63748292 10.46610185 10.30501659 10.1843549\n",
      " 10.12650817 10.1424526  10.22994581 10.37389529 10.54883802 10.72312402\n",
      " 10.8641163  10.94355306]\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.92976322 10.78534702 10.5296871  10.20695402  9.87529179  9.5925822\n",
      "  9.40227327  9.32274176  9.34279382  9.42440566]\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": "iVBORw0KGgoAAAANSUhEUgAAA6MAAAHWCAYAAACCFl9HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABrnklEQVR4nO3dd3RU1cKG8eckBAg1tFBC7yAgCIiCoIIaG4jYe+8Nr6Ji7w2s1957u4KooKKIKFYEAWkioFhCFyItQEjO90eAT5DOJGeSPL+1siBnZs688c4d8s7eZ+8gDEMkSZIkSSpICVEHkCRJkiQVP5ZRSZIkSVKBs4xKkiRJkgqcZVSSJEmSVOAso5IkSZKkAmcZlSRJkiQVuBJRB6hatWpYv379qGNIkiRJkvLBuHHjFoVhWG3T45GX0fr16zN27NioY0iSJEmS8kEQBL9t7rjTdCVJkiRJBc4yKkmSJEkqcJZRSZIkSVKBs4xKkiRJkgqcZVSSJEmSVOAiX013W5YuXcqCBQvIzs6OOoqKkKSkJFJTU6lQoULUUSRJkqRiKa7L6NKlS5k/fz5paWkkJycTBEHUkVQEhGFIVlYWGRkZABZSSZIkKQJxPU13wYIFpKWlUaZMGYuoYiYIAsqUKUNaWhoLFiyIOo4kSZJULMV1Gc3OziY5OTnqGCqikpOTnf4tSZIkRSSuyyjgiKjyja8tSZIkKTpxX0YlSZIkSUWPZVSSJEmSVOAso/ng9NNPJwgCgiDYsIXI/vvvz6OPPrpD1yiOGjWKIAhYtGhRPqaVJEmSpIJnGc0nBxxwAHPnzmX27Nl8/PHH9OzZk5tuuomuXbuyYsWKqONJkiRJUqSKRRkdMj6DLnePpME1w+hy90iGjM/I9+csVaoUNWrUIC0tjbZt2/Kf//yHUaNG8cMPP3DvvfcC8Morr9CxY0fKly9PamoqxxxzzIa9L2fPns3+++8PQLVq1QiCgNNPPx2Ajz76iK5du1KpUiUqV65Meno606ZNy/efSZIkSZJipciX0SHjM+g/eBIZmVmEQEZmFv0HTyqQQrqpVq1acfDBBzNo0CAA1qxZwy233MLEiRMZOnQoixYt4oQTTgCgTp06G+43ZcoU5s6dy0MPPQTAihUr6Nu3L2PGjGHUqFFUrFiRnj17smbNmgL/mSRJ+qcoPgCWJBVOJaIOkN8GDJ9OVnbORseysnMYMHw6vdulFXieli1bMmLECADOPPPMDccbNmzI448/TosWLfjzzz+pXbs2lStXBiA1NZWqVatuuO9RRx210Tmff/55KlSowJgxY9hnn30K4KeQJOnf1n8AvP7f3fUfAAOR/JsrSYpvRX5kdE5m1g4dz29hGG7Y3/KHH37giCOOoF69epQvX54OHToA8Pvvv2/1HLNmzeLEE0+kUaNGVKhQgerVq5Obm7vNx0mSlJ+29gGwJEmbKvIjo7VSksnYTPGslZIcQRqYOnUqDRs2ZMWKFaSnp3PAAQfw8ssvk5qayqJFi+jates2p9v27NmTtLQ0nnzySdLS0ihRogQtW7Z0mq4kKVLx9gGwlF+GjM9gwPDpzMnMolZKMv3Smzn6L+2EIj8y2i+9GclJiRsdS05KpF96swLPMnnyZD766COOPvpofvrpJxYtWsSdd95Jt27daN68OQsWLNjo/iVLlgQgJ+f/P2X+66+/mDZtGtdeey0HHHAALVq0YNmyZaxdu7ZAfxZJkja1pQ96o/oAWMoP8bQeiVTYFfky2rtdGnf1aU1aSjIBkJaSzF19Wuf7p1erV69m3rx5zJkzh4kTJ3L//fez33770b59e6688krq1q1LqVKleOSRR/jll18YNmwYN9xww0bnqFevHkEQMGzYMBYuXMjy5cupVKkSVatW5emnn2bmzJl8/vnnnH/++ZQoUeQHuSVJcS6ePgCW8ovT0aXYKRYNpne7tAKfOjFixAhq1qxJYmIiKSkptGrViptuuonzzjuPkiVLUrZsWV588UWuvfZaHn30Udq0acP999/PwQcfvOEcaWlp3HLLLVx33XWcffbZnHrqqbzwwgu8+eabXHrppbRq1YrGjRtz3333/WtRI0mS/qkgphWuP5/TF1WUOR1dip0gDMNIA3To0CEcO3bsZm+bNm0aLVq0KOBEKk58jUkqDjZd5RbyRiwLYqaQVNR0uXvkZtcjSUtJ5qtrukeQSIp/QRCMC8Oww6bHt2uabhAE3YIgeC8IgowgCMIgCE7f5PY+QRAMD4Jg4brb94tJakmStMucVijFjtPRpdjZ3mtGywGTgcuAzc1BKAt8DfwnRrkkSVKMOK1Qip2o1iORiqLtumY0DMMPgA8AgiB4YTO3v7zutqqxDCdJknZdvG1zJhV2UaxHIhVFRX41XUmSijunFUqS4lEkq+kGQXAucC5A3bp1o4ggSVKx4Sq3kqR4FEkZDcPwKeApyFtNN4oMkiQVJ04rlGJn1SoYOhRq1oQuXaJOIxVexWKfUUmSJGlX5ObCl1/CKy+HzHh9LIeueIuqjTOhZ3ko/4+vcuU2/v6fX2XLQoJXyUnrWUYlSZKkLZg+HV5+Gb58fgb7znmNfsGrNAlnkJtUkmB5ZXhqGaxYsf0nLFv2/8vp7rvDTTdBq1b59wNIcWy7ymgQBOWAxuu+TQDqBkHQFlgchuHvQRBUBuoCKevu0zgIgkxgXhiG82KaWJIkxaUVK+DTT+HDd1axeOwvJKZWIblOVarXSqR6dahRY+OvChUgCKJOLf3bggXwxhsw7Ln5NJ/4BifxKrfzPWEQkLvPvnDqVSQcdRRUqpT3gNxcWL4cli37/z+39LX+9qVL4aOPYNAgOOkkuOUWaNgw2h9cKmDbOzLaAfjsH9/fsu7rReB0oBfw/D9uf/of97t5lxJKkqS49eefedfOffn2PMp9Poz0tUMZwCeUI2+kKIcEFlGV+VRnPtX5nRp8v+7vi5Oqk12pOmFqdRLTalC+QVWOOq4E++5rSVXBy8qCd9+FQS8so+zH73BC+CoXMYJEcslu1RZOG0Bw/PEk1q797wcnJOR9ulKhwo496eLFcM898PDDee333HPh+uvzLkaVioEgDKNdP6hDhw7h2LFjN3vbtGnTaNGiRQEn2jXBNv71PO2003jhhRcKJoy2qTC+xiQpSrm5MG4cvP9eyIy3xtPk56EczlD25HsAVlWrQ1KfniR27Qx//w3z5hHOm8+aP+eTM2c+wfz5JC2eR4k1/973NJsSvMORfNjoErpduw8nnBhQunRB/4QqbjIy4Nbr15D55kf0yXqVI4L3KB2uYk1afUqedmLeqGXLlvkbYs4cuO02eOYZSEqCSy+Fq6/+/5FXqZALgmBcGIYd/nXcMhpb8+b9/6zkoUOHcs455zB37twNx5KTk6lYseKG77Ozs0lKSirQjPp/hfE1JkkFbcUKGDEChr+zkuXvfkqXzLwCmsYcwiAgq00nko/pSdDzcGjdetvDmmGYN1Vx/vyNvtZO/om1L7xM6ZVLGE9bni9/KVUuPJ5zLk2mVq2C+VlVfOTmwpNPhEy8/FluXXMtqSxkSamKLO7Vm0Z9z4W99y74IfpZs/KuIX3ttbxR1quugssuy7vOVCrEtlRGXc4rxmrUqLHhKyUlZaNjq1atIiUlhddff53u3buTnJzMk08+yQsvvEC5cuU2Os+oUaMIgoBFixZtOPb111+z7777UqZMGdLS0rjgggtYunRpQf54kqRiZOlS+M9pi7g65QkSex/OfS9W4aXMXpxZ+jWq9twbXniBYN48ykz4huC6a6FNm+375T0I8hZvadw4b1+MPn3gggso8ehDlF74J+GTT9GkwVoeXnYmF91Th5frXMulR/7B99/n/8+s4mHaNDipw3SaX9SdJ9acw681qnPG0TfR8ZIXObzZ8QxJrhfNXPFGjeCVV2DCBOjWDa67Lu/YI4/AmjUFn0fKZ5bRCPTv358LL7yQqVOn0rt37+16zKRJkzjooIPo1asXEydOZPDgwUyYMIEzzzwzf8NKkoqlzz6Daxr9j2tfas4jay+gR82plLzoXPj4Y5IyF1HqvbfhtNMgNTUmzzdkfAZd7h5Jg1s/Y5/FjRjx9ofw2WeUPrAb/cJ7uH9IA2bveQwXth7NW2+GrF0bk6dVMbNmDdx502rean0rL4xvwx5JP3B1+iUcf+qdfNaoI2sTS5CVncOA4dOjDdqmDbz3Hnz1FTRvDpdcAs2awUsvQU5OtNmkGCp0W7v07Zv3YVFBatsWHnwwdue75JJLOProo3foMQMGDOC4447jiiuu2HDs8ccfp127dixYsIDUGP0yIEkq3rKy4La+f9H6qYt5jDdY3qIjvPIxye3a5dtI0ZDxGfQfPIms7LxfsjMys+j/zmTo05reHw+G2bPJfuAxej7zDMdMfpvxx7flmkqXUOPyEzjjwmSqVMmXWCpivvkGHj9xNP1nn0sLfmJV7+M5IO0wFpb793WZczL/fU1zJDp3zvtk6OOP4dpr8z4AuvdeuP12OOIIV/pSoefIaAQ6dPjXdOltGjduHK+88grlypXb8NWlSxcAZs2aFeuIkqRiaMwYuLzJUC55qhXHJAxizU23U+7Hr2GPPfL1l94Bw6dvKKLrbTQ6Vb8+pR66l9IL/yT3iadoVC+HgUvO4vQb6/Bcjf7cedl8Vq/Ot3gq5JYtg6vPXcLkzufy0uxu1EvNgg8+oPQ7r1Oy9uYvRq6VklzAKbciCCA9HcaOhf/9D9auhSOPhL32yiuqUiFW6EZGYzlCGZWym1yEnpCQwKYLSWVnZ2/0fW5uLmeffTaXX375v86XlpYW+5CSpGJjzRoYcP3f1Bp4OU+Ez7OsYRtKDPowb2pQAdjSKNS/jpcpQ8J551Dh3LPh889Juv1hrhh5L5kPP8lNg//LsUNOZI/2jhTp/w0bGvLB6W9yw199qRYsYs0lV1Lmzps3LAjUL73ZRqPyAMlJifRLbxZR4q0IAjj6aOjdO2+67s03Q/fuedvCXHJJ1OmkneLIaByoVq0aK1eu3GgxogmbzEXeY489mDJlCo0bN/7XV3JyHH16J0kqVCZPhktbjuCUAa05jRdZdcW1lJ86psCKKGx5FGqLo1NBAPvtR8URg0mYOoWgeXPu/vNk/uzYm4FXzHWdF7FgAVzW61cSeh7Ko3+dQLnmdUj8YSwlHxqw0cq0vdulcVef1qSlJBMAaSnJ3NWnNb3bxfEH/SVKwJlnwvTpeSOkl14Kt96at0q1VMhYRuNAp06dKFu2LP3792fmzJkMGjSIxx57bKP7XH311YwZM4bzzz+f8ePHM3PmTIYOHcp5550XUWpJUmGWkwMP3rGCL3e/iCdmHUiltDIkfPM1pQfeAaVKFWiWfunNSE5K3OjYdo9ONW9OpcmjWXnbfRyc8DFn3r8btzR5hR8n+ot5cTNkfAad7xpJjUPHMjDtHu56fzd6lBzN2oEPUm7yt1v8gKV3uzS+uqY7v959GF9d0z2+i+g/JSfDW2/B6afnbQdz+eV5+9VIhYhlNA5UrlyZV199lU8++YTWrVvz1FNPcdttt210nzZt2vDFF18we/Zs9t13X3bffXf69+9P9erVI0otSSqsfvkFLmn3JYdfvzvn5j7OyvMup/yM8dCpUyR5dnl0KjGRMtf/h5JTJxI2b8Edv5/C7+2O4KGr57jqbjExZHwGV70+leyn4MMPz+HetdfwVb09+GzIZ5S44jJITNz2SQqjEiXg2WfzVvh86CE46yx80aswCTa9VrGgdejQIRw7duxmb5s2bRotWrQo4EQqTnyNSSpOwhCefXQVy/9zA5dm38eKavUp978XCPbtFnW02MnJYfld/yXp5mtZmVOKB+s/xHFDT6Hlbl5LWpTtccU3tH5qCo8tv4zlScnceNh5DG+6N2mVyvDVNd2jjpf/wjBvhd0bb8y7pvT116F06ahTSRsEQTAuDMN/reLqyKgkScXAypVwxb5j6XzJHvTNHsjKk8+j/C8/Fq0iCpCYSLnr+1Jq2kTWNm/FLbNPY3brnjx2XYbbMxYiG/advWYYXe4eyZDxGVu87/uD1nDmg6/w4vJz+LFaUw4/9wGGN+sMQRA/W7TktyCAG27IW8xoyBA4/PC8ZYSlOGcZlSSpiFu5Em7bZzh3jt6HOhWXkvvBR5R7+XEoVy7qaPmnSROqTfmcZbc9SI+EkZx4527c2fQFpv/ktaTxbv2+sxmZWYSs23d28KR/FdKcHBjY90+qHL0fF+c+xpNtjuaU025mYbnKG+4TV1u0FIRLLslbaXfUKDjgAFi8OOpE0lZZRiVJKsKysuC2Lh9x0/gjWFG3BeVnTiDhkPSoYxWMhATKX38ZJaf9yJrmbbjhlzP4ZbfDeeqmDNd5iWPb3HcWWLIEruv8Gac81J49SvzI17c/wYO9zmJt4v/vWhi3W7Tkt1NOgUGDYOJE6NYN5syJOpG0RZZRSZKKqKwsuL3zB9w84QhW1GtJlR9GQNWqUceKme2dyhk0aUzqlFH8fdvD7B+M4thbd+Oe3V5i3rwCDqztsq19Z3+cGPJU43u5Y8wBlKxRmVITv6fzdecVvi1a8tMRR8CHH8Jvv8E+++StWibFIRcwUrHma0xSUZWVBXd0HsYNE/qwrF4rqv7wCVSuvO0HFhLrp3L+cwQtOSlxmwUknDmL+YedQY2fR/Na8pnUe+8RuhyQP1M5h4zPYMDw6czJzKJWSjL90pttsxztzGOKmi53jyRjM4U0LSWZy6p2oOR5Z9A7dzCLuh9D1SHPQvnyEaQsJMaMgUMOyduu6eOPoVWrqBOpmHIBI0mSiomsLLiz89C8Ilq/dZErorB9Uzk3J2jciBpTP2PBuddzYtZzlDtwL56/dgax/mx+e6973NXHFEWb23e2dEIJ2o4uz+7ndOTw3HdZetN9VB3xpkV0W/bcE774Im+Bo27d4Lvvok4kbcQyKklSEZKVBXd1fp8bJvRhaf02RbKIwrancm5VYiKpT97Giv99QMOkP+lzVwfu6zI4pouP7kxZ3tmCXdRsuu9stcQK9HxjPvd+sB81y/wNn46kws3/yStY2rbddoMvv8x7H+jRA0aMiDqRtEGJbd9FkiQVBllZcE/nd7l+wjH83aAt1X74GFJSoo6VL2qlJG92Kue2Vk/ddBrsDe98SJsLLuHKb47i+Xr/odNnd9Ny96StPmZ7ps7uTFne2YJdFKf29m6XRu92aYz5cg2TDu7HWSseZkHTLqR+9hbUqhV1vMKnQQMYPRrS0+Gww+CNN+DII6NOJTkyKklSUbBqFdzbeQjXTjiGvxu0K9JFFDY/lXNbq6dubhrs5d8uZuz/XufPIy/mjCX3k7nH/rz7WMZWH7M9U2e3VIq3VpZ35jFFeWrvqwPmkN2te14RPakvqZM/s4juipo187Z82WMPOPpoeOGFqBNJllHtuosvvpj99ttvw/enn346hx9++C6d8+abb6aVF9lL0mZtuorsW99mcO/e73DthGPIbLgH1cYX7SIK/57KuT2rp25pGuw9n82m9uD/svixN2jLRPa+qB2PHjmCNWt2furszpTlnXlMUZzau3o1DDj8c3pctQd7JExg+TNvkPrKA5CUtO0Ha+sqV86bptujB5xxBjz4YNSJVMxZRvNJRkYG5557LrVr16ZkyZKkpaVxzjnn8Oeff250v20Vt4kTJ3LEEUdQo0YNSpcuTd26dTnqqKP47bff8vtH2GkPPfQQr7zyynbdd/bs2QRBwKYrKl955ZV8/vnn+RFPkgq1TUfC/ly0mnePGE7/CceypGEHUsd/DBUrRh2zQPRul8ZX13Tn17sP46truu/y1NnKFxxH0oTvyamSygVDDuLFxrfxx+9rd+hc/8y2o2V5Zx6zS9fOxqE//wh5rOkDXD6sB4mVK1Lyh+8od9ZxUccqWsqWhfffh6OOgssvh8ceizqRijGvGc0Hv/76K507d6ZBgwa8+OKLNGnShFmzZnHdddfRsWNHvvnmG+rXr7/N8yxcuJAePXqQnp7OsGHDqFKlCr/99hvDhg1j6dKlMc28Zs0aSpYsGZNzVYzBL0HlypWjXLlyMUgjSUXLP0fCwrUJdHl5Di8uOI8fK7dmj/HDoUKFiBPGr+25zjSpdXNq/vYdvx18Hud8eSN1njqIvr3+w6qma7f4mC1Zf93jjtjRx+zstbPxaPSHy1l85Flcvvot5nQ6klofv+DrOb+UKpV33WifPnDppdCwIRx8cNSpVAw5MpoPLrroIhISEhgxYgQ9evSgbt267L///owYMYKEhAQuuuii7TrPV199xZIlS3j++edp37499evXZ9999+Xee++ldevWW3zc+tHW22+/nerVq1OuXDnOOOMMsrL+/x+r/fbbjwsuuIArr7ySatWq0aVLFwCmTp3KYYcdRvny5UlNTeWEE05g3j92Bc/JyeHKK6+kUqVKVKpUib59+5KTk7PZ518vDEPuu+8+mjRpQqlSpahduzb9+/cHoEGDBgB07NiRIAg2TPfddJpubm4ut912G3Xq1KFUqVK0bt2ad999d8Pt60dYBw0axIEHHkiZMmVo2bIln3zyyXb9t5akwmL9iFe4NoF9Xs7gpQXnMTalNaeeeq2/uG/Ddk+DLVuWel+8zNybnmD/nFF8/M45NP542YbtX7Y1dbYg7czU3ngThvDitdOpcuieHL76beZfcQ+1vhnk6zm/lSgBr70GrVvDscfC5MlRJ1IxZBmNscWLF/PRRx9x0UUXUaZMmY1uK1OmDBdeeCEffvghS5Ys2ea5atSoQW5uLm+//TbhDm6A9vnnnzNx4kQ+/fRTBg0axMcff8zVV1+90X1eeeUVwjBk9OjRvPTSS8ydO5du3brRqlUrxowZw4gRI1i+fDm9evUiNzcXgPvuu4+nn36aJ598km+++YacnBxeffXVrWa59tprue222+jfvz9Tpkzhf//7H3Xq1AFgzJgxAHz00UfMnTuXwYMHb/YcDz30EAMGDOCee+5h0qRJHHnkkfTp04cJEyZsdL/rrruOSy+9lIkTJ9KxY0eOP/54li9fvkP/7SQpntVKSSbMCej28h+8uOB8vq/UhrNOv5aK1atEHS3u7dA02CCg5s3nsfbzrylROoEPxp/KUc9PIDWp3DanzhaknZnaG09WroSH9h3MkXd1JK3kQla99wnVB17lti0FpVy5vCm75crB4YfD/PlRJ1IxU/im6fbtC5sUkHzXtu12X+A9Y8YMwjCkRYsWm729ZcuWhGHIjBkz2HPPPbd6rr322otrr72W0047jYsuuoiOHTuy3377cdJJJ1GvXr2tPjYxMZHnn3+ecuXK0apVK+655x7OOuss7rrrLsqWLQvkjUred999Gx5z4403svvuu3PPPfdsOPbSSy9RuXJlxo4dy5577smDDz7IVVddxbHHHgvklcThw4dvMcfy5ct54IEHePDBBznzzDMBaNy4MXvvvTcA1apVA6BKlSrUqFFji+cZOHAgV155JSeeeCIAt956K1988QUDBw7c6PrUyy+/nJ49ewJw55138tJLLzFhwgT22Wefrf73kqTC4sqDmvHK0SN5fsGFjKm0O2ef3p+wbLlCNRIWpR2dBlu2W3vKZIzn126ncf+U6/l4YC/K7/kstMvHkDtoZ6YDx4NfZ6xl1D7X0XfBvWTU3pOaX75NQr06UccqfmrXziuk3brBEUfAZ59BcuGb5q3CyZHRfBJs4RO99SOcW7p9U3fccQfz5s3jqaeeonXr1jz77LO0bNmSTz/9dKuPa9OmzUbXXO69996sWbOGWbNmbTjWvn37jR4zbtw4vvjiiw3Xa5YrV27DCOasWbP4+++/mTt37oYiCZCQkECnTp22mGPq1KmsXr2aHj16bNfPuzlLly5lzpw5G6YSr7fPPvswderUjY61adNmw99rrVv+fcGCBTv93JIUb5a9mcnzv17E9HJNOef0/lRKrZxvI2GbrtpbFLYL2RlB5Uo0/HEIv1/+APut+oh6vdrw8umfsslVKtoBo95awO8t0zljwb38fuj5pM38wiIapfbt4ZVXYMyYvFV2182Ik/Jb4RsZjfMlqJs0aUIQBEyZMoXevXv/6/Zp06YRBAGNGjXa7nNWqVKFY445hmOOOYa77rqLdu3acdttt+1SwQM2jJCul5uby2GHHcbAgQP/dd/q1atvmKq7I3Z0evHWbK7Ab3os6R/Lvq+/bWdyS1I8+vTNRexzb09ySpWhxeThTMnHX97Xr9q7frGk9ftXAoVyFG6XJSRQ9/6+LO+zH+FhJ3DSiwfy+sir2HfUrdRuGJsFAIuDMISXLvqO7o8fTbVgEQvueZ66V50edSwBHHkk3HMPXHUVNGkCt90WdSIVA46MxljlypVJT0/nscceY+XKlRvdtnLlSh599FEOOeQQKleuvFPnL1myJI0aNdrmdZCTJk1ixYoVG77/9ttvNzx2S/bYYw+mTJlCvXr1aNy48UZf5cuXp2LFitSsWZNvv/12w2PCMNxw3efmtGzZklKlSm1xJHf9Cr6bLoL0TxUqVKBWrVp8+eWXGx3/8ssvadmy5RYfJ0lFybSJayh9Uh/SgjmU+vDdfB9FKor7V8ZCuX3akjZnLDP3P4eT/riHhU27MOLxGVHHKhSWLQ15ao8nOOHxrpQsm0T45dekWkTjy5VXwllnwe23w0svRZ1GxYBlNB888sgjrF27lgMOOICRI0fyxx9/MGrUKA488EDCMOSRRx7Z6P5Lly5lwoQJG33Nnj2boUOHcvLJJzN06FB+/vlnpk+fzsCBA/nggw848sgjt5ph7dq1nHnmmUyZMoVPPvmEa665hnPOOedfo6H/dNFFF/H3339z3HHH8d133/HLL78wYsQIzj33XJYtWwbAZZddxr333svbb7/N9OnT6du3L3Pnzt3iOcuXL89ll11G//79ef7555k1axZjxozh8ccfByA1NZXk5GSGDx/O/Pnz+fvvvzd7nn79+jFw4EBef/11fv75Z2688UZGjx7NFVdcsdX/DpJUFPy1KGRy1/PpkjOaZQ89T5n9t3x5RKwUtf0rY6psWZqOfJI5/x1EQ2ax14XtePmAF1mVFbvZQEXNzxOzGFH3DM6bcAEZzQ8g9bexJHeOowtvlScI8vYd3X9/OPtsGD066kQq4grfNN1CoFGjRowdO5Zbb72VU045hQULFlCtWjUOPfRQ3nzzTWrXrr3R/UePHk27dhu/IR911FHce++9lCtXjiuvvJI//viDEiVK0KBBAwYOHMhll1221Qz77rsvu+22G/vvvz8rV67ccL6tqVWrFl999RX9+/fn4IMPZtWqVdStW5eDDjqIUqVKAXDFFVcwb948zj77bABOOeUUTjrpJKZNm7bF8951111UqlSJ2267jT///JPq1atz6qmnAlCiRAkefvhhbr31Vm655Ra6du3KqFGj/nWOSy+9lGXLlnHVVVcxf/58mjVrxqBBg2jbtu1WfyZJKuyys+F/e9/P+cue548zbqTOJScUyPMWpf0r80uti/uw+uCO/LH/KZzy6el8WPMjGg5/nGadUqKOFjdyc+Gtq8ex2/1nckTuJH497WYaPHcDJDgeErdKloRBg2DvvfOm7n77LTRuHHUqFVFBLK/p2xkdOnQIx44du9nbpk2btsVVabVlp59+OosWLWLo0KFRR4l7vsYkxbvHDxvKeR/04veOR1P/2zcK7Jf4Ta8Zhbz9KwvTtiEFJieH6WfeQ6OXbiQjqM3Eq16j512di/3uJD9PWMnYw27iuDn3k1mqOrlPPUu1Uw+JOpa218yZ0KkTVK2aV0grVYo6kQqxIAjGhWHYYdPjfiwlSVKcevO6Hzn5gxOYU30P6o96oUBHkwr7/pUFKjGRZi9eS+b7X1KydAKH3dOV/7W+lcxFa6NOFom1a+H1sz+lRLvWnDhnILP2PYvKc6daRAubxo1hyBD49Vc46ihYsybqRCqCnKYrSVIcGv32fPa6sydrSlWg5ph3oUyZAs9QWPevjErVw/cid84Epu1/IcdOuInva39Cwmuv0L7P1vcGL0omj17CzCOv5IS/nmNu2cYsfukzmvbZL+pY2lldu8Izz8Bpp8GFF8LTT1Psh/wVU46MFkEvvPCCU3QlqRCbOXkVJY/vQ2qwkFLD3yOxroWwsEhIqcBu419h5k0v03zNRBodtTtv7/swCzOK9qjS6lUhbxz9NlW7teDwv15keu+rqbnwRypbRAu/U0+F666DZ5+FzWz/J+0Ky6gkSXEkc0nI5C7n0inna5b+9yXK7ds+6kjaCY1vPhl+GM+8tPYc/cVlLK/TnHeOfZ2Vy4ve3tfj3p/D6NQ+HD/oGFZVrsWKkWNo9s7dkOxiV0XGrbfCscfC1VfDO+9EnUZFSNyX0agXWFLR5WtLUrxZuxYGd7qH3ktf5tczb6P6RUdHHUm7oHzbRjT/YwR/PP0RlK/Akf87kVmVO/Dhfz5hbRG4nHTl8lzeOuApGvdqwT7LPuKnM+6h/vwxVNx/j6ijKdYSEuCFF2DPPeGkk2DcuKgTqYiI6zKalJREVpb7mSl/ZGVlkZSUFHUMSdrgpSPf4cwZ/Zm55wk0eOa6qOMoFoKAOmen02DJD/x84ytUSVjCIQ8cxJiUA/nigXEU1s9Fv315Bj9W686xn57H/Fp7sPaHH2n+3FVQwuVIiqzkZHj3XUhNhZ494c8/o06kIiCuy2hqaioZGRmsXLnSUSzFTBiGrFy5koyMDFJTU6OOI0kAvHPjeI4bejKzq3ei8ahnXSSkqElIoOktJ1Ez8ycmnvEgLbLG0+0/HRhZ/QQmDJoVdbrt9suEpQze8y7antqalmsmMP3Kp2n650jKtWsSdTQVhOrVYehQWL4cDj88709pF8T1PqMAS5cuZcGCBWRnZxdgKhV1SUlJpKamUqFChaijSBLfvjOXtD57Uqp0QOUZYyhRu0bUkZTPshf9zcRTB9Lyw/tJYg0jGp9P81duoEGn+PuQdPWqkNH3fkP2E8/Qbe6blGUlkxofSeOPHiG5Ua2o4ykKH30Ehx0Ghx6at/1LYmLUiRTntrTPaNyXUUmSirLZ07L4q81+NM+ZQu7nX1K+a9uoI6kALZ8xl2kn3Eq7cU+TRTKj97ySjq//h2oNy0cdjZnfLGRq/5dp9uUzNMuZxvKgHD+3P4HaN55Fas9OUcdT1B59FC6+GP7zH7jvvqjTKM5ZRiVJijMrV4SMrHMqhy55lbmPDibtwt5RR1JEFn71M7+feh3tf3mbBUEqP7Q5g+TDe9DynC5Uqxf7PWaHjM9gwPDpzMnMolZKMv3Sm9G7XRpZK3L5+tYRJDz3DF0WDaEk2UyvvDfZp51Ny5uPJaFCuZhnUSF2ySXwyCPw3HNwxhlRp1Ecs4xKkhRHwhBe2fNhThl7GT+ffCtNX74h6kiKA7Pf/I7Mvjex27xPSWIta0jix7J781eb7lQ4sgetz9qTcpVL7tJzDBmfQf/Bk8jKztlwrNYfq7lw0ni6TXmZurm/sSShMjP3PpUGt59F1f1a7eqPpaJq7Vo4+GD44gv47DPo0iXqRIpTllFJkuLI4Ms+p9fDPZjZ/HCaTxmct3WCtM7azOXMeOErMgd9SuWJI2my7AcSCFlBGSaldGVp++5UObYHrU5uS6ky23+9XhhC55s/Y+kfS6nwdxZNf53D8VM/4qDVn5JILj9WP4CEc86m5bW9SUgulY8/oYqMxYuhUydYuhTGjIF69aJOpDhkGZUkKU6MGfQH9Y9uz5qylan15xgSUlxMTVuXNWcJPz/9OSveH0mNqZ/SMGsqAEtIYWrqfqzs1J0ybRqzZv4S1s7/CxYvJliymBJLF1N6xV8kr1pMuTWLScn5i0osIYH///0vI6jJ63V78c6+nfnqxVOj+hFVmP30E+y1F9SvD19+CeWczq2NWUYlSYoDc39dxbxm3Wiy9idyvx1DhT2bRx1JhdDf0+cx46nPWPPRSGrPGEnd7F/+fZ+gIstKVGZ56SpkJVcmu1xlcipWZszSgPlJZVhatixzq1RifJv65CYmkpaSzFfXdI/gp1GRsH6F3SOOgLffdraHNrKlMurOxJIkFZA1q0N+2OtCDsv+nt8efId6FlHtpIrNatDhvhPgvhMAWDBmNpnT5lK+XmUqNKhCmVopVEwqQcXNPHb++Aye2OSa0eSkRPqlNyug9CqSDj4YBg7MW1335pvh1lujTqRCwDIqSVIBeSf9CY5b8DxT+tzAbpf1jjqOipDUPeuTumf97bpv73ZpAJtdTVfaJX37wqRJcNttsNtucNxxUSdSnHOariRJBeDD67+ixx37M7PBgbSc+b5T2CQVTatXQ48eMG4cjB4NHf41M1PF0Jam6W7Xv4RBEHQLguC9IAgygiAIgyA4fZPbgyAIbg6CYE4QBFlBEIwKgmC3GGWXJKlQm/TxHNrecTQLkuvRdMyrFlFJRVepUjB4MKSmQu/eMHdu1IkUx7b3X8NywGTgMiBrM7dfBVwBXAJ0BBYAnwRBUD4WISVJKqz+mruGNb2OpkKwjOQP36FE1ZSoI0lS/kpNhffeg8zMvEKatbn6IG1nGQ3D8IMwDK8Nw/BtIPeftwVBEAB9gbvDMBwUhuFk4DSgPHBijPNKklRo5OTAN50uo/3qb5hz+/NU2bdV1JEkqWDsvju8/HLe3qPnnJO3ya20iVjME2oA1AA+Xn8gDMMs4AugcwzOL0lSofRer2c4/I8nmHDw1TS59pio40hSwTryyLzFjF59Fe69N+o0ikOxKKM11v05f5Pj8/9x20aCIDg3CIKxQRCMXbhwYQwiSJIUXz6/9zsO/eAipqQdSNuhd0QdR5Kicd11cPzx0L9/3tRd6R9iuYLCpmPvwWaO5d0xDJ8Kw7BDGIYdqlWrFsMIkiRFb+ZX82l8zVH8VbIWjb57HRITo44kSdEIAnjuOWjfHk46KW/rF2mdWJTReev+3HQUNJV/j5ZKklSkLVucTeaBx1ApXEzCu+9QOq1K1JEkKVrJyTBkCJQvD716gTMjtU4syuiv5BXSA9cfCIKgNNAV+DoG55ckqVAIQxi95xV0yBrNr9c9Q42D2271/kPGZ9Dl7pE0uGYYXe4eyZDxGQUTVNpJvma109LS8grp3Llw9NGwZk3UiRQHtnef0XJBELQNgqDtusfUXfd93TAMQ+BB4JogCPoEQdAKeAFYDryWP7ElSYo/w457iUNn/Zex3S5nt9u3vqD8kPEZ9B88iYzMLEIgIzOL/oMn+cu94pavWe2yPffMm7L7xRdw8cWusKvtHhntAIxf95UM3LLu77euu/1e4H7gUWAsUBM4KAzDZTFNK0lSnBp952gO/N85TE7dn/Yjtr1q5IDh08nKztnoWFZ2DgOGT8+viNIu2dnXrKOp2siJJ+YtZvT00/DII1GnUcRKbM+dwjAcRd6CRFu6PQRuXvclSVKxMuXdmbS47kjmla5Pox/eJkja9j+vczI3vwn8lo5LUduZ1+z60dT1JXb9aCpA73ZpsQ+pwuH222HqVLj8cmjeHA48cNuPUZEUy9V0JUkqEnZkJGfO5MWUPuowEhKgzMhhJKdV3q7nqJWSvEPHpajtzGvWGQDarIQEePllaNkS+vSBb76JOpEiYhmVJOkfduS6uBVL1jBn76OonTObxc++Q7W9G2/38/RLb0Zy0sZbviQnJdIvvdmu/ghSvtiZ16wzALRF5cvDRx9BjRpw8MEwdmzUiRQBy6gkSf+wvSM5uTkh3+1+Lh2Wj2LaFc/S+PSuO/Q8vdulcVef1qSlJBMAaSnJ3NWntVMXFbd25jXrDABtVa1aMHIkVK4MBx0EEydGnUgFbLuuGZUkqbjY3pGcj7vfxcF/vMi36Tex18CTd+q5erdLs3yqUNnR12y/9GYbXTMKzgDQJurUySuk3brBAQfA55/nTd9VseDIqCRJ/7A9IzmfXfAWB39xHWOanEinD24CXDFU2hxnAGi7NGiQV0hLlIAePeDnn6NOpAIShBHv79OhQ4dwrHPEJUlxYtPVPyFvJGf9L9A/PPYtLS/aj5kpHWj2+wiSypfe5mMkSdth6lTYbz8oWTJvL9KGDaNOpBgJgmBcGIYdNj3uyKgkSf+wtZGcX0f+Sp2Le7EgKY06379DUvnSgCuGSlJMtGwJI0ZAVlbeCOnvv0edSPnMa0YlSUXakPEZDBg+nTmZWdRKSaZferNtjlZu7rq4Jb9msvaQw0kim1XDhlGxcbUNt7liqCTFSJs28PHH0L173tcXX+QtdKQiyZFRSVKRtSPbtGzNmhXZzGx/LPXX/MyfDw6izoHNN7rdFUMlKYbat4fhw2H+/LwR0vnzo06kfGIZlSQVWbGYPhvmhnzd/mI6LvmE7895ilaXdv/XfdwzVJJibK+94IMP8qbqHnAALFoUdSLlA8uoJKnIisX02VG97me/6U/xRedr6PzUGZu9jyuGSlI+6NoV3nsPZszI24c0MzPqRIoxrxmVJBVZtVKSydhM8dze6bPf9R/CvsP68W3to9nn8zu2el/3DJWkfNCjB7zzDhxxBKSnwyefQIUKUadSjDgyKkkqsnZl+uy0V8fR+u4TmVq2I7tPfImEEv6TKUmROOQQ+N//4Icf4NBDYfnyqBMpRvyXVZIK0JDxGXS5eyQNrhlGl7tH7vBCOtoxOzt9dtZ7U6hy6mH8lZhK6jfvkVzZhYgkKVJHHAGvvQbffAO9euVt/6JCz2m6klRA1q/sun5BnfUruwJO78xHOzp9dsabP1D5hIPIDkqy+p0PqdO6ej6mkyRtt2OOgdWr4dRT4cgj4d13oVSpqFNpFzgyKkkFJBYruyp//fT8N6Se0J2shLJkffQFjXu2iDqSJOmfTj4Znn46b+uXY46BNWuiTqRdYBmVpAISi5VdlX+mPDqK2mceyJLEquSMGk2DAxtHHUmStDlnnQWPPgrvvw/HHQdLl0adSDvJMipJBWRLK7hu78quyj+TBg6n4cWHMC+pLklff0G9fepGHUmStDUXXggPPZS39cvuu8Po0VEn0k6wjEpSAdmVlV2Vfybc8i5N+/VidqnmlP/hc9I61oo6kiRpe1x6KXzxBSQkwL77wjXX5F1TqkLDMipJBWRnV3ZV/vnhqjdodfNRTC/Tjqo/jqR6q2pRR5Ik7YguXWDCBDj7bLjnHujUCSZPjjqVtlMQhmGkATp06BCOHTs20gySpOJn3MXP0+7RsxhfrisNJg+lcr3yUUeSJO2K99/PK6WZmXDnnXD55XmjpopcEATjwjDssOlx/9eRJBU7Y894lPaPnsnYigfQ+OcPLaKSVBT07Jk3KnrooXDlldCjB/z2W9SptBWWUUlSsfL9cQPp8MLFfFWlJy1mvEfFmmWijiRJipVq1WDwYHjuORg7Ftq0gZdegohng2rzLKOSpOIhDPm+5610fKsfn1c/lnazBlG+WumoU0mSYi0I4Iwz4Mcf88roaafl7Um6aFHUybQJy6gkqegLQ8YecA0dh97Ep7VPo9Os1yhTMSnqVJKk/NSgAYwalbew0XvvQevW8OGHUafSP1hGJWknDRmfQZe7R9LgmmF0uXskQ8ZnRB1JmxGuzWFcl0vpMPJehjc4n64znqN02cRtP1CSVPglJsJVV8H330PVqnnXk15wAaxYEXUyYRmVpJ0yZHwG/QdPIiMzixDIyMyi/+BJFtI4s/DzqcxI7UL7bx7hg+b/ocf0xyhZ2n/6JKnY2X33vEJ65ZXw5JPQrh18913UqXZebi6MGQM33ghnnhl1mp1WIuoAklQYDRg+nazsnI2OZWXnMGD49M3uG5qTA9OmwZjPs1gwfDxBbg4la1WlTN2qlK9XmWo1EklNzVt3oVo1SHIG6S4Js9cy8eQBtHjrZhIox/vHv8qhL59AYokg6miSpKiULg0DBsBhh+VdR9qlC1xzDVx0EdSsGXW6bVu2DD75BIYOhWHDYMGCvK1runaF7OxC+cuDZVSSdsKczKytHl+8OO8D1ynD/yRr5NdUnv4NHdZ8zSn8QBJrN3pMLgGZpLCIqsymKmOpyrKSVckqW5XsilUJK1chqFaVsrs3pt1JLWnVOiCwU23RwpGT+LvPGbT9exwjKx9FnXcfpec+1aOOJUmKF/vtl7e40aWXwh135H21a5c3hfeww2DPPfOm98aDX3/NK59Dh+Zd/7pmDVSsCIccAocfDgcfDFWqRJ1ypwVhxMscd+jQIRw7dmykGSRpR3W5eyQZ64pnmAvZf5Un54+y7P7HAvZd/DONFnxDZ76mLn8AsCaxNIsb70nJfTtT6ZC9CMqWIVy4iKw/FrHy90WsmbOInIV/ESxaRIm/F1Fq2SLKZi2iZO7qjZ53LjX4pnR3Fu9xAFWO7UHn4+tS3Z4FQLgmm4kn3EXLwbeTSQrfnfoohz53TNz8PiFJikOTJ+cVvQ8+gK+/zpvKVLlyXsk79FBIT8+71rSgrF0L336bl+n992Hq1LzjzZrllc+ePaFz50I3ChoEwbgwDDv867hlVJJ23JDxGfR98DfaDp/N3gvHs1fuGDryPWXIK6h/V6zDqj06k3JIZ0rt3znvWpUd/YcjDGHlSoaPmsxbH/5ApSmz6DrrJ7rOH0PlNQsBmEFjxlc5gNVdelD3tP3pdGgVShfD3UoWfDye5cecQcOlE/mk6vE0fP9hGu1VLepYkqTCZMkS+PjjvGL64YewcGHeNjF77ZVXTA89NG8ENdbTkzIzYfjw/y/FixdDiRKw7755BfSww6BJk9g+ZwGzjEoqFoaMz2DA8OnMycyiVkoy/dKbbfYazl2xYAE8dtY4Dh56EXvxHWuCJKZVag5d9qLNyQcSdN4bateOyXOtXyjpn9enJpdI4L+tS9JszCRWffApaTNGUSZnObkETAjaMbPeASQe1INmZ+3Dbh3LFOkpveHqNfx47O20fO8uFlOFMWc+waFP9XY0VJK0a3JzYezYvHL4wQd5ix9B3rWlhxySV0wPOCBvyuym1q7NK7Z//ZVXLP/66/+/Nv3+r79gypS8Edn1q/0efjgcdNDmz11IWUYlFXmbLW5JidzVp3VMCunatfD8fYsJbryeM9c8wYqyqSTddw+lTz0WkpN3+fyb88/pwP+UlpLMV9d0z/smO5uVn3/PHy98SuKoEdTL+IYksllNSb4vuQ9zO/Skxrm96HRCQ0qWzJeYkZg/bCxZx59B/eWT+ajaKTQZ9iCNOlaOOpYkqSiaPx8++iivmA4fDn//nTd62bkzlCmzcbn8++8tnycxMW8acJUq//9nq1Z5BbRTp/i5VjXGLKOSirztKm476esvc/n4xBe46I+rqcxi/j7lEir/95Z8/9SywTXD2Ny7dAD8evdhm3/QihUsHDyaea9+SoVvPqLe0skATEloxcwWvSh/8hF0uqgDZcsXzi1OwqxVTDrqZnb7cADzgxqMPftJDn/icBIK548jSSpssrPhm2/yiumnn+Ydq1Jl44K56d/Xf1+hAsXxHyzLqKQib6eK2zYsWACPnf0D6e9fxN58y8JmXaj6xqMEbXffpazbKxYFe9XUX5h5/3skDnuXJvNGU4Ic5lKTifV6UuLIXrS7sgdV0uL/QtPcv5bwy0PvU/L+u6i74ieGVT+Tlh/eR4N2KVFHkyRJW7GlMlr8armkIqtWyuanym7p+Nbk5MCzA5fwbp2LueH9jrQq8wurnnyRatNGF1gRBeiX3ozkpI2n7CQnJdIvvdl2n6N0y4a0eqYvLeZ+BvMXMLX/y8xrvA/7/P4aBzx4OKVqV2V0tT58cvKLZExcFOsfYZfkzF/Ez1c/y9R6h7C2anUa33Yaa1eu4b0LP+KQOc9aRCVJKsQcGZVUZMTqmtFvv87loxNf4sLfrqIKf/H3SRdR+ZFbISUlH1Jv284syrQ9jwlXrWbGU5+x5MX3qPfje9RYm0EOCfxYrjOLW+xDqS4dqHtUR+p0rkOQUHCrIK3NmM+Me9+BQYNokvEZJcjhVxowofHRlD75aPa+tCMplYrwqkySJBUxTtOVVCzsymq6CxfCY+dO4IAhF9GFr1nUtDNV3niUoF3b/A0dYztVysOQ2YN/IOOxd6k85kMaLZ9ISbIBWBhUY3bVjqxs2YHy3TvS4JgOVGpRI6aZ18yew4x7B5P4zts0nfcFCYT8HDRlcrOjKXva0XS+sC3lK1hAJUkqjCyjkrQVw99exuyTr+Ps1Y+SVaYKiffdS/K5pxbKRQZicZ1p9rJVzBoyiYUffA9jx1L9j+9ptHoqieQCMLdEbTJqdGDN7h2pfFAHGhzdnlI1K+ct6rBmTd6f2dnkrs5mbVY2OVlryFmVvdFX7qo1LPh4PKXef5umC78GYGrCbkxreTQVzzqazufsRpmyFlBJkgo7y6gkbcGHr2dS5aR02odj+fuEC6j86G1QqVLUsXZafizkBLB07gpmvDWeJZ+MJWni99SeO5ZGOT/v9PnW+zGxLTPaHE3lc45i7zOaUzr+11KSJEk7YEtltEQUYSQpXnzwWibVTj6Itkwg69XBVD7xiKgj7bJaKcmbHRndmYWc/qlCzbK0v2wfuGyfDcfmTPubWW+NY9nn42DFSsKkJCiRRFgiCZKSCJNKQskkgqS870lKgpIlCUomEZRMokLreux5fEPaFKH9TyVJ0vaxjEoqtj54dQmppxzE7kxkzWuDKHd8z6gjxUS/9GabvWZ0R1bg3V61WlSk1k3dgV3bx1WSJBU/llFJxdIHryymxikH0iqYzJrXB1P2uMOjjhQz6xcp2tmFnCRJkgqCZVRSsTPs5cXUOvUAdgumkP3mO5Q95tCoI8Vc73ZpcV0+d2XVY0mSVDRYRiUVK8Ne+ou00w6gRTCN7LffpWyfg6OOVOxsuvVMRmYW/QdPArCQSpJUjBS+PQskaScNfWERdU7rTotgGmstopEZMHz6RtezAmRl5zBg+PSIEkmSpCg4MiopbsVyKuf7zy2k3lk9aJYwg7Vvv0fZIw/Kl+fRts3ZzEq/WzsuSZKKJsuopLgUy6mc7z+7gPpn96Bpwkxy3nmfsr0OyJfn0fbJr61nJElS4RKzabpBEJQPguDBIAh+C4IgKwiCr4Mg6Bir80sqXmI1lfO9p+fT8Oz9aZIwi5whQynzjyIay+fR9uuX3ozkpMSNjuXX1jOSJCl+xXJk9BmgDXAa8CdwMjAiCIKWYRhmxPB5JBUDsZjK+d5T82h8XncaJswm971hlDls/3x5Hu0Yt56RJEkQozIaBEEycBRwVBiGo9YdvjkIgp7ABcD1sXgeScXHrk7lfPeJuTS9oDv1E34n9/0PKHPofvnyPNo58b71jCRJyn+xmqZbAkgEVm1yPAvYJ0bPIakY2ZWpnMOemUuzC/anXsIfhMM+3GIR3dXnkSRJ0s6LychoGIbLgiD4Brg+CILJwDzgBGBvYOam9w+C4FzgXIC6devGIoKkImZnp3L+NGEVaeceRt2EP+GDDymT3jVfnkcFz1WPJUkqWoIwDGNzoiBoBDwHdANygB+An4E9wjBsuaXHdejQIRw7dmxMMkgq3rKy4N3aF3H84sf46/n3qHJ6z6gjKUY2XfUY8kaw7+rT2kIqSVKcC4JgXBiGHTY9HrPVdMMwnBWG4b5AOaBOGIZ7AknAr7F6Dknampd6vc3xix/j16OusIgWMa56LElS0ROzMrpeGIYrwjCcGwRBJSAdeDfWzyFJmxr60CyOH3EWv9fqRIPX74o6jmLMVY8lSSp6Yra1SxAE6eSV25+AxsAAYDrwfKyeQ5I2Z+aU1dT6z3EkJCZQ6/M3ICkp6kiKMVc9liSp6InlyGhF4BHyyuhLwJfAQWEYZsfwOSRpI6tXw5j9r2KP3HGsfuJ5SjSuH3Uk5QNXPZYkqeiJ2choGIZvAW/F6nyStD1e7vMOZy98mF969aXh2b2jjlMsFcQqt656LElS0ROz1XR3lqvpStpZw5/4lU4XtGNp9abU/f1LKFky6kjFjqvcSpKkbcn31XQlqSDN/nkNVS8+jsREqPH5mxbRiLjKrSRJ2lmWUUmFTnY2fLPvNbTP+Z4V/32Oks0aRB2p2HKVW0mStLMso5IKnVePe48T5j3AjIMvpsYFfaKOU6xtaTVbV7mVJEnbYhmVVKiMfOE3er1zOr9X3YMmQwZGHafYc5VbSZK0s2K2mq4k5bc/f82m/DnHk5SQQ/XP34JSpaKOVOy5yq0kSdpZllFJhcLatTC667WcsPZbMh58i7SWjaKOpHV6t0uzfEqSpB1mGZVUKLx+8jBOyRjI9B4X0OyyY6KOU2QVxJ6hkiRJYBmVVAiMfu0PDn3zVH6r3JZmQ++POk6RtemeoRmZWfQfPAnAQipJkmLOBYwkxbV5f2RT8vQTKJ2whmoj34LSpaOOVGS5Z6gkSSpIllFJcSsMYVS3G+mU/RWL736KMrs3iTpSkeaeoZIkqSBZRiXFrWE3fsexs+9hapezqdPvhKjjFHnuGSpJkgqSZVRSXFo0N5v6d53LopK1aD70vqjjFAvuGSpJkgqSCxhJikuf9nyA43J+5Pf73yEhpULUcYoF9wyVJEkFKQjDMNIAHTp0CMeOHRtpBknx5ZvXfmX3k3ZjdtN0Wk5/J+o4kiRJ2gVBEIwLw7DDpsedpisprqzKCll7zgXkBCVoMPS/UceRJElSPrGMSoorQ09+g64rh/PnBXeQ3KR21HEkSZKUTyyjkuLGz98upuvgvsyssictHr4w6jiSJEnKR5ZRSXEhDOHnI6+mCn9R6a2nIDFx2w+SJElSoeVqupIKxJDxGVtdpfXD60Zz+LxnmHhQP3bvvnuESSVJklQQLKOS8t2Q8Rn0HzyJrOwcADIys+g/eBKQt53IoozVNL73XOaUqk/rt2+KMqokSZIKiNN0JeW7AcOnbyii62Vl5zBg+HQAvjj8Hprm/ET2Q4+TUL5sFBElSZJUwCyjkvLdnMysLR7/9sXpHDbhDia2PJ565x1cwMkkSZIUFcuopHxXKyV5s8drlClDcMH5rEooQ9NhDxZsKEmSJEXKMiop3/VLb0Zy0sar4yYnJXLqyGl0yhrFn5fcS3L96hGlkyRJUhRcwEhSvlu/au4/V9M9s0Z1et9+Ej9V24fd7j8r4oSSJEkqaJZRSQWid7u0DaU0DOGTmqdQnmXkvv0UJDhJQ5IkqbjxN0BJBe6Tqz7hoPmvMOmwa6jWrUXUcSRJkhQBy6ikArXojywa338Bv5duSru3ro06jiRJkiJiGZVUoL497DYa5s4i59EnSChTOuo4kiRJiohlVFKB+f65SaRPGsC41qfT4Mz9o44jSZKkCFlGJRWI1Vm5JF50HksTUmgxbGDUcSRJkhQxy6ikAvHpyc+zx6pvyLj8PsrUqRJ1HEmSJEXMMiop382dlknHwf2ZWrkLbQacEnUcSZIkxQHLqKR8N7HPLVRhEeWf/y8EQdRxJEmSFAcso5Ly1cTXp3LAT//l+7bnUqdXu6jjSJIkKU5YRiXlm9yckFXnX8byoDy7vXN71HEkSZIURyyjkvLNZ33fpdPSEfx80q2Uq1816jiSJEmKI5ZRSfli6fwsmjx2ObOSd6PDsxdEHUeSJElxxjIqKV98ffR91M2dTfbAh0koWSLqOJIkSYozllFJMffL53/Q7cs7+b7+0TS/sHvUcSRJkhSHLKOSYu7PE/oREFLvrYFRR5EkSVKcsoxKiqlv7/2CbnPfZPyBV5PasV7UcSRJkhSnLKOSYmbNyrWk3HAJGSXq0uGtq6KOI0mSpDhmGZUUM1+c+jTN1/zIvCvvo2RKmajjSJIkKY5ZRiXFxIJpf7HH4OuZWHl/2t95VNRxJEmSFOcso5JiYvJRN1Ih/JsKzz8EQRB1HEmSJMU5y6ikXTb19YnsO+0Jvm17AQ16tY46jiRJkgqBmJTRIAgSgyC4LQiCX4MgWLXuz9uDIHCne6mIC3NDVp1/GZlBJdq8c0vUcSRJklRIxKosXg1cBJwGTALaAC8Cq4HbYvQckuLQl5f9j65LP+erU56gS/3KUceRJElSIRGrMtoZeD8Mw/fXfT87CIL3gE4xOr+kODFkfAYDhk9nTmYW9RITeeWxK5me3Ja9nz076miSJEkqRGJ1zeiXwP5BEDQHCIKgJdAd+CBG55cUB4aMz6D/4ElkZGYRAoe8OJTauX/wfd9bSUhKjDqeJEmSCpFYldF7gJeBqUEQZANTgBfDMHxsc3cOguDcIAjGBkEwduHChTGKICm/DRg+nazsHABq/J7JZRnP8FaFI3i8QtmIk0mSJKmwiVUZPQ44FTgR2GPd3y8MguCszd05DMOnwjDsEIZhh2rVqsUogqT8Nicza8Pf+7//KjkkMqD3cRsdlyRJkrZHrK4ZHQAMDMPwjXXfTwqCoB7QH3g2Rs8hKWK1UpLJyMyi47ezOGL5h9xWvy8La1YgLSU56miSJEkqZGI1MloGyNnkWE4Mzy8pDvRLb0b5MOS2Lx9lVkIDXj5iP5KTEumX3izqaJIkSSpkYjUy+j5wTRAEv5J3vWg74D/ASzE6v6Q40LtdGqv7PkfznJ85ufN9VKtRnn7pzejdLi3qaJIkSSpkYlVGLyFvP9HHgFRgLvA0cGuMzi8pDiyYMIdDv7iX76odzstf/ocgiDqRJEmSCquYlNEwDJcBfdd9SSqiZh51Fe1ZQ+prD1pEJUmStEu8plPSdpn82Bd0/uVVvup8FQ0OaBR1HEmSJBVyllFJ25Szei2lrriYPxPrsuc7/aOOI0mSpCLAMippm7497XGarJrE7EsfoFxqmajjSJIkqQiwjEraqiU/zafVWzfwfaUD6TLwyKjjSJIkqYiwjEraqp+O7E9yuJIKL/yXIMFViyRJkhQbllFJW/TzS9+y90/P88Uel9OsV7Oo40iSJKkIsYxK2qxwbQ7hhRcxN6EW7d+5Puo4kiRJKmIso5I2a8y5z9BsxQ9MO+s+KtUtH3UcSZIkFTGWUUn/smz2XzR54VrGld+X/R4/Luo4kiRJKoIso5L+ZXLv66kQ/k3SE4+QkOiiRZIkSYo9y6ikjfw6aBydJj7JyN0uoc2JraKOI0mSpCLKMippgzAnl6yzLmZhkEq7ITdHHUeSJElFmGVU0gbjL3+Jln9/y4QT7qVa44pRx5EkSVIRZhmVBEDW3EzqPnoV48t0psfzJ0cdR5IkSUWcZVQSAD8eeROVcv8i58FHKFHStwZJkiTlL3/jlETGhz/S4btHGNH4fDqc0y7qOJIkSSoGSkQdQFLEwpDMUy6mFJVoNfi2qNNIkiSpmHBkVCrmJl37Orv9NZoxve8irXXlqONIkiSpmLCMSsXYyt8XUWtAXyaW6kiP186KOo4kSZKKEcuoVIz9lH4p5XMyWfXIs5RK9u1AkiRJBcffPqViavq977LHT6/zUYfr6XR266jjSJIkqZixjErF0Oq5i6l87flMSdqd/T7qH3UcSZIkFUOupisVEUPGZzBg+HTmZGZRKyWZfunN6N0ubbP3nZp+Oa1yFvHzgx+wW5WkAk4qSZIkWUalImHI+Az6D55EVnYOABmZWfQfPAngX4V01sPDaDfpJd5rcz29LnZPUUmSJEXDabpSETBg+PQNRXS9rOwcBgyfvtGx7EV/U+7K8/ipxG50HX59QUaUJEmSNuLIqFQEzMnM2q7jkw++gjbZc/lpwDs0r1Fqh6b2SpIkSbHkyKhUBNRKSd7m8dlPfUy7cc8ytEU/9r2y44apvRmZWYT8/9TeIeMzCii1JEmSijPLqFQE9EtvRnJS4kbHkpMS6ZfeDIC1S5ZR6pJz+DmxOZ2H3wxs/9ReSZIkKT84TVcqAtZPrd3SlNvJh15FmzV/MOXWr2hapzSw/VN7JUmSpPxgGZWKiN7t0jZ7vecfL46k7bdP8G6j/9Dr+r03HK+VkkzGZornlqb8SpIkSbHkNF2pCMtdupzE885mVkJj9hx+G0Hw/7dta2qvJEmSlJ8so1IRNqnXtdRYPZufr36Omo3KbHRb73Zp3NWnNWkpyQRAWkoyd/Vp7Wq6kiRJKhBO05WKqDlvjmb3z//Lu3UvodcdXTd7ny1N7ZUkSZLymyOjUhEUrlhJ7hln8mvQgD2G37XR9FxJkiQpHlhGpSJo8pE3UDtrJpP7Pkud5mWjjiNJkiT9i2VUKmIWvPsNu33yAO/VOp/DBu4fdRxJkiRpsyyjUhESZq1i9Uln8mdQh1Yf3EuC/w+XJElSnPJXVakImXLMzdRZ8RPjznuahruXjzqOJEmStEWWUamImPvch7Qcdi/vp55Fr0cOijqOJEmStFWWUakIWDFhBmXPOYEpiW1o9elDJCZGnUiSJEnaOsuoVMiFS5exeN/erMktQebzQ2jQytVzJUmSFP8so1JhlpvLjL1PpebS6Xx6zpt0PaV+1IkkSZKk7WIZlQqxGaffQdOpQ3it3UCOfbJH1HEkSZKk7WYZlQqpOU++R5OXb+T9lFM46ovLCIKoE0mSJEnbr0TUAST925DxGQwYPp05mVnUSkmmX3ozerdL23D7inE/UeHCk5mQ2J7WXz9J2XI2UUmSJBUullEpzgwZn0H/wZPIys4BICMzi/6DJwHQu10aYebf/L3fEZTILc2yV9+hbYvkKONKkiRJO8VpulKcGTB8+oYiul5Wdg4Dhk+H3Fxm7n0y1Zb/wsgL3qbriXUiSilJkiTtGsuoFGfmZGZt8fjMk2+myU9DeaX9gxz3aLcCTiZJkiTFTkzKaBAEs4MgCDfzNSwW55eKk1opm592e9S0H2n8+m28U/lMjvv8QhcskiRJUqEWq5HRjkDNf3ztAYTAWzE6v1Rs9EtvRnJS4kbHdpufwS3v38HYxE60++pRypS1iUqSJKlwi8kCRmEYLvzn90EQnAUsBf4Xi/NL8WBbK9zGyvpzrn+upiXX8uSb97A0LE/WK4Oo37x0zJ9TkiRJKmgxX003CIIAOAt4JQzDlbE+vxSFba1wG2u926XlnTcnh1ktDiM16w8GXTyKE46P/XNJkiRJUciPBYwOBBoAz+TDuaVIbHWF23w08/jraDRjOC/t+SjHP9w5X59LkiRJKkj5UUbPAb4Pw3DClu4QBMG5QRCMDYJg7MKFC7d0NylubG2F2/wy96G3aPz2Pfyvyvmc+Nk5LlgkSZKkIiWmZTQIglTgCODprd0vDMOnwjDsEIZhh2rVqsUygpQvtrTC7ZaO76plQz+n4uVn8G2JLnT8+iHKlMmXp5EkSZIiE+uR0dOB1cAbMT6vFKnNrXCbnJRIv/RmMX+uv54dQsle6fwW1iP7tbep37RkzJ9DkiRJilrMyui6hYvOBt4Iw3BZrM4rxYPe7dK4q09r0lKSCYC0lGTu6tM65osXzb/rOVLOPoofg7YsHDyarsfUiOn5JUmSpHgRy9V09wOaACfH8JxS3Niwwm1+CEMy+g4g7eGr+SzpICqOGES3buXy57kkSZKkOBCzMhqG4WeAS6xIOyo3l9+Pv4q6/7uPd8scT4vvXqRpK6fmSpIkqWiL+T6jknZAdja/pZ9Dvc9e5LVKF7HvxIdJq5Mfi1xLkiRJ8cUyKkUlK4vf9zqWej8O5enaN3PUhBupXMXJBZIkSSoeLKNSBMIlmfy5R09qz/6KR1o+yhljLqRs2ahTSZIkSQXHMqpCb8j4DAYMn86czCxqpSTTL71Z/i00FAO5GXOZ1+5gqi+cxn+7vMEFI4+lpJeISpIkqZixjKpQGzI+g/6DJ5GVnQNARmYW/QdPAojLQpo9bSZL9jyICssX8PQRw7hk8IEkeImoJEmSiiF/DVahNmD49A1FdL2s7BwGDJ8eUaItW/XtBJa324eE5Ut589yRXPiORVSSJEnFl78Kq1Cbk5m1Q8ejsmzY56zdZ1+Wr07ikxtGc9aTexK4VpEkSZKKMcuoCrVaKck7dLzA5eSw+I7HKdkznT9zajHh0a854dYWUaeSJEmSImcZVaHWL70ZyUmJGx1LTkqkX3qziBL9v9Wffc3cOh2pfP2FfJPQhXlvf0nPC+tEHUuSJEmKC5ZRFWq926VxV5/WpKUkEwBpKcnc1ad1pIsXhXPn8Uf30yjVvQs5cxfwQKc3aDBzBPsdVSWyTJIkSVK8cTVdFXq926XFx8q52dksuPERyg24ieo5q3iqSn+avnAtlx9eLupkkiRJUtyxjEoxsHLoSJaefgk1/prKJ4kHk9H/Ic64pSlJSVEnkyRJkuKT03SlXRD+/ge/dTqWMj17sPKvLB4+4F12z/iA0++0iEqSJElbYxmVdsbq1cy5+E5WNWhO6pj3eSLtVv76fAqXftKL1Oru2SJJkiRti9N0pR209PVhrDq/L7WWzmRoyT6svO0+zr2yPgl+tCNJkiRtN8uotD3+/JP5jw9mzStvUuf3r8mgOW/3/pgTnz+QlJSow0mSJEmFj2VUxdKQ8RkMGD6dOZlZ1EpJpl96s3+tyBv+9jtz/juInDffpu6fX1Md+JHWDG36AN3euJAL25WMJrwkSZJUBFhGVewMGZ9B/8GTyMrOASAjM4v+gycB0KtiNn888DYMept6c78jDRhPWz6ufzulTjqa/c5rxgV1IgwvSZIkFRGWURU7A4ZP31BEAer8NY/078dR/4GrSVjxI/WAcUF7Pm16F2VPOYr9z21Cu9Sde67tGYGVJEmSiiPLqIqX3FxK/DKLHr/NpcXvf3DA7G9ou2oKAGPowIu73UvFM49ivzMb0j5l155qayOwFlJJkiQVd5ZRFVk5c+Yzf8QkFn8+idyJkyj36yRqZk7l89yVAOQS8F2wJ9ek3sjw3TtSqkt5vr1p35g9/6YjsABZ2TkMGD7dMipJkqRizzKqQm/t4qXMHTmNRaMmsXb8ZMrMmkTNvyZRee1CagG1gPmk8nOp1kyofS5/1mrCB4mV+bl2NXJqryEokUtyUiI39Woc01xzMrN26LgkSZJUnFhGVSiEuSGLpi5g3mfTWPrdNJg2jbJ/TKPGkmnUWJtBHaAOsJyyzEjajTHVe5HVuDWlOrSm2v6taLpPKl0r/v/56m64ljM3367lrJWSTMZmimetlOSYPo8kSZJUGFlGtU0FuQjP2uyQGZ/MZtHoaayeMI0SM6dRad406q6YRrVwCdXW3W8Z5ZhdugXTavZgfMMWJLVuQZX9WtOoR33apSRs83l6t0vL96my/dKbbXTNKEByUiL90pvl6/NKkiRJhYFlVFtVEIvwzJu1gikPf0o4dBgtfv2AFuGfG25bmJBKRoUWTGpxHLnNWlCmfQtS921B7U5ptE4KYvL8+WX9fx9X05UkSZL+LQjDMNIAHTp0CMeOHRtpBm1Zl7tHbnaqaVpKMl9d032nzpmdDePfnsX854ZR5bthtF82ilKsYXlQjp/qHkRu9wOp1K01tXq0oGydyrv4E0iSJEmKUhAE48Iw7LDpcUdGtVWxWoTnj1lr+PHR0eS+P4zmvwxjz9yfAfg9uRkT97mYKqccSsPTutKhVMldzixJkiQp/llGtVW7sgjPpM8WMe2ed6n89TA6LfuEw1jOakoyo/b+TEy/iIYXH0bdto2omx/BJUmSJMU1y6i2amcW4flt2kq+PeFBDpl4N61ZxoJStfml04lUOvkw6pzeg1blym72cQW5UJIkSZKkaFlGtVU7sgjP4oU5fHzqK+zz0XUcRwZTmx5BnaduJLVbO1KDrS82VBALJUmSJEmKH5ZRbdO2tkFZtQreu3QEzZ/rx/E5E/ilakcWPP4aLY/utt3PMWD49I1GXwGysnMYMHy6ZVSSJEkqgra9IaO0Bbm5MPTuyXxT6VCOffpAqpdcwu93v0bD+d+SugNFFGK3UJIkSZKkwsGRUe2U0W/NZeGFN3LEX8+xMrE8M84dQJOHLobSpXfqfLuyUJIkSZKkwseRUe2Qyd8u5+Umt9DuuCYc/teLzDj4UsrOnUWTJ6/c6SIKeQslJSclbnRsWwslSZIkSSq8HBnVdvnztxw+OfF5Dv76Bk5hHj+3OZp6r91F890ax+T8O7JQkiRJkqTCzzKqbfr6fxmUPb4nZ+SOZ3atvVn67CCaHtw55s+zrYWSJEmSJBUdllFt1ccP/0Szy9KpkrCE+Q+9Qf1LjoVtbNMiSZIkSdtiGdUWDbrqO/YdcBgJSYnkDB9F9f33iDqSJEmSpCLCBYz0L2EIL5/4IQcP6E52mYok//A1FS2ikiRJkmLIMqqNrF0Lz3d/meNf78WiKs2o9vPXJLdqFHUsSZIkSUWM03S1QVYWvNlxIGdO6ces+t1pOOEdgooVdupcQ8ZnuDKuJEmSpC2yjAqAzMW5fNTmKk7PuI8Z7Y6hyTcvQ6lSO3WuIeMz6D94ElnZOQBkZGbRf/AkAAupJEmSJMBpugIyZmfzRaPTOT7jPmakX0ST71/f6SIKeXuFri+i62Vl5zBg+PRdjSpJkiSpiLCMFnPTf1jBzy2OoFfmy8w6/TaafPhfSEzcpXPOyczaoeOSJEmSih/LaDE2bvgilnfqTrdVw/ntuqdo9Pz1MdlDtFZK8g4dlyRJklT8WEaLqc9e+I1yh+xDq5yJLHh8EPVuPydm5+6X3ozkpI1HV5OTEumX3ixmzyFJkiSpcHMBo2LovTsn0/66dCokrmDF2x9Ts3e3mJ5//SJFrqYrSZIkaUsso8XMxw9Ooet1XckpmUzC56OpvFfrfHme3u3SLJ+SJEmStsgyWozM+GEZ9a84itykUpT/8WtKNasfdSRJkiRJxVTMrhkNgqBmEAQvBkGwMAiCVUEQTA2CYN9YnV+7ZuWKkF+6n0Wj3BmsffkNi6gkSZKkSMWkjAZBkAJ8BQTAYUAL4BJgQSzOr10ThjB4/4dJ//t/zDrzTqoft1/UkSRJkiQVc7GapnsVMDcMw1P/cezXGJ1bu+i9a77muO+v5KdmvWj+zFVRx5EkSZKkmE3T7Q18FwTBm0EQLAiCYEIQBBcHQQw2rdQumfjJAtrfeywLk+vS5KsXY7KPqCRJkiTtqliV0YbAhcAvQDrwEHA3cNHm7hwEwblBEIwNgmDswoULYxRBm1qyKIdlvU6kavAXycMGkVglJepIkiRJkgTErowmAD+EYdg/DMPxYRg+DzzMFspoGIZPhWHYIQzDDtWqVYtRBP1Tbi4M73wz+6z6lIz+j1Jp/7ZRR5IkSZKkDWJVRucCUzc5Ng2oG6Pzawe9feYHHD/jdqbsdSaN7jgz6jiSJEmStJFYldGvgGabHGsK/Baj82sHfPP6bA548WRmp7Sl5aePRB1HkiRJkv4lVmX0AWCvIAiuC4KgcRAExwCXAo/G6PzaTnN+WUXyqUdTIiGXap+/TVAmOepIkiRJkvQvMdnaJQzD74Mg6A3cCdwA/L7uz8dicX7925DxGQwYPp05mVnUSkmmX3ozDmuVxpgufem9dhy//3cIdds0ijqmJEmSJG1WrPYZJQzDYcCwWJ1PWzZkfAb9B08iKzsHgIzMLPoPnsTsS4bRd96TTOl5NbtdfETEKSVJkiRpy2I1TVcFaMDw6RuK6HqpX2dy7ld9+bnWfuw2+PaIkkmSJEnS9rGMFkJzMrM2+r7UvJCnR95AZkJF6n39OpSI2YC3JEmSJOULW0shVCslmYx1hTR3TcBdrz9FQ37lvOMe5tl6Nbb62M1da9q7XVpBxJYkSZKkDRwZLYT6pTcjOSkRgBPe+I4+a97j7rYX0bPf1q8TXX+taUZmFiH/f63pkPEZBZBakiRJkv6fZbQQ6t0ujbv6tGbP7zO5Y+4dfFjtIFo922+bI5ybu9Y0KzuHAcOn52dcSZIkSfoXp+kWUvumpdL087tZWKo2B/30FomVK27zMZtea7qt45IkSZKUXxwZLaRGHvUoLXMmk3Xng9tVRCHvWtMdOS5JkiRJ+cUyWghNGTGXA7+8kSl1D6Hh5du/n+g/rzVdLzkpkX7pzWIdUZIkSZK2ymm6hUwYwp8nXkVjVlN78MMQBNv92PXXlLqariRJkqSoWUYLmY+v/4L0ha8w4bDraNu+8Q4/vne7NMunJEmSpMg5TbcQWbp4LXXuvZi5JevS5o1ro44jSZIkSTvNMlqIfNrnUVquncSK2x8koVyZqONIkiRJ0k6zjBYSP42aR/fPb2RK7XQaX9k76jiSJEmStEsso4VAGMJvx19FMlnUfHvHFi2SJEmSpHhkGS0EPr3lS9Lnv8zkg/tRuVPTqONIkiRJ0i6zjMa5ZUvWUuuOi5iXVIfd33TRIkmSJElFg2U0zn16zOO0XPsjS295gMQKZaOOI0mSJEkxYRmNYz+Pns/+n17P5FoH0fSaPlHHkSRJkqSYsYzGqTCE2cddvW7Rov+6aJEkSZKkIsUyGqdG3fEVB819kUkHXkGVvV20SJIkSVLRYhmNQyv+XkvqrRcxN6kOu//v+qjjSJIkSVLMWUbj0KfHPMFu2RPJvOF+SlR00SJJkiRJRU+JqAMIhozPYMDw6czJzKLh3zDok+uZXOMAWl1/VNTRJEmSJClfODIasSHjM+g/eBIZmVnkhnDaK69QhpX8fPsNLlokSZIkqciyjEZswPDpZGXnALDbt/M5dfkbPFb3VO5buDbiZJIkSZKUf5ymG7E5mVl5f1kdcseXD/FHQhpP9j6EVeuPS5IkSVIR5MhoxGqlJANw1Hvf0y73R27b+3yykktvOC5JkiRJRZFlNGL90ptRcXnI1b88zufJe/Nxl7YkJyXSL71Z1NEkSZIkKd84TTdivdulseybZ0llIWf3uJm0SmXol96M3u3Soo4mSZIkSfnGMhqxBTOXcugPDzGh1iG89/6FUceRJEmSpALhNN2IjTvtYaqwmEoP3xp1FEmSJEkqMJbRCM37KZO9vr6PH2r3ot5RHaKOI0mSJEkFxjIaoQmnPUAlMqn66C1RR5EkSZKkAmUZjcicyYvpPOYBxtY/irq92kYdR5IkSZIKlGU0IpNOv49yLKfG4zdHHUWSJEmSCpxlNAIZExbSZdxDjGt4LLUPbhV1HEmSJEkqcJbRCEw5fQDJZFHrqZujjiJJkiRJkbCMFrA/xs5nn4mPMLbpiaT1aB51HEmSJEmKhGW0gE0/425KsoZ6z9wYdRRJkiRJioxltAD99s0c9pn8ON+3OJUaXZtEHUeSJEmSImMZLUAzzrqLRHJo+NwNUUeRJEmSpEhZRgvI7C9+p+u0p/i+1ZlU36tB1HEkSZIkKVKW0QLyyzl3EhDS+Pnroo4iSZIkSZGzjBaAWZ/OpuvPz/L97ueQ2qFu1HEkSZIkKXKW0QLwx7m3kUMizV68NuookiRJkhQXLKP5bOZHM9nnlxcZ2/58qu6eFnUcSZIkSYoLltF8Nuf8W1lDSVq8eE3UUSRJkiQpblhG89HP7/1El99eZVyni6iyW42o40iSJElS3IhJGQ2C4OYgCMJNvubF4tyF2fyLbyWLZFq9dFXUUSRJkiQprsRyZHQ6UPMfX61jeO5CZ/qgyXT54w3Gd7mESk2rRR1HkiRJkuJKiRiea20YhsV+NHS9vy69hZqUo81LV0YdRZIkSZLiTixHRhsGQZARBMGvQRC8EQRBwxieu1CZ9sZEOs95mwn79qViwypRx5EkSZKkuBOrMvodcDpwCHAOUAP4OgiCzTaxIAjODYJgbBAEYxcuXBijCPHj77438XdQkd1fvDzqKJIkSZIUl2JSRsMw/DAMw7fCMPwxDMMRwOHrzn3aFu7/VBiGHcIw7FCtWtG6nnL625PYa/67TNjvcirWqxR1HEmSJEmKS/mytUsYhsuBKUCT/Dh/PFvY716WU5bdn7kk6iiSJEmSFLfypYwGQVAaaA7MzY/zx6s/v5zNXrNfZ1z780hpWDnqOJIkSZIUt2K1z+jAIAj2DYKgQRAEnYC3gbLAi7E4f2Ex68L7yCWBpo97ragkSZIkbU2sRkZrA6+Tt9foYGA1sFcYhr/F6Pxxb9G0hXSc9CzfNTmZmh1rRx1HkiRJkuJaTPYZDcPw+FicpzCbfN5/6cYqat3fL+ookiRJkhT38uWa0eJm+dxl7P7lI4yp2ZtGh7eIOo4kSZIkxT3LaAyMu+BpKoVLKHfb1VFHkSRJkqRCwTK6i7JXrKHp0PsZX3E/Wp3VKeo4kiRJklQoWEZ30ZjLXqVmTgbZV1wTdRRJkiRJKjRisoBRcZW7NpeaL9/DtNJt6XjdQQAMGZ/BgOHTmZOZRa2UZPqlN6N3u7SIk0qSJElSfHFkdBeMu/FdGq6ZzqKzriFICBgyPoP+gyeRkZlFCGRkZtF/8CSGjM+IOqokSZIkxRXL6M4KQ8r+925+K9GQvQYcBcCA4dPJys7Z6G5Z2TkMGD49ioSSJEmSFLcsoztp0iOf03L5GH7p04+k5LzZznMyszZ73y0dlyRJkqTiyjK6k7Jvu5v5QXX2fOz0DcdqpSRv9r5bOi5JkiRJxZVldCfMens8eywczqTufSlbpfSG4/3Sm5GclLjRfZOTEumX3qygI0qSJElSXHM13Z2wqN89VKM8ezx1/kbH16+a62q6kiRJkrR1ltEdNGf0LDrM/h+ftb+SAxqm/Ov23u3SLJ+SJEmStA1O091Bv148kLWUoPnjfaOOIkmSJEmFlmV0ByyZNo/2Pz7P141Po3bHmlHHkSRJkqRCyzK6A6ae9xAlWUOtB/pFHUWSJEmSCjXL6HZaOfdvWn/5GF/VPJpmhzeJOo4kSZIkFWqW0e004YInqRAupdytV0cdRZIkSZIKPcvodli7fBWNhz7AmIoH0Pas9lHHkSRJkqRCzzK6HX7o+xKpOfPI/s81BEHUaSRJkiSp8LOMbkO4NocaLw9gUukO7H1d96jjSJIkSVKRYBndhgk3DabumpksPPMaEhIdFpUkSZKkWLCMbsWQH/4kGHgHPweNuaF6RYaMz4g6kiRJkiQVCZbRLRgyPoO3rn2Htmsm8mjTU8hYuZr+gydZSCVJkiQpBiyjWzBg+HS+bdSQ85rew7AD2wKQlZ3DgOHTow0mSZIkSUVAiagDxKs5mVmE5WH4kbv967gkSZIkadc4MroFtVKSd+i4JEmSJGn7WUa3oF96M5KTEjc6lpyUSL/0ZhElkiRJkqSiw2m6W9C7XRqQd+3onMwsaqUk0y+92YbjkiRJkqSdZxndit7t0iyfkiRJkpQPnKYrSZIkSSpwllFJkiRJUoGzjEqSJEmSCpxlVJIkSZJU4CyjkiRJkqQCZxmVJEmSJBU4y6gkSZIkqcBZRiVJkiRJBc4yKkmSJEkqcJZRSZIkSVKBs4xKkiRJkgqcZVSSJEmSVOAso5IkSZKkAmcZlSRJkiQVOMuoJEmSJKnABWEYRhsgCBYCv0UaYtuqAouiDqFiz9eh4oGvQ8ULX4uKB74OFQ8Kw+uwXhiG1TY9GHkZLQyCIBgbhmGHqHOoePN1qHjg61Dxwtei4oGvQ8WDwvw6dJquJEmSJKnAWUYlSZIkSQXOMrp9noo6gISvQ8UHX4eKF74WFQ98HSoeFNrXodeMSpIkSZIKnCOjkiRJkqQCZxmVJEmSJBU4y+hWBEFwYRAEvwZBsCoIgnFBEHSNOpOKjyAIbg6CINzka17UuVT0BUHQLQiC94IgyFj3ujt9k9uDda/POUEQZAVBMCoIgt0iiqsiajtehy9s5j3y24jiqogKgqB/EATfB0GwNAiChUEQvB8EQatN7uN7ovLVdr4OC+V7omV0C4IgOA54CLgTaAd8DXwYBEHdSIOpuJkO1PzHV+to46iYKAdMBi4DsjZz+1XAFcAlQEdgAfBJEATlCyyhioNtvQ4BRrDxe+ShBRNNxch+wGNAZ6A7sBYYEQRB5X/cx/dE5bf92PbrEArhe6ILGG1BEATfAT+GYXjOP47NAN4Ow7B/dMlUXARBcDNwdBiGrbZ1Xym/BEGwHLg4DMMX1n0fAHOAR8IwvGPdsWTyfvm6MgzDJ6PKqqJr09fhumMvAFXDMDw8qlwqfoIgKAf8DfQOw/B93xMVhU1fh+uOvUAhfE90ZHQzgiAoCbQHPt7kpo/J+0RCKigN101R+zUIgjeCIGgYdSAVew2AGvzj/TEMwyzgC3x/VMHbJwiCBUEQ/BwEwdNBEKRGHUhFXnnyfn9esu573xMVhU1fh+sVuvdEy+jmVQUSgfmbHJ9P3huOVBC+A04HDgHOIe+193UQBFWiDKVib/17oO+PitpHwKlAD/KmSO4JjAyCoFSkqVTUPQRMAL5Z973viYrCpq9DKKTviSWiDhDnNp3DHGzmmJQvwjD88J/fr7sI/RfgNOD+SEJJ/8/3R0UqDMM3/vHtpCAIxgG/AYcBg6NJpaIsCIL7gX2AfcIwzNnkZt8TVSC29DosrO+Jjoxu3iIgh39/opXKvz/5kgpEGIbLgSlAk6izqFhbv6Kz74+KK2EYzgH+xPdI5YMgCB4ATgC6h2H4yz9u8j1RBWYrr8N/KSzviZbRzQjDcA0wDjhwk5sOJG9VXanABUFQGmgOzI06i4q1X8n75WvD++O612ZXfH9UhIIgqAqk4XukYiwIgoeAE8krAD9tcrPviSoQ23gdbu7+heI90Wm6W3Y/8HIQBGOAr4DzgVrAE5GmUrERBMFA4H3gd/I+Yb0BKAu8GGUuFX3rVulrvO7bBKBuEARtgcVhGP4eBMGDwHVBEPwE/AxcDywHXosgroqorb0O133dDAwi7xet+sBd5K1g+k4BR1URFgTBo8ApQG9gSRAE60dAl4dhuDwMw9D3ROW3bb0O171f3kwhfE90a5etCILgQvL2jqpJ3l5nl4dh+EW0qVRcBEHwBtCNvAW1FgLfAjeEYTg10mAq8oIg2A/4bDM3vRiG4enrtjK4CTgPqETeYlsXhWE4ucBCqsjb2usQuAAYQt4+4Cnk/fL1GXnvkX8USEAVC0EQbOkX5VvCMLx53X18T1S+2tbrcN12QkMohO+JllFJkiRJUoHzmlFJkiRJUoGzjEqSJEmSCpxlVJIkSZJU4CyjkiRJkqQCZxmVJEmSJBU4y6gkSZIkqcBZRiVJkiRJBc4yKkmSJEkqcJZRSZIkSVKB+z9IyFdCnvH4MgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x1, y, 'o', label=\"Data\")\n",
    "ax.plot(x1, y_true, 'b-', label=\"True\")\n",
    "ax.plot(np.hstack((x1, x1n)), np.hstack((ypred, ynewpred)), 'r', label=\"OLS prediction\")\n",
    "ax.legend(loc=\"best\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting with Formulas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using formulas can make both estimation and prediction a lot easier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.formula.api import ols\n",
    "\n",
    "data = {\"x1\" : x1, \"y\" : y}\n",
    "\n",
    "res = ols(\"y ~ x1 + np.sin(x1) + I((x1-5)**2)\", data=data).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the `I` to indicate use of the Identity transform. Ie., we do not want any expansion magic from using `**2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept           5.002138\n",
       "x1                  0.499515\n",
       "np.sin(x1)          0.494250\n",
       "I((x1 - 5) ** 2)   -0.020000\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.929763\n",
       "1    10.785347\n",
       "2    10.529687\n",
       "3    10.206954\n",
       "4     9.875292\n",
       "5     9.592582\n",
       "6     9.402273\n",
       "7     9.322742\n",
       "8     9.342794\n",
       "9     9.424406\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.6rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
