Skip to content

Instantly share code, notes, and snippets.

@ChadFulton
Last active March 27, 2019 15:40
Show Gist options
  • Save ChadFulton/fac61cc2b2cb07eddfa3e87babb6e9cc to your computer and use it in GitHub Desktop.
Save ChadFulton/fac61cc2b2cb07eddfa3e87babb6e9cc to your computer and use it in GitHub Desktop.
Different types of SARIMAX
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import pandas as pd\n",
"import statsmodels.api as sm\n",
"import matplotlib.pyplot as plt\n",
"\n",
"np.set_printoptions(suppress=True, precision=5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Regression with AR(1) errors:**\n",
"\n",
"This is what Statsmodels' SARIMAX implements.\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"y_t & = x_t \\beta + \\varepsilon_t \\\\\n",
"\\varepsilon_t & = \\phi \\varepsilon_{t-1} + \\eta_t, \\quad \\eta_t \\sim N(0, \\sigma^2)\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"**AR(1) with exogenous intercept:**\n",
"\n",
"This is what you want.\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"y_t & = \\phi y_{t-1} + x_t \\beta + \\eta_t, \\quad \\eta_t \\sim N(0, \\sigma^2)\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here I've simulated some data that matches what you gave in your example; it's clear that this is the AR(1) with exogenous intercept case."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6cAAAFgCAYAAABHbg9VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl83FW9//H3ZyaTTLoladKW7jtdgFJoactO2UVZRFQQVLwo13tBvdfrvYrL1Yv3qlfvdfvJRUBxRVAW7y0IItCyiLQlhbZQ2tKNrrRJ920mmeX8/phJmNYuaZvMme93Xs/HI4/Od5vvZzqZk/OZs5lzTgAAAAAA+BTxHQAAAAAAACSnAAAAAADvSE4BAAAAAN6RnAIAAAAAvCM5BQAAAAB4R3IKAAAAAPCO5BQAAAAA4B3JKQAAAADAO5JTAAAAAIB3JKcAABwjMxtgZg+bWbOZrTKzT+f3P25m/11w3m/N7N7844iZfdnMVptZk5n90sxqCs79SP7YFjP7ipm9ZWYXFv/VAQBQHCSnAAAcAzOLSHpU0gJJAyVdIOkfzOwSSX8j6cNmdr6ZXS/pNEmfyV96Y/5nuqQRknpI+lH+OcdL+h9J10vqL6km/9wAAISWOed8xwAAQGCZ2VRJDzrnhhTsu03S8c65j5nZ1ZJ+IKla0lXOuT/nz3lG0sPOuf/Jb4+R9Hr+vC9KGuecuy5/rJuk7ZIuc849XbxXBwBA8VT4DgAAgIAbKmmAmW0v2BeV9EL+8WPKtYgubUtM8wZIWl2wvVq5v8v98sfWth1wzu01sy1dEDsAACWDbr0AABybtZJWOedqC356Oucuyx//D0mLJfU3s+sKrtugXGLbZoiktKRNkt6WNKjtgJlVS6rvyhcBAIBvtJwCAHBs5kraaWafl/RDSa2SxinXPbda0scknSxpuKT/NbPnnXPrJd0v6fNm9oSkZknfkPRb51zazB6SNNvMzpDUKOnfJFmRXxcAAEVFyykAAMfAOZeRdLmkiZJWSdos6SfKTWT0S0m3OufW57v0/lTSz8zMJN0r6VeSns9fl5T0qfxzLso/fkC5VtRdkpoktRTvlQEAUFxMiAQAQIkzsx7KTYg02jm3ync8AAB0BVpOAQAoQWZ2uZl1M7Pukv5L0muS3vIbFQAAXYfkFACA0nSlcpMmbZA0WtK1ju5OAIAQo1svAAAAAMA7Wk4BAAAAAN6V3FIyDQ0NbtiwYb7DAAAAAAAco3nz5m12zvXpyLkll5wOGzZMjY2NvsMAAAAAABwjM1vd0XPp1gsAAAAA8O6wyamZ3WtmTWb2+kGOm5n90MyWm9lCMzu14NhHzWxZ/uejnRk4AAAAACA8OtJy+nNJlx7i+LuUm+J+tKSbJd0pSWbWW9JXJU2VNEXSV82s7liCBQAAAACE02HHnDrnnjezYYc45UpJv8yvvTbbzGrNrL+k8yQ95ZzbKklm9pRySe79xxo0UO7ufHaF7nx2ue8wUIYG1XXTjFvPVEWUUSEoDXta0rr8//1Zm3e3+A4FALwY3Lub/vDps32H0Sk6Y0KkgZLWFmyvy+872P6/YmY3K9fqqiFDhnRCSEC4zVu9VbFoRJefPMB3KCgjSzbu1OyVW7UrmVZd90rf4QCSpLd3JLVy8x6de3wfDW/o7jscACi63iH6m9wZyakdYJ87xP6/3unc3ZLulqTJkycf8BwA70imshpa301fu+IE36GgjDwwd41mr9yqZDrjOxSgXTKV+3380NQhuuSE4zxHAwA4Fp3RL2udpMEF24MkbTjEfgDHKJHKqLoy6jsMlJm237lEK8kpSkdbclodo0wEgKDrjOR0hqSP5GftnSZph3PubUlPSrrYzOryEyFdnN8H4BglUxnFK6iIobiq8r9zyVTWcyTAO9p+H+MkpwAQeIft1mtm9ys3uVGDma1TbgbemCQ5534s6XFJl0laLmmvpI/lj201s69Lejn/VLe3TY4E4NgkUhnFaTlFkbW3nKZoOUXpSNByCgCh0ZHZeq87zHEn6ZaDHLtX0r1HFxqAg2lJZamIoejafudaSE5RQtq79VYygzQABB0lORBAiVRG8RgfXxRX2+8cLacoJW2/j1UMdQCAwKN2CwRQojVDyymKru13juQUpeSdllPKRAAIOpJTIGCcc0qmM0z+gaJr+51jQiSUkrbklDIRAIKP5BQImJZ0Vs5REUPxxWk5RQlKtOZn662gSgMAQUdJDgRMS77Vim69KLa2MadMiIRSkkxnFIuaKqJUaQAg6CjJgYBJ0IUNnrS3nLaSnKJ0JFoZ5gAAYUFyCgRMgmUT4EksGlEsanTrRUlpSTNBHACEBbVbIGDaJ/9g2QR4EK+IMiESSgotpwAQHiSnQMC0d+tl2QR4EK+M0nKKkpJI0XIKAGFBcgoEDC2n8CkeizAhEkpKMpVtn6wLABBslOZAwLDgPHyqjtFyitKSSNGtFwDCguQUCJi2Nf3oxgYf4iSnKDEtJKcAEBokp0DAtHfrpRsbPIjHou2/g0ApYMwpAIQHtVsgYNqXkqEyBg9y3XqZrRelI5nKMswBAEKC5BQImLZWqyqSU3jAhEgoNbkxp1RnACAMKM2BgEnScgqPmBAJpSbJOqcAEBokp0DAJFNZRSOmWNR8h4IyxJhTlJpkmuQUAMKC5BQImEQqo3hFRGYkpyi+eCyqRCvJKUpDOpNVKuPoSQIAIdGh5NTMLjWzpWa23My+cIDj3zOz+fmfN81se8GxTMGxGZ0ZPFCOEqkMk3/Am+rKqJJMiIQSkUyztBYAhEnF4U4ws6ikOyRdJGmdpJfNbIZz7o22c5xz/1hw/qcknVLwFAnn3MTOCxkob8lURlUVVMTgR7wiqtZMVpmsUzRC6z38amvFZ0IkAAiHjpTmUyQtd86tdM61SnpA0pWHOP86Sfd3RnAA/lqSllN4VF2Z+7PBuFOUgnfWfaZMBIAw6EhyOlDS2oLtdfl9f8XMhkoaLmlmwe64mTWa2Wwzu+og192cP6exubm5g6ED5SmZytKFDd60JQEkpygFJKcAEC4dSU4P1G/LHeTcayU95JwrrLUMcc5NlvQhSd83s5F/9WTO3e2cm+ycm9ynT58OhASUr0Qra/rBn7YkgOVkUAoSLK0FAKHSkRruOkmDC7YHSdpwkHOv1X5dep1zG/L/rpT0rPYdjwrgCOUWnKciBj9oOUUpaZucizIRAMKhI8npy5JGm9lwM6tULgH9q1l3zWyMpDpJLxXsqzOzqvzjBklnSnpj/2sBdFyS5BQeVbcnp8zYC//aW04r6U0CAGFw2Nl6nXNpM7tV0pOSopLudc4tMrPbJTU659oS1eskPeCcK+zyO07SXWaWVS4R/lbhLL8AjlwylaELG7ypplsvSghjTgEgXA6bnEqSc+5xSY/vt+9f99v+2gGu+4ukk44hPgD7SaayjDmFN22/e3TrRSkgOQWAcKGGCwRMgpZTeNQ+IVIrySn8a/s9pEwEgHAgOQUCJpnKKM46p/CkfUKkNGNO4R8tpwAQLiSnQIBks04t6aziFVTE4Ed1/ouRJC2nKAGJ/MRctJwCQDiQnAIBkky3zUxJRQx+xCtyfzaYEAmloK3ltKqC6gwAhAGlORAg7Wv6URGDJ+0tpySnKAHJVEZVFRFFIuY7FABAJ6CGCwTIO2v60XIKP9q6lNNyilKQTGUoDwEgREhOgQBh8g/4FomYKisi7a34gE+JVIYx+AAQIiSnQIC0LZtAcgqfqmNRuvWiJCRSWVpOASBESE6BAGlLCJiZEj7FYxHWOUVJaBtzCgAIB0p0IEDaJ0QiOYVH1bFo+8zRgE+MOQWAcCE5BQIkQcspSkA8FqXlFCUhyZhTAAgVklMgQN6ZEImPLvyJx6JKppkQCf4laDkFgFChhgsESILZelECqmNRJWk5RQlItGboSQIAIUJyCgRIknVOUQLisQjrnKIkJFNZVdGTBABCgxIdCBDWOUUpqK5kKRmUhmSKllMACBOSUyBAEq352XpZOgEexSuitJyiJCRTGb6sA4AQoYYLBEgynVEsaqqI8tGFP/HKaPuyRoAvzrnchEgkpwAQGtRwgQBJtNJKAP+qY3TrhX+pjFPWMQYfAMKkQ8mpmV1qZkvNbLmZfeEAx280s2Yzm5//+XjBsY+a2bL8z0c7M3ig3LSkaSWAf/FYhOQU3rV1La9imAMAhEbF4U4ws6ikOyRdJGmdpJfNbIZz7o39Tv2tc+7W/a7tLemrkiZLcpLm5a/d1inRA2WGllOUgupYVOmsUyqTVYwu5vCE2csBIHw6UquYImm5c26lc65V0gOSruzg818i6Snn3NZ8QvqUpEuPLlQAjK9CKWj7goRJkeBT++zlFZSJABAWHUlOB0paW7C9Lr9vf+8zs4Vm9pCZDT6Sa83sZjNrNLPG5ubmDoYOlJ9kKqs4a/rBs7bklK698ClByykAhE5Harl2gH1uv+1HJQ1zzk2Q9LSkXxzBtXLO3e2cm+ycm9ynT58OhASUpwTLJqAEtCenrczYC3/aZozmCzsACI+OlOjrJA0u2B4kaUPhCc65Lc65lvzmPZImdfRaAB3XksrQSgDv2rqWJ9O0nMKfRGu+Wy9f2AFAaHQkOX1Z0mgzG25mlZKulTSj8AQz61+weYWkxfnHT0q62MzqzKxO0sX5fQCOQiKVYXwVvKuuzP3paEsOAB/avhxhHD4AhMdhZ+t1zqXN7FblksqopHudc4vM7HZJjc65GZI+bWZXSEpL2irpxvy1W83s68oluJJ0u3Nuaxe8DqAsJGg5RQlo+4KECZHgU5KWUwAIncMmp5LknHtc0uP77fvXgse3SbrtINfeK+neY4gRQB4TIqEUxCuZEAn+tU+IRHIKAKFBLRcIkCTrnKIEtLWckpzCp3cmRKJMBICwIDkFAiSZJjmFf9XtLafM1gt/aDkFgPAhOQUCIpXJKpVxVMTgXVvXcsacwqe2lvsqhjoAQGhQogMBkaSVACWi7XeQ2XrhUzKVkZlUVUFVBgDCghIdCAgWnEepiLPOKUpAojWj6lhUZuY7FABAJ6GWCwREW8spY07hW1VFRGbvLOUB+MAYfAAIH5JTICBITlEqzEzxiqiSaSZEgj+J1izDHAAgZEhOgYBgZkqUkngswphTeJVMZ5gMCQBChlIdCIi2Madty3gAPlXHoqxzCq+S+TGnAIDwIDkFAiLR3q2Xjy38i1dGWUoGXiXTJKcAEDbUcoGAaOtCyZhTlIJ4BS2n8CvRyoRIABA2JKdAQLSkSU5ROqoro+1dzQEfEqks5SEAhAzJKRAQbS2ndGNDKYjHInTrhVctqQzDHAAgZCjVgYBIMlsvSggTIsG3RIoxpwAQNiSnQEAk8l0o6caGUlAVY0Ik+JVMMeYUAMKG5BQIiLZEoKqCjy38q45FlWSdU3iUSGVYWgsAQoZaLhAQLamMqioiikTMdyhALjlNMyES/HDOKcmESAAQOiSnQEDQSoBSEo9F2ifpAoqtJd02zIFqDACESYdKdTO71MyWmtlyM/vCAY5/1szeMLOFZvaMmQ0tOJYxs/n5nxmdGTxQTpKpjOIVJKcoDbmW04ycc75DQRli9nIACKeKw51gZlFJd0i6SNI6SS+b2Qzn3BsFp70qabJzbq+Z/Z2kb0v6YP5Ywjk3sZPjBspOIpWl5RQloyoWlXO5Fiy6VqLYkqz7DACh1JGW0ymSljvnVjrnWiU9IOnKwhOcc7Occ3vzm7MlDercMAEwMyVKSVuLVUuKcacoPlpOASCcOpKcDpS0tmB7XX7fwdwk6YmC7biZNZrZbDO76kAXmNnN+XMam5ubOxASUH6SLDiPEtL2RQnLycCHZIoxpwAQRoft1ivpQFODHnCQkZndIGmypHMLdg9xzm0wsxGSZprZa865Ffs8mXN3S7pbkiZPnswAJuAAEq0sOI/SUV2ZSwpITuFD2+8dvUkAIFw68pXjOkmDC7YHSdqw/0lmdqGkL0m6wjnX0rbfObch/+9KSc9KOuUY4gXKVjJNt16UjrYvSpIkp/CgJUW3XgAIo44kpy9LGm1mw82sUtK1kvaZddfMTpF0l3KJaVPB/jozq8o/bpB0pqTCiZQAdBAtpyglVXTrhUe0nAJAOB22W69zLm1mt0p6UlJU0r3OuUVmdrukRufcDEnfkdRD0oNmJklrnHNXSBon6S4zyyqXCH9rv1l+AXQQC86jlNByCp/aklNmMAeAcOnImFM55x6X9Ph++/614PGFB7nuL5JOOpYAAeQwIRJKSZzkFB61T4jE2s8AECrUdIGASKTo1ovS0fa7mGhlKRkUX3u33kqqMQAQJpTqQAA451jnFCWFbr3wiQmRACCcSE6BAGjNZJV1jK9C6WjrYs6ESPAh0cqESAAQRiSnQAC0ja+qquAji9IQr6TlFP4kUhlVREyxKGUiAIQJpToQAElmpkSJaZuIhuQUPjB7OQCEE8kpEABtXdgYX4VSEYuaohGjWy+8SDAGHwBCieQUCIBkmvFVKC1mpnhFpL3LOVBMLSytBQChRMkOBAAtpyhF1ZVRWk7hBUtrAUA4kZwCAdA+IRItBSgh8ViUMafwIpnKMAYfAEKImi4QAEnW9EMJIjmFL4lUpn1SLgBAeJCcAgHAbL0oRdWxKGNO4UUilW1fzggAEB4kp0AAtI3ro6UApSQei7SPhwaKqSWVUZx1nwEgdCjZgQBI0HKKEhSPMSES/Egw5hQAQonkFAiAtq6TtJyilDDmFL4kGXMKAKFEcgoEQFsCEK/kI4vSUU1yCk8SrbScAkAYUdMFAiCZyihiUmWUjyxKBxMiwZdkOqs4s5cDQOhQ0wUCINGaUTwWlZn5DgVoF49FGHOKostknVrTWcVZ9xkAQoeSHQiARCrDGqcoOfFKJkRC8bHuMwCEV4eSUzO71MyWmtlyM/vCAY5Xmdlv88fnmNmwgmO35fcvNbNLOi90oHwkU3RhQ+mJV0TVms4qm3W+Q0EZaR+DT5kIAKFz2OTUzKKS7pD0LknjJV1nZuP3O+0mSducc6MkfU/Sf+avHS/pWkknSLpU0v/knw/AEUimMnRhQ8lpm5Ammab1FMWToOUUAEKrogPnTJG03Dm3UpLM7AFJV0p6o+CcKyV9Lf/4IUk/stzguCslPeCca5G0ysyW55/vpc4Jv/i+9cQSZqdE0S1cv1294jHfYQD7iFfkvjD5+mOLVVXBlycojp2JlCSpii/sACB0OpKcDpS0tmB7naSpBzvHOZc2sx2S6vP7Z+937cD9b2BmN0u6WZKGDBnS0di9eHTBBu1KpnyHgTJ0wdh+vkMA9nHCwBo19KjUHxZu8B0KykzfnlUac1xP32EAADpZR5LTA00Puv8Ao4Od05Fr5Zy7W9LdkjR58uSSHrz04hfO9x0CAJSE04b1VuOXL/IdBgAACImO9IlZJ2lwwfYgSft/Td5+jplVSKqRtLWD1wIAAAAAylxHktOXJY02s+FmVqncBEcz9jtnhqSP5h9fI2mmc87l91+bn813uKTRkuZ2TugAAAAAgLA4bLfe/BjSWyU9KSkq6V7n3CIzu11So3NuhqSfSvpVfsKjrcolsMqf9zvlJk9KS7rFOcdsQgAAAACAfViugbN0TJ482TU2NvoOAwAAAABwjMxsnnNucofOLbXk1Mx2SVp6BJfUSNrRReEcTIOkzUW8n4/XWOx7BuF+x/q+B+E1Bu2eR3q/oL2HPu4Zxte4//vO/2nw73ege3b13+Zy/D8t9fsdzXsetNdY6vcrxj19l+E+7hmG+xW+b0Odc306dJVzrqR+lOsqfCTn313qMXbC/Xy8xqLeMwj3O9b3PQivMWj3PNL7Be09DML/aRDuuf/7zv9p8O93oHt29d/mcvw/LfX7Hc17HrTXWOr3K8Y9fZfh5fA+dsX9jrZMDsMK1o/6DqAIfLzGYt8z7PfzcU9eY/Dv5+OevMbg38/HPXmNwb+fj3vyGoN/Px/35DUG/34HVYrdehtdB/sk+xKEGNH5eN+Dj/ewPPG+lwfe5/LDe14eeJ+D6Wjft1JsOb3bdwAdEIQY0fl434OP97A88b6XB97n8sN7Xh54n4PpqN63kms5BQAAAACUn1JsOQUAAAAAlBmS00Mws92+Y0DxmFnGzOYX/Aw7xLnnmdljxYsOHWFmzsx+VbBdYWbNvFfhZ2bvzb//Y33Hgs7F5xrUx8rL4d5vM3vWzBiDGlIkp8A7Es65iQU/b/kOCEdsj6QTzaw6v32RpPVH8gRmVtHpUaEYrpP0Z0nXHslFZhbtmnDQiY75cw0ACAaS08Mwsx5m9oyZvWJmr5nZlfn9w8xssZndY2aLzOxPBX84ERJmFjWz75jZy2a20Mz+tuBwLzP7vZm9YWY/NjM+T6XhCUnvzj++TtL9bQfMbIqZ/cXMXs3/Oya//0Yze9DMHpX0p+KHjGNhZj0knSnpJuWT03zvhucP9Bk1s91mdruZzZF0ur/IcQSO5nP9gplNLDjvRTObUNSo0Wn277FkZj8ysxvzj98ys38rqKvRgyLgDvV+I9yoTB9eUtJ7nXOnSpou6b/NzPLHRku6wzl3gqTtkt7nKUZ0juqCLr2/z++7SdIO59xpkk6T9AkzG54/NkXSP0k6SdJISVcXPWIcyAOSrjWzuKQJkuYUHFsi6Rzn3CmS/lXSNwqOnS7po86584sWKTrLVZL+6Jx7U9JWMzs1v/9gn9Hukl53zk11zv256NHiaBzN5/onkm6UJDM7XlKVc25h0SJGsW3O19XulPQ538EAODp0Xzs8k/QNMztHUlbSQEn98sdWOefm5x/PkzSs+OGhEyWccxP323expAlmdk1+u0a5LyVaJc11zq2UJDO7X9JZkh4qVrA4MOfcwvx44eskPb7f4RpJvzCz0ZKcpFjBsaecc1uLEiQ623WSvp9//EB++w86+Gc0I+lhD3HiKB3l5/pBSV8xs3+W9DeSfl6UYOHLI/l/54kvi4HAIjk9vOsl9ZE0yTmXMrO3JMXzx1oKzstIoltv+JikTznnntxnp9l5ylWCCrEuU+mYIem/JJ0nqb5g/9clzXLOvTdf0X224NieIsWGTmRm9ZLOV25MopMUVe6z+LgO/hlNOucyxYsSneSIPtfOub1m9pSkKyV9QBITqARbWvv2+Ivvd7ytTpYR9dswONz7jZCiW+/h1Uhqyiem0yUN9R0QiupJSX9nZjEp1zXMzLrnj00xs+H5cWwfVG4yFpSGeyXd7px7bb/9NXpnIpUbixoRuso1kn7pnBvqnBvmnBssaZVyraR8RsPlaD7XP5H0Q0kv0zMi8FZLGm9mVWZWI+kC3wGhS/F+lymS04PIz9jZIuk+SZPNrFG5VtQlXgNDsf1E0huSXjGz1yXdpXe+kX1J0rckva5cZfj3B3wGFJ1zbp1z7gcHOPRtSd80sxeVa2FD8F2nv/7sPSzpQ+IzGipH87l2zs2TtFPSz4oQIrpAW33MObdW0u8kLVSubvaq18DQJXi/Yc7RE/FAzOxkSfc456b4jgUAcGTyXe8/55x7j+9Y4I+ZDVCum+9Y51zWczg4CtTHygvvN2g5PQAz+6Ry09R/2XcsAADgyJnZR5Sb1fdLJKbBRH2svPB+Q6LlFAAAAABQAmg5BQAAAAB4R3IqycwGm9ksM1tsZovM7DP5/b3N7CkzW5b/ty6/f6yZvWRmLWb2uf2e6x/zz/G6md2fXzAcAAAAAHAIJKc5aUn/5JwbJ2mapFvMbLykL0h6xjk3WtIz+W1J2irp08qtt9bOzAbm9092zp2o3KyB1xbnJQAAAABAcJGcSnLOve2ceyX/eJekxZIGKrdw9y/yp/1C0lX5c5qccy9LSh3g6SokVeenwu4maUMXhw8AAAAAgUdyuh8zGybpFOVm+OvnnHtbyiWwkvoe6lrn3HrlWlPXSHpb0g7n3J+6Ml4AAAAACAOS0wJm1kO5xdv/wTm38yiur1OutXW4pAGSupvZDZ0bJQAAAACED8lpnpnFlEtM73POPZLfvcnM+ueP95fUdJinuVDSKudcs3MuJekRSWd0VcwAAAAAEBYkp5LMzCT9VNJi59x3Cw7NkPTR/OOPSvq/wzzVGknTzKxb/jkvUG78KgAAAADgEMw55zsG78zsLEkvSHpNUja/+4vKjTv9naQhyiWe73fObTWz4yQ1SuqVP3+3pPHOuZ1m9m+SPqjcDMCvSvq4c66lmK8HAAAAAIKG5BQAAAAA4B3degEAAAAA3pGcAgAAAAC8IzkFAAAAAHhHcgoAAAAA8I7kFAAAAADgHckpAAAAAMA7klMAAAAAgHckpwAAAAAA70hOAQAAAADekZwCAAAAALwjOQUAAAAAeEdyCgAAAADwjuQUAABPzOznZvbvvuMAAKAUkJwCAAAAALwjOQUAAAAAeEdyCgDAETKzAWb2sJk1m9kqM/t0fv/XzOx3ZvZLM9tlZovMbHLBdaeY2Sv5Y7+VFN/veT9hZsvNbKuZzTCzAQXHLjazpWa2w8z+x8yeM7OPF+1FAwDQxUhOAQA4AmYWkfSopAWSBkq6QNI/mNkl+VOukPSApFpJMyT9KH9dpaT/lfQrSb0lPSjpfQXPe76kb0r6gKT+klbnn0dm1iDpIUm3SaqXtFTSGV34MgEAKDqSUwAAjsxpkvo45253zrU651ZKukfStfnjf3bOPe6cyyiXiJ6c3z9NUkzS951zKefcQ5JeLnje6yXd65x7xTnXolwierqZDZN0maRFzrlHnHNpST+UtLFrXyYAAMVV4TsAAAACZqikAWa2vWBfVNILyrV2FiaNeyXFzaxC0gBJ651zruD46oLHAyS90rbhnNttZluUa50dIGltwTFnZus66fUAAFASaDkFAODIrJW0yjlXW/DT0zl32WGue1vSQDOzgn1DCh5vUC7xlSSZWXfluvCuz187qOCYFW4DABAGJKcAAByZuZJ2mtnnzazazKJmdqKZnXaY616SlJb0aTOrMLOrJU0pOP4bSR8zs4lmViXpG5LmOOfekvQHSSeZ2VX5VthbJB3X2S8MAACfSE4BADgC+bGkl0uaKGmVpM2SfiKp5jDXtUq6WtKNkrZJ+qCkRwqOPyPpK5IeVq6ldKTy41idc5slvV/StyVwPqkKAAAgAElEQVRtkTReUqOklk57YQAAeGb7Dn0BAAClLj9j8DpJ1zvnZvmOBwCAzkDLKQAAAWBml5hZbb7L7xclmaTZnsMCAKDTkJwCABAMp0taoVw34sslXeWcS/gNCQCAzkO3XgAAAACAd7ScAgAAAAC8q/AdwP4aGhrcsGHDfIcBAAAAADhG8+bN2+yc69ORc0suOR02bJgaGxt9hwEAAAAAOEZmtrqj59KtFwAAAADgHckpAAAAAMC7kuvWCwDoOs45ZbJOGeeUzUqZ/HY265TOOmXbjucfZ52UdU7OOTkn1XSLqW/PuO+XAQDHzOXLuExh2edy5eH+5WS2oFx0knKLXeTKxbbtrGvbdmpbDKNw2xXcc/9rXcHztp9/gGsPusbGIRbfOMRVOtSiHYc8dtBrDnGvgz/dIe91sCuPJr7D3yuYulVGNX1sX99hdAqSUwDoIs45taSzSrRmtDeVUaI197O3Nd2+vbc1o0RrOvdvKqOWdFapdFapTFatGZf7N7/dtq81nVFqv2OtmaxSaad0NluQXGrfypZzx/xHOR6LaM4XL1RNdaxz/pMAlI1M1mlvazpXFqZy5d/egnJx331pJVO58q2loAxMZZxa0/kyr+0n7dSSeafsbD8vU1Ae5svAwkQ0G8IkBeVpaH03klMAKCeJ1ow2727R5t0t2rK7NffvnlbtTKS0M5nSjkRKOxPpgscp7UymlTnC2k8saopFI6qsiOT+jUYUi1r7dtux6lhUveIV+50XUUXUFI2YIpb7953HUtRMkYi982/hY1Pu3H32mSIRKWImM9ParXv1nSeX6uVVW3Xh+H5d9D8NoNRlsk7b97Zq8+5WbdndouZ8ubh9b6t2JtMFZeC+ZePe1swR3Sdi2qeMe6ccfKecbCv7aipjqtyvrIxF9y/n7IDlXDSifY/vd140Xw62HTPLxWdmMklmkskUMeWPWX7ffudY4fa+17adL8uXuW3ntD/Pvs97IGYHO3Lwa9T+3Ae77hDPeagnLdK9ihl7KYtFwzNSk+QUQNnb3ZLW+m0Jrd++V+u2JbR+W0Lrtif09vZEe+Vrz0EqVVUVEfWqjqmmOqZe8Qr17l6p4Q3d1SseU6/qCnWrrFC3yqi6VUYVj0Xbt6vz+7rFKlSd366ORRWNlO5fzGQqox88s0yzV24hOQVCKpN1atqVzJWD2xJav/2df5t2JrV5d6u27mk5YKujmdSzqqKgTIxpeEP39sc94m3lX4W6xaIFZWFFQRn5zv7KaOSQCReA8CE5BVAWMlmn1Vv26M1Nu7Vs0y692bRbK5t3a/32hLbvTe1zbmU0ooF11epfE9cpQ2pV371KDT0r1ZD/N7ddpfrulYrHop5eUfHFY1GdOqRWs1dt8R0KgGO0I5HKlYWbduvNTbu0vGm31mzdqw3bE0rvl3nWd6/UwLpqDe7dTacMqVVDj1z5V9+jSg09qtTQo1INPapUUx1TpIS/YANQ+khOAYTO9r2tmr92uxZt2Nle+VrevFut6Wz7OYPqqjWqbw+dMqRWA2u7aWBdtQbVVWtQbbUaelRRwTqIaSPq9YNnlmlHIsW4UyAAMlmnZU27tHDtDi3ZuEvLmnbpzU27tGlnS/s53SqjGt23hyYOrtW7J/TXwNp8eVhXrQG11epWSXURQHFQ2gAINOecVjTv0eyVWzRv9TbNX7tdqzbvaT/evyau4/v11Jmj6jW6X0+N6ddTo/r2UPcqir+jMW1Evb7/9DLGnQIlak9LWo2rt2nOyi16dc12LVy3vX1YQjwW0ai+PXTmyIZceXhcD43u21MDa6v5Qg5ASaB2BiBwmne1aNbSJv152WbNXrlFTbtyLQB9elZp4uBaXTNpkE4ZXKsTBtbQutfJJg6uVWVFhHGnQInIZJ1eXbNNs5Y26aUVW7Rw3Q6ls04VEdP4Ab109amDNHFwrSYOqdWw+u4lPa4dAEhOAZQ855yWbNylp97YpGeWNGnB2u2Scsno6SPqdfrIek0bUa9h9d2YPKOLMe4U8G9PS1ozlzRp5pImPbu0Sdv2phSNmE4eVKObzxmh00fWa9LQOrrjAgicopRaZlYr6SeSTlRuXdy/cc69VIx7AwiuNVv2asaC9fq/+Ru0rGm3zKSTB9Xqny46XheM66dx/XuSjHrAuFOg+FrSGT3/5mb93/z1enrxJiVTWfXuXqnpY/vqgrH9dPbxDeoV5/MIINiK9ZXaDyT90Tl3jZlVSupWpPsCCJhkKqPHX3tb981Zo3mrt0mSThtWp69feYIuPbG/+vSs8hwhGHcKFM/ypl369ew1+v2r67UjkVLv7pW6ZtIgXXHyQE0aWkc3XQCh0uXJqZn1knSOpBslyTnXKqm1q+8LIFjWbNmr++as1u8a12rb3pRGNHTX5y8dqysmDtDA2mrf4aEA406BrpXOZPXkok361ey3NHvlVsWipktP7K+rTx2os0Y1KBaN+A4RALpEMVpOR0hqlvQzMztZ0jxJn3HOtU+naWY3S7pZkoYMGVKEkACUimWbdumOWcs1Y8EGmZkuHt9PN0wbqjNG1tNlt0Qx7hToGi3pjB6et153Prdca7cmNKiuWv9y6Rh9YPJgNfSg1wiA8CtGcloh6VRJn3LOzTGzH0j6gqSvtJ3gnLtb0t2SNHnyZHfAZwEQKos27NCPZi7XHxdtVHUsqk+cPUIfO3O4jquJ+w4NHTBtRL1+yLhToFMkUxndP3eN7npupTbuTOrkQTX68rvH68Jx/ei2C6CsFCM5XSdpnXNuTn77IeWSUwBlaNPOpL79x6V6+JV16hmv0Kemj9LHzhyuuu6VvkPDEWgbd9r41lZdMI6uvcDRcM7psYVv61tPLNH67QlNGd5b33n/BJ01qoGeIwDKUpcnp865jWa21szGOOeWSrpA0htdfV8ApSWZyuie51fqzudWKJ1x+ttzR+jvzxtFq1tAFY47JTkFjtyCtdv19cfeUOPqbRrXv5e+c80EnTGqwXdYAOBVsWbr/ZSk+/Iz9a6U9LEi3RdACfjzss36/MMLtX57Qu868Tjd9q5xGlLPpN1B1j7udOVW36EAgbK7Ja1vPL5Yv5mzRg09qvSf7ztJ10waTPddAFCRklPn3HxJk4txLwClo7ASNqJPdz1w8zRNG1HvOyx0kvZxp3tTqulGCzhwOC8u36x/eWihNuxI6ONnDddnLhytnqxNCgDtitVyCqDM/GX5Zv1zvhJ28zkj9NmLjlc8FvUdFjrR2aMb9P2nl+n5Zc26/OQBvsMBStaelrS++cRi/Xr2Go1o6K6HPnm6Jg3t7TssACg5JKcAOlU26/TDmcv0/aeXaTiVsFCbOLhOdd1imrWkieQUOIjlTbv1yV/P04rm3fr4WcP1uUvG8EUdABwEySmATrMjkdJnfztfzyxp0tWnDtR/XHWSqiuphIVVNGI6b0xfPftmszJZx5g5YD9/fH2jPvfgAlVVRHTfTVOZ8AgADoPkFECnWLJxp/72V/O0fltCX7/yBN0wbShLIZSB6WP76vevrtf8tds1aWid73CAkpDJOv3Xn5bqzmdX6OTBtbrz+lM1oLbad1gAUPJITgEcs2eXNunvfv2KesYr9Nu/nUY33jJy7ug+ikZMs5Y0kZwCyi2bdct9r+iZJU26bsoQfe2K8aqqoAcJAHRExHcAAILt0QUb9IlfNmp4Q3c99qmzSEzLTE23mCYNrdPMJU2+QwG825lM6SM/nauZS5v09StP0DevPonEFACOAMkpgKN235zV+vQDr+qUwXV64G+nqW+vuO+Q4MH5Y/vqjbd36u0dCd+hAN4072rRtXfN1itrtumH156iD58+zHdIABA4JKcAjphzTnfMWq4v/f51nT+mr3550xT1Yq2+snX+2L6SpFlLmj1HAvixbttefeCul7Ry82795KOTmb0aAI4SySmAI/bj51bqO08u1VUTB+jHH57EsghlbnTfHhpYW03XXpSlpp1JXXfPbG3Z3aJf3zRV543p6zskAAgsklMAR+S3L6/Rf/5xia44eYC++4GJikUpRsqdmemCcX314vLNSqYyvsMBimZHIqWP3DtXW3a36pc3TdXkYYy5B4BjQa0SQIc9uWijbnvkNZ1zfB/91/tPVoR1LZE3fWxfJVIZzV65xXcoQFEkWjP6+C9e1orm3brrw5M0cXCt75AAIPBITgF0yOyVW/Sp+1/VhEG1+vENp6qyguID7zh9RL3isYhm0bUXZSCVyerW37yixtXb9L0PTtTZo/v4DgkAQoHaJYDDenPTLn3iF40a0rubfnbjaepWyRLJ2Fc8FtVZoxr0zJImOed8hwN0GeecvvK/r+uZJU26/coT9Z4JTH4EAJ2F5BTAIe1KpvTJX81TvDKqX900RXXdK32HhBI1fWxfrduW0NJNu3yHAnSZ++eu1QMvr9Ut00fqw9OG+g4HAEKF5BTAQTnn9LkHF2j11r2640Onqn9Nte+QUMIuHn+cIiY9tuBt36EAXWLB2u362oxFOuf4PvrsRWN8hwMAoUNyCuCg7np+pZ5ctEm3vWuspgxnFkocWp+eVTpzVINmLNhA116EztY9rfr7+15Rn55V+sEHJyrKhHAA0OlITgEc0F9WbNa3/7hE757QXzedNdx3OAiIyycM0Jqte7Vw3Q7foQCdJpN1+swDr6p5d4t+fMMkhjcAQBchOQXwVzbtTOrT97+q4Q3d9Z/vmyAzWgjQMZeceJxiUdOMBRt8hwJ0mh88/aZeWLZZX7/yBJ00qMZ3OAAQWiSnAPbhnNNtj7ym3S1p3fXhSepRxcy86Lia6pjOPb6vHlu4QdksXXsRfPPXbtePZi3XNZMG6YOnDfEdDgCEGskpgH088sp6zVzSpH+5ZKxG9e3pOxwE0BUTB2jTzhbNfWur71CAY9KSzuifH1yg43rF9dXLx/sOBwBCj+QUQLumnUn926OLNHlonW48Y5jvcBBQF47rq+pYVI/StRcB98NnlmlZ0259830T1DMe8x0OAIQeySkASbnuvF/8/etqSWf17WsmKMJMlDhK3SordOH4fnr8tbeVymR9hwMclYXrtuvHz63UByYP0rnH9/EdDgCUBZJTAJKkGQs26OnFm/S5i8doRJ8evsNBwF0+ob+27U3pxeWbfYcCHLFcd96FauhRqS+9m+68AFAsJKcA1LyrRV+dsUinDKnV37BsDDrBuWP6qGe8gll7EUg/mrlcSzft0jevPkk11XTnBYBiKVpyamZRM3vVzB4r1j0BdMx3nlyiPS1pfeeaCSwsj05RVRHVpSccpz8t2qRkKuM7HKDDVm/Zo7ueW6n3njJQ54/t5zscACgrxWw5/YykxUW8H4AOWLRhhx6ct04fO3M4s/OiU105caB2t6T15KKNvkMBOuxbTyxRRdT0hXeN9R0KAJSdoiSnZjZI0rsl/aQY9wPQMc45/ccfFqu2OqZbpo/yHQ5C5oyR9RrSu5vum7PGdyhAh8xdtVVPvL5Rnzx3pPr1ivsOBwDKTrFaTr8v6V8kHXDaRjO72cwazayxubm5SCEBeGZxk/6yYov+8aLjGVeFTheJmD40dYjmrtqqZZt2+Q4HOKRs1unf//CGjusV1yfOHuE7HAAoS12enJrZeyQ1OefmHewc59zdzrnJzrnJffowXTtQDKlMVt94fLFG9umu66YM8R0OQur9kwYpFjVaT1Hy/m/Bei1ct0P/cukYVVdGfYcDAGWpGC2nZ0q6wszekvSApPPN7NdFuC+AQ7hv9mqt3LxHX3r3OMWiTNyNrlHfo0rvOrG/Hn5lnRKtTIyE0pRozejbf1yqCYNqdNXEgb7DAYCy1eU1Uufcbc65Qc65YZKulTTTOXdDV98XwMHt2JvS959ZprNGNWj6mL6+w0HIXT91iHYl03p0IcvKoDTd88JKvb0jqS+/e7wizFgOAN7QXAKUobueX6EdiZS+eNk4mVERQ9eaMry3RvftQddelKTte1t113MrdOkJx2nK8N6+wwGAslbU5NQ596xz7j3FvCeAfe1IpPSrl1brspP6a/yAXr7DQRkwM10/dYgWrN2u19fv8B0OsI+f/+Ut7WnN6B8uGu07FAAoe7ScAmXmVy+9pV0taf39eSN9h4Iy8t5TBykei+i+Oat9hwK0292S1s9efEsXjuunscfxZR0A+EZyCpSRva1p3fviW5o+po9OGFDjOxyUkZrqmK44eYD+99UN2plM+Q4HkCT9Zs5q7UikdMt0vqwDgFJAcgqUkQfmrtXWPa269fxRvkNBGbph2lAlUhk9MJexp/AvmcronhdW6cxR9TplSJ3vcAAAIjkFykZLOqO7n1+pqcN7a9JQJv1A8U0YVKszR9XrnhdWKZliWRn49eC8dWre1aJbzuPLOgAoFSSnQJn4/SvrtXFnUrdMpyIGf26ZPkrNu1r04Lx1vkNBGUtlsvrxsyt0ypBanT6y3nc4AIA8klOgDKQzWd353ApNGFSjs0c3+A4HZez0EfU6ZUitfvzsCqUyWd/hoEzNmL9B67cndOv0USynBQAlhOQUKAOPv75Rq7fs1d+fR0UMfpmZbp0+Suu3JzRj/gbf4aAMOed053MrNPa4njp/bF/f4QAACpCcAmXgl395S8Pqu+ni8f18hwLo/LF9Nfa4nvqfZ5crm3W+w0GZ+cuKLVretFs3nzOCL+sAoMSQnAIht2TjTjWu3qbrpw5VJEJFDP6ZmW6ZPkormvfoj4s2+g4HZea+OatV2y2my07q7zsUAMB+SE6BkPvNnDWqrIjofZMG+Q4FaHfZSf01vKG77pi1XM7ReoriaNqZ1J8WbdL7Jw1SPBb1HQ4AYD8kp0CI7WlJ65FX1uvdJ/VX7+6VvsMB2kUjpr87d6QWbdipZxY3+Q4HZeJ3jWuVzjpdN2WI71AAAAdAcgqE2KMLNmh3S1rXT6UihtLz3lMHakRDd33zicXM3Isul8k63T93rc4cVa8RfXr4DgcAcAAkp0CI3Tdnjcb066lJQ+t8hwL8lVg0otsuG6cVzXt0/9w1vsNByD33ZpPWb0/o+qlDfYcCADgIklMgpBau267X1u/Q9dOGMCMlStaF4/rq9BH1+t5Tb2pHIuU7HITYr2evUZ+eVbqIWcsBoGSRnAIhdd/sNepWGdV7TxnoOxTgoMxMX3r3OG1PpHTHrOW+w0FIrdu2V7OWNuna0wYrFqXqAwClihIaCKEdiZRmLNigKycOUM94zHc4wCGdOLBG15w6SD9/8S2t3rLHdzgIoQfmrpVJupaJkACgpJGcAiH0v6+uVyKV0YemMLYKwfC5S8YoGjF964klvkNByKQzWf22ca2mj+mrgbXVvsMBABwCySkQQr9/db3G9e+lkwbV+A4F6JB+veL65Lkj9cTrGzV75Rbf4SBE/rJii5p3tej9kwf7DgUAcBgkp0DIrNmyV/PXbteVEwf4DgU4IjefM0KD6qp12yOvKdGa8R0OQmLGgg3qWVWh88b08R0KAOAwSE6BkHl04QZJ0nsm9PccCXBkqiuj+vb7JmjV5j367z8t9R0OQiCZyujJ1zfqkhOPUzwW9R0OAOAwSE6BkHl0wQZNGlqnQXXdfIcCHLEzRjXo+qlD9NMXV2ne6q2+w0HAPfdms3a1pHX5yfQkAYAgIDkFQuTNTbu0ZOMuXU6rKQLstsvGaUBNtf75oYVKpujei6M3Y8EG1Xev1Jkj632HAgDoAJJTIEQeW7BBEZMuIzlFgPWoqtC33neSVjbv0feeetN3OAioPS1pPbN4ky47qb8qWNsUAAKhy0trMxtsZrPMbLGZLTKzz3T1PYFy5JzTjAUbdPrIevXtGfcdDnBMzh7dR9dNGax7Xlipeau3+Q4HAfT04k1KprJ06QWAACnGV4lpSf/knBsnaZqkW8xsfBHuC5SV19fv1Ftb9uoKKmIIiS9eNk4Daqv1qd+8oi27W3yHg4CZMX+D+tfENXlone9QAAAd1OXJqXPubefcK/nHuyQtljSwq+8LlJsZC9YrFjVdegJdehEOPeMx3Xn9JG3e06rPPDBfmazzHRICYvveVj2/rFnvmdBfkYj5DgcA0EFFHYRhZsMknSJpzn77bzazRjNrbG5uLmZIQChks06PLXxb5x7fRzXdYr7DATrNSYNq9O9Xnqg/L9+s7z7F8jLomD++vlGpjNMVJ/NdOAAESdGSUzPrIelhSf/gnNtZeMw5d7dzbrJzbnKfPiySDRypxtXb9PaOJGOrEEofOG2wrpsyWHfMWqE/LdroOxwEwKMLN2hYfTedOLCX71AAAEegKMmpmcWUS0zvc849Uox7AuXkDws3KB6L6MJx/XyHAnSJr15+gk4aWKN/+t0Crdq8x3c4KGFbdrfopRVbdPnJA2RGl14ACJJizNZrkn4qabFz7rtdfT+g3Djn9PTiJp01qo+6V1X4DgfoEvFYVHfecKoqoqaP3jtXTTuTvkNCiXp2abOyTrp4/HG+QwEAHKFitJyeKenDks43s/n5n8uKcF+gLCxr2q312xO6YFxf36EAXWpQXTf97GNTtHl3iz5y71ztSKR8h4QSNHNJk/r2rNIJA+jSCwBBU4zZev/snDPn3ATn3MT8z+NdfV+gXDyzuEmSNH0MySnCb+LgWt314Ula0bxbn/hFo5KpjO+QUEJSmayef7NZ08f0ZZZeAAigos7WC6DzzVrSpPH9e+m4mrjvUICiOHt0H333AxP18uqtuvU3ryidyfoOCSWi8a1t2tWS1vSxfFkHAEFEcgoE2I69Kc1bs03nUxFDmbn85AG6/YoT9PTiJn3mt/PVmiZBhTRzySbFoqazRjf4DgUAcBSYPQUIsOeWNSuTdbQSoCx9+PRhSqQy+sbjS7Q7mdadN5yqbpX8WStnM5c0adqIevVgcjgACCRaToEAm7WkSb27V2ri4FrfoQBe3HzOSH3r6pP0wrJmffinc7VjL5Mklas1W/ZqRfMext8DQICRnAIBlck6Pbu0Sece30dRJv5AGbt2yhDd8aFT9dq6Hfrg3S+xzEyZmrlkkyQxzAEAAozkFAio+Wu3adveFBUxQNK7Tuqve288TWu27tUVP3pRr6zZ5jskFNkzS5o0oqG7hjV09x0KAOAokZwCATVzSZOiEdM5x/fxHQpQEs4a3aAHP3m6YhWmD971ku6bs1rOOd9hoQj2tKQ1Z+VWvqwDgIAjOQUCauaSZk0aWqea6pjvUICSccKAGj1661k6Y2SDvvT71/X5hxeyFmoZeHH5ZrVmsiSnABBwJKdAAL29I6HFb++kIgYcQG23St1742n69Pmj9LvGdbrqjhf1+vodvsNCF5q5pEk9qio0eVhv36EAAI4BySkQQLOWNEti4g/gYKIR02cvHqOf3Xiatu5p1VV3vKjvPvUm66GGkHNOs5Y26ezRDaqsoFoDAEFGKQ4E0LNLmzSwtlqj+/bwHQpQ0qaP7as//eM5uvzkAfrhM8t0Ja2oobN00y5t2tnCEjIAEAIkp0DAZLNOc1Zt1Zmj6mXGEjLA4dR2q9T3PjhR93xksjbvbtHlP/qzvvDwQjXtYsmZMJi9Yosk6fSR9Z4jAQAcK5JTIGCWbNylHYmUpo2gIgYciYvG99PTnz1XN505XA+/sk7Tv/Os7pi1nAmTAm72yq0aVFetwb27+Q4FAHCMSE6BgJm9MtdKMJXkFDhiNdUxffk94/WnfzxXZ4xq0HeeXKrp//Wsfv7iKpLUAMr1JNnCl3UAEBIkp0DAzF65RUN6d9PA2mrfoQCBNbyhu+75yGT95hNTNbium7726Bs66z9n6sfPrdDulrTv8NBBbzbt0ra99CQBgLCo8B0AgI5rG2968fh+vkMBQuGMkQ06Y2SD5qzcoh/NWq5vPbFEd8xcrqtPHagbpg3V6H49fYeIQ2gbbzp1OEvIAEAYkJwCAcJ4U6BrTB1Rr6kj6jV/7Xb9/MVVun/uWv3ipdWaOry3PjR1iC4a30/dKvmTWWoYbwoA4cJfWiBA3hlvSisB0BUmDq7V9689RV95T4senLdO981Zrc88MF/dKqO6aHw/XXHyAJ09ug/raZaAtvGmF4yjJwkAhAXJKRAgc1Zt0eDe1RpURysB0JXqe1Tpk+eO1M1nj9Dct7ZqxoINevy1t/V/8zeoV7xC547pqwvG9tW5x/dRXfdK3+GWJcabAkD4kJwCAdE23vQiWgmAoolETNNG1GvaiHp97fIT9MKyZj25aKNmLmnWows2KGLSqUPqdPrIep0+ol6nDq1TPBb1HXZZYLwpAIQPySkQEEs37dJ2WgkAbyorIrpgXD9dMK6fslmnhet3aOaSJj33ZrPumLVc/2/mclVGIzp5cI1OGVKniYNrNXFwrfrXxGVmvsMPHcabAkD4kJwCAcF4U6B0RCLWnnx+9qLjtSuZUuNb2/TSyi2au2qrfv7iW2rNZCVJfXpWaexxPTW6b08d36+Hjj+u5/9v785j5KzrOI6/PzM7e/TYpRQKsvSKVmtFBKwVhBBOY6IJ4AVVIkQ84x0xMV6JR5B4RYxERTzQKAQVYj0BiUZAlHKFU1EuaSFAC7TddrvHzNc/nqe703Xb7S47zzPPzOeVTOaZ3/PMfL/P/PLMPN95fs8zrFg0j/ndlZzXorh8vqmZWWtycWpWEH9/aDOHLvD5pmbNaH53hRNXLuLElYsAGBqtct/jW7nzsee4e8MWHnhqGz+/5VF2jtTGnnNIX/dYobo4/e/iQxfMoX9BD/O6/PW8Nz7f1MysNfnbz6wAdp1veoqPEpgVQldHmSOXLODIJQvG2qq1YMOzO3jgyQEeeHJbehvgbw9uZni0ttvz+3oq9O/XQ/+CHvr362FRbxcHzO3igPmdLJzbxQHzu1g4t7Ntz2/1+aZmZq0pk+JU0uuAi4AycGlEXJhFXLNW4fNNzYqvXBJLF85l6cK5nLpq/IemWi3YNDDEhucG2fjsIBvT+w3P7uDRzdu5+cHNDAyNTvqa8y34Ge0AAAm/SURBVLs6WDivk745nfR2d9DXU6G3p0JvdyWd7qC3u8K8rg56OsvMSW89nR3MqZTp6SzT1VEq3DmxPt/UzKw1Nbw4lVQGLgZOBTYA6yWti4j7Gh3brFWMnW/qowRmLadUEot6u1nU281RdUda6w0OV9k0MMSmgSE2Dwwn99uHeXpbcr9lcIQtgyNsfHZwbHq0FvsWX9BTSQvWtHjtqpTpKpeodIhKuUSlXKKzXKKzo0SlPN7W1VEam650aGyZjlKJcglKEuVScts1Pd5WN1+iNGG5pI2xaUmUBJJ8vqmZWYvK4sjpGuA/EfEQgKQrgNOAQhanq7/0J7buHMk7DWszo9Ua/fv5KIFZu+rpLLN4/zn7/BkQEewcqY0VqtuHRxkcrrJjuMrgSJXB4VF27Ho81j7eNjRaY2S0xtBIjW07RxkerTFSrTFSDUaqNYZHawxXx9uq+1gIz6ZjPJLEzKzlZFGc9gOP1T3eALy6fgFJ7wHeA7BkyZIMUpq5s49estsFLcyycuyLvCNmZvtGEj2dybDdg/u6Gx6vWkuL1mpS1I5Ug2oEtVpSuI5NR/K4VmN8eqwtJrSx2/wgKbprEVTKJV676uCGr5eZmWUri+J0shNZdvuJNSIuAS4BWL16dfY/v07DR095cd4pmJmZNZVkmG65bS/QZGZms6OUQYwNwOK6x4cCj2cQ18zMzMzMzAoii+J0PbBC0nJJncBZwLoM4pqZmZmZmVlBNHxYb0SMSvogcA3JX8n8MCLubXRcMzMzMzMzKw5FNNcpnpK2Af+axlP6gC0NSmdPDgA2ZRgvj3XMOmYR4j3ffi/COhYt5nTjFa0P84jZius4sd/9nhY/3mQxG/3d3I7vabPHm0mfF20dmz1eFjHz/gzPI2YrxKvvt6URceA+PSsimuoG3DrN5S9p9hxnIV4e65hpzCLEe779XoR1LFrM6cYrWh8W4T0tQsyJ/e73tPjxJovZ6O/mdnxPmz3eTPq8aOvY7PGyiJn3Z3g79GMj4s30MzmLc04b7Td5J5CBPNYx65itHi+PmF7H4sfLI6bXsfjx8ojpdSx+vDxieh2LHy+PmF7H4sfbo2Yc1ntrRKzOO4+9KUKONvvc78XnPmxP7vf24H5uP+7z9uB+LqaZ9lszHjm9JO8E9kERcrTZ534vPvdhe3K/twf3c/txn7cH93Mxzajfmu7IqZmZmZmZmbWfZjxyamZmZmZmZm3GxamZmZmZmZnlzsXpXkgayDsHy46kqqQ7627L9rLsCZJ+m112ti8khaSf1j3ukPS0+6r1SToj7f+Veedis8vbtXl/rL1M1d+S/iLJF0hqUS5OzcYNRsQRdbdH8k7Ipm07cJiknvTxqcDG6byApI5Zz8qysBa4EThrOk+SVG5MOjaLnvd2bWZmxeDidAqS5km6XtLtku6WdFravkzS/ZK+L+leSdfWfXFai5BUlvRVSesl3SXpvXWzeyVdLek+Sd+V5O2pOfwBeH06vRa4fNcMSWsk/U3SHen9S9L2cyX9QtJvgGuzT9meD0nzgGOB80iL03R0w18n20YlDUj6gqR/AMfkl7lNw0y26xskHVG33E2SDs80a5s1E0csSfq2pHPT6Uckfb5uX80jKApub/1trc0701PbCZwREUcBJwJfl6R03grg4oh4GfAc8KaccrTZ0VM3pPfqtO08YEtEvAp4FfBuScvTeWuAjwMvB14IvDHzjG0yVwBnSeoGDgf+UTfvn8DxEXEk8Dnggrp5xwDnRMRJmWVqs+V04I8R8QDwjKSj0vY9baNzgXsi4tURcWPm2dpMzGS7vhQ4F0DSi4GuiLgrs4wta5vSfbXvAOfnnYyZzYyHr01NwAWSjgdqQD9wUDrv4Yi4M52+DViWfXo2iwYj4ogJba8FDpf05vRxH8mPEsPALRHxEICky4HjgF9mlaxNLiLuSs8XXgv8fsLsPuAySSuAACp1866LiGcySdJm21rgm+n0Fenj37HnbbQK/CqHPG2GZrhd/wL4rKRPAO8EfpxJspaXq9L72/CPxWaF5eJ0am8HDgReGREjkh4ButN5Q3XLVQEP6209Aj4UEdfs1iidQLITVM9/Gtw81gFfA04AFta1fxH4c0Scke7o/qVu3vaMcrNZJGkhcBLJOYkBlEm2xd+z5210Z0RUs8vSZsm0tuuI2CHpOuA04K2AL6BSbKPsPuKve8L8XftkVbx/2wqm6m9rUR7WO7U+4Km0MD0RWJp3Qpapa4D3S6pAMjRM0tx03hpJy9Pz2M4kuRiLNYcfAl+IiLsntPcxfiGVczPNyBrlzcBPImJpRCyLiMXAwyRHSb2NtpaZbNeXAt8C1ntkROE9CqyS1CWpDzg574SsodzfbcrF6R6kV+wcAn4GrJZ0K8lR1H/mmphl7VLgPuB2SfcA32P8F9mbgQuBe0h2hq+e9BUscxGxISIummTWV4AvS7qJ5AibFd9a/n/b+xXwNryNtpSZbNcRcRuwFfhRBilaA+zaH4uIx4ArgbtI9s3uyDUxawj3tynCIxEnI+kVwPcjYk3euZiZ2fSkQ+/Pj4g35J2L5UfSISTDfFdGRC3ndGwGvD/WXtzf5iOnk5D0PpLL1H8m71zMzMxs+iS9g+Sqvp92YVpM3h9rL+5vAx85NTMzMzMzsybgI6dmZmZmZmaWOxengKTFkv4s6X5J90r6SNq+v6TrJP07vV+Qtq+UdLOkIUnnT3itj6WvcY+ky9M/DDczMzMzM7O9cHGaGAU+HhEvBY4GPiBpFfBJ4PqIWAFcnz4GeAb4MMn/rY2R1J+2r46Iw0iuGnhWNqtgZmZmZmZWXC5OgYh4IiJuT6e3AfcD/SR/3H1ZuthlwOnpMk9FxHpgZJKX6wB60kthzwEeb3D6ZmZmZmZmhefidAJJy4AjSa7wd1BEPAFJAQss2ttzI2IjydHU/wJPAFsi4tpG5mtmZmZmZtYKXJzWkTSP5M/bPxoRW2fw/AUkR1uXA4cAcyWdPbtZmpmZmZmZtR4XpylJFZLC9GcRcVXa/KSkF6TzXwA8NcXLnAI8HBFPR8QIcBXwmkblbGZmZmZm1ipcnAKSBPwAuD8ivlE3ax1wTjp9DvDrKV7qv8DRkuakr3kyyfmrZmZmZmZmtheKiLxzyJ2k44AbgLuBWtr8KZLzTq8ElpAUnm+JiGckHQzcCvSmyw8AqyJiq6TPA2eSXAH4DuBdETGU5fqYmZmZmZkVjYtTMzMzMzMzy52H9ZqZmZmZmVnuXJyamZmZmZlZ7lycmpmZmZmZWe5cnJqZmZmZmVnuXJyamZmZmZlZ7lycmpmZmZmZWe5cnJqZmZmZmVnu/gc/sDp0YLz1lgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 936x360 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from scipy.signal import lfilter\n",
"\n",
"exog = np.zeros(200)\n",
"exog[40:80] = 1\n",
"exog[120:] = 1\n",
"eps = np.zeros_like(exog)\n",
"eta = np.random.normal(scale=0.1, size=len(exog)) * 0\n",
"endog = lfilter([1], [1, -0.85], exog + eta)\n",
"\n",
"index = pd.PeriodIndex(start='2018-01-01', periods=len(endog), freq='D')\n",
"endog = pd.Series(endog, index=index)\n",
"exog = pd.Series(exog, index=index)\n",
"\n",
"fig, axes = plt.subplots(2, figsize=(13, 5))\n",
"\n",
"exog.plot(ax=axes[0], title='exog')\n",
"endog.plot(ax=axes[1], title='endog')\n",
"\n",
"fig.tight_layout();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although Statsmodels doesn't offer a built-in class with the model you're looking for, I can set one up to demonstrate how it would work.\n",
"\n",
"In the state space setup, the key here is that the `exog` array needs to affect the `state_intercept` and not (as in the `SARIMAX` class) the `obs_intercept`.\n",
"\n",
"Another important aspect here is that we need to use diffuse initialization here, since the `exog` array is non-stationary, and it is determining a time-varying mean in the state process."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"class AR1X(sm.tsa.statespace.MLEModel):\n",
" _start_params = [0., 0., 1]\n",
" _param_names = ['phi', 'beta', 'sigma2']\n",
"\n",
" def __init__(self, endog, exog):\n",
" super(AR1X, self).__init__(\n",
" endog, k_states=1, exog=exog, initialization='diffuse')\n",
" \n",
" # For demonstration purposes, it's easier to just support k_exog=1\n",
" if not self.exog.shape == (self.nobs,):\n",
" raise NotImplementedError\n",
"\n",
" self['design', 0, 0] = 1.\n",
" self['state_intercept'] = np.zeros((1, self.nobs))\n",
" self['selection', 0, 0] = 1.\n",
" \n",
" def transform_params(self, params):\n",
" return np.r_[params[0:2], params[-1]**2]\n",
" \n",
" def untransform_params(self, params):\n",
" return np.r_[params[0:2], params[-1]**0.5]\n",
"\n",
" def update(self, params, **kwargs):\n",
" params = super(AR1X, self).update(params, **kwargs)\n",
" \n",
" self['state_intercept', 0, :] = self.exog * params[1]\n",
" self['transition', 0, 0] = params[0]\n",
" self['state_cov', 0, 0] = params[-1]\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can fit the model via MLE:\n",
"\n",
"- Parameter estimates are close to what we specified.\n",
"- Dynamic forecasts work pretty well."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Statespace Model Results \n",
"==============================================================================\n",
"Dep. Variable: y No. Observations: 200\n",
"Model: AR1X Log Likelihood 134.405\n",
"Date: Wed, 15 Aug 2018 AIC -260.811\n",
"Time: 19:53:38 BIC -247.617\n",
"Sample: 01-01-2018 HQIC -255.471\n",
" - 07-19-2018 \n",
"Covariance Type: opg \n",
"==============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"phi 0.8479 0.137 6.200 0.000 0.580 1.116\n",
"beta 1.0044 0.911 1.103 0.270 -0.781 2.790\n",
"sigma2 0.0150 0.000 46.568 0.000 0.014 0.016\n",
"===================================================================================\n",
"Ljung-Box (Q): 111.96 Jarque-Bera (JB): 33081.18\n",
"Prob(Q): 0.00 Prob(JB): 0.00\n",
"Heteroskedasticity (H): 0.01 Skew: 2.58\n",
"Prob(H) (two-sided): 0.00 Kurtosis: 65.95\n",
"===================================================================================\n",
"\n",
"Warnings:\n",
"[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwIAAAEiCAYAAABKj2WnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8XFd9///XmX20r7ZkybZkO5YdrwHHJBCysSUkJJQ9pRDW8IWGQku/hfRbCOUHpbRAKW2/paE0EEIIS1gSEpbwIwmEJE6cxXYWy5tkW7Jl7dJoNOu95/vHjGR5k2VHoxlJ7+fjMY9Z7rlzP6OHPXM/95zzOcZai4iIiIiIzC+efAcgIiIiIiIzT4mAiIiIiMg8pERARERERGQeUiIgIiIiIjIPKREQEREREZmHlAiIiIiIiMxDSgRERGYJY0yTMcYaY3xTaPseY8zDZ3mcs95XRERmDyUCIiI5YIxpN8YkjTE1x73+TPZkvik/kYmIiGQoERARyZ024LqxJ8aYdUA4f+FMv6n0ToiISGFSIiAikjvfBd494fn1wG0TGxhjyo0xtxljeowx+40xf2eM8WS3eY0xXzbG9Bpj9gFXnWTfbxljDhtjOo0xnzfGeE8WiDFmlTHmfmNMvzGm1Rjztgnbqo0xdxtjho0xjwPLT/WBJgxPer8x5gDwu+zrFxhjHjHGDBpjthljLp2wz3uMMfuMMRFjTJsx5p0TXv+jMebfjDFDxpidxphXTdhvUTaufmPMHmPMByds+6wx5ofZv13EGPOcMWbThO2fzP5NItnP+6rs6x5jzKeMMXuNMX3Z96g61ecVEZnLlAiIiOTOY0CZMWZ19gT97cDtx7X5N6AcWAZcQiZxeG922weBq4HzgE3AW47b9ztAGliRbfNa4APHB2GMKQbuB+4AFpDppfi/xpg12Sb/AcSBeuB92dvpXAKsBl5njGkA7gU+D1QBfw3cZYypzR7768CV1tpS4OXAMxPe52XAPqAGuBn4yYQT8+8DHcCi7Gf/h4mJAnANcCdQAdwN/Hv287YANwLnZ4/5OqA9u89fAG/Mxr8IGMh+fhGReUeJgIhIbo31CrwG2Al0jm2YkBzcZK2NWGvbga8A78o2eRvwNWvtQWttP/DFCfsuBK4EPm6tjVpru4F/Ad5xkhiuBtqttbdaa9PW2qeAu4C3ZGN4M/CZ7Ps8SybBOJ3PZtvHgD8D7rPW3metda219wNbgddn27rAWmNM2Fp72Fr73IT36c5+xpS19gdAK3CVMWYxcBHwSWtt3Fr7DPDfE/42AA9nj+lk/84bsq87QBA41xjjt9a2W2v3Zrd9CPg/1toOa20C+Gz276AhTiIy7ygREBHJre8Cfwq8h+OGBZG5Ch4A9k94bT/QkH28CDh43LYxSwE/cDg7HGcQ+C8yV/yPtxR42Vi7bNt3AnVALeCb5DinMrH9UuCtx73/RUC9tTZKJtn5X9lY7zXGrJqwb6e11h537EXZW7+1NnLctoYJz7smPB4FQsYYn7V2D/BxMif53caYO40xiybE+tMJcb5AJnFYOIXPLCIypygREBHJIWvtfjKThl8P/OS4zb1AiszJ6ZglHO01OAwsPm7bmINAAqix1lZkb2XW2jWc6CDw0IR2FdbaEmvth4EeMsOLTnWcU360497/u8e9f7G19h8BrLW/tta+hszQo53ANyfs22CMMccd+1D2VmWMKT1uWydTYK29w1p7EZm/rQW+NCHWK4+LNWStndL7iojMJUoERERy7/3A5dmr4+OyQ1p+CHzBGFNqjFkK/BVH5xH8EPgLY0yjMaYS+NSEfQ8DvwG+Yowpy06CXW6MueQkx/8FsNIY8y5jjD97O98Yszobw0+Azxpjiowx55KZ1HwmbgfeYIx5XXaCc8gYc2k27oXGmGuycwUSwAiZK/BjFmQ/o98Y81Yy8w7us9YeBB4Bvph9v/XZv+P3TheMMabFGHO5MSZIZu5DbMIxv0Hm770027bWGHPtGX5eEZE5QYmAiEiOWWv3Wmu3nmLzR4EomQmzD5OZ0Ps/2W3fBH4NbAOe4sQehXeTGVr0PJlJrz8mc9X9+ONHyEwkfgeZK+1dZK6QB7NNbgRKsq9/G7j1DD/fQeBa4G/J9DAcBP43md8YD/CJ7HH7yUzS/ciE3bcA55DpHfkC8BZrbV9223VAU3bfnwI3Z+cfnE4Q+Mfse3aRSTb+NrvtX8lMLP6NMSZCZkL3y87k84qIzBXm2KGZIiIiM8MY8x7gA9khPCIiMsPUIyAiIiIiMg8pERARERERmYc0NEhEREREZB5Sj4CIiIiIyDykREBEREREZB4qqCXVa2pqbFNTU77DEBERERGZtZ588slea23t6doVVCLQ1NTE1q2nKrUtIiIiIiKnY4zZP5V2GhokIiIiIjIPKREQEREREZmHlAiIiIiIiMxDSgREREREROYhJQIiIiIiIvOQEgERERERkXkoZ4mAMabFGPPMhNuwMebjuTqeiIiIiIhMXc7WEbDWtgIbAYwxXqAT+GmujiciIiIyE6y1mXvsCa8d//rx207V5lTtprLfmZjKMeTFMcYQ8oXyHcaUzNSCYq8C9lprp7S4gYiIiMhUWWuxWFzrYm32fsJzix1//VQ3YPzxxH0AXLLb3cw95uhxTfaJxWKMOeH147dNjHlim8nanuzznmzf0zlZEnG6Y8mZ8+BhRfWKfIcxJTOVCLwD+P7JNhhjbgBuAFiyZMkMhSMiIiKFzFqLYx0c18GxDq51STtpUm6KtJsev421GT8xHjuvtRNOqm3mdYPBGDPeduzx2Mnw+PaTvD7WXuR0RhIj+Q5hynKeCBhjAsA1wE0n226tvQW4BWDTpk3qrxIREZknXOuOn9CnnBQJJ0EinSDlpEjb9PgJ/NgVcGMMHuMZv/k8PgImoBN0kbM0Ez0CVwJPWWuPzMCxREREpACl3TRJJ0kynSSWjhFLxUi76fHtYyf5Po+PgC9AyMyOMdYis9lMJALXcYphQSIiIjL3WGtJuSniqTjRVJTR1CiO64wPz9HJvkhhyGkiYIwpAl4DfCiXxxEREZH8clyHeDpOJBFhJDUyPm7f5/UR9AXxGC1dJFJocpoIWGtHgepcHkNERETyw3EdYukYg7FBoqkoBoPX4yXkC+nEX2QWmKmqQSIiIjIHWGuJp+MMxAcy1VEMBLwBSoOl+Q5NRM6QEgERERE5Lcd1GEmO0DfaR8pN4ff6KQ4Uq2KPyCymREBEZjXXtcRSDsm0SyLtZu8dEsc9H9s+9loy7ZB2LWnX4riWtGNJu+7485Tj4riWmpIgH718hU52ZN5Ku2mGE8P0jfZhrSXkDxHyz55Jvq61JFIuKceSSFuS6bHHmfvM90H2fsLzlJP5Lhi/WcYfuxMeOy7HtHOPa5dZnAysJXtvJzyesP3453Ys/gkrGWf3c7MPTraPJftmJzHZqsST1W+fbDHiSeu+n+V+k61+nIvjTbcVCwPcev3KGTzi2VMiICJ55biW4ViK/tEkA9Ek/dEkA6NJ+qMphmIpRhIpogmHSDxNNJEmmkwzEk8zkhh77kxbLB4DPo8Hr8fg8xgcaxlNOlyzYRFNNcXTdhyR2cBxHYYSQ/RGezHGEPaH8zLu37WWkbjDUCzN4KjD0GiawVia4ViaoZjDaMIhmnQYTbiMJh2i2fvRpEs04RBLutN6Eugx4PUYvB5zzGPvxMceg8cDHmMwZmyNM4sHwGRO+41xsxtczNh9drvBjrfzmOw+1mbeyxgMFjO+3WQWToOj+xvD0VPfiaf/9uh6ayecbFuOXu84w7/YMcc7mZNvm+zyymTbrDn1sSZ9zxm6oFNaFJ6R40wHJQIikjPWWvqjSdp6o3QMxOgajtM1lLkdHo5zZChOdySOe4rv9IDXQ0nIR3HQS3HAR2nIR1VxgMVVRZQGfRRnb0UBL0Gfh6Avcx/weTLP/V4CXg9Bv4eA10PIn2kT8GWe+7zmmBN/j+fYH4lH9/Zx3Tcfo2MgpkRA5g1rLSPJEbqj3TiuQ1GgKKcJwHAszYH+BIcGE/RGUvREUnQPZ+57Ikl6R1I47sn39XrIfj94KQp6CPs9FAUTlJck8PpieL0JvN4EmDjWJHBsHJc4Dklcm8CxSRybwCWVeewmcWwax6ZIZx+n3BSOmyZt05n77C3hOqRt5vHYwmiu62ZWOnYc3LSLazM3x3UmvRovc8v5i84H3pTvMKZEiYCIvGiua2nri/LC4WH29URp642yrzdKW88Iw/H0MW1Lgz4WloeoKwtxzjk1LCwLUlMSpKo4QGVRIHNfHKCqKEA44M3TJ8porMxc1ekYGM1rHCIzJekk6RrpYjQ5SlGgiLB/eq5sWms5NJhk15EYB/viHOhPcKAvwcH+OEOxY3v1wn4PtWV+FpT62bAkSDA4ivEN4jBIigES7iBxZ5BoepDhxABDiSEOx4cYig8xPDKMa0+RNZyEx3gI+UIEvAFC3hB+r5+AN0DAGxh/7Pf4CfmK8Hv8+Lw+fMaH1+PF7/Hj9XjxGm/mPvvYYzzHPD7+lrmin6mu5OHo87EF1YwxJ7wOjD+f+PiY7ZzY7ujdie2Ofzy276m2nardMa9Pdj1+mi/GT3qsPCvyFeU7hClTIiAiZ8RaS+dgjO0dQ9nbIDs6h4hMOOFfVB6iubaYazYuormmhGW1xSyuDFNXHqYkOHu+durLQ3g9ho6BWL5DEckpay2D8UG6o934vX7KQmUv6v16hpO8cHh0/Lbz8CjDE074a0v9LK4KcklLKSUlAxjfIUbdToaShzkcPUjHcAfPDx+m70jfSd+/IlhBZbiSynAlNeEallcupzxYTmmwlLJgGSWBEkoDpRQHiin2F1PkL6LIn0lswr4wIV+IkC9z4i8y3UYSI/kOYcpmzy+yiOTNwf5R/rinlz/u7ePRvX30jiQA8HsNq+rKeMOGRWxoLGfNonKW15bk/Ur+dPF5PdSXhzioHgGZw1JOKtMLkBqlOFB8VsOAuoeTbG2P8GT7CE+2R+iOpADwGlhWG+aSlgoaq2O4vn30p/awd6CVF3pe4L49bSScxPj7hHwhFpctZnHZYjbWbaSupI76knoWFi+ktriWBcULqApX4fPo9EVkOuh/koicwHUt2zoG+c3zR7j/+SPs6c5c3agtDfKKFdVsWlrJusYKVtWVEvLPjZP+U2msDKtHQOasWCpG53AnxpgzWgfAWsuurhh/2D3EH3YNsftI5v9IRZGPly4tYU1jmHBRB13xZ3i6ays/P/wMB1sPju/fUNrAqppVXNZ8Gcsrl7OschnNlc3UFtWqQpfIDFIiICLjdh+J8JOnO/n5050cGorj9Rhe1lzFdZuXcPE5NaxYUDLvfqQbK4v4w+6efIchMq3GhgIdGTlC2B+e8hCZg/1xfv3sAL9+tp/OgSQGWL+4mI9cvogV9TH2RR7hgbYH+MnWRxiMDwKZk/7z6s/jPRvfw4aFGzi39lzKQ+U5/HQiMlVKBETmuXjK4e5th7j9sf1s7xjC6zFcfE4N//uKFi5vWUh50fweQ9tYGebIcIJE2iHom9u9HzI/WGvpjnYzEB+gJFhy2qFAKcflgZ2D/GRrL9s7ohjgpU0lXP/yOprrhnnwwH3csftetj26DYD6knquWH4FFy6+kAsbL6ShrGEGPpWInA0lAiLzVHckzrf/2M73Hz/AwGiKlQtL+MzV5/KGDYuoLQ3mO7yCsbgyU/3h0GCcZpUQlVnOtS5HRo4wnBimLDj5hODhWJofb+3hZ0/10juSprEyyEcuX8RFKwP8sfM+vvn8D9n6u60AnFd3Hp98xSd59bJXs7pm9bzrORSZrZQIiMwzPZEE//XQXm7fsp9k2uU15y7k+pc3ceGyav14n8RYCdGD/aNKBGRWc63L4chhoqnopPMBhmNp7tzSzQ+f6GE06XLBslI+dVUtVRVdfHfb1/n89+8ikozQUt3Cpy76FNesvIalFUtn8JOIyHRRIiAyT0QTaf7vg3v41sNtJNMubzyvgY9efo5Obk+jsSrTI6AJwzKbudalc7iTeDpOSaDkpG0SaZc7t3Rz+6NHiCZcLltVwfteWceIu4uvb/kiv9r7K4LeIFetvIp3b3g3m+o36eKByCynREBkjrPWcs/2w/zDvS/QNRzn2o2L+NirzmFZ7clPBuRYdWUhfB6jRcVk1rLWcmTkCLFUjJLgyf/f/3H3EF+7v4POgSSvXFnODZfUEzf7uPkPH+SB9gcoD5bzlxf8Je87731Uhatm+BOISK4oERCZww72j/I3P97Oo/v6WLOojH//0/PY1KQf8TPh9RgWVaiEqMxO1lp6R3sZTgyfdDhQ30iKL913gId3D7O0OsjXrltO88IE//zIzdyx4w4qQhXcdNFNXL/h+jMqLyois4MSAZE5yFrLXU918tm7nwPg829cy3Wbl+D1qBv/bDRWhrWomMxKA/EB+mP9Jx0O9ODOQb503wFiKZcbX7WIN2+q5vbt3+Fd9/0T8XScD770g3zsZR+jIlSRh8gLg7U2c4+d9PnE145//fhtp2ozpXhO8j5T3vcsjidzX04TAWNMBfDfwFrAAu+z1j6ay2OKzHdDoyk+edd2fvVcF5ubqvjK2zawODvOXc5OY2WYB1q1loDMLpFEhO5oN6WB0mPG8seSDl/5dQf3be9nVV2Yz1zbBL4u/vSut/NY52Nc3nQ5n73ssyyvXJ6/4M+Aa11c62KtzdxjsdZiyTzHAobT3ltrMZhjHns8mdKqHjL3Y3/HsZKrxpjMPhNeO/51AIM5YT7F8W0mazuV/abqdO8vL57XM3tKTee6R+BfgV9Za99ijAkAOhsRyaE93RE+eNuTdAyMctOVq/jAK5epF2AaNFYW0RNJEE85c34lZZkbkk6Sw5HDFPuLjznpOzyU5FM/2seeIzHe84qFvPeiOr7/3O187qHP4fP4+Orrvsrbzn1b3k8UrbU41sFxHVzr4lhn/PXjT4J9Hh8e48Hr8RLwBvAYz/hrY6+PnfxOdg+c8FhkrstZImCMKQMuBt4DYK1NAslcHU9kvntgZzd/8f2nCfo93PHBCzhfcwGmzeKqTAnRzsEYyzXJWgqca10ODR/C5/Udc2Xy6QMj/J+72kg7lq+8Yzkblvj5xP0f564X7uKSpZfw5dd+mUWli2YsTmstaTdN2k0fcyV/7Eq83+Mn5AsR8AbweXzHnPCPneSfbjE0EZlcLnsElgE9wK3GmA3Ak8DHrLXRiY2MMTcANwAsWbIkh+GIzF23/rGNz/3iec6tL+OWd2+ioSKc75DmlMbsomIH+0eVCEjB64n2kHJTFAeOlgb+5Y5+/uEX+2moDPKlty7D4+vh2jvfz3M9z/GJCz/Bxy/4eE5Pql3rknJSpNzUMSf7QW+QsmAZIV8InyeTuIyd8ItI7uUyEfABLwE+aq3dYoz5V+BTwKcnNrLW3gLcArBp0ybNZBE5Q//xwB7++detvG7NQr729vMIBzR0ZbqNLSqmykFS6CKJCAPxgWNWDf7ZU7380y8PsqmphH948zJa+5/h+h9cj2Mdvv3Gb/OaZa+Z9jgc1yHpJEm76fGT/mJ/MVX+KgLeAH6vH59H9UpE8i2X/ws7gA5r7Zbs8x+TSQREZBpYa/nKb3bx7w/s4U/Oa+Cf37Ien1dX0XJhQWkIv9coEZCClnbTdI10Uew/2hNw55Zuvv7bTl6xoozPv7mZLZ0P8/67309tUS23v+l2llUum7bjJ50kyXQSi8Xv8VMaLKXYXzx+4i8ihSdniYC1tssYc9AY02KtbQVeBTyfq+OJzCfWWv7hvhf45h/auG7zYr7wxnV4NCk4Z7weQ0NFWIuKSUHrifZgjBmfF3DbI11844HDXL66gpuvXcpv9/2aj9z3EZZXLud7b/oeC0sWvuhjJp0kSSeJtTZzxb+kipA/M65fRApfrvvlPgp8L1sxaB/w3hwfT2ReuOX3+/jmH9p494VL+ftr1qi6xQxorCzioHoEpECNJEYYig9RFsoMCfrFM31844HDvHZNJX93zVLu3XU3N/7yRjbWbeS2N95GZbjyrI/lWpdYKoZrXcL+MHXFdRQFijTUR2QWyun/WmvtM8CmXB5DZL75+TOdfPGXO7l6fT2ffYOSgJnSWBnmty8cyXcYIicYGxJUFMhMan90zxBfuu8Am5eV8ndvWMofD/yej/3qY2xetJnb/uS2YyYRn4mUkyKejuPz+KguqqY0UKohPyKznNJ3kVnkkb29/PWPtvGy5sxCYRoONHMaK8P0jiSJJR1NyJaC0jfaB2Tq6e88PMrf/aSd5QvCfOFNzTzbs40P3PMBzqk+h1vfeOtZJQFJJ0k8HSfkDdFY1kjYH1ZVH5E5Qv+TRWaJfT0jfOi7T9JUXcwt79pE0KeT0Zk0tjpz56DmCUjhSKQTDMYHCfvD9Awn+esf7KW8yMeX376cw9E23vXTd1Edrub2P7n9mEpCU5F20wwnhjEYlpQvYWnFUooDxUoCROYQ9QiIzALxlMNHvvcUPo/h1veeT3mRuuNn2lgJ0YP9MVYsKM1zNCIZPdEe/F4/joXP/KydWNLl3965Ar9/lOt/eD0Gwx1vvuOMJga71mU0NYrP+Fhctpgif5GGIIrMUUoERGaBv7/neXZ2Rbj1veePL24lM2vs736gXz0CUhhGU6OMJEcoC5XxjQcOse1glJuvXUpTTYgP3H0jHZEOfvzWH59RidBYKobjOtQW11IeKtfVf5E5TomASIH7+TOdfP/xA3z40uVc1rIg3+HMWwtKgxQFvLT3RU/fWCTHrLV0j3QT8od4dM8Qtz1yhDdsrOZ1a6v4zyf+k1/t/RU3X3Iz5zecP6X3c63LSGKE4kAxdeV1mgQsMk8oERApYHt7RrjpJzs4v6mST7xmZb7DmdeMMTRVF9PWq0RA8i+SiJBwEsQSQT53935WLAjxV69t5LGOx/jiw1/kqnOu4oMv+eCU3iuRTpB0ktSX1lMWLNMwIJF5RH1+IgXKcS1//aNtBH0e/u26l2jV4ALQXKtEQPLPtS49oz2EfWG+dN9BEmnL//emZqLpQT5874dZWrGUr7z2K1M6oR9NZoa6NVc2Ux4qVxIgMs/ozEKkQN3+2H6ePjDIZ95wLnXloXyHI8CymmIO9o+STLv5DkXmsZHECGk3zQM7Izy6d5gPXVrP0uoQNz9wMwOxAb5x9TcoDU4+od1ay3B8mOJAMUvKl2glYJF5SomASAE6NBjjn361k4tX1vLGjQ35DkeymmuKca0mDEv+WGvpHe0lmfLzL7/pYHV9EW/ZVMv9++7nJzt/wl+87C9YU7tm0vdwrUskGWFByQLqSurwelSKWGS+UiIgUmCstXz6Z8/iWvjCG9eqq76ANNdkFmPS8CDJl9HUKCk3xX8+cIRIPM1NVy0hmorwqd9+itU1q7lx842T7u+4DiOJEepL6qkKV+n7RWSeUyIgUmDu3XGY/39nN5947crxRaykMBxNBEbyHInMR9ZaeqI97OhIcd/2ft55wUJWLAzz+d9/nu5oN19+7ZcnHeKTdtNEU1EayhooD5XPYOQiUqhUNUikgIwm03zunudZ11DOe17elO9w5DgVRQGqigPqEZC8iKVjRJNxvn5/F4urgrz3lXU8cvARvrfje3x404fZWLfxlPs6rkMsFWNJ+RKK/LrAICIZ6hEQKSD/83Ab3ZEEN7/hXFUJKlDNNcXs61EiIDOvb7SP3z4X5UBfghtf1YDfC5998LM0ljXyiQs/ccr9XOsSTUZpKG1QEiAix9CZhkiB6BtJ8I2H9vGacxeyqakq3+HIKTTXqISozLx4Ok5fdITv/LGXDYuLueicMu564S6e63mOmy66ibA/fNL9rLWZOQGl9ZQES2Y4ahEpdEoERArEv/1uD6PJNJ+8oiXfocgkmmuK6Y4kGEmk8x2KzCPDiWF++uQw/dE0H7m8gXg6zpce/hIbFm7gmpZrTrnfWHUgzQkQkZNRIiBSAPb3Rfnelv28/fwlrFgwef1vya9l2QnD7eoVkBniuA7tfX388Il+Lm0pZ11jMd986pscHjnMpy/+NB5z8p/ykcQIlaFKqsLqYRSRk1MiIFIA/vnXrfg8Hv7y1efkOxQ5jebaTCKwT4mAzJBoMsrtj/aRTLl86LJF9I728h9P/AevW/46Llx84Un3iafjBLwBaotrZzhaEZlNclo1yBjTDkQAB0hbazfl8ngis9ELh4f5xfbDfPTyFSwo0wrCha6pOltCVBOGZYa8cKSLe7cNcc15NSytDvF3v/s8sVSMv33l3560fdpNk3bSNFU2nbK3QEQEZqZ86GXW2t4ZOI7IrHTL7/dRFPDygYuW5TsUmYKQ30tDRVhrCciMSKQT3Pl4NwDvevlCuqPd3LHjDt6+5u2sqFpxQntrLaOpURpLGyddU0BEBDQ0SCSvOgZGuXvbIa7bvITyIn++w5EpaqopUuUgmREHBvv45fZhXrOmkrryAN966luk3BQfPv/DJ20/mhqlJlyjCkEiMiW5TgQs8BtjzJPGmBtyfCyRWed/Hm7HAO+/qDnfocgZGCshaq3Ndygyhzmuw/ceO0A8ZXnnBQsZTgzznW3f4apzrmJZ5Yk9iEknidfjpapIk4NFZGpynQi8wlr7EuBK4M+NMRcf38AYc4MxZqsxZmtPT0+OwxEpHIOjSe584gDXbFzEooqT1wCXwtRcU8JwPE1/NJnvUGQO6x+N8NMnB7lweRnLF4S5bdttRJIR/vz8Pz+hrbWWeCpOfUm95gWIyJTl9NvCWnsoe98N/BTYfJI2t1hrN1lrN9XWqrqBzB/ffXQ/o0mHGy7W3IDZZqyEqIYHSS7d+UQbQzGHP7twAbFUjP9+6r+5ZOklrFu47oS2o6lRqouqT7mwmIjIyeQsETDGFBtjSsceA68Fns3V8URmk3jK4duPtHNZSy2r6sryHY6coeYalRCV3EqkU9yxpYc1i4rYuKSEHz3/I3pGe7hx840ntE05KbzGq/UCROSM5bJq0ELgp8aYsePcYa39VQ6PJzJr/PTpTvqiST50yfJ8hyJnobEyjM9j1CMgOXPv9oMcHkxx46saca3LN7Z+g/PqzuPCxhPXDYilYiypWILX481oxWhIAAAgAElEQVRDpCIym+UsEbDW7gM25Or9RWazO7YcYFVdKS9r1hW82cjn9bCkuoh9PSohKrlx5xMHWVDm5+KV5Ty0/wH2D+3nplfeRPbi2rhYKkZpsJQif1GeIhWR2UwzikRm2LOdQ+zoHOK6zUtO+FGX2WNFbQm7jygRkOnX1jvME20jXL2hGq/H8L3t36O2qJYrll9xTDtrLSk3RU1RTZ4iFZHZTomAyAz7wRMHCfg8vHFjQ75DkRdhVV0p7X1R4ikn36HIHHPnE+0AXLWhmq6RLu7fdz9vW/M2/N5j1xqJpWJUhasI+oJ5iFJE5gIlAiIzKJZ0+Nkznbx+bZ0WEJvlWurKcC3s6VavgEwfx7X85Kkuzl9WQn15gB889wMc63Dd2uuOaedaF9e6miAsIi+KEgGRGXTfjsNE4mnesXlJvkORF6mlLrNya2tXJM+RyFzyu52H6YmkuPa8Wlzr8v0d3+eiJRfRXHnsooOjyVFqi2vxeXJZ80NE5jolAiIz6AdPHKSpukiThOeApupiAj4PrUeUCMj0uePx/VQUebnonDL+sP8PHBw+yDvXvfOYNo7r4DEeykPleYpSROYKJQIiM2RP9wiPt/fz9vM1SXgu8Hk9rKgtYad6BGSadA/Heai1nyvXVeL3erh9x+1Uh6u5YsWxk4RjqRjVRdVaQVhEXjR9i4jMkB9uPYjPY3jzSzVJeK5oqStllxIBmSY/fPIAroVrz6ulO9rNb/b+hree+1YC3sB4G9e6GGMoC2ohQhF58ZQIiMwA17X87OlOLlu1gAWloXyHI9Okpa6UruE4Q6OpfIcic8DPnz7EmoYQS6pD/HTnT0m7aa5bd+wk4bFKQVo8TESmgxIBkRmwdf8A3ZEEV6+vz3coMo1a6koBNE9AXrQ93RF2d0e5bHVm3P89rfewbsE6VlStGG9jrcW1ruYGiMi0USIgMgPu23GYgM/Dq1YvzHcoMo1aFmYTga7hPEcis9292w9jgFevruHg0EGe7nqaN6x8wzFtYukYleFKVQoSkWmjREAkx1zXct+Ow1zWUktJUD/gc0l9eYjSkE8ThuVF+8WOQ6xpDLGgLMAvdv0CgDe0HE0ErLU4rkNFqCJfIYrIHKREQCTHxoYFvX6dhgXNNcYYVtWVsktDg+RF2NMdYfeRKJe0ZHqY7tl1DxsXbmRJ+dH1RuLpOGXBsmMmDouIvFhKBERy7L4dhwlqWNCctXJhKTu7Ilhr8x2KzFL3bu/CAK86t5r2wXa2Hdl2TG8AQNpNqzdARKadEgGRHBobFnSphgXNWavqSonE03QNx/MdisxSv9h+iLWNIerKwuPDgq5eefX49rSbxu/1E/Kp4piITC8lAiI5NDYs6Kr1i/IdiuRIS12mnrvmCcjZ2H0kwu7uEV65sgSAu1vv5iX1L6GxrHG8TTwVpzpcrYUIRWTaKREQyaHxYUGrFuQ7FMmRo5WDlAjImbt3R6Za0GWrK9g7sJfnep47plqQtRaLpThQnL8gRWTOUiIgkiMThwUVa1jQnFVe5KeuLKQVhuWs3LfjMGsbwywqL+Ge1nsAuGrlVePb4+k45cFylQwVkZxQIiCSI9s7h+iOJLhyraoFzXUtdaUaGiRn7GD/KLuOjPDyc4oxxvCbvb/hJfUvoaG0YbxN2k1rATERyZmcJwLGGK8x5mljzC9yfSyRQvJQaw/GwMUra/MdiuRYS10pe3pGSDtuvkORWeTB1m4ALlhWQu9oL9uObOPy5svHtzuug9+jScIikjsz0SPwMeCFGTiOSEF5cFc36xsrqCpW3e+5bs2iMpJpl93dI/kORWaRh3b1sKgiwNLqMA+2PwjA5U1HE4FYKkZ1kSYJi0ju5DQRMMY0AlcB/53L44gUmoFokm0HB7lUvQHzwrqGzNCN7R2DeY5EZotE2uGRvX28tClM0BfkgbYHqCmqYd3CdeNtLJYif1EeoxSRuS7XPQJfA/4GOGV/uTHmBmPMVmPM1p6enhyHIzIz/rCnF9fCJS1KBOaDpupiSkM+tncM5TsUmSWeaBtgNOlwfnMRrnV5cP+DXNp0KR6T+VlOOSnC/jB+rz/PkYrIXJazRMAYczXQba19crJ21tpbrLWbrLWbamt10iRzw0OtPVQU+dnQqJVA5wOPx7CuoVyJgEzZg63d+L2GDYuLeLrraQbjg8cMC0qkE1QE9f0hIrmVyx6BVwDXGGPagTuBy40xt+fweCIFwXUtD+3q4ZXn1OL1aGzvfLG+sYKdXcMk0k6+Q5FZ4KFdPWxcUkJZODMsyGM8XLz04vHtFkvYH85jhCIyH+QsEbDW3mStbbTWNgHvAH5nrf2zXB1PpFA8f3iY3pGE5gfMM+sby0k5lp2HVUZUJtc5GGN39wgvbQrh9/p5oP0BXlL/EirDlQAknaSGBYnIjJhSImCM+dhUXhORzJU+UNnQ+WZ9Y3bCcKeGB8nkxsqGvrQpTN9oH9uObOOypsvGtyfTSQ0LEpEZMdUegetP8tp7pnoQa+2D1tqrp9peZDZ7qLWHtQ1l1JYG8x2KzKCGijBVxQG2H1TlIJncQ6091JcHWVzl58H9DwIcs36AxVIUULUgEcm9SdcsN8ZcB/wp0GyMuXvCplKgL5eBicxGQ7EUTx4Y4MOXLM93KDLDjDGsbyxnh3oEZBLJtMsf9/Ty2jWVBHyB8bKhaxeszWzPDgvyeSb9eRYRmRan+6Z5BDgM1ABfmfB6BNieq6BEZqtH9vTiuFZlQ+ep9Q3l/H5XD6PJNEUBncjJiZ7cP0A06bCxKYDP4+Oh/Q/x6mWvHi8bmnSS1BXX5TlKEZkvJv2lstbuB/YDF85MOCKz25a2fooCXjYu1vje+Wh9YwWuhecPDbOpqSrf4UgB2tLWhzGwtiFIa18rg/FBXrnklce00bAgEZkpU50sHDHGDGdvcWOMY4wZznVwIrPNlrZ+XrKkEr8312v1SSFal50wvE3rCcgpPNHeT8vCEkrDPrZ0bAHggsYLgMwiYkFvUMOCRGTGTOlsxVpbaq0ty95CwJuBf89taCKzy1Asxc6uYTY360rwfLWwLMTCsiA7OjRhWE6UTLs8uX+ADUuK8Rovj3U8RkNpA41ljZntTpKyYFmeoxSR+eSsLltaa38GXH7ahiLzyJP7+7EWJQLz3PrGCq0wLCf17KEh4imX1fV+fB4fWzq38LLGl41vd61LkV/DgkRk5kyp/9EY86YJTz3AJsDmJCKRWWpLWz8Br0fzA+a59Q3l3P/8EYbjKcpCWhBKjnq8rR+AVYt87B/aT+9oLxc0ZIYFudbF6/ES8AbyGaKIzDNTHYj4hgmP00A7cO20RyMyiz3e1s/6xnJCfm++Q5E8Wp9NBJ/tGOLlK2ryHI0Uksfb+mmuKaKqxM/Pd2bmB4z1CCTSCUr8JRhj8hmiiMwzU0oErLXvzXUgIrPZaDLNjo4hbrh4Wb5DkTzbkJ0w/OT+ASUCMs5xLU+09/PaNTVg4bHOx6gtqmV55fLsdofSYGmeoxSR+WaqVYOWGWPuMcb0GGO6jTE/N8bojEck65kDg6Rdq/kBQkVRgFV1pTze3p/vUKSAtHZFiMTTrGsM4/V42dKRmR8w1gNgsQR9Wo1cRGbWVCcL3wH8EKgHFgE/Ar6fq6BEZpstbf14DLx0aWW+Q5ECsLm5iif3D5B23HyHIgXi8bY+AFrqvXRHu+mMdI7PD0g5KUK+kMqGisiMm2oiYKy137XWprO329FkYZFxj7f1c+6iMko1OVTIJAKjSYfnDmm5Fcl4vL2fRRUhqksMTxx6ApgwP8BJqGyoiOTFVBOBB4wxnzLGNBljlhpj/ga41xhTZYzRWAiZ15Jpl6cPDrC5qTrfoUiB2JxdVXisSozMb9ZaHm8bYNPSCiyWxzoeozxYzqqaVePbVTZURPJhqv2Qb8/ef+i4199HpmdA8wVk3trRmakNvrlZw4IkY0FZiKbqIra09fNBTSCf99p6o/SOJNi4pASP8fBYx2NsbtiMx3hUNlRE8mqqicBqa2184gvGmNDxr4nMR2NXfc9vUueYHLW5uYpfP3cE17V4PCoJOZ+NfUecuyhIX6yLtsE2/mz9nwGZ1YRVNlRE8mWqQ4MemeJrIvPO1vZ+ltcWU12iih9y1ObmaoZiKXZ3j+Q7FMmzJ9oHqC4OUFvmsK1rGwCbGzYDmYnCJcGSfIYnIvPYpD0Cxpg6oAEIG2POA8YuWZQBGtAo8561lm0dQ1zaUpvvUKTAvKx5bJ5AHy11qg8/n23vGGR9YzkODju6d+D3+FlTu2Z8u4YFiUi+nG5o0OuA9wCNwFcnvB4B/nayHY0xIeD3QDB7nB9ba28+60hFCtDhoTi9IwnWZxeREhnTWBmmrizElrZ+3nVhU77DkTyJJtLs6RnhtWtqsNbyTNczrK5dTdAXHJ8f4Peo2piI5MekiYC19jvAd4wxb7bW3nWG750ALrfWjhhj/MDDxphfWmsfO9tgRQrN9o4hANY1KBGQYxlj2NxcxWP7+rDWagz4PPVs5xDWwqr6IiwJdnTv4NqWa4HMsKBif7H+bYhI3kx1svBaY8ya41+01n7uVDtYay0wNjjWn71p7QGZU7Z3DOLzGFbXqwa4nGhzcxV3bzvE/r5RmmqK8x2O5MGOzszFgmW1XjqHOxlODLNh4QYgkwhUh1V2WETyZ6qThUeAaPbmAFcCTafbyRjjNcY8A3QD91trt5xlnCIFaUfnEC11pYT83nyHIgVofJ5Au9YTmK+2dQyxqDxEOJjm+Z7nAdhQt2F8e9CnIgMikj9TSgSstV+ZcPsCcCmZScSn28+x1m4kM8dgszFm7fFtjDE3GGO2GmO29vT0nGH4IvljrWV7xxDrGyvyHYoUqBULSqgs8vPYvr58hyJ5sqNjkLUN5aTcFNu7txPyhVhZvXJ8uJgmCotIPk21R+B4RZzBImLW2kHgQeCKk2y7xVq7yVq7qbZWlVdk9jjQP8pQLKWJwnJKxhhevqKGP+zuxXU1MnK+GRpN0d43ytqGzDoB27q2sXbBWnweH0knSZG/SPMDRCSvppQIGGN2GGO2Z2/PAq3A10+zT60xpiL7OAy8Gtj5YgMWKRSaKCxTcenKWnoiCV7oGs53KDLDxuYHrKovIuWkeLb72fH5AWk3TUlA6weISH5NdbLw1UAl8EqgArjPWvvkafapJ1NxyEsm4fihtfYXZx2pSIHZ3jFIwOdRjXiZ1CUrMz2dD7b2sGaRksb5ZFvHIAArFgTY2b+fWDo2nghYrOYHiEjeTXVo0LXAd4EaMtV/bjXGfHSyHay1262151lr11tr105WYUhkNtreMcS59WX4vWc7wk7mgwVlIc6tL+OhXZoDNd/s6BhiaXURfn+KZ3ueBTIThTNF9bSQmIjk31TPYD4AXGCtvdla+xngQuCDuQtLpLA5ruXZziE2aH6ATMGlLbU8uX+A4Xgq36HIDMqsKFxBwknwXPdzlAZKWVa5jJSbIuwL4zG6iCAi+TXVbyFDpmzoGCf7msi81NY7QjTpsE4Vg2QKLm1ZgONa/ri7N9+hyAzpiSQ4NBRnzaJSHNdh25FtrFu4Do/xkHJSlAY0pFBE8m+qicCtwBZjzGeNMZ8FHgO+lbOoRArctoOZSYDqEZCpOG9JBaVBn4YHzSM7OjPzA85dlJko/HzP82xcuBEA17qaHyAiBWFKk4WttV81xjwIXESmJ+C91tqncxmYSCHb0TlEUcDLslpV/ZDT83s9XHRODQ+29ozXj5e5bXvHEMbAOQtCPNq5i5SbGl9ITOsHiEihmGrVIKy1TwFP5TAWkVlje8cgaxeV4/XohE6m5tKWWn75bBetRyKsqivLdziSY9s7hlhRW4LXlz46UXjhBhzXwWd8eD1ajVxE8k8zlUTOUNpxee7QMOs0LEjOwMXZMqIPtWp40Hywo3OIdQ3lxFIxnu95nspQJY1ljZmJwv5wvsMTEQGUCIicsfa+KIm0y5pFuqorU1dfHmZVXSkPKhGY8/pGEvREEqyuLyXhJGjta+Xc2nMxxuC4DsX+4nyHKCICKBEQOWOtXSMArFyoqh9yZi5tWcAT7f0MjaqM6FzW2hUBYMXCYlzrsqtvF6tqVgGZicIBn+YHiEhhUCIgcoZau4bxGFixQBOF5cxcsbaOtGv5zfNd+Q5Fcqj1SCYRWFYbpGOog9HU6HgiYDD4Pf58hiciMk6JgMgZaj0SoammmJBfk/3kzGxoLKehIsx9Ow7nOxTJodauCJVFfsrDht39uwFYVbMKx3Xwe/2aKCwiBUOJgMgZau2KsKpOw4LkzBljuGp9PQ/v6dXwoDms9UiElrpSYukYewb2ALCyeiVJJ0mRvyjP0YmIHKVEQOQMxJIO+/tHNT9AztpV6+pJORoeNFe5rmVXV4SWhaXE03F29+1mSfkSSgIlOK6jREBECooSAZEzsLs7grWoR0DO2vrGchorw9yr4UFzUudgjGjSYWVdCUknSWtf6/j8AAxaSExECooSAZEzMFYNRD0CcraMMVy1rp6Hd2t40Fw0XjFoQREpJ8W+gX20VLdkNlrwezVRWEQKhxIBkTPQ2hUh6POwtFp1wOXsvX5dPWnX8msND5pzxioGNdeE2Dewj7SbZnXNatJumoA3gMfoZ1dECoe+kUTOQOuRCOcsLMHrMfkORWaxseFBqh4097R2RWioCBP0W3b17wKgpaaFlJOiOKALCCJSWJQIiJyB1q4ILQu1orC8OBoeNHe1dmUqBiWcBHv69+D3+FleuRzHOoT94XyHJyJyjJwlAsaYxcaYB4wxLxhjnjPGfCxXxxKZCQPRJN2RBC11WkhMXryr1y8i7Vru2X4o36HINEmmXfb2jNBSd7Ri0IqqFZl5ARYtJCYiBSeXPQJp4BPW2tXABcCfG2POzeHxRHJqbOxvS516BOTFW9tQxqq6Un7wxMF8hyLTpK03Stq1rKorJekk2dW/6+hEYaOJwiJSeHKWCFhrD1trn8o+jgAvAA25Op5Iro1VA1HpUJkOxhjecf5idnQO8WznUL7DkWmws2sYgBW1xQwnhukY7mBVbXZFYY9fE4VFpODMyLeSMaYJOA/YMhPHE8mF1iMRysN+FpQG8x2KzBF/cl4jAZ9HvQJzxK4jEXwew5LqILv6shOFq1tIuSnCPs0PEJHCk/NEwBhTAtwFfNxaO3yS7TcYY7YaY7b29PTkOhyRszY2CdAYVQyS6VFe5Of1a+v42TOdxJJOvsORF6m1K8Ky2mI8Hnc8EVhdsxrH1URhESlMOU0EjDF+MknA96y1PzlZG2vtLdbaTdbaTbW1tbkMR+SsWWvZ1RWhRQuJyTR7x+YlROJplRKdA3Z2RVi5MDM/YHffbor9xTSWNWKxWlFYRApSLqsGGeBbwAvW2q/m6jgiM+HQUJxIIk2L5gfINHtZcxXNNcXc+cSBfIciL8JIIk3HQIxVdaUk0gl29++mpaYFYwzWWnweX75DFBE5QS57BF4BvAu43BjzTPb2+hweTyRndnWNVQxSIiDTyxjD289fzBPtA+zpjuQ7HDlLu7NVxVYuzKwhsLt/N6uqV2GtxWM8SgREpCDlsmrQw9ZaY61db63dmL3dl6vjieTS3p4RAFbUag0BmX5vfkkjPo/hji2aNDxb7euJArC8tpjuaDf9sX6WVy0n7aYJeAOaWyQiBUm1zESmoK03SkWRn8pijfOV6VdbGuSq9fX84IkDWml4lmrrjeL1GBZVBmkfaAdgWeUy0m6aIn9RfoMTETkFJQIiU9DWG6W5pjjfYcgcdsPFy4gmHW7fsj/fochZaOuNsrgyjDEubUNtQCYRcFyHkC+U5+hERE5OiYDIFCgRkFxbs6icV55Tw61/bCeeUinR2WZf9jsi7aZpH2zHa7wsKV8CaEVhESlcSgRETmM0mebwUJxlSgQkxz58yXJ6RxL89OnOfIciZ8BaS3tvlOaaEpJOkgNDB1hctpiAN4BFFYNEpHApERA5jfbeUQCaazRRWHLrwuXVrGso55bf78Nxbb7DkSk6MpwglnJori0mkU7QPtjOssplWGvxerxKBESkYCkREDmNtt5MNRANDZJcM8bwvy5ZTltvlPuf78p3ODJF+3ozVcWW1RQTS8VoG2yjubKZlJsi5NX8ABEpXEoERE6jLfsj31Sjyh+Se1esrWNJVRH/+eBerFWvwGwwdrGgqbqIzkgno6nR8YnCYX84z9GJiJyaEgGR02jrHaWuLERRQN37kntej+Ejly5nW8cQv37uSL7DkSlo64kS8nuoKfHRNjihYpB1CPqCeY5OROTUlAiInEZb74iGBcmMestLG1mxoIR/+tVOUo6b73DkNNp6ozRVF+Pi0D7YDsDyyuUA+D2qGCQihUuJgMhptPVGaa5VIiAzx+f18MkrVrGvN8oPntBqw4WurTfKstpM6dD9g/sJeUPUl9aDRROFRaSgKREQmcRANMnAaEqlQ2XGvXr1As5vquRrv91NNJHOdzhyCinH5UD/KE3VxSSdJPuH9tNU0QRkkgCvx5vfAEVEJqFEQGQSbX2qGCT5YYzhU1eupnckwX//oS3f4cgpdAzESLuW5ppMIjBWOjTlpDQ/QEQKnhIBkUm09SgRkPx56dJKrlhTxy2/30t3JA5A2k3TMdRBIp3Ic3QCR6uKLastZjQ1yoGhAyyrXIZrXUI+lQ4VkcKmREBkEm29Ubwew+IqlQ6V/PjklatIOZa/v+d5ILOKbSQZYf/Qfobjw3mOTvaNXywooX2wnZSbGi8dqh4BESl0SgREJtHWG2VJVRF+r/6rSH401xTz0ctXcO/2w/z2+Uw5Ua/HS5G/iEMjhzgycgTXqrJQvrT1RikP+ykLedg3uA/IlA4FTRQWkcKnsxuRSWTKAqo3QPLrQ5csp2VhKZ/++bNE4pmJwx7joSxYxnBimAODB0g6yTxHOT+19UZprinGxaV9oB1QIiAis4cSAZFTsNZmf+RL8h2KzHMBn4d/fPM6uobjfPU3u4/ZVhzInITuH9xPNBnNU4TzV1tvlGU1xTiuw/6h/ZQHy6kKV2GMwWtUMUhECpsSAZFTODKcIJZytIaAFITzllRy/YVL+c8n/5P7ntt+zLaQL0TQF+TA0AF6o71Ya/MU5fwSSzocHorTXJNZQ6B9sJ3mimYc6+Dz+DDG5DtEEZFJ5SwRMMb8jzGm2xjzbK6OIZJL+8aqgahikBSI91+8gIj/B9z0wCfoGj52orDP46MsWEZfrI+O4Q5STipPUc4f7WPlhWuPLR3quI4qBonIrJDLHoFvA1fk8P1FcqqtV6VDpbAsrljAv7z2m8TtQd56x1+dcOXfGENpsHT8pFRDhXJr4nfEUHyIQ5FDmUTAKhEQkdkhZ4mAtfb3QH+u3l8k19p6ooT8HurK9IMuheOGzddweeN72Tf6S/73L/7npG3C/jBBX5CDQwc5MnIEx3VmOMr5YSwRaKouZnf/bix2fA2BgDeQ5+hERE4v73MEjDE3GGO2GmO29vT05DsckXHtfaMsrSrG49E4XyksX7/mr1gY3MCdu77AL1/YftI2Po+PslAZkUSEtoE29Q7kwP6+KDUlQYqDPtoGM6s/LylfAqhikIjMDnlPBKy1t1hrN1lrN9XW1uY7HJFxHQOjLK4K5zsMkRP4vX5++PZb8JoQN/7yQ+zt6T1l26JAEQFfgINDB+kc7lSZ0WnUMRBjcVUYx3XoGO4AjiYCXo8qBolI4ct7IiBSiKy1dA7EaKzUGgJSmFbUNPJPr/o34raTa+64nv5o7JRtx3oH4uk4bQNt9EZ7NVxoGnQMxFhcWYRjHTqHOwn7wlSGKjEY9QiIyKygREDkJIZjaSKJNI2V6hGQwvX2Da/iY5s+z6DzDK//9odJpic/uQ/7w5QESuiP9dM20MZAbECrEp8lx7UcGozRWJnpEeiMdLKkfAkumh8gIrNHLsuHfh94FGgxxnQYY96fq2OJTLeDA6MA6hGQguP1ePHgGT+B/5tL3s0bl9/IwcT9/Mltn8Y9zRoCxhhKgiWE/CF6RnvY17+PofiQegjOUNdwnLRraawsIu2m6RzupLGskbSbVsUgEZk1clk16Dprbb211m+tbbTWfitXxxKZbh3jiYB6BKSweIyH6qJqRlOj46/9+zWf4vwFb+SZwe/wlts+T9o9/YJiHuOhJFBC0BfkSPQI+wb20Rvt1foDU9TRf/Q7Iukk6Yx0srhssdYQEJFZRUODRE6iYyAz3nqxegSkAJWHyvHgGb+Kb4zhR9d9jQ1VV7Kl7xtce+unSaWnNuTH6/FSEigh7A8zEB9g38A+OoY6iCajGjY0ifHviKoiekd7GUoMsbh8MdZa/F5/nqMTEZkaJQIiJ9ExEKM06KMsrAl/Ung8xkNtcS2x9NEJwn6vn3ve/V+cv+Banhm6lSv/52+Ip6Y+3MdjPBQHiikNlpJyU3QMd7Cvfx/d0W5iqdgJi5fNd2PDBxdVhNjbvxeAxrLG/9fenUdJVpZ3HP8+VbequnqdhYEZZnoGZREQ2QKMiAsocoyiBFcQEjgSlkQFEjEnJmqiOSjHaNwAIyKJGgOyRjQg24GIoCwizgyLyDbTDQ0zPd3TW3VX1/Lkj3u7qenpWeiZ2n+fc/pU1b233nr63qq696l3AyBuGjFIROqDEgGRWfQOZlg6P42Z5hCQ2tSR6iBu8c3a9sdjcW487VLevORDPDF6NcdecS4vDY2+6rJTQYqOVActiRZGsiOsG1rHM4PPsH5sPZlcRv0JCH8s2KMzRSqIs3ZoLQDdnd2A5hAQkfqhREBkFuH44GoWJLUrZjEWtS7arK/A1PJrTvk6J+19Hj3ZW3nb90/m/mfXzfk10ol0mBQEYVLQM9TDMwPPsG5oHYPjg4znxpuyCVHvYIZl81vDOQRGwjkElnYuJRaLaQ4BEakbSgREZnB3egYy6igsNa891U4iltiig6+ZcflJn+Pzx3yTDM9w6jQyX1YAABQ9SURBVE3v5bv337tTr1WaFLSn2nF3+jP99Az18PTGp3l+8Hk2jG1gdHKUbD7b8MlB72A0dGg0h0BropXOZCepeKraoYmI7DDVX4rMsCmTY2yyoKFDpebFLMaSjiWsHVpLEAu2aMp27soPctDivTnjpo/xxV+fzi+eOoerPnQR89t2/mI1EU9s1ik2X8wznB1mcGIQHBwnGU/SErTQErSQjCcJYgFBLKj7X8zzhSJ9QxN0lwwd2t3ZTcELtAVt1Q5PRGSHKREQmWFqNBDVCEg9SCfSLEwvZNPEJtqSW16EHrPiMO4/+05Ov/ZCHhy4nDd+916+dsLXOfGgA3ZpHFMX+aUKxQLj+XFGJ0enawjMDMNIxBKkghSJWJhQTCUIMYsRt/C2Vvvo9A1NUCj6ZpOJLetcRsELqhEQkbqiREBkBs0hIPVmQXoBI9kRcoXcrENX7t62kNvO/CHfuv/HfO2BL3Debe/lqIc/xqUnnc+e89vLFlc8Fp/11393p+AFJvITZDxDwaPOx1MDE0XX/1NJwVSSMH3fwnLNLEwYsOkEYyqBKF0G7NKk4pUfC1qn5xBYuXQlOARxnVZFpH7oG0tkBs0qLPUmHouHTYQ2zd5ECMIL4QuOOZ0T938LZ/3P3/HAwGUcc9VPOf2Az/DZE95HKqhclzEzI7AtaxBmcneKXqToRfL5PI5PL/OprMEBC7c1LEwiomWb3QKxWIxY1DVuOokoSSCALW5L708lF0++vB6AjtZJXhjpZzg7HM4hgGvoUBGpK0oERGboHRynsyWgK61JgaR+pBNpFrYuZHBikPbk1n/l33vhCu456yf8ZNWt/NM9/8xVT36c6//wfT528Plc8LZ3kKxgQrA9Zhb++s+uubh29+lkApi+7/h0rUTputLnlXqufwgDkskx/jjwPBAOHWpmGjpUROpK7Xzji9SIcDQQ1QZI/VnYupB0kN5iSNHZfOTgP+X3f/1Lzjzo00zwLN/4/Zkc8u2TuPiOWxgZz1cg2sqbqgWYarIUxILpTs/JeJJkPEkqSJEKUtOdnFuCFtKJ9GZ/G0acRR0Julra6M/0A9Dd1b1ZzYKISD3QN5bIDOH44OofIPVnahShGDGy+ex2t08FKS4+4ULWfPzB6YTg8tVnc8jlx3LG1d/mib6BCkRdf/qGJlkyLwlAz3APAHu277nVZlkiIrVKiYBICXdXjYDUtSAWsKxrGblCjnxxx37Zb0u2cfEJF/L4Jx7mU0d+ifZUwJ19l/DO/17JykvP5st33MLG0e0nFs2ib1OWxV1RIjDUQ1uijfZkO8l4ssqRiYi8OkoEREoMjE2SmSyoRkDqWjKeZFnXMjK5DIViYYefl06k+du3nMHvP343P3jfdRy2+3H0Td7NpavP5rDvHsVbv3M+F9/+c3oGx8oYfW3LF5wNIzn2nBcOE9o73Et3Z9hRWEOHiki9Ua8mkRJTwwJ2L1CNgNS31kQryzqW0TvSS2ui9VV1YjUzjt/nTRy/z5vI5DL86JGfc/XqG3lm+GdcvuYGvrO6jYXB4Ry6+9G8+3Vv490HHERHujlOJy8PT1J0WNL1StOgZV3RHAKBEgERqS/N8c0tsoM0mZg0kvZUO8tjy+kZ6iGdSM9pRJvWRCvnrvww5678MOO5cW587C6uX3Mra/p/zZ1993Jn31f49N3zmRccyN5dh3Dk0sN4+96HcVj3HqQSjVfp3Dc0CcDiqI9A73Dv9BwC9T5jsog0HyUCIiWmJhNbqkRAGkRropXurm56h3unR8aZq3QizWmHnshph56Iu/P0xue4fs3d3Lv2Nzw7tIoHB+7jwQG4bDUExT2Yn9yX7o592Hv+Phy0x34cuWw/9ttjAelk/V4wv7QpTASWdCUZmhhiODvMss5lmkNAROpSWRMBM3sX8E0gDlzp7peU8/VEdlbv4Dhd6QSdLZpDQBpHa6KV5V3LeWH4BcaL46QTO5/omhn77vZaPnPsa/kMZwEwMD7APc8+zK+ef5Q1Lz/GupE/8Mjg/TwyWOS6Z8PnxXweaVvCgtRSFqYXs7h9T7q7lrJi/hL2Wbgn++62mEUdLcRjtTn6Tt/QJDGDPTqTPLnxj4DmEBCR+lW2by0ziwOXAe8EeoGHzOxmd3+8XK8psrN6NHSoNKiWoIUV81bQN9LH6OQobYm2XT7U5YL0At7/+hN4/+tPmF6WzWf5Q/+zPNTzJKtffornBtfRN9pDf3YNvdm7eXRTITxDTHEjRgcJ6yIVm0dbMJ+2RBedyS66WrroSnXR1dLJvJYO5rV0sqC1k/npdha0dbCgtZ0Fra20pQJiZRrGs28oy6KOBEHc6BkKhw7VHAIiUq/K+fPFUcDT7v4sgJldA5wE1GQiMJErsG5g+5PwSGN7vn+M/Rd3VjsMkbIIYgHLOpfRn+ln4/hG0kGaRLy8tV+pIMXBiw/g4MUHbLGu6EVeHl3PHzb08PTGF1k7+BIvDr/Ehkw/Q9kBRiY3MZpfy8bcCLmxYZwdGw7VPIWRIk6KmCWnb4NYkrglCCxJEEsQxKZuAxKxBEEsMT3RWBBLEFiceCwgEQuIx2IEsYA1L4zTnkryX6se5Te9vwE0h4CI1K9yJgJLgZ6Sx73AyjK+3k556uUR3nfpfdUOQ2rAuw5aUu0QRMrGzFjUtoj2ZDt9I31MFiZpTbRW5SI2nABtMUs6FnPsa7e9rbuTyWXoz2xiw+gwG8aGGMgMsWlijE3jIwxPjDI6mWEsl2FsMsNEfpzJwgSTxSyThSz5Yo5cMUu+OMGED1PM5yh4+Ff0PE6eInmcAu4FsG0kHTlYdWd4d3H7Ys0hICJ1q5yJwGxnFd9iI7NzgHMAli9fXsZwtm35glYu++jhVXt9qQ0xgzfts1u1wxApu3QizYp5KxgYH2BjZiPJIElL0FLtsLbKzGhLttGWbGPFvKUVec2iF8kX80wWcmRzeXKFPJOFAqmEU/QiRS8yr2Uejtf0vhMR2ZpyJgK9QHfJ42XAizM3cvcrgCsAjjjiiC0ShUqZ15rkPQfrl2ARaR7xWJxFbYvoTHWyfmw9wxPDpBPlby5UL2IWmx5pqX0bP/iPTo6qRkBE6lI5ezY9BOxrZq8xsyRwCnBzGV9PRETmIBWk6O7qZvm85bg7w9lhsvlstcOqH5pDQETqVNlqBNw9b2afAG4jHD70Knd/rFyvJyIiO6c10cqKeSsYz4/Tn+lnJDtCEAtoCVrUEXYbNIeAiNSrsg567O63ALeU8zVERGTXMbPpeQcm8hMMZ4cZmhii6EVSQUpNYGahOQREpF7pm0tERGbVErTQErSwML2QTC7DpolNjGRHMDMSsQTJeLLpawrcHcPUNEhE6pISARER2aZ4LE5HqoOOVAf5Yp6J/ARDE0OM5cZw9+lOtc04ln7Ri6oNEJG6pW8vERHZYUEsoD3ZTnuynaIXmSxMMp4bZ2RyhLHcWLhR1Hk2EU80/EVyvpgnFaSqHYaIyJw09je0iIiUTcxi082H5qfnTycGuUKOTC5DJpdhIj8x3XzGLGxCE7d4w9QeFL1IKq5EQETqkxIBERHZJUoTg45UB/DKpFz5Yp5cIcdEfoJsIct4fpxisRhOPemAgWHELLbFXy0nDAUvqAO1iNQtJQIiIlI2pZNykYAuuqbXFb1IoVig4AUKxcJ0wjBZmJy+X/ACRS9OP2eqdsFxzGy6j4KZYdGE9mZhQgFM10SU3k5ts0toDgERqWNKBEREpCpiFiMWj5Fg2zMZuztFL07/OeHjmcsLxQJFotvSbYvhrbtTpBjWRABM5QLOZolFaaIxtT688c0TCNccAiJS35QIiIhITTMz4hYnzq694HZ3PLrKd49uZzze3jJATYNEpG4pERARkaZU2lSI2u2GICJSNrFqByAiIiIiIpWnREBEREREpAkpERARERERaUJKBEREREREmpASARERERGRJqREQERERESkCSkREBERERFpQjUzj4CZnQP0m9na7WzaBQyVMZTdgP4ylV3O2Bu97Lkel1qIvV7Lnss+L/fnsxb2Sy2UXXps6inuSpZfqbJ39TmjEfZJrZS/o8emXvdLPZU981jUU+z1VHbpfl6xQ89w95r4Ax7ewe2uqIU45lh22WJv9LLnelxqIfZ6LXsu+7wCn8+q75daKLv02NRT3PW8z7dW9q4+ZzTCPqmV8mvhukJlz34s6in2eip7Lt9H9dg06GfVDmAnlDN2lV358lV2Zcsud/kqu7Jll7t8ld04ZZe7fJVd2bLLXb7KfhUsyiCqzswedvcjFIfMRsel8rTPa5eOTe3QsahdOja1Q8eiMuayn2upRuCKagcQqZU4ZHM6LpWnfV67dGxqh45F7dKxqR06FpXxqvdzzdQIiIiIiIhI5dRSjYCIiIiIiFSIEgERERERkSbUlImAmY1WOwbZnJkVzOzRkr+9trHtsWb288pF15jMzM3sRyWPAzPboH1bG8zs5OgY7V/tWJqVPiP1Qef02rO9Y2Jm95iZOg/XgKZMBKQmjbv7oSV/z1c7oCYwBhxkZuno8TuBF15NAWZWM5MSNqBTgV8Bp7yaJ5lZvDzhNKWd/oyIiNSypk0EzKzdzO4ys0fMbLWZnRQt38vMnjCz75nZY2Z2e8lJQCrIzOJm9q9m9pCZrTKzc0tWd5rZTWb2uJn9u5k17Xt5J90KvCe6fypw9dQKMzvKzO43s99Ft6+Llp9pZteZ2c+A2ysfcuMzs3bgGOAsokQgqgn75WzvezMbNbMvmtkDwNHVi7whzeUzcq+ZHVqy3X1mdnBFo24yM2uKzexSMzszuv+8mX2h5HyvWrYK2NYxkdrRzBdPE8DJ7n44cBzwNTOzaN2+wGXu/npgE/CBKsXYTNIlzYJuipadBQy5+5HAkcDZZvaaaN1RwKeANwB7A++veMSN4RrgFDNrAQ4GHihZ9yTwVnc/DPg88KWSdUcDZ7j72ysWaXP5M+AX7v4UMGBmh0fLt/a+bwPWuPtKd/9VxaNtbHP5jFwJnAlgZvsBKXdfVbGIZTb90fn+O8BF1Q5GpFY0c7W+AV8ys7cCRWApsEe07jl3fzS6/1tgr8qH13TG3f3QGctOAA42sw9Gj7sIk7RJ4EF3fxbAzK4G3gxcX6lgG4W7r4r6Y5wK3DJjdRfwAzPbF3AgUbLuDncfqEiQzelU4BvR/Wuix//L1t/3BeCGKsTZ8Ob4GbkO+JyZfRr4GPCfFQlWtuXG6Pa36IcjkWnNnAicBiwC/sTdc2b2PNASrcuWbFcA1DSoOgz4pLvfttlCs2MJT7qlNCHG3N0MfBU4FlhYsvxfgLvd/eToQuieknVjFYqt6ZjZQuDthG3THYgTvr9vYevv+wl3L1Quyqbzqj4j7p4xszuAk4APA+oUWX55Nm/l0DJj/dR5vUBzX/tU0vaOidSAZm4a1AWsj5KA44AV1Q5ItnAb8FdmloCwit3M2qJ1R5nZa6I20h8h7FQpc3MV8EV3Xz1jeRevdIw8s6IRNbcPAj909xXuvpe7dwPPEf76r/d9dczlM3Il8C3gIdWeVcRa4EAzS5lZF/COagckOib1oOkSgWiUkyzwY+AIM3uYsHbgyaoGJrO5EngceMTM1gDf5ZVfcn4NXAKsIbxIumnWEmS73L3X3b85y6qvAF82s/sIf5WWyjiVLd/PNwAfRe/7qpjLZ8TdfwsMA/9RgRCb1tQ53d17gGuBVYTn999VNbAmpmNSX8y9uVpUmNkhwPfc/ahqxyIisqOiJnEXufuJ1Y5Fts/M9iRsKrS/uxerHE7D0jm99uiY1JemqhEws/MIh377bLVjERGRxmRmf0E4utA/KgkoH53Ta4+OSf1puhoBERERERFpghoBM+s2s7ujScIeM7MLouULzOwOM/tjdDs/Wr6/mf3azLJmdtGMsv4mKmONmV0djSstIiIiIlJ3Gj4RIBy+6lPufgDwRuDjZnYg8PfAXe6+L3BX9BhgADifcKi4aWa2NFp+hLsfRNgx7JTK/AsiIiIiIrtWwycC7t7n7o9E90eAJwgnDzsJ+EG02Q8IZ/LE3de7+0NAbpbiAsIZcAOgFXixzOGLiIiIiJRFwycCpaIJXw4j7MS1h7v3QZgsALtv67nu/gJhLcE6oA8YcvfbyxmviIiIiEi5NE0iYGbthGNxX+juw3N4/nzCWoTXAHsCbWZ2+q6NUkRERESkMpoiEYhmpr0B+LG73xgtftnMlkTrlwDrt1PM8cBz7r7B3XPAjcCbyhWziIiIiEg5NXwiYGYGfB94wt3/rWTVzcAZ0f0zgJ9up6h1wBvNrDUq8x2E/Q1EREREROpOw88jYGZvBu4FVgNTE7v8A2E/gWuB5YQX+R9y9wEzWww8DHRG248CB7r7sJl9AfgI4UhEvwP+0t2zlfx/RERERER2hYZPBEREREREZEsN3zRIRERERES2pERARERERKQJKREQEREREWlCSgRERERERJqQEgERERERkSYUVDsAERGpHjNbCNwVPVwMFIAN0eOMu2viRBGRBqXhQ0VEBAAz+2dg1N2/Wu1YRESk/NQ0SEREZmVmo9HtsWb2f2Z2rZk9ZWaXmNlpZvagma02s72j7RaZ2Q1m9lD0d0x1/wMREdkWJQIiIrIjDgEuAN4A/Dmwn7sfBVwJfDLa5pvA1939SOAD0ToREalR6iMgIiI74iF37wMws2eA26Plq4HjovvHAwea2dRzOs2sw91HKhqpiIjsECUCIiKyI7Il94slj4u8ci6JAUe7+3glAxMRkblR0yAREdlVbgc+MfXAzA6tYiwiIrIdSgRERGRXOR84wsxWmdnjwHnVDkhERLZOw4eKiIiIiDQh1QiIiIiIiDQhJQIiIiIiIk1IiYCIiIiISBNSIiAiIiIi0oSUCIiIiIiINCElAiIiIiIiTUiJgIiIiIhIE1IiICIiIiLShP4fmH7fnV9I/RsAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 936x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"mod = AR1X(endog, exog=exog)\n",
"res = mod.fit(maxiter=1000)\n",
"print(res.summary())\n",
"\n",
"startPredDate = '2018-04-01'\n",
"predict_dy = res.get_prediction(dynamic=startPredDate)\n",
"predict_dy_ci = predict_dy.conf_int()\n",
"\n",
"# Graph\n",
"fig, ax = plt.subplots(figsize=(13, 4))\n",
"ax.set(title='Modeled response', xlabel='Time', ylabel='output')\n",
"# Plot data points\n",
"endog.plot(ax=ax, label='Observed')\n",
"\n",
"# Plot predictions\n",
"predict_dy.predicted_mean.loc[startPredDate:].plot(ax=ax, style='g', label='Dynamic forecast')\n",
"ci = predict_dy_ci.loc[startPredDate:]\n",
"ax.fill_between(ci.index, ci.iloc[:,0], ci.iloc[:,1], color='g', alpha=0.1);"
]
}
],
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@bobbejaan
Copy link

I get the following error:
ValueError: Invalid state space initialization method.
I guess it is caused by:
super(AR1X, self).__init__(endog, k_states=1, exog=exog, initialization='diffuse')

The following does work:
super(AR1X, self).__init__(endog, k_states=1, exog=exog, initialization='approximate_diffuse')

@ChadFulton
Copy link
Author

Ah, yes, for diffuse you need the version of Statsmodels in master, but approximate_diffuse will work in v0.9

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment