cstwMPC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Distribution of Wealth and the MPC\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook presents a selection of results from the paper [The Distribution of Wealth and the Marginal Propensity to Consume](http://econ.jhu.edu/people/ccarroll/papers/cstwMPC), using the [Econ-ARK/HARK](https://github.com/econ-ark/HARK) toolkit.  It sketches the steps that would need to be taken to solve the model using the [dolARK](https://github.com/EconForge/dolARK) modeling system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# This cell does some standard python setup\n",
    "\n",
    "# Tools for navigating the filesystem\n",
    "import sys\n",
    "import os\n",
    "\n",
    "# Import related generic python packages\n",
    "import numpy as np\n",
    "from time import clock\n",
    "mystr = lambda number : \"{:.4f}\".format(number)\n",
    "from copy import copy, deepcopy\n",
    "\n",
    "# Plotting tools\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import plot, draw, show\n",
    "\n",
    "# iPython gives us some interactive and graphical tools\n",
    "from IPython import get_ipython # In case it was run from python instead of ipython\n",
    "\n",
    "# The warnings package allows us to ignore some harmless but alarming warning messages\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Abstract\n",
    "\n",
    "In a model calibrated to match micro- and macroeconomic evidence on household income dynamics, this paper shows that a modest degree of heterogeneity in household preferences or beliefs is sufficient to match empirical measures of wealth inequality in the U.S. The hegerogeneity-augmented model's predictions are consistent with microeconomic evidence that suggests that the annual marginal propensity to consume (MPC) is much larger than the roughly 0.04 implied by commonly-used macroeconomic models (even ones including some heterogeneity).  The high MPC arises because many consumers hold little wealth despite having a strong precautionary motive.  The model also plausibly predicts that the aggregate MPC can differ greatly depending on how the shock is distributed across households (depending, e.g., on their wealth, or employment status)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notation For the Core Model\n",
    "\n",
    "We define the following notation.\n",
    "\n",
    "| Exogenous Random Variable | Description | Code | Value |\n",
    "| :---: | ---         | ---  | :---: |\n",
    "| $\\newcommand{\\tshk}{\\zeta}\\tshk$ | Transitory Income | $\\texttt{tshk}$ |  |\n",
    "| $\\newcommand{\\pshk}{\\psi}\\pshk$ | Permanent Shock | $\\texttt{pshk}$ |  |\n",
    "\n",
    "| Parameter | Description | Code | Value |\n",
    "| :---: | ---         | ---  | :---: |\n",
    "| $\\newcommand{\\PLives}{\\Lambda} \\PLives$ | Probability of living | $\\texttt{PLives}$ | 0.99375\n",
    "| $\\newcommand{\\Discount}{\\beta}\\Discount$ | Time Preference Factor | $\\texttt{Discount}$ | 0.96 |\n",
    "| $\\newcommand{\\CRRA}{\\rho}\\CRRA$ | Coeﬃcient of Relative Risk Aversion| $\\texttt{CRRA}$ | 1 |\n",
    "| $\\sigma_\\tshk$ | Transitory Income Standard Deviation | $\\texttt{PermShkStd}$ | 0.1  |\n",
    "| $\\sigma_\\pshk$ | Permanent Shock Standard Deviation |$\\texttt{TranShkStd}$ | 0.1 |\n",
    "\n",
    "\n",
    "| Variable | Description | Code | \n",
    "| :---: | ---         | ---  | \n",
    "| $\\newcommand{\\aRat}{a}\\aRat$ | Assets | $\\texttt{aRat}$ |  \n",
    "| $\\newcommand{\\mRat}{m}\\mRat$ | Market resources | $\\texttt{mRt}$ |  \n",
    "| $\\newcommand{\\KLev}{K}\\KLev$ | Capital Aggregate | $\\texttt{KLev}$\n",
    "| $\\newcommand{\\kapShare}{\\alpha}\\kapShare$ | Capital share | $\\texttt{kapShare}$\n",
    "| $\\newcommand{\\LLev}{L}\\LLev$ | Labor Aggregate | $\\texttt{LLev}$\n",
    "| $\\newcommand{\\labor}{\\ell}\\labor$ | Labor share | $\\texttt{labor}$ \n",
    "| $\\newcommand{\\kRat}{k}\\kRat$ | $K/P$| $\\texttt{kRat}$ |  \n",
    "| $\\newcommand{\\pRat}{p}\\pRat$ | Permanent Income | |  \n",
    "| $\\mathbf{P}$ | | | \n",
    "| $\\newcommand{\\rProd}{r}\\rProd$ | Interest rate | $\\texttt{rProd}$ |  \n",
    "| $\\newcommand{\\yLev}{y} \\yLev$ | Income | $\\texttt{yLev}$ \n",
    "| $\\newcommand{\\Wage}{W}\\Wage$ | Aggregate Wage Rate | $\\texttt{Wage}$ | \n",
    "\n",
    "| Functions | Description | Code | Value |\n",
    "| :---: | ---         | ---  | :---: |\n",
    "| $\\newcommand{\\cFunc}{\\mathrm{c}}\\cFunc$ | Consumption | $\\texttt{cFunc}$ |  |\n",
    "| $\\newcommand{\\valfn}{\\mathrm{v}} \\valfn$ | Value | $\\texttt{valfn}$\n",
    "| $\\newcommand{\\uFunc}{{\\mathrm{u}}}\\uFunc$ | Utility | $\\texttt{uFunc}$ |  |\n",
    "\n",
    "$\\newcommand{\\cRat}{c}$\n",
    "$\\newcommand{\\Ex}{\\mathbb{E}}$\n",
    "$\\newcommand{\\PDies}{\\mathsf{P}}$\n",
    "$\\newcommand{\\ptyLev}{a}$\n",
    "$\\newcommand{\\YLev}{Y}$\n",
    "$\\newcommand{\\wEndRat}{\\aRat}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The consumer has a standard Constant Relative Risk Aversion utility function $$u(c)=\\frac{c^{1-\\rho}}{1-\\rho}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idiosyncratic (household) income process is logarithmic Friedman:\n",
    "\\begin{eqnarray*}\n",
    "\\yLev_{t+1}&=&\\pRat_{t+1}\\tshk_{t+1}\\Wage\\\\\n",
    "\\pRat_{t+1}&=&\\pRat_{t}\\pshk_{t+1}\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bellman form of the value function for households is:\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\valfn(\\mRat_{t})&=&\\underset{\\cFunc_{t}}{\\max } ~~ \\uFunc(\\cFunc_{t}(\\mRat_t))+\\Discount \\PLives \\Ex_{t}\\left[ \\pshk_{t+1}^{1-\\CRRA}\\valfn(\\mRat_{t+1})\n",
    "\\right]   \\\\\n",
    "\\notag &\\text{s.t.}&\\\\\n",
    "\\wEndRat_{t} &=&\\mRat_{t}-\\cRat_{t},\\\\\n",
    "\\wEndRat_{t} &\\geq &0, \\\\\n",
    "\\kRat_{t+1} &=&\\wEndRat_{t}/(\\PLives \\pshk_{t+1}),\n",
    "\\\\\n",
    "\\mRat_{t+1} &=&(\\daleth +\\rProd_t)\\kRat_{t+1}+\\tshk_{t+1},\\\\\n",
    "\\rProd &=&\\kapShare\\ptyLev(\\KLev/\\labor\\LLev)^{\\kapShare-1}\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "This will run the absolute minimum amount of work that actually produces\n",
    "relevant output-- no aggregate shocks, perpetual youth, matching net worth.\n",
    "Will run both beta-point and beta-dist versions.\n",
    "'''\n",
    "import os\n",
    "\n",
    "'''\n",
    "Copied here from do_min.py.\n",
    "Design decisions about whether to include this code explicitly,\n",
    "or import it, or execute it as is here, TBD.\n",
    "'''\n",
    "\n",
    "here = os.path.dirname(os.path.realpath(\"cstwMPC.ipynb\"))\n",
    "my_path = os.path.join(here,'')\n",
    "path_to_models = os.path.join(my_path,'Code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# For speed here, use the \"tractable\" version of the model\n",
    "# This is not the \"right\" model, but illustrates the key point\n",
    "'''\n",
    "This options file specifies parameter heterogeneity, making the choice in the paper:\n",
    "uniformly distributed discount factors.\n",
    "'''\n",
    "param_name = 'DiscFac'        # Which parameter to introduce heterogeneity in\n",
    "dist_type = 'uniform'         # Which type of distribution to use\n",
    "\n",
    "'''\n",
    "This options file specifies the \"standard\" work options for cstwMPC, estimating the model only.\n",
    "'''\n",
    "run_estimation = True         # Runs the estimation if True\n",
    "run_sensitivity = [False, False, False, False, False, False, False, False] # Choose which sensitivity analyses to run: rho, xi_sigma, psi_sigma, mu, urate, mortality, g, R\n",
    "find_beta_vs_KY = False       # Computes K/Y ratio for a wide range of beta; should have do_beta_dist = False\n",
    "do_tractable = True          # Uses a \"tractable consumer\" rather than solving full model when True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beginning an estimation with the specification name PYbetaPointIndNetWorth...\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.95 and got -9.737467406398238\n",
      "The given type violates the Absolute Impatience Condition with the supplied parameter values; the AIF is 1.00495  \n",
      "The given parameter values violate the Growth Impatience Condition for this consumer type; the GIF is: 1.0050 \n",
      "\n",
      "\n",
      "The given parameter values violate the Growth Impatience Condition for this consumer type; the GIF is: 1.0050 \n",
      "getKYratioDifference tried center = 0.995 and got 65.56174461486972\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.9558192645251606 and got -9.682962081909512\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.9754096322625803 and got -9.30238283304393\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.9852048161312901 and got -8.202435342294791\n",
      "The given type violates the Absolute Impatience Condition with the supplied parameter values; the AIF is 1.00000  \n",
      "The given parameter values violate the Growth Impatience Condition for this consumer type; the GIF is: 1.0000 \n",
      "\n",
      "\n",
      "The given parameter values violate the Growth Impatience Condition for this consumer type; the GIF is: 1.0000 \n",
      "getKYratioDifference tried center = 0.9901024080656451 and got 5.330284130112025\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.9881733385525523 and got -3.9914530499444387\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.9889993423588024 and got -0.667776237950152\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.989146822569607 and got 0.034770717735451484\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.9891395234236353 and got -0.000711323657323959\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "getKYratioDifference tried center = 0.9891400234236357 and got 0.0017168379920207144\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "Estimate is center=0.9891395234236353, spread=0.0, took 31.892662 seconds.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAENCAYAAADgwHn9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZzNZfvA8c81g8HYt599J2tJkpTIli2UJKVHCm1IFOpp97TJk6VFSUg9lUiMZEmLZDw0JFmSyTbD2A0zxixmrt8f58w8h2Y4xpzzPTNzvV+v85rvfl/HGeea+76/3/sWVcUYY4wBCHI6AGOMMYHDkoIxxpgMlhSMMcZksKRgjDEmgyUFY4wxGSwpGGOMyeCXpCAiM0XksIhsyWK/iMhUEYkUkc0i0twfcRljjDmXv2oKs4EuF9jfFajnfg0FpvkhJmOMMefxS1JQ1Z+A4xc4pBcwR13+C5QSkUr+iM0YY8z/FHA6ALcqQJTHerR7W8z5B4rIUFy1CUJDQ69p0KCBXwI0xphAcPr0af744w8qV65MpUrZ+9t5w4YNR1W1fGb7AiUpSCbbMh1/Q1WnA9MBWrRooREREb6MyxhjAsrEiRN58sknGTZsGE899VS2riEie7PaFyh3H0UD1TzWqwIHHIrFGGMCVvofwvXr1/fJ9QMlKYQB/3DfhdQKOKmqf2s6MsaY/G7z5s0AlCxZ0ifX90vzkYh8BrQDyolINPA8UBBAVd8DvgG6AZFAAjDIH3EZY0xukpSUxM6dOwEIDQ31SRl+SQqq2v8i+xV41B+xGGNMbrV3716KFClCXFwcxYoV80kZgdJ8ZIwx5iLq16/PlClTACwpGGOMcd2SCrm8+cgYY8zl6969O8HBwYDVFIwxJl87ceIE33zzDfHx8YgIRYoU8Uk5lhSMMSYXWL9+PQBlypQhNDQUkcye+b18lhSMMSYXWLduHSJCsWLFfNZ0BJYUjDEmV1i/fj0NGzYkJSXFZ53MYEnBGGNyhfr169OnTx/i4+N9WlOwu4+MMSYXePPNNwHo0KGD1RSMMSY/O3DgAGlpaYDrOQXrUzDGmHzslltuoW/fvgDEx8dbTcEYY/KrrVu3smXLFm6++WbAagrGGJOvzZ07l6CgIO644w4An3c0W1IwxpgAparMnTuXtm3bUrFiRcBVU7DmI2OMyYc2bdrEn3/+yV133QVAamoqZ86csVtSjTEmP2rcuDGLFy+mdevWgO9HSAVLCsYYE7AKFSpEjx49MtbTk4L1KRhjTD7zyy+/8Oyzz3Ls2LGMbfHx8YAlBWOMyXcmT57M1KlTKVy4cMY2fzQfWVIwxpgAc+TIEebPn8/AgQPPSQBWUzDGmHxo5syZJCcn89BDD52z3WoKxhiTz6SmpvLee+/Rrl07GjVqdM4+f9QU7O4jY4wJICdPnuTaa6/NeDbBkyUFY4zJZ8qUKcMXX3yR6T5rPjLGmHwkJiaG7du3Z7nfOpqNMSYfeeWVV2jWrBlHjx7NdH96TaFIkSI+i8GSgjHGBIAjR47w4Ycfcs8991CuXLlMj0mfSyEoyHdf3ZYUjDEmALz11lskJiby5JNPZnmMryfYAUsKxhjjuPj4eN5++2169epFw4YNszzO1xPswCUkBREpKCJtRKSfez1URHybsowxJh/45ZdfSExMZOzYsRc8ztcT7ICXt6SKSFMgDEgCqgJzgbbAQKCfz6Izxph84Oabb+bAgQOUKlXqgsf5eoId8L6mMA14TlUbACnubauAG30SlTHG5BMxMTGo6kUTAvinpuBtUmgMfOJeVgBVPQ347r4oY4zJ45KSkrjuuusYNmyYV8cHUk1hD3CN5wYRaQlE5nRAxhiTX3zwwQdERUXRu3dvr44PpJrCs8ASEXkRKCQiTwHzgGe8LUhEuojIDhGJFJFxmeyvLiI/iMivIrJZRLp5e21jjMltzpw5wyuvvEKbNm3o2LGjV+cETEezqn4tIl2Bwbj6EmoAt6vqBm/OF5Fg4B2gExAN/CIiYaq6zeOwZ4AvVHWaiDQCvgFqev1OjDEmF5k2bRoxMTF8/vnniIhX5/ij+eiiScH9hT4TGKqqj2SznJZApKrucl/zc6AX4JkUFCjhXi4JHMhmWcYYE9BUlY8//piOHTty0003eXVOWloaCQkJztcUVDVVRDoDaZdRThUgymM9GrjuvGNeAFaIyHAgFMi0PiUiQ4GhANWrV7+MkIwxxhkiQnh4+DnzL19MQkIC4NsRUsH7PoVJwIsiUjCb5WRWN9Lz1vsDs1W1KtAN+FhE/hafqk5X1Raq2qJ8+fLZDMcYY5xx8uRJkpKSKFKkCFWrVvX6PH+MkAreJ4XhwJNAnIhEici+9JeX50cD1TzWq/L35qEHgC8AVHUtUBjIfFQoY4zJpcaOHctVV11FcnLyJZ2XPkKq481HbgMus5xfgHoiUgvYD9wF3H3eMfuADsBsEWmIKykcucxyjTEmYGzbto0PPviARx99lEKFCl3Suek1Bcc7mgFUddXlFKKqZ0VkGLAcCAZmqupWEXkJiFDVMGA08IGIPI6raek+VT2/ickYY3KtsWPHUqxYMZ577rlLPtdfzUdeT8cpIs2ANriadDL6CFTVq3enqt/gus3Uc9tzHsvbgBu8jccYY3KT77//nq+//prXX389y/kSLsQfU3GCl30K7jt+1gDtgbFAU1x/2df1XWjGGJN3zJkzh5o1azJixIhsnR9oNYUxQBdVXS0iJ1T1NvfDbHf5MDZjjMkzZs6cSXR0NIULF87W+QFVUwAqqOpq93KaiASp6lLgVh/FZYwxecLRo0c5fPgwQUFBl/VsVaDdkhotIjXdy38CvUSkDXBp91QZY0w+M2bMGJo2bZrxl352+euWVG+TwgQgfY64l3ANo/098KIvgjLGmLxg3bp1zJo1i/vuu++ym33SawpFixbNidCy5O0tqbM9lpeKSGmgkKrG+yowY4zJzVJSUhg6dChVqlThmWe8HlA6S/Hx8RQtWpSgIK9nUc4Wr29JBRCRCkAxz/X0Qe6MMcb8z+TJk9m8eTMLFiygePHil309f4yQCt7P0dwF+BCoyLnjGCmuh9GMMcZ42LFjB7169eK2227Lkev5Yy4F8L6m8A4wHvhIVc/4MB5jjMkTZsyYccnjG13I6dOn/ZIUvG2cKg28bwnBGGMubOXKlWzduhXgksc3upD4+Hi/NB95mxQ+BAb5MhBjjMntjhw5Qv/+/Rk6dCg5PXSbv2oKWTYfichq/jfngQCPuedWPuh5nKp6N22QMcbkccOGDePUqVO8//77Xk+x6a34+Hj8MYfMhfoUZlxk3RhjjNuCBQv44osv+Ne//kWTJk1y/PqOdzSr6kc+L90YY/KAY8eO8cgjj3D11VczZswYn5QRUB3NItLfPfENInKFiKwSke9FpIFvwzPGmMAXGhrKoEGDmDVrFgULZnfW4gvzV0ezt7ek/gto7V6eiGsmtXjgXVzDaRtjTL6kqhQuXJhXX33VZ2WkpaWRkJAQODUFoLyqHhKRwsCNwD9xjYHUzGeRGWNMgNu/fz/XXHMN69ev92k5Z86cQVUDqqZwRETq4ppc5xdVTRKRopz7dLMxxuQbaWlpDBo0iB07dlCmTBmfluWvYbPB+6QwHtgApAL93Ns6AL/5IihjjAl077zzDt9++y3Tpk2jbl3fTkJ56tQpIICSgqrOFpEv3MsJ7s3rsJnXjDH50ObNm3nyySfp3r07Dz74oM/L27NnDwA1atTweVlej5LqkQzS1w/nfDjGGBP43nvvPUqXLs2sWbNy/CG1zOzcuROAevXq+bws3w7MbYwxedDbb7/NmjVr/PKEMUBkZCRFixalUqVKPi/LkoIxxnhp1apV7N+/n6CgIGrXru23cnfu3EndunX9UiuxpGCMMV7YvXs3vXr14oEHHvB72ZGRkT7vzE7ndZ+CiJQErsBj5jUAVf0+p4MyxphAkpSUxJ133gnAtGnT/Fp2amoqu3btolevXn4pz9uZ1+7DNdFOPODZ4ayA/+pQxhjjgCeffJKIiAgWLFhArVq1/Fp2VFQUycnJAVdTeBm4Q1WX+jIYY4wJNIsXL+att95i5MiROTa15qXw551H4H2fQgFghS8DMcaYQNS2bVuefvppXn/9dUfKj4yMBPBbTcHbpPA68IyIWMe0MSZfiI+P58yZM5QoUYKXX345R6fWvBQ7d+6kSJEiVK5c2S/lXWjmtSjOnXmtIjBGRI55Hqeq1X0XnjHG+J+qcv/997Nnzx7WrFnjs+GwvZF+55E/bkeFC/cpDPBLBMYYE2DefPNN5s2bx+uvv+5oQgBXTaFRo0Z+K+9CM6+tSl8Wkb6qOu/8Y0TkDl8FZowxTvj2228ZM2YMffr04cknn3Q0lvTbUXv27Om3Mr3tI/gwi+3TcyoQY4xx2l9//UW/fv1o1KgRs2fP9luTTVbSb0f1151HcJGkICK1RaQ2ECQitdLX3a+OQKK3BYlIFxHZISKRIjIui2PuFJFtIrJVRD69tLdijDGXR0Ro2rQpixYt8ssw1ReTfjuqv+48gos/pxCJq7NZgL/O23cQeMGbQkQkGNfDb52AaOAXEQlT1W0ex9QDngJuUNUTIlLBq3dgjDGXKS0tDRGhdu3arFq16uIn+En67agBU1NQ1SBVDQZWu5c9X5VV1dvmo5ZApKruUtVk4HPg/Ge2hwDvqOoJd9k2NLcxxi+efvpp+vfvT0pKitOhnCP9dlR/jI6azqs+BVVte5nlVAGiPNaj3ds81Qfqi8gaEfmviHS5zDKNMeaiZs+ezeuvv06pUqUoUMDr4eD8Iv121KAg/z0idqHnFFbzv+cUsqSqN3lRTma9NedfuwBQD2gHVAVWi0gTVY09L66hwFCA6tXtEQljTPatXr2aoUOH0qFDB9566y3HO5bPt3PnTho2bOjXMi+UFmfkYDnRQDWP9arAgUyO+a+qpgC7RWQHriTxi+dB7iar6QAtWrS4aNIyxpjMREZGctttt1GrVi3mzZvn+PMI53PidlS48HMKH+VgOb8A9USkFrAf19zOd593zEKgPzBbRMrhak7alYMxGGNMhkOHDlGmTBmWLFlC6dKlnQ7nb/w9Omq6S5lP4f9wdRiXw6M5SFVnXuxcVT0rIsOA5UAwMFNVt4rIS0CEqoa593UWkW1AKvCkqh7L+qrGGHPpVBUR4YYbbmDbtm0B14+QbuPGjQA0bdrUr+V6O59Cb+ATYCfQGNgKNAF+Bi6aFABU9Rvgm/O2PeexrMAo98sYY3JcWload999N82aNWPcuHEBmxAAwsPDCQkJoXnz5n4t19su7X8Bg1T1auC0++dQYIPPIjPGmBykqowePZq5c+f69W6e7AoPD6dFixZ+H53V23+Z6pmMffQR8I8cjscYY3zizTffZPLkyYwYMcLxMY0uJjExkQ0bNtC6dWu/l+1tUjjs7lMA2CMi1wN1cPUPGGNMQPv000954okn6Nu3L5MmTQq4W0/Pt3HjRpKTkwM6KXwA3OhengT8APwGvOuLoIwxJiclJCTQvn175syZk2uajgCuv/56v5ctrv7dSzxJpDoQqqrbcz4k77Vo0UIjIiKcDMEYE8ASExMpXLgw4Opkzg0JAaBPnz789ttvGWMf5TQR2aCqLTLb5/W/kIgUFJE2ItJPVfcB+0QkNMeiNMaYHLR161bq1q3LypUrAXJNQlBVwsPDHWk6Ai+Tgog0Bf7E1YyUPrdCW7y8HdUYY/xp165ddO7cGVWlTp06TodzSfbs2cPBgwcDOykA04DnVLUBkD6M4Cr+189gjDEBYf/+/XTo0IHExERWrFhBrVq1nA7pkjjZnwDeP9HcGNfDa+AeyE5VT4tIEZ9EZYwx2XDixAk6duzIsWPH+P7772ncuLHTIV2y8PBwihUrRpMmTRwp39uawh7gGs8NItIS1yQ8xhgTEEqUKEGnTp1YsmQJLVpk2o8a8MLDw2nVqhXBwc7c8e9tUngWWCIiLwKFROQpYB7wjM8iM8YYL8XGxhIVFUVwcDBTp06lTZs2ToeULXFxcWzevNmx/gTwfpKdr4GuQHlcfQk1gNtVdYUPYzPGmIuKi4uja9eudOjQIeBmTrtUP//8M2lpadxwww2OxeD1aFCquhF4xIexGGPMJYmPj6dbt25EREQE5JwIlyosLIzQ0FBuusmbuct8w9tRUo8BP+GqJawCNml2nnozxpgckp4Q1q5dy2effUbv3r2dDumypKWlERYWRpcuXTIeuHOCt30K1wFhQDNgAXBcRBaLyBM+i8wYYy5g3LhxhIeH8+mnn9K3b1+nw7lsGzZs4MCBA/Tq1cvROLyqKahqJK47jWaJSH1co6MOB24BJvouPGOMydz48ePp0aMHXbp0cTqUHLFo0SKCg4Pp3r27o3F4+0TzQyLymYjsA+YAhXBNnVnOl8EZY4ynkydP8uSTT5KYmEjp0qXzTEIAV1Jo06YNZcqUcTQObzua3wX+AsYDX6tqjO9CMsaYvzt+/Di33HILv/32Gz179sy1t51mZteuXWzZsoVJkyY5HYrXSaEqrrGObgIeE5GCuDueVfWTC55pjDGX6ciRI3Tu3Jnt27fz1Vdf5amEAK5aAkDPnj0djiQbQ2eLyNXAHcAwoJiqOjbRjg2dbUzet3//fjp27MjevXtZuHAhnTt3djqkHNeuXTuOHTvG77//7pfyLnvobBF5XEQWichxYBFQHXgCaJhzYRpjzN/FxsaSlJTE8uXL82RCOHbsGKtXr3b8rqN03jYfXQUsBB5X1V0+jMcYYwCIiYmhYsWKNG7cmB07duT6B9OyMn/+fNLS0rjtttucDgXwfpiL+1R1liUEY4w/rF+/niZNmvDmm28C5NmEADBjxgyuvPJKmjdv7nQowCXMvGaMMf6wcuVK2rdvT6lSpQLmr2df2bRpExEREQwePBgRcTocwJKCMSaAzJs3j+7du1O7dm1+/vlnateu7XRIPvXhhx8SEhLCPffc43QoGSwpGGMCwr59+7jnnnu49tprWbVqFZUqVXI6JJ86c+YMn3zyCX369HH8gTVPXo+SCiAi1YAqqvpfH8VjjMmnqlevztdff02bNm0oUiTvT+q4YMECYmNjGTx4sNOhnMPbW1Kri8ga4A9gpXvbHSIyw5fBGWPytpSUFIYMGUJYWBgAnTt3zhcJAVwdzHXq1KFt27ZOh3IOb5uP3geWAMWB9FksvgU6+SIoY0zeFxcXx6233sqMGTPYsmWL0+H41c6dO/nxxx+5//77CQoKrFZ8b5uPWgLdVTVNRBRAVU+KSEnfhWaMyasOHDhAjx492Lx5MzNmzOCBBx5wOiS/mjJlCgULFuS+++5zOpS/8TYpHALqAn+mbxCRRsA+XwRljMm7jhw5QqtWrTh+/DhhYWF069bN6ZD86siRI8ycOZN7772XypUrOx3O33hbb5kIfC0ig4ACItIfmAu87rPIjDF5Urly5bjvvvtYvXp1vksIAG+//TZnzpzhiScCc44yrwfEE5HewFCgBq4awvuqutCHsV2UDYhnTO4xY8YMWrVqRZMmTZwOxTGnT5+mevXqtGnThoULnfv6zIkB8YJVdaGqdlPVxqra1emEYIzJHVJTU3n88ccZMmQIU6dOdTocR3344YccP36cMWPGOB1KlrxtPjooIu+KyA3ZLUhEuojIDhGJFJFxFzjuDhFREck0ixljco+4uDh69erF5MmTGTFiBO+++67TITkmJSWFf//739x44420bt3a6XCy5G1S6AzEA5+JyB4ReVVEmnpbiIgEA+8AXYFGQH93R/X5xxUHRgDrvL22MSYwxcTE0Lp1a5YtW8Y777zDlClTKFDgkp6XzVM+++wz9u3bF9C1BPB+lNRfVXWMqlYHBgKlge9EZLOX5bQEIlV1l6omA58DmQ0ePh6YACR6eV1jTIAqW7YsderUYfny5TzyyCNOh+Oo5ORknn/+eZo3b0737t2dDueCsvPUxA5gOxAF1PTynCru49NFu7dlcM/oVk1Vv77QhURkqIhEiEjEkSNHvA7aGON7qsrs2bM5evQohQoVYuHChXTo0MHpsBw3ffp09uzZwyuvvBJwD6udz9uO5lIi8oCIfAf8BbTDdTtqBS/LyWxM2IzbnkQkCJgEjL7YhVR1uqq2UNUW5cuX97J4Y4yvJScn8+CDDzJo0KB836HsKT4+nvHjx9OuXbtcMXOctw18B4Bw4FPgdlU9eYnlRAPVPNaruq+ZrjjQBPjRPaZ4RSBMRHqqqt1zakyAO3jwIH369CE8PJxx48bx/PPPOx1SwJgyZQqHDx9m0aJFATNnwoV4mxTqqGrMZZTzC1BPRGoB+4G7gLvTd7qTTLn0dRH5EXjCEoIxgW/Tpk1069aN2NhYPv/8c/r16+d0SAHj2LFjTJgwgZ49e9KqVSunw/FKlklBRG5S1Z/cqw1FpGFmx6nq9xcrRFXPisgwYDkQDMxU1a0i8hIQoaph2YjdGBMAqlSpQoMGDZg8eTJXXnml0+EElPHjxxMXF8e//vUvp0PxWpZPNIvIFlVt4l7encX5qqqOTY1kTzQb44yEhASmTJnCE088kafnT74c27Zt48orr2Tw4MG89957Todzjgs90ZxlTSE9IbiXa/kiMGNM7hMZGUmfPn34/fffad68ObfccovTIQUcVeWxxx6jePHiuaqWAN7ffbQoi+0LcjYcY0wgW7hwIS1atCA6OpolS5ZYQsjCokWLWLlyJS+99BLlypW7+AkBxNsbZm/OYnu7HIrDGBPgJk6cyG233Ua9evXYsGEDXbt2dTqkgJSYmMioUaNo3LgxDz/8sNPhXLIL3n3k7ggGKOSxnK42sNcnURljAk6nTp0YOXIkr732GiEhIU6HE7AmTJjA7t27WblyZa4c1uNiNYVq7leQx3I1XM8ZRAF9fRqdMcZRYWFhGWP1XHXVVUyaNMkSwgVs376dl19+mbvuuivXPsl9wTSmqoMARCRcVT/wT0jGGKclJyczduxYJk+ezNVXX83p06cJDQ11OqyAlpaWxpAhQyhWrBhTpkxxOpxs86puk54Q3KOYlsNj2ApV3eWb0IwxTti5cyf9+/dnw4YNDB8+nDfeeMNqB16YPn06a9asYdasWVSo4O0IQIHHq6TgfnDtU+AqXGMWCf8buyjYN6EZY/wtKSmJdu3acebMGb766it69+7tdEi5wv79+xk7diwdOnRg4MCBTodzWbztBZkG/IDrLqTduEZHfRXXeEjGmFwuLi6OYsWKERISwuzZs2nQoAHVqlW7+IkGVWXIkCGkpKTw/vvv54rxjS7E21tSrwLGqmosrqegTwJP4pr/wBiTi/388880bdo046nbTp06WUK4BNOmTWPp0qW88cYb1KlTx+lwLpu3SSERSH+W/aiIVHefW9YnURljfC4lJYVnn32Wtm3bEhwczNVXX+10SLnOjh07eOKJJ+jSpUuemUjI2+aj1cCdwGxgPrAUSAIuOhieMSbw7NixgwEDBhAREcHAgQOZOnUqJUqUcDqsXCUlJYUBAwZQtGhRZs6cmeubjdJ5e/fRnR6rTwNbgWLAHF8EZYzxrb1797J7927mz59Pnz59nA4nV3ruueeIiIjgyy+/pFKlSk6Hk2OyHCU1N7BRUo3xXlRUFKtWrWLAgAGAq3O5ePHiDkeVO3399dfceuutDBkyhOnTpzsdziW70CipFxo6+2M8pszMiqr+4/LCyz5LCsZcnKryySefMHz4cAB2795N6dKlHY4q99qzZw/NmzenZs2ahIeHU7hwYadDumTZGjobiPRRPMYYPzl48CAPPvggYWFh3HDDDXz00UeWEC5DUlISffv2JS0tjXnz5uXKhHAxF5pP4UV/BmKMyVkJCQlcffXVxMbG8uabbzJixAiCg+1Z0+xSVYYNG0ZERARfffVVnrj9NDPePtHcPqt93kzHaYzxn9jYWEqVKkXRokWZMGECLVu25IorrnA6rFxv0qRJzJgxg6effjpPP+ntVUdzJtNxlgcKAdE2HacxgUFVmT17NqNGjWLOnDnceuutToeUZyxevJhevXpx++2388UXXxAU5O0jXoEpu30KGc6fjlNEgoFngLjLD88Yc7l27drFgw8+yMqVK7npppto0KCB0yHlGb/99hv9+/enefPmzJkzJ9cnhIvJ1rtT1VTgZWBMzoZjjLlU06dPp2nTpqxbt45p06bxww8/UK9ePafDyhP27t1Lt27dKFWqFGFhYRQtWtTpkHzucqYF6gSk5VQgxpjsKVSoEO3bt2fatGlUrVrV6XDyjGPHjnHLLbdw+vRpVq9eTeXKlZ0OyS+87WiO4txnFooChYG8MdiHMblIfHw8zz//PHXr1uXhhx9m4MCBDBw4MM8MsxAITp8+TY8ePdizZw8rVqygadOmTofkN97WFAact34a+FNVT+VwPMaYC1i4cCHDhw8nOjqa0aNHA1gyyGHJycn07duX9evXM3/+fG666SanQ/IrbzuaV/k6EGNM1vbu3cuIESMICwujadOmzJ07l9atWzsdVp5z9uxZ7r77bpYuXcr06dO57bbbnA7J77xtPioJjACuxjUQXgZV7eyDuIwxHnbv3s13333HhAkTGDlyJAULFrz4SeaSpKWlcf/99/Pll18yadIkhgwZ4nRIjvC2+Wgermk3vwLO+C4cY0y6FStWsGXLFkaNGkW7du3Yt28fZcqUcTqsPCktLY2HH36Yjz/+mPHjxzNy5EinQ3KMt0mhFVBWVVN8GYwxxtVUNGrUKBYsWECDBg149NFHCQkJsYTgI6mpqQwZMoRZs2bx1FNP8c9//tPpkBzl7XMKPwMNfRmIMfldQkICL7zwAg0aNGDZsmW88sorbNq0iZCQEKdDy7POnj3LwIEDmTVrFi+88AIvv/xyvu+497amcB/wjYisAw557lDVl3I6KGPyowMHDvDqq69y2223MWHCBKpXr+50SHlaUlIS99xzD19++SWvvPIKTz31lNMhBQRvk8LLQDVgD+A5Z1/unaHHmACwceNGvvrqK8aPH0/dunX5888/qVGjhtNh5XlxcXH07t2b77//nsmTJ/PYY485HVLA8DYp3AXUV9UYXwZjTH4RE+B5M30AABtNSURBVBPDM888w6xZsyhbtiyPPvooFStWtITgB0eOHKFr165s2rSJOXPmcO+99zodUkDxtk9hF2CdzMZcpoSEBF566SXq1avHxx9/zOjRo4mMjKRixYpOh5Yv7Ny5k9atW7Nt2zYWLVpkCSET3tYUPgbCROQt/t6nYPMpGOOl5ORkpk6dSteuXXnttdfy7EQtgWj16tX07t2boKAgVq5caQ//ZcHbpPCo++cr521XwKv5FESkCzAF1/MOM1T1tfP2jwIGA2eBI8D9qrrXy/iMCUiqyrJly5g9ezaffvoppUqVYvv27ZQvX97p0PKV//znP9x///3UqlWLJUuWWDK+AK+aj1S1VhYvbxNCMPAO0BVoBPQXkUbnHfYr0EJVrwTmAxO8fxvGBJ7169fToUMHunXrxoYNG4iKigKwhOBHqampjBkzhgEDBnD99dcTHh5uCeEi/DVbREsgUlV3qWoy8DnQy/MAVf1BVRPcq/8FbAxgkysdP36cvn37ct1117FlyxamTp3Ktm3bqFmzptOh5SsnTpygW7duvPHGGzzyyCOsWLHCHgD0QnaHzs6gqt7cTF0FiPJYjwauu8DxDwBLs4hlKDAUsPu4TUBJSkoiJCSEEiVKsHv3bp5//nlGjx5N8eLFnQ4t39m4cSN9+/YlKiqKDz74gMGDBzsdUq6R3aGzKwGP4fqL3xuZPSKYaZIRkQFAC6BtZvtVdTowHVxzNHtZvjE+c+jQIV599VUWLFjA1q1bKV68OOvXr8/z0zYGIlXlvffeY+TIkVSoUIFVq1Zx/fXXOx1WrpLtobNF5EdgGa7O44uJxvXwW7qqwIFMrtkR+CfQVlWTvInNGKccO3aMN954g7feeoukpCQGDhxIUlISxYsXt4TggNjYWB566CHmzp1L165dmTNnDuXKlXM6rFzncqbjTAJqeXnsL0A9EakF7Mf1MNzdngeIyNXA+0AXVT18GXEZ43P79u2jcePGnD59mv79+/PCCy/YvMgO+vHHH/nHP/5BTEwMr776KmPGjLHEnE3e9imcP75RUaAbWbT7n09Vz4rIMGA5rltSZ6rqVvd1I1Q1DHgD11wN89wDUu1T1Z7evQ1jfO/48eOsXbuW7t27U716dZ544gnuuOMOGjdu7HRo+VZiYiLPP/88b7zxBnXr1iU8PJxrr73W6bByNVG9eLO8iMw6b9NpYBPwsZPNPC1atNCIiAinijf5xNGjR5k0aRJvvfUWKSkpxMTEUKpUKafDyvfWrVvHoEGD2L59O0OGDGHSpEmEhoY6HVauICIbVLVFZvu87VMYlLMhGRP4Dh8+zIQJE5g2bRpnzpyhb9++PPvss5YQHJaQkMCLL77IxIkTqVy5MkuXLqVLly5Oh5VnXLDRTURuEJHXs9j3moi08k1YxjgnLS0NcHVcTp06ldtvv50tW7Ywd+5cmjRp4nB0+dvSpUtp0qQJEyZMYNCgQWzZssUSQg67WE3haeDdLPb9iOtOoVtzMiBjnLJ161Zef/11zpw5w7x586hfvz7R0dFUqFDB6dDyvaioKEaPHs28efO44oor+OGHH2jXrp3TYeVJF+ueb4brttPMrASuydlwjPEvVWX16tX06NGDJk2a8OWXX1KjRo2M2oIlBGclJiby8ssv06BBAxYvXsxLL73Eb7/9ZgnBhy5WUygBFALOZLKvIGCPappcbfr06Tz00EOUK1eOF198kUcffZSyZcs6HVa+p6p8+eWXjB07ll27dtGnTx8mTpxoQ4X4wcWSwh9AZ2BRJvs6u/cbk2vEx8cze/Zs6tevT+fOnbn99ttJSUnh/vvvp2jRok6HZ4C1a9cyevRo1q5dS9OmTfn222/p2LGj02HlGxdrPpoEvC8it4tIEICIBInI7cB7wJu+DtCYnBAdHc24ceOoVq0aw4cPZ8GCBYBrxNJhw4ZZQggAf/zxB7fffjutW7dm9+7dzJgxg19//dUSgp9dsKagqp+KSEXgIyBERI4C5YBE4HlV/cwPMRpzWcaOHcu///1vVJXbb7+dUaNG2Xg4AWTfvn289NJLzJo1i9DQUF588UVGjRpFsWLFnA4tX7rocwqq+qaIzACuB8oCx4C1qnrK18EZkx2JiYl88cUX9OnTh9DQUBo0aMDIkSMZNmyYtUkHkOjoaF555RVmzJiBiDBixAiefvppm2/CYd4+vHYK1xAVxgSsPXv28P777zNjxgyOHj1KgQIFuPvuuxk0yJ69DCS7du1iwoQJzJo1C1Vl8ODBPPXUU1SrVu3iJxufu5wB8YwJCKdPn+auu+5iyZIliAg9e/Zk+PDh3HzzzU6HZjxs2rSJN954g88//5wCBQowaNAgnnrqKWrUqOF0aMaDJQWTK8XExBAREcGtt95KaGgoIsI///lPhg4dan9xBhBVZfny5UycOJHvvvuO0NBQRo0axeOPP07lypWdDs9kwpKCyTVSU1NZtmwZM2bMYPHixRQqVIjDhw9TrFgxwsLCnA7PeIiPj+fjjz9m6tSp/PHHH1SuXJnXXnuNoUOHUrp0aafDMxdgScHkCsuXL2fw4MFER0dTvnx5Ro0axeDBg+0OlQDz559/8t577zFr1ixiY2Np0aIFH3/8MXfeeSeFChVyOjzjBUsKJiDFx8ezYMECGjRoQMuWLalevTpNmjRh0qRJ9OzZ075gAkhKSgqLFy/mvffe49tvv6VAgQL06dOHESNGcP311+OeH8XkEpYUTMBIS0vjxx9/ZM6cOcyfP5/Tp08zfPhwWrZsScOGDVm61Ks5nYyf7Ny5k5kzZzJr1iwOHTpE1apVGT9+PIMHD6ZixYpOh2eyyZKCCRg33ngja9eupXjx4vTr149BgwZxww03OB2W8XDq1Cnmz5/PzJkzWbNmDcHBwXTv3p2hQ4fSpUsXgoODnQ7RXCZLCsYRe/fuZe7cuaxYsYJly5ZRoEABHnzwQUaMGEHPnj1t2IkAkpyczPLly/nkk08ICwsjMTGRK664gtdee417773X7iLKYywpGL85cuQIc+fO5bPPPiM8PByAVq1acejQIapUqcLAgQMdjtCkS05OZuXKlXzxxRcsWrSI2NhYypYty/3338+AAQNo1aqV9RXkUZYUjE/FxMSgqlSuXJnNmzczfPhwmjZtyssvv8xdd91F7dq1nQ7RuJ06dYply5bx1Vdf8c0333Dq1ClKlixJ79696du3L507d6ZgwYJOh2l8zJKCyXF79+5l4cKFzJ8/nzVr1jBq1CgmTpxI27Zt2bp1K40aNXI6ROMWGRnJkiVL+Prrr1m1ahUpKSmUL1+eO++8k969e9OpUye70yufsaRgcoyq0rFjR77//nsArrzySl544QXuvPNOAAoUKGAJwWHx8fGsWrWKZcuWsXTpUv766y8AGjZsyOOPP06PHj1o3bq1dRjnY5YUTLYkJyezevVqFi9ezI4dO1i6dCkiQtu2benSpQu9e/emXr16ToeZ76WmprJhwwa+++47vv32W37++WdSUlIoUqQI7du3Z+TIkXTt2pU6deo4HaoJEJYUzCX56aefeOutt1i+fDlxcXGEhITQsWNHEhISKFq0KM8995zTIeZrqsrWrVv5/vvv+eGHH/jxxx+JjY0FXDW3kSNH0rlzZ2688UYKFy7scLQmEFlSMFlKTk4mPDw8Y4iJOnXqEBUVRXh4OP3796dHjx60b9+e0NBQp0PNt1JTU9myZQurVq1i1apV/PTTTxw9ehSA2rVrc8cdd9ChQwfat29PhQoVHI7W5AaWFMw5Tp48yezZs/n2229ZtWoV8fHxFChQgKuvvpo6derQr18/7r77brsd0SEnT55k/fr1rF27ljVr1vDf//6XU6dc813VqFGDbt260bZtW9q3b28TCplssaSQj6kqf/31Fz/++CMVKlSgZ8+epKWlMWrUKGrXrs2AAQPo0qULN998MyVKlABcncXGP5KTk9m6dSvr169n3bp1rF+/nm3btqGqiAiNGzfm7rvv5oYbbuDGG2+0JGByhP0Pz4dmz57NihUrWLVqFQcOHACgb9++9OzZk9KlSxMdHU2lSpUcjjJ/SUhIYMuWLWzatIlff/2VDRs28Ntvv5GcnAxAuXLlaNmyJf369aNVq1a0bNmSkiVLOhy1yYssKeRhCQkJREREEB4ezqFDh5g0aRLgSgo7duygXbt2tG3blrZt29KgQYOM8ywh+E5qaiq7d+9my5YtbNmyhc2bN7N582Z27txJWloaACVLlqR58+Y89thjXHPNNVx77bXUqlXLmuyMX1hSyCNSU1Mz7i2fMWMG7777Lps3byY1NRWAxo0bZxyzaNEiSpQoYV8yPpSYmEhkZCQ7duzgjz/+YNu2bWzfvp0//viDM2fOZBxXp04dmjZtSr9+/WjWrBnNmjWjZs2a9tkYx1hSyIVSUlLYsWNHRjPDxo0b2bhxI3/++SeVK1fm7NmzlClThnHjxnHddddx/fXXU65cuYzzrdkhZyQnJ7Nnzx4iIyOJjIxk586d7Ny5kz///JO9e/dm/OUPUL16dRo1akS7du1o0qQJTZo0oVGjRjZJkAk4oqpOx5BtLVq00IiICKfD8BlV5dChQ/z+++/8/vvv9OrVizp16vDJJ59w7733AlCkSBGaNWvGNddcw9ixY6latarDUecdZ8+e5cCBA+zdu5e9e/eyZ88edu/eze7du9m1axdRUVHnfPEXL16c+vXrZ7yuuOIKrrjiCurXr29f/iagiMgGVW2R2T6rKQSAs2fPsmfPHkJDQ6lUqRKRkZHcd999bNu2jRMnTmQcV758eerUqcPNN9/MJ598QrNmzWjQoIENSZANiYmJxMTEEBMTw/79+9m/fz/R0dFER0cTFRVFVFQUBw4cyGh+S1epUiVq1qxJmzZtqF27NrVr16ZevXrUrVuX8uXLW7OPyfX8lhREpAswBQgGZqjqa+ftDwHmANcAx4B+qrrHX/H5kqoSGxtLSkoKFSpU4PTp04wbN46//vqLnTt3smfPHs6ePctLL73Es88+S8mSJQkKCqJfv340atSIxo0b07RpU8qXLw9AlSpVuOeeexx+V4Hl7NmzHD9+nGPHjnHkyBGOHj3KkSNHOHLkCIcPH+bw4cMcOnSIgwcPcvDgwYynfD2FhIRQtWpVqlWrRrt27ahevTo1atSgRo0aGctFihRx4N0Z4z9+aT4SkWDgT6ATEA38AvRX1W0exzwCXKmqD4nIXcBtqtrvQtd1uvlIVTl9+jSHDh3i0KFDFC5cmObNmwPw8MMPZzQxREdHExcXx5AhQ5g+fTqpqalUqVKFypUrU7du3Yy/NFu3bs0VV1zh2Ptxkqpy5swZTp06xalTp4iNjeXkyZOcPHmS2NhYYmNjOXHixDmv9CRw7NgxTp48meW1S5YsSYUKFahQoQKVKlWiYsWKVKxYkcqVK1OpUiUqV65MlSpVKFOmjP2lb/KFQGg+aglEquoud0CfA72AbR7H9AJecC/PB94WEVEfZK39+/ezdu1aChQokPFQVnJyMiEhIdx8882uAObPZ/fu3cTFxWV8UVWpUoXx48cD0KZNG3755ReSkpIyrtulS5eMeYQ3btyIqtKwYUM6depEjRo1aNmyJQDBwcEcPHgwp99WjlFVUlNTSUlJISUlheTkZJKTk0lKSsr46flKTEzM+HnmzJlzXgkJCZw5c4bTp0+f84qLiyM+Pp64uLiM5fObas4XFBRE6dKlM15lypShXr16lC1bljJlylC+fHnKli1LuXLlKF++POXKlaNcuXKEhIT46V/OmNzPX0mhChDlsR4NXJfVMap6VkROAmWBozkdzH/+8x/Gjh37t+2FChWiUaNGFCxYkB07dmQMH1CgQAEKFChA6dKl2bp1K0FBQcTGxlKrVi1CQkIoXLgwhQsXJjg4mAEDBgCcM0Jo+l+zGzduZNq0aRnb0/Odqv5t2fOVlpb2t+W0tLS/Lae/UlNTSU1NPWfZ83X27NmMn2fPns348vdczimFChWiSJEihIaGZryKFStG2bJlqV69OsWLF894lShRghIlSlC8eHFKlSpFyZIlKVmyJKVLl6ZUqVIUK1bM/pI3xsf8lRQy+598fg3Am2MQkaHAUHDd5pcdd955J8HBwezbty+jmSIlJYXU1FRCQkIyhhY+e/Zsxhdr+nJkZCSpqakZX9IJCQkkJCSc88V+MZ5fbOnLIvK35fRXUFDQOcvp6+cvBwcHExQURIECBQgODs5YT18ODg7OSHDBwcEULFgwY1vBggUzfnq+ChUqRKFChShYsCAhISEZ64ULFyYkJOScpBgSEkKRIkUoUqQIhQsXpmjRotYJbkwu46+kEA1U81ivChzI4phoESkAlASOn38hVZ0OTAdXn0J2gqlZsyajR4/OzqnGGJOnBfmpnF+AeiJSS0QKAXcBYecdEwakz9x+B/C9L/oTjDHGZM0vNQV3H8EwYDmuW1JnqupWEXkJiFDVMOBD4GMRicRVQ7jLH7EZY4z5n1z9RLOIHAH2ZvP0cvigEzvA2XvOH+w95w+X855rqGr5zHbk6qRwOUQkIqv7dPMqe8/5g73n/MFX79lffQrGGGNyAUsKxhhjMuTnpDDd6QAcYO85f7D3nD/45D3n2z4FY4wxf5efawrGGGPOY0nBGGNMhnyZFESki4jsEJFIERnndDy+ICLVROQHEdkuIltF5DH39jIi8q2I7HT/LO10rDlJRIJF5FcR+dq9XktE1rnf71z3E/V5hoiUEpH5IvKH+7O+Ph98xo+7f6e3iMhnIlI4r33OIjJTRA6LyBaPbZl+ruIy1f19tllEml9O2fkuKbjndngH6Ao0AvqLSCNno/KJs8BoVW0ItAIedb/PccB3qloP+M69npc8Bmz3WH8dmOR+vyeABxyJynemAMtUtQFwFa73nmc/YxGpAowAWqhqE1wjJNxF3vucZwNdztuW1efaFajnfg0FpnEZ8l1SwGNuB1VNBtLndshTVDVGVTe6l+NwfVlUwfVeP3If9hHQ25kIc56IVAW6AzPc6wK0xzU/B+S991sCuAnXEDGoarKqxpKHP2O3AkAR98CZRYEY8tjnrKo/8fcBQbP6XHsBc9Tlv0ApEamU3bLzY1LIbG6HKg7F4hciUhO4GlgH/J+qxoArcQAVnIssx00GxgBp7vWyQKyqnnWv57XPujZwBJjlbjKbISKh5OHPWFX3AxOBfbiSwUlgA3n7c06X1eeao99p+TEpeDVvQ14hIsWAL4GRqnrK6Xh8RUR6AIdVdYPn5kwOzUufdQGgOTBNVa8GTpOHmooy425H7wXUAioDobiaT86Xlz7ni8nR3/P8mBS8mdshTxCRgrgSwn9UdYF786H0qqX752Gn4sthNwA9RWQPribB9rhqDqXczQyQ9z7raCBaVde51+fjShJ59TMG6AjsVtUjqpoCLABak7c/53RZfa45+p2WH5OCN3M75Hru9vQPge2q+qbHLs95KwYCi/wdmy+o6lOqWlVVa+L6TL9X1XuAH3DNzwF56P0CqOpBIEpErnBv6oBr3vM8+Rm77QNaiUhR9+94+nvOs5+zh6w+1zDgH+67kFoBJ9ObmbIjXz7RLCLdcP0VmT63w8sOh5TjRORGYDXwO/9rY38aV7/CF0B1XP/B+qrq32a4y81EpB3whKr2EJHauGoOZYBfgQGqmuRkfDlJRJrh6lgvBOwCBuH6Yy/PfsYi8iLQD9cddr8Cg3G1oeeZz1lEPgPa4Roe+xDwPLCQTD5Xd3J8G9fdSgnAIFWNyHbZ+TEpGGOMyVx+bD4yxhiTBUsKxhhjMlhSMMYYk8GSgjHGmAyWFIwxxmSwpGBMPiUiS0VkoHv5PhH52emYjPMsKRjHicgeEenodBx5mYi8ICKfeG5T1a6q+lFW55j8yZKCMQHCY5gGYxxjScEElPRmDBGZKCInRGS3iHT12F9GRGaJyAH3/oUe+4a4Jxo5LiJhIlLZY5+KyCPuCUriRGS8iNQRkbUickpEvvCcmEVEeojIJhGJFZFwEbnyAjGriIwQkV0iclRE3hCRII/994trApwTIrJcRGqcd+6jIrIT2One1tg9icpxETkkIk+7tweJyDgR+UtEjrljLuPeV9N9rYEiss8dxz/d+7rgepq9n4jEi8hv7u0/isjgLN5TA48YdojInV5/iCZ3U1V72cvRF7AH6Ohevg9IAYbgGobkYVyDe6U/fb8EmAuUBgoCbd3b2wNHcQ0IFwK8BfzkUYbiGiOmBNAYSMI1UUltoCSu8XMGuo9tjmuwsevcMQx0xxiSRfyKa+ydMriGIPgTGOze1xuIBBriGtX0GSD8vHO/dZ9bBCiOa0jo0UBh9/p17mNHAv/FNeBZCPA+8Jl7X033tT5wX+cq93ts6N7/AvDJeXH/6BHnfcDP7uVQXEMxD+J/I7EeBRo7/btiLz/8f3Q6AHvZK5OkEOmxr6j7y64iUAnXOE6lM7nGh8AEj/Vi7uRS072uwA0e+zcAYz3W/w1Mdi9PA8afd/0d6Qkok7IV6OKx/giuGbIAlgIPeOwLwjU+TQ2Pc9t77O8P/JpFOduBDh7rldzvsYBHUqjqsX89cJd7+VKSQj9g9XnHvg887/Tvir18/7LmIxOIDqYvqGqCe7EYruGBj6vqiUzOqQzs9TgvHjjGuZONHPJYPpPJejH3cg1gtLvpKFZEYt1lVyZrnpOc7PU4tgYwxeM6x3GNf18li3OrAX9lUUYN4CuPa20HUoH/8zjmoMdygsd7uhQ1gOvOe//34ErMJo+zji2Tm0QBZUSklLqmnfR0ANeXGQDimoGsLLA/m+W8rJc2em41YKt7uTr/G88+/Vr/ucC5nqNSRuGqLWQV1/2quub8HeKaXe9CLmXkyyhglap2uoRzTB5hNQWTa6hrjPilwLsiUlpECorITe7dnwKDRKSZiIQArwDrVHVPNor6AHhIRK5zj1EfKiLdRaT4Bc550h1TNeAxXP0eAO8BT4lIYwARKSkifS9wna+BiiIyUkRCRKS4iFznca2X0zuqRaS8iHg7v/ghoKZnB/hFYqgvIve6/40Lisi1ItLQy7JMLmZJweQ29+JqR/8DV2fwSABV/Q54FtdMczFAHVyT7VwydY1FPwTXGPUncHUU33eR0xbh6qfYhKsz/EP3tb4CXgc+F5FTwBYynz4yvew4oBNwK66moJ3Aze7dU3B1lq8QkThcnc7XZXadTMxz/zwmIhsvdKA7hs64/v0OuON4HVfntsnjbD4FYy6TiChQT1UjnY7FmMtlNQVjjDEZLCkYY4zJYM1HxhhjMlhNwRhjTAZLCsYYYzJYUjDGGJPBkoIxxpgMlhSMMcZk+H/ALGSpz0QzGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimate is center=0.9891395234236353, spread=0.0\n",
      "Lorenz distance is 42.691225404510476\n",
      "Average MPC for all consumers is 0.099\n",
      "Average MPC in the top percentile of W/Y is 0.068\n",
      "Average MPC in the top decile of W/Y is 0.071\n",
      "Average MPC in the top quintile of W/Y is 0.072\n",
      "Average MPC in the second quintile of W/Y is 0.074\n",
      "Average MPC in the middle quintile of W/Y is 0.074\n",
      "Average MPC in the fourth quintile of W/Y is 0.076\n",
      "Average MPC in the bottom quintile of W/Y is 0.199\n",
      "Average MPC in the top percentile of y is 0.075\n",
      "Average MPC in the top decile of y is 0.078\n",
      "Average MPC in the top quintile of y is 0.084\n",
      "Average MPC in the second quintile of y is 0.114\n",
      "Average MPC in the middle quintile of y is 0.123\n",
      "Average MPC in the fourth quintile of y is 0.082\n",
      "Average MPC in the bottom quintile of y is 0.093\n",
      "Average MPC for the employed is 0.097\n",
      "Average MPC for the unemployed is 0.131\n",
      "Average MPC for the retired is nan\n",
      "Of the population with the 1/3 highest MPCs...\n",
      "61.057% are in the bottom wealth quintile,\n",
      "28.655% are in the second wealth quintile,\n",
      "7.285% are in the third wealth quintile,\n",
      "2.207% are in the fourth wealth quintile,\n",
      "and 0.796% are in the top wealth quintile.\n",
      "\n",
      "These results have been saved to ./Code/Results/PYbetaPointIndNetWorth.txt\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Solve for the $\\beta-Point$ (do_param_dist=False) for speed\n",
    "'''\n",
    "This options file establishes the second simplest model specification possible:\n",
    "with heterogeneity, no aggregate shocks, perpetual youth model, matching net worth.\n",
    "'''\n",
    "do_param_dist = False          # Do param-dist version if True, param-point if False\n",
    "do_lifecycle = False          # Use lifecycle model if True, perpetual youth if False\n",
    "do_agg_shocks = False         # Solve the FBS aggregate shocks version of the model\n",
    "do_liquid = False            # Matches liquid assets data when True, net worth data when False\n",
    "\n",
    "\n",
    "os.chdir(path_to_models)\n",
    "exec(open('cstwMPC_MAIN.py').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAFNCAYAAABmLCa9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3gUZdfH8e9JCC006UpHegkBQgdFpIOISJUiKKCCiD4WFBVQUbE9isJrwUcQQ0dRpJmEJl2qiPQSIDQhtFDSz/vHJBggJJtkN7Ob3J/rykV2Z3bmlyiHKfecW1QVwzAM48687A5gGIbh7kyhNAzDSIUplIZhGKkwhdIwDCMVplAahmGkwhRKwzCMVJhCabiUiJQXERWRHCmsoyJSyYn7bCEi+5y4vaUi8njC9wNFZK0Tt91XRIKctT3DNUyhzIZE5DURWXLLewfu8F5vJ+97lYgMzsDnx4lIjIhEJHztF5FJInJ34jqqukZVqzq4rcDU1lPVDqr6fXozJ9nfbf9oqOoMVW2b0W0brmUKZfb0O9BMRLwBRKQk4APUu+W9Sgnrups5qpofKAw8ApQEtiYtls4gFvN3xDCFMpvajFUY/RNe3wesBPbd8t4hVT0JICLVRCRYRM6LyD4R6Zm4MRHpJCLbReSyiBwXkXHJ7VRE3gVaAJNE5IqITEqyuHXCEewFEZksIpLaD6GqMar6N9ALOAu8mLCfliISlmS/o0TkRMIR6D4ReVBE2gOjgV4JWf5MWHeViLwrIuuAa0DFZI6CRUS+EJFLIrJXRB5MsiBURFoneZ30qDXxH52LCftscuupvIg0FZHNCdveLCJNkyxbJSLviMi6hJ8lSESKpvZ7MjLOFMpsSFWjgU1YxZCEP9cAa29573cAEfEFgoGZQHGgD/B/IlIzYd2rwACgENAJeEZEuiaz39cT9vOsquZT1WeTLO4MNADqAD2Bdmn4eeKAX7CK8E1EpCrwLNAg4Si0HRCqqsuA97COTvOpap0kH+sPDAXyA0eT2WUj4DBQFBgL/CQihR2Imvi7LZSwzw23ZC0MLAY+B4oA/wUWi0iRJKs9BgzC+u+QE3jJgf0aGWQKZfa1mn//4rbAKmBrbnlvdcL3nbGKy1RVjVXVbcCPQHcAVV2lqn+paryq7gRmAfenMc8EVb2oqsewjm79U/vALU5inYrfKg7IBdQQER9VDVXVQ6lsa5qq/p3ws8Yks/wf4LOEI9o5WEfindKYNzmdgAOq+kPCvmcBe4GHkqwzVVX3q+p1YC5p/z0Z6WAKZfb1O9BcRO4CiqnqAWA90DThvVr8e6pYDmgkIhcTv4C+WNcGEZFGIrJSRM6KyCXgaayjrbQ4neT7a0C+NH6+FHD+1jdV9SDwPDAO+EdEZovIPals63gqy0/ozd1kjgKpbdMR93D7EexRrJ8tUUZ/T0Y6mEKZfW0ACmKdYq4DUNXLWEdmQ4GTqnokYd3jwGpVLZTkK5+qPpOwfCawECijqgWBr4A7XWN0eruqhBsuD2EdEd++Q9WZqtocq+Ar8EEqWVLLWOqWa6hlsX5vYF2GyJtkWck0bPdkQsakygInUvmc4WKmUGZTCaduW4D/cHOBWZvwXtK73YuAKiLSX0R8Er4aiEj1hOX5gfOqGikiDbGuo93JGaCiM36GhBzVsU71S2Jd07t1naoi0kpEcgGRwHWs0/HELOXTcWe7OPBcwv57ANWBxKFVO4DeCcsCSLg8keAsEM+df/4lWL/nx0Qkh4j0Ampg/f4NG5lCmb2txvpLn3QA9ZqE924USlWNANoCvbGOek5jHZXlSlhlGPC2iEQAY7Cund3JRKB7wt3tz9OZu5eIXAEuYh3JhgP1E+/Q3yIXMAE4l5C7ONbdboB5CX+Gi8i2NOx/E1A5YZvvAt1VNTxh2ZvAvcAF4C2so20AVPVawvrrEi5hNE660YRtdMa6ex8OvAJ0VtVzachmuICYxr2GYRgpM0eUhmEYqXBZoRSR70TkHxHZdYflIiKfi8hBEdkpIvVclcUwDCMjXHlEOQ1on8LyDljXeSpj3WX90oVZDMMw0s1lhVJVfyeZcW1JPAxMV8tGoJCzn9U1DMNwBjuvUZbi5oG9Ydw8sNYwDMMt3LFHYCZIbkBysrfgRWQo1uk5vr6+9atVq+bKXIZheKg4jWP/uf3kzZmXcgVvHbufsq1bt55T1WLJLbOzUIYBZZK8Ls2/TzfcRFW/Ab4BCAgI0C1btrg+nWEYHuV6zHXaz2hP9PFo5vaaS6cqaXv8XkSSa4AC2HvqvRAYkHD3uzFwSVVP2ZjHMAwPFRMXQ495PVhzdA2fNJvOseWdOJ/SHZI0cuXwoFlYzxNXFZEwEXlSRJ4WkacTVlmC1arqIDAF6+kOwzCMNFFVBv4ykMUHFvNlpy/Rv/owbBhcuuS8fbjs1FtV+6SyXIHhrtq/YRjZg4jQslxL/Ir78VTAU3QeB5UqQYUKztuHndconSYmJoawsDAiIyPtjpLl5M6dm9KlS+Pj42N3FMO4TejFUMoXKs+Q+kMAiI6GVatgwADn7idLFMqwsDDy589P+fLlcWAGAcNBqkp4eDhhYWFUcOY/z4bhBFO2TuG5Zc+xafAm/Er4AbBxI1y9Cm3aOHdfWaJQRkZGmiLpAiJCkSJFOHv2rN1RDOM2Xat15fjl49QsVvPGe0FB4OUFDzzg3H1lmaYYpki6hvm9Gu5m3bF1RMdFU8y3GG8/8DbeXt43lgUHQ6NGUKiQc/eZZQqlO3j33XepWbMmfn5++Pv7s2nTJmJiYnj11VepXLkytWrVomHDhixduhSA8uXLU7t2bfz9/fH392f9+vU2/wSG4d6WHFhCy+9b8vbqt29bduECbNni/NNuyCKn3u5gw4YNLFq0iG3btpErVy7OnTtHdHQ0b775JqdOnWLXrl3kypWLM2fOsHr16hufW7lyJUWLmhlHDSM1a46u4dG5j+JXwo9Xmr1y2/IVKyA+3hRKt3bq1CmKFi1KrlxW0++iRYty7do1pkyZwpEjR268X6JECXr27JnSpgzDuMW2U9voPKsz5QuVZ1nfZRTIVeC2dYKDIX9+69Tb2cypt5O0bduW48ePU6VKFYYNG8bq1as5ePAgZcuWpUCB2/+jJnrggQfw9/enkSv+6xpGFrD33F7aBbajUO5CBPULophvso9jExwMLVuCK0ayZckjypYtW972Xs+ePRk2bBjXrl2jY8eOty0fOHAgAwcO5Ny5c3Tv3v2mZatWrUp1n/ny5WPr1q2sWbOGlStX0qtXL0aPHp3q58ypt2Hc2dGLR2nzQxu8xZuQ/iGUKVgm2fUOHYLDh+GFF1yTI0sWSrt4e3vTsmVLWrZsSe3atfn66685duwYERER5M+f3+54huFRzlw5Q+sfWnMl+gqrHl9F5SKV77hucLD1pyuuT0IWLZQpHQHmzZs3xeVFixZ16AjyVvv27cPLy4vKla3/mDt27KBq1arUrVuX5557jq+//pqcOXNy6tQpli9fTr9+/dK8D8PITl4KfomTEScJ7h9MnZJ1Ulw3OBjKlIEqVVyTJUsWSjtcuXKFESNGcPHiRXLkyEGlSpX45ptvKFCgAG+88QY1atQgd+7c+Pr68vbbtw9tMAzjZl90+IKn6z9N0zJNU1wvLs66492tG7hq2K/HTVebXD/KPXv2UL16dZsSZX3m92tklqjYKD5c9yEvNX2JPD55HPrMpk3QuDHMmgW9e6d/3yKyVVUDkltm7nobhuE2VoWuYtzqcSw/stzhzwQFWUeSrVu7Lpc59TYMw220q9SOPcP3UKWI4xcbg4Ohbl1w5eARc0RpGIatVJXXl79O0KEggDQVyYgI2LDBdXe7E5lCaRiGrd5d8y7vrX2PZQeXpfmzq1dDbKwplIZhZGGT/pjEmyvfZECdAXzc9uM0fz44GHLnhmbNXBAuCVMoDcOwxQ9//sCIpSPoWq0r/+vyP7wk7eUoOBjuv98qlq5kCqWTeHt74+/vT506dahXr16qLdMuXrzI//3f/6V5P9OmTaNYsWL4+/tTo0YNpkyZkuL6CxcuZMKECSmuExoaysyZM9OcxTDS65e9vzDol0E8WOFBZj06ixxeab+vHBYGe/a4/rQbsC6ketJX/fr19Va7d+++7b3M5uvre+P7ZcuW6X333Zfi+keOHNGaNWumeT9Tp07V4cOHq6rqmTNntGjRonr69Ok0byeplStXaqdOne643B1+v0bWsfzwcs35Tk5tNKWRRkRFpHs7332nCqp//umcXMAWvUPdMUeULnD58mXuuuuuG68/+ugjGjRogJ+fH2PHjgXg1Vdf5dChQ/j7+/Pyyy9z5coVHnzwQerVq0ft2rX55ZdfUt1P8eLFuffeezl69Cjnz5+na9eu+Pn50bhxY3bu3AlYR6DPPvssYDX+eO6552jatCkVK1Zk/vz5N7KsWbMGf39/Pv30U2f/Ogzjhk1hm+gyqwtVilRhSd8l5MuZL93bCg6GEiWgdm0nBrwDM47SSa5fv46/vz+RkZGcOnWKFStWABAUFMSBAwf4448/UFW6dOnC77//zoQJE9i1axc7duwAIDY2lgULFlCgQAHOnTtH48aN6dKlS4pTMRw+fJjDhw9TqVIlxo4dS926dfn5559ZsWIFAwYMuLHtpE6dOsXatWvZu3cvXbp0oXv37kyYMIGPP/6YRYsWueaXYxgJjlw8QukCpQnqF0ThPIXTvZ34eAgJgbZtXffYYlJZrlA+/zwkUx8yxN8fPvss5XXy5MlzozBt2LCBAQMGsGvXLoKCgggKCqJu3bqA9Uz4gQMHKFu27E2fV1VGjx7N77//jpeXFydOnODMmTOULFnytn3NmTOHtWvXkitXLr7++msKFy7M2rVr+fHHHwFo1aoV4eHhXEpmBviuXbvi5eVFjRo1OHPmTHp+HYaRZrHxseTwykHvWr3pVr0bOb1zZmh7O3fC2bNWocwMWa5QuoMmTZpw7tw5zp49i6ry2muv8dRTT920Tmho6E2vZ8yYwdmzZ9m6dSs+Pj6UL1+eyMhIJk+efOOGzZIlSwDo1asXkyZNuunzmswz+8kdjSZ2Wr/TZwzD2U5GnOTB6Q/yQesP6FK1S4aLJFiPLYJrH1tMKssVytSO/DLD3r17iYuLo0iRIrRr144333yTvn37ki9fPk6cOIGPjw/58+cnIiLixmcuXbpE8eLF8fHxYeXKlRw9ehSA4cOHM3z48FT3ed999zFjxgzefPNNVq1aRdGiRVPsrJ7UrVkMw5ly58hN2YJluTvf3U7bZnAw1KwJ99xz+7JPP/2UK1eu8Oabbzptf1muUNol8RolWEdq33//Pd7e3rRt25Y9e/bQpEkTwOqEHhgYyL333kuzZs2oVasWHTp0YNSoUTz00EMEBATg7+9PtWrV0rT/cePGMWjQIPz8/MibNy/ff/+9w5/18/MjR44c1KlTh4EDB/KCq9pEG9lKRFQEPt4+FM5TmN/6/ea07V6/DmvWwDPPJL/8xx9/JEeOHE4tlKbNmpEq8/s10ioyNpKOMzqS0zsnS/suder88MHB1rXJxYshmVldKFWqFG3btmXq1Klp2q5ps2YYRqaJiYuh1/xerApdRX+//k4tkmAVSh8f64mcW0VGRnLy5EkqVKjg1H2aU2/DMJwmXuN5YuETLNy3kMkdJ9PXr6/T9xEUZD3b7et7+7Jjx44BOL1QmiNKwzCcQlV5bulzBO4MZPwD4xnWYJjT93HmDPz5550fW7xw4QIlS5akYsWKTt2vOaI0DMMpxqwcw+TNk3mpyUuMbpH6VM3psTyh8fmdCmWjRo04deqU0/drjigNw8iw/274L+PXjGdw3cF82OZDp1+XTBQcDHfdBfXquWTzd2QKpWEYGbJgzwJeDHqRHjV68FXnr1xWJFWtQtm6NXh7J7/Oa6+9xogRI5y+b1MonSguLo66devSuXPnG+8NHDiQChUq4O/vj7+/f7LPX69atYqCBQtSt25dqlatyn333WeeuzY8RuuKrXmjxRsEdgvE2+sOFcwJ9u6FEydSbqu2fPly9u/f7/R9m2uUTjRx4kSqV6/O5cuXb3r/o48+onv37il+tkWLFjeK444dO+jatSt58uThwQcfzFCmG22ivMy/iYZzbQzbSO3itcmfKz/vtHrH5ftLfGwxpUJ55MgR6rngvNz87XGSsLAwFi9ezODBgzO8LX9/f8aMGXPjee6zZ8/y6KOP0qBBAxo0aMC6detuvN+mTRvq1avHU089Rbly5Th37hyhoaFUr16dYcOGUa9ePY4fP05QUBBNmjShXr169OjRgytXrgCwdetW7r//furXr0+7du1cciHcyHrCr4XT9oe2PL/s+UzbZ3AwVKoE5csnv/zKlSucO3fO6UODwBRKp3n++ef58MMPkz1ye/311/Hz8+OFF14gKirKoe3Vq1ePvXv3AjBy5EheeOEFNm/ezI8//nijGL/11lu0atWKbdu28cgjj9wYQwawb98+BgwYwPbt2/H19WX8+PGEhISwbds2AgIC+O9//0tMTAwjRoxg/vz5bN26lSeeeILXX3/dCb8NI6srkrcIs7vPZnyr8Zmyv+hoWLUq9aNJgPJ3qqQZkCVPvVtOa5nqOp2rdOalpi/dWH+g/0AG+g/k3LVzdJ9782nyqoGrUtzWokWLKF68OPXr12fVqpvXff/99ylZsiTR0dEMHTqUDz74gDFjxqSaL+mjpSEhIezevfvG68uXLxMREcHatWtZsGABAO3bt7+pWXC5cuVo3LgxABs3bmT37t00S5iBKTo6miZNmrBv3z527dpFm4T/++Li4rj7buc1LjCynv3h+zkQfoBOVTrRsXIyzw+6yMaNcPVqyoUyOjqaBg0aUKWK49PdOipLFsrMtm7dOhYuXMiSJUuIjIzk8uXL9OvXj8DAwBuFJ1euXAwaNIiPP3Zsprnt27ffeL46Pj6eDRs2kCdPnpvWSek5fd8kjy2oKm3atGHWrFk3rfPXX39Rs2ZNNmzY4FAmI3s7dukYrae3Jk7jODDiAHl98mbavoODrTvdDzxw53Xq16/PH3/84ZoAd5ojwl2/3HXOnES3zj9z8uRJVVWNj4/XkSNH6qhRo1L9zJ9//qnly5fXkJAQVVXt06ePfvjhhzeWb9++XVVVhw0bphMmTFBV1d9++00BPXv27G3z8fzzzz9apkwZPXDggKqqXr16Vfft26dRUVF677336vr161VVNTo6Wnft2nVbPnf6/Rr2OHPljFb5oooWeL+Abju5LdP337ChapMmrt0HZs4c+/Tt25fatWtTu3Ztzp07xxtvvJHsemvWrLkxPGj48OF8/vnnN+54f/7552zZsgU/Pz9q1KjBV199BcDYsWMJCgqiXr16LF26lLvvvpv8+fPftu1ixYoxbdo0+vTpc2NOnb1795IzZ07mz5/PqFGjqFOnDv7+/qnOHmlkP5ciL9E+sD3HLx1n8WOLqXt33Uzd/4ULsGVL6rMtDhw4kN69e7smxJ0qqLt+ufsRZWaKjIzUmJgYVVVdv3691qlTxyX7ya6/X0P1avRVbf5dc/V520eXHlhqS4b5863ZFtesSXk9Pz8/7dy5c7r3QwpHlOYapQc7duwYPXv2JD4+npw5c6Y6x7dhpEV0XDSPzn2U9cfXM/vR2bSv1N6WHMHBkD8/NGp053VUldDQUO5PrveaE5hC6cEqV67M9u3b7Y5hZEFx8XH0X9CfZQeXMeWhKfSo2cO2LMHB1k0cH587r3PhwgUuX77skqFBYMZRGoaRjG2ntrFgzwI+bvMxg+tl/CGK9Dp0CA4fTv36ZOIYSlcMNgcXF0oRaS8i+0TkoIi8mszysiKyUkS2i8hOEUn3wCz1sCktPIX5vWZPDUo1YPfw3bzY9EVbcwQHW3+mVih9fHx49NFHqVGjhktyuKxQiog3MBnoANQA+ojIrT/FG8BcVa0L9Ab+Lz37yp07N+Hh4eYvtZOpKuHh4eTOndvuKEYm+WDtB8zeNRuASoUr2ZzGKpRlykBqY8j9/PyYP38+VatWdUkOV16jbAgcVNXDACIyG3gY2J1kHQUS51QtCJxMz45Kly5NWFgYZ8+ezUBcIzm5c+emdOnSdscwMkFsfCxLDy6lTMEy9K7lomE2aRAXBytWQLdukFrntri4OLzv1HvNCVxZKEsBx5O8DgNuvW81DggSkRGAL5DsdOYiMhQYClC2bNnblvv4+Ljs2oRhZAeqSg6vHCztu5QcXu5xj3fLFrh40ZpxMTVdunQhKiqKkJAQl2Rx5TXK5P4NuPXcuA8wTVVLAx2BH0Tktkyq+o2qBqhqQLFixVwQ1TCyr1/3/cqD0x/kUuQl8vjkwcc7hdvLmSg42DqSdKTT4OHDhylYsKDLsriyUIYBZZK8Ls3tp9ZPAnMBVHUDkBso6sJMhmEksSp0FT3m9SAiOsJlncnTKygI6taFoqlUhMQxlK4aGgSuLZSbgcoiUkFEcmLdrFl4yzrHgAcBRKQ6VqE0FxoNIxNsObmFLrO6cG/he1nadykFchVI/UOZJCICNmxI/W43wJkzZ4iMjHTp5TeXFUpVjQWeBX4D9mDd3f5bRN4WkS4Jq70IDBGRP4FZwEA1t64Nw+V2n91N+8D2FMlbhKB+QRTN614ncqtXQ2ysY4XS1WMowcVP5qjqEmDJLe+NSfL9bqCZKzMYhnGzIxeO0OaHNvh4+xDSP4RSBUrZHek2wcGQJw80c6A6FC5cmJEjR1KrVi2X5XGP21uGYWSKUxGnaPNDG67HXGf1wNXcW/heuyMlKzgY7rsPHBnCW7VqVT777DOX5jGPMBpGNnHh+gXaBbbj9JXTLO27lNolatsdKVlhYbBnj2On3QD//PMPkZGRLs1kCqVhZBO5cuTi3sL38kvvX2hUOoVWPDZz9LHFRI899hgtW7Z0WR4wp96GkeVFxkYSHRdNgVwFWNBrgd1xUhUcDCVKQG0HD3iPHDlCw4YNXZrJHFEaRhb3+M+P8+D0B4mJi7E7Sqri4yEkxDqadGRYZ1xcHMeOHXPpGEowR5SGkeUN8BtA6MVQt3niJiU7d8LZs46fdp84cYLY2FiXP8JsCqVhZEGqyo7TO6h7d106VelkdxyHJV6fbJ1s14fbZcYYSjCn3oaRJb21+i0CpgSwMWyj3VHSJCgIataEe+5xbP3y5cvz8ccf4+fn59JcplAaRhbz2cbPeGv1WwyoM4CGpVx7k8OZrl+HNWscP+0GKFeuHC+++CIlSpRwXTAcOPUWkWLAEKB80vVV9QnXxTIMIz2mbp/KC7+9QLfq3Zjy0BS8bm/G5bbWroWoKMfaqiXavXs3vr6+lCtXznXBcOwa5S/AGiAEiHNpGsMw0u2nPT8x+NfBtKnYhpndZrpNX0lHzZxpPbZ4332Of+bpp58G4Pfff3dRKosjv8m8qjrKpSkMw8iQkMMh9PmxD41KNWJBrwXkypHL7khpcuwYBAbCsGHg6+v4544cOUKrVq1cFyyBI8flizIy6ZdhGK61MWwjXWd3pVrRaix+bDG+OdNQadzEJ59Yf76YhrnMoqKiOHHiRKbMbnDHI0oRicDqSC7AaBGJAmISXququk/zOsPIpq7HXKfbnG7cnf9ufuv3G3flucvuSGl29ixMmQL9+kEyM73c0fHjx1FVlw82hxQKparmd/neDcPIkDw+eZj16CzKFSpHyXwl7Y6TLhMnQmQkjErjBb7MGkMJDpx6i8hyR94zDCPzhF0OuzGt7P3l76d8ofL2Bkqny5dh0iR45BGoVi1tn/Xz82POnDkuH0MJKZ9658aaGbGoiNzFv5OFFQAcHA5qGIYrvLfmPWb+NZPWFVu7XXfytPjqK7h0CV57Le2fLVGiBD179nR+qGTInWZeEJGRwPNYRTHppGCXgSmqOsn18W4XEBCgW7ZssWPXhuE2omKj2Be+D78Srj+acpXISKhQweoSFBSU9s+vXr2aPHnyOK1zkIhsVdWA5Jbd8dRbVSeqagXgJVWtkOSrjl1F0jCys+sx13lu6XOEXwsnV45cHl0kAaZNg9On03c0CTBq1ChGjx7t1Ex3ktKpd7eEb08k+f4GVf3JZakMw7hJTFwMPeb1YMmBJbSu2JouVbuk/iE3FhsLH34IjRpBenvuHjlyhC5dMuf3kNKA84dSWKaAKZSGkQni4uN4/OfHWXxgMV92+tLjiyTAnDlw5Ah8+qljfSdvde3aNf75559MueMNKQ8PGpQpCQzDuCNV5dklzzJr1yzef/B9ng542u5IGRYfDxMmWF2CHkrpcCwFoaGhAJkyhhIc7EcpIp2AmsCNOdFU9W1XhTIMw/L6itf5autXjGo2ilebv2p3HKdYvBh27YIffgCvdPbsyMwxlODYOMqvgF7ACKwhQj0A17bqMAyDj9Z9xPtr3+ep+k/x/oPv2x3HKVThvfegfHno3Tv927nvvvtYv359poyhBMee9W6qqgOAC6r6FtAEKOPaWIaRvU3ZOoVXQl6hd63eTO44GUnPhTw3tHo1bNwIL78MOTLQ3Ch//vw0adIE37R00MgARwrl9YQ/r4nIPVjPe2fO8a5hZEOqyrJDy+hYuSPTu07H28vb7khO8/77ULw4DMrgHZB58+axaNEi54RygCM1fZGIFAI+ArZh3fH+1qWpDCObUlVEhDnd5xATF+MRE4I5autWa2D5hAlW38mMeP/99ylZsiSdO3d2TrhUpHpEqarvqOpFVf0R69pkNVV90/XRDCN7WXN0DU3+14TTV06TwysHeXwyWE3czIQJULAgPPNMxrd15MiRTLuRA47dzMkrIm+KyBRVjQKKi0jmlHHDyEai46KJ0zi8Jeucaifatw9+/BGGD4cCGWzQePHiRS5evJiphdKRU++pwFasmzgAYcA8IPMuEBhGFnY1+iq+OX15sOKDbKqwyaPmuXHUBx9ArlwwcmTGt5XZQ4PAsZs596rqh1g3cVDV6/zbScgwjAw4evEoNf6vBv/b9j+ALFkkjx+3xkwOHmzdyMmoxEKZWYPNwbEjymgRyYN1EwcRuReIcmkqw8gGzlw5Q5sf2nA56jIB9yTbtCZLSJzm4aWXnLO9rl27cvz4cYo7o+o6yJFCOQ5YBkDeHq4AACAASURBVJQRkRlAM2CgCzMZRpZ3MfIi7QLbcSLiBMH9g6lTso7dkVzi3Dlrmoe+fcFZM8p6eXlRunRp52zMQakWSlUNEpGtQGOsU+6RqnrO5ckMI4u6Gn2VTjM7sfvsbn7t8ytNyzS1O5LLfP45XL+e9mkeUjJp0iR8fX0ZlNHBmGngyF3vH4BuwCFVXWSKpGGkX1RsFN3mdmNj2EZmPjqTdpXa2R3JZS5fhi++gK5doXp15233yy+/ZOHChc7boAMcuXI8Fbgb+EJEDonIjwndzw3DSIO4+Dj6LehH0KEgpjw0he41utsdyaW+/houXkx/Y97kqCqhoaGZeiMHHDv1XiEiq4EGwAPA01idhCa6OJthZCkjl41k/u75fNL2E56o+4TdcVwqMhL++19o3RoaNHDeds+ePcu1a9cydWgQOFAoE2Zc9AU2AGuABqr6j6uDGUZW81jtxyiVvxT/afIfu6O43PffW9M8zJjh3O3aMYYSHDv13glEA7UAP6BWwnAhwzAcsO3UNgCalmnKay2ceB7qphKneWjYEB54wLnbPnPmDN7e3pl+6u3Is94vqOp9wCNAONY1y4uuDmYYWcGv+36l/jf1WbQ/+zzINncuHD5sXZt0dne4Ll26EBkZSc2aNZ274VQ4cur9LNACqA8cBb7DOgU3DCMV7Su1Z2L7ibSv1N7uKJlC1Wp+UaMGuGrerxwZaWSZ3n06sE4e4L/AVlWNdXEew8gSQg6H4FfCj+K+xXmu0XN2x8k0ixfDX3/B9Onpn+YhJa+//joFCxbklVdecf7GU+DIqfdHqrrJFEnDcMyKIyvoNLMTL/z2gt1RMpWq1Zi3XLmMTfOQknnz5rFlyxbXbDwFWe8JfMOw0R8n/qDLrC5UKVKFLzp8YXecTLVmDaxfb03z4OOCfsPx8fEcPXo00+94QwqFUkRyZWYQw/B0u/7ZRYcZHSiRrwRB/YIonKew3ZEy1XvvWd2BnnDRENGTJ08SHR3tXoUSa9xk4iOM6SIi7UVkn4gcFJFk59oUkZ4isltE/haRmendl2HY6fCFw7T9oS25vHMR0j+Eu/PfbXekTLVtG/z2G7zwQsanebiTzJ7LO6mUbubkFJHHgaYi0u3Whar6U0obFhFvYDLQBqvZ72YRWaiqu5OsUxl4DWimqhdEJPP6JhmGk5yMOEnr6a2Jiovi94G/U+Gu7Df33oQJVudyZ0zzcCdXr16lTJkyVKxY0XU7uYOUCuXTQF+gEPDQLcsUSLFQAg2Bg6p6GEBEZgMPA7uTrDMEmKyqFwDMEz+Gpwm/Fk7bH9py9tpZVgxYQc3imTu+zx3s3w/z58Orr1pz4rhKu3btOHbsmOt2kII7FkpVXQusFZEtqvq/dGy7FHA8yeswoNEt61QBEJF1gDcwTlWXpWNfhmGLL/74goPnD7K071IalHLiQ80e5K23rGkenn/e7iSu48g4yh9E5DngvoTXq4GvVDUmlc8lNyZfk9l/ZaAlUBpYIyK1VPWmJ39EZCgwFKBs2bIORDaMzPHmfW/ycNWHqXt3Xbuj2GLFCpg5E954wznTPKTkySefpGjRonzwwQeu3VEyHBke9H9YT+X8X8JXPeBLBz4XBpRJ8ro0cDKZdX5R1RhVPQLswyqcN1HVb1Q1QFUDihUr5sCuDcN1YuJieGHZCxy/dBxvL+9sWySjo2HYMKhYEUaPdv3+li9fzokTJ1y/o2Q4UigbqOrjqroi4WsQVsu11GwGKotIBRHJCfQGbu22+TNW6zZEpCjWqfhhx+MbRuY7cP4A3+34jpDDIXZHsdUnn1jT0E6a5Lo73YliYmI4fvy4LUODwLFT7zgRuVdVDwGISEUgLrUPqWpswnPiv2Fdf/xOVf8WkbeBLaq6MGFZWxHZnbDNl1U1PL0/jGFkhhrFarD/2f2UyFfC7ii2OXIE3nkHunWDDh1cv7+wsDDi4+PdulC+DKwUkcNY1x3LAQ5NVqGqS4Alt7w3Jsn3Cvwn4csw3NqYlWPInzM/Lzd7OVsXSbDm5/bygs8+y5z92dWHMpEjHc6XJ4x3rIpVKPeqqpmu1shW/rvhv7zz+zsMrjsYVUWc3T/Mg/zyC/z6K3z0EZQpk/r6zuDl5UXTpk2pVKlS5uzwFmId1HmOgIAAteOheCP7+t+2/zH418H0qNGDWY/OwtvL2+5Itrl61WqhVqCA9TSOK57ptouIbFXVZCdYz/zGbobhQebvns/QRUNpd287ArsFZusiCdZ1yWPHrAYYWalIpsZ0DzKMO/jt4G889uNjNCndhB97/khO75x2R7LV7t3Wne6BA6F588zdd8eOHRk6dGjm7jQJR+b1FhHpJyJjEl6XFZGGro9mGPZZd2wdj8x5hJrFa7LosUX45vS1O5KtVK0xk/nzW/PhZLYdO3YQF5fqYBuXcXTAeROgT8LrCKxmF4aRJe04vYNOMztRpmAZfuv3G4VyF7I7ku0CA2H1aqv5RWY/83H9+nVOnTplS9egRI4UykaqOhyIBEhoYJG9z0GMLC02PpaKd1UkuH8wxX1NQ6sLF+Cll6BRIxg8OPP3f/ToUcC+oUHg2M2cmISWaQogIsWAeJemMgwbXIm+Qr6c+ah/d322Dt2arYcAJfX663DuHCxb5pp5cFJj9xhKcOyI8nNgAVBcRN4F1gLvuTSVYWSy89fPE/BNAO+veR/AFMkEmzfDV1/BiBFQ16ZH2u+66y569epl2xhKcHAcpYhUAx7EGnC+XFX3uDrYnZhxlIYrxMbHMnLpSHrX6k2Lci3sjuMW4uKs0+2TJ2HvXmvsZFaW7nGUIuIF7FTVWsBeV4QzDDtdi7nGpchL3J3/biZ3Mvcok/rqK9i6FWbPtrdIRkVFkSuXvVN4pXjqrarxwJ8iYppAGllOdFw03ed2p8XUFkTGRtodx62cPm21TmvdGnr2tDdLs2bN6N69u60ZHLmZczfwt4j8AVxNfFNVu7gslWG4WFx8HP0X9GfpwaV80/kbcufIbXckt/LSSxAZCZMng92Xa0NDQ2nQwN7u8Y4UyrdcnsIwMpGq8sziZ5j791w+bP0hQ+oPsTuSW1m5EmbMgDffhCpV7M0SERFBeHi4rWMowbHuQaszI4hhZAZVZVTIKKZsm8Lo5qN5udnLdkdyK0m7lr/2mt1p3GNoEDj2CGNjEdksIldEJFpE4kTkcmaEMwxnm7B2Ah+t/4hhAcMY32q83XHcziefWHe4M6NruSM8plACk7AeXzwA5AEGJ7xnGB7ly81fMnrFaB6r/RhfdPzCjJW8RWZ3LXdEhQoVeOmll2wdQwkOjKNMmK42QER2qqpfwnvrVbVppiS8hRlHaaTH+uPraf5dczpV6cRPPX/Cxzsb9QhzUJcu1qyKe/ZkXkNed5LRfpTXEiYH2yEiHwKngOzdSsXwOI1LN+bzDp/zZN0nTZFMhh1dyx1x/PhxChcujK+vvSXHkSPKcsAZrEYYLwAFgf9T1YOuj3c7c0RppMXaY2splb8UFe6y9xqXO3PnruV+fn6UL1+ehQtvncDV+TJ0RKmqRxO+jcQMFTI8SHRcNP1+6kflIpUJ7h9sdxy35a5dy1WV0NBQWrZsaXeU1AuliDQDxmHNvnhjfVWt6LpYhpFxOb1zsuixRdyV+y67o7itxK7lgwZlftfy1Jw/f56IiAjb73iDY9co/4d1yr0VB+bzNgy7HblwhJ/3/szzjZ+nVvFadsdxW0m7ln/wgd1pbucuQ4PAsUJ5SVWXujyJYTjBqYhTtPmhDeHXw+lVqxf35L/H7khuK7Fr+TffZH7XckeEhoYC2P5UDqRQKEWkXsK3K0XkI+An4MZ83qq6zcXZDCNNzl8/T9vAtpy+cpqQASGmSKbg6FF47jlo0gSefNLuNMnz9/dn4sSJto+hhBTueovIyhQ+p6rayjWRUmbuehvJuRJ9hdbTW7P99HYWP7aY1hVb2x3JbcXEwP33w99/w/bt1uOKRjrveqvqA66LZBjOExUbRdfZXdl8cjPze8w3RTIVY8fChg1Wn0l3LpJbt26lSJEibnHq7ciz3iNFpEDCtLXfisg2EWmbGeEMIzWx8bH0+bEPy48s57su3/FI9UfsjuTWgoOtmRQHD4ZevexOk7K+ffvy4osv2h0DcOxZ7ydU9TLQFigODAImuDSVYTggXuMZvHAwC/Yu4LN2n/G4/+N2R3JrZ85A//5QvTpMnGh3mpSpKkePHnWLO97g2F3vxM4BHYGpqvqnmG4Chhs4GXGSZQeXMe7+cYxsPNLuOG4tPh4GDIBLlyAkBPLmtTtRyk6fPk1kZKRHFcqtIhIEVABeE5H8mOlqDTdQukBpdj6zk2J53XBsi5v5+GMICoKvv4ZaHjC01J3GUIJjhfJJwB84rKrXRKQI1um3Ydhi4saJhF0O48M2H1Lct7jdcdzexo3W3Nw9esAQD2nm7nGFMmGCsW1JXocD4a4MZRgpOXD+AKeunCJe4/EWb7vjuLWLF6FPHyhd2hpY7ikXzVq1asVPP/1ERTe5Le/QvN7uxIyjzL5i4mLw8fZBVYmNjzXt0lKhah1F/vKL1fCicWO7E7m3lMZROnLX2zBsF3I4hOqTq3Mg/AAiYoqkA775Bn78Ed591/OK5OLFi9m8ebPdMW5wqFCKSHMRGZTwfTERcY8LB0a2sOH4BrrO7kpen7wUzVvU7jge4a+/4PnnoV07a+pZTzNixAg+/fRTu2Pc4MiA87HAKCBxTjYfINCVoQwj0c4zO+k4syMl85UkqH8Qd+UxLdNSc/WqNZi8UCGYPh28POy8MTY2lmPHjrnNjRxw7K73I0BdEm7oqOrJhCFChuFSB88fpO0PbfH18SVkQAgl85W0O5JHeP55aybFoCAo7oGDAk6cOEFcXJxbPLqYyJFCGa2qKiIKICJmvhzD5U5cPkHr6a2JjY9l5eMrKV+ovN2RPMLs2fDttzB6NLT20Efe3W1oEDh2jXKuiHwNFBKRIUAIMMW1sYzs7Ny1c7T5oQ3nr59nWb9lVC9W3e5IHuHQIRg6FJo2hXHj7E6Tfu5YKB0ZR/mxiLQBLgNVgTGqaiYgMVziWsw1OszowJGLR1jWdxkB9yQ7WsO4RXQ09O4N3t4wc6Z7zX2TVt27d6d27dqUK1fO7ig3ODJnzgvAPFMcjcyQJ0ce2lZsy9j7x3J/+fvtjuMxRo+GLVus4UBuVF/SJX/+/AQEuNc/kI5coywA/CYi54HZwHxVPePaWEZ2ExMXw8mIk5QrVI53H3zX7jgeZckSa4Kw4cOhWze702Tct99+S6lSpejQoYPdUW5I9Rqlqr6lqjWB4cA9wGoRCXF5MiNbeeG3F2j4bUPOXz9vdxSPcvIkPP44+PlZjS+ygnHjxjFnzhy7Y9zEkSPKRP8Ap7Ge8/bAQQeGOxvRcAQ1itWgcJ7CdkfxGHFx0K8fXLsGc+ZA7tx2J8q4qKgoTp486VY3csCxAefPiMgqYDlQFBiiqn6ObFxE2ovIPhE5KCKvprBedxFREXGvCxOGywUdCkJVqVq0KsMaDLM7jkd5/31YuRImT4Zq1exO4xzHjh1DVd1qDCU4NjyoHPC8qtZU1bGqutuRDYuINzAZ6ADUAPqISI1k1ssPPAdscjy2kRV8tO4j2gW2Y/au2XZH8Thr1lhz3/Tta516ZxXuODQIUiiUIlIg4dsPgWMiUjjplwPbbggcVNXDqhqNdSPo4WTWeydhH5FpzG54sClbp/BKyCv0qtmLnjV72h3Ho4SHw2OPWRODffml57ROc0TiXN4eUyiBmQl/bgW2JPy5Ncnr1JQCjid5HZbw3g0iUhcoo6qLHA1seL45u+bw1KKn6FCpA9MfmY63l+kp6ai4OOsI8swZ6ymc/FnsYeIhQ4Zw+vRpSpUqlfrKmSil6Wo7J/yZ3tKe3L9zN5pfiogX8CkwMNUNiQwFhgKULVs2nXEMd7D0wFL6LehH87LNmd9zPjm9c9odyaO89BIsXmxdl6xf3+40zicilChRwu4Yt3HkZs5yR95LRhhQJsnr0sDJJK/zA7WAVSISCjQGFiZ3Q0dVv1HVAFUNKFbMzI/iqdYcXcOjcx/Fr4Qfv/b5lbw+bj7DlZuZNAk++wxGjoRhWfS+1zvvvMOMGTPsjnGblK5R5k64FllURO5Kcn2yPNZ4ytRsBiqLSAURyQn0BhYmLlTVS6paVFXLq2p5YCPQRVVN+/IsaNupbXSe1ZlyhcqxrO8yCuYuaHckj7J4sVUgu3SxBpdnVV988QWrV6+2O8ZtUhpH+RTwPFZR3Mq/p9KXse5mp0hVY0XkWeA3wBv4TlX/FpG3gS2qujDlLRhZhaoyeOFgCuUuRFC/IIr5mrOCtNi+3eov6e9vPcftnUUv6V65coWzZ8+63Y0cSPka5URgooiMUNUv0rNxVV0CLLnlvTF3WLdlevZhuD8R4adePxETF0OZgmVS/4BxQ1gYdO4MhQvDokXgm4WbHCbe8Xa3MZTgWPegL0SkFtZYyNxJ3p/uymCG5ztz5QzfbvuW11q8ZvpJpkNEhFUkIyJg3Tq4+267E7mWuw4NAsengvgi4esBrDGPXVycy8gCZu+azXtr32N/+H67o3ic2FjrdHvXLpg3D2rXtjuR64WHh5MzZ063LJSpTlcrIn8BdYDtqlpHREoA36rqQ5kR8FZmulrPoaocuXiEine5x9zMnkLV6gT05Zfw9ddWM97sIj4+HhFBbBhFn9Hpaq+rajwQm/C0zj+A+T/fSFZ0XDSP//w4f535CxExRTIdPv3UKpKvvJJ9imR4eDiqipeXly1FMjWOFMotIlIIa/qHrViTjP3h0lSGR4qLj6PvT32Z/ud0tp3aZnccj7RggTWovHt3q+lFdhAXF0eLFi0YMmSI3VHuyJF+lMNU9aKqfgW0AR5X1UGuj2Z4ElXlqUVPMX/3fD5p+wmP+2ehTg2Z5I8/rCYXDRt65jSz6TVv3jz27NlDmzZt7I5yR3e8Riki9VL6oKracshgrlG6H1Xl5eCX+WTDJ7zR4g3eafWO3ZE8TmgoNGoEefPCpk2eOc1sesTFxVG7dm1EhL/++gsvG/91SOkaZUrDg1Ia/69AqwylMrKM99a8xycbPuHZBs/y9gNv2x3H41y8CJ06WROErVqVfYok/Hs0OWfOHFuLZGpSGnD+QGYGMTzTpD8m8cbKN+jv15+JHSa65YV4dxYTAz16wP79EBQE1bPZzLxfffUVNWvWpHv37nZHSZEjszAOSO59M+DcCNwZyIilI3i46sN89/B3eIn7HhG4I1V45hkICYGpU+GBbHhosnjxYo4fP+7WR5Pg2Jw5DZJ8nxt4EOvOtymU2VzoxVBaVWjF7O6zyeGVlumXDIAJE+B//4M33oCBA+1Ok7ni4uJQVXx9fanmAfNYpDrg/LYPiBQEflBVW57OMTdz7BcdF32jj2RMXAw+3j42J/I8c+ZA795Wp/LAwKzVpdwRM2fOZOzYsSxfvtxtesxmdMD5ra4BlTMWyfBUm09spsoXVW6MkzRFMu3Wr7e6lDdvDt99l/2KZFxcHG+//Ta5c+emdOnSdsdxiCPXKH/l387kXljNMea6MpThvorkLUKVIlUoma+k3VE80qFD8PDDULYs/Pwz5Mpld6LMN2fOHPbt28e8efPc/tpkIkee9b4/yctY4Kiqhrk0VQrMqbc9/rn6D0XzFjU3bDLg/Hlo0gTOnYONG6FyNjwvi4uLo2bNmuTMmZMdO3a4VaFM7zhKAFR1dcJGCiSuLyKFVfW8U1MabutkxEmaf9ech6s+zKftP7U7jke6dg26drUGli9fnj2LJMDChQvZt28f8+fPd6simRpHTr2HYk0pex2Ix+p0rpjGGNlC+LVw2v7QlrPXztKndh+743ika9esvpLr1lkdyps3tzuRfR5++GEWLlxIp06d7I6SJo6M6XgZqKmq51wdxnAvEVERdJzZkYPnD7K071IalmpodySPk1gkV6+G77+3ekxmV4ndgR56yJYOjRniyLHvIaw73UY2EhkbycOzH2brya3M7TGXBypkw9HQGXTtGjz00L9Fsl8/uxPZJzY2lqZNmzJ16lS7o6SLI4XyNWC9iHwtIp8nfrk6mGGf2PhYes/vzcrQlUzrOo0uVU1D+7RKLJIrV5oiCTB79mw2btxIwYKeOfumI6feXwMrgL+wrlEaWVi8xvPEL0/wy75fmNRhEv38svnf8HS4ds2aVnblSqtdWnYvkrGxsbz99tvUqVOHrl272h0nXRwplLGq+h+XJzHcwvjfx/PDzh8Y/8B4hjccbnccj5NYJFesMEUy0axZszhw4AA//fSTR93pTsqRcZTvAkeBX4GoxPftGh5kxlG61onLJ5i9azb/afIf0wkojZIWye+/h/797U5kv/j4eKpVq0bevHnZtm2bWxfKlMZROlIojyTztqqqLcODTKF0jeBDwbSq0ApvL2+7o3ika9esJ26WLzdF8lbr168nPj6e5m4+LiqjA87db+5Iw6m2ndpG28C2fNbuM0Y2Hml3HI9z/fq/RXLaNFMkb9W0aVO7I2SY6UdpUO/uesztPtfc3U6H69et0+3EIjkg2b8t2dPMmTNZs2YNH3/8Mb6+vnbHyRDTjzIbCz4UTOE8hal/T3161OxhdxyPk7RITp1qimRSsbGxjB07lnz58pE3b16742SYI6feI5K+TuxH6bJERqZYf3w9Xed0JeCeAFY9vsrcuEmjpKfbU6dabdOMf82YMYODBw/y888/Z4n/t0w/ymzoz9N/0nFGR0rlL8Xc7nOzxP/ImSmxSIaEWP0kTZG8WWxsLO+88w7+/v506ZI1LueYfpTZzP7w/bQNbEuBXAUIGRBCiXwl7I7kUW4tktltCgdHzJgxg0OHDmWZo0lw7Brlx0m+t70fpZF+xy8dp80PbVBVgvsHU7age7Tg9xTXr1ut0kyRTFmLFi0YM2ZMljmahBTGUYpIJaCEqq675f0WwElVPZQJ+W5jxlGmz9mrZ2kxtQWnrpxi1eOrqHt3XbsjeZTEIhkcbE0INmiQ3YkMZ0vvnDmfARHJvH89YZnhIS5FXqJdYDuOXTrGoj6LTJFMo8hIeOQRUyRTExMTw5AhQ/jzzz/tjuJ0KRXK8qq689Y3VXULUN5liQynW3JgCbv+2cWPPX+kRbkWdsfxKFevWtckg4Lg229NkUxJYGAg3377LceOHbM7itOldOp9UFUrpXWZq5lT7/Q5cuEIFe4yD1mlxalTVqu07dtNkUxNTEwM1apVo1ChQmzZssUjb+Kk99R7s4gMSWZjTwJbnRXOcI24+DiG/jqUtcfWApgimUZ//w2NG8PevfDLL6ZIpiYwMJDDhw8zbtw4jyySqUnprvfzwAIR6cu/hTEAyAk84upgRsacv36e1UdXU61oNZqXde9mBO5m+XLo1g18feH336FePbsTubeYmBjGjx9P/fr16dy5s91xXOKOhVJVzwBNReQBoFbC24tVdUWmJDPSLV7jKeZbjG1Dt+Gb07Ofsc1s06bBkCFQrRosXmzNv22kLDY2lv79+9OkSZMseTQJDrRZczfmGmXKJqydwJ9n/mR61+n4ePvYHcdjqMLYsfDOO9C6NcyfDx46a0Gmio+Pd+sek2mR3muUhof5cvOXvLb8NbzEy/SVTIOoKKuhxTvvwBNPwJIlpkg64tSpU9SvX5+1a9faHcXlHHkyx/AAs/6axfAlw+lcpTPTHp6Gl5h/Ax1x4YI1RnL1ahg/HkaPhix69uhUkZGRPPLII+zfv598+fLZHcflTKHMAhbtX8SAnwdwX7n7mNt9rjnldtCRI9CxIxw+DDNmwGOP2Z3IM6gqQ4YMYdOmTfz000/4+/vbHcnlTKH0cKtDV9NjXg/8S/qzsM9C8vjksTuSR9i0yeolGRNjPXFz3312J/IcH3zwAYGBgYwfP55HHskeA2DM+ZkH23JyCw/NeogKhSqwtO9SCuQqYHckj7BgAbRsaQ3/2bDBFMm0UFW2b99Onz59GD16tN1xMo05ovRQe8/tpX1ge4rkLUJw/2CK5i1qdyS3pwqffQYvvggNG8LChVC8uN2pPIuIMGvWLGJiYrLsUKDkmCNKD3VX7rsIuCeA4P7BlCpQyu44bi8uDp57Dv7zH+vmzcqVpkimxdmzZ3nooYcIDQ3Fy8uLXLly2R0pU7m0UIpIexHZJyIHReTVZJb/R0R2i8hOEVkuIuVcmScrCL8WTkxcDCXylWBZv2VUKmzLI/ce5epVqzhOmmQdTc6bB3nMpVyHRUVF0a1bN0JCQjh37pzdcWzhskIpIt7AZKADVlf0PiJS45bVtgMBquoHzAc+dFWerCAqNopW01vx+M9m7gFHnToF999vPWUzaRJ8/DFkkfHRmUJVeeaZZ1i7di3Tpk0jICDZ8dhZniuvUTYEDqrqYQARmQ08DOxOXEFVVyZZfyPQz4V5PF6uHLl4uv7TVC5ipixyxN9/W8N/zp2zGltk0ceQXeqzzz5j6tSpjBkzhl69etkdxzau/Le1FHA8yeuwhPfu5ElgaXILRGSoiGwRkS1nz551YkTPEBUbxc4zVmvQZxo8Q+uKrW1O5P5+/hmaNoXoaKuxhSmSaRcTE0NgYCCPPvooY8eOtTuOrVxZKJO7JZbsg+Ui0g+rM9FHyS1X1W9UNUBVA4oVK+bEiO4vNj6WPj/2odl3zTh95bTdcdxedPS/N2wqV7bGS9avb3cqz+Tj48Pq1av5/vvvs8zz3Onlyp8+DCiT5HVp4OStK4lIa+B1oIuqRrkwj8eJ13iG/DqEBXsXMP6B8ZTMV9LuSG4tNBRatIBPP4Vnn4V160z3n/QIDw9nxIgRXLlyhXz58uHrazpQubJQbgYqi0gFEckJ9AYWJl1BROoCX2MVyX9cmMXjqCr/+e0/TNsxjbH3j2Vk45F2R3JrCxdC3bpWo9158+CLLyCbjWBxipiYGHr06ME333zD3r17KsdmCgAAGQ9JREFU7Y7jNlxWKFU1FngW+A3YA8xV1b9F5G0RSZzH8iMgHzBPRHaIyMI7bC7beXv120zcNJGRjUYy9v7sfX0oJTEx8NJL1rw2FSvCtm3QvbvdqTzXyJEjWblyJd9++222vcOdLFX1qK/69etrVvfZhs+UcejAnwdqXHyc3XHc1tGjqo0bq4LqsGGq16/bncizTZ48WQEdNWqU3VFsAWzRO9Qd07jXzUzbMY1BvwzikWqPMLfHXHJ4madMk7NoETz+uHVEOWUKZOORK05x9epVKleuTEBAAAsWLMDbO/v1M02pca/5W+hGQi+GMuTXIbSu2JpZj84yRTIZMTHw+uvw0Ufg7w9z51p3t42M8fX1ZcOGDRQuXDhbFsnUZO97/m6mfKHy/NTzJxb0WkCuHOZOxK2OH7e6/nz0ETz9tNX5xxTJjLl48SITJ04kPj6ecuXKkT9/frsjuSVTKN3AxrCNrDhizdn2UNWHyJcz63eMTqulS6272jt3wsyZ8OWXkDu33ak8W2xsLL169eLll1/m77//tjuOWzOF0maqyqiQUYxcNpK4+Di747id2Fh47TXrUcRSpWDLFujTx+5Uni86OpohQ4YQFBTEl19+Se3ate2O5NbMRTCbiQgLei0gIirCTAh2ixMnoHdvWLsWhg61ekmarj8Zd/LkSbp3786GDRsYM2YMTz75pN2R3J45orTJicsnGLFkBFGxURTOU5hyhUyHuaR++826WbN9OwQGwtdfmyLpLIcPH2bPnj3MmTOHt956y+44HsEUShucu3aOtoFt+f7P7zl4/qDdcdxKbCy88Qa0bw8lS1qn2n372p3K86kqmzZtAqB58+YcPXqUnj172pzKc5hCmckuR12mw4wOHDp/iIV9FlKzeE27I7mNbdusKRrefReefNJqaFGtmt2pPN/169cZNGgQjRs3Zs2aNQAUKGDmV0oLUygz0fWY6zw8+2F2nN7B/J7zaVm+pd2R3MK1a/Dyy9CggdVod/58+PZbyJvX7mSeLzQ0lGbNmvH9998zbtw4mjVrZnckj2Ru5mSSmLgYes3vxerQ1QR2C6RzFdMgEWD5cutGzeHDMGQIfPghFCpkd6qsISQkhN69exMbG8uvv/5KZ9OUM93MEWUmiNd4Bv4ykF/3/8rkjpN5rPZjdkey3fnz8MQT0Lq1NTXDihXwzTemSDrT4cOHKVmyJJs3bzZFMoNMoXQxVWXEkhHM/Gsm77V6j2caPGN3JFupWo8dVq8O06fDq69ag8gfeMDuZFlDRETEjeuQQ4cOZcuWLVQ2jy9lmDn1drGrMVf54+QfvNL0FV5tfttElNlKWBj8f3tnHl9Vde3x78rAGIghIYLWEGkYVDRAGC3aRAQqnwKiCFinNOJAlNantJW2jyp5Vn22WrG0joTJ54ChDBoMg2FUIERQk8jYEgXCFCpjCCR3vT/2SbzQhHuT3ORm2N/P53xyzj075/zWueesu8/ee62dnAxLlpis42VDgCy+YceOHYwePZq9e/eyZ88ewsLCaGHDl3yCdZS1iEtdhDQLYXXialoGtWxSE8a743KZcZC/+Y0Z/vOnP8EvfwlB9u7zGUuWLOHuu+8mODiYBQsWEBYW5m9JjQr76l1LvJH9BsPfHs7pc6dpFdyqyTrJr7+GG280Ncn+/SEnx8ytbZ2kb1BV/vCHPzBy5Ei6dOlCdnY2gwcP9resRod1lLVEUEAQzQKbNdlUaWfPQkqKebXOy4PUVFi2zGQht/gOEaGwsJDExETWrl1Lp042wqs2sIl7fcyJ4hO0aW5SValqk6xJbtgAEyaYebXHjYOXX4ZLL/W3qsbFV199RWlpKT179qS0tJSAgIAmea/5kosl7rU1Sh+yNn8t0S9Hl6dMa2o37smT8NhjZj7tY8fMhF/vvmudpC9RVebMmcOAAQNITk5GVQkMDGxy91pdYx2lj9hSsIWfvvNT2rdqT4/IHv6WU6eompkPr7nG1B6Tk01tcsQIfytrXKxbt46BAwdy33330atXL9LS0qyDrCOso/QB249sZ9i8YVzS4hKW37OcyNaR/pZUZ6xcaeKzx46FNm1MSrS//hVsKLFvWbFiBTfccAPffvstqamprF69mo4dO/pbVpPBOsoa8s2xbxgydwgiwvJ7lnNF6BX+llQnfP45DB1qImsOHjSdNV98ATaU2HccPnyYVatWAZCQkMCMGTPYuXMniYmJdl6bOsY6yhpw6NQhhswdwvHi4yy7exldw7v6W1Kts2uXSaYbFwfZ2fDnP8OOHZCYCPbZ9Q1FRUU899xzxMTEMG7cOIqLiwkMDCQ5OZlWNlOIX7COspp8d+Y7hs0bxrfHvuWjn31EbIdYf0uqVQ4cMG2PV11lImt+9zuTyOLxx+3cNb7C5XIxd+5cunXrxpQpU4iPj2f16tU0b24nmvM3TXOQXw1xqYtR744i91AuS+5cwo+iGu/75rFjZtbDl16C4mKT6ee//xts85jv2bhxI/feey9xcXHMmTOH+Ph4f0uyONhxlNXk/dz3CZAAxlw9xt9SaoUzZ+Bvf4M//hEKC814yJQUOz2sr8nLy2PDhg0kJSUB8MknnxAfH09AgH3Zq2suNo4SVW1QS1xcnPqLktISzdqX5bfz1wUlJaqpqapRUaqgOmSI6ubN/lbV+CgoKNAHH3xQAwICtH379nry5El/S2ryAJu1Er9jf7aqwPPrn2fgWwPZUbjD31J8jqoZIB4bCz//ObRvD8uXm7DDuDh/q2s8nDp1imnTphETE8PMmTN59NFHycvLo3Xr1v6WZrkIto2yCjzS9xE6hHRodL3b69aZvJDr10NMDLz3HowZYxLqWnzL/v37SUlJYdSoUTz77LM2V2QDwTpKL3gv5z1GdBtBaItQknol+VuOTygthQ8/hOnTTXbxDh3g7383k3oFB/tbXeNAVcnKyuKtt97iyJEjpKWl0aVLF3bu3El0dLS/5VmqgK0zeGDGphmMTxvPyxte9rcUn/Ddd/Dii6ZT5tZbzRjI55834yMfftg6SV9w9OhRpk+fTmxsLP3792fevHmEhoZSWloKYJ1kA8TWKC/CvC/n8ejSRxnZbSSTr5/sbzk1Yts2eOUVmD0bTp2CQYOMgxw92uaG9AUulwuXy0VQUBCpqalMnjyZvn378uqrrzJ+/HhCQ0P9LdFSEyrr5amvS131ei/atkgDnw7UhFkJWnSuqE7O6WtKS1U/+kh16FDTg92smWpiomp2tr+VNR727dunzzzzjHbu3Flnz56tqqqFhYW6detWPyuzVBUu0utt6xIVkPmvTMbOH0vcZXEsGr+IFkENK/Tk+HGYNcvUIHftMoPDU1LMYPHIppOvo9ZQVT788EPeeOMN0tPTKS0tJT4+vjxJRbt27WjXrp2fVVp8iXWUF7Bp3yZGvjuSmHYxpP8svTwJb0Ng506TuSc1FU6cgAEDjIO87TZo1szf6ho+hYWFhIeHAzB16lQOHDjAr3/9a5KSkoiJifGzOkttYh2lG7mHcrnl7VuIbB3JsnuWEd4q3N+SPKJqxjtOnw7p6aa9cdw4mDTJpD+z1IwjR46wdOlSUlNTycrKYt++fbRt25aFCxdy+eWXE2QbeJsE9lt24+nVT9M8sDnL71nOZW0u87eci3LiBMybZxzktm0mi/jUqabnukMHf6tr+KxZs4bJkyezefNmVJUrr7ySJ598EnVCfu3cNE0L6yjdSB2Vyv4T++kcVj9nwDp61GTuSUszETPFxdCnD8ydC3fcATbJTPXYs2cPGRkZZGRkMGHCBIYPH07btm0JDg7mqaeeYtiwYfTt29fGXzdhmryjPFp0lN+u/C0vDHmBNs3b0CW8fkVKHDgAixYZ55iZaebFjoqCiRNNXsh+/cDOBlB1ioqKePLJJ8nIyGD79u0AREVFcdtttwHQs2dP1q9f70+JlnpEk3eUG/du5O2v3ua+2PsYeMVAf8sB4JtvYMECs6xbZ9ohu3SByZPh9ttN7LV1jt6jquTk5JCRkYGI8MQTT9CiRQs+/vhjOnfuzMMPP8ywYcPo3r27nYPGUiE2zRpw5PQRIlpF+PSYVWXnTlNrXLAAsrLMZ9deaxzj7bebibvsM1w1li9fzjvvvENGRgb79+8HYPDgwaxYsQIwg8Tt67SljIulWWuSNcoSVwl3LbiLsVeP5farb/eLk1SFnBzjHNPSzDpA377w3HNmSI/Nl+CZkpIStm3bxpYtW/j888/Jyclh6dKlBAUF8Y9//IOFCxdy8803M2zYMIYOHcoVV3w/p5F1khZvaXKO0qUukhYl8X7u+9wYdWOdnvvUKTPPTHq6cY67dpla4qBB8Je/mHDCqKg6ldSgOHPmDDk5OVx11VW0bt2aN998k0mTJnHmzBkAWrZsSWxsLIWFhVx66aU888wzvPLKK3YiLkuNaVKOUlV57OPHmPvlXFISUnik3yO1dq6SEjO39aZN3y85OeBymUm4brrJtDmOGmWH81RGQUEB8+fPL68t5uXlUVJSwooVKxg8eDDXXHMNycnJ9OrVi969e9O1a9fzxjWGhYX5Ub2lMdGk2iinZk4lZU0KTwx8gheGvOCzhntVyM8/3ylmZ8Pp02Z/WJjpnS5brr8ebIQbnD17li1btpCfn8+ePXvKlwceeIDRo0eTnZ1Nnz59iIyMpHfv3vTu3ZtevXqRkJBQHiFjsfgKv7VRishPgJeBQOBNVX3ugv3NgTlAHFAIjFPVPbWh5aXPXiJlTQr397q/xk7y6FHT4eLuGA8dMvuaN4devWDCBOjf3zjGH/6waXbEuFwuNm3aVO4AyxziiBEjSE5O5tixYwwYMKC8fLt27YiOjqaoqAiA6667jn379tGxY0fbG23xK7XmKEUkEJgBDAH2AlkislhV89yK3Q/8W1VjRGQ88DwwztdaZm6ZyePLHmfM1WN47aeveXzoTp82ju/QITh82Pw9eNC8Om/caNoWjY3QvTvccsv3TvHaaxtvXLWqcvz4cYqLi4l0smvMnz+fffv2UVhYyMGDB8nPz6dfv36kpKQgIiQkJJS3IYaHhxMdHY3L5QIgIiKCxYsXEx0dTadOnWjbtu155wsODuayy+p3hJSlaVCbNcp+wC5V/SeAiLwLjALcHeUo4Cln/QPgryIi6sP2gLS8NB5Y8gA/vnwov+o8j5UrAs9zgO5L2WenTlV8rI4djUNMSjJ/4+KgIaUZVFWKi4s5e/ZsuVPavXs3Bw4coLCwkMLCQo4cOUJISAgTJ04EIDExkaysrPL9JSUl3HTTTaxcuRKAKVOmsHv3bgICAoiIiKBTp040c34pRIT09HQiIyPp1KkTISEh5+kREUaMGFGHV8BiqR616SgvB751294L9K+sjKqWiMgxIBw44gsB587BPXcF4LountX/s4D+586P8QsOhrCwc8AhQkKKCA09Q1TUWcLCzjJoUDe6dQtD5AiFhV/ToUMgZc+5qhIbG0tISAgFBQXs3Lmz/PMyH9+vXz9atWpFfn4+u3btOi+3ncvlIiEhgebNm5Obm0tubi4lJSXnLUlJSQQFBbFq1SqysrIoLS0t3+dyuZg2bRoAs2fPJjMzk9OnT1NUVMTp06dp1qwZS5cuBeChhx4iLS2NoqIiioqKUFWioqLIz88HIDk5mWXLlp13XXr06FHuKENDQ+nevTvh4eFEREQQHh5+3jwvmZmZtG7dmksuuaTC4TYJCQk1/BYtFv9Ta505InIHMExVJzjb9wD9VHWSW5lcp8xeZ3u3U6bwgmM9CDzobHYDtldRTgQ+cr5+prHYAdaW+kpTtqWTqravaEdt1ij3Ale4bf8A2F9Jmb0iEgSEAkcvPJCqvg68Xl0hIrK5st6shkRjsQOsLfUVa0vF1GZoQhbQRUSuFJFmwHhg8QVlFgP3OetjgE982T5psVgsvqDWapROm+OjQAZmeNBMVc0VkWmYuSkWA28Bc0VkF6YmOb629FgsFkt1qdVxlKqaDqRf8NlUt/UzwB21qcGh2q/t9YzGYgdYW+or1pYKaHCRORaLxVLX2PQpFovF4oFG4yhF5Ccisl1EdonIkxXsby4i7zn7N4pIdN2r9A4vbHlcRPJE5EsRWSki9XYCF0+2uJUbIyIqIvW2x9UbW0RkrPPd5IrI/9W1Rm/x4h6LEpFMEdni3GfD/aHTEyIyU0QOiUhOJftFRKY7dn4pIr2rdaLKJvxuSAums2g30BloBnwBXH1BmWTgVWd9PPCev3XXwJYEoJWzPrEh2+KUawOsATYAffytuwbfSxdgCxDmbEf6W3cNbHkdmOisXw3s8bfuSmy5EegN5FSyfziwFBBgALCxOudpLDXK8nBJVT0LlIVLujMKmO2sfwAMlvqZacGjLaqaqapObiI2YMao1ke8+V4AUoD/Bc7Upbgq4o0tDwAzVPXfAKp6qI41eos3tihQFnwfyn+Oga4XqOoaKhh77cYoYI4aNgCXiEjHqp6nsTjKisIlL6+sjKqWAGXhkvUNb2xx537ML2Z9xKMtItILuEJVP6xLYdXAm++lK9BVRNaLyAYne1Z9xBtbngLuFpG9mJErk2iYVPV5qpDGkri3oprhhd353pSpD3itU0TuBvoAP65VRdXnoraISADwEpBYV4JqgDffSxDm9TseU8tfKyI9VPW7WtZWVbyx5U5glqr+WUQGYsY791BVV+3L8yk+ee4bS42yKuGSXCxcsh7gjS2IyM3A74CRqlpcR9qqiidb2gA9gFUisgfThrS4nnboeHuPLVLVc6r6L0xOgvo485E3ttwPvA+gqp8BLTCx0w0Nr54nTzQWR9mYwiU92uK8rr6GcZL1tR0MPNiiqsdUNUJVo1U1GtPeOlJVfTvNpm/w5h5biOloQ0QiMK/i/6xTld7hjS3fAIMBROQqjKM8XKcqfcNi4F6n93sAcExVC6p8FH/3Wvmw92s4sAPTm/c757NpmAcPzBc9H9gFbAI6+1tzDWxZARwEtjrLYn9rrq4tF5RdRT3t9fbyexHgRUzO1a+A8f7WXANbrgbWY3rEtwJD/a25EjveAQqAc5ja4/3Aw8DDbt/JDMfOr6p7f9nIHIvFYvFAY3n1tlgsllrDOkqLxWLxgHWUFovF4gHrKC0Wi8UD1lFaLBaLB6yjbOKISKmIbBWRHBGZLyKt/KTjMfdzi0i6iFzirJ/00TnuEJGvRSTTB8dKFJFanXRcROJF5Hq37VkiMqY2z2mpGOsoLUWq2lNVewBnMWPQvEJEAn2o4zGg3FGq6nD1fejf/UCyqvpiDt1EoNYcpRM9Fg9c76GopQ6wjtLizlogBkwcuYhscmqbr5U5RRE5KSLTRGQjMFBE+orIpyLyhVO+jYgEisgLIpLl5AB8yPnfeBFZJSIfiMg2EXnbiZj4BcbpZJbV9kRkjxPdch4i8iu34z5dkREicqeIfOXUkp93PpsKDAJeFZEXLihfoS5nX5yIrBaRbBHJEJGOTq2uD/C2c31auh0rUkSynfVYMTk2o5zt3SLSSkQ6ickjWpZPtGz/LBF50bkG72F+tP7LOccNziludK73P23tsg7x98h6u/h3AU46f4OARZj8llcBS4BgZ9/fgHuddQXGOuvNMCF6fZ3tts5xHgR+73zWHNgMXImpIR3DxNsGAJ8Bg5xye4AIN13l224ah2LyJIrz/x8CN15gz2WY8Lv2jpZPgFudfauoIDKjMl1AMPAp0N4pNw4zSV6lx3L25TrX4lFMuOBdQCfgM2f/EuA+Zz0JWOisz3JsCnS2nwImux13Fia6LAATObPL3/dPU1kaS/YgS/VpKSJbnfW1mJkxHwTigCynYtUSKIspLwXSnPVuQIGqZgGo6nEAERkKXOdW4wnFJIc4C2xS1b1Oua1ANLDOS61DnWWLsx3iHHeNW5m+wCpVPeyc421McteFHo5dka7vMEk7ljvXIRATLueJT4EfOef9I/ATjHNf6+wfCNzmrM/F5OIsY76qll7k2AvVZPDJE5FLvdBi8QHWUVqKVLWn+wfOa+dsVZ1SQfkzbg+yUHHKKgEmqWrGBceNB9wzHZVStXtQgGdV9TUPZapDRboEyFXVgVU81lrgBkwtchHwG8x1qiznpvs1PFUFnfUx8XSjxLZRWipiJTBGRCIBRKSdVDwvzzbgMhHp65Rr43RCZAATRSTY+byriLT2cM4TmLRrFyMDSBKREOe4l5dpdGMj8GMRiXDaVe8EVns4bmVsB9qLyceIiASLyDVe6F0D3A3sdGp/RzFJKNY7+z/l+zns76LyGrU318RSB9gapeU/UNU8Efk9sExMct1zwCNA/gXlzorIOOAVp0OjCLgZeBPz6vq5Uzs9DNzq4bSvA0tFpEAr6ZVW1WViUn595rwKn8Q4pENuZQpEZAqQialxpavqoipdgPPtGwNMF5FQzPPyF0wb5CxMx1ARMFBVi9z+b4+jr6xJYB3wA3WmiAB+AcwUkV9hrs3PK5GwBPhAREbRcDOMNwps9iCLxWLxgH31tlgsFg9YR2mxWCwesI7SYrFYPGAdpcVisXjAOkqLxWLxgHWUFovF4gHrKC0Wi8UD1lFaLBaLB/4f1WNsMRmAR2MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get some tools for plotting simulated vs actual wealth distributions\n",
    "from HARK.utilities import getLorenzShares, getPercentiles\n",
    "\n",
    "# The cstwMPC model conveniently has data on the wealth distribution \n",
    "# from the U.S. Survey of Consumer Finances\n",
    "from HARK.cstwMPC.SetupParamsCSTW import SCF_wealth, SCF_weights\n",
    "\n",
    "# Construct the Lorenz curves and plot them\n",
    "\n",
    "pctiles = np.linspace(0.001,0.999,15)\n",
    "SCF_Lorenz_points = getLorenzShares(SCF_wealth,weights=SCF_weights,percentiles=pctiles)\n",
    "\n",
    "sim_wealth = EstimationEconomy.aLvlNow[0]\n",
    "sim_Lorenz_points = getLorenzShares(sim_wealth,percentiles=pctiles)\n",
    "\n",
    "# Plot \n",
    "plt.figure(figsize=(5,5))\n",
    "plt.title('Wealth Distribution')\n",
    "plt.plot(pctiles,SCF_Lorenz_points,'--k',label='SCF')\n",
    "plt.plot(pctiles,sim_Lorenz_points,'-b',label='Beta-Point')\n",
    "plt.plot(pctiles,pctiles,'g-.',label='45 Degree')\n",
    "plt.xlabel('Percentile of net worth')\n",
    "plt.ylabel('Cumulative share of wealth')\n",
    "plt.legend(loc=2)\n",
    "plt.ylim([0,1])\n",
    "plt.show('wealth_distribution_1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time Preference Heterogeneneity\n",
    "\n",
    "Our specific approach is to replace the assumption that all households have the same time\n",
    "preference factor with an assumption that, for some dispersion $\\nabla$, time\n",
    "preference factors are distributed uniformly in the population between\n",
    "$\\grave{\\Discount}-\\nabla$ and $\\grave{\\Discount}+\\nabla$ (for this reason, the model is referred to as the $\\Discount$-Dist model).  Then,\n",
    "using simulations, we search for the values of $\\grave{\\Discount}$ and\n",
    "$\\nabla$ for which the model best matches the fraction of net worth held by the top $20$, $40$, $60$, and $80$ percent of the population, while at the same time matching\n",
    "the aggregate capital-to-output ratio from the perfect foresight\n",
    "model. Specifically, defining $w_{i}$ and $\\omega _{i}$ as the proportion of total aggregate net worth held by the top $i$ percent in our model and in the data, respectively, we solve the following minimization problem:\n",
    "\n",
    "$$  \\{\\grave{\\Discount}, \\nabla\\}= \\underset{\\{{\\Discount}, \\nabla\\}}{\\text{argmin} }\\Big(\\sum_\\text{i=20, 40, 60, 80}\n",
    "  \\big(w_{i}({\\Discount}, \\nabla)-\\omega _{i}\\big)^{2}\\Big)^{1/2}$$\n",
    "  subject to the constraint that the aggregate wealth (net worth)-to-output ratio in the model matches the aggregate\n",
    "capital-to-output ratio from the perfect foresight model ($\\KLev_{PF}/\\YLev_{PF}$). When solving the problem for the FBS specification we shut down the aggregate shocks (practically, this does not affect the estimates given their small size).\n",
    " \n",
    "$$\\KLev / \\YLev = \\KLev_{PF} / \\YLev_{PF}$$\n",
    "\n",
    "The solution to this problem is $\\{\\grave{\\Discount}, \\nabla\\}=\\{0.9867, 0.0067\\}$\n",
    ", so that the discount factors are evenly spread roughly between 0.98 and 0.99. We call the optimal value of the objective function the 'Lorenz distance' and use it as a measure of fit of the models.\n",
    "\n",
    "The introduction of even such a relatively modest amount of time\n",
    "preference heterogeneity sharply improves the model's fit to the targeted\n",
    "proportions of wealth holdings, bringing it reasonably in line with the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beginning an estimation with the specification name PYbetaDistIndLiquid...\n",
      "getKYratioDifference tried center = 0.95 and got -6.075744713972091\n",
      "getKYratioDifference tried center = 0.995 and got 111.76737279427269\n",
      "getKYratioDifference tried center = 0.9523201058993506 and got -6.055195625199659\n",
      "getKYratioDifference tried center = 0.9736600529496753 and got -5.684198676919181\n",
      "getKYratioDifference tried center = 0.9843300264748376 and got -3.0115687152149486\n",
      "getKYratioDifference tried center = 0.9896650132374187 and got 16.765594594165385\n",
      "getKYratioDifference tried center = 0.9851424119003377 and got -1.7495265290069613\n",
      "getKYratioDifference tried center = 0.9861621794459002 and got 0.464734381802705\n",
      "getKYratioDifference tried center = 0.985948148191699 and got -0.06632048718483752\n",
      "getKYratioDifference tried center = 0.9859748773639189 and got -0.00206330128011345\n",
      "getKYratioDifference tried center = 0.9859757318452261 and got -6.473122304839762e-07\n",
      "getKYratioDifference tried center = 0.9859762318452265 and got 0.0012065867285953047\n",
      "findLorenzDistanceAtTargetKY tried spread = 0.006 and got 19.749358529784313\n",
      "getKYratioDifference tried center = 0.95 and got -6.074329349907804\n",
      "getKYratioDifference tried center = 0.995 and got 153.76485470647563\n",
      "getKYratioDifference tried center = 0.9517101239746658 and got -6.059198957121332\n",
      "getKYratioDifference tried center = 0.9733550619873329 and got -5.677994182319397\n",
      "getKYratioDifference tried center = 0.9841775309936664 and got -1.0649288044380345\n",
      "getKYratioDifference tried center = 0.9865895842994721 and got 5.91923800548585\n",
      "getKYratioDifference tried center = 0.9845453150276249 and got -0.3090101531308731\n",
      "getKYratioDifference tried center = 0.9846882012829227 and got 0.009287781264188943\n",
      "getKYratioDifference tried center = 0.9846840319304342 and got -0.000199768947886092\n",
      "getKYratioDifference tried center = 0.9846845319304346 and got 0.0009373473997769111\n",
      "findLorenzDistanceAtTargetKY tried spread = 0.008 and got 14.267990340611645\n",
      "getKYratioDifference tried center = 0.95 and got -6.07111352958761\n",
      "getKYratioDifference tried center = 0.995 and got 287.31581963388754\n",
      "getKYratioDifference tried center = 0.9509311938534059 and got -6.062792341529118\n",
      "getKYratioDifference tried center = 0.9729655969267029 and got -5.647989738263023\n",
      "getKYratioDifference tried center = 0.9839827984633515 and got 4.3791970758657435\n"
     ]
    }
   ],
   "source": [
    "'''\n",
    "This options file establishes the second simplest model specification possible:\n",
    "with heterogeneity, no aggregate shocks, perpetual youth model, matching net worth.\n",
    "'''\n",
    "do_param_dist = True          # Do param-dist version if True, param-point if False\n",
    "do_lifecycle = False          # Use lifecycle model if True, perpetual youth if False\n",
    "do_agg_shocks = False         # Solve the FBS aggregate shocks version of the model\n",
    "do_liquid = True              # Matches liquid assets data when True, net worth data when False\n",
    "do_tractable = True           # \n",
    "os.chdir(path_to_models)\n",
    "exec(open('cstwMPC_MAIN.py').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot \n",
    "plt.figure(figsize=(5,5))\n",
    "plt.title('Wealth Distribution')\n",
    "plt.plot(pctiles,SCF_Lorenz_points,'--k',label='SCF')\n",
    "plt.plot(pctiles,sim_Lorenz_points,'-b',label='Beta-Dist')\n",
    "plt.plot(pctiles,pctiles,'g-.',label='45 Degree')\n",
    "plt.xlabel('Percentile of net worth')\n",
    "plt.ylabel('Cumulative share of wealth')\n",
    "plt.legend(loc=2)\n",
    "plt.ylim([0,1])\n",
    "plt.show('wealth_distribution_2')"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-",
   "formats": "ipynb,py:percent"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": false,
   "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": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}