{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fifteen-treaty",
   "metadata": {},
   "source": [
    "# Perfect Foresight CRRA Model - Saving Rate\n",
    "\n",
    "[![badge](https://img.shields.io/badge/Launch%20using%20-Econ--ARK-blue)](https://econ-ark.org/materials/perfforesightcrra-savingrate#launch)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "employed-classroom",
   "metadata": {},
   "source": [
    "This notebook demonstrates how to obtain the saving rate of a forward-looking consumer using HARK."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "threaded-apollo",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preliminaries\n",
    "\n",
    "from HARK.ConsumptionSaving.ConsIndShockModel import (\n",
    "    PerfForesightConsumerType,\n",
    ")  # Import the consumer type\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def mystr(number):\n",
    "    return \"{:.4f}\".format(number)\n",
    "\n",
    "\n",
    "# These last two will make our charts look nice\n",
    "plt.style.use(\"seaborn-v0_8-darkgrid\")\n",
    "palette = plt.get_cmap(\"Dark2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "broke-myrtle",
   "metadata": {},
   "source": [
    "## 1. Creating an agent\n",
    "\n",
    "We now import the HARK class that represents a perfect-foresight consumer with CRRA utility, and create an instance of it with the parameters we'd like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "thorough-booth",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up a HARK Perfect Foresight Consumer called PFsavrate\n",
    "\n",
    "# Now we need to \"fill\" our consumer with parameters that allow us to solve the consumer's problem\n",
    "\n",
    "# First we need to set out a dictionary\n",
    "ρ = CRRA = 2.0  # Coefficient of relative risk aversion\n",
    "Rfree = 1.03  # Interest factor on assets\n",
    "β = DiscFac = 0.97  # Intertemporal discount factor\n",
    "LivPrb = [1.0]  # Survival probability\n",
    "Γ = PermGroFac = [1.01]  # Permanent income growth factor\n",
    "# Number of agents of this type (only matters for simulation# Number of periods in the cycle for this agent type\n",
    "AgentCount = 1\n",
    "cycles = 0  # Agent is infinitely lived\n",
    "T_cycle = 1  # Every period is the same\n",
    "\n",
    "# Make a dictionary to specify a perfect foresight consumer type\n",
    "dict_wealth = {\n",
    "    \"CRRA\": CRRA,\n",
    "    \"Rfree\": Rfree,\n",
    "    \"DiscFac\": DiscFac,\n",
    "    \"LivPrb\": LivPrb,\n",
    "    \"PermGroFac\": PermGroFac,\n",
    "    \"AgentCount\": AgentCount,\n",
    "    \"cycles\": cycles,\n",
    "    \"T_cycle\": T_cycle,\n",
    "}\n",
    "\n",
    "# Now lets pass our dictionary to our consumer class\n",
    "PFsavrate = PerfForesightConsumerType(**dict_wealth)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "golden-syndication",
   "metadata": {},
   "source": [
    "## 2. Solving the consumer's problem and finding the saving rate.\n",
    "\n",
    "The next step is to ask HARK to solve the agent's optimal consumption problem. The result, or what we call \"solution\", is a function `cFunc` that tells us what is the optimal amount to consume at each level of normalized market resources.\n",
    "\n",
    "The next cell solves the agent and tests the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "political-country",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "When normalized market resources are m=4, the agent's normalized consumption should be c=1.611193200855211\n"
     ]
    }
   ],
   "source": [
    "PFsavrate.solve()\n",
    "test_m = 4\n",
    "test_c = PFsavrate.solution[0].cFunc(test_m)\n",
    "\n",
    "print(\n",
    "    \"When normalized market resources are m={}, the agent's normalized consumption should be c={}\".format(\n",
    "        test_m, test_c\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "informative-moment",
   "metadata": {},
   "source": [
    "We are now ready to think about the saving rate of the consumer.\n",
    "\n",
    "Remember that the saving rate of an agent at time $t$ is\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Saving Rate}_t = \\frac{\\text{Capital Income}_t + \\text{Labor Income}_t - \\text{Consumption_t}}{\\text{Capital Income}_t + \\text{Labor Income}_t}\n",
    "\\end{equation}\n",
    "\n",
    "Dividing both the numerator and denominator by the agent's permanent income $\\textbf{P}_t$ and adopting the notation in [the course's lecture notes](https://www.econ2.jhu.edu/people/ccarroll/public/LectureNotes/Consumption/PerfForesightCRRA/) we have\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\text{Saving Rate}_t &= \\frac{a_{t-1}*r/\\Gamma + 1 - c(m_t)}{a_{t-1}*r/\\Gamma + 1}\\\\\n",
    "&=\\frac{a_{t-1}*r/\\Gamma + 1 - c(a_{t-1}*R/\\Gamma + 1)}{a_{t-1}*r/\\Gamma + 1}\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "We now have an expression for the saving rate as a function of the previous period's end-of-period normalized assets $a_t$. The cell below calculates and plots this relationship for different levels of $a_t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "future-finder",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Saving Rate')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAGqCAYAAAALLU5kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQqElEQVR4nO3dd3hUZfrG8e8kmfRQUiD0QCAFUKSrWBBEaQEBXVwVdFdcCyuwa0BAJWAUF40FXFxd6+rq6k8FIYKA2AuIFBEhCSR0CCENSJlJm/P7I5qVJWiCmTkzyf25Lq7LnHPm5JknQ7h9zznvazEMw0BERESkCfAyuwARERERV1HwERERkSZDwUdERESaDAUfERERaTIUfERERKTJUPARERGRJkPBR0RERJoMBR8RERFpMhR8REREpMnwMbsAd5ObW+S0c4eGBlFQUOK080s19dk11GfXUJ9dQ312HWf1OiIipE7HacTHRSwW8Pb2wmIxu5LGTX12DfXZNdRn11CfXccdeq3gIyIiIk2Ggo+IiIg0GQo+IiIi0mQo+IiIiEiToeAjIiIiTYaCj4iIiDQZCj4iIiLSZCj4iIiISJOh4CMiIiJNhoKPiIiINBkKPiIiItJkKPiIiIhIk6HgIyIiIk2Ggo+IiIg43c5jRaR8nMmB/BJT6/Ax9buLiIhIo2avqOIfX+3nza1HcBgQ1boZ1/ZoZVo9Cj4iIiLiFN8eLOThdXs4ctIOwIj4Vtx0YSfsRTbTalLwERERkQZVZK9k8ed7WbHjGACtgn2ZOyyGS6JDCfbzwV5kXm1ueY9PaWkpc+bMYeDAgfTt25dZs2ZRUvLr1wSPHz/OxRdfzLJly2q2jRo1it69e5/2JzY2lueee86Zb0FERKRJ+iwzj9+9srkm9Fzbqw1v3dKPQV1CTa6smluO+CQnJ5Odnc3atWupqqpixowZpKSkkJSUdNbXOBwOEhMTKSwsPG37qlWrTvv6qaee4tNPP+Wmm25ySu0iIiJNUX5JOSkfZ7F+dy4AHVsGcN9V3ejTvoW5hf0Ptws+NpuN1NRUXn31VVq0aAFAYmIikydPZtasWQQEBNT6uqVLlxIZGUmbNm3Oeu6NGzfyr3/9i+XLlxMUFOSM8kVERJoUwzD4IO04T3ySxUl7Jd4WuKl/B6Zc2BF/q7fZ5Z3BlOBjt9vJycmpdZ/NZqOiooKYmJiabdHR0djtdvbv3098fPwZr9m4cSOrVq3i3XffJSEhodbzVlVVkZSUxJ133klUVNQv1mex1P291NVP53TGueW/1GfXUJ9dQ312DfX53B07ZWfhh3v4el/11ZaYVkHMuzqGuNYhtR7vDr02Jfhs376dyZMn17pv+vTpAAQGBtZs+2mUp7b7fPLz85k7dy5Lliz5xVGc1NRUSktLz/p9fxIaGoS3t/NufQoLq/3DIA1LfXYN9dk11GfXUJ/rzuEweH3TQf62Oo2S8ip8vb2YfmU3/nRZF6x1+DfUzF6bEnwGDhxIRkZGrft27drF4sWLsdlsNUHGZqt+7C04OPi0Yw3DYNasWUyaNImePXv+4vf8v//7PyZOnIi/v/8vHldQUOK0EZ+wsBDy84swjIY/v1RTn11DfXYN9dk11Of6OVho46G1u9l6+CQA57dtxryrY4gKC+Rk4S8/iOTMXoeH1y1Mud09Pp07d8ZqtZKZmUmvXr0AyMrKwmq1nnGJKjs7m02bNrF9+3aWLl0KQHFxMQsWLGDt2rU1T27l5eWxdetWFi1aVKcanPnBNwznnl+qqc+uoT67hvrsGurzL6t0GPxny2Ge+/oAZZUO/H28mHppZ667oC3eXpZ69c7MXrtd8AkICGDEiBGkpKSwePFiAFJSUhg9evQZozVt27Zlx44dp20bMmQIf/7znxk/fnzNtq1bt9KqVSs6dOjg/DcgIiLSyGTmlpC8bje7jlVPwDOgYwvmXtWNds1rf+DInbnlPD5JSUlERUWRkJDA8OHDad++PfPmzavZP2rUKJ599tk6n+/QoUO0bt3aGaWKiIg0WhVVDv759X4m/Xsru44VEeznzQNXxfD3a8/zyNADYDEMDez9XG6uc6aTtFiqrz/m5ekasjOpz66hPruG+uwa6nPtdh0rInntbjLzqu/buTw6jHuv7EpEsN85n9OZvY6I8NB7fERERMQ89ooqnt9wgH9vPozDgBYBVmYOiWZYbASWRvDMv4KPiIiIAPDd4ZMkr9vNwcLqp6mvjovgniuiaRnoa3JlDUfBR0REpIkrLa/imS/38X/bjmIA4UG+zL6yG5d3DTO7tAan4CMiItKEbTpQyMPrdnP0VBkAY3q2Zsbl0YT4N86I0DjflYiIiPyi4rJKFn+2l/d+XEU9MsSP+67qxoVR7rGKurMo+IiIiDQxX+0tYOGHuzleXA7AdRe0ZeqlUQT5Nv5Y0PjfoYiIiABw0lbBk59msWrXcQDat/Dn/qti6NuhhbmFuZCCj4iISBPw6Z48/vZRJvkl5ViA3/dtx52DovC3eptdmksp+IiIiDRiJ0oreOzjTNZl5AIQFRrAA1fHcn7bZiZXZg4FHxERkUbIMAzW787jsY8yKbRV4GWBSf07cNtFnfDzccsVq1xCwUdERKSRyS8pZ9FHmXyyJw+A6PBA5l0dS/fIui3r0Jgp+IiIiDQShmGwNj2XlI8zOWmvxNvLwh8GdOCPF3bE6t10R3l+TsFHRESkEcgrLuNv6zP5LCsfgJiIIOYNjyW2VbDJlbkXBR8REREPZhgGH6Qd5/FPsjhlr8THy8KUizpyc/8O+GiU5wwKPiIiIh7qeFEZj6zfw5d7CwCIbx3MvKtj6RoRZHJl7kvBR0RExMMYhsH7O3N44tMsisuqsHpbuO2iTkzq3wEfL4vZ5bk1BR8REREPklNUxsIPd/P1vkIAekSG8MDVMUSHa5SnLhR8REREPIBhGKT+UD3KU1Jeha+3hdsvjuKGfu01ylMPCj4iIiJu7tgpOws/3MOG/dWjPOe1CeGBq2PpHBZocmWeR8FHRETETRmGwcofjvHkp3trRnnuGBTFDX3b461RnnOi4CMiIuKGjp2y8/CHe9j4s1GeeVfHEqVRnt9EwUdERMSNGIbBih3HeOqz6lEePx8v7hgUxe/7tNMoTwNQ8BEREXETx07ZeXjdHjYe+GmUpxnzhscQFapRnoai4CMiImKy/72XR6M8zqPgIyIiYqKcojIeXrf7Z09saZTHmRR8RERETGAYBqk7c3jikyw9seVCCj4iIiIudryojIUf7uGrfdVrbPVsE0KSnthyCQUfERERFzEMg1W7cnj8k+o1tn6affnGfhrlcRUFHxERERfILa4e5flpJfUekSEkDdfsy66m4CMiIuJEhmGwJv04KR9nccpeidXbwp8u6sRNWkndFAo+IiIiTpJfUs7f1u/h08x8AOJbB5M0PFYrqZtIwUdERMQJPszIZdH6PZy0V+LjZeG2izoxeYBGecym4CMiItKATpRWsOijPazfnQdATEQQ80fE0i0i2OTKBBR8REREGswne/L42/o9FJRW4O1l4Y8DO/DHgR3x8fYyuzT5kYKPiIjIb3TSVsFjH2eyNj0XgOjwQOYPjyWudYjJlcn/UvARERH5Db7IyufhD/eQX1KOlwVuHtCBKRd2wtdHozzuSMFHRETkHBSXVfLEJ1mk7swBICo0gPnDY+nRppnJlckvUfARERGpp2/2F5K8bjc5RWVYgBv7teeOQVH4aZTH7bnlT6i0tJQ5c+YwcOBA+vbty6xZsygpKfnV1x0/fpyLL76YZcuW1WxzOBw8+eSTXHbZZfTt25ff/e53bNq0yZnli4hII1VaXsXf1u/hz+/uIKeojPYt/PnnxF5Mv7yLQo+HcMufUnJyMtnZ2axdu5Z169aRnZ1NSkrKL77G4XCQmJhIYWHhadvffPNN1q9fz9tvv823337LyJEjuf322ykrK3PmWxARkUZm6+ET3PDqFt7dng3A7y5oyxuT+3JB++YmVyb14XbBx2azkZqayrRp02jRogVhYWEkJiaybNkybDbbWV+3dOlSIiMjadOmzWnb9+7di8PhwOFwYBgGFosFf39/Z78NERFpJOwVVTz5aRZ3vPU9R07aiQzxY+m15zFzaFcCrN5mlyf1ZMo9Pna7nZycnFr32Ww2KioqiImJqdkWHR2N3W5n//79xMfHn/GajRs3smrVKt59910SEhJO23f99dfz0UcfMXjwYLy9vfHz8+Of//wnfn5+DfumRESk0fkh+xTzP8jgQGH1/3iP7RnJjMFdCPbTLbKeypSf3Pbt25k8eXKt+6ZPnw5AYOB/V6sNCAgAqPU+n/z8fObOncuSJUsICjpz7ZOKigoGDBjA7bffTtu2bXnxxReZNm0aK1euJCIiotYaLE6YTfynczrj3PJf6rNrqM+uoT67Rm19rqhy8M+vD/CvTYdwGBAe5Mv9V8dwSZdQc4psJNzhM21K8Bk4cCAZGRm17tu1axeLFy/GZrPVBJmfLnEFB58+3bdhGMyaNYtJkybRs2fPWs83a9Ys7rjjDrp06QLA1KlTWbFiBWvWrGHSpElnHB8aGoS3E2fYDAvTZFauoD67hvrsGuqza/zU511HT/HX/9tO+rEiAK65oC3zx/SgRaCvmeU1KmZ+pt1urK5z585YrVYyMzPp1asXAFlZWVitVqKiok47Njs7m02bNrF9+3aWLl0KQHFxMQsWLGDt2rU899xzHD16lPLy8tNe5+Pjg9VqrfX7FxSUOG3EJywshPz8Igyj4c8v1dRn11CfXUN9do2f+pyTe4p/fXOIf359gEqHQYsAK7Ov7MqVsRFUlpaRV6qHYn4rZ36mw8PrFqbcLvgEBAQwYsQIUlJSWLx4MQApKSmMHj36jJuS27Zty44dO07bNmTIEP785z8zfvz4mq//8Y9/0K9fPyIjI3njjTfIzc3liiuuOGsNzvwFYxjOPb9UU59dQ312DfXZ+bJyi5n+xnf8kF09ynN5dBhzhnUjLMhXvXcCMz/TbvdUF0BSUhJRUVEkJCQwfPhw2rdvz7x582r2jxo1imeffbZO55o/fz6XXXYZN954IxdffDEffvghL774Iq1bt3ZW+SIi4iEchsF/thxh5OIv+CG7iGA/b+YPj+Wxsd0JC9KlrcbIYhjKsj+Xm1vklPNaLNXDcHl5GrJ2JvXZNdRn11CfnevoSTsPrs1gy6GTAAzs1IL7r4ohspmmPHEWZ36mIyI89FKXiIiIMxmGwcofjvHkp3spKa/C38eL+0Z3Z3h0S0CP0DV2Cj4iItJk5JWU8/C63Xy5twCAXm2bsWBkLBd0baWRtSZCwUdERJqE9Rm5/G39Hk7aK7F6W7hzUBQ39G2Pj7dGeZoSBR8REWnUTtoqeOzjTNam5wIQ2yqY+SNi6Rp+5qS30vgp+IiISKP19b4CHlq3m9zicrwtcPPAjky5sCNWJ05UK+5NwUdERBqd0vIqFn+2l2XfV6+k3qllAAtGxNKjTTOTKxOzKfiIiEijsv3ISZI+yODISTsA1/dpx9RLovDXSuqCgo+IiDQS5ZUO/rnhAK99W72waOsQP5KGx9C/Y0uzSxM3ouAjIiIeb09uMUkfZLAntwSA0T1ac88V0QT76Z85OZ0+ESIi4rGqHAb/3nyYZ7/aT6XDoGWAlbnDujG4W7jZpYmbUvARERGPdPiEjfkfZLD96CmgemHRuVd1IzRQa2zJ2Sn4iIiIRzEMg+U7jvHUp1nYKhwE+XpzzxXRjO7RGotFkxHKL1PwERERj5FXXEbyut18va8QgL4dmpM0PJY2WlhU6kjBR0REPMJHu3N55MPqJSd8vS1MvbQz1/dph5dGeaQeFHxERMStFdkreezjTD5IOw5ULzmxYEQs0VpyQs6Bgo+IiLitbw8WsmDNbnKKyvCywC0DOjDlok5ackLOmYKPiIi4HXtFFc98uZ//bD0CQPsW/iwYEcf5bbXkhPw2Cj4iIuJW0nOKmPdBBvvySwEYf34bpl/ehUBfLTkhv52Cj4iIuIVKh8Fr3x7iua8PUOUwCAvy5YGrYhjUJdTs0qQRUfARERHTHT5hI+mDDL7/cTLCId3CmXNlN1oEWk2uTBobBR8RETGNYRis2HGMJ342GeHMIV0Z2b2VJiMUp1DwERERU+SXlPPwut18sbcAgD7tmzN/hCYjFOdS8BEREZf7LDOfh9ftptBWgdXbwp2Dorihb3u8vTTKI86l4CMiIi5TUl7Jk5/sZcUPxwDoGh7EgyNj6RYRbHJl0lQo+IiIiEtsP3KSpA8yOHLSjgW4sV977hwUha+PJiMU11HwERERp6qscvD8hgO8sukQDgMiQ/yYPyKWvh1amF2aNEEKPiIi4jT780uZ90E6aTnFAIzs3oqZQ7oS7Kd/fsQc+uSJiEiDMwyDt7/LZsnneymrdNDM34c5V3bjytgIs0uTJk7BR0REGlRecRkPrt3Nhv2FAAzs1IJ5V8fSKsTP5MpEFHxERKQBfbInj4fX7eakvRJfbwt3X9aF3/Vui5cmIxQ3oeAjIiK/WUl5JY9/nEXqzhwAYiKCeHBkHNHhQSZXJnI6BR8REflN/vcx9Un9O3D7xZ30mLq4JQUfERE5J5VVDp7feJBXvjmIw4A2zaofU+/TvoXZpYmclYKPiIjU24GCUuZ9kMGuY0WAHlMXz6FPqIiI1JlhGCzfcYwnP8nCXukgxM+H2Vd25aq4VmaXJlInCj4iIlInhaXlPLRuD59n5QPQr0NzkobHEqnV1MWDKPiIiMiv+mpvAQ+uzaCgtHo19bsu6cwNfdvpMXXxOAo+IiJyVvaKKhZ/tpd3tmcD0CUskOSRccS00mrq4pkUfEREpFYZOcXcvzqN/QU2AK7v046pl0Thb/U2uTKRc+eWkyyUlpYyZ84cBg4cSN++fZk1axYlJSW/+rrjx49z8cUXs2zZsppthmHw/PPPM2TIEPr06cMtt9zC7t27nVm+iIhHcxgGr246xC1vbGN/gY3wIF+entCTe66IVugRj+eWwSc5OZns7GzWrl3LunXryM7OJiUl5Rdf43A4SExMpLCw8LTtr732Gi+88AIpKSls2rSJoUOHMnnyZAoKCpz5FkREPNKxU3amvv09T3+xj0qHweCuYfxncl8ujAo1uzSRBuF2wcdms5Gamsq0adNo0aIFYWFhJCYmsmzZMmw221lft3TpUiIjI2nTps1p299//30mTZpEnz598PHxYdKkSbRs2ZI1a9Y4+62IiHiUDzNyueHVrWw+dBJ/Hy/uG9aNR8d0p0Wg1ezSRBqMKff42O12cnJyat1ns9moqKggJiamZlt0dDR2u539+/cTHx9/xms2btzIqlWrePfdd0lISDhtX1VVFYGBgadt8/LyYu/evWetzxkPKfx0Tj0A4Vzqs2uoz67hqj4Xl1WS8nEW7/+4zlb3yBAeGhVHx5YBzv3GbkKfZ9dxh16bEny2b9/O5MmTa903ffp0gNPCSkBA9V++2u7zyc/PZ+7cuSxZsoSgoDMXw7v66qt57bXXuOiii+jatSvvvPMO+/bto0+fPrV+/9DQILy9nTcQFhYW4rRzy3+pz66hPruGM/u85UABM976jkMFNrwscNfgrky/shtWJ/4edFf6PLuOmb02JfgMHDiQjIyMWvft2rWLxYsXY7PZaoLMT5e4goNPf3zSMAxmzZrFpEmT6NmzZ63n++Mf/4jNZmPq1KmUl5czYsQILrnkEpo1a1br8QUFJU4b8QkLCyE/vwjDaPjzSzX12TXUZ9dwZp8rHQYvbjjAixv/u85W8sg4LmjfnJOFv/4wSWOiz7PrOLPX4eF1C1Nu9zh7586dsVqtZGZm0qtXLwCysrKwWq1ERUWddmx2djabNm1i+/btLF26FIDi4mIWLFjA2rVree6558jJyeHaa6+tGUmqrKxkyJAhjBs37qw1OPODbxjOPb9UU59dQ312jYbu8+ETNuatzmBH9ikARsS3YtbQ6nW2mvLPU59n1zGz124XfAICAhgxYgQpKSksXrwYgJSUFEaPHo2//+nTordt25YdO3actm3IkCH8+c9/Zvz48QCsWrWKVatW8corr+Dr68vf//53fH19GTJkiGvekIiImzAMgw/SjvPoR5mUlFcR7OfN7KHduDpe62xJ0+F2wQcgKSmJRYsWkZCQQEVFBUOHDuWBBx6o2T9q1CgSEhK44447fvVcf/jDH8jOzmbkyJFUVFTQr18/XnnlFfz8/Jz5FkRE3EqRvZK/rd/DuoxcAHq3a8aCkXG00Tpb0sRYDEMDez+Xm1vklPNaLNXXH/PydA3ZmdRn11CfXaOh+rz18AmSVmdwrKgMbwv86eIobh7QAW8vPcYE+jy7kjN7HRHhoff4iIhIw6iscvD8hgO8sukQDgPat/AneWQcPdvU/nCHSFOg4CMi0ggdLLTxwOp0dh2rHsVO6NGae4ZEE+SrX/vStOlvgIhII2IYBqk7c0j5OBNbhYMQPx/mDuvGlbERZpcm4hYUfEREGolT9goe+TCT9burb2Du0745C0bEEqkbmEVqKPiIiDQCWw+fYN7qDHKKyvD2snDHxZ2Y1F83MIv8LwUfEREPVlnl4PmNB3nlm+oZmDu08Cd5VDw9IrX8gkhtFHxERDzU4RPVNzD/kF19A/PoHq1J1A3MIr9IfztERDxMbTMwzx0WwzDdwCzyqxR8REQ8SHFZ9QzMa9Orb2C+oF0zHtQMzCJ1puAjIuIhvj96igdWpXH0VPUMzLdd3IlbBnTUDcwi9aDgIyLi5qocBi9tPMgLGw5QZUDb5tUzMJ/fVjMwi9SXgo+IiBs7csLGn/9vO9sOnwLg6rgIZl/ZjWA//foWORf6myMi4qbWZ+Sy8MM9nLJXEmj15t4ruzKye2uzyxLxaAo+IiJuxlZRxeMfZ7Hih2MA9IgM4aFRcbRvEWByZSKeT8FHRMSNpOcUcd+qdA4W2rAAdw6OZnLvNnh7eZldmkijoOAjIuIGHIbBf7Yc4e9f7KPSYRAR7MuDI2MZ0acjeXlFGIbZFYo0Dgo+IiImyy8pZ/6aDDbuLwTg8ugw7r86hpaBVpMrE2l8FHxEREy0YX8B8z/IoKC0Aj8fL2Zc3oUJvdpgsWhuHhFnUPARETFBeaWDpV/u440tRwCIDg/koVHxdA0PMrkykcZNwUdExMUOFJRy/6p00o8XA3BtrzZMv7wL/lZvkysTafwUfEREXMQwDN7fmcNjH2diq3DQ3N+HB66O4fKu4WaXJtJkKPiIiLjA/y4u2rdDcxaMiKN1iJ/JlYk0LecUfHbt2sXhw4cZPHgwRUVFhIWFNXRdIiKNxg/Zp7hvVTpHT9rxtsCfLo7i5gEdtLioiAnqFXzy8/OZOnUqP/zwA1arlXfeeYdrr72Wl156id69ezurRhERj+QwDF7ddIhnvz5AlcOgTTM/HhoVr8VFRUxUr6lAFy5cSExMDN9++y0+Pj5ER0fzpz/9iUcffdRZ9YmIeKS84jLufmcHS7/cT5XD4MqYCF6f1FehR8Rk9Rrx2bhxI+vXrycgIKBmjokpU6bw0ksvOaU4ERFP9NW+AhZ8kEGhrXpunplDohnTM1Jz84i4gXoFH6vVit1uJyAgAOPH+dNLSkoICtK8EyIi/zs3T7eIIB4eFU/nsECTKxORn9TrUteQIUOYOXMm+/fvx2KxkJ+fz4IFC7jsssucVZ+IiEc4WGjj1v98VxN6JvZuy8s39FboEXEz9Qo+99xzD4GBgQwfPpxTp05xySWXYLPZSExMdFZ9IiJub/WuHCa9tpX048U09/chZWwPEod0xc9HK6qLuJt6XeoKCgpiyZIl5Ofnc+TIESIjI2nVqhXFxcXOqk9ExG2VlFfy6EeZrN51HIA+7ZuTPDKOVpqbR8Rt1Sv4DBgwgE2bNhEWFnba3D2DBw9m8+bNDV6ciIi7Ss8p4r5V6RwstOFlgdsu6sQfBnbU3Dwibu5Xg8+BAweYN28ehmFQXFzM5MmTT9tfXFxMs2Z6PFNEmgbDMHhz21GWfLaXSodB6xA/HhoZxwXtm5tdmojUwa8Gn06dOnHVVVdRWFjI1q1bGTBgwGn7fX19GTJkiNMKFBFxFydKK1iwNoMv9xYAMLhrGPdfFUPzAKvJlYlIXdXpUteNN94IQPv27bnmmmucWY+IiFvacugED6xOJ7e4HF9vC9Mvj+a6C9pobh4RD1Ove3yuueYaNm7cSE5OTs08PhUVFWRkZHD//fc7pUARETNVOgxe3HCAFzcexAA6tQxg4eh4YloFm12aiJyDegWfhx56iDfffLNmwsKqqipKSkq49NJLnVKciIiZcorKeGBVGtuOnAJgTM/WJA7pSoDV2+TKRORc1Sv4fPDBB/z73//GZrOxcuVKFi5cyKJFiygtLXVWfSIipvg8K58H12Rw0l5JkK83s6/sxvD4VmaXJSK/Ub2Cj81m44ILLiA3N5edO3disVj485//zMiRI51Vn4iIS5VXOnj6i328ubV6Bub41sE8PCqeDi0DTK5MRBpCvaYVjYyMJD8/n4iICI4dO0ZFRQX+/v4NPoFhaWkpc+bMYeDAgfTt25dZs2ZRUlJy1uOTkpLo2bMnvXv3rvnz1ltv1exfvnw5w4YN44ILLmD8+PFs27atQesVkcbh0I/LTvwUem7o244Xf3+BQo9II1KvEZ/LL7+cW265hX/961/079+fuXPn4ufnR1RUVIMWlZycTHZ2NmvXrqWqqooZM2aQkpJCUlJSrcfv2LGD5ORkxo0bd8a+b775huTkZJ5//nnOP/98Xn/9de68804++eQTAgL0y0xEqq1NO84j6/dQUl5Fc38fkobHcml02K+/UEQ8Sr1GfP76178yduxYrFYr8+bN48SJE2RmZpKcnNxgBdlsNlJTU5k2bRotWrQgLCyMxMREli1bhs1mO+P48vJydu/eTc+ePWs939tvv82oUaPo27cvVquVW265hZYtW7J69eoGq1lEPJetoorktRncvzqdkvIqerdrxuuT+yr0iDRS9RrxsVqtTJkyBYCQkBCef/55AA4fPlyvb2q328nJyal1n81mo6KigpiYmJpt0dHR2O129u/fT3x8/GnHp6enU1lZyZIlS9iyZQshISFMmDCBKVOm4OXlRWZmJhMmTDjtNV27diU9Pf2s9TljWo6fzqkpP5xLfXaNxtLnzNwS5qSmsa+gFAsw5aKO3HpRJ3zcZNmJxtJnd6c+u4479LrOwefdd98lLS2Nyy+//LTH19977z0eeuiheq3VtX379jOWvvjJ9OnTAQgMDKzZ9tMlqdru8ykqKmLAgAFMmjSJJ554grS0NKZOnYqXlxdTpkyhpKTkjEta/v7+Z30SLTQ0CG9v562oHBYW4rRzy3+pz67hqX02DIM3vz3E/JU7Kat00CrEj6euv4CLo8PNLq1WntpnT6M+u46Zva5T8HniiSd4+eWXiYmJ4c0332Tp0qVccsklzJ07lxUrVnDTTTfV65sOHDiQjIyMWvft2rWLxYsXY7PZauYL+ukSV3DwmROGDRo0iEGDBtV8ff7553PzzTezevVqpkyZQkBAAHa7/bTX2O12WrZsWev3LygocdqIT1hYCPn5Rfw496M4gfrsGp7c5+KySh5et4cPM3IBuDiqJfNHxBIa5EteXpHJ1Z3Ok/vsSdRn13Fmr8PD6xam6hR83n//fZ5++mkGDx7MBx98wH/+8x/effddNm/ezHPPPcfll1/+m4r9uc6dO2O1WsnMzKRXr14AZGVlYbVaa72Jev369eTl5XH99dfXbCsvL8ff3x+Abt26sWfPntNek5mZyWWXXXbWGpz5wTcM555fqqnPruFpfU7LKWLu+2kcPmHH28vC1EuiuLFfe7wsFrd+H57WZ0+lPruOmb2u0zWdwsLCmnBz5ZVX8uWXX3LgwAGWL1/eoKEHqi9rjRgxgpSUFAoKCigoKCAlJYXRo0fXhJmfMwyDRx55hA0bNmAYBtu2bePVV19l4sSJAFx77bWkpqayceNGKioqeOWVV8jPz2fYsGENWreIuC/DMHhz6xH++MZ3HD5hp00zP56f2ItJ/TvgpRs7RJqUOo34WCyWmoX4rFYrXl5ePPnkk7Ru3dopRSUlJbFo0SISEhKoqKhg6NChPPDAAzX7R40aRUJCAnfccQfDhg1jzpw5zJ8/n5ycHMLDw7n77rsZO3YsABdddBFJSUk1+7t27crzzz9PixYtnFK7iLiXk7YKktfu5rOsfKB6RfUHro6hmb9WVBdpiiyG8euDTX369GHr1q01X/fv359vv/3WqYWZJTfXOdf4LZbq6495ebqG7Ezqs2t4Sp+3HznJ/avSOVZUhtXbwozLu3DdBW09ZkV1T+mzp1OfXceZvY6IaMB7fP6Xp/zSEJGmyWEYvPbtYf7x5T6qDOjQwp+Fo+OJa62ndkSaujoFH5vNxtChQ2u+LioqOu1rgI8++qhhKxMROQeFpeUkfZDBhv2FAFwdF8HsK7sR7HdO/58nIo1MnX4TLFy40Nl1iIj8ZlsOneCB1enkFpfj5+NF4hXRjD0vUqPUIlKjTsGntjWwRETcRZXD4OVvDvL8hgM4DOgcGsjChHi6hgeZXZqIuBmN/YqIR8svKeeB1el8e/AEAKN6tObeoV0JsHqbW5iIuCUFHxHxWJsOFPLA6nQKSivw9/Hi3iu7MrpHpNlliYgbU/AREY9T5TB4YcMBXtx4EAPoEhbIIwnxdAnTpS0R+WUKPiLiUXKLy3hgdTpbDp0EYGzPSBKHROOvS1siUgf1Cj6TJk2q9ekIq9VKaGgoV1xxBSNHjmyw4kREfm7j/gLmrc6g0FZBgNWLOcO6MSLeOTPIi0jjVKe1un7Sq1cv0tLSOO+88xg5ciQXXHABGRkZhIaGEh4ezsMPP8xrr73mrFpFpImqdBj846v9THv3BwptFXSLCOLVm/oo9IhIvdVrxGfr1q384x//oF+/fjXbhg4dymOPPcZjjz3G2LFjmT59OpMmTWrwQkWkacotLuP+VelsPVx9aWvc+ZH8dbAubYnIualX8Nm9ezd9+vQ5bdt5553Hrl27AIiLiyM3N7fhqhORJu2b/dVPbRXaKgi0ejNnWDeGx7cyuywR8WD1utTVoUMH3n333dO2paam0rZtWwB27txJREREw1UnIk1SlcPg2a/2c/e7O2oubf3rpt4KPSLym9VrxGfmzJnceeedvPvuu7Rr146jR4+Snp7OkiVLSEtL46abbuK+++5zVq0i0gTkFZdx/8+e2tKlLRFpSPUKPhdffDGrVq0iNTWVY8eOccUVV/DUU0/RunVrjh07xhtvvEF8fLyzahWRRu6bA4XM+3FCQl3aEhFnqPc8Pu3bt+fOO+88Y3tkZCSRkZoxVUTq738nJOwWEcTC0fFEhQaaXZqINDL1Cj579uzh0UcfZf/+/TgcjtP2ffTRRw1amIg0DXk/rrW1+ce1tq45L5J7rtClLRFxjnoFn3nz5hEQEMCf/vQnfHw06bOI/DabD57g/tXp5JeUa0JCEXGJeqWXjIwMPv/8c4KDg51Vj4g0AQ7D4KWNB3l+wwEcBkSHB/K30d2JCtOlLRFxrnoFn1atWlFeXu6sWkSkCSgsLWfe6gw2HigEIKFHa2YN7apLWyLiEvUKPjfddBNTp05l8uTJhIeHn7avf//+DVqYiDQ+2w6f5L5VaeQWl+Pn48XsK7syuoceihAR16lX8HnooYcA2LZt22nbLRYLaWlpDVeViDQqDsPgtW8P848v91FlQOfQQB5JiCc6PMjs0kSkialX8ElPT3dWHSLSSJ2wVTD/gwy+2lcAwMjurbh3aDcCfXVpS0Rcr07B59ixY0RGRnL06NGzHvPTshUiIj/5IfsUc1LTOFZUhp+PF4lXRDP2vEgsFovZpYlIE1Wn4DNy5Ei2bt3KkCFDsFgsGIYBUPPfutQlIj9nGAZvbjvKks/2Uukw6NgygEdGxxPTSk+Eioi56hR8Vq1aBWiSQhH5dcVllTy4djef7MkD4MqYcO67KoZgP839JSLmq9NvojZt2gDw0ksvcd111xEXF+fUokTEM2XkFDP7/V0cPmHHx8vCXwZ34boL2urSloi4jXr9L1h+fj4TJ04kOjqa6667jtGjRxMSEuKs2kTEQxiGwfIdx3j840zKqwzaNPPjkYTu9IjU7wcRcS/1Cj5PPfUURUVFpKamsnz5chYtWsRVV13Fddddp3l8RJqo0vIqHlm/hzVpxwG4tEsoScNjaR5gNbkyEZEz1fuie0hICDfccAM33HADGzZs4L777iM1NVU3N4s0QXvzSrh3ZRr7CkrxtsDUSztzU7/2urQlIm6r3sGnpKSENWvW8N577/H9998zePBgkpOTnVGbiLixFd8d4d53vsde6SAi2JeHR8XTu31zs8sSEflF9Qo+99xzDx9//DGRkZFcd911LF68mNDQUGfVJiJuqKzSwZOfZvHu9mwA+ndswUOj4ggN9DW5MhGRX1ev4OPj48Pzzz9Pv379nFWPiLixIydtzF6ZRvrxYiwWuHVgR6Zc1AlvL13aEhHPUK/gs2jRojO2VVZWsnv3brp3795gRYmI+/ksM58FazIoKqukeYAPS37fhx6h/vw4n6mIiEeoV/D57LPPmD9/Pjk5OTWzN0P1SNCOHTsavDgRMV+lw+AfX+7j1W8PA3BemxAeSYinZ5cI8vKKTK5ORKR+6hV8HnvsMa666iqaNWtGRkYGo0ePZunSpVx77bXOqk9ETJRXXMbcVelsO3wSgN/3acfdl3XG18fL5MpERM5NvYLPoUOHmDlzJocPH2bjxo1cddVVdOnShb/85S9MmjTJWTWKiAk2HzzBfavSKCitIMjXm3lXxzAkJsLsskREfpN6BZ/Q0FC8vLxo27YtWVlZAHTt2pVjx445pTgRcT2HYfCvTYd49qv9OAzoFhHE3xK607FlgNmliYj8ZvUar46NjWXx4sUAhIWF8dlnn/HNN9/g5+fnlOJExLVO2Su4572dPPNldegZ3aM1L/3+AoUeEWk06hV8Zs6cyfr168nNzWXatGncdddd3HLLLdx6660NWlRpaSlz5sxh4MCB9O3bl1mzZlFSUnLW45OSkujZsye9e/eu+fPWW2+dcdzLL7+sS3IiZ5GWU8Sk17by5d4CfL0t3H9VN+ZdHYO/1dvs0kREGky9LnVFR0ezatUqANq1a8cnn3xCSUkJnTt3btCikpOTyc7OZu3atVRVVTFjxgxSUlJISkqq9fgdO3aQnJzMuHHjat1fWlrKkiVLePnllxkwYECD1iri6QzDYPn32aR8kkVFlUG75v4sSuhObOtgs0sTEWlwdR7xcTgcFBYW1ny9cePGmhDUkGw2G6mpqUybNo0WLVoQFhZGYmIiy5Ytw2aznXF8eXk5u3fvpmfPnmc959ixY8nNzeX3v/99g9cr4slsFVXMX5PBI+szqagyuCw6jNdu6qPQIyKNVp1GfHJycvjjH//I+eefzyOPPEJqair33nsvcXFxLF26lJdffpnzzjuvzt/UbreTk5NT6z6bzUZFRQUxMTE126Kjo7Hb7ezfv5/4+PjTjk9PT6eyspIlS5awZcsWQkJCmDBhAlOmTMHLqzrXvfbaa0RGRvL000/X3JT9S5yxvuJP59Tajc6lPtfd/oJS7l25i6y8Urx+XGB0cv+6LTCqPruG+uwa6rPruEOv6xR8nnzySWJjY0lMTATg6aef5rbbbuMvf/kLK1eu5Omnn+af//xnnb/p9u3bmTx5cq37pk+fDkBgYGDNtoCA6hsra7vPp6ioiAEDBjBp0iSeeOIJ0tLSmDp1Kl5eXkyZMgWAyMjIOtcWGhqEt7fz5igJCwtx2rnlv9TnX7Z6RzYz395OSXkV4cF+/P2G3lzYJaze51GfXUN9dg312XXM7HWdgs9XX33FihUrCA0N5ejRoxw8eJAxY8YAMHToUB566KF6fdOBAweSkZFR675du3axePFibDYbQUFBADWXuIKDzxx+HzRoEIMGDar5+vzzz+fmm29m9erVNcGnPgoKSpw24hMWFkJ+fpGm+Hci9fmXVVY5ePrzfby+5QgAfdo3Z+HoOMKDfes1C7P67Brqs2uoz67jzF6Hh9ctTNUp+BQXF9eswr59+3aaNWtGdHQ0AH5+flRUVJxjmWfq3LkzVquVzMxMevXqBUBWVhZWq5WoqKgzjl+/fj15eXlcf/31NdvKy8vx9/c/5xqc+cE3DOeeX6qpz2fKLS5j7vtpfHfkFACT+rXnrks74+NlOedeqc+uoT67hvrsOmb2uk7XdJo3b05BQQEAmzZtok+fPjX79u7dS8uWLRusoICAAEaMGEFKSgoFBQUUFBSQkpLC6NGjaw0zhmHwyCOPsGHDBgzDYNu2bbz66qtMnDixwWoS8XRbDp3gpte28t2RUwT5evPomO5Mu7wLPlpVXUSamDoFnyuuuILk5GRWr15Namoqo0aNAuDUqVMsXryYSy+9tEGLSkpKIioqioSEBIYPH0779u2ZN29ezf5Ro0bx7LPPAjBs2DDmzJnD/Pnz6d27NzNnzuTuu+9m7NixDVqTiCcyfpyF+a63v6egtIKu4UG8elMfrugWbnZpIiKmsBjGrw82nTp1ihkzZrB161ZGjRrFww8/DEDv3r2JiIjgjTfeIDy8cfwizc11zmrTFkv19ce8PF1Ddib1+b+K7JUsWJPBZ1n5AIzq3orZV3ZrkAkJ1WfXUJ9dQ312HWf2OiKiAe/xadasGS+99NIZ259++mn69++vJStE3Mzu48Xcm7qLwyfsWL0tJA7pyrjzIuv0qLqISGNWr5mb/9cll1zSUHWISANZtTOHR9bvoazSQZtmfvwtoTvdI/WYrogI/MbgIyLuo7zSwROfZvHu9mwALopqyYMj42gRYDW5MhER96HgI9IIHDtl597UNHYdK8ICTLmoI1Mu6oSXLm2JiJxGwUfEw32zv5D7VqVx0l5JM38fHhwZx6DOoWaXJSLilhR8RDyUwzB45ZtDPPvVfgwgrlUwi8Z0p23zc5+8U0SksVPwEfFAp+wVzP8ggy/2Vk8sOva8SGYO6Yqfj/PWmRMRaQwUfEQ8TMbxYu5duYsjJ+34elu4d2g3xpxX94V4RUSaMgUfEQ/y80fV2zbzY9GY7sS11qPqIiJ1peAj4gH+91H1izu35MERcTTXo+oiIvWi4CPi5nKKypiduosfsqsfVb/tok7celFHPaouInIOFHxE3Ni3Bwu57/10Cm0VelRdRKQBKPiIuCHDMHjt28Ms/XIfDgNiIoJYNKY77VsEmF2aiIhHU/ARcTPFZdWrqn+a+eOq6j1aM3to1wZZVV1EpKlT8BFxI1l5JcxauYuDhbbqVdWviGbc+W20qrqISANR8BFxEx9m5JK8NgNbhYNWwb48OqY7Pdo0M7ssEZFGRcFHxGSVDoOnP9/LG1uOANCvYwsWjoqjZaCvyZWJiDQ+Cj4iJiooLWfu+2lsOXQSgMn9O3DnJVH4eOnSloiIMyj4iJhkx9FTzE7dxfHicgKt3iQNj2FITITZZYmINGoKPiIuZhgGy7/P5rGPs6h0GESFBvDomB50Dgs0uzQRkUZPwUfEhewVVTz6USapO3MAuKJbOPOujiHYT38VRURcQb9tRVwk+5SdWSt2kX68GC8L3HVJZyb3b69H1UVEXEjBR8QFvjlQyH3vp3HSXklzfx8eHh3PwE4tzS5LRKTJUfARcSLDMHj128M88+PSE/Gtg1k0pjttmvmbXZqISJOk4CPiJCXllSSv3c1Hu/MAGNOzNbOGdsPPx8vkykREmi4FHxEnOFBQysyVu9iXX4qPl4WZQ7T0hIiIO1DwEWlgn2Xmk/RBOiXlVUQE+7IooTvntdXSEyIi7kDBR6SBOAyD578+wAsbDwLQu10zFiZ0JzxIS0+IiLgLBR+RBnDKXsG81Rl8ta8AgIm92zLj8i74eOt+HhERd6LgI/IbZeaWMHPlTg6fsOPn48XcYd0Y2b212WWJiEgtFHxEfoMPM3J5cE0G9koHbZv58eiYHsS2Dja7LBEROQsFH5FzUOkweOaLfby2+TAAAzu14KFR8bQIsJpcmYiI/BIFH5F6OlFawdxVaXx78AQAk/t34K5LovD20qPqIiLuTsFHpB4ycoqZuXIn2afKCLB6Me/qWK6MjTC7LBERqSMFH5E6Wr0rh4Uf7qGs0kH7Fv48NrYHXcODzC5LRETqQcFH5FdUVjlY/Pk+3tx6BIBBnUN5cGQszfx1P4+IiKdR8BH5BQWl5cxJTWPr4ZMA3HphR/50cSe8tPSEiIhHUvAROYtdx4qYuWInx4vLCfL1Zv7wWAZ3Cze7LBER+Q3cclrZ0tJS5syZw8CBA+nbty+zZs2ipKTkrMcnJSXRs2dPevfuXfPnrbfeAsAwDJYuXcqQIUPo06cPCQkJrFmzxlVvRTzU+zuPcdub33G8uJyOLQN4+YbeCj0iIo2AW474JCcnk52dzdq1a6mqqmLGjBmkpKSQlJRU6/E7duwgOTmZcePGnbHvX//6F8uWLeP555+nS5cufPLJJ8yYMYO2bdty/vnnO/utiIeprHLw1Gd7eWvbUQAu7RLKgyPjCPZzy78qIiJST2434mOz2UhNTWXatGm0aNGCsLAwEhMTWbZsGTab7Yzjy8vL2b17Nz179qz1fKdOnWLq1KlER0djsVgYMmQI0dHRbN261dlvRTxMQWk5U9/ZURN6bruoIynX9FDoERFpREz5jW6328nJyal1n81mo6KigpiYmJpt0dHR2O129u/fT3x8/GnHp6enU1lZyZIlS9iyZQshISFMmDCBKVOm4OXlxbRp0047Pisriz179tCjR4+z1ueM+1Z/OqfuiXWuc+1z2rEiElfsIqeojCBfbxaM0P08v0SfZ9dQn11DfXYdd+i1KcFn+/btTJ48udZ906dPByAwMLBmW0BAAECt9/kUFRUxYMAAJk2axBNPPEFaWhpTp07Fy8uLKVOmnHbsvn37uO222xgzZgz9+/ev9fuHhgbh7cQVtcPCQpx2bvmv+vT53S2HmbN8B+WVDrqEB/HPyX3p2ko/p7rQ59k11GfXUJ9dx8xeWwzDMEz77rXYtWsX48aNY+vWrQQFVU8OV1xcTN++fVmxYgVxcXG/eo4XXniB1atXs2zZspptH3/8MbNnz2b8+PHce++9WM4SN3Nzi5w24hMWFkJ+fhHu1fHGpT59rqxysPizffznx/l5Lu0SSvIo3c9TF/o8u4b67Brqs+s4s9fh4XULU273G75z585YrVYyMzPp1asXUH15ymq1EhUVdcbx69evJy8vj+uvv75mW3l5Of7+/jVfL126lBdeeIEHH3yQhISEX63BmR98w3Du+aXar/W5sLScOe+nseVQ9fw8Uy7syG0/zs+jn0/d6fPsGuqza6jPrmNmr93u5uaAgABGjBhBSkoKBQUFFBQUkJKSwujRo08LMz8xDINHHnmEDRs2YBgG27Zt49VXX2XixIkAvPzyy7z88su8/vrrdQo90vhl5BQz+d/b2HLoJEG+3jw2pju3D4rSpIQiIk2A213qgupLW4sWLeLjjz+moqKCoUOH8sADD9Tc9zNq1CgSEhK44447AHjzzTd5+eWXycnJITw8nD/84Q/ceOONGIZB//79sdls+Pr6nvY9br/99prX/1xubpFT3pPFUj0Ml5enoVRn+rU+r0k7zkPrdlNW6aBjywAeG9udLmFab6u+9Hl2DfXZNdRn13FmryMi6napyy2Dj5kUfDzb2fpc6TBY+sU+/r35MAAXd27JQyPjCfF3u6u9HkGfZ9dQn11DfXYddwg++q0vjd5JWwX3rUrjmwMnAPjDwA7cfnEU3l66tCUi0tQo+EijlplbQuKKnRw5acffx4uk4bFcGRthdlkiImISBR9ptD7ancuCNRnYKhy0be5PytjudIsINrssERExkYKPNDoOh8EzX+7jpY2HAOjfsQULR8fTIsBqcmUiImI2BR9pVIrLKrn31c18lH4cgBv6tuPuy7rgo/t5REQEBR9pRPYXlJK4YicHCmz4elu476oYRnZvbXZZIiLiRhR8pFH4am8B961Ko6S8ijbN/VmUEE98a627IyIip1PwEY9mGAb/2nSIZ77cjwH0ateM52/pj1dZhebjEBGRMyj4iMeyV1Tx4NrdfJiRC8C48yOZNbQrrUL8ySurMLk6ERFxRwo+4pGyT9lJfG8nu3NL8PayMHNINBN6tUXLbYmIyC9R8BGPs+XQCWanpnHCVkHLACuLxnSnd/vmZpclIiIeQMFHPIZhGLyzPZvHP8miymEQ1yqYx8Z2J7KZv9mliYiIh1DwEY9QUeXg0Y8yeW/HMQCujovg/qti8Ld6m1yZiIh4EgUfcXsFpeXcu3IX3x05hQW4+7LO3NSvPRbd0CMiIvWk4CNuLSOnmHtW7CSnqIxgP28eGhXPoM6hZpclIiIeSsFH3Na69OM8uHY3ZZUOOrYM4PFrehAVGmh2WSIi4sEUfMTtOAyDf3y5n1c2VS8yenHnljw0Mp4Qf31cRUTkt9G/JOJWissqeWB1Ol/uLQBgcv/23HVJZ7y1yKiIiDQABR9xG4cKbdzz3k72FZTi5+PFfVd1Y0S8FhkVEZGGo+AjbuGb/YXMeT+NorJKWgX78tjYHnSP1CKjIiLSsBR8xFSGYfDWtqM89WkWVQac16YZj47tTniQr9mliYhII6TgI6Ypr6yelHDFD9WTEo7u0Zo5V3bD18fL5MpERKSxUvARU+SXlDNr5S6+P3oKLwtMv7wLv+/TTpMSioiIUyn4iMul5xRxz3s7OV5cToifDwtHx3FhlCYlFBER51PwEZf6MCOXBWsyKKt00OnHSQk7aVJCERFxEQUfcQmHYfDc1wd4aeNBoHpSwodHxRPsp4+giIi4jv7VEacrLa8i6YN0Ps3MB2BSv/ZMvVSTEoqIiOsp+IhTHT1p5573dpKZV4Kvt4X7rophZHdNSigiIuZQ8BGn2Xr4BPeuTOOErYKwIF8eG9Od89o2M7ssERFpwhR8xCmWf5/Noo8yqXIYxLUKJuWaHrQO8TO7LBERaeIUfKRBVToMnvo0i7e2HQVgWGwE866Owd/qbXJlIiIiCj7SgE7ZK5j7fhrfHDgBwB2DOvHHgR01KaGIiLgNBR9pEPvzS7lnxU4OFtoIsHqxYEQcV3QLN7ssERGR0yj4yG/29b4C7luVRnFZFZEhfjx+TQ9iWgWbXZaIiMgZFHzknBmGwX+2HmHxZ3txGHBBu2YsGtOd0ECtrC4iIu5JwUfOSUWVg0UfZbJiR/XK6gk9WjNbK6uLiIibU/CReissLefelbvYdkQrq4uIiGdR8JF6ycwr4Z7lP3D0VBlBvt48PDqeQZ21srqIiHgGt7wuUVpaypw5cxg4cCB9+/Zl1qxZlJSUnPX4pKQkevbsSe/evWv+vPXWWwCUlZWRnJzMoEGD6N27N7/73e/YsGGDq95Ko/JFVj63vvEdR0+V0a65Py/dcIFCj4iIeBS3DD7JyclkZ2ezdu1a1q1bR3Z2NikpKWc9fseOHSQnJ7Nt27aaPxMnTgTgySefZPv27bz33nts2bKFsWPHctddd/1ikJLTGYbBa98e4p73dlJaUUXfDs155cbedAkLMrs0ERGRenG7S102m43U1FReffVVWrRoAUBiYiKTJ09m1qxZBAQEnHZ8eXk5u3fvpmfPnrWeb+bMmZSXlxMQEEBpaSknTpwgJCQEq9Xq7LfSKJRXOli4fg+rduYAMO78SGYO6YrV2y0zs4iIyC8yJfjY7XZycnJq3Wez2aioqCAmJqZmW3R0NHa7nf379xMfH3/a8enp6VRWVrJkyRK2bNlCSEgIEyZMYMqUKXh5eeHt7U1AQABvvfUWSUlJ+Pj4kJKSgq/v2R+5dsY9uj+d05Pu/y0oKSdxxS6+P1p9E/Nfr4hmYu+2bn0Tsyf22ROpz66hPruG+uw67tBrU4LP9u3bmTx5cq37pk+fDkBgYGDNtp9GeWq7PFVUVMSAAQOYNGkSTzzxBGlpaUydOhUvLy+mTJlSc9w111zD+PHjWbduHYmJiURERNC3b98zzhcaGoS3E0czwsJCnHbuhpSWfYop/9nOkRM2Qvx9WHpDHy6LiTC7rDrzlD57OvXZNdRn11CfXcfMXlsMwzBM++612LVrF+PGjWPr1q0EBVXfQ1JcXEzfvn1ZsWIFcXFxv3qOF154gdWrV7Ns2bJa999222106tSJ+++//4x9ublFThvxCQsLIT+/CPfq+Jk+z8rn/vfTKa2ookMLf54c15OosMBff6Eb8KQ+ezL12TXUZ9dQn13Hmb0OD69bmHK7e3w6d+6M1WolMzOTXr16AZCVlYXVaiUqKuqM49evX09eXh7XX399zbby8nL8/f0BmDFjBhdccAG33HLLaft/un+oNs784BuGc8//WxiGwb83H+bpz/dhAP06NOdvCd1pHmB125rPxp373Jioz66hPruG+uw6Zvba7e5QDQgIYMSIEaSkpFBQUEBBQQEpKSmMHj26Jsz8nGEYPPLII2zYsAHDMNi2bRuvvvpqzVNdvXv35vnnnycjI4PKykrefvttduzYwZgxY1z91txaeaWDB9fuZsmPoWf8+W14esJ5NA/QTeAiItJ4uN2ID1TPy7No0SISEhKoqKhg6NChPPDAAzX7R40aRUJCAnfccQfDhg1jzpw5zJ8/n5ycHMLDw7n77rsZO3YsAJMnT6asrIw777yToqIi4uLieOWVV+jYsaNZb8/t/O9MzH8dHM3v3PwmZhERkXPhdvf4mC03t8gp57VYqq8/5uW51zXk/52JeeHoeC724EkJ3bXPjY367Brqs2uoz67jzF5HRHjoPT7iOl/tLeC+VWmUlFfRrnn1TcydPeQmZhERkXOh4NMEGYbBf7YeYfFne3EY0Kd9cxaN6U4L3c8jIiKNnIJPE1NZ5WDRR5m8t+MYAGN7RnLvlZqJWUREmgYFnybkpK2C2e+nsfngCSzA9Mu7cEPfdrqJWUREmgwFnybiQEEpf31vJwcLbQRavUkeFcdl0WFmlyUiIuJSCj5NwLcHC5mdmsYpeyWRIX48Ma4H3SKCzS5LRETE5RR8Grnl32ez6KNMqhwG57UJ4bGxPQgLOvsCrSIiIo2Zgk8jVeUwWPL5Xt7YcgSAq+MieODqWPx8dBOziIg0XQo+jVBJeSX3r0rny70FANx+cSduvbCjbmIWEZEmT8Gnkck+Zeevy3eSmVeCn48XScNjGRYbYXZZIiIibkHBpxH5IfsU97y3k4LSCkIDrTxxTQ96tGlmdlkiIiJuQ8GnkViXfpwFazIorzLoFhHEE9f0ILLZmavZi4iINGUKPh7OMAxe3HiQ574+AMAlXUJ5aFQcQb760YqIiPwv/evowcoqHTy0bjdr0o4DcEPfdky7rAveXrqJWUREpDYKPh6qsLScxBW7+P7oKby9LMwa2pXx57cxuywRERG3puDjgfbml/CX5Ts5etJOsJ83ixK6M6BTS7PLEhERcXsKPh5m4/4CZqemUVJeRfsW/jx5TU+iwgLNLktERMQjKPh4kHe+O0rKx5lUGdC7XTMeHdODFoFWs8sSERHxGAo+HuB/l58Y1aM1c6/shq+WnxAREakXBR83V1pexf2r0vjix+Un7rokilsGdNDyEyIiIudAwceN5RSV8dflP7A7V8tPiIiINAQFHzeVnlPEX9/bSW5xOaGBVh6/pgc9tfyEiIjIb6Lg44Y+y8zj/lXp2CsddAkL5MlxPWnbXMtPiIiI/FYKPm7EMAxe33KEJZ/txQAu7NSSRxLiCfbTj0lERKQh6F9UN1FZ5eCxj7NY9n02ABN6tSFxSFd8tPyEiIhIg1HwcQPFZZXMTt3FNwdOYAFmDO7C7/u005NbIiIiDUzBx2RHT9r5y/If2JtfSoDVi4dGxXNZdJjZZYmIiDRKCj4m2pl9ir++t5OC0goign158pqexLYONrssERGRRkvBxyQf785l3gcZlFU66BYRxJPjetI6xM/sskRERBo1BR8XMwyD1749zJLP9wFwSZdQHhoVR5CvfhQiIiLOpn9tXaiiysHCD/ew/PtjAEzs3ZYZg6P15JaIiIiLKPi4SHFZJX955Vu+2JOHlwX+Mjia6/u0M7ssERGRJkXBxwVO2ir401vb9eSWiIiIyRR8XGDL4ZPszS+ldTM/Hh/bg9hWenJLRETEDAo+LnBZl1AeG9udK85rC/ZyDMPsikRERJomL7MLaAp8vL24ols44cF6XF1ERMRMCj4iIiLSZCj4iIiISJOh4CMiIiJNhlsGn9LSUubMmcPAgQPp27cvs2bNoqSk5KzHJyUl0bNnT3r37l3z56233jrjuK+++or4+HgOHz7szPJFRETETbnlU13JyclkZ2ezdu1aqqqqmDFjBikpKSQlJdV6/I4dO0hOTmbcuHFnPWdubi733nsvDofDWWWLiIiIm3O7ER+bzUZqairTpk2jRYsWhIWFkZiYyLJly7DZbGccX15ezu7du+nZs+dZz+lwOEhMTOS6665zZukiIiLi5kwZ8bHb7eTk5NS6z2azUVFRQUxMTM226Oho7HY7+/fvJz4+/rTj09PTqaysZMmSJWzZsoWQkBAmTJjAlClT8PKqznXPPPMMYWFhTJgwgWeeeeZX67M4Yemsn87pjHPLf6nPrqE+u4b67Brqs+u4Q69NCT7bt29n8uTJte6bPn06AIGBgTXbAgICAGq9z6eoqIgBAwYwadIknnjiCdLS0pg6dSpeXl5MmTKFTZs2sXLlSpYtW8aJEyd+tbbQ0CC8vZ03EBYWFuK0c8t/qc+uoT67hvrsGuqz65jZa1OCz8CBA8nIyKh1365du1i8eDE2m42goCCAmktcwcFnLvUwaNAgBg0aVPP1+eefz80338zq1asZP348s2fP5sknnyQ4OLhOwaegoMRpIz5hYSHk5xdp5mYnUp9dQ312DfXZNdRn13Fmr8PD6xam3O7m5s6dO2O1WsnMzKRXr14AZGVlYbVaiYqKOuP49evXk5eXx/XXX1+zrby8HH9/f7744gvy8/O59dZbAWpubB4zZgx33HEHf/rTn2qtwZkffMNw7vmlmvrsGuqza6jPrqE+u46ZvXa7m5sDAgIYMWIEKSkpFBQUUFBQQEpKCqNHj8bf3/+M4w3D4JFHHmHDhg0YhsG2bdt49dVXmThxImPHjmX79u1s3ryZzZs3s3LlSgBWrlx51tAjIiIijZfbBR+onpcnKiqKhIQEhg8fTvv27Zk3b17N/lGjRvHss88CMGzYMObMmcP8+fPp3bs3M2fO5O6772bs2LFmlS8iIiJuymIYGtj7udzcIqec12Kpvv6Yl6dryM6kPruG+uwa6rNrqM+u48xeR0TU7R4ftxzxEREREXEGjfiIiIhIk6ERHxEREWkyFHxERESkyVDwERERkSZDwUdERESaDAUfE9hsNiZOnMiyZcvMLqVRyc/P56677qJfv34MHDiQhx9+mMrKSrPLarQKCgoYNmwY33zzjdmlNErp6en84Q9/YMCAAQwaNIhZs2ZRUFBgdlmNzoYNG7juuuvo06cPgwYNIjk5GbvdbnZZjVZVVRWTJk1i9uzZptWg4ONie/bs4cYbb+S7774zu5RGZ8aMGQQGBvLFF1/wzjvvsGHDBl555RWzy2qUtmzZwsSJEzl48KDZpTRKdrudKVOm0Lt3b7788kvef/99Tpw4wdy5c80urVEpKCjg9ttv5/e//z2bN29m+fLlbNq0iX/+859ml9Zo/f3vf2fz5s2m1qDg40IbNmzg5ptvZty4cbRt29bschqVAwcOsGnTJmbOnElAQAAdOnTgrrvu4vXXXze7tEZn+fLlJCYm8pe//MXsUhqto0ePEhcXx9SpU/H19aVly5ZMnDiRb7/91uzSGpXQ0FC+/vprxo8fj8Vi4cSJE5SVlREaGmp2aY3Shg0bWLduHVdddZWpdbjdIqWezG63k5OTU+u+iIgI4uLi+OSTT/Dz8+Pll192cXWN2549e2jRogWtW7eu2RYdHc3Ro0c5deoUzZo1M7G6xuWSSy4hISEBHx8fhR8n6dKlCy+88MJp29auXUuPHj1MqqjxCg4OBuDyyy8nJyeHfv36MX78eJOranzy8/O57777eOaZZ0wfiVfwaUDbt29n8uTJte5bunQpV155pYsrajpKSkoICAg4bdtPX5eWlir4NKCIiAizS2hSDMPgqaee4pNPPuHf//632eU0WuvWrePkyZMkJiYybdq0M4KnnDuHw8HMmTP5wx/+QFxcnNnlKPg0pIEDB5KRkWF2GU1SYGAgNpvttG0/fR0UFGRGSSK/WXFxMXPmzGHnzp38+9//JjY21uySGi1/f3/8/f2ZOXMm1113HSdPnqR58+Zml9UoPPfcc/j6+jJp0iSzSwEUfKSR6NatGydOnCAvL4/w8HAAsrKyiIyMJCSkbgvXibiTgwcPctttt9G2bVveeecd3XfiBFu3bmXu3LmsXLkSX19fAMrLy7FarWeMIMu5W7FiBcePH6dfv34ANU/NrV+/3pQbnXVzszQKUVFR9O3bl4ULF1JcXMyhQ4d45plnuPbaa80uTaTeTp48yc0330yfPn148cUXFXqcJDY2FrvdzuOPP055eTlHjhxh0aJFXHvttTVBSH67NWvWsHXrVjZv3szmzZsZPXo0o0ePNu3pLo34SKOxZMkSHnzwQYYOHYqXlxfXXHMNd911l9llidTbsmXLOHr0KB988AFr1qw5bd+2bdtMqqrxCQoK4oUXXmDhwoUMGjSIkJAQEhISmDp1qtmliRNpdXYRERFpMnSpS0RERJoMBR8RERFpMhR8REREpMlQ8BEREZEmQ8FHREREmgwFHxEREWkyFHxERESkyVDwERERkSZDwUdERESaDC1ZISIe6e9//ztr167l0KFDBAYGMmzYMO6//36sVqvZpYmIG1PwERGPYxgGDoeDBQsW0Lp1azIzM7n33nuJjY3lhhtuMLs8EXFjCj4i4nEsFgvTpk2r+bpdu3YMGjSIvXv3AnDo0CH27NnDkCFDzun8U6dOZdOmTVx00UUsWbKkQWoWEfeg4CMiHufIkSO8+OKLbNq0iZycHCorKykrK+Ovf/0rAJ9//jk2m+2cg8+kSZOYMGEC7733XgNWLSLuQDc3i4hHKSgo4LrrrqOwsJDZs2fzxhtv8Oabb+Ll5UV8fDybNm3iqaee4u233+aaa67BbrfX+3tceOGFBAUFOaF6ETGbRnxExKN8/vnnVFVV8cQTT2CxWAB4/fXXqaioID4+ntDQULp168bjjz9OmzZtTK5WRNyNgo+IeJTmzZtTXFzMRx99RNeuXfnkk0947rnnaN26NaGhoQAcO3as1tAzfvx4ysvLz9j+4osv0rp1a6fXLiLmU/AREY8yePBgJkyYwKxZs/D392fMmDGMGDGCo0ePAtWh52whZtmyZa4sVUTckIKPiHgUi8XCgw8+yIMPPljr/sOHD9OqVSsXVyUinkI3N4tIo9KtWzcOHDhAQkICWVlZ53SOW2+9lenTp/PZZ59x2WWX8f333zdwlSJiFothGIbZRYiIiIi4gkZ8REREpMlQ8BEREZEmQ8FHREREmgwFHxEREWkyFHxERESkyVDwERERkSZDwUdERESaDAUfERERaTIUfERERKTJUPARERGRJkPBR0RERJoMBR8RERFpMv4fUs0a/mCwodsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create a grid for a\n",
    "a_tm1 = np.linspace(-1, 4)\n",
    "\n",
    "# Compute income components at every a\n",
    "cap_income_t = a_tm1 * (Rfree - 1) / Γ\n",
    "lab_income_t = 1\n",
    "\n",
    "# and market resources\n",
    "m_t = a_tm1 * Rfree / Γ + 1\n",
    "\n",
    "# Consumption\n",
    "c_t = PFsavrate.solution[0].cFunc(m_t)\n",
    "\n",
    "# and finally the saving rate\n",
    "ς_t = sav_rate_t = (cap_income_t + lab_income_t - c_t) / (cap_income_t + lab_income_t)\n",
    "\n",
    "# And now the plot\n",
    "plt.plot(a_tm1, sav_rate_t)\n",
    "plt.xlabel(r\"$a_{t-1}$\")\n",
    "plt.ylabel(\"Saving Rate\")"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-",
   "formats": "py:percent,ipynb",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}