Code should execute sequentially if run in a Jupyter notebook

- See the set up page to install Jupyter, Julia (0.6+) and all necessary libraries
- Please direct feedback to contact@quantecon.org or the discourse forum

# Optimal Savings III: Occasionally Binding Constraints¶

Contents

## Overview¶

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

- For further details on economic intuition, implication and models, see [LS18]
- Proofs of all mathematical results stated below can be found in
`this paper`

To solve the model we will use Euler equation based time iteration, similar to this lecture

This method turns out to be

- Globally convergent under mild assumptions, even when utility is unbounded (both above and below)
- More efficient numerically than value function iteration

## The Optimal Savings Problem¶

Let’s write down the model and then discuss how to solve it

### Set Up¶

Consider a household that chooses a state-contingent consumption plan \(\{c_t\}_{t \geq 0}\) to maximize

subject to

Here

- \(\beta \in (0,1)\) is the discount factor
- \(a_t\) is asset holdings at time \(t\), with ad-hoc borrowing constraint \(a_t \geq -b\)
- \(c_t\) is consumption
- \(z_t\) is non-capital income (wages, unemployment compensation, etc.)
- \(R := 1 + r\), where \(r > 0\) is the interest rate on savings

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

We further assume that

- \(r > 0\) and \(\beta R < 1\)
- \(u\) is smooth, strictly increasing and strictly concave with \(\lim_{c \to 0} u'(c) = \infty\) and \(\lim_{c \to \infty} u'(c) = 0\)

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

- \((a_0, z_0) = (a, z)\)
- the feasibility constraints in (1), and
- measurability of \(c_t\) w.r.t. the filtration generated by \(\{z_1, \ldots, z_t\}\)

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

### Value Function and Euler Equation¶

The *value function* \(V \colon S \to \mathbb{R}\) is defined by

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

and

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

### Optimality Results¶

Given our assumptions, it is `known`

that

- For each \((a,z) \in S\), a unique optimal consumption path from \((a,z)\) exists
- This path is the unique feasible path from \((a,z)\) satisfying the Euler equality (5) and the transversality condition

Moreover, there exists an *optimal consumption function* \(c^* \colon S \to [0, \infty)\) such that the path from \((a,z)\) generated by

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^*\)

## Computation¶

There are two standard ways to solve for \(c^*\)

- Time iteration (TI) using the Euler equality
- Value function iteration (VFI)

Let’s look at these in turn

### Time Iteration¶

We can rewrite (5) to make it a statement about functions rather than random variables

In particular, consider the functional equation

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

- each \(c \in \mathscr{C}\) is continuous and (weakly) increasing
- \(\min Z \leq c(a,z) \leq Ra + z + b\) for all \((a,z) \in S\)

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

where

We refer to \(K\) as Coleman’s policy function operator [Col90]

It is `known`

that

- \(K\) is a contraction mapping on \(\mathscr{C}\) under the metric

- The metric \(\rho\) is complete on \(\mathscr{C}\)
- Convergence in \(\rho\) implies uniform convergence on compacts

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)

### Value Function Iteration¶

The Bellman operator for this problem is given by

We have to be careful with VFI (i.e., iterating with \(T\)) in this setting because \(u\) is not assumed to be bounded

- In fact typically unbounded both above and below — e.g. \(u(c) = \log c\)
- In which case, the standard DP theory does not apply
- \(T^n v\) is not guaranteed to converge to the value function for arbitrary continous bounded \(v\)

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

### Implementation¶

Here’s the code for a type called `ConsumerProblem`

that stores primitives, as well as

- a
`bellman_operator`

function, which implements the Bellman operator \(T\) specified above - a
`coleman_operator`

function, which implements the Coleman operator \(K\) specified above - an
`initialize`

, which generates suitable initial conditions for iteration

