odu.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='odu-v2'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Job Search V: Search with Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Job Search V: Search with Learning](#Job-Search-V:-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",
    "  - [Another Functional Equation](#Another-Functional-Equation)  \n",
    "  - [Solving the RWFE](#Solving-the-RWFE)  \n",
    "  - [Implementation](#Implementation)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  \n",
    "  - [Appendix A](#Appendix-A)  \n",
    "  - [Appendix B](#Appendix-B)  \n",
    "  - [Examples](#Examples)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to what’s in Anaconda, this lecture deploys the libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "  !pip install --upgrade quantecon\n",
    "  !pip install interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this lecture, we consider an extension of the [previously\n",
    "studied](https://python.quantecon.org/mccall_model.html)  job search model of McCall\n",
    "[[McC70]](https://python.quantecon.org/zreferences.html#mccall1970).\n",
    "\n",
    "We’ll build on a model of Bayesian learning discussed in [this\n",
    "lecture](https://python.quantecon.org/exchangeable.html) on the topic of exchangeability and its relationship to\n",
    "the concept of IID (identically and independently distributed) random variables and to  Bayesian updating.\n",
    "\n",
    "In the McCall model, an unemployed worker decides when to accept a\n",
    "permanent job at a specific fixed wage, given\n",
    "\n",
    "- his or her discount factor  \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\n",
    "must be learned.\n",
    "\n",
    "- The following is based on the presentation in\n",
    "  [[LS18]](https://python.quantecon.org/zreferences.html#ljungqvist2012), section 6.6.  \n",
    "\n",
    "\n",
    "Let’s start with some imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from numba import njit, prange, vectorize\n",
    "from interpolation import mlinterp, interp\n",
    "from math import gamma\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from matplotlib import cm\n",
    "import scipy.optimize as op\n",
    "from scipy.stats import cumfreq, beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Features\n",
    "\n",
    "- Infinite horizon dynamic programming with two states and one binary\n",
    "  control.  \n",
    "- Bayesian updating to learn the unknown distribution.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model\n",
    "\n",
    "Let’s first review the basic McCall model\n",
    "[[McC70]](https://python.quantecon.org/zreferences.html#mccall1970) and then add the variation we\n",
    "want to consider."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Basic McCall Model\n",
    "\n",
    "Recall that, in the baseline model <mccall_model>, an\n",
    "unemployed worker is presented in each period with a permanent job offer\n",
    "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\n",
    "  reconsiders next period  \n",
    "\n",
    "\n",
    "The wage sequence $ {W_t} $ is IID and generated from known density $ q $.\n",
    "\n",
    "The worker aims to maximize the expected discounted sum of earnings\n",
    "$ \\mathbb{E} \\sum_{t=0}^{\\infty}\\beta^t y_t $ 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')q(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 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\n",
    "[[LS18]](https://python.quantecon.org/zreferences.html#ljungqvist2012), section 6.6.\n",
    "\n",
    "The model is as above, apart from the fact that\n",
    "\n",
    "- the density $ q $ is unknown  \n",
    "- the worker learns about $ q $ by starting with a prior and updating\n",
    "  based on wage offers that he/she observes  \n",
    "\n",
    "\n",
    "The worker knows there are two possible distributions $ F $ and $ G $ —\n",
    "with densities $ f $ and $ g $.\n",
    "\n",
    "At the start of time, “nature” selects $ q $ to be either $ f $ or $ g $\n",
    "— the wage distribution from which the entire sequence $ {W_t} $ will be\n",
    "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 $,\n",
    "where $ \\pi_t $ updates via\n",
    "\n",
    "\n",
    "<a id='equation-odu-pi-rec-2'></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}\\{q = f \\,|\\, W = w\\}\n",
    "= \\frac{\\mathbb{P}\\{W = w \\,|\\, q = f\\}\\mathbb{P}\\{q = f\\}}\n",
    "{\\mathbb{P}\\{W = w\\}}\n",
    "\\quad \\text{and} \\quad\n",
    "\\mathbb{P}\\{W = w\\} = \\sum_{\\omega \\in \\{f, g\\}} \\mathbb{P}\\{W = w \\,|\\, q = \\omega\\} \\mathbb{P}\\{q = \\omega\\}\n",
    "$$\n",
    "\n",
    "The fact that [(2)](#equation-odu-pi-rec-2) is recursive allows us to\n",
    "progress to a recursive solution method.\n",
    "\n",
    "Letting\n",
    "\n",
    "$$\n",
    "q_{\\pi}(w) := \\pi f(w) + (1 - \\pi) g(w)\n",
    "\\quad \\text{and} \\quad\n",
    "\\kappa(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\n",
    "as 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') \\, q_{\\pi}(w') \\, dw'\n",
    "\\right\\}\n",
    "\\quad \\text{where} \\quad\n",
    "\\pi' = \\kappa(w', \\pi) \\tag{3}\n",
    "$$\n",
    "\n",
    "Notice that the current guess $ \\pi $ is a state variable,\n",
    "since it affects the worker’s perception of probabilities for future\n",
    "rewards."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameterization\n",
    "\n",
    "Following section 6.6 of [[LS18]](https://python.quantecon.org/zreferences.html#ljungqvist2012),\n",
    "our baseline parameterization will be\n",
    "\n",
    "- $ f $ is $ \\operatorname{Beta}(1, 1) $  \n",
    "- $ g $ is $ \\operatorname{Beta}(3, 1.2) $  \n",
    "- $ \\beta = 0.95 $ and $ c = 0.3 $  \n",
    "\n",
    "\n",
    "The densities $ f $ and $ g $ have the following shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@vectorize\n",
    "def p(x, a, b):\n",
    "    r = gamma(a + b) / (gamma(a) * gamma(b))\n",
    "    return r * x**(a-1) * (1 - x)**(b-1)\n",
    "\n",
    "\n",
    "x_grid = np.linspace(0, 1, 100)\n",
    "f = lambda x: p(x, 1, 1)\n",
    "g = lambda x: p(x, 3, 1.2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 8))\n",
    "ax.plot(x_grid, f(x_grid), label='$f$', lw=2)\n",
    "ax.plot(x_grid, g(x_grid), label='$g$', lw=2)\n",
    "\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looking Forward\n",
    "\n",
    "What kind of optimal policy might result from\n",
    "[(3)](#equation-odu-mvf) and the parameterization specified above?\n",
    "\n",
    "Intuitively, if we accept at $ w_a $ and $ w_a\\leq w_b $,\n",
    "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\n",
    "threshold value $ \\bar w $.\n",
    "\n",
    "But $ \\bar w $ should depend on $ \\pi $ — in\n",
    "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\n",
    "  $ g $  \n",
    "\n",
    "\n",
    "Thus larger $ \\pi $ depresses the worker’s assessment of\n",
    "her future prospects, and relatively low current offers become more\n",
    "attractive.\n",
    "\n",
    "**Summary:** We conjecture that the optimal policy is of the form\n",
    "$ \\mathbb 1{w\\geq \\bar w(\\pi) } $ for some\n",
    "decreasing function $ \\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\n",
    "intuition.\n",
    "\n",
    "We begin by solving via value function iteration (VFI), which is natural\n",
    "but ultimately turns out to be second best.\n",
    "\n",
    "The class `SearchProblem` is used to store parameters and methods\n",
    "needed to compute optimal actions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "class SearchProblem:\n",
    "    \"\"\"\n",
    "    A class to store a given parameterization of the \"offer distribution\n",
    "    unknown\" model.\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self,\n",
    "                 β=0.95,            # Discount factor\n",
    "                 c=0.3,             # Unemployment compensation\n",
    "                 F_a=1,\n",
    "                 F_b=1,\n",
    "                 G_a=3,\n",
    "                 G_b=1.2,\n",
    "                 w_max=1,           # Maximum wage possible\n",
    "                 w_grid_size=100,\n",
    "                 π_grid_size=100,\n",
    "                 mc_size=500):\n",
    "\n",
    "        self.β, self.c, self.w_max = β, c, w_max\n",
    "\n",
    "        self.f = njit(lambda x: p(x, F_a, F_b))\n",
    "        self.g = njit(lambda x: p(x, G_a, G_b))\n",
    "\n",
    "        self.π_min, self.π_max = 1e-3, 1-1e-3    # Avoids instability\n",
    "        self.w_grid = np.linspace(0, w_max, w_grid_size)\n",
    "        self.π_grid = np.linspace(self.π_min, self.π_max, π_grid_size)\n",
    "\n",
    "        self.mc_size = mc_size\n",
    "\n",
    "        self.w_f = np.random.beta(F_a, F_b, mc_size)\n",
    "        self.w_g = np.random.beta(G_a, G_b, mc_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes an instance of this class and returns\n",
    "jitted versions of the Bellman operator `T`, and a `get_greedy()`\n",
    "function to compute the approximate optimal policy from a guess `v` of\n",
    "the value function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def operator_factory(sp, parallel_flag=True):\n",
    "\n",
    "    f, g = sp.f, sp.g\n",
    "    w_f, w_g = sp.w_f, sp.w_g\n",
    "    β, c = sp.β, sp.c\n",
    "    mc_size = sp.mc_size\n",
    "    w_grid, π_grid = sp.w_grid, sp.π_grid\n",
    "\n",
    "    @njit\n",
    "    def κ(w, π):\n",
    "        \"\"\"\n",
    "        Updates π using Bayes' rule and the current wage observation w.\n",
    "        \"\"\"\n",
    "        pf, pg = π * f(w), (1 - π) * g(w)\n",
    "        π_new = pf / (pf + pg)\n",
    "\n",
    "        return π_new\n",
    "\n",
    "    @njit(parallel=parallel_flag)\n",
    "    def T(v):\n",
    "        \"\"\"\n",
    "        The Bellman operator.\n",
    "\n",
    "        \"\"\"\n",
    "        v_func = lambda x, y: mlinterp((w_grid, π_grid), v, (x, y))\n",
    "        v_new = np.empty_like(v)\n",
    "\n",
    "        for i in prange(len(w_grid)):\n",
    "            for j in prange(len(π_grid)):\n",
    "                w = w_grid[i]\n",
    "                π = π_grid[j]\n",
    "\n",
    "                v_1 = w / (1 - β)\n",
    "\n",
    "                integral_f, integral_g = 0.0, 0.0\n",
    "                for m in prange(mc_size):\n",
    "                    integral_f += v_func(w_f[m], κ(w_f[m], π))\n",
    "                    integral_g += v_func(w_g[m], κ(w_g[m], π))\n",
    "                integral = (π * integral_f + (1 - π) * integral_g) / mc_size\n",
    "\n",
    "                v_2 = c + β * integral\n",
    "                v_new[i, j] = max(v_1, v_2)\n",
    "\n",
    "        return v_new\n",
    "\n",
    "    @njit(parallel=parallel_flag)\n",
    "    def get_greedy(v):\n",
    "        \"\"\"\"\n",
    "        Compute optimal actions taking v as the value function.\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        v_func = lambda x, y: mlinterp((w_grid, π_grid), v, (x, y))\n",
    "        σ = np.empty_like(v)\n",
    "\n",
    "        for i in prange(len(w_grid)):\n",
    "            for j in prange(len(π_grid)):\n",
    "                w = w_grid[i]\n",
    "                π = π_grid[j]\n",
    "\n",
    "                v_1 = w / (1 - β)\n",
    "\n",
    "                integral_f, integral_g = 0.0, 0.0\n",
    "                for m in prange(mc_size):\n",
    "                    integral_f += v_func(w_f[m], κ(w_f[m], π))\n",
    "                    integral_g += v_func(w_g[m], κ(w_g[m], π))\n",
    "                integral = (π * integral_f + (1 - π) * integral_g) / mc_size\n",
    "\n",
    "                v_2 = c + β * integral\n",
    "\n",
    "                σ[i, j] = v_1 > v_2  # Evaluates to 1 or 0\n",
    "\n",
    "        return σ\n",
    "\n",
    "    return T, get_greedy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will omit a detailed discussion of the code because there is a more\n",
    "efficient solution method that we will use later.\n",
    "\n",
    "To solve the model we will use the following function that iterates\n",
    "using T to find a fixed point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def solve_model(sp,\n",
    "                use_parallel=True,\n",
    "                tol=1e-4,\n",
    "                max_iter=1000,\n",
    "                verbose=True,\n",
    "                print_skip=5):\n",
    "\n",
    "    \"\"\"\n",
    "    Solves for the value function\n",
    "\n",
    "    * sp is an instance of SearchProblem\n",
    "    \"\"\"\n",
    "\n",
    "    T, _ = operator_factory(sp, use_parallel)\n",
    "\n",
    "    # Set up loop\n",
    "    i = 0\n",
    "    error = tol + 1\n",
    "    m, n = len(sp.w_grid), len(sp.π_grid)\n",
    "\n",
    "    # Initialize v\n",
    "    v = np.zeros((m, n)) + sp.c / (1 - sp.β)\n",
    "\n",
    "    while i < max_iter and error > tol:\n",
    "        v_new = T(v)\n",
    "        error = np.max(np.abs(v - v_new))\n",
    "        i += 1\n",
    "        if verbose and i % print_skip == 0:\n",
    "            print(f\"Error at iteration {i} is {error}.\")\n",
    "        v = v_new\n",
    "\n",
    "    if i == max_iter:\n",
    "        print(\"Failed to converge!\")\n",
    "\n",
    "    if verbose and i < max_iter:\n",
    "        print(f\"\\nConverged in {i} iterations.\")\n",
    "\n",
    "\n",
    "    return v_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s look at solutions computed from value function iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "sp = SearchProblem()\n",
    "v_star = solve_model(sp)\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "ax.contourf(sp.π_grid, sp.w_grid, v_star, 12, alpha=0.6, cmap=cm.jet)\n",
    "cs = ax.contour(sp.π_grid, sp.w_grid, v_star, 12, colors=\"black\")\n",
    "ax.clabel(cs, inline=1, fontsize=10)\n",
    "ax.set(xlabel='$\\pi$', ylabel='$w$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also plot the optimal policy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "T, get_greedy = operator_factory(sp)\n",
    "σ_star = get_greedy(v_star)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "ax.contourf(sp.π_grid, sp.w_grid, σ_star, 1, alpha=0.6, cmap=cm.jet)\n",
    "ax.contour(sp.π_grid, sp.w_grid, σ_star, 1, colors=\"black\")\n",
    "ax.set(xlabel='$\\pi$', ylabel='$w$')\n",
    "\n",
    "ax.text(0.5, 0.6, 'reject')\n",
    "ax.text(0.7, 0.9, 'accept')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results fit well with our intuition from section [looking\n",
    "forward](#looking-forward).\n",
    "\n",
    "- The black line in the figure above corresponds to the function\n",
    "  $ \\bar w(\\pi) $ introduced there.  \n",
    "- It is decreasing as expected.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take 2: A More Efficient Method\n",
    "\n",
    "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\n",
    "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\n",
    "bit of 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') \\, q_{\\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",
    "\\, q_{\\pi}(w') \\, dw'\n",
    "$$\n",
    "\n",
    "Multiplying by $ 1 - \\beta $, substituting in $ \\pi' = \\kappa(w', \\pi) $\n",
    "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 \\kappa(w', \\pi) \\right\\} \\, q_{\\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",
    "- $ \\| \\omega \\| := \\sup_{x \\in [0,1]} | \\omega(x) | $  \n",
    "\n",
    "\n",
    "Consider the operator $ Q $ mapping $ \\omega \\in b[0,1] $ into $ Q\\omega \\in b[0,1] $ via\n",
    "\n",
    "\n",
    "<a id='equation-odu-dq'></a>\n",
    "$$\n",
    "(Q \\omega)(\\pi)\n",
    "= (1 - \\beta) c +\n",
    "\\beta \\int \\max \\left\\{ w', \\omega \\circ \\kappa(w', \\pi) \\right\\} \\, q_{\\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 $ \\omega, \\omega' \\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 \\omega)(\\pi) - (Q \\omega')(\\pi)|\n",
    "\\leq \\beta \\int\n",
    "\\left|\n",
    "\\max \\left\\{w', \\omega \\circ \\kappa(w', \\pi) \\right\\} -\n",
    "\\max \\left\\{w', \\omega' \\circ \\kappa(w', \\pi) \\right\\}\n",
    "\\right|\n",
    "\\, q_{\\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 \\omega)(\\pi) - (Q \\omega')(\\pi)|\n",
    "\\leq \\beta \\int\n",
    "\\left| \\omega \\circ \\kappa(w', \\pi) -  \\omega' \\circ \\kappa(w', \\pi) \\right|\n",
    "\\, q_{\\pi}(w') \\, dw'\n",
    "\\leq \\beta \\| \\omega - \\omega' \\| \\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 \\omega - Q \\omega'\\|\n",
    "\\leq \\beta \\| \\omega - \\omega' \\| \\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 \\omega \\to \\bar w $ uniformly as $ k \\to \\infty $, for any $ \\omega \\in b[0,1] $.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "The following function takes an instance of `SearchProblem` and\n",
    "returns the operator `Q`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def Q_factory(sp, parallel_flag=True):\n",
    "\n",
    "    f, g = sp.f, sp.g\n",
    "    w_f, w_g = sp.w_f, sp.w_g\n",
    "    β, c = sp.β, sp.c\n",
    "    mc_size = sp.mc_size\n",
    "    w_grid, π_grid = sp.w_grid, sp.π_grid\n",
    "\n",
    "    @njit\n",
    "    def κ(w, π):\n",
    "        \"\"\"\n",
    "        Updates π using Bayes' rule and the current wage observation w.\n",
    "        \"\"\"\n",
    "        pf, pg = π * f(w), (1 - π) * g(w)\n",
    "        π_new = pf / (pf + pg)\n",
    "\n",
    "        return π_new\n",
    "\n",
    "    @njit(parallel=parallel_flag)\n",
    "    def Q(ω):\n",
    "        \"\"\"\n",
    "\n",
    "        Updates the reservation wage function guess ω via the operator\n",
    "        Q.\n",
    "\n",
    "        \"\"\"\n",
    "        ω_func = lambda p: interp(π_grid, ω, p)\n",
    "        ω_new = np.empty_like(ω)\n",
    "\n",
    "        for i in prange(len(π_grid)):\n",
    "            π = π_grid[i]\n",
    "            integral_f, integral_g = 0.0, 0.0\n",
    "\n",
    "            for m in prange(mc_size):\n",
    "                integral_f += max(w_f[m], ω_func(κ(w_f[m], π)))\n",
    "                integral_g += max(w_g[m], ω_func(κ(w_g[m], π)))\n",
    "            integral = (π * integral_f + (1 - π) * integral_g) / mc_size\n",
    "\n",
    "            ω_new[i] = (1 - β) * c + β * integral\n",
    "\n",
    "        return ω_new\n",
    "\n",
    "    return Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next exercise, you are asked to compute an approximation to\n",
    "$ \\bar w $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Use the default parameters and `Q_factory` to compute an optimal\n",
    "policy.\n",
    "\n",
    "Your result should coincide closely with the figure for the optimal\n",
    "policy [shown above](#odu-pol-vfi).\n",
    "\n",
    "Try experimenting with different parameters, and confirm that the change\n",
    "in 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.\n",
    "\n",
    "You should find that the run time is shorter than that of the value\n",
    "function approach.\n",
    "\n",
    "Similar to above, we set up a function to iterate with `Q` to find the\n",
    "fixed point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def solve_wbar(sp,\n",
    "               use_parallel=True,\n",
    "               tol=1e-4,\n",
    "               max_iter=1000,\n",
    "               verbose=True,\n",
    "               print_skip=5):\n",
    "\n",
    "    Q = Q_factory(sp, use_parallel)\n",
    "\n",
    "    # Set up loop\n",
    "    i = 0\n",
    "    error = tol + 1\n",
    "    m, n = len(sp.w_grid), len(sp.π_grid)\n",
    "\n",
    "    # Initialize w\n",
    "    w = np.ones_like(sp.π_grid)\n",
    "\n",
    "    while i < max_iter and error > tol:\n",
    "        w_new = Q(w)\n",
    "        error = np.max(np.abs(w - w_new))\n",
    "        i += 1\n",
    "        if verbose and i % print_skip == 0:\n",
    "            print(f\"Error at iteration {i} is {error}.\")\n",
    "        w = w_new\n",
    "\n",
    "    if i == max_iter:\n",
    "        print(\"Failed to converge!\")\n",
    "\n",
    "    if verbose and i < max_iter:\n",
    "        print(f\"\\nConverged in {i} iterations.\")\n",
    "\n",
    "    return w_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution can be plotted as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "sp = SearchProblem()\n",
    "w_bar = solve_wbar(sp)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 7))\n",
    "\n",
    "ax.plot(sp.π_grid, w_bar, color='k')\n",
    "ax.fill_between(sp.π_grid, 0, w_bar, color='blue', alpha=0.15)\n",
    "ax.fill_between(sp.π_grid, w_bar, sp.w_max, color='green', alpha=0.15)\n",
    "ax.text(0.5, 0.6, 'reject')\n",
    "ax.text(0.7, 0.9, 'accept')\n",
    "ax.set(xlabel='$\\pi$', ylabel='$w$')\n",
    "ax.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix A\n",
    "\n",
    "The next piece of code generates a fun simulation to see what the effect\n",
    "of a change in the underlying distribution on the unemployment rate is.\n",
    "\n",
    "At a point in the simulation, the distribution becomes significantly\n",
    "worse.\n",
    "\n",
    "It takes a while for agents to learn this, and in the meantime, they are\n",
    "too optimistic and turn down too many jobs.\n",
    "\n",
    "As a result, the unemployment rate spikes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F_a, F_b, G_a, G_b = 1, 1, 3, 1.2\n",
    "\n",
    "sp = SearchProblem(F_a=F_a, F_b=F_b, G_a=G_a, G_b=G_b)\n",
    "f, g = sp.f, sp.g\n",
    "\n",
    "# Solve for reservation wage\n",
    "w_bar = solve_wbar(sp, verbose=False)\n",
    "\n",
    "# Interpolate reservation wage function\n",
    "π_grid = sp.π_grid\n",
    "w_func = njit(lambda x: interp(π_grid, w_bar, x))\n",
    "\n",
    "@njit\n",
    "def update(a, b, e, π):\n",
    "    \"Update e and π by drawing wage offer from beta distribution with parameters a and b\"\n",
    "\n",
    "    if e == False:\n",
    "        w = np.random.beta(a, b)       # Draw random wage\n",
    "        if w >= w_func(π):\n",
    "            e = True                   # Take new job\n",
    "        else:\n",
    "            π = 1 / (1 + ((1 - π) * g(w)) / (π * f(w)))\n",
    "\n",
    "    return e, π\n",
    "\n",
    "@njit\n",
    "def simulate_path(F_a=F_a,\n",
    "                  F_b=F_b,\n",
    "                  G_a=G_a,\n",
    "                  G_b=G_b,\n",
    "                  N=5000,       # Number of agents\n",
    "                  T=600,        # Simulation length\n",
    "                  d=200,        # Change date\n",
    "                  s=0.025):     # Separation rate\n",
    "\n",
    "    \"\"\"Simulates path of employment for N number of works over T periods\"\"\"\n",
    "\n",
    "    e = np.ones((N, T+1))\n",
    "    π = np.ones((N, T+1)) * 1e-3\n",
    "\n",
    "    a, b = G_a, G_b   # Initial distribution parameters\n",
    "\n",
    "    for t in range(T+1):\n",
    "\n",
    "        if t == d:\n",
    "            a, b = F_a, F_b  # Change distribution parameters\n",
    "\n",
    "        # Update each agent\n",
    "        for n in range(N):\n",
    "            if e[n, t] == 1:                    # If agent is currently employment\n",
    "                p = np.random.uniform(0, 1)\n",
    "                if p <= s:                      # Randomly separate with probability s\n",
    "                    e[n, t] = 0\n",
    "\n",
    "            new_e, new_π = update(a, b, e[n, t], π[n, t])\n",
    "            e[n, t+1] = new_e\n",
    "            π[n, t+1] = new_π\n",
    "\n",
    "    return e[:, 1:]\n",
    "\n",
    "d = 200  # Change distribution at time d\n",
    "unemployment_rate = 1 - simulate_path(d=d).mean(axis=0)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "ax.plot(unemployment_rate)\n",
    "ax.axvline(d, color='r', alpha=0.6, label='Change date')\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_title('Unemployment rate')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix B\n",
    "\n",
    "In this appendix we provide more details about how Bayes’ Law contributes to the workings of the model.\n",
    "\n",
    "We present some graphs that bring out additional insights about how learning works.\n",
    "\n",
    "We build on graphs proposed in [this lecture](https://python.quantecon.org/exchangeable.html).\n",
    "\n",
    "In particular, we’ll add actions of  our searching worker to a key graph\n",
    "presented in that lecture.\n",
    "\n",
    "To begin, we first define two functions for computing the\n",
    "empirical distributions of unemployment duration and π at the time of\n",
    "employment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def empirical_dist(F_a, F_b, G_a, G_b, w_bar, π_grid,\n",
    "                   N=10000, T=600):\n",
    "    \"\"\"\n",
    "    Simulates population for computing empirical cumulative\n",
    "    distribution of unemployment duration and π at time when\n",
    "    the worker accepts the wage offer. For each job searching\n",
    "    problem, we simulate for two cases that either f or g is\n",
    "    the true offer distribution.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "\n",
    "    F_a, F_b, G_a, G_b : parameters of beta distributions F and G.\n",
    "    w_bar : the reservation wage\n",
    "    π_grid : grid points of π, for interpolation\n",
    "    N : number of workers for simulation, optional\n",
    "    T : maximum of time periods for simulation, optional\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    accpet_t : 2 by N ndarray. the empirical distribution of\n",
    "               unemployment duration when f or g generates offers.\n",
    "    accept_π : 2 by N ndarray. the empirical distribution of\n",
    "               π at the time of employment when f or g generates offers.\n",
    "    \"\"\"\n",
    "\n",
    "    accept_t = np.empty((2, N))\n",
    "    accept_π = np.empty((2, N))\n",
    "\n",
    "    # f or g generates offers\n",
    "    for i, (a, b) in enumerate([(F_a, F_b), (G_a, G_b)]):\n",
    "        # update each agent\n",
    "        for n in range(N):\n",
    "\n",
    "            # initial priori\n",
    "            π = 0.5\n",
    "\n",
    "            for t in range(T+1):\n",
    "\n",
    "                # Draw random wage\n",
    "                w = np.random.beta(a, b)\n",
    "                lw = p(w, F_a, F_b) / p(w, G_a, G_b)\n",
    "                π = π * lw / (π * lw + 1 - π)\n",
    "\n",
    "                # move to next agent if accepts\n",
    "                if w >= interp(π_grid, w_bar, π):\n",
    "                    break\n",
    "\n",
    "            # record the unemployment duration\n",
    "            # and π at the time of acceptance\n",
    "            accept_t[i, n] = t\n",
    "            accept_π[i, n] = π\n",
    "\n",
    "    return accept_t, accept_π\n",
    "\n",
    "def cumfreq_x(res):\n",
    "    \"\"\"\n",
    "    A helper function for calculating the x grids of\n",
    "    the cumulative frequency histogram.\n",
    "    \"\"\"\n",
    "\n",
    "    cumcount = res.cumcount\n",
    "    lowerlimit, binsize = res.lowerlimit, res.binsize\n",
    "\n",
    "    x = lowerlimit + np.linspace(0, binsize*cumcount.size, cumcount.size)\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define a wrapper function for analyzing job search models with\n",
    "learning under different parameterizations.\n",
    "\n",
    "The wrapper takes parameters of beta distributions and  unemployment\n",
    "compensation as inputs  and then displays various things we want to know\n",
    "to interpret the solution of our search model.\n",
    "\n",
    "In addition, it computes empirical cumulative distributions of two key objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def job_search_example(F_a=1, F_b=1, G_a=3, G_b=1.2, c=0.3):\n",
    "    \"\"\"\n",
    "    Given the parameters that specify F and G distributions,\n",
    "    calculate and display the rejection and acceptance area,\n",
    "    the evolution of belief π, and the probability of accepting\n",
    "    an offer at different π level, and simulate and calculate\n",
    "    the empirical cumulative distribution of the duration of\n",
    "    unemployment and π at the time the worker accepts the offer.\n",
    "    \"\"\"\n",
    "\n",
    "    # construct a search problem\n",
    "    sp = SearchProblem(F_a=F_a, F_b=F_b, G_a=G_a, G_b=G_b, c=c)\n",
    "    f, g = sp.f, sp.g\n",
    "    π_grid = sp.π_grid\n",
    "\n",
    "    # Solve for reservation wage\n",
    "    w_bar = solve_wbar(sp, verbose=False)\n",
    "\n",
    "    # l(w) = f(w) / g(w)\n",
    "    l = lambda w: f(w) / g(w)\n",
    "    # objective function for solving l(w) = 1\n",
    "    obj = lambda w: l(w) - 1.\n",
    "\n",
    "    # the mode of beta distribution\n",
    "    # use this to divide w into two intervals for root finding\n",
    "    G_mode = (G_a - 1) / (G_a + G_b - 2)\n",
    "    roots = np.empty(2)\n",
    "    roots[0] = op.root_scalar(obj, bracket=[1e-10, G_mode]).root\n",
    "    roots[1] = op.root_scalar(obj, bracket=[G_mode, 1-1e-10]).root\n",
    "\n",
    "    fig, axs = plt.subplots(2, 2, figsize=(12, 9))\n",
    "\n",
    "    # part 1: display the details of the model settings and some results\n",
    "    w_grid = np.linspace(1e-12, 1-1e-12, 100)\n",
    "\n",
    "    axs[0, 0].plot(l(w_grid), w_grid, label='$l$', lw=2)\n",
    "    axs[0, 0].vlines(1., 0., 1., linestyle=\"--\")\n",
    "    axs[0, 0].hlines(roots, 0., 2., linestyle=\"--\")\n",
    "    axs[0, 0].set_xlim([0., 2.])\n",
    "    axs[0, 0].legend(loc=4)\n",
    "    axs[0, 0].set(xlabel='$l(w)=f(w)/g(w)$', ylabel='$w$')\n",
    "\n",
    "    axs[0, 1].plot(sp.π_grid, w_bar, color='k')\n",
    "    axs[0, 1].fill_between(sp.π_grid, 0, w_bar, color='blue', alpha=0.15)\n",
    "    axs[0, 1].fill_between(sp.π_grid, w_bar, sp.w_max, color='green', alpha=0.15)\n",
    "    axs[0, 1].text(0.5, 0.6, 'reject')\n",
    "    axs[0, 1].text(0.7, 0.9, 'accept')\n",
    "\n",
    "    W = np.arange(0.01, 0.99, 0.08)\n",
    "    Π = np.arange(0.01, 0.99, 0.08)\n",
    "\n",
    "    ΔW = np.zeros((len(W), len(Π)))\n",
    "    ΔΠ = np.empty((len(W), len(Π)))\n",
    "    for i, w in enumerate(W):\n",
    "        for j, π in enumerate(Π):\n",
    "            lw = l(w)\n",
    "            ΔΠ[i, j] = π * (lw / (π * lw + 1 - π) - 1)\n",
    "\n",
    "    q = axs[0, 1].quiver(Π, W, ΔΠ, ΔW, scale=2, color='r', alpha=0.8)\n",
    "\n",
    "    axs[0, 1].hlines(roots, 0., 1., linestyle=\"--\")\n",
    "    axs[0, 1].set(xlabel='$\\pi$', ylabel='$w$')\n",
    "    axs[0, 1].grid()\n",
    "\n",
    "    axs[1, 0].plot(f(x_grid), x_grid, label='$f$', lw=2)\n",
    "    axs[1, 0].plot(g(x_grid), x_grid, label='$g$', lw=2)\n",
    "    axs[1, 0].vlines(1., 0., 1., linestyle=\"--\")\n",
    "    axs[1, 0].hlines(roots, 0., 2., linestyle=\"--\")\n",
    "    axs[1, 0].legend(loc=4)\n",
    "    axs[1, 0].set(xlabel='$f(w), g(w)$', ylabel='$w$')\n",
    "\n",
    "    axs[1, 1].plot(sp.π_grid, 1 - beta.cdf(w_bar, F_a, F_b), label='$f$')\n",
    "    axs[1, 1].plot(sp.π_grid, 1 - beta.cdf(w_bar, G_a, G_b), label='$g$')\n",
    "    axs[1, 1].set_ylim([0., 1.])\n",
    "    axs[1, 1].grid()\n",
    "    axs[1, 1].legend(loc=4)\n",
    "    axs[1, 1].set(xlabel='$\\pi$', ylabel='$\\mathbb{P}\\{w > \\overline{w} (\\pi)\\}$')\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "    # part 2: simulate empirical cumulative distribution\n",
    "    accept_t, accept_π = empirical_dist(F_a, F_b, G_a, G_b, w_bar, π_grid)\n",
    "    N = accept_t.shape[1]\n",
    "\n",
    "    cfq_t_F = cumfreq(accept_t[0, :], numbins=100)\n",
    "    cfq_π_F = cumfreq(accept_π[0, :], numbins=100)\n",
    "\n",
    "    cfq_t_G = cumfreq(accept_t[1, :], numbins=100)\n",
    "    cfq_π_G = cumfreq(accept_π[1, :], numbins=100)\n",
    "\n",
    "    fig, axs = plt.subplots(2, 1, figsize=(12, 9))\n",
    "\n",
    "    axs[0].plot(cumfreq_x(cfq_t_F), cfq_t_F.cumcount/N, label=\"f generates\")\n",
    "    axs[0].plot(cumfreq_x(cfq_t_G), cfq_t_G.cumcount/N, label=\"g generates\")\n",
    "    axs[0].grid(linestyle='--')\n",
    "    axs[0].legend(loc=4)\n",
    "    axs[0].title.set_text('CDF of duration of unemployment')\n",
    "    axs[0].set(xlabel='time', ylabel='Prob(time)')\n",
    "\n",
    "    axs[1].plot(cumfreq_x(cfq_π_F), cfq_π_F.cumcount/N, label=\"f generates\")\n",
    "    axs[1].plot(cumfreq_x(cfq_π_G), cfq_π_G.cumcount/N, label=\"g generates\")\n",
    "    axs[1].grid(linestyle='--')\n",
    "    axs[1].legend(loc=4)\n",
    "    axs[1].title.set_text('CDF of π at time worker accepts wage and leaves unemployment')\n",
    "    axs[1].set(xlabel='π', ylabel='Prob(π)')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now provide some examples that provide insights about how the model works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1 (Baseline)\n",
    "\n",
    "$ F $ ~ Beta(1, 1), $ G $ ~ Beta(3, 1.2), $ c $=0.3.\n",
    "\n",
    "In the graphs below, the red arrows in the upper right figure show how $ \\pi_t $ is updated in response to the\n",
    "new information $ w_t $.\n",
    "\n",
    "Recall the following formula from [this lecture](https://python.quantecon.org/exchangeable.html)\n",
    "\n",
    "$$\n",
    "\\frac{\\pi_{t+1}}{\\pi_{t}}=\\frac{l\\left(w_{t+1}\\right)}{\\pi_{t}l\\left(w_{t+1}\\right)+\\left(1-\\pi_{t}\\right)}\\begin{cases}\n",
    ">1 & \\text{if }l\\left(w_{t+1}\\right)>1\\\\\n",
    "\\leq1 & \\text{if }l\\left(w_{t+1}\\right)\\leq1\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "The formula implies that  the direction of motion of $ \\pi_t $ is determined by the relationship between\n",
    "$ l(w_t) $ and $ 1 $.\n",
    "\n",
    "The magnitude is small if\n",
    "\n",
    "- $ l(w) $ is close to $ 1 $, which means the new $ w $ is\n",
    "  not very informative for distinguishing two distributions,  \n",
    "- $ \\pi_{t-1} $ is close to either $ 0 $ or $ 1 $, which\n",
    "  means the priori is strong.  \n",
    "\n",
    "\n",
    "Will an unemployed  worker accept an offer earlier or\n",
    "not, when the actual ruling distribution is $ g $ instead of\n",
    "$ f $?\n",
    "\n",
    "Two countervailing effects are at work.\n",
    "\n",
    "- if f generates successive wage offers, then $ w $ is more likely to be low, but\n",
    "  $ \\pi $ is moving up toward to 1, which lowers the reservation wage,\n",
    "  i.e., the worker becomes  less selective the longer he or she remains unemployed.  \n",
    "- if g generates wage offers, then $ w $ is more likely to be high, but\n",
    "  $ \\pi $ is moving downward toward 0, increasing the reservation wage, i.e., the worker becomes  more selective\n",
    "  the longer he or she remains unemployed.  \n",
    "\n",
    "\n",
    "Quantitatively, the lower right figure sheds light on which effect dominates in this example.\n",
    "\n",
    "It shows the probability that a previously unemployed  worker\n",
    "accepts an offer at different values of $ \\pi $ when $ f $ or $ g $ generates\n",
    "wage offers.\n",
    "\n",
    "That graph shows that for the particular $ f $ and $ g $ in this example, the\n",
    "worker is always more likely to accept an offer when $ f $ generates the data  even when $ \\pi $ is close to zero so\n",
    "that  the worker believes the true distribution is $ g $ and therefore is relatively\n",
    "more selective.\n",
    "\n",
    "The empirical cumulative distribution of the duration of\n",
    "unemployment verifies our conjecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "job_search_example()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2\n",
    "\n",
    "$ F $ ~ Beta(1, 1), $ G $ ~ Beta(1.2, 1.2), $ c $=0.3.\n",
    "\n",
    "Now $ G $ has the same mean as $ F $ with a smaller variance.\n",
    "\n",
    "Since the unemployment compensation $ c $ serves as a lower bound\n",
    "for bad wage offers, $ G $ is now an “inferior” distribution to\n",
    "$ F $.\n",
    "\n",
    "Consequently, we observe that the optimal policy\n",
    "$ \\overline{w}(\\pi) $ is increasing in $ \\pi $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "job_search_example(1, 1, 1.2, 1.2, 0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 3\n",
    "\n",
    "$ F $ ~ Beta(1, 1), $ G $ ~ Beta(2, 2), $ c $=0.3.\n",
    "\n",
    "If the variance of $ G $ is smaller, we observe in the result that\n",
    "$ G $ is even more “inferior” and the slope of\n",
    "$ \\overline{w}(\\pi) $ is larger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "job_search_example(1, 1, 2, 2, 0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 4\n",
    "\n",
    "$ F $ ~ Beta(1, 1), $ G $ ~ Beta(3, 1.2), and $ c $=0.8.\n",
    "\n",
    "In this example, we keep the parameters of beta distributions to be the\n",
    "same with the baseline case but increase the unemployment compensation\n",
    "$ c $.\n",
    "\n",
    "Comparing outcomes to the baseline case (example 1) in which\n",
    "unemployment compensation if low ($ c $=0.3), now the worker can\n",
    "afford a longer  learning period.\n",
    "\n",
    "As a result, the worker tends to accept\n",
    "wage offers much later.\n",
    "\n",
    "Furthermore, at the time of accepting employment, the belief\n",
    "$ \\pi $ is closer to either $ 0 $ or $ 1 $.\n",
    "\n",
    "That means that\n",
    "the worker has a better  idea about what the true distribution is\n",
    "when he eventually chooses to accept a wage offer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "job_search_example(1, 1, 3, 1.2, c=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 5\n",
    "\n",
    "$ F $ ~ Beta(1, 1), $ G $ ~ Beta(3, 1.2), and $ c $=0.1.\n",
    "\n",
    "As expected, a smaller $ c $ makes an unemployed worker  accept wage offers earlier\n",
    "after having acquired less information about the wage distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "job_search_example(1, 1, 3, 1.2, c=0.1)"
   ]
  }
 ],
 "metadata": {
  "date": 1584334748.8340664,
  "filename": "odu.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Job Search V: Search with Learning"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}