aiyagari.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "44072df1",
   "metadata": {},
   "source": [
    "\n",
    "<a id='aiyagari'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d16262e",
   "metadata": {},
   "source": [
    "# The Aiyagari Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5279fc1d",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [The Aiyagari Model](#The-Aiyagari-Model)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [The Economy](#The-Economy)  \n",
    "  - [Firms](#Firms)  \n",
    "  - [Code](#Code)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b3b714b",
   "metadata": {},
   "source": [
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc344a7b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6f5888d",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this lecture, we describe the structure of a class of models that build on work by Truman Bewley [[Bewley, 1977](https://python.quantecon.org/zreferences.html#id177)].\n",
    "\n",
    "We begin by discussing an example of a Bewley model due to Rao Aiyagari [[Aiyagari, 1994](https://python.quantecon.org/zreferences.html#id141)].\n",
    "\n",
    "The model features\n",
    "\n",
    "- Heterogeneous agents  \n",
    "- A single exogenous vehicle for borrowing and lending  \n",
    "- Limits on amounts individual agents may borrow  \n",
    "\n",
    "\n",
    "The Aiyagari model has been used to investigate many topics, including\n",
    "\n",
    "- precautionary savings and the effect of liquidity constraints [[Aiyagari, 1994](https://python.quantecon.org/zreferences.html#id141)]  \n",
    "- risk sharing and asset pricing [[Heaton and Lucas, 1996](https://python.quantecon.org/zreferences.html#id133)]  \n",
    "- the shape of the wealth distribution [[Benhabib *et al.*, 2015](https://python.quantecon.org/zreferences.html#id134)]  \n",
    "- etc., etc., etc.  \n",
    "\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "234965c8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)  #set default figure size\n",
    "import numpy as np\n",
    "from quantecon.markov import DiscreteDP\n",
    "from numba import jit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97851b18",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "The primary reference for this lecture is [[Aiyagari, 1994](https://python.quantecon.org/zreferences.html#id141)].\n",
    "\n",
    "A textbook treatment is available in chapter 18 of [[Ljungqvist and Sargent, 2018](https://python.quantecon.org/zreferences.html#id186)].\n",
    "\n",
    "A continuous time version of the model by SeHyoun Ahn and Benjamin Moll can be found [here](https://nbviewer.org/github/QuantEcon/QuantEcon.notebooks/blob/master/aiyagari_continuous_time.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d0ba9f",
   "metadata": {},
   "source": [
    "## The Economy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab6f5ba4",
   "metadata": {},
   "source": [
    "### Households\n",
    "\n",
    "Infinitely lived households / consumers face idiosyncratic income shocks.\n",
    "\n",
    "A unit interval of  *ex-ante* identical households face a common borrowing constraint.\n",
    "\n",
    "The savings problem faced by a typical  household is\n",
    "\n",
    "$$\n",
    "\\max \\mathbb E \\sum_{t=0}^{\\infty} \\beta^t u(c_t)\n",
    "$$\n",
    "\n",
    "subject to\n",
    "\n",
    "$$\n",
    "a_{t+1} + c_t \\leq w z_t + (1 + r) a_t\n",
    "\\quad\n",
    "c_t \\geq 0,\n",
    "\\quad \\text{and} \\quad\n",
    "a_t \\geq -B\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ c_t $ is current consumption  \n",
    "- $ a_t $ is assets  \n",
    "- $ z_t $ is an exogenous component of labor income capturing stochastic unemployment risk, etc.  \n",
    "- $ w $ is a wage rate  \n",
    "- $ r $ is a net interest rate  \n",
    "- $ B $ is the maximum amount that the agent is allowed to borrow  \n",
    "\n",
    "\n",
    "The exogenous process $ \\{z_t\\} $ follows a finite state Markov chain with given stochastic matrix $ P $.\n",
    "\n",
    "The wage and interest rate are fixed over time.\n",
    "\n",
    "In this simple version of the model, households supply labor  inelastically because they do not value leisure."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5433965b",
   "metadata": {},
   "source": [
    "## Firms\n",
    "\n",
    "Firms produce output by hiring capital and labor.\n",
    "\n",
    "Firms act competitively and face constant returns to scale.\n",
    "\n",
    "Since returns to scale are constant the number of firms does not matter.\n",
    "\n",
    "Hence we can consider a single (but nonetheless competitive) representative firm.\n",
    "\n",
    "The firm’s output is\n",
    "\n",
    "$$\n",
    "Y_t = A K_t^{\\alpha} N^{1 - \\alpha}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ A $ and $ \\alpha $ are parameters with $ A > 0 $ and $ \\alpha \\in (0, 1) $  \n",
    "- $ K_t $ is aggregate capital  \n",
    "- $ N $ is total labor supply (which is constant in this simple version of the model)  \n",
    "\n",
    "\n",
    "The firm’s problem is\n",
    "\n",
    "$$\n",
    "max_{K, N} \\left\\{ A K_t^{\\alpha} N^{1 - \\alpha} - (r + \\delta) K - w N \\right\\}\n",
    "$$\n",
    "\n",
    "The parameter $ \\delta $ is the depreciation rate.\n",
    "\n",
    "From the first-order condition with respect to capital, the firm’s inverse demand for capital is\n",
    "\n",
    "\n",
    "<a id='equation-aiy-rgk'></a>\n",
    "$$\n",
    "r = A \\alpha  \\left( \\frac{N}{K} \\right)^{1 - \\alpha} - \\delta \\tag{67.1}\n",
    "$$\n",
    "\n",
    "Using this expression and the firm’s first-order condition for labor, we can pin down\n",
    "the equilibrium wage rate as a function of $ r $ as\n",
    "\n",
    "\n",
    "<a id='equation-aiy-wgr'></a>\n",
    "$$\n",
    "w(r) = A  (1 - \\alpha)  (A \\alpha / (r + \\delta))^{\\alpha / (1 - \\alpha)} \\tag{67.2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41523b76",
   "metadata": {},
   "source": [
    "### Equilibrium\n",
    "\n",
    "We construct  a *stationary rational expectations equilibrium* (SREE).\n",
    "\n",
    "In such an equilibrium\n",
    "\n",
    "- prices induce behavior that generates aggregate quantities consistent with the prices  \n",
    "- aggregate quantities and prices are constant over time  \n",
    "\n",
    "\n",
    "In more detail, an SREE lists a set of prices, savings and production policies such that\n",
    "\n",
    "- households want to choose the specified savings policies taking the prices as given  \n",
    "- firms maximize profits taking the same prices as given  \n",
    "- the resulting aggregate quantities are consistent with the prices; in particular, the demand for capital equals the supply  \n",
    "- aggregate quantities (defined as cross-sectional averages) are constant  \n",
    "\n",
    "\n",
    "In practice, once parameter values are set, we can check for an SREE by the following steps\n",
    "\n",
    "1. pick a proposed quantity $ K $ for aggregate capital  \n",
    "1. determine corresponding prices, with interest rate $ r $ determined by [(67.1)](#equation-aiy-rgk) and a wage rate $ w(r) $ as given in [(67.2)](#equation-aiy-wgr)  \n",
    "1. determine the common optimal savings policy of the households given these prices  \n",
    "1. compute aggregate capital as the mean of steady state capital given this savings policy  \n",
    "\n",
    "\n",
    "If this final quantity agrees with $ K $ then we have a SREE."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "395e779e",
   "metadata": {},
   "source": [
    "## Code\n",
    "\n",
    "Let’s look at how we might compute such an equilibrium in practice.\n",
    "\n",
    "To solve the household’s dynamic programming problem we’ll use the [DiscreteDP](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/markov/ddp.py) class from [QuantEcon.py](https://quantecon.org/quantecon-py/).\n",
    "\n",
    "Our first task is the least exciting one: write code that maps parameters for a household problem into the `R` and `Q` matrices needed to generate an instance of `DiscreteDP`.\n",
    "\n",
    "Below is a piece of boilerplate code that does just this.\n",
    "\n",
    "In reading the code, the following information will be helpful\n",
    "\n",
    "- `R` needs to be a matrix where `R[s, a]` is the reward at state `s` under action `a`.  \n",
    "- `Q` needs to be a three-dimensional array where `Q[s, a, s']` is the probability of transitioning to state `s'` when the current state is `s` and the current action is `a`.  \n",
    "\n",
    "\n",
    "(A more detailed discussion of `DiscreteDP` is available in the [Discrete State Dynamic Programming](https://python-advanced.quantecon.org/discrete_dp.html) lecture in the [Advanced\n",
    "Quantitative Economics with Python](https://python-advanced.quantecon.org) lecture series.)\n",
    "\n",
    "Here we take the state to be $ s_t := (a_t, z_t) $, where $ a_t $ is assets and $ z_t $ is the shock.\n",
    "\n",
    "The action is the choice of next period asset level $ a_{t+1} $.\n",
    "\n",
    "We use Numba to speed up the loops so we can update the matrices efficiently\n",
    "when the parameters change.\n",
    "\n",
    "The class also includes a default set of parameters that we’ll adopt unless otherwise specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9701614d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "class Household:\n",
    "    \"\"\"\n",
    "    This class takes the parameters that define a household asset accumulation\n",
    "    problem and computes the corresponding reward and transition matrices R\n",
    "    and Q required to generate an instance of DiscreteDP, and thereby solve\n",
    "    for the optimal policy.\n",
    "\n",
    "    Comments on indexing: We need to enumerate the state space S as a sequence\n",
    "    S = {0, ..., n}.  To this end, (a_i, z_i) index pairs are mapped to s_i\n",
    "    indices according to the rule\n",
    "\n",
    "        s_i = a_i * z_size + z_i\n",
    "\n",
    "    To invert this map, use\n",
    "\n",
    "        a_i = s_i // z_size  (integer division)\n",
    "        z_i = s_i % z_size\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    def __init__(self,\n",
    "                r=0.01,                      # Interest rate\n",
    "                w=1.0,                       # Wages\n",
    "                β=0.96,                      # Discount factor\n",
    "                a_min=1e-10,\n",
    "                Π=[[0.9, 0.1], [0.1, 0.9]],  # Markov chain\n",
    "                z_vals=[0.1, 1.0],           # Exogenous states\n",
    "                a_max=18,\n",
    "                a_size=200):\n",
    "\n",
    "        # Store values, set up grids over a and z\n",
    "        self.r, self.w, self.β = r, w, β\n",
    "        self.a_min, self.a_max, self.a_size = a_min, a_max, a_size\n",
    "\n",
    "        self.Π = np.asarray(Π)\n",
    "        self.z_vals = np.asarray(z_vals)\n",
    "        self.z_size = len(z_vals)\n",
    "\n",
    "        self.a_vals = np.linspace(a_min, a_max, a_size)\n",
    "        self.n = a_size * self.z_size\n",
    "\n",
    "        # Build the array Q\n",
    "        self.Q = np.zeros((self.n, a_size, self.n))\n",
    "        self.build_Q()\n",
    "\n",
    "        # Build the array R\n",
    "        self.R = np.empty((self.n, a_size))\n",
    "        self.build_R()\n",
    "\n",
    "    def set_prices(self, r, w):\n",
    "        \"\"\"\n",
    "        Use this method to reset prices. Calling the method will trigger a\n",
    "        re-build of R.\n",
    "        \"\"\"\n",
    "        self.r, self.w = r, w\n",
    "        self.build_R()\n",
    "\n",
    "    def build_Q(self):\n",
    "        populate_Q(self.Q, self.a_size, self.z_size, self.Π)\n",
    "\n",
    "    def build_R(self):\n",
    "        self.R.fill(-np.inf)\n",
    "        populate_R(self.R,\n",
    "                self.a_size,\n",
    "                self.z_size,\n",
    "                self.a_vals,\n",
    "                self.z_vals,\n",
    "                self.r,\n",
    "                self.w)\n",
    "\n",
    "\n",
    "# Do the hard work using JIT-ed functions\n",
    "\n",
    "@jit(nopython=True)\n",
    "def populate_R(R, a_size, z_size, a_vals, z_vals, r, w):\n",
    "    n = a_size * z_size\n",
    "    for s_i in range(n):\n",
    "        a_i = s_i // z_size\n",
    "        z_i = s_i % z_size\n",
    "        a = a_vals[a_i]\n",
    "        z = z_vals[z_i]\n",
    "        for new_a_i in range(a_size):\n",
    "            a_new = a_vals[new_a_i]\n",
    "            c = w * z + (1 + r) * a - a_new\n",
    "            if c > 0:\n",
    "                R[s_i, new_a_i] = np.log(c)  # Utility\n",
    "\n",
    "@jit(nopython=True)\n",
    "def populate_Q(Q, a_size, z_size, Π):\n",
    "    n = a_size * z_size\n",
    "    for s_i in range(n):\n",
    "        z_i = s_i % z_size\n",
    "        for a_i in range(a_size):\n",
    "            for next_z_i in range(z_size):\n",
    "                Q[s_i, a_i, a_i*z_size + next_z_i] = Π[z_i, next_z_i]\n",
    "\n",
    "\n",
    "@jit(nopython=True)\n",
    "def asset_marginal(s_probs, a_size, z_size):\n",
    "    a_probs = np.zeros(a_size)\n",
    "    for a_i in range(a_size):\n",
    "        for z_i in range(z_size):\n",
    "            a_probs[a_i] += s_probs[a_i*z_size + z_i]\n",
    "    return a_probs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4f3b13f",
   "metadata": {},
   "source": [
    "As a first example of what we can do, let’s compute and plot an optimal accumulation policy at fixed prices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6883a6ba",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Example prices\n",
    "r = 0.03\n",
    "w = 0.956\n",
    "\n",
    "# Create an instance of Household\n",
    "am = Household(a_max=20, r=r, w=w)\n",
    "\n",
    "# Use the instance to build a discrete dynamic program\n",
    "am_ddp = DiscreteDP(am.R, am.Q, am.β)\n",
    "\n",
    "# Solve using policy function iteration\n",
    "results = am_ddp.solve(method='policy_iteration')\n",
    "\n",
    "# Simplify names\n",
    "z_size, a_size = am.z_size, am.a_size\n",
    "z_vals, a_vals = am.z_vals, am.a_vals\n",
    "n = a_size * z_size\n",
    "\n",
    "# Get all optimal actions across the set of a indices with z fixed in each row\n",
    "a_star = np.empty((z_size, a_size))\n",
    "for s_i in range(n):\n",
    "    a_i = s_i // z_size\n",
    "    z_i = s_i % z_size\n",
    "    a_star[z_i, a_i] = a_vals[results.sigma[s_i]]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 9))\n",
    "ax.plot(a_vals, a_vals, 'k--')  # 45 degrees\n",
    "for i in range(z_size):\n",
    "    lb = f'$z = {z_vals[i]:.2}$'\n",
    "    ax.plot(a_vals, a_star[i, :], lw=2, alpha=0.6, label=lb)\n",
    "    ax.set_xlabel('current assets')\n",
    "    ax.set_ylabel('next period assets')\n",
    "ax.legend(loc='upper left')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82ffb4c1",
   "metadata": {},
   "source": [
    "The plot shows asset accumulation policies at different values of the exogenous state.\n",
    "\n",
    "Now we want to calculate the equilibrium.\n",
    "\n",
    "Let’s do this visually as a first pass.\n",
    "\n",
    "The following code draws aggregate supply and demand curves.\n",
    "\n",
    "The intersection gives equilibrium interest rates and capital."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76dd0661",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A = 1.0\n",
    "N = 1.0\n",
    "α = 0.33\n",
    "β = 0.96\n",
    "δ = 0.05\n",
    "\n",
    "\n",
    "def r_to_w(r):\n",
    "    \"\"\"\n",
    "    Equilibrium wages associated with a given interest rate r.\n",
    "    \"\"\"\n",
    "    return A * (1 - α) * (A * α / (r + δ))**(α / (1 - α))\n",
    "\n",
    "def rd(K):\n",
    "    \"\"\"\n",
    "    Inverse demand curve for capital.  The interest rate associated with a\n",
    "    given demand for capital K.\n",
    "    \"\"\"\n",
    "    return A * α * (N / K)**(1 - α) - δ\n",
    "\n",
    "\n",
    "def prices_to_capital_stock(am, r):\n",
    "    \"\"\"\n",
    "    Map prices to the induced level of capital stock.\n",
    "\n",
    "    Parameters:\n",
    "    ----------\n",
    "\n",
    "    am : Household\n",
    "        An instance of an aiyagari_household.Household\n",
    "    r : float\n",
    "        The interest rate\n",
    "    \"\"\"\n",
    "    w = r_to_w(r)\n",
    "    am.set_prices(r, w)\n",
    "    aiyagari_ddp = DiscreteDP(am.R, am.Q, β)\n",
    "    # Compute the optimal policy\n",
    "    results = aiyagari_ddp.solve(method='policy_iteration')\n",
    "    # Compute the stationary distribution\n",
    "    stationary_probs = results.mc.stationary_distributions[0]\n",
    "    # Extract the marginal distribution for assets\n",
    "    asset_probs = asset_marginal(stationary_probs, am.a_size, am.z_size)\n",
    "    # Return K\n",
    "    return np.sum(asset_probs * am.a_vals)\n",
    "\n",
    "\n",
    "# Create an instance of Household\n",
    "am = Household(a_max=20)\n",
    "\n",
    "# Use the instance to build a discrete dynamic program\n",
    "am_ddp = DiscreteDP(am.R, am.Q, am.β)\n",
    "\n",
    "# Create a grid of r values at which to compute demand and supply of capital\n",
    "num_points = 20\n",
    "r_vals = np.linspace(0.005, 0.04, num_points)\n",
    "\n",
    "# Compute supply of capital\n",
    "k_vals = np.empty(num_points)\n",
    "for i, r in enumerate(r_vals):\n",
    "    k_vals[i] = prices_to_capital_stock(am, r)\n",
    "\n",
    "# Plot against demand for capital by firms\n",
    "fig, ax = plt.subplots(figsize=(11, 8))\n",
    "ax.plot(k_vals, r_vals, lw=2, alpha=0.6, label='supply of capital')\n",
    "ax.plot(k_vals, rd(k_vals), lw=2, alpha=0.6, label='demand for capital')\n",
    "ax.grid()\n",
    "ax.set_xlabel('capital')\n",
    "ax.set_ylabel('interest rate')\n",
    "ax.legend(loc='upper right')\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "date": 1727927324.0410063,
  "filename": "aiyagari.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "The Aiyagari Model"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}