Hi guys,
I’m using Julia to compute a dynamic programming model. There is one particular part of the model that takes a long time to run (about 15 seconds, which is a lot, given that there is an iterative procedure involved). I ran the @time macro on it and the output is:
15.600687 seconds (75.16 M allocations: 31.067 GiB, 43.40% gc time)
This part of the code is as follows:
function solv_active(pa::Param,pr::Prices,q::Array{<:Real},vn_ret::Array{<:Real},vh_ret::Array{<:Real},vd_ret::Array{<:Real})
# vn_ret is the value function of retired renters in the first period
# vh_ret is the value function of retired owners in the first period
# vd_ret is the value function of retired defaulters in the first period
@unpack Jret, J, β, ϑ, H_num, na, nϵ, nb, b_grid, ϵ_chain, b_min, M_num, infm = pa
@unpack q_b, w = pr
ϵ_grid = ϵ_chain.state_values # idiosyncratic prod states
nw = Jret-1 # number of working periods
# Renters - array pre-allocation
vn_act = Array{Float64}(undef,nϵ,na,nb,nw) # value function
cnf_act = similar(vn_act) # consumption function
bnf_act = similar(vn_act) # asset function
htnf_act = similar(vn_act) # rent decision
hnf_act = similar(vn_act) # house decision
mnf_act = similar(vn_act) # mortgage decision
go_act = Array{Int}(undef,nϵ,na,nb,nw) # decision to own
# Homeowners - array pre-allocation
vh_act = Array{Float64}(undef,nϵ,na,M_num,H_num,nb,nw) # value function
chf_act = similar(vh_act) # consumption function
bhf_act = similar(vh_act) # asset function
hthf_act = similar(vh_act) # rent decision
hhf_act = similar(vh_act) # house decision
mhf_act = similar(vh_act) # mortgage decision
πf_act = similar(vh_act) # mortgage payment function
gh_act = Array{Int}(undef,nϵ,na,M_num,H_num,nb,nw) # decision
# Defaulters - array pre-allocation
vd_act = Array{Float64}(undef,nϵ,na,M_num,nb,nw) # value function
cdf_act = similar(vd_act) # consumption function
bdf_act = similar(vd_act) # asset function
htdf_act = similar(vd_act) # rent decision
πdf_act = similar(vd_act) # mortgage balance
# last year of active life
Threads.@threads for ib = 1:nb
# renter
for ia = 1:na, iϵ = 1:nϵ
bn,vn,cn,htn,hn,mn,gn = solve_act_last_rent(pa,pr,ib,iϵ,ia,q[iϵ,ia,:,:,:,Jret-1],vn_ret[iϵ,ia,:],vh_ret[iϵ,ia,:,:,:])
vn_act[iϵ,ia,ib,end] = vn # value function
cnf_act[iϵ,ia,ib,end] = cn # consumption function
bnf_act[iϵ,ia,ib,end] = bn # bond function
htnf_act[iϵ,ia,ib,end] = htn # decision to rent
hnf_act[iϵ,ia,ib,end] = hn # own house
mnf_act[iϵ,ia,ib,end] = mn # mortage decision
go_act[iϵ,ia,ib,end] = gn # decision to own
end
end
Threads.@threads for ib = 1:nb
# homeowner
for ih = 1:H_num, im = 1:M_num, ia = 1:na, iϵ = 1:nϵ
bh,vh,ch,hth,hh,mh,πh,gh = solve_act_last_home(pa,pr,ib,iϵ,ih,im,ia,q[iϵ,ia,:,:,:,Jret-1],vn_ret[iϵ,ia,:],vh_ret[iϵ,ia,:,:,:],vd_ret[iϵ,ia,:,:])
vh_act[iϵ,ia,im,ih,ib,end] = vh # value function
chf_act[iϵ,ia,im,ih,ib,end] = ch # consumption function
bhf_act[iϵ,ia,im,ih,ib,end] = bh # bond function
hthf_act[iϵ,ia,im,ih,ib,end] = hth # decision to rent
hhf_act[iϵ,ia,im,ih,ib,end] = hh # decision to own home
mhf_act[iϵ,ia,im,ih,ib,end] = mh # mortgage decision
πf_act[iϵ,ia,im,ih,ib,end] = πh # mortgage payment function
gh_act[iϵ,ia,im,ih,ib,end] = gh # decision to own, sell or rent
end
end
Threads.@threads for ib = 1:nb
# defaulter
for im = 1:M_num, ia = 1:na, iϵ = 1:nϵ
bd,vd,cd,htd,πd = solve_act_last_def(pa,pr,ib,iϵ,im,ia,vn_ret[iϵ,ia,:],vd_ret[iϵ,ia,:,:])
vd_act[iϵ,ia,im,ib,end] = vd # value function
cdf_act[iϵ,ia,im,ib,end] = cd # consumption function
bdf_act[iϵ,ia,im,ib,end] = bd # bond function
htdf_act[iϵ,ia,im,ib,end] = htd # rent function
πdf_act[iϵ,ia,im,ib,end] = πd # payment function
end
end
# Remaining years of active life
for ij = 1:nw-1
Threads.@threads for ib = 1:nb
# Renter
for ia = 1:na, iϵ = 1:nϵ
bn,vn,cn,htn,hn,mn,gn = solve_act_rent(pa,pr,ib,iϵ,ia,nw-ij,q[iϵ,ia,:,:,:,Jret-1-ij],vn_act[:,ia,:,end-ij+1],vh_act[:,ia,:,:,:,end-ij+1])
vn_act[iϵ,ia,ib,end-ij] = vn # value function
cnf_act[iϵ,ia,ib,end-ij] = cn # consumption function
bnf_act[iϵ,ia,ib,end-ij] = bn # bond function
htnf_act[iϵ,ia,ib,end-ij] = htn # decision to rent
hnf_act[iϵ,ia,ib,end-ij] = hn # decision to own house
mnf_act[iϵ,ia,ib,end-ij] = mn # mortgage decision
go_act[iϵ,ia,ib,end-ij] = gn # decision to own
end
end
Threads.@threads for ib = 1:nb
# Homeowner
for ih = 1:H_num, im = 1:M_num, ia = 1:na, iϵ = 1:nϵ
bh,vh,ch,hth,hh,mh,πh,gh = solve_act_home(pa,pr,ib,iϵ,ih,im,ia,nw-ij,q[iϵ,ia,:,:,:,Jret-1-ij],vn_act[:,ia,:,end-ij+1],vh_act[:,ia,:,:,:,end-ij+1],vd_act[:,ia,:,:,end-ij+1])
vh_act[iϵ,ia,im,ih,ib,end-ij] = vh # value function
chf_act[iϵ,ia,im,ih,ib,end-ij] = ch # consumption function
bhf_act[iϵ,ia,im,ih,ib,end-ij] = bh # bond function
hthf_act[iϵ,ia,im,ih,ib,end-ij] = hth # decision to rent
hhf_act[iϵ,ia,im,ih,ib,end-ij] = hh # decision to own house
mhf_act[iϵ,ia,im,ih,ib,end-ij] = mh # mortgage decision
πf_act[iϵ,ia,im,ih,ib,end-ij] = πh # mortgage payment function
gh_act[iϵ,ia,im,ih,ib,end-ij] = gh # decision to keep owning, selling and renting or buying another house
end
end
Threads.@threads for ib = 1:nb
# Defaulter
for im = 1:M_num, ia = 1:na, iϵ = 1:nϵ
bd,vd,cd,htd,πd = solve_act_def(pa,pr,ib,iϵ,im,ia,nw-ij,vn_act[:,ia,:,end-ij+1],vd_act[:,ia,:,:,end-ij+1])
vd_act[iϵ,ia,im,ib,end-ij] = vd # value function
cdf_act[iϵ,ia,im,ib,end-ij] = cd # consumption function
bdf_act[iϵ,ia,im,ib,end-ij] = bd # bond function
htdf_act[iϵ,ia,im,ib,end-ij] = htd # rent function
πdf_act[iϵ,ia,im,ib,end-ij] = πd # payment function
end
end
end
return vn_act, cnf_act, bnf_act, htnf_act, hnf_act, mnf_act, go_act, # renter functions
vh_act, chf_act, bhf_act, hthf_act, hhf_act, mhf_act, πf_act, gh_act, # homeowner
vd_act, cdf_act, bdf_act, htdf_act, πdf_act # defaulter
end
So what I’m doing here is: for each point in the state space and for each agent type, I solve a dynamic programming problem. The solution, the maximizer and other info are stored in a matrix. This function is type stable (I ran @code_warntype on it) and all sub-functions are also type stable.
From the @time output, the biggest issue seems to be memory allocations and garbage collection, although I’m not entirely sure how I’m going to solve that. Perhaps there is a better way to organize this function that I’m overlooking?