28_zurcher.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Foundations of Computational Economics #28\n",
    "\n",
    "by Fedor Iskhakov, ANU\n",
    "\n",
    "<img src=\"_static/img/dag3logo.png\" style=\"width:256px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Rust model of bus engine replacement\n",
    "\n",
    "<img src=\"_static/img/lecture.png\" style=\"width:64px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"_static/img/youtube.png\" style=\"width:65px;\">\n",
    "\n",
    "[https://youtu.be/Yo69tlKKplY](https://youtu.be/Yo69tlKKplY)\n",
    "\n",
    "Description: Model background and formulation. Mileage process. Optimal replacement choice with and without EV taste shocks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Empirical Model of Harold Zurcher\n",
    "\n",
    "📖 **John Rust (1987, Econometrica) “Optimal Replacement of GMC Bus Engines: An Empirical Model of Harold Zurcher”**\n",
    "\n",
    "- Foundational paper in dynamic structural econometrics  \n",
    "- Develops the framework used extensively in many fields today  \n",
    "\n",
    "\n",
    "Ingredients:\n",
    "\n",
    "1. Simple dynamic model: binary discrete state regenerative optimal stopping problem  \n",
    "1. Smart solution method: polyalgorithm with fast Newton-Kantorovich iterations  \n",
    "1. Coherent econometrics specification: unobserved states ensuring not non-degenerate likelihood  \n",
    "1. New maximum likelihood estimator: nested fixed point (NFXP) estimator  \n",
    "1. Very real application with actual data collected by a real person named Harold Zurcher  \n",
    "\n",
    "\n",
    "**Still very powerful method for both solving and estimating dynamic discrete choice models**\n",
    "\n",
    "📖 Iskhakov, Lee, Rust, Schjerning and Seo (2016, Econometrica) “Comment on “Constrained Optimization Approaches to Estimation of Structural Models””"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Components of the dynamic model\n",
    "\n",
    "- **State variables** — vector of variables that describe all relevant\n",
    "  information about the modeled decision process, $ x_t $  \n",
    "- **Decision variables** — vector of variables describing the choices, $ d_t $  \n",
    "- **Instantaneous payoff** — utility function, $ u(x_t,d_t) $, with\n",
    "  time separable discounted utility  \n",
    "- **Motion rules** — agent’s beliefs of how state variable evolve\n",
    "  through time, conditional on choices, $ x_{t+1} \\sim F(x_t,d_t) $  \n",
    "\n",
    "\n",
    "Solution is given by:\n",
    "\n",
    "- **Value function** — maximum attainable utility $ V(x_t) $  \n",
    "- **Policy function** — mapping from state space to action space that\n",
    "  returns the optimal choice, $ d^{\\star}(x_t) $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Zurcher optimization problem\n",
    "\n",
    "- choice set: $ \\{ \\text{keep} , \\text{replace} \\} = \\{0,1\\} $  \n",
    "\n",
    "\n",
    "Each bus comes in for repair once a month and Zurcher chooses between ordinary maintenance $ (d_{t}=0) $ and overhaul/engine replacement $ (d_{t}=1) $.\n",
    "\n",
    "- state space: mileage at time t since last engine overhaul $ x_t $  \n",
    "\n",
    "\n",
    "Harold observes many different attributes of the buses which come into the shop, but we focus on the main one for now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Zurcher’s preferences\n",
    "\n",
    "Instantaneous payoffs are given by the cost function that depends on the choice\n",
    "\n",
    "$$\n",
    "u(x_{t},d_t,\\theta_1)=\\left \\{\n",
    "\\begin{array}{ll}\n",
    "    -RC-c(0,\\theta_1) & \\text{if }d_{t}=\\text{replace}=1 \\\\\n",
    "    -c(x_{t},\\theta_1) & \\text{if }d_{t}=\\text{keep}=0\n",
    "\\end{array} \\right.\n",
    "$$\n",
    "\n",
    "- $ RC $ = replacement cost  \n",
    "- $ c(x,\\theta_1) $ = cost of maintenance with preference parameters $ \\theta_1 $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Motion rules\n",
    "\n",
    "- First of all, mileage is continuous. How to deal with the *continuous* state space?  \n",
    "- Rust discretized the range of travelled miles into $ n=175 $ bins, indexed with $ i $, $ \\hat{X} = \\{\\hat{x}_1,...,\\hat{x}_n\\} $ with $ \\hat{x}_1=0 $  \n",
    "\n",
    "\n",
    "Mileage transition probability: for $ j = 0,...,J $\n",
    "\n",
    "$$\n",
    "p(x'|\\hat{x}_k, d,\\theta_2)=\n",
    "\\begin{cases}\n",
    "Pr\\{x'=\\hat{x}_{k+j}|\\theta_2\\}= \\theta_{2j} \\text{ if }  d=0  \\\\\n",
    "Pr\\{x'=\\hat{x}_{1+j}|\\theta_2\\}= \\theta_{2j}  \\text{ if }  d=1\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "- Mileage in the next period $ x' $ can move up at most $ J $ grid points  \n",
    "- $ J $ is determined from the observed distribution of mileage  \n",
    "- Effectively, the probabilities of increase of mileage from any existing mileage are the same  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Transition matrix for mileage\n",
    "\n",
    "If not replacing ($ d=0) $\n",
    "\n",
    "$$\n",
    "\\Pi(d=0)_{n x n} =\n",
    "\\begin{pmatrix}\n",
    "\\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & \\cdot & \\cdot & 0 \\\\\n",
    "0 & \\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & \\cdot & 0 \\\\\n",
    "0 & 0 &\\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & 0 \\\\\n",
    "\\cdot & \\cdot & \\cdot & \\cdot & \\cdot & \\cdot & \\cdot & \\cdot \\\\\n",
    "0 & \\cdot & \\cdot & 0 & \\theta_{20} & \\theta_{21} & \\theta_{22} & 0 \\\\\n",
    "0 & \\cdot & \\cdot & \\cdot & 0 & \\theta_{20} & \\theta_{21} & \\theta_{22} \\\\\n",
    "0 & \\cdot & \\cdot & \\cdot & \\cdot  & 0 & \\theta_{20} & 1-\\theta_{20} \\\\\n",
    "0 & \\cdot & \\cdot & \\cdot & \\cdot & \\cdot  & 0 & 1\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Transition matrix for mileage\n",
    "\n",
    "If replacing ($ d=1) $\n",
    "\n",
    "$$\n",
    "\\Pi(d=1)_{n x n} =\n",
    "\\begin{pmatrix}\n",
    "\\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & \\cdot & \\cdot & 0 \\\\\n",
    "\\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & \\cdot & \\cdot & 0 \\\\\n",
    "\\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & \\cdot & \\cdot & 0 \\\\\n",
    "\\cdot & \\cdot & \\cdot & \\cdot & \\cdot & \\cdot & \\cdot & \\cdot \\\\\n",
    "\\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & \\cdot & \\cdot & 0 \\\\\n",
    "\\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & \\cdot & \\cdot & 0 \\\\\n",
    "\\theta_{20} & \\theta_{21} & \\theta_{22} & 0 & \\cdot & \\cdot & \\cdot & 0 \\\\\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Dynamic optimization problem\n",
    "\n",
    "To minimize the discounted expected value of the costs, Zurcher should find such policy function $ f(x_t):X\\rightarrow \\{\\text{keep},\\text{replace}\\} $ such that $ d_t=f(x_t) $ maximizes\n",
    "\n",
    "$$\n",
    "\\mathbb{E}\\sum_{t=0}^{\\infty} \\beta^t u(x_t,d_t) \\longrightarrow \\max\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Value function\n",
    "\n",
    "The function $ V(x_t) $ denotes the maximum attainable value at period t\n",
    "\n",
    "$$\n",
    "V(x_t) = max_{\\Pi} \\mathbb{E} \\sum_{j=t}^{\\infty} \\beta^{j-t} u(x_t,d_t)\n",
    "$$\n",
    "\n",
    "where $ \\Pi $ is a space of policy functions $ f(x_t):X\\rightarrow \\{\\text{keep},\\text{replace}\\} $, and $ d_t = f(x_t) $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Recursive form of the maximization problem (Bellman equation)\n",
    "\n",
    "Using Bellman Principle of Optimality, we can show that the value function $ V(x_t) $ constitutes the solution of the following functional equation\n",
    "\n",
    "$$\n",
    "V(x) = \\max_{d\\in \\{\\text{keep},\\text{replace}\\}} \\big\\{ u(x,d) + \\beta \\mathbb{E}\\big[ V(x')\\big|x,d\\big] \\big\\}\n",
    "$$\n",
    "\n",
    "where expectation is taken over the next period values of state $ x' $ given the motion rule of the problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Bellman equation and Bellman operator\n",
    "\n",
    "$$\n",
    "V(x) = \\max_{d\\in \\{\\text{keep},\\text{replace}\\}} \\big\\{ u(x,d) + \\beta \\mathbb{E}\\big[ V(x')\\big|x,d\\big] \\big\\}\n",
    "$$\n",
    "\n",
    "can be written as a fixed point equation of the **Bellman operator** in the functional space\n",
    "\n",
    "$$\n",
    "T(V)(x) \\equiv \\max_{d \\in \\{\\text{keep},\\text{replace}\\}} \\big\\{ u(x,d) + \\beta \\mathbb{E}\\big[ V(x') \\big|x,d\\big] \\big\\}\n",
    "$$\n",
    "\n",
    "The Bellman equations is then $ V(x) = T(V)(x) $, with the\n",
    "solution given by the fixed point $ T(V) = V $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the **finite horizon** models there is last period $ T $\n",
    "\n",
    "- for $ t=T $ Bellman equation does not have the second (discounted expectation) term  \n",
    "- for $ t=T $ can typically be solved analytically  \n",
    "- value and policy functions are **time specific**, i.e. differ in time  \n",
    "- can be computed with **backward induction** with finite number of steps  \n",
    "- all the examples in previous video  \n",
    "\n",
    "\n",
    "In the **infinite horizon** models there is no last period, so $ T=\\infty $\n",
    "\n",
    "- value and policy functions are **time invariant**  \n",
    "- constitute the fixed point of the functional Bellman operator  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Contraction Mapping Theorem (Banach Fixed Point Theorem)\n",
    "\n",
    "Let $ (S,\\rho) $ be a complete metric space with a contraction mapping $ T: S \\rightarrow S $.\n",
    "\n",
    "Then\n",
    "1. $ T $ admits a unique fixed-point $ V^{\\star} \\in S $, i.e. $ T(V^{\\star}) = V^{\\star} $.\n",
    "2. $ V^{\\star} $ can be found by repeated application of the operator $ T $, i.e. $ T^n(V) \\rightarrow V^{\\star} $ as $ n\\rightarrow \\infty $.\n",
    "\n",
    "We will soon see that Bellman operator in a contraction mapping!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Value function iterations solution algorithm\n",
    "\n",
    "Algorithm:\n",
    "\n",
    "1. Initialize the value function $ V $ with arbitrary values  \n",
    "1. Solve Bellman equation to compute $ T(V) $  \n",
    "1. Repeat until convergence  \n",
    "\n",
    "\n",
    "- Unique fixed point $ \\Leftrightarrow $ unique solution to the Bellman equation  \n",
    "- The fixed point algorithm can start from **arbitrary initial guess**! VFI algorithm converges globally  \n",
    "- Direct resemblance with backward induction, only have to stop after unspecified number of iterations  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Making the model suitable for empirical work\n",
    "\n",
    "- Remember that Zurcher observes many different attributes of the busses that come into the shop  \n",
    "- But we as an econometrician do not!  \n",
    "- Yet, these are likely to be the reason for observing different behavior in same states (for the same observed mileage)  \n",
    "\n",
    "\n",
    "Need to include **error terms** into the model, denote $ \\epsilon $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Three independence assumptions\n",
    "\n",
    "1. Error terms are **independent across observations** due to random sampling  \n",
    "1. Error terms come in pairs, one for each decision $ d=0 $ and $ d=1 $, and are **independent across choices**  \n",
    "1. Conditional on $ x $, there is no serial correlation in error terms **across time**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Updating the Bellman equation\n",
    "\n",
    "$ \\varepsilon $ is a new (vector) state variable\n",
    "\n",
    "$$\n",
    "V(x,\\varepsilon) = \\max_{d\\in \\{0,1\\}} \\big\\{ u(x,\\varepsilon_d,d) + \\beta \\mathbb{E}\\big[ V(x',\\varepsilon')\\big|x,\\varepsilon,d\\big] \\big\\}\n",
    "$$\n",
    "\n",
    "$$\n",
    "V(x,\\varepsilon) = \\max_{d\\in \\{0,1\\}} \\big\\{ u(x,\\varepsilon_d,d) + \\beta\n",
    "\\int_{X} \\int_{\\Omega} V(x',\\varepsilon') p(x',\\varepsilon'|x,\\varepsilon,d) dx'd\\varepsilon' \\big\\}\n",
    "$$\n",
    "\n",
    "where $ \\varepsilon_d $ is the component of vector $ \\varepsilon \\in \\mathbb{R}^2 $ which corresponds to $ d $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Rust assumptions\n",
    "\n",
    "**(AS)** Additive separability in preferences\n",
    "\n",
    "$$\n",
    "u(x,\\varepsilon_d,d) = u(x,d) + \\varepsilon_d,\n",
    "$$\n",
    "\n",
    "**(CI)** Conditional independence\n",
    "\n",
    "$$\n",
    "p(x',\\varepsilon'|x,\\varepsilon,d) = q(\\varepsilon'|x')\\cdot \\pi(x'|x,d)\n",
    "$$\n",
    "\n",
    "**(EV)** Extreme value Type I (EV1) distribution of $ \\varepsilon $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### What Rust assumptions allow:\n",
    "\n",
    "$$\n",
    "V(x,\\varepsilon) = \\max_{d\\in \\{0,1\\}} \\big\\{ u(x,d) + \\varepsilon_d + \\beta\n",
    "\\int_{X} \\int_{\\Omega} V(x',\\varepsilon') \\pi(x'|x,d) q(\\varepsilon'|x') dx' d\\varepsilon' \\big\\}\n",
    "$$\n",
    "\n",
    "1. Separate out the deterministic part of **choice specific value function** $ v(x,d) $ (SA)  \n",
    "1. Compute the expectation by part (CI)  \n",
    "1. Use max-stability of EV1 to compute expectation w.r.t. $ \\varepsilon' $ (EV)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$$\n",
    "V(x,\\varepsilon) = \\max_{d\\in \\{0,1\\}} \\big\\{ \\underbrace{u(x,d) + \\beta\n",
    "\\int_{X} \\Big( \\int_{\\Omega} V(x',\\varepsilon') q(\\varepsilon'|x') d\\varepsilon'\\Big)\n",
    "\\pi(x'|x,d) dx'}_{v(x,d)}\n",
    "+ \\varepsilon_d \\big\\}\n",
    "$$\n",
    "\n",
    "$$\n",
    "V(x',\\varepsilon') = \\max_{d\\in \\{0,1\\}} \\big\\{ v(x',d) + \\varepsilon'_d \\big\\}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbb{E}\\big[ V(x',\\varepsilon')\\big|x,d\\big] =\n",
    "\\int_{X} \\log \\big( \\exp[v(x',0)] + \\exp[v(x',1)] \\big) \\pi(x'|x,d) dx'\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Expected value function\n",
    "\n",
    "Let $ \\mathbb{E}\\big[ V(x',\\varepsilon')\\big|x,d\\big] = EV(x,d) $, then\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "EV(x,d) &=& \\int_{X} \\log \\big( \\exp[v(x',0)] + \\exp[v(x',1)] \\big) \\pi(x'|x,d) dx' \\\\\n",
    "v(x,d) &=& u(x,d) + \\beta EV(x,d)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "- this is Bellman equation *in expected value function space*  \n",
    "- when the state space is discrete the integral is, of course, a simple sum over future values  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Bellman equation and Bellman operator in expected value function space\n",
    "\n",
    "$$\n",
    "EV(x,d) = \\sum_{X} \\log \\big( \\exp[u(x',0) + \\beta EV(x',0)] + \\exp[u(x',1) + \\beta EV(x',1)] \\big) \\pi(x'|x,d)\n",
    "$$\n",
    "\n",
    "$$\n",
    "T^*(EV)(x,d) \\equiv \\sum_{X} \\log \\big( \\exp[u(x',0) + \\beta EV(x',0)] + \\exp[u(x',1) + \\beta EV(x',1)] \\big) \\pi(x'|x,d)\n",
    "$$\n",
    "\n",
    "Solution to the Bellman functional equation $ EV(x,d) $ is also a fixed point of $ T^* $ operator, $ T^*(EV)(x,d)=EV(x,d) $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Benefits of Rust’s approach\n",
    "\n",
    "- Bellman operator in expected terms is also a contraction mapping $ \\Rightarrow $ VFI work!  \n",
    "- Dimentionality of this fixed point problem is smaller that the one in value function terms (because $ \\varepsilon $ is not included)  \n",
    "- It is also numerically easier to work with smooth expected values $ EV(x,d) $ rather than $ V(x,\\varepsilon) $  \n",
    "- Later we’ll also see a very nice numerical optimization possibilities as well  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Choice probabilities\n",
    "\n",
    "Once the fixed point is found, the *optimal* choice probability $ P(d|x) $ is given by the Logit structure (assumption EV):\n",
    "\n",
    "$$\n",
    "P(d|x) = \\frac{\\exp[v(x,d)]}{\\sum_{d'\\in \\{0,1\\}} \\exp[v(x,d')]}\n",
    "$$\n",
    "\n",
    "The choice probability serve as the bases for forming the likelihood function.  Will continue with this when talking about structural estimation!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Further learning resources\n",
    "\n",
    "- John Rust [https://en.wikipedia.org/wiki/John_Rust](https://en.wikipedia.org/wiki/John_Rust)  \n",
    "- Optimal Replacement of GMC Bus Engines: An Empirical Model of Harold Zurcher [https://www.jstor.org/stable/1911259](https://www.jstor.org/stable/1911259)  \n",
    "- Google scholar citing Rust (1987) [https://scholar.google.com/scholar?oi=bibs&hl=en&cites=16527795233338248687](https://scholar.google.com/scholar?oi=bibs&hl=en&cites=16527795233338248687)  "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "date": 1612589586.005288,
  "download_nb": false,
  "filename": "28_zurcher.rst",
  "filename_with_path": "28_zurcher",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "title": "Foundations of Computational Economics #28"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}