dynamic_programming/odu.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='odu'></a>\n",
    "<div id=\"qe-notebook-header\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Job Search III: Search with Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Job Search III: Search with Learning](#Job-Search-III:-Search-with-Learning)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Model](#Model)  \n",
    "  - [Take 1: Solution by VFI](#Take-1:-Solution-by-VFI)  \n",
    "  - [Take 2: A More Efficient Method](#Take-2:-A-More-Efficient-Method)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this lecture we consider an extension of the [previously studied](https://lectures.quantecon.org/mccall_model.html) job search model of McCall [[McC70]](https://lectures.quantecon.org/zreferences.html#mccall1970)\n",
    "\n",
    "In the McCall model, an unemployed worker decides when to accept a permanent position at a specified wage, given\n",
    "\n",
    "- his or her discount rate  \n",
    "- the level of unemployment compensation  \n",
    "- the distribution from which wage offers are drawn  \n",
    "\n",
    "\n",
    "In the version considered below, the wage distribution is unknown and must be learned\n",
    "\n",
    "- The following is based on the presentation in [[LS18]](https://lectures.quantecon.org/zreferences.html#ljungqvist2012), section 6.6  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model features\n",
    "\n",
    "- Infinite horizon dynamic programming with two states and one binary control  \n",
    "- Bayesian updating to learn the unknown distribution  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide-output": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1mActivated\u001b[0m /home/qebuild/repos/lecture-source-jl/_build/website/jupyter/Project.toml\u001b[39m\n",
      "\u001b[36m\u001b[1mInfo\u001b[0m quantecon-notebooks-julia 0.1.0 activated, 0.2.0 requested\u001b[39m\n"
     ]
    }
   ],
   "source": [
    "using InstantiateFromURL\n",
    "github_project(\"QuantEcon/quantecon-notebooks-julia\", version = \"0.2.0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model\n",
    "\n",
    "\n",
    "<a id='index-0'></a>\n",
    "Let’s first review the basic McCall model [[McC70]](https://lectures.quantecon.org/zreferences.html#mccall1970) and then add the variation we want to consider"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Basic McCall Model\n",
    "\n",
    "\n",
    "<a id='index-1'></a>\n",
    "Recall that, [in the baseline model](https://lectures.quantecon.org/mccall_model.html), an unemployed worker is presented in each period with a\n",
    "permanent job offer at wage $ W_t $\n",
    "\n",
    "At time $ t $, our worker either\n",
    "\n",
    "1. accepts the offer and works permanently at constant wage $ W_t $  \n",
    "1. rejects the offer, receives unemployment compensation $ c $ and reconsiders next period  \n",
    "\n",
    "\n",
    "The wage sequence $ \\{W_t\\} $ is iid and generated from known density $ h $\n",
    "\n",
    "The worker aims to maximize the expected discounted sum of earnings $ \\mathbb{E} \\sum_{t=0}^{\\infty} \\beta^t y_t $\n",
    "The function $ V $ satisfies the recursion\n",
    "\n",
    "\n",
    "<a id='equation-odu-odu-pv'></a>\n",
    "$$\n",
    "V(w)\n",
    "= \\max \\left\\{\n",
    "\\frac{w}{1 - \\beta}, \\, c + \\beta \\int V(w')h(w') dw'\n",
    "\\right\\} \\tag{1}\n",
    "$$\n",
    "\n",
    "The optimal policy has the form $ \\mathbf{1}\\{w \\geq \\bar w\\} $, where\n",
    "$ \\bar w $ is a constant depending called the *reservation wage*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Offer Distribution Unknown\n",
    "\n",
    "Now let’s extend the model by considering the variation presented in [[LS18]](https://lectures.quantecon.org/zreferences.html#ljungqvist2012), section 6.6\n",
    "\n",
    "The model is as above, apart from the fact that\n",
    "\n",
    "- the density $ h $ is unknown  \n",
    "- the worker learns about $ h $ by starting with a prior and updating based on wage offers that he/she observes  \n",
    "\n",
    "\n",
    "The worker knows there are two possible distributions $ F $ and $ G $ — with densities $ f $ and $ g $\n",
    "\n",
    "At the start of time, “nature” selects $ h $ to be either $ f $ or\n",
    "$ g $ — the wage distribution from which the entire sequence $ \\{W_t\\} $ will be drawn\n",
    "\n",
    "This choice is not observed by the worker, who puts prior probability $ \\pi_0 $ on $ f $ being chosen\n",
    "\n",
    "Update rule: worker’s time $ t $ estimate of the distribution is $ \\pi_t f + (1 - \\pi_t) g $, where $ \\pi_t $ updates via\n",
    "\n",
    "\n",
    "<a id='equation-odu-pi-rec'></a>\n",
    "$$\n",
    "\\pi_{t+1}\n",
    "= \\frac{\\pi_t f(w_{t+1})}{\\pi_t f(w_{t+1}) + (1 - \\pi_t) g(w_{t+1})} \\tag{2}\n",
    "$$\n",
    "\n",
    "This last expression follows from Bayes’ rule, which tells us that\n",
    "\n",
    "$$\n",
    "\\mathbb{P}\\{h = f \\,|\\, W = w\\}\n",
    "= \\frac{\\mathbb{P}\\{W = w \\,|\\, h = f\\}\\mathbb{P}\\{h = f\\}}\n",
    "{\\mathbb{P}\\{W = w\\}}\n",
    "\\quad \\text{and} \\quad\n",
    "\\mathbb{P}\\{W = w\\} = \\sum_{\\psi \\in \\{f, g\\}} \\mathbb{P}\\{W = w \\,|\\, h = \\psi\\} \\mathbb{P}\\{h = \\psi\\}\n",
    "$$\n",
    "\n",
    "The fact that [(2)](#equation-odu-pi-rec) is recursive allows us to progress to a recursive solution method\n",
    "\n",
    "Letting\n",
    "\n",
    "$$\n",
    "h_{\\pi}(w) := \\pi f(w) + (1 - \\pi) g(w)\n",
    "\\quad \\text{and} \\quad\n",
    "q(w, \\pi) := \\frac{\\pi f(w)}{\\pi f(w) + (1 - \\pi) g(w)}\n",
    "$$\n",
    "\n",
    "we can express the value function for the unemployed worker recursively as\n",
    "follows\n",
    "\n",
    "\n",
    "<a id='equation-odu-mvf'></a>\n",
    "$$\n",
    "V(w, \\pi)\n",
    "= \\max \\left\\{\n",
    "\\frac{w}{1 - \\beta}, \\, c + \\beta \\int V(w', \\pi') \\, h_{\\pi}(w') \\, dw'\n",
    "\\right\\}\n",
    "\\quad \\text{where} \\quad\n",
    "\\pi' = q(w', \\pi) \\tag{3}\n",
    "$$\n",
    "\n",
    "Notice that the current guess $ \\pi $ is a state variable, since it affects the worker’s perception of probabilities for future rewards"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameterization\n",
    "\n",
    "Following  section 6.6 of [[LS18]](https://lectures.quantecon.org/zreferences.html#ljungqvist2012), our baseline parameterization will be\n",
    "\n",
    "- $ f $ is $ \\operatorname{Beta}(1, 1) $ scaled (i.e., draws are multiplied by) some factor $ w_m $  \n",
    "- $ g $ is $ \\operatorname{Beta}(3, 1.2) $ scaled (i.e., draws are multiplied by) the same factor $ w_m $  \n",
    "- $ \\beta = 0.95 $ and $ c = 0.6 $  \n",
    "\n",
    "\n",
    "With $ w_m = 2 $, the densities $ f $ and $ g $ have the following shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXwU9f0/8PdnZjbn5k4gJCThDPcNCsgVriCQBMSrnoj1qE09vlWxiLS/Sr8IbbW1pbVUo4iI+kWuyB1iuC8RBQIE5QwhIQlsSHaTbHZmPr8/ooicYdndmdl5Pf/iMdnsvtndz7zyOWY+jHNOAAAAZiVoXQAAAICWEIQAAGBqCEIAADA1BCEAAJgaghAAAEwNQQgAAKaGIAQAAFNDEAIAgKkhCAEAwNQQhAAAYGq+DkLO+b59+5r+eFVVFUXxXj3QFLIs41Z82uKcy7KsdRVmpyiKqqpaV2F23jgd+ToIGxoa+vXr1/TH19bW1tbWeq8eaIqqqiqchbWlKIrNZtO6CrPD6UgPzp8/7/E/RzA0CgAApoYgBAAAU0MQAgCAqSEIAQDA1BCEAABgaghCAAAwNQQhAACYGoIQAABMTdK6AAAAg1E4ldbyEzV03snPOelCAzkvuf9VoEjRgRQdyFqEUIqVxQZpVyg0DYIQAOAGjtfw3RV8Z6lUVC0U1cjFDh4TyFqFUUwgiwmiyAAKFH96cJ1MX1fSead6ppZO2blTpU6RrFsU6xHDBjZnPaOZhJE4nUEQAgBcxeEqvr6Ebyzjm8vUAIH1i2Ndw+iB1kqfFgGtw1hAk8Os2kWFNn7gPP/6HH+3SD1RwwfFsztbCmOTWNtw5s3/ATQVghAA4AcKp4JSvvSEuqqYq5xGt2QTU9jf+0uJoYyIampkxpjVenPpFW6hAc3YgGY//FZVA20oUdec5rO+VVqEsPvbCL9oy1qGIhG1hCAEALPjRFvL+Effq5+fUFuHsUmthNzRQpcor4RTZABNai1Mak0qFzeW8U+Oqj2WKIPihac6CmNaMgGBqAUEIQCYV2kt5RxRc4rUYIkeaS98NUFKuckOn9sERmktWFoL8c3+4qfH1Bl7lBd30kvdhQfbCU0fdwWPQBACgBltKuNvF6r5Z9R7WwufjRD7xGrWFwuVaEqqMCVVyD/D5+xT/vC1+vvewqPtBRG9Q19BEAKAicgqfXJMfeuAWivTs12E94dYwixa1/Sj4QlseIK0o5z/brfy533qn28TxycjDH0BQQgApuBU6L0i9S/71RQr/amvmN6S6TNk+jdjX46T1p7mz+9Q/nOY/tZfxOJSb8NQNAD4uQaV/n1Ibf+ZvPq0ujBN/HKcNEavKXhReku27y5pSLzQf4X8532qwrUuyK8hCAHAb6mcFnyvpn4mf3FKXTJKzB0tXbyMQf8sAr3UXdiVJa05rd6RKx+uQhh6C4ZGAcA/rS/hL+5UrBb6aJg4KN4w+XeZ1mEsb6z0ziF1yBfyG7eJU1JvuvfSq1cvbxTme2lpaW+++aY3nhlBCAD+pugCf3GnUnSB5twmTEgx/LgXI/pVJ2FoC3Z/vpJXwucNEq03s8Dnm2+++frrr71WnY8UFBRs3LjRS0+OIAQA/2F30et7lfePqK/0ED8f6VcX5HWOZDszpd9sVwbmystGiW3CbqKP6wedwpMnT3ovCP3oawIA5rbkhNp5sVxWR/snWf6nm1+lYKNgid4dLD7dUbhjhbzhDKYMPQY9QgAwvGIHz96mHq3mC9PEwYadDmyiZzoLXaPZfRvkWf3EyTc/ZQhXwpsIAAbGif51UO2zVO4Xy76eKPl9CjYaEs8Kxkuv71Vf36tqXYs/QI8QAIzqaDV/fLPiUmlzhtQhwhQReFGHCLY1Uxq/Vi6v428PxO3Ybgl6hABgPJzo34fU/ivkCSnC5vGmS8FG8cGUP07ae44/uVlRMWN4CxCEAGAwZ2r5nWvkD46oWzKk57sKZt66KNxCa++UjtXwRzcquPuM2xCEAGAky0+qfZbKA5sLW803HHpVoRJ9MVoqq+NPbkYUuglBCADGUCfT01uU3+5Ul46SZvQSJJy9fhQs0bJR0pFq/uw2RetaDAlfJQAwgENV/Lblsl2mvROl/sa5X6jPhEq0Ml3aWcFf/QpZeNMQhACgd/O/U4d+Ib/QTfhomKif7QP1JtxCq9Klxcf5vw/hmoqbg8snAEC/6hXK3qZsO8sLxkudI9ERvIHYIFo9RhycqySGUKbxb7LqM3inAECnjtfwASvkWpl2ZSEFm6pNGFs+WvzlZmVPpU6XzlRXVz/++OPx8fEdOnSYP38+08HWkOgRAoAerT3NJ2+UX+0lZnfG3+s3p28smzdYvCtP2ZklxQdf5QGfHFU/OeajmLy/Dbu/7c8+wZdeeomIjh07JgjCc88955syrg9BCAD6wone+Fade1D9vxGScfcR1NaEFGH/ebprvfzluKuc5HvGsCBfnfs7XnGJy7Jly/bv3x8SEkJEf/zjH+fNm+ejUq4NQQgAOlIr05RNykk735UlJoQgBd03vZew/zx/ZutVFpF2jGQdtRtqVhTl4nCoKIpalXEpjDkAgF6csvNBuXKQSAXjJKTgLWJE7w8Vd1bobqZw/Pjx06ZNq6urq6+vnzFjhtblECEIAUAndpTzgbnKw+2FD4aKgbroJxheqESLR+jurXzrrbfq6+uTkpL69OkzcODAiIgIrStCEAKADnx8VM1aL88bJL7QFSclT9JwCPRaNm/ePHPmzMrKysLCwn79+jVr1kzrihCEAKApTvTHveprX6n5Y6WxSbo7a4PHbdmyJTs7u7y8vKSkZOrUqXfffbfWFSEIAUA7DSpN3qisPKVuy5S6RCEFTWHGjBmRkZGpqal9+/Zt3rz5tGnTtK4Iq0YBQCNVDTRxvRwTxArGScE4FZmG1WpdsGCB1lX8DHqEAKCBYgcfnCv3jGGfDReRgqAtBCEA+Nr+8/yOFcqUDsJb/UUzb6sLOoEgBACf2lTGR62W/3K7gAWipvXaa6/Fx8drXcVPMCQBAL6z9IT69FZlUZo0PAE9QfP697//ffjwYa2r+AmCEAB85L0idcYedc0YqVcMUtDUzp07Fxsbq3UVP8HQBAD4wpx96p++UQvGiUhBk5swYQIR9ezZU+tCfoIeIQB43e92K7mn+ObxYmIoUlB79Yd21x/6yjevFdSpb1CnfpceWbZsGWPsm2++8U0BTeFmENpstkceeWTr1q2DBg2aP39+VFRU4/GNGzc+99xzR48ebdu27dtvvz1kyBDPlQoAxsOJXtihbCnjBeOk2CCtqwEiIhJCwqTYBJ+9lm9e6Fa4GYSzZ89OSUlZvHjxb3/72zlz5syaNavx+EMPPfTmm29OmDBh6dKlDz300KlTpzxXKgAYjKzSY5uUYgf/cpwUZtG6GvhRQErHgJSOWlehI24G4dKlS5cvXx4YGJidnZ2VlXUxCMPDwy9cuGC322tqaqxW61V/V1XV11577crjbdq0ueeeey476HA4OOcXN68CTTgcDlEULRacyTQjy7LD4QgOvtp243rVoNKU7ZY6hf7vDhdzkt2pdUG3zOFw4FykIUVR7Ha7w+EIDAxs+kaGISEhgnCD1TBuBmFJSUlKSgoRpaSklJaWXjw+f/78fv36PfHEE0S0e/fua/16Q0PDlQddLhfnl2+dxX/kXp3gEfgUNGe4j8Cp0GPbA2TOPxroChLJOIVfT+P7b6BPwf/wSzT9V274GDeD8GIvjXOuKD/tgDx16tSXX375+eeff+utt1555ZW8vLwrf1cQhNmzZzfxhRhjnPOwMAOMMvsxp9MZFhaGHqGGZFmWZdkoDaFWpgfWy7FB7MOhoiT41cQgY+xaY13QdO79MSGKYlhYWF1dXVhYmGe3tnczCBMSEoqLi9u3b19SUpKYmHjx+M6dOxcuXBgfHz916tRWrVp5pkYAMI5amTLXyc2CG1NQ62oAmsDN72lGRkZOTg7nPCcnJysri4gKCgqIqHv37u+9957dbv/www979OjhwUIBQP9qXDRmjdw6jH00DCkIhuHmV3XGjBn79u1LSkoqLCycPn06EaWlpRFRTk7OqlWrWrRosXjx4nfffdeTlQKAvtW46M41cpcoNm8wbqUNRuLm0GhkZOTKlSsvPdI45tuxY8etW7d6oC4AMBSHTBnr5M5RbO5AESEIxoLBCwC4VXYXpa+WO0Wy/wxCXxCMB0EIALekVqaMdXLHSPQFwahwr1EAcF+tTOPWym3D2X8HIwXBqNAjBAA31SuUtV5uFcbmDUIKgoEhCAHAHQ0q3bNBjgxg/8W8IBgchkYB4Ka5VLpvgxIkskVpuF7QAP76179qXcKtOnDggPeeHEEIADdH4fToRkXh9BlS0AhWrFjReMOTJso/w0MtdHucvrr50dHRQ4cO9dKTIwgB4CZwoqe2KBX1PHe0ZEEKGkFGRkZGRkbTH3+ihvddJufea4kK9F5R+oIvMgDchOe3K4er+LJRUpAnb3oMOtIqjGWmCP88qGpdiO8gCAGgqWbsUbac5SvTpVCMJfm1qT2Efx5UHLLWdfgKghAAmuTtQvXTY3xVuhQRoHUp4GUdItiQeOG/h83SKUQQAsCNvVuk/v2Amj9WbB6sdSngE7/rIfx1v9pgjihEEALADSw+rv7ha3XtnWJiqL5WEoL39I5lXaJowXemSEIEIQBcz7oSnr1NWZUutgtHCprLKz3Ev+5X3dlL3mgQhABwTTvL+cMF8pKRUvdopKDpDGvBgiVaXez/UYggBICrO1jFJ6yX3x8iDWyOFDSpF7oKb+5XtK7C6xCEAHAVp+z8zjXKX24XxyYhBc3rvjZC0QX69ryfdwoRhABwucp6Sl+j/Lab8GA7nCJMzSJQdmfhrf1+vmQG33IA+BmHTOPXyRNT2LNdcH4AerKjkHtKLa3Vug5vwhcdAH7iUumeDXKXSPanfriFGhARRQXS/W2FeX59cT2CEAB+wIme3KJIjP0HG+3CJZ7pJPy3SHX5bxQiCAHgB6/uVg7a+KLh2FwJfqZLFGsfTstO+m0S4vsOAEREcw+qS0/yVWNwQ224il93Fub6734UCEIAoKUn1FnfqqvTxRjTbEEHN2VCinC0mg7Y/PM6CgQhgNntquC/2qrkjhZbhWFmEK5OEuiJjsK//LRTiCAEMLWiCzxrnTx/qNQrBikI1/NEB+HTY6rdpXUdXoAgBDCv8joau0Z54zYxvSVSEG6gRQgNiRc+PeaHnUIEIYBJ1ck0Yb08OVV4tD3OA9Akv+wovFeEIAQAv6Bwui9f6RTJXuuFkwA01ZiWrNjhh0tm0AYAzOjZ7YpT4e8Mwu1j4CaIjB5LZe/6XacQQQhgOn/Zr24s5Z+OkCw4AcBN+mUH4ePv1Xr/2poJ7QDAXD4/rr59QF09RowM0LoUMKBkK+sVw5ae8KtOIYIQwER2lvNntikrRotJoVgmCm56LFWY/x2CEAAM6EQNvytP+e8gsScuGYRbMKGVsLuClzj8Z8kMghDAFKoaaNxaZVpPITMFrR5uSZBId7USPj6KIAQA42jcZXB0S/brzmjy4AGPtBcWfO8/o6NoFQD+79ntSoBAf7kdF0uAZwyKZ7Uy7T3nJ51CBCGAn5v9rbqrgn82QsJmu+ApjOihduxDf1kygyAE8GdLTqhzD6orRonYZRA86+F2wqKjquwXUYggBPBbX1Xyp7coy0aJibhYAjytbThrE8byzvjD6CiCEMA/nXbwieuVeYPF3rFIQfCK+9sKnxz1hy4hghDADzlkylqvPNdVmICLJcBr7mkt5J7yh9utoZEA+BtO9PgmpXs0e7EbGjh4UYsQ6hnDVhcbvlOIdgLgb6btVsrq+H+wswR43y/aCp8cM/w0IYIQwK8s+F5dfJx/PlIKQOMG77urlbD2tGp3aV3HrUFbAfAf287y3+5Qlo4SYwK1LgXMITqQBsez5SeNPTqKIATwEyft/J4NyoJhUtcoLBMF37m/jfDZcWOPjiIIAfyB3UWZ65SXuwvpLZGC4FPjk4WNpWqNkUdHEYQAhqdyeqhA6RPLnuuKFg2+FhFAdzRnq4y8dhTNBsDwXv1KOefk72CZKGhkUmvhcyOPjiIIAYzt46PqZ8f4UiwTBe1kJgvrS9Q6Wes63IWmA2BgX1fyF3Yoy0aJsUFalwImFhtEfWLZ2hKjjo4iCAGMqrSWJuYp/x0kdovGAhnQ2F2thCWGHR1FEAIYUr1CE9bLT3cSMnE3UdCBia2ElcVqgzH7hGhCAIb0xGalbTh7pQeaMOhCixDqFMm+NOauTGhFAMYz+1u10MbfHYw950FHMpKFFacM2SVEEAIYzOpi/o+D6orRYgg2nQc9yUxhuSe5EbuECEIAIym6wB/dKC9KE1ti03nQmU6RLEiib88ZLwoRhACGUe2iu9Yr/9tPHByPFAQ9Gp/EVpxCEAKAdyic7tsgj27JftkBzRZ0KiNFyDXgThRoUQDGMPNgkMzpz7fhPmqgX4Obs2M1/EytwTqFCEIAA/j4KM89E7AoTZLQZEHHJIHGtBS+MNroKFoVgN7tqeQvfUUf9XfgPmqgfxnJLNdoF1EgCAF0rbyOJuUpc/tTxzBF61oAbmxMkrCplNcb6tvqZhDabLaMjIzo6OjMzEybzXbxuCzLzzzzTFxc3B133FFSUuKhIgFMyqXSvfny5FQ2IQXLRMEYIgOoewzbVGak0VE3g3D27NkpKSmlpaXJyclz5sy5ePxvf/tbdXX1yZMnBw4c+Pvf/95DRQKY1LPblYgANqMXFsiAkYxpKaw9baTRUTeDcOnSpdnZ2YGBgdnZ2UuWLLl4/OOPP37ppZdCQkJmzJjxq1/9ykNFApjRu0XqxlK+YJgooDcIhjKmJVtdbKQeoZv3aCopKUlJSSGixn7hxeMnT55ctGhRWlpamzZt3n///av+rqqqTz/99JXHO3ToMGXKlMsOOhwOzrmqGumPC/9TXV1NRBaLRetCTGTXOeGVXYGrhtXzWn6BSJbl6urqgIAAresyNbvdzhhTFENNf2mhrYXO1QfvO1OdEur5OKyurpYkSRSbOkxitVpv+GA3g5Bzzhhr/MelX4vq6mrOeWFh4dy5c5944okdO3Zc+buMsXbt2l15PDEx8cpTrSRJhFOw1iw/0roQszhbz6bskOYNULrG/NBCGWP4CDTX+P7jU2iKEfF8Y2XA45Ge78M0NoSmB2FjVN3gMdytW6S2b99+1apV7du3/+6778aNG3fkyJHG4wkJCXv27GnRokVZWVm7du3sdvtlv+h0OsPDw51OZxNfyG63c87DwsLcKBI8pbKyMiIiAu3fN5wKDVspZyQL03r+NHMhy7LNZouLi9OwMKipqWGMWa1WrQsxgEVH1c+O8aWjPD+9XV5eHhMT0/QgbAo35wgzMjJycnI45zk5OVlZWURUUFBAROnp6R988IHT6Zw3b17fvn09WCiASfxmu9IihP2uJy5tAgMblSgUlBpmn143G9uMGTP27duXlJRUWFg4ffp0IkpLSyOiWbNm5efnN2/efMOGDe+++64nKwUwgX8dVLeU8Q+GYqNBMLbYIEqNYNvOGmPJjJtzhJGRkStXrrz0SOMQa3x8/Pr16z1QF4D5bCnjf9yrbM2QwjEIDcY3piVbe1od1sIAF/9g+AVAF0oc/L58Zf5QqW04eoPgD0YmCnklxugRIggBtOdS6b585dkuQnpLpCD4idubsSMX+PmmrozUEoIQQHu/2qokhLCXe6A9gv8IEGhgc1ZQaoAFM2h4ABp755C6s5znDMECGfA3IxOFDWcMMDqKIATQ0vZy/oevlaWjRCsWyIDfGZHANhhhmhBBCKCZsjq6b4OSM0RqhwUy4I96xLDzTl7s0HsWIggBtNGg0qQ8+alOwtgkpCD4J0aUliDk6350FEEIoI1ntylxQWwa7iADfs0Qo6NohAAa+PA7taCUz8cdZMDfjUhgeWdUnSchghDA13ZV8Jd2KctHixHYVQn8XdtwFiiyIxd0HYUIQgCfOltHd+cp/xoodohAbxBMYWg821SKIAQAIiKSVbovX57SgU1qjaYHZjGkBdtYhiAEACIienGXYpVoRi8D3IYYwFOGxrON6BECABF9fFRdVcw/SpMEjImCmTTeR/5YjX6zEEEI4Av7zvMXdihLRoqRWCAD5jNY39OECEIAr7M5aVKe8o8BYtcodAbBjIa2YJt0PE2IIATwLpXTgwXyXa3YvW3Q3MCkhuh7mhAtE8C7XtujOBX6U18skAHz6hjJHDI/ZddpFiIIAbxo+Ul14ff8k+GShKYGJsaIBscLm/U6OorWCeAth6v4k1uUz0eKcUFalwKgtaE6vpoQQQjgFXYXTcpTZvUT+8RigQwA3dGcbUUQApgHJ5qySRmWwKakookBEBF1j2Yltfy8U+s6rgatFMDz3vhWLXbwN2/HAhmAH4iM+sayHeV67BQiCAE8bMMZ/s9CdfEIMRA5CHCJgc3Z9nJV6yquAkEI4Emn7PzhAvnjNDExFFODAD8zoJmw7Sx6hAB+rV6hSXnKKz3EoS2QggCXG9Cc7a7gsv76hAhCAI/59VYlNYI92wXNCuAqIgMo2cr2ndddpxAtFsAz5h5Uv6rk/x2MiUGAaxrYnG3T33oZBCGAB+wo5zP3KktGiiGS1qUA6NiAZmy7/qYJEYQAt6q0lu7ZoOQMkRr3XQOAa7kDPUIA/+NS6f58+amOwp1JSEGAG2gfwRwuXuLQVxYyzn1akNPpDA8PdzqbencBu91eu/w/FgGBraX6+vqAgAABn8LV7CjndheNTPRuCqqq2tDQEBSEm5ZqyeVyEZHFYtG6EGPLK+HtIqiV1c0m44xNihs+URQ9ORlvgAkNKSk1IAC7emupwW63BAd79pvnH3aU87Uin95HCPByS1JV1eVwBISFefdl4LrU+nrGWEBgoNaFGJsq8v0qT23l5h/WTinYs/WQIYLQ0mtYKNq/puoqK4MjIvCH8GX2necPHJXzH5Rivb/vvCzL9TZbaFyct18IrkOtqWGMhVqtWhdibJEt+Zx9yiP93UwfR3m5Z+shzBECuMfmpLvylH8MELt4PwUB/Em/OPZVBVf1NEuIIAS4aSqnX3wpT2rF7m2DFgRwc6IDqVkwK7qgoyREMwa4aa/tUWSV/rcfJk0B3NEvju2uQBACGNaKk+rC7/mi4ZKIMVEAtyAIAQys6AJ/Yovy+UgxDhcyALjrtji2C0EIYEQ1Lpq4Xnmjn9gnFp1BAPf1imGFNt6gm20oEIQATcKJHt+kpCWwx1LRagBuSYhE7SN0tA0FmjRAk7zxrVrs4G/ejgUyAB7QL47t0s1NRw1wQT2A5taX8LkH1V1ZYiByEMAT+sXq6O7b6BEC3MCJGv5IgbwoTUwIwdQggGfc1kxHC0cRhADXU6/Q3RuUV3uJg+ORggAe0yWSHa/hdpfWdRARghDg+p7crHSOZNmd0VIAPEkSqFMk22/TRacQzRvgmt4uVPfb+H8GYWIQwPN6xbBvzukiCLFYBuDqtp3ls75RtmZKwWglAF7QK4bt1UcQokcIcBWltXRfvvL+UKlNGKYGAbyiVyzbW4kgBNAll0r35su/6SKMaYkUBPCWHtHsUBV36eD+MghCgMv9epsSG8he6o7WAeBFIRIlWdnhKu07hWjqAD8z77C6uYzPH4a9JQC8TifThAhCgJ/sLOev7VGWjxLDLVqXAmACCEIAfTlbR/fmK+8NllIj0BsE8AUEIYCOuFS6L19+vIMwPhkpCOAjvWLZN+e45kmIIAQgIvqfHUpEAJveEy0CwHdiAincwo7XaByFaPYAtOB7dV0J/3CoKKA3COBberiaEEEIZre7gr+4U1k2SowI0LoUAPPpGU3far1DL4IQTK28ju7eoMwbJHaKRGcQQAPdotn+8xrXgCAE83KpdM8GeXJ7lpWChgCgjW7R7IDWe1Cg/YN5NS6Q+X1vbC4BoJl24ay0ljtkLWtAEIJJzf9OXV/CFwzDAhkALYmMUiPYQU07hQhCMKPdFfzlXcpSLJAB0IGuURqPjiIIwXTO1tGkPGXuQCyQAdCFrtGsEEEI4DMX7yBzd2t8+QF0oWsU26/pFRTunAtsNltGRkZ0dHRmZqbNZrvyAQcOHAgNDb3l2gA87/kdSkQAe60XUhBAL7pF037D9Qhnz56dkpJSWlqanJw8Z86cy35aVVX16KOP1tbWeqI8AE/68Dt1A+4gA6AzSaGsTqZzTs0KcCcIly5dmp2dHRgYmJ2dvWTJkkt/pKrqo48+Om3aNA+VB+AxO8r5b3cqS7BABkB/ukSxA9qNjkpu/E5JSUlKSgoRNfYLL/3RrFmzUlNTJ02adJ1fV1X1rrvuuvJ4t27dnn/++csOOhwOzrksa3qNielVVVWpqmqxGHiPvrP1wqT84H/2cbbg8tWG8/VOluWqqipJcqfBgqfY7XbGmMvl0roQP5QaGrizRO0edOP3tqqqShAEUWzq5b/h4eE3fLA77Ypzzhhr/IeiKBeP5+XlbdiwYd26ddf/dcbYxIkTrzweHx8fHBx82UFFUTjnVx4HXwoODg4ODjZuELpU+uVm4emOlNXaQmTI/4Usy06nEw1BW7IsM8bwKXhDtxh26IIUHHzjSGo8HTU9CBvT6vrcCcKEhITi4uL27duXlJQkJiZePJ6Xl/fll19ePF0yxjZv3jxo0KAry3r44Yeb+FqyLHPOg4KC3KgTPCUwMDAoKMi4QfibzUrzEJreRzTuzKAsy42fgtaFmJrL5WKM4VPwht7N+JJiJSjoxvMWjQ2h6UHYFO7MEWZkZOTk5HDOc3JysrKyiKigoICI3njjDf4jIuKcX5mCAD72z4Pqrgr+4TADpyCA3+sazQ6c12yHXneCcMaMGfv27UtKSiosLJw+fToRpaWlebowAA8oKOV/2qssGyWGYnINQMdiAilYohKHNlHozukhMjJy5cqVlx7hVwT5lUcAfOyknT/wpbwwTWodht4ggN51imSHq6ilFpeg47Ji8E8OmbLWKa/0EIcnIAUBDKBjJDt8QZseFIIQ/BAnemyj0juWPdsF33AAY+gQwYqqEIQAHjJzr3rawf99BzYaBDCMDhGsSFjWNTEAABaBSURBVKMeIZYQgL9ZflL9z2F1Z6YYiBwEMI6OkXS4SpuXRo8Q/EqhjT+5RVk6UkwMxdQggJEkW9l5J7drcd8eBCH4j/NOmrBeeau/2C8OKQhgMIyofQQ7osXoKIIQ/ISs0j0b5Ltbswfa4lsNYEhaLRzFKQP8xPM7lGCR/tQXE4MARtUxgh3WYuEoghD8wftH1A1n+MI0CRsNAhhXhwgquqDB6yIIwfC2lPFpu5Xc0dhoEMDYOkaiRwhw847X8Hvz5Y/SpHbh6AwCGFtqBPu+mqs+j0IEIRhYjYuy1ivTeoojcB81AOMLkSguiJ20+zoJEYRgVCqnhwqUfrEsuzO+xgB+QpPL6nEGAaOa9pVS3cDfGYRlogD+Q5MrKHCLNTCkBd+rnx/nO7IkC/6WA/AjHSLYt+cwNApwI9vL+Ys7lRWjxZhArUsBAI9K1eLW2whCMJhTdn53njJ/qNQpEgtkAPxN2zA6Wu3rF0UQgpHYXZS5TnmpuzCmJVIQwA8lWVlFPa9XfPqiCEIwjMZlon1i2fNd8b0F8E8ioyQrO1Hj09FRnFDAMH63G8tEAfxf2zA6WuPTV8SqUTCGnCPq0pN8RyaWiQL4ubbh7Gg1J/Ld9AeCEAxgcxl/ZZfy5TgpGstEAfxd23B2DEOjAJf6vprfu0FeNFzqEoUFMgD+r20YHa1GEAL8yOak8WuV/9cHdxMFMIu24czHV1AgCEG/XCrds0Eel8ye7IgvKoBZtA5jJ+w+3YMC5xfQr2e2KiES+/NtWCYKYCIhEkUFsJJa3yUhghB06s/71D2V/OM0EZvOA5hN23Cf3l8GQQh69MUp/nahuny0aLVoXQoA+FybMJ8uHMXlE6A7uyv445vlNWOkpFB0BgHM6MdLCX0EPULQl5N2PjFPeW+w1CsGKQhgUhgaBfOqdlHmOmVqd2F8MlIQwLza+nZoFEEIeuFS6e48OS2B/aYLvpYApoahUTCpp7coQSL76+24WALA7GKDSOFkc/ro5RCEoAsz96rfnOcfp4kixkQBwLcLRxGEoL2Pj6rvHVG/GC3hYgkAaNTGh6OjuHwCNLapjP/PDmXDWKlFiNalAIBupFip2OGj10KPELR0qIrfu0H+aBh2lgCAn0kKZcUODI2Cvyuro3Frlb/cLo5MRAoCwM8khdIpu49eC0EI2rC7aPxaeUqq8FA7fAkB4HLJVlZsR48Q/JfC6aECpWcMm94L30AAuIokKzuFoVHwY09vUVwqf+cOXDIIAFfXPJiqG6he8cVrIQjB1/64V/3mHP90uCTh2wcA18CIEkN9NDqKUxH41HtF6vwjam46LhkEgBtIDvXRFRS4jhB8Z1Uxf22PsnGcFB+sdSkAoHvJVnbKzom8vqocQQg+squCT94o546W2kfgYgkAuLEkK53ySY8QQ6PgC0er+YT18rzB4u3NkIIA0CRJmCMEv1FWR6NXKzP7ihNS8H0DgKZKtvro5jI4MYF3VbvozjXylA7ClFR82QDgJvjs5jI4N4EXNah0d548oBl7tSe+aQBwc1LQIwSjUzk9+KUSFcj+ORAXzgPATbNaSGJ03vvb8yIIwVuytylVDfzDoaKA9TEA4BbfTBMiCMErZuxRvqrkS0ZKgegNAoC7kkLplPcXjuI6QvC8fxSqnx3jmzOkMNw+BgBuQbKV+WC9DIIQPGzh9+pf9qubxotxQVqXAgAGl2RlpzE0CsayvoS/vEtdPUZMsWJiEABuVbJPrqBAjxA8ZksZf6hAzh0tdY5ECgKAByRZ2SmH6u1XQY8QPGPvOX73BnlRmnRbHFIQADwjOZSKvd8jRBCCB3x3gY9fq/xzoDg8ASkIAB7TIoSV1XFvTxIiCOFWnajhI1crb9wm3N0aXycA8KRAkUIkqvLyNfU4c8EtKa+jMWuU57oID7fDdwkAPK95MDtb590+IU5e4L5zThqxSn40VfifbvgiAYBXxAdTWZ13XwLnL3DThQYas1oen8x+1wPfIgDwFvQIQadqZcpcJ/dvxmb1wy3UAMCLmgfTWfQIQW/qFcpaL7cLZ29jWwkA8DL0CEF3GrcYjAti8waLuFQCALxNpz1Cm82WkZERHR2dmZlps9ku/dGSJUu6du0aGRk5ZMiQI0eOeKhI0AuF08MFSoDIPhwqIgYBwAfiQ1hZrf56hLNnz05JSSktLU1OTp4zZ87F48ePH588eXJOTk5paWlGRsZjjz3muTpBewqnB79UamX+yXBRwlACAPiETnuES5cuzc7ODgwMzM7OXrJkycXjx44de+CBB2677bbg4ODJkycXFRV5rk7QmMJp8kbF5uT/N0IKQAoCgK/4IAgZv/mb11it1oqKiuDg4Lq6uubNm1dXV1/2AEVRsrOzBUGYO3fuZT9yOp1Wq3XAgAFXPm2fPn1ee+21yw46HA7OudVqvdkiwYMqKs/NONbibIP48YD6INHrW6LAlWRZrqqqio2N1boQU7Pb7Yyx0NBQrQsxF6dCrb8IK8mqaZyNqaysjIqKEsWmrtSLiIi44YPd2X2Cc84Ya/yHoiiX/XTdunVTp05NT0+fOXPmVX+dMTZt2rQrj0dHR1/5DeOcc87xzdOQyukPxyxlDZZlIyhECtG6HJOSZdnlcqEhaEtVVQSh74UShUjUIIVGBxIR1dbWhoaGNj0IBeHGQ1juBGFCQkJxcXH79u1LSkoSExMvHuecT506dfv27Z9++mlqauq1fp0xNmbMmCa+lsvl4pwHBga6USfcOk70663K8Tp1dboQFYL95jUjimJAQAAagrYaGhoYY/gUfK95sFylBrQIZETU2BCaHoRN4c5sT0ZGRk5ODuc8JycnKyuLiAoKCoho06ZNK1asyM3NTUhIsNvtdrv3N88Ab1I5PbVFKbTxzwbWWhGCAKCReC9PE7rTI5wxY8aDDz6YlJTUu3fvBQsWEFFaWhrnvKCgoKioKCoq6uIj3ZiABJ1oTMEjF/jKdKn+Aj5HANDMj9fUe+uaLXcWy9wKp9MZHh7udDZ1Uw273c45DwsL82pVcJnGEdH9Nr4qXQqzUGVlZUREhMWCXqFmZFm22WxxcXFaF2JqNTU1jDGs3fO957YrbcPZs10EIiovL4+JifHs0Kg7PULwbwqnX25WTtTwNWOkUHxBAEBrzYO9e009rgiDn1E4TdmkHK/huaORggCgC/Eh+psjBH8lq/TwRsXm5KvTpWB8NQBAH5oHs7N1qveeH2c7+EGDSvfnKy6VLx8lBWJXCQDQDW/fXAZBCEREToXuy1dUTotHIAUBQF+8vUk95giBHDKNWyuHSrRkpIgUBAC9aRbMKuq9eIUDgtDsLjTQ6NVyQgibPxR7SgCAHgUIFCKRramX3d00nPlMraKehq+Sb4tj84chBQFAv+K9uU89Tn7mddrBh3whj0tib/XHJrsAoGteXS+DxTImdbyGj16tPNVJeLEb/hgCAL378Zp6r/zRjiA0o2/O8XFrlT/2ER7vgBQEAAPw6jX1CELT2VzG794g/2ugOKk1UhAAjCEuiFXUe2uOEEFoLitOqk9sUT5Ok0YkYFoQAAwjzIIeIXhCzhF1+lfKynSpbyxSEACMxGohh8tbT44gNItZ36rvHlY3jpPaRyAFAcBgQiWyy956cgSh/1M4Pbdd2VLGt2RILUK0rgYA4OZZLczu8tZ9txGEfs6p0CMblYo6vnG8FBGgdTUAAG4JlcjhtR4h1g36s3NOGrlalhituRMpCAAGZrWQ3WtzhAhCv3W8hg/KlXvHsAXDxAB8zgBgZOgRwk3bWc4H5SrPdxX+PkAUsDgGAAzOqz1CzBH6oc+Pq89sU94fIo1NQgYCgD+wSszuwgX10DRz9qlzD6rr7pR6RCMFAcBPWC1eHBpFEPqPBpV+tUX55jzfliEmhiIFAcB/WARiRE7FK0+OIPQT55w0KU+ODmSbxkuh+FQBwO9YLd66ph6LZfxBoY3fvlwe0IwtHiEiBQHAL1ktzOGdaUKcNQ1vzWn+6EZ5dj9xcir+rAEAv9V4BUWQF54ZQWhgnGj2t+o/D6rLR0n9m2FSEAD8WeMVFDFeeGYEoVHVK/TkZqWwim/LEJOtSEEA8HM/3Hdb9PwzYzDNkI7V8P7LZUmgbRkSUhAAzMBqIS9dSoggNJ41p/mAFfIv2go5Q8RAL/xxBACgQ6ES89KWhBgaNRKV05++UecdVpeOlAY2R0cQAEzEe5dPIAgNo9pFkzcqlfV8Vxa2FQQA0/HeJvUYGjWG3RW85xK5bRjlj0UKAoAZWb22ST16hAYw77A6/Stl7h3iPa3xhwsAmFSohdU04IJ687E56ZeblVN2vjNLah2GSUEAMC+rRGccXnlm9DD0a3cF77dcbhlKWzKQggBgdqEWqsXQqHkonGZ/q75dqMwbJGam4I8VAACySlSDyydMotjBHy5QBKKvJkgtsZsSAAARNd50W1a98czobejLp8fUvsvkMS2FvLFIQQCAnzTea9Qb0CPUi6oGyt6mfF3JV6VLfWIRgQAAPxPqtcsn0CPUhXUlvPvnclQA7ZmAFAQAuArvXVCPHqHGalz08i5ldTF/f6g4IgERCABwdegR+qf1Jbz7Etml0r5JElIQAOA6rBbmpd0n0CPURlUDvbRTWV/C/zNITG+JCAQAuAGrRA70CP3GkhNq18/lAJH2TZKQggAATSEJJDJyqp4/Z6JH6FOn7Pw329XvLvBPh4t3YB8lAICbYZWoVvH8mRM9Qh9xqfTX/WqfZXK/WPbNXRJSEADgZlktzBujo+gR+sLGUp69TUkMpe2ZUrtwRCAAgDtCJXLIGBo1mtMO/vIudetZ/lZ/4a5W6H8DALjPasHQqKHUyTRzr9prqdwunA7dLSEFAQBuEXqEhqFyWnhUfXW32r8Z250ltcIOSgAAnoA5QmPIK+FTdysWgT4ZLg7EihgAAM/x0tAogtBjdlfwV79STtppZl/h7tYCMhAAwLMwNKpf+8/zGXvUryr5qz2FxzsIFswGAgB4gdVCDvQI9Wb/ef76XnVzmfpyD3HRcDFI1LogAAD/ZbVQbR2CUDd2V/BZ36rbz6ovdhffH2oJxRsJAOBloRI7o3j+aXH+vjmcaH0J//M+5cgFerGb8NEwSwjeQgAAn7BaqBZzhBpyKvTJMfXN/Son+m034YG2mAsEAPApq4Q5Qo2UOPg7h9X/HlZ7xrDZt4npLRlWhAIA+F6oBatGfUvllHeGv3NI3VSq/qKtsHG81CECCQgAoBmrxGoxR+gbx2r4h9+p7x/hzYLoyY7CgmFYCwMAoD0reoTeds5Ji4+rH32vfneB399W+GK02C0aXUAAAL3ABfXeUlFPK06qi4+rO8r5mCRhancxvSXDQhgAAL3BLdY87FAVX1nMvzilfnuOj24pTOkgfD5SwLUQAAC6hR6hB1TU05dn1LwzfH0J55zGJrGXu4vDExjuCAMAoH9e2n3CzRFAm82WkZERHR2dmZlps9lueFwrnKjoAv/wO/VXW5Uui+XUz1wffc+7RrGV6eKJ+6V/3SGOTUIKAgAYQ6ikp415Z8+enZKSUlpampycPGfOnBse9xmHTHsq+fzv1Bd3KsNXytEfusasUdac5p0i2UdpYuXDlhWjxWe7CJ0jsQoGAMBgJIEkxus83SlknHM3fq1Dhw7Lly/v2LHj4cOHs7KyioqKrn/8IqfTGR4e7nQ6m/hCdrudcx4WFvazgy6qauAV9VRaS2fr+Ek7P2mnEzX8u2o6V887RbLOUaxbFOsVy3rHsNggN/5/8DOVlZUREREWi0XrQsxLlmWbzRYXF6d1IaZWU1PDGLNarVoXYmoxHzoPTpKah3pyKM/NOcKSkpKUlBQiauz/3fD4pRRF6dy585XH+/fvf2Un0uFw3LUr9pzLSUQyZw6FHDKzMB4RQDEBanyQGhvIk0LUXlaeGae2DlUSgvnPdgK0U6Xdvf8i/OT8+fMulwtBqCFZli9cuMBwUyNN2e12xlh9fb3WhZja6Bh2/rwqNnkPisjISEm6QdK5GYSc88Y2yTlXFOWGxy8lCMJHH3105fGwsLCIiIjLDoqi+Nc+roSoACKSGFklCrVQwA9ZJxJhfs8XXC4XeoTakmVZVdUrGwj4kiAI6BFq7i89y2NiYkSxqSf/pjzSzSBMSEgoLi5u3759SUlJYmLiDY9fijHWu3fvJr6QxWLpE8fDwnAK1pLlR1oXYl6MMXwEmrNYLI0fhNaFmFpjQ2h6EDaFm4tlMjIycnJyOOc5OTlZWVlEVFBQcNXjAAAAeuZmEM6YMWPfvn1JSUmFhYXTp08norS0tKseBwAA0DM3gzAyMnLlypWnT59esWJF47xF4+rTK4/foi1btuTl5d3688CtmDdvXllZmdZVmFp5efk777yjdRVml5+fv2nTJq2rMLu5c+eeP3/es8+p91tq7tq1a/PmzVpXYXYLFy4sLy/XugpTq6ysXLhwodZVmN2WLVt27typdRVm98EHH1RVVXn2OfUehAAAAF6FIAQAAFNDEAIAgKm5eYs1tzU0NLRp02b48OFNfPzhw4ddLle3bt28WhVcX15eXr9+/XA1t4aqq6t37tw5atQorQsxtQMHDoii2KlTJ60LMbU1a9YMHjw4NDS0iY+fOXNmcnLy9R/j6yAkor179+7fv9/HLwoAACaUkZERFRV1/cdoEIQAAAD6gTlCAAAwNQQhAACYGoIQAABMTV9BaLPZMjIyoqOjMzMzbTbbDY+DZ13nfR44cCD70dNPP61VheahKErHjh0vPYJW4EtXvv+EVuBzS5Ys6dq1a2Rk5JAhQ44cOXLxuMfbgr6CcPbs2Y07+iYnJ1+6Se+1joNnXet95pwXFRWVlpbW1NTU1NT87W9/07BIM/j73/8+cODAoqKiSw+iFfjMVd9/tAIfO378+OTJk3NyckpLSzMyMh577LGLP/J8W+B6kpqaeujQIc75oUOHUlNTb3gcPOta73NJSUl4eHjv3r2tVmtWVtbZs2e1q9EU8vPzc3NzL2ueaAU+c9X3H63Ax/Ly8p566qnGf5eXl8fExFz8kcfbgr4un7BarRUVFcHBwXV1dc2bN6+urr7+cfCsa73Pe/bseemll958883k5OQXXnihoaFh0aJF2pZqBoz9rHmiFfjYZe8/WoFWFEXJzs4WBGHu3LmNRzzeFvQVhKGhoefOnQsKCqqtrY2Li3M4HNc/Dp7VlPe5tLS0S5cuHt8GBa502YkYrcDHLnv/L4VW4DPr1q2bOnVqenr6zJkzJUlqPOjxtqCvOcKEhITi4mIiKikpSUxMvOFx8Kxrvc979uzZtm1b478DAgICAwO1qc/c0Aq0hVbgY5zzl19++fXXX//000/feOONiylIXmgL+grCjIyMnJwcznlOTk5WVhYRFRQUXPU4eMO13n+HwzFx4sRDhw41NDS8/vrrEyZM0LhQk0Er0BZagSY2bdq0YsWK3NzchIQEu91ut9vJe23h1qcZPchms40dOzYxMTEjI6Oqqopz3ljhlcfBG671/quq+o9//KNt27axsbGPPPLIhQsXtK7UFC42T7QCTVz2/qMV+Ngf/vCHK9PKS21BX3OEAAAAPqavoVEAAAAfQxACAICpIQgBAMDUEIQAAGBqCEIAADA1BCEAAJgaghAAAEwNQQgAAKaGIAQAAFNDEAIAgKkhCAEAwNQQhAAAYGoIQgAAMDUEIQAAmNr/ByrJfdDoODxCAAAAAElFTkSuQmCC"
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra, Statistics\n",
    "using Distributions, Plots, QuantEcon, Interpolations, Parameters\n",
    "\n",
    "gr(fmt=:png);\n",
    "\n",
    "w_max = 2\n",
    "x = range(0,  w_max, length = 200)\n",
    "\n",
    "G = Beta(3, 1.6)\n",
    "F = Beta(1, 1)\n",
    "plot(x, pdf.(G, x/w_max)/w_max, label=\"g\")\n",
    "plot!(x, pdf.(F, x/w_max)/w_max, label=\"f\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='looking-forward'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looking Forward\n",
    "\n",
    "What kind of optimal policy might result from [(3)](#equation-odu-mvf) and the parameterization specified above?\n",
    "\n",
    "Intuitively, if we accept at $ w_a $ and $ w_a \\leq w_b $, then — all other things being given — we should also accept at $ w_b $\n",
    "\n",
    "This suggests a policy of accepting whenever $ w $ exceeds some threshold value $ \\bar w $\n",
    "\n",
    "But $ \\bar w $ should depend on $ \\pi $ — in fact it should be decreasing in $ \\pi $ because\n",
    "\n",
    "- $ f $ is a less attractive offer distribution than $ g $  \n",
    "- larger $ \\pi $ means more weight on $ f $ and less on $ g $  \n",
    "\n",
    "\n",
    "Thus larger $ \\pi $ depresses the worker’s assessment of her future prospects, and relatively low current offers become more attractive\n",
    "\n",
    "**Summary:**  We conjecture that the optimal policy is of the form\n",
    "$ \\mathbb 1\\{w \\geq \\bar w(\\pi) \\} $ for some decreasing function\n",
    "$ \\bar w $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take 1: Solution by VFI\n",
    "\n",
    "Let’s set about solving the model and see how our results match with our intuition\n",
    "\n",
    "We begin by solving via value function iteration (VFI), which is natural but ultimately turns out to be second best\n",
    "\n",
    "The code is as follows\n",
    "\n",
    "\n",
    "<a id='odu-vfi-code'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "res_wage_operator (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# use key word argment\n",
    "function SearchProblem(;β = 0.95, c = 0.6, F_a = 1, F_b = 1,\n",
    "                       G_a = 3, G_b = 1.2, w_max = 2.0,\n",
    "                       w_grid_size = 40, π_grid_size = 40)\n",
    "\n",
    "    F = Beta(F_a, F_b)\n",
    "    G = Beta(G_a, G_b)\n",
    "\n",
    "    # scaled pdfs\n",
    "    f(x) = pdf.(F, x/w_max)/w_max\n",
    "    g(x) = pdf.(G, x/w_max)/w_max\n",
    "\n",
    "    π_min = 1e-3  # avoids instability\n",
    "    π_max = 1 - π_min\n",
    "\n",
    "    w_grid = range(0,  w_max, length = w_grid_size)\n",
    "    π_grid = range(π_min,  π_max, length = π_grid_size)\n",
    "\n",
    "    nodes, weights = qnwlege(21, 0.0, w_max)\n",
    "\n",
    "    return (β = β, c = c, F = F, G = G, f = f,\n",
    "            g = g, n_w = w_grid_size, w_max = w_max,\n",
    "            w_grid = w_grid, n_π = π_grid_size, π_min = π_min,\n",
    "            π_max = π_max, π_grid = π_grid, quad_nodes = nodes,\n",
    "            quad_weights = weights)\n",
    "end\n",
    "\n",
    "function q(sp, w, π_val)\n",
    "    new_π = 1.0 / (1 + ((1 - π_val) * sp.g(w)) / (π_val * sp.f(w)))\n",
    "\n",
    "    # Return new_π when in [π_min, π_max] and else end points\n",
    "    return clamp(new_π, sp.π_min, sp.π_max)\n",
    "end\n",
    "\n",
    "function T!(sp, v, out;\n",
    "                           ret_policy = false)\n",
    "    # simplify names\n",
    "    @unpack f, g, β, c = sp\n",
    "    nodes, weights = sp.quad_nodes, sp.quad_weights\n",
    "\n",
    "    vf = extrapolate(interpolate((sp.w_grid, sp.π_grid), v,\n",
    "                    Gridded(Linear())), Flat())\n",
    "\n",
    "    # set up quadrature nodes/weights\n",
    "    # q_nodes, q_weights = qnwlege(21, 0.0, sp.w_max)\n",
    "\n",
    "    for (w_i, w) in enumerate(sp.w_grid)\n",
    "        # calculate v1\n",
    "        v1 = w / (1 - β)\n",
    "\n",
    "        for (π_j, _π) in enumerate(sp.π_grid)\n",
    "            # calculate v2\n",
    "            integrand(m) = [vf(m[i], q.(Ref(sp), m[i], _π)) *\n",
    "                        (_π * f(m[i]) + (1 - _π) * g(m[i])) for i in 1:length(m)]\n",
    "            integral = do_quad(integrand, nodes, weights)\n",
    "            # integral = do_quad(integrand, q_nodes, q_weights)\n",
    "            v2 = c + β * integral\n",
    "\n",
    "            # return policy if asked for, otherwise return max of values\n",
    "            out[w_i, π_j] = ret_policy ? v1 > v2 : max(v1, v2)\n",
    "        end\n",
    "    end\n",
    "    return out\n",
    "end\n",
    "\n",
    "function T(sp, v;\n",
    "                          ret_policy = false)\n",
    "    out_type = ret_policy ? Bool : Float64\n",
    "    out = zeros(out_type, sp.n_w, sp.n_π)\n",
    "    T!(sp, v, out, ret_policy=ret_policy)\n",
    "end\n",
    "\n",
    "\n",
    "get_greedy!(sp, v, out) = T!(sp, v, out, ret_policy = true)\n",
    "\n",
    "get_greedy(sp, v) = T(sp, v, ret_policy = true)\n",
    "\n",
    "function res_wage_operator!(sp, ϕ, out)\n",
    "    # simplify name\n",
    "    @unpack f, g, β, c = sp\n",
    "\n",
    "    # Construct interpolator over π_grid, given ϕ\n",
    "    ϕ_f = LinearInterpolation(sp.π_grid, ϕ, extrapolation_bc = Line())\n",
    "\n",
    "    # set up quadrature nodes/weights\n",
    "    q_nodes, q_weights = qnwlege(7, 0.0, sp.w_max)\n",
    "\n",
    "    for (i, _π) in enumerate(sp.π_grid)\n",
    "        integrand(x) = max.(x, ϕ_f.(q.(Ref(sp), x, _π))) .* (_π * f(x) + (1 - _π) * g(x))\n",
    "        integral = do_quad(integrand, q_nodes, q_weights)\n",
    "        out[i] = (1 - β) * c + β * integral\n",
    "    end\n",
    "end\n",
    "\n",
    "function res_wage_operator(sp, ϕ)\n",
    "    out = similar(ϕ)\n",
    "    res_wage_operator!(sp, ϕ, out)\n",
    "    return out\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The type `SearchProblem` is used to store parameters and methods needed to compute optimal actions\n",
    "\n",
    "The Bellman operator is implemented as the method `T()`, while `get_greedy()`\n",
    "computes an approximate optimal policy from a guess `v` of the value function\n",
    "\n",
    "We will omit a detailed discussion of the code because there is a more efficient solution method\n",
    "\n",
    "These ideas are implemented in the `.res_wage_operator()` method\n",
    "\n",
    "Before explaining it let’s look at solutions computed from value function iteration\n",
    "\n",
    "Here’s the value function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute iterate 10 with error 0.19801710153283736\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute iterate 20 with error 0.007608221868107279\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute iterate 30 with error 0.0002901698734376623\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converged in 34 steps\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAGQCAIAAADZeHv+AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de1RU16E/8H3mcQaYYWYAjTIEuYmJQTQqJrSpDxCaEkt4aK+/1d7CXZGQ/JblxwopBqUm0VSJiiap1NrbBEEbVlKxiVPz+9kQn7BE2pSQlYcWMLeKD8Qnj+EMMI8z+/fHVMpVZxzIzJzZ4/ez8sdkz5zdb1L5Zp/NmXM4SikBAAh4MqkDAAB4BG0FAGxAWwEAG9BWAMAGtBUAsAFtBQBsQFsBgM+dPHlSrVY7X/f29mZlZUVGRmZnZ/f29no+if/aav/+/TNnztTr9cnJyadPnx4ZH3d0AGBCf3//8uXLBwcHnX9bUVERFxfX3d09ZcqULVu2eD6Pn9rq/PnzeXl5VVVV3d3d2dnZ+fn5I2+NOzoABD5K6fLly8vKykZGjEZjUVGRSqUqKirat2+f51Nx/rmWvaGh4b333quqqiKEXLt2bfr06devX3e+9cgjj+zfvz8+Pr69vT0nJ6ejo8MPeQDAPzZv3nz9+vU33niD4/7ZNhqN5tq1a6GhoUNDQ5MmTTKZTB5O5ae2GiGKYlFRkUwm27Fjh3PEffT8n/6v9z8YQ/sCwLeRlvy9jw833eGNvj5itbo5sP2bb8oqKojsn6drGzZsePTRR48dO/bLX/7y0KFDSqVypK3UavWNGzdCQkIGBwcnTpxoNps9zKYY2z/Kt3P48OFVq1alp6eXl5ePDFJKOY5zvhBF8ZZD3v9gn9XmcH4AAHztk6PNd37jxAnS2enmwBtnznT8/e+b3njD+bfR0dGEkCNHjjQ2NvI87xzkOO748eMGg+HChQsPP/xwV1dXTEyM59n81FaU0jVr1pw4cWLPnj3Tpk0b/dZdo3Mc53j6af/kBLiXaT7+eNDhGPfhkXr9kiVLRo+Ul5ePLE1G1lZZWVk1NTUbN26sqanJycnxfH4/7bI3NzcbjcaPPvrIYDAIgiAIAiGkoaGB3IxOKR1rdABg0dq1a7/66qvY2NhTp0698sornh/op7VVQ0NDR0dHRETEyAilNDU1lVK6du3a3Nzc2NjYuXPn1tbW+icPAPjZyBa5Xq8/cODAOGbw9y77WKl4udXmeDA8TOogAMHvrDBICOe448nggQPu961OnDmzqrHxxGef+Sgb8fMu+/hwHHco9wmpUwAEv5k7G4fF8e9b+VrgtxVHCH1Qj7UVgM/JAvt374H+PUFcugDgN7LA/nkL/LUVIYT0DtukjgAQ/AJ7E5uBtuIopVPfOSp1DIDgZ7bZA/l0JtDbSibjOI7rqXpK6iAAwU//vz8xDdqlTuFSoO9bAYDfBO6yihCCtgIAVgT6maDT4ZPXpY4AEPzsYkBvswd6W1GOUEorDv+31EEAgt+w7daboASUgG8r57XsO3EtO4DPRSz4pF/ALjsAwLcT6G0lk+FOfAB+opAH9A9boJ8JEkIopRU1/5A6BUDwG7Jg3+pbcH4V4Cpxd0NoAPAKR0D/SjDg24rjOI7j1r82XeogAMGv9r0Lw1bcMWa8OOxbAfiLzM0tYyJ5Yg93d/BAmK/vOBPobQUAAYGzEseAuw/QQUJ8eybJQFs5HPT/vPiV1CkAgp9pIHAvtiIstJWM47hJM/RSxwAIfgpVN4d9q3HjOI6Tkey8KVIHAQh+u7d9Yw7g5VWgXx0KAOCEtgIANgT6mSAlxCHS5U8dlzoIQPDruWaROoI7DLQVJ+N+tGGW1EEAgt9/Pfc3oSdwvzcS6G1FCOE4MmWmTuoUAMFPrgzoraGADkcIkSnkXIA/khEgWCiUAf3dkUBvKwAAp4A/E3TI7FZH0UP/T+ocAPcEmdzLK5j6+vqSkpJLly4ZDIZt27alp6dfvXq1pKTkyJEjcrk8IyPjzTffDA93+w3EmwK9rThOJlfJ81t+KHUQgOC3N7tBuDjkxQkdDkdubu7evXtTU1P37duXn5/f1dX17LPPJiYmnjt3zuFwrFu37rXXXnvzzTc9mS3QzwQ5mSyAz6MBgopMznl338put9fW1qalpZnNZpVKpdfrCSGNjY0vvfQSz/MhISFlZWUffvihh7MF+toKANjF83xGRoYgCFqtluO4pqYmQsjjjz++efPm1atXW63W119/vbu728PZGGgrSun1v/dLnQIg+Inf7k7H1270vPvuu87X6enpkydPdr7WaDSCIFRWVhYXF7e0tOzevbuwsDA2Nva+++574YUXIiMjPZw/4NuKyhxW2rjulNQ5AILf4FXLt9kd6us3HT582Pk6KSlp8uTJnZ2dO3bs2Lp1q1qtLigo2LhxIyFErVYbjUae5wkhjY2NjzzyiIfzB3xbOTgZL5+/40dS5wAIfo3PfTB40Tzuwx9+8N9G1lZOBoOhuro6MzMzOTm5rq4uMTGRELJq1SpKaWVlpdlsLisrKy0t9XD+QN9lBwB28TxvNBpLSkqioqL27NlTVVVFCNm6dWtPT09sbOz3v//9Z5999kc/8nQtEuhrK5lM5uubPQOAk0whI97+HXxKSkpra+vokaioqP37949jqkBvK04mp6Lj/IF2qYMABD/rgNXlFQwqnmjdXsMZFkbIvf0UCc4hJyLp+zxw72cIEDTEQcJRF7tDaiuxuX2KRNggkd3bT5HgCJHx/Kyf5kkdBCD4mc5tMV+6JnUKlwJ9l12hUMnkcqlTANwT5CpeJgvcFUygtxUAgFPg9ugIu8XyVe0eqVMABD/z5cA9DSSB31YKGiKXh+jV86QOAhD8+lSX7bZeqVO4FOhtJZMp5MrQKXOypQ4CEPyu/KPZbDZJncIl7FsBABvQVgDAhkA/EySEWgf7/vr+S1LHAAh+pusdvEwpdQqXAr2tOEoUCu1DD7wsdRCA4NcxuME+dE7qFC4FelsRQmQy1YRJqVKnAAh+Z//7vwK5rQJ930opD5HLQ6ROAXBPkCvUCrlK6hQuBXpbyWQKjmNgAQgQBGQyZSB/8yZwk41CLcNXrZbrVmuPzdpjtfbYrL12a7/N1ifaBuw2wW4fsNsHbLYB0W4WRW8+XwjgnmIXB6O090udwqVAbysZcQgDp48f/K6KjwrhJ4QoI1R8ZIgyQqPUq0JjlFqtUq7hFeG8UqtUaBRytUIeJnVkAFb95eRLwtBnUqdwKdDbKlQZOjnyiex5R6QOAhD8QlX3iaJHj02WRKDvW3GE8/UNCQGACYG+tgLXqMXm0WMW7eKg6LD4Og0EAYs9oB/cibaSjOgY7r7R1Cd0DFquDA1fttn77XbBYuunRLTYTIQQuzgsOoadH7aJQ7c1Dhei1HnyP6SUh8plgftraQgcQ7aeuIkPSZ3CJbSVBLpvHD955tddN45P1s64X/foJNVE7YTvhijCVQp1iEIr4+QhSi0hRCFTKeWhzkOU8lCFjJc0NQS/A3//5XXrp1KncAlt5W89A6eOfP6fT09/9T/nvBWm1EsdB8AzOgUJcfsL94sqrz/d6xZoK7+i1HGwZdmSGeWzDDlSZwG4lULGK+UulvAKO7ENujuYsxBybz/zJshQIgpD3agq94Zs/ZQ6hu0mBxUtdsE54vxJGLaZqI9/JO5l/cOXtOERUqdwCW0FvuWgdsFybcBy3Wy9Zrb2Dtn6hqx9Q7Y+i91ksQvDNtOwrd/qGLaJg0O2AZs4bHNYQhXhMk4eptRyHBemCCeEhCm1HOE4QkKV4TJc0eIz182dD034N6lTuIS2Am86c6P5ykBHz+DZ/sGLvUMX+4e7zdbecFWUlo/SqaK0fIRGqdMptYbQCWrlg6EKTZgyPEwRrpKHqOShamW4Uqbi8SV26fzfMzWXZF9KncIltBV40x+/KHo06rsPaB+ZGDFzQmjMhNDJOj6K4wL9ImRwilDdRyJjpE7hEtoKvOynj7wYFRotdQoYD4VMySsD99I8tBXAmIlUFKx9ZpvJbDeZbaYhuzBkEyyOYYs4RAixi1aLOCx1xvG4ZD47675p3p2zvr6+pKTk0qVLBoNh27Zt6enpjY2NxcXF//jHP6ZOnfrrX/86OTnZw6nQVgCe+vCbHY0XPuy39g2JgzqlLlyp1SjDw5XhGoVGrQgLlYeGykM4QsLkyhCeyd03hV13/wMTvTihw+HIzc3du3dvamrqvn378vPzu7q68vLy3nrrrSVLlhiNxry8vPPnz3saz4vJAILb9aGuZXFL/33Kv4crA/dGBd/GN6ZvrPdbvTih3W6vra1NS0sTBEGlUun1ekKIVqvt7+8XBGFgYECj0Xg+G9oKYAz0Sn2wVpUv8DyfkZEhCIJWq+U4rqmpiRDy+9//Pikp6fnnnyeEtLS0eD4bflkDAN7hcFDzTQ6HY2Rco9EIglBeXl5cXEwIWb169apVqy5dulRaWlpWVub5/GgrAPCOv332xX03HT9+nBDS2dlZWlpKCFGr1QUFBW1tbYSQTz/99Oc//3l0dPTq1as//XQMX6JGWwGAdzzxncSRtVVKSgohxGAwVFdXNzY2Ukrr6uoSExMJIbNmzaqurhYE4d133509e7bn86OtAMBXeJ43Go0lJSVRUVF79uypqqoihNTU1Pz5z3+Ojo7+4IMPdu7c6fls2GUHAB9KSUlpbW0dPRIfH3/ixIlxTIW1FQCwAW0FAGzwX1uJohgfH3/L4Pz587mbVqxY4bcwAMAcP+1bVVZWvv/++x0dHaMHKaXt7e0XL17U6XSEEIUCm2gA4JKf1lazZs169dVXbxm8cuWK1WrNycmJjo7Oy8szmUz+CQMALPJTW6WmpmZmZt4yePny5aSkpJ07d547d06n0zmvc70FHXVFLAD4miiKdxynShVVad39xauD+SkSc+bMOXr0qPP15s2bZ8yYcftnOB//8wPAaDLZnVcwdpXFoXB39mMPMfv6KRJS/k7w888/b25udr7meV6lutNtwNBWAH4UyOsDadqqoaGBEGI2m5cuXdrW1ma1Wjds2LBkyRJJwgAAE6Rpq9TUVELIggUL1q1bl5WVFRMT09vbW1FRIUkYAGCCX/etKKWjX3AcV1hYWFhY6M8MAMAoXMsOAGxAWwEAG9BWAMAGtBUAsAFtBQBsQFsBABvQVgDABrQVALABbQUAbEBbAQAb0FYAwAa0FQCwAW0FAGxAWwEAG9BWAMAGtBUAsAGP8AOAuzOJimExzM0H+hwqBwneZ94AACusxD7kGHTzAYvDEszPvAEA8BzaCgB8qL6+PiEhQa/XJyQkHDx4kLuN51PhTBAAfMXhcOTm5u7duzc1NXXfvn35+fkDAwMj727ZssVqtXo+G9oKAHzFbrfX1tampaUJgqBSqfR6vUajcb719ddfNzc319fXez4b2goAfIXn+YyMDEEQtFotx3FNTU3OcavV+txzz+3atUuhGEMFYd8KALyjre2bxTd9+eWXI+MajUYQhPLy8uLiYufIm2+++Z3vfCchIWFM82NtBQDeER096cUXX3S+jo2NJYR0dnbu2LFj69atarW6oKBg48aNhBBRFH/3u98dOXJkrPOjrQDAO/R67eLFi0ePGAyG6urqzMzM5OTkurq6xMREQsjRo0djY2Mfeuihsc6PtgIAX+F53mg0lpSUnD17Nj4+vqamhhCye/fu9PT0ccyGtgIAH0pJSWltbR098t57741vKuyyAwAb0FYAwAa0FQCwAW0FAGxAWwEAG9BWAMAGtBUAsAFtBQBsQFsBABtwLTsA3J1IVZTTuvsEp6Z4igQASK5XtNywmNx84JrNTPEUCQAAgrYCAFagrQCADWgrAGAD2goA2IC2AgA2oK0AgA1oKwBgA9oKANiAtgIANqCtAIANaCsAYAPaCgDYgLYCADagrQCADWgrAGAD2goA2IC2AgA2oK0AwIfq6+sTEhL0en1CQsLBgwcJIXa7vbCwcOLEifPnz+/q6vJ8KrQVAPiKw+HIzc3dvn17T0/P+vXr8/PzCSHbtm0zmUznzp2bN2/eunXrPJ8NT5EAgLvrFxU3RJW7DziUtz/zxm6319bWpqWlCYKgUqn0ej0h5P3339+1a1dYWNjatWtPnz7teQa0FQDcnc1hHxQtbj5gEW2U3vrMG57nMzIyBEHQarUcxzU1NRFCzp0794c//CE1NfXBBx/ctWuX5xlwJggA3nH96vV3brp06dLIuEajEQShvLy8uLiYEGIymSilp06dWrx48fPPP+/5/GgrAPAOs2D+7KaBgQFCSGdnZ2lpKSFErVYXFBS0tbURQiZOnPjiiy9GR0cXFRWdPHnS8/lxJggA3hH3YNw777wzesRgMFRXV2dmZiYnJ9fV1SUmJhJCnnrqqd27d5eUlLzzzjuPP/645/NjbQUAvsLzvNFoLCkpiYqK2rNnT1VVFSFk06ZNR48enTRp0pEjR3bu3On5bFhbAYAPpaSktLa2jh6ZPHnyoUOHxjEV1lYAwAa0FQCwAW0FAGxAWwEAG9BWAMAGtBUAsAFtBQBsQFsBABvG3FYvvfTSkSNHLBZ338YGAPC6MbfVb37zmyeffDIqKio7O/u3v/3tmTNnfBELAOAWY26rvr6+xsbG1atXDw0NlZaWTp069ZFHHnnxxRfr6+t9kQ8AwGnMbRUSEpKcnPzqq68eOnSor6/v448/njRpUmVl5Q9/+MO7HiuKYnx8/OiR3t7erKysyMjI7Ozs3t7esYYBgHvHeHbZr1y58uGHH77wwgtJSUkZGRknTpx47LHHVq5c6f6oysrKefPmdXR0jB6sqKiIi4vr7u6eMmXKli1bxhEGAO4RY74Hw/Tp09vb2+Vy+dy5c9PT019//fUFCxbodLq7Hjhr1qypU6dmZWWNHjQajfv371epVEVFRTk5OZs2bRprHgC4R4y5rdrb2wkhCxcuzMjIWLhw4dy5c3me9+TA1NTU2we7urri4uIIIc4V1u0fcIjiWBMCwLiJ9jv/xMlpWBgX5eZAFXed3PYUCe8ac1t1d3c3NTU1NTXt2bOnrKyM5/knnnhiwYIFCxcuTE9PH+tslFKO45wvxDsVk0wmH+ucADBucsWdf+IGhsnlAZubA3sE8baHSHjZmNtq8uTJy5YtW7ZsGSFkYGDg2LFjb7zxRnl5OSHk9ide3JXBYLhw4cLDDz/c1dUVExNzh0/4tqwBgBnjuXeoxWL529/+duzYsYaGhubmZovFMnv27MWLF49pkoaGhkWLFmVlZdXU1GzcuLGmpiYnJ2ccYQDgHjHmtnryySdPnDgxPDwcFRX1gx/84O23305PT4+Ojh7rPKmpqZTStWvX5ubmxsbGzp07t7a2dqyTAMC9Y8xtNTQ0tGbNmqeeeuqxxx6Ty8e8qTRytuh8odfrDxw4MNZJAOAeNOa2OnHihC9yAAC4h3swAAAb0FYAwAa0FQCwAW0FAGxAWwEAG9BWAMAGtBUAsAFtBQBsQFsBgA/V19cnJCTo9fqEhISDBw86B+fPn8/dtGLFCg+nGs+3mgEAPOFwOHJzc/fu3Zuamrpv3778/Pyuri5KaXt7+8WLF5138VQoPG0hrK0AwFfsdnttbW1aWprZbFapVHq9nhBy5coVq9Wak5MTHR2dl5dnMpk8nA1tBQC+wvN8RkaG2WzWarU5OTlVVVWEkMuXLyclJe3cufPcuXM6na64uNjD2XAmCADeIYrijRs3nK91Ot3IKZ5GoxEEobKysri4uKWlZc6cOUePHnW+tXnz5hkzZng4P9ZWAOAdJ7/4etpNf/nLXwghnZ2dpaWlhBC1Wl1QUNDW1kYI+fzzz5ubm52H8DyvUqk8nB9tBQDeMfuxOTduWrhwISHEYDBUV1c3NjZSSuvq6hITEwkhZrN56dKlbW1tVqt1w4YNS5Ys8XB+nAkCwN3JBiYoet09JUE+0E/orasfnueNRmNJScnZs2fj4+NramoIIQsWLFi3bl1WVlZ/f39GRsb27ds9zIC2AoC7c1Di4tld/yQ67jyekpLS2to6eoTjuMLCwsLCwrFmwJkgALABbQUAbEBbAQAb0FYAwAa0FQCwAW0FAGxAWwEAG9BWAMAGtBUAsAFtBQBsQFsBABvQVgDABrQVALABbQUAbEBbAQAb0FYAwAa0FQCwAW0FAGxAWwEAG3BfdgC4O5UjXC+jbj6gkWk54u4xE98e2goA7s4yoOq7EubmA0KPilLfthXOBAGADWgrAGAD2goA2IC2AgA2oK0AgA1oKwBgA9oKANiAtgIANqCtAIANaCsA8KH6+vqEhAS9Xp+QkHDw4MGR8ZMnT6rV6jFNhbYCAF9xOBy5ubnbt2/v6elZv359fn6+c7y/v3/58uWDg4Njmg1tBQC+Yrfba2tr09LSzGazSqXS6/WEEErp8uXLy8rKxjobvtUMAL7C83xGRoYgCFqtluO4pqYmQkhFRcXUqVOXLVs21tnQVgDgHR3tpxYtWuR8/atf/SoxMdH5WqPRCIJQWVlZXFy8ZcuW+vr6Q4cOjWN+tBUAeEdMTOxrr73mfP3AAw8QQjo7O3fs2LF161a1Wl1QULBx48YjR440NjbyPO/8GMdxx48fX7BggSfzo60AwDs04dqRtZWTwWCorq7OzMxMTk6uq6tLTEwsLy8vLy93vstxHKXu7vB3C7QVAPgKz/NGo7GkpOTs2bPx8fE1NTXfZja0FQD4UEpKSmtr6x3fGtPCiuAKBgBgBdoKANiAM0EA8IBpArke6e4D/RcJ9e3qB2srAGAD2goA2IC2AgA2oK0AgA1oKwBgA9oKANiAtgIANqCtAIANaCsAYAPaCgDYgLYCADagrQCADWgrAGAD2goA2IC2AgA2oK0AgA1oKwBgg5/aqre3NysrKzIyMjs7u7e3d/Rb8+fP525asWKFf/IAAHP81FYVFRVxcXHd3d1TpkzZsmXLyDiltL29/eLFiwMDAwMDA9u2bfNPHgBgjp/aymg0FhUVqVSqoqKiffv2jYxfuXLFarXm5ORER0fn5eWZTCb/5AEA5viprbq6uuLi4gghzhXWyPjly5eTkpJ27tx57tw5nU5XXFx8y4EOUfRPQgAghIj2O//EaZTcfVp3f0WEEc7H2fz0zBtKKcdxzhfiqAKaM2fO0aNHna83b948Y8aMWw6UyeX+SQgAhBC54s4/cXSIWG64O9BqImRsDzMdMz+trQwGw4ULFwghXV1dMTExI+Off/55c3Oz8zXP8yqVyj95AIA5fmqrrKysmpoaSmlNTU1OTg4hpKGhgRBiNpuXLl3a1tZmtVo3bNiwZMkS/+QBAOb4qa3Wrl371VdfxcbGnjp16pVXXiGEpKamEkIWLFiwbt26rKysmJiY3t7eiooK/+QBAOb4ad9Kr9cfOHBg9AillBDCcVxhYWFhYaF/YgAAu3AtOwD4UH19fUJCgl6vT0hIOHjwoJvBu0JbAYCvOByO3Nzc7du39/T0rF+/Pj8/39WgJ9BWAOArdru9trY2LS3NbDarVCq9Xu9q0BN+2rcCgHsQz/MZGRmCIGi1Wo7jmpqaXA16Am0FAN5x9Wr39u3bna+XLl16//33O19rNBpBECorK4uLi1taWtwMuoczQQDwDotl+PRNQ0NDhJDOzs7S0lJCiFqtLigoaGtrczXoCaytAMA7YmMfGFlbORkMhurq6szMzOTk5Lq6usTERFeDnkBbAYCv8DxvNBpLSkrOnj0bHx9fU1PjatATaCsA8KGUlJTW1lZPBu8K+1YAwAa0FQCwAW0FAGxAWwEAG9BWAMAGtBUAsAFXMADA3clkROG2LeTyYHmKBAAwLdxhi7Db3H1AtHI+fowEzgQBgA1oKwBgA9oKANiAtgIANqCtAIANaCsAYAPaCgDYgLYCADagrQCADWgrAGAD2goA2IC2AgA2oK0AgA1oKwBgA9oKANiAtgIANqCtAIANaCsAYAPaCgDYgLYCADbgKRIAcHdqlRilc/eQCJ2acj5+6A3WVgBwd3KLyPVb3P1ltt3xkTf19fUJCQl6vT4hIeHgwYPOwf3798+cOVOv1ycnJ58+fdrDDGgrAPAVh8ORm5u7ffv2np6e9evX5+fnE0LOnz+fl5dXVVXV3d2dnZ3tHPQE2goAfMVut9fW1qalpZnNZpVKpdfrCSFnzpz5yU9+8r3vfS80NPSZZ57p6OjwcDbsWwGAr/A8n5GRIQiCVqvlOK6pqYkQsmjRokWLFhFCRFFcu3btj3/8Yw9nw9oKALzDarN23WS1WkfGNRqNIAjl5eXFxcUjg4cPH05KStLpdJWVlR7Oj7YCAO841fb1d25qaWkhhHR2dpaWlhJC1Gp1QUFBW1sbIYRS+otf/GL9+vV79uzZvHmzQuHpGR7OBAHAOxJnPfbXz5pHjxgMhurq6szMzOTk5Lq6usTEREJIc3Oz0Wj861//qlAoBEEghGg0Gk/mR1sBgK/wPG80GktKSs6ePRsfH19TU0MIaWho6OjoiIiIGPkYpe6u5BqBM0EA8KGUlJTW1taenp7m5ub4+HhCyMsvv0z/Jw+nQlsBABvQVgDABrQVALABbQUAbEBbAQAb0FYAwAa0FQCwAW0FAGxAWwEAG9BWAMAGtBUAsAHfagaAu5MrOIXK3VMi5EqO+PgpEmgrALg7td0iswy5+UC4zcJ5/P3k8cGZIACwAW0FAGxAWwEAG9BWAMAGtBUAsAFtBQBsQFsBABvQVgDABrQVALABbQUAbEBbAQAb0FYAwAa0FQCwAW0FAGxAWwEAG9BWAMAGtBUAsMF/bdXb25uVlRUZGZmdnd3b23vXcQCA0fzXVhUVFXFxcd3d3VOmTNmyZctdxwEARvPffdmNRuP+/ftVKlVRUVFOTs6mTZvcj/8Tpebhy1+f2eG3nD5FiUiI44suo9RBfMUqDnX0faEbvCB1EJ/oHb7aZ43uGuySOoivXB2+qraqpU7hEkd9fOP3ERqN5tq1a6GhoUNDQ5MmTTKZTO7HnQrzVn+4/18/2xzn46dq+B5HHOF8uNQp/sX5/7+3/r1SKup4tYzz7ZpddDjkMgm2XB3EIZPbQxQhYz3QLto5wsnlcl+k8q6frfnZf0zpgtAAAAUpSURBVDz/H7ePW7otoll0c+CN3hunLp9anLXYZ9H82FZqtfrGjRshISGDg4MTJ040m83uxwEARvPff6AMBsOFCxcIIV1dXTExMXcdBwAYzX9tlZWVVVNTQymtqanJyckhhDQ0NNxxHADgdv47E+zr68vNzf3yyy/nzp1bW1ur0+k4jqOU3j7unzwAwBgK49LT05OZmRkREZGVldXT0zP6rT/96U8zZszQ6XQLFy7s6OiQKuHt3GR2+vrrr8PCwvwfzA1XmW02289+9rMJEybMmzfv4sWLEia8havADQ0Ns2fP1mg0s2fPbmxslDAhu3At+zi5ukzs/PnzeXl5VVVV3d3d2dnZ+fn5Eoa8hftL2/r7+5cvXz44OChJNldcZd62bZvJZDp37ty8efPWrVsnYcJbuAqcl5f38ssv9/T0rFmzJi8vT8KEDJO6Llk1bdq0trY2SmlbW9u0adNGxo8dO/bcc885X1+9ejUqKkqafHfiKjOl1OFwLFmy5I9//GOg/ZFwlTkxMfGLL76glJpMps8++0yyfLdxFTghIaGqqqqnp2fnzp3Tp0+XLiDDAuuPJkPUavXg4CCldHBwMDw8/PYP2O32FStWFBYW+j2aS24yb9q0aeXKlZTSQGsrV5kjIyNXr14dERHx2GOPffXVV9IFvJWrwC0tLSNLhJaWFukCMgxnguNEKXVeqkopFcVbr5o7fPhwUlKSTqerrKyUIt2ducp87Nix+vr6W79FEBhcZTaZTJTSU6dOLV68+Pnnn5cu4K1cBV69evWqVasuXbpUWlpaVlYmXUCWSdaTjHvooYdOnz5NKT19+vTDDz88Mu5wOMrKygJtf93JVeaXX375lj8Vx48fly7m/+Aqc3R09KVLlyil3d3darVasny3cRVYrVZ3d3dTSq9fv67RaCTLxzKsrcbJ1eVjzc3NRqPxo48+MhgMgiAIgiBx0FFcZS4vLx/5A0EIoZQuWLBA2qgjXGV+6qmndu/ebbFY3nnnnccff1zilKO4Cjxr1qzq6mpBEN59993Zs2dLnJJR0pQk+3p7ezMyMmJiYrKysvr6+ujNn/Py8vKA/TfsKvNoARWYus7c3d395JNP6nS65OTkb775RuqY/+IqcFtb27x58zQazbx585zb8DBW/rs6FADg28CZIACwAW0FAGxAWwEAG9BWAMAGtBUAsAFtBQBsQFuBb8XHxwfB3fQhEKCtAIANuDoUfMtsNlNKNRqN1EGAeWgrAGADzgTBC6ZNm8Zx3OXLl5cuXRoREZGQkLBy5UrnN7qxbwXegrYCL7Db7YSQp59+WqvVlpWVxcbGvvXWW/Pnz7fZbFJHg+DhvyfLQxBz3nbumWeeeeGFFwghq1atKigo2LVr19tvvy11NAgeWFuBF4iiqFQqV6xY4fxbjuNWrlxJbt7aCcAr0FbgBaIoRkdH8zw/MvLAAw8QQs6cOSNdKAg2aCvwAlEUu7u7rVbryMj58+cJIVOmTJEuFAQbtBV4gd1ut9lsI7tUlNI33niDEJKSkiJpLggquN4KvECn0zkcDkLIsmXLpk+ffvTo0U8++WTmzJmtra2zZs1yPlBD6ozAPKytwAtEUYyJiTl27NiZM2fKy8s7OzuLi4ubm5tH72QBfEtYW4EXhIaGxsXFtbe3Sx0EghnWVuAFtz//FcDr0FbgBc5r2QF8CteygxesWbNmwoQJUqeAIId9KwBgA84EAYANaCsAYAPaCgDYgLYCADagrQCADWgrAGAD2goA2PD/AVkEr208DVhNAAAAAElFTkSuQmCC"
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GKS: invalid bitmap size\n"
     ]
    }
   ],
   "source": [
    "# Set up the problem and initial guess, solve by VFI\n",
    "sp = SearchProblem(;w_grid_size=100, π_grid_size=100)\n",
    "v_init = fill(sp.c / (1 - sp.β), sp.n_w, sp.n_π)\n",
    "f(x) = T(sp, x)\n",
    "v = compute_fixed_point(f, v_init)\n",
    "policy = get_greedy(sp, v)\n",
    "\n",
    "# Make functions for the linear interpolants of these\n",
    "vf = extrapolate(interpolate((sp.w_grid, sp.π_grid), v, Gridded(Linear())),\n",
    "                Flat())\n",
    "pf = extrapolate(interpolate((sp.w_grid, sp.π_grid), policy,\n",
    "                Gridded(Linear())), Flat())\n",
    "\n",
    "function plot_value_function(;w_plot_grid_size = 100,\n",
    "                            π_plot_grid_size = 100)\n",
    "  π_plot_grid = range(0.001,  0.99, length =  π_plot_grid_size)\n",
    "  w_plot_grid = range(0,  sp.w_max, length = w_plot_grid_size)\n",
    "  Z = [vf(w_plot_grid[j], π_plot_grid[i])\n",
    "          for j in 1:w_plot_grid_size, i in 1:π_plot_grid_size]\n",
    "  p = contour(π_plot_grid, w_plot_grid, Z, levels=15, alpha=0.6,\n",
    "              fill=true, size=(400, 400), c=:lightrainbow)\n",
    "  plot!(xlabel=\"pi\", ylabel=\"w\", xguidefont=font(12))\n",
    "  return p\n",
    "end\n",
    "\n",
    "plot_value_function()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='odu-pol-vfi'></a>\n",
    "The optimal policy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAGQCAIAAADZeHv+AAAABmJLR0QA/wD/AP+gvaeTAAAfi0lEQVR4nO3deXhU5b3A8d97JitJyIZiEpMoyL5EQayyFbyKiA6RXmpdEKu2FWkstsgiCvRRXIJeL1httdDgUi3WXhC81KUiwUpuBaTKIgkgsoWwB7Jv57z3j4EYgSwwIWde8v08/jE5c87kdx6TL2fOnJkorbUAQMCz3B4AAJqEWgEwA7UCYAZqBcAM1AqAGagVADNQKwBmcKdWS5Ys6dmzZ0xMzODBg7ds2VK7vLCw0Ov1xsXFjRw5srCw0JXZAAQmF2q1a9euMWPGzJs3r6CgYOTIkffcc0/tXZmZmampqQUFBSkpKbNnz2752QAELNXy17JnZ2e/+eab8+bNE5GDBw9269bt0KFDvru6dOmyZMmSrl275ubmpqen5+XltfBsAAKWC7WqZdt2RkaGZVkvvfSSb0lkZOTBgwfDw8PLy8vbt29fVFRUd/1J4x/+vxWf1X7ZsWPHyzpe1qITNwctWkSUKLcH8ZejbUt53J7Cf9rR2lLGn8A1/efqzgfu7tC1Q8PrBLXMKKf6+OOPJ0+ePGzYsFmzZtUu1ForpXw3bNs+aZO/Zi3oVXW0zYkvy/PWbGihYQGcQ5u1rqwon/VKIyd/XKiV1nratGmrVq1auHBh586d696VmJi4e/fuTp065efnJyUlnbylkruV58KWmxRAS3hDHFGNHxW6cACck5OzePHipUuXJiYmlpSUlJSUiEh2draIeL3erKwsrXVWVlZ6enrLzwYgYLlwbJWdnZ2XlxcbG1u7RGs9dOhQrfWMGTPuvPPO5OTkPn36vPHGGy0/G4CA5eZZ9jOVGh7/XFUxzwSB88wb2rnoFxNnvZzZ8GrGvxQCoJWgVgDMQK0AmIFaATADtQJgBmoFwAzUCoAZqBUAM1ArAGagVgDMQK0AmIFaATADtQJgBmoFwAzUCoAZqBUAM1ArAGagVgDMQK0AmIFaATADtQJgBmoFmKFapMrtGdxFrQAz3OvUDHOq3Z7CTdQKgBmoFQAzuPCX5YGAtUP0n7WzXeu9ooNF2ou6Vlk/VlZwnXU+1c574mzROlRUJ6VuEtVfWU28d6V2/lecrVrbIqmiblLWDcry3X2XU7Nb9D+s4Ne1/YHWJaI7ibpOWV5lKZEhJ54D+m5kf2+i1oJaAcct186T2nZELhU1UFnlIv/Wzh+1XSB6ovL41vmDtt/WTqRIL6W0yGrtrBKZIDJKWY3e+4K2F2mnnajeyqoWvVHrTG1/LnqG8tT27Cltf6adbqKUqM2i12t7o+hpyjNBeV7XdqHIhBOTtELUCjjuNe04IiOUNUl5lIiIHFTWHU7NCu34avW16Le1c5moZ62gWBER2an0eKdmvrZvUta2Bu/9WutF2klXVoby+I6LSpU8o+1s7Vwp6uYTx1//1s4frKDLRIlIvuhJjv2RdoaLNUpZi7RTKHqUar1nb1rvngMnuV1Zk5TnbmWpE0suEBUvquTElwu1IyK/VJ7YE0tSRf1EeVJEHRTd8L3viBMiUpsqEYkQeVh5ROQf2qmdYazy+FIlIkmifqksEfkf+W6F1oxjK+C4G08ctpSK7BC9Reu1oveJrl1hh9bBImlK1d1qrLLGKqsp99oi9zo1J31TJfJtnW9x9fc3/4GyPNrepbUoAbUCjqsQma/tz7TeJ9oSuURUF6XaihSdWGGf6Dip95lYw/ceFG2L7K4TplrldW7Hfz9LQSLxovafbqtWiFoBx83R9gfauVqp8eLpp6xwERG5S+uiE7GIF1UoWsv3iuKI+JY0fG97UeUif7Ma+Y07LDqxzgM4IkdEX8iRlYhw3gqotVI7bUWeUkGDT6RKi5TWOa5JUVIhsvn7RzqvaPs/nOoNohu7Vx0WfdJRUoHo57W9rM55q8/191ZYo3WNSAqxEhFqBdSKEFVZ53lfjUiWto+IiIgtIiI/Eo+IzHXsoyfW2SN6mXbaiHQT1fC9PxZLizzh2MdO3FsmkqntpdoJrzPDa9r+5kTR8kX/TtsiMkq+u2qhNb9VkGeCwHHDlfqz1j93avorVSXyldaWSIqoXaKf0vZdyvqBUjeK9b527nKqeypliVqjnSqRR5UnRKThey9XaqRYS7Vzh1PdXalwUV9qp0jkemUNqXOyK1nUOKemuyilZLPWlSLXKesqpUTEF7WZuiZJVEarvOqKWgHH/VR52oh6XzsfaOdiUT9U6m7l+VLrTF3zL+0MF+sSJZOVJ03UB+Js1DpIdE+l7hLPFSdeyGv43t8oTy9R74uTp7WIThY1UlnXnbiW3ed5K2i+tj/X+oA+fi17+okrKu5Vnhe0vUbrXNGts1ZKa2NebkgNj3+uqvhCt8cAzgXfO29a51tq3tDORb+YOOvlzIZX47wVADNQKwBmoFYAzMBZdiAg/M4K4t2ADaNWQECIcXuAwMczQQBmoFYAzECtAJiBWgEwA7UCYAZqBcAM1AqAGQy73mq76KONr/WdeFHx52oWAC3KpFodC2rztBVmqTOYuabiwH/pms58UCxgPpNqVaGCLhj0elB4QtM3KS/4x6/XTpytq3sQLMBwJp23Co5ItIIiz2iT8ITrI3/w0iQV/DV/NQQwnEnHVmcn/MJBctWLD30+PsIKtsQJVUGWCrzjLC2j7bJb3J4CCGTnf61EJLz94NAbc5yqYyLi1BQ7OuDe6+7UlM3/99SvKw9Pcqpb42dHAk3QKmolIlZwWyu4rdtTNCRk6P+uXvPgw4fXPutUh7g9DBCAWkutAp8V1Cbm6j/uWv3grQc+S7aCU5zqRKcq8J6yfkeJXKtUe16+QEuhVoFEeaKverG65NsDpbvzy3ZJ2V5L17g9U710dcnr+cvu1fZ/mvViDYxFrQKMsoKjOgZHdQxvfFX3hXW+/89rf/NRyfbnneoot4fBeY9a4ewFR14aO2TRkfWzRu9YeLknZEhN+ZVKhQXwc0MlcmaXwCCQUCv4SUX1nt6my/hv9q/cmv/3uUe+cBzb7ZHqpbWd4gkZW1M+6Pt/cxRGoFZoBp7Q+MiUH0nKjwL6ZVcR0XZxwfL/zntxbunO+52qYRKAl96hXvwDg9ZEedokDosZujS0/4IXw5N+ZoV8xZsczMGxFVqj0Lg+odd9VLTnvWkbn7nQrmgTwMdYlkg3u7K3trsrK87tYdxFrdBaKSsiOb1N0ojyw2vLAvkIy67ac3TD8oM5JUe/HqDkMae61T4holZo1ZQVHHbBNW5P0Yjwi4ZI1wcj7YovVt39UFHuTKe6dX5qG7UCzKA8YTGD3srPfWHstgU/d6pjAvl4UCRcVJySGJF4Uc11MEitAGMo5Ynq9uuQ9kNf2/Ky0oF7pYiIVlXHaqoOVVYe7Sv2406NpzkelFoBhgmNuzz06pfdnqJJtFO1MeeeGUc3/Nap8f/DRVrtCTsA55yyQmKuydrU7poMK6TI70fj2ArAOaQ8oTE/+MOhTc/85NuF5XZ5fatNj4pu9KGoFYBzTFlRPadF9ZxW3/3H8n6vYhq/mIxnggBcFhzV8eKUDo2uRq0AuEx5wkNDwxpdjVoBMIM7tbJtu2vXrictHDBggDph3LhxrgwGIGC5cJZ97ty5b731Vl5eXt2FWuvc3Nw9e/ZER0eLSFAQp/8BfI8Lx1a9e/eePn36SQv3799fVVWVnp6ekJAwZsyYoqLTXJyhdUC/1QDAWWvK381zoVZDhw69+eabT1q4b9++fv36zZ8/f+fOndHR0RMmTGj5wQAEskB5wnX55Zd/8sknvtvPPPNMjx49Tl1HBfCHEAHwh9WE9z4HymuC69aty8nJ8d0OCQkJDQ11dx4Agcb9WmVnZ4tIaWnpqFGjNm/eXFVV9cQTT9xyyy1uzwUgsLhfq6FDh4rIwIEDZ86c6fV6k5KSCgsLMzMz3Z4LQGBx7bxV7Qt8vhtKqfHjx48fP96teQAEOPePrQCgKagVADNQKwBmoFYAzECtAJiBWgEwA7UCYAZqBcAM1AqAGagVADNQKwBmoFYAzECtAJiBWgEwA7UCYAZqBcAM1AqAGagVADNQKwBmoFYAzECtAJiBWgEwA7UCYAZqBcAM1AqAGagVADNQKwBmoFYAzECtAJiBWgEwA7UCYAZqBcAM1AqAGagVADNQKwBmoFYAzECtAJiBWgEwA7UCYAZqBcAM1AqAGagVADNQKwBmoFYAzECtAJiBWgEwA7UCYAZqBcAM1AqAGagVADNQKwBmoFYAzECtAJiBWgEwA7UCYAZqBcAM1AqAGagVADNQKwBm8LdW27Ztu+OOOy6++OLw8PDc3Nzf/va3q1evbpbJAKAuv2q1devWfv36vffee3379q2oqBCRtWvXDhkyZO3atc00HgAc51etZsyYERUV9fXXXy9ZssS35J133unevfuTTz7ZHLMBwHf8qtXy5csnT56cnJxcuyQ8PPzOO+/817/+5fdgAPA9ftWquro6ISHhpIVJSUmlpaX+PCwAnMqvWl1zzTV///vfT1r44YcfpqWl+fOwAHCqIH82njNnztVXXx0WFnbbbbeJyMaNG996661XX311xYoVzTQeABzn17FV586dN2zY4PF4HnjggTZt2jz44INr1qzJyckZPHhwc80HAD5+HVuJSFJS0gsvvNAsowBAA/w6tlKnExwc3L59+27dut17772ff/55cw0KoJXzq1Y5OTnJycl33XXXsmXL1q1b9/777999992XXnrp3/72tzlz5sTExAwePPijjz5qrlkBtGZ+PRN8/vnnR4wY8fLLL9cuGT58+Lhx47KyshYsWHDDDTfExcU9/vjjw4YN83tOuCn/42HVJd9ecsvWZnk07VSJ1soT2iyPhtbD36tDTy3R9ddfv2zZstrb69ev9+db4Pyz95Obd77X0+0pYB6/anXxxRef+pbANWvWtGvXznd748aNF154oT/fAoHggqteTBiyyO0p0Nr59UzwkUceGTt2bHh4+JgxYxISEvbt2/fmm28+++yzWVlZpaWlCxYsmD59+rRp05prVrglpG1nt0cA/KvV7bffHhER8fjjj8+cOVNrLSI9evT4y1/+cuutt27duvW111578sknx40bd9ptbdvu0aNHbm5u7ZLCwsKxY8euWrVq4MCBr732WmxsrD+z4bRqz0AdzXupaNufoi65PbbHpNK9H5TseLvq2NfasYPbXhaV+uPI5FGirFO3ql3S6CZlez8s3rGw6ugm5QkLiekRmfrjNhddKyI73u3kW8F3o7nOhaE18PfzrUaOHLl27dry8vJvv/22rKxs48aNt956q4h06tRpzZo148ePt6zTfIu5c+f2798/Ly+v7sLMzMzU1NSCgoKUlJTZs2f7ORgaULzzr0c3zxFRVkjMkfVPHFz9YFXR1tD4K8Pi+1QXbTu0burBtQ+JdurbvNFNCjdmHlidUVm4PjTu8uC2Xcr3rzzwr/uLt/9ZROJ6z/CExvtuxPWe0TL7i/ODv1eH+oSGhl5yySW+22VlZd98802vXr0aWL93794dO3b0er11Fy5evHjJkiWhoaEZGRnp6elPP/30yZvpZhkWUrhxdvsBr4VfcE3FodX7Ns2OuvSOuF6PKStYRJyakkNfTCnNfz/sggFRl/zk1G0rDn1etP31BjapPPLlsW3zQ6K7te+f5QltJyLVxd8UrBxd+PV/Rab+uG2Hu4q3v2FXHm7b4a4W3msEMt2EX29/a1VRUbFz507f00Cf5cuXP/LII0VFRQ1sNXTo0FMX5ufnp6amiojvCOvUFZz6/7XHGYnt/pvwC/qLSNG2BcoTGtfrUV93RMQKioy/YlZZwUele5aetlaNbnJs23wRies1zZcqEQmO6ti208/K962oqdgXHJHaEnsI0zh247/dftXq/fffHz16dFlZWd2FwcHBv/jFL87i0bTWSinfDdu2T13htE8qcRbCLrjGd6O6eJs4NXs/ufmUVVRV0ZbTbtvoJtXFW5UVHBZ/Vd37Yrr8MqbLL/2fHOcrj8fT6Dp+/f4/8sgjI0aMKCsre/jhh++//37HcTZs2JCcnPzAAw+cxaMlJibu3r1bRPLz85OSkvwZDA0LanP8AxRrKvZpbVeXfHvSfyJa15z+Q8oa3aSmLN8TdmHdM+5As/Dr2CovL2/KlCnh4eG33Xbbfffdp5Tq2bPnpEmTfv3rX5/RG26ys7OHDBni9XqzsrKeeuqprKys9PR0fwZDw5R1/P97UHiiU1OaPHxV07dtdBNP6AV25WERLaK+W6odLVqJomI4a3796LRr1853gqlXr15btmzZu3eviHTo0OFM38zsO401Y8aM9evXJycnb9q06bHHHvNnMDRRcFRHu+JgTVl+3YU1ZXsOfzm9ZOc7Z7dJcFRHbZdXHvmq7gqFm57duaRrxWH+vAjOnl+1GjBgwJw5c957772QkJArrrhizpw5JSUlb7/9dmJiYlM2rz0377sRExOzbNmyPXv2LF26NDo62p/B0ERtO/5URB9c+xu7qtC3xKkpPbRuavGOhSqozdlt4nux78j6x+3KI74Vqkt2FO/8qxUUERr73YfKarvy3OwTzlt+PRPMzMwcMmTI3LlzvV7v73//+xtvvPHZZ5+NjY1duHBhc82Hcyqs3Q+iLr29+Nu/5P/jP0JjL1dB4RWHPneqjkUkj4xIvPHsNglv/8PIlP8s2fU/+R8PC4vvI2KVH/xM25Xt+j7neyezCooQkYNrHgyKSInrxUE0msqvWv3ud797+eWXe/ToISJpaWn5+fn5+fnt27cPDg5upvFwzsWnPR4Wd2Xxrr9VHl0vIsGRHaIuvT3i4pENnGBqdJN2fZ4Oa3dVya5FFYfXKSsoNK5vTJcHwtpd7bs3ttuEw+sfLz/wTys4hlqh6VTdS6XOVFhYWGVlZURExLXXXjt8+PDhw4d36NChGYc7SVTCgLi+L1rBUefuW+C0mvcTY4CTlO//dPZDl429tZFPSPfrvNXRo0dXrlw5ZcqU8vLySZMmdezYsUuXLg899NAHH3zgz8Mi0Ginyu0RAP+eCYaFhQ0ePNj3NyOqq6uXL1/+1FNPzZ07d+7cuf4csiGAaLvi8JqasnzlCXd7FLR2/r7zZv/+/Z999tnKlSs//fTT9evXK6X69u07ZMiQ5pgN7juwOqOs4GMRCW/P3zGCy/yqVbdu3XJzcz0eT58+fYYNG/bkk08OHDiQiw/OJ0FtLg6NuyIsvm9057N5fwLQjPyqle/TqQYNGjRixIhBgwb16dMnJCSkmQZDQIjr9ajbIwDH+XWWvaCg4J133klLS1u4cOGAAQOio6OHDh06ffp0/s4NgGbn17HVRRddNHr06NGjR4tIcXHxihUrnnvuuVmzZkmd69QBoFn4e5a9srJy9erVK1asyM7OzsnJqaysTEtLGz58eLMMBwC1/KrVddddt2rVqoqKivj4+Ouvv/6VV14ZNmxYQkJCcw0HALX8qlV5efm0adNuuOGGvn37NuXDtADgrPlVq1WrzuBzkQDAH3w0GgAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYwYVaFRYWer3euLi4kSNHFhYW1r1rwIAB6oRx48a1/GwAApYLtcrMzExNTS0oKEhJSZk9e3btcq11bm7unj17iouLi4uL58yZ0/KzAQhYLtRq8eLFGRkZoaGhGRkZixYtql2+f//+qqqq9PT0hISEMWPGFBUVtfxsAAKWC7XKz89PTU0VEd8RVu3yffv29evXb/78+Tt37oyOjp4wYcJJGzqO06KDAmgptm03uk5QC8xxEq21Usp3o+6Il19++SeffOK7/cwzz/To0eOkDS3FawLA+cnyNP7b7cLvf2Ji4u7du0UkPz8/KSmpdvm6detycnJ8t0NCQkJDQ0/eUrXUiABalmrCr7cLtfJ6vVlZWVrrrKys9PR0EcnOzhaR0tLSUaNGbd68uaqq6oknnrjllltafjYAAcuFWs2YMWP9+vXJycmbNm167LHHRGTo0KEiMnDgwJkzZ3q93qSkpMLCwszMzJafDUDAUlprt2doqqiEAXF9X7SCo9weBEBzKt//6eyHLht76+CGV+O8NQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZnCnVoWFhV6vNy4ubuTIkYWFhY0uBwB3apWZmZmamlpQUJCSkjJ79uxGlwOAO7VavHhxRkZGaGhoRkbGokWLGl0OAEGufNf8/PzU1FQR8R1JNbrcJyjIKtrykmUF+75UlmVZ5p13s2tqDJ38JNXVVcHBIW5P4S+tHdu2g4KC3R7EX0b/XNWU7y8pHtPoau7USmutlPLdsG270eU+r/9+8rvvvitS5fsyJSUlJSWlpUYGcA7dfPM1ja7jTq0SExN3797dqVOn/Pz8pKSkRpf7eL1er9fbspMCCBTuHDd6vd6srCytdVZWVnp6uohkZ2efdjkA+Citdct/16NHj955551fffVVnz593njjjejoaKWU1vrU5S0/G4AApXFuHDly5Oabb46NjfV6vUeOHKl717vvvtujR4/o6OhBgwbl5eW5NWETNbAjPhs2bGjTpk3LD3YW6tuX6urqBx54oF27dv3799+zZ4+LEzZFfXuRnZ2dlpYWGRmZlpa2cuVKFyc8R4x8BcEI9V07tmvXrjFjxsybN6+goGDkyJH33HOPi0M2RcMXwR07duynP/1pWVmZK7Odqfr2Zc6cOUVFRTt37uzfv//MmTNdnLAp6tuLMWPGPProo0eOHJk2bdqYMY2/xGYet3N53urcufPmzZu11ps3b+7cuXPt8hUrVvzsZz/z3T5w4EB8fLw78zVZfTuitXYc55ZbbnnnnXdM+UGqb1+uuOKKL7/8UmtdVFS0du1a1+Zrmvr2onv37vPmzTty5Mj8+fO7devm3oDnihk/ZCaKiIgoKyvTWpeVlUVFRZ26Qk1Nzbhx48aPH9/io52ZBnbk6aefnjhxotbalFrVty9xcXFTpkyJjY3t27fv+vXr3RuwSerbizVr1tQehaxZs8a9Ac8VngmeK7rBa8c+/vjjfv36RUdHz507143pzkB9O7JixYoPPvjg6aefdm+0M1bfvhQVFWmtN23aNHz48J///OfuDdgk9e3FlClTJk+evHfv3kmTJk2dOtW9Ac8Z1zp5vrvsssu2bNmitd6yZUunTp1qlzuOM3XqVCPOr/vUtyOPPvroST9L//znP90bs0nq25eEhIS9e/dqrQsKCiIiIlybr2nq24uIiIiCggKt9aFDhyIjI12b75zh2Opcqe+aspycnMWLFy9dujQxMbGkpKSkpMTlQRtT347MmjWr9sdIRLTWAwcOdHfURtW3LzfccMOrr75aWVn5xz/+8corr3R5ysbUtxe9e/f+05/+VFJS8vrrr6elpbk85bngTiRbgcLCwhEjRiQlJXm93qNHj+oTv9KzZs0y639BfTtSV+DvhU99+1JQUHDddddFR0cPHjx469atbo/ZiPr2YvPmzf3794+MjOzfv7/vNPx5xp2rQwHgTPFMEIAZqBUAM1ArAGagVgDMQK0AmIFaATADtQJgBmqFQNe1a1ff2+LQylErAGbgWnYEutLSUq11ZGSk24PAZdQKgBl4JoiA0LlzZ6XUvn37Ro0aFRsb271794kTJ/o+oILzVvChVggINTU1InLTTTe1bdt26tSpycnJzz///IABA6qrq90eDYHCnb9+CpzE9xmYd999969+9SsRmTx58n333bdgwYJXXnnF7dEQKDi2QkCwbTs4OHjcuHG+L5VSEydOlBOfMwcItUKAsG07ISEhJCSkdsmll14qItu3b3dvKAQWaoWAYNt2QUFBVVVV7ZJdu3aJSEpKintDIbBQKwSEmpqa6urq2rNUWuvnnntORH74wx+6OhcCCNdbISBER0c7jiMio0eP7tat2yeffPLhhx/27Nnziy++6N27t+/vA7k9I1zGsRUCgm3bSUlJK1as2L59+6xZs3bs2DFhwoScnJy6Z7LQynFshYAQHh6empqam5vr9iAIXBxbISCc+uesgZNQKwQE37XsQAO4lh0BYdq0ae3atXN7CgQ0zlsBMAPPBAGYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCYgVoBMAO1AmAGagXADNQKgBmoFQAzUCsAZqBWAMxArQCY4f8Bn3el+qeSH0gAAAAASUVORK5CYII="
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function plot_policy_function(;w_plot_grid_size = 100,\n",
    "                              π_plot_grid_size = 100)\n",
    "    π_plot_grid = range(0.001,  0.99, length = π_plot_grid_size)\n",
    "    w_plot_grid = range(0,  sp.w_max, length = w_plot_grid_size)\n",
    "    Z = [pf(w_plot_grid[j], π_plot_grid[i])\n",
    "            for j in 1:w_plot_grid_size, i in 1:π_plot_grid_size]\n",
    "    p = contour(π_plot_grid, w_plot_grid, Z, levels=1, alpha=0.6, fill=true,\n",
    "                size=(400, 400), c=:coolwarm)\n",
    "    plot!(xlabel=\"pi\", ylabel=\"wage\", xguidefont=font(12), cbar=false)\n",
    "    annotate!(0.4, 1.0, \"reject\")\n",
    "    annotate!(0.7, 1.8, \"accept\")\n",
    "    return p\n",
    "end\n",
    "\n",
    "plot_policy_function()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code takes several minutes to run\n",
    "\n",
    "The results fit well with our intuition from section [looking forward](#looking-forward)\n",
    "\n",
    "- The black line in the figure above corresponds to the function $ \\bar w(\\pi) $ introduced there  \n",
    "- It is decreasing as expected  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take 2: A More Efficient Method\n",
    "\n",
    "Our implementation of VFI can be optimized to some degree\n",
    "\n",
    "But instead of pursuing that, let’s consider another method to solve for the optimal policy\n",
    "\n",
    "We will use iteration with an operator that has the same contraction rate as the Bellman operator, but\n",
    "\n",
    "- one dimensional rather than two dimensional  \n",
    "- no maximization step  \n",
    "\n",
    "\n",
    "As a consequence, the algorithm is orders of magnitude faster than VFI\n",
    "\n",
    "This section illustrates the point that when it comes to programming, a bit of\n",
    "mathematical analysis goes a long way"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Another Functional Equation\n",
    "\n",
    "To begin, note that when $ w = \\bar w(\\pi) $, the worker is indifferent\n",
    "between accepting and rejecting\n",
    "\n",
    "Hence the two choices on the right-hand side of [(3)](#equation-odu-mvf) have equal value:\n",
    "\n",
    "\n",
    "<a id='equation-odu-mvf2'></a>\n",
    "$$\n",
    "\\frac{\\bar w(\\pi)}{1 - \\beta}\n",
    "= c + \\beta \\int V(w', \\pi') \\, h_{\\pi}(w') \\, dw' \\tag{4}\n",
    "$$\n",
    "\n",
    "Together, [(3)](#equation-odu-mvf) and [(4)](#equation-odu-mvf2) give\n",
    "\n",
    "\n",
    "<a id='equation-odu-mvf3'></a>\n",
    "$$\n",
    "V(w, \\pi) =\n",
    "\\max\n",
    "\\left\\{\n",
    "    \\frac{w}{1 - \\beta} ,\\, \\frac{\\bar w(\\pi)}{1 - \\beta}\n",
    "\\right\\} \\tag{5}\n",
    "$$\n",
    "\n",
    "Combining [(4)](#equation-odu-mvf2) and [(5)](#equation-odu-mvf3), we obtain\n",
    "\n",
    "$$\n",
    "\\frac{\\bar w(\\pi)}{1 - \\beta}\n",
    "= c + \\beta \\int \\max \\left\\{\n",
    "    \\frac{w'}{1 - \\beta} ,\\, \\frac{\\bar w(\\pi')}{1 - \\beta}\n",
    "\\right\\}\n",
    "\\, h_{\\pi}(w') \\, dw'\n",
    "$$\n",
    "\n",
    "Multiplying by $ 1 - \\beta $, substituting in $ \\pi' = q(w', \\pi) $ and using $ \\circ $ for composition of functions yields\n",
    "\n",
    "\n",
    "<a id='equation-odu-mvf4'></a>\n",
    "$$\n",
    "\\bar w(\\pi)\n",
    "= (1 - \\beta) c +\n",
    "\\beta \\int \\max \\left\\{ w', \\bar w \\circ q(w', \\pi) \\right\\} \\, h_{\\pi}(w') \\, dw' \\tag{6}\n",
    "$$\n",
    "\n",
    "Equation [(6)](#equation-odu-mvf4) can be understood as a functional equation, where $ \\bar w $ is the unknown function\n",
    "\n",
    "- Let’s call it the *reservation wage functional equation* (RWFE)  \n",
    "- The solution $ \\bar w $ to the RWFE is the object that we wish to compute  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the RWFE\n",
    "\n",
    "To solve the RWFE, we will first show that its solution is the\n",
    "fixed point of a [contraction mapping](https://en.wikipedia.org/wiki/Contraction_mapping)\n",
    "\n",
    "To this end, let\n",
    "\n",
    "- $ b[0,1] $ be the bounded real-valued functions on $ [0,1] $  \n",
    "- $ \\| \\psi \\| := \\sup_{x \\in [0,1]} | \\psi(x) | $  \n",
    "\n",
    "\n",
    "Consider the operator $ Q $ mapping $ \\psi \\in b[0,1] $ into $ Q\\psi \\in b[0,1] $ via\n",
    "\n",
    "\n",
    "<a id='equation-odu-dq'></a>\n",
    "$$\n",
    "(Q \\psi)(\\pi)\n",
    "= (1 - \\beta) c +\n",
    "\\beta \\int \\max \\left\\{ w', \\psi \\circ q(w', \\pi) \\right\\} \\, h_{\\pi}(w') \\, dw' \\tag{7}\n",
    "$$\n",
    "\n",
    "Comparing [(6)](#equation-odu-mvf4) and [(7)](#equation-odu-dq), we see that the set of fixed points of $ Q $ exactly coincides with the set of solutions to the RWFE\n",
    "\n",
    "- If $ Q \\bar w = \\bar w $ then $ \\bar w $ solves [(6)](#equation-odu-mvf4) and vice versa  \n",
    "\n",
    "\n",
    "Moreover, for any $ \\psi, \\phi \\in b[0,1] $, basic algebra and the\n",
    "triangle inequality for integrals tells us that\n",
    "\n",
    "\n",
    "<a id='equation-odu-nt'></a>\n",
    "$$\n",
    "|(Q \\psi)(\\pi) - (Q \\phi)(\\pi)|\n",
    "\\leq \\beta \\int\n",
    "\\left|\n",
    "\\max \\left\\{w', \\psi \\circ q(w', \\pi) \\right\\} -\n",
    "\\max \\left\\{w', \\phi \\circ q(w', \\pi) \\right\\}\n",
    "\\right|\n",
    "\\, h_{\\pi}(w') \\, dw' \\tag{8}\n",
    "$$\n",
    "\n",
    "Working case by case, it is easy to check that for real numbers $ a, b, c $ we always have\n",
    "\n",
    "\n",
    "<a id='equation-odu-nt2'></a>\n",
    "$$\n",
    "| \\max\\{a, b\\} - \\max\\{a, c\\}| \\leq | b - c| \\tag{9}\n",
    "$$\n",
    "\n",
    "Combining [(8)](#equation-odu-nt) and [(9)](#equation-odu-nt2) yields\n",
    "\n",
    "\n",
    "<a id='equation-odu-nt3'></a>\n",
    "$$\n",
    "|(Q \\psi)(\\pi) - (Q \\phi)(\\pi)|\n",
    "\\leq \\beta \\int\n",
    "\\left| \\psi \\circ q(w', \\pi) -  \\phi \\circ q(w', \\pi) \\right|\n",
    "\\, h_{\\pi}(w') \\, dw'\n",
    "\\leq \\beta \\| \\psi - \\phi \\| \\tag{10}\n",
    "$$\n",
    "\n",
    "Taking the supremum over $ \\pi $ now gives us\n",
    "\n",
    "\n",
    "<a id='equation-odu-rwc'></a>\n",
    "$$\n",
    "\\|Q \\psi - Q \\phi\\|\n",
    "\\leq \\beta \\| \\psi - \\phi \\| \\tag{11}\n",
    "$$\n",
    "\n",
    "In other words, $ Q $ is a contraction of modulus $ \\beta $ on the\n",
    "complete metric space $ (b[0,1], \\| \\cdot \\|) $\n",
    "\n",
    "Hence\n",
    "\n",
    "- A unique solution $ \\bar w $ to the RWFE exists in $ b[0,1] $  \n",
    "- $ Q^k \\psi \\to \\bar w $ uniformly as $ k \\to \\infty $, for any $ \\psi \\in b[0,1] $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementation\n",
    "\n",
    "These ideas are implemented in the `.res_wage_operator()` method from `odu.jl` as shown above\n",
    "\n",
    "The method corresponds to action of the operator $ Q $\n",
    "\n",
    "The following exercise asks you to exploit these facts to compute an approximation to $ \\bar w $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<a id='odu-ex1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Use the default parameters and the `.res_wage_operator()` method to compute an optimal policy\n",
    "\n",
    "Your result should coincide closely with the figure for the optimal policy [shown above](#odu-pol-vfi)\n",
    "\n",
    "Try experimenting with different parameters, and confirm that the change in\n",
    "the optimal policy coincides with your intuition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "This code solves the “Offer Distribution Unknown” model by iterating on\n",
    "a guess of the reservation wage function. You should find that the run\n",
    "time is much shorter than that of the value function approach in\n",
    "`examples/odu_vfi_plots.jl`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute iterate 10 with error 0.007194437603255555\n",
      "Compute iterate 20 with error 0.0004348703417873523\n",
      "Converged in 26 steps\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXxU5cHo8Wf2TGYyCWEriSRSyo5QAmKbvlCjrfiKAazVilqBLhZrXOhVuEUBKVw1WBW13kv72lhfK+7QilqXlqZ+6lI1orxFDFoBMQ4xkMnsOWfm5Nw/RofpZE9mcmKe3/ev4SzDkzOOvzznzBlMuq4LAABkZTZ6AAAAGIkQAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgtV6HcMeOHdOnTy8oKJg/f/6BAweSy30+X2VlZWFh4aJFi3w+X0YHCQBAtvQuhAcPHly+fHlNTY3X662srFyxYkVyVXV1dWlpqdfrLSkp2bJlS6bHCQBAVph69c0yf/nLXx5//PFt27YJIZqamqZMmXLs2LHEqkmTJv3xj3+cPHnye++9t3jx4vr6+qyMFwCAjOpdCJM0TauqqjKbzffee29iidvtbmpqcjqd0Wh09OjRgUAgbRdd12+++eZ//OMf7Z9t+vTpGzduTFsYi8W0Tw+3abE+DA+Z0qa1mUwmk9lk9ECkput6m9ZmsVqMHojstLhmtphNJt4ORrK4C6z5oyyWnr4drFZrty+ZtQ/jeOGFF9asWbNgwYLNmzcnF+q6nvjLdF3XNK39Xqqqbty48aGHHmq/avjw4e13OXbsmPbJe3qMEBqptbXVbDbb7XajByK1tra2SCTidruNHojswuGw0+k0m/mMoZHM+V9yCkd+fn4Pt7dYLBkOoa7ra9asefXVVx999NGJEyemrioqKjpy5MiECRMaGhqKi4s73N1kMl1wwQU9/LscDoc6bY4u+E5wQ4VCwmy25uYaPQ6pxeNxze+3Dh9u9EBkpzU3mz0eq7Uv8wdkirXNmmPOycnJyeBz9u5Xm5deeumpp57atWtXUVFRKBQKhUJCiNraWiFEZWVlTU2Nrus1NTWLFy/O4BABAMie3v1qU1tbW19fP2zYsOQSXdcrKip0XV+/fv0ll1wyduzYsrKyBx98MNPjBAAgK/r4YZk+UBTF4/EoitLD7RsbG9UclVOjxgqFQmazOZdTo4aKx+N+v384p0aN1tzc7OHUqNGsbVaX2dXza4Q9wVVfAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhADQFVVRlVbF6FEgiwghAHRlwX8smPiliUaPAllECAEAUiOEAACpEUIAxtu/b3/VD6q+9bVvTfrSpFnjZy351pLtv9uu63rqNk/vfPqSJZfMHDfzGzO+sfKyla+9/FrP1+5+fvfyC5fPnTx36klTF5+5+KH7H9I0LbGq4tSK0oJSVVFvXn/znElzppdMX7po6TN/eCaxtrSg9MP3P0w8mPXlWVk8BDAOIQRgsJ2P7vzP//jPXTt2WcyWb5/z7Wkzpv3PO//z82t/ft+99yW32bR205Urrty7Z++sU2eVjit9/unnl1Yuff7p53uy9ub1N6/43oo9b+yZOmPqad847dCHh9auWrvy+yuTLRRCXP3jq+//9f3jJ4yfOGXi66++/tPlP91661YhxPpb1g8fMTzx4Lobrxu4g4IBZEr7nSt7FEXxeDyK0tMPXzU2Nqo5qi4GaHjoUCgUMpvNubm5Rg9EavF43O/3Dx8+3OiBZEvFnIoPP/jw0h9cuvn2zSaTSQjxxmtvfPfs786aM+sPf/6DEKLuH3XfWfCdqadMffDJB0eMGiGEePO1N7937vdKTi7565t/7Xpt4qkqzqr41W9/5c5zCyHCofDqq1Y/vfPpzbdv/v4Pv19xasWH7384cvTI7X/cPnHyRCHEO2+9s+y7y4KB4O43dpeOK01scLjlcHNzs8fjsVqtRh4s6VnbrC6zKz8/P4PPyYwQgMGuWHVF9d3VV/6vKxMVFELMnDVTCOH3+xN//M2vfiOE2HDLhkTnhBBzvjZn6bKlLrcrEol0s/buz9YmKiiEcLld1XdXWyyWp554KjmGa1Zfk6igEGJm2cwrf3ZlPB5/6P6HsvyjY1DgVxsABrvwkgsTD8KhcP279e+89c4Lz76QusH79e/b7La55XNTF26+fXNP1h5474AQ4gff+4Ewpf+979e/n3x8+rdOT1115tlnbr5x88EPDvbxR8IXCiEEYLBoNHrbL2578U8vfnToI6vNOmXalJllM1956ZXkBh9/9PGo0aPM5o7PYHW91vuJVwjx4Qcftl8VCoWSj0eNHpW6akzxGCHEx0c+7v1Pgy8eQgjAYDesuuHJR548Y8EZ625eN/+M+Tk5OUKI39f8PrnBqNGjjn16TNf15LlTIURbW1tbW5vFYul6bfHY4g/f//Ddj991uV1djOHTxk/Hlo498cejnwohRo8ZncEfE4MW1wgBGOy5Xc8NKxx23/b7zjrnrEQFQ8FQ6gbjJ4yPRqNv172duvDay68dP2J8w5GGrtdOmDhBCFH3el3q2qPeo2tXrX3kvx9JLqn9c23qBruf3y2EGP+V8Rn5ATHIEUIABsvz5EUj0YA/kPhjJBJZf/16IYQW/+z2huWXLxdC3LTmpubjzYklb9e9/exTz5aOKz2p5KSu1/6o6kdCiLWr1h7812cX/KLR6Jqr1jx0/0NOpzM5hru23PX+e59dMty7Z+/dt91tsVguXn5xcgO+bnQI49QoAIOdv/T8e++4d+E3F5654MxQKPRy7ctFY4sKhxcePnh4w5oN1665tuKsivMvOv/JR548ffbps0+bHY/FX3npFU3Tbtl6ixCi67Vzvz73h1f88Lf/77dnlZ81a/asgsKCN1978/ix4wuXLKw8vzI5htJxpQtPXzh77uyYGttTtycei1+9+urxE8cLIVwulxDiiuVXjB4zetNtmww6SMgi7iNEV7iPcDAY8vcRxmPxbXdve2L7E94G74TJExacu+CKa67Y/sD26o3Vdrv9mb89U3RSka7rjz/0+BPbn9i/b78jxzFtxrSf/fxnM8tmJp6h67VCiOd2Pff7mt+/+893VUX98le+vHTZ0gsuvsBqswohErcJHjh64Nabbq39S21TY9O0GdMu+9FlC5csTOz71xf+un71+k8aPsnPz3/9vde5j9BY2biPkBCiK4RwMBjyITRW8n75brfkhvrBgBvqAQDIMEIIAJAaIQQASI2T3QCk9tjTj8W1uNGjgJEIIQCpjRw90ughwGCcGgUASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqD+vaJf/3rX1ab1eVy2Ww2j8dj9HAAAEPQoA7hRRcsbWw6KhL/6LRJ2O12p9PpcrnsdnteXp7D4cjJyXG73Xa73e12O51Ou93u8XhsNpvL5crNzU38scO9DP7BAACDxqAOYVH+OHssNxwOxmJqMOBXhaqaVL/JL4QQJpEMZO/+aBJCiLy8PLvd7nK5kvm02+25ubnJ2Wcin8nK5uTkOByO5F65ubk2my2zX38OADDEoA7ho4/+XddPXMVUlFZFiYZCgVhMDYUCra1RRWkNBv2qqkQioWg0HIupfr8vFlOj0fBn+Qz6FaW1tTWxlxIKBVujEVVVgoGgEOK46Xiv8tn+jzabLZHPxHSzfS+7mJWmxTj7hxMA0IFBHcI0DkeOw5Hj8Qzr/1MFAi2xmBqJhCKRUCymBgItqqq0tkbC4aCqKqmVjcXUcDjY2hpRVSVZ2c/28rfERMxv8vd6ktrR8j7PStNi3P+DAwBS+SKFMIM8ngIhxPDho/r5PMlJqqoq4XAwGo2oqhII+GIxNRL5t1lp+82CwZZYTA2HQ5FwKB6PBQIBIfoyK037Y98unaZVlkkqAHlIGsJMydQkVdf1YLBFVZVoNJLMZ2trVFVPzErTKpsyl02e+1VDoYBiUhShZGSSmp+fb7Vac3NzE1ltf1K3s1lpWoz7eXAAIKsI4aBgMpkycspXCJHIZ8q533AspgYCvrTK9mSS6g/4RdqV1D5ltT+z0tTKWq385wog8/g/y1CTk+PMyXFmapLa3Hw8FlPjcfXzDyglZqWJyWs4EgmpqhoMtqR8IkkNhQKKElWUEzGO+qNRU/TEJFX08QywyWzyeDxp080eXjpN7uVwOFwuVz8PDoChhBCiY4lJqtlss1gsTqezP0/V/gO9qqqGQic+0JuYlSYr2/5zTIm5rKbF/f52p3wHcJKatpnb7bZYLP05MgAGA0KIrHM6XU6nq/+TVE3TwuHuP9Dbo0lqIBoV0e6vpHZXWYvF0udZaVqM+3lwAPQZIcQXhsVi8XiGZeRianKSmnLbzIlJameVTf0ck6K0hoJ+rU3z+/z9P+srTKKLWanNZktcLu3w0mniZptkZZmkAr1FCCGjzE1S46kf6FXVlNtMAy2dfAXEZ59ISjljrESjkUggEhGRFlOLEH069yuYpAJ9RAiBvrNYrJmapKZ9oDf1tpnW1mhT01Gr1drtV0BkdpKa9t0OXVw6TVQ2uZnD4UiNsdnMv3KDQY0QAoOCy5UnhCgoGN5+laZpfr+/sLCwJ8+TmKSm3DaTvHSqpN4205NJajgQFkL0fZL6+QOr1drnWWlajPt4cIEuEUJgSBmYSWrirlNVVRLf7ZCobOK7HVI2C6uqGvD74iLub+58ktqbc789v+s07dJp6s02Ho/HZDL1//hgyCCEADrWxSS1Vz7/qFFIVZXU22Y6/EBvZ18BkahsOBAOi7DP5OvtpdO0B32bpLa2to4YMSKxPZPUoYQQAsgum81us9kzMklN+0Bv8rsdwuFQLKamfk9h4tJporLZnqSmXRNlkvqFQwgBfGHk5eULIQoLR/bzefo2SY1EwrGY8nl0/22S2mxqzuwktbN/h5wrqdlACAFIp2+TVJ/Pl5eXl/adtz2cpKZV1pBJag+/JV/CSSohBIC+G5hJarsvIAwnP6/U2ZXUgZmkpp0K7uzfVe3nwck2QggAxhv4K6ntJ6mqqkQiocEwSU1GN1HTbP8TNIQQAIaUQTVJTVS2v5PUlD9evvzyLRu39PNHS0MIAQAdMHCSmvz3xuPxWOJO1tbWaDgU1LR4W8TR//GkIYQAgOzK1CS1ra3NlPiqo4ziOwABAF8MZrM5JyfzH70hhAAAqRFCAIDUCCEAQGqEEAAgtV6HUNO0yZMnpy0sLy83fW7lypUZGhsAAFnXu9sn7rrrru3bt9fX16cu1HW9vr7e6/W63W4hRJbu/AcAIBt6NyOcMWPGunXr0hZ6vd54PL5w4cIxY8ZceumlgUAgc8MDACC7TLqu93of07/tVVdXd/31199xxx0lJSWrVq1SVfXhhx9uv5eiKG63e968ee1XTZ8+/aabbkpb2NTUFAyO0nWuYhopHA6bzebB/525Q1s8Hg8Gg8OGZeALPtAfLS0tiW+/NHogUrPbI/n5isfj6eH2+fn5Foul620yEMJUXq932rRpzc3N7VcpipKXl7djx472q4YNGzZnzpy0hZ9++mkgMJIQGisUClksFkJoLE3TAoEAITRcS0uL2+3m6o+x7PZIQYGan5/f4+3t3f57Uhl4Revq6hRFKS8vT/yVDkenXwRnMpnOPffcHj6t3W632+2E0Fh2u91isdjtdqMHIjVN02w2G6+C4RKvAiE0lsOhORymLkLTB/3KTG1trRAiHA6fd955+/fvV1V106ZNS5YsyczQAADIvn6FsKKiQggxb968devWVVZWFhcX+3y+6urqDI0NAICs68s1wr5RFMXj8SiK0sPtGxsbuUZoOK4RDgaapvn9/sLCQqMHIjufz5eXl8epUWM5HNFeXSPsCTIDAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQTSnX325EmTTJl6NlVVFKU1U88GIOMIIZBdixbNnDHDafQoAHTKavQAgEHnzjsfYQ4HyIMQAummTPmq0UMAMHA4NYqhLHG1T9f1u+/eUFbm+c1vbhVC1NY+c/nlC+fNKy4r81xwwWmPPPJrTdPa75W6pNtd/vSnx1as+PbcucPPOGPc1Vd/9/XX/5ZYPmmS6eDB+sSDDF53BJBBhBBD32OP/de99/7CZDK53Z7bblv9k5+c+/bbr02Z8tVTT51/+PAHGzasvOqq89PClqrbXW655WfXXvu9f/7zza9+9WslJeNffHHnsmVnvPjiTiHE2rV3FhaOTDxYu/bOgfl5AfQKp0Yx9N1++88feOAvp51W8dZbL2/ceOU3v3nOHXc87HZ7hBDhcPCGG370pz899thjv1m69Ir2+9bV/f2++27rYpc9e1753e/unDLlq/fd99yIEaOFEG+99fKll57+y1/+729/+7xly659+OFtzc1Ny5ZdO8A/NYAeYkaIoW/Vqv/zta+dYTKZfvvbXwohbrhhayJpQgiXK2/z5v+yWCxPP/1wh/t2u0tNze1CiLVr70xUUAhRVvaNCy/8scuVF42Gs/yTAcgAZoQY+r7+9TMSDz74YJ8Q4vLLzzWZ0i/XffDBux3u2+0uH3zwrs1mnzNnfuqqm276v5kYOICBQAgx9J100rjEg6NHPxZCHDp0oP024XCww3273aWh4dDIkWPMZk6uAF9UhBBDn9VqSzwoKio9eLD+rbcCLldeD/ftdpeRI8ccP96o63rqlLGtrU3X28xmS/t5JIDBhl9jIZHx46cIIfbseTV1YWNjw4YNKx9//L6+7TJu3KRoNLJ37+upG1x//aVTp9o++eRwZscPIBsIISSyYsXPhBAbNqw8fPj9xJJoNHLDDT965JFf5+Tk9m2X73//KiHE5s1X+3zHEhvs3fv6888/WVIyvrj45OTz8FU1wKDFqVFIZM6cecuWXfvAA1vPPfeUmTNPKygYXlf39+bmprPPvmDhwov6tss3v3nOkiWX/eEP/33WWRPLysrj8dirr+5ua9N+8YtfJ54hN9cthLjmmgvGjh1/ww1bB+yHBdBDJl3XB+ZvUhTF4/EoitLD7RsbGwOBkbrOnNVIoVDIYrE4nV/U74w+++zJBw/W19f/23/kL7yw45FHtr333juqqpx88sQLL/zxd76zPHkdscO9ut5F1/UdO+7fufOB+vq9dnvO1Kmzrr564ymnnJpY+7e/Pbtp01Ve75GCgsKXXz7ah59C0zS/319YWNiHfZFBPp8vLy/PamX+YCSHI1pQoObn52fwOQkhuvJFD2HfdJhPAxHCQYIQDgbZCCGZAdKpak9/XQMwBPCrDXCCpml1dX9vaDjkdLqMHguAAUIIgROuu+6SZ599VAgxf/7ZRo8FwAAhhMAJxcUnz5r19bKyb6xcudbosQAYIIQQOOG66241eggABhoflgEASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCk1pcQapo2efLk1CU+n6+ysrKwsHDRokU+ny9DYwMAIOt6HcK77rqrvLy8vr4+dWF1dXVpaanX6y0pKdmyZUvmhgcAQHb1OoQzZsxYt25d2sKdO3dWVVU5HI6qqqodO3ZkaGwAAGSdSdf1vuxm+rcd3W53U1OT0+mMRqOjR48OBALtd1EUxe12n3rqqe1XnXLKKbfcckvawqampmBwlK5zFdNI4XDYbDY7nU6jByK1eDweDAaHDRtm9EBk19LS4nK5bDab0QORmt0eyc9XPB5PD7fPz8+3WCxdb2Pt96iEEELXdZPJlHigaVpnm5lMpk2bNrVfnp+f73K50haGw+F43EkIjaVpmsViIYTG0jQtHo/zKhhOUZTc3FyrNTP/20Tf2O16bq6lfTI6YzZ3H5HMvKJFRUVHjhyZMGFCQ0NDcXFxZ5uZTKYzzzyzh89pt9vtdjshNJbdbrdYLHa73eiBSE3TNJvNxqtguMSrQAiN5XBoDofJ4XBk8Dn7m5na2lohRGVlZU1Nja7rNTU1ixcvzsC4AAAYEP0NYUVFhRBi/fr1e/fuHTt27L59+2688cZMDAwAgIHQxzl+8pMyiQcFBQXPPPNMxgYFAMBA4QocAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1HoXQp/PV1lZWVhYuGjRIp/Pl7qqvLzc9LmVK1dmdJAAAGRL70JYXV1dWlrq9XpLSkq2bNmSXK7ren19vdfrDQaDwWBw69atmR4nAABZ0bsQ7ty5s6qqyuFwVFVV7dixI7nc6/XG4/GFCxeOGTPm0ksvDQQCmR4nAABZYe3V1g0NDaWlpUKIxLwwudzr9c6ePfuOO+4oKSlZtWrVNddc8/DDD7ffva2traysrP3ymTNn3nbbbWkLm5ubg0GrrnMV00jhcNhsNre2tho9EKlpmhYMBk0mk9EDkZ3f74/H41Zr7/63icxyOKKa1hqLxXq4fUFBQbcvmUnX9Z6PwOVyHT9+PCcnJxKJjBw5MhwOt9/G6/VOmzatubk5bbmiKHl5ebW1te138Xg8kyZNSlvY2NgYDI4ihMYKhUJmszk3N9fogUgtHo8HAoHCwkKjByI7n8/ndrttNpvRA5Ga3R4pKFDz8/N7uL3Vau32l8je/WpTVFR05MiRCRMmNDQ0FBcXJ5fX1dUpilJeXi6EsNvtDoejw91NJlNim56w2WxWKzNCg1mtVovFwq/AxjKZTFarlVfBcFarNfH/JaMHIjWbzWaz6Zn9daR3mamsrKypqdF1vaamZvHixUKIxAwvHA6fd955+/fvV1V106ZNS5YsyeAQAQDInt6FcP369Xv37h07duy+fftuvPFGIURFRYUQYt68eevWrausrCwuLvb5fNXV1VkZLAAAmda7a4T9oSiKx+NRFKWH2zc2NgYCIzk1aqxQKGSxWJxOp9EDkZqmaX6/n2uEhvP5fHl5eZwaNZbDEe3VNcKeIDMAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEiNEAIApEYIAQBSI4QAAKkRQgCA1AghAEBqhBAAIDVCCACQGiEEAEht8IZw1apV0WjU6FHI7plnnnnppZeMHoXsjh49es899xg9Coht27Z99NFHRo9Cdi+//PIDDzyQ2eccvCF89tlnY7GY0aOQ3YEDBw4dOmT0KGQXCoVeffVVo0cB8cYbb7S0tBg9CtkdOnRoz549mX3OXofQ5/NVVlYWFhYuWrTI5/N1uxwAgMGs1yGsrq4uLS31er0lJSVbtmzpdjkAAINZr0O4c+fOqqoqh8NRVVW1Y8eObpcDADCYmXRd79UObre7qanJ6XRGo9HRo0cHAoGulyfF4/GTTz752LFj7Z/T6XQWFRWlLfz4449HjTrFbLb2anjILL/fb7FY3G630QORmqqqLS0to0aNMnogsmtqasrPz7fb7UYPRGqRyLF4/HhhYWEPt3/iiSemTJnS9Ta9zoyu6yaTKfFA07Rul5/4m6zW55577sCBA+1X5ebmFhcX93YkAAB0bdy4cd1u0+sQFhUVHTlyZMKECQ0NDan16mx5qunTp0+fPr23fyMAANnT62uElZWVNTU1uq7X1NQsXrxYCFFbW9vhcgAABr9eXyNsaWm55JJL3nnnnbKysgcffDA/P99kMum63n55lkYMAEAG9TqEAAAMJUZ+s0zP783nbv0s6eLA7tixY/r06QUFBfPnz09+xKm8vNz0uZUrVxox5KGpixei/THn7ZAlnR1YUzu8EQaApmmTJ09OXZK9NBgZwp7fm8/d+lnS2YE9ePDg8uXLa2pqvF5vZWXlihUrhBC6rtfX13u93mAwGAwGt27datzAh5rOXogOjzlvhyzp7MAGU6xbt+Retc4AAAKtSURBVG716tW8EbLtrrvuKi8vr6+vT12YxTToxpk4ceL+/ft1Xd+/f//EiRO7WN7Zluinzg7sn//855/85CeJx59++unw4cN1XW9oaPB4PGVlZW63e/HixY2NjYaMeUjq7IXo8JjzdsiSbg/s3r17zzzzzMOHD/NGyLbdu3fv2rUrrVDZS4ORIXS5XJFIRNf1SCSSl5fXxfLOtkQ/dXtg4/H4ypUrf/rTn+q6/uabb1ZUVOzZs+f48eOXXXbZRRddNNDDHbo6eyE6POa8HbKk6wOrKMrcuXP37dvHG2HApIUwe2kw8ntb9B7fm9/Zluinrg/sCy+8sGbNmgULFmzevFkIMXv27N27dydW3XrrrdOmTRvg0Q5hnb0QHR5z3g5Z0vWBvf322+fOnTt16lQhBG8EQ2QvDUZeI0zcgy+E6PDe/NTlnW2JfurswOq6vnr16k2bNj366KO33nqr1WoVQtTV1b3yyiuJDex2u8PhMGTMQ1JnL0SHx5y3Q5Z0cWA1Tdu2bds111wjeCMYJ3tpMDKEPb83n7v1s6Szl+Cll1566qmndu3aVVRUFAqFQqGQECIcDp933nn79+9XVXXTpk1LliwxdvBDSWcvRIfHnLdDlnT2Kgghdu/ePXbs2K985SuCN4IRsp6GPp9U7T+fz3fOOecUFxdXVla2tLTon58Rbr+8/RJkRGcvwU033dT+v5O2trZ77rln/PjxI0aMuOyyy/x+v8GjH0I6eyE6POa8HbKks1dB1/WLL75448aNice8EQZM8vhnOw3cUA8AkJqRp0YBADAcIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASI0QAgCkRggBAFIjhAAAqRFCAIDUCCEAQGqEEAAgNUIIAJAaIQQASO3/A9ZjWHidLjN0AAAAAElFTkSuQmCC"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp = SearchProblem(π_grid_size = 50)\n",
    "\n",
    "ϕ_init = ones(sp.n_π)\n",
    "f_ex1(x) = res_wage_operator(sp, x)\n",
    "w̄ = compute_fixed_point(f_ex1, ϕ_init)\n",
    "\n",
    "plot(sp.π_grid, w̄, linewidth = 2, color=:black,\n",
    "     fill_between = 0, fillalpha = 0.15, fillcolor = :blue)\n",
    "plot!(sp.π_grid, 2 * ones(length(w̄)), linewidth = 0, fill_between = w̄,\n",
    "      fillalpha = 0.12, fillcolor = :green, legend = :none)\n",
    "plot!(ylims = (0, 2), annotations = [(0.42, 1.2, \"reject\"),\n",
    "                                     (0.7, 1.8, \"accept\")])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next piece of code is not one of the exercises from QuantEcon – it’s\n",
    "just a fun simulation to see what the effect of a change in the\n",
    "underlying distribution on the unemployment rate is\n",
    "\n",
    "At a point in the simulation, the distribution becomes significantly\n",
    "worse. It takes a while for agents to learn this, and in the meantime\n",
    "they are too optimistic, and turn down too many jobs. As a result, the\n",
    "unemployment rate spikes\n",
    "\n",
    "The code takes a few minutes to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute iterate 10 with error 0.007194437603255555\n",
      "Compute iterate 20 with error 0.0004348703417873523\n",
      "Converged in 26 steps\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 20\n",
      "date = 40\n",
      "date = 60\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 80\n",
      "date = 100\n",
      "date = 120\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 140\n",
      "date = 160\n",
      "date = 180\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 200\n",
      "date = 220\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 240\n",
      "date = 260\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 280\n",
      "date = 300\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 320\n",
      "date = 340\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 360\n",
      "date = 380\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 400\n",
      "date = 420\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 440\n",
      "date = 460\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 480\n",
      "date = 500\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 520\n",
      "date = 540\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 560\n",
      "date = 580\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "date = 600\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOy9eXwb9Z3///qMTtuS79vxkQtygJNwlgCF0FAKJQlH6bbf0mO/u7Tfttku7W7LbqF0t8uXbNl+F/ot/VFYCG350WtbwlEopUASSigkHMHkPp04vmXZlqxbM5/vH5+Z0eiwLDty5LHfz0f+UEaj0UdjaV7zvhnnHARBEAQxV5EKvQCCIAiCKCQkhARBEMSchoSQIAiCmNOQEBIEQRBzGhJCgiAIYk5DQkgQBEHMaUgICYIgiDkNCSFBEAQxpyEhJAiCIOY0JIQEQRDEnGZGCOHIyMiJEydy3Dkej8/ltnCc83g8XuhVFBJ53z7u9xd6FQWDvgB0BZjjXwBZlhVFye8xZ4QQvvTSS9/85jdz3HlkZGQufw+i0ajP5yv0KgoJ/+u/jr/7bqFXUTBisdjo6GihV1FIfD5fNBot9CoKRjweHxkZKfQqCsl0fAFmhBASBEEQRKEgISQIgiDmNCSEBEEQxJyGhJCYzQyECr0CgiBmPCSExGzmV8fynF1GEMTsg4SQmM38oUvxhAu9CIIgZjYkhMRsZjSK358ko5AgiGyQEBKzmUAcW07M3eJrgiBygYSQmM0E4njplOKPFXodBEHMYEgIidlMMM7DMv54iryjBEGMCwkhMZsZiwHAlk7yjhIEMS4khMRsJhgHgBe6FIWkkCCIcSAhJGYtoThkDgAjURz3kxISBJEZEkJi1hIwDCl530tCSBBEZkgIiVlLIJ4Qvw4SQoIgxoGEkJi1GC3CDm/h1kEQxMyGhJCYtQRiAFBkBYCDI2QREgSRGRJCYtYyFgeAlhIGwBshISQIIjMkhMSsRViEzS4AGI4Wdi0EQcxcrIVeAEFMFyJZpsLOiqw8FEcwjmL6vs86tm/ffvvttxd6FWcOzrmiKBaLpdALKQBr1qz5z//8z+k4Ml0YiNlJMI4d/RyAy4YKOwvF+UiUF1tZoddF5JnR0dGysrL777+/0Ashppdt27Zt3759mg5OQkjMTl48pYgEmRIrKhzoCWI4gsbiQi+LmAZKS0tXrVpV6FUQ08uJEyemTwgpRkjMTrZ08q4AAJTYUOEAgOFIYVdEEMQMhSxCYhYSU/B8lxJXAKDEyirsALg3wgFyjRIEkQoJITEL2drLdfvPpVuElDhKEEQmyDVKzEKeOZEYQFhsJdcoQRDZICEkZhscePZEonzeZUWFnQEYppp6giAyQa5RYlbx3Eml0sFOBRKaV2xFJblGCYIYHxJCYlZx61Y5pWreZWPlDg5ghFyjBEFkgoSQmD2E4vDF4IslbawtgkgW7Rwj1yhBEBmgGCExe+gPqVJ3aV2iTKKhiK2qYgx418PjyjivJAizwRjVAuUNEkJi9tAXAoBl5ezHl1okBgA2CdVOVDsx380CceylYUwEMSXWrl07TTvPBEgIidmDsAgXlbEVlazMDgD1RUwo4kW1DMCuQRJCgpgKr7zyyjTtPBMgISRmD/0hAKh1AkClgwGo15qLnl/NAHzgJSGcExwa5buHztC/UDz13Y1OS/GYMfarX/1qxYoVVVVVDzzwAIC+vr5bbrmltrZ2wYIFn/vc53p7e8Vu3/72txsaGr73ve/9y7/8y1lnnVVWVnbvvfeKp+67777a2trLL7+8s7NTP77X6/3c5z7X0NDQ2Nj4+c9/3uv1fvGLX9RbkH/pS1+67777shx2vGWkrPaGG24AsHLlSuPnevLJJ2tqap599tmVK1eWl5c3NDT84Ac/MO6c8eAzE0qWIWYPwjUqxK/SgaNAY7F6SVrgBoDOsUItjTijfOpV+b2hM3TT884N1vOqJw7XnTx5cvfu3Vu3br3++utvv/3222677dZbb/35z38eiUR++MMf3nbbbb///e8BnHvuuS+99FJ7e/v999+/Z8+e7du3r1+//tvf/jaAsbGxnp6eu+++++tf//qWLVvEYW+//Xa73X7s2DEAX/va1/7hH/7h05/+9D333PP1r389Go0+9dRT77777h133DHeYcdbRspqn376acbY7t27jZ9ox44dW7du/cxnPvPZz37261//ekdHx4c+9KF//Md/1Hdet25dxoPPQEgIidmDcI3WFTFArR1s0CzC+W4G4LifLMI5wdnlTDlTf+och1x++ctfZoxdddVVoVAIwNatW42qUFNTIx7cfPPNNpsNwMaNG61W69q1a8PhsHjqC1/4gtVq/cY3vnH22WfrL3zhhRf27dtXVFQE4J577mlvb3/kkUc+85nP9Pf3v/322+ecc05zc3OWw463jJTVZuRf//Vfa2pq3nvvvV27dv3sZz/bvn17NJpUqzvewWcgJITE7EG4RuuKANU1yuuL1Fv1NhcD0ElCODf45ZoZMbfW7/frj91ut/GpioqKjo6OBQsWAAgEAkNDQ2K73W4XD6xWKzKlhkqSJMuycYu+D2NMlmWbzbZu3bpnn312+/btt956a/bDjreMlNVmRAjbJz/5Sbvd/ulPf3rTpk0///nPc/mMMxCKERKzh4EQB1DrZAA+v1j69wst17doP3gHyu3wx+AJF3KFxFzA4XBs3bqVc/7QQw+Nt89NN920adOmYDA4MDCwYcOGTZs2TXjYxx9/PB6P33///R/+8If1jddee+2dd94ZDodDodCdd9553XXXAfjEJz7x5JNPvvjiizfffHP2Y05qGbFYLH3jn/70pzvvvPP6669/8cUXAcTjcX3nKXzGQkFCSMwejBbhNfPYHSuk8w3BmzY3ldUTZ4J77rnn5ptvbm9vr6ury7KPLMvz589fvnx5W1ubSDPJTiwWq6+vf/XVV3/4wx/qGx944IFQKNTW1rZgwYJoNCpyWz7ykY/s3r37yiuvLC8vn3CpOS7juuuuW7hwYfr2e++998orrzz33HOHhoauueaaL3zhC/rOU/iMBYPPAH7961/fcsstOe48ODgYjUandT0zmXA4PDQ0VOhVFJLYRRdFt23L+FT1E1H8V3QglPmFN7wUx39Ff3FEnsbFTT+RSMTj8RR6FYVkaGgoHA7r/33mmWfWrVtXwPWcGSZ7rV69evWWLVumaTEFYcuWLRs2bOCce73eUGicH/lUIYuQmCVwYCQKBlTYM+9wcS0DcO9uJUr9ZYjZSywW271799GjR6+99tpCr8U0kBASswR/DHEFLhus43ypbz9HWuBme4b5633kHSVMxi9/+csc93zuuec++tGPPvjggw6HY1qXNJugrFFiliDGDYo6+ow4Lbiolh3z8/4QF224CcIsfOpTn8pxz5tuuummm26a1sXMPsgiJGYJYgB9RdabYJFH0zduZRRBEHMREkJiliDm7o4XIBSIWnt9SAVBEARICIlZg3CNlo/vGoVmEYoqi5iWMmMcZ08QxByEhJCYJfzqKMdkLMK9w1yM8H3xVKoQRuS0VxIEMXshISRmA10BvuWEAq3F6HgYLcLdQ7wvyENxbO9NFcIOLycjkSDmDiSExGzgJ/sVMX0+d9foB8O8N4iDo/yoL1X1eoK8N0hSSBBzBSqfIExPRMajB9WIX3bXaK2TMWAgxDnQ4eUXVvPeYIama54wOv1oLM54DGJmYbfbDxw48LWvfa3QCyGml2PHjomO4dMBCSFhep4+oQxoFRHZyyccFpQ7MByBN4IOL+8NYTTK+4IIy3AaxhUMhiEPc4cFxVYsLaeKwxnNxz72se9973sDAwOFXsgZQpblcDhcUlJS6IWcaRYtWpSx2WleICEkTM8RX+JxRVbXKIBqJxuO8H3DfCCE3iA/MQYOdPr5EoPgecL85W7+pdf50nL26set9UXTtHAiP+RebD4LiMVio6Oj1dXVhV7IrIJihITpGQpzaDN4s7tGAZTbAeC1Pg6gN4gDIxxA5xjChkxRTxjvDnEAg2H+zIlJdCY9Y8NgCYLIIySEhOkZigDAFxZLmMg1Cl0IexUA3QF+aJQD6PTzV3oSIuYJc9GnZjiCrT1cHD8X9o+QEhKE+SAhJEyPsAhX10nLytmErtFyOwPwxgAH8LaHC0Owc4xv6UxYfh5N+WSOl7uVnpwr7t/oJyEkCPNBQkiYHqFb1U7c2MYmdI0Kk3EsBgB+beB2px87+nlvUP3voKEZ6VAE3UHkiNBXgiDMBQkhYXq8EQCodOCzi6XxZjDplGdSyuN+3hvkb3tUo3AokqRn3TlbhDsHuExSSBBmg4SQMD3CNVrlZGeXTVzqkLHifv8IH41i5yAHEFXgiyY9m7tFOBjmh0dJCQnCZJAQEuZG5hiJQmIT54sKMlqEwke6c4ADGAojRcpytwjH4ujwkhAShMkgISTMzXAECkeFHVJuhe8pQmiX0FyivnKXh3PAE05VshyFUOYIkxAShAkhISTMjYjnVTlz7f8iskYBVTjbK9mCUvWp4QiOjHJPWrGEcI16whOIXCAODnR4c105QRAzBBJCwtwMhQGgaqLyQR290HBZOQNwUS1rKk6I6M5BPpg2tldYhD1Bfsmz8V8cHbe+PhADQBYhQZiP0xLC4eHhdevWVVZWrl+/fnh4WGyUZXnJkiXZ9yGIfCGGC9YU5W4Rqg8uqWMALqphTYaujTsHuSec+hJPGBEZviiCcXxmq/yfH2TWwrEYB3ByjI9EMz5PEMQM5bSE8Pvf/35ra2tvb29LS8t9990H4Ic//OHq1asPHjyYZR+CyCOdYwDQ6sp1fz1r9JJaTQgNFuG7Q6pr9O+WS1uuVvtwc6AnqE7xBXDn2/KhTKmhgbi6MxmFBGEuTqvp9pYtW5555hmHw7Fx48YNGzZs2rSpvb194cKF69aty7JPxkMdO3bs3nvvTd++bt26+fPnG7cEAgGLxWKz2U5n5eYlEokEg0G7PbcUydmITZYjoVBkbEz897DXClgabNGxsZzmytsVfG+FVWJY6Iy6rbYmS6DKIgHq16ljSFnujgOWZkd0kVNhsC9w86N+dsQTGghC7BaW8bP9kX8+J55y5EGfepxdveHzXNM15D4ajQYCAYcjZ1/wrCMQCHDOY7HYxLvORmKxWCAQcDqdhV5IwQgEArIsx+OpP8DxKC4ulqQJTL7TEsLu7u7W1lYAwuYDsGbNmgn3yUgkEsn4bCgU4snTwrnG6azcvMzxjy8wnoHOMQagpUTJ8ZzYGG5fEgPQE8R5VZyBNxSpL7RJ8MXYe14GoMrBK+18voufU86P+tlhP8TgXwZwwBMB57wnmDSzMKD9MPeMsOn7A9EXYI6fgTn+8TH5M5DLnqclhJxzxph4IMuZb4Fz2QfA0qVLf/SjH+XyppFIxO12z1mL0G63S5LkdrsLvZCCEbdYiouLbdoZOBWOA3xZbbHbPbnBgQtLcGmD7HY7FjMOxCWG1XVsey9/d1gCMK+8qLGCXdEoX1zLnumSf3HCfkObBMgNxawnyEdlm9vt/O8jyl2rEneaskcBZABHA9Y9IYdwvQJQeK6lHbkQjUYBzOUvQCwWKykpmbM2cSwWUxRlLn8B4vF4UVFRfm3i04oRNjY2dnV1Aeju7m5qapryPgQxZTr9HECba9JSY5Xw8WYJQEMxkxhqnLigmgGq5TffBQDXNbNbF0nlduzo57sGOQBRa+EJ86EInjL06e4LYVTz1fUE8eiBxFM/PTyJQU4EQZx5TksI161bt3nzZs755s2bN2zYAGDbtm0T7kMQ+cIThj+GMvvE05cy8qFaBsAmodaJxmJ2ZYP6c2gqYYvLGICrm6QSKy6tYwCEEC50M/G+vz+p9BkKLU74+R4tR6YnyH9zXIkqAKBw3LlL9s/ReBZBmIPTEsK77767o6Ojubl57969d911FzLFCNP3IYh8ccTHASwsnaLnUfdYNpawphKsbWJuGwBc2aA+UWYHgGYXA3BiTFiEDMBgGFs6ucfQjK0/xN/2qP8LxTEWU2cTnhzjfSH8bKpGIU36JYgzwGnFCMvLy59//nnjFhGWNAYn0/chiHwhxGZJDr22s9NUzBqL4bTg2mbpN8eUNQ1JBxQ92IQmLXADgCfMX+rmMQWjUbUwsT+E94aSVKvDy1dUsgOjAPDjfcpXl00lUPjcSWVDK3W9IIjphX5jhIk5IISw/LSFsARNJQzAja0MBotQ0GwoUqxxMrcNUQWhOGBoTNofUmcc6nQMcX2FB0b4K92TNu68ETxxhExCgph2SAgJEyPsraXlp3ucpmK1v8x1LdLiMpbia20pSfy31I5qQ19TvQ1Nf1pjNlFWL2xWAL/rnLR39AMv396rCEt0jKKMBDFtkBASJiafFmExA1Bqw7+dn/qjMFqEZXZUG9K2DUKoPtC7ngohPKAJYV8Qo5Nsvfa+l3vC2DPMx2LUrYYgphESQsKscOCYnzNg8VSTZXSaShIdRz+5IPVHMa+E6fG9UluKEOquUfXBIi1g2RfCQChhEQ6G+ZZJGoVC/Lb28KN+fmqiUVBbe0kpCWKKkBASZiUcR1yB0wqH5XQP1VQMveNouqjapUSdYqmd1TgZtOEVg5pF2KdZhEZVfrVH0XcYCCPL5IqMqELYy4+McjEKKguPHlCoSIMgpgYJIWFWQjIAFJ22CgKY72bZKxHXtTIAEoPLimonFrjZZxZJyGgRljIArS4GqFUT4nF3gG/t5YFc+yNC5tg7zAG81qcc8k08HPjAKH+hiyr3CWIqkBASZiUU5wCKrHloX1Y8URnRja0SAJcVEkO1k31lmVRXBGgxwrCciP8tKgWAb5wrAfhTNwewqorZJATjiCvoDebqwDzi48E4AAxH8NRxZUKLsNPPt3SSd5QgpgIJIWFWhE7kxSKckMvqWY0TpXbR4Bv/8yxJRArFzCZjyujiMlbtxOcXSxYGmQPAfHcirNg7kZ7pGLNj3vHw8SzC3xxT/mmX/Kdu7o3gD11KZKpDL474eChna5UgZhkkhIRZEa7RCY25vGBhWNcildoA4Ob5UoVDLaIQrlE9ZRRAi4tdWieV2XFetWqqznezWq3iIneL8MBI4jEHugOpOwTiOObnjx1Uvv++8vABBYAvhld6pmgUvuPhkw1hEsSsgYSQMCvCgik6I0II4MY2qdQOAE4LANXIE65Ro0XotmF9KwOgt6dpc6OmSH22zyCZ2Tnu59Ca2gDoCabOknm6U1n46/hL3RzAn/tUDcsxMTXdMD3iw4/3kRAScxQSQsKs5DFZJhfWNrFGwyz7SgcDMJRsETKg2IrrmyUAegvvNtdULEIxVeOSOvWFYRlD4aQdjHk3A9oCRGfwCTk4ysXxdY74+HtDfEc/RRmJuQgJIWFWzrBF6LTgrxYYWszYAEBULOhCWGIDA2qLAOCyemaTAKDNzXSLMPcYYecYAOhDDQH0JItoxgTUkdxq9r0RvnMwVQgBPEhGITEnISEkzEpQZI1a8jf0diLWG/pfOyxwWBBVEJZV12i5HSUGVXbbcEE1q3bCbcNkLcK4gq4AlxguNghhSpgwkKlqcDSa0/G9EaQK4SgH8NRxRQzZIIg5BQkhYVbOZLKMwJnshnVrRqGQqA83SC5bkipf2cBEJX6NljWaS4xQ5jgV4HEFDcXMOHC4O9UizKBYvlhOk5tShHAsphq1UQVf2ZFIPP3JfuX/fKB4IxMfkCBMDQkhYVbOsGs0HbeNAfBF+btDHMAt81lJ8mLWNEptbgagtghLyxlyswhvfln+2WEOYL4btUXQLd5UizDZNWphsEpQOHLpL+ON8Hc9PK4gpsATxlE/58C8Euaw4IUuLo7AgX94S/7Ht+TDo2QjErMcEkLCrJzhZJl0RJjwqB+dfl5qw4ZWSdiIOqvr2NllAFDjZJ9cwKwSPGHEsobhIjJePKX8y7sygMvrmIWpEUekWYTGeRRVDvzhY9aGIgZgJAfvqDeCQBx7R/iuQf62hwu/6KoqVukAAH+MA+gJqBX91LmNmPWQEBJmJVhwi9AOAC93KwDOr2Zum2r26ZRY1RbetUX4SKNU62Qc6DPUWqR7Sju8XBTFt7rYt1daANQXqcdMqanXLcJWF3v7BuvVTazMDiCnGRfC27lzgG/t5R1efsQHAItKUarauAAgNgIYy+SDJYjZBAkhYVZCZzxZJgVhEW7r5QAuqmUALqhJXUx7JQPQUMwurmX1xUBy9b1e/6fzx1McwIpK9pPLLC4bADQUQ8y+ODmWtGcgpurTedVMOGDLHUDWxFG9d8xwhAPYOci39SodXi4SZBaUMrchFVbkkYIsQmIOQEJImBXVNVroGKEYtDTfzQBcUJ1ZlUussEsQ7Ul1IQzFk8r+wjL++jX5O+/IAL60VPrYPPVQDcXsklpmlXBgVPVVCnSLsFJrF15mA8ZPHFU4Ln42fnCUQ7MIX+/nb/Tz94d4V4ADaCmB6BjgSxHCSY5RJAjTQUJImBVh35zJrNEUhGyIWJ0QuXMrs5mndUUMhjY0XQF+3J949u535J8eUg3EiwyWZUMx2ivZORUsruC9oYTICSFcWs7OqVB3LncwJLtGhwwJn6NRfODlFz0d39KpiNr8AyM8GMfBUX54FACaXUxIu4gRHtVco7lbhD6yHQlzQkJImJWCJ8sYU2OEyNmz/p5SLMKuAIztXURS6BUNrNiqOlQFDUVsUSm7uIYBSTUPoo7wt2stt5+jvquIEequ0a4Af2sgsb83wgH4Yrj5ZdmYdxNTIMzEFhcTzl5fFIPhRP2+P5ZrjFCES2cCMkd3gHvCE+9JECAhJMxLwZNlSg1Vg3VFWXbU92EA+jSB6Rrjx/1c1lRmIMwBfHul5Y31Vpvhd9lQjEWlagzy9ycTSjMWB5BUwl+uJcsIY3HXINfdm9DcoQA4MtQaFltR5VDTf/wxPHJA0Zvg5G4RvnhqpqTVPH9SecfD/+GtqQ7jMAn9ObeuJbJDQkiYldAZHMOUEaNFqOd2ZqE+2SI8GcBQBFc9H983wqH1C611YkWyf7WhmC0qY+tapEoHXu3hT51Qt4vGOiWGcYxCmEej/Fs7ZVEyny6E+vzhcnvS+ltcTP9Evhh+sl/pDU06WebNAT48M6rvf3JAOTCKJw4rL3fPFG2eDow3RsTpQEJImJWQnLfBvFNDxAgBuG05hSpTY4RjYgA9v/iZ+G+PK4NhDkMPGp3GYixwsxon7lplAfBsl/p5RWwyySJ0AMDL3fwn+5X9I3zngCqEwgQUrlFdZRuLmRiaKGguATQp3TvMTwUS4wlzFEKZ4/Ao75wBHdr8Mbx0iv+5T+HA/9ohh+KQz/iihFEeU/DcSeXpE9MlV8+dTJ1JQkwNEkLCrBQ8WUa3qOpyMAeRIUaoXsTGYrhjp+IJgwE1aYdqcTHR2u2qRgZglwcAZI6wDAuD0/Dxqx0A8N4QVzj2DvO3PWqB4B4vH4qoFuHScuawAEClE19dJn1usSTmJop5T+ITvd6XdHXNMUZ43M/Dsjo9qrD0h7jM8fxJkfLD/9itPHfGLae3B8VYEr7+Jfnml+Xo9Lx/h5frJ7zw593MkBASZmUGJMuoopVLgBBpFqGxLrAnyGMKyuzZ0m2Wl7MSK476MBxlAe0mwCib61ql5VoG6dMnFH8MnX4eV7C1l3vCXAhhlRNNxQxAhZ0tKWf3XCDtWGf9/GKp2QVoNq7ed1tIZo4WoZgk3OnPtk/GiRl55MVTvD+k3mrowjAQwotd/GBao7hQPFFYmXfeHFBTkwAofBLjt3JHdGbv8KpHFvWsxNQgISTMygxIllEf5GgRVjrgsGA4oiZ2dhk6xYRlAKjNehyrhPOqGQf+2GcNpPlFAdglPHKZRVTfv3SKA4gpODHGt/VyT1gtoq90sBaXuhgAjcXMacFPr7BsXGZBctQTwAK3qKbI5cPhwCgHkO4aDcbx16/Jol3OL48qD+5TPr9d3j7Jq7bMEyL61pDlq2+yxw4mGVkc+N+7lY//MX5ghO8eSjq4J4wOL0/Xib0jvDc0LeIRU/COh8NQgpnSJzYvdAd5XMH7Q+p/Nx9UwrM8N2gaISEkzErBk2UuqWN3rJBYzhahxHBxDePAa72KN5LULFSQHiBMQUy9//t3i773ngygxJYqnKvr2I8usQCJqNihUWzvVXSLsNKhekH1MnxBtRNIzoMFcFbZJIRw/zAHUl2jx/189bPxnx5ShDg9ckB56rjy88PK65OcALx/hOuWdMeI5bHD7LmTSUf4QYdy19uywrFzkN+7O0kjB8L8g2G+tSf1HTu8PPfxkJNiKILDvoRFiLQ+sROSywgRUYSqW4SHfXx2ZwZNKySEhFnRkmUKtgALw79faPnNRywLS3NN2FnTyABs6+VdgQzXrOwWIYDvnmf5zgpw4Cf7FRhMUiN/u0QyVl881al4I/CEoQvhhlb2xSXS1U0ZfvspFuFZZQDgn6iLd1eAH/PzD4Y50lyj33hTed+rtnOLyHhviO/o5wBEZlDu7BzkPs268sfVoR/GHU5p5/PBfYpeASnixzsH+FgM23pTxeX9IT4dHksAnjA/6uMwhFcnaxG+7514YaIIVRfC3iC2dFIS6RQp3FWEIE4PzSIsWNao4BPzpdyvplc2SP8K5fkufnl9hhdNaBFKDHet4H/qVt70WAC0ujJ8druEBW6mh8R+dUwB4ImoWaOVDnZpHbt5fubj63mwDOCaRTiWNZCmcCz7bVy3bo2u0ZNjXM9SeWuQxxQlrl2oB3IugIspsEnYOcAXurV02TiQ1sVGN1tPGhbwoVr2ag/f5RF5K9g3zPUYKoAOL1+U8x3MpBgMwx9Dfwi+hGt0cor7p26+qmqCtYlTfczPx2IosaEvyJ87yWVuKfQPwpSQRUiYlYL3GtXJ/cpzSS1bWMoOjfI7385w817lSN+Wgfkl6msXlWbeYWk5W1LO/vcFFmhVFsd8fCAMpHlEUxAWIQPuWCFJDG0u5rQgImfrnXbcz40+3rEY9H4uPzvMZa52yfnFEcVY3p67RSicqDsHud5D1R9jGF8IdRhwWR0DoKtvSpiwwzs5i5ADOSZ/esIcwBEf79X0vnuSPtiXTk38TsI1qnDsGeZDYUQVDGyWwRgAACAASURBVIbx575psXFnPSSEhClROCIyJKZmNpoFhwX/cZEErVW3oFyzw0Sz0AlJCGFZ5v2/tFR6a4P1rxYknn30oCLm61ZkfYtSO3Pb8Nu1lk0XWs6tYKLNKYB3PZkvr1xzzRkPqocJhafu386Xiqypyf05WoShOLb2KKE49gzzUU3qxjTX6Dse/o5HH5qYusK1Tezy+qTrmzFD51SAeyPonUxnlhN+fnAkJ5kRtwJHfFxvHtszSR/sodEk0zYjen++B/YoenO7XLyj7w2RWKZCQkiYEmEOOi2TsMZmCMaG2lYJVgm6g05XxOzML1FNq/E8ex+bx0ptWFjKRHKNTouLZbcIXVbs3GC9qU0CcGUjayhmF2k9TnuCPD0pcY+Xd3gB4JMLErX5wmXX6efvDfFSG66ZJ52X5uXL0SJ8b4j3BvHeEI8pCTejsAV9Mdz0snzB03E1PzPNIvzackvKVKxTBv+kWPakLMIOLz+QVoORESGEr/bwQ9r+3kk23BmJcmNf2Yx0auU3vz6mPHlU1b9cutw9eSSXXJy5BQkhYUoKnjI6ZRpLmN4EYHUtayxWB+rC0P8sOwtc6nVsPNeozsblST/wjpus2duCSwxLtNnCH5snVThwoSaEbw7wf3svSQnHYvjLgFrHtr6VXd+iJumcCgCaS3Ntk+Sw4DurLGJAsc5gOKcC8J2D3BNWW43rUzWEJzYUVwNvIvYpNhr95B9uYOV21WEgSkqMcznEsieVNdrhxf6RnPYUrtHfHk8YZ5Oa6ShzjMWQXQhjCroD3MJwRQODpusAjvh4ejZyCjv6eca+Bzy3bNXxMPXsERJCwpSEC91fbcowQGSZSgwfnSc1l8Bl0y3CnD5OW4kCwGFRCyGysL5F0vcptqIsN4tTsLaRMahCuHuIH/HhPzqUDwzZjP0h3uHlolHquRVs4zLp/GoGYCDEofl+z6kAgGvmqdMzdOIKbvvzxFVvOwe5J6LaRj7N+SmyRqGViAhjSyhNe6XaN6fIqqbUVjsZgBtaJSSPLBZC2DeZOsL3vfxAbq7RwTCg1bkKfBNl3hrxRcExgRC+3M1ljnkl7FMLJBhShETIMMsLxQ4dmbJSt/acVhrtXyZZEjOjICEkTEmo0NX0p4Nwaba62LJytLiYXhefo2u0zMarnZjvZtJEummV8L+Wqr/xhuLJ3TRYJQBoLGYAPGF+xMdjCr74uqwbDf0hvO/lPQEOoNnFPtLEbp4vQZMB0WhGty9XVDEALlvClf3YQWVCO2lXRoswnvRBNCHkAFZWsr85W4KhxUG1Awz455USkkcWvz/EAXjCiOVccdCRsxB60hy/k7IIR6IcwDsenqVF6sY3ZABfXSYZM43F9yGjyOmILFPdgjTy432K5zR6pu/oN3HxBgkhYUoK3l/tdBAuzSXlaHOz5hK4tOq9HIUQwKJSNqFfVPD350jCSGoonuwyAaDUDguDLwqRJ/LmAP+jFoXqD/HdQ9wXg8OCMjsYsLgU0BJhhGYs1YRQJN1cUM2WlDPdPWscjpHOUARHffzEGD8mitMTQpi0m3CNCqVZVc1uO1tqLGb1WouDmiJc3cTOr2ZWCYG4KnthGYd8HIDCEx3vshOM46iPHxzlWZyHegw1fQ5iTJlEOzch+WMx7BvHtospOObnVglfP0cy1p6KiovsQiieTd+nK8CfPamczgTHfSOYIbNHpgAJIWFKZoFFuKSMzXezZqNFmFvWKFQhzGnnEiuqHAxAQ2594FJgQLkd3JBqqHve+kNqZK6uiIlDi4vyYJjHFBzxcYmplYgAapxoKMaScramkZ1TycQf7mgmIdTlZNcg50AorkYTU8ondAbD6PDyiAybhGvnsZVVrL3SYBE62cblEtOaDwiNedfD9ZqKHMOEe4a5zBGMZ0vmfEabMpFRTnI3CnUX7q5MybphGYE4A+Cywiol1Z6KtuyHs2b0jCeED+1T4koGWzYjYmLlIweUvxgmP3cHZsTskalBQkiYEpNbhAzA0nJWbse5FWyqFmGuwiaulVOzCKFVXOjXcd09qI/R0DvMiTcaCKHTz6MKWlzMOBukvZItKWNXNbCHL7N8dakE4LAvw9u9eEoRLsGUIJkon5A5QnLSB3+hiz9zggNw29Spiu2VTF/ShdXsumYJQJldHdaI5ILC3iD3Riae06Qrx4HRcff53XEuEkY8kaTDieXmOMQDBsnfOZDhJX8ZQFAGtAZ7RotwZRUDINybGW8yoKXVHPVxYwP0iIxHDyoYR8LT2TPMv/S6/KXX5UcPJNyh3YFERUd2UozjyAxokUpCSJgSk1uEgBY/u7iWieG6LhtsOf8cF5WOWzuRjugjWj/JGKGOnssqXq+XEOhORV11xEV5IKwW7zcmS297JVtawa5rli6oZqICMuPF+s99XExq3DWYFHMSxpw/lppu2h3gIjqlDwMxCuGXlqp1HeImQxhbW3sVaLLdG8J7Q/z9iUrrdCHcP36YsD/EXzipjMVSL/Ti/OeeVKlbhBnfa3t/wiIEUOFQvzbFVoiMJCFmL43Td1R8EIVjj/aJFI7vvCOLyG6OZS26T3sgnDhIX4gfzzp7ROebO2VjIemXd8gyTypuOfOQEBKmRG00as52UvNKmNOixs+cFjVGmGPKqCD3GCGAKxqktU1sZeVUhVCzU0XCy/4MFqF65FIbnBaMxVT/ocjY1GmvZEvK1HsXoeIZY4S7BtUxiikWoU8VwgwvETOPdMN6RRXTl6S7nUXG7GgUUQVv9HMGiNSe3iDfP8K3TjQNI2ERji+EI1E8fYLr3sXaItUz3FQyid7lMKQFGes9BHEFr/cjKDMAJVobICG051Wzei2zSeHI2IDbH8Mx7Zzr7UyfPqH8R0c2p246R7SbIV3P+kOIKRlmj6TzRj9/aL+ie333DvPHDym/OKpc90e5gOWNJISEKTG1RSgxXFjDqrXojrhY5+4XBXBWGcvYaDQj31kl/ela67XNU7YI1RdeVscqHBiOqBKYbhECqHEmxLI6uSxydR1r1tYsVPxImms0ruAdDz/i473B1O4zwmHoM2hDuybtQjD0juFnl7HmktQji/uMkSh/a4AH4zinUu0V8Otj/P0hvrVngozHXCzCkShe6FJENzULwxX10opKJjHVKe1LU7XxSAhhmnb2BPmh0SSLEECtkwG4qIaVWFFkRTCO/SOJWn4jx/yJifb6XYgI9IqMqixCaCw+0f9wg9r+onVO9mmUAKJa4rFu7j+4TwHwN6/JH3jVvu0FgYSQMCWmjhECWNeS+OmpFmFu1fSCCoda3nAG0F2ji8vY4tKES7MvzSIEUFsEAPtGAKAmeTrVAjfT95tXwqwSeoM8pXvn3hEejOOIL0Ohmy8GnuxgFJUSOroQ2iSkNFeDwSIUxt+aBvY/FkpLy9mBEf7Tw8qf+xLpMym87eFPdSrDEbhtYMA7ntQ16wxHuD+GXx1VACyrYGsaWXslc1nV8OQUYoSjadWHJ8fQF8JAJGERQjvPF1YzANUOBuDPfTxjUo8xq1MXNlHvOK9EtSbHW9Xf/yURyku4RrWbIdHcIGOdvpFfHlX2DnMgURLz/x9RADWbN31U1hmDhJAwJaa2CAFc35IQj4tq2MOXWf55xQxVdd01uqiUtbkZtCZqGS1CUbcgpg+muEaNWBjqi1j66HaRHnLElyGtUeEYiyXe1CrhM4sk492A2zBMMb1Hjy6E23sVAFc2MIcFn18sAYgr8MXw7jhhwvc8/PPbZQAX1bDlFSwYz5zDEldU0+qJIwqA86rYtfNYe6XavhWZYoTjRcV0iTIGRN8c4Ff/IS7KB/eOWgC4tG4SqkVYy6C5Sf/cx30xjEZTqyRFhaJ4mW53intKYUNntAijCvpCeOKw8pb2wYUQMiAQV3VU2METukZPam3hdg1yAK/1KcY+OOnDk88YJISEKTG7RagX2AFocbEvLpGum6rrcrrRXaOLSjHfDQDH/QjEE4OFjRbhWYYsmOqsNm5TCZA2qE8YCrpFqJ8lrfhBHaXrsODD9azKkRBppA1TTEFLllEDkOdVMyT3qBvvKvy+V21a1l7JxDjJjAHF0Zhe5qEev83NVlQxt01dWHqMcHsvD2YqLtQlKq4gYEjWfbmbi8DeB6MSki1C0WABmhCK/nYnxxJFn4KRCKDduIxoqa1iDSLhdjCTEP7TTnlbj8I1N2ZfCN4IXDYIR7cwCoVFaJw9khE9GccTxjE/F97vC6qZ6H/0Wp8i85x67+UdEkLClJjdIjQRwrqyMLS5WZuLAej0c2MdutEiXGIQ+CwWIbSeNd1B/o6hWk4I4TE/F/MUP1zPABRb1XyT0Sj6ghzAQjdubJMAVDqYbvyltDNNQSufSMxlRPL4Dj1MmOL51G3TFVXs0jqGcWZxjCSXTIgAZIuLNZeg1JYQDCM7B9XMmpTsWaNHVH/cb4iY7h21IjlGeGGN6nYWMVrhF+0K4JkTirFMQtiawqzXA5ChOAfQ4gKAoQh/c4Cn1DM8fED57rsKgP8+rgyE8PxJBcCH61mtE9C0s0crx8xuFIqdhR2/a5CL/17VyLZfb2l1seEIjvr4W5kM7umGhJAwJSJr1GnOrFFzIayuFhezS9Bdo+K6LHp6GS3CpUlCmO2wTcUA0B3AIwcU0SZtNAoRQArF1ZyUs8uZOL7u2BSDk84uww2tDEClA18/xwLgrxZMYFKL4OXJMYzFYJdUK22hIWz5er9q+b01wHWbjAN6e9X2Siba1HkjfGsvF6qgY8wlsUtqSR8DLq2XxHv9oEN5O1lBRUtxAPfsTtJe46F8MXWeorEtal84OUboVAOESD7nXQF+YJQbZXs4wqEJ4YgWLxQWYYWduWyIyLh/j/IbQ7vw0SiCcYjUm4iMXxxVnj6hALixTRKxSVEqo08ezh4mFHcDq2sZgF2DXPy3tojNdzOhxG97+GuFGKlIQkiYEvNOnzAdwjUqvIi6a1RYhOdXM7uUFJAzWoQ12YWwhEF09jqh7B/h61+Sm34R0wvbhSZ9uJ79+4WWu1dJQgh9MTWmuH6eIpI75pWwry2X6orww0sm+CoIW1YEAvUFu2yo14odx2JqeO9UgL/WpyrByTE+EkWlA3VFWFauDrHyRvD3f5E/+YpsFDahXmKH9krm1JZzWR0rtQMAB368LyEwUQW7h7gofn/qeFI3cyHDuvbfv0dBskWoLl6PERbholr1Sm60wrvG+JHRpHFOYpGtLnFkdbuIMhRbVWvy6U7lwb2JdaYEcQ+N8ld7OID1LZKITb7WqwxHEpOHsyeOiluH61skADs1i1Br+MAA/PKo0qe945nUQxJCwpTMnPH0sx4hG3qjcAZ0jXExa2llFVtewYyGWI0zYZRkd42KGOFTnbwvhO++Iz93UgnEcVYZu8IwQ3FZObtjhfSFsyTh2PzZYeXEmFiMusOnF7IyO3631lpXlPYGyQgJF61PKg2t7MTnEiMYHzuoAOgNJuKFwi96XjX73Vqrw6Lq3FCED4Y4AH+M65E/4Rq9skFqLGYib0VwcS3Tg5e/Oppo5tnh5WEZnjCPyPDF1PCbeqio8FUyAB94+ffek6NKhp6oukVYW8R0i9AY9dw3gv6QmpaiHRnQ7glSskaLrOofLqpg5yDfNcjFsHs9N1hI+xv9PBhHYzGrLVKN7Ps6lO19yngWYUoppIgRXtfMRP6tUFnRh0EUmfzxlOpsODzK0ytVxsvXPX1ICAlTQjHCM4ZwjQrBcFrQ7GJRBU8dVwAsK2cfnZeqdl9bbmGAXUJp1srIpmI13AjglR4O4CeXWT642XpVo3pRclgSf99PL2TFVvzqqGo56e3ihG0hQnfZaShWJzRBs9sE4nPdsUKyS9g9xANx9IW4nscvZmicU6FGB4WCeiNqGzNfFC90qddmoSsVDmxoZcbZyyVWVDjYfRdZPtrEwjJe6lYARBVVnwZDasXCk0cUcYSeIB8KwyqpdtsPPlDGYvjAyzNZhOqD9kpWq90HXNss6T3Nt/cqPLkvgXiLxmLYJYRltbeZ+CkVW3Gl4Rbkx/uUu96WuwKqVjkt+HiLBO3OYGk5AJxXxYRaH/ElZDUl3mnMxVU4PGEw4KwytqCUBeNqJwTVIixiAGKa5P9/+5WhsLo8PYXn9yenSwlJCAlTYvasURNhdI0CWFIGANv7OICLatjVTanXkO+skn671tLmZtnVaWk5S4nw3twm2SXogyOM0xNvaJU+peXCMKBWMzpzb0onMdXGQrJFuLCUWRjaK9m5lUzmwkzB+14uwpbiotyodadzWlBsRURW43b+GJ4+oV7oh4UQ2nFjm3RR8vDFq5vYN9slEe8cCmMshk27FaFPnojqHQ3E8fghBcDD+xWZ48ZW1fEognM7B3m6Rah30ikx3A6W2yFSW6ENqDru57qQCLO13KHGXIV6BeMcQJGFfVlrRwfgv48rr/XxnQNqmu5fnyXdvlyC5q4UDvBPLZT+51kStGoZkQIjknJf7eFiYuLbBhkWPV1FT7izytRPDS18q9/c9IfUsyEimgdG+flb4uK+4b+PZ86zPX1ICAlTIlLdzDiY13SISUx6guXSCgZA4bBJWKUlUqZwU5v02vUTWOsNxVjfmrj+6D5VPbJYltxzTu8jU+NQ7FO6bs13qQ+MFuHiUrS5mdOi5nnuHOS9Qa5wiDChMYilvTaxKl+M7x5Sp73rGnNlAzu7LGnlzPCm3ggGw3zT+/IfT6lNrnVn6X3vy8MR/OY4B/DlZZLxPmBHvyrMRkrGOcE3tKaeHd0oFMpXbk9qQa7HCOe72XXNUmMxa3ExoTfibABodbGLaxPzL/VIsIiwik6t51QwC8OJMR5VcPUf4lf/IT4cgTElWDQmFY5QY6dcEZvU52X2h/gTh5XRqHpvcdTHuwL88t/Hf3+S9wW5OG95h4SQMCWULHPGYEClAwvc6nVqiXaVP7eSFVnhHOdPMGHQDsBXlyWuP/q1VR+nkNJzTnQ6BXBl7RSnFbS5dYswsXFhKRM2rhBCYRFC63KipzXq+xtf64+hO6Bmo+gaY5OQcWCyUNB/f19e+VQ8IqvjnzzhRDOXvhD+7i/y4VFuk7C6lhkdyy92KXofTt1oc9ky3wVuaJUklrSGnQPKOx7+crdqVpbb1XM7qlqEgBZl2Lhc2tDK1mg+0p2DXKTpNhTDJiUMcT03uKEIAESty7wSNLuYzHFgRJ3aOBJNqo0ZDAHaXYUuhKI/LZDwBAxF8KO9CrQ+OMLEjMh45oTii2FL57Tk0JAQEqZETK0b7ypM5JdzDWmQumJ9euHpXj2uamTLNAepfliDRZi6BvHg+qYpesdWVelCmJQsIy7roiC9J8B7Qxxa1XwmizDx+FQA/hhShHA8xAvDclKLmcGw6rdc28QcFjx5RJE5zi5jDotqtFU54LRgyGAOXtWgPhjPImwoxuX17KNNic/4h1P8q2/IV/8hLgKNFcmuUT1GCODqJvZ3yyXdufqOh4ssGGGuiUTfUht036/YLvrX1DjVRvB6BuzJMRz3c/Es1xJQ65MtQn0oim4RKhz7Rji0Yg+9ndsRH/dF8fsuJTYNNiEJIWFKwpQ1ega50BD0ErKxsoqJoNHpwIB/XildM4/BYGTo5leKa7TKgXklrNiKq+qmaBHqySBGMSuz4/J6Bi05tjekGiJ7vNwbUdP9aw3WrVFExSQKVQg11+h4716R6SndIry0jn2rXT2fwgksmun8zdmSrt+ia8+GFritHIYYYTrfbLfcooVUnRbsGuS6OFkluGzqOoVrNGgItzNgaTkTJ8oqwR/DW4NCvQBNCD9/lqTnwRqHXNYWqfKmtyA46uMyV6v7X+7mfaHES/SQ8ycXqJ+uyokUj7dwjRqEEL4YH47gDU/+f/YkhIQpCZFFeAa5oDpxEa8rQrUTj1xmyUvX71sXSeK6r3tcy+zqBTHdumqvxEebpCLLFJ1jS8qZuApXJvd+W9skQRNCfUoRB7oDXKT71xjqQIyvFZksYt69KKSrGr+rXMqbCtelJ6zW1Fc72deWW8T3WTiBy+ywMHx5qaTr99VNDECbCx+qimN8ixDAdc1sXYtkYbBJ+OISCZr/U6yQaZnA23vV9gXQLEJBq4td1chubJXEs0yTwKZiMOArSxN/+LqiROCwxsnEn3JHvyaEfg6gcwwAfrJfUS3CYgagzc0sDFYJX1qSyIFKGZlpdI0C6A5wseX7+22ngsgvJISEKaEY4ZnkwuQ0yB+vtqRsOR2ELSgy8gEwoKZIVYIU2ivZjW2n9b5XNkgAKpMLHIUGVNgZkDSqvnOMh2W4bEkiYdQz4ZYIxvF6P//Ayx0WrBh/6GOKEP7t2RIDhjSjs9qBaif+aoEEzSIss7OPt0htbrZGqycR3s5mFy6tkTF+jBDiHDpxWT1rc7P25CUJD+TG5ZJNwkP7lUOjPJgp7+yxyy0f1/rCr6pWOwk0FrO1TczYM8EmJbS/tkj9nuiVi0d9ANDp58f8/KVuRYRFRVjRLqHFxW5olURjBG1tSTNVhiM8GEdPgNslLHAzrvlg3/RYtvXlOUuOhJAwJWHRYo2yRs8IKbMPs3f1nCxLylmxNVHbAK06oixtUvH51ez6ltN6a2FdVWay21y2VNfcwVEg2RyE5hpNcUX8oEOWOS6uYVl89ZXJrtHzq1m5A3FFdf2J5gMblwshBIBSGzYukwCsrmN2CRLDFQ2SVcK8YlxWPYFFKLixVVpUmjqIQ9jEKyrZzfMlmWN7L1ezRpM/UZubXd+iDve4UUtDbSpRl2SkSVMyYRE6LYmyd1FTeNzPn+7kYzH1k+qxwEWlqUerL2IfaUycpeEojvs5B9rc7OzyxG5X1MqfWZDnlJmpf6uGh4fXrVtXWVm5fv364eHhLNu3b9++cuVKt9u9cuXK1157LQ+rJuY8YaojnC0sKWdnlzFjlqNITilLi4Gtb5UyaljuZBFCpI2EFCHA2uT010oHzqtmz1+TpEIvdHEgUcCXkQo7jE9XO1XxEymXom/nBdXs+hYmyhaXlLO1TQxAiRUX1bL5blZmx1mlzG3DOWVylWPiAPkNbWxxKUuJTTZoIdhLahmANwd4OA4GONOOVuVQm57f0JZ4ycfTbkSu0j51rTPRZFUghLBzDFtOKNBih3pY8cY2ydhFCMCKSujtFAAMRyAaGDWXqKlMAJwW/J9V4bzf/05dCL///e+3trb29va2tLTcd999Wbbfeuutd955p9fr/fa3v33rrbfmYdXEnEe4RilGOAsoteEjyRLy8Rbpi0ukC9K8r1MrHzRyVhlrKmEVabamQN8uVPmNfo7kluIArm5iO9ZZz032Nwrz5JLabOuzJrfaqXYyMaZKRO/0kVV609Qyg3Be2aB6OMU5sTBc1zzxuWh1sVvmSxXJHmZdh0Tm52t9nANFVmQ8Ize0SovL2DkV6pNL0nogwGAvijsGvcwG2gjGnYNcnElRPq9bhLednfoRPrVQElk54lDDEd4d5ACaSph+fs4qYwtc+a+gmPo3a8uWLRs3bnQ4HBs3bnzqqaeybC8tLR0dHR0bG/P7/S6Xa/xDEkROKEBUgcTgICGcFdzQlnQh+tpy6eHLLLk0TpsCH2lk6dFHge5FXFjKAHUChhhzoSOq70szZWy2jx8gFBi9ozXORC/WYmvisVFIdNY0SMJfqmctGQc7Z+GyepbiGtUTUlZWMbukuivH86zc0MZumigou7qOiZpR4UNO76t3eJTrRZDGmGJ6stXSciZuO0Tr1OGoOq6yqRhV2vnJ3rdvykw9D7W7u7u1tRWAsP+ybP/Zz3524YUX3nbbbQB27dqV8WivvPJKVVVV+vYHH3zwqquuMm7xer3hcNhmyzqFc/YSjUaDwWAsljblc85QGot5hkYAOCTe399f6OWcaWKx2NjYWDw+PZ2mCsQCjtz/kiMjI4FAwG6f4hVxTbltcCDzz6cYLsAOYFVp5PCoHUCFTbnKNdzfn8EEsbGKGGcAiiw8JLMquyL5+/uzzl64traow2F93WMDwP2easlZZrPXOfhD54+NeLLVhCwE67PZ+vujCy3W/v6Q3+9f5VYyriqdWJwBFfp/iyIj/f3qxz/LXbZn1ALAwZSMPyU78Nk6qb9/gsK9j9aWPNNjH/H0A7BEi4BxmynU2JXBgWx/aVvIApS12YNFFmdIZnsGQoC9VBmzhhXABcCpRAYHB51Op8ORq5e8qqrKap1A6aYuhJxzxph4IMtylu133HHHt771rdtvv/3+++//p3/6p5dffjn9aJdffvnDDz+cvr2srCzlG88YKysrm7NCGIlEAoFAZWVloRdSMGSr1ekqhxfFVlZTU1Po5ZxpotGow+HIeNc4R7BYLCUlJblfB1O4pWLc6FqdS8EAB7CmxfmbLgXAJXVSc311xp3LHLInLHJY2IunsKJKmvDb+KMrsHeYt29RJIbFjVWPNuFR9ZmJf8631MHCcEUFrNxlt9urqzOvKp1qwCrJcU3Lzq4vq9F8zqtqlD2jHIDLYRlv8bn8wD69BH8ZVsQR6ss5MK5wzi8b940ES1wA5Pb6kso+3h3gR0IOgJ9d53JaGCADqHY5qqqqioqKnM6sU74MSNLEjs+pC2FjY2NXV9fixYu7u7ubmpqybH/rrbeefPLJ+vr6O+64o62tLePR7HZ7XV1dLu8raUx55aZmjn98AApjwhpyWtkcPA/0BTjNM1CSpfmLkwO8tgiLyyRxNV9WMe4bNRUrMue1TnZBDXvxlLKiKqclNblwy3wUW2GzTG79krb4WGzSH7/cLnvCaHGxcyswv1SStMSkFVV44ogMoMSak1qMx0easLhMEUcos2czH+e7J/jN1hZBYvLicqnCLncHVO/0vJLES8rsbDp+AlM/1rp16zZv3sw537x584YNGwBs27Yt4/b29vbHHntsbGzs5z//+YoVK/K0cmLuEqGUUWIaEOG0pdoAXiTPGU5h901W72dtB26xcQZsHAAAIABJREFUitDghAFCQaUDv/mI5adXnNHvrkgcXdvIfn+NtclQt6ev+brm0wrH2iV8YbEqJRljeHqPujb3BIeySqh2Qq/6EDlxTSWo0o4wTTHCqQvh3Xff3dHR0dzcvHfv3rvuugvAmjVrMm7fvHnzCy+80NDQ8Nvf/vbRRx+d4LgEMRERaitDTAMia3RJWU5CqDMpISwIInE0XUJE+b9Vwv9aerrWlZ7ulDGN6BPz1WfbXBOfpVYXaypmerKuVUJdEavS8oxKx28jcDpM3TVaXl7+/PPPG7dwzjNuX7JkyY4dO6b8RgSRQkQBqNEokW+EFbKknOnpnbkI4aJSVmrDsooZLIQOIJNE1RahvgiX1ie1d5kaemVLaVppCgP+aoH00H4FhgEgWVhdxySGNq28oKWEWRhK7bBJiCkzzyIkiEIhWleQRUjkF9016rbBLqHGma13qI6F4cY2aSZ/G4VrNF2iALRXsvRmMaeDLrfG4YV6lf38iVyjAETZjF5BL9pzM60NQkaL8/Shm2rCfERlkbZe6HUQswvVNVoOABUOpMzXzcKti2a0RfHlpdKaBnZZfYaP8z8WJZp65wXdYltUyg6N8goH1jSyMjucFsQUNOdgeq6uZTDYjvpE6CoH6w/xGVdHSBCFIqpQo1Ei/1Q4UKJ1Pa10sFz8ooLszdUKzofr2YczqSCAzy3Os4TrMbzzqtmhUb68gl3TxADUFTEO2HJ4N5HOo6fVLNQU8ep5rLEkte1tviAhJMwHZY0S00GFA2eXM3GhrXTkFCAUpDceMwt5X7husa2qYr86ihonu7ZZAlBXNLmgfpuLMYAbJhc+8CH1B2/obJ03ZrRFTxAZoRghMR1U2Jk+FrHSkXhM5E6xVe2dtryCOSyocqhWYF0RyyVl1Hgc0XF00Rn5K5AQEuZDWITF5M4g8orLpg7FBVDpSIxIJHKHAWJ+vduGxaWsWqv/qyuauIgwhTY3kxjmT48vNAUSQsJ8RJUMM+EI4vT5UK162W0qSRqRSOSOCBMWWbCknOnNsuuKEomgOdLmYvNKso14zCN0U02Yjwi5RonpQR/v8KHapBGJRO58Z5XkCaPVxZaWJ0Y/TtY1CmC+GwOh/C8vIySEhPmI0nh6YnrQ/e2r68hbNkX+Rhs0uKSclWjnc2qu0aFIXlc2PiSEhPmIyICFskaJaSSXUnoiO0vKWEibS9TsYpPtX9PmYqPR/K8qIySEhPmIcAZKliGImc2SctYVUIcmrqjMMN0+O21ujJ2puat0LSHMRzAGWNXkNIIgZibFViwqZfrjydLqYtFs44rzCfnBCfMRiAPTNpCFIIh8cTqtBpyWM9fKnISQMB/BOMe0DWQhCGKGcMZa9pAQEuYjGGcAuUYJgsgPJISE+QjKALlGCYLIEySEhPkIxDjIIiQIIk+QEBLmIygzUIyQIIg8QUJImAzOEZEhMZSQRUgQRD4gISRMhgxwwGXN/yg1giDmJiSEhMlQOACU2kkHCYLIDySEhMmQFYAyZQiCyB8khITJEE2XqHaCIIh8QUJImAyZUzU9QRD5hISQMBmyiBFS7QRBEHmChJAwGUIIySIkCCJfkBASJkMhISQIIq+QEBImQwghVdMTBJEvSAgJk6EAAIrO2IAWgiBmOySEhMkQFmHR5AdeEwRBZISEkDAZQgiLSQgJgsgTJISEyRB1hCSEBEHkCxJCwmSIGGGxpcDLIAhi1kBCSJgMco0SBJFfSAgJk6EJIWWNEgSRH0gICZOhlk+QRUgQRJ4gISRMBrlGCYLILySEhMkgISQIIr+QEBImQ6HyCYIg8goJIWEy1PIJSpYhCCJPkBASJkN1jVIdIUEQeYKEkDATCofCwQAnuUYJgsgTJISEmQjJAOCwgByjBEHkCxJCwkwE4wDgIL8oQRD5g4SQMBOhOAfgoK8tQRD5g64ohJkQFqGTLEKCIPIHCSFhJsg1ShBE3iEhJMwEWYQEQeQdEkLCTIisUbvEC70QgiBmDySEhJmIyABgo68tQRD5g64ohJmIKhyATaIyQoIg8gYJIWEmYgoA2EgHCYLIHySEhJmIKgBgpa8tQRD5g64ohJmIUoyQIIh8Q1cUwkyQRUgQRN6hKwphJkSM0MqofIIgiLxBQkiYCdUipGQZgiDyBwkhYSbUGCF1liEIIn+QEBJmQtQRkkVIEEQemboQDg8Pr1u3rrKycv369cPDw1m2x+Pxr3zlKzU1NZdeeml3d3ceVk3MVdQYIRXUEwSRP6YuhN///vdbW1t7e3tbWlruu+++LNsfeOABn8934sSJ1atXf/e7383Dqom5inCNkkVIEEQesU75lVu2bHnmmWccDsfGjRs3bNiwadOm8bb/4he/ePzxx4uLi+++++5Dhw5lPFpXV9dDDz2Uvv2qq65qbm42bgmFQjabzWazTXnlpiYSiYRCoWAwWOiFFIZQzAIAcmzOnoFoNDqXvwAAQqEQY0yW5UIvpDDEYjH6AnDOFUXJcX+n0ylJE5h8UxfC7u7u1tZWAML+y7L9xIkTv/zlL9esWbNgwYLHH38849GGh4ffeOON9O3nnntufX29cUssFovFYlNettmJaRR6IYUhHGMALFDm7BmY418AaGdgwkvbbIW+ALFYzGKxWCy5psw5HI4J95m6EHLOGWPigfHuLH27z+fjnO/du/fHP/7xbbfd9uabb6Yfrb29/YknnsjlfWOxWFlZ2Vy2CK1Wa1lZWaEXUiCsMoBip33OnoFoNCpJ0pz9+ABkWS4pKcnl6jYrERI4l78AiqIUFRU5nc48HnPqd1WNjY1dXV0Auru7m5qasmyvqam5/fbbGxoaNm7cuGfPntNeMzF3oRghQRB5Z+pCuG7dus2bN3PON2/evGHDBgDbtm3LuP2aa6756U9/GolEHnnkkQsuuCBPKyfmIqKgnnqNEgSRR6Z+Rbn77rs7Ojqam5v37t171113AVizZk3G7Zs2bXr11Vfr6upeeeWVRx99NF9LJ+YgMRJCgiDyzdRjhOXl5c8//7xxC+c84/b6+vo//elPU34jgtBRC+qpjpAgiPxBt9aEmaAxTARB5B26ohBmgsYwEQSRd+iKQpgJVQhBY5gIgsgbJISEmYiRRUgQRL6hKwphJihGSBBE3qErCmEmKEZIEETeoSsKYSZoQj1BEHmHhJAwE1RQTxBE3qErCmEmojIHYKOCeoIg8secEMKeII/nOruKmNEI16iFUfkEQRB5Y04I4bsenArk+dJ5aJSuxQWAYoQEQeSdOSGEHV7eOZbnY/6hi4SwAFD5BEEQeWdOXFE6vPy4P8+69cducrYWACqoJwgi78yJK8oHXt6ZVyHsCvA93jwez2QUyhaOKeAAYyDPKEEQeWT2C2FYxiEfP+7P5zF3DvChyBxyjSocR308Iqshuv0jhfnswhwkFSQIIr/MfiHcN8zjCjrHUq/dUQWj0Skec+cgD8YRip/u2kyBN4LrX4r/371Kf4gf8XEAbw7wgnz2KAkhQRDTwOwXwg4vB9CZZhHuGuQn0tQxR3YOcgBGo1CZpfahN4ILn47/oYsf8fH+EA6McADHfBOfuuk4ISJTRqLaCYIg8spcEcLuII8mZ7ds7eG9wakcUOF4xyOEMLHxg2EekU9jlQa+skP+0uvy/XsmnYwzEMrPAowc8fFjfg7giA/9IewfAYDOMWTJwu0Ngj0acz4eE1o4GsXhPJWaiPH0FCIkCCK/zE4h9ITxjofvGebQhFDhOKkZMX0hANjaq/QGp3KB3j/C/TEAGAonNnb6+elEzvTXesJ4aL/yyAHlySO5CmFXgAMIy/jbP+dJig3op6jTz7uDXFiEnf5sWbieMAcQUzAQBoB/3iUfz1PtivDHMrIICYLIK7NTCB87qFzwdPzc38W/+LoshBAG7+im3XJExl8GeO+UTCjhF0Wya7Q7AP2NptDF5sVTfDgCaKqGZHMzO5t2KwBOBfhzJ5XfHMtzUYduNEcV7BpUhfC4H53+cS3gQe3+oDvA3+jnDx9QArH8SJcnAgA2MggJgsgrs1MIdTn5rwOKfl3WjZhfH1MeO6iE4uibkkWYEEKDRdgT5O9rQvhqL+/KoZGNbNhl56Cah9I1ph88p7X1BvHTw0pMwYkxAPj7v8gjU80AGuf4iWXs6OcHRnlYRl+Id45h86HMouvRVn7cz7/4uqxwBE47s0b87QZDHFRESBBEvpmdF5WTYwAwryTJdnjiiNLh5TLHYBjf2inDYO5koVuTtJh22f/Am8Fo6w4mLMLuAH94/8SW2W+PJ/bZOaAJYUDd4o8hmnaM7jR9feaEEorjjX71HftCuO/9fDpIhdHssADAwRE+FsOOfq5wvDfE/+/e8YRQffCv7yp7hzmQByH89/eVz2yVRWCSLEKCIPLL7BRCYZDd1KZeMltcDMCf+/gvjyqDYeg2Sm9oAqtrS6fyu051n/eG+Oe3y6E4hGIh2WjrDnBdCPtD+K+DSvbcGW8E/3VAFRJPGMf9/IgvsXLt+EkvefEUf+JI6oK3dCoA7t0tP3NCPdqxvFZMCovwvCoGrY7+D10KgMOj/MBI5rCoR7s/EDFaAGOxcY/fk5tR3h3gvziqfPMtGYBVohghQRD5ZJYK4RgH8In56qf7xjnSV5ZJALb28n6D+GW3CH0x/N1flIFQ4mr+88PKh56N92uRxSSLMICBkJqG0x/iAyH85ng2o3DnIH9jQA2z7RrkHKq+njTklRhjkME4vrJDHkhW7pEotvVyAP+vvXsPb6u88wT+fY8ky5JsS5bt2JbjS26OkyZAriwhAQIMoUydkNKd9GF3uDzdThkmO+1seZppy4ZuyTMk9DLw0NmHZVlvt9NOG9oJJTSlpIYkhEshXEpIcALkAomjxE5sx3ZsS7LOu3+8uhxZ8iWyHOVY388fPPaJJL/nWOjr33s7O1vjc2I7o886O5DOaOUQ6hItLI1XYX84GW/Ds8cT2qPacDapU3e4ivBcAE8fHlsQ9sVfnxUhEWXWJAzCC4M4F4DDimvLhTsPANbN0LYssdg0vHNWfnLeGIQjfQp/e1+49YKMDTH2huJzUJWEirAv/q8qKX9yUD+bWNJ9aKif3mqX/YP4U5tEdNAx2jUqAeRp6vXjz33onfCxnnhjlB2f6cndp53RMcJ3z8kjF7mxXF9SYqkgvKokHj4HO41BmPDjd3ym/7FVDjlrAMNNlnn6kH5ibEs5jYUjxwiJKLMm1YfKQBiIloNTXUITWFImagpEhQMFNiwpE4M6th1PKLOSN5dRgdTSJZ9s0WFYnNdrCIkpDgB4sVV6/zW05o/h3lDkdaJBGIm3aVtDvzbUhT/ar8dKtH3tOhCp595q1wEciUyWAYA5HgFDRXikWz52UDUmITZeO5MiRTqjder+jsgkz7H718Q1G7pE24AUwJXe1FXYO2cTpgW1dOFvXwt/lpRtKSvCsMSTh4b+rZBSIIx2w/xeVoRElFmTKgg3vhPuCkbmm9S4AGBJqVhaFvngXFAiALx6RgIotEWeoopCY8Hx0ikZCOPj81KtB2+Pln3GgS71moM6OgM41iNbo0/ffy5eEVY60RvCf94V3uOPFEkfdsnfRkfyVBW4y68D2NceeVZXEKf6pACu8AoYKsIXT0ZuLDykIkyZc7Gu0f3n5KGuUa6Y8cQ7A/hp4kTQ7hAGdRTlDZ12FCOBl1qNQSiPdEsVz49fY6lwoNwBABcGcaxH/uO+8E8+jL/+7z7Tj/fIs2PYstXfJ40P4hghEWXWpArCV0/Lb70VVlPt1QSZpVPEkmgQVrsEEFlW/7UGTVU5p/sxqOOqbYMvnYp8vDa3yl1+GRv/a4t3jQLRfrmry+LX7eyAbI1O9XzfUBH+ea3tbxq0oI4bdgz+/RthAN1BqCQ43iNVofmnNvlhV7zD8/UzMqSjNB8+F2AYg1SFo7ExyqHzkdpXfVtggwC6gpFZLfs7Rlnj3z6A3xrq49+d0IdsAXM+KAEU2URZPrRh6jDjWkxjMK+tE2/fbr23XgNwIYQDnXLL+/o3/xSesXXwxx/oQORSjKUiVAOExfbIt7wZIRFl1qT6UPmkWz59SFf9ezOKBIClZfGKsLog/sh5XjHbIwD4+2T7ANoHsOqFwSdbdAnsOqXv9uux2za1xybLDALALVUCwJIyERupOjeAWEV4qEv2D+LcACwCpfl4dKnlc8UC0WVw3SHs8csPOmRsJWIgDDUTUnnltA6gyiVK7KoilAAksNuvxxqjntl6QXYF4e+Dy4rbqsUcj3BZUVsgCm0IS3QHEdTR0iUPjbi32R6/roYttx3Xz/Tj2ePyXADdhsJX9fe682DVUGIHooOXRs99qv/jvvD5ICRw+LwEML1QTC8UpXZR5RIrKjQA3aFITRzUcbRHHu2RZwcipWT7aBN3EV00srxcs2oQvD09EWXa5AnC80G0D0ACr56WAGYWAYDPKZaVJ1SESrlDVDoB4LPeSOdnWOJ/tuhHu+W5AFq6sDPa43c+CDW3U834WFGh1RaIejemFUZeLajjw+j8kaCON9qkBFQJ5c7DH261AFCx2h2UAP7lQ/3dcxLAolIB4PcnJBCpt/b4JYAqJ0ryAWDvGXk+iCPdsn0APqcoskXumHE2oH1lb1iVX7M94qoS8aVpYnmFqClAsV0A6AzIQ10yqKfuO43Z5ZcHO+W394W/1Bx+/ED4xZN6rKmxcwegJhyVOwSAO2dqP7jaouJQXc832+SW9/UTF+TJC7InhCkOHFlnPbLO6rACgMsKAK+dSdjWtSOAkxciid4ZTNhYICVVEU4rxIxC4c4btjYlIkrP5AnC2PI+ZWZR5PMyVsTUGCrCcgeuLhMA/u2IHpsOc6hLvn1WAjAWbTK6HkBVhIU2/NV0UVMgGtzxz+P3DVNJ9/h1RGMDgM8p8jSc6ceFwcgr/O/D+v/9SAfw3au0a6ZEHqa+ePesBFDlEktKhdeON9vkjz6IbBF3ZQmmOASAswP47n7762ekevAcj7jCK9bWaTdUatUuofoPO4ORaTvngyMt1Nt9Su72y83v6xL4wQe6mjJqDEJVHaogrHBGrtsD87W9jVYAN1fFr8DZAajxSDXNJ8ZlA4DeEF44GR8d7AjEc1GX6BhtMzk1+6nKJRo8oszBGCSiDJv8QRjjcwpL9FiFQ9wxTatwYH+HjC1FD+nYelQCON4jjXcRUiNzaoywwIb/+jnNItDgiT/g/XNAtEJS43lqkggATaCmQEjgQIfUJTQBXUZmol7hFbGVjl+o0RBdJ1flEvO94snlFgD72iPr9K/0irJ8APi3Y+I3J6w9ITx5SAdws08sKhULSsTKSlHtEsV5ALDzpHzvXKy3NvXlOtOfcH/d2HRW4x2MI2OEeQLRaFf/XVImXFbc5DMGoVRlsfHvA0QrQgCvGya4dgQStjJIXnc4xJEeAJheiDkeTMkf+bFERBdtMgUhAHjyAKAsP1LHGNk0VDoFgNluMcWBPA23VmsA9pyOfxDv+ExHdAsVAPkWAGjrR1jildMSQIE1EngNhtJHVV1fqBEA3mxXQRj/12mFQHRlRaVTzPUIAIU2TC8Sa+sEALslobqqcgLRGvGtdqlS9gqvUBXhlg8ij/ygQ1o1fKFGU41cVCo+VxzpGv3u2+HY/mfDzZfZ7Y9n/RyPKLbDaQUQv4PxZ70y0jVqA4Cf32CR/8X2D/M0AAKY4xHXlidUhC+26gCunpIQhAXR2bnGPyw6Ajht6Ckddb6MWvo5o0g0eERZPitCIsqwSRKEA2Ec7pIA/tNMDanKQaXaBQE8udyiSkNVt6knqo/sIevT690CQPuA3PFZZPpMQXQV25DSx6rhi3UaEBlQjFWEAOoKBaJBWGSDCr95xUIA0wrFVSVirkdMLzQEoUsAmOoSlU50BNB8SkckCIHErbqvLRel0QrJquHz1ZrqGg3LeIU33DDhrlPx4/+0ROv4a9v/u14NZ0YOfu9dXW1KnvwnBYC/rBEVznibj3TLl09Ji0BjTcI7yplqZktHQBq3MmhPqgj//o3w4t8O/rFVAtAl1A0RZxaJBnfkIhARZdAkCcJnjuq7/BLA7bVaaf7wQVgg7q3XbqiM/KsqL1T43VKV4lLMLRYATvUhtgAuVuIMGQyb7RbGztISQ+FSV2AIwjysrdMAXBndq+X2Wu0Kr/Da472Ic6Kvs7RMA9AbQpENs90iuVfwRl9Cm/Mt8NqHPubwMBNHjYvxZxQKRANbVYTtA/jlEX3f2XjX6BCra7R8C4qiV+PXx2QgjOUV8WBWCqzJT8X5IE5eAKLle3JF+MYZ+c7ZyA2BW/vkQBjlDhTa0OBhEBJR5k2SIPzJh/rJCxJAdQGur9RmDBOEi0vFD662xL41fqreXieSF6hdXyEAPHNUb45OIo0FYbE9oey7witKDeFXYgikKhcQXfNXZMPCUlFbIK6I7tWyti7ytVr4OL1Q1BZE/umWaH/p8gph1SKzTGGY/rOycuhpqtUaAOrd4oH5GhLH/IxiO8KI6FITFdjHemRXEE8f1gfCeKtt2IpwQakAEJu68mlvZObOkIfZLSl2RNNlpMN2njcy/WcItSOBmkSjerzVXzbuvGH3uCEiSttkCMI32+S+6CTPSqdYWSnU2olk/zBPM9ZMxgGnaQXxpfdKoQ3LKwSAd87GdzYxljgNHjHLLVQsXeEV+ZZ4TJYYCqMKh0B0q7aiPCGA2+viQXiFV6yuFQCqXQBwgyHbbq8T6puVPg3A9ZWaWjmwvGwQgNM6dEAO0S1vAMwrFt9fZBHAZ70yeX2CcW85n0uo0cHSfFQ4cD6IGVtDjx0II7qiP1b2Gakfc+tUsciwH3dpqpksrlRFodoEdX6xCsKE9snojgS/PqoP6pEBwliJrxYmEhFl0GT4WIn1W7qsKLLhhkoxXNfokOrE2Nk4xTG0wCrJF3M8ojAxBgoMO102uMVfTRO7v2D1OSPBFisK1Yp4RS1YVFSofLFOm2+obFRrVUV4g2Eqps8pVNSphnntka3XVk4JO61YVi6Sl7fPdkf2GW/wwGFFhRNBPcUKCjVEV+mERcD4R8MNPg1ARyC+wyqGqQiVnyyz/K/l8Qq71J7isrtSDROqIUx1UwtV88V0BSKd1e+dk/98QFeV9OxorVnOrlEiyjTTB2FbP2IbW6tJoXOLxfyxdaCVGT5VpzjEDZUJV6PEDouIdEjGIqfAkItzPKLCKa6ZIt5Za/0PUwSA0mi5aawIKw2TSoryAGBFhUgus6oL4pkXc3ut5smL3/xB/evsIn22W6ysTPG70wQWlwpEeylVb+fxpN5RdbuougJRV5jwR0Pspxsb4U41RhhjrAJTV4TRMy204eaqyCIQAE4rVk0VAPae1o01q/E+Wd97N6w6pecYxl+JiDLL9EHYfCp+C1xVewlEVgKMKtY1Oq1QePKG1lhqnG9pmahw4JvzNQCaSHjlBo+odABAhSMySyWWrMaKsMQe3x5T5V/KYKl2YWaRGLK99do6cV2lFlv+qDpOZxWE53jESl/qfFIdvGp1h9r+5ljSzZiiFaGY40FyEM4oEjs/H6/jRqgIAUzJF7FHlqZa26C6Rme5xYIS8cfPW2+ZGrkWUxyi2iVmFImuIP58Lt7CM4ZitG8wMsloyBxdIqIMMn0QGvfuMtZeYxGb96gCwGnFbdWRDbxmu4Wa+bm0TDx2jUV14rmsCRnW4Bn6E2N9g8aKUBORYUIMMwNTqSlIkW31bvF3c+O/o+sqtUIbpjrl4jKxuDT1Sy0tEyKaHHWFAHC8d+hj1EWrdKLBnTCeOsstprrE387Rbq6KR3LKMcIY1QGrlKWqCF+6zXpknfXhRZqKSV/ig9WV/+ab4djdo84k7T5q0zDc7CciovGbBEEY/9w0jsaNkVqlHhuZ+2/zNY8dABprhKoIb6vW1k3X1NKCgsRIUAWN8YjqGyy0Dd2cOhYVI4RKtSthpkzMLYa19p48fKlWagJ/PVMb7iYMS8vEVJdQTVUTUD8dviJs8AwdT72tWqhbRqjVIF+s06oLRgmh2Cuk7Br12jG9UKyaqqnhPV/0Twc1QHv/XK3SiT3+yP4+YZkitmcWpZjQS0SUKWPrQ7yMGSvCiovfiLIsH590x8fGVlSI6yu0ExfkHdM0te+22jx6frE49B+tDkvCczUxdO6GKnpKknoIK51C7VczYhCKghFrL+Vv6iWGiRylyiVuimanqupaEyfLvHxK/vbTyGSZBvfQLN+02KK6eT9XLG70aRuuHH2P65lFYu9p1aphH+vJw19UCQBTXZEj6k+QBSXim/MtD7wZfqtd+lxy1QsJ9/C9dap48aRsSFqVQUSUQZMgCCWAqS5x8oJMryKc5U648eyzfxGJuz7DZ7LdgtljGKZS+VSStKo91h84QteowwrHGO4wtKAEFy6M8pi7ZmnGn9ua+PitR3W1jK/SIRaVCntiuse6N/9pscUxtneHqgid1lGGZm+dqsFQEcZ+kBrUfKtd3ugbWrn+ZY0W0vUG95iaQUSUHtN3Oakbw66oELj4MUIAZfmpOySByCLCi1JTIKYXxu+AGPM/FlrUTxl5vC1Tro+2vCpVRahurPHfF2jJKWg0xhRE9I5XI5SDxhesSqwIASwqFRaB/R3SuGZDqSsQfzdXm1PMipCIJtAkqQhXVIhfHkHlxS8ym+LAvGE+Z5NX6Y3q1qniyLoUl3SKA3/8vHXDvvDIMzAzJdabWZoPuwWdAfQNRsq1vkEc6JQ2Dd++cqwF36hURThCb61RpVOoW3DEKkKXFQ0ecbAzfusrABuu1AAsKEGFU2u9MLRSJCLKIHMHodohJd8S6V5LoyKckj90+eAEsWr40dWWS/yJLgCfUxzrka0X5Cy3APDuWTmoY2GpyFQKIrpDW8opo8nyNJTl40x/vCIE4HPiYGf8RlpFNmxeEq9Va0abrUNENB7mDkJVDlY4xcwikaclLFoYoxUVIo2RxbRd+k9cqi9dAAAP/0lEQVT0KheO9aC1D7PcAPDuOQlguKUX6XHnoSx/9K7RmA1XWs70y7mGNfKl+QKI7LINcLEEEV1S5g7Cw+cBoK4AnjzM84o0Pj4XZjQSLkNVTgHI1gtSpfDHaseyTK9Pn1k09L4TI1A3NTRSz43dpHC4HfKIiCaCuSfLqL22Vb/o8nJ+eqagJqd85+3I3TlU9+Nwm5KnbWaRSLnR6BipajLWb5zx5hERjcDcQfhWu47oLRdWXPwkz1zwzfnakjLxWa9UW5NH7mqU8YrQfREVYbIh44sZbx4R0QhMHIQdAbzZJhENwuW8QU8qPqf4yTILgP9zWO8N4dNeqQlML8x41+hYJ8ukNCRE2TVKRJeSicPj5t8PngugrlCoWYUVvEHPMJaWiTkecXYAO07oIR1VTpE//PLB9MwsEmOfLJNsyHPZNUpEl5KJg1DdY/2FVZn+UJ+M1NqJPf4JGSDERU6WSVZq2ItnqkuksQyGiChtJp41OhAGgCH3LaKUalxAdG6RbwKumNeOfEtmKsK1denM/iUiSpuJK0IVhBnv5ZuU1J8L6t5+E3ST9zHeAzKl0vz4Css76kz8niQiMzLrh86gjkEdVg1Ws57BJVVTAABBHQDKL/4eHRPNpkWi2mnFMi6DIaJLy6wxwnLwolQbukMnqCIcpx9crQEosoG3HiSiSyzNT53Ozs7Gxkav17t69erOzs5Rjx84cMDlcqV6pTT1hwHAwSAcG1URKmnctfESWDdda6zR3MPfpoqIaIKkGYRbtmypra31+/01NTWPPvroyMe7urruvvvuvr6+YV4sHYGwBGAfxwSNnOJzilgf8pTLsiIE8MQy7VLu+0pEpKQZhM8+++z69evtdvv69eu3bds2wnFd1+++++7vfOc7mWlvVP8gwK7RMbNqWFAS+aPh8uwaBVBbIH54NX+jRHSppTnVr7W1tba2FoCq/0Y4/sgjj9TX199xxx0jv+CePXumTZuWfPyHP/zhihUrjEc6OjoCgcCp/nygyIbBtraO9E7BpILBYF9f3+Dg4MU+canbsa89XxNAd3tb70Q0LQOqgba2UR5TGAr1dnXJUR83SYVCod7e3nA4nO2GZE1XV1dfX19e3iW5t+flJxQK9fT06Lqe7YZkTVdXl8PhsNvtoz8UAOD1eq3WUZIuzSCUUgoh1BfG/yeHHG9ubn7ppZd27tw56gsuXrz4xz/+cfLxysrKIYOLUkq3233yvA3QnTZLSUlJeqdgUoFAwG63e73ei33i56fhXz7RvXaUl5n7iulWa1FRkTXHfu8xwWDQZrPl2tveSAjhcrnG/jk4yYRCIavVmstvAE3THA5Hfv5Yt/DQtNE7PtMMQp/Pd+LEiVmzZrW2tlZVVQ13vLm5edeuXTabTf2rEGLv3r3Lly9PfkGXyzV79uyx/GiLxWKxWIJSA3SnTVgsudWZZom62Cde74NN0yscpr9iUghN08x+FmlL+w0waeT4FdB1PZdPHxPzBkhzjLCxsbGpqUlK2dTUtGbNGgC7d+9OPr5582YZBUBKmTIF08DlExerwIZFpeKyHSAkIsqWNINw48aN+/fvr66uPnjw4IMPPghg5cqVKY9PkP5BCYxrW68ctLJSVHAbTyKiRGl2jXo8nh07dhiPqJov+fiQB2QKK8I0rPRpfziZu2PsREQpmXUbDwZhGq4tF9Xco5yIKJFZgzCggtDEN8/IAqcVX6hhEBIRJTBrEPazIkwLb/5ORDSEWYNwgHuNEhFRJpg7CFkREhHROJk1CLnpNhERZYRZg5AVIRERZQSDkIiIcpppg1DdhonLJ4iIaHzMGoRcPkFERBlh1iBk1ygREWWEeYOQm24TEVEGmDUIA6wIiYgoE8wahP1qsgyDkIiIxsesQcgxQiIiyghzB6GDyyeIiGh8zB2ErAiJiGiczBuEnDVKREQZYN4gBFgREhHRuJk1CDlrlIiIMsKUQSij6wjtDEIiIhofUwZhIAwJ2C3QOERIRETjY8og5AAhERFlCoOQiIhymimDMLrRKDtGiYhovEwZhNxWhoiIMsXEQciuUSIiGj8GIRER5TQGIRER5TQGIRER5TRTByFnjRIR0XiZOQg5a5SIiMbNzEHIrlEiIho3BiEREeU0Ewehg0FIRETjZuIgZEVIRETjZ8ogDOoAYOesUSIiGjdTBiErQiIiyhRTBmH/IMAgJCKiTDBlEAa4jpCIiDLEnEGoxghN2XYiIrq8mDJMVEWYx65RIiIaN1MGoZo1mmfKthMR0eXFlGESDAOAnRUhERGNmymDMCQBIE/jOkIiIhovUwZhkGOERESUIeYMQo4REhFRhpgyTBiERESUKaYME3aNEhFRppgyCAOsCImIKENMGSaRitCUbSciosuLKcMkMkbIrlEiIho3UwahWkdo5zpCIiIaN1MGISfLEBFRppgzCDlZhoiIMsSUYcIgJCKiTDFfmOgSgzo0Aav52k5ERJcd84VJSAqwHCQiogwxX56wX5SIiDIozTzp7OxsbGz0er2rV6/u7Owc+fi2bdvmzZvn8Xiuu+66jz76aJwtDuoCnDJKREQZkmYQbtmypba21u/319TUPProoyMcP3bs2D333NPU1OT3+xsbG++9995xtjgYluDNCImIKEOs6T3t2Weffe655+x2+/r169esWfPII48Md/zo0aN33nnn0qVLAdxzzz1btmxJ+YJ+v/+Xv/xl8vFrrrmmoqLCeKR3IAggT5MDAwPpNd7UAoFAIBDIzXNXNF0PBYPhXL0CwWAwx98AgUDAYrFIKbPdkOwIhUJ8AwhxEYVQXl6epo1S8qUZhK2trbW1tQBU/TfC8Ztuuummm24CEA6HN27cuG7dupQvOFwQlpeXu91u4xEVhFbo/f396TXe1AKBQH9/f35+frYbkjUOXQ8EAiInf/sAgsFgf39/br75lf7+fiGEruvZbkh2hEIhvgEAjP0vIZvNNupj0gxCKaXKZCllOBwe9fjOnTs3bNiwatWqTZs2pXzBhQsXPvPMM2P50fmdAODIsxQXF6fXeFMLBAJ5eXm5ee7KoNVaWFhoy9UrEAwGrVZrLr8BpJQul8tut2e7IdkRCoU0TcvlNwAAh8OR2WIgzSD0+XwnTpyYNWtWa2trVVXVCMellBs2bHjjjTe2bt1aX18//hZz1igREWVQmnnS2NjY1NQkpWxqalqzZg2A3bt3pzz+yiuvbN++/fnnn/f5fL29vb29veNsMYOQiIgyKM082bhx4/79+6urqw8ePPjggw8CWLlyZcrju3fvPnz4cHFxcWHUOFsc4vIJIiLKnDS7Rj0ez44dO4xH1NBl8vGHHnrooYceSrt9yXh7eiIiyiDz5UmkIjRfw4mI6HJkvjwJ6hJAnoUL6omIKAPMF4SqIrRzjJCIiDLBfEHIWaNERJRB5ssT3oaJiIgyyHx5EgwDXD5BREQZYr4gZEVIREQZZL48qXWG76jFVSWcNUpERBmQ5oL6LLq1cnBdA2w280U4ERFdhhgnRESU0xiERESU0xiERESU08wXhE899dTp06ez3YqsOXz48K9+9atstyKb2tra2trast2KrDl27NjPfvazbLcim37xi1988skn2W5F1pw6derpp5/Odiuy6Zlnnvnwww8z+5rmC8Kf//zn7e3t2W5F1hw5cuS5557Ldiuy6ezZs7kchCdOnPjNb36T7VZk07Zt244fP57tVmSN3+/P8T+Ff/e733388ceZfU3zBSEREVEGMQiJiCinMQiJiCinCXVn+ez6/e9/v3Hjxrlz547lwc3NzUuWLHG73RPdqsuT3+8/evTotddem+2GZM2X//3f/7BgQdf06dluSHa0tbW1tLRcf/312W5I1uzdu7e+vr68vDzbDcmOzs7O995778Ybb8x2Q7LmjTfeqK6unjp16hgfv2nTppqampEfc1kEIYDt27d3dXVluxVERDSpNDY2FhcXj/yYyyUIiYiIsoJjhERElNMYhERElNMYhERElNNME4SdnZ2NjY1er3f16tWdnZ3Zbs6lEA6HGxoa1NfJpz+5L8i2bdvmzZvn8Xiuu+66jz76CLl3BV544YW5c+d6PJ65c+fu3LkTuXcFABw4cMDlcqmvc+30ly1bJqLuu+8+5NgVGBwcvP/++8vKyq699trW1lZM8OmbJgi3bNlSW1vr9/tramoeffTRbDdnwj3++OPLli07fPiw+jb59CfxBTl27Ng999zT1NTk9/sbGxvvvfde5NgVCIfDd9555xNPPNHR0fH9738/B68AgK6urrvvvruvr099m1OnL6U8fPiw3+/v6enp6el57LHHkGNX4LHHHuvu7v7000+XLVv20EMPYaJPX5pEfX19S0uLlLKlpaW+vj7bzZlwL7/88vPPPx/7BSWf/iS+IM3NzV/72tfU121tbSUlJTLHrsDAwMCOHTt0Xe/u7t6+ffvcuXNljl2BcDi8evVqtauqOpJTp9/a2lpUVLRw4cKCgoI1a9acOXNG5tgVWLBgwZ///GcpZXd399tvvy0n+PRNE4Qul6uvr09K2dfXV1hYmO3mXCKxT4Hk08+FCzI4OHjffffdf//9MievQE9PDwAhxGuvvSZz7Aps2rTpgQcekLn6v8Dbb7+9cuXK995779y5c3fdddeXv/xlmWNXwOv1btiwobi4eNGiRfv375cTfPqm6RqVUgoh1BfhcDjbzbnUkk9/0l+QnTt3Ll682O12P/7448jJK1BQUNDb27tp06avf/3ryKUr0Nzc/NJLLz3yyCPGg7lz+gAWLVr08ssvX3XVVV6vd/PmzS+++CJy7Ap0d3dLKQ8ePHjrrbd+9atfxQSfvmmC0OfznThxAkBra2tVVVW2m3OpJZ/+JL4gUspvfetbDz/88NatWzdv3my1WpFjV+Do0aMbNmwA4HK5vvKVr7S0tCCXrkBzc/OuXbtsNpv6mBNCvPrqq7lz+gDeeeed119/XX2dl5dnt9uRS28AAGVlZd/4xjcqKyvXr19/4MABTPDpmyYIGxsbm5qapJRNTU1r1qzJdnMuteTTn8QX5JVXXtm+ffvzzz/v8/l6e3t7e3uRY1fA5/M99dRTe/fulVJu3bp1wYIFyKUrsHnz5lifFQAp5fLly3Pn9AFcuHBh7dq1LS0twWDw4Ycfvv3225FLbwAAq1at+ulPfxoIBJ566qnFixdjok9/PP2ql1JnZ+dtt91WVVXV2NjY1dWV7eZcIrFfUPLpT+IL8r3vfS/5XZpTV0BKuWvXroULFxYXF19zzTVqRkCuXQElN/8X0HX9iSeemDFjRmlp6V133XX+/HmZY1fA7/fffPPNbrf7uuuu+/jjj+UEnz73GiUiopxmmq5RIiKiicAgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinMYgJCKinPb/AVm0zwjsncZLAAAAAElFTkSuQmCC"
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determinism and random objects.\n",
    "using Random\n",
    "Random.seed!(42)\n",
    "\n",
    "# Set up model and compute the function w̄\n",
    "sp = SearchProblem(π_grid_size = 50, F_a = 1, F_b = 1)\n",
    "ϕ_init = ones(sp.n_π)\n",
    "g(x) = res_wage_operator(sp, x)\n",
    "w̄_vals = compute_fixed_point(g, ϕ_init)\n",
    "w̄ = extrapolate(interpolate((sp.π_grid, ), w̄_vals,\n",
    "                    Gridded(Linear())), Flat())\n",
    "\n",
    "# Holds the employment state and beliefs of an individual agent.\n",
    "mutable struct Agent{TF <: AbstractFloat, TI <: Integer}\n",
    "    _π::TF\n",
    "    employed::TI\n",
    "end\n",
    "\n",
    "Agent(_π=1e-3) = Agent(_π, 1)\n",
    "\n",
    "function update!(ag, H)\n",
    "    if ag.employed == 0\n",
    "        w = rand(H) * 2   # account for scale in julia\n",
    "        if w ≥ w̄(ag._π)\n",
    "            ag.employed = 1\n",
    "        else\n",
    "            ag._π = 1.0 ./ (1 .+ ((1 - ag._π) .* sp.g(w)) ./ (ag._π * sp.f(w)))\n",
    "        end\n",
    "    end\n",
    "    nothing\n",
    "end\n",
    "\n",
    "num_agents = 5000\n",
    "separation_rate = 0.025  # Fraction of jobs that end in each period\n",
    "separation_num = round(Int, num_agents * separation_rate)\n",
    "agent_indices = collect(1:num_agents)\n",
    "agents = [Agent() for i=1:num_agents]\n",
    "sim_length = 600\n",
    "H = sp.G                 # Start with distribution G\n",
    "change_date = 200        # Change to F after this many periods\n",
    "unempl_rate = zeros(sim_length)\n",
    "\n",
    "for i in 1:sim_length\n",
    "    if i % 20 == 0\n",
    "        println(\"date = $i\")\n",
    "    end\n",
    "\n",
    "    if i == change_date\n",
    "        H = sp.F\n",
    "    end\n",
    "\n",
    "    # Randomly select separation_num agents and set employment status to 0\n",
    "    shuffle!(agent_indices)\n",
    "    separation_list = agent_indices[1:separation_num]\n",
    "\n",
    "    for agent in agents[separation_list]\n",
    "        agent.employed = 0\n",
    "    end\n",
    "\n",
    "    # update agents\n",
    "    for agent in agents\n",
    "        update!(agent, H)\n",
    "    end\n",
    "    employed = Int[agent.employed for agent in agents]\n",
    "    unempl_rate[i] = 1.0 - mean(employed)\n",
    "end\n",
    "\n",
    "plot(unempl_rate, linewidth = 2, label = \"unemployment rate\")\n",
    "vline!([change_date], color = :red, label = \"\")"
   ]
  }
 ],
 "metadata": {
  "filename": "odu.rst",
  "kernelspec": {
   "display_name": "Julia 1.2",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  },
  "title": "Job Search III: Search with Learning"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}