alanlujan91/HAOS.ipynb

## HAOS.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heterogeneous Agent Optimal Savings I\n",
    "## Linear Interpolation\n",
    "### By Alan Lujan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Huggett, 1993 model is an ideal starting point for developing heterogeneous agent solution algorithms. This notebook presents standard solution algorithms with a few extensions. \n",
    "\n",
    "First, we import relevant python libraries that will be useful in our calculations. \n",
    "- `interpolation:` provides \"jit\" linear interpolation\n",
    "- `matplotlib:` graph plotting library\n",
    "- `numba`: provides just-in-time compilation for speedy calculation\n",
    "- `numpy:` numerical library for working with arrays and vectorized functions\n",
    "- `quantecon:` provides \"jit\" maximum search\n",
    "- `scipy:` provides root finding algorithm\n",
    "- `time:` enables time keeping "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt  # plotting library\n",
    "import numpy as np  # numerical library\n",
    "from interpolation import interp  # linear interpolation\n",
    "from numba import njit  # just-in-time compiling nopython\n",
    "from numba import prange  # parallel range\n",
    "from quantecon.optimize.scalar_maximization import brent_max as argmax  # finding maximum\n",
    "from scipy.optimize import toms748 as root  # finding roots\n",
    "from time import time  # time keeping\n",
    "\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "The most commonly used function for these models is the constant relative risk aversion (CRRA) specification, where the coefficient of risk aversion is $\\sigma$. When $\\sigma=1$, the limit of the function is $u(c)=\\log(c)$, so in our definition we allow for different cases depending on the parameter. We use a function factory to create a jitted function for the model depending on the value of $\\sigma$.\n",
    "\n",
    "\\begin{equation}\n",
    "\t\t\\begin{gathered}\n",
    "\t\t\\mathbb{E}_0 \\left( \\sum_{t=0}^{\\infty} \\beta^t u(c_t) \\right) \\quad \\text{ where } \\beta \\in (0, 1), \\text{ and}  \\\\\n",
    "\t\tu(c) = \\lim_{\\nu \\to \\sigma} \\frac{c^{1-\\nu}-1}{1-\\nu} \\quad \\text{ where } \\sigma \\ge 1.\n",
    "\t\t\\end{gathered}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def utility(σ):\n",
    "\n",
    "    @njit\n",
    "    def log(c):\n",
    "        return np.log(c)\n",
    "\n",
    "    @njit\n",
    "    def iso(c):\n",
    "        return (c**(1 - σ) - 1) / (1 - σ)\n",
    "\n",
    "    return log if σ == 1.0 else iso\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The environment is composed of ex-ante identical households with a labor income endowment $w_i$ and a net bond balance of $b$ (savings if positive, debt if negative). The household can then use their total period budget to choose consumption of a non-storable good and next period debt or savings in the form of discount bonds at the price of $q$. Thus the following recursive value function summarizes the household problem:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "    \\begin{gathered}\n",
    "    V(w_i, b;q) = \\underset{(c,b')\\in \\Gamma(w_i,b;q)}{\\max} u(c) + \\beta \\sum_{i=1}^{Nz} \\pi_{ij} V(w_j,b';q) \\\\\n",
    "    \\text{where} \\\\\n",
    "    \\Gamma(w_i,b;q) = \\{(c,b'): c+qb' \\le w_i + b; c \\ge 0; b' \\ge \\underline{b} \\}\n",
    "    \\end{gathered}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in solving the model is creating an object that contains the model parameters. I take advantage of python's object oriented programming to create a class of moodels named HAOS (heterogeneous agent optimal savings) that takes in the following parameters:\n",
    "\n",
    "- $\\sigma \\ge 1$: CRRA parameter of relative risk aversion \n",
    "- $u$: jitted function to override the CRRA utility form, useful for different utility specifications\n",
    "- $\\beta$: household discount factor for future utility of consumption \n",
    "- $\\omega$: the set of possible labor income endowments\n",
    "- $\\pi$: the transition probability matrix for the Markov chain process of labor income\n",
    "- $\\texttt{b\\_bounds}$: the range of bonds over which we solve the value, policy, and distribution functions, where the lower bound serves as the credit limit\n",
    "- $\\texttt{Nb}$: number of grid points used to solve the value function\n",
    "- $\\texttt{N}\\mu$: number of grid points used to solve the distribution function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heterogeneous agent optimal savings\n",
    "class HAOS(object):\n",
    "\n",
    "    def __init__(\n",
    "            self,\n",
    "            σ=1.5,  # CRRA parameter\n",
    "            u=None,  # user defined utility override\n",
    "            β=0.99322,  # discount factor\n",
    "            ω=[0.1, 1.0],  # labor earnings endowment\n",
    "            π=[[0.5, 0.5], [0.075, 0.925]],  # earnings transition matrix\n",
    "            b_bounds=[-4.0, 4.0],  # credit limits\n",
    "            Nb=100,  # number of value function grid points\n",
    "            Nμ=1000):  # number of distribution grid points\n",
    "\n",
    "        self.σ = σ if u is None else None\n",
    "        self.u = utility(σ) if u is None else u\n",
    "        self.β = β\n",
    "        self.ω, self.π = np.array(ω), np.array(π)\n",
    "        blow, bhigh = b_bounds\n",
    "        self.blow, self.bhigh = b_bounds\n",
    "        self.b_bounds = np.array(b_bounds)\n",
    "        Nz = self.ω.size\n",
    "        self.Nz, self.Nb, self.Nμ = Nz, Nb, Nμ\n",
    "\n",
    "        # grid for value and policy functions\n",
    "        self.bgrid = np.geomspace(1, bhigh - blow + 1, Nb) + blow - 1\n",
    "        # grid for distribution\n",
    "        self.μgrid = np.linspace(blow, bhigh, Nμ)\n",
    "\n",
    "        self.init_value()\n",
    "        self.init_distribution()\n",
    "\n",
    "        self.T, self.G, self.A = None, None, None\n",
    "\n",
    "    def init_value(self):\n",
    "        Nz, Nb = self.Nz, self.Nb\n",
    "        self.v, self.g = _init_value(Nz, Nb)\n",
    "\n",
    "    def init_distribution(self):\n",
    "        Nz, Nμ = self.Nz, self.Nμ\n",
    "        self.μ, self.ind, self.w = _init_distribution(Nz, Nμ)\n",
    "\n",
    "    def set_value(self, v, g):\n",
    "        self.v, self.g = v, g\n",
    "\n",
    "    def set_distribution(self, μ, ind, w):\n",
    "        self.μ, self.ind, self.w = μ, ind, w\n",
    "\n",
    "    def set_solution(self, qstar, Bstar):\n",
    "        self.qstar, self.B = qstar, Bstar\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "@njit\n",
    "def _init_value(Nz, Nb):\n",
    "    v = np.zeros((Nz, Nb))\n",
    "    g = np.zeros((Nz, Nb))\n",
    "    return v, g\n",
    "\n",
    "@njit\n",
    "def _init_distribution(Nz, Nμ):\n",
    "    μ = np.zeros((Nz, Nμ))\n",
    "    mid = int(Nμ / 2)\n",
    "    val = 1.0 / Nz\n",
    "    μ[:, mid] = np.ones(Nz) * val\n",
    "    ind = np.empty((Nz, Nμ), dtype=np.int64)\n",
    "    w = np.empty((Nz, Nμ)) * np.nan\n",
    "    return μ, ind, w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I use a log-spaced grid to solve the value and policy functions, expecting more concavity in the value function near the lower bound, or the credit limit. To solve for the stationary distribution, I use a finer linear grid. Additional helper methods such as $\\texttt{init\\_*}$ and $\\texttt{set\\_*}$ are used to store and clear previous solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Value Function Iteration and the $T$ operator\n",
    "\n",
    "To solve the value function, I use value function iteration (VFI) and the $T$ operator, the properties of which are well established in dynamic programming. To avoid having to pass all the model parameters every time we use the T operator, we can use a function factory as described below. \n",
    "\n",
    "The first component is the $\\texttt{objective}$ function, which evaluates the value function at a given proposed next period bond choice, given a price of bonds $q$ and a total period budget $y$. The $T$ operator then uses a golden section search algorithm with inverse parabolic interpolation (Brent, 1973) to evaluate the value function at a range of next period bond choices to find a maximum. \n",
    "\n",
    "Thus, the $T$ operator takes an initial value function defined over the labor income and the bond grid points to produce a new value function. This operator is designed to allow skipping of the maximization search as described in (Judd, 1998). If we provide an already solved policy function $g$, the operator skips the max search and just evaluates the value of the current state. This option will be used later to speed up VFI. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "    (Tv)(\\omega_i, b) = \\underset{(c,b')\\in \\Gamma(\\omega_i,b;q)}{\\max} u(c) + \\beta \\sum_{i=1}^{Nz} \\pi_{ij} v(\\omega_j,b';q)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def T_operator(model, use_parallel=True, max_tol=1e-5, max_iter=500):\n",
    "\n",
    "    u, β = model.u, model.β\n",
    "    ω, π = model.ω, model.π\n",
    "    blow, bhigh = model.blow, model.bhigh\n",
    "    Nz, Nb = model.Nz, model.Nb\n",
    "    bgrid = model.bgrid\n",
    "\n",
    "    v, g = model.v, model.g\n",
    "\n",
    "    @njit\n",
    "    def objective(bp, q, y, Evz):\n",
    "        c = y - q * bp\n",
    "        if c <= 0.0:\n",
    "            return -1e16\n",
    "        # Ev = np.interp(bp, bgrid, Evz)\n",
    "        Ev = interp(bgrid, Evz, bp)\n",
    "        return u(c) + β * Ev\n",
    "\n",
    "    @njit(parallel=use_parallel)\n",
    "    def T(v, q, g, solve=True):\n",
    "\n",
    "        v_new = np.zeros_like(v)\n",
    "        g_new = np.zeros_like(g) if solve else g\n",
    "\n",
    "        Ev = π @ v\n",
    "\n",
    "        for bi in prange(Nb):\n",
    "            for zi in prange(Nz):\n",
    "                y = bgrid[bi] + ω[zi]\n",
    "                if solve:\n",
    "                    bphigh = np.maximum(blow, np.minimum(bhigh, y / q))\n",
    "                    if bphigh <= blow:\n",
    "                        b_max = blow\n",
    "                        v_max = objective(b_max, q, y, Ev[zi])\n",
    "                    else:\n",
    "                        b_max, v_max, _ = argmax(objective,\n",
    "                                                 blow,\n",
    "                                                 bphigh,\n",
    "                                                 args=(q, y, Ev[zi]),\n",
    "                                                 xtol=max_tol,\n",
    "                                                 maxiter=max_iter)\n",
    "                    v_new[zi, bi] = v_max\n",
    "                    g_new[zi, bi] = b_max\n",
    "                else:\n",
    "                    v_new[zi, bi] = objective(g[zi, bi], q, y, Ev[zi])\n",
    "\n",
    "        return v_new, g_new\n",
    "\n",
    "    T(v, 1.0, g, solve=True)\n",
    "\n",
    "    return T\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is an important extension embedded in the $T$ operator above. Whereas in Huggett 1993 the credit limit is picked such that $\\underline{b} > \\frac{-w_1}{1-q}$ to ensure that any debt can be repaid feasibly, this $T$ operator allows for arbitrarily large credit limits. While this in theory means that agents can borrow more than the maximum amount they can repay, the value function is highly negative when consumption is less than or equal to 0. This induces endogenous and idyosincratic borrowing constraints, such that agents self enforce their borrowing limits even if credit limits are large. \n",
    "\n",
    "The $\\texttt{solve\\_model}$ subroutine solves the value function by VFI. To do so, it repeatedly uses the $T$ operator defined above to derive the optimal policy, and avoids the costly optimal policy computation when the optimal policy function has reached some level of numerical convergence. The $T$ operator was also designed to take advantage of parallelized computation, so this is an option for solving that we use to compare the efficiency of parallel computing. \n",
    "\n",
    "- $\\texttt{use\\_parallel}$: if $\\texttt{True}$, computations are performed in parallel\n",
    "- $\\texttt{use\\_existing}$: if $\\texttt{True}$, the value function uses a pre-solved value function, which might be useful if we have solved the value function before\n",
    "- $\\texttt{tol}$: the subroutine stops when $|Tv-v|\\le \\texttt{tol}$\n",
    "- $\\texttt{max\\_iter}$: the subroutine stops at a maximum number of loops to avoid infinite looping\n",
    "- $\\texttt{skip\\_policy}$: if $\\texttt{True}$, the subroutine skips policy evaluation when $|g_{t+1}-g_{t}| \\le \\texttt{tol}*0.1$\n",
    "- $\\texttt{solve\\_skip}$: when the subroutine skips policy evaluation, it re-evaluates policy after a set number of iterations\n",
    "- $\\texttt{verbose}$: optional print argument to see if VFI error is decreasing (debug)\n",
    "- $\\texttt{print\\_skip}$: avoid printing every single iteration, instead only print error every $\\texttt{print\\_skip}$ iterations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_value(model,\n",
    "                T,\n",
    "                q,\n",
    "                use_existing=False,\n",
    "                tol=1e-4,\n",
    "                max_iter=2000,\n",
    "                skip_policy=False,\n",
    "                solve_skip=100,\n",
    "                verbose=True,\n",
    "                print_skip=100):\n",
    "\n",
    "    start = time()\n",
    "\n",
    "    if not use_existing:\n",
    "        model.init_value()\n",
    "\n",
    "    v, g = model.v, model.g\n",
    "\n",
    "    i = 0\n",
    "    skipped = 0\n",
    "    solved = 0\n",
    "    error = tol + 1.0\n",
    "    g_err = tol + 1.0\n",
    "    solve = True\n",
    "\n",
    "    while i < max_iter and (error > tol or g_err > tol * 0.1):\n",
    "\n",
    "        v_new, g_new = T(v, q, g, solve)\n",
    "        error = np.max(np.abs(v - v_new))\n",
    "        g_err = np.max(np.abs(g - g_new))\n",
    "        v, g = v_new, g_new\n",
    "        i += 1\n",
    "\n",
    "        if solve:\n",
    "            solved += 1\n",
    "\n",
    "        if skip_policy:\n",
    "\n",
    "            if g_err < tol * 0.1 and skipped < solve_skip and error > 2 * tol:\n",
    "                solve = False\n",
    "                skipped += 1\n",
    "            else:\n",
    "                solve = True\n",
    "                skipped = 0\n",
    "\n",
    "        if verbose and i % print_skip == 0:\n",
    "            print(f\"Value error at iteration {i} is {error:.4e}.\")\n",
    "            if solve:\n",
    "                print(f\"-- Policy error at iteration {solved} is {g_err:.4e}.\")\n",
    "\n",
    "    end = time()\n",
    "\n",
    "    runtime = end - start\n",
    "\n",
    "    if i == max_iter:\n",
    "        print(\"Failed to converge!\")\n",
    "\n",
    "    if verbose and i < max_iter:\n",
    "        print(f\"Value error at iteration {i} is {error:.4e}.\")\n",
    "        print(f\"-- Policy error at iteration {solved} is {g_err:.4e}.\\n\")\n",
    "        print(f\"Value converged in {i} iterations and {runtime:.4f} seconds.\\n\")\n",
    "\n",
    "    model.set_value(v, g)\n",
    "\n",
    "    return v, g\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start calculation, we instanciate a HAOS object called model with default parameters. Notice that the HAOS class is flexible an allows an override of parameters to experiment with different specifications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = HAOS()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to solve for the value function, we create a $T$ operator for the model and use parallel computation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tp = T_operator(model, use_parallel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now solve the value function given a price of bonds. For simplicity, lets use $q=1$ and use the $\\texttt{skip\\_policy}$ option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Value error at iteration 100 is 8.4790e-02.\n",
      "-- Policy error at iteration 100 is 1.2247e-04.\n",
      "Value error at iteration 200 is 4.2858e-02.\n",
      "Value error at iteration 300 is 2.1705e-02.\n",
      "Value error at iteration 400 is 1.0993e-02.\n",
      "Value error at iteration 500 is 5.5675e-03.\n",
      "Value error at iteration 600 is 2.8197e-03.\n",
      "Value error at iteration 700 is 1.4281e-03.\n",
      "Value error at iteration 800 is 7.2326e-04.\n",
      "Value error at iteration 900 is 3.6630e-04.\n",
      "Value error at iteration 1000 is 1.8552e-04.\n",
      "-- Policy error at iteration 166 is 2.6750e-10.\n",
      "Value error at iteration 1091 is 9.9891e-05.\n",
      "-- Policy error at iteration 257 is 2.1690e-07.\n",
      "\n",
      "Value converged in 1091 iterations and 0.1107 seconds.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "v, g = solve_value(model, Tp, 1.0, skip_policy=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxV1bn/8c+TMIR5ykQgYZ6nAJFBVERQAVEc61S1Wmu1eq+9t71ay7V6a3t/tvbW/lpbW2tbO2CttSooIqOzooQ5QJinkJCEBJKQkOGc8/z+WCclP24Yz0n2yTnP+/XKi5x99j7ncWB/915r7bVEVTHGGBO74rwuwBhjjLcsCIwxJsZZEBhjTIyzIDDGmBhnQWCMMTGuldcFnI/ExETt27ev12UYY0yLsmbNmsOqmnTy9hYZBH379iU7O9vrMowxpkURkX2NbbemIWOMiXEWBMYYE+MsCIwxJsZZEBhjTIyzIDDGmBhnQWCMMTHOgsAYY2KcBYExxkQ6Xw3s/gCWPwnlBWH/+Bb5QJkxxkQ1VSjcDLvfg13vwb5PwXcc4lpB+iTo3DOsX2dBYIwxkaDiEOx+H3atdH8eK3TbEwfDuDthwDToexG07RT2r7YgMMYYL9RWwf5P3RX/rvegaLPb3j4R+l/qTvz9p0GXXk1eigWBMcY0h0AACnPcFf+ulbB/FfhrIL4t9JkMo/8LBlwGKSMhrnm7by0IjDGmqVQcCl7xr3Tt/ZXFbnvyCJjwNXfVn3EhtGnvaZkWBMYYEy511cHmnpUuAApz3Pb2ie6kP2C6a/YJc2dvqCwIjDHmfKlC8TbYtQJ2roB9n4CvGuJaQ8YkmPFksLlnVLM395wLCwJjjDkXVaXB0T0r3FV/+UG3vccgGP8Vd9Xfdwq06eBllefEgsAYY07H74ODa05c9eevBQ1AQpfg6J5H3FV/1wyvKz1vFgTGGHOysjx30t+1wl39V5eBxEGv8XDJIzBwOqSNg/joOIVGxz+FMcaEor6Td2fwqr94q9veKQ2GXX2ik7d9dy+rbDIhBYGI3AQ8CQwDJqhqdnD77cB/NNh1NDBOVdefdPyTwNeA4Jgqvquq74RSkzHGnJEqlOyCncvdz96P3RQO8W2gzxQYezsMnAFJQ0HE62qbXKh3BDnA9cBvGm5U1fnAfAARGQUsODkEGnhWVX8SYh3GGHN6Ncdgz4cnTv5Hg+u49xjopnAYOKPFdfKGS0hBoKpbAeT0iXkr8NdQvscYY86ZKhRthZ3L3Il/32cQqIPWHaD/VJjyr67Jp3s/ryv1XHP0EdwMzD3N+w+JyJ1ANvAtVT3S2E4ich9wH0BGRsvtnTfGNKHqMte5u3O5a+uvH9qZPBwm3Q8DL3fj+1u19bTMSHPGIBCR5UBqI2/NU9UFZzh2IlClqjmn2OV54ClAg3/+D3BPYzuq6gvACwBZWVl6prqNMTGgfrrmHUvdyf/A5xDwQdvOrnN36qOuyacZJm5ryc4YBKo6I4TPv4XTNAupamH97yLyW+DtEL7LGBMLqsvdVf+Ope6qvyLfbU8ZBRf+i7vqT58A8a09LTNc6vwBsvce4f3tRbyfW8yzN2cyPK1zWL+jyZqGRCQOuAm45DT79FTV+uV2rsN1PhtjzAn1bf31V/37Pztx1T9gmjvxD5wRcfP3hKKwvJoPthXz3rYiPt5xmIoaH63jhQn9ulPj84f9+0IdPnod8AsgCVgkIutV9crg25cAeaq6+6RjXgR+HRxq+mMRycQ1De0Fvh5KPcaYKFFzDPZ84E7+O5ZDeZ7bnjIyKq/6/QFl3f4jvLetiPdyi9lSUA5AaucE5ozpyaVDkpkyMJGObZvm2l1UW15ze1ZWlmZnZ3tdhjEmXFShZGfwxL/ULc3or4U2nWDApSeu+qOorb+0spYPtrsT/4c7ijlaVUd8nDA+oxvThiYzbWgSQ1I6nWlU5jkRkTWqmnXydnuy2Bjjjbpq9yDXjqWwYwkc2eu2Jw2FiV+HQVe49XlbtfG0zHBRVTbnl7Myt4j3thWx/sBRVCGxYxumD01h2tAkLh6URJd2zX+XY0FgjGk+ZXmwfQnsWOaafuqqoFU76HcxTH7Infy79fG6yrA5VuPj4x2HeS948i+qqEEERvfuysPTBzFtSDKjenUhLs7bp5ctCIwxTcfvg7zV7op/+9IT6/J2zYDM22HwlW5B9tbtvK0zjPYermRlbhErc4v4fE8JdX6lU9tWXDI4iWlDk7l0SBKJHSPrOQYLAmNMeFWVumGdO5a4UT7Hj0BcK8iYDJc/5U7+iYOjZg6fOn+A1XtLeS+3iBW5RewurgRgQFIHvnJhX6YNTeaCvt1pHW8L0xhjopUqFOe6Jp/tS+DAKjdff/tEGDwLBl/h5utP6OJ1pWFTWlnL+9vcif/D7cVUVPtoEx/HxP7duXNSHy4bmkJGD2/XIT4XFgTGmHPnq3EdvdvfdT9H97vtqaPh4m/BoCvd3P0RvDzjuVBVthceY0VuISu2FrF2/xFUIalTW2aP7Mllw5K5aGAiHZpoeGdTa5lVG2Oa37HiYFv/u26JxtpjrqO3/1S46N9dR28UDe+s8flZtbuUlVsLWZFbRN6R4wCM6tWFf71sENOHJTMyzfuO3nCwIDDGNK5+Hp/ti2Hbu265RhQ694LRX3LNPv0ujqqO3pJjNazMLWLF1iI+2lFMZa2fhNZxXDQwiQenDeSyocmkdE7wusywsyAwxpzgq4E9H7mT//YlUHbAbe81HqbNgyEz3dO9UdLRq6rsKDrG8q2FLN9SyLrg2P7UzgnMHduLGcOSuXBAIgmt470utUlZEBgT6yoPBzt6F8POlVBXCa3buw7eqY+49v5OKV5XGTZ1/gCr95SyfGsRy7cWsr+0CnBNPg9PH8SMYSmMSOsc1id6I50FgTGxRhUOb4dt78C2xXDgC0Dd+rxjbg42+VwCraOnCaTseB0fbC9m+ZZC3ttW5Eb5tIpjyoAefH1qf6YPTSG1S/T8854rCwJjYoHf54Z1blvsAqA0OBdkzzFuzv4hs9zvUXQVfKC0ihVbC1m2tZDPd5fiCyg9OrRh5ohUZgxP4eJBibRvY6dAsCAwJnrVVLgHu7a945p+qo+6xdn7TXXTOQyeGVWjfFSVnIPlLNtyiKVbCsk9VAHAwOSO3Htxfy4fnkxmejfio2CUT7hZEBgTTcoLgk0+77iF2v210K6bu+IfMsu1+7ft5HWVYVPrC7BqdwnLthSyfGshBWXVxAlk9enOvNnDmDE8hX6JsbcY/bmyIDCmJatftGXbIsh9B/LXuu3d+sGE+2DIbEifCPHR81e9vLqO97cVs2xLIe/nFlFR46Nd63guGZzIt64YwmVDk+neITpmLG0u0fN/hzGxIuB3a/PmLoLct09M39wrCy57HIZe5aZyjqL2/qLyapZuKWTplkI+23WYOr9r7589qieXD0/hokHRP8SzKVkQGNMS1B13T/PmLnLDPKtKTrT3T3nYjfSJoqUaAXYVH2Pp5kKWbjnEuv1HAejToz13T+nH5cNTGJdh7f3hYkFgTKSqKnWLtmx9C3atdHP3t+3iJnEbepVbsSuK2vtVlY15ZSzZ7Dp7dxYdA2B07y58+4rBXDEilUHJHWNqfH9zsSAwJpKUHQw2+bwFez8B9bvx/Zm3wdA50GdK1KzYBe7hri/2lLI0ePIvKKsmPk6Y2K87d0zqw+XDU0jrGj1TWESqkINARJ4BrgZqgV3A3ap6NPjeY8BXAT/wr6q6pJHj+wGvAN2BtcAdqlobal3GtBjF29xVf+7bkL/ObUsc4pp8hs6BtLFRM4snQHWdnw+3F7NkcyErcgs5WlVHQus4LhmUxLevGML0Ycl0bR89YdcShOOOYBnwmKr6RORHwGPAoyIyHLgFGAGkActFZLCq+k86/kfAs6r6ioj8Ghccz4ehLmMik6ob3bP1bRcAJTvc9l7jYfoTMOxqSBzkbY1hVlFdx8rcIt7NOcT724o5Xuenc0IrZgxL4YoRqUwdnES7NtbZ65WQg0BVlzZ4uQq4Mfj7XOAVVa0B9ojITmAC8Fn9zuIa+y4Dbgtu+iPwJBYEJtr4fbD/M3fVv/VtKM8DiXfLNE78umvz75zmdZVhVVpZy7Ith3g35xCf7Cyh1h8gqVNbrh/Xi5kjU5nUv0dEr9oVS8LdR3AP8Lfg771wwVAvL7itoR7AUVX1nWYfAETkPuA+gIyMjHDVa0zT8dW4h7q2LnTt/lUl0CrBPdR12Tz3ZG/77l5XGVaF5dUs2XyIxZsO8fmeEgIKvbu1487JfZg1KpWx6d2iYv7+aHNWQSAiy4HURt6ap6oLgvvMA3zA/PrDGtlfT/7os9jHbVR9AXgBICsrq9F9jPFcbZVbp3frQjetQ005tOnkRvoMuyY40qej11WG1YHSKnfyzznEmn1HALde7zcuHcjMkakxN5NnS3RWQaCqM073vojcBcwBpqtq/Uk6D0hvsFtvIP+kQw8DXUWkVfCuoLF9jIls1eVumOeWBS4E6qqgXXcYfo07+febGlUzeQLsPVzJOzkFvJtziI15ZQAM79mZf798MLNGpjIoJXqGtcaCcIwamgk8CkxV1aoGby0EXhaRn+I6iwcBXzQ8VlVVRN7D9Su8AtwFLAi1JmOa3PEjbibPLQth1wo3p0/HFBhzqwuAPhdF1bQOADuLjrF4UwHv5Bxia0E5AGPSu/KdWUOZNTKVPj1sTp+WKhz/pz4HtAWWBW//Vqnq/aq6WUReBbbgmowerB8xJCLvAPeqaj4uRF4RkR8A64DfhaEmY8KvssTN6bP5TdjzAQR80Lk3XHCvu/JPnxhVwzwBdhRWsGhTAe9sKmB7oXvAK6tPNx6fM5yZI1PpZWP8o4KcaMlpObKysjQ7O9vrMkwsOFbsHu7assAt4ah+6NbXnfiHXwu9xkXVnD4A2wsreHujO/nvLDqGCFzQtzuzR6Yyc2TPmF7ApaUTkTWqmnXy9ui6dzUmHI4Vuc7ezW/Cvk9AA9C9v3vAa8S1kDo6Kk/+izYWsCh48o8TmNCvO3dOHsHMEakkR+GC7eYECwJjACoK3cl/y4ITJ/8eg+Dib7kr/5QRUXfy31nkrvwXbSxgR/DKf2K/7tw1eQRXjkwluZOd/GOFBYGJXceKg1f+b5w4+ScOhkv+w538k4dF3cl/d/ExFm0s4O2NBWwrrPhns89Tc+3kH8ssCExsqSw5cfLf+1GDK/9vw4jrovLkf6C0irc25vP2hgK2BEf7ZPXpxpNXD2f2qJ7W7GMsCEwMqCp1UztsfgN2f+A6fLsPcM0+I66D5OFRd/IvKDvOoo0FvLUhnw3Bcf6Z6V35z6uGcdXonvTsYqN9zAkWBCY6VZe5pRs3v+7m8g/43GifKQ/DyOshZWTUnfwPH6th8aYCFm7IZ/Ve94TviLTOfGfWUK4a1ZP07u09rtBEKgsCEz1qK920Djn/gB3LwF8DXdJh0gMw4no3nXOUnfzLjtexZPMh3tqQz6e7SvAHlIHJHfn3ywczZ3RP+idF13QWpmlYEJiWzVcDO1e4k/+2xVBXCR1TIesed+Xf+4KoO/kfr/WzIreQhevzeX9bMbX+AOnd2/H1S/pzTWYaQ1I62dw+5pxYEJiWJ+B3s3rmvObm868uc3P7jL4JRt4IfS6EuOia277OH+CTnYdZuD6fJZsPUVnrJ6lTW26flME1Y9LITO9qJ39z3iwITMugCnnZ7uSf8zpUFkGbjm4Fr1E3Qv9LIb6111WGlaqydv8RFqzPZ9HGAkoqa+mc0Io5o9OYm5nGxP49bPF2ExYWBCayFeXCpr+7ADiyF+LbuimdR94Ig6+E1tE3+mVnUQVvrstnwYaDHCg9TttWccwYnsLcMWlMHZJE21bRdbdjvGdBYCJP2UF34t/0dzi0CSTOTeV8ySMwbA4kdPG6wrArLK9m4fp83lx/kM355cQJTBmYyMPTB3PliBQ6JUTX3Y6JLBYEJjIcP+oe9Nr4Kuz9GFC3hu/MH7mx/p1SvK4w7CprfLybc4g31x/kk52HCSiM7t2F780ZzpwxPe0pX9NsLAiMd3w1bpjnxr+5YZ/+Gveg16WPuXb/HgO8rjDsfP4AH+88zBvrDrJ0cyHH6/ykd2/HQ9MGMndsLwbYcE/jAQsC07xU4cDnsOEV96Rv9VHokOSGe46+CdKib1pngC355by+No8FG/IprqihS7vWXDeuF9eP7cX4Pt1sxI/xlAWBaR4lu9yV/4ZX4Og+aNXOtfePvhn6T4u61bwAisqrWbA+n3+szSP3UAWt44VpQ5K5flxvpg21Tl8TOaLvb5+JHFWl7qp/wyuQ9wUg0H+qa/oZNgfaRt+6ttV1fpZtKeQfa/P4cHsxAXVz/Dw1dwRzRqfRrUMbr0s05n+xIDDh5at1C7hv+Ctsf9et5Zs8HC7/Poy6CTqneV1h2Kkq6w4c5bU1eby1IZ+Kah9pXRJ44NIBXD+ut7X7m4hnQWBCpwoFG9zJf9PfoarEtftfcK9bzD11VFS2+x8qq+b1dXm8tiaP3cWVJLSOY/bIntwwvjeT+/cgzh72Mi1ESEEgIs8AVwO1wC7gblU9KiKXA08DbYLv/Yeqrmzk+CeBrwHFwU3fVdV3QqnJNKNjRW645/qXoWgzxLeBIbMh8zYYcFnUPekLUOPzs3xLEa9mH+CjHa7pZ0Lf7tx/yQBmjUq18f6mRQr1jmAZ8Jiq+kTkR8BjwKPAYeBqVc0XkZHAEqDXKT7jWVX9SYh1mObiq3VNPutfhh1L3dz+vcbDVT91k7y16+Z1hU1iS345r2Yf4M31BzlaVUfPLgl849KB3Di+N30TO3hdnjEhCSkIVHVpg5ergBuD29c12L4ZSBCRtqpaE8r3GQ8d2gTr5ruRP8dL3QyfFz4EmbdD0hCvq2sSZcfrWLj+IH/LPkDOwXLaxMdxxYgUvpSVzpSBiTbPj4ka4ewjuAf4WyPbbwDWnSYEHhKRO4Fs4FuqeqSxnUTkPuA+gIyMjDCUa87o+BHY9Bqs+wsUrA82/cyCsXdE7ZBPVWXV7lL+tno/i3MOUeMLMLxnZ568ejjXju1F1/Y26sdEH1HV0+8gshxIbeSteaq6ILjPPCALuF4bfKCIjAAWAleo6q5GPjsF14ykwFNAT1W950xFZ2VlaXZ29pl2M+cjEHBr+a77s5vi2VcNKaNg3B1u1E/77l5X2CSKKqr5x5qD/G31fvaWVNEpoRVzM9O45YIMRvaKvrmNTGwSkTWqmnXy9jNe0qnqjDN88F3AHGD6SSHQG3gDuLOxEAh+dmGD/X8LvH2mekwTKTvo2v3X/dk98JXQBcZ+2V39p2V6XV2T8AeUj3YU89cv9rNiaxG+gDKhX3f+dfogZo3sSbs29sCXiQ2hjhqaiescnqqqVQ22dwUW4TqSPznN8T1VtSD48jogJ5R6zDny17k5ftb+CXYuAw1A34vhssfdA19ROMUzuJk+X119gFdWH+Dg0eP06NCGey7qx80XpNuYfxOTQm3kfQ5oCywLzpWySlXvBx4CBgKPi8jjwX2vUNUiEXkR+LWqZgM/FpFMXNPQXuDrIdZjzkbpHnflv24+HDvkOn6nfNM1/3Tv73V1TSIQUD7cUczLn+9nRW4R/oBy0cBEvjt7GJcPT6FNqzivSzTGM2fsI4hE1kdwHvx1kLsI1rwEu99zc/wPugLG3eX+jMKOX4DDx2p4NfsAL3++n7wj7ur/pqx0bp2QTp8eNuzTxJbz7iMwLdyRvbDmj27kT2URdO7t5voZewd0OdWjHS2bqvLFnlL+8vl+3s0poM6vTOrfnUdmDuXKESk22ZsxJ7EgiEZ+n3vYK/t3sHOFm95h0JWQdTcMnBF1C7vXq6iu4811B/nzqn1sLzxG54RW3DGpL7dNzGBgsrX9G3MqFgTRpOKQ6/hd8xKUH4ROPWHqo67tv0tvr6trMjsKK/jTZ/t4fW0elbV+RvXqwo9vGM3VY9Js5I8xZ8GCoKVTdUs7rv6t6wMI+Nw8P7N+DINnRm3bv88fYPnWQv746T4+211Cm1ZxzBndkzsn9yUzvavX5RnTokTnWSIWVJe76R6++C0c3gYJXWHi/W6lryhc4rFeaWUtr6zez/xV+zl49Di9urbjkZlDuDkrnR4d23pdnjEtkgVBS1O8Hb54wU35XHsMembC3F+5Cd+idNw/wNaCcl76ZC9vrj9IjS/A5P49+N7Vw5kxLMXm/DEmRBYELUHA7zp/P/+NG/oZ3wZGXA8T7oPe472ursn4A8rK3CJ+//EePttdQkLrOG4Y35u7JvdlSGr0rW5mjFcsCCJZdZkb9vnFC24YaKc0uOw/YdxXoGOS19U1mWM1Pl7LPsAfPt3LvpIq0rok8J1ZQ7nlgnSb9M2YJmBBEIlKdsHnv3Zz/9Qeg/RJMONJGDonKhd7qVdQdpyXPtnLy1/sp6Lax7iMrjxypRv73yrenvw1pqlYEEQKVdjzIaz6lZv/J741jLzBdQBH6aRv9XIOlvHbj3azaGMBAVVmjerJvRf1Y2xGdC5yY0yksSDwmq8Wcv4Bn/0SCjdB+0SY+ghkfRU6pXhdXZNRVd7fXsxvP9zNp7tK6Ni2FXdd2JevXNiX9O7tvS7PmJhiQeCVqlJY8wf4/AU38VvSMLjmFzDqS9A6wevqmkydP8BbG/J54cPd5B6qILVzAt+dPZRbJmTQ2db7NcYTFgTN7cg+1/yz9s9QV+lW+rr2lzBgupsKIkpV1fp45YsDvPjRbvLLqhmc0pH/uWkMV49Js5k/jfGYBUFzKdgAn/wcNr/hTvgjb4QL/wVSR3pdWZM6WlXLS5/u5aVP93K0qo4J/brzg+tGMm1IMhLFwWdMS2JB0JRU3bKPHz8Lu1ZCm04w+Rsw8YGonfmzXlF5NS9+vIf5q/ZRWetnxrBkHrh0AOP7ROdSl8a0ZBYETSEQgO2L4aP/gYNroEMyTH/CTf/QLrrnwck7UsVvPtjN37IP4PMHuGZMGg9cOtAeADMmglkQhJPfB5tfh49+CsVboWsfuOqnkHl7VHcAA+wvqeKX7+3kH2vzEIEbx/fm/qkDbPEXY1oAC4Jw8NXCxldcABzZ40YAXf8ijLguamf/rLf3cCXPvbeTN9YdJD5OuH1iBl+fOoC0rtE775Ex0SbUxeufAa4GaoFdwN2qelRE+gJbgW3BXevXMj75+O7A34C+uDWLv6SqR0KpqVn5atzavx//DMoOuAngbv4LDLkK4qJ7JMy+kkp+sdIFQKs44a7Jfbl/an+SO0f3nY8x0SjUy9VlwGOq6hORHwGPAY8G39ulqmd6JPY7wApVfVpEvhN8/egZjvFefQB89FO3AEzvC2DOs271rygfCZN3pIpfrNjJa2vzTgTApf1J7mQBYExLFVIQqOrSBi9XATee40fMBS4N/v5H4H0iOQh8tbB+Pnz4EyjPc3MAzX3OPQsQ5QFQVF7NL1bu5JXV+xGEOyb14RuXDrA7AGOiQDgbsO/BNfPU6yci64By4D9V9aNGjklR1QIAVS0QkeQw1hM+fp/rA/jgR3B0v7sDmPuLmAiAo1W1PP/BLv746V58fuVLF6Tz0LSB1gdgTBQ5YxCIyHIgtZG35qnqguA+8wAfMD/4XgGQoaolIjIeeFNERqhq+fkWKiL3AfcBZGRknO/HnJtAALa8Ce/9N5TscH0AV/00JpqAqmp9/OGTvfz6g10cq/FxbWYvvjljkI0CMiYKnTEIVHXG6d4XkbuAOcB0VdXgMTVATfD3NSKyCxgMZJ90eKGI9AzeDfQEik5TxwvACwBZWVl6prpDouoeAFvxX+6J4KRhrhN46JyoDwCfP8Cr2Xn8bPl2iipqmDEshW9fOZihqZ29Ls0Y00RCHTU0E9emP1VVqxpsTwJKVdUvIv2BQcDuRj5iIXAX8HTwzwWh1BMWB9fC8ifclNBdM+C638ComyAu3uvKmpSqsnxrEU8v3squ4krG9+nGr24fR1ZfexLYmGgXah/Bc0BbYFlw3pj6YaKXAN8XER/gB+5X1VIAEXkR+LWqZuMC4FUR+SqwH7gpxHrO35F9sOL7kPMatO8BM38EWXdDq+hfEH1TXhk/WLSFz/eU0j+pA7+5YzxXDE+xuYCMiRESbM1pUbKysjQ7++RWpvNUXeamglj1PEg8TH4QpjwMCdHfFHKorJofL8nl9bUH6dGhDd+8fDC3XJBOa1sNzJioJCJrVDXr5O3R/djr6QT8sPZPsPIHUFUCY26F6Y9D5zSvK2ty1XV+fvvhbn71/i78AeX+qQN4cNoAOtl6AMbEpNgMgr2fwOJH3YpgGRfCzP+GtLFeV9XkVJV3cw7xg0VbOXj0OLNHpfLYrGG2IpgxMS62gqA8H5Y+7voBuqTDjX9w8wHFQFv4zqIKnli4mU92ljA0tRN//dokJg/o4XVZxpgIEFtBsPy/YOtbMPVRmPJNaBP9V8KVNT5+vmIHv/t4D+3bxPP9uSO4bUIGrawfwBgTFFtBMONJuPRR6N7f60qanKqydEshTy7cTEFZNV/K6s2jM4fSo2P0j4Iyxpyb2AqCzj29rqBZ5B89zvcW5LB8axFDUzvx3G1jbWUwY8wpxVYQRDl/QPnTZ3v5yZJtBBS+O3sod0/pZ8NBjTGnZUEQJXYWVfDIaxtZu/8olwxO4ofXjrTRQMaYs2JB0ML5/AFe+Gg3P1u2g/Zt43n25jFcm9nLngo2xpw1C4IWbEdhBd/++wY25JUxe1Qq3587kkTrDDbGnCMLghYoEFB+/8kefrxkGx3axPPL28Zx1ejY6Ag3xoSfBUELk3/0OP/+6npW7S5lxrAU/s/1o0jqZHcBxpjzZ0HQgry9MZ/HXt+EP6D86IZRfCkr3foCjDEhsyBoAapqfTy5cDOvZueRmd6V/3tLpq0UZowJGwuCCJd7qJyHXl7HruJjPDhtAN+cMdieCzDGhJUFQQR7dfUBHl+QQ+d2rfnLVycyZWCi1yUZY6KQBUEEqq7z8/ibOfx9TR5TBvbgZzePtQ5hY0yTsSCIMAdKq3hg/hpyDpbzL5cN5JszBhMfZx3CxpimY0EQQT7bVcI35q/BF1B+d1cW04eleF2SMSYGhNTrKCLPiEiuiGwUkTdEpBFHRQsAAA5HSURBVGtw++0isr7BT0BEMhs5/kkROdhgv9mh1NOS/XnVPr78u8/p0bEtCx+6yELAGNNsQh1+sgwYqaqjge3AYwCqOl9VM1U1E7gD2Kuq60/xGc/W76uq74RYT4vj8wd4YkEOj7+Zw9TBSbzxjQvpl2hDQ40xzSekpiFVXdrg5SrgxkZ2uxX4ayjfE60qa3z8y1/XsTK3iK9d3I/vzBpm/QHGmGYXzgHp9wCLG9l+M6cPgoeCTUu/F5Fup9pJRO4TkWwRyS4uLg61Vs8VVVRz8wuf8cH2Yn5w7UjmXTXcQsAY44kzBoGILBeRnEZ+5jbYZx7gA+afdOxEoEpVc07x8c8DA4BMoAD4n1PVoaovqGqWqmYlJSWd+Z8sgu05XMkNz3/KrqJKXrwziy9P6uN1ScaYGHbGpiFVnXG690XkLmAOMF1V9aS3b+E0dwOqWtjgc34LvH2melq6nINlfOUPX+APKH+9bxKZ6V29LskYE+NC6iMQkZnAo8BUVa066b044CbgktMc31NVC4IvrwNOdecQFdbsK+Urv19Np4RW/Om+iQxM7uh1ScYYE3IfwXNAJ2BZcPjnrxu8dwmQp6q7Gx4gIi+KSFbw5Y9FZJOIbASmAf8WYj0R67NdJdzxuy9I7NSWvz9woYWAMSZihDpqaOBp3nsfmNTI9nsb/H5HKN/fUnyy8zBf/eNq0ru1Z/7XJpLcKcHrkowx5p/syeIm9umuw9zz0mr69ujA/K9NtKUkjTERx4KgCa3eW8pXX8qmT4/2vPy1ifSwEDDGRCCb2L6J5Bws454/rKZnlwTm3zvJQsAYE7EsCJrA7uJj3PX7L9w6AvdOtCmkjTERzYIgzIrKq7nz918A8Jd7J5LWtZ3HFRljzOlZH0EYVdb4uOePqymtrOWV+ybZ5HHGmBbB7gjCJBBQHn5lPVvyy3nutrGM7m1PDBtjWgYLgjB5Zuk2lm8t5ImrR3DZUFtLwBjTclgQhMFbG/J5/v1d3DYxgzsn2wRyxpiWxYIgRNsOVfDIaxvJ6tONJ68egYhNJW2MaVksCEJwrMbHA/PX0DGhFb+6fRxtWtm/TmNMy2NnrvOkqnz39U3sPVzJz28ZS3Jnmz/IGNMyWRCcp7+vyWPhhny+OWMwkwf08LocY4w5bxYE52FfSSVPLtzMxH7deXDaKSdgNcaYFsGC4Bz5A8q3Xt1AfJzw7M2Zts6wMabFsyeLz9HvP95D9r4jPHvzGJs+whgTFeyO4BzsOVzJT5ZuY8awFK7N7OV1OcYYExYWBGdJVXns9Y20aRXHD68bac8LGGOihgXBWXptTR6rdpfy2KxhpNhQUWNMFAk5CETkKRHZGFy8fqmIpAW3i4j8XER2Bt8fd4rjxwcXsN8Z3D/iLrXLqup4enEu4/t045YL0r0uxxhjwiocdwTPqOpoVc0E3ga+F9w+CxgU/LkPeP4Uxz8ffL9+35lhqCmsnl2+nSNVtXx/7gjibJSQMSbKhBwEqlre4GUHQIO/zwX+pM4qoKuI9Gx4bPB1Z1X9TFUV+BNwbag1hdOu4mP8edU+bp2QwYi0Ll6XY4wxYReW4aMi8kPgTqAMmBbc3As40GC3vOC2ggbbegW3n7xPY99xH+7OgYyMjHCUfVaeXpxLu9bx/Nvlg5vtO40xpjmd1R2BiCwXkZxGfuYCqOo8VU0H5gMP1R/WyEfpSa/PZh+C3/GCqmapalZSUtLZlB2yNftKWbalkAcuHUCiLT5vjIlSZ3VHoKozzvLzXgYWAU/gru4b9qz2BvJP2j8vuP10+3jmmSXbSOzYlrun9PW6FGOMaTLhGDU0qMHLa4Dc4O8LgTuDo4cmAWWq2rBZiODrChGZFBwtdCewINSawmHV7hJW7S7lwWkDaN/GHsA2xkSvcJzhnhaRIUAA2AfcH9z+DjAb2AlUAXfXHyAi64OjjAAeAF4C2gGLgz+e+/mKHSR1asutE5qvP8IYY7wQchCo6g2n2K7Ag6d4L7PB79nAyFDrCKd1+4/w6a4S5s0eRkLreK/LMcaYJmVPFjfiNx/spku71tw20e4GjDHRz4LgJHsPV7JkyyG+PCmDDm2tb8AYE/0sCE7y51X7iBfhrsl9vS7FGGOahQVBA1W1Pl7NPsCsUT1tDWJjTMywIGjgrQ35VFT7uHNyH69LMcaYZmNB0MDLn+9ncEpHsvp087oUY4xpNhYEQVvyy9mQV8atEzJs0RljTEyxIAj6+5oDtImP47qxtgSlMSa2WBAAtb4AC9bnc/nwFLq2b+N1OcYY06wsCICPdhRTWlnL9ePsbsAYE3ssCIB3Nh2ic0IrLh7UPNNbG2NMJIn5IKj1BVi25RCXD0+lTauY/9dhjIlBMX/mW7W7hPJqH7NGpnpdijHGeCLmg2D51kISWsdx0aBEr0sxxhhPxHQQqCrLtxRy8aAkm27aGBOzYjoIthVWkF9WzYxhyV6XYowxnonpIPh4x2EAGy1kjIlpsR0EOw/TP6kDaV3beV2KMcZ4JmaDoNYX4PPdpVw00DqJjTGxLaQgEJGnRGSjiKwXkaUikhbcfntw+0YR+VRExpzi+JdEZE/w+PUiktnYfk0hJ7+M43V+Jvfv0VxfaYwxESnUO4JnVHV0cDH6t4HvBbfvAaaq6mjgKeCF03zGf6hqZvBnfYj1nLXVe0oByOrbvbm+0hhjIlJIi/KqanmDlx0ADW7/tMH2VUDvUL6nKazeW0q/xA4kdWrrdSnGGOOpkPsIROSHInIAuJ0TdwQNfRVYfJqP+GGwCelZETnlWVlE7hORbBHJLi4uDqlmVSV73xFbgMYYYziLIBCR5SKS08jPXABVnaeq6cB84KGTjp2GC4JHT/HxjwFDgQuA7qfZD1V9QVWzVDUrKSm04Z57S6o4WlXHeAsCY4w5c9OQqs44y896GVgEPAEgIqOBF4FZqlpyis8uCP5aIyJ/AL59lt8VkvUHjgAwJr1rc3ydMcZEtFBHDQ1q8PIaIDe4PQN4HbhDVbef5viewT8FuBbICaWes7XhQBntWsczKLljc3ydMcZEtJA6i4GnRWQIEAD2AfcHt38P6AH8Krj+r09VswBE5B3gXlXNB+aLSBIgwPoGxzep9QeOMqp3F1rFx+xjFMYY80+hjhq64RTb7wXuPcV7sxv8flko338+AgFl26EKbpmQ3txfbYwxESnmLokPHj3O8To/g1M6eV2KMcZEhJgLgu2FFQAMTrH+AWOMgZgMgmMADEy2OwJjjIEYDIIdRRWkdG5Ll3atvS7FGGMiQswFwZ7DlQxIsmYhY4ypF3NBcKD0OOnd2ntdhjHGRIyYCoLqOj+Hj9WQ3t0WojHGmHoxFQR5R44D0NvuCIwx5p9iKggOHKkCsDsCY4xpIKaCwO4IjDHmf4uxIKiiTas4kjraYjTGGFMvpoKgf2IHrsvsRVyceF2KMcZEjFBnH21Rbr4gg5svyPC6DGOMiSgxdUdgjDHmf7MgMMaYGGdBYIwxMc6CwBhjYpwFgTHGxDgLAmOMiXEWBMYYE+MsCIwxJsaJqnpdwzkTkWJg33kenggcDmM54WJ1nRur69xYXecmUuuC0Grro6pJJ29skUEQChHJVtUsr+s4mdV1bqyuc2N1nZtIrQuapjZrGjLGmBhnQWCMMTEuFoPgBa8LOAWr69xYXefG6jo3kVoXNEFtMddHYIwx5v8Xi3cExhhjGrAgMMaYGBfTQSAi3xYRFZFEr2sBEJGnRGSjiKwXkaUikuZ1TQAi8oyI5AZre0NEunpdE4CI3CQim0UkICKeD/UTkZkisk1EdorId7yuB0BEfi8iRSKS43UtDYlIuoi8JyJbg/8NH/a6JgARSRCRL0RkQ7Cu//K6poZEJF5E1onI2+H83JgNAhFJBy4H9ntdSwPPqOpoVc0E3ga+53VBQcuAkao6GtgOPOZxPfVygOuBD70uRETigV8Cs4DhwK0iMtzbqgB4CZjpdRGN8AHfUtVhwCTgwQj591UDXKaqY4BMYKaITPK4poYeBraG+0NjNgiAZ4FHgIjpLVfV8gYvOxAhtanqUlX1BV+uAnp7WU89Vd2qqtu8riNoArBTVXerai3wCjDX45pQ1Q+BUq/rOJmqFqjq2uDvFbiTWy9vqwJ1jgVftg7+RMTfQxHpDVwFvBjuz47JIBCRa4CDqrrB61pOJiI/FJEDwO1Ezh1BQ/cAi70uIgL1Ag40eJ1HBJzYWgIR6QuMBT73thIn2PyyHigClqlqRNQF/Ax38RoI9wdH7eL1IrIcSG3krXnAd4Ermrci53R1qeoCVZ0HzBORx4CHgCcioa7gPvNwt/Tzm6Oms60rQkgj2yLiSjKSiUhH4B/AN0+6I/aMqvqBzGBf2BsiMlJVPe1jEZE5QJGqrhGRS8P9+VEbBKo6o7HtIjIK6AdsEBFwzRxrRWSCqh7yqq5GvAwsopmC4Ex1ichdwBxgujbjwyfn8O/La3lAeoPXvYF8j2ppEUSkNS4E5qvq617XczJVPSoi7+P6WLzubJ8CXCMis4EEoLOI/EVVvxyOD4+5piFV3aSqyaraV1X74v4Cj2uOEDgTERnU4OU1QK5XtTQkIjOBR4FrVLXK63oi1GpgkIj0E5E2wC3AQo9riljirsJ+B2xV1Z96XU89EUmqHxUnIu2AGUTA30NVfUxVewfPWbcAK8MVAhCDQRDhnhaRHBHZiGu6ioghdcBzQCdgWXBo66+9LghARK4TkTxgMrBIRJZ4VUuwM/0hYAmu4/NVVd3sVT31ROSvwGfAEBHJE5Gvel1T0BTgDuCy4P9T64NXu17rCbwX/Du4GtdHENahmpHIppgwxpgYZ3cExhgT4ywIjDEmxlkQGGNMjLMgMMaYGGdBYIwxMc6CwBhjYpwFgTHGxLj/B4bfAQmsNKbvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(model.bgrid, v.T);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The value function converges quickly, and we can plot the conditional value function based on labor income draw above. \n",
    "\n",
    "We have assumed that parallel computation is more efficient, but this is not always the case, as parallel computing usually requires some overhead computation time and memory, as well as good programming decisions. To check if indeed parallel computation is desirable, we can create a $T$ operator that does not use parallelization, and compare computation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tnp = T_operator(model, use_parallel=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "74.7 ms ± 1.39 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "solve_value(model, Tp, 1.0, skip_policy=True, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "78.9 ms ± 314 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "solve_value(model, Tnp, 1.0, skip_policy=True, verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see above, parallel computation is on average a few microseconds faster, and I believe in more complex models the advantages are more significant. This result varies depending on your hardware and other processes going on in your system.\n",
    "\n",
    "We have now developed a sub-routine that is able to quickly solve the value and policy functions, given a price of bonds. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The stationary distribution and the $G$ operator\n",
    "\n",
    "\\begin{equation}\n",
    "    (G\\mu)(B) = \\int_S P(s,B)d\\mu \\text{ for } B \\in \\beta_S\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def G_operator(model, parallel_flag=True):\n",
    "\n",
    "    π = model.π\n",
    "    Nz, Nμ = model.Nz, model.Nμ\n",
    "    bgrid, μgrid = model.bgrid, model.μgrid\n",
    "\n",
    "    μ, g = model.μ, model.g\n",
    "    ind, w = model.ind, model.w\n",
    "\n",
    "    @njit\n",
    "    def interpgrid(g):\n",
    "\n",
    "        grid = np.empty((Nz, Nμ))\n",
    "\n",
    "        for zi in range(Nz):\n",
    "            grid[zi] = interp(bgrid, g[zi], μgrid)\n",
    "\n",
    "        return grid\n",
    "\n",
    "    @njit\n",
    "    def weight(x, ind):\n",
    "\n",
    "        next = np.take(μgrid, ind)\n",
    "        prev = np.take(μgrid, ind - 1)\n",
    "\n",
    "        return (x - prev) / (next - prev)\n",
    "\n",
    "    @njit\n",
    "    def get_weights(g):\n",
    "\n",
    "        gμ = interpgrid(g)\n",
    "        ind = np.searchsorted(μgrid, gμ)\n",
    "        w = weight(gμ, ind)\n",
    "\n",
    "        return ind, w\n",
    "\n",
    "    @njit\n",
    "    def G(μ, g, ind, w):\n",
    "\n",
    "        μ_new = np.zeros_like(μ)\n",
    "\n",
    "        if np.isnan(w[0, 0]):\n",
    "            ind, w = get_weights(g)\n",
    "\n",
    "        for zi in range(Nz):\n",
    "            for bi in range(Nμ):\n",
    "                if μ[zi, bi] > 0.0:\n",
    "\n",
    "                    bj = ind[zi, bi]\n",
    "                    wb = w[zi, bi]\n",
    "\n",
    "                    for zj in range(Nz):\n",
    "                        μ_new[zj, bj] += wb * π[zi, zj] * μ[zi, bi]\n",
    "                        μ_new[zj, bj - 1] += (1 - wb) * π[zi, zj] * μ[zi, bi]\n",
    "\n",
    "        return μ_new, ind, w\n",
    "\n",
    "    G(μ, g, ind, w)\n",
    "\n",
    "    return G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- $\\texttt{use\\_existing}$:\n",
    "- $\\texttt{tol}$:\n",
    "- $\\texttt{max\\_iter}$:\n",
    "- $\\texttt{verbose}$:\n",
    "- $\\texttt{print\\_skip}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_distribution(model,\n",
    "                       G,\n",
    "                       use_existing=False,\n",
    "                       tol=1e-5,\n",
    "                       max_iter=2000,\n",
    "                       verbose=True,\n",
    "                       print_skip=25):\n",
    "\n",
    "    start = time()\n",
    "\n",
    "    if not use_existing:\n",
    "        model.init_distribution()\n",
    "\n",
    "    μ, g = model.μ, model.g\n",
    "    ind, w = model.ind, model.w\n",
    "\n",
    "    i = 0\n",
    "    error = tol + 1\n",
    "\n",
    "    while i < max_iter and error > tol:\n",
    "        μ_new, ind, w = G(μ, g, ind, w)\n",
    "        error = np.max(np.abs(μ_new - μ))\n",
    "        μ = μ_new / np.sum(μ_new)\n",
    "        i += 1\n",
    "\n",
    "        if verbose and i % print_skip == 0:\n",
    "            print(f\"Distribution error at iteration {i} is {error:.4e}.\")\n",
    "\n",
    "    end = time()\n",
    "\n",
    "    runtime = end - start\n",
    "\n",
    "    if i == max_iter:\n",
    "        print(\"Failed to converge!\")\n",
    "\n",
    "    if verbose and i < max_iter:\n",
    "        print(f\"Distribution error at iteration {i} is {error:.4e}.\\n\")\n",
    "        print(\n",
    "            f\"Distribution converged in {i} iterations and {runtime:.4f} seconds.\\n\"\n",
    "        )\n",
    "\n",
    "    model.set_distribution(μ, ind, w)\n",
    "    return μ\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = G_operator(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distribution error at iteration 25 is 2.0038e-02.\n",
      "Distribution error at iteration 50 is 1.3703e-03.\n",
      "Distribution error at iteration 75 is 1.0302e-04.\n",
      "Distribution error at iteration 91 is 9.8976e-06.\n",
      "\n",
      "Distribution converged in 91 iterations and 0.0070 seconds.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "μ = solve_distribution(model, G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD4CAYAAADo30HgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3xW5d348c83O5BFBiGQQMLeokYUB1oVxQXWasVaa1ue2qEdT+3QDmsdj9ra6s8q9rHVVlvno7ZSt4iKiiIRUdmEHWYgg0B2cv3+uE64R+4kd0KSc4/v+/Xidc65znVOvrdCvvc51xJjDEoppVSbGLcDUEopFVo0MSillPKhiUEppZQPTQxKKaV8aGJQSinlI87tAHpDdna2KSwsdDsMpZQKKx9//PF+Y0yOf3lEJIbCwkJKSkrcDkMppcKKiGwLVK6vkpRSSvkIKjGIyGwRWS8ipSJyQ4DziSLytHN+mYgUep270SlfLyLnOmUFIvKWiKwVkdUi8kOv+jeLyE4RWen8Of/oP6ZSSqlgdfkqSURigQeAWUAZsFxEFhpj1nhVmw9UGmNGi8g84C7gchGZCMwDJgFDgUUiMhZoBq43xqwQkVTgYxF5w+ue9xhj7u6tD6mUUip4wTwxTAdKjTGbjTGNwFPAXL86c4FHnf1ngbNERJzyp4wxDcaYLUApMN0Ys9sYswLAGFMDrAWGHf3HUUopdbSCSQzDgB1ex2W0/yV+pI4xphmoBrKCudZ57XQssMyr+DoR+UxEHhGRQYGCEpFrRKRERErKy8uD+BhKKaWCEUxikABl/jPvdVSn02tFJAV4DviRMeagU/wgMAqYBuwG/hAoKGPMQ8aYYmNMcU5Ou95WSimleiiYxFAGFHgd5wO7OqojInFAOlDR2bUiEo9NCo8bY55vq2CM2WuMaTHGtAJ/wb7KUkop1U+CSQzLgTEiUiQiCdjG5IV+dRYCVzv7lwKLjZ3PeyEwz+m1VASMAT5y2h8eBtYaY/7ofSMRyfM6/CKwqrsfSql+YQx88jg01bkdiVK9qsteScaYZhG5DngNiAUeMcasFpFbgBJjzELsL/l/iEgp9klhnnPtahF5BliD7Yl0rTGmRUROBa4CPheRlc6P+oUx5mXgdyIyDfvKaSvw7V78vEr1nvWvwAvfg/0bYNZvfc+1tkJLA8QnuxObUkchqJHPzi/sl/3KbvLarwcu6+Da24Hb/creI3D7A8aYq4KJSSnXVW2328ZDsOdz+GABzPkTxMbBm7+F9++FX+6F+CR341Sqm3Tks1I9VV9tt0kZ8MzV8OkTULkFFv7AJgWAhoMdX69UiNLEoFRP1VfZbVI6tDbZ/bpKWPGop87i22DH8v6PTamjoIlBqZ6qcxKDxEBthd0/tNe3zopH4eGz4bGL+zc2pY6CJgalumPVc/DYXKg/aF8dAax70bYzACz5feDrNr8Fj3/Z9mSq2AJL7++feJXqgYiYdlupPvXCtfDJP+Hmanj2m7bsXa+pvLZ/4Nnf/WnH99n4Giy6GVY9D9XbYdpXYEBmn4Ss1NHQJwaluvLJP+22dJGnbPuHvnVyJ3v2C0/z7BfN9K33/r02KQA01fZejEr1Ik0MSgXrnd959ncsg+EnQ5wzTmHoNM+54Sd59gdk223GiPb32/VJ78eoVC/QxKBUZza87tlvbfY9lzMOmp1Rz0OmesoLvBKDabXb037c/t5Pf7V3YlSql2liUKozT3iN29z5MYw+23OcM96z750Yhh3nbI/3JIaEFM/5r73g2b+ryI6gViqEaOOzUh0x/pMIY58G2toacsZ5yodMhvEXQlyibVD+znuQNRrKSmDtQhhxMkz+km2LaG3xXFdXAZvfgXHn9e1nUaobNDGo6FOzB0rfhGOvDHx+4xu20bhtZLO3Ycd69nPGw5m/tqObE1Nh3uOec0Om2G3RabY3E8Clj9htcyMcdzUccwX8bTYsexDO/BUkej1VKOUifZWkos/z37KT31Vua3+ufAM8fincNhjWOK98ir/pOZ8zwbOfOgRm/gRm3dK9nx+XAHPugxEzPGV3DIOSv9n5l25Ot11alXKJJgYVfdpGLFduaX+uxmupkZd/Yrd5x3jK0obC9z6Eb74GEnAeyO750eeQOdLur/6XZxzEqueO/t5K9ZAmBhV9EtPstnqn3b7+a1jkTJtds8e3blyy77gEERg8wbdL6tHIGA4/+ASOvQp2rvBMqREb3zv3V6oHNDGo6GOcxt+2mU+X3gfv/RH2roGDfosTDj3WawxCLzwhdGT6t6CxBl663h7HJvTdz1KqC5oYVORrrIUFJ3tGK7c4M6E21PjWe3AGfPqkfaI44Vu2LGO4XV/hovvgu0v7Lsa8Y2CaV2N4jD4xKPdoYlCRb+8q2LcaXr3RHrclhPrq9l1S92+A1DwYkGWPBzojl4+/GnIn9m2cFy+AbzoD6uKT7RTeLc2dX6NUH9DEoCJf2yCzNnWVdttwsP1TA0BGgac8JbdvY/M3/ERIy7cN43cVwlu39e/PVwpNDCoaNDhTYovYJ4Q6Z+2Ehhqo3d++fsFJMON7dsDa8Vf3X5xtUod4BtGVvtn/P19FPU0MKvK1rbSG2AV12uY8qj/YvhcS2PEF6fl2wFpSer+FecTc+2H6NXbfe6oNpfqJJgYV+doSw6G9noV0ElLsq6Sa3fb4zF976g/u47aErgyeAOf/3g6mO5LUlOo/mhhU5Gub2qJ6h51+YmAOjDzDjnzetRIk1vMNHSB5kBtRtpecAVvfgzdvCTw9h1J9RBODinx1ft+6M0ZAXBIc3mfHMAw7HpLSPOd7Y0Rzb8idbJ8Y3v0DfPq029GoKKKT6KnI5/9tW2J8eyMVnmq3P1jpGfQWCmbfYZ9snr4ycCO5Un1EnxhU5KuvguyxMOM6e2xaYdZvPecHFdptZpHvvEhui42HCRfaQXYVm21ZS7Onl5VSfUQTg4p8h8rtgLWJc+1x7iTbwJs12h7391iF7socaQfetTTDv75tZ2JVqg9pYlCRa/uH8NYddtRz1mjIPwHmPQFnOCOgpzirs6UOcS/GYOROtrOu3poFq561Zd6L/SjVy7SNQUWmDxbAazd6jvNPsI3K4y/wlJ3+cxg7G4ZO6//4uuP0n8O+NbBpsaesqdYuDqRUH9AnBhV5Du/3TQrgO0FdG5HQTwpge0xd/rhvWWOtO7GoqKCJQUWe/Rt9j2Pi7Ayp4SxhAKQO9Rw3aWJQfUcTg4o8FZt8j1sjZIbSweM9+5oYVB/SxKAiz4FN9inh7JvdjqR3ne3VxXbnCvtHqT4QVGIQkdkisl5ESkXkhgDnE0Xkaef8MhEp9Dp3o1O+XkTOdcoKROQtEVkrIqtF5Ide9TNF5A0R2ehsQ2R+AhXyDmyC9+6FDa/ZsQmTLnE7ot6VNxW+/pLdX3gd/OUL7sajIlaXiUFEYoEHgPOAicAVIuI/y9h8oNIYMxq4B7jLuXYiMA+YBMwGFjj3awauN8ZMAE4CrvW65w3Am8aYMcCbzrFSXXvoDFj0G9s9NXOUZ7GdSJKa53YEKgoE88QwHSg1xmw2xjQCTwFz/erMBR519p8FzhIRccqfMsY0GGO2AKXAdGPMbmPMCgBjTA2wFhgW4F6PAhf37KOpqOM9ncXAbEgY6F4sfeXI+tNK9Z1gEsMwYIfXcRmeX+Lt6hhjmoFqICuYa53XTscCy5yiXGPMbudeu4HBgYISkWtEpERESsrLy4P4GCpiNdXDY37fVQpPs91Rz70D5i9yJ66+EBsH+dPdjkJFuGASQ6CpJk2QdTq9VkRSgOeAHxljujV7mTHmIWNMsTGmOCcnpzuXqkiz+W37x9sx8+x2xveg4IT+jqhvffU5T3JoaXI3FhWRgkkMZUCB13E+sKujOiISB6QDFZ1dKyLx2KTwuDHmea86e0Ukz6mTB+wL9sOoKLV/g+9xUnroTJ3dF5LSYLLTsL7hNZ0eQ/W6YBLDcmCMiBSJSAK2MXmhX52FQNviuJcCi40xximf5/RaKgLGAB857Q8PA2uNMX/s5F5XAy9090OpKFO1zffY/3k2ErU1Qj99Jbz//9yNRUWcLhOD02ZwHfAatpH4GWPMahG5RUTmONUeBrJEpBT4MU5PImPMauAZYA3wKnCtMaYFOAW4CjhTRFY6f8537nUnMEtENgKznGOlOlaxGfKmwcyfOQVRkBm8p/IoW+5eHCoiBTVPgDHmZeBlv7KbvPbrgcs6uPZ24Ha/svcI3P6AMeYAcFYwcako99JP7GujA6X2nfvI02HJ78BEQWIYVGiXI/3oIbcjUREozCeQUVFt+V+cHYFjvgLxA+xhS6NrIfWr838P5evspIFK9SKdEkNFAGPXW2hbeKelwd1w+tOAbKg94HYUKsJoYlDhp2IzvOI3IH5Qoe2tM+osuPAeV8JyRWqenTTw0TnadVX1Gk0MKvy8dD0se9C3LN0ZN3nV81D8zf6PyS0jz7DbLe/Ans/cjERFEE0MKvzUVfoeS2zor9vcV8bMgpk/tftVOzqvq1SQNDGo8FOxxfc4dQjExLoTi9tEYMZ1dr9aE4PqHZoYVHiprYD6Kt+ytKGB60aL5AxITIeq7W5HoiKEJgYVXg44q7OdfgNkj7X70Z4YwDa8f/QQLPw+lDzidjQqzGliUOGlYrPdTv4SjD3X7qfluxdPqBh5ut2ueAxe/G93Y1FhTxODCi8Vm0BiYNAIO9026BMDwEX3QcFJnuNoGP2t+owmBhVeDmyC9HyIS4QDG21Zuv/yIFEoJtbz1ADQVOteLCrsaWJQ4eGJeXD7UPsqKXOULTv2KrsdpVNrAXaQX5v6atfCUOFPE4MKbetfgf/Jhw2vQNNh+yopy0kMUy6Fm6ttw6uCQUWe/bqqjusp1QVNDCq0bXsfGms8x/XVkDnSvXhCWX6xZ//VG+Cg/3paSgVHE4MKbYcCrOedrr2QAoqNh++vsPtb3oHnr3E3HhW2NDGo0HZob/uygYP7P45w4f06qaGm43pKdUITgwpt1WWexuY2KZoYOhTj9U86NsG9OFRY08SgQldrC1RuhQkXwddf8pRrYujc2Nl2q43yqoc0MajQVV0GrU22F1LCQE95Qop7MYWDuQsgdaiuz6B6TBODCl0VzrxImSMhIdVTLgGXC1dtBmZBwXQ4uNPtSFSY0sSgQlfbvEiZIyFRnxK6JT0fDpTCYxe7HYkKQ5oYVOjY/Sncnge7VtrjA5shLtkuX5noPDHoGIbgZDq9kza/ZdtqlOoGTQwqdLz7RzvHz0POnD8Vm20iELFtDF/+B8x/w90Yw8XwGZ59/xXvlOqCJgYVOhK92hGaG2xiyPJ6Qpg4BwZm939c4Sh3Epx/t90/vN/dWFTY0cSgQkdthWf/tsGwf73vgC3VPTnj7LZmt7txqLCjiUGFDv8lO0HXWjgabe0x/7gYGg+7G4sKK5oYVOioq4Tscb5lA3PciSUSpHol1R3L3ItDhR1NDCp01FXBsON8yzQx9FxMDMx7wu7v3+huLCqsaGJQoaG1BWr32+kuTr/BU67TXxydcedDXBJU73A7EhVGNDGo0FC5FVoaIXssnHa9p1yfGI6OiG2n0bUZVDdoYlChYd9au82ZYNcVaJOc6U48kWRANqx6Dqp1igwVnKASg4jMFpH1IlIqIjcEOJ8oIk8755eJSKHXuRud8vUicq5X+SMisk9EVvnd62YR2SkiK50/5/f846mQZoxnv7wtMYzznQspRr+7HLWRzoDBpfe5G4cKG13+qxORWOAB4DxgInCFiEz0qzYfqDTGjAbuAe5yrp0IzAMmAbOBBc79AP7ulAVyjzFmmvPn5e59JBUWNr8Nv82AF//bHu9bB+nDfedEGn6yK6FFnJk/te0M5evcjkSFibgg6kwHSo0xmwFE5ClgLrDGq85c4GZn/1ngfhERp/wpY0wDsEVESp37fWCMWeL9ZKGizGNz7bbkETsf0r41MHi85/wvdutCM70lLhHGngt7V7sdiQoTwTynDwO8uzSUOWUB6xhjmoFqICvIawO5TkQ+c143DQpUQUSuEZESESkpLw+wLrAKHx8+YBNDjldiSBgAscF8b1FBSS+Aqh06oZ4KSjCJIdDk9ybIOsFc6+9BYBQwDdgN/CFQJWPMQ8aYYmNMcU6O9lwJO/ED2pcNntD/cUSLIVOgpQE++ovbkagwEExiKAMKvI7zAf++b0fqiEgckA5UBHmtD2PMXmNMizGmFfgL9tWTiiQtzXYW1fiBvuUpue7EEw1GnWm3nz7hbhwqLASTGJYDY0SkSEQSsI3JC/3qLASudvYvBRYbY4xTPs/ptVQEjAE+6uyHiUie1+EXgVUd1VVhqrHGbnP9+jAkZfR/LNEiZTCc9hO75sWfjte5k1SnukwMTpvBdcBrwFrgGWPMahG5RUTmONUeBrKcxuUfAzc4164GnsE2VL8KXGuMaQEQkSeBD4BxIlImIvOde/1ORD4Xkc+ALwD/3UufVYWKBicx+C+6k5Te/7FEk+yxdnugFMrXuxuLCmlBte45XUZf9iu7yWu/Hrisg2tvB24PUH5FB/WvCiYmFcaqy+w2x2/CvGR9YuhTOWM9+7putuqEjh5S/e9v59lt7hTf8sS0/o8lmmSN8ew3HHIvDhXyNDEo9wzxSwxxOm6hTyWmwEnX2v1GTQyqY5oYVN+qq4LFt8H+Uk/ZgGz77TUtD079sXuxRaPjv263+sSgOqGJQfWtt++AJb+H+4+Hmj22rLkexpxj9091+haMu8Cd+KJN25Qjz/+XjoRWHdLEoPpWU51nv7bCHjcegmRnQHtSGly/AS7/pzvxRRvvacxL/uZeHCqkaWJQfSsu0bNfXwWbFtv9odM85am5Ootqf4mNh2HFdj+joPO6Kmrpv0bVt5obPPt1VZ4ZPofPcCceBV97wW69n+aU8qKJQfWt+mrPLKl1FVC1HQZk+U6vrfpXYoqdjqT+oNuRqBCliUH1rcZDdsStxNikcGgfpOZ1fZ3qW0npULHZ7ShUiNLEoPrOiz+G0kX2CWFQIZQth8PlOpAtFLQ2w4ZXYPW/3I5EhSBNDKpvNNVBycN2PzEVUofahuey5bYnknLXHGeZz81vuxqGCk2aGFTf8H5N0dIEdZWeY50sz33jzrMLI9VWuB2JCkGaGFTfaBvMBnbtBe/EoK+SQsOALKg94HYUKgRpYlC9wxj4YAEcKoc1L8Cbv7XlKUPgwnt85//XWVRDw8Bs2Pa+bxJXCk0MqrfsWgGv3Qgv/gie+ZpdEAbg+yWQPQaOvdJTV1dqCw0T59rt8ofdjUOFHE0MqnccKrfb+mqvQvEs33nObZA+3O6nDO7X0FQHJn/Jdgo4UNp1XRVVNDGo3lHjLOW99V1PWcJAz1QXMbGQf7zdH1TYr6GpTuRO1MSg2glqBTelutS2XKe3hIG+x2f9BiZdAnnH9E9MqmuZo2DbUttGpKu6KYc+Maje0Vjbvsw/MWQWwcQ57esp9+RNtb3GVj3ndiQqhGhiUL0j0IpgMfH9H4fqnqnzbM+xtf9xOxIVQjQxqKO3/lX44P725TH6pjLkxcbB0GNh/wa3I1EhRBODOnpPXh64XEc4h4esUXakekuT25GoEKGJQR2d1lbPftZo+NZbnoVgYvVVUlgoONEut1pW4nYkKkRoYlBHZ5/XusFN9TDsOJj3uD0u/qY7ManuKZppp0Xf/JbbkagQoYlBdU/lVlj5pOfYezqF1ma7TR0Cv6mCSRf3a2iqh5IzYOhxsEkTg7I0MajueXQO/Ps7nmUhvUc6J6Z69rVPfHgpmmmnNQnU7VhFHU0MqnsOO1Nf/Os79tVRg7M85Izr4IonO75OhbbhM+wTn/fIdRW1NDGo7olPtts1/4aV//SsG/yFX9jJ8lR4GnmG7UWm4xkUmhhUd8Ule/bbnhgkFuIHuBeTOnpxCTDiFDsNt4p6OgJJdU+M13eJ139pt8mDtE0hEhSeCutfhoO7IG2o29EoF+kTg+oeiW1fpiuyRYYRp9jt1vfcjUO5LqjEICKzRWS9iJSKyA0BzieKyNPO+WUiUuh17kanfL2InOtV/oiI7BORVX73yhSRN0Rko7Md1POPp3qdBPgrk6SJISIMmQJJGbD5HbcjUS7rMjGISCzwAHAeMBG4QkQm+lWbD1QaY0YD9wB3OddOBOYBk4DZwALnfgB/d8r83QC8aYwZA7zpHKtQ0VzfvixRp76ICDGxttvqyn96Fl5SUSmYJ4bpQKkxZrMxphF4CpjrV2cu8Kiz/yxwloiIU/6UMabBGLMFKHXuhzFmCVAR4Od53+tRQEdJhZKGALOo6hND5Ciaabd/GOtuHMpVwSSGYcAOr+MypyxgHWNMM1ANZAV5rb9cY8xu5167gYDrQIrINSJSIiIl5eX67abPrXwS/lQMDdUw/RoYcarn3IBM9+JSvWus8xBvWjuvpyJaMIkhUHcTE2SdYK7tEWPMQ8aYYmNMcU5OTm/cUnlrafY9fvUGOLDR7qfmwddegLkP2OOs0f0bm+o7GQUw7nzInex2JMpFwSSGMqDA6zgf2NVRHRGJA9Kxr4mCudbfXhHJc+6VB+wLIkbVm7Yvg1uzYNsHnjLv1dgSU+08/tOuhC//A078Tv/HqPrOgEzYu0qnx4hiwSSG5cAYESkSkQRsY/JCvzoLgaud/UuBxcYY45TPc3otFQFjgI+6+Hne97oaeCGIGFVvWv5Xu93r1WHMe9GdticEEbtUZ7zXoDcV/rLH2e3iW92NQ7mmy8TgtBlcB7wGrAWeMcasFpFbRKRtAd+HgSwRKQV+jNOTyBizGngGWAO8ClxrjGkBEJEngQ+AcSJSJiLznXvdCcwSkY3ALOdY9adWZ8GWuCRPWV0VZIywTwlt/d1VZJr+LbstX+9uHMo1QY18Nsa8DLzsV3aT1349cFkH194O3B6g/IoO6h8AzgomLtVH2tZqbnJeJRhj13Q+YT6c/Rv34lL9Iz7ZJv/mBrcjUS7RKTFUe22vjRoPwb618PmzYFogQedDihopubBtqZ0PKz6p6/oqouiUGKpjjYdhwUnw7t32OH5g5/VV5Dj2q3BoD6x4zO1IlAs0MUSrZf8LdxUGPtfivELw75WiTwzRY/RZdnLE8nVuR6JcoK+SotUrP7PbliaIjfc91/ZuedmDvuX6xBBd0gugYrPbUSgX6BNDtGs83L6so0ZHfWKILsNPgi3vwMHdbkei+pkmhmjX5Pe6qLbC8yrJny7GE12O+5qdGmOLzrYabTQxRKsjPY+8EsPHf4ffFcHeNfZ49Czfa1Jy+yU0FSIGT7RjWfZ87nYkqp9pYohWbYmhyetVUskjdlu7H8acA199Fn5T5Tmfnt9/8Sn3xcTaOZM2v23HsqiooYkhWrUti9F4GOoq7chm74nzkp31kbyX7NTptaPPlMvs1CiVW92ORPUj7ZUUrbxfJbV1W82f7jmfO8mzf/7dkNbVbOkqIo08w263vguZRW5GovqRPjFEqxjnicH7VVKr1xPDuAs8+9O/BePP75+4VGjJGQcDB8OWJW5HovqRJoZoFajxuWY3TLoEfrYFsnWNBYV9lVh0Gmx5V9sZoogmhmjlPR9Sm5rdth1BV2RT3opm2ukx9m90OxLVTzQxRKrGw7Cjk6Uv2hLDyz/xLU/K6LuYVHhqWwdaxzNEDU0Mkerf34OHZ8GhDhbAi+ngf33qkL6LSYWnQUWQPlwTQxTRxBCpdn9qt/UHA59vaQpcrolB+ROxTw1b3oXWVrejUf1AE0Okik2w25bGwOc7mg8pdWjfxKPC28jTob4K9nzqdiSqH2hiiFRtiaG5LvD55gb7LfD6DfDtdz3lmSP7PjYVfkaeYbelb7oZheonmhgiVWyA7qhg50N64yY7Ud7Q4yA1F4ZM8ZwfmN1vIaowkjLY/j3Z9Jbbkah+oCOfI1Xbus3+02r/54ee/bhEuxWBaV+FzELfKTCU8jbqLPjgfmiogcRUt6NRfUifGCJV26sk73EK/g2HbYkB4OIHYOZP+z4uFb5Gn2VHx299z+1IVB/TxBCpjrxKOmxHrL53D5Sv9auT2P46pTpScKJdk0PbGSKevkqKVG1PDE21UL4eFt1s13n2FqeJQXVDXCIUngabNDFEOn1iiFRH2hgOeXom1fgt0aiJQXXX6LPsOtAVW9yORPUhTQyRbvFtHS/oPuqs/o1Fhb9RZ9rtpsXuxqH6lCaGSOU9sO3f1/qem7sAfvAJpOsaC6qbskbb6TE0MUQ0TQzhasndUPZxx+dbvEY2ew9yi0uCaV/RgWyqZ0Rg1Bdg3YtwYJPb0ag+ookhHBkDi2+Fv57ZcZ2OprzIGK5jFdTRmTjXbj9c4G4cqs9oYghHHU2A562pDlJy7bKc3jJG9E1MKnqMPgvSC3R9hgimiSEcNdd3XaepDkacbJflvOJpz7QX+gpJ9Yb8YjsNd0fTuquwpokhHHWWGF7/FXz2jK0Tl2zLxs2GQmexFe95kZTqqWlX2u2+Ne7GofpEUIlBRGaLyHoRKRWRGwKcTxSRp53zy0Sk0OvcjU75ehE5t6t7isjfRWSLiKx0/kw7uo8YgZo6mDG1rgqW/gme/5Yd2Baf7Dl3wnyY+TOYcln/xKgiW940kFhY/W+3I1F9oMuRzyISCzwAzALKgOUistAY4/1VYT5QaYwZLSLzgLuAy0VkIjAPmAQMBRaJyFjnms7u+VNjzLO98Pkik3/Dcn01fPgg5Iz3lDXV+yaGrFFw5i/7Jz4V+VJyYMw5Om9ShArmiWE6UGqM2WyMaQSeAub61ZkLPOrsPwucJSLilD9ljGkwxmwBSp37BXNP1RH/NRbe+R28fQe8/mtPWdNh2zVVqb6SXwwHSnUUdAQKJjEMA3Z4HZc5ZQHrGGOagWogq5Nru7rn7SLymYjcIyI6b4M//yeGtjaH6u2+5d5PDEr1timXAQY2vu52JKqXBZMYAnV6N0HW6W45wI3AeOAEIBP4ecCgRK4RkRIRKSkvLw9UpUsV+3ayc/Pariv2p8pt8Oatna+t693GULMXypZ7nfT6Tzvs+F4PT6kjBo2AtHxdvCcCBZMYyoACr+N8YFdHdUQkDkgHKjq5tsN7GmN2G6sB+Bv2tVM7xpiHjDHFxgNjs+gAABjBSURBVJjinJycID5Gexuf/gXJj53To2v7zAvXwrt3w+6VHdfx7pX0h7Gw22sd3vEX2G16gR2hqlRfmjjXzrbacKjruipsBJMYlgNjRKRIRBKwjckL/eosBK529i8FFhtjjFM+z+m1VASMAT7q7J4ikudsBbgYWHU0H7BTEoO0e/hxWUys3dZVdFyns+6qOePg6y/Df+nUyKofjD3Xzsu19d2u66qw0WWvJGNMs4hcB7wGxAKPGGNWi8gtQIkxZiHwMPAPESnFPinMc65dLSLPAGuAZuBaY0wLQKB7Oj/ycRHJwb4TWQl8p/c+rh+JIYZOXtm4IXmQ3dZVdVynqZPEkDEcCk/p3ZiU6sjwGZCQYtsZxp3ndjSqlwS1UI8x5mXgZb+ym7z264GAHeSNMbcDtwdzT6e8kwmAepc5mieG+oO2cTc2vneDSkyz24aDHddp65V0/Ndh3AXwhNd/+rxjejcepToTlwAjz4ANr9s5vHQerogQ3SOfJYZY08MnhjsL4Jmru67XXW2vkjqaBM/73Jk3wdhz4NcH4NtL4Au/hCFTez8mpToz/gI4WAa7PnE7EtVLojwxxB5dG8P6l4KvW7oIGmqCigkI3I5QVwWfP+vplRTvjFOIjbNPCqf/zJNYlOovY2fbv7frXnQ7EtVLojoxGJH+aWOo2Qv//BI8+82u67Y9ijfVw+EDULnV/mlpgrtGwHPzPbNa6gA2FQoGZELhqbD2P25HonpJUG0MEUtiiOnJE0NnYwwCaXtdteOjruu2TandXA9/nOBZcOeSv3rqVG2zazrr04EKFRMugpd/AuXrbc84Fdai+omhx91Vbceq4LU22219ddd1vROD9yps29737Fdu1VHNKrSMv9Bu1/r3ZFfhKOoTQ49eJbX2MDEEk4Ta1mqu9+uV5D26uXpHcO0VSvWXtDzIP0FfJ0WIqE8MsWIw3X01dOQXfR/Ub3tK+PQJ3/K9q3y7omaN7l4MSvW1CRfZUfiV29yORB2lKE8M9h29HaTdDX2aGDpZtjNjhPPILjBfJy5TIabtddK6bvTWUyEpyhOD7QHU2u1XQ92sH8wazW38xy/c4DVjauZIuPyfcNMB2xNEqVCSNQpyJ+vrpAgQ1b2SxHliaA2lV0mNhz3733gVktI9x7mTbTIT7Y0USGurYXtFLfXNLeyuqictOZ7Rg1NYWrqfv763hTGDU7j42GGs3nWQ51eUcc/l0xibm0p5TQOxMULmwASf+9XUN3G4oYUh6dotOGjjL4R37rJrQacMdjsa1UNRnRhMT58YetoryVtLk21A9v/m33jIrst84ndhxAxbljXaLogS4XMg1TY2s+1ALRPy0tqdq2tsob6phb019fz13S0s23KAkdkpDE5NJD05npJtlZRV1rL/UGOH9/94WyVPLfcsA3Lpg0spLsxk+dYKauqbGT8klYS4GHZV1dHUYqius0968bHCmMGp3HrxJBLjYslJTSQ3zSYLYwyi00B4TLgI3rnTvk4q/obb0ageiurEIGLfpPV74/OWJfDCdXY8wq/KYWcJxCXa9RMaDsLwk+HYKz31/2sRVO2AtKHd+7khrLq2ifg4ITk+lh0Vdeytqefb//iYisON/Hz2eLIGJvDh5gOICK+s2k1toycZx8YILa2GHRWedSkmDU0jLz35SGK4+7JjWLv7II8v28ZXpo/gp+eO4+NtlTy3oowTCjOZkJfKT5/9jMXr9jFpaBq7qupYt8f29EpJjGNYRjLVdU0kxMWQlhTHmt0H+dKDHwA2Ucwck8N7pftpaG7l6hkjGDckjeyUBGZNzI3uRJE7CQYV2W6rmhjCVlQnBmJsYmhp6cNXQ9C+jeHRizz7u1fC35xZKW+qtE8Riam+9ZMHeWZdjQAtrYaz/vgO+w/5tqckxdv/H3e9us6nfNLQNIZnDmBT+SFSEuNYcOXx7D/UQOm+Q2QMiGdUTgoFmQMA+LysmsFpnm/0v75w4pH7nDomm1PHZB85fv1HM6mpbyZ9gJ0IcfWuamobWzhu+KAjySc2RqhrbOG/n17JwMQ4JuSl8ud3NvH+pv2kJMbR0NzIox94euEMy0hmQEIs+2oaGD8klcuKCzhv8hAGJMRGR8IQsWs0LP0THN4PA7O7vkaFnOhODD1uYziKV0n+SWjjG579t++wC54kpnTv/iGiqraRuNgYUhLjaGk17Ky0TwJ3vLyWFduruPOSKeSmJ3Hjc5+3SwoAi358OtkpifxfyQ5WbK8iPTmeOdOGckx+BrExvr9Uh6QnMXlYert7TMlvX9aRmBg5khQAJg31vbbtZyYnxPLnqzyr4X3zlCKaWw0JcTGs2llNWWUtk4elc/tLa48kq437DrFsSwXLtlTwk//7lKLsgXxh3GDOmZTL9MJMYmIiOElMuRTevxfWvAAnzHc7GtUDUZ4Yevoq6SgSw6E9vud2lnj2Sx6G1ibfBucwYIyh4nAjx9+2CLC/ON/dWM7Gfb6ret3w/OeA/VJ5wdQ87r18Gk0trSTHx9JqPL+Ir5pRyFUz+vczdEdMjJDgxDp5WPqRBPXgVz3J4631+4gRYcOeGm5/eS07K+t45P0tPPL+liN1Lj0+n2tmjmRsrt8TYrjLnQzZY2HVc5oYwpQmBsAczS/67tb/86m+5zYthoRUSB8G5c4rlIwR3bt/PyqrrOXb//iYu740lYl5aTy7ooyH393C+r2ekdjev/y+PXMklxXnIyL8+e1NNLa0ctOFE8lKSQQgPtb+P4iNsC/QXxhne+ScPjaHb80cCcDS0v0s2bif0n01LFq7j2c/LuPZj8sYMziFmWNzOGV0FjPH5FDf3EpKYhj/0xSBqV+GxbfB3tW23UGFlTD+23f0xGlj6P44hqNoY6irbH9+3HmQWQTvOInBpX9IOypq2VlVxx9eX8+vLpjIyh1VnDwqi9GDU7h/cSlPl+xg78F6mloMF/7pvXbXX3p8Pt86bSQrd1RScbiJM8cPZtwQz7fh318W3YsInTw6m5NH23fura2GlWVVLNtcwYuf7eLxZdt4+D1PQr16xghOH5fDKaOzSYwLw+7JxfPh/fvgvXvhS39xOxrVTVGdGNqeGPq1jSGQwRPgpO/Zx+8RJ/dL76Paxmb2HmygKHsgn5dVExsjXPfECjbvt+Mo5j7wfrtrUhPjaGrxjBIfkTWAaQUZ/PrCiWQkxxPnfPv3TgYqsJgY4bjhgzhu+CC+e8YomlpaeXPtXj7cXMHfl27l0Q+28egH2xiansScacOYc8xQxg9JDZ+2iQGZMPkS+OwZOzYnYaDbEalu0MRAL75KWvgDe8+L7g2ufptx59tFd6Zc2r04eqB03yEqDjfy6Adbeemz3Sy94Uwuur/9t/+0pDgO1vvGveKmWdQ1tZCW1MvLmSriY2OYPTmP2ZPz+M1FEymrrOPtDeW89Nku/vruZv78zibSk+P5wrgcUpPi+e2cSaGfJKZ8GT7+O6x/pV/+bqveE9WJQWLaeiX10gC3FY/a7UX3worHYNRZtu2gxW/Q1ZXPwsgv2PUU+qgLY9v8T5W1TWQOTKB0Xw2FWQM5+4/v+NQ7+c7FPsejB6dwxyVTOKHQDryrbWxmc/lhMgcmEB8bc6RNQPUdEaEgcwBXnTSCq04awb6D9by5bh9vrt3Lv1fuAuClz3czeVg6N543nrG5qe16bYWE4TMgbZh9atDEEFaiOzH01gC37ct8xyYc2gcLv2/357/h265w7v/AmFk9iLZ7Fn66ix8+tRKwfflvfXENJxYFnl/p3Em5FGYNJCslga+fXERCnOeX/4CEuIDdQlX/GZyWxBXTh3PF9OHUN7Xw9PId/O39LSzZUM6SDeWMH5LKORNz+cYpRQzym9bDVTExMPlL8OECuxrhwCy3I1JBiurEYHraxuD9BGAMvHmL76I6ez7z7D88C07+vl2G83sfwqDCngfsp6a+ic/LqmlqNSwt3c+XTyhgxbZKdlTWcd+bG4/Uu/XFNQAs21LBoAHxXHJcPgCXFeczMjvFJxGo0JYUH8vVJxdy9cmFbN1/mNfX7OGZkjLuW1zKfYtLSUuK44KpQ7n+nLFkDUxwf1Dd1C/D0vtgzb+162oYierE0OMnhibPVAxs/wCqtvue37XS93jfOhiYY3seHYWG5hbW76lh0tB0FrxVyh/e2OBz/n+XbPY5josRpuSnUzBoAFPz0xmbm8r0okyS4sOwl4tqpzB7INfMHMU1M0fxWVkVr67aw4K3N/HkR9t58qPtnDd5CJcV53Pm+Fz3gsydDDnj4dMnNTGEkehODDGBG58PVh2g7IE5pF3+v+SPntz+Qu/E8LfzIMbvP+PiW32PS9+A0WcHFdOW/YfZsLeGmWNyWLxuH1sPHOZwQzNnT8zlkgVLAThueAYrtlcduWb04BT2H2qgqraJvPQkbpk7GWMMJ47MIj1ZG4qjwdT8DKbmZ/Cjs8fyzoZy7nxlLa+s2sMrq/ZwTEEG3zylkNmTh/R/11cROO5qeO1Gu4hPXnR3WQ4XUZ0Y2nollS192icBbFjyDMVNq/j4hV+Tf/2/2l/XVOt7HKjXUVIGzPkTPHOVPb6k677cra2Gqx/5iO0Vte3OLXh705H9FdurmFaQwaPfmE5achwidl6f/YcajswRpKJTQlwMsybmMmtiLvsPNfDsx2U89dF2fvjUSlIT4zhldDZFOQP577PH9t8rxGlfsV+Wlj8Mc+7rn5+pjkpUJ4a2Xkknbb4P8HzLj4m3v1xjWjuYwrmpvuubX/AHOwXxsV+FIcfw7s4WkuIrOKEwk0+2V7Jx3yHW7a6htrGZz8qqKT/UQHlN+/mDvP303HG0ttreRtedOdrn/XFsjGhSUD6yUxL5zumjmH9qEUs3HeA/n+7i2Y/LAHjxs13816kjmTttKBkD+rjBOjnDdsle+x8459awm/IlGkV1YvBe8GbDircZe9wZAMTE2ekaYloC/6JuaTxMhw/kZ/6a2tHns+DzWL6UV0vR3AeoqW/iqpvtUpxzjhnKwk93dRrWoh/PZNXOg4zPSyUvLZnkhFjqm3X8gOqZ+NgYTh+bw+ljc7jjkik89sE2bn9pDb9ZuJrbX1rLqWOyOWdiLudNzvOZVLBXnfRdWP08/OeHMO4C233V7YZx1aGoTgzi1fd77MK5VOSvYcPztxOTMwaAWGOnsqiu3E9CQiLJA+2I3vrDhxgIzGv8FU9dfwkL3/2Ivy4r54ux77G/4jSW/ruST7ZXcf9bpZw6Opv3Svcf+TltSSEhNobfXzaVysONJMbHcuzwDGLFzvY5ODWJ0YN9Rw9rzyHVG+JjY5h/ahFXnjicDzcf4M21+1i0di+L1+3jV/9exYxRWZw+NofiwkyOyU/vvV5N+cUw8WKbHFb/y472HxKg/U6FhKhODP5LZJa++3+ctOdxqvbYaa/jW+ppbKgn4d4JNEk8e2KyaTj3bhIqKskzCXzYOpHC368D0oA0PmseBR/u97mnd1IoHjGIzIEJXHzsMM6fktfXn06pDiXFx3LGuMGcMW4wt8ydxKdl1by6ag+vrtrNbS+tBaAgM5k5xwzlnIlDmNobSWLaV2xiACj7SBNDCIvqxOD/F72lyi77mGYOg8CAloOsfPkvTJdGkmkkrfUwK969nwP1DSTR8ZoJJ4/K4oGvHEdsrNgngrhYslMSjswlpFQoERGmFWQwrSCDn88ex4a9h3hr/T7eWLOXBW9v4oG3NjEkLYmzJgzm7Am5zBiV1bMuz2NmwbeXwP/OhIrNXddXrpG2qRPCWXFxsSkpKem6op+Vi55k2nvf6dHPXCGTSPvua/xlyRa+fEI+x48IPKpYqXBWVdvIG2v28ubafSzZWE5tYwtJ8TEUj8jkpJGZnDNpCGMGp3TvaeL+EyBnHFz+z74LXAVFRD42xhT7l0f1E0NrkEt6roubwLhfLGXLmuVUvP47sg9tIGbGdYwenMpdl07t4yiVck/GgAQuKy7gsuIC6pta+HDzAd5at48lG/dz9+sbuPv1DWSnJHLamGxOHZ3NiSMzGZaR3HmiGFQIlVv76yOoHggqMYjIbOD/AbHAX40xd/qdTwQeA44HDgCXG2O2OuduBOYDLcAPjDGvdXZPESkCngIygRXAVcaYDvqNHp2mg3uP7JeknU3riFOJ2fYexQcXsUOGkjD/RSr3bGV88VkAjJx8IiMnP9cXoSgV8rzbJQB2VtXx3sZy3i89wJIN5fzrk50ADElL4vjCQRSPGMTxIwYxfkiab+eJQYWw/UM7nYz2TApJXb5KEpFYYAMwCygDlgNXGGPWeNX5HjDVGPMdEZkHfNEYc7mITASeBKYDQ4FFwFjnsoD3FJFngOeNMU+JyJ+BT40xD3YWY09fJdUdrmHdgnkMO7yGxB98RHpmDg31tVSW7yIrt4D4hMRu31OpaNTaali3p4aSbRWUbK2kZGsFu6rteJ+EuBgm5KUxLT+d0bmpTCh7huJVt8EPPoHMkS5HHt06epUUTGKYAdxsjDnXOb4RwBhzh1ed15w6H4hIHLAHyAFu8K7bVs+5rN09gTuBcmCIMabZ/2d3pKeJQSnVd3ZW1fHJ9ko+K6vm0x1VfL6zmtrGFgplN4sTfkKNpFAVk3HUPyf8W0mPTt25f2DCiZ3+iuzQ0bQxDAN2eB2XASd2VMf5hV4NZDnlH/pdO8zZD3TPLKDKGNMcoL7/B7oGuAZg+PDhQXwMpVR/GpaRzLCMZC6calckbGk17D1Yz+7qOv5v8WFGVbyNdLS2SZAk6tMCpCV13EOyp4JJDIFeAvr/3+ioTkflgfptdla/faExDwEPgX1iCFRHKRU6YmOEoRnJDM1I5vhv/BD4odshqQ4E07G+DCjwOs4H/Od0OFLHeZWUDlR0cm1H5fuBDOceHf0spZRSfSiYxLAcGCMiRSKSAMwDFvrVWQhc7exfCiw2tvFiITBPRBKd3kZjgI86uqdzzVvOPXDu+ULPP55SSqnu6vJVktNmcB3wGrZr6SPGmNUicgtQYoxZCDwM/ENESrFPCvOca1c7vYzWAM3AtcbYl4qB7un8yJ8DT4nIbcAnzr2VUkr1k6ge+ayUUtGso15JOnmPUkopH5oYlFJK+dDEoJRSyocmBqWUUj4iovFZRMqBbT28PBs7fiLUaFzdo3F1j8bVPaEaFxxdbCOMMTn+hRGRGI6GiJQEapV3m8bVPRpX92hc3ROqcUHfxKavkpRSSvnQxKCUUsqHJgZnIr4QpHF1j8bVPRpX94RqXNAHsUV9G4NSSilf+sSglFLKhyYGpZRSPjQxeBGRn4iIEZFst2MBEJFbReQzEVkpIq+LyFC3YwIQkd+LyDontn+JyNGvz9gLROQyEVktIq0i4nrXQhGZLSLrRaRURG5wOx4AEXlERPaJyCq3Y/EmIgUi8paIrHX+H4bEKj4ikiQiH4nIp05cv3U7Jm8iEisin4jIi715X00MDhEpAGYB292OxcvvjTFTjTHTgBeBm9wOyPEGMNkYMxXYANzYRf3+sgq4BFjidiAiEgs8AJwHTASuEJGJ7kYFwN+B2W4HEUAzcL0xZgJwEnBtiPz3agDONMYcA0wDZovISS7H5O2HwNrevqkmBo97gJ8RQmuLG2MOeh0OJERiM8a87rUu94fYlfZcZ4xZa4xZ73YcjulAqTFmszGmEXgKmOtyTBhjlmDXTAkpxpjdxpgVzn4N9pddwPXe+5OxDjmH8c6fkPh3KCL5wAXAX3v73poYABGZA+w0xnzqdiz+ROR2EdkBXEnoPDF4+ybwittBhKBhwA6v4zJC4BddOBCRQuBYYJm7kVjO65qVwD7gDWNMSMQF3Iv9Mtva2zfucgW3SCEii4AhAU79EvgFcE7/RmR1Fpcx5gVjzC+BX4rIjcB1wG9CIS6nzi+xrwAe74+Ygo0rREiAspD4phnKRCQFeA74kd8Ts2ucVSenOW1p/xKRycYYV9toRORCYJ8x5mMROaO37x81icEYc3agchGZAhQBn4oI2NciK0RkujFmj1txBfAE8BL9lBi6iktErgYuBM4y/TgYphv/vdxWBhR4HecDu1yKJSyISDw2KTxujHne7Xj8GWOqRORtbBuN2433pwBzROR8IAlIE5F/GmO+2hs3j/pXScaYz40xg40xhcaYQuw/6OP6Iyl0RUTGeB3OAda5FYs3EZmNXZt7jjGm1u14QtRyYIyIFIlIAnYd9IUuxxSyxH4rexhYa4z5o9vxtBGRnLZedyKSDJxNCPw7NMbcaIzJd35nzQMW91ZSAE0Moe5OEVklIp9hX3WFRBc+4H4gFXjD6Ur7Z7cDAhCRL4pIGTADeElEXnMrFqdx/jrgNWxD6jPGmNVuxdNGRJ4EPgDGiUiZiMx3OybHKcBVwJnO36mVzrdht+UBbzn/Bpdj2xh6tWtoKNIpMZRSSvnQJwallFI+NDEopZTyoYlBKaWUD00MSimlfGhiUEop5UMTg1JKKR+aGJRSSvn4/+6pt/1jB/SkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(model.μgrid, μ.T);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def A_operator(model):\n",
    "\n",
    "    μ, μgrid = model.μ, model.μgrid\n",
    "\n",
    "    @njit\n",
    "    def A(μ):\n",
    "        return np.sum(μ @ μgrid)\n",
    "\n",
    "    A(μ)\n",
    "    \n",
    "    return A\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_model(model,\n",
    "                T,\n",
    "                G,\n",
    "                A,\n",
    "                q_bounds=None,\n",
    "                tol=1e-4,\n",
    "                max_iter=2000,\n",
    "                reset=True,\n",
    "                use_existing=False,\n",
    "                skip_policy=False,\n",
    "                solve_skip=100,\n",
    "                verbose='all',\n",
    "                print_skip=100):\n",
    "\n",
    "    start = time()\n",
    "\n",
    "    if q_bounds is None:\n",
    "        q_bounds = [model.β, 1 / model.β]\n",
    "\n",
    "    verb = verbose == 'all'\n",
    "\n",
    "    qlow, qhigh = q_bounds\n",
    "\n",
    "    if reset:\n",
    "        model.init_value()\n",
    "        model.init_distribution()\n",
    "\n",
    "    def objective(q):\n",
    "\n",
    "        if verbose == 'all' or verbose == 'summary':\n",
    "            print(f\"Evaluating solution at q = {q:.6f}.\\n\")\n",
    "\n",
    "        v, g = solve_value(model,\n",
    "                           T,\n",
    "                           q,\n",
    "                           tol=tol,\n",
    "                           max_iter=max_iter,\n",
    "                           use_existing=use_existing,\n",
    "                           skip_policy=skip_policy,\n",
    "                           solve_skip=solve_skip,\n",
    "                           verbose=verb,\n",
    "                           print_skip=print_skip)\n",
    "\n",
    "        μ = solve_distribution(model,\n",
    "                               G,\n",
    "                               tol=tol * 0.1,\n",
    "                               max_iter=max_iter,\n",
    "                               use_existing=False,\n",
    "                               verbose=verb,\n",
    "                               print_skip=int(print_skip / 4))\n",
    "        B = A(μ)\n",
    "\n",
    "        if verbose == 'all' or verbose == 'summary':\n",
    "            print(f\"The aggregate balance for q = {q:.6f} is B = {B:.6f}.\\n\")\n",
    "            print(\"-------------------------------------------------------\\n\")\n",
    "\n",
    "        return B\n",
    "\n",
    "    qstar, result = root(objective,\n",
    "                         qlow,\n",
    "                         qhigh,\n",
    "                         xtol=tol * 0.1,\n",
    "                         full_output=True)\n",
    "\n",
    "    Bstar = objective(qstar)\n",
    "\n",
    "    iters = result.function_calls + 1\n",
    "\n",
    "    end = time()\n",
    "\n",
    "    runtime = end - start\n",
    "\n",
    "    print(f\"Model converged in {iters} iterations and {runtime:.4f} seconds.\")\n",
    "\n",
    "    model.set_solution(qstar, Bstar)\n",
    "\n",
    "    return qstar, Bstar\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = A_operator(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluating solution at q = 0.993220.\n",
      "\n",
      "The aggregate balance for q = 0.993220 is B = 2.052284.\n",
      "\n",
      "-------------------------------------------------------\n",
      "\n",
      "Evaluating solution at q = 1.006826.\n",
      "\n",
      "The aggregate balance for q = 1.006826 is B = -1.600514.\n",
      "\n",
      "-------------------------------------------------------\n",
      "\n",
      "Evaluating solution at q = 1.000865.\n",
      "\n",
      "The aggregate balance for q = 1.000865 is B = -0.765006.\n",
      "\n",
      "-------------------------------------------------------\n",
      "\n",
      "Evaluating solution at q = 0.998150.\n",
      "\n",
      "The aggregate balance for q = 0.998150 is B = -0.062666.\n",
      "\n",
      "-------------------------------------------------------\n",
      "\n",
      "Evaluating solution at q = 0.997858.\n",
      "\n",
      "The aggregate balance for q = 0.997858 is B = 0.040579.\n",
      "\n",
      "-------------------------------------------------------\n",
      "\n",
      "Evaluating solution at q = 0.997970.\n",
      "\n",
      "The aggregate balance for q = 0.997970 is B = 0.000295.\n",
      "\n",
      "-------------------------------------------------------\n",
      "\n",
      "Evaluating solution at q = 0.997972.\n",
      "\n",
      "The aggregate balance for q = 0.997972 is B = -0.000296.\n",
      "\n",
      "-------------------------------------------------------\n",
      "\n",
      "Evaluating solution at q = 0.997971.\n",
      "\n",
      "The aggregate balance for q = 0.997971 is B = 0.000000.\n",
      "\n",
      "-------------------------------------------------------\n",
      "\n",
      "Model converged in 8 iterations and 1.1749 seconds.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.997971053979755, 1.3600828877025073e-07)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve_model(model, Tp, G, A, verbose='summary')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD4CAYAAADo30HgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXhV1bn48e+beU7IxAwJEIYwCBhxAGdUVBRnwZ9Wrb3aVq3V9latt621w9XWqq1Dna/WCXGqOFQcUHAAIcwzhDmMIQmZ56zfH2uHcxIChJBkn+H9PE+etfY6++zzHkjOe/baa68lxhiUUkqpJiFuB6CUUsq3aGJQSinVjCYGpZRSzWhiUEop1YwmBqWUUs2EuR1AR0hNTTUZGRluh6GUUn5l0aJF+4wxaS3bAyIxZGRkkJub63YYSinlV0Rka2vt2pWklFKqGU0MSimlmtHEoJRSqhlNDEoppZrRxKCUUqoZTQxKKaWa0cSglFKqGU0MSvmSle9CZZHbUaggp4lBKV9RtgfevhHevM7tSFSQ08SglK+or7Zl0SZbrv0Ylr/lXjwqaLUpMYjIJBFZJyJ5InJPK49HisibzuPfi0iG12P3Ou3rROQ8p62viHwpImtEZJWI3OG1f7KIfCYiG5yy27G/TaX8QF2lLU2DLadPg3d/ZOtLXoMv/xcq9sErl8L+7e7EqILCERODiIQCTwLnA9nANBHJbrHbTUCxMWYQ8CjwkPPcbGAqMByYBDzlHK8e+IUxZhhwEnCr1zHvAb4wxmQBXzjbSgW+2gpbNjYc/Nj7P4U5D8LCF2DjbHhsRNfGpoJKW84YxgF5xphNxphaYDowpcU+U4CXnfrbwNkiIk77dGNMjTFmM5AHjDPG7DLGLAYwxpQBa4DerRzrZeCS9r01pfxMbbktK/dBVbGnfddyr/pST333SnhjGuzb0DXxqaDRlsTQG/A+b83H8yF+0D7GmHqgBEhpy3OdbqcxwPdOU3djzC7nWLuA9NaCEpGbRSRXRHILCgra8DaU8lE7l8KSV6GmzNP2eI6n/sypnvq6jz31p8fb7SdywJjmx1z/KVSXdk68KuC1JTFIK22mjfsc9rkiEge8A/zcGHNUv8XGmGeNMTnGmJy0tIOmE1fKP5TugmdPh/dvhTevtW1J/e1ZA8DUNyAq0dZPvs2WWedCSHjz4yx9Df7vAqgohD2r4PUrYda9XfMeVMBpy3oM+UBfr+0+wM5D7JMvImFAIlB0uOeKSDg2KbxmjHnXa589ItLTGLNLRHoCe4/i/Sjl+4o2wUuTISIWqksOfvz2RTDvSdu1NOR8uH0JFOZB33Ew9nqITgIJgWdOg0Fnw+J/2cQC8PlvYcCZtq73Q6h2aktiWAhkiUgmsAN7MfmaFvvMBK4H5gFXALONMUZEZgKvi8gjQC8gC1jgXH94AVhjjHnkEMd60Cnfb9c7U8rX7FkNL06CPjlQusPTftw0uPRp+O5xyDgVQsNhws89j8em2B+AtMGe9jtXgQjEpMI3zp/Rmg+gsdHWE/t07vtRAeuIicEYUy8itwGzgFDgRWPMKhF5AMg1xszEfsi/IiJ52DOFqc5zV4nIDGA1diTSrcaYBhGZAFwHrBCRpqtpvzbGfIxNCDNE5CZgG3BlR75hpVyx7E1472Zb3/gFpGfbawol2+Gkn9r2U24/umOK01N71m+g+3B7H8T7t8Ky1217XZXtWqopheTMjnkfKiiIaXnRyg/l5OQYXdpT+ayqYngoo3nbHcugWwbUVUN4VMe8TmMDPHsG7HZGMQ272L72lq/h7i0QrbcEqeZEZJExJqdle0Cs+ayUT9v89cFt3TJs2VFJASAkFK77N+S+CJu+sqOdSrbZxwo3QZ/jO+61VEDTKTGU6gz1tZD3BWz9Dr5+GOJ7wZAL7GPXf9h5rxubAqf/N4z7L09SAKjQId2q7fSMQanO8MZUey2hyUV/hzE/ANMIoV3wZzf8Esi/DfI+h4K19mL3u7fAkEkw/NLOf33l1/SMQamOUrAOFjxn+/q9kwLY/v6QkK5JCk3O+xPcPAdCI+CLB2D5dJjz1657feW39IxBqY7y7JlQVwFhkXa711jYuRgiEyEm2Z2YwqNg7A9g4fN2u6HGnTiUX9EzBqU6QkWhTQoAXz8CacPgDGf+x5QB7sUFcOZ99vpGt0zYv82e0dRrglCHpolBqY7wxtWeevFmyJoImafByKvsB7ObYpJh2hsw4U5oqIXXr4KHMqGm3N7roFQL2pWkVEco3Nh8e9RUCI+Gy59zJ57W9D/Flnmf2/KRbKgpgftbmZZDBTU9Y1CqLRY8B+/fZrthlr/V/Jv23jVQVQSJ/TxtaUO7PsYjSc2C8V5TbdQ4CaG19R9UUNPEoFRbfPxLWPIKbJ5rV1V7eoK9kcwY+ORee1fxFS969u/K0UdH45zfw693QYhXfCW6GpxqThODUkfjrettWZgHH97p3GX8pZ31tLef3FkcEWPnamqyfpZ7sSifpIlBqaPRcprsj+6y5fE32PsUxt8BV7/a5WEdtbE/8NR3LrUL+2iXknJoYlDqSCr2HXmfJOf6wjkPwLCLOjeejjDmOhh3s+1SWva6Xdjn27+7HZXyEZoYlDqc+lpY/4mtn9/iruFsr+XIQ0K7LqaOEB4FF/wVxt3iadu9wr14lE/RxKBUa774A6z9CD77jV3jQEJgtNf6VMMvgwtbrjHlh7xnXC3XxRKV5aNDJ5Ry2dcP27LnaFumDYPIOOgxCnoeB1OesO1n/w5SB7d+DH8wdLKnvnuFHWUlrS3VroKJJgalWvK+R6Fp3qOpzgXlH7dYW+HUu7omps4SFgk/mQcr3rLLg+5bb6fOCItwOzLlIu1KUqqlqmJPffv3EJsGyS7Pd9SZumfDMOfM4clxMOMHh99fBTxNDEoBlO603SjQPDGAvQAd6Jrd1/AfaGx0LxblOk0MSu1cAo8MgyVOd9HeNbY8/W5b1gTBXELh0TBoome7cIN7sSjXaWJQal+eLTfOtuXaDyE2HUZdfejnBKKpb8A1M2w9P9fdWJSrNDEoFRpuy4Za25206SvIOge6ZbgZVdcLi7BnDRHxkPsCrP3Y7YiUSzQxKNU0PLNsN+xaaq8x9DvJ/25a6wghodB7DOxYBNOnQV212xEpF2hiUKq20pY7cuHT39j6kAtt+eNv4bYg61bpe6Knnr/QvTiUazQxKNW0JCfAlq8hNNKzRnOPEXYdg2By2q/gkn/a+o5F7saiXKGJQanaiubbMcnBffdvWISd/qNbpj2LUkFHE4NSe1Y1345OdicOX9P7eFjzAdyfCM+d5XY0qgtpYlDBrWIfLH/T1s/6H1uGR7kXjy/xvtagXUpBRRODCm4Fa205dLLn7t+acvfi8SWjpzWfWryhzr1YVJfSxKCC2+6Vtjzvz5A6xNZ1AjkrMh6uetmzXV3qXiyqS2liUMFp7l9h6Ruw+GU7lXZSPztR3oS74PIX3I7Ot1z6jC2r97sbh+oyOu22Cj75i2D2Hz3bZ/2PHYUkAhN/515cvioqyZaPj4W71kJCT3fjUZ1OzxhU8Hm+xQibpg8+1bq+4zz1R4bqtYYgoIlBBZemqbW9RcR2fRz+JCYZ7vQa0luxz71YVJfQxKCCS8u1FgAi4ro+Dn+T2MdzraGy0N1YVKfTxKCCS9EmW06b7mnTM4a2SehtS00MAU8TgwoeOxbB61fZeorX/EeR8e7E429i02z5r4th/3Z3Y1GdShODCh7PnWW/7fY9CVIHeS46x+somzZJHQwDzrD1uX91MxLVyTQxqMDX2Agl+Z7tnqNsOW06TH4Mkvq6E5e/CQmBH7wPI6+CZdN1rYYA1qbEICKTRGSdiOSJyD2tPB4pIm86j38vIhlej93rtK8TkfO82l8Ukb0isrLFse4XkR0istT5uaD9b08FvYY6ePtGeHS4py3rXFv2PxlybnQnLn827CJoqIG9q468r/JLR0wMIhIKPAmcD2QD00Qku8VuNwHFxphBwKPAQ85zs4GpwHBgEvCUczyAl5y21jxqjBnt/Oj6gqr9PrwTVv/b1uN62JlTs85xNyZ/lzzAlvu3uRuH6jRtOWMYB+QZYzYZY2qB6cCUFvtMAZomVXkbOFtExGmfboypMcZsBvKc42GMmQsUdcB7UOrQVs/01Mt3Q9oQ92IJFEn9bPnWDXY5VBVw2pIYegPeQxDynbZW9zHG1AMlQEobn9ua20RkudPd1K21HUTkZhHJFZHcgoKCNhxSBaWW6zaH6ZTaxywqAc51phSZ+7C7sahO0ZbE0NpSVi1vHz3UPm15bkv/BAYCo4FdwN9a28kY86wxJscYk5OWlnaEQ6qgJS1+xcOj3Ykj0JxyO4yaCstnQGOD29GoDtaWxJAPeA/b6APsPNQ+IhIGJGK7idry3GaMMXuMMQ3GmEbgOZyuJ6WO2u4VULkPzvg1jLratukZQ8fJmAA1JbB/q9uRqA7WlsSwEMgSkUwRicBeTJ7ZYp+ZwPVO/QpgtjHGOO1TnVFLmUAWsOBwLyYi3oPKLwVWHmpfpQ5r/tO2zJoIMSm2Hh7jXjyBpul6zc6l7sahOtwRE4NzzeA2YBawBphhjFklIg+IyMXObi8AKSKSB9wF3OM8dxUwA1gNfALcaoxpABCRN4B5wBARyReRm5xj/UVEVojIcuBM4M4Oeq8q2BTmQY+Rdu3iprubddnOjtNzNEQmwHu3QIVOkxFIxLQ226SfycnJMbm5uW6HodxWvhc2z7XLdZ76S/jbYBhxBUx+BPaugadOghs+hozxbkcaOFbPhBnXwdm/g1PvcjsadZREZJExJqdluy7UowLHa1fArmW2npIF1SWQPsxupw+D+0vciy1QZV8M/SfA8jc1MQQQnRJDBY49qz31fOdSVnKmO7EEk4Fn2rO01qY0V35JE4MKHN5DURc+b8v4Xu7EEkz6nGDLHYvcjUN1GE0MKnC0NhQ1QRNDp+s91i529Olvob7W7WhUB9DEoAJHy7ucw2MgKtGdWIJJZDxMedJOqrf2Q7ejUR1AE4MKDMbYUUne4nuCtHbzvepwwy6GxL52JtumbjzltzQxqMDwyb1gGuCs/4Gzf2vbYlPdjSmYhIR4RoB99At3Y1HHTBODCgzrnNnZB0+CCOdmtvge7sUTjIZe6HYEqoNoYlD+qbYCcv8PVrxtJ3Er2wXj77B3OjdNC515ursxBpvjb4D04ZAyyO1I1DHSG9yUf3o8B8qc+Rh7jIKGWntTG8CQSfDLPO1KcsOoK+Hz+6FwI6QMdDsa1U56xqD8U5nXJL07nOlQUrM8bXFpeuHZDcdNA8Qmh/oat6NR7aSJQfmflvN77V1jy6YuJOWe+B4QGgFrZsK8J9yORrWTJgblf+qqmm/v2wAIxKa7Eo5q4aLHbFmS724cqt00MSj/U91iMrx96yEuHUL1kplPGH0NJPWH6lK3I1HtpIlB+Z+WiaFoow5N9TXdR8DW76Cx0e1IVDtoYlD+56kTD26L73lwm3LPiMvsAIFt89yORLWDJgblX7wnafvRbAiNtPW47u7Eo1o3+Dxbvnq5TqznhzQxKP9Sst2WSf2g1xjPvQppQ92LSR0sMh7G3QL1VbD1G7ejUUdJE4PyL8WbbXnpM3Z+nvP+bNdcyDrH3bjUwU77pS0L1rsbhzpqmhiUb9uXB/cnwrb5drt4iy27Zdhy+CXwizXNb25TviE2DeJ62HW4lV/RxKB825JXbPnieVBRCLtX2EVh4nQUks8TgZFXwIZZ9v9O+Q1NDMq3ed/M9u1jsO176D/ediMp33fcNGishzkPuh2JOgr616X8R2UhlO+GxD5uR6LaqscIGDoZlr0JDfVuR6PaSBOD8m21FZ562W6oKrZ3OSv/MeoqqCmB/AVuR6LaSBOD8m21ZZ76Pmd0S2yaO7Go9hlwBoSEwYbP3I5EtZEmBuW7lr4Oq9/3bDfdw6BnDP4lKhH6ngjrZx08M67ySZoYlO/690889cGTPHWdRdX/DL8U9q6yyUH5PE0Mync1rd3cOwfGXOdpj9OuJL9z/A0Qkwor33E7EtUGmhiUbzAG8nObdzVEJ9ny6ldh4Jme9oTeXRubOnah4ZB5Gmz5WruT/IAmBuUbNnwGz58NC5/3tFUWwUm3QkJPiIiFny2B2xZBWKR7car2yzwVynZ5BhEon6WJQfmGukpbrplpy5IdUFcBSX09+yQPgNRBXR+b6hhZ5wECy6a7HYk6Ak0Myjc01NmyscGW6/9jy4wJ7sSjOl5ibzu31YJnD15sSfkUTQzKN9Q5N7IZZ8Wv/Fw7a2r3Ee7FpDreuJuhthw2zXE7EnUYmhiUb6h1upKaEkN1KcQk24nYVODocwJEJkDe525Hog5DE4Ny367lMOteW9/+PbzzI6gptR8gKrCEhsOA0yHvCx2d5MM0MSj3vX5V8+0Vb9m7nKM0MQSk7EugNB/Wfex2JOoQNDEo97X2zbF4i54xBKrsSyCxHyx6ye1I1CFoYlDui4j11K97z1OP7tb1sajOFxpml2Ld8Cmseu/I+6sup4lBuS++p6fe72RPPVHvcA5YTcOQ37rB1TBU69qUGERkkoisE5E8EbmnlccjReRN5/HvRSTD67F7nfZ1InKeV/uLIrJXRFa2OFayiHwmIhucUr82BpKqYnj7h/biY40zpXZDrS2nvgHh0dDvFLutk+UFruGXQnq2nQ9LF/DxOUdMDCISCjwJnA9kA9NEJLvFbjcBxcaYQcCjwEPOc7OBqcBwYBLwlHM8gJectpbuAb4wxmQBXzjbKlCs+8ROpPbqZfCPsbatfDeMuhqGXmC3r34Vcm6CrHPdi1N1LhE49Rd2vY3dy92ORrXQljOGcUCeMWaTMaYWmA5MabHPFOBlp/42cLaIiNM+3RhTY4zZDOQ5x8MYMxcoauX1vI/1MnDJUbwf5evqqz31ir1QVw2lu5p3J8WmwORHbKkCV//xtlz2hrtxqIO0JTH0BrZ7bec7ba3uY4ypB0qAlDY+t6XuxphdzrF2AdqfEEjqa5pvL30VGuugW4Yr4SgXJThfBhY8C+UF7saimmlLYmjt1tOW4wsPtU9bntsuInKziOSKSG5Bgf5S+Y36KltOcyZS+/oRO/po1NXuxaTcM/4OW26f724cqpm2JIZ8wGuKS/oAOw+1j4iEAYnYbqK2PLelPSLS0zlWT2BvazsZY541xuQYY3LS0nThFr/RdMbQ7yQIi4bSHdBrDETEuBuXcseZ90FoJGzTxOBL2pIYFgJZIpIpIhHYi8kzW+wzE7jeqV8BzDbGGKd9qjNqKRPIAhYc4fW8j3U98P5h9lX+pr7afhBEd4ORl9u2WE3sQSssEnqP1cTgY46YGJxrBrcBs4A1wAxjzCoReUBELnZ2ewFIEZE84C6ckUTGmFXADGA18AlwqzGmAUBE3gDmAUNEJF9EbnKO9SBwjohsAM5xtlWgqKuGsChbP+lWW2ac6l48yn39x8OOXNix2O1IlENMAExklZOTY3Jzc90OQx1JfS38ZYD9lvirjbat6fdPZ1ENXqU74bGRMOgcuEYX8elKIrLIGJPTsl3vfFZd58s/2nHrlfs8bSKaFIJdQi8Y/f9g63eehZqUqzQxqK6zb4PbEShflXEq1JTA7hVuR6LQxKA6kzFwfyLMus9uV5fa8pav3YtJ+aYM52a3DZ+5G4cCNDGoztLYAHucabDmPWHLkm0w4groOcq9uJRvSugFgybCd//wfIFQrtHEoDrHN4/C084MmlFJdg3n/dugh67hrA5h/B125b5t89yOJOhpYlCdY+t3nnpiX1jxtq0PucCdeJTv63MChITD1m/djiToaWJQnaPpprWYFKgth72roddYSBviblzKd4VHQ98TYe1Huh60yzQxqM5Rvgd6Hw+n3A7Fm2HzHJ0oTx3Z2OugMA82feV2JEFNE4PqHNUldtqL0dfaswZoPrW2Uq3JvsT+vix5xe1IgpomBtU5asshIg7i0uxqXWDrSh1OeBQMPh/yPteV3VykiUF1jppyiIyz9RN/DCMuh6GT3Y1J+YfB59ozzi16v4tbNDGojjfzdijbCZEJdjs1C6540ZZKHUnWeRCTCt8/43YkQUsTg+pYdVWw+F+2HhHnbizKP4VHwehrbHeS3uzmCk0MqmOV7fbUQ8Lci0P5t8Hn2SVfdXSSKzQxqPYzxvYFeyvfY8teY+zQQ6Xao++JEJ0Mq951O5KgpIlBtd+SV+HBfnaivLkP27birbac8qSd/0ap9ggNh5FXwtqPoarY7WiCjiYG1X4Faz312X+A4i3w3s12O66HKyGpADJ6GjTUwKp/ux1J0NHEoNovNLz59mav4YUxyV0biwo8PUdDyiBY+Y7bkQQdTQyq/eqq7JDU2521er92upOyztVV2dSxE7H3v2z5xi7/qbqMJgbVfnVVduKzlIEw8irblQRwyT9dDUsFkFFXAwaWz3A7kqCiiUG1X10VhEXZ+vifedqju7kTjwo8KQPtdNwr33Y7kqCiiUG1X30VhMfYenevBXhCQt2JRwWmEZfbtaAL1rsdSdDQxKDaZ+5fYc0H9i5VsP3B138A/zXb3bhU4Bl+GUgoLHze7UiChiYGdfQa6mH2H229qSsJIPM0uwaDUh0pvjsMOR8WPgf78tyOJihoYlBts2+Dp+59w1FdZdfHooLPuX8A0wjr/+N2JEFBE4M6PGNg5bvwRA6sft+2VRZ6Hq+rdicuFVySB0BiP8hf6HYkQUETgzq8bx6Bt2+09R2LbNmUGKK76dBU1XUGnmGnyChY53YkAU8Tgzq8prMEgGXTbdk04+UNH0EfvaagusipvwTTAEtfczuSgKeJQR1e6hBPvXyPvXdh7l/sdmJfd2JSwalbf8g4FdZ94nYkAU8Tgzq8mrLm202jkQCiEro2FqWGnA/71sGe1W5HEtA0MajW1ZTB7D9B6Q7IPB1+tx+GXwrznrCPn3ybu/Gp4DT8Ursy4Ld/dzuSgKaJQbVuwXO2y2j3cohNtTewnX6P5/Ghk92LTQWv+B6QfQms+w/U17odTcDSxKBaFxbpqcek2DJ9qKctSa8vKJcMuwhqSmDzXLcjCViaGFTraso99abEADDmWlvG9+zaeJRqMvBM+zu54Fm3IwlYmhhU6yr3eepx3T31ix6He7brRHnKPWGRMPr/wcbZUF3qdjQBSRODcu5ufge2L/C0VXglhjSvIashIToaSblvyPnQWAfrZ7kdSUAKczsA5bKq/fBQf8/2XWugoQ5WvQtpw2Dcf0G/k92LT6nW9D3RTpMx/ykYdaXb0QQcPWMIdk2rrjV5aTLM+IGtRyXCCTfpMp3K94SEwok/hp2LYc8qt6MJOJoYgl3LtXSLNsKupbaeMrDr41GqrYZfBhICK95yO5KAo4kh2JXvPvRjF/6t6+JQ6mjFpcGgc2Dp67b7U3WYNiUGEZkkIutEJE9E7mnl8UgRedN5/HsRyfB67F6nfZ2InHekY4rISyKyWUSWOj+jj+0tqsOqLvHUpzwJo6baeo9REB7tTkxKtVXOD+0cXut0nYaOdMTEICKhwJPA+UA2ME1EslvsdhNQbIwZBDwKPOQ8NxuYCgwHJgFPiUhoG47538aY0c7P0mN6h6p1dVXw0S+ar6M7aCL0O9HWvUclKeWrss6BhN6w+GW3IwkobTljGAfkGWM2GWNqgenAlBb7TAGa/mfeBs4WEXHapxtjaowxm4E853htOabqTOv+Y9fQXfa63b76NTvdQP/xdjupn3uxKdVWIaH2psu8L6B4q9vRBIy2JIbewHav7XynrdV9jDH1QAmQcpjnHumYfxKR5SLyqIh4zc3gISI3i0iuiOQWFBS04W2oZrynvEgeCMOcuY/ShsAP3oerX3EnLqWO1pjr7Mi5xf9yO5KA0ZbE0NpYRdPGfY62HeBeYChwApAM3N1aUMaYZ40xOcaYnLS0tNZ2UYfjPZ12eEzzxwacAXHpXRmNUu2X1BeyzoUlr+hF6A7SlsSQD3jPmNYH2HmofUQkDEgEig7z3EMe0xizy1g1wP9hu51UR/OeSkDvZFb+7sBF6I/djiQgtCUxLASyRCRTRCKwF5NntthnJnC9U78CmG2MMU77VGfUUiaQBSw43DFFpKdTCnAJsPJY3mDQ278NPvg5VBXDrPugfK9t9x6NNPxSd2JTqqMMmmhXFFz4vNuRBIQjTolhjKkXkduAWUAo8KIxZpWIPADkGmNmAi8Ar4hIHvZMYarz3FUiMgNYDdQDtxpjGgBaO6bzkq+JSBq2u2kp8OOOe7tB6K0bYUcuhEbAgmfst6rLn4fKQoiIh7u3QKjOjKL8XEioPWv44vewdw2kD3M7Ir/Wpk8EY8zHwMct2n7rVa8GWp2wxBjzJ+BPbTmm035WW2JSbdTgLGay4BlbluywZdkuSOipSUEFjrHXw5yH7HTckx91Oxq/pnc+B7rIFtcPtn0Hu5bDtvl2eKpSgSI2BUZeAcum265T1W6aGAJdfZWnPvkxW7482U6F0WOUOzEp1VnG3QJ1lbDkNbcj8WuaGAKd9+ijnBvhpJ96Ljzn/NCdmJTqLD1HQb9TYOFz0NjgdjR+SxNDoKsphZBwOPePdrv7cM9j3TLdiUmpznTizXY6+Q2fuh2J39LEEKgKN8Ljx9tRSKfcBqfcbtsHnGHLHqPsamxKBZqhk+38Sd8/7XYkfks/GfzZtu/htatav9vz279DYZ6tx3rdGZ7YB2752k57oVQgCg2HcTfDpq9gzYduR+OXNDH4s3d/BBtmQdGm5u21FbD6357t2BZThvQcBTHJnR+fUm45+Va79OfC59yOxC9pYvBnEXG23L8d3vkRrHVuC9k2r/mdzXHduz42pdwUGg7DLobNcyE/1+1o/I4mBn8WEWvLnYvt8obTp9ntzXPtBeemM4VeY9yJTyk3nXoXRCXB/KfcjsTvaGLwZyHhtvzS68bywo22X7XPCfCT7+D6D3SSPBWcohJhxGWw6j3Ys9rtaPyKJgZ/Vlt2cNvjY6FoIww43U6dnXla18ellK848z4IjbTTZKg208Tgz7zXVGh5F/OgiV0bi1K+KCYZRl1pF/Ep3Oh2NH5DZ1DzZ96J4ZKn7EXomGS7nnOfHB9Qty0AABf+SURBVPfiUsqXnPk/sPwtmPswXPpPt6PxC0GdGKoqyqitqSYx2cdXgDPGLl3ovQ02MYy4HLIvgR4j7Y9Sqrn47nDCTfYi9Gm/hJSBbkfk84K6K2n5i7fR+I+xbodxsMKN8MhwKFhn679P8gxFBXj3Zniwv51SOz0bsi92L9YAUFxRizEtV6tVAeWUn9lrDXP/6nYkfiGozxiQUEJo50Rb1SVQX9M5ayN/+xiU5sP6TyAs2rat+QCGXgC7V8KKGZ59W968ptrkidkb6J4QxckDU5jw0JekxEYw+xdnkBgTTnFFLU/P3cjx/brRPyWWwooa1uwqY/3uMsZlJlPb0Ej/5BjCw0IYmBZHRU09K3aUUFBWw+Du8Zw8MAVjDMZASIg90yssryEiLIT4qHCX33mQ8j5rmHAXpA12OyKfFtSJwYSEEmIa2/fkv4+GqiK4v+TI+x6tpm+voZFQuc/Wl70Ok/4M6//TfN/UrI5//QDS2GjYX1VHcmwEjY2GRmOYv6mIhz9d32y/wopaxv35c04ZmMKCzUVU1Lb+heHN3O1HfM2eiVHsKqkmIjSEHolRNDQaduy3058PSI21yaW+kYHpcaTFRZIQHc6ukiomjejB2t1ljO6TRLfYiGN/86q58T+HRS/DZ7+Fa6a7HY1PC+rEgIQSSjsTQ1XR0e3fUAef3w8T7oTY1MPvGxrheY2KAk/713+DLd9C95Gwx1kKu/uIo4sjgH25bi+fr97D5cf3IbtnAm8tyuelbzezsaACgAFpsWxy6k3b1bUN3HzaAEqr65mzvoA56ws4ISOZ28/KYk9pNblbiygoq+XnE7PonxLD6p2lpCdEsX5PGSEirN9Thgh0i4kgKz2Ov8xaR0RoCFce34eSqjpytxaTFh9JbUMjYSHCpn0VFJTXEBoivLtkR7P4f/+BZ6z9GUPSiI8KJz4qjNMHp5EaF0lidBipcZGEhYYQFxmGMQbxvvakDi8uDU77hf073PSVZ0JJdRAJhL7VnJwck5t79Le9z3vmVsbsfJOo3+87+he9P9Ep23jGsOo9eOsGGHMtTHny8Pu+fysseRX6nWynt2gSHmMXITn3j5B1LjTWN59GO0jU1Ddw/YsLmL+piAmDUrnn/KH8ZdY65q4vOGjf1LhI9pXXECLQ6Pyq/+6ibM4f0ZMeiVEH7V9V20BUeEiXfOBW1tZTWF5LUUUtJVV1rN5Vypx1BfTpFs2y/P1s3ldBXcOh/z5T4yK5Y2IWp2el0S8lptPjDQh11fDECfbmt1vm2LWig5iILDLGHDSEUc8Y2nvGcLQOzF3Uhg+cOmfVNe+kkHk6bJ5j633GQdqQDg3P1zQ2Gj5asYst+yoorKilR2IUQ3vEk7ulmCe+zDuw3zd5+5j8+DckRoczcVg6PROjWbiliMKKWm4/axDXndQfEcEY253TKzH6QL9/a6Ijuu6DIiYijJjkMPom2w/10wan8ePTm4+Y2V9Zy/aiKnbsr2R7URWrdpYQGRbK52v2sK+8ht/8eyVhIcI52d1JjYtkQlYqJ2YmU13X2GriC3rhUXDO/fD2D2Hp6zD2Orcj8knBnRhCujAx1DpdGE0T3x1238rm23etsfMfNSWG9KEdG1sXqmto5OfTl3LD+AxOyPDM8FpSWcdjX6zng2U7SYuPYlB6HB8s29nqMU7NSiWnfzJnDEmjrLqeBVuKuHRMbzJTYw/5uiJCn27+9606KSaCpJgIRvZJPOixqtoG1u0p4/Xvt/JtXiHFlbW8Mn/rgccvGNmDCYPSSImL4Nzs7trt1GT4ZfD9s/D572DohTrTcCuCOzFIKCFiMI2NSGcvWlPnfNiHO6OMyvfaGSCju7W+b3wvKNtpH0/oBUPOt4+lDbWnwX6mqraBx75YT3hICB+t2MVHK3YxcVg6e8tqWJ7fvDtuX3kta3aVMi7TfvifPbQ7MRGh7NhfRd/kGHonRTfbf0LWEa7ZBKjoiFBG901idN8kwCbdhZuLWLJ9P6t2ljB77V4+XrEbgIFpsfRKiuaMIelcNKon3WIjCA8N0tHqIjD5EXjmNPjsN0fu2g1CwZ0YnP7FhoZ6wkI6eRRIQ70tjTPa5eEse83gvl12u7EB6qvtjKl1Vbar6Kz7bLcR2GRw91a/Wce2sdGwu7SauKgwEqLCefHbzTwzx7NuRGxEKGt2lR0YrQNw04RMLhjZk6raBjYWlDN1XF8iwzxdO01dLqp14aEhnDIolVMG2URZU99AYXkts9fuZdaq3eworuIPH67mDx+uJiYilFMGptIvOYZbTh9A94Qg63bqPhxOvs0ODT/uGsgY73ZEPkUTA05iCO/kxFDvfAB6dxPVedU//Lmdz2XgWbB/G/Q70V6o9had1LkxtsPmfRWs3FHCRcf1whjDf1bupqq2gf9+e9mBi72ZqbFs3ldBTv9uXHRcL5JiwpkyujfGGPaU1rB0ezGnDU4jJsLz6xisZwEdKTIslF5J0Vx7Un+uPak/xhgWbS1mxY4ScrcW89Fy+6XklflbiA4P5eoT+jJldG/6JseQGB0E91ucfjesetf+7d3ytb3+oIAgTwwiTmKob2VpzI5WX2PLukrPxWWA6lKoKrZJAWDjbFt2y+z8mI7R9qJKznz4KwAWbS0mOiKUf37VfKIy27dtx+//7qLhzUbPiAg9EqOYlNizK8MOWiJCTkYyORnJ3Dg+k4evaKCgrIbnv9nE/E2FPPf1Zp77ejPR4aFcdFxPhvdK5PyRPUiPD9APzIgYmPwovHo5zP4DnPenIz8nSAR1YvCcMRxD90xjQ9uGvNVX27J0B5Tt8rQ/2Lf1/bPObX9MHaCh0TB3fQFnDEkjv7iKn7y2iAcvG8Wwngk8PWcjr8zbyu7S6gP7/2velgNnCOdkd+fG8Rkc379bs64g5VuiI0LplxLDA1PsvTDbiypZuKWIb/MK+feSnczIzef3H6xiUHocx/dP5qdnDKRPt+jAuog9aCKc8F8w7wnIOkfvbXBoYgAajyUxNNS1MTE4Zwx710Jp66NtuPIlqK+1SaQL+jw3FpSTX1zF6YPttBq19Y28/N0WsnslcO+7K9hWVEnvpOgD1wEmP/5Ns+eP7J3I3ZOGMiErldLqOpZvL6F/SoxeC/BTfZPt/91lY/vw8JWj2FhQzsylNkG8sWAbbyzYRnp8JFNP6MvE7O6M6uN7XZvtcs4DdsTfez+Bn3yro5QI9sQgTYmhvv3HaGzjc5vOGMp3w0sX2vqZ/wM1pfDdPyCxHwy/tP1xHCVjDGf/zQ5/7ZUYxdRx/Vi3p+xAv3MT74vD157UjxARdpdUk90rgZ+dlXXgnoCEqHC9LhBARIRB6fHcde4Q7jp3CCt3lDB/UyHTF27niS/z+MfsPIb2iOfUrFQuOq4X2T0TCPPXUU4RMXDZc/D8RHjzOrjuXQiLdDsqVwV1YhCvi8/t1tjG6xP1NXb92er9djuuO5x8K2CgaBOMvKL9MbRBSWUdVXUNlFbXkRQdzrg/f3HgsZ0l1TzymWfuoMiwELK6x/HAlBEUlddy6uBU7RIKciN6JzKidyI/OnUAJZV1vJm7jbnr9/HSd1t47uvNpMZFcP6InowflMrJA1JIjPGzi9e9RsOFD8MHd8Cs+2w9iAV1YiDEvv1jO2NwuqGKt9jb7Vu7+axsjz1jSBkIOxbZGVN/sc6zxsLU19r/+odRUVPPvvIabnxpYbM5gry9cH0O3WIjeG7uJsZlJnPdSf3995uf6hKJMeHcfNpAbj5tIHvLqpm3sZBPV+3hrUXbeWX+VuKjwjhraDpnDknntMFpJPvLhIDH3wB7VtllQHuMsNtBKqgTg4R0QFfS6n/DCT+Cvx9nt3+3337g11baC1pfOiMdopPt2OmfzLOPd8AFvNLqOhZtKSY+Kox1e8ooqarj7KHd+dU7y+mTFM1HK5p3C00cls7na/YSFiIkx0bwwvUnHLij9p/XHn/M8ajgkx4fxZTRvZkyujcVNfV8sGwnCzYXMWd9Ae8v3UloiHD20HQmDuvOucO7kxTj40nivP+1a6B89Et7M2m/k9yOyBVBnRg65OLzR7+AMV7zrZTuhMTe8Pjx9s7lJlVFdu2G7tltPnRDoyFvbzmDu8cdGAlSUlXH72eu4rqT+/Ojl3MprKht9py/fLIOgGXb9x9o++3kbC4b2/vAH2V1XQORYV0zUZwKHrGRYUwd14+p4/rR2GhYsaOEj1bsYubSnXy6eg93vwtj+iYxaUQPzs3uQcZhpjBxTWgYXPECPHcWvDEVrv8gKFdGDOrE0CFnDAAl+Z76o9me6SxaypjQ5kM+PWcjX63by/xNdnrvn54xkFU7S5njzCDacsrmaeP60S0mnIqaetITorhwZE/6JsdQ39h40PWBqHC9XqA6V0iIcFzfJI7rm8S95w9lxY4SPl+zl89W7+HPH6/lzx+vZVB6HGcPTefEAcmcMjDVd34vo7vBte/AS5Ph5Yud5BBc09sH9bTbuR88Q86iX7Htmjn0Gzz66J58fxvnK7pzNexdYxfauejvEBkP2G/+BWXVNDRCUkw4u0qqeW9xPkN6JPDh8p18t7HwkIe8cGRPSqvrCA8N4YlrxhAi4jt/VEodQX5xJZ+u2sOX6/Yyf1MhdQ2G+Mgwxg9K5ayh6ZwxJI10X5iio3CjTQ4NNTY5BOAU9zrtdisk1L79ytJDfwi3qrEBJARarP7WGBZNyKgrITYdKvbyavJtzPtoD7X1Sfz1iqd5e0E+STElzMjdzoLNR17o58bxGQxMiyM+KozVu0r5wckZ9EqM0i4g5df6dIvhhxMy+eGETHudbGsxs1buZu76Aj5ZZSf9y+6ZwKlZqUzISuWEjGR3vvikDIQbPrTDy1++CH7wftB0KwX1GcPiT15i7Pw7AGj8TREhoZ5fvt3b84iJ70ZCUsrBT6wsgr9k0nj8DwlprKOk5ymc8G4MtYSz8L6JpMVH8tHyXdz6+uLDvv6I3gnsKa1BgL1lNfRw1iAODxV+OCGToT0Sjvo9KeWvjDGs2VXGl+v2Mnd9AYu3FVPXYIgIC+GEjG6My0hh/KAURvZJ7Nrh001nDtUlcPlzdqruAHGoM4agTgxLPn2VMd/dCsDG0EwG/mYpC99/ir5jzqHHizlsC+lNv9+uZsFj02hMG0rWxJuoefpMKkb/iKwlf+aXjbdhRlzFO4vzmx03LjKM8pqDr1ukx0dyXN8kNhaUM3lkT+48Z7B++1fqECpq6lmwuYivN+zjm7wCNuwtxxiIiQhlbL9unDwwhZMGpDCydyIRYZ08xLp0F0y/BnYuhtPvgdP+216o9nOaGFqx9PM3GP3Njw+7T+7Yh8hZfHerj11fezdzGo87sB0VHkJtfSNj+nVj0dZiLj6uF3+4ZARxkWHs3F9FWnykXgtQqp1KKuuYt2kf320sZMHmItbuLgPs393wXomM6JXA8F6JDO+dQFZ6fMcni7oq+PAue72w70lwyVO2u8mPaWJoxbLZ0zlu7i1t2rfahBMlddSYcDZHZFHVGEr+2U+xsTKKkqo6fnZWFt28buSpb2jUG8WU6kT7ymvI3VLEgs3FrNxRwqqdJVTU2qHnEaEhDO4Rx9AeCQzpHk9W9zgGpccdcWnXNlk+ww5Tr6uCk34MKYPsqnDz/wl9jrcT8/mJY0oMIjIJ+DsQCjxvjHmwxeORwL+A44FC4GpjzBbnsXuBm4AG4GfGmFmHO6aIZALTgWRgMXCdMab5YP0W2psYCnZuIe3Z45q1VZgoYsXOa/T9sF9z4po/s2Dk/Yy7/E4K9+STlNKD0DD/P4VUKtA0Nhq2FlWyckcJK3eWsHpnKWt3l1FQVnNgn8iwEDJTY8lMjSUjNZbeSdH0SoqiV1I0PROjSYgKa1v3btlu+M/d9gZXgJhUqNxnl+69e4tdndEPtDsxiF20YD1wDpAPLASmGWNWe+3zU2CUMebHIjIVuNQYc7WIZANvAOOAXsDnwGDnaa0eU0RmAO8aY6aLyNPAMmPMPw8XY3sTA0B1ZTnr5n9IWuZokrv3oa6ultLCXRTl5zHytCntOqZSyncUV9Syfk8ZG/aWs2VfBVsKK9i0r4JthZXUNzb//IuLDCM9IZJuMREkRoeTFB1OQnQ4STHhdtspE6MjSIwKJW33HOJX/ouQDZ96DnL5C50+91lHOZbEcDJwvzHmPGf7XgBjzP967TPL2WeeiIQBu4E04B7vfZv2c5520DGBB4ECoIcxpr7lax/KsSQGpVRwamg0FJTVsLOkip37m36q2VtWTUlVHfsr6yipqqOkso6yVgaTNEmhhBmRf+TlhnO5IXQWA2QnW0MOsc7KEbSnY7/qvL8x7MTDfkQe0rHcx9Ab2O61nQ+ceKh9nA/0EiDFaZ/f4rm9nXprx0wB9htj6lvZv+Ubuhm4GaBfv35teBtKKeURGmJXEOyRGMXYft0Ou299QyOl1fU2UVTVsb+y9kC9pLKO16veJiE8lK/KTmfvtjeIaKw+7PHa7+DUkRAV1+Gv0pbE0FqHW8voDrXPodpbuyp7uP0PbjTmWeBZsGcMre2jlFIdISw0hOTYiDbMFDsE6Lp1VTpLW4bN5APe50V9gJYTAR3Yx+lKSgSKDvPcQ7XvA5KcYxzqtZRSSnWitiSGhUCWiGSKSAQwFZjZYp+ZwPVO/QpgtrEXL2YCU0Uk0hltlAUsONQxned86RwD55jvt//tKaWUOlpH7EpyrhncBszCDi190RizSkQeAHKNMTOBF4BXRCQPe6Yw1XnuKmeU0WqgHrjVGNMA0NoxnZe8G5guIn8EljjHVkop1UWC+gY3pZQKZocalaS35iqllGpGE4NSSqlmNDEopZRqRhODUkqpZgLi4rOIFABb2/n0VOz9E75G4zo6GtfR0biOjq/GBccWW39jTFrLxoBIDMdCRHJbuyrvNo3r6GhcR0fjOjq+Ghd0TmzalaSUUqoZTQxKKaWa0cTgTMTngzSuo6NxHR2N6+j4alzQCbEF/TUGpZRSzekZg1JKqWY0MSillGpGE4MXEfmliBgRSXU7FgAR+YOILBeRpSLyqYj0cjsmABH5q4isdWJ7T0SS3I4JQESuFJFVItIoIq4PLRSRSSKyTkTyROQet+MBEJEXRWSviKx0OxZvItJXRL4UkTXO/+EdbscEICJRIrJARJY5cf3e7Zi8iUioiCwRkQ878riaGBwi0hc4B9jmdixe/mqMGWWMGQ18CPzW7YAcnwEjjDGjgPXAvUfYv6usBC4D5rodiIiEAk8C5wPZwDQRyXY3KgBeAia5HUQr6oFfGGOGAScBt/rIv1cNcJYx5jhgNDBJRE5yOSZvdwBrOvqgmhg8HgV+RfvW4+4UxphSr81YfCQ2Y8ynXutyz8eutOc6Y8waY8w6t+NwjAPyjDGbjDG1wHRgissxYYyZi10zxacYY3YZYxY79TLsh12r6713JWOVO5vhzo9P/B2KSB/gQuD5jj62JgZARC4GdhhjlrkdS0si8icR2Q78P3znjMHbD4H/uB2ED+oNbPfazscHPuj8gYhkAGOA792NxHK6a5YCe4HPjDE+ERfwGPbLbGNHH/iIK7gFChH5HOjRykP3Ab8Gzu3aiKzDxWWMed8Ycx9wn4jcC9wG/M4X4nL2uQ/bBfBaV8TU1rh8hLTS5hPfNH2ZiMQB7wA/b3HG7Bpn1cnRzrW090RkhDHG1Ws0IjIZ2GuMWSQiZ3T08YMmMRhjJrbWLiIjgUxgmYiA7RZZLCLjjDG73YqrFa8DH9FFieFIcYnI9cBk4GzThTfDHMW/l9vygb5e232AnS7F4hdEJBybFF4zxrzrdjwtGWP2i8hX2Gs0bl+8Hw9cLCIXAFFAgoi8aoy5tiMOHvRdScaYFcaYdGNMhjEmA/sHPbYrksKRiEiW1+bFwFq3YvEmIpOwa3NfbIypdDseH7UQyBKRTBGJwK6DPtPlmHyW2G9lLwBrjDGPuB1PExFJaxp1JyLRwER84O/QGHOvMaaP85k1FZjdUUkBNDH4ugdFZKWILMd2dfnEED7gCSAe+MwZSvu02wEBiMilIpIPnAx8JCKz3IrFuTh/GzALeyF1hjFmlVvxNBGRN4B5wBARyReRm9yOyTEeuA44y/mdWup8G3ZbT+BL529wIfYaQ4cODfVFOiWGUkqpZvSMQSmlVDOaGJRSSjWjiUEppVQzmhiUUko1o4lBKaVUM5oYlFJKNaOJQSmlVDP/H85/r9ZKBFgTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(model.μgrid, model.μ.T);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model converged in 8 iterations and 0.7779 seconds.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.99797105445102, -3.121781626258535e-08)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve_model(model, Tp, G, A, skip_policy=True, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model converged in 10 iterations and 0.3600 seconds.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.997971371090888, -0.001244835981068601)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve_model(model, Tp, G, A, skip_policy=True, verbose=False, use_existing=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "# References\n",
    "\n",
    "Huggett, M. (1993). The risk-free rate in heterogeneous-agent incomplete-insurance economies. Journal of Economic Dynamics and Control, 17(5-6):953-969"
   ]
  }
 ],
 "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": 4
}