{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Persistent Shock Model and Income Expectations\n",
    "\n",
    "[![badge](https://img.shields.io/badge/Launch%20using%20-Econ--ARK-blue)](https://econ-ark.org/materials/incexpectationexample#launch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This module creates an example application extending $\\texttt{PersistentShockConsumerType}$ from $\\texttt{ConsGenIndShockModel}$. It uses the HARK tool $\\texttt{GenIncProcessModel}$ (whose documentation you can find [here](https://github.com/econ-ark/DemARK/blob/master/notebooks/IncExpectationExample.ipynb).)\n",
    "\n",
    "Most simply, it solves a consumption-saving model with shocks that are neither necessarily fully transitory nor fully permanent. Persistent income is tracked as a state variable and follows an AR(1) process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  What if beliefs about persistent income differ from actual ?\n",
    "\n",
    "What if the consumer has beliefs about the persistence of his/her income process which differ from the *actual* persistence?\n",
    "\n",
    "We can use the class $\\texttt{PersistentShockConsumerType}$ to solve the problem of a consumer with a given set of beliefs, but then simulate a population of consumers for whom that actual persistence differs from what they believe.\n",
    "\n",
    "(This thought experiment is motivated by an interesting recennt paper presented at the NBER Summer Institute's _Behavioral Macroeconomics Conference_ <cite data-cite=\"undefined\"></cite>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial imports and notebook setup, click arrow to show\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from copy import copy\n",
    "import numpy as np\n",
    "\n",
    "from HARK.ConsumptionSaving.ConsGenIncProcessModel import PersistentShockConsumerType\n",
    "from HARK.distribution import Uniform\n",
    "from HARK.utilities import get_lorenz_shares, calc_subpop_avg\n",
    "\n",
    "\n",
    "def mystr(number):\n",
    "    return \"{:.4f}\".format(number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell makes a subclass of PersistentShockConsumerType including the MPC\n",
    "\n",
    "\n",
    "class PersistentShockConsumerTypeX(PersistentShockConsumerType):\n",
    "    def getControls(self):\n",
    "        cLvl = np.zeros(self.AgentCount) + np.nan\n",
    "        MPCnow = np.zeros(self.AgentCount) + np.nan\n",
    "        for t in range(self.T_cycle):\n",
    "            these = t == self.t_cycle\n",
    "            cLvl[these] = self.solution[t].cFunc(\n",
    "                self.state_now[\"mLvlNow\"][these], self.state_now[\"pLvl\"][these]\n",
    "            )\n",
    "            MPCnow[these] = self.solution[t].cFunc.derivativeX(\n",
    "                self.state_now[\"mLvlNow\"][these], self.state_now[\"pLvl\"][these]\n",
    "            )\n",
    "        self.controls[\"cLvl\"] = cLvl\n",
    "        self.MPCnow = MPCnow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell defines a dictionary to make an instance of \"persistent shocks\" consumer for the infinite horizon model\n",
    "BaselineDict = {\n",
    "    \"CRRA\": 2.0,  # Coefficient of relative risk aversion\n",
    "    \"Rfree\": 1.01 / (1.0 - 1.0 / 160.0),  # Interest factor on assets\n",
    "    \"DiscFac\": 0.97,  # Default intertemporal discount factor\n",
    "    \"LivPrb\": [1.0 - 1.0 / 160.0],  # Survival probability\n",
    "    # Number of agents of this type (only matters for simulation)\n",
    "    \"AgentCount\": 10000,\n",
    "    \"IndL\": 10.0 / 9.0,  # Labor supply per individual (constant)\n",
    "    \"aNrmInitMean\": np.log(\n",
    "        0.00001\n",
    "    ),  # Mean of log initial assets (only matters for simulation)\n",
    "    # Standard deviation of log initial assets (only for simulation)\n",
    "    \"aNrmInitStd\": 0.0,\n",
    "    # Mean of log initial permanent income (only matters for simulation)\n",
    "    \"pLvlInitMean\": 0.0,\n",
    "    # Standard deviation of log initial permanent income (only matters for simulation)\n",
    "    \"pLvlInitStd\": 0.0,\n",
    "    # Aggregate permanent income growth factor (only matters for simulation)\n",
    "    \"PermGroFacAgg\": 1.0,\n",
    "    \"T_age\": 400,  # Age after which simulated agents are automatically killed\n",
    "    \"T_cycle\": 1,  # Number of periods in the cycle for this agent type\n",
    "    # Number of periods to simulate (idiosyncratic shocks model, perpetual youth)\n",
    "    \"T_sim\": 1200,\n",
    "    \"aXtraMin\": 0.001,  # Minimum end-of-period \"assets above minimum\" value\n",
    "    \"aXtraMax\": 30,  # Maximum end-of-period \"assets above minimum\" value\n",
    "    \"aXtraExtra\": [\n",
    "        0.005,\n",
    "        0.01,\n",
    "    ],  # Some other value of \"assets above minimum\" to add to the grid\n",
    "    # Exponential nesting factor when constructing \"assets above minimum\" grid\n",
    "    \"aXtraNestFac\": 3,\n",
    "    \"aXtraCount\": 48,  # Number of points in the grid of \"assets above minimum\"\n",
    "    \"PermShkCount\": 7,  # Number of points in discrete approximation to permanent income shocks\n",
    "    \"TranShkCount\": 7,  # Number of points in discrete approximation to transitory income shocks\n",
    "    \"PermShkStd\": [\n",
    "        (0.01 * 4 / 11) ** 0.5\n",
    "    ],  # Standard deviation of permanent shocks to income\n",
    "    \"TranShkStd\": [\n",
    "        (0.01 * 4) ** 0.5\n",
    "    ],  # Standard deviation of transitory shocks to income\n",
    "    \"UnempPrb\": 0.05,  # Probability of unemployment while working\n",
    "    \"UnempPrbRet\": 0.005,  # Probability of \"unemployment\" while retired\n",
    "    \"IncUnemp\": 0.3,  # Unemployment benefits replacement rate\n",
    "    \"IncUnempRet\": 0.0,  # \"Unemployment\" benefits when retired\n",
    "    \"tax_rate\": 0.0,  # Flat income tax rate\n",
    "    \"T_retire\": 0,  # Period of retirement (0 --> no retirement)\n",
    "    # Artificial borrowing constraint; imposed minimum level of end-of period assets\n",
    "    \"BoroCnstArt\": 0.0,\n",
    "    # Use cubic spline interpolation when True, linear interpolation when False\n",
    "    \"CubicBool\": False,\n",
    "    \"vFuncBool\": True,  # Whether to calculate the value function during solution\n",
    "    \"cycles\": 0,  # Make this type have an infinite horizon\n",
    "    \"pLvlPctiles\": np.concatenate(\n",
    "        (\n",
    "            [0.001, 0.005, 0.01, 0.03],\n",
    "            np.linspace(0.05, 0.95, num=19),\n",
    "            [0.97, 0.99, 0.995, 0.999],\n",
    "        )\n",
    "    ),\n",
    "    # Permanent income growth factor (no perm growth)\n",
    "    \"PermGroFac\": [1.000**0.25],\n",
    "    \"PrstIncCorr\": 0.99,  # Serial correlation coefficient for persistence of income\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell defines a function to solve and simulate a consumer misperceiving the correlation of persistent income shocks\n",
    "\n",
    "\n",
    "def runRoszypalSchlaffmanExperiment(CorrAct, CorrPcvd, DiscFac_center, DiscFac_spread):\n",
    "    \"\"\"\n",
    "    Solve and simulate a consumer type who misperceives the extent of serial correlation\n",
    "    of persistent shocks to income.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    CorrAct : float\n",
    "        Serial correlation coefficient for *actual* persistent income.\n",
    "    CorrPcvd : float\n",
    "        List or array of *perceived* persistent income.\n",
    "    DiscFac_center : float\n",
    "        A measure of centrality for the distribution of the beta parameter, DiscFac.\n",
    "    DiscFac_spread : float\n",
    "        A measure of spread or diffusion for the distribution of the beta parameter.\n",
    "\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    AggWealthRatio: float\n",
    "        Ratio of Aggregate wealth to income.\n",
    "    Lorenz: numpy.array\n",
    "        A list of two 1D array reprensenting the Lorenz curve for assets in the most recent simulated period.\n",
    "    Gini: float\n",
    "        Gini coefficient for assets in the most recent simulated period.\n",
    "    Avg_MPC: numpy.array\n",
    "        Average marginal propensity to consume by income quintile in the latest simulated period.\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    # Make a dictionary to construct our consumer type\n",
    "    ThisDict = copy(BaselineDict)\n",
    "    ThisDict[\"PrstIncCorr\"] = CorrAct\n",
    "\n",
    "    # Make a 7 point approximation to a uniform distribution of DiscFac\n",
    "    DiscFac_list = (\n",
    "        Uniform(\n",
    "            bot=DiscFac_center - DiscFac_spread, top=DiscFac_center + DiscFac_spread\n",
    "        )\n",
    "        .discretize(N=7)\n",
    "        .atoms\n",
    "    )\n",
    "\n",
    "    type_list = []\n",
    "    # Make a PersistentShockConsumerTypeX for each value of beta saved in DiscFac_list\n",
    "    for i in range(len(DiscFac_list)):\n",
    "        ThisDict[\"DiscFac\"] = DiscFac_list.flatten()[i]\n",
    "        ThisType = PersistentShockConsumerTypeX(**ThisDict)\n",
    "\n",
    "        # Make the consumer type *believe* he will face a different level of persistence\n",
    "        ThisType.PrstIncCorr = CorrPcvd\n",
    "        # Now he *thinks* E[p_{t+1}] as a function of p_t is different than it is\n",
    "        ThisType.update_pLvlNextFunc()\n",
    "\n",
    "        # Solve the consumer's problem with *perceived* persistence\n",
    "        ThisType.solve()\n",
    "\n",
    "        # Make the consumer type experience the true level of persistence during simulation\n",
    "        ThisType.PrstIncCorr = CorrAct\n",
    "        ThisType.update_pLvlNextFunc()\n",
    "\n",
    "        # Simulate the agents for many periods\n",
    "        ThisType.T_sim = 100\n",
    "        # ThisType.track_vars = ['cLvl','aLvl','pLvl','MPCnow']\n",
    "        ThisType.initialize_sim()\n",
    "        ThisType.simulate()\n",
    "        type_list.append(ThisType)\n",
    "\n",
    "    # Get the most recent simulated values of X = cLvl, MPCnow, aLvl, pLvl for all types\n",
    "    np.concatenate([ThisType.controls[\"cLvl\"] for ThisType in type_list])\n",
    "    aLvl_all = np.concatenate([ThisType.state_now[\"aLvl\"] for ThisType in type_list])\n",
    "    MPC_all = np.concatenate([ThisType.MPCnow for ThisType in type_list])\n",
    "    pLvl_all = np.concatenate([ThisType.state_now[\"pLvl\"] for ThisType in type_list])\n",
    "\n",
    "    # The ratio of aggregate assets over the income\n",
    "    AggWealthRatio = np.mean(aLvl_all) / np.mean(pLvl_all)\n",
    "\n",
    "    # first 1D array: Create points in the range (0,1)\n",
    "    wealth_percentile = np.linspace(0.001, 0.999, 201)\n",
    "\n",
    "    # second 1D array: Compute Lorenz shares for the created points\n",
    "    Lorenz_init = get_lorenz_shares(aLvl_all, percentiles=wealth_percentile)\n",
    "\n",
    "    # Stick 0 and 1 at the boundaries of both arrays to make it inclusive on the range [0,1]\n",
    "    Lorenz_init = np.concatenate([[0], Lorenz_init, [1]])\n",
    "    wealth_percentile = np.concatenate([[0], wealth_percentile, [1]])\n",
    "\n",
    "    # Create a list of wealth_percentile 1D array and Lorenz Shares 1D array\n",
    "    Lorenz = np.stack((wealth_percentile, Lorenz_init))\n",
    "\n",
    "    # Compute the Gini coefficient\n",
    "    Gini = 1.0 - 2.0 * np.mean(Lorenz_init[1])\n",
    "\n",
    "    # Compute the average MPC by income quintile in the latest simulated period\n",
    "    Avg_MPC = calc_subpop_avg(\n",
    "        MPC_all,\n",
    "        pLvl_all,\n",
    "        cutoffs=[(0.0, 0.2), (0.2, 0.4), (0.4, 0.6), (0.6, 0.8), (0.8, 1.0)],\n",
    "    )\n",
    "\n",
    "    return AggWealthRatio, Lorenz, Gini, Avg_MPC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The user needs to call the function $\\texttt{runRoszypalSchlaffmanExperiment}$ with specific values for $\\texttt{CorrAct}$, $\\texttt{CorrPcvd}$ and estimates of $\\texttt{DiscFac\\_center}$, $\\texttt{DiscFac\\_spread}$ to solve the model accordingly.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Lorenz curve for assests is\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG2CAYAAACTTOmSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXAElEQVR4nO3deVxU5eIG8GfY91H2fRHcURRQFK+i5l6aVmpa7i1eKxey0muu1c9WM0ut3NpMbTPtSiaZKIorgqigoiCIDiCLMyyyzMz5/YHOjUBlcIYzMM/38+GTnDkz8+BJ5/Gc97yvRBAEAURERERGyETsAERERERiYREiIiIio8UiREREREaLRYiIiIiMFosQERERGS0WISIiIjJaLEJERERktFiEiIiIyGixCBEREZHRYhEiIiIioyVqETp06BBGjhwJT09PSCQS/Prrrw98zsGDBxEWFgYrKyu0adMGn3/+uf6DEhERUYskahEqKytDSEgIPvvsswbtn5mZiREjRqBv375ISkrCf/7zH8yePRs///yznpMSERFRSyQxlEVXJRIJdu7cidGjR99znzfeeAO7d+9GWlqaZtvMmTNx5swZHD16tAlSEhERUUtiJnYAbRw9ehRDhgyptW3o0KHYtGkTqqurYW5uXuc5lZWVqKys1HyvVqtRVFQEJycnSCQSvWcmIiKihycIAkpKSuDp6QkTE91d0GpWRSg3Nxdubm61trm5uUGpVKKgoAAeHh51nrNy5UosX768qSISERGRHl27dg3e3t46e71mVYQA1DmLc/fK3r3O7ixcuBDR0dGa7+VyOXx9fXHt2jU4ODjoLygRERHpzMCVvyPx3adhb2+v09dtVkXI3d0dubm5tbbl5+fDzMwMTk5O9T7H0tISlpaWdbY7ODiwCBERETUDt6tUyK+ouRym62EtzWoeod69eyM2NrbWtn379iE8PLze8UFERETU/KXKFHp7bVGLUGlpKZKTk5GcnAyg5vb45ORkZGdnA6i5rDV58mTN/jNnzkRWVhaio6ORlpaGzZs3Y9OmTZg/f74Y8YmIiKgJJGYV6e21Rb00durUKQwYMEDz/d2xPFOmTMFXX30FmUymKUUAEBAQgJiYGMybNw9r166Fp6cn1qxZgyeffLLJsxMREVHTOHW1WG+vbTDzCDUVhUIBqVQKuVzOMUJEREQGThAEhL/9J24W3cK11eN0/vndrMYIERERkXG5WliOwrIqmJvqp7KwCBEREZHBOnm1ZnxQsKd+ruKwCBEREZHBSrwzPqibXyu9vD6LEBERERmsU3fuGOvu01ovr88iRERERAapuKwKV26WAQC6+bTSy3uwCBEREZFBSsyquSzWxsUWjrYWenkPFiEiIiIySKfuFKFwP/1cFgNYhIiIiMhAHbp0EwAQEVD/eqK6wCJEREREBidXXoFUmQISCRDV3kVv78MiRERERAYn7mI+AKCrdys421nq7X1YhIiIiMjgHLhThAa2d9Xr+7AIERERkUGpVKpwOL0AADCgg/4uiwEsQkRERGRgTl0tRlmVCs52lgj2lOr1vViEiIiIyKAcuFBzWax/exeYmEj0+l4sQkRERGRQ/ro7PqiDfscHASxCREREZEAyC8qQcbMMZiYS/Kuts97fj0WIiIiIDMbv52QAgN6BTnCwMtf7+7EIERERkcHYey4XADAs2L1J3o9FiIiIiAzC9Vu3kZIjh0QCDOnEIkRERERG5O7ZoB7+jnCx199s0n/HIkREREQGYe+d8UHDm+iyGMAiRERERAYgv6QCp7KKATTd+CCARYiIiIgMwB/n8yAIQDefVvCQWjfZ+7IIERERkej+e+YGgKa9LAawCBEREZHIbty6jRNXiwAAI0M8m/S9WYSIiIhIVL+duQFBAHoGOMKzVdNdFgNYhIiIiEhku5JrLouN7ubV5O/NIkRERESiSc8rQapMAXNTSZOPDwJYhIiIiEhEu+8Mko5q54LWthZN/v4sQkRERCQKQRA0l8VGiXBZDGARIiIiIpEkZhUju6gcNhamGNTRVZQMLEJEREQkih9P5QAAHu3iARsLM1EysAgRERFRkyuvUmLP2Zq1xZ4K8xYtB4sQERERNbk/zueitFIJX0cb9AxwFC0HixARERE1uZ8Say6LPRXmDYlEIloOFiEiIiJqUjnF5Ui4UggAeCJUnLvF7mIRIiIioib1c+J1CAIQGegE79Y2omZhESIiIqImo1IL2HEyGwAwLtxH5DQsQkRERNSEDl7Kxw15BVrbmGOYCEtq/BOLEBERETWZ74/XnA16MtQbVuamIqdhESIiIqImcuPWbfx1IR8A8HRPX5HT1GARIiIioibxw6lrUAtARIAjglztxI4DgEWIiIiImoBSpcaOk9cAABMjDONsEMAiRERERE3gz7R8yOQVcLS1MIhB0nexCBEREZHefXP0KgBgQk8fWJqJP0j6LhYhIiIi0qv0vBIkXCmEiQSYGOEndpxaWISIiIhIr745mgUAGNzJDV6trEVOUxuLEBEREemNoqIaP5+uWWB1Sm9/ccPUg0WIiIiI9ObnxByUV6kQ5GqH3oFOYsepg0WIiIiI9EKlFrD5SCYAYEqkPyQSiciJ6mIRIiIiIr3443wurhXdRmsbczwV6i12nHqxCBEREZFebIjPAAA828sP1haGc8v837EIERERkc4lZhUhKfsWLExNMKm3Yd0y/3csQkRERKRzGw7VjA0a3d0TrvZWIqe5NxYhIiIi0qmswjL8kZoLAHiubxuR09wfixARERHp1JYjVyEIQFQ7F7Rzsxc7zn2xCBEREZHOyMur8cOpmlXmnzfws0EAixARERHp0NYTWSivUqGDuz36BBneBIr/xCJEREREOlGlVOPrhKsAas4GGeIEiv/EIkREREQ6sSv5OvIUlXBzsMTIEE+x4zQIixARERE9NJVawPq4KwCAaX0CYGHWPCpG80hJREREBm3PWRkyCsrQysYcz/Yy3AkU/4lFiIiIiB6KWi3gs7/SAQDT+wTAztJM5EQNJ3oRWrduHQICAmBlZYWwsDDEx8ffd/+tW7ciJCQENjY28PDwwLRp01BYWNhEaYmIiOif9qXm4lJeKewtzTAl0l/sOFoRtQjt2LEDc+fOxaJFi5CUlIS+ffti+PDhyM7Ornf/w4cPY/LkyZgxYwbOnz+PH3/8ESdPnsRzzz3XxMmJiIgIAARBwKd/XQYATO3jD6m1uciJtCNqEVq1ahVmzJiB5557Dh07dsTq1avh4+OD9evX17v/sWPH4O/vj9mzZyMgIAD/+te/8OKLL+LUqVNNnJyIiIgA4MDFfJy/oYCNhSmm9QkQO47WRCtCVVVVSExMxJAhQ2ptHzJkCBISEup9TmRkJHJychATEwNBEJCXl4effvoJjz766D3fp7KyEgqFotYXERERPTxBELBmf83ZoEm9/OBoayFyIu2JVoQKCgqgUqng5uZWa7ubmxtyc3PrfU5kZCS2bt2K8ePHw8LCAu7u7mjVqhU+/fTTe77PypUrIZVKNV8+Pj46/TmIiIiM1ZHLhUi+dguWZiYGv7jqvYg+WPqfs04KgnDPmShTU1Mxe/ZsLFmyBImJidi7dy8yMzMxc+bMe77+woULIZfLNV/Xrl3TaX4iIiJjtebOnWITevrCxd5S5DSNI9r9bc7OzjA1Na1z9ic/P7/OWaK7Vq5ciT59+uC1114DAHTt2hW2trbo27cv3n77bXh4eNR5jqWlJSwtm+fBISIiMlTHMwpxIrMIFqYmmBkVKHacRhPtjJCFhQXCwsIQGxtba3tsbCwiIyPrfU55eTlMTGpHNjU1BVBzJomIiIj0TxAEfLjvIgBgXA9vuEutRE7UeKJeGouOjsbGjRuxefNmpKWlYd68ecjOztZc6lq4cCEmT56s2X/kyJH45ZdfsH79emRkZODIkSOYPXs2evbsCU/P5rGmCRERUXMXd+kmTl4thqWZCV4e0FbsOA9F1Kkfx48fj8LCQqxYsQIymQzBwcGIiYmBn1/N1NwymazWnEJTp05FSUkJPvvsM7z66qto1aoVBg4ciPfee0+sH4GIiMioqNUCPvyj5mzQlEj/Zn02CAAkgpFdU1IoFJBKpZDL5XBwcBA7DhERUbOyJ0WGl74/DTtLMxx6fUCT3TKvr89v0e8aIyIiouZBqVLjo9ias0HP9Q1olvMG/ROLEBERETXIL6evI+NmGRxtLZrtvEH/xCJERERED1SpVGH1n5cAALP6BzarFebvh0WIiIiIHuj749m4Ia+Au4MVnu3lJ3YcnWERIiIiovsqq1Ri7YGaNcVmP9IWVuamIifSHRYhIiIiuq8N8RkoKK2Cn5MNxoZ7ix1Hp1iEiIiI6J7yFRX44mAGAOC1oe1hbtqyqkPL+mmIiIhIpz7adwm3q1UI9W2FR7vUXdOzuWMRIiIionqlyRT4IfEaAGDRo50gkUhETqR7LEJERERUr/+LSYMgAI928UCYX2ux4+gFixARERHVEXcxH/HpBbAwNcEbwzqIHUdvWISIiIioFqVKjf+LSQMATIn0g6+TjciJ9IdFiIiIiGr5MTEHl/JK0crGHC8PaCt2HL1iESIiIiKNskolPtpXs5TG7IFtIbUxFzmRfrEIERERkcZnBy6joLQS/k42LWopjXthESIiIiIAQGZBGTbG10ye+OajnWBh1vJrQsv/CYmIiKhB3vpvKqpVAvq3d8EjHV3FjtMkWISIiIgIf13Iw18X8mFuKsGSx1rm5In1YREiIiIycpVKFVb8lgoAmP6vALRxsRM5UdNpVBG6desWNm7ciIULF6KoqAgAcPr0aVy/fl2n4YiIiEj/Nh3OxNXCcrjaW+KVgS37dvl/MtP2CSkpKRg0aBCkUimuXr2K559/Ho6Ojti5cyeysrLwzTff6CMnERER6UGuvAKf/XUZALBwRAfYWWpdDZo1rc8IRUdHY+rUqUhPT4eVlZVm+/Dhw3Ho0CGdhiMiIiL9Wvl7GsqrVAjza43R3bzEjtPktC5CJ0+exIsvvlhnu5eXF3Jzc3USioiIiPTvWEYhdiXfgEQCLB/V2WgGSP+d1kXIysoKCoWizvaLFy/CxcVFJ6GIiIhIv6qUarz56zkAwISevgj2koqcSBxaF6HHH38cK1asQHV1NQBAIpEgOzsbCxYswJNPPqnzgERERKR7G+IzcDm/FM52FnhjaMtdXf5BtC5CH374IW7evAlXV1fcvn0bUVFRCAoKgr29Pd555x19ZCQiIiIdulZUjk//SgcALHq0Y4tfT+x+tB4a7uDggMOHD+Ovv/7C6dOnoVarERoaikGDBukjHxEREemQIAhYsuscKqrV6N3GySgHSP+dVkVIqVTCysoKycnJGDhwIAYOHKivXERERKQHf5zPxYGLN2FuKsFbo4ONcoD032l1aczMzAx+fn5QqVT6ykNERER6UlqpxPI7M0jPjApEkKvxzCB9L1qPEXrzzTdrzShNREREzcPq2EuQySvg62iDlwYEiR3HIGg9RmjNmjW4fPkyPD094efnB1tb21qPnz59WmfhiIiISDfOXZdjS8JVAMCKxzvDytxU3EAGQusiNHr0aD3EICIiIn2pVqnx+k8pUKkFPNrVA/3bu4odyWBoXYSWLl2qjxxERESkJ18eykCqTIFWNuZYNrKz2HEMSqNWnyciIqLm4XJ+KT7ZXzNn0JLHOsHF3lLkRIZF6zNCKpUKH3/8MX744QdkZ2ejqqqq1uMcRE1ERGQY1GoBC35OQZVSjah2LhjT3bjnDKqP1meEli9fjlWrVmHcuHGQy+WIjo7GE088ARMTEyxbtkwPEYmIiKgxvjuehVNZxbC1MMU7YzhnUH20LkJbt27Fhg0bMH/+fJiZmWHChAnYuHEjlixZgmPHjukjIxEREWkpp7gc7/1+AQDwxvAO8G5tI3Iiw6R1EcrNzUWXLl0AAHZ2dpDL5QCAxx57DHv27NFtOiIiItKaIAhYtPMcyqpUCPdrjWcj/MSOZLC0LkLe3t6QyWQAgKCgIOzbtw8AcPLkSVhacgAWERGR2H44dQ0HL92EhZkJ3nuqK0xMeEnsXrQuQmPGjMH+/fsBAHPmzMHixYvRtm1bTJ48GdOnT9d5QCIiImq4a0XlWHFnGY1XB7dDoAuX0bgfiSAIwsO8wLFjx5CQkICgoCCMGjVKV7n0RqFQQCqVQi6Xw8HBQew4REREOqNWC5i48RiOZRShh39rbH+hN0xbyNkgfX1+a337/D/16tULvXr10kUWIiIieghfJVzFsYwiWJub4sOxIS2mBOlTo4rQpUuXEBcXh/z8fKjV6lqPLVmyRCfBiIiIqOGu3CzFe3tr7hL7z6Md4edk+4BnENCIIrRhwwb8+9//hrOzM9zd3WvNSSCRSFiEiIiImphSpUb0D2dQqVSjb1tnPBvhK3akZkPrIvT222/jnXfewRtvvKGPPERERKSlLw5l4My1W7C3MsP7T3XlxIla0PquseLiYowdO1YfWYiIiEhLZ3PkWP3nJQDA8lGd4SG1FjlR86J1ERo7dqxm7iAiIiIST1mlErO3J6FaJWBYZ3euJdYIDbo0tmbNGs2vg4KCsHjxYhw7dgxdunSBubl5rX1nz56t24RERERUrxW/pSKzoAzuDlZ498kuvCTWCA2aRyggIKBhLyaRICMj46FD6RPnESIiopYg5qwMs7aehkQCfP9cL/QOdBI7kl6JOo9QZmamzt6QiIiIHs6NW7ex4OcUAMC/owJbfAnSJ63HCP2TSqVCcnIyiouLdZGHiIiI7kOlFjBvRzIUFUqEeEsxb3A7sSM1a1oXoblz52LTpk0AakpQv379EBoaCh8fH8TFxek6HxEREf3N5wev4HhmEWwtTPHJ091hbvrQ5zSMmta/ez/99BNCQkIAAL/99huuXr2KCxcuYO7cuVi0aJHOAxIREVGNpOxirIq9c6v848Hwd+bs0Q9L6yJUUFAAd3d3AEBMTAzGjh2Ldu3aYcaMGTh79qzOAxIRERFQWqnEnO3JUKkFPNbVA0+G8lZ5XdC6CLm5uSE1NRUqlQp79+7FoEGDAADl5eUwNTXVeUAiIiIClu46j+yicni1ssY7Y3irvK5ovcTGtGnTMG7cOHh4eEAikWDw4MEAgOPHj6NDhw46D0hERGTsdiVfx8+nc2AiAVY/3Q1Sa/MHP4kaROsitGzZMgQHB+PatWsYO3YsLC0tAQCmpqZYsGCBzgMSEREZs8v5pVj4S83Qk1cGtkUPf0eRE7UsDZpQsSXhhIpERNRc3K5SYfTaI7iYV4LebZzw3XMRMDUxzkti+vr85j13REREBmrJrnO4mFcCF3tLfDKhm9GWIH1iESIiIjJAP5y6hh8Ta8YFrXm6O1ztrcSO1CKxCBERERmYC7kKLNl1DgAQPbgdl9DQIxYhIiIiA1JaqcSsradRUa1Gv3YumNU/SOxILZrWd40BgFqtxuXLl5Gfnw+1Wl3rsX79+ukkGBERkbERBAH/+eUsMm6Wwd3BCqvHd4MJxwXpldZnhI4dO4agoCB07NgR/fr1Q//+/TVfAwYM0DrAunXrEBAQACsrK4SFhSE+Pv6++1dWVmLRokXw8/ODpaUlAgMDsXnzZq3fl4iIyNB8fyIbu8/cgJmJBGuf6Q5HWwuxI7V4Wp8RmjlzJsLDw7Fnzx7NpIqNtWPHDsydOxfr1q1Dnz598MUXX2D48OFITU2Fr69vvc8ZN24c8vLysGnTJgQFBSE/Px9KpbLRGYiIiAzBuetyLN+dCgB4Y1gHhPlxvqCmoPU8Qra2tjhz5gyCgh7+mmVERARCQ0Oxfv16zbaOHTti9OjRWLlyZZ399+7di6effhoZGRlwdGzc/yCcR4iIiAzNrfIqjPrsCLKLyjGooxs2TA7jEhr/YDDzCEVERODy5csP/cZVVVVITEzEkCFDam0fMmQIEhIS6n3O7t27ER4ejvfffx9eXl5o164d5s+fj9u3b9/zfSorK6FQKGp9ERERGQqVWsDs7cnILiqHj6M1PhobwhLUhBp0aSwlJUXz61deeQWvvvoqcnNz0aVLF5ib117vpGvXrg1644KCAqhUKri5udXa7ubmhtzc3Hqfk5GRgcOHD8PKygo7d+5EQUEBZs2ahaKionuOE1q5ciWWL1/eoExERERNbVXsRRy6dBNW5ib44tlwSG24jlhTalAR6tatGyQSCf5+FW369OmaX999TCKRQKVSaRXgn6337uvUR61WQyKRYOvWrZBKpQCAVatW4amnnsLatWthbW1d5zkLFy5EdHS05nuFQgEfHx+tMhIREenD3nMyrD1wBQDw3pNd0cmTQzaaWoOKUGZmps7f2NnZGaampnXO/uTn59c5S3SXh4cHvLy8NCUIqBlTJAgCcnJy0LZt2zrPsbS01CwMS0REZCjS80rw6g9nAADP/SsAj3fzEjmRcWrQGCE/Pz/NV1ZWFry8vGpt8/Pzg5eXF7Kyshr8xhYWFggLC0NsbGyt7bGxsYiMjKz3OX369MGNGzdQWlqq2Xbp0iWYmJjA29u7we9NREQkJkVFNV74NhFlVSr0buOEBcM7iB3JaGk9WHrAgAEoKiqqs10ul2s9j1B0dDQ2btyIzZs3Iy0tDfPmzUN2djZmzpwJoOay1uTJkzX7T5w4EU5OTpg2bRpSU1Nx6NAhvPbaa5g+fXq9l8WIiIgMjVotIHpHMjILyuAptcJnE7vDzJQLPYhF63mE7jWGp7CwELa2tlq91vjx41FYWIgVK1ZAJpMhODgYMTEx8PPzAwDIZDJkZ2dr9rezs0NsbCxeeeUVhIeHw8nJCePGjcPbb7+t7Y9BREQkijV/pePPtHxYmJngi0nhcLLj8A0xNXgeoSeeeAIAsGvXLgwbNqzWuBuVSoWUlBS0b98ee/fu1U9SHeE8QkREJJb9aXmY8fUpAMAHT3XF2HDevNNQ+vr8bvAZobsDlAVBgL29fa1LURYWFujVqxeef/55nQUjIiJqSa7cLMXc7ckAgCm9/ViCDESDi9CWLVsAAP7+/pg/f77Wl8GIiIiM1a3yKsz46iRKKpXo4d8abz7WSexIdIfWY4SWLl2qjxxEREQtUrVKjX9/dxpXC8vh1coa658NgzkHRxuMBhWh7t27N3i679OnTz9UICIiopZCEAQs2XUeRzMKYWthik1Tw+HMwdEGpUFFaPTo0XqOQURE1PJ8lXAV205kQyIB1kzojg7uvEnH0DSoCPFyGBERkXbiLubjrf+mAgD+M7wjHulY/6oJJC5epCQiItKx9LwSvPJ9EtQCMC7cG8/1DRA7Et2D1oOlVSoVPv74Y/zwww/Izs5GVVVVrcfrm3WaiIjIWBSVVWHG16dQUqlEzwBHvD26S4PH2VLT0/qM0PLly7Fq1SqMGzcOcrkc0dHReOKJJ2BiYoJly5bpISIREVHzUKVUY+Z3icguKoePozU+fzYMFma8+GLItD46W7duxYYNGzB//nyYmZlhwoQJ2LhxI5YsWYJjx47pIyMREZHBEwQBi3aexYnMIthbmmHzlB5wtLUQOxY9gNZFKDc3F126dAFQs/aXXC4HADz22GPYs2ePbtMRERE1E2v2X8aPiTkwkQBrJnZHWzd7sSNRA2hdhLy9vSGTyQAAQUFB2LdvHwDg5MmTtdYfIyIiMhY/Jebg4z8vAQDeGh2MAe1dRU5EDaV1ERozZgz2798PAJgzZw4WL16Mtm3bYvLkyZg+fbrOAxIRERmyw+kFWPBzCgDg3/0D8UyEn8iJSBsNXn3+Xo4dO4aEhAQEBQVh1KhRusqlN1x9noiIdCVNpsDYz4+itFKJx7t54uNx3WBiwjvE9EH01efvpVevXujVq5cushARETUbMvltTNtyEqWVSkQEOOL9p7qyBDVDjbqn79tvv0WfPn3g6emJrKwsAMDq1auxa9cunYYjIiIyRIqKakzbchK5igoEudrhy0nhsDQzFTsWNYLWRWj9+vWIjo7GiBEjcOvWLahUKgBAq1atsHr1al3nIyIiMijVKjVmfXcaF3JL4GJvia+m9YDUxlzsWNRIWhehTz/9FBs2bMCiRYtgavq/9hseHo6zZ8/qNBwREZEhEQQBC385i8OXC2BjYYotU3vAu7WN2LHoIWhdhDIzM9G9e/c62y0tLVFWVqaTUERERIbo/T8u4qfEHJiaSLD2mVAEe0nFjkQPSesiFBAQgOTk5Drbf//9d3Tq1EkXmYiIiAzOhkMZWB93BQDwf2M4V1BLofVdY6+99hpeeuklVFRUQBAEnDhxAtu2bcPKlSuxceNGfWQkIiIS1U+JOXgnJg0A8MawDhjfw1fkRKQrWhehadOmQalU4vXXX0d5eTkmTpwILy8vfPLJJ3j66af1kZGIiEg0f6bm4Y07EyY+3zcAM6PaiJyIdOmhJlQsKCiAWq2Gq2vzOT3ICRWJiKihjmcUYvLmE6hUqvFkqDc+HNsVEgnnChKDvj6/tR4jtGHDBqSnpwMAnJ2dm1UJIiIiaqjzN+R47utTqFSqMaijG957sgtLUAukdRH66KOP0L59e3h6emLChAn44osvcOHCBX1kIyIiEkVWYRmmbD6JkkolegY44rOJ3WFm2qg5iMnAaX1UL1y4gBs3buCjjz6CVCrFxx9/jM6dO8Pd3Z1jhIiIqNnLV1Tg2U3HUVBaiY4eDtg4JRxW5pw1uqV6qDFCZWVlOHz4MLZv347vvvsOgiBAqVTqMp/OcYwQERHdS2FpJZ7+8hjS80vh52SDH2f2hqu9ldixCAa06Orvv/+OgwcPIi4uDmfOnEHnzp3Rr18//Pzzz+jbt6/OghERETWlW+VVmLTpBNLzS+HuYIVvp0ewBBkBrYvQo48+ChcXF7z66qv4448/IJVyVk0iImreFBXVmLL5BFJlCjjbWWLr8xHwdeLSGcZA6zFCq1atQp8+ffDBBx+gffv2GD9+PNavX4+0tDR95CMiItKrskolpm85iTM5crS2McfW5yIQ6GIndixqIg81Rujs2bM4ePAgDhw4gN9++w1OTk6QyWS6zKdzHCNERER3VVSrMG3LSRzNKISDlRm+f74X1w8zUAYzRuiupKQkxMXF4cCBA4iPj4darYa3t7fOghEREelTpVKFF75NxNGMQthZmuGbGREsQUZI60tjo0aNgqOjI3r06IGtW7eiXbt2+Pbbb1FUVISTJ0/qIyMREZFOVavUeGlrEg5duglrc1NsmdYD3XxaiR2LRKD1GaF27drhhRdeQL9+/XhpiYiImp1qlRpztyfjz7Q8WJqZYNOUcPTwdxQ7FolE6yL04Ycf6iMHERGR3lWr1JizPQkxZ3NhYWqCLyaFITLIWexYJKJGjxEiIiJqTqpVarzyfRL2nq8pQeufDUX/9lwv09ixCBERUYtXpVTj5e9PY19qnuZM0IAOLEHEIkRERC1clVKNWVtP48+0PFiYmeDLSWE8E0QaLEJERNRiVSpVmPXdaey/kA9LMxNsmByOfu1cxI5FBkTr2+cB4MqVK3jzzTcxYcIE5OfnAwD27t2L8+fP6zQcERFRY1VUqzDz20RNCdo0pQdLENWhdRE6ePAgunTpguPHj+OXX35BaWkpACAlJQVLly7VeUAiIiJtVVSrMPO7RBy4eBNW5ibYPLUH/tWWd4dRXVoXoQULFuDtt99GbGwsLCwsNNsHDBiAo0eP6jQcERGRtsqrlJjx9UnEXayZLHHz1B7ow1vk6R60LkJnz57FmDFj6mx3cXFBYWGhTkIRERE1hqKiGpM3ncCRy4WwtaiZMToykCWI7k3rItSqVat6F1ZNSkqCl5eXTkIRERFpq6isChM3HMOprGI4WJnhu+ci0KuNk9ixyMBpXYQmTpyIN954A7m5uZBIJFCr1Thy5Ajmz5+PyZMn6yMjERHRfeUrKvD0l0dx7roCTrYW2P5Cb3T3bS12LGoGtC5C77zzDnx9feHl5YXS0lJ06tQJ/fr1Q2RkJN588019ZCQiIrqnnOJyjPviKC7llcLNwRI7XuyNTp5cC5MaRiIIgtCYJ165cgVJSUlQq9Xo3r072rZtq+tseqFQKCCVSiGXy7loLBFRM5dZUIZnNhzDDXkFvFtb4/vnesHXyUbsWKQH+vr81npCxYMHDyIqKgqBgYEIDAzUWRAiIiJtXMorwTMbj+NmSSXauNhi63MR8JBaix2LmhmtL40NHjwYvr6+WLBgAc6dO6ePTERERPd1OrsY4744ipsllejgbo8dL/RmCaJG0boI3bhxA6+//jri4+PRtWtXdO3aFe+//z5ycnL0kY+IiKiWg5du4pkNx3GrvBrdfFph+wu94GJvKXYsaqYaPUYIADIzM/H9999j27ZtuHDhAvr164e//vpLl/l0jmOEiIiar91nbuDVH5JRrRLQt60zPn82DLaWXDbTGOjr8/uhihAAqFQq/P7771i8eDFSUlKgUql0lU0vWISIiJqnb45exdLd5yEIwMgQT3w0NgQWZo1aMpOaIX19fjf6/6AjR45g1qxZ8PDwwMSJE9G5c2f897//1VkwIiIiABAEAR/HXsKSXTUlaHJvP3wyvhtLEOmE1ucT//Of/2Dbtm24ceMGBg0ahNWrV2P06NGwseHtikREpFsqtYBlu8/j22NZAIC5g9piziNtIZFIRE5GLYXWRSguLg7z58/H+PHj4ezM9VuIiEg/qpRqRP+QjP+myCCRACtGdcak3v5ix6IWRusilJCQoI8cREREGoqKasz8NhEJVwphbirBqnHdMDLEU+xY1AI1qAjt3r0bw4cPh7m5OXbv3n3ffUeNGqWTYEREZJxy5RWYuuUELuSWwNbCFJ9PCkPfti5ix6IWqkF3jZmYmCA3Nxeurq4wMbn34DSJRMK7xoiIqNEu5ZVg6uYTuCGvgIu9JbZM7YFgL6nYscgAiLrEhlqtrvfXREREunI8oxDPf3MKigol2rjY4utpPeHjyBtxSL+0vvfwm2++QWVlZZ3tVVVV+Oabb3QSioiIjMueFBkmbToBRYUSYX6t8fPMSJYgahJaT6hoamoKmUwGV1fXWtsLCwvh6urKS2NERKSVzYcz8daeVAgCMLSzGz55ujuszE3FjkUGxmBWnxcEod75G3JyciCV8jouERE1jEotYGVMGjYezgRQM1Hi0pGdYWrCOYKo6TS4CHXv3h0SiQQSiQSPPPIIzMz+91SVSoXMzEwMGzZMLyGJiKhlKa9SYs72ZMSm5gEA3hjWATOj2nCiRGpyDS5Co0ePBgAkJydj6NChsLOz0zxmYWEBf39/PPnkk1oHWLduHT744APIZDJ07twZq1evRt++fR/4vCNHjiAqKgrBwcFITk7W+n2JiEgceYoKPPf1KZy9LoeFmQk+HBuCUZwjiETS4CK0dOlSAIC/vz/Gjx8PKyurh37zHTt2YO7cuVi3bh369OmDL774AsOHD0dqaip8fX3v+Ty5XI7JkyfjkUceQV5e3kPnICKippF6Q4EZX5+ETF4BR1sLbJgchjA/R7FjkRF76NXnH0ZERARCQ0Oxfv16zbaOHTti9OjRWLly5T2f9/TTT6Nt27YwNTXFr7/+qtUZIQ6WJiISx4EL+Xj5+9Moq1Ih0MUWW6b2hK8T7wyjhjGY1edVKhU+/PBD9OzZE+7u7nB0dKz11VBVVVVITEzEkCFDam0fMmTIfZfx2LJlC65cuaI5Q/UglZWVUCgUtb6IiKhpfZ1wFTO+PomyKhUiA53wy6w+LEFkELQuQsuXL8eqVaswbtw4yOVyREdH44knnoCJiQmWLVvW4NcpKCiASqWCm5tbre1ubm7Izc2t9znp6elYsGABtm7dWmuw9v2sXLkSUqlU8+Xj49PgjERE9HCUKjWW7T6PpbvPQy0A48N98PX0npBam4sdjQhAI4rQ1q1bsWHDBsyfPx9mZmaYMGECNm7ciCVLluDYsWNaB/jnHQL3uj1fpVJh4sSJWL58Odq1a9fg11+4cCHkcrnm69q1a1pnJCIi7cnLqzHtq5P4KuEqgJo7w959sgvMTbX+6CHSG63nEcrNzUWXLl0AAHZ2dpDL5QCAxx57DIsXL27w6zg7O8PU1LTO2Z/8/Pw6Z4kAoKSkBKdOnUJSUhJefvllADXLfQiCADMzM+zbtw8DBw6s8zxLS0tYWlo2OBcRET28y/mleP6bU8gsKIO1uSlWjQvB8C4eYsciqkPrWu7t7Q2ZTAYACAoKwr59+wAAJ0+e1KpwWFhYICwsDLGxsbW2x8bGIjIyss7+Dg4OOHv2LJKTkzVfM2fORPv27ZGcnIyIiAhtfxQiItKDAxfzMWbtEWQWlMGrlTV+/nckSxAZLK3PCI0ZMwb79+9HREQE5syZgwkTJmDTpk3Izs7GvHnztHqt6OhoTJo0CeHh4ejduze+/PJLZGdnY+bMmQBqLmtdv34d33zzDUxMTBAcHFzr+a6urrCysqqznYiImp4gCNgQn4GVv1+AIAA9/R2x7tlQONvxrDwZLq2L0Lvvvqv59VNPPQVvb28kJCQgKCgIo0aN0uq1xo8fj8LCQqxYsQIymQzBwcGIiYmBn58fAEAmkyE7O1vbiERE1MQqqlX4zy9n8UvSdQDAhJ4+WD4qGBZmHA9Ehk3UeYTEwHmEiIh0K19RgRe+TUTytVswNZFgyWOdMLm3H5fLIJ0SddHV3bt3N/gFtT0rREREzdeZa7fwwrenkKeohNTaHOueCUWfIGexYxE1WIOK0N11xh5EIpFApVI9TB4iImomdiVfx+s/paBSqUZbVztsnBIOPydbsWMRaaVBRUitVus7BxERNRPVKjXe/f0CNh3OBAA80sEVq5/uBnsrTpJIzY/Wg6WJiMh43SypxMvfn8bxzCIAwEsDAhE9uD1MTTgeiJonrYvQihUr7vv4kiVLGh2GiIgM1+nsYsz67jRyFRWwszTDh2NDMCzYXexYRA9F6yK0c+fOWt9XV1cjMzMTZmZmCAwMZBEiImphBEHA9yeysWz3eVSrBAS62OKLSeEIcrUTOxrRQ9O6CCUlJdXZplAoMHXqVIwZM0YnoYiIyDBUVKuwZNc5/HAqBwAwPNgdH4wNgZ0lR1ZQy6CzeYTOnTuHxx57DFevXtXFy+kN5xEiImqYnOJy/Pu70zh7XQ4TCfD6sA54sV8bzg9EohB1HqGGuHXrlmYBViIiat4OpxfglW2nUVxejdY25vhsIucHopZJ6yK0Zs2aWt8LggCZTIZvv/0Ww4YN01kwIiJqemq1gM8PXcGHf1yEWgC6eEmx/tlQeLe2ETsakV5oXYQ+/vjjWt+bmJjAxcUFU6ZMwcKFC3UWjIiImlZxWRVe/fEM/rqQDwAYG+aNt0YHw8rcVORkRPqjdRHKzMzURw4iIhJRYlYxXvn+NG7IK2BhZoLlozrj6R4+HA9ELR6H/RMRGTFBELAxPhPv7b0ApVpAgLMt1k4MRSdP3kxCxkHrIlRRUYFPP/0UBw4cQH5+fp3lN06fPq2zcEREpD/y8mq8+uMZ/JmWBwB4rKsHVj7RhUtlkFHRughNnz4dsbGxeOqpp9CzZ0+eNiUiaoaSr93CS1tP4/qt27AwNcHikZ3wbIQv/04no6N1EdqzZw9iYmLQp08ffeQhIiI9EgQBXyVcxf/FpKFaJcDX0QbrnglFsJdU7GhEotC6CHl5ecHe3l4fWYiISI/kt6vxxk8p2Hs+F0DNLNHvPdUVDrwURkbMRNsnfPTRR3jjjTeQlZWljzxERKQHiVnFeOzTeOw9nwtzUwmWjeyEdc+EsgSR0dP6jFB4eDgqKirQpk0b2NjYwNy89h+ioqIinYUjIqKHo1ILWHfgMlbvT4dKLcC7tTXWTgxFiE8rsaMRGQSti9CECRNw/fp1/N///R/c3Nw4sI6IyEDduHUbc3ck40RmzT9QH+/mibdGB/MsENHfaF2EEhIScPToUYSEhOgjDxER6UDMWRkW/JwCRYUSthameGt0MJ4I9RY7FpHB0boIdejQAbdv39ZHFiIiekjlVUqs+C0V209eAwCE+LTCmqe7wc/JVuRkRIZJ6yL07rvv4tVXX8U777yDLl261Bkj5ODA2UiJiMRw7rocs7cnIeNmGSQSYFb/QMwd1A7mplrfF0NkNCSCIAjaPMHEpOYP1D/HBgmCAIlEApVKpbt0eqBQKCCVSiGXy1naiKhFUKsFbD5Ss0xGtUqAu4MVVo0PQWSgs9jRiHRGX5/fWp8ROnDggM7enIiIHo5Mfhuv/5SC+PQCAMCQTm5478muaG1rIXIyouZB6yIUFRWljxxERKQFQRCw+8wNLP71HBQVSliZm2DxY50wsSeXySDShtZF6NChQ/d9vF+/fo0OQ0RED3arvAqLfj2HPSkyAECItxSrxndDoIudyMmImh+ti1D//v3rbPv7vz4MfYwQEVFzFncxH6//lIL8kkqYmkgwe2BbvDQgEGYcEE3UKFoXoeLi4lrfV1dXIykpCYsXL8Y777yjs2BERPQ/5VVK/F9MGr47lg0ACHSxxcfju6GrdytxgxE1c1oXIam07grFgwcPhqWlJebNm4fExESdBCMiohqns4sRvSMZVwvLAQBTI/2xYHgHWJmbipyMqPnTugjdi4uLCy5evKirlyMiMnpVSjXW7E/HurjLUAuAh9QKHzwVgn+15W3xRLqidRFKSUmp9b0gCJDJZHj33Xe57AYRkY6cuy7Haz+lIE2mAACM6e6FZaM6Q2rNdcKIdEnrItStWzdIJBL8cx7GXr16YfPmzToLRkRkjKqUanz2VzrWxV2BUi2gtY053h7dBY929RA7GlGLpHURyszMrPW9iYkJXFxcYGVlpbNQRETGKCXnFl77MQUX80oAACO6uGP5qGC42FuKnIyo5dK6CPn5+ekjBxGR0aqoVmHN/nR8cSgDKrUAJ1sLrHg8mGeBiJpAgyee+Ouvv9CpUycoFIo6j8nlcnTu3Bnx8fE6DUdE1NIlZRfjsU8PY13cFajUAkaGeGLfvH4sQURNpMFnhFavXo3nn3++3oXOpFIpXnzxRaxatQp9+/bVaUAiopaoolqFj2MvYUN8BtQC4GxnibdHB2NYsLvY0YiMSoPPCJ05cwbDhg275+NDhgzhHEJERA1w8moRRnwSjy8O1ZSg0d08ETuvH0sQkQgafEYoLy8P5ub3vm3TzMwMN2/e1EkoIqKWSFFRjfd+v4Ctx2tmh3a1t8Q7Y7pgcCc3kZMRGa8GFyEvLy+cPXsWQUFB9T6ekpICDw9e0yYiqs/eczIs2XUe+SWVAIBx4d5YNKITpDacF4hITA0uQiNGjMCSJUswfPjwOrfK3759G0uXLsVjjz2m84BERM1ZrrwCS3adw77UPABAgLMt3hkTjMhAzg5NZAgkwj9nRryHvLw8hIaGwtTUFC+//DLat28PiUSCtLQ0rF27FiqVCqdPn4abm2Gf4lUoFJBKpZDL5fUO/CYi0gW1WsDW41l4b+9FlFYqYWYiwYtRbfDKwLZcI4yoEfT1+d3gM0Jubm5ISEjAv//9byxcuFAzs7REIsHQoUOxbt06gy9BRERNIT2vBAt+OYvErGIAQDefVnj3yS7o4M5/fBEZGq0mVPTz80NMTAyKi4tx+fJlCIKAtm3bonXr1vrKR0TUbFRUq7Au7grWx11GtUqArYUpXhvaHpN6+8PURCJ2PCKqR6NWn2/dujV69Oih6yxERM3WoUs3sXT3eWQWlAEAHungirdGB8OzlbXIyYjofhpVhIiIqEauvAJv7UnFnhQZgJpb4peM7IRHu3hAIuFZICJDxyJERNQISpUaXyVcxcexl1BWpYKJBJgaGYB5g9vC3oq3xBM1FyxCRERaOnW1CG/+eg4XcmtWiQ/1bYW3Rgejs6dU5GREpC0WISKiBioqq8K7v6fhh1M5AIBWNuZYOLwDxob5wISDoYmaJRYhIqIHUKsF/HDqGt7dewG3yqsBAOPDffDG8A5wtLUQOR0RPQwWISKi+zidXYxlu88jJUcOAOjgbo93xgQjzM9R5GREpAssQkRE9chXVOC9vRfx8+may2B2lmaYO6gtpkb6w8zUROR0RKQrLEJERH9TpVRjy5FMrNmfjrIqFYCaBVJfG9oBLvaWIqcjIl1jESIiuuPAxXy89VsqMu5MitjNpxWWjeqMbj6txA1GRHrDIkRERu9qQRne+m8q9l/IBwA421liwfAOeKK7F+8GI2rhWISIyGiVVFRj7YEr2Hw4E1UqNcxNJZjeJwAvDwzipIhERoJFiIiMjlKlxo5T1/Bx7CUUlFYBAKLauWDJyE4IdLETOR0RNSUWISIyKgcv3cQ7e1JxKa8UANDG2RaLHu2IgR1cuTYYkRFiESIio3AxtwTvxKTh0KWbAGpmhZ77SFs808sP5rwdnshosQgRUYt2s6QSH/95CdtPZEMtAOamEkyN9MfLA9pCasNxQETGjkWIiFqkimoVNh3OxPq4KyitVAIAhge7Y8HwDvBzshU5HREZChYhImpRVGoBvyZdx6rYS7h+6zYAoKu3FG8+2gk9A7gsBhHVxiJERC2CIAiIu3gT7+29gAu5JQAAD6kV3hjWAaNCPDkfEBHVS/QRguvWrUNAQACsrKwQFhaG+Pj4e+77yy+/YPDgwXBxcYGDgwN69+6NP/74ownTEpEhSr52CxM2HMO0r07iQm4JHKzMsGB4BxyY3x+jOSkiEd2HqEVox44dmDt3LhYtWoSkpCT07dsXw4cPR3Z2dr37Hzp0CIMHD0ZMTAwSExMxYMAAjBw5EklJSU2cnIgMQWZBGWZtTcTotUdwLKMIFmYmeLFfGxx6fQBmRgXCytxU7IhEZOAkgiAIYr15REQEQkNDsX79es22jh07YvTo0Vi5cmWDXqNz584YP348lixZ0qD9FQoFpFIp5HI5HBwcGpWbiMSVX1KBNfvTse3ENajUAiQS4KlQb8wb3A6erazFjkdEeqCvz2/RxghVVVUhMTERCxYsqLV9yJAhSEhIaNBrqNVqlJSUwNHx3gMgKysrUVlZqfleoVA0LjARiU5RUY0NhzKwMT4Tt6trVoZ/pIMrXh/WAe3d7UVOR0TNkWhFqKCgACqVCm5ubrW2u7m5ITc3t0Gv8dFHH6GsrAzjxo275z4rV67E8uXLHyorEYmrvEqJLUeu4stDGZDfrgYAdPdthQXDOiCijZPI6YioORP9rrF/TmkvCEKDprnftm0bli1bhl27dsHV1fWe+y1cuBDR0dGa7xUKBXx8fBofmIiaTEW1Ct8fz8a6uMuaNcHautrh1SHtMLSzO5fEIKKHJloRcnZ2hqmpaZ2zP/n5+XXOEv3Tjh07MGPGDPz4448YNGjQffe1tLSEpaXlQ+cloqZTrVLjx1M5+PSvdMjkFQAAPycbzB3UFqNCvGDKu8CISEdEK0IWFhYICwtDbGwsxowZo9keGxuLxx9//J7P27ZtG6ZPn45t27bh0UcfbYqoRNREVGoBu5KvY/Wf6cguKgdQMxfQ7Efa4qkwb64JRkQ6J+qlsejoaEyaNAnh4eHo3bs3vvzyS2RnZ2PmzJkAai5rXb9+Hd988w2AmhI0efJkfPLJJ+jVq5fmbJK1tTWkUqloPwcRPRy1WsDe87n4OPYS0vNrVoV3trPArP5BmBjhy9vgiUhvRC1C48ePR2FhIVasWAGZTIbg4GDExMTAz88PACCTyWrNKfTFF19AqVTipZdewksvvaTZPmXKFHz11VdNHZ+IHpJaLSDmnAyf7r+Mi3k1s0FLrc3xYlQbTI30h42F6MMYiaiFE3UeITFwHiEi8anUAvacleHT/emaM0D2lmaY1scfM/q2gdSaq8ITUW0tbh4hIjI+KrWA/6bcwKd/XcbluwXIygzT+wRgep8ASG1YgIioabEIEZHeKVVq/HanAGXcLAMAOFiZYca/2mBqH3+eASIi0bAIEZHeKFVq7Eq+gc8OXEZmQU0Bklqb47l/BWBKH384WLEAEZG4WISISOeq7xagv9JxtbDmNvjWNuZ4rm8bTO7tB3sWICIyECxCRKQzFdUq/HDqGr44mIHrt24DABxtLfB83zaY1NsPdpb8K4eIDAv/ViKih6aoqMa3R7Ow5UimZikMZztLPN83AM/28oMtCxARGSj+7UREjVZQWonNhzPx7dEslFQqAQDera3xYr82GBvuw4kQicjgsQgRkdZyisux4VAGtp+8hkqlGkDNYqizBgTisa6eXAqDiJoNFiEiarDL+SVYF3cFu5NvQKmumYs1xKcVXuofiEEd3WDCxVCJqJlhESKiB0rKLsbnB69gX2oe7s5F/68gZ8zqH4jegU6QSFiAiKh5YhEionqp1QL+TMvDhvgMnLxarNk+tLMbZvUPQohPK/HCERHpCIsQEdVSUa3CT4k52HQ4UzMJormpBKNCvDAzqg3autmLnJCISHdYhIgIQM0dYN8ezcK3x7JQVFZzC7yDlRme7eWHKZH+cHOwEjkhEZHusQgRGbkrN0uxMT4TP5/OQdWdO8C8W1tjxr8CMC7ch3MAEVGLxr/hiIyQIAg4kVmEDfEZ+DMtX7M9xFuKF/oFYmhnN5jxFngiMgIsQkRGpFKpwp4UGb5KuIqUHDkAQCIBBnV0w/N926CHf2veAUZERoVFiMgI5Csq8N3xbHx/PEuzBIalmQmeDPPGjH8FINDFTuSERETiYBEiaqEEQUDStVv46shVxJyVaSZAdHewwqTefni6hw+c7CxFTklEJC4WIaIWplKpQsxZGb46chVn7lz+AoBwv9aY2scfQzu7cwkMIqI7WISIWoj6Ln9ZmJpgVDdPTI30R7CXVOSERESGh0WIqBkTBAGJWcX49lgW9qTUvvz1bC9fPN3TF868/EVEdE8sQkTNkKKiGjtPX8fW41m4lFeq2c7LX0RE2mERImpGUnJuYeuxbOw+cwO3q1UAACtzE4zs6onJvf3RxZuXv4iItMEiRGTgyquU2J18A1uPZ+Ps9f8Nfm7raodnInwxJtQbUmtzERMSETVfLEJEBupCrgLfH8/GztPXUVKpBFAz+HlEF3c808sP4X6c/JCI6GGxCBEZkPIqJX4/m4ttJ7JxKqtYs93fyQYTI3zxVJgPHG0tRExIRNSysAgRiUwQBCRfu4UfTuXgtzM3UHrn7I+ZiQSDO7nhmQg/RAY6wcSEZ3+IiHSNRYhIJIWlldiZdB0/nLpW684vX0cbjAv3xrhwH7g6WImYkIio5WMRImpCKrWAQ5du4odT1/BnWh6qVTXz/liamWBEFw+MC/dBRIAjz/4QETURFiGiJpBVWIYfTl3Dz4nXkauo0GwP8ZZibLgPRnXzhIMV7/wiImpqLEJEelJSUY3fz+bi59M5OJ5ZpNne2sYcY7p7Y1wPb3RwdxAxIRERsQgR6ZBSpUZ8egF+Pp2D2NQ8VCrVAACJBOjX1gXje/jgkY6usDQzFTkpEREBLEJED00QBJy7rsAvSTV3fd1d8BQAglzt8ESoF0Z384JnK2sRUxIRUX1YhIga6cat2/g1+Tp2nr6O9Pz/3fXlZGuBUd088UR3bwR7OXDSQyIiA8YiRKSFkopq7D2Xi51J13E0oxBCzU1fsDAzwZBObngi1At927pwwVMiomaCRYjoASqqVdiflo/fztzAXxfzUXVn3A8A9AxwxJOhXhjexYN3fRERNUMsQkT1qFKqcfjyTexOvoHY1DyUVak0j7VxscUT3b3weDcv+DjaiJiSiIgeFosQ0R0qtYDjGYXYfeYGfj+XC/ntas1jXq2sMTLEEyNDPNDJg+N+iIhaChYhMmpKlRonMovw+7lc7D2fi5sllZrHXOwt8WgXD4wM8USobyuWHyKiFohFiIxOlVKNhCsF2HsuF/tS81BU9r/b3aXW5hjRxR0ju3oioo0TTLnUBRFRi8YiREaholqF+PQC/H5Whj/T8qCoUGoea21jjsGd3DA82AN9gpxhYcY7voiIjAWLELVY5VVKxF28iZizMhy4kF9rwLOznSWGBdeUn4gAR5jxdnciIqPEIkQtys2SShy4kI8/0/JwKP0mKqr/d6u7p9QKQ4PdMTzYA2F+rXnZi4iIWISoeRMEAen5pfgzLQ9/puYh6dotzSSHAODraIPhXWrKT4i3lAOeiYioFhYhanaqVWqcvFqEP1NrzvxkF5XXeryrtxSDOrrhkY6uvNWdiIjui0WImgX57WocvHQTf6bmIe5ifq3BzhZmJugT6IRBndzwSAc3uEutRExKRETNCYsQGSS1WkCqTIGDl27i4MWbSMwuhkr9v2teTrYWGNjBFYM6ueFfQc6wteT/ykREpD1+epDBKC6rQvzlAsRdzMehSwUoKK2s9Xigiy0Gd3LH4E6u6ObDwc5ERPTwWIRINCq1gJScWzVnfS7dxJlrt/C3kz6wsTBFZKAzotq7oH87F67rRUREOsciRE0qp7gcCVcKcTi9APHpN1FcXl3r8Q7u9ohq54Kodi4I828NSzNTkZISEZExYBEivSosrUTClcI7XwXIKqx9h5e9lRn6tnVGVDsX9GvnAg+ptUhJiYjIGLEIkU6VVFTjRGYRjlyuKT4XcktqPW5qIkGItxR9gpzRr50Luvu04qzOREQkGhYheijlVUokZd/C0TtnfM7kyGvd3QXUXO7qE+SMyEAn9AxwhL2VuUhpiYiIamMRIq3Iy6tx8moRTl4twvHMIpy7LofyH8XHz8kGkYFOiAx0Ru9AJzjbWYqUloiI6P5YhOi+8hQVOJFZhBOZNeXnYl5JrSUsAMBDaoWeAY7oE+iMyCAneLfm3V1ERNQ8sAiRRrVKjTSZAqezinE6+xaSrhXjWtHtOvu1cbZFzwBH9PB3RM8AR3i3tuYyFkRE1CyxCBmxPEUFkrLvlJ7sYqTkyFGpVNfaRyIBOnk4oIe/IyICHBHu7wgXe17qIiKiloFFyEjIy6tx7oYcKTlynL1+C2euyXH9Vt2zPVJrc3T3bYVQ39bo7tsKIT6t4MDBzURE1EKxCLVAJRXVOHddgbPXb90pPvI68/cAgIkEaOdmj1C/1pri08bZlpe5iIjIaLAINWOCIOCGvAIXZApcyC1BmkyB1BsKZBSU1bu/r6MNunhL0cVLiq5eUnT1aQU7LlZKRERGjJ+CzUR5lRKX8kqRJlPggkyBtNwSXJApoKhQ1ru/VytrdPWWoou3FF29WiHYywGtbCyaODUREZFhYxEyMEVlVbicX4orN0s1/71ysxQ5xbfr3LYOAGYmEgS52qGjhwM6uNujg4cDunhJ4WjL0kNERPQgLEIiKK1UIruwHNeKy5FdWF6r9PxzEdK/c7G3RAd3e3T0cEBHD3t0cHdAoIsdLMy4RAUREVFjiF6E1q1bhw8++AAymQydO3fG6tWr0bdv33vuf/DgQURHR+P8+fPw9PTE66+/jpkzZzZh4gerVKqQJ69ETnE5sov+93Wt+DauFZWjqKzqvs/3bm2NQBc7BLna/e2/tnDiDM1EREQ6JWoR2rFjB+bOnYt169ahT58++OKLLzB8+HCkpqbC19e3zv6ZmZkYMWIEnn/+eXz33Xc4cuQIZs2aBRcXFzz55JN6zSoIAkorlSguq0ZhWSWKyqpQUFqJXHklchUVyJXfRq6iEnmKigcWHQBobWMOX0cbeDvaINClpugEudqhjbMdrC1M9fqzEBERUQ2JINQ38qRpREREIDQ0FOvXr9ds69ixI0aPHo2VK1fW2f+NN97A7t27kZaWptk2c+ZMnDlzBkePHm3QeyoUCkilUuw4cgGWNnZQqgSUVSlRWqlEacWd//7t18Xl1Sguq0JRWRWqVOoHv8EdlmYm8G5tDR9HG/je+fJuXfNfH0drLjxKRESkhbuf33K5HA4ODjp7XdHOCFVVVSExMRELFiyotX3IkCFISEio9zlHjx7FkCFDam0bOnQoNm3ahOrqapib1y0XlZWVqKys1Hwvl8sBANHfHYOJpfZrYlmZm6C1jQVa21jA0c4cbvZWcHOwgqu9JVylVnCzt4SrvRVa2Zjfcz4eoeo2FFV1JzMkIiKi+ikUCgA1V2h0SbQiVFBQAJVKBTc3t1rb3dzckJubW+9zcnNz691fqVSioKAAHh4edZ6zcuVKLF++vM726+unNj48ERERiaKwsBBSqVRnryf6YOl/njURBOG+MxvXt3992+9auHAhoqOjNd/funULfn5+yM7O1ulvJDWOQqGAj48Prl27ptNTnaQ9HgvDwWNhOHgsDIdcLoevry8cHR11+rqiFSFnZ2eYmprWOfuTn59f56zPXe7u7vXub2ZmBicnp3qfY2lpCUvLundbSaVS/k9tQBwcHHg8DASPheHgsTAcPBaGw8REt1PGiDYBjYWFBcLCwhAbG1tre2xsLCIjI+t9Tu/evevsv2/fPoSHh9c7PoiIiIjofkSdiS86OhobN27E5s2bkZaWhnnz5iE7O1szL9DChQsxefJkzf4zZ85EVlYWoqOjkZaWhs2bN2PTpk2YP3++WD8CERERNWOijhEaP348CgsLsWLFCshkMgQHByMmJgZ+fn4AAJlMhuzsbM3+AQEBiImJwbx587B27Vp4enpizZo1Ws0hZGlpiaVLl9Z7uYyaHo+H4eCxMBw8FoaDx8Jw6OtYiDqPEBEREZGYuEgVERERGS0WISIiIjJaLEJERERktFiEiIiIyGi1yCK0bt06BAQEwMrKCmFhYYiPj7/v/gcPHkRYWBisrKzQpk0bfP75502UtOXT5lj88ssvGDx4MFxcXODg4IDevXvjjz/+aMK0LZ+2fzbuOnLkCMzMzNCtWzf9BjQi2h6LyspKLFq0CH5+frC0tERgYCA2b97cRGlbNm2PxdatWxESEgIbGxt4eHhg2rRpKCwsbKK0LdehQ4cwcuRIeHp6QiKR4Ndff33gc3Ty+S20MNu3bxfMzc2FDRs2CKmpqcKcOXMEW1tbISsrq979MzIyBBsbG2HOnDlCamqqsGHDBsHc3Fz46aefmjh5y6PtsZgzZ47w3nvvCSdOnBAuXbokLFy4UDA3NxdOnz7dxMlbJm2Px123bt0S2rRpIwwZMkQICQlpmrAtXGOOxahRo4SIiAghNjZWyMzMFI4fPy4cOXKkCVO3TNoei/j4eMHExET45JNPhIyMDCE+Pl7o3LmzMHr06CZO3vLExMQIixYtEn7++WcBgLBz58777q+rz+8WV4R69uwpzJw5s9a2Dh06CAsWLKh3/9dff13o0KFDrW0vvvii0KtXL71lNBbaHov6dOrUSVi+fLmuoxmlxh6P8ePHC2+++aawdOlSFiEd0fZY/P7774JUKhUKCwubIp5R0fZYfPDBB0KbNm1qbVuzZo3g7e2tt4zGqCFFSFef3y3q0lhVVRUSExMxZMiQWtuHDBmChISEep9z9OjROvsPHToUp06dQnV1td6ytnSNORb/pFarUVJSovMF9oxRY4/Hli1bcOXKFSxdulTfEY1GY47F7t27ER4ejvfffx9eXl5o164d5s+fj9u3bzdF5BarMcciMjISOTk5iImJgSAIyMvLw08//YRHH320KSLT3+jq81v01ed1qaCgACqVqs6irW5ubnUWa70rNze33v2VSiUKCgrg4eGht7wtWWOOxT999NFHKCsrw7hx4/QR0ag05nikp6djwYIFiI+Ph5lZi/qrQlSNORYZGRk4fPgwrKyssHPnThQUFGDWrFkoKiriOKGH0JhjERkZia1bt2L8+PGoqKiAUqnEqFGj8OmnnzZFZPobXX1+t6gzQndJJJJa3wuCUGfbg/avbztpT9tjcde2bduwbNky7NixA66urvqKZ3QaejxUKhUmTpyI5cuXo127dk0Vz6ho82dDrVZDIpFg69at6NmzJ0aMGIFVq1bhq6++4lkhHdDmWKSmpmL27NlYsmQJEhMTsXfvXmRmZmrWyKSmpYvP7xb1zzxnZ2eYmprWafL5+fl1WuNd7u7u9e5vZmYGJycnvWVt6RpzLO7asWMHZsyYgR9//BGDBg3SZ0yjoe3xKCkpwalTp5CUlISXX34ZQM2HsSAIMDMzw759+zBw4MAmyd7SNObPhoeHB7y8vCCVSjXbOnbsCEEQkJOTg7Zt2+o1c0vVmGOxcuVK9OnTB6+99hoAoGvXrrC1tUXfvn3x9ttv8ypCE9LV53eLOiNkYWGBsLAwxMbG1toeGxuLyMjIep/Tu3fvOvvv27cP4eHhMDc311vWlq4xxwKoORM0depUfP/997zmrkPaHg8HBwecPXsWycnJmq+ZM2eiffv2SE5ORkRERFNFb3Ea82ejT58+uHHjBkpLSzXbLl26BBMTE3h7e+s1b0vWmGNRXl4OE5PaH52mpqYA/nc2gpqGzj6/tRpa3QzcvRVy06ZNQmpqqjB37lzB1tZWuHr1qiAIgrBgwQJh0qRJmv3v3n43b948ITU1Vdi0aRNvn9cRbY/F999/L5iZmQlr164VZDKZ5uvWrVti/QgtirbH459415juaHssSkpKBG9vb+Gpp54Szp8/Lxw8eFBo27at8Nxzz4n1I7QY2h6LLVu2CGZmZsK6deuEK1euCIcPHxbCw8OFnj17ivUjtBglJSVCUlKSkJSUJAAQVq1aJSQlJWmmMtDX53eLK0KCIAhr164V/Pz8BAsLCyE0NFQ4ePCg5rEpU6YIUVFRtfaPi4sTunfvLlhYWAj+/v7C+vXrmzhxy6XNsYiKihIA1PmaMmVK0wdvobT9s/F3LEK6pe2xSEtLEwYNGiRYW1sL3t7eQnR0tFBeXt7EqVsmbY/FmjVrhE6dOgnW1taCh4eH8Mwzzwg5OTlNnLrlOXDgwH0/A/T1+S0RBJ7LIyIiIuPUosYIEREREWmDRYiIiIiMFosQERERGS0WISIiIjJaLEJERERktFiEiIiIyGixCBEREZHRYhEiIr3o378/5s6de999vvrqK7Rq1apJ8hg6iUSCX3/9FQBw9epVSCQSJCcni5qJyBiwCBG1YJ9//jns7e2hVCo120pLS2Fubo6+ffvW2jc+Ph4SiQSXLl3SSxZ/f3+sXr1aL6/dnCxbtgzdunWrs10mk2H48OFNH4jIyLEIEbVgAwYMQGlpKU6dOqXZFh8fD3d3d5w8eRLl5eWa7XFxcfD09ES7du3EiCq6qqoqUd/f3d0dlpaWomYgMkYsQkQtWPv27eHp6Ym4uDjNtri4ODz++OMIDAxEQkJCre0DBgwAUFMKXn/9dXh5ecHW1hYRERG1XqOwsBATJkyAt7c3bGxs0KVLF2zbtu2eOfr374+srCzMmzcPEokEEomk1uN//PEHOnbsCDs7OwwbNgwymeyerxUXFweJRII9e/YgJCQEVlZWiIiIwNmzZ2vtl5CQgH79+sHa2ho+Pj6YPXs2ysrKNI/7+/vj7bffxtSpUyGVSvH8888DAI4cOYKoqCjY2NigdevWGDp0KIqLiwHUrC7+/vvvo02bNrC2tkZISAh++umnOtn279+P8PBw2NjYIDIyEhcvXgRQcylw+fLlOHPmjOb34auvvgJQ+9JYfVJTUzFixAjY2dnBzc0NkyZNQkFBwT33J6KGYREiauH69++PAwcOaL4/cOAA+vfvj6ioKM32qqoqHD16VFOEpk2bhiNHjmD79u1ISUnB2LFjMWzYMKSnpwMAKioqEBYWhv/+9784d+4cXnjhBUyaNAnHjx+vN8Mvv/wCb29vrFixAjKZrFbRKS8vx4cffohvv/0Whw4dQnZ2NubPn//An+u1117Dhx9+iJMnT8LV1RWjRo1CdXU1AODs2bMYOnQonnjiCaSkpGDHjh04fPgwXn755Vqv8cEHHyA4OBiJiYlYvHgxkpOT8cgjj6Bz5844evQoDh8+jJEjR0KlUgEA3nzzTWzZsgXr16/H+fPnMW/ePDz77LM4ePBgrdddtGgRPvroI5w6dQpmZmaYPn06AGD8+PF49dVX0blzZ83vw/jx4x/4s8pkMkRFRaFbt244deoU9u7di7y8PIwbN+6BzyWiB3jY1WKJyLB9+eWXgq2trVBdXS0oFArBzMxMyMvLE7Zv3y5ERkYKgiAIBw8eFAAIV65cES5fvixIJBLh+vXrtV7nkUceERYuXHjP9xkxYoTw6quvar6PiooS5syZo/nez89P+Pjjj2s9Z8uWLQIA4fLly5pta9euFdzc3O75PndXqN6+fbtmW2FhoWBtbS3s2LFDEARBmDRpkvDCCy/Uel58fLxgYmIi3L59W5Nn9OjRtfaZMGGC0KdPn3rft7S0VLCyshISEhJqbZ8xY4YwYcKEWtn+/PNPzeN79uwRAGjed+nSpUJISEid1wcg7Ny5UxAEQcjMzBQACElJSYIgCMLixYuFIUOG1Nr/2rVrAgDh4sWL9eYlooYxE6+CEVFTGDBgAMrKynDy5EkUFxejXbt2cHV1RVRUFCZNmoSysjLExcXB19cXbdq0wY8//ghBEOqMFaqsrISTkxMAQKVS4d1338WOHTtw/fp1VFZWorKyEra2tlrns7GxQWBgoOZ7Dw8P5OfnP/B5vXv31vza0dER7du3R1paGgAgMTERly9fxtatWzX7CIIAtVqNzMxMdOzYEQAQHh5e6zWTk5MxduzYet8vNTUVFRUVGDx4cK3tVVVV6N69e61tXbt2rfXzAEB+fj58fX0f+HPVJzExEQcOHICdnV2dx65cuWK047qIdIFFiKiFCwoKgre3Nw4cOIDi4mJERUUBqBmcGxAQgCNHjuDAgQMYOHAgAECtVsPU1BSJiYkwNTWt9Vp3P4g/+ugjfPzxx1i9ejW6dOkCW1tbzJ07t1EDjs3NzWt9L5FIIAhCY35UzdgjtVqNF198EbNnz66zz9/LyD+Lm7W19T1fW61WAwD27NkDLy+vWo/9c5Dz33+mv2dqLLVajZEjR+K9996r89jdokVEjcMiRGQEBgwYgLi4OBQXF+O1117TbI+KisIff/yBY8eOYdq0aQCA7t27Q6VSIT8/v84t9nfFx8fj8ccfx7PPPgug5oM6PT1dc6alPhYWFpqxNrpw7NgxTakpLi7GpUuX0KFDBwBAaGgozp8/j6CgIK1es2vXrti/fz+WL19e57FOnTrB0tIS2dnZmjLZGI35fQgNDcXPP/8Mf39/mJnxr20iXeJgaSIjMGDAABw+fBjJycm1PsSjoqKwYcMGVFRUaAZKt2vXDs888wwmT56MX375BZmZmTh58iTee+89xMTEAKg5yxQbG4uEhASkpaXhxRdfRG5u7n0z+Pv749ChQ7h+/bpO7nZasWIF9u/fj3PnzmHq1KlwdnbG6NGjAQBvvPEGjh49ipdeegnJyclIT0/H7t278corr9z3NRcuXIiTJ09i1qxZSElJwYULF7B+/XoUFBTA3t4e8+fPx7x58/D111/jypUrSEpKwtq1a/H11183OLe/vz8yMzORnJyMgoICVFZWPvA5L730EoqKijBhwgScOHECGRkZ2LdvH6ZPn67TcklkjFiEiIzAgAEDcPv2bQQFBcHNzU2zPSoqCiUlJQgMDISPj49m+5YtWzB58mS8+uqraN++PUaNGoXjx49r9lm8eDFCQ0MxdOhQ9O/fH+7u7poSci8rVqzA1atXERgYCBcXl4f+md59913MmTMHYWFhkMlk2L17NywsLADUnNk5ePAg0tPT0bdvX3Tv3h2LFy9+4GWkdu3aYd++fThz5gx69uyJ3r17Y9euXZqzMG+99RaWLFmClStXomPHjhg6dCh+++03BAQENDj3k08+iWHDhmHAgAFwcXG577QDd3l6euLIkSNQqVQYOnQogoODMWfOHEilUpiY8K9xoochERp7MZ6ISAR35zsqLi7m8hxE9ND4TwkiIiIyWixCREREZLR4aYyIiIiMFs8IERERkdFiESIiIiKjxSJERERERotFiIiIiIwWixAREREZLRYhIiIiMlosQkRERGS0WISIiIjIaLEIERERkdH6f2EFA11QeE+kAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The aggregate wealth to income ratio is 1.9763887031145013\n",
      "The Gini Coefficient for assests is 1.0\n",
      "The average MPC by income quintile is [0.10692187291548673, 0.09327176601929044, 0.08810286571670196, 0.07414283841778725, 0.046026708398195525]\n"
     ]
    }
   ],
   "source": [
    "# Call the function with test values for (CorrAct, CorrPcvd, DiscFac_center, DiscFac_spread)\n",
    "AggWealthRatio, Lorenz, Gini, Avg_MPC = runRoszypalSchlaffmanExperiment(\n",
    "    0.97, 0.9831, 0.9867, 0.0067\n",
    ")\n",
    "\n",
    "# Plot the Lorenz curve\n",
    "print(\"The Lorenz curve for assests is\")\n",
    "plt.plot(Lorenz[0], Lorenz[1])\n",
    "plt.xlabel(\"Wealth percentile\")\n",
    "plt.ylabel(\"Cumulative wealth share\")\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.0])\n",
    "plt.show()\n",
    "\n",
    "print(\"The aggregate wealth to income ratio is \" + str(AggWealthRatio))\n",
    "print(\"The Gini Coefficient for assests is \" + str(Gini))\n",
    "print(\"The average MPC by income quintile is \" + str(Avg_MPC))"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "collapsed,code_folding",
   "cell_metadata_json": true,
   "formats": "ipynb,py:percent",
   "notebook_metadata_filter": "all"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}