{ "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:22.644808Z", "iopub.status.busy": "2021-02-02T06:51:22.644040Z", "iopub.status.idle": "2021-02-02T06:51:23.666459Z", "shell.execute_reply": "2021-02-02T06:51:23.667523Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:23.672592Z", "iopub.status.busy": "2021-02-02T06:51:23.671137Z", "iopub.status.idle": "2021-02-02T06:51:23.918040Z", "shell.execute_reply": "2021-02-02T06:51:23.917634Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAEvCAYAAADW/SmEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACVxklEQVR4nOzdd3hb53k28PvF3gD3EqlNiZJsDUvykC3vkem4SeOs1tl7J23Tpm3atP26MpvRZjnbdpomTmLHe1uWbVnWntQgRYl7Ym+c749zDohNkARJgLp/16VLEgiCEHUInOc8S0iSBCIiIiIiIipPmoV+AkRERERERJQfgzYiIiIiIqIyxqCNiIiIiIiojDFoIyIiIiIiKmMM2oiIiIiIiMoYgzYiIiIiIqIypluIL1pbWystW7ZsIb40ERERERHRgnv11VdHJEmqK+a+CxK0LVu2DHv37l2IL01ERERERLTghBDnir0vyyOJiIiIiIjKGIM2IiIiIiKiMsagjYiIiIiIqIwVHbQJIe4WQgwJIY6k3PYPQoheIcQB5ddr5+ZpEhERERERXZymk2n7CYDbctz+dUmSNim/HirN0yIiIiIiIiJgGkGbJEnPARibw+dCREREREREGUrR0/ZxIcQhpXyyKt+dhBAfFELsFULsHR4eLsGXJSIiIiIiWvxmG7T9N4CVADYB6Afw1Xx3lCTp+5IkbZUkaWtdXVE75IiIiIiIiC56swraJEkalCQpLklSAsAPAGwvzdMiIiIiIiIiANDN5pOFEE2SJPUrf70DwJFC9yciorlz4PwEzgz5YDXqYDPqYDVq0d5gh9U4q5d6IiIiWmBFv5MLIe4FcB2AWiHEBQBfAnCdEGITAAlAN4APlf4pEhFRMT7ws70Y9obTbrthbT3ufve2BXpGREREVApFB22SJL09x80/KuFzISKiGfKEohj2hvHha1fiDRub4A/H8e+PnEDfRHChnxoRERHNEmtmiIgWgXMjAQDAplYX1jc7AQDLaqx46ezoQj4tIiIiKoFSjPwnIqIF1jXqBwAsr7Umb3OYdfAEowv1lIiIiKhEGLQRES0C3SNy0La0xpK8zWHSwxuOIZ6QFuppERERUQkwaCMiWgS6R/xocppg0muTtznMegCALxRbqKdFREREJcCgjYhoEege9WNZjTXtNodJblv2hFgiSUREVMkYtBERLQLdowEsq80I2pRMm5t9bURERBWNQRsRUYVzB6MY80ewLKWfDZB72gBm2oiIiCodgzYiogp3TpkcmZ1pU8ojg+xpIyIiqmQM2ohoXrkDUYSi8YV+GotK10j2uH+AmTYiIqLFgkEbEc2rP/nvF/C1xzsX+mksKt3KYu226ozySKWnjbvaiIiIKhuDNiKaN9F4AmdH/Dg54F3op1KRBj0hSFL2zrXuUT+aM8b9A4DdqIMQgIcj/4mIiCoagzYimjfD3jAkCeh3Bxf6qVScEwMeXPmvT+LxY4NZH+se9Wf1swGARiNgM+qYaSMiIqpwDNqIaN4MeEIAgP6J0AI/k8rz8OEBJCTgqRNDWR/rHvFjaU120AbIfW3saSMiIqpsDNqIaN4MuOVgzRuOwctAYlqeOC5n2F44M5J2uzsQxXggiuW1llyfBodZz+mRREREFY5BGxHNGzVoy/wzFdbvDuJonwdLayw4PxZEz2gg+bFuddx/3kybjpk2IiKiCsegjYjmzaBnMlDrY9BWtCeOyyWRf/PaDgDp2bbuPDvaVHKmjUEbERFRJWPQRkTzZsATgkEnv+z0T3AYSbGeODaI5bVW3LKuAQ0OI3adngzaukb8ECJ73L/KYdLDy+mRREREFY1BGxHNmwF3COuaHBAC6GemrSi+cAwvnhnFjWvrIYTAjpW1ePHMKBIJefT/udEAmp3mrHH/KoeZ0yOJiIgqHYM2Ipo3g54QWqstqLUZOfa/SM93DiMST+CmdQ0AgB2rajHmj+CEsuuua8SPpTW5s2yAkmkLxxBPZO93IyIiosrAoI2I5oUkSeh3h9DoMKLZaWKmrUiPHx+E06zH1qVVAOSgDQBeUEok8+1oUznMegCAjyWSREREFYtBGxEV7cSAB1959CQkafpZG3cwinAsgQaHCU1OM/rY0zaleELC0yeGcMPaeui08st1o9OElXVWvHBmBBOBCCYCUSzPMzkSkKdHAuAESSIiogrGoI2Iivbw4QF8++nTuDA+/YBLXazd6DShySVn2mYS/F1M9vWMYzwQxU0dDWm371hVi5fPjuHUkA8ACpdHKpk2N/vaiIiIKhaDNiIqmj8sl9gd6/dM+3PVvWyNDhOanCYEInF4WLJX0BPHBqHXCuxsr027/aqVtQhG47h/fy8AYHmh8kiTHLQx00ZERFS5GLQRUdH8kTgA4Gjf9IO2wdRMm9MMABxGMoXHjw/iihU1sCuBl+rKFTXQCOD3+3shBNCaZ9w/IE+PBABPkAEyERFRpWLQRkRFS2baZhC0DbjDAIB6uwnNLhMAoH+Cw0jyOTvsw9lhP25e15D1MadFj0tanPBH4gXH/QPMtBERES0GDNqIqGhq0HZ8JuWRniBqbQYYdJpkpq2Pmba8Dl6YACBn1XJRp0guq82fZQMme9q4q42IiKhyMWgjoqL5lKCtdyKIcX9kWp874A6hwSFn2OrtRmgEM22FDHrkzGSTy5zz48mgrcDkSACwG3UQAuwfJCIiqmAM2oioaIFIHCa9/LIx3WzbgCeMRiVo02k1qLdzV1shg54QbEYdbEZdzo9ftrQKy2utuDxPJk6l0QjYjDpm2oiIiCoYgzYiKpo/HMPmVnnJ83QnSA56QmhwmpJ/l8f+szwynyFPGPUOY96Pm/RaPP356/DGjc1TPpbDpGdPGxERUQVj0EZERfOFY1haY0GDwzitCZLhWBxj/kgy0wYAzU4zM20FDHhCaLCbpr5jEZxmfVHTI08MeHDr15/DiC9ckq9LREREpcGgjYiKFojEYTHosK7JMa0JkkNKf1ZjaqbNaULfRJALtvMY9ITSvl+z4TDrisq0PXSoHycHvTjc6y7J1yUiIqLSYNBGREWRJAn+SAw2oxbrmh04PexDKBov6nP7UxZrq5pcZoRjCYwHWLaXSZKkKcsjp8Nh0hfV0/Zy1xgA4MJYoCRfl4iIiEqDQRsRFSUQiUOSAKtRh/XNTsQTEjoHvUV97kDKYm1Vk/Jn9rVlGw9EEYknSlYe6TDr4Z1iemQoGsf+8xMAgB4GbURERGWFQRsRFcUfkU/6LUa5PBIofsn2oJJpa3DkCNo49j/LYI4gdzaKybQdPD+BSCwBADg/xkCaiIionDBoI6Ki+MNyKaTNqEVbtQU2o67oCZIDnhDMei0cpsnx9c3K/jFm2rKpQVtDqcojzTp4wzHEE/n7B1/uGoMQwJY2F86PM9NGRERUThi0EVFR/MpibatBB41GoKPJXnSmbUAZqiGESN5WazNCpxHo4wTJLGrQVl+q8kiTHgDgK1Ai+XLXKNY02LGhxYnzLI8kIiIqKwzaiKgoatCmLnte1+TA8X4PEinZmyFvCB/95atZJ/2D7lDaEBIA0GoEGhwmDDBoyzKoTNss2SASsxy05ZsgGYkl8Oq5cVyxogatVRZ4QjG4OSCGiIiobDBoI6KipPa0AcD6Zif8kTjOKQGaJEn4wm8O46HDA/j5S+fSPncgz/h6dew/pRv0hFBtNcCo05bk8dSyVHeevrbDvW6EoglcvrwardUWAGCJJBERURlh0EZERfGl9LQBwLrm9GEk971yHk+dGILLoscfDvQl+6cSCQmDnlDaEBJVk4sLtnMZ9IRQby9Nlg2YOtP2ctcoAGDb8mq0Vsu9hiyRJCIiKh8M2oioKMmeNiXTtrrBBp1G4GifGz2jAfzTg8dw1coafPn2DRjwhPDyWTkQGAtEEI1LaMxR6tfslMsjEwUGZFyMBj3hnEHuTKk9bZ5g7p62l8+OYVW9DbU245SZtu4Rf9H7+YiIiKg0GLQRUVHUoM1ikIM2o06LVfU2HO5143O/PgCtEPjPP92IW9Y1wGbU4XcHegEg2bOWrzwyEk9g1B+Zp39FZRj0ZPcAzobDLP+f5cq0xeIJ7O0ew+XLq+X7mvRwmvU5d7X5wjHc+o3n8JPd3SV7bkRERDQ1Bm1EVBR15L/VMNlnta7ZgedPjeCV7nH84+3r0eIyw6TX4tb1jXj48ABC0XjK+PrsIKTRKZficRjJpFg8gRFfuGTj/oGU8sgcPW3H+j3wR+K4fEVN8ra2akvOXW3H+jwIxxI40usu2XMjIiKiqTFoI6Ki+CMxmPQa6LSTLxvqku3b1jfijs0tydvv2NwCbziGp04MYaDAouhml3xbH3e1JY34IkhIQH0JM202gw5CAJ4cI/9fPjsGALhCybQBQGu1OWd5pBqsnR7yzfi5nBjw4B8fOApJYkksERFRsYoO2oQQdwshhoQQR1JuqxZCPC6EOKX8XjU3T5OIFpo/HIPVoEu77db1jbh9UzP+5Y4NaTvYrlxZg3q7Eb/b34tBdwgaAdTZsjNHTUqmrZ8TJJMKZSZnSqMRsBt1OTNtL3eNYnmtNS1IbK2y4MJ4MKvX8EifHLSdHfEjFk/M6Lk8cWwQP36hG+NcKUBERFS06WTafgLgtozbvgDgSUmSVgN4Uvk7ES1C/nAsOYRE1VptwTffthk1GQGZViPwho3NePrkEE4MeFFnN6Zl6FQ1VgMMWg0nSKZQg7ZS9rQBcolkZk9bPCFhT9cYti+rTrt9SbUFkVgCQ95w2u1Hez3QCHmv2/nxmQXa6hRStUeSiIiIplZ00CZJ0nMAxjJuvh3AT5U//xTAm0rztIio3PjC8aygrZA7NrcgGpfwxPHBvAGIRiPQ6DShj0Fb0mSmrXQ9bYA8YCRzeuSJAQ88oRguX5EetLXlmCAZjMRxasiLHatqAcy8RDKg7PtT9/4RERHR1Gbb09YgSVI/ACi/18/+KRFROQpEYmlDSKayvtmBlXVWJKTCpX6NDhMGGbQlDXrC0GpEVvZythxmXVambU+XfB0udQgJALRWZe9qOz7gQUICbt8k9y6eGvLO6Hn4lAwbM21ERETFm7dBJEKIDwoh9goh9g4PD8/XlyWiEslVHlmIEAJvUk7wcw0hUdXYDBj1h/N+fLE41ufB7jMjU95v0BNCnc0IrUZMed/pkDNt6UHb3nPjaHaa0OIyp93eUmWGEEgb+39UGUJy5coaNDpMOD04s0ybPxm0cdcbERFRsWYbtA0KIZoAQPl9KN8dJUn6viRJWyVJ2lpXVzfLL0tE880XjsE2jaANkLMyQkyW2+VSbTVg7CLY0/bXvz2Ed/zgZXzwZ3vRW2DwyoAnVPLSSEDuafNmTI/cd24cl2X0swHyDr4Guylt7P+RXg+qLHo0O01Y3WDD6eGZBm3saSMiIpqu2QZtfwBwl/LnuwD8fpaPR0RlKhCJwzKN8kgAaKux4IGPX413XN6W9z41NiMmgtEZTyOsBImEhJODXqxttOP5UyO46avP4r+fOYNILPvfPOQJl3Tcv8ph0sOdkmnrnQii3x3CZW2unPdvq7ak9bQd6XNjQ4sTQgisrLPh9JAva7pkMdReNh+DNiIioqJNZ+T/vQBeBLBGCHFBCPE+AP8G4GYhxCkANyt/J6JFyDfN8kjVhhYnLIb8n1djNUCSsKhHwJ8fDyAUTeA9O5bh8c/uxDWra/Hvj5zAu3+8J+u+g965yrTp4AvHksHxq+fGAQBbc2TaAGBJtRkXlPLIcCyOzkEv1jc7AQCrG2wIROIz2q+nZtgCEZZHEhERFWs60yPfLklSkyRJekmSlkiS9CNJkkYlSbpRkqTVyu+Z0yWJaBGQJAn+GZRHFqPaagCARV0i2an0f61usGNJlQXf//Ot+NSNq7H7zGjasI9QNI6JQLTk4/4BOdMGTGa49p0bh8WgxdpGe877t1ZZ0O8JIRyL49SgD9G4hA0t8jL11fXy58xkgqRaHslMGxERUfHmbRAJEVWuUDSBhARYjNMrjyxGjRK0LeZhJJ2D8qTF1fW25G1v2NgMAHju1ORgpiGP/D2Yk/JIsxy0qWP/954bw6ZWV879eYC8g0+SgL6JEI4oQ0g2KJm2Vcq/YyZBmy+ZaWPQRkREVCwGbUQ0JbUPaU4ybbaLIdPmRYvLDLuS7QKAlXVWtLjMePbkZNA26FV3tM1Fpk3+v/OEovCHYzje78VlS6vy3j+5q20sgCN9btiNuuRt1VYDaqwGnJrmBEk1YwtweiQREdF0MGgjoqQfPn8WP3mhK+t29UTbWqA3baYuhvLIkwNetDfY0m4TQuC6NXV44fRIciCJulh7Tsojk5m2KA6en0A8IRUM2lqrlV1t4wEc6fVgXbMDmpQ1BKvqpz9BMhJPIKYML2F5JBERUfEYtBFR0v+9egH37+/Nul09wbbOQXlktUUpj/QtzqAtFk/g7LAf7Q3ZvWPXttfBH4knh4IMuNVM2xwMIlGyfJ5QNPn1NrflD9oa7CYYtBp0DftxvN+DDS3OtI+vbrDh1KAXklT8BMnU7BrLI4mILk737unB+CK+UDtXGLQRUdKoP4KRHMGTOulvJtMjp6LTauCy6BdtT1v3aACReCJn0HbVqlroNALPdsolkkPeMAw6DZxmfdZ9Z8thVsojgzG82jOO9gZbwa+j0Qi0VJnxTOcwwrFEcgiJalWdDZ5QDMPe4v/fUnez+VgeSUR00embCOKvf3sYDx7uX+inUnEYtBERAHmX2Jg/krNMcTLTVvqgDVjcC7bVISS5gjabUYety6qSQdugslhbCJF139lSyyMnghF5qfbS3KP+U7VWW5LDRtQhJKrVDdOfIJlaEsnl2kREFx9vSO1r5nvAdDFoIyIAgDsYRTwhIRiNZ5WuqS+uczGIBJAnSC7W8sjOQS+EmJy4mOm6NfU43u/BoCeEAXdoTvrZAMBm0EEIYN+5CXhCsYL9bKrWKrmvzaTXYEVd+vNXJ2GemkbQph5HJr2Gb9hERBchX1jeyRrge8C0MWgjIgDpI/czA6iAUspmMZS+pw0ov0zboCeE3+67UJLH6hz0oq3aAnOe79217XUAgGc7hzHkDc/JuH9ALne0G3XYdXoEALC1mKBNmRa5rskBrSY9+1dnN8Ju0uHUkLfo5+BXymzr7abkRFKiTHfv6sI3nzi10E+DiOaAmmlT2y6oeAzaiAgA0nrZMgMo3xxn2qqtxrIK2n75cg8++78HS9Io3Tnoy1kaqVrbaEeDw4hnTw7L5ZH2uQnaALlE0heOocZqwNIay5T3V0f8Zw4hAeTpl6vrbdMqj1Sza/V2I0f+U14PHOrDHw5mD0Qiosqnvvb7GbRNG4M2IgKQnl3LDKDUk23LHIz8B4BamwHjgQgSieInEc6l3vEgAKDPHZzV44RjcXSN+LPG/acSQuDa9jo8fXIIgUgcjc7ST45UqRMkL1taVVTf3LIaKwDgkhxBGwCsrrfPqKetwWFieSTlNeQJYzwQXeinQURzQC2PDFZYtYUkSYjFEwv6HBi0ERGA9PLIEV/6REB/JA6DVgODbm5eMqqtBiQkYCJYHidqvRMBAED/RGhWj9M14kc8IRXMtAHAte31yVKRuVisrVInSBbTzwYAHU12fP/PLsPtm1pyfnxVvQ0jvkjRGUk1UKuzGxGOJRb8DZDKTyIhYdATwkQggniZXMQhotJJDiKpsEzbL146h2v/85kFfQ4M2ogIQOHySH84Nic72lTqgu1RX3mM/e+dkDNs/bPMtJ0cyD85MtXVq2qhtozVz2V5pJJp27qsuKBNCIFb1jfmDdZXKRnEYpdsZwamLJGkTGOBCGIJCQlJHo5ERIuLWnERrLCgrXPQh96J4IJeTGLQRkQA5ICp2mqAQafJE7TNTWkkANRY5ZLA0TLoa4snpGSGrc89u0zbqUEftBqBFXXWgvdzWvTYoiy6novF2smvY9bDoNVgfXPucsfpSk6QHCwuaPOFY9BrBaoscvDIYSSUadAz+TNXTn2uRFQavmSmrbJe/9WLSOHYwgWbc3cWRkQVZdQXQa3NAKNOk7Vg2xeOwTpH/WzAZKatHE7Shr1hxJQraf0Ts8y0DXqxrMYCo27qLOXN6xpwtM+DJqd5Vl+zkLuuWoYdq2ph0pcma9rsNMOs1xY9QVIN/tULAOxro0xDnslsezm8HhBRaanBWqVl2tT2jVA0AYthYZ4DgzYiAiD3sdVYjUqmLb1MMRCJz2l5ZI1NKY8sg5M0tZ9NI0qRafNiXbOjqPu+7+rleMPG5ryrAUphQ4sz5yTImdJoBFqrzegrMrhVg3/1WPIxaKMMzLQRLW7eSs20BeTXo4XMtLE8kogAyAFTjc2Qc/y+b47LI6uUy1ZjZbBg+4IyOXJto2NWPW3BSBznxgJYXV+4n02l02rQ7Jq7LNtcqbebMOgprhdR7Y1Us7bc00OZBplpI1rUKrWnzZ2SaVsoDNqICICcaau1GVFjNWRlvPxzXB5p0GngMOmyMnwLQR1CsnVZFQbcoRmvITgz7IMkAWsaiwvaKlW9w4hhb3H/b3LGdrI8kpk2yjToDcGuHB/jAQZtRItNsqetwgZRTZRBTxuDNiJCOBaHNyQvXa6xGtJ2tgGTJ9tzqcZmxEgZXFnvmwjCZdFjVb0N0biEkRkGkpOTI/PvaFsM6u0mDHlDkKSpg1tfOAYbe9qogCFPCEuqLbAatFmvQ0RU+ZKZtmi8bHazTiWRkJhpI6LyoJYh1diMqLYZEIzG00oX5JPtueu1AuRhJOVQHtk7HkSLy5wcCNI3w11tnUNeGLQaLK0pPDmy0tXbjYjGpaKWIfvDMVgM2mRPW6Xt6aG5N+gJo8Ehvw4x00a0+KRWWASjlfEe4A3HoF6XDC3gc2bQRkTJK9o1NjnTBkwu25YkST7ZnuNMW7XVUBY9LL0TatAm7xKb6QTJzgEvVtRZodcu7pdZdefakHfq4NYfljO2NmbaKI9BTwgNdhOqLdll2kRU+XzhGLTKYtJK6Wt2p1yUDMeYaSOiBTSiLLWutRkmd6YpgVwknkAsISVPtOdKrl66+SZJEnrHg2h2mZNDQWYyQXLEF8b+8xOLvp8NkHvaABQ1jEQtjzTrtRCCQRuli8UTGPEpmTarAeMM2ogWFUmS4AvFUGeT3zcCFTJBUi2NBJhpI6IFpgZotUp5JDBZMqk2C1vncBQ9IGf5xgORrBr33x/oxY5/eyrtRXOueIIx+CNxLKkyo8qih1GnmXamLRZP4BP37EcwEscHrlkxR8+0fDTYlUybZ+rgNhCRp5AKIWA16CquEZ3m1qg/goQE1DtMqCqTzDsRlU44Jl8EVi/2VUqmbSI4+VrETBsRLSi1FLJGmR4p36YGbfKVsLkvjzQinpDgCaUHZ891jqB3Ioj79vTM6dcHgAvKjrYWlxlCCDS7zOifZqbtPx89iRfPjuJf7rikpDvRypX65js0xQTJcCyOaHwyY2s1aplpozTqjrYGhwk1DNqIFh21n63eXlmZtokAM21EVCZGfREYdRpYDVrU2NTySPkkXF2AOR/lkQAwkjGM5Hi/BwDw4xe6EZnjK1y9yo62liq5NLLJaULfNHa1PXS4H9977izedUUb3nLZkjl5juXGpNfCbtJNmWlTs2oWJWNrNegqbrkqzS21xLbBYUSVNXsgEhFVNnXcf51SoVEpmbbUSp8wgzYiWkgjvghqbUalbE0Lg06TUh4pv8jO9cj/amt6WSYAROMJnB7yoaPJgQFPCH883Denz6FPKYVU+9manGb0Fzk98vSQF3/x64PY3ObC379+/Zw9x3LU4DBNmWnLPI6sRh0zbZQmM9MGAGOcIEm0aKiZtgalQqNSSuTTgjaWRxLRQhrxhVGj9LIJIVCbMhTEN089bZNB2+TJ/9lhPyLxBD5wzXKsrrfhB891FbUPbKZ6J4Iw6TXJE8Zml7yDLBYv/CIdT0j48C/2wWzQ4r/feRkMuovrpbXebkyecOeTmbG1GLQV84ZN82PIE4JGyFn3KovyelAGa0CIqDS8IbU8Us60BaOVceHOHYxCr5UnXrI8kogW1Kg/nAxUAKDaZpgsj5ynTJsaNKZOkDwxIJdGrmt24P3XLMexfg92nxmds+fQOyFPjhRCfnFucpqRkIDBKbJIp4d8OD3kw1/cugaNyqqAi8lMMm02I8sjKd2gJ4xamxE6rSb5esBMG9Hi4c/oaauUC3cTgQiqrQZoBDNtRLTARn2RZC8bIA8FySyPnOuetmSmLeXK+rF+D/RagZV1Nty+qQW1NiN+8PzZOXsO6mJtVZOruF1tB89PAAC2Lques+dWzurtRgx5wgWzoGrGVl3SzvJIyjToDSX3/iUzbf6pV0kQUWVIDiJRyiMrpWd1IhCFy2yASa9lpo2IFo4kSUrQNplpSy2PnK9Mm1Gnhd2oS8+09Xuxqt4OvVYDk16Lu65cimdODqNz0Dsnz6F3IpQWtDU7i9vVduDCBBwmHZbXWOfkeZW7OrsRkXii4FqG5BRSw+T0SF+FXGWl+THoCSd7XdR9kWP+uV/1QUTzwxtOL4+slGoLdzAKp1mvBG3MtBEtqD8e6r9ox0t7wzFE4gnUWlMzbYbk7jZ/JH3q31yqtqWP+T7e70FH0+SC6nddsRQmvQY/nINsWygax4gvPKNM24GeCWxsdUGjESV/XpVAzY4UKpH0ZWRsrQZdxYx7ptkpdsfikCeEeuVYspt00GoEM21Ei4g6PdJplvegVkqmzR2Mwqnsbg3HmGkjWjB9E0F87J59+NZTpxb6qSwINThLzbRV2ybHbfvDMeg0AsZ5GK5RbTUkd8aN+sIY8obR0ehIfrzKasBbt7bi/v29+PvfH8GTxwdLVmKnTo5Ux/0DgMOkh92oK7irLRiJ4+SgFxuXuEryPCqR2p9QaBhJIMf0yEAknrVMnRaXIU8I2/75CTx5fLDg/SKxBEb9keSydo1GoMqiZ6aNaBHxhaPQCMCk18gl8hVy4c4djMLFTBvRwlP3gD16ZGBOJxOWK3XgSGpPm5p1G/WH4Q/HYDXqksM55lJNSobvxIBcAtnR5Ei7z8dvWIXr19Tj13sv4H0/3YtNX34M7/rhyxiY5hLsTL0Z4/5VTS5TMqDL5WifG/GEhE2trll9/UqWzLR58mdF1IytNdnTplVur4w3bZqZrhF5Auzec+MF7zfsm9zRpqq2GphpI1pE/OE4bMr5hFmvrZg9bROBaDI7yEwb0QJSg7Y+dwgHL7iL+pye0QAePNQ35ULhSqAus06bHqn8edQXgS8cn/Nx/6lfVy2PVP9fUssjAbkW/vt/vhUHvnQz7nn/5XjbtjbsOj2C3WdG8j5u14gfx/s9BbM6yUxbZtDmNBfMtB1QhpBc2urM/w9b5NSm8kFv/u+TT8nYGrTy246acauUN22amQHlNbJzoHAfauqONlWVxYBxZtqIFg1vKAa7SQ9AvnAXKKO+5hMDHnSP+LNuD8fiCEbjcFn0MC5wpm1uJwsQVYDjA17U2oxwByN4+HB/URmTz//6IPZ0jwEAVtRZceWKGty6vhE72+vm+NmWnlqOWGdPucJtm1x0HYjE5nwIiarGZsR4IAJJknC834s6uzEtA5jKqNPiqlW12LDEiZ+/dA4jvvxX5N/1w5fROxFEtdWAy5dX46qVNXjNJU2oTXns3vEgNAJZI/ubXSYc7csfzB84P4EWlznZWH0xshh0sBl1hTNtGRlbtbfNF46hYV6eJS0ENRg7OcXwIPUCWH1Kpq3GZsDJKYI9IqocvnA0+dpvNugQWMBJjJk+fd8BNLvMuPvd29JuV3tynRYDTDoNp0cSLaTj/R5sbnNhx6paPHSkf8oSyVODXuzpHsO7r1qGL762A0urLfjd/l68+8d74A1V3lVhtRxRHbENpJZHRuALz2PQZjUgGpfgCcVwYsCTVRqZi92og0GnSWYMM8UTEvrdQVzbXofr19Tj0AU3/u73R/GOH7yEeErm7cJEEI0OE/Ta9JfFJqcZI75I3pKIgxcmsPEizrKp6h1GDE8xiCR1bYQ6RZJj/yvLrlMj+NrjnUXff8AtHxMXxoPJYTS57ycHbY2ZmbZA5b2mElFuvnAMNpM6jEqb7HVeaImEhK4RPy6MB7I+5lZeg5xmOdPGPW1ECyQYiaN7xI+OJgdes6ER58eCONrnKfg59+45D71W4OM3rMIHdq7Aj9+zHd96x2YkJFTkVeFRXxhOsx6GlEEjaqZt1Kf2tM1feSQgX3U/NehDR6N9is8AhBCosxkxkidgGPNHkJCAGzvq8dW3bsSuv7oeX79zIzoHfXjwUF/yfr3jwax+NgBoUjJvuXrmRn1hnB8LXtT9bKp6u3GKQSTxtOMo2dNWRuUxNLWf7O7Gfz15ChNFLr1OLZkttKpj0BuGXivSLh7VWA0YD0TSLq4QUeXyhSYvAlsM5dPTNuQNIxxL5GyFUDNtLrOemTaihdQ56EVCAjoa7bh5XSO0GoGHj/TnvX8oGsdv9l3Aresb00rr1IyQ2odVSUb86TvaAPkKmEGnwZg/An84DqthfjJtatD2Svc4IvFEUZk2AKi1GZKDDDKp2Z865f9LCIHbN7ZgTYMd//XkqeQJYZ87mDY5UqUGcr05hpEcvDABABf15EhVg8NUcOS/P6PMVs26MdNWWQ4px/y+nsKDRVSD7sndh4X62gY9IdTbTWlrM6qsBkhS8SsDiKi8+cIx2JNBW/msfekZkzNs3lAsqyJgIiPTFmGmjWhhnBhQh104UG014MoVNXjocP4pkg8f6Yc7GMU7trel3d7oMMFp1uNYf+Vl2ka84bQdbYAc2KgLtv2R9LK2uaQu1H1BGSqytmnqTBsA1NqMecsj1WAutWdPoxH41E2rcWbYjwcO9skllBmLtVVqpq1/IvsK3IHzbmgEsKGF5ZFqpi3fz44vHEsL/pPlkWXypk1TG/SEkoH53u7igrYBTwhbl1XBYtAW7Gsb8oTT+tmAyYs4F+sOTaLFJrVMvpwybWrQBgAD7vQLtBNqps3CTBvRgjre74XFoEVbtQUA8JpLGtE14s97cnHPyz1YVmPBlStr0m4XQqCjyV6RmbbRHJk2YHLRtT8cg2W+yiOV57H79AgMWg1W1tmK+jw5aCucaavNGGhy2/pGrG2Us2397iBiCSlnpq3JKd/W786RaTs/gfYG+7z1/JWzersJ4VgCnlDuICyzzHYy01Yeb9o0tYPKpFSLQTvlCH8AkCQJQ54wGp0mrG6wFywfH/SEkjvaVAzaiBYXX2iyp03OtJXH63/P6OTUyMwSycnySAOMeg1CzLQRlc4DB/vw2m8+X9QujWP9HqxptCdLcm5Z1wghgIcPD2Td99SgF690j+Pt29ty7izraHLg5IC34pYFj/rCuYM2q1HpaYvP6yASABgPRLGq3pY1FCSfWrscYOb63o/kyLQBcrbt0zetxtkRP/7n2TMAsne0AYDZoEWVRY++jBdySZJw8MIE+9kUapZkOM/Y/8zjaLKnjZm2SnHoghtajcDtm1pw8PzElGVC44EoIvEEGuwmrGmwFe5p84TSdrQBDNqIFpNEQoI/Es/ItMXKYj9uz1gguY4mK2gLRCAEYDfpYNJpEWamjah0nj45hGP9Huw+PVrwfpIk4UR/+oTCOrsR25dV5+xru2dPD/RagbdctiTn43U0OhCMxnFuLHv6ULmKxRMYD0STZYmpaq0GDHrCiMQTsM1TT5tJr03uhCu2NBKQs2jxhITxHMMRhr1hWAzanIHnLesa0dHkwC9f7gEALMkRtAHKrraMnrZzowFMBKLYyKANAJIrDwbzjP3PLLNVyyMLTRSk8nKo1432BjuuXlWLcCxRcBUGkDIR0mlCe4MdI75Izox4MBKHJxRDvYOZNqLFSi2FTwZtRi0SEhZ0GqOqZyyAS5bIbQ6ZQ8cmglE4THpoNAImvZaZNqJSUq/mFhooAsjLtD2hWNaEwtde0oTOQR9OD01eFQ5F4/jtvl7cur4x796wShxGMqYEObU5M22G5GJcyzyW/6klkuuKHEICTJY+5uprG/aGs7JsKo1G4FM3roZ6oS9XeSQg72rLvPqmDiFhpk2mZkmG8mba0geRaDUCZr22bBrRqTBJknDowgQubXFi67IqAMCrU5RIpi7MXtso/zznyrapx0xDRtCmTpLMdTGGiCqLeoEuWR6ply/QlkOJZM9YEKvrbaixGnKWRzrN8kJwo06DeEJCLL4wgRuDNlpU4gkJpwZ9AIDHjw0W/ME60T85hCTVbRsaAQCf+dVB/OvDx3H//gu4+4UueQDJ5W1Zj6Na3WCDViMqKmhTd7TlCkSrUwI52zz1tAFyWSaA5EleMSaDtuyr+MPecFY/W6pb1zdgXZMDNVZDMvuTqclpRl9Gpm1/zwTMei1W1xfXd7fYqVmSXJm2cCyOaFxKZlFVVqMOPva0VYQL40FMBKK4tNWJBocJS6rMUw4jUYO2RqcJ7Y3yz0muvjb1mMksj1Qz76N5hgwRUeXwhTIzbfLvC33hzh+OYcQXRluNBU0uU1b/+kQgCpdFDtpMSqC5UNk2ds/TonJu1I9wLIGbOurxxPEhvNw1hh2ranPeVw2u1mRk2hocJvzNa9fi/v19uHtXF6JxOQ2zvNaKK1fUZD2OyqTXYkWtFccraIJkMmizZmfaUidKzuegDfW5dEyjPLLOLn9OzqDNF8aqAgNNhBD4zju3ZJU/pmpymeAJxdKyRQcvTOCSFid0RfbdLXY2ow4WgxZDOYI2ddhI5nFkNWrZ01Yh1MzypS0uAMC2ZdV4/tQIJEnK2eMLIJmpr7MZlR1s+pyZttSMXKZqm4GZNqJFwJuZaTOUR6btvLJQu63agkaHOWvBdlqmTS+/34ei8Xmbqp2qJF9RCNENwAsgDiAmSdLWUjwu0XSpJwQfuGYFXjg9ioeP9OcP2ga8aK02w27SZ33sgztX4oM7VyIaT6BrxI8TA16sabDnPTlRrW1yYF8RU9XKxahfPsHOmWlLCeTma08bIL9wtlab85ah5qJm0oZz7Akb8YULBtuAHJAvr7Xm/bi6CuB1//U8OpocWNvowNE+D+66cmnRz/FiIO9qyy6PVAOzrKCtjPb0UGGHLrhh0GqSF7kuW1qF+/f3omcsgKU1uX92Bj0h1NoMMOjkE532PBMkk0GbPUfQZpFXjxBRZVMzbeqeNvW8YqGDtp7RyaCtyWnCK91jaR93B6NYorROmHRyoLlQfXilvER8vSRJmxiw0UI6OeCDEMClS1y4bk0dHj06mHea4/F+DzqmKMHTazVob7DjjRubszJyuXQ02dE7Ecy5DLYcJiRlUnvA6qYoj5zPTNvnb12D//vwVdP6HKdZD71WZPW0hWNxTASieXvainVTRwM+e3M71jTKax2+8WQnIrEELl9eOBi82NTZjbkzbRkN6CqbUTcng0gisQQmFjg7E09IWQ3tlezQhQl0NNmTAZja11aoRHLAHUrLnq1ttKNz0Jf1WjjkDcOo08Bhzn6dqbYaMM6gjajiZV68M6uZtgWutlB3tC2ttqLRaYI7GE27mDgRiCTLI1MzbQuB5ZG0qJwc9GBptQVmgxa3bWjEw0cG8GrPOLYtq067XzASR/eIH6+/tLmkX1/tjzs54MX25ZNf88ywD7d/+wVUWfVY2+hAR6MdHU0OXL+2PlkjvRBGfWHoNCLnyVJ6eeT8PUebUTftsgMhBGqs2bva1PLPQj1txbAadfjkjauTfw9EYuh3h7CiQHbuYtTgMOGwUkaXSn2ztmT0tFmM2jmZDPgPDxzF7/f34tcfvgrrmovvjSyle/f04J//eAz7/u7mvL2SlSKRkHCk14M7Nrckb2uvt8Nu0mHvuXG8Oc9E3UFPOLmcHgDaG+3whWPonQhiSZUl5X5ycJerkqHKakCn0qdMRJUrWR5Zbpm2sQAcJh2cFj2aXfLrVb87hJV1NiQSEtzBKFxm+SK2Uc20RSs70yYBeEwI8aoQ4oMlekyqAD9/6Rx+8+qFhX4aSScHvGhvkDNiN6yth0GrwUOHs6dIdg56kZCAddPomyqGmrnLHEbyi5fOIRyL49IWF84O+/Dtp0/jI7/ch3uUUfO5nBjw4FSBvUalMOqTF2vnOllKy7RVwElnrd2QFbSp5ZKzzbRlshh0WFlnm7Jc9mJTbzdi0BPOyqSow0Yyg3HrHGTaQtE4HjjQB38kjvf99BUMeRYm27WvZxyhaALePMvGK8nZET984RguVUZiA/Lk1S1tVXj13Fjezxv0hNLG+K9RXpsz+9py7WhT1VgNHPlPtAgkyyNN6Zk2/wKXyPeMBdBWI19EanTIZZBqlYQvEkNCQrKnzaRm2orYAzwXShW07ZAkaQuA1wD4mBBiZ+YdhBAfFELsFULsHR4eLtGXpYX23adP4ye7u/N+fE/XGL7wm0PzUhoYisbRPRpIljHaTXpcs7oWjx4ZyPr6JwbkoGo6EwqL0eAwosqiTwvaQtE4fvPqBdy2oQnfeecWPPm563Dsy7fBadbj7Ej+K8h/9ZvD+MvfHCrp88s04gvn3NEGAFaDNlkKNZ/lkTNVZ8vOtOVbrE1zo95uRDAazwrE8ve0lX4QyTMnh+ENx/CF16yFOxjF+366tyR9cz97sRvPdhb/3qUOJFroq8ilcEgdQrLElXb71qVV6Bz05SxFDcfiGPVH0JgStK1WgraTA5Ove55QFKcGfWh05l63UWU1IBiNI7gIvo9EFzNfxvuAWsGz0D/bPaMBLK2Wq2bUygB17L87ILe6ONXySCXTtlDlkSUJ2iRJ6lN+HwJwP4DtOe7zfUmStkqStLWurq4UX5YW2Jg/gn53COdG/Xnv8+ChPtz3yvlkWnwunR32I56Q0nrPXnNJE/rcIRy8kL4E9ni/F1aDFm3VlsyHmRUhBDqaHGlB2x8P9cMTiuHt21uTt5n0WrRWm3F+LP/Ewu4RP471eeZ0H8iIX8605SKEQK0yjGQhpiRNV63NiBFv+snjXGXaKDe1f2koYyCMP5y7p81q1CFQ4pH/DxzqQ7XVgPdfvRzfevtmHO1z49P3Hcjb21oMSZLwH4+cxI92dRV1/2g8kdzzuBgGrRy64IZZr8WqjPUWlyl9bft6svva1J+9Rufkz57TrEez05SWafuH3x/FRDCK9+5YlvNrVyu72sY4QZKoovnDMRh1GuiVicsWvfx+4F/AoC2ekHBhPIhW5VywUQnaBpSx/+p8gsxMW8UOIhFCWIUQdvXPAG4BcGS2j0vl71ifHJh4QrHk1YhM3cpUnlxT/Wbqnx48hq89djLr9pODygj/hsmg7eaOBug0ImvR9vF+D9Y02qHRlL68raPJgZODXsSVk8R79/RgRY51Aa1VluSo2UzuQBTuYBThWAJnhvMHxbM17AkV7PeqthmgEZMvVOWs1m7EqD+9NE897nKtNKDSq1eC48GMksR8mTabUQd/JFayTLw/HMOTxwfx2ksaodNqcGNHA/7u9evw2LFB/PsjJ2b8uMPeMHzhGDpzTD7M5eywP7kqpNRXZGcTfE4lFI3jK4+exO8P9KbdfujCBDa0OKDNeL3c1OqCViNyDiPJN8a/vdGOE8r38cFDffjt/l58/PpV2NxWlfM5qVNsx7irjaiiecOxZGkkMFkeGVzAC1uDnhAi8UTyAr5Jr0V1yoLtCeXc1mVO39MWruBMWwOAXUKIgwD2APijJEmPlOBxqcwd65/MXp0byx1YqFm4kRIFbROBCH72YjfufqEb4Yya4pMDPui1AstShkM4LXpcubIGjxwZwIhPPqGXJAnH+z1Y2zQ3AwrWNtoRiibQPepH56AXe8+N4+3b27L6n1qrLbgwHsx5EpYazB3pdWd9vBT294yjzx3CplZX3vtUW42wGnQV0btVazMiGpfSJncO+8JwmHQLOuzlYqL2L2VepFGvpGYNIjHokJCAYIneAJ84PohQNIE3bpwcmPGeHcvx9u1t+N5zZ7MWpBfr7Ij8OjbgCeW9QJVKLb8GSlce2TMawBu+tQt/fveeOSk3Pz3kw5u+8wK+/fRpfOq+A/jly+cAyFnDo32erNJIQP7/29DswN4ca04G3OrC7PSgbU2DHWeGfLgwHsAX7z+Cja0ufPyGVXmfVzJoY6aNqKL5QrG0aguDTgO9VsxJpi2RkPDIkYHkxfN81MmRqVVXjQ7TZHmkcj7hsqiDSCo80yZJ0llJkjYqv9ZLkvQvpXhiVP6O9XmgXng9N5qdMYrGE7gwLp8kDedYejwTjxwZQDQuwReO4YXTI2kf6xz0YmWdLZl6V73ukiacGw1g6z8/gfVfehQ3fe1ZeEKx5KTHUlMf93i/B/e83AODVpNzulprlRmRWCLn90Z9IQGAI31zE7T9+IVu2I26vJPfAKDJYUJVhWSpam3ZC7ZHfGGWRs6jekfuTJsvHINOI5JveCqb0tPgL1GJ5AMH+9HoMGHr0vSszTsvbwMAvNw1OqPHPZuS7e4cmjrbdiylPLoU/RrPnBzCG769C8f7Pdh1egRPnRia9WOm+r9XL+AN39qFIW8Y3/uzy3Dj2np88f4j+PlL53Bq0IdwLJE2hCTVZUurcfD8BCIZJzHqYu3GzExbgx2ReALv+fEriMQS+Madm7Jes1MlgzZ/6ao1iGj++cKx5GJtlcWgm5OetkePDuDDv3gVTx4fLHg/dUfb0prJoK3JORm0TQTli0XOjExbRfe00cXpWL8nOdY+NchQ9Y4Hk1c5SpVpe+BQH9qqLbAbdXjkyEDax04OeHPuUnvLZUtw97u34u9fvw5v29aGFXU2bF9Wjeva56a3cnWDDVqNwP6eCfx23wXctqExbVG1aolyZed8ju+d+v1sb7DhaK8n6+Oz1e8O4qHD/XjrttaC/Wqfu6Ud3/uzy0r+9edCXXLB9uQV+WEvg7b5ZDfqYNJrsna1BcIxWI3ZGVu1XLIUw0jcgSie7RzC6y9tyip77mhywGHS4eWz+ScdFtI1Iu9/BJBzOXSmE/3eZFZxNlnERELCt586hff85BU0OU145NM7sbTGgv989GRJyiR94Rg++6sD+PyvD2JjqxMPffIa3Lq+Ed991xbc1FGPv/vdEfzLQ8cAZA8hUW1bVoVwLJF1cWnIE4JBp0nuN1Kpr9Gnhnz4u9evK7jUHkgN2qbOcBJR+crMtAFy9UWph1EBwOPH5GBtX89Ewfv1jAWg1Yi01SSNTlNWT1tyT5tO3dO2MJm28p8uQGUpFI3jzLAft65vxOkhf/JqRarulAElmUuPZ2LIG8KLZ0bx8etXoWcsgMePDSIWT0Cn1cAbiqJ3Ioh3NLRlfZ5Oq8ENaxtm/fWLZdRpsbLOil++fA6haAJv3579nAC5pw2QXzS2ZuyR6xkLoNpqwBUravDbfb1IJKSS9t/9/MVzSEgS3n3VsoL3q3eY0kZ2l7NaJThLzVwOe8O4JM/JJpWeEAINDlPWIBJfOJ7z4oC6v6wUY/8fPSZn4d+wMXv3olYjsH15DV46O/NM25oGOy6MB7PG1edyYsCDjUtcePHs6KyuIv/N/Ydx3yvncfumZvzbn1wKs0GLz97cjk/ddwAPHu7HG3P8W4t1tM+NT9yzH92jfnz6ptX4xA2rkz1rRp0W333nZfjoL/fhieODsJt0WFaTe2iT+tq1p2sMW1L60gaUMf6ZgfqqehuMOg2uXlWbNpwpH4dJD61GMNNGVOF84VhyD5rKYtAiUOKsVSyewFMn5WqEA+ezS7dT9YwF0OIyQ5eS7W92mTEeiCIUjcMdiMKg0yQzbMmetgof+U8XmU5l0Ma6JgeW1lhy9rSpJZMGnaYkg0gePjyAhAS8YWMzbtvQiPFAFHu6xpTnI4+QTh1CspA6mhwIRRNYUWvFFSuqc95nSZU84jrXBMnzYwG0VluwvtkBXziGczmycTMVjMRxz54e3LyuITkxaTFQB6qkZnWHveFk2STNjwa7Cf3u9GPaH47lXNCuBnKl6Pt64GAfltZY8pbxXbGiGt2jgeT+nenoGvFjRZ0V7Q22KTNtY/4IBj1hbFnqAjDzf1soGsdv9l3AW7cuwTfu3JRs2n/Dpc1Y22jH1x47iWjGZNlRXxjPnyq8lkCSJPzsxW7c8d3d8Edi+OX7r8Cnb2rPGjJi0Gnw3XduwVu3LsHbtrXm7WutsxuxotaKV7rSs5gD7lBWaSQgn/T88ZNX4zvv3FJUr6xGI1Bl0TPTRvOuc9A7p4N/Lja+cK5Mmw6BEmfa9p4bx0QgihaXGYcvuAv2tZ0bC6SVRgKTJd0D7hAmAtHkEBJg4TNtDNpoRtTJkeuaHWirtuTMtJ0bDcBi0GJlnS1rf9ZMPHCwD2sb7VjdYMfO9jqY9Bo8clQukVSvfucqj1wIal9brgEkKpNei3q7MecEyZ6xAFqrzFjfLJ+AlnIYye8O9GIiEMV7dywv2WOWA5dZviKvHmuBSAz+SJzlkfNsbZMdRzNWVfgjsWRWLZU12dM2uzftEV8YL5wewRsubc7783aFMr11un1t0XgCPWMBLK+1Yk2jHZ2D3oKDQE4o/Wxq1mmm5ZFH+zyIxiXc2NGQ9m/SaAQ+d8sadI8G8JtXLyRv398zjtf91y782Y/2YPeZkVwPiUgsgY/+ch/+/vdHcdXKGjz0yWtw5cqanPcF5MDtP96yEV983bqCz3X78mq80j2WdoI75A1nDSFRraq3T2s4ULXVwEwbzau7d3Xhlq8/h1/u6Vnop1JxDl2YwJ3fezFr3UnunjZtyXdZPn5sEAatBh+5biX8kThOFehDVi+Qp1JLJfvcQbiD0bQSb51WA51GMNNGleVYvwd2ow6tVRa0VVvQ7wllHcTnRv1oq7ag3m6c9SCS3okg9p4bT5Y+WQw6XNdej0eODCCRkHByQN671uLKvaB1vt28rgHXr6nDn27NP+QDkCdIZva0xeIJ9I4H0VZtQXuDHXqtKNkwEkmScPeuLqxrciT7ERcLjUagxmpIBm3qzra6AisNqPQuW1qFQCSeHOsO5L7CCkz2tE2nPFKSJPx0dze+/MAxfO/ZM7h//wV879kzySx8Ph1NDthNummXSPaMBRBLSFhRa0N7gx3jgWjB17Pjyr97Y6sLQsx8EMl+ZffZ5jZX1sdu6qjHplYXvvnkKYSicdzzcg/u/N5L0OsEGhxG/McjJ3MGlj/a1YWHjwzgL29bg7vv2oaaEv1sbFtWDU8ohpPKxTNJkvJm2maiymLAODNtNE9+t78XX35Q7uV8aopBFpTtkSMDeLlrDEf70vvx5Z629B7XUgdtkiThieODuGpVDXasqgUAHMjT1+YNRTHmj2Tt653c1RbCRDCSHEKiMum1zLRRZTnW50FHkwMajcDSGgskCclJkaruUT+W1ViVpcezC9r+eKgPgFwapLptQyOGvGHsPz+BkwNerG6Ym71rM7GyzoYfv2d7ckxsPm3K2P9U/e4QYgkJbdUWGHQarGm0JzObs7Xr9AhODfnw3quXV8QY/+mqtRmT/ZPDPrkMjpm2+aVmmFIXLgfC8ZzlkdZkeWRxQVs8IeFv7j+ML/3hKO7Zcw7/+vAJfOZXB/GD57uwttFeMNOu1QhsX1Y97WEkXcrkyOV11mT5daESyeP9HtTZjai1GWHRz/yEZH/PBJZUmVFvzw58hBD4y1vXoN8dwh3f3Y2/uf8wrlhZgwc+fjU+c1M7DpyfwGPH0k82+91BfOupU7ipowEfvW5VSV8r1QtAr3TL31tPKIZgNJ430zZdNTYDRplpo3nwzMkhfP7XB3HFimq8fXsrXjo7ljerEorGs0qUCclgLfV1MhyLIxJPpO1pAwCLUVfU6/9XHzuJT9+3f8r7nRry4dxoADd1NGBZjQVOsx4Hzk/kvK868G1pVqZNvvjf7w7BHYzBaU4/jzPqNJweSZUjkZD3nK1rlksA1Xrg1BLJeELC+bEgltZaUGeXT6Rns1vogYP92NjqQltK7fENHfXQawUeOdKPzkEv1pZJaeR0tFaZ0e8Opr3wn8/YG7K+yYkjve6S7Gb68QvdqLUZ8IaNTbN+rHJUazcmM21qH2Wh5eFUenKgYcSrKbu7fMr0yEy25CCSqd8Ao/EEPnXffty75zw+fv0qHP/ybTjyj7fiyc9di3vefznufve2KR/jihU1ODviz1pJIEkSnjoxmPPk7OyI3C+7otaK9sapg7YTA57ka5HZoJ1xeeS+nvG8C6cB4KpVtdixqgbH+z34+PWr8ON3b4PLYsBbLluCFbVWfOXRk2m9HP/voROIJST8/esLlzrOxJIqMxodpmSPcXKxtrN0mbZRP/e00dza3zOOj/xiH9Y02vGDP9+KmzoaEIzGcy6PlyQJb/3ei/ji/YcX4JmWN3XlyamUoU3qWhdr5q7OIi5sJRIS7t1zHi93TX3BTZ0aeZNSVr6x1ZU3aFPPtTLLI80GLVwWPQbcIbgDkawJuCa9tnL3tNHF59xYAP5IHOuUvq22anlk87mUaZH97iAi8YSSaTMgEk/AE5xZ30rXiB+He914w6XpgYbDpMeOVbX47b5ejPojaC+TISTTsaTagoSEtKW/PRkvJBtaHBgPRNE3gwEKqbpG/HjqxBDeeflSGHWLc9l0rc2QzOqqQVs9M23zSgiBy5ZWpQVt/kgM1hw9bZYie9pC0Tg+9PNX8eChfvz1a9bi87eugRACNqMOK+tsuGpVLZqLKI2+XBkKlFki+cChfrz3J3vxm1d7sz6na8SPaqsBLosBtTYjaqyGvBMkY/EEOgd9yZ5Ws0GLYJFZxFT97iD63SFsyVEamerbb9+CBz5+NT5/65rkIBGdVoPP37oGp4Z8uH+//O958cwoHjjYhw9fuzLtwlepCCGwfXk19nSNQZKkZNBWqvLIJVUWTASi8IRYIklzo28iiPf85BXUO4z4yXu2w27S44oVNdBrBZ7LMdzneL8Xhy64iwokLiZD3lDyvfdkyuukLyS/DtpM6QGQ1aibMmg71u/BiC+MicDUP/+PHxvEpUucyRLHTa0udA56c77HJBdr53hNlBdsBzERjGaVRzLTRhUldQgJIJ8oWwzatAmH51IWFtblGMVeyIA7hPNjAfjCMUiShAcO9kEI4PWXZvervGZDY/IKbLkMIZkOdex/6gTJnrEAdCl7Q9a3lGYYyc9fPAe9VuCdV+ReQbAY1Nkms7rDvgiEQM4deTS3LltahQvjweTJuz9Ppk2v1cCg08BfILCRJAkf/PmrePrkEP7ljg340LUrZ/y81jU5YDfq0k60IrEEvvrYSQDACzkGeJwZ9mNFyi6xNY12nFSm1WbqGvEjEkskM20WvW5GmTa1B6NQpg0AqqwGXJJjWuZrNjTikhYnvv54J/zhGL70hyNYUmXGR6+b+fduKtuWV2PIG0bP2OSEzlIFbeout+6R7CnFNDuSJLHED/LgEW8ohp+8Z3vynMVq1OGypVV4rjP7deH3B+ULIudGA/Au8osJ7mAUZ4d9mAhEppymqZZGLq+1Jqd6A4A3LH+PMnubzQbtlOWRT5+Qx/cHo/GCwdKQN4QD5ydwc8fkiqfNrS4kJODQhezzp56xAKosejgyAklAHkbSMxZAIBJPmx4JAEZm2qiSHOt3Q6cRWFVvAyBfZW3LGKih7mhbVmNNWXo8ddA25A1h538+jWv+42ls+NKjaP/bh/Gtp05h27Lq5JWTVDd1NEBtzajETFtrtTL2P2WCZM9YAC1Vk3tDOhod0AhkNfVOhz8cw69fPY/XbGjK2SOzWNTajHJWNxTDsDeMGqshbf8KzY/Llip9befGEY7FEY1LsOXoaQPkN/FCmTZPMIbnOofx0etW4p2XL53V89JpNdi2vDot03bfKz04NyqPfX7pzGhWGXLXiD9tAXR7gx2n8owCV4eQqJk20wyb7Pf1jMOg0ySrGaZLCIG/um2tvLvyBy+hc1BeZD2diY3TtT1lX5sarNc7SpPlVr//XQzaSkqSJPzVbw7h+q88U3As+mLnD8fwq73n8ZoNjVnL3ne21+F4vwdD3slKl0RCwoMH++FQ+rOO90+9u7HY53Hz156d8T7JufL277+EG776LDZ9+XGs+uJD2PJPj+Nvf5e7LFS9qH/7pmaM+SPJdgU105bZ02Y1aBGNS4gUCIKe6ZzMdHqC+QPkJ4/Lwd1N6yaDto2tLgDIWSJ5bjSQNYRE1eQy46zSz+zMKo8sTaZtxBfGvz9yYlqfw7MZmrZjfR6sqrelnQC0VVuS2TVA/mEw6DRodJiSS4+LGfv/6JEBRGIJfPG1Hfib167F+65egbdctgR/eeuanPevsRlx+fIa1NoMFTlwoslphk4j0gLe82PpLyRmZW3C0Vlk2n53oBfeUAx3XTW7k95yV2uXs2ojvrCyo63yjonFYH2zEwadBq+eG0dA7WXIkWkD5Olh/gI9bSPKAIrV9aW5KHP58mqcHfZjyBuCPxzDfz15Cpcvr8bHrl+FUX8k/epwKIphbxgr6mzJ29Y02hGIxNE7kb1f8US/BzqNwErl/ha9dkZv7vt7JnBJi/w9nKmrV8s9bwcvuLGzvQ63pJzIzIXV9TY4zXrs6RrDgCcEl0VfsiBxaY0FQgDdI6XbV7kYSJKEp08OzfgE8scvdON/917AhfFgclrpxei3+y7AG4rhPTnW4OxcXQcAeD4l2/Zqzzh6J4L42PWrAADHSjTd+eywH6eGfHiyjCZWnh324Vi/B3dubcXfvX4dPnb9KiyvteJXr5zPmSE71udBW7UFW5fKF3E6lQtZajVFdqZN/nu+KbsTgQj294xjZZ0cTE8UCNoePzaIJVXmtPkG1VYDltZYci7ZzjXuX9XkMCGmXMjIVR4ZnsX0yAF3CF9+4Biu/ven8D/PnpnW5zJoo2k71u/JugK8tMaCnrFA8upz94gfS6st0GjEtDJtDx0ewKp6Gz6wcwU+uHMlvvCatfjXP7kUW5flH0//z3dswLffsWUW/6KFo9UINLvMOD+eXh6Z+UKyocU547H/kiThZ7vPYX2zIznZb7FKXbA97AtXZCC/GBh0Gmxc4sSrPePJcf75grapMm2jyjTQUgXgyX1tZ8fww+e7MOKL4AuvWYsrldtfTCmRVDM7mZk2IPcwkhMDXqyqtyWDLfMMMm2RWAKHet1T9rMV44uvXYfty6vxj29cP+fTYjUagW3L5H1tg55wyUojAbnxv9lpRtdI7rLUUhvxhfHxe/bhlq8/i6MlOiGfC3e/0I33/PgV/OKlc9P+3N1nRvAvDx3Hte110GkEnlCyFItVvzuIZzuze9MSCQk/2d2NS5c4c/7MrWtyoNZmSOtr+/2BXpj0GrzriqWosRpmVQWTSr0QVKrHKwV1sMcnblyF9129HJ+7ZQ0+deNqRONScvBQqqN9bqxvdqC9Ub5wpfa1eUO53wfUwSSBaO73gOdOjSAhAXdsbgGAvH1tgUgMu06PJAeQpNqUYxjJoQsT6B4NJNt8MqVWdmVOAZcHkczsQsk3nziFnf/xNH76Yjded0kzHv/MtdP6fAZtNC0jvjAGPeGsA72t2oJwLIEhJTCTy43kEx2nWQ9dytLjQo/9ctcoXruhcVrPaWWdLXkiVolaq83JTJs3FMV4IJqVsl/f7MCgJ1xU4JtpT9cYTg56cdeVyxblmP9UyaDNF8GIN8wdbQtoy9IqHOl1Y0zpOc01iASQ38QL9bSprxs1ttL0Jq5vdsBm1OGhw/34/nNn8JoNjdjcVoXWagtaq814MaU0SS2PUa/yAkB7Q/rJSKrj/Z5kaSQws+mRx/s9iMQSU/azFWNdswP/+6Ers0q+5sr25VXoHg3gaK+7ZOP+VctqLegazZ9pOz1UeOl5MdQe6pu/9iweOzqIMX8Uf/Ld3fj13vOzetxiuQNRHM7Re5PL3u4x/OtDxwEAT52YXsB1YTyAj9+zH8tqLPj2Ozbj8hXVZZXdmQtffawTd929B48cGUi7/fnTIzgz7Md7duR+f9RoBK5eVYvnT40gkZD7//54qB83r2uE1ajDumZHyYIsdSBZqaZFl8JjxwaxvtmBJVWT5yTbllXDoNVg95n0Mk5vKCoHQk0O1NmMcFn0ycoF9eJdZnmk2aAOo8r9OvnMySFUWfTY2S5nPCcCuafI7j49ikgsgZtzVBRsanVh0BNGv1v+/kqShH9+8DhqbQb82RW5q4/Usf9A7kzbTPa0haJxfPPJTlyxsgbPfP46fPWtG5NtRsVi0EbTcrw/fQiJqk0J0HrGApAkCefG/FimTOTRaARqbIYpA45Hjw4gIQGvvXRxjqPPp7XKggtKT5s6kCQ7aJOHDczkqu/PXjwHl0WPN27Kv3h4sZgM2phpW2iXtVUhGpeS/Rm59rQBU5dHjvpKu7pBp9Vg67IqPHxkAKFYAp9PKb2+ckUNXjo7lqwYODvih0akTxezm/RocZmzJkhOBCLod4fSSnPMeu20l2ur++0qMSu+TamI6HOH0FCifjbV8loruoZ9OU9md50awU1few4f/eW+nIvaz48F8On79uM9P96DT923H3/3uyP4yqMn8f3nzuC+PT14+HA/nuscxod/8So+ce9+tNVY8cdPXo1HPn0NtrRV4S/+7xD++reH53xi3Kd+tR9v+u4LaZNXcxnxhfGxe/ahpcqMd1zehj1dY0UPwwhF4/jwL15FNJbA9/98K+wmPW5c26Dst1qcPYOSJOE5Jcv2+V8fxOmhyYztT17oQq3NiNdekv+8Y2d7Hcb8ERzr92DX6RGMB6K4faP8frqu2YFTQ96CPVnFUoM2TyiWtb91IQx7w9jXM45b1qVfSDcbtNiy1IVdp9IHtKi9fetbHBBCoL3BnnydTE6PzMq05S+PTCTk/7drVtclB4rly7Spg/DW58icbVL72pQBT48eHcCe7jF85uZ22HMMIQGAJldKpi3Hcu2ZZNrODvuRkIA/vWxJ3rLMqTBoo2lJTo7MLI9UDsBzo34MecMIRRNYmnJ1ty5lf1Y+Dx3ux4rayQW2F4vWagtGfBEEIrHJEbQZP9BqkDzdK3oD7hAeOTqAO7e2zukQgnJRbTVAI+Q6/EgswZ62BbRFGUbyvPLGnvlmrZqqPFKdAlplyf3mOhNqZv7Oba3J/jMAuHJlDdzBaHLP0NlhH5ZUWbJWZLQ32LLKI08of1+b8tpomUF55P6eCTQ5TTkHL5W7DS1OmJXXmVKWRwLyUCtPKIbxHCdtapDz6NEBvOk7L+DMsHxSnkhI+PlL53DrN57D48cGMeqP4MD5CTx4qA///ewZ/L+HTuALvz2Mj/xyH/787j14+uQwvvCatfjNh6/E6gY7am1G/Px92/GR61bi3j09eOv3XkxeYCuGJElFZ0wOnp/AMyflwOKz/3sg789EPCHhU/ftx0Qgiu++cwtu39iMWELKOoHO558ePIYjvR58422bksf+Tcq0vVwlkpIkYW/3GA5dmEib6lxJTg56MeQN4zM3tcOg0+DDv3gVvnAMZ4d9ePrkMN51RVvBNTjXKH1tz3YO4w8H+uA0T2Z+1jU5EI1LODU0+2Ek/e5QcnVHsRdonz45hAcP9c36a+fy5PFBSBJwy/rs7NWOlbU41u9JVlIAk7196kXmNQ12dA7IGXB/OAYh5NfEVOrfc1VbHOlzY8QXwfVr65IlihPB3Jm2cX8EWo3IOQlyXbMDBq0GB85PIBJL4F8fPoH2Bhvu3Nqa99+e+vqVmWkz6bQzyrSdVl6XVjdML7uWKve7KFEex/o9aHGZs2p8W6rM0GoEesYCybHMy1KuTtfajAVH/o/6wnjp7Bg+cu3KRV/Cl2lJlZyGvzAezLvs0WnWY2mNBft7JuAPx2AxaIv6Pt3z8jkkJAnvylMCsNhoNQLVVmNyih8zbQun1mbEshoL9nTLfQ/5etqsU/a0hVFlKe0U0Ndd0oS93WP49E2r026/ckUtAHmP24YWZ9bkSFV7ox0vnB5FNJ6AXnleahVCR2amrUB2RpKkrJ/j/efHKzLLBsgrHLYsdeGF06MlW6ytSp0gmbnG42ifGyvqrPjn2zfg4/fux+3ffgFffF0Hfn+gFy+dHcM1q2vxr39ySVqJlyRJ8Efi8ATl/W/uQBRLqi1oydj3p9Nq8Fe3rcWmVhc+/78H8fpv7cJ/vW1z8qQ9n1A0juu/8gzef80KvO/q7AEXmb711Gk4zXp8485NeO9PX8E/PXgM//bmS7Pu9/XHO/HC6VH8x5svxfpmJ2LxBBwmHZ46MYTXFMgWAfJFiHv39ODdVy3DjSlj0dtqLGhvsOHJ44NZz/U3+3rx+V8fTLvNrNfi7ndvw5UrK6MtQR0i8tZtS7BtWRXe9aOX8Zf/dxB1NiP0WoF3XF54DU6d3Yh1TQ48dmwQpwa9uH1Tc7JvVQ1QjvV5kn+eqd6JILa0ubCvZwJH+zy4bUPh/88njw/igz9/FS6zHq+7pGlG507H+z341H378Re3rs0qLXwsx2AP1VWravHVxzux+8xIch3T0T4Pam2G5G7U9gYbvOEY+t0heMMx2Iy6rOdoMebPtD1zchhCyMNgrAYtdBqRN9M26o+gyqKHRpP9PTDqtOhodmD/+Qn87MVunBsN4Cfv2VbwPcVq1MFh0sETisGRNfJfg9AMMm2nB73QCMyqXJ1BG03L0b70ng2VXqtBs8uEc6MBtFYpO9qqUzJtNiNOFBiL+/ixQcQTEl5zyfT62RYDNUDrGQ2gZywAp1mfdWUHAC5pceLBQ/1Y/6VHIYRcVmA1amE16mAz6pS/62BLue3/Xr2AG9fWzzgVX4lqbYbkCTSDtoW1ZWkVuvfJ+4zyZdqsBi38BbJRI74wakvUz6Zqrbbgh3dty7q90WnCilorXjwzivddvRxdI35sX549BGlNgx2ReALnRv1YVS+vAPjBc2fR4jKnHXNmgxaRWALxhJS8gq4a90dw6zeew0euW5mcWjfsDeP8WBB3XbmspP/e+bRtWbUctJV4tUjqrjZ1pYTqaJ8HW5ZW4apVtXjgE1fjI794FX/928OwG3X4tz+5BHdua806WVSXs9uMOjRj6sXst65vRPsn7PjIL17FXT/eg8/c1I6PX78q50kiII8f73eH8MuXz+G9efqlVEd63Xji+CA+e3M7rl9bjw9fuxL//cwZXL+2Hreul98T/eEYvvnkKXz/ubN469YleOs2OUug02qws70OT58cRiIh5X0+APD1J07BqNMmpx6murGjAT947izcKcuEw7E4vv54Jza0OPCpG9sxHohg3B/B1x7vxJPHB8sqaPvp7m40OEy4LUdP/HOnhrG63oYmpxlNTjP+8ra1+LeHT0AI4I5NLUWtwbmmvRbfe/YsAOCNG1uSty+vtcKs1+Jonwd/Ost/Q99EENevqYcnGJtyL+tLZ0fx0V/ug1Gnwag/gp6xyTkCxeoe8ePPfrQHI74w/u53R7BjVQ0sSrmiPywP9njX5UtzHrsblzhhM+rwwunRtKCto8mRvL86tKlz0AtfKJbzPaBQpu3pk0O4tMWJGqVixmXR550eOeYPF9zJurnVhV+9ch4n+j3Y2V6H69bU572vSu5rC2a9dpv02hlNjzw97ENbdXblxnSwPLKMlKImei4FInI5Qb5pO23V8gTJ7lE/dBqB5pSa4Fq7EaP+cN7FjH883I+lNZYZ7yWqZMkF2+Ny0JZvb8jfvLYD//ymDfjr16zFJ65fhbdubcV17fXoaHSgymJANJ7AhfEAXu0Zx8NHBvDj3d1wB6N4/zUr5vOfs+Dq7MbkpCoGbQtLHfsMZJfFqNRMW76Sq1FfBDXW+ft/vGJlDfZ0jaF3IohAJJ62WFs1OUHShz1dY3jzf+9GNCHh+39+WdoJjvpvzpVtuzAexJA3jH984Bj++xl57LM6dr0UQ0gWys3rGlBnN6Ijz/vETLVWW6DViKxdbeP+CHongslelhaXGf/7oSvxL3dswKOf2Ym3bW8rWfXG8lor7v/oDrxpUwu+9ngn3v+zvXn73O7fL1+sODvsn7Ks/dtPnYbdpMNdVy0DAHzmpnasb3bgr397GEPeEB450o+bvvYsvv/cWdy5tRVfvn1D2uffsLYeI75wwQnDx/s9eOBgH969Y1nO18WbOhoQS0hpExbv23MevRNB/NVtchbmrVtb8aFrV+KSFmey97IcDHlC+KcHj+Ef/nAUsYxF4aFoHHu6xpIljgDwoZ0r8JoNjZAk4N07lhX1Na5VPr/RYUq7kKPVCKxtsidLqmcqHItjyBtGs8uM9VMMNzl0YQLv/+letFVb8MO7tgLAtP8/+t1BvPOHLyOeSODf33wJBjwhfPfpyfHzz3UOIxJL5CyNBOSLBVesqMZuZdpuJJbAqSFvWrYxLWgLFw7aMsvIx5VS5tTgymnWw50n0zbuj6LKkj9o29TqQjAahy8cwxdf25H3fqmaXKasHW2AMogkFp92mfDpIR9WzXJ1TVll2iRJwr/88TieOzWMRz61s+AVo0J+u+8C1jc7sSZHSrdcPds5jA///FU88xfXlXzqVqkc7/ciIckZn1zaqq149OgAWlxmtFZb0lLPdTYjonEJ7mAUVRlXQ8b9Eew+M4oP7lxx0ZVGAnJmyKzX4vyYXB6ZK5MJAM0u87TLHKe68roYpfaxsadtYaVmRAqVR8YSEsKxRM6+y1F/BBvyvObMhStX1OCel3vwwMF+AEjb0aZaVW+DRgA/fbEbB85PYEmVGT99z/asjLba3xWIZJ+wqAMz1jTY8e+PnEA4FkcomoBeK3I201eK9c1OvPLFm0r+uHqtBq1VZnRlDMtQT5ZTv2cmvXbWi9jzMRu0+NpbN2LjEif+4YFjuHdPT9Z+rzF/BM+cHMJbty7B/ft78bv9vXmP4ZMDXjxydACfvGFVMsNl0GnwjTs34fXf2oXbvvE8xvwRrG2049vv2IzLlmZnfq9tr4MQ8hTJS5e4cn6drz/eCbtRhw/tzH0Rb1OrCzVWA548Pog3bmxGIBLDt546jStWVOPqVbVp992ytAo/eaEb4Vi8qKzBd54+jWFvGP/wxvVT3ncm7nvlPGIJCQOeEJ4+OZxW5renawzhWAI72yf/DUIIfONtm/DJYX/e99tMly2rQrXVgDdf1pKVeVnf7MDv9/fN6v120C23jzS7TLCZdPjt/l4MeUKozzgfPDXoxV1374HLosfP33c56uxGWA1a7Ds3gTs2Lynqa435I/izH+2BOxjFvR+4ApcsceKls2P4/nNn8adbl2BpjRWPHRtElUWPrUvzX0DasaoWTxwfwvmxANzBKKJxKe3nsMoql0qeHPDJQZspV9Am3xbIKJF/7tQwJAm4bs1ksO2yGPL2tI36wwXP+dWS8zu3tRUdG3zwmhXJieipTHotJAmIxiUYdMX9f8fiCXSN+HHD2tntyyybTJskSfjnPx7HD3d1oXPQhwFPaOpPyiEaT+Av/+8QvvFEZ4mf4dx66ewogtE49k0xNWohqY2xG1pyv8gtrbFgzB/BkT43ltakn7wUWrCtlka+dor67cVKCIHWajN6xvy4MB4saSnjxRawAUiW0uk0ImvqE82v1fU22I066DQCxjyLoq15rrSqRrxh1BQoeyk1dUjJPXvk3Ve5+g9Mei2W1Vqxp2sMG5od+M2Hr8r5c6sujg1Fsqso1KDt399yKd68ZQm+8cQp/OzFbqxrdl4UQ4NmYlmtNdkzrTqaMfxgPggh8O4dy7FtWRV++HwXohnZnT8e6kMsIeHdVy3HdWvq8cChPsTzVJl866lTsBq0eG9GL9nqBju+9Ib1SEgS/vZ1HXjwE1fnDNgAoMZmxKZWV97R/4cuTOCxY4N4/zUrsvrRVVqNwPVr6/H0iSFE4wn8ZHc3Rnxh/MWta7Iupm5pcyESTxQ1GEuSJPzsxW7c83LPtCepFiMWT+Cel3tw1coaNDiMuHdPT9rHnz81DINWg8uXp5dyGnXaogM29f5Pfe5afOam9qyPrWtywhue3cRHdUdbi5JpA7IHj0mShE/cux9ajQa/eN/laHSaoNUIbGx1FZ1pC0biuOvuPTg/FsAP79qKS5bIPzdfeM1a6LQC//TgcUTjCTx5fBA3rG0o2Pe1Qwnmd58ZSV48yazEam+w49RQEZm2jIz1MyeHUW01pF2EqLLo8/a0jQeiBcsj22os+MX7Lsffvb64LBsg9+29aXNL1u3qe9l0+trOjQUQjUvTHvGfqSyCNkmS8G+PnMCPdnXhciXtnDqSdTrOjQYQS0jYfWY074tkOVInkR2eoo55Lg24Q3jTd17IKj9RHel1o8ZqyDsVbHKCZADLMmqr1RPpXGP/HzrSj9Zqc95g8GLQWmXBq+fGEYkn8pZHUnHU7FqtzXhRBq3lRKMR2Ly0CtYcDegqNQOXaxhJKBqHNxyb1zLXOrsR7Q02nB8LwqzX5n29e/u2NrxtWyt++f4rsqoHVMlMW47Fsb6wfPLhNOvxn2+5FG/f3oZAJF6SpdqL1bIaK7pG/GllSUf7PGhymgqesM2VD+1cid6JIB463J92+/37e7GmwY6OJjtu39SMQU8YL6fs/1OdHvLhj4f7cddVy3IGU++4vA37/+5mvP+aFVMO4rlhTT0OXXBjyJt9wfurj3WiyqLHe69eVvAxbupogCcUw1MnhvA/z5zBjWvrcwaKavluMReZzwz7MegJIxJPYO+57GXM+RRbevbE8SEMeEJ491XL8NatrXjm5FAyAALk6bXbllcl94HNhivPQKTJIGvm52/quP9mlzllWnT64x3t8+DEgBefuXk1lqVcTNrSVoUTA14ECuy7VD12bACHe9342ls3pe23bXCY8IkbVuOJ44P4yqMn4QnF8pZGqlbX21BnN2LX6VEc6/PAYtBieca5nzr23xOM5gzajDoNNAIIZKx9eaV7DFeuqEnLajrNhpxBWzwhYTwQQXWB8kgAuHp1bTKzNxtG5XV9OitA1JhmdaUHbZIk4SuPncT3nj2LP7tiKb79ji0AZh60nVVGarqD0eR4+kowVdDmDUVx+f97Ao8eHcj58VJ48FAfDpyfyFo+qTrS68H6Fmfek6/UK82ZmTZ1mlDmBEl3IIoXTo/gtRtmNvlosWittiRHWbdWT90UT/mpQRv72crDh3auwMeuX5n342rQpvYhphpVxknPZ6YNkEskATmzky/w/8DOFfi3N19a8GQw2dOWI8PgU05SbEYdNBqB/3fHBvzX2zfnHBJBshV1VgQi8bSLf0f7PAtWTnrD2nqsrrfhf549mwwyzo36sa9nAndsaYEQAjd1NMBq0OL3B9LHskuShP989ARMOm3B6ZLFvi9ev1bu/VHXBqhe6R6T2y+uXZl3J5XqmtW1MGg1+Mv/OwRPKIbP3bIm5/0aHCa0uMzYf35iyuf1wmm550kjgBdOZweuuZwY8ODSf3isqOzRL146h2anCTesrced21ohAfjVK/Iy9EFPCCcGvNi5uvCkz9la02iHViNm1demBm2NThMcJj2W1VhwpDf98X63vxd6rcDrMqaEblnqQjwh4VARi9lfPDMKu0mXc2DLe69ehmU1FnzvubMw6TVTft+EENixsga7T4/gSK8bHU2OrNfLNY02hKIJdI8GcgZtQghYDbq0SgtPKIoL48GsrJ3Looc7xyCSiUAEkoR5u3BjUjJt0xlGosY0Kys9aPvZi+fwnafP4O3b2/CPb1yPWpsBTrM+uWclUyIh4fcHerOaTVVnU7JEu04Xt7dkoXlCUfROyBNqjvS6c15hevXcOAY9Yeyew3/TY0cHAQB7u7OvhoVjcXQOerGhwJtjaqCWnWlTgraMTNue7jFE41La+OGLkTr2H8je0UbTo5bilnriIM3MjlW1+ODO/EGbGlznyhCoi7Vr5rk3UZ2Kt6Ju5qOZASTLHHMGbRnLZoUQeOPGZvZhFqC+r6jVIMFIXBmONX+lkak0GoEP7lyB4/0ePKfsSfvd/j4IAbxRWb5s0mtx64ZGPHSkP+3K/L17zuPRo4P45I2rS3J8r292oMFhxNMpJZKnh7z42/uPoNZmxJ8XMZHUatQldxW+YWNz3qFjALC5zYX9RWTadp0eQVu1BZctrUoGcFN54fQovOEY/u3hEwUzbmeHfdh1egTvuLwNOq0GS6osuLa9Dr96pQexeCK5I/KaOQ7aTHotVtZZp71HNVWfO4hamyH5mrG+2Ymj/ZNBWDwh4Q8H+3DdmvqsrOzmViXzWUSQ+8KZEVyRkcFSGXVa/P0b1gGQv2fFZCd3rKrFqD+CV3vGc148Wa0MI4knpJw9bYDcJ5qaJexUkhgdTem9Zy6zHr5wLKsceTwgX9zLV/FQamqmbToLtk8P+dDsNOWdolysBQ/a/nioHxtaHPiXN22ARiMghMCqelveTNvzp0fwqfsO4Injgzk/fnbYlyxvUafalFogEsPt395VsgBKzbJd116H8UA0LbWvUpeHdg7OLAM5lRFfGHvPjUGvFdh7bjxrymPngA+xhFRwIIDdpE9e6cjMtDnNeui1AiO+9CbS/T3j0GkELl2yMG+65ULNUmqEXB5BM6cGa8y0VQZ1L1au1z21B3a+A/DLl9dApxFon+Wkr3yT0QC5PFKrETDpF/xtuGKk7moDgOMDHiQkLOjglts3taDBYcT3nj0DSZLwuwO9uGJ5Tdrr+Js2tcAbiiWzYCcGPPjHB47imtW1eQeDTJcQAtevqcfzp0bgD8fwjSc68dpv7sKAJ4T/eMslRZcHvmFjM0x6DT6TscMw05a2KvS5Qxhw558/EIsn8NKZUexYVYsdq2pxpM+NiUDuQRKp1HH3e7rGkoFXLr98uQd6rUiuPwCAd2xvw6AnjKdPDuP5U8OotRlz7hkrtfXNzqzqrm89eQp3fu/Foko9eydCacfM+hYHzo8Fk9MSXzwziiFvGG/alN1jVWU1YEWtFfvOTRT8GufHAjg/FsRVBVY13LC2AX//+nX41I2F//9Val+bJCHn9O/UckB7gWFUqa+RJ5Rz4jWN2Zk2AFnZtlGfWpExP+/5aqZtOgu2Tw15Z51lAxY4aJMkCZ1DXly6xJWWUl1VZ8ubaVN/mI/l2fl1ZtiPFbVWXLWyFnu6xqZVc1qsl86O4uAFNx7PEzgCcolEX46TkFzUA/Qtl8mTf3Lt53hFyX6dGsq/62w2njo+hIQE/NkVy+AORtGZ8XXUUcIbprii2VZtgUYgbYkpIL+h1NqMWYNI9vdMoKPJcdE33qtj/5td5uSyXpqZOpZHVpQGhwk6jUBvjiZ+9SLPfGefqqwG/O5jO/C+a6ZeilxIoZH/6t6ii7ksfLqaXWYYtJrkBEk1s7GQQZtBp8H7rl6O3WdG8YuXzqFrxI87MoYXXLWyBrU2I35/oBeBSAwf++U+OMx6fP3OTSXtu71+bT184Riu/8oz+MYTp3DbhkY8+blrpzWx7s1bWvDKF2/KOTU11ZalU2d3DvW64Q3HsGNVDXasqoUkycHHVA73urGzvQ4tLjP+89GTOYOeYCSOX+89j9s2NKXtWbthbT0aHEb84qVzeP7UCK5ZXTsvvc3rmhwY8ISS1QE/f+kcvvp4J17uGitqQEn/RBDNzpSgTTnXUrNtvzvQC7tRhxs7cu8X29TmwoHz4wUDRDWRsSNjGmim9169vOiJvc0uc/JiSq5hQHaTPnlhLm+mTZ+eaTsx4IHdpEOzM72f2KlkGDMD/8lM2/wMHptupi2RkHBmyD/rISTAAgdtw74wJgLRrMa8lfVWjPgiOa/IqI2ZJ/LUDp8d9mFFnQ1Xr6pFOJaYk10i6pWf4wXqlz/8i3345L37i3q8k8oBev3aeug0IqsuORJL4MD5CVgMWoz4IhjzT32lKpdXusfwD384mnNX2mPH5FH971b2xLzSlV4ieaTXDbtJN2W/VUeTHavqbTDkmBRXazOmlUfKNdgT2MzG++T3laWRs1drM+L2Tc24Ye3UyzNp4Wk1Ao1OU86LXMkrqAtQ6rqhxTnrUpaC5ZHh+Kwf/2Kj1Qi01ViSEySP9bnhNE+eFC6Ut29vg92kwz88cAwGnQa3XZLeL6TTavD6S5vw5Ikh/MWvD+HsiB/fvHNTyS9GXL2qFjajDnqtBj9+9zb819s3T/trCCGm7H0D5CDFoNMUHEbygnKudNXKWmxqdcFq0OKFKSqg/OEYzgz7sKXNhc/c3I7Dve6cvfwPHOyDJxTDuy5vS7tdp9Xgzq2teLZzGGP+SNqo/7mkXjg41u/Bo0cH8KXfH0lmnqbq/ZMkCX0TwfRMmzqMpNeDUDSOR44M4LYNjXkvcG9pq8KIL4LzY/kDxN1nRlFrM856GEamnatrYdRpsLoh9+O2K7fnX/uiTc+09XvR0ejIuqClToPOHEYy2ftcnpm2PncQwWgcq2dZuQEscNB2Sin1UxfwqdRoNFe2Tb2ydnwgO2Aa80cwHohiZZ0Vl6+ohlYjiq6hng71MY/3e3Ne1QhEYjg54MGrPePJqy6FnOj3Ym2jHSa9Fqsb7FnDSI72uRGKJnD7JrlGvnNw+tm2eELCF+8/jJ/s7sZDR9InXQUiMTx/agQ3r2tAa7UZDQ4j9nSnvxAf6fNgQ3P+ISSqL75uHX75/ityfqzOnp5pOzXkhT8Sx6ZW17T/PYuN3aRHo8NU8hfTi5FGI/DNt+XeZ0TlqdllzlseaTFoSzLxayEUzLSFc09To8KW11qT5ZHqEJKFzlbaTXq864qliCck3NzRAEeOoOdNm1sQiSXwx8P9+MQNq3HVFNmOmbAadXj8szvxxGevTQ4mmSsGnWbKJdu7To9gfbMD1VYD9FoNti+vxu4phpEc6/dAUvbB3rG5BavqbfjKY51p08BHfGH8aFcX2htsaYuuVW/d1gr1kJgqq1Qqav/fPS/34JP37selS1y470NXwKTXYP8UyQNPMAZ/JI5m12RmqdZmRKPDhKN9bjx5fAi+cCzn+HmVuocs3/+HJMlT1a9aWVPyn5fP3rIGv/nIVXkDynalPDXf653ZoINfCdokScLJAW/OXWpqeWRm0Dbmm99Mm2mambZTSrtXxWfa1OAjMzpfVSf/Z2X2tXlCUZwbDcBp1uP8WBDeUPp/nDo5ckWdFXaTHhuXOIueVlSsQU8InYM+LK2xwB2Moj9HPffxfrnOXpKyJzllkiQJJwcnD9BLWhxZw0j2KgHUO7bLy0JPzSBoe/BQHzoHfbAZdfjGE6fSXgCf6xxBOCZvvhdCYNuyarzSNZZ8DtF4Asf7PUWN5LcZdXnL0mpthrRM24GeCQCT44Mvdv/7oSvx2TzTuogWsyUuc87yyFFfeEGybKViLtjTlnvZLBW2vNaK7tEAIrEETgx4y2YR+Xt2LMOaBjvuUqpVMm1c4sT6Zgd2rKrBJ2+YuwmhTU5zScbbF2NLmwtH+jw5T14DkRj290ykLebesaoWZ0cKt44cViqNLmlxQqsR+NzN7Tg95MPv9vcikZDwi5fO4YavPIOzIz585qb2nAHIkioLbl3XiC1trrTSybnkshjQ4jLj4SMDaHaZcfe7t8Fh0uPSFhcOTJFpS93RlmpDiwNH+jy4f38vGhzGtBH9mdY02uUl23mCttNDPgx7w9ixKv9jzJTTrC9YTrlGSczY87zeWQ1aBJXyyN6JILzhGNY25QjazEp5ZEZP21ggAptRV9Si91Iw6qeXaTuzeII2H1wWfbIHRdVSZYZBp8GZ4fR9YWqTZ76M01nl/iuVWuwdq2px6MJEzhGhM7VLSfe/XxnRm2utgPqiYzPq8i67VPW5Q/CGYlirNFxe0uLMGkay99wYltZYsKHFAbtRN+1hJNF4Al9/vBMdTQ78659cgtNDPjx4aHL88GPHBuA067F9mXzFavvyagx4Qsk67DPDPkRiiaJrnPOpsxsx6o8kyzP390zAZZFH25K8/NHJZdB0EWqpMmPAE8qaCjbii1T0NEWTrlCmjeWRM7GsxopILIFdp4cRiSXmdal2IfV2Ex79zM6cmR9ALjv8zUeuwi/ed/mUO9cqxZa2KkRiiZznQa90y3tHd2QEbQAKVkAd6XWj3m5EvbIf8bYNjbikxYmvPd6JO/57N/72d0ewvtmJhz91DV6TMfo+1Tffvilv1c9c2dzmQq3NiJ++Z3tyKNvmNheO9uYObFWpO9pSrWt24uywD892DuGNG5tzTnxUqUu29ysXwzOp3/OrVs5P5jHVVStrsaXNhXVNuX9WzQYt/MoKlBPKvIq1jdkXY5zJTFt6i9CYPzKvexrV1/ViM22nh3yosRpK8hwXuDzSi/Z6e9aVEq1GYEWtNSvTppZGvnmLPLAjcxjJmREfDMrYV0B+gUhIyLnUcqZ2nR5BjdWA25U0da6+tiN9HtTaDHjdJU14rlN+Y8nnpFLmqU43UgMjdRiJJEnY2z2OrUurIYTA6gYbTk4z0/bbfRfQPRrA525ux+suacLaRju+8cQpxOIJxOIJPHl8CDeurU++kWxTgrc9Sl+buitktm+OtTZjcgkiABw4P4FNra4FL20hooXV4jIjIcmVDKlGfOF561OYCxplOmQwx9JbX4jlkTOhDj148JBc5l8umbZimPTaRfV+l1yynSNQeOH0CAxaTfJ8ApAzLrU2Q8Gg7XCvG5ekXCAWQuDzt65B70QQveNBfOPOTbjnA5dj1RT9QUaddt4yjqr/eMuleOKzO9GWciF6U6sLkXgCx/MMzwPknicAaHKlZwU3NDuQkIBoXCpYGqna3ObC8X5Pzh7a3WdG0VptTtunO18anSb89qM70OjMnfW0GnTJC1snlHPiXOWRDpMOWo3ISsSM+SPzNu4fmH6m7dSQrySTI4EFDNokSULnoDdv42Kusf9HlSswly5xwmHSZQ0jOTvsx9IaS/JqxOY2F0x6Tcn62iRJwq7TI7hqVS0cJj3aqi05e+uO9LqxocWJGzvq4Q3Hcu49U6k/yGrNb0eTA1qNSPa1dY34MeqPYNsy+cWxvcGOU4O5e+lyCcfi+K8nT2Njqws3dtRDoxH49E3t6Brx43cH+rCnewzuYBS3rJ+cLrWmwQ6HSYe959SgzS1vuq+d3c4i9Yr5iC8Cb0ieUKnuFyGii5d6hTmzRHLEF0GdvXLLIwHAknJCksoXjjFomwH1fejxo4Mw6TVTTjmkudPoNKHZacpZkrfr1AguW1qVFjhpNAJXrqzFC2dG884DODPsy6rquba9Dr/64BV48nPX4k2bW8o28LUYdNk71JTAtlBfW+9EEAatBrUZF6jWK9+H1fW2nOP0M21pq0JMGfCWKp6Q8NLZUexYgCxbMSwGLfxh+cLWiQEvWqvNeRdxO8367J42fwQ1C5BpK2Y6vSRJOD3kK0lpJLCAQduQNwxPKJY1hES1ss6G8+OBtG9KatPx2iZHclS+6sywL20ZqlGnxfblNXhhihGzp4e8uHtX15SB0MlBL4a9YVyjpPjXNTmyrp6EonGcGvJhQ7MTO1bVwqDT4Inj+UskTw540eIyJxuXTXotVtfbcFjJbqn9bFuVq1XtDXaMB6JZ+87yuW/PefROBPEXt6xJvtDdur4BG1oc+K8nT+Ghw/0w6jTY2T65fFKjEdi6rDqZaTva58Y6JZicDbXXbcQXxuELbkiSPKaWiC5uLVXZu9oSCQlj/srOtAHqOOs8I//Z0zZtDQ4jzHqt3PfSOPv3JZqdzUurspZsj/rCONbvydk/dfWqGgx7w8nhDKmO9cnzAC7J0Ypx+YqaimwfaHSa0OgwFexr65sIocllylpN0Ow04coVNfjANSuKClTzZT6P9LrhCcVwZYH9bAvJYtAhHEsgnpBwYsCbszRS5TLrs3va5rs8MjmIZOpM24gvAncwe0r+TC1Y0JZvCIlqVb0NkjS5RDMUjeN0yhWYjkY7Tg54k/1R0XgCPaOBZD+basfKGpwe8hVcAPmjXV348oPH8OjR/HvXgMl+tqtXy0FbR5MD3aP+tP0Sx/o9iCtLqK1GHa5cUYMnTwzmDQhPDnizFj9e0uJMDiN5pXsMVRY9VirBqBrkFjOMJBiJ49tPn8bly6vTXjyFEPjsze3oGQvgnpd7cM3q2qzpbNuWVePMsB/D3jCO9nlm3c8GTGbahr3h5AjcTUtcs35cIqpsLTkybRPBKBLS/C/WLjWzQZtVrpRISPBH2NM2E0IILEvuhaqc0sjFKteS7d3KhfJckxvVnqpcFVBqhdElS8qjT7FUNhXoNQPknrbUHW0qIQTu/eAVacvDC6m2GrC81pqV+VT/Pxain60Y6pTd8UAEZ4d9BZehOy36tJ42SZLmPWgzJkf+T51pU3crV3ymrTPPuH+V+g9USyRPDngRT0jJF+m1TQ74wrHksIzzYwHEElJWqYT6orG7wG6Qg+flF4p//uOxgv8Ju06PYEWdNVnK09FkhyQhLeN3NONF56aOepwbDeDsiD/r8SKxBM4M+7Jqdy9Z4sSYP4I+dwh7z43jMqWfDZjcd1HM2P/vP3cWw94wPpeSZVNdv6Yem1pdSEjAzeuyF29uXy5fsfn1q+cRiMRL8uaYmmnb3zOBFXXWZGMpEV28THotam2GZG8HgOR6kJoKHkQCyCckmeWRfuVCH4O2mVleK/fllMsQkovZFqVa5levnMezncN4tnMYfzjYB7tJlzNj1lptwdIaS96grc5uRINjfiY+zpfNbS70jAXyroDK3NE2q6/V6sJLZ0bxxLHJZMHuMyNob7Dlney90CxGOWg7dGECCSn3EBKVy6xP62kLROIIxxLzGrRpNAIGraaoTFspJ0cCCxi0nRr0otpqyDsZbHmtFUJM7mo7oizVVl+kO5T6XrWnTJ0cmVoeCcgljFUWPXbl6WsLRuI4OejF9uXVuDAexPefO5vzfuFYHC+fHUuWRqY+h9TJSYd73ai2GpKb3NVdKU/lKJE8M+xDLCFlBW1qVuuZk0PoGvEn+9kAOfBxmvXozFFakOpnL3bj60904vWXNuWcZiWEwN++rgOXtDhxy7rGrI9vaHHCoNPg5y+eS3tOs+Ew6WDQajDsDePA+XH2sxFRUrPLnLwIB6QGbZWdaTPlKI/0Kf0bLI+cmWU1zLSVi3XNDtiMOnz9iU7cdfce3HX3Hjx+bBA7V9flnZK5Y1Utdp8ZzZoCeCRjCMlioe6izVUiGYsnMOgJocVVmkD1Q9euRJ3DiPf/bC/e+cOXceD8BF7pHivbLBswmWnbd24CAHKO+1e5LIbkMDtALo0EgGrL/L5PGHWaIjNt8qqtxhJdiFiwd4zOQW/BGk+TXovWKksy03a0zwOHSYclSu9De4MNQsjjQW9d35gM7lbWpj+mRiNw+fIavHw29zCQo31uxBMSPnDNCtTaDPjuM6fx5suWZO3L2HduAsFoPC3dv6TKDLtJlzZB8nCvXEqoZraWVFmwttGOJ08M4gM7V6Q95kklQ9eR0WCq9o/9dHc3gMl+NkAOttobbAXLI3+6uxtf+sNR3LyuAV9766a899u6rBoPfOLqnB8z6rTY1OrCnq4xGHSaklwlEEKg1mbA/vMTGPFF2M9GREktLnPaZFy1bzdzJUylsRi0yRMLldp0z0zbzNy0rkHufSlwckfzw6jT4pFPX4NBT3oWKdf0P9WfX7kU9+7pwf88exZfeM1aAPIQktNDPty2If8Y/0p1yRJ559yB8xO4sSO9smnQG0ZCyh73P1NrGu149NM7ce+eHnz98U686TsvAJi/JeMzobbn7OsZh1GnSV6UySVzEEkyaJvHTBsAGPXaoqZHnlYmR5ZqeM4CZtp8eUsjVakTJOUhJJPBkMWgw7IaazJgOjvsR63NkLPc7vIV1eidCOLCeCDrY+qVj41LnPib13ZAkoD/99DxrPvtOj0MrUbgipXpvWEdjY7kcwhF4zg16MWGjKt/N6ytxyvd41ljSk8MeKHXiqypjOowks5BH4w6TdZS69UNdnQO+nL2yf34hS586Q9Hccu6BnznHVtg0M38v1jd29bRaIe+RHtl6uzG5DTNzcrVJyKiFpcZfRPB5Ova6CIpj8w1iMQbYtA2G1vaqnD3u7fN2zJdKmxJlQWXLa1K+1Xo2F7b6MDtG5vxk91dGFLWfBzvzz+EpNJZDDqsabDn7GvLt6NtNvRaDf78ymV45i+ux4d2rsCWNlfZDiEBJjNtB85PoL3BXnC4kMuihzcUQ0zZ6ZkM2ua5IsOk1xS1p+30kK9kQ0iABQraovEEvOFYsj8rn1X1NnSN+BGJJXCi35NVCtHRZE/udDg74sOK2tyPp5YHqtMQUx264EaT04R6hwlLqiz4yHUr8cdD/XgxY+LkrlMj2NTqSk55VK1rlqdYJhISTg54EUtIWS86N3bUI56Q8FzncNrtJwY8WFlnyxkQqeWIG5e4st6Y2uttcAejGPKmX9n6+Uvn8I8PHMOt6xvw7VkGbACwTfm+rS/hi2itzYiEJB/whZpNieji0lJlRiiaSL4Jj/oi0GoEXBU4MS5VrkEkLI+ki92nb2pHLC7h20+fBgAcvqDMA1iEQRsg97UdPD+RHJ6nmougTeU06/HXr+3Abz+6o6wvEKmZtkAkPuV5ofp+4FEufI0uYHlkeIpMmzckn6dnDkicjQUJ2tSU4lTLEVfWWRGOJfBc5zDCsQTWZ2Sc1jY6cG4sAH84hjPD/qx+ttT7OUy6nEHbwQsT2JgywfDD165Ei8uMv7n/MP79kRP490dO4F8fPo5DvW5cnSO93NFkRyASR89YIDn5KLP/a1NrFaosejx2LH06Za7JkSr1hWvrsuy+LzVDmTqMZNgbxr/88Riuba8rScAGAJctrUJrtRnXpawDmC21EfbSFlfeenciuvgkd7UpJzEjvjCqrYasMdiVxqzPHkTiY6aNLnLLaq1467ZW3LunB+fHAjjc60GtzYgGR2Vn1vPZ1OqCNxxLtvKoepNB2+IavjIdlpRdfoXKagEkl2ir/ZDjC5Zp006ZaVOr60q5Q25BzprVf2gxmTYA+N2BXgDZk6LWNsrTG/d0jWHMH8kbtGk1AttS9o6pxv0RnBsNYGNKmZ5Jr8U/37EBY/4IfvR8F370fBd+vKsbNqMOt23IHtiRHEbS78GRXjdcFn2y7y7167/u0iY8cLAPn7x3P8b8EbgDUfS7Q1ibZ2Hi9uXVEAJp+9NU6iJudQInAPz3M2cQjUv4hzeuL1kpo82ow/N/eQNuWZ/9754pdfAM+9mIKFXm2P8RXyTvoKpKYimUaWPQRhexT96wGkIIfOOJU8oQEkfZLs6ereSS7YxhJH0TQbgs+qy1SxeT1KAtc8ZDJnVXn7qrbdQfgV4rYJ/n11J5EEnhTJs6qCR1wfxsLchREoom0GI1TNmroKYUnzg+CKNOgxUZvV/qf+4fD/cDQN7ySEDua3vyxBCGPCHUK1NcDipb4ze2pgeD16+px8Ev3VLUv6W9wQ6NkOuxj/S5sSGl7y7Vl96wHnU2E7799CnsPjOCt1wm793Id1Who8mBvV+8Kef3qNZmRLXVkBxGMuAO4Rcvn8OfbG7J6o8rN2qmjf1sRJRqSVV2pq3Sd7QBgNmgQzAaRyIhJbOGDNqI5MXTd125FD/a1QUAuHV99vqhxWJFrRV2kw77eybw1q2Te9f6JkI5d7RdTFID1inLI5UySLcyjGTcH0GVxTDvwb5Jr51yemQwIgd1Zn3pgrYFKo+M512qncplMaDWZkAomkBHkyOrnG5JlRk2ow6PHR0AAKws0Oy3fbnchLmnezLbduiCG0LMrobapNdiRZ0NBy+4cXLAm3c0vl6rwaduWo0/fPxqNDpN+J9nzwAofIAWCmrlQSVy0PbdZ04jkZDwyRtXz/jfMV82t7mwotaKy1eUb1MsEc0/p1kPq0GbDNpG/eGSlpUsFPUNO5RSSqOWR1oZtNFF7iPXrYJZr0VCKs1qoXKl0QhsanVljf0v5Y62SqVm2mptximTOa5kpk3pfZ7nxdoqo27qPW0BZR9nKTNtC1QemZhycqRKzbbl2scihMDaRjs8oRj0WoHWqvwH/oZmBywGbVqJ5MHzE1hZZ4PdNLtG93VNDrxwegTRePYQkkwdTQ7c/9Ed+Itb1+BPL1sy490N7Q12nBr0oXciiPv2nMefbm1Fa7VlRo81ny5d4sJTn79uQX7IiKh8CSHQ7DJPlkd6F095JIC0EklfJAajTlOS3mOiSlZtNeCDO1dCpwQ1i9nmVhdODnjweMri696JYMl2tFUq9cJWRxErPFzKhHh17P+YP7wg55NFZdrmoDxyQd4xEpKE1UUGbWpfW2Y/m0otkWyrthQcbKHTanDZ0qrkvjZJkrKGkMxUR5MDcWUiUDFZO71Wg49dvwr/+acbZ5zSbW+wwRuO4W/vPwwA+PgNq2b0OERE5aKlyozeiSACkRiC0XjFj/sHJk9IUsf++0Ix2Dk5kggA8IkbVuHJz12bbF1ZrP50ayuW11rxgZ/txTt+8DJePjsKbyh20WfaNBqBZqcJly3NHryXyW7SQ4jJoG08EF2woG2qTFuyp63cyiOFELcJIU4KIU4LIb5QzOe0F7m3YDJoy92cqC7XXFHESM0rVtTg5KAX4/4IeieC8oLn1tmn49WrAw6TDq3V8/PDpwa9T58cxtu2t2YtAyciqjQtLjloG1UWay+OnjalPDLlqqwvHGNpJJFCoxFYWmCh8mLRWm3BI5/eiS/fvh4nBjy48/svAZibcf+V5uFP7cTHrp86+aDVCDhM+uRkxlHfwmTa5EEkhTNt6oW6UgZts37XEEJoAXwHwM0ALgB4RQjxB0mSjhX6vGLLI+/Y3AKpwMLFtY1yMFfMHoTkvrbusWRmbGMJ0vHrlGzfhpbcQ0jmgvr9M+g0RR3oRETlrqXKjIlAFD1jAQBYFOWR+TJtHEJCdPFRF1/fvqkF3336NB481F+Siq9K57QU36bksugxHoggGk/AE4qVfXmkpcymR24HcFqSpLMAIIS4D8DtAPIGbTqNSO5amIrLYsB7r16e9+PrmhzY2OrCtUXsErt0iRNGnQZ7usag0wgYtJpk0DcbdXYjNrQ4cMPa+lk/VrGqrQasa3Lgxo56NCzykgIiujioFQPqZN/FELQle9oyMm0M2oguXuri679+bcdCP5WK4zLrMRGIYlzZ1Vaug0jUPmZTmQVtLQDOp/z9AoDLC32CqYSpQrNBi99/bEdR9zXqtNjc5sKerjFYDFp0NDtK0gguhMCDn7hm1o8zXX/85NXz/jWJiOaKGrQdOu8GANQsovLItEEk4diMh1AREV3MnBYDJoJRjPkXMGhTetokScpbYRecg/LIUvS05Xq2UtadhPigEGKvEGKvPh4qwZedme3La3C0z41DF9zYtKSyx8sKIRbtIkoiuvi0VKVn2hbDlFk1aAtkBG02DiIhIpq2Kose7kBkQYM2k14Onwpl24LROPRaAX2BIYnTVYpHugCgNeXvSwD0Zd5JkqTvS5K0VZKkrUsaFm5H1xXLq5GQ5G9mKfrZiIioNOrtJug0Av3uEOwmXUmrMhaKRS8HZ6nlkX6WRxIRzYjLrF/4TJtOfm8KRwsHbaV+DytF0PYKgNVCiOVCCAOAtwH4Qwked05sbquCXitnpy5l8ycRUdnQagQanXLZ4GLoZwMAk0F+mw0qi1YBwBtipo2IaCacFgPcwShGvGEAC5tpC8XyDyMJRuIlLY0EShC0SZIUA/BxAI8COA7gfyVJOjrbx50rZoMWly5xwW7UYUXt4h8xS0RUSdS+tsUw7h8ALIb0TFsklkA4loDNwKCNiGi6XGY9JAk4p0wZrrKUb6atlJMjgdIMIoEkSQ8BeKgUjzUfPndzO/rdIWg07AcjIionLVVmoAuosS6OTFvmyH9/WM64MdNGRDR9LmU9wNlhPxwmXUl7xopVbKat1OWRF+W7xlWrahf6KRARUQ5L1EybfXFk2rQaAYNOk8y0+dSgjT1tRETTpgZtXSN+1CxQGb2pyEybucSZtvkPT4mIiPJoVoK2xZJpA+Rdber4ZzVoszPTRkQ0bU6zfEHvwngAVdNYyl1KxiIzbaUuj2TQRkREZUMd+79YetoAuUQyM2izMtNGRDRtaqYtIQHVC3RxTy17DEULBG3RMhxEQkREVCprGu2wGXVY1+xY6KdSMmaDFgGWRxIRzZrLPJldq7YuUKZNp+xpK1QeyZ42IiJazOrtJhz5x1sX+mmUVFqmLcTySCKimXKmBW0LnGkrVB45B9MjmWkjIiKaQ7l62lgeSUQ0fTqtJnnRq2YBdrQBRWbaWB5JRERUWUz6lPLIEMsjiYhmQ+1rq1qgoK2YTFsgEoeJmTYiIqLKYTFoEcrMtHG5NhHRjLiUCZILlWlTR/6H8mTa4gkJkVgCFn1pX+cZtBEREc0hs16LQFQO1nzhGGxGHTQascDPioioMi10pk0d+R/Ok2lTp0qaDaUNsxi0ERERzSGzQYdgRL4i6wvFYDWWtmSGiOhiog4jWeietnyZtoBSWcGeNiIiogoiDyJRMm2RGPvZiIhmQc20VS9Q0CaEgEGnmTLTVuqR/wzaiIiI5pBZr0UwGockSfCFYrCZFma3EBHRYrCm0YHWanPJR+pPh0mnyTs9MqgEbZYS9y7zch8REdEcMhu0SEhAOJaALxyDnZk2IqIZ+7MrluJdl7dBiIXrDTbptcmMWqZkeSR72oiIiCqH2tcQjMTZ00ZEVAILGbAB8jCScCxPpi3C8kgiIqKKo5bwBKNxZXokyyOJiCqZSZc/0xaao/JIBm1ERERzyKwEbYGIHLTZTSyPJCKqZIUybZweSUREVIHSyiPDLI8kIqp0hTJt6iASBm1EREQVRC2RGQ9EEE9ILI8kIqpwhQaRJIO2Ek+3ZNBGREQ0h9QJYsPeMADAxvJIIqKKZtQVGkQi7+Vk0EZERFRBzHo5SBtSgjaO/CciqmwFM20ROZhjeSQREVEFUa+2qpk2K4M2IqKKVjDTFo3DoNNAqyntWgIGbURERHNIHfk/7FPKIxm0ERFVNKNei1A0f3lkqbNsAIM2IiKiOaUuWB3yhACAI/+JiCqcUadBuMAgEgZtREREFYaZNiKixcWk1xYoj0wkX/dLiUEbERHRHNJrNdBrBXvaiIgWCZNeg0g8gXhCyvpYMBJLVliUEoM2IiKiOWbSa+ENyWOgWR5JRFTZjDo5KIvkyLYFo/GSj/sHGLQRERHNObVURqcRMOr41ktEVMlMevl1PNfY/2AkzvJIIiKiSqQ2pVuNOghR2jHQREQ0v9RMWyiWHbQFInGWRxIREVUis0EuieQQEiKiymc1ykGZPxzL+liI0yOJiIgqk1kppWE/GxFR5XNZDACAiUA062PBKMsjiYiIKpKFmTYiokXDZdYDANzB7KCN5ZFEREQVypTS00ZERJXNZZGDtlyZthCnRxIREVUmtVTGxvJIIqKK5zIr5ZEZmbZoPIFoXIKFmTYiIqLKowZtdmbaiIgqnt2kgxCAOxBJu11dAcBMGxERUQVSyyPZ00ZEVPk0GgGHSZ+VaQtG5KCNPW1EREQVSM20saeNiGhxcFn0WT1tQSXTxumRREREFUjd2cOR/0REi4PLnCPTppZHMtNGRERUedT+BpZHEhEtDk6LIWvkf0Atj2SmjYiIqPKYWR5JRLSouMz67EEkStDG6ZFEREQViCP/iYgWF5elQHkkM21ERESVx6yXgzWO/CciWhxcZj3cwSgSCSl5m1oeyZ42IiKiCrRtWRXevGUJ1jU7FvqpEBFRCTjMekgS4A3Fkrcx00ZERFTBamxGfPWtG2ExMNNGRLQYuCwGAMBEcLKvLcTpkUREREREROXBZdYDQNqutmR5ZLll2oQQ/yCE6BVCHFB+vbZUT4yIiIiIiKgcuSxy0JY69j+ojvzXlT5oK0WdxtclSfpKCR6HiIiIiIio7KlBW+oEyVA0DpNeA41GlPzrsTySiIiIiIhoGpxmuactdVdbIBKfk342oDRB28eFEIeEEHcLIary3UkI8UEhxF4hxN7h4eESfFkiIiIiIqL558zR0xaMxuds4NSUQZsQ4gkhxJEcv24H8N8AVgLYBKAfwFfzPY4kSd+XJGmrJElb6+rqSvX8iYiIiIiI5pVBp4HFoE0rjwwq5ZFzYcpQUJKkm4p5ICHEDwA8OOtnREREREREVOZcZn16pi0Sn5PJkcDsp0c2pfz1DgBHZvd0iIiIiIiIyp/TYoA7ZU9bMBKHRT835ZGzfdT/EEJsAiAB6Abwodk+ISIiIiIionLnMuvTR/5H43AovW6lNqugTZKkPyvVEyEiIiIiIqoULosep4d8yb8HI3E0OIxz8rU48p+IiIiIiGiaXBZ91iCSBZseSUREREREROkcZj3cgSgkSQKgTo8sw0EkREREREREFyOX2YBIPIFgNA5AmR7JoI2IiIiIiKg8uCyTC7YlSUIwGofZMDfhFYM2IiIiIiKiaXKZJ4O2aFxCPCGxp42IiIiIiKhcOJVMmzsYRTAil0iyp42IiIiIiKhMuMwGAIA7GEn2tbGnjYiIiIiIqEyk9rSpQZvFwKCNiIiIiIioLDjVnrZgFIFIDADLI4mIiIiIiMqGxaCFXiswEYgipJZHMtNGRERERERUHoQQcJoNck9bJAGA5ZFERERERERlxWXRYyIwWR7JQSRERERERERlxGXWyyP/oxz5T0REREREVHbUTFuI0yOJiIiIiIjKj0PJtAUi3NNGRERERERUdlxmAyYCKcu1mWkjIiIiIiIqHy6LHv5IHJ5gDEIARt3chFcM2oiIiIiIiGbAZZEXbA+4gzDrtRBCzMnXYdBGREREREQ0A06zHLT1u0Nz1s8GMGgjIiIiIiKaEZfFAAAY9ITmrJ8NYNBGREREREQ0Iy5m2oiIiIiIiMqXWh4ZjiWYaSMiIiIiIio36iASYO52tAEM2oiIiIiIiGbEbtJDHRjJTBsREREREVGZ0WoEHCY528ZMGxERERERURlSSySZaSMiIiIiIipD6gRJZtqIiIiIiIjKkFPZ1cagjYiIiIiIqAypY/8tLI8kIiIiIiIqP2p5pIlBGxERERERUflJDiJheSQREREREVH5YXkkERERERFRGXMpg0hMzLQRERERERGVH478JyIiIiIiKmNLaywQAmh2mefsa+jm7JGJiIiIiIgWudUNduz94k2osRnn7Gsw00ZERERERDQLcxmwAQzaiIiIiIiIyhqDNiIiIiIiojLGoI2IiIiIiKiMMWgjIiIiIiIqYwzaiIiIiIiIyhiDNiIiIiIiojLGoI2IiIiIiKiMMWgjIiIiIiIqYwzaiIiIiIiIyhiDNiIiIiIiojImJEma/y8qhBfAyXn/wouTE4B7oZ/EIlELYGShn8QiweOyNHhMlg6PydLhcVk6PC5Lh8dl6fC4LJ2pjsulkiTVFfNAutI8n2k7KUnS1gX62ouKEOL7kiR9cKGfx2IghNjL47I0eFyWBo/J0uExWTo8LkuHx2Xp8LgsHR6XpVPK45LlkZXvgYV+AkQ58LikcsNjksoRj0sqRzwuyxCDtgonSRJ/sKjs8LikcsNjksoRj0sqRzwuy9NCBW3fX6CvS1QIj0sqNzwmqRzxuKRyxOOSylHJjssFGURCRERERERExWF5JBERERERURkrSdAmhLhbCDEkhDiScttGIcSLQojDQogHhBAO5fZlQoigEOKA8ut/Uj7nTiHEISHEUSHEf5TiudHFazrHpfKxS5WPHVU+blJu53FJJTPN18t3prxWHhBCJIQQm5SP8bikkpnmcakXQvxUuf24EOKvUz6HxyWVxDSPSYMQ4sfK7QeFENelfA6PSSoZIUSrEOJp5bXvqBDiU8rt1UKIx4UQp5Tfq1I+56+FEKeFECeFELem3D69Y1OSpFn/ArATwBYAR1JuewXAtcqf3wvgn5Q/L0u9X8r9awD0AKhT/v5TADeW4vnx18X5a5rHpQ7AIQAblb/XANDyuOSvUv+aznGZ8XmXADir/JnHJX+V9Nc0Xy/fAeA+5c8WAN3KezuPS/4q2a9pHpMfA/Bj5c/1AF6FnJjgMclfJf0FoAnAFuXPdgCdANYB+A8AX1Bu/wKAf1f+vA7AQQBGAMsBnJnp+WVJMm2SJD0HYCzj5jUAnlP+/DiAN0/xMCsAdEqSNKz8/YkiPocor2kel7cAOCRJ0kHlc0clSYqDxyWV2CxeL98O4F7lzzwuqaSmeVxKAKxCCB0AM4AIAA94XFIJTfOYXAfgSeXzhgBMANgKHpNUYpIk9UuStE/5sxfAcQAtAG6HHHhB+f1Nyp9vh3yRKyxJUheA0wC2YwbH5lz2tB0B8Eblz38KoDXlY8uFEPuFEM8KIa5RbjsNYK1SPqmD/I9N/RyiUsh3XLYDkIQQjwoh9gkh/lK5ncclzYdCr5eqOzEZtPG4pPmQ77j8PwB+AP2QrxR/RZKkMfC4pLmX75g8COB2IYROCLEcwGXKx3hM0pwRQiwDsBnAywAaJEnqB+TADnLGF5ADuvMpn3ZBuW3ax+ZcBm3vBfAxIcSrkNOHEeX2fgBtkiRtBvBZAPcIIRySJI0D+AiAXwF4HnK5RWwOnx9dnPIdlzoAVwN4p/L7HUKIG3lc0jzJd1wCAIQQlwMISJJ0BAB4XNI8yXdcbgcQB9AMudznc0KIFTwuaR7kOybvhnwyvBfANwDsBhDjMUlzRQhhA/AbAJ+WJMlT6K45bpNmcmzqZvZUpyZJ0gnIJWcQQrQDeJ1yexhAWPnzq0KIM5CzHHsleZnfA8rnfBDymwJRyeQ7LiG/2D8rSdKI8rGHINfSP8njkuZageNS9TZMZtnUz+FxSXOqwHH5DgCPSJIUBTAkhHgBcinaWR6XNJcKnFvGAHxGvZ8QYjeAU8rHeExSSQkh9JADtl9KkvRb5eZBIUSTJEn9QogmAEPK7ReQnkFbAqAPmP6xOWeZNiFEvfK7BsDfAvgf5e91Qgit8ucVAFYDOJvxOVUAPgrgh3P1/OjilO+4BPAogEuFEBYlTX0tgGMZn8PjkuZEgeNSve1PAdyX53N4XNKcKHBc9gC4QcisAK4AcCLjc3hcUskVOLe0KMcihBA3Q86y8T2cSk4IIQD8CMBxSZK+lvKhPwC4S/nzXQB+n3L724QQRqV0dzWAPcpjTevYLEmmTQhxL4DrANQKIS4A+BIAmxDiY8pdfgvgx8qfdwL4shAiBjmi/LBSCw8A3xRCbFT+/GVJkjpL8fzo4jSd41KSpHEhxNcgT6aSADwkSdIflfvxuKSSmebrJSC/Zl6QJOlsxkPxuKSSmeZx+R3lz0cgl/78WJKkQ8rHeFxSSUzzmKz//+3coQ1CQRBF0Td1UQcShaQDKqCW3wFdkBBKoAvCIvYLLAGSEefIFavG3N1kkpyr6pnknmT3dpWZ5Jc2mfN1rarLenZMckqyVNU+82FrmyRjjFtVLZkfAY8kh3XRXfLhbNa6ZhIAAICG/rmIBAAAgC+JNgAAgMZEGwAAQGOiDQAAoDHRBgAA0JhoAwAAaEy0AQAANCbaAAAAGnsBD8quIkNv1l0AAAAASUVORK5CYII=\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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:23.925654Z", "iopub.status.busy": "2021-02-02T06:51:23.924811Z", "iopub.status.idle": "2021-02-02T06:51:23.988163Z", "shell.execute_reply": "2021-02-02T06:51:23.988577Z" } }, "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, 02 Feb 2021 AIC 951.427\n", "Time: 06:51:23 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 (L1) (Q): 8.43 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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:23.998223Z", "iopub.status.busy": "2021-02-02T06:51:23.994602Z", "iopub.status.idle": "2021-02-02T06:51:24.000262Z", "shell.execute_reply": "2021-02-02T06:51:24.000610Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.013417Z", "iopub.status.busy": "2021-02-02T06:51:24.010140Z", "iopub.status.idle": "2021-02-02T06:51:24.015640Z", "shell.execute_reply": "2021-02-02T06:51:24.016041Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.024942Z", "iopub.status.busy": "2021-02-02T06:51:24.022725Z", "iopub.status.idle": "2021-02-02T06:51:24.026925Z", "shell.execute_reply": "2021-02-02T06:51:24.027266Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2009Q4 3.689210\n", "2010Q1 3.772434\n", "Freq: Q-DEC, Name: predicted_mean, dtype: float64\n" ] } ], "source": [ "print(res.forecast(steps=2))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.048509Z", "iopub.status.busy": "2021-02-02T06:51:24.030346Z", "iopub.status.idle": "2021-02-02T06:51:24.050564Z", "shell.execute_reply": "2021-02-02T06:51:24.049127Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.057080Z", "iopub.status.busy": "2021-02-02T06:51:24.056399Z", "iopub.status.idle": "2021-02-02T06:51:24.065332Z", "shell.execute_reply": "2021-02-02T06:51:24.065992Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2009Q4 3.689210\n", "2010Q1 3.772434\n", "2010Q2 3.826039\n", "Freq: Q-DEC, Name: predicted_mean, dtype: float64\n" ] } ], "source": [ "print(res.forecast('2010Q2'))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.083157Z", "iopub.status.busy": "2021-02-02T06:51:24.073596Z", "iopub.status.idle": "2021-02-02T06:51:24.088356Z", "shell.execute_reply": "2021-02-02T06:51:24.089228Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.108553Z", "iopub.status.busy": "2021-02-02T06:51:24.106925Z", "iopub.status.idle": "2021-02-02T06:51:24.436093Z", "shell.execute_reply": "2021-02-02T06:51:24.437070Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3IAAAEvCAYAAAAAWPPhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABo/klEQVR4nO3dd1icVdoG8PvMDMNQh6FDqGmkF0ISk1iTWNeY6FriunbXtrq2b13bqquuuq7r2teNfe09iRp7N1GTQCAVSEgIEHofGGDa+f6YGSQJfd5pcP+uay6Gt56Bl2Ge95zzPEJKCSIiIiIiIgocKl83gIiIiIiIiIaGgRwREREREVGAYSBHREREREQUYBjIERERERERBRgGckRERERERAGGgRwREREREVGA0fi6Af2JjY2VGRkZvm4GERERERGRT+Tm5tZLKeMOXe7XgVxGRgY2b97s62YQERERERH5hBBif2/LObSSiIiIiIgowDCQIyIiIiIiCjAM5IiIiIiIiAIMAzkiIiIiIqIAw0COiIiIiIgowDCQIyIiIiIiCjAM5IiIiIiIiAIMAzkiIiIiIqIAw0COiIiIiIgowDCQIyIiIiIiCjCawW4ohHgBwKkAaqWU05zLogG8BSADQCmAs6WUTb3sexKAxwCoATwnpXzQ7ZYTEREREdGoJ6Xs/up69Py+t+0OXTaY9b19DwB2u31Qy3oeo7fj9NMO0dt2gw7kALwE4EkA/+ux7BYAX0kpHxRC3OL8/i89dxJCqAE8BeB4ABUANgkh1kopdw7h3EREREREFOCklLDb7d1fD31ut9thtVoP+grgoPWuAMf1fW/nEKLX2KfP9QPtM9B2g102mHU9dXV1AX2Mohx0ICel/F4IkXHI4uUAjnU+fxnAtzgkkAMwD8AeKeVeABBCvOncj4EcEREREdEIYLPZYLPZYLVaYbVa0dXV1f3cZrN1B2Yu/QVEKpUjblGpVN3bCCEghIBGozno+8EGRIGqv567ofTI9SZBSlnlPEmVECK+l23GACjv8X0FgPlunpeIiIiIiLxEStkdqNlsNpjNZpjNZnR2dsJiscBms0EI0R2gqdVqqFSq7odGo0FwcPCID7y8yd1AbjB6+231GVoKIS4HcDkApKWleapNRERERETUC4vFgq6uLnR1daGzsxNmsxkWiwXArz1prmBNo9EgJCSEAZoPuBvI1Qghkpy9cUkAanvZpgJAao/vUwBU9nVAKeUqAKsAICcnp+++RCIiIiIiUoTZbEZHRweamprQ1dV1UK+aqzeN/Iu7gdxaABcCeND5dU0v22wCMEEIkQngAICVAH7n5nmJiIiIiGiYpJTo6upCe3s7WlpaYLFYoFKpoNVqERER4evm0SAMpfzAG3AkNokVQlQAuAuOAO5tIcSlAMoAnOXcNhmOMgOnSCmtQohrAHwGR/mBF6SUO5R9GURERERE1B+73Y7Ozk60tbWhtbUVdrsdKpUKwcHB0Ol0vm4eDdFQslae28eqJb1sWwnglB7frwOwbsitIyIiIiKiYbPZbOjs7ERrayuMRiOklNBoNNDpdN3ZISkweSPZCREREREReYndbu/udWtvb+9O2x8WFsakJCMIAzkiIiIiohHCYrHgwIED6OrqQnBwMOe7jWAM5IiIiIiIRgCTyYQDBw5ArVYzgBsFGMgREREREQUwKSWam5tRU1OD0NBQaDT8iD8a8LdMRERERBSgbDYbamtr0dLSgvDwcCYwGUX4myYiIhqGbRUtuHvtDkgpfd0UIhqlzGYzysvL0dbWhsjISAZxowx/20RERMPwXl4FXtpQiurWTl83hYhGofb2dpSWlsJutyMsLMzXzSEfYCBHREQ0DEXVRgDAvvp2H7eEiEYTKSUaGxtRVlYGnU7HQt6jGAM5IiKiYSiucQRypfUmH7eEiEYLm82GyspK1NXVITIykklNRjn+9omIiIaoztiFhnYzAGBffZuPW0NEo0FXVxcOHDgAu93O0gIEgIEcERHRkLl644QA9rFHjog8zGg0orKyElqtFqGhob5uDvkJBnJERERDVOicH5edZkBpA+fIEZFnSCnR0NCA+vp6hIWFQa1W+7pJ5Ec4R46IiGiIiquNiAnTIifdgLIGE2x2liAgImVZrVZUVlaioaEBERERDOLoMAzkiIiIhqiwxoiJCRHIiA2D2WZHZXOHr5tE5LZb39+Kz3dU+7oZBMd8uLKyMnR2diIiIgJCCF83ifwQAzkiIqIhsNsldtcYkZUYgYwYR+0mDq+kQFdr7MQbG8vx2Y4aXzeFADQ1NUFKiZCQEF83hfwYAzkiIqIhqGjqgMlsQ1ZiBDJjnYEca8lRgNtS1gwAqG5l77I/kFJCpeLHdOofrxAiIqIhKHJmrMxKjEBCZDBCgtTMXEkBL6+sCQBQ1dLp45YQ0WAxkCMiIhqCoupWAMDEBMe8lfSYUA6tpIC3ZX8zAKC6pRNSMnkPUSBgIEdERDQERTVtSDGEIDzYUcFnbFwYh1ZSQLPY7Nh6oBnBGhVMZhuMXVZfN4mIBoGBHBER0RAUVxuRlRDR/X1GTBjKGk2w2uw+bBXR8BVWGdFpsePYrDgAjl45IvJ/DOSIiIgGyWy1o6SuDVmJPQK52DBY7RIHWIKAApRrftwp05MAcJ4cUaBgIEdERDRI++rbYbXLgwI5V+bKvRxeSQEqr6wJCZHByE4zAABqGMgRBQQGckRERINU6Ex0clCPXAxLEFBgyytrQnaaAQmROgDskSMKFG4HckKILCFEfo9HqxDi+kO2OVYI0dJjmzvdPS8REZG3FdcYoVEJjI0N714WG65FeLCGgRwFpDpjF8obOzA7LQpajQqx4VrWkiMKEBp3DyClLAIwCwCEEGoABwB80MumP0gpT3X3fERERL5SVG1EZmwYtJpf74MKIZARG4p9DawlR4Fni3N+nGtYZaJex2QnRAFC6aGVSwCUSCn3K3xcIiIinyuqMR40rNIlI4YlCCgw5ZU1I0gtMG2MHgCQGBnCoZVEAULpQG4lgDf6WLdACFEghPhECDFV4fMSERF5VFuXFeWNHQeVHnDJjA1DRZMJZitLEFBgyStrwpRkPXRBagBAoj4Y1a0M5IgCgWKBnBBCC+A0AO/0sjoPQLqUciaAJwCs7uc4lwshNgshNtfV1SnVPCIiIrfsrjECQJ89cnYJlDdxeCUFDqvNjq0VzZidGtW9LEkfgmaTBZ0Wm+8aRkSDomSP3MkA8qSUNYeukFK2SinbnM/XAQgSQsT2dhAp5SopZY6UMicuLk7B5hEREQ1fUXU/gVwsM1dS4CmsdhQCz043dC9LdGau5Dw5Iv+nZCB3LvoYVimESBRCCOfzec7zNih4biIiIo8qqjEiJEiNVEPoYetcteT2MZCjAJLXnegkqntZkp4lCIgChdtZKwFACBEK4HgAV/RYdiUASCmfAXAmgKuEEFYAHQBWSimlEucmIiLyhqJqIyYmhEOlEoetM4QGIVKnQWkDAzkKHHn7mxAfEYwxUSHdyxKcgRxLEBD5P0UCOSmlCUDMIcue6fH8SQBPKnEuIiIiXyiuMWLxpPhe1wkhkBkXjtJ6zpGjwLGlvBmz06LgHDQF4NehleyRI/J/SmetJCIiGnHq27pQ32bGxF4yVrpkxoRyaCUFjPq2LuxvMHXXj3MJC9YgUqdBDQM5Ir/HQI6IiGgAxc5EJ5MSI/vcJiM2DJUtHcz2RwFhS1kzAByU6MQlUa9jjxxRAGAgR0RENIAiZ+mBiYnhfW6TGRsGKYGyRg6vJP+XV9YEjUpgurMQeE+J+hDWkiMKAAzkiIiIBlBUbYQhNAhx4cF9bpMRw8yVFDi2lDVhanJkdyHwnpIidSw/QBQAGMgRERENoKjGiKzEiIOSQhyKteQoUFhtdhSUt2B22uHDKgHH0Mq6ti5YbHYvt4yIhoKBHBERUT/sdoniaiOy+kl0AgD6kCBEh2lZgoD8XmG1ER0WG2b3qB/XU6JeBymBWmOXdxtGREPCQI6IiKgfB5o70G62IaufRCcuGcxcSQFgS3ch8L575ACguoW15Ij8GQM5IiKifhQ5M1Zm9ZPoxCUjNoy15MjvbSlrRlxEMFIMIb2uT+oO5NgjR+TPGMgRERH1oztj5QBDKwEgMyYM1a2d6DCzBAH5r7yyJsxOjepzzmdSpCPAq2KPHJFfYyBHRETUj6JqI8ZEhSBCFzTgtt0JTzhPjoboz+8U4M412z1+noa2LpQ2mHqtH+cSGaKBLkjFzJVEfo6BHBERUT+KnRkrByOTmStpGErq2vBObgXe3FiOFpPFo+fqLgTex/w4ABBCIEkfgirWkiPyawzkiIiI+mCx2VFS1zaoYZXArz1y+9gjR0Pwvw2lEAIw2+xYt73Ko+faUt53IfCeEiN1qGGPHJFfYyBHRETUh3317bDYJCYNskcuPFiD2PBg9sjRoLV2WvBubgVOnz0GY2PDsHrLAY+eL29/MyYnRSJEe3gh8J4S9TpUMZAj8msM5IiIiPrgylg52B45ABgbG8YSBDRo72yuQLvZhosXZmL5rDH4ZV8jKps9k2TEarOjoKIZ2X3Uj+spUa9DTWsn7HbpkbYQkfsYyBEREfWhqNoItUpgXHzYoPfJiA3FPpYgoEGw2SVe3lCKnHQDpqfosXxWMgBgbUGlR85XVGOEyWzrN9GJS5JeB6tdoqHd7JG2EJH7GMgRERH1oajGiMzYMARr+h+G1lNGbBjq27pg7PRs0goKfF8X1qKs0YSLF2UCcFw7s1KjPDa8Mm8QiU5cEiNdteQ4vJLIXzGQIyIi6kNRtRFZQxhWCThqyQHA/gb2ylH/XtqwD0l6HU6cmtC9bMWsZBRWG7uH9SppS1kTYsO1fRYC7ynRWRScteSI/BcDOSIiol6YzFaUNZoGXXrApTtzJefJUT+Kqo1Yv6cB5y9Ih0b968exU2cmQ60SWJ2vfK/clrJmzE4z9FkIvCdXIFfNEgREfouBHBERUS+Ka9oADC3RCQBkxLCWHA3spQ2lCNaocO7ctIOWx4YH48jxsVibX6loopHGdjP21bcPalglAMSGBUOjEhxaSeTHGMgRERH1otg5tG2wpQdcQrRqJEbqWEuO+tRsMuODLY6SA4Yw7WHrV8xOxoHmDmze36TYObeUOY41mIyVAKBSCSRE6hjIEfkxBnJERES9KKw2QhekQmp06JD3zYgNZY8c9emNjeXotNhx0aKMXtefMCURIUFqRYdXbilrhlolMD2l/0LgPbGWHJF/YyBHRETUi+IaIyYmRECtGng+0aEyY8NQymQn1AurzY5XfirFwnExmJQY2es2YcEaHD8lAeu2VcFstSty3ryyJkxOikCoVjPofRL1Os6RI/JjDOSIiIh6UeQM5IYjIyYMje1mtHSwBAEd7POdNahs6cRFCzP63W7F7GQ0myz4rrjO7XPa7BIF5c2Dnh/nkuQcWikli4IT+SMGckRERIdobDejztg15PlxLq7MlRxeSYd6aX0pUqNDsGRyQr/bHTUhDtFhWkWGVxZVG9Futg05kEvU69BhsaG1w+p2G4hIeYoEckKIUiHENiFEvhBicy/rhRDicSHEHiHEViFEthLnJSIi8gRXDa/h9shlugI5JjyhHrYfaMHG0kZcuCBjwCG7QWoVfjM9CV/urHG7uPyWcleik6EHcgBQ1cpackT+SMkeueOklLOklDm9rDsZwATn43IA/1HwvERERIoqqm4FgCHXkHNJiw6FEKwlRwd7cX0pQrVqnJWTOqjtV8xORpfVjs921Lh13rz9zYgJ0yI1euBC4D0luWrJMeEJkV/y1tDK5QD+Jx1+BhAlhEjy0rmJiIiGpKimDVGhQYiPCB7W/rogNZL1IQzkqFt9Wxc+LKjEmXNSoA8JGtQ+2WkGpEaHYI2bwyu3lDUNuhB4T4l6R+DHQI7IPykVyEkAnwshcoUQl/eyfgyA8h7fVziXERER+Z2i6lZMTIgY8gffnjJjwzhHjrq9/ksZzDY7LhwgyUlPQggsnzkG6/fUo9Y4vGCqqd2MvfXtyE6PGvK+8RHBEAIsQUDkp5QK5BZJKbPhGEL5RyHE0Yes7+0/Ya8pkIQQlwshNgshNtfVuZ+piYiIaCiklCiuaRt2ohOXjNhQ7KtvZ8Y/gtlqx6s/78cxE+MwLi58SPuumJ0MuwQ+LKga1rnzy5sBDH1+HOCYpxcbHsweOSI/pUggJ6WsdH6tBfABgHmHbFIBoOeA8BQAlX0ca5WUMkdKmRMXF6dE84iIiAbtQHMH2rqsw0504pIRE4bWTiuaTCxBMNp9sr0KtcYuXNxHAfD+jI+PwNTkyGEPr8wra4JaJTBjCIXAe0piLTkiv+V2ICeECBNCRLieAzgBwPZDNlsL4AJn9sojALRIKYd3a4mIiMiDimscGSvd7ZFzZa7kPDl6YX0pxsaG4egJw7tBvWLWGGytaMHeurYh75tX1oRJiUMrBN5TorOWHBH5HyV65BIA/CiEKACwEcDHUspPhRBXCiGudG6zDsBeAHsAPAvgagXOS0REpLhCZ+mBCe72yLGWHMGRaKSgvBkXLcqAaoCSA305bVYyhABW5/c6mKlPNrtEftnQC4H3lKjXoaqF5QeI/NHwbs/0IKXcC2BmL8uf6fFcAviju+ciIiLytOJqI5L1ukFnFuxLqiEUKsFacp5msdnx3+9KcNykeExNHt7wQU96cX0pIoI1OCM7ZdjHSIjUYeG4GKzJP4Ablk4YdBKe3bXOQuDDSHTikqjXobXTCpPZOuxePSLyDG+VHyAiIgoIhdVGTHRzWCUAaDUqpBhCObTSgyw2O657cwse/rwYl/8v1+3C2Uqrae3Eum1VOHtuKsKD3QuCls8ag/0Npu7kJYORt9+x7ezU4ffIsZYckf9iIEdERORksdmxt6592IXAD5URG8YeOQ+x2uy4/q18rNtWjd/NT0NVSwfu+XCnr5t1kFd/3g+blLhwQYbbxzppWiK0GhXWDGF4ZV5ZE6LDtEiPCR32eRMiGcgR+SsGckRERE77G9phttmR5eb8OJfMmFCU1ptYgkBhVpsdN75dgI+3VuH2Uybj/tOn4+pjx+Od3Ap8vqPa180DAHRabHj9lzIsmZSANDcCKZdIXRCWTo7HR1srYbXZB7VPXlkTstOi3KqHmOQsCs5ackT+h4EcERGRkyvRibulB1wyYsPQ1mVFfZtZkeORI4HH/71TgLUFlbjl5En4w9FjAQB/WjIBU5Mjcev721Bn7PJxK4EPCyrR0G7GJcMoOdCX5bPGoL7NjB/31A+4bbPJjL117ZjtRqITwJG1EgBLEBD5IQZyI5TdLrGptBGPfF6E8kaTr5tDRKQ4u135Xq7iaiNUAhgfP7SizX3pzlzJ4ZWKsNkl/vxOAVbnV+LPJ2bhymPGda/TalR49JxZMHZZcev7W33aCyqlxIvrS5GVEIEF42IUO+6xWXGI1GkGNbxyi3Mu3ey0KLfOGaJVIyo0iEMrifwQA7kRREqJvLIm3PPhTix88Guc9cxPePzrPbjghY1oaPP93UkiIqV8U1iLefd/ha0VzYoet7DaiIzYMOiC1IocLzPGWUuujoGcu+x2ib+8txXvbzmAm46fiD8eN/6wbSYkRODmE7Pw5a5avL253AetdNhU2oSdVa24aFGGW8MaDxWsUeM3M5Lw2Y5qmMzWfrfdsr8JKgHMTIly+7yJkToOrSTyQwzkApyUEtsqWvDAul048h/f4IynN+DVn/dj2phIPHrOLPzvknmobO7AJS9vHvBNn4goEJjMVtyxejvq27pw30e7FO15Ka4xul0IvKcUQwg0KoF97JFzi90uccv7W/FubgWuXzoB1y6Z0Oe2lyzKxIKxMbjnw50oa/DNiJQX1+9DVGgQVswao/ixl88aA5PZhi921vS7XV5ZMyYlRiLMzWyZgKMEQXUra8kR+RsGcgFISoldVa3452eFOO7hb7HsyR/x/I/7MCEhHA+fNROb7liK5y6cixWzx+DoiXF4/NzZ2FbRjGtf3zLoCdJERP7qsa9240BzB86ak4KNpY34cletIsc1ma3Y32hSbH4cAGjUKqRFh/ptUfDyRhNufDsf9X48asNul7h99Ta8vbkCf1o8Htcvndjv9iqVwMNnz4RKCNz0Tj5sHhiC25+KJhM+21GNlXPTEKJVpme3p3kZ0UjW67B6y4E+t7HZJfLLm92qH9dTkl6H6hb/vUaIApmUsvsBADabDRaLBRaLBWazGVZr3x0xrOwYQPbUGvFhQRU+2lqJkrp2qASwcFwsrjxmHE6cmghDmLbX/U6cmoi/LZ+Gv67ejr+u2Y77T5+u6FCPgdr86s9luH7pBESF9t4+IqLBKqo24vkf9uHsnBTcf/p05JY14cFPduG4rDho1O7dm9xT2wYpoWiPHOCYJ+ePteTqjF04//lfUNpgwryMaKycl+brJh1GSom/rtmONzaW45rjxuOG4/sP4lzGRIXg7tOm4qZ3CvDsD3sPmkvnaa/8vB9CCJy/IN0jx1epBJbNSsZzP+xDQ1sXYsKDD9tmT20b2rqsyHYz0YlLYmQI6tu6YLbaodWwD2C0k1KitbUVZrMZZrMZXV1dMJvNMBgMSEhIQFdXFzZs2NAdjFitVlitVkydOhWTJk1Ca2sr3nnnne7lrsdxxx2H7OxsVFVV4emnn4bFYoHdbofNZoPNZsN5552HuXPnYs+ePXjooYcOWme323Hddddh7ty5yMvLw3333de9TkoJu92Ov//978jOzsZ3332He+65B3a7vfshpcQzzzyDadOmYe3atbj33nu7l7u+vvfeexg/fjxeffVVPPDAAwDQvU5Kia+//hpjxozBU089hUceeeSgAE1Kifz8fBgMBjzwwAN48sknD/u5lpaWIigoCH/961/x8ssvdy8///zz+/xdMJDzc6X17fhoayU+2lqFwmojhHDcjbtoUSZOnpaI2F7ewHtz/hHpqG7pwFPflCAhUjfgHU0lbNzXiMte3oTWTiu6rDY8cMYMj5+TiEYuu13ijtXbEKHT4JaTJ0OjVuGWkybh8ldy8eamcvz+CPc+OCudsdIlIyYMP5U0QErptZtoA2npsOCCFzaiprUL4cEabN7f5HeBnJQSd67Zgdd+KcNVx47DTSdMHNLP74zsMfhiZw0e+bwYx0yMw+SkSA+21sFktuLNjeU4cWoCxkSFeOw8K2aNwX+/24uPt1Xhgl5q1OWVNQGA2xkrXRL1js8aNa2dSI12v5QCed6+fftgNBphMplgMpnQ0dGB2NhYzJ8/HwDw9NNPo7GxER0dHd3bzJs3D5deeikAYOnSpWhrazsoUDvvvPNw9913o6urC1OmTDnsnNdccw1uvfVWtLe34/e///1h6//yl79g0qRJaG5uxp133nnYeoPBgOzsbBiNRrz33nvQaDTQaDRQqVRQqVRYunQpAKCzsxPFxcVQq9VQqVRQq9VQq9Uwmx3ZgYUQUKlUCAoK6l4nhIBW6+hQCA0NxdixYyGE6N5WpVIhNNRxbSckJOCoo46CSqXq3kYIgbAwx5znsWPHYsWKFd3rXed07T9jxgxcfPHFB+0LADqdIwPswoULodFoDtrX1Q7Xzz4+Pr57/aRJk/DKK6/0+ntmIOenWjosuGP1dnxY4MhMNSfdgLuWTcEp05O6i3MO1f+dkIXqli48+uVuJEbqPPpP++OtVbjh7XykGEKwdEoC3thYjjPnpGJOujL/VIho9Hk3rwKbSpvw0G9nINo5AuH4KQmYlxGNR78sxorZYxDuxnyg4mojgjUqpDsTlCglMzYUHRYbalq7kKgf3vu3kjrMNvzh5c3YU2vEcxfOxSs/7Ufe/iZfN+sgUkrcvXYHXvl5P644eixuPjFryEGwEAL3nzEdJ/z7e9zwVj7WXLMIwRrlhzr2tHpLJVo6LLh4UaZHzzM5KRJZCRFYveVA74Hcfkch8AwF6tcBQKKzllw1Azmvqa+vR2VlJTo6OtDc3IyWlhbodDosW7YMAPCPf/wDhYWFaGlpQUtLC5qbmzF9+nS89NJLAIDzzjsP+/fvP+iYxx9/fHcg9+yzz6K1tRUhISEIDQ1FSEgIMjN/vW4nT54MwBF8BAUFQavVdu8bHByMu+66C1qtFsHBwdBqtdBqtZgwwTF3NTIyEmvWrIFWq0VQUFB3QGYwOD4DjhkzBtu2bTtonSsoA4CJEydi586dff5spk2bhm+//bbP9bNnz8a7777b5/q5c+di7ty5fa6fP39+92vtzcKFC7Fw4cI+1x911FE46qij+lx/zDHH4Jhjjulz/eLFi7F48eLu702mvuf6MpDzQ5tKG3H9m/mobu3EtYvHY+W8NEXu7Akh8OBvp6OurQu3r96OuIhgLJmcoECLD/b8j/tw38c7MSfNgOcuzEGQWoUNexocgek1i9we/kRE7jOZrVCrBIJUKqhU/tFL1J+mdjMeWLcLOekGnDknpXu5EAK3/WYyVjy1Hqu+K8GNJ2QN+xxFNUZMSAiHWuGfh6sEwb76dp8HchabHde8nodN+xvx+MrZOGZiHHZWtuLLXTVobDd3B8i+JKXEPR/txMs/7cdlR2bilpMnDbsnMzpMi4fOnI5LXtqMR74oxq0nT1a4tb+SUuKlDfswbUwkcrxw03L57GQ89GkRyhpMhxUczytrwuxU9wqB95TkvG5ZgsB7rr76avz4448HLZs4cWJ3IFdSUoKKigro9XpkZGRAr9d3B18AcO+998JqtSIkJARhYWEIDQ3tDqQAIDc3tztw6s0TTzzR5zohBC6//PI+12s0GuTk5PS5Xq1WIzo6us/1NHgM5PyI1WbH41/vwZNf70aKIRTvXrlAsWERLkFqFf5zXjZWrvoZf3w9D2/84QjFzmG3S/x93S48/+M+nDQ1EY+unNWdwvuuZVNw1Wt5ePmn/bj0SM/eqRyJals7cev723DnsimK9xbQ6GJzpnB/N7eie5lGJRCkViFILaDVqJzPHd8HqVU9lgloNWpo1QKhWg2uXTweExQehtiXBz8phLHTivtOn3ZY4DkrNQqnzkjCsz/sw3lHpA971EJRtRFHTYhTorkHyYj5tZackjXFhspul7j53a34qrAW962YhmUzkwEAORmO/wG5+5tw/BTlb+4NhZQSf/94F15cX4qLF2Xg9t9MdjsYWTwpAefOS8Oq7/diyaQEzMv0zAfI9XsaUFzThofPmumVIbSnzXQEcmvyDxyUxbPFZEFJXTvOyE7pZ++hcf1NMZDrX0uHBXvr2hw3bSJ1WDg+dtjHuvrqq3H22WcjNjYWer0eer0eUVFR3etXrVrV7/5Llizpd31/QRwFDgZyfqK80YTr38pH7v4mnJE9Bn87bSoidEEeOVdYsAYvXDQXv/3PBlz68ma8d9VCZMa6Fxx0Wmy46e0CfLytChctzMBfT51y0F3tk6Yl4tisODzyeRF+Mz3J53elA83rG8vwVWEtOiw2vHbZfL+ZZ0OBxWqz46Z3CrAmvxK/m+/o6bfY7M6HhNlqh9lmh8XaY1n3ejssVonWDgssNjvKGk3YXNqI1dcsQnyEZ/+eN5c24q3N5bji6LGYlNj7PKebT5yEz3ZU45HPi/GPM4c+H7ep3YxaYxeyEpUpBN5TclQItGqVTzNXunq5PthyAP93wsSD5hNOH6NHkFr4PJCTUuKBTwrx3I/7cNHCDNx56hTF3uvu+M1kbCipx41v5+OT647yyP/XlzbsQ2y4FstmJil+7N6kGEIxLyMaq/MP4JrF47t/VlvKnfPjUqMUO1ekToNQrZq15AB0WW0oazBhb3079ta1Y199m/NrOxrazd3bnTYz2a1A7thjj0VnZyeCgweXC4FGJwZyfmBtQSVuf38bAOCxlbOw3AN1Zw4VFxGMly+Zh9/+ZwMueOEXvH/VIsRFDO/NosVkwR9e2YyN+xpx+ymTcdlRmYf98xVC4G+nTcUJ//4e9368E0/9LluJlzEq2O0S7+VVIEKnwYaSBnyw5YCid1ppdLDa7Ljx7QKsLajEn0/M6rWY8lBsP9CCs575CX94eTPevHyBR9KsA46hgHes3o5kvQ5/6qd2WFpMKC5YkIEX1+/DJUdmImuImSeLajyT6AQA1CqBtJhQn2aufOLrPXhpQykuPTLzsN+9LkiNqcl65O5v9FHrHEHcPz4twqrv9+L8I9Jx1zLlgjjAcQPzkbNn4qxnfsK9H+3EQ2fOVOzYJrMVq77fi68Ka3Ht4gken4fX0/LZybj9g+3YUdmKaWP0ABz141QCmKlgICeEGFW15Ox2ierWTuyrb8feurYeQVs7KppM6FnRIjY8GGPjwnD8lARkxoZhbFw4MmPDkMa5hOQFDOT6sLvGiGd/2IslkxOweFI8gjwwr6uty4o712zH+3kHkJ0WhcdWzvbqJOLM2DC8cNFcnLvqZ1z80ka8efmCIScKqGgy4aIXN6GswYTHz52N05xDdXqTHhOGPx43Ho98UYxzcupw9ETlhzCNRL/sa0R5YwceOXsmXvl5P+77eBeOy4rvs9wE0aGsNjtueLsAHxZU4uaTsnD1se4FcQAwbYwej62chStezcVN7+TjyXOzPTLX7sX1+1BYbcSq8+cMWNj42sXj8c7mcjzwyS68dPG8IZ2n2BnI9dXj566MmDCU+qgo+Cs/leKRL4rx2+wU3H5K70MVc9IN+N/P+32SXl5KiYc/L8Iz35XgvPlp+NtpUz0y6mBOejSuPGYcnv62BEsnJ+CEqYluHc9ml3hnczke+aIYtcYunDI9EZcd5d2pA7+ZnoS71+7A6i0HugO5LWVNyFKoEHhPjlpyI79HrqC8GRe+uBHNJkv3spAgNTJjwzAjRY8Vs5K7g7XMuDBEemj0FNFgMJDrRVVLBy54YSOqWjrx9uYKxEcE46ycFJyTk3bYhOLhyi9vxnVvbkF5owl/WjIBf1o83idJQGalRuGp82bjD//LxdWv5eF5Z3KSwdhR2YKLX9yEDosNL18yb1BzP644Ziw+2HIAd67Zjk+vP7p7Dh317d3cCoQHa3DytCRMSY7EqY//iAc/KRzW8DEafXoGcbecPEnRelonTE3EbSdPxt/X7cLDMUW4+aRJih0bACqbO/Dol7uxdHL8oD50R4Vqcc3i8bh/XSHW76nHoiEMayqsNiJSp0FCpGeGMWXGhuL73XWw26VXk8usyT+AO9fuwNLJCfjHb6f3ee456QY89+M+7KhsUXxu9kDWbavGU9+U4Nx5qbh3+eFzIJV0/dKJ+LaoDre+vw3Z6YZBl/DpSUqJb4vr8OC6QhTVGJGdFoX//D4bc9K9n7whKlSLYybGY21BJW49ZTIEgPyyZpw2q++bqsOVEKnDzyUNih/Xn0jpmOuvUalw34ppGOsM1hIjdZzSQH6JMx0P0dppwcUvboKx04q11yzCqvPnYNoYPf7zbQmO/uc3+P1zv+CjrZUwW+3DOr7NLvHUN3tw5n82wGK1483LF+DG4yf6NJPj4kkJ+PuKafi+uA63vLetu7J8f37YXYdz/vsz1CqB965aOOgJ/MEaNe5dPg2lDSY8812Ju00f8dq6rFi3rQqnzkhCiFaNSYmRuPSoTLy1uRwb9/luGBQFBqvNjuvfyseHBZW4VeEgzuWyozJx7rxUPP1tCd7ZXK7osf/24Q7YpcRdy6YOep8LFmRgTFQI7l+3C3b7wO9lLsXVRkxKjPTYh7WM2DCYrXZUtnhvaNo3RbW46e0CzMuIxpO/m93v/xlXaZhcH5Qh+HxnNWLDtfj7ir4DTaVoNSr8+5xZMHZacev7g/t/19OOyhac//xGXPziJnRabfjPedl476qFPgniXFbMTkatsQs/723Anro2GLusHgnGk/Q61Bi7YBvC31Wg+WF3PTbua8SflozH749Ix8LxsUjShzCII7/FQK4Hs9WOK/6Xiz21bfjP77MxIyUKJ0xNxAsXzcWPf1mMG5ZOxL76dlzz+hYseOAr3L9uF0rq2gZ9/KqWDvz+uV/wz8+KcOLURHxy3dEey541VCvnpeH6pRPwXl4FHv68qN9t38+rwMUvbkKKIQQfXL1oyHNKjpwQi2Uzk/H0tyU+nfwfCNZtq0KHxYazcn6dE3fdkglIMYTgtg+2DfuGgq9Ut3TCagusNgcqi82O697Mx0dbq3DbKZNwhQeCOMAxd+ae5dOwaHwMbvtgG37eq8wd+6921eCzHTW4bsnEIQ051wWpcfNJWdhR2YrV+QcGtY+UEkU1Rkz0QKITl0xX5sr6vusBKSl3fyOuejUXWYkRePbCnAFHP8RH6pAaHeL1QE5KiQ0lDVgwLtZrPZVZiRH484lZ+GJnDd7pkb21P1UtHbjp7QKc+sSP2F7ZgruWTcEXNxyDk6cn+fxD/tLJCQgP1mD1lgPd9QCz06IUP0+iPgQ2u0RDW5fix/YHriG+Y6JCsHKu5+rsEimJgZyT3S7x53cL8NPeBjx05ozDUlAnR4XguqUT8P3Nx+Gli+dibkY0XvhxH5b86zuc/d+f8MGWCnRabH0e/9Pt1Tj5sR9QUNGMh86cgSd/Nxv6UP8aV33dkgk4d14qnvqmBK/8VHrYeikdvYk3vl2A+WOj8faVC4adffKO30yGVq3CnWt3DPmO6Gjybm4FxsaGIbvH3dVQrQb3rpiGPbVtWPW9//dq7q1rw5Nf78bJj/2AIx74Cn95b5uvm+RTJXVt+GRbVb/vF+5yBHFb8PG2Ktx+ymRcfrRngjiXILUKT583B2nRobjilVzsHcINrt50mG24a+0OTIgPH1a5kmUzkjFtTCQe/qxoUD/nqpZOGDutyPLQ/DgAyIxz1pLzwjy5XVWtuPjFTUjSh+DlS+YNeg7PnDQDNu9v8up7ckldG+qMXVjk5bIMlx6ZifmZ0bjnw50ob+w7uDZ2WvDPzwpx7D+/xYdbK3H50WPx3Z+Pw8WLMr0+l7AvuiA1TpyaiE+3V+OnvQ0whAa5nYm6N4nOEgQjNXPlZztqsLWiBdctneA3v1uigfBKdfrHZ4VYk+/I5tZfRkC1SuDYrHg8c/4cbLh1MW4+KQs1rZ244a0CzL//K9y9dgcKq1u7t+8w23DbB9tw5au5SDWE4qNrj8TZOak+v4PXGyEE7l0+DUsmxePOtTvw6fbq7nVWZ+a4f35WhBWzkvHiRYP/cNCbhEgdbjphIr4vrsO6bdUD7zAK7W9ox8Z9jfjtnJTDrpfjsuLxm+lJeOLrPX7Zq1lS14YnvtqNkx79Hov/9R0e/rwYIUEqnDQ1Ee/lVWBtQaWvm+hVe2rb8Ljz57HkX9/hqtfycOw/v8UbG8sU76G02Oz40xtbsG5bNe74zWT84eixih6/L/qQILxw0VyoVQKXvrwZzSbzwDv14Ymvd6OiqQP3rZg2rA9UKpXAbadMRmVLJ17aUDrg9kXVjkQnWR6siZcQoYMuyPMlCPY3tOOCFzYiVKvBK5fOG9IcsDkZ0agzdqGiyXvDP9fvcfTgLhw3/DTtw6FSCfzrbEfmypveKThsuKDFZscrP5Xi2H9+i6e+KcHJ0xLx9U3H4NaTJ0Mf4l83YQHH8EpjlxUfFlRidprBI58xXEXBR2IgZ7NLPPJFEcbGheGM2Z7PHE6kFL9OdtLSYRl4IwW8vKEU//1uL86bn4arjx38nev4CB2uPnY8rjx6HH7e24A3NpXj9V/K8NKGUsxOi8KyGcl47Zf9KKlrxxXHjMVNx2f5/V0ejVqFJ343G7979hdc9+YWvHbZfExJjsSf3tiCL3fV4qpjx+HPJ2QpMgTm/CPS8c7mCtzz0Q4ckxU35IyZI917uRVQCeCM7N7/qdy5bAq+L67DX9dsx/8umefzmwN7ao34eGs11m2r6k7lPifdgL+eOgUnT0tEclQIrDY7zv7vT7j9g23ITotCimHkpmc+9OchhCMz4F3LpiA9JhRPfr0Ht76/Dau+34ubTpiIU6Yluf135QriPtnuCOIuO8o7QZxLekwYVp0/B7979hdc8UouXrl0/pDf83bXGLHq+734bXYK5o8dfi/NwnGxWDwpHk99swfn5KT2m+XVdb16MpBTqYQjc6UHA7na1k6c//xGWGx2vH7FgiH/fc1J+3WenLcyKG8oqUeKIUSxRGJDkWIIxd2nTcX/vVOA53/ci8uPHgcpJb7YWYMHPynE3vp2zM+Mxou/mYwZKVFeb99QLBwXi7iIYNQZuxStH9eTawROtRfneXrLhwWVKK5pG3AuKZG/8etPzmWNJtzz4U7cesokj6T/BxxDHu/+0JHR657l04b1YVilElg4PhYLx8eisd2M9/Mq8Oamctzz0U7ERwTj1Uvn48gJ3r3b6I5QrQbPX5iDM5/5CZe+vBnpMaHYfqAF9y6fivMXZCh2Ho1ahb+fPg1n/GcD/v1FMf566hTFjh3oHLXjDuDICXFI0of0uk1CpA43n5SFv67ZgbUFlV6pP3io3TVGfLytCuu2VaG4pu2gYOXkaYcXfteoVXhs5Wyc/NgPuOGtfLzxhyNG1D/N3n4ec9OjcfeyKTh5ehISIn/9eRyXFY8vd9Xi4c+KcM3rWzA1uQQ3nzQJR0+IHdb7kMVmx7Wvb8GnO6rx11OnDGtIohJyMqLx0JkzcP1b+bjtg23455kzBv16pJS4Y/V2hAVrcNsp7mfAvPXkSTjx0e/x+Ne7+02YUlRtRGKkzuPD3TNiwlBca/TIsVtMFlzwwkbUt3Xh9T8cgQnDCEqzEiMQHqzB5v2NWOGFXgmbXeKnkgacNM29MgDu+G32GHyxsxoPf1aM+AgdXt9Yho37GjEuLgzPXZCDJZPjfX6TbDDUKoFlM5Lxwvp9yE73TNbR6FAttGoVqltH1hw5i82OR74oxuSkSJwyzTvF3ImU4nYgJ4RIBfA/AIkA7ABWSSkfO2SbYwGsAbDPueh9KeU9Ax07JlyLF9bvw9aKZjz5u+xhz8fqS+7+Rlz35hbMTInCE+fOhlqBXqboMC0uO2osLj0yE4XVRiRHhfjlMIyBxIQH4+WL5+GM/2xAcY0Rz/x+jts1d3ozO82AlXPT8NKGUvw2OwVTkj03RyWQ/LS3AQeaO/CXk/v/MPu7+el4N+8A7v1oJ46dGO+VeZfFNUZ8vNURrOyu7T9Y6U1qdCjuXTEVN7xVgKe/Lem3yLO/k1KiuKatO3jb4/p5ZETjb6dNxUnTEvv8eQghcPwUR53KNfkH8MgXxbjwhY04Ymw0bj5p0kHzIgdittpx7Rt5+GxHDe48dQou8VEQ57Ji9hjsrW/H41/txti4sEHXrXs/7wB+2deIB86YjphhpIU/1ISECJwzNxWv/rwfFy3MQHpM7/OGiqqNQy4gPhwZsWH4qrAGVptd0RsYHWYbLnl5E/bWteOFi+Zi1jB7ZNQqgdlpUcjd36xY2/qzs7IVrZ3WIZWJUJoQAvefPh0nPvoDrn8rH7HhWty3YhpWzk0NuJtMlx2VCbuUyMnwTCCnUgnERwaPuB65dzZXoKzRhBcuyvFqaRAiJSjRI2cFcJOUMk8IEQEgVwjxhZRy5yHb/SClPHUoB07Wh+Cuc2fjL+9txalP/IDHz52t2Dj6kro2XPryZiTpdXj+whyEaJWtZyaEwOSkwA5K0mJCsfaaReiy2j0ycdrlLydl4fMd1bhj9Ta8e+VCvpECeGdzOSJ0GpwwJaHf7dQqgftPn4bTnlyPBz8txANnTPdYm37YXYe/fbhzSMFKX06fnYJvi+rw2Fe7sWh8bHfq80DgynC4bmsVPt5WhZK6dggBzMuIxgXLp+KkqYmIH8LPQ60SOCM7BafOSMYbG8vwxNd7cMbTG3D8lAT83wlZAwYYZqsdf3w9D1/srMFdy6bg4kW+DeJcblg6AaX17Xjo0yJkxIThlOn93+luNpnx93W7kJ0WhXNyUhVsx0Ssya/EQ58W4anzsg9bb7XZsaeuzSujJjJjQ2GxSVQ2dyo2lNBsteOq13KxpawJT/0u2+3XkZ1mwBNf74ax04IIDxc63lBSDwBY4MYQWiXEhAdj1QVzsGlfI847Ij1gh/knR4Xg7tMGX6pjOJL0uhE1R67TYsPjX+1GdloUjsuK93VziIbM7dtNUsoqKWWe87kRwC4Aio3JWDYzGWuvWQR9SBB+/9wveOa7ErczatUaO3HhCxuhFgIvXzJPkTu/I1VyVIhHgzjAUdD01lMmI6+sGW8rXIcqELV2WvDpjmqcNjN5UAXTpybrcemRmXhjYxk2l3qmttybG8tw0YubAAD3LJ+KX25dgrevWIALF2YMOYhzuXfFNCTpdbj+rS0wdnpnPqy7LDY7Ln15M0569Ac8+c0exEUE497lU/HLbUvw1hULcMGCjCEFcT1pNSpcuDAD3/35WPzfCRPxc0kDTnrse9z4Vn6fWfV6BnF3+1EQBzhuZj105gxkp0XhhrfykV/e3O/2//i0CC0dFtyncC2x+Egd/nDUWHy8rQp5ZYen1i9tMMFstQ+5jMpwZMQon7nyvo934tuiOvz99Ok4eYBgeTByMgywSwz4+1LC+pIGjI8PH/bfjJKy0wy44phxARvEeUuiPgTVrSMnkHv15/2obu3E/52YFRBDaIkOpei4ASFEBoDZAH7pZfUCIUSBEOITIcSQbhmNj4/AmmuOxMnTk/DgJ4W44pVctA7zg197lxWXvLQJDW1mvHDR3D6H2pB3/TZ7DOZlROPBTwvR2D78bHcjwcdbq9BpseOsIfRKXL90AsZEheD2D7bDomAWRCkl/vlZIW55fxsWjY/FB1cvdCtY6SlSF4RHz5mFA00duGvNDgVa61lSSty5Zge+LqzFjcdPxC+3LcWbly/A+QsyEB+h3AfRsGANrlk8AT/85ThcfrQjAFn8r29x99odqDP+OjfFbLXj6tccQdzfTpuKi/woiHPRBamx6oIcxEUE47KXN+NAc+9DsnL3N+GNjWW4eGGGR4ZXX370WMRFBOP+j3cddiOw2JnoZJIXhla6booplfBkR2ULXnEOGz13njJ1r2alRkEIzxcGN1vt2LSv0etlB8g9SXodqls6R0TZoLYuK/7zbQmOHB/r9aypREpRLJATQoQDeA/A9VLK1kNW5wFIl1LOBPAEgNX9HOdyIcRmIcTmurq67uXhwRo8ee5s3HnqFHxdWIvTnvgRu6oOPU3/LDbHB59dVUY8dd5szPRQZicaOiEE7l0xDW2dVjz4yS5fN8en3s2twPj4cMxM0Q96n1CtBvcsn4qiGiOe/WGvIu3ostpw/Vv5eOqbEpw7LxXPX5ij+FCrnIxoXLt4At7fcgBrBlm82VdeXF+KNzaW4cpjxuFPSyYgLsKzPflRoVrcevJkfPfn43DmnFS88vN+HPPPb/Cvz4vQ0NaFq1/LxZe7anDP8qm4cGGGR9vijtjwYLxw0Vx0WWy49KVNh/W+ukqbJEbqcP3xEz3ShrBgDW5YOhGb9zfhsx01B60rrDZCJYDx8Z4rBu4SFxGMMK0a+xQI5KSUuOfDnTCEanHDUuV+bhG6IGQlRHg8kMsvb0aHxYYF/AAdUBIideiy2tFsCoxRFP158cd9aGg34/9OzPJ1U4iGTZFATggRBEcQ95qU8v1D10spW6WUbc7n6wAECSF6ffeWUq6SUuZIKXPi4g4uyi2EwCVHZuLNy49Ah8WG059ej/fzKgbVRiklbv9gG74rrsN9K6Zh8aT+5x6R92UlRuDSIzPx9uYKjw0R9Hd769qQu78JZ/VSO24gSyYn4KSpiXjsy90oa+i7wO1gNJvMOP+5jViTX4mbT8rC/adP91jm2GsXj8ecdAPu+GB7v4V5fembolrc9/FOnDAlATd7+Z9+ol6HB86Yji9vPAZLJifgia/3YN79X+HLXbW4d8U0XKBgJllPmZgQgafOy8bu2jb86Y0tB9XOe2lDKXZVteLu06Z4dFjb2TkpGB8fjn98WnhQr3VxtREZMWGDGsbsLiEE0mPCUKrA0MpPt1fjl32NuPH4iYonOcrJMGBLWfNhtdWUtKGkHkL4fn4cDc1IqSXXbDJj1Q97cfyUhGEnByLyB25/MhOOT5vPA9glpXykj20SndtBCDHPed6G4Z4zJyMaH117FGalRuHGtwtw+wfb0GW19bvPo1/uxtubK/CnxeMVG4JCyvvTkglI1utwx2plhgh2WmxYt60KL63fFxBDQd511o47fZipv+8+bSqC1CrcsWb7sF9vWYMJZ/xnA/LLm/H4ubNx9bHjPTp3QKNW4dFzZgEArn8rX/EC2e4qrjHi2te3ICsxEv8+Z5bPkvFkxobhiXNn46Nrj8RJ0xLx0JkzcP4R6T5py3AcPTEOd582Fd8U1eG+jx297lUtHfj3F8U4LisOJ3ogK25PGrUKt5w0Cfvq2/HGxrLu5UU1Rq/Mj3PJjA1zu0eu02LD39ftwqTECKycq1xiGJc56Qa0dVm7C6V7woaSBkxL1nsl0y4px5U9vCbA58n99/u9aOuy4qYTPDMKgMhblLjFvgjA+QAWCyHynY9ThBBXCiGudG5zJoDtQogCAI8DWCnd/FQd56zPduUx4/DaL2U465mfUNHU+938tzaV4bGvduOsOSm4wUNDd0gZYcEa3LlsKgqrjXh5Q+mwjmGx2fFNYS1ueCsfc+79Ale/loe7P9yJ134pG3hnH7LZJd7PO4BjJsYNew5aol6H/zthIr4vrsNHW6uGvP+Wsiac/vR6NLab8epl83HazORhtWOoUqNDcd/p05C7vwlPfrPHK+ccjIa2Llz68iaEaNV4/sIchPlBIoRpY/R46nfZOFvBzI7ecv4R6bhkUSZe2lCK//1Uins+3AmrXeJvpw2vhudQLZkcj/mZ0XjsS0dWxk6LDaUN7V4pPeCSERuKiqYOt25UPf/jPlQ0deDOU6d4JEV+Tno0ACC3l+QwSjCZrdhS1oSF49kbF2gSIwO/R67W2ImX1pdi2YxkTEoM7OziREpkrfxRSimklDOklLOcj3VSymeklM84t3lSSjlVSjlTSnmElHKD+0133mE9eRL+e/4c7Ktrx6lP/IjviusO2uabwlrc9sF2HD0xDvefMZ1ZiQLAiVMdtbX+/UUxqgZZr8Zul/h5bwNu+2Ab5v39S1z80iZ8tasGp85IxmuXzcfRE+Nw38c7UVLX5uHWD9+Pe+pR3do5pCQnvTl/QQZmpOhxz0c70dIx+HkMn26vxspVPyMsWIP3rlqIeZnRbrVjqJbPGoPTZ4/B41/tRu5+3w+t7bLacOWruaht7cKzF+QgOar3wuw0NLf/ZjKWTIrHXWt34JPt1fjTkgmKpeIfiBACt/9mMhrazfjvd3uxu6YNUsK7gVxMGGx2OexhxDWtnXjqmz04cWoCFnqo/lqKIQRxEcHI89A8uc2lTbDYJBNMBKC4iGCoBAK6ltzT35TAbLPzxj6NCIFV7bIPJ05NxIfXHonESB0uenEjHv2yGHa7xNaKZlz9Wh4mJ0Xg6fOyPTbHh5QlhMDdy6bCandM5u+LlI7f8X0f7cTCB7/GylU/44O8AzhyQhyevSAHm+5Yin+cOQOLxsfin2fOgC5IjRveylc0q6OS3s2tQFRoEJZMdq+WjaO23HQ0tHXhn58VDmqf53/ch6tey8XkpEi8f/VCjIvzfOKH3tyzfCrGGEJw3Zv5w85MqwQpJW57fzs2lTbh4bNmcg6FgtQqgcfPnY2pyZGYlBiBPxw11qvnn5EShdNmJuO5H/fi+92OG3/eDOTGxjkzVw5zntw/Pi2E1SZx+ylTlGzWQYQQmJNmwGYP3VBZX1KPILXAXA8VribPCVKrEBcRHLA9cgeaO/D6L2U4a06Kx0srEXmD78cJKSQjNgwfXL0It6/ehke/3I3c/U3YVdWKmHAtXrhoLmvDBJi0mFBcu3g8Hv68GN8U1R5UqHNPrRFr8yuxtqASpQ0mBKkFjpkYh1tPmYTjpyQgVHv47zohUocHTp+Oq17LwxNf7caNJ/hXlqoWkwWf7ajGuXNTEaxxP+nCtDF6XLwoEy+s34czslOQndb7ByabXeLej3bipQ2lOGlqIh5dOcsrSR/6EqELwqPnzMbZ//0Jd67ejkdXzvZJO/77/V68l1eB65ZMwDIvDS8dTcKCNVh99SJY7RJajfdvsP35xCx8ur0aj325G1qNCunR3ukRBHrUkqsfeo9cfnkz3s87gCuPGefxXsycDAM+3VGN2tZOxeu8/VTSgNmphl7fq8n/BXItuce/3A0AuHbJBB+3hEgZI6qLKkSrxr/Omon7T5+OX/Y2wmqXePmSeYrWeCLv+cPRYzE2Lgx3rdmBkro2/OfbEpz82A9Y+sj3eOKbPUiOCsGDZ0zHptuX4rkL52L5rDH9fjA4eXoSfpudgie/2ePx1NpD9eHWSpitdpw5R7l5TzcePxGJkTrc9v62XnshTWYrrnglFy9tKMVlR2bi6fOyfRrEucxJN+BPiydgdX4lVm/xfkmCz3ZU4x+fFuLUGUm4fin/2XuKRq3y2fWWGh2KCxemw2yzY0J8uEfmmfUlOkyLCJ1myLXkHOUGdiAuIhjXLB7vodb9KjvdcfNH6ffKFpMF2w+0YAHrxwWsxMhgVAdgj9zeuja8m1eB845IwxgOlacRYkQFcoBjSMjv5qfhk+uPwpo/LvLZEDFyX7BGjfuWT0NZowlL/vUd/vFpIYI1Ktx56hT8cusSvP6HI7ByXhqiQrWDPubdp01BclQIbnw7H+1dVg+2fmjeya3ApMQITBuj3MTrsGAN/naaI3HMCz/uO2hdrbETK1f9jK8LHcWk7zh1is+yMfbmj8eNQ066AXes9m5Jgh2VLbj+zXzMGKPHw2fN5JzaEeya4ybAEBqEGUOo16gEIQQyY4degmBtQSXyyprx5xOzvDLCZFqyHlqNSvFA7ud9DbBLYJGH5veR5yXpQwIykPv3l7uhVatw9bGevxFC5C0jLpBzGRcXjvQYjn8OdAvHx+LOU6fgzydm4fs/H4fVf1yES47MHPZQnwhdEB45exbKGk2496O+59950+4aIwrKm3HmMGrHDeSEqYk4YUoC/v1lcXdAtKfWiDOe3oDdNW1YdX6OXxaT1qhV+Pc5syAAXPfmFq+UJKg1duKylzcjKjQIz16Q4xe9k+Q5+tAgfHLd0bjtlMleP3dGzNBKEJjMVjywrhDTx+hxZnaKB1v2K61GhZkpesUzV/5U0gBdkIrzTgNYol4HY5cVbX50M3QgOytb8WFBJS5elIG4iGBfN4dIMSM2kKOR45IjM/HH48YrNidkXmY0rjxmHN7cVI4vdtYockx3vJtXAY1KYMUwa8cN5O7TpkItBO5csx0bSupxxtMb0Gmx460rjsDSKQkeOacSXCUJ8sqa8cTXni1J0Gmx4Q//y0WzyYJnL8hRfE4Q+adEvQ4ROu/XMcuIDUNlc8eA9U9dnvluL6pbO3HnMu/2nGenG7D9QAs6LYNr52Cs31OPuRnRPpkbScpwFQUPpF65R74oQoROgyuOHufrphApiu+kNCrdsHQipiZH4pb3tqLO2OWzdlhtdryfdwDHZsUjNtwzdwmTo0Jw4wlZ+KaoDr9/7hfER+rwwdULMSMlyiPnU9LyWWNwxuwxeOLr3dhc6pkMelJK/PndrSgob8ajK2dh2hjvDrWj0SczNhR2iUENGz7Q3IH/fleCZTOTMTfDuyVBctKjYbFJbDvQosjxao2d2F3bxmGVAS4hMrACubyyJny5qxZXHD2WBehpxGEgR6OSVqPCo+fMQluXFX95byvcrE8/bD/srkedsQtn5Xh2uNSFC9IxPzMaR02Iw3tXLUSqF7P0uetvy6cixRDqsZIEj3+1Bx8WVOLmk7Jw4tRExY9PdKihZK58YN0uCAHccvIkTzfrMNlpUQAcdd+U8FNJAwBgIROdBDRXj9xg67z62sOfFSEmTIuLF2X6uilEimMgR6PWhIQI3HLyJHxdWIvXN5b5pA3v5JYjOkx7UHkFT9CoVXjz8iPw8iXzoA8JrDuSEbogPLZyFqpbO3HHB9sVDbo/2lqJf39ZjDOyx+CqYzjkhrzDVb9qoMyVm0ob8dHWKlx+9DifZNmLCQ/G2NgwxRKebNjTgEidBlOT2esdyFw9cjUBUIJgw556bChpwNXHjUcYy1DRCMRAjka1Cxdk4KgJsbjvo13YW9fm1XM3tZvx5c5arJg1xivzRQI5A+PsNAOuXzIBawsq8eovZYrM2Skob8ZNbxcgJ92AB86YHtA/HwosUaFaRIUGYW8/gZzdLvG3D3cgMVKHK4/xbtH0nrLTDcgra1LkBsqGvfU4YmwM1H6UIZeGThekhiE0yO+Lgksp8c/Pi5Ck1+G8+Wm+bg6RRzCQo1FNpRL455kzodWocMPbBb3WW/OUtQWVMNvsOHOOd7LQBbqrjxuPeRnR+Ovq7Zhy56dY/PC3uOrVXPz7i2J8sq0KJXVtsNkH92GzqqUDf/jfZsRFBOO/589RpAg70VBkxIT12yP3bm4Fth9oxa2nTPJp4ew56QY0tpuHlGWzN+WNJpQ3dnB+3AiRGAAlCL7aVYstZc3405IJzEJMIxb7mWnUS9TrcP/p0/HH1/Pw5Nd7cMPxE71y3ndzKzAlKRJTkpWrHTeSqVUCL10yF98U1qGoxoii6lYUVhvx6Y5quDoLgjUqjI8PR1ZiBLISIhxfEyOQGKnr7nEzma247OXNMJltePWy+YjxUJIZov5kxobh570Nva4zdlrw0GdFyE6Lwmkzk73csoPl9CgMPtaNuqwbSuoBcH7cSJGk1/l1j5zdLvHw50VIjwnlzVIa0RjIEQH4zYwkfLVrDJ78Zg+OzYrD7DSDR89XWN2KbQdacNeyKR49z0gTqtXgNzOS8BskdS/rMNuwp7atO7grqmnD+j31eD/vQPc2kTpNd1BXWm/CrqpWPH/RXExMiPDFyyBCRkwYPthyAB1mG0K0B/cWPPVNCerbuvD8hTk+H/I7Li4ckToNcvc34ayc1GEfZ0NJA+IigjE+fvjBIPmPRL0OBeXNvm5Gnz7aVoXCaiMePWcWgtQcfEYjFwM5Iqe7l0/FL/saccNb+fj4T0d5dGL0u5srEKQWWD7LM7XjRpMQrRrTU/SYnnJwAoVmkxlF1UYU1xhR6Py6Jr8S7V1W/PXUKR5PMEPUn8w4R8KT/Y3tmJT4a6/8/oZ2vPDjPvw2OwUz/aBotkolkJ1ucCvhiZQSG0oasHBcjM8DU1JGYqQODe1mdFpsfjds0Wqz49EvipGVEIFlPu7RJvI0BnJETpG6IDxy9kysfPZn3PfxLjxwxnSPnMdis2N1/gEsmZSA6DCtR85BjoQS88fGYP7YX4dySSnRYbH5dM4REQBkxvyaubJnIPf3j3dBoxa4+aQsXzXtMDnpBnxbVIdmkxlRoUN/z9pT24Y6YxeHVY4gic4SBLWtXUiL8a9yNu/nHcDe+nb89/w5TKxDIx77m4l6mD82BpcfPRZvbCzDlztrPHKOb4vqUN9m5rh9HxBCMIgjv5AR6/jw27OW3IY99fh8Zw3+eNz47hTv/iDbOU9uS1nzsPZfv8c1P46JTkYKf64lt6bgACbEh+OEKQm+bgqRxzGQIzrEjcdPxOSkSNzy/lbUt3Upfvx3c8sRGx6MY7LiFD82EQWGCF0QYsO13ZkrrTY77vloJ1IMIbj0SP8qXDwrNQpqlRj28MoNJQ1IjQ5BarR/9dzQ8LkCuWo/rCW3v8GEyUmRHMZLowIDOaJDBGvUeGzlLLR2WnHLe1sVLUDd0NaFr3bV4vTZyZyATTTKZcSEYV+DI5B7Y1M5CquNuP2UyX435yhUq8GUpEhs3t845H1tdomf9zZgEXvjRhRXj7G/lSCw2uyoaulEGm8a0CjBT5JEvZiYEIG/nDQJX+6qxRsbyxU77pr8SljtEmfOGX72NyIaGTJiHbXkWkwWPPJ5EeZnRuOkaYm+blav5qQbUFDeMuRamzsqW9DaacUCzo8bUSJ0QQgP1vhdCYKqlk7Y7BKp0SG+bgqRVzCQI+rDxQszsGh8DO79aKfbxXBd3smtwIwUPbISmfaeaLTLjA1DrbELD366C80dFty5bIrfDgebk25Ah8WGXVWtQ9pvQ4mjVh4DuZEnUa9DjZ8NrSxrdMw5TTWwR45GBwZyRH1QqQQePmsmgtQCN7yVP+Q70YfaUdmCXVWtOItJTogIjqGVAPDGxnKsnJuGqcn6AfbwnTk9CoMPxfo99ZiYEI74CP9J3kLKSIz0v6Lg5a5AjkMraZRgIEfUjyR9CP5++nTklzdjxt2fY/mTP+Lmdwvw3A978ePuetQaOwc9h+6dzRXQqlWsa0NEAH7NXBkRrMFNJ0z0cWv6lxwVgmS9bkiBnNlqx6bSRmarHKES9Tq/myNX3mSCWiW6k7EQjXTMw000gGUzk6HVqPDL3kYU1bTi68I6vL25ont9dJgWExPCMSkxEhMTIpCVGIGJCeGI0AV1b2O22rEm/wCOn5owrDpMRDTyjI0NR2x4MK5bMh6x4cG+bs6AhloYPL+8GZ0WO+vHjVBJeh1qjZ2w2uzQ+EnyrrLGDiRH6fymPUSexkCOaBBOnJqIE6f+moSgoa0LRTVGFFUbUVxjRGG1Ee9sLke72da9zZioEGdQFwEpJZpMFtaOI6JuIVo1Nt62BKoAKVo8J92Aj7ZWobK5A8lRAyeTWL+nHirhqM9JI0+iXge7BOrbzN0Fwn2tvNHE+XE0qigSyAkhTgLwGAA1gOeklA8esl44158CwATgIillnhLnJvKFmPBgLAwPPmjIkJQSFU0d3YFdsTPQ+2F3HSw2iSS9DkdPYO04IvpVoARxAJCTHg0A2Ly/CacNIpD7qaQB08booQ8JGnBbCjyJkb8WBfeXQK6iyYSlk1kInEYPtwM5IYQawFMAjgdQAWCTEGKtlHJnj81OBjDB+ZgP4D/Or0QjhhACqdGhSI0OxZIe/0gsNjv21bcjQqeBOoA+tBER9TQpKQIhQWrk7W/CaQPM9TWZrdhS3oRLjxzrpdaRt7mCN3+ZJ9feZUV9m5mJTmhUUWIQ8TwAe6SUe6WUZgBvAlh+yDbLAfxPOvwMIEoIkaTAuYn8XpBahYkJEUjSs64NEQWuILUKs1KjBjVPblNpEyw2yflxI5jrf5q/ZK6saOoAAKQY+L+WRg8lArkxAHpWTK5wLhvqNgAAIcTlQojNQojNdXV1CjSPiIiIlDAn3YCdVa1o77L2u92GknoEqQXmZkR7qWXkbYbQIGg1Kr+pJecqPZDGHjkaRZQI5HobK3ZoPvbBbONYKOUqKWWOlDInLo7ziYiIiPzFnHQDbHaJgormfrfbsKcBs9MMCNGqvdMw8johhF/VkitvYg05Gn2UCOQqAKT2+D4FQOUwtiEiIiI/lp3mLAxe2vfwyhaTBdsrW7CI9eNGPH+qJVfWaEJIkBoxYSzxQ6OHEoHcJgAThBCZQggtgJUA1h6yzVoAFwiHIwC0SCmrFDg3EREReYk+NAgT4sORW9Z3IPfT3gZICSwcz/lxI12SXodqvxla2YG06FA4EqUTjQ5uZ62UUlqFENcA+AyO8gMvSCl3CCGudK5/BsA6OEoP7IGj/MDF7p6XiIiIvG9OugHrtlXBbpe9lk/4qaQeIUFqzEyJ8n7jyKsSIx09clJKnwdQFU0mpEYz0QmNLorUkZNSroMjWOu57JkezyWAPypxLiIiIvKdOekGvLmpHHvq2jAxIeKw9RtKGjAvMxpajRKDfsifJep1MNvsaGw3IyY82GftkFKirNGEI1h8nkYZvssSERHRoM1Jd86T66UMQW1rJ3bXtrHswCiRpHcVBfft8MrGdjNMZhszVtKow0COiIiIBi0zNgzRYVps7iXhyU97GwAAi8Yz0clokOisJefrEgTlzhpyzFhJow0DOSIiIho0IQSy0wzI6yXhyfo99dCHBGFyUqQPWkbelhjpHz1yZY2u0gOcI0ejCwM5IiIiGpI56Qbsq29HQ1vXQcs3lDRgwdgYqHtJgkIjT1xEMNQq4fMSBK5i4KkG9sjR6MJAjoiIiIYkJ+PweXLljSZUNHWw7MAoolYJxEcE+7xHrqLJhJgwLcKCFcnhRxQwGMgRERHRkEwfo0eQWhxUT279nnoAYKKTUSZRr/P9HLnGDqRwfhyNQgzkiIiIaEh0QWpMG6NHbo+EJxtKGhAfEYxxceE+bBl5W2KkDlUtHT5tQ1mjCakGzo+j0YeBHBEREQ3ZnDQDth5oQZfVBiklNpQ0YOG4GJ8XhibvStTrUOUsCu4LNrtEZXMHSw/QqMRAjoiIiIZsTroBZqsdOypbsbu2DfVtXVjIsgOjTpJeB5PZBmOX1Sfnr2rpgNUuWXqARiXOCiUiIqIh6y4MXtqEILWjF47z40afBGcJgpqWTkTqgrx+/rIRnLHSV72cFDgYyBEREdGQxUfqkBodgtz9TbBJibToUKSMwA/T1L8kZ1HwqpZOTEiI8Pr5Kxod8/NG2tDK8PBwtLa2QgiB4OBgXzeH/BSHVhIREdGw5KRHY/P+Rvy8twGLWHZgVErSO3rkfFVLrrzJBJUAkqJ0Pjm/p0RERCAjIwMqlQpGoxF2u93XTSI/xECOiIiIhiU73YD6NjOMnVYsGMf5caNRfKSjt6jaRyUIyhpNSNKHIEg98j7SBgcHIy0tDfHx8Whvb0dnp2/LPJD/GXlXPREREXnFnDRD9/MFY9kjNxoFa9SICdP6rCh4eaNpxA2r7EkIAYPBgIyMDKjVavbO0UEYyBEREdGwZCVGIDxYg6yECMRFcB7PaJWo16HaR7Xkyps6kBo98mvIsXeOesNkJ0RERDQsapXA/50wEXERI2t+Eg1Nkl6HiibvB3IdZhvqjF0jMmNlb1y9c2FhYaiurobRaERYWBhUKvbLjFYM5IiIiGjYLlqU6esmkI8lROqQu7/J6+etaHKWHhjBQyt7o9VqkZqaiubmZtTV1UGj0UCn482U0YghPBERERENW5JehyaTBZ0Wm1fPWz5KAzng4LlzQUFBnDs3SjGQIyIiIqJhS3TWkvN2CYJyZw250TBHri9arRYpKSlITExEe3s7Ojp8M1eRfIOBHBERERENm6uWnLczV5Y1mqALUiEufHQn2hFCQK/XIzMzE1qtlr1zowgDOSIiIiIatoRIRyBX4+VacuWNJqQaQiGE8Op5/VXP3jmTyQSTyQSr1errZpEHMdkJEREREQ1boo965BylB0bf/Lj+uHrnQkND0dLSgvb2drS1tUFKCSEEgoKCoNVqGfyOEAzkiIiIiGjYwoM1iNBpvFpLTkqJ8kYT5mUYBt54FAoKCkJsbCxiY2Nht9thNpvR1dWFtrY2mEwm2O12CCGgVquh1WqhVqt93WQaBrcCOSHEPwEsA2AGUALgYillcy/blQIwArABsEopc9w5LxERERH5jyS9DtVeHFrZbLKgrcvKHrlBUKlU0Ol00Ol00Ov1kFLCYrHAbDbDZDKhra2tO0mKSqVCUFAQNBoNe+0CgLs9cl8AuFVKaRVC/APArQD+0se2x0kp6908HxERERH5mYRInVezVo7m0gPuEkJAq9VCq9UiPDwc8fHxsFqtsFgs6OzsRFtbG9rb2w/bT6VSdffiqVSq7gf5jluBnJTy8x7f/gzgTPeaQ0RERESBJkmvQ1G10WvnK2t0BnIGBnJK0Gg00Gg0CAkJgcFggN1uh8Vigd1uh91uh81m6w72rFYrrFYrzGZzv8lUhBAHBXpCiO5evsE+p/4pOUfuEgBv9bFOAvhcCCEB/FdKuUrB8xIRERGRDyXqQ1DX1gWLzY4gted7aVhDzrNUKhWCgwdX1sEV7PUM+ux2O6xWa/eyntv1fC6lhJQSNpsNUsrDtnOHK8GLJ7na7KLk+VzHdn6VvW0zYCAnhPgSQGIvq26XUq5xbnM7ACuA1/o4zCIpZaUQIh7AF0KIQinl932c73IAlwNAWlraQM0jIiIiIh9LjNRBSqDO2IXkKM8HV+VNJhhCgxChC/L4uah/nhxieWigdOj3fS3rb7mneDho7DWyHTCQk1Iu7W+9EOJCAKcCWCL7+IlJKSudX2uFEB8AmAeg10DO2Vu3CgBycnK8+xsgIiIioiHrWRTcK4Fco4nz40aBQ4MjDrc8mFvhsxDiJDiSm5wmpTT1sU2YECLC9RzACQC2u3NeIiIiIvIfrlpy3kp44ioGTjSaudsP+iSACDiGS+YLIZ4BACFEshBinXObBAA/CiEKAGwE8LGU8lM3z0tEREREfuLXHjnP15Kz2SUONLMYOJG7WSvH97G8EsApzud7Acx05zxERERE5L/0IUEI1qhQ44VacjWtnbDYJBOd0KjH4g9ERERE5BYhBJL0OlR5YWglSw8QOTCQIyIiIiK3Jeq9UxS83BnIpXFoJY1yDOSIiIiIyG1J+hBUe2FoZXlTB4SAV7JjEvkzBnJERERE5LbkKEePXKfF5tHzlDeakBSpg1bDj7E0uvEvgIiIiIjcNiMlCla7xI7KFo+ehzXkiBwYyBERERGR27LTDACAvP3NHj1PeRMDOSKAgRwRERERKSAuIhhp0aHI3d/ksXN0Wmyoae1ixkoiMJAjIiIiIoVkp0Uht6wJUkqPHL+iyVFwPC2GiU6IGMgRERERkSLmpBtQZ+zqDriUVt7EGnJELgzkiIiIiEgRs13z5Mo8M7yywlUMnHPkiBjIEREREZEyJiVGIFSrRp6H5smVNZqg1agQFx7skeMTBRIGckRERESkCI1ahVmpjnlynlDe2IFUQwhUKuGR4xMFEgZyRERERKSY7DQDdlUZYTJbFT82Sw8Q/YqBHBEREREpZk66ATa7REG58oXByxpNTHRC5MRAjoiIiIgUMzstCoDyCU9aTBYYO61IY48cEQAGckRERESkoKhQLcbFhSme8KS79EA0a8gRAQzkiIiIiEhh2WkG5ClcGLzMWXoghUMriQAwkCMiIiIihc1JN6DJZMG++nbFjlnuDOTSYhjIEQEM5IiIiIhIYdnpjsLguQoOryxvMkEfEoRIXZBixyQKZAzkiIiIiEhR4+PCEaHTIK+sWbFjljd2cH4cUQ8M5IiIiIhIUSqVcMyTU7JHrtHEjJVEPTCQIyIiIiLFZacZUFxrRGunxe1j2e0SFU0drCFH1AMDOSIiIiJS3Jx0A6QE8hUYXllr7ILZZkcKe+SIujGQIyIiIiLFzUzVQwhlCoO7Sg+kGjhHjsjFrUBOCHG3EOKAECLf+Tilj+1OEkIUCSH2CCFuceecREREROT/InRByEqIUCRzZXfpAfbIEXVTokfu31LKWc7HukNXCiHUAJ4CcDKAKQDOFUJMUeC8REREROTHstMNyC9rht3uXmHw8iYThADGsEeOqJs3hlbOA7BHSrlXSmkG8CaA5V44LxERERH50Jw0A4xdVuyubXPrOGWNJiRE6BCsUSvUMqLAp0Qgd40QYqsQ4gUhhKGX9WMAlPf4vsK5rFdCiMuFEJuFEJvr6uoUaB4RERER+cIchQqDVzR2cFgl0SEGDOSEEF8KIbb38lgO4D8AxgGYBaAKwL96O0Qvy/rsX5dSrpJS5kgpc+Li4gb3KoiIiIjI76THhCI6TOt2wpPyJhNSWAyc6CCagTaQUi4dzIGEEM8C+KiXVRUAUnt8nwKgclCtIyIiIqKAJYT7hcG7rDZUt3ayhhzRIdzNWpnU49vTAWzvZbNNACYIITKFEFoAKwGsdee8RERERBQYstOjsLe+HY3t5mHtf6CpA1IyYyXRodydI/eQEGKbEGIrgOMA3AAAQohkIcQ6AJBSWgFcA+AzALsAvC2l3OHmeYmIiIgoAMxJc8yT2zLM4ZXlTR0AgFQGckQHGXBoZX+klOf3sbwSwCk9vl8H4LDSBEREREQ0ss1IiYJGJZBX1oQlkxOGvL+rhlwq58gRHcQb5QeIiIiIaJQK0aoxJTly2JkryxtN0KpVSIjQKdwyosDGQI6IiIiIPCo7zYCC8hZYbfYh71veZEKKIQQqVW+J0IlGLwZyRERERORR2ekGdFhsKKw2Dnnf8sYOpHB+HNFhGMgRERERkUe5Uxi8rNGEVAPnxxEdioEcEREREXlUsl6HhMjgIRcGb+20oKXDwtIDRL1gIEdEREREHiWEwJx0w5B75H7NWMlAjuhQDOSIiIiIyOOy0wyoaOpAbWvnoPfpDuQMDOSIDsVAjoiIiIg8Lts5T24owyvLGx3FwDm0kuhwDOSIiIiIyOOmJkdCq1YNaXhleZMJEToN9KFBHmwZUWBiIEdEREREHhesUWN6ih55Zc2D3qe80cRhlUR9YCBHRERERF6RnRaFbRUt6LLaBrV9WaOJwyqJ+sBAjoiIiIi8Yk66AWabHTsqWwfcVkqJiqYOpEazhhxRbxjIEREREZFXZKc5E54MYp5cnbELXVY7Sw8Q9YGBHBERERF5RXykDimGkEFlrixj6QGifjGQIyIiIiKvcRUGl1L2u115E4uBE/WHgRwREREReU12mgE1rV2obOm/MLirhlyKgXPkiHrDQI6IiIiIvGaOszD4QPXkyhpNiI8Ihi5I7Y1mEQUcBnJERERE5DWTEiMQEqQeMOFJOUsPEPWLgRwREREReY1GrcLMVP2ACU8cpQcYyBH1hYEcEREREXnVnHQDdla2osPce2Fws9WOqpYOpHJ+HFGfGMgRERERkVdlpxlgtUtsrWjudX1lcwfskhkrifrDQI6IiIiIvGq2szB4bh/DK1l6gGhgDOSIiIiIyKuiw7QYGxuGvP3Nva53lR5gIEfUN407Owsh3gKQ5fw2CkCzlHJWL9uVAjACsAGwSilz3DkvEREREQW27HQDvi6shZQSQoiD1pU1mhCkFkiM1PmodUT+z61ATkp5juu5EOJfAFr62fw4KWW9O+cjIiIiopEhO82Ad3MrUNpgQmZs2EHryptMGBMVArVK9LE3ESkytFI4bqOcDeANJY5HRERERCObqzB4b/XkKhpNHFZJNACl5sgdBaBGSrm7j/USwOdCiFwhxOUKnZOIiIiIAtSE+HBEBGt6TXhS1mhCioGBHFF/BhxaKYT4EkBiL6tul1KucT4/F/33xi2SUlYKIeIBfCGEKJRSft/H+S4HcDkApKWlDdQ8IiIiIgpAKpXArLSow3rk2rqsaDJZkMYeOaJ+DRjISSmX9rdeCKEBcAaAOf0co9L5tVYI8QGAeQB6DeSklKsArAKAnJwcOVD7iIiIiCgwzUk34LGvdsPYaUGELggAUN7oKj3AYuBE/VFiaOVSAIVSyoreVgohwoQQEa7nAE4AsF2B8xIRERFRAMtOM0BKoKD813x5Za5AjkMrifqlRCC3EocMqxRCJAsh1jm/TQDwoxCiAMBGAB9LKT9V4LxEREREFMBmpUVBCCC3x/BKV48ch1YS9c+t8gMAIKW8qJdllQBOcT7fC2Cmu+chIiIiopElUheEifERyOuR8KSiqQPhwRpEhQb5sGVE/k+prJVEREREREOWnW5AXlkT7HZHaoTyRhNSDCGHFQknooMxkCMiIiIin8lOi4Kx04o9dW0AHHPkOKySaGAM5IiIiIjIZ3oWBpdSoqKpg8XAiQaBgRwRERER+UxmbBgMoUHI3d+E+jYzOiw2pBpYeoBoIAzkiIiIiMhnhBDITnPMk3OVHkiLYY8c0UAYyBERERGRT2WnG1BS147tBxz15FhDjmhgDOSIiIiIyKey0xzz5NYWVAIAUhjIEQ2IgRwRERER+dTMVD3UKoHc/U2IDQ9GiFbt6yYR+T0GckRERETkU6FaDSYnRQAA0qKZ6IRoMBjIEREREZHPzXEOr2TpAaLBYSBHRERERD6X7awnx0QnRIPDQI6IiIiIfG5eZjS0ahWmJEf6uilEAUHj6wYQERERESXpQ7Dh1sWICdP6uilEAYGBHBERERH5hdjwYF83gShgcGglERERERFRgGEgR0REREREFGAYyBEREREREQUYBnJEREREREQBhoEcERERERFRgGEgR0REREREFGAYyBEREREREQUYBnJEREREREQBhoEcERERERFRgGEgR0REREREFGCElNLXbeiTEMIIoMjX7fAxPYAWXzfCD8QCqPd1I/wArwdeCy68FngtuPBacOD1wGvBhdeCA6+HkXMtpEsp4w5dqPFFS4agSEqZ4+tG+JIQYpWU8nJft8PXhBCbR/u1APB6AHgtuPBa4LXgwmvBgdcDrwUXXgsOvB5G/rXAoZX+70NfN4D8Cq8HcuG1QC68FsiF1wL1xOthhGMg5+eklPwjpG68HsiF1wK58FogF14L1BOvh5HP3wO5Vb5uAPkNXgvkwmuBXHgtUE+8HsiF1wK5jOhrwa+TnRAREREREdHh/L1HjoiIiIiIiA7h1UBOCPGCEKJWCLG9x7KZQoifhBDbhBAfCiEincu1QogXncsLhBDH9tjnHCHEViHEDiHEQ958DaQMIUSqEOIbIcQu5+/xOufyaCHEF0KI3c6vhh773CqE2COEKBJCnNhj+d+FEOVCiDZfvBZyj8LXwqfO94sdQohnhBBqX7wmGh6Fr4VvncvynY94X7wmGj6lrgchRESP6yBfCFEvhHjURy+LhkHh9wZ+hgxgQ70WhBAxzu3bhBBPHnKswP/8KKX02gPA0QCyAWzvsWwTgGOczy8BcK/z+R8BvOh8Hg8gF47AMwZAGYA457qXASzx5uvgQ5FrIQlAtvN5BIBiAFMAPATgFufyWwD8w/l8CoACAMEAMgGUAFA71x3hPF6br18XHz6/FiKdXwWA9wCs9PXr48Nn18K3AHJ8/Zr48I/r4ZDj5gI42tevjw/vXwv8DBn4j2FcC2EAjgRwJYAnDzlWwH9+9GqPnJTyewCNhyzOAvC98/kXAH7rfD4FwFfO/WoBNAPIATAWQLGUss653Zc99qEAIaWsklLmOZ8bAewCMAbAcjjeWOH8usL5fDmAN6WUXVLKfQD2AJjn3P9nKWWVF5tPClL4Wmh1bqMBoAXAScABRMlrgQKfJ64HIcQEOG4O/+DxF0CKUfBa4GfIADfUa0FK2S6l/BFAZy/HCvjPj/4wR247gNOcz88CkOp8XgBguRBCI4TIBDDHuW4PgElCiAwhhAaOX1QqKGAJITIAzAbwC4AE1x+V86trONQYAOU9dqtwLqMRRIlrQQjxGYBaAEYA73q+1eQJCr0vvOgcSvdXIYTwfKvJUxT8P3EugLek83Y8BR43rwV+hhxBBnktjGj+EMhdAuCPQohcOLpIzc7lL8Dxh7cZwKMANgCwSimbAFwF4C047qiVArB6t8mkFCFEOBxD4K7v0ZvS66a9LOM/4hFEqWtBSnkiHEMlggEsVrSR5BUKXQvnSSmnAzjK+Thf2VaStyj8f2IlgDeUaht5l7vXAj9DjhxDuBZGNJ8HclLKQinlCVLKOXC8uZY4l1ullDdIKWdJKZcDiAKw27nuQynlfCnlAgBFruUUWIQQQXD8Eb4mpXzfubhGCJHkXJ8ER88K4Ajqe941SwFQ6a22kmcpfS1IKTsBrIVjqAUFEKWuBSnlAedXI4DXwSGXAUnJ9wYhxEwAGillrscbTopT8L2BnyED3BCvhRHN54GcK5OYEEIF4A4Azzi/DxVChDmfHw9Hb9zOQ/YxALgawHM+aDq5wTnM6XkAu6SUj/RYtRbAhc7nFwJY02P5SiFEsHOo7QQAG73VXvIcpa4FIUR4jzdxDYBTABR64zWQMhS8FjRCiFjnMYMAnArHMH4KIB74P3Eu2BsXkJS8FvgZMrAN41oY2byZWQWON9AqABY47pZcCuA6ODLOFAN4EL8WKc+A407JLjgmo6Yfcpydzgez0gXgA44MQhLAVgD5zscpcGSU+gqOO2RfAYjusc/tcPTYFgE4ucfyh5zXk9359W5fvz4+vH8tAEiAIwvuVgA7ADwBx913n79GPrx+LYTBkZnQdS08hl6yF/Lh3w8l/0841+0FMMnXr4sP314L/AwZ2I9hXgulcCRbbHN+TpziXB7wnx9dQRMREREREREFCJ8PrSQiIiIiIqKhYSBHREREREQUYBjIERERERERBRgGckRERERERAGGgRwREREREVGAYSBHREREREQUYBjIERERERERBRgGckRERERERAHm/wHBjV7AXiBCTQAAAABJRU5ErkJggg==\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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.442103Z", "iopub.status.busy": "2021-02-02T06:51:24.440643Z", "iopub.status.idle": "2021-02-02T06:51:24.519114Z", "shell.execute_reply": "2021-02-02T06:51:24.520133Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.524405Z", "iopub.status.busy": "2021-02-02T06:51:24.523078Z", "iopub.status.idle": "2021-02-02T06:51:24.543215Z", "shell.execute_reply": "2021-02-02T06:51:24.544179Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.548608Z", "iopub.status.busy": "2021-02-02T06:51:24.547266Z", "iopub.status.idle": "2021-02-02T06:51:24.625742Z", "shell.execute_reply": "2021-02-02T06:51:24.626728Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.631241Z", "iopub.status.busy": "2021-02-02T06:51:24.629837Z", "iopub.status.idle": "2021-02-02T06:51:24.648995Z", "shell.execute_reply": "2021-02-02T06:51:24.649953Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:24.654257Z", "iopub.status.busy": "2021-02-02T06:51:24.652927Z", "iopub.status.idle": "2021-02-02T06:51:25.678150Z", "shell.execute_reply": "2021-02-02T06:51:25.679242Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 1999Q2 1999Q3 1999Q4 2000Q1 2000Q2\n", "1999Q3 2.552620 NaN NaN NaN NaN\n", "1999Q4 3.010790 3.588286 NaN NaN NaN\n", "2000Q1 3.342616 3.760863 3.226165 NaN NaN\n", "2000Q2 NaN 3.885850 3.498599 3.885225 NaN\n", "2000Q3 NaN NaN 3.695908 3.975918 4.196649\n" ] } ], "source": [ "# Setup forecasts\n", "nforecasts = 3\n", "forecasts = {}\n", "\n", "# Get the number of initial training observations\n", "nobs = len(endog)\n", "n_init_training = int(nobs * 0.8)\n", "\n", "# Create model for initial training sample, fit parameters\n", "init_training_endog = endog.iloc[:n_init_training]\n", "mod = sm.tsa.SARIMAX(training_endog, order=(1, 0, 0), trend='c')\n", "res = mod.fit()\n", "\n", "# Save initial forecast\n", "forecasts[training_endog.index[-1]] = res.forecast(steps=nforecasts)\n", "\n", "# Step through the rest of the sample\n", "for t in range(n_init_training, nobs):\n", " # Update the results by appending the next observation\n", " updated_endog = endog.iloc[t:t+1]\n", " res = res.append(updated_endog, refit=False)\n", " \n", " # Save the new set of forecasts\n", " forecasts[updated_endog.index[0]] = res.forecast(steps=nforecasts)\n", "\n", "# Combine all forecasts into a dataframe\n", "forecasts = pd.concat(forecasts, axis=1)\n", "\n", "print(forecasts.iloc[:5, :5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have a set of three forecasts made at each point in time from 1999Q2 through 2009Q3. We can construct the forecast errors by subtracting each forecast from the actual value of `endog` at that point." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "execution": { "iopub.execute_input": "2021-02-02T06:51:25.684199Z", "iopub.status.busy": "2021-02-02T06:51:25.682739Z", "iopub.status.idle": "2021-02-02T06:51:25.742628Z", "shell.execute_reply": "2021-02-02T06:51:25.742231Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 1999Q2 1999Q3 1999Q4 2000Q1 2000Q2\n", "1999Q3 0.797380 NaN NaN NaN NaN\n", "1999Q4 -0.160790 -0.738286 NaN NaN NaN\n", "2000Q1 0.417384 -0.000863 0.533835 NaN NaN\n", "2000Q2 NaN 0.304150 0.691401 0.304775 NaN\n", "2000Q3 NaN NaN -0.925908 -1.205918 -1.426649\n" ] } ], "source": [ "# Construct the forecast errors\n", "forecast_errors = forecasts.apply(lambda column: endog - column).reindex(forecasts.index)\n", "\n", "print(forecast_errors.iloc[:5, :5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To evaluate our forecasts, we often want to look at a summary value like the root mean square error. Here we can compute that for each horizon by first flattening the forecast errors so that they are indexed by horizon and then computing the root mean square error fore each horizon." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "execution": { "iopub.execute_input": "2021-02-02T06:51:25.746134Z", "iopub.status.busy": "2021-02-02T06:51:25.745427Z", "iopub.status.idle": "2021-02-02T06:51:25.788014Z", "shell.execute_reply": "2021-02-02T06:51:25.788621Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 1999Q2 1999Q3 1999Q4 2000Q1 2000Q2\n", "horizon \n", "1 0.797380 -0.738286 0.533835 0.304775 -1.426649\n", "2 -0.160790 -0.000863 0.691401 -1.205918 -0.311464\n", "3 0.417384 0.304150 -0.925908 -0.151602 -2.384952\n" ] } ], "source": [ "# Reindex the forecasts by horizon rather than by date\n", "def flatten(column):\n", " return column.dropna().reset_index(drop=True)\n", "\n", "flattened = forecast_errors.apply(flatten)\n", "flattened.index = (flattened.index + 1).rename('horizon')\n", "\n", "print(flattened.iloc[:3, :5])" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "execution": { "iopub.execute_input": "2021-02-02T06:51:25.793021Z", "iopub.status.busy": "2021-02-02T06:51:25.791374Z", "iopub.status.idle": "2021-02-02T06:51:25.798479Z", "shell.execute_reply": "2021-02-02T06:51:25.798825Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "horizon\n", "1 3.292700\n", "2 3.421808\n", "3 3.280012\n", "dtype: float64\n" ] } ], "source": [ "# Compute the root mean square error\n", "rmse = (flattened**2).mean(axis=1)**0.5\n", "\n", "print(rmse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Using `extend`\n", "\n", "We can check that we get similar forecasts if we instead use the `extend` method, but that they are not exactly the same as when we use `append` with the `refit=True` argument. This is because `extend` does not re-estimate the parameters given the new observation." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "execution": { "iopub.execute_input": "2021-02-02T06:51:25.807552Z", "iopub.status.busy": "2021-02-02T06:51:25.805873Z", "iopub.status.idle": "2021-02-02T06:51:26.509498Z", "shell.execute_reply": "2021-02-02T06:51:26.510450Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.514512Z", "iopub.status.busy": "2021-02-02T06:51:26.513364Z", "iopub.status.idle": "2021-02-02T06:51:26.565980Z", "shell.execute_reply": "2021-02-02T06:51:26.566865Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.570426Z", "iopub.status.busy": "2021-02-02T06:51:26.569617Z", "iopub.status.idle": "2021-02-02T06:51:26.602415Z", "shell.execute_reply": "2021-02-02T06:51:26.603086Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.605976Z", "iopub.status.busy": "2021-02-02T06:51:26.605312Z", "iopub.status.idle": "2021-02-02T06:51:26.611328Z", "shell.execute_reply": "2021-02-02T06:51:26.611769Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.616585Z", "iopub.status.busy": "2021-02-02T06:51:26.615569Z", "iopub.status.idle": "2021-02-02T06:51:26.619943Z", "shell.execute_reply": "2021-02-02T06:51:26.620333Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.623089Z", "iopub.status.busy": "2021-02-02T06:51:26.622426Z", "iopub.status.idle": "2021-02-02T06:51:26.627755Z", "shell.execute_reply": "2021-02-02T06:51:26.628385Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.633534Z", "iopub.status.busy": "2021-02-02T06:51:26.633136Z", "iopub.status.idle": "2021-02-02T06:51:26.636427Z", "shell.execute_reply": "2021-02-02T06:51:26.636798Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.639518Z", "iopub.status.busy": "2021-02-02T06:51:26.638840Z", "iopub.status.idle": "2021-02-02T06:51:26.644126Z", "shell.execute_reply": "2021-02-02T06:51:26.644565Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.647278Z", "iopub.status.busy": "2021-02-02T06:51:26.646599Z", "iopub.status.idle": "2021-02-02T06:51:26.654036Z", "shell.execute_reply": "2021-02-02T06:51:26.654592Z" } }, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.659593Z", "iopub.status.busy": "2021-02-02T06:51:26.658958Z", "iopub.status.idle": "2021-02-02T06:51:26.700219Z", "shell.execute_reply": "2021-02-02T06:51:26.699818Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/base/tsa_model.py:583: 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:583: 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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.709461Z", "iopub.status.busy": "2021-02-02T06:51:26.708788Z", "iopub.status.idle": "2021-02-02T06:51:26.718493Z", "shell.execute_reply": "2021-02-02T06:51:26.719341Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/base/tsa_model.py:379: 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": 1, "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": { "execution": { "iopub.execute_input": "2021-02-02T06:51:26.723168Z", "iopub.status.busy": "2021-02-02T06:51:26.721978Z", "iopub.status.idle": "2021-02-02T06:51:26.730044Z", "shell.execute_reply": "2021-02-02T06:51:26.730882Z" }, "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.9" } }, "nbformat": 4, "nbformat_minor": 2 }