using LinearAlgebra
using PlotsLecture 4: The Aiyagari Model
In this lecture, I will introduce the basic set-up of the workhorse incomplete market model – the Aiyagari (1996) model and discuss how to compute it numerically. When writing this lecture, I largely follow from my answer for problem set 3 and 4 for Manuel Amador’s Macro class.
Useful notes: Jesus Fernandez-Villaverde’s notes on heterogeneous agent models: part 1, part 2, part 3.
1.General set up of the Aiyagari model
1.1 Environment:
Household:
Continuum of household with mass 1:
Household has identical preferences: \[ \sum_{t=0}^\infty \beta^t u(c_t) \],
subject to a budget constraint: \(c_t + a_{t+1} \leq w_t l_t(s^t) + R_t a_t, \forall t\),
and a borrowing limit: \(a_{t + 1} \geq \phi, \forall t\),
with \(a_0, l_0 \geq 0\) given.
Let labor endowments for each household follow a markov chain with a transition matrix \(\pi\), \(\pi(s'|s) > 0, \forall s, s'\). The invariant distribution of labor supply is a probability distribution \(\lambda\) such that \(\lambda \pi = \pi\).
By law of large number, the aggregate labor supply \(L = \sum_{s \in S} \lambda(s) l(s)\) does not depend on \(s\) (no aggregate uncertainty). This is a key feature of the problem.
Firm:
A perfectly competitive firm with Neo-classical technology: \[ Y_t = F(K_t, L_t). \]
Aggregate resource constraint:
\[ C_t + K_{t+1} - (1 - \delta) K_t = F(K_t, L) \] with depreciation rate \(\delta \in (0,1)\).
2. Recursive formulation of the household problem
Assume utility is strictly increasing. Taking the interest rate \(R\) and \(w\) as given, the household problem takes the recursive form:
\(V(a, s) = \max_{a'\geq \phi} u(Ra + w l(s) - a') + \beta E[V(a', s')]\)
where \(s\) is the realization of the labor endowment today.
2.1 Side notes on alternative specification of the household problem:
Alternatively we can write the problem as \(V(x, s) = \max_{a'\geq \phi} u(x - a') + \beta E[V(Ra' + w l(s'), s')]\), where \(x \equiv Ra + w l(s)\) is the cash in hand today. Note that the value function will be different as it is now a function of \(x\). The key here is that one need to make sure that the state variables on both sides are consistent.
On top of this, Aiyagari used a trick to simplify the problem. Define \(\hat{a} = a + \phi\), \(\hat{x} = x + \phi\). Then, the household problem can be rewritten as:
\(v(\hat{x}, s) = \max_{\hat{a}' \geq 0} u(\hat{x} - \hat{a}') + \beta E[v(R\hat{a}' + w l(s') - (R - 1) \phi), s']\)
Now let’s try to solve the household problem using \(a\) as the first state variable given particular \(R\) and \(w\). This is essentially a partial equilibrium and is often referred to as a Bewley/Hugget model.
2.2 Computing the household problem
Note that here only the product of \(w\) and \(l\) matters here and there’s no interesting mechanism for determining \(w\). For simplicity, we define \(y(s) \equiv wl(s)\) such that \(V(a, s) = \max_{a'\geq \phi} u(Ra + y(s) - a') + \beta E[V(a', s')]\). Further assume \(s\) follows a first-order markov chain governed by states \(y(s)\) and transition matrix \(P\).
We will construct grids using \(a\). Let \(a\) start from \(\phi\) so the agent cannot over-borrow.
The following codes are adopted from problem set 2 and 3 from Manuel Amador’s Macro class.
As before, Base.@kwdef struct is used to store the parameters. However here instead of some specific types like Float64, we use parametric types like R1 and T1. These can really be any types dependent on the input. Check the documentation for parametric types here.
Base.@kwdef struct Household{T1, T2, T3, R1, S}
β::R1 = 0.95
ρ::R1 = 2.0
ϕ::R1 = 0.0 # The borrowing limit is zero
P::T1 = [0.5 0.5; 0.2 0.8] # transition matrix for state s
l::T2 = [0.5, 1.0] # labor endowment: l(s) (normalize w to 1.0 so that y(s) = l(s))
a_min::R1 = ϕ
a_max::R1 = 5.0
points::S = 10_000
a_grid::T3 = range(a_min, a_max, points)
end
u(c, m) = c ^ (1 - m.ρ) / (1 - m.ρ)
uprime(c, m) = c ^ (-m.ρ) uprime (generic function with 1 method)Now solve the value function as before. Since we have two states now, we need a matrix to store the value and policy for each different combination of \(a\) and \(y\).
There are mainly two tricks in the code. The first is to pre-compute the expected discounted value for each \(a'\) and \(s'\) conditional on \(s\) within the outer loop. The second is to take use of the fact of strictly increasing policy function as we have seen in Lecture 2.
function solve_household(h; R = 0.9, w = 1.0, v0 = zeros(length(h.a_grid), length(h.l)), tol = 1e-6)
# unpacking
(; a_grid, l, β, P) = h
v1 = similar(v0)
pol = similar(v0, Int)
βv = similar(v1)
n = length(a_grid)
iter = 0
while true
distance = zero(eltype(v0))
iter += 1
# Precompute βE[V(a', s')] for better performance
for s in eachindex(l)
for i in eachindex(a_grid)
accum = zero(eltype(v0))
for sprime in eachindex(l)
accum += β * P[s, sprime] * v0[i, sprime]
end
βv[i, s] = accum
end
end
# or simply βv = β .* v0 * P'
for s in eachindex(l)
pol_i = 1
for i in eachindex(a_grid)
just_started = true
vmax = zero(eltype(v0))
for j in pol_i:n #take use of strictly increasing policy function
c = R * a_grid[i] + w*l[s] - a_grid[j]
if c > 0.0
v_tmp = u(c, h) + βv[j, s]
if just_started
vmax = v_tmp
pol_i = j
just_started = false
elseif v_tmp > vmax
vmax = v_tmp
pol_i = j
else
break
end
end
end
v1[i, s] = vmax
pol[i, s] = pol_i
dis = abs(vmax - v0[i, s])
if dis > distance
distance = dis
end
end
end
if distance < tol
break
else
v0 .= v1
end
end
return (v = v1, pol = pol, h = h, R = R, w = w)
end
solve_household (generic function with 1 method)Explanations over how the code works (GPT generated, which is actually pretty good):
This function solve_household solves the household problem defined earlier using value function iteration. The function takes in the model parameters as input, par, and some optional arguments such as an initial guess for the value function v0 and a tolerance level tol.
The function initializes some variables such as v1 and pol as arrays of zeros with the same shape as v0, which is the value function from the previous iteration. βv is also initialized as an array of zeros, which will be used to store the expected value of the next iteration.
The function then performs the value function iteration using a nested loop over ytilde and ahat_grid. The inner loop first computes the expected value of the next iteration using the Bellman equation and stores it in βv.
The inner loop then maximizes the Bellman equation over the choice of next period asset holding ahat_grid[j] using a nested loop that starts at pol_i. This is done to take advantage of the fact that the optimal choice of ahat_grid[j] is increasing in j. The loop computes the value of the Bellman equation for each ahat_grid[j] and stops once the first ahat_grid[j] is found that yields a negative consumption level c. The function then stores the optimal value function vmax and the corresponding choice of j in pol_i and moves on to the next ahat_grid[i].
The function then computes the distance between the value function from the current iteration v0 and the value function from the previous iteration v1. If the distance is less than the tolerance level tol, the function exits the loop and returns the value function v1, the policy function pol, and the model parameters m.
Overall, the function uses a nested loop to solve the household problem using value function iteration, and it takes advantage of the fact that the problem has a recursive structure and a simple functional form for the utility function to compute the optimal value and policy functions efficiently.
h = Household() # struct with default parameter values
sol_1 = solve_household(h, R = 0.658) #input some random R. (v = [-26.71326843693425 -25.314667038332857; -26.711953170852723 -25.314338113656607; … ; -24.019524603316466 -23.72455416253825; -24.019409824729976 -23.72445885993552], pol = [1 1; 1 1; … ; 4194 4864; 4194 4864], h = Household{Matrix{Float64}, Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Int64}(0.95, 2.0, 0.0, [0.5 0.5; 0.2 0.8], [0.5, 1.0], 0.0, 5.0, 10000, 0.0:0.0005000500050005:5.0), R = 0.658, w = 1.0)Let’s take a look at the value function and policy functions
function do_v_plot(sol_1)
income = sol_1.w * sol_1.h.l
p1 = Plots.plot(xlabel = "Current saving, a", ylabel = "Value, v")
for i in eachindex(income)
temp = income[i]
plot!(p1, sol_1.h.a_grid, sol_1.v[:,i], label = "y = $temp")
end
return p1
end
do_v_plot(sol_1)function do_pol_plot(sol_1)
(; pol, h)= sol_1
p1 = Plots.plot(xlabel = "Current saving, a", ylabel = "Saving next period, a' ")
for s in eachindex(h.l)
temp = sol_1.w * h.l[s]
yvals = [[h.a_grid[pol[i, s]] for i in eachindex(h.a_grid)]]
plot!(p1, sol_1.h.a_grid, yvals, label = "y = $temp")
end
return p1
end
do_pol_plot(sol_1)Note that when \(y\) is low. The policy function for \(a'\) is flat. The intuition is that the agent has incentive to borrow for consumption smoothing when endowment is low (why?) but is borrowing constrained due to the \(0\) borrowing limit. As a result, the agent consumes all her endowment. This is sometimes referred to as “Hand-to-mouth” consumption.
3. Recursive competitive equilibrium (Aiyagari 1994)
With some understanding of the partial equilibrium, we now try to build towards the general equilibrium. With household optimality, firm optimality and market clearing conditions, the equilibrium can be formally defined by a recursive competitive equilibrium.
The key difference is that every household now needs to know the distribution of everyone else’s asset holdings to infer future interest rate.
Let \(a \in A \equiv [\phi, \infty]\), \(s \in S\). From the recursive form of the household problem, for given R we can solve for a value function \(v(a, s): A \times S \to R\) and a policy function for saving \(g(a, s): A \times S \to A\).
3.1 Prices in general equilibrium
Now for simplicity, let the firm technology be characterized by a CRS Cobb-Douglas production function, i.e. \[ Y_t = A K_t^\alpha L_t^{1 - \alpha}. \]
The firm’s problem takes the form: (why?)
\[\max_{K_t, L_t} A K_t^\alpha L_t^{1 - \alpha} - (R_t -1 +\delta) K_t - w_t L_t\]
Firm’s first order conditions: \[ \begin{cases} R_t = A \alpha K_t^{\alpha - 1} L_t^{1 - \alpha} +1 - \delta \\ w_t = A (1 - \alpha) K_t^\alpha L_t^{-\alpha} \end{cases} \]
Note that \(w\) can be pinned down by \(R\) in equilibrium: \[ w_t(R_t) = A^{\frac{1}{1-\alpha}} (1 - \alpha) \alpha^{\frac{\alpha}{1 - \alpha}} \left(R_t - 1 +\delta \right)^{\frac{\alpha}{1 - \alpha}} \]
This is very handy when it comes to computing the equilibrium prices.
3.2 Set up of a (recursive) stationary competitive equilibrium
A stationary competitive equilibrium is defined by a system of constant prices (interest rate and wage) and allocations such that individuals optimize and markets clear.
Define \(\lambda_t(a, s): A \times S \to [0,1]\) to be the distribution of household over current asset \(a\) and state \(s\). Imagine it to be a matrix with each entry filled with the fraction of household with each combination of possible \(a\) and \(s\). Hypothetically, this should evolve over time following a law of motion: \[ \lambda_{t+1}(a', s') = \sum_{s \in S} \sum_{a: g(a,s) = a'} \pi(s'|s) \lambda_t(a,s), \quad \forall a', s' \]
where \(\pi(s'|s)\) is the transition probability and \(\sum_{a: g(a,s) = a'}\) reads that all \(a\) such that \(g(a,s) = a'\). Intuitively, the fraction of households that falls into the state \(a'\) and \(s'\) next period will be the fraction of households (with different \(a\) and \(s\)) that choose \(a'\) times the conditional probability that the household moves into \(s'\) next period.
Does the economy converge to some kind of steady state (stationary equilibrium)? Here, in a stationary equilibrium, \(\lambda_t = \lambda, \forall t\) and all aggregate real quantities and prices stay constant. Aiyari (1994) proved that under some conditions, with no aggregate shock, a stationary competitive equilibrium exists for the economy. For clarity, I will formally define the stationary competitive equilibrium as follows.
A stationary competitive equilibrium is scalars \(\{R, w, K, L\}\) and functions \(\{v, g, \lambda\}\) such that:
- The value function, \(v\) satisfies the Household’s problem given \(R, w\) and \(l(s)\), and \(g\) is an optimal policy to the problem.
- \(\lambda\) is the stationary distribution that arises from \(g\) and \(\pi\) (from above).
- Asset market clears, i.e. \(K = \sum_{s \in S} \sum_{a \in A} \lambda(a,s) g(a,s)\).
- Labor market clears, i.e. \(L = \sum_{s \in S} \sum_{a \in A} \lambda(a,s) l(s)\).
- Firm’s optimality conditions are satisfied.
Good’s market will clear with a Walras’ Law-ish argument.
Now let’s build towards computing the model. Add some more parameters in the struct.
Base.@kwdef struct Household_GE{T1, T2, T3, R1, S}
β::R1 = 0.7
ρ::R1 = 2.0
ϕ::R1 = 0.0 # The borrowing limit is zero
P::T1 = [0.5 0.5; 0.2 0.8] # transition matrix for state s
l::T2 = [1.0, 5.0] # l(s)
a_min::R1 = ϕ
a_max::R1 = 5.0
points::S = 10_000
a_grid::T3 = range(a_min, a_max, points)
A::R1 = 1.20
α::R1 = 0.7
δ::R1 = 1.0 #full depreciation
end
u(c, m) = c ^ (1 - m.ρ) / (1 - m.ρ)
uprime(c, m) = c ^ (-m.ρ) uprime (generic function with 1 method)h_GE = Household_GE()Household_GE{Matrix{Float64}, Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Int64}(0.7, 2.0, 0.0, [0.5 0.5; 0.2 0.8], [1.0, 5.0], 0.0, 5.0, 10000, 0.0:0.0005000500050005:5.0, 1.2, 0.7, 1.0)It is useful to note that at stationary equilibrium, \(L\) is a constant. This is because if you solve for the invariant distribution of \(P\) (the transition probability \(\pi\) ), \(L\) can be written as \(\sum_{s \in S} \pi(s) l(s)\). Here we will numerically compute \(L\).
function compute_L(P, l)
# approximate the invariant distribution of labor endowment (inefficient)
temp = P^100
L = 0
for i in eachindex(l)
L = L + l[i] * temp[1,i]
end
return L
end
compute_L(h_GE.P, h_GE.l)3.857142857142877First let’s compute the invariant distribution given some \(R\). For simplicity let \(A\) be 1 (not the correct way if you want to calibrate the model). Since we want to include firm’s optimality conditions, modify the function for computing \(v\) to include firm’s optimality condition. Now we have \(K\) as input and then solve for \(R\) and \(w\) endogenously.
function solve_household_GE(h_GE; K = 0.5, v0 = zeros(length(h_GE.a_grid), length(h_GE.l)), tol = 1e-6)
# unpacking
(; a_grid, l, β, P, A, α, δ) = h_GE
L = compute_L(P, l)
R = A * α * K^(α - 1) * L^(1 - α) + 1 - δ# Get R from the firm's first order condition
w = A * (1 - α) * K^α * L^(- α) # Get w from the firm's first order condition
#Given K
v1 = similar(v0)
pol = similar(v0, Int)
βv = similar(v1)
n = length(a_grid)
iter = 0
while true
distance = zero(eltype(v0))
iter += 1
# Precompute βE[V(a', s')] for better performance
for s in eachindex(l)
for i in eachindex(a_grid)
accum = zero(eltype(v0))
for sprime in eachindex(l)
accum += β * P[s, sprime] * v0[i, sprime]
end
βv[i, s] = accum
end
end
for s in eachindex(l)
pol_i = 1
for i in eachindex(a_grid)
just_started = true
vmax = zero(eltype(v0))
for j in pol_i:n #take use of strictly increasing policy function
c = R * a_grid[i] + w * l[s] - a_grid[j]
if c > 0.0
v_tmp = u(c, h_GE) + βv[j, s]
if just_started
vmax = v_tmp
pol_i = j
just_started = false
elseif v_tmp > vmax
vmax = v_tmp
pol_i = j
else
break
end
end
end
v1[i, s] = vmax
pol[i, s] = pol_i
dis = abs(vmax - v0[i, s])
if dis > distance
distance = dis
end
end
end
if distance < tol
break
else
v0 .= v1
end
end
return (v = v1, pol = pol, h = h_GE, R = R, w = w)
end
solve_household_GE (generic function with 1 method)sol_2 = solve_household_GE(h_GE, K = 0.75)(v = [-17.561774071593458 -7.641618063043207; -17.509638148903143 -7.637543522427641; … ; -1.5204371786915847 -1.4164179272473385; -1.5203069683874395 -1.4163050174215155], pol = [1 324; 1 325; … ; 9364 9940; 9365 9941], h = Household_GE{Matrix{Float64}, Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Int64}(0.7, 2.0, 0.0, [0.5 0.5; 0.2 0.8], [1.0, 5.0], 0.0, 5.0, 10000, 0.0:0.0005000500050005:5.0, 1.2, 0.7, 1.0), R = 1.3729054349841805, w = 0.11440878624868113)function compute_stationary_distribution(sol; tol = 1e-8, pdf_0 = fill(1.0 / prod(size(sol.v)), size(sol.v)) )
(; pol, h) = sol
pdf_1 = similar(pdf_0)
while true
fill!(pdf_1, zero(eltype(pdf_0)))
for i in eachindex(h.a_grid)
for s in eachindex(h.l)
for sprime in eachindex(h.l)
pdf_1[pol[i, s], sprime] += h.P[s, sprime] * pdf_0[i, s]
end
end
end
distance = zero(eltype(pdf_0))
for (a, b) in zip(pdf_0, pdf_1)
distance = max(abs(a - b), distance)
end
(distance < tol) && break
pdf_0 .= pdf_1
end
return pdf_1
end compute_stationary_distribution (generic function with 1 method)pdfA = compute_stationary_distribution(sol_2)
plot(sol_2.h.a_grid, sum(pdfA, dims = 2)[:, 1], label = "pdf")The PDF is not super smooth due to computational errors. You can try adding more grids (doesn’t work very well), but the better method is to use continuous optimization. Keep in mind the numerical errors from here. It will come back later.
Of course you can also plot the CDF.
plot(sol_2.h.a_grid, cumsum(sum(pdfA, dims = 2)[:, 1]), label = "cdf")We can also compute the aggregate asset supply under the invariant distribution, which is simply:
# function for computing A under the invariant distribution
function compute_A(sol, pdf)
a = similar(sol.v)
a = sol.h.a_grid[sol.pol]
return a⋅pdf
endcompute_A (generic function with 1 method)@show compute_A(sol_2, pdfA)compute_A(sol_2, pdfA) = 1.0478295961721261.0478295961721263.3 The capital demand and supply curve
Since the aggregate labor supply \(L\) is a constant as we have argued, the only equilibrium object that has to been determined endogenously is the market clearing capital (or asset). Now, the equilibrium interest rate \(R = AαK^{α - 1} L^{1 - α} + 1 - δ\) is a one-to-one function of \(K\). The inverse function \(K(R)\) gives the capital demand curve, which is continuous and downward sloping. From the household problem, given each R, we can endogenously solve out the asset supply curve with the invariant distribution of household, which will be continuous and upward sloping. I denote the supply as \(A(R)\). The unique intersection of \(K(R)\) and \(A(R)\) gives out the equilibrium capital and the corresponding interest rate.
The following function solves out the aggregate asset supply and capital demand over a grid of possible equilibrium interest rate \(R\).
function K_demand_supply(h_GE; R_grids = range(0.95, 1.0/h_GE.β, 100), K_demand = similar(R_grids), A_supply = similar(R_grids), w_grids = similar(R_grids))
(; a_grid, l, β, P, A, α, δ) = h_GE
L = compute_L(P, l)
# get the demand function
K_demand = @. ((R_grids - 1 + δ)/(A * α * L^(1 - α)))^(1/(α - 1))
# equilibrium wages
w_grids = @. A^(1/(1-α)) * (1 - α) * α^(α/(1 - α)) / (R_grids - 1 + δ)^(α/(1-α))
# get the supply function (use the original VFI function)
for i in eachindex(R_grids)
sol = solve_household(h_GE, R = R_grids[i], w = w_grids[i])
pdf = compute_stationary_distribution(sol)
A_supply[i] = compute_A(sol, pdf)
end
return R_grids, K_demand, A_supply
end
R_grids, K_demand, A_supply = K_demand_supply(h_GE)(0.95:0.004834054834054834:1.4285714285714286, [2.5592811230762695, 2.516345927319649, 2.474342418197531, 2.4332458368040855, 2.3930322016718937, 2.353678280645201, 2.3151615639038146, 2.277460238085242, 2.2405531614551712, 2.2044198400790145 … 0.7282763869637299, 0.719867466600011, 0.7115843300746376, 0.7034246696389707, 0.695386227650787, 0.6874667953188908, 0.6796642114833356, 0.6719763614301456, 0.6644011757394387, 0.6569366291659031], [0.3935551678730847, 0.39275776490277464, 0.39199031525591366, 0.3911147352032359, 0.39089303063609615, 0.39115584409938453, 0.39160472443612426, 0.39218564845568704, 0.3928700139055018, 0.39364469678411473 … 1.2179868137367416, 1.3090272943476338, 1.4205468127468952, 1.5608866788486613, 1.7443162211259264, 1.9941783002778637, 2.3511808837319235, 2.847716210221898, 3.433878827587586, 3.958195215474183])function K_plot(K_demand, A_supply, h_GE; R_grids = range(0.95, 1.0/h_GE.β, 100))
plot(R_grids, K_demand, label = "K(R)")
plot!(R_grids, A_supply, label = "A(R)")
end
K_plot(K_demand, A_supply, h_GE)We can look at when \(A(R)\) first exceeds \(K(R)\) to have a guess on the equilibrium capital and interest rate.
function find_first(K_demand, A_supply)
idx = 0
for i in eachindex(K_demand)
if A_supply[i] > K_demand[i]
idx = i
break
end
end
return idx
end
i = find_first(K_demand, A_supply)
K1 = K_demand[i]
K2 = K_demand[i - 1]
R1 = R_grids[i - 1]
R2 = R_grids[i]
print("The equilibrium capital is between $K1 and $K2. The equilibrium interest rate is between $R1 and $R2.")The equilibrium capital is between 0.8003704925496153 and 0.8100242546324644. The equilibrium interest rate is between 1.3415584415584416 and 1.3463924963924965.Challenge for solving the equilibrium:
3.4 Solving the stationary competitive equilibrium
Now we have everything ready to compute the stationary competitive equilibrium.
The algorithm goes as follows:
- Start from a given \(K_0\)
- For each \(K_t\) and obtain \(R\) and \(w\)
- Solve the household problem given \(R\) and \(w\)
- Compute invariant distribution \(\lambda(a, s)\)
- Get \(A \equiv \sum_{s \in S} \sum_{a \in A} \lambda(a,s) g(a,s)\).
- Guess \(K_{t+1} = ϵ A + (1 - ϵ) K_t\). Go back to the second step and iterate until convergence.
compute_A(sol_2, pdfA)1.047829596172126function compute_stationary_equilibrium(h_GE; K0 = 0.70, tol = 1e-5, iterate = 0, ϵ = 0.3)
sol = solve_household_GE(h_GE, K = K0)
pdf = compute_stationary_distribution(sol)
A = compute_A(sol, pdf)
K = K0
K1 = ϵ * A + (1 - ϵ) * K
while true
distance = zero(eltype(K))
sol = solve_household_GE(h_GE, K = K)
pdf = compute_stationary_distribution(sol)
A = compute_A(sol, pdf)
K1 = ϵ * A + (1 - ϵ) * K # for speed
distance = abs(K1 - K)
(distance < tol || iterate == 100) && break
iterate += 1
K = copy(K1)
end
return (K_rce = K, sol_rce = sol, iterate = iterate)
endcompute_stationary_equilibrium (generic function with 1 method)Warning: this can get slow as the parameter values are casually set and the algorithms are not optimized. If you play around with the parameters you might notice that sometimes it doesn’t converge. I will discuss this with more detail.
@time begin
sol_GE = compute_stationary_equilibrium(h_GE)
end 0.706224 seconds (24.39 k allocations: 12.320 MiB, 2.37% compilation time)(K_rce = 0.807696820287375, sol_rce = (v = [-16.717896304757282 -7.336975047166834; -16.67191258870082 -7.3334192227477315; … ; -1.606486499114749 -1.490479825567775; -1.6063520123876334 -1.4903641858039804], pol = [1 337; 1 338; … ; 9198 9812; 9199 9813], h = Household_GE{Matrix{Float64}, Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Int64}(0.7, 2.0, 0.0, [0.5 0.5; 0.2 0.8], [1.0, 5.0], 0.0, 5.0, 10000, 0.0:0.0005000500050005:5.0, 1.2, 0.7, 1.0), R = 1.342717011889535, w = 0.12050091789432643), iterate = 7)# equilibrium interest rate
@show sol_GE.sol_rce.Rsol_GE.sol_rce.R = 1.3427170118895351.342717011889535# equilibrium wage
@show sol_GE.sol_rce.wsol_GE.sol_rce.w = 0.120500917894326430.12050091789432643function do_pol_plot_GE(sol_1)
(; pol, h)= sol_1
p1 = Plots.plot(xlabel = "Current saving, a", ylabel = "Saving next period, a' ")
for s in eachindex(h.l)
temp = h.l[s]
yvals = [[h.a_grid[pol[i, s]] for i in eachindex(h.a_grid)]]
plot!(p1, sol_1.h.a_grid, yvals, label = "l = $temp")
end
return p1
end
do_pol_plot_GE(sol_GE.sol_rce)3.5 What properties are important for solving the stationary equilibrium?
Originally the value that labor endowment process takes are 2 and 4, with same transition probability (now it’s 1 and 5). We can see that the algorithm for solving the stationary equilibrium doesn’t converge.
h_GE1 = Household_GE(l = [2.0, 4.0])Household_GE{Matrix{Float64}, Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Int64}(0.7, 2.0, 0.0, [0.5 0.5; 0.2 0.8], [2.0, 4.0], 0.0, 5.0, 10000, 0.0:0.0005000500050005:5.0, 1.2, 0.7, 1.0)Let’s look at the capital demand and supply curve. Do you notice anything different?
R_grids1, K_demand1, A_supply1 = K_demand_supply(h_GE1)
K_plot(K_demand1, A_supply1, h_GE1)The main difference is steeper asset supply curve, especially close to where the equilibrium asset and interest rate lies. We can check if we can solve for the stationary equilibrium.
sol_GE1 = compute_stationary_equilibrium(h_GE1)(K_rce = 0.689482173249647, sol_rce = (v = [-10.938373253202755 -8.118259358389057; -10.926025926801701 -8.114221970562886; … ; -1.5501109947084142 -1.4943803850038575; -1.5499807113138844 -1.4942593212902437], pol = [1 118; 1 119; … ; 9491 9788; 9492 9789], h = Household_GE{Matrix{Float64}, Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Int64}(0.7, 2.0, 0.0, [0.5 0.5; 0.2 0.8], [2.0, 4.0], 0.0, 5.0, 10000, 0.0:0.0005000500050005:5.0, 1.2, 0.7, 1.0), R = 1.3591154017758917, w = 0.11713573011418807), iterate = 100)@show sol_GE1.iteratesol_GE1.iterate = 100100Recall that 100 is the upper bound we set for the total iterations. You can change it to higher value but in short you don’t get convergence. It is also often the case with other parameter choices as long as the asset supply curve is steep around where equilibrium lies around, i.e. the elasticity for interest rate is high.
Why this is troublesome? Now because the aggregate asset supply rises sharply where \(R\) is close to the equilibrium, to find an equilibrium we need good approximations of the policy function around where it matters. However recall that we used even grids for the assets, which causes inaccuracy of the policy function as we simply don’t have enough points where it truly matters. Again this is just due to approximation errors. A naive way to resolve it is to put more grids around the equilibrium point. A better alternative is to simply obtain better precisions with continuous optimization method (you are asked to perform this in the assignment 4). (For enthusiasts, another possible solution is endogenous grid methods).
For our purpose, we can just make the asset supply less elastic. That’s why I make the possible values of labor endowments wider. The reason this works goes back to the very nature of the “incompleteness” of the model. Household can’t insure against all future risks, so they have to save more than in a complete market due to strict concavity of utility. This is called precautionary saving motivation. When the labor endowments get more fluctuant, either through wider values or smaller autocorrelation, the household saves more due to this motive. Fixing other parameters, households’ saving behavior is more driven by the uncertainty of the income shock.
In general, we should try to obtain better precisions through continuous optimizations (good exercise) and other computational techniques. An important lesson here is that precision does matter and you should always keep that in mind.
For clarity, the key here lies in the comparison between the incentives to save due to precautionary saving and higher interest rate. Just for intuition, we enhanced the first channel and made the asset supply less elastic of the second.
A final comment is that don’t be afraid when your code doesn’t work. Of course mostly it is some random mistakes and as long as you correct them, it is perfectly fine. More importantly, there are cases where something that doesn’t work actually teaches you a lesson, and often a most valuable one. Digging into some seemingly weird situations can provide valuable insights, as is the case here.
4. Transition dynamics for incomplete market (extended reading)
Economists are often interested in the short-run response of agents given an unexpected one-time shock (often referred to as the MIT shock). Now we have the tools to compute the stationary competitive equilibrium. However as you can see, even under stationarity, this is not very simple.
Assuming that households don’t have perfect foresight over prices, this implies that household needs to know everyone else’s state today (knowing the distribution is sufficient) to estimate the prices next period. In dynamic programming term, you need the whole distribution as an additional state variable. Note that the distribution is of very high dimensionality (here \(A \times s\)), taking conditional expectation over this gigantic object is merely impossible.
Krussel Smith (1998) is a very important contribution to the computation technique of the model. They are able to show that first order approximations over prices provide very accurate results because the policy function is close to linear when \(a\) is not very low.
Some latest developments of the field is the sequence Jacobian. Auclert et al. provides a computation technique that is crazy fast and the purtabation method by Bhandari et al.
5. Conclusion
Incomplete market models are the workhorse models in modern Macroeconomic research, particularly in Macro-Labor and Macro-Finance. Knowing how to compute the basic model is of vital importance. In this lecture I started with the partial equilibrium version of the heteorgeneous household saving model (known as a Bewley/Hugget model), then I built towards the general equilibrium version of the model (known as an Aiyagari model). This lecture serves more as an introduction to allude readers to dive deeper themselves.