Cap08-07-StatsModels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font color='blue'>Python Fundamentos - Capítulo 8</font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Versão da Linguagem Python\n",
    "from platform import python_version\n",
    "print('Versão da Linguagem Python Usada Neste Jupyter Notebook:', python_version())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statsmodels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Regression Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Para visualização de gráficos\n",
    "from pylab import *\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels as st\n",
    "import sys\n",
    "import warnings\n",
    "if not sys.warnoptions:\n",
    "    warnings.simplefilter(\"ignore\")\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "warnings.filterwarnings(\"ignore\", category=FutureWarning)\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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.18.2'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.0.3'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.11.1'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "st.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = sm.add_constant(X)\n",
    "y = np.dot(X, beta) + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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:                Mon, 23 Mar 2020   Prob (F-statistic):          2.83e-239\n",
      "Time:                        20:09:10   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": 10,
   "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": 11,
   "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:                Mon, 23 Mar 2020   Prob (F-statistic):           6.30e-27\n",
      "Time:                        20:09:10   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": 12,
   "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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc4d1028fd0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFlCAYAAAAzqTv+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzde3zP9RfA8ddn313dr4m5q1xyDWUki1zKNbqhFF1+SklKKEVUiC7kEuWSqFDMrUI0ZHO/5K6Q2pDrxmyz2/v3xzHXmV2+1+08H489tu/t83l/N77n876dYxljUEoppZRzebm6AUoppVRupAFYKaWUcgENwEoppZQLaABWSimlXEADsFJKKeUCGoCVUkopF/B25smKFStmypcv78xTKqWUUi6zefPmk8aY4mk95tQAXL58eTZt2uTMUyqllFIuY1nW4Rs9pkPQSimllAtoAFZKKaVcQAOwUkop5QJOnQNOS2JiIhEREcTHx7u6KQ7j7+9P6dKl8fHxcXVTlFJKuQmXB+CIiAjy589P+fLlsSzL1c2xO2MMp06dIiIiggoVKri6OUoppdyEy4eg4+PjKVq0aI4MvgCWZVG0aNEc3cNXSimVeS4PwECODb6pcvr7U0oplXluEYDdyZAhQxg9evQNHw8JCWH37t1ObJFSSqmcyOMCcMjWSBqNWEmFAUtoNGIlIVsjnXt+DcBKKaXswKMCcMjWSAbO20FkVBwGiIyKY+C8HdkOwh988AF33HEH9957L/v27QPgyy+/pH79+tSqVYtOnToRGxtLWFgYCxcupF+/ftSuXZsDBw6k+TyllFLqZjwqAI9auo+4xOSr7otLTGbU0n1ZPubmzZv5/vvv2bZtGz/99BMbN24EoGPHjmzcuJHt27dTtWpVpkyZQsOGDWnXrh2jRo1i27ZtVKpUKc3nKaWUUjfj8m1ImXEkKi5T92fEmjVrePjhh8mTJw8A7dq1A2Dnzp0MGjSIqKgoYmJiaNmyZZqvz+jzlFJKubnoaChY0Gmn86gecKlCAZm6PzueeeYZxo0bx44dOxg8ePANtxFl9HlKKaXc3KBBEOm8dUUeFYD7taxMgI/tqvsCfGz0a1k5y8e87777CAkJIS4ujnPnzrFo0SIAzp07R8mSJUlMTGTWrFmXnp8/f37OnTt36faNnqeUUsrNbdoEnTrB+vVye8AAcGLGQo8KwB3qBDK8Yw0CCwVgAYGFAhjesQYd6gRm+Zh33XUXjz/+OLVq1eLBBx+kfv36AAwbNox77rmHRo0aUaVKlUvPf+KJJxg1ahR16tThwIEDN3yeUkopN2QMhIZCixZQvz6sWAEHD8pjgYFwyy1Oa4pljHHayerVq2eurQe8Z88eqlat6rQ2uEpueZ9KKeXW2rSBJUugRAno2xd69oQCBRx2OsuyNhtj6qX12E0XYVmWNRVoAxw3xlS/eN8ooC2QABwAuhtjouzXZKWUUspOLlwAPz/5uVUreOgh6N4dAuy/figzMjIEPR1odc19y4HqxpiawH5goJ3bpZRSSmXfli1QsybMnCm3X34ZXnrJ5cEXMhCAjTGrgdPX3LfMGJN08eY6oLQD2qaUUkplTUoKjB4NDRrA+fNQ2v3ClD32AfcAZtvhOEoppVT2HTkCTz8Nv/4KHTvC5MlQtKirW3WdbK2CtizrbSAJuOH+G8uyXrAsa5NlWZtOnDiRndMppZRSN7d+PYSFSeD94Qe3DL6QjR6wZVnPIIuzmpl0llIbYyYDk0FWQWf1fEopZU8hWyMZtXQfR6LiKFUogH4tK2d8S6MxsHMnLFgAf/8NTz0FTZpIJqUVK+QDv0gR+V606OUFQMpx4uNh7Vpo1gwefhgOHIBbb3V1q9KVpQBsWVYr4E2giTHGo6sPnDp1imbNmgFw7NgxbDYbxYsXB2DDhg34+vq6snlKKQdILeySmls+tbALkH4QPncOBg+GkBA4dEjuK1FCgi/A3r2S2OFKvr4waxY88oi934ZKdfSoBN1t2+TvUrKk2wdfyNg2pO+AYKCYZVkRwGBk1bMfsPxisfl1xpieDmynwxQtWpRt27YBUgs4X758vPHGG5ceT0pKwtvbo1JmK6VuIr3CLlcF4PPnYelSiI2FJ5+EPHngxx+henXJmtS2rXzYp6pRA7ZuhdOn4dQp+b5zpywEAli2TLIvPfMMlCrl+DeaG2zeDO3bw5kz8O23V/893NxNI4sxpnMad+fokj/PPPMM/v7+bN26lUaNGlGgQIGrAnP16tVZvHgx5cuXZ+bMmYwdO5aEhATuueceJkyYgM1mu8kZlFKulKHCLnPnwosvSiCtXVsCsM0mQ5s3uijPk0eeeyMrV8LIkfDuu9C6NTz3HDz44I2Pp9I3Z45czBQrJsPP6f3u3ZBb/dX79JERBHuqXRs++yzzr4uIiCAsLAybzcaQIUPSfM6ePXuYPXs2a9euxcfHh5deeolZs2bRrVu37DVaKeVQpQoFEJlGEC5VKEB6rS+/DN99J6kK586Fxo0vPyk7wXLECHj2WZgyBaZPh4ULITgYfvst68fMzf74A+rUgXnzZCrAw7hVAHYnjz766E17sitWrGDz5s2X8kfHxcVxixPziCqlsqZfy8pXzQHDFYVd9u+H+fNh2DAZZrZ37/T22yUQDxsGixfLgi6AhATpgdx9t33Pl9OcPy/zvNWrw9ChkJjosYvc3CoAZ6Wn6ih58+a99LO3tzcpKSmXbqeWHDTG8PTTTzN8+HCnt08plXWp87ypq6BvCzB86HeI+nVaAYGystnRPSofH1k4lGrSJOjdW/avjhzpkT06hzt8WOZ7jx+Hv/6SIX8PDb7gYdWQXKV8+fJs2bIFgC1btnDo4urHZs2a8cMPP3D8+HEATp8+zeHDh13WTqVUxnWoE8jaAU051Covy6e/Qv1Br1y9stnZuneXHve338Idd8CYMZCUdPPX5Ra//y6jA4cOwdSpEnw9nAbgDOjUqROnT5/mzjvvZNy4cdxxxx0AVKtWjffff58WLVpQs2ZNmjdvztGjR13cWqVUhiQkwGuvwf33S290zRqoUMF17cmXD4YPhx07ZNV0nz7SG87tjIEJE+TvVLCgJNlodW15As+k5QidJLe8T6U8gjGyL3fePFlwNWIEXDHt5HLGyF7jUqXgnnvg7FmZ63TTjE4OZYwM1SclSUGFQoVc3aJMyVY5QqWUynEsSz7UmzWTyjjuJrV9qfr3l4VhEyZIbuPc4N9/JfiWLSvD8v7+4JWzBm1z1rtRSqn0xMbKXCLIvl53DL5p6dlTesOdOsFjj8kipJxs1SqoW/fyEHyePDku+IIGYKVUbhETI8kvWrSAY8dc3ZrMqVVL5j4/+EDyT1erJjmncxpjZDtMs2Yy3D5xot0OHbI1kkYjVlJhwBIajVhJyNZIux07qzQAK6VyvrNnZeHO6tXw1VcekSf4Oj4+8NZbkuqydm2oWNHVLbKv2FgpavHaa5Lic/16qFLFLodOzf0dGRWH4XLub1cHYQ3ASqmcLSoKWraUD/Tvv4cuXVzdouypVk3q3FaoID3GLl1g2jQpQO/JUlIks9WwYZJvu0ABux06vdzfrqQBWCmVs02ZIgn7586FRx91dWvs6+xZiIiAHj2gUSPYuNHVLcqc2Fh47z2pMpUvn7R/0CC7z/dmKPe3C2gARvI+t2/fnttvv51KlSrx6quvkpCQQGhoKG3atLnu+YsXL6ZOnTrUqlWLatWqMWnSJBe0WimVIX37ygd7hw6ubon9FSwIoaGSV/rvvyVRRffuks/a3a1cKdWjhgyBn36S+xyU1apUoYBM3e8suT4AG2Po2LEjHTp04M8//2T//v3ExMTw9ttvp/n8xMREXnjhBRYtWsT27dvZunUrwcHBzm20Uip9sbHQuTPs2ydbemrVcnWLHMfLS1YL79sHb74pgc3Hx9WturGoKKkC1ayZtD00FB5/3KGn7NeyMgE+V+f2v5T724U8MwCHh0vGmPDwbB9q5cqV+Pv70717dwBsNhuffvopU6dOJTY29rrnnzt3jqSkJIpe3BDv5+dH5cqu/SMqpa6QnAxdu8Ls2bB3r6tb4zwFCkgO6X37IH9+SdzRvv3l3qW7ePFF6bH37y9zvk2aOPyUHeoEMrxjDQILBWABgYUCGN6xxtW1n13A/RJxpNWbfOwx2a8XGyvzHH/8IRP2Xl5Qsya8+qrUhDx5UrLbXCk0NN3T7dq1i7p16151X4ECBShbtix//fXXdc8vUqQI7dq1o1y5cjRr1ow2bdrQuXNnvHLgHjWlPI4xsoo2JERyKbdv7+oWOZ+/v3yPiJALkNat4aGHJOA1biwjAs62caNcIFSuLFup3nhD9vk6UYc6gS4PuNfyvKgRHX15tV9Kitx2sq+++ooVK1Zw9913M3r0aHr06OH0Niil0vDZZ/D55xKEe/d2dWtcq0IFySs9ejSEhUlPs2pVOHLEOeePi5PV2fXry9z0Rx/J/RUrOj34ui1jjNO+6tata661e/fu6+5LV1iYMQEBxths8j0sLHOvv8by5ctN48aNr7ovOjraFClSxCxZssS0bt063defOHHC5MuX76bnyfT7VEplTmKiMffea0ynTsYkJ7u6Ne7l/Hljpk83pnNnY1JS5L5vvjFm1arLt+1p6FBjChc2BoypVs2YceOMiY62/3k8ALDJ3CAmel4POChIMsAMGybfg4KydbhmzZoRGxvLjBkzAEhOTub111/nmWeeIU8a5a5iYmIIvWJYe9u2bZQrVy5bbVBK2YG3NyxbBt98kyPTFmZLnjyyUOvbb2UIOiUF3n33cq/4k09g507ptWaWMTL99/PPl0cnz5+HBx6A336T4/bqZdd9vTmFVkMC/v33X1566SX27t1LSkoKDz30EKNHjyY8PJwHH3zw0oIrgO+++47hw4dz4MABAgICyJs3L2PGjKFevTSLXVziDu9TqRzpr7/g7bdh8mTZlqMyJjZW9kZPmnR5Qeu778q+3Kgoybp1xx1w++3y3d9fKhHlzy/7qkeOhAMH5Ct1KvDXX2V1szGumWt2Q1oN6SbKlCnDokWLrrs/ODiYuDSuCBs3buyMZimlbubkSXjwQThzBk6c0ACcGam94qeflsVa27ZJli2QBVzffSeB+EqzZ8ui2Ph4eX6lSjIKWamSvDZ1RbMHBt+QrZGMWrqPI1FxlCoUQL+WlR2+aEsDsFLKM8XFQbt2UrZu5Uq47TZXt8hzValydd7l6tUlmcepU7B/v3xduCALqkB2o+zf75q2OkBqrujUdJWpuaIBhwZhDcBKKc+TkgLdusG6dTKM2rChq1uU81gWFCsmXzn895termgNwEopdaUjRyT4jholNXJdwBVDlsoxXJUr2i0CsDEGywPnDDLKmQvdlMoVSpeWPa4umvN11ZClcoxShQKIjIrjrsg9NPhnB+vK1mBLYFWH54p2+Vp9f39/Tp06lWODlDGGU6dO4Z+anUYplXUrV0pxheRkWZHrogt3dy1vp7KmX8vK3Be5k7mz+tN3zUxmff82Qf/td3iuaJf3gEuXLk1ERAQnTpxwdVMcxt/fn9KlS7u6GUp5tr17Zbg5MFAq6LhwX6m7lrdTmXTgAKxfT4cuXbi9UBSWScEGkJzEoHzHuTOnr4L28fGhQoUKrm6GUsqdnTgh+Yx9fWHxYpcndUgdskzrfuXGwsMlgZNlSZ2AX3+FgABo25Y7u7SDKWMgIQFvX1+57WAuD8BKKZWu+Hip5Xv0qHxoli/v6hbRr2Xlq+aAwT3K26l0hIdLsZ+EBLldogQMHQo9ekhykdQsi6Gh8rxsZlnMCA3ASin3tnWrJH345hu45x5Xtwa4vNBKV0F7iMRECazJFy+YvLzglVckg9qVgoKcEnhTaQBWSrm3oCA4eFB6LG7EHcvbqWucPQtvvgl//in1A3x9pQfs6wtNm156WkqKdJBnz5afx41zTvM0ACul3NPMmZJ96dln3S74pjp9WkbGb7kFihQBm83VLVKX/PQT/O9/sme8Tx+oV++qIWbTIIiNGyTozp0rCdX8/KSkvLNSWWsAVkq5n2XLZG7u3nvhmWfcJrIlJ0tt+V9+gaVLwWtdGE34jd9oygavIIoXl2B8yy1yzVA/KZx6MaFU/l8wxds5b2gzV4uKkuHlmTPhzjvhhx8uTV2YBkFs8w9i9myY0xUOHQIfH2jVCoYPl8ym+fM7r6kagJVS7iUsDB5+WJL7//ijy4Pv0aMSbH/5BZYvh5jTF7ifUD4sPIlgQsCySPL+gKX3j6Th2lGcj83HuYP5SE5MplrCNgyQ8JMfA5usoOnbQTRrptUSHSI8XHq39etLlrTBg6Wik68vIAue+/aV/C02m1RLfOcdWd9XuLBrmqwBWCnlPrZvl+1GgYES9Vz1yXixKRO7hVP4j1BCCeZs8dv4qcDL3BX3Mz5x5yDGBzBgDD4pCbS5Yz+UakaRmBiIieH89p14HU3BAgwJ1A4bx6EW03gqsCt1ejfmmR5eFCvmsreXs4SGQosWMoHr6yu1iS9WZoqIgNdfhzlzoGJFqb7YsSNu8bvXAKyUch+rVske319/ddm8b3y8rNdZM+J3fk1pijdJ4OcPPyzF68Wd8OQTMlaZL59cLKQu6unS5dIK2pCtkcweO4epMwfgk5xEks1GQo1zPL0jhBciv+Sf/mWY/lZnTrfqyhNtz1PjVCjW/cHXrcDNLfmms/U+Dx+GJ5+Ulc4gf4+wMBKCmvDZZ7LTKDlZvvfrJ2WN3YXlzBSQ9erVM5s2bXLa+ZRSHuLKVS/R0S7L8bx6NTz/PFj79/J73pYUO/+PPGCzSVQeOPDqF6QOe16zb7TRiJVp5haulAdW3BbNuS9mkef3pfxNOUqao/iRgOXvi9fKFVcF8bT2Gg/vWMOuQdjVQT5b73P5cujcWUpTJiVJpPX1ZePIFXSbEMTevXKt9Nln4Kp8T5ZlbTbG1EvrMZ2JUEq51vHjUu4uPFxuuyD4RkdDz57QrEkiz50Yzi6f2hTzOiM9W5tNvgcHX//CoCAJytf0XFNTUm4JrMqEoMfYElgVgIOxQJcu5F+9BNt/RynzwkP4WwnYSIb4OLb3/ZqkJDmGM/JNpwa/yKg4DJeLSoRsjbTbOW4my+9zyhRZPVWypOwTDw0l+o1hDApawd29g0hIgEWLYMEC1wXfm9EhaKWU60RFQcuWsG/f5SQJGWDPXltICPTqBceOwWdPbOSV79+CRx7h5/+9zfwf13Db7k38Va0eD/mXpUMGj5mhVJXFi+P7dBeYMQVzIQErJZla6yax9JYoSs0Yabd80+n9rlxVB/dKWX6fjRrJCvmxYzF58vL12tt5eVwQycnw3nuy/dedhpvTogFYKeUa589DmzawaxcsXChbjjLAXqUAz56VLcYLf7hAjwqh9FjXkvr1G8KbWwjhFjlHwYosC6oIwJpMnCPDqSovpj+0QkMxd9/D7i9Wcd8Po7DahtCn9HN81ulBzDVBJDP5pm/2u3KHohIZzqudminjzBmYPh2qVIEpU4iJgZeelkRpwcHSMa5Y0SlNzzYdglZKOV9CglQ2Cg+Hb7+VocQMssfQ7MmTkgipyLwvOZW3LBMOP0T9YofkwTp1sn2ODnUCGd6xBoGFArCAwEIBN57TvDiMbTVrSrW575GwYx/bb3+E2yMOEjG9KXduOMZL4XO4K3JPpvNN3+x93CiYO7OoRL+WlQnwuXqr2XXvMzWP85gxMGOGjC0jK9Xr1oVZs6TX++uvnhN8QXvASilXsNkgTx748ktJPZQJ2e21RUbKjpVu+97izZThWOeROd5jxy5NFtqjZ5jVVJUFq5fhnv0zCVuVSOsu6/jut5fxJ54Eb1/WfTmXJpk45s3ehzsUlbhpXu2kJBlPTi2iYLNhdu7ii8h2vPaaZCBbsSLtKXp3pwFYKeU8Bw9KsCtdWjIUZSEjRXZKAR48KAkYHo/8hP7Jwy8/kJwsq5kvLqZyh3KDDZv48MOLv+P1TgIW4JuUQO2Jn0PXByV9Uwbc7H24S1GJdC9WHn0Ufv8dvL3BGIyvL4N+DebD32TgZMYMKF7cqc21Gx2CVko5x/r10KABdOsmt7OYDipDQ5Zp2LVLppmjo6FnL5uMQQcEpLnKOavnsDfvZsF4BfhhvGykYKPghl/5r/zdpGzZlqHXZ+R9dKgTyNoBTTk0ojVrBzR1v33Gzz8vE7yrVxPxv2E8UmgFI1cHMXIkLFniucEXdB+wUsoZfvxRkiWUKiVJ8itnL5BldhX0pk3wYMsUqtj+ZOLKylSvjuw9XrfuhvVfXb0/9pKLe43P3hXMpPeO8VT4i2ws3ZH6Gydw6603f7nbvI+MMgYmTpQh5z59Lt312WfQv7/sOvr+e6dWDcyW9PYBawBWSjmOMfDJJ5KC6J57ZLWzk7ssq1dD+9ZJfGmeoyM/4rVrJ5Qr59Q22Isx8M2Y07wx0Afy52f+u1tpFLBF9lI7qYi8Q8XGwosvyrhy+/Ywfz6nTlu0ejiOTWsCyHP7Map33s/ADpXc+yLiCukFYJ0DVspNeFxPJSPi42XLSMeOMowY4Lw5VJDOdueOF/jBtzPNY+ZLPsKyZZ3aBnuyLOjWpwj1WsATT4DtlRcxrAfLwvL3l9VInhqEf/hBqhgdOyZLmgcNYm2YRYdHkjh1wo/CD+wi/11/898FsrTtzB3pHLBSbsAdMhLZVUyMpAcMCJAh3jlznB58V48MZ2Ob91hra0zzc/NlC8s77zin0KuDVasmU+qn67XEYGEZg4mLk6FbT7R8uSy2OnYMfH1JadacER950aQJnEtIpMSTaylQ9+9Lfzp7ZwRzFQ3ASrkBZ6QddJpFi6BmTclSBFC0qNPr7/0xOZx6A5rxjnmPO2M3wttvQ+/eTm2DowUEwENjW5Hi608SXhgs+OYbkuf+6OqmZdzJk/J906ZLF0YmOZmvu4cycKAMnJTothq/W89e99Irt1iFbI2k0YiVVBiwhEYjVnrMhasGYKUcKTYWNm+G3bvl9smTMkT4yCPw2mvw8ccwZw7ehw6m+XJnZiTKtoMHoW1byX4fECD5HV1g3z5Y0CcUXxIkLHl5Qd68LmmLwwUF4R26gri33uftRqv4H1/QfGw7Dh5E/s2lVghyN8nJ8OmnMhe/fLnMX/v7k+JlIz7Zl+l/BzNxoiS+Kl0i7S1XqVupPHn0SOeAlbK3IUMkOfzOnRKUjJHe4LRp0hssXFj2xPzyi6RjBLq06MHwwh255dwp5s3sx4Yyd7K2XG0O1LzHpW8lw0JCpCqNtzeMHi29zQzuVbWn//6D1q2S+TB5N15+PpDEjQsp5BRBQeQPCuJDAzNmNOb73tCgxnn+9m5GQOkiWBMnyt/iBqu9nW73bskBum6dpCKtVo34ooF888QKDk0P5UCZYMYuCqJWLXn6zZKFuEM+66zSVdBK2UNkpBSRB6hRQwqD33knVK8uX3XqXF+SxRgpRhARwS9HLvDa7ycpcPo/3lk5hQb//EGx2Gh5XuXKkjGqcWPnvqeMOH9eepdHjsBbb8EHH1z+PThZTAwENzE898cr9EwaL3tWChZ0j6DjRP/8Az16QMCKRXwV8Aol4g7LXmdjwM/PdQu1wsPh/fdh2TL5u4wdC507s/xXi1694M8/4amnYPx4yJ//6pemt0CxwoAlpBXFLODQiNYOf1s3o6uglXKU3bsl8CxfLp8gpUpJgtqMzHlalvSGCxemVQ2IvyWSUUt9eaV9fwIL+DG0kqFp5A75wEzd8DlzJnz+OTRvLlWEGjRwSU+Tgwehb184fRpWrZL3PX2689txUVISPPYYNN86kp5mPLzxBowY4bL2uFLZshLjJkxoy539mjLfuzX3Jq3CAtlbe0XGL6e4cEGC70MPyap4Ly/4+muO1GlN384yzHz77dLm5s3TPkR6mbLcIWtZlhljnPZVt25do1SO8M8/xnTvboyXlzEFChjz/vvGxMQ4/rzz5hnTsKExNpsxYEz+/Ma0b29MXJzjz334sDEffWRMgwZy7jx5jBkxwpjERMefOx0pKcY895wxT/G1tKtLF2OSk13aJnexb58xz90ZZmLxN8l4mSS/AGPCwoz55BNjvv/esX+7w4eNGTjQmOLFjWnX7tK/2RSbzfze+kOTP78xfn7GvPde9v75zt8SYaoM+tmU67/40leVQT+b+Vsi7PdesgHYZG4QE3UIWqnMOnECypeXbtfLL0tB9mLFnNuGqCj47TdYuhT+/lvmkwFefVV6OffeC1WryvB1VhcgGSO9+TJlZO566lSZu6tbV5andusmOZ2v4ez9zMOGwXvvJvFPifqUql5UNv/6+jrsfO4kI7/rpCRYODCcfZNDWXg2mDzB97Ag4i7y/XXxb9u7t6xa37w568P1F7N1ERwsCw/HjZOkKyDzvK1bQ58+mIQELqT4cr9ZQcGWQYwbB7fdlt3fgnvvoddMWErZw4kTl7M4ffWVjJe5W0albt1g/nyZEE3VubOU/AMZ7ytVSgJnbCycOyfvqVIlGR6cPFnuO3ZMEu0eOiQfpr16SQHdM2fSfc/X1p8FWTBzw1J82TR9OnTvLm97+qdnsGxeMr+YC2T2dx0fD5MmwfDhcPy/FN6ps4TXrU8osCVUnuDlJXPEixbJ37hixctTKVcG2CsDtDHw889yQZaUJBc+lSrB0aOSw7lnTyhXjj/+gCWDwjm3KJSdxYLpNjGITp1yxJbsm9IArFR2LVgAXbtKtp5M1K51iYQE2L8f9u6FPXukl/PMM7IlJW/e67em9OkjW0JiYy/3lgMC5MO2Y0dJCZjB9JGNRqxMcz4usFAAawc0zd77usayZdDroUN8eutIWuz+DN8C/jd/UQ6S1d/1+fMwYQKMHAmnTsGi8q/Q+u/xWBhZrNWli2Qty5dPFhSWKCEXYykpEmDffVfm/Q8flq/Y2MsHt9lke92wYZw458+338pF0rZtslTh5ZclydW1i6xyMl2EpVR2TJsGzz0H9evLl7vz9b28+vpK3t6yUGzPHumh5M0LBQpIjwUk6J48KZ+OWRzCtUcd3YzYuRO+aP8zG+hGoegLWP+9DgVut+s5nCE7Q6dZ/V3nzSupuXv2lIXIY0d0oSlT8CWBFMuX7WU7cXYMUeAAACAASURBVOeY+8jz1x+wfTvJP/2M7eJFW1L8Bf77dQ2BZ47LFEerVtILnjABkpMxvr6sKtKRTx/356efpFNct+6lBc9On6lxdxqAlUrPqFFSDLxFC6noky+fq1uUdZYlQ4s3GkK2LJnrzQZnrEg9fhyGPRDK3Pg2eJGClegnFw63e1YAvnYIOTWBBGQsx3F2f9f580uCsF69gpg7YAXxv4QyKzKYNR8EYVlQqxaUrhqD1XAuc9b0xDsliUSbN68Vb8X9rzxK3ZK3cuaMLIT3tT2Gz9pQJuwOZulbQdx6qwysPP309deB6jIdglbqRlasuFi9/XGpzpJLFvZkh6PngC9cgGZNDRPCa1PT/CF32myyEmvgwGwf35myO1zviN91XJzkx1i9WkaZV61JJiXJRgPCud+2gt+Sm7GOtBdp+fnJbMXTT8v1qvcV3Tt3XiTlaDoErVRWNG0Kc+fCww/Lh7y6qdQPVUd82Boj63qOhB2iSsAhSPK5PC/pgZmusjtc74jfdUAA3H+/fAGUf2MZ8ccKsvefouyK7YqXXyKFA3Zh809kWs/aFC4MRYpw6Xta16jZ7ennZDcNwJZlTQXaAMeNMdUv3lcEmA2UB/4GHjPGnHFcM5Vykvh4eOkleP11yWT1yCOubpHHSS9pQnaMGCFrg957ryK+L+yXZCCrVnlspit7DNc76nedKrCYH5HeZ/AvffXHe2ChANq0ydgxPDlVpKNlpBjDdODaZZ8DgBXGmNuBFRdvK+XZzp2TbD3Tpsk4nHIb8+bBT2+t4fuaH/LOICOZwRo2lGFnDwy+IDmOA3yuHlm5MsexO7BHG521MM8T3TQAG2NWA6evubs98PXFn78GOti5XUo5V0qKLNNcvVrSPT77rKtbpC7asgUGd/2LRd4P82j8DKzzMTd/kQfoUCeQ4R1rEFgoAAvpVTpqv3RW2aONN+rRu12qyNOnZaGlE2VoEZZlWeWBxVcMQUcZYwpd/NkCzqTeTuO1LwAvAJQtW7bu4cOH7dNypezp3XdlIc/48TIErdzCkSPwQN0zLDrZgAoFTuG1fp19Uicpp3F2cpZMOXFCJrC9vSV/+Mcfyx56O66oT28RVrbrAV/MdXnDKG6MmWyMqWeMqVc8g5v5lXKq5GTYsEFKyLz4oqtboy6KjYWObRKYeKITFay/8QqZn+Hg66kF2nMit+vpHz0q+5abNZOpjNWr5f5evWDjRqde4GV1FfR/lmWVNMYctSyrJHDcno1SyqlsNsn0k5SUO3LjeYCUFPigTTjPbp1KY6/VeE2dnuFyjLrq1v04erFYhhw5InukVqyQJfWVK0sls4oV5fEKFa4vGepgWe0BLwSevvjz08AC+zRHKSeKipIkwkePShD283N1ixTy2Ti2czhv/9aMZ61pePn6XM7WlQHprbpVuVBKinwvWlQ2Or/7rqRS27NHpp3Kl3dZ024agC3L+g4IBypblhVhWdazwAiguWVZfwIPXLytlOdISYEnn4TvvoMDB1zdGnWFUaPg1jmf4c8FvEyy5K4ODc3w63XVrQJkROvzz6F2bUmA7ecHa9bAkCGyxdANRrtuOgRtjOl8g4ea2bktSjnPkCEy7DxunJTuUxnmyKxGM2bAb/1/ZhE/YnkBli3TiTacVqD9woXLoyabNsHWrVJxwMdH2uzjI2X4fHzse151cytWSGnOXbukallUlCTBdoOgeyXNhKVyn/nzZeipe3dd8ZxJjpxf/eUX+LJHOMu9OuFVsybWqI9kUUwmE230a1k5zVW3dtlfe/y4XLgtXCjlmHbvltzaixbB0KHXP//sWQnAo0dLSaAOHaSAgSfnFHdn58/LPO+PP8p8bkgItGvndoE3leaCVrlLSopUNPL2lixK/vYrYZcb8t06qtzghg3wUpNdrEhsTL5yRbGF/S5l8LLI7n+LP/6QFfLh4TJJXaaMfLD36ycBOCYGoqOlFGRionwlJEhFAy8vCc5jx0r9P39/6ZV17ixfyn6MgbZtJUlL3752/f+dVVoPWKkrRUXJYoySJe12SLfe63iF7AamCgOWpLnn0AIOjWidpTbt3w+NGsFbSUPp7TsR27owp69GTdPatfKBfu+98N9/Mpzctq1UHKhVK/O9qqQk+P13GYEJCYG77pKfATZvlttu2lNza1FRkjp2yBC5MDLGrX6PDt0HrJRHMEZSTMbHQ6FCdg2+4Bkrb1MvEiKj4jBcHj7OzB5Ze2c1OnoUWraUz8u2G97B9sc21wffM2fghRck8A4eLPeVKCHzvIMHy6KerHzAe3vLcPqYMfD331KpHmDHDhmVadxYgr6Hceme640b5cJlxozLvzs3Cr43owFY5Q4zZ0qijRkzHHJ4T1h5a4+LBHvmL46Oho4tYhjzb0dWfr6L2263sjXsnG3GyKr4KlVg6lTJjLRwoWPOZVlQsKD8XLUqfPGFFJe4916ZJ9692zHntTN7XNRliTHw6acydJKcLMk0nnjCsed0AA3AKueLiIBXXpEPNwflePaEfLf2uEiwV1aj+Hh4pF0CQ3d1pI1ZSPW8hzL1eodYtAi6dJF9oZs2yX6ovHkdf15vb+lx//knfPAB/PYbNGggxUHcnMtGfj79VOZ4H3pIVp97aEEOXQWtcjZjJOgmJcmQn4Pq+jp05a2d2Gt7TnazGkVHw9vBaxm97UVqsQOmTiPDte3sLSFBepu1a0sbvvsOHn3UNfWf8+aVzEwvvCBDq/nzy7/fL76QC4PUHrMbcfrIT1KSXLA89xwUKCD/tz1oyPla2gNWOdukSbJdZPToTGVTyiy3y3ebBncof3fsGPStu4rPtjWR4OvjIykBXWH/fpk/bNpUrgq8vGQY0xXB90rFisGDD8rPW7bIVrmqVR03HJ4NThv5MUZGJBo0kAWUBQpIEPbg4AsagFVO16SJDFX9738OP1WHOoGsHdCUQyNas3ZAU7cKvuD6i4S//pIpu3KHV2Pj4khBSkqmslzZzcaN0pj//pN1AW7YuwSgbl1pa/Hisvq6Sxc4edLVrbrEKRd10dHQsSO8+aYs0EtKst+xXUy3Iamcyc22IuR2W7bA0y2OEpfix4KP9nFn72Yy/OvrK1mLnDmHt3QpdOoEt9wiP9ux9JzDJCTA8OHw/vuSRnHrVrf59+3Q/e9//CF/q7//lh7wq6+6zfvOKN0HrHKfjz+WnsP06W6xGT83W7EC+rXbx4ILLSkSVIW8a36RhBahoZnOcmUXPXpIAPv5ZylH50n++EOSedx/vyT7OHNGLiRyImMkocbhwzBnjsemjNUArHKX3btlbu/BB2HevAxdMeeGLFauMGcOfN51HQtNGwoUtmH75ScZVnWF6GgZak5MvDyP6Mk+/FAuNMeOlaFpD+sZ3lB8vAwz58sHhw5Bnjyu3Z6WTZqIQ+UeiYmSCzZ/flk9msHg65K9jDnc+PHwzeOLWZ7SlAJlC0mGK1cE35QUyZRUrx6cPi0Lvzw9+AI8/DDccYdU9WrfXrbbebq//5aEJM89J7crVPDo4HszGoBVzjJihOzhnDgxw/9xPSGLlSc5dw4+fiScoy+/z1T/nvjWvlOCrwNXod9QQoLUfP7kExkRKVTI+W1wlKpVJbXl6NHw668yN/zjj65uVdb9dHF05M8/c02ObN0HrHKO8+dhwgQZjnvkkQy/zJF7GVNSYN8+mY7euBFiV4RT+Wgou28J5liFIIoVkwWuxYpx6edbb5URdE+bujZGUhx/+/xKvj7VBj8rAS98sD762jXzlDExsoBn2TIZrh0wIOcM06ay2aR336GDFIsoW9bVLcq86GjZqTB1KtSsKRcRt93m6lY5hQZglXPkzSuLa3x9M/Uye9aPPX5csuJt3CgVfjZvvpzQqL3fL8xNaI/NJJIc7cO4uBGEJjfmt7N3cCRGhkQbEE4woQz0CcZ2bxBNm8p6m/r1M/a2XDWXffiwJBtj0UJmenUjgHgsYyAR+UU0c0H58L59ZQXY1KlSejInq1RJLjRS9ekjI0BvvOH+9Yjj4mDxYhg4EN591/OuPLNBF2GpnGHjRhm+8sr8rIo9KhmdOwcjR8Lvo8JpmLCSf2wVqVg2mSb5t3Ch9cOU69qYqp/9D6+vJl//4kWLuNC8DXEffkrBoX0BSLG8+aLEu7x/7HmOcSt588oi0NSAfNdd1+eLcEVFpqQkqS0wedA/jE7sTdvkBZgKFbGORMqDrthmlOrUKVi/XtIV5ibJyZJQ5IcfpGrTV1/J/Lc7OX1a1mgMGCD/Z8+dk3UbOZCuglY5265dUKcODBokV9BZkNWeY3Ky7HQaNAjqHFvCYqsdlknh0kCnv7/MP774omy7adVKApOPj0SukiXhnntkiPbll2UI/Zr/k79+vJ0Fh2qyZ+k/7PvTojQRtM4TSkqTYKp2D6JFC1nc66havTeyfr3kNym/PYTvbE/i52vwem+I9L42bcrWNqMs9+RTq1517Qp+fpk+b44yfz706iXJRvr2lf8b7hDkFiyAnj3hxAmZw27QwNUtcigNwCrnSk6WjEYHDsj2o+LFnXbqlSvlc23v9njuCvJnds0PKDNpkDzo5SUB9eOPJXdtqvT2v4aHy1BtaoKKzz+X+bHeveUYL78M48djsDBAIj60ZjGrvJtz332wld0EVDqOT5HzVx02O7V6rxUfL+9707hwLvy8kh3FmtJraAlarOyPNXq0FKfPpmz15IcMgffeg8mT4fnns90WjxcVJRmkZs6Uaku33uq6JDUnT8q/5e++k575tGly4ZzDaQBWOdeYMdLjmjVLFl85wb590K8fhC06ycj8H/CE74/kObQba+eOqwNoVoZe0wvQe/bIiZcsuXRXUt6CvPPyGRYvsTizM4IjlKJR/t9olm8p4WVrsKNqecrdlkT42/dn+f2ePCmnXLgQVv0SxzOx4xjBQEknGRCAlcn3ebPebZZ78hMnSt7k7t1hypSct+AqO44eldEWY2RIvnZtuXp04gUrLVpIpad33pGh50yu1fBUGoBVznToEFSvLsFq8WKHf+DGx8MXT4cTPXcZZWxH6Gr7Ht/EGKwePWT7U9Gijs/wdGUv2dtbPszefhuMIeaW0iRExVAg6TwWKSTiSwuWss73PurUtqhbl0tfVarI4ZKSZBAhKenqn6OjZU3PggVS57xZyjJe8/+CpklL8UuKxSA9a2w2GDZMFtBkQEZ6txUGLCGtT6V0e/Jz58Ljj0tFo3nzrh51UJfFxkomsDlzICBAhoLfeEOCsz0ZI8PLX3wB48ZB4cJyQVq8uKx0zkU0AKucafNmGWZcsADKlHHoqWJiYOD94Xy0qSn+xEvwue8++YCpWtWh575OWkE+KQm+/Zbo90dQ4M89l+agD5W9m/GPrmf7xgTqr/+c9Rdq40MCddhGKMGsQ15vkUIeYvElAT8uEMxvPMksFlXoTYmnWvLc2U8oPfcTrPbtJXfyW29lqaefkd5tpnvA58/LKuDbbpOrhjx5MtSWXG3vXtma9e23crHy88+yui+7YmPlmOPGwfbtsu86JESKouRSGoBVzuWE+azTp2XUrtmG4Qyz3sErJVnmeN9/P8M9P6e5sodss8kFQvfusHMn1Khx6Wmp/+vXNh/CxgcHU/bwGjqNue/64/n6SrC/6y75OfV3ncWefkZ6t1maA969W3pxhQtnuC0KWTsxZoyM4OTJI0F5yRJZGBUUJN9Ll077tYmJMrSdkCAXP6dPy8XZ6dPSy33lFZkWyuUXROkFYB2nUZ7n2DFZoPT22w7/z330KLRvFkPffS9QqU8bvL7wvdzzCw526LmzJChIeqTXBsfq1WXVaf/+MG2a7NG1LO69/Tj3vgZEVoTSo+R9rVwJixZJFpHkZDnWtUE2KChLQ+wZ2XOdGmRvugr6r79k6uHVV6FatUy3RSEjB2PHXr5dvLhcZI0fL6v3QeYrdu+W+3v3hnXrJO3lsWNyAdy0qfybK1JE1mMEB8ueOZ2DvyntASvP8+ijEiC2b3doMfdDh+DxpieY8E9r6rIZa8YMqFjRdVV87OHaldZpDR9n5DlZZLe9yidPSu8sOlp69zk4X7BLJCTI/6916+DsWbnYBcnRHBEBgYHSMy5dWv4P3pfG6IkCdAha5SQhIZKE/oMPZB7SQXbvhmfvP8jMky0p7xOJbe5saNvWYedzqowMHztwMVm2s3VduAAPPCDJV0JDc/w+UuXZNACrnCEqSoYab7lFPnwdlGJv40Z4ucV+Fp9tTOH8SXj/vNgze7s5kTFS7eqbb+D772Xls1JuTMsRqpxh4EDJ6jNlisOCb2ioTGlFFyyLf9sWeK9bq8HXnWzaJHu+hw7V4Ks8ni7CUp6jd2/JnOOgmrIbPw/nSJ/xPFH6KYasbUn+wG8cch6VDfXryxBFLsigpHI+HYJW7i8pyeGJFY7OC6dopyb4kojx9sZavVp7vu5k/XoprpDbCisoj6dD0Mqz9eolBbpTUhxy+Lg4CHt+Gj4kAsgWndBQh5xLZcHhw9CunWxxSUhwdWuUshsNwMq9LVsmifXLlMlSqcGbMQbe7nKIZqdny75Fm8199/jmRmfPSnrJCxck41kuyR+scgedA1buKzpa9h1WqSKLbhxg/HhoEfIifv5eWF9/L5mBPHWPb06TlCR1bffsgV9+cX7KT6UcTAOwcl+vvw6RkRAWJnV17WzNGnjtNejSfBot3jkIjRvZ/RwqG+bOlRzFkybJvl+lchgNwMo9nTghSTfefFMK1tvZkSMwsd3PVCrfgrFzS+JV0M7VYFT2PfGE1K+1R5EApdyQBmDlnooXh127pJqKnSUkwKQm3/JtVFeOvvQ5BQu+bPdzOEq2s0h5gvnzZdqhalUNvipH00VYyv388osUAShRAvz87H74TztvYOBfPThR9T5KDn7B7sd3lNQ8ypFRcRggMiqOgfN2ELI10tVNs5+VKyXBhgPTjCrlLjQAK/eyeDE8+CB8+aVDDj/7k0iemteB2IIlKb76R49aVTtq6b6rihgAxCUmM2rpPhe1yM62b4cOHeCOO2DaNFe3RimH0yFo5T5On4YXXpC6td272/3we6aEcffrXSniFYV36AYoVszu53CkI2mU8Uvvfo9y+LBceBUsKCMgDph6UMrdaABW7uPVV2Xx1eLFdh96jlsZTvnnH8CXC3h5e2PFnbPr8Z0hI7V0PdawYZIR5fffb1wAXqkcRoeglXuYPx9mzpS6o3fdZffDr+6/GB+TgI0UrNQi8x6mX8vKBPjYrrovwMdGv5aOq4nsNOPGwapVcOedrm6JUk6jPWDlHkqWhPbtHbL4Zt28IzTY9Llk0rLw2ExXqaudHb0K2mkrrZOTpef76qtQuDDUrGn/cyjlxrQYg3ItYyQFpIPEnjdsKv4g9eNXY02bhv+Rg5rpKh2pK62vXOwV4GNjeMca9g3CxkiO74kT4euvoVs3+x1bKTeSXjEG7QEr13rtNcm/PHq0QwLxkjYTeTRuKfv7TOCOp7V+7M2kt9LabgHYGHjjDQm+b76pwVflWjoHrFznxx9hzBgZinRA8N00ax+tQ99gd9lW3PFJT7sfPydy+ErrpCR49ln45BN4+WUYPtw+x1XKA2kAVq7x11/QowfcfTd89JHdDx8bC1/2288p7xKUWzHVocPcOcmNVlTbbaX16dOSbGPwYBg71iEVrpTyFDoErZwvPh4efVSGnufMcUgyjLffhslH29JleSvK3OZjl2PmhjSQ/VpWTnMOONsrrWNiICAAbrlFEm4ULJjNlirl+TQAK+fbvBn274fZs6FcObsc8srg2GD/SfLPv0Cvl56nyQP2C75XBqbUNJBAjgrCDllpffIkPPQQ1KsHEyZo8FXqIg3AyvkaNYJDh6Q3ZAdXBke/2ATeW/A+/lY86zsEA3fY5RxOWZzkJjrUCbTfe4qIgBYt4OBBeOcd+xxTqRxCJ2CU8+zZA998Iz/bKfjC1cGx3/c/clvKAfoED+SzzRF2O0eOTgPpKPv3y8VWRAQsXQpt27q6RUq5Fe0BK+c4f17mfY8flw9iO+b6TQ2CTy/7hedOzGJO4bZsubsclh2DY45OA+kICQnQqpWklwwNdUh2M6U8nfaAleMZAz17wu7dMGuW3RPtlyoUwD0HdzB46zgM0PbcUu6K3GPX4Gi3NJCJiXDsGJw7J9uvcprYWHlfvr4weTKsWaPBV6kb0ACsHMsYSbYwcyYMHQrNm9v9FP1aVqb6hpMYvLAAn+Qk7o3cdVVwDNkaSaMRK6kwYAmNRqzMdA3dDnUCGd6xBoGFArCAwEIBmcsO9cMP0KSJLEAqWRIKFABvbwlYICkZq1WD+vXldzRgACxc6FlBeulSqF5dEmwAPPAAVM4BeaqVchAdglaOtXGjZLnq1Uv2BjlAVVsxPo7oyOvWeHxJIMnbh7rdHqbJxeBorxXMGVqcdOgQLFsG4eEQFgZLlsDtt8s2nPh4Kbd4++0yNHv+vGzNAakAdOedct/x45Ko4quvpDoUwJQpksSiYUN5njvtnz1+HPr2ldGNKlWgdm1Xt0gpj6C5oJXjrVwp+ZcdEDRMzHn+veUuxvIqb82tQ5E/Qq/L9dxoxMo0528DCwWwdkBT+zTk5Em5wJgyRXqtxYtLsHz/fekVZlZcnATzatXkdqNGEtBBes/33w9PPCFfrjR/Pjz3nAypv/UWDBxo91KSSnkyzQWtnO+HH2Sl8333QVM7Bbk07H1kEFXj9nPfGzUo0joIWl9fZMEpK5iNkdSaL70Er7wCt912XfatTCXyCAi4HHxB6uQeOCA96zVrZLi3WDEJwMZIoG/SRC48fOyz9zldiYlynkKFpJ2TJl3dXqXUTWkPWNlf6paT4GD52UFpIKOXriN/q4bML/EiDx8Zf8MOtkN6wAkJEnSWLoVFi+Q9xsRAvnxpPt3uVYaMkV5ynjyyx7ZyZRmiLlhQ5pAfeABat7ZfcfvkZNiwQd7r4sUS6CdNutwWTfWpVJrS6wG70USSyhHCwqBjR5mnnDPHcR/MFy5w/olniaA0lecNT3d0266F7FNSJINXtWrQu7csojpzRh67QfCF9BN5ZIllSfAFqFgRTp2SHvgjj8jfoGdP2LJFHt+6VeZoZ8+WYe3MXnQPGCALxxo2lLzdRYvCPfdc3RalVKbpELSynz/+kF5XYCD88ovdtxtdaeektVSO2s+Mjgt4tmGBdJ9rt/SK//4rAW7DBike//PP0LJlhgKQw4fBCxSQC5+OHSXAXplpbNcuWZn86adyu1gxmZf+8UcoUkTq8X777eX3YVkyjzt/vvxsjPSo27aVvb2FC9unzUrlchqAlf1MmQJ588Ly5VCiRIZfltkiB4mJ0PnLpuQrdYBfZ5TN0Dnskl6xUCF5f19/DV27SjGJDHJqIg/Lkl5xqiefhMcfhx075OJhwwbJUpXqwgWIjpafjZEvmw2OHoVSpWDkSPu3USmVvTlgy7JeA54DDLAD6G6Mib/R83UOOIdKTpYP7JQU+dAOzHigy/TcaHIy374cRtcvGrNgAbRrZ483cBMrVkCDBhJ8szjfafc5YKWUR3DIHLBlWYFAb6CeMaY6YANcvCdCOd2iRVCrFvz3n2wzykTwhczPjZ4eMpYuX9xH//vCnRN8x42TYgLDhsntLM53ZjuRh1Iqx8nuELQ3EGBZViKQBziS/SYpj5CSIltfBg+WVIOJiVk6TGbmRs2Bg+T98G2W2NrS65sGWTpfhqWkQL9+khCjfXu7VPKxa5UhpZTHy3IP2BgTCYwG/gGOAtHGmGX2aphyY2fPymKfwYPhqadkj2oWt7vcaA70uvvDwjjfqDlJKRaRb02gTFkHrryNi4PHHpPg27u3LFbKm9dx51NK5UrZGYIuDLQHKgClgLyWZT2ZxvNesCxrk2VZm06kptVTnq1/f9kL+tlnsiApIOsLiTK0RSg8HBMcTL7/DuJrJdKj+b9ZPl+G/PefXFR89hmMGZOpxVZKKZVR2dkH/ABwyBhzwhiTCMwDGl77JGPMZGNMPWNMveLFi2fjdMrlkpLk+wcfSHrJV1/N9h7QDM2NhoZiLs4Te1speP8emq1z3lDqSuDy5WWV8KuvOuY8SilF9uaA/wEaWJaVB4gDmgG6xDknSp3vXb5cVgQXKSIpJu3kZnOjs+Jr8jB++JBAss2b9SXvpIndzn7RmTOy0vmRR+QCo0D6e4uVUiq7sjMHvB74AdiCbEHyAibbqV3KHRgjySbq1ZP53goVJBg70foPJrD20wO08P+JTxo+RZfH36fnAd9MlxNMV0KCzGn//bckmlBKKSfQXNAqbUeOSKL/NWsk8A4bBl26ODft4NGjRJetwo6kGjzS/hP8q1xeQ2C3SkbGwLPPwrRp8M03krRCKaXsRKshqYw7e1aGX4sXl+A0fryUm/P1dW47jOH80y/im5TAS+VG4lf56gV8dkvhOGKEBN9339Xgq5RyKg3ASuzfL8PMq1fDn39Kov81a1zWHPPd9+RdvoB+tpGcfige72s63nZL4VilivSAhwyxz/GUUiqDNADnZmfOQEiIbClasEAS8L/2mtPnea9z7hwJ/3uFrdxDVJ9nyR+wkbgr8nxkuZLRlWJj5SLj4YflSymlnEzLEeYmKSmSiH/vXrn911/QowesWycJJw4elNXO6ZTVc4bTifl52jaTT6pP44uRRe2fwvHvv+GOO+C77+zVZKWUyjTtAed0GzbAvn3w66+yovnECakVO3Ei1K0rtWJr1XKfmq6xsfTvn4cfYlqxeabkwLBrCsfoaGjTBs6fhzp17HNMpZTKAg3Anu7PP6X268GDcOCAfA8MhLFj5fHHH5ceX5EissXmoYekhi1I8YTatV3W9OucPMmFarXhxGD69nueWrXsfPzEREkxuW8fLF0q879KKeUiGoDd3eLFEB4Ox45d/ipQAH77TR5//nlYtUp+9vWVOrBXZhz77jsoWFCGXN05pWJ4OCn/exGv4t7uOgAAIABJREFUE8f4t1QDPhvsgHO8/josWwZffQVN7bCFSSmlskEDsDs4eRK2b7/8deiQBFXLgjlzJIiWKCFft94Kt912+bXDh0vPrkIF6fl6XTOt38DBVYPsITwcgoPxSkggCW+GvBFj/9oHxkDRovDGG7LqWSmlXEwDsKu99ZYE0VSlSsmc7Pnzshjq889ln+qNeq9BQc5ppyMtWIBJSMACbJahQXwoYOf3ZVmyzUoppdyEroJ2ti1bZNh461a5/cQT8NFHkmf5+HGIjISffrq8ErlgQfceOraDZP88ACRhw8vPF4KD7XdwY6B7dxl6VkopN6I9YGeIi5Oh5IkTYf16Kd/XuLGswq1ZU75ysQ+83mUNDRj9+GZqvRps3179uHEwfTrUrw8tWtjvuEoplU0agB0tORnuvFPmdatUkfqy3bpBoUKubpnr/f47+zdEMXRoGx7v0oJas+wcIHfuhH79oHVrePFF+x5bKaWySQOwo/z7L5QuLcPH77wjNWaDg91nv62rnTxJyuNP4H0iL2VLtmT8eB/7Hj8+XopHFCwIU6bo710p5XZ0DtjejIHJk2Xbz8yZcl/37nD//RoEUl2cl00+doJHEr9nygwf+w8IzJwJO3bIArYSJex8cKWUyj7tAdvT2bPwwgswezY0b65zjjcyZgwsXkxfxtL09Trcf78DzvHss1C5ssy1K6WUG9IAbC+bN1/OOvXhh9C///V7chUcOIB5802W+rVj1e0vs/EDOx//5Ek4d072RWvwVUq5MQ3A9vLvv5CQIAk0GjVydWvclqlQkTHVv2LkztYs+9bCz8+eBzdSuzgsTFJyuriohFJKpUcDcHZERUnAbd8eOnSQHMsBdqpTm9OEhcGiRSy22vHa1m6MHg01atj5HF99JWUVP/5Yg69Syu1ZxhinnaxevXpm06ZNTjufQ8XFScDdvBn++UfSHKo0rZq+kIbPdcI7OYkL+PFclaXM2NXEviP0Bw9KRG/YUAot6PC/UsoNWJa12RhTL63H9FMqK5KToWtX+P13WWWrwfeGQrZG8t+YiXgnJ0mqSZKoXXgmC7dH2u8kxsjiN5sNpk7V4KuU8gg6BJ1ZxkCvXjB/vqzmfewxV7fI5UK2RjJq6T6ORMVRqlAA/VpWvlS/98evf+aLXaEYLJLwItHLm7BKlZm/dJ/9avxeuCCLrh55BMqUsc8xlVLKwTQAZ9bSpTBpkqxy7t3b1a1xuZCtkQyct4O4xGQAIqPiGDhvBwAdbvVi1JT+nLXlp2viLGrdspEdLYqyJbAqVlSc/Rrh7w9ffmm/4ymllBPoWF1mtWwJixZdXcEoFxu1dN+l4JsqLjGZUUv3QfHiLK/UkpZJS/m5ZFOmPXkfWwKrAlCqkB0WqxkDAwZATllXoJTKVTQAZ9TPP0tuYcuCNm00q9VFR9LoyQYkxJMQcYSIY970OfINe/JU4ZaOm/DySZHHfWz0a1k5+yf/8UcYORJWrsz+sZRSysk0AGdEWBh07Aivv+7qlrida3uy3slJjFs4knnf9qdT63gS422M/vIsZUt7YQGBhQIY3rFG9ud/T52Sufi6daFv3+wdSymlXEDngG9mzx7p8ZYpczm3s7qkX8vKl+eAjeHDpeNodmAjw28by6ad/ixZAq1aleBV7JyPuW9fOH1a6vx66z9jpZTn0U+u9EREyJyvn58svipe3NUtcjupPdlRS/fRZcEXPLbjV76rPpC3dr7C+PHQqpUDTvrbbzBjBgwaBLVqOeAESinleBqA0/Pee5LtavVq2eai0tShTiAd5k+CdXM5WKM9XXZ8QJ8+8NJLDjpho0Ywdqzs/VVKKQ+lmbDSk5gIe/c6IGdiDhMeDs2aYeLiicOfdxutYOSqIGw2B5wrIQF8fR1wYKWUsj/NhJVZ+/fL/KKPjwbf9Jw7BwMHwvLlmAsJWBh8SeCDB0IdE3x//x0qVYLt2x1wcKWUci4NwNeKj4eHH4YHH5R9piptf/0FDRrAqFGciM9PvPElERte/r74tQy2//ni46XSkc0mQVgppTyczgFfa9Ag2L1bFl3pXt+0/fILdO4MNhs7Ri8leFgz7i7QgClPhVKqSzAEBdn/nEOHwr598nfRSkdKqRxAA/CVVq2CTz6R1UMtWri6Ne5p6lTpidasycLu83n0zQqULw/jfgqiVCUHBF6Abdvgo4/gmWf076KUyjE0AKc6e1Y+4CtVkg97dVl4OISGQnAw3HsvpsezjAr8jP598nLffVKXokgRB55/xgwoVkzq/CqlVA6hATjVhQtQvTq89Rbkzevq1thNepWKMiQ8HO6/X1Yf+/uT+MsKXkj+kulD4ckn4auvZJu0Q338MfTp4+Aor5RSzqUBOFXx4lJkIQdJt1JRRoLwjh3w/PNycQKYhAS+7h7K9INBDBkC777r4Gnygwdl0VW5clC2rANPpJRSzqeroE+cgEcfhcOHXd0Su0u3UlF6jh2TOsc1a8KhQ+DtjbHZiE/x5evDwcyYAYMHOzj4pqTIlEDjxrIfWymlcpjc3QM2Bv73P1iyRLpzV8j20K0bSKtS0bX3h2yN5KfJ87ht9yYO31GL5i89Toc7CsHWrTBoEKbPa/w+ZR+rh4ayygpm2KIggoOd0PhJk2DNGln05ePjhBMqpZRz5e4APHOmrCD66KOrEm5ke+jWTZQqFEBkGkE4tYJRyNZI5n72HVNnvYVPciKstngqIQX6dKbD3r1s2mrjjU6walUQVasGMW8eVKnihIb/8w+8+SY88ID0gpVSKgfKvUPQx49D795w773XlbPL8tDtNUK2RtJoxEoqDFhCoxErCdkame1mZ0a/lpUJ8Lk6JdWlWryrV1O062NMmfUWfsmJF/8hGOoe3Mqw7w/TtZuN+vVlS/SECZJ8yinB1xjo2VOGoCdP1r3YSqkcK/f2gEeMgJgY+PJLrs2bmJGh25txh150hzqB2OJi2TRmGkFbf6NISgLn33uf++sEwsLNlDl6iNUV7iL44Ga8TAqJNh9+OtWBLR/fw24fWRDevz8UKOCU5oqEBCn9OHy4FsBQSuVouTcADxkCTZum2a272dBtRqTXi04NwBmaZ75yD+61Gabi4mT+evlyKFUKAgIkPeR998GBA3D33bQ9fZq2V77m3+3Ag9C2LV375yMyKo46/+6lzvp/WfpvR8L2NKb4XcfYHFKSMmUy/Hbtx89P5n+VUiqHy30BOCUFkpOlW9emTZpPuarI/EWXhm4z6Ga96Az1kJctg3btpFdos8kFQ6dOUoYvJgby57/+BO+8IwG4RAlZyXzwIPz6q7xvm+1SfuvosxZB1ObLxfEsPNCckHhf/MudpHzzMD59qfyl4OvUxWjvvAPt20O9NAuHKKVUjpL75oBnzZIi7keO3PApHeoEMrxjDQILBWABgYUCGN6xRqYCz416y6n3p/aQ74rcw0vhc6j37y4qRuxn2Rdz5YkpKRJ8L1yQoJmUBBs2wNGj8ni+fNCy5eU5UptN9gYNHXr58YkTpafv5wc2Gyk+vsz+L5gHHriYWGpgEYgsQdFqpyjx2Hr+396dx0dV3nsc/zzsCEJALBhcALFUubaiqRrX2Iq4luVWC4qiFxSL4nIVgVJBxKooLq316lXRVqvIS4qyK2sQFBUUUBCQkCsqYEDZXBBI8tw/fhMIIQkhOTNnzuT7fr3yyjBz5pzn4WTOb57lPL8OfT/msX6t9mmhDx7/Ceu27sCz90tCXMayJ06E++6DadOC37eISBKqXvmAv/sOfv5zW9RhwQKoEb/vHyVbuGCt6KJA3nrQFM5ds4hnx4+gVqFt44BPmh/LiV/n2BuGDrUZ2vn5lgN31qx9u6FjeXj35Mgt9vrmzZbKeOVK+H7GAgpnZzN2Yxbvkcnxx1tsv+wy67EuK3XgmQ/OLrUrvmVafd4Z9JtA/p8AyMuzWejp6fYlQ/l+RSRFlJcPuHp1Qd93ny0yMWFCXIMv7O1GLqv7Nj2tPoOzn6d2LPgWApN+cTYvXtaXcUU7ufdeS4tYxhjw7oxMNrwwix+nZrO4cRZzXshk5UALups27d2uTp1MzjgjkysGwUuXQdu2FatDEJPRDsh7uO46+3L0yisKviJSbVSfALxqFTz2mF3sTz01IYfs0qHlvt3W774L3frDww8zoFM7nn2vO/dNfoxahQXsrlmLMad3pWePrH32UXBqJqubZJKTAzmPw+rVloo3J8cW7yooyAQsMDdrZnPKOne230U/rVqV3cotTxCT0Q5ozBjrdn7ySTjhhOD2KyKS5KpPAH7ySZsl/MADiT3u/Pk2Frt0KSxfDk2awIoVdLn0Uhh2C7cccSRtP11EzgkZ/OGGbnQ+qSWffWbzpmbOhDlzYOvWvbtr3BiOO86+Q1x5pbVm27aFdu0sAAcpiMloB3T55TYprmfP4PYpIhIB1WcMOD/fAuCvfpW4Y86ebas5eW+TpW6/3bqVS2Rbysuz4duioPvll/b8McdAx462HHK7dhZomzZN7NoUcZsFvWMH/PBD8N8aRESSSPUeA965E3780VqeiQy+AO+/v/dxjRoWbGLB13sLusOHWyMZrIi/+Y0tgNGxI7RpU36wTcQtQvt1owdl4EAYP96+FDVuHPz+RUSSXOrfhvT44zbz+euvE3O8H3+E3r2t7zgrC+rVswHYOnUoymIwb56l2O3YET7/HP7yF1i40CZOjRtnKzEee+yBg2/CbhEK2tSp8MQT1v2s4Csi1VRqt4DXr4cRI6wbuEWL+B9vzRpbKGPpUptQdMcd1syNzWL+oGYmd3ey9TVatIC//c3S7dard/CHqshKW0kpL88mwv3yl4kfjxcRSSKpHYAHDrSx30cfjf+xJk6Ea66xruYpU+Dii+35zEyW1M9k6FCYNMl6oR9+GPr1g0MOqfzhEnKLUNCKbjnavt3GxyvzzUNEJEWkbgBesMDSDQ4ZYoOp8TpGdjakpVlEPeUU60Nu1QqALVvgj3+EsWNtk/vuswRMpa0gebAScotQ0H74wfrVR42C9u3DLo2ISKhSNwBPmmT9vIMGxWf/JVehGjDAZjjHWnU5ObbUdG4u/PnP1hudlhbc4RNyi1DQGjaEyZPDLoWISFJI3UlY998PS5bYRT8eZsywW2kKCiwIN2myJ/jOm2dLPG7aZLcVjRgRbPCFYNarTpgtW+Daa23GmXPK8SsiQiq2gHfvhq++slyyzZvH5xgFBTa5CmzMt9gM55degj59rBd6ypSKL/tYGXG7RShIO3dC1662Cti11+7pnhcRqe5SrwU8erStWrF8eXz27z3cdhu8/bb9vu8+mDWLwtMyuftum4d15pnw3nvxDb6RUFhok67mzoV//GPPlxQREUm1FvB331lKvtNPj9+6wqNGwd//boO6o0YB1hN93ZU22ap3b/if/wkmp0BCc/HGw5AhttbzAw/YupkiIrJHagXgRx6BjRvtlqB4jDPm51u/8hVXWJpA7LbWLl1s0auRI20uVhCHLpnOsGihDSAaQfj77+089O1rt4OJiMg+UicAb9hgLdLLL4fTTovPMWrVgjfftMc1apCTY2t8bNwI//63DXUGJbILbRRp2NDGfRs00KQrEZFSVGkM2DmX5pwb55xb6Zxb4ZzLPPC74mTuXBtzvP/+4Pe9bJllsN+82WY616vHtm12m9H339twcJDBFyK60AbYmprXXgs//WTLTNZKne94IiJBqurV8a/Am9773zvn6gBVWNupirp3t8WVDzss2P1+9RVcdJEF9++/h6ZNKSiAHj1s5clZsyCj1DwXVRPJhTZyc+1bySGHwLZtWulKRKQclW4BO+caA+cAowG897u891vLf1ecfPqp/Q46+M6YYdF182ZLIHD00YBlK5o2DY6+dCW9pk7hzAdnB54EYUCndtSvXXOf55J6oY1vv7UvKrt3239OvG4BExFJEVXpgm4NbAJecM4tds4955xrUHIj59wNzrlFzrlFmzZtqsLhyjB/vi1rOHZssPt991248EKbZVVQYFmOsNUtH3oI0k5ZS0G7NXHLRBSphTa2b4fOnWHtWpgwAX7xi7BLJCKS9KrSBV0LOBno771/3zn3V2AQcHfxjbz3zwDPAGRkZPgqHG9/3tu04/R0uOyyQHfN5MnW7Qw2+zk7mw9qZtKnDzRqs4VG5+17n3E8JkhFYqENsAlwK1bYKiRnnx12aUREIqEqLeCvgK+890VZ58dhATlxxo+3FS/uvbdqqYVKc9llUL/+nly+m9pn0bUrHHEENL5kIa7m/t8lkn6CVNA+/ti+BLVrZ+O/l18edolERCKj0gHYe/818KVzrmhQ8rfAp4GUqiJ277ZEC+3bQ69ewe23sBCefho6dLAZViNGsHPqLC79SybbttmtrUell95xkNQTpILkveVU7NDBWr1gM55FRKTCqjoLuj/wcmwGdC5wXdWLVEHLlsE339igbCVudSlzlamnnoKbb7aA0qMH/vRMru8FH3xgDe4TT4QB+RHMRBSUH3+0xa7HjLEFSf7zP8MukYhIJFUpAHvvlwBxuAmnAjp0sOw6jRod9FvLWmWqwRe5dLzrLpvN2707YItrvfSS9XIX3etbNC4b6WUiK+OLL2zZryVL7H7rQYO0yIaISCVFe5WESnZ7lrbK1M6du2h+651Qty489xw4x5tv2iqKl19uOX2Li8wEqSAtWwb/93+Wa/mSS8IujYhIpEU7AFdSaZOl+ix8g1+uXQ4vvwzp6WzZYol82reHF16oxg29DRtg9my46iq4+GILwEEnNxYRqYaqZQAubZWpOW0ySHe7uLZHD8Bavhs3Wu6FBvvd3VwNrF1rNzyPHm2Tri64AA4/XMFXRCQgqZcPuAKKrzLlvN3r+1V6G9IeGQnOMXcuPPss3H47nJzYG6vCt3695VRs29b+E3r1snt8Dz887JKJiKSUatkCLj6J6opJz3H8D3nseGY0nTu05KefLINeq1YwfHi45UyoXbv2JjEePx769YM774Sjjgq3XCIiKapaBmCITaL66E1Y8KotOXlqK8Am965aBW+9leJdz99+a8ttvvuuLedZu7aN9aanw7p1wS9sIiIi+6iWXdAAzJsH119v45vZ2bBgAcuXw4MPQs+eNuSZMgoLYfXqvf++6SZo1sxSLD7yiC1qcu65tuQmKPiKiCRAtW0BM3KkBV+AXbsonJPN9ZMzadQIHn003KJVWGEhbNpk47YbNsA550DDhtZ8f+45+Ppr+1m/3hbQyMuDn/3M0jYefTSccYZle6pfTVbwEhFJItUzAOfnw4cfQo0adn9RnTq8vjmLBQvgn/9MovlGhYWwfDksXWrdwn/4gw1OT55sY7QbNuxttQIsWgSnnGJBedkyaNECfv1r+33iiXvz83bpEkp1RERkr+oZgGvVsqA2fz6sWsXGE7K47upMzj8frr467MJh3cUDB8Lbb9tYbZHjj7cA3KIFnHcetGxpP+npliWiKA1gz572IyIiSav6BeAtW2wFrZ/9DLp1w3u4oas1JJ9+OoQFNz77zG42zs62Mdneva0beckSy8iUlQWnnWazkYtmhWVkWFNdREQiq3oFYO/h97+3Mc/JkwG742bCBBsSPvbYBJZl2TIYMQJee83K1batzcYGa83m5iawMCIikmjVKwBPmGC32vz97wBs3Qr9+8NJJ8F//3eCy3L99RaEBw2y8dwjj0xwAUREJEzVJwDv3Al33GGLO/ftC1jsy8uzHL+VyGh4cBYvtmb2E0/YLK/nn4fmzaFp0zgfWEREklHKBuCS+X6f2TCT9rm5MH061KrF0qXwzDNwyy02pBo3H31kS2pNnGhjz0uW2G1Axx8fx4OKiEiyS8mFOIry/a7bugMPrN/yAzXGjmXDuRdY8AMGD7Z4OHRonArhvbV4MzJsNvPw4Za/OHZ8ERGp3lKyBVwy3693NejScxRt6xUyBZgzB6ZNs2Q/ce0BXrwYrrgC/vd/K527WEREUlNKBuDi+X7Tt29kc/1G/FS7Hp/mW8N04ECb83TzzXE4+Nq1UFAAbdrYrUJ16lTjZMIiIlKWlOyCTk+LLa3oPX+b+DCvjhkM3pOeVp9x42DhQrj33jiswDh3rnU5X321Rfq6dRV8RUSkVCkZgIvy/d7yzhgy1q1g/tEnUb9OLW7/TTv+9CebCH3NNQEe0Ht48kk4/3xLcvDCCwq8IiJSrpTsgu7SoSXNFr7Lme+8ggf6fDSRE2/syacftiQnByZNgpo1AzrYzp2WXWj0aFu56l//gkaNAtq5iIikqpRsAQOcNX40DnBAvcJ8Tv98OcOHw9lnwyWXBHigggK7tejPf4Y33lDwFRGRCknJFjAA3323T7ajl9dlkZcHr78eUO/wjh3WjD7kEHjnHRvvFRERqaDUDcDz58OsWbBwIZt/mcWt3TPp1g0yMwPYd0EBXHUVbN9uuXcVfEVE5CClXgD+5hsLkM2b26So88/nnluswXr//QEdY8AAa0r/9a8BDiaLiEh1knpjwMOGWV7c7dsBWLPG0gz27g3t2gWw/yeegMceg1tvtXUsRUREKiG1AnBuri3w3L37nslQd99tiRaGDQtg/5MmwW23QefO8MgjAexQRESqq9QKwMOGQe3aFnWxPAhjxliqwfT0APbfujV07Qovv6yuZxERqZLUCcCffGKBsX//PdF24EA47DAbsq2SbdtssY3/+A8YNw4aNKh6eUVEpFpLnQA8cyakpVnUBWbPtqeGDKliHoStW23q9ODBwZRTRESEVArAt98OOTnQtCnewz33WEP4j3+swj537YJu3Wy/nToFVVIREZEUCMDeW4CEPbkFs7Nh3jwYNAjq1avCvm+91XIXjh4N551X5aKKiIgUiX4Anj4dfv5zWxAjZvhwOOIIuP76Kuz3zTft/qU777TsRiIiIgGK9kIchYU2NnvMMXtaqHPn2s/jj1ex9fvTT3DWWTBiRDBlFRERKSbaAXjcOFi8GF580RLfY3l+W7SAG26o4r67dLH7fZVWUERE4iC6XdDz5lkawNat4corAVv+efZsuOsuqF+/kvudMsWWmCwsVPAVEZG4iWYAXrAAOna0dZ/XrYMPPgBs7Ld5c+jbt5L73bwZ+vSxSVf5+cGVV0REpIRodkFnZ+8NkAUFkJ3Nuz6TmTNh1CjLEFgp/ftbUJ82bU+XtoiISDxEswWclWUBsmZN+52VxfDhcPjhcOONldzn+PHwyiu2jOVJJwVZWhERkf1EswWcmWm5frOzISuL91wm06fDyJGVXCVyxw7o1w86dNCKVyIikhDRDMBgQTgzE4DhF0GzZhZDK6V+fXj1VdtJ7drBlVFERKQM0Q3AMR98YGtmPPAANGxYiR1s326pC7Oygi6aiIhImaI5BlzM8OG2AuVNN1XizXl5cNxx8NRTgZdLRESkPJEOwAsXwtSpcMcdcOihldhBv36WalCtXxERSbBIBuA3Fq/jzAdnc26PPGrV302rs9cf/E6mT7eZz/fcA8cfH3gZRUREyhO5APzG4nUMHv8JuStrs2NNcxpk5DJi+se8sXhdxXeSn2/N5jZtLI2hiIhIgkUuAD/81ip27C7gxxXp1Ki7m0anfM6O3QU8/Naqiu/ko49g9Wp46CGoWzd+hRURESlD5GZBr9+6A4C0rJU07LCWGnXz93m+Qk49FdasgfT0eBRRRETkgCLXAk5PsywLzkHttB37PX9An30G3kPLlkq2ICIioYlcAB7QqR31a9fc57n6tWsyoFO7A785NxdOPBEefTROpRMREamYyHVBd+nQErCx4PVbd5CeVp8Bndrteb5cAwdCrVrQo0ecSykiIlK+yAVgsCBcoYBb3Lx5MG6crdyhsV8REQlZ5LqgK6Ww0G43atkS7rwz7NKIiIhEswV80HJz4csvq5gsWEREJDjVIwC3bQs5OZXMVSgiIhK81O+CXrTIVr469FCokfrVFRGRaEjtiLRuHZx7Ltx1V9glERER2UdqB+A//clav/37h10SERGRfaRuAF68GF580WY/t24ddmlERET2kboBeOhQaNIEBg8OuyQiIiL7qXIAds7VdM4tds5NDqJAgfjuO1i71u75bdw47NKIiIjsJ4jbkG4FVgCNAthXMA49FJYssfFfERGRJFSlFrBz7kjgEuC5YIoTgNxc2LrVbjmqUyfs0oiIiJSqql3QjwN3AYUBlCUYffrAGWdYykEREZEkVekA7Jy7FNjovf/wANvd4Jxb5JxbtGnTpsoermLmzLGfvn2V61dERJKa85VsKTrnHgCuBvKBetgY8Hjvfc+y3pORkeEXLVpUqeMdkPdwzjnWBb1mDdSrF5/jiIiIVJBz7kPvfUZpr1W6Bey9H+y9P9J73wroDswuL/jG3YwZMH8+DBmi4CsiIkkvde4DnjkTjj4aevcOuyQiIiIHFEgA9t5ne+8vDWJflfbQQ7b6Vd26oRZDRESkIqLfAvYevvjCHjdtGm5ZREREKij6Afj11+HYY2HBgrBLIiIiUmHRDsCFhTBsGLRpA7/+ddilERERqbAglqIMz2uvwbJl8MorUCvaVRERkeolui3gggK45x5o3x6uuCLs0oiIiByU6DYbFy2CnBwYMwZq1gy7NCIiIgclugH4tNMsAB91VNglEREROWjRDcAAxxwTdglEREQqJbpjwCIiIhGmACwiIhICBWAREZEQKACLiIiEQAFYREQkBArAIiIiIVAAFhERCYECsIiISAgUgEVEREKgACwiIhICBWAREZEQKACLiIiEQAFYREQkBM57n7iDObcJWBvgLpsB3wS4vzCpLsknVeoBqkuySpW6pEo9IPi6HOO9P7y0FxIagIPmnFvkvc8IuxxBUF2ST6rUA1SXZJUqdUmVekBi66IuaBERkRAoAIuIiIQg6gH4mbALECDVJfmkSj1AdUlWqVKXVKkHJLAukR4DFhERiaqot4BFREQiKRIB2Dl3oXNulXMuxzmgygi2AAAFGklEQVQ3qJTX6zrnxsZef9851yrxpTww59xRzrk5zrlPnXPLnXO3lrJNlnNum3NuSexnaBhlrQjn3OfOuU9i5VxUyuvOOfe32Hn52Dl3chjlLI9zrl2x/+slzrntzrnbSmyTtOfEOfe8c26jc25ZseeaOudmOOdWx343KeO9vWLbrHbO9UpcqUtXRl0eds6tjP39vO6cSyvjveX+LSZaGXW5xzm3rtjf0cVlvLfc610ilVGPscXq8LlzbkkZ7022c1Lq9TfUz4v3Pql/gJrAGqANUAdYCpxQYpt+wNOxx92BsWGXu4y6HAGcHHt8KPBZKXXJAiaHXdYK1udzoFk5r18MTAMccDrwfthlPkB9agJfY/ftReKcAOcAJwPLij33EDAo9ngQMLKU9zUFcmO/m8QeN0nCulwA1Io9HllaXWKvlfu3mCR1uQe48wDvO+D1Lux6lHj9EWBoRM5JqdffMD8vUWgBnwrkeO9zvfe7gFeBziW26Qz8M/Z4HPBb55xLYBkrxHu/wXv/Uezxd8AKoGW4pYqrzsCL3rwHpDnnjgi7UOX4LbDGex/kYjFx5b1/G9hc4unin4d/Al1KeWsnYIb3frP3fgswA7gwbgWtgNLq4r2f7r3Pj/3zPeDIhBesEso4LxVRketdwpRXj9g19gpgTEILVUnlXH9D+7xEIQC3BL4s9u+v2D9o7dkm9mHdBhyWkNJVUqybvAPwfikvZzrnljrnpjnn2ie0YAfHA9Odcx86524o5fWKnLtk0p2yLyZROScAzb33G2KPvwaal7JN1M4NwH9hPSqlOdDfYrK4Odad/nwZXZ1ROi9nA3ne+9VlvJ6056TE9Te0z0sUAnDKcc41BP4N3Oa9317i5Y+wLtBfAU8AbyS6fAfhLO/9ycBFwE3OuXPCLlBlOefqAL8DXivl5Sidk3146z+L/K0OzrkhQD7wchmbROFv8SngWOAkYAPWfRtlPSi/9ZuU56S862+iPy9RCMDrgKOK/fvI2HOlbuOcqwU0Br5NSOkOknOuNnbyX/bejy/5uvd+u/f++9jjqUBt51yzBBezQrz362K/NwKvY91nxVXk3CWLi4CPvPd5JV+I0jmJySvq6o/93ljKNpE5N865a4FLgatiF8j9VOBvMXTe+zzvfYH3vhB4ltLLGInzErvOdgPGlrVNMp6TMq6/oX1eohCAFwLHOedax1op3YGJJbaZCBTNSvs9MLusD2qYYmMmo4EV3vtHy9imRdH4tXPuVOwcJd2XCedcA+fcoUWPsckyy0psNhG4xpnTgW3FunqSTZnf5qNyToop/nnoBUwoZZu3gAucc01iXaEXxJ5LKs65C4G7gN95738sY5uK/C2GrsT8h66UXsaKXO+SwfnASu/9V6W9mIznpJzrb3ifl7BnplXkB5tN+xk2O3BI7Ll7sQ8lQD2s6zAH+ABoE3aZy6jHWVj3xsfAktjPxcCNwI2xbW4GlmOzH98Dzgi73GXUpU2sjEtj5S06L8Xr4oAnY+ftEyAj7HKXUZcGWEBtXOy5SJwT7EvDBmA3Ni7VG5v/MAtYDcwEmsa2zQCeK/be/4p9ZnKA65K0LjnY2FvR56Xobod0YGp5f4tJWJeXYp+Dj7GL/hEl6xL7937Xu2SqR+z5fxR9Poptm+znpKzrb2ifF62EJSIiEoIodEGLiIikHAVgERGRECgAi4iIhEABWEREJAQKwCIiIiFQABYREQmBArCIiEgIFIBFRERC8P+dNESQXRETaAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.arima_process import arma_generate_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": 16,
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dmpm/opt/anaconda3/lib/python3.7/site-packages/statsmodels/tsa/base/tsa_model.py:162: ValueWarning: No frequency information was provided, so inferred frequency M will be used.\n",
      "  % freq, ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "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": 18,
   "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:                Mon, 23 Mar 2020   AIC                            501.773\n",
      "Time:                        20:09:11   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.037       0.052       1.630\n",
      "ar.L2.y       -0.2693      0.247     -1.092      0.275      -0.753       0.214\n",
      "ma.L1.y        0.5352      0.412      1.299      0.194      -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.1793           +0.0000j           32.1793            0.5000\n",
      "-----------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "print(arma_res.summary())"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}