lake_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='lake-model'></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": [
    "# A Lake Model of Employment and Unemployment\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [A Lake Model of Employment and Unemployment](#A-Lake-Model-of-Employment-and-Unemployment)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [The Model](#The-Model)  \n",
    "  - [Implementation](#Implementation)  \n",
    "  - [Dynamics of an Individual Worker](#Dynamics-of-an-Individual-Worker)  \n",
    "  - [Endogenous Job Finding Rate](#Endogenous-Job-Finding-Rate)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  \n",
    "  - [Lake Model Solutions](#Lake-Model-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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture describes what has come to be called a *lake model*.\n",
    "\n",
    "The lake model is a basic tool for modeling unemployment.\n",
    "\n",
    "It allows us to analyze\n",
    "\n",
    "- flows between unemployment and employment.  \n",
    "- how these flows influence steady state employment and unemployment rates.  \n",
    "\n",
    "\n",
    "It is a good model for interpreting monthly labor department reports on gross and net jobs created and jobs destroyed.\n",
    "\n",
    "The “lakes” in the model are the pools of employed and unemployed.\n",
    "\n",
    "The “flows” between the lakes are caused by\n",
    "\n",
    "- firing and hiring  \n",
    "- entry and exit from the labor force  \n",
    "\n",
    "\n",
    "For the first part of this lecture, the parameters governing transitions into\n",
    "and out of unemployment and employment are exogenous.\n",
    "\n",
    "Later, we’ll determine some of these transition rates endogenously using the [McCall search model](https://python.quantecon.org/mccall_model.html).\n",
    "\n",
    "We’ll also use some nifty concepts like ergodicity, which provides a fundamental link between *cross-sectional* and *long run time series* distributions.\n",
    "\n",
    "These concepts will help us build an equilibrium model of ex-ante homogeneous workers whose different luck generates variations in their ex post experiences.\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from quantecon import MarkovChain\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy.stats import norm\n",
    "from scipy.optimize import brentq\n",
    "from quantecon.distributions import BetaBinomial\n",
    "from numba import jit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prerequisites\n",
    "\n",
    "Before working through what follows, we recommend you read the [lecture\n",
    "on finite Markov chains](https://python.quantecon.org/finite_markov.html).\n",
    "\n",
    "You will also need some basic [linear algebra](https://python.quantecon.org/linear_algebra.html) and probability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Model\n",
    "\n",
    "The economy is inhabited by a very large number of ex-ante identical workers.\n",
    "\n",
    "The workers live forever, spending their lives moving between unemployment and employment.\n",
    "\n",
    "Their rates of  transition between employment and unemployment are  governed by the following parameters:\n",
    "\n",
    "- $ \\lambda $, the job finding rate for currently unemployed workers  \n",
    "- $ \\alpha $, the dismissal rate for currently employed workers  \n",
    "- $ b $, the entry rate into the labor force  \n",
    "- $ d $, the exit rate from the labor force  \n",
    "\n",
    "\n",
    "The growth rate of the labor force evidently equals $ g=b-d $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregate Variables\n",
    "\n",
    "We want to derive the dynamics of the following aggregates\n",
    "\n",
    "- $ E_t $, the total number of employed workers at date $ t $  \n",
    "- $ U_t $, the total number of unemployed workers at $ t $  \n",
    "- $ N_t $, the number of workers in the labor force at $ t $  \n",
    "\n",
    "\n",
    "We also want to know the values of the following objects\n",
    "\n",
    "- The employment rate $ e_t := E_t/N_t $.  \n",
    "- The unemployment rate $ u_t := U_t/N_t $.  \n",
    "\n",
    "\n",
    "(Here and below, capital letters represent stocks and lowercase letters represent flows)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laws of Motion for Stock Variables\n",
    "\n",
    "We begin by constructing laws of motion for the aggregate variables $ E_t,U_t, N_t $.\n",
    "\n",
    "Of the mass of workers $ E_t $ who are employed at date $ t $,\n",
    "\n",
    "- $ (1-d)E_t $ will remain in the labor force  \n",
    "- of these, $ (1-\\alpha)(1-d)E_t $ will remain employed  \n",
    "\n",
    "\n",
    "Of the mass of workers $ U_t $ workers who are currently unemployed,\n",
    "\n",
    "- $ (1-d)U_t $ will remain in the labor force  \n",
    "- of these, $ (1-d) \\lambda U_t $ will become employed  \n",
    "\n",
    "\n",
    "Therefore,  the number of workers who will be employed at date $ t+1 $ will be\n",
    "\n",
    "$$\n",
    "E_{t+1} = (1-d)(1-\\alpha)E_t + (1-d)\\lambda U_t\n",
    "$$\n",
    "\n",
    "A similar analysis implies\n",
    "\n",
    "$$\n",
    "U_{t+1} = (1-d)\\alpha E_t + (1-d)(1-\\lambda)U_t + b (E_t+U_t)\n",
    "$$\n",
    "\n",
    "The value $ b(E_t+U_t) $ is the mass of new workers entering the labor force unemployed.\n",
    "\n",
    "The total stock of workers $ N_t=E_t+U_t $ evolves as\n",
    "\n",
    "$$\n",
    "N_{t+1} = (1+b-d)N_t = (1+g)N_t\n",
    "$$\n",
    "\n",
    "Letting $ X_t := \\left(\\begin{matrix}U_t\\\\E_t\\end{matrix}\\right) $, the law of motion for $ X $  is\n",
    "\n",
    "$$\n",
    "X_{t+1} = A X_t\n",
    "\\quad \\text{where} \\quad\n",
    "A :=\n",
    "\\begin{pmatrix}\n",
    "    (1-d)(1-\\lambda) + b & (1-d)\\alpha + b  \\\\\n",
    "    (1-d)\\lambda & (1-d)(1-\\alpha)\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "This law tells us how total employment and unemployment evolve over time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laws of Motion for Rates\n",
    "\n",
    "Now let’s derive the law of motion for rates.\n",
    "\n",
    "To get these we can divide both sides of $ X_{t+1} = A X_t $ by  $ N_{t+1} $ to get\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "    U_{t+1}/N_{t+1} \\\\\n",
    "    E_{t+1}/N_{t+1}\n",
    "\\end{pmatrix} =\n",
    "\\frac1{1+g} A\n",
    "\\begin{pmatrix}\n",
    "    U_{t}/N_{t}\n",
    "    \\\\\n",
    "    E_{t}/N_{t}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Letting\n",
    "\n",
    "$$\n",
    "x_t :=\n",
    "\\left(\\begin{matrix}\n",
    "    u_t\\\\ e_t\n",
    "\\end{matrix}\\right) =\n",
    "\\left(\\begin{matrix}\n",
    "    U_t/N_t\\\\ E_t/N_t\n",
    "\\end{matrix}\\right)\n",
    "$$\n",
    "\n",
    "we can also write this as\n",
    "\n",
    "$$\n",
    "x_{t+1} = \\hat A x_t\n",
    "\\quad \\text{where} \\quad\n",
    "\\hat A := \\frac{1}{1 + g} A\n",
    "$$\n",
    "\n",
    "You can check that $ e_t + u_t = 1 $ implies that $ e_{t+1}+u_{t+1} = 1 $.\n",
    "\n",
    "This follows from the fact that the columns of $ \\hat A $ sum to 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "Let’s code up these equations.\n",
    "\n",
    "To do this we’re going to use a class that we’ll call LakeModel.\n",
    "\n",
    "This class will\n",
    "\n",
    "1. store the primitives $ \\alpha, \\lambda, b, d $  \n",
    "1. compute and store the implied objects $ g, A, \\hat A $  \n",
    "1. provide methods to simulate dynamics of the stocks and rates  \n",
    "1. provide a method to compute the state state of the rate  \n",
    "\n",
    "\n",
    "To write a nice implementation, there’s an issue we have to address.\n",
    "\n",
    "Derived data such as $ A $ depend on the primitives like $ \\alpha $\n",
    "and $ \\lambda $.\n",
    "\n",
    "If a user alters these primitives, we would ideally like derived data to\n",
    "update automatically.\n",
    "\n",
    "(For example, if a user changes the value of $ b $ for a given instance of the class, we would like $ g = b - d $ to update automatically)\n",
    "\n",
    "To achieve this outcome, we’re going to use descriptors and decorators such as @property.\n",
    "\n",
    "If you need to refresh your understanding of how these work, consult [this lecture](https://python.quantecon.org/python_advanced_features.html).\n",
    "\n",
    "Here’s the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "class LakeModel:\n",
    "    \"\"\"\n",
    "    Solves the lake model and computes dynamics of unemployment stocks and\n",
    "    rates.\n",
    "\n",
    "    Parameters:\n",
    "    ------------\n",
    "    λ : scalar\n",
    "        The job finding rate for currently unemployed workers\n",
    "    α : scalar\n",
    "        The dismissal rate for currently employed workers\n",
    "    b : scalar\n",
    "        Entry rate into the labor force\n",
    "    d : scalar\n",
    "        Exit rate from the labor force\n",
    "\n",
    "    \"\"\"\n",
    "    def __init__(self, λ=0.283, α=0.013, b=0.0124, d=0.00822):\n",
    "        self._λ, self._α, self._b, self._d = λ, α, b, d\n",
    "        self.compute_derived_values()\n",
    "\n",
    "    def compute_derived_values(self):\n",
    "        # Unpack names to simplify expression\n",
    "        λ, α, b, d = self._λ, self._α, self._b, self._d\n",
    "\n",
    "        self._g = b - d\n",
    "        self._A = np.array([[(1-d) * (1-λ) + b,      (1 - d) * α + b],\n",
    "                            [        (1-d) * λ,   (1 - d) * (1 - α)]])\n",
    "\n",
    "        self._A_hat = self._A / (1 + self._g)\n",
    "\n",
    "    @property\n",
    "    def g(self):\n",
    "        return self._g\n",
    "\n",
    "    @property\n",
    "    def A(self):\n",
    "        return self._A\n",
    "\n",
    "    @property\n",
    "    def A_hat(self):\n",
    "        return self._A_hat\n",
    "\n",
    "    @property\n",
    "    def λ(self):\n",
    "        return self._λ\n",
    "\n",
    "    @λ.setter\n",
    "    def λ(self, new_value):\n",
    "        self._α = new_value\n",
    "        self.compute_derived_values()\n",
    "\n",
    "    @property\n",
    "    def α(self):\n",
    "        return self._α\n",
    "\n",
    "    @α.setter\n",
    "    def α(self, new_value):\n",
    "        self._α = new_value\n",
    "        self.compute_derived_values()\n",
    "\n",
    "    @property\n",
    "    def b(self):\n",
    "        return self._b\n",
    "\n",
    "    @b.setter\n",
    "    def b(self, new_value):\n",
    "        self._b = new_value\n",
    "        self.compute_derived_values()\n",
    "\n",
    "    @property\n",
    "    def d(self):\n",
    "        return self._d\n",
    "\n",
    "    @d.setter\n",
    "    def d(self, new_value):\n",
    "        self._d = new_value\n",
    "        self.compute_derived_values()\n",
    "\n",
    "\n",
    "    def rate_steady_state(self, tol=1e-6):\n",
    "        \"\"\"\n",
    "        Finds the steady state of the system :math:`x_{t+1} = \\hat A x_{t}`\n",
    "\n",
    "        Returns\n",
    "        --------\n",
    "        xbar : steady state vector of employment and unemployment rates\n",
    "        \"\"\"\n",
    "        x = 0.5 * np.ones(2)\n",
    "        error = tol + 1\n",
    "        while error > tol:\n",
    "            new_x = self.A_hat @ x\n",
    "            error = np.max(np.abs(new_x - x))\n",
    "            x = new_x\n",
    "        return x\n",
    "\n",
    "    def simulate_stock_path(self, X0, T):\n",
    "        \"\"\"\n",
    "        Simulates the sequence of Employment and Unemployment stocks\n",
    "\n",
    "        Parameters\n",
    "        ------------\n",
    "        X0 : array\n",
    "            Contains initial values (E0, U0)\n",
    "        T : int\n",
    "            Number of periods to simulate\n",
    "\n",
    "        Returns\n",
    "        ---------\n",
    "        X : iterator\n",
    "            Contains sequence of employment and unemployment stocks\n",
    "        \"\"\"\n",
    "\n",
    "        X = np.atleast_1d(X0)  # Recast as array just in case\n",
    "        for t in range(T):\n",
    "            yield X\n",
    "            X = self.A @ X\n",
    "\n",
    "    def simulate_rate_path(self, x0, T):\n",
    "        \"\"\"\n",
    "        Simulates the sequence of employment and unemployment rates\n",
    "\n",
    "        Parameters\n",
    "        ------------\n",
    "        x0 : array\n",
    "            Contains initial values (e0,u0)\n",
    "        T : int\n",
    "            Number of periods to simulate\n",
    "\n",
    "        Returns\n",
    "        ---------\n",
    "        x : iterator\n",
    "            Contains sequence of employment and unemployment rates\n",
    "\n",
    "        \"\"\"\n",
    "        x = np.atleast_1d(x0)  # Recast as array just in case\n",
    "        for t in range(T):\n",
    "            yield x\n",
    "            x = self.A_hat @ x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As desired, if we create an instance and update a primitive like\n",
    "$ \\alpha $, derived objects like $ A $ will also change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm = LakeModel()\n",
    "lm.α"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm.A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm.α = 2\n",
    "lm.A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregate Dynamics\n",
    "\n",
    "Let’s run a simulation under the default parameters (see above) starting from $ X_0 = (12, 138) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm = LakeModel()\n",
    "N_0 = 150      # Population\n",
    "e_0 = 0.92     # Initial employment rate\n",
    "u_0 = 1 - e_0  # Initial unemployment rate\n",
    "T = 50         # Simulation length\n",
    "\n",
    "U_0 = u_0 * N_0\n",
    "E_0 = e_0 * N_0\n",
    "\n",
    "fig, axes = plt.subplots(3, 1, figsize=(10, 8))\n",
    "X_0 = (U_0, E_0)\n",
    "X_path = np.vstack(tuple(lm.simulate_stock_path(X_0, T)))\n",
    "\n",
    "axes[0].plot(X_path[:, 0], lw=2)\n",
    "axes[0].set_title('Unemployment')\n",
    "\n",
    "axes[1].plot(X_path[:, 1], lw=2)\n",
    "axes[1].set_title('Employment')\n",
    "\n",
    "axes[2].plot(X_path.sum(1), lw=2)\n",
    "axes[2].set_title('Labor force')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The aggregates $ E_t $ and $ U_t $ don’t converge because  their sum $ E_t + U_t $ grows at rate $ g $.\n",
    "\n",
    "On the other hand, the vector of employment and unemployment rates $ x_t $ can be in a steady state $ \\bar x $ if\n",
    "there exists an $ \\bar x $  such that\n",
    "\n",
    "- $ \\bar x = \\hat A \\bar x $  \n",
    "- the components satisfy $ \\bar e + \\bar u = 1 $  \n",
    "\n",
    "\n",
    "This equation tells us that a steady state level $ \\bar x $ is an  eigenvector of $ \\hat A $ associated with a unit eigenvalue.\n",
    "\n",
    "We also have $ x_t \\to \\bar x $ as $ t \\to \\infty $ provided that the remaining eigenvalue of $ \\hat A $ has modulus less that 1.\n",
    "\n",
    "This is the case for our default parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm = LakeModel()\n",
    "e, f = np.linalg.eigvals(lm.A_hat)\n",
    "abs(e), abs(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s look at the convergence of the unemployment and employment rate to steady state levels (dashed red line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm = LakeModel()\n",
    "e_0 = 0.92     # Initial employment rate\n",
    "u_0 = 1 - e_0  # Initial unemployment rate\n",
    "T = 50         # Simulation length\n",
    "\n",
    "xbar = lm.rate_steady_state()\n",
    "\n",
    "fig, axes = plt.subplots(2, 1, figsize=(10, 8))\n",
    "x_0 = (u_0, e_0)\n",
    "x_path = np.vstack(tuple(lm.simulate_rate_path(x_0, T)))\n",
    "\n",
    "titles = ['Unemployment rate', 'Employment rate']\n",
    "\n",
    "for i, title in enumerate(titles):\n",
    "    axes[i].plot(x_path[:, i], lw=2, alpha=0.5)\n",
    "    axes[i].hlines(xbar[i], 0, T, 'r', '--')\n",
    "    axes[i].set_title(title)\n",
    "    axes[i].grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dynamics of an Individual Worker\n",
    "\n",
    "An individual worker’s employment dynamics are governed by a [finite state Markov process](https://python.quantecon.org/finite_markov.html).\n",
    "\n",
    "The worker can be in one of two states:\n",
    "\n",
    "- $ s_t=0 $ means unemployed  \n",
    "- $ s_t=1 $ means employed  \n",
    "\n",
    "\n",
    "Let’s start off under the assumption that $ b = d = 0 $.\n",
    "\n",
    "The associated transition matrix is then\n",
    "\n",
    "$$\n",
    "P = \\left(\n",
    "        \\begin{matrix}\n",
    "            1 - \\lambda & \\lambda \\\\\n",
    "            \\alpha & 1 - \\alpha\n",
    "        \\end{matrix}\n",
    "    \\right)\n",
    "$$\n",
    "\n",
    "Let $ \\psi_t $ denote the [marginal distribution](https://python.quantecon.org/finite_markov.html#mc-md) over employment/unemployment states for the worker at time $ t $.\n",
    "\n",
    "As usual, we regard it as a row vector.\n",
    "\n",
    "We know [from an earlier discussion](https://python.quantecon.org/finite_markov.html#mc-md) that $ \\psi_t $ follows the law of motion\n",
    "\n",
    "$$\n",
    "\\psi_{t+1} = \\psi_t P\n",
    "$$\n",
    "\n",
    "We also know from the [lecture on finite Markov chains](https://python.quantecon.org/finite_markov.html)\n",
    "that if $ \\alpha \\in (0, 1) $ and $ \\lambda \\in (0, 1) $, then\n",
    "$ P $ has a unique stationary distribution, denoted here by $ \\psi^* $.\n",
    "\n",
    "The unique stationary distribution satisfies\n",
    "\n",
    "$$\n",
    "\\psi^*[0] = \\frac{\\alpha}{\\alpha + \\lambda}\n",
    "$$\n",
    "\n",
    "Not surprisingly, probability mass on the unemployment state increases with\n",
    "the dismissal rate and falls with the job finding rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ergodicity\n",
    "\n",
    "Let’s look at a typical lifetime of employment-unemployment spells.\n",
    "\n",
    "We want to compute the average amounts of time an infinitely lived worker would spend employed and unemployed.\n",
    "\n",
    "Let\n",
    "\n",
    "$$\n",
    "\\bar s_{u,T} := \\frac1{T} \\sum_{t=1}^T \\mathbb 1\\{s_t = 0\\}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\bar s_{e,T} := \\frac1{T} \\sum_{t=1}^T \\mathbb 1\\{s_t = 1\\}\n",
    "$$\n",
    "\n",
    "(As usual, $ \\mathbb 1\\{Q\\} = 1 $ if statement $ Q $ is true and 0 otherwise)\n",
    "\n",
    "These are the fraction of time a worker spends unemployed and employed, respectively, up until period $ T $.\n",
    "\n",
    "If $ \\alpha \\in (0, 1) $ and $ \\lambda \\in (0, 1) $, then $ P $ is [ergodic](https://python.quantecon.org/finite_markov.html#ergodicity), and hence we have\n",
    "\n",
    "$$\n",
    "\\lim_{T \\to \\infty} \\bar s_{u, T} = \\psi^*[0]\n",
    "\\quad \\text{and} \\quad\n",
    "\\lim_{T \\to \\infty} \\bar s_{e, T} = \\psi^*[1]\n",
    "$$\n",
    "\n",
    "with probability one.\n",
    "\n",
    "Inspection tells us that $ P $ is exactly the transpose of $ \\hat A $ under the assumption $ b=d=0 $.\n",
    "\n",
    "Thus, the percentages of time that an  infinitely lived worker  spends employed and unemployed equal the fractions of workers employed and unemployed in the steady state distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergence Rate\n",
    "\n",
    "How long does it take for time series sample averages to converge to cross-sectional averages?\n",
    "\n",
    "We can use [QuantEcon.py’s](http://quantecon.org/quantecon-py)\n",
    "MarkovChain class to investigate this.\n",
    "\n",
    "Let’s plot the path of the sample averages over 5,000 periods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm = LakeModel(d=0, b=0)\n",
    "T = 5000  # Simulation length\n",
    "\n",
    "α, λ = lm.α, lm.λ\n",
    "\n",
    "P = [[1 - λ,        λ],\n",
    "    [    α,    1 - α]]\n",
    "\n",
    "mc = MarkovChain(P)\n",
    "\n",
    "xbar = lm.rate_steady_state()\n",
    "\n",
    "fig, axes = plt.subplots(2, 1, figsize=(10, 8))\n",
    "s_path = mc.simulate(T, init=1)\n",
    "s_bar_e = s_path.cumsum() / range(1, T+1)\n",
    "s_bar_u = 1 - s_bar_e\n",
    "\n",
    "to_plot = [s_bar_u, s_bar_e]\n",
    "titles = ['Percent of time unemployed', 'Percent of time employed']\n",
    "\n",
    "for i, plot in enumerate(to_plot):\n",
    "    axes[i].plot(plot, lw=2, alpha=0.5)\n",
    "    axes[i].hlines(xbar[i], 0, T, 'r', '--')\n",
    "    axes[i].set_title(titles[i])\n",
    "    axes[i].grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stationary probabilities are given by the dashed red line.\n",
    "\n",
    "In this case it takes much of the sample for these two objects to converge.\n",
    "\n",
    "This is largely due to the high persistence in the Markov chain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Endogenous Job Finding Rate\n",
    "\n",
    "We now make the hiring rate endogenous.\n",
    "\n",
    "The transition rate from unemployment to employment will be determined by the McCall search model [[McC70]](https://python.quantecon.org/zreferences.html#mccall1970).\n",
    "\n",
    "All details relevant to the following discussion can be found in [our treatment](https://python.quantecon.org/mccall_model.html) of that model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reservation Wage\n",
    "\n",
    "The most important thing to remember about the model is that optimal decisions\n",
    "are characterized by a reservation wage $ \\bar w $\n",
    "\n",
    "- If the wage offer $ w $ in hand is greater than or equal to $ \\bar w $, then the worker accepts.  \n",
    "- Otherwise, the worker rejects.  \n",
    "\n",
    "\n",
    "As we saw in [our discussion of the model](https://python.quantecon.org/mccall_model.html), the reservation wage depends on the wage offer distribution and the parameters\n",
    "\n",
    "- $ \\alpha $, the separation rate  \n",
    "- $ \\beta $, the discount factor  \n",
    "- $ \\gamma $, the offer arrival rate  \n",
    "- $ c $, unemployment compensation  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linking the McCall Search Model to the Lake Model\n",
    "\n",
    "Suppose that  all workers inside a lake model behave according to the McCall search model.\n",
    "\n",
    "The exogenous probability of leaving employment remains $ \\alpha $.\n",
    "\n",
    "But their optimal decision rules determine the probability $ \\lambda $ of leaving unemployment.\n",
    "\n",
    "This is now\n",
    "\n",
    "\n",
    "<a id='equation-lake-lamda'></a>\n",
    "$$\n",
    "\\lambda\n",
    "= \\gamma \\mathbb P \\{ w_t \\geq \\bar w\\}\n",
    "= \\gamma \\sum_{w' \\geq \\bar w} p(w') \\tag{1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fiscal Policy\n",
    "\n",
    "We can use the McCall search version of the Lake Model  to find an optimal level of unemployment insurance.\n",
    "\n",
    "We assume that  the government sets unemployment compensation $ c $.\n",
    "\n",
    "The government imposes a lump-sum tax $ \\tau $ sufficient to finance total unemployment payments.\n",
    "\n",
    "To attain a balanced budget at a steady state, taxes, the steady state unemployment rate $ u $, and the unemployment compensation rate must satisfy\n",
    "\n",
    "$$\n",
    "\\tau = u c\n",
    "$$\n",
    "\n",
    "The lump-sum tax applies to everyone, including unemployed workers.\n",
    "\n",
    "Thus, the post-tax income of an employed worker with wage $ w $ is $ w - \\tau $.\n",
    "\n",
    "The post-tax income of an unemployed worker is $ c - \\tau $.\n",
    "\n",
    "For each specification $ (c, \\tau) $ of government policy, we can solve for the worker’s optimal reservation wage.\n",
    "\n",
    "This determines $ \\lambda $ via [(1)](#equation-lake-lamda) evaluated at post tax wages, which in turn determines a steady state unemployment rate $ u(c, \\tau) $.\n",
    "\n",
    "For a given level of unemployment benefit $ c $, we can solve for a tax that balances the budget in the steady state\n",
    "\n",
    "$$\n",
    "\\tau = u(c, \\tau) c\n",
    "$$\n",
    "\n",
    "To evaluate alternative government tax-unemployment compensation pairs, we require a welfare criterion.\n",
    "\n",
    "We use a steady state welfare criterion\n",
    "\n",
    "$$\n",
    "W := e \\,  {\\mathbb E} [V \\, | \\,  \\text{employed}] + u \\,  U\n",
    "$$\n",
    "\n",
    "where the notation $ V $ and $ U $ is as defined in the [McCall search model lecture](https://python.quantecon.org/mccall_model.html).\n",
    "\n",
    "The wage offer distribution will be a discretized version of the lognormal distribution $ LN(\\log(20),1) $, as shown in the next figure\n",
    "\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/python.quantecon.org/_static/lecture_specific/lake_model/lake_distribution_wages.png\" style=\"\">\n",
    "\n",
    "  \n",
    "We take a period to be a month.\n",
    "\n",
    "We set $ b $ and $ d $ to match monthly [birth](http://www.cdc.gov/nchs/fastats/births.htm) and [death rates](http://www.cdc.gov/nchs/fastats/deaths.htm), respectively, in the U.S. population\n",
    "\n",
    "- $ b = 0.0124 $  \n",
    "- $ d = 0.00822 $  \n",
    "\n",
    "\n",
    "Following [[DFH06]](https://python.quantecon.org/zreferences.html#davis2006flow), we set $ \\alpha $, the hazard rate of leaving employment, to\n",
    "\n",
    "- $ \\alpha = 0.013 $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fiscal Policy Code\n",
    "\n",
    "We will make use of techniques from the [McCall model lecture](https://python.quantecon.org/mccall_model.html)\n",
    "\n",
    "The first piece of code implements value function iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "html-class": "collapse"
   },
   "outputs": [],
   "source": [
    "# A default utility function\n",
    "\n",
    "@jit\n",
    "def u(c, σ):\n",
    "    if c > 0:\n",
    "        return (c**(1 - σ) - 1) / (1 - σ)\n",
    "    else:\n",
    "        return -10e6\n",
    "\n",
    "\n",
    "class McCallModel:\n",
    "    \"\"\"\n",
    "    Stores the parameters and functions associated with a given model.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self,\n",
    "                α=0.2,       # Job separation rate\n",
    "                β=0.98,      # Discount rate\n",
    "                γ=0.7,       # Job offer rate\n",
    "                c=6.0,       # Unemployment compensation\n",
    "                σ=2.0,       # Utility parameter\n",
    "                w_vec=None,  # Possible wage values\n",
    "                p_vec=None): # Probabilities over w_vec\n",
    "\n",
    "        self.α, self.β, self.γ, self.c = α, β, γ, c\n",
    "        self.σ = σ\n",
    "\n",
    "        # Add a default wage vector and probabilities over the vector using\n",
    "        # the beta-binomial distribution\n",
    "        if w_vec is None:\n",
    "            n = 60  # Number of possible outcomes for wage\n",
    "            # Wages between 10 and 20\n",
    "            self.w_vec = np.linspace(10, 20, n)\n",
    "            a, b = 600, 400  # Shape parameters\n",
    "            dist = BetaBinomial(n-1, a, b)\n",
    "            self.p_vec = dist.pdf()\n",
    "        else:\n",
    "            self.w_vec = w_vec\n",
    "            self.p_vec = p_vec\n",
    "\n",
    "@jit\n",
    "def _update_bellman(α, β, γ, c, σ, w_vec, p_vec, V, V_new, U):\n",
    "    \"\"\"\n",
    "    A jitted function to update the Bellman equations.  Note that V_new is\n",
    "    modified in place (i.e, modified by this function).  The new value of U\n",
    "    is returned.\n",
    "\n",
    "    \"\"\"\n",
    "    for w_idx, w in enumerate(w_vec):\n",
    "        # w_idx indexes the vector of possible wages\n",
    "        V_new[w_idx] = u(w, σ) + β * ((1 - α) * V[w_idx] + α * U)\n",
    "\n",
    "    U_new = u(c, σ) + β * (1 - γ) * U + \\\n",
    "                    β * γ * np.sum(np.maximum(U, V) * p_vec)\n",
    "\n",
    "    return U_new\n",
    "\n",
    "\n",
    "def solve_mccall_model(mcm, tol=1e-5, max_iter=2000):\n",
    "    \"\"\"\n",
    "    Iterates to convergence on the Bellman equations\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    mcm : an instance of McCallModel\n",
    "    tol : float\n",
    "        error tolerance\n",
    "    max_iter : int\n",
    "        the maximum number of iterations\n",
    "    \"\"\"\n",
    "\n",
    "    V = np.ones(len(mcm.w_vec))  # Initial guess of V\n",
    "    V_new = np.empty_like(V)     # To store updates to V\n",
    "    U = 1                        # Initial guess of U\n",
    "    i = 0\n",
    "    error = tol + 1\n",
    "\n",
    "    while error > tol and i < max_iter:\n",
    "        U_new = _update_bellman(mcm.α, mcm.β, mcm.γ,\n",
    "                mcm.c, mcm.σ, mcm.w_vec, mcm.p_vec, V, V_new, U)\n",
    "        error_1 = np.max(np.abs(V_new - V))\n",
    "        error_2 = np.abs(U_new - U)\n",
    "        error = max(error_1, error_2)\n",
    "        V[:] = V_new\n",
    "        U = U_new\n",
    "        i += 1\n",
    "\n",
    "    return V, U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second piece of code is used to complete the reservation wage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "html-class": "collapse"
   },
   "outputs": [],
   "source": [
    "def compute_reservation_wage(mcm, return_values=False):\n",
    "    \"\"\"\n",
    "    Computes the reservation wage of an instance of the McCall model\n",
    "    by finding the smallest w such that V(w) > U.\n",
    "\n",
    "    If V(w) > U for all w, then the reservation wage w_bar is set to\n",
    "    the lowest wage in mcm.w_vec.\n",
    "\n",
    "    If v(w) < U for all w, then w_bar is set to np.inf.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    mcm : an instance of McCallModel\n",
    "    return_values : bool (optional, default=False)\n",
    "        Return the value functions as well\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    w_bar : scalar\n",
    "        The reservation wage\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    V, U = solve_mccall_model(mcm)\n",
    "    w_idx = np.searchsorted(V - U, 0)\n",
    "\n",
    "    if w_idx == len(V):\n",
    "        w_bar = np.inf\n",
    "    else:\n",
    "        w_bar = mcm.w_vec[w_idx]\n",
    "\n",
    "    if return_values == False:\n",
    "        return w_bar\n",
    "    else:\n",
    "        return w_bar, V, U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let’s compute and plot welfare, employment, unemployment, and tax revenue as a\n",
    "function of the unemployment compensation rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Some global variables that will stay constant\n",
    "α = 0.013\n",
    "α_q = (1-(1-α)**3)   # Quarterly (α is monthly)\n",
    "b = 0.0124\n",
    "d = 0.00822\n",
    "β = 0.98\n",
    "γ = 1.0\n",
    "σ = 2.0\n",
    "\n",
    "# The default wage distribution --- a discretized lognormal\n",
    "log_wage_mean, wage_grid_size, max_wage = 20, 200, 170\n",
    "logw_dist = norm(np.log(log_wage_mean), 1)\n",
    "w_vec = np.linspace(1e-8, max_wage, wage_grid_size + 1)\n",
    "cdf = logw_dist.cdf(np.log(w_vec))\n",
    "pdf = cdf[1:] - cdf[:-1]\n",
    "p_vec = pdf / pdf.sum()\n",
    "w_vec = (w_vec[1:] + w_vec[:-1]) / 2\n",
    "\n",
    "\n",
    "def compute_optimal_quantities(c, τ):\n",
    "    \"\"\"\n",
    "    Compute the reservation wage, job finding rate and value functions\n",
    "    of the workers given c and τ.\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    mcm = McCallModel(α=α_q,\n",
    "                    β=β,\n",
    "                    γ=γ,\n",
    "                    c=c-τ,          # Post tax compensation\n",
    "                    σ=σ,\n",
    "                    w_vec=w_vec-τ,  # Post tax wages\n",
    "                    p_vec=p_vec)\n",
    "\n",
    "    w_bar, V, U = compute_reservation_wage(mcm, return_values=True)\n",
    "    λ = γ * np.sum(p_vec[w_vec - τ > w_bar])\n",
    "    return w_bar, λ, V, U\n",
    "\n",
    "def compute_steady_state_quantities(c, τ):\n",
    "    \"\"\"\n",
    "    Compute the steady state unemployment rate given c and τ using optimal\n",
    "    quantities from the McCall model and computing corresponding steady\n",
    "    state quantities\n",
    "\n",
    "    \"\"\"\n",
    "    w_bar, λ, V, U = compute_optimal_quantities(c, τ)\n",
    "\n",
    "    # Compute steady state employment and unemployment rates\n",
    "    lm = LakeModel(α=α_q, λ=λ, b=b, d=d)\n",
    "    x = lm.rate_steady_state()\n",
    "    u, e = x\n",
    "\n",
    "    # Compute steady state welfare\n",
    "    w = np.sum(V * p_vec * (w_vec - τ > w_bar)) / np.sum(p_vec * (w_vec -\n",
    "            τ > w_bar))\n",
    "    welfare = e * w + u * U\n",
    "\n",
    "    return e, u, welfare\n",
    "\n",
    "\n",
    "def find_balanced_budget_tax(c):\n",
    "    \"\"\"\n",
    "    Find the tax level that will induce a balanced budget.\n",
    "\n",
    "    \"\"\"\n",
    "    def steady_state_budget(t):\n",
    "        e, u, w = compute_steady_state_quantities(c, t)\n",
    "        return t - u * c\n",
    "\n",
    "    τ = brentq(steady_state_budget, 0.0, 0.9 * c)\n",
    "    return τ\n",
    "\n",
    "\n",
    "# Levels of unemployment insurance we wish to study\n",
    "c_vec = np.linspace(5, 140, 60)\n",
    "\n",
    "tax_vec = []\n",
    "unempl_vec = []\n",
    "empl_vec = []\n",
    "welfare_vec = []\n",
    "\n",
    "for c in c_vec:\n",
    "    t = find_balanced_budget_tax(c)\n",
    "    e_rate, u_rate, welfare = compute_steady_state_quantities(c, t)\n",
    "    tax_vec.append(t)\n",
    "    unempl_vec.append(u_rate)\n",
    "    empl_vec.append(e_rate)\n",
    "    welfare_vec.append(welfare)\n",
    "\n",
    "fig, axes = plt.subplots(2, 2, figsize=(12, 10))\n",
    "\n",
    "plots = [unempl_vec, empl_vec, tax_vec, welfare_vec]\n",
    "titles = ['Unemployment', 'Employment', 'Tax', 'Welfare']\n",
    "\n",
    "for ax, plot, title in zip(axes.flatten(), plots, titles):\n",
    "    ax.plot(c_vec, plot, lw=2, alpha=0.7)\n",
    "    ax.set_title(title)\n",
    "    ax.grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Welfare first increases and then decreases as unemployment benefits rise.\n",
    "\n",
    "The level that maximizes steady state welfare is approximately 62."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Consider an economy with an initial stock  of workers $ N_0 = 100 $ at the\n",
    "steady state level of employment in the baseline parameterization\n",
    "\n",
    "- $ \\alpha = 0.013 $  \n",
    "- $ \\lambda = 0.283 $  \n",
    "- $ b = 0.0124 $  \n",
    "- $ d = 0.00822 $  \n",
    "\n",
    "\n",
    "(The values for $ \\alpha $ and $ \\lambda $ follow [[DFH06]](https://python.quantecon.org/zreferences.html#davis2006flow))\n",
    "\n",
    "Suppose that in response to new legislation the hiring rate reduces to $ \\lambda = 0.2 $.\n",
    "\n",
    "Plot the transition dynamics of the unemployment and employment stocks for 50 periods.\n",
    "\n",
    "Plot the transition dynamics for the rates.\n",
    "\n",
    "How long does the economy take to converge to its new steady state?\n",
    "\n",
    "What is the new steady state level of employment?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "Consider an economy with an initial stock  of workers $ N_0 = 100 $ at the\n",
    "steady state level of employment in the baseline parameterization.\n",
    "\n",
    "Suppose that for 20 periods the birth rate was temporarily high ($ b = 0.0025 $) and then returned to its original level.\n",
    "\n",
    "Plot the transition dynamics of the unemployment and employment stocks for 50 periods.\n",
    "\n",
    "Plot the transition dynamics for the rates.\n",
    "\n",
    "How long does the economy take to return to its original steady state?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lake Model Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "We begin by constructing the class containing the default parameters and assigning the\n",
    "steady state values to x0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm = LakeModel()\n",
    "x0 = lm.rate_steady_state()\n",
    "print(f\"Initial Steady State: {x0}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the simulation values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N0 = 100\n",
    "T = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "New legislation changes $ \\lambda $ to $ 0.2 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm.lmda = 0.2\n",
    "\n",
    "xbar = lm.rate_steady_state()  # new steady state\n",
    "X_path = np.vstack(tuple(lm.simulate_stock_path(x0 * N0, T)))\n",
    "x_path = np.vstack(tuple(lm.simulate_rate_path(x0, T)))\n",
    "print(f\"New Steady State: {xbar}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot stocks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 1, figsize=[10, 9])\n",
    "\n",
    "axes[0].plot(X_path[:, 0])\n",
    "axes[0].set_title('Unemployment')\n",
    "\n",
    "axes[1].plot(X_path[:, 1])\n",
    "axes[1].set_title('Employment')\n",
    "\n",
    "axes[2].plot(X_path.sum(1))\n",
    "axes[2].set_title('Labor force')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And how the rates evolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=(10, 8))\n",
    "\n",
    "titles = ['Unemployment rate', 'Employment rate']\n",
    "\n",
    "for i, title in enumerate(titles):\n",
    "    axes[i].plot(x_path[:, i])\n",
    "    axes[i].hlines(xbar[i], 0, T, 'r', '--')\n",
    "    axes[i].set_title(title)\n",
    "    axes[i].grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that it takes 20 periods for the economy to converge to its new\n",
    "steady state levels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "This next exercise has the economy experiencing a boom in entrances to\n",
    "the labor market and then later returning to the original levels.\n",
    "\n",
    "For 20 periods the economy has a new entry rate into the labor market.\n",
    "\n",
    "Let’s start off at the baseline parameterization and record the steady\n",
    "state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm = LakeModel()\n",
    "x0 = lm.rate_steady_state()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the other parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "b_hat = 0.003\n",
    "T_hat = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s increase $ b $ to the new value and simulate for 20 periods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm.b = b_hat\n",
    "# Simulate stocks\n",
    "X_path1 = np.vstack(tuple(lm.simulate_stock_path(x0 * N0, T_hat)))\n",
    "# Simulate rates\n",
    "x_path1 = np.vstack(tuple(lm.simulate_rate_path(x0, T_hat)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we reset $ b $ to the original value and then, using the state\n",
    "after 20 periods for the new initial conditions, we simulate for the\n",
    "additional 30 periods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lm.b = 0.0124\n",
    "# Simulate stocks\n",
    "X_path2 = np.vstack(tuple(lm.simulate_stock_path(X_path1[-1, :2], T-T_hat+1)))\n",
    "# Simulate rates\n",
    "x_path2 = np.vstack(tuple(lm.simulate_rate_path(x_path1[-1, :2], T-T_hat+1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we combine these two paths and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# note [1:] to avoid doubling period 20\n",
    "x_path = np.vstack([x_path1, x_path2[1:]])\n",
    "X_path = np.vstack([X_path1, X_path2[1:]])\n",
    "\n",
    "fig, axes = plt.subplots(3, 1, figsize=[10, 9])\n",
    "\n",
    "axes[0].plot(X_path[:, 0])\n",
    "axes[0].set_title('Unemployment')\n",
    "\n",
    "axes[1].plot(X_path[:, 1])\n",
    "axes[1].set_title('Employment')\n",
    "\n",
    "axes[2].plot(X_path.sum(1))\n",
    "axes[2].set_title('Labor force')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=[10, 6])\n",
    "\n",
    "titles = ['Unemployment rate', 'Employment rate']\n",
    "\n",
    "for i, title in enumerate(titles):\n",
    "    axes[i].plot(x_path[:, i])\n",
    "    axes[i].hlines(x0[i], 0, T, 'r', '--')\n",
    "    axes[i].set_title(title)\n",
    "    axes[i].grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "date": 1584334742.2079854,
  "filename": "lake_model.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "A Lake Model of Employment and Unemployment"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}