Next we study an optimal savings problem for an infinitely lived consumer—the “common ancestor” described in [LS18], section 1.3.
This is an essential sub-problem for many representative macroeconomic models
It is related to the decision problem in the stochastic optimal growth model and yet differs in important ways.
For example, the choice problem for the agent includes an additive income term that leads to an occasionally binding constraint.
Our presentation of the model will be relatively brief.
To solve the model we will use Euler equation based time iteration, similar to this lecture.
This method turns out to be
Let’s write down the model and then discuss how to solve it.
Consider a household that chooses a state-contingent consumption plan $ \{c_t\}_{t \geq 0} $ to maximize
$$ \mathbb{E} \, \sum_{t=0}^{\infty} \beta^t u(c_t) $$subject to
$$ c_t + a_{t+1} \leq Ra_t + z_t, \qquad c_t \geq 0, \qquad a_t \geq -b \qquad t = 0, 1, \ldots \tag{1} $$
Here
Non-capital income $ \{z_t\} $ is assumed to be a Markov process taking values in $ Z\subset (0,\infty) $ with stochastic kernel $ \Pi $.
This means that $ \Pi(z, B) $ is the probability that $ z_{t+1} \in B $ given $ z_t = z $.
The expectation of $ f(z_{t+1}) $ given $ z_t = z $ is written as
$$ \int f( \acute z) \, \Pi(z, d \acute z) $$We further assume that
The asset space is $ [-b, \infty) $ and the state is the pair $ (a,z) \in S := [-b,\infty) \times Z $.
A feasible consumption path from $ (a,z) \in S $ is a consumption sequence $ \{c_t\} $ such that $ \{c_t\} $ and its induced asset path $ \{a_t\} $ satisfy
The meaning of the third point is just that consumption at time $ t $ can only be a function of outcomes that have already been observed.
The value function $ V \colon S \to \mathbb{R} $ is defined by
$$ V(a, z) := \sup \, \mathbb{E} \left\{ \sum_{t=0}^{\infty} \beta^t u(c_t) \right\} \tag{2} $$
where the supremum is over all feasible consumption paths from $ (a,z) $.
An optimal consumption path from $ (a,z) $ is a feasible consumption path from $ (a,z) $ that attains the supremum in (2).
To pin down such paths we can use a version of the Euler equation, which in the present setting is
$$ u' (c_t) \geq \beta R \, \mathbb{E}_t [ u'(c_{t+1}) ] \tag{3} $$
and
$$ u' (c_t) = \beta R \, \mathbb{E}_t [ u'(c_{t+1}) ] \quad \text{whenever } c_t < Ra_t + z_t + b \tag{4} $$
In essence, this says that the natural “arbitrage” relation $ u' (c_t) = \beta R \, \mathbb{E}_t [ u'(c_{t+1}) ] $ holds when the choice of current consumption is interior.
Interiority means that $ c_t $ is strictly less than its upper bound $ Ra_t + z_t + b $.
(The lower boundary case $ c_t = 0 $ never arises at the optimum because $ u'(0) = \infty $)
When $ c_t $ does hit the upper bound $ Ra_t + z_t + b $, the strict inequality $ u' (c_t) > \beta R \, \mathbb{E}_t [ u'(c_{t+1}) ] $ can occur because $ c_t $ cannot increase sufficiently to attain equality.
With some thought and effort, one can show that (3) and (4) are equivalent to
$$ u' (c_t) = \max \left\{ \beta R \, \mathbb{E}_t [ u'(c_{t+1}) ] \,,\; u'(Ra_t + z_t + b) \right\} \tag{5} $$
Given our assumptions, it is known that
Euler equality (5) and the transversality condition
$$ \lim_{t \to \infty} \beta^t \, \mathbb{E} \, [ u'(c_t) a_{t+1} ] = 0. \tag{6} $$
Moreover, there exists an optimal consumption function $ c^* \colon S \to [0, \infty) $ such that the path from $ (a,z) $ generated by
$$ (a_0, z_0) = (a, z), \quad z_{t+1} \sim \Pi(z_t, dy), \quad c_t = c^*(a_t, z_t) \quad \text{and} \quad a_{t+1} = R a_t + z_t - c_t $$satisfies both (5) and (6), and hence is the unique optimal path from $ (a,z) $.
In summary, to solve the optimization problem, we need to compute $ c^* $.
There are two standard ways to solve for $ c^* $
Let’s look at these in turn.
We can rewrite (5) to make it a statement about functions rather than random variables.
In particular, consider the functional equation
$$ u' \circ c \, (a, z) = \max \left\{ \gamma \int u' \circ c \, \{R a + z - c(a, z), \, \acute z\} \, \Pi(z,d \acute z) \, , \; u'(Ra + z + b) \right\} \tag{7} $$
where $ \gamma := \beta R $ and $ u' \circ c(s) := u'(c(s)) $.
Equation (7) is a functional equation in $ c $.
In order to identify a solution, let $ \mathscr{C} $ be the set of candidate consumption functions $ c \colon S \to \mathbb R $ such that
In addition, let $ K \colon \mathscr{C} \to \mathscr{C} $ be defined as follows:
For given $ c\in \mathscr{C} $, the value $ Kc(a,z) $ is the unique $ t \in J(a,z) $ that solves
$$ u'(t) = \max \left\{ \gamma \int u' \circ c \, \{R a + z - t, \, \acute z\} \, \Pi(z,d \acute z) \, , \; u'(Ra + z + b) \right\} \tag{8} $$
where
$$ J(a,z) := \{t \in \mathbb{R} \,:\, \min Z \leq t \leq Ra+ z + b\} \tag{9} $$
We refer to $ K $ as Coleman’s policy function operator [Col90].
It is known that
In consequence, $ K $ has a unique fixed point $ c^* \in \mathscr{C} $ and $ K^n c \to c^* $ as $ n \to \infty $ for any $ c \in \mathscr{C} $.
By the definition of $ K $, the fixed points of $ K $ in $ \mathscr{C} $ coincide with the solutions to (7) in $ \mathscr{C} $.
In particular, it can be shown that the path $ \{c_t\} $ generated from $ (a_0,z_0) \in S $ using policy function $ c^* $ is the unique optimal path from $ (a_0,z_0) \in S $.
TL;DR The unique optimal policy can be computed by picking any $ c \in \mathscr{C} $ and iterating with the operator $ K $ defined in (8).
The Bellman operator for this problem is given by
$$ Tv(a, z) = \max_{0 \leq c \leq Ra + z + b} \left\{ u(c) + \beta \int v(Ra + z - c, \acute z) \Pi(z, d \acute z) \right\} \tag{10} $$
We have to be careful with VFI (i.e., iterating with $ T $) in this setting because $ u $ is not assumed to be bounded
Nonetheless, we can always try the popular strategy “iterate and hope”.
We can then check the outcome by comparing with that produced by TI.
The latter is known to converge, as described above.
Here’s the code for a named-tuple constructor called ConsumerProblem
that stores primitives, as well as
T
function, which implements the Bellman operator $ T $ specified aboveK
function, which implements the Coleman operator $ K $ specified aboveinitialize
, which generates suitable initial conditions for iterationusing InstantiateFromURL
# optionally add arguments to force installation: instantiate = true, precompile = true
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.8.0")
using LinearAlgebra, Statistics
using BenchmarkTools, Optim, Parameters, Plots, QuantEcon, Random
using Optim: converged, maximum, maximizer, minimizer, iterations
gr(fmt = :png);
# utility and marginal utility functions
u(x) = log(x)
du(x) = 1 / x
# model
function ConsumerProblem(;r = 0.01,
β = 0.96,
Π = [0.6 0.4; 0.05 0.95],
z_vals = [0.5, 1.0],
b = 0.0,
grid_max = 16,
grid_size = 50)
R = 1 + r
asset_grid = range(-b, grid_max, length = grid_size)
return (r = r, R = R, β = β, b = b, Π = Π, z_vals = z_vals, asset_grid = asset_grid)
end
function T!(cp, V, out; ret_policy = false)
# unpack input, set up arrays
@unpack R, Π, β, b, asset_grid, z_vals = cp
z_idx = 1:length(z_vals)
# value function when the shock index is z_i
vf = interp(asset_grid, V)
opt_lb = 1e-8
# solve for RHS of Bellman equation
for (i_z, z) in enumerate(z_vals)
for (i_a, a) in enumerate(asset_grid)
function obj(c)
EV = dot(vf.(R * a + z - c, z_idx), Π[i_z, :]) # compute expectation
return u(c) + β * EV
end
res = maximize(obj, opt_lb, R .* a .+ z .+ b)
converged(res) || error("Didn't converge") # important to check
if ret_policy
out[i_a, i_z] = maximizer(res)
else
out[i_a, i_z] = maximum(res)
end
end
end
out
end
T(cp, V; ret_policy = false) =
T!(cp, V, similar(V); ret_policy = ret_policy)
get_greedy!(cp, V, out) =
update_bellman!(cp, V, out, ret_policy = true)
get_greedy(cp, V) =
update_bellman(cp, V, ret_policy = true)
function K!(cp, c, out)
# simplify names, set up arrays
@unpack R, Π, β, b, asset_grid, z_vals = cp
z_idx = 1:length(z_vals)
gam = R * β
# policy function when the shock index is z_i
cf = interp(asset_grid, c)
# compute lower_bound for optimization
opt_lb = 1e-8
for (i_z, z) in enumerate(z_vals)
for (i_a, a) in enumerate(asset_grid)
function h(t)
cps = cf.(R * a + z - t, z_idx) # c' for each z'
expectation = dot(du.(cps), Π[i_z, :])
return abs(du(t) - max(gam * expectation, du(R * a + z + b)))
end
opt_ub = R*a + z + b # addresses issue #8 on github
res = optimize(h, min(opt_lb, opt_ub - 1e-2), opt_ub,
method = Optim.Brent())
out[i_a, i_z] = minimizer(res)
end
end
return out
end
K(cp, c) = K!(cp, c, similar(c))
function initialize(cp)
# simplify names, set up arrays
@unpack R, β, b, asset_grid, z_vals = cp
shape = length(asset_grid), length(z_vals)
V, c = zeros(shape...), zeros(shape...)
# populate V and c
for (i_z, z) in enumerate(z_vals)
for (i_a, a) in enumerate(asset_grid)
c_max = R * a + z + b
c[i_a, i_z] = c_max
V[i_a, i_z] = u(c_max) / (1 - β)
end
end
return V, c
end
initialize (generic function with 1 method)
Both T
and K
use linear interpolation along the asset grid to approximate the value and consumption functions.
The following exercises walk you through several applications where policy functions are computed.
In exercise 1 you will see that while VFI and TI produce similar results, the latter is much faster.
Intuition behind this fact was provided in a previous lecture on time iteration.
The first exercise is to replicate the following figure, which compares TI and VFI as solution methods
The figure shows consumption policies computed by iteration of $ K $ and $ T $ respectively
Consumption is shown as a function of assets with income $ z $ held fixed at its smallest value.
The following details are needed to replicate the figure
consumerProblem
.initialize(cp)
.When you run your code you will observe that iteration with $ K $ is faster than iteration with $ T $.
In the Julia console, a comparison of the operators can be made as follows
cp = ConsumerProblem()
v, c, = initialize(cp)
([-17.328679513998615 0.0; -4.664387247760648 7.125637140580436; … ; 69.82541010621486 70.57937835479346; 70.32526591846735 71.06452735149534], [0.5 1.0; 0.8297959183673469 1.329795918367347; … ; 16.33020408163265 16.83020408163265; 16.66 17.16])
@btime T(cp, v);
492.517 μs (4676 allocations: 428.20 KiB)
@btime K(cp, c);
695.109 μs (7925 allocations: 742.17 KiB)
Next let’s consider how the interest rate affects consumption.
Reproduce the following figure, which shows (approximately) optimal consumption policies for different interest rates
The figure shows that higher interest rates boost savings and hence suppress consumption.
Now let’s consider the long run asset levels held by households.
We’ll take r = 0.03 and otherwise use default parameters.
The following figure is a 45 degree diagram showing the law of motion for assets when consumption is optimal
# solve for optimal consumption
m = ConsumerProblem(r = 0.03, grid_max = 4)
v_init, c_init = initialize(m)
c = compute_fixed_point(c -> K(m, c),
c_init,
max_iter = 150,
verbose = false)
a = m.asset_grid
R, z_vals = m.R, m.z_vals
# generate savings plot
plot(a, R * a .+ z_vals[1] - c[:, 1], label = "Low income")
plot!(xlabel = "Current assets", ylabel = "Next period assets")
plot!(a, R * a .+ z_vals[2] - c[:, 2], label = "High income")
plot!(xlabel = "Current assets", ylabel = "Next period assets")
plot!(a, a, linestyle = :dash, color = "black", label = "")
plot!(xlabel = "Current assets", ylabel = "Next period assets")
plot!(legend = :topleft)
The blue line and orange line represent the function
$$ a' = h(a, z) := R a + z - c^*(a, z) $$when income $ z $ takes its high and low values respectively.
The dashed line is the 45 degree line.
We can see from the figure that the dynamics will be stable — assets do not diverge.
In fact there is a unique stationary distribution of assets that we can calculate by simulation
Ergodicity is valid here, so stationary probabilities can be calculated by averaging over a single long time series.
Your task is to replicate the figure
MarkovChain
type from quantecon
Following on from exercises 2 and 3, let’s look at how savings and aggregate asset holdings vary with the interest rate
For a given parameterization of the model, the mean of the stationary distribution can be interpreted as aggregate capital in an economy with a unit mass of ex-ante identical households facing idiosyncratic shocks.
Let’s look at how this measure of aggregate capital varies with the interest rate and borrowing constraint.
The next figure plots aggregate capital against the interest rate for b in (1, 3)
As is traditional, the price (interest rate) is on the vertical axis.
The horizontal axis is aggregate capital computed as the mean of the stationary distribution.
Exercise 4 is to replicate the figure, making use of code from previous exercises.
Try to explain why the measure of aggregate capital is equal to $ -b $ when $ r=0 $ for both cases shown here.
cp = ConsumerProblem()
N = 80
V, c = initialize(cp)
println("Starting value function iteration")
for i in 1:N
V = T(cp, V)
end
c1 = T(cp, V, ret_policy=true)
V2, c2 = initialize(cp)
println("Starting policy function iteration")
for i in 1:N
c2 = K(cp, c2)
end
plot(cp.asset_grid, c1[:, 1], label = "value function iteration")
plot!(cp.asset_grid, c2[:, 1], label = "policy function iteration")
plot!(xlabel = "asset level", ylabel = "Consumption (low income)")
plot!(legend = :topleft)
Starting value function iteration Starting policy function iteration
r_vals = range(0, 0.04, length = 4)
traces = []
legends = []
for r_val in r_vals
cp = ConsumerProblem(r = r_val)
v_init, c_init = initialize(cp)
c = compute_fixed_point(x -> K(cp, x),
c_init,
max_iter = 150,
verbose = false)
traces = push!(traces, c[:, 1])
legends = push!(legends, "r = $(round(r_val, digits = 3))")
end
plot(traces, label = reshape(legends, 1, length(legends)))
plot!(xlabel = "asset level", ylabel = "Consumption (low income)")
plot!(legend = :topleft)
function compute_asset_series(cp, T = 500_000; verbose = false)
@unpack Π, z_vals, R = cp # simplify names
z_idx = 1:length(z_vals)
v_init, c_init = initialize(cp)
c = compute_fixed_point(x -> K(cp, x), c_init,
max_iter = 150, verbose = false)
cf = interp(cp.asset_grid, c)
a = zeros(T + 1)
z_seq = simulate(MarkovChain(Π), T)
for t in 1:T
i_z = z_seq[t]
a[t+1] = R * a[t] + z_vals[i_z] - cf(a[t], i_z)
end
return a
end
cp = ConsumerProblem(r = 0.03, grid_max = 4)
Random.seed!(42) # for reproducibility
a = compute_asset_series(cp)
histogram(a, nbins = 20, leg = false, normed = true, xlabel = "assets")
M = 25
r_vals = range(0, 0.04, length = M)
xs = []
ys = []
legends = []
for b in [1.0, 3.0]
asset_mean = zeros(M)
for (i, r_val) in enumerate(r_vals)
cp = ConsumerProblem(r = r_val, b = b)
the_mean = mean(compute_asset_series(cp, 250_000))
asset_mean[i] = the_mean
end
xs = push!(xs, asset_mean)
ys = push!(ys, r_vals)
legends = push!(legends, "b = $b")
println("Finished iteration b = $b")
end
plot(xs, ys, label = reshape(legends, 1, length(legends)))
plot!(xlabel = "capital", ylabel = "interest rate", yticks = ([0, 0.045]))
plot!(legend = :bottomright)
Finished iteration b = 1.0 Finished iteration b = 3.0