```
using QuantEcon
using Optim
# utility and marginal utility functions
u(x) = log(x)
du(x) = 1 / x
"""
Income fluctuation problem
##### Fields
- `r::Real` : Strictly positive interest rate
- `R::Real` : The interest rate plus 1 (strictly greater than 1)
- `β::Real` : Discount rate in (0, 1)
- `b::Real` : The borrowing constraint
- `Π::Matrix{T} where T<:Real` : Transition matrix for `z`
- `z_vals::Vector{T} where T<:Real` : Levels of productivity
- `asset_grid::AbstractArray` : Grid of asset values
"""
struct ConsumerProblem{T <: Real}
r::T
R::T
β::T
b::T
Π::Matrix{T}
z_vals::Vector{T}
asset_grid::AbstractArray
end
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 = linspace(-b, grid_max, grid_size)
ConsumerProblem(r, R, β, b, Π, z_vals, asset_grid)
end
"""
Apply the Bellman operator for a given model and initial value.
##### Arguments
- `cp::ConsumerProblem` : Instance of `ConsumerProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output
- `;ret_policy::Bool(false)`: Toggles return of value or policy functions
##### Returns
None, `out` is updated in place. If `ret_policy == true` out is filled with the
policy function, otherwise the value function is stored in `out`.
"""
function bellman_operator!(cp::ConsumerProblem,
V::Matrix,
out::Matrix;
ret_policy::Bool=false)
# simplify names, set up arrays
R, Π, β, b = cp.R, cp.Π, cp.β, cp.b
asset_grid, z_vals = cp.asset_grid, cp.z_vals
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 = optimize(obj, opt_lb, R .* a .+ z .+ b)
c_star = Optim.minimizer(res)
if ret_policy
out[i_a, i_z] = c_star
else
out[i_a, i_z] = - Optim.minimum(res)
end
end
end
out
end
bellman_operator(cp::ConsumerProblem, V::Matrix; ret_policy=false) =
bellman_operator!(cp, V, similar(V); ret_policy=ret_policy)
"""
Extract the greedy policy (policy function) of the model.
##### Arguments
- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output
##### Returns
None, `out` is updated in place to hold the policy function
"""
get_greedy!(cp::ConsumerProblem, V::Matrix, out::Matrix) =
update_bellman!(cp, V, out, ret_policy=true)
get_greedy(cp::ConsumerProblem, V::Matrix) =
update_bellman(cp, V, ret_policy=true)
"""
The approximate Coleman operator.
Iteration with this operator corresponds to policy function
iteration. Computes and returns the updated consumption policy
c. The array c is replaced with a function cf that implements
univariate linear interpolation over the asset grid for each
possible value of z.
##### Arguments
- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `c::Matrix`: Current guess for the policy function
- `out::Matrix` : Storage for output
##### Returns
None, `out` is updated in place to hold the policy function
"""
function coleman_operator!(cp::ConsumerProblem, c::Matrix, out::Matrix)
# simplify names, set up arrays
R, Π, β, b = cp.R, cp.Π, cp.β, cp.b
asset_grid, z_vals = cp.asset_grid, cp.z_vals
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] = Optim.minimizer(res)
end
end
out
end
"""
Apply the Coleman operator for a given model and initial value
See the specific methods of the mutating version of this function for more
details on arguments
"""
coleman_operator(cp::ConsumerProblem, c::Matrix) =
coleman_operator!(cp, c, similar(c))
function initialize(cp::ConsumerProblem)
# simplify names, set up arrays
R, β, b = cp.R, cp.β, cp.b
asset_grid, z_vals = cp.asset_grid, cp.z_vals
shape = length(asset_grid), length(z_vals)
V, c = Array{Float64}(shape...), Array{Float64}(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
```

Both `bellman_operator`

and `coleman_operator`

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

## Exercises¶

### Exercise 1¶

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

- In the case of iteration with \(T\), the final value function is used to compute the observed policy

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

- The parameters are the default parameters in the definition of
`consumerProblem`

- The initial conditions are the default ones from
`initialize(cp)`

- Both operators are iterated 80 times

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)
```

```
@time bellman_operator(cp, v);
```

```
elapsed time: 0.095017748 seconds (24212168 bytes allocated, 30.48% gc time)
```

```
@time coleman_operator(cp, c);
```

```
elapsed time: 0.0696242 seconds (23937576 bytes allocated)
```

### Exercise 2¶

Next let’s consider how the interest rate affects consumption

Reproduce the following figure, which shows (approximately) optimal consumption policies for different interest rates

- Other than r, all parameters are at their default values
- r steps through linspace(0, 0.04, 4)
- Consumption is plotted against assets for income shock fixed at the smallest value

The figure shows that higher interest rates boost savings and hence suppress consumption

### Exercise 3¶

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

```
using Plots
pyplot()
# === solve for optimal consumption === #
m = ConsumerProblem(r=0.03, grid_max=4)
v_init, c_init = initialize(m)
c = compute_fixed_point(c -> coleman_operator(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")
```

The blue line and orange line represent the function

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

- Can be proved via theorem 2 of [HP92]
- Represents the long run dispersion of assets across households when households have idiosyncratic shocks

Ergodicity is valid here, so stationary probabilities can be calculated by averaging over a single long time series

- Hence to approximate the stationary distribution we can simulate a long time series for assets and histogram, as in the following figure

Your task is to replicate the figure

- Parameters are as discussed above
- The histogram in the figure used a single time series \(\{a_t\}\) of length 500,000
- Given the length of this time series, the initial condition \((a_0, z_0)\) will not matter
- You might find it helpful to use the
`MarkovChain`

type from`quantecon`

### Exercise 4¶

Following on from exercises 2 and 3, let’s look at how savings and aggregate asset holdings vary with the interest rate

- Note: [LS18] section 18.6 can be consulted for more background on the topic treated in this exercise

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

## Solutions¶

### Exercise 1¶

```
cp = ConsumerProblem()
K = 80
V, c = initialize(cp)
println("Starting value function iteration")
for i=1:K
V = bellman_operator(cp, V)
end
c1 = bellman_operator(cp, V, ret_policy=true)
V2, c2 = initialize(cp)
println("Starting policy function iteration")
for i=1:K
c2 = coleman_operator(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)")
```

```
Starting value function iteration
Starting policy function iteration
```

### Exercise 2¶

```
r_vals = linspace(0, 0.04, 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 -> coleman_operator(cp, x),
c_init,
max_iter=150,
verbose=false)
traces = push!(traces, c[:, 1])
legends = push!(legends, "r = $(round(r_val, 3))")
end
plot(traces, label=reshape(legends, 1, length(legends)))
plot!(xlabel="asset level", ylabel="Consumption (low income)")
```

### Exercise 3¶

```
function compute_asset_series(cp, T=500000; verbose=false)
Π, z_vals, R = cp.Π, cp.z_vals, cp.R # Simplify names
z_idx = 1:length(z_vals)
v_init, c_init = initialize(cp)
c = compute_fixed_point(x -> coleman_operator(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=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)
a = compute_asset_series(cp)
histogram(a, nbins=20, leg=false, normed=true, xlabel="assets")
```

### Exercise 4¶

```
M = 25
r_vals = linspace(0, 0.04, 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, 250000))
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]))
```

```
Finished iteration b=1.0
Finished iteration b=3.0
```