ChadFulton/statespace_forecasting.ipynb

## statespace_forecasting.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forecasting in Statsmodels\n",
    "\n",
    "This notebook describes forecasting using time series models in Statsmodels.\n",
    "\n",
    "**Note**: this notebook applies only to the state space model classes, which are:\n",
    "\n",
    "- `sm.tsa.SARIMAX`\n",
    "- `sm.tsa.UnobservedComponents`\n",
    "- `sm.tsa.VARMAX`\n",
    "- `sm.tsa.DynamicFactor`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "macrodata = sm.datasets.macrodata.load_pandas().data\n",
    "macrodata.index = pd.period_range('1959Q1', '2009Q3', freq='Q')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic example\n",
    "\n",
    "A simple example is to use an AR(1) model to forecast inflation. Before forecasting, let's take a look at the series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1c21a29438>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAEyCAYAAABQ0omaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4ZAd5L/7vmd5HGtXRrtr2aq+9zQ1sgjE2MeAQQoxvgIATLhCSkAsh9xcukHYTCCGhhISEmNAuxNhgjIkLbri3XW/f1WqbumYkjaTpc6ad3x9nztFUaSTNSDPa7+d5eNY7GmkOWunMec/bBEmSQERERERERPVFs9oHQERERERERIvHYI6IiIiIiKgOMZgjIiIiIiKqQwzmiIiIiIiI6hCDOSIiIiIiojrEYI6IiIiIiKgOMZgjIiIiIiKqQwzmiIiIiIiI6hCDOSIiIiIiojqkW+0DyNbc3Cz19PSs9mEQERERERGtisOHD09JktRSznNrKpjr6enBoUOHVvswiIiIiIiIVoUgCIPlPpdllkRERERERHWIwRwREREREVEdYjBHRERERERUhyoSzAmC8G1BECYEQTiZ9dhfCIIwKgjC0cz/3laJ1yIiIiIiIqLKZea+A+DWIo//kyRJezL/e7hCr0VERERERHTZq0gwJ0nSswCmK/G1iIiIiIiIaGHV7pn7uCAIxzNlmI3FniAIwocFQTgkCMKhycnJKh8OERERERHR2lDNYO5fAWwEsAfAOIAvF3uSJEn/LknSPkmS9rW0lLUbj4iIiIiI6LJXtWBOkiSvJEkpSZLSAL4F4EC1XouIiIiIiOhyo6vWFxYEwS1J0njmr78B4OR8zyciouoYnY3ixfNTsBp1sBp1sBm16Gy0oNVhWu1DIyIiomWoSDAnCMKPANwEoFkQhBEAnwdwkyAIewBIAAYA/M9KvBYRES3OFx7pw0PHxnIesxl1OPb5W6DVCKt0VERERLRcFQnmJEl6b5GH76nE1yYiouW5MBHCtRua8Pl37EBYTOJnR8bw/ZcHEYol4bToV/vwiIiIaImqVmZJRESrT5IkDPrC+K19ndjW7gAAXJwMAwACsQSDOSIiojpW7dUERES0iiZDIsLxFHqbrepjDrMcwPmjidU6LCIiIqoABnNERGvYwFQEANCTHcyZ5GAuEGMwR0REVM8YzBERrWEDU3JJZU+TRX3MYZYr7APR5KocExEREVUGgzkiojVswBeGTiNgXYNZfYyZOSIiorWBwRwR0Ro24Aujy2WBTjt3uld65gLsmSMiIqprDOaIiNawS1MRdGeVWAKA3aiDIACBGMssiYiI6hmDOSKiNUpZS5A9/AQANBoBNqOOmTkiIqI6x2COiGpCLJHiqPwKmwyKiOStJVA4THr2zBEREdU5BnNEVBP++hen8b57Xlntw1hTLmUmWXY3FQnmzHpOsyQiIqpzDOaIqCac84Zw1hOEJEmrfSh1ZyYch5hMFTw+4JODud5iwZxJx8wcERFRnWMwR0Q1wROIQUymMRNhgLEYiVQab/3Ks/jiI2cLPjbgi0CvFdDRYCr4mJyZ4/eaiIionjGYI6JVJ0kSPIEYAGBsNrrKR1NfDg3MYCIo4qk+b8HHBqbC6GzMXUugcJj0CHKaJRERUV1jMEdEq24mkkA8mQYAePyxVT6a+vLEGTmIG/BFMDITyfnYpanCSZYKh5nTLImIiOodgzkiWnXZAdy4n5m5ckmShCfOeNU9ci+e9+V8bNAXQU+Rfjkgk5kTk0il2aNIRERUrxjMEdGq8wbmgrkxZubKdmEyhEFfBL/3hg1othnxwoUp9WMTQRHRRAo9zZain+sw6wEAIZZaEhER1S0Gc0S06pR+OYNOg3H2zJXt8dMTAICbt7fi+k1NeOG8T50GqqwlKJ2Z0wEAJ1oSERHVMQZzRLTqPP4YBAHY4XZgnJm5sj1xxotd6xxwO824fmMzpkIi+r0hAMCgspagZM+cnJnjonYiIqL6xWCOiFadNxBDs82ITpeFwVyZpkIiXh+awc3b2wAA129uBgC8cF4utbw0Ja8lcDsL1xIAcs8cwMwcERFRPWMwR0SrbtwfQ7vDhA6nCR5/DGkO5VjQU30TkCSowdy6BjN6mixqMDcwFUanq/haAkCeZgkAgSh75oiIiOoVgzkiWrZIPInP/uwkZiPxJX2+NxBDm8MEt9OEeCoNX3hpX+dy8sRpL9xOE3Z2ONTHrt/UjFcuTSORSmPAF0ZviX45gJk5IiKitYDBHBEt27FhP77/8iAeP124uLocnkAM7U4j3A1mAFxPsJBYIoXnzk3h5u1tEARBffz6Tc0IiUkcG57FgC+M7vmCuUzPHHfNERER1S8Gc0S0bGFRLtU7PR5Y9OfGEinMRhJoz2TmALBvbgEvXfAhmkjh5h1tOY9fu6EJggD89MgoYok0ekusJQAAu1EHQQACXE1ARERUtxjMEdGyheNyQHBqbPHBnLJjrt1phtuZycxxPcG8Hj/jhdWgxTUbXDmPN1oN2NnhwINHRgEAPSUmWQKARiPAZtQxM0dERFTHGMwR0bKFxRQA4MxYQN1zVi5PJgvX7jChyWqAQathZm4ekiThyTNe3Li1BUadtuDj129sRjgu/3uU2jGncJj07JkjIiKqYwzmiGjZlDLLoJjEyMzismoeNTNnhEYjoN1pwhiDuZI8gRi8ARHXbmgq+vHrN8krCgxaDToyPYilOMx6TrMkIiKqYwzmiGjZQuJcQHBqzL+oz1Uyc20OuV/O7TSxzHIe3oAIACUDtf09Lhi0GnS6zNBqhKLPUThMOmbmiIiI6hiDOSJatkg8CYNWA40AnF5k35wnEIPVoIU9Myrf7TSxzHIeSo+hEvzmMxu0uHVXu5qhm4+cmWMwR0REVK90q30ARFT/QmIKDrMeDRb9oidaegMxtDnnAhN3gxnewDhSaWnBzNLlaCITzLU6jCWf87X3XlXW13KY9AhymiUREVHdYmaOiJYtLCZhM2qxs8Ox6ImWHn8M7VlZpg6nCcm0hKmQWOnDXBM8gRi0GgFN1tLBXLmcZWbmYokU7vjGC3ju3OSyX5OIiIgqh8EcES1bJJ6ExaDDDrcD4/4YpsPxsj/XGxDRnp2Zy6wnGGPfXFHegIhWu7EiWUuHWYegmEQqPf8E0tcHZ3B0eBbPn59a9msSERFR5TCYI6JlC4lJ2Iw67OhwAADOlFlqmU5L8AZyM3PuBi4On483EENriX65xXJk+hRDC5RavnJpGgAwMs0Am4iIqJYwmCOiZQuLKViNWuxwy8FcuRMtp8IikmmpaGaOwVxx3kAMbfbll1gC8gAUAAtOtHzlkg8AMDQdqcjrEhERUWUwmCOiZQvHk7AYdWiyGdHuMJU90dLrl/visiczNlr0MOo0XE9QQn5Z6nI4TPIMLP88fXNiMoUjQ7MAgOEZBnNERES1hMEcES1bWEzCZpADg50djrInWqoLw7OCOUEQ0NFgZmauiFgiBX80UXItwWKVk5k7NuyHmExjX3cjZiMJBLmXjoiIqGYwmCOiZZPLLOVgbkeHAxcmw4glUgt+nhrM5WWa3E4TxvzMzOVTdsy1VqrMMtMzF4iW7pl75aJcYvmuq9cDAIbZN0dERFQzGMwR0bJIkoRwXF5NAAA73A6k0hLOeoLqc1JpCZ/92Uk8ecab87levzxmv9mWG5y4nWZ4mJkr4A0UlqUuh8MsB+DzZeZeHZjGtnY7dq9zAmDfHBERUS1hMEdEyxJNpCBJgMWolFnKF/3ZpZb3PH8R3395EF998lzO53oCsaJj9t1OE7yBGJKpdJWPvr54S2Qyl0otsyzRM5dIpXF4cAYHe13oclkAACPsmyMiIqoZDOaIaFlColyip5RZrm80w27UqUNQ+jwB/MNj/Wiw6HF8xI8LkyH1cz3+WNEsk7vBhLQETAS5ODybEsy12SsTzNkMOggCECixmuDEqB+ReAoHepvgtOhhN+kwzMwcERFRzWAwR0TLEhbl3jilzFKjEbDd7cCpMT/iyTT+5N5jcJh1+OHvXQONADx4ZFT9XE/ejjlFh7qegP1Z2byBGIw6jVoeuVwajQC7UVcyM/fKRXm/3IFeFwCgs9GC4Zni/yaTQRG+EINvIiKilcRgjoiWJZzJzFkMcwHGjg4H+jxB/OPj/TgzHsDfvesK7Ohw4LqNzfjZ0TFIkgRA7pkrVjKoLA4fm2XfXDZlLYEgCAs/uUwOs75kz9wrl3zY2GJFS2bgSqfLXLJn7qM/OIxP3nesYsdFREREC2MwR0TLogRzNmNuMBeJp/DNZy7gPfvW4y072gAA79zTgaHpCI4MzyIsJhEUk8XLLB1yZo5DUHLJC8MrU2KpcJj0RadZptISDg3M4OCGJvWxLpcFIzMRNRhXJFJpHB/x4+RoeSspiIiIqDIYzBHRsoTjuT1zgDzREgDWNZjx2dt3qI/fuqsdRp0GPzsymrWWoHDMvsOsg8Wg5XqCPN5ADK2OyqwlUDjMuqKZudNjAYTEJA5mSiwBoNNlQSyRxmReOeU5bwjxVBpTIRGzkfiSjiMaT+FT9x3DFEs1iYiIysZgjoiWRemZsxq06mPb2u149971+Oe7roI9s8sMAOwmPW7e0YZfHB/HSKb3qj2ThcsmCALcThPGWWapkiQJ3oBYsbUECjkzVxjMvXJJ3i93sHcuM9fZKE+0zN81d3LMr/73+YkQluLUmB/3Hx7Bq5eml/T5REREl6OKBHOCIHxbEIQJQRBOZj3mEgThcUEQzmX+bKzEaxFRbQmLhZk5nVaDf/itK3FVV+Gv/R171mE6HMd9h4YBlB6z39Fg5gCULEExiWgiVXRgzHI4zHoEi0yzfPniNLqbLDn/Pp0uOfDOn2h5atQPZbvEuSUGc8pUVOXniYiIiBZWqczcdwDcmvfY/wbwpCRJmwE8mfk7Ea0x+asJFnLjlhY0WPR4+MQ4AJQMTtxOE8bYM6fyZr4XFS+zLJKZS6clvDYwnVNiCQDr1cxcbjB3YtSPq7saYdJrlpyZi8TlDC+DOSIiovJVJJiTJOlZAPm1Me8E8N3Mf38XwB2VeC0iqi3KRXh2meV8DDoNfn23G2kJcJh0MJf4vHaHCVMhkYvDM7wBuZes4mWWZh2CYhKp9NxQk7PeIPzRRE6JJQCY9Fq02o0YzlocnkpLOD0ewO71TmxssS0/M5f5eSIiIqKFVbNnrk2SpHEAyPzZWuxJgiB8WBCEQ4IgHJqcnKzi4RBRNYTFJIw6DXTa8k8nd1y1DkDpEksAaLIZIUnATKT42Py1Yiok4hfHx5BYIGhVFoZXvMwy09MYyiq1PDQ4AwDY3+MqeH6ny5KznuDiZAixRBq7OpzY3GrDeW9wSccRZpklERHRoq36ABRJkv5dkqR9kiTta2lpWe3DIaJFConJnLUE5djb1YgulwVdLkvJ57isBgDAdHhp0xHrxbeevYiP//AI3v7153FooPTwD2X6Z+WnWcrBXPZEy9cHZ9BiN6o9ctk6G805A1CU4Se71jmxuc2OMX9MzbItBoM5IiKixatmMOcVBMENAJk/J6r4WkS0SiLxFCzG8kosFRqNgB/cfRB/c8fuks9pssnBnC+8tkfV93mCaHMYEYwl8e5vvoQ/ve9Y0QB2IhCD3aTLWc5eCQ6T/PX8WX1zhwansbersehy8i6XBeP+qJpJPDkagFGnwcYWKza22AAAF5ZQaqmUV4ZEllkSERGVq5rB3M8BfCDz3x8A8GAVX4uIVklITMK6hACjK29SYr4mq5yB8oXWdmau3xvEdRub8fj/eiM+cuNGPHBkFLf80zPw5e1bq8ZaAiArM5cJ5iYCMQxPR7Gvp/gA4vUuC9IS1LURJ0f92O52QKfVYHObHMwtpW9OychF4szMERERlatSqwl+BOAlAFsFQRgRBOFuAF8A8BZBEM4BeEvm70S0xoSXUGZZjsuhzDIQS2DcH8PmNhssBh3+923b8P27D2IqFMfTZ3N7iL3BWMX75YC5njmlzPL1Iblf7uru4sGcsmtuaDqCdFrC6bEAdq2Tl8R3uyzQa4UlTbRUSjOXUqJJRER0uarIFZgkSe8t8aE3V+LrE1HtCotJOC2Gin/dRoscZPjWcDB3LjMsZGubXX3smg0utNiNeKZ/Eu/eu1593OuP4ZqNTQVfY7kcZvltIBCVg6hDAzMw6DTY1eEs+nx119xMBEPTZgTFpPpcnVaD3mYrzk8sfgjKXGaOZZZERETlqvztdCK6rITjKaxrXFzPXDl0Wg0aLHpMr+GeuX6vnMHakhXMCYKAN25uwZN9XqTSErQaAem0hIlglcssM5m5w0MzuHK9EwZd8cINt9MMnUbA8HQE9ky/3a51c4Hf5la7OhRlMcIi98wREREt1qpPsySi2vfi+Sn82f3Hi34svMSeuXK4rIY1XWZ51hOExaDFuobcqZE3bW3BbCSBYyOzAIDpSBzJtFSVMkubQQdBkHvmYokUTo76sbe7cCWBQqsRsK7RjOGZKE6OBqDXCmqvHABsarVheDqCWGJxGbZwnGWWREREi8VgjogW9MvTXtx7aLjocIqQmIS1Cj1zANBsNa7pASj93iA2t9qg0eROjbxhUzM0AvBMpm/O45eHjbRVeC0BIE8WtRt1CMSSODHqRyIlYW+JfjlFZ6O8a+7UmB9b2uww6uYys5vbbEhLwMXJ8KKOg2WWRESXt1+dncDJ0cVXdlzuGMwR0YKUvrX8wEqSJETiKVgXuZqgXC6rYU33zPV7QzkllopGqwFXdjbgmX45mJsIKjvmKp+ZA+RSy0A0gcOZZeFXdzXM+/xOlxkj0xGcHPUX9NZtalUmWi6ub04ps2Rmjojo8vS5B0/hX5+5sNqHUXcYzBHRgpQx+fklj2IyjVRaqlpmzmVbu2WW0+E4pkJi0WAOAG7c0oJjI7OYCcfhDcjf/2r0zAHyRMtALIFDAzPY0GxFk23+DGCnywJfOI6ZSEKdZKnobbZCIyx+15wSxMWTaXWHHRERXT6CsQT7ppeAwRwRLUjJyOUHVsoFeDVWEwBAk9WAmUgcqbRUla+/mvozkyy3tBcP5m7a2gpJAp49N6mWWbbaK19mCcgTLf3RBF4fmim5kiCbsp4AAHauy83MGXVa9DRZF7VrTpIkhMUkTHr5LSnCxeFERJcVSZIQEpM8/y8BgzkiWpAvM1FyKm+RtXLStVRxAIokAbOR2sjOJVNpfOeFS4se7lGMGsxlDQ/JtnudE40WPZ7pn8REMIZmmwF6bXVO2Q6THmfGg5gOx7GvnGDOJQdzGgHY3u4o+PjGVtuigrl4Ko1kWkKrXc48hrg4nIp45aIPH/n+4TV5c4focicm00ikJEQSPP8vFoM5IppXKi2pGbnSmbnq9cwVe93V8urANP7iodP41dmJZX+tfm8QdpOu5IRKrUbAGza34Nn+KYz7Y2qgUw0Os179t1xo+AkAdGWCuU2tNpgNhf/2m1ttGJgKl10uqfTLKZlHltlQMb/qn8Sjpzw1cz4gospRh2AxM7doDOaIaF6zkTiUG+H5F1HKOPmqTbPM9G7VyhCU0Zmo/OdsbNlfq98jDz8RBKHkc27a2oKpkIhXL02j3VnFYM4k75pzmvXY2FI8U5it0aKHw6TD7nXFB6VsbrMhmZYw6CtvoqXyJq70BDKYo2K8Afn3bqZGMvVEVDmhOp5ovNp93gzmiGhe2YHUVN40S+Wiu5pllkDtZOZGZ+Vgbjzz51JJkoT+iWDJ4SeKN2xuASC/uVVjLYHCYZb//a7uaihYk1CMIAj4zw/ux6feuqXoxze1yP+/znnLK7VU3sRb1Mxc/b2ZU/UpwVytnA+IqHKCMfl9IFxnZfZ9ngB2fO5RDEwtbh1PJTGYI6J5ZffJTYdze+aUi+5qDkAB5qZprjYlMzfuX15mbjIoYjaSKNkvp2ixG9VpkVUts8xk5vb1lF4Wnm9vtwtup7noxza2WgEA58vsm1P2FyqZOa4noGKUqa4M5ojWHuW8H42nIEn10xd7YSKMRErC2DJv8i4HgzkimpcyybLDaSossxSVMsvq9Mw1KsFcjVy8KZm5Mf/yTtr9mYzV1gUycwBw05ZWANVbSwDI5ZUAcHXXwv1y5bAYdFjfaC57CEooc1NAyT4WW05PxMwc0doVymTmkmkJ8TpaT+OPJgAAseTqVZQwmCOieSlZsc1t9oIyS+VOmrVKZZZ6rQYOk65mLt7G1DLL5WXmzmYmWW4uI5i7eUcbAGBDi3VZrzmfN29vxadv3Yr9PZUJ5gB5OEq5wZxyU0DJPrJnjvJF4km1DKtWzgdEVDnZ5ZXROuqbm43K56NYYvUCUAZzRDSvqVAcGgHY2GIruIiKVHkACgA02Yw1kZlLpyWMzcagEYCJYAzJZdw5POcNwmU1oNlmWPC5ezob8Nyn34SDveWXQC5Wg8WAj920CboKrj7oclnKLjtRbgq0ZjJzIfbMUZ6JQHa59+qfD4iospSbNQAQrqNgzh+RM3MiM3NEVKt8YREuqxHNdgOiiVTOHbOQmIJBq4FBV71TSZPVgOnQ6l+8TYVExFNpbGt3IC0B3uDS+/jOeoPY3Gqbd5Jltk6Xpezn1opWuxH+aKKsnXxKJq7JaoBGYJklFVJKLAEGc0RrUXavdLSO3gPUMktm5oioVk2F4mi2GeaGkWQNQQmLSViq1C+ncFkNNXHxNpLJMimliEudaClJEs55Q9javnCJZT1rzfT4TZYR9CqjqG0mHawGHQegUAHl5ondqONqAqI1KJSdmauj6oxZJTNXxo3LamEwR0Tz8oVENNkMaLJmdr5lZcnC8WTV+uUUTTZDTgC5WpSSQWXi49gSJ1qO+WMIicmy+uXqmbIAPDujUkpITEKvFWDUaWExatkzRwUmMj9H29z2nHMQEa0N2Tfx6mnXnNozl2RmjohqlC8cR5PVCJetcOdbWExWbS2BwmU1YCaSQDq9uqOKlbUE+zKZuaWOIe7PDD8pZ5JlPVOGmUyUkZkLi0l1V6HVqKurfglaGd5ADCa9Bt1NVmbmiNag3GCufm7o+aPysZbTUlAtDOaIaF6+UDyTmStcExAWUytQZmlEKi2pdemrZXQ2CodJB7fTDLtRt+Qyy36PHMwttGOu3ilrBibKzMwpNwVsRh0zc1TAGxDR5jChyWqALxyvqz1URLSwUCwJrUbuDa+nzJw/c3NJZGaOiGpRLJFCSEyi2WZEk00ps8zqmYtXPzNXLIhcDWOzUXQ0yEuy3Q2mJZVZptISnuqbQJvDiAbLwpMs61mjxQCdRihrUExYTKq7Ci0GlllSIW8ghja7CS6rAfFkuq4u9ohoYSExiRZb/e0anRuAwswcEdUgJYBqthlgNWhh0GkKyixXomcOKJxgN+gLY89f/RKHB6er+vqKkZko1jdmgjmnGeNLWBz+1SfP4ZVL0/ijN2+u9OHVHI1GQKvdmDNSvpRIPKWut5Azc7xQp1wTQRGtDiMarcXPB0RU34JiUl1PUy83a+LJtNoWwMwcEdUkJQvXZDVCEAS1xEmxMmWWysVbblDw6qVpzEYS+Oenzlf19RWjs1Gsy2TmOhpMi14c/uQZL7725Dm8e+963HWgqxqHWHNaHCZMBBdXZin3zNXPXVmqPkmS5MxcpswSYDBHtNaExaQ6OKtegrns9g9m5oioJilT45TsWJPNsApllvLJfSpvgt2Zcbn37OmzkziXGSpSLYFYAsFYEuuyMnO+cLzsk/fAVBifuPcodq1z4G/u2FV3O+OWqtzMnDwARSmzZGaOcoXEJCLxFNqYmSNas0KxJBotBmg1Qt2UWWYHcyL3zBFRLZrKBG7NmTp2l9VYWGZZ5WCu0aoHUHjx1ucJoLfZCpNeg/947lJVj0GZXKn2zDnlSY2eMvrmIvEkPvKDw9BqBPzr/9gLk766mcxa0uYwlpWZC4vZZZbsmaNc3swNAWbmiNaukJiE3aTP9E3Xxw09f3TuPCQmmZkjoho0lZeZa84qsxSTKSRSEqyG6gYnRp0WdqMu5+JNkiScGQ/gYK8L7967Hg8cGS0raFgqZS3BXJml/OdYGX1zf/PfZ3DWG8TX7rwKnS5L1Y6xFrXaTZiJJBZ8k8vO8FoMOkQTKaRWeRUF1Q5lImqr3cTMHNEalE5LmXJ7LSwGLaJ1Vmap1wqIMTNHRLXIFxJh1mvVHWAuq0EtvVTunFU7MwcALltur95EUMRMJIFt7XbcfcMGJNJpfP+lwaq9/mgmMzdXZiln5hbqm5MkCY+cGMcde9bhjVtaqnZ8tUpZTzC5wETL7AyvEtTVS5kNVZ83c6OmzWGE3aiDXitgmrvmiNaMSKZlwWbSwWqon77p2YgczLXaTczMEVFt8oXjalYOkIOqaCKFaDyllsKtRDDXZDXkDEA5PR4AAGx3O9DbbMVbtrfh+y8PVu1u3uhMFAadBs2Z/j23Uw7qFppoOTwdxUwkoS4av9woi8O98/TNKRne7AEoAOqmzIaqT/n5aXWYIAgCGi0GTIcYzBGtFaGYfD1hM+phrqPMnBrMOYzMzBFRbZoKiep+OQBqMOMLi+qds2oPQAHkXj1f1sVbX2b4ybZ2BwDg99+4AbORBO4/PFyV1x+djaLDaYIms9DUbNCi0aJfcNfc0ZFZAMCezoaqHFeta7ErmbnS3yclaFMGoCj75kLsm6MMbyAGm1GnnmtcVgMzc0RrSEiUg6J6y8wpZZatdiOnWRLVqsOD0zg9Fljtw1g1vlAczdaszJyywDsUVzNzlir3zAFKZm7u4u3MeADrGsxwWuThKPu6G7GnswH3PH+pKr1Wo7NRtcRS4XaaMT47f2bu6NAsTHoNtrTZK35M9aDNIWfmJuYps8zP8Cp7C1lmufZF4knEy9jNNBEQ1f1TQCaYY88c0ZoRzGTm7EZdXWXm/NEE7CYdLAYd98wR1ao/ufcYPnnfsdU+jFXjC4sFZZaAPHxAyaisSGbOJl+8SZIcqPV5AtjWPhcgCYKAD79xAwZ8EXz0B4fxk8MjFR2IMjozt2NO0dFgxvgCmbljI7PY1eGEXnt5nmqbrPKYaW9gnsxcXoZXCeqYmVv73vUvL+JLj/Ut+DxvIIa2TMkuADRaDZhhMEe0ZoSybupZjVp1EXeWFQczAAAgAElEQVSt80cTaLDoYdJrVjUzV/2rMKI6FYwlMDQdAQAM+sLobrKu8hGtLEmS4AvFS5RZxtUplivVM5dMSwhEkzDqNbgwGcYtO9pznvPWne344PU9eOjYOH552gsA2OF24P/cvh3XbWxe8muLyRQmgqI6wVLR0WDCq5d8JT8vkUrj5Kgf77ume8mvXe80GgEttvl3zRVk5jJlluyZW9tSaQnnJkIwl5HZ9wZj2Ns113faZM0diERE9U15H7AZdTDrdXWTmZuNxOE062HUaZmZI6pFZz1zi6gfOekp63Mi8SQeODKC8xNBNYtUrwLRJJJpSd3rBMxl5nwhce5OmmEleuYyrxsWcX4ihFRawna3I+c5Wo2Az799J1798zfjv//oBnz61q0Y8IXxsyOjJb+uLyTi8ODMvFOolF1y+Zk5t9OMQCxZcifaWU8QYjKNKy/TfjlFq8MI7zxlliFlKmrezQGWWa5tvpCIVFpCv2f+c6UkSfAGRLVkFwAaLQb4owkkU6t38URElaOWWZqUzFztnP8ngjEcHpwp+rHZaAINZgOMzMwR1aYzmWCuw2nCIyfG8ZEbNy74Od/81QV87anzAOThD9dsaMIbNjXjXVevg67OSu2mMtMjlSEWgHzBbdBpMB2Oq8uvlUxKNSnZwelwHJemwgCAbe7ifWgajYCdHU7s7HDiF8fG1V15xfztw334yesjMOk12NvdiGs3NOGtO9uxOavHTd0x11iYmQPkiZabWguP5ejw5T38RNFqN2JkpnRvYX5mzsYyy8uCJ1N6G46nMDobxfrG4jsY/dEE4sk0WrOCOaX0eyaSyDk/EVF9CmVn5gxaRGooM/cvT1/Ajw8N49RfvhWCIOR8zB9NoKPBDFMmMydJUsFzVkJ9XV0SraAz4wE4TDq879oeHBvxY2QmMu/zk6k07j00jGs3NOGLv7kb129swisXffj0T47j2XOTK3TUlaNMj2yyzl0sCYKgLg4PrfBqAkAu7+zzBGHSa9BTRtlrs92IqVDprNDITAQbmq2460A3fKE4/uGX/fj1rz+f8289khlysr4h92JTWU8wVmLX3LHhWbisBqzPCwIvN60O07x75rLfxIG5gTqlMp5Um8b9UfzpfcfK/nfzZPWbZldBFDwvE/S152XmAGCGEy2J1gRlNYHVKE+zjCfTNZN5vzgVRiSewkxmDUE2fyQhl1nq5XBqtUotGcwRldA3HsA2twNv2y33Zj26QKnlU30T8AZEfPD6Hvz2/i585c6r8NSnbgKAupyI6csEQdkDUIDMAu+QiLCYhFYjwKir/mlEKbOcDsdxZjyArW12aDUL3/1qthkwNU8gMRUSsc1tx+fevgOPfuKN+NWnbgIk4BtPX1CfMzoThSAA7U5Tzueqi8NL7Jo7OjyLPZ0Nq3KXrpa02o3wheMlpxZGSkyzZM9cffnv4+O47/AInj8/Vdbzs0tvz3pLB3PKjrm2rGmWTVlTdYmo/oXEJAw6DQw6jXpDL7KKZYvZhnxyNVD+e70kSfIAFLMeJp18zOIq7ZpjMEdURDotoc8TxPZ2O7qbrNjhdizYN/fDV4fQ5jDi17a1qo/ZjDp0N1lwZrz0xUqtmsoMGCgI5qzGzDTLJKwG7YoEK3MrEUScGQ8U9MuV0mIzYioUL9mTMxkU0ZI14KWn2Yo7D3TivkPDGM4MvxmbjaLVboQhL2htd5ogCMBokcxcMJbA+ckQrlx/eZdYAnPrCUplSJWpZUq5rkYjwGLQMjNXZ46P+AGgZG9JPq8/Bq1GQLvDhP55MnPKJNScnjkrM3NEa0lITMKuVmdk+qZr4IZeKi2pbQKevOnV4XgKybSUl5lbnWNmMEdUxPBMBJF4Sg0a3ra7HYcHZwp+mRUjMxE80z+J397XWdAbt63djjPj9ZeZUzJaLktuMKeUWYbjqRVZSwAAJr0WVoMWZzxBzEQSOWsJ5tNsMyKeSiMQKwwMYokUArFkQc/NR2/aCI0g4BtPy72Po7OFawkAQK/VoMVmLLpr7sSoH5IEXNnpLOs417LWzPe31HqCkJiETiPAkPV7Y6mjpbEkOzEqB3OHBqbLer4nEEOLzYjtbjvOekMlnzeR+bnJ/j1tysrUE1H9C4lJ2Ey5pfa1MARr3B9FMrO7dizv+m82czOpwTKXmYsxM0dUO5TgSwnmbtvtBgA8enK86PPvfW0YAPDbB7oKPrbd7cAlX7gmTkyL4QuLaLToC4JTZWFvWEzCskLBHCCXd76YKeEqNzPXbJcv+oplhZTHmm25wZzbacZ7D3Ti/sMjGJ6OZBaGFx/O4C6xa+7YsHxhe7kPPwGAVvv8i8PDYhJWoy4nw2szallmWUf80QQuTYVhMWhxcjRQ1lQ3byCGNqcJW9rtuDARQqJEf4w3IGb2OM0NWmqwMJgjWktCsWRB33QtDEEZ8s31z3vyyiz9UbmHzpmZZgkAMWbmiKpryBfB9V94Cv3z9GcoTo8HoRGALZmphhtbbNjSZitaaplMpXHva8O4aUtL0QzOdrcDkgT0z3P3uRbl75hTuGwGROIpTIXEFRl+omiyGtUG5G3lBnOZ4y/WN6dMuSw2De9jb9oEjUbA1586h/HZmDq5Ml+H04SxIj1zR4dn0NNkUS86L2dKr1PpYK4ww2s16lhmWUdOZEos3713PeKptJqlm4+8CNyIrW12xFNpDGb6Uoo/L/f3z6DTwG7SMZgjWiOCYnYwp6ynqYFgLtNuYdBqCm7c+iNKMMeeOaIV8/IlH0Zno/jJ4ZEFn9s3HkBPszVnoe1tu9x4dWC6YDLfk30TmAiKuOtg8eXQ29vlwKPeSi19oXjOjjmFsjh8eDoK2wqsJVAox7KuwQynWV/W56jBXJFBCcq/Y7Fgrs1hwl0HunDf4RHEU2msLxKkA3IWb3w2VtCTd2zYf9nvl1M02YzQCHPlcvnkzFzuz5HVoONqgjpyfFRew/G71/UAAA4NLNw35/HH0O40qTfMznqK3+zyBkW0OorcVMpUCBBR/QtnB3PG2imzHJqOQKcRsL3DUdBmM5vJzGVXDjAzR1RlSpP9Iyc9Cy70PuMJqEGY4m273ZAk4LFTudm5H706hHaHCW/a2lL0a61vNMNm1NVdMDcVFgtKEIG5YSSeQEy9g7YSlNfdXmK/XDFzwVxhVmi+YA6Qe+eUPq78HXOKjgYToomUWm4ByBepnkCMJZYZWo2AZpsRE4FSA1CSBRleq7G29gzR/I4P+9HdZMGGFht6m604PDh/31w0LvertjlM2NRqg0YoPdFyIhDLGX6icFkNHIBCtEYU75lb/feAoekI1jeasb7BXBDMzZVZZg1AYWaOqLqUi4Wh6QhOzxNYBWMJDE9HC4KGLW02bGix4htPn8dnf3YS/++VQTxx2otn+ifxnv2Fg08UGo1Ql0NQ5DLLwsycK+uxlRqAkv2629rLK7EE5As+jTB/MJe9Ry9bm8OEuw7KPZCdpXrmiuyaU5aFMzM3p9VhhDdYegCK1cAyy3p2YtSPKzKTW/d2N+Lw4My8N8y8WbvjTHotepqtOOspPD+m0xImgmLOWgKFy2LgagKiNSK7Z85aY2WWXU1WuDMtFdnntdlIVmZOHYDCzBxRVZ31BPGmrS3QCPPvjFMW2OYHDYIg4HO370BHgxkPHBnFZx44id/73iEIAH57f+e8r73d7UDfeHDBjGCtiCfT8EcTRQOd5qzH8svjqqlJzcyVH8xpNQJcVkPxYC4UQ4NFX7ByINsnb9mKr965B5tabUU/7m4o3DV3bGQWOo2AHYs4zrWu1W4qnZljmWVdmwqJGJ2N4op18uTW/T2NmIkkcGGyeA8cMLcIXMm4bW2zF+0p9oXjSKUlZuaI1rhgVmbOXEPTLIemI+hymdHuNCGWSOdU4fijCRi0Gpj12lUfgLJyt9WJVtFMOI6JoIjf29iEaCKFR0568MlbthZ97plMMLe9o/Bi/KatrbhpayskSd49ctYThFGvKTr4JNs2tx3Bl5MYmYmi01U8y1NLlIukhTJz+RmVaupyWaARgN3rFjfuv9lmxGSw8KJvKhjP2TFXjM2owzv3rCv5caWX7s9+chw7O5zY1m7Hr85OYrvbkTN973LX5jCqe8jyhcVUkTJLXU3claWFHR+RM9FXrJd/L/d2uwAAhwenS94EUTNzTvn3b0ubHY+e8iCWSOX83ijPa7UXD+Z8YXmH5ErsuiSi6ogn04gn0+qeuVrJzPmjCcxGEuhyWdQqnHF/TB1s5o/G4TDrIQjC2h+AIgjCgCAIJwRBOCoIwqFqvx5RMcoEyy1tdty2y43zEyGcnyjeo3FmPACHSYcOZ/EJhoCcpet0WXDzjja8YXPxXrlsSjapWKllLWbrSo3tBwCrQatms1ZymuUtO9rx1CdvQlfT4oLhZpuxRGZOLNkvV65WhwlfeNduvHFzCyaCIv7zhQGc9QZxsNe1rK+71rTYTfCFRSSLjJ8Px5MF5bo2oxbheLIqvxsTJco9V9LYbLQmf++X4viIH4IA7MzcZNnYYkWjRT/vEBSl90TJuG1rt0OSgPMTudk55d+qaJml1YB4Mr3qF3xEtDxKSb1yPWHSayAIQGSVqzOGM5Msu1xWtGeuB7P75mYjCTRY5GFspsskM/cmSZKmVui1iAoo/XLb2h3Y1u7A539+Co+c8OAP31w4TKNvPIBtbkdF7/ZubbNDEIA+TxC37GxXH0+k0njbV5/DTCSObe0ObG23Y1u7HTdtbV12oLEcSi9Kc5HMnCAIaLYaMOaPrWgwp9EI6Gm2Lvrzmm0GDBQZez4ZFCsypOTOA124M7NfMJFKY2g6smCm9nLT5jBCkuSpou15N0nCYrJgkI7FqIMkAdFEqqJDdn5+bAx/9KMj+PvfvALvWaA0ulrOeYO45SvP4oe/dw2u3di0KsdQScdH/NjUYlMDckEQ1L65UrwBEVaDFnaTfCG0pV0+D/d5gtiVlXn3Zkpzi5VZNmYtDl/J8xARVZZSUp99DrHoV38I1pAazFnQaJXPVdmriPzRBBoyk7WNaz0zR5enVy768HcPn1ntw1Cd9QThMOnQ5jCi3WnC1V0NeLhI31w6LaHPE6x4v5PVqEO3y1KQmXuqbwLnJkLYvc4JfzSBH7w8iD+9/zg+88CJkl9rKiTipQu+ih5fPl84MxykRBmiUmppNdR+KaGSmcvPhEwGl5+Zy6fXarCxxcYSyzxKmZw3bz2BmEwhkZIKVlwoF+eV7pu7P7OW5M8fOIEXL6zO/cUjw7OQpNrIEC6XJEk4PjI3/ESxt9uFi1Nh+IpkxIHM7risAK3bZYFBpynYAar8vBT7PW2ycnE40VoQjMnnebtp7qaM2aBDuFaCuSYLWjIrdvIzc8qapNXumVuJYE4C8EtBEA4LgvDh/A8KgvBhQRAOCYJwaHJycgUOh1bCf702jH979mLJi7HpcBwf/t6hkm/2ldbvDWJru13Ntt22y40z44GCRbXDMxFE4ilsay9//H25trsdBcHcD1+R1xp86/378NAf3oDTf3Urfm1bKy5OlR4e8K3nLuJ37nmlqtP+poKle+YAwJUZglIPd8Rb7EbEEumcN4awmEQ0kVrV7OflpNVefHF4WJT/TQp65jI3CZSPV4IvJOKF81N43zXd6G224iPfP1xQ1rcUz5+bwn++cKns5yvngNW+61wJ4/4YpkKi2i+n2NfTCAA4VCI758kL5nRaDTa12NThU4AcKL4+NIsWuxH6IpOCGxnMEa0Jc5m5uf2xVqMW0VUegDLoi6DJaoDNqINOq0Gr3ZSzONwfTcBpUTJzmWBuDWfmrpck6WoAtwH4A0EQ3pj9QUmS/l2SpH2SJO1raVm494jqw+kx+YJlyBcp+vGXL/rwy9PeeUtxKkWSJJz1yMGc4tZdcqnjI3nZOeVCazETE8u13e3A4HREDcKGpyN49lzuWgOtRsCGZitGZiIle2oGpsJIpaWqrjqYCoswaDVqQ3K+5syF1EquJlgqdddcViCh7phbYAAKVYZy4Z6fjcrvlVAof6/kDYuHT3qQSku462AXvv27+2HQafCh77y27GDgOy9ewhcf7UMqXV4PXN+4HLCshWAuf/iJYvc6JwxaTcnzuzcQKyi33dZuz8nM3XdoBM/2T+LDb9hQ9Gu4LAzmiNaCufeBuQoNs1676pm54elIzsC6dqcpJzPnj85l5gRBgFGngbhWM3OSJI1l/pwA8ACAA9V+TVpdsUQK5yflO95Kmjqf0sM0WcHM3M+PjeED334V6byLKk8ghkAsia1tc8Fcp8uC3eucRYK5IDSCPCil0ra7HZCkuf69e18bhgDgzrzenU6XBbFEuuT3Zmhartk+NVa9YG4yIKLJZijZN6gs8K6HzFyzvXBxuPK9bWZmbkU02wwQhLkeKEU4ntsrobBVIZh76NgYNrXasK3djk6XBd96/z54AzF8+HuHlvUGfHEyjFgirTbLz0eSJPRl9qlVeh9R/nmvkiRJwo9fG8Y//vJsTtB6fMQPnUYouPll0muxa50DhwYKl4dLkoSJgFjQB7el3Y5xfwz+SAKDvjD+8qFTuHZDE+6+obfoMSml3gzmiOpbUCwss7QadYjWQJllV1Yw19FgUtcQJVJphMQkGsxz1UsmvXZt9swJgmAVBMGu/DeAWwCcrOZr0urr9wbVN/yh6eLlgoNT8oXPVJGR8Uv1H89dxDP9kziauVusUEp38gO0W3e149jwLE6N+dXjPTMeQE+zVd1zUklK6eaZ8QASqTTuPTSMN21tRUfesIxOl/z34elowdeQJEm9aDw5WnzU+3JF4yk82Tcx73AQ5ULKUhc9c/Kx5gRzzMytKJ1WgyarEZMlMnP5P0fK38MVKrMZ90fx2sA03nFlh3qD4qquRvz9u6/AocEZPHKi9N7J+SgDb4C5mzTzmQiKmMksmq3UDqVIPIlP/NcRXPuFJ6tSth6IJfDxHx7Bp39yHF976jw+dd8x9Xx5fMSPre32oj2i+3tcODkaKAhap8NxxFPpggmVys220+MB/Mm9R6HRCPjye66ERlP8hpLdqINeK2Cau+aI6looVlhmaTFoK3b+z/fSBV/RCdfZkqk0RmejOcFcu8OMcX8MkiQhEJ1bGK5Yy5m5NgDPC4JwDMCrAP5bkqRHq/yatMqUEkuNINccFzOXmavMEIBLU2F1j1X+QnCldGdrXh/c23a7AQC//rXnse2zj+DGLz2N585NVaXEEgDWN5phN+lwZjyAJ894MRkU8d7MFMRsnY3yyWNkpvB7NxNJqPXlJ6uUmXvgyCj80QQ+eH3xO+IA4M6USCl9K7VMCdgmQ3MXfcqJnD1zK6fVbizIzIUyPXGlM3OVeWP87+PjkCTg9ivcOY/ffkUHHCbdkgcKDU9HkMwENv2ehYO501ml0dH48u/gDkyF8RvfeBEPHhvDRFDEN5+5sOyvme3o8Cx+/WvP4dFTHvzZrdvwqVu24IEjo/jkj48imUrj+MhswfATxd7uRsRTaZzIu+mkLAxvL5KZA4DPPXgSrw/N4m/u2FVwoyubIAhotBgwHWIwR1TPQqIcGNmyMnMWg7YqmbnJoIj/8R8v41+env9cOTYbQyot5axCcjtNiMRTCMSSmM0Ec0qZJSBn5larZ66qNVKSJF0EcGU1X4Nqz+nxAGxGHXqbrSXLLJUgr1KZuV8cGwMA7HA78OhJD/6/27apd+D7PEG0OYzqokdFb7MVP/3YdegbD2J4JoLh6Qg8thjumGdJ9HIIgoDt7Q6cGQ9iaDoKt9OEm7YW9omuzwRzxcq2lO/n1jY7znmDEJMpdSRuJUiShG+/cAk7OxzYnxliUMyv7+5Am8NUFyP4XVa5xC+/Z04jzJWLUvW1OowFPXORFeqZe+j4OHatc2BDS+4Sa61GwIHeJrxyaWnB3KXMoCJBKC8zp/TLWQxaRBPL+//2VJ8Xf/xfR6HVCPjOBw/g50fH8N2XBvGhG3rVBbdLlU5LuOd5uRewzWHCj//ntdjbLZ8PBEHAlx47i6lQHIFYsqBfTrGvR961+NrANPb3zO1dnFDWDeT1zHU4TbAbdTg3EcI7ruzAO8s4D7usBmbmiOpcKJaEIACWrAy/xaCrSmbu6b4JpCXg9aH55zVkryVQZO+aU47NmZeZq3T5fLm4moAq7vRYANvddnQ3WYoGc9F4Sr07u1Cqu1wPHR/DgR4X3n9tN4amIzl3wPu9wZI9cFd3NeKug134s1u34Z/vuhr3f/Q6vGVHW0WOqZhtbjtOjPrx3LlJvGff3OCTbGaDFs02Y9HvnfLYbbvbkUxL6PcsfxpftufOTeH8RAgfur533j17Bp0G121sruhrV4tOq0GjxZDTgzgZFNFkM0JbooSLKq/NblIv5BX5+4UUVkPlVhMM+SI4NjyLt1/RUfTj12xwYcAXUXshFuPipBzM7e1qLBirX0yfJ4B1DWY024zLuuv809dHcPd3D6Gz0YKHPn4DbtzSgk/cvBmSJOHrT51f8tcF5DLIu7/7Gv7vw2fw5u2tePiP3qAGcgDwB2/ahD+7dRuePy+vdigVzLmsBmxqteHVS7l9c8q5P79nThAE7FzngNtpwl+/c1dZx+qyGtgzR1TnQmIKVoMup6S6Wpm5X572ApCvU+criSwWzHU0yOescX8U/kjxzJyYXIM9c3T5SWemLO5wO9DdZMHoTBTJVO4Pt/JLYtBpKjIA5awniH5vCG+/0o2bd7RBIwCPZUotU2kJ57yhnOEnq2m724F4Mg0BwG/Ps7S402Uu2jOnZOuUaZwnxyrbN/ftFy6h2WbE7Ve6F35yHWm2GQoyc83sl1tRbQ553188682u9DRL+Q5tJSY+PnRcztrffmWpYE5e3P3KxcJhHQu5OBWGy2rAgV4XLk6Gc/6/FXNmXL7RZV7mQtz/em0Ym1tt+OnHrlOnrXW6LHjvgS78+LVhDOStNkmk0nj8tHfBfo6XL/pw21efxQvnffjLd+zEN39nb86dZ8VHb9qIz96+A2/Y3DzvsKj9PS4cHpjJGZri8ccgCHPrKrJ97c6r8OAfXF/0NYtptBoww2COVtjYbBT+TJkdLV9ITBTc0LMYtBVdTQPIiYTnz09iXYMZ8VQaZ8ZL34AbnA7DoNXklIO3ZyoePP4YZqPyeafBzMwcrUFD0xGE4yns6HCgy2VBMi1hbDa3tErZ7XblemfOBfZSPXRsDFqNgNt2u9FsM2J/jwuPnvKoxyMm0wX9cqtF6ccrNvgkW2ejBcNFeuaGfBE024zY2maH3aSr6BCUC5Mh/OrsJN53TXdFSzdrgbI4XDEVqvzCcJrfNrcDaQk4lXUDQhk9nT8ARafVwKjTVKTM8qFjY9jX3ViyJHi72wG7SbekUsuLkyH0Nluxtd2OZFpSyy6LEZMpXJgMY1u7A2aDFtElvukrvWrXb2ouGDzy8Tdtgk4r4CtP9KuPTQRjuOtbL+P3v3cI35inT+TfnrmAu771MiwGHX76sevwget65s3O331DL75/98GiO+AUB3tdCIpJdYKncjxN1uK741odJrTmZezm02Q1wMdgjlbQ4cEZ/NqXf4U/f+DEah9K3QnGEnjPv72U8x4AyBUY2f1ygFxmGU2kKjql9/nzU4gl0vjEzZsBAEfnKbUcno5gvcucky1stRshCMBYZuougJz2HWbmaM1Qyht3djjR5bICkO9wZFP65fZ2uxCOp5Y11U2SJDx0fAzXbWxSMy237WpHvzeEC5MhnM1cRNROMGfHW3a04eO/tmne53W65KlJxbKaXS6zXJLU4ajoEJTvvDAAg1aDuw4WDmWpd3IwN3fRNxkUOclyhSmletm7x0JiEjqNoC5czWY16hZdZvlM/yQ+9+BJfP3Jc/jxa8O479Aw+jxBvL1EVg6Q++YO9rrw8hIzcxuarWp2ar6+ufMTIaTSErZlMnNLLSHq8wQRS6RxdVdhT2urw4Tfva4XDx4bQ58ngMODM7j9a8/j5GgAOzscuOe5i0VL2w8PzuDvHunDW3e246E/vAG71hUvnVys/b1yr1x2qaXHH0O7szK/e40WA/zRRMF5kqga+r1BfOg7ryGWSOPZs5NI8OduUY4MzeLVS9N46sxEzuPBWLJoZg7Akm96FfPEaS/sRh3euWcd2hxGHB2eLfnc/LUEAKDXatBiM8Ljj6oDUBxZQahJz8wcrRGnxwLQaQRsarWhOzMFKL/3a8AXRoNFj40tcrC3nCEoJ0b9GPRFci7W3popQXz0pAdnPSEIArC5tTaCOaNOi2+9fx+uKnIhlq3LZUEqLWHcn5vVzD7B7Opwom88UJELGX8kgfsPj+AdezrWZMYqOzMnSRImmZlbccrAnOzG84iYhNWoK5oBshoXV4p436FhfPA/X8W9rw3jy4/349M/OY4/vf84dBoBt+1un/dzD/Y24dJUGN5A+dN1g7EEJoMielus2NBihVYjqDePilFKera7HbAYll5meSTz/buqq/gUyY/cuAE2gw4f/+ER3PnvL8Fs0OKnH7sOX73zKkQTKXzj6dyeulRawud/fhJtDiO+9FtXFlxULce6BjPWNZjxWta+OU9ARJu9/OzbfJoya0eUdQ9E1TI6G8X773kVRp0Gn3nbdgTFZMlgIJlKr/qOtFqk7MbNv+kVFpM5O+YAwJI5Dy10nvzV2Qnc/I/PLBhEpdISnuzz4satLTDoNNjT2VDy30+SJAz6IujOC+YAwN0g32j3RxOwG3U5cw+MOi2DOVobTo8HsKnVBpNei3aHCQadBkN56wkGfRF0N1nVi+nlrCd46NgY9FoBb905d7Hmdpqxp7MBj570oN8bRLfLUpW9cdXUWWSiZTyZxrh/bu/JznUOiMk0LkyWLu0q172HhhBNpPDB63uW/bVqUbPdgEgmC+yPJpBISer+OVo5e7sbcXhwBpIkl86ExFTJ4MFqKD8z958vXMKf3n8c129qxpHPvQV9f30rnvv0m9+AkaQAACAASURBVHD/R67Fgx+/Hq0LBA9K39zLFwtLLV8fmoHHX3iOUkoqNzTbYNRp0dtsxdl5BhL1jQdg1GnQ02SFybD0N/3Xh2bRYjeWLBttsBjw4TduwPmJEG7Y1Iyf/8EN2O52YFOrDb+1txP/7+WhnLUn//XaEE6OBvDnb9te0UBOsb+nEa9emvs39wZiBZMsl6oxU+LkC1d+vx6RYjocx/vueQWReBLfu/sA3rO/E1qNgGf7J4s+/+8fO4u3fuVZ9WeeZErl1jlv7nkyJCbVoVcKZbLlQpVb9x0ewfmJ0II34o4Oz2IqFFcH3F3Z2YABX6Roz60/mkAwllT7kbO5HSZ4MmWW+b29Rr2GZZa0Npwa82NHpi9MoxHQ2Wgu2DU34Aujp8milkVOLjEzl05L+MXxcdy4pTVnohAgDwg5MerHSxd98zbo1yrlJJLdNzc2G0VamvvYrg65FGq5fXOptITvvjiIg70u7OyoTHlVrVF+1qaC8bmF4czMrbi93Y3wBkSMzsrDfcJisuTieatRt2DPnCRJ+PqT5/CXD53GW3e24T8+sA8Wgw4mvRadLgv29ZT3M72jwwG7UVdQajkyE8Gd//Yy/vbhMwWfowZzmQqDrW32eSda9nmC2Npuh1YjwLKMAShHhmZwdVfDvP1sH71pI374+wdxzwf251xw/PHNmwEB+MoT5wAAM+E4vvTYWRzsdeEd85SiLseB3iZMhURcmgpDTKYwHY4X7JhbqvWNckBbbFgUUSWk0hLu/u5rGJ2J4p7f3Y9t7Q44zXrs6WzAs+emCp6fSKVx/+ERDE1H1MmtJFN65S5MhnKGRYVihT1z5QzBSqbSeC4TUM8ukJ1//LQXOo2Am7a2AgD2dMqVDUdHCrNzxSZZKtqdJoz7Y5iNJgquO+XMHIM5qnNTIRHegIgdHXNLt7tcFgxmZZfEZApjs9GczFy56wn80QQuToYwG4kjnZZwaHAG4/4Y3l5k8uKtmUzddDheM/1yi+F2mqDVCDkXKfknmA0tNpj0mmVPtHzyjBejs9E1m5UDsheHi+oEVQZzKy+/by4cTxZMslRYjTp1QEop//rMBXz58X686+p1+MZdVy95cI9WI2B/r6tgCMo/Pt6PeCqNFy/4Cu6yX5gMQxCglpNvbbdjaDpS9E6yJMlTfrdlzkWWJQ5AmQ7HMeCLLFimrdPKq0M0eas3OhrM+MC13fjp6yPo9wbxpV+eRTCWxF++c+e8weFyHOiVj/W1gWl1NUWlgrneZjmQvjRV2RUtJGNPmLyX7EhmiX32vsQ3bG7G8ZHZgszO8+en1HUZpyvY016L4sk0zk+EMBkUF/xZCYtJXJoKo7fZimRawoBvrqIoKBb2zJkNSpll6Rt6R4ZnEYjJH59dYLroE2e8OLjBpQZgV6xvgCAAR4dKB3PdTdaCj7mdJoTEJEZnomjIy8yZ9JoFJwZXC4M5qpgzmRS6kpkD5F+G4emIeiE0MiNnl3qaLOoy58kyJ1q+999fxq99+Rns+avHsekzD+N997wCk16Dm7cX7oXrabaqF071mJnTaTVwO005mTk1mMtcPGo1Ara7HWod+lJ976VBuJ2mot/HtULNzIVE9eet2Gh0qq5t7XZYDFq8ngnmQkXexBU2o3bBzNzPj45hf08j/uHdVxbd2bgY12yQ1wtMZO6mnxkP4IEjo+husmAqJOL8RG7AcGkqjPWNZjWAVM4z+SVEgHwTwReOq9NsTUvcoaT0yxUbflKuj920CVaDDv/rx0fxo1eH8P5ru7Gt3bHwJy7RxhYbXFYDXr00o5ZCtToq87vXYDGg0aLHpanCyb+0PPcdGsauzz+GC5OXd6D8nRcH0O4w4Y6rcpfYv3FLCyQJ6r5FxUNHx2DPnNOW+96c7Q9/dKSg33W1feGRPtz8j89g//99Aps/8wh2/8Vj+I1/eaFoCXmfJwBJAu7YI38flSoGSZKK9sxZDQtn5n51dm6QymykdIXXpakwzk+Ecq5xbEYdtrTai/bNKdVkna7CUnZ3prz9wmSoIDNn0mshViAzJyZT+MHLg4v6HAZzVDHKXajt7tzMXEhMqneqlLUE3U1W6DPLnMvJzJ3zBnF6PID3HujC527fgT940yb81r71+Kt37ip5Z/+2Xe7M8dRfMAdk1hNkZTWHpyMwaDU5wwN2dThxeiyw5PG95yeCeP78FH7nmu5lXwzXsma73FuTHcxxz9zK02nlxvPDmaAkIqbUcpp8FsPCZZa+cBwbW2wFGailONib2TeXmbz4pcfOwm7U4Rt3XQ0AePFCbtbu4mQIG5pt6t+VCoBiEy37MsNPlKDJotchnkovenjRkaFZ6DQCdi9j2mSjVe6pOzkaQJPVgE/cvGXJX6scgiBgX3cjXh3wqWVn7RXqmQPk7Fz+Xj0Cjg3PYmx2aeWnR4Zm8JkHTkJMpvHIifEKH1n96PfK74/vu7a7YJXGlesb4DTrc/rmovEUHjvlwe1XutHTZKlYZk6SJDx1xotHM/tza0E6LeHhE+PY39OIv3rnTvzJzVtw45YWHBmaxaGBwpH/yvfiHXs6oBGAfo98TowmUkhLKJKZk98X5ts196uzk+ogvfn2/j1+Wv6+5d+w3tPZgGMjswVVF8PT8gooi6Hw2tKdOXcl0xKc5ty+e6NOg3gqnbNXczGi8RTuef4S3vj3T+P//Ozkoj537V690Yo7PR5Ah9OERuvcD7hSgqSUWg5k7qD2ZB5vsRnLysw9fMIDQQD+5ObN+NANvfjkLVvxN3fsxnv2lV68ffcbevG1916FTTUyyXKxOl1mDM/kllnm7z3Ztc6BkJjMKWVdjO+/NAiDVjPvAvO1oMma1TMXEqHXCgV31Whl7O1uxJnxIMJiUm58L5mZmz+YS6clTIfjFQvKd3Y4YDPq8PJFH16+6MNTfRP42Js2Ydc6J9Y1mPFSVjAnSZJaMqTocllg1GnUi5Rsyp415caS2SC/9S621PL1oRlsdzuWPdDpQzf04g2bm/G3v7F7RX4PDvS6MDwdxfERuSS8UmWWgFyFMd9+v0qKJ9P4p8f7ce3fPYkHj46uyGsuxZGhGbz7my/icw+eWvTnTgRj+MgPDqPNKe8zfSJvjPxaIyZT+PmxsaJlgt95cQBGnQbvPVC4rkerEXDDpmY8e25SDQae7PMiHE/h7Vd2YGeHE6fGK7MHNhBNIhxP4awnWDOlrydG/fAEYrhzfxfef20P/vjmzfjib14BnUYoyFYCcpay0aJHT5MFPc1W9aZXKFMmmf8+oAxEiSaKvwdMBGI4NRZQM33z9cw9cXoC29rtBQNN9nQ1YDaSwEDWXIewmMSTfRM57ULZss9dhWWW8nk5voQhKE+c9uKGLz6Fv/7FafQ2W/GDuw8u6vMZzFHFnB4LFPwCKP1dykTLQV8YdqMOrkzA12wvLzP38Ilx7O92LWqhrM2oq1pT/0robLRgMiiqJQvDM4V7T3YuYwhKSEziJ6+P4vYr3Gs+S2XQaeA06zEVEjEVjKPFZqxajxDN7+ruRqTSEo6NzMo9c0XufgJyA3w4nio5EW42mkAqLanj6ZdLp9VgX08jXrrowxcf7UO7w4Tfva4HAHDdxia8fMmnZsC9ARGReEq9KwzIF3eb22xFM3NnxoNwO03qglmzeqFSfjCXSks4NjxbciXBYliNOnz/7oO4Zef8Kxsq5UBm39wvjo3BmPldrJTeJis8gVjJstXh6UhFxoWfHPXjHf8/e/cdF8d95w3889u+CwtLbwIBEqgX1KzmKrc4cRwncWKn2XFydhKn59LuSe7yPHeXu5RLufhSnMSO05zmJHZ8jm25yZZcZHWBJIoKQoDosCzL9nn+mJ1hl51dFljKos/79eJlRF3wMDvf+bb79+L7zzVDJwQ+9bsj+NrjDVO6cJusYEjCa2f6kqrA6B/x4b7fHII/KGFfS++kfnZfIISP/foQhkb9+Mn7NuEta0tw9MIguocX7iCPvx3txCcfORwz5GjQ7cOfD13A29aXqdcr411ek48upxdN4dLqx450oCjLjMuq8rCyNAtt/aMJM0bJUgZG+YIhzTLuufDMiYvQ6wSuWV6ovi3DbEBdhQOvnNYO5laWZkEIgdpCu/pzDIdv2MWsJpggM/diOCO6a0UR7GZD3GBuyO3HgdZ+dYplJHUISttYJvEne06jZ9irLhYfryjLAuXSIXYAihxSTeV886M9p2Ez6/HHj2zD7+7Zhp01+ZP6fAZzlBIefxCne1xYOW5ynHInROn3OtfnxuJ8m3ohnZ9pVgdSxNPS7UJj1zBummBX1EKj/O6UMeLn+2KDudoiO4x6MaXa/L8cugCXN4APhC9YF7r8TPnGAXfMza0N5XK/16HWAYwkyMzZTAYEQ1LcUc99rtSXy26tzsOZnhEcPj+Iz15Xq95p3bYkD4NuP06Fs25nwgM3qiLKLAH571FromXk8BMAsIa/7mT65pq6hjHiC06rX26urCzJQoZJj44hT/hiKHU3UqrCAXXkQAXFsMeP6767B7fcv08ze+f2BfAffz+J9//8ddz320P48p+P4z//fgo/fLEFv36tFY8f7cCeph58++lG3PI/+9A/4sPPPrAJL37+Kty9owq/eOUc7vjpa5qrK1Lppy+fwe0PvIaf7T2T8OOCIQmf/v0R9Lp8+Nx1tRj1BzXXbcTz/55owIHWAXzzneuwsjQLu1YUQZLkISALlVIm+dC+c1HZ1t+/0QaPP4S7EgwGu6K2AADwcnMPhtx+7Gnswc1rS6HXCfXGtjJLYDoiy2UbpjnwLFV2n+jClsrcqEosANixNB/H24eietj8wRAaLw6rN59ri+041zcCjz+oZuZiloaH/x3vHLmnsQdFWWasKLEj22aM2zN3YdCNkCRXXoxXWyT3cCtDUDqHRvHAy2dw87rSuOdZk0GnPuc4NHrmAEx6PYEkSWjuGsYVNQVRQ3Ymg8EcpUTjxWGEpOjhJwDUfXNKQ+n5fnfUhKCCTDN6h30J97EoNfs3ro6dWrmQKc23bf2jGHL74fQEYoI5k0GH2iL7pE/wkiTh4VdbsW5Rtnp3aqErsJvVnrmFnomcz7JtRtQUZuK1M/3wByVkxumZU57c45VaKjeBUpWZA8b2zdUUZuLtG8YGHmxbIr9dueN8pid6LYFiWZEdXU5v1IWFLxDC6R4XlkecG21JNPePdzh8wZGKzNxsM+h12BCeZJrKEksAqMxTJlrGBmsnOpzw+OXf/1vv34vnTnap73vldC9u+N5L+MmeMxh0+3Gy04ndJ7rw4N6z+OZTjfjKX+vxyUcO484H9+P+F1rwtvVl2P2ZK3HtyiIY9Tr8880rcf976nCy04m3/OBlzWxEPJIkJb2DzO0L4KcvnYFeJ/CtpxsT9mH94PlmvNTUg39560r8wxXVsBh1SQdifzvagV+/dh73XlGtVrSsKLGjzGGNW2rZ1DWMg60DONs7ok6ZTiehkIS9Lb1467pSbK7MwZcePY5TF50IBEP45aut2FqdGzUDYLxShxVLCzOxp6kHTzV0whcM4a3r5d/dqvDnpWIISueQHMzpdcnfuD3d48KPXjw9I7vuzvaOoKnLhetXxWa7dizNhyRF7+w83eOCLxhSA6plRXaEJPlGvXJ+j+mZCwdGIxrTLAPBEF5q7sFVtYUQQsBhM8adZjkwIr89NyP2OV8f7j9WhqB86+lGhCTgCzcsS/jzK31zsQNQppaZ63F54fQEUFOYOfEHx5H6DaF0SVKWQWrd/ajIs+F8/wgCwRDa+t1RGbZ8uxmj/iBGfPGXBz9ZfxGbFuektGk+HaiLwwfcavChtcRydWk2njlxEUOjfmSaDdAnMQzi1dN9aOl24b9uW5faBz2P5Wea0dDhxLAngHWLFuY+vXSxcXEO/nxYvgueaDUBIJfZ5Gk8x/W55IAplYH56tIsvGVtCd4/biBQSbYVlXk2vHamDx++vBpne0dgMepiApPacPatqcullhae7nHBH5S0M3NxnvQlSYrJXh0+P4C8DJPm7qN0sKUyFy8396ZsYbiiMj9+MKdc+P7xI9vw1cfq8aGHD+AT1yzFoNuPX73Wiso8G/5w7zb1/5XC4w/COeqH0+PH0GgAVqNes4fmLWtLsazIjo/8+iDe97PX8YUbl+PeK6onzDze91u5DPKnH9g04c/329fPo2/Eh59+YBP+6S/H8enfH8bjH9+pZgEUe5p68P3nmvH2ujK8Z0sFhBDYsSQfzzd242sax1MkfzCEbz59CitLsvD5iAtZIQR2rSjEHw9cgMcfjPqeZ3tHcNP3X0YgIoDT6wQ+d30tPnbV0gl/rvmgocOJ/hEfrlleiO1L8vDmH+zFvb86iPuuWor2wVF89S0rJ/waV9QU4Devt8LtC6IqP0MdTlSYZUF+pjklQ1DaBz0w6XVYXZaV1I3b831u3PHAa+ge9uLaFYWomcJEb6fHj4/++iDWlzvw+RuWR71PGSiiVbq4bpEDNpMee1t61RvwDe3R14e1RfIJvalrWD3Pj98zp9cJWIw6zczcofODGPYEcNUyOTPqsJriZub6RuSbfrkZ2qXd6ysceHDvWRw4148/H2rHR65conmdFak4y4JjGIpdGh6ebOyZ5HqClnDJ6XTmOzCYo5Ro6BiC3WxQl7hGWpxrw56mHnQMehAISVicG52ZA4DeYa9mMHe2dwQnO51JnVQXmgK7GWaDDuf73OoAD60LuTWLsvH7A21Y93+fASDfHco0G5BhNiDDZAi/rkeG2aC+/Y1z/cjNMOHNay+dbGd+phndTg9G/UGWWc6xDYtz8Ls32gAkCOZM8e/MAmP7KVMZzBn0Otwfnl453rYl+XjiWAeCIQlnelyoyo+dormsaGyi5ZaqXDg9fvzrEycAICoDrgwwiVdCdPsDr6HUYcU337lWnaJ36PwA6iZYFj6fbQ4HTEUp/tvLNBtQaDdrTrSs7xhCgd2Muooc/Okj2/F//lKPHzzfAiGAD+2swj9ev0xzmIzFqIfFqE+qR7umyI7HPr4TX/zTMfzn30/hyPlBfOu2tbBbtC8eu50ePFV/ERKAi0OehDcpPf4gfvLSGWxfkofrVhbBqBe466E38M2nGvHPN8vPiaGQhD8ebMO//e9J1Bba8e+3rlGPkauXF+K5U9043eNKeKH4hwNtaOsfxYN3rYqZarxrRRF++WorXjndi2uWj128f3d3E4x6He5/z3q4fQEMuP14ZP95PHm8c14Fc/taenGkbRD3XR37mF5qlkssdyzNR4HdjB+9dwNuf+A1fPHPx1DmsGoGK+NdXpuPB/edxcHWAXxqV03U3+eq0uSCr4l0DI6ixGHBmrJs/OngBYRCUtwJvt1OD97389fhCme8Dp0fmHQwN+oL4sO/OID95/rxyuk+vGl1CVZHTNDdfaILK0uysCgn9nrEZNDhsqpcvNIylplr6HDCYtSpZemV+Rkw6gWaulxqNkrr+s9mMmie/19o7IZBJ7Aj3FeWbTOiY0h7cqsySV0rMwcAdeUO+IMSPvHIYeRlmPCxq5doflwkJTPnGDfNUsnMTXY9QUuPEsxNPTPHMss0MBsN1tPV0OHEinBz63gVuTZ0D3txMjzRTZlwCciZOSD+4vAnwyWWb1p9afXLAfJd0UU5VrQNuNWeQ607Ru/YsAjfeudafOXNK/Dpa2vw/q2Lcd3KYqxb5EBJttys2+PyoqHDiedOdeOR/edx7MIQPrSzKubu7kJWYDdjxCePQWYwN7c2LR7rR4g/ACVxmWWfywe9TsT0LcyUbUvyMOwJoKFjCGd7R1Cdr71Q1m42oOniMLqcHrzrx69i/9l+fOdd66LKy20TBHP17UP4y+F2fOw3h+ANBDHk9uN0z8iEy8Lns/XlDlQXZKiL41Mp3kTLEx1ONRtgMerx7dvW4ofv3YA/f3Q7vvqWldOeCqrINBtw/3vq8JU3r8Duk1146/371DU84z1+tAMhCZAk4IljHQm/7u/2n0fPsBef3CUPY7hqWSE+sG0xHtx3Fnube3Giw4l3/vgVfPHR41hRnIWf3bkp6me6Ojyc4vkEpZYefxA/eK4FdRUOXL2sMOb9W6tzkWHSR5Vanuhw4vGjHbh7ZyVuXF2Mt29YhA/trMKNq4pxsnM44aLn2SRJEv7f307gW083xuyJBOR+uZUlWerzwabKXHzlzSsgScBd2yuTqnLZWpUHU3jwhVJiqVhZmoWWbte0F0l3DI6iNNuKVaXZGPEFNftDAXloy/t/vh99Li9+8+HL4LAZcag1do9aIr5ACB/9zUG80dqPr9+6Brk2E772eINartnr8uJA64BmiaVix9J8nOkdUXv9GjqGsLw4S/19GvU6LCnIRFPXsBp0agdzes1S9Bcbe7BxcQ6ywjdMHFYjhuIMQBkY8UGI2JJIxfpwD3fnkAefua5W/ZqJKLvm4mbmJllm2dLtgt1sQNE09m+mTWbuiWMd+NrjJ/Cnj2xTyyom641z/Rhy+3FtEndb5ov+ER8u/8bz+O6718/a5LHJCgRDONnpxHsvW6z5fmXJ9d5muacg8v+fkpmLt57g7/WdqKtwoNQRm/G7FJTn2tDWP4rcDDPyMkyaJzyrSY/bEqxo0JLozt5ClR/RW8WeublVlZ+BHJsRA25/3D1zSjDnihfMjXiRm2GateN4a7WcWdrT2IO2gVHcrDEpVwiB2mI79p3uxfOnujHo9uGhD27G5TUFUR+nlFm6NZ70QyEJI74gaosysftEF+791UF1NHo69sspLEY9nv/cVTPytavzM/BsRD8cIF9QNXe7sGvFWIAihMBNa2amGkEIgQ9fXo21ixz40MNv4OtPnsRP3h9bRvnXI+1YGy7z/uuRdnz48mrNr+cNBPHjPWewpTJX7eUEgC+/aQX2tfTio785CLcviGyrEd++bR3esaEs5mZqmcOK5cV2PH+qG/dcoZ1x+O3r53HR6cF33rVO82as2aDHFbUFeO5kF6S3rYYQAt/Z3YgsiwH3XB79NTcsdiAYknDswlDUY45nb3Mv7n+hGQ/fvUW9EE6lN84NqNNlH9l/PqrCx+UN4ND5AXxoZ/Tv/87tldhUmRvT/x+P1aTHFTUFGHDLOy8jrSrNQiAkobnLFZXZmqyOwVFsW5KPVWXyY6rvcKJ63PdyeQO486E3cLZvBL+4azPqKnJQV+7AofOxO9/iCYYkfPYPR/BiYw/+4+1rcMeWCuh1wBcfPY7HjnTgbXVl8nEgAdevjH89umOpnDHb19KLd25chBOdzpjJ4rVFdhxsHVBv7owvswTCwdy4aZYXhzw42enEF28cK/1Ueua0ytP7RnzIsZniBubF2RaUZluQYTbg9iRXNN28rhRefwil47LqamZukgmY5i4XlhRmTqvqIi0yc38/3olP/e4Iel3eSR2Y43376UZ84dFjadWke7x9CCO+IPZp7O2YL870jsDjD2F1mfbJT7kj/VJzDyxGHQojsiKRy5zHO9/nRn27EzddYoNPIlXk2tA24EZbv3vCOu7JuNQCOSA6gGNmbm4JIcaexOP2zCUeEtIz7ENenJHhM6HQbsHSwkz84WAbgiEpasdcpGXFdpzpGYE3EMLv790WE8gBY2WWHo2fTSkrum1jOf7j7Wuwp6kHn/n9EeiE3I9CsSrzM9Dr8sHpGbs739Q1jGBIUifozZYtVbm4a3slnjnRhdM90dmg5q5h1LfLu7FuWV+G+nanZsYIAP544AIuOj1qVk5hNenx/dvrYNAJvHtzOZ7/3JV458ZFcS8Er15eiDfODWiOyHf7Avjhiy3YVp2H7Uvjj0LftaIIXU4v6tudONg6gGdPduPeK5fEZCbqlEm1SV6n/fFgG14704+DrVO/rkvkV6+1wm4xYNfyQjx66EJUxuT1M33wByVcMW4EvBACq8uyJ/Uc+YM76vDLu7fEvF0JCKfTNxcIhnDR6UGZw4KawvD0ao1VRN9+uhH17UO4/4469f9lXUUOmrtdSa9H+NrjDXjiWCe+/Kbl6g2k2zaWY+2ibHz9yZNweQN4pqELZQ6rujNTy7IiO/IyTHjldB8uDIxi2BOI6TmtLcpE++AoupxyP6BWMG8zGWJueO1pkjPEVy8fO6/m2EwIhiTNG38Dbl/c1RKKn925GQ/etTmmxDieMocVn7q2JuZvbsqZuR7XtEosgTQI5p5puIhPPHIY6xZlw6ATcU98yTjdM4L+EZ9a7pcOGsOP9fgU9oiliiRJ+IdfHsDfjmqXhCg7zlbHedJcHA5CWvvcqMzLiPoDyLWZ5DJAjczck/XhEstLbCVBpPIcm1rala6DD+aLqGCOmbk5p0w3jN8zN3FmbraD8u1L8tDWL5cOjb8zrrh5bSmuqC3AXz62Pe7deDUzp1GOppYdWQy4Y0sFvv3OdfD4g1hWnBX3d3WpUyZaRvbNKcNPtIZyzbQ7t1fCpNfhZy9HrxL4y+F26HUCN68rxc1rS6ATwOMay8d9gRB+9OJpbKhwYMfS2AzX6rJsHPrqdfj6rWvU3YXxXLO8EMGQhJfD/WGRHn6lFb0uH/7xhtqEX+PqZQXQCeDZk1341tOnkJ9pwgc1RvbnZJhQlZ+RVGmfJEnYF+6rmszN6mSnM3YPe/BUfSdu21iOu3dWYdDtx1P1F9X3v9zcC6tRj42V0y/7tZr0mn+blXkZsJn00+qb6xr2IiTJkzNNBh2WFdtjJlr6AiH89Ug7blpTElXBpYzXP9o28f+PLqcHv3qtFe/fuhj3XjmWcdXpBP7l5lXoHvbiW0+dwsstvbh+VVHCLJJOJ7B9aT72tvSqP/v4myq14T6+w+cH41ZnyJm56HPkG+cGkJ9pUvuTgbESSq1dc30uH3In+BtZWZqVkpvl6jTLSWTmhtx+9Ax7pzXJEpjnwdxzJ7tw328PYXVZNh6+ewsq8zOmHMw5PX41+xPZmDnfKXuNTnTK43K1fPCh/WqT/Uxo6JBHNv/pB+VPtwAAIABJREFU4AXN99e3y82t8S5uHDYj7OETXWS/HCAPHMjLMKHHFTuJ6O/HO7FuUbZmk+2lQllPMOD2q6/T1OTbmZmbT26tK8P7tlbElCYp1DJLT/wBKLOZmQOAbRGlY/Eyc9uW5OGXd29JeHFgU5eGx57Tx4/qfsfGRfjdPdvwrXeunfLjXuiUFRFno4I5eShX+Rw8f+RnmnHbpkV49GC7unA7FJLw2JEOXF4jD9sozLJg+5J8PHa0IyZA+c3rrWgfHMUndsXe/VckW5JVV+5AttUY0zfn9Pjx4z2ncfWyAmxcnHi3VV6mGRsqcvCLV87htTP9+PjVS9VjOOb7VThwpG1gwqCrqcuFXpcXOgE1qJuIPxjCzm+8gAf3np3wY//wRhv8QQnv3VqBbdV5qMyz4bf7z6vvf6mpB1urc2ekvFOh0wmsKMlSp31PhdJ3pvRprS7NRn3HUNTv96WmHgy6/XjbuJ69deXZECK5TOmrp+X/B+/SaNnYuDgHb99QhodfbYUvEEpYYqnYsSQPPcNe/PVwB/Q6ETXNF5ArGAD52larxBIIZ+bGVS+cuujEipLo+QzKDQ2tDGT/yMSZuVRRjiXvJDJzLT3yNf6CzcwduzCIj/76EFaUZOHhu7fAbjFiSUFGTNlCpJeaetDt1F7eqewFAoC987hkcbzGi8PQ60R4V05s0+uIN4CXmnvVfrSZ8EyDfDfrUOsAgholqvUdQ1hZkhW3JlkIofbNVebFXgDlZ5pjMnNDo34cvTCEXSvSp79xJkQGsszMTY9y4W81at9FpdlVkm3Fv71tjTo8YDyH1QijXqA7Tj9tn8uHvFnOsF4WDubyM01xG+qTYQ7/zKMamblhjSW6W6pyp9Vzs9BV5NogBHCu162+TRnKNVcl5R/eWY1AKISH9p0DIPfstw+O4ta6sf2Fb11fitY+t7rnCpBL8v7j76dwZW0BrqqNLdGdLINehytrC7CnsUdtMRly+/HlR49jaNSPz16XeKeWYteKIgyN+lHmsOKOyyriftyGihz0unxqBjse5Trs1rpFOHZhMKlSwKauYbQPjuK7zzbFHXgByL1fv339PHYuzceSAnnq7O1bKrD/bD9auofR1u/Gmd4RzRLoVFtVmoUTHc4pt/cowVyZw6J+vUG3Hx0Ri+r/eqQdOTajusRcYbcYsazIjkPnJ87M7WvpRbbVqLmCAwC+dONyZJj0cNiM2JxENlPpm3v6xEUsKciIGbZWnmODxahDMCQh06x9LpUHoIydIwPBEJq6XDG7/xy2+Jm5AbcvZrH5TDFPITOnJKhqprGWAJjHwdyzJ7oQlCQ8/MEt6pPm0sJMtPa54dfIULl9AXzwF2/ghy+e1vx6Z8JB4M6l+dh/tn/GJkR+4U9H8T8vtKTkawWCITR3u9QTulap5ZG2QXlUdq9L8/eSCs+c6IJRLzDsDeDUuBLVUEjCiQ7nhBcaSkZusUYwpyxzjqSUBczE5LN0Enl3P5U9c5cii1EPu8XArFya0OkESrKtaB+MvSh0+wJw+4KzPsgmN8OE1WVZ03/i1QlYjdqT2iLLLCk5FqMepdlWnO2Vn+eDIQmnOofnpMRSUZmfgTetLsGvX2vFsMePvx5ph82kjxp3f+PqYpgMOjx2RG5hGPEG8PHfHkKOzRh3IMlUXLO8EH0jPhy9MIgnjnVg13f24KmGi/jcdbVYk+TOzRtXF8OoF/jHG2oTZrOU0r6JskH7WnpRnZ+B2zYtQkiSe9gmorR0DHsC+MlL2td6gDy9s2PIg/dtHRvK9s6Ni2DUCzyyv00NJK+ojd8nmCqrSrMw4guqE6kB4MXGblzzXy+qgVoiyvmvJFvOzK0KX2uN/S782H2iC29ZW6quMYlUV5GDw+cHEgaTkiThldN92FqdG/emfGGWBf99Rx3+/W1rkuotK8+1oSLXBkmC5jAZnU6o51F7gr7pyHPkuT43fIFQVIklAHWi8eBodIVXKCRhwO2ftQoOJWCdTGauucsFs0GHMo21XpMxb4O5pi4XFufZoiLqpYWZCIQkzZG/JzudCIakuOns0z0uGHQCt28px6g/iMPTGKQSj9Pjx6OH2vHEsc64HzPk9ifdDHuubwS+QAg3ri6GzaRX/3gjvXGuHwDgD2r/XqartW8Epy4O4/1bK+Xvd7Y/+v39bri8gbj9coqK8G65yrzYgCQ/MzaYO3x+EEJAnfp1qcq2GpEVvqhjZm76CjLNDObSSJnDivYBd8zblYXheZmzW2YJAD9670Z8+13rpv11bCa95tJwl0ZmjiZWlZ+Bs33ysXK214VRf3DWh5+Md88V1Rj2BPDwK+fwxLFO3LiqOKo8MctixLUrCvHEsQ4EgiF89bF6nOsbwffeXZfSrPOVtXLP272/OoiP//YwSrIteOy+HfjEuOEqiVTlZ+DgV6/DrXWLEn7csmI7bCZ9wmDOHwzhtTN92L40D3UVDliMuqT65o63y6WzN68rxUP7zqklrOP96rVWFGdZcG3EJNP8TDOuX1WMRw9dwO4TXSjJtsQt8U6llSXyMaj0uR1tk6vOzvSMJFUl1jnogcNmVKtJVhRnQSfGvt7TDV3wBkJ4W0TGN1JdhQPDngDO9Mavajvf70b74KiaTYtn14qiSe2mVfo94/0dKn1z8W5cWY3RZZZKMmH5uOEryiCegXGZOafHj2BImr3MnGHy0yxbelyoLshMag1GIvM3mOsejmkIVP7wWrpjgxblwD7V6dSs1T7TM4KKXBsuXyqf1PadTn3f3Kun+xAMSWjpHo6b+fvec0249Yf74jb1R1L65VaWZmFVaRaOXYhNlR84N6DuLGrqmlo/Yc+wF5/7w1HNE+PuE/K45w/uqERptgVvnIs+QSsB5qo4kywVK0rsMOiEZl1wfqYJPcPeqP9vR9oGUFOYGXfx6qWkPNcGQzhLQdPzzk2LYvoKaP4qy7GiYzD2vKTc/JmLQTbluTaUpWBVisWo19wzl2jvEsVXmW/D2R4XJEma0+EnkdaVO7CtOg/ffbYZw56A5gX3W9eVodflwxcePYY/H2rHJ3fVYNuSicf6T0ZOhgmXVcl7Er/y5hUJh/MkkswOLr1OYN2ixCPxj7QNwu0LYufSfJgNemypykvqmux4uxOryrLwuetq4QuG8MMXYrNz53pH8FJTD+7YUhGTQXrvlgoMuv14/lQ3rqgpSFnmM5GaokwYdAInOodwrncEd//iDeRlmmA3G6LKa+NRdswprCY9lhRkqhMtHzvSjopcGzbEWVuiZkoTDKV5Jfy7374ktZnKK8JlrPFuytcWydeDcYdgmeUyS+Xa8FSn3HY0/jpSqd4bckdn5vrCC8NnKzNn0usgxOSmWbZ0u6Y9/ASYp8GcNxBEa59bjdoVSjCn1TfX0C6fvJ2eQFQtseJMzwiqCzKQbTNiTVn2jIz6V76mPyjF7e070jYIbyCEvRqTpcaLPHBXl2XHDEEJBEM4dH4Ab1lbAiHkevKp+P5zTXj00AX893PNMe97pqELy4vtKM+1YXNVLvaf648Kuuo7hmDS6yYsO7p5bSle+MerUJhliXlfgd0MbyCkXsRIkoQjbYNYX84x3ABQU5iJJSm4c0PAx65aivdvq5zrh0FJKnVY0TXsibk51juHmblUiZuZYzA3JVX5mXB6Ahhw+9HQ4YTJoJv2UIFU+MhVSxAMSSiwm7FdI0i7enkB7BYD/nyoHVurc/GJa5LPlk3Gj9+3EXu/eDU+fHl10iPYp2rDYkfC5eF7m3shBLCtWg4edizJQ0u3C11xZh4AcjbvZKcTa8qyUZmfgXdtKsdvXm/FhYjMvS8Qwvefa1arsMbbGh6EAgCXz0KJJSDftFlamIlXTvfhzof2IyRJ+OXdW7C+woHDSfSytQ+OxuzZXVWahYYOJ7qdHuxr6cUt60vjBqbV+RnIthoTBtf7WnpRaDdjScHUdjjHc+PqYvzunq3YUqU9YKc2PAQl3rnOatIjJI1luk5ddKI6PyOmzNds0MNm0sf0zPWHg7nZGoAihIDFoE86M+f2BXBhYDQl56l5Gcyd6RlBMCShZlwwl2E2oDTbojnRsr5jSI3OT40rtQyGJJztG1GnLe5Ymo8jbYMY9iS3eyNZe5t71d6wkxrlnspybQB47mR3zPvHO3VxGFXhA3dNWXbMEBT5ZBnE5TUFqMi1oXkKmbm2fjd+/0YbMs0G/P6NtqgTY5/LiwOt/eqo282VuegZ9qK1L6LJvN2JZcX2uEMMFDqdiNvzlT9ucXhrnxsDbj/qKi7tfjnFP9+8Cg9+cPNcPwyiWbfIYYUkyYtiI/WFM3OzPQAllaymOD1z4TJLDumZnKp8+fnlbO8IGjqGsKzIrtlDNNuuqMnHruWF+PDOKs0gymzQ412bylFoN+P7t9fN2E27bJtx1v5eNlTkIBiScPyC9kj+V073Ym1ZtloeF7lkOp7mLhd8gZCaUfzkrqUQQuD7zzarX/NN338Jfzncjru2V6JI48axTidw1/ZKZJj02JHiLFQiK0uzcPj8ILqcHvz8rs2oLshEXbkDjRedcQNeRcfgqDr8RLG6LBsXnR489Mo5hCTglvXaJZaA/DPXVcTPlEqShFdP92HH0vyUZyqFENhanRf36yq9b/Y4ZZbKehrlPHnq4jCWx1nm7rDKi8MjzXYwB8hDUJLNzCmDGRdsMKdkmJQUbKQlhZkxWS9fIISmrmHcvE6u5VXKExXtA6PwBULqXYcdS/MRDEnYP67/azraB0dxpncE77tsMUwGnWZfXEuPCx5/CJlmA15o7J5wulFjl1Md57omfAKLHIJyoFV+/Jsqc1BTaJ9SZu77zzVDJwR+9aEtEBBRw1ueO9WNkARcH27YVu6u7A/36UmShPqOobjLwpOl9DApd9sPt8knnbo4ZQOXmtwMU0rKuojSjdIUfmEwum9OKbOc7dUEqWQ1xsnM+QIwG3QT3iCjaMqkZDmYc855iaVCCIGf37U5anfXeP/nphV46QtXawYg6ahOHYISm3lyeQM4fH4wqj9rZUkWcmzGhCsKlJYO5VqoJNuKD2xdjEcPXcC9vzqA9/z0dfiDEh66azO+8paVcb/Ondsr8cqXd81aHxUgB7c6Adx/xwa17LGuIgchCTgWJ+AF5OEmTk8gJjOnTJx8cO9ZrCnLnjAY2BBeHu7USGA0dg2jb8SX8tLeZJRkW3D9yiJsrdbO3FnDLUQj3gCcHj8uDIzGrDhQZNtMc56ZAyBn5jRWzmgZm2S5QIO55i4X9DqhucdnSUEmTne7okr9mruH4Q9KuKwqD+W51pghKKfDjZ9KZm7j4hyYDbqkd5skQymbvHJZAZYV2TUXk9eHS0Hv2l6JXpcPxxIsAnd5A2jrHztwqwsyY4agHDg3gDKHFSXZVtQWZeJs78ikpnSe7nHhz4cu4P1bF6OuIgd3bCnHHw9cwPlw5u2Zhi6UOazqk+LSgkw4bEZ1CEr74CgG3f5pN5mPz8wdOT+IDJN+2hPjiCi9KTcxxvfN9bp8sJsNMeOu04nVFKdnzhNgieUUlOfaoNcJvHK6N/y8ND+CuWTodCKtj+XxcjNMqMyzaWaD9p/tQyAkYWdEMKfTCWxfko99Lb1x99Mdbx9CptkQtd7oo1ctgdWoxwunevDJa5bimc9cgauXF2p+vkIIMa21IlNxx5YKvPrlXbg2YpLpunAbSaK+uc5wRUJJTJmlfM2VaPBJpLoKByRJe3m4snd5ouEnM0EIgQc+sAnXLNdeQZWh7uMMoimcpIkXzOXYjBgaN81yzjJzgeQyc83dciuV1pT3yZqXwVxT1zAq82ya42+XFmZixBdUD3JgrF9uVWkWlhdnxZRZKqnM6nBwaDHqsakyJ6V9cy83yzXHNYWZWFFix8nO4ZiTUn37EGwmPT64oxI6ATx/sivu12sMH7jLiuUnJL1OYGVJlpqZkyQJb5zrV/d91BbZEQhJUUtTJ/K9Z5thMerxkavkO4Yfu3op9DqB/36+GW5fAC839+C6lUVqilynE9i0OBcHWgfCP4/8e57u/iMlmFPuth9uG8TaRQ72iBFd4oqz5UxF+0D0CO9elzdqCXw6StQzx7UEk2fU61CeY8XuBvl5deUcT7K81G0Ij8Qffx20t7kPZoMOG8atHdq+NA8XnR6ciXMNc7x9CKvG7Q3MyzTjzx/bgWc/eyU+e/2yeRsQ63UiJuuqBLyJJqu3j9sxp8i2GlGea4VOQK1IS2R9uUNeHq4xBOWV072ozEvNUKdUs0Vk5pSKu7hlljajZmbOZtLP6nFhMeiTLrNs6Zan9qeiCmNeBnPN3a6Y4ScKrSEoDR1DyDDpUZmXgRUlWTjbOxL1yzzd40K21RgVne9Ymo/GruG4o20B+QD63rNNE/bWhULyjo6d4ZrjlSVZ6B/xxSy7Pd4uL9fOyzRj4+IcPJugb65R4y7E6rJsnOiQVzC09Y+ie9iLTZVyelr5fSVbanmy04m/He3A3Tuq1GCqKMuC921djD8fuoCHX2mFNxBSSywVW6pycLZ3BN3DHjR0DEGvE3HvlCQrN8MEnZAv0Dz+IE50OLGeJZZElzyLUY8Cuxnt48os+1y+tC6xBOSx28zMpVZlfgaGvQEIIU9QprlTt1h7efi+ll5srsyNucDemaBvLhAx/GS8ZcV2VGisPEoH68vlISjxspHKHrrxZZYA8NZ1peFey4lLc+0WI2oL7TGZ0kAwhNfP9GPbLPYPToYSzI36gjh10Qm7RZ6boSXbatLsmZvNrBwAWIy6pAegpGqSJTAPgzmPP4jWvpGY4ScKpTY4cghKQ4cTK8N3bFYU2xGSEDUM5EyPC0sKMqKaMJXm11cTjMN99mQXvvdsM/7rmaaEj/lEpxP9Iz7srJG/prKdPrLcMzhuufY1y4twotOJziHtpZGNF53INBuwKGKR4JqybIz6gzjd41L3y20OB3PVBRnQCaA5yWDuO7ubYLcY8A+XV0e9/SNXLoHZoMe3nj6FbKsRm8dNIVK+34FzA6hvH0JNYea073rodQK5GWb0DHvR0DGEQEhCHSdZEhHCu+bGLdftG/HO+sLwVLOadJrDD1xeBnNTpbRmVOdnRO1zo9mnjMpXeuABoHvYg8auYc2Svorwyg+tYK652wVvIJT0gvN0sb7cge5hb1SlWaSOwVHodUIzYPv8Dcvxn+9Ym/T32rDYEbM8/Hj7EIa9AXUf3Hyj/A2P+II41TmM5cX2uMNUHDYjhtz+qMB4LoI5c5KZOV8ghHN97pRN3J13wdzpHhdCkvbwE0DeSZZtNarBXCi8KFypIVZSsJHTJOW1BNFfb3VZNrIshoSllkfb5JLGX73WqmbKtCiLH5U7S9qPQV5iqtxZUpZZPn9KOzt38uIwaosyow5c5UR2/MIQDrT2I8tiUKN6i1HOTCaza+61M33YfaIL91xerU6TUhTYzfjA9sUIScA1ywtjpoGtLsuGxajD/rP9qO9wpmwpa4FdXhyujOplZo6IAO1dc70uX1qvJQDkC5W4ZZYM5qZECebmelk4yZMKbSY9nqq/iD1NPdjT1INfv9oKAFH9cgohBHYuzVf39UZS2kum29Ix3yiDYuL1zXUMelCcZUlJy0ldeQ6cngB+/NJpdbaCsl9uW/U8DebMcqLA7Qug8eIwlhfH74N1WI3wBUNR59Q5CeaSzMy19slT+xdsMKdk1OKVWQoh711TyizP9o3A7Quq030W59pgNerVASTDHj+6h72oHrc/Q68T2LYkD/ta+uKmuI9eGMSyIjsyzQb83781xP24vc29WFZkV3eoZVuNKHNYoyZaKicjJSBbWpiJ8lwrntcotZQkCY0Xh9V+OcWSgkxYjXocbx/CgXMD2Lg4J6p+vKYoE03diTNzB1sH8OGHD2Bxng0f3Fml+TH3XrEE68oduH1z7J4Wo16HuvIcPN1wET3D3mlPslQoi8MPtw2izGFNqnSAiBY+JTOn3FEOBEMYcPvSei0BIN+A8/hDMVON2TM3dcpwjHQafrJQGfQ6bKrMxd/rL+LOB/fjzgf347+fb0F+plm9XhtvR00+nJ5AzE32+vDwk6oUDIqYT1aUZMFk0MXtm5PXEqSml+2mtSW4orYA33yqEdd/dw+eqr+IfS29WF5sn7fnUqXMsrnLhWFvAMsTlE47womJgYi+uf4RH3Jtc5GZmziYa1YnWaamHHzePWM0dQ3DoBNRE4vGW1KQgedPydMjGzrGhp8A8pCOZcV2nOqUg5qx4Sex0e+26jw83dCFCwOjMTvQ/MEQ6tuH8L6ti7E4z4Z/fqwBT9VfxJvWRDebevxB7D/Xj/ddtjjq7StKsqIyc8fbh2Ax6tSePyEEdi0vwiP7z2PUF1RHsAJAl9OLoVF/TM2/XiewsjQLe1t60dLtipliVFtkx7Mnu+ENBDWHxxxs7cedD76B/EwTHrlna9y7v7kZJjx23w7N9wHA5qpcvHpGvqOTqjtlBZlmnOkZQa/Lx6wcEanKHFb4AiH0jnhRaLeg3+2DJAEFaZ+Zk8/RnkAwqiRwhJm5KVtf4cCOpXm4bqX2dDyaXT+4oy5mL3CZwxo303T9yiKUZFvwnd1NuLxmbO/Z8fYhtZVmITEZdFhdmhU/Mzc0io0p2rebaTbgl3dvwYuN3fj6kyfxkV8fBAB8KM5N/flAOS8qvX6J5jNkW+Xng0G3Tw2A56xnLokyS+XvYnyiaarmXWauqcuFqvyMhNNdlhZmotflxZDbj4aOIZj0uqjodkWJvBpAkiScCa8lWFoY+wu7LJxa1to313hxGN5ACOvKHXjPlgosL7bj3/73ZEzD+hvn+uELhHB5TXTZwMoSe9QgloZ2J1aWZEWdxK5ZXghvIIRXz0TfhToVziou08hOrinLVg8CpX9NUVNkRzAkqQFspAPn+vGBn+9Hgd2M392zDSXZU7/bsyX8feUm89TcAS2wm9E5NIr2wVH2yxGRSnliViZa9oX3Uc7Xu8nJshqVEqLo55RhDkCZsiyLEb/58NaYtgqaG9lWIzYuzol6KY4zwAKQs9Wf2lWDI22D6oC4RMNPFoL15Tk4dmEI/mB0NicYknBxyKM5/GQ6rlpWiCc/eTm+fusarCjJwi3rS1P69VNJueGlBLvxKvaAsczcUDgzN+oLYtQfRO4s3/SzGPVJlVm2dLuwKMeast7eeRfMNXcPJ/wfBkQMQelxoaHdidrizKjgb0VJFgbdfnQ5vTjTMwK9TqAiNzaYW1ZkR7bViNfPxg5BURY5rluUDYNeh3+5eRXaB0fxk5dOR33c3uZeGPUCl41beriyNAshSQ4KQyEJDR1DMSejy6pzYTPp8dy4Ukt1BKtGfbCSCTPqBdaOawZW+gzHT7Q80jaIOx/cj6IsCx75h60JT6bJqKtwqHsAU3XRkZ9phlJtxGXhRKRQFocrfXNKMJf+A1DGJrUpfIEQvIEQgzm6ZL1j4yJU5Wfg2083IhSS0NLjgscfWrDBXF2FA95AKGYuQ6/LC39QSnkwB8glsO+5rAJ//9TlWLto/l5vGfU6mPQ6uH1BlOdaYbfE3w+oBHPKRMu+EXma/OyXWeqSGoByuselVuqlwrx6xghJwPl+N962PvESRHU9QbcLDR1DuH5lcdT7lSDo5EUnTve4UJ5j1cz06XQCmytzNTNzR9sG4bAZUREuv9y2JA9vXluCH714Gh5/CMpckieOdWJDRU5MdL0iYghKhtmAEV8wpiTRbNDj8pp8PH+qG4FgCIbwsJHGi8MozrLEDCcBoJ7Q1pRlx0yRrMrPgF4nooI5SZLw1b/WI9tqxCP3bI3ZdTIVGWYDrl1RmJJFh4qC8M4oo16weZ2IVMrFjLKeQNlHme4DUJTMXGTD/ohXnm7Jnjm6VBn1Onzmulp88pHD+NuxDnVYx0IbfqJYH65EOnx+IOpnbFfXElza8wOsJj18oyEsK0pcBZZjU8os5WBuYET+7+yXWSaXmRsa9WtW303VvMrMeQNBSFLiVCoALMqRl+y91NyDAbcfq8YN4VgWrqs92enUnGQZ6bKqXJzrc6PLGT0t7eiFQaxb5IiaJvlPN61AcbYFD+49i5+/LL/0urwxvWsAUJ5jQ4ZJjxOdTtSPG34S6eZ1pegc8uDWH76illeeujgct9FzSUEG8jNNuKK2IOZ9ZoMeVfnREy13n+jC8fYhfPq62pQEcoqfvH8T/ummFSn7espd9hUlWfN28ScRzb5sqxF2s0Ets1SCuXTPzNk0MnMuJZhjZo4uYW9ZU4LlxXZ8Z3cTDrcNIsOkR3X+whp+oliUY0V+phmHx/XNJdoxdylRzpMT7Y3MtiqZOblyQ83MzfpqguQycx5/9KyM6ZpXzxhK02C8tQQKvU6gOj8Dz57sAhA7uSpymuTZ3hHNMbgKpTzy9bP9eOs6uXbY7QugqWsY16+KzviVOazY8/mrk/pZdDqB5eEhKCa9DmaDDks1gsq3rC2FTgh89a/1uPkHe3Hf1UtxutuFK2q1H7NBr8Nzn71KHdk6Xm1RpjpFMxSS8J3dTajKz8DbNQLO+UTJzLFfjojGK8sZ2zXX6/LBqBfISvPslfJE7mYwRxRFpxP4/A3L8KGHD6BjsA115TkLbviJQgiB9eUOHDnPYE6LEswlWksAyBkxi1Gn9swNuOWgbi4yc4GQFFVtp2XUF1SrM1JhXmXmPIEQjHqByiTuwCwpzFTLHbWGcKwoseOlph54AyEsSbDHYWVJFjLNBuyP6Jurb3ciJAHry6eX1l9RIk/VPNY+hBUlWXH/x960pgS7P3slblxdgu892wxfMJR4ao/NGLP/TVFTaEdrvxsefxBPNVzEqYvD+NSumoQH1XxQnmvF0sJM3LC6eOIPJqJLiryeQOmZ8yIvwxx3eWy6GCuzHFsc7mKZJREAeUDchgoH/EFpwZZYKuoqHDjTO4LBcAACyD3CdrMBWQn6xC4FSgvTsgTXxAqH1aSWWaqDsjJmt4JmQnFbAAAXy0lEQVTDHG7pSlRqKUkS3CnOzM2rK3yvP4iq/Iy4gUokJctVnZ+hOQ1meXEWnJ6A+jHxGPQ6bFycE9U3dzSc7p5uY+jKkmwMewM42DowYfNuboYJP7ijDj9+30ZcvawAl9fEllEmo7bIDkmS93J8d3cTlhZm4uZ183dakcJmMuDZz16J7UviZ1GJ6NJU6rCifWCsZy7fnt79csDYRcqob+xJn5k5IpkQAp+/YTkAYFNlasbzz1dKRdLPXj6r9gi2D45e8lk5QM7MmQ06VObZJvxYh82olln2j/ig1wnYZ/nGmNImlKjU0hsIQZKwcMssPf4QapJsCFQmWsYblhGZrZtoTPGWqlx86+lGdSfFkQvy4urp9mQoNb7BkJT0JKYbVxfjxmlkp5QS1e/sbkRztwv3v6cu7k4XIqJ0UJZjhdMTwLDHj74R36zfbZ0JY6sJIjJz4RuQs30BQjQfbVuShz2fvwrlORNfyKezTZW5uHZFIe5/oQV/O9aBL924HB2Do5f88BNArsoQAklVl2VbjWMDUNw+5NhMs16eazFOnJlTAr20KrMUQtwohGgUQrQIIb6U6GN9wRBqk9yGPhbMadfRKgNE7BYD8ieYerY13DenZOeOtg2qE4amY1mxXZ16OVtlApX5GTDqBV5o7MHyYjtuWl0y8ScREc1j6q65wVH0uXxpP/wEGLsrG3kHV8nMZTAzRwQAWJyXsWD75RQmgw4/u3MzfvHBzTAbdPjobw6hocPJzByAf791DX525+akPtZhM2JodKzMMjdj9ktUzYaJM3NKn3TaBHNCCD2A/wHwJgArAdwhhFiZ6HMmGn6iWF5sx7/cvBLv2lSu+f7KvAyYDTosKcicsLdiTZkDFqMOr5/tQ5/LiwsDo1g3zX45QC6jqcqTF6DXJPlzTZdRr0NVuKz009fWLviTIBEtfMquufaBUfS4vBPeoEsHmgNQPCyzJLpURS70Xpxnw9bqvLl+SHPOatInfT50WE3q4JMBt2/Wh58AY5k5jz9+Zk5ZR5NOZZZbALRIknQGAIQQvwNwC4AT8T4h2TJLIQQ+uKMq7vv1OoG3rS9LapiKyaDDhgq5b+5YjbxGIFWLFK9cVoALA6NJ9QGmys6lBXBYTbhhVdGsfU8iopmiZOYau4bhC4QWRmZOY8+cmpnT6AMnooVPWej9nssq5vqhpB2HbazMsm/El3CQ4ExRMnPeQPzM3OgMZOZm+hmjDEBbxL8vALgs3gcLIKkmx2R9451rk/7YLVW5+P5zzXi5uRc6gaR73CbyLzevSsnXmYx/vnklJElK+2lvREQAUJBphkmvw7E2+WZbui8MB+QbjmaDLmbPXIZJz4oKIqJJyrYZ4Q2E4PEH1RkYs808R5m5mU4XaT0jSVEfIMQ9QogDQogDdkNwzkbob6nKhSQBfzjQhppCe9r3LDCQI6KFQqcTKHFYcPSCPGk4bwFk5gD5yXx8mSXXEhARTV6OTQ7eel1eDI36kTsHg7KUaZbJZOZsaRTMXQAQ2dS2CEBH5AdIkvSAJEmbJEnatLgod4YfTnwbKnJg1Au4vIGU9MsREVHqlDms6BySd80thJ45ALAZ9dFllr4A++WIiKbAYZUHnrT2uSFJQK5tLgagJJ+Zs6TLABQAbwCoEUJUCSFMAG4H8PgMf88psRj1WBfuk0tVvxwREaVG5GS3hdAzBwAWkz66zNITQOYlviSYiGgqssPB25keFwAgdw6eJyaTmUubaZaSJAUAfBzA0wBOAviDJEkNM/k9p+Oy8IqCVKwlICKi1CmLCObmohdiJthM+pgBKJnm1D3BExFdKhxW+XnhTO8IACDXNgc9c+HMnDeJzJwthYOuZryeQ5KkJwE8OdPfJxXee9liGHS6qIXjREQ095T1BA6bcVanA88kq1EfszQ8P3NhL0gmIpoJDjUzFw7m5mQ1QXjP3ELKzKWbUocVn7muFnpOEiMimlcWhTNzC6XEEgCsJgNGI+7gypk5llkSEU2WEsydDWfm5mLqsVpmmUzPnCl1IRiDOSIimveUnrm8BVJiCYQHoERm5rwB2DnNkoho0qxGPUx6HS4MuAGMBXezaWwASuLMnF4nYEphhQmDOSIimvdKHBYACy0zN9YzJ0mSvGeOPXNERJMmhEC2zYiQBNjNBnWB92wy6nXQ60TiMkt/EFajPqUrxBjMERHRvGc26LF2UTZWlS2cnmZrxDRLbyCEYEhimSUR0RQp6wly5rCCw2zQJSyzdPuCKV1LAMzCABQiIqJUePzjO+f6IaSUPABFDuaGPXK5JZeGExFNjbI4fC4nHluM+oSZOY8/mNKF4QAzc0RERHNCWU2glFgC4GoCIqIpUnbNzWVv9USZuVFfMKWTLAEGc0RERHPCYtRDkuQSS5eSmWOZJRHRlMyHMks5M5egzNIfhIWZOSIiovSnlNqM+oIRmTmWWRIRTYVjnmTmEk2z9PiCsDEzR0RElP6UUhu3fyyY42oCIqKpcYR75uZ0AIpRD2+CzNyoPwgrM3NERETpzxqVmfMDADKYmSMimpLscJnlnA5AmSAz5/YF2DNHRES0ENhMcuAmB3Pykz/LLImIpmZelFlOkJnz+EMpX03AYI6IiGgOKHdnR/1BdQAKyyyJiKamptAOk0GHJQWZc/YYLAYdvAkyc6MzsJqAzxpERERzQCmzdPsCcHn9MOgEzAbeYyUimoplxXY0/uuNEELM2WOwGPUTl1myZ46IiCj9qZk5n5yZyzAb5vQihIgo3c31OdRs0MUtswyFJJZZEhERLRTqagK/3DPHfjkiovSWKDOnBHmpLrNkMEdERDQHxsos5WmW7JcjIkpviTJzbp/cG81plkRERAuAEsx5wnvmuJaAiCi9KZk5SZJi3jcaztgxmCMiIloA1KXhPpZZEhEtBBajDiEJ8Adjgzml/JIDUIiIiBYAo14Ho16EVxP4kckySyKitGY2yIGaNxDbN+f2MTNHRES0oFiN+vDS8ADszMwREaU1i1EOrTz+2L65UR8zc0RERAuK1aSPWk1ARETpK1FmbpRllkRERAuLzWSAyxfAiI89c0RE6c6cTGaOZZZEREQLg8WoR++wFwC4moCIKM0pmTmtXXOcZklERLTA2Ex69LjkYI6ZOSKi9Kb0zGntmlOCOS4NJyIiWiBsJj16wpk59swREaU3Szjr5tXKzIXLLC0M5oiIiBYGi1GPYU8AALiagIgozZkNCTJz7JkjIiJaWCLLbbiagIgovSmZuXg9c0a9gFGf2vCLwRwREdEcibxDyzJLIqL0pmTmPHGWhltSnJUDGMwRERHNmch9QxyAQkSU3pTzuMsbG8x5/MGUl1gCDOaIiIjmTOQTO1cTEBGltyyrEQAw5PbFvG/UH0z5JEuAwRwREdGciXxiZ5klEVF6sxj1sBr1GBr1x7yPZZZEREQLjPLEbjboUt4UT0REs89hM2LQHRvMefzBqNL6VOEzBxER0RyxmeRsHEssiYgWhmyrEYMamblRH8ssiYiIFhTliZ3DT4iIFgaHzYghjczcKAegEBERLSxKmSUXhhMRLQxyZk5jAAp75oiIiBYWJTOXYWIwR0S0EDisJs2eOU6zJCIiWmCUZnj2zBERLQwOW5yeOZZZEhERLSzKEzt75oiIFoZsmxG+QAgef/TicLcvCAszc0RERAuHkpnjjjkiooXBYTUBQFSpZTAkwRcIwWZM/bmewRwREdEcUadZssySiGhBcNiMABA1BEXJ0llNqQ+9GMwRERHNEeUurZ2ZOSKiBcFhDQdzEZk5ty8czLFnjoiIaOHIshpw1/ZK7FpRNNcPhYiIUiBLI5gby8yl/sYdbwUSERHNESEEvvbWVXP9MIiIKEWUMsuhiDLLUT8zc0RERERERPOawxY7AEUts0ynnjkhxNeEEO1CiCPhl5tm6nsRERERERHNtQyTHgadwFDErrlRtWcu/cosvytJ0rdn+HsQERERERHNOSFEzOLwsZ45llkSERERERHNW9lWI4YWyDTLjwshjgkhHhRC5Gh9gBDiHiHEASHEgZ6enhl+OERERERERDMn22qM2jOnDECxzbfMnBDiWSFEvcbLLQB+BGAJgPUAOgH8l9bXkCTpAUmSNkmStKmgoGA6D4eIiIiIiGhOOWymqAEoSjBnmYHM3LR65iRJujaZjxNC/BTAE9P5XkRERERERPOdw2pE48Vh9d+jvgCANOuZE0KURPzzVgD1M/W9iIiIiIiI5oNsm3HcNMsQgJnpmZvJaZbfFEKsByABOAfg3hn8XkRERERERHPOYTXB5Q3AHwzBqNdh1B+EyaCDXidS/r1mLJiTJOn9M/W1iYiIiIiI5iOHzQgAcI76kZdpxqgvMCNZOYCrCYiIiIiIiFJGCeaUXXOj/uCMTLIEGMwRERERERGlTJY1HMy5lWAuxMwcERERERHRfOcIB3ND4V1zo77AjKwlABjMERERERERpYzDZgIQmZkLzshaAoDBHBERERERUco4xpdZ+tgzR0RERERENO9lqWWWcjDn9gVZZklERERERDTf6XUCWRaDGsx5/EEOQCEiIiIiIkoHDpsJg+7wABSuJiAiIiIiIkoP2VajumeOZZZERERERERpwmEzqgNQPJxmSURERERElB6yrUYMjfrhD4bgD0qwMTNHREREREQ0/8mZOR9G/UEAYGaOiIiIiIgoHTisJgyN+uH2ysEce+aIiIiIiIjSgMNmREgCel1eAOA0SyIiIiIionSgLA7vGBwFAO6ZIyIiIiIiSgeOcDB30ekBAFiYmSMiIiIiIpr/HDYTAKBzSA7mOM2SiIiIiIgoDThscmauUymzZGaOiIiIiIho/lPKLJXMHHvmiIiIiIiI0oAyAKUr3DPHzBwREREREVEasBj1sBr1zMwRERERERGlm2yrEd5ACAAzc0RERERERGlDGYICABYDgzkiIiIiIqK0kB3um7MYddDpxIx8DwZzREREREREKaZk5maqXw5gMEdERERERJRyDqu8ONxmMszY92AwR0RERERElGJKZs5inLmQi8EcERERERFRimUrZZYzNMkSYDBHRERERESUcsoAFJuRZZZERERERERpQ+mZszAzR0RERERElD7GplmyZ46IiIiIiChtqGWWnGZJRERERESUPsamWbLMkoiIiIiIKG04bHLPHJeGExERERERpZEMkx75mWaUOiwz9j1mroCTiIiIiIjoEiWEwHOfvRI288xl5hjMERERERERzQBlcfhMYZklERERERFRGmIwR0RERERElIYYzBEREREREaUhBnNERERERERpiMEcERERERFRGmIwR0RERERElIYYzBEREREREaWhaQVzQojbhBANQoiQEGLTuPd9WQjRIoRoFELcML2HSURERERERJGmuzS8HsDbAfwk8o1CiJUAbgewCkApgGeFELWSJAWn+f2IiIiIiIgI08zMSZJ0UpKkRo133QLgd5IkeSVJOgugBcCW6XwvIiIiIiIiGjPdzFw8ZQBei/j3hfDbYggh7gFwT/iffiHEsRl6TJeabABDc/0gFogKAOfn+kEsEDwuU4fHZerwuEwNHpOpw2MydXhcpg6Py9SZ6LhcnOwXmjCYE0I8C6BY413/R5Kkx+J9msbbJK0PlCTpAQAPhL9XjyRJm7Q+jiZHCPGAJEn3TPyRNBEel6nD4zJ1eFymDo/L1OAxmTo8JlOHx2Xq8LhMnVQelxMGc5IkXTuFr3sBQHnEvxcB6Eji8wan8L1I29/m+gEsIDwuU4fHZerwuEwdHpepwWMydXhMpg6Py9ThcZk6KTsuZ2o1weMAbhdCmIUQVQBqAOxP4vOYuk0RSZL4B5c6PC5ThMdlSvG4TBEelynDYzJFeEymFI/LFOFxmVIpOy6nu5rgViHEBQDbAPyvEOJpAJAkqQHAHwCcAPAUgPuSnGT5wHQeD9EM4XFJ8xGPS5pveEzSfMTjkuajlB2XQpI0W9mIiIiIiIhoHpupMksiIiIiIiKaQQzmiIiIiIiI0tCMB3NCiAeFEN1CiPqIt60TQrwqhDguhPibECIr/PZKIcSoEOJI+OXHEZ/zbiHEMSFEgxDimzP9uGnhmswxGX7f2vD7GsLvt4TfzmOSUmaS58r3RpwnjwghQkKI9eH38biklJnkcWkUQjwcfvtJIcSXIz7nU0KI+vBx+em5+Flo4ZjkcWkSQjwUfvtRIcRVEZ/D8yWlhBCiXAjxQvjc1yCE+FT47blCiN1CiObwf3PCbxdCiP8WQrSEj8ENEV/rG+HzZb0Q4t0TfnNJkmb0BcAVADYAqI942xsArgy/fjeAfw2/Xhn5cREfnwd5sV5B+N8PA9g104+dLwvzZZLHpAHAMQDrwv/OA6DnMcmXVL9M5rgc93lrAJwJv87jki8pfZnk+fI9AH4Xft0G4Fz4eX01gPrw2wwAngVQM9c/G1/S92WSx+V9AB4Kv14I4CDkZAbPl3xJ2QuAEgAbwq/bATQBWAngmwC+FH77lwB8I/z6TQD+Dnk391YAr4ff/mYAu8PnygwABwBkJfreM56ZkyTpJQD94968DMBL4dd3A3jHBF+mGkCTJEk94X8/m8TnEGma5DF5PYBjkiQdDX9unyRPZuUxSSk1jXPlHQAeCb/O45JSapLHpQQgQwhhAGAF4APgBLACwGuSJLklSQoA2APg1pl+7LRwTfK4XAngufDndUPe77UJPF9SCkmS1ClJ0qHw68MATgIoA3AL5BsFCP/3beHXbwHwS0n2GgCHEKIE8vG6R5KkgCRJIwCOArgx0feeq565egBvDb9+G6IXjFcJIQ4LIfYIIS4Pv60FwPJwGaYB8i8i8nOIpiveMVkLQBJCPC2EOCSE+EL47TwmaTYkOlcq3o2xYI7HJc2GeMflnwCMAOiEnPH4tiRJ/eGPv0IIkSeEsEG+I83jklIt3nF5FMAtQgiDkHcfbwy/j+dLmhFCiEoAdQBeB1AkSVInIAd8kLPDgBzotUV82oXw244CeJMQwiaEyAdwNSY4LucqmLsbwH1CiIOQU5G+8Ns7AVRIklQH4LMAfiuEyJIkaQDARwH8HsDLkEs3ArP+qGkhi3dMGgDsBPDe8H9vFULs4jFJsyTecQkAEEJcBsAtSVI9APC4pFkS77jcAiAIoBRAFYDPCSGqJUk6CeAbkLMlT0G+WOFxSakW77h8EPKF8gEA3wPwCoAAz5c0E4QQmQAeBfBpSZKciT5U422SJEnPAHgS8nH6CIBXMcFxaZjiY50WSZJOQS5fgxCiFnJ9KCRJ8gLwhl8/KIQ4DTkzckCSt87/Lfw590B+wiBKiXjHJOQngD2SJPWG3/ck5Dr953hM0kxLcFwqbsdYVk75HB6XNKMSHJfvAfCUJEl+AN1CiH2Qy9nOSJL0c/z/du7fNYogCuD49xnTqI2oAUEsBNvYCoKK4D8gpFOD2IipbEUQAhGrgIVgoaTyB4IBhQRiK1iLiopKComKNqKFjcFnsS8QD+7gzJ3hwvfT3NzMzt4Wj7l9M7MLt6rPFZqxVeqZDveWy8CFleMi4inwrtocL9UzETFMk8jdzszZqv4SEbsz83Nto/xa9Uv8veK2B/gEkJlTwFSd8w4Vr+2sy8pcRIzU5ybgEnCjvu+KiKEq7wP2A4stfbYD54Gb///KtVG1i0lgARit5e7NwBHgVUsfY1J90SEuV+rGgHtt+hiX6osOcfkBOFZvadtK81D/m5Y+e4ETtExCSGvV4d5yS8UjEXGcZlXO/3H1VEQEzYTV68ycXtX0CBiv8jjwcFX96RovDwLfK+Ebiogddc5RYBR43Om3+74yFxF3gaPAzohYAi4D2yJiog6ZBWaqfBiYjIhlmtmRc7XfHuBaRByo8mRmvu33tWtj6iYmM/NbREzTvCUrgfnMnKvjjEn1TJdjJTTj5VJmLracyrhUz3QZl9er/JJmC9FMZj6vtgd1g/ILmKgtbtI/6TIuR4CFiPgNfAROrTqV46V65RBNbL2IiGdVdxG4CtyPiLM0E15j1TZP8/zwe+AncKbqh4EnTW7ID+BkrS63FfUaTEmSJEnSAFmvF6BIkiRJktbAZE6SJEmSBpDJnCRJkiQNIJM5SZIkSRpAJnOSJEmSNIBM5iRJkiRpAJnMSZIkSdIA+gOOS73skmVAjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "endog = macrodata['infl']\n",
    "endog.plot(figsize=(15, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constructing and estimating the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to formulate the econometric model that we want to use for forecasting. In this case, we will use an AR(1) model via the `SARIMAX` class in Statsmodels.\n",
    "\n",
    "After constructing the model, we need to estimate its parameters. This is done using the `fit` method. The `summary` method produces several convinient tables showing the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Statespace Model Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                   infl   No. Observations:                  203\n",
      "Model:               SARIMAX(1, 0, 0)   Log Likelihood                -472.714\n",
      "Date:                Sun, 07 Jul 2019   AIC                            951.427\n",
      "Time:                        20:47:31   BIC                            961.367\n",
      "Sample:                    03-31-1959   HQIC                           955.449\n",
      "                         - 09-30-2009                                         \n",
      "Covariance Type:                  opg                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "intercept      1.3962      0.254      5.488      0.000       0.898       1.895\n",
      "ar.L1          0.6441      0.039     16.482      0.000       0.568       0.721\n",
      "sigma2         6.1519      0.397     15.487      0.000       5.373       6.930\n",
      "===================================================================================\n",
      "Ljung-Box (Q):                       81.33   Jarque-Bera (JB):                68.45\n",
      "Prob(Q):                              0.00   Prob(JB):                         0.00\n",
      "Heteroskedasticity (H):               1.47   Skew:                            -0.22\n",
      "Prob(H) (two-sided):                  0.12   Kurtosis:                         5.81\n",
      "===================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "# Construct the model\n",
    "mod = sm.tsa.SARIMAX(endog, order=(1, 0, 0), trend='c')\n",
    "# Estimate the parameters\n",
    "res = mod.fit()\n",
    "\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forecasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Out-of-sample forecasts are produced using the `forecast` or `get_forecast` methods from the results object.\n",
    "\n",
    "The `forecast` method gives only point forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2009Q4    3.68921\n",
      "Freq: Q-DEC, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "# The default is to get a one-step-ahead forecast:\n",
    "print(res.forecast())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `get_forecast` method is more general, and also allows constructing confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "infl       mean   mean_se  mean_ci_lower  mean_ci_upper\n",
      "2009Q4  3.68921  2.480302      -0.390523       7.768943\n"
     ]
    }
   ],
   "source": [
    "# Here we construct a more complete results object.\n",
    "fcast_res1 = res.get_forecast()\n",
    "\n",
    "# Most results are collected in the `summary_frame` attribute.\n",
    "# Here we specify that we want a confidence level of 90%\n",
    "print(fcast_res1.summary_frame(alpha=0.10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default confidence level is 95%, but this can be controlled by setting the `alpha` parameter, where the confidence level is defined as $(1 - \\alpha) \\times 100\\%$. In the example above, we specified a confidence level of 90%, using `alpha=0.10`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specifying the number of forecasts\n",
    "\n",
    "Both of the functions `forecast` and `get_forecast` accept a single argument indicating how many forecasting steps are desired. One option for this argument is always to provide an integer describing the number of steps ahead you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2009Q4    3.689210\n",
      "2010Q1    3.772434\n",
      "Freq: Q-DEC, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(res.forecast(steps=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "infl        mean   mean_se  mean_ci_lower  mean_ci_upper\n",
      "2009Q4  3.689210  2.480302      -1.172092       8.550512\n",
      "2010Q1  3.772434  2.950274      -2.009996       9.554865\n"
     ]
    }
   ],
   "source": [
    "fcast_res2 = res.get_forecast(steps=2)\n",
    "# Note: since we didn't specify the alpha parameter, the\n",
    "# confidence level is at the default, 95%\n",
    "print(fcast_res2.summary_frame())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, **if your data included a Pandas index with a defined frequency** (see the section at the end on Indexes for more information), then you can alternatively specify the date through which you want forecasts to be produced:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2009Q4    3.689210\n",
      "2010Q1    3.772434\n",
      "2010Q2    3.826039\n",
      "Freq: Q-DEC, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(res.forecast('2010Q2'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "infl        mean   mean_se  mean_ci_lower  mean_ci_upper\n",
      "2009Q4  3.689210  2.480302      -1.172092       8.550512\n",
      "2010Q1  3.772434  2.950274      -2.009996       9.554865\n",
      "2010Q2  3.826039  3.124571      -2.298008       9.950087\n"
     ]
    }
   ],
   "source": [
    "fcast_res3 = res.get_forecast('2010Q2')\n",
    "print(fcast_res3.summary_frame())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the data, forecasts, and confidence intervals\n",
    "\n",
    "Often it is useful to plot the data, the forecasts, and the confidence intervals. There are many ways to do this, but here's one example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3YAAAEyCAYAAAC2+0LeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VNXWwOHfmUw6aSQhDUJCSyChJAQQUFFAAaXZUa+9e+1+em1XvfZ2vfaCimBXFCmKvQMikEZNAiGFkN4ndTIz+/sjxQBJSJlJMrDe58mT5Jwz52wmYTLr7LXX0pRSCCGEEEIIIYSwX7q+HoAQQgghhBBCiJ6RwE4IIYQQQggh7JwEdkIIIYQQQghh5ySwE0IIIYQQQgg7J4GdEEIIIYQQQtg5CeyEEEIIIYQQws5JYCeEEEIIIYQQdk4COyGEEEIIIYSwcxLYCSGEEEIIIYSd0/f1ADri5+enwsLC+noYQgghhBBCCNEn4uPji5VS/kc7rl8HdmFhYWzbtq2vhyGEEEIIIYQQfULTtKzOHCepmEIIIYQQQghh5ySwE0IIIYQQQgg7J4GdEEIIIYQQQtg5CeyEEEIIIYQQws5JYCeEEEIIIYQQdk4COyGEEEIIIYSwcxLYCSGEEEIIIYSdk8BOCCGEEEIIIeycBHZCCCGEEEIIYecksBNCCCGEEEIIO6fv6wEIIYQQQgghRHcppVo+N3+0/r6t4w7f1pn9bX0PYLFYOrWt9TnaOk8H49DaPbiVLgV2mqYtA+YDhUqp6KZtA4FPgTAgEzhfKVXWxmMvAx5o+vYxpdSKrlxbCCGEEEIIYd+UUlgslpbPh39tsVgwmUyHfAYO2d8c8DR/39Y1NK39WKit/Ud7zNGO6+y2zuxrrb6+HjqZZdnVGbvlwCvAe6223QP8pJR6StO0e5q+/1frBzUFfw8BcYAC4jVNW9tWACiEEEIIIYSwL2azGbPZjMlkwmQyUV9f3/K12WxuCdSadRQg6XSNcYxOp2s5RtM0NE1Dr9cf8n1nAyR71dHM3uG6FNgppX7XNC3ssM2LgFOavl4B/MphgR0wB/hBKVUKoGnaD8Bc4OOuXF8IIYQQQgjR+5RSLYGb2WzGaDRiNBqpq6ujoaEBs9mMpmktAZuDgwM6na7lQ6/X4+zsfMwHYn3JGmvsApRSeQBKqTxN0wa1cUwIcKDV9zlN246gadq1wLUAoaGhVhieEEIIIYQQorMaGhqor6+nvr6euro6jEYjDQ0NwN8zbc3Bm16vx9XVVQK2fqC3iqe09ZNuc15RKbUUWAoQFxfX+blHIYQQQgghRLcYjUZqa2spKyujvr7+kFm35tk20b9ZI7Ar0DQtqGm2LggobOOYHP5O1wQYTGPKphBCCCGEEKKXKaWor6+nurqaiooKGhoa0Ol0ODk54eHh0dfDE91gjcBuLXAZ8FTT5zVtHPMd8ISmaT5N358O3GuFawshhBBCCCE6wWKxUFdXR1VVFZWVlVgsFnQ6Hc7Ozri4uPT18EQPdbXdwcc0zrz5aZqWQ2Oly6eAzzRNuwrIBs5rOjYOuF4pdbVSqlTTtEeBrU2neqS5kIoQQgghhBDCNsxmM3V1dVRWVmIwGFBKodfrcXFxaak+KY4NXa2KeWE7u2a1cew24OpW3y8DlnVpdEIIIYQQQogusVgsLbNy1dXVLW0C3N3dpcjJMay3iqcIIYQQQgghbKyhoYGDBw9SX1+Ps7OzrJc7jkhgJ4QQQgghxDGgpqaGgwcP4uDgIAHdcUgCOyGEEEIIIeyYUory8nIKCgpwc3NDr5e3+Mcj+akLIYQQQghhp8xmM4WFhVRUVDBgwAApiHIck5+8EEII0Q07cip4eO0ulFJ9PRQhxHHKaDRy4MABqqqq8PT0lKDuOCc/fSGEEKIbvkjIYfmmTPIr6/p6KEKI41B1dTWZmZlYLBbc3d37ejiiH5DATgghhOiG1HwDABnF1X08EiHE8UQpRWlpKdnZ2bi4uEhjcdFCAjshhBCiG9IKGgO7zOKaPh6JEOJ4YTabyc3NpaioCE9PTymSIg4hvw1CCCFEFxUZ6impNgKQUVzVx6MRQhwP6uvrOXjwIBaLRVoZiDZJYCeEEEJ0UfNsnaZBhszYCSFszGAwkJubi5OTE25ubn09HNFPSWAnhBBCdFFK0/q62FAfMktkjZ0QwjaUUpSUlFBcXIy7uzsODg59PSTRj8kaOyGEEKKL0vIN+Lo7ETfUh+ySGswWaXkghLAuk8lEbm4uJSUleHh4SFAnjkoCOyGEEKKLUgoMjArwIMzPHaPZQm55bV8PSYgeu3fVdr7fld/XwxA0rqfLzs6mrq4ODw8PNE3r6yEJOyCBnRBCCNEFFotib4GBiEAPwnwbe0dJOqawd4WGOj7ecoDvdhX09VAEUFZWhlIKV1fXvh6KsCMS2AkhhBBdkFNWS43RTESgB+F+TYGd9LITdi4xuxyA/EqZfe4PlFLodPI2XXSN/MYIIYQQXZDaVBEzItCDAE9nXB0dpDKmsHsJ2WUA5FXU9fFIhBDdJYGdEEII0QWp+ZUAjApoXPcy1NdNUjGF3UvMapqxq6hDKSkGJIQ9ksBOCCGE6ILUgioG+7gywLmxY9Awf3dJxRR2rcFsYfvBcpz1OmqMZgz1pr4ekhCiGySwE0IIIbogLd9ARIBHy/dhvu5kl9ZgMlv6cFRCdF9KnoG6BgunRPgDjbN2Qgj7I4GdEEII0UlGk4X0oioiAlsFdn7umCyKg9LyQNip5vV1Z4wNAmSdnRD2SgI7IYQQopMyiqsxWdQhgV1zZcz9ko4p7FRCdhkBns7EhvoAUCCBnRB2SQI7IYQQopNSmgqnHDJj5ystD4R9S8guIzbUhwBPF0Bm7ISwVxLYCSGEEJ2UVmBAr9MY5jegZZvfACcGOOslsBN2qchQz4HSWmJCvXHS6/Ab4CS97ISwUz0O7DRNi9A0LanVR6WmabcddswpmqZVtDrmwZ5eVwghhOhtqfkGwv3ccdL//edT0zTC/NzIKJFedsL+JDatr2tOwwz0cpHiKULYKX1PT6CUSgUmAGia5gAcBL5s49A/lFLze3o9IYQQoq+kFhgYP9j7iO1hvu5sz6nogxEJ0TMJ2eU4OmhEh3gBEOjpSk6Z3KQQwh5ZOxVzFpCulMqy8nmFEEKIPlVVb+JAae0hrQ6ahfu5k1NWg9EkLQ+EfUnILmNMsBcujg4ABHo5k18pM3ZC2CNrB3ZLgI/b2TdV07RkTdO+0TQtqr0TaJp2raZp2zRN21ZUVGTl4QkhhBDds7fAABxaOKVZmK87FgUHZKZD2BGT2cL2nHJihvw9Cx3k5Up5TQN1DeY+HJkQojusFthpmuYELARWtrE7ARiqlBoPvAysbu88SqmlSqk4pVScv7+/tYYnhBBC9EhqfgeBnZ9UxhT2JyW/sTF57FCflm2BTZUxZZ2dEPbHmjN284AEpVTB4TuUUpVKqaqmr9cDjpqm+Vnx2kIIIYRNpRYYcHV0YIiP2xH7mnvZZUhgJ+xIQkvhlNYzdtLyQAh7Zc3A7kLaScPUNC1Q0zSt6evJTdctseK1hRBCCJtKzTcwKmAAOp12xD4fN0c8XfRklkhgJ+xHQlYZgzycCfF2bdkW0BTYScsDIexPj6tiAmia5gacBlzXatv1AEqpN4BzgRs0TTMBtcASpZSyxrWFEEKI3pBWYGBm5KA292maRrj/ADKLZY2dsB+JB8qJCfWm6d478HcqpszYCWF/rBLYKaVqAN/Dtr3R6utXgFescS0hhBCitxVX1VNcZWRUGxUxm4X7urE1s6wXRyVE9xVX1ZNVUsNFk0MP2e7urMfTRU+BBHZC2B1rV8UUQgghjjlpTYVTIgM92z0mzM+d3IpaqSYo7EJidjnAIYVTmgV6uciMnRB2SAI7IYQQ4ihSm1odjAoc0O4x4X7uKAXZpZKOKfq/hOwy9DqNsU2NyVsL9HKVXnZC2CEJ7IQQQoijSM034OPmiP8A53aPCfOVypjCfiRmlxEV7NnSmLy1IE8XaXcghB2SwE4IIYQ4itQCAxGBHocUmTic9LIT9sJktpB8oIKY0CPTMKExFbOoqp4Gs6WXRyaE6AkJ7IQQQogOWCyKtHwDER0UTgHwcnVkoLuTtDwQ/V5KvoHaBjMxrfrXtRbo5YJSUGio7+WRCSF6QgI7IYQQogMHy2upNpqJ6KBwSrMwXzdJxRT9XmJLY/L2Z+wA8iukl50Q9kQCOyGEEKIDqU0VMSM6KJzSLMzPXXrZiX4vMbscfw9nBvu4trk/qCWwkxk7IeyJBHZCCCFEB1oqYh4lFRMg3Ned/Mo6ao3S8kD0XwnZZcQM8W53zWiQZ2PAlyczdkLYFQnshBBCiA6k5hsI8XbFw8XxqMe2FFCRdXaii+5amcyDa3ba/DolVfVkltS02b+umaerHhdHnVTGFMLOSGAnhBBCdCCtqSJmZ4RLZUzRDelFVayMz+GTLQeoqGmw6bVaGpO3s74OQNM0grxcyZNedkLYFQnshBBCiHY0mC2kF1V1Kg0T/p6xy5AZO9EF723KRNPAaLawfmeeTa+VeKD9xuStBXq6UCAzdkLYFQnshBBCiHZkFFfTYFZEdnLGboCzHr8BzjJjJzqtsq6Bz+NzOCsmhGF+7qxOPGjT6yVklTM6yBNXpyMbk7cW6OVCngR2QtgVCeyEEEKIdjRXxOzsjB3AMD93aXkgOm3lthyqjWaumBbOogkh/JVRSm65bYqWmMwWknPKiW2nf11rgV4uFFTWYbEom4xFCGF9EtgJIYQQ7UjNN+Cg0xg+yL3TjwnzcyNDWh6ITjBbFCs2ZRI31Iexg71YNCEYgLXJuTa5XmqBgRqjucPCKc2CvFwwWRQl1UabjEUIYX0S2AkhhBDtSC0wEO7njrO+47S11sL83CmuqsdQZ9siGML+/ZxSSHZpDVdMDwcaf3cmDPG2WTpmQicKpzQL9GzuZSfpmELYCwnshBBCiHak5huI6EIaJjT2sgPIKpFZO9Gx5ZsyCPJyYU5UQMu2xROCSck3tKQBW1Nidhl+A5zabUzeWmBTk3LpZSeE/ZDATgghhGhDjdFEdmlNp1sdNGupjCnr7EQHUvMNbNxXwiVTh6J3+Pvt2PzxwTjoNFYnWX/WLjG7nJhQn3Ybk7fWHNjlS8sDIeyGBHZCCCFEG9IKqoCuFU4BCPOVXnbi6JZvysRZr+PCSaGHbPcb4MyJI/xYm5Rr1cIlpdVGMoqrO5WGCeDn7oxep0kqphB2RAI7IYQQog1pTalwnW110MzVyYFATxfpZSfaVV5j5MvExhYHPu5OR+xfHBPMwfJatmWVWe2aidmN5+pMRUwAnU4jwNNFAjsh7IgEdkIIIUQbUvINuDjqGDLQrcuPDfNzkxk70a6PtxygrsHC5dPD2tx/+phAXB0drJqOmZhdjoNOY+zgjhuTtya97ISwLxLYCSGEEG1IKzAwKsADB93R1yMdLtzPnUwpniLaYDJbeP/PTKYN9yUy0LPNY9yd9Zw2JoD1O/IwmixWuW5Cdhmjgzxwc9J3+jGBXi6yxk4IOyKBnRBCCNGG1KbArjvCfN0prTZSUSstD8Shvt9dQG5FHZdPC+vwuMUxwZTXNPBbWlGPr2m2KJIPlHd6fV2zoKZUTKWkSbkQ9kACOyGEEOIwpdVGigz1XV5f16y5MqakY4rDLd+YyZCBrswaHdDhcSeN9Gegu5NV0jFT8w1UG81dDuwCvVyobTBTWWvq8RiEELZntcBO07RMTdN2aJqWpGnatjb2a5qmvaRp2j5N07ZrmhZrrWsLIYQQ1tTcQ6y7M3bhzYGdFFARrew8WMGWzFIumxp21BRfRwcdZ44N4sfdBT1udp94oLlwStcDO4C8SullJ4Q9sPaM3alKqQlKqbg29s0DRjZ9XAu8buVrCyGEEFaRml8J0OUeds1CB7qhadLLThzq3Y2ZuDk5cF7ckE4dvzgmmHqThe92FfTouglZ5fi6OzFk4NEbk7cW1NzLTgqoCGEXejMVcxHwnmq0GfDWNC2oF68vhBBCdEpqQRXebo4M8nDu1uNdHB0I9nKVwE60KK6qZ11yLudOHIyXq2OnHhMb6sOQga6s6WE6ZmJ2Wacbk7cW6NUYCEpgJ4R9sGZgp4DvNU2L1zTt2jb2hwAHWn2f07RNCCGE6FdS8ysZFeDR5TfCrYX7ucsaO9Hio7+yMZotXHaUoimtaZrGovEhbNxXTKGhe8FVWbWR/cXVxA7tXP+61gZ5OKNpSMsDIeyENQO76UqpWBpTLv+padrJh+1v66/jEWWWNE27VtO0bZqmbSsq6nklKCGEEKIrlFKkFVR1u3BKszA/NzKKq6WioMBosvDB5ixmjPJnuP+ALj12cUwwFgXrkvO6de2kA+VA19fXQeM6P78BzjJjJ4SdsFpgp5TKbfpcCHwJTD7skBygdVL5YCC3jfMsVUrFKaXi/P39rTU8IYQQolMOltdSVW/qduGUZmG+7lTWmSirkZYHx7tvduZRaKjninYakndkxCAPooI9u52OmZBdhoNOY1wXGpO3FiS97ISwG1YJ7DRNc9c0zaP5a+B0YOdhh60FLm2qjnkCUKGU6t7tJyGEEMJG0goaK2L2dMauuTKmrLMTyzZmMszPnZNHdu+G9eIJIWzPqWB/UVWXH5uQXUZkYNcak7cW2NTLTgjR/1lrxi4A2KBpWjKwBfhaKfWtpmnXa5p2fdMx64H9wD7gLeBGK11bCCGEsJqUplYHI3s6Yye97ASNhUuSD5Rz+fQwdEdpcdCehROC0TRYnXREolOHzBZFUnbXG5O3FujlQl6FtDsQwh507/bNYZRS+4HxbWx/o9XXCvinNa4nhBBC2EpavoFgL5dOVy5szxAfN3Sa9LKztQazhTd/S+fUyEFEBXcv3dCW3t2YiYeznrNjB3f7HAGeLkwb7suapIPcPntkp4v67C1sakzejcIpzQK9XKisM1FjNHV71k8I0Tt6s92BEEII0e+l5BsY1cM0TAAnvY7BPm6SimlDDWYLt36SyHPfp3Hte/E9buRtbQWVdazfkcf5k4YwwLlnQdGiCSFkldS0FEPpjISsxmNjhnR/xk562QlhPySwE0IIIZo0mC3sL6rudmPyw4X5ucuMnY2YzBZu+zSJ9TvyuWhKKHkVtTyybndfD+sQH2zOwqwUl00N6/G55kYH4qTXsaYL6ZgJ2WUMdHdiqK9bt68b4CmBnRD2QgI7IYQQoklWSTVGs4WIHq6vaxbu60ZmcY20PLAyk9nCHZ8l8/X2PO4/YzRPnDWWG08Zwcr4HL7fld/XwwOgrsHMR39lMysygNAeBFbNPF0cmT16EF9tz8VktnTqMQnZZcSGeveoH2NQU5Ny6WUnRP8ngZ0QQgjRpLlwSk9bHTQL83Onqt5EcZXRKucTjQVB/m9lMmuTc7lnXiTXnDwMgFtmjSQq2JN7V+2gyFDfx6OEdcm5lFQbubIbLQ7as2hCCMVVRjbsKz7qseU1RvYXVRPTg8Ip0FgVE5CWB0LYAQnsjlEWi2JrZinPf5/KgdKavh6OEEJYncVi/VmwtHwDOg1GDOpaE+n2tFTGlHRMqzBbFHetTGZ1Ui53zYng+hnDW/Y56XW8cMEEDPUm7l21vU9nSZVSvLsxk4gAD6YO97XaeU+J8MfTRd+pdMzEprV4MaHdL5wC4OrkgLebo6RiCmEHJLA7hiilSMgu45F1u5n21M+c98afvPTzPi5dtoWSqr6/eymEENbyS0ohk5/4ie05nS8k0Rkp+QbC/NxxcXSwyvnCfZt62RVJYNdTFoviX19sZ1XiQe48bRT/PHXEEceMDPDg7jkR/LinkM+2HeiDUTbamlnG7rxKLp8e1qM0yMM56x04c1wQ3+3Kp8Zo6vDYxKwydBqMH9yzwA4aZ+0kFVOI/k8COzunlGJHTgVPrt/DiU//wtmvbeKDzVlEh3jywgUTeO/KyeSW13Llim1H/SMghBD2oMZo4oHVOymuquexr/ZYdWYmrcDQ48bkrQ32cUWv08iQGbsesVgU96zazufxOdw2eyQ3zxrZ7rFXTg9n6jBfHlm3m+ySvslYeXdjBt5ujiyeEGL1cy+aEEKN0cwPuws6PC4hu5zIQE/ce1iNExpbHuRXSi87Ifo7CezskFKKPXmVPPtdCqc+9ysLXtnAOxsyGBkwgOfOG8/WB2bz9mWTWBwTwsmj/Hnpwhh25JRz80eJnV5wLYQQ/dWLP+3lYHkt500czJbMUn7cU2iV89YYTWSV1lhtfR2A3kFH6EC3ftuk/EBpDXd8lkRxP87qsFgU96/ewWfbcrhl5ghumz2qw+N1Oo3nzh+PTtO4c2USZhuk7HYkp6yG73bls2RSKK5O1pn5bW1y2ECCvVxYnXiw3WPMFkXSgfIe9a9rLcjLhfyK/vs7IoQ9U0q1fACYzWYaGhpoaGjAaDRiMnV+YkY6TdqRfYUG1iXn8dX2XNKLqtFpMG24H9fPGM6cqEB83J3afNycqED+syiaf6/eyb/X7OSJs8ZaNTXkaGP+YHM2t80eibdb2+MTQojOSs038M4fGZwfN5gnzhpLfHYZT32zh1Mj/NE79Oxe5b7CKpTCqjN20LjOrj/2sisy1HPJO3+RWVLD5LCBLJkc2tdDOoJSin+v2cnHWw5w06kjuP20joO6ZiHerjy8MIo7Vybz1h/7D1mLZ2vvb85C0zQumTrUJufX6TQWTAjm7T8yKKmqx3eA8xHH7CusoqreRGwPC6c0C/R0pbiqHqPJgpNe5gSOd0opKisrMRqNGI1G6uvrMRqN+Pj4EBAQQH19PZs2bWoJTkwmEyaTiaioKCIjI6msrGTlypUt25s/Tj31VGJjY8nLy+O1116joaEBi8WC2WzGbDZz8cUXM2nSJPbt28czzzxzyD6LxcKtt97KpEmTSEhI4LHHHmvZp5TCYrHw+OOPExsby2+//cYjjzyCxWJp+VBK8cYbbxAdHc3atWt59NFHW7Y3f/7iiy8YMWIEH3zwAU8++SRAyz6lFD///DMhISG8+uqrPP/884cEbEopkpKS8PHx4cknn+SVV1454nnNzMzE0dGRf//736xYsaJl+yWXXNLpn40Edv1cZnE1X23P5avteaTkG9C0xrt1l08PZ150IH5tvKC35ZIThpJfUcurv6QT4Oly1Due1rAlo5SrV2ylss5EvcnMk2ePs/k1hRDHLotF8cDqHXi46Lln3mj0DjrumRvJte/H88nWA/zjhJ69kbZ2RcxmYb7u/JleglKq126qHU1FbQOXLttCQWU9A5z1bMsq63eBnVKKB9fs4sO/srnhlOHcefqoLj1/Z8eG8MPuAp7/Po0Zo/wZHeRpw9E2qjGa+GTLAeZEBRDi7Wqz6yyeEMKbv+3n6x15XNpGj7yE7DKAHlfEbBbo1fheo6CyjiEDe966QdheRkYGBoOBmpoaampqqK2txc/PjylTpgDw2muvUVpaSm1tbcsxkydP5qqrrgJg9uzZVFVVHRK4XXzxxTz88MPU19czZsyYI6550003ce+991JdXc0//vGPI/b/61//IjIykvLych588MEj9vv4+BAbG4vBYOCLL75Ar9ej1+vR6XTodDpmz54NQF1dHWlpaTg4OKDT6XBwcMDBwQGjsbH6sKZp6HQ6HB0dW/ZpmoaTU+MEg5ubG8OGDUPTtJZjdTodbm6Nv9sBAQGcdNJJ6HS6lmM0TcPdvXHN9LBhw1i8eHHL/uZrNj9+3LhxXHHFFYc8FsDFpbHC7LRp09Dr9Yc8tnkczc/9oEGDWvZHRkby/vvvd+rnLoFdP1VR28ADq3eyLrmx8tXEoT48tGAMZ4wNamkW2lX/d3oE+RX1vPDjXgI9XWz6R/zr7Xnc/lkSg31cmT0mgI+3HODciUOYONQ6f2SEEMefzxNy2JpZxjPnjGNgU4bCaWMCmBw2kBd+TGNxTAgDerCeKC3fgLNex9CmgifWEu7nRm2DmYLKegK9uvf6bU21RjPXrNjGvkIDb182iff/zCIhq6yvh3UIpRQPr93F+5uzuO7kYdw9J6LLQbGmaTxx9lhO/9/v3P5pEmtumo6z3vqpka2tTsyloraBK6aH2/Q6o4M8iQjwYHXiwbYDu6zGxuRhVuifBxDY1MsuXwK7XlNcXExubi61tbWUl5dTUVGBi4sLCxYsAODpp58mJSWFiooKKioqKC8vZ+zYsSxfvhyAiy++mKysrEPOedppp7UEdm+99RaVlZW4urri5uaGq6sr4eF//96OHj0aaAxGHB0dcXJyanmss7MzDz30EE5OTjg7O+Pk5ISTkxMjRzauffX09GTNmjU4OTnh6OjYEqD5+DS+BwwJCWHHjh2H7GsO0gBGjRrF7t27231uoqOj+fXXX9vdHxMTw+eff97u/kmTJjFp0qR290+ZMqXl39qWadOmMW3atHb3n3TSSZx00knt7p8xYwYzZsxod//MmTOZOXNmy/c1NZ1fKyyBXT+0NbOU2z5JIr+yjptnjmDJ5FCr3PnTNI2nzhlLUVU996/eib+HM7NGB1hhxId6Z0MGj329m4mhPrx9WRyODjo27StpDFRvmt7jdCkhRM/VGE046DQcdTp0uv4xi9SRsmojT67fQ9xQH86dOLhlu6Zp3HfmaBa/upGlv6Vzx+kR3b5GaoGBkQEDcLDy89Hc8iCjuLrPA7sGs4WbPkpga1YpLy2JYcYof3bnVvLjngJKq40tAXNfUkrxyFe7WfFnFlefGM498yK7PdM50N2JZ84dy5XLt/H8D2ncO2+0lUf7N6UUyzdlEB3iSVwv3MRcFBPMM9+mkl1Sc0QD9ITsMmKG9KwxeWtBTb+30vKg99x4441s2LDhkG2jRo1qCezS09PJycnBy8uLsLAwvLy8WoIxgEcffRSTyYSrqyvu7u64ubm1BFYA8fHxLYFUW15++eV292maxrXXXtvufr1eT1xcXLv7HRwcGDhwYLv7RfdJYNePmMwWXvp5H6/8vJfBPm58fv1Uq6VRNHN00PH6xbEsWbqZf36UwMfXnGC1a1gsisfX7+GdDRnMjQrkhSUTWkqGP7RgDDd8mMCKP7O46kTb3sk8FhVW1nHvqh08uGCM1We2AILVAAAgAElEQVQTxPHF3FQy/vP4nJZtep2Go4MORwcNJ72u6evG7x0ddK22aTjpHXBy0HBz0nPzzBGMtHLaYnue+iYFQ52Jx86KPiIQnTDEm/njgnjrjwwuPmFot7MaUvMNnDTS3xrDPUSY79+97KzZ06yrLBbF3Z9v56eUQh5bHM2C8cEAxIU1/g2IzyrjtDHWv9nXFUopHv96D+9uzOSK6WHcf+boHgcnMyMDuHByKEt/38+syAAmh9vmDeXGfSWkFVTx3HnjeyXlduH4xsBuTdLBQ6qEVtQ0kF5Uzdmxgzt4dNc0/5+SwK5jFbUN7C+qaryJ4+nCtBF+3T7XjTfeyPnnn4+fnx9eXl54eXnh7f13MZylS5d2+PhZs2Z1uL+joE7YLwns+okDpTXc9mkS8VllnB0bwn8WRuHh4miTa7k761l2+STOeX0TV63Yxhc3TCPcr2fBQl2DmTs/S+brHXlcPi2Mf88fc8hd77nRgZwS4c/z36dy5tigPr9rbW8+2pLNTymF1DaY+fDqKf1mnY6wLyazhTtXJrMmKZeLpjRmAjSYLU0fCqPJgtFsocHUalvLfgsNJkVlbQMNZgvZpTVsyyxl9U3TGeRh2//P2zJL+XTbAa47eRiRgW2vk7p7TiTf7crn+e/TePrcrq/nLas2UmioJyLQOo3JWwv2dsXJQdenlTGbZ8G+TDzI/50+6pD1iGNDvHB00Po8sFNK8eQ3Kby9IYPLp4Xx4PwxVnute+DM0WxKL+aOz5L45taTbPL3dfmmDPwGOLFgfJDVz92WwT5uTA4byOqkg9w0c0TLc5V4oGl93RDrVMQE8HTR4+bkIL3sgHqTmeySGvYXV7O/qJqM4qqmz9WUVBtbjls4PrhHgd0pp5xCXV0dzs6dq6UgBEhg1y+sTc7l/lU7AHhxyQQW2aDvzeH8PZxZceVkznl9E5cu+4tVN0zH36N7Lx4VNQ1c8/42tmSUcv8Zo7n6pPAj/hhrmsZ/FkZx+v9+59Gvd/PqRbHW+GccFywWxRcJOXi46NmUXsKXiQeteidWHB9MZgt3fJbM2uRc7poT0WZz567YebCC8974k2tWbOOTa6fapKw7NKYOPrB6J8FeLtzSQe+yUF83Lp0axrsbM7jyxHAiuljZMrXANoVTABx0GqG+bn1aGfPln/exfFMmV50YfsTP3sXRgahgL+KzSvtodI1B3dPfprL09/1ccsJQHlpgvaAOGm9oPn/+eM57408e/Wo3z5w73mrnrjGaWPr7fn5KKeTmmSNtvo6vtUUxwdz/5U525VYSHeIFNPav02kw3oqBnaZpx1UvO4tFkV9ZR0ZxNfuLqloFcdXklNXQuoOG3wBnhvm7c9qYAML93BnmP4BwP3dCZS2i6AMS2LVjb4GBt/7Yz6zRAcyMHISjDdaFVdWbeHDNTlYlHCQ21JsXl8T06qLkcD93ll0+iQuXbuaK5Vv45NqpXS48kFNWw+XvbiW7pIaXLoxhYVNqT1uG+rrzz1NH8PwPaVwQV8TJo6yf8nQs+iujlAOltTx//nje35zFY1/v4dSIQe22txDicCazhds/S2Zdci53z43gxlN6FtQBRId48eKSCVz3QTx3rkzilQtjbbJW792NGaTkG1h6ycSjNlq+eeYIVm47wJPf7GH5FZO7dJ20psCuvRnBngrzdSezj5qUv/9nJs//kMY5sYO5/4y2Uxvjhvrw3uasPilnr5Tiue9TeeO3dC6eEsp/FkbZJCth4tCBXD9jOK/9ms7s0QGcHhXYo/OZLYqV2w7w/A9pFBrqOWNsIFef1LtLDc4cG8TDa3exOvFgS2CXmF1GhJUak7fW2Mvu2J+xSz5QzmXvbqG8pqFlm6ujA+F+7owb7MXiCcEtwVu4vzueNsquEqI7JLBrQ15FLZcu20JeRR2fbcthkIcz58UN5oK40CMWKHdX0oFybv0kkQOlNdwyayS3zBzRJ0VFJgzx5tWLY7jmvXhu/DCBd5qKnXTGrtwKrnh3K7UNZlZcOblTa0eumzGMLxMP8uCanXx728kta/BE+z6Pz2GAs5550UGMCfZk/ksbeOqblG6lm4njT+ug7p55kVbt53V6VCD3zRvN4+v38JxvKnfPjbTauQFyy2t54ce9zB49qFNvwr3dnLhp5gieWJ/Cxn3FTO9CGlRKvgFPFz0BnrZJewr3c+P3vUVYLKpXi9WsSTrIg2t3MXt0AE+fM7bda08c6sPbGzLYlVth9bXdR7N+Rz6v/pLOhZOH8OiiI9dQWtNts0fxa2oR967aQexQn063DGpNKcWvaUU8tT6F1AIDsaHevP6PWCYO7f1iEN5uTswYNYi1ybnce8ZoNCApu5yFE9q/ydpdAZ4ubE4vsfp5+xOlGmsF6HU6HlsczbCm4C3Q00WWQAi7ICsnD1NZ18AV727FUGdi7U3TWXrJRKJDvHj913ROfvYX/vH2X3y1PRejydKt85stild/2ce5r2+iwWThk2uncsdpo/q0UuTMyAAeXxzN72lF3PPFDpRSR33MH3uLuODNzTjoNL64YVqnCwI46x14dFE0mSU1vPFbek+Hfsyrqjexfkce88cF4erkQGSgJ1edFM6n2w6wJaPv0qaEfTCZLdz2aRLrknO518pBXbOrTwrnwslDeO3XdFZuO2DVc/9n3S4sSvHQgqhOP+bSqWGEeLvyxPo9WCxHfy1rlpZvIDLQ02Zv3sL83DGaLORW9F4q2y+phdz5WTKTwwbyykUxHf6daW5FE98HbQ++352P3wAnHl/cfuBpLU56Hf+7YAKGOhP3rurc37vWduVWcMk7W7ji3a3Umcy8fnEsX9wwrU+CumaLY4IpNNSzeX8J+4qqMNSbbBKcB3m5UGCox9yF/1f25o+9xWzJKOWWWSP4xwlDmTbCjyAvVwnqhN2QwK4Vo8nCde/Fs6+witf/Ecu4wd6cHhXIsssnseFfM7l99igyiqu56aNEpj75E0+s30N6UVWnz59XUcs/3v6LZ79LZU5UIN/cerLNqnN11ZLJodw2eyRfJOTw3PepHR67KiGHK97dymAfV768cXqX16ScONKPBeODee3X9D4tJmAP1u/Io7bBzHlxf6+pu3XWSAb7uHLflzu6fYOhr+RX1GEy29eY7VWD2cKtnyTx1fY87jsjkutsENRB49qbRxZFM32EL/d9uYPN+61zR/+nPQV8t6uAW2eN6lKKuoujA3fPjWBXbiWrkw526jFKKVILDIyyQeGUZuHNlTGLO9+PqCfis0q54YN4IgI9eOuyuKNmRwzydGHIQNdeD+yUUmxKL2HqcL9em8mMCPTgrjkR/LC7gJWtqsN2JK+iljs/S2b+yxvYmVvBQwvG8MPtM5g3NqjP3/TPHh3AAGc9qxMPtvQjjA213vq6ZoFerpgtipKqequfuz9oTgkO8XZlySTb9fkVwpYksGtisSju+jyZP/eX8My5444oeR3s7cqts0fy+92nsvyKSUwKG8iyDRnM+u9vnP/mn3yZmENdg7nd83+7M595L/5Bck45z5w7jlcuisHLrX/lZd86ayQXTh7Cq7+k8/6fmUfsV6pxtvGOz5KZMmwgn10/tdvVLR84czRODjoeXLury3dMjyefx+cwzM+d2FZ3X92c9Dy6OJp9hVUs/b3/z3ruL6rilZ/3Mu/FPzjhyZ/41xc7+npIfSq9qIpvduR1+HrRU41BXSJf78jj/jNGc+3Jtgnqmjk66Hjt4omEDnTjuvfj2d+FG15tqTWaeWjtLkYOGtCt9igLxgUTHeLJc9+ldup5zquow1BnIsJG6+sAwv2betn1wjq7PXmVXPHuVoK8XFlx5eROrwGaGOrDtqyyXn1NTi+qoshQz/RebgNx1YnhTAkfyCPrdnOgtP1g21DXwLPfpXDKs7+ybnsu1548jN/uOpUrpof3+lrE9rg4OjAnKpBvd+bz5/4SfNwce1zpui2BTS0PjtXKmN/tKmB7TgW3zh7Zb362QnSV/OY2efq7FNYkNVaL66jioINO45SIQbxxyUQ23TuTu+dGUFBZx+2fJjPliZ94eO0uUvIrW46vNZq578sdXP9BPEN83Pjq5hM5P25In9/ha4umaTy6KJpZkYN4cO0uvt2Z37LP1FSZ7tnvUlk8IZh3L+/8m4W2BHi6cOfpo/g9rYj1O/KP/oDjUFZJNVsySjln4uAjfl9OjRjEmWODePnnff1y1jO9qIqXf9rL3Bd+Z+Z/f+O579NwddQxNyqQLxJyWJuc29dD7FX7Cqt4qen5mPXf37jhwwROefZXPt6SbfUZzAazhVs+TmT9jnweOHM015w8zKrnb4+XqyPLLp+Eg07jqhXbKK8xHv1B7Xj5573klNXy2OLobr3B0uk07jtjNLkVdSzflHnU41PzGwunRNiwJ1+AhwsujrZveZBVUs2ly7bg5qTn/asmd2kN2cSwgRQZ6skp67100Y37Gmd4pw3vfln47tDpNP57fmNlzDtXJh+RXthgtvD+n5mc8uyvvPpLOvOiA/n5zhncO280Xq7966YsNKZjGupNrEvOJSbUxybvMZqblB+LgZ3Zonj+h1SG+btzdoztK5MLYSv9unhKRW3D0Q+yghWbMnnzt/1cPCWUG0/p/J3tQR4u3HjKCK4/eTib95fw8dYDfPRXNss3ZRIT6s2CccF8+FcW6UXVXDdjGHeeFtHv7wLpHXS8fFEMF731F7d+ksiHV09hTLAnt3ycyI97CrnhlOHcdXqEVVJmLjlhKCu35fDIV7uYEeHf5Yqcx7ov4nPQaXB2bNt/ZB5cMIbf04r495qdvHfl5D6/WbCv0MDX2/NZvyOvpXT8xKE+/Hv+GOZFBxLs7YrJbOH8N//k/i93EBvqzWCfY7cc9OHPh6Y1Vh58aMEYhvq68crP+7h31Q6W/r6fO08fxRnRQT3+f9Uc1H2zszGou/qk3gnqmg31dWfpJRO56K2/uO79eN6/akqXX/P2FhhY+vt+zokdzJRh3Z/FmTbcj5mRg3j1l31cEDekwyqyzb+vtgzsdDqtsTKmDQO7wso6LnlnCw1mCx9dN7XL/78mhv69zq63KjRvSi9msI+r1QqTdcVgHzceXhjF/61M5p0N+7n25OEopfhhdwFPfZPC/uJqpoQP5N0zRzNusPVTG61p2nA//D2cKTLUW7V/XWvNGTr5vbhOtLesS84lraDqqGtRhejv+vU76ezSGh5Zt5t7z4i0SbsBaEyRfHhdY8WwRxZFd+vNsU6nMW2EH9NG+FFabWRVQg6fbD3AI1/tZpCHMx9cNYUTR/bu3ciecHPS885lcZz7xp9ctWIbQ33d2HmwgkcXRXHJ1DCrXUfvoOPxs6I5+/VN/O+HNP49f4zVzm3vGnvXHeTEkf4Eebm2eUyApwt3z43g32t2sTY5t1f6Hx5ub4GBr3fksX5HHmkFVYcEL/Oij2xEr3fQ8eKSGOa9+Ae3f5rEx9eccEz9EW3r+Zg0dCAPLxjDvLFBBHj+/XycGjGIH/cU8tx3qdz0USJRwencPTeSk0f6det1qMFs4eaPEvl2Vz7/nj+mWymM1hAXNpBnzh3HbZ8mcd+XO3j23HGd/vcopXhg9U7cnfXcd0bPK2zeOy+SOS/8zks/7+2wAEtqvoFATxebp8eH+bqTVmiwybkrahq4dNkWiqvq+eiaExjZjSA1ItCDAc56tmWVsrgXZi3MFsWf6SXMje5Z24GeOCc2hB925/Pcd2kM8nDhoy3ZbMkoZbi/O29fGses0YP6/KZZZzjoNBaMC2bZxgxih9qmqulANyecHHTkVx5ba+wazBae/yGN0UGenBHdO83lhbCVfh3Y+Q5wYtnGDLbnlPPKRbHdXs/VnvisUm79JJHxg715+cIYHKwwCzXQ3YmrTxrGVSeGk5JvINjbtV+mbRyN7wBnVlwxmbNf30RagYE3/jGxxz1/2hIT6sOSSaEs35TJObGDGRNsuzUu9uTP/SUcLK/lX/M6fnN70ZShfJ5wkEe/2s0powb1yrrNtAIDX29vDF72FnYcvLRlyEA3Hl0cxe2fJvPar+kdNp3u75RSpBVUtQRz+5qfj7CB/GdhFHOjA9t9PjRN47QxjX0y1yQd5Pkf0rhs2RZOGDaQu+dGHrKu8miMJgs3f5zAd7sKeHD+GK7so6Cu2eKYEPYXV/PST3sZ5u/e6b55qxIO8ldGKU+ePRbfbpShP9zIAA8umDSEDzZncfm0MIb6tr3uKDXf0OWG5t0R5ufOTykFmMwWq97QqDWauXLFVvYXVbPs8klM6OaMjYNOIybUm/iscquNrSO7cyuprDN1qS2FtWmaxhNnjWXOC39w26dJ+A1w4rHF0SyZNMTubjpdfVI4FqWIC7NNYKfTaQzydD7mZuxWbsshu7SGZZfH9WorEiFsoceBnaZpQ4D3gEDAAixVSr142DGnAGuAjKZNq5RSjxzt3MFerjx0YQz/+mI781/+g5cujLFaHn56URVXrdhGkJcL71wWh6uTdfupaZrG6CD7DlJCfd1Ye9N06k0WmyzEbvavuRF8vyufB1bv4PPrp8kLK7By2wE8XPScPiagw+McdBpPnBXNwlc28tS3KTx59libjemPvUX8Z93uLgUv7TkrZjC/phbx4k97mT7Cr6XUuj1orqC4fnseX+/II72oGk2DyWEDuXRRFHOjAhnUhefDQadxduxg5o8L5uMt2bz88z7Ofm0Tp40J4P9OjzhqwGE0WfjnRwn8sLuAhxaM4YrpfRvUNbt99kgyi6t55ttUwnzdOWNsx3fCy2uMPL5+D7Gh3lwQN8SK4xjFmqRcnvk2lVcvjj1iv8lsYV9RVa9kVYT7udFgVuSW11kt9dBosnDDh/EkZpfx6kWxPf53xIb68PLPezHUNeBh48bLm9KLAZjag5Rba/Ad4MzSSyeyNaOUi08YarfLAoK9XXl4Yedbg3RHkJfLMbXGrq7BzEs/7SU21JtTIwb19XCE6DFrvHqZgDuVUgmapnkA8Zqm/aCU2n3YcX8opeZ39eQLxgczOsiD696P5x9v/8XdcyO57uRhPUqNKDTUcdmyLThoGiuunGyVO8PHqmDvttMArcnbzYl7zxjN/61M5rNtB1gy+fguM1xZ18C3u/I5J3Zwpxq4RwV7cdWJ4U3rkkKIC7N+C41PtmRz/+qdhPu580g3gpe2PLo4mvisMm77NJH1t5xk8zeR1tBgtnDd+/H8nFKIToPJ4QO5fFoYc6IDGeTRs+fDSa/jsmlhnDtxMO9uzODN3/Yz98XfOWtCCLef1nbJ/9ZB3cMLxnB5PwnqoPHm1jPnjiOnrIbbP00i2Nu1w5mkp79NpaK2gces3MtskKcL15w0jBd/2stV2WVHzIRmltRgNFm63LalO8J8/66Maa3A7rGvd/NrahFPnj2WeUcJnjsjLswHi4KkA+VHVIe2to3pJYwYNKDHryXWEBvq06VZ8uNVoJcr23N6Z0a3N3ywOYv8yjqev2C8XaTcCnE0Pc4zUErlKaUSmr42AHsAqybnjxjkwZqbTmTe2CCe+iaF696Pp7Kue4VVqutNXLl8KyVVRpZdPqnd1BzRu86JDWFy2ECe+jaF0uruV9M7Fny9PY+6BgvndWHW4rbZIwnxduX+L3fSYMUqi0opnv0uhXtW7WD6CD++vHEal04Ns8obMU8XR164YAIHy2p5aM0uK4zWtpRSPLhmFz+nFHLHaaP4677ZfHLtVC6ZGtbjoK41d2c9N80cyR//OpVrTx7G1zvymPnfX3l47S6KDH+vbTGaLNz4YWNQ95+FUf0qqGvm4ujA0kvj8Pdw5uoV2zhY3nYKV3xWGR9vyeaKaWE2Sce+9uRh+Hs488TXe44o5Z/WVDglshdSMZszH6xVQGVXbgXvN6WZXmilG2IThnijabZvVG40WdiaUdrrbQ5EzwR5uZBfUXdMtCmqqjfx+q/pnDjCr9ersgphK1ZNINc0LQyIAf5qY/dUTdOSNU37RtO0dnMFNE27VtO0bZqmbSsqKmrZPsBZzysXxvDg/DH8nFLIwpc3sCevsr3TtKnB3PhGaE+egVcvjmG8jSpHia7TNI1HF0dTVWfiqW/29PVw+tTn8TmMGDSA8YO9Ov0YNyc9jyyKIrXAwFt/7LfKOOpNZm77NIlXf0nnwslDeOeyOKvPqsWFDeTmmSNZlXiQNZ1sJt1X3t2Yycdbsrl+xnBumTUSfw/bzvR7uzlx77zR/HbXqZw7cQjvb85ixrO/8N/vUympqufGD+P5cU8BjyyK4rJpYTYdS0/4DXBm2eWTqG8wc9XyrRgOuynX3Eol0NOF204bZZMxuDvruX32KLZllfHdroJD9qXkG9BpMGKQ7ZqTN/P3cMbdyYEMKwR2SikeWbcbHzcnbp9tvefNw8WRiAAPmwd2SQfKqW0wM1XeUNuVAE8X6k0Wymt6p2q5Lb27IYOSaiP/Nyeir4cihNVYLbDTNG0A8AVwm1Lq8IgrARiqlBoPvAysbu88SqmlSqk4pVScv/+haSCapnHlieF8cu0J1DaYOeu1jaxKyOnU+JRS3P/lDn5LK+KxxdHMjOx47ZLofRGBHlx1YjifbcthW2ZpXw+nT+wvqiI+q4zz2uhddzSzRgcwNyqQF3/cS3ZJ+w13O6O8xsglb29hTVIud8+N4ImzxtqsMu3NM0cwcagPD3y5s8NGwX3pl9RCHvt6N6ePCeDuXn4TEOjlwpNnj+XHO2Ywa3QAL/+8j8lP/MSPewp5dHE0l1qxUq2tjArw4NWLY9lbWMUtHyce0rtv+aZM9uRV8vDCMTZd23R+3GBGDBrA09+mHDKrnZZvIMzXvVNpzz2laRpDfd3JtEKT8m935vNXRil3nDbK6kWT4sJ8SMwuP6K3mzVtSi9G0/p+fZ3ommOll115jZGlf+zntDEB3S42JER/ZJV3apqmOdIY1H2olFp1+H6lVKVSqqrp6/WAo6Zp3b5NFxc2kK9uPokJQ7y547Nk7v9yB/Umc4ePeeHHvXy2LYdbZo6wWsqKsL5bZo0k2MuFB1ZbJ6WwrsHM+h15LN+YYRepI5839a47q5ulxh9eGIWjg44H1uzs9r83u6SGs1/fRNKBcl66MIYbTxlh07UHegcdL1wwAYDbPk2yesPunkorMHDzR4lEBHryvwsm9Flxn3A/d16+MIavbj6RudGBPHPuOC45YWifjKU7Th7lz8MLo/gltYjHvm6clc+rqOV/P6RxaoQ/c2xQdbc1vYOOe+ZGklFczcdbslu2pxYYemV9XbNwP/cez9jVNZh5fP0eIgM9WDLJeoVmmk0c6kNVvamlcbstbEovITrYq1cq+Qrraa5OXlBp34Hdm7/vp6rexJ2n2yZLQIi+0uPATmt8x/cOsEcp9Xw7xwQ2HYemaZObrlvSk+v6N/WHu37GcD78K5vz3viTnLK27/Z/ujWbF3/ay3kTB3O7jVJ9hHW4O+t5cEEUKfkGVmzK7NY5GswWfkkp5PZPk5j46A/c+GECD6/bzYd/ZR/9wX3IbFGsSjjIjFH+3V7DFujlwv+dPorf04r4antelx+fmF3GWa9tpLTayAdXT2Hh+OBujaOrhgx047GzGoupvPLLvl65ZmeUVNVz1YqtuDo58M5lcbj3g2p50SFevHpRLOdbsXJkb7nkhKFcOT2c5Zsyee/PTB5ZtxuTRfGfhd3rIdpVs0YPYkr4QF78sbHqY12DmcyS6l5pddAszM+NnLLaHt24emdDBjlltTw4f4xNSvLHDW0swBSfbZt0zBqjicTsMqaNkNk6exPoaf8zdoWGOpZvzGTBuGAiA+27erkQh7PGX4TpwCXATE3Tkpo+ztA07XpN065vOuZcYKemacnAS8ASZYXpE72DjnvmRfLmJRPJKKpm/ssb+C2t6JBjfkkp5L4vd3LyKH+eOHusVD2yA3OiGnt7/e+HNPI62S/HYlFs3l/CfV/uYPLjP3LF8q38tKeA+eOC+fDqKZw8yp/Hvt5NelGVjUfffRv2FZNfWdeloiltuWRqGOMGe/HIV7upqO38Oohvd+azZOlm3J31fHHDNCaHW7+6ZkcWTQjhrJgQXvppL/FZfZ+KW28yc/0H8RRW1vPWpXG9UiH2eHD/maOZFTmIh9bu4pud+dwya6TVKkQejaZp3H/maEqqjbz52372FlShFL0b2Pm6Y7aobqcdF1TW8eov+5gTFcA0G/V/G+zjir+HMwk2Wme3LbOMBrOSghV2yN/DGZ2GXfeye+2XdIxmi9zoF8cka1TF3KCU0pRS45RSE5o+1iul3lBKvdF0zCtKqSil1Hil1AlKqU09H/rf5kQFsu7mEwn0dOHyd7fwwo9pWCyK7Tnl3PhhAqODPHjt4librRES1qVpGg8viMJkaSwO0B6lGn/Gj321m2lP/cySpZv5MuEgJ470561L49j6wGyePncc00f48ey543BxdOD2T5OsWjXSmj6Pz8HbzZFZo3vWS6ext91YSqrqefa7lE495p0NGdzwYTyjgzxZdeM0hvvbvpBEWx5ZFEWIjyu3fpLU7cq31qCU4r5VO9maWcZz542XNRhW5KDTeOnCGKKCPYkM9OCak4b16vXHDfZm4fhg3t6wn9/3Nt4I7M3Abph/U2XMbq6ze/rbFExmxf1njLHmsA6haRoTQ33YZqMbLBvTi3F00Jhko0bawnYcHXT4ezjb7YzdwfJaPvorm/MmDrZpf14h+krf5xVZSZifO1/eOJ37V+/ghR/3Ep9Vxp68SnwHOLHs8kl223D0eBXq68bNM0fw3Pdp/JJaeEjj0H2FBtYm5bI2OZfMkhocHTRmjPLn3jMiOW1MAG5OR/6sAzxdePKssdzwYQIv/7SXO07vX1WwKmoa+G5XPhdOGoKzvudFHKJDvLhiejjLNmZwduzgdvszmS2KR7/azfJNmcyNCuSFJRN6pYhEezxcHHnhghjOf/NPHly9kxeWxPTJON78fT9fJJWqiGgAACAASURBVORw66yRLOildNTjibuzntU3TsdkUTjpe/+G211zIvh2Zz4v/rgXJ72OoW30CLSVll52xV2fsUs6UM6qhINcP2O4zWc548J8+HZXPoWVdVbvM/dnegkxQ3zafK0W/V+glyv5drrG7qUf9wJw86yRfTwSIWzjmJrCcnVy4L/njeeJs8by1/5STBbFiisnW7XHlOg915w8jGH+7jy0ZhfpRVW8/ms68178g9nP/87Lv+wj2NuVp84ey9b7Z/P2ZZNYNCGkwzcK88YGcU7sYF75ZZ/NS3l31brtuRhNFs6daL11U3ecNopATxfuW7WjzVnKGqOJ696PZ/mmTK4+MZzXLo7t06Cu2cShPtwycySrk3JZndj7LRC+25XP09+mMH9cELfNlj/+tqJ30PXZ79uQgW5cNm0oRrOFkYMG2GSdWnsGujvh4aLvci+7xvYGu/D3cOammSNsNLq/xQ5tvBlk7dfKipoGdh6sYKr0r7NbgZ7O5NvhjN3+oio+T8jh4hNCCZHUenGMOqYCO2hMIbloSijf3HYSa/45vc9SykTPOesdeGxRNNmlNcz67288/W0KznodD84fw1/3zuKja05gyeRQvN2cOn3OhxeOIdjblTs+S6K63mTD0XfNyvgcIgM9iA6x3kJud2c9/1nYWIhm2YaMQ/YVGupYsnQzP6c0Nrd+YP6YPqv22JZ/njqcuKE+PLC6d1sg7Mqt4LZPkhgX4sVz542XNbnHsJtOHYmPmyPjutAv0ho0TSPcr+stD9Ym55KQXc5dcyJ6JQMlOtgLJ73O6oHd5owSLAqm22h9oLC9IC9Xuwzs/vfjXpwcdNx4iu1vjAjRV465wK7ZcP8BDPWV/Gl7N22EHw/OH8NdcyL4/a5TWf3P6Vx5Yni3U4M8XBx5/vwJZJfW8OhX7a/f6017CwwkHyjn3G70rjua06MCOX1MAP/7Ma0lQNpXaODs1zaxt6CKpZfE9cvm1noHHf+7YAIacOsnib3SAqHQUMfVK7bh7ebIW5fG9YvZS2E7Xm6OfHPrydx3xuhev3aYb9daHtQYTTy5PoWxIV6cGzvYhiP7m5Nex/jBXlavjPlnegkujjpZt2rHAr1cMNSbqOpHN0ePZnduJeuSc7liehj+Hs59PRwhbOaYDezEsePKE8P556kjrLamZHL4QK6fMZxPth7gh90FVjlnT3yekINep7G4m73rjubhhVE4aBoPrtnJpvRizn5tE3UNFj697gRmjwmwyTWtobkFQkJ2OS//bNsWCHUNZq55L57ymgbeujTO6muKRP8U6OWCh0vv91EL83Mnt7z2qP1Xm73x237yK+t4cEHvzqzHDvVh58EK6ho6N87O+P/27j06rrO89/jvmRmNRtfRSLYl2dLIdi5OHBI7kpOQpFzSQBpCQwolTbrO6UmBnrRAV9u0PW0oLW1pOS2UUuih65QU6KGrbbgEaEJIITcoECcktrFzs53YJtbI8kWyZFnS6Daa9/wxM45i6z57bprvZ61Z2tqzZ+93269GeuZ93+d54kC/rljfWJC1lfBGpkh5KY3afeqR/aoLBfTrbzyv0E0Bcop3VpSlu95yoS5ZW6+7v/6s+oYnCtaOxHRS39h1RG/etEaranPzKeLahir97g2b9L39ffrvn/+x1tSH9M0PXKPL2or/E/Nbtq7Tuy5fp//z+Mva8UpuMvQ55/S/7ntWe2Kn9Onbt+p16/I7NQ/lZ8OqaiWdFjXN+MipMX3uvw7q5i1rdcX6/JYg2dbRqKlpp+eODHlyvhPD43r5xAjTMEtcc31pBXa7ugf16N4T+vU3blS4Ov8f5AD5RGCHshQM+PTp27ZqZCKhP/z6s/KgrOKy/PDlfvUNT+jWbbmdXnXH1R26akOj3nDBan39/deoPY9ZALP157dcorZIdc5KIPz9Ywf0rT29+oMbN+nnLmnx/PzA2ZaSGfOvHtorM+nut12U62adozOa+vBnxyveTMd88uBJSdI1JE4paZkRu8XWmS20T353v5pqgnrPtRsK3RQg5wjsULYuaK7T3W+7SI/vO6F/f7q7IG342s6YGmuCrynnkAsBv09fvvP1+tJ7r1S4qrQ+sawLVegzt2/VsdPj+uNvPu9pEP7gs736u0df0rs61+n9b2KKDvIjUz9rocyYz7wyoAefPao733heQbL4NdVWauOqGs8SqGw/cFL1oYAuWcuoeCnLjNgdL4GSB9sP9Gv7wZP6wHXnq4ayVygDBHYoa3dcvV5vuGCV/vLBvTrUN5LXaw+OTurRF0/oF7auy8t6k1LO8Hh5NKLfuf4CPbCnV//6425P1vzsiZ3S7311j7Z1RPRX77q0pP99UFoaqoNqqK7QoXkCu2TS6c+/9YJa6kP6jTflt4j7TJ0dEe3qHvTkA5Xth/r1+o1N8hdRBl4sXajCr0h1RdEXKXfO6W8e3q/WcEj/7apooZsD5AWBHcqaz2f6m3dvUTDg011f3TNrvbdceWBPryank3p3V36y3JW6D1x3vq5c36g/+Y/ntfkj39HPfvL7ev+/7tTfPfKS/vO5ozrYN6Lp5OL++Dw6NKb/+S87tLquUp/7lS5PisIDS7G+qWbeEbv7dvbo+SOn9aGbLipoIe+ujogGRieXlMVzNrGBuGIDY6yvWyFaSqDkwWN7T+gn3af0W9dfQJZjlA3GpVH2WsIh/e93XqoP/vsuffbxA7rrrRfm5br37ezR5tZ6bV7rXe26lczvM/2/916h7+3r0/7jw9p/7LT2HRvWd144psxgQmXAp/PX1GpTS502NdelvrbUqaU+dGZELj6Z0K99aYfik9P611+7Sk05SloDzGfDqho9dejkrM8Nj0/pE9/dr85og96xZW2eW/Za22YUKt+YRV3Y7Qf7JbG+bqVoDYeKesQumXT65MP71dFUzYenKCsEdoCkt1/Wqsf2rtNnv3dAb960WpdHIzm93r5jp/XckSH96c2bc3qdlaY6GNDbL2vV29V6Zt/Y5LQOnBg5E+ztPz6iJw706xu7jpw5pj4UOBPkvdIf196jp/WFX71CFzbXFeI2AK1vqtE3f3JEY5PTqgq+djThH753UP0jE/rCHdsKPkX4vNW1qg8FtPPwoG7d1r7s82w/eFKr6yp1/prlB4coHi3hkPbEThW6GXN68Lmj2ndsWJ++basq/ExOQ/kgsAPS/uyWS/Tjnw7orq/s1rd/6w05XWh9344eVfhNt2zNTe26clIV9OvStrAubXttQoZT8UntPzasl44Pa1/66/27ezU6kdCf/PzmnCesAeazYXUqgcrhgVFd1PLqqP3hk6P64o9+ql/sbNOWIiji7fOZOjsiWSVQcc5p+8GTuua8poIHqvBGS31IJ0cnNT41XXTTHBPTSX36kZe0qblONxd4xBvINwI7IK0+VKFP/dIW3f5PT+kvv71Xf/WuS3NynanppP5j9xFdf1GzGmuCObkGUgkqrtrYpKs2vjr1yzmnsanpgq5ZAiRpQ9OrmTFnBnYf+/ZeBfymP7hxU6Gado5tHRF9f3+fTsUn1VC99PesAydG1Dc8wTTMFaQlXfLgxOkJRZuKq3zON3Yd0aH+UX3uV7pI1IOyw/g0MMNVG5t05xs36t6nu/Xoi8dzco3v7+9T/8gk8/4LwMwI6lAU1q9K/TE8s5bd9gP9evjF4/rgdeefSSlfDDrT6+x+0r28qXdPHMisryNxykpRzLXs7t9zRBesqdUNm5sL3RQg7wjsgLP87lsv1MWt9br7G8+qf2TC8/PftzOmVbWVetOm1Z6fG0BpqAtVaFVt8ExmzMR0Uh998EW1Rar0vp8prkLKW9sb5PfZsqdjbj94Uu2NVWpvLK6RHSxfJrA7VoS17A6fjOvi1nqm/aIsEdgBZ6kM+PWZ27fq9HhCd3/9WU8LYp8cmdBje0/onZevZUE3UObWN9XopydTgd29z8S079iwPnzTxUW3Zqk6GNDm1nrtODyw5NdOJ52eOnRS1zJat6JkRpSLreRBYjqpo0PjivIhAsoUf1kCs7iwuU5/eONFenTvCd37dMyz896/u1eJpNO7u5afXQ7AyrB+VaqW3VB8Sp96eL+u2tCoG1/XUuhmzaqrI6I9saEl1/p8oXdIp8cTupr1dStKXahCtZWBoit5cHRoXNNJp/bGqkI3BSgIAjtgDu+5Zr2uPb9Jf/Hgi1kX58342s4eXdYW1qYW0uwD5W7DqhqdGJ7QX39nr06NTekjN28u2uljXR0RjU1Na+/R00t63faDqVp9BHYrT0s4pONFNhWzeyC1ZrU9wogdyhOBHTAHn8/0yVu3qMJvuusru5f8SfXZXugd0t6jp3UrSVMAKDUVU5LufTqm26+I6pK14QVeUThdMwqVL8UTB/p1YXOt1tQVTzIYeKOlvviKlMcygR1TMVGmCOyAebSGq/Sxd16q3bFTuuzPHtYtn/2R/uC+Pfr8Dw/pRy/368Tw+KLX4H1tR4+Cfh91dQBIejUzZl1lQL93w4UFbs381jZUaW04tKTAbjKR1DOvDJANc4VqCYeKbo1dbDAuv8/OJHcByg15v4EF3LxlrYIBn358aED7j5/W4/v69NUdPWeeb6wJ6sLmWl3UUq8Lm+u0qaVOFzbXqi5UceaYyURS9+8+orde0rysOlAAVp6Nq2q1qrZSv339+VpVW1no5ixoqYXKd8dOaXwqSf26Fao1HNKJ4XElppMKFEkysO6BMa1tCBVNe4B8I7ADFuHnLmnRz13yalKDkyMT2n98WPuPDeul48Pad2xYX9sR0+jk9Jlj1jVUpYO8OjnnNBifonYdgDOqgn49/UfXy1ciRZS7OiJ68Nmj6j01prUNCyeneOJAv3yWqg+KlaclHFLSSf0jk2cKlhdabCDO+jqUNU8COzO7UdJnJPklfd4599dnPV8p6V8kdUk6Kek259wrXlwbKISm2kpdU1v5milGzjn1DI6dCfReSgd+P3y5T1PTTq3hkN54AbXrALyqVII6SdrW0ShJ2nF4UO9YRGD35MGTet26sMJVFQsei9LTUv9qkfJiCex6BuN6y8UUJkf5yjqwMzO/pH+Q9FZJPZKeMbMHnHMvzjjsfZIGnXPnm9ntkj4u6bZsrw0UEzNTe2O12hurdf2MXyxT00n9tH9UdaGA/CX0RxwAzHRRa52qKvzadXhQ71hgrXB8MqGfxAb1vp/ZmKfWId8ywVyxrLMbnUiof2SSxCkoa15MQr5S0gHn3CHn3KSkL0u65axjbpH0pfT2fZKut2LN6Qx4rMLv04XNdWoNU1cHQOmq8Pu0tb1hUevsnnllUFPTjvV1K1jmd1qxZMbsGRyTJLVF+F2L8uVFYLdO0swKzj3pfbMe45xLSBqSNOu7vZndaWY7zGxHX1+fB80DAABe6OqI6MWjpzU6kZj3uO0H+1XhN12xvjFPLUO+RaorFAz4iqaWXabUQZQRO5QxLwK72Ubezs7/vphjUjudu8c5t805t231atYjAQBQLLo6IppOOu3pOTXvcdsPnNTl0Yiqgv48tQz5ZmZFVcsuNkgNO8CLwK5HUvuM79sk9c51jJkFJIUlDXhwbQAAkCed0XSh8lfmno45FJ/S871Dupb6dSteMdWy6x6Iq6rCr6YaSgqhfHkR2D0j6QIz22BmQUm3S3rgrGMekHRHevvdkh53i63qDAAAikK4ukIXrKnVzu65A7snD52Uc9I157O+bqVrDYd0rGimYo4p2lgtUjignGUd2KXXzP2mpO9K2ivpq865F8zso2b2jvRhX5DUZGYHJP2upLuzvS4AAMi/ro6Idh0eVDI5++ezTx7sV1WFX1vaGvLcMuRbS31qxK4YPqvvGYyrvZHEKShvntSxc849JOmhs/Z9ZMb2uKRbvbgWAAAonK6OiL78TEwH+kZ0YXPdOc9vP3hSV25oVDDgxaQgFLOWcEiT00kNjE6qqbayYO1wzql7IK7Xb2SUGOWNd10AALBoXR3pdXazlD04cXpcL58YocxBmWgNZ4qUF3Y65sDopOKT02TERNkjsAMAAIu2YVWNGmuC2jFLApUnD52UJF17PolTykFLupZdoUsexNI17MiIiXJHYAcAABbNzNQZjWjXLAlUnjjQr3BVhS5urS9Ay5BvLfXFMWLXPZApdcAaO5Q3AjsAALAkXR0R/bR/VCdHJl6zf/vBk7p6Y5P8PjITloPVdZXy+6zgJQ8yxcnbI4zYobwR2AEAgCXZtv7cdXaxgbh6Bscoc1BG/D7TmrrKgo/Y9QzG1VQTVE2lJzkBgZJFYAcAAJbk0nVhVfjtNfXsnjjQL0kkTikzLeFQ4dfYDYypjfV1AIEdAABYmlCFX69bF9bOGQlUth88qTV1lTpvdW0BW4Z8a6kP6ejQWEHb0D0QV3uE9XUAgR0AAFiyrmhEzx4Z0kRiWs45bT94Utec1yQz1teVk5ZwSEcLWKR8OunUe2qMUgeACOwAAMAydHVENJlI6oXe03r5xIj6RyZ0DWUOyk5rOKT45LSGJxIFuf7RoTElko5SB4AkVpkCAIAlO1Oo/JVBVfhTo3Ssrys/zemSB8eHxlUfqsj79btXcEbMQo2ConQR2AEAgCVbUx9Se2OVdh4e1LRzijZWq20F/nGN+bWmi5QfHRrXBc11eb9+z0Bqfd9Km4pZW1ur06dPy8xUWVlZ6OagRDAVEwAALMu2jkbtODygpw6d1LWUOShLreHUiF2hatnFBuPymdTaECrI9XOlrq5O69evl8/n0/DwsJLJZKGbhBJAYAcAAJalsyOi/pFJDY8ndPV5rK8rR2vqU6NJxwpU8qB7IK7WcJUq/CvvT9rKykpFo1GtWbNGo6OjGh8vbFkJFL+V91MAAADyoisaObN99UZG7MpRZcCvpppgwYqUxwbiK24a5kxmpkgkovXr18vv9zN6h3kR2AEAgGXZ1FKn2sqANjXXaXUd64DKVUs4pGMFqmUXGxxTe+PKr2HH6B0Wg+QpAABgWfw+0+/fcKFW162s9U1YmtZwSD2D+Q/sxian1Tc8sSIzYs4mM3pXU1OjY8eOaXh4WDU1NfL5GKdBCoEdAABYtl+9dkOhm4ACa64Paefhwbxft2cwXepgBU/FnE0wGFR7e7tOnTqlvr4+BQIBhUJ8uAKmYgIAACALreGQBuNTGp+azut1Y2Ua2EmvXXtXUVHB2jtIIrADAABAFlrStezyXfIglq5hVw5r7OYSDAbV1tamlpYWjY6OamysMGsdURwI7AAAALBsmVp2+c6M2T0QV6jCp9W15Z24x8wUDoe1YcMGBYNBRu/KGIEdAAAAlq25PhXYHc9zLbvYQFztkWqZWV6vW6xmjt7F43HF43ElEolCNwt5RPIUAAAALFtLgUbsUqUOym993Xwyo3fV1dUaGhrS6OioRkZG5JyTmamiokLBYJBgeIUisAMAAMCy1VYGVBcK5LWWnXNOsYG4rlwfyds1S0lFRYVWrVqlVatWKZlManJyUhMTExoZGVE8HlcymZSZye/3KxgMyu/3F7rJ8ACBHQAAALLSGg7pWB6nYp6KT2lkIsGI3SL4fD6FQiGFQiGFw2E55zQ1NaXJyUnF43GNjIycSbri8/lUUVGhQCDAqF4JyiqwM7O/kXSzpElJByW9xzl3apbjXpE0LGlaUsI5ty2b6wIAAKB4NNeH8poVs5xLHWTLzBQMBhUMBlVbW6s1a9YokUhoampK4+PjGhkZ0ejo6Dmv8/l8Z0b5fD7fmQeKR7Yjdo9I+pBzLmFmH5f0IUl/OMex1znn+rO8HgAAAIpMazik/ceG83a97oF0YBchsPNCIBBQIBBQVVWVIpGIksmkpqamlEwmlUwmNT09fSb4SyQSSiQSmpycnDc5i5m9JvAzszOjgIvdxtJkFdg55x6e8e1Tkt6dXXMAAABQalrCVeobmdDUdFIV/tyP4lDDLrd8Pp8qKxdXRiIT/M0MApPJpBKJxJl9M4+bue2ck3NO09PTcs6dc1w2MgljcinT5gwvr5c5d/qrm/fgNC/X2L1X0lfmeM5JetjMnKTPOefumeskZnanpDslKRqNetg8AAAA5EJLfUjOSX3DE1rbkPtgKzYYV6S6QnWhipxfC/PL5ZTMswOns7+fa998+3Mlx0HkoiLdBQM7M3tUUsssT33YOXd/+pgPS0pI+rc5TnOtc67XzNZIesTM9jnnfjDbgemg7x5J2rZtW37/RwAAALBkM4uU5yWwG4izvq4MnB0sMT1zfgsGds65t8z3vJndIennJV3v5giNnXO96a8nzOybkq6UNGtgBwAAgNKSqWWXrwQqsYG4Llkbzsu1gFKR1bipmd2oVLKUdzjn4nMcU2NmdZltSTdIej6b6wIAAKB4vDpil/tadtNJpyOnKE4OnC3bCbGflVSn1PTK3Wb2j5JkZmvN7KH0Mc2SfmRmeyQ9LenbzrnvZHldAAAAFIlwVYUqAz4dz0Mtu+OnxzU17UicApwl26yY58+xv1fSTentQ5K2ZHMdAAAAFC8zU2s4pKN5mIpJqQNgdlQVBAAAQNZawvkpUh5LB3ZRpmICr0FgBwAAgKy1hqt0LA9TMWODYzJTXrJvAqWEwA4AAABZW9uQGrEbn5rO6XViA3G11ocUDPBnLDATPxEAAADI2mVtDUoknV7oHcrpdahhB8yOwA4AAABZ64xGJEm7Dp/K6XVigwR2wGwI7AAAAJC11XWVijZWa+fhwZxdY3xqWsdPT5ARE5gFgR0AAAA80Rlt0M7uQTnncnL+nsFUAfRoE4lTgLMR2AEAAMATXR0R9Q1PnAnAvBYbpIYdMBcCOwAAAHji8sw6u+7cTMfsyRQnZ40dcA4COwAAAHjiopY6VQf92pWjdXbdA3EFAz6trq3MyfmBUkZgBwAAAE8E/D5tbU+ts8uF2MCY2iNV8vksJ+cHShmBHQAAADzTGY1o79FhxScTnp+bUgfA3AjsAAAA4Jmujoimk057Yt4XKu8eiJM4BZgDgR0AAAA8c3m0QZL3CVSG4lMaHk8oyogdMCsCOwAAAHimoTqo81bXeJ5A5Uypg0Zq2AGzIbADAACApzqjEe3yuFB5d7rUQRtTMYFZEdgBAADAU10dEQ3Gp/TT/lHPzhlLB3bRJgI7YDYEdgAAAPBUZ0eqUPlOD6djxgbjCldVqD5U4dk5gZWEwA4AAACeOn91repCAe3qPuXZOWMDY6yvA+ZBYAcAAABP+XyWWmfn5YjdQJyMmMA8COwAAADguc5oRC+dGNbp8amsz5VMOvUMjlHDDpgHgR0AAAA819URkXPSbg+mY54YntDkdFJtjNgBcyKwAwAAgOe2tIdl5k2h8kypg/YIa+yAuWQV2JnZn5nZETPbnX7cNMdxN5rZfjM7YGZ3Z3NNAAAAFL+6UIU2Ndd5khnzTKkDRuyAOXkxYvd3zrmt6cdDZz9pZn5J/yDpbZI2S/plM9vswXUBAABQxDo7ItrdfUrJZHaFymODcZlJ6xixA+aUj6mYV0o64Jw75JyblPRlSbfk4boAAAAooK5oRMMTCb18YiSr83QPxNVcF1JlwO9Ry4CVx4vA7jfN7Fkz+6KZRWZ5fp2k2Izve9L7ZmVmd5rZDjPb0dfX50HzAAAAUAhdHhUq7xkYYxomsIAFAzsze9TMnp/lcYuk/yvpPElbJR2V9LeznWKWfXOOxzvn7nHObXPObVu9evUibwMAAADFpqOpWo01wawTqMQG42qjODkwr8BCBzjn3rKYE5nZP0l6cJaneiS1z/i+TVLvoloHAACAkmWWfaHyicS0jp0ep4YdsIBss2K2zvj2nZKen+WwZyRdYGYbzCwo6XZJD2RzXQAAAJSGzo4GHeof1cDo5LJef2RwTM6RERNYSLZr7D5hZs+Z2bOSrpN0lySZ2Voze0iSnHMJSb8p6buS9kr6qnPuhSyvCwAAgBLQFU2ts/vJMqdjxgbHJEntBHbAvBacijkf59yvzLG/V9JNM75/SNI5pRAAAACwsl3W1qCAz7Sre1DXX9y85Ndnati1s8YOmFc+yh0AAACgTFUF/dq8tn7ZmTFjA3EF/T4114U8bhmwshDYAQAAIKc6oxHtiQ0pMZ1c8mtjg3G1Rark882WaB1ABoEdAAAAcqqzI6KxqWntOza85NfGBsbUxvo6YEEEdgAAAMipbAqVdw/E1R5hfR2wEAI7AAAA5NTacEjN9ZVLLlR+enxKQ2NTlDoAFoHADgAAADllZurqiCx5xO7VjJgEdsBCCOwAAACQc53RiHoGx3Ti9PiiX3MmsIsQ2AELIbADAABAznWm19ktZTpmbCBVnJypmMDCCOwAAACQc5esrVfQ71vSdMzYYFx1oYDC1RU5bBmwMhDYAQAAIOcqA35d2hbWru5Ti35NbCDONExgkQjsAAAAkBed0QY91zOkicT0oo7vHogzDRNYJAI7AAAA5EVXR0ST00m90Ht6wWOdc+oZHFN7IzXsgMUgsAMAAEBedEbTCVQWsc6ub3hCE4kkpQ6ARSKwAwAAQF6sqQ+pLVK1qMyY3ZQ6AJaEwA4AAAB5kylU7pyb97jYIMXJgaUgsAMAAEDedEYjOn56Qr1D8xcqz9Swa4uwxg5YDAI7AAAA5E1XulD5QvXsugfiWlNXqVCFPx/NAkoegR0AAADy5qKWOlVV+BdMoBKj1AGwJAR2AAAAyJuA36ct7eEFE6ikSh0Q2AGLRWAHAACAvOrqiOjF3tMam5y9UPlkIqmjQ2NqZ30dsGgEdgAAAMirzmhEiaTTsz2nZn2+99SYko6MmMBSENgBAAAgry5PFyrfOcd0TEodAEtHYAcAAIC8aqwJauOqGu06PPuIXabUAYEdsHiBbF5sZl+RtCn9bYOkU865rbMc94qkYUnTkhLOuW3ZXBcAAAClrbMjosf3nZBzTmb2mue6B+Kq8Jta6kMFah1QerIK7Jxzt2W2zexvJQ3Nc/h1zrn+bK4HAACAlaEzBNjBKwAACclJREFUGtF9O3v0ysm4Nqyqec1zscG41jVUye+zOV4N4GyeTMW01McsvyTpXi/OBwAAgJUtU6h8tnp2PQNxpmECS+TVGrs3SDrunHt5juedpIfNbKeZ3enRNQEAAFCiLlhTq7rKwKwJVLoH4mqLENgBS7HgVEwze1RSyyxPfdg5d396+5c1/2jdtc65XjNbI+kRM9vnnPvBHNe7U9KdkhSNRhdqHgAAAEqQz2faGm04Z8RuZCKhwfiUoozYAUuyYGDnnHvLfM+bWUDSuyR1zXOO3vTXE2b2TUlXSpo1sHPO3SPpHknatm2bW6h9AAAAKE1dHRF95rGXNTw+pbpQhSQpNpApdUBxcmApvJiK+RZJ+5xzPbM9aWY1ZlaX2ZZ0g6TnPbguAAAASlhnNCLnpD2xV/PvdWcCO6ZiAkviRWB3u86ahmlma83sofS3zZJ+ZGZ7JD0t6dvOue94cF0AAACUsK3RBplJO2dMx8yM2DEVE1iarModSJJz7ldn2dcr6ab09iFJW7K9DgAAAFaW+lCFLlxTp10zEqj0DI6ptjKghuqKArYMKD1eZcUEAAAAlqyzI6Jd3YNKJlOpFWIDcbVFqs4pWg5gfgR2AAAAKJjOaIOGxxM60DciKbXGjmmYwNIR2AEAAKBgZhYqd86pZ3CM4uTAMhDYAQAAoGA2rKpRpLpCOw8Pqn9kUmNT02qPUOoAWCoCOwAAABSMmakzmlpnlyl1EG1ixA5YKgI7AAAAFFRnR0QH+0b1/JFUPTtq2AFLR2AHAACAguqMptbZPbCnV5LURmAHLBmBHQAAAApqS3tYfp9p5+FBraqtVFXQX+gmASWHwA4AAAAFVR0M6OLWOklStJHEKcByENgBAACg4LrS0zEpdQAsD4EdAAAACq4zXc+OxCnA8hDYAQAAoOCu3NCooN+nzWvrC90UoCQFCt0AAAAAoDVcpe0f+lk11QQL3RSgJBHYAQAAoCisqq0sdBOAksVUTAAAAAAocQR2AAAAAFDiCOwAAAAAoMQR2AEAAABAiSOwAwAAAIASR2AHAAAAACWOwA4AAAAAShyBHQAAAACUOAI7AAAAAChxBHYAAAAAUOLMOVfoNszJzMYkvVDodhRYWNJQoRtRBKKSugvdiCJAf6AvZNAX6AsZ9IUU+gN9IYO+kEJ/WDl9ocM5t3qhg4o9sOtbzE2sZGZ2j3PuzkK3o9DoCyn0B/pCBn2BvpBBX0ihP9AXMugLKfSH8usLxT4V81ShG1AEvlXoBhQJ+kIK/YG+kEFfoC9k0BdS6A/0hQz6Qgr9ocz6QrEHduU+fCznHD+UKWXfFyT6Qxp9QfSFNPqC6AszlH1/oC+cUfZ9QaI/pJVVXyj2wO6eQjcARYO+gAz6AjLoC5iJ/oAM+gIyyqovFPUaOwAAAADAwop9xA4AAAAAsAACOwAAAAAocXkN7Mzsi2Z2wsyen7Fvi5k9aWbPmdm3zKw+vT9oZv+c3r/HzN484zW3mdmzZvaCmX0in/cAb5hZu5l9z8z2pv8ffzu9v9HMHjGzl9NfI+n9ZmZ/b2YH0v/3nTPO9R0zO2VmDxbqfrB8XvUFM+sws51mtjt9nt8o5H1heTx+b5hO94fdZvZAoe4Jy+Phe8N1M/rBbjMbN7NfKOS9YWk8fl/4uJk9n37cVqh7wvIsoy9cZKk4Y8LMfv+sc50Tl5Q851zeHpLeKKlT0vMz9j0j6U3p7fdK+ov09gcl/XN6e42knUoFok1KFRpcnX7uS5Kuz+d98PCkL7RK6kxv10l6SdJmSZ+QdHd6/92SPp7evknSf0oySa+X9OMZ57pe0s2SHiz0ffEoXF+QFJRUmd6ulfSKpLWFvj8ehekP6edGCn0/PIqjL8w4Z6OkAUnVhb4/HvnvC5LeLukRSQFJNZJ2SKov9P3xyGlfWCPpCkkfk/T7Z53rnLik1B95HbFzzv1AqTfUmTZJ+kF6+xFJv5je3izpsfTrTihVh2KbpI2SXnLO9aWPe3TGa1AinHNHnXO70tvDkvZKWifpFqWCdaW/Zj5VvUXSv7iUpyQ1mFlr+vWPSRrOZ/vhHa/6gnNu0jk3kT6mUkw1L0levjegtOWoL7xb0n865+I5vwF4xsO+sFnSfznnEs65UUl7JN2Yx1tBlpbaF5xzJ5xzz0iamuVcs8UlJa0Y/vB5XtI70tu3SmpPb++RdIuZBcxsg6Su9HMHJF1kZuvNLKDUf1y7ULLMbL2kyyX9WFKzc+6olPrhVeqTFin1Qxub8bKe9D6sINn2hfQUjWfTz3/cOdebn5YjFzx4bwiZ2Q4ze4qpd6XNw98Tt0u6N5dtRW5l2Rf2SHqbmVWb2SpJ14m/IUvWIvtCWQkUugFKTb/8ezP7iKQHJE2m939R0sVKDZMflrRdUsI5N2hm75f0FUnJ9P6NeW81PGFmtZK+Lul3nHOnzWzOQ2fZR62OFcSLvuCci0m6zMzWSvoPM7vPOXc8Jw1GTnn03hB1zvWa2UZJj5vZc865gzloLnLIq98T6RGbSyV91/NGIi+y7QvOuYfN7Aql/nbsk/SkpEROGoucWkJfKCsFH7Fzzu1zzt3gnOtS6lO0g+n9CefcXc65rc65WyQ1SHo5/dy3nHNXOeeulrQ/sx+lxcwqlPqh/Dfn3DfSu49nps6kv55I7+/Raz9Va5PEaMwK4XVfSI/UvSDpDblsN3LDq/6QGbF1zh2S9H2lPtlFCfH4veGXJH3TOXfOlCwUPw/fFz6W/tvyrUoFgPwNWWKW2BfKSsEDOzNbk/7qk/THkv4x/X21mdWkt9+q1Gjdi2e9JiLpA5I+X4CmIwuW+mjlC5L2Ouc+NeOpByTdkd6+Q9L9M/b/j3Smq9dLGsoMuaO0edUXzKzNzKrS54xIulapD35QQjzsDxEzq0yfc5VS/eHFvNwEPJGD3xO/LKZhliQP3xf8ZtaUPudlki6T9HBebgKeWEZfKCvmXP5ms5nZvZLeLGmVpOOS/lSp7HUfTB/yDUkfcs659LzZ7yo13fKIpPc55w7POM+W9Gs+6pz7cp5uAR4xs5+R9ENJzyn1fyxJf6TUPOmvSooqlf30VufcQPoH+bNKLXKOS3qPc25H+lw/lHSRUn3ppFJ9hak2JcKrvpD+AOhvlZp6ZZI+65y7J683g6x52B+ukfS59Dl8kj7tnPtCXm8GWfH498R6SU9IanfOJYWS4uH7QkjSrvTrT0v6Defc7vzdCbK1jL7QonT20/TxI5I2p6dvnhOXlPrvibwGdgAAAAAA7xV8KiYAAAAAIDsEdgAAAABQ4gjsAAAAAKDEEdgBAAAAQIkjsAMAAACAEkdgBwAAAAAljsAOAAAAAErc/weK9Chg4a9PwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 5))\n",
    "\n",
    "# Plot the data (here we are subsetting it to get a better look at the forecasts)\n",
    "endog.loc['1999':].plot(ax=ax)\n",
    "\n",
    "# Construct the forecasts\n",
    "fcast = res.get_forecast('2011Q4').summary_frame()\n",
    "fcast['mean'].plot(ax=ax, style='k--')\n",
    "ax.fill_between(fcast.index, fcast['mean_ci_lower'], fcast['mean_ci_upper'], color='k', alpha=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note on what to expect from forecasts\n",
    "\n",
    "The forecast above may not look very impressive, as it is almost a straight line. This is because this is a very simple, univariate forecasting model. Nonetheless, keep in mind that these simple forecasting models can be extremely competitive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prediction vs Forecasting\n",
    "\n",
    "The results objects also contain two methods that all for both in-sample fitted values and out-of-sample forecasting. They are `predict` and `get_prediction`. The `predict` method only returns point predictions (similar to `forecast`), while the `get_prediction` method also returns additional results (similar to `get_forecast`).\n",
    "\n",
    "In general, if your interest is out-of-sample forecasting, it is easier to stick to the `forecast` and `get_forecast` methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross validation\n",
    "\n",
    "**Note**: some of the functions used in this section were first introduced in Statsmodels v0.11.0.\n",
    "\n",
    "A common use case is to cross-validate forecasting methods by performing h-step-ahead forecasts recursively using the following process:\n",
    "\n",
    "1. Fit model parameters on a training sample\n",
    "2. Produce h-step-ahead forecasts from the end of that sample\n",
    "3. Compare forecasts against test dataset to compute error rate\n",
    "4. Expand the sample to include the next observation, and repeat\n",
    "\n",
    "Economists sometimes call this a pseudo-out-of-sample forecast evaluation exercise, or time-series cross-validation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will conduct a very simple exercise of this sort using the inflation dataset above. The full dataset contains 203 observations, and for expositional purposes we'll use the first 80% as our training sample and only consider one-step-ahead forecasts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A single iteration of the above procedure looks like the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intercept    1.162076\n",
      "ar.L1        0.724242\n",
      "sigma2       5.051600\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "# Step 1: fit model parameters w/ training sample\n",
    "training_obs = int(len(endog) * 0.8)\n",
    "\n",
    "training_endog = endog[:training_obs]\n",
    "training_mod = sm.tsa.SARIMAX(\n",
    "    training_endog, order=(1, 0, 0), trend='c')\n",
    "training_res = training_mod.fit()\n",
    "\n",
    "# Print the estimated parameters\n",
    "print(training_res.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        true  forecast    error\n",
      "1999Q3  3.35   2.55262  0.79738\n"
     ]
    }
   ],
   "source": [
    "# Step 2: produce one-step-ahead forecasts\n",
    "fcast = training_res.forecast()\n",
    "\n",
    "# Step 3: compute root mean square forecasting error\n",
    "true = endog.reindex(fcast.index)\n",
    "error = true - fcast\n",
    "\n",
    "# Print out the results\n",
    "print(pd.concat([true.rename('true'),\n",
    "                 fcast.rename('forecast'),\n",
    "                 error.rename('error')], axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To add on another observation, we can use the `append` or `extend` results methods. Either method will produce the same forecasts, but they differ in the other results that are available:\n",
    "\n",
    "- `append` is the more complete method. It always stores results for all training observations, and it optionally allows refitting the model parameters given the new observations (note that the default is *not* to refit the parameters).\n",
    "- `extend` is a faster method that may be useful if the training sample is very large. It *only* stores results for the new observations, and it does not allow refitting the model parameters (i.e. you have to use the parameters estimated on the previous sample).\n",
    "\n",
    "If your training sample is relatively small (less than a few thousand observations, for example) or if you want to compute the best possible forecasts, then you should use the `append` method. However, if that method is infeasible (for example, becuase you have a very large training sample) or if you are okay with slightly suboptimal forecasts (because the parameter estimates will be slightly stale), then you can consider the `extend` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A second iteration, using the `append` method, would go as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intercept    1.171544\n",
      "ar.L1        0.723152\n",
      "sigma2       5.024580\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "# Step 1: append a new observation to the sample and refit the parameters\n",
    "append_res = training_res.append(endog[training_obs:training_obs + 1], refit=True)\n",
    "\n",
    "# Print the re-estimated parameters\n",
    "print(append_res.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that these estimated parameters are slightly different than those we originally estimated. With the new results object, `append_res`, we can compute forecasts starting from one observation further than the previous call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        true  forecast     error\n",
      "1999Q4  2.85  3.594102 -0.744102\n"
     ]
    }
   ],
   "source": [
    "# Step 2: produce one-step-ahead forecasts\n",
    "fcast = append_res.forecast()\n",
    "\n",
    "# Step 3: compute root mean square forecasting error\n",
    "true = endog.reindex(fcast.index)\n",
    "error = true - fcast\n",
    "\n",
    "# Print out the results\n",
    "print(pd.concat([true.rename('true'),\n",
    "                 fcast.rename('forecast'),\n",
    "                 error.rename('error')], axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting it altogether, we can perform the recursive forecast evaluation exercise as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          1999Q2    1999Q3    1999Q4    2000Q1    2000Q2\n",
      "1999Q3  2.552620       NaN       NaN       NaN       NaN\n",
      "1999Q4  3.010790  3.588286       NaN       NaN       NaN\n",
      "2000Q1  3.342616  3.760863  3.226165       NaN       NaN\n",
      "2000Q2       NaN  3.885850  3.498599  3.885225       NaN\n",
      "2000Q3       NaN       NaN  3.695908  3.975918  4.196649\n"
     ]
    }
   ],
   "source": [
    "# Setup forecasts\n",
    "nforecasts = 3\n",
    "forecasts = {}\n",
    "\n",
    "# Get the number of initial training observations\n",
    "nobs = len(endog)\n",
    "n_init_training = int(nobs * 0.8)\n",
    "\n",
    "# Create model for initial training sample, fit parameters\n",
    "init_training_endog = endog.iloc[:n_init_training]\n",
    "mod = sm.tsa.SARIMAX(training_endog, order=(1, 0, 0), trend='c')\n",
    "res = mod.fit()\n",
    "\n",
    "# Save initial forecast\n",
    "forecasts[training_endog.index[-1]] = res.forecast(steps=nforecasts)\n",
    "\n",
    "# Step through the rest of the sample\n",
    "for t in range(n_init_training, nobs):\n",
    "    # Update the results by appending the next observation\n",
    "    updated_endog = endog.iloc[t:t+1]\n",
    "    res = res.append(updated_endog, refit=False)\n",
    "    \n",
    "    # Save the new set of forecasts\n",
    "    forecasts[updated_endog.index[0]] = res.forecast(steps=nforecasts)\n",
    "\n",
    "# Combine all forecasts into a dataframe\n",
    "forecasts = pd.concat(forecasts, axis=1)\n",
    "\n",
    "print(forecasts.iloc[:5, :5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a set of three forecasts made at each point in time from 1999Q2 through 2009Q3. We can construct the forecast errors by subtracting each forecast from the actual value of `endog` at that point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          1999Q2    1999Q3    1999Q4    2000Q1    2000Q2\n",
      "1999Q3  0.797380       NaN       NaN       NaN       NaN\n",
      "1999Q4 -0.160790 -0.738286       NaN       NaN       NaN\n",
      "2000Q1  0.417384 -0.000863  0.533835       NaN       NaN\n",
      "2000Q2       NaN  0.304150  0.691401  0.304775       NaN\n",
      "2000Q3       NaN       NaN -0.925908 -1.205918 -1.426649\n"
     ]
    }
   ],
   "source": [
    "# Construct the forecast errors\n",
    "forecast_errors = forecasts.apply(lambda column: endog - column).reindex(forecasts.index)\n",
    "\n",
    "print(forecast_errors.iloc[:5, :5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To evaluate our forecasts, we often want to look at a summary value like the root mean square error. Here we can compute that for each horizon by first flattening the forecast errors so that they are indexed by horizon and then computing the root mean square error fore each horizon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           1999Q2    1999Q3    1999Q4    2000Q1    2000Q2\n",
      "horizon                                                  \n",
      "1        0.797380 -0.738286  0.533835  0.304775 -1.426649\n",
      "2       -0.160790 -0.000863  0.691401 -1.205918 -0.311464\n",
      "3        0.417384  0.304150 -0.925908 -0.151602 -2.384952\n"
     ]
    }
   ],
   "source": [
    "# Reindex the forecasts by horizon rather than by date\n",
    "def flatten(column):\n",
    "    return column.dropna().reset_index(drop=True)\n",
    "\n",
    "flattened = forecast_errors.apply(flatten)\n",
    "flattened.index = (flattened.index + 1).rename('horizon')\n",
    "\n",
    "print(flattened.iloc[:3, :5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "horizon\n",
      "1    3.292700\n",
      "2    3.421808\n",
      "3    3.280012\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "# Compute the root mean square error\n",
    "rmse = (flattened**2).mean(axis=1)**0.5\n",
    "\n",
    "print(rmse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using `extend`\n",
    "\n",
    "We can check that we get similar forecasts if we instead use the `extend` method, but that they are not exactly the same as when we use `append` with the `refit=True` argument. This is because `extend` does not re-estimate the parameters given the new observation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          1999Q2    1999Q3    1999Q4    2000Q1    2000Q2\n",
      "1999Q3  2.552620       NaN       NaN       NaN       NaN\n",
      "1999Q4  3.010790  3.588286       NaN       NaN       NaN\n",
      "2000Q1  3.342616  3.760863  3.226165       NaN       NaN\n",
      "2000Q2       NaN  3.885850  3.498599  3.885225       NaN\n",
      "2000Q3       NaN       NaN  3.695908  3.975918  4.196649\n"
     ]
    }
   ],
   "source": [
    "# Setup forecasts\n",
    "nforecasts = 3\n",
    "forecasts = {}\n",
    "\n",
    "# Get the number of initial training observations\n",
    "nobs = len(endog)\n",
    "n_init_training = int(nobs * 0.8)\n",
    "\n",
    "# Create model for initial training sample, fit parameters\n",
    "init_training_endog = endog.iloc[:n_init_training]\n",
    "mod = sm.tsa.SARIMAX(training_endog, order=(1, 0, 0), trend='c')\n",
    "res = mod.fit()\n",
    "\n",
    "# Save initial forecast\n",
    "forecasts[training_endog.index[-1]] = res.forecast(steps=nforecasts)\n",
    "\n",
    "# Step through the rest of the sample\n",
    "for t in range(n_init_training, nobs):\n",
    "    # Update the results by appending the next observation\n",
    "    updated_endog = endog.iloc[t:t+1]\n",
    "    res = res.extend(updated_endog)\n",
    "    \n",
    "    # Save the new set of forecasts\n",
    "    forecasts[updated_endog.index[0]] = res.forecast(steps=nforecasts)\n",
    "\n",
    "# Combine all forecasts into a dataframe\n",
    "forecasts = pd.concat(forecasts, axis=1)\n",
    "\n",
    "print(forecasts.iloc[:5, :5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          1999Q2    1999Q3    1999Q4    2000Q1    2000Q2\n",
      "1999Q3  0.797380       NaN       NaN       NaN       NaN\n",
      "1999Q4 -0.160790 -0.738286       NaN       NaN       NaN\n",
      "2000Q1  0.417384 -0.000863  0.533835       NaN       NaN\n",
      "2000Q2       NaN  0.304150  0.691401  0.304775       NaN\n",
      "2000Q3       NaN       NaN -0.925908 -1.205918 -1.426649\n"
     ]
    }
   ],
   "source": [
    "# Construct the forecast errors\n",
    "forecast_errors = forecasts.apply(lambda column: endog - column).reindex(forecasts.index)\n",
    "\n",
    "print(forecast_errors.iloc[:5, :5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           1999Q2    1999Q3    1999Q4    2000Q1    2000Q2\n",
      "horizon                                                  \n",
      "1        0.797380 -0.738286  0.533835  0.304775 -1.426649\n",
      "2       -0.160790 -0.000863  0.691401 -1.205918 -0.311464\n",
      "3        0.417384  0.304150 -0.925908 -0.151602 -2.384952\n"
     ]
    }
   ],
   "source": [
    "# Reindex the forecasts by horizon rather than by date\n",
    "def flatten(column):\n",
    "    return column.dropna().reset_index(drop=True)\n",
    "\n",
    "flattened = forecast_errors.apply(flatten)\n",
    "flattened.index = (flattened.index + 1).rename('horizon')\n",
    "\n",
    "print(flattened.iloc[:3, :5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "horizon\n",
      "1    3.292700\n",
      "2    3.421808\n",
      "3    3.280012\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "# Compute the root mean square error\n",
    "rmse = (flattened**2).mean(axis=1)**0.5\n",
    "\n",
    "print(rmse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By not re-estimating the parameters, our forecasts are slightly worse (the root mean square error is higher at each horizon). However, the process is quite a bit faster, even with only 200 datapoints. Using the `%%timeit` cell magic on the cells above, we found a runtime of 570ms using `extend` versus 1.7s using `append` with `refit=True`. (Note that using `extend` is also faster than using `append` with `refit=False`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexes\n",
    "\n",
    "Throughout this notebook, we have been making use of Pandas date indexes with an associated frequency. As you can see, this index marks our data as at a quarterly frequency, between 1959Q1 and 2009Q3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PeriodIndex(['1959Q1', '1959Q2', '1959Q3', '1959Q4', '1960Q1', '1960Q2',\n",
      "             '1960Q3', '1960Q4', '1961Q1', '1961Q2',\n",
      "             ...\n",
      "             '2007Q2', '2007Q3', '2007Q4', '2008Q1', '2008Q2', '2008Q3',\n",
      "             '2008Q4', '2009Q1', '2009Q2', '2009Q3'],\n",
      "            dtype='period[Q-DEC]', length=203, freq='Q-DEC')\n"
     ]
    }
   ],
   "source": [
    "print(endog.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In most cases, if your data has an associated data/time index with a defined frequency (like quarterly, monthly, etc.), then it is best to make sure your data is a Pandas series with the appropriate index. Here are three examples of this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PeriodIndex(['2000', '2001', '2002', '2003'], dtype='period[A-DEC]', freq='A-DEC')\n"
     ]
    }
   ],
   "source": [
    "# Annual frequency, using a PeriodIndex\n",
    "index = pd.period_range(start='2000', periods=4, freq='A')\n",
    "endog1 = pd.Series([1, 2, 3, 4], index=index)\n",
    "print(endog1.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DatetimeIndex(['2000-01-01', '2000-04-01', '2000-07-01', '2000-10-01'], dtype='datetime64[ns]', freq='QS-JAN')\n"
     ]
    }
   ],
   "source": [
    "# Quarterly frequency, using a DatetimeIndex\n",
    "index = pd.date_range(start='2000', periods=4, freq='QS')\n",
    "endog2 = pd.Series([1, 2, 3, 4], index=index)\n",
    "print(endog2.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DatetimeIndex(['2000-01-31', '2000-02-29', '2000-03-31', '2000-04-30'], dtype='datetime64[ns]', freq='M')\n"
     ]
    }
   ],
   "source": [
    "# Monthly frequency, using a DatetimeIndex\n",
    "index = pd.date_range(start='2000', periods=4, freq='M')\n",
    "endog3 = pd.Series([1, 2, 3, 4], index=index)\n",
    "print(endog3.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, if your data has an associated date/time index, it is best to use that even if does not have a defined frequency. An example of that kind of index is as follows - notice that it has `freq=None`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DatetimeIndex(['2000-01-01 10:08:00', '2000-01-01 11:32:00',\n",
      "               '2000-01-01 17:32:00', '2000-01-02 06:15:00'],\n",
      "              dtype='datetime64[ns]', freq=None)\n"
     ]
    }
   ],
   "source": [
    "index = pd.DatetimeIndex([\n",
    "    '2000-01-01 10:08am', '2000-01-01 11:32am',\n",
    "    '2000-01-01 5:32pm', '2000-01-02 6:15am'])\n",
    "endog4 = pd.Series([0.2, 0.5, -0.1, 0.1], index=index)\n",
    "print(endog4.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can still pass this data to Statsmodels' model classes, but you will get the following warning, that no frequency data was found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/fulton/projects/statsmodels/statsmodels/tsa/base/tsa_model.py:219: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  ' ignored when e.g. forecasting.', ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "mod = sm.tsa.SARIMAX(endog4)\n",
    "res = mod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this means is that you can't specify forecasting steps by dates, and the output of the `forecast` and `get_forecast` methods won't have associated dates. The reason is that without a given frequency, there is no way to determine what date each forecast should be assigned to. In the example above, there is no pattern to the date/time stamps of the index, so there is no way to determine what the next date/time should be (should it be in the morning of 2000-01-02? the afternoon? or maybe not until 2000-01-03?).\n",
    "\n",
    "For example, if we forecast one-step-ahead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/fulton/projects/statsmodels/statsmodels/tsa/base/tsa_model.py:576: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.\n",
      "  ValueWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4    0.011866\n",
       "dtype: float64"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.forecast(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The index associated with the new forecast is `4`, because if the given data had an integer index, that would be the next value. A warning is given letting the user know that the index is not a date/time index.\n",
    "\n",
    "If we try to specify the steps of the forecast using a date, we will get the following exception:\n",
    "\n",
    "    KeyError: 'The `end` argument could not be matched to a location related to the index of the data.'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'The `end` argument could not be matched to a location related to the index of the data.'\n"
     ]
    }
   ],
   "source": [
    "# Here we'll catch the exception to prevent printing too much of\n",
    "# the exception trace output in this notebook\n",
    "try:\n",
    "    res.forecast('2000-01-03')\n",
    "except KeyError as e:\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ultimately there is nothing wrong with using data that does not have an associated date/time frequency, or even using data that has no index at all, like a Numpy array. However, if you can use a Pandas series with an associated frequency, you'll have more options for specifying your forecasts and get back results with a more useful index."
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}