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:

- a scalar corresponding to the interpolated value (truncated if needed);
- an index
`np`

which is such that the interpolation is performed between indices`np`

and`np+1`

.

In [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:

- previous period savings
`agrid`

(in`params`

), - current savings policy function
`ga`

(in`solution`

), - model parameters in
`params`

.

In [2]:

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

In [3]:

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

The 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:

- either the difference between the policy function
`ga`

and its update is below the threshold`tol`

; - or the number of iterations is above the number of repetitions
`maxiter`

. The output message is different in both cases.

In [4]:

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

- the sparse matrix is of size
`na*ny`

$\times$`na*ny`

; - a matrix element corresponds to the probability to switch from a pair of (savings, productivity level) to another pair of (savings, productivity level).

In [5]:

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

.

In [6]:

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

relies 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.

In [7]:

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

- the stationary distribution sums to 1;
- the capital is correctly computed;
- the aggregate savings are consistent with saving policy function;
- the aggregate savings are consistent with the asset grid;
- the consumption is properly computed.

In [8]:

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

- the Gini coefficient
`Gini`

, - the debt-to-GDP
`B/Y`

, - the public spending-to-GDP
`G/Y`

, - total consumption-to-GDP
`C/Y`

, - investment-to-GDP
`I/Y`

, - lump-sum transfer-to-GDP
`Transfers/Y`

, - total labor supply
`L`

, - Marginal propensity to consume
`MPC`

, - the proportion of constrained agents in the population
`Share of constrained agents`

.

In [9]:

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