Go back to the Main.ipynb notebook.
The first function interpLinear is a linear interpolation function, that is more specialized and hence slightly faster than the interpolation function of the LinearInterpolation package.
The function specification is interpLinear(x, xs, ys, n, lb) where:
x is a scalar where the interpolation is computed;xs and ys are sorted vectors such that (xs[i], ys[i]) for i=1,...,n is a set of n known points;n is the common length of vectors xs and ys;lb is the lower bound for the truncation of the interpolation. If the interpolation result is lower than lb, the interpolation is trauncated to lb.The function returns a pair containing:
np which is such that the interpolation is performed between indices np and np+1.function interpLinear(x::T,
xs::Vector{T},
ys::Vector{T},
n::I, lb_y::T)::Tuple{T,I} where{T<:Real,I<:Integer}
nx = searchsortedlast(xs, x)
#np is the largest index such that xs[np]≤x (and hence xs[np+1]>x). Returns 0 if x≤xs[1]. xs sorted.
#Adjust nx if x falls out of bounds of xs
if nx == 0
nx += 1
elseif (nx==n)
nx -= 1
end
# Linear interpolation in x between xs[nx] and xs[nx+1]
x_l,x_h = xs[nx],xs[nx+1]
y_l,y_h = ys[nx],ys[nx+1]
y = y_l + (y_h-y_l)/(x_h-x_l)*(x-x_l)
# returns interpolated value and corresponding index
return ((y<lb_y) ? lb_y : y), nx
end;
The second function is update_cEGM!(solution, params), where:
solution::AiyagariSolution is a mutable struct AiyagariSolution which is updated by the function;params::Params is a immutable struct Params which contains economy parameters.The function updates the consumption policy function gc (in solution) using:
agrid (in params), ga (in solution), params. function update_cEGM!(solution::AiyagariSolution,
params::Params)::Nothing
@unpack Tt,u′,na,a_min,aGrid,ny,ys = params
@unpack ga,R,w = solution
cs = similar(ga)
for iy = 1:ny
for ia = 1:na
cs[ia,iy] = find_zero(c -> c-(R*aGrid[ia] +
w*ys[iy]
-interpLinear(aGrid[ia], ga[:,iy], aGrid, na, a_min)[1] + Tt),
one(w))
# we solve for the agent's budget constraint (a fsolve is needed because
# of the labor supply Euler equation). We use a linear interpolation for the
# saving policy function -- which is somewhat 'inverted' because of the EGM
# solution.
end
end
solution.gc = cs
return nothing
end;
The third function is EulerBackEGM!(solution, params), where (as before):
solution::AiyagariSolution is a mutable struct AiyagariSolution which is updated by the function;params::Params is a immutable struct Params which contains economy parameters.The function updates the policy functions for savings ga, for consumption gc, and for labor supply gl (all of them in solution) using the Euler equations for labor and consumption.
function EulerBackEGM!(solution::AiyagariSolution,
params::Params)::Nothing
@unpack β,Tt,u′,inv_u′,na,aGrid,ny,ys,Πy = params
update_cEGM!(solution, params)
@unpack gc,R,w = solution
u′s_next = similar(gc)
u′s_next .= u′.(gc)
Eu′s = similar(gc)
Eu′s .= β*R*u′s_next*Πy'
cs = similar(gc)
as = similar(gc)
cs .= inv_u′.(Eu′s)
for ia = 1:na
for iy = 1:ny
as[ia,iy] = (aGrid[ia] + cs[ia,iy] -Tt - w*ys[iy])/R
end
end
# as is a(a')
# cs is c(a')
# ls is l(a')
# updates policy function in solution
solution.gc .= cs
solution.ga .= as
return nothing
end;
The final function (the one actually computing policy functions) is SolveEGM!(solution, params; tol::Float64=1e-6, maxiter::Int64=10000), where (as before):
solution::AiyagariSolution is a mutable struct AiyagariSolution which is updated by the function;params::Params is a immutable struct Params which contains economy parameters;tol::Float64 is a precision criterion to stop the convergence process;maxiter::Int64 is a number of maximal repetitions (in case of non-convergence of policy function). The meThe function computes the policy functions for savings ga, for consumption gc, and for labor supply gl (all of them in solution) by iterating over the function EulerBackEGM!.
The function stops when:
ga and its update is below the threshold tol;maxiter.
The output message is different in both cases.function SolveEGM!(solution::AiyagariSolution,
params::Params;
tol::T=1e-8, maxiter::I=10000)::Nothing where {T<:Real,I<:Integer}
# ierates on policy functions until convergence
as = similar(solution.ga)
as .= solution.ga
i = 0
while true
i += 1
EulerBackEGM!(solution, params) #updates policy functions once
if i%50 == 0
# The convergence is only checked by blocks of 50 iterations
test = maximum(abs.(solution.ga .- as) ./ (
abs.(solution.ga) .+ abs.(as))) #computation of the convergence criterion
# println("iteration: ", i , " ", maximum(test))
# flush(stdout) #avoids that the previous line be printed once at the end of the program
if test < tol
# convergence is reached
break
end
if i > maxiter
# maximal nb of iterations is reached (but no convergence)
println("Convergence not reached after ",i,"iterations. Criterion = ", test)
break
end
end
as .= solution.ga
end
return nothing
end;
There is one helper function transitionMat(ga,params) where:
ga::Matrix{T} is the saving policy function (a na$\times$ ny matrix defined on the product grid of assets $\times$ productivity);
params::Param is the collection of model parameters.
NB: For types, we have: T<:Real and I<:Integer.
The function returns a sparse matrix (of type SparseMatrixCSC{T,I} from the package SparseArrays) such that:
na*ny$\times$na*ny;function transitionMat(ga::Matrix{T}, params::Params{T,T1,T2,T3,T4,I})::SparseMatrixCSC{T,I} where {
T<:Real,T1<:Function,T2<:Function,
T3<:Function,T4<:Function,I<:Int64}
@unpack na,a_min,aGrid,ny,Πy = params
trans = spzeros(T, I, na*ny, na*ny)
p = zero(T)
a_val = zero(T)
ia_val = zero(I)
i_mat = zero(I)
j_mat = zero(I)
for ia = 1:na
for iy = 1:ny
a_val, ia_val = interpLinear(aGrid[ia], ga[:,iy], aGrid, na, a_min)
p = (aGrid[ia] - ga[ia_val,iy])/(
ga[ia_val+1,iy] - ga[ia_val,iy])
p = min(max(p,0.0),1.0)
j_mat = (iy-1)*na
for jy = 1:ny
i_mat = (jy-1)*na
trans[i_mat+ia_val+1,j_mat+ia] = p * Πy[iy,jy]
trans[i_mat+ia_val,j_mat+ia] = (one(T)-p) * Πy[iy,jy]
end
end
end
return trans
end;
The function actually computing the stationary distribution is stationaryDist(M; tol::Float64=1e-16, maxiter::Int64=100000), where:
M::SparseMatrixCSC{T,I} is a (sparse) transition matrix that results from function transitionMat;tol::Float64 = 1e-6 is a precision criterion to stop the convergence process;maxiter::Int64=10000 is a number of maximal repetitions (in case of non-convergence of the computation).
The function returns a vector $\Pi$ of size na*ny that verifies $\Pi M=\Pi$ and is stationary distribution -- that is known to exist. It is computed as the normalised eigenvector corresponding to the largest eigen value of the transition matrix -- which is $1$. To do so, we rely on the function powm! of the package IterativeSolvers.
function stationaryDist(M::SparseMatrixCSC{T,I}; tol::T=1e-16, maxiter::I=100000
)::Vector{T} where {T<:Real,I<:Integer}
nM = size(M)[1]
_, x = powm!(Matrix(M), ones(T, nM),
maxiter = maxiter,tol = tol)
# returns the approximate largest eigenvalue λ of M and one of its eigenvector
return x/sum(x)
end;
The function steady(params; tolEGM::T=1e-8, maxiterEGM::I=100000,tolSD::T=1e-16, maxiterSD::I=100000, tolGE::T=1e-8, maxiterGE::I=50, Rmin::T=0.995) computes the steady-state solution of the Aiyigary model, where (as before):
params::Params is a immutable struct Params which contains economy parameters;tolEGM::T=1e-8 is a precision criterion to stop the convergence process of the EGM procedure SolveEGM!;maxiterEGM::I=100000 is a number of maximal repetitions (in case of non-convergence of computations) for the EGM procedure SolveEGM!;tolSD::T=1e-16 is a precision criterion to stop the convergence of the computation of the stationary distribution stationaryDist;maxiterSD::I=100000 is a number of maximal repetitions (in case of non-convergence of computations) for the computation of the stationary distribution stationaryDist;tolGE::T=1e-8 is a precision criterion to find the general equilibrium; the criterion concerns the supply-demand gap for capital.maxiterGE::I=50 is a number of maximal repetitions for finding the GE.Rmin::T=-1.0 is an initial guess for the lower bound on interest rate. If $R_{\min} - 1 \le 0$ (as it is the case for the default value), we set Rmin = 1.The function returns the steady-state allocation under the form of a mutable struct of type AiyagariSolution{T,I}.
The function steadyrelies on previous functions, in particular SolveEGM! for computing steady-state policy functions and stationaryDist to compute the stationnary distribution.
It uses a dichotomy algorithm to solve for the GE.
function steady(params::Params;
tolEGM::T=1e-10, maxiterEGM::I=100000,
tolSD::T=1e-16, maxiterSD::I=100000,
tolGE::T=1e-6, maxiterGE::I=100, Rmin::T=-1.0, Rmax::T=-1.0)::AiyagariSolution{T,I} where {T<:Real,I<:Integer}
@unpack β,α,δ,Tt,u′,na,a_min,aGrid,ny,ys = params
# We initialize the solution
solution = AiyagariSolution(params)
# We set the initial values for the dichotomy
if Rmax > 1.0/β || Rmax < zero(Rmax)
R1 = 1.0/β
else
R1 = Rmax
end
R0 = ((Rmin-1 + δ) ≤ zero(Rmin)) ? 1.0 : Rmin
Rm = 0.5*(R1+R0)
Km = ((Rm - 1 + δ)/(α*1^(1-α)))^(1 /(α-1))
# Start of the dichotomy
for kR = 1:maxiterGE
solution.R = Rm
solution.w = (1-α)*((Rm - (1-δ))/α)^(α/(α-1))
SolveEGM!(solution, params, tol=tolEGM, maxiter=maxiterEGM)
@unpack ga = solution
solution.transitMat = transitionMat(ga,params)
solution.stationaryDist = reshape(stationaryDist(solution.transitMat,tol=tolSD,maxiter=maxiterSD), na, ny)
Am = sum(repeat(aGrid,ny).*vec(solution.stationaryDist))
if abs(Km - Am) < tolGE
# The dichotomoy has converged: we fill the solution structure
println("Markets clear! in: ",kR," iterations")
as = similar(ga) #policy rules as a function of beginning of period asset
cs = similar(ga) #policy rules as a function of beginning of period asset
# we 'invert' policy functions (to obtain policy rules as a function of beginning of period asset)
for ia = 1:na, iy = 1:ny
as[ia,iy] = interpLinear(aGrid[ia], ga[:,iy], aGrid, na, a_min)[1]
cs[ia,iy] = Rm*aGrid[ia] - as[ia,iy] + solution.w*ys[iy] + Tt
end
solution.ga = as
solution.gc = cs
# We compute aggregate quantities
solution.A = sum(solution.stationaryDist.*as) #aggregate savings
solution.C = sum(solution.stationaryDist.*cs) #aggregate consumption
solution.L = 1.0
solution.K = solution.L*(((Rm-1)+δ)/α)^(one(T)/(α-one(T))) #aggregare capital
solution.Y = (solution.K)^α*(solution.L)^(1-α)
solution.B = (solution.A)-(solution.K)
solution.G = solution.Y- δ*solution.K - solution.C
# We check Euler equations by computing their residuals
solution.residEuler = u′.(cs) - β*Rm*reshape(
solution.transitMat'*reshape(u′.(cs),na*ny,1),na,ny)
# Returning the result (successful case)
return solution
break
end
# We update the interval [R0, Rm] based on the comparison
# between supply and demand for savings.
if (Am > Km)
R1 = Rm
else
R0 = Rm
end
Rm = 0.5*(R1+R0)
Km = ((Rm-1 + δ)/(α*1^(1-α)))^(1 /(α-1))
println("ir: ",fmt.(100*([R0,R1,Rm].-1.0),6),", supply: ",fmt(Am),
", demand-supply gap: ",Km-Am)
flush(stdout)
end
println("Markets did not clear")
# Returning the result (failure case)
return solution
end;
The function check_solution(solution::AiyagariSolution, params::Params) checks some consistency requirements for the Aiyagari solution solution corresponding parameters Params.
It returns a Boolean. If true, all requirements are fullfilled. If false, it returns an error message indicating the specific error.
We check the following elements:
function check_solution(solution::AiyagariSolution, params::Params)::Bool
@unpack β,α,δ,Tt,u′,na,a_min,aGrid,ny,ys = params
@unpack ga,gc,R,w,A,K,C,L,stationaryDist = solution
if !(sum(stationaryDist) ≈ 1.0)
println("error in stationary distribution: ",
sum(stationaryDist))
return false
end
if !(L*((R - (1-δ))/α)^(1/(α-1)) ≈ K)
println("error in aggregate savings",L*((1/β - (1-δ))/α)^(1/(α-1)) - K )
return false
end
if !(A ≈ sum(stationaryDist.*ga))
println("error in aggregate savings:", A,sum(stationaryDist.*ga))
return false
end
if abs(A-sum(stationaryDist.*repeat(params.aGrid,1,ny))) > na*ny*1e-10
println("difference in savings: ",
A, sum(stationaryDist.*repeat(params.aGrid,1,ny)),
A-sum(stationaryDist.*repeat(params.aGrid,1,ny)))
return false
end
C_ = w -A + Tt + R*sum(stationaryDist.*repeat(params.aGrid,1,ny))
if !(C ≈ C_)
println("error in aggregate consumption: ", C_, C, C_-C)
return false
end
return true
end;
The function describe_solution(solution::AiyagariSolution, params::Params) returns a dictionary that characterizes several elements of the solution:
Gini,B/Y,G/Y,C/Y,I/Y,Transfers/Y,L,MPC,Share of constrained agents.function describe_solution(solution::AiyagariSolution, params::Params)
@unpack β,α,δ,Tt,u′,na,a_min,aGrid,ny,ys = params
@unpack ga,gc,R,w,A,C,K,L,G,Y,B,stationaryDist,residEuler = solution
#computing MPC
diff_gc = gc[2:end,:] - gc[1:end-1,:]
diff_ga = aGrid[2:end,:] - aGrid[1:end-1,:]
mpc = sum(diff_gc.*stationaryDist[1:end-1,:]./diff_ga)
p1 = plot(repeat(aGrid,1,ny),residEuler,legend=:none) # by bins
p2 = plot(repeat(aGrid,1,ny),ga .- repeat(aGrid,1,ny),legend=:none)
p3 = plot(repeat(aGrid,1,ny),gc,legend=:none)
layout = @layout [a ; b c]
p = plot(p1, p2, p3, layout=layout,
title=["Resid Euler" "Labor supply" "Net saving" "Consumption Rule"])
return Dict("Gini" => Gini(ga, stationaryDist),
"B/Y" => B/(4*Y),
"G/Y" => G/Y,
"C/Y" => C/Y,
"I/Y" => δ*K/Y,
"Transfers/Y" => Tt/Y,
"L" => L,
"MPC" => mpc,
"Share of constrained agents" => sum(stationaryDist[1,:]),
"K/Y" => K/(4*Y)
), p
end;