{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color='blue'>Data Science Academy - Python Fundamentos - Capítulo 8</font>\n",
    "\n",
    "## Download: http://github.com/dsacademybr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statsmodels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Regression Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Para visualização de gráficos\n",
    "from pylab import *\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
    "np.random.seed(9876789)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Criando dados artificiais\n",
    "nsample = 100\n",
    "x = np.linspace(0, 10, 100)\n",
    "X = np.column_stack((x, x**2))\n",
    "beta = np.array([1, 0.1, 10])\n",
    "e = np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = sm.add_constant(X)\n",
    "y = np.dot(X, beta) + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       1.000\n",
      "Model:                            OLS   Adj. R-squared:                  1.000\n",
      "Method:                 Least Squares   F-statistic:                 4.020e+06\n",
      "Date:                Thu, 02 Aug 2018   Prob (F-statistic):          2.83e-239\n",
      "Time:                        12:16:44   Log-Likelihood:                -146.51\n",
      "No. Observations:                 100   AIC:                             299.0\n",
      "Df Residuals:                      97   BIC:                             306.8\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.3423      0.313      4.292      0.000       0.722       1.963\n",
      "x1            -0.0402      0.145     -0.278      0.781      -0.327       0.247\n",
      "x2            10.0103      0.014    715.745      0.000       9.982      10.038\n",
      "==============================================================================\n",
      "Omnibus:                        2.042   Durbin-Watson:                   2.274\n",
      "Prob(Omnibus):                  0.360   Jarque-Bera (JB):                1.875\n",
      "Skew:                           0.234   Prob(JB):                        0.392\n",
      "Kurtosis:                       2.519   Cond. No.                         144.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "model = sm.OLS(y, X)\n",
    "results = model.fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [ 1.34233516 -0.04024948 10.01025357]\n",
      "R2:  0.9999879365025871\n"
     ]
    }
   ],
   "source": [
    "print('Parameters: ', results.params)\n",
    "print('R2: ', results.rsquared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.933\n",
      "Model:                            OLS   Adj. R-squared:                  0.928\n",
      "Method:                 Least Squares   F-statistic:                     211.8\n",
      "Date:                Thu, 02 Aug 2018   Prob (F-statistic):           6.30e-27\n",
      "Time:                        12:16:44   Log-Likelihood:                -34.438\n",
      "No. Observations:                  50   AIC:                             76.88\n",
      "Df Residuals:                      46   BIC:                             84.52\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             0.4687      0.026     17.751      0.000       0.416       0.522\n",
      "x2             0.4836      0.104      4.659      0.000       0.275       0.693\n",
      "x3            -0.0174      0.002     -7.507      0.000      -0.022      -0.013\n",
      "const          5.2058      0.171     30.405      0.000       4.861       5.550\n",
      "==============================================================================\n",
      "Omnibus:                        0.655   Durbin-Watson:                   2.896\n",
      "Prob(Omnibus):                  0.721   Jarque-Bera (JB):                0.360\n",
      "Skew:                           0.207   Prob(JB):                        0.835\n",
      "Kurtosis:                       3.026   Cond. No.                         221.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "nsample = 50\n",
    "sig = 0.5\n",
    "x = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x, np.sin(x), (x-5)**2, np.ones(nsample)))\n",
    "beta = [0.5, 0.5, -0.02, 5.]\n",
    "\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + sig * np.random.normal(size=nsample)\n",
    "\n",
    "res = sm.OLS(y, X).fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [ 0.46872448  0.48360119 -0.01740479  5.20584496]\n",
      "Standard errors:  [0.02640602 0.10380518 0.00231847 0.17121765]\n",
      "Predicted values:  [ 4.77072516  5.22213464  5.63620761  5.98658823  6.25643234  6.44117491\n",
      "  6.54928009  6.60085051  6.62432454  6.6518039   6.71377946  6.83412169\n",
      "  7.02615877  7.29048685  7.61487206  7.97626054  8.34456611  8.68761335\n",
      "  8.97642389  9.18997755  9.31866582  9.36587056  9.34740836  9.28893189\n",
      "  9.22171529  9.17751587  9.1833565   9.25708583  9.40444579  9.61812821\n",
      "  9.87897556 10.15912843 10.42660281 10.65054491 10.8063004  10.87946503\n",
      " 10.86825119 10.78378163 10.64826203 10.49133265 10.34519853 10.23933827\n",
      " 10.19566084 10.22490593 10.32487947 10.48081414 10.66779556 10.85485568\n",
      " 11.01006072 11.10575781]\n"
     ]
    }
   ],
   "source": [
    "print('Parameters: ', res.params)\n",
    "print('Standard errors: ', res.bse)\n",
    "print('Predicted values: ', res.predict())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c1d6ba470>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAFpCAYAAABNgFv/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Wd0VFUXgOH3ZlLpVYTQUSlSBZSISAQpShWwADawfCiKiCKgKCoqIFhAioAUESygEJpKNYAQepGOVCWA1AAhCSlzvh+bUAOk3ClJ9rNWVjKTmXvPpMy+p+1tGWNQSimllGf5eLoBSimllNKArJRSSnkFDchKKaWUF9CArJRSSnkBDchKKaWUF9CArJRSSnkBDchKKaWUF9CArJRSSnkBDchKKaWUF9CArJRSSnkBX3eerFChQqZ06dLuPKVSSinlMevWrTtujCmcmse6NSCXLl2atWvXuvOUSimllMdYlnUgtY/VIWullFLKC2hAVkoppbyABmSllFLKC7h1DjklCQkJHDx4kLi4OE83xaUCAwMpXrw4fn5+nm6KUkopL+TxgHzw4EFy585N6dKlsSzL081xCWMMJ06c4ODBg5QpU8bTzVFKKeWFPD5kHRcXR8GCBbNsMAawLIuCBQtm+VEApZRS6efxgAxk6WCcLDu8RqWUUunnFQHZm7z//vsMGTLkut8PCwtj27ZtbmyRUkqp7CDTBeSwDZHUHbiYMr3nUnfgYsI2RLr3/BqQlVJKuUCmCshhGyLpM30zkVGxGCAyKpY+0zdnOCh//PHHlC9fngcffJCdO3cCMHbsWGrXrk21atVo27YtMTExrFixglmzZtGzZ0+qV6/Onj17UnycUkoplVaZKiAPnreT2ISkK+6LTUhi8Lyd6T7munXr+PHHH9mwYQPTp09nzZo1ALRp04Y1a9awadMmKlasyLhx47j33ntp2bIlgwcPZuPGjZQrVy7FxymllFJp5fFtT2lxKCo2TfenxrJly3jkkUfIkSMHAC1btgRgy5Yt9O3bl6ioKKKjo2nSpEmKz0/t45RSSnm506chb16PnT5T9ZCL5QtK0/2pldIK6GeffZbhw4ezefNm+vXrd90tS6l9nFJKKS/Xty9Eundd0uUyVUDu2aQ8QX6OK+4L8nPQs0n5dB/z/vvvZ8aMGcTGxnL27Flmz54NwNmzZylatCgJCQlMmTLl4uNz587N2bNnL96+3uOUUkp5ubVroW1bWLVKbvfuDR7MppipAnLrGsEMaFOF4HxBWEBwviAGtKlC6xrB6T7mXXfdxeOPP0716tVp27Yt9erVA6B///7cc889NGrUiAoVKlx8/BNPPMHgwYOpUaMGe/bsue7jlFJKeSFjIDwcGjeG2rVh0SLYu1e+FxwMt9zisaZZxhi3naxWrVrm6nrI27dvp2LFim5rgydlp9eqlFJeqXlzmDsXihSBHj2gSxfIk8dlp7Msa50xplZqHnvTRV2WZY0HmgNHjTGVL9w3GGgBxAN7gE7GmKj0N1kppZRykfPnISBAvm7aFB5+GDp1gqCMrT+yW2qGrCcCTa+6bwFQ2RhTFdgF9LG5XUoppVTGrV8PVavC5Mly+5VX4OWXvS4YQyoCsjFmKXDyqvvmG2MSL9xcCRR3QduUUkqp9HE6YcgQqFMHzp2D4t4fpuzYh9wZ+MmG4yillFIZd+gQPPMMLFwIbdrAmDFQsKCnW3VTGVplbVnWO0AicN39PpZlvWhZ1lrLstYeO3YsI6dTSimlbm7VKlixQgLxzz9nimAMGeghW5b1DLLYq6G5wVJtY8wYYAzIKuv0nk8ppewUtiGSwfN2cigqlmL5gujZpHzqt1AaA1u2wMyZsH8/PPUU1K8vmZ4WLZIAUKCAfC5Y8NKCIuU6cXGwfDk0bAiPPAJ79sCtt3q6VWmSroBsWVZToBdQ3xiTqaspnDhxgoYNGwJw5MgRHA4HhQsXBmD16tX4+/t7snlKKRdILlSTnBs/uVANcOOgfPYs9OsHYWGwb5/cV6SIBGOAHTsk0cTl/P1hyhRo187ul6GSHT4sQXjjRvm9FC2a6YIxpG7b0w9AKFDIsqyDQD9kVXUAsOBC2smVxpguLmynyxQsWJCNGzcCUgs5V65cvPnmm1c8xhiDMQYfn0yVR0UpdR03KlRzRUA+dw7mzYOYGHjySciRA375BSpXlqxOLVrIm3+yKlVgwwY4eRJOnJDPW7bIwiKA+fMlO9Szz0KxYq5/odnBunXQqhWcOgXff3/l7yOTuWlANsa0T+HuLF/SaPfu3bRu3Zr77ruPVatWERYWRrVq1YiKku3WP/74IwsXLuSbb77hv//+46WXXuKff/7Bx8eHYcOGUSf5H1Ap5XVSVahm2jR46SUJrNWrS0B2OGQo1Pc6b505cshjr2fxYhg0CN57D5o1g+efh4ceuv7x1I1NnSoXN4UKyXD1jX72mYBX/RV07y4jDnaqXh2+/DJ9z922bRsTJkzg66+/JjEx8bqP69atG2+99RZ16tRh//79NG/enC1btqSzxUopVyuWL4jIFIJysXxB0qt95RX44QdJrThtGlxIqQtkLHgOHAjPPQfjxsHEiTBrFoSGwh9/pP+Y2dlff0GNGjB9ukwdZHJeFZC9Tbly5ahdu/ZNH7dw4UJ27rxUk/nUqVPExsYS5IUbz5VSUqjm8jlkuKxQza5dMGMG9O8vw9J2915vv10Cc//+MGeOLBADiI+XHsndd9t7vqzm3DmZJ65cGT78EBISssyiOa8KyOntybpKzpw5L37t4+PD5YvJLy+zaIzRBWBKZSLJ88TJq6xvCzJ8ErCP2jWaAsGyctrVPS4/P1mIlGz0aOjWTfbPDhqUJXp8tjtwQOaLjx6F3btliiCLBGPIZNWePMnHx4f8+fPz999/43Q6mTFjxsXvPfjgg4wYMeLi7Y12j7srpWzXukYwy3s3YF/TnCyY+Cq1+7565cppd+vUSXrk338Pd9wBQ4fCDabKsp0//5TRg337YPx4CcZZjAbkNBg0aBBNmzalYcOGFL8sDduIESNYvnw5VatWpVKlSowdO9aDrVRKpUp8PLz+OjzwgPRWly2DMmU8155cuWDAANi8WVZld+8uveXszhgYOVJ+T3nzStKPpleXV8gatPyiG2Wn16qUVzNG9gVPny4LuAYOhMumqDzOGNnrXKwY3HMPnDkjc6WZJOOUrYyRof3ERCkQkS+fp1uUJraWX1RKqSzHsuRNvmFDqfzjbZLbl6xXL1loNnKk5GbODv79V4JxyZIyjB8YCFk8F0TWfnVKKXW5mBiZiwTZV+yNwTglXbpIb7ltW3jsMVnUlJUtWQI1a14ass+RI8sHY9CArJTKLqKjJRlH48Zw5IinW5M21arJ3OnHH0v+7EqVJGd2VmOMbLdp2FCG50eNsu3QYRsiqTtwMWV6z6XuwMWEbYi07dh20YCslMr6zpyRhUBLl8I332TKPMf4+cHbb0tqzurVoWxZT7fIXjExUqTj9dclJemqVVChgi2HTs5dHhkVi+FS7nJvC8oakJVSWVtUFDRpIm/wP/4IHTp4ukUZU6mS1PktU0Z6lB06wIQJ4HR6umUZ43RK5q3+/SVfeJ48th36RrnLvYkGZKVU1jZunBQgmDYNHn3U062x15kzcPAgdO4MdevCmjWeblHaxMTABx9IFa1cuaT9ffvaPl+cqtzlXkADMnDw4EFatWrF7bffTrly5XjttdeIj48nPDyc5s2bX/P4OXPmUKNGDapVq0alSpUYPXq0B1qtlEqVHj3kjb51a0+3xH5580J4uOTF3r9fEmd06iT5uL3d4sVSHev99+HXX+U+F2XdKpYv5TTG17vfU7J9QDbG0KZNG1q3bs3ff//Nrl27iI6O5p133knx8QkJCbz44ovMnj2bTZs2sWHDBkJDQ93baKXUjcXEQPv2sHOnbCGqVs3TLXIdHx9ZjbxzJ7z1lgQ6Pz9Pt+r6oqKkylXDhtL28HB4/HGXnrJnk/IE+TmuuO9i7nIvkjkDckSEZLSJiMjwoRYvXkxgYCCdOnUCwOFw8MUXXzB+/HhiYmKuefzZs2dJTEyk4IUN+gEBAZQv712/VKWytaQk6NgRfvoJduzwdGvcJ08eyYG9cyfkzi2JRFq1utT79BYvvSQ9+l69ZM64fn2Xn7J1jWAGtKlCcL4gLCA4XxAD2lS5sva1F/C+xCAp9TYfe0z2C8bEyDzJX3/JAgAfH6haFV57TWpiHj8u2XcuFx5+w9Nt3bqVmjVrXnFfnjx5KFmyJLt3777m8QUKFKBly5aUKlWKhg0b0rx5c9q3b49PNtgjp5TXM0ZW6YaFSS7oVq083SL3CwyUzwcPygVJs2bw8MMSAOvVkxEDd1uzRi4YypeXrVtvvin7jN2odY1grwvAV8t8UeT06UurCZ1OuZ0BxhisFP5Ar3c/wDfffMOiRYu4++67GTJkCJ07d85QG5RSNvnyS/jqKwnK3bp5ujWeVaaM5MUeMgRWrJCeaMWKcOiQe84fGyurv2vXlrntTz+V+8uWdXswzjSMMW77qFmzprnatm3brrnvhlasMCYoyBiHQz6vWJG2519lwYIFpl69elfcd/r0aVOgQAEzd+5c06xZsxs+/9ixYyZXrlypOleaX6tSKvUSEoy57z5j2rY1JinJ063xLufOGTNxojHt2xvjdMp9331nzJIll27b6cMPjcmf3xgwplIlY4YPN+b0afvPkwkAa00qY2Tm6yGHhEiGmv795XNISIYO17BhQ2JiYpg0aRIASUlJvPHGGzz77LPkSKG8V3R0NOGXDYNv3LiRUqVKZagNSikb+PrC/Pnw3XfZIs1imuTIIQu/vv9ehqydTnjvvUu95s8/hy1bpFebVsbIdOFvv10avTx3Dh58EP74Q47btaut+4qzKq32BPz777+8/PLL7NixA6fTycMPP8yQIUOIiIjgoYceuriAC+CHH35gwIAB7Nmzh6CgIHLmzMnQoUOpVevmxTy84bUqleXs3g3vvANjxsg2IJU6MTGyN3v06EsLZN97T/YFR0VJVrA77oDbb5fPgYFSaSl3btnXPWgQ7NkjH8lThwsXyuppYzwzV+2FtNpTGpUoUYLZs2dfc39oaCixKVwx1qtXzx3NUkrdzPHj8NBDcOoUHDumATktknvNzzwji782bpQsYCALwn74QQLz5X76SRbZxsXJ48uVk1HKcuXkuckrpjNhMA7bEMngeTs5FBVLsXxB9GxS3u2LwDQgK6Uyp9hYaNlSyvQtXgy33ebpFmVeFSpcmTe6cmVJLnLiBOzaJR/nz8sCLZDdLrt2eaatLpCc6zo5vWZyrmvArUFZA7JSKvNxOuHpp2HlShl2vfdeT7co67EsKFRIPrL4z/dGua41ICul1I0cOiTBePBgqRHsAd4wxKns4S25rr0iIJsb7PnNKty5eE6pLK94cdlj66E5Y28Z4lT2KJYviMioWO6K3E6dfzazsmQV1gdXdHuua4/vDQgMDOTEiRNZOmAZYzhx4gSByRl0lFLps3ixFItISpIVvx66kM8s5fxU6vRsUp77I7cwbUoveiybzJQf3yHkv11uz3Xt8R5y8eLFOXjwIMeOHfN0U1wqMDCQ4sWLe7oZSmVeO3bI8HRwsFQI8uC+Vm8Z4lQZtGcPrFpF6w4duD1fFJZx4gBISqRvrqPcmd1WWfv5+VGmTBlPN0Mp5c2OHZN8zP7+MGeOx5NMJA9xpnS/8mIREZJQyrKkzsHChRAUBC1acGeHljBuKMTH4+vvL7fdzOMBWSmlbiguTmoZHz4sb6KlS3u6RfRsUv6KOWTwznJ+6jIREVK8KD5ebhcpAh9+CJ07S7KT5CyQ4eHyuAxmgUwPDchKKe+2YYMkofjuO7jnHk+3Bri0cEtXWWcSCQkSaJMuXED5+MCrr0qGt8uFhHgkECfTgKyU8m4hIbB3r/RovEhmKOeX7Z05A2+9BX//LfUP/P2lh+zvDw0aXHyY0ykd6J9+kq+HD/dMczUgK6W80+TJkh3quee8LhgnO3lSRtJvuQUKFACHw9MtUhf9+iv873+yZ717d6hV64ohaVMnhDWrJQhPmyYJ3wICoF07z6Xi1oCslPI+8+fL3N5998Gzz3pNpEtKgjVr4PffYd488Fm5gvr8wR80YLVPCIULS3C+5Ra5hqidGEGt6HDK/y+Uwi09NxSarURFyXD05Mlw553w888XpzpMnRA2Bobw008wtSPs2wd+ftC0KQwYIJlYc+f2XNM1ICulvMuKFfDII1Ks4JdfPB6MDx+W4Pv777BgAUSfPM8DhPNJ/tGEEgaWRaLvx8x7YBD3Lh/MuZhcnN2bi6SEJCrFb8QA8b8G0Kf+Ihq8E0LDhlod0iUiIqT3W7u2ZHHr108qVvn7A7KgukcPySfjcEh1yHfflfWC+fN7tunJNCArpbzHpk2yvSk4WKKgB98pN22CUU9HkP+vcMIJ5Uzh2/g1zyvcFfsbfrFnIdoPMGAMfs54mt+xC4o1pEB0NERHc27TFnwOO7EAQzzVVwxnX+MJPBXckRrd6vFsZx8KFfLYy8tawsOhcWOZAPb3l9rMFypPHTwIb7wBU6dC2bJSbbJNG7zyZ68BWSnlPZYskT3GCxd6bN44Lk7W/ywb+CcLnQ3wJRECAuHnefi8tAWefELGNnPlkouH5EVCHTpcXKEbtiGSn4ZNZfzk3vglJZLocBBf5SzPbA7jxcix/NOrBBPfbs/Jph15osU5qpwIx3og9JoVvtklX3aGXueBA/Dkk7KSGuT3sWIF8SH1+fJL2dmUlCSfe/aUss7eynJnyspatWqZtWvXuu18SqlM4vJVNKdPeyxH9dKl8MILYO3awZ85m1Do3D/yDYdDonSfPlc+IXmY9Kp9q3UHLk4xN3K5HLDottOc/XoKOf6cx35KUdQcJoB4rEB/fBYvuiKop7TXeUCbKrYGZU8H/Qy9zgULoH17KcWZmCiR19+fNYMW8fTIEHbskGunL78ET+WfsixrnTGmVmoeqzMZSinPOnpUyvtFRMhtDwTj06ehSxdoWD+B548NYKtfdQr5nJKer8Mhn0NDr31iSIgE6at6tskpNNcHV2RkyGOsD64IwN4YoEMHci+di+O/w5R48WECrXgcJEFcLJt6fEtiohzDHfmyk4NhZFQshktFMsI2RNp2jptJ9+scN05WYxUtKvvUw8M5/WZ/+oYs4u5uIcTHw+zZMHOm54JxWumQtVLKc6KioEkT2LnzUtKGVLCzVxcWBl27wpEj8OUTa3j1x7ehXTt++987zPhlGbdtW8vuSrV4OLAkrVN5zFSl1ixcGP9nOsCkcZjz8VjOJKqtHM28W6IoNmmQbfmyb/Sz8oY6wOl+nXXrygr8YcMwOXLy7fLbeWV4CElJ8MEHsv3Ym4enU6IBWSnlGefOQfPmsHUrzJolW5xSwa7Sh2fOyBbnWT+fp3OZcDqvbELt2vfCW+sJ4xY5R96yzA8pC8CyNJwj1ak1L6RrtMLDMXffw7avl3D/z4OxWoTRvfjzfNn2IcxVQSUt+bJv9rPyhiIZqc4Lnpy549QpmDgRKlSAceOIjoaXn5FEbqGh0nEuW9YtTbedDlkrpdwvPl4qN0VEwPffy9BjKtkxlHv8uCRqKjB9LCdylmTkgYepXWiffLNGjQyfo3WNYAa0qUJwviAsIDhf0PXnRC8Me1sNG1Bp2gfEb97JptvbcfvBvRyc2IA7Vx/h5Yip3BW5Pc35sm/2Oq4X3N1ZJKNnk/IE+V25te2a15mch3roUJg0ScaikZXwNWvClCnSK164MPMGY9AeslLKExwOyJEDxo6V1EhpkNFeXWSk7JB5eufbvOUcgHUOmSM+cuTiZKMdPcf0ptbMW7kE9+yazIolCTTrsJIf/niFQOKI9/Vn5dhp1E/DMW/2OryhSMZN84InJsr4c3JRCIcDs2UrX0e25PXXJUPaokUpT/FnNhqQlVLus3evBL/ixSWDUjoyZGSk9OHevZIQ4vHIz+mVNODSN5KSZLX0hcVZ3lBe8d76fvz80p/4vBuPBfgnxlN91FfQ8SFJL5UKN3sd3lIk44YXL48+Cn/+Cb6+YAzG35++C0P55A8ZWJk0CQoXdmtzXUaHrJVS7rFqFdSpA08/LbfTma4qVUOcKdi6VaapT5+GLl0dMmYdFJTiKur0nsNuvg1D8QkKwPg4cOIg7+qF/Ff6bpzrN6bq+al5Ha1rBLO8dwP2DWzG8t4NvG+f8wsvyATx0qUc/F9/2uVbxKClIQwaBHPnZp1gDLoPWSnlDr/8IskbihWTpP/lMxbY0rrKeu1aeKiJkwqOvxm1uDyVKyN7n1euvG79W0/vz73owl7nM3eFMvqDIzwV8RJrireh9pqR3HrrzZ/uNa8jtYyBUaNkiLp794t3ffkl9Oolu5x+/NGjVRLTJC37kDUgK6Vcxxj4/HNJkXTPPbKa2s1dmqVLoVWzRMaa52nDL/hs3QKlSrm1DXYxBr4bepI3+/hB7tzMeG8DdYPWy17uFC4qMp2YGHjpJRmHbtUKZszgxEmLpo/EsnZZEDluP0Ll9rvo07qcd19UXCYtAVnnkJXyEpmuJ5MacXGyRaVNGxl2DHLfHCxIZ7x9m/P87N+eRtEzJH9iyZJubYOdLAue7l6AWo3hiSfA8epLGFaBZWEFBsrqpswalH/+Wao0HTkiS6b79mX5CovW7RI5cSyA/A9uJfdd+/nvPOna5pYZ6ByyUl7AGzIm2So6WtIZBgXJkPDUqW4PxksHRbCm+Qcsd9Sj0dkZsmXm3Xc9U+jWZpUqyZT8yVpNMFhYxmBiY2WoNzNasEAWbx05Av7+OBs2YuCnPtSvD2fjEyjy5HLy1Nx/8Vdnd8Yyb6EBWSkv4I40iW4zezZUrSpZlAAKFnR7vcG/xkRQq3dD3jUfcGfMGnjnHejWza1tcLWgIHh4WFOc/oEk4oPBgu++I2naL55uWuodPy6f1669eKFkkpL4tlM4ffrIwEqRp5cScOuZa556+ZausA2R1B24mDK951J34OJMeyGrAVkpV4qJgXXrYNs2uX38uAwptmsHr78On30GU6fiu29vik93Z8akDNu7F1q0kGz+QUGSj9IDdu6Emd3D8SdewpSPD+TM6ZG2uFxICL7hi4h9+yPeqbuE//E1jYa1ZO9e5G8uuQKSt0lKgi++kLn8BQtk/jswEKePg7gkfybuD2XUKEnMVbxIylu8krduZaXRJZ1DVspu778vye63bJEgZYz0FidMkN5i/vyyB+f33yV9JNChcWcG5G/DLWdPMH1yT1aXuJPlpaqzp+o9Hn0pqRYWJlV3fH1hyBDpjaZyr6yd/vsPmjVN4pOkbfgE+EEi1y8MkVWEhJA7JIRPDEyaVI8fu0GdKufY79uQoOIFsEaNkt/FdVaTu922bZKzdOVKSZ1aqRJxBYP57olF7JsYzp4SoQybHUK1avLwmyUv8YZ83HbRVdZK2SEyEoIv/PNXqSKF0u+8EypXlo8aNa4tOWOMFFc4eJDfD53n9T+Pk+fkf7y7eBx1/vmLQjGn5XHly0tGq3r13PuaUuPcOel9HjoEb78NH3986efgZtHREFrf8Pxfr9IlcYTskcmb1zuCkBv98w907gxBi2bzTdCrFIk9IHutjYGAAM8t/IqIgI8+gvnz5fcybBi0b8+ChRZdu8Lff8NTT8GIEZA795VPvdGCxzK955JSFLOAfQObufxl3YyuslbKXbZtk0C0YIG8oxQrJgl2UzNnalnSW86fn6ZVIO6WSAbP8+fVVr0IzhPAh+UMDSI3yxto8obTyZPhq6+gUSOpklSnjkd6ouzdCz16wMmTsGSJvO6JE93fjgsSE+Gxx6DRhkF0MSPgzTdh4ECPtceTSpaUmDdyZAvu7NmAGb7NuC9xCRbI3t7LMpK5xfnzEowfflhW3fv4wLffcqhGM3q0l2Hp22+XNjdqlPIhbpTJyxuyqtnGGOO2j5o1axqlsoR//jGmUydjfHyMyZPHmI8+MiY62vXnnT7dmHvvNcbhMAaMyZ3bmFatjImNdf25Dxww5tNPjalTR86dI4cxAwcak5Dg+nPfgNNpzPPPG/MU30q7OnQwJinJo23yFjt3GvP8nStMDIEmCR+TGBBkzIoVxnz+uTE//uja392BA8b06WNM4cLGtGx58W/W6XCYP5t9YnLnNiYgwJgPPsjYn++M9QdNhb6/mVK95lz8qND3NzNj/UH7XksGAGtNKmOkDlkrlVbHjkHp0tIte+UVKVBfqJB72xAVBX/8AfPmwf79Mh8N8Npr0gu67z6oWFGGu9O7oMkY6e2XKCFz3+PHy9xfzZqy/PXppyUn9VXcvZ+6f3/44L1E/ilSm2KVC8rmY39/l53Pm6TmZ52YCLP6RLBzTDizzoSSI/QeZh68i1y7L/xuu3WTVfHr1qV/eP9CNjFCQ2Uh4/DhkgQGZJ64WTPo3h0TH895pz8PmEXkbRLC8OFw220Z/Sl49x5+zdSllCscO3Ypy9Q338j4mrdlfHr6aZgxQyZUk7VvLyUOQcYHixWTQBoTA2fPymsqV06GE8eMkfuOHJFEwfv2yZtr165SQPjUqRu+5qvr74IswLlu6cEMmjgROnWSlz3xi1NYDh+Zn8wG0vqzjouD0aNhwAA4+p+Td2vM5Q3rc/KsD5cH+PjIHPPs2fI7Llv20tTL5QH38oBtDPz2m1ygJSbKhVC5cnD4sOSg7tIFSpXir79gbt8Izs4OZ0uhUJ4eFULbtlliS/hNaUBWym4zZ0LHjpJNKA21ez0iPh527YIdO2D7dukFPfusbIHJmfParTDdu8sWlJiYS73poCB5823TRlIYpjLdZd2Bi1OczwvOF8Ty3g0y9rquMn8+dH14H1/cOojG277EP0+grcf3dun9WZ87ByNHwqBBcOIEzC79Ks32j8DCyOKvDh0kq1quXLJAsUgRuThzOiXgvveerBs4cEA+YmIuHdzhkO18/ftz7Gwg338vF00bN8pSh1dekSRcVy/aysp0UZdSdpowAZ5/HmrXlg9v5+9/aXX35Xx9ZeHZ9u3Sg8mZE/LkkR4NSBA+flzeLdM55GtHHeHU2LIFvm71G6t5mnynz2P99wY380/xAAAgAElEQVTkud3Wc7hDRoZa0/uzzplTUot36SILnYcN7EADxuFPPE7Ln00l23Ln0PvJsfsv2LSJpF9/w3HhIi4x7jz/LVxG8KmjMiXStKn0kkeOhKQkjL8/Swq04YvHA/n1V+k016x5cUG122d2MhsNyErdyODBUhy9cWOpWJQrl6dblH6WJUOR1xtytiyZK84Ad6x4PXoU+j8YzrS45vjgxEoIkAuJ2zNXQL56yDk5oQWkLkdzRn/WuXNLArOuXUOY1nsRcb+HMyUylGUfh2BZUK0aFK8YjXXvNKYu64KvM5EEhy+vF27KA68+Ss2it3LqlCy093c8ht/ycEZuC2Xe2yHceqsMvDzzzLXXher6dMhaqetZtOhCNfvHpfpMNlkolBGunkM+fx4aNjCMjKhOVfOX3OlwyMquPn0yfHx3yujwvit+1rGxkq9j6VIZlV6yLAlnooM6RPCAYxF/JDVkJSkv+goIkNmNZ56R61ffy7p73rzoytVsHbK2LGs80Bw4aoypfOG+AsBPQGlgP/CYMeZUehuslFdq0ACmTYNHHpE3fXVTyW+yrnjzNUbWCR1asY8KQfsg0e/SvGYmzMSV0eF9V/ysg4LggQfkA6D0m/OJO5KXHf8UZGtMR3wCEsgftBVHYAITulQnf34oUICLn1O6Zs3oSEB2kpoh64nAcGDSZff1BhYZYwZaltX7wu1e9jdPKTeLi4OXX4Y33pBMW+3aebpFmc6NkjhkxMCBstbogw/K4v/iLklOsmRJps3EZcfwvqt+1smCCwUQ6XuKwOJX9reC8wXRvHnqjpGVUlu62k3TCRljlgInr7q7FfDtha+/BVrb3C6l3O/sWckmNGGCjNsprzF9Ovz69jJ+rPoJ7/Y1krns3ntlmDoTBmOQHM1BfleOvFyeo9kb2NFGdy30ywrSW+2piDHmMMCFz7fY1ySlPMDplGWgS5dKesrnnvN0i9QF69dDv467me37CI/GTcI6F33zJ2UCrWsEM6BNFYLzBWEhvU5X7ddOLzvaeL0ev9eltjx5UhZuepDLV1lblvUi8CJAyZIlXX06pdLn/fdlr+WIEbLfWHmFQ4fgyWanmJ3YjDx5wGfunCy1idXVQ852yGgbb1atyaOOHZMJcF9f+OQTKYe6a5fHVuynt4f8n2VZRQEufD56vQcaY8YYY2oZY2oVTmVyAaXcKikJVq+WEjkvveTp1qgLYmKgTfN4Rh1rSxlrPz5hM1KdZzGrFKzPCrxuJODwYdk33bChTH0sXSr3d+0Ka9bYk8szndLbQ54FPAMMvPB5pm0tUsrdHA7pHScmZo9cfpmA0wkfN4/guQ3jqeezFJ/xE1NdflJX9XofrxgJOHRI9mQtWiRL9suXl0ptZcvK98uUubZEqpvdtIdsWdYPQARQ3rKsg5ZlPYcE4kaWZf0NNLpwW6nMJSpKkiAfPixBOSDA0y1SyHvlsPYRvPNHQ56zJuDj73cpm1gq3GhVr8qGnE75XLCgbLR+7z1J9bZ9u+xfL13ao8273E17yMaY9tf5VkOb26KU+zid8OSTUi3pxRehaFFPt0hdMHgwlJz6JYGcx8c4IYE01fDVVb0KkBGvUaNg7FgpjpEzJyxb5tWjYJo6U2VPyYu4hg+XUoUq1VyZdWnSJPij12/M5hcsH8BypDnxh9sK1p8/f2lUZe1a2LBBKij4+Umb/fyk7KCfn73nVTe3aJGUIt26VaqyRUVJQPbiYAwakFV2NGOGDFV16iRJQFSquXJ+9vffYWznCBb4tMWnalWswZ/KIps0Jv5w6areo0flQm7WLCk3tW2b5AafPRs+/PDax585IwF5yBApedS6tRRkyMw50b3ZuXMyT/zLLzIfHBYGLVt6fSBOprmsVfbidErFJl9fyfIUaF/JvuyQr9dV5RVXr4aX629lUUI9cpUqiGPFn1L2L51s/1389ZeswI+IkEnuEiXkjb5nTwnI0dFw+rSUvkxIkI/4eKnQ4OMjwXrYMKl3GBgovbb27eVD2ccYaNFCksb06GHr/3d6aT1kpW4kKkoWd9g4b+zqogp2yWigKtN7Lim9Y1jAvoHN0tWmXbugbl14O/FDuvmPwrFyhcdXuwKwfLm8wd93H/z3nww/t2ghFRSqVUt7rysxEf78U0ZowsLgrrvka4B16+R2JunJeZWoKEl1+/77cqFkjFf9HNMSkNO7D1mpzMUYSYkZFwf58tm+iCszrOxNvmiIjIrFcGm4OS17dO3OunT4MDRpIu+fLVa/i+OvjZ4PxqdOyUK/++6Dfv3kviJFZJ64Xz+oXj19b/i+vjL8PnQo7N8PEyfK/Zs3y6hNvXpyEZDJeHTP95o1ciEzadKln50XBeO00oCssofJkyXxx6RJN39sOmSGlb12XDTYmX/59Glo0ziaof+2YfFXW7ntditDw9QZZgz88ANUqADjx8Obb8pcsStYFuTNK19XrAhffy3FMu67T+aZt21zzXltZsdFXroYA198IUMrSUmS3OOJJ1x7TjfQgKyyvoMH4dVX5c3ORTmqM0O+XjsuGuzKuhQXB+1axvPh1jY0N7OonHNfmp7vErNnQ4cOsi917VrZf5Uzp+vP6+srPfK//4aPP4Y//oA6daTYiZfz2MjQF1/IHPHDD8vq9kxaYORquspaZW3GSBBOTJQhQhfVNfbqfL0X2LUdKKNZl06fhndClzNk40tUYzOMn0Cqa/nZLT5eeqPVq0sbfvgBHn3UM/Wvc+aUzFEvvihDsblzy9/v11/LhUJyj9qLuH1kKDFRLmCefx7y5JH/7Uw8RH017SGrrG30aNmeMmRImrI9pZXX5etNgTeU+ztyBHrUXMKXG+tLMPbzkxSGnrBrl8w/NmggVwk+PjLs6YlgfLlCheChh+Tr9etla17Fiq4bPs8At40MGSMjFnXqyILMPHkkKGehYAwakFVWV7++DG39738uP1XrGsEs792AfQObsbx3A68KxuD5i4bdu2XKr9SBpTi4MJLgdEoWLndbs0Ya899/sq7AC3ufANSsKW0tXFhWd3foAMePe7pVF7nlIu/0aWjTBt56Sxb8JSbad2wvo9ueVNbkZVsfsrv16+GZxoeJdQYw89Od3NmtoQwX+/tLViV3zgHOmwdt28Itt8jXHiq1lybx8TBgAHz0Edx5p8ybesnft0v33//1l/yu9u+XHvJrr3nN604t3Yes1GefSc9i4kSvSA6QnS1aBD1b7mTm+SYUCKlAzmW/S4KN8PA0Z+GyRefOEtB++03K72Umf/0lyUUeeECSj5w6JRcWWZExkuDjwAGYOjXTprjVgKyyt23bZG7woYdg+vRUXVFnhyxbnjB1KnzVcSWzTHPy5Hfg+P1XGYb1hNOnZWg6IeHSPGRm9skncuE5bJgMZWeynuN1xcXJsHSuXLBvH+TI4dntcBmkiUFU9pWQILlsc+eW1ampDMYe2UuZxY0YAd89PocFzgbkKZlPMnB5Ihg7nZLJqVYtOHlSFpJl9mAM8MgjcMcdUrWsVSvZ3pfZ7d8vCVKef15ulymTqYNxWmlAVlnLwIGyh3TUqFT/I2eGLFuZydmz8Fm7CA6/8hHjA7vgX/1OCcYuXOV+XfHxUvP6889lxCRfPve3wVUqVpRUnEOGwMKFMrf8yy+eblX6/Xph9OTvv7Ntjm/dh6yyjnPnYORIGb5r1y7VT3PlXkqnE3bulOnsNWsgZlEE5Q+Hs+2WUI6UCaFQIVlAW6gQF7++9VYZcc9sU9/GSIrm719YzLcnmhNgxeODH9an33pmnjM6WhYEzZ8vw7u9e2edYd1kDof0/lu3luIXJUt6ukVpd/q07IQYPx6qVpWLittu83SrPEIDsso6cuaUxTr+/ml6mp31c48elSx+a9ZIBaN16y4lXGoV8DvT4lvhMAkknfZjeOxAwpPq8ceZOzgULUOodYgglHD6+IXiuC+EBg1k/U7t2ql7WZ6aCz9wQJKhMXsWk32eJog4LGMgAflBNGzo8jZco0cPWVE2fryU2szKypWTC49k3bvLCNGbb3p/PebYWJgzB/r0gffey3xXojbSRV0qa1izRoa7fNI+C2NHpaazZ2HQIPhzcAT3xi/mH0dZypZMon7u9Zxv9gilOtaj4pf/w+ebMdc+efZszjdqTuwnX5D3wx4AOC1fvi7yHh8deYEj3ErOnLLINDlA33XXtfkrPFFxKjFRaiWM6fsPQxK60SJpJqZMWaxDkfJNT2xrSnbiBKxaJekVs5OkJElw8vPPUpXqm29k/tybnDwpazx695b/2bNnZd1HFqSrrFX2snUr1KgBffvKFXY6pLdnmZQkO6v69oUaR+Yyx2qJZZxcHBgNDJT5y5dekm0+TZtKoPLzk0hWtCjcc48M6b7yigy5X/U/ufCzTczcV5Xt8/5h598WxTlIsxzhOOuHUrFTCI0by+JhV9Uqvp5VqyTfSulNYfzgeJIAf4PPB+9L72zt2gxta0p3Tz+5qlfHjhAQkObzZikzZkDXrpL8pEcP+d/whqA3cyZ06QLHjskceJ06nm6RS2lAVtlHUpJkXNqzR7Y7FS7stlMvXizvczs2xXFXSCA/Vf2YEqP7yjd9fCTAfvaZ5N5NdqP9txERMrSbnDDjq69kfq1bNznGK6/AiBEYLAyQgB/NmMMS30bcfz9sYBtB5Y7iV+DcFYfNSK3iq8XFyeteOzyC878tZnOhBnT9sAiNF/fCGjIESpXK8Dky1NN//3344AMYMwZeeCHDbcn0oqIkw9XkyVJN6tZbPZc05/hx+Vv+4QfpuU+YIBfSWZwGZJV9DB0qPbIpU2Qxlxvs3Ak9e8KK2ccZlPtjnvD/hRz7tmFt2XxlQE3PUO2NAvb27XLiuXMv3pWYMy/vvnKKOXMtTm05yCGKUTf3HzTMNY+IklXYXLE0pW5LJOKdB9L9eo8fl1POmgVLfo/l2ZjhDKSPpL8MCsJK4+u8We833T39UaMk73OnTjBuXNZbwJURhw/LaIwxMoRfvbpcTbrxApbGjaWS1bvvylB1Gtd6ZFYakFX2sG8fVK4swWvOHJe/AcfFwdfPRHB62nxKOA7R0fEj/gnRWJ07y3arggVdn4Hq8l60r6+8ub3zDhhD9C3FiY+KJk/iOSycJOBPY+ax0v9+alS3qFmTix8VKsjhEhNlkCEx8cqvT5+WNUIzZ0rd94bO+bwe+DUNEucRkBiDQXreOBzQv78syEmF1PR+y/SeS0rvSjfs6U+bBo8/LhWbpk+/clRCXRITI5nKpk6FoCAZOn7zTQnWdjJGhqO//hqGD4f8+eUCtXBhWUmdjWhAVtnDunUyLDlzJpQo4dJTRUdDnwci+HRtAwKJk2B0//3yhlOxokvPfY2Ugn5iInz/Pac/Gkiev7dfnMPeV/JuRjy6ik1r4qm96itWna+OH/HUYCPhhLISeb6FkxzE4E88AZwnlD94kinMLtONIk814fkzn1N82udYrVpJ7ue3307XSEBqer9p7iGfOyerjG+7Ta4icuRIVVuytR07ZCvY99/Lxctvv8lqwYyKiZFjDh8OmzbJvu+wMCnykk1pQFbZhxvmw06elFG+hqsH0N96Fx9nkswRf/RRqnuGbnN5D9rhkAuGTp1gyxaoUuXiw5L/65c3ep81D/Wj5IFltB16/7XH8/eX4H/XXfJ18s86nSMBqen9pmsOeds26eXlz5/qtihk7cXQoTLCkyOHBOm5c2WhVUiIfC5ePOXnJiTIUHh8vFwMnTwpF2snT0ov+NVXZRopm18gpSUg67iOynyOHJEFT++84/J/9sOHoVXDaHrsfJFy3Zvj87X/pZ5haKhLz50uISHSY706WFauLKtae/WCCRNkj7Blcd/tR7nvdSCyLBQfLK9r8WKYPVuymiQlybGuDrohIekakk/Nnu/koHvTVda7d8tUxWuvQaVKaW6LQkYWhg27dLtwYbnoGjFCdgeAzG9s2yb3d+sGK1dKms4jR+SCuEED+ZsrUEDWc4SGyh49ncNPM+0hq8zn0UclYGza5NLi9vv2weMNjjHyn2bUZB3WpElQtqznqhTZ4eqV3CkNN6fmMelk217p48el93b6tPT+s1G+Y7eIj5f/r5Ur4cwZufgFyTF98CAEB0vPuXhx+R+8P4XRFQXokLXKysLCJKn+xx/LPKaLbNsGzz2wl8nHm1DaLxLHtJ+gRQuXnc+tUjPc7MLFaRnOJnb+PDz4oCSDCQ/P8vtYVeamAVllTVFRMjR5yy3yZuyilIBr1sArjXcx50w98udOxPe3OZmzN5wVGSPVvL77Dn78UVZWK+XFtPyiypr69JGsQ+PGuSwYh4fLlNjpvCUJbNEY35XLNRh7k7VrZc/5hx9qMFZZji7qUplHt26S2cdFNXXXfBXBoe4jeKL4U7y/vAm5g79zyXlUBtSuLUMY2SDDk8p+dMhaeb/ERJcnejg8PYKCbevjTwLG1xdr6VLtGXuTVaukWER2KxShMj0dslZZS9euUrDc6XTJ4WNjYcULE/AjAUC2BIWHu+RcKh0OHICWLWVLTXy8p1ujlMtoQFbebf58KRRQokS6SivejDHwTod9NDz5k+ybdDi8d49xdnTmjKTDPH9eMrJlk/zHKnvSOWTlvU6fln2PFSrIIh4XGDECGoe9RECgD9a3P0rmosy6xzirSUyUur7bt8Pvv7s/RalSbqYBWXmvN96AyEhYsULqCtts2TJ4/XXo0GgCjd/dC/Xq2n4OlQHTpkmO5dGjZd+xUlmcBmTlnY4dkyQgb70F99xj++EPHYJRLX+jXOnGDJtWFJ+8Nle7URn3xBNSv9eOogdKZQIakJV3KlwYtm6VajE2i4+H0fW/5/uojhx++Svy5n3F9nO4SoazXGUGM2bINEXFihqMVbaii7qU9/n9dylqUKQIBATYfvgv2q+mz+7OHKt4P0X7vWj78V0lOQ90ZFQsBoiMiqXP9M2EbYj0dNPss3ixJPxwYVpUpbyVBmTlXebMgYcegrFjXXL4nz6P5KnprYnJW5TCS3/JVKt2B8/beUVRBoDYhCQGz9vpoRbZbNMmaN0a7rgDJkzwdGuUcjsdslbe4+RJePFFqdvbqZPth98+bgV3v9GRAj5R+IavhkKFbD+HKx1KoWzhje7PVA4ckAuxvHllhMQFUxVKeTsNyMp7vPaaLOaaM8f2oerYxRGUfuFB/DmPj68vVuxZW4/vDqmpJZxp9e8vGVr+/FNK+imVDemQtfIOM2bA5MlSd/Wuu2w//NJec/Az8ThwYiUlZcpMXD2blCfIz3HFfUF+Dno2cV1NaLcZPhyWLIE77/R0S5TyGO0hK+9QtCi0auWSxTwrpx+iztqvJNOXRabNxJW8mtrVq6zdtpI7KUl6xq+9BvnzQ9Wq9p9DqUxEi0sozzJGUla6SMw5w9rCD1E7binWhAkEHtqrmbhuIHkl9+WLx4L8HAxoU8XeoGyM5CgfNQq+/Raeftq+YyvlRdJSXEJ7yMqzXn9d8kcPGeKSwDy3+SgejZ3Hru4jueMZrZ97MzdayW1bQDYG3nxTgvFbb2kwVuoCnUNWnvPLLzB0qAxduiAYr52yk2bhb7KtZFPu+LyL7cfPily+kjsxEZ57Dj7/HF55BQYMsOe4SmUBGpCVZ+zeDZ07w913w6ef2n74mBgY23MXJ3yLUGrReJcOi2cl11uxbdtK7pMnJflHv34wbJhLKngplVnpkLVyv7g4ePRRGaqeOtUlyTneeQfGHG5BhwVNKXGbny3HzA5pK3s2KZ/iHHKGV3JHR0NQENxyiyQAyZs3gy1VKuvRgKzcb9062LULfvoJSpWy5ZCXB8s6u46Te8Z5ur78AvUftC8YXx6oktNWAlkqKLtkJffx4/Dww1CrFowcqcFYqevQgKzcr25d2LdPeks2uDxYBsTE88HMjwi04ljVOhS4w5ZzuGWxk5doXSPYvtd08CA0bgx798K779pzTKWyKJ3AUe6zfTt89518bVMwhiuDZc8ff+E25x66h/bhy3UHbTtHlk5b6Sq7dsnF18GDMG8etGjh6RYp5dW0h6zc49w5mTc+elTemG3MVZwcFJ+Z/zvPH5vC1PwtWH93KSwbg2WWTlvpCvHx0LSppMMMD3dJ9jWlshrtISvXMwa6dIFt22DKFNsLBxTLF8Q9ezfTb8NwDNDi7Dzuitxua7C0LW1lQgIcOQJnz8p2r6wmJkZel78/jBkDy5ZpMFYqlTQgK9cyRpI/TJ4MH34IjRrZfoqeTcpTefVxDD5YgF9SIvdFbr0iWIZtiKTuwMWU6T2XugMXp7mGcOsawQxoU4XgfEFYQHC+oLRlr/r5Z6hfXxY0FS0KefKAr68EMJAUkpUqQe3a8jPq3RtmzcpcQXvePKhcWRJ+ADz4IJTPAnm2lXITHbJWrrVmjWTh6tpV9iK5QEVHIT472IY3rBH4E0+irx81n36E+heCpV0rpFO12GnfPpg/HyIiYMUKmDsXbr9dtv3ExUl5ydtvl6Hcc+dkKxBIhaM775T7jh6VxBnffCPVrwDGjZOkGvfeK4/zpv27R49Cjx4y+lGhAlSv7ukWKZUpaS5r5XqLF0v+aBcEERN9jn9vuYthvMbb02pQ4K/wa3JV1x24OMX53+B8QSzv3cCehhw/Lhcc48ZJr7ZwYQmeH30kvca0io2V4F6pktyuW1cCPEjv+oEH4Ikn5MOTZsyA55+XIfi334Y+fWwvnalUZqa5rJXn/fyzrKS+/35oYFPQS8GOdn2pGLuL+9+sQoFmIdDs2qIRblkhbYykAn35ZXj1Vbjttmuyg6UpsUhQ0KVgDFIneM8e6XkvWybDw4UKSUA2RgJ//fpyIeJnz97rG0pIkPPkyyftHD36yvYqpdJMe8jKfslbXEJD5WsXpa08PW8luZvey4wiL/HIoRHX7YC7pIccHy9BaN48mD1bXmN0NOTKleLDba+iZIz0onPkkD2+5cvLkHbevDIH/eCD0KyZDIXbISkJVq+W1zpnjgT+0aMvtUVTkyqVorT0kL1oIkplCStWQJs2Ms85darr3qjPn+fcE89xkOKUnz7ghqPhtq2QBnA6JcNYpUrQrZssyjp1Sr53nWAMN04ski6WJcEYoGxZOHFCeujt2snvoEsXWL9evr9hg8zx/vSTDIOn9SK8d29ZiHbvvZJ3vGBBuOeeK9uilMowHbJW9vnrL+mVBQfD77/bvr3pcltGL6d81C4mtZnJc/fmueFjbUsH+e+/EvBWr4aqVeG336BJk1QFJJcPm+fJIxdCbdpIwL08E9rWrbLy+Ysv5HahQjKv/csvUKCA1CP+/vtLr8OyZB54xgz52hjpcbdoIXuL8+e3p81KqStoQFb2GTcOcuaEBQugSJFUPy2tRRsSEqD92AbkKraHhZNKpuoctqSDzJdPXt+330LHjlIcI5XcmljEsqTXnOzJJ+Hxx2HzZrmYWL1asmglO38eTp+Wr42RD4cDDh+GYsVg0CD726iUukaG5pAty3odeB4wwGagkzEm7nqP1znkLCopSd7AnU55Ew9OfeBL89xqUhLfv7KCjl/XY+ZMaNnSjhdwE4sWQZ06EozTOV9q+xyyUipTcMscsmVZwUA3oJYxpjLgADy8B0O53ezZUK0a/PefbGtKQzCGtM+tnnx/GB2+vp9e90e4JxgPHy7FEfr3l9vpnC/NcGIRpVSWl9Eha18gyLKsBCAHcCjjTVKZgtMpW2369ZPUiAkJ6TpMWuZWzZ695PzkHeY6WtD1uzrpOl+qOZ3Qs6ck6GjVypZKRbZWUVJKZTnp7iEbYyKBIcA/wGHgtDFmvl0NU17szBlZPNSvHzz1lOyRTef2muvNoV5z/4oVnKvbiESnReTbIylR0oUre2Nj4bHHJBh36yaLn3LmdN35lFKKjA1Z5wdaAWWAYkBOy7KeTOFxL1qWtdayrLXHktMAqsytVy/Zi/rll7LAKSj9C5NStSUpIgITGkqu//bibyXQudG/6T5fqvz3n1xkfPklDB2apsVbSimVXhnZh/wgsM8Yc8wYkwBMB+69+kHGmDHGmFrGmFqFCxfOwOmUxyUmyuePP5Z0mK+9luE9qKmaWw0Px1yYZ/a1nPj+GZ6hc15X8krj0qVlFfJrr7nmPEoplYKMzCH/A9SxLCsHEAs0BHQJdVaUPF+8YIGsOC5QQFJi2uRmc6tT4qryCAH4EU+Sw5dVRe+kvm1nv+DUKVlJ3a6dXHDkufHeZqWUsltG5pBXAT8D65EtTz7AGJvapbyBMZL8olYtmS8uU0aCsxut+ngky7/YQ+PAX/n83qfo8PhHdNnjn+byiTcUHy9z4vv3S+ILpZTyAM1lrVJ26JAULli2TAJx//7QoYN70yQePszpkhXYnFiFdq0+J7DCpTUItlVqMgaeew4mTIDvvpMkGkopZROt9qTS78wZGa4tXFiC1YgRUl7P39+97TCGc8+8hH9iPC+XGkRA+SsXBNqWcnLgQAnG772nwVgp5VEakJXYtUuGpZcuhb//lsIFy5Z5rDnmhx/JuWAmPR2DOPlwHL5XdcxtSzlZoYL0kN9/357jKaVUOmlAzs5OnYKwMNnCNHOmFBR4/XW3zxNf4+xZ4v/3Khu4h6juz5E7aA2xl+UdSXelpsvFxMhFxyOPyIdSSnmYll/MTpxOKSywY4fc3r0bOneGlSslAcbevbKa+gZlBN3hZEJunnFM5vPKE/h6UEH7U07u3w933AE//GBXk5VSKsO0h5zVrV4NO3fCwoWyYvrYMamVO2oU1KwptXKrVfOemrYxMfTqlYOfo5uybrLk5LA15eTp09C8OZw7BzVq2HNMpZSygQbkzO7vv6X27d69sGePfA4OhmHD5PuPPy49wgIFZEvPww9LDV+QYhDVq3us6dc4fpzzlarDsX706PkC1arZfPyEBEmJuXMnzJsn88dKKeUlNCB7uzlzICICjhy59JEnD/zxh3z/hRdgyRL52t9f6nGoaIEAACAASURBVOBenhHthx8gb14ZovXmFJARETj/9xI+x47wb7E6fNnPBed44w2YPx+++QYa2LBlSimlbKQB2RscPw6bNl362LdPgqxlwdSpElSLFJGPW2+F22679NwBA6TnV6aM9Ix9rloWUMfFVZHsEBEBoaH4xMeTiC/vvxltfy0HY6BgQXjzTVlVrZRSXkYDsqe9/bYE1WTFismc7rlzsrjqq69kn+z1erchIe5ppyvNnImJj8cCHJahTlw4YPPrsizZ1qWUUl5KV1m72/r1Msy8YYPcfuIJ+PRTyRN99ChERsKvv15a6Zw3r3cPNdsgKTAHAIk48Anwh9BQ+w5uDHTqJEPVSinlxbSH7A6xsTL0PGoUrFol5Qrr1ZNVvlWrykc29rHPeyyjDkMeX0e110Lt7fUPHw4TJ0Lt2tC4sX3HVUopm2lAdrWkJLjzTpkXrlBB6us+/TTky+fplnnen3+ya3UUH37YnMc7NKbaFJsD5pYt0LMnNGsGL71k77GVUspmGpBd5d9/oXhxGW5+912psRsa6j37fT3t+HGcjz+B77GclCzahBEj/Ow9flycFMPImxfGjdOfu1LK6+kcst2MgTFjZJvR5MlyX6dO8MADGhSSXZjXTTpyjHYJPzJukp/9AwaTJ8PmzbIgrkgRmw+ulFL20x6ync6cgRdfhJ9+gkaNdM7yeoYOhTlz6MEwGrxRgwcecME5nnsOypeXuXqllMoENCDbZd26S1mxPvkEevW6dk+wgj17MG+9xbyAliy5/RXWfGzz8Y8fh7NnZV+2BmOlVCaiAdku//4L8fGS0KNuXU+3xmuZMmUZWvkbBm1pxvzvLQIC7Dy4kdrNK1ZIClEPF8lQSqm00ICcEVFREoBbtYLWrSVHdJBNdXqzmhUrYPZs5lgteX3D0wwZAlWq2HyOb76RMpKffabBWCmV6VjGGLedrFatWmbt2rVuO59LxcZKAF63Dv75R9IyqhQtmTiLe59vi29SIucJ4PkK85i0tb69I/p790qEv/deKRyh0wVKKS9gWdY6Y0yt1DxW37XSIykJOnaEP/+UVbwajK8rbEMk/w0dhW9SoqTGJJHq+Scza1OkfScxRhbTORwwfrwGY6VUpqRD1mllDHTtCjNmyGrhxx7zdIs8LmxDJIPn7eRQVCzF8gXRs0n5i/WLf/n2N77eGo7BIhEfEnx8WVGuPDPm7bSvxvH587KIq107KFHCnmMqpZSbaUBOq3nzYPRoWUXdrZunW+NxYRsi6TN9M7EJSQBERsXSZ/pmAFrf6sPgcb0448hNx4QpVLtlDZsbF2R9cEWsqFj7GhEYCGPH2nc8pZTyAB3bS6smTWD27CsrNGVjg+ftvBiMk8UmJDF43k4oXJgF5ZrQJHEevxVtwIQn72d9cEUAiuWzYfGbMdC7N2SVdQlKqWxNA3Jq/fab5Ea2LGjeXLNuXXAohZ5uUHwc8QcPcfCIL90Pfcf2HBW4pc1afPyc8n0/Bz2blM/4yX/5BQYNgsWLM34spZTyMA3IqbFiBbRpA2+84emWeJ2re7q+SYkMnzWI6d/3om2zOBLiHAwZe4aSxX2wgOB8QQxoUyXj88cnTshcfs2a0KNHxo6llFJeQOeQb2b7dukRlyhxKTe1uqhnk/KX5pCN4ZN5w2m4Zw0DbhvG2i2BzJ0LTZsW4TVszifdowecPCl1jn31z1gplfnpO9mNHDwoc8YBAbKYq3BhT7fI6yT3dAfP20mHmV/z2OaF/FC5D29veZURI6BpUxec9I8/YNIk6NsXqlVzwQmUUsr9NCDfyAcfSDaupUtlW41KUesawbSeMRpWTmNvlVZ02Pwx3bvDyy+76IR168KwYbL3WCmlsgjN1HUjCQmwY4cLcjxmMRER0LAhJjaOWAJ5r+4iBi0JweFwwbni48Hf3wUHVkop+2mmrozatUvmJ/38NBjfyNmz0KcPLFiAOR+PhcGfeD5+MNw1wfjPP6FcOdi0yQUHV0opz9KAfLW4OHjkEXjoIdnnqlK2ezfUqQODB3MsLjdxxp8EHPgE+hPQJNT+88XFSSUnh0OCslJKZTE6h3y1vn1h2zZZxKV7jVP2++/Qvj04HGweMo/Q/g25O08dxj0VTrEOoRASYv85P/wQdu6U34tWclJKZUEakC+3ZAl8/rmsRmrc2NOt8U7jx0tPtWpVZnWawaNvlaF0aRj+awjFyrkgEANs3AiffgrPPqu/F6VUlqUBOdmZM/KGX66cvPmrSyIiIDwcQkPhvvswnZ9jcPCX9Oqek/vvlzobBQq48PyTJkGhQlLnWCmlsigNyMnOn4fKleHttyFnTk+3xjY3qsSUKhER8MADsro5MJCE3xfxYtJYJn4ITz4J33wj27Rd6rPPoHt3F0d9pZTyLA3IyQoXlqIRWcgNKzGlJihv3gwvvCAXK4CJj+fbTuFM3BvC++/De++5eJp9715ZxFWqFJQs6cITKaWU5+kq62PH4NFH4cABT7fEdjesxHQjR45IneeqVWHfPvD1xTgcxDn9+fZAKJMmQb9+Lg7GTqdMIdSrJ/vBlVIqi8vePWRj4H//g7lzpbt3mQwP9XqBlCoxXX1/2IZIfh0zndu2reXAHdVo9PLjtL4jH2zYAH37Yrq/zp/jdrL0w3CWWKH0nx1CaKgbGj96NCxbJovI/PzccEKllPKs7B2QJ0+WFUmffnpFApAMD/V6iWL5gohMISgnV2gK2xDJtC9/YPyUt/FLSoClFk/FO6F7e1rv2MHaDQ7ebAtLloRQsWII06dDhQpuaPg//8Bbb8GDD0ovWSmlsoHsO2R99Ch06wb33XdN+b50D/VeJWxDJHUHLqZM77nUHbiYsA2RGW52WvRsUp4gvytTZl2sRbx0KQU7Psa4KW8TkJRw4Q/BUHPvBvr/eICOTzuoXVu2ZI8cKcmx3BKMjYEuXWTIeswY3QuulMo2sm8PeeBAiI6GsWO5Os9jaoZ6b8YbetmtawTjiI1h7dAJhGz4gwLOeM598BEP1AiGWesocXgfS8vcRejedfgYJwkOP3490Zr1n93DNj9ZcN6rF+TJ45bmivh4KXU5YIAW9FBKZSvZNyC//z40aJBit+9mQ72pcaNednJATtU89eV7gK/OgBUbK/PfCxZAsWIQFCTpLO+/H/bsgbvvpsXJk7S4/Dn/bgIeghYt6NgrF5FRsdT4dwc1Vv3LvH/bsGJ7PQrfdYR1YUUpUSLVL9c+AQEyf6yUUtlM9gvITickJUm3r3nzFB/Ss0n5K3q3cNlQbyrdrJedqh70/PnQsqX0Gh0OuYBo21bKDkZHQ+7c157g3XclIBcpIiul9+6FhQvldTscF/Nznz5jEUJ1xs6JY9aeRoTF+RNY6jilG63gi5dLXwzGbl3c9u670KoV1EpVYRSllMpSst8c8pQpUtT+0KHrPqR1jWAGtKlCcL4gLCA4XxAD2lRJUyC6Xm86+f7kHvRdkdt5OWIqtf7dStmDu5j/9TR5oNMpwfj8eQmiiYmwejUcPizfz/X/9u48Pqrq7uP45xBWRQiIgsEFEUuV2oqmarRqbEVcC9pKQXErKBbF5VEESkUR64pLa60+Ktrqo8hLigoIyhoERQUFFEQlpKIsBpTNBYEk5/njN4EQMskwubPcyff9es0rw8yde8/hzsxvznJ/pyl07bpzjDUry65FuuOOnc8/9pj1BDRqBFlZlDVoyJjifE4/PZL4akhLWNWafY/8htY93qVzvw95qH+7XVrwQ8Z9xKqNW/Ds/NGQkLHw8ePhzjth8uTg9y0iEgJ1az3kb7+Fn/zEkkzMnQv1Evd7pHILGKyVXR7YDx38Gqcun8+T40ZQv8y2ccBHrQ/jqK8K7QXDhtkM8JISWwN4+vRdu60j6xDvWCO4wvPr19tSzp98At9NnUvZjALGrM3nHfI44giL9eedZz3c0ZZKPOmeGVV23bfNbsJbg38dyP8TAMXFNss9J8d+dGi9YxHJEHuyHnLd6rK+805LevHqqwkNxrCz2zlad29OdhOGFDxNg0gwLgMm/PRknj2vH2PLd3LHHbYMZJQx5O25eax5Zjo/TCpgQfN8Zj6TxyeDLAivW7dzu4YN8zjxxDx6DIbnzoMOHWKrQxCT22rkPVxxhf1YeuEFBWMRqbPqTkD+9FN46CH78j/uuKQcsnvntrt2c7/9NlwwAO6/n4FdO/LkOz25c+JD1C8rZXtWfUafcD69e+Xvso/S4/JY1iKPwkIofBiWLbOliAsLLblYaWkeYIG6VSubo9atm/0tv7VrF70VXJ0gJrfVaPRo66Z+9FE48sjg9isiEjJ1JyA/+qjNQr777uQed84cG8tdtAiWLIEWLWDpUrqfey7cdh3XHXAgHT6eT+GRufzhqgvodnRbPvvM5mFNmwYzZ8LGjTt317w5HH64/aa46CJr7XboAB07WkAOUhCT22p04YU2ya537+D2KSISQnVnDLmkxALiL36RvGPOmGHZpry3yVc33mjd0JVWkyoutuHf8iD85Zf2+CGHQJculs65Y0cLvC1bJjdXRsJmWW/ZAt9/H/yvCBGRNKIx5Iq2boUffrCWaTKDMcC77+68X6+eBZ9IMPbegvDw4daIBivir39tCTm6dIH27asPvsm4JGm3bvegDBoE48bZj6TmzYPfv4hIyGT+ZU8PP2wzq7/6KjnH++EH6NPH+prz86FxYxvAbdiQ8lUZZs+2JYa7dIHPP4e//hXmzbOJWGPHWubIww6rORgn7ZKkoE2aBI88Yt3VCsYiIkCmt5BXr4YRI6zbuE2bxB9v+XJL3LFokU1QuukmawZHZkm/l5XHrV0t30ebNvD3v9tyw40b7/mhYskElpaKi21i3c9/nvzxfBGRNJbZAXnQIBs7fvDBxB9r/Hi49FLrmn7tNTj7bHs8L4+FTfIYNgwmTLBe6/vvh/79Ya+94j9cUi5JClr5JU6bN9v4ejy/REREMlTmBuS5c215xaFDbTA2UccoKIDsbIuwxx5rfc7t2gGwYQP86U8wZoxtcuedtsBUVRkv91RSLkkK2vffWz/8yJHQqVOqSyMiklYyNyBPmGD9woMHJ2b/lbNkDRxoM6gjrb7CQkuVXVQEf/mL9V5nZwd3+KRckhS0pk1h4sRUl0JEJC1l7qSuu+6ChQstCCTC1Kl26U5pqQXlFi12BOPZsy0l5bp1dhnTiBHBBmMIJt920mzYAJdfbjPYnNMaxyIiVci8FvL27bBypa2l27p1Yo5RWmqTtcDGjCvMoH7uOejb13qtX3st9jSV8UjYJUlB2roVzj/fspRdfvmO7nwREdlV5rWQR42yLBpLliRm/97DDTfAm2/a3zvvhOnTKTs+j1tvtXldJ50E77yT2GAcCmVlNolr1iz41792/GgREZHdZVYL+dtvbQnCE05IXF7kkSPhH/+wQeGRIwHrub7iIpu81acP/POfwayRkNS1iBNh6FDLVX333ZbnU0REosqsgPzAA7B2rV2ClIhxypIS64fu0cOWRcQuq+3e3ZJy3Xuvze0K4tCVl28sT/wBhCMof/ednYd+/ezyMxERqVatArJzLht4CvgZ4IE/eu/nBlGwPbZmjbVYL7wQjj8+MceoXx9ef93u16tHYaHlHFm7Fv7zHxsqDUpoE3+Ua9rUxo333luTuEREYlDbMeS/Aa97738K/AJYWvsixWnWLBuzvOuu4Pe9eDH89rewfr3NpG7cmE2b7LKm776z4eQggzGENPEHWA7Qyy+HH3+0tJj1M6sTRkQkUeL+tnTONQNOAS4H8N5vA7YFU6w49OxpyaH33TfY/a5cCWedZcH+u++gZUtKS6FXL8uUOX065Ma0jseeCWXij6Ii+5Wy116waZMycYmI7IHatJDbA+uAZ5xzC5xzTznn9q7pRQnx8cf2N+hgPHWqRdv1621BhIMPBmw1psmT4eBzP+GySa9x0j0zAl/UYWDXjjRpkLXLY2md+OObb+yHy/bt9p+TqEvOREQyVG0Ccn3gGOAx731n4Htgt7RYzrmrnHPznXPz161bV4vDRTFnjqVhHDMm2P2+/TaceabN2iottVWcsGyc990H2ceuoLTj8oSttBSqxB+bN0O3brBiBbz6Kvz0p6kukYhI6NRmgG8lsNJ7X77o71iqCMje+yeAJwByc3N9LY63O+9tWnNODpx3XqC7ZuJE66YGm11dUMB7WXn07QvN2m+g2Wm7XueciAlXoUj8ATahbulSy4py8smpLo2ISCjF3UL23n8FfOmcK+9D/Q3wcSClitW4cZaB4447ard0UlXOOw+aNNmxlvG6Tvmcfz4ccAA0P2ceLmv33xZpP+EqaB9+aD+KOna08eMLL0x1iUREQqu2s6wHAM875z4EjgYSMMU5iu3bbeGITp3gssuC229ZGTz+OHTubDO2Roxg66TpnPvXPDZtsktrD8qpumMhrSdcBcl7W0Oyc2drFYPNqBYRkbjV6poU7/1CIAFzjGOweDF8/bUN6sZxaU3ULFiPPQbXXmsBplcv/Al5XHkZvPeeNciPOgoGloRwpaWg/PCDJesePdoSpPzud6kukYhIRgjvRaKdO9vqQc2a7fFLo2XB2vuLIrrccovNFu7ZE7DkX889Z73i5dcal4/rhjqtZTy++MLSki1caNd7Dx6spB8iIgEJb0CGuLtJq8qCtXXrNlpffzM0agRPPQXO8frrlvXxwgttTeOKQjPhKkiLF8N//2trTZ9zTqpLIyKSUcIdkONU1eSrvvNe4ecrlsDzz0NODhs22EJFnTrBM8/U4YbgmjUwYwZcfDGcfbYF5KAXdxYRkboZkKvKgjWzfS45bhuX9+oFWMt47VpbS2Lv1KQ7Sa0VK+yC61GjbBLXGWfAfvspGIuIJEjmrYccg4pZsJy3a41X5rQn+4F7wTlmzYInn4Qbb4RjjkllSVNg9WpbQ7JDB/tPuOwyu8Z4v/1SXTIRkYxWJ1vIFSdl9ZjwFEd8X8yWJ0bRrXNbfvzRVgxs1w6GD09tOZNq27adiziPGwf9+8PNN8NBB6W2XCIidUSdDMgQmZT1wesw90VLkXlcO8AmD3/6KbzxRoZ3VX/zjaUHffttSz/aoIGNFefkwKpVwSdaERGRatXJLmsAZs+GK6+08dGCApg7lyVL4J57oHdvGzLNGGVlsGzZzn9fcw20amVLSj7wgCVZOfVUSxEKCsYiIilQZ1vI3HuvBWOAbdsom1nAlRPzaNYMHnwwtUWLWVkZrFtn475r1sApp0DTpta8f+op+Ooru61ebQk9ioth//1tmcqDD4YTT7TVrJrUkQxjIiJprG4G5JISeP99qFfPrmdq2JCX1+czdy78+99pNH+prAyWLIFFi6wb+Q9/sMHtiRNtjHfNmp2tWoD58+HYYy1IL14MbdrAL39pf486auf6xN27p6Q6IiISXd0MyPXrW5CbMwc+/ZS1R+ZzxSV5nH46XHJJqguHdS8PGgRvvmljveWOOMICcps2cNpp0Lat3XJybNWL8mUPe/e2m4iIhEbdC8gbNliGr/33hwsuwHu46nxraD7+eAoSgHz2mV3sXFBgY7p9+li388KFtuJUfj4cf7zNdi6fZZaba015ERHJGHUrIHsPv/+9jZlOnAjYFT6vvmpDyocdlsSyLF4MI0bASy9ZuTp0sNneYK3doqIkFkZERFKtbgXkV1+1S3v+8Q8ANm6EAQPg6KPhf/4nyWW58koLyoMH23jwgQcmuQAiIpJO6k5A3roVbrrJklP36wdYLCwutjWO41jBcc8sWGDN8EcesVljTz8NrVtDy5YJPrCIiIRBxgbkyusdP7FmGp2KimDKFKhfn0WL4Ikn4LrrbEg2YT74wFJ+jR9vY9cLF9plR0cckcCDiohI2GRkYpDy9Y5XbdyCB1Zv+J56Y8aw5tQzLBgCQ4ZYfBw2LEGF8N5axLm5Nlt6+HBbvzlyfBERkYoysoVceb1j7+rRvfdIOjQu4zVg5kyYPNkWM0poj/GCBdCjB/zv/8a9drOIiNQNGRmQK653nLN5LeubNOPHBo35uMQaroMG2Ryqa69NwMFXrIDSUmjf3i5NatiwDi+mLCIiscrILuuc7EgqSO/5+/j7eXH0EPCenOwmjB0L8+bBHXckIGPkrFnWRX3JJRb5GzVSMBYRkZhkZEAuX+/4urdGk7tqKXMOPpomDetz46878uc/20TrSy8N8IDew6OPwumn26INzzyjQCwiInskI7usu3duS6t5b3PSWy/ggb4fjOeoq3vz8fttKSyECRMgKyugg23daqsnjRplmbX+7/+gWbOAdi4iInVFRraQAX41bhQOcEDjshJO+HwJw4fDySfDOecEeKDSUruU6S9/gVdeUTAWEZG4ZGQLGYBvv91lNafnV+VTXAwvvxxQb/KWLdbM3msveOstGy8WERGJU+YG5DlzYPp0mDeP9T/P5/qeeVxwAeTlBbDv0lK4+GLYvNnWHlYwFhGRWsq8gPz11xYwW7e2SVann87t11mD9q67AjrGwIHW1P7b3wIcjBYRkbos88aQb7vN1gXevBmA5cttWcU+faBjxwD2/8gj8NBDcP31lndTREQkAJkVkIuKLEF1z547JlfdeqstHHHbbQHsf8IEuOEG6NYNHngggB2KiIiYzArIt90GDRpYFMbWdRg92pZWzMkJYP+HHgrnnw/PP6+uahERCVTmBOSPPrJAOWDAjug7aBDsu68N+dbKpk2W/ONnP4OxY2HvvWtfXhERkQoyJyBPmwbZ2RaFgRkz7KGhQ2u5rsPGjTY1e8iQYMopIiJShcwJyDfeCIWF0LIl3sPtt1tD+U9/qsU+t22DCy6w/XbtGlRJRUREdhP+gOy9BUzYsZZiQQHMng2DB0PjxrXY9/XX21qNo0bBaafVuqgiIiLRhD8gT5kCP/mJJeiIGD4cDjgArryyFvt9/XW7Xurmm231JhERkQQKd2KQsjIb2z3kkB0t2Fmz7Pbww7VsHf/4I/zqVzBiRDBlFRERqUa4A/LYsbBgATz7LDRsCNg6x23awFVX1XLf3bvb9cZaRlFERJIgvF3Ws2fbsoeHHgoXXQRY+uoZM+CWW6BJkzj3+9prlhKzrEzBWEREkiacAXnuXOjSxfJWr1oF770H2Nhx69bQr1+c+12/Hvr2tUlcJSXBlVdERKQG4eyyLijYGTBLS6GggLd9HtOmwciRtiJiXAYMsCA/efKOLnAREZFkCGcLOT/fAmZWlv3Nz2f4cNhvP7j66jj3OW4cvPCCpd08+uggSysiIlKjcLaQ8/JsreOCAsjP5x2Xx5QpcO+9cWa13LIF+veHzp2VkUtERFIinAEZLCjn5QEw/Cxo1cpialyaNIEXX7SdNGgQXBlFRERiFN6AHPHee5bD4+67oWnTOHawebMt1ZifH3TRREREYhbOMeQKhg+3jJnXXBPHi4uL4fDD4bHHAi+XiIjIngh1QJ43DyZNgptugn32iWMH/fvb0opqHYuISIqFMiC/smAVJ90zg1N7FVO/yXbanbx6z3cyZYrNrL79djjiiMDLKCIisidCF5BfWbCKIeM+ouiTBmxZ3pq9c4sYMeVDXlmwKvadlJRYs7p9e1u2UUREJMVCF5Dvf+NTtmwv5YelOdRrtJ1mx37Olu2l3P/Gp7Hv5IMPYNkyuO8+aNQocYUVERGJUehmWa/euAWA7PxPaNp5BfUalezyeEyOOw6WL4ecnEQUUUREZI+FroWck22rRjgHDbK37PZ4jT77DLyHtm21eISIiKSN0AXkgV070qRB1i6PNWmQxcCuHWt+cVERHHUUPPhggkonIiISn9B1WXfv3BawseTVG7eQk92EgV077ni8WoMGQf360KtXgkspIiKyZ0IXkMGCckwBuKLZs2HsWMskorFjERFJM6Hrso5LWZld3tS2Ldx8c6pLIyIisptQtpD3WFERfPllLRdLFhERSZy6EZA7dIDCwjjXZhQREUm8zO+ynj/fMnPtsw/Uy/zqiohIOGV2hFq1Ck49FW65JdUlERERqVZmB+Q//9laxwMGpLokIiIi1crcgLxgATz7rM2uPvTQVJdGRESkWpkbkIcNgxYtYMiQVJdERESkRrUOyM65LOfcAufcxCAKFIhvv4UVK+ya4+bNU10aERGRGgVx2dP1wFKgWQD7CsY++8DChTZ+LCIiEgK1aiE75w4EzgGeCqY4ASgqgo0b7RKnhg1TXRoREZGY1LbL+mHgFqAsgLIEo29fOPFEW2JRREQkJOIOyM65c4G13vv3a9juKufcfOfc/HXr1sV7uNjMnGm3fv201rGIiISK83G2JJ1zdwOXACVAY2wMeZz3vne01+Tm5vr58+fHdbwaeQ+nnGJd1suXQ+PGiTmOiIhIjJxz73vvc2PZNu4Wsvd+iPf+QO99O6AnMKO6YJxwU6fCnDkwdKiCsYiIhE7mXIc8bRocfDD06ZPqkoiIiOyxQAKy977Ae39uEPuK2333WXauRo1SWgwREZF4hL+F7D188YXdb9kytWURERGJU/gD8ssvw2GHwdy5qS6JiIhI3MIdkMvK4LbboH17+OUvU10aERGRuAWROjN1XnoJFi+GF16A+uGuioiI1G3hbSGXlsLtt0OnTtCjR6pLIyIiUivhbVbOnw+FhTB6NGRlpbo0IiIitRLegHz88RaQDzoo1SURERGptfAGZIBDDkl1CURERAIR3jFkERGRDKKALCIikgYUkEVERNKAArKIiEgaUEAWERFJAwrIIiIiaUABWUREJA0oIIuIiKQBBWQREZE0oIAsIiKSBhSQRURE0oACsoiISBpQQBYREUkDznufvIM5tw5YEeAuWwFfB7i/VFJd0k+m1ANUl3SVKXXJlHpA8HU5xHu/XywbJjUgB805N997n5vqcgRBdUk/mVIPUF3SVabUJVPqAamti7qsRURE0oACsoiISBoIe0B+ItUFCJDqkn4ypR6guqSrTKlLtmklLAAABU5JREFUptQDUliXUI8hi4iIZIqwt5BFREQyQigCsnPuTOfcp865Qufc4Cqeb+ScGxN5/l3nXLvkl7JmzrmDnHMznXNLnXNLnHPXV7FNvnNuk3NuYeQ2LBVljYVz7nPn3EeRcs6v4nnnnPt75Lx86Jw7JhXlrI5zrmOF/+uFzrnNzrkbKm2TtufEOfe0c26tc25xhcdaOuemOueWRf62iPLayyLbLHPOXZa8UlctSl3ud859Enn/vOycy47y2mrfi8kWpS63O+dWVXgfnR3ltdV+3yVTlHqMqVCHz51zC6O8Nt3OSZXfv2n1efHep/UNyAKWA+2BhsAi4MhK2/QHHo/c7wmMSXW5o9TlAOCYyP19gM+qqEs+MDHVZY2xPp8Drap5/mxgMuCAE4B3U13mGuqTBXyFXTcYinMCnAIcAyyu8Nh9wODI/cHAvVW8riVQFPnbInK/RRrW5QygfuT+vVXVJfJcte/FNKnL7cDNNbyuxu+7VNej0vMPAMNCck6q/P5Np89LGFrIxwGF3vsi7/024EWgW6VtugH/jtwfC/zGOeeSWMaYeO/XeO8/iNz/FlgKtE1tqRKqG/CsN+8A2c65A1JdqGr8BljuvQ8yeU1Cee/fBNZXerji5+HfQPcqXtoVmOq9X++93wBMBc5MWEFjUFVdvPdTvPclkX++AxyY9ILFIcp5iUUs33dJU109It+xPYDRSS1UnKr5/k2bz0sYAnJb4MsK/17J7kFsxzaRD+8mYN+klC5OkW71zsC7VTyd55xb5Jyb7JzrlNSC7RkPTHHOve+cu6qK52M5d+mkJ9G/XMJyTgBae+/XgH0JAftXsU3Yzg3AH7Eel6rU9F5MF9dGut+fjtI1GqbzcjJQ7L1fFuX5tD0nlb5/0+bzEoaAXFVLt/LU8Fi2SRvOuabAf4AbvPebKz39AdZl+gvgEeCVZJdvD5zkvT8GOAu4xjl3SqXnQ3NenHMNgd8CL1XxdJjOSaxCc24AnHNDgRLg+Sib1PReTAePAYcBRwNrsO7eysJ0XnpRfes4Lc9JDd+/UV9WxWOBn5cwBOSVwEEV/n0gsDraNs65+kBz4usuSjjnXAPszfC8935c5ee995u9999F7k8CGjjnWiW5mDHx3q+O/F0LvIx1t1UUy7lLF2cBH3jviys/EaZzElFcPjQQ+bu2im1Cc24iE2jOBS72kQG9ymJ4L6ac977Ye1/qvS8DnqTqMobivES+Zy8AxkTbJh3PSZTv37T5vIQhIM8DDnfOHRppxfQExlfaZjxQPuvt98CMaB/cVIqMuYwClnrvH4yyTZvy8W/n3HHYOfomeaWMjXNub+fcPuX3sck3iyttNh641JkTgE3lXUNpKOqv/bCckwoqfh4uA16tYps3gDOccy0iXadnRB5LK865M4FBwG+99z9E2SaW92LKVZo/cT5VlzGW77t0cDrwifd+ZVVPpuM5qeb7N30+L6me+RbLDZut+xk2+3Bo5LE7sA8pQGOsq7EQeA9on+oyR6nHr7Bujg+BhZHb2cDVwNWRba4FlmCzK98BTkx1uaPUpX2kjIsi5S0/LxXr4oBHI+ftIyA31eWOUpe9sADbvMJjoTgn2I+INcB27Fd8H2z+xHRgWeRvy8i2ucBTFV77x8hnphC4Ik3rUoiN3ZV/XsqvpsgBJlX3XkzDujwX+Rx8iAWBAyrXJfLv3b7v0qkekcf/Vf75qLBtup+TaN+/afN5UaYuERGRNBCGLmsREZGMp4AsIiKSBhSQRURE0oACsoiISBpQQBYREUkDCsgiIiJpQAFZREQkDSggi4iIpIH/B55jbJLrQVQhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(x, y, 'o', label=\"data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "ax.plot(x, res.fittedvalues, 'r--.', label=\"OLS\")\n",
    "ax.plot(x, iv_u, 'r--')\n",
    "ax.plot(x, iv_l, 'r--')\n",
    "ax.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time-Series Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.arima_process import arma_generate_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gerando dados\n",
    "np.random.seed(12345)\n",
    "arparams = np.array([.75, -.25])\n",
    "maparams = np.array([.65, .35])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parâmetros\n",
    "arparams = np.r_[1, -arparams]\n",
    "maparam = np.r_[1, maparams]\n",
    "nobs = 250\n",
    "y = arma_generate_sample(arparams, maparams, nobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dmpm/anaconda3/lib/python3.6/site-packages/statsmodels/tsa/base/tsa_model.py:171: ValueWarning: No frequency information was provided, so inferred frequency M will be used.\n",
      "  % freq, ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "dates = sm.tsa.datetools.dates_from_range('1980m1', length=nobs)\n",
    "y = pd.Series(y, index=dates)\n",
    "arma_mod = sm.tsa.ARMA(y, order=(2,2))\n",
    "arma_res = arma_mod.fit(trend='nc', disp=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                              ARMA Model Results                              \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                  250\n",
      "Model:                     ARMA(2, 2)   Log Likelihood                -245.887\n",
      "Method:                       css-mle   S.D. of innovations              0.645\n",
      "Date:                Thu, 02 Aug 2018   AIC                            501.773\n",
      "Time:                        12:16:44   BIC                            519.381\n",
      "Sample:                    01-31-1980   HQIC                           508.860\n",
      "                         - 10-31-2000                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "ar.L1.y        0.8411      0.403      2.089      0.038       0.052       1.630\n",
      "ar.L2.y       -0.2693      0.247     -1.092      0.276      -0.753       0.214\n",
      "ma.L1.y        0.5352      0.412      1.299      0.195      -0.273       1.343\n",
      "ma.L2.y        0.0157      0.306      0.051      0.959      -0.585       0.616\n",
      "                                    Roots                                    \n",
      "=============================================================================\n",
      "                  Real          Imaginary           Modulus         Frequency\n",
      "-----------------------------------------------------------------------------\n",
      "AR.1            1.5618           -1.1289j            1.9271           -0.0996\n",
      "AR.2            1.5618           +1.1289j            1.9271            0.0996\n",
      "MA.1           -1.9835           +0.0000j            1.9835            0.5000\n",
      "MA.2          -32.1800           +0.0000j           32.1800            0.5000\n",
      "-----------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "print(arma_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conheça a Formação Cientista de Dados, um programa completo, 100% online e 100% em português, com 340 horas, mais de 1.200 aulas em vídeos e 26 projetos, que vão ajudá-lo a se tornar um dos profissionais mais cobiçados do mercado de análise de dados. Clique no link abaixo, faça sua inscrição, comece hoje mesmo e aumente sua empregabilidade:\n",
    "\n",
    "https://www.datascienceacademy.com.br/pages/formacao-cientista-de-dados"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Obrigado - Data Science Academy - <a href=\"http://facebook.com/dsacademybr\">facebook.com/dsacademybr</a>"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}