A Simple Ramsey Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### [**NICOLAS CACHANOSKY**](http://www.ncachanosky.com) | Department of Economics | Metropolitan State University of Denver | ncachano@msudenver.edu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A SIMPLE RAMSEY MODEL\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This note illustrates how to code a simple Ramsey model in Python. The purpose of the note is to walk through Python applications, not to offer a detailed discussion of Ramsey model or to show best coding practices. The note also assumes familiarity with Ramsey model and a beginner experience with Python.\n",
    "\n",
    "For a more complete and detailed discussion of Python applications see the material in [Quant Econ](https://quantecon.org/).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TABLE OF CONTENTS\n",
    "1. [Description of the model](#1.-DESCRIPTION-OF-THE-MODEL)\n",
    "2. [The respresentative household](#2.-THE-REPRESENTATIVE-HOUSEHOLD)\n",
    "3. [The Representative Firm](#3.-THE-REPRESENTATIVE-FIRM)\n",
    "4. [Equilibrium](#4.-EQUILIBRIUM)\n",
    "5. [The code](#5.-THE-CODE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. DESCRIPTION OF THE MODEL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A known issue with the Solow model is that the savings rate is exogenous (and constant). This means the ratio of consumption to income $Y/C$ is also constant. While the Solow model can be useful to highlight some features of long-run growth, this assumption also limits the applicability of this model. Yet, the Solow model can be interpreted as a specIal case of Ramsey model (or Ramsey-Kass-Koopmans model). Ramsey model can be considered as the building block of macro models with a household that optimizes its utility and, therefore, endogenously sets the savings rate.\n",
    "\n",
    "This note presents a \"simple\" Ramsey model (it offers a simplified discussion, does not discuss shocks, and does not add other features). For a more complete and detailed discussion see an advanced macroeconomics textbook. \n",
    "\n",
    "The model is populated by a representative household that maximizes its utility and by a representative firm that maximizes its profits. The model assumes a competitive market (no monopoly power, agents are price takers, etc.) The household provides labor services to the firm in exchange of a wage. The wage can be used to consume goods or save and accumulate assets. The model assumes a closed economy with no government. Finally, the model is presented in continuous (rather than discrete) time. A \"dot\" on top of a variable denotes its instantaneous time-change $(\\dot x (t) = \\frac{\\partial x}{\\partial t})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. THE REPRESENTATIVE HOUSEHOLD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume a large number of identicial infinitely-lived households. Because these households are identical, we can use  a **representative** household to populate the model. The assumptions that this households are infinitely lived captures a continuum of generations. The representative household has $L(t)$ individuals with a population growth rate of $n$. Therefore, $L(t) = L(0) \\cdot e^{nt}$\n",
    "\n",
    "## 2.1. THE HOUSEHOLD UTILITY\n",
    "\n",
    "The household utility $(U)$ is the present value of the utility (felicity function) of all household individuals $(u(c))$.\n",
    "\n",
    "\\begin{equation}\n",
    "    U = \\int_{0}^{\\infty} u[c(t)] \\cdot e^{nt} \\cdot \\underbrace{e^{-\\rho t}}_{L(t)} \\cdot dt = \\int_{0}^{\\infty} u[c(t)] \\cdot e^{(n-\\rho)t} \\cdot dt\n",
    "\\end{equation}\n",
    "\n",
    "where $\\rho > 0$ is the time preference of the household individual and where $n -\\rho >0$. Also, $u(c)$ is concave $(u^{'}(c) > 0, u^{''}(c) < 0)$ and satisfies the Inada conditions $(u^{'}(c) \\to \\infty \\text{ as } c \\to 0 \\text{, and } u^{'}(c) \\to 0 \\text{ as } c \\to \\infty$.\n",
    "\n",
    "The parameter $\\rho$ can be interpreted as the time preference of the household or as capturing a discount of the future generation by the present generations. If $\\rho = 1$, then present generation assigns equal \"value\" to the next generatio (i.e. parents with respect to their kids). If $\\rho < 0$, then the present generation is altruistic and values the utility of the next generation more than their own.\n",
    "\n",
    "## 2.2. THE HOUSEHOLD BUDGET CONSTRAINT\n",
    "\n",
    "The household budget constrained is composed of assets $(\\Lambda)$ and wages $(w)$. The non-consumed wage becames an increase in assets. At some period in time, assets can be negative, which means that in that period the household has a net debt. Assets have a rate of return of $r$. Therefore, the stock of assets increase by the return of the assets plus wages, and decrease by the amount of consumption. Since population growths at rate $n$, the change in household's assets can be represented as follows.\n",
    "\n",
    "\\begin{equation}\n",
    "    \\frac{d\\Lambda}{dt} = r(t) \\cdot \\Lambda(t) + w(t)L(t) - C(t)\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Dividing by $L(t)$ we get the change in assets per capita $(a)$ (and dropping the $t$ variable for notation simplicity):\n",
    "\n",
    "\\begin{align}\n",
    "    \\dot{\\left(\\frac{A}{L} \\right)} &= \\frac{\\dot{A}L - A \\dot{L}}{L^2}                \\\\\n",
    "    \\dot{a} &= \\frac{\\dot{A}}{L} - \\frac{A}{L}n                                        \\\\\n",
    "    \\dot{a} &= ra + w - c - na                                                         \\\\\n",
    "    \\dot{a} &= (r - n)a + w - c\n",
    "\\end{align}\n",
    "\n",
    "To avoid the situation where the household borrows money indifenetely, the model also imposses a \"no Ponzi game condition\". Intuitively, at the end of the household life $(T)$, assets have to be at least zero $a(T) \\geq 0$. A household may have debt at any point in time, but at the end of its lifetime all debts are paid. The present value of the assets have to be asymptotically nonnegative.\n",
    "\n",
    "\\begin{equation}\n",
    "    \\lim\\limits_{t \\to \\infty} a(t) e^{-\\tilde{r}(t)} \\geq 0\n",
    "\\end{equation}\n",
    "\n",
    "where $\\tilde{r}(t) = \\int_{0}^{t} (r(\\tau) - n) d\\tau$. This expression means that any moment $\\tau \\in (0, \\infty)$, the present value of $a$ cannot be negative. Note that $\\tilde{r}$ is taking into consideration the growht rate of the household as well $n$. In a more simple notation, this condition can be written as $\\lim\\limits_{t \\to \\infty} a(t) \\geq 0$.\n",
    "\n",
    "## 2.3. THE HOUSEHOLD OPTIMIZATION PROBLEM\n",
    "\n",
    "The household is  facing a dynamic optimization problem. It needs to set the path of consumption that maximizes the present value of its life utility. The optimization problem can be described with the following equations (where $n - \\rho > 0$).\n",
    "\n",
    "\\begin{align}\n",
    "    \\max_{c(t)} U &= \\int_{t=0}^{\\infty} u[c(t)] \\cdot e^{-(n-\\rho)t}    \\\\\n",
    "    \\text{subject to:}                                                   \\\\\n",
    "    \\dot{a} &= (r(t) - n)a(t) + w(t) - c(t)                              \\\\\n",
    "    a(0) &= a_{0}                                                        \\\\\n",
    "    \\lim\\limits_{t \\to \\infty}a(t) &\\geq 0                               \\\\\n",
    "    \\text{transversality condition:}                                     \\\\\n",
    "    \\lim\\limits_{t \\to \\infty} \\mu(t) \\cdot a(t) = 0\n",
    "\\end{align}\n",
    "\n",
    "The optimization problem is given an initial asset value of $a_0$ and a terminal condition for which assets cannot be negative. Let $\\mathscr{H}$ be the present-value Hamiltonian:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathscr{H} = e^{-(\\rho - n)t} \\cdot u[c(t)] + \\mu(t) \\cdot \\left[r(t) - n)a(t) + w(t) - c(t) \\right] \n",
    "\\end{equation}\n",
    "\n",
    "where $\\mu(t)$ is the present-value shadow price. Since the household faces an optimization problem at each moment $t$, the shadow price may vary over time as well. The \"transversality condition\" states that in the limit, either the shadow price or the assets have to be zero. \n",
    "\n",
    "The first-order-conditions (FOC) of the hamiltonian are:\n",
    "\n",
    "\\begin{align}\n",
    "    \\frac{\\partial \\mathscr{H}}{\\partial c} &= e^{-(\\rho - n)t} \\cdot u'_{c} - \\mu(t) = 0       \\\\\n",
    "    \\frac{\\partial\\mathscr{H}}{\\partial\\mu} &= (r(t) - n)a(t) + w(t) - c(t) = \\dot{a}           \\\\\n",
    "    \\frac{\\partial \\mathscr{H}}{\\partial a} &= - (r(t) - n)\\cdot \\mu (t) = \\dot{\\mu}(t)\n",
    "\\end{align}\n",
    "\n",
    "Now we proceed to use the first equation to get $\\dot{\\mu}(t)$.\n",
    "\n",
    "\\begin{align}\n",
    "    \\mu(t)       &= e^{-(\\rho - n)t} \\cdot u'_{c} \\\\\n",
    "    \\dot{\\mu}(t) &= -(\\rho - n) \\cdot e^{-(\\rho - n)t} \\cdot  u'_{c} + e^{-(\\rho - n)t} \\cdot u''_c \\cdot \\dot{c}\n",
    "\\end{align}    \n",
    "\n",
    "Make now two substitutions: (1) The first FOC in the third FOC, and (2) this last derivative in the third FOC as well.\n",
    "\n",
    "\\begin{align}\n",
    "    -[r(t) -n] \\cdot \\mu(t) &= -(\\rho - n) \\cdot e^{-(\\rho - n)t} \\cdot u'(c) + e^{-(\\rho - n)t} \\cdot u''(c) \\cdot \\dot{c}        \\\\[10pt]\n",
    "    -[r(t) -n] \\cdot e^{-(\\rho - n)t} \\cdot u'(c) &= e^{-(\\rho - n)t} \\left[-(\\rho - n) \\cdot u'(c) + u''(c) \\cdot \\dot{c} \\right] \\\\[10pt]\n",
    "    -[r(t) -n] \\cdot u'(c) &= -(\\rho - n) \\cdot u'(c) + u''(c) \\cdot \\dot{c}                                                       \\\\[10pt]\n",
    "    r(t) &= \\rho - \\frac{u''(c)\\cdot c}{u'(c)} \\cdot \\frac{\\dot{c}}{c}\n",
    "\\end{align}\n",
    "    \n",
    "The last equation is the **Euler equation**, which defines the path for $c$ given the rate of return $r(t)$. Recall that $u''(c)<0$ and therefore $\\frac{u''(c)c}{u'(c)}>0$. The right-hand side of the equation can be interpreted as the \"rate of return of consumption.\" Consumers would like to smooth consumption in time. Therefore, the rate of return $(r)$ has to be equal to the time preference $(\\rho)$ minus the change of consumption over time. Or, put it differently, if $r = \\rho$, then consumption will be flat. There is no opportunity to increase utility by moving consumption over time (extend a loan and postpone consumption or take a debt and advance consumption). Consumers will be willing to postpone consumption if $r$ is sufficiently above $\\rho$. The required compensation to postpone consumption is the second term in the righ-hand side of the equation.\n",
    "\n",
    "\\begin{align}\n",
    "    r &= \\rho \\to \\frac{\\dot{c}}{c} = 0 \\\\\n",
    "    r &> \\rho \\to \\frac{\\dot{c}}{c} > 0 \\text{    (saving today: postpone present consumption)} \\\\\n",
    "    r &< \\rho \\to \\frac{\\dot{c}}{c} < 0 \\text{    (consuming today: advance future consumption)}\n",
    "\\end{align}\n",
    "\n",
    "The Euler equation also includes the elasticity of marginal utility: $\\frac{u''(c)\\cdot c}{u'(c)}$. For $r$ and $\\frac{\\dot{c}}{c}$ to be constant, this elasticity of marginal utility has to be constant. A common practice is to assume a CRRA (constant relative risk aversion) utility function that depicts constant intertemporal elasticity of substitution (CIES).\n",
    "\n",
    "\\begin{equation}\n",
    "    u(c) = \\frac{c^{1-\\theta}-1}{1-\\theta}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\theta > 0$ is the coefficient of relative risk aversion. The $1 - \\theta$ denominator assures the function will keep its desired properties of the function $(u'>0, u''<0)$ even if $\\theta > 1$ (see below). The inclution of the term $-1/(1-\\theta)$ is for convenience; when $\\theta \\to 1$ then the CRRA function becomes $(u(c) = log(c))$. Consider now the first and second derivative of the CRRA function and the elasticity of substitution $(\\sigma)$.\n",
    "\n",
    "\\begin{align}\n",
    "    u'(c)  &= c^{-\\theta} > 0                      \\\\\n",
    "    u''(c) &= -\\theta \\cdot c^{-(\\theta + 1)} < 0  \\\\\n",
    "    \\sigma &= \\left[ \\frac{u''(c) \\cdot c}{u'(c)} \\right]^{-1} = \\left[\\frac{\\theta \\cdot c^{-(\\theta + 1)} \\cdot c}{c^{-\\theta}} \\right]^{-1} = \\frac{1}{\\theta}\n",
    "\\end{align}\n",
    "\n",
    "Then, the Euler equation from the household maximization problem becomes\n",
    "\n",
    "\\begin{align}\n",
    "    r(t) &= \\rho - \\frac{u''(c)\\cdot c}{u'(c)} \\cdot \\frac{\\dot{c}}{c} \\\\\n",
    "    r(t) &= \\rho - \\theta \\cdot \\frac{\\dot{c}}{c}                      \\\\\n",
    "\\end{align}\n",
    "\n",
    "We can rewrite the Euler equation to more conveniently depict the consumption path that maximizes the household utility.\n",
    "\n",
    "\\begin{equation}\n",
    "    \\frac{\\dot{c}}{c} = \\frac{1}{\\theta} \\cdot \\left(r(t) - \\rho \\right)\n",
    "\\end{equation}\n",
    "\n",
    "This is one of the two equations that describes Ramsey model. The second equation is derived from the firms profit maximization problem (next section).\n",
    "\n",
    "---\n",
    "\n",
    "The following code plots the utility function for different levels of $\\theta$. The function also allows for $\\theta = 1$, in such case the utility function becomes the log of consumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGKCAYAAADQeD9lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXxcZ33v8c8zm6x1vMR2YsexnFUTEAQSkrAHEUKoWW8LBVqqQoH6opT93qq3FASvC3VvCyVtdSu4QOv0lnApBUoj1iAgARKSQEKGZOTEsWXHlu3Em6zNWmae+8dzjudoNJsWSzP29/16ndc5c86jozNOwN88z3N+j7HWIiIiIiKLJ7TcDyAiIiJytlHAEhEREVlkClgiIiIii0wBS0RERGSRKWCJiIiILLKKCljGmB8v9zOIiIiILFRFBSwRERGRs4ECloiIiMgiU8ASERERWWQKWCIiIiKLTAFLREREZJEtKGAZY2qNMT8xxoQLXI8ZY+4yxkQW8ntEREREqslCe7DeAXzdWpvOd9FaOwn8EPjdBf4eERERkapRVsAyxtxrjGn2jjcaYx7wLv0e8B+BdhuMMf9ujHnQGNNvjLkW+KbXTkREROScUHLozhhjgIuAvd6pZwFJY0wMuNhaO+C1iwDfAf7cWnuHMaYOCANjwPPOwLOLiIiIVKRy5kZdCuyx1lrv87OAJHAecCLQ7vVAylp7B4C1dsy/YIyZNMY0WmuHF+exRURERCpXOUOErbhA5bsGeBgYB1YEzl8F3FvgHjXAqfk8oIiIiEi1KSdgrcaFKYwxCWArkLTWHgfCxhg/ZB0CnuH/kDFmrbdfAzxtrZ1azAdfasaYsDHmVmPMI8aYpDHm4nne52ZjzE5jzC5jTGee65uMMT8yxqS83/W+nOsD3u9/KDAXTkRERCpIOQHre8DLjTFfBd4IHLXWHvaufR94kXf8z8B6LxQ8BDzfO/8y4NuL98jL5s+A3dbaZwB/B7xnrjfwyll0A68CrgTeYoy5MqfZNPAha20CuB7oyNPmZdbaq6y118z1GUREROTMKzkHy1r7JG7ele8TgeN/AD4I3GmtHQFem+cWb8WFk6pljKkH3mCtvdo7tQfXkzdX1wK7rLW7vft+BXgd8KjfwFp7EDjoHQ8bY1LAxmAbERERqWwLKgBqrX3QG84K56uF5b1p+E1r7c6F/B5fc2fvh4ArFuNeATsHtm/9dIk2NwKbvJ45cMOmdwYbGGPuBhrz/OyHrbV+243Ak4Fr+4HrCv1SrzTGc4BfBE5b4PvGGAt8zlr7+RLPLiIiIktswRXWrbVfKnJtErhtob+jAlwFfNRa2wNgjPkCbqL/adbaF5dxH5PnnM1zDmNMA/DvwPuttScDl15orR00xqwDfmCM6bfW3lXOlxAREZGlUVVL2JTR03SmrMINC/r1vm4CPhlsUGYP1n5gU+DahcBg7g8YY6K4cPWv1tqvB69Zawe9/VPGmG/ghh0VsERERCpIVQWsZfQYbsL5bcAHgF5r7Z5ggzJ7sO4HLjPGbAEOAG/GzVE7zSvs+kVcTbHP5FyrB0Le3Kx6XNALzokTERGRCrDQtQjPFbcDzzXG7MJN+P/gfG5irZ0GbsG9mZkCvmqtfQTAGPNtY8wG4IXA24A2rxTDQ8aY3/JusR74qTHm18B9uKD33YV8MREREVl8JlugffkZY35srb1huZ9DREREZCHUgyUiIiKyyBSwRERERBaZApaIiIjIIlPAEhEREVlkClgiIiIii0wBS0RERGSRKWCJiIiILDIFLBEREZFFpoBVJmNM2BhzqzHmEWNM0hhz8Tzvc7MxZqcxZpcxpjPP9SsCFdwfMsacNMa8P3B9wPv9DxljHljIdxIREZEzQ2sRlu/PgN3W2mcYY94FvAf48FxuYIwJA93AK3ALP99vjPmWtfZRv421didwVaD9AeAbObd6mbX2yLy/iYiIiJxR6sEqg7ew8hustbd6p/YAl87jVtcCu6y1u621k8BXgNcVaf9y4Alr7d55/C4RERFZJtXVg9UV/xBwxSLfdSddQ58u0eZGYJMx5iHv82rgzmADY8zdQGOen/2wtdZvuxF4MnBtP3Bdkd/7ZtxC00EW+L4xxgKfs9Z+vsSzi4iIyBKrroC1fK4CPmqt7QEwxnwBeDjYwFr74jLuY/Kcy7vatjEmBrwWNzQZ9EJr7aAxZh3wA2NMv7X2rjJ+t4iIiCyR6gpYpXuazpRVuGFBjDER4Cbgk8EGZfZg7Qc2Ba5dCAwW+J2vAn5lrT0cPGmtHfT2TxljvoEbdlTAEhERqSDVFbCWz2PA9cBtwAeAXmvtnmCDMnuw7gcuM8ZswU1efzPw1gJt30LO8KA3FyxkrR32jm8CPjGXLyIiIiJnnia5l+d24LnGmF3As4APzucm1tpp4Bbge0AK+Kq19hEAY8y3jTEbvOM63JuGX8+5xXrgp8aYXwP34YLed+fzLCIiInLmGGvzTgFaFsaYH1trb1ju5xARERFZCPVgiYiIiCwyBSwRERGRRaaAJSIiIrLIFLBEREREFpkCloiIiMgiU8ASERERWWQKWCIiIiKLTAFLREREZJEpYJXJGBM2xtxqjHnEGJM0xlw8z/t8yRjzlDHmN0Xa3GyM2WmM2WWM6Sz3moiIiFQGVXIvkzHmI8CwtfZWY8y7gCustR+ex31eAowAt1lrn5nnehi39uErcItD3w+8xVr7aLFr8/1eIiIiUlxzZ28jsBlo9vabgc0D27f+bqGf0WLPZfAWVn6DtfZq79QeYOt87mWtvcsY01ykybXALmvtbu93fwV4HfBoiWsiIiIyT82dvSFgAzODVLO3rQ40zeA6OQaK3a+qAlbrjtYPAVcs8m13JtuTny7R5kZgkzHmIe/zauDOYANjzN1AY56f/bC19s485wvZCDwZ+LwfuK6MayIiIlJCc2dvA9leqObAtgmIBpqexIWon3r7vd7+wMD2rdOlfk9VBaxldBXwUWttD4Ax5gvAw8EG1toXL9LvMnnO2TKuiYiICKd7o85nZoDyA9WaQNNgb9RPyYaovQPbt55YyDNUVcAqo6fpTFmFGxbEGBMBbgI+GWywiD1Y+3Ep2nchMFjGNRERkXNKc2dvPTN7ovyeqYuAWKCp3xv1c2b3Rk2diWerqoC1jB4DrgduAz4A9Fpr9wQbLGIP1v3AZcaYLcAB4M3AW8u4JiIictbJ6Y3KnR91XqBpBvd34wBwLzOD1ImB7VuXdMRHAas8twPfMcbsAu4B3j3fGxljbgduAM4zxuwHPmat/aIx5tvAO621g8aYW4DvAWHgS9baRwCstdOFromIiFSz5s7eOmYHKL9HKtgbNYwLTvcSGNID9g9s3zq5ZA9cgso0iIiIyJLweqPWkz9IrQs0zeCmwAQD1IC3HV/q3qj5UA+WiIiILCqvN+oi8gepmkDTUVxoup+ZIaqieqPmQwFLRERE5qy5s9fg3shrBrbk7IO9UZZsb9QDZEPUXuBYNfRGzYcCloiIiBTU3NkbxtVhbCYbovwg1RBoOoYLTn6I8nuknqz23qj5UMASERERmjt7a8kO5QVDVG4BzqO40kXf9fYD3v7ps7U3aj4UsERERM4R3rDeKmb2Rvn78wNNcwtw+kFqYGD71uEletyqpoAlIiJylvHe1ruAmT1R/r4p0PQULjg9xMzeqKqfZL7cFLBERESqVHNnbw0zK5n7ISq3dtQxXHi6k2yIGgAOD2zfmlmapz23KGCJiIhUuObO3jize6K24Hqp/HVq/bf19gC/8PZ7cMN6J5f2iUUBq0zGmDDwGeBG3Nj066y1u+dxny8BrwaestY+s0CbAVyl2jQwba29JnDtZuBWXCX3L1hrt8/1GUREpPIUWKDYD1OrAk0ncb1PvwHuINsbtVfDepVDldzLZIz5CDBsrb3VGPMu4Apr7YfncZ+XACPAbSUC1jXW2iM558O4dRFfgZt8eD/wFmvto3N9DhERWR7Nnb0R3Jt5/vyoi8kO660INB1i5rwof39Iw3qVTz1YZTDG1ANvsNZe7Z3aA2ydz72stXcZY5rn+SjXArv8njNjzFeA1wEKWCIiFaa5szdGtifqYmYGqXCg6UHc3yu/JBui9gxs33piyR5WFl1VBaxUS+JDwBWLfNudif7Up0u0uRHYZIx5yPu8GjdR8DRjzN1AY56f/bC19s4854uxwPeNMRb4nLX28975jcCTgXb7gevmeG8REVlE3rIwzWR7o/yeqY1AyGvmlz3YA/zE2+/GDeuNL/EjyxKoqoC1jK4CPmqt7QEwxnwBeDjYwFr74kX8fS+01g4aY9YBPzDG9Ftr7yI7kXHGr17E3ysiIgU0d/Y2ke2NCg7vBetHTQP7gJ3Ad3A9UruBfZofdW6pqoBVRk/TmbIK918bGGMiwE3AJ4MNFrMHy1o76O2fMsZ8Azc0eBfuv342BZpeiHtjREREFkGgEGcwQPn7NYGmk7i/Fx4k+7bebuDAwPat00v5zFKZqipgLaPHgOuB24APAL3W2j3BBovVg+XN9wpZa4e945uAT3iX7wcuM8ZsAQ4Abwbeuhi/V0TkXOIFqbXMDlFbgHig6RguOP2cbIjaAxzURHMpRgGrPLcD3zHG7ALuAd493xsZY24HbgDOM8bsBz5mrf2iMebbwDtxb5B8wxgD7p/Pl6213wWw1k4bY24BvoebIPkla+0j8/9aIiJnt0DpA3+SebBnqi7Q9CQuPP2QmT1SWl9P5kVlGkREpOo1d/aGcdMmckNUM1ATaHqUbC9UsEfquIKULCb1YImISNUI1JC6GLgEF6QuAS5i5t9ph3Dh6ZfefjeqaC5LSAFLREQqTqBH6hKyw3uX4GpI+X93+UvD7AbuJtsbNTCwfevYUj+zSJACloiILBtvjlRukPKH9qKBpn6Q+hnwhHe8Z2D71lNL+bwi5VLAEhGRM84LUhvJ9kQFg1Qs0PQgLjzdw8yhPfVISVVRwBIRkUXjBakNzOyN8udKBYPUYVxP1H24EPUErkdKQUqWR1c8hqtnGQeacrbcc+5z19DrC91OAUtEROYsUP7A743y91uY+dbeU7jw9IC33wPsHti+dXRJH1jODV1xgyu/EWeuYWnmQtu5MsAIbgHuk8AJZi5dN4sCloiIFJSnjlSwnlRtoOnTuJ6ofyc7tLd7YPvWkSV9YDl7dMVX4ILPSrIhKO599o9ztyZmLqSdaxIXkoa9/SCQwoWm4BZscxIYpWtoToVlFbBERCRY2fxSZvZK5QapI7jw9E1cb5Q/tKfyB5JfVzxC4UDkh6KVeT7H8t3Oc4psT9IQ7t/DocCWPyx1DU0s7pcrTAGrTMaYMPAZ4EZcV+HrrLW753iPTbjlds737vF5a+2tedrdDNyKS+FfsNZuL+eaiEg5mjt747jwdCnZQHUJ0BBodgz3l9a3yL61t1tB6hzWFQ/h/h0ptzfJ3+qL3DXNzGB0AHg051xucDpB11DFL5ytSu5lMsZ8BBi21t5qjHkXcIW19sNzvMcFwAXW2l8ZYxpxBfBeb619NNAmjFv78BW4xZ3vB95irX202LVF+IoicpZp7uytxfVABYPUpcxctHgY2IULUae3ge1bh5b2aWVJubDk9xQFt3iBz/6cplCBO1pmDqkN4XqX/ON8n4eAMbqGKieILCL1YJXBW3T5Ddbaq71Te4Ctc72PtfYg7hVkvMWcU7jXloMB6Vpgl987Zoz5CvA6r02xayJyjmru7I3iKpnn9kptAIzXbBIXnu7x9n6o0lp71c5N7G5gdlgqtjWR/Xcj1yRwnGwo2snMyd1DeT4Pz3WO0tmuqgJW97a+DwFXLPJtd3b0tH26RJsbgU3GmIe8z6uBO4MNjDF349J9rg9ba+/MPWmMaQaeA/wi59JGZr6ZsB+4roxrInKW8yacX8DsIBWsbp4B9uL+w+s/yfZK7R/YvlV/AVY6F5ZqmVtYilN4Yvc0Liz5Qehx77jw1jWk4q2LoKoC1jK6CviotbYHwBjzBeDhYANr7YvLvZkxpgH3ps37rbW58xny/ReFLeOaiJwlvAnna8jOjfLD1MXMfJV8EBee7ibbI7V3YPvWip+fcs7oiodxAWg1sCqwFQtM0bz3cuE5GIYGKBWWYPxsHYKrdFUVsMroaTpTVuGGBTHGRICbgE8GG5Tbg2WMieLC1b9aa7+ep/1+3EKmvgtx/yda6pqIVKHmzt5Gsm/sBXul4oFmx3AB6hvMnCelopxLLRuYVpENTSuZHaD8rdBQnD9nyQ9CB3HlAoqFpRENwy2N1h2tYVw9rdrA3j8+fT7Znry90D2qKmAto8eA63FvAH4A6LXW7gk2KKcHyxhjgC8CKWvtZwo0ux+4zBizBfc2xZuBt5ZxTUQqWHNnb4zshPNgmFoXaDaGC1J9ZHukdg9s33psaZ/2HDIzMPmhaSUzA1S5gWkINxx3HPfP7jguHJ8IHGeH67qG0mfqa51LvDBUjws99WRDUKFtBTlBiZkBqpbiJSKCFLAW6HbgO8aYXbgJou+e531eCLwNSAbmc/0Pa+23jTHfBt5prR00xtwCfA83pv4la+0jANba6ULXRKQyeMN7F5ANUJd520Vk38CawvWK+9XN/UnnhzXhfIHcHKYm3BBrbkjK1+NULDCdJBuK/MCUu/kB6oR6l8rTuqM1igsy5YSi3Gv5fqbcMATuP2LGvc0/PolbumkscC5fu3zXC1KZBhGReWru7G1gdpC6hJl1fw7gwtPj3n4X8OTA9q3qvSiXKymwEheaVgf2wS14Lt+Ebz8wBUORH57ybeph8niBqD6wFQtD+QJT7vlyO3cywCguyIySDTjBrdT1YJtTyfbkkoVg9WCJiJTQ3Nkbxr2pFwxTl+J6qnx+PaleskFql+ZJFeCqe+cLSPnC00ry11+awoWko7gK8zu9z/45//icC0ytO1r9txHrc7YGskEn+Dm3TV3gc7k9RFPMDjZjuGWUCoWfYkFpMtmerJxeoDlSwBIRCWju7F1NtjfKD1JbyP4lk8aVQXgY98KK3zv11Dk/vNcVj5E/KOU711TgLqfIBqQDuD/nYzmbH55GzrY35Fp3tEaYHXhyw5EfkAqFI/96oTpXQZO4RYz9kDMKHAoc+5vfJhicZoSjZHtyaiHf/WyjgCUi56Tmzt4a3Nt7fojye6ZWBZodwYWn+8gO8Q2cU2UQZs5pOi+wrQns/eN8b1KD+8vZD0a7cS/s5PYyHQWO0zVUtT1+rTtaY7iA42+NZXzODUfl9BZZZgYif3sa92edez5fWBoFxhSKzhwFLBE5q3nFOc8nG6T8nqlNZIedJsjWk3rc23YNbN96YskfeKl0xaPMDEfFAlS+vytO4QJoMDT5n2eGpypYN87rOWpgbgEpdysnHPkBZwQ3rHwU1yM6I/jkfM4NR0s6l0jmRwFLRM4azZ299cwe3rsUN1zi24/rifo+2V6ps6PKuettaiR/cMoNUIWG6I6TDUp7yM5v8vf+VlEFLL3eo0Zva8o5LicwrZh911nGyQajEdxE+f1kA1NwG87zeUzB6NyhgCUiVSdQCuFyb7vM228MNBvGBag7yM6Tqt7inF3xWmBtGVu+XpRJssFoL26h+dzg5Hqcuoamz+j3KCAwKTtfQMo9zneupsSv8OcaBYPPU8wOQoXC0WiyPXnOTJKXhVPAEpGK5s2VuoSZYeoyXK8DuPkoT+KqYP8HLkg9RrVMOndDdefhwtF5uMKja/N8rs/z0+O4eTdPA0lv7wepYHBaksngrTtaQ7h/Lk1ke47842KByd8KracH7p+zH3aGcSUXBgLHw0WOR5LtyYofppSziwKWiFSEwPp7lzMzTG0mO1fKr3T+XVyI8udKFS34tyxc7aZVZAOSH6JyP6/K89PTZIPTbuDewOcjuJ4X9+r7GQhOeYJS3NtKHTeSv5xC8HudxAWlk2SH2E7mnM93PKLhNakmClgisuSaO3sjQDOzw1QwbBzEBagf4sLUY8BgRcyVcjWczgPWe9s6b1sf2J/H7LBhcRO/n8a9Cv8b79gPTEe8/dBiVAX3glI9cwtJ5QQlP/wMefvBwLG/93uPgj1Kp6q5rpHIXChgicgZ1dzZ28TMIHU5rq5U1GsySfYNPj9I7RrYvvXk0j8tfi0nPzzlhib/eA2zawyN45bbeApX1sEPTcHt6HyLXQZ6lVaWuZUTlEZxgcgPRQeZGZTyHQ8n25PLMk9LJJ9USyKEyzPRwD6a83ku5yM51wseJ/pTf1rouRSwRGRReOUQLmRmj9QVzFzM+Biu2va9ZMPU3iVbNsaFp2BvU74QtTrPT46SDU+7vGP/s78ve55TYEJ3uWHJ3wqFpUn89fDc1s/sYOQPyfnHJxWUpBAvtMTIBo/gcTBoLDTMLMb5Yv8RsVBp3ND2lLfPPS5IaxGKyJw1d/bGcBPPrwhsl+FCA7g1xAZwYcqfdP7YwPatx87YQ7k5T6txNa/87YLA8XpcSMk1zMygFAxP7nPX0GixX+31LjWRXUA4uKhw7uafj+a9mfuzOzGH7TgaeqtqqZZEmPxhJpZznBt0Cl2P4UJIsev5flfw+pkKLVPMDitTZ+j8VJ5rxc7nhqfpRH9q3v+7UsASkaK8BY0vB1rIhqktZN/4GsMFqZ1ke6V2L3q18654DS4k5QYn/3g9s0PLGG7Y6zBuzlNuiHq6UOXw1h2tdcwMS6uZHaBKrZUH2Z6jcoLSCVw5gOWfZ3aO8HpqanChwt/HipzLPV/Oz/j7QqFnMcPMNK5Hcyqwn8o5V+713PP5gslcwk96IYGl2ihgiQhw+i2+83ABKhiogrWljuKCVD/ZULXwieeuQOZKZoen4JY7dGfJThbP3Q56+9PDdq07WsPMDkfFglOhwpNjzFxEOLjPPTek2knl8YLOClwQ8fc1OeeCx+WEmxUlrtdQvDREOTK4lQAmA/tJXKX7yZxtLkGm3OvTweNEf0rhvEIoYImcgwLzpa7I2YIh5kmyvVL9wM6B7VuPzusXugC1CtiQZ/MDVW6hyFNkg9KhnONDwFOtWy5K496SC1YuD66PFzxfqJdpmvLC0jHgeLI9OTGvP4Mq5A1d5Qab3BCUG4bKuZ4vQBUaMi3FMjvclAo8+a4XOpd7fsb1RH9KAVryUsASOcs1d/ZGcYsaX0G2V+oyssvHpHFv8e0MbI8NbN9adN7RDNklWjbiAlNw7wep3AB1Avd6/0FygtRdtSuOfnDdeZGJUCh3yZd8x/kql0/jetuCxTb9tfGOMDM4DVfr/KVUS8Lg/lxrcUEluA9uuefKbTvfF6GmcIHGDyinAvtTec5N5Fwv9LO5e7/Xpir/+cnZTW8RipxFvLX4LmfmMN/FZP+3PkZ2+Rh/mG9PWfOluuL15O+B8rfcSuPDuAA1APwcGDwZMod76+vHvtzUaAdi0UZmViq/nNLr5J0gG5YeChwHw9QRKjA0pVoSEVyozd3qcWGmPs+1UiFoBbPLRZRyCldSYjzn+OkC5/MFnqIhSMNUIgpYIlXLqy/V4m0Jb78p0OQ4LkDdQ3beVOFFjd1beOtwQ4fBbaO35YaecVyAOjANDw5Eo8d/uaJm/Be1K6YfqonxdCTSRDZAXUO2anluIEiTrRM1ADxA/uB0PNmenCrzj2fBvEBUH9jyhaBi4Si3Tb6etnwsLgiPe3s/7IzgL7I8OwQVCka5xwo/57jubX0h3FB5CDf/LBw4zt2Hy2ibe24+bUNz2MK4/w+Zy34u95/T1tHTdmOhP2sFLJEq0NzZG8eFKD9IJXC9Rr5BXIC6g+ww35FZa/F1xWN0sZHZIepC737BeTDTwOAkHNgdi+57NBYbf6QmNt0fizEQjUROhkKNGLMWN9z4AvL//4lftfxp4FFmFt18ChcYTizmW3PesFktrihnAy7gBPeFzuVeK7V4sC+Nq5M1lrMdy/k8igs5+doGr5/SkNfi6N7WZ3B/yfp1m8I5WyRnn3u+VLvwItzbP7cUQWahE/qXUgb3HxvpPPvMIm3TBc7n+52F9gVpDpZIhWnu7F3FzCDVgpvP5DuAW9jY33YObN86dPpqV7yJ/AHqQlwv0ukepDSMHQ6Hn94Tiw49Ho2OPR6LTj0RizIYicSOh0L1GLMeN9cpd3L4KDOXd8lbtXw+PU7e22QNzFwMOHfx4EZmByL/uD7P8+bye4lGvO8ykrONBvb+li8UjSX6U2fVIsJeKMmtXJ1vK1b1upzP871fobBS6Hg5+X9Z+8Uq8x3n/kUfPE7n3CN4rtrapnN+Lt+9/M129LRVTjiZJwUskWXU3Nm7mmzPlL8FK5/vw/VM9ZMNUyfpitcBF+Vsm3FDhE3gVasMhcKDkcjovmhkaCAaHdsTjUzvj0bM4XA4eiIcrp+CVRiT+5fQKbJ1og4Fjk+fS7Yn89aO8nlvn+WGo0Zmh6Z85+opPq9oGje/KzcM5Qaj3NAUPD++XENlXoDxC0sG6yL5x/kqZ+crMFnqON+1Qr8vGGDOZFVsn/8Xc3CbmsPnNKWDS267dIGfD7Ytdu9i7fLdO9PR06bh2HOYApbIEmnu7F3LzF6pBK5HCVyPyj5ciOoHUltD9+zujv39Klxo2szMMHVeBjgeCoUPRSLRfdHIyEA0MjIQjU49GYmEnoqEY0OhUOOEMcYaE/wf+RTZQpt5wxM5E8RTLYka3Lp2wTXu4nmOg58bSvxxTFJ4QeDguXznFzyE5s1DiTK7cGQ5tZTyFZosdo+anPtFmfvE9FKKFY7M9zlYPym3unWhbS4BqOj9FDzkXKCAJbLIvIKda5kZpBK4oTZwYWoASEWY3vn68M+OfCDytYmN5uh6AkEqA+cfC4eih8KR6KFIOHYgEpnYG42c2h+NpA+FI6Fj4VDdKWPspDGTgRDlL9ibW/7gMHDojXenh974UxtndlHNQmFpJcXnIo2RrUTur3sXPM4XmoYLDat5vTt+qFnoVlvg/EJqLvnSlFdjKV/tpFLBp9Rx7uf02TCcInK2UcASWSBvztSVwDO8/ZVkC3ZmgIE4I4+/NPTro78Vvm/iRaFkuMGcuhBoHjNm82AkUjcYCUcPRiKxwUiEvdHIxMFIxBwJh6KnjElPGjM5acxkxpgMLrScDk8mYwfXn+DYdcqrpG8AACAASURBVDvt+EuTmakLj7KCwhXK/SCVr/fE4gJQsbA049rR1Vee/PWzOqLMLivgvzWXr+RAvmvB0gPzKTsAM0sHjOd8zi0r4IehckNR8NwkMNHR06bikiJSlAKWyBx46/IlmBmozgcwWHseJ/ZfE3rsqeeHHh1+XqifLWawbiRiLzkYjlw4GI1ED4XDsUORcPTJaGTqQCTCiVA4NGWY8ENUGk5gzIFw2g5uOMaJ5sN2/NJBO3XxIZve9DThukn80gfnkQ1NhV7/HyVbTPOotz+WDkWPj9WtOzVSf+HEcOOmyZNNzenhhk3YUCTfG3XByePBQFRH4aVk8pkkOzF8nJkTxXPLDBTaCl2fVA+OiFQaBSyRApo7e1fgCnb6vVLPAC4KkeECjkZbQvuGrgo9cfxZZtfkqujB+lPR0XWHouGVByPh2OFwOHYoEgkNRsJ2NBRKTxozOYWZmLJMNI7w9LoTHNtwjJGNR+34hmNMrz9uWT1MtGaaOC485asXBXDCwpGpaP3QRGzl6MSKVeOnalZPjdeeN31qxer0qZrVZqImHpqKNkRsKBIsVRB8266cYDTJ7MnhwVICfjAKfs4NTaePO3rapuf4xy8iUtUUsEQ4vZzMpWTD1JUGe8k6jscuDQ3WXGH2TVwefmK8oeZATTp6cuXRSCZ8KBKOHYqEY4ORSGg0FJqeskysGGe6YdRMrjnBsbXHOHHBMSbXDVm7epho4zi1NZM0hIJlEkJRMxVtCE3EmkYnY/HRiZqV4xM1KycmYk3TkzXxzGS00U5FG0LTkdpIOhyrwYT8N+2Kvek1gZvrlO/tuVLnRoDRjp62s6r0gIjIUlPAknOOt9BxM4E5U6s5mbjEDDZsDh2oWR8bCMej+62JnVhxPJKOPB3FHI6EoqPTYepHzHTDCOnGEdJrhhhZc5zh84aYXjlKuGGcmlAmmsmEG+10tCE8Fa2PTEUbMhM18bGJmpXjE7H49GSsyQWmaF0kHa5ZkTGRiDWhacyszip/TtQwbs7TyZxtKOea//lkR0/bObMYsYhIpVIldzmreW/0XQA8E3hGE6OtzzUHnn1RaDC+KravvjZ2OBaOHTVD0YnIUQsnJkKRqRHsiSHSDSPhzJqT4clLTjDeNFpro+nGUCbUmLHhpkw63MhUrDEyUbOyfmxNPH3igiYzFa2PZMIxkzHhaWvCU5lQeAxCGS88hXFB6HhgO5Hz2T93AhjRq+wiItVLAUvOKt4k9CtDZJ652Ry69mXhJ69eHdu3oSE2WBuOHI1l7Fh4ws0OCseOW1t/0qQbxxrtxePrMjXTcZsOxcfTkZWZdCTOdKSJyVg8PbS+IXo8FE1nTHjKmvBkJhSZsiY8ZU1oGGOOMHOR4aNkJ5b7+xO4eUiV010sIiJnlIYIpWo1d/ZGgEsbGHt2a2jXi1bW7L2+NnZ4c609XFszMVzDWNpER8OmaWQlDeMr7YrJOBHblElH4unpcDyTjjQxFY2nrYlOeD1Ok9ZEpqwJDVsTOpAJRfdjjL8UTDBAHQGOdPS0Fa1mLiIi5y4FLKkK3lDf+ZeZ/c+7LProK2Oxg9c1Th/d0jB6oj42Ho2Ex+KmYTxuaifiYFaSiazMTEdW2XS4wVoTnrQmMmFNZBwTHsPaw5lQeF86XLvHhsKDzK5mPqzeJhERWQgFLKlIzZ299YnQrudtrnnov6yZHLyhcWzqovrxWH10Ih6OTawmll5t0uHVTEdWWkzYQiidCYUnLNGRcMYeNNY+mQlFdk1FG/oz4dheXHHOQeC4wpOIiJxpCliy7Jo7e0NXhx9qbTa/ftuakYmbG06t2FwzUV8XmVppol6QSodrMRhrjUlbY06ZzMRT4bTtD9vobyajjcmJFasfBw4Ah1RzSURElpsCliy5t3+kY+2WsUMfjI9HXlszVb85MhWvDaVXhzBryYS8Ze+MsTAxaRh9KpS2D4fT9T+biq25/9SKNf3AQS1VIiIilUxvEcoZ9ca/fEfoOYdO/VH9ROydsan6lsj06oYXZS4P2dA1p9tYM5XJRI6PwIF94XTmzujU+m+M1m/55R9/8TXDy/joIiIi86aAJYvqs3/SflVkMvzx6FTddeHpVatvts+O2lDT6evGHs/Y8JGRdGjvbgNfD09d+s8T9S37NC9KRETOJpU2RHgCeGi5n0PKt7KhfkO8tm5DNBSrDVETgmjg6jSWycx0ZmpyfOrUsSMnh/dOZ9JagkVERM4KxaY1qQdLymaMCZ0Xb2xuiNWujZhYjSFmsv8KWbATpBmfOpWeOHlsZGTf+MTkyWV9YBERkWVSaQHrIU1yrxxv/0jH2svGDv1FfKTp9bGJteeHMpuiNtQAgLGjZEKHJtLh4wemI+NfC01FP/aeL33u1DI/soiISEWotIAly+iVf/7pjc9N/+yW809E31x7at3GF09vjmbC/mT0CZsJD45PRVKPT8dG/+7BDeu+9E8f+evKGV8WERGpIJU2B0tlGpZIc2dvaDUnL752xQ9/78LR/W9adWLlJbXjF9VgtpAJrQAywOHxdPjoTmOnb105tOm2N339w1p8WEREpAwKWOeI5s7eWITpxAuj921dF02+fu2xY1fGhzbWrjjVEpqObcYSwpqJacyhvZHM5OfJXNn99n/5vdHlfm4REZFqpIB1lmru7I0Bz9xkDr7oWSseeG1D+InnnHf0WFPTkYvD9RNXmIkVW7ChSAbGToYzQ3eG0mtvHWnc8vOOnjb1UomIiCyQ5mCdJZo7e2uAZ8WYuvracPKm32545Nq60O5VdU+PrVi790JTN/kcc6rucpsJx6ZP1WUORqYzXxuvWf+lTDj2qGpQiYiILC4FrCrV3NlbCzwLuPp8jj7/NeHktevqH147GRuszxxJh5sfqWfV6LPMWMOzM5OxVVOjNbHDken0NyZja3fYUDj5TvVUiYiInDEKWFWiubO3Dng2cLUh89xLzMHnXht5cOWq+t+sOr7iSGxsiPDKx2DT4WY7UXft1Ej9pelja+qPGRP62WS04TZMqO9d/9Q2sdzfQ0REpBoYYzYAfw9cDNQCf2Ctva/cn1fAqlDeHKpW4FrgeWHSz3iuebzx2sivG5sakjW768Zq96VDkdonSF+3MzwV5Vmnjq96vj2yfsN4OhzbNx2p/b+Y0H909LTtX+avIiIiUlWMMRHgO8CfW2vvMMbUAeG53EMBq0I0d/aGgBbgebhQdVWUqRXXhfprbwr/fLKhvn/Nw3Xpul/URGJrD4amrv5laOTVe1dMD8efH3pq3QvsRM3KselIza8w4X8BftDR06YlaURERObn9UDKWnsHgLV2bK43UMBaJs2dvQZoxgWq5wHXAI0xpswNoYeG3hC+a3RF3RNr7qsz8f+srakdS4fHn90fPvHuXzKxcrTOHFn3Eh5LvCA6EWsas6HoD4F/AR7UhHUREZEFuwq4dyE3UMBaQs2dveuA68iGqrUAEaYP/Xb47ifeEu6rXRN58tIfNdRc/NW62ob90brJ2ChHbv6RefL6lF1bM7UivHfzKyceuvwF9VPR+imM+SHwhY6etseW8WuJiIicbQ7h5j0DYIxZa619ei43UB2sM8grnfBc4PnA9biJcgDHDZn7fzf848PvjXx93Wpz7Pp76mov/F59Xf0vV9QwHAoNX3CIB//wBxl76SCXYkJm36abDu676OXrpyN1EeD7wBc7etqeWKavJiIictYyxjQAXwYuAaaAj1prvzWneyhgLR5v2O9iXJh6Pi5cxYBJ4CHgnneF7xjojNyeCBt786Fw+OJvNdQ33dFQb56OhCdGjdnTOmB//iffyqyNj/F8C6cObHzpb57Y8trN6ciKdcAvgc929LSllus7ioiISGkaIlyg5s7eOG7Y73pvW+dd2gN8Dbj3A5F/e+R9kW+8CHhDBt73YE1N7ZebGiZ/UbuibjgUOpYx5icvfyjz83d9J311CF4LjJ9ouvibDz37ls2ZcM3zgAHgA8BPNcdKRESk8ilgzZHXS3U58GLgRcAzAAOcBO7DTYq7d2D71kN0xS8H3gB8agIa/7OxfvK2pqb0k9HI1LQxx4F/+q37Mj/8wx+mXwf8KXBqOrziX+657mOZqVjTW4Fp4K+Br3X0tKWX/MuKiIjIvGiIsAxe1fTn4QLVi3C9VBZ4BPg5cA/wyMD2rRm64nXAzbhglRg1ZvoLK5sO/ltjw7qhUKgGYx4DvvL7fekfv/YX9k1AO24Y8asPPvu9Pzm+6or/hhvzvRP4dEdP25wm1YmIiMjyU8AqoLmz9wKyvVTX4ELQGC5M/RT42cD2rcdO/0BXfCPwJuB1QMNT4fC+v1yz6tiP6movT7sCZQ8A//zGu9O/eONP7StwQ35rgb6JaOP//tkLt98IvAs4DvzPjp62ny7VdxUREZHFpSFCj1foMwG8DBesLvEuPYmbS/VT4MGB7VunTv9QV9zgwtebgZcAmcFw+GfvW792qj8WfRHGbAJ+BOxIticfSbUkLsSV3b8eSAF/1ndD91PAX+KGGr8H/FVHT9vJM/6FRURE5Iw5p3uwmjt7w7g3/V4G3IAb+ksDDwJ3A3cPbN+6b9YPdsXDwI3AHwKXAUMHw+E73rZhfeRwJPIaYAXwXeBLyfbkQKolEQXeBrwTN6+qG/i3vhu6Xwh8HAgBn+roafv+mfu2IiIislTOuR4sb42/63Gh6iVAHJjADf39A/DTge1b8/cgdcVjwKtx86Y2AntOhEKfumnThlXjodDvA43AD4DPJ9uTewBSLYnLgU8AlwJ9wN/03dB9FHgPLqDtBP5UawaKiIicPc6JgOWFqhcCN3n7OmAEuAv4MXDPwPat4wVv4Cau/zbw+8Aa4NFp+NvnNW+KThvzJ8AGXI/XPybbk48BpFoSIeAPgG3AEPDBRH/qru5tffXA3wIvAL4O/I3WDRQRETm7nLUByxv+uxZ4Ja63qh43gfy7uJ6kX86YT5WP67H6HeDtwCrgfuAvntO8aWTamA/i1ip6HPivyfbk/f6PpVoS5wP/07v+Q+AvE/2pE93b+s4HPosrRvqpjp62ry/eNxYREZFKcVYFLG+i+rNxPVU34kLRCC7kfA94YGD71tL1pLriEeA1uLf61uHqW/W0brloN3ALrjfrOC5EfSvZnsz4P5pqSbzAOx8BPgp8J9Gfst3b+q7E9VzVALd09LTdtxjfWURERCrPWRGwmjt7N+DmRr0GuAA3p+ouXKj6+cD2reUNwbm3Al8MvB+4CEgCH2vdctEDuF6w/wWsBm4HPpdsT476P+oNCf4x8A5gF/Cnif7UPoDubX3PxfVcnQD+a0dP2+4FfmURERGpYFUbsLyFlF+GW1rmWlzhz18A/xv4ycD2rWNzumFX/BLgQ969BoAPAne3brloLfBp3IT4x4APJtuTjwZ/NNWSqMeVWngB8C3grxL9qQmA7m191wGfAQaB96hwqIiIyNmv6gJWc2fvJlxBz1fj3tobBHqAOwa2bz005xt2xZtwE9F/BxgF/gb4Gl1D0607Wm8GOnF/Tp8Fbk+2J2cMMaZaEhu8a5uBTyX6U6fnVXVv63sJ8Fe4wNbR0dN2DBERETnrVUXA8uZWXYcr6PkCXK2qHwLfxE1WzxT58fzccODLgf8OrMQVE/0cXUNDrTtam9jR2omby/Uw8NFke3JWGYVUS6IV17sVA/4k0Z86Pa+qe1vf83FDio/h5lypeKiIiMg5oqIDlhesXomb17QFOAb8H+DrA9u3Hpn3jbvi63A9Uy8B+oFb6Bp6DKB1R2srsB04DzfcuCO31wog1ZJ4Pq636wjw7kR/asC/1r2t75m4cLUbNyw4Mu9nFRERkapTkQHLC1Y3A3+EG3p7AvdG3p1lT1jPx/VavQF4H+673wp8ma6hdOuOVoMbevwA8BTw9ty5Vr5US6IN+BQuQN2S6E+dHvrr3ta3xbvvMeBPFK5ERETOPRUXsJo7e1uA/wFciasx9d+BH89rGDCoK74SF9Jegqtn9Um6hvYDtO5orQM+ghsSvAvoSrYn8w7ppVoSvwV0Ab8B3pfoTw3717q39a3DVYNP4+ZcHV3QM4uIiEhVqqiAFW4873zgNlw5g78AvrfgYAXQFb8Wt1xNHPdG31foGsoAtO5oXYerT3UZLhzdFqxrFeSFq48DDwAfSvSnTr+p2L2tb4V3n0bgnVr6RkRE5NxVUQHLTp0aw002/8eB7VuHS7UvqSsews3f+mNgL/Bef64VQOuO1stxbwA2AB9Itid/VuhWqZbES3E9Vw8A7/fLMAB0b+szuN6xy4EPdPS0PZb3JiIiInJOqKiAlTk1cnJg+9b/tSg3c+sHfhxXK+vbwKfoGjrlX27d0foC3GT2EeCd/hqC+aRaEtfg6lylcD1XEzlN/gA3vPj3HT1tP12U5xcREZGqVVEBa9F0xS/EDQU2e/vb6Rqy/uXWHa3+JPUngPcn25MFi3+mWhJXevfYD7w3OCwI0L2t7zlAB/B93PCmiIiInOPOvoDVFb8CN5cqjCu/MGPNP6946Cdwk9Tfm2xPFnzLz1u02V/ipiPRnxoKXu/e1rcS+CRwAPhkR0+bnX0XEREROdecXQGrK341bqL5SaCDrqG9wcutO1pfi5s8/yCu56rgcjqplkQdrueqBtiW6E/N6OUKzLtaBby9o6dtdPZdRERE5FwUWu4HWDRd8ZcCfw8cAt6RJ1y9CheIfoHruSoWrgxu/talwJ8l+lP5Fmd+Da7kw9919LT1L86XEBERkbPB2RGwXLj6a9yyNO+ia+ip4OXWHa0vxL0B+CvcYs2nZt1jprfhJsd/NtGf+nnuxe5tfetxC0P/Cvh/C35+EREROatU/xBhV/x63NuAjwLvoWtoRs9U645Wf9max3GlGIpWgk+1JJ4D3AL0AbfnXveGBv8CF04/3tHTtvA6XSIiInJWqe4erK74M3GLLe/B1bjKDVfn4+ZRHcENCxadJ5VqSazEvV14APhEoj+Vb9L6q4DrcSUZDiz8S4iIiMjZpnp7sLriwfB0C11DM5a28Za/OT1JPdmePDb7JlnevKs/w01ab0/0p2a9Xdi9ra8ReD/wCK4gqoiIiMgs1RmwXBHR0+GJrqEZ4clbuPkjuEnq70u2J/NNUs/1SuDlwD8k+lM7C7T5Y1wAe5+GBkVERKSQ6gtYXXGDWwz6UuB9dA3lC0//BVdZ/R+S7cl7St0y1ZJYB/wp8DAFioV2b+u7HHgT8LWOnrbUPJ9eREREzgHVOAfrt4Cbgc/RNTQrPLXuaL0C94bfzym/svp/A2LAxxL9qUI9U+8FhoF/nPMTi4iIyDmlugJWV3wTrqfpV8A/5V5u3dEaw1VWPwF8LNmeLDmMl2pJvARXkuH/JPpTT+Zr072t71rcxPYvdvS0nczXRkRERMRXPUOEXfEQrpbVNPAXdA3lC09/jFt/8JZke/J4qVumWhK1uN6r3cC/5mvTva0vhOu9OogmtouIiEgZqqkH61XAs4HP0jV0OPdi647WK3EFQv8j2Z68t8x7vh24APjLRH9qqkCbVwAtwD929LQVraElIiIiAtUSsLriDcD7cAs035F7uXVHawT4GK5kw9+Wc8tUS2I98PvAtxP9qQfztfF6r94N7AK+O69nFxERkXNOdQSsbHmEvyowNPg7wCXAXyXbk7PqVxXwHm9fbNJ6G7AZN/dKZRlERESkLJUfsLriG3DlEb5B19Cs8gitO1rjuF6m+4C7yrllqiVxOe5txK8k+lMH87XxlsR5B7AX+OH8Hl5ERETORZUfsFzISQNfKHD93UAD8JlkezLf0jb5vBc4SZ43EQNeCFwO/JN6r0RERGQuKjtgdcXXAa8BvknX0FO5l1t3tF4EvBH4erI9uaucW6ZaEq24kgs7Ev2p4SJN3wEMorlXIiIiMkeVHbDc3CoD/N8C198BTAGfn8M9/wgYokjJhe5tfVcCzwK+3NHTNj2He4uIiIhUcMDqisdwS97cRdfQYO7l1h2tG3DzqP691ELOvlRL4grgRcC/JvpTY0WavgkYA/5zzs8tIiIi57zKDVjwUmAl8G8Frv8hbm7Wv8zhnn8EjABfLdSge1vfKtw6hnd09LSNzuHeIiIiIkBlB6xX4epa3Z97oXVHqz836z+S7cmny7lZqiVxEa7swv9L9KeKlXJ4PW5dwkLBTkRERKSoygxYXfEm4AXA9wrUvXozEKb8xZzBTYafpnjvVRg37+u+jp62PXO4t4iIiMhplRmw4CW4dRK/l3uhdUdrFHg1cFeyPTlrblY+qZZEHfBa4AeJ/tTRIk2vB9ajNQdFRERkASo1YF0PHAP681x7CbAa+MYc7rcVqKdI75Xn1bg3DO+ew71FREREZqi8gNUVDwHXAfcWGB58A3AYuKec26VaEgb3VmAKt5ZhXt3b+ppwE+u/29HTVmjhZxEREZGSKi9gwWW4dQfvzb3glWa4Dje5vdzq6s8FtgBfTfSnilV6vxE3ub13bo8rIiIiMlMlBqyrvf19ea69ztt/aw73ezWuptX3y2i3G9fTJSIiIjJvlRiwrgCO0DV0JHiydUerwZVu+EWyPXmonBulWhIrgJcDdyb6UxOF2nVv67sIV7n9jo6etnLXMxQRERHJqxID1mXAY3nOXwFsAH4wh3vdANRRetjvFd5e6w6KiIjIglVUwDJu3cGLgcfzXH45kAF+PIdbbgUOAQ+WaNcGPNzR0zZrQWkRERGRuaqogFUfZQWu/lW+Hqw24IFke3KonHulWhLn4SbEfzvRnyo4Ib57W99GXO9Y39yfWERERGS2igpYdTFT6x3OCFitO1o3ApuBu+Zwu5tw36/U8ODLvP2P5nBvERERkYIqK2BFqQMmgX05l17o7X8+h9u1AY8l+lN7y2nX0dN2YA73FhERESmoogLW4LA9DNxC11A659ILgP3J9mRu8Mor1ZJYhXsr8CfF2nVv61vrtfvhPB5XREREJK+KCliTaaboGvpV8Jy39uA1zK336sW471Y0YOHeMgQND4qIiMgiqqiAVUALsAJ4YA4/81Lccjo7S7R7EfBkR0/b7nk+m4iIiMgs1RCwnuPtHyqnsVdc9HrgJ8WWxune1hfDVY0va01DERERkXJVQ8B6LjCQbE8eK7P9tUANpYcHr8L1jClgiYiIyKKq6IDVuqM1hAtCpQqFBt0AjAK/KtHuemAa+OW8Hk5ERESkgIoOWEAz0AA8XE7jVEvC4ILTvYn+1FSJ5tcDv+7oaRtb0BOKiIiI5Kj0gNXi7R8ts/0mYB1wX7FG3dv61gCXo+FBEREROQOqIWBNAANltr/G25d64/A6b3/vPJ5JREREpKhqCFiPJduTuYVHC7kGeJrZleBzPR84Tv41D0VEREQWpGIDljfBvQVIldPem391DfBAifIMfrv7O3raCi4CLSIiIjJfFRuwgAuBOqC/zPZbgNXA/SXaXQCspfRbhiIiIiLzUskBq9nbl1tlvdz5V1d5+1/P9YFEREREylENAWtvme2vAQYT/anBEu2uAkaAJ+b5XCIiIiJFVXrAOpZsTw6XaphqSYTw5l+Vcd+rgIc1/0pERETOlEoOWJspvzzDRUATJdYr7N7W1wRcXKqdiIiIyEJUcsBqpvzhwVZvnyzR7lneXgFLREREzpiKDFitO1rjQJzye7BacfOqSgWyq3DrD5ZbGV5ERERkzioyYOGGB6H8HqxnAr9J9KdKzau6Cujv6Gk7Ne8nExERESmhUgPWRm9/oFTDVEuiDriUEsOD3dv6osCVqDyDiIiInGGVGrDWevvDZbS9HPc9Sg37XQLEKD1PS0RERGRBKjVgrQdGku3JsTLatnj7UhXfy20nIiIisiCVGrDWAU+V2bYFOAYcKdHuCmAUKFWIVERERGRBKjVgrae84UFwAau/2ALPgXY7VWBUREREzrRKDVjrgKdLNUq1JGpwhUOLDvt1b+sL4+ZqaXhQREREzriKC1itO1ojwBrK68G6BPcdSgWnzUBNGe1EREREFqziAhYuXBnKm4N1qbd/vEQ7TXAXERGRJVOJAWu9ty+nB+tSYILS9bJavHblFi4VERERmbdKDFjrvH05PViXALvLqODeAjzW0dOWXtCTiYiIiJShEgPWed6+VNkFcAHriWINurf1hXAlGnYu8LlEREREylKJAavR2w8Xa5RqScRxYaxowAI2APVo/pWIiIgskUoNWGPJ9mSp4bxLvH2pgHW5t39sQU8lIiIiUqZKDVhFe688W7z97jLb7Zn3E4mIiIjMQTUHrM24NwNLTYZvBg519LSNL/C5RERERMpSqQHrZBntLgL2lvEG4RbUeyUiIiJLqFID1kgZ7TYD+4o18N4gbAYGFvxUIiIiImWq1IBV6g3CKLCR0oVD1wMrUA+WiIiILKFKDVilhgg34p69VMBq9vYDC3skERERkfJVYsCqp/Qk94u8famApTcIRUREZMlVVsAKEcYt9FxqDtZmb/9kiXbNuN6wEwt7MBEREZHyVVTAMmET9g5L9WBtBE4m+lOlhhK3AHs6etrsgh9OREREpEyVFbBCpwNWqeC0ARgs45bNaHhQREREllhFBSzCRLyjcnqwigas7m19cWAVmuAuIiIiS6yiAlagB6vgHKxUSyIEXAAcKHE7TXAXERGRZVGpAatYD9YaIEbpIUIFLBEREVkWFRWwAkOExeZgbfD2pXqwNgOTwKGFPpaIiIjIXFRUwPJ6sDJAsYWZN3r7Uj1YG4EDHT1tpdYqFBEREVlUFRWwMpOZcaA32Z4sFor8HqxSAetCYP+iPJiIiIjIHFRUwEoPp48n25MfL9FsA3A00Z+aLNSge1ufwevBWsznExERESlHRQWsMq0DDpdoswqoQz1YIiIisgzO1oDlz9NSD5aIiIgsuWoNWE+XaHOht1cPloiIiCy5qgpYqZZEHdBA+T1Y5SynIyIiIrKoqipg4XqvAJ4q0W4j8HRHT9vEGX4eERER5gLGkAAAB/VJREFUkVmqNWCV6sHaABw8w88iIiIikle1Baz13r5UD9YFKGCJiIjIMqm2gFVyiLB7W18IF8QUsERERGRZVGPAOlGsyChwHhBGAUtERESWSTUGrFLDg+d7ewUsERERWRbVFrDWAEdKtLnA2x86w88iIiIiklc1BqxjJdooYImIiMiyqpqAlWpJGGA15QWskx09bWNn/qlEREREZquagIWr4B4FjpZodz6afyUiIiLLqJoC1hpvX6oHaz2lC5GKiIiInDHVFLBWe/tSPVjlvGkoIiIicsZUU8Aq2YPVva2vBmhCAUtERESWUTUFrHJ6sNZ6+6fP8LOIiIiIFFRNAWsNkAGGirQpuZSOiIiIyJlWTQFrNXA80Z/KFGnjByz1YImIiMiyqaaAtYbyJriD3iIUERGRZVRtAatUiYa1wFhHT9voEjyPiIiISF7VFLBWAsdLtFmHhgdFRERkmVVTwIpTfII7uB4sTXAXERGRZVUVASvVkogA9cDJEk3Xoh4sERERWWZVEbBwxUMBThRq0L2tz+DmaR1ZkicSERERKaBaAlbc2xcbImwAYpR+01BERETkjKqWgOX3YBUbIvSX0lHAEhERkWVVLQGrnB4sBSwRERGpCNUWsArOwUIBS0RERCpEtQUsDRGKiIhIxaumgJUGilVoXwNMA8NL8kQiIiIiBVRTwBpK9KdskTZrgGMdPW3FFoMWEREROeOqKWCVKjKqGlgiIiJSEaolYDVRepmcchaDFhERETnjqiVglbMO4Ro0wV1EREQqQLUErJUUCVjd2/pCwGoUsERERKQCVEvAaqT4HKxG3Hc5vjSPIyIiIlJYxQesVEsiDNQCI0WarfT2xQqRioiIiCyJig9YQL23L1YDa5W3Vw+WiIiILLtqCFgN3l49WCIiIlIVqiFg+T1YxQKWerBERESkYlRDwCqnB8sPWOrBEhERkWVXDQGrnDlYK4Hxjp62iSV4HhEREZGiqiFglTsHS8ODIiIiUhGqIWCV24Ol4UERERGpCNUQsMqdg6WAJSIiIhWhWgJWBig2v0o9WCIiIlIxqiFg1QMjif6ULdJmFZqDJSIiIhWiGgJWA0WGB7u39a0AVqCAJSIiIhWi6gMWEPf2Q0vwLCIiIiIlVUvAKvYGoQKWiIiIVJRqCFj1FO/BavL2ClgiIiJSEaohYJXqwfID1vASPIuIiIhISdUQsNSDJSIiIlWlogNWqiVhKL8H6+SZfyIRERGR0io6YAFRIELptwingFNL8kQiIiIiJVR6wPLXIRwr0qYRGO7oaStWiFRERERkyVR6wKr19uNF2sTR/CsRERGpIJUesFZ4+2LDf01o/pWIiIhUkGoJWMV6sBSwREREpKJUesDyhwjVgyUiIiJVo9IDVjk9WJqDJSIiIhWl0gNW0Unu3dv6IkAdquIuIiIiFaRaAlahIcJGb68eLBEREakYlR6wSg0Rqoq7iIiIVJxKD1ilerAUsERERKTiVHrAKlUHyx8iLLaUjoiIiMiSqvSAVQtMJvpTmQLXG7y9ApaIiIhUjEoPWCsoXqJBAUtEREQqTqUHrFoUsERERKTKVEPAKlbFvQHIUDyEiYiIiCypSg9Y5QwRjnb0tNkleh4RERGRkio9YJUzRKjhQREREakolR6wVlB6iFABS0RERCpKtQesehSwREREpMJUesDSEKGIiIhUnWoIWBoiFBERkapS6QGrnLcIFbBERESkolRswEq1JEIUCVjd2/oMClgiIiL/v717yYkjhsIwelGHbjpsJGyD1URiPUishm2QjZCEZ2eAS5ngpiKVit/RORMPPPHw03U9CBQbWFW1a2vvinBXVZsSWABAmOTAOmtr74rQb3IAgEjJgbVva2+CJbAAgEjJgTVNsAQWADCU5MD66BmsKbDuVzgLAMBsyYF12tanzv55W02wAIAoyYE1TbAeOvtTYJlgAQBRkgPrownW17YKLAAgSnJgbdv62NmfAuvYl94BAFY3emA9X91c9iZcAACfYvTA+rnSWQAAZhshsI49gyWwAIA4IwRW7y1CgQUARBohsEywAIChJAfW9JkGz2ABAENJDqxdVb1W1UtnX2ABAJGSA+u0qp4uftwdOvsCCwCIlBxY2+pfD1YJLAAglMACAFjYkIF1/f120/YFFgAQJzmwTuvjHz0LLAAgTnJgbev4R0arBBYAECg9sEywAIDhpAfWsY+MVlX9WuksAACzjR5Y9yudBQBgtlEDa99WEywAIM6ogXXeVhMsACBOcmAd+0yDCRYAECs5sHbV/0zDNMHyFiEAECc5sEywAIAhJQfWsWewzqrq8erm8nXF8wAAzDJyYP1e8SwAALNFBtbdt4uTersiFFgAwHAiA6ve4qqq/wyWwAIAYqUG1ratvbcIBRYAECs9sEywAIDhpAbWdEXoGSwAYDipgbVrq8ACAIaTGlhzJlg+MgoAREoNrOkZLBMsAGA4AgsAYGGjBta+BBYAECo1sLofGr3+frtp+wILAIg0XGDV3zcMBRYAECk1sDZtfX5n76yt3iIEACKlBtaXtr68szcFlgkWABBJYAEALCw1sOZcEQosACDSiIG1b6vAAgAipQaWK0IAYFipgeUtQgBgWOmBZYIFAAwnNbBcEQIAw0oPLG8RAgDDSQ2sOVeEDyudBQDgn5wcDofPPgMAwH8ldYIFADAsgQUAsDCBBQCwMIEFALAwgQUAsDCBBQCwsD9ZA/0uaAqnNgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"1|IMPORT PACKAGES\"\n",
    "import numpy as np              # Package for scientific computing with Python\n",
    "import matplotlib.pyplot as plt # Matplotlib is a 2D plotting library\n",
    "\n",
    "\n",
    "\"2|DEFINE PARAMETERS AND ARRAYS\"\n",
    "# PARAMETERS\n",
    "size    = 5        # Model domain\n",
    "steps   = size*100    # Number of \"dots\" in the domain\n",
    "\n",
    "# ARRAY\n",
    "c = np.linspace(0, size, steps) # Create array of consumption\n",
    "\n",
    "\n",
    "\"3|DEFINE UTILITY FUNCTION\"\n",
    "# Production function\n",
    "def u(x, CRRA):\n",
    "    if CRRA == 1:\n",
    "        u = np.log(x)\n",
    "    else:\n",
    "        u = ((x)**(1-CRRA)-1)/(1-CRRA)\n",
    "    return u\n",
    "\n",
    "\n",
    "\"4|CALCULATE UTILITY FUNCTIONS FOR DIFFERENT VALUES OF THETA\"\n",
    "theta = [0.25, 0.75, 1, 1.5, 2]\n",
    "\n",
    "# Ignore c[0] = 0 to avoid dividing by zero\n",
    "u1 = u(c[1:], theta[0])\n",
    "u2 = u(c[1:], theta[1])\n",
    "u3 = u(c[1:], theta[2])\n",
    "u4 = u(c[1:], theta[3])\n",
    "u5 = u(c[1:], theta[4])\n",
    "\n",
    "\n",
    "\"5|PLOT UTILITY FUNCTIONS\"\n",
    "### AXIS RANGE\n",
    "v = [0, size, -10, 5]    \n",
    "### BUILD PLOT AND POPULIATE WITH LOCI LINES\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "ax.plot(c[1:], u1, alpha = 0.8, label = r'$\\theta=%.2f$' %theta[0])\n",
    "ax.plot(c[1:], u2, alpha = 0.8, label = r'$\\theta=%.2f$' %theta[1])\n",
    "ax.plot(c[1:], u3, alpha = 0.8, label = r'$\\theta=%.2f$' %theta[2])\n",
    "ax.plot(c[1:], u4, alpha = 0.8, label = r'$\\theta=%.2f$' %theta[3])\n",
    "ax.plot(c[1:], u5, alpha = 0.8, label = r'$\\theta=%.2f$' %theta[4])\n",
    "### AXIS\n",
    "ax.axhline(0, color = 'k') # Add horizontal axis\n",
    "ax.axvline(0, color = 'k') # Add vertical axis\n",
    "ax.yaxis.set_major_locator(plt.NullLocator())                # Hide ticks\n",
    "ax.xaxis.set_major_locator(plt.NullLocator())                # Hide ticks\n",
    "### TEXT \n",
    "plt.text(size-0.10, -0.5, r'$c$', color = 'k') # x-axis label\n",
    "plt.text(    -0.25,  4.5, r'$u(c)$')           # y-axis label\n",
    "# SETTINGS\n",
    "plt.box(False)                                               # Hide axis\n",
    "plt.legend(loc=0, frameon=False)\n",
    "plt.axis(v)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. THE REPRESENTATIVE FIRM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The representative firm uses a labor-technology augmenting production function $(Y = F(K; AL)$ with constant returns to scale (CRS). The production function depicts the typical properties of diminishing marginal returns and the Inada conditions. Technology $(A)$ grows at rage $g$: $A(t) = A(0)\\cdot e^{gt}$.\n",
    "\n",
    "In terms of effective labor $(AL)$, and  normalizing the price of the good to one $(p = 1)$, the profit of the firm becomes (where the $\"\\hat{}\"$ symbol denotes the variable is in per effective labor units):\n",
    "\n",
    "\\begin{equation}\n",
    "    \\hat{\\pi} = f(\\hat{k}) - R \\hat{k} - \\hat{w}\n",
    "\\end{equation}\n",
    "\n",
    "The rate of return on capital equals the gross rate of return $(R)$ minus the depreciation rate $(\\delta)$: $r = R - \\delta$. In this model there is no intertermporal optimization. The firm rents capital from the household one period at the time.\n",
    "\n",
    "The FOC for the firm with respect to capital per effective unit of labor utilization is: \n",
    "\n",
    "\\begin{equation}\n",
    "    \\frac{\\partial \\pi}{\\partial \\hat{k}} = f'(\\hat{k}) - R = 0\n",
    "\\end{equation}\n",
    "\n",
    "In equilibrium the firm's profit is zero, therefore:\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{w} &= f(\\hat{w}) - R\\hat{k}       \\\\\n",
    "    \\hat{w} &= f(\\hat{w}) - (r + \\delta)k\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. EQUILIBRIUM\n",
    "\n",
    "In equilibrium, the assets owned by the household coincide with the capital per unit of labor. Recall that in this model households own the capital that is rented to the firms. Also, note that assets $(a)$ are expressed in per labor terms, while $(k)$ is expressed in per effective labor. Therefore: $a = Ak$. Then:\n",
    "\n",
    "\\begin{align}\n",
    "    \\dot{a}                      &= ra + w - c - na                          \\\\[5pt]\n",
    "    \\dot{(Ak)}                   &= r \\cdot Ak + Aw - Ac - n \\cdot Ak        \\\\[5pt]\n",
    "    \\dot{A}k + A \\dot{k}         &= r \\cdot Ak + Aw - Ac - n \\cdot Ak        \\\\[5pt]\n",
    "    \\frac{\\dot{A}}{A}k + \\dot{k} &= rk + w - c - nk                          \\\\[5pt]\n",
    "    gk + \\dot{k}                 &= rk + w - c - n    k                      \\\\[5pt]\n",
    "    \\dot{k}                      &= rk + f(k) - (r + \\delta)k - c - (n + g)k \\\\[5pt]\n",
    "    \\dot{k}                      &= f(k) - c - (n + g + \\delta)k\n",
    "\\end{align}\n",
    "\n",
    "Knowing that $R = r + \\delta$ and that $f'(k) = R + \\delta$, then the two motion equations that model Ramsey model are:\n",
    "\n",
    "\\begin{cases}\n",
    "    \\frac{\\dot{c}}{c} &= \\frac{1}{\\theta} \\cdot (f'(k) - \\delta - \\rho) \\\\[10pt]\n",
    "    \\dot{k}           &= f(k) - c - (n + g + \\delta)\n",
    "\\end{cases}\n",
    "\n",
    "In equilibrium, $\\hat{c} = 0$ and $\\hat{k} = 0$. Assuming the equilibrium conditions, we can rewrite the first of Ramsey model equations as an equilibrium **value** of $\\hat{k}$, and the second equation can be represented as $\\hat{c}$ in terms of $\\hat{k}$. First, assume the production following, $Y = A \\cdot \\left[K^{\\alpha} (AL)^{1 - \\alpha} \\right] (\\alpha \\in (0, 1))$. In terms of effective labor: $\\hat{y} = \\hat{k}^{\\alpha}$.\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{k} &= (\\delta + \\rho)^{1/\\alpha}             \\\\\n",
    "    \\hat{c} &= \\hat{k}^{\\alpha} - (n + g + \\delta) \\cdot \\hat{k}\n",
    "\\end{align}\n",
    "\n",
    "From the second equation we can also derive the value of $\\hat{k}$ that maximizes $\\hat{c}$; the *golden-rule*.\n",
    "\n",
    "\\begin{align}\n",
    "    \\hat{c}                                   &= \\hat{k}^{\\alpha} - (n + g + \\delta) \\cdot \\hat{k}                   \\\\[10pt]\n",
    "    \\frac{\\partial \\hat{c}}{\\partial \\hat{k}} &= \\alpha \\cdot \\hat{k}^{-(1 - \\alpha)} - (n + g + \\delta) \\hat{k} = 0 \\\\[8pt]\n",
    "    \\hat{k}^{*}_{g}                           &= \\left(\\frac{\\alpha}{n + g + \\delta} \\right)^{\\frac{1}{1-\\alpha}}\n",
    "\\end{align}\n",
    "\n",
    "The model has one stable saddle path that converges to equilibrium. This saddle-path is unstable, if the model is outside this path by the smallest distance, then the model diverges away from equilibirum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. THE CODE\n",
    "\n",
    "The Python code has the following structure. Following the above discussion, after importing the required packages and defining parameters and arrays, the equilibrium, gold values, and _loci_ lines are calculated. The code plots two graphs. \n",
    "\n",
    "The first one shows the _loci_ lines with five sample paths (in red). Each one of these sample paths depicts the unstable dynamics of the model. The second graph plots the stable-path (in blue). \n",
    "\n",
    "## 5.1 ESTIMATING THE STABLE-PATH\n",
    "\n",
    "A common way to estimate the stable-path is with a \"shooting\" methodology. This approach shoots a number of paths until it finds one that is close enough the actual stable path. Let $\\tilde{k}$ and $\\tilde{c}$ denote estimated values of the path produced by the \"shooting\" algorithm.\n",
    "\n",
    "This note uses a **simple** code and methodology to illustrate this type of numerical method. It is not intended to be the more efficient one. The stable path is estimated in two halves. One estimation for the stable path section to the left of the equilirium point, and another estimation for the stable path to the right of the equilibrium point. These two estimations are coded as two different functions: (1) backward shoot and (2) forward shoot respectively. Alternatively, one simple \"shoot\" function could include both initial locations. \n",
    "\n",
    "These type of iterations can be time consuming depending on how many shoots the code is going to perform (which depends on how big a correction is when a shoot misses the steady-state), what the tolerance level is, etc. A more sofisticated code could reduce the number of attempts and save computational time.\n",
    "\n",
    "Since we know that is very unlikely that our \"shoots\" will exactly hit the steady state $(k^*, c(k^*)$, we need to define a success criteria. The success criteria is when the shoots hits close enough to the target. We need to define circle around the equilibrium that plays the role of the tolerance level. Once we find a \"shoot\" that hits inside this circle, the shoots is considered a success and an approximation good enough for the estimation purposes.\n",
    "\n",
    "The distance of a point to the steady-state is: $distance = \\sqrt[2]{(\\tilde{k} - k^*)^2 + (\\tilde{c} - c^*)^2}$, where $\\tilde{k}, \\tilde{c}$ are the shoot estimations.\n",
    "\n",
    "The shooting functions have two inputs. The initial level of $k_0$ and a \"step\" size. The \"step\" size is how big the correction of the initial shoot will be after it misses the tolerance level (the target). Given $k_0$, the function estimates a level of $c(k_0)$ (the first shoot). If the shoot misses the target, then the algorithm performs a correction to the initial level of $c_0$ and shoots again. This sequence repeats until the shoot hits the target.  \n",
    "\n",
    "### 5.1.1 \"FORWARD\" SHOOT\n",
    "The first shoot of the forward shoot function is located one step below the _loci_ function. We know is very likely that this initial shoot is too high, and therefore the path will go above $c^*$. The loop has two break points to define a succesfull shoot. The first one is if the path produces a point $(k, c)$ inside the tolerance distance from the steady-state (the above discussion). The second one is the path produces a value of $\\tilde{k} > k_*$. This break point is in place because the step size (how far the new shoot starting point moves) is a fixed amount without a correction or feedback mechanism. It is possible, then, that after missing the target the next iteration starts at a value of $c$ that is too low. Such path may still be an estimation good enough even if it does not satisfy the tolerance level (for the purposes of this notebook, this is enough). \n",
    "\n",
    "There is a third break, which is a cap on the number of iterations the loop can perform. When building a code like this, this type of break can help avoid a mistake that sends the code into an infinite loop.\n",
    "\n",
    "If the path produces a value of $\\tilde{c} > c^*$, then the loop resets and starts again from a new initial $c_0$ that is located one \"step size\" below the previous shoot.\n",
    "\n",
    "### 5.1.2 \"BACKWARD\" SHOOT\n",
    "\n",
    "The backward shoot function follow a similar structure than the forward shoot. Yet, because the starting point is different, the break points have to be modified. However, there is a difference in the starting point. Rather than starting at the _loci_ level given the initial on step above $c_0 = c(k_0)$, the code starts with $c_0 = c(k_0)\\cdot 1.5$. This is an example of an easy way to save some computational power. Looking at the first figure (the unstable paths) we know that the level of consumption of the stable path will be several steps above $c(k_0)$. Therefore, we can save computation time by avoiding shoot attempts between $c(k_0)$ and $c(k_0)\\cdot 1.5$. \n",
    "\n",
    "The forward shoot loop has two stop points of success. The first one is if the shoot produces a point inside the tolerance circle. The second one is if the path produces a value of $\\tilde{k} < k^*$ (for similar reasons than the forward shoot break). There is also the third break built on the number of iterations the loop will perform. Since the starting point is farther away from the equilibrium point, this break has a larger limit of iterations than the forward shoot loop.\n",
    "\n",
    "If the shoot produces a value of $\\tilde{c} < c^*$, then the code resets and calculates a new shoot with a new $c_0$ one step above the previous shoot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAGKCAYAAADHWq+sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5jdVbXG8e9LAoReDV2KonQEpSkiIoIgYr1cUAFBhCsCKogiVRARCxGRokhvAiIKUqRIsyBF6UFaQi+hJoEQ0tb9Y/2GDCHTz5z9O2fez/PMc5hTF5PkzDp7r72WIgIzMzMz67s5SgdgZmZm1qqcSJmZmZn1kxMpMzMzs35yImVmZmbWT06kzMzMzPqpSCI1h3QjANLZSHuViKGvJJ0s6eTScZiZmVl9DC/xogEdPRcuAZ4qEYOZmZnZQBVJpJaCdwAQcVGJ1zczMzNrhCJbe4vBYgBIw5CWKRGDmZmZ2UAVSaSGwTAkAdsBJyGNKBGHmZmZ2UAUSaQEAuYFxua3LF8iDjMzM7OBKNn+YCEykQJYsWAcZmZmZv3Sq0RK0pqSnpW0RiNetDqytwAwDphEiyZSjf65mJmZWWvp7YrUgcAHq8sBuwvuJeIBIoJclWrJRIoG/1zMzMystShzmSa/qHRDRGxafbMOEETc2fRA+qCjGWdE7N7g5/0E8EtgGHBKRBzdyOc3MzOzwVOkRmp5WAbpQwBE3FH3JGqwSBoGnABsBawG7CBptbJRmZmZWW/1tkZqaUl/kHSHpP9KWn8gL7pQFpqvVj35cKT3tWI/qQb8XNYHHo6IMRExBTgf+HTjIzUzM7PB0GNnc0nDgSuBgyLiMknzkttQ/TYNppHJFGQydzjwe+CcgTxvM/X0c5H0N7KgflbfiYhrq/9eBnii021PAhsMUshmZmbWYL0ZEfMZ4P6IuAwgIiYN9EWnw3Q6EqmIKUhP0XoF593+XCLiw714Ds3muuYXrZmZmVm/9CaReh/wr0a+6CwrUpAn91Zt5Gs0Qbc/l16uSD0JLNfptmWBpxsWoZmZmQ2q3iRSzwJrd3wj6R0R8fxAXnQqTCWTqQ5jgE2QFiBi4kCeu4m6/bn0ckXqNmBlSSsCTwHbA19sdKBmZmY2OHpTbH4GsISk+yTdCWw00Bd9HJ4mYt9OV3V0OF9hoM/dRGcwwJ9LREwD9gKuAu4HLoyI+xoapZmZmQ2a8n2k8ooRZOH140RMbXpAvTBYfaTMzMysdRXpIzU/zIt0GNJIACImE/FIXZMoMzMzs9kpkkgNy9f9ALDYm1dKayN9tkQ8ZmZmZv1RJJGaNrPQvPPJvXWAncj+TGZmZma1VySRmpp9pODtLRCGky0AzMzMzGqvVCI1uxWpjpN7rdaY08zMzIaoIolUZPfuR4E3Ol39FNlfyomUmZmZtYRy9UgRe8/y/XSkx4AlygRkZmZm1jd1K+z+PhGTSwdhZmZm1htFtvYAkL6E9L23XOckyszMzFpIuUQqe0it9pZrpMWQ9kdavUxIZmZmZr1XMpEaDyyIpE7XTQY2YdYEy8zMzKyGSidSw4F537wm4jVgHD65Z2ZmZi2gdCIFb+0lBTAGWKnJsZiZmZn1WclE6jngvtnEMBZYGmlE80MyMzMz672SfaRGAwfM5paHyGRqIbJmyszMzKyW6tZHCiJuA24rHYaZmZlZT0r2kRqOdBLStsViMDMzMxuAcolUxDSyl9TbR8JIuyMd1vSYzMzMzPqgZLE55Mm9WU/tAQhYY5YeU2ZmZma1UtdEaiwwAliyueGYmZmZ9V6dEylwY04zMzOrsdKJ1L3A6Nlc/xgQOJEyMzOzGivb/iDij11cPwXpKuDp5gZkZmZm1nv16yPVIeKE0iGYmZmZdafs1p60EdKFSMt0cfsIpPome2ZmZjakla6RmgLMAyz4tlukNYALgVWbHJOZmZlZr5ROpF6pLmd3cu8psp+UC87NzMyslkonUuOry7cnUhEvk4nWSs0MyMzMzKy3SidSE6rL2a1IQfaT8oqUmZmZ1VLZRCpiCnAV2TdqdsYC73TBuZmZmdVR+QQl4vhubv0n8DwwDJjWnIDMzMzMeqd8IpWDiYcTMfVtt0U8ADzQ9JjMzMzMeqF0jRTAAcAxXd4qLYm0QrOCMTMzM+utOiRSr9F1sTnA94CvNikWMzMzs16rQyI1Hliw2uKbnbHASt3cbmZmZlZEHRKpV8harXm7uH0s2fl8kaZFZGZmZtYLdUiketNLCtxPyszMzGqmDonUI8Dvgcld3O5EyszMzGqpfPuDiMeBs7q5/TWkI8iEy8zMzKw2yidSWUQ+PzCDiNdme5+I25oak5mZmVkv1GFrbzhwHvDJLu8hjUTaBmnupkVlZmZm1oPyiVR2NJ9E972k3gXsASzflJjMzMzMeqF8IpXG030iNaa6dMG5mZmZ1UZdEqkJdJ9IjSNXrVbq9lmkPZHOb2BcZmZmZl0qX2yeXgFGdnlrRCA9Ss8rUnsDKyMtSMSEHu5rZmZmNiB1WZG6Brikh/uMBZbvclSMtBywAvA63RWum5mZmTVIPRKpiFuI+GsP9/od8BUioovbPw8E2Uphp0aGZ2ZmZjY79UikpBFIKyB1vdUYMZ6I17t5lp2Bear/3hSpq9l9ZmZmZg1Rj0QKNgJ+RXd1UpKQvoS06WxuWwJYtdM1U4AtGhuimZmZ2VvVJZEaX112fXIvt/Q+BHx4Nrd+epbvFwR2bEhkZmZmZl2oSyL1SnXZXQsEyILz2Z3cexK4r9P3/wHubkBcZmZmZl2qSyLV0aqgN4nUO5Dmf8u1EVcQ8X7gb9X37yfi8EYHaWZmZtZZXRKpnrf20tjqsvvGnGZmZmZNUI+GnBFTkY4DHurhng8CrwELD35QZmZmZt2rRyIFEHFNL+4zEWmHbnpJmZmZmTVNXbb2QFoWaeUe7+ckyszMzGqiPokU7ELOyuue9B6kXyH1NHfPzMzMbFDVKZEaT8/F5pCtElYAVh/UaMzMzMx6UL9EqquhxDM9D7wArDb4IZmZmZl1rW6J1DCg+xl5WSM1Gli9F0mXmZmZ2aCpUyLV0ZRzwV7c9z5gUWCJwQvHzMzMrHt1SqTuBg4HXu7Ffe8BbiJXsMzMzMyKqFMfqRfI2qfe3PcJ4GeDGo+ZmZlZD+qzIiXNibQB0jK9vL+QFh/kqMzMzMy6VJ9EKrfpDgY26OX9twFOR+pNy4T2I83hYnszM7Oy6pRIvQFMpXe9pAAeri6HXhsEaSTwD+D3pUMxMzMbyuqTSGVbg/H07tQeZCI1laGWSEkbku0f1gOWG+BzbYJ0IdLt1eUmjQjRzMxsqKhPIpUm0NtEKmIq8CBDpcN51oR9E7gOWIyBnljMpGkUMBJ4troc5WTKzMys9+pzai/1PpFK9wFfQBpBxORBiqkuvgf8AJi703XrIzViiPNTwB3Vf+9FtpYwMzOzHtQtkToFmNGH+98EPNbHx7Sq64D9yZWojj+3SeTPrD+2B14j5xYuQyZSrwIrDShKMzOzIaReiVTEY/24f98e06oibkVaHbiU3M6cF7iXiG/26/mkpYBlgXcyc5twfmDMwIM1MzMbGupVIyUth7Q1Uu8TvHzMBwcxqvqIeBb4IHAyA1+FOx7o3IdrATKROn6Az2tmZjZk1CuRgjWBr5O/1HtrK2DfPiVfrSxiGhHfBj4P/HwAz3MTcE2na8YB+1bXm5mZWS/ULfnoGFy8EL2buQfZCuBTDLXanog/Dejx2cxzPmAisDAR2zUiLDMzs6GkbitSHYlUX07uja4uh1Y/qYHKvl37kwXrZmZm1g+tn0hFvET2QRoa/aQaKeJ5YHrpMMzMzFpV3bb2xleXfVmRguwn9f4Gx9LepG3I9gdmZmbWT3VbkRoP7Alc38fHnQ3s0fhw2lTWR32eHDNjZmZm/VSvFamIGcAT/XjciwBIDQ6obS1Jtj64t3QgZmZmraxuK1IgbYa0cT8ety05g856tkZ1eU/RKMzMzFpc/RKp7Au1JdCxBdVbqwOLDkZAbWhNchv1ydKBmJmZtbI6JlITgIWQFgauQTq3l48bDcw1eGG1lfmBu6oWCGZmZtZP9aqRShPIIujRwEjgP7183HxkggDShcDx7tLdhYgjkOqYRJuZmbWUOv4yXRPYGliKmcN0uydtAuwLdGwFjgRGVdfb7GRhv5mZmQ1AvVakpG8BO/LWBGo9pL5uQU2sLvcCvCrVmfR1YF4ijikdipmZWaurVyIFtwGvk4lU52Tqgh4et0X1ODGzfcKrDLX5ez3J4v31gAdKh2JmZtYO6pVIRfwDaQ3gCuDdwDzAbURs3+3jsiZqJDNXoiDrpcYMUqStaiTwDuCi0oGYmZm1gzrWSE0CjgMuofdz4I4nE6cFyFWpBarvjx+MAFvYmtWlG3GamZk1QB0TqbmBTwO/BXamN8lQns7bFxhHdu0eB+zrU3tvsyZ5KrLv3ePNzMzsbeq1tZcmVJcLEdHbHlIdyZQTp+49DDzn/lFmZmaNUcdEahIwDViwdCBtJ+LPpUMwMzNrJ/Xb2svVkgk4kWosaVGkEaXDMDMzayf1S6TSBGDe0kG0mZ2B3/RxfqGZmZl1o45bewDfJmJa6SDazHuBB10fZWZm1jj1XJFyEtVY0pzA0sCjhSMxMzNrK/VMpKQPIX2jdBhtZBne2vXdzMzMGqCeiRQsD2yJVNf4Ws1y1eXjRaMwMzNrM3VNVCYws0O5DdxDwEnA06UDMTMzayd1TqTALRAaI+JZIq4gYkrpUMzMzNpJ3ROphYpG0S6kdZDeUToMMzOzdlPXROoVYCIwV+lAWp40HDgU2Lp0KGZmZu2mnn2kIh4Fvlg6jDaxJPnn7EJzMzOzBqvripQ1zjuryyeLRmFmZtaG6ptISd9C+mTpMNpAR+sD95AyMzNrsPomUrAGsErpINrAcsA4IiaXDsTMzKzd1LNGKo3Hp/Ya4SxgkdJBmJmZtaM6J1ITcAIwcBHjgHGlwzAzM2tHdd7aG48bcg6MtAjSp5AWKx2KmZlZO6pzIvUM8HLpIFrcysDugJtxmpmZDYL6bu1FXABcUDqMFtf1iT1pQbKYf5nq+/WBMUS80KzgzMzMWl2dV6Rs4JYDXiLitdncdihwU6fvbwROaUpUZmZmbaK+iZT0XqSfIb2z5ztbF95J1x3NbwLeYObfgRnAX5sRlJmZWbuobyKVsa0CLF46kJYkidy266oR5zW8dZbhHMDFgx2WmZlZO6lvjVS2PwCf3OufiEDaGZizi9tfR7qOmcOMxxDxVLPCMzMzawd1XpHqSKTclLO/IiYTMbGbe5wFTAQmV/9tZmZmfVDnROpVsm7HK1L9Ia2EtCtSd01Nr2Dm9t5FTYjKzMysrdQ3kYoI4G7gldKhtKgVgc/y1jqot8rVqpuBp4l4pElxmZmZtY0610hBxCGlQ2hh81WXs2t90NneeNXPzMysX+qdSNlAdCRSk7q9V8S9gx+KmZlZe6rv1h6AtAvSD0uH0aLmAyYRMaN0IGZmZu2q3okUzEvW+ljfjaDnbT0zMzMbgLpv7Y0HFkRSVXxuvRVxPNKw0mGYmZm1s7qvSE0ABMxfOpCWFDG9dAhmZmbtrO6J1Pjq0qfKZkNpDknDJM0xy407IX2iUGhmZmZDQt239p4BbiEbc7Y1SXORs/GWA5YAFiPnDC4GLELWi3V8zUP+2WmW55hKDiKefAms+gC89F1pO7IX18uzfL0APAs8GxFTBv1/0MzMrA3VO5GKeBA4snQYjaQcJrwUsCqwGvAeYHlgSd66QhjMTHheBsaRrQwmAa8DU8kEM6rLOYC5ySLzuZeE5R7Kx0yuXm81MiF7W92UpBeAp8nE6kngUWAs8GhEvN6w/3kzM7M2U+9Eqg1UidNKwAbA+sCazJwfOA14BLgHuBx4AnicXIl7Jfpb45SvucH6cMF+EefOEsv8ZEK1OJlgLV1dLgWsDmxOp4RO0rNkUvUw8ADwX+DxcFsFMzOzmidS0pzAqcAfifhj6XB6q9qm+xDwUTKBWqy66XHgRuA+4H7g4UHaVpuX3PZ7tfOVkScfJ1Zfj3cR+5zk9uIKZOuJFclE8APAnNXdJkt6kEyq7iNH+TwZPllpZmZDTL0TqYipSPORKyi1ViUgGwBbAJuSycx4ssbrFuDWiHimSeGMIE88TuzrAyNiKjCm+nqTpOFkUrUK8N7qchtgu+ouL0m6m0yq7gLud+2VmZm1u3onUmk6NT5dKGlp4HPAp8mEbyJwLXA1cFu/t+cGIuJF4EuNfcqYBjxUff0ZoDopuCKwNrBWdblp9ZBpkkYDtwO3Anc7sTIzs3bTColULUlaC/gK8GGy4PtvwCXAzdWqTtur6qQeqb4uBpC0KFkHtjawDvkz2hWYIulOMqm6jVyxcp2VmZm1NCdSfSRpA2A3MkkYD5wGXBwRzxUNrDNpVeCzwClEjGvmS0fES2Qd2I0ZiuYD1gXWI4vt96ruOlHSbcDfgb9XjzMzM2sprZBI3URuJ5W2MHmqbV2yrcAxwJ9q2h5gSWAj4IzCcRARr5GrdX+DN1esPkDWk30Q2Ky6fjRVUgX816tVZmbWCuqfSEUcX/LlJS0GfIs8gfcG8ENyBarO9T4dsc1VNIrZqFaergaurtoxrAxsXH19DdidLFz/O3AD8K+a/6zNzGwIq38iVUhVSP0FYE+y0eUDwIMRcX7RwHpncnU5omgUPajaJTxYfZ0maWFylWpjcqVqW2BSlVRdD/wjIiaVitfMzGxW9U+kpHOBa4k4vXkvqWWBI8iTaLcCRwMHN+v1G6Aj2ZinaBR9FBGvAFcAV1TtJD5ArgRuSraVmCLpX2RSdWNETCgVq5mZGbRCIpVNIN821mQwVFtNnwb2I9suHApcGRGRN7WMV4HnyNOELak6+XgzcLOko8lTgJtVX5sAUyX9E7gKuCkiJnf5ZGZmZoOkFRKpppA0L5k4bU72PjqsVifx+iLiCfJkYVuoCs/vAO6QNIqcG7gl8HHgI8Drkm4A/gLcUvW8MjMzG3ROpABJ7wR+To5F+RVwtk+N1VNVV3UfcJ+kY8lTlFsCHwO2AsZL+iuZVN3pP0czMxtMQz6RkrQR8GNyK2+viLi1cEgDJw0DDgeuJ+KvpcMZLFWSdDtwu6Sfki0ftgS2JrvNPy3pMuCyiHi6XKRmZtauWiGRuooc8Ntwkj4JHAY8DHynbX7ZRkxHWp08DTckVDVVNwE3SZqHLFD/FFVLBUn/Bi4Frqtp7y8zM2tB9U+kIk4djKeVtCPwTXJcyXeqxpHtZDItdmqvUapE6UrgSklLAp8kk6rDge9JupacF3hntVVoZmbWL/VPpKThwAwaVOtSncz7Ojn/7Vrg0DZt+Pg6QzSR6iwingVOlXQa8D4yodqc7FH1pKQ/An/2iBozM+uPOUoH0AvnkUlPo+xWPd+fgAPbNIkCJ1JvEemOiDgC+AS5pTsO2JvsW/VjSeupxfpcmJlZWfVfkWogSTsBewCXA0e1+YmuMcDE0kHUUdUd/XLgckkrkgOetyHbKTzeaZXqlYJhmplZCxgyiZSkzwH7ANcAR7R5EgURx5QOoRVExFhglKQTyBYKnydr5/aUdB3wB+AO11KZmdnsDIlEStKGwAHAP4BDImJ64ZCsZiLiDWaOp3kX2T7hk2Q7hYclXUB2uXcHdTMze1Mr1EgNSLV18xNyq+v7Q6brtfS/SIeXDqMVRcQjEfEzspbqh8AM4CAyydpH0tJFAzQzs9pohUTqEuCu/jxQ0sLAscAbwLer2pihYhFg5dJBtLKImBwRlwBfIvtR3Qp8GfiTpGNcnG5mZvXf2os4tz8PkzQHcCTwDmCPiHimoXHV38vAAkgj8HbUgFT1UR2z/pYg66g+R875GyPpfHLbz40+zcyGmPqvSEkLIM3dj0d+GdgQOCYi7mlwVK3giepyuaJRtJmIeC4iTiTH0BwOTAEOZOa238iiAZqZWVPVP5GCU4Ed+/IA5XiUbwDXARcPRlAt4LHq8p1Fo2hTETElIv5M/t3cFfgXmbz/WdIRkt5TNEAzM2uK+m/t9ZGk+YCjgOeBI4fwsfVngDvJxpw2SKq/X3cDd1dF6NsDnwG2lnQrcA5w8xD+e2hm1tZaYUWqr74NLAUcFBETSgdTTMQMIg4h4p+lQxkqIuLpiBhFtk04Dlixujxf0raS5ioaoJmZNVxbJVKS1iVXA86JiH6d9Gs70rDSIQw1ETExIs4i5/kdRrZPOJTc9ttV0oJFAzQzs4Zpm0Sq+rR/EPA08NvC4dSDtAXwe3K705osIqZGxOXAF8mavYeAPcnRNN92YbqZWetrhRqp88lmmj35CrA8sI+Pob/pZWBOsuD8/sKxDFlVfdQtwC2SViYL1HcAtpN0GXBWRDzR3XOYmVk91X9FKuJiIu7s7i6Slgd2Aa4O1wR19nh16ZN7NRERD0XEoeQW9CVkPdUfJB3lk35mZq2n/omUtDjS/D3c65tkP5+fNyGiVjIOmEyu1FmNVIXpRwOfIk/2bQycJ+mXkt5XNjozM+ut+idScCJ5pHy2ql86mwBnRsRLTYuqFeSW0uN4Raq2IuLFiDiOXJk6EVgdOEXSbyV90CNozMzqrRUSqS5Vv2T2Bl4Eflc4nLq6Gri5dBDWveqk32nANuTK6tJk64SzJG3ihMrMrJ5aOpEit0PWBk52gXkXIq4iT45ZC6gGJZ9P1lAdASwAjALOkbSpEyozs3pp2USqGkq8F7l1dUnhcOpLEtKivagzsxqpWidcSg5I/gEwL7lSda6kzaq//2ZmVlgrvxlvCbwLOCkippUOpsYWBc4k68isxUTE9Ii4DPgC2dRzBPBTsjB9cydUZmZltcKb8OnkQNg3Vb88dgEeBv5aIqgW8hIwEXhv6UCs/6qE6goyoTqY7AF3NDl+ZgsnVGZmZdT/zTfiSiLuneXajYGVgDMiYkaBqFpHnty7A3g/rq9peRExIyL+AmwHHFhdfRRwgaRPOKEyM2uu+r/pSsshLTLLtV8hR8Fc0/yAWtJtwELAyqUDscaoEqqrydYg3wcCOJJcoXINlZlZk7TCm+0xZMEtAJLWANYCzo2I6cWiai3/Jn/Rrlc6EGusKqG6hpkJ1RxkDdXZkj7kU35mZoOrFRKpWf0vMAn4c+lAWkbERODHwJWlQ7HB0Smh2g44jGyb8EvgVElOoM3MBklLJVKSFgM+DlwaEZNKx9NSIm7Gnd/bXpVQXQ58jqydWhI4SdJJktYqG52ZWftpqUQK+DR5WunC0oG0HGk40ifwHLchISKmRcTFZGPPY8hWIadVs/xWKRudmVn7aJlEqiqe/TRwe0Q8XjqeFjSd3BbdunQg1jwRMSUifkf+2zkeWJPskv5TSSuVjc7MrPW1QiJ1EvA3YF1gGdzFvH+yDcLtwDpIc5YOx5orIl6PiDOAbYGTgQ3Jlgk/kLRk0eDMzFpY/ROpiOuJeIDcopgIXFc4olZ2O9kZe7XSgVgZEfFqRJwMfAo4F9gCuFjStyUtVDY6M7PWU/9ESlr5cml54KPAVRHxRumQWthdwFTcBmHIi4jxEXEsWZT+F2AH4FJJu0qap2x0Zmato/6JFPxoKuwLzI2P7w9MxGTgbmDx0qFYPUTEsxFxBNmH6t/AnsCfJH1O0vCy0ZmZ1V8rJFL8F1YHngXuKR1LGziSiKNLB2H1EhFjImJf4KvAk+T4mQslfcxNPc3MulbvREraZDp8cCfY+C5YJnLGng1ExDQAz92z2YmIu4DdyFXgacBPgDMlfaBoYGZmNVXfREraBBg1FeafCDNWhCnAqOp6Gwhpd+CQ0mFYPUW6idzuOxxYDPi1pF9Jek/Z6MzM6qW+iRTsBbw6GYYPh6nzwzjg1ep6G5jxwHpIy5QOxOqr6pL+Z+CzwLHkFvt5ko5wywQzs1TnRGqlqfD6WJgBPFftQ70KuIngwF1NbttsVToQq7+qqec5ZFPPM4DNyZYJ35A0X9HgzMwKq3MiNeYVWOplmLFgFpoDzA+MKRlUW4h4GfgnsDnSiNLhWGuIiIkRcTzZMuFaYBfyhN8XJA0rG52ZWRl1TqSOfx1GLgSxGEwgp9nPT465sIG7HJgP+EjpQKy1VC0TDgV2AsYCB5Bd0j/sE35mNtTUNpES/H1veDrg6Tlygv04YF+yCNYG7n7gFLJ3kFmfRcRoYA9gP/K95BfASZLeWzQwM7MmqnPDvfdcClwK+0fEFaWDaTs5e89zC21AIv8e3SjpH+SW3x7kUOQrgBMj4rmiAZqZDbLarkiRQ1UBbikaRbuT1kH6bOkwrLVFxLSIuJCciXk2OcPvj5L2lDRv2ejMzAZPnROpDYCHIuLF0oG0ufWAnfDAWmuAqiD9OODzwPXArswcOVPn9xszs36p5Rub8iTZ+/BqVDNcTm7xblk6EGsfEfF0RBwE7Aw8To6cOdcd0s2s3dQykQLWBebEidTgi3gKuAvYCh9htwaLiPuAr5En++YnO6T/TG4Ga2Ztoq6J1DpkI847SwcyRFwOLE5u85k1VDVy5lrgC8CJZP3jRZL2dv2UmbW6uiZSawP/jYjXSwcyRNwKjKbepzitxUXEGxFxGnm672py2+9PkrZ1/ZSZtaravXlJmgtYA69GNU/EdCK+R8TfS4di7S8ino+Iw8hE6kngUOBMSe8rG5mZWd/VLpECVgHmwolU80nDkdYqHYYNDVX91FeBg4HFgFMkHeWByGbWSuqYSHV8KnUi1XzbAUciLVc6EBsaqvqpv5DtEk4mRxZdLOn/JM1TNjozs57VNZF6IiJeKh3IEHQ5MJnccjFrmoh4PSJOJuunrgd2A/4gaUvP7zOzOqtVIlW9Ya4B3F06liEpYjzwB2ADpFVLh2NDT0Q8V/Wf2hV4EfgR8BtJK5eNzMxs9mqVSJFH8BclB+paGZcALwG74JUAKyQi7iZXRn8EvIts5rm/pAXLRmZm9lZ1S6RWqS7/WzSKoSxiMvA7YG5ggcLR2BAWETMi4o/AZ4GLgP8h66c+43YJZlYXdXszWhUI4MHSgQxxVwPfImJC6UDMImJCRPwU+DLwKHnK7wxJqxcNzMyM+iVSqwCPRcSk0oEMaREziAikBZDeWzocM4CIeJAcN3MwMJLsPXWopMssSuAAACAASURBVEXLRmZmQ1kdEylv69XHfsCB5BBps+I6tUv4HHAWsDW53be9PCvSzAqoTSJVfaociQvN6+QCsvh/29KBmHUWEZMi4jhge+Be4DvAeZLeXzYyMxtqapNIAR1bSA8UjcJmirgfuAX4PD4tZTUUEY8Ce5OJ1Dxkq4SjJI0sGpiZDRl1SqRWqi4fLhqFzepM8hfU/5YOxGx2qu2+G8hTfScDm5LbfTtJmrNkbGbW/uqUSK0IvBIRr5QOxDqJeII8xbew+0pZnUXEG1V39C+QK6n7kP2n1i0bmZm1s+GlA+hkBWBs6SBstk4iYnrpIMx6IyKeBvaT9GFgf+BkSVcAv4yIF8tGZ2btpm4rUk6k6qgjiZKWRVqncDRmvRIRfyMHcZ8KfJyc3bedm3maWSPV4g1F0iLAQmSzPauvbwD7k39eZrUXEZMj4iSyxm808F3gLDfzNLNGqUUiRa5GgROpujsRGAHs6XopayUR8Tj5QeD7wGJkZ/QDPbvPzAaqLonUCtWlt/bqLAvPzwE2BDYpHI1Zn1Sn+64hi9F/B3yGPN23rbf7zKy/6vLmsQIwGXiucBzWsz+Rvb7+D4/msBYUEa9FxCjgS8DjwKHAbyWtXDYyM2tFdUmkVgQejYgZpQOxHuSf0S+Aa4BXC0dj1m8R8RCwG3AEsDzZKmFfSfOWjczMWkldEqnlyE+G1goiniLiNCKmuFbKWllEzIiIS8nZfX8CdiBP920m/902s14onkhVtQlLAM+UjsX6KLdCjvEWn7W6iJgQEUcBuwAvAz8FRklaqmxkZlZ3xRMp8gTNnDiRakWvkVsie3llytpBRNwL7EhuX68H/F7SlyUNKxuZmdVVHRKpjk98TqRaTXaQPov8hbNZ4WjMGiIipkfEueTsvluBbwFnS1qjbGRmVkd1SKSWrC6fLRqF9delZKPDryEtVjoYs0aJiGeA/cgxMwsDp0v6rqT5y0ZmZnVSh0TKK1KtLCKAY8nt2W0KR2PWUFXvqevJ1akLqsuLJG3uYnQzg/okUhMiYlLpQKyf8pP7/sDZpUMxGwxV76mfAzsDLwJHA8dKWrpsZGZWWl0SKW/rtbqIMUTMQFoM15JYm4qI0cBOwChgXbIYfSdJw8tGZmal1CWR8rZe+9gbOBhpmdKBmA2Gqhj9PHLUzL+Afchi9DXLRmZmJRRNpKoagyVxItVOTgRmAAfhDtHWxiLiuYjYD/gOsBBwmqQDJC1QODQza6LSK1ILAPPiRKp9RIwj60eWAfZ1fylrdxFxAzMHIX8OF6ObDSmlEymf2GtHEXcDpwIbAJ8qHI3ZoIuISdUg5J2B58kPE8dIGlk2MjMbbKUTqSWqy+eKRmGD4c/AycBfSwdi1iwRcT+ZTP2S/CBxkaQvVKOwzKwNlf7HvXB1+XLRKKzxIoKIPxPxGtJc+JO5DRFVMfrZwP8C9wIHAL+VtGLZyMxsMDiRsmY4ADgCab7SgZg1S0Q8CXwDOBxYEThP0m6S5iwbmZk1UulEahFgCjC5cBw2uC4mT2fu5+JzG0qqzuh/JovRrwf+DzhX0lplIzOzRimdSC0MvBw5ZsTaVcS9ZL3UesCXCkdj1nQR8VJEHAh8G5gfOFXS/nKLELOWVzqRWgR4pXAM1hxXAtcA/4u0UelgzEqIiL+Rq1O/B7YjO6NvXDYqMxuIOiRSro8aCnLV8STgJmBs4WjMiqlaJfwU+CrwGjmz7yhJixYOzcz6oXQitTBekRo6IqYS8TMinkUS0lI9P8isPUX2W/sy8BtgM7JVwjZu5GnWWkonUl6RGrq+CByL9K7SgZiVEhFTIuK35L+HR4EfACfIsyrNWkaxRErSXOR4GK9IDU1XkdsaRyAtVzoYs5IiYgywG/ATYA3gQkk7ShpWNjIz60nJFSn3kBrKIl4ADgKmA0ciLVk4IrOiImJGRPwe+B/gFuCbwOnyqq1ZrdUhkfKK1FAV8QxwCDAncBj+9G1GRDwH7AccCCxN9p36mht5mtWTEykrK+Ix4DDgN0RMLx2OWR1UjTyvJlenrgP2AM6StGrZyMxsViUTqUWqS2/tDXURDxFxJwDSBh4lY5Yi4uWqked+5HvmmZK+UdWYmlkNeEXK6kNaDPgeuc03onQ4ZnURETeSDTwvA3Yh5/Z5zIxZDZRekQpgfMEYrE4iXgR+DqwCHIQ/dZu9KSImRMQRwN7ACHLMzL6S5ikcmtmQVnpFakJEzCgYg9VNxD+BXwLvA76LNLxwRGa1EhE3k6tTF5H9p86X9IGyUZkNXaVXpFwfZW8X8Vey2/MGwPqFozGrnWrMzE+A3YEZwK8lHSjXF5o1XekVKddH2exFXAbsV61QgcdmmL1NRPwH2AE4B/gMOQT5g2WjMhtaSiZS8wCTCr6+1V3EgwBIKwE/RxpZNiCz+omIyRFxLDOHIB8n6QeSFiwcmtmQUDKRGkZ2tTbryQhgGeCnSMuXDsasjiLiHuBLwKnA1uQQ5I+Wjcqs/ZVMpIYD0wq+vrWKiNHAAYCAnyCtVjgis1qqhiCfBOwIvAD8TNKPJS3Sw0PNrJ+8ImWtIeJRYH+yru6HSCuXDcisviLiAWAn4ERgU3II8mZFgzJrU6UTKa9IWe9FjCMbdl4OjC0cjVmtRcS0iDgN+DLwLPBTSUdJWriHh5pZH5Te2vOKlPVNxHgiTiNiGtJCSFv5RJ9Z1yLiEbIb+knAZuTqlGunzBqk/itS0ue9jWNd2ArYE9jNyZRZ16rVqVPJ2qnnydqpH0laqHBoZi2v3itS0txkf5T9mxGQtZwLgEuAbYHvuAu6Wfci4iFgZ7Lh7eZk36mPlI3KrLWVTqR6WpH6GHlS63+QSsZqdRQR5FHvM4FNgMNx7xyzblWrU79l5sm+YyQd4b5TZv1T90Tqy8Bc5DagR4XY20UEERcBo8hu+Z7daNYLkQ1vdwZOBrYkV6c2KRuVWespXSPV9dZebtN8ilyRmhfYvjlhWUuKuB7Yh4hXkebEYzLMehQRUyPiZLJVwkvAKEmHe3XKrPfqvCK1CTNXF4YB27ug2LoV0ZGYfwL4PtI+SHOVDMmsFXTqO3UKeYjjAkkbl43KrDXUd0UqB3HO3+n7+YG1BzUiK0P6QfXVqL+PlwMXAh8HjkZ6R4Oe16xtVatTvya3+yYAx0o6TNIChUMzq7XSBdzdrUitD7zR6X5zA+sOekRDibQJ0oVIt1eXpeojtgUOBK6hEaMsImYQcTZwJDmj71ik1Qf8vGZDQETcTxainwZ8klyd8la5WRdKJVIdW3TdrUh9jJmJ09XAmsBZgxnUkJJJ0yhgJNn1eCQwqmAyNSfwIWA00joNecaIW4B9gWfIT9hm1gvVzL4Tga8ArwLHSTpU0vzdP9Js6CnVd6cjkep6RSriBeAFsizqGSL+24S4hpK9yDfIucgTOwICuBZpcpNj6dg6mBtYEvg70g5EXDrgZ454Cml/IqKqsdsKuJ6I1wf83GZtLiJGS/oy8DVyy29DSUdExL8Kh2ZWG6UTKY+IKWclciVqU97+92DOpkfzVh0nNRsj+00BrADsAXwK6SginmjYa5i1qYiYApwg6QbgcOB4SRcBvwx/IDErvrXnocXljCEL+Dv+DtwBjCO3wX4KLEeEmvKVrw359+FlYAsizm/4/3HEWOBgcgVsFD6VZNZrEXEf8CXgXODzwO8k+QCQDXlOpIau48lEquPPYgLwFHAcmWhkTZG0RJNaCEwF7gJWI+Lvg/YqEfcA3wIeA76H9NVBey2zNhMRb0TEL8iV3TmA30raR24zYkOYE6mhKuImshC7o1fXOOBbRBwNfIOICVVN0feAU5G2Y/COQV8P/BrYiIhnB+k1Zsr6uwOA84F7B/31zNpMRPyHbJJ8Cdl/6mxJ7y0blVkZmlk+0sQXlW4m63AOi4jLe7hzzlOL2K0ZsXUdhk4GiIjdS8bRcNLjdGzjvf02kaclP0+eoJxMnqC8hIhxzQxz0ElfABYFziTijZ7ubmZJ0oeAQ4BFyHEzZ8TM5rhmba/OfaSstJxjdzcRhwF7AzeTfWWyLUV7dZpfgBxJ9Ev8ydqs1yLiH8B2wF+BrwOnS1qhZExmzVR6a8+fWlpFxKNEjAJ2A66rrt0a6UdI67V8UhVxOnAQ2Q7iZ0g7VvMezawHETEhIg4Evg8sC5wnaQc1blqBWW05kbK+iXiBPA4N2Xl+GeBQ4ASkLVt6tl3E3eTK21+B/wFWLBuQWWuJiGvI1albgf2AkyQtXTYqs8FVOpHy1l4ri7iWXKEaRZ6624ss4m5dEa8R8Uvg60Q8BID0vgbOATRra5GHOb4NHAGsCpwv6dNq9VVrsy6UTqS8ItXqIqYRcT3ZUuAg4CIApIWQvo60TMnw+i3iKQCkFYEfksOPlyoak1mLiHQpebLvfrIY/ReSFi8bmVnjlU6kWm9FSpoDaXGkd3e67r1I70daueq7NG/L1wz11czC9NHVNe8FPg78GukQpDVa9GfyKHAM8E7gV0jbIg0rG5JZa4iIp8kC9GPIQfQXSvp42ajMGqvkiJigjolU/pJcFBhJdvIFacuz4QOLwzzAxcAwYDLSdtX4kc8DG83yTM+R214g7QQsBUwkG11OAJ6rhuq2p4hbkXYhT/l9Evgx8BDS94iYWja4Psg/3xuQ7iW3Lr8GbFz9fzS/d4hZi4mIGWQX9JvJ7b4fS/oo8JOIGF82OrOBK5JICRR1SqSkjwEfA5YAFqdjpU7anojXgPnmgmGPwHgykRpXfXU4ubp+QfIY/YLMbHRJdd0KnW4X8CBwS/U6B1av+Qg5uuUR4MWW/0Wdb5LnkXO5NgOWfTOJkj4I3EPExIIR9l7EC0iHAxsC83Yagrwg/mVg1qOIeFT54eorwO7A+6sByP8oG5nZwBRpyDlMunNGJlE7R8eqT1cGqyGnNBJ4nYiJSOuTJ02eYWaS9BxwX8cv/oY15MxfvvMDcxHxYnXdPmRR5jLM3Pa8vmo3ANKGwBPA0w1PrrpryDlYsk7iNLJA/XrgMiIebdrrN0p+qv46OXvsciLq8cHArOaqLuhHAO8C/gAc6wHI1qqKrUhV/9n8XzzZKO5zwEeA3wPnEHEreVx38GUiNHGW646rYhtBHrlfiY4VL2lBsogbcjuxY8XqRiIeaErMjZarO/uQDTA/CmxZbZ2dSMQTZYPrkweA0eQW7hZIJxHhkTNmPYiIByTtCOxJDkJeT9Kh4X8/1oKKFJtHiVN70qpIhwK/IuuZLgX+0rTX742IyUTcT8TlRNxWXfsa8E1ymPC15M/u42TxM0hLIe2P9AmkpVumoDsbfP6KXOY/gxwv0TEoebkqgay3LKQ9HDgSGAH8GGnXskGZtYaImBIRx5KrunMDp0naXW6Eay2myF/YGVkfBYO9IiWp01bYZ4FVmLkN0yq1OdPJuqkxb16XPY06kuCRwBrAJtX3LyLdTa601X8eXv45/AHp4k5/VnsB70a6kdz2G9P1ExSWMd+CdCfwBeBhgDe7onu7z6xbEXG7pO2B/cnaqY0lHRwRjxcOzaxXSg0tvhuYAnw2etrK6U+NVP4S+zC5hXc0EU8hLQpMImJyP2Ou79DiXIVaGlir+loD+AYRE8ijxqsD9wB3E/H8LI9tfo1UT3L79ZNkgfpc5PbZBeTE+daQvxg2AX5ddUw3sx5UrRG+T65Q/QL4Q7T6oRtre0VWpOaGud7IRKqxn9YzodiCLBwfCTxOnpR7ioiXGvpadZJvNE9VX1fOshK3MNm/5WMASM8Ad5C/4Ov5BpWF5ycgnUluY25DJor/qUbQzNMCJ+UeBjYHfoT0T+BsIp4sHJNZrUXENZLuIsdOHQBsUp3se6FwaGZdKrIiNY/0wOQsuN46etp+6suKlLQdsCPwX+BC4PZGJQu1XpHqSSaYK5CrVWsCcxBxRHXbS2R90lrA/bXciurYyoyYhvQJYA/gn8BVZAuFeiaEmfR9juwzNjf59/iSskGZ1V81TuZ/yIkJrwM/iojrun+UWRnteGrveuAXtf3lWkL+LMZWXzN/keebVUeX7qPIU4F3AVd1KnYvLxv6dfTluge4ktz22wR4Gukq4I+1+zPP4c7nI11J1k9l13dp4er2V4rFZlZj1XbehZJuJUc0/VTSZcDPI+LVstGZvVWRU3tzzHzdxiRS0twARFyIk6jey59TxxbZkWQSuiKwLADS/Ei7V+Nv5i4T5CwiniLiZGAncljyK8CH3vwzzxN/9an3gmxMGnHqm0OQM/ZTkHZEmq9kaGZ1FrnNvwtwCrA12SF93aJBmc2i5IgYaET7A2kdYF+kHxDxiJOofspxNbfMskq1ArAl2e9pKtJ9wL/JZqFla5Ryted64Pqq/xZICwC/BF5Auhq4tqarPn8gt/q2Az5ZdX6/rL8HIczaWWS5wa+VtYZHAL+RdDbw68j3AbOiiqxITZu5TTOwFansjnsguSrx7ADDMugYPjyt+u97gR3Iws8rgMWAr5Kd2TuGNX8Qad4ywVZmJiCTyUTqRWBn4AykA5CWKxbb7OSq2s/I/mCjyVh3KhuUWb1Fnn79IvBH8t/Lmeo8PN6skCIrUm/MXInq/4pU/nI8DHgZOLSaiWeNlp/47qi+Osa7vFjduiV5qm460mhytep2Ih4rECnVOJ8bgRuRlqni24yOxF1aCnijNic4sz/WEUirAtmWQnoXuRJ4fVUbZmaViJgEHCXpb8AhwNmSTgDOC/97sUJK9ZF6gDy1t16PPUJmd2pPWgw4hlxR+y4Rg74a1dKn9rozkD5S2a9rFeD9wAfIBOApIv6vun1l4ElKztCShr+5wiYdQHa1z+1JuIW6bQ1Ie5DtHp4ELiJHAdXvJKVZYZIWAQ4mx339GzgsmvC7wGxWRRKpRaQnXoFnImL9Hu88+0RqTnKswKXNGnbrRKpXz7U4sBgRDyANA84hR6eMBm4n3+yeKFbHJi1NrqB9lNymnET+HTq3SDyzkzVqG5FbqiuQq3/nE1GvcUZmNVC1SfgU8B1yYsZPgCvdxNOaqdSpPdGf+ihpXqQFiJhKxHHNSqLakrQL0kPActX3D1WzCPsv4oVOg5SDbKlwKbAQsCtwApkggDTszSLxZol4mogzq1gOBv4FTK3iGY60XbX9V07WqP0T2Ifcun6KbKqa/bSyQ7+ZkW0SIuJS8n3lYbIY/Si1wqxOaxtFaqSq9gd9S6SyueFBwAJI365m0DWWtDfwVyJGN/y562cROgYfp5WAJRr27FmvcE/1dXq1WvV+oKMFwKrADzvVVv0HeKwpq1UZ213VV4d3A18GdkT6L7n197diMxnz5/Afspt7x2rhRsB3kK4n+2Z1P17JbIiIiKckfY0cgr4HsLakQyLi32Ujs6GgyNbeEtKz42B0RGzW7R2lTcji4RfIU3kTgIOIuGFQApMmk4XJXyXid2+9qc229qTlyQ7wHatCE4Gtifh7k15/aeATwLrA8tW1LwHfJ+LpWcbcNEfW3n2EHKfzTjLZ/1ax4vlZSUuQndI3J2cQ3gL8gYj7i8ZlViOSViP74i0HnEW2SZhaNiprZ0USqaWkcc/CXRHx8S7vlEnUKHIV42UywXmZTHJuGpTAMpGam6ydORvYp6MYue0SKYAs+n9P9d14YNEiJ8VytWodcnzNcdUomK+Qw5f/Q65YPdS02HIFaEXgg8C5RATSDuTcxhuBB4r2K5MWIoc6b0MmwP/n/mlmM0maB9gX+Cz5gfHgcCmIDZIiidSC0lMT4e6I2KrLO0kXkoOHP0KuDLwAvEEmORcPUmgHdfrvSeQ21JZEPNemidRB5BHiuYAzidilcEQz5Uy9zclET8Cr5FbbiYXi2YMciD0X+XfxH+SJuoe6fdzgxjQCWIKIx6rO8z8EbgBuII+Jmw1pkjYl3+NGkB/ML3YhujVaqfYHY4H7ImKbbu50O7mdtzFZrFzCdOBjRNzYponUe8j+UNOBzxFxbeGI3i67la9DbgFOqsbDgHQI2SLg3+Sw5cFfus9xLusDHyJXSq8n4rhqBevdwMOFTyR+F3gX2Zj0BuDKqleV2ZAl6R3kwY0NgZuAH0bEy2WjsnZSJJEaLj0+He6IiE93eaeZK1Kdi30XAMYRsd2gBDZza286uSK1Xcex87ZMpACkx4DFgYWbkow0Qq7EHAKsRh6YmALcS7YyaE5xaSZVI4h4EWklsqP68+RK1d+BBwvUeAlYmaw9+wi5era3T7faUCdpDmB7YG/yd8oPIk/HNuv11wSuATaPnBhhbaRI+4NlMyHq6dTe8eQokgXIrZ0Fqu+PH9zomExu6a01RHr3jAJ+0zJJFORImIiDyCPPPwSuAt7BzNE1yyB9F2kLpJGDFMNrRHR0eH+G/DmOJeuWfg6cirTCoLx21zEFEQ8ScRwzhzpnobz01WoA9Tu7ewqzdhQRMyLiPPLfxcvAcZK+o+YNYz+QrLk8sEmvZ01UZEVqRemVR+Hq6GllKQvO9yKP5o8Bjh+0QvN8vSfIY++7zzpAtm1XpNpJx0m/HGT9LaCj59IzZKuD3w36eJhcqdqA3EY4hog3qnqvpcnVquavVGVce5GnEYcD9wF/Af5Zu87uZoNM2Upnb/KD2BjgoChZ62gtr0gitZI0fixcHhFfbPqLd0eas6uVGSdSLSa3uZYF3gesDawO7EbEa0gfJdsb3EnWVw1uMiHtRp6yG04O2L6NTGJuH9TXfXscC5LJ1FbAUmQvqtOaGoNZTUjaCPgBeRr3V8D5ntdn/VEkkXq3NOER+FNEtMzEeydSLU6a4832CZnYbAMMI+ur7gf+TcQfB/H15yPnEa5fXd5HxBHVbRuTCd2LXT9BQ2MRmVw+R8QzSGuRWx43kCcjxzclDrPCZpnXdwtZO/V8g19jc2DniNhxgM/zCbIWcxhwSkQc3Yj4bOCKJFIrS68+DBdGxK5Nf/F+ciLVZrJgfQ1mrlhNqOquOlodjCe7sj/Y8PqxHPa8IBEvkW/kZ1W3PEy+md8KjG3aFqC0LrAzuYU+g+zddQPwDw9MtnZXzev7LNl36g3gyIi4voHPvx85zWbUAJ5jGPAgOSv0SXJVe4cYGlM4aq/IiJiX8y9r40e8mPVW1sDdXn11JDdUw5ZXIdsZAEypRsZc2bCu75mcdNRqvQJ8g1yp2gD4IvAl8pPntdVopBjUwwARHaNolgc2JT+df408fUg1f3DcoIxlMius6it1saR/kx3RfybpT8CoaEw/trWB06vC9t8AT5N1WX35oLQ+8HBU7UwknQ98mhwIb4UVSaRezO0Uf9K1+uhYeclk4dtV/6rVyG7ra5KzCUFamOzXdC+5YvXAgGqs8s308erroqpr+XpkfyyAjwK7Id1JrhTdQcSz/X697mN5DDgT6SxgJBHTq23Ao4DhSDeRK1Xl+mWZDZKIeEzSrsDu5My+dSUd3IBVn7WBceTp4lMi4pyOGyT9jTyRPqvvxFv7+i0DdJ6t+ST5wctqoEgiNQKGTfaKlNVZDiu+pfrqqCuCTKjmIXvS7ABMrUbt/JaIMQOeEZj1SZ3fQMcA15HJ1YZVLM+Q/aHeGJSZhPl8z1XfCTiZXKnaGtgWeArpdCJuaejrmhVWzeQ7QdLNZGuV0yX9BjijP4XokuYEVgB+B+wRETfP8nof7u1TzS7cvsZjg6NIIrUSLDDaK1LWSjqSlYix5IrVfLx1xer16p5bIG1LLrnfX3092+9kJ49lP4T0a7KFwrrA0kS8Ud3jwGoVK1ercrWocR9S8pfHzcDNSPOTvXA2JWupQFqOHOVzM6VnEJo1SET8Rzlf8wBgT2AjSYdGxDN9fKrVyHqmRZnN4kEfVqSeJIcwd1iW3CK0GihSbL6GNPk+OCEi9mv6i/eTi82tV6QNyPYCqwDzVde+DHytWkFaFJjYsJonaTtyperd5KfW14DLiTi7Ic/f8+tvBuxDniR6iVzBuxm42zVV1uqqQvStyIRqBnBURFzdh8fvRP7bPAf4I7BZRDzX/aNm+zzDyWLzjwFPkcnZFyPivr4+lzVekRUp5Ru+V6Ss/eR21y3VVuDyZEK1VKcVpL2A9yE9RK5a/ZdsfTChn693IXBhVdO1Nrlile0Lsrh1FNmA817g3oY3JI24DukWcutxI2Az8s3+S8B0pGWAF2dtcGvWCqqC8CuUNYo/Ao6StCHw814Woq8N3BIRD0r6HnChpM2jjx+kImKasqnuVeSHltOcRNVHkRWptaQ37oFjIqJl2uV7RcoaQvoA+ea6KjlgeDiZSH23un0Tsj5p7IAbhUqLk9sSawIjqmufBk4m4t+DUl+VpwxXIOLB6vtRZEJ5J/Avsl/X4HaXNxsE1arQ14BdycLvgyLi/rJRWR14RcqsmbKbeUfLhbnIIcNzVN//f3t3H1x1feVx/H1ugAAiKD6gVCwWKq5aY/GpxW1rO7p92LXTmd0dVy3VrSO21rViu1tBXbdia21ZO9IqKmIZVlTsUCyuum53tu7oYkW3IvWhgiCgYEwg5fkpCWf/ON9f7kViSC4hv9zk85r5zY353dz7TQbh5HzP95y+wETi/8tmzFYS6fxncV9SxnutA25OLR0+QvTNOpksYwWfwOwyiicQXyHaHOxPsfyutObML4hM1SeJI9xg9hju96aP+6hXlVQCjz+n0y0ysLcQheh3AnPUEb13yyUjNcxscx38xLPOzhVAGSnpEmaHAccTAVb2+Aju81JR+XVEoLIsPdaXHfiYnUz0ojmJYsHrOmAi7htS09KdnZK1iq3OkcTW4zu4P5/qxe4BlhDF8r+n48W8Il3OYtzSDcRW9vPATR6/uEgvlEtGqg4aUfsDkb3FmJjn0pUFIH3T3SFAP6IFQfb/7kbMpuK+OJ0kHAisa1fw4x51U/EexxLZqlEUM1bfAM5I7R3+mK5luG9v7eX28V4OvJWuTBXR2uE0itmqd4Gfoi0T6cbcfVOqefoK8F3gYTP7vrs/k/PSJAe5BFKDoc8mXwienQAADlpJREFUbe2J7FsEILvSx6uB76Qu7MdRzFrVpWePI07QbcZsBdGDagXwuzaLveM9VqWr1AtEr5oTiGJyiNqQKwEwO5XIYK0pK2sVM82mp0DuaCJbNRaoT69/LlG4nm09/vGAD5gWaadUiD6/pBD9p2Y2F7jD9ee0V8lla+80s6bfw3XuPrXL37xM2tqTbi9GuXycyCp9hNhK6wNcjPum1KrgBCK4Wg6sandgEj2kxgB9cf9dCn7mEFuCW4A30vV/LYXm+//9fBY4n2Jrh6b0Hjfg3nRAiuVFymBR73gVMeLpTWByNs5Fer68is1BGSmRzhX1RcUao8hcDS9prXA08GmiLw7AbsyWA9/B3dOsvW20tjXovoXi2JrM94jAbEx6HEt0fV+aCum/TfyjshxYjvvWDn4/vwV+m7Ys/4xsVE+xOP16olYly1i9rjYLkoeUgbo9FaL/C/CAxYnVeR2cqScVKJeM1BlmzS/CNe7+sy5/8zIpIyU9QmSShhEZq1FAf9xnpHs/JgKWbRS3+l5LAU17XnsgkbHaiNlw4mTTESXPeBe4D/dFKdCqTqN4yv1eLiC2HEcT9VbNwFO4T0/3hwJ/UtZKupLFgZHvE41ynwameIx+kh5KGSmR3iSCitp0LXzf3ZlEgPXhdP05MdoiAimzqcQonFUl1+qWLFBpg0L3tcDX00nDUekaTbGQ/ePADZjVUZq1glfbnVVynwvMTacLs4xVbVprNdF6YVNJsfwbRLG8slZywLj7ejO7mpjF+Q9EIfqNHq1PpAfKJSN1llnzIvimZ71kKoAyUtLrRPaqGvcd6eOrKQZZ/dKznsB9eupVdRkxvuKddDV8YDYo6rnGEcHVKGLbEeAbuK9Jo3ZqiGBtJVHP1f4AKIKrzxHbjmOAD6U7M3F/NG0JnkG0kFiD+gDJAWBmY4AfEqdiZwH3uPqm9Ti5ZKRWxykk/WES6c4iCNpR8vEdAJgViO3BkcSpPYBDieHFA0peYQdm9+L+m1TnVEMEWO+meq55Lc+M+8dRrPEaAZxHsSM7mNUCV+LemOq5mtNr7d1KJYKuJ9JFGqFzPHHqEKLVwzXp450lpxznU8YsNJHWuPsbZvZVokXC3wNnmtn17v5OzkuTTpRLRsrMNgPfcvfZXf7mZVJGSmQfImt1KDGZPruexf01zE4hjohDtFSoJYKmB3Bflk4FDiE6qzeWvN6RRMA2EhhaUv90I9F7qoliFmwp7r9K96vaHJocrz2CYkYsu76N+1rMzgO+RHHLcTmwUu0XpFwW7TyuJ+r5fuTuT+S8JOkkuQRSR5ptr4fL3P3BLn/zMimQEtkPUVw+ggiuRgDDie28O3F/M7U6uJYIsuqJmYDvAnNxX4/ZAKC5JZCJjNRoYsske80G3K9L928nArO3KW41LsN9eRtrjPLNOME4jjjdOIpi1/fdwAVpq/NUYBCx7dh6VkzkfczsKOIQxqnAk0RA1bHTrNLt5LK1NxqqN8RfSiLSG0QAlGV2WvMKcDvFAOto4FPAw+n++cB4zBqIoc7vEY1IH8R9ZzoxWJoteoYItI4hitD7Ac8CtwFgNgXYRARrWdD2TsspQveFwMIUXB1BFOEfXVKndT5ZN3Zowuxtov1CljEbDGxR7ZWUcvdaM7uCGHx8OVCTtvr+kPPSZD/kEkgBDIrfPEVEsi7ne7dZyLJE8DIx3HlYuk4kitXnpPuXAl/AbB3FQGs17j9Jr3EcWV1m9NdyombqU7QcJGY+cH/Knk2kGGC9R9RPvVCystuIIG0kxQL8Q0ruTwFGYLaG2HpcQwRaOrnVy3lkL2eY2SIiOzXTzO4CZmv4cWXKLZAarEBKRPYlqz1wzzqnF0UdVPYPz0KitUIWaNUQmai5aatuPHA6ZhuJAvl6Ys7fL9PzTwXWptOHWcuGcUTwlrkP+DVmRxBH2+uIIOt5YAHQUPLcR4ng7RgimzUO+F/gxbT2acBmisFaLbAC99qyfk5Scdz9ZTO7iKibuoooRL/RY96mVJDc+kgdsudfUCIiHVNal+S+GFi8x/0IijJPEoHY4ekaDvRLXdLXYHYdkV3aDfyJCIoWAk8RgdZngF2pNuswYtDy0Pet6FZiO/CjRJ+sdUTg9CSwAVhfsq4VRJB1NsUarEeBmSkjdgvFgvys71dx61F6BHffbGaTiGD8H4GHzOwmd38u56VJB+SWkRoWJxdERA6MPQOtRcCiNp49gwiuDicCpaFAoWXsjtk/Ee0dMo1Ea4UFwFHAeOCkVBt1FJGB6s+emfdriQDtzPQea4jRNluJ3y+zrcNBxDbkx4BzKG49zgAWEAXLVxMZsXqKmbEVaZSPVJCS4ccvAz8CfmZms4Hpnp1glW4tl0BqOewcUOxwLCKSL/clwJI2njGRYoCVPb6VmofWEVmrUe/7mvlElulDwM3AhFQsfzgxm3A78QtlXyLoWpSyVZ8naq4aiNOGzURAls06HJHe/8NENisLtLKM2MnAFURGbB2RCVsPvID7hg79XKTLuPsKiy3oa4GvAaeZ2WR3X5Pz0mQfcgmk3ovftrbt84kiIt1BNOlsvVFnZA0uxKwvUV91KFF4Xp9aNzQCz6XPH5ue0wzMBh4jxttMA36eXrEvcBCwJT0eRRTGv0704Ps8EYhtS6+zPX38avr6sURgdywRoFUTGbQriS1G6abcfSdwaxp+/M/Ag2b2A3f/z5yXJm3IpY/UcLMdH4W//h/3x7v8zcukPlIi0mki81RIXdoPImqqhpRcg4HHcX8FsxOJIbhZl/cB6foVsa1XA/wFsU24hTjReCTRQ2s7MJD4pXkqsRNwIdHSoZaoB9tKbAs+ktZ2FpEp25KurcDmPWYpygFnxcHfpwC/Bqa6+/Z8VyWtySWQOtts99/AFROzqfMVQIGUiOQqMl4HE0HWwcT8wU2YjQA+S2SvBqXrYGK22yYig3VR+hgiwzWAKKTfTYzLGUkEYtuB04mgbRmRyRpCnDD8Zbr/d0SgtSx9vgpYnV5vB/Dp9LgxPX8HkZ17L30fA4Ed6rG1bxatOiYQ42VWApPc/c1cFyV7ya3YfEikr0VEpD1iC7GBPdssgPvbxDZh6+KXwFlEoDWw5PH11KV9DBE8DUifX0kET3OIIOkLwCeALxN/bx9LNDh9lQjqxhLbll9O73g6exbPjwC2YbYE2An8JbAUs+3pv3cC/477vDRselL63K6S+4twX5yCsHNL7jWm6y3c69KJx2PS55qAXVRwO4E04PguM3uR6E0226Jr/zzPIwsircolIzXEzKvhzfo4tVIpjk+PS3NdhYhITgysAFYFVbsiWKE/9OsP1QUoFKDQB6qqoKoBNlRB4VA4pBqqd0NzAQoD4aD1sD69VqEAhU2wuQE29oGqUTCyAIXsvkGhHupqYV1/6HdC8e/iFmthbR00DIT+x0dHewCaofkP8HpX/owOoD5E5nAIsSW7iqiRky7g7ud80L1cAqm+Zs0GSxor6+SeAikRkZz1gSqLYC4LtmwXNDVCUxUUBsMgA7PIivn6yvp3pj2GESdBG4l+ZJrV1wXaCqRy2dprgmfaWlR3pBopERHpDiwOINxKzKScDczSeJn8qLu4iIhIBXH314gDBL8h2lrcZTG6SHKgQEpERKTCuPtW4Aai2evJxHiZs/NdVe+kQEpERKQCeVhAjCiqB+4ws4kWrTKkiyiQEhERqWDu/hZwCfAIcDFwv0V/MWknM/uYmdVajFjqEAVSIiIiFc7dd7n7j4HvEqf65pjZF3NeViWZTAwbn9zRL1QgJSIi0kO4+9NEIfpSYIqZ3WRmA/JdVffn7he6+wp3v6ijX6tASkREpAdx91rgCuA+4K+Ijuij2/4qKZcCKRERkR7G3Zvd/W7gW0Q39Nlm9hUzs5yX1u2Y2blm9m/lfr0CKRERkR7K3RcRg6YXE+0SbrGYWShFNcBL5X6xAikREZEezN0bgKuAu4DziEL0MfmuqlupAV4ys2ozm2VmP+xI5k6BlIiISA/n7rvd/X6idqoamGVmf6utPiACqTrgKeC/3H2yd2AQsQIpERGRXsLdXyJO9b0AfA+4zcwOzndV+UnNS0cCDwGT3P2Bjr6GAikREZFexN03ANcA04DPEFt9J+W7qtycSASVTUBzOS+gQEpERKSXSVt9s4HLiVhgppld1Au3+mqAhURB/i/MbFhHX0CBlIiISC/l7kuIrb5ngWuBfzWzwfmuqkvVAK+4+1Jiq/ORjs4qtA7UU3UaM3va3c/p8jfeD2Z2L4C7T8h7LSIiIp0pZaIuILb81gOT3f3lfFdVGZSREhER6eU8PAx8nagXmmFml5qZ4oR90A9IREREAHD314CLgf8mek/dYWZD811V96ZASkRERFq4+xZgEnArcBrwoJmNzXdV3ZcCKREREdlD2uqbB1wKbAPuNrPLtdW3N/1AREREpFXpNNtXgf8guqLfaWaH57uq7kWBlIiIiHwgd98G3ATcDJxCbPWdle+qug8FUiIiItKmtNW3ABgPbAB+bmZXmllVzkvLnQIpERERaRd3XwF8DVhAtEq4x8yOzHdV+VIgJSIiIu3m7jvcfQpwIzAGeMjMxuW8rNwokBIREZEOc/cniUL0OmCamV3VG7f6FEiJiIhIWdx9FdEiYX56vLu3bfUpkBIREZGyuftOd/8BsdV3AnGq75M5L6vLKJASERGR/Za2+sYTQ4+nmdk3e0MDzx7/DYqIiEjXcPeVwCXAY8BlwPSe3sBTgZSIiIh0mnSq72aiiedJxFbfmTkv64BRICUiIiKdzt0fJ3pObSBGy0zoiVt9Pe4bEhERke4hNfC8BHgCmEAEVIflu6rOZe6e9xpEREREKpIyUiIiIiJlUiAlIiIiUiYFUiIiIiJlUiAlIiIiUiYFUiIiIiJlUiAlIiIiUqb/B3GvEWPPiOpyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAGKCAYAAADHWq+sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7zlY93/8dfHmZJzTil0S44jiZKkCUWSurvlTOXwuwsldeeU3Z6K0q2kIuSQcyoKoajIKYc7h3KWJMmhEBqMMZ/fH9dXJsbM3muvta91eD0fj/1Ymdlrrffs2Ps91/f6fq7ITCRJkjR6s9UOIEmS1KssUpIkSS2ySEmSJLXIIiVJktQii5QkSVKLqhSpiLi4xvuORUQcHRFH184hSZK6hytSkiRJLbJISZIktcgiJUmS1CKLlCRJUossUpIkSS2ySEmSJLVoREUqIlaLiPsjYtVOB+olfl0kSRpsI12R2g9Yt3nU8/y6SJI0wEZUpDJz68y8KzO36XSgXtKOr0tEvDsibouIOyNin3bmkyRJneUeqYoiYnbg28AmwMrA1hGxct1UkiRppEa6R2qpiPhRRFwXEbdGxNqdDtYL2vB1WRu4s1nVmgKcDryv/UklSVInzDGrT4iIOYDzgf0z89yImA+YvePJutysvi4RcSkw/wye+unMvKj530sDf57u9+4F1ulQZEmS1GazLFLAFsAtmXkuQGZO7myknjHTr0tmvm0ErxEz+LVsQzZJkjQORlKk1gB+0+kgPWimX5cRrkjdCywz3e+9CrivbQklSVJHjaRI3Q9MeO4fImKxzHyoc5F6xky/LiNckboGWCEilgP+AmwFeGekJEk9YiSbzU8AFo+ImyLieuAtnY3UM05gjF+XzJwK7A78DLgFOCMzb2prSkmS1DGzXJHKzCeAzcchS09p19clM88Dzht7IkmSNN6cIyVJktQii5QkSVKLLFKSJEktskhJkiS1yCIlSZLUIouUJElSiyxSkiRJLbJISZIktcgiJUmS1CKLlCRJUossUpIkSS2ySEmSJLXIIiVJktQii5QkSVKLLFKSJEktskhJkiS1yCIlSZLUIotUH4tgiwheVTuHJEn9ao7aAdQZESwLfJRSlu+tm0aSpP7kilT/2gyYAlxYO4gkSf3KItWHIpgfeAdwcSaP184jSVK/skj1p42AuYBzageRJKmfWaT60zzAtZncXTuIJEn9zM3mfSiTUyOI2jkkSep3rkj1mQheHUFkkrWzSJLU7yxSfSSCZYBvAxvXziJJ0iCwSPWXzYBngN/UDiJJ0iCwSPWJCF4GTAR+nck/aueRJGkQWKT6xzspd+udWzuIJEmDwiLVP9YHbs3kztpBJEkaFI4/6B/7AwvXDiFJ0iCxSPWBZtzB08Bfa2eRJGmQeGmvx0WwFHBkBCvWziJJ0qCxSPW+9wCLAw/WDiJJ0qCxSPWwCOYBNgSuyOSR2nkkSRo0FqneNhGYDzindhBJkgaRRapHNYcSbwbcCdxWOY4kSQPJu/Z628nAUx5QLElSHRapHtWUpytq55AkaZB5aa8HRbBYBFtFMH/tLJIkDTKLVG96N7AN5Ww9SZJUiUWqx0QwJ/Au4OpMHqqdR5KkQWaR6j1vARYAzqsdRJKkQWeR6j3voZypd13tIJIkDTrv2ushEcwFPEaZZO7IA0mSKrNI9ZBMpgBfqp1DkiQVXtrrERHMG8GStXNIkqTnWaR6x0Tg6AiWrh1EkiQVFqke0JyrtylwRyZ/qZ1HkiQVFqnesArwahx5IElSV7FI9YZNgSeAS2sHkSRJz7NIdbkI5gHeBFyUydO180iSpOc5/qDLZfJUBDvXziFJkl7MItXFIohMMpN/1M4iSZJezEt73W2dCA6JYJHaQSRJ0otZpLrbe4DFgEdrB5EkSS9mkepSzeDNNYDzM3m2dh5JkvRiFqnutQkwFbiwdhBJkjRjFqkuFMHcwIbAFZk8UjuPJEmaMe/a604JnAjcUTuIJEl6aRapLpTJFDwORpKkruelvS4TwTIRbNJMNJckSV3MItV9NgN2AeasHUSSJM2cRaqLNKtQ7wAuy+Tx2nkkSdLMWaS6y9uAeYHzaweRJEmzZpHqLpsA9wC31g4iSZJmzSLVJSKYD5gGXJBJ1s4jSZJmzfEHXSKTycCnI4jaWSRJ0si4ItUFIpg7gvkBXI2SJKl3WKS6wzuA70WwZO0gkiRp5CxS3eHdwF+A+2sHkSRJI2eRqiyC/wBei5vMJUnqORap+t4NPA1cXDmHJEkaJYtURRHMBbwd+HUm/6ydR5IkjY7jDyrKZEoEn6DMj5IkST3GIlVZJvfVziBJklrjpb1KIlghgv0iWKx2FkmS1BqLVD2bAG8A90ZJktSrLFIVNOfqrU/ZZD65dh5JktQai1Qd6wNzAz+rHUSSJLXOIlXHxsDdwB2Vc0iSpDHwrr1xFsHswLXAX51kLklSb7NIjbNMngVOrZ1DkiSNnZf2xlEEc0WwToQFVpKkfmCRGl9vBg4AVq0dRJIkjZ1FanxtDDwI3FA7iCRJGjuL1DiJYHFgAnChm8wlSeoPFqnxsxGQwEW1g0iSpPawSI2fVYDfZvK32kEkSVJ7ePfY+NkPeHntEJIkqX0sUuMggtkymQY8XjuLJElqHy/tdVgECwInRLBO7SySJKm9LFKdNxFYCPhL7SCSJKm9LFIdFEFQZkfdksm9tfNIkqT2skh11krA0sDPageRJEntZ5HqrI2BJ4HLaweRJEnt5117nXU+8LtMnqodRJIktZ9FqoMyuQ24rXYOSZLUGV7a65AI3h/BcrVzSJKkzrFIdUAESwEfAdasnUWSJHWORaoz3kk5oPhXtYNIkqTOsUi1WQSzUYZw/jaTh2vnkSRJnWORar/VgUWBX9QOIkmSOssi1X6LAg8CV9UOIkmSOsvxB22WyUUR/CKTrJ1FkiR1litSbRTB/BGEJUqSpMFgkWqvA4D9a4eQJEnjwyLVJhEsCawM3Fo7iyRJGh8WqfaZiLOjJEkaKBapNoggKEM4r8vk77XzSJKk8WGRao/VgMVwdpQkSQPF8QftcStwCM6OkiRpoFik2iCTKcCltXNIkqTx5aW9MYrgTRFsGcGctbNIkqTxZZEau/cCGwNTaweRJEnjyyI1BhEsBKwBXOI0c0mSBo9FamzWAwK4uHIOSZJUgUVqbDYA7srkz7WDSJKk8eddey1qNpf/A7ihdhZJklSHRapFmTwDTKqdQ5Ik1eOlvRZEEBEsUjuHJEmqyyLVmuWB4yNYp3YQSZJUj0WqNW8HngVurh1EkiTV4x6pUYoggLcBv83k8bpZIijjFwLIzJxWM48kSYNmTEUqIlYG1gZ+ATyamVWLxTh5LbAocHI7XzQi5gKWBpYBFgcWad5nEWAhYL7pPual/H8XL3iNZ4Cngaeme/wn8CjwyAs+/gbcD9yfmVPa+WeRJGlQjHVFak5gD+AVwDfHHqcnrAMkcE0rT25WkZYEVgJWBl4HvAZYgn+/1Jo8X3geAR4EJjcfTwLPANOaz5vWPHduYJ4XPL68eb+VKYVs9hlk+htwH6VY3QvcDfwRuDszn2zlzylJ0iAYa5FaBjgeuAuYH3hszIm63/mUIZwj+rM2xWl5SgFbG1gNWKD57anAH4DfAT8F/gzcA/yVssL3bDuDN1leTilUi1IK1lLN45LAKsCGTFfoIuJ+Sqm6E7gNuBW4x8uIkiSNsEhFxFKUFaflKZeVdsjMqzPz3IhYIjPv72TIbpLJw8CVM/uc5jLdW4F3UArUc6MS7gEuAW4CbgHuHM/LapmZwOPNxz0z+pyImJNSkJcFlms+lgfWoqxAAjwVEbdTStVNwI3Avc3rS5I0MGZZpCJiDsoqzP5NcZqP6S4PDVKJgg0XjmAT4MJMpk7/O00BWQfYmHJ0zHyUyedXNR9XZ+Zfxzfv6GXmM5QVxrum//Xm34PlgNcDKzaPmwFbNp/ycETcSClVNwC3uPdKktTvRrIitQXlh+K5AJk5ubORutlmSwPbAhc89yvNat0HgPdRLpk9DlwE/By4pt2X52rJzKnAHc3HOQARMRulXE0AVm8eN2ieMjUibgauBa4GbrRYSZL6zUiK1BrAbzodpDf8xwLArZlkRKwO7EQZhZDApcBPgCubVZ2+1+yT+kPzcSZARCxM2Qc2AXgD5Wv0EWBKRFxPKVXXUMq5+6wkST1tJEXqfsoPRQAiYrHMfKhzkbpVAIvMCxfOGbHxMZSS8A/gOODMzHygarwukZkPU/aBXQIQES8D1gTeRNlsv3vzqY9HxDXAZcBlzfMkSeopIylSJwCnRsRNlFvuDwTO7mSo7rTIovDEonDCNpTLW4cCP3Y8wMxl5j8pq3WXwr9WrNai7CdbF5jY/PrNNKUKuNXVKklSL4gaN1rttttujy611FKfHBoaOmF4eHhO4ELgu0NDQycPDw/PB5wHHDk0NPT94eHhBSiXzA4fGho6c3h4eFHgh8ChQ0ND5wwPDy8BnA58eWho6ILh4eFlgJOALw4NDV00PDy8PGXVaGhoaOiS4eHhFYGjgP2GhoauGB4eXhX4FvCZoaGha4aHh9cADgM+OTQ0dP3OO++84XzzvezYM8748fwPPDAx4YkvwA+PcL/P2DXjGFYA1ms+VqMs/T1MKVQXA7/xay1J6lYeEfMSnnnmmYiILZdddo39Jk7cdsmIXR6Ez/08Mw+rna1fNOMSbm8+jouIBSmrVOtRVqo2ByZHxGXAr4DLB/tmB0lSt6myIhURF2fmBuP+xiMUEa8CJgGrw2cfhv1ng1tXgD2uzPzNVrXzDYJmnMRalFlcGwALA1MoNz78CrgkMwdhAKwkqYu5IjWd5lLT+4C9gWdh9gPh4C0gnoDvJ9zpasg4ae58vBK4MiK+TLnhYWLzsT7wTERcAfwM+HVmPlUtrCRpYLki1WgGjR5IOSLlWmAI8lnKpvKTILYEyMxd66VUU3ZXBt4FbAQsRjl78GLKfK+rmplXkiR1nEUKiIhXA/9LORbl28BJz901FsFswGwQR4BFqps0A0HXpJSqd1IOz/4H8AtKqbreu/8kSZ008Jf2IuItwMHAs8DumXl1+fVycG8m04BpEfUyasaaknQtcG1EHAK8hVKqNqVMm78vIs4Fzs3M++ollST1q9lqB6gpIt4DfAO4D9j+uRLVeANwfASvrhJOo5KZz2TmrzNzf8olv88B9wK7AGdHxFER8Z6ImLdqUElSXxnYFamI2B74BOW4kk83gyOn93ZgbqDrDxrWv2uGpJ4PnB8RSwDvAd4LDAOfjYiLKOcFXp81rm1LkvrGwBWpZrPyf1POf7sIOPCFAx8jmItymejSTAbi3Lx+lZn3A8dGxHGUcyPfS7mhYHPg3og4CzjHI2okSa0YxEt7O1NK1I+B/V5iavbawDw058Wp92VxXWZOAt4NDAEPAnsA50XEwRHxpqZoS5I0IgO1IhUROwC7AT8FDprJHV3rA48AvxuvbBo/zXT0nwI/jYjlgPcDm1H2Vt0z3SrVoxVjSpJ6wMCMP4iIDwD7Uc71OyAzn33pz2UCsEAmv57u+UeD4w/6VUTMTRmh8J+U4Z/PAL8EfgRc514qSdKMDMSKVES8GdgHuBz43MxKFEAmN4xLMHWNzHyaclj2eRHxWsr4hPdQxincGRHfB853grokaXp9v0equXTzFeAuYN9ZTb2OYKMIXjMu4dSVMvMPmflVyl6qLwDTgP0pJWvPiFiqakBJUtfo6yIVEQsChwFPA3s1e2Nm8vm8AtidckiuBlxmPpWZPwG2pcyjuhrYDvhxRBzq5nRJUt9e2muOD/ki5Sy23TJzJPOg1qWUy1/P6hM1OJr9UdcB10XE4pR9VB+gzBq7KyJOp1z2e7JiTElSBf28IrUd8Gbg0Mwc6d1361OmYd/dqVDqbZn5QGYeQTmGZhiYQrmJ4bnLfq+sGlCSNK76skhFxCrAxyl3XZ05suewMLAqZQind2hppjJzSmaeA2xPmUv2G0p5PyciJkXE66oGlCSNi767tBcRLwMOAh4CvjiK29aXo6wuXNqpbOo/zb9fNwI3NpvQtwK2ADaNiKuBk4ErHZ8gSf2pH1ek9gKWBPbPzMdG+qRM/g/YJpM/dyyZ+lpm3peZX6OMTTicUs4PB06PiM0jYq6qASVJbddXRSoi1qSsBpycmSOeBRVBAGQyo+NipFHJzMcz80TKeX5DlPEJB1Iu+30kIl5RNaAkqW36pkg1f9vfH7gPOGaUT988gsMjmLf9yTSoMvOZzPwpsA1lz94dwMcoR9Ps5cZ0Sep9/bRHaifgNcCeLdyG/jYgM/H2dbVdsz/qKuCqiFiBskF9a2DLiDgXODEzvaQsST2oL1akIuI1wIeBn2fmFaN7Lq8EVgQu60Q2aXqZeUdmHki5BP0Tyn6qH0XEQd7pJ0m9py+KFPAJyh13/9vCc9/aPFqkNG6ajelfBt5LubNvPeDUiPhGRKxRN50kaaR6vkg1P3TWB76XmQ+38BLrAX/IZCSTz6W2ysy/Z+bhlJWpI4BVgO9GxDERsa5H0EhSd+vpItX8kNkD+Dtw2uifTwDnAT9oczRpVJo7/Y4DNqOsrC5FGZ1wYkSsb6GSpO7U00WKspo0ATi6lXPOMslMfpHJ5e2PJo1ec1Dy6ZQ9VJOA+YGvASdHxAYWKknqLj1bpJpDiXcH7qFs2m3hNVg/gkXaGkxqg2Z0wtmUA5I/D8xHWak6JSImNv/+S5Iq6+Vvxu8CXgscmZlTR/vkCBYFPgO8s93BpHbJzGcz81zgg5ShnvMAh1A2pm9ooZKkunrym3Dzw+PDwJ3AL1p8mefu1vOynrpeU6jOoxSqAygz4L5MOX5mYwuVJNXRq9981wOWB07IzGktvsZbgbsz+Uv7YkmdlZnTMvMCYEtgv+aXDwK+HxHvtlBJ0vjq1W+6O1GOgrmwlSdHsDCwEq5GqUc1hernwFbAvkACX6SsULmHSpLGSc99s42IVYHVgVMy89kWX2Ylyg8ei5R6WlOoLuT5QjUbZQ/VSRHxVu/yk6TO6rkiBXwImAyc0+oLNOMOts/E883UF6YrVFsCQ5SxCd8Ajo2IN1UNJ0l9rKeKVEQsAmwEnJ2Zk8fyWpn8oz2ppO7RFKqfAh+g7J1aAjgyIo6MiNXrppOk/tNTRQp4H+VupTNafYEIJkYwKYKXtS+W1F0yc2pmnkkZ7HkoZVTIcc1Zfq+vm06S+kfPFKlm8+z7gGsz854xvNT6lOM3xrSiJfWCzJySmadR/tv5FrAaZUr6IRGxfN10ktT7eqZIAWsCS9PiFHOAZhVqAnBFJtmuYFK3y8wnM/MEYHPgaODNlJEJn4+IJaqGk6Qe1ktFagvgceCXY3iNN1IuDV7ZlkRSj8nMJzLzaOC9wCnAxsCZEbFXRCxQN50k9Z6eKFIRMR/wDuBnmfn0GF5qXeAR4Na2BJN6VGb+IzMPo2xKvwDYGjg7Ij4SEfPWTSdJvaMnihSwATA3cP4YX+cm4Cwv60lFZt6fmZMoc6j+D/gY8OOI+EBEzFE3nSR1v14pUu8C7gd+N5YXyeScTM5qTySpf2TmXZn5KeCjwL2U42fOiIh3OtRTkl5a1xepZt/Gm4Gfj+FcPSJYIYJ52pdM6j+ZeQOwM/ApYCrwFeB7EbFW1WCS1KW6vkhR9kbNDvys1ReIYHZgmHLZQtJMZPFryuW+YWAR4DsR8c2IeF3ddJLUXXqhSG1AOaD49jG8xkqUIzOuakcgaRA0U9LPAd4PHAasApwaEZMcmSBJRVcXqeZuvbWBizNzLBvE3ww8A/y2LcGkAdIM9TyZMtTzBGBDysiEj0eEJwRIGmhdXaQoBWgu4JJWXyCCaF7nhkyebFcwadBk5uOZ+S3KyISLgA9T7vD7YETMXjedJNXR7UXq7cBjwPVjeI1XAYsDv2lLImnANSMTDgR2AP4I7EOZkv427/CTNGi6tkg1Z+u9FbgsM58dw0vdC+wBXN6WYJIAyMybgd2AvSnfS74OHBkRK1YNJknjqGuLFPA6YEHGuJKUSWZydyZPtCeWpOc0d/hdAmwJHAKsQDkUeTgiFq+bTpI6r5uL1Jubx5bvtItg0Qg+GcHSbcokaQYyc2pmnkE5E/Mkyhl+Z0XEx5qbRiSpL3VzkVoHuCMz/z6G11gbeCfd/eeU+kazIf1w4D+BXwEf4fkjZ/zvUFLf6cpvbBExD7AGY5/7tDbwV8o+KUnjJDPvy8z9gR2BeyhHzpzihHRJ/aYrixSwJjAnY7usNw8wAbjaQ4qlOjLzJmAXyp19L6dMSP9qRHi5XVJf6NYi9QZgGmMbe/AGYA6cZi5V1WxIvwj4IHAEZf/jDyNiD/dPSep13VqkJgC3ZuZYBmjOBdwF3NKeSJLGIjOfzszjKAM9f0657PfjiNjc/VOSelXXffOKiLmAVRnbahSZXJLJJzKZ2p5kktohMx/KzCFKkboXOBD4XkSsUTeZJI1e1xUp4PWU1aSWi1QE8zRHw0jqUs3+qY8CBwCLAN+NiIM8EFlSL+nGIvXc30rHsiL1X8AJEczVhjySOqTZP3UBZVzC0ZRjoc6MiP8XEfPWTSdJs9atRerPmfnwGF5jLeCvmUxpUyZJHZSZT2bm0ZT9U78CdgZ+FBHv8vw+Sd2sq4pU8w1zVeDG1l+DRYHlgWvblUvS+MjMB5r5Ux8B/g58CTgqIlaom0ySZqyrihSwKLAwY7vT7rmBf1ePPY6kGjLzRspm9C8Br6UM8/xMRLyibjJJ+nfdVqRe3zzeOobXeBPwIPDnsceRVEtmTsvMs4D3Az+k7H08MyK2cFyCpG7Rbd+MVgISuH0Mr/ET4DinmUv9ITMfy8xDgO2Auyl3+Z0QEatUDSZJdF+Rej3wp8yc3OoLZHJjJpe3MZOkLpCZt1OOmzkAeCVl9tSBEbFw3WSSBlk3FqmWL+tFsFYEK7Yxj6QuMt24hA8AJwKbUi73bRURs9dNJ2kQdU2Rav5W+UrGttF8Z2Dr9iSS1K0yc3JmHg5sBfwe+DRwakS8sW4ySYOma4oU/Gsl6bZWnhzBksDSwP+1LZGkrpaZdwN7UIrUvJRRCQdFxCurBpM0MLqpSC3fPN7Z4vOfG3vg/ChpgDSX+y6m3NV3NLAB5XLfDhExZ81skvpfNxWp5YBHM/PRFp//RuC+TP7axkySekRmPt1MR/8gcBWwJ2X+1Jp1k0nqZ91UpJYF/tjKEyOYHVgBL+tJAy8z78vMvYG9gHmAoyNiUkQsUjmapD7UTUVqOVosUpk8C+wEnNbOQJJ6V2ZeCmwJHAtsRDm7b0uHeUpqp674hhIRCwELUIbttSSTZzJ5vG2hJPW8zHwqM48EPgTcDPwPcKLDPCW1S1cUKcpqFLRYpCL4bAQbtS+OpH6SmfcAHwf2BRahTEbfz7P7JI1VtxSpZZvHUV/ai2BRYD3g5e0MJKm/NHf3XUjZjH4asAXl7r7NvdwnqVXd8s1jWeAp4IEWnvuG5vG6tqWR1Lcy85+Z+TVgW+Ae4EDgmIhYoW4ySb2oW4rUcsDdmTmtheeuCTwM/Km9kST1s8y8g3IawiTgNZRRCZ+KiPnqJpPUS7qlSC1D+ZvhqEQwG7AGcF0m2fZUkvpaZk7LzLMpZ/f9mHLE1I8iYmJERN10knpB9SLV7E1YHFoapPky4EbgyraGkjRQMvOxzDwI+DDwCHAI8LWIWLJuMkndrnqRotxBMyctFKlMHs/k4Eyuan8sSYMmM38PbA98HXgT8IOI2C4iZq+bTFK36oYi9dzf+EZdpCLw1mVJbZWZz2bmKZSz+64GPgmcFBGr1k0mqRt1Q5Faonm8fzRPimAe4HsR/Ff7I0kadJn5V2Bv4DPAgsDxEfE/EeGoFUn/0g1FqtUVqZWBOYA72xtHkopm9tSvKKtT328efxgRG7oZXRJ0T5F6LDMnj/J5awBTKcc+SFLHNLOn/hfYEfg78GXgsIhYqm4ySbV1S5Ea1WW9xhrAzZk83eY8kjRDmXkzsAPwNcoMux9ExA4RMUfdZJJq6ZYiNarLehEsQBnieX1HEknSS2g2o59KOWrmN8CelM3oq9VNJqmGqkWq2WOwBKPfHzWFcnvy5W0PJUkjkJkPZObewKeBBYDjImKfiJi/cjRJ46j2itT8wHyMskhl8mQmv8zkvs7EkqSRycyLef4g5A/gZnRpoNQuUi3dsRfBOyN4ZQfySNKoZebk5iDkHYGHKJvRD40Iv09Jfa52kVq8eXxgpE+IYDHKgLy3dCSRJLUoM2+hlKlvAOtQVqc+2ByFJakP1f6Pe8Hm8ZFRPOe5DZ03tjmLJI1Zsxn9JOBDwO+BfYBjImK5uskkdUKvFqnHgbvbnkaS2iQz7wU+DgxT7jI+NSJ2jog56yaT1E61i9RClDvwnhrFc1YHfp9JdiaSJLVHMxn9HMpm9F8B/w84JSJWr5tMUrvULlILAo9k5ohKUQQLA4sBv+toKklqo8x8ODP3A/YCXg4cGxGfiYj5KkeTNEa1i9RCwKMj/eRMHga2BX7ZsUSS1CGZeSlldeoHwJaUyejr1U0laSy6oUiNZn8UmTyeyT87lEeSOqoZlXAI8FHgn5Qz+w6KiIUrR5PUgtpFakFGuCIVQUSwdwTrdDiTJHVcZt4IbAccBUykjErYzEGeUm+pXaRGsyK1OLABsGjH0kjSOMrMKZl5DLAN5U7kzwPfjoila+aSNHLVilREzEU5Hmake6RWaR5/35lEklRHZt4F7Ax8BVgVOCMito+I2esmkzQrNVekRjtDamXKfoJ7OhNHkurJzGmZ+QPgv4CrgE8Ax0fEa+smkzQz3VCkRroitTJws/OjJPWzzHwA2BvYD1iKMndqFwd5St2pJ4pUBHMA/wBu6GgiSeoCzSDPn1NWp34J7AacGBEr1U0m6YVqFqmFmsdZXtrLZGom+2Tykw5nkqSukZmPNIM896Z8z/xeRHy82WMqqQv0yoqUtwNLGliZeQllgOe5wIcp5/Z5zIzUBWqvSCXlkt2sfDmCj3U4jyR1rcx8LDMnAXsA81COmflURMxbOZo00GqvSD2WmdNm9kkRzAOsCDwxLqkkqYtl5pWU1akfUuZPnR4Ra9VNJQ2u2itSIxl98HpOnAIAAA49SURBVDpgduDmzsaRpN7QHDPzFWBXYBrwnYjYLyJeVjmaNHBqr0iNZPTBypRLgLd0No4k9ZbM/C2wNXAysAXlEOR166aSBkvNIjUvMHkEn7cK8CcPKpakF8vMpzLzMJ4/BPnwiPh8RLyicjRpINQsUrMDz47g864EftrhLJLU0zLzd8C2wLHAppRDkN9RN5XU/2oWqTmAqbP6pEzOy+SCccgjST2tOQT5SGB74G/AVyPi4IhYaBZPldSiXliRkiSNQmbeBuwAHAFsQDkEeWLVUFKfql2kZrkiJUkavcycmpnHAdsB9wOHRMRBEbHgLJ4qaRRqX9pzRUqSOigz/0CZhn4kMJGyOuXeKalNXJGSpD7XrE4dS9k79RBl79SXImKBytGknueKlCQNiMy8A9gROArYkDJ36u11U0m9rXaRckVKksZRszp1DM/f2XdoRExy7pTUGouUJA2gzLydsjp1NPAuyurU+nVTSb2n9h4pL+1JUiWZ+UxmHk0ZlfAw8LWIGHZ1Sho5V6QkacBNN3fqu8AmwPcjYr26qaTe4IqUJOm51anvUC73PQYcFhFDETF/5WhSV6tZpMAVKUnqKpl5C2Uj+nHAeyirU+vWTSV1r1pFKppHV6Qkqcs0Z/YdAewEPAEcHhEHRsTL6yaTuk/tIuWKlCR1qcy8mXLEzPHAZpSp6G+um0rqLrWLlCtSktTFmtWpb1OOmZkMfCsi9omIeStHk7pC7SLlipQk9YDMvAnYFjgF+E/gtIiYUDeVVJ9FSpI0Ipn5dGZ+HdiN8vPjmIjYMyLmqhxNqsYiJUkalcz8LbAV8BPK/KmTImLFuqmkOmoXKfdISVIPyszJmfkl4BPAAsCJEfHRiJi9cjRpXDlHSpLUssy8HNgS+AXw38DxEbFszUzSeHJFSpI0Jpn5WGbuB+wLvAo4NSK2jojaf1mXOs4iJUlqi8y8kLI6dTWwN3BkRCxVN5XUWbWLlJf2JKmPZObfgL2AScBKwOkR8b6IiJk/U+pNtYuUK1KS1GeyOJtyZ98twOeAr0fEonWTSe1Xu0i5IiVJfSoz76NsQD8UWJtyxMxGdVNJ7WWRkiR1TGZOy8zTgG2Ae4GDI+KgiFigcjSpLSxSkqSOy8y7Kef1HQm8k7I69daqoaQ2qF2k3CMlSQMiM5/NzGMp09AfBb4REft6ALJ6We0ZH65ISdKAyczbgO2Bk4EPAKdExKp1U0mtcUVKkjTuMnNKZh5G2Yw+N3BcROwaEXNUjiaNSu0i5YqUJA2wzLyWMibhAmBXSqF6dd1U0sjVLlKuSEnSgMvMxzPzQJ4/Yua0iPigQzzVC2oXKVekJEnAv46Y2Qq4DtiHshndIZ7qarWLlCtSkqR/ycwHgT2AQ4C1gO9HxMS6qaSXVrtIuSIlSfo3zREzZ1CGeN4HHBIRn4+Il1eOJr2IRUqS1JWmG+L5XWBTyt6pNauGkl6gdpHy0p4k6SVl5tTM/A6wM+VnxlERsWdEzFU5mgQ4kFOS1AMy80bKpb6zKJPRvxcR/1E3leSKlCSpR2Tm5Mw8CNgLWAQ4KSK2i4jaiwIaYBYpSVJPycxLgQ8BVwCfBI6MiCXqptKgqlmkpmVmVnp/SVIPy8xHgE8Dk4CVKGMSNnWIp8ZbzSLl/ihJUsuaMQlnA1sDd1JK1UER8Yq6yTRILFKSpJ6WmX8BdgGOACYCp0fEG+um0qCwSEmSel5mTsvM4yhzp54CvhMRe0TEnJWjqc/VLFJuNJcktVVm3gxsC/wY2BE4PiKWrZlJ/c0VKUlSX8nMJzPzS5TN6EsCp0TEf7oRXZ1Qc/aGK1KSpI7JzIuBrYDrgH2BQyNioaqh1He8tCdJ6luZ+RCwJ/A14C2UMQnrjmeGiFgtIu6PiFXH8301Pry0J0nqa81G9FMpR8s8AhweEZ+OiLnHKcJ+wLrNo/qMK1KSpIGQmXcA2wOnUS75nRQRK4zD+26dmXdl5jadfi+NP1ekJEkDIzOnZOahwB7AAsCJEbGN5/WpVRYpSdLAycwrKatSVwCfAr4ZEYu1+30iYsOIOKkNr/PuiLgtIu6MiH3akU3tYZGSJA2k6c7rOwiYQJmI/o42v80Eyl2DLYuI2YFvA5sAKwNbR8TKbcimNnCPlCRpYDXn9Z1JGeJ5H/DViDggIuZr01tMAK6LiLkj4oSIOKiFeVZrA3c2+6ymAKcD72tTPo2RK1KSpIGXmX8CPgIcTykpJ7dp1WcC8CDwM+CizNwvMxMgIi6NiOtn8LHhC15jaeDP0/3zvc2vqQvMUel9XZGSJHWVzHwG+HZEXAl8gXK8zFHACZk5bbSv15zztyzlLsHdmn1Z07/f20b6UjOKO9o86oyadym4IiVJ6jqZ+Vtga+AXwMeAoyJiyRZeamXgGsrPuxctHoxiRepeYJnp/vlVlMuQ6gKuSEmS9AKZ+VhE7A9cBuwDnBYRB2Xmz0fxMhModwWeDJwVERMz84Hp3mOkK1LXACtExHLAXyh3GzqTqku4R0qSpBloNqKfRykufwQOiogDR7ERfQLw+8y8HfgscEZzuW+0OaYCu1P2Wd0CnJGZN432ddQZ0ex5G983jXgIOCYze2ZcfkQcDZCZu9bOIkkaXxExB7ALZUP6n4H9M/OWuqnUDVyRkiRpFjJzamYeCewGzEPZiL69E9FlkZIkaYSajehbAb8GPkGZiL5o3VSqyYGckiSNQmY+Rtnz9CVgDcpE9JFuHFefcUVKkqRRajainwVsRxm4+fWI+ExEzFU5msaZK1KSJLUoM/8I7AScCnwIODEilq8aSuPKgZySJI1BZk7JzK9R9kwtQjle5oMtnKmnHuSKlCRJbZCZl1M2ov+WMsTzqxGxQN1U6jT3SEmS1CaZ+XdgT+DrwHqUjehr1U2lTnJFSpKkNsrMaZl5CrAjMBk4MiI+3gz1VJ9xj5QkSR2QmbdR7uo7G/gwcGxEvKpuKrWbRUqSpA7JzCcz8wuUPVOvBk6NiE0rx1IbWaQkSeqwzLwI2Bq4HZgUEV+IiJdVjqU2qFmk3CMlSRoYmXk/5ay+o4B3AadFxGp1U2msXJGSJGmcZOazmXkMsAvlxqtjI2InDz/uXa5ISZI0zjLzBmAb4JfA7sC3ImKRuqnUCouUJEkVZObjwL48f/jxaRHxlrqpNFpe2pMkqZIXHH78KPDNiNgzIuasHE0j5IqUJEmVZeZdwPbAj4AdgO9GxNJ1U2kkXJGSJKkLZObTmXkw8D/AaygzpzauHEuzYJGSJKmLZOYvKRvR/wAcFBGfi4h5K8fSS7BISZLUZTLzPmBX4Dhgc+DEiPiPuqk0I+6RkiSpC2Xm1Mw8Avg48ApKmfpgRETlaJpOZOb4v2nEVOB64Ilxf/PWva55vL1qCknSIJoDWBZYAHgE+BMuSIybzNzgpX6v5orU+Dc4SZJ601TgTuBeYEFgZcCz+rrAHJXe97KZtbtuFBFHA2TmrrWzSJIGV0SsDBwMLAmcCJyQmdPqphpcnu0jSVIPycybKXf1XQh8DDgiIharm2pwWaQkSeoxmflP4ABgErAq5XiZt9ZNNZgsUpIk9aDmeJmzKRPRHwK+ERF7ebzM+LJISZLUwzLzj8COwBnAtsBxEbFM3VS9JSJWi4j7I2LV0T7XIiVJUo/LzCmZeQjwaWBp4JSI2KRyrF6yH7Bu8zgqFilJkvpEZl5M2Yh+O/CFiBjyeJlZy8ytM/OuzNxmtM+1SHVARKwcETtFxDIRMX/tPJKkwZGZ9wO7Ad8FNsPjZTrKIjUGEbFhRJw0g9+aE9gDeD+9Nb1dktQHMvPZzPwO5XiZBShlaguPl3mxmfwsHxGL1NhMAK6bwa8vAxxPmULripQkqYrMvBrYinIs2wHAFyNivrqpus5L/SwfEYvU2EwArouIuSPihIg4KCIiM88FfpiZ52XmY7VDSpIGV2Y+DOwOHAFsRNmIvmLdVF1lhj/LR/pki9TYTAAeBH4GXJSZ+2VzCnRzjVqSpOoyc1pmHkfZOzU3cEJE/JeX+oCZ/CwfCYtUi5qBZ8sCpwH7ZubJdRNJkjRzmXkd5a6+a4DPAl8Z5Jui2vGz3CLVupUp/yJOBZ6tnEWSpBHJzEeBTwKHA2+nXOpbpW6qasb8s9wi1boJwBWUTXzHR8TilfNIkjQizaW+E4FdKF3g2IjYZgAv9Y35Z7lFqnUTgN9n5u2U5dEzPN9IktRLMvNGyqW+y4BPAYdGxCvqphpXY/5ZHqPYT9U2EXFxZm4w7m88BhFxNEBm7lo7iyRJ7dSsRH2Icsnv78B+mXlD3VS9wRUpSZIGXBanAx+h7Bc6pjmhw54wC36BJEkSAJl5M7At8EvK7KlvRMTCdVN1N4uUJEn6l8x8AtgXOBh4I3BqRKxZN1X3skhJkqR/01zq+xGwEzAZ+E5E7OKlvhfzCyJJkmaouZttO+ACylT0b0fEonVTdReLlCRJekmZORkYAiYBq1Mu9a1TN1X3sEhJkqSZai71nQ1sDzwKfCsiPhYRs1eOVp1FSpIkjUhm3gXsAJxNGZVwVES8sm6quixSkiRpxDLzqcz8AvA5YEXgtIhYt3KsaixSkiRp1DLzfMpG9AeBwyNi90G81GeRkiRJLcnMP1FGJJzVPH5n0C71WaQkSVLLMvPpzPwS5VLf6yl39b2lcqxxY5GSJElj1lzq255y6PHhEfHfgzDAs+//gJIkaXxk5t3AjsA5wEeBI/t9gKdFSpIktU1zV98kyhDPVSiX+tauHKtjLFKSJKntMvOnlJlTj1KOltm1Hy/19d0fSJIkdYdmgOeOwHnArpRCtUjdVO0VmVk7gyRJUk9yRUqSJKlFFilJkqQWWaQkSZJaZJGSJElqkUVKkiSpRRYpSZKkFv1/AI8sIBqzpSsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"1|IMPORT PACKAGES\"\n",
    "import numpy as np              # Package for scientific computing with Python\n",
    "import matplotlib.pyplot as plt # Matplotlib is a 2D plotting library\n",
    "\n",
    "\"2|DEFINE PARAMETERS AND ARRAYS\"\n",
    "# PARAMETERS\n",
    "k_size  =  150        # Model domain\n",
    "steps   =  k_size*100 # Number of \"dots\" in the domain\n",
    "alpha   =  0.30       # Output elasticity of capital\n",
    "delta   =  0.35       # Depreciation rate\n",
    "rho     =  0.35       # Time preference\n",
    "n       =  0.05       # Population growth rate\n",
    "g       =  0.05       # TFP growth rate\n",
    "theta   =  0.8        # Coefficient or Relative Risk Aversion\n",
    "# ARRAY\n",
    "k_array = np.linspace(0, k_size, steps) # Create array of k\n",
    "\n",
    "\n",
    "\"3|DEFINE FUNCTIONS\"\n",
    "# PRODUCTION FUNCTION\n",
    "def y(k):\n",
    "    y = k**alpha\n",
    "    return y\n",
    "\n",
    "# MARGINAL PRODUCT OF CAPITAL\n",
    "def y_prime(k):\n",
    "    y_prime = alpha*k**(alpha-1)\n",
    "    return y_prime\n",
    "\n",
    "# CONSUMPTION\n",
    "def consumption(k):\n",
    "    c = y(k) - (n + g + delta)*k\n",
    "    return c\n",
    "\n",
    "# MOTION FUNCTION OF CONSUMPTION\n",
    "def cdot(k, c):\n",
    "    cdot = 1/theta * (y_prime(k) - delta - rho) * c\n",
    "    return cdot\n",
    "\n",
    "# MOTION FUNCTION OF CAPITAL\n",
    "def kdot(k, c):\n",
    "    kdot = y(k) - c - (n + g + delta)*k \n",
    "    return kdot\n",
    "    \n",
    "\n",
    "\"4|CALCULATE STEADY STATE AND GOLD VALUES\"\n",
    "# STEADY STATE VALUES\n",
    "k_star = (delta + rho)**(1/alpha)\n",
    "c_star = consumption(k_star)\n",
    "y_star = y(k_star)\n",
    "s_star = (y_star - c_star)/y_star\n",
    "\n",
    "# GOLD VALUES\n",
    "k_gold = (alpha/(n+g+delta))**(1/(1-alpha))\n",
    "c_gold = consumption(k_gold)\n",
    "y_gold = y(k_gold)\n",
    "s_gold = (y_gold - c_gold)/y_gold\n",
    "\n",
    "\n",
    "\"5|CALCULATE LOCI FUNCTIONS\"\n",
    "k_loci = k_star\n",
    "c_loci = consumption(k_array)\n",
    "\n",
    "\n",
    "\"6|CALCULATE SAMPLE PATHS\"\n",
    "k00, k01, k02 = k_star*0.25, k_star*0.50, k_star*2.2\n",
    "c0 = [consumption(k00)*0.25, consumption(k00), c_star*1.25,\n",
    "      consumption(k01)*0.40, consumption(k01)*1.50]\n",
    "\n",
    "def sample_path(k0, c0, n):\n",
    "    path = np.zeros(shape=(n, 2))\n",
    "    path[0, 0] = k0\n",
    "    path[0, 1] = c0\n",
    "    \n",
    "    for j in range(n-1):\n",
    "        if path[j,0] < 0: # Stop if motion goes off-chart\n",
    "            break\n",
    "        else:\n",
    "            path[j+1, 0] = path[j, 0] + kdot(path[j, 0], path[j, 1])*0.25\n",
    "            path[j+1, 1] = path[j, 1] + cdot(path[j, 0], path[j, 1])*0.25\n",
    "    return path\n",
    "\n",
    "path1 = sample_path(k00, c0[0], 30)\n",
    "path2 = sample_path(k00, c0[1], 30)\n",
    "path3 = sample_path(k01, c0[2], 30)\n",
    "path4 = sample_path(k02, c0[3], 30)\n",
    "path5 = sample_path(k02, c0[4], 30)\n",
    "\n",
    "\n",
    "\"7|PLOT RAMSEY MODEL WITH SAMPLE PATHS\"\n",
    "# Value of k such that c = 0\n",
    "k_zero = (1/(n + g + delta))**(1/(1-alpha))\n",
    "# Axis range\n",
    "y_max = np.max(c_loci)*1.7\n",
    "x_max = k_zero*1.2\n",
    "v = [0, x_max, -0.1, y_max]                        \n",
    "\n",
    "### BUILD PLOT AND POPULIATE WITH LOCI LINES\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "ax.plot(k_array, c_loci, color=\"K\", alpha = 0.8)\n",
    "ax.axvline(k_loci, color=\"K\", alpha = 0.8)\n",
    "### SAMPLE PATHS\n",
    "ax.plot(path1[:,0], path1[:,1], \"r--\", alpha = 0.7)\n",
    "ax.plot(path2[:,0], path2[:,1], \"r--\", alpha = 0.7)\n",
    "ax.plot(path3[:,0], path3[:,1], \"r--\", alpha = 0.7)\n",
    "ax.plot(path4[:,0], path4[:,1], \"r--\", alpha = 0.7)\n",
    "ax.plot(path5[:,0], path5[:,1], \"r--\", alpha = 0.7)\n",
    "### ADD DOTS AT BEGINNING OF SAMPLE PATHS\n",
    "plt.plot(k00, c0[0], 'ro', alpha = 0.7)\n",
    "plt.plot(k00, c0[1], 'ro', alpha = 0.7)\n",
    "plt.plot(k01, c0[2], 'ro', alpha = 0.7)\n",
    "plt.plot(k02, c0[3], 'ro', alpha = 0.7)\n",
    "plt.plot(k02, c0[4], 'ro', alpha = 0.7)\n",
    "### ADD MOTION ARROWS\n",
    "plt.arrow(k00, c0[0], 0, 0.10, head_width=0.03,head_length=0.02,fc='r',ec='r')\n",
    "plt.arrow(k00, c0[0], 0.18, 0, head_width=0.02,head_length=0.03,fc='r',ec='r')\n",
    "plt.arrow(k01, c0[2], 0, 0.10, head_width=0.03,head_length=0.02,fc='r',ec='r')\n",
    "plt.arrow(k01, c0[2],-0.08, 0, head_width=0.02,head_length=0.03,fc='r',ec='r')\n",
    "plt.arrow(k02, c0[3], 0,-0.08, head_width=0.03,head_length=0.02,fc='r',ec='r')\n",
    "plt.arrow(k02, c0[3], 0.22, 0, head_width=0.02,head_length=0.03,fc='r',ec='r')\n",
    "plt.arrow(k02, c0[4], 0,-0.12, head_width=0.03,head_length=0.02,fc='r',ec='r')\n",
    "plt.arrow(k02, c0[4],-0.12, 0, head_width=0.02,head_length=0.03,fc='r',ec='r')\n",
    "### AXIS\n",
    "ax.axvline(0, color=\"k\")                                     # y-axis\n",
    "ax.axhline(0, color=\"k\")                                     # x-axis\n",
    "ax.yaxis.set_major_locator(plt.NullLocator())                # Hide ticks\n",
    "ax.xaxis.set_major_locator(plt.NullLocator())                # Hide ticks\n",
    "### TEXT \n",
    "plt.text(x_max*0.98, -0.05      , r'$\\hat{k}$', color = 'k') # x-axis label\n",
    "plt.text(-0.15 , y_max*0.95     , r'$\\hat{c}$', color = 'k') # y-axis label\n",
    "plt.text(k_star*1.05, y_max*0.95, r'$\\hat{c}=0$', color = \"k\")\n",
    "plt.text(k_zero*0.95, 0.07      , r'$\\hat{k}=0$', color = \"k\")\n",
    "# SETTINGS\n",
    "plt.axis(v)\n",
    "plt.box(False)                                               # Hide axis\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\"8|STABLE-PATH: FORWARD AND BACKWARD SHOOTING\"\n",
    "def forward_shoot(k0, step):\n",
    "    # Set conditions for initial shoot\n",
    "    tol   = 1.0e-10\n",
    "    c0    = consumption(k0) - step\n",
    "    error = np.abs(((k0 - k_star)**2 + (c0 - c_star)**2)**0.5)\n",
    "    path      = np.zeros(shape=(1, 2))\n",
    "    path[0,0] = k0\n",
    "    path[0,1] = c0\n",
    "\n",
    "    # Start loop (forward shoot)\n",
    "    count = 0\n",
    "    while 1:\n",
    "        k1 = path[count, 0]\n",
    "        c1 = path[count, 1]\n",
    "        k2 = k1 + kdot(k1, c1)*0.1\n",
    "        c2 = c1 + cdot(k1, c1)*0.1\n",
    "        error = np.abs(((k2 - k_star)**2 + (c2 - c_star)**2)**0.5)\n",
    "        path = np.append(path, [[k2, c2]], axis = 0)\n",
    "        count = count + 1\n",
    "\n",
    "        # Check if this is the stable-path\n",
    "        if error < tol:\n",
    "            break\n",
    "        \n",
    "        # Set up code break and reset the shoot\n",
    "        if c2 > c_star:\n",
    "            path = np.zeros(shape=(1, 2))\n",
    "            c0 = c0 - step\n",
    "            path[0,0] = k0\n",
    "            path[0,1] = c0\n",
    "            count = 0\n",
    "        \n",
    "        # Break code when k < k*\n",
    "        if k2 > k_star:\n",
    "            break\n",
    "        \n",
    "        # Put a limit of iterations to the code\n",
    "        if count > 100:\n",
    "            break\n",
    "       \n",
    "    return path\n",
    "\n",
    "\n",
    "def backward_shoot(k0, step):\n",
    "    # Set conditions for initial shoot\n",
    "    tol   = 1.0e-10\n",
    "    c0    = consumption(k_gold) * 1.5\n",
    "    error = np.abs(((k0 - k_star)**2 + (c0 - c_star)**2)**0.5)\n",
    "    path      = np.zeros(shape=(1, 2))\n",
    "    path[0,0] = k0\n",
    "    path[0,1] = c0\n",
    "\n",
    "    # Start loop (forward shoot)\n",
    "    count = 0\n",
    "    while 1: # Infinite loop with break conditions inside\n",
    "        k1 = path[count, 0]\n",
    "        c1 = path[count, 1]\n",
    "        k2 = k1 + kdot(k1, c1)*0.01\n",
    "        c2 = c1 + cdot(k1, c1)*0.01\n",
    "        error = np.abs(((k2 - k_star)**2 + (c2 - c_star)**2)**0.5)\n",
    "        path = np.append(path, [[k2, c2]], axis = 0)\n",
    "        count = count + 1\n",
    "\n",
    "        # Check if this is the stable-path\n",
    "        if error < tol:\n",
    "            break\n",
    "        \n",
    "        # Set up code break and reset the shoot\n",
    "        if c2 < c_star:\n",
    "            path = np.zeros(shape=(1, 2))\n",
    "            c0 = c0 + step\n",
    "            path[0,0] = k0\n",
    "            path[0,1] = c0\n",
    "            count = 0\n",
    "        \n",
    "        # Break code when k < k*\n",
    "        if k2 < k_star:\n",
    "            break\n",
    "        \n",
    "        # Put a limit of iterations to the code\n",
    "        if count > 400:\n",
    "            break\n",
    "       \n",
    "    return path\n",
    "\n",
    "stable_L = forward_shoot(k00, 0.000001)\n",
    "stable_R = backward_shoot(k_star*2.00, 0.000001)\n",
    "\n",
    "\n",
    "\n",
    "\"9|PLOT RAMSEY MODEL WITH STABLE PATH\"\n",
    "# Value of k such that c = 0\n",
    "k_zero = (1/(n + g + delta))**(1/(1-alpha))\n",
    "# Axis range\n",
    "y_max = np.max(c_loci)*1.7\n",
    "x_max = k_zero*1.2\n",
    "v = [0, x_max, -0.1, y_max]                        \n",
    "\n",
    "### BUILD PLOT AND POPULIATE WITH LOCI LINES\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "ax.plot(k_array, c_loci, color=\"K\", alpha = 0.8)\n",
    "ax.axvline(k_loci, color=\"K\", alpha = 0.8)\n",
    "### STABLE PATH\n",
    "ax.plot(stable_L[:,0], stable_L[:,1], \"b--\", alpha = 0.7)\n",
    "ax.plot(stable_R[:,0], stable_R[:,1], \"b--\", alpha = 0.7)\n",
    "### AXIS\n",
    "ax.axvline(0, color=\"k\")                                     # y-axis\n",
    "ax.axhline(0, color=\"k\")                                     # x-axis\n",
    "ax.yaxis.set_major_locator(plt.NullLocator())                # Hide ticks\n",
    "ax.xaxis.set_major_locator(plt.NullLocator())                # Hide ticks\n",
    "### TEXT \n",
    "plt.text(x_max*0.98, -0.05      , r'$\\hat{k}$', color = 'k') # x-axis label\n",
    "plt.text(-0.15 , y_max*0.95     , r'$\\hat{c}$', color = 'k') # y-axis label\n",
    "plt.text(k_star*1.05, y_max*0.95, r'$\\hat{c}=0$', color = \"k\")\n",
    "plt.text(k_zero*0.95, 0.07      , r'$\\hat{k}=0$', color = \"k\")\n",
    "plt.text(k_star*1.05, -0.05     , r'$k^*$'      , color = 'k')\n",
    "plt.text(-0.15, c_star          , r'$c^*$'      , color = 'k')\n",
    "ax.axhline(c_star, 0, k_star/v[1], linestyle=\":\", color = \"grey\")\n",
    "# SETTINGS\n",
    "plt.axis(v)\n",
    "plt.box(False)                                               # Hide axis\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}