{
 "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 0x7ff99cec43a0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAEvCAYAAADW/SmEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3wbh303/s9hb5AASICUOCRRlCjJlmzJ8oxHvN24rvvLk9RJ2sTNaJq4TdrkadPkSdOmza8zSduk6ZOktZ3lOE0dx7HjEa94yLZk7UGKpCRuAiABEBs4rHv+ONwR4wCCJEAA1Pf9evllGwTBE0Ue7nvfxXAcB0IIIYQQQgghjUlW7wMghBBCCCGEEFIaBW2EEEIIIYQQ0sAoaCOEEEIIIYSQBkZBGyGEEEIIIYQ0MAraCCGEEEIIIaSBUdBGCCGEEEIIIQ1MUY8varPZuN7e3np8aUIIIYQQQgipuyNHjng4jmur5Ll1Cdp6e3tx+PDhenxpQgghhBBCCKk7hmEmKn0ulUcSQgghhBBCSAOjoI0QQgghhBBCGhgFbYQQQgghhBDSwCoO2hiGeZBhmDmGYU7nPPZXDMPMMAxzPPvPXbU5TEIIIYQQQgi5OC0n0/YwgDskHv86x3F7sv88XZ3DIoQQQgghhBACLCNo4zjuVQC+Gh4LIYQQQgghhJAC1ehpe4BhmJPZ8snWUk9iGOZjDMMcZhjm8Pz8fBW+LCGEEEIIIYSsf6sN2v4DwBYAewA4AXy11BM5jvsOx3H7OI7b19ZW0Q45QgghhBBCCLnorSpo4zjOzXFcmuO4DIDvAthfncMihBBCCCGEEAIAitV8MsMwHRzHObP/ey+A0+WeTwghpHaOT/lxfi4MvVoBg1oBvVqOfrsRevWqTvWEEEIIqbOK38kZhvkxgBsB2BiGmQbwJQA3MgyzBwAHYBzAH9TgGAkhhFTgo98/jPkQm/fYO7e348EPXVGnIyKEEEJINVQctHEcd5/Ew/9VxWMhhBCyQsF4EvMhFh+/YQvu3t2BCJvGPzx7FrP+WL0PjRBCCCGrRDUzhBCyDkx4ogCAPV0t2NlpBgD0WvV464K3nodFCCGEkCqoxsh/QgghdTbmjQAANtn04mMmrQLBWLJeh0QIIYSQKqGgjRBC1oFxDx+09Vh14mMmjRIhNoV0hqvXYRFCCCGkCihoI4SQdWDcE0GHWQONUi4+ZtIqAQDheKpeh0UIIYSQKqCgjRBC1oFxbwS9Vn3eYyYN37YcjFOJJCGEENLMKGgjhJB1YNwbRa+tIGjLZtoC1NdGCCGENDUK2gghpMkFYkn4Ign05vSzAXxPG0CZNkIIIaTZUdBGCCFNbiI7ObI405Ytj4xRTxshhBDSzChoI4SsqUA0iXgyXe/DWFfGPMXj/gHKtBFCCCHrBQVthJA19dv/cQBfe36k3oexroxnF2t3WwrKI7M9bbSrjRBCCGluFLQRQtZMMp3BBU8Ew65QvQ+lKbmDcXBc8c61cW8EnQXj/gHAqFaAYYAgjfwnhBBCmhoFbYSQNTMfYsFxgDMQq/ehNJ2zriCu/rsX8fygu+hj495IUT8bAMhkDAxqBWXaCCGEkCZHQRshZM24gnEAgNMfr/ORNJ9nTrmQ4YCXzs4VfWzcE0GPtThoA/i+NuppI4QQQpobBW2EkDXjCvDBWohNIUSBxLK8MMRn2A6c9+Q9HogmsRBNYpNNJ/VpMGmVND2SEEIIaXIUtBFC1owQtBX+NynPGYjhzGwQPVYdpnwxTHqj4sfGhXH/JTNtCsq0EUIIIU2OgjZCyJpxBxcDtVkK2ir2whBfEvn5uwYA5GfbxkvsaBPwmTYK2gghhJBmRkEbIWTNuIJxqBT8acfpp2EklXph0I1NNj1u22GH3aTG6+cWg7YxTwQMUzzuX2DSKBGi6ZGEEEJIU6OgjRCyZlyBOHZ0mMAwgJMybRUJsym8ed6Lm7e3g2EYXLvFhjfPe5HJ8KP/J7xRdJq1ReP+BSYtTY8khBBCmh0FbYSQNeMOxtFl0cFmUNPY/wq9NjKPRDqDW3bYAQDX9tngiyRwNrvrbswTQY9VOssGZDNtbArpTPF+N0IIIYQ0BwraCCFrguM4OANxOExqdJo1lGmr0PNDbpi1SuzraQXAB20AcCBbIllqR5vApFUCAMJUIkkIIYQ0LQraCCEVO+sK4p+fGwbHLT9rE4glwaYysJs06DBrMUs9bUtKZzi8fHYO79zeDoWcP107zBpsadPjwHkP/NEE/NEkNpWYHAnw0yMB0ARJQgghpIlR0EYIqdgzp1z45svnML2w/IBLWKztMGvQ0cJn2lYS/F1Mjk4uYCGaxC0D9rzHr+2z4eAFH0bnwgBQvjwym2kLUF8bIYQQ0rQoaCOEVCzC8iV2g87gsj9X2MvmMGnQYdYgmkgjSCV7Zb0w6IZSzuD6flve49dssSGWTOPxYzMAgE3lyiM1fNBGmTZCCCGkeVHQRgipWCSRBgCcmV1+0ObOzbSZtQBAw0iW8PyQG1dttsKYDbwEV2+2QsYATxybAcMAXSXG/QP89EgACMYoQCaEEEKaFQVthJCKiZm2FQRtrgALAGg3atDZogEAOP00jKSUC/NhXJiP4NYd9qKPmXVKXLLBjEgiXXbcP0CZNkIIIWQ9oKCNEFIxIWgbWkl5ZDAGm0EFlUImZtpmKdNW0olpPwA+qyZFmCLZayudZQMWe9poVxshhBDSvChoI4RULJwN2mb8MSxEEsv6XFcgDruJz7C1G9WQMZRpK8cd5DOTHS1ayY+LQVuZyZEAYFQrwDCg/kFCCCGkiVHQRgipWDSRhkbJnzaWm21zBVk4skGbQi5Du5F2tZXjDsZhUCtgUCskP763pxWbbHpcWSITJ5DJGBjUCsq0EUIIIU2MgjZCSMUibAqXdfFLnpc7QdIdjMNu1oj/z4/9p/LIUuaCLNpN6pIf1yjlePmzN+I3d3cu+VomjZJ62gghhJAmRkEbIaRiYTaFHqsOdpN6WRMk2VQavkhCzLQBQKdZS5m2MlzBOOxGzdJPrIBZq6xoeuRZVxC3f/1VeMJsVb4uIYQQQqqDgjZCSMWiiTR0KgV2dJiWNUFyLtuf5cjNtJk1mPXHaMF2Ce5gPO/7tRomraKiTNvTJ50YdodwaiZQla9LCCGEkOqgoI0QUhGO4xBJpGBQy7Gj04Rz82HEk+mKPteZs1hb0NGiBZvKYCFKZXuFOI5bsjxyOUwaZUU9bQfHfACAaV+0Kl+XEEIIIdVBQRshpCLRRBocB+jVCuzsNCOd4TDiDlX0ua6cxdqCjux/U19bsYVoEol0pmrlkSatEqElpkfGk2kcm+LXDExS0EYIIYQ0FAraCCEViST4i36dmi+PBCpfsu3OZtrsJomgjcb+F3FLBLmrUUmm7cSUH4lUBgAw5aNAmhBCCGkkFLQRQioSYflSSINajm6LDga1ouIJkq5gHFqlHCbN4vj6zuz+Mcq0FROCNnu1yiO1CoTYFNKZ0v2DB8d8YBjg8u4WTC1Qpo0QQghpJBS0EUIqEsku1tarFJDJGAx0GCvOtLmyQzUYhhEfsxnUUMgYzNIEySJC0NZerfJIjRIAEC5TInlwzIttdiN2bTBjisojCSGEkIZCQRshpCJC0CYse97RYcKQM4hMTvZmLhTHJ350pOii3x2I5w0hAQC5jIHdpIGLgrYi7uy0zaoNItHyQVupCZKJVAZHJhZw1WYrulp1CMZTCNCAGEIIIaRhUNBGCKlIbk8bAOzsNCOSSGMiG6BxHIfPPXYKT59y4QdvTeR9rqvE+Hph7D/J5w7GYdGroFbIq/J6QllqoERf26mZAOLJDK7cZEGXRQcAVCJJCCGENBAK2gghFQnn9LQBwI7O/GEkj749hZfOzqFFp8Qvjs+K/VOZDAd3MJ43hETQ0UILtqW4g3G0G6uTZQOWzrQdHPMCAK7YZEGXhe81pBJJQgghpHFQ0EYIqYjY05bNtG21G6CQMTgzG8CkN4q/eWoQ12yx4sv37IIrGMfBC3wg4IsmkExzcEiU+nWa+fLITJkBGRcjd5CVDHJXSuhpC8ake9oOXvChr90Am0G9ZKZt3BOpeD8fIYQQQqqDgjZCSEWEoE2n4oM2tUKOvnYDTs0E8JmfHoecYfBP/2s3btthh0GtwM+PzwCA2LNWqjwykc7AG0ms0Z+iObiDxT2Aq2HS8n9nUpm2VDqDw+M+XLnJwj9Xo4RZq5Tc1RZmU7j9X17Fw2+MV+3YCCGEELI0CtoIIRURRv7rVYt9Vjs6TXht1IO3xxfw1/fsxIYWLTRKOW7f6cAzp1yIJ9M54+uLgxCHmS/Fo2Eki1LpDDxhtmrj/oGc8kiJnrZBZxCRRBpXbraKj3VbdJK72gZng2BTGZyeCVTt2AghhBCyNAraCCEViSRS0ChlUMgXTxvCku07djpw72UbxMfvvWwDQmwKL52dg6vMoujOFv6xWdrVJvKEE8hwQHsVM20GlQIMAwQlRv4fvOADAFyVzbQBQJdFK1keKQRr5+bCKz6Ws64g/vrJM+A4KoklhBBCKlVx0MYwzIMMw8wxDHM65zELwzDPMwwzmv13a20OkxBSbxE2Bb1KkffY7TsduGdPJ75y7668HWxXb7Gi3ajGz4/NwB2IQ8YAbYbizFFHNtPmpAmSonKZyZWSyRgY1QrJTNvBMS822fR5QWJXqw7TC7GiXsPTs3zQdsETQSqdWdGxvDDoxkMHxrFAKwUIIYSQii0n0/YwgDsKHvscgBc5jtsK4MXs/xNC1qEImxKHkAi6LDr86+9cBmtBQCaXMbh7dydeHp7DWVcIbUZ1XoZOYNWroJLLaIJkDiFoq2ZPG8CXSBb2tKUzHA6N+bC/15L3+EaLDolUBnMhNu/xMzNByBh+r9vUwsoCbWEKqdAjSQghhJClVRy0cRz3KgBfwcP3APhe9r+/B+C3qnRchJAGE2bTRUFbOfdetgHJNIcXhtwlAxCZjIHDrMEsBW2ixUxb9XraAH7ASOH0yLOuIILxFK7cnB+0dUtMkIwl0hidC+HaPhuAlZdIRrP7/oS9f4QQQghZ2mp72uwcxzkBIPvv9tUfEiGkEUUTqbwhJEvZ2WnCljY9Mlz5Uj+HSQM3BW0id5CFXMYUZS9Xy6RVFGXaDo3x9+Fyh5AAQFdr8a62IVcQGQ64Zw/fuzg6F1rRcYSzGTbKtBFCCCGVW7NBJAzDfIxhmMMMwxyen59fqy9LCKkSqfLIchiGwW9lL/ClhpAIrAYVvBG25MfXi8HZIN4471nyee5gHG0GNeQyZsnnLgefacsP2g5PLKDTrMGGFm3e4xtatWAY5I39P5MdQnL1FiscJg3OuVeWaYuIQRvteiOEEEIqtdqgzc0wTAcAZP89V+qJHMd9h+O4fRzH7Wtra1vllyWErLUwm4JhGUEbwGdlGGax3E6KRa+C7yLY0/YXPzuJ9333ID72/cOYKTN4xRWMV700EuB72kIF0yOPTixgb0E/G8Dv4LMbNXlj/0/PBNGqU6LTrMFWuwHn5lcatFFPGyGEELJcqw3afgHgg9n//iCAJ1b5eoSQBhVNpKFbRnkkAHRbdXjygevwviu7Sz7HalDDH0uueBphM8hkOAy7Q9juMOK1UQ9u+eor+I9fn0ciVfxnnguyVR33LzBplAjkZNpm/DE4A3Hs7W6RfH63RZfX03Z6NoBdG8xgGAZb2gw4Nxcumi5ZCaGXLUxBGyGEEFKx5Yz8/zGANwFsYxhmmmGYDwP4ewC3MgwzCuDW7P8TQtah8DLLIwW7NpihU5X+PKteBY7Duh4BP7UQRTyZwf3X9uL5P70e79hqwz88exYfeuhQ0XPdoVpl2hQIsykxOD4ysQAA2CeRaQOAjRYtprPlkWwqjRF3CDs7zQCArXYDoon0ivbrCRm2aILKIwkhhJBKLWd65H0cx3VwHKfkOG4jx3H/xXGcl+O4mzmO25r9d+F0SULIOsBxHCIrKI+shEWvAoB1XSI5ku3/2mo3YmOrDt/5vX341M1b8cZ5b96wj3gyDX80WfVx/wCfaQMWM1xHJxagU8mx3WGUfH5Xqw7OYBxsKo1RdxjJNIddG/hl6lvb+c9ZyQRJoTySMm2EEEJI5dZsEAkhpHnFkxlkOECnXl55ZCWs2aBtPQ8jGXHzkxa3thvEx+7e3QkAeHV0cTDTXJD/HtSkPFLLB23C2P/DEz7s6WqR3J8H8Dv4OA6Y9cdxOjuEZFc209aX/XOsJGgLi5k2CtoIIYSQSlHQRghZktCHVJNMm+FiyLSFsKFFC2M22wUAW9r02NCixSvDi0GbOyTsaKtFpo3/uwvGk4iwKQw5Q9jb01ry+eKuNl8Up2cDMKoV4mMWvQpWvQqjy5wgKWRsAZoeSQghhCwHBW2EENF/vnYBDx8YK3pcuNDWl+lNW6mLoTxy2BVCv92Q9xjDMLhxWxsOnPOIA0mExdo1KY8UM21JnJjyI53hygZtXZbsrraFKE7PBLGj0wRZzhqCvvblT5BMpDNIZYeXUHkkIYQQUjkK2gghov85Mo3Hj80UPS5cYOtrUB5p0WXLI8PrM2hLpTO4MB9Bv724d+yG/jZEEmlxKIgrIGTaajCIJJvlC8aT4te7rLt00GY3aqCSyzA2H8GQM4hdG8x5H99qN2DUHQLHVT5BMje7RuWRhBBycfrxoUksrOMbtbVCQRshROSNJOCRCJ6ESX8rmR65FIVchhadct32tI17o0ikM5JB2zV9NihkDF4Z4Usk50IsVAoZzFpl0XNXy6TNlkfGUjgyuYB+u6Hs15HJGGxo1eLXI/NgUxlxCImgr82AYDyF+VDlf2+5u9nCVB5JCCEXnVl/DH/xs1N46pSz3ofSdChoI4QA4HeJ+SIJyTLFxUxb9YM2YH0v2BaGkEgFbQa1Avt6W8WgzZ1drM0wTNFzV0soj/THEvxS7R7pUf+5uiw6cdiIMIREsNW+/AmSuSWRtFybEEIuPqG40NdM7wHLRUEbIQQAEIglkc5wiCXTRaVrwsm1FoNIAH6C5Hotjxxxh8AwixMXC924rR1DziDcwThcgXhN+tkAwKBSgGGAoxN+BOOpsv1sgq5Wvq9No5Rhc1v+8QuTMEeXEbQJP0capYzesAkh5CIUZvmdrFF6D1g2CtoIIQDyR+4XBlDRbCmbTlX9njag8TJt7mAcPzs6XZXXGnGH0G3RQVvie3dDfxsA4JWRecyF2JqM+wf4ckejWoHXz3kAAPsqCdqy0yJ3dJggl+Vn/9qMahg1CozOhSo+hki2zLbdqBEnkhJS6MHXx/CvL4zW+zAIITUgZNqEtgtSOQraCCEAkNfLVhhAhWucabPo1Q0VtP3o4CT+9L9PVKVResQdliyNFGx3GGE3qfHK8DxfHmmsTdAG8CWSYTYFq16FHqtuyecLI/4Lh5AA/PTLre2GZZVHCtm1dqOaRv6Tkp48OYtfnCgeiEQIaX7CuT9CQduyUdBGCAGQn10rDKCEi21dDUb+A4DNoMJCNIFMpvJJhLU0sxADAMwGYqt6HTaVxpgnUjTuPxfDMLihvw0vD88hmkjDYa7+5EiBMEFyb09rRX1zvVY9AOASiaANALa2G1fU02Y3aag8kpQ0F2SxEE3W+zAIITUglEfGmqzaguM4pNKZuh4DBW2EEAD55ZGecP5EwEgiDZVcBpWiNqcMi16FDAf4Y41xoTbjjwIAnP74ql5nzBNBOsOVzbQBwA397WKpSC0WawuECZKV9LMBwECHEd/53b24Z88GyY/3tRvgCScqzkgKgVqbUQ02lan7GyBpPJkMB3cwDn80gXSD3MQhhFSPOIikyTJtP3xrAjf806/regwUtBFCAJQvj4ywqZrsaBMIC7a94cYY+z/j5zNszlVm2oZdpSdH5rquzwahZay9luWR2Uzbvt7KgjaGYXDbTkfJYL0vm0GsdMl2YWBKJZKkkC+aQCrDIcPxw5EIIeuLUHERa7KgbcQdxow/VtebSRS0EUIA8AGTRa+CSiErEbTVpjQSAKx6viTQ2wB9bekMJ2bYZgOry7SNusOQyxhsbtOXfZ5Zp8Tl2UXXtVisLX4drRIquQw7O6XLHZdLnCDprixoC7MpKOUMWnV88EjDSEghd3Dxd66R+lwJIdURFjNtzXX+F24isan6BZu1uwojhDQVbzgBm0EFtUJWtGA7zKagr1E/G7CYaWuEi7T5EItU9k6a07/KTJs7hF6rDmrF0lnKW3fYcWY2iA6zdlVfs5wPXtOLa/ts0CirkzXtNGuhVcorniApBP/CDQDqayOF5oKL2fZGOB8QQqpLCNaaLdMmtG/EkxnoVPU5BgraCCEA+D42q16dzbTllylGE+malkdaDdnyyAa4SBP62WRMNTJtIezoNFX03A9ftwl37+4suRqgGnZtMEtOglwpmYxBl0WL2QqDWyH4F36WwhS0kQKUaSNkfQs1a6Ytyp+P6plpo/JIQggAPmCyGlSS4/fDNS6PbM3etvI1wILt6ezkyO0O06p62mKJNCZ8UWxtL9/PJlDIZehsqV2WrVbajRq4g5X1Igq9kULWlvb0kEJuyrQRsq41a09bICfTVi8UtBFCAPCZNptBDateVZTxitS4PFKlkMGkURRl+OpBGEKyr7cVrkB8xWsIzs+HwXHANkdlQVuzajepMR+q7O+Nz9gulkdSpo0UcofiMGZ/PhaiFLQRst6IPW1NNojK3wA9bRS0EULAptIIxfmly1a9Km9nG7B4sV1LVoManga4sz7rj6FFp0RfuwHJNAfPCgPJxcmRpXe0rQftRg3mQnFw3NLBbZhNwUA9baSMuWAcGy066FXyovMQIaT5iZm2ZLphdrMuJZPhKNNGCGkMQhmS1aCGxaBCLJnOK13gL7Zr12sF8MNIGqE8cmYhhg0tWnEgyOwKd7WNzIWgksvQYy0/ObLZtRvVSKa5ipYhR9gUdCq52NPWbHt6SO25gyzsJv48RJk2Qtaf3AqLWLI53gNCbArCfcl4HY+ZgjZCiHhH22rgM23A4rJtjuP4i+0aZ9oselVD9LDM+IWgjd8lttIJkiOuEDa36aGUr+/TrLBzbS60dHAbYfmMrYEybaQEdzAOu1EDi664TJsQ0vzCbAry7GLSZulrDuTclGRTlGkjhNSRJ7vU2mZQLe5MywZyiXQGqQwnXmjXilQv3VrjOA4zCzF0tmjFoSArmSDpCbM4NuVf9/1sAN/TBqCiYSRCeaRWKQfDUNBG8qXSGXjC2UybXoUFCtoIWVc4jkM4nkKbgX/fiDbJBEmhNBKgTBshpM6EAM2WLY8EFksmhWZhfQ1H0QN8lm8hmiiqcX/i+Ayu/fuX8k6atRKMpRBJpLGxVYtWnRJqhWzZmbZUOoM/euQYYok0PvqOzTU60sZhN2YzbcGlg9togp9CyjAM9CpF0zWik9ryRhLIcEC7SYPWBsm8E0Kqh03xN4GFm33NkmnzxxbPRZRpI4TUlVAKac1Oj+QfE4I2/k5Y7csj1UhnOATj+cHZqyMezPhjePTQZE2/PgBMZ3e0bWjRgmEYdLZo4Vxmpu2fnhvGmxe8+Mq9l1R1J1qjEt5855aYIMmm0kimFzO2erWcMm0kj7CjzW7SwEpBGyHrjtDP1m5srkybP0qZNkJIg/CGE1ArZNCr5LAahPJI/iJcWIC5FuWRAOApGEYy5AwCAB46MI5Eje9wzWR3tG1o5UsjO8wazC5jV9vTp5z49qsX8IGruvHuvRtrcoyNRqOUw6hRLJlpE7JqumzGVq9SNN1yVVJbQomt3aRGq754IBIhpLkJ4/7bshUazZJpy630YSloI4TUkyecgM2gzpatyaFSyHLKI/mTbK1H/lv0+WWZAJBMZ3BuLoyBDhNcwTh+eWq2pscwmy2FFPrZOsxaOCucHnluLoT//dMTuKy7BX/5rp01O8ZGZDdplsy0Ff4c6dUKyrSRPIWZNgDw0QRJQtYNIdNmz1ZoNEuJfF7QRuWRhJB68oRZWLO9bAzDwJYzFCS8Rj1ti0Hb4sX/hfkIEukMPvqOTdjabsB3Xx2raB/YSs34Y9AoZeIFY2cLv4MslS5/kk5nOHz8h0ehVcnxH+/fC5Xi4jq1thvV4gV3KYUZW51K3jRv2GRtzAXjkDF81r1Vlz0fNMAaEEJIdYTiQnkkn2mLJZvjxl0gloRSzk+8pPJIQkhdeSOsGKgAgMWgWiyPXKNMmxA05k6QPOviSyN3dJrwkXdswqAziDfOe2t2DDN+fnIkw/An5w6zFhkOcC+RRTo3F8a5uTD+9+3b4MiuCriYrCTTZlBTeSTJ5w6ysBnUUMhl4vmAMm2ErB+Rgp62Zrlx548mYNGrIGMo00YIqTNvOCH2sgH8UJDC8sha97SJmbacO+uDziCUcgZb2gy4Z88G2AxqfPe1CzU7BmGxtqCjpbJdbSem/ACAfb2Wmh1bI2s3qjEXZMtmQYWMrbCkncojSSF3KC7u/RMzbZGlV0kQQpqDOIgkWx7ZLD2r/mgSLVoVNEo5ZdoIIfXDcVw2aFvMtOWWR65Vpk2tkMOoVuRn2pwh9LUboZTLoFHK8cGre/Dr4XmMuEM1OYYZfzwvaOs0V7ar7fi0HyaNApus+pocV6NrM6qRSGfKrmUQp5CqFqdHhpvkLitZG+4gK/a6CPsifZHar/oghKyNEJtfHtks1RaBWBJmrTIbtFGmjZC6+uVJ50U7XjrEppBIZ2DT52baVOLutkgif+pfLVkM+WO+h5xBDHQsLqj+wFU90Chl+M8aZNviyTQ8YXZFmbbjk37s7mqBTMZU/biagZAdKVciGS7I2OpViqYZ90xWp9Idi3PBONqzP0tGjQJyGUOZNkLWEWF6pFnL70FtlkxbIJaEObu7lU1Rpo2Qupn1x/DJR47iGy+N1vtQ6kIIznIzbRbD4rjtCJuCQsZAvQbDNSx6lbgzzhtmMRdiMeAwiR9v1avwnn1dePzYDP7yidN4cchdtRI7YXKkMO4fAEwaJYxqRdldbbFEGsPuEHZvbKnKcTQjoT+h3DCSqMT0yGgiXbRMnawvc8E4rvjbF/DikLvs8xKpDLyRhLisXSZj0KpTUqaNkHUkzCYhYwCNUsaXyDfJjbtALIkWyrQRUn/CHijHvFsAACAASURBVLDnTrtqOpmwUQkDR3J72oSsmzfCIsKmoFcrxOEctWTNyfCddfElkAMdprznPPDOPty0rR0/PTyND3/vMPZ8+Vf4wH8ehGuZS7ALzRSM+xd0tGjEgE7KmdkA0hkOe7ou3qBNzLQFS2dFhIytXuxpk2cfb443bbIyYx5+AuzhiYWyz5sPL+5oE1j0Ksq0EbKORNg0DNnrCa1S3jR72vzRpJgdpEwbIXUkBG2zgThOTAcq+pxJbxRPnZxdcqFwMxCWWedNj8z+tzecQJhN13zcf+7XFcojhb+X3PJIgK+F/87v7cPxL92KRz5yJX7nim68fs6DN857Sr7umCeCIWewbFZHzLQVBm1mbdlM2/HsEJJLu8xl/mTrm9BU7g6V/j6FsxlblZx/2xEybs3ypk1WxpU9R464yveh5u5oE7TqVFigTBsh60YonoJRowTA37iLNlBf81lXEOOeSNHjbCqNWDKNFp0S6jpn2mo7WYCQJjDkCsFmUCMQS+CZU86KMiaf/ekJHBr3AQA2t+lx9WYrbt/pwPX9bbU+3KoTyhHbjDl3uA2Li66jiVTNh5AIrAY1FqIJcByHIWcIbUZ1XgYwl1ohxzV9NuzaaMYP3pqAJ1z6jvwH/vMgZvwxWPQqXLnJgmu2WHHnJR2w5bz2zEIMMgZFI/s7WzQ4M1s6mD8+5ceGFq3YWH0x0qkUMKgV5TNtBRlbobctzKZgX5OjJPUgBGPDSwwPEm6Atedk2qwGFYaXCPYIIc0jzCbFc79WpUC0jpMYC3360ePobNHiwQ9dkfe40JNr1qmgUchoeiQh9TTkDOKy7hZc22fD06edS5ZIjrpDODTuw4eu6cUX7hpAj0WHnx+bwYceOoRQvPnuCgvliMKIbSC3PDKBMLuGQZtehWSaQzCewllXsKg0UopRrYBKIRMzhoXSGQ7OQAw39Lfhpm3tODkdwBefOIP3ffctpHMyb9P+GBwmDZTy/NNih1kLTzhRsiTixLQfuy/iLJug3aTG/BKDSHLXRghTJGnsf3N5fdSDrz0/UvHzXQH+Z2J6ISYOo5F+Hh+0OQozbdHmO6cSQqSF2RQMGmEYlVzsda63TIbDmCeC6YVo0ccC2XOQWctn2mhPGyF1EkukMe6JYKDDhDt3OTDli+HMbLDs5/z40BSUcgYPvLMPH71+Mx66fz++8b7LkOHQlHeFvWEWZq0SqpxBI0KmzRsWetrWrjwS4O+6j7rDGHAYl/gMgGEYtBnU8JQIGHyRBDIccPNAO776nt14/c9vwtffuxsj7jCeOjkrPm9mIVbUzwYAHdnMm1TPnDfMYsoXu6j72QTtRvUSg0jSeT9HYk9bA5XHkKU9/MY4/u3FUfgrXHqdWzJbblWHO8RCKWfybh5Z9SosRBN5N1cIIc0rHF+8CaxTNU5P21yIBZvKSLZCCJm2Fq2SMm2E1NOIO4QMBww4jLh1hwNyGYNnTjtLPj+eTOOxo9O4facjr7ROyAgJfVjNxBPJ39EG8HfAVAoZfJEEImwaetXaZNqEoO3t8QUk0pmKMm0AYDOoxEEGhYTsT1v274thGNyzewO22Y34txdHxQvC2UAsb3KkQAjkZiSGkZyY5vvZLubJkQK7SVN25H+koMxWyLpRpq25nMz+zB+dLD9YROAOLO4+LNfX5g7G0W7U5K3NaNWrwHGVrwwghDS2MJuCUQzaGmfty6SPz7CF4qmiigB/QaYtQZk2QurjrEsYdmGCRa/C1ZutePpU6SmSz5x2IhBL4n37u/Med5g0MGuVGHQ2X6bNE2LzdrQBfGAjLNiOJPLL2mpJWKh7IDtUZHvH0pk2ALAZ1CXLI4VgLrdnTyZj8KlbtuL8fARPnpjlSygLFmsLhEyb0198B+74VAAyBti1gcojhUxbqd+dMJvKC/7F8sgGedMmS3MH42Jgfni8sqDNFYxjX28rdCp52b62uSCb188GLN7EuVh3aBKy3uSWyTdSpk0I2gDAFci/QesXMm06yrQRUldDzhB0Kjm6LToAwJ2XODDmiZS8uHjk4CR6rTpcvcWa9zjDMBjoMDZlps0rkWkDFhddR9gUdGtVHpk9jjfOeaCSy7ClzVDR5/FBW/lMm61goMkdOx3Y7uCzbc5ADKkMJ5lp6zDzjzkDEpm2KT/67cY16/lrZO1GDdhUBsG4dBBWWGa7mGlrjDdtsrQT2UmpOpV8yRH+AMBxHOaCLBxmDbbajWXLx93BuLijTUBBGyHrSzi+2NPGZ9oa4/w/6V2cGllYIrlYHqmCWilDnDJthFTPkydmcde/vlbRLo1BZxDbHEaxJOe2HQ4wDPDMKVfRc0fdIbw9voD79ndL7iwb6DBh2BVqumXB3jArHbTp1dmetvSaDiIBgIVoEn3thqKhIKXYjHyAKfW990hk2gA+2/bpW7bigieC//vKeQDFO9oAQKuSo1WnxGzBiZzjOJyY9lM/W5aQJZkvMfa/8OdosaeNMm3N4uR0AHIZg3v2bMCJKf+SZUIL0SQS6QzsRg222Q3le9qC8bwdbQAFbYSsJ5kMh0giXZBpSzXEftxJX1RcR1MUtEUTYBjAqFFAo5CDpUwbIdXz8vAcBp1BvHHOW/Z5HMfhrDN/QmGbUY39vRbJvrZHDk1CKWfw7r0bJV9vwGFCLJnGhK94+lCjSqUzWIgmxbLEXDa9Cu4gi0Q6A8Ma9bRplHJxJ1ylpZEAn0VLZzgsSAxHmA+x0KnkkoHnbTscGOgw4UcHJwEAGyWCNiC7q62gp23CG4U/msRuCtoAQFx54C4x9r+wzFYojyw3UZA0lpMzAfTbjbiuzwY2lSm7CgPImQhp1qDfboQnnJDMiMcSaQTjKbSbKNNGyHollMKLQZtajgyHuk5jFEz6orhkI9/mUDh0zB9LwqRRQiZjoFHKKdNGSDUJd3PLDRQB+GXawXiqaELhXZd0YMQdxrm5xbvC8WQaPzs6g9t3OkruDWvGYSS+bJBjk8y0qcTFuLo1LP8TSiR3VDiEBFgsfZTqa5sPsUVZNoFMxuBTN2+FcKNPqjwS4He1Fd59E4aQUKaNJ2RJ5kpm2vIHkchlDLRKecM0opPyOI7DyWk/Lt1gxr7eVgDAkSVKJHMXZm938L/PUtk24WfGXhC0CZMkpW7GEEKai3CDTiyPVPI3aBuhRHLSF8PWdgOsepVkeaRZyy8EVytkSGc4pNL1CdwoaCPrSjrDYdQdBgA8P+gu+4t11rk4hCTXHbscAIA/+ckJ/N0zQ3j82DQePDDGDyC5srvodQRb7QbIZUxTBW3CjjapQNSSE8gZ1qinDeDLMgGIF3mVWAzaiu/iz4fYon62XLfvtGNHhwlWvUrM/hTqMGsxW5BpOzbph1Ypx9b2yvru1jshSyKVaWNTaSTTnJhFFejVCoSpp60pTC/E4I8mcWmXGXaTBhtbtUsOIxGCNodZg34H/3si1dcm/MwUlkcKmXdviSFDhJDmEY4XZtr4f9f7xl2ETcETZtFt1aGjRVPUv+6PJtGi44M2TTbQrFe2jbrnyboy4Y2ATWVwy0A7Xhiaw8ExH67ts0k+VwiuthVk2uwmDT5/13Y8fmwWD74+hmSaT8Nssulx9WZr0esINEo5Ntv0GGqiCZJi0KYvzrTlTpRcy0EbwrEMLKM8ss3If45k0BZm0VdmoAnDMPj3919eVP6Yq6NFg2A8lZctOjHtxyUbzFBU2He33hnUCuhUcsxJBG3CsJHCnyO9Wk49bU1CyCxfuoHPLF/Ra8Frox5wHCfZ4wtAzNS3GdTZHWxKyUxbbkaukMWgokwbIetAqDDTpmqMTNtUdqF2t0UHh0lbtGA7L9Om5N/v48n0mk3VzlWVr8gwzDiAEIA0gBTHcfuq8bqELJdwQfDRd2zGgXNePHPaWTpoc4XQZdHCqFEWfexj12/Bx67fgmQ6gzFPBGddIWyzG0tenAi2d5hwtIKpao3CG+EvsCUzbTmB3FrtaQP4E2eXRVuyDFWKkEmbl9gT5gmzZYNtgA/IN9n0JT8urAL4jX97DQMdJmx3mHBmNogPXt1T8TFeDPhdbcXlkUJgVhS0NdCeHlLeyekAVHKZeJNrb08rHj82g0lfFD1W6d8ddzAOm0EFlYK/0OkvMUFSDNqMEkGbjl89QghpbkKmTdjTJlxX1Dtom/QuBm0dZg3eHvflfTwQS2JjtnVCo+ADzXr14VXzFvFNHMftoYCN1NOwKwyGAS7d2IIbt7XhuTPuktMch5xBDCxRgqeUy9BvN+I3d3cWZeSkDHQYMeOPSS6DbYQJSYWEHrC2Jcoj1zLT9tnbt+F/Pn7Nsj7HrFVCKWeKetrYVBr+aLJkT1ulbhmw409v7cc2B7/W4V9eHEEilcGVm8oHgxebNqNaOtNW0IAuMKgVNRlEkkhl4K9zdiad4Yoa2pvZyWk/BjqMYgAm9LWVK5F0BeJ52bPtDiNG3OGic+FciIVaIYNJW3yesehVWKCgjZCmV3jzTitk2upcbSHsaOux6OEwaxCIJfNuJvqjCbE8MjfTVg9UHknWlWF3ED0WHbQqOe7Y5cAzp104MrmAK3otec+LJdIY90Twrks7q/r1hf64YVcI+zctfs3z82Hc880DaNUrsd1hwoDDiIEOE27a3i7WSNeDN8xCIWMkL5byyyPX7hgNasWyyw4YhoFVX7yrTSj/LNfTVgm9WoE/vnmr+P/RRArOQByby2TnLkZ2kwansmV0uYQ3a11BT5tOLa/JZMC/evIMnjg2g59+/Brs6Ky8N7KafnxoEn/7y0Ec/eKtJXslm0Umw+H0TBD3XrZBfKy/3QijRoHDEwv4/0pM1HUHWXE5PQD0O4wIsynM+GPY2KrLeR4f3ElVMrTqVRjJ9ikTQpqXWB7ZaJk2XxQmjQJmnRKdLfz5yhmIY0ubAZkMh0AsiRYtfxNbLWTaks2daeMA/IphmCMMw3ysSq9JmsAP3prAY0em630YomFXCP12PiP2zu3tUMllePpU8RTJEXcIGQ7YsYy+qUoImbvCYSQ/fGsCbCqNSze04MJ8GN98+Rz+8EdH8Uh21LyUs64gRsvsNaoGb5hfrC11sZSXaWuCi06bUVUUtAnlkqvNtBXSqRTY0mZYslz2YtNuVMMdZIsyKcKwkcJgXF+DTFs8mcaTx2cRSaTx4e+9jblgfbJdRycXEE9mECqxbLyZXPBEEGZTuDQ7EhvgJ69e3t2KIxO+kp/nDsbzxvhvy56bC/vapHa0Cax6FY38J2QdEMsjNfmZtkidS+QnfVF0W/mbSA4TXwYpVEmEEylkOIg9bRoh01bBHuBaqFbQdi3HcZcDuBPAJxmGub7wCQzDfIxhmMMMwxyen5+v0pcl9fatl8/h4TfGS3780JgPn3vs5JqUBsaTaYx7o2IZo1GjxDu22vDcaVfR1z/r4oOq5UworITdpEarTpkXtMWTaTx2ZBp37OrAv7//crz4mRsx+OU7YNYqccFT+g7ynz92Cn/22MmqHl8hT5iV3NEGAHqVXCyFWsvyyJVqMxRn2kot1ia10W5UI5ZMFwVipXvaqj+I5NfD8wixKXzuzu0IxJL48PcOV6Vv7vtvjuOVkcrfu4SBRPW+i1wNJ4UhJBvz11vs62nFiDssWYrKptLwRhJw5ARtW7NB27Br8bwXjCcx6g7DYZZet9GqVyGWTCO2Dr6PhFzMwgXvA0IFT71/tye9UfRY+KoZoTJAGPsfiPKtLmahPDKbaatXeWRVgjaO42az/54D8DiA/RLP+Q7Hcfs4jtvX1tZWjS9L6swXScAZiGPCGyn5nKdOzuLRt6fEtHgtXZiPIJ3h8nrP7rykA7OBOE5M5y+BHXKGoFfJ0W3RFb7MqjAMg4EOU17Q9suTTgTjKdy3v0t8TKOUo8uixZSv9MTCcU8Eg7PBmu4D8UT4TJsUhmFgyw4jqceUpOWyGdTwhPIvHmuVaSPShP6luYKBMBFWuqdNr1YgWuWR/0+enIVFr8JHrtuEb9x3Gc7MBvDpR4+X7G2tBMdx+Mdnh/Ffr49V9PxkOiPueVwPg1ZOTgegVcrRV7DeYm+2r+3oZHFfm/C75zAv/u6ZtUp0mjV5mba/euIM/LEkfv/aXsmvbcnuavPRBElCmlqETUGtkEGZnbisU/LvB5E6Bm3pDIfphRi6steCjmzQ5sqO/RfmExRm2pp2EAnDMHqGYYzCfwO4DcDp1b4uaXyDs3xgEoynxLsRhcazU3mkpvqt1N88NYiv/Wq46PFhd3aEv30xaLt1wA6FjClatD3kDGKbwwiZrPrlbQMdJgy7Q0hnLxJ/fGgSmyXWBXS16sRRs4UC0SQCsSTYVAbn50sHxas1H4yX7feyGFSQMYsnqkZmM6rhjeSX5gk/d1IrDUj1tWeDY3dBSWKpTJtBrUAkkapaJj7CpvDikBt3XeKAQi7DzQN2fPFdO/CrQTf+4dmzK37d+RCLMJvCiMTkQykX5iPiqpBq35FdTfC5lHgyjX9+bhhPHJ/Je/zktB+7NpggLzhf7ulqgVzGSA4jKTXGv99hxNns9/Gpk7P42bEZPHBTHy7rbpU8JmGKrY92tRHS1EJsSiyNBBbLI2N1vLHlDsaRSGfEG/gapRyWnAXb/uy1bYs2f08b28SZNjuA1xmGOQHgEIBfchz3bBVelzS4Qedi9mrCJx1YCFk4T5WCNn80ge+/OY4HD4yDLagpHnaFoZQz6M0ZDmHWKXH1FiuePe2CJ8xf0HMchyFnENs7ajOgYLvDiHgyg3FvBCPuEA5PLOC+/d1F/U9dFh2mF2KSF2G5wdzpmUDRx6vh2OQCZgNx7OlqKfkci14NvUrRFL1bNoMayTSXN7lzPszCpFHUddjLxUToXyq8SSPcSS0aRKJSIMMBsSq9Ab4w5EY8mcFv7l4cmHH/tZtw3/5ufPvVC0UL0it1wcOfx1zBeMkbVLmE8mugeuWRk94o7v7G6/i9Bw/VpNz83FwYv/XvB/DNl8/hU48ex48OTgDgs4ZnZoNFpZEA//e3q9OEwxJrTlwBYWF2ftC2zW7E+bkwphei+MLjp7G7qwUPvLOv5HGJQRtl2ghpauF4Kq/aQqWQQSlnapJpy2Q4PHvaJd48L0WYHJlbdeUwaRbLI7PXEy06YRBJk2faOI67wHHc7uw/OzmO+0o1Dow0vsHZIIQbrxPe4oxRMp3B9AJ/kTQvsfR4JZ497UIyzSHMpnDgnCfvYyPuELa0GcTUu+A3LunAhDeKfX/7AnZ+6Tnc8rVXEIynxEmP1Sa87pAziEcOTkIll0lOV+tq1SKRykh+b4QTCQCcnq1N0PbQgXEY1YqSk98AoMOkQWuTZKlshuIF254wS6WRa6jdJJ1pC7MpKGSM+IYnMGR7GiJVKpF88oQTDpMG+3ryszbvv7IbAHBwzLui172Qk+0emVs62zaYUx5djX6NXw/P4e5vvo4hZxCvn/PgpbNzq37NXP9zZBp3f+N1zIVYfPt39+Lm7e34wuOn8YO3JjDqDoNNZfKGkOTa22PBiSk/EgUXMcJibUdhps1uRCKdwf0PvY1EKoN/ee+eonN2LjFoi1SvWoMQsvbCbEpcrC3QqRQ16Wl77owLH//hEbw45C77PGFHW491MWjrMC8Gbf4Yf7PIXJBpa+qeNnJxGnQGxbH2uUGGYGYhJt7lqFam7cmTs+i26GBUK/DsaVfex4ZdIcldau/euxEPfmgf/vJdO/A7V3Rjc5sB+3stuLG/Nr2VW+0GyGUMjk368bOj07hjlyNvUbVgY/bOzpTE9074fvbbDTgzEyz6+Go5AzE8fcqJ91zRVbZf7TO39ePbv7u36l+/FtrEBduLd+TnQxS0rSWjWgGNUla0qy3KpqBXF2dshXLJagwjCUSTeGVkDu+6tKOo7HmgwwSTRoGDF0pPOixnzMPvfwQguRy60FlnSMwqriaLmMlw+OZLo7j/4bfRYdbg2U9fjx6rDv/03HBVyiTDbAp/+pPj+OxPT2B3lxlP//E7cPtOB771gctxy0A7vvjz0/jK04MAioeQCK7obQWbyhTdXJoLxqFSyMT9RgLhHD06F8YX37Wj7FJ7IDdoWzrDSQhpXIWZNoCvvqj2MCoAeH6QD9aOThavoMk16YtCLmPyVpM4zJqinjZxT5tC2NNWn0xb408XIA0pnkzj/HwEt+904NxcRLxbkWs8Z0BJ4dLjlZgLxfHmeS8euKkPk74onh90I5XOQCGXIRRPYsYfw/vs3UWfp5DL8M7t9lV//UqpFXJsadPjRwcnEE9mcN/+4mMC+J42gD9p7CvYIzfpi8KiV+GqzVb87OgMMhmuqv13P3hzAhmOw4eu6S37vHaTJm9kdyOzZYOz3MzlfIjFJSUuNkn1MQwDu0lTNIgkzKYlbw4I+8uqMfb/uUE+C3/37uLdi3IZg/2brHjrwsozbdvsRkwvxIrG1Us56wpi98YWvHnBu6q7yJ9//BQefXsK9+zpxN//9qXQquT401v78alHj+OpU078psSftVJnZgP4o0eOYdwbwadv2Yo/eudWsWdNrZDjW+/fi0/86CheGHLDqFGg1yo9tEk4dx0a8+HynL40V3aMf2Gg3tdugFohw3V9trzhTKWYNErIZQxl2ghpcmE2Je5BE+hUckSrnLVKpTN4aZivRjg+VVy6nWvSF8WGFi0UOdn+zhYtFqJJxJNpBKJJqBQyMcMm9rQ1+ch/cpEZyQ7a2NFhQo9VJ9nTJpRMqhSyqgwieeaUCxkOuHt3J+7Y5cBCNIlDY77s8fAjpHOHkNTTQIcJ8WQGm216XLXZIvmcja38iGupCZJTvii6LDrs7DQhzKYwIZGNW6lYIo1HDk3i1h12cWLSeiAMVMnN6s6HWLFskqwNu1EDZyD/ZzrCpiQXtAuBXDX6vp48MYseq65kGd9Vmy0Y90bF/TvLMeaJYHObHv12w5KZNl8kAXeQxeU9/M2Clf7Z4sk0Hjs6jffs24h/ee8esWn/7ks7sd1hxNd+NYxkwWRZb5jFa6Pl1xJwHIfvvzmOe7/1BiKJFH70kavw6Vv6i4aMqBQyfOv9l+M9+zbid67oKtnX2mZUY7NNj7fH8rOYrkC8qDQS4C96fvnH1+Hf3395Rb2yMhmDVp2SMm1kzY24QzUd/HOxCbNSmTYFolXOtB2eWIA/msSGFi1OTQfK9rVN+KJ5pZHAYkm3KxCHP5oUh5AA9c+0UdBGVkSYHLmj04Rui04y0zbhjUKnkmNLm6Fof9ZKPHliFtsdRmy1G3F9fxs0ShmePcOXSAp3v6XKI+tB6GuTGkAi0CjlaDeqJSdITvqi6GrVYmcnfwFazWEkPz8+A380id+/dlPVXrMRtGj5O/LCz1o0kUIkkabyyDW2vcOIMwWrKiKJlJhVy6UXe9pW96btCbM4cM6Duy/tLPn7dlV2euty+9qS6QwmfVFssumxzWHEiDtUdhDI2Ww/m5B1Wml55JnZIJJpDjcP2PP+TDIZg8/ctg3j3igeOzItPn5scgG/8W+v43f/6xDeOO+RekkkUhl84kdH8ZdPnME1W6x4+o/fgau3WCWfC/CB2z++eze+8Bs7yh7r/k0WvD3uy7vAnQuxRUNIBH3txmUNB7LoVZRpI2vqwdfHcNvXX8WPDk3W+1CazslpP9777TeL1p1I97TJq77L8vlBN1RyGf7wxi2IJNIYLdOHLNwgzyWUSs4GYgjEknkl3gq5DAoZQ5k20lwGnUEY1Qp0terQbdHBGYwX/RBPeCPotujQblSvehDJjD+GwxMLYumTTqXAjf3tePa0C5kMh2EXv3dtQ4v0gta1dusOO27a1ob/ta/0kA+AnyBZ2NOWSmcwsxBDt0WHfrsRSjlTtWEkHMfhwdfHsKPDJPYjrhcyGQOrXiUGbcLOtrYyKw1I9e3taUU0kRbHugPSd1iBxZ625ZRHchyH770xji8/OYhvv3Iejx+bxrdfOS9m4UsZ6DDBqFEsu0Ry0hdFKsNhs82AfrsRC9Fk2fPZUPbPvburBQyz8kEkx7K7zy7rLi7vvWWgHXu6WvCvL44inkzjkYOTeO+334JSwcBuUuMfnx2WDCz/6/UxPHPahT+7Yxse/OAVsFbpd+OKXguC8RSGszfPOI4rmWlbiVadCguUaSNr5OfHZvDlp/hezpeWGGRBij172oWDYz6cmc3vx+d72vJ7XKsdtHEchxeG3Limz4pr+2wAgOMl+tpC8SR8kUTRvt7FXW1x+GMJcQiJQKOUU6aNNJfB2SAGOkyQyRj0WHXgOIiTIgXj3gh6rfrs0uPVBW2/PDkLgC8NEtyxy4G5EItjU34Mu0LYaq/N3rWV2NJmwEP37xfHxJbSnR37n8sZiCOV4dBt0UGlkGGbwyhmNlfr9XMejM6F8fvXbWqKMf7LZTOoxf7J+TBfBkeZtrUlZJhyFy5H2bRkeaReLI+sLGhLZzh8/vFT+NIvzuCRQxP4u2fO4k9+cgLffW0M2x3Gspl2uYzB/l7LsoeRjGUnR25q04vl1+VKJIecQbQZ1bAZ1NApV35BcmzSj42tWrQbiwMfhmHwZ7dvgzMQx73fegOff/wUrtpixZMPXIc/uaUfx6f8+NVg/sWmMxDDN14axS0Ddnzixr6qniuFG0Bvj/Pf22A8hVgyXTLTtlxWgwpeyrSRNfDr4Tl89qcncNVmC+7b34W3LvhKZlXiyXRRiTKBGKzlnifZVBqJdCZvTxsA6NSKis7/X/3VMD796LElnzc6F8aEN4pbBuzotepg1ipxfEo6aBMGvvUUZdr4m//OQByBWApmbf51nFoho+mRpHlkMvyesx2dfAmgUA+cWyKZznCY8sXQY9OhzchfSK9mt9CTJ5zY3dWC7pza43cOtEMpZ/DsaSdG3CFsb5DSyOXoatXCGYjlnfinCvaG7Oww4/RMoCq7mR46MA6bQYW7d3es+rUakc2oFjNtQh9lueXhpPr4QEONIzm7u8LZ6ZGFDOIgkqXf8ovQawAAIABJREFUAJPpDD716DH8+NAUHripD0NfvgOn//p2vPiZG/DIR67Egx+6YsnXuGqzFRc8kaKVBBzH4aWzbsmLswsevl92s02PfsfSQdtZV1A8F2lV8hWXRx6dXCi5cBoArumz4do+K4acQTxwUx8e+tAVaNGp8O69G7HZpsc/Pzec18vx/z99FqkMh798V/lSx5XY2KqFw6QRe4zFxdrm6mXavBHa00Zq69jkAv7wh0exzWHEd39vH24ZsCOWTEsuj+c4Du/59pv4wuOn6nCkjU1YeTKaM7RJWOuiL9zVWcGNrUyGw48PTeHg2NI33ISpkbdky8p3d7WUDNqEa63C8kitSo4WnRKuQByBaKJoAq5GKW/ePW3k4jPhiyKSSGNHtm+r28KPbJ7ImRbpDMSQSGeymTYVEukMgrGV9a2MeSI4NRPA3ZfmBxomjRLX9tnws6Mz8EYS6G+QISTLsdGiQ4ZD3tLfyYITya4NJixEk5hdwQCFXGOeCF46O4f3X9kDtWJ9Lpu2GVRiVlcI2top07amGIbB3p7WvKAtkkhBL9HTpquwpy2eTOMPfnAET5104i/u3I7P3r4NDMPAoFZgS5sB1/TZ0FlBafSV2aFAhSWST5504vcfPozHjswUfc6YJwKLXoUWnQo2gxpWvarkBMlUOoMRd1jsadWq5IhVmEXM5QzE4AzEcblEaWSub953OZ584Dp89vZt4iARhVyGz96+DaNzYTx+jP/zvHneiydPzOLjN2zJu/FVLQzDYP8mCw6N+cBxnBi0Vas8cmOrDv5oEsE4lUiS2pj1x3D/w2+j3aTGw/fvh1GjxFWbrVDKGbwqMdxnyBnCyelARYHExWQuFBffe4dzzpPhOH8eNGjyAyC9WrFk0DboDMITZuGPLv37//ygG5duNIsljnu6WjDiDkm+x4iLtSXOifyC7Rj8sWRReSRl2khTyR1CAvAXyjqVPG/C4UTOwsI2iVHs5bgCcUz5ogizKXAchydPzIJhgHddWtyvcucuh3gHtlGGkCyHMPY/d4LkpC8KRc7ekJ0bqjOM5AdvTkApZ/D+q6RXEKwHbYbFrO58OAGGgeSOPFJbe3taMb0QEy/eIyUybUq5DCqFDJEygQ3HcfjYD47g5eE5fOXeXfiDG7as+Lh2dJhgVCvyLrQSqQy++qthAMABiQEe5+cj2JyzS2ybw4jh7LTaQmOeCBKpjJhp0ykVK8q0CT0Y5TJtANCqV+ESiWmZd+5y4JINZnz9+RFE2BS+9IvT2NiqxSduXPn3bilXbLJgLsRi0rc4obNaQZuwy23cUzylmKwOx3FU4gd+8EgonsLD9+8Xr1n0agX29rTi1ZHi88ITJ/gbIhPeKELr/GZCIJbEhfkw/NHEktM0hdLITTa9ONUbAEIs/z0q7G3WquRLlke+fJYf3x9LpssGS3OhOI5P+XHrwOKKp8u6WpDhgJPTxddPk74oWnVKmAoCSYAfRjLpiyKaSOdNjwQANWXaSDMZdAagkDHoazcA4O+ydhcM1BB2tPVa9TlLj5cO2uZCcVz/Ty/jHf/4MnZ96Tn0/59n8I2XRnFFr0W8c5LrlgE7hNaMZsy0dVmyY/9zJkhO+qLY0Lq4N2TAYYKMQVFT73JE2BR+emQKd+7qkOyRWS9sBjWf1Y2nMB9iYdWr8vavkLWxtyfb1zaxADaVRjLNwSDR0wbwb+LlMm3BWAqvjszjEzduwfuv7FnVcSnkMlyxyZKXaXv07UlMePmxz2+d9xaVIY95InkLoPvtRoyWGAUuDCERMm2aFTbZH51cgEohE6sZlothGPz5Hdv53ZXffQsjbn6R9XImNi7X/px9bUKw3m6qTpZb+P6PUdBWVRzH4c8fO4mb/vnXZceir3cRNoWfHJ7CnbscRcver+9vw5AziLnQYqVLJsPhqRNOmLL9WUPOpXc3Vnoct37tlRXvk6yV+77zFt751Vew58vPo+8LT+Pyv3ke/+fn0mWhwk39e/Z0whdJiO0KQqatsKdNr5IjmeaQKBME/XpkMdMZjJUOkF8c4oO7W3YsBm27u/hqBakSyQlvtGgIiaCjRYsL2X5mc1F5ZHUybZ4wi3949uyyPoeuZsiyDc4G0dduyLsA6LboxOwawP8yqBQyOEwacelxJWP/nzvtQiKVwRfuGsDn79qOD1+3Ge/euxF/dvs2yedbDWpcuckKm0HVlAMnOsxaKGRMXsA75cs/kWizaxPOrCLT9vPjMwjFU/jgNau76G10NiOfVfOE2eyOtub7mVgPdnaaoVLIcGRiAVGhl0Ei0wbw08MiZXraPNkBFFvbq3NT5spNFlyYj2AuFEeETeHfXhzFlZss+ORNffBGEvl3h+NJzIdYbG4ziI9tcxgRTaQx4y/er3jWGYRCxmBL9vk6pXxFb+7HJv24ZAP/PVyp67byPW8npgO4vr8Nt+VcyNTC1nYDzFolDo354ArG0aJTVi1I7LHqwDDAuKd6+yrXA47j8PLw3IovIB86MI7/PjyN6YWYOK30YvSzo9MIxVO4X2INzvVb2wAAr+Vk245MLmDGH8Mnb+oDAAxWabrzhfkIRufCeLGBJlZemA9j0BnEe/d14Yvv2oFP3tSHTTY9fvL2lGSGbHA2iG6LDvt6+Js4I9kbWUI1RXGmjf//UlN2/dEEjk0uYEsbH0z7ywRtzw+6sbFVmzffwKJXoceqk1yyLTXuX9Bh0iCVvZEhVR7JrmJ6pCsQx5efHMR1//AS/u8r55f1uRS0kWUbdAaL7gD3WHWY9EXFu8/jngh6LDrIZMyyMm1Pn3Khr92Aj16/GR+7fgs+d+d2/N1vX4p9vaXH0//tvbvwzfddvoo/Uf3IZQw6W7SYWsgvjyw8kezaYF7x2H+O4/D9Nyaws9MkTvZbr3IXbM+H2aYM5NcDlUKG3RvNODK5II7zLxW0LZVp82angVYrABf3tV3w4T9fG4MnnMDn7tyOq7OPv5lTIilkdgozbYD0MJKzrhD62g1isKVdQaYtkcrg5ExgyX62Snzhrh3Yv8mCv/7NnTWfFiuTMbiil9/X5g6yVSuNBPjG/06zFmMe6bLUavOEWTzwyFHc9vVXcKZKF+S18OCBcdz/0Nv44VsTy/7cN8578JWnh3BDfxsUMgYvZLMU65UzEMMrI8W9aZkMh4ffGMelG82Sv3M7OkywGVR5fW1PHJ+BRinDB67qgVWvWlUVTC7hRlC1Xq8ahMEef3RzHz583SZ85rZt+NTNW5FMc+LgoVxnZgPY2WlCv4O/cSX0tYXi0u8DwmCSaFL6PeDVUQ8yHHDvZRsAoGRfWzSRwuvnPOIAklx7JIaRnJz2Y9wbFdt8CuVWdhVOAecHkazsRsm/vjCK6//xZXzvzXH8xiWdeP5PbljW51PQRpbFE2bhDrJFP+jdFh3YVAZz2cCMLzfiL3TMWiUUOUuPy732wTEv7trlWNYxbWkziBdizajLohUzbaF4EgvRZFHKfmenCe4gW1HgW+jQmA/D7hA+eHXvuhzzn0sM2sIJeEIs7Wiro8t7WnF6JgBftudUahAJwL+Jl+tpE84bVkN1ehN3dppgUCvw9CknvvPqedy5y4HLulvRZdGhy6LFmzmlSUJ5jHCXFwD67fkXI7mGnEGxNBJY2fTIIWcQiVRmyX62SuzoNOG//+DqopKvWtm/qRXj3ijOzASqNu5f0GvTYcxbOtN2bq780vNKCD3Ut37tFfzqjBu+SBK//a038NPDU6t63UoFokmckui9kXJ43Ie/e3oIAPDS2eUFXNMLUTzwyLH/x959h8d1lmkDv9/pM5omadSLZdmWJXc7LkkcO52EsiTZQAJZINTQCaEtLOwHu8vuAksJS89CaIGEXUoghRSHJI7tJI4T9ybLsi1bvWt6Pd8f55zRjKaojaSRfP+uS5flkTQaycdn5jlPQ12xBd+/Yz221BflVXZnJnzrqWbcef9ePHGkK+n2F1r6cLrXi/dsTf/8qNEIXLHUhRdO9SEWk/v/HjvUietXlKPAqMOKSnvOgix1IFmupkXnwlPHurGy0o7qwtHXJJvqimDQarDndHIZpzsQlgOhCjtKrEY4Lfp45YJ68W5seaTZoA6jSn+efO5kDwotemxvkDOeQ770U2T3tPQjFInh+jQVBetqnOgeCaJzWP79SpKErz56HC6rAe+8NH31kTr2H0ifaZvKnrZAOIrvPtOMS5cU47nPXIVv3bY23mY0UQzaaFKOdyYPIVHVKgFa24APkiTh3IAXdcpEHo1GoNhqGDfgePJoF2IS8IY1C3McfSY1hRZcUHra1IEkqUGbPGxgKld9f/XiOTgterx5XebFwwvFaNDGTNtcu6S2EOGoFO/PSLenDRi/PLLfk9vVDTqtBhvrCvHXI10IRGL4TELp9WX1xXipdSBeMdDa54VGJE8Xs5n0qHKaUyZIDvlC6BwOJJXmmPXaSS/XVvfbzces+CalIqJjOICyHPWzqRa7CnCm15P2xeyuU3247ts78ZHfvJZ2Ufv5AR8++dB+vOfne3H3Q/vxzw8fwTefPIn7dp7GQ3vb8NfDndjZ3IsPPfAqPv7gftQWF+CxT1yBJz65DRtqC/HZ3x/CF/54eMYnxt39u/24+Ye7kyavptPnCeKjv30NVYVm3LGlFnvPDEx4GEYgHMWHHngV4UgM971rI2wmPa5tLFP2Wy3MnkFJkrBTybJ95v8OoqVnNGP7i91n4LIa8YbVmV93bG8owYA3hGOdI9jV0odBXxg3rZWfT1dU2nGqx521J2ui1KBtJBBJ2d86F3rdQbzWNojXrUi+kG42aLFhkRO7TiUPaFF7+1ZW2SGEQEOZLX6ejE+PTMm0ZS6PjMXkf7dty0riA8UyZdrUQXgr02TO1ql9bcqApyePdmHv2QHcc30DbGmGkABAhTMh05ZmufZUMm2tvV7EJOCtl1RnLMscD4M2mpT45Mix5ZHKAXiu34sedxCBcAyLEq7uliTsz8rk8cOdqHeNLrC9WNQUWdDnCcEXioyOoB3zH1oNkid7Ra9rOIAnjnbh9o01MzqEIF8UFRigEXIdfigSY0/bHNqgDCN5QXliH/tkrRqvPFKdAlpoSf/kOhVqZv72TTXx/jMAuGxJMYb94fieodZeD6oLLSkrMhrKrCnlkSeUvzcmnBstUyiP3N82hAqHKe3gpXy3qsoBs3KeyWV5JCAPtRoJRDCY5kWbGuQ8ebQLN/9gN073yi/KYzEJv37pHG64dyeePtaNfm8IB84P4dFDHfjR86fxH4+fwOf/eBgf/s1reNf9e/HsyV58/vWN+MOHLsOyMhtcViN+/b7N+PBVS/Dg3jbc9pMX4xfYJkKSpAlnTA6eH8JzJ+XA4lP/eyDj/4loTMLdD+3HkC+MH/7DBty0thKRmJTyAjqTf3v0GI60j+Det62LH/vXKdP20pVISpKEfWcHcOjCUNJU5/nkZLcbPe4g7rmuAQadBh964FV4ghG09nrw7MlevOPS2qxrcLYpfW3PN/fiLwc64DCPZn5WVNgRjko41TP9YSSdw4H46o6JXqB99mQPHj3UMe3vnc4zx7shScDrVqZmr7YuceFY50i8kgIY7e1TLzIvL7OhuUvOgHuDEQghnxMTqX9PV21xpGMYfZ4Qrm4siZcoDvnTZ9oGvSFoNSLtJMgVlXYYtBocOD+EUCSG//zrCTSUWXH7xpqMP3vi+Wtsps2k004p09ainJeWlU0uu5Yo/bMoUQbHOkdQ5TSn1PhWFZqh1Qi0DfjiY5nrEq5Ou6zGrCP/+z1BvNQ6gA9fuWTBl/CNVV0op+EvDPozLnt0mPVYVGzB/rYheIMRWAzaCf2efvvyOcQkCe/IUAKw0Gg1AkUFxvgUP2ba5o7LakRdsQV7z8p9D5l62grG7WkLotCS2ymgb1xdgX1nB/DJ65Yl3X5ZvQuAvMdtVZUjZXKkqqHcht0t/QhHY9Arj0utQmgam2nLkp2RJCnl//H+84PzMssGyCscNixyYndLf84Wa6sSJ0iOXeNxtGMY9SUF+OpNq/CxB/fjpu/vxhff2IQ/H2jHS60D2LbMhf/8+9VJJV6SJMEbimLEL+9/G/aFUV1kQdWYfX86rQb/eGMj1tU48Zn/PYg3fW8X/vtt6+Mv2jMJhKO4+pvP4f3b6vG+K1IHXIz1vb+1wGHW497b1+G9v3wF//boMXzt1jUpn/edp5uxu6Uf37h1DVZWOhCJxmA36fC3Ez14fZZsESBfhHhwbxvefXkdrk0Yi15bbEFDmRXPHO9Oeax/eK0dn/m/g0m3mfVa3P/uTbhsyfxoS1CHiNy2qRqb6grxjp+9jM/9/iBKrEbotQJ3bMm+BqfEZsSKCjueOtaNU91u3LSuMt63qgYoxzpG4u9PVfuQHxtqnXitbQhHO0Zw46rs/57PHO/GXb9+FU6zHm9cXTGl107HO0dw90P78dkbGlNKC59KM9hDdflSF771dDP2nO6Lr2M62jECl9UQ343aUGaFOxhB53AA7mAEVqMu5TFajJkzbc+d7IUQ8jCYAoMWOo3ImGnr94ZQaNFDo0n9HRh1WjRV2rH//BB+9eJZnOv34Rfv2ZT1OaXAqIPdpMNIIAJ7ysh/DQJTyLS1dLuhEZhWuTqDNpqUox3JPRsqvVaDSqcJ5/p9qClUdrQVJWTarEacyDIW9+lj3YjGJLx+9eT62RYCNUBr6/ehbcAHh1mfcmUHAFZXOfDooU6s/PKTEEIuKygwalFg1MFq1Cl/18GacNvvX72AaxtLp5yKn49cVkP8BTSDtrm1YVEhzr4m7zPKlGkrMGjhzZKN6vME4cpRP5uqpsiCn965KeX2cocJ9a4CvHi6H++7YjHO9HmxeXHqEKTlZTaEojGc6/diaam8AuB/draiymlOOubMBi1CkRiiMSl+BV016A3hhnt34sNXLYlPret1B3F+wI87L6vL6c87mzbVFclBW45XiyTualNXSqiOdoxgw6JCXL7UhUc+fgU+/MCr+MIfD8Nm1OFrf78at2+qSXmxqC5ntxp1qMT4i9lvWFmOho/b8OEHXsWdP9+Le65rwMeuXpr2RSIgjx/vHA7gNy+fw3sz9EupjrQPY8fxbnzq+gZc3ViKD125BD967jSubizFDSvl50RvMILvPnMK9+1sxW0bq3HbJjlLoNNqsL2hBM+e7EUsJmV8PADwnR2nYNRp41MPE13bVIb/2dmK4YRlwsFIFN95uhmrquy4+9oGDPpCGPSG8O2nm/HM8e68Ctp+uecsyuwm3JimJ37nqV4sK7WiwmFGhcOMz93YiK/99QSEAG5ZVzWhNTjbGlz4yfOtAIA3r62K377YVQCzXoujHSN46zR/ho4hP65eXooRf2TcvawvtfbjI795DUadBv3eENoGRucITNTZPi/e+bO96PME8c8PH8HWpcWwKOWK3qA82OMdWxalPXbXVjtgNeqwu6U/KWhrqrDHP18d2tTc7YYnEEn7HJAt0/bsyR6sqXKgWKmYcVr0GadHDniDWXeyrq9x4nevnMeJzhFsbyjBVctLM36uSu5r86ecu0167ZSmR7b0elBblFq5MRksj8wjuaiJnkm+kFxOkGnaTm2RPEHybL8XOo1AZUJNsMtmRL83mHEx42OHO7Go2DLlvUTzWXzB9qActGXaG/JPb2jCV29ehS+8vhEfv3opbttYg6saStFUbkehxYBwNIYLgz682jaIvx7pws/3nMWwP4z3b6ufzR9nzpXYjPFJVQza5pY69hlILYtRqZm2TCVX/Z4Qigtm79/x0iXF2HtmAO1DfvhC0aTF2qrRCZIe7D0zgFt/tAfhmIT73nVJ0gsc9WdOl227MOhHjzuIf3nkGH70nDz2WR27noshJHPl+hVlKLEZ0ZTheWKqaoos0GpEyq62QW8I7UP+eC9LldOM//3gZfj3W1bhyXu2422ba3NWvbHYVYA/fWQrbl5XhW8/3Yz3/2pfxj63P+2XL1a09nrHLWv//t9aYDPpcOfldQCAe65rwMpKO77wx8PocQfwxJFOXPft53HfzlbcvrEG/3rTqqSvv6axFH2eYNYJw8c7R/DIwQ68e2td2vPidU1liMSkpAmLD+09j/YhP/7xRjkLc9vGGnzwyiVYXeWI917mg56RAP7t0WP4yl+OIjJmUXggHMXeMwPxEkcA+OD2erx+VTkkCXj31roJfY8rla8vt5uSLuRoNQKNFbZ4SfVUBSNR9LiDqHSasXKc4SaHLgzh/b/ch9oiC35650YAmPS/R+ewH//w05cRjcXw9VtXo2skgB8+Ozp+fmdzL0KRWNrSSEC+WHBpfRH2KNN2Q5EYTvW4k7KNSUFbMHvQNraMfFApZU4MrhxmPYYzZNoGvWEUWjIHbetqnPCHo/AEI/jiG5oyfl6iCqcpZUcboAwiiUQnXSbc0uPB0mmursmrTJskSfj3x45j56lePHH39qxXjLL542sXsLLSgeVpUrr56vnmXnzo16/iuc9elfOpW7lyvNONmCRnfNKpLSrAk0e7UOU0o6bIkpR6LrEaEY5KGPaHUTjmasigN4Q9p/tx1/b6i640EpAzQ2a9FucH5PLIdJlMAKh0midd5jjeldeFKLGPjT1tcysxI5KtPDISkxCMxNL2XfZ7Q1iV4ZwzEy6rL8ZvX27DIwc7ASBpR5tqaakVGgH88sWzOHB+CNWFZvzyPZtTMtpqf5cvlPqCRR2YsbzMhq8/cQLBSBSBcAx6rUjbTD9frKx04JUvXpfz+9VrNagpNOPMmGEZ6ovlxN+ZSa+d9iL2TMwGLb5921qsrXbgK48cw4N721L2ew14Q3juZA9u21iNP+1vx8P72zMewye73HjiaBc+cc3SeIbLoNPg3tvX4U3f24Ub730BA94QGstt+P4d63HJotTM75UNJRBCniK5pjr9qojvPN0Mm1GHD25PfxFvXY0TxQUGPHO8G29eWwlfKILv/a0Fl9YX4YqlrqTP3bCoEL/YfRbBSHRCWYMfPNuCXncQX3nzynE/dyoeeuU8IjEJXSMBPHuyN6nMb++ZAQQjMWxvGP0ZhBC4923r8Ileb8bn27EuqStEUYEBt15SlZJ5WVlpx5/3d0zr+bZ7WG4fqXSaYDXp8Mf97egZCaB0zOvBU91u3Hn/Xjgtevz6fVtQYjOiwKDFa+eGcMv66gl9rwFvCO/82V4M+8N48AOXYnW1Ay+1DuC+na1468ZqLCouwFPHulFo0WPjoswXkLYudWHH8R6cH/Bh2B9GOCol/T8sLJBLJU92eeSgzZQuaJNv840pkd95qheSBFy1fDTYdloMGXva+r3BrK/51ZLz2zfVTjg2uGtbfXwieiKTXgtJAsJRCQbdxP69I9EYzvR5cU3j9PZl5k2mTZIkfPWx4/jprjNo7vagayQw/helEY7G8LnfH8K9O5pz/Ahn1kut/fCHo3htnKlRc0ltjF1Vlf4kt6jYggFvCEc6hrGoOPnFS7YF22pp5BvGqd9eqIQQqCkyo23AiwuD/pyWMl5sARuAeCmdTiNSpj7R7FpWaoXNqINOI2DMsCi6IMOVVlWfO4jiLGUvuaYOKfntXnn3Vbr+A5NeizpXAfaeGcCqSjv+8KHL0/6/VRfHBkKpVRRq0Pb1t6zBrRuqce+OU/jVi2exotJxUQwNmoo6V0G8Z1p1dMzwg9kghMC7ty7GprpC/PSFMwiPye48dqgDkZiEd1++GFctL8UjhzoQzVBl8r2/nUKBQYv3juklW1Zmw5f/biVikoQvvbEJj378irQBGwAUW41YV+PMOPr/0IUhPHWsG+/fVp/Sj67SagSubizFsyd6EI7G8Is9Z9HnCeKzNyxPuZi6odaJUDQ2ocFYkiThVy+exW9fbpv0JNWJiERj+O3Lbbh8STHK7EY8uLct6eMvnOqFQavBlsXJpZxGnXbCAZv6+X/79JW457qGlI+tqHDAHZzexEd1R1uVkmkDUgePSZKEjz+4H1qNBg+8bwvKHSZoNQJra5wTzrT5Q1Hcef9enB/w4ad3bsTqavn/zedf3widVuDfHj2OcDSGZ45345rGsqx9X1uVYH7P6b74xZOxlVgNZTac6plApm1Mxvq5k70oKjAkXYQotOgz9rQN+sJZyyNriy144H1b8M9vmliWDZD79m5eX5Vyu/pcNpm+tnMDPoSj0qRH/I+VF0GbJEn42hMn8LNdZ7BFSTsnjmSdjHP9PkRiEvac7s94ksxH6iSyw+PUMc+kruEAbv7B7pTyE9WR9mEUFxgyTgUbnSDpQ92Y2mr1hXS6sf+PH+lETZE5YzB4MagptODVc4MIRWMZyyNpYtTsmstqvCiD1nyi0QisX1SIgjQN6Co1A5duGEkgHIU7GJnVMtcSmxENZVacH/DDrNdmPN+9fVMt3rapBr95/6Up1QOqeKYtzeJYT1B+8eEw6/Ffb1mDt2+uhS8UzclS7YWqrrgAZ/q8SWVJRztGUOEwZX3BNlM+uH0J2of8ePxwZ9Ltf9rfjuVlNjRV2HDTukp0jwTxcmt/yte39Hjw2OFO3Hl5Xdpg6o4ttdj/z9fj/dvqxx3Ec83yUhy6MIwed+oF72891YxCix7vvaIu631c11SGkUAEfzvRgx8/dxrXNpamDRTV8t2JXGQ+3etF90gQoWgM+86lLmPOZKKlZzuO96BrJIB3X16H2zbW4LmTPfEACJCn125aXBjfBzYdzgwDkUaDrKm/flPH/Vc6zQnTopPv72jHCE50uXHP9ctQl3AxaUNtIU50ueHLsu9S9dSxLhxuH8a3b1uXtN+2zG7Cx69Zhh3Hu/HNJ09iJBDJWBqpWlZqRYnNiF0t/TjWMQKLQYvFY177qWP/R/zhtEGbUaeBRgC+MWtfXjk7gMvqi5Oymg6zIW3QFo1JGPSFUJSlPBIArljmimf2psOonNcnswJEjWmWzfegTZIkfPOpk/jJ861456WL8P07NgCYetDWqozRR3lVAAAgAElEQVTUHPaH4+Pp54PxgjZ3IIwt/7EDTx7tSvvxXHj0UAcOnB9KWT6pOtI+gpVVjowvvhKvNI/NtKnThMZOkBz2hbG7pQ9vWDW1yUcLRU2RJT7KuqZo/KZ4ykwN2tjPlh8+uL0eH716ScaPq0Gb2oeYqF8ZJz2bmTZALpEE5MxOpsD/A9vr8bVb12R9MRjvaUuTYfAoL1KsRh00GoH/uGUV/vvt69MOiSBZfUkBfKFo0sW/ox0jc1ZOek1jKZaVWvHj51vjQca5fi9eaxvCLRuqIITAdU1lKDBo8ecDyWPZJUnCfz15AiadNut0yYk+L17dKPf+qGsDVK+cHZDbL65cknEnlWrbMhcMWg0+9/tDGAlE8OnXLU/7eWV2E6qcZuw/PzTu49rdIvc8aQSwuyU1cE3nRNcI1nzlqQlljx546RwqHSZc01iK2zfVQALwu1fkZejdIwGc6HJj+7Lskz6na3m5DVqNmFZfmxq0lTtMsJv0qCu24Eh78v09vL8deq3AG8dMCd2wyIloTMKhCSxmf/F0P2wmXdqBLe+9og51xRb8ZGcrTHrNuL83IQS2LinGnpY+HGkfRlOFPeV8ubzcikA4hrP9vrRBmxACBQZdUqXFSCCMC4P+lKyd06LHcJpBJEO+ECQJs3bhxqRk2iYzjESNaZbM96DtVy+eww+ePY23b67Fv7x5JVxWAxxmfXzPylixmIQ/H2hPaTZVtSZkiXa1TGxvyVwbCYTRPiRPqDnSPpz2CtOr5wbRPRLEnhn8mZ462g0A2Hc29WpYMBJFc7cbq7I8OSYGaqmZNiVoG5Np23t2AOGolDR++GKkjv0HUne00eSopbi5njhIU7N1qQt3bc8ctKnBdboMgbpYu3iWexPVqXj1JVMfzQwgXuaYNmgbs2xWCIE3r61kH2YW6vOKWg3iD0WV4VizVxqZSKMRuGt7PY53jmCnsift4f0dEAJ4s7J82aTX4oZV5Xj8SGfSlfkH957Hk0e78Ylrl+Xk+F5ZaUeZ3YhnE0okW3rc+NKfjsBlNeJdE5hIWmDUxXcV/t3ayoxDxwBgfa0T+yeQadvV0ofaIgsuWVQYD+DGs7ulH+5gBF/764msGbfWXg92tfThji210Gk1qC604MqGEvzulTZEorH4jshtMxy0mfRaLCkpmPQe1UQdw364rIb4OWNlpQNHO0eDsGhMwl8OduCq5aUpWdn1NUrmcwJB7u7Tfbh0TAZLZdRp8f/+bgUA+Xc2kezk1qUu9HtDeLVtMO3Fk2XKMJJoTErb0wbIfaKJWcJmJYnRVJHce+Y06+EJRlLKkQd98sW9TBUPuaZm2iazYLulx4NKhynjFOWJmvOg7bFDnVhVZce/37wKGo2AEAJLS60ZM20vtPTh7ocOYMfx7rQfb+31xMtb1Kk2ueYLRXDT93flLIBSs2xXNZRg0BdOSu2r1OWhzd1Ty0COp88TxL5zA9BrBfadG0yZ8tjc5UEkJmUdCGAz6eNXOsZm2hxmPfRagT5PchPp/rZB6DQCa6rn5kk3X6hZSo2QyyNo6tRgjZm2+UHdi5XuvKf2wM52AL5lcTF0GoGGaU76yjQZDZDLI7UaAZN+zp+G543EXW0AcLxrBDEJczq45aZ1VSizG/GT509DkiQ8fKAdly4uTjqP37yuCu5AJJ4FO9E1gn955Ci2LXNlHAwyWUIIXL28FC+c6oM3GMG9O5rxhu/uQtdIAN94y+oJlwf+3dpKmPQa3DNmh+FYG2oL0TEcQNdw5vkDkWgML53ux9alLmxd6sKRjmEM+dIPkkikjrvfe2YgHnil85uX26DXivj6AwC4Y3MtukeCePZkL1441QuX1Zh2z1iurax0pFR3fe+ZU7j9Jy9OqNSzfSiQdMysrLLj/IA/Pi3xxdP96HEHcfO61B6rwgID6l0FeO1c9szn+QEfzg/4cXmWVQ3XNJbh/71pBe6+Nvu/v0rta5MkpJ3+nVgOaMsyjCrxHHlCeU28vDw10wYgJdvW71ErMmbnOV/NtE1mwfapHve0s2zAHAdtkiShuceNNdXOpJTq0hJrxkyb+p/5WIadX6d7vah3FeDyJS7sPTMwqZrTiXqptR8HLwzj6QyBIyCXSHSkeRGSjnqAvuUSefJPuv0cryjZr1M9mXedTcffjvcgJgHvvLQOw/4wmsd8H3WU8KpxrmjWFlmgEUhaYgrITyguqzFlEMn+tiE0Vdgv+sZ7dex/pdMcX9ZLU1PC8sh5pcxugk4j0J6miV+9yDPb2afCAgMe/uhWvG/b+EuRs8k28l/dW3Qxl4VPVqXTDINWE58gqWY25jJoM+g0eN8Vi7HndD8eeOkczvR5ccuY4QWXLymGy2rEnw+0wxeK4KO/eQ12sx7fuX1dTvtur24shScYwdXffA737jiFG1eV45lPXzmpiXW3bqjCK1+8Lu3U1EQbFo2f3TnUPgx3MIKtS4uxdakLkiQHH+M53D6M7Q0lqHKa8V9Pnkwb9PhDUfzfvvO4cVVF0p61axpLUWY34oGXzuGFU33Ytsw1K73NKyrs6BoJxKsDfv3SOXzr6Wa8fGZgQgNKOof8qHQkBG3Kay012/bwgXbYjDpc25R+v9i6WicOnB/MGiCqiYytY6aBjvXeKxZPeGJvpdMcv5iSbhiQzaSPX5jLmGnTJ2faTnSNwGbSodKR3E/sUDKMYwP/0Uzb7Awem2ymLRaTcLrHO+0hJMAcB229niCGfOGUxrwlpQXo84TSXpFRGzNPZKgdbu31oL7EiiuWuhCMxGZkl4h65ed4lvrlDz3wGj7x4P4J3d9J5QC9urEUOo1IqUsORWI4cH4IFoMWfZ4QBrzjX6lK55WzA/jKX46m3ZX21DF5VP+7lT0xr5xJLpE80j4Mm0k3br9VU4UNS0utMKSZFOeyGpPKI+Ua7CGsZ+N9/PfK0sjpc1mNuGldJa5pHH95Js09rUag3GFKe5ErfgV1DkpdV1U5pl3KkrU8Mhid9v1fbLQagdpiS3yC5LGOYTjMoy8K58rbN9fCZtLhK48cg0GnwY2rk/uFdFoN3rSmAs+c6MFn/+8QWvu8+O7t63J+MeKKpS5YjTrotRr8/N2b8N9vXz/p7yGEGLf3DZCDFINOk3UYyW7ltdLlS1xYV+NEgUGL3eNUQHmDEZzu9WBDrRP3XN+Aw+3DaXv5HznYgZFABO/YUpt0u06rwe0ba/B8cy8GvKGkUf8zSb1wcKxzBE8e7cKX/3wknnkar/dPkiR0DPmTM23qMJL2EQTCUTxxpAs3rirPeIF7Q20h+jwhnB/IHCDuOd0Pl9U47WEYY21f5oJRp8GysvT326DcnnntizY509bpRlO5PeWCljoNeuwwktHe5/zMtHUM++EPR7FsmpUbwBwHbaeUUj91AZ9KjUbTZdvUK2vHu1IDpgFvCIO+MJaUFGBLfRG0GjHhGurJUO/zeKc77VUNXyiCk10jeLVtMH7VJZsTnW40lttg0muxrMyWMozkaMcwAuEYblon18g3d08+2xaNSfjinw7jF3vO4vEjyZOufKEIXjjVh+tXlKGmyIwyuxF7zyafiI90jGBVZeYhJKovvnEFfvP+S9N+rMSWnGk71eOGNxTFuhoGbTaTHuV2U85PphcjjUbgu29Lv8+I8lOl05yxPNJi0OZk4tdcyJppC6afpkbZLXYVxMsj1SEkc52ttJn0eMelixCNSbi+qQz2NEHPzeurEIrE8NjhTnz8mmW4fJxsx1QUGHV4+lPbseNTV8YHk8wUg04z7pLtXS19WFlpR1GBAXqtBpsXF2HPOMNIjnWOQFL2wd6yvgpLS6345lPNSdPA+zxB/GzXGTSUWZMWXatu21QD9ZAYL6uUK2r/329fbsMnHtyPNdVOPPTBS2HSa7B/nOTBiD8CbyiKSudoZsllNaLcbsLRjmE8c7wHnmAk7fh5lbqHLNO/hyTJU9UvX1Kc8/8vn3rdcvzhw5dnDCgblPLUTOc7s0EHrxK0SZKEk13utLvU1PLIsUHbgGd2M22mSWbaTintXvM+06YGH2Oj86Ul8j/W2L62kUAY5/p9cJj1OD/ghzuQ/A+nTo6sLymAzaTH2mrHhKcVTVT3SADN3R4sKrZg2B9GZ5p67uOdcp29JKVOchpLkiSc7B49QFdX2VOGkexTAqg7NsvLQk9NIWh79FAHmrs9sBp1uHfHqaQT4M7mPgQj8uZ7IQQ21RXhlTMD8ccQjsZwvHNkQiP5rUZdxrI0l9WQlGk70CZffVLHB1/s/veDl+FTGaZ1ES1k1U5z2vLIfk9wTrJsuWLO2tOWftksZbfYVYCz/T6EIjGc6HLnzSLy92ytw/IyG+5UqlXGWlvtwMpKO7YuLcYnrpm5CaEVDnNOxttPxIZaJ450jKR98eoLRbC/bShpMffWpS609mVvHTmsVBqtrnJAqxH49PUNaOnx4OH97YjFJDzw0jlc883n0NrnwT3XNaQNQKoLLbhhRTk21DqTSidnktNiQJXTjL8e6UKl04z7370JdpMea6qcODBOpi1xR1uiVVV2HOkYwZ/2t6PMbkwa0T/W8nKbvGQ7Q9DW0uNBrzuIrUsz38dUOcz6rOWUy5XEjC3D+a7AoIVfKY9sH/LDHYygsSJN0GZWyiPH9LQN+EKwGnUTWvSeC0b95DJtpxdO0OaB06KP96CoqgrNMOg0ON2bvC9MbfLMlHFqVT5/iVKLvXWpC4cuDKUdETpVu5R0//uVEb3p1gqoJx2rUZdx2aWqYzgAdyCCRqXhcnWVI2UYyb5zA1hUbMGqKjtsRt2kh5GEozF85+lmNFXY8Z9/vxotPR48emh0/PBTx7rgMOuxuU6+YrV5cRG6RgLxOuzTvR6EIrEJ1zhnUmIzot8bipdn7m8bgtMij7Ylefmjg8ug6SJUVWhG10ggZSpYnyc0r6cpmnTZMm0sj5yKuuIChCIx7GrpRSgSm9Wl2tmU2kx48p7taTM/gFx2+IcPX44H3rdl3J1r88WG2kKEIrG0r4NeOSvvHd06JmgDkLUC6kj7MEptRpQq+xFvXFWO1VUOfPvpZtzyoz340sNHsLLSgb/evQ2vHzP6PtF3374uY9XPTFlf64TLasQv37M5PpRtfa0TR9vTB7aqxB1tiVZUOtDa68HzzT1489rKtBMfVeqS7f1t6QNE9Xd++ZLZyTwmunyJCxtqnVhRkf7/qtmghVdZgXJCmVfRWJ56McYRz7QltwgNeEOzuqdRPa9PNNPW0uNBcYEhJ49xjssj3WgotaVcKdFqBOpdBSmZNrU08tYN8sCOscNITvd5YFDGvgLyCSImIe1Sy6na1dKH4gIDblLS1On62o50jMBlNeCNqyuws1l+YsnkpFLmqU43UgMjdRiJJEnYd3YQGxcVQQiBZWVWnJxkpu2Pr13A2X4fPn19A964ugKN5Tbcu+MUItEYItEYnjneg2sbS+NPJJuU4G2v0tem7gqZ7pOjy2qML0EEgAPnh7CuxjnnpS1ENLeqnGbEJLmSIVGfJzhrfQozQaNMh/SnWXrrCbA8cirUoQePHpLL/PMl0zYRJr12QT3fxZdspwkUdrf0waDVxF9PAHLGxWU1ZA3aDrcPY3XCBWIhBD5zw3K0D/nRPujHvbevw28/sAVLx+kPMuq0s5ZxVH3jLWuw41PbUZtwIXpdjROhaAzHMwzPA+SeJwCocCZnBVdV2hGTgHBUyloaqVpf68TxzpG0PbR7TvejpsictE93tpQ7TPjjR7ai3JE+61lg0MUvbJ1QXhOnK4+0m3TQakRKImbAG5q1cf/A5DNtp3o8OZkcCcxh0CZJEpq73RkbF9ON/T+qXIFZU+2A3aRLGUbS2uvFomJL/GrE+lonTHpNzvraJEnCrpY+XL7UBbtJj9oiS9reuiPtw1hV5cC1TaVwByNp956p1P/Ias1vU4UdWo2I97Wd6fOi3xvCpjr55NhQZsOp7vS9dOkEI1H89zMtWFvjxLVNpdBoBD55XQPO9Hnx8IEO7D07gGF/GK9bOTpdanmZDXaTDvvOqUHbsLzp3jW9nUXqFfM+TwjugDyhUt0vQkQXL/UK89gSyT5PCCW2+VseCQCWhBckiTzBCIO2KVCfh54+2g2TXjPulEOaOeUOEyodprQlebtO9eGSRYVJgZNGI3DZEhd2n+7POA/gdK8nparnyoYS/O6uS/HMp6/Ezeur8jbwtRh0qTvUlMA2W19b+5AfBq0GrjEXqFYqv4dlpda04/TH2lBbiIgy4C1RNCbhpdZ+bJ2DLNtEWAxaeIPyha0TXW7UFJkzLuJ2mPWpPW3eEIrnINM2ken0kiShpceTk9JIYA6Dth53ECOBSMoQEtWSEivOD/qSfimJTceNFfb4qHzV6V5P0jJUo06LzYuLsXucEbMtPW7cv+vMuIHQyW43et1BbFNS/Csq7ClXTwLhKE71eLCq0oGtS10w6DTYcTxzieTJLjeqnOZ447JJr8WyUisOK9kttZ9to3K1qqHMhkFfOGXfWSYP7T2P9iE/Pvu65fET3Q0ry7Cqyo7/fuYUHj/cCaNOg+0No8snNRqBjXVF8Uzb0Y5hrFCCyelQe936PEEcvjAMSZLH1BLRxa2qMHVXWywmYcA7vzNtgDrOOsPIf/a0TVqZ3QizXiv3vZRP/3mJpmf9osKUJdv9niCOdY6k7Z+6Ymkxet3B+HCGRMc65HkAq9O0YmypL56X7QPlDhPK7aasfW0dQwFUOE0pqwkqHSZcVl+MD2yrn1CgminzeaR9GCOBCC7Lsp9tLlkMOgQjMURjEk50udOWRqqcZn1qT9tsl0fGB5GMn2nr84Qw7E+dkj9Vcxa0ZRpColpaaoUkjS7RDISjaEm4AtNUbsPJLne8PyocjaGt3xfvZ1NtXVKMlh5P1gWQP9t1Bv/66DE8eTTz3jVgtJ/timVy0NZUYcfZfm/SfoljnSOIKkuoC4w6XFZfjGdOdGcMCE92uVMWP66ucsSHkbxydgCFFj2WKMGoGuROZBiJPxTF959twZbFRUknTyEEPnV9A9oGfPjty23YtsyVMp1tU10RTvd60esO4mjHyLT72YDRTFuvOxgfgbuumkEb0cWuKk2mbcgfRkya/cXauWY2aFPKlWIxCd4Qe9qmQgiBuvheqPlTGrlQpVuyvUe5UJ5ucqPaU5WuAkqtMFpdnR99irmyLkuvGSD3tCXuaFMJIfDgXZcmLQ/PpqjAgMWugpTMp/rvMRf9bBOhTtkd9IXQ2uvJugzdYdEn9bRJkjTrQZsxPvJ//Eybult53mfamjOM+1epP6BaInmyy41oTIqfpBsr7PAEI/FhGecHfIjEpJRSCfWksSfLbpCD5+UTxVcfO5b1H2FXSx/qSwripTxNFTZIEpIyfkfHnHSuayrFuX4fWvu8KfcXisRwuteTUru7utqBAW8IHcMB7Ds3iEuUfjZgdN/FRMb+37ezFb3uID6dkGVTXb28FOtqnIhJwPUrUhdvbl4sX7H5v1fPwxeK5uTJMTHTtr9tCPUlBfHGUiK6eJn0WrishnhvB4D4epDieTyIBJBfkIwtj/QqF/oYtE3NYpfcl5MvQ0guZhuUapnfvXIezzf34vnmXvzlYAdsJl3ajFlNkQWLii0Zg7YSmxFl9tmZ+Dhb1tc60Tbgy7gCauyOtml9rxonXjrdjx3HRpMFe073oaHMmnGy91yzGOWg7dCFIcSk9ENIVE6zPqmnzReKIhiJzWrQptEIGLSaCWXacjk5EpjDoO1UtxtFBYaMk8EWuwogxOiutiPKUm31JN2k1PeqPWXq5MjE8khALmEstOixK0Nfmz8UxcluNzYvLsKFQT/u29ma9vOCkShebh2Il0YmPobEyUmH24dRVGCIb3JXd6X8LU2J5OleDyIxKSVoU7Naz53swZk+b7yfDZADH4dZj+Y0pQWJfvXiWXxnRzPetKYi7TQrIQS+9MYmrK5y4HUrylM+vqrKAYNOg1+/eC7pMU2H3aSDQatBrzuIA+cH2c9GRHGVTnP8IhyQGLTN70ybKU15pEfp32B55NTUFTPTli9WVNphNerwnR3NuPP+vbjz/r14+lg3ti8ryTglc+tSF/ac7k+ZAnhkzBCShULdRZuuRDISjaF7JIAqZ24C1Q9euQQldiPe/6t9+IefvowD54fwytmBvM2yAaOZttfOyb+fdOP+VU6LIT7MDpBLIwGgyDK7zxNGnWaCmTZ51VZ5ji5EzNkzRnO3O2uNp0mvRU2hJZ5pO9oxArtJh2ql96GhzAoh5PGgN6wsjwd3S1zJ96nRCGxZXIyXW9MPAznaMYxoTMIHttXDZTXgh8+14NZLqlP2Zbx2bgj+cDQp3V9daIbNpEuaIHm4XS4lVDNb1YUWNJbb8MyJbnxge33SfZ5UMnRNYxpM1f6xX+45C2C0nw2Qg62GMmvW8shf7jmLL//lKK5fUYZv37Yu4+dtrCvCIx+/Iu3HjDot1tU4sffMAAw6TU6uEggh4LIasP/8EPo8IfazEVFcldOcNBlX7dsduxJmvrEYtPEXFiq16Z6Ztqm5bkWZ3PuS5cUdzQ6jTosnPrkN3SPJWaR00/9U77psER7c24YfP9+Kz7++EYA8hKSlx4MbV2Ue4z9fra6Wd84dOD+Ea5uSK5u63UHEpNRx/1O1vNyGJz+5HQ/ubcN3nm7GzT/YDWD2loxPhdqe81rbIIw6TfyiTDpjB5HEg7ZZzLQBgFGvndD0yBZlcmSuhufMYabNk7E0UpU4QVIeQjIaDFkMOtQVF8QDptZeL1xWQ9pyuy31RWgf8uPCoC/lY+qVj7XVDvzTG5ogScB/PH485fN2tfRCqxG4dElyb1hTuT3+GALhKE51u7FqzNW/axpL8crZwZQxpSe63NBrRcpURnUYSXO3B0adJmWp9bIyG5q7PWn75H6++wy+/JejeN2KMvzgjg0w6Kb+T6zubWsqt0Gfo70yJTZjfJrm+hoGbUQkq3Ka0THkj5/X+hdIeWS6QSTuAIO26dhQW4j7371p1pbpUnbVhRZcsqgw6S3bsd1YbsdNayvxiz1n0KOs+TjemXkIyXxnMeiwvMyWtq8t04626dBrNXjXZXV47rNX44Pb67Gh1pm3Q0iA0UzbgfNDaCizZR0u5LTo4Q5EEFF2esaDtlmuyDDpNRPa09bS48nZEBJgjoK2cDQGdzAS78/KZGmpFWf6vAhFYjjROZJSCtFUYYvvdGjt86Delf7+1PJAdRpiokMXhlHhMKHUbkJ1oQUfvmoJHjvUiRfHTJzcdaoP62qc8SmPqhWV8hTLWEzCyS43IjEp5aRzbVMpojEJO5t7k24/0TWCJSXWtAGRWo64ttqZ8sTUUGrFsD+MHnfyla1fv3QO//LIMdywsgzfn2bABgCblN/byhyeRF1WI2KSfMBnazYlootLVaEZgXAs/iTc7wlBqxFwzsOJcYnSDSJheSRd7D55XQMiUQnff7YFAHD4gjIPYAEGbYDc13bw/FB8eJ5qJoI2lcOsxxfe0IQ/fmRrXl8gUjNtvlB03NeF6vPBiHLhq38OyyOD42Ta3AH5dfrYAYnTMSdBm5pSHG854pKSAgQjMexs7kUwEsPKMRmnxnI7zg344A1GcLrXm9LPlvh5dpMubdB28MIQ1iZMMPzQlUtQ5TTjn/50GF9/4gS+/sQJ/Odfj+NQ+zCuSJNebqqwwReKom3AF598NLb/a11NIQotejx1LHk6ZbrJkSr1xLWxLrXvS81QJg4j6XUH8e+PHcOVDSU5CdgA4JJFhagpMuOqhHUA06U2wq6pcmasdyeii098V5vyIqbPE0RRgSFlDPZ8Y9anDiLxMNNGF7k6VwFu21SDB/e24fyAD4fbR+CyGlFmn9+Z9UzW1TjhDkbirTyq9njQtrCGr0yGJWGXX7ayWgDxJdpqP+TgnGXatONm2tTqulzukJuTV83qDzqRTBsAPHygHUDqpKjGcnl6494zAxjwhjIGbVqNwKaEvWOqQW8I5/p9WJtQpmfSa/HVW1ZhwBvCz144g5+9cAY/33UWVqMON65KHdgRH0bSOYIj7cNwWvTxvrvE7//GNRV45GAHPvHgfgx4Qxj2hdE5HEBjhoWJmxcXQQgk7U9TqYu41QmcAPCj504jHJXwlTevzFkpo9WowwufuwavW5n6c0+VOniG/WxElGjs2P8+TyjjoKr5xJIt08agjS5in7hmGYQQuHfHKWUIiT1vF2dPV3zJ9phhJB1Dfjgt+pS1SxeTxKBt7IyHsdRdfequtn5vCHqtgG2Wz6XyIJLsmTZ1UEnigvnpmpOjJBCOoarAMG6vgppS3HG8G0adBvVjer/Uf9zHDncCQMbySEDua3vmRA96RgIoVaa4HFS2xq+tSQ4Gr15eioNfft2EfpaGMhs0Qq7HPtIxjFUJfXeJvvx3K1FiNeH7z57CntN9eMsl8t6NTFcVmirs2PfF69L+jlxWI4oKDPFhJF3DATzw8jn8/fqqlP64fKNm2tjPRkSJqgtTM23zfUcbAJgNOvjDUcRiUjxryKCNSF48fedli/CzXWcAADesTF0/tFDUuwpgM+mwv20It20c3bvWMRRIu6PtYpIYsI5bHqmUQQ4rw0gGvSEUWgyzHuyb9Npxp0f6Q3JQZ9bnLmibo/LIaMal2omcFgNcVgMC4RiaKuwp5XTVhWZYjTo8dbQLALAkS7Pf5sVyE+bes6PZtkMXhiHE9GqoTXot6kusOHhhGCe73BlH4+u1Gtx93TL85WNXoNxhwo+fPw0g+wGaLaiVB5XIQdsPn2tBLCbhE9cum/LPMVvW1zpR7yrAlvr8bYolotnnMOtRYNDGg7Z+bzCnZSVzRX3CDiSU0qjlkQUM2ugi9+GrlsKs1yIm5Wa1UL7SaATW1ThTxv7nckfbfKVm2lxW47jJHGc806b0Ps/yYm2VUTf+njafso8zl5m2OSqPjI07OVKlZtvS7WMRQqCx3LwB8NwAABXtSURBVIaRQAR6rUBNYeYDf1WlHRaDNqlE8uD5ISwpscJmml6j+4oKO3a39CEcTR1CMlZThR1/+shWfPaG5XjrJdVT3t3QUGbDqW4P2of8eGjvebx1Yw1qiixTuq/ZtKbaib995qo5+U9GRPlLCIFKp3m0PNK9cMojASSVSHpCERh1mpz0HhPNZ0UFBty1fQl0SlCzkK2vceJk1wieTlh83T7kz9mOtvlKvbDVNIEVHk5lQrw69n/AG5yT15MTyrTNQHnknDxjxCQJyyYYtKl9bWP72VRqiWRtkSXrYAudVoNLFhXG97VJkpQyhGSqmirsiCoTgSaStdNrNfjo1UvxX29dO+WUbkOZFe5gBF/602EAwMeuWTql+yEiyhdVhWa0D/nhC0XgD0fn/bh/YPQFSeLYf08gAhsnRxIBAD5+zVI88+kr460rC9VbN9ZgsasAH/jVPtzxPy/j5dZ+uAORiz7TptEIVDpMuGRR6uC9sWwmPYQYDdoGfeE5C9rGy7TFe9ryrTxSCHGjEOKkEKJFCPH5iXxNwwT3FowGbembE9XlmvUTGKl5aX0xTna7MegNoX3ILy94rpl+Ol69OmA36VBTNDv/+dSg99mTvXjb5pqUZeBERPNNlVMO2vqVxdoLo6dNKY9MuCrrCUZYGkmk0GgEFmVZqLxQ1BRZ8MQnt+Nfb1qJE10juP2+lwDMzLj/+eavd2/HR68eP/mg1QjYTfr4ZMZ+z9xk2uRBJNkzbeqFulwGbdN+1hBCaAH8AMD1AC4AeEUI8RdJko5l+7qJlkfesr4KUpaFi43lcjA3kT0I8X1tZwfimbG1OUjHr1Cyfauq0g8hmQnq78+g00zoQCciyndVhWYM+cJoG/ABwIIoj8yUaeMQEqKLj7r4+qZ1Vfjhsy149FBnTiq+5juHZeJtSk6LHoO+EMLRGEYCkbwvj7Tk2fTIzQBaJElqBQAhxEMAbgKQMWjTaUR818J4nBYD3nvF4owfX1Fhx9oaJ66cwC6xNdUOGHUa7D0zAJ1GwKDVxIO+6SixGbGqyo5rGkunfV8TVVRgwIoKO65tKkXZAi8pIKKLg1oxoE72XQhBW7ynbUymjUEb0cVLXXz9hTc0zfVDmXecZj2GfGEMKrva8nUQidrHbMqzoK0KwPmEv18AsCXbF5hymCo0G7T480e3TuhzjTot1tc6sffMACwGLZoq7TlpBBdC4NGPb5v2/UzWY5+4Yta/JxHRTFGDtkPnhwEAxQuoPDJpEEkwMuUhVEREFzOHxYAhfxgD3jkM2pSeNkmSMlbY+WegPDIXPW3pHq2U8klC3CWE2CeE2KePBnLwbadm8+JiHO0YxqELw1hXPb/HywohFuwiSiK6+FQVJmfaFsKUWTVo840J2qwcREJENGmFFj2GfaE5DdpMejl8ypZt84ej0GsF9FmGJE5WLu7pAoCahL9XA+gY+0mSJN0nSdJGSZI2VpfN3Y6uSxcXISbJv8xc9LMREVFulNpM0GkEOocDsJl0Oa3KmCsWvRycJZZHelkeSUQ0JU6zfu4zbTr5uSkYzh605fo5LBdB2ysAlgkhFgshDADeBuAvObjfGbG+thB6rZydWsPmTyKivKHVCJQ75LLBhdDPBgAmg/w061cWrQKAO8BMGxHRVDgsBgz7w+hzBwHMbaYtEMk8jMQfiua0NBLIQdAmSVIEwMcAPAngOID/lSTp6HTvd6aYDVqsqXbCZtSh3rXwR8wSEc0nal/bQhj3DwAWQ3KmLRSJIRiJwWpg0EZENFlOsx6SBJxTpgwXWvI305bLyZFAbgaRQJKkxwE8nov7mg2fvr4BncMBaDTsByMiyidVhWbgDFBcsDAybWNH/nuDcsaNmTYioslzKusBWnu9sJt0Oe0Zm6iJZtpyXR55UT5rXL7UNdcPgYiI0qhWM222hZFp02oEDDpNPNPmUYM29rQREU2aGrSd6fOieI7K6E0TzLSZc5xpm/3wlIiIKINKJWhbKJk2QN7Vpo5/VoM2GzNtREST5jDLF/QuDPpQOIml3LlknGCmLdflkQzaiIgob6hj/xdKTxsgl0iODdoKmGkjIpo0NdMWk4CiObq4p5Y9BsJZgrZwHg4iISIiypXl5TZYjTqsqLTP9UPJGbNBCx/LI4mIps1pHs2uFRXMUaZNp+xpy1YeyZ42IiJayEptJhz5lxvm+mHkVFKmLcDySCKiqXIkBW1znGnLVh45A9MjmWkjIiKaQel62lgeSUQ0eTqtJn7Rq3gOdrQBE8y0sTySiIhofjHpE8ojAyyPJCKaDrWvrXCOgraJZNp8oShMzLQRERHNHxaDFoGxmTYu1yYimhKnMkFyrjJt6sj/QIZMWzQmIRSJwaLP7XmeQRsREdEMMuu18IXlYM0TjMBq1EGjEXP8qIiI5qe5zrSpI/+DGTJt6lRJsyG3YRaDNiIiohlkNujgD8lXZD2BCAqMuS2ZISK6mKjDSOa6py1Tps2nVFawp42IiGgekQeRKJm2UIT9bERE06Bm2ormKGgTQsCg04ybacv1yH8GbURERDPIrNfCH45CkiR4AhFYTXOzW4iIaCFYXm5HTZE55yP1J8Ok02ScHulXgjZLjnuXebmPiIhoBpkNWsQkIBiJwROMwMZMGxHRlL3z0kV4x5ZaCDF3vcEmvTaeURsrXh7JnjYiIqL5Q+1r8Iei7GkjIsqBuQzYAHkYSTCSIdMWYnkkERHRvKOW8PjDUWV6JMsjiYjmM5Muc6YtMEPlkQzaiIiIZpBZCdp8ITlos5lYHklENJ9ly7RxeiQREdE8lFQeGWR5JBHRfJct06YOImHQRkRENI+oJTKDvhCiMYnlkURE81y2QSTxoC3H0y0ZtBEREc0gdYJYrzsIALCyPJKIaF4z6rINIpH3cjJoIyIimkfMejlI61GCNo78JyKa37Jm2kJyMMfySCIionlEvdqqZtoKGLQREc1rWTNt4SgMOg20mtyuJWDQRkRENIPUkf+9HqU8kkEbEdG8ZtRrEQhnLo/MdZYNYNBGREQ0o9QFqz0jAQDgyH8ionnOqNMgmGUQCYM2IiKieYaZNiKihcWk12Ypj4zFz/u5xKCNiIhoBum1Gui1gj1tREQLhEmvQSgaQzQmpXzMH4rEKyxyiUEbERHRDDPptXAH5DHQLI8kIprfjDo5KAulybb5w9Gcj/sHGLQRERHNOLVURqcRMOr41EtENJ+Z9PJ5PN3Yf38oyvJIIiKi+UhtSi8w6iBEbsdAExHR7FIzbYFIatDmC0VZHklERDQfmQ1ySSSHkBARzX8FRjko8wYjKR8LcHokERHR/GRWSmnYz0ZENP85LQYAwJAvnPIxf5jlkURERPOShZk2IqIFw2nWAwCG/alBG8sjiYiI5ilTQk8bERHNb06LHLSly7QFOD2SiIhoflJLZawsjyQimvecZqU8ckymLRyNIRyVYGGmjYiIaP5RgzYbM21ERPOezaSDEMCwL5R0u7oCgJk2IiKieUgtj2RPGxHR/KfRCNhN+pRMmz8kB23saSMiIpqH1Ewbe9qIiBYGp0Wf0tPmVzJtnB5JREQ0D6k7ezjyn4hoYXCa02Ta1PJIZtqIiIjmH7W/geWRREQLg8NiSBn571PLI5lpIyIimn/MLI8kIlpQnGZ96iASJWjj9EgiIqJ5iCP/iYgWFqclS3kkM21ERETzj1kvB2sc+U9EtDA4zXoM+8OIxaT4bWp5JHvaiIiI5qFNdYW4dUM1VlTa5/qhEBFRDtjNekgS4A5E4rcx00ZERDSPFVuN+NZta2ExMNNGRLQQOC0GAMCQf7SvLcDpkURERERERPnBadYDQNKutnh5ZL5l2oQQXxFCtAshDihvb8jVAyMiIiIiIspHTosctCWO/ferI/91uQ/aclGn8R1Jkr6Zg/shIiIiIiLKe2rQljhBMhCOwqTXQKMROf9+LI8kIiIiIiKaBIdZ7mlL3NXmC0VnpJ8NyE3Q9jEhxCEhxP1CiMJMnySEuEsIsU8Isa+3tzcH35aIiIiIiGj2OdL0tPnD0RkbODVu0CaE2CGEOJLm7SYAPwKwBMA6AJ0AvpXpfiRJuk+SpI2SJG0sKSnJ2Q9AREREREQ0mww6DSwGbVJ5pF8pj5wJ44aCkiRdN5E7EkL8D4BHp/2IiIiIiIiI8pzTrE/OtIWiMzI5Epj+9MiKhL/eAuDI9B4OERERERFR/nNYDBhO2NPmD0Vh0c9MeeR07/UbQoh1ACQAZwF8cNqPiIiIiIiIKM85zfrkkf/hKOxKr1uuTStokyTpnbl6IERERERERPOF06JHS48n/nd/KIoyu3FGvhdH/hMREREREU2S06JPGUQyZ9MjiYiIiIiIKJndrMewLwxJkgCo0yPzcBAJERERERHRxchpNiAUjcEfjgJQpkcyaCMiIiIiIsoPTsvogm1JkuAPR2E2zEx4xaCNiIiIiIhokpzm0aAtHJUQjUnsaSMiIiIiIsoXDiXTNuwPwx+SSyTZ00ZERERERJQnnGYDAGDYH4r3tbGnjYiIiIiIKE8k9rSpQZvFwKCNiIiIiIgoLzjUnjZ/GL5QBADLI4mIiIiIiPKGxaCFXisw5AsjoJZHMtNGRERERESUH4QQcJgNck9bKAaA5ZFERERERER5xWnRY8g3Wh7JQSRERERERER5xGnWyyP/wxz5T0RERERElHfUTFuA0yOJiIiIiIjyj13JtPlC3NNGRERERESUd5xmA4Z8Ccu1mWkjIiIiIiLKH06LHt5QFCP+CIQAjLqZCa8YtBEREREREU2B0yIv2O4a9sOs10IIMSPfh0EbERERERHRFDjMctDWORyYsX42gEEbERERERHRlDgtBgBA90hgxvrZAAZtREREREREU+Jkpo2IiIiIiCh/qeWRwUiMmTYiIiIiIqJ8ow4iAWZuRxvAoI2IiIiIiGhKbCY91IGRzLQRERERERHlGa1GwG6Ss23MtBEREREREeUhtUSSmTYiIiIiIqI8pE6QZKaNiIiIiIgoDzmUXW0M2oiIiIiIiPKQOvbfwvJIIiIiIiKi/KOWR5oYtBEREREREeWf+CASlkcSERERERHlH5ZHEhERERER5TGnMojExEwbERERERFR/uHIfyIiIiIiojy2qNgCIYBKp3nGvoduxu6ZiIiIiIhogVtWZsO+L16HYqtxxr4HM21ERERERETTMJMBG8CgjYiIiIiIKK8xaCMiIiIiIspjDNqIiIiIiIjyGIM2IiIiIiKiPMagjYiIiIiIKI8xaCMiIiIiIspjDNqIiIiIiIjyGIM2IiIiIiKiPMagjYiIiIiIKI8xaCMiIiIiIspjQpKk2f+mQrgBnJz1b7wwOQAMz/WDWCBcAPrm+kEsEDwuc4PHZO7wmMwdHpe5w+Myd3hc5g6Py9wZ77hcJElSyUTuSJebxzNpJyVJ2jhH33tBEULcJ0nSXXP9OBYCIcQ+Hpe5weMyN3hM5g6PydzhcZk7PC5zh8dl7vC4zJ1cHpcsj5z/HpnrB0CUBo9Lyjc8Jikf8bikfMTjMg8xaJvnJEnifyzKOzwuKd/wmKR8xOOS8hGPy/w0V0HbfXP0fYmy4XFJ+YbHJOUjHpeUj3hcUj7K2XE5J4NIiIiIiIiIaGJYHklERERERJTHchK0CSHuF0L0CCGOJNy2VgjxohDisBDiESGEXbm9TgjhF0IcUN5+nPA1twshDgkhjgohvpGLx0YXr8kcl8rH1igfO6p83KTczuOScmaS58t/SDhXHhBCxIQQ65SP8biknJnkcakXQvxSuf24EOILCV/D45JyYpLHpEEI8XPl9oNCiKsSvobHJOWMEKJGCPGscu47KoS4W7m9SAjxtBDilPJnYcLXfEEI0SKEOCmEuCHh9skdm5IkTfsNwHYAGwAcSbjtFQBXKu+/F8C/Ke/XJX5ewucXA2gDUKL8/ZcArs3F4+Pbxfk2yeNSB+AQgLXK34sBaHlc8i3Xb5M5Lsd83WoArcr7PC75ltO3SZ4v7wDwkPK+BcBZ5bmdxyXfcvY2yWPyowB+rrxfCuBVyIkJHpN8y+kbgAoAG5T3bQCaAawA8A0An1du/zyAryvvrwBwEIARwGIAp6f6+jInmTZJknYCGBhz83IAO5X3nwZw6zh3Uw+gWZKkXuXvOybwNUQZTfK4fB2AQ5IkHVS+tl+SpCh4XFKOTeN8+XYADyrv87iknJrkcSkBKBBC6ACYAYQAjIDHJeXQJI/JFQCeUb6uB8AQgI3gMUk5JklSpyRJrynvuwEcB1AF4CbIgReUP29W3r8J8kWuoCRJZwC0ANiMKRybM9nTdgTAm5X33wqgJuFji4UQ+4UQzwshtim3tQBoVMondZB/2MSvIcqFTMdlAwBJCPGkEOI1IcTnlNt5XNJsyHa+VN2O0aCNxyXNhkzH5e8BeAF0Qr5S/E1JkgbA45JmXqZj8iCAm4QQOiHEYgCXKB/jMUkzRghRB2A9gJcBlEmS1AnIgR3kjC8gB3TnE77sgnLbpI/NmQza3gvgo0KIVyGnD0PK7Z0AaiVJWg/gUwB+K4SwS5I0CODDAH4H4AXI5RaRGXx8dHHKdFzqAFwB4B+UP28RQlzL45JmSabjEgAghNgCwCdJ0hEA4HFJsyTTcbkZQBRAJeRyn08LIep5XNIsyHRM3g/5xfA+APcC2AMgwmOSZooQwgrgDwA+KUnSSLZPTXObNJVjUze1hzo+SZJOQC45gxCiAcAblduDAILK+68KIU5DznLsk+Rlfo8oX3MX5CcFopzJdFxCPtk/L0lSn/KxxyHX0j/D45JmWpbjUvU2jGbZ1K/hcUkzKstxeQeAJyRJCgPoEULshlyK1srjkmZSlteWEQD3qJ8nhNgD4JTyMR6TlFNCCD3kgO03kiT9Ubm5WwhRIUlSpxCiAkCPcvsFJGfQqgF0AJM/Nmcs0yaEKFX+1AD4EoAfK38vEUJolffrASwD0DrmawoBfATAT2fq8dHFKdNxCeBJAGuEEBYlTX0lgGNjvobHJc2ILMelettbATyU4Wt4XNKMyHJctgG4RsgKAFwK4MSYr+FxSTmX5bWlRTkWIYS4HnKWjc/hlHNCCAHgZwCOS5L07YQP/QXAncr7dwL4c8LtbxNCGJXS3WUA9ir3NaljMyeZNiHEgwCuAuASQlwA8GUAViHER5VP+SOAnyvvbwfwr0KICOSI8kNKLTwAfFcIsVZ5/18lSWrOxeOji9NkjktJkgaFEN+GPJlKAvC4JEmPKZ/H45JyZpLnS0A+Z16QJKl1zF3xuKScmeRx+QPl/SOQS39+LknSIeVjPC4pJyZ5TJYCeFIIEQPQDuCdCXfFY5JyaSvk4+uwEOKActs/Afga8P/buYMTAGEgiKKzRVqBRzuzH7sR4iFpQETYw3s15PI3MDmras88bG1JMsa4qurM/Ai4kxxr6C55+TZrzUwCAADQ0J9DJAAAAHwk2gAAABoTbQAAAI2JNgAAgMZEGwAAQGOiDQAAoDHRBgAA0JhoAwAAaOwBD8quIqr3rTMAAAAASUVORK5CYII=\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 convenient tables showing the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                               SARIMAX Results                                \n",
      "==============================================================================\n",
      "Dep. Variable:                   infl   No. Observations:                  203\n",
      "Model:               SARIMAX(1, 0, 0)   Log Likelihood                -472.714\n",
      "Date:                Mon, 24 Feb 2020   AIC                            951.427\n",
      "Time:                        22:49:06   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 did not 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": "iVBORw0KGgoAAAANSUhEUgAAA3IAAAEvCAYAAAAAWPPhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd1gU5/bA8e/sLl06SFOKBVCwgKhREzVqEk1s6Sa5pvebnpvcVJOb3n656ck1idH0xMRYEtO7GhvNCihSpHdY6rK78/uDElRAyi6wej7PwwPMzM68Lrjsmfe85yiqqiKEEEIIIYQQwnZo+nsAQgghhBBCCCG6RwI5IYQQQgghhLAxEsgJIYQQQgghhI2RQE4IIYQQQgghbIwEckIIIYQQQghhYySQE0IIIYQQQggbo+vvAXTGx8dHDQ0N7e9hCCGEEEIIIUS/iI+PL1FV1ffo7QM6kAsNDWXnzp39PQwhhBBCCCGE6BeKomS1t11SK4UQQgghhBDCxkggJ4QQQgghhBA2RgI5IYQQQgghhLAxEsgJIYQQQgghhI2RQE4IIYQQQgghbIwEckIIIYQQQghhYySQE0IIIYQQQggbI4GcEEIIIYQQQtgYCeSEEEIIIYQQwsZIICeEEEIIIYQQNkbX1QMVRVkBzAeKVFWNbt7mBXwGhAKZwEWqqpa389i5wMuAFnhHVdVnej1yIYQQQgghxElPVdXWzy0fbb9v77ijt3Vlf3vfA5jN5i5ta3uO9s7TyTiU9o7rciAHrAReA95vs+0+4GdVVZ9RFOW+5u//3fZBiqJogdeBM4AcYIeiKOtVVd3XjWsLIYQQQgghbJyqqpjN5tbPR39tNpsxGo1HfAaO2N8S4LR83941FKXd2KfD/cd7zPGO6+q2ruxrq6GhATrIouxyIKeq6h+KooQetXkRMLP561XAbxwVyAGTgIOqqh4CUBTl0+bHSSAnhBBCCCHECcBkMmEymTAajRiNRhoaGlq/NplMrYFZi84CIo2mKW7RaDStxyiKgqIo6HS6I77vakBkqzqbuevOjFx7/FRVzW++SL6iKIPbOSYIONzm+xxgci+vK4QQQgghhOgjqqq2BmomkwmDwYDBYKC+vp7GxkZMJhOKorQGaFqtFo1G0/qh0+lwcHA44QOvvtTbQK4r2vtpdRhaKopyPXA9QHBwsLXGJIQQQgghhGhHY2MjDQ0NNDQ0UF9fj8FgoLGxEfh7Jq0lWNPpdDg5OUmA1g96G8gVKooS0DwbFwAUtXNMDjC0zfdDgLyOTqiq6nJgOUBcXFzHc4lCCCGEEEIIizAYDNTV1VFeXk5DQ8MRs2ots2liYOltILceuAJ4pvnzunaO2QGMVBQlDMgFlgCX9vK6QgghhBBCiB5SVZWGhgZqamqorKyksbERjUaDvb09rq6u/T080QXdaT/wCU2FTXwURckBHqEpgPtcUZRrgGzgwuZjA2lqM3C2qqpGRVFuAb6nqf3AClVV91r2nyGEEEIIIYTojNlspr6+nurqaqqqqjCbzWg0GhwcHHB0dOzv4Ylu6k7Vyks62DW7nWPzgLPbfL8R2Njt0QkhhBBCCCF6zGQyUV9fT1VVFXq9HlVV0el0ODo6tlaHFLapL4qdCCGEEEIIIfqI2WxunXWrqalpLdvv4uIiRUlOIBLICSGEEEIIcYJobGwkNzeXhoYGHBwcZL3bCUwCOSGEEEIIIU4AtbW15ObmotVqJYA7CUggJ4QQQgghhA1TVZWKigoKCwtxdnZGp5O3+CcD+SkLIYQQQghho0wmE0VFRVRWVjJo0CApYHISkZ+0EEII0QO7cyp5dP1eVFXt76EIIU5SBoOBw4cPU11djZubmwRxJxn5aQshhBA98GVCDiu3ZFJQVd/fQxFCnIRqamrIzMzEbDbj4uLS38MR/UACOSGEEKIHUgv0AGSU1PTzSIQQJxNVVSkrKyM7OxtHR0dp5H0Sk0BOCCGE6IG0wqZALrOktp9HIoQ4WZhMJvLy8iguLsbNzU2Kmpzk5KcvhBBCdFOxvoHSGgMAGSXV/TwaIcTJoKGhgdzcXMxms7QWEIAEckIIIUS3tczGKQpkyIycEMLK9Ho9eXl52Nvb4+zs3N/DEQOEBHJCCCFEN6U0r4+LDfYks1TWyAkhrENVVUpLSykpKcHFxQWtVtvfQxIDiKyRE0IIIboprUCPt4s9cSGeZJfWYjJLCwIhhGUZjUby8vIoLS3F1dVVgjhxDAnkhBBCiG5KKdQT7udKqI8LBpOZvIq6/h6SEL12/5pd/LC3oL+HIWhaD5ednU19fT2urq4oitLfQxIDkARyQgghRDeYzSoHCvVE+LsS6t3Uu0nSK4WtK9LX88n2w3y/t7C/hyKA8vJyVFXFycmpv4ciBjAJ5IQQQohuyCmvo9ZgIsLflTCf5kBOeskJG5eYXQFAQZXMLg8Eqqqi0cjbdNE5+Q0RQgghuiG1uWJlhL8rfm4OONlppXKlsHkJ2eUA5FfW9/NIhBBdJYGcEEII0Q2pBVUAhPs1rVsJ8XaW1Eph8xKzmmfkKutRVSneI4QtkEBOCCGE6IbUwmqGeDoxyKGpg88wXxdJrRQ2rdFkZlduBQ46DbUGE/oGY38PSQjRBRLICSGEEN2QVqAnws+19ftQbxeyy2oxmsz9OCohei4lX099o5mZEb5A06ycEGLgk0BOCCGE6CKD0Ux6cTUR/m0COR8XjGaVXGlBIGxUy/q4s8cEALJOTghbIYGcEEII0UUZJTUYzeoRgVxL5cpDkl4pbFRCdjl+bg7EBnsCUCiBnBA2QQI5IYQQootSmgudHDEj5y0tCIRtS8guJzbYEz83R0Bm5ISwFb0O5BRFiVAUJanNR5WiKHccdcxMRVEq2xyzrLfXFUIIIfpaWqEenUZhmM+g1m0+g+wZ5KCTQE7YpGJ9A4fL6ogJ9sBep8FnkL30khPCRuh6ewJVVVOB8QCKomiBXOCrdg79U1XV+b29nhBCCNFfUgv0hPm4YK/7+z6ooiiE+jiTUSq95ITtSWxeH9eSVunv7ijFToSwEZZOrZwNpKuqmmXh8wohhBD9LrVQf0RaZYtQb2lBIGxTQnYFdlqF6CB3APzdnCS1UggbYelAbgnwSQf7piiKkqwoyreKokRZ+LpCCCGEVVU3GDlcVndE64EWYT4u5JTXYjBKCwJhWxKyyxkd6I6jnRYAf3cHCqokkBPCFlgskFMUxR5YCKxuZ3cCEKKq6jjgVWBtJ+e5XlGUnYqi7CwuLrbU8IQQQoheOVCoB+hwRs6swuFySa8UtsNoMrMrp4KYoR6t2wLcnaiobaS+0dSPIxNCdIUlZ+TmAQmqqhYevUNV1SpVVaubv94I2CmK4tPeSVRVXa6qapyqqnG+vr4WHJ4QQgjRc6kFnQRyPlK5UtielIKmRuCxIZ6t2/ybK1fKOjkhBj5LBnKX0EFapaIo/oqiKM1fT2q+bqkFry2EEEJYVWqhHic7LUM9nY/Z19JLLkMCOWFDEloLnbSdkZMWBELYil5XrQRQFMUZOAO4oc22GwFUVX0LuAC4SVEUI1AHLFFVVbXEtYUQQoi+kFqgJ9xvEBqNcsw+T2c73Bx1ZJZKICdsR0JWOYNdHQjycGrd5tccyEkLAiEGPosEcqqq1gLeR217q83XrwGvWeJaQgghRH9IK9QzK3Jwu/sURSHMdxCZJbJGTtiOxMMVxAR70Jw0BfydWikzckIMfJauWimEEEKccEqqGyipNhDeTsXKFmHezpJaKWxGSXUDWaW1rf3jWrg46HBz1FEogZwQA54EckIIIcRxpDUXOon0d+vwmFAfF/Iq66Tan7AJidkVAEcUOmnh7+4oM3JC2AAJ5IQQQojjSG1uPRDuP6jDY8J8XFBVyC6T9Eox8CVkl6PTKIxpbgTelr+7k/SSE8IGSCAnhBBCHEdqgR5PZzt8Bzl0eEyot1SuFLYjMbucqEC31kbgbQW4OUr7ASFsgARyQgghxHGkFuqJ8Hc9oijE0aSXnLAVRpOZ5MOVxAQfm1YJTamVxdUNNJrMfTwyIUR3SCAnhBBCdMJsVkkr0BPRSaETAHcnO7xc7KUFgRjwUgr01DWaiGnTP64tf3dHVBWK9A19PDIhRHdIICeEEEJ0IreijhqDiYhOCp20CJXKlcIGJLY2Au94Rg6goFJ6yQkxkEkgJ4QQQnQitbliZUQnhU5ahPq4SC85MeAlZlfg6+rAEE+ndvcHtAZyMiMnxEAmgZwQQgjRidaKlcdJrQQI83ahoKqeOoO0IBADV0J2OTFDPTpc8xng1hTg5cuMnBADmgRyQgghRCdSC/QEeTjh6mh33GNbC57IOjnRTfesTmbZuj1Wv05pdQOZpbXt9o9r4eakw9FOI5UrhRjgJJATQgghOpHWXLGyK8KkcqXogfTialbH5/Dp9sNU1jZa9VqtjcA7WB8HoCgKAe5O5EsvOSEGNAnkhBBCiA40msykF1d3Ka0S/p6Ry5AZOdEN72/JRFHAYDKzcU++Va+VeLjjRuBt+bs5UigzckIMaBLICSGEEB3IKKmh0aQS2cUZuUEOOnwGOciMnOiyqvpGvojP4dyYIIb5uLA2Mdeq10vIqmBUgBtO9sc2Am/L392RfAnkhBjQJJATQgghOtBSsbKrM3IAw3xcpAWB6LLVO3OoMZi4amoYi8YHsS2jjLwK6xQZMZrMJOdUENtB/7i2/N0dKayqx2xWrTIWIUTvSSAnhBBCdCC1QI9WozB8sEuXHxPq40yGtCAQXWAyq6zakklciCdjhrizaHwgAOuT86xyvdRCPbUGU6eFTloEuDtiNKuU1hisMhYhRO9JICeEEEJ0ILVQT5iPCw66ztPQ2gr1caGkugF9vXWLVgjb90tKEdlltVw1LQxo+t0ZP9TDaumVCV0odNLC362ll5ykVwoxUEkgJ4QQQnQgtUBPRDfSKqGplxxAVqnMyonOrdySQYC7I2dF+bVuWzw+kJQCfWtaryUlZpfjM8i+w0bgbfk3NwWXXnJCDFwSyAkhhBDtqDUYyS6r7XLrgRatlStlnZzoRGqBns0HS1k6JQSd9u+3Y/PHBaLVKKxNsvysXGJ2BTHBnh02Am+rJZArkBYEQgxYEsgJIYQQ7UgrrAa6V+gEINRbesmJ41u5JRMHnYZLJgYfsd1nkAOnjvBhfVKeRQuNlNUYyCip6VJaJYCPiwM6jSKplUIMYBLICSGEEO1Ia05t62rrgRZO9lr83Ryll5zoUEWtga8Sm1oOeLrYH7N/cUwguRV17Mwqt9g1E7ObztWVipUAGo2Cn5ujBHJCDGASyAkhhBDtSCnQ42inYaiXc7cfG+rjLDNyokOfbD9MfaOZK6eFtrv/zNH+ONlpLZpemZhdgVajMGZI543A25JeckIMbBLICSGEEO1IK9QT7ueKVnP89URHC/NxIVOKnYh2GE1mPvgrk6nDvYn0d2v3GBcHHWeM9mPj7nwMRrNFrpuQXc6oAFec7XVdfoy/u6OskRNiAJNATgghhGhHanMg1xOh3i6U1RiorJMWBOJIP+wrJK+yniunhnZ63OKYQCpqG/k9rbjX1zSZVZIPV3R5fVyLgObUSlWVpuBCDEQSyAkhhBBHKasxUKxv6Pb6uBYtlSslvVIcbeXmTIZ6OTF7lF+nx5020hcvF3uLpFemFuipMZi6Hcj5uztS12iiqs7Y6zEIISzPIoGcoiiZiqLsVhQlSVGUne3sVxRFeUVRlIOKouxSFCXWEtcVQgghrKGlh1dPZ+TCWgI5KXgi2tiTW8n2zDKumBJ63JRdO62Gc8YE8NO+wl43l0883FLopPuBHEB+lfSSE2IgsuSM3Omqqo5XVTWunX3zgJHNH9cDb1rwukIIIYRFpRZUAXS7h1yLYC9nFEV6yYkjvbc5E2d7LRfGDe3S8YtjAmkwmvl+b2GvrpuQVYG3iz1DvY7fCLytgJZeclLwRIgBqa9SKxcB76tNtgIeiqIE9NG1hRBCiG5JLazGw9mOwa4OPXq8o52WQHcnCeREq5LqBjYk53HBhCG4O9l16TGxwZ4M9XJiXS/TKxOzy7vcCLwtf/emwE8COSEGJksFcirwg6Io8YqiXN/O/iDgcJvvc5q3CSGEEANOakEV4X6u3X7j21aYj4uskROtPt6WjcFk5orjFDlpS1EUFo0LYvPBEor0PQumymsMHCqpITaka/3j2hrs6oCiIC0IhBigLBXITVNVNZamFMp/Kooy/aj97f0lbLcEkqIo1yuKslNRlJ3Fxb2v1CSEEEJ0h6qqpBVW97jQSYtQH2cySmqk4p/AYDTz4dYsZoT7Mtx3ULceuzgmELMKG5Lze3TtpMMVQPfXx0HTOj2fQQ4yIyfEAGWRQE5V1bzmz0XAV8Ckow7JAdomhA8B8jo413JVVeNUVY3z9fW1xPCEEEKILsutqKO6wdjjQictQr1dqKo3Ul4rLQhOdt/uyadI38BVHTQA78yIwa5EBbr1OL0yIbscrUZhbDcagbcVIL3khBiweh3IKYrioiiKa8vXwJnAnqMOWw9c3ly98hSgUlXVnt1aEkIIIaworbCpYmVvZ+RaKlfKOjmxYnMmw3xcmD6yZzeoF48PYldOJYeKq7v92ITsciL9u9cIvC3/5l5yQoiBxxIzcn7AJkVRkoHtwDeqqn6nKMqNiqLc2HzMRuAQcBB4G7jZAtcVQgghLC6lufXAyN7OyEkvOUFToZHkwxVcOS0UzXFaDnRk4fhAFAXWJrWbzNQhk1klKbv7jcDb8nd3JL9S2g8IMRD17PZMG6qqHgLGtbP9rTZfq8A/e3stIYQQwtrSCvQEujt2ubJgR4Z6OqNRpJectTWazPzv93ROjxxMVGDP0get6b3Nmbg66DgvdkiPz+Hn5sjU4d6sS8rlzjkju1yE50BRcyPwHhQ6aeHv7khVvZFag7HHs3pCCOvoq/YDQgghhE1IKdAT3su0SgB7nYYhns6SWmlFjSYzt3+ayAs/pHH9+/G9bpxtaYVV9Wzcnc9FE4cyyKF3QdCi8UFklda2Fi/pioSspmNjhvZ8Rk56yQkxcEkgJ4QQQjRrNJk5VFzT40bgRwv1cZEZOSsxmszc8VkSG3cXcOnkYPIr63hsw77+HtYRPtyahUlVuWJKaK/PNTfaH3udhnXdSK9MyC7Hy8WeEG/nHl/Xz00COSEGKgnkhBBCiGZZpTUYTGYierk+rkWYtzOZJbXSgsDCjCYzd32ezDe78nnw7FE8de4Ybp45gtXxOfywt6C/hwdAfaOJj7dlMzvSj+BeBFIt3BztmDNqMF/vysNoMnfpMQnZ5cQGe/SqH2JAc1Nw6SUnxMAjgZwQQgjRrKXQSW9bD7QI9XGhusFISbXBIucTTQU8/rU6mfXJedw3L5Lrpg8D4LbZI4kKdOP+Nbsp1jf08yhhQ3IepTUGru5By4GOLBofREm1gU0HS457bEWtgUPFNcT0otAJNFWtBKQFgRADkARyJyizWWVHZhkv/pDK4bLa/h6OEEJYnNls+VmutAI9GgVGDO5e0+aOtFaulPRKizCZVe5ZnczapDzuOSuCG2cMb91nr9Pw0sXj0TcYuX/Nrn6dBVVVlfc2ZxLh58qU4d4WO+/MCF/cHHVdSq9MbF5LFxPc80InAE72Wjyc7SS1UogBSAK5E4iqqiRkl/PYhn1MfeYXLnzrL1755SCXr9hOaXX/350UQghL+TWliElP/cyunK4XfuiKlAI9oT4uONppLXK+MO/mXnLFEsj1ltms8u8vd7EmMZe7zwjnn6ePOOaYkX6u3HtWBD/tL+LznYf7YZRNdmSWsy+/iiunhfYqrfFoDjot54wN4Pu9BdQajJ0em5hVjkaBcUN6F8hB06ycpFYKMfBIIGfjVFVld04lT2/cz6nP/sp5b2zhw61ZRAe58dLF43n/6knkVdRx9aqdx33RF0IIW1BrMPLQ2j2UVDfwxNf7LTrzklao73Uj8LaGeDqh0yhkyIxcr5jNKvet2cUX8TncMWckt84e2eGxV08LY8owbx7bsI/s0v7JSHlvcwYeznYsHh9k8XMvGh9ErcHEj/sKOz0uIbuCSH83XHpZLROaWhAUVEkvOSEGGgnkbJCqquzPr+L571M4/YXfWPDaJt7dlMFIv0G8cOE4djw0h3eumMjimCCmh/vyyiUx7M6p4NaPE7u8QFoIIQaql38+QG5FHRdOGML2zDJ+2l9kkfPWGoxkldVabH0cgE6rIdjLecA2BT9cVstdnydRMoCzNsxmlQfX7ubznTncNmsEd8wJ7/R4jUbhhYvGoVEU7l6dhMkKKbidySmv5fu9BSyZGIyTvWVmdtuaFOpFoLsjaxNzOzzGZFZJOlzRq/5xbQW4O1JQOXB/R4SwZaqqtn4AmEwmGhsbaWxsxGAwYDR2PBEjnR1tyMEiPRuS8/l6Vx7pxTVoFJg63IcbZwznrCh/PF3s233cWVH+/GdRNA+v3cPD6/bw1LljLJrqcbwxf7g1mzvmjMTDuf3xCSFEV6UW6Hn3zwwuihvCU+eOIT67nGe+3c/pEb7otL27N3mwqBpVxaIzctC0Tm4g9pIr1jew9N1tZJbWMinUiyWTgvt7SMdQVZWH1+3hk+2HueX0Edx5RudBXIsgDyceXRjF3auTefvPQ0espbO2D7ZmoSgKS6eEWOX8Go3CgvGBvPNnBqXVDXgPcjjmmINF1VQ3GIntZaGTFv5uTpRUN2AwmrHXyRzAyU5VVaqqqjAYDBgMBhoaGjAYDHh6euLn50dDQwNbtmxpDUaMRiNGo5GoqCgiIyOpqqpi9erVrdtbPk4//XRiY2PJz8/njTfeoLGxEbPZjMlkwmQycdlllzFx4kQOHjzIc889d8Q+s9nM7bffzsSJE0lISOCJJ55o3aeqKmazmSeffJLY2Fh+//13HnvsMcxmc+uHqqq89dZbREdHs379eh5//PHW7S2fv/zyS0aMGMGHH37I008/DdC6T1VVfvnlF4KCgnj99dd58cUXjwjQVFUlKSkJT09Pnn76aV577bVjntfMzEzs7Ox4+OGHWbVqVev2pUuXdvizkEBugMssqeHrXXl8vSuflAI9itJ0N+7KaWHMi/bHp50X8PYsPSWEgso6Xv81HT83x+Pe0bSE7RllXLtqB1X1RhqMJp4+b6zVrymEOHGZzSoPrd2Nq6OO++aNQqfVcN/cSK7/IJ5PdxzmH6f07o2zpStWtgj1duGv9FJUVe2zm2jHU1nXyOUrtlNY1cAgBx07s8oHXCCnqirL1u3lo23Z3DRzOHefGd6t5++82CB+3FfIiz+kMSPcl1EBblYcbZNag5FPtx/mrCg/gjycrHadxeOD+N/vh/hmdz6Xt9OjLiG7HKDXFStb+Ls3vdcorKpnqFfvWykI68vIyECv11NbW0ttbS11dXX4+PgwefJkAN544w3Kysqoq6trPWbSpElcc801AMyZM4fq6uojArXLLruMRx99lIaGBkaPHn3MNW+55Rbuv/9+ampq+Mc//nHM/n//+99ERkZSUVHBsmXLjtnv6elJbGwser2eL7/8Ep1Oh06nQ6PRoNFomDNnDgD19fWkpaWh1WrRaDRotVq0Wi0GQ1N1YEVR0Gg02NnZte5TFAV7+6YJBWdnZ4YNG4aiKK3HajQanJ2bfrf9/Pw47bTT0Gg0rccoioKLS9Oa52HDhrF48eLW/S3XbHn82LFjueqqq454LICjY1MF2KlTp6LT6Y54bMs4Wp77wYMHt+6PjIzkgw8+aPfnLIHcAFVZ18hDa/ewIbmpMtWEEE8eWTCas8cEtDbn7K5/nRlBQWUDL/10AH83R6v+0f5mVz53fp7EEE8n5oz245Pth7lgwlAmhFjmj4oQ4uTzRUIOOzLLee78sXg1ZyCcMdqPSaFevPRTGotjghjUi/VAaQV6HHQaQpoLlFhKmI8zdY0mCqsa8Hfv2eu3JdUZTFy3aicHi/S8c8VEPvgri4Ss8v4e1hFUVeXR9Xv5YGsWN0wfxr1nRXQ7CFYUhafOG8OZ//2DOz9LYt0t03DQWT7Vsa21iXlU1jVy1bQwq15nVIAbEX6urE3MbT+Qy2pqBB5qgf51AP7NveQKJJDrMyUlJeTl5VFXV0dFRQWVlZU4OjqyYMECAJ599llSUlKorKyksrKSiooKxowZw8qVKwG47LLLyMrKOuKcZ5xxRmsg9/bbb1NVVYWTkxPOzs44OTkRFvb37+2oUaOApuDDzs4Oe3v71sc6ODjwyCOPYG9vj4ODA/b29tjb2zNyZNPaVTc3N9atW4e9vT12dnatAZmnZ9N7wKCgIHbv3n3EvpagDCA8PJx9+/Z1+NxER0fz22+/dbg/JiaGL774osP9EydOZOLEiR3unzx5cuu/tT1Tp05l6tSpHe4/7bTTOO200zrcP2PGDGbMmNHh/lmzZjFr1qzW72trO17rK4HcALQjs4w7Pk2ioKqeW2eNYMmkYIvc2VMUhWfOH0NxdQMPrt2Dr6sDs0f5WWDER3p3UwZPfLOPCcGevHNFHHZaDVsOljYFprdM63X6kxCi92oNRrQaBTuNBo1mYMwSdaa8xsDTG/cTF+LJBROGtG5XFIUHzhnF4tc3s/z3dO46M6LH10gt1DPSbxBaCz8fLS0IMkpq+j2QazSZueXjBHZklfHKkhhmhPuyL6+Kn/YXUlZjaA2Q+5Oqqjz29T5W/ZXFtaeGcd+8yB7PZHq52PPcBWO4euVOXvwxjfvnjbLwaP+mqiort2QQHeRGXB/ctFwUE8hz36WSXVp7TMPxhOxyYob2rhF4WwHNv7fSgqDv3HzzzWzatOmIbeHh4a2BXHp6Ojk5Obi7uxMaGoq7u3tr8AXw+OOPYzQacXJywsXFBWdn59ZACiA+Pr41cGrPq6++2uE+RVG4/vrrO9yv0+mIi4vrcL9Wq8XLy6vD/aLrJJAbQIwmM6/8cpDXfjnAEE9nvrhxisXSIlrYaTW8eVksS5Zv5Z8fJ/DJdadY7Bpms8qTG/fz7qYM5kb589KS8a0lvB9ZMJqbPkpg1V9ZXHOqde9UnoiKquq5f81uli0YbfHZAnFyMTWXcP8iPqd1m06jYKfVYKdVsNdpmr9u+t5Oq2mzTcFep8Veq+Bsr3y4pcEAACAASURBVOPWWSMYaeE0xI48820K+nojT5wbfUzgOX6oB/PHBvD2nxlcdkpIj7MWUgv0nDbS1xLDPUKo99+95CzZU6y7zGaVe7/Yxc8pRTyxOJoF4wIBiAtt+hsQn1XOGaMtf3OvO1RV5clv9vPe5kyumhbKg+eM6nUwMivSj0smBbP8j0PMjvRjUph13kBuPlhKWmE1L1w4rk9SaBeOawrk1iXlHlHFs7K2kfTiGs6LHdLJo7un5f+UBHKdq6xr5FBxddNNGzdHpo7w6fG5br75Zi666CJ8fHxwd3fH3d0dD4+/i9csX76808fPnj270/2dBXHCdkggN0AcLqvljs+SiM8q57zYIP6zMApXRzurXMvFQceKKydy/ptbuGbVTr68aSphPr0LDuobTdz9eTLf7M7nyqmhPDx/9BF3tedG+zMzwpcXf0jlnDEB/X5X2tZ8vD2bn1OKqGs08dG1kwfMOhthW4wmM3evTmZdUh6XTm6a6W80mZs/VAxGMwaTmUZjm22t+800GlWq6hppNJnJLqtlZ2YZa2+ZxmBX6/5/3plZxmc7D3PD9GFE+re/zunesyL5fm8BL/6QxrMXdH89bnmNgSJ9AxH+lmkE3laghxP2Wk2/Vq5smeX6KjGXf50ZfsR6wjFB7thplX4P5FRV5elvU3hnUwZXTg1l2fzRFnute+icUWxJL+Guz5P49vbTrPL3deWWDHwG2bNgXIDFz92eIZ7OTAr1Ym1SLrfMGtH6XCUebl4fN9QyFSsB3Bx1ONtrpZcc0GA0kV1ay6GSGg4V15BRUt38uYbSGkPrcQvHBfYqkJs5cyb19fU4OHStFoI4OUkgNwCsT87jwTW7AXh5yXgWWaHvzNF8XR1YdfUkzn9zC5ev2Maam6bh69qzF4vK2kau+2An2zPKePDsUVx7Wtgxf3wVReE/C6M4879/8Pg3+3j90lhL/DNOCmazypcJObg66tiSXspXibkWvdMqTg5Gk5m7Pk9mfXIe95wV0W4z5e7Yk1vJhW/9xXWrdvLp9VOsUmYdmlIBH1q7h0B3R27rpHdYsLczl08J5b3NGVx9ahgR3aw8mVponUInAFqNQrC3c79Wrnz1l4Os3JLJNaeGHfOzd7TTEhXoTnxWWT+NrimIe/a7VJb/cYilp4TwyALLBXHQdAPzxYvGceFbf/H41/t47oJxFjt3rcHI8j8O8XNKEbfOGmn1dXhtLYoJ5MGv9rA3r4roIHegqX+cRoFxFgzkFEU5qXrJmc0qBVX1ZJTUcKi4uk3QVkNOeS1tO1r4DHJgmK8LZ4z2I8zHhWG+gwjzcSFY1hKKPiCBXAcOFOp5+89DzB7lx6zIwdhZYV1XdYORZev2sCYhl9hgD15eEtOni4jDfFxYceVELlm+latWbufT66d0u1BATnktV763g+zSWl65JIaFzak67QnxduGfp4/gxR/TuDiumOnhlk9hOhFtyyjjcFkdL140jg+2ZvHEN/s5PWJwh+0mhDia0WTmzs+T2ZCcx71zI7h5Zu+COIDoIHdeXjKeGz6M5+7VSbx2SaxV1tq9tzmDlAI9y5dOOG5j41tnjWD1zsM8/e1+Vl41qVvXSWsO5Dqa8eutUG8XMvupKfgHf2Xy4o9pnB87hAfPbj9VMS7Ek/e3ZvVLeXlVVXnhh1Te+j2dyyYH85+FUVbJOpgQ4sWNM4bzxm/pzBnlx5lR/r06n8mssnrnYV78MY0ifQNnj/Hn2tP6dunAOWMCeHT9XtYm5rYGconZ5URYqBF4W0295E78GbnkwxVc8d52KmobW7c52WkJ83Fh7BB3Fo8PbA3WwnxdcLNS9pQQXSGBXDvyK+u4fMV28ivr+XxnDoNdHbgwbggXxwUfs6C4p5IOV3D7p4kcLqvlttkjuW3WiH4pAjJ+qAevXxbDde/Hc/NHCbzbXJykK/bmVXLVezuoazSx6upJXVr7ccOMYXyVmMuydXv47o7prWvoRMe+iM9hkIOOedEBjA50Y/4rm3jm25QepY+Jk0/bIO6+eZEW7ad1ZpQ/D8wbxZMb9/OCdyr3zo202LkB8irqeOmnA8wZNbhLb7o9nO25ZdYIntqYwuaDJUzrRlpTSoEeN0cdfm7WSWMK83HmjwPFmM1qnxaXWZeUy7L1e5kzyo9nzx/T4bUnhHjyzqYM9uZVWnxt9vFs3F3A67+mc8mkoTy+6Ng1kJZ0x5xwfkst5v41u4kN8exyC5+2VFXlt7RintmYQmqhnthgD978RywTQvq+eIOHsz0zwgezPjmP+88ehQIkZVewcHzHN1V7ys/Nka3ppRY/70Ciqk1r/XUaDU8sjmZYc7Dm7+YoSxrEgCQrHY9SVd/IVe/tQF9vZP0t01i+dALRQe68+Vs605//lX+8s42vd+VhMJp7dH6TWeX1Xw9ywZtbaDSa+fT6Kdx1Rni/VnKcFenHk4uj+SOtmPu+3N3aWb4zfx4o5uL/bUWrUfjypqldXsDvoNPy+KJoMktreev39N4O/YRX3WBk4+585o8NwMleS6S/G9ecFsZnOw+zPaP/0qCEbTCazNzxWRIbkvO438JBXItrTwvjkklDeeO3dFbvPGzRc/9nw17MqsojC6K6/JjLp4QS5OHEUxv3YzYf/7WsRVqBnkh/N6u9WQv1ccFgNJNX2Xepab+mFnH358lMCvXitUtjOv0709IaJr4f2hD8sK8An0H2PLm440DTUux1Gv578Xj09UbuX9O1v3dt7c2rZOm727nqvR3UG028eVksX940tV+CuBaLYwIp0jew9VApB4ur0TcYrRKMB7g7UqhvwNSN/1e25s8DJWzPKOO22SP4xykhTB3hQ4C7kwRxYsCSQK4Ng9HMDe/Hc7Comjf/EcvYIR6cGeXPiisnsunfs7hzTjgZJTXc8nEiU57+mac27ie9uLrL58+vrOMf72zj+e9TOSvKn29vn2616lndtWRSMHfMGcmXCTm88ENqp8euScjhqvd2MMTTia9untbtNSWnjvRhwbhA3vgtvV8X/9uCjbvzqWs0cWHc32vibp89kiGeTjzw1e4e31DoLwWV9RhNtjVmW9VoMnP7p0l8vSufB86O5AYrBHHQtHbmsUXRTBvhzQNf7WbrIcvcsf95fyHf7y3k9tnh3Uo5d7TTcu/cCPbmVbE2KbdLj1FVldRCPeFWKHTSIqylcmVJx/2ALCk+q4ybPownwt+Vt6+IO272w2A3R4Z6OfV5IKeqKlvSS5ky3KfPZioj/F2556wIftxXyOo21Vs7k19Zx92fJzP/1U3syavkkQWj+fHOGcwbE9Dvb/LnjPJjkIOOtYm5rf0AY4Mttz6uhb+7EyazSml1g8XPPRC0pPgGeTixZKL1+uwKYUkSyDUzm1Xu+SKZvw6V8twFY48pQR3o4cTtc0byx72ns/KqiUwM9WLFpgxm/9/vXPS/v/gqMYf6RlOH5/9uTwHzXv6T5JwKnrtgLK9dGoO788DKq7599kgumTSU139N54O/Mo/Zr6pNs4l3fZ7M5GFefH7jlB5Xn3zonFHYazUsW7+323dETyZfxOcwzMeF2DZ3V53tdTy+OJqDRdUs/2Pgz2oeKq7mtV8OMO/lPznl6Z/595e7+3tI/Sq9uJpvd+d3+nrRW01BXCLf7M7nwbNHcf106wRxLey0Gt64bALBXs7c8EE8h7pxg6s9dQYTj6zfy8jBg3rUrmTB2ECig9x44fvULj3P+ZX16OuNRFhpfRxAmG9zL7k+WCe3P7+Kq97bQYC7E6uuntTlNTwTgj3ZmVXep6/J6cXVFOsbmNbHbRmuOTWMyWFePLZhH4fLOg6u9fWNPP99CjOf/40Nu/K4fvowfr/ndK6aFtbnawk74min5awof77bU8Bfh0rxdLbrdSXq9vg3tyA4UStXfr+3kF05ldw+Z+SA+dkKcTzym9rs2e9TWJfUVM2ts4qAWo3CzIjBvLV0Alvun8W9cyMorKrnzs+SmfzUzzy6fi8pBVWtx9cZTDzw1W5u/DCeoZ7OfH3rqVwUN7Tf7+C1R1EUHl8UzezIwSxbv5fv9hS07jM2V457/vtUFo8P5L0ru/7moD1+bo7cfWY4f6QVs3F3wfEfcBLKKq1he0YZ508Ycszvy+kRgzlnTACv/nJwQM5qphdX8+rPB5j70h/M+r/feeGHNJzsNMyN8ufLhBzWJ+f19xD71MGial5pfj5m/9/v3PRRAjOf/41PtmdbfIay0WTmtk8S2bi7gIfOGcV104dZ9PwdcXeyY8WVE9FqFK5ZtZOKWsPxH9SBV385QE55HU8sju7RGyqNRuGBs0eRV1nPyi2Zxz0+taCp0EmEFXvi+bk64mhn/RYEWaU1XL5iO872Oj64ZlK31oBNCPWiWN9ATnnfpX9uPtg0gzt1eM/LtPeERqPwfxc1Va68e3XyMemCjSYzH/yVycznf+P1X9OZF+3PL3fP4P55o3B3Glg3YaEpvVLfYGRDch4xwZ5WeY/R0hT8RAzkTGaVF39MZZivC+fFWL9yuBCWMqCLnVTWNR7/IAtYtSWT//1+iMsmB3PzzK7fuR7s6sjNM0dw4/ThbD1Uyic7DvPxtmxWbskkJtiDBWMD+WhbFunFNdwwYxh3nxEx4O/y6LQaXr00hkvf3sbtnyby0bWTGR3oxm2fJPLT/iJumjmce86MsEgKzNJTQli9M4fHvt7LjAjfblfMPNF9GZ+DRoHzYtv/o7JswWj+SCvm4XV7eP/qSf1+c+BgkZ5vdhWwcXd+ayn3CSGePDx/NPOi/Qn0cMJoMnPR//7iwa92ExvswRDPE7c889HPh6I0VQZ8ZMFoQrydee2Xg9y/ZjfL/zjE3WeGc3Z0QK//X7UEcd/uaQrirj2tb4K4FiHeLixfOoFL397GDR/E88E1k7v9mnegUM/yPw5xfuwQJg/r+SzN1OE+zIoczOu/HuTiuKGdVnlt+X21ZiCn0ShNlSutGMgVVdWz9N3tNJrMfHzDlG7//5oQ/Pc6ub6qoLwlvYQhnk4WKyTWHUM8nXl0YRT/Wp3Mu5sOcf304aiqyo/7Cnnm2xQOldQwOcyL984Zxdghlk9VtKSpw33wdXWgWN9g0f5xbbVk4BT04TrPvrIhOY+0wurjriUVYqAZ0O+cs8tqeWzDPu4/O9Iq5f+hKeXx0Q1NFb0eWxTdozfDGo3C1BE+TB3hQ1mNgTUJOXy64zCPfb2Pwa4OfHjNZE4d2bd3G3vD2V7Hu1fEccFbf3HNqp2EeDuzJ7eSxxdFsXRKqMWuo9NqePLcaM57cwv//TGNh+ePtti5bV1T77hcTh3pS4C7U7vH+Lk5cu/cCB5et5f1yXl90n/waAcK9XyzO5+Nu/NJK6w+IliZF31s43edVsPLS2KY9/Kf3PlZEp9cd8oJ9UezvedjYogXjy4YzbwxAfi5/f18nB4xmJ/2F/HC96nc8nEiUYHp3Ds3kukjfXr0OtRoMnPrx4l8t7eAh+eP7lFKoiXEhXrx3AVjueOzJB74ajfPXzC2y/8eVVV5aO0eXBx0PHB27ytg3j8vkrNe+oNXfjnQacGU1AI9/m6OVk93D/V2Ia1Ib5VzV9Y2cvmK7ZRUN/DxdacwsgdBaYS/K4McdOzMKmNxH8xKmMwqf6WXMje6d20AeuP82CB+3FfAC9+nMdjVkY+3Z7M9o4zhvi68c3kcs0cN7vebZF2h1SgsGBvIis0ZxIZYp+qol7M99loNBVUn1hq5RpOZF39MY1SAG2dH900zdyEspdeBnKIoQ4H3AX/ADCxXVfXlo46ZCawDMpo3rVFV9bHjndt7kD0rNmewK6eC1y6N7fF6rI7EZ5Vx+6eJjBviwauXxKC1wCyTl4s91542jGtODSOlQE+gh9OATMM4Hu9BDqy6ahLnvbmFtEI9b/1jQq977rQnJtiTJRODWbklk/NjhzA60HprVGzJX4dKya2o49/zOn8ze+nkEL5IyOXxr/cxM3xwn6y7TCvU882upmDlQFHnwUp7hno58/jiKO78LJk3fkvvtMnzQKeqKmmF1a3B28GW5yPUi/8sjGJutH+Hz4eiKJwxuqlP5bqkXF78MY0rVmznlGFe3Ds38oh1kcdjMJq59ZMEvt9byLL5o7m6n4K4FotjgjhUUsMrPx9gmK9Ll/vWrUnIZVtGGU+fNwbvHpSFP9pIP1cunjiUD7dmceXUUEK82183lFqg73YD8Z4I9XHh55RCjCazRW9g1BlMXL1qB4eKa1hx5UTG93BGRqtRiAn2ID6rwmJj68y+vCqq6o3dahNhaYqi8NS5YzjrpT+547MkfAbZ88TiaJZMHGpzN5muPS0Ms6oSF2qdQE6jURjs5nDCzcit3plDdlktK66M69PWIEJYgiVm5IzA3aqqJiiK4grEK4ryo6qq+4467k9VVed358SB7k48ckkM//5yF/Nf/ZNXLomxWB59enE116zaSYC7I+9eEYeTvWX7mSmKwqgA2w5Kgr2dWX/LNBqMZqssnG7x77kR/LC3gIfW7uaLG6fKCymweudhXB11nDnar9PjtBqFp86NZuFrm3nmuxSePm+M1cb054Fi/rNhX7eClY6cGzOE31KLefnnA0wb4dNa+twWtFQ43Lgrn29255NeXIOiwKRQLy5fFMXcKH8Gd+P50GoUzosdwvyxgXyyPZtXfznIeW9s4YzRfvzrzIjjBhgGo5l/fpzAj/sKeWTBaK6a1r9BXIs754wks6SG575LJdTbhbPHdH6nu6LWwJMb9xMb7MHFcUMtOI5w1iXl8dx3qbx+Wewx+40mMweLq/skayLMx5lGk0peRb3FUgkNRjM3fRRPYnY5r18a2+t/R2ywJ6/+cgB9fSOuVm50vCW9BIApvUihtQTvQQ4sv3wCOzLKuOyUEJtN8w/0cOLRhV1v1dETAe6OJ9QaufpGE6/8fIDYYA9Ojxjc38MRott6fbtJVdV8VVUTmr/WA/sBi+VkLBgXyPpbpuHuZMc/3tnGW7+n97qiVpG+nitWbEerKKy6epJF7vyeqAI9nKwaxEFTQ9P7zx5FQnYFn1u4D5Utqqpv5Lu9BSwcF9ilhulRge5cc2oYn2zPZmemdXrLfbo9myvf2wHAY4ui2Hb/bD6/YQpXTA3tdhDX4vHF0QS4O3LHZ4no6/tmPWxvNZrMXLNqJ3Nf+pPXfj2Ir6sDjy+KYtsDs/nshilcPiW0W0FcW/Y6DVdMDeX3e2byrzPD2ZpeytyX/+Cuz5I6rKrXNoh7dAAFcdB0M+u5C8YSG+zBnZ8lkXS481meZ79LpbKukScs3EtssJsj1502jG9255OQfWxp/czSWgxGc7fbqPREqLflK1c+8c0+fkst5slzxzDvOMFyV8SFemJWOe7PyxI2p5cyYvCgHv+fsaTYYE9umDHcZoO4vuLv7kRB1YkTyH24NYuCqnr+dVaETaTQCnE0i+YNKIoSCsQA29rZPUVRlGRFUb5VFKVbt4xGDHZl3S2nMm9MAM98m8INH8RT1cM3fjUNRq5euYPSagMrrpzYYaqN6FvnxwYxKdSLZ75Loaym59XuTgTf7MqnvtHMhd2YlbhjzkiCPJx48Ks9NFqwCqKqqjz/fQr3rdnNtBE+fHXz1F4FK225Odrx0sXjyS2v45F1ey0wWutSVZVl6/byS0oRd50RzrYH5vDp9VNYOiWUwa6WeyPq4qDjllkj+fPfp3P99KYAZNb//caj6/dSrP97bYrBaObmj5qCuP8sjOLKARTEtXC007L88jh8XR24dtVOcivaT8mKzyrnk+3ZXDU11Crp1ddPH4avqwNPfbP/mBuBac2FTiL7ILWy5aaYpQqe7M2r5IPmtNFLJlmm79X4oR4oivUbgxuMZnZklPV52wHROwHujhRU1p8QbYOqG4y8+Vs6p47w6fOqqUJYisUCOUVRBgFfAneoqlp11O4EIERV1XHAq8DaTs5zvaIoOxVF2VlcXNy6fZCDjtcuiWHZ/NH8klLEwlc3sT//6Mt0rtHU9MZnf76e1y+LYZyVKjuJ7lMUhccXR1Ndb+SZb/f393D61RfxOYwYPIhxQ9y7/Bhnex2PLYoitVDP238essg4Gowm7vgsidd/TeeSSUN594o4i6daxYV6ceuskaxJzGVdF5s395f3NmfyyfZsbpwxnNtmj8TX1boz+R7O9tw/bxS/33M6F0wYygdbs5jx/K/83w+plFY3cPNH8fy0v5DHFkVxxdRQq46lN3wGObDiyok0NJq4ZuWOY2ZfW1qb+Ls5cscZ4VYZg4uDjjvnhLMzq5zv9xYesS+lQI9GgRGDrdcMvIWvqwMu9loyLBDIqarKYxv24elsz51zLPe8uTraEeHnavVALulwBXWNJqbIG2ib4ufmSIPRTEWtbWRRdOa9TRmU1hj411kR/T0UIXrMIoGcoih2NAVxH6mquubo/aqqVqmqWt389UbATlGUdl+9VVVdrqpqnKqqcb6+RzblVhSFq08N49PrT6Gu0cS5b2xmTUJOl8aoqioPfrWb39OKeWJxNLMiO197JPpehL8r15waxuc7c6yWIjjQHSquJj6rnAvb6R13PLNH+TE3yp+XfzpAdmnHDW67oqLWwNJ3trMuKY9750bw1LljrFY59tZZI5gQ4slDX+3ptDFvf/o1tYgnvtnHmaP9uLeP/+j7uzvy9Hlj+OmuGcwe5cervxxk0lM/89P+Ih5fHM3lFqwkay3hfq68flksB4qque2TxCN6563cksn+/CoeXTjaqmltF8UNYcTgQTz7XcoRs9ZpBXpCvV26lMbcW4qiEOLtQqYFUiu/21PAtowy7joj3OJFjuJCPUnMrjimt5olbUkvQVH6f32c6J4TpZdcRa2B5X8e4ozRfj0uDiTEQNDrd2ZK07vNd4H9qqq+2MEx/s3HoSjKpObrlvb0mnGhXnx962mMH+rBXZ8n8+BXu2kwmjp9zEs/HeDznTncNmuExVJQhOXdNnskge6OPLTWMimC9Y0mNu7OZ+XmDJtIBfmiuXfcuT0s/f3owijstBoeWrenx//e7NJazntzC0mHK3jlkhhunjnCqmsHdFoNL108HoA7PkuyeIPs3kor1HPrx4lE+Lvx34vH91sxnjAfF169JIavbz2VudH+PHfBWJaeEtIvY+mJ6eG+PLowil9Ti3nim6ZZ9/zKOv77YxqnR/hylhWq4ral02q4b24kGSU1fLI9u3V7aqG+T9bHtQjzcen1jFx9o4knN+4n0t+VJRMtVximxYQQT6objK2N0q1hS3op0YHufVJpV1hOS/XwQhtfJ/e/Pw5R3WDk7jOtkwUgRF+xxC32acBSYJaiKEnNH2crinKjoig3Nh9zAbBHUZRk4BVgidrLd9W+zf3ZbpwxnI+2ZXPhW3+RU97+3fzPdmTz8s8HuHDCEO60UuqOsAwXBx3LFkSRUqBn1ZbMHp2j0WTm15Qi7vwsiQmP/8jNHyXw6IZ9fLQt+/gP7kcms8qahFxmhPv2eA2av7sj/zoznD/Sivl6V363H5+YXc65b2ymrMbAh9dOZuG4wB6No7uGejnzxLnRxGeV89qvB/vkml1RWt3ANat24GSv5d0r4nAZAIUQooPcef3SWC6yYGXHvrL0lBCunhbGyi2ZvP9XJo9t2IfRrPKfhT3r4dlds0cNZnKYFy//1FSVsb7RRGZpTZ+0HmgR6uNMTnldr25Uvbspg5zyOpbNH22VEvlxIV4AxLdTHMYSag1GErPLmTpCZuNsjb+b7c/IFenrWbk5kwVjA4n0t+3q4kJYomrlJlVVFVVVx6qqOr75Y6Oqqm+pqvpW8zGvqaoaparqOFVVT1FVdUvvh958h3VeJP9bOoGM4hrmv7qJ39OKjzjm15QiHvhqD9PDfXnqvDFSlcgGnBXV1Fvrvz+mkd/FfjVms8rWQ6U88NVuJj35E1et3MHP+wuZPzaQj66dzPRwX574Zh/pxdVWHn3PbTpYQkFVfbeKnLRn6ZRQxg5x57Gv91FZ1/V1DN/tKWDJ8q24OOj48qapTArz6tU4umvR+CDOjQnilZ8PEJ/V/6m1DUYTN34YT1FVA29fHkegR/uN2UX3PHjOKGZHDuaR9Xv5dk8Bt80eabFS/MejKAoPnjOK0hoD//v9EAcKq1FV+jaQ83bBZFZ7nEZcWFXP678e5KwoP6Zaqf/aEE8nfF0dSLDSOrmdmeU0mlQpMGGDfF0d0CjYdC+5N35Nx2Ayy419cUKwrW6XHTgryp8Nt56Kv5sjV763nZd+SsNsVtmVU8HNHyUwKsCVNy6LtdoaH2FZiqLw6IIojOamxfwdUdWmn/ETX+9j6jO/sGT5Vr5KyOXUkb68fXkcOx6aw7MXjGXaCB+ev2AsjnZa7vwsyaJVHS3pi/gcPJztmD2qd71smnrLjaG0uoHnv0/p0mPe3ZTBTR/FMyrAjTU3T2W4r/ULP7TnsUVRBHk6cfunST2uTGsJqqrywJo97Mgs54ULx8kaCgvSahReuSSGqEA3Iv1due60YX16/bFDPFg4LpB3Nh3ijwNNN/76MpAb5ttcubKH6+Se/S4Fo0nlwbNHW3JYR1AUhQnBnuy00g2Vzekl2GkVJlqpcbWwHjutBl9XB5udkcutqOPjbdlcOGGI1VsrCdEX+j9PyEJCfVz46uZpPLh2Ny/9dID4rHL251fhPcieFVdOlN4wNibY25lbZ43ghR/S+DW16IhGnQeL9KxPymN9ch6ZpbXYaRVmhPty/9mRnDHaD2f7Y3/Wfm6OPH3uGG76KIFXfz7AXWcOrCpVlbWNfL+3gEsmDsVB1/uiC9FB7lw1LYwVmzM4L3YIscHtv2EymVUe/3ofK7dkMjfKn5eWjO+Tog8dcXW046WLY7jof3+xbO0eXloS0y/j+N8fh/gyIYfbZ49kQR+ll55MXBx0rL15Gkazir2u72+w3XNWBN/tKeDlnw5gr9MQ4tU3M4LQppdcSfdn5JIOV7AmbbBkjAAAIABJREFUIZcbZwy3+ixmXKgn3+0toKiq3uJ93v5KLyVmqGe7r9Vi4LPlXnKv/HQAgFtnj+znkQhhGSfUFJWTvZb/u3AcT507hm2HyjCaVVZdPcmiPZ5E37lu+jCG+brwyLq9pBdX8+Zv6cx7+U/mvPgHr/56kEAPJ545bww7HpzDO1dMZNH4oE7fGMwbE8D5sUN47deDVi+t3V0bduVhMJq5YILl1j3ddUY4/m6OPLBmd7uzkLUGIzd8EM/KLZlce2oYb1wW269BXIsJIZ7cNmska5PyWJvY9y0Jvt9bwLPfpTB/bAB3zJE/9tai02r67fdtqJczV0wNwWAyM3LwIKusM+uIl4s9ro66bveSa2o3sBdfVwdumTXCSqP7W2xI080fS79WVtY2sie3kinSP85m+bs5UGCDM3KHiqv5IiGHy04JJkhS5cUJ4oQK5KApJeTSycF8e8dprPvntH5LERO956DT8sSiaLLLapn9f7/z7HcpOOg0LJs/mm33z+bj605hyaRgPJztu3zORxeOJtDDibs+T6KmwWjF0XfP6vgcIv1diQ6y3MJrFwcd/1nYVDhmxaaMI/YV6etZsnwrv6Q0NZN+aP7ofqvG2J5/nj6cuBBPHlrbty0J9uZVcsenSYwNcueFC8fJmtoT2C2nj8TT2Y6x3ejXaAmKohDm0/0WBOuT80jIruCesyL6JMMkOtAde53G4oHc1oxSzCpMs9L6PmF9Ae5ONhnI/fenA9hrNdw80/o3QoToKydcINdiuO8gQrwl/9nWTR3hw7L5o7nnrAj+uOd01v5zGlefGtbjVB9XRztevGg82WW1PP51x+vv+tKBQj3Jhyu4oAe9447nzCh/zhztx39/SmsNiA4W6TnvjS0cKKxm+dK4AdlMWqfV8N+Lx6MAt3+a2CctCYr09Vy7aicezna8fXncgJidFNbj7mzHt7dP54GzR/X5tUO9u9eCoNZg5OmNKYwJcueC2CFWHNnf7HUaxg1xt3jlyr/SS3G008i6Uxvm7+6IvsFI9QC6GXo8+/Kq2JCcx1XTQvF1dejv4QhhMSdsICdOHFefGsY/Tx9hsTUhk8K8uHHGcD7dcZgf9xVa5Jy98UVCDjqNwuIe9o47nkcXRqFVFJat28OW9BLOe2ML9Y1mPrvhFOaM9rPKNS2hpSVBQnYFr/5i3ZYE9Y0mrns/noraRt6+PM7ia4LEwOTv7oirY9/3MQv1cSGvou64/U9bvPX7IQqq6lm2oG9nzmNDPNmTW0l9Y9fG2RWbD5YwMdSrX9ZGCstoaQpuS7NyL/6YiqujjhumD+/voQhhUfJKKk5Kd84JJyrQjfu+3EWxvqHfxmE0mVmTkMvMiMH4DLLOXcJAj/9v796j47zrO49/vnPT6DoaybYkWxfbiePEIbEjmYQkBUoDIaSFFAolnC7Lpd1sgZ62tLttWNouvW1boBS29BRSSrc9bbkFaELIQi7AcjEhsY3tXGwnthNrZNmWZF0saXQbzW//mBlFVnTXMzfN+3XOnHn0PM88z2+Sn6X5zu/3+37L9Tu37tR3j/fqP33uJ9pUE9bX33+Trm0u/G/E79izRW+5bov+9jvPaf8L2cmg55zTf7/3iA7HBvXJO/foZVtyO9UOpWfbhgolnZY1bfjM4Jg++/9O6o27N+vlW3NbEmRvW52mpp2ePDPkyfV6hsf1XM8I0yqLXENNcQVyBzsH9MjRHv3XV22nAD3WHQI5lKRQwKdPvn2PRiYS+v2vHtEa69Ov2g+e61Pv8ITetje706XedWObbthWp1fu2Kivvu8mteQwS99a/fEdV6s5WpG1kgT/+9ET+sbhbv3ebTv1+qsbPb8+MNdKMlf+xYNHZSbd/YYrs92sl2hvTX3Zs/8Fb6ZX/vjkBUnSTSQ6KWqZEbnl1nnNt49/+7jqK0N6z83b8t0UwHMEcihZOxqqdfcbrtR3jvXo3x/vzEsbvnIgprrK0CXlFbIh4Pfpi3e9Qv/83usVKS+ubySrw0F96s49OndxXH/w9ac8DbofONKtv3nkWb2lfYve92qm3CA3MvWrlspc+cQL/XrgyFnd9arL8pJlr76qTNs3VHqW8GTfiQuqCQd09WZGvYtZZkTufBGUINh3ok/7Tl7Q+19zuSopQ4V1iEAOJe1dN27VK3ds0J89cFSnekdyeu+B0Uk98kyPfnHPlpysFynmDIzXtUb127fs0P2Hu/WvP+n0ZM3O4digfvfLh7W3Laq/eMs1Rf3fB8WltiKk2oqgTi0SyCWTTn/8jafVWBPWr786t0XTZ2tvi+pg54AnX6DsO9WnV2yvl7+AMuRi5cJBv6IVwYIvCu6c08ceOq6mSFi/ckNrvpsDZAWBHEqaz2f62Ft3KxTw6YNfPjxvvbVsuf9wtyank3prR26y0BW797/mcl2/tU5/+B9PadcffUs/9/Hv6X3/ekB/8/Cz+r9PntXJ3hFNJ5f3YfPs0Jj+y7/s18bqMn32nR2eFGEHVmJrfeWiI3L3HujSU2cu6kO3X5nXwtkdbVH1j06uKMvmfGL9ccX6x1gft040FkEJgkeP9uinnYP6zVt2kIUY6xbjzCh5jZGw/tebr9EH/v2gPv2dE/rg667IyX3vPdClXU012rXZu9px65nfZ/o/7325vnusV8fPD+v4uYs6dm5Y33r6nDKDBWUBny7fVKWdjdXa2VCdem6sVmNNeGbELT6Z0K/9837FJ6f1r792g+qzlGQGWMy2DZV67NSFeY8Nj0/po98+rvbWWr1p9+Yct+xSe2cVBt++hrqs+072SWJ93HrRFAkX9IhcMun08YeOq62+gi9Lsa4RyAGSfv7aJj16dIs+/d0T+tmdG3VdazSr9zt27qKePDOk//nGXVm9z3pTEQro569t0s+raWbf2OS0TvSMzAR3x8+P6Ecn+vS1g2dmzqkJB2aCuhf64jp69qL+8d0v1xUN1fl4G4C21lfq6z89o7HJaZWHLh0t+LvvnlTfyIT+8V178z7l97KNVaoJB3Tg9IDetrdl1dfZd/KCNlaX6fJNqw8GUTgaI2Edjg3muxkLeuDJszp2bliffPseBf1MPsP6RSAHpH3kjqv1k+f79cEvHdI3f/OVWV0Yfe/+LgX9pjv2ZKd2XCkpD/l1TXNE1zRfmkBhMD6p4+eG9ez5YR1LP993qFujEwn94S/synqCGWAx2zamEp6c7h/VlY0vjsqfvjCqz//wef1Se7N2F0DRbJ/P1N4WXVPCE+ec9p28oJsuq897YApvNNaEdWF0UuNT0wU3bTExndQnH35WOxuq9cY8j2gD2UYgB6TVhIP6xC/v1p3/8Jj+7JtH9RdvuSYr95maTuo/Dp3RLVc2qK4ylJV7IJVQ4obt9bph+4tTuZxzGpuazuuaI0CSttW/mLlydiD35988qoDf9Hu37cxX015ib1tU3zveq8H4pGorVv4760TPiHqHJ5hWuY40pksQ9FycUGt9YZWz+drBMzrVN6rPvrODxDpY9xhvBma5YXu97nrVdn3h8U498sz5rNzje8d71Tcyybz9PDAzgjgUhK0bUh9+Z9eS23eiTw89c14feM3lMyneC0F7ep3cTztXN5XuRycy6+NIdLJeFHItufsOn9GOTVW6dVdDvpsCZB2BHDDH77zuCl3VVKO7v3ZEfSMTnl//3gMxbagq06t3bvT82gCKQ3U4qA1VoZnMlYnppP7kgWfUHC3Xr/5MYRUu3tNSK7/PVj29ct/JC2qpK1dLXWGN3GD1MoHcuQKsJXf6QlxXNdUwjRclgUAOmKMs4Nen7tyji+MJ3f3VI54WoL4wMqFHj/bozddtZgE2UOK21lfq+QupQO4LT8R07NywPnz7VQW35qgiFNCuphrtP92/4tdOJ50eO3VBNzMat65kRowLrQRBYjqps0PjauVLA5QIPkkC87iioVq/f9uVeuRoj77weMyz6953qFuJpNNbO1af/Q3A+rB1Q6qW3FB8Sp946Lhu2Fan217WmO9mzaujLarDsaEV19p8untIF8cTupH1cetKdTioqrJAwZUgODs0rumkU0tdeb6bAuQEgRywgPfctFU3X16vP33gmTUXw834yoEuXdsc0c5G0t4DpW7bhkr1DE/oL791VINjU/qjN+4q2OlgHW1RjU1N6+jZiyt63b6TqVp5BHLrT2MkrPMFNrWysz+15rQlyogcSgOBHLAAn8/08bftVtBv+uCXDq34m+i5nu4e0tGzF/U2kpwAUGpqpSR94fGY7nx5q67eHFniFfnTMasw+Er86ESfrmio0qbqwkneAm801hReUfBYJpBjaiVKBIEcsIimSLn+/M3X6FBsUNd+5CHd8ekf6vfuPazP/eCUfvhcn3qGx5e9hu4r+7sU8vuoawNA0ouZK6vLAvrdW6/Ic2sWt7m2XJsj4RUFcpOJpJ54oZ9sletUYyRccGvkYgNx+X02k4wFWO/Iww0s4Y27NysU8Oknp/p1/PxFfedYr768v2vmeF1lSFc0VOnKxhpd0VCtnY3VuqKhStXh4Mw5k4mk7jt0Rq+7umFVdZgArD/bN1RpQ1WZfuuWy7WhqizfzVnSSguDH4oNanwqSf24daopElbP8LgS00kFCiR5V2f/mDbXhgumPUC2EcgBy/D6qxv1+qtfTEJwYWRCx88P6/i5YT17fljHzg3rK/tjGp2cnjlnS215OqirlnNOA/EpascBmFEe8uvx/3GLfEVStLijLaoHjpxV9+CYNtcunUziRyf65LNUfU6sP42RsJJO6huZnCkQnm+x/jjr41BSPAnkzOw2SZ+S5Jf0OefcX845bunjt0uKS3q3c+6gF/cG8qG+qkw3VZVdMmXIOaeugbGZwO7ZdKD3g+d6NTXt1BQJ61U7qB0H4EXFEsRJ0t62OknS/tMDetMyArkfn7ygl22JKFIeXPJcFJ/GmheLghdKINc1ENdrr6IQOErHmgM5M/NL+jtJr5PUJekJM7vfOffMrNPeIGlH+nGDpL9PPwPrhpmppa5CLXUVumXWH5Kp6aSe7xtVdTggfxF9aAOA2a5sqlZ50K+Dpwf0piXW+sYnE/ppbEC/+jPbc9Q65FomeCuUdXKjEwn1jUyS6AQlxYtJxNdLOuGcO+Wcm5T0RUl3zDnnDkn/4lIek1RrZk0e3BsoeEG/T1c0VKspQl0bAMUr6PdpT0vtstbJPfHCgKamHevj1rHM37RCyVzZNTAmSWqO8rcWpcOLQG6LpNkVk7vS+1Z6jiTJzO4ys/1mtr+3t9eD5gEAAC90tEX1zNmLGp1ILHrevpN9CvpNL99al6OWIdeiFUGFAr6CqSWXKT3QyogcSogXgdx8c8Xm5mNfzjmpnc7d45zb65zbu3Ej64kAACgUHW1RTSedDncNLnrevhMXdF1rVOUhf45ahlwzs4KqJRcboIYcSo8XgVyXpJZZPzdL6l7FOQAAoIC1t6YLg7+w8PTKofiUnuoe0s3Uj1v3CqmWXGd/XOVBv+orKfGD0uFFIPeEpB1mts3MQpLulHT/nHPul/SfLeUVkoacc2c9uDcAAMiRSEVQOzZV6UDnwoHcj09dkHPSTZezPm69a4qEda5gplaOqbWuQqlE6UBpWHPWSudcwsx+Q9K3lSo/8Hnn3NNm9uvp45+R9KBSpQdOKFV+4D1rvS8AAMi9jraoHnzyrJJJN2/5hB+f7FN50K/dzbV5aB1yqbEmNSLnnMt7ANU1EFdLHYlOUFo8qSPnnHtQqWBt9r7PzNp2kj7gxb0AAED+dLRF9cUnYjrRO6IrGqpfcnzfyQu6fludQgEvJv2gkDVGwpqcTqp/dFL1VWV5a4dzTp39cb2C4vMoMfyWBQAAy9bRll4nN08Zgp6L43quZ4SyAyWiKZIpCp7f6ZX9o5OKT06TsRIlh0AOAAAs27YNlaqrDGn/PAlPfnzqgiTp5stJdFIKGtO15PJdgiCWriFHxkqUGgI5AACwbGam9taoDs6T8ORHJ/oUKQ/qqqaaPLQMudZYUxgjcp39mdIDrJFDaSGQAwAAK9LRFtXzfaO6MDJxyf59Jy/oxu318s+TBAXrz8bqMvl9lvcSBJli4C1RRuRQWgjkAADAiuzd+tJ1crH+uLoGxig7UEL8PtOm6rK8j8h1DcRVXxlSZZknOfyAokEgBwAAVuSaLREF/XZJPbkfneiTJBKdlJjGSDj/a+T6x9TM+jiUIAI5AACwIuGgXy/bEtGBWQlP9p28oE3VZbpsY1UeW4Zca6wJ6+zQWF7b0NkfV0uU9XEoPQRyAABgxTpaozpyZkgTiWk557Tv5AXddFl93gtDI7caI2GdTRcFz4fppFP34BilB1CSCOQAAMCKdbRFNZlI6unui3quZ0R9IxO6ibIDJacpElZ8clrDE4m83P/s0JgSSUfpAZQkVoUCAIAVmykM/sKAgv7UKBzr40pPQ7oEwfmhcdWEgzm/f+c6zliZr1FOFA8COQAAsGKbasJqqSvXgdMDmnZOrXUVal6HH6axuKZ0UfCzQ+Pa0VCd8/t39afW5623qZVVVVW6ePGizExlZWX5bg4KFFMrAQDAquxtq9P+0/167NQF3UzZgZLUFEmNyOWrllxsIC6fSU214bzcP1uqq6u1detW+Xw+DQ8PK5lM5rtJKEAEcgAAYFXa26LqG5nU8HhCN17G+rhStKkmNVp0Lk8lCDr742qKlCvoX38facvKytTa2qpNmzZpdHRU4+P5LfOAwrP+ej0AAMiJjtbozPaN2xmRK0VlAb/qK0N5Kwoe64+vu2mVs5mZotGotm7dKr/fz+gcLkEgBwAAVmVnY7WqygLa2VCtjdWs4ylVjZGwzuWpllxsYEwtdeu/hhyjc5gPyU4AAMCq+H2m/3brFdpYvb7WJ2FlmiJhdQ3kPpAbm5xW7/DEusxYOZ/M6FxlZaXOnTun4eFhVVZWyudjXKZUEcgBAIBVe/fN2/LdBORZQ01YB04P5Py+XQPp0gPreGrlfEKhkFpaWjQ4OKje3l4FAgGFw3yZUooI4QEAALBqTZGwBuJTGp+azul9YyUayEmXrp0LBoOsnStRBHIAAABYtcZ0LblclyCIpWvIlcIauYWEQiE1NzersbFRo6OjGhvLz1pF5AeBHAAAAFYtU0su15krO/vjCgd92lhV2ol2zEyRSETbtm1TKBRidK6EEMgBAABg1RpqUoHc+RzXkov1x9USrZCZ5fS+hWr26Fw8Hlc8Hlcikch3s5BFJDsBAADAqjXmaUQuVXqg9NbHLSYzOldRUaGhoSGNjo5qZGREzjmZmYLBoEKhEMHvOkEgBwAAgFWrKguoOhzIaS0555xi/XFdvzW69MklKBgMasOGDdqwYYOSyaQmJyc1MTGhkZERxeNxJZNJmZn8fr9CoZD8fn++m4xVWFMgZ2Yfk/RGSZOSTkp6j3NucJ7zXpA0LGlaUsI5t3ct9wUAAEDhaIqEdS6HUysH41MamUgwIrcMPp9P4XBY4XBYkUhEzjlNTU1pcnJS8XhcIyMjM0lSfD6fgsGgAoEAo3ZFYK0jcg9L+pBzLmFmfyXpQ5J+f4FzX+Oc61vj/QAAAFBgGmrCOc1aWcqlB9bKzBQKhRQKhVRVVaVNmzYpkUhoampK4+PjGhkZ0ejo6Ete5/P5ZkbxfD7fzAP5s6ZAzjn30KwfH5P01rU1BwAAAMWmKRLW8XPDObtfZ386kIsSyHkhEAgoEAiovLxc0WhUyWRSU1NTSiaTSiaTmp6engn2EomEEomEJicnF02mYmaXBHpmNjPKt9xtLM7LNXLvlfSlBY45SQ+ZmZP0WefcPR7eFwAAAHnUGClX78iEpqaTCvqzP0pDDbns8vl8KitbXlmHTLA3O+hLJpNKJBIz+2afN3vbOSfnnKanp+Wce8l5a5FJ8JJNmTZneHm/zLXTz26+c5YM5MzsEUmN8xz6sHPuvvQ5H5aUkPRvC1zmZudct5ltkvSwmR1zzn1/gfvdJekuSWptbV2qeQAAAMizxpqwnJN6hye0uTb7wVVsIK5oRVDV4WDW74XFZXOK5dxAae7PC+1bbH+2ZDlonDeyXTKQc869drHjZvYuSb8g6Ra3wH8x51x3+rnHzL4u6XpJ8wZy6dG6eyRp7969uf0/AAAAgBWbXRQ8J4Fcf5z1cSVgbnDEdMtLrSl8NrPblEpu8ibnXHyBcyrNrDqzLelWSU+t5b4AAAAoHJlacrlKeJIpBg6UsrWOg35aUrVS0yUPmdlnJMnMNpvZg+lzGiT90MwOS3pc0jedc99a430BAABQIF4ckct+LbnppNOZQYqBA2vNWnn5Avu7Jd2e3j4lafda7gMAAIDCFSkPqizg0/kc1JI7f3FcU9OORCcoeRR/AAAAwJqYmZoiYZ3NwdRKSg8AKQRyAAAAWLPGSG6KgsfSgVwrUytR4gjkAAAAsGZNkXKdy8HUytjAmMyUk+yYQCEjkAMAAMCaba5NjciNT01n9T6x/riaasIKBfgYi9LGvwAAAACs2bXNtUoknZ7uHsrqfaghB6QQyAEAAGDN2lujkqSDpwezep/YAIEcIBHIAQAAwAMbq8vUWlehA6cHsnaP8alpnb84QcZKQARyAAAA8Eh7a60OdA7IOZeV63cNpAqOt9aT6AQgkAMAAIAnOtqi6h2emAm4vBYboIYckEEgBwAAAE9cl1kn15md6ZVdmWLgrJEDCOQAAADgjSsbq1UR8utgltbJdfbHFQr4tLGqLCvXB4oJgRwAAAA8EfD7tKcltU4uG2L9Y2qJlsvns6xcHygmBHIAAADwTHtrVEfPDis+mfD82pQeAF5EIAcAAADPdLRFNZ10OhzzvjB4Z3+cRCdAGoEcAAAAPHNda60k7xOeDMWnNDyeUCsjcoAkAjkAAAB4qLYipMs2Vnqe8GSm9EAdNeQAiUAOAAAAHmtvjeqgx4XBO9OlB5qZWglIIpADAACAxzraohqIT+n5vlHPrhlLB3Kt9QRygEQgBwAAAI+1t6UKgx/wcHplbCCuSHlQNeGgZ9cEihmBHAAAADx1+cYqVYcDOtg56Nk1Y/1jrI8DZiGQAwAAgKd8Pkutk/NyRK4/TsZKYBYCOQAAAHiuvTWqZ3uGdXF8as3XSiadugbGqCEHzEIgBwAAAM91tEXlnHTIg+mVPcMTmpxOqpkROWAGgRwAAAA8t7slIjNvCoNnSg+0RFkjB2SsKZAzs4+Y2RkzO5R+3L7AebeZ2XEzO2Fmd6/lngAAACh81eGgdjZUe5K5cqb0ACNywAwvRuT+xjm3J/14cO5BM/NL+jtJb5C0S9I7zGyXB/cFAABAAWtvi+pQ56CSybUVBo8NxGUmbWFEDpiRi6mV10s64Zw75ZyblPRFSXfk4L4AAADIo47WqIYnEnquZ2RN1+nsj6uhOqyygN+jlgHFz4tA7jfM7IiZfd7MovMc3yIpNuvnrvS+eZnZXWa238z29/b2etA8AAAA5EOHR4XBu/rHmFYJzLFkIGdmj5jZU/M87pD095Iuk7RH0llJfz3fJebZt+D4unPuHufcXufc3o0bNy7zbQAAAKDQtNVXqK4ytOaEJ7GBuJopBg5cIrDUCc651y7nQmb2D5IemOdQl6SWWT83S+peVusAAABQtMzWXhh8IjGtcxfHqSEHzLHWrJVNs358s6Sn5jntCUk7zGybmYUk3Snp/rXcFwAAAMWhva1Wp/pG1T86uarXnxkYk3NkrATmWusauY+a2ZNmdkTSayR9UJLMbLOZPShJzrmEpN+Q9G1JRyV92Tn39BrvCwAAgCLQ0ZpaJ/fTVU6vjA2MSZJaCOSASyw5tXIxzrl3LrC/W9Lts35+UNJLShMAAABgfbu2uVYBn+lg54Buuaphxa/P1JBrYY0ccIlclB8AAABAiSoP+bVrc82qM1fG+uMK+X1qqA573DKguBHIAQAAIKvaW6M6HBtSYjq54tfGBuJqjpbL55svETpQugjkAAAAkFXtbVGNTU3r2LnhFb821j+mZtbHAS9BIAcAAICsWkth8M7+uFqirI8D5iKQAwAAQFZtjoTVUFO24sLgF8enNDQ2RekBYB4EcgAAAMgqM1NHW3TFI3IvZqwkkAPmIpADAABA1rW3RtU1MKaei+PLfs1MIBclkAPmIpADAABA1rWn18mtZHplrD9VDJyplcBLEcgBAAAg667eXKOQ37ei6ZWxgbiqwwFFKoJZbBlQnAjkAAAAkHVlAb+uaY7oYOfgsl8T648zrRJYAIEcAAAAcqK9tVZPdg1pIjG9rPM7++NMqwQWQCAHAACAnOhoi2pyOqmnuy8uea5zTl0DY2qpo4YcMB8COQAAAOREe2s64cky1sn1Dk9oIpGk9ACwAAI5AAAA5MSmmrCao+XLylzZSekBYFEEcgAAAMiZTGFw59yi58UGKAYOLIZADgAAADnT3hrV+YsT6h5avDB4poZcc5Q1csB8COQAAACQMx3pwuBL1ZPr7I9rU3WZwkF/LpoFFB0COQAAAOTMlY3VKg/6l0x4EqP0ALAoAjkAAADkTMDv0+6WyJIJT1KlBwjkgIUQyAEAACCnOtqieqb7osYm5y8MPplI6uzQmFpYHwcsiEAOAAAAOdXeGlUi6XSka3De492DY0o6MlYCiyGQAwAAQE5dly4MfmCB6ZWUHgCWRiAHAACAnKqrDGn7hkodPD3/iFym9ACBHLCwwFpebGZfkrQz/WOtpEHn3J55zntB0rCkaUkJ59zetdwXAAAAxa29LarvHOuRc05mdsmxzv64gn5TY004T60DCt+aAjnn3Nsz22b215KGFjn9Nc65vrXcDwAAAOtDe2tU9x7o0gsX4tq2ofKSY7GBuLbUlsvvswVeDcCTqZWW+hrllyV9wYvrAQAAYH3LFAafr55cV3+caZXAErxaI/dKSeedc88tcNxJesjMDpjZXR7dEwAAAEVqx6YqVZcF5k140tkfV3OUQA5YzJJTK80PqthTAAAIrklEQVTsEUmN8xz6sHPuvvT2O7T4aNzNzrluM9sk6WEzO+ac+/4C97tL0l2S1NraulTzAAAAUIR8PtOe1tqXjMiNTCQ0EJ9SKyNywKKWDOScc69d7LiZBSS9RVLHItfoTj/3mNnXJV0vad5Azjl3j6R7JGnv3r1uqfYBAACgOHW0RfWpR5/T8PiUqsNBSVKsP1N6gGLgwGK8mFr5WknHnHNd8x00s0ozq85sS7pV0lMe3BcAAABFrL01Kuekw7EX8+V1ZgI5plYCi/IikLtTc6ZVmtlmM3sw/WODpB+a2WFJj0v6pnPuWx7cFwAAAEVsT2utzKQDs6ZXZkbkmFoJLG5N5QckyTn37nn2dUu6Pb19StLutd4HAAAA60tNOKgrNlXr4KyEJ10DY6oqC6i2IpjHlgGFz6uslQAAAMCKtbdFdbBzQMlkKjVCrD+u5mj5S4qEA7gUgRwAAADypr21VsPjCZ3oHZGUWiPHtEpgaQRyAAAAyJvZhcGdc+oaGKMYOLAMBHIAAADIm20bKhWtCOrA6QH1jUxqbGpaLVFKDwBLIZADAABA3piZ2ltT6+QypQda6xmRA5ZCIAcAAIC8am+L6mTvqJ46k6onRw05YGkEcgAAAMir9tbUOrn7D3dLkpoJ5IAlEcgBAAAgr3a3ROT3mQ6cHtCGqjKVh/z5bhJQ8AjkAAAAkFcVoYCuaqqWJLXWkegEWA4COQAAAORdR3p6JaUHgOUhkAMAAEDetafryZHoBFgeAjkAAADk3fXb6hTy+7Rrc02+mwIUhUC+GwAAAAA0Rcq170M/p/rKUL6bAhQFAjkAAAAUhA1VZfluAlA0mFoJAAAAAEWGQA4AAAAAigyBHAAAAAAUGQI5AAAAACgyBHIAAAAAUGQI5AAAAACgyBDIAQAAAECRIZADAAAAgCJDIAcAAAAARYZADgAAAACKjDnn8t2GBZnZsKTj+W5HnkUkDeW7EQVgg6S+fDeiANAf6AsZ9AX6QgZ9IYX+QF/IoC+k0B/WT19oc85tnLszkI+WrMBx59zefDcin8zsHufcXfluR76Z2f5S7wsS/UGiL2TQF+gLGfSFFPoDfSGDvpBCf1j/fYGplYXvG/luAAoK/QEZ9AVk0BeQQV/AbPSHdY5ArsA55/hHiBn0B2TQF5BBX0AGfQGz0R/Wv0IP5O7JdwNQMOgLyKAvIIO+gNnoD8igLyBjXfeFgk52AgAAAAB4qUIfkQMAAAAAzJHTQM7MPm9mPWb21Kx9u83sx2b2pJl9w8xq0vtDZvZP6f2HzexnZ73m7WZ2xMyeNrOP5vI9wBtm1mJm3zWzo+n/j7+V3l9nZg+b2XPp5+is13zIzE6Y2XEze/2s/X9uZjEzG8nHe8HaeNwXvpX+ffG0mX3GzPz5eE9YHY/7wvfS+w6lH5vy8Z6wel71BzOrntUPDplZn5l9Ml/vCyvn8e8GPkMWsZX2BTOrT58/YmafnnOt4v/86JzL2UPSqyS1S3pq1r4nJL06vf1eSX+a3v6ApH9Kb2+SdECpwLNeUqekjelj/yzplly+Dx6e9IUmSe3p7WpJz0raJemjku5O779b0l+lt3dJOiypTNI2SScl+dPHXpG+3ki+3xePvPeFmvSzSfqqpDvz/f545K0vfE/S3ny/Jx6F0R/mXPeApFfl+/3xyH1f4DNk8T9W0RcqJf2MpF+X9Ok51yr6z485HZFzzn1fUv+c3TslfT+9/bCkX0pv75L0aPp1PZIGJe2VtF3Ss8653vR5j8x6DYqEc+6sc+5gentY0lFJWyTdodQvVqWffzG9fYekLzrnJpxzz0s6Ien69Osfc86dzWX74R2P+8LF9DkBSSFJLAIuIl72BRS/bPQHM9uh1JfDP8j+O4BXPOwLfIYscivtC865UefcDyWNz3Otov/8WAhr5J6S9Kb09tsktaS3D0u6w8wCZrZNUkf62AlJV5rZVjMLKPU/qkUoWma2VdJ1kn4iqSHzjyr9nJkOtUVSbNbLutL7sI540RfM7NuSeiQNS7o3641GVnj0e+Gf0lPp/tDMLOuNRtZ4+HfiHZK+5NJfx6P4rLEv8BlyHVlmX1jXCiGQe6+kD5jZAaWGSCfT+z+v1D+8/ZI+KWmfpIRzbkDS+yR9Salv1F6QlMhxm+ERM6tSagrcb88aTZn31Hn28Yd4HfGqLzjnXq/UVIkyST/naSOREx71hV9xzl0j6ZXpxzu9bSVyxeO/E3dK+oJXbUNurbUv8Bly/VhBX1jX8h7IOeeOOedudc51KPXL9WR6f8I590Hn3B7n3B2SaiU9lz72DefcDc65GyUdz+xHcTGzoFL/CP/NOfe19O7zZtaUPt6k1MiKlArqZ39r1iypO1dtRXZ53Recc+OS7ldqqgWKiFd9wTl3Jv08LOnfxZTLouTl7wYz2y0p4Jw7kPWGw3Me/m7gM2SRW2FfWNfyHshlMomZmU/SH0j6TPrnCjOrTG+/TqnRuGfmvCYq6f2SPpeHpmMN0tOc/lHSUefcJ2Ydul/Su9Lb75J036z9d5pZWXqq7Q5Jj+eqvcger/qCmVXN+iUekHS7pGO5eA/whod9IWBmG9LXDEr6BaWm8aOIZOHvxDvEaFxR8rIv8BmyuK2iL6xvucysotQv0LOSppT6tuRXJf2WUhlnnpX0l3qxSPlWpb4pOarUYtS2Odd5Jv0gK10RPpTKIOQkHZF0KP24XamMUo8q9Q3Zo5LqZr3mw0qN2B6X9IZZ+z+a7k/J9PNH8v3+eOS+L0hqUCoL7hFJT0v6W6W+fc/7e+SR875QqVRmwkxf+JTmyV7Io7AfXv6dSB87JenKfL8vHvntC3yGLO7HKvvCC0olWxxJf07cld5f9J8fM0ETAAAAAKBI5H1qJQAAAABgZQjkAAAAAKDIEMgBAAAAQJEhkAMAAACAIkMgBwAAAABFhkAOAAAAAIoMgRwAAAAAFBkCOQAAAAAoMv8fwY1ewLtxotoAAAAASUVORK5CYII=\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": [
    "### 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 can 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, because 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 and refitting the parameters, would go as follows (note again that the default for `append` does not refit the parameters, but we have overridden that with the `refit=True` argument):"
   ]
  },
  {
   "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.760863       NaN       NaN\n",
      "2000Q2       NaN  3.885850  3.885850  3.885850       NaN\n",
      "2000Q3       NaN       NaN  3.976371  3.976371  3.976371\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.000863       NaN       NaN\n",
      "2000Q2       NaN  0.304150  0.304150  0.304150       NaN\n",
      "2000Q3       NaN       NaN -1.206371 -1.206371 -1.206371\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.000863  0.304150 -1.206371\n",
      "2       -0.160790 -0.000863  0.304150 -1.206371 -0.151930\n",
      "3        0.417384  0.304150 -1.206371 -0.151930 -2.269410\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.227973\n",
      "2    3.263653\n",
      "3    3.305805\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 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": [
      "/usr/lib/python3/dist-packages/statsmodels/tsa/base/tsa_model.py:216: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  warnings.warn('A date index has been provided, but it has no'\n",
      "/usr/lib/python3/dist-packages/statsmodels/tsa/base/tsa_model.py:216: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  warnings.warn('A date index has been provided, but it has no'\n"
     ]
    }
   ],
   "source": [
    "mod = sm.tsa.SARIMAX(endog4)\n",
    "res = mod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this means is that you cannot specify forecasting steps by dates, and the output of the `forecast` and `get_forecast` methods will not 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": [
      "/usr/lib/python3/dist-packages/statsmodels/tsa/base/tsa_model.py:580: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.\n",
      "  warnings.warn('No supported index is available.'\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.8.2rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
