Go back to the Main.ipynb notebook.
The function convertBasisE10(vE, E, N) converts a number expressed in base-E (and represented as the vector vE) into its decimal representation p (as an integer).
Convention: vE(1) represents the idiosyncratic state in the present period and vE(N) represents the past idosyncratic state N-1 periods
function convertBasisE10(vE::Vector{I}, E::I)::I where {I<:Integer}
# converts the a vector vE into an integer p
@assert (minimum(vE)≥1)&&(maximum(vE)≤E)
N = length(vE)
p = one(I) #otherwise p lieas between 0 and (E+1)^N - 1
for k ∈ eachindex(vE)
p += (vE[k]-1)*E^(N-k)
end
return p
end;
The function convertBasis10E(p, E, N) converts a number expressed in base-10 (and represented as the integer p) into its base-E representation vE (as a vector of length N). This is the inverse of convertBasisE10.
Convention: vE(1) represents the idiosyncratic state in the present period and vE(N) represents the past idosyncratic state N-1 periods
function convertBasis10E(p::I, E::I, N::I)::Vector{I} where{I<:Integer}
vE = zeros(I, N)
@assert (p≥1)&&(p≤E^N)
ptemp = p-1;
for k ∈ eachindex(vE)
ptemp, r = divrem(ptemp,E)
vE[N-k+1] = 1+r
end
return vE
end;
The function historySizes(N::I, Πy::Matrix{T};
maxiter::I=1000000, tol::T=1e-16) computes the size of (positive-size) truncated histories. It takes as input:
N: the truncation length;Πy: the transition matrix of the productivity process;maxiter and tol: controls for the convergence of the function powm! to compute the eigenvector of Πy.The function returns a t-uple with three elements:
ind_h corresponding to positive-size histories;S_h corresponding to the positive sizes of histories; note that the pair (ind_h,S_h) allows one to reconstruct the sparse vector of history sizes). y0h[ind_h] of positive-size histories.function historySizes(N::I, Πy::Matrix{T};
maxiter::I=1000000, tol::T=1e-16) where{I<:Integer, T<:Real}
ny = size(Πy,1) # nb of producitivty states
_,Sy = powm!(Πy', ones(typeof(Πy[1]), ny),maxiter = maxiter,tol = tol)
Sy /= sum(Sy) # distribution over producitivity levels
Ntot = ny^N # total nb of truncated histories
Sh = zeros(typeof(Πy[1]), Ntot)
# initialization of size of truncated histories
y0h = zeros(I, Ntot)
# initialization of current productivity indices of history
# We loop over all truncated histories
for h ∈ eachindex(Sh)
vE = convertBasis10E(h, ny, N)
# vector of producitvity levels for
# truncated history h
Sh[h] = Sy[vE[N]]
# distribution according to the
# terminal productity level
y0h[h] = vE[1]
# current prod. level of h
for j = N-1:-1:1
Sh[h] *= Πy[vE[j+1],vE[j]]
# We move backwards and compute the probability to move
# from productivity vE[j] to productivity vE[j+1]
end
end
ind_h, S_h = findnz(sparsevec(Sh))
return ind_h, S_h, y0h[ind_h]
end;
The function historyDist(ny::I, N::I, stationaryDist::Matrix{T},
transitMat::SparseMatrixCSC{T,I}) computes the distribution of (positive-size) truncated histories. It takes as input:
ny: the nb of productivity states;N: the truncation length;stationaryDist: the stationary distribution over truncated histories; transitMat: the transition matrix over the product grid asset $\times$ productivity. The function returns a matrix $((na \times ny), ny^N)$ containing the distribution of truncated histories over the product grid asset $\times$ productivity. More precisely, if statDist_h = historyDist(ny, N, stationaryDist, transitMat), then statDist_h[:,h] is the distribution of truncated history h over the product grid asset $\times$ productivity.
function historyDist(ny::I, N::I, stationaryDist::Matrix{T},
transitMat::SparseMatrixCSC{T,I})::Matrix{T} where{I<:Integer, T<:Real}
na = div(size(transitMat,1),ny) # size of asset grid
Ntot = ny^N # nb of truncated histories
# initialization of the stat distribution: matrix of size (na*ny, Ntot)
statDist_h = zeros(T, length(stationaryDist), Ntot)
temp_v = zeros(T, na)
# Loop over all truncated histories
for h ∈ 1:Ntot
vE = convertBasis10E(h, ny, N)
shift0 = (vE[N]-1)*na # shift for the index of vE[N] in the product grid
statDist_h[1+shift0:na+shift0,h] = stationaryDist[:,vE[N]]
for j=N-1:-1:1
shift1 = (vE[j+1]-1)*na
shift0 = (vE[j]-1)*na
temp_v .= statDist_h[1+shift1:na+shift1,h]
statDist_h[:,h] .= zero(T)
statDist_h[1+shift0:na+shift0,h] .= (
transitMat[1+shift0:na+shift0,1+shift1:na+shift1] * temp_v)
end
end
return statDist_h
end;
The function historyTrans(N::I,ind_h::Vector{I},Πy::Matrix{T}) computes the transition matrix for truncated histories. It takes as input:
N: the truncation length;ind_h: the index of non-zero size truncated histories;Πy: the productivity transition matrix.The function returns a matrix ((na $\times$ ny), ny^N) containing the distribution of truncated histories over the product grid asset $\times$ productivity. More precisely, if statDist_h = historyDist(ny, N, stationaryDist, transitMat), then statDist_h[:,h] is the distribution of truncated history h over the product grid asset $\times$ productivity.
function historyTrans(N::I,ind_h::Vector{I},Πy::Matrix{T}) where{I<:Integer,T<:Real}
ny = length(Πy[:,1])
Ntot = length(ind_h)
Πh = spzeros(T,Ntot,Ntot)
for (i,h) ∈ enumerate(ind_h)
vE_h = convertBasis10E(h, ny, N)
for (j,ht) ∈ enumerate(ind_h)
vE_ht = convertBasis10E(ht, ny, N)
if all(vE_ht[2:end] .== vE_h[1:end-1])
Πh[i,j] = Πy[vE_h[1],vE_ht[1]]
end
end
end
return Πh
end;
The function credit_constrained_h(shareCC::T, S_h::Vector{T}, x_h::Vector{T};method::String="first larger") identifies credit-constraint histories. It takes as input:
shareCC: the target share of credit-constrained agents;S_h: the distribution of truncated histories;x_h: a vector of the same size as S_h;method: a string caharacterizing the selction mechanism.The function returns a pair:
i_cc indicating the number of credit-constrained truncated histories;The credit-constrained histories are computed such that based the first i_cc truncated histories (ordered along x_h) have a size similar to the target shareCC. The meaning of "similar" is specified by the method: either the first larger or the closest.
function credit_constrained_h(shareCC::T, S_h::Vector{T}, x_h::Vector{T};method::String="first larger") where{T<:Real}
ind_x_indiv = sortperm(x_h,rev = true)
S_h_sorted = S_h[ind_x_indiv]
if method=="closest"
i_cc = argmin(abs.(cumsum(S_h_sorted) .- shareCC)) #closest distance
else #(method=="first larger")
i_cc = findlast(x->x<0.0, cumsum(S_h_sorted) .- shareCC)#at least as large
end
return i_cc, ind_x_indiv[1:i_cc]
end;
The function Truncate(N::Integer, # length of the truncation
solution::AiyagariSolution,
params::Params;maxiter=1000000,tol=1e-16) computes the the truncated model. It takes as input:
N: the truncation length;solution: the solution of the Aiyagary model (see SolveAiyagari for the resolution and Parameters for the details of the structures);params: the economy parameters.It returns a Truncationstructure containing:
N: the truncation length;Ntot: the nb of truncated histories;ind_h: the indices of non-zero histories;allocation_trunc: the truncated allocation;ξs: the residual heterogeneity parameters $\xi$s.function Truncate(N::Integer, # length of the truncation
solution::AiyagariSolution,
params::Params;maxiter=1000000,tol=1e-16)
@unpack β,α,δ,γ,u,u′,u′′,na,a_min,aGrid,ny,ys,Πy = params
@unpack ga,gc,R,w, A,K,L,transitMat,stationaryDist,residEuler = solution
# Identifying positive-size truncated histories
(ind_h, S_h, ind_y0_h) = historySizes(N, Πy)
y0_h = [ys[i] for i in ind_y0_h]
T = typeof(Πy[1,1])
# Distribution of truncated histories
statDist_h = historyDist(ny, N, stationaryDist, transitMat)[:,ind_h]
Ntot = length(ind_h)#nb of history with positive size
Π_h = historyTrans(N,ind_h,Πy)
#Compute truncated allocations
c_h = sum(statDist_h .* repeat(gc[:],1,Ntot), dims=1)[:]./S_h
a_beg_h = sum(statDist_h .* repeat(aGrid[:],ny,Ntot), dims=1)[:]./S_h
a_end_h = sum(statDist_h .* repeat(ga[:],1,Ntot), dims=1)[:]./S_h
u_h = sum(statDist_h .* repeat(u.(gc[:]),1,Ntot), dims=1)[:]./S_h
u′_h = sum(statDist_h .* repeat(u′.(gc[:]),1,Ntot), dims=1)[:]./S_h
u′′_h = sum(statDist_h .* repeat(u′′.(gc[:]),1,Ntot), dims=1)[:]./S_h
resid_E_h = sum(statDist_h .* repeat(residEuler[:],1,Ntot), dims=1)[:]./S_h
# Indentify credit constrained histories
share_cc = sum(stationaryDist[1,:]) #Share of credit constrained agents
nb_cc_h, ind_cc_h = credit_constrained_h(share_cc, S_h, resid_E_h;method="first larger")
# The truncated allocation structure
allocation_trunc = Allocation_trunc(
S_h=S_h,
Π_h=Π_h,
y0_h=y0_h,
a_beg_h=a_beg_h,
a_end_h=a_end_h,
c_h=c_h,
u_h=u_h,
u′_h=u′_h,
u′′_h=u′′_h,
resid_E_h=resid_E_h,
nb_cc_h=nb_cc_h,
ind_cc_h=ind_cc_h)
# The ξs allocation structure
ξu0 = u_h./u.(c_h)
ξu1 = u′_h./u′.(c_h)
ξu2 = u′′_h./u′′.(c_h)
ξs = ξs_struct(
ξu0=ξu0,
ξu1=ξu1,
ξu2=ξu2)
# The Lagrange mutiplier is constructed in the function
# LagrangeMult defined below
lagrange_mult,θ = LagrangeMult(Ntot,
allocation_trunc,solution,
params,fit_theta=true)
# The Truncation structure
return Truncation(
N = N,
Ntot = Ntot,
ind_h = ind_h,
allocation_trunc = allocation_trunc,
ξs = ξs,
lagrange_mult = lagrange_mult,
θ = θ)
end;
The function LagrangeMult(trunc::Truncation,
solution::AiyagariSolution,
params::Params; fit_theta::Bool=true) computes the Lagrange mutipliers associated to the Ramsey program in the truncated economy. These Lagrange multipliers are needed for the model simulation in the presence of aggregate shocks. The function takes as input:
trunc: the Truncation outcome (as a Truncation structure);solution: solution of the Aiyagari model;params: the economy parameters;fit_theta: a Boolean wich is true when the value of the θ parameter must be set so as to verify the planner's FOC on lump-sum transfer. If no, the default θ value is taken.The function check_truncated_model(trunc::Truncation,
solution::AiyagariSolution, tol::T=1e-6) checks the consistency of the truncectiion result. It takes as input:
trunc: the Truncation outcome (as a Truncation structure);solution: solution of the Aiyagari model;tol: the tolerance for difference in numerical checks.The function returns a Boolean and checks three elements:
An error message is displayed when an error occurs.
function LagrangeMult(Ntot::Int,allocation_trunc::Allocation_trunc,solution::AiyagariSolution,
params::Params; fit_theta::Bool=true) where{Int<:Integer}
@unpack β,α,δ,θ,Tt,u,u′,u′′,na,a_min,aGrid,ny,ys = params
@unpack R,L,K = solution
@unpack S_h,Π_h,y0_h,a_beg_h,a_end_h,c_h,u_h,u′_h,u′′_h,nb_cc_h,ind_cc_h = allocation_trunc
T = typeof(β)
Π_h_bar = spdiagm(one(T)./S_h)*Π_h'*spdiagm(S_h)
Id = sparse(I,Ntot,Ntot)
Pc = sparse(ind_cc_h,ind_cc_h,ones(T,nb_cc_h),Ntot,Ntot)
P = Id - Pc
FKK = α*(α-1)*K^(α-2)*L^(1-α)
FLK = α*(1-α)*K^(α-1)*L^(-α)
Γ = β*sparse(R*Π_h +
repeat((S_h .* (FKK*a_beg_h + FLK*y0_h))',Ntot,1))
M = Pc + P*sparse(
(β*FKK*repeat((S_h.*u′_h)',Ntot,1)*Π_h_bar) +
(Id-Γ)*diagm(u′′_h)*(Id-R*Π_h_bar))
λ = sparse(M\(P*(Id-Γ)*u′_h))
ψh= u′_h - spdiagm(u′′_h)*(Id-R*Π_h_bar)*λ
if fit_theta
F_θ(x) = sum(S_h.*ψh) .- x.*(-Tt)^(x-1)
θ_f = fzero(F_θ, zero(T), one(T))
else
θ_f = θ
end
λt = Π_h_bar*λ
Φ = (β*FKK*repeat((S_h.*u′_h)',Ntot,1)*Π_h_bar)*λ
return Lagrange_mult(Vector(λ),Vector(λt),Vector(ψh),Vector(Φ)),θ_f
end;
function check_truncated_model(trunc::Truncation,
solution::AiyagariSolution, params::Params,
tol::T=1e-6)::Bool where{T<:Real}
@unpack α,β,Tt = params
@unpack A,R,K,L = solution
@unpack Ntot, allocation_trunc, lagrange_mult, θ = trunc
@unpack S_h,Π_h,y0_h,a_beg_h,a_end_h,u′_h,u′′_h,nb_cc_h,ind_cc_h = allocation_trunc
@unpack λ,λt,ψh,Φ = lagrange_mult
if !(abs.(sum(S_h.*a_beg_h) - sum(S_h.*a_end_h)) < tol)
println("error in end and beginning distributions")
@show (sum(S_h.*a_beg_h) - sum(S_h.*a_end_h))
return false
end
if !(maximum(abs.(S_h.*a_beg_h - Π_h'*(S_h.*a_end_h))) < tol)
println("error in transition matrix")
@show (maximum(abs.(S_h.*a_beg_h - Π_h*(S_h.*a_end_h))))
return false
end
if !(abs(sum(S_h.*a_beg_h)-A) < tol)
println("error in aggregate savings")
@show (sum.(S_h.*a_beg_h)-A)
return false
end
Id = sparse(I,Ntot,Ntot)
Pc = sparse(ind_cc_h,ind_cc_h,ones(T,nb_cc_h),Ntot,Ntot)
P = Id - Pc
FKK = α*(α-1)*K^(α-2)*L^(1-α)
FLK = α*(1-α)*K^(α-1)*L^(-α)
Γ = β*sparse(R*Π_h +
repeat((S_h .* (FKK*a_beg_h + FLK*y0_h))',Ntot,1))
if !(maximum(abs.(P*ψh - P*Γ*ψh - P*Φ))<tol)
println("error in ψ for unconstrained histories")
@show maximum(abs.(P*ψh-P*Γ*ψh-P*Φ))
return false
end
if !(maximum(Pc*λ)<tol)
println("error in λ for credit-constrained histories")
@show maximum(Pc*λ)
return false
end
if !(maximum(ψh - (u′_h-diagm(u′′_h)*(λ - R*λt)))<tol)
println("error in ψ (general expression)")
@show maximum(ψh - (u′_h-diagm(u′′_h)*(λ - R*λt)))
return false
end
if !(sum(S_h.*ψh) - θ*(-Tt)^(θ-1) < tol)
println("error in θ")
@show sum(S_h.*ψh) - θ*(-Tt)^(θ-1)
return false
end
return true
end;