{ "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": "2023-12-14T14:43:51.552351Z", "iopub.status.busy": "2023-12-14T14:43:51.552106Z", "iopub.status.idle": "2023-12-14T14:43:54.645457Z", "shell.execute_reply": "2023-12-14T14:43:54.644713Z" }, "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": "2023-12-14T14:43:54.649569Z", "iopub.status.busy": "2023-12-14T14:43:54.649072Z", "iopub.status.idle": "2023-12-14T14:43:54.661821Z", "shell.execute_reply": "2023-12-14T14:43:54.661104Z" }, "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": "2023-12-14T14:43:54.665375Z", "iopub.status.busy": "2023-12-14T14:43:54.665059Z", "iopub.status.idle": "2023-12-14T14:43:57.673580Z", "shell.execute_reply": "2023-12-14T14:43:57.672721Z" }, "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": "2023-12-14T14:43:57.677796Z", "iopub.status.busy": "2023-12-14T14:43:57.677168Z", "iopub.status.idle": "2023-12-14T14:43:57.778219Z", "shell.execute_reply": "2023-12-14T14:43:57.777449Z" }, "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: Thu, 14 Dec 2023 AIC -49.020\n", "Time: 14:43:57 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.405e-11 0.000 1.6e-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": "2023-12-14T14:43:57.783601Z", "iopub.status.busy": "2023-12-14T14:43:57.782313Z", "iopub.status.idle": "2023-12-14T14:43:57.797797Z", "shell.execute_reply": "2023-12-14T14:43:57.796988Z" }, "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": "2023-12-14T14:43:57.805179Z", "iopub.status.busy": "2023-12-14T14:43:57.803244Z", "iopub.status.idle": "2023-12-14T14:43:58.593644Z", "shell.execute_reply": "2023-12-14T14:43:58.592779Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAFlCAYAAAAj9p2/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACIuUlEQVR4nO3dd3gUVdsG8Hu2Z9MT0hsh9F7FgAooGhAL+iqKKCJYgVewIPLZwIYNe8UCFhQLgr6KIkWKgAgICEgJJCE9lJCebTPz/bHZJUvabrKbsrl/17XGzM7OnDnZhHn2Oec5gizLMoiIiIiIiLyIoqUbQERERERE5G4MdIiIiIiIyOsw0CEiIiIiIq/DQIeIiIiIiLwOAx0iIiIiIvI6DHSIiIiIiMjrMNAhIiIiIiKvw0CHiIiIiIi8DgMdIiIiIiLyOgx0iIiIiIjI67gU6IiiiCeeeAKJiYnw8fFBUlISnnnmGciyXO/rNm7ciIEDB0Kr1aJz585YunRpU9pMRERERERUL5UrO7/44ot477338Omnn6JXr17YtWsX7rjjDgQGBuL++++v9TXp6ekYN24c7r33Xixbtgzr16/HnXfeiaioKKSkpLjlIoiIiIiIiKoT5IbSMdVcddVViIiIwMcff2zf9p///Ac+Pj744osvan3N3Llz8fPPP+PAgQP2bTfffDOKiorw66+/NqHpREREREREtXMpozNs2DAsXrwYR48eRdeuXbFv3z788ccfePXVV+t8zfbt2zF69GiHbSkpKZg9e3adrzEajTAajfbvJUlCYWEhQkNDIQiCK00mIiIiIiIvIssySktLER0dDYWi7pk4LgU6jz76KEpKStC9e3colUqIoojnnnsOkyZNqvM1+fn5iIiIcNgWERGBkpISVFZWwsfHp8ZrFi5ciAULFrjSNCIiIiIiakeysrIQGxtb5/MuBTrffPMNli1bhi+//BK9evXC3r17MXv2bERHR+P2229vcmNt5s2bhwcffND+fXFxMeLj45GVlYWAgAC3nYeIiIiIiNqWkpISxMXFwd/fv979XAp05syZg0cffRQ333wzAKBPnz44ceIEFi5cWGegExkZiYKCAodtBQUFCAgIqDWbAwBarRZarbbG9oCAAAY6RERERETU4JQWl8pLV1RU1BgHp1QqIUlSna9JTk7G+vXrHbatXbsWycnJrpyaiIiIiIjIaS4FOldffTWee+45/Pzzz8jIyMDKlSvx6quv4rrrrrPvM2/ePEyePNn+/b333ou0tDQ88sgjOHz4MN5991188803eOCBB9x3FURERERERNW4NHTtrbfewhNPPIHp06fj5MmTiI6Oxj333IMnn3zSvk9eXh4yMzPt3ycmJuLnn3/GAw88gDfeeAOxsbH46KOPuIYOERERERF5jEvr6LSUkpISBAYGori4mHN0iIiIiIjaMWdjA5eGrhEREREREbUFDHSIiIiIiMjrMNAhIiIiIiKvw0CHiIiIiIi8DgMdIiIiIiLyOgx0iIiIiIjI6zDQISIiIiIir8NAh4iIiIiIvA4DHSIiIiIi8joMdIiIiIiIyOsw0CEiIiIiIq/DQIeIiIiIiLwOAx0iIiIiIvI6DHSIiIiIiMjrMNAhIiIiIiKvw0CHiIiIiIi8DgMdIiIiIiLyOgx0iIiIiIjI6zDQISIiIiIir8NAh4iIiIiIvA4DHSIiIiIi8joMdIiIiIiIyOsw0CEiIiIiIq/DQIeIiIiIiLwOAx0iIiIiIvI6DHSIiIiIiMjrMNAhIiIiIiKvw0CHiIiIiIi8jkuBTseOHSEIQo3HjBkzat1/6dKlNfbV6XRuaTgREREREVFdVK7svHPnToiiaP/+wIEDuPzyy3HjjTfW+ZqAgAAcOXLE/r0gCI1oJhERERERkfNcCnTCwsIcvn/hhReQlJSEESNG1PkaQRAQGRnZuNYRERERERE1QqPn6JhMJnzxxReYOnVqvVmasrIyJCQkIC4uDtdeey0OHjzY4LGNRiNKSkocHkRERERERM5qdKCzatUqFBUVYcqUKXXu061bN3zyySf44Ycf8MUXX0CSJAwbNgzZ2dn1HnvhwoUIDAy0P+Li4hrbTCIiIiIiaocEWZblxrwwJSUFGo0G//vf/5x+jdlsRo8ePTBx4kQ888wzde5nNBphNBrt35eUlCAuLg7FxcUICAhoTHOJiIiIiMgLlJSUIDAwsMHYwKU5OjYnTpzAunXr8P3337v0OrVajQEDBuDYsWP17qfVaqHVahvTNCIiIiIiosYNXVuyZAnCw8Mxbtw4l14niiL279+PqKioxpyWiIiIiIjIKS4HOpIkYcmSJbj99tuhUjkmhCZPnox58+bZv3/66afx22+/IS0tDX///TduvfVWnDhxAnfeeWfTW05ERERERFQHl4eurVu3DpmZmZg6dWqN5zIzM6FQnIudzp49i7vuugv5+fkIDg7GoEGDsG3bNvTs2bNprSYiIiIiIqpHo4sRNCdnJxwREREREZF3czY2aHR56fZKlmWcLTeh3GiBRZRaujlERERERFSLRlVda+9MogSTKKHMCCgEAVq1AhqlAlqVot7FU4mIiIiIqHkw0GkiSZZRaRJRCRECALVSYQ98VEomzIiIiIiIWgIDHTeScS7bAwBKhQCNyprp0SiZ7SEiIiIiai4MdDxIlKqyPSZrtkejUlQFPkooFQx6iIiIiIg8hYFOM5EBGC0SjBYJpbBApRDgq1VBp1a2dNOIiIiIiLwOJ5G0EIsko7jSjNNlRlSaRLSBKt9ERERERG0GA50WJkoySgxmnC4zocJkYcBDREREROQGDHRaCUmWUWqw4FSZEeVGBjxERERERE3BOTqtjCwDZUYLyk0W6DUq6NVKKFi4gIiIiIjIJQx0WilZBsqNFlQYLfDRKOGrUTHgISIiIiJyEgOdVk4GUFFVolpXFfCwNDURERERUf0Y6LQRMoBKkwiDSYRWrYSvRgmVklOsiIiIiIhqw0CnjZEBGMwiDGYROpUSfjpmeIiIiIiIzseUQBtmsIgoqjCxQhsRERER0XkY6LRxFklGSaWlpZtBRERERNSqMNDxAgaLiHIjgx0iIiIiIhsGOl6izGiBwSy2dDOIiIiIiFoFBjpepKTSDIsotXQziIiIiIhaHAMdLyIDKKo0Q5JYnICIiIiI2jcGOi4yWSQUVZggtdJKZ6Iko8RgbulmEBERERG1KK6j44If9+XiiVUHUFxphkIAgvQahOg1CPZVI1ivQYivBsG+tWzTa+CjUTZbO40WCaUGM/x16mY7JxERERFRa8JAxwnFlWY89cMBrNqba98myUBhuQmF5SbgVMPH8FErHYKfEF8NLu7SAcOSOnhkwc8Kkwi1UgGduvkCLCIiIiKi1oKBTgO2Hz+Dh77Zi9xiAxQCMGNUZ9wwKBZlRos90DlbYcLZcjPOVtT+vdEiodIsorJIRG6RwX7sH/bmIjJAh+sGxODqflEI9dO6te0llWYoFQLUSo5QJCIiIqL2RZDlVjrZpJqSkhIEBgaiuLgYAQEBzXJOo0XEq78dxeItaZBlID5Ej9du6o9wRRn+2vcvOiUlITomtsHjyLKMSrOIs+VmFNoCoXIT0k+XY/WBPPtinyqFgJHdwnDDoFj0jwuCILgny6MQBIT6aqDwQNaIiIiIiKi5ORsbtJtAJzs7G6mpqejSpQtiY+sPUI7kl2L213txKK8EAHDT4Dg8cXVPfP3Fp7j77rshSRIUCgVeeeMd3DJ5SqPaAwAGs4gNh09ixd/ZOJBTYt/eqYMvrh8Yg7G9o+Cna3rSTa1UIFivdlvwRERERETUUhjoVPPxxx87BCiLFy/GtGnTauwnSTKWbMvAi78ehskiIcRXg4XX90FKr0hkZ2cjISEBknRunRqlUomd+w87ldlpyNGCUqzYnY1fD+bDYLaew0etREqvCPxnUCy6Rvg36fg+GiUCWJyAiIiIiNo4Z2MDlyZvdOzYEYIg1HjMmDGjztd8++236N69O3Q6Hfr06YPVq1e7csomy87Otgc5ACBJEu655x5kZ2c77JdXXInbPtmBZ376FyaLhFHdwvDr7IuR0isSAJCamuoQ5ACAKIpIT0tzSzu7Rvhj3pU98PN/L8bDV3RFYgdfVJpFrNqbi9s+/gt3froLq/fnwWgRG3X8SpOISlPjXktERERE1Na4NC5q586dEMVzN8sHDhzA5ZdfjhtvvLHW/bdt24aJEydi4cKFuOqqq/Dll19i/Pjx+Pvvv9G7d++mtdxJdQUox44dsw9h++mfXDy20lo2WqdW4LFxPXHr0HiHoV5dunSBQqGokdFJ7NTJre3106lw4+A43DAoFnsyi7Di72z8fuQU9ucUY39OMV5fl4qr+0XhugExiA3Wu3TsUoO1OIFGxeIEREREROTdmjR0bfbs2fjpp5+Qmppa6/yPm266CeXl5fjpp5/s2y688EL0798f77//vtPnacrQtbqGnGVkZCCgQwSe+uEgVu7JAQD0jQ3Eazf1R1KYX63H+vjjj3HPPfdAFEUolUq8/PrbjZ6jk5uTjbTjx50qanCmzIgf9+Vi5Z4cFJQY7dsv7BSCiRfE48JOoU6fVxCAUF+tR0paExERERF5mkeGrlVnMpnwxRdfYOrUqXVOct++fTtGjx7tsC0lJQXbt2+v99hGoxElJSUOj8aKjY3F4sWLoVRa15NRKpX44IMPkGPywdjXt2DlnhwoBOC/l3bGivuG1RnkAMC0adOQnp6OFT+twc79hxsd5Hz52VIM7t0NN1w9BoN7d8OXny2td/9QPy3uGJ6IldOH4+Ub+iK5UygEAH+mFWLW8r34377cel9fnSwDRRUmtIGpWUREREREjdboQGfVqlUoKirClClT6twnPz8fERERDtsiIiKQn59f77EXLlyIwMBA+yMuLg4AkJOT06i2Tps2DRkZGfj9999x9HgaTkYPw80f/omcokrEh+jx7b3JeOiKbk6tNxMbG4vhF1/S6AIEuTnZeHjWDIc5Q3Nmz0RuTnYDrwSUCgGXdA3D6zf3x4r7huGqvlEAgBd+OYw9mWedboNFku1lrYmIiIiIvFGjA52PP/4YY8eORXR0tDvbAwCYN28eiouL7Y+srCwAQK9evfDxxx836pixsbGI7jEIM1edwAebrGvjTBgci9WzLsaghBB3Nr9eacePu6WoQUywDx4b1wOXdQ+HRZIxd8V+5JytdPr1BouIcqN7gh1ZlmERJZhFqeGdiYiIiIiaQaMWaTlx4gTWrVuH77//vt79IiMjUVBQ4LCtoKAAkZGR9b5Oq9VCq9XW2C7LMu655x6kpKQ0uBbO+TYdPYW7PtsFk0VCsF6NF/7T115RrTl1SkpyW1EDhSDgyat7IqeoEofzS/HQt/vw0e2D4ad17sdaZrRApRSgVSmd2l+UZFgkCaIk2x+Wqq82WpUCfloVVE5kx4iIiIiIPKVRd6NLlixBeHg4xo0bV+9+ycnJWL9+vcO2tWvXIjk5uTGnBXCuYpqrlm5Nh8kiYVhSKNbMvqRFghwAiI6JxStvvOMwZ+jl199u9FA4nVqJl2/siw5+GqSfLscTqw44BB4NKa40w1ItEyPLMsyiBINZRJnRguIKM86UGXGyxIDTZUYUVZhRarCgwiTCaJFqnMtokXCm3IQSgxmSC+0gIiIiInInlzM6kiRhyZIluP3226FSOb588uTJiImJwcKFCwEAs2bNwogRI7Bo0SKMGzcOy5cvx65du7B48eJGN1ipVKJz584uv66iag2ZW4bGIzxA1+jzu8Mtk6dg5GWjkZ6WhsROnZq84Gi4vw6v3NgP93y+G9uOn8HbG45h1uguTr1WloGzFWaoFAIskgzJTUUKKk0iDCYReq0KvhplnQUriIiIiIg8weWMzrp165CZmYmpU6fWeC4zMxN5eXn274cNG4Yvv/wSixcvRr9+/fDdd99h1apVjV5DR6FQ4IMPPnB52BoAGCzWrIWzw7Q8LTqmaUUNztcjKgBPXtUTAPDlX5n4Ya/zhRskWYZJlNwW5NjIAMqNFpwuM3GxUiIiIiJqVk1aR6e52Gpl//vvv+jRo0ejjjHm9c04nF+Kz6ddgIu7hDW6LbIs42SpseEdW8hHW9Lw4ZZ0KBUC3po4AIMSglu6SXZKhQA/rQo6desINomIiIio7fH4OjotISYmptGvNbayjI6nTLsoEZf3jIAoyXj0+3+QfbaipZtkJ0oyiivNOFtuYoU2IiIiIvKoNhXoNIXRbB06pVN79yULgoDHx/VAz6gAlFRa8NA3+1BmaF1r5phECYXlJhRXmF0qnEBERERE5CzvvuuvprXN0fEknVqJl27oizB/LTLOVOCxVfthkVpfBsVgEXGmzIhSVmgjIiIiIjdrN4FOa8noCAKgUgjQKBXwZB2yMH8tXrmxL7QqBf5MK8Sb610vyd0cZFgr4p0uN6LcaEEbmDJGRERERG1A+wl0miGjI8A64V6jVECnVsJXq0KATo0gvRqhvhqE+2sR7q9DqJ8Wwb4ahPlrEeijhlblmaCne2QA5l/TCwDw9c4sfP93tgfO4h6ybF3A9HSZCQYzK7QRERERUdO4vI5OW2QRJViqhkZpVU2P7bQqBRQKAUpBgFIhQCEIUCkEKBSuhSuCIECnVkKnVkKSZBgt1oU6TW6cqH9p93DcO6IT3t+Uhld+O4r4ED0Gdwxx2Cc3Jxtpx4+jU1KS28pdN5YkWwsWVJpE+OtUUCnbTSxORERERG7ULu4ibdkcAE0ubSwIAoL0GgTo1PCtKpWsqQp8mkKhEOCjUVozPX5aBOjU0LjpJn/KsI5I6WWtxDbv+/3ILDxXie3Lz5ZicO9uuOHqMRjcuxu+/GypW87ZVLaCBWUczkZEREREjdAuAp3qQ6HckdHxtPODHn+dCuomBD2CIOD/ruyBXtEBKDFY8PA3+1BqMCM3JxsPz5oBqapQgSRJmDN7JnJzWscQt+oLjhotHM5GRERERM5r/Xf9bmDL6GiUTc+8NDeFQoBeo0KIrwYd/LTw06qgasQ16NRKvHxDX0QEaHGisAL/t/IAUo8dswc5NqIoIj0tzV3NdwtJllFUYUZRhYnlqImIiIjIKe0q0GkL2Zz6KBUCfLUqhPpp0cFPi+JT+di6ZZPTGZhQPy1eubEfdGoF/kovxMYz/lAoHPtEqVQisVMnl9plkSSkny7HP9lFHi1jbbRIOFPG6mxERERE1LB2UYzANnRN60WLhS5d8gnuvvtuSJIEhUKBV954B7dMntLg67pG+OPpa3rjkRX/YPXhIkx+bik+f/wOiKIIpVKJl19/u86CBJIsI7/YgLRT5Th+qsz+NeNMOcyiNfDoGKrHfSOTMKJrGATB/dkzGdbqbJVma7GC9rAuEhERERG5TpDbwEfjJSUlCAwMRHFxMQICAlx+/d6sIox/Zytignyw9dFLPdDC5pWdnY2EhASHYWdKpRI79x92umrap9sy8O7G41AKAp4YHYMAQwESO3VCdEwsZFlGYbnpXEBzuhzHTpYh/XQ5Kky1z5XxUSuhUADlRuvzvWMCMGNkZwxMCG76BddDp1LCX6dqc0MSiYiIiKhxnI0NmNFpg1JTU+ucW+NsoDM5OQHpp8vxy4F8LNqSj7su7oTtB8px/PfdSDtVjqJKc62vUykEdOzgi6QwX3QK80NSmC+SwvwQGahDhVHEF3+ewFc7M3EgpwT3Lfsbw5JCMX1UErqE+zf5umtjsIgwlonw06mg17SLtzNRs7GIEsqNImTI8NEomUElIqI2pV3cGdrm6Oi85B/pLl26QKFQ1Mjo9OjWxeljCIKAeVd2R/bZSuzPKcara486PK8QgNhgPTpVBTK2wCYu2KfOtW38dCrcOzIJNwyOxSd/pGPV3lxsO34G24+fQUrvSNxzSSdEB/k07qLrIQMoNViq1t5RQ9PG52IRtTRJklFmssBgEmFL+RstElQKC3y1Kusixx4YmkpERORO7SPQ8bKMTmxsLBYvXox77rnHPrfmgw8+QI/OiSiqMDmsG1QfrUqJF//TB8/+fAgyYM/OdArzRcdQX5fWHDp/0dFHxnTHxAvi8f6m41h36CR+PZCPdf8W4PqBMZg6PBHBvppGXn3dLJKMsxUm6NRK+Gs5nI3IVbIso9wkosJoQW1jmi2SdUFfhSBAr1FWDVnl7xkREbVO7WKOzo/7cnH/V3twYacQLL872QMtbBnZ2dk4duwYOnfujNhY65A12/waSzOWYf7ys6X29XhqK4xwKK8E7248jr/SCwEAeo0Sk4bGY+IF8fDVeibWFgQgQKdu8gKxRO2BLMuoNItVC/Q6/zoBgE6jhF6trDPTS0RE5G7OxgbtItD5dlcW5nz3D0Z2C8PSOy7wQAtbF1GyBjtSM/xoc3OyMbh3N6cKI/yVXoh3fj+Gw/mlAIBgvRpThyfiuoExTVoQtT46tRIBOhWH2RDVwVAV4DR1jSqtSgG9RsWho0RE5HHOxgbt4l8kg5eso+MspUJAkF6N5ri1Tzt+3OlFRy9IDMGSO4bgufG9ERvsg7MVZixaexQ3ffAn1hzM90hgZjCLOFNugln03Po+RG2R0SLiTJkRxZVmtyzEa7RIOFthwpkyo70ADBERUUtqF3f+tjk67WkYk1qpQKBe7fHzdEpKcmnRUYUgYHTPCHx994WYO6YbQn01yCmqxJM/HMTtn/yF7cfPuH0xUFGScbbchHKjxa3HJWqLzKKEs+UmFFWYPTLE1TaP51SpdXFfqRmH0RIREVXXPgKddpbRsdFWrTHjSdExsXjljXegVFqDyIYWHbVRKRW4fmAsVtw3DPeNSIKvVomjBWWY/fVezPxyD9JOlbm1nbaFRgvLTW759JqorRElGcUVZhSWm2BqhgynJMsoM1pwusyIEoMZFmZViYiombWrqmvtKaNjo9eoIEpynQt9usMtk6dg5GWjkZ6WZl901Fk+GiWmDO+I6wbEYOm2DHy7Owu7TpzFrR/9hZuGxGHaxYnwc2PBArMo4Uy5kYUKqN2orVR0c5IBVJpEVJpECAKgUiigFAQolQJUCgFKhfUr59EREZG7tY9Ap51mdGz8dWqIkux02enGiI6JdSnAOV+gXo1Zo7vgxsGxeH1dKjYdPYUv/8rEmoP5+O9lnTGmV6TbboRkGSiuNMNolhDgw0IF5J0yM7Nw4NARRMd3RFQTfjfdSZatHzaYAeC8kaQKoSrwcQiAFFCyfDURETVSu7jzt02Mbc+regf6qD1W2cydooN88NINffH6Tf0RF+KDM+UmzP/xX9z7xd9IPVnq1nMZLCJOl5lg8mAASNTcJEnGO+8vRmJiR4wbczkG9e6GLz9b2tLNapAkyzCJEipNIkoNFhRVmHG6zIiTJQZr0YQKM8qMFhjMIoefEhGRU1r/na8b2DIZOi9ZMLQxBEFAkI8aijaSvUhOCsWXd16I+0YmQadWYG9WEW7/eCcW/XYEpQaz284jydZFRss8UKhAlmWYLBLKjRbOTyCPEyUZJQYz/jlyHPfPuM9eDVGSJMyZPRO5Odkt3MLGkWEtcGCwiCg3WlBcaQ2ATpVag59yowUmi+T2IiZERNT2tYs7f2Z0rBQKAcF6NVpzrJObk40/Nm9Cbk42NCoFpgzriK/vTsaobmEQZRnf7MrGje9vx8//5Lm1HHV5VaGCpgQkkiTDYBZRarBO+D5VarQHUYWs+kYeYhElFFeacabMiEqTiOMulHxvyyTZGvyUGS04W2H9fTtTVfjAYBb54QIREbWPQIcZnXNUSgUCfZpnjR1XffnZUgzu3Q03XD0Gg6sNt4kM1OGF//TFmxP7IyFEj7MVZjz907+45/PdOFrgvuFsZlFCYbkJlU4WbhCrAhv7J8xVa5JUmESYRclh4nf1qm+8ASN3MIsSiipMOFNugsF8rtCAqyXfvYUt81Npsv5Onik34WSpAUVVHzYYLSKzPkRE7Uy7uPM/V4ygfWd0bKxlp927xo4ANGlYXG5ONh6eNaPe4TZDE0Ox7K6hmDmqM3zUSvyTXYzbP/kLL685gpJK9wxnkwGUGMworjDXWP/DLEqoMFlQXGFdI+R0VWDjypwBWzBVYWJ2hxrHZLGug1NYbqq1wEhjS76fT5JlrPu3ALd+tAM3L/4TX+/ManNZSVm2/v0vN1rn/Jys9nvL9X2IiLyfy4FOTk4Obr31VoSGhsLHxwd9+vTBrl276tx/48aNEAShxiM/P79JDXeFfegaMzp2Phol9JrGBX6CAGiUCug1SgT6qBHqq0F4gA5h/loE+qgbVSUpzcnhNmqlArclJ+Drey7E6B7hkGTgu93W4Ww/7st123A2g0XEmXITSg1mnC034WSJAYXlJpQaLDBYxCadRwZQarDgLNf08TqVJhHFFWaUGsyoNIkwWSS33VAbzCIKy004W9HwOji3TJ6CnfsPY8VPa7Bz/2HcMnmKS+falVGIqUt34rFVB5B6sgzpp8vx6tqjuObtrXhzfSryiw1NuJKWZcvEFlbw94+IyNu5VF767NmzGD58OEaNGoVffvkFYWFhSE1NRXBwcIOvPXLkCAICAuzfh4eHu97aRmJGp3b+OjUkyXpTXxeFIECtFKBSKqBSCFAr6y/3qlMroVMrUWmyjp13NiCwDbepHuzUN9wmIkCH567rg+sGFOLlNUeQcaYCz/18CKv25GBOSjf0iAqo9XWukGTPrj9kEiWcKTPCT6eCXtMuKr3bybK13Lkky9CplFC08RLCBrP1/e5443zuvWNfP0bhuHaMyolKiAazdRK+xcWb8saUfD9aUIp3fz+O7WlnAAB6jRKThsYjWK/B8p1ZyCyswLIdmVj+VxYu7RGOWy6IR8/opv+utQRRklFYbkKQvm1UpCQiIte5dHf14osvIi4uDkuWLLFvS0xMdOq14eHhCAoKcqlx7sKMTt0CfFQQK2SYRQlKhQC1QgGVUoBKaf3/xt6A+miU0KkVqDCJKDdZ0FC8YxtuM2f2TIii6PRwm8EdQ7DszqH4elcWPtqSjoO5JbhjyU6MHxCD+0YmIdDHvUP03M2W3bGu6dO4bFhbYRYlmCwSjBYJlmpzmMpggVZtzTC2tRtOZ4MQ+/ox58XNAqxFQs5fN0apEKwV+0znB0+ekVdciQ82peHXA/mQASgVAq4fEIOpFyUixFcDALhuYAy2HTuDL//KxO4TZ7H23wKs/bcA/WIDccvQeFzcJazNvX9tVReDfDTQtNN11oiIvJkguzA7s2fPnkhJSUF2djY2bdqEmJgYTJ8+HXfddVedr9m4cSNGjRqFhIQEGI1G9O7dG/Pnz8fw4cOdbmRJSQkCAwNRXFzskBVy1uhXN+HYyTJ8ddeFSE4Kdfn13s72FvDUwpmyLKPcJKLCiYAnNycb6WlpSOzUyeVPo0+VGvHWhlSsOVgAAAjz0+LJq3vigsSQxja9WQmwZtl8GjmksLWRJOu6KEaLNcBxJrunrhoSqVUpWvVCrkaLiHKjteiEJ+TmZCPt+HF0Skpq0kK8DSmuMGPptgx8uzsLZtH68xndIxz3jkhCXIi+ztcdLSjFV39l4reDBfYgLybIBzcNicNVfaPgq21bGUoBQICPGjq1d/zuERF5O2djA5cCHZ1OBwB48MEHceONN2Lnzp2YNWsW3n//fdx+++21vubIkSPYuHEjBg8eDKPRiI8++giff/45duzYgYEDB9b6GqPRCKPR6HAxcXFxjQ50Ln5pA7IKK/H99GEYGN/wMDvyDEmSUW6yoNIkwpOfUf994iwW/nIYmYUVAICbh8Rh+qgkl4cuNtfN5vk0SkWbze6YLBJMVZmbpgQBggDoNSr4qJWtqh/MooQyg6XBOTJN8eVnS+2FORQKBV554x2X59g0xGAW8fXOLHy2/YR9DanBCcGYeWnnGsM+6/s9OFVqxHe7s/H9nmyUVFqP46dVYfyAaEwYHIeIAJ1b2+1p/u1wGCkRUVvkkUBHo9Fg8ODB2LZtm33b/fffj507d2L79u1ON27EiBGIj4/H559/Xuvz8+fPx4IFC2psb2ygM+S5dThVasTP91+EXtGBLr+e3EuSZJSZLDB4MOCpNIl4a0MqVvydAwDo1MEXC67tha4R/k69vjluNusjCIC/tvVndyRJtmdsjKLYYMauMXQqJXw0yhYdWmQRpaoSxZ4tDZ6bk43BvbvVmKu2c/9htwTbFknC6n/ysXhLGk6VWj9M6hLuh5mXdsbQxJAaWTRnfw8qTSJW78+zz+MBAKUgtMl5PL5aFfzaWEaKiKi98Uigk5CQgMsvvxwfffSRfdt7772HZ599Fjk5OU43bs6cOfjjjz/qDI7cndHpO38NSgwWrH9oBJLC/Fx+PXmGKMkoM1rsc6g84Y9jp/HsT//ibIUZaqWAe0ck4Zah8fWWwvb0zaYrtCoF/HUtn92RZRkWSYYoyZBk61eTRXJ5gnxTqBQC9BoVdOrmG9ZmESWUG8V6C3a40x+bN+GGq8fU2L7ipzUYfvEljT6uLMvYknoa7248jvTT5QCAqEAd7hnRCSm9Imv9fWjM74Ekyw7zeGz6xQYipbMvAg0F6Nmtc7P/HrnKR6NEgJtL8BMRkfs4G+i49LHV8OHDceTIEYdtR48eRUJCgkuN27t3L6Kioup8XqvVQqvVunTM+hjsVdc42bQ1USoEBPqo4atReuxm8qLOHfDVXRfi+V8OYfPR03hrwzFsPXYa86/pVeewmvpKXTf3DZrRIsFUbkSAzrPzByRJhiifC2QskmzdVrW9NayzaJFklBjMKDUCPmol9BqVxwJAsWqYpSezjrVxtfqgM/7JLsLbG45hX3YxAGsBkqnDE/GfgbH1Zska83ugEARc1KUDLurSwWEez77sYvv5xXXbkBDqiyE9k5AQqkd8iPURFaSDStE6/kZXmkRIkmxdXLkVzxUjIqL6uRToPPDAAxg2bBief/55TJgwAX/99RcWL16MxYsX2/eZN28ecnJy8NlnnwEAXn/9dSQmJqJXr14wGAz46KOPsGHDBvz222/uvZI6yLL1k2cAnGjaSqmUCgTqFdCLSpS7aXjQ+fMKXvpPX/y4LxevrU3F35lFuOXDHXhkTDek9Iqs8VpP3Gw2hSwDxZVmVJhEKARAQNWNV7X7L+G8Tbabs+q3aLZ9ZBkQ5WqBjCQ36818U8kyUGESUWESoVUp4KNRQq1QnOuDRt6YZmdn48iRo4hJSERIeFSL9Eljqw/W5vjJMnywOQ2bjp4CYP2g5+YL4jD5wo7w0zX8p7+pvwddI/zx1NW98J/uvrjxwYXQ974UKv8OUPqFItsIZO9xHAWgUgiIDfZBQqivQwCUEKpHkF7jwpW7h9Ei4WyFGUE+6jZf/pyIqL1yKdAZMmQIVq5ciXnz5uHpp59GYmIiXn/9dUyaNMm+T15eHjIzM+3fm0wmPPTQQ8jJyYFer0ffvn2xbt06jBo1yn1XUY/qN83M6LRuaqUCQXoNjBYRpYbGl9Wta17Btf1jMDA+GE/9eBAHc0vw5A8H8UfqaTwyphv8qw1TcefNpjt5qsJXY7RUoYbzGavKVZ9PqPqPAAGC4BgACkDVNusTggB8tnQJZs24r8XmZFV3y+QpGHnZ6EZXHzyQU4yl2zKwJfU0AEAhAFf3i8adFyci3N/54gDu+j0oyjuBs5s/w9nNn0HQ6KEOjYU6JAYT75sDybcDThRWIKuwAkaLhIwzFcg4U1HjGAE+KiSE+CI+VI/EDr7oFRWAHlEBHp/DZhYla/lpvabFh48SEZHrXJqj01KaUl66uMKMfk9bs0epz41tc+t0tGcVJgvKjA2XpK7OmXkFFknC0q0Z+OSPDIiyjIgALZ68qicGdwypcazG3mzayLKMghIjgn3VXrNgbUsXanC31jQnq7FkWcaujLNYui0Du6rmxggALu0ejrsu6YTEDr6NPnZTfw+c6V9JllFQYsCJMxXIPFOBE4W2r+UoKDHWelylICAp3Be9ogPROyYAvaMDER+qr3f+XWMpBAHBerVTC7wSEZHneWSOTltkrJr3oRCsQyOo7dBrVNCplCh1oWCBM/MKVAoF7ry4Ey7sFIqnfjyI7LOVmPnlHtwyNB73jkiyz1tozMryoiTj2Mky7M0qsj8Ky00I1qsxd0x3jOoe7tLxWpvcnGx7kAMAkiRhzuyZGHnZ6DYTFJyvNc3JcpVUVWTg020ZOJhbAsA6921s70hMTk5AQmjjAxybxvwenP/6hjJDCkFAVKAPogJ9cGEnx7XOKk0iss5W4MuVq7H8xzVQhXWENrob4N8BRwvKcLSgDCurhsH5aVXoFR2A3jGB1q/RgQjUN72ogCTLKOTCokREbY7XBzoGs60QgZKTStsgRVXBAr1GiVKDpcHhW67MK+gdE4jPp12AN9alYtXeXCzbkYkd6YV4+ppeSAp3rjqfySLhUF4J9lQFNf9kF6HcWDMoO1thxqPf78flPSMw54pubrn5agltOSioizvnZDXXkD6LJGHdvyfx6bYMpFVVUdOqFLi2fzQmDU1AZGDrWr+mKcPxfDRK+FmKsXjeVIefkSYwHK8t/xU5lSoczC3BobwSlBkt2JFeiB3phfb94kJ80Ds60B4AdQn3a1RmRpaBogoTAvXek50lIvJ2Xh/o2DI6OjU/hWvL1EoFQnw1qDSJKDWa6xzO5uq8Ar1GhXlX9sDwzh3w/OpDOHayDFOW7MT0UUm4aUhcjWEw5UYL9ucUW7M1mUX4N6+kxhwRvUaJfrFB6B8fhP5xQegc7ofPt5/A59tPYO2/BdiVUYhHx3bHyG5tL7vT2go1uIO75qI0x5A+k0XCT//k4os/M5FTVAkA8NUqccOgWNw8JB4hvs0/ad9ZTckM1RZgm4pPIlI8hf9cZi27bRElHDtVhoM5JTiQW4wDOSXILKxAVmElsgor8cuBfADWgLBbpD/G94/BuL51V/+sjQzrcOgAHxa3ISJqC7x+js7+7GJc/fYfiAzQ4c//u8xDLaTmJMvW9Xcq6yn925h5BWfKjHhu9SFsPXYGAHBBxxDMHt0FWWcr7MPQjuaXQTzvVyZYr0b/uCAMiA+2Bza1TVw+lFeCBf/7176OSUqvCDx0efNmd9yRcfjys6U1goK2PEfHpilzUXJzsjFk0ED4DboGmg4JMBflQSzKw7uvvoiB3RIQ4qtpUka5wmTBqj25WLbjBE6XmQAAQT5qTLwgHv8ZFONQTMMbNXYeVXGlGf/mluBATjEO5loDoFKDxf787cMScN+IpEb9bPy0KvhyYVEiohbhkQVDW0pTAp1dGYW44f3t6Biqx8Y5zVPpjZqHRZRQarDA5MZqZLIsY+WeHLy+LrXOMtfRQTr0j7NmawbEBSMuxMfpGyWTRcJHf6Th8+0nIMlAiK8Gj47tjhFdw9x2DXVxZ8ahqRPUzaKEfVlFOFVmxLCkDgj0abs36hUmC176dgt+PlICha72IY9+WpW1XHJV2eSEqv+PC9bXWzmsuNKMb3dl4etdWSiptN6gh/trceuFCbi2f7THswoqhQA/nQoCBFSY3FP6vbHcEWBLsoyswgqs3p+PpdsyAABX9Y3CvCu7N2oNH51KCZVSgFIhQCEIUCkElqImImoGDHSqbD12GpM+2oFuEf5Y80DjVxan1stgtpajltz4Vs48U4EFPx3EgZwSdOrgiwHxQehXFdzUtdBoXWrLohzMLcbT//vXXkp3TK9IPHhF1zpv+JuaiWkNlcXyiw3Ydvw0tqedwa6Ms6gwWYeValUKXNErAjcOikO3SP9maYs7mCwSVu7JwZKt6ThbYbZuO5WB8gProfQPgyY0Bgl9LsTJMnO9a/KE+2sd1o2JD9UjzF+LXw/k4/u/c+z9FBfig8nJHTG2d6THq0cqBAH+OlWNQKqlFlK1cUclRJsf9+Zi4S+HIMnWhYWfu663WwJHAda5hbagRylYAyFl1f8zECIiajoGOlXWHyrAtE93oW9sIH6ceZGHWkgtTZZllJtEVBgtbr0BM4tSk24q68uiGC0iPtqSji/+tGZ3QquyO5ecl91xRybmj82bcMPVY2psX/HTGgy/2DMfAJgsEvZmFWF72hlsP37GPmTPJlivRpBe47C9T0wgbhwci0u7h7faUvAWScLq/fn4eEs68ksMAIDYYB/0ErLx8f9NgyhaHDIORouInLOV1tLJhedKJ2cWVqC40tzg+bqE+2HKsI4Y1T3c42u5CALgq1FBr6m/eIssy/ZFW935AUNz23z0FB5fdQBGi4Q+MYFYdGM/jw8lrS0QUisVrOZGROQCBjpVVu/Pw/Rlf+OCjiH45t5kD7WQWgtRklFmsMBgca4ctSc5m0U5kFOMZ36qlt3pHYkHL7dmd9yViWmujE5uUSW2Hz9jz9pUVisLrhCsle6SO4UiOSkU/mIx0o8fhykgGhtPmLD+8En7IrHBejXGD4jBdQNiGsygNVelM0mW8fvhk/hgUxpOFFp/VmH+Wky7KBFX942CSqlwOeNQXGGuCn7K7YFQ5pkK5BZXomuEP24f1hHDk0I9XjFSgLW6ma9G5XLGwWC2BjytaUFbV+zLKsLD3+5DicGCjqF6vDlxgMtZW3dQCAJ8NEroVAqu10NE1ACuo1PFtv6KllXX2gWlQkCgXg2dRYGSSvcOZ3OVs6WYe8cE4rNpF+DDzelYtuMEfj2Qb6/MJuS5p5yzuyqLnc9oEbEnswh/VmVtzl/VPtRXgwuTQpHcKRQXJIbYh+bVlqWaPXMiVu3Nxcq/c3CqzIglWzPw2bYTuKRrB9w4OA4D44Nq3PA3R6UzWZaxPe0M3t+YhiMFpQCAQB81bh+WgP8MjHUY7uRqZbFAvRp99IHoExvo1ja7QqdWwk+ranS2SKdWQqdWwmSRUGkSW8WHDK7oFxeED24bhFnL9yLjTAWmfboLb97cH53CnCsx7y6SLKPcaEG50To3ytavns7iERF5M6/P6Hz1Vybmfb8fo3tE4KPbB3uohdQaybKMEoPzi426W2OyKPtzivFstezOyKQAfDFrHMTKUqeP0VCbGjvHwWAWkVdsQG5RJTILK7AzoxC7T5y1r1UFWFer7x0TgGFJHZCcFIouEX41SnQ31C8WUcKmo6fw3e5s/J1ZZN+nUwdf3DAoFmN6R8JXq2qWLNWezLN4b+Nx7MsuBmAtHT5paDxuviAefm284pZWpYCfVuX27IEoyagwWVBpFussA98aFZQYcP9Xe5BxpgL+OhUW3dgP/eKCWrpZUCsV8FEroVUpOL+HiKgKMzpVjMzotFuCYF1sVKtSoMRQ99o7ntKYLEqfmEB8OvUCfLglDV/uyMTG4yXoNvtznPjueZSn/tXkTEx9GQeTRUJ+iTWQsQU01b8WlptqfV0HPw0u7BSKYUnWrE1DpY4bynSplApc1iMCl/WIwPGTZfhudzZ+OZCPtNPleGnNEbz9+zFc1TcKHaU8jy1eeji/BO9vSsP249ZS41qVAjcMisXk5AQE6VvvWjXOUCutAY6n5oQoFQL8dWr4aVWoNIsoN7aNeTwRATosnjwYD32zD/tzivHfr/bg2fG9a8yZa25mUYJZlCAA0KgU0FUFPVwAm4ioYV6f0Xl/03G88MthXD8wBq9O6O+ZBlKrJ0kySgzmFimP29gsyv7sYjz907/IrJoPMihcwMDEMIQEB0NRvZJTVWnbc1Wd4Fjutuo528TnCpMFuUUG5BZXIs/2tdiA06XGBgs56DVKRAf5ICpQhz4xgdasTbifSzddjcnElBks+Hl/Hr7bnW3vDwAwZOxFyd8/oTJtNyCam5zRyThdjsWb07D+8ElruxQCru0XjTsu6ohw/+aft+FOSoUAP23NSmrNwWAWUVk1j6e1/4NjMIt4bOUB/HHsNBQC8OjY7ri2f0xLN8uBAECrVkKnVkCr4sKlRNT+sBhBlTfWpeK1dUdxy9B4PH9dHw+1kNqKSpOIUkP95X5bE4NZxAeb0/DVjsxmabNOrUBUoA+ig3S1fg3QqdzySXJj10SRZBk7Mwrx7a5sbD12GtJ5naKEBL1OA01VFSv7V5UCWtW5/9corTeIaqVQ9ZwSp0qN+O3ffEiy9UYypVck7rw4EXEh+iZfb0sSBOs6PnpN60jgi5IMiyRBkqwV7Kzfy/ZCFK2BRZLwwi+H8b99eQCAey7phDuGd2yVWRRBsM6T8lErW22lQiKiWtlCEHsoIjv9/yVl5QgMDefQNdvEWC1LdxKslaXUSgElBkubqBKlUysx67IuuLRbOFbuzYHBJEKUZIiy9cZQqvpqe0gy7M9L0rl9LNX21aqUdQYywXp1s9zM3TJ5CkZeNtrlTJdCEDA0MRRDE0ORW1SJlXtysOrvbJQYrb/nIhQoNVia1LZLunbAPZckoXN4805GdzdBAPQaFXwbKBXd3KxZSFsW4lw2Qq56f9qCnnNfpWYfdqpSKPDYlT0Q6qvF0m0Z+GBzGs6Um/Dg5V1bXXEAWbZ+gFNpEqFUCNBrlNCplJzPQ0T1cwgy5JpfnX2uxn6oY1ttr28Ck9Gp3bw+0DFWTZRuieEa1DqplAqE+GpQZrS4fd0dT+kT27KVuTzB1QplNV4f5IMZozrjvpFJqDCJMFkk+8MonvveaLHOcbD9/7l9JIfXAMClPcLRJ6Zt97NCsN7sNrQWTmsjCAJUSgG1jcSSpHPBulmUYBFljw+DEwQB941MQqivBq+uPYrvdmejsNyE+df0bLXDxURJRqnBgjJYoFUpodM0/9A2WZZhtEhQCALXBiJqLFkGZKnuQMOpr2j4+XbA6wMdZnSoLn5aFbQqBYorza1q2Ay5RiFY555A29ItaVkKQYCv1jqEqS0FOM5QKARoqjIU1T+0sgU9pqoJ+87+Hruy9tKEIXEI8dXgqR8PYsPhkyiqMOHlG/rBT9d6//mUYf23z2ARoRAs8NFY3xeeykZZRAkmUYLRLDkEoApBgFatgE6lbLGgxyxaM4ItHXRZxHMftCiq5laqqs+lZAau7ZLlc4EJzg9QagtWJMdAxP669heENIfW+5faTZjRofqolQqE+mpQarSg0uTeMtQqhQB11XwQo1lqc+uLUNvQkkUGWppaqYBaCfhUDYGTZVvQI9tvvs+/Z2jM2kuje0Yg0EeNR1b8g78zi3DvF7vx6IgIFOWd8PhCtU11bn0eCzRKBXw0Ta/aZutn2417XQGmJMv2YXXNFfSIkuyQ2bX9/AUB0Kqs195cVessogSDRYLRLMJSvY9q+adAgDWgtwU91YMgpULwug8vWo3qmRNZOveoK2CpM4Ch1srrixHM/PJv/PRPHp66uifuGJ7ooRaSNzBaxCYtMmq96RKgVlqDm/M/oTOYxRYpc03eSaUQ4NtOAxxXiJJ1mJtJlHAiMwv9e3Rp9NpLh/NLMHv5XpytMMNSVICCb56AVJzvkYVqPUkQAJ+qAgbOrqNkCyCMFuuw0Kb8GXNn0GMbKmcS6w+6qrOV6rYFPu7MptQZ3DRR9cqaSqXQ5tfxcjtZBiSxlkBFOi+QqS3zQm1RSYURgVGJLEZgW8ywtY6pptZDq1Ii1Nc6mb2h7IsA61wfW9UujbLhTwh1aiU0SueOT1QXtVIBX62Sf9OcZCt+oFMrcTrnRJPWXuoeGYDnUmJx1yfboA6OQuStL6Py2F94+ustkDsOQc+O0YgK9EGYv7bVFS2oTpaBCpOICpNoX5BUp675N6x6YOPOm/bqmR5b1ThXgp7qc+4sjZirJQMwVr0esP5O2TI9jVlA1xbcGMyix4ZBS7IMSZRhBiBY4P2BTvUMS/UAps4HAxaqnZf/plg/pQesZXOJGqJQCAjUq6E1Oy4yagtsNKqq4MaJwMaV45N3sJWwFqtVDpMk2S2T5jVKBXw9uNBne9ClSxcoFIoaGZ3ETp2cPkblqSzkfzEH4TfOhzayM/z6Xg4AePOPfOCPfOsxFQIiArSICrSuORUVqLOvP9XaAiHbgqSlBuvaPBqlosawL0+qXjXOFvRoVY5FFOoajuYutj4oM1p/drbz1/e71hzBjdeSxKrg5fyvsvX/GbiQG3l/oMOMDjWCrmpNCoNZtA9Jc+cYaWZ3zo1J94abhIaCEIso2csli7IMUTxXRawhOpUSei3XSHGH2NhYLF68GPfcc499DafX3nrHpTk2nZKSAEMJCpbNhU+XoVAFRUIdFIlhKeNRaJCRX2yARZKti/IWGWo9xvmBUJdwP1zaI7xFF6WVYR1eazC33N8jx6DHDK1SCbPkfJEJdxAl2Z7tOn9ejyjJDG4aUiOIkar+Xzq3jagZeX+gw4wONZKyag6Ep1TP7pQaGj83qC2xfVpafbifySKh0iTCaBHbRKnv6pwdRqZSKmotmyzL560bI1oDIYskQatSwlfj/BwKcs60adOQkpKCY8eOoXPnzoiNjUWJwex0MZLomFi88sY7mDN7JioOba624G0yAOuN8ukyI/KKDcgrrkRekeHc/xcb6gyEXl+XioEJwRjTKxKjuofBX6f2yPW3FbKMFv8QSJbPBX8C0Ob+PnmMJAKSBRDN1q+2wIaZGGqFvL4YwZjXN+Nwfim+mDYUF3Xp4KEWEjWNJMkoNVpa9NNUT7BN+rVN/K1vuI4kyagwWz/Nbe1BX3uudOatyqoqkzkrNyfb5QVvgZqB0OHMAvyZdgYZJef2USsFDE/qgCt6RWB45w58n5EDAUB4QDNk/2oLaCQzgxlqFViMoIrtxlHLjA61YgqFgEAfNbSqtp/dUSqEqsDGuSINNoqq4MFPq4KhKuAxiVLDL2xGtjV7fDS88fQ2floVFAJQanAu2GnsgrfWYWs6RATo8O/vq/ByValrdVAkbpjzCvI1MUg7XY6NR09h49FT8NUqMbJbOFJ6RWBwQohH5/a4sr4QeRFJrBbMVHu04X+HiGy8PtCxVVXRcY4OtQH2uTttKLvjStbGWTq1tUqWRZRQYRZhMLXssDZBsN4Ie+NinHSOXqOCQhBQUmn2+PstNyfbvp4PAJiL8vHNk7dj5/7DKFcH4reDBVhzMB8FJUb8/E8efv4nD6G+GozuGYExvSLRI8rfre/Fxqwv5GmtKfBqTW1xG3MlYChmQENezevTHMzoUFtjy+4E6dVQtMKbagHWNVx8NEoE6dUI89ciSK+BXqNy+6fNKqUCATrrOfx17j9+Q2wBTpifFnqNikFOO6BTKxGk18DTP+q048frLHXdJdwfM0Z1xqoZw/H+rQNx/YAYBPiocKbchK93ZuGOpTtxw/vbsXhzGjLPVDS5LecHXZIkYc7smcjNyW7ysRvry8+WYnDvbrjh6jEY3LsbvvxsKdvibpxTQ+2A18/R6fnkr6gwidg8ZxTiQ/UeaiGRZ7TU3B3bity2FbqV1Vbpdufieo1htIgwmCSPTlQWAPholPDVqFr8eqllmEUJRRVmjw0jzc3JxuDe3ZxevNQsStiRVog1B/OxOfWUfY04AOge6Y8xvSNxRc8IhPppXW7LH5s34Yarx9TYvuKnNRh+8SUuH6+pXO2b9tIWG7fN0TGVA4aShvcjaoU4R6eKbegaMzrUFtmyOzr1eXN3zrv3cvVWTBAAlUJhX2VbKQhQKKq2tfIbe2u5VyX8JBUqzSIqTBa3fSgpANBplPBjgNPuqZUKhPhqcLbC5JFSwtWrt9lKXb/8+tt13jyrlQpc1KUDLurSARUmCzYfPY01B/OxI60Qh/NLcTi/FG+tP4aLunTANf2icWFSCFQK5/7d65SU1OT1hdypvmxXcwcX7myLVw5/I2rlXA50cnJyMHfuXPzyyy+oqKhA586dsWTJEgwePLjO12zcuBEPPvggDh48iLi4ODz++OOYMmVKU9rtFLN4rv6+lovsURumVSmh9XN+ntn5idrq3woCvGIIlq3yma9GCUk+d80yqtadg+xw3bZt5/7f8TUC4JHhd9R2KRUCQvQaFFWaYfZAYYxbJk/ByMtGu1y9Ta9RYUzvSIzpHYmz5SasP3wSvxzIw4GcEmw6egqbjp5CmJ8W4/pG4ep+UYgNrn80g6tBV13MooR9WUXYkV4IX60Kl3YLb9RIitYUeLmrLa1xDhRRe+DS0LWzZ89iwIABGDVqFO677z6EhYUhNTUVSUlJSEpKqvU16enp6N27N+69917ceeedWL9+PWbPno2ff/4ZKSkpTp23sUPXyowW9H5qDQDg8DNjWKKTiIhcJssyiirMra4KYHW5OdnYtv8YjhgCsDm9DEWVZvtzgxKCcW3/aIzsFlbvmk+NKZl9psyIbcfPYOux09iRXoiK89Yj6hzuh8u6h+PS7uHo2MHX6ev58rOlNQIvVwIDWZaRfrocO9ILUWkSkZwUiu6RjSvg0NS27D+ajv/c9yh8ug6HJqITxLJCiMUFGD/mUnSO6YCYIB/EBPkgOkjn1PpJHLpG5PzQNZcCnUcffRRbt27Fli1bnG7I3Llz8fPPP+PAgQP2bTfffDOKiorw66+/OnWMxgY6p8uMGPzsOgBA2vNXcigKERE1WnGluVVWQzw/W/Di6+8g+oKx+HFfLnakFdqHtvrrVBjTKxLX9I9G1wj/Rp1LkmUcyivB1mNnsO34aRzKK3V4PlivRr8oH+QXliK1SIZY7Q4jKcwXl3YPx2U9IpDoRNDjauBVWG7CzoxC7EgrxF/phThVZnR4Psxfi0u6dMAlXcMwKCEYahcW43W1LSdLDdhw6CTWHz6Jf7KLnT5PgE6F6CAfRFcLfmz/Hxmog1qpYKBDBA8FOj179kRKSgqys7OxadMmxMTEYPr06bjrrrvqfM0ll1yCgQMH4vXXX7dvW7JkCWbPno3i4tp/+Y1GI4zGc3+gSkpKEBcX53Kgk1NUieEvbIBGpcDRZ8c6/ToiIqLalBrMNbIWLamhyfL5xQb89E8u/rcvD/klBvs+3SP9cU2/aFzRK6LBLEKpwYy/0gvtwc3ZCrPD8z2i/DE8qQOGd+6APWu/x5zZ1qBLpQ/A7U++jfKQLtiZXghLtblOnTr44rIe1kxPpzC/Rl270SJiX1YxdqSfwV/phThaUObwvFalQM9wHSzGSqQWAwbLufPrNUoMSwrFJV3DMCwp1KlMSkMKSgzYcPgkNtQS3BhzDqH88BYYMvZBoQ+AJjgKd8+ZjxJRjdziSuScrazRr+dTCEC4vw4xQTokhfvh8at6IqAp7WagQ22YR4oRpKWl4b333sODDz6I//u//8POnTtx//33Q6PR4Pbbb6/1Nfn5+YiIiHDYFhERgZKSElRWVsLHx6fGaxYuXIgFCxa40rRa2UtLc34OERG5gb/OWva9zOjcwqKe1tBk+chAHe68uBPuGJ6InRmF+N++XGw8cqqqgMERvLE+FZd2D8e1/aPRPy4IgiDYh31tPX4G246dxr7sYoeCDL5aJYYmhmJ451Akdwq1V3rLzcm2BzkAYKkowdLH7sDO/Yfhd00vbEk9jfWHC7AjrRBpp8uRtiUdH25JR8dQPS7rEYHLuoejU5hvncPLZFnGsVNl+CvdmrXZm1VkLzhk0zXCD0MTQzE0MQSHNv+IR6vao1BrMeO59yHE9cOW1NMoLDdh3aGTWHfoJJQKAQPjg3BJlzBc0jUMkYHOZ0tswc36QyexP8cxuOkbG4jLuodjVPdwrF+VhTlf/Wwf/rZwznTccvUAh/0rTBbkFhmQW1SJnKJK5BZVIrfIYP9/o0VCfokB+SUG7MkqwnPX9XG6nUTtlUsZHY1Gg8GDB2Pbtm32bffffz927tyJ7du31/qarl274o477sC8efPs21avXo1x48ahoqKi1kDHXRmdf3NLcOWbWxDmr8XOx0Y7/ToiIqL6GMxisyws2pDGlD8uqjDhlwP5+HFvLtJOl9u3x4X4oF9sEHafOIu8YoPDazqG6jGscwcMTwpFv7igWod9OVumutRgtgY9h05iR/oZmKuNb+sYqsel3cNxaY9wdA7zQ2G5CTvSrUPR/kovxJlyk8Oxw/y0uKBTCIYmhmBIxxCE+Goa7JfI6BgczC3B5qOnsPnoKWSctxZR1wg/e9DTNcKvRuBVV3AjoCq46RGBUd3DEO7vGDA1Zg6UjSzLKCw32YMeiwRMuyjRpWPUwIwOtWEeyehERUWhZ8+eDtt69OiBFStW1PmayMhIFBQUOGwrKChAQEBArUEOAGi1Wmi1rq8FcD6jhRkdIiJyP51aCZVCgEWSIUoyRFmGJMmwSNavzRUANaZiWpBeg4kXxOPmIXE4mFuCH/flYu2/BcgqrERWYSUAQKNUYFBCMIZ3DsWwpA6ICa793+vqnK1Q5q9T48o+UbiyTxTKDBZsOXYK6w+dxJ9pZ5BxpgKfbM3AJ1szEOKrQeF5gY1OrcDA+GBckGgNbhI71J4BaijT1ScmEH1iAjFjVGdkFlbYg579OcU4WlCGowVl+OiPdEQG6HBxlw4Y1jkUGacrsP5wAQ7knAsOBAD94oLsmZsw/7rvXaJjYhtdVloQBIT6aRHqp0WYUIaz+VnIzlYjNpZlqonq41KgM3z4cBw5csRh29GjR5GQkFDna5KTk7F69WqHbWvXrkVycrIrp24U24JqDHSIiMjdVEoF6ipiJlUFP2IzBEKNLVMtCAJ6xwSid0wgZo/ugnWHTuLEmXL0jwvC4IQQ+Ghcq1TamKDLT6fC2N5RGNs7CmVGC/6oGt725/FCFJabIADoFumPoZ1CMDQxFH1iAqFx4t90V8pCx4foceuFCbj1wgScLTfhj2OnsTn1FHakFSK/xIBvd2fj293Z5/oNQP+4IFzWIxwju9Uf3Ljb+YUnFi9ejGnTpjXb+YnaGpeGru3cuRPDhg3DggULMGHCBPz111+46667sHjxYkyaNAkAMG/ePOTk5OCzzz4DcK689IwZMzB16lRs2LAB999/f7OUl9545CSmLNmJXtEB+Pn+i51+HRERkSdVD4TMogSDWTq3IHAb15QhWjblRguOFpQisYMvgvSaRh2jqWWhDWYRf6UXYnPqKezKOIvIAB0u62HN3HTwa77gxqau4XgZGRmNy+xw6Bq1YR4ZujZkyBCsXLkS8+bNw9NPP43ExES8/vrr9iAHAPLy8pCZmWn/PjExET///DMeeOABvPHGG4iNjcVHH33kdJDTFMzoEBFRa6RQCFBAgFppHQbnr7MOtzaYJRgtItpyzNOUIVo2vloVBsQHN+kYjc102ejUSlzS1TpXpzWoazjesWPHXA90SkuB8iJA7/K68URtisvv8KuuugpXXXVVnc8vXbq0xraRI0diz549rp6qyWxzdLhQKBERtXZalRJalRKyrILRIsFoC3paumFtmDuCrtairuF4nTt3dv1gX34J3HsvEBgIJCZUPTo6fo2OAhT8oJjaNq8O5W1lJ5nRISKitkIQBOjUSujU54KeSpMIkyg1/GLyWrXNgfrggw8aN2zNViSquBjY+4/1cT6tFkiIAzrWEgglxAE6NyxaSuRh3h3o2NfRYUaHiIjanupBjyTJMFQNbzMz6GmXbpk8BaMuG42igmx07ty58VXXnnwSuP8+4PABID0DSD9x7mtGJpCZBRiNwNFj1sf5BAGIjrQGPR2rBUG2oCikacMOidzFuwOdqoyOTs2MDhERtW0KhQC9RgW9BhAlGQaziEqz6LCYJ3m/6JhY9O/RiOFq59PrgZ7drY/ziSKQnVszCErPADJOAKVlQE6e9fFHLesockgctRJeHegYmNEhIiIvpFQI8NWq4KtVVVVtE72qchu1MKXSOjwtIQ4YeV7VWlkGzhSeC3rSTzgGQ/kFzg2JS+wIdIznkDjyKK8OdJjRISIib6dWKqBWKuCvQ1XAI9r//SNyO0EAOoRaH0MG1Xy+osI6/K1GNuhE44fE2b4GB3nyysgLtYtAR8uqa0RE1A7Y5vO0xNA2AYBGpYBWpYQkW89v4bA6B4KANl063Cn1DYmzWKzD3TgkjpqJVwc654au8Y1PRETtR/Whbaaqqm2eKFWtVAhVwY0CGqUCgiDYn/PVquwBV3MHPQKsmS6NSgFBAIxVBRxaIsZQKQRo1UroVAqolArIVQvF2haMdXjIsncHQiqVc0Pizg+C0k8ABSedHxJXvTACh8S1a14d6BjNtqFrzOgQEVH7pFFZb/glSQWDRUSlqfFBhy2A0KqtgY1KWf8HidUDLlvQY7R4pmqcLejSKK2BV/WgS68BZFk+tz6R6NlFWc8PbqoTBAEqpVDnDZh0fhAkyxBF61fJmzNkDQ2JK6+omhOUUXNeUFY2h8RRrbw60DFYmNEhIiICqldtsxYwqKzKtDR0w68QzmVtzg8gXHEu6IFbgh5BALRKpT2QUyrqb1f1Ut2AGkaLaA983FHEob7gxhUKhQAFBNT1GW27rbLnqwd69bA+zmexAFk5VUPgqs0Pyqj6WlZe/5C4oKCqLFA8h8R5Ga8OdGwZHc7RISIiOsdewEBb+4KktmFfWpV1P3c7P+gxOrE+UPXhaBo3tEurUlqrsuoAsyhVBT2uZbvcFdy4oqGArl1Sqc7N3TmfLAOnz9Q9JO7kKaCoCNhTBOzZV/P15w+JswVBHRM4JK4N8OpAhxkdIiKiulXPcoiSDLMoQaNUQNGMN9PK89YHqh701DcHyJ1sgZ9f1RA7o0Wsc16PSmHtM20zBjfUBIIAhHWwPi4YXPN5Donzal4d6HCODhERkXOUCgFKRcv+e1k96JFl2WOBjbNtkCQZJtE6vE2lFBjceCNnh8TZq8Nluj4krrYqcVGRHBLXDLw70GFGh4iIqE1qiSDnfAqFAJ1CyQ9M26uWGBKXmAAkxFufpybz6kDHYJujw0CHiIiIiNyloSFxZeXAiczGD4mLiap7SFxQoEcvzZt4daBjy+jwkxgiIiKiagR+COxRfr6uDYlLPwGcOHFuSFx2rvWxZVvN1wcHn1cYIZ5D4urg1YEOMzpEREREtVD7AAoVYCoDLEZ490qlrUxTh8SdPWt9/L235ut1uqohcbVkguLj2t2QOK8OdIwWFiMgIiIiqpVSDfgEA5IImCsAUwUgu38xV3KBM0PiMqoyPxnnDYnLzAIMBuBIqvVR27Fjo63D4drJkDgvD3RYjICIiIioXgoloPUHNH6AuRIwlQOSpaVbRbXx8wV697Q+zlfXkDhb+ezyCuvzWTnODYmrPjSujQ6J8+5AhwuGEhERETlHEACN3vqwGK0Bj8XY0q0iZ3FIXA1eG+jYat8DgI4ZHSIiIiLnqbTWh2gBzOXWTA/n8bRdrgyJc8gEuTIkruO5QKtj6xgS57WBjm1+DsCMDhEREVGjKFWAMhDQ+Fvn8ZgrrHN6yLu4bUjc1pqvr2tIXGICEBnh0SFxXhzonPslZEaHiIiIqAkUCkDrZ32YK62FC0RTS7eKmkNDQ+JOnbYGPBmZjRsS1zHesUS2G4fEeXGgY83oKBUCVEoGOkRERERuofaxPiwm67A2lqduvwQBCA+zPoYOqfl8bUPibEGRbUjc4aPWR23HPn9InD0TFO1U87w20DGYWXGNiIiIyGNUGutDkgBLpTXTI5pbulXUmnhySJwTvDbQ4Ro6RERERM1AoQA0vtaHaK6ay2PgmjxUP2eHxNVWJe7UaedO4c72tibM6BARERE1M6XaWrxAGwBYDNYsD0tUk6saGhJ38gzQpU+Dh/HaQIcZHSIiIqIWIgjn5vJIYlWWp5IV28g9/Pyc2s1rAx1mdIiIiIhaAYUS0PpbHxajNehhAQNqBi5FAfPnz4cgCA6P7t2717n/0qVLa+yv0+ma3GhnGM3WjA7X0CEiIiJqJVRawCcY8A0HdAHWoW5EHuJyRqdXr15Yt27duQOo6j9EQEAAjhw5Yv9eEARXT9kotqFrzOgQERERtTIsYEDNwOVAR6VSITIy0un9BUFwaX934dA1IiIiojbAVsBAF8ihbeRWLkcBqampiI6ORqdOnTBp0iRkZmbWu39ZWRkSEhIQFxeHa6+9FgcPHmzwHEajESUlJQ4PV7EYAREREVEbYxva5hcB+ARZvydqJJcCnaFDh2Lp0qX49ddf8d577yE9PR0XX3wxSktLa92/W7du+OSTT/DDDz/giy++gCRJGDZsGLKzs+s9z8KFCxEYGGh/xMXFudJMAMzoEBEREbVZtqpt+hBr0KMLAJSalm4VtTGCLDc+L1hUVISEhAS8+uqrmDZtWoP7m81m9OjRAxMnTsQzzzxT535GoxFG47ma6yUlJYiLi0NxcTECAgKcatt7G4/jxV8P44ZBsXjlxn5OvYaIiIiIWjHRAlgqrfN5JEtLt4ZaSEmFEYFRiQ3GBk0qLx0UFISuXbvi2LFjTu2vVqsxYMCABvfXarXQapuWqjRamNEhIiIi8ipKFaC0lao2nQt6WMSAatGkKKCsrAzHjx9HVFSUU/uLooj9+/c7vX9TGGzlpVWco0NERETkdVQaawED/wjrvB61zjrkjaiKS4HOww8/jE2bNiEjIwPbtm3DddddB6VSiYkTJwIAJk+ejHnz5tn3f/rpp/Hbb78hLS0Nf//9N2699VacOHECd955p3uvoha2jI5OzYwOERERkVdT6xyLGDDoIbg4dC07OxsTJ07EmTNnEBYWhosuugh//vknwsLCAACZmZlQKM4FFmfPnsVdd92F/Px8BAcHY9CgQdi2bRt69uzp3quoBTM6RERERO2MrYiB2sdantpiAMyVgGhiuep2yKVAZ/ny5fU+v3HjRofvX3vtNbz22msuN8odmNEhIiIiaseqBz2SZA16LAYGPe1Ik4oRtGa2dXRYjICIiIionVMoAI3e+qge9FiMDb+W2izvDXRs6+hwwVAiIiIisqkR9FRVbhNNLd0ycjPvDXSqMjocukZEREREtVIoAI2v9SGJVXN6GPR4C68NdAy2jA6LERARERFRQxRKBj1exmsDHWZ0iIiIiKhRHIKequFtFiPn9LQxXhvoMKNDRERERE3mMLyN1dvaEq8NdJjRISIiIiK3qqt6G4OeVsl7Ax0uGEpEREREnlI96OHipK2S1wY6Bott6BozOkRERETkQdUXJ7UFPbZ1ehj0tBivDXRsGR0d19EhIiIiouZSI+gxVgt6pJZuXbvilYGOLMvM6BARERFRyxIEQK2zPgDHoEcSW7Zt7YBXBjpmUbZnCbXM6BARERFRa6DSWh8AIJqtc3osRkCytGy7vJRXBjpGy7kImRkdIiIiImp1lGrrAwBEy7m1ekRzy7bLi3hloGMwnxv/yECHiIiIiFo1pQpQ+gNaf+uQNosBMFeVraZG88pAx1htfo4gCC3cGiIiIiIiJymUXKDUTbwy0DHY19BhNoeIiIiI2qgaa/UYq4a4mVjBzQleGejYMjosLU1EREREXoEV3FzmpYFOVUZHzYwOEREREXmh8yu42YIeFjOw88pAx2C2zdFhRoeIiIiIvJytglv1YgYWo/XRjnlloGPL6OiY0SEiIiKi9oTFDOy8M9BhRoeIiIiI2rtaixlUBT3tYF6PdwY6zOgQEREREZ1To5iBCRCN1vV6JEvLts1DvDLQ4RwdIiIiIqJ6qDTWhxfP6/HKQIcZHSIiIiIiJ50/r0esVrq6Dc/r8c5Ax75gKDM6REREREROUygAhQ+g9rEGOaKpza7X45WBzrmha8zoEBERERE1iiC06fV6vDLQOTd0jRkdIiIiIiK3qLFej7FVl652KeUxf/58CILg8OjevXu9r/n222/RvXt36HQ69OnTB6tXr25Sg53BjA4RERERkQcplNay1foQwC8C8Am2DndTtJ5Eg8uRQK9evZCXl2d//PHHH3Xuu23bNkycOBHTpk3Dnj17MH78eIwfPx4HDhxoUqMbYsvoaJnRISIiIiLyLFvpap8gwC8c8O0AaP2s2Z8W5PLQNZVKhcjISKf2feONNzBmzBjMmTMHAPDMM89g7dq1ePvtt/H++++7emqnGS3M6BARERERtYgaQ9xs83qad4iby5FAamoqoqOj0alTJ0yaNAmZmZl17rt9+3aMHj3aYVtKSgq2b99e7zmMRiNKSkocHq4w2KuuMdAhIiIiImoxttLV5w9xEzx/n+7SGYYOHYqlS5fi119/xXvvvYf09HRcfPHFKC0trXX//Px8REREOGyLiIhAfn5+vedZuHAhAgMD7Y+4uDhXmmnP6LAYARERERFRK1F9iJt/hMeHuLkU6IwdOxY33ngj+vbti5SUFKxevRpFRUX45ptv3NqoefPmobi42P7Iyspy6fXM6BARERERtXK24W2+Haxze3SB1lLWguCWwzepvHRQUBC6du2KY8eO1fp8ZGQkCgoKHLYVFBQ0OMdHq9VCq9W61BZJkmAymQAAvkoRMf5K+KlkGAwGl45DRO6n0WigUPCDByIiIqqDrYqbRm+dx2MxAqKxSQuVNinQKSsrw/Hjx3HbbbfV+nxycjLWr1+P2bNn27etXbsWycnJTTltDSaTCenp6ZAkaybnph4+uK6LFqGKEqSnl7v1XETkOoVCgcTERGg0mpZuChEREbV2tiFuap31e/tCpSZrQQMnuRToPPzww7j66quRkJCA3NxcPPXUU1AqlZg4cSIAYPLkyYiJicHChQsBALNmzcKIESOwaNEijBs3DsuXL8euXbuwePFiV05bL1mWkZeXB6VSibi4OOunxqfLYbKIiA32ga+2ZcvaEbV3kiQhNzcXeXl5iI+Ph+CmdDQRERG1E/YqbgAkCbAUOvUylwKd7OxsTJw4EWfOnEFYWBguuugi/PnnnwgLCwMAZGZmOgxPGTZsGL788ks8/vjj+L//+z906dIFq1atQu/evV05bb0sFgsqKioQHR0NvV4PAFCozBAgQqvzgU7bpKQVEblBWFgYcnNzYbFYoFbzwwciIiJqJIUCUDk3QkSQ5WYsZt1IJSUlCAwMRHFxMQICAhyeMxgMSE9PR8eOHeHj4wMAOJxXApMooXO4H/QaBjpELa2yshIZGRlITEyETqdr6eYQERFRG1ZfbFCd18wOrj4cRpJrbiOilsPfRSIiImpuXhPoVGdLUnnlxRERERERUYO8MhaQqr56y6fIHTt2xOuvv97SzXCbjRs3QhAEFBUVtXRTiIiIiMhLeV2gI8vyuYxOG4hzsrKyMHXqVERHR0Oj0SAhIQGzZs3CmTNnWrppbjFy5EiH8uKAtUhFXl4eAgMDW6ZRREREROT1vC7QkaqVVmjtGZ20tDQMHjwYqamp+Oqrr3Ds2DG8//77WL9+PZKTk1FY6FzpPHcTRdG+JpEnaDQaREZGtvqfDxERERG1XV4X6FQvIudqRic7Oxu///47srOz3dyq2s2YMQMajQa//fYbRowYgfj4eIwdOxbr1q1DTk4OHnvsMfu+paWlmDhxInx9fRETE4N33nnH/pwsy5g/fz7i4+Oh1WoRHR2N+++/3/680WjEww8/jJiYGPj6+mLo0KHYuHGj/fmlS5ciKCgIP/74I3r27AmtVouPPvoIOp2uxvCyWbNm4dJLLwUAnDlzBhMnTkRMTAz0ej369OmDr776yr7vlClTsGnTJrzxxhsQBAGCICAjI6PWoWsrVqxAr169oNVq0bFjRyxatMjhvB07dsTzzz+PqVOnwt/fH/Hx8Q7rMZlMJsycORNRUVHQ6XRISEiwr+dERERERO2P1wU69oprEFzKGHz88cdISEjApZdeioSEBHz88cceaqFVYWEh1qxZg+nTp9vLYttERkZi0qRJ+Prrr+2B28svv4x+/fphz549ePTRRzFr1iysXbsWgDVIeO211/DBBx8gNTUVq1atQp8+fezHmzlzJrZv347ly5fjn3/+wY033ogxY8YgNTXVvk9FRQVefPFFfPTRRzh48CAmTZqEoKAgrFixwr6PKIr4+uuvMWnSJADW0t6DBg3Czz//jAMHDuDuu+/Gbbfdhr/++gsA8MYbbyA5ORl33XUX8vLykJeXh7i4uBp9sXv3bkyYMAE333wz9u/fj/nz5+OJJ57A0qVLHfZbtGgRBg8ejD179mD69Om47777cOTIEQDAm2++iR9//BHffPMNjhw5gmXLlqFjx46N/OkQERERUZsntwHFxcUyALm4uLjGc5WVlfK///4rV1ZWyrIsywaTRd6XdVY+kF3k9PGzsrJkhUIhA7A/lEqlnJWV5bZrON+ff/4pA5BXrlxZ6/OvvvqqDEAuKCiQExIS5DFjxjg8f9NNN8ljx46VZVmWFy1aJHft2lU2mUw1jnPixAlZqVTKOTk5Dtsvu+wyed68ebIsy/KSJUtkAPLevXsd9pk1a5Z86aWX2r9fs2aNrNVq5bNnz9Z5XePGjZMfeugh+/cjRoyQZ82a5bDP77//LgOwH+eWW26RL7/8cod95syZI/fs2dP+fUJCgnzrrbfav5ckSQ4PD5ffe+89WZZl+b///a986aWXypIk1dk2ajnn/54SERERNVZ9sUF1XpfRsQ1ccyWbk5qaWmNOiiiKOHbsmBtbVjvZyfVak5OTa3x/6NAhAMCNN96IyspKdOrUCXfddRdWrlwJi8UCANi/fz9EUUTXrl3h5+dnf2zatAnHjx+3H0+j0aBv374O55g0aRI2btyI3NxcAMCyZcswbtw4BAUFAbD20TPPPIM+ffogJCQEfn5+WLNmDTIzM13qg0OHDmH48OEO24YPH47U1FSIomjfVr19giAgMjISJ0+eBGAdJrd3715069YN999/P3777TeX2kBERERE3sXrAh2pKnBwZZ57ly5doFA4doVSqUTnzp3d2TQHnTt3hiAI9mDlfIcOHUJwcDDCwsIaPFZcXByOHDmCd999Fz4+Ppg+fTouueQSmM1mlJWVQalUYvfu3di7d6/9cejQIbzxxhv2Y/j4+NQIDocMGYKkpCQsX74clZWVWLlypX3YGmAdTvfGG29g7ty5+P3337F3716kpKTAZDI1slfqp1arHb4XBMEeoA4cOBDp6el45plnUFlZiQkTJuCGG27wSDuIiIiIqPXzukDHliBRuBDpxMbGYvHixVAqlQCsQc4HH3yA2NhYTzQRABAaGorLL78c7777LiorKx2ey8/Px7Jly3DTTTfZg48///zTYZ8///wTPXr0sH/v4+ODq6++Gm+++SY2btyI7du3Y//+/RgwYABEUcTJkyfRuXNnh0dkZGSD7Zw0aRKWLVuG//3vf1AoFBg3bpz9ua1bt+Laa6/Frbfein79+qFTp044evSow+s1Go1DVqY2PXr0wNatWx22bd26FV27drX/TJwREBCAm266CR9++CG+/vprrFixosUq1xERERFRy1K1dAPcrTEZHQCYNm0aUlJScOzYMXTu3NmjQY7N22+/jWHDhiElJQXPPvssEhMTcfDgQcyZMwcxMTF47rnn7Ptu3boVL730EsaPH4+1a9fi22+/xc8//wzAWjVNFEUMHToUer0eX3zxBXx8fJCQkIDQ0FBMmjQJkydPxqJFizBgwACcOnUK69evR9++fR0Cl9pMmjQJ8+fPx3PPPYcbbrgBWq3W/lyXLl3w3XffYdu2bQgODsarr76KgoIC9OzZ075Px44dsWPHDmRkZMDPzw8hISE1zvHQQw9hyJAheOaZZ3DTTTdh+/btePvtt/Huu+863ZevvvoqoqKiMGDAACgUCnz77beIjIy0D7MjIiIiovaFGZ1qYmNjMXLkyGYJcgBroLBr1y506tQJEyZMQFJSEu6++26MGjUK27dvdwgKHnroIezatQsDBgzAs88+i1dffRUpKSkAgKCgIHz44YcYPnw4+vbti3Xr1uF///sfQkNDAQBLlizB5MmT8dBDD6Fbt24YP348du7cifj4+Abb2LlzZ1xwwQX4559/HIatAcDjjz+OgQMHIiUlBSNHjkRkZCTGjx/vsM/DDz8MpVKJnj17IiwsrNb5OwMHDsQ333yD5cuXo3fv3njyySfx9NNPY8qUKU73pb+/P1566SUMHjwYQ4YMQUZGBlavXl1jSCIRERERtQ+C7Oxs+BZUUlKCwMBAFBcXIyAgwOE5g8GA9PR0JCYmWtd9qTAhs7ACvloVksL8WqjFRFTd+b+nRERERI1VX2xQndd93C01IaNDRERERETewesCHVuCimEOEREREVH75XWBDjM6RERERETkdYGOjMZVXSMiIiIiIu/hdYHOuYxOy7aDiIiIiIhajtcFOvY5OkzpEBERERG1W14Y6Fi/MqNDRERERNR+eV2gIzGjQ0RERETU7nldoMOMDhEREREReV2gw4yOd9q4cSMEQUBRUVGzn3vp0qUICgpq9vPWRRAErFq1CgCQkZEBQRCwd+/eRh/PHccgIiIiam28LtBpSxmdrKwsTJ06FdHR0dBoNEhISMCsWbNw5syZlm4apkyZgvHjx7d0M9oMQRDsj8DAQAwfPhwbNmzw+Hnj4uKQl5eH3r17O7V/bT9XV49BRERE1BZ4XaDTVjI6aWlpGDx4MFJTU/HVV1/h2LFjeP/997F+/XokJyejsLCwpZtILlqyZAny8vKwdetWdOjQAVdddRXS0tJq3ddsNrvlnEqlEpGRkVCpVC16DCIiIqLWxusCHXtGp2Wb0aAZM2ZAo9Hgt99+w4gRIxAfH4+xY8di3bp1yMnJwWOPPWbft2PHjnj++ecxdepU+Pv7Iz4+HosXL3Y4XlZWFiZMmICgoCCEhITg2muvRUZGRr1t+O6779CnTx/4+PggNDQUo0ePRnl5OebPn49PP/0UP/zwgz1LsXHjRqfOY8sYLFiwAGFhYQgICMC9994Lk8lUb1s+//xzDB48GP7+/oiMjMQtt9yCkydP1thv9+7dGDx4MPR6PYYNG4YjR444PP/DDz9g4MCB0Ol06NSpExYsWACLxWJ//tVXX0WfPn3g6+uLuLg4TJ8+HWVlZQ7HWLp0KeLj46HX63Hdddc5nWELCgpCZGQkevfujffeew+VlZVYu3YtAGvg/d577+Gaa66Br68vnnvuOafam5qaiksuuQQ6nQ49e/a0H8+mtmFnBw8exFVXXYWAgAD4+/vj4osvxvHjx+v8udZ2jE2bNuGCCy6AVqtFVFQUHn30UYd2jRw5Evfffz8eeeQRhISEIDIyEvPnz3eqn4iIiIiaQ5PigRdeeAGCIGD27Nl17rN06VKHYT2CIECn0zXltPUSJQkGs4hKs4gKk6VZH7Y1fBpSWFiINWvWYPr06fDx8XF4LjIyEpMmTcLXX3/tcLxFixZh8ODB2LNnD6ZPn4777rvPfpNvNpuRkpICf39/bNmyBVu3boWfnx/GjBlTZ4CRl5eHiRMnYurUqTh06BA2btyI66+/HrIs4+GHH8aECRMwZswY5OXlIS8vD8OGDXP6POvXr7cf86uvvsL333+PBQsW1NsnZrMZzzzzDPbt24dVq1YhIyMDU6ZMqbHfY489hkWLFmHXrl1QqVSYOnWq/bktW7Zg8uTJmDVrFv7991988MEHWLp0qT2oAACFQoE333wTBw8exKeffooNGzbgkUcesT+/Y8cOTJs2DTNnzsTevXsxatQoPPvss/W2vTa2n2v1fpk/fz6uu+467N+/H1OnTm2wvZIk4frrr4dGo8GOHTvw/vvvY+7cufWeNycnB5dccgm0Wi02bNiA3bt3Y+rUqbBYLHX+XGs7xpVXXokhQ4Zg3759eO+99/Dxxx/X6IdPP/0Uvr6+2LFjB1566SU8/fTTNQIxIiIiopbS6LEqO3fuxAcffIC+ffs2uG9AQIDDJ++eHFZWaZYw4YM/PXb8+vz7dAr0moa7NDU1FbIso0ePHrU+36NHD5w9exanTp1CeHg4AODKK6/E9OnTAQBz587Fa6+9ht9//x3dunXD119/DUmS8NFHH9n7dsmSJQgKCsLGjRtxxRVX1DhHXl4eLBYLrr/+eiQkJAAA+vTpY3/ex8cHRqMRkZGR9m1ffPGFU+fRaDT45JNPoNfr0atXLzz99NOYM2cOnnnmGSgUtcfW1QOWTp064c0338SQIUNQVlYGPz8/+3PPPfccRowYAQB49NFHMW7cOBgMBuh0OixYsACPPvoobr/9dvtxnnnmGTzyyCN46qmnAMAhKO/YsSOeffZZ3HvvvXj33XcBAG+88QbGjBljD366du2Kbdu24ddff6213bWpqKjA448/DqVSaW8rANxyyy244447HK65vvauW7cOhw8fxpo1axAdHQ0AeP755zF27Ng6z/3OO+8gMDAQy5cvh1qttl+DTW0/1/O9++67iIuLw9tvvw1BENC9e3fk5uZi7ty5ePLJJ+0/w759+9r7tUuXLnj77bexfv16XH755U73FREREZGnNCqjU1ZWhkmTJuHDDz9EcHBwg/sLgoDIyEj7IyIiojGndYoM57IqrYGzGSAADgGlrT9tQ7v27duHY8eOwd/fH35+fvDz80NISAgMBgOOHz+OLVu22Lf7+flh2bJl6NevHy677DL06dMHN954Iz788EOcPXu23jY0dB6bfv36Qa/X279PTk5GWVkZsrKysGzZMoe2bNmyBYB1SNrVV1+N+Ph4+Pv72wOEzMzMOvshKioKABz64emnn3Y4/l133YW8vDxUVFQAANatW4fLLrsMMTEx8Pf3x2233YYzZ87Ynz906BCGDh3qcM7k5OSGfjwAgIkTJ8LPzw/+/v5YsWIFPv74Y4f2Dh48uEZ/1tfeQ4cOIS4uzh7kONOWvXv34uKLL7YHOY1x6NAhJCcnO3wgMXz4cJSVlSE7O9u+7fwPOaKiomodbkhERETUEhqV0ZkxYwbGjRuH0aNHOzWsp6ysDAkJCZAkCQMHDsTzzz+PXr16NebUDdIoFfjmnguRFOYLHyeyK+7ko1Y6tV/nzp0hCAIOHTqE6667rsbzhw4dQnBwMMLCwuzbzr9xFQQBkiQBsPbvoEGDsGzZshrHCgsLg0ajcZh/ERERAaVSibVr12Lbtm347bff8NZbb+Gxxx7Djh07kJiYWGu7GzqPM6655hqHQCImJgbl5eVISUlBSkoKli1bhrCwMGRmZiIlJaXG0Lvq/WC7Ea/eDwsWLMD1119f47w6nQ4ZGRm46qqrcN999+G5555DSEgI/vjjD0ybNg0mk8khOGuM1157DaNHj0ZgYGCt/eHr6+vwfUPtbYzzh0J6Un3vSSIiIqKW5nIksHz5cvz999/YuXOnU/t369YNn3zyCfr27Yvi4mK88sorGDZsGA4ePIjY2NhaX2M0GmE0Gu3fl5SUON0+GYBOrYSvVg2dk4FHcwsNDcXll1+Od999Fw888IDDzWl+fj6WLVuGyZMnOz3Eb+DAgfj6668RHh6OgICAWvfp3LlzjW2CIGD48OEYPnw4nnzySSQkJGDlypV48MEHodFoIIqiy+cBrJmKyspK+3X9+eef8PPzQ1xcHBQKBfz9/R323717N86cOYMXXngBcXFxAIBdu3Y5de3nt+/IkSO1XqvtPJIkYdGiRfbhV998843DPj169MCOHTsctv35p3NDISMjI+s8d2Pa26NHD2RlZSEvL8+evWqoLX379sWnn34Ks9lca1antp9rbeddsWIFZFm2vwe3bt0Kf3//On9niYiIiFobl4auZWVlYdasWVi2bJnTnzgnJydj8uTJ6N+/P0aMGIHvv/8eYWFh+OCDD+p8zcKFCxEYGGh/2G5+ndFW1tF5++23YTQakZKSgs2bNyMrKwu//vorLr/8csTExDhMoG/IpEmT0KFDB1x77bXYsmUL0tPTsXHjRtx///0OQ42q27FjB55//nns2rULmZmZ+P7773Hq1Cn7vKGOHTvin3/+wZEjR3D69GmYzWanz2MymTBt2jT8+++/WL16NZ566inMnDmzzvk58fHx0Gg0eOutt5CWloYff/wRzzzzjAu9afXkk0/is88+w4IFC3Dw4EEcOnQIy5cvx+OPPw7AGuyZzWb7eT7//HO8//77Dse4//778euvv+KVV15Bamoq3n77bZfm57izvaNHj0bXrl1x++23Y9++fdiyZYtDNb7azJw5EyUlJbj55puxa9cupKam4vPPP7fPkavt53q+6dOnIysrC//9739x+PBh/PDDD3jqqafw4IMP1vkzJCIiImptXLpr2b17N06ePImBAwdCpVJBpVJh06ZNePPNN6FSqRr8pBiwDncZMGAAjh07Vuc+8+bNQ3Fxsf2RlZXlVPtkWW4z6+h06dIFu3btQqdOnTBhwgQkJSXh7rvvxqhRo7B9+3aEhIQ4fSy9Xo/NmzcjPj4e119/PXr06IFp06bBYDDUmXkJCAjA5s2bceWVV6Jr1654/PHHsWjRIvtE97vuugvdunXD4MGDERYWhq1btzp9nssuuwxdunTBJZdcgptuugnXXHNNvaWHw8LCsHTpUnz77bfo2bMnXnjhBbzyyitOX79NSkoKfvrpJ/z2228YMmQILrzwQrz22mv2Ygv9+vXDq6++ihdffBG9e/fGsmXLsHDhQodjXHjhhfjwww/xxhtvoF+/fvjtt9/sgYe7NdRehUKBlStXorKyEhdccAHuvPPOBgPg0NBQbNiwAWVlZRgxYgQGDRqEDz/80J7dqe3ner6YmBisXr0af/31F/r164d7770X06ZN81g/EBEREXmCILswI760tBQnTpxw2HbHHXege/fumDt3rlMrq4uiiF69euHKK6/Eq6++6tR5S0pKEBgYiOLi4ho37gaDAenp6UhMTIRGq8WBnGIAQK/oACj56XOzmzJlCoqKirBq1aqWbgq1ItV/Tz1ZXp6IiIi8X32xQXUuzdHx9/evEcz4+voiNDTUvn3y5MmIiYmxf1L+9NNP48ILL0Tnzp1RVFSEl19+GSdOnMCdd97p6jU1SKoWs7X2jA4REREREXmO28uSZWZmOozjP3v2LO666y7k5+cjODgYgwYNwrZt29CzZ093nxrVc1MMc4iIiIiI2i+Xhq61FGeHrilUahzOL4VCENA7JrCFWktE5+PQNSIiInIXZ4euedUkFqkqZOOoNSIiIiKi9s2rAh1bckrBSIeIiIiIqF3zqkCHGR0iIiIiIgK8LNBhRoeIiIiIiAAvC3TsGZ2WbQYREREREbUwrwp0mNEhIiIiIiLAywIdqeor4xwiIiIiovbNqwKdtpTRmTJlCgRBqPE4duxYSzetUZYuXYqgoKCWbgYREREREQBA1dINcKe2VnVtzJgxWLJkicO2sLAwl49jMpmg0Wjc1SwiIiIiojaPGZ0WpNVqERkZ6fBQKpXYtGkTLrjgAmi1WkRFReHRRx+FxWKxv27kyJGYOXMmZs+ejQ4dOiAlJQUAcODAAYwdOxZ+fn6IiIjAbbfdhtOnT9tfJ0kSXnrpJXTu3BlarRbx8fF47rnn7M/PnTsXXbt2hV6vR6dOnfDEE0/AbDbbn9+3bx9GjRoFf39/BAQEYNCgQdi1axc2btyIO+64A8XFxfbM1Pz58z3fgUREREREdfC+jI4sQ1FRDmjl5m+AXt/kdFJOTg6uvPJKTJkyBZ999hkOHz6Mu+66CzqdziF4+PTTT3Hfffdh69atAICioiJceumluPPOO/Haa6+hsrISc+fOxYQJE7BhwwYAwLx58/Dhhx/itddew0UXXYS8vDwcPnzYfkx/f38sXboU0dHR2L9/P+666y74+/vjkUceAQBMmjQJAwYMwHvvvQelUom9e/dCrVZj2LBheP311/Hkk0/iyJEjAAA/P78m9QMRERERUVMIsi0N0oqVlJQgMDAQxcXFCAgIcHjOYDAgPT0diYmJKDLKOFVQiD7dYlumoWVlgK+vU7tOmTIFX3zxBXQ6nX3b2LFj0bVrV6xYsQKHDh2CUBU0vfvuu5g7dy6Ki4uhUCgwcuRIlJSU4O+//7a/9tlnn8WWLVuwZs0a+7bs7GzExcXhyJEjiIqKQlhYGN5++23ceeedTrXxlVdewfLly7Fr1y4AQEBAAN566y3cfvvtNfZdunQpZs+ejaKiIqeOTe1L9d/T6u95IiIiIlfVFxtU530ZnTZk1KhReO+99+zf+/r6YsaMGUhOTrYHOQAwfPhwlJWVITs7G/Hx8QCAQYMGORxr3759+P3332vNpBw/fhxFRUUwGo247LLL6mzP119/jTfffBPHjx9HWVkZLBaLw5vnwQcfxJ133onPP/8co0ePxo033oikpKRGXz8RERERkad4VaAjyzJkHz0Kck8jIqAFPjXW613a3dfXF507d27UqXzPyxyVlZXh6quvxosvvlhj36ioKKSlpdV7vO3bt2PSpElYsGABUlJSEBgYiOXLl2PRokX2febPn49bbrkFP//8M3755Rc89dRTWL58Oa677rpGXQMRERERkad4VaAjyQAEAYKfL+DbNofH9OjRAytWrIAsy/asztatW+Hv74/Y2LqH5A0cOBArVqxAx44doVLV/LF26dIFPj4+WL9+fa1D17Zt24aEhAQ89thj9m0nTpyosV/Xrl3RtWtXPPDAA5g4cSKWLFmC6667DhqNBqIoNuaSiYiIiIjczsuqrlm/KtA2qq7VZvr06cjKysJ///tfHD58GD/88AOeeuopPPjgg1Ao6v5xzZgxA4WFhZg4cSJ27tyJ48ePY82aNbjjjjsgiiJ0Oh3mzp2LRx55BJ999hmOHz+OP//8Ex9//DEAayCUmZmJ5cuX4/jx43jzzTexcuVK+/ErKysxc+ZMbNy4ESdOnMDWrVuxc+dO9OjRAwDQsWNHlJWVYf369Th9+jQqKio821FERERERPXwqkBHqop02kh16VrFxMRg9erV+Ouvv9CvXz/ce++9mDZtGh5//PF6XxcdHY2tW7dCFEVcccUV6NOnD2bPno2goCB7gPTEE0/goYcewpNPPokePXrgpptuwsmTJwEA11xzDR544AHMnDkT/fv3x7Zt2/DEE0/Yj69UKnHmzBlMnjwZXbt2xYQJEzB27FgsWLAAADBs2DDce++9uOmmmxAWFoaXXnrJQz1ERERERNQwr6q6llcmotRgRmywHiG+XECTqLVg1TUiIiJyF2errnllRkfRhjM6RERERETUdF4V6Njn6LTlsWtERERERNRkXhXoeMMcHSIiIiIiajqvCnSY0SEiIiIiIsDLAh1mdIiIiIiICPCiQEeWZWZ0iFqpNlDckYiIiLyMqqUb0FRqtRqCIODUqVOwSDrIAEwGAyAqW7ppRARrkHPq1CkIggC1Wt3SzSEiIqJ2os0HOkqlErGxscjOzkZB4WnrtnIdlKwxTdRqCIKA2NhYKJX8AIKIiIiaR5sPdADAz88PiZ2ScMd366EQgO/vG4ZAPRcMJWot1Go1gxwiIiJqVl4R6ACARRaQVyYCAPz99NBpvObSiIiIiIjIRU0qRvDCCy9AEATMnj273v2+/fZbdO/eHTqdDn369MHq1aubctpaGS2S/f+1Kn5yTERERETUnjU60Nm5cyc++OAD9O3bt979tm3bhokTJ2LatGnYs2cPxo8fj/Hjx+PAgQONPXWtDGZrNketFDg/h4iIiIionWtUoFNWVoZJkybhww8/RHBwcL37vvHGGxgzZgzmzJmDHj164JlnnsHAgQPx9ttvN6rBdbFldHTM5hARERERtXuNmsgyY8YMjBs3DqNHj8azzz5b777bt2/Hgw8+6LAtJSUFq1atqvM1RqMRRqPR/n1xcTEAoKSkpM7XnDlbAslYAZVKXe9+RERERETUdtnu9Rtap8/lQGf58uX4+++/sXPnTqf2z8/PR0REhMO2iIgI5Ofn1/mahQsXYsGCBTW2x8XFNXi+LACBzznVNCIiIiIiaqNKS0sRGBhY5/MuBTpZWVmYNWsW1q5dC51O1+TG1WXevHkOWaCioiIkJCQgMzOz3otpyJAhQ5wO0NrKcdxxjJKSEsTFxSErKwsBAQEt2pbWdhz2r2eP05r6113tYf96tj3sX8+2pzX1rzuOw/717HHYv549Dvu3brIsY9CgQYiOjq53P5cCnd27d+PkyZMYOHCgfZsoiti8eTPefvttGI3GGmtlREZGoqCgwGFbQUEBIiMj6zyPVquFVqutsT0wMLBJP2ilUtnkN0prO4672gIAAQEB7F8PtQVg/3qyLUDT+9dd7WH/erY97F/Ptqc19a87j8P+9exx2L+ePQ77t3YajQYKRf3lBlwqRnDZZZdh//792Lt3r/0xePBgTJo0CXv37q11QcDk5GSsX7/eYdvatWuRnJzsyqndYsaMGV53HHe1xR1aU7+46zjsX88epzX1L9C6rqk1tcVdWtM1taa2uEtruqbWdhx3YP96FvvXs9pr/wpyQ7N4GjBy5Ej0798fr7/+OgBg8uTJiImJwcKFCwFYy0uPGDECL7zwAsaNG4fly5fj+eefx99//43evXs7dY6SkhIEBgaiuLjYbZ+e0TnsX89i/3oW+9ez2L+exf71LPavZ7F/PYv923RNWjC0NpmZmcjLy7N/P2zYMHz55ZdYvHgx+vXrh++++w6rVq1yOsgBrEPZnnrqqVqHs1HTsX89i/3rWexfz2L/ehb717PYv57F/vUs9m/TNTmjQ0RERERE1Nq4PaNDRERERETU0hjoEBERERGR12GgQ0REREREXoeBDhEREREReZ1mC3Q2b96Mq6++GtHR0RAEAatWrXJ4vqCgAFOmTEF0dDT0ej3GjBmD1NRU+/MZGRkQBKHWx7fffmvfLzMzE+PGjYNer0d4eDjmzJkDi8XSXJfZYpqrf2t7fvny5c11mS2mqf0LAPn5+bjtttsQGRkJX19fDBw4ECtWrHDYp7CwEJMmTUJAQACCgoIwbdo0lJWVefryWlxz9W/Hjh1rvH9feOEFT19ei3NH/x4/fhzXXXcdwsLCEBAQgAkTJtRYDJrvX8/2b3t9/y5cuBBDhgyBv78/wsPDMX78eBw5csRhH4PBgBkzZiA0NBR+fn74z3/+U6P/nLk/2LhxIwYOHAitVovOnTtj6dKlnr68Ftdc/btx48Za7yHy8/Ob5Tpbirv69/7778egQYOg1WrRv3//Ws/1zz//4OKLL4ZOp0NcXBxeeuklT11Wm9FsgU55eTn69euHd955p8Zzsixj/PjxSEtLww8//IA9e/YgISEBo0ePRnl5OQAgLi4OeXl5Do8FCxbAz88PY8eOBQCIoohx48bBZDJh27Zt+PTTT7F06VI8+eSTzXWZLaY5+tdmyZIlDvuNHz++OS6xRTW1fwHrGlNHjhzBjz/+iP379+P666/HhAkTsGfPHvs+kyZNwsGDB7F27Vr89NNP2Lx5M+6+++5mucaW1Fz9CwBPP/20w/v3v//9r8evr6U1tX/Ly8txxRVXQBAEbNiwAVu3boXJZMLVV18NSZLsx+L717P9C7TP9++mTZswY8YM/Pnnn1i7di3MZjOuuOIKh9//Bx54AP/73//w7bffYtOmTcjNzcX1119vf96Z+4P09HSMGzcOo0aNwt69ezF79mzceeedWLNmTbNeb3Nrrv61OXLkiMN7ODw8vFmus6W4o39tpk6diptuuqnW85SUlOCKK65AQkICdu/ejZdffhnz58/H4sWLPXZtbYLcAgDIK1eutH9/5MgRGYB84MAB+zZRFOWwsDD5ww8/rPM4/fv3l6dOnWr/fvXq1bJCoZDz8/Pt29577z05ICBANhqN7r2IVsxT/Vvbsdujxvavr6+v/NlnnzkcKyQkxL7Pv//+KwOQd+7caX/+l19+kQVBkHNycjx0Na2Pp/pXlmU5ISFBfu211zzW9ragMf27Zs0aWaFQyMXFxfZ9ioqKZEEQ5LVr18qyzPevjaf6V5b5/rU5efKkDEDetGmTLMvWvlKr1fK3335r3+fQoUMyAHn79u2yLDt3f/DII4/IvXr1cjjXTTfdJKekpHj6kloVT/Xv77//LgOQz54923wX0wo1pn+re+qpp+R+/frV2P7uu+/KwcHBDve7c+fOlbt16+b+i2hDWsUcHaPRCADQ6XT2bQqFAlqtFn/88Uetr9m9ezf27t2LadOm2bdt374dffr0QUREhH1bSkoKSkpKcPDgQQ+1vvVzV//azJgxAx06dMAFF1yATz75BHI7X4rJ2f4dNmwYvv76axQWFkKSJCxfvhwGgwEjR44EYH3/BgUFYfDgwfbXjB49GgqFAjt27Giei2mF3NW/Ni+88AJCQ0MxYMAAvPzyy+1iaGt9nOlfo9EIQRAcFq3T6XRQKBT2ffj+rZ27+teG71+guLgYABASEgLA+u+V2WzG6NGj7ft0794d8fHx2L59OwDn7g+2b9/ucAzbPrZjtBee6l+b/v37IyoqCpdffjm2bt3q6ctpdRrTv87Yvn07LrnkEmg0Gvu2lJQUHDlyBGfPnnVT69ueVhHo2H6g8+bNw9mzZ2EymfDiiy8iOzsbeXl5tb7m448/Ro8ePTBs2DD7tvz8fIdfMgD27719DGh93NW/gHXYxDfffIO1a9fiP//5D6ZPn4633nqrOS6j1XK2f7/55huYzWaEhoZCq9XinnvuwcqVK9G5c2cA1vfo+Sl8lUqFkJAQvn/d0L+AdYzz8uXL8fvvv+Oee+7B888/j0ceeaQlLqvVcKZ/L7zwQvj6+mLu3LmoqKhAeXk5Hn74YYiiaN+H79/auat/Ab5/AUCSJMyePRvDhw9H7969AVjfexqNBkFBQQ77RkRE2N97ztwf1LVPSUkJKisrPXE5rY4n+zcqKgrvv/8+VqxYgRUrViAuLg4jR47E33//7eGraj0a27/O4D1w7VpFoKNWq/H999/j6NGjCAkJgV6vx++//46xY8dCoajZxMrKSnz55Ze1ZhuoJnf27xNPPIHhw4djwIABmDt3Lh555BG8/PLLzXEZrZaz/fvEE0+gqKgI69atw65du/Dggw9iwoQJ2L9/fwu2vvVzZ/8++OCDGDlyJPr27Yt7770XixYtwltvvWX/1L09cqZ/w8LC8O233+J///sf/Pz8EBgYiKKiIgwcOLDWvyF0jjv7l+9f64iCAwcOtIsiOC3Bk/3brVs33HPPPRg0aBCGDRuGTz75BMOGDcNrr73m9nO1Vnz/Nj9VSzfAZtCgQdi7dy+Ki4thMpkQFhaGoUOHOgyDsPnuu+9QUVGByZMnO2yPjIzEX3/95bDNVrUiMjLSc41vA9zRv7UZOnQonnnmGRiNRodhF+1NQ/17/PhxvP322zhw4AB69eoFAOjXrx+2bNmCd955B++//z4iIyNx8uRJh+NaLBYUFhby/euG/q3N0KFDYbFYkJGRgW7dujXb9bQ2zvx9uOKKK3D8+HGcPn0aKpUKQUFBiIyMRKdOnQCA7996uKN/a9Pe3r8zZ860F7mIjY21b4+MjITJZEJRUZHDp+IFBQX2954z9weRkZE1Kl0VFBQgICAAPj4+nrikVsXT/VubCy64oM4h9N6mKf3rjLrev7bn2qtW91FcYGAgwsLCkJqail27duHaa6+tsc/HH3+Ma665BmFhYQ7bk5OTsX//fod/bNeuXYuAgAD07NnT421vC5rSv7XZu3cvgoOD23WQU11d/VtRUQEANT79ViqV9qpKycnJKCoqwu7du+3Pb9iwAZIkYejQoc10Ba1bU/q3Nnv37oVCofD6qj/OcubvQ4cOHRAUFIQNGzbg5MmTuOaaawDw/euMpvRvbdrL+1eWZcycORMrV67Ehg0bkJiY6PD8oEGDoFarsX79evu2I0eOIDMzE8nJyQCcuz9ITk52OIZtH9sxvFVz9W9t9u7di6ioKDdfUevijv51RnJyMjZv3gyz2WzftnbtWnTr1g3BwcFNv5C2qrmqHpSWlsp79uyR9+zZIwOQX331VXnPnj3yiRMnZFmW5W+++Ub+/fff5ePHj8urVq2SExIS5Ouvv77GcVJTU2VBEORffvmlxnMWi0Xu3bu3fMUVV8h79+6Vf/31VzksLEyeN2+ex6+vpTVH//7444/yhx9+KO/fv19OTU2V3333XVmv18tPPvmkx6+vpTW1f00mk9y5c2f54osvlnfs2CEfO3ZMfuWVV2RBEOSff/7Zvt+YMWPkAQMGyDt27JD/+OMPuUuXLvLEiROb/XqbW3P077Zt2+TXXntN3rt3r3z8+HH5iy++kMPCwuTJkye3yDU3J3f8ffjkk0/k7du3y8eOHZM///xzOSQkRH7wwQcd9uH713P9257fv/fdd58cGBgob9y4Uc7Ly7M/Kioq7Pvce++9cnx8vLxhwwZ5165dcnJyspycnGx/3pn7g7S0NFmv18tz5syRDx06JL/zzjuyUqmUf/3112a93ubWXP372muvyatWrZJTU1Pl/fv3y7NmzZIVCoW8bt26Zr3e5uaO/pVl6/3Znj175HvuuUfu2rWr/W+OrcpaUVGRHBERId92223ygQMH5OXLl8t6vV7+4IMPmvV6W5tmC3RsZQXPf9x+++2yLMvyG2+8IcfGxspqtVqOj4+XH3/88VpLQs+bN0+Oi4uTRVGs9TwZGRny2LFjZR8fH7lDhw7yQw89JJvNZk9eWqvQHP37yy+/yP3795f9/PxkX19fuV+/fvL7779f58/Cm7ijf48ePSpff/31cnh4uKzX6+W+ffvWKId85swZeeLEibKfn58cEBAg33HHHXJpaWlzXWaLaY7+3b17tzx06FA5MDBQ1ul0co8ePeTnn39eNhgMzXmpLcId/Tt37lw5IiJCVqvVcpcuXeRFixbJkiQ57MP3r+f6tz2/f2vrWwDykiVL7PtUVlbK06dPl4ODg2W9Xi9fd911cl5ensNxnLk/+P333+X+/fvLGo1G7tSpk8M5vFVz9e+LL74oJyUlyTqdTg4JCZFHjhwpb9iwobkus8W4q39HjBhR63HS09Pt++zbt0++6KKLZK1WK8fExMgvvPBCM11l6yXIcjuvDUxERERERF6n1c3RISIiIiIiaioGOkRERERE5HUY6BARERERkddhoENERERERF6HgQ4REREREXkdBjpEREREROR1GOgQEREREZHXYaBDREREREReh4EOERERERF5HQY6RERERETkdRjoEBERERGR12GgQ0REREREXuf/AV8E24pGnUVdAAAAAElFTkSuQmCC", "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.13" } }, "nbformat": 4, "nbformat_minor": 4 }