{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Numerical Dynamic Programming\n",
" \n",
"Florian Oswald, Sciences Po, 2019\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"## Intro\n",
"\n",
"* Numerical Dynamic Programming (DP) is widely used to solve dynamic models.\n",
"* You are familiar with the technique from your core macro course.\n",
"* We will illustrate some ways to solve dynamic programs.\n",
" 1. Models with one discrete or continuous choice variable\n",
" 1. Models with several choice variables\n",
"\t1. Models with a discrete-continuous choice combination\n",
"* We will go through:\n",
"\t1. Value Function Iteration (VFI)\n",
" 1. Policy function iteration (PFI)\n",
" 1. Projection Methods\n",
"\t1. Endogenous Grid Method (EGM)\n",
"\t1. Discrete Choice Endogenous Grid Method (DCEGM)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"## Dynamic Programming Theory\n",
"\n",
"* Payoffs over time are \n",
"\t$$U=\\sum_{t=1}^{\\infty}\\beta^{t}u\\left(s_{t},c_{t}\\right) $$\n",
"\twhere $\\beta<1$ is a discount factor, $s_{t}$ is the state, $c_{t}$ is the control.\n",
"\n",
"* The state (vector) evolves as $s_{t+1}=h(s_{t},c_{t})$.\n",
"* All past decisions are contained in $s_{t}$.\n",
"\n",
"### Assumptions\n",
"\n",
"* Let $c_{t}\\in C(s_{t}),s_{t}\\in S$ and assume $u$ is bounded in $(c,s)\\in C\\times S$.\n",
"* Stationarity: neither payoff $u$ nor transition $h$ depend on time.\n",
"* Write the problem as \n",
"\t$$ v(s)=\\max_{s'\\in\\Gamma(s)}u(s,s')+\\beta v(s') $$\n",
"* $\\Gamma(s)$ is the constraint set (or feasible set) for $s'$ when the current state is $s$\n",
"\n",
"### Existence\n",
"\n",
"**Theorem.** Assume that $u(s,s')$ is real-valued, continuous, and bounded, that $\\beta\\in(0,1)$, and that the constraint set $\\Gamma(s)$ is nonempty, compact, and continuous. Then there exists a unique function $v(s)$ that solves the above functional equation.\n",
"\n",
"**Proof.** [@stokeylucas] theoreom 4.6."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Solution Methods\n",
"\n",
"## Value Function Iteration (VFI)\n",
"\n",
"* Find the fix point of the functional equation by iterating on it until the distance between consecutive iterations becomes small.\n",
"* Motivated by the Bellman Operator, and it's characterization in the Continuous Mapping Theorem.\n",
"\n",
"## Discrete DP VFI\n",
"\n",
"* Represents and solves the functional problem in $\\mathbb{R}$ on a finite set of grid points only.\n",
"* Widely used method.\n",
"\t* Simple (+)\n",
"\t* Robust (+)\n",
"\t* Slow (-)\n",
"\t* Imprecise (-)\n",
"* Precision depends on number of discretization points used. \n",
"* High-dimensional problems are difficult to tackle with this method because of the curse of dimensionality.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Deterministic growth model with Discrete VFI\n",
"\n",
"* We have this theoretical model:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" V(k) &= \\max_{0\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Bellman Operator\n",
"# inputs\n",
"# `grid`: grid of values of state variable\n",
"# `v0`: current guess of value function\n",
"\n",
"# output\n",
"# `v1`: next guess of value function\n",
"# `pol`: corresponding policy function \n",
"\n",
"#takes a grid of state variables and computes the next iterate of the value function.\n",
"function bellman_operator(grid,v0)\n",
" \n",
" v1 = zeros(n) # next guess\n",
" pol = zeros(Int,n) # policy function\n",
" w = zeros(n) # temporary vector \n",
"\n",
" # loop over current states\n",
" # current capital\n",
" for (i,k) in enumerate(grid)\n",
"\n",
" # loop over all possible kprime choices\n",
" for (iprime,kprime) in enumerate(grid)\n",
" if f(k) - kprime < 0 #check for negative consumption\n",
" w[iprime] = -Inf\n",
" else\n",
" w[iprime] = ufun(f(k) - kprime) + beta * v0[iprime]\n",
" end\n",
" end\n",
" # find maximal choice\n",
" v1[i], pol[i] = findmax(w) # stores Value und policy (index of optimal choice)\n",
" end\n",
" return (v1,pol) # return both value and policy function\n",
"end\n",
"\n",
"\n",
"\n",
"# VFI iterator\n",
"#\n",
"## input\n",
"# `n`: number of grid points\n",
"# output\n",
"# `v_next`: tuple with value and policy functions after `n` iterations.\n",
"function VFI()\n",
" v_init = zeros(n) # initial guess\n",
" for iter in 1:N_iter\n",
" v_next = bellman_operator(kgrid,v_init) # returns a tuple: (v1,pol)\n",
" # check convergence\n",
" if maximum(abs,v_init.-v_next[1]) < tol\n",
" verrors = maximum(abs,v_next[1].-v_star(kgrid))\n",
" perrors = maximum(abs,kgrid[v_next[2]].-k_star(kgrid))\n",
" println(\"Found solution after $iter iterations\")\n",
" println(\"maximal value function error = $verrors\")\n",
" println(\"maximal policy function error = $perrors\")\n",
" return v_next\n",
" elseif iter==N_iter\n",
" warn(\"No solution found after $iter iterations\")\n",
" return v_next\n",
" end\n",
" v_init = v_next[1] # update guess \n",
" end\n",
"end\n",
"\n",
"# plot\n",
"using Plots\n",
"function plotVFI()\n",
" v = VFI()\n",
" p = Any[]\n",
" \n",
" # value and policy functions\n",
" push!(p,plot(kgrid,v[1],\n",
" lab=\"V\",\n",
" ylim=(-50,-30),legend=:bottomright),\n",
" plot(kgrid,kgrid[v[2]],\n",
" lab=\"policy\",legend=:bottomright))\n",
" \n",
" # errors of both\n",
" push!(p,plot(kgrid,v[1].-v_star(kgrid),\n",
" lab=\"V error\",legend=:bottomright),\n",
" plot(kgrid,kgrid[v[2]].-k_star(kgrid),\n",
" lab=\"policy error\",legend=:bottomright))\n",
"\n",
" plot(p...,layout=grid(2,2) )\n",
" \n",
"end\n",
"\n",
"plotVFI()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Exercise 2: Discretizing only the state space (not control space)\n",
"\n",
"* Same exercise, but now use a continuous solver for choice of $k'$.\n",
"* in other words, employ the following numerical approximation:\n",
"\t$$ V(k_i) = \\max_{k'\\in[0,\\bar{k}]} \\ln(f(k_i) - k') + \\beta V(k') $$\n",
"* To do this, you need to be able to evaluate $V(k')$ where $k'$ is potentially off the `kgrid`.\n",
"* use `Interpolations.jl` to linearly interpolate V.\n",
" * the relevant object is setup with function `interpolate((grid,),v,Gridded(Linear()))`\n",
"* use `Optim::optimize()` to perform the maximization.\n",
" * you have to define an ojbective function for each $k_i$\n",
" * do something like `optimize(objective, lb,ub)` "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.01:0.013355704697986578:2.0"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kgrid"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"attributes": {
"classes": [
"julia"
],
"id": ""
},
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"continuous VFI:\n",
"Found solution after 418 iterations\n",
"maximal value function error = 0.04828453368161689\n",
"maximal policy function error = 0.004602693711777683\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Interpolations\n",
"using Optim\n",
"function bellman_operator2(grid,v0)\n",
" \n",
" v1 = zeros(n) # next guess\n",
" pol = zeros(n) # consumption policy function\n",
"\n",
" Interp = interpolate((collect(grid),), v0, Gridded(Linear()) ) \n",
" Interp = extrapolate(Interp,Interpolations.Flat())\n",
"\n",
" # loop over current states\n",
" # of current capital\n",
" for (i,k) in enumerate(grid)\n",
"\n",
" objective(c) = - (log.(c) + beta * Interp(f(k) - c))\n",
" # find max of ojbective between [0,k^alpha]\n",
" res = optimize(objective, 1e-6, f(k)) # Optim.jl\n",
" pol[i] = f(k) - res.minimizer # k'\n",
" v1[i] = -res.minimum\n",
" end\n",
" return (v1,pol) # return both value and policy function\n",
"end\n",
"\n",
"function VFI2()\n",
" v_init = zeros(n) # initial guess\n",
" for iter in 1:N_iter\n",
" v_next = bellman_operator2(kgrid,v_init) # returns a tuple: (v1,pol)\n",
" # check convergence\n",
" if maximum(abs,v_init.-v_next[1]) < tol\n",
" verrors = maximum(abs,v_next[1].-v_star(kgrid))\n",
" perrors = maximum(abs,v_next[2].-k_star(kgrid))\n",
" println(\"continuous VFI:\")\n",
" println(\"Found solution after $iter iterations\")\n",
" println(\"maximal value function error = $verrors\")\n",
" println(\"maximal policy function error = $perrors\")\n",
" return v_next\n",
" elseif iter==N_iter\n",
" warn(\"No solution found after $iter iterations\")\n",
" return v_next\n",
" end\n",
" v_init = v_next[1] # update guess \n",
" end\n",
" return nothing\n",
"end\n",
"\n",
"function plotVFI2()\n",
" v = VFI2()\n",
" p = Any[]\n",
" \n",
" # value and policy functions\n",
" push!(p,plot(kgrid,v[1],\n",
" lab=\"V\",\n",
" ylim=(-50,-30),legend=:bottomright),\n",
" plot(kgrid,v[2],\n",
" lab=\"policy\",legend=:bottomright))\n",
" \n",
" # errors of both\n",
" push!(p,plot(kgrid,v[1].-v_star(kgrid),\n",
" lab=\"V error\",legend=:bottomright),\n",
" plot(kgrid,v[2].-k_star(kgrid),\n",
" lab=\"policy error\",legend=:bottomright))\n",
"\n",
" plot(p...,layout=grid(2,2) )\n",
" \n",
"end\n",
"\n",
"plotVFI2()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Policy Function Iteration\n",
"\n",
"* This is similar to VFI but we now guess successive *policy* functions\n",
"* The idea is to choose a new policy $p^*$ in each iteration so as to satisfy an optimality condition. In our example, that would be the Euler Equation.\n",
"* We know that the solution to the above problem is a function $c^*(k)$ such that\n",
"\n",
"$$ \n",
" c^*(k) = \\arg \\max_z u(z) + \\beta V(f(k)-z) \\text{ }\\forall k\\in[0,\\infty]\n",
"$$\n",
"\n",
"* We **don't** directly solve the maximiation problem outlined above, but it's first order condition:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"u'(c^*(k_t)) & = \\beta u'(c^*(k_{t+1})) f'(k_{t+1}) \\\\\n",
" & = \\beta u'[c^*(f(k_{t})-c^*(k_t))] f'(f(k_{t})-c^*(k_t)) \n",
"\\end{aligned}\n",
"$$\n",
"\n",
"* In practice, we have to find the zeros of\n",
"\n",
"$$\n",
"g(k_t) = u'(c^*(k_t)) - \\beta u'[c^*(f(k_{t})-c^*(k_t))] f'(f(k_{t})-c^*(k_t)) \n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PFI:\n",
"Found solution after 39 iterations\n",
"max policy function error = 7.301895796647112e-5\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your turn!\n",
"\n",
"\n",
"using Roots\n",
"function policy_iter(grid,c0,u_prime,f_prime)\n",
" \n",
" c1 = zeros(length(grid)) # next guess\n",
" pol_fun = extrapolate(interpolate((collect(grid),), c0, Gridded(Linear()) ) , Interpolations.Flat())\n",
" \n",
" \n",
" # loop over current states\n",
" # of current capital\n",
" for (i,k) in enumerate(grid)\n",
" objective(c) = u_prime(c) - beta * u_prime(pol_fun(f(k)-c)) * f_prime(f(k)-c)\n",
" c1[i] = fzero(objective, 1e-10, f(k)-1e-10) \n",
" end\n",
" return c1\n",
"end\n",
"\n",
"uprime(x) = 1.0 ./ x\n",
"fprime(x) = alpha * x.^(alpha-1)\n",
"\n",
"function PFI()\n",
" c_init = kgrid\n",
" for iter in 1:N_iter\n",
" c_next = policy_iter(kgrid,c_init,uprime,fprime) \n",
" # check convergence\n",
" if maximum(abs,c_init.-c_next) < tol\n",
" perrors = maximum(abs,c_next.-c_star(kgrid))\n",
" println(\"PFI:\")\n",
" println(\"Found solution after $iter iterations\")\n",
" println(\"max policy function error = $perrors\")\n",
" return c_next\n",
" elseif iter==N_iter\n",
" warn(\"No solution found after $iter iterations\")\n",
" return c_next\n",
" end\n",
" c_init = c_next # update guess \n",
" end\n",
"end\n",
"function plotPFI()\n",
" v = PFI()\n",
" plot(kgrid,[v v.-c_star(kgrid)],\n",
" lab=[\"policy\" \"error\"],\n",
" legend=:bottomright,\n",
" layout = 2)\n",
"end\n",
"plotPFI()\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Projection Methods\n",
"\n",
"* Many applications require us to solve for an *unknown function*\n",
" * ODEs, PDEs\n",
" * Pricing functions in asset pricing models\n",
" * Consumption/Investment policy functions\n",
"* Projection methods find approximations to those functions that set an error function close to zero."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Example: Growth, again\n",
"\n",
"* We stick to our working example from above.\n",
"* We encountered the Euler Equation $g$ for optimality.\n",
"* At the true consumption function $c^*$, $g(k) = 0$.\n",
"* We define the following function operator:\n",
"\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"0 & = u'(c^*(k)) - \\beta u'[c^*(f(k)-c^*(k))] f'(f(k)-c^*(k)) \\\\\n",
" & \\equiv (\\mathcal{N(c^*)})(k)\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"* The Equilibrium solves the operator equation\n",
"$$\n",
"0 = \\mathcal{N}(c^*)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Projection Method example\n",
"\n",
"1. create an approximation to $c^*$: \n",
" find \n",
" $$\n",
" \\bar{c} \\equiv \\sum_{i=0}^n a_i k^i\n",
" $$\n",
" \n",
" which nearly solves \n",
" \n",
" $$\\mathcal{N}(c^*)=0\n",
" $$\n",
"2. Compute Euler equation error function:\n",
" $$ g(k;a) = u'(\\bar{c}(k)) - \\beta u'[\\bar{c}(f(k)-\\bar{c}(k))] f'(f(k)-\\bar{c}(k))$$\n",
"3. Choose $a$ to make $g(k;a)$ small in some sense"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"What's *small in some sense*?\n",
"\n",
"* Least-squares: minimize sum of squared errors\n",
"$$\n",
"\\min_a \\int g(k;a)^2 dk\n",
"$$\n",
"* Galerkin: zero out weighted averages of Euler errors\n",
"* Collocation: zero out Euler equation errors at grid $k\\in\\{k_1,\\dots,k_n\\}$:\n",
"$$\n",
"P_i(a) \\equiv g(k_i;a) = 0, i=1,\\dots,n\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### General Projection Method\n",
"\n",
"1. Express solution in terms of unknown function\n",
"$$\n",
" \\mathcal{N}(h)=0\n",
" $$\n",
" where $h(x)$ is the equilibrium function at state $x$\n",
"1. Choose a space for appximation\n",
"1. Find $\\bar{h}$ which nearly solves $$ \\mathcal{N}(\\bar{h})=0$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Projection method exercise\n",
"\n",
"* suppose we want to find effective supply of an oligopolistic firm in cournot competition.\n",
"* We want to know $q = S(p)$, how much is supplied at each price $p$.\n",
"* This function is characterized as\n",
"\n",
"$$\n",
"p + \\frac{S(p)}{D'(p)} - MC(S(p)) = 0,\\forall p>0\n",
"$$\n",
"\n",
"* Take $D(p) = p^{-\\eta}$ and $MC(q) = \\alpha \\sqrt{q} + q^2$.\n",
"* Our task is to solve for $S(p)$ in \n",
"\n",
"$$\n",
"p - \\frac{S(p)p^{\\eta+1}}{\\eta} - \\alpha \\sqrt{S(p)} - S(p)^2 = 0,\\forall p>0\n",
"$$\n",
"\n",
"* No closed form solution. But collocation works!\n",
"\n",
"#### TASK\n",
"\n",
"1. solve for $S(p)$ by collocation\n",
"2. Plot residual function\n",
"3. Plot resulting $mS(p)$ together with market demand and $m=1,10,20$ for market size."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Results of Nonlinear Solver Algorithm\n",
" * Algorithm: Trust-region with dogleg and autoscaling\n",
" * Starting Point: [0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3]\n",
" * Zero: [0.248768, 0.0838916, -0.13965, 0.0447411, 0.00701804, -0.0135233, 0.00715223, -0.00229524, 0.000329096, 0.000224355, -0.000290968, 0.000210507, -0.000107876, 3.32945e-5, 4.04856e-6, -1.48411e-5, 1.28265e-5, -7.34842e-6, 2.74595e-6, -1.36709e-7, -8.28018e-7, 8.69542e-7, -5.79429e-7, 2.87391e-7, -1.02456e-7]\n",
" * Inf-norm of residuals: 0.000000\n",
" * Iterations: 9\n",
" * Convergence: true\n",
" * |x - x'| < 0.0e+00: false\n",
" * |f(x)| < 1.0e-08: true\n",
" * Function Calls (f): 8\n",
" * Jacobian Calls (df/dx): 7\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using CompEcon\n",
"using Plots\n",
"using NLsolve\n",
"function proj(n=25)\n",
"\n",
" alpha = 1.0\n",
" eta = 1.5\n",
" a = 0.1\n",
" b = 3.0\n",
" basis = fundefn(:cheb,n,a,b)\n",
" p = funnode(basis)[1] # collocation points\n",
"\n",
" c0 = ones(n)*0.3\n",
" function resid!(c::Vector,result::Vector,p,basis,alpha,eta)\n",
" # your turn!\n",
" q = funeval(c,basis,p)[1]\n",
" q2 = similar(q)\n",
" for i in eachindex(q2)\n",
" if q[i] < 0\n",
" q2[i] = -20.0\n",
" else\n",
" q2[i] = sqrt(q[i])\n",
" end\n",
" end\n",
" result[:] = p.+ q .*((-1/eta)*p.^(eta+1)) .- alpha*q2 .- q.^2\n",
" end\n",
" f_closure(r::Vector,x::Vector) = resid!(x,r,p,basis,alpha,eta)\n",
" res = nlsolve(f_closure,c0)\n",
" println(res)\n",
"\n",
" # plot residual function\n",
" x = collect(range(a, stop = b, length = 501))\n",
" y = similar(x)\n",
" resid!(res.zero,y,x,basis,alpha,eta);\n",
" y = funeval(res.zero,basis,x)[1]\n",
" pl = Any[]\n",
" push!(pl,plot(x,y,title=\"residual function\"))\n",
" \n",
" # plot supply functions at levels 1,10,20\n",
" \n",
" # plot demand function\n",
" y = funeval(res.zero,basis,x)[1]\n",
" p2 = plot(y,x,label=\"supply 1\")\n",
" plot!(10*y,x,label=\"supply 10\")\n",
" plot!(20*y,x,label=\"supply 20\")\n",
" d = x.^(-eta)\n",
" plot!(d,x,label=\"Demand\")\n",
"\n",
" push!(pl,p2)\n",
" \n",
" plot(pl...,layout=2)\n",
" \n",
"end\n",
"proj()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Endogenous Grid Method (EGM)\n",
"\n",
"* Fast, elegant and precise method to solve consumption/savings problems\n",
"* One continuous state variable\n",
"* One continuous control variable\n",
" $$V(M_t) = \\max_{0 introduced this method.\n",
"* The idea is as follows:\n",
" * Instead of using non-linear root finding for optimal $c$ (see above)\n",
" * fix a grid of possible end-of-period asset levels $A_t$\n",
" * use structure of model to find implied beginning of period cash in hand.\n",
" * We use euler equation and envelope condition to connect $M_{t+1}$ with $c_t$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Recall Traditional Methods: VFI and Euler Equation \n",
"\n",
"* Just to be clear, let us repeat what we did in the beginning of this lecture, using the $M_t$ notation.\n",
" $$\n",
" \\begin{aligned}\n",
" V(M_t) &= \\max_{0\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# minimal EGM implementation, go here: https://github.com/floswald/DCEGM.jl/blob/master/src/dc_algo.jl#L4\n",
"# try out: \n",
"# ] dev https://github.com/floswald/DCEGM.jl\n",
"using DCEGM\n",
"DCEGM.minimal_EGM(dplot = true);"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
"\n",
"## Discrete Choice EGM\n",
"\n",
"* This is a method developed by Fedor Iskhakov, Thomas Jorgensen, John Rust and Bertel Schjerning.\n",
"* Reference: [@iskhakovRust2014] \n",
"* Suppose we have several discrete choices (like \"work/retire\"), combined with a continuous choice in each case (like \"how much to consume given work/retire\").\n",
"* Let $d=0$ mean to retire.\n",
"* Write the problem of a worker as\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"V_t(M_t) & = \\max \\left[ v_t(M_t|d_t=0), v_t(M_t|d_t=1) \\right] \\\\\n",
" & \\text{with}\\\\\n",
"v_t(M_t|d_t=0) & = \\max_{0show that there will be a kink point $\\bar{M}$ such that \n",
" $$ v_t(\\bar{M}|d_t=0) = v_t(\\bar{M}|d_t=1) $$\n",
" * We call any such point a **primary kink** (because it refers to a discrete choice in the **current period**)\n",
"* $V$ is not differentiable at $\\bar{M}$.\n",
"* However, it can be shown that both left and right derivatives exist, with\n",
" $$ V^-(\\bar{M}) < V^+(\\bar{M}) $$\n",
"* Given that the value of the derivative changes discretely at $\\bar{M_t}$, the value function in $t-1$ will exhibit a discontinuity as well:\n",
" * $v_{t-1}$ depends on $V_t$.\n",
" * Tracing out the optimal choice of $c_{t-1}$ implies next period cash on hand $M_t$, and as that hits $\\bar{M_t}$, the derivative jumps.\n",
" * The derivative of the value function determines optimal behaviour via the Euler Equation.\n",
" * We call a discontinuity in $v_{t-1}$ arising from a kink in $V_t$ a **secondary kink**.\n",
"* The kinks propagate backwards. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"\n",
"* [@iskhakovRust2014] provide an analytic example where one can compute the actual number of kinks in period 1 of T.\n",
"* Figure 1 in [@clausenenvelope]:\n",
"\n",
"\n",
"\n",
"![[@iskhakovRust2014] figure 1](../assets/figs/fedor-1.png)\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Kinks\n",
"\n",
"* Refer back to the work/retirement model from before.\n",
"* 6 period implementation of the DC-EGM method:\n",
"\n",
"![github/floswald](../assets/figs/Dchoice_condV.png)\n",
"\n",
"* [Iskhakov @ cemmap 2015: Value functions in T-1](http://www.cemmap.ac.uk/event/id/1213)\n",
"\n",
"\n",
"\n",
"* [Iskhakov @ cemmap 2015: Value functions in T-2](http://www.cemmap.ac.uk/event/id/1213)\n",
"\n",
"\n",
"* [Iskhakov @ cemmap 2015: Consumption function in T-2](http://www.cemmap.ac.uk/event/id/1213)\n",
"\n",
"\n",
"* Optimal consumption in 6 period model:\n",
"![github/floswald](../assets/figs/Dchoice_envC.png)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### The Problem with Kinks\n",
"\n",
"* Relying on fast methods that rely on first order conditions (like euler equation) will fail.\n",
"* There are multiple zeros in the Euler Equation, and a standard Euler Equation approach is not guaranteed to find the right one.\n",
"* picture from Fedor Iskhakov's master class at [cemmap 2015](http://www.cemmap.ac.uk/event/id/1213):\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"### DC-EGM Algorithm\n",
"\n",
"1. Do the EGM step for each discrete choice $d$\n",
"1. Compute $d$-specific consumption and value functions\n",
"1. compare $d$-specific value functions to find optimal switch points\n",
"1. Build envelope over $d$-specific consumption functions with knowledge of which optimal $d$ applies where."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### But EGM relies on the Euler Equation?!\n",
"\n",
"* Yes.\n",
"* An important result in [@clausenenvelope] is that the Euler Equation is still the necessary condition for optimal consumption\n",
" * Intuition: marginal utility differs greatly at $\\epsilon+\\bar{M}$. \n",
" * No economic agent would ever locate **at** $\\bar{M}$.\n",
"* This is different from saying that a proceedure that tries to find the zeros of the Euler Equation would still work.\n",
" * this will pick the wrong solution some times.\n",
"* EGM finds **all** solutions. \n",
" * There is a proceedure to discard the \"wrong ones\". Proof in [@iskhakovRust2014]\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"### Adding Shocks\n",
"\n",
"* This problem is hard to solve with standard methods.\n",
"* It is hard, because the only reliable method is VFI, and this is not feasible in large problems.\n",
"* Adding shocks to non-smooth problems is a widely used remedy.\n",
" * think of \"convexifying\" in game theoretic models\n",
" * (Add a lottery)\n",
" * Also used a lot in macro\n",
"* Adding shocks does indeed help in the current model.\n",
" * We add idiosyncratic taste shocks: Type 1 EV.\n",
" * Income uncertainty: \n",
" * In general, the more shocks, the more smoothing.\n",
"* The problem becomes\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"V_t(M_t) & = \\max \\left[ v_t(M_t|d_t=0) + \\sigma_\\epsilon \\epsilon_t(0), v_t(M_t|d_t=1) + \\sigma_\\epsilon \\epsilon_t(1)\\right] \\\\\n",
"v_t(M_t|d_t=1) & = \\max_{0