{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to Solve a Basic Dynamic Programming Problem\n", "\n", "This notebook will demonstrate how to solve a class of problems commonly encountered in economics. The application will be a simple income fluctuation problem that is widely known and forms the basis of many heterogeneous agents models. Hopefully, this code is straightforward to extend to new economic setups.\n", "\n", "There are many ways of solving these types of problems. The method we will discuss is the most intuitive and a great starting point for learning about more advanced techniques. Moreover, it is particularly well-suited for Julia. It will highlight the following features:\n", "\n", "* Fast loops\n", "* The type system\n", "* Some useful packages\n", "* Easy parallelization\n", "\n", "## Economic Environment\n", "\n", "Consider an infinitely-lived household who seeks to choose a plan for consumption $\\{c_t\\}_{t=1}^{\\infty}$ to maximize:\n", "$$\n", "\\mathbb{E} \\displaystyle \\sum_{t=0}^{\\infty} \\beta^t u(c_t)\n", "$$\n", "\n", "subject to \n", "\\begin{align*}\n", "c_t + a_{t+1} &\\leq (1 + r) a_t + w y_t \\\\\n", "c_t &\\geq 0 \\\\\n", "a_t &\\geq 0\n", "\\end{align*}\n", "\n", "The budget constraint says that the household needs to allocate resources to consumption and savings, $a_{t+1}$. The resources it has at its disposal are its savings plus interest, $(1+r) a_t$ and its labor income, $w y_t$. Furthermore, savings are subject to a borrowing limit -- here it's zero, meaning that households can't go into debt.\n", "\n", "The income process $y_t$ is assumed to follow a Markov chain with transition matrix $\\Pi$. For example, $\\Pi(y_j | y_i)$ is the probability that tomorrow's state is $y_j$ if today's state if $y_i$.\n", "\n", "## Dynamic Programming Setup\n", "\n", "We are looking for consumption plans and saving plans for the household under any possible state that can arise. What is a state here? It consists of:\n", "\n", "1. Income $y$: today's income determines how much cash-in-hand the household has and is used to forecast tomorrow's income. This is the **exogenous** state, and it is not chosen by the household.\n", "2. Current assets $a$: also determines how much resources are available. This is chosen by the household each period and gets carried over to next period. This is the **endogenous** state.\n", "\n", "Let's formulate the household's problem recursively. The Bellman equation is:\n", "\n", "$$\n", "V(a, y) = \\max_{c, a'} u(c) + \\beta \\displaystyle \\sum_{y' \\in Y} \\Pi(y' | y) V(a', y')\n", "$$\n", "\n", "subject to \n", "\\begin{align*}\n", "c + a' &\\leq (1 + r) a + w y \\\\\n", "c &\\geq 0 \\\\\n", "a &\\geq 0\n", "\\end{align*}\n", "\n", "We call $V(a, y)$ the **value function**. Each period, the household chooses assets and consumption to maximize the payoff from being in state $(a, y)$. This is the sum of two components:\n", "\n", "1. The instantaneous utility it gets in the current period\n", "2. The continuation value, which represents the expected value of tomorrow's state. This depends on the level of assets the household chooses today, and an expectation is taken over the exogenous income process.\n", "\n", "The solution to this problem will be a value function and associated policy functions for consumption and assets that satisfy the household's problem and the constraints.\n", "\n", "## Solution Strategy\n", "\n", "This solving this problem requires us to find an unknown *function*. To compute it, we will use an iterative method called **value function iteration**. We will start with some guess for the value function, compute the right-hand-side of the Bellman equation, and then use the new resulting function as our next guess. We will continue iterating in this fashion until our guesses converge. \n", "\n", "Since a computer can't store functions, we will need to approximate our guesses for value function. There are many ways to do this. Our approach here will be to record the value at a finite set of asset and income points, and linearly interpolate between those points. \n", "\n", "We will be more specific about the details of each step as we walk through the code. We'll first discuss how the algorithm works and then show how to parallelize it.\n", "\n", "## Implementation\n", "\n", "We first import the packages needed." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Plots.PyPlotBackend()" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using QuantEcon\n", "using Interpolations\n", "using Optim\n", "using KernelDensity\n", "using StatsBase\n", "using Plots\n", "pyplot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting Up the Infrastructure\n", "\n", "We are first going to set up a type, `Household`, that will store all of the parameters needed to solve the household's problem. In general, creating custom types is helpful for organizing and structuring large programs. Here, it is a clean way to pass all of the parameters to the functions we will define." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "mutable struct Household\n", "\n", " # User-inputted\n", " β::Float64 # discount factor\n", " y_chain::MarkovChain{Float64,Matrix{Float64},Vector{Float64}} # income process\n", " r::Float64 # interest rate\n", " w::Float64 # wage\n", " amax::Float64 # top of asset grid\n", " Na::Int64 # number of points on asset grid\n", " curv::Float64 # curvature of asset grid\n", "\n", " # Created in constructor\n", " a_grid::Vector{Float64} # asset grid\n", " y_grid::Vector{Float64} # income grid\n", " ay_grid::Matrix{Float64} # combined asset/income grid\n", " ayi_grid::Matrix{Float64} # combined asset/income index grid\n", " Ny::Int64 # number of points on income grid\n", " V::Matrix{Float64} # current guess for value function on grid points\n", " ap::Matrix{Float64} # current guess for asset policy function on grid points\n", " c::Matrix{Float64} # current guess for consumption policy function on grid points\n", "\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `Household` struct will store all of the exogenous parameters of the model, the grids over which we will approximate the value function, and the current guesses for the value function, asset policy function, and consumption policy function. Each of these comes paired with a *type declaration* (after the `::`). Julia does not require these, but including them avoids the need for the compiler to infer the types of each field. In many situations, this can dramatically improve the speed of the code.\n", "\n", "The first block of fields are the parameters the user needs to specify. The second block is created by the code in the constructor. Before describing that, let's define a few functions that will be useful throughout the algorithm" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "budget_constraint (generic function with 2 methods)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u(c::Float64) = log(c)\n", "\n", "function budget_constraint(a::Float64, ap::Float64, y::Float64, r::Float64, w::Float64)\n", " c = (1 + r)*a + w*y - ap\n", " return c\n", "end\n", "\n", "function budget_constraint(a::Float64, ap::Float64, y::Float64, h::Household)\n", " budget_constraint(a, ap, y, h.r, h.w)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We've set up the utility function, and also defined a function that backs out consumption from the budget constraint given a state (`a` and `y`) and a choice of assets tomorrow `ap`. The `budget_constraint` function has two methods, one that also accepts values for the interest rate and the wage, and another that accepts a `Household` type and takes the interest rate and wage from there. Based on the types of the arguments, Julia determines the appropriate method to execute. This is an example of *multiple dispatch*, one of the major features of Julia. It's useful tool for organizing different implementations of the same concept.\n", "\n", "Next, we define the constructor function for the `Household`:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Household" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function Household(;β::Float64=0.96,\n", " y_chain::MarkovChain{Float64,Matrix{Float64},Vector{Float64}}=MarkovChain([0.5 0.5; 0.04 0.96], [0.25; 1.0]),\n", " r::Float64=0.038, \n", " w::Float64=1.09,\n", " amax::Float64=30.0,\n", " Na::Int64=100,\n", " curv::Float64=0.4)\n", "\n", " # Set up asset grid\n", " a_grid = linspace(0, amax^curv, Na).^(1/curv)\n", "\n", " # Parameters of income grid\n", " y_grid = y_chain.state_values\n", " Ny = length(y_grid)\n", "\n", " # Combined grids\n", " ay_grid = gridmake(a_grid, y_grid)\n", " ayi_grid = gridmake(1:Na, 1:Ny)\n", "\n", " # Set up initial guess for value function\n", " V = zeros(Na, Ny)\n", " c = zeros(Na, Ny)\n", " for (yi, y) in enumerate(y_grid)\n", " for (ai, a) in enumerate(a_grid)\n", " c_max = budget_constraint(a, 0.0, y, r, w)\n", " c[ai, yi] = c_max\n", " V[ai, yi] = u(c_max)/(1 - β)\n", " end\n", " end\n", "\n", " # Corresponding initial policy function is all zeros\n", " ap = zeros(Na, Ny)\n", "\n", " return Household(β, y_chain, r, w, amax, Na, curv, a_grid, y_grid, ay_grid, ayi_grid, Ny, V, ap, c)\n", "\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notable features:\n", "* All inputs (parameters) are optional, with the given defaults\n", "* The asset grid is designed with more points closer to the borrowing constraint. In models with incomplete markets and a borrowing constraint, we expect there to be more curvature in this area. We want to have better precision in places where the policy is not as close to linear.\n", "* The initial guess used for the value function is the present discounted value of consuming the maximum posible amount in the given state forever. In principle any guess should work, but this is a good one because it gives the value function the right shape.\n", "\n", "### Applying the Bellman Operator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to be able to compute the right-hand side of the Bellman equation:\n", "\n", "$$\n", "V(a, y) = \\max_{c, a'} u(c) + \\beta \\displaystyle \\sum_{y' \\in Y} \\Pi(y' | y) V(a', y')\n", "$$\n", "\n", "subject to \n", "\\begin{align*}\n", "c + a' &\\leq (1 + r) a + w y \\\\\n", "c &\\geq 0 \\\\\n", "a &\\geq 0\n", "\\end{align*}\n", "\n", "For every grid point, we need to be able to compute this maximum (we will eliminate $c$ using the budget constraint and maximize over $a'$). This will involve using an optimization package. Because the optimal asset policy will not necessarily be on the grid, we need to be able to evaluate our guess for the value function at points off this grid. This will involve using an interpolation package.\n", "\n", "Let's start off by writing a function to compute the flow value and continuation value given the indices of a set of states, `ai` and `yi`, and a choice for tomorrow's asset level, `ap`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "value (generic function with 1 method)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function value(h::Household, itp_V::Interpolations.GriddedInterpolation,\n", " ap::Float64, ai::Int64, yi::Int64)\n", " \n", " # Interpolate value function at a', for each possible income level\n", " Vp = zeros(h.Ny)\n", " for yii = 1:h.Ny\n", " Vp[yii] = itp_V[ap, yii]\n", " end\n", "\n", " # Take expectations\n", " continuation = h.β*dot(h.y_chain.p[yi, :], Vp)\n", "\n", " # Compute today's consumption and utility\n", " c = budget_constraint(h.a_grid[ai], ap, h.y_grid[yi], h)\n", " flow = u(c)\n", "\n", " RHS = flow + continuation # This is what we want to maximize\n", " return RHS\n", "\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* The argument `itp_V` is an interpolation object that contains all the information necessary for interpolating the value function. We will show how it's created later, but the expression `itp_V[ap, yii]` returns the interpolated value at asset level `ap` and income index `yii`.\n", "* We take expectations by first creating a vector with the interpolated values at each possible income level with tomorrow's choice of assets. Then we take the dot product of that vector with the transition probabilities from the approprate row of the Markov transition matrix.\n", "\n", "Next, we need to compute the optimal asset policy. This amounts to maximizing the function above with respect to `ap`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "opt_value (generic function with 1 method)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function opt_value(h::Household, itp_V::Interpolations.GriddedInterpolation,\n", " ai::Int64, yi::Int64)\n", "\n", " f(ap::Float64) = -value(h, itp_V, ap, ai, yi)\n", " ub = (1 + h.r)*h.a_grid[ai] + h.w*h.y_grid[yi]\n", " lb = 0.0\n", " res = optimize(f, lb, ub)\n", " return res.minimizer, -res.minimum\n", "\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* The default optimizer from the `Optim` package works well here. \n", "* We define `f` based off of `value` with the current interpolant and the given state pair. \n", "* The lower bound on assets is the zero borrowing constraint, while the upper bound is the amount the household would save if it consumed nothing.\n", "* The function returns the argmax (asset policy), `res.minimizer`, and the maximum (the value), `-res.minimum`\n", "\n", "Next, we just need to be able to compute this max at every grid point in the state space. The function `bellman_iterator!` does this, as well as creates the interpolation object based off of the current value function guess." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "bellman_iterator! (generic function with 1 method)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function bellman_iterator!(h::Household)\n", "\n", " # Set up interpolant for current value function guess\n", " knots = (h.a_grid, [1, 2])\n", " itp_V = interpolate(knots, h.V, (Gridded(Linear()), NoInterp()))\n", "\n", " # Loop through the grid to update value and policy functions\n", " V = zeros(h.Na, h.Ny)\n", " ap = zeros(h.Na, h.Ny)\n", " c = zeros(h.Na, h.Ny)\n", " for (yi, y) in enumerate(h.y_grid)\n", " for (ai, a) in enumerate(h.a_grid)\n", " ap[ai, yi], V[ai, yi] = opt_value(h, itp_V, ai, yi)\n", " c[ai, yi] = budget_constraint(a, ap[ai, yi], y, h)\n", " end\n", " end\n", "\n", " # Update solution\n", " h.V = V;\n", " h.ap = ap;\n", " h.c = c;\n", " \n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* To create the interpolation object, we need to provide\n", " * The knots, which are the points for which we are providing the function value. Here, these are the grid points for assets and the indices of the income points (these are used as opposed to the points themselves because we are using the `NoInterp()` option, described below\n", " * The values of the function at the knots. This is the current guess for the value function `h.V`\n", " * The type of interpolation: here, we are using linear interpolation between the asset points. Since income only assumes values in a finite set, we use the `NoInterp()` option, which means the interpolation object won't have to save any information for that dimension\n", "* Once we have the interpolation object, we can loop through all the grid points and solve the optimization problem at each point. There is no need to vectorize the algorithm; Julia handles loops very well.\n", "* The exclamation point in the function name means that the function will modify its arugments; in this case, the `Household` object\n", "\n", "### Iteration\n", "\n", "Finally, we just need to build the infrastructure for iterating on the value function, updating the guess, and checking convergence." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "vfi! (generic function with 1 method)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function vfi!(h::Household; tol::Float64=1e-8, maxit::Int64=3000)\n", "\n", " dist = 1.0\n", " i = 1\n", " while (tol < dist) & (i < maxit)\n", " V_old = h.V\n", " bellman_iterator!(h)\n", " dist = maximum(abs, h.V - V_old)\n", "\n", " if i % 50 == 0\n", " println(\"Iteration $(i), distance $(dist)\")\n", " end\n", "\n", " i = i + 1\n", " end\n", "\n", " println(\"Converged in $(i) iterations!\")\n", "\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use the maximum absolute difference in the value functions (the sup-norm) to check for convergence. This norm is guaranteed to converge monotonically to the true value function.\n", "\n", "## Running the Code\n", "\n", "To run the code, all we need to do is create a `Household` object and run the `vfi!` function. Below, we do that with the default parameters:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 50, distance 0.38175976506264675\n", "Iteration 100, distance 0.07087936221186553\n", "Iteration 150, distance 0.014663718899711853\n", "Iteration 200, distance 0.0008171100387741603\n", "Iteration 250, distance 9.04814686748523e-5\n", "Iteration 300, distance 1.1324895226039189e-5\n", "Iteration 350, distance 1.450412696613057e-6\n", "Iteration 400, distance 1.876855009186329e-7\n", "Iteration 450, distance 2.4358325134699044e-8\n", "Converged in 473 iterations!\n", " 0.957675 seconds (6.58 M allocations: 418.709 MiB, 13.70% gc time)\n" ] } ], "source": [ "h = Household();\n", "@time vfi!(h)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function plot_policies(h::Household)\n", " \n", " labs = reshape([\"Unemployed\", \"Employed\"], 1, 2)\n", " plt_v = plot(h.a_grid, h.V, label=labs, lw=2, title = \"Value Function\", xlabel=\"Assets Today\")\n", " plt_ap = plot(h.a_grid, h.ap, label=labs, lw=2, title = \"Asset Policy Function\", xlabel=\"Assets Today\")\n", " plot!(h.a_grid, h.a_grid, color=:black, linestyle=:dash, lw=0.5, label=\"\")\n", " plt_c = plot(h.a_grid, h.c, label=labs, lw=2, title = \"Consumption Policy Function\", xlabel=\"Assets Today\")\n", " plot(plt_v, plt_ap, plt_c, layout=(1, 3), size=(1000, 400))\n", " \n", "end\n", "\n", "plot_policies(h)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 50, distance 0.38149422466877425\n", "Iteration 100, distance 0.06966792637230412\n", "Iteration 150, distance 0.013525481521678984\n", "Iteration 200, distance 0.000739719665276084\n", "Iteration 250, distance 8.977807224042067e-5\n", "Iteration 300, distance 1.1276834968043659e-5\n", "Iteration 350, distance 1.449225191407777e-6\n", "Iteration 400, distance 1.8781478061669077e-7\n", "Iteration 450, distance 2.438638802004789e-8\n", "Converged in 473 iterations!\n", " 71.028029 seconds (459.38 M allocations: 28.539 GiB, 10.45% gc time)\n" ] } ], "source": [ "h = Household(;Na=7000);\n", "@time vfi!(h)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparitive Statics Exercise: Part 1\n", "\n", "We are going to lower the level of unemployment benefits from 0.25 to 0.15 and examine how households' decision rules change. Another way to think of this is that we are increasing the variance of the income process. Using the functions above:\n", "* Create a new `Household` object with the new parameters (500 points on the asset grid is more than enough)\n", "* Solve the problem\n", "* Plot the consumption policy functions against the ones from the original setup\n", "\n", "\n", "