{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# State space modeling: Local Linear Trends" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook describes how to extend the statsmodels statespace classes to create and estimate a custom model. Here we develop a local linear trend model.\n", "\n", "The Local Linear Trend model has the form (see Durbin and Koopman 2012, Chapter 3.2 for all notation and details):\n", "\n", "$$\n", "\\begin{align}\n", "y_t & = \\mu_t + \\varepsilon_t \\qquad & \\varepsilon_t \\sim\n", " N(0, \\sigma_\\varepsilon^2) \\\\\n", "\\mu_{t+1} & = \\mu_t + \\nu_t + \\xi_t & \\xi_t \\sim N(0, \\sigma_\\xi^2) \\\\\n", "\\nu_{t+1} & = \\nu_t + \\zeta_t & \\zeta_t \\sim N(0, \\sigma_\\zeta^2)\n", "\\end{align}\n", "$$\n", "\n", "It is easy to see that this can be cast into state space form as:\n", "\n", "$$\n", "\\begin{align}\n", "y_t & = \\begin{pmatrix} 1 & 0 \\end{pmatrix} \\begin{pmatrix} \\mu_t \\\\ \\nu_t \\end{pmatrix} + \\varepsilon_t \\\\\n", "\\begin{pmatrix} \\mu_{t+1} \\\\ \\nu_{t+1} \\end{pmatrix} & = \\begin{bmatrix} 1 & 1 \\\\ 0 & 1 \\end{bmatrix} \\begin{pmatrix} \\mu_t \\\\ \\nu_t \\end{pmatrix} + \\begin{pmatrix} \\xi_t \\\\ \\zeta_t \\end{pmatrix}\n", "\\end{align}\n", "$$\n", "\n", "Notice that much of the state space representation is composed of known values; in fact the only parts in which parameters to be estimated appear are in the variance / covariance matrices:\n", "\n", "$$\n", "\\begin{align}\n", "H_t & = \\begin{bmatrix} \\sigma_\\varepsilon^2 \\end{bmatrix} \\\\\n", "Q_t & = \\begin{bmatrix} \\sigma_\\xi^2 & 0 \\\\ 0 & \\sigma_\\zeta^2 \\end{bmatrix}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-11-02T17:02:05.295532Z", "iopub.status.busy": "2022-11-02T17:02:05.294794Z", "iopub.status.idle": "2022-11-02T17:02:08.539274Z", "shell.execute_reply": "2022-11-02T17:02:08.538504Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import numpy as np\n", "import pandas as pd\n", "from scipy.stats import norm\n", "import statsmodels.api as sm\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To take advantage of the existing infrastructure, including Kalman filtering and maximum likelihood estimation, we create a new class which extends from `statsmodels.tsa.statespace.MLEModel`. There are a number of things that must be specified:\n", "\n", "1. **k_states**, **k_posdef**: These two parameters must be provided to the base classes in initialization. The inform the statespace model about the size of, respectively, the state vector, above $\\begin{pmatrix} \\mu_t & \\nu_t \\end{pmatrix}'$, and the state error vector, above $\\begin{pmatrix} \\xi_t & \\zeta_t \\end{pmatrix}'$. Note that the dimension of the endogenous vector does not have to be specified, since it can be inferred from the `endog` array.\n", "2. **update**: The method `update`, with argument `params`, must be specified (it is used when `fit()` is called to calculate the MLE). It takes the parameters and fills them into the appropriate state space matrices. For example, below, the `params` vector contains variance parameters $\\begin{pmatrix} \\sigma_\\varepsilon^2 & \\sigma_\\xi^2 & \\sigma_\\zeta^2\\end{pmatrix}$, and the `update` method must place them in the observation and state covariance matrices. More generally, the parameter vector might be mapped into many different places in all of the statespace matrices.\n", "3. **statespace matrices**: by default, all state space matrices (`obs_intercept, design, obs_cov, state_intercept, transition, selection, state_cov`) are set to zeros. Values that are fixed (like the ones in the design and transition matrices here) can be set in initialization, whereas values that vary with the parameters should be set in the `update` method. Note that it is easy to forget to set the selection matrix, which is often just the identity matrix (as it is here), but not setting it will lead to a very different model (one where there is not a stochastic component to the transition equation).\n", "4. **start params**: start parameters must be set, even if it is just a vector of zeros, although often good start parameters can be found from the data. Maximum likelihood estimation by gradient methods (as employed here) can be sensitive to the starting parameters, so it is important to select good ones if possible. Here it does not matter too much (although as variances, they should't be set zero).\n", "5. **initialization**: in addition to defined state space matrices, all state space models must be initialized with the mean and variance for the initial distribution of the state vector. If the distribution is known, `initialize_known(initial_state, initial_state_cov)` can be called, or if the model is stationary (e.g. an ARMA model), `initialize_stationary` can be used. Otherwise, `initialize_approximate_diffuse` is a reasonable generic initialization (exact diffuse initialization is not yet available). Since the local linear trend model is not stationary (it is composed of random walks) and since the distribution is not generally known, we use `initialize_approximate_diffuse` below.\n", "\n", "The above are the minimum necessary for a successful model. There are also a number of things that do not have to be set, but which may be helpful or important for some applications:\n", "\n", "1. **transform / untransform**: when `fit` is called, the optimizer in the background will use gradient methods to select the parameters that maximize the likelihood function. By default it uses unbounded optimization, which means that it may select any parameter value. In many cases, that is not the desired behavior; variances, for example, cannot be negative. To get around this, the `transform` method takes the unconstrained vector of parameters provided by the optimizer and returns a constrained vector of parameters used in likelihood evaluation. `untransform` provides the reverse operation.\n", "2. **param_names**: this internal method can be used to set names for the estimated parameters so that e.g. the summary provides meaningful names. If not present, parameters are named `param0`, `param1`, etc." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-11-02T17:02:08.543632Z", "iopub.status.busy": "2022-11-02T17:02:08.543259Z", "iopub.status.idle": "2022-11-02T17:02:08.553224Z", "shell.execute_reply": "2022-11-02T17:02:08.552589Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "\"\"\"\n", "Univariate Local Linear Trend Model\n", "\"\"\"\n", "class LocalLinearTrend(sm.tsa.statespace.MLEModel):\n", " def __init__(self, endog):\n", " # Model order\n", " k_states = k_posdef = 2\n", "\n", " # Initialize the statespace\n", " super(LocalLinearTrend, self).__init__(\n", " endog, k_states=k_states, k_posdef=k_posdef,\n", " initialization='approximate_diffuse',\n", " loglikelihood_burn=k_states\n", " )\n", "\n", " # Initialize the matrices\n", " self.ssm['design'] = np.array([1, 0])\n", " self.ssm['transition'] = np.array([[1, 1],\n", " [0, 1]])\n", " self.ssm['selection'] = np.eye(k_states)\n", "\n", " # Cache some indices\n", " self._state_cov_idx = ('state_cov',) + np.diag_indices(k_posdef)\n", "\n", " @property\n", " def param_names(self):\n", " return ['sigma2.measurement', 'sigma2.level', 'sigma2.trend']\n", "\n", " @property\n", " def start_params(self):\n", " return [np.std(self.endog)]*3\n", "\n", " def transform_params(self, unconstrained):\n", " return unconstrained**2\n", "\n", " def untransform_params(self, constrained):\n", " return constrained**0.5\n", "\n", " def update(self, params, *args, **kwargs):\n", " params = super(LocalLinearTrend, self).update(params, *args, **kwargs)\n", " \n", " # Observation covariance\n", " self.ssm['obs_cov',0,0] = params[0]\n", "\n", " # State covariance\n", " self.ssm[self._state_cov_idx] = params[1:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using this simple model, we can estimate the parameters from a local linear trend model. The following example is from Commandeur and Koopman (2007), section 3.4., modeling motor vehicle fatalities in Finland." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-11-02T17:02:08.556395Z", "iopub.status.busy": "2022-11-02T17:02:08.556156Z", "iopub.status.idle": "2022-11-02T17:02:10.010649Z", "shell.execute_reply": "2022-11-02T17:02:10.009868Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import requests\n", "from io import BytesIO\n", "from zipfile import ZipFile\n", " \n", "# Download the dataset\n", "ck = requests.get('http://staff.feweb.vu.nl/koopman/projects/ckbook/OxCodeAll.zip').content\n", "zipped = ZipFile(BytesIO(ck))\n", "df = pd.read_table(\n", " BytesIO(zipped.read('OxCodeIntroStateSpaceBook/Chapter_2/NorwayFinland.txt')),\n", " skiprows=1, header=None, sep='\\s+', engine='python',\n", " names=['date','nf', 'ff']\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we defined the local linear trend model as extending from `MLEModel`, the `fit()` method is immediately available, just as in other statsmodels maximum likelihood classes. Similarly, the returned results class supports many of the same post-estimation results, like the `summary` method.\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-11-02T17:02:10.014866Z", "iopub.status.busy": "2022-11-02T17:02:10.014607Z", "iopub.status.idle": "2022-11-02T17:02:10.057368Z", "shell.execute_reply": "2022-11-02T17:02:10.056556Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Statespace Model Results \n", "==============================================================================\n", "Dep. Variable: lff No. Observations: 34\n", "Model: LocalLinearTrend Log Likelihood 27.510\n", "Date: Wed, 02 Nov 2022 AIC -49.020\n", "Time: 17:02:10 BIC -44.623\n", "Sample: 01-01-1970 HQIC -47.563\n", " - 01-01-2003 \n", "Covariance Type: opg \n", "======================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "--------------------------------------------------------------------------------------\n", "sigma2.measurement 0.0010 0.003 0.346 0.730 -0.005 0.007\n", "sigma2.level 0.0074 0.005 1.564 0.118 -0.002 0.017\n", "sigma2.trend 2.501e-11 0.000 1.67e-07 1.000 -0.000 0.000\n", "===================================================================================\n", "Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 0.68\n", "Prob(Q): 0.95 Prob(JB): 0.71\n", "Heteroskedasticity (H): 0.75 Skew: -0.02\n", "Prob(H) (two-sided): 0.64 Kurtosis: 2.29\n", "===================================================================================\n", "\n", "Warnings:\n", "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n" ] } ], "source": [ "# Load Dataset\n", "df.index = pd.date_range(start='%d-01-01' % df.date[0], end='%d-01-01' % df.iloc[-1, 0], freq='AS')\n", "\n", "# Log transform\n", "df['lff'] = np.log(df['ff'])\n", "\n", "# Setup the model\n", "mod = LocalLinearTrend(df['lff'])\n", "\n", "# Fit it using MLE (recall that we are fitting the three variance parameters)\n", "res = mod.fit(disp=False)\n", "print(res.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can do post-estimation prediction and forecasting. Notice that the end period can be specified as a date." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-11-02T17:02:10.061466Z", "iopub.status.busy": "2022-11-02T17:02:10.061198Z", "iopub.status.idle": "2022-11-02T17:02:10.068441Z", "shell.execute_reply": "2022-11-02T17:02:10.067776Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Perform prediction and forecasting\n", "predict = res.get_prediction()\n", "forecast = res.get_forecast('2014')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-11-02T17:02:10.071581Z", "iopub.status.busy": "2022-11-02T17:02:10.071334Z", "iopub.status.idle": "2022-11-02T17:02:10.377653Z", "shell.execute_reply": "2022-11-02T17:02:10.376738Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAFlCAYAAAAj9p2/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAACIx0lEQVR4nO3dd3gUVdsG8Hu2Z9MT0hshofcmBlRA0YBY0FdRRBFBLMArWBD5bGDDhhUbFrCgWBD0VRQpUgREQEBASkghPUBCerbOfH9sdsmStpvspmzu33UtIbOzM2dONjDPPuc8R5AkSQIREREREZEHkbV2A4iIiIiIiFyNgQ4REREREXkcBjpERERERORxGOgQEREREZHHYaBDREREREQeh4EOERERERF5HAY6RERERETkcRjoEBERERGRx2GgQ0REREREHoeBDhEREREReRynAh2z2Ywnn3wS8fHx8PLyQkJCAp599llIktTg67Zs2YJBgwZBrVYjMTERK1asaE6biYiIiIiIGqRwZueXXnoJ7733Hj799FP07t0be/fuxV133QV/f3888MADdb4mPT0d48ePx3333YeVK1di06ZNuPvuuxEREYHk5GSXXAQREREREVFNgtRYOqaGa665BmFhYfj4449t2/7zn//Ay8sLX3zxRZ2vmT9/Pn7++WccPnzYtu3WW29FcXExfv3112Y0nYiIiIiIqG5OZXSGDx+OZcuW4cSJE+jWrRsOHjyIP/74A6+99lq9r9m1axfGjBljty05ORlz586t9zV6vR56vd72vSiKKCoqQnBwMARBcKbJRERERETkQSRJQllZGSIjIyGT1T8Tx6lA57HHHkNpaSl69OgBuVwOs9mM559/HpMnT673Nfn5+QgLC7PbFhYWhtLSUlRVVcHLy6vWaxYvXoxFixY50zQiIiIiIupAsrKyEB0dXe/zTgU633zzDVauXIkvv/wSvXv3xoEDBzB37lxERkbizjvvbHZjrRYsWICHHnrI9n1JSQliY2ORlZUFPz8/l52HiIiIiIjal9LSUsTExMDX17fB/ZwKdObNm4fHHnsMt956KwCgb9++OHXqFBYvXlxvoBMeHo6CggK7bQUFBfDz86szmwMAarUaarW61nY/Pz8GOkRERERE1OiUFqfKS1dWVtYaByeXyyGKYr2vSUpKwqZNm+y2bdiwAUlJSc6cmoiIiIiIyGFOBTrXXnstnn/+efz888/IyMjAmjVr8Nprr+GGG26w7bNgwQJMmTLF9v19992HtLQ0PProozh27BjeffddfPPNN3jwwQdddxVEREREREQ1ODV07e2338aTTz6JmTNn4vTp04iMjMS9996Lp556yrZPXl4eMjMzbd/Hx8fj559/xoMPPog333wT0dHR+Oijj7iGDhERERERuY1T6+i0ltLSUvj7+6OkpIRzdIiIiIiIOjBHYwOnhq4RERERERG1Bwx0iIiIiIjI4zDQISIiIiIij8NAh4iIiIiIPA4DHSIiIiIi8jgMdIiIiIiIyOMw0CEiIiIiIo/DQIeIiIiIiDwOAx0iIiIiIvI4DHSIiIiIiMjjMNAhIiIiIiKPw0CHiIiIiIg8DgMdIiIiIiLyOAx0iIiIiIjI4zDQISIiIiIij8NAh4iIiIiIPA4DHSIiIiIi8jgMdIiIiIiIyOMw0CEiIiIiIo/DQIeIiIiIiDwOAx0iIiIiIvI4DHSIiIiIiMjjMNAhIiIiIiKPw0CHiIiIiIg8DgMdIiIiIiLyOAx0iIiIiIjI4zDQISIiIiIij8NAh4iIiIiIPI5TgU7nzp0hCEKtx6xZs+rcf8WKFbX21Wg0Lmk4ERERERFRfRTO7Lxnzx6YzWbb94cPH8aVV16Jm2++ud7X+Pn54fjx47bvBUFoQjOJiIiIiIgc51SgExISYvf9iy++iISEBIwcObLe1wiCgPDw8Ka1joiIiIiIqAmaPEfHYDDgiy++wLRp0xrM0pSXlyMuLg4xMTG4/vrrceTIkUaPrdfrUVpaavcgIiIiIiJyVJMDnbVr16K4uBhTp06td5/u3bvjk08+wQ8//IAvvvgCoihi+PDhyM7ObvDYixcvhr+/v+0RExPT1GYSEREREVEHJEiSJDXlhcnJyVCpVPjf//7n8GuMRiN69uyJSZMm4dlnn613P71eD71eb/u+tLQUMTExKCkpgZ+fX1OaS0REREREHqC0tBT+/v6NxgZOzdGxOnXqFDZu3Ijvv//eqdcplUoMHDgQJ0+ebHA/tVoNtVrdlKYRERERERE1beja8uXLERoaivHjxzv1OrPZjEOHDiEiIqIppyUiIiIiInKI04GOKIpYvnw57rzzTigU9gmhKVOmYMGCBbbvn3nmGfz2229IS0vD33//jdtvvx2nTp3C3Xff3fyWExERERER1cPpoWsbN25EZmYmpk2bVuu5zMxMyGTnY6dz585hxowZyM/PR2BgIAYPHoydO3eiV69ezWs1ERERERFRA5pcjKAlOTrhiIiIiIiIPJujsUGTy0t3VJIk4VyFARV6E0xmsbWbQ0REREREdWhS1bWOzmAWYTCLKNcDMkGAWimDSi6DWiFrcPFUIiIiIiJqGQx0mkmUJFQZzKiCGQIApVxmC3wUcibMiIiIiIhaAwMdF5JwPtsDAHKZAJXCkulRyZntISIiIiJqKQx03MgsVmd7DJZsj0ohqw585JDLGPQQEREREbkLA50WIgHQm0ToTSLKYIJCJsBbrYBGKW/tphEREREReRxOImklJlFCSZURZ8v1qDKY0Q6qfBMRERERtRsMdFqZWZRQqjPibLkBlQYTAx4iIiIiIhdgoNNGiJKEMp0JZ8r1qNAz4CEiIiIiag7O0WljJAko15tQYTBBq1JAq5RDxsIFREREREROYaDTRkkSUKE3oVJvgpdKDm+VggEPEREREZGDGOi0cRKAyuoS1ZrqgIelqYmIiIiIGsZAp52QAFQZzNAZzFAr5fBWyaGQc4oVEREREVFdGOi0MxIAndEMndEMjUIOHw0zPEREREREF2JKoB3TmcworjSwQhsRERER0QUY6LRzJlFCaZWptZtBRERERNSmMNDxADqTGRV6BjtERERERFYMdDxEud4EndHc2s0gIiIiImoTGOh4kNIqI0xmsbWbQURERETU6hjoeBAJQHGVEaLI4gRERERE1LEx0HGSwSSiuNIAsY1WOjOLEkp1xtZuBhERERFRq+I6Ok748WAunlx7GCVVRsgEIECrQpBWhUBvJQK1KgR5qxDoXcc2rQpeKnmLtVNvElGmM8JXo2yxcxIRERERtSUMdBxQUmXE0z8cxtoDubZtogQUVRhQVGEAzjR+DC+l3C74CfJW4dKunTA8oZNbFvysNJihlMugUbZcgEVERERE1FYw0GnErtRCPPzNAeSW6CATgFmjE3HT4GiU6022QOdcpQHnKow4V1n393qTiCqjGVXFZuQW62zH/uFALsL9NLhhYBSu7R+BYB+1S9teWmWEXCZAKecIRSIiIiLqWARJaqOTTWooLS2Fv78/SkpK4Ofn1yLn1JvMeO23E1i2PQ2SBMQGafH6LQMQKivHXwf/RZeEBERGRTd6HEmSUGU041yFEUXWQKjCgPSzFVh3OM+22KdCJmBU9xDcNDgaA2ICIAiuyfLIBAHB3irI3JA1IiIiIiJqaY7GBh0m0MnOzkZKSgq6du2K6OiGA5Tj+WWY+/UBHM0rBQDcMiQGT17bC19/8SnuueceiKIImUyGV998B7dNmdqk9gCAzmjG5mOnsfrvbBzOKbVt79LJGzcOisK4PhHw0TQ/6aaUyxCoVboseCIiIiIiai0MdGr4+OOP7QKUZcuWYfr06bX2E0UJy3dm4KVfj8FgEhHkrcLiG/siuXc4srOzERcXB1E8v06NXC7HnkPHHMrsNOZEQRlW78vGr0fyoTNazuGllCO5dxj+Mzga3cJ8m3V8L5UcfixOQERERETtnKOxgVOTNzp37gxBEGo9Zs2aVe9rvv32W/To0QMajQZ9+/bFunXrnDlls2VnZ9uCHAAQRRH33nsvsrOz7fbLK6nCHZ/sxrM//QuDScTo7iH4de6lSO4dDgBISUmxC3IAwGw2Iz0tzSXt7BbmiwVX98TP/70Uj1zVDfGdvFFlNGPtgVzc8fFfuPvTvVh3KA96k7lJx68ymFFlaNpriYiIiIjaG6fGRe3Zswdm8/mb5cOHD+PKK6/EzTffXOf+O3fuxKRJk7B48WJcc801+PLLLzFhwgT8/fff6NOnT/Na7qD6ApSTJ0/ahrD99E8uHl9jKRutUcrw+PheuH1YrN1Qr65du0Imk9XK6MR36eLS9vpoFLh5SAxuGhyN/ZnFWP13Nn4/fgaHckpwKKcEb2xMwbX9I3DDwChEB2qdOnaZzlKcQKVgcQIiIiIi8mzNGro2d+5c/PTTT0hJSalz/sctt9yCiooK/PTTT7ZtF198MQYMGID333/f4fM0Z+hafUPOMjIy4NcpDE//cARr9ucAAPpF++P1WwYgIcSnzmN9/PHHuPfee2E2myGXy/HKG0ubPEcnNycbaampDhU1KCzX48eDuVizPwcFpXrb9ou7BGHSRbG4uEuww+cVBCDYW+2WktZERERERO7mlqFrNRkMBnzxxReYNm1avZPcd+3ahTFjxthtS05Oxq5duxo8tl6vR2lpqd2jqaKjo7Fs2TLI5Zb1ZORyOT744APkGLww7o3tWLM/BzIB+O/liVh9//B6gxwAmD59OtLT07H6p/XYc+hYk4OcLz9bgSF9uuOma8diSJ/u+PKzFQ3uH+yjxl0j4rFm5gi8clM/JHUJhgDgz7QizFl1AP87mNvg62uSJKC40oB2MDWLiIiIiKjJmhzorF27FsXFxZg6dWq9++Tn5yMsLMxuW1hYGPLz8xs89uLFi+Hv7297xMTEAABycnKa1Nbp06cjIyMDv//+O06kpuF05HDc+uGfyCmuQmyQFt/el4SHr+ru0Hoz0dHRGHHpZU0uQJCbk41H5syymzM0b+5s5OZkN/JKQC4TcFm3ELxx6wCsvn84rukXAQB48Zdj2J95zuE2mETJVtaaiIiIiMgTNTnQ+fjjjzFu3DhERka6sj0AgAULFqCkpMT2yMrKAgD07t0bH3/8cZOOGR0djciegzF77Sl8sNWyNs7EIdFYN+dSDI4LcmXzG5SWmuqSogZRgV54fHxPXNEjFCZRwvzVh5Bzrsrh1+tMZlToXRPsSJIEk1mE0Sw2vjMRERERUQto0iItp06dwsaNG/H99983uF94eDgKCgrsthUUFCA8PLzB16nVaqjV6lrbJUnCvffei+Tk5EbXwrnQ1hNnMOOzvTCYRARqlXjxP/1sFdVaUpeEBJcVNZAJAp66thdyiqtwLL8MD397EB/dOQQ+asd+rOV6ExRyAWqF3KH9zaIEkyjCLEq2h6n6q5VaIYOPWgGFA9kxIiIiIiJ3adLd6PLlyxEaGorx48c3uF9SUhI2bdpkt23Dhg1ISkpqymkBnK+Y5qwVO9JhMIkYnhCM9XMva5UgBwAio6Lx6pvv2M0ZeuWNpU0eCqdRyvHKzf3QyUeF9LMVeHLtYbvAozElVUaYamRiJEmC0SxCZzSjXG9CSaURheV6nC7V4Wy5HsWVRpTpTKg0mKE3ibXOpTeJKKwwoFRnhOhEO4iIiIiIXMnpjI4oili+fDnuvPNOKBT2L58yZQqioqKwePFiAMCcOXMwcuRILFmyBOPHj8eqVauwd+9eLFu2rMkNlsvlSExMdPp1ldVryNw2LBahfpomn98VbpsyFaOuGIP0tDTEd+nS7AVHQ301ePXm/rj3833YmVqIpZtPYs6Yrg69VpKAc5VGKGQCTKIE0UVFCqoMZugMZmjVCnir5PUWrCAiIiIicgenMzobN25EZmYmpk2bVuu5zMxM5OXl2b4fPnw4vvzySyxbtgz9+/fHd999h7Vr1zZ5DR2ZTIYPPvjA6WFrAKAzWbIWjg7TcrfIqOYVNbhQzwg/PHVNLwDAl39l4ocDjhduECUJBrPosiDHSgJQoTfhbLmBi5USERERUYtq1jo6LcVaK/vff/9Fz549m3SMsW9sw7H8Mnw+/SJc2jWkyW2RJAmny/SN79hKPtqehg+3p0MuE/D2pIEYHBfY2k2ykcsE+KgV0CjbRrBJRERERO2P29fRaQ1RUVFNfq2+jWV03GX6JfG4slcYzKKEx77/B9nnKlu7STZmUUJJlRHnKgys0EZEREREbtWuAp3m0BstQ6c0Ss++ZEEQ8MT4nugV4YfSKhMe/uYgynVta80cg1lEUYUBJZVGpwonEBERERE5yrPv+mtoa3N03EmjlOPlm/ohxFeNjMJKPL72EExi28ug6ExmFJbrUcYKbURERETkYh0m0GkrGR1BABQyASq5DO6sQxbiq8arN/eDWiHDn2lFeGuT8yW5W4IES0W8sxV6VOhNaAdTxoiIiIioHeg4gU4LZHQEWCbcq+QyaJRyeKsV8NMoEaBVIthbhVBfNUJ9NQj2USPQW4UQXzX8vZRQK9wT9PQI98PC63oDAL7ek4Xv/852w1lcQ5IsC5ieLTdAZ2SFNiIiIiJqHqfX0WmPTGYRpuqhUWpF82M7tUIGmUyAXBAglwmQCQIUMgEymXPhiiAI0Cjl0CjlEEUJepNloU6DCyfqX94jFPeN7IL3t6bh1d9OIDZIiyGdg+z2yc3JRlpqKrokJLis3HVTiZKlYEGVwQxfjQIKeYeJxYmIiIjIhTrEXaQ1mwOg2aWNBUFAgFYFP40S3tWlklXVgU9zyGQCvFRyS6bHRw0/jRIqF93kTx3eGcm9LZXYFnx/CJlF5yuxffnZCgzp0x03XTsWQ/p0x5efrXDJOZvLWrCgnMPZiIiIiKgJOkSgU3MolCsyOu52YdDjq1FA2YygRxAE/N/VPdE70g+lOhMe+eYgynRG5OZk45E5syBWFyoQRRHz5s5Gbk7bGOJWc8FRvYnD2YiIiIjIcW3/rt8FrBkdlbz5mZeWJpMJ0KoUCPJWoZOPGj5qBRRNuAaNUo5XbuqHMD81ThVV4v/WHEbKyZO2IMfKbDYjPS3NVc13CVGSUFxpRHGlgeWoiYiIiMghHSrQaQ/ZnIbIZQK81QoE+6jRyUeNkjP52LF9q8MZmGAfNV69uT80Shn+Si/ClkJfyGT2fSKXyxHfpYtT7TKJItLPVuCf7GK3lrHWm0QUlrM6GxERERE1rkMUI7AOXVN70GKhK5Z/gnvuuQeiKEImk+HVN9/BbVOmNvq6bmG+eOa6Pnh09T9Yd6wYU55fgc+fuAtmsxlyuRyvvLG03oIEoiQhv0SHtDMVSD1TbvuaUVgBo9kSeHQO1uL+UQkY2S0EguD67JkES3W2KqOlWEFHWBeJiIiIiJwnSO3go/HS0lL4+/ujpKQEfn5+Tr/+QFYxJryzA1EBXtjx2OVuaGHLys7ORlxcnN2wM7lcjj2HjjlcNe3TnRl4d0sq5IKAJ8dEwU9XgPguXRAZFQ1JklBUYTgf0JytwMnT5Ug/W4FKQ91zZbyUcshkQIXe8nyfKD/MGpWIQXGBzb/gBmgUcvhqFO1uSCIRERERNY2jsQEzOu1QSkpKvXNrHA10piTFIf1sBX45nI8l2/Mx49Iu2HW4Aqm/70PamQoUVxnrfJ1CJqBzJ28khHijS4gPEkK8kRDig3B/DSr1Znzx5yl8tScTh3NKcf/KvzE8IRgzRyega6hvs6+7LjqTGfpyM3w0CmhVHeLtTNRiTGYRFXozJEjwUsmZQSUionalQ9wZWufoaDzkP+muXbtCJpPVyuj07N7V4WMIgoAFV/dA9rkqHMopwWsbTtg9LxOA6EAtulQHMtbAJibQq961bXw0Ctw3KgE3DYnGJ3+kY+2BXOxMLcSu1EIk9wnHvZd1QWSAV9MuugESgDKdqXrtHSVU7XwuFlFrE0UJ5QYTdAYzrCl/vUmEQmaCt1phWeTYDUNTiYiIXKljBDoeltGJjo7GsmXLcO+999rm1nzwwQfomRiP4kqD3bpBDVEr5HjpP33x3M9HIQG27EyXEG90DvZ2as2hCxcdfXRsD0y6KBbvb03FxqOn8evhfGz8twA3DorCtBHxCPRWNfHq62cSJZyrNECjlMNXzeFsRM6SJAkVBjMq9SbUNabZJFoW9JUJArQqefWQVf6eERFR29Qh5uj8eDAXD3y1Hxd3CcKqe5Lc0MLWkZ2djZMnTyIxMRHR0ZYha9b5NaYWLMP85WcrbOvx1FUY4WheKd7dkoq/0osAAFqVHJOHxWLSRbHwVrsn1hYEwE+jbPYCsUQdgSRJqDKaqxfodfx1AgCNSg6tUl5vppeIiMjVHI0NOkSg8+3eLMz77h+M6h6CFXdd5IYWti1m0RLsiC3wo83NycaQPt0dKozwV3oR3vn9JI7llwEAArVKTBsRjxsGRTVrQdSGaJRy+GkUHGZDVA9ddYDT3DWq1AoZtCoFh44SEZHbORobdIj/kXQeso6Oo+QyAQFaJVri1j4tNdXhRUcvig/C8ruG4vkJfRAd6IVzlUYs2XACt3zwJ9YfyXdLYKYzmlFYYYDR7L71fYjaI73JjMJyPUqqjC5ZiFdvEnGu0oDCcr2tAAwREVFr6hB3/tY5Oh1pGJNSLoO/Vun283RJSHBq0VGZIGBMrzB8fc/FmD+2O4K9VcgprsJTPxzBnZ/8hV2phS5fDNQsSjhXYUCF3uTS4xK1R0aziHMVBhRXGt0yxNU6j+dMmWVxX7EFh9ESERHV1DECnQ6W0bFSV68x406RUdF49c13IJdbgsjGFh21UshluHFQNFbfPxz3j0yAt1qOEwXlmPv1Acz+cj/SzpS7tJ3WhUaLKgwu+fSaqL0xixJKKo0oqjDA0AIZTlGSUK434Wy5HqU6I0zMqhIRUQvrUFXXOlJGx0qrUsAsSvUu9OkKt02ZilFXjEF6Wppt0VFHeankmDqiM24YGIUVOzPw7b4s7D11Drd/9BduGRqD6ZfGw8eFBQuMZhGFFXoWKqAOo65S0S1JAlBlMKPKYIYgAAqZDHJBgFwuQCETIJdZvnIeHRERuVrHCHQ6aEbHylejhFmUHC473RSRUdFOBTgX8tcqMWdMV9w8JBpvbEzB1hNn8OVfmVh/JB//vSIRY3uHu+xGSJKAkioj9EYRfl4sVECeKTMzC4ePHkdkbGdENON305UkyfJhgxEALhhJKhOqAx+7AEgGOctXExFRE3WIO3/rxNiOvKq3v5fSbZXNXCkywAsv39QPb9wyADFBXiisMGDhj//ivi/+RsrpMpeeS2cy42y5AQY3BoBELU0UJbzz/jLEx3fG+LFXYnCf7vjysxWt3axGiZIEg1lElcGMMp0JxZVGnC3X43SpzlI0odKIcr0JOqOZw0+JiMghbf/O1wWsmQyNhywY2hSCICDASwlZO8leJCUE48u7L8b9oxKgUcpwIKsYd368B0t+O44yndFl5xElyyKj5W4oVCBJEgwmERV6E+cnkNuZRQmlOiP+OZ6KB2bdb6uGKIoi5s2djdyc7FZuYdNIsBQ40JnMqNCbUFJlCYDOlOlRXGkpMmIwiS4vYkJERO1fh7jzZ0bHQiYTEKhVoi3HOrk52fhj21bk5mRDpZBh6vDO+PqeJIzuHgKzJOGbvdm4+f1d+PmfPJeWo66oLlTQnIBEFCXojGaU6SwTvs+U6W1BVBGrvpGbmMwiSqqMKCzXo8pgRqoTJd/bM1GyDMct15twrtLy+1ZYXfhAZzTzwwUiIuoYgQ4zOucp5DL4e7XMGjvO+vKzFRjSpztuunYshtQYbhPur8GL/+mHtyYNQFyQFucqjXjmp39x7+f7cKLAdcPZjGYRRRUGVDlYuMFcHdjYPmGuXpOk0mCG0SzaTfyuWfWNN2DkCkaziOJKAworDNAZzxcacLbku6ewZn6qDJbfycIKA06X6VBc/WGD3mRm1oeIqIPpEHf+54sRdOyMjpWl7LRr19gRgGYNi8vNycYjc2Y1ONxmWHwwVs4YhtmjE+GllOOf7BLc+clfeGX9cZRWuWY4mwSgVGdESaWx1vofRrOISoMJJZWWNULOVgc2zswZsAZTlQZmd6hpDCbLOjhFFYY6C4w0teT7hURJwsZ/C3D7R7tx67I/8fWerHaXlZQky7//FXrLnJ/TNX5vub4PEZHnczrQycnJwe23347g4GB4eXmhb9++2Lt3b737b9myBYIg1Hrk5+c3q+HOsA1dY0bHxkslh1bVtMBPEACVXAatSg5/LyWCvVUI9dMgxFcNfy9lk6okpTk43EYpl+GOpDh8fe/FGNMzFKIEfLfPMpztx4O5LhvOpjOZUVhhQJnOiHMVBpwu1aGowoAynQk6k7lZ55EAlOlMOMc1fTxOlcGMkkojynRGVBnMMJhEl91Q64xmFFUYcK6y8XVwbpsyFXsOHcPqn9Zjz6FjuG3KVKfOtTejCNNW7MHjaw8j5XQ50s9W4LUNJ3Dd0h14a1MK8kt0zbiS1mXNxBZV8vePiMjTOVVe+ty5cxgxYgRGjx6NX375BSEhIUhJSUFgYGCjrz1+/Dj8/Pxs34eGhjrf2iZiRqduvholRNFyU18fmSBAKRegkMugkAlQyhsu96pRyqFRylFlMKNcb3I4ILAOt6kZ7DQ03CbMT4Pnb+iLGwYW4ZX1x5FRWInnfz6KtftzMC+5O3pG+NX5OmeIknvXHzKYRRSW6+GjUUCr6hCV3m2k6vkVoiRBo5BD1s5LCOuMlve7/Y3z+feObf0Ymf3aMXIH1o/RGS2T8E1O3pQ3peT7iYIyvPt7KnalFQIAtCo5Jg+LRaBWhVV7spBZVImVuzOx6q8sXN4zFLddFItekc3/XWsNZlFCUYUBAdr2UZGSiIic59Td1UsvvYSYmBgsX77cti0+Pt6h14aGhiIgIMCpxrkKMzr18/NSwFwpwWgWIZcJUMpkUMgFKOSWvzf1BtRLJYdGKUOlwYwKgwmNxTvW4Tbz5s6G2Wx2eLjNkM5BWHn3MHy9NwsfbU/HkdxS3LV8DyYMjML9oxLg7+XaIXquZs3uWNb0aVo2rL0wmkUYTCL0JhGmGnOYymGCWmnJMLa3G05HgxDb+jF1xM32wY8lGJLLBEvFPsOFwZN75JVU4YOtafj1cD6k6jbdODAK0y6JR5C3CgBww6Ao7DxZiC//ysS+U+ew4d8CbPi3AP2j/XHbsFhc2jWk3b1/rVUXA7xUUHXQddaIiDyZIDkxO7NXr15ITk5GdnY2tm7diqioKMycORMzZsyo9zVbtmzB6NGjERcXB71ejz59+mDhwoUYMWKEw40sLS2Fv78/SkpK7LJCjhrz2lacPF2Or2ZcjKSEYKdf7+msbwF3LZwpSRIqDGZUOhDw5OZkIz0tDfFdujj9afSZMj3e3pyC9UcKAAAhPmo8dW0vXBQf1NSmtygBliybVxOHFLY1omhZF0VvsgQ4jmT3lNVDItUKWZteyFVvMqNCbyk64Q65OdlIS01Fl4SEZi3E25iSSiNW7MzAt/uyYDRbfj5jeobivpEJiAnS1vu6EwVl+OqvTPx2pMAW5EUFeOGWoTG4pl8EvNXtK0MpAPDzUkKj9IzfPSIiT+dobOBUoKPRaAAADz30EG6++Wbs2bMHc+bMwfvvv48777yzztccP34cW7ZswZAhQ6DX6/HRRx/h888/x+7duzFo0KA6X6PX66HX6+0uJiYmpsmBzqUvb0ZWURW+nzkcg2IbH2ZH7iGKEioMJlQZzHDnZ9R/nzqHxb8cQ2ZRJQDg1qExmDk6wemhiy11s3khlVzWbrM7BpMIQ3XmpjlBgCAAWpUCXkp5m+oH67pIjc2RaY4vP1thK8whk8nw6pvvOD3HpjE6oxlf78nCZ7tO2daQGhIXiNmXJ9Ya9tnQ78GZMj2+25eN7/dno7TKchwftQITBkZi4pAYhPlpXNpud/PtgMNIiYjaI7cEOiqVCkOGDMHOnTtt2x544AHs2bMHu3btcrhxI0eORGxsLD7//PM6n1+4cCEWLVpUa3tTA52hz2/EmTI9fn7gEvSO9Hf69eRaoiih3GCCzo0BT5XBjLc3p2D13zkAgC6dvLHo+t7oFubr0Otb4mazIYIA+KrbfnZHFCVbxkZvNjeasWsKjUIOL5W8VYcWmcxidYli95YGz83JxpA+3WvNVdtz6JhLgm2TKGLdP/lYtj0NZ8osHyZ1DfXB7MsTMSw+qFYWzdHfgyqDGesO5dnm8QCAXBDa5Tweb7UCPu0sI0VE1NG4JdCJi4vDlVdeiY8++si27b333sNzzz2HnJwchxs3b948/PHHH/UGR67O6PRbuB6lOhM2PTwSCSE+Tr+e3MMsSijXm2xzqNzhj5Nn8dxP/+JcpRFKuYD7RibgtmGxDZbCdvfNpjPUChl8Na2f3ZEkCSZRglmUIEqWrwaT6PQE+eZQyARoVQpolC03rM1kFlGhNzdYsMOV/ti2FTddO7bW9tU/rceISy9r8nElScL2lLN4d0sq0s9WAAAi/DW4d2QXJPcOr/P3oSm/B6Ik2c3jseof7Y/kRG/46wrQq3tii/8eOctLJYefi0vwExGR6zga6Dj1sdWIESNw/Phxu20nTpxAXFycU407cOAAIiIi6n1erVZDrVY7dcyG6GxV1zjZtC2RywT4eynhrZK77WbyksRO+GrGxXjhl6PYduIs3t58EjtOnsXC63rXO6ymoVLXLX2DpjeJMFTo4adx7/wBUZRgls4HMiZRsmyr3t4W1lk0iRJKdUaU6QEvpRxeSjkUbipe0BJBeF2crT7oiH+yi7F080kczC4BYClAMm1EPP4zKLrBLFlTfg9kgoBLunbCJV072c3jOZhdYju/eeNOxAV7Y2ivBMQFaxEbZHlEBGigkLWNf6OrDGaIomRZXLkNzxUjIqKGORXoPPjggxg+fDheeOEFTJw4EX/99ReWLVuGZcuW2fZZsGABcnJy8NlnnwEA3njjDcTHx6N3797Q6XT46KOPsHnzZvz222+uvZJ6SJLlk2cAnGjaRinkMvhrZdCa5ahw0fCgC+cVvPyffvjxYC5e35CCvzOLcduHu/Ho2O5I7h1e67XuuNlsDkkCSqqMqDSYIRMAAdU3XjXuv4QLNllvzmreoln3kSTALNUIZETJrXOmXE2SgEqDGZUGM9QKGbxUcihlsvN90MQb0+zsbBw/fgJRcfEICo1olT5pavXBuqSeLscH29Kw9cQZAJYPem69KAZTLu4MH03j//Q39/egW5gvnr62N/7Twxs3P7QY2j6XQ+HbCXKfYGTrgez99qMAFDIB0YFeiAv2tguA4oK1CNCqnLhy19CbRJyrNCLAS9nuy58TEXVUTgU6Q4cOxZo1a7BgwQI888wziI+PxxtvvIHJkyfb9snLy0NmZqbte4PBgIcffhg5OTnQarXo168fNm7ciNGjR7vuKhpQ86aZGZ22TSmXIUCrgt5kRpmu6WV165tXcP2AKAyKDcTTPx7BkdxSPPXDEfyRchaPju0O3xrDVFx5s+lK7qrw1RStVajhQvrqctUXEqr/ECDYBYGCIFR/rQ4YBcvfP1uxHHNm3d9qc7Jqum3KVIy6YkyTqw8ezinBip0Z2J5yFgAgE4Br+0fi7kvjEerreHEAV/0eFOedwrltn+Hcts8gqLRQBkdDGRSFSffPg+jdCaeKKpFVVAm9SURGYSUyCitrHcPPS4G4IG/EBmsR38kbvSP80DPCz+1z2Ixm0VJ+Wqtq9eGjRETkPKfm6LSW5pSXLqk0ov8zluxRyvPj2t06HR1ZpcGEcn3jJalrcmRegUkUsWJHBj75IwNmSUKYnxpPXdMLQzoH1TpWU282rSRJQkGpHoHeSo9ZsLa1CzW4Wluak9VUkiRhb8Y5rNiZgb3Vc2MEAJf3CMWMy7ogvpN3k4/d3N8DR/pXlCQUlOpwqrASmYWVOFVk/VqBglJ9nceVCwISQr3RO9IffaL80CfSH7HB2gbn3zWVTBAQqFW6bagkERE5xy1zdNojffW8D5lgGRpB7YdWpYBGIUeZE3MlHJlXoJDJcPelXXBxl2A8/eMRZJ+rwuwv9+O2YbG4b2SCbd5CU1aWN4sSTp4ux4GsYtujqMKAQK0S88f2wOgeoU4dr63Jzcm2BTkAIIoi5s2djVFXjGk3QcGF2tKcLGeJ1UUGPt2ZgSO5pQAsc9/G9QnHlKQ4xAU3PcCxasrvwYWvbywzJBMERPh7IcLfCxd3sV/rrMpgRta5Sny5Zh1W/bgeipDOUEd2B3w74URBOU4UlGNN9TA4H7UCvSP90CfK3/I10h/+2uYXFRAlCUVcWJSIqN3x+EBHZ7QWIpBzUmk7JKsuWKBVyVGmMzU6fMuZeQV9ovzx+fSL8ObGFKw9kIuVuzOxO70Iz1zXGwmhjlXnM5hEHM0rxf7qoOaf7GJU6GsHZecqjXjs+0O4slcY5l3V3SU3X62hPQcF9XHlnKyWGtJnEkVs/Pc0Pt2ZgbTqKmpqhQzXD4jE5GFxCPdvW+vXNGc4npdKDh9TCZYtmGb3M1L5h+L1Vb8ip0qBI7mlOJpXinK9CbvTi7A7vci2X0yQF/pE+tsCoK6hPk3KzEgSUFxpgL/Wc7KzRESezuMDHWtGR6Pkp3DtmVIuQ5C3ClUGM8r0xnqHszk7r0CrUmDB1T0xIrETXlh3FCdPl2Pq8j2YOToBtwyNqTUMpkJvwqGcEku2JrMY/+aV1pojolXJ0T86AANiAzAgJgCJoT74fNcpfL7rFDb8W4C9GUV4bFwPjOre/rI7ba1Qgyu4ai5KSwzpM5hE/PRPLr74MxM5xVUAAG+1HDcNjsatQ2MR5N3yk/Yd1ZzMUF0BtqHkNMLNZ/CfKyxlt01mESfPlONITikO55bgcE4pMosqkVVUhayiKvxyOB+AJSDsHu6LCQOiML5f/dU/6yLBMhzaz4vFbYiI2gOPn6NzKLsE1y79A+F+Gvz5f1e4qYXUkiTJUvq3qoEFR5syr6CwXI/n1x3FjpOFAICLOgdh7piuyDpXaRuGdiK/HOYLfmUCtUoMiAnAwNhAW2BT18Tlo3mlWPS/f23rmCT3DsPDV7ZsdscVGYcvP1tRKyhoz3N0rJozFyU3JxtDBw+Cz+DroOoUB2NxHszFeXj3tZcwqHscgrxVzcooVxpMWLs/Fyt3n8LZcgMAIMBLiUkXxeI/g6Psiml4oqbOoyqpMuLf3FIczinBkVxLAFSmM9mev3N4HO4fmdCkn42PWgFvLixKRNQq3LJgaGtpTqCzN6MIN72/C52Dtdgyr2UqvVHLMJlFlOlMMLiwGpkkSVizPwdvbEypt8x1ZIAGA2Is2ZqBMYGICfJy+EbJYBLx0R9p+HzXKYgSEOStwmPjemBktxCXXUN9XJlxaO4EdaNZxMGsYpwp12N4Qif4e7XfG/VKgwkvf7sdPx8vhUxT95BHH7XCUi65umxyXPXfYwK1DVYOK6ky4tu9Wfh6bxZKqyw36KG+atx+cRyuHxDp9qyCQibAR6OAAAGVBteUfm8qVwTYoiQhq6gS6w7lY8XODADANf0isODqHk1aw0ejkEMhFyCXCZAJAhQygaWoiYhaAAOdajtOnsXkj3aje5gv1j/Y9JXFqe3SGS3lqEUXvpUzCyux6KcjOJxTii6dvDEwNgD9q4Ob+hYarU9dWZQjuSV45n//2krpju0djoeu6lbvDX9zMzFtobJYfokOO1PPYldaIfZmnEOlwTKsVK2Q4areYbh5cAy6h/u2SFtcwWASsWZ/DpbvSMe5SqNl25kMVBzeBLlvCFTBUYjrezFOlxsbXJMn1Fdtt25MbLAWIb5q/Ho4H9//nWPrp5ggL0xJ6oxxfcLdXj1SJgjw1ShqBVJmUUKFwQRdA9lUd3JFJUSrHw/kYvEvRyFKloWFn7+hj0sCRwGWuYXWoEcuWAIhefXfGQgRETUfA51qm44WYPqne9Ev2h8/zr7ETS2k1iZJEioMZlTqTS69ATOaxWbdVDaURdGbzPhoezq++NOS3Qmuzu5cdkF2xxWZmD+2bcVN146ttX31T+sx4lL3fABgMIk4kFWMXWmF2JVaaBuyZxWoVSJAq7Lb3jfKHzcPicblPULbbCl4kyhi3aF8fLw9HfmlOgBAdKAXegvZ+Pj/psNsNtllHPQmM3LOVVlKJxedL52cWVSJkipjo+frGuqDqcM7Y3SPULev5SIIgLdKAa2q4eItkiTZFm115QcMLW3biTN4Yu1h6E0i+kb5Y8nN/d0+lLSuQEgpl7GaGxGRExjoVFt3KA8zV/6NizoH4Zv7ktzUQmorzKKEcp0JOpNj5ajdydEsyuGcEjz7U43sTp9wPHSlJbvjqkxMS2V0coursCu10Ja1qapRFlwmWCrdJXUJRlJCMHzNJUhPTYXBLxJbThmw6dhp2yKxgVolJgyMwg0DoxrNoLVUpTNRkvD7sdP4YGsaThVZflYhvmpMvyQe1/aLgEIuczrjUFJprA5+KmyBUGZhJXJLqtAtzBd3Du+MEQnBbq8YKcBS3cxbpXA646AzWgKetrSgrTMOZhXjkW8PolRnQudgLd6aNNDprK0ryAQBXio5NAoZ1+shImoE19GpZl1/Rc2qax2CXCbAX6uExiRDaZVrh7M5y9FSzH2i/PHZ9Ivw4bZ0rNx9Cr8ezrdVZhPyXFPO2VWVxS6kN5mxP7MYf1ZnbS5c1T7YW4WLE4KR1CUYF8UH2Ybm1ZWlmjt7EtYeyMWav3NwplyP5Tsy8NnOU7isWyfcPCQGg2IDat3wt0SlM0mSsCutEO9vScPxgjIAgL+XEncOj8N/BkXbDXdytrKYv1aJvlp/9I32d2mbnaFRyuGjVjQ5W6RRyqFRymEwiagymNvEhwzO6B8TgA/uGIw5qw4go7AS0z/di7duHYAuIY6VmHcVUZJQoTehQm+ZG2XtV3dn8YiIPJnHZ3S++isTC74/hDE9w/DRnUPc1EJqiyRJQqnO8cVGXa0pWZRDOSV4rkZ2Z1SCH76YMx7mqjKHj9FYm5o6x0FnNCOvRIfc4ipkFlViT0YR9p06Z1urCrCsVt8nyg/DEzohKSEYXcN8apXobqxfTGYRW0+cwXf7svF3ZrFtny6dvHHT4GiM7RMOb7WiRbJU+zPP4b0tqTiYXQLAUjp88rBY3HpRLHzaecUttUIGH7XC5dkDsyih0mBCldFcbxn4tqigVIcHvtqPjMJK+GoUWHJzf/SPCWjtZkEpl8FLKYdaIeP8HiKiaszoVNMzo9NhCYJlsVG1QoZSXf1r77hLU7IofaP88em0i/Dh9jR8uTsTW1JL0X3u5zj13QuoSPmr2ZmYhjIOBpOI/FJLIGMNaGp+Laow1Pm6Tj4qXNwlGMMTLFmbxkodN5bpUshluKJnGK7oGYbU0+X4bl82fjmcj7SzFXh5/XEs/f0krukXgc5intsWLz2WX4r3t6ZhV6ql1LhaIcNNg6MxJSkOAdq2u1aNI5RyS4DjrjkhcpkAX40SPmoFqoxmVOjbxzyeMD8Nlk0Zgoe/OYhDOSX471f78dyEPrXmzLU0o1mE0SxCAKBSyKCpDnq4ADYRUeM8PqPz/tZUvPjLMdw4KAqvTRzgngZSmyeKEkp1xlYpj9vULMqh7BI889O/yKyeDzI4VMCg+BAEBQZCVrOSU3Vp2/NVnWBf7rb6OevE50qDCbnFOuSWVCHP+rVEh7Nl+kYLOWhVckQGeCHCX4O+Uf6WrE2oj1M3XU3JxJTrTPj5UB6+25dt6w8A0GUcQOnfP6EqbR9gNjY7o5NxtgLLtqVh07HTlnbJBFzfPxJ3XdIZob4tP2/DleQyAT7q2pXUWoLOaEZV9Tyetv4fjs5oxuNrDuOPk2chE4DHxvXA9QOiWrtZdgQAaqUcGqUMagUXLiWijofFCKq9uTEFr288gduGxeKFG/q6qYXUXlQZzCjTNVzuty3RGc34YFsavtqd2SJt1ihliPD3QmSAps6vfhqFSz5JbuqaKKIkYU9GEb7dm40dJ89CvKBT5BCh1aigqq5iZfuqkEGtOP93ldxyg6iUC9XPyXGmTI/f/s2HKFluJJN7h+PuS+MRE6Rt9vW2JkGwrOOjVbWNBL5ZlGASRYiipYKd5XvJVoiiLTCJIl785Rj+dzAPAHDvZV1w14jObTKLIgiWeVJeSnmbrVRIRFQnawhiC0Ukh/9eWl4B/+BQDl2zToxVs3QnwVJZSikXUKoztYsqURqlHHOu6IrLu4dizYEc6AxmmEUJZslyYyhWf7U+RAm250Xx/D6mGvuqFfJ6A5lArbJFbuZumzIVo64Y43SmSyYIGBYfjGHxwcgtrsKa/TlY+3c2SvWW33MzZCjTmZrVtsu6dcK9lyUgMbRlJ6O7miAAWpUC3o2Uim5pliykNQtxPhshVb8/rUHP+a9iiw87VchkePzqngj2VmPFzgx8sC0NhRUGPHRltzZXHECSLB/gVBnMkMsEaFVyaBRyzuchoobZBRlS7a+OPldrP9Szra7XN4NB79BuHh/o6KsnSrfGcA1qmxRyGYK8VSjXm1y+7o679I1u3cpc7uBshbJarw/wwqzRibh/VAIqDWYYTKLtoTef/15vssxxsP79/D6i3WsA4PKeoegb1b77WSZYbnYbWwunrREEAQq5gLpGYoni+WDdaBZhMktuHwYnCALuH5WAYG8VXttwAt/ty0ZRhQELr+vVZoeLmUUJZToTymGCWiGHRtXyQ9skSYLeJEImCFwbiKipJAmQxPoDDYe+ovHnOwCPD3SY0aH6+KgVUCtkKKkytqlhM+QcmWCZewJ1a7ekdckEAd5qyxCm9hTgOEImE6CqzlDU/NDKGvQYqifsO/p77MzaSxOHxiDIW4WnfzyCzcdOo7jSgFdu6g8fTdv971OC5f8+nckMmWCCl8ryvnBXNspkFmEwi9AbRbsAVCYIUCtl0CjkrRb0GM2WjGBrB10m8/kPWmTVcysVNedSMgPXfknS+cAEFwYodQUron0gYntdxwtCWkLb/ZfaRZjRoYYo5TIEe6tQpjehyuDaMtQKmQBl9XwQvVFsd+uLUPvQmkUGWptSLoNSDnhVD4GTJGvQI9luvi+8Z2jK2ktjeoXB30uJR1f/g78zi3HfF/vw2MgwFOedcvtCtc11fn0eE1RyGbxUza/aZu1n6417fQGmKEm2YXUtFfSYRckus2v9+QsCoFZYrr2lqtaZzCJ0JhF6oxmmmn1Ux38FAiwBvTXoqRkEyWWCx3140WbUzJxI4vlHfQFLvQEMtVUeX4xg9pd/46d/8vD0tb1w14h4N7WQPIHeZG7WIqOWmy4BSrkluLnwEzqd0dwqZa7JMylkArw7aIDjDLNoGeZmMIs4lZmFAT27NnntpWP5pZi76gDOVRphKi5AwTdPQizJd8tCte4kCIBXdQEDR9dRsgYQepNlWGhz/hlzZdBjHSpnMDccdNVkLdVtDXxcmU2pN7hpppqVNeVyod2v4+VykgSI5joCFfGCQKauzAu1R6WVevhHxLMYgXUxw7Y6ppraDrVCjmBvy2T2xrIvAixzfaxVu1Tyxj8h1CjlUMkdOz5RfZRyGbzVcv6b5iBr8QONUo6zOaeatfZSj3A/PJ8cjRmf7IQyMALht7+CqpN/4Zmvt0PqPBS9Okciwt8LIb7qNle0oCZJAioNZlQazLYFSTXK2v+G1QxsXHnTXjPTY60ap1Y4Pp/IUCOwMTVhrpYEQF89Zw+w/E5ZMz1NWUDXGtzojGa3DYMWJQmiWYIRgGCC5wc6NTMsNQOYeh8MWKhuHv6bYvmUHrCUzSVqjEwmwF+rhNpov8ioNbBRKaqDGwcCG2eOT57BWsLaXKNymChKLpk0r5LL4O3GhT47gq5du0Imk9XK6MR36eLwMarOZCH/i3kIvXkh1OGJ8Ol3JQDgrT/ygT/yLceUCQjzUyPC37LmVIS/xrb+VFsLhKwLkpbpLGvzqOSyWsO+3Klm1ThBMNYZ9NQ3HM1VrH1Qrrf87Kznb+h3rSWCG48lmquDlwu/Spa/M3AhF/L8QIcZHWoCTfWaFDqj2TYkzZVjpJndOT8m3RNuEhoLQkxm0VYu2SxJMJvPVxFrjEYhh1bNNVJcITo6GsuWLcO9995rW8Pp9bffcWqOTZeEBEBXioKV8+HVdRgUAeFQBoRjePIEFOkk5JfoYBIly6K8xbo6j3FhINQ11AeX9wxt1UVpJViG1+qMrffv0YVBj1ouh1F0vMiEK5hFyZbtunBej1mUGNw0plYQI1b/XTy/jagFeX6gw4wONZG8eg6Eu9TM7pTpmj43qD2xflpac7ifwSSiymCG3mRuF6W+a3J0GJlCLquzbLIkXbBujNkSCJlEEWqFHN4qx+dQkGOmT5+O5ORknDx5EomJiYiOjkapzuhwMZLIqGi8+uY7mDd3NiqPbqux4G0SAMuN8tlyPfJKdMgrqUJese7830t09QZCb2xMwaC4QIztHY7RPULgq1G65frbC0lCq38IJEnngz8BaHf/PrmNaAZEE2A2Wr5aAxtmYqgN8vhiBGPf2IZj+WX4YvowXNK1k5taSNQ8oiihTG9q1U9T3cE66dc68beh4TqiKKHSaPk0t60HfR250pmnKq+uTOao3Jxspxe8BWoHQscyC/BnWiEySs/vo5QLGJHQCVf1DsOIxE58n5EdAUCoXwtk/+oKaEQjgxlqE1iMoJr1xlHNjA61YTKZAH8vJdSK9p/dkcuE6sDGsSINVrLq4MFHrYCuOuAxmMXGX9iCrGv2eKl44+lpfNQKyASgTOdYsNPUBW8tw9Y0CPPT4N/f1+KV6lLXyoBw3DTvVeSropB2tgJbTpzBlhNn4K2WY1T3UCT3DsOQuCC3zu1xZn0h8iCiuUYwU+PRjv8fIrLy+EDHWlVFwzk61A7Y5u60o+yOM1kbR2mUlipZJrOISqMZOkPrDmsTBMuNsCcuxknnaVUKyAQBpVVGt7/fcnOybev5AICxOB/fPHUn9hw6hgqlP347UoD1R/JRUKrHz//k4ed/8hDsrcKYXmEY2zscPSN8XfpebMr6Qu7WlgKvttQWlzFWAboSBjTk0Tw+zcGMDrU31uxOgFYJWRu8qRZgWcPFSyVHgFaJEF81ArQqaFUKl3/arJDL4KexnMNX4/rjN8Ya4IT4qKFVKRjkdAAapRwBWhXc/aNOS02tt9R111BfzBqdiLWzRuD92wfhxoFR8PNSoLDCgK/3ZOGuFXtw0/u7sGxbGjILK5vdlguDLlEUMW/ubOTmZDf72E315WcrMKRPd9x07VgM6dMdX362gm1xNc6poQ7A4+fo9HrqV1QazNg2bzRig7VuaiGRe7TW3B3ritzWFbrlNVbpduXiek2hN5mhM4hunagsAPBSyeGtUrT69VLrMJpFFFca3TaMNDcnG0P6dHd48VKjWcTutCKsP5KPbSlnbGvEAUCPcF+M7ROOq3qFIdhH7XRb/ti2FTddO7bW9tU/rceISy9z+njN5WzfdJS2WLlsjo6hAtCVNr4fURvEOTrVrEPXmNGh9sia3dEoL5i7c8G9l7O3YoIAKGQy2yrbckGATFa9rY3f2FvKvcrhIypQZTSj0mBy2YeSAgCNSg4fBjgdnlIuQ5C3CucqDW4pJVyzepu11PUrbyyt9+ZZKZfhkq6dcEnXTqg0mLDtxFmsP5KP3WlFOJZfhmP5ZXh700lc0rUTrusfiYsTgqCQOfb/XpeEhGavL+RKDWW7Wjq4cGVbPHL4G1Eb53Sgk5OTg/nz5+OXX35BZWUlEhMTsXz5cgwZMqTe12zZsgUPPfQQjhw5gpiYGDzxxBOYOnVqc9rtEKP5fP19NRfZo3ZMrZBD7eP4PLMLE7U1vxUEeMQQLGvlM2+VHKJ0/polVK87B8nuuq3bzv/d/jUC4Jbhd9R+yWUCgrQqFFcZYXRDYYzbpkzFqCvGOF29TatSYGyfcIztE45zFQZsOnYavxzOw+GcUmw9cQZbT5xBiI8a4/tF4Nr+EYgObHg0g7NBV32MZhEHs4qxO70I3moFLu8e2qSRFG0p8HJVW9riHCiijsCpoWvnzp3DwIEDMXr0aNx///0ICQlBSkoKEhISkJCQUOdr0tPT0adPH9x33324++67sWnTJsydOxc///wzkpOTHTpvU4euletN6PP0egDAsWfHskQnERE5TZIkFFca21wVwJpyc7Kx89BJHNf5YVt6OYqrjLbnBscF4voBkRjVPaTBNZ+aUjK7sFyPnamF2HHyLHanF6HygvWIEkN9cEWPUFzeIxSdO3k7fD1ffraiVuDlTGAgSRLSz1Zgd3oRqgxmJCUEo0d40wo4NLcth06k4z/3PwavbiOgCusCc3kRzCUFmDD2ciRGdUJUgBeiArwQGaBxaP0kDl0jcnzomlOBzmOPPYYdO3Zg+/btDjdk/vz5+Pnnn3H48GHbtltvvRXFxcX49ddfHTpGUwOds+V6DHluIwAg7YWrORSFiIiarKTK2CarIV6YLXjpjXcQedE4/HgwF7vTimxDW301CoztHY7rBkSiW5hvk84lShKO5pVix8lC7Ew9i6N5ZXbPB2qV6B/hhfyiMqQUSzDXuMNICPHG5T1CcUXPMMQ7EPQ4G3gVVRiwJ6MIu9OK8Fd6Ec6U6+2eD/FV47KunXBZtxAMjguE0onFeJ1ty+kyHTYfPY1Nx07jn+wSh8/jp1EgMsALkTWCH+vfw/01UMplDHSI4KZAp1evXkhOTkZ2dja2bt2KqKgozJw5EzNmzKj3NZdddhkGDRqEN954w7Zt+fLlmDt3LkpK6v7l1+v10OvP/wNVWlqKmJgYpwOdnOIqjHhxM1QKGU48N87h1xEREdWlTGeslbVoTY1Nls8v0eGnf3Lxv4N5yC/V2fbpEe6L6/pH4qreYY1mEcp0RvyVXmQLbs5VGu2e7xnhixEJnTAisRP2b/ge8+Zagi6F1g93PrUUFUFdsSe9CKYac526dPLGFT0tmZ4uIT5Nuna9yYyDWSXYnV6Iv9KLcKKg3O55tUKGXqEamPRVSCkBdKbz59eq5BieEIzLuoVgeEKwQ5mUxhSU6rD52GlsriO40eccRcWx7dBlHIRM6wdVYATumbcQpWYlckuqkHOuqla/XkgmAKG+GkQFaJAQ6oMnrukFv+a0m4EOtWNuKUaQlpaG9957Dw899BD+7//+D3v27MEDDzwAlUqFO++8s87X5OfnIywszG5bWFgYSktLUVVVBS8vr1qvWbx4MRYtWuRM0+pkKy3N+TlEROQCvhpL2fdyvWMLi7pbY5Plw/01uPvSLrhrRDz2ZBThfwdzseX4meoCBsfx5qYUXN4jFNcPiMSAmAAIgmAb9rUjtRA7T57FwewSu4IM3mo5hsUHY0RiMJK6BNsqveXmZNuCHAAwVZZixeN3Yc+hY/C5rje2p5zFpmMF2J1WhLSzFUjbno4Pt6ejc7AWV/QMwxU9QtElxLve4WWSJOHkmXL8lW7J2hzIKrYVHLLqFuaDYfHBGBYfhKPbfsRj1e2RKdWY9fz7EGL6Y3vKWRRVGLDx6GlsPHoacpmAQbEBuKxrCC7rFoJwf8ezJdbgZtPR0ziUYx/c9Iv2xxU9QjG6Ryg2rc3CvK9+tg1/WzxvJm67dqDd/pUGE3KLdcgtrkJOcRVyi6uQW6yz/V1vEpFfqkN+qQ77s4rx/A19HW4nUUflVEZHpVJhyJAh2Llzp23bAw88gD179mDXrl11vqZbt2646667sGDBAtu2devWYfz48aisrKwz0HFVRuff3FJc/dZ2hPiqsefxMQ6/joiIqCE6o7lFFhZtTFPKHxdXGvDL4Xz8eCAXaWcrbNtjgrzQPzoA+06dQ16Jzu41nYO1GJ7YCSMSgtE/JqDOYV+Olqku0xktQc/R09idXghjjfFtnYO1uLxHKC7vGYrEEB8UVRiwO90yFO2v9CIUVhjsjh3io8ZFXYIwLD4IQzsHIchb1Wi/hEdG4UhuKbadOINtJ84g44K1iLqF+diCnm5hPrUCr/qCGwHVwU3PMIzuEYJQX/uAqSlzoKwkSUJRhcEW9JhEYPol8U4doxZmdKgdc0tGJyIiAr169bLb1rNnT6xevbre14SHh6OgoMBuW0FBAfz8/OoMcgBArVZDrXZ+LYAL6U3M6BARketplHIoZAJMogSzKMEsSRBFCSbR8rWlAqCmVEwL0Kow6aJY3Do0BkdyS/HjwVxs+LcAWUVVyCqqAgCo5DIMjgvEiMRgDE/ohKjAuv+/rsnRCmW+GiWu7huBq/tGoFxnwvaTZ7Dp6Gn8mVaIjMJKfLIjA5/syECQtwpFFwQ2GqUMg2IDcVG8JbiJ71R3BqixTFffKH/0jfLHrNGJyCyqtAU9h3JKcKKgHCcKyvHRH+kI99Pg0q6dMDwxGBlnK7HpWAEO55wPDgQA/WMCbJmbEN/6710io6KbXFZaEAQE+6gR7KNGiFCOc/lZyM5WIjqaZaqJGuJUoDNixAgcP37cbtuJEycQFxdX72uSkpKwbt06u20bNmxAUlKSM6duEuuCagx0iIjI1RRyGeorYiZWBz/mFgiEmlqmWhAE9InyR58of8wd0xUbj57GqcIKDIgJwJC4IHipnKtU2pSgy0ejwLg+ERjXJwLlehP+qB7e9mdqEYoqDBAAdA/3xbAuQRgWH4y+Uf5QOfB/ujNloWODtLj94jjcfnEczlUY8MfJs9iWcga704qQX6rDt/uy8e2+7PP9BmBATACu6BmKUd0bDm5c7cLCE8uWLcP06dNb7PxE7Y1TQ9f27NmD4cOHY9GiRZg4cSL++usvzJgxA8uWLcPkyZMBAAsWLEBOTg4+++wzAOfLS8+aNQvTpk3D5s2b8cADD7RIeektx09j6vI96B3ph58fuNTh1xEREblTzUDIaBahM4rnFwRu55ozRMuqQm/CiYIyxHfyRoBW1aRjNLcstM5oxl/pRdiWcgZ7M84h3E+DK3paMjedfFouuLGqbzheRkZG0zI7HLpG7Zhbhq4NHToUa9aswYIFC/DMM88gPj4eb7zxhi3IAYC8vDxkZmbavo+Pj8fPP/+MBx98EG+++Saio6Px0UcfORzkNAczOkRE1BbJZAJkEKCUW4bB+Wosw611RhF6kxntOeZpzhAtK2+1AgNjA5t1jKZmuqw0Sjku62aZq9MW1Dcc7+TJk84HOmVlQEUxoHV63XiidsXpd/g111yDa665pt7nV6xYUWvbqFGjsH//fmdP1WzWOTpcKJSIiNo6tUIOtUIOSVJAbxKhtwY9rd2wdswVQVdbUd9wvMTEROcP9uWXwH33Af7+QHxc9aOz/dfICEDGD4qpffPoUN5adpIZHSIiai8EQYBGKYdGeT7oqTKYYTCLjb+YPFZdc6A++OCDpg1bsxaJKikBDvxjeVxIrQbiYoDOdQRCcTGAxgWLlhK5mWcHOrZ1dJjRISKi9qdm0COKEnTVw9uMDHo6pNumTMXoK8aguCAbiYmJTa+69tRTwAP3A8cOA+kZQPqp818zMoHMLECvB06ctDwuJAhAZLgl6OlcIwiyBkVBzRt2SOQqnh3oVGd0NEpmdIiIqH2TyQRoVQpoVYBZlKAzmlFlNNst5kmeLzIqGgN6NmG42oW0WqBXD8vjQmYzkJ1bOwhKzwAyTgFl5UBOnuXxRx3rKHJIHLURHh3o6JjRISIiDySXCfBWK+CtVlRXbTN7VOU2amVyuWV4WlwMMOqCqrWSBBQWnQ960k/ZB0P5BY4NiYvvDHSO5ZA4ciuPDnSY0SEiIk+nlMuglMvgq0F1wGO2/f9H5HKCAHQKtjyGDq79fGWlZfhbrWzQqaYPibN+DQxw55WRB+oQgY6aVdeIiKgDsM7naY2hbQIAlUIGtUIOUbKc38RhdXYEAe26dLhDGhoSZzJZhrtxSBy1EI8OdM4PXeMbn4iIOo6aQ9sM1VXb3FGqWi4TqoMbGVRyGQRBsD3nrVbYAq6WDnoEWDJdKoUMggDoqws4tEaMoZAJUCvl0ChkUMhlkKoXirUuGGv3kCTPDoQUCseGxF0YBKWfAgpOOz4krmZhBA6J69A8OtDRG61D15jRISKijkmlsNzwi6ICOpMZVYamBx3WAEKttAQ2CnnDHyTWDLisQY/e5J6qcdagSyW3BF41gy6tCpAk6fz6RGb3Lsp6YXBTkyAIUMiFem/AxAuDIEmC2Wz5KnpyhqyxIXEVldVzgjJqzwvKyuaQOKqTRwc6OhMzOkREREDNqm2WAgZV1ZmWxm74ZcL5rM2FAYQzzgc9cEnQIwiAWi63BXJyWcPtqlmqG1BCbzLbAh9XFHFoKLhxhkwmQAYB9X1G22Gr7Hlrgd49LY8LmUxAVk71ELga84Myqr+WVzQ8JC4goDoLFMshcR7GowMda0aHc3SIiIjOsxUwUNe9IKl12JdaYdnP1S4MevQOrA9UcziaygXtUivklqqsGsBoFquDHueyXa4KbpzRWEDXISkU5+fuXEiSgLOF9Q+JO30GKC4G9hcD+w/Wfv2FQ+KsQVDnOA6Jawc8OtBhRoeIiKh+NbMcZlGC0SxCJZdB1oI30/IL1geqGfQ0NAfIlayBn0/1EDu9yVzvvB6FzNJn6hYMbqgZBAEI6WR5XDSk9vMcEufRPDrQ4RwdIiIix8hlAuSy1v3/smbQI0mS2wIbR9sgihIMZsvwNoVcYHDjiRwdEmerDpfp/JC4uqrERYRzSFwL8OxAhxkdIiKidqk1gpwLyWQCNDI5PzDtqFpjSFx8HBAXa3mems2jAx2ddY4OAx0iIiIicpXGhsSVVwCnMps+JC4qov4hcQH+br00T+LRgY41o8NPYoiIiIhqEPghsFv5eDs3JC79FHDq1Pkhcdm5lsf2nbVfHxh4QWGEWA6Jq4dHBzrM6BARERHVQekFyBSAoRww6eHZK5W2Mc0dEnfunOXx94Har9doqofE1ZEJio3pcEPiPDrQ0ZtYjICIiIioTnIl4BUIiGbAWAkYKgHJ9Yu5khMcGRKXUZ35ybhgSFxmFqDTAcdTLI+6jh0daRkO10GGxHl4oMNiBEREREQNkskBtS+g8gGMVYChAhBNrd0qqouPN9Cnl+VxofqGxFnLZ1dUWp7PynFsSFzNoXHtdEicZwc6XDCUiIiIyDGCAKi0lodJbwl4TPrWbhU5ikPiavHYQMda+x4ANMzoEBERETlOobY8zCbAWGHJ9HAeT/vlzJA4u0yQM0PiOp8PtDq3jSFxHhvoWOfnAMzoEBERETWJXAHI/QGVr2Uej7HSMqeHPIvLhsTtqP36+obExccB4WFuHRLnwYHO+V9CZnSIiIiImkEmA9Q+loexylK4wGxo7VZRS2hsSNyZs5aAJyOzaUPiOsfal8h24ZA4Dw50LBkduUyAQs5Ah4iIiMgllF6Wh8lgGdbG8tQdlyAAoSGWx7ChtZ+va0icNSiyDok7dsLyqOvYFw6Js2WCIh1qnscGOjojK64RERERuY1CZXmIImCqsmR6zMbWbhW1Je4cEucAjw10uIYOERERUQuQyQCVt+VhNlbP5dFxTR5qmKND4uqqEnfmrGOncGV72xJmdIiIiIhamFxpKV6g9gNMOkuWhyWqyVmNDYk7XQh07dvoYTw20GFGh4iIiKiVCML5uTyiuTrLU8WKbeQaPj4O7eaxgQ4zOkRERERtgEwOqH0tD5PeEvSwgAG1AKeigIULF0IQBLtHjx496t1/xYoVtfbXaDTNbrQj9EZLRodr6BARERG1EQo14BUIeIcCGj/LUDciN3E6o9O7d29s3Ljx/AEUDR/Cz88Px48ft30vCIKzp2wS69A1ZnSIiIiI2hgWMKAW4HSgo1AoEB4e7vD+giA4tb+rcOgaERERUTtgLWCg8efQNnIpp6OAlJQUREZGokuXLpg8eTIyMzMb3L+8vBxxcXGIiYnB9ddfjyNHjjR6Dr1ej9LSUruHs1iMgIiIiKidsQ5t8wkDvAIs3xM1kVOBzrBhw7BixQr8+uuveO+995Ceno5LL70UZWVlde7fvXt3fPLJJ/jhhx/wxRdfQBRFDB8+HNnZ2Q2eZ/HixfD397c9YmJinGkmAGZ0iIiIiNota9U2bZAl6NH4AXJVa7eK2hlBkpqeFywuLkZcXBxee+01TJ8+vdH9jUYjevbsiUmTJuHZZ5+tdz+9Xg+9/nzN9dLSUsTExKCkpAR+fn4Ote29Lal46ddjuGlwNF69ub9DryEiIiKiNsxsAkxVlvk8oqm1W0OtpLRSD/+I+EZjg2aVlw4ICEC3bt1w8uRJh/ZXKpUYOHBgo/ur1Wqo1c1LVepNzOgQEREReRS5ApBbS1Ubzgc9LGJAdWhWFFBeXo7U1FREREQ4tL/ZbMahQ4cc3r85dNby0grO0SEiIiLyOAqVpYCBb5hlXo9SYxnyRlTNqUDnkUcewdatW5GRkYGdO3fihhtugFwux6RJkwAAU6ZMwYIFC2z7P/PMM/jtt9+QlpaGv//+G7fffjtOnTqFu+++27VXUQdrRkejZEaHiIiIyKMpNfZFDBj0EJwcupadnY1JkyahsLAQISEhuOSSS/Dnn38iJCQEAJCZmQmZ7Hxgce7cOcyYMQP5+fkIDAzE4MGDsXPnTvTq1cu1V1EHZnSIiIiIOhhrEQOll6U8tUkHGKsAs4HlqjsgpwKdVatWNfj8li1b7L5//fXX8frrrzvdKFdgRoeIiIioA6sZ9IiiJegx6Rj0dCDNKkbQllnX0WExAiIiIqIOTiYDVFrLo2bQY9I3/lpqtzw30LGuo8MFQ4mIiIjIqlbQU125zWxo7ZaRi3luoFOd0eHQNSIiIiKqk0wGqLwtD9FcPaeHQY+n8NhAR2fN6LAYARERERE1RiZn0ONhPDbQYUaHiIiIiJrELuipHt5m0nNOTzvjsYEOMzpERERE1Gx2w9tYva098dhAhxkdIiIiInKp+qq3Mehpkzw30OGCoURERETkLjWDHi5O2iZ5bKCjM1mHrjGjQ0RERERuVHNxUmvQY12nh0FPq/HYQMea0dFwHR0iIiIiaim1gh59jaBHbO3WdSgeGehIksSMDhERERG1LkEAlBrLA7APekRz67atA/DIQMdolmxZQjUzOkRERETUFijUlgcAmI2WOT0mPSCaWrddHsojAx296XyEzIwOEREREbU5cqXlAQBm0/m1eszG1m2XB/HIQEdnPD/+kYEOEREREbVpcgUg9wXUvpYhbSYdYKwuW01N5pGBjr7G/BxBEFq5NUREREREDpLJuUCpi3hkoKOzraHDbA4RERERtVO11urRVw9xM7CCmwM8MtCxZnRYWpqIiIiIPAIruDnNQwOd6oyOkhkdIiIiIvJAF1ZwswY9LGZg45GBjs5onaPDjA4REREReThrBTdrMQNr2eoOXszAIwMda0ZHw4wOEREREXUkMjmg9rE8OngxA88MdJjRISIiIqKOrs5iBtVD3DpAMQPPDHSY0SEiIiIiOq9WMQNDjWIGptZtm5t4ZKDDOTpERERERA1QqCwP4PwipSa95eEhPDLQYUaHiIiIiMhBFy5Saq45xK39zuvxzEDHtmAoMzpERERERA6TyQCZF6D0sgQ5ZkO7Xa/HIwOd80PXmNEhIiIiImoSQWjX6/V4ZKBzfugaMzpERERERC5x4Xo91ipubbR0tVMpj4ULF0IQBLtHjx49GnzNt99+ix49ekCj0aBv375Yt25dsxrsCGZ0iIiIiIjcSCa3lK3WBgE+YYBXoGW4m6ztJBqcjgR69+6NvLw82+OPP/6od9+dO3di0qRJmD59Ovbv348JEyZgwoQJOHz4cLMa3RhrRkfNjA4RERERkXtZS1d7BQA+oYB3J8uCpXJlqzbL6aFrCoUC4eHhDu375ptvYuzYsZg3bx4A4Nlnn8WGDRuwdOlSvP/++86e2mF6EzM6REREREStotYQN+u8npYd4uZ0JJCSkoLIyEh06dIFkydPRmZmZr377tq1C2PGjLHblpycjF27djV4Dr1ej9LSUruHM3S2qmsMdIiIiIiIWo21dPWFQ9wE99+nO3WGYcOGYcWKFfj111/x3nvvIT09HZdeeinKysrq3D8/Px9hYWF228LCwpCfn9/geRYvXgx/f3/bIyYmxplm2jI6LEZARERERNRG1Bzi5hvm9iFuTgU648aNw80334x+/fohOTkZ69atQ3FxMb755huXNmrBggUoKSmxPbKyspx6PTM6RERERERtnHV4m3cny9wejb+llLUguOTwzSovHRAQgG7duuHkyZN1Ph8eHo6CggK7bQUFBY3O8VGr1VCr1U61RRRFGAwGAIC33IwoXzl8FBJ0Op1TxyEi11OpVJDJ+MEDERER1cNaxU2ltczjMekBs75ZC5U2K9ApLy9Hamoq7rjjjjqfT0pKwqZNmzB37lzbtg0bNiApKak5p63FYDAgPT0domjJ5NzS0ws3dFUjWFaK9PQKl56LiJwnk8kQHx8PlUrV2k0hIiKits46xE2psXxvW6jUYClo4CCnAp1HHnkE1157LeLi4pCbm4unn34acrkckyZNAgBMmTIFUVFRWLx4MQBgzpw5GDlyJJYsWYLx48dj1apV2Lt3L5YtW+bMaRskSRLy8vIgl8sRExNj+dT4bAUMJjOiA73grW7dsnZEHZ0oisjNzUVeXh5iY2MhuCgdTURERB2ErYobAFEETEUOvcypQCc7OxuTJk1CYWEhQkJCcMkll+DPP/9ESEgIACAzM9NueMrw4cPx5Zdf4oknnsD//d//oWvXrli7di369OnjzGkbZDKZUFlZicjISGi1WgCATGGEADPUGi9o1M1KWhGRC4SEhCA3NxcmkwlKJT98ICIioiaSyQCFYyNEBElqwWLWTVRaWgp/f3+UlJTAz8/P7jmdTof09HR07twZXl5eAIBjeaUwmEUkhvpAq2KgQ9TaqqqqkJGRgfj4eGg0mtZuDhEREbVjDcUGNXnM7OCaw2FEqfY2Imo9/F0kIiKiluYxgU5N1iSVR14cERERERE1yiNjAbH6q6d8ity5c2e88cYbrd0Ml9myZQsEQUBxcXFrN4WIiIiIPJTHBTqSJJ3P6LSDOCcrKwvTpk1DZGQkVCoV4uLiMGfOHBQWFrZ201xi1KhRduXFAUuRiry8PPj7+7dOo4iIiIjI43lcoCPWKK3Q1jM6aWlpGDJkCFJSUvDVV1/h5MmTeP/997Fp0yYkJSWhqMix0nmuZjabbWsSuYNKpUJ4eHib//kQERERUfvlcYFOzSJyzmZ0srOz8fvvvyM7O9vFrarbrFmzoFKp8Ntvv2HkyJGIjY3FuHHjsHHjRuTk5ODxxx+37VtWVoZJkybB29sbUVFReOedd2zPSZKEhQsXIjY2Fmq1GpGRkXjggQdsz+v1ejzyyCOIioqCt7c3hg0bhi1bttieX7FiBQICAvDjjz+iV69eUKvV+Oijj6DRaGoNL5szZw4uv/xyAEBhYSEmTZqEqKgoaLVa9O3bF1999ZVt36lTp2Lr1q148803IQgCBEFARkZGnUPXVq9ejd69e0OtVqNz585YsmSJ3Xk7d+6MF154AdOmTYOvry9iY2Pt1mMyGAyYPXs2IiIioNFoEBcXZ1vPiYiIiIg6Ho8LdGwV1yA4lTH4+OOPERcXh8svvxxxcXH4+OOP3dRCi6KiIqxfvx4zZ860lcW2Cg8Px+TJk/H111/bArdXXnkF/fv3x/79+/HYY49hzpw52LBhAwBLkPD666/jgw8+QEpKCtauXYu+ffvajjd79mzs2rULq1atwj///IObb74ZY8eORUpKim2fyspKvPTSS/joo49w5MgRTJ48GQEBAVi9erVtH7PZjK+//hqTJ08GYCntPXjwYPz88884fPgw7rnnHtxxxx3466+/AABvvvkmkpKSMGPGDOTl5SEvLw8xMTG1+mLfvn2YOHEibr31Vhw6dAgLFy7Ek08+iRUrVtjtt2TJEgwZMgT79+/HzJkzcf/99+P48eMAgLfeegs//vgjvvnmGxw/fhwrV65E586dm/jTISIiIqJ2T2oHSkpKJABSSUlJreeqqqqkf//9V6qqqpIkSZJ0BpN0MOucdDi72OHjZ2VlSTKZTAJge8jlcikrK8tl13ChP//8UwIgrVmzps7nX3vtNQmAVFBQIMXFxUljx461e/6WW26Rxo0bJ0mSJC1ZskTq1q2bZDAYah3n1KlTklwul3Jycuy2X3HFFdKCBQskSZKk5cuXSwCkAwcO2O0zZ84c6fLLL7d9v379ekmtVkvnzp2r97rGjx8vPfzww7bvR44cKc2ZM8dun99//10CYDvObbfdJl155ZV2+8ybN0/q1auX7fu4uDjp9ttvt30viqIUGhoqvffee5IkSdJ///tf6fLLL5dEUay3bdR6Lvw9JSIiImqqhmKDmjwuo2MduOZMNiclJaXWnBSz2YyTJ0+6sGV1kxxcrzUpKanW90ePHgUA3HzzzaiqqkKXLl0wY8YMrFmzBiaTCQBw6NAhmM1mdOvWDT4+PrbH1q1bkZqaajueSqVCv3797M4xefJkbNmyBbm5uQCAlStXYvz48QgICABg6aNnn30Wffv2RVBQEHx8fLB+/XpkZmY61QdHjx7FiBEj7LaNGDECKSkpMJvNtm012ycIAsLDw3H69GkAlmFyBw4cQPfu3fHAAw/gt99+c6oNRERERORZPC7QEasDB2fmuXft2hUymX1XyOVyJCYmurJpdhITEyEIgi1YudDRo0cRGBiIkJCQRo8VExOD48eP491334WXlxdmzpyJyy67DEajEeXl5ZDL5di3bx8OHDhgexw9ehRvvvmm7RheXl61gsOhQ4ciISEBq1atQlVVFdasWWMbtgZYhtO9+eabmD9/Pn7//XccOHAAycnJMBgMTeyVhimVSrvvBUGwBaiDBg1Ceno6nn32WVRVVWHixIm46aab3NIOIiIiImr7PC7QsSZIZE5EOtHR0Vi2bBnkcjkAS5DzwQcfIDo62h1NBAAEBwfjyiuvxLvvvouqqiq75/Lz87Fy5UrccssttuDjzz//tNvnzz//RM+ePW3fe3l54dprr8Vbb72FLVu2YNeuXTh06BAGDhwIs9mM06dPIzEx0e4RHh7eaDsnT56MlStX4n//+x9kMhnGjx9ve27Hjh24/vrrcfvtt6N///7o0qULTpw4Yfd6lUpll5WpS8+ePbFjxw67bTt27EC3bt1sPxNH+Pn54ZZbbsGHH36Ir7/+GqtXr261ynVERERE1LoUrd0AV2tKRgcApk+fjuTkZJw8eRKJiYluDXKsli5diuHDhyM5ORnPPfcc4uPjceTIEcybNw9RUVF4/vnnbfvu2LEDL7/8MiZMmIANGzbg22+/xc8//wzAUjXNbDZj2LBh0Gq1+OKLL+Dl5YW4uDgEBwdj8uTJmDJlCpYsWYKBAwfizJkz2LRpE/r162cXuNRl8uTJWLhwIZ5//nncdNNNUKvVtue6du2K7777Djt37kRgYCBee+01FBQUoFevXrZ9OnfujN27dyMjIwM+Pj4ICgqqdY6HH34YQ4cOxbPPPotbbrkFu3btwtKlS/Huu+863JevvfYaIiIiMHDgQMhkMnz77bcIDw+3DbMjIiIioo6FGZ0aoqOjMWrUqBYJcgBLoLB371506dIFEydOREJCAu655x6MHj0au3btsgsKHn74YezduxcDBw7Ec889h9deew3JyckAgICAAHz44YcYMWIE+vXrh40bN+J///sfgoODAQDLly/HlClT8PDDD6N79+6YMGEC9uzZg9jY2EbbmJiYiIsuugj//POP3bA1AHjiiScwaNAgJCcnY9SoUQgPD8eECRPs9nnkkUcgl8vRq1cvhISE1Dl/Z9CgQfjmm2+watUq9OnTB0899RSeeeYZTJ061eG+9PX1xcsvv4whQ4Zg6NChyMjIwLp162oNSSQiIiKijkGQHJ0N34pKS0vh7++PkpIS+Pn52T2n0+mQnp6O+Ph4y7ovlQZkFlXCW61AQohPK7WYiGq68PeUiIiIqKkaig1q8riPu8VmZHSIiIiIiMgzeFygY01QMcwhIiIiIuq4PC7QYUaHiIiIiIg8LtCR0LSqa0RERERE5Dk8LtA5n9Fp3XYQEREREVHr8bhAxzZHhykdIiIiIqIOywMDHctXZnSIiIiIiDoujwt0RGZ0iIiIiIg6PI8LdJjRISIiIiIijwt0mNHxTFu2bIEgCCguLm7xc69YsQIBAQEtft76CIKAtWvXAgAyMjIgCAIOHDjQ5OO54hhEREREbY3HBTrtKaOTlZWFadOmITIyEiqVCnFxcZgzZw4KCwtbu2mYOnUqJkyY0NrNaDcEQbA9/P39MWLECGzevNnt542JiUFeXh769Onj0P51/VydPQYRERFRe+BxgU57yeikpaVhyJAhSElJwVdffYWTJ0/i/fffx6ZNm5CUlISioqLWbiI5afny5cjLy8OOHTvQqVMnXHPNNUhLS6tzX6PR6JJzyuVyhIeHQ6FQtOoxiIiIiNoajwt0bBmd1m1Go2bNmgWVSoXffvsNI0eORGxsLMaNG4eNGzciJycHjz/+uG3fzp0744UXXsC0adPg6+uL2NhYLFu2zO54WVlZmDhxIgICAhAUFITrr78eGRkZDbbhu+++Q9++feHl5YXg4GCMGTMGFRUVWLhwIT799FP88MMPtizFli1bHDqPNWOwaNEihISEwM/PD/fddx8MBkODbfn8888xZMgQ+Pr6Ijw8HLfddhtOnz5da799+/ZhyJAh0Gq1GD58OI4fP273/A8//IBBgwZBo9GgS5cuWLRoEUwmk+351157DX379oW3tzdiYmIwc+ZMlJeX2x1jxYoViI2NhVarxQ033OBwhi0gIADh4eHo06cP3nvvPVRVVWHDhg0ALIH3e++9h+uuuw7e3t54/vnnHWpvSkoKLrvsMmg0GvTq1ct2PKu6hp0dOXIE11xzDfz8/ODr64tLL70Uqamp9f5c6zrG1q1bcdFFF0GtViMiIgKPPfaYXbtGjRqFBx54AI8++iiCgoIQHh6OhQsXOtRPRERERC2hWfHAiy++CEEQMHfu3Hr3WbFihd2wHkEQoNFomnPaBplFETqjGVVGMyoNphZ9WNfwaUxRURHWr1+PmTNnwsvLy+658PBwTJ48GV9//bXd8ZYsWYIhQ4Zg//79mDlzJu6//37bTb7RaERycjJ8fX2xfft27NixAz4+Phg7dmy9AUZeXh4mTZqEadOm4ejRo9iyZQtuvPFGSJKERx55BBMnTsTYsWORl5eHvLw8DB8+3OHzbNq0yXbMr776Ct9//z0WLVrUYJ8YjUY8++yzOHjwINauXYuMjAxMnTq11n6PP/44lixZgr1790KhUGDatGm257Zv344pU6Zgzpw5+Pfff/HBBx9gxYoVtqACAGQyGd566y0cOXIEn376KTZv3oxHH33U9vzu3bsxffp0zJ49GwcOHMDo0aPx3HPPNdj2ulh/rjX7ZeHChbjhhhtw6NAhTJs2rdH2iqKIG2+8ESqVCrt378b777+P+fPnN3jenJwcXHbZZVCr1di8eTP27duHadOmwWQy1ftzresYV199NYYOHYqDBw/ivffew8cff1yrHz799FN4e3tj9+7dePnll/HMM8/UCsSIiIiIWkuTx6rs2bMHH3zwAfr169fovn5+fnafvLtzWFmVUcTED/502/Eb8u8zydCqGu/SlJQUSJKEnj171vl8z549ce7cOZw5cwahoaEAgKuvvhozZ84EAMyfPx+vv/46fv/9d3Tv3h1ff/01RFHERx99ZOvb5cuXIyAgAFu2bMFVV11V6xx5eXkwmUy48cYbERcXBwDo27ev7XkvLy/o9XqEh4fbtn3xxRcOnUelUuGTTz6BVqtF79698cwzz2DevHl49tlnIZPVHVvXDFi6dOmCt956C0OHDkV5eTl8fHxszz3//PMYOXIkAOCxxx7D+PHjodPpoNFosGjRIjz22GO48847bcd59tln8eijj+Lpp58GALugvHPnznjuuedw33334d133wUAvPnmmxg7dqwt+OnWrRt27tyJX3/9tc5216WyshJPPPEE5HK5ra0AcNttt+Guu+6yu+aG2rtx40YcO3YM69evR2RkJADghRdewLhx4+o99zvvvAN/f3+sWrUKSqXSdg1Wdf1cL/Tuu+8iJiYGS5cuhSAI6NGjB3JzczF//nw89dRTtp9hv379bP3atWtXLF26FJs2bcKVV17pcF8RERERuUuTMjrl5eWYPHkyPvzwQwQGBja6vyAICA8Ptz3CwsKaclqHSHAsq9IWOJoBAmAXUFr70zq06+DBgzh58iR8fX3h4+MDHx8fBAUFQafTITU1Fdu3b7dt9/HxwcqVK9G/f39cccUV6Nu3L26++WZ8+OGHOHfuXINtaOw8Vv3794dWq7V9n5SUhPLycmRlZWHlypV2bdm+fTsAy5C0a6+9FrGxsfD19bUFCJmZmfX2Q0REBADY9cMzzzxjd/wZM2YgLy8PlZWVAICNGzfiiiuuQFRUFHx9fXHHHXegsLDQ9vzRo0cxbNgwu3MmJSU19uMBAEyaNAk+Pj7w9fXF6tWr8fHHH9u1d8iQIbX6s6H2Hj16FDExMbYgx5G2HDhwAJdeeqktyGmKo0ePIikpye4DiREjRqC8vBzZ2dm2bRd+yBEREVHncEMiIiKi1tCkjM6sWbMwfvx4jBkzxqFhPeXl5YiLi4Moihg0aBBeeOEF9O7duymnbpRKLsM3916MhBBveDmQXXElL6Xcof0SExMhCAKOHj2KG264odbzR48eRWBgIEJCQmzbLrxxFQQBoigCsPTv4MGDsXLlylrHCgkJgUqlspt/ERYWBrlcjg0bNmDnzp347bff8Pbbb+Pxxx/H7t27ER8fX2e7GzuPI6677jq7QCIqKgoVFRVITk5GcnIyVq5ciZCQEGRmZiI5ObnW0Lua/WC9Ea/ZD4sWLcKNN95Y67wajQYZGRm45pprcP/99+P5559HUFAQ/vjjD0yfPh0Gg8EuOGuK119/HWPGjIG/v3+d/eHt7W33fWPtbYoLh0K6U0PvSSIiIqLW5nQksGrVKvz999/Ys2ePQ/t3794dn3zyCfr164eSkhK8+uqrGD58OI4cOYLo6Og6X6PX66HX623fl5aWOtw+CYBGKYe3WgmNg4FHSwsODsaVV16Jd999Fw8++KDdzWl+fj5WrlyJKVOmODzEb9CgQfj6668RGhoKPz+/OvdJTEystU0QBIwYMQIjRozAU089hbi4OKxZswYPPfQQVCoVzGaz0+cBLJmKqqoq23X9+eef8PHxQUxMDGQyGXx9fe3237dvHwoLC/Hiiy8iJiYGALB3716Hrv3C9h0/frzOa7WeRxRFLFmyxDb86ptvvrHbp2fPnti9e7fdtj//dGwoZHh4eL3nbkp7e/bsiaysLOTl5dmyV421pV+/fvj0009hNBrrzOrU9XOt67yrV6+GJEm29+COHTvg6+tb7+8sERERUVvj1NC1rKwszJkzBytXrnT4E+ekpCRMmTIFAwYMwMiRI/H9998jJCQEH3zwQb2vWbx4Mfz9/W0P682vI9rLOjpLly6FXq9HcnIytm3bhqysLPz666+48sorERUVZTeBvjGTJ09Gp06dcP3112P79u1IT0/Hli1b8MADD9gNNapp9+7deOGFF7B3715kZmbi+++/x5kzZ2zzhjp37ox//vkHx48fx9mzZ2E0Gh0+j8FgwPTp0/Hvv/9i3bp1ePrppzF79ux65+fExsZCpVLh7bffRlpaGn788Uc8++yzTvSmxVNPPYXPPvsMixYtwpEjR3D06FGsWrUKTzzxBABLsGc0Gm3n+fzzz/H+++/bHeOBBx7Ar7/+ildffRUpKSlYunSpU/NzXNneMWPGoFu3brjzzjtx8OBBbN++3a4aX11mz56N0tJS3Hrrrdi7dy9SUlLw+eef2+bI1fVzvdDMmTORlZWF//73vzh27Bh++OEHPP3003jooYfq/RkSERERtTVO3bXs27cPp0+fxqBBg6BQKKBQKLB161a89dZbUCgUjX5SDFiGuwwcOBAnT56sd58FCxagpKTE9sjKynKofZIktZt1dLp27Yq9e/eiS5cumDhxIhISEnDPPfdg9OjR2LVrF4KCghw+llarxbZt2xAbG4sbb7wRPXv2xPTp06HT6erNvPj5+WHbtm24+uqr0a1bNzzxxBNYsmSJbaL7jBkz0L17dwwZMgQhISHYsWOHw+e54oor0LVrV1x22WW45ZZbcN111zVYejgkJAQrVqzAt99+i169euHFF1/Eq6++6vD1WyUnJ+Onn37Cb7/9hqFDh+Liiy/G66+/biu20L9/f7z22mt46aWX0KdPH6xcuRKLFy+2O8bFF1+MDz/8EG+++Sb69++P3377zRZ4uFpj7ZXJZFizZg2qqqpw0UUX4e677240AA4ODsbmzZtRXl6OkSNHYvDgwfjwww9t2Z26fq4XioqKwrp16/DXX3+hf//+uO+++zB9+nS39QMRERGROwiSEzPiy8rKcOrUKbttd911F3r06IH58+c7tLK62WxG7969cfXVV+O1115z6LylpaXw9/dHSUlJrRt3nU6H9PR0xMfHQ6VW43BOCQCgd6Qf5Pz0ucVNnToVxcXFWLt2bWs3hdqQmr+n7iwvT0RERJ6vodigJqfm6Pj6+tYKZry9vREcHGzbPmXKFERFRdk+KX/mmWdw8cUXIzExEcXFxXjllVdw6tQp3H333c5eU6PEGjFbW8/oEBERERGR+7i8LFlmZqbdOP5z585hxowZyM/PR2BgIAYPHoydO3eiV69erj41auamGOYQEREREXVcTg1day2ODl2TKZQ4ll8GmSCgT5R/K7WWiC7EoWtERETkKo4OXfOoSSxidcjGUWtERERERB2bRwU61uSUjJEOEREREVGH5lGBDjM6REREREQEeFigw4wOEREREREBHhbo2DI6rdsMIiIiIiJqZR4V6DCjQ0REREREgIcFOmL1V8Y5REREREQdm0cFOu0pozN16lQIglDrcfLkydZuWpOsWLECAQEBrd0MIiIiIiIAgKK1G+BK7a3q2tixY7F8+XK7bSEhIU4fx2AwQKVSuapZRERERETtHjM6rUitViM8PNzuIZfLsXXrVlx00UVQq9WIiIjAY489BpPJZHvdqFGjMHv2bMydOxedOnVCcnIyAODw4cMYN24cfHx8EBYWhjvuuANnz561vU4URbz88stITEyEWq1GbGwsnn/+edvz8+fPR7du3aDVatGlSxc8+eSTMBqNtucPHjyI0aNHw9fXF35+fhg8eDD27t2LLVu24K677kJJSYktM7Vw4UL3dyARERERUT08L6MjSZBVVgBqqeUboNU2O52Uk5ODq6++GlOnTsVnn32GY8eOYcaMGdBoNHbBw6effor7778fO3bsAAAUFxfj8ssvx913343XX38dVVVVmD9/PiZOnIjNmzcDABYsWIAPP/wQr7/+Oi655BLk5eXh2LFjtmP6+vpixYoViIyMxKFDhzBjxgz4+vri0UcfBQBMnjwZAwcOxHvvvQe5XI4DBw5AqVRi+PDheOONN/DUU0/h+PHjAAAfH59m9QMRERERUXMIkjUN0oaVlpbC398fJSUl8PPzs3tOp9MhPT0d8fHxKNZLOFNQhL7do1unoeXlgLe3Q7tOnToVX3zxBTQajW3buHHj0K1bN6xevRpHjx6FUB00vfvuu5g/fz5KSkogk8kwatQolJaW4u+//7a99rnnnsP27duxfv1627bs7GzExMTg+PHjiIiIQEhICJYuXYq7777boTa++uqrWLVqFfbu3QsA8PPzw9tvv40777yz1r4rVqzA3LlzUVxc7NCxqWOp+Xta8z1PRERE5KyGYoOaPC+j046MHj0a7733nu17b29vzJo1C0lJSbYgBwBGjBiB8vJyZGdnIzY2FgAwePBgu2MdPHgQv//+e52ZlNTUVBQXF0Ov1+OKK66otz1ff/013nrrLaSmpqK8vBwmk8nuzfPQQw/h7rvvxueff44xY8bg5ptvRkJCQpOvn4iIiIjIXTwq0JEkCZKXFgW5ZxHm1wqfGmu1Tu3u7e2NxMTEJp3K+4LMUXl5Oa699lq89NJLtfaNiIhAWlpag8fbtWsXJk+ejEWLFiE5ORn+/v5YtWoVlixZYttn4cKFuO222/Dzzz/jl19+wdNPP41Vq1bhhhtuaNI1EBERERG5i0cFOqIEQBAg+HgD3u1zeEzPnj2xevVqSJJky+rs2LEDvr6+iI6uf0jeoEGDsHr1anTu3BkKRe0fa9euXeHl5YVNmzbVOXRt586diIuLw+OPP27bdurUqVr7devWDd26dcODDz6ISZMmYfny5bjhhhugUqlgNpubcslERERERC7nYVXXLF9laB9V1+oyc+ZMZGVl4b///S+OHTuGH374AU8//TQeeughyGT1/7hmzZqFoqIiTJo0CXv27EFqairWr1+Pu+66C2azGRqNBvPnz8ejjz6Kzz77DKmpqfjzzz/x8ccfA7AEQpmZmVi1ahVSU1Px1ltvYc2aNbbjV1VVYfbs2diyZQtOnTqFHTt2YM+ePejZsycAoHPnzigvL8emTZtw9uxZVFZWurejiIiIiIga4FGBjlgd6bST6tJ1ioqKwrp16/DXX3+hf//+uO+++zB9+nQ88cQTDb4uMjISO3bsgNlsxlVXXYW+ffti7ty5CAgIsAVITz75JB5++GE89dRT6NmzJ2655RacPn0aAHDdddfhwQcfxOzZszFgwADs3LkTTz75pO34crkchYWFmDJlCrp164aJEydi3LhxWLRoEQBg+PDhuO+++3DLLbcgJCQEL7/8spt6iIiIiIiocR5VdS2v3IwynRHRgVoEeXMBTaK2glXXiIiIyFUcrbrmkRkdWTvO6BARERERUfN5VKBjm6PTnseuERERERFRs3lUoOMJc3SIiIiIiKj5PCrQYUaHiIiIiIgADwt0mNEhIiIiIiLAgwIdSZKY0SFqo9pBcUciIiLyMIrWbkBzKZVKCIKAM2fOwCRqIAEw6HSAWd7aTSMiWIKcM2fOQBAEKJXK1m4OERERdRDtPtCRy+WIjo5GdnY2CorOWrZVaCBnjWmiNkMQBERHR0Mu5wcQRERE1DLafaADAD4+PojvkoC7vtsEmQB8f/9w+Gu5YChRW6FUKhnkEBERUYvyiEAHAEySgLxyMwDA10cLjcpjLo2IiIiIiJzUrGIEL774IgRBwNy5cxvc79tvv0WPHj2g0WjQt29frFu3rjmnrZPeJNr+rlbwk2MiIiIioo6syYHOnj178MEHH6Bfv34N7rdz505MmjQJ06dPx/79+zFhwgRMmDABhw8fbuqp66QzWrI5SrnA+TlERERERB1ckwKd8vJyTJ48GR9++CECAwMb3PfNN9/E2LFjMW/ePPTs2RPPPvssBg0ahKVLlzapwfWxZnQ0zOYQEREREXV4TZrIMmvWLIwfPx5jxozBc8891+C+u3btwkMPPWS3LTk5GWvXrq33NXq9Hnq93vZ9SUkJAKC0tLTe1xSeK4Wor4RCoWxwPyIiIiIiar+s9/qNrdPndKCzatUq/P3339izZ49D++fn5yMsLMxuW1hYGPLz8+t9zeLFi7Fo0aJa22NiYho9XxYA/+cdahoREREREbVTZWVl8Pf3r/d5pwKdrKwszJkzBxs2bIBGo2l24+qzYMECuyxQcXEx4uLikJmZ2eDFNGbo0KEOB2jt5TiuOEZpaSliYmKQlZUFPz+/Vm1LWzsO+9e9x2lL/euq9rB/3dse9q9729OW+tcVx2H/uvc47F/3Hof9Wz9JkjB48GBERkY2uJ9Tgc6+fftw+vRpDBo0yLbNbDZj27ZtWLp0KfR6fa21MsLDw1FQUGC3raCgAOHh4fWeR61WQ61W19ru7+/frB+0XC5v9hulrR3HVW0BAD8/P/avm9oCsH/d2Rag+f3rqvawf93bHvave9vTlvrXlcdh/7r3OOxf9x6H/Vs3lUoFmazhcgNOFSO44oorcOjQIRw4cMD2GDJkCCZPnowDBw7UuSBgUlISNm3aZLdtw4YNSEpKcubULjFr1iyPO46r2uIKbalfXHUc9q97j9OW+hdoW9fUltriKm3pmtpSW1ylLV1TWzuOK7B/3Yv9614dtX8FqbFZPI0YNWoUBgwYgDfeeAMAMGXKFERFRWHx4sUALOWlR44ciRdffBHjx4/HqlWr8MILL+Dvv/9Gnz59HDpHaWkp/P39UVJS4rJPz+g89q97sX/di/3rXuxf92L/uhf7173Yv+7F/m2+Zi0YWpfMzEzk5eXZvh8+fDi+/PJLLFu2DP3798d3332HtWvXOhzkAJahbE8//XSdw9mo+di/7sX+dS/2r3uxf92L/ete7F/3Yv+6F/u3+Zqd0SEiIiIiImprXJ7RISIiIiIiam0MdIiIiIiIyOMw0CEiIiIiIo/DQIeIiIiIiDxOiwU627Ztw7XXXovIyEgIgoC1a9faPV9QUICpU6ciMjISWq0WY8eORUpKiu35jIwMCIJQ5+Pbb7+17ZeZmYnx48dDq9UiNDQU8+bNg8lkaqnLbDUt1b91Pb9q1aqWusxW09z+BYD8/HzccccdCA8Ph7e3NwYNGoTVq1fb7VNUVITJkyfDz88PAQEBmD59OsrLy919ea2upfq3c+fOtd6/L774orsvr9W5on9TU1Nxww03ICQkBH5+fpg4cWKtxaD5/nVv/3bU9+/ixYsxdOhQ+Pr6IjQ0FBMmTMDx48ft9tHpdJg1axaCg4Ph4+OD//znP7X6z5H7gy1btmDQoEFQq9VITEzEihUr3H15ra6l+nfLli113kPk5+e3yHW2Flf17wMPPIDBgwdDrVZjwIABdZ7rn3/+waWXXgqNRoOYmBi8/PLL7rqsdqPFAp2Kigr0798f77zzTq3nJEnChAkTkJaWhh9++AH79+9HXFwcxowZg4qKCgBATEwM8vLy7B6LFi2Cj48Pxo0bBwAwm80YP348DAYDdu7ciU8//RQrVqzAU0891VKX2Wpaon+tli9fbrffhAkTWuISW1Vz+xewrDF1/Phx/Pjjjzh06BBuvPFGTJw4Efv377ftM3nyZBw5cgQbNmzATz/9hG3btuGee+5pkWtsTS3VvwDwzDPP2L1///vf/7r9+lpbc/u3oqICV111FQRBwObNm7Fjxw4YDAZce+21EEXRdiy+f93bv0DHfP9u3boVs2bNwp9//okNGzbAaDTiqquusvv9f/DBB/G///0P3377LbZu3Yrc3FzceOONtucduT9IT0/H+PHjMXr0aBw4cABz587F3XffjfXr17fo9ba0lupfq+PHj9u9h0NDQ1vkOluLK/rXatq0abjlllvqPE9paSmuuuoqxMXFYd++fXjllVewcOFCLFu2zG3X1i5IrQCAtGbNGtv3x48flwBIhw8ftm0zm81SSEiI9OGHH9Z7nAEDBkjTpk2zfb9u3TpJJpNJ+fn5tm3vvfee5OfnJ+n1etdeRBvmrv6t69gdUVP719vbW/rss8/sjhUUFGTb599//5UASHv27LE9/8svv0iCIEg5OTluupq2x139K0mSFBcXJ73++utua3t70JT+Xb9+vSSTyaSSkhLbPsXFxZIgCNKGDRskSeL718pd/StJfP9anT59WgIgbd26VZIkS18plUrp22+/te1z9OhRCYC0a9cuSZIcuz949NFHpd69e9ud65ZbbpGSk5PdfUltirv69/fff5cASOfOnWu5i2mDmtK/NT399NNS//79a21/9913pcDAQLv73fnz50vdu3d3/UW0I21ijo5erwcAaDQa2zaZTAa1Wo0//vijztfs27cPBw4cwPTp023bdu3ahb59+yIsLMy2LTk5GaWlpThy5IibWt/2uap/rWbNmoVOnTrhoosuwieffAKpgy/F5Gj/Dh8+HF9//TWKioogiiJWrVoFnU6HUaNGAbC8fwMCAjBkyBDba8aMGQOZTIbdu3e3zMW0Qa7qX6sXX3wRwcHBGDhwIF555ZUOMbS1IY70r16vhyAIdovWaTQayGQy2z58/9bNVf1rxfcvUFJSAgAICgoCYPn/ymg0YsyYMbZ9evTogdjYWOzatQuAY/cHu3btsjuGdR/rMToKd/Wv1YABAxAREYErr7wSO3bscPfltDlN6V9H7Nq1C5dddhlUKpVtW3JyMo4fP45z5865qPXtT5sIdKw/0AULFuDcuXMwGAx46aWXkJ2djby8vDpf8/HHH6Nnz54YPny4bVt+fr7dLxkA2/eePga0Ia7qX8AybOKbb77Bhg0b8J///AczZ87E22+/3RKX0WY52r/ffPMNjEYjgoODoVarce+992LNmjVITEwEYHmPXpjCVygUCAoK4vvXBf0LWMY4r1q1Cr///jvuvfdevPDCC3j00Udb47LaDEf69+KLL4a3tzfmz5+PyspKVFRU4JFHHoHZbLbtw/dv3VzVvwDfvwAgiiLmzp2LESNGoE+fPgAs7z2VSoWAgAC7fcPCwmzvPUfuD+rbp7S0FFVVVe64nDbHnf0bERGB999/H6tXr8bq1asRExODUaNG4e+//3bzVbUdTe1fR/AeuG5tItBRKpX4/vvvceLECQQFBUGr1eL333/HuHHjIJPVbmJVVRW+/PLLOrMNVJsr+/fJJ5/EiBEjMHDgQMyfPx+PPvooXnnllZa4jDbL0f598sknUVxcjI0bN2Lv3r146KGHMHHiRBw6dKgVW9/2ubJ/H3roIYwaNQr9+vXDfffdhyVLluDtt9+2fereETnSvyEhIfj222/xv//9Dz4+PvD390dxcTEGDRpU578hdJ4r+5fvX8uIgsOHD3eIIjitwZ392717d9x7770YPHgwhg8fjk8++QTDhw/H66+/7vJztVV8/7Y8RWs3wGrw4ME4cOAASkpKYDAYEBISgmHDhtkNg7D67rvvUFlZiSlTpthtDw8Px19//WW3zVq1Ijw83H2Nbwdc0b91GTZsGJ599lno9Xq7YRcdTWP9m5qaiqVLl+Lw4cPo3bs3AKB///7Yvn073nnnHbz//vsIDw/H6dOn7Y5rMplQVFTE968L+rcuw4YNg8lkQkZGBrp3795i19PWOPLvw1VXXYXU1FScPXsWCoUCAQEBCA8PR5cuXQCA798GuKJ/69LR3r+zZ8+2FbmIjo62bQ8PD4fBYEBxcbHdp+IFBQW2954j9wfh4eG1Kl0VFBTAz88PXl5e7rikNsXd/VuXiy66qN4h9J6mOf3riPrev9bnOqo291Gcv78/QkJCkJKSgr179+L666+vtc/HH3+M6667DiEhIXbbk5KScOjQIbv/bDds2AA/Pz/06tXL7W1vD5rTv3U5cOAAAgMDO3SQU1N9/VtZWQkAtT79lsvltqpKSUlJKC4uxr59+2zPb968GaIoYtiwYS10BW1bc/q3LgcOHIBMJvP4qj+OcuTfh06dOiEgIACbN2/G6dOncd111wHg+9cRzenfunSU968kSZg9ezbWrFmDzZs3Iz4+3u75wYMHQ6lUYtOmTbZtx48fR2ZmJpKSkgA4dn+QlJRkdwzrPtZjeKqW6t+6HDhwABERES6+orbFFf3riKSkJGzbtg1Go9G2bcOGDejevTsCAwObfyHtVUtVPSgrK5P2798v7d+/XwIgvfbaa9L+/fulU6dOSZIkSd988430+++/S6mpqdLatWuluLg46cYbb6x1nJSUFEkQBOmXX36p9ZzJZJL69OkjXXXVVdKBAwekX3/9VQoJCZEWLFjg9utrbS3Rvz/++KP04YcfSocOHZJSUlKkd999V9JqtdJTTz3l9utrbc3tX4PBICUmJkqXXnqptHv3bunkyZPSq6++KgmCIP3888+2/caOHSsNHDhQ2r17t/THH39IXbt2lSZNmtTi19vSWqJ/d+7cKb3++uvSgQMHpNTUVOmLL76QQkJCpClTprTKNbckV/z78Mknn0i7du2STp48KX3++edSUFCQ9NBDD9ntw/ev+/q3I79/77//fsnf31/asmWLlJeXZ3tUVlba9rnvvvuk2NhYafPmzdLevXulpKQkKSkpyfa8I/cHaWlpklarlebNmycdPXpUeueddyS5XC79+uuvLXq9La2l+vf111+X1q5dK6WkpEiHDh2S5syZI8lkMmnjxo0ter0tzRX9K0mW+7P9+/dL9957r9StWzfbvznWKmvFxcVSWFiYdMcdd0iHDx+WVq1aJWm1WumDDz5o0etta1os0LGWFbzwceedd0qSJElvvvmmFB0dLSmVSik2NlZ64okn6iwJvWDBAikmJkYym811nicjI0MaN26c5OXlJXXq1El6+OGHJaPR6M5LaxNaon9/+eUXacCAAZKPj4/k7e0t9e/fX3r//ffr/Vl4Elf074kTJ6Qbb7xRCg0NlbRardSvX79a5ZALCwulSZMmST4+PpKfn5901113SWVlZS11ma2mJfp337590rBhwyR/f39Jo9FIPXv2lF544QVJp9O15KW2Clf07/z586WwsDBJqVRKXbt2lZYsWSKJomi3D9+/7uvfjvz+ratvAUjLly+37VNVVSXNnDlTCgwMlLRarXTDDTdIeXl5dsdx5P7g999/lwYMGCCpVCqpS5cudufwVC3Vvy+99JKUkJAgaTQaKSgoSBo1apS0efPmlrrMVuOq/h05cmSdx0lPT7ftc/DgQemSSy6R1Gq1FBUVJb344ostdJVtlyBJHbw2MBEREREReZw2N0eHiIiIiIiouRjoEBERERGRx2GgQ0REREREHoeBDhEREREReRwGOkRERERE5HEY6BARERERkcdhoENERERERB6HgQ4REREREXkcBjpERERERORxGOgQEREREZHHYaBDREREREQeh4EOERERERF5nP8H0sLbdQJcHM0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(10,4))\n", "\n", "# Plot the results\n", "df['lff'].plot(ax=ax, style='k.', label='Observations')\n", "predict.predicted_mean.plot(ax=ax, label='One-step-ahead Prediction')\n", "predict_ci = predict.conf_int(alpha=0.05)\n", "predict_index = np.arange(len(predict_ci))\n", "ax.fill_between(predict_index[2:], predict_ci.iloc[2:, 0], predict_ci.iloc[2:, 1], alpha=0.1)\n", "\n", "forecast.predicted_mean.plot(ax=ax, style='r', label='Forecast')\n", "forecast_ci = forecast.conf_int()\n", "forecast_index = np.arange(len(predict_ci), len(predict_ci) + len(forecast_ci))\n", "ax.fill_between(forecast_index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], alpha=0.1)\n", "\n", "# Cleanup the image\n", "ax.set_ylim((4, 8));\n", "legend = ax.legend(loc='lower left');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### References\n", "\n", " Commandeur, Jacques J. F., and Siem Jan Koopman. 2007.\n", " An Introduction to State Space Time Series Analysis.\n", " Oxford ; New York: Oxford University Press.\n", "\n", " Durbin, James, and Siem Jan Koopman. 2012.\n", " Time Series Analysis by State Space Methods: Second Edition.\n", " Oxford University Press." ] } ], "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.10.8" } }, "nbformat": 4, "nbformat_minor": 4 }