{ "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": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAEvCAYAAADW/SmEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3wbh303/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/pM5omadSLbdmWJXc7duzEsdNJCCxJNpBAFgg1dEJoCwv7we7CLrCU0CELoQUSdikhvTgkcWwncZy4yU2W5areNb2e749zzmhGU9RG0ki+f9ely/JIGo3k4zPznKdhYbEFP7p9HTbVFeVVdmc6fOfpZtxx3x482dSZdPuLLb042ePFe7ekf37UaAQuW+LCiyd6EYvJ/X+PHezAtcvLUWDUYXmlPWdBljqQLFfTonPh6SNdWFFpR3XhyGuSjQuLYNBqsPtkchmnOxCWA6EKO0qsRjgt+njlgnrxbnR5pNmgDqNKf558/ng3Ci16bKuXM56DvvRTZHe39CEUieHaNBUFa2uc6BoOomNI/v1KkoSvPXoULqsB79qcvvpIHfsPpM+0TWZPWyAcxfefbcbmxcV4/rNX4Du3rom3GY0XgzaakKMdyUNIVLVKgHa23wdJknCm34uFykQejUag2GoYM+B46nAnYhJww+r5OY4+k5pCC84rPW3qQJLUoE0eNjCZq76/fekMnBY93rI28+Lh+WIkaGOmbbZdVFuIcFSK92ek29MGjF0e2efJ7eoGnVaDDQsL8URTJwKRGD6bUHp9SV0xXm7tj1cMtPZ6oRHJ08VsJj2qnOaUCZKDvhA6hgJJpTlmvXbCy7XV/XZzMSu+UamIaB8KoCxH/WyqRa4CnOrxpH0xu/NEL6757g589Pevp13Ufq7fh089uA/v/dUe3PXgPvzrQ0349lPHce+Ok3hwz1k8cagDO5p78OH7X8MnHtiH2uICPPbJy/Dkp7ZifW0hPveng/jiXw5N+8S4u/64Dzf9ZFfS5NV0ej1BfOwPr6Oq0IzbN9Viz6n+cQ/DCISj+PD9ryEcieHed2+AzaTH1Q1lyn6r+dkzKEkSdihZts/+3wG0dI9kbH+96xRcViNuWJX5dce2+hL0e0M40jGMnS29GPCFceMa+fl0eaUdJ7rdWXuyxksN2oYDkZT9rbOhxx3E62cH8IblyRfSzQYt1i9wYueJ5AEtam/fiio7hBCoL7PFz5Px6ZEpmbbM5ZGxmPzvtnVpSXygWKZMmzoIb0WazNlata9NGfD01OFO7Dndj7uvrYctzRASAKhwJmTa0izXnkymrbXHi5gEvO2i6oxlmWNh0EYTEp8cObo8UjkAz/R50e0OIhCOYUHC1d2ShP1ZmTx+qAN1rpEFtheKmiILej0h+EKRkRG0o/5Dq0HyRK/odQ4F8OThTty2oWZahxDki6ICAzRCrsMPRWLsaZtF65VhJC8qT+yjn6xVY5VHqlNACy3pn1wnQ83M37axJt5/BgCXLC7GkD8c3zPU2uNBdaElZUVGfZk1pTzymPL3hoRzo2US5ZH7zg6iwmFKO3gp362scsCsnGdyWR4JyEOthgMRDKR50aYGOU8d7sRNP96Fkz3yi/JYTMLvXj6D6+7ZgWeOdKHPG8L+c4N49GA7fvrCSfzn48fwhb8cwkd+/zrefd8ePHe8B194YwP+/OFLsLTMBpfViN+9/2J85IrFeGDPWdz685fiF9jGQ5KkcWdMDpwbxPPH5cDi0/+7P+P/iWhMwl0P7sOgL4yf/NN63LimEpGYlPICOpP/ePQImtqGcc/b18aP/WuUaXvpSiQlScLe0/04eH4waarzXHK8y41udxB3X1MPg06DD9//GjzBCFp7PHjueA/eubk26xqcrUpf2wvNPXh4fzsc5pHMz/IKO8JRCSe6pz6MpGMoEF/dMd4LtM8d78ajB9un/L3TefZoFyQJeMOK1OzVlsUuHOkYjldSACO9fepF5mVlNjR3yhlwbzACIeRzYiL17+mqLZrah9DrCeHKhpJ4ieKgP32mbcAbglYj0k6CXF5ph0Grwf5zgwhFYvivJ46hvsyK2zbUZPzZE89fozNtJp12Upm2FuW8tLRsYtm1ROmfRYkyONIxjCqnOaXGt6rQDK1G4Gy/Lz6WeWHC1WmX1Zh15H+fJ4iXW/vxkcsXz/sSvtGqC+U0/PkBf8Zljw6zHguKLdh3dhDeYAQWg3Zcv6c/vHIGMUnCOzOUAMw3Wo1AUYExPsWPmbbZ47IasbDYgj2n5b6HTD1tBWP2tAVRaMntFNA3rarA3tP9+NQ1S5Nuv6TOBUDe47ayypEyOVJVX27DrpY+hKMx6JXHpVYhNI7OtGXJzkiSlPL/eN+5gTmZZQPkFQ7rFzixq6UvZ4u1VYkTJEev8TjcPoS6kgJ87caV+PgD+3Djj3bhS29qxN/2t+Hl1n5sXerCf/3jqqQSL0mS4A1FMeyX978N+cKoLrKgatS+P51Wg3++vgFra5z47P8ewJt/uBM/ePu6+Iv2TALhKK789vP4wNY6vP+y1AEXo/3w7y1wmPW457a1eN9vXsV/PHoE37hldcrnfe+ZZuxq6cO3blmNFZUORKIx2E06/P1YN96YJVsEyBchHthzFu+5dCGuThiLXltsQX2ZFc8e7Up5rH9+vQ2f/b8DSbeZ9Vrc956NuGTx3GhLUIeI3LqxGhsXFuKdv3wFn//TAZRYjdBrBW7flH0NTonNiOUVdjx9pAsnuty4cW1lvG9VDVCOtA/H35+stkE/1tc68frZQRxuH8b1K7P/ez57tAt3/u41OM16vGlVxaReOx3tGMZdD+7D565rSCktfDrNYA/VpUtc+M4zzdh9sje+julw+zBcVkN8N2p9mRXuYAQdQwG4gxFYjbqUx2gxZs60PX+8B0LIw2AKDFroNCJjpq3PG0KhRQ+NJvV3YNRp0Vhpx75zg/jtS6dxps+HX793Y9bnlAKjDnaTDsOBCOwpI/81CEwi09bS5YZGYErl6gzaaEIOtyf3bKj0Wg0qnSac6fOhplDZ0VaUkGmzGnEsy1jcZ450IRqT8MZVE+tnmw/UAO1snw9n+31wmPUpV3YAYFWVA48e7MCKrzwFIeSyggKjFgVGHaxGnfJ3HawJt/3ptfO4uqF00qn4uchlNcRfQDNom13rFxTi9OvyPqNMmbYCgxbeLNmoXk8Qrhz1s6lqiiz4xR0bU24vd5hQ5yrASyf78P7LFuFUrxcXL0odgrSszIZQNIYzfV4sKZVXAPzPjlZUOc1Jx5zZoEUoEkM0JsWvoKsGvCFcd88OfOSKxfGpdT3uIM71+3HHJQtz+vPOpI0Li+SgLcerRRJ3takrJVSH24exfkEhLl3iwiOfuAwfuf81fPEvh2Az6vCNf1yF2zbWpLxYVJezW406VGLsxezXrShH/Sds+Mj9r+GOX+3B3dfU4+NXLkn7IhGQx493DAXw+1fO4H0Z+qVUTW1D2H60C5++th5XNpTiw5cvxk+fP4krG0px3Qr5OdEbjOD7z57AvTtaceuGaty6Uc4S6LQabKsvwXPHexCLSRkfDwB8b/sJGHXa+NTDRFc3luF/drRiKGGZcDASxfeeacbKKjvuuroeA74QBrwhfPeZZjx7tCuvgrbf7D6NMrsJ16fpid9xogdLS62ocJhR4TDj89c34BtPHIMQwM1rq8a1BmdrvQs/f6EVAPCWNVXx2xe5CmDWa3G4fRhvm+LP0D7ox5XLSjHsj4y5l/Xl1j589Pevw6jToM8bwtn+kTkC43W614t3/XIPej1B/OtDTdiypBgWpVzRG5QHe7xz04K0x+6aagesRh12tfQlBW2NFfb456tDm5q73PAEImmfA7Jl2p473o3VVQ4UKxUzTos+4/TIfm8w607WdTVO/PHVczjWMYxt9SW4Yllpxs9VyX1t/pRzt0mvndT0yJYeD2qLUis3JoLlkXkkFzXR08kXkssJMk3bqS2SJ0ie7vNCpxGoTKgJdtmM6PMGMy5mfOxQBxYUWya9l2guiy/YHpCDtkx7Q/7lhkZ87aaV+OIbG/CJK5fg1g01uKK+FI3ldhRaDAhHYzg/4MNrZwfwRFMnfrX7NIb8YXxga91M/jizrsRmjE+qYtA2u9Sxz0BqWYxKzbRlKrnq84RQXDBz/46bFxdjz6l+tA364QtFkxZrq0YmSHqw51Q/bvnpboRjEu5990VJL3DUnzldtu38gB/d7iD+7ZEj+Onz8thndex6LoaQzJZrl5ehxGZEY4bnicmqKbJAqxEpu9oGvCG0DfrjvSxVTjP+90OX4Os3r8RTd2/D2y+uzVn1xiJXAf760S24aW0VvvtMMz7w270Z+9z+uk++WNHa4x2zrP1Hf2+BzaTDHZcuBADcfU09VlTa8cW/HEK3O4AnmzpwzXdfwL07WnHbhhr8+40rk77+qoZS9HqCWScMH+0YxiMH2vGeLQvTnhevaSxDJCYlTVh8cM85tA368c/Xy1mYWzfU4EOXL8aqKke89zIfdA8H8B+PHsFXHz6MyKhF4YFwFHtO9cdLHAHgQ9vq8MaV5ZAk4D1bFo7re1yufH253ZR0IUerEWiosMVLqicrGImi2x1EpdOMFWMMNzl4fhAf+M1e1BZZ8Is7NgDAhP89Oob8+KdfvIJoLIZv3rIKncMB/OS5kfHzO5p7EIrE0pZGAvLFgs11RditTNsNRWI40e1OyjYmBW3B7EHb6DLyAaWUOTG4cpj1GMqQaRvwhlFoyRy0ra1xwh+OwhOM4Es3NGb8vEQVTlPKjjZAGUQSiU64TLil24MlU1xdk1eZNkmS8PXHjmLHiR48ede2rFeMsvnL6+exotKBZWlSuvnqheYefPh3r+H5z12R86lbuXK0w42YJGd80qktKsBThztR5TSjpsiSlHousRoRjkoY8odROOpqyIA3hN0n+3DntroLrjQSkDNDZr0W5/rl8sh0mUwAqHSaJ1zmONaV1/kosY+NPW2zKzEjkq08MhKTEIzE0vZd9nlDWJnhnDMdLqkrxh9eOYtHDnQAQNKONtWSUis0AvjNS6ex/9wgqgvN+M17L07JaKv9Xb5Q6gsWdWDGsjIbvvnkMQQjUQTCMei1Im0z/VyxotKBV790Tc7vV6/VoKbQjFOjhmWoL5YTf2cmvXbKi9gzMRu0+O6ta7Cm2oGvPnIED+w5m7Lfq98bwvPHu3Hrhmr8dV8bHtrXlvEYPt7pxpOHO/HJq5bEM1wGnQb33LYWb/7hTlx/z4vo94bQUG7Dj25fh4sWpGZ+L68vgRDyFMnV1elXRXzvmWbYjDp8aFv6i3hra5woLjDg2aNdeMuaSvhCEfzw7y3YXFeEy5a4kj53/YJC/HrXaQQj0XFlDX78XAt63EF89S0rxvzcyXjw1XOIxCR0Dgfw3PGepDK/Paf6EYzEsK1+5GcQQuCet6/FJ3u8GZ9vR7toYSGKCgy45aKqlMzLiko7/ravfUrPt11DcvtIpdMEq0mHv+xrQ/dwAKWjXg+e6HLjjvv2wGnR43fv34QSmxEFBi1ePzOIm9dVj+t79XtDeNcv92DIH8YDH9yMVdUOvNzaj3t3tOJtG6qxoLgATx/pQqFFjw0LMl9A2rLEhe1Hu3Gu34chfxjhqJT0/7CwQC6VPN7pkYM2U7qgTb7NN6pEfseJHkgScMWykWDbaTFk7Gnr8wazvuZXS85v21g77tjgzq118YnoiUx6LSQJCEclGHTj+/eORGM41evFVQ1T25eZN5k2SZLwtceO4hc7T6G5y4PO4cDYX5RGOBrD5/90EPdsb87xI5xeL7f2wR+O4vUxpkbNJrUxdmVV+pPcgmIL+r0hNLUPYUFx8ouXbAu21dLIG8ao356vhBCoKTLjbL8X5wf8OS1lvNACNgDxUjqdRqRMfaKZtbTUCptRB51GwJhhUXRBhiutql53EMVZyl5yTR1S8oc98u6rdP0HJr0WC10F2HOqHysr7fjzhy9N+/9WXRwbCKVWUahB2zffuhq3rK/GPdtP4LcvncbySscFMTRoMha6CuI906rDo4YfzAQhBN6zZRE2LizEL148hfCo7M5jB9sRiUl4z6WLcMWyUjxysB3RDFUmP/z7CRQYtHjfqF6ypWU2fOUfViAmSfjymxrx6CcuSxuwAUCx1Yi1Nc6Mo/8Pnh/E00e68IGtdSn96CqtRuDKhlI8d6wb4WgMv959Gr2eID533bKUi6nra50IRWPjGowlSRJ++9Jp/OGVsxOepDoekWgMf3jlLC5dXIwyuxEP7Dmb9PEXT/TAoNVg06LkUk6jTjvugE39/L9/5nLcfU19yseWVzjgDk5t4qO6o61KybQBqYPHJEnCJx7YB61Gg/vfvwnlDhO0GoE1Nc5xZ9r8oSjuuG8PzvX78Is7NmBVtfz/5gtvbIBOK/Afjx5FOBrDs0e7cFVDWda+ry1KML/7ZG/84snoSqz6MhtOdI8j0zYqY/388R4UFRiSLkIUWvQZe9oGfOGs5ZG1xRbc//5N+Nc3jy/LBsh9ezetq0q5XX0um0hf25l+H8JRacIj/kfLi6BNkiR848lj+OXOU9ikpJ0TR7JOxJk+HyIxCbtP9mU8SeYjdRLZoTHqmKdT51AAN/14V0r5iaqpbQjFBYaMU8FGJkj6sHBUbbX6Qjrd2P/HmzpQU2TOGAxeCGoKLXjtzABC0VjG8kgaHzW75rIaL8igNZ9oNALrFhSiIE0DukrNwKUbRhIIR+EORma0zLXEZkR9mRXn+v0w67UZz3fv2FiLt2+swe8/sDmlekAVz7SlWRzrCcovPhxmPf77ravxjotr4QtFc7JUe75aWFyAU73epLKkw+3DqHCYsr5gmy4f2rYYbYN+PH6oI+n2v+5rw7IyGxorbLhxbSW6hoN4pbUv5etbuj147FAH7rh0Ydpg6vZNtdj3r9fiA1vrxhzEc9WyUhw8P4Rud+oF7+883YxCix7vu2xh1vu4prEMw4EI/n6sGz97/iSubihNGyiq5bvjuch8sseLruEgQtEY9p5JXcacyXhLz7Yf7UbncADvuXQhbt1Qg+ePd8cDIECeXrtxUWF8H9hUODMMRBoJsib/+k0d91/pNCdMi06+v8PtwzjW6cbd1y7FwoSLSetrC3Gs0w1fln2XqqePdOJQ2xC+e+vapP22ZXYTPnHVUmw/2oVvP3Ucw4FIxtJI1dJSK0psRuxs6cOR9mFYDFosGvXaTx37P+wPpw3ajDoNNALwjVr78urpflxSV5yU1XSYDWmDtmhMwoAvhKIs5ZEAcNlSVzyzNxVG5bw+kRUgakyzdK4HbZIk4dtPH8fPX2jFuzYvwI9uXw9g8kFbqzI+AcDvAAAgAElEQVRSc8gfjo+nnwvGCtrcgTA2/ed2PHW4M+3Hc+HRg+3Yf24wZfmkqqltGCuqHBlffCVeaR6daVOnCY2eIDnkC2NXSy9uWDm5yUfzRU2RJT7KuqZo7KZ4ykwN2tjPlh8+tK0OH7tyccaPq0Gb2oeYqE8ZJz2TmTZALpEE5MxOpsD/g9vq8I1bVmd9MRjvaUuTYfAoL1KsRh00GoH/vHklfvCOdWmHRJCsrqQAvlA06eLf4fbhWSsnvaqhFEtLrfjZC63xIONMnxevnx3EzeurIITANY1lKDBo8bf9yWPZJUnCfz91DCadNut0yfE+L17ZIPf+qGsDVK+e7pfbLy5fnHEnlWrrUhcMWg0+/6eDGA5E8Jk3LEv7eWV2E6qcZuw7Nzjm49rVIvc8aQSwqyU1cE3nWOcwVn/16XFlj+5/+QwqHSZc1VCK2zbWQALwx1flZehdwwEc63Rj29Lskz6nalm5DVqNmFJfmxq0lTtMsJv0WFhsQVNb8v09tK8Neq3Am0ZNCV2/wIloTMLBcSxmf+lkH2wmXdqBLe+7bCEWFlvw8x2tMOk1Y/7ehBDYsrgYu1t60dQ2hMYKe8r5clm5FYFwDKf7fGmDNiEECgy6pEqL4UAY5wf8KVk7p0WPoTSDSAZ9IUgSZuzCjUnJtE1kGIka0yye60Hbb186gx8/dxLvuLgW//aWFXBZDXCY9fE9K6PFYhL+tr8tpdlU1ZqQJdrZMr69JbNtOBBG26A8oaapbSjtFabXzgygaziI3dP4Mz19uAsAsPd06tWwYCSK5i43VmZ5ckwM1FIzbUrQNirTtud0P8JRKWn88IVIHfsPpO5oo4lRS3FzPXGQJmfLEhfu3JY5aFOD63QZAnWxdvEM9yaqU/HqSiY/mhlAvMwxbdA2atmsEAJvWVPJPsws1OcVtRrEH4oqw7FmrjQykUYjcOe2OhztGMYOZU/aQ/vaIQTwFmX5skmvxXUry/F4U0fSlfkH9pzDU4e78Mmrl+bk+F5RaUeZ3YjnEkokW7rd+PJfm+CyGvHucUwkLTDq4rsK/2FNZcahYwCwrtaJfePItO1s6UVtkQUXLSiMB3Bj2dXSB3cwgm88cSxrxq21x4OdLb24fVMtdFoNqgstuLy+BH989Swi0Vh8R+TWaQ7aTHotFpcUTHiPaqL2IT9cVkP8nLGi0oHDHSNBWDQm4eED7bhiWWlKVnZdjZL5HEeQu+tkLzaPymCpjDot/t8/LAcg/87Gk53cssSFPm8Ir50dSHvxZKkyjCQak9L2tAFyn2hilrBZSWI0ViT3njnNeniCkZRy5AGffHEvU8VDrqmZtoks2G7p9qDSYco4RXm8Zj1oe+xgB1ZW2fH1m1ZCoxEQQmBJqTVjpu3Fll7c9eB+bD/alfbjrT2eeHmLOtUm13yhCG780c6cBVBqlu2K+hIM+MJJqX2Vujy0uWtyGcix9HqC2HumH3qtwN4zAylTHps7PYjEpKwDAWwmffxKx+hMm8Osh14r0OtJbiLdd3YAOo3A6urZedLNF2qWUiPk8giaPDVYY6ZtblD3YqU776k9sDMdgG9aVAydRqB+ipO+Mk1GA+TySK1GwKSf9afhOSNxVxsAHO0cRkzCrA5uuXFtFcrsRvz8hZOQJAkP7W/D5kXFSefxm9ZWwR2IxLNgxzqH8W+PHMbWpa6Mg0EmSgiBK5eV4sUTvfAGI7hnezNu+P5OdA4H8K23rhp3eeA/rKmESa/B3aN2GI62vrYQ7UMBdA5lnj8Qicbw8sk+bFniwpYlLjS1D2HQl36QRCJ13P2eU/3xwCud379yFnqtiK8/AIDbL65F13AQzx3vwYsneuCyGtPuGcu1FZWOlOquHz57Arf9/KVxlXq2DQaSjpkVVXac6/fHpyW+dLIP3e4gblqb2mNVWGBAnasAr5/Jnvk81+/DuX4/Ls2yquGqhjL8vzcvx11XZ//3V6l9bZKEtNO/E8sBbVmGUSWeI48pr4mXladm2gCkZNv6PGpFxsw856uZtoks2D7R7Z5ylg2Y5aBNkiQ0d7uxutqZlFJdUmLNmGlT/zMfybDz62SPF3WuAly62IU9p/onVHM6Xi+39uHA+SE8kyFwBOQSifY0L0LSUQ/Qt14kT/5Jt5/jVSX7daI7866zqfj70W7EJOBdmxdiyB9G86jvo44SXjnGFc3aIgs0AklLTAH5CcVlNaYMItl3dhCNFfYLvvFeHftf6TTHl/XS5JSwPHJOKbOboNMItKVp4lcv8sx09qmwwICHPrYF79869lLkbLKN/Ff3Fl3IZeETVek0w6DVxCdIqpmN2QzaDDoN3n/ZIuw+2Yf7Xz6DU71e3DxqeMGli4vhshrxt/1t8IUi+NjvX4fdrMf3blub077bKxtK4QlGcOW3n8c920/g+pXlePYzl09oYt0t66vw6peuSTs1NdH6BWNndw62DcEdjGDLkmJsWeKCJMnBx1gOtQ1hW30Jqpxm/PdTx9MGPf5QFP+39xyuX1mRtGftqoZSlNmNuP/lM3jxRC+2LnXNSG/z8go7OocD8eqA3718Bt95phmvnOof14CSjkE/Kh0JQZvyWkvNtj20vw02ow5XN6bfL7a21on95wayBohqImPLqGmgo73vskXjnthb6TTHL6akGwZkM+njF+YyZtr0yZm2Y53DsJl0qHQk9xM7lAzj6MB/JNM2M4PHJpppi8UknOz2TnkICTDLQVuPJ4hBXzilMW9xaQF6PaG0V2TUxsxjGWqHW3s8qCux4rIlLgQjsWnZJaJe+TmapX75w/e/jk8+sG9c93dcOUCvbCiFTiNS6pJDkRj2nxuExaBFryeEfu/YV6rSefV0P7768OG0u9KePiKP6n+Psifm1VPJJZJNbUOwmXRj9ls1VtiwpNQKQ5pJcS6rMak8Uq7BHsQ6Nt7Hf68sjZw6l9WIG9dW4qqGsZdn0uzTagTKHaa0F7niV1BnodR1ZZVjyqUsWcsjg9Ep3/+FRqsRqC22xCdIHmkfgsM88qJwtrzj4lrYTDp89ZEjMOg0uH5Vcr+QTqvBm1dX4Nlj3fjc/x1Ea68X379tbc4vRly2xAWrUQe9VoNfvWcjfvCOdRP+HkKIMXvfADlIMeg0WYeR7FJeK1262IW1NU4UGLTYNUYFlDcYwckeD9bXOnH3tfU41DaUtpf/kQPtGA5E8M5NtUm367Qa3LahBi8096DfG0oa9T+d1AsHRzqG8dThTnzlb03xzNNYvX+SJKF90J+caVOHkbQNIxCO4smmTly/sjzjBe71tYXo9YRwrj9zgLj7ZB9cVuOUh2GMtm2pC0adBkvL0t9vvXJ75rUv2uRMW4cbjeX2lAta6jTo0cNIRnqf8zPT1j7khz8cxdIpVm4Asxy0nVBK/dQFfCo1Gk2XbVOvrB3tTA2Y+r0hDPjCWFxSgE11RdBqxLhrqCdCvc+jHe60VzV8oQiOdw7jtbMD8asu2RzrcKOh3AaTXoulZbaUYSSH24cQCMdw41q5Rr65a+LZtmhMwpf+egi/3n0ajzclT7ryhSJ48UQvrl1ehpoiM8rsRuw5nXwibmofxsrKzENIVF9603L8/gOb036sxJacaTvR7YY3FMXaGgZtNpMe5XZTzk+mFyKNRuD7b0+/z4jyU6XTnLE80mLQ5mTi12zImmkLpp+mRtktchXEyyPVISSzna20mfR45+YFiMYkXNtYBnuaoOemdVUIRWJ47FAHPnHVUlw6RrZjMgqMOjzz6W3Y/unL44NJpotBpxlzyfbOll6sqLSjqMAAvVaDixcVYfcYw0iOdAxDUvbB3ryuCktKrfj2081J08B7PUH8cucp1JdZkxZdq27dWAP1kBgrq5Qrav/fH145i08+sA+rq5148EObYdJrsG+M5MGwPwJvKIpK50hmyWU1otxuwuH2ITx7tBueYCTt+HmVuocs07+HJMlT1S9dXJzz/y+ffsMy/Pkjl2YMKOuV8tRM5zuzQQevErRJkoTjne60u9TU8sjRQVu/Z2YzbaYJZtpOKO1ecz7TpgYfo6PzJSXyP9bovrbhQBhn+nxwmPU41++HO5D8D6dOjqwrKYDNpMeaase4pxWNV9dwAM1dHiwotmDIH0ZHmnruox1ynb0kpU5yGk2SJBzvGjlAV1XZU4aR7FUCqNsvlpeFnphE0PbowXY0d3lgNepwz/YTSSfAHc29CEbkzfdCCGxcWIRXT/XHH0M4GsPRjuFxjeS3GnUZy9JcVkNSpm3/Wfnqkzo++EL3vx+6BJ/OMK2LaD6rdprTlkf2eYKzkmXLFXPWnrb0y2Ypu0WuApzu8yEUieFYpztvFpG/d8tCLCuz4Q6lWmW0NdUOrKi0Y8uSYnzyqumbEFrhMOdkvP14rK91oql9OO2LV18ogn1nB5MWc29Z4kJrb/bWkUNKpdGqKge0GoHPXFuPlm4PHtrXhlhMwv0vn8FV334erb0e3H1NfdoApLrQguuWl2N9rTOpdHI6OS0GVDnNeKKpE5VOM+57z0bYTXqsrnJi/xiZtsQdbYlWVtnR1D6Mv+5rQ5ndmDSif7Rl5TZ5yXaGoK2l24MedxBblmS+j8lymPVZyymXKYkZW4bzXYFBC79SHtk26Ic7GEFDRZqgzayUR47qaev3hWA16sa16D0XjPqJZdpOzp+gzQOnRR/vQVFVFZph0Glwsid5X5ja5Jkp49SqfP5ipRZ7yxIXDp4fTDsidLJ2Kun+DygjetOtFVBPOlajLuOyS1X7UADuQAQNSsPlqipHyjCSvWf6saDYgpVVdtiMugkPIwlHY/jeM81orLDjv/5xFVq6PXj04Mj44aePdMJh1uPihfIVq4sXFaFzOBCvwz7Z40EoEht3jXMmJTYj+ryheHnmvrODcFrk0bYkL390cBk0XYCqCs3oHA6kTAXr9YTm9DRFky5bpo3lkZOxsLgAoUgMO1t6EIrEZnSpdjalNhOeuntb2swPIJcd/vkjl+L+928ac+faXLG+thChSCzt66BXT8t7R7eMCtoAZK2AamobQqnNiFJlP+L1K8uxqsqB7z7TjJt/uhtffqgJKyodeOKurXjjqNH3ib7/jrUZq36my7paJ1xWI37z3ovjQ9nW1TpxuC19YKtK3NGWaHmlA609HrzQ3I23rKlMO/FRpS7Z3nc2fYCo/s4vXTwzmcdEly52YX2tE8sr0v9fNRu08CorUI4p8yoaylMvxjjimbbkFqF+b2hG9zSq5/XxZtpauj0oLjDk5DHOcnmkG/WltpQrJVqNQJ2rICXTppZG3rJeHtgxehjJyV4PDMrYV0A+QcQkpF1qOVk7W3pRXGDAjUqaOl1fW1P7MFxWA960qgI7muUnlkyOK2We6nQjNTBSh5FIkoS9pwewYUERhBBYWmbF8Qlm2v7y+nmc7vPhM9fW402rKtBQbsM9208gEo0hEo3h2aPduLqhNP5EslEJ3vYofW3qrpCpPjm6rMb4EkQA2H9uEGtrnLNe2kJEs6vKaUZMkisZEvV6gjPWpzAdNMp0SH+apbeeAMsjJ0MdevDoQbnMP18ybeNh0mvn1fNdfMl2mkBhV0svDFpN/PUEIGdcXFZD1qDtUNsQViVcIBZC4LPXLUPboB9tA37cc9ta/OGDm7BkjP4go047YxlH1bfeuhrbP70NtQkXotfWOBGKxnA0w/A8QO55AoAKZ3JWcGWlHTEJCEelrKWRqnW1ThztGE7bQ7v7ZB9qisxJ+3RnSrnDhL98dAvKHemzngUGXfzC1jHlNXG68ki7SQetRqQkYvq9oRkb9w9MPNN2otuTk8mRwCwGbZIkobnLnbFxMd3Y/8PKFZjV1Q7YTbqUYSStPV4sKLbEr0asq3XCpNfkrK9NkiTsbOnFpUtcsJv0qC2ypO2ta2obwsoqB65uLIU7GEm790yl/kdWa34bK+zQakS8r+1Urxd93hA2LpRPjvVlNpzoSt9Ll04wEsUPnm3Bmhonrm4shUYj8Klr6nGq14uH9rdjz+l+DPnDeMOKkelSy8pssJt02HtGDdqG5E33rqntLFKvmPd6QnAH5AmV6n4RIrpwqVeYR5dI9npCKLHN3fJIALAkvCBJ5AlGGLRNgvo89MzhLpj0mjGnHNL0KXeYUOkwpS3J23miFxctKEwKnDQagUsWu7DrZF/GeQAnezwpVT2X15fgj3duxrOfuRw3ravK28DXYtCl7lBTAttsfW1tg34YtBq4Rl2gWqH8HpaWWtOO0x9tfW0hIsqAt0TRmISXW/uwZRaybONhMWjhDcoXto51ulFTZM64iNth1qf2tHlDKJ6FTNt4ptNLkoSWbk9OSiOBWQzaut1BDAciKUNIVItLrDg34Ev6pSQ2HTdU2OOj8lUnezxJy1CNOi0uXlSMXWOMmG3pduO+nafGDISOd7nR4w5iq5LiX15hT7l6EghHcaLbg5WVDmxZ4oJBp8H2o5lLJI93ulHlNMcbl016LZaWWnFIyW6p/WwblKtV9WU2DPjCKfvOMnlwzzm0DfrxuTcsi5/orltRhpVVdvzg2RN4/FAHjDoNttWPLJ/UaAQ2LCyKZ9oOtw9huRJMToXa69brCeLQ+SFIkjymlogubFWFqbvaYjEJ/d65nWkD1HHWGUb+s6dtwsrsRpj1WrnvpXzqz0s0NesWFKYs2e7zBHGkYzht/9RlS4rR4w7GhzMkOtIuzwNYlaYVY1Nd8ZxsHyh3mFBuN2Xta2sfDKDCaUpZTVDpMOGSumJ8cGvduALVTJnPprYhDAciuCTLfrbZZDHoEIzEEI1JONbpTlsaqXKa9ak9bTNdHhkfRDJ2pq3XE8KQP3VK/mTNWtCWaQiJakmpFZI0skQzEI6iJeEKTGO5Dcc73fH+qHA0hrN9vng/m2rL4mK0dHuyLoD85c5T+PdHj+Cpw5n3rgEj/WyXLZWDtsYKO073eZP2SxzpGEZUWUJdYNThkrpiPHusK2NAeLzTnbL4cVWVIz6M5NXT/Si06LFYCUbVIHc8w0j8oSh+9FwLNi0qSjp5CiHw6Wvrcbbfhz+8chZbl7pSprNtXFiEkz1e9LiDONw+POV+NmAk09bjDsZH4K6tZtBGdKGrSpNpG/SHEZNmfrF2rpkN2pRypVhMgjfEnrbJEEJgYXwv1NwpjZyv0i3Z3q1cKE83uVHtqUpXAaVWGK2qzo8+xVxZm6XXDJB72hJ3tKmEEHjgzs1Jy8OzKSowYJGrICXzqf57zEY/23ioU3YHfCG09niyLkN3WPRJPW2SJM140GaMj/wfO9Om7lae85m25gzj/lXqD6iWSB7vdCMak+In6YYKOzzBSHxYxrl+HyIxKaVUQj1p7M6yG+TAOflE8bXHjmT9R9jZ0ou6koJ4KU9jhQ2ShKSM3+FRJ51rGktxps+H1l5vyv2FIjGc7PGk1O6uqnag3xtC+1AAe88M4CKlnw0Y2XcxnrH/9+5oRY87iM8kZNlUVy4rxdoaJ2IScO3y1MWbFy+Sr9j832vn4AtFc/LkmJhp23d2EHUlBfHGUiK6cJn0WrishnhvB4D4epDiOTyIBJBfkIwuj/QqF/oYtE3OIpfcl5MvQ0guZOuVapk/vnoOLzT34IXmHjx8oB02ky5txqymyIIFxZaMQVuJzYgy+8xMfJwp62qdONvvy7gCavSOtil9rxonXj7Zh+1HRpIFu0/2or7MmnGy92yzGOWg7eD5QcSk9ENIVE6zPqmnzReKIhiJzWjQptEIGLSacWXacjk5EpjFoO1ElxtFBYaMk8EWuQogxMiutiZlqbZ6km5U6nvVnjJ1cmRieSQglzAWWvTYmaGvzR+K4niXGxcvKsL5AT/u3dGa9vOCkSheae2Pl0YmPobEyUmH2oZQVGCIb3JXd6X8PU2J5MkeDyIxKSVoU7Nazx/vxqleb7yfDZADH4dZj+Y0pQWJfvvSaXxvezPevLoi7TQrIQS+/KZGrKpy4A3Ly1M+vrLKAYNOg9+9dCbpMU2F3aSDQatBjzuI/ecG2M9GRHGVTnP8IhyQGLTN7UybKU15pEfp32B55OQsLGamLV8sr7TDatThe9ubccd9e3DHfXvwzJEubFtaknFK5pYlLuw+2ZcyBbBp1BCS+ULdRZuuRDISjaFrOIAqZ24C1Q9dvhgldiM+8Nu9+KdfvIL95wbx6un+vM2yASOZttfPyL+fdOP+VU6LIT7MDpBLIwGgyDKzzxNGnWacmTZ51VZ5ji5EzNozRnOXO2uNp0mvRU2hJZ5pO9w+DLtJh2ql96G+zAoh5PGg160ojwd3i13J96nRCGxaVIxXWtMPAzncPoRoTMIHt9bBZTXgJ8+34JaLqlP2Zbx+ZhD+cDQp3V9daIbNpEuaIHmoTS4lVDNb1YUWNJTb8OyxLnxwW13SfR5XMnSNoxpM1f6x3+w+DWCknw2Qg636MmvW8sjf7D6Nrzx8GNcuL8N3b12b8fM2LCzCI5+4LO3HjDot1tY4sedUPww6TU6uEggh4LIasO/cIHo9IfazEVFcldOcNBlX7dsdvRJmrrEYtPEXFiq16Z6Ztsm5ZnmZ3PuS5cUdzQyjTosnP7UVXcPJWaR00/9U775kAR7YcxY/e6EVX3hjAwB5CElLtwfXr8w8xn+uWlUt75zbf24QVzcmVzZ1uYOISanj/idrWbkNT31qGx7Ycxbfe6YZN/14F4CZWzI+GWp7zutnB2DUaeIXZdIZPYgkHrTNYKYNAIx67bimR7YokyNzNTxnFjNtnoylkarECZLyEJKRYMhi0GFhcUE8YGrt8cJlNaQtt9tUV4S2QT/OD/hSPqZe+VhT7cC/3NAISQL+8/GjKZ+3s6UHWo3A5sXJvWGN5fb4YwiEozjR5cbKUVf/rmooxaunB1LGlB7rdEOvFSlTGdVhJM1dHhh1mpSl1kvLbGju8qTtk/vVrlP4ysOH8YblZfjx7eth0E3+n1jd29ZYboM+R3tlSmzG+DTNdTUM2ohIVuU0o33QHz+v9c2T8sh0g0jcAQZtU7G+thD3vWfjjC3TpeyqCy24aEFh0lu2Y7uh3I4b11Ti17tPoVtZ83G0I/MQkrnOYtBhWZktbV9bph1tU6HXavDuSxbi+c9diQ9tq8P6WmfeDiEBRjJt+88Nor7MlnW4kNOihzsQQUTZ6RkP2ma4IsOk14xrT1tLtydnQ0iAWQrawtEY3MFIvD8rkyWlVpzq9SIUieFYx3BKKURjhS2+06G114M6V/r7U8sD1WmIiQ6eH0KFw4RSuwnVhRZ85IrFeOxgB14aNXFy54lerK1xxqc8qpZXylMsYzEJxzvdiMSklJPO1Y2liMYk7GjuSbr9WOcwFpdY0wZEajnimmpnyhNTfakVQ/4wut3JV7Z+9/IZ/NsjR3DdijL8aIoBGwBsVH5vK3J4EnVZjYhJ8gGfrdmUiC4sVYVmBMKx+JNwnycErUbAOQcnxiVKN4iE5ZF0ofvUNfWIRCX86LkWAMCh88o8gHkYtAFyX9uBc4Px4Xmq6QjaVA6zHl+8oRF/+eiWvL5ApGbafKHomK8L1eeDYeXCV98slkcGx8i0uQPy6/TRAxKnYlaCNjWlONZyxMUlBQhGYtjR3INgJIYVozJODeV2nOn3wRuM4GSPN6WfLfHz7CZd2qDtwPlBrEmYYPjhyxejymnGv/z1EL755DF888lj+K8njuJg2xAuS5NebqywwReK4my/Lz75aHT/19qaQhRa9Hj6SPJ0ynSTI1XqiWvDwtS+LzVDmTiMpMcdxNcfO4LL60tyErABwEULClFTZMYVCesApkpthF1d5cxY705EF574rjblRUyvJ4iiAkPKGOy5xqxPHUTiYaaNLnALXQW4dWMNHthzFuf6fTjUNgyX1Ygy+9zOrGeytsYJdzASb+VRtcWDtvk1fGUiLAm7/LKV1QKIL9FW+yEHZi3Tph0z06ZW1+Vyh9ysvGpWf9DxZNoA4KH9bQBSJ0U1lMvTG/ec6ke/N5QxaNNqBDYm7B1TDXhDONPnw5qEMj2TXouv3bwS/d4QfvniKfzyxVP41c7TsBp1uH5l6sCO+DCSjmE0tQ3BadHH++4Sv/+bVlfgkQPt+OQD+9DvDWHIF0bHUAANGRYmXryoCEIgaX+aSl3ErU7gBICfPn8S4aiEr75lRc5KGa1GHV78/FV4w4rUn3uy1MEz7GcjokSjx/73ekIZB1XNJZZsmTYGbXQB++RVSyGEwD3bTyhDSOx5uzh7quJLtkcNI2kf9MNp0aesXbqQJAZto2c8jKbu6lN3tfV5Q9BrBWwzfC6VB5Fkz7Spg0oSF8xP1awcJYFwDFUFhjF7FdSU4vajXTDqNKgb1ful/uM+dqgDADKWRwJyX9uzx7rRPRxAqTLF5YCyNX5NTXIweOWyUhz4yhvG9bPUl9mgEXI9dlP7EFYm9N0l+so/rECJ1YQfPXcCu0/24q0XyXs3Ml1VaKywY++Xrkn7O3JZjSgqMMSHkXQOBXD/K2fwj+uqUvrj8o2aaWM/GxElqi5MzbTN9R1tAGA26OAPRxGLSfGsIYM2Innx9B2XLMAvd54CAFy3InX90HxR5yqAzaTDvrODuHXDyN619sFA2h1tF5LEgHXM8kilDHJIGUYy4A2h0GKY8WDfpNeOOT3SH5KDOrM+d0HbLJVHRjMu1U7ktBjgshoQCMfQWGFPKaerLjTDatTh6cOdAIDFWZr9Ll4kN2HuOT2SbTt4fghCTK2G2qTXoq7EigPnh3C8051xNL5eq8Fd1yzFwx+/DOUOE372wkkA2Q/QbEGtPKhEDtp+8nwLYjEJn7x66aR/jpmyrtaJOlcBNtXlb1MsEc08h1mPAoM2HrT1eYM5LSuZLeoTdiChlEYtjyxg0EYXuI9csQRmvRYxKTerhfKVRiOwtsaZMvY/lzva5io10+ayGsdM5jjjmTal93mGF2urjLqx97T5lH2cucy0zVJ5ZGzMyZEqNduWbh+LEAIN5eKbBDkAABaySURBVDYMByLQawVqCjMf+Csr7bAYtEklkgfODWJxiRU209Qa3ZdX2LGrpRfhaOoQktEaK+z460e34HPXLcPbLqqe9O6G+jIbTnR50Dbox4N7zuFtG2pQU2SZ1H3NpNXVTvz9s1fMyn8yIspfQghUOs0j5ZHu+VMeCSCpRNITisCo0+Sk95hoLisqMODObYuhU4Ka+WxdjRPHO4fxTMLi67ZBf852tM1V6oWtxnGs8HAqE+LVsf/93uCsvJ4cV6ZtGsojZ+UZIyZJWDrOoE3taxvdz6ZSSyRriyxZB1votBpctKAwvq9NkqSUISST1VhhR1SZCDSerJ1eq8HHrlyC/37bmkmndOvLrHAHI/jyXw8BAD5+1ZJJ3Q8RUb6oKjSjbdAPXygCfzg658f9AyMvSBLH/nsCEdg4OZIIAPCJq5bg2c9cHm9dma/etqEGi1wF+OBv9+L2/3kFr7T2wR2IXPCZNo1GoNJhwkULUgfvjWYz6SHESNA24AvPWtA2VqYt3tOWb+WRQojrhRDHhRAtQogvjOdr6se5t2AkaEvfnKgu16wbx0jNzXXFON7lxoA3hLZBv7zguWbq6Xj16oDdpENN0cz851OD3ueO9+DtF9ekLAMnIpprqpxy0NanLNaeHz1tSnlkwlVZTzDC0kgihUYjsCDLQuX5oqbIgic/tQ3/fuMKHOscxm33vgxgesb9zzVP3LUNH7ty7OSDViNgN+njkxn7PLOTaZMHkWTPtKkX6nIZtE35WUMIoQXwYwDXAjgP4FUhxMOSJB3J9nXjLY+8eV0VpCwLFxvK5WBuPHsQ4vvaTvfHM2NrcpCOX65k+1ZWpR9CMh3U359BpxnXgU5ElO+qCs0Y9IVxtt8HAPOiPDJTpo1DSIguPOri6xvXVuEnz7Xg0YMdOan4musclvG3KTktegz4QghHYxgORPK+PNKSZ9MjLwbQIklSKwAIIR4EcCOAjEGbTiPiuxbG4rQY8L7LFmX8+PIKO9bUOHH5OHaJra52wKjTYM+pfug0AgatJh70TUWJzYiVVXZc1VA65fsar6ICA5ZX2HF1YynK5nlJARFdGNSKAXWy73wI2uI9baMybQzaiC5c6uLrL97QONsPZc5xmvUY9IUxoOxqy9dBJGofsynPgrYqAOcS/n4ewKZsX2DKYarQbNDibx/bMq7PNeq0WFfrxJ5T/bAYtGistOekEVwIgUc/sXXK9zNRj33yshn/nkRE00UN2g6eGwIAFM+j8sikQSTByKSHUBERXcgcFgMG/WH0e2cxaFN62iRJylhh55+G8shc9LSle7RSyicJcacQYq8QYq8+GsjBt52cixcV43D7EA6eH8La6rk9XlYIMW8XURLRhaeqMDnTNh+mzKpBm29U0GblIBIiogkrtOgx5AvNatBm0svhU7Zsmz8chV4roM8yJHGicnFP5wHUJPy9GkD76E+SJOleSZI2SJK0obps9nZ0bV5UhJgk/zJz0c9GRES5UWozQacR6BgKwGbS5bQqY7ZY9HJwllge6WV5JBHRpDjN+tnPtOnk56ZgOHvQluvnsFwEba8CWCqEWCSEMAB4O4CHc3C/02JdbSH0Wjk7tZrNn0REeUOrESh3yGWD86GfDQBMBvlp1q8sWgUAd4CZNiKiyXBYDBjyh9HrDgKY3UxbIJJ5GIk/FM1paSSQg6BNkqQIgI8DeArAUQD/K0nS4ane73QxG7RYXe2EzahDnWv+j5glIppL1L62+TDuHwAshuRMWygSQzASg9XAoI2IaKKcZj0kCTijTBkutORvpi2XkyOB3AwigSRJjwN4PBf3NRM+c209OoYC0GjYD0ZElE+qCs3AKaC4YH5k2kaP/PcG5YwbM21ERBPnVNYDtPZ4YTfpctozNl7jzbTlujzygnzWuHSJa7YfAhERpVGtZtps8yPTptUIGHSaeKbNowZt7GkjIpowNWg71etF8SyV0ZvGmWkz5zjTNvPhKRERUQaVStA2XzJtgLyrTR3/rAZtNmbaiIgmzGGWL+idH/ChcAJLuXPJOM5MW67LIxm0ERFR3lDH/s+XnjZALpEcHbQVMNNGRDRhaqYtJgFFs3RxTy17DISzBG3hPBxEQkRElCvLym2wGnVYXmmf7YeSM2aDFj6WRxIRTZnTPJJdKyqYpUybTtnTlq08kj1tREQ0n5XaTGj6t+tm+2HkVFKmLcDySCKiyXIkBW2znGnLVh45DdMjmWkjIiKaRul62lgeSUQ0cTqtJn7Rq3gWdrQB48y0sTySiIhobjHpE8ojAyyPJCKaCrWvrXCWgrbxZNp8oShMzLQRERHNHRaDFoHRmTYu1yYimhSnMkFytjJt6sj/QIZMWzQmIRSJwaLP7XmeQRsREdE0Muu18IXlYM0TjMBq1EGjEbP8qIiI5qbZzrSpI/+DGTJt6lRJsyG3YRaDNiIiomlkNujgD8lXZD2BCAqMuS2ZISK6kKjDSGa7py1Tps2nVFawp42IiGgOkQeRKJm2UIT9bEREU6Bm2opmKWgTQsCg04yZacv1yH8GbURERNPIrNfCH45CkiR4AhFYTbOzW4iIaD5YVm5HTZE55yP1J8Kk02ScHulXgjZLjnuXebmPiIhoGpkNWsQkIBiJwROMwMZMGxHRpL1r8wK8c1MthJi93mCTXhvPqI0WL49kTxsREdHcofY1+ENR9rQREeXAbAZsgDyMJBjJkGkLsTySiIhozlFLePzhqDI9kuWRRERzmUmXOdMWmKbySAZtRERE08isBG2+kBy02UwsjyQimsuyZdo4PZKIiGgOSiqPDLI8kohorsuWaVMHkTBoIyIimkPUEpkBXwjRmMTySCKiOS7bIJJ40Jbj6ZYM2oiIiKaROkGsxx0EAFhZHklENKcZddkGkch7ORm0ERERzSFmvRykdStBG0f+ExHNbVkzbSE5mGN5JBER0RyiXm1VM20FDNqIiOa0rJm2cBQGnQZaTW7XEjBoIyIimkbqyP8ej1IeyaCNiGhOM+q1CIQzl0fmOssGMGgjIiKaVuqC1e7hAABw5D8R0Rxn1GkQzDKIhEEbERHRHMNMGxHR/GLSa7OUR8bi5/1cYtBGREQ0jfRaDfRawZ42IqJ5wqTXIBSNIRqTUj7mD0XiFRa5xKCNiIhompn0WrgD8hholkcSEc1tRp0clIXSZNv84WjOx/0DDNqIiIimnVoqo9MIGHV86iUimstMevk8nm7svz8UZXkkERHRXKQ2pRcYdRAit2OgiYhoZqmZtkAkNWjzhaIsjyQiIpqLzAa5JJJDSIiI5r4CoxyUeYORlI8FOD2SiIhobjIrpTTsZyMimvucFgMAYNAXTvmYP8zySCIiojnJwkwbEdG84TTrAQBD/tSgjeWRREREc5QpoaeNiIjmNqdFDtrSZdoCnB5JREQ0N6mlMlaWRxIRzXlOs1IeOSrTFo7GEI5KsDDTRkRENPeoQZuNmTYiojnPZtJBCGDIF0q6XV0BwEwbERHRHKSWR7KnjYho7tNoBOwmfUqmzR+Sgzb2tBEREc1BaqaNPW1ERPOD06JP6WnzK5k2To8kIiKag9SdPRz5T0Q0PzjNaTJtankkM21ERERzj9rfwPJIIqL5wWExpIz896nlkcy0ERERzT1mlkcSEc0rTrM+dRCJErRxeiQREdEcxJH/RETzi9OSpTySmTYiIqK5x6yXgzWO/Ccimh+cZj2G/GHEYlL8NrU8kj1tREREc9DGhYW4ZX01llfaZ/uhEBFRDtjNekgS4A5E4rcx00ZERDSHFVuN+M6ta2AxMNNGRDQfOC0GAMCgf6SvLcDpkURERERERPnBadYDQNKutnh5ZL5l2oQQXxVCtAkh9itvN+TqgREREREREeUjp0UO2hLH/vvVkf+63AdtuajT+J4kSd/Owf0QERERERHlPTVoS5wgGQhHYdJroNGInH8/lkcSERERERFNgMMs97Ql7mrzhaLT0s8G5CZo+7gQ4qAQ4j4hRGGmTxJC3CmE2CuE2NvT05ODb0tERERERDTzHGl62vzh6LQNnBozaBNCbBdCNKV5uxHATwEsBrAWQAeA72S6H0mS7pUkaYMkSRtKSkpy9gMQERERERHNJINOA4tBm1Qe6VfKI6fDmKGgJEnXjOeOhBD/A+DRKT8iIiIiIiKiPOc065MzbaHotEyOBKY+PbIi4a83A2ia2sMhIiIiIiLKfw6LAUMJe9r8oSgs+ukpj5zqvX5LCLEWgATgNIAPTfkRERERERER5TmnWZ888j8chV3pdcu1KQVtkiS9K1cPhIiIiIiIaK5wWvRo6fbE/+4PRVFmN07L9+LIfyIiIiIioglyWvQpg0hmbXokERERERERJbOb9RjyhSFJEgB1emQeDiIhIiIiIiK6EDnNBoSiMfjDUQDK9EgGbURERERERPnBaRlZsC1JEvzhKMyG6QmvGLQRERERERFNkNM8ErSFoxKiMYk9bURERERERPnCoWTahvxh+ENyiSR72oiIiIiIiPKE02wAAAz5Q/G+Nva0ERERERER5YnEnjY1aLMYGLQRERERERHlBYfa0+YPwxeKAGB5JBERERERUd6wGLTQawUGfWEE1PJIZtqIiIiIiIjygxACDrNB7mkLxQCwPJKIiIiIiCivOC16DPpGyiM5iISIiIiIiCiPOM16eeR/mCP/iYiIiIiI8o6aaQtweiQREREREVH+sSuZNl+Ie9qIiIiIiIjyjtNswKAvYbk2M21ERERERET5w2nRwxuKYtgfgRCAUTc94RWDNiIiIiIioklwWuQF251Dfpj1WgghpuX7MGgjIiIiIiKaBIdZDto6hgLT1s8GMGgjIiIiIiKaFKfFAADoGg5MWz8bwKCNiIiIiIhoUpzMtBEREREREeUvtTwyGIkx00ZERERERJRv1EEkwPTtaAMYtBEREREREU2KzaSHOjCSmTYiIiIiIqI8o9UI2E1yto2ZNiIiIiIiojyklkgy00ZERERERJSH1AmSzLQRERERERHlIYeyq41BGxERERERUR5Sx/5bWB5JRERERESUf9TySBODNiIiIiIiovwTH0TC8kgiIiIiIqL8w/JIIiIiIiKiPOZUBpGYmGkjIiIiIiLKPxz5T0RERERElMcWFFsgBFDpNE/b99BN2z0TERERERHNc0vLbNj7pWtQbDVO2/dgpo2IiIiIiGgKpjNgAxi0ERERERER5TUGbURERERERHmMQRsREREREVEeY9BGRERERESUxxi0ERERERER5TEGbURERERERHmMQRsREREREVEeY9BGRERERESUxxi0ERERERER5TEGbURERERERHlMSJI0899UCDeA4zP+jecnB4Ch2X4Q84QLQO9sP4h5gsdlbvCYzB0ek7nD4zJ3eFzmDo/L3OFxmTtjHZcLJEkqGc8d6XLzeCbsuCRJG2bpe88rQoh7JUm6c7Yfx3wghNjL4zI3eFzmBo/J3OExmTs8LnOHx2Xu8LjMHR6XuZPL45LlkXPfI7P9AIjS4HFJ+YbHJOUjHpeUj3hc5iEGbXOcJEn8j0V5h8cl5Rsek5SPeFxSPuJxmZ9mK2i7d5a+L1E2PC4p3/CYpHzE45LyEY9Lykc5Oy5nZRAJERERERERjQ/LI4mIiIiIiPJYToI2IcR9QohuIURTwm1rhBAvCSEOCSEeEULYldsXCiH8Qoj9ytvPEr7mNiHEQSHEYSHEt3Lx2OjCNZHjUvnYauVjh5WPm5TbeVxSzkzwfPlPCefK/UKImBBirfIxHpeUMxM8LvVCiN8otx8VQnwx4WvuEkI0Kcflp2bjZ6H5YYLHpEEI8Svl9gNCiCsSvobnSsoZIUSNEOI55dx3WAhxl3J7kRDiGSHECeXPQuV2IYT4gRCiRTkO1yfc1zeV82WTEOK2Mb+5JElTfgOwDcB6AE0Jt70K4HLl/fcB+A/l/YWJn5fw+cUAzgIoUf7+GwBX5+Lx8e3CfJvgcakDcBDAGuXvxQC0PC75luu3iRyXo75uFYBW5X0el3zL6dsEz5e3A3hQed8C4LTy3L4SQJNymw7AdgBLZ/tn49vcfJvgMfkxAL9S3i8F8BrkxATPlXzL6RuACgDrlfdtAJoBLAfwLQBfUG7/AoBvKu/fAOAJAALAZgCvKLe/CcAzyrmyAMBeAPZs3zsnmTZJknYA6B918zIAO5T3nwFwyxh3UwegWZKkHuXv28fxNUQZTfC4fAOAg5IkHVC+tk+SpCh4XFKOTeF8+Q4ADyjv87iknJrgcSkBKBBC6ACYAYQADANoBPCyJEk+SZIiAF4AcPN0P3aanyZ4TC4H8Kzydd0ABgFsAM+VlGOSJHVIkvS68r4bwFEAVQBuhHxRAMqfNynv3wjgt5LsZQBOIUQF5GP2BUmSIpIkeQEcAHB9tu89nT1tTQDeorz/NgA1CR9bJITYJ4R4QQixVbmtBUCDUj6pg/zDJn4NUS5kOi7rAUhCiKeEEK8LIT6v3M7jkmZCtvOl6jaMBG08LmkmZDou/wTAC6ADchbj25Ik9Sufv00IUSyEsEC+wszjknIp0zF5AMCNQgidEGIRgIuUj/FcSdNGCLEQwDoArwAokySpA5ADO8gZX0AO6M4lfNl55bYDAN4ohLAIIVwArsQYx+Z0Bm3vA/AxIcRrkNOHIeX2DgC1kiStA/BpAH8QQtglSRoA8BEAfwTwIuRyi8g0Pj66MGU6LnUALgPwT8qfNwshruZxSTMk03EJABBCbALgkySpCQB4XNIMyXRcXgwgCqASwCIAnxFC1EmSdBTANyFnQJ6E/KKExyXlUqZj8j7IL4b3ArgHwG4AEZ4raboIIawA/gzgU5IkDWf71DS3SZIkPQ3gccjH6gMAXsIYx6Zuko91TJIkHYNccgYhRD3k2k1IkhQEEFTef00IcRJylmOvJC/ze0T5mjshPykQ5Uym4xLyyf4FSZJ6lY89DrmW/lkelzTdshyXqrdjJMumfg2PS5pWWY7L2wE8KUlSGEC3EGIX5FK0VkmSfgngl8rX/CfkcytRTmR5bRkBcLf6eUKI3QBOKB/juZJySgihhxyw/V6SpL8oN3cJISokSepQyh+7ldvPIzmDVg2gHQAkSfo6gK8r9/kHKMdsJtOWaRNClCp/agB8GcDPlL+XCCG0yvt1AJYCaB31NYUAPgrgF9P1+OjClOm4BPAUgNVKmloH4HIAR0Z9DY9LmhZZjkv1trcBeDDD1/C4pGmR5bg8C+AqZSpaAeTm+mOjvqYWwD9i1MUGoqnI8trSohyLEEJcCznLxudwyjkhhIB8YeqoJEnfTfjQwwDuUN6/A8DfEm5/t3K+3Azg/7dzhyoVBFEcxr+D2GyKj3F9AEF8CZsGsQj3DSzCTSbBYDYKgjYvaPUFRE0ipmsWDRbFY5gJIigsrjDg90vDzjK74TDsf3d2nmqwm4qI2TrmABgA5z9du5cvbRFxCCwDcxExAbaBmYgY1lNOgIPaXgJGEfFGeduxWdfCA+xFxEJtjzLzto/70//UpS4z8zEidik7UyUwzszTep51qd50nC+hzJmTzLz/MpR1qd50rMv92r6hLP05yMyr2ndcH0RegWFdniZ11rEm54GziHgHHoC1T0M5V6pPi5T6uo6Iy3psC9gBjiJig/Jia6X2jSn/994BL8B6PT4NXJQMyDOwWr8YfyvqtpOSJEmSpAb95UYkkiRJkqRfMrRJkiRJUsMMbZIkSZLUMEObJEmSJDXM0CZJkiRJDTO0SZIkSVLDDG2SJEmS1DBDmyRJkiQ17APyarAttLEi8gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "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: Tue, 17 Dec 2019 AIC 951.427\n", "Time: 23:40:41 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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd1gU5/bA8e8sS5cO0pRiARQsIGrURI0aS2JLN8X0ftNzk5tqctPbLze9mMRoemJiLInpMUWNjWYFFCnSOyx12d35/UEJKiBlF1g9n+fhAWZmZ14XXPbMe95zFFVVEUIIIYQQQghhPTR9PQAhhBBCCCGEEF0jgZwQQgghhBBCWBkJ5IQQQgghhBDCykggJ4QQQgghhBBWRgI5IYQQQgghhLAyEsgJIYQQQgghhJXR9vUAOuLt7a2GhIT09TCEEEIIIYQQok/ExcUVq6rqc+z2fh3IhYSEsGvXrr4ehhBCCCGEEEL0CUVRMtvaLqmVQgghhBBCCGFlJJATQgghhBBCCCsjgZwQQgghhBBCWBkJ5IQQQgghhBDCykggJ4QQQgghhBBWRgI5IYQQQgghhLAyEsgJIYQQQgghhJWRQE4IIYQQQgghrIwEckIIIYQQQghhZSSQE0IIIYQQQggro+3KwYqirADmA4WqqkY1bfMEvgBCgAzgIlVVy9p47JXAw03fPqmq6qruD1sIIYQQQgghQFXVls/NH62/b+u4Y7d1Zn9b3wOYTKZObWt9jrbO08E4lLaO61IgB6wEXgc+bLXtfuBXVVWfVRTl/qbv/9P6QU3B3qNALKACcYqirG8r4BNCCCGEEEKcnFRVxWQytXw+9muTyYTBYDjqM3DU/uYAp/n7tq6hKG3GPu3uP9FjTnRcZ7d1Zl9r9fX10E4WZZcCOVVV/1QUJeSYzYuA6U1frwJ+55hADpgD/KyqaimAoig/A3OBz7pyfSGEEEIIIUT/YzQaMRqNGAwGDAYD9fX1LV8bjcaWwKxZRwGRRtMYt2g0mpZjFEVBURS0Wu1R33c2ILJWHc3cdXVGri2+qqrmNV0oT1GUgW0cEwgcafV9dtM2IYQQQgghRD+nqmpLoGY0GtHr9ej1eurq6mhoaMBoNKIoSkuAZmNjg0ajafnQarXY29uf9IFXbzJHINcZbf3E2gwvFUW5AbgBICgoyJJjEkIIIYQQQhyjoaGB+vp66uvrqaurQ6/X09DQAPwzk9YcrGm1WhwdHSVA6wPmCOQKFEXxb5qN8wcK2zgmm3/SLwEG0ZiCeRxVVZcDywFiY2Pbn0sUQgghhBBCmIVer6e2tpaysjLq6+uPmlVrnk0T/Ys5Arn1wJXAs02f17VxzI/A04qieDR9Pxt4wAzXFkIIIYQQQnSRqqrU19dTXV1NRUUFDQ0NaDQa7OzscHFx6evhiU7oavuBz2icWfNWFCWbxkqUzwJfKopyLZAFXNh0bCxwk6qq16mqWqooyhPAzqZTPd5c+EQIIYQQQghheSaTibq6OqqqqqisrMRkMqHRaLC3t8fBwaGvhye6qKtVKy9pZ9fMNo7dBVzX6vsVwIoujU4IIYQQQgjRbUajkbq6OiorK9HpdKiqilarxcHBoaU6pLBOvVXsRAghhBBCCNELTCZTy6xbdXV1S9l+Z2dnKUpyEpFATgghhBBCiJNEQ0MDOTk51NfXY29vL+vdTmISyAkhhBBCCHESqKmpIScnBxsbGwngTgESyAkhhBBCCGHFVFWlvLycgoICnJyc0GrlLf6pQH7KQgghhBBCWCmj0UhhYSEVFRUMGDBACpicQuQnLYQQQnTDnuwKHlu/D1VV+3ooQohTlF6v58iRI1RVVeHq6ipB3ClGftpCCCFEN3wdn83KrRnkV9b19VCEEKeg6upqMjIyMJlMODs79/VwRB+QQE4IIYTohpR8HQDpxdV9PBIhxKlEVVVKS0vJysrCwcFBGnmfwiSQE0IIIbohtaAxkMsorunjkQghThVGo5Hc3FyKiopwdXWVoianOPnpCyGEEF1UpKunpFoPQHpxVR+PRghxKqivrycnJweTySStBQQggZwQQgjRZc2zcYoC6TIjJ4SwMJ1OR25uLnZ2djg5OfX1cEQ/IYGcEEII0UXJTevjYoI8yCiRNXJCCMtQVZWSkhKKi4txdnbGxsamr4ck+hFZIyeEEEJ0UWq+Di9nO2KDPcgqqcFokhYEQgjzMhgM5ObmUlJSgouLiwRx4jgSyAkhhBBdlFygI8zXhRBvZ/RGE7nltX09JCF67IE1u/lpX35fD0PQuB4uKyuLuro6XFxcUBSlr4ck+iEJ5IQQQoguMJlUDhboCPdzIcSrsXeTpFcKa1eoq+OzHUf4cV9BXw9FAGVlZaiqiqOjY18PRfRjEsgJIYQQXZBdVkuN3ki4nwuh3k2BnPSSE1YuIascgPxKmV3uD1RVRaORt+miY/IbIoQQQnRBSlPFynA/F3xd7XG0tZHKlcLqxWeVAZBXUdfHIxFCdJYEckIIIUQXpORXAhDm27huJdjLSVIrhdVLyGyakauoQ1WleI8Q1kACOSGEEKILUgqqGOThyAD7xg4+Q3ycJbVSWLUGo4ndOeXYazXU6I3o6g19PSQhRCdIICeEEEJ0QWq+jnBfl5bvQ7ycySqtwWA09eGohOi+5DwddQ0mpof7AI2zckKI/k8COSGEEKKT9AYTaUVVhPu1CuS8nTGYVHKkBYGwUs3r484e5Q/IOjkhrIUEckIIIUQnpRdXYzCpRwVyzZUrD0t6pbBS8Vll+LraExPkAUCBBHJCWAUJ5IQQQohOSm4qdHLUjJyXtCAQ1i0+q4yYIA98XR0AmZETwlr0OJBTFCVcUZTEVh+ViqLcecwx0xVFqWh1zLKeXlcIIYTobakFOrQahSHeA1q2eQ+wY4C9VgI5YZWKdPUcKa0lOsgdO60G7wF20ktOCCuh7ekJVFVNAcYCKIpiA+QA37Rx6F+qqs7v6fWEEEKIvpKSryPU2xk77T/3QRVFIcTbifQS6SUnrE9C0/q45rRKPzcHKXYihJUwd2rlTCBNVdVMM59XCCGE6HMpBbqj0iqbhXhJCwJhneKzyrG1UYgKdAPAz9VRUiuFsBLmDuSWAJ+1s2+SoihJiqJ8ryhKpJmvK4QQQlhUVb2BI6W1R7UeaBbq7Ux2WQ16g7QgENYlPquMkQFuONjaAODnZk9+pQRyQlgDswVyiqLYAQuB1W3sjgeCVVUdA7wGrO3gPDcoirJLUZRdRUVF5hqeEEII0SMHC3QA7c7ImVQ4UibplcJ6GIwmdmeXEz3YvWWbv5sj5TUN1DUY+3BkQojOMOeM3DwgXlXVgmN3qKpaqapqVdPXGwFbRVG82zqJqqrLVVWNVVU11sfHx4zDE0IIIbovJb+DQM5bKlcK65Oc39gIPCbYo2WbX1PlSlknJ0T/Z85A7hLaSatUFMVPURSl6esJTdctMeO1hRBCCItKKdDhaGvDYA+n4/Y195JLl0BOWJH4lkInrWfkpAWBENaix1UrARRFcQLOAm5ste0mAFVV3wYuAG5WFMUA1AJLVFVVzXFtIYQQojek5OsI8x2ARqMct8/DyRZXBy0ZJRLICesRn1nGQBd7At0dW7b5NgVy0oJAiP7PLIGcqqo1gNcx295u9fXrwOvmuJYQQgjRF1ILdMyIGNjmPkVRCPUZQEaxrJET1iPhSDnRQe40JU0B/6RWyoycEP2fuatWCiGEECed4qp6iqv0hLVRsbJZqJeTpFYKq1FcVU9mSU1L/7hmzvZaXB20FEggJ0S/J4GcEEIIcQKpTYVOIvxc2z0mxNuZ3IpaqfYnrEJCVjnAUYVOmvm5OciMnBBWQAI5IYQQ4gRSmloPhPkNaPeYUG9nVBWySiW9UvR/8VllaDUKo5oagbfm5+YoveSEsAISyAkhhBAnkJKvw8PJFp8B9u0eE+IllSuF9UjIKiMywLWlEXhr/q4O0n5ACCsggZwQQghxAikFOsL9XI4qCnEs6SUnrIXBaCLpSAXRQcenVUJjamVRVT0NRlMvj0wI0RUSyAkhhBAdMJlUUvN1hHdQ6ATAzdEWT2c7aUEg+r3kfB21DUaiW/WPa83PzQFVhUJdfS+PTAjRFRLICSGEEB3IKa+lWm8kvINCJ81CpHKlsAIJLY3A25+RA8ivkF5yQvRnEsgJIYQQHUhpqlgZ3kGhk2Yh3s7SS070ewlZ5fi42DPIw7HN/f4tgZzMyAnRn0kgJ4QQQnSgpWLlCVIrAUK9nMmvrKNWLy0IRP8Vn1VG9GD3dtd8+rs2Bnh5MiMnRL8mgZwQQgjRgZR8HYHujrg42J7w2JaCJ7JOTnTRvauTWLZur8WvU1JVT0ZJTZv945q5OmpxsNVI5Uoh+jkJ5IQQQogOpDZVrOyMUKlcKbohraiK1XHZfL7jCBU1DRa9Vksj8HbWxwEoioK/myN50ktOiH5NAjkhhBCiHQ1GE2lFVZ1Kq4R/ZuTSZUZOdMGHWzNQFNAbTWzcm2fRayUcab8ReGt+rg4UyIycEP2aBHJCCCFEO9KLq2kwqkR0ckZugL0W7wH2MiMnOq2yroGv4rI5NzqQId7OrE3Isej14jPLGeHviqPd8Y3AW/NzcyBPAjkh+jUJ5IQQQoh2NFes7OyMHMAQb2dpQSA6bfWubKr1Rq6eHMqisYFsTy8lt9wyRUYMRhNJ2eXEtNM/rjU/NwcKKuswmVSLjEUI0XMSyAkhhBDtSMnXYaNRGDrQudOPCfF2Il1aEIhOMJpUVm3NIDbYg1GD3Fg0NgCA9Um5FrleSoGOGr2xw0InzfzdHDCYVEqq9RYZixCi5ySQE0IIIdqRUqAj1NsZe23HaWithXg7U1xVj67OskUrhPX7LbmQrNIarp4SCjT+7owd7G6x9Mr4ThQ6aebn2txLTtIrheivJJATQggh2pGSryO8C2mV0NhLDiCzRGblRMdWbk3H382BOZG+LdsWjw0gOV/XktZrTglZZXgPsGu3EXhrfk1NwaWXnBD9lwRyQgghRBtq9AaySms63XqgWUvlSlknJzqQkq9jy6ESlk4KRmvzz9ux+WMCsNEorE00/6xcQlY50UEe7TYCb605kMuXFgRC9FsSyAkhhBBtSC2oArpW6AQgxEt6yYkTW7k1A3uthkvGBx213XuAPacP82Z9Yq5ZC42UVutJL67uVFolgLezPVqNIqmVQvRjEsgJIYQQbUhtSm3rbOuBZo52Nvi5OkgvOdGu8ho93yQ0thzwcLY7bv/i6AByymvZlVlmtmsmZDWeqzMVKwE0GgVfVwcJ5IToxySQE0IIIdqQnK/DwVbDYE+nLj82xNtJZuREuz7bcYS6BhNXTQlpc//skX442tqYNb0yIascG43CqEEdNwJvTXrJCdG/SSAnhBBCtCG1QEeYrws2mhOvJzpWqLczGVLsRLTBYDTx0d8ZTB7qRYSfa5vHONtrOWukLxv35KE3mMxy3fisMkb4u+Bkp+30Y/zcHGSNnBD9mARyQgghRBtSmgK57gjxcqa0Wk9FrbQgEEf7aX8BuRV1XDU5pMPjFkcHUF7TwB+pRT2+ptGkknSkvNPr45r5N6VWqqo0BReiP5JATgghhDhGabWeIl19l9fHNWuuXCnpleJYK7dkMNjTkZkjfDs87ozhPng625klvTIlX0e13tjlQM7PzYHaBiOVtYYej0EIYX5mC+QURclQFGWPoiiJiqLsamO/oijKq4qiHFIUZbeiKDHmurYQQghhTs09vLo7IxfaHMhJwRPRyt6cCnZklHLlpJATpuza2mg4Z5Q/v+wv6HFz+YQjzYVOuh7IAeRVSi85Ifojc8/Inamq6lhVVWPb2DcPGN70cQPwlpmvLYQQQphFSn4lQJd7yDUL8nRCUaSXnDjaB1sycLKz4cLYwZ06fnF0APUGEz/uK+jRdeMzy/FytmOw54kbgbfm39xLTgqeCNEv9WZq5SLgQ7XRNsBdURT/Xry+EEII0SkpBVW4O9ky0MW+W493sLUhwM1RAjnRoriqng1JuVwwbhBujradekxMkAeDPR1Z18P0yoSssk43Am/Nz60x8JNAToj+yZyBnAr8pChKnKIoN7SxPxA40ur77KZtQgghRL+Skl9JmK9Ll9/4thbq7Sxr5ESLT7dnoTeauPIERU5aUxSFRWMC2XKomEJd94Kpsmo9h4uriQnuXP+41ga62KMoSAsCIfopcwZyU1RVjaExhfJfiqJMPWZ/W38NjyuDpCjKDYqi7FIUZVdRUc8rNQkhhBBdoaoqqQVV3S500izE24n04mqp+CfQG0x8vC2TaWE+DPUZ0KXHLo4OwKTChqS8bl078Ug50PX1cdC4Ts97gL3MyAnRT5ktkFNVNbfpcyHwDTDhmEOygdZJ4YOA3DbOs1xV1VhVVWN9fHzMNTwhhBCiU3LKa6mqN3S70EmzEC9nKusMlNVIC4JT3fd78yjU1XN1Ow3AOzJsoAuRAa7dTq+MzyrDRqMwuguNwFvzl15yQvRbZgnkFEVxVhTFpflrYDaw95jD1gNXNFWvPA2oUFW1e7eXhBBCCAtJLWisWNnTGbnmypWyTk6s2JLBEG9npg7v3g3qxWMD2Z1dweGiqi4/Nj6rjAi/rjUCb82vqZecEKL/MdeMnC+wWVGUJGAH8J2qqj8oinKToig3NR2zETgMHALeBW4x07WFEEIIs0luaj0wvKczctJLTtBYaCTpSDlXTQlBc4KWA+1ZODYARYG1icclMnXIaFJJzOp6I/DW/NwcyKuQ9gNC9Efduz1zDFVVDwNj2tj+dquvVeBf5rieEEIIYSmp+ToC3Bw6XVmwPYM9nNAo0kvO0hqMJt75I40zIwYSGdC99EFL+mBLBi72Ws6LGdTtc/i6OjB5qBfrEnO4a9bwThfhOVjY1Ai8G4VOmvm5OVBZZ6BGb+j2rJ4QwjJ6s/2AEEII0e8l5+sI62FaJYCdVsMgDydJrbSgBqOJOz5P4MWfUrnhw7geN842t4LKOjbuyeOi8YMZYN+zIGjR2EAyS2paipd0Rnxm47HRg7s/Iye95ITovySQE0IIIZo0GE0cLqrudiPwY4V4O8uMnIUYjCbu/CKRjXvyuXRiEHkVtTy+YX9fD+soH2/LxKiqXDkppMfnmhvlh51Ww7oupFfGZ5Xh6WxHsJdTt6/r6yqBnBD9lQRyQgghRJPMkmr0RhPhPVwf1yzUy4mM4hppQWBmBqOJu79M4rvdeTx09giePncUt0wfxuq4bH7al9/XwwOgrsHIp9uzmBnhS1APAqlmrg62zBoxkG9352Iwmjr1mPisMmKC3HvUD9G/qSm49JITov+RQE4IIYRo0lzopKetB5qFeDtTVW+guEpvlvOJxgIe/16dxPqkXO6fF8H1U4cAcPvM4UQGuPLAmj0U6er7eJSwISmXkmo913Sj5UB7Fo0NpLhKz+ZDxSc8trxGz+GiaqJ7UOgEGqtWAtKCQIh+SAK5k5TJpLIzo5SXfkrhSGlNXw9HCCHMzmQy/yxXar4OjQLDBnataXN7WipXSnqlWRhNKveuTmJtYi73zgnnpmlDW/bZaTW8fPFYdPUGHlizu09nQVVV5YMtGYT7ujBpqJfZzjs93AdXB22n0isTmtbSRQd1v9AJgKOdDe5OtpJaKUQ/JIHcSURVVeKzynh8w34mP/sbF779N6/+dogrVuygpKrv704KIYS5bEouZMLTv7I7u/OFHzojOV9HiLczDrY2ZjlfqFdTL7kiCeR6ymRS+c/Xu1mTkMM9Z4XxrzOHHXfMcF8X7psTzi8HCvly15E+GGWjnRll7M+r5KopIT1KazyWvdaGc0b78+O+fGr0hg6PTcgsQ6PAmEE9C+SgcVZOUiuF6H8kkLNyqqqyJ7uCZzYe4PTnNnHem1v5eFsmUYGuvHzxWD68ZgK55bVcs2rXCV/0hRDCGtToDTy8di/FVfU8+e0Bs868pBboetwIvLVBHo5oNQrpMiPXIyaTyv1rdvNVXDZ3zhrObTOHt3vsNVNCmTTEi8c37CerpG8yUj7Yko67ky2Lxwaa/dyLxgZSozfy8/6CDo+Lzyonws8V5x5Wy4TGFgT5ldJLToj+RgI5K6SqKgfyKnnhx2TOfPF3Fry+mfc3pzPcdwAvXjiGnQ/P4r0rx7M4OpCpYT68ekk0e7LLue3ThE4vkBZCiP7qlV8PklNey4XjBrEjo5RfDhSa5bw1egOZpTVmWx8HoLXREOTp1G+bgh8preHuLxMp7sdZGyaTykNr9/DlrmxunzGMO2eFdXi8RqPw4kVj0CgK96xOxGiBFNyOZJfV8OO+fJaMD8LRzjwzu61NCPEkwM2BtQk57R5jNKkkHinvUf+41vzdHMiv6L+/I0JYM1VVWz4AjEYjDQ0NNDQ0oNfrMRjan4iRzo5W5FChjg1JeXy7O5e0omo0Ckwe6s1N04YyJ9IPD2e7Nh83J9KP/y6K4pG1e3lk3V6ePneUWVM9TjTmj7dlcees4bg7tT0+IYTorJR8He//lc5FsYN4+txRxGWV8ez3Bzgz3AetTc/uTR4qrEJVMeuMHDSuk+uPveSKdPUsfX87GSU1TAjxZMmEoL4e0nFUVeWRdXv5bMcRbj1zGHed1XEQ1yzQ3ZHHFkZyz+ok3v3r8FFr6Szto22ZKIrC0knBFjm/RqOwYGwA7/2VTklVPV4D7I875lBhFVX1BmJ6WOikmZ+rI8VV9egNJuy0MgdwqlNVlcrKSvR6PXq9nvr6evR6PR4eHvj6+lJfX8/WrVtbghGDwYDBYCAyMpKIiAgqKytZvXp1y/bmjzPPPJOYmBjy8vJ48803aWhowGQyYTQaMRqNXHbZZYwfP55Dhw7x/PPPH7XPZDJxxx13MH78eOLj43nyySdb9qmqislk4qmnniImJoY//viDxx9/HJPJ1PKhqipvv/02UVFRrF+/nieeeKJle/Pnr7/+mmHDhvHxxx/zzDPPALTsU1WV3377jcDAQN544w1eeumlowI0VVVJTEzEw8ODZ555htdff/245zUjIwNbW1seeeQRVq1a1bJ96dKl7f4sJJDr5zKKq/l2dy7f7s4jOV+HojTejbtqSijzovzwbuMFvC1LTwsmv6KWNzal4evqcMI7muawI72U61btpLLOQL3ByDPnjbb4NYUQJy+TSeXhtXtwcdBy/7wRaG003D83ghs+iuPznUe4/LSevXE2d8XKZiFezvydVoKqqr12E+1EKmobuGLFDgoq6xlgr2VXZlm/C+RUVWXZun18sj2Lm6cP5Z7ZYV16/s6LCeTn/QW89FMq08J8GOHvasHRNqrRG/h8xxHmRPoS6O5osessHhvIO38c5rs9eVzRRo+6+KwygB5XrGzm59b4XqOgso7Bnj1vpSAsLz09HZ1OR01NDTU1NdTW1uLt7c3EiRMBePPNNyktLaW2trblmAkTJnDttdcCMGvWLKqqqo4K1C677DIee+wx6uvrGTly5HHXvPXWW3nggQeorq7m8ssvP27/f/7zHyIiIigvL2fZsmXH7ffw8CAmJgadTsfXX3+NVqtFq9Wi0WjQaDTMmjULgLq6OlJTU7GxsUGj0WBjY4ONjQ16fWN1YEVR0Gg02NratuxTFAU7u8YJBScnJ4YMGYKiKC3HajQanJwaf7d9fX0544wz0Gg0LccoioKzc+Oa5yFDhrB48eKW/c3XbH786NGjufrqq496LICDQ2MF2MmTJ6PVao96bPM4mp/7gQMHtuyPiIjgo48+avPnLIFcP1VR28DDa/eyIamxMtW4YA8eXTCSs0f5tzTn7Kp/zw4nv6Kel385iJ+rg0X/aH+3O4+7vkxkkIcjs0b68tmOI1wwbjDjgs3zR0UIcer5Kj6bnRllPH/+aDybMhDOGunLhBBPXv4llcXRgQzowXqg1Hwd9loNwU0FSswl1NuJ2gYjBZX1+Ll17/XbnGr1Rq5ftYtDhTreu3I8H/2dSXxmWV8P6yiqqvLY+n18tC2TG6cO4b454V0OghVF4enzRjH7f39y1xeJrLt1CvZa86c6trY2IZeK2gaunhJq0euM8Hcl3NeFtQk5bQdymY2NwEPM0L8OwK+pl1y+BHK9pri4mNzcXGpraykvL6eiogIHBwcWLFgAwHPPPUdycjIVFRVUVFRQXl7OqFGjWLlyJQCXXXYZmZmZR53zrLPOagnk3n33XSorK3F0dMTJyQlHR0dCQ//5vR0xYgTQGHzY2tpiZ2fX8lh7e3seffRR7OzssLe3x87ODjs7O4YPb1y76urqyrp167Czs8PW1rYlIPPwaHwPGBgYyJ49e47a1xyUAYSFhbF///52n5uoqCh+//33dvdHR0fz1Vdftbt//PjxjB8/vt39EydObPm3tmXy5MlMnjy53f1nnHEGZ5xxRrv7p02bxrRp09rdP2PGDGbMmNHyfU1N+2t9JZDrh3ZmlHLn54nkV9Zx24xhLJkQZJY7e4qi8Oz5oyiqquehtXvxcbFn5ghfM4z4aO9vTufJ7/YzLsiD966MxdZGw9ZDJY2B6a1Tepz+JITouRq9ARuNgq1Gg0bTP2aJOlJWreeZjQeIDfbggnGDWrYrisKD54xg8RtbWP5HGnfPDu/2NVIKdAz3HYCNmZ+P5hYE6cXVfR7INRhN3PppPDszS3l1STTTwnzYn1vJLwcKKK3WtwTIfUlVVR7/dj+r/s7kutNDuX9eRLdnMj2d7Xj+glFcs3IXL/2cygPzRph5tP9QVZWVW9OJCnQlthduWi6KDuD5H1LIKqk5ruF4fFYZ0YN71gi8Nf+m31tpQdB7brnlFjZv3nzUtrCwsJZALi0tjezsbNzc3AgJCcHNza0l+AJ44oknMBgMODo64uzsjJOTU0sgBRAXF9cSOLXltddea3efoijccMMN7e7XarXExsa2u9/GxgZPT89294vOk0CuHzEYTbz62yFe/+0ggzyc+OqmSWZLi2hma6PhrctiWLJ8G//6NJ7Prj/NbNcwmVSe2niA9zenMzfSj5eXjG0p4f3ogpHc/Ek8q/7O5NrTLXun8mRUWFnHA2v2sGzBSLPPFohTi7GphPtXcdkt27QaBVsbDbY2CnZaTdPXjd/b2mhabVOw09pgZ6PgZKflthnDGCr9YFcAACAASURBVG7mNMT2PPt9Mro6A0+eG3Vc4Dl2sDvzR/vz7l/pXHZacLezFlLydZwx3Mccwz1KiNc/veTM2VOsq0wmlfu+2s2vyYU8uTiKBWMCAIgNafwbEJdZxlkjzX9zrytUVeWp7w7wwZYMrp4SwkPnjOhxMDIjwpdLJgSx/M/DzIzwZUKoZd5AbjlUQmpBFS9eOKZXUmgXjmkM5NYl5hxVxbOipoG0omrOixnUwaO7pvn/lARyHauobeBwUVXjTRtXByYP8+72uW655RYuuugivL29cXNzw83NDXf3f4rXLF++vMPHz5w5s8P9HQVxwnpIINdPHCmt4c4vEonLLOO8mED+uzASFwdbi1zL2V7LiqvGc/5bW7l21S6+vnkyod49Cw7qGozc82US3+3J46rJITwyf+RRd7XnRvkxPdyHl35K4ZxR/n1+V9rafLoji1+TC6ltMPLJdRP7zTobYV0MRhP3rE5iXWIul05snOlvMJqaPlT0BhN6o4kGQ6ttLftNNBhUKmsbaDCayCqtYVdGKWtvncJAF8v+f96VUcoXu45w49QhRPi1vc7pvjkR/Lgvn5d+SuW5C7q+HresWk+hrp5wP/M0Am8twN0ROxtNn1aubJ7l+iYhh3/PDjtqPeGoQDdsbZQ+D+RUVeWZ75N5b3M6V00OYdn8kWZ7rXv4nBFsTSvm7i8T+f6OMyzy93Xl1nS8B9ixYIy/2c/dlkEeTkwI8WRtYg63zhjW8lwlHGlaHzfYPBUrAVwdtDjZ2UgvOaDeYCSrpIbDxdUcLqomvbiq6XM1JdX6luMWjgnoUSA3ffp06urqsLfvXC0EcWqSQK4fWJ+Uy0Nr9gDwypKxLLJA35lj+bjYs+qaCZz/1lauWLGdNTdPwceley8WFTUNXP/RLnakl/LQ2SO47ozQ4/74KorCfxdGMvt/f/LEd/t549IYc/wzTgkmk8rX8dm4OGjZmlbCNwk5Zr3TKk4NBqOJu79MYn1SLvfOCW+zmXJX7M2p4MK3/+b6Vbv4/IZJFimzDo2pgA+v3UuAmwO3d9A7LMjLiSsmhfDBlnSuOT2U8C5WnkwpsEyhEwAbjUKQl1OfVq587bdDrNyawbWnhx73s3ewtSEywI24zNI+Gl1jEPfcDyks//MwS08L5tEF5gvioPEG5ksXjeHCt//miW/38/wFY8x27hq9geV/HubX5EJumzHc4uvwWlsUHcBD3+xlX24lUYFuQGP/OI0CY8wYyCmKckr1kjOZVPIr60gvruZwUVWroK2a7LIaWne08B5gzxAfZ84a6UuotzNDfAYQ6u1MkKwlFL1AArl2HCzQ8e5fh5k5wpcZEQOxtcC6rqp6A8vW7WVNfA4xQe68siS6VxcRh3o7s+Kq8VyyfBtXr9zB5zdM6nKhgOyyGq76YCdZJTW8ekk0C5tSddoS7OXMv84cxks/p3JxbBFTw8yfwnQy2p5eypHSWl66aAwfbcvkye8OcGb4wHbbTQhxLIPRxF1fJrEhKZf75oZzy/SeBXEAUYFuvLJkLDd+HMc9qxN5/ZIYi6y1+2BLOsn5OpYvHXfCxsa3zRjG6l1HeOb7A6y8ekKXrpPaFMi1N+PXUyFezmT0UVPwj/7O4KWfUzk/ZhAPnd12qmJssAcfbsvsk/Lyqqry4k8pvP1HGpdNDOK/CyMtknUwLtiTm6YN5c3f05g1wpfZkX49Op/RpLJ61xFe+jmVQl09Z4/y47ozenfpwDmj/Hls/T7WJuS0BHIJWWWEm6kReGuNveRO/hm5pCPlXPnBDsprGlq2OdraEOrtzOhBbiweG9ASrIX6OONqoewpITpDArk25FXUcsWKHeRV1PHlrmwGuthzYewgLo4NOm5BcXclHinnjs8TOFJaw+0zh3P7jGF9UgRk7GB33rgsmus/jOOWT+J5v6k4SWfsy63g6g92UttgZNU1Ezq19uPGaUP4JiGHZev28sOdU1vW0In2fRWXzQB7LfOi/BkZ4Mr8Vzfz7PfJ3UofE6ee1kHc/fMizNpPa3akHw/OG8FTGw/wolcK982NMNu5AXLLa3n5l4PMGjGwU2+63Z3suHXGMJ7emMyWQ8VM6UJaU3K+DlcHLb6ulkljCvV24s+DRZhMaq8Wl1mXmMOy9fuYNcKX584f1e61xwV78N7mdPblVph9bfaJbNyTzxub0rhkwmCeWHT8GkhzunNWGL+nFPHAmj3EBHt0uoVPa6qq8ntqEc9uTCalQEdMkDtvXR7DuODeL97g7mTHtLCBrE/K5YGzR6AAiVnlLBzb/k3V7vJ1dWBbWonZz9ufqGrjWn+tRsOTi6MY0hSs+bk6yJIG0S/JSsdjVNY1cPUHO9HVGVh/6xSWLx1HVKAbb/2extQXNnH5e9v5dncueoOpW+c3mlTe2HSIC97aSoPBxOc3TOLus8L6tJLjjAhfnlocxZ+pRdz/9Z6WzvId+etgERe/sw0bjcLXN0/u9AJ+e60NTyyKIqOkhrf/SOvp0E96VfUGNu7JY/5ofxztbIjwc+XaM0L5YtcRdqT3XRqUsA4Go4k7v0hkQ1IuD5g5iGt23RmhXDJhMG/+nsbqXUfMeu7/btiHSVV5dEFkpx9zxaQQAt0deXrjAUymE7+WNUvN1xHh52qxN2sh3s7oDSZyK3ovNW1TSiH3fJnEhBBPXr80usO/M82tYeL6oA3BT/vz8R5gx1OL2w80zcVOq+F/F49FV2fggTWd+3vX2r7cCpa+v4OrP9hJncHIW5fF8PXNk/skiGu2ODqAQl092w6XcKioCl29wSLBuL+bAwW6eoxd+H9lbf46WMyO9FJunzmMy08LZvIwb/zdHCWIE/2WBHKt6A0mbvwwjkOFVbx1eQyjB7kzO9KPFVeNZ/N/ZnDXrDDSi6u59dMEJj3zK09vPEBaUVWnz59XUcvl723nhR9TmBPpx/d3TLVY9ayuWjIhiDtnDefr+Gxe/Cmlw2PXxGdz9Qc7GeThyDe3TOnympLTh3uzYEwAb/6e1qeL/63Bxj151DYYuTD2nzVxd8wcziAPRx78Zk+3byj0lfyKOgxG6xqztWowmrjj80S+3Z3Hg2dHcKMFgjhoXDvz+KIopgzz4sFv9rDtsHnu2P96oIAf9xVwx8ywLqWcO9jacN/ccPblVrI2MadTj1FVlZQCHWEWKHTSLLS5cmVx+/2AzCkus5SbP44j3M+Fd6+MPWH2w0BXBwZ7OvZ6IKeqKlvTSpg01LvXZirD/Vy4d044P+8vYHWr6q0dyauo5Z4vk5j/2mb25lbw6IKR/HzXNOaN8u/zN/mzRvgywF7L2oScln6AMUHmWx/XzM/NEaNJpaSq3uzn7g+aU3wD3R1ZMt5yfXaFMCcJ5JqYTCr3fpXE34dLeP6C0ceVoA5wd+SOWcP5874zWXn1eMaHeLJiczoz/+8PLnrnb75JyKauwdju+X/Ym8+8V/4iKbuc5y8YzeuXRuPm1L/yqu+YOZxLJgzmjU1pfPR3xnH7VbVxNvHuL5OYOMSTL2+a1O3qkw+fMwI7Gw3L1u/r8h3RU8lXcdkM8XYmptXdVSc7LU8sjuJQYRXL/+z/s5qHi6p4/beDzHvlL0575lf+8/Wevh5Sn0orquL7PXkdvl70VGMQl8B3e/J46OwR3DDVMkFcM1sbDW9eNo4gTydu/CiOw124wdWWWr2RR9fvY/jAAd1qV7JgdABRga68+GNKp57nvIo6dHUGwi20Pg4g1Kepl1wvrJM7kFfJ1R/sxN/NkVXXTOj0Gp5xQR7syizr1dfktKIqinT1TOnltgzXnh7KxFBPHt+wnyOl7QfXuroGXvgxmekv/M6G3bncMHUIf9x7JldPCe31tYTtcbC1YU6kHz/szefvwyV4ONn2uBJ1W/yaWhCcrJUrf9xXwO7sCu6YNbzf/GyFOBH5TW3y3I/JrEtsrObWUUVAG43C9PCBvL10HFsfmMF9c8MpqKzjri+SmPj0rzy2fh/J+ZUtx9fqjTz4zR5u+jiOwR5OfHvb6VwUO7jP7+C1RVEUnlgUxcyIgSxbv48f9ua37DM0VY574ccUFo8N4IOrOv/moC2+rg7cMzuMP1OL2Lgn/8QPOAVlllSzI72U88cNOu735czwgZwzyp/XfjvUL2c104qqeO3Xg8x9+U9m/N8fvPhTKo62GuZG+vF1fDbrk3L7eoi96lBhFa82PR8z/+8Pbv4knukv/M5nO7LMPkPZYDRx+2cJbNyTz8PnjOD6qUPMev72uDnasuKq8dhoFK5dtYvyGv2JH9SO1347SHZZLU8ujurWGyqNRuHBs0eQW1HHyq0ZJzw+Jb+x0Em4BXvi+bo44GBr+RYEmSXVXLFiB052Wj66dkKX1oCNC/GkSFdPdlnvpX9uOdQ4gzt5aPfLtHeHRqPwfxc1Vq68Z3XScemCDUYTH/2dwfQXfueNTWnMi/Ljt3um8cC8Ebg59q+bsNCYXqmrN7AhKZfoIA+LvMdobgp+MgZyRpPKSz+nMMTHmfOiLV85XAhz6dfFTipqG058kBms2prBO38c5rKJQdwyvfN3rge6OHDL9GHcNHUo2w6X8NnOI3y6PYuVWzOIDnJnwegAPtmeSVpRNTdOG8I9Z4X3+7s8WhsNr10azaXvbueOzxP45LqJjAxw5fbPEvjlQCE3Tx/KvbPDzZICs/S0YFbvyubxb/cxLdynyxUzT3Zfx2WjUeC8mLb/qCxbMJI/U4t4ZN1ePrxmQp/fHDhUqOO73fls3JPXUsp9XLAHj8wfybwoPwLcHTEYTVz0zt889M0eYoLcGeRx8pZnPvb5UJTGyoCPLhhJsJcTr/92iAfW7GH5n4e5Z3YYZ0f59/j/VXMQ9/3exiDuujN6J4hrFuzlzPKl47j03e3c+FEcH107scuveQcLdCz/8zDnxwxi4pDuz9JMHurNjIiBvLHpEBfHDu6wymvz76slAzmNRmmsXGnBQK6wso6l7++gwWji0xsndfn/17igf9bJ9VYF5a1pxQzycDRbIbGuGOThxGMLI/n36iTe33yYG6YORVVVft5fwLPfJ3O4uJqJoZ58cM4IRg8yf6qiOU0e6o2Piz1Funqz9o9rrTkDJ78X13n2lg1JuaQWVJ1wLakQ/U2/fuecVVrD4xv288DZERYp/w+NKY+PbWis6PX4oqhuvRnWaBQmD/Nm8jBvSqv1rInP5vOdR3j82/0MdLHn42sncvrw3r3b2BNOdlrevzKWC97+m2tX7SLYy4m9ORU8sSiSpZNCzHYdrY2Gp86N4ry3tvK/n1N5ZP5Is53b2jX2jsvh9OE++Ls5tnmMr6sD980N55F1+1iflNsr/QePdbBAx3d78ti4J4/UgqqjgpV5Ucc3ftfaaHhlSTTzXvmLu75I5LPrTzup/mi29XyMD/bksQUjmTfKH1/Xf56PM8MH8suBQl78MYVbP00gMiCN++ZGMHW4d7dehxqMJm77NIEf9uXzyPyR3UpJNIfYEE+ev2A0d36RyIPf7OGFC0Z3+t+jqioPr92Ls72WB8/ueQXMB+ZFMOflP3n1t4MdFkxJydfh5+pg8XT3EC9nUgt1Fjl3RU0DV6zYQXFVPZ9efxrDuxGUhvu5MMBey67MUhb3wqyE0aTyd1oJc6N61gagJ86PCeTn/fm8+GMqA10c+HRHFjvSSxnq48x7V8Qyc8TAPr9J1hk2GoUFowNYsSWdmGDLVB31dLLDzkZDfuXJtUauwWjipZ9TGeHvytlRvdPMXQhz6XEgpyjKYOBDwA8wActVVX3lmGOmA+uA9KZNa1RVffxE5/YaYMeKLenszi7n9Utjur0eqz1xmaXc8XkCYwa589ol0diYYZbJ09mO684YwrWnh5KcryPA3bFfpmGciNcAe1ZdPYHz3tpKaoGOty8f1+OeO22JDvJgyfggVm7N4PyYQYwMsNwaFWvy9+EScspr+c+8jt/MXjoxmK/ic3ji2/1MDxvYK+suUwt0fLe7MVg5WNhxsNKWwZ5OPLE4kru+SOLN39M6bPLc36mqSmpBVUvwdqj5+Qjx5L8LI5kb5dfu86EoCmeNbOxTuS4xh5d+TuXKFTs4bYgn982NOGpd5InoDSZu+yyeH/cVsGz+SK7poyCu2eLoQA4XV/PqrwcZ4uPc6b51a+Jz2J5eyjPnjcKrG2XhjzXc14WLxw/m422ZXDU5hGCvttcNpeTrutxAvDtCvJ35NbkAg9Fk1hsYtXoj16zayeGialZcNZ6x3ZyRsdEoRAe5E5dZbraxdWR/biWVdYYutYkwN0VRePrcUcx5+S/u/CIR7wF2PLk4iiXjB1vdTabrzgjFpKrEhlgmkNNoFAa62p90M3Krd2WTVVrDiqtie7U1iBDmYI4ZOQNwj6qq8YqiuABxiqL8rKrq/mOO+0tV1fldOXGAmyOPXhLNf77ezfzX/uLVS6LNlkefVlTFtat24e/mwPtXxuJoZ95+ZoqiMMLfuoOSIC8n1t86hXqDySILp5v9Z244P+3L5+G1e/jqpsnyQgqs3nUEFwcts0f6dnicjUbh6XOjWPj6Fp79IZlnzhtlsTH9dbCI/27Y36VgpT3nRg/i95QiXvn1IFOGebeUPrcGzRUON+7O47s9eaQVVaMoMCHEkysWRTI30o+BXXg+bDQK58UMYv7oAD7bkcVrvx3ivDe3ctZIX/49O/yEAYbeYOJfn8bz8/4CHl0wkqun9G0Q1+yuWcPJKK7m+R9SCPFy5uxRHd/pLq/R89TGA8QEuXNx7GAzjiOMdYm5PP9DCm9cFnPcfoPRxKGiql7Jmgj1dqLBqJJbXme2VEK9wcTNn8SRkFXGG5fG9PjfERPkwWu/HURX14CLhRsdb00rBmBSD1JozcFrgD3LrxjHzvRSLjst2GrT/APcHXlsYedbdXSHv5vDSbVGrq7ByKu/HiQmyJ0zwwf29XCE6LIe325SVTVPVdX4pq91wAHAbDkZC8YEsP7WKbg52nL5e9t5+4+0HlfUKtTVceWKHdgoCquumWCWO78nqwB3R4sGcdDY0PSBs0cQn1XOl2buQ2WNKusa+GFfPgvHBHSqYXpkgBvXnh7KZzuy2JVhmd5yn+/I4qoPdgLw+KJItj8wky9vnMSVk0O6HMQ1e2JxFP5uDtz5RQK6ut5ZD9tTDUYT167axdyX/+L1TYfwcbHniUWRbH9wJl/cOIkrJoV0KYhrzU6r4crJIfxx73T+PTuMbWklzH3lT+7+IrHdqnqtg7jH+lEQB403s56/YDQxQe7c9UUiiUc6nuV57ocUKmobeNLMvcQGujpw/RlD+G5PHvFZx5fWzyipQW8wdbmNSneEeJm/cuWT3+3n95Qinjp3FPNOECx3RmyIByaVE/68zGFLWgnDBg7o9v8Zc4oJ8uDGaUOtNojrLX5ujuRXnjyB3MfbMsmvrOPfc8KtIoVWiGOZNW9AUZQQIBrY3sbuSYqiJCmK8r2iKF26ZTRsoAvrbj2deaP8efb7ZG78KI7Kbr7xq643cM3KnZRU6Vlx1fh2U21E7zo/JpAJIZ48+0MypdXdr3Z3Mvhudx51DSYu7MKsxJ2zhhPo7shD3+ylwYxVEFVV5YUfk7l/zR6mDPPmm1sm9yhYac3VwZaXLx5LTlktj67bZ4bRWpaqqixbt4/fkgu5+6wwtj84i89vmMTSSSEMdDHfG1Fney23zhjOX/85kxumNgYgM/7vdx5bv48i3T9rU/QGE7d80hjE/XdhJFf1oyCumYOtDcuviMXHxZ7rVu0ip7ztlKy4zDI+25HF1ZNDLJJefcPUIfi42PP0dweOuxGY2lToJKIXUiubb4qZq+DJvtwKPmpKG71kgnn6Xo0d7I6iWL4xuN5gYmd6aa+3HRA94+/mQH5F3UnRNqiq3sBbv6dx+jDvXq+aKoS5mC2QUxRlAPA1cKeqqpXH7I4HglVVHQO8Bqzt4Dw3KIqyS1GUXUVFRS3bB9href2SaJbNH8lvyYUsfG0zB/KOvUzHGoyNb3wO5Ol447JoxliospPoOkVReGJxFFV1Bp79/kBfD6dPfRWXzbCBAxgzyK3Tj3Gy0/L4okhSCnS8+9dhs4yj3mDkzi8SeWNTGpdMGMz7V8aaPdUqNsST22YMZ01CDus62by5r3ywJYPPdmRx07Sh3D5zOD4ulp3Jd3ey44F5I/jj3jO5YNxgPtqWybQXNvF/P6VQUlXPLZ/E8cuBAh5fFMmVk0MsOpae8B5gz4qrxlPfYOTalTuPm31tbm3i5+rAnWeFWWQMzvZa7poVxq7MMn7cV3DUvuR8HRoFhg20XDPwZj4u9jjb2ZBuhkBOVVUe37AfDyc77pplvufNxcGWcF8XiwdyiUfKqW0wMkneQFsVX1cH6g0mymusI4uiIx9sTqekWs+/54T39VCE6DazBHKKotjSGMR9oqrqmmP3q6paqapqVdPXGwFbRVHafPVWVXW5qqqxqqrG+vgc3ZRbURSuOT2Uz284jdoGI+e+uYU18dmdGqOqqjz0zR7+SC3iycVRzIjoeO2R6H3hfi5ce3ooX+7KtliKYH93uKiKuMwyLmyjd9yJzBzhy9xIP1755SBZJe03uO2M8ho9S9/bwbrEXO6bG87T546yWOXY22YMY1ywBw9/s7fDxrx9aVNKIU9+t5/ZI325r5f/6Pu5OfDMeaP45e5pzBzhy2u/HWLC07/yy4FCnlgcxRVmrCRrKWG+LrxxWQwHC6u4/bOEo3rnrdyawYG8Sh5bONKiaW0XxQ5i2MABPPdD8lGz1qn5OkK8nDuVxtxTiqIQ7OVMhhlSK3/Ym8/29FLuPivM7EWOYkM8SMgqP663mjltTStGUfp+fZzompOll1x5jZ7lfx3mrJG+3S4OJER/0ON3Zkrju833gQOqqr7UzjF+TcehKMqEpuuWdPeasSGefHvbGYwd7M7dXybx0Dd7qDcYO3zMy78c5Mtd2dw+Y5jZUlCE+d0+czgBbg48vNY8KYJ1DUY27slj5ZZ0q0gF+aqpd9y53Sz9/djCSGxtNDy8bm+3/71ZJTWc99ZWEo+U8+ol0dwyfZhF1w5obTS8fPFYAO78ItHsDbJ7KrVAx22fJhDu58r/Lh7bZ8V4Qr2dee2SaL697XTmRvnx/AWjWXpacJ+MpTumhvnw2MJINqUU8eR3jbPueRW1/O/nVM4M92GOBaritqa10XD/3AjSi6v5bEdWy/aUAl2vrI9rFurt3OMZuboGI09tPECEnwtLxpuvMEyzccEeVNUbWhqlW8LWtBKiAtx6pdKuMJ/m6uEFVr5O7p0/D1NVb+Ce2ZbJAhCit5jjFvsUYCkwQ1GUxKaPsxVFuUlRlJuajrkA2KsoShLwKrBE7eG7ap+m/mw3TRvKJ9uzuPDtv8kua/tu/hc7s3jl14NcOG4Qd1kodUeYh7O9lmULIknO17Fqa0a3ztFgNLEpuZC7vkhk3BM/c8sn8Ty2YT+fbM868YP7kNGksiY+h2lhPt1eg+bn5sC/Z4fxZ2oR3+7O6/LjE7LKOPfNLZRW6/n4uoksHBPQrXF01WBPJ548N4q4zDJe33SoV67ZGSVV9Vy7aieOdja8f2Uszv2gEEJUoBtvXBrDRWas7Nhblp4WzDVTQlm5NYMP/87g8Q37MZhU/ruwez08u2rmiIFMDPXklV8aqzLWNRjJKKnuldYDzUK8ncguq+3Rjar3N6eTXVbLsvkjLVIiPzbYE4C4NorDmEON3kBCVhmTh8lsnLXxc7X+GblCXR0rt2SwYHQAEX7WXV1cCHNUrdysqqqiqupoVVXHNn1sVFX1bVVV32465nVVVSNVVR2jquppqqpu7fnQm+6wzovgnaXjSC+qZv5rm/kjteioYzYlF/LgN3uZGubD0+eNkqpEVmBOZGNvrf/9nEpeJ/vVmEwq2w6X8OA3e5jw1C9cvXInvx4oYP7oAD65biJTw3x48rv9pBVVWXj03bf5UDH5lXVdKnLSlqWTQhg9yI3Hv91PRW3n1zH8sDefJcu34Wyv5eubJzMh1LNH4+iqRWMDOTc6kFd/PUhcZt+n1tYbjNz0cRyFlfW8e0UsAe5tN2YXXfPQOSOYGTGQR9fv4/u9+dw+c7jZSvGfiKIoPHTOCEqq9bzzx2EOFlShqvRuIOfljNGkdjuNuKCyjjc2HWJOpC+TLdR/bZCHIz4u9sRbaJ3crowyGoyqFJiwQj4u9mgUrLqX3Jub0tAbTXJjX5wUrKvbZTvmRPqx4bbT8XN14KoPdvDyL6mYTCq7s8u55ZN4Rvi78OZlMRZb4yPMS1EUHlsQicHUuJi/Para+DN+8tv9TH72N5Ys38Y38TmcPtyHd6+IZefDs3jugtFMGebNCxeMxsHWhru+SDRrVUdz+iouG3cnW2aO6Fkvm8becqMoqarnhR+TO/WY9zenc/MncYzwd2XNLZMZ6mP5wg9teXxRJIEejtzxeWK3K9Oag6qqPLhmLzszynjxwjGyhsKMbDQKr14STWSAKxF+Llx/xpBevf7oQe4sHBPAe5sP8+fBxht/vRnIDfFpqlzZzXVyz/2QjMGo8tDZI805rKMoisK4IA92WeiGypa0YmxtFMZbqHG1sBxbGw0+LvZWOyOXU17Lp9uzuHDcIIu3VhKiN/R9npCZhHg7880tU3ho7R5e/uUgcZllHMirxGuAHSuuGi+9YaxMkJcTt80Yxos/pbIppfCoRp2HCnWsT8xlfVIuGSU12NooTAvz4YGzIzhrpC9Odsf/rH1dHXjm3FHc/Ek8r/16kLtn968qVRU1Dfy4L59Lxg/GXtvzogtRgW5cPSWUFVvSOS9mEDFBbb9hMppUnvh2Pyu3ZjA30o+Xl4ztlaIP7XFxsOXli6O56J2/WbZ2Ly8vie6Tcbzz52G+js/mjpnDWdBL6aWnEmd7LWtvmYLBpGKn7f0bbPfOCeeHvfm88stB7LQagj17Z0YQWvWSK+76jFzikXLWxOdw07ShFp/Fh/FOfgAAIABJREFUjA3x4Id9+RRW1pm9z9vfaSVED/Zo87Va9H/W3Evu1V8OAnDbzOF9PBIhzOOkmqJytLPh/y4cw9PnjmL74VIMJpVV10wwa48n0XuunzqEIT7OPLpuH2lFVbz1exrzXvmLWS/9yWubDhHg7siz541i50OzeO/K8SwaG9jhG4N5o/w5P2YQr286ZPHS2l21YXcueoOJC8aZb93T3WeF4efqwINr9rQ5C1mjN3DjR3Gs3JrBdaeH8uZlMX0axDUbF+zB7TOGszYxl7UJvd+S4Md9+Tz3QzLzR/tz5yz5Y28pWhtNn/2+DfZ04srJweiNJoYPHGCRdWbt8XS2w8VB2+Veco3tBvbh42LPrTOGWWh0/4gJbrz5Y+7XyoqaBvbmVDBJ+sdZLT9Xe/KtcEbucFEVX8Vnc9lpQQRKqrw4SZxUgRw0poRcOjGI7+88g3X/mtJnKWKi5+y1Njy5KIqs0hpm/t8fPPdDMvZaDcvmj2T7AzP59PrTWDIhCHcnu06f87GFIwlwd+TuLxOprjdYcPRdszoumwg/F6ICzbfw2tley38XNhaOWbE5/ah9hbo6lizfxm/Jjc2kH54/ss+qMbblX2cOJTbYg4fX9m5Lgn25Fdz5eSKjA9148cIxsqb2JHbrmcPxcLJldBf6NZqDoiiEene9BcH6pFzis8q5d054r2SYRAW4YafVmD2Q25ZegkmFKRZa3ycsz9/N0SoDuf/9chA7Gw23TLf8jRAhestJF8g1G+ozgGAvyX+2dpOHebNs/kjunRPOn/eeydp/TeGa00O7nerj4mDLSxeNJau0hie+bX/9XW86WKAj6Ug5F3Sjd9yJzI70Y/ZIX/73S2pLQHSoUMd5b27lYEEVy5fG9stm0lobDf+7eCwKcMfnCb3SkqBQV8d1q3bh7mTLu1fE9ovZSWE5bk62fH/HVB48e0SvXzvEq2stCGr0Bp7ZmMyoQDcuiBlkwZH9w06rYcwgN7NXrvw7rQQHW42sO7Vifm4O6OoNVPWjm6Ensj+3kg1JuVw9JQQfF/u+Ho4QZnPSBnLi5HHN6aH868xhZlsTMiHUk5umDeXznUf4eX+BWc7ZE1/FZ6PVKCzuZu+4E3lsYSQ2isKydXvZmlbMeW9upa7BxBc3nsaskb4WuaY5NLckiM8q57XfLNuSoK7ByPUfxlFe08C7V8SafU2Q6J/83Bxwcej9PmYh3s7klteesP9ps7f/OEx+ZR3LFvTuzHlMsAd7cyqoa+jcODtjy6Fixod49snaSGEezU3BrWlW7qWfU3Bx0HLj1KF9PRQhzEpeScUp6a5ZYUQGuHL/17sp0tX32TgMRhNr4nOYHj4Q7wGWuUsY4O7I3bPD2ZRSxOXvbWegq8P/t3fv0XGd5b3Hf8/MaDS6ayTZkmxpZDt3h8SO5CQkKdA0QENoSKGkSddZPSnQk1Ogqy1tTxtKS1taTgulLfTQdUpa6KGrbbgEaEJIITcoECcktrFzs53YJtbIsi3Jumt0G817/pgZRbZ1156b5vtZa5a29t6z9zv2q5Geed/3efSND1yvK1vy/xPx23Zu1ruu2qz/88Qr2vNqZjLoOef0v+5/Tgeig/r0nTv1us3ZnWqH4rO1oVwJp2VNGz4xOK7P/ddR3bpjk67ekt2SILva6jQ94/T8iSFPrtczMqFXekaZVlngGqsLK5Db1zmgxw726H++cRsF6LHuEMihKAUDPn36jp0anYzr97/2nNZYn37VfvBKn3pHJnX7rsxOl7rrujZdu7VOb7hog772/uvVmsUsfWv1p7ddrpZwecZKEvzd40f0zQPd+r2bL9HPXt7k+fWBc60kc+VfPHxQZtI9b7s00806T3sk+WHPnle9mV751NEzkqTrSXRS0NIjcsut85prn/rOYdVXBPWeG7bmuimA5wjkULQuaqzSPW+7VE8c6tG/P9OZkzZ8dW9UdRXBs8orZELA79OX7n69vvjea1RTVlifSFaFSvSZO3fq1PCE/vAbL3gadD/0XLf+9rGX9a72zXr/m5hyg+xI169aKnPls6/266HnTuruN16Qkyx79ZWl2tZQ4VnCk91Hzqg6FNDlmxj1LmTpEbnTBVCCYPeRPu0+ekYfuPFCVVCGCusQgRyK2l3XbdEbLmrQnz90UMd6R7N674GxKT32Uo9+fufmrKwXKeQMjFdFwvqtmy7Sgwe69a8/6vRkzc6B6KB+5ysHtKstrL941xUF/e+DwlJbHlRteYmOLRLIJRJOf/rNF9VUHdKvvSm7RdPnam8La1/ngCcfoOw+1qfXb6uXP48y5GLlQiV+hctL8r4ouHNOf/XIYTXXhPTfro3kujlARhDIoaj5fKa/evcOBQM+fegrB+att5YpDx7o1tRMQu/uyE4WukL3gRsv1DVb6vRH//GCtn/02/qZT31P7//XvfrbR1/Wfz5/Ukd7RzWTWN4fmyeHxvU//mWPNlSV6nO/3OFJEXZgJbbUVyw6Inf/3i69cGJYH77l0pwWzu5oC6t/bGpFWTbnE+2PKdo/zvq4daKpAEoQPH6wRz/uHNRv3HQRWYixbjHOjKLXVBPS/37nFfrgv+/TZ584og+95eKs3Pf+vV3a3lyt7Zu8qx23nvl9pv/33qv13UO9Onx6RIdPDevQqRF9+8VTSg8WlAZ8unBjpS5pqtIljVXJr01VaqoOzY64xabi+tUv7lFsakb/+qvXqj5DSWaAxWxtqNDTx87Me2xkYlqf/M5htUdq9Y4dm7LcsrPtmlMYfNsa6rLuPtonifVx60VzTSivR+QSCadPPXJYbfXlfFiKdY1ADpD09iub9fjBzfrsd4/opy/ZoKsi4Yze79CpYT1/Ykh/fOv2jN5nvSkPBvT2K5v1djXP7hufmtGRntHZ4O7w6VE9eaRPX993Yvac6lBgNqh7tS+mgyeH9flfuVoXN1bl4mUA2lJfoW/8+ITGp2ZUFjx7tODvv3tUfaOT+vxdu3I+5feCDZWqDgW09/iAbt/Vuurr7D56RhuqSnXhxtUHg8gfTTUhHYgO5roZC3ro+ZM6dGpEn75jp0r8TD7D+kUgB6T8yW2X60c/6deHvrxf3/qNN2R0YfT9e7pU4jfdtjMzteOKSVnQrytaanRFy9kJFAZjUzp8akQvnx7RodTXB/Z3a2wyrj/6ue0ZTzADLGbrhmTCk+P9Y7q06bVR+eNnxvSFH/5Ev9Deoh15UDTb5zO1t4XXlPDEOafdR8/o+gvqcx6YwhtN1SGdGZvSxPRM3k1bjM8k9OlHX9YljVW6Nccj2kCmEcgBKdWhEv3NL+7Qnf/4tP78Wwf1F++6IiP3mZ5J6D/2n9BNlzaqriKYkXsgmVDi2m31unbba1O5nHMan57J6ZojQJK21r+WuXJuIPfxbx1UwG/6vZsvyVXTzrOrLazvHe7VYGxKteUrf8860jOq3pFJplWuI02pEgQ9w5OK1OdXOZuv7zuhY31j+twvd5BYB+se483AHNduq9fdb9ym+57p1GMvnc7IPb53uFd9o1PM288BMyOIQ17Y0pD843duLbndR/r0yEun9cEbL5xN8Z4P2lPr5H7cubqpdE8eSa+PI9HJepHPteQeOHBCF22s1Fu3N+a6KUDGEcgB5/jtt1ysy5qrdc/Xn1Pf6KTn179/b1QNlaV60yUbPL82gMJQFSpRQ2VwNnNlfCahjz30klrCZXrfT+VX4eKdrbXy+2zV0yt3Hz2j1roytdbl18gNVi8dyJ3Kw1pyx8/EdFlzNdN4URQI5IBzlAb8+sydOzU8Edc9X3vO0wLUZ0Yn9fjBHr3zqk0swAaK3Jb6Cv3kTDKQu+/ZqA6dGtFHbrks79YclQcD2t5crT3H+1f83JmE09PHzugGRuPWlfSIcb6VIIjPJHRyaEIRPjRAkeAvSWAeFzdW6fdvvlSPHezRfc9EPbvuA/u7FU84vbtj9dnfAKwPWxqSteSGYtP6m0cO69qtdbr5dU25bta8OtrCOhAdWnGtzRe7hzQ8Edd1rI9bV6pCJaosDeRdCYKTQxOaSTi11pXluilAVhDIAQt4z/VbdMOF9fqzh15aczHctK/u7dKVLTW6pIm090Cx29pQoZ6RSf3ltw9qcHxaH711e95OB+toC2t8ekYHTw6v6Hm7jyZr5RHIrT9NNSGdzrOplZ39yTWnrWFG5FAcCOSABfh8pk/dvkMlftOHvrx/xZ9En+vF7iEdPDms20lyAkDJqZWSdN8zUd15dUSXb6pZ4hm50zGnMPhKPHmkTxc3VmpjVf4kb4E3mqrzryh4NB3IMbUSRYJADlhEc02ZPv7OK7Q/Oqgr/+QR3fbZH+r37j+gf/rBMf3wlT71jEwsew3dV/d0Kej3UdcGgKTXMldWlQb0O2+9OMetWdym2jJtqgmtKJCbiif07Kv9ZKtcp5pqQnm3Ri46EJPfZ7PJWID1jjzcwBJu3bFJwYBPPzrWr8Onh/XEoV59ZU/X7PG6iqAubqzUpU3VurixSpc0VenixkpVhUpmz5mKJ/TA/hN6y+WNq6rDBGD92dZQqYbKUv3mTReqobI0181Z0koLg++PDmpiOkH9uHWquSaknpEJxWcSCuRJ8q7O/nFtqg3lTXuATCOQA5bhZy9v0s9e/loSgjOjkzp8ekSHT43o5dMjOnRqRF/dE9XY1MzsOZtry1JBXZWccxqITVM7DsCssqBfz/zBTfIVSNHijrawHnrupLoHx7WpdulkEk8e6ZPPkvU5sf401YSUcFLf6NRsgfBci/bHWB+HouJJIGdmN0v6jCS/pH9yzv3lOcdLJf2LpA5JZyTd4Zx71Yt7A7lQX1mq6ytLz5oy5JxT18D4bGD3cirQ+8ErvZqecWquCemNF1E7DsBrCiWIk6RdbXWSpD3HB/SOZQRyTx09o9dtrlFNWcmS56LwNFW/VhQ8XwK5roGY3nwZhcBRPNYcyJmZX9LfS3qLpC5Jz5rZg865l+ac9j5JA865C83sTkmfkHTHWu8N5BMzU2tduVrrynXTnF8k0zMJ/aRvTFWhgPwF9EcbAMx1aXOVykr82nd8QO9YYq1vbCquH0cH9L6f2pal1iHb0sFbvqyTG5uMq290ikQnKCpeTCK+RtIR59wx59yUpC9Juu2cc26T9MXU9v2SbrJ8zbEMeKzE79PFjVVqrqGuDYDCVeL3aWdr7bLWyT376oCmZxzr49ax9O+0fMlc2TUwLklqCfO7FsXDi0Bus6S5FZO7UvvmPcc5F5c0JGned3czu9vM9pjZnt7eXg+aBwAAvNDRFtZLJ4c1Nhlf9LzdR/tU4jddvaUuSy1DtoXLSxQM+PKmlly69ECEETkUES8CuflG1s7Nx76cc5I7nbvXObfLObdrwwbWEwEAkC862sKaSTgd6Bpc9LzdR87oqkhYZUF/llqGbDOzvKolFx2ghhyKjxeBXJek1jnft0jqXugcMwtIqpHU78G9AQBAlrRHUoXBX114euVQbFovdA/pBurHrXv5VEuusz+mshK/6iso8YPi4UUg96yki8xsq5kFJd0p6cFzznlQ0l2p7XdLesItt4oyAADICzXlJbpoY6X2di4cyD117Iyck66/kPVx611zTUin8mZq5bgideUiBQOKyZoDudSat1+X9B1JByV9xTn3opl9zMzekTrt85LqzeyIpN+WdM9a7wsAALKvoy2sfccHlEjM/3nsU0f7VFbi146W2iy3DNnWVJ0ckcuHz+a7BmJqrSPRCYqLJ3XknHMPS3r4nH0fnbM9Iel2L+4FAAByp6MtrC89G9WR3lFd3Fh13vHdR8/omq11Cga8mPSDfNZUE9LUTEL9Y1OqryzNWTucc+rsj+n1FJ9HkeFdFgAALFtHW2qd3DxlCHqGJ/RKzyhlB4pEc026KHhup1f2j00pNjVDxkoUHQI5AACwbFsbKlRXEdSeeRKePHXsjCTphgtJdFIMmlK15HJdgiCaqiFHxkoUGwI5AACwbGam9khY++ZJePLkkT7VlJXosubqHLQM2dZUnR8jcp396dIDrJFDcSGQAwAAK9LRFtZP+sZ0ZnTyrP27j57Rddvq5feRObAYbKgqld9nOS9BkC4G3hpmRA7FhUAOAACsyK4t56+Ti/bH1DUwTtmBIuL3mTZWleZ8RK5rIKb6iqAqSj3J4QcUDAI5AACwIldsrlGJ386qJ/fkkT5JItFJkWmqCeV+jVz/uFpYH4ciRCAHAABWJFTi1+s212jvnIQnu4+e0caqUl2woTKHLUO2NVWHdHJoPKdt6OyPqTXM+jgUHwI5AACwYh2RsJ47MaTJ+Iycc9p99Iyuv6BeZqyPKyZNNSGdzGFR8JmEU/fgOKUHUJQI5AAAwIp1tIU1FU/oxe5hvdIzqr7RSV1P2YGi01wTUmxqRiOT8Zzc/+TQuOIJR+kBFCVWhQIAgBWbLQz+6oBK/MlRONbHFZ/GVAmC00MTqg6VZP3+nes4Y2WuRjlROAjkAADAim2sDqm1rkx7jw9oxjlF6srVsg7/mMbimlNFwU8OTeiixqqs37+rP7k+b71NraysrNTw8LDMTKWlpbluDvIUUysBAMCq7Gqr057j/Xr62BndQNmBotRckxyRy1UtuehATD6TmmtDObl/plRVVWnLli3y+XwaGRlRIpHIdZOQhwjkAADAqrS3hdU3OqWRibiuu4D1ccVoY3VytOhUjkoQdPbH1FxTphL/+vuTtrS0VJFIRBs3btTY2JgmJnJb5gH5Z/31egAAkBUdkfDs9nXbGJErRqUBv+orgjkrCh7tj627aZVzmZnC4bC2bNkiv9/P6BzOQiAHAABW5ZKmKlWWBnRJY5U2VLGOp1g11YR0Kke15KID42qtW/815Bidw3xIdgIAAFbF7zP97lsv1oaq9bU+CSvTXBNS10D2A7nxqRn1jkyuy4yV80mPzlVUVOjUqVMaGRlRRUWFfD7GZYoVgRwAAFi1X7lha66bgBxrrA5p7/GBrN+3ayBVemAdT62cTzAYVGtrqwYHB9Xb26tAIKBQiA9TihEhPAAAAFatuSakgdi0JqZnsnrfaJEGctLZa+dKSkpYO1ekCOQAAACwak2pWnLZLkEQTdWQK4Y1cgsJBoNqaWlRU1OTxsbGND6em7WKyA0COQAAAKxaupZctjNXdvbHFCrxaUNlcSfaMTPV1NRo69atCgaDjM4VEQI5AAAArFpjdTKQO53lWnLR/phaw+Uys6zeN1/NHZ2LxWKKxWKKx+O5bhYyiGQnAAAAWLWmHI3IJUsPFN/6uMWkR+fKy8s1NDSksbExjY6OyjknM1NJSYmCwSDB7zpBIAcAAIBVqywNqCoUyGotOeecov0xXbMlvPTJRaikpEQNDQ1qaGhQIpHQ1NSUJicnNTo6qlgspkQiITOT3+9XMBiU3+/PdZOxCmsK5MzsryTdKmlK0lFJ73HODc5z3quSRiTNSIo753at5b4AAADIH801IZ3K4tTKwdi0RifjjMgtg8/nUygUUigUUk1NjZxzmp6e1tTUlGKxmEZHR2eTpPh8PpWUlCgQCDBqVwDWOiL3qKQPO+fiZvYJSR+W9PsLnHujc65vjfcDAABAnmmsDmU1a2Uxlx5YKzNTMBhUMBhUZWWlNm7cqHg8runpaU1MTGh0dFRjY2PnPc/n882O4vl8vtkHcmdNgZxz7pE53z4t6d1raw4AAAAKTXNNSIdPjWTtfp39qUAuTCDnhUAgoEAgoLKyMoXDYSUSCU1PTyuRSCiRSGhmZmY22IvH44rH45qamlo0mYqZnRXomdnsKN9yt7E4L9fIvVfSlxc45iQ9YmZO0uecc/d6eF8AAADkUFNNmXpHJzU9k1CJP/OjNNSQyyyfz6fS0uWVdUgHe3ODvkQioXg8Prtv7nlzt51zcs5pZmZGzrnzzluLdIKXTEq3Oc3L+6Wvnfrq5jtnyUDOzB6T1DTPoY845x5InfMRSXFJ/7bAZW5wznWb2UZJj5rZIefc9xe4392S7pakSCSyVPMAAACQY03VITkn9Y5MalNt5oOr6EBM4fISVYVKMn4vLC6TUyzPDZTO/X6hfYvtz5QMB43zRrZLBnLOuTcvdtzM7pL0c5Jucgv8iznnulNfe8zsG5KukTRvIJcarbtXknbt2pXd/wEAAACs2Nyi4FkJ5PpjrI8rAucGR0y3PNuawmczu1nJ5CbvcM7FFjinwsyq0tuS3irphbXcFwAAAPkjXUsuWwlP0sXAgWK21nHQz0qqUnK65H4z+wdJMrNNZvZw6pxGST80swOSnpH0Lefct9d4XwAAAOSJ10bkMl9LbibhdGKQYuDAWrNWXrjA/m5Jt6S2j0nasZb7AAAAIH/VlJWoNODT6SzUkjs9PKHpGUeiExQ9ij8AAABgTcxMzTUhnczC1EpKDwBJBHIAAABYs6aa7BQFj6YCuQhTK1HkCOQAAACwZs01ZTqVhamV0YFxmSkr2TGBfEYgBwAAgDXbVJsckZuYnsnofaL9MTVXhxQM8Gcsihs/AQAAAFizK1tqFU84vdg9lNH7UEMOSCKQAwAAwJq1R8KSpH3HBzN6n+gAgRwgEcgBAADAAxuqShWpK9fe4wMZu8fE9IxOD0+SsRIQgRwAAAA80h6p1d7OATnnMnL9roFkwfFIPYlOAAI5AAAAeKKjLazekcnZgMtr0QFqyAFpBHIAAADwxFXpdXKdmZle2ZUuBs4aOYBADgAAAN64tKlK5UG/9mVonVxnf0zBgE8bKkszcn2gkBDIAQAAwBMBv087W5Pr5DIh2j+u1nCZfD7LyPWBQkIgBwAAAM+0R8I6eHJEsam459em9ADwGgI5AAAAeKajLayZhNOBqPeFwTv7YyQ6AVII5AAAAOCZqyK1krxPeDIUm9bIRFwRRuQASQRyAAAA8FBteVAXbKjwPOHJbOmBOmrIARKBHAAAADzWHglrn8eFwTtTpQdamFoJSCKQAwAAgMc62sIaiE3rJ31jnl0zmgrkIvUEcoBEIAcAAACPtbclC4Pv9XB6ZXQgppqyElWHSjy7JlDICOQAAADgqQs3VKoqFNC+zkHPrhntH2d9HDAHgRwAAAA85fNZcp2clyNy/TEyVgJzEMgBAADAc+2RsF7uGdHwxPSar5VIOHUNjFNDDpiDQA4AAACe62gLyzlpvwfTK3tGJjU1k1ALI3LALAI5AAAAeG5Ha43MvCkMni490BpmjRyQtqZAzsz+xMxOmNn+1OOWBc672cwOm9kRM7tnLfcEAABA/qsKleiSxipPMlfOlh5gRA6Y5cWI3N8653amHg+fe9DM/JL+XtLbJG2X9Etmtt2D+wIAACCPtbeFtb9zUInE2gqDRwdiMpM2MyIHzMrG1MprJB1xzh1zzk1J+pKk27JwXwAAAORQRySskcm4XukZXdN1OvtjaqwKqTTg96hlQOHzIpD7dTN7zsy+YGbheY5vlhSd831Xat+8zOxuM9tjZnt6e3s9aB4AAAByocOjwuBd/eNMqwTOsWQgZ2aPmdkL8zxuk/R/JV0gaaekk5L+er5LzLNvwfF159y9zrldzrldGzZsWObLAAAAQL5pqy9XXUVwzQlPogMxtVAMHDhLYKkTnHNvXs6FzOwfJT00z6EuSa1zvm+R1L2s1gEAAKBgma29MPhkfEanhieoIQecY61ZK5vnfPtOSS/Mc9qzki4ys61mFpR0p6QH13JfAAAAFIb2tlod6xtT/9jUqp5/YmBczpGxEjjXWtfIfdLMnjez5yTdKOlDkmRmm8zsYUlyzsUl/bqk70g6KOkrzrkX13hfAAAAFICOSHKd3I9XOb0yOjAuSWolkAPOsuTUysU45355gf3dkm6Z8/3Dks4rTQAAAID17cqWWgV8pn2dA7rpssYVPz9dQ66VNXLAWbJRfgAAAABFqizo1/ZN1avOXBntjyno96mxKuRxy4DCRiAHAACAjGqPhHUgOqT4TGLFz40OxNQSLpPPN18idKB4EcgBAAAgo9rbwhqfntGhUyMrfm60f1wtrI8DzkMgBwAAgIxaS2Hwzv6YWsOsjwPORSAHAACAjNpUE1JjdemKC4MPT0xraHya0gPAPAjkAAAAkFFmpo628IpH5F7LWEkgB5yLQA4AAAAZ1x4Jq2tgXD3DE8t+zmwgFyaQA85FIAcAAICMa0+tk1vJ9Mpof7IYOFMrgfMRyAEAACDjLt9UraDft6LpldGBmKpCAdWUl2SwZUBhIpADAABAxpUG/LqipUb7OgeX/Zxof4xplcACCOQAAACQFe2RWj3fNaTJ+Myyzu/sjzGtElgAgRwAAACyoqMtrKmZhF7sHl7yXOecugbG1VpHDTlgPgRyAAAAyIr2SCrhyTLWyfWOTGoynqD0ALAAAjkAAABkxcbqkFrCZcvKXNlJ6QFgUQRyAAAAyJp0YXDn3KLnRQcoBg4shkAOAAAAWdMeCev08KS6hxYvDJ6uIdcSZo0cMB8COQAAAGRNR6ow+FL15Dr7Y9pYVapQiT8bzQIKDoEcAAAAsubSpiqVlfiXTHgSpfQAsCgCOQAAAGRNwO/TjtaaJROeJEsPEMgBCyGQAwAAQFZ1tIX1UvewxqfmLww+FU/o5NC4WlkfByyIQA4AAABZ1R4JK55weq5rcN7j3YPjSjgyVgKLIZADAABAVl2VKgy+d4HplZQeAJZGIAcAAICsqqsIaltDhfYdn39ELl16gEAOWFhgLU82sy9LuiT1ba2kQefcznnOe1XSiKQZSXHn3K613BcAAACFrb0trCcO9cg5JzM761hnf0wlflNTdShHrQPy35oCOefcHeltM/trSUOLnH6jc65vLfcDAADA+tAeCev+vV169UxMWxsqzjoWHYhpc22Z/D5b4NkAPJlaacmPUX5R0n1eXA8AAADrW7ow+Hz15Lr6Y0yrBJbg1Rq5N0g67Zx7ZYHjTtIjZrbXzO726J4AAAAoUBdtrFRVaWDehCed/TG1hAkFaZKFAAAJRElEQVTkgMUsObXSzB6T1DTPoY845x5Ibf+SFh+Nu8E5121mGyU9amaHnHPfX+B+d0u6W5IikchSzQMAAEAB8vlMOyO1543IjU7GNRCbVoQROWBRSwZyzrk3L3bczAKS3iWpY5FrdKe+9pjZNyRdI2neQM45d6+keyVp165dbqn2AQAAoDB1tIX1mcdf0cjEtKpCJZKkaH+69ADFwIHFeDG18s2SDjnnuuY7aGYVZlaV3pb0VkkveHBfAAAAFLD2SFjOSQeir+XL60wHckytBBblRSB3p86ZVmlmm8zs4dS3jZJ+aGYHJD0j6VvOuW97cF8AAAAUsJ2RWplJe+dMr0yPyDG1EljcmsoPSJJz7lfm2dct6ZbU9jFJO9Z6HwAAAKwv1aESXbyxSvvmJDzpGhhXZWlAteUlOWwZkP+8yloJAAAArFh7W1j7OgeUSCRTI0T7Y2oJl51XJBzA2QjkAAAAkDPtkVqNTMR1pHdUUnKNHNMqgaURyAEAACBn5hYGd86pa2CcYuDAMhDIAQAAIGe2NlQoXF6ivccH1Dc6pfHpGbWGKT0ALIVADgAAADljZmqPJNfJpUsPROoZkQOWQiAHAACAnGpvC+to75heOJGsJ0cNOWBpBHIAAADIqfZIcp3cgwe6JUktBHLAkgjkAAAAkFM7Wmvk95n2Hh9QQ2WpyoL+XDcJyHsEcgAAAMip8mBAlzVXSZIidSQ6AZaDQA4AAAA515GaXknpAWB5COQAAACQc+2penIkOgGWh0AOAAAAOXfN1joF/T5t31Sd66YABSGQ6wYAAAAAzTVl2v3hn1F9RTDXTQEKAoEcAAAA8kJDZWmumwAUDKZWAgAAAECBIZADAAAAgAJDIAcAAAAABYZADgAAAAAKDIEcAAAAABQYAjkAAAAAKDAEcgAAAABQYAjkAAAAAKDAEMgBAAAAQIEhkAMAAACAAmPOuVy3YUFmNiLpcK7bkWM1koZy3Yg80CCpL9eNyAP0B/pCGn2BvpBGX0iiP9AX0ugLSfSH9dMX2pxzG87dGchFS1bgsHNuV64bkUtmdq9z7u5ctyPXzGxPsfcFif4g0RfS6Av0hTT6QhL9gb6QRl9Ioj+s/77A1Mr8981cNwB5hf6ANPoC0ugLSKMvYC76wzpHIJfnnHP8EGIW/QFp9AWk0ReQRl/AXPSH9S/fA7l7c90A5A36AtLoC0ijL2Au+gPS6AtIW9d9Ia+TnQAAAAAAzpfvI3IAAAAAgHNkNZAzsy+YWY+ZvTBn3w4ze8rMnjezb5pZdWp/0Mz+ObX/gJn99Jzn3GFmz5nZi2b2yWy+BnjDzFrN7LtmdjD1//ibqf11Zvaomb2S+hpO7Tcz+zszO5L6v2+fc61vm9mgmT2Uq9eD1fOqL5hZm5ntNbP9qev8Wi5fF1bH4/eGmVR/2G9mD+bqNWF1PHxvuHFOP9hvZhNm9vO5fG1YGY/fFz5hZi+kHnfk6jVhdVbRFy61ZJwxaWa/e861zotLCo5zLmsPSW+U1C7phTn7npX0ptT2eyX9WWr7g5L+ObW9UdJeJQPPekmdkjakjn1R0k3ZfB08POkLzZLaU9tVkl6WtF3SJyXdk9p/j6RPpLZvkfSfkkzS6yX9aM61bpJ0q6SHcv26eOSuL0gKSipNbVdKelXSply/Ph656Q+pY6O5fj088qMvzLlmnaR+SeW5fn08st8XJL1d0qNKlt+qkLRHUnWuXx+PjPaFjZKulvRxSb97zrXOi0sK7ZHVETnn3PeVfAOd6xJJ309tPyrpF1Lb2yU9nnpej6RBSbskbZP0snOuN3XeY3OegwLhnDvpnNuX2h6RdFDSZkm3KRmcK/U1/anpbZL+xSU9LanWzJpTz39c0kg22w/veNUXnHNTzrnJ1DmlYup4QfLyvQGFLUN94d2S/tM5F8v4C4BnPOwL2yX9l3Mu7pwbk3RA0s1ZfClYo5X2Bedcj3PuWUnT81xrvrikoOTDHzovSHpHavt2Sa2p7QOSbjOzgJltldSROnZE0qVmtsXMAkr+R7UKBcvMtki6StKPJDU6505KyR9WJT9JkZI/pNE5T+tK7cM6sta+kJpy8Vzq+Cecc93ZaTkywYP3hpCZ7TGzp5lKV9g8/D1xp6T7MtlWZNYa+8IBSW8zs3Iza5B0o/gbsmAtsy+sa4FcN0DJ6ZR/Z2YflfSgpKnU/i9IukzJYe/jknZLijvnBszs/ZK+LCmR2r8t662GJ8ysUtLXJP2Wc27YzBY8dZ59pFxdR7zoC865qKQrzWyTpP8ws/udc6cz0mBklEfvDRHnXLeZbZP0hJk975w7moHmIoO8+j2RGpG5QtJ3PG8ksmKtfcE594iZXa3k3469kp6SFM9IY5FRK+gL61rOR+Scc4ecc291znUo+SnZ0dT+uHPuQ865nc652yTVSnoldeybzrlrnXPXSTqc3o/CYmYlSv4Q/ptz7uup3afTU2FSX3tS+7t09qdmLZIYbVknvO4LqZG4FyW9IZPtRmZ41R/SI7LOuWOSvqfkJ7coIB6/N/yipG84586bYoX85+H7wsdTf1u+RcmAj78hC8wK+8K6lvNAzsw2pr76JP2hpH9IfV9uZhWp7bcoORr30jnPCUv6gKR/ykHTsQaW/Ojk85IOOuf+Zs6hByXdldq+S9IDc/b/91QmqtdLGkoPoaOwedUXzKzFzMpS1wxLukHJD3pQQDzsD2EzK01ds0HJ/vBSVl4EPJGB3xO/JKZVFiQP3xf8ZlafuuaVkq6U9EhWXgQ8sYq+sK5ltSC4md0n6aclNUg6LemPlcwu98HUKV+X9GHnnEvNe/2OktMnT0h6n3Pu+Jzr7Eg952POuS9l6SXAI2b2U5J+IOl5Jf+PJekPlJzn/BVJESWzk97unOtP/eB+VslFyTFJ73HO7Uld6weSLlWyL51Rsq8wdaZAeNUXUh/4/LWSU6lM0medc/dm9cVgzTzsD9dL+lzqGj5Jn3bOfT6rLwZr4vHviS2SnpTU6pxLCAXFw/eFkKR9qecPS/o159z+7L0SrNUq+kKTUtlJU+ePStqemo55XlxSaL8nshrIAQAAAADWLudTKwEAAAAAK0MgBwAAAAAFhkAOAAAAAAoMgRwAAAAAFBgCOQAAAAAoMARyAAAAAFBgCOQAAAAAoMAQyAEAAABAgfn/RfEZJED5JfIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "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": [ "/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/base/tsa_model.py:218: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n", " ' ignored when e.g. forecasting.', ValueWarning)\n", "/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/base/tsa_model.py:218: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n", " ' ignored when e.g. forecasting.', ValueWarning)\n" ] } ], "source": [ "mod = sm.tsa.SARIMAX(endog4)\n", "res = mod.fit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What this means is that you 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": [ "/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/base/tsa_model.py:583: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.\n", " ValueWarning)\n" ] }, { "data": { "text/plain": [ "4 0.011866\n", "dtype: float64" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res.forecast(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The index associated with the new forecast is `4`, because if the given data had an integer index, that would be the next value. A warning is given letting the user know that the index is not a date/time index.\n", "\n", "If we try to specify the steps of the forecast using a date, we will get the following exception:\n", "\n", " KeyError: 'The `end` argument could not be matched to a location related to the index of the data.'\n" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "'The `end` argument could not be matched to a location related to the index of the data.'\n" ] } ], "source": [ "# Here we'll catch the exception to prevent printing too much of\n", "# the exception trace output in this notebook\n", "try:\n", " res.forecast('2000-01-03')\n", "except KeyError as e:\n", " print(e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ultimately there is nothing wrong with using data that does not have an associated date/time frequency, or even using data that has no index at all, like a Numpy array. However, if you can use a Pandas series with an associated frequency, you'll have more options for specifying your forecasts and get back results with a more useful index." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.5" } }, "nbformat": 4, "nbformat_minor": 2 }