mccall_fitted_vfi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<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 III: Fitted Value Function Iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Job Search III: Fitted Value Function Iteration](#Job-Search-III:-Fitted-Value-Function-Iteration)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [The Algorithm](#The-Algorithm)  \n",
    "  - [Implementation](#Implementation)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon\n",
    "!pip install --upgrade interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this lecture we again study the [McCall job search model with separation](https://python.quantecon.org/mccall_model_with_separation.html), but now with a continuous wage distribution.\n",
    "\n",
    "While we already considered continuous wage distributions briefly in the\n",
    "exercises of the [first job search lecture](https://python.quantecon.org/mccall_model.html),\n",
    "the change was relatively trivial in that case.\n",
    "\n",
    "This is because we were able to reduce the problem to solving for a single\n",
    "scalar value (the continuation value).\n",
    "\n",
    "Here, with separation, the change is less trivial, since a continuous wage distribution leads to an uncountably infinite state space.\n",
    "\n",
    "The infinite state space leads to additional challenges, particularly when it\n",
    "comes to applying value function iteration (VFI).\n",
    "\n",
    "These challenges will lead us to modify VFI by adding an interpolation step.\n",
    "\n",
    "The combination of VFI and this interpolation step is called **fitted value function iteration** (fitted VFI).\n",
    "\n",
    "Fitted VFI is very common in practice, so we will take some time to work through the details.\n",
    "\n",
    "We will use the following imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import quantecon as qe\n",
    "from interpolation import interp\n",
    "from numpy.random import randn\n",
    "from numba import njit, jitclass, prange, float64, int32"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Algorithm\n",
    "\n",
    "The model is the same as the McCall model with job separation we [studied before](https://python.quantecon.org/mccall_model_with_separation.html), except that the wage offer distribution is continuous.\n",
    "\n",
    "We are going to start with the two Bellman equations we obtained for the model with job separation after [a simplifying transformation](https://python.quantecon.org/mccall_model_with_separation.html#ast-mcm).\n",
    "\n",
    "Modified to accommodate continuous wage draws, they take the following form:\n",
    "\n",
    "\n",
    "<a id='equation-bell1mcmc'></a>\n",
    "$$\n",
    "d = \\int \\max \\left\\{ v(w'), \\,  u(c) + \\beta d \\right\\} q(w') d w' \\tag{1}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-bell2mcmc'></a>\n",
    "$$\n",
    "v(w) = u(w) + \\beta\n",
    "    \\left[\n",
    "        (1-\\alpha)v(w) + \\alpha d\n",
    "    \\right] \\tag{2}\n",
    "$$\n",
    "\n",
    "The unknowns here are the function $ v $ and the scalar $ d $.\n",
    "\n",
    "The difference between these and the pair of Bellman equations we previously worked on are\n",
    "\n",
    "1. in [(1)](#equation-bell1mcmc), what used to be a sum over a finite number of wage values is an integral over an infinite set.  \n",
    "1. The function $ v $ in [(2)](#equation-bell2mcmc) is defined over all $ w \\in \\mathbb R_+ $.  \n",
    "\n",
    "\n",
    "The function $ q $ in [(1)](#equation-bell1mcmc) is the density of the wage offer distribution.\n",
    "\n",
    "Its support is taken as equal to $ \\mathbb R_+ $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Value Function Iteration\n",
    "\n",
    "In theory, we should now proceed as follows:\n",
    "\n",
    "1. Begin with a guess $ v, d $ for the solutions to [(1)](#equation-bell1mcmc)–[(2)](#equation-bell2mcmc).  \n",
    "1. Plug $ v, d $ into the right hand side of [(1)](#equation-bell1mcmc)–[(2)](#equation-bell2mcmc) and\n",
    "  compute the left hand side to obtain updates $ v', d' $  \n",
    "1. Unless some stopping condition is satisfied, set $ (v, d) = (v', d') $\n",
    "  and go to step 2.  \n",
    "\n",
    "\n",
    "However, there is a problem we must confront before we implement this procedure:\n",
    "The iterates of the value function can neither be calculated exactly nor stored on a computer.\n",
    "\n",
    "To see the issue, consider [(2)](#equation-bell2mcmc).\n",
    "\n",
    "Even if $ v $ is a known function,  the only way to store its update $ v' $\n",
    "is to record its value $ v'(w) $ for every $ w \\in \\mathbb R_+ $.\n",
    "\n",
    "Clearly, this is impossible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitted Value Function Iteration\n",
    "\n",
    "What we will do instead is use **fitted value function iteration**.\n",
    "\n",
    "The procedure is as follows:\n",
    "\n",
    "Let a current guess $ v $ be given.\n",
    "\n",
    "Now we record the value of the function $ v' $ at only\n",
    "finitely many “grid” points $ w_1 < w_2 < \\cdots < w_I $ and then reconstruct $ v' $ from this information when required.\n",
    "\n",
    "More precisely, the algorithm will be\n",
    "\n",
    "\n",
    "<a id='fvi-alg'></a>\n",
    "1. Begin with an array $ \\mathbf v $ representing the values of an initial guess of the value function on some grid points $ \\{w_i\\} $.  \n",
    "1. Build a function $ v $ on the state space $ \\mathbb R_+ $ by interpolation or approximation, based on $ \\mathbf v $ and $ \\{ w_i\\} $.  \n",
    "1. Obtain and record the samples of the updated function $ v'(w_i) $ on each grid point $ w_i $.  \n",
    "1. Unless some stopping condition is satisfied, take this as the new array and go to step 1.  \n",
    "\n",
    "\n",
    "How should we go about step 2?\n",
    "\n",
    "This is a problem of function approximation, and there are many ways to approach it.\n",
    "\n",
    "What’s important here is that the function approximation scheme must not only\n",
    "produce a good approximation to each $ v $, but also that it combines well with the broader iteration algorithm described above.\n",
    "\n",
    "One good choice from both respects is continuous piecewise linear interpolation.\n",
    "\n",
    "This method\n",
    "\n",
    "1. combines well with value function iteration (see., e.g.,\n",
    "  [gordon1995stable] or [stachurski2008continuous]) and  \n",
    "1. preserves useful shape properties such as monotonicity and concavity/convexity.  \n",
    "\n",
    "\n",
    "Linear interpolation will be implemented using a JIT-aware Python interpolation library called [interpolation.py](https://github.com/EconForge/interpolation.py).\n",
    "\n",
    "The next figure illustrates piecewise linear interpolation of an arbitrary\n",
    "function on grid points $ 0, 0.2, 0.4, 0.6, 0.8, 1 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    y1 = 2 * np.cos(6 * x) + np.sin(14 * x)\n",
    "    return y1 + 2.5\n",
    "\n",
    "c_grid = np.linspace(0, 1, 6)\n",
    "f_grid = np.linspace(0, 1, 150)\n",
    "\n",
    "def Af(x):\n",
    "    return interp(c_grid, f(c_grid), x)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(f_grid, f(f_grid), 'b-', label='true function')\n",
    "ax.plot(f_grid, Af(f_grid), 'g-', label='linear approximation')\n",
    "ax.vlines(c_grid, c_grid * 0, f(c_grid), linestyle='dashed', alpha=0.5)\n",
    "\n",
    "ax.legend(loc=\"upper center\")\n",
    "\n",
    "ax.set(xlim=(0, 1), ylim=(0, 6))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "The first step is to build a jitted class for the McCall model with separation and\n",
    "a continuous wage offer distribution.\n",
    "\n",
    "We will take the utility function to be the log function for this application, with $ u(c) = \\ln c $.\n",
    "\n",
    "We will adopt the lognormal distribution for wages, with $ w = \\exp(\\mu + \\sigma z) $\n",
    "when $ z $ is standard normal and $ \\mu, \\sigma $ are parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def lognormal_draws(n=1000, μ=2.5, σ=0.5, seed=1234):\n",
    "    np.random.seed(seed)\n",
    "    z = np.random.randn(n)\n",
    "    w_draws = np.exp(μ + σ * z)\n",
    "    return w_draws"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here’s our class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mccall_data_continuous = [\n",
    "    ('c', float64),          # unemployment compensation\n",
    "    ('α', float64),          # job separation rate\n",
    "    ('β', float64),          # discount factor\n",
    "    ('σ', float64),          # scale parameter in lognormal distribution\n",
    "    ('μ', float64),          # location parameter in lognormal distribution\n",
    "    ('w_grid', float64[:]),  # grid of points for fitted VFI\n",
    "    ('w_draws', float64[:])  # draws of wages for Monte Carlo\n",
    "]\n",
    "\n",
    "@jitclass(mccall_data_continuous)\n",
    "class McCallModelContinuous:\n",
    "\n",
    "    def __init__(self,\n",
    "                 c=1,\n",
    "                 α=0.1,\n",
    "                 β=0.96,\n",
    "                 grid_min=1e-10,\n",
    "                 grid_max=5,\n",
    "                 grid_size=100,\n",
    "                 w_draws=lognormal_draws()):\n",
    "\n",
    "        self.c, self.α, self.β = c, α, β\n",
    "\n",
    "        self.w_grid = np.linspace(grid_min, grid_max, grid_size)\n",
    "        self.w_draws = w_draws\n",
    "\n",
    "    def update(self, v, d):\n",
    "\n",
    "        # Simplify names\n",
    "        c, α, β, σ, μ = self.c, self.α, self.β, self.σ, self.μ\n",
    "        w = self.w_grid\n",
    "        u = lambda x: np.log(x)\n",
    "\n",
    "        # Interpolate array represented value function\n",
    "        vf = lambda x: interp(w, v, x)\n",
    "\n",
    "        # Update d using Monte Carlo to evaluate integral\n",
    "        d_new = np.mean(np.maximum(vf(self.w_draws), u(c) + β * d))\n",
    "\n",
    "        # Update v\n",
    "        v_new = u(w) + β * ((1 - α) * v + α * d)\n",
    "\n",
    "        return v_new, d_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then return the current iterate as an approximate solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def solve_model(mcm, tol=1e-5, max_iter=2000):\n",
    "    \"\"\"\n",
    "    Iterates to convergence on the Bellman equations\n",
    "\n",
    "    * mcm is an instance of McCallModel\n",
    "    \"\"\"\n",
    "\n",
    "    v = np.ones_like(mcm.w_grid)    # Initial guess of v\n",
    "    d = 1                           # Initial guess of d\n",
    "    i = 0\n",
    "    error = tol + 1\n",
    "\n",
    "    while error > tol and i < max_iter:\n",
    "        v_new, d_new = mcm.update(v, d)\n",
    "        error_1 = np.max(np.abs(v_new - v))\n",
    "        error_2 = np.abs(d_new - d)\n",
    "        error = max(error_1, error_2)\n",
    "        v = v_new\n",
    "        d = d_new\n",
    "        i += 1\n",
    "\n",
    "    return v, d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here’s a function `compute_reservation_wage` that takes an instance of `McCallModelContinuous`\n",
    "and returns the associated reservation wage.\n",
    "\n",
    "If $ v(w) < h $ for all $ w $, then the function returns np.inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def compute_reservation_wage(mcm):\n",
    "    \"\"\"\n",
    "    Computes the reservation wage of an instance of the McCall model\n",
    "    by finding the smallest w such that v(w) >= h.\n",
    "\n",
    "    If no such w exists, then w_bar is set to np.inf.\n",
    "    \"\"\"\n",
    "    u = lambda x: np.log(x)\n",
    "\n",
    "    v, d = solve_model(mcm)\n",
    "    h = u(mcm.c) + mcm.β * d\n",
    "\n",
    "    w_bar = np.inf\n",
    "    for i, wage in enumerate(mcm.w_grid):\n",
    "        if v[i] > h:\n",
    "            w_bar = wage\n",
    "            break\n",
    "\n",
    "    return w_bar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exercises ask you to explore the solution and how it changes with parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Use the code above to explore what happens to the reservation wage when the wage parameter $ \\mu $\n",
    "changes.\n",
    "\n",
    "Use the default parameters and $ \\mu $ in `mu_vals = np.linspace(0.0, 2.0, 15)`\n",
    "\n",
    "Is the impact on the reservation wage as you expected?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "Let us now consider how the agent responds to an increase in volatility.\n",
    "\n",
    "To try to understand this, compute the reservation wage when the wage offer\n",
    "distribution is uniform on $ (m - s, m + s) $ and $ s $ varies.\n",
    "\n",
    "The idea here is that we are holding the mean constant and spreading the\n",
    "support.\n",
    "\n",
    "(This is a form of *mean-preserving spread*.)\n",
    "\n",
    "Use `s_vals = np.linspace(1.0, 2.0, 15)` and `m = 2.0`.\n",
    "\n",
    "State how you expect the reservation wage vary with $ s $.\n",
    "\n",
    "Now compute it.  Is this as you expected?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Here is one solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mcm = McCallModelContinuous()\n",
    "mu_vals = np.linspace(0.0, 2.0, 15)\n",
    "w_bar_vals = np.empty_like(mu_vals)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "for i, m in enumerate(mu_vals):\n",
    "    mcm.w_draws = lognormal_draws(μ=m)\n",
    "    w_bar = compute_reservation_wage(mcm)\n",
    "    w_bar_vals[i] = w_bar\n",
    "\n",
    "ax.set(xlabel='mean', ylabel='reservation wage')\n",
    "ax.plot(mu_vals, w_bar_vals, label=r'$\\bar w$ as a function of $\\mu$')\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not surprisingly, the agent is more inclined to wait when the distribution of\n",
    "offers shifts to the right."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "Here is one solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mcm = McCallModelContinuous()\n",
    "s_vals = np.linspace(1.0, 2.0, 15)\n",
    "m = 2.0\n",
    "w_bar_vals = np.empty_like(s_vals)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "for i, s in enumerate(s_vals):\n",
    "    a, b = m - s, m + s\n",
    "    mcm.w_draws = np.random.uniform(low=a, high=b, size=10_000)\n",
    "    w_bar = compute_reservation_wage(mcm)\n",
    "    w_bar_vals[i] = w_bar\n",
    "\n",
    "ax.set(xlabel='volatility', ylabel='reservation wage')\n",
    "ax.plot(s_vals, w_bar_vals, label=r'$\\bar w$ as a function of wage volatility')\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reservation wage increases with volatility.\n",
    "\n",
    "One might think that higher volatility would make the agent more inclined to\n",
    "take a given offer, since doing so represents certainty and waiting represents\n",
    "risk.\n",
    "\n",
    "But job search is like holding an option: the worker is only exposed to upside risk (since, in a free market, no one can force them to take a bad offer).\n",
    "\n",
    "More volatility means higher upside potential, which encourages the agent to wait."
   ]
  }
 ],
 "metadata": {
  "date": 1584334746.2663426,
  "filename": "mccall_fitted_vfi.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Job Search III: Fitted Value Function Iteration"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}