This post will discuss various approaches to designing solvers for intertemporal optimization problems (in economics and hopefully other domains).
Consider a fairly generic intertemporal optimization problem with 1-state variable and 1-control variable.
To make things precise, this post will study the deterministic, infinite horizon T=\infty, Neoclassical growth model with log utility and full depreciation \delta =1 (so we can compare the numeric solution to the closed form solution).
To use the methods in CLMM (or any perturbation/projection method) the user can solve for the relevant equations (Euler Equations, Envelope conditions etc):
Note: there are several ways to do this
- a user can solve for these equations symbolically by hand
- a user can solve for these equations symbolically with a CAS (like
Symbolics.jl)
The user only needs to supply the objective u(s,c), transition T(s,c), & constraint g(s,c)\leq 0 - a user can solve some (or all) of these equations numerically
Next consider 3 different Fixed Point methods (inspired by CLMM):
Below is generic (1-state/1-control) Julia code for all 3 methods.
(NOT written for performance, only to discuss DSLs in economics)
using Interpolations, Plots, LinearAlgebra, SpeedMapping
##########################################################
# 1: choose approximation method. approx(gr, val0)
##########################################################
@inline approx(gr, val0) = LinearInterpolation(gr, val0, extrapolation_bc = Line())
d(x,y) = norm(x - y, Inf)
##########################################################
# 2: Define MDP. γ/S/C(s)/sp=T(s,c)/R(s,c)
##########################################################
struct MDP γ; S; C; T; R end;
RHS(P,V,s,c) = P.R(s,c) + P.γ*V(P.T(s,c))
##########################################################
# 3: solvers
##########################################################
# VFI Vanilla iteration
@inline function MDP_VFI(Vin, P::MDP; tol=1e-10, itmax::Int64=100)
POL = similar(Vin); #POLp = similar(Vin);
V = copy(Vin); Vp = similar(V);
rhs = []
#rhs = zeros(length(P1.A(P1.S[1])), length(P.S));
dist=1.0; iter = 0;
@inbounds while(dist > tol && iter<itmax)
println(iter)
v̂ = approx(P.S, V)
#if mod(iter,1) == 1
#= BEGIN VFI =#
rhs = [[RHS(P, v̂, s, c) for c ∈ P.C(s)] for s ∈ P.S]
rhs = hcat(rhs...)
@views Vp .= maximum(rhs, dims=1) |> vec
#= END VFI =#
dist = d(V, Vp) #
iter += 1
V .= Vp
#plot!(Vp) |> display
end
pol_ix = argmax(rhs, dims=1) |> vec
POL = [P.C.(P.S[j])[getindex.(pol_ix, 1)[j]] for j in eachindex(P.S)]
return V, POL
end
# VFI SpeedMapping => returns V(s)
function map_VFI!(Vout, Vin, P)
v̂ = approx(P.S, Vin)
rhs = [[RHS(P, v̂, s, c) for c ∈ P.C(s)] for s ∈ P.S]
rhs = hcat(rhs...)
Vout .= maximum(rhs, dims=1) |> vec
end
# EC DVF SpeedMapping => returns Vs(s)
function map_ECDVF!(Vsout, Vsin, P)
v̂s = approx(P.S, Vsin)
#Vout .= [β*df(s)*v̂s(T(s,df(s)/v̂s(s))) for s in S]
#Vout .= [β*df(s)*v̂s(T(s,FOC(s,v̂s))) for s in S]
Vsout .= [EC(s, v̂s) for s in P.S]
end
# EE Iteration SpeedMapping => returns POL(s)
function map_EE!(Πout, Πin, P)
Πh = approx((P.S), Πin)
sp = P.T.(P.S, Πin)
cp = Πh.(sp)
Πout .= EE1.(P.S, sp, cp) # Π Improvement/Update
#plot!(P.S, Πout) |> display
end
##########################################################
# Example 1: NGM
##########################################################
α=.33; δ=1.0; z=1.0; β=.80;
u(s,c) = c > 1e-10 ? log(c) : -1e10
TC(t,s,c;β=β) = (β^t)*(1.0/c)*s # Terminal Condition
f(s; δ=δ,z=z,α=α) = (1.0 -δ)*s + z*(s^α)
df(s;δ=δ,z=z,α=α) = (1.0 -δ) + α*z*(s^(α-1))
# sp = f(s) - c #sp>=0 => c <=f(s)
c_max(s) = f(s) - eps(); # => sp>0
c_min(s) = 0.0 + eps(); # => c>0
#
FOC(s,Vs) = clamp((df(s)/Vs(s)), c_min(s), c_max(s)) # => c(s,Vs)
EC(s, Vs;β=β) = β*df(s)*Vs(T(s,FOC(s,Vs)))
#
dVEC(s,c) = df(s)/c(s) # takes policy c(s), returns Vs(s)
#
c1_EE1(s1,s1p,c1p; β=β) = (c1p) * ((β*df(s1p))^(-1))
c1_min(s1) = 0.01 # c1 >= 0.0
c1_max(s1) = f(s1) # c1 <= f(s1) => s1p>=0
EE1(s1,s1p,c1p) = clamp(c1_EE1(s1,s1p,c1p), c1_min(s1), c1_max(s1)) # clamp(x, lo, hi)
#
s_ss = (((1/β) -1 +δ)/(α*z))^(1/(α-1))
c_ss = z*(s_ss)^(α) - δ*s_ss
v_ss = (1/(1-β)) * u(s_ss, c_ss)
vs_ss = (1/c_ss)*df(s_ss)
#
c_s_ss(s) = f(s) - s; # c(s), on sp=s Locus.
#Y_ss = z*(s_ss)^(α) #K_ss = Y_ss +(1-δ)*K_ss - C_ss
Ω1 = α/(1-α*β)
ω = (β*Ω1)^(-1)
Ω0 = (1/(1-β))*( log(ω) -(1+β*Ω1)*log(1+ω) + (1+β*Ω1)*log(z))
c_cf(s) = (1-α*β)*f(s) # (1-α*β)*z*(s^α)
v_cf(s) = Ω0 + Ω1*log(s)
vs_cf(s) = Ω1/s
# MDP
γ = β
T(s,c) = f(s) - c #(1.0 -δ)*s + z*(s^α) - c
R(s,c) = u(s,c)
##########################################
" MDP VFI "
##########################################
#S = [(0.5*s_ss : .0005 : 1.5*s_ss);]
#S = [range(0.5*s_ss, 1.5*s_ss, length= 3);]
S = [range(0.5*s_ss, 1.5*s_ss, length= 7);]
C(s;c_ss=c_ss) = [range(0.25*c_ss, f(s), length= 5*length(S));]
P1 = MDP(γ, S, C, T, R)
V0 = [(1/(1-β)) * u(s, c_s_ss(s)) for s in P1.S]
@time V_VFI, Pol_VFI = MDP_VFI(V0, P1; tol=1e-18, itmax=1_000)
#
V̂_VFI = approx(P1.S, V_VFI)
S = [range(0.5*s_ss, 1.5*s_ss, length= 200);]
C(s;c_ss=c_ss) = [range(0.25*c_ss, f(s), length= 5*length(S));]
P1 = MDP(γ, S, C, T, R)
V0 = [V̂_VFI(s) for s in P1.S]
@time V_VFI, Pol_VFI = MDP_VFI(V0, P1; tol=1e-18, itmax=1_000)
##########################################
" MDP VFI SpeedMapping "
##########################################
S = [range(0.5*s_ss, 1.5*s_ss, length= 7);]
C(s;c_ss=c_ss) = [range(0.25*c_ss, f(s), length= 5*length(S));]
P1 = MDP(γ, S, C, T, R)
V0 = [(1/(1-β)) * u(s, c_s_ss(s)) for s in P1.S]
@time V_SM = speedmapping(V0; (m!)=(Vout, Vin) -> map_VFI!(Vout, Vin, P1), tol=1e-18, Lp=Inf).minimizer
#
V̂_SM = approx(P1.S, V_SM)
S = [range(0.5*s_ss, 1.5*s_ss, length= 200);]
C(s;c_ss=c_ss) = [range(0.25*c_ss, f(s), length= 5*length(S));]
P1 = MDP(γ, S, C, T, R)
V0 = [V̂_SM(s) for s in P1.S]
@time V_SM = speedmapping(V0; (m!)=(Vout, Vin) -> map_VFI!(Vout, Vin, P1), tol=1e-18, Lp=Inf).minimizer
##########################################
" EC DVF SpeedMapping "
##########################################
S = [range(0.5*s_ss, 1.5*s_ss, length= 7);]
P1 = MDP(γ, S, C, T, R)
DV0 = [df(s) * 1.0/(c_s_ss(s)) for s in P1.S]
@time DV_ECDVF = speedmapping(DV0; (m!)=(Vout, Vin) -> map_ECDVF!(Vout, Vin, P1), tol=1e-12, Lp=Inf).minimizer
#
DV̂_ECDVF = approx(P1.S, DV_ECDVF)
S = [range(0.5*s_ss, 1.5*s_ss, length= 200);]
P1 = MDP(γ, S, C, T, R)
DV0 = [DV̂_ECDVF(s) for s in P1.S]
@time DV_ECDVF = speedmapping(DV0; (m!)=(Vout, Vin) -> map_ECDVF!(Vout, Vin, P1), tol=1e-12, Lp=Inf).minimizer
##########################################
" EE SpeedMapping "
##########################################
S = [range(0.5*s_ss, 1.5*s_ss, length= 7);]
P1 = MDP(γ, S, C, T, R)
Π0 = [c_s_ss(s) for s in P1.S]
@time Pol_EE = speedmapping(Π0; (m!)=(Πout, Πin) -> map_EE!(Πout, Πin, P1), tol=1e-12, Lp=Inf).minimizer
#
P̂ol_EE = approx(P1.S, Pol_EE)
S = [range(0.5*s_ss, 1.5*s_ss, length= 200);]
P1 = MDP(γ, S, C, T, R)
Π0 = [P̂ol_EE(s) for s in P1.S]
@time Pol_EE = speedmapping(Π0; (m!)=(Πout, Πin) -> map_EE!(Πout, Πin, P1), tol=1e-12, Lp=Inf).minimizer
###########################################
norm(V_VFI - v_cf.(S), Inf)
norm(V_SM - v_cf.(S), Inf)
norm(Pol_VFI - c_cf.(S), Inf) # inaccurate bc method minimized norm(V-Vp)
norm(DV_ECDVF - vs_cf.(S), Inf)
norm(Pol_EE - c_cf.(S), Inf)
ds = (maximum(S) - minimum(S))/(length(S)-1)
# V_VFI/dv̂_VFI(s)/POL_VFI(s)
dV_VFI = (V_VFI[2:end] - V_VFI[1:(end-1)])/ds
dv̂_VFI = approx(S[2:end], dV_VFI)
POL_VFI(s) = FOC(s,dv̂_VFI) #df(s)/v̂s_SM(s)
# V_SM/dv̂_SM(s)/POL_SM(s)
dV_SM = (V_SM[2:end] - V_SM[1:(end-1)])/ds
dv̂_SM = approx(S[2:end], dV_SM)
POL_SM(s) = FOC(s,dv̂_SM) #df(s)/v̂s_SM(s)
# V_ECDVF/dV̂_ECDVF/POL_ECDVF
dV̂_ECDVF = approx(S, DV_ECDVF)
POL_ECDVF(s) = FOC(s, dV̂_ECDVF)
function V_ECDVF(s;β=β)
T_sim = 1_000;
c_path = zeros(T_sim)
u_path = zeros(T_sim)
s_path = zeros(T_sim+1)
s_path[1] = s;
for tt =1:(T_sim)
c_path[tt] = POL_ECDVF(s_path[tt])
u_path[tt] = β^(tt-1) * u(s_path[tt], c_path[tt])
s_path[tt+1] = T(s_path[tt], c_path[tt])
end
return sum(u_path)
end
# V_EE(s)/dV_EE(s)/P̂OL_EE(s)
P̂OL_EE = approx(S, Pol_EE)
dV_EE(s) = dVEC(s, P̂OL_EE)
function V_EE(s;β=β)
T_sim = 1_000;
c_path = zeros(T_sim)
u_path = zeros(T_sim)
s_path = zeros(T_sim+1)
s_path[1] = s;
for tt =1:(T_sim)
c_path[tt] = P̂OL_EE(s_path[tt])
u_path[tt] = β^(tt-1) * u(s_path[tt], c_path[tt])
s_path[tt+1] = T(s_path[tt], c_path[tt])
end
return sum(u_path)
end
#
plot(legend=:topleft)
plot!(S, V_VFI, lab="V: VFI")
plot!(S, V_SM, lab="V: VFI SM")
plot!(S, V_ECDVF, lab="V: ECDVF")
plot!(S, V_EE, lab="V: EE")
plot!(S, v_cf, lab="V: closed form")
plot!(S, x -> (1/(1-β)) * u(x, c_s_ss(x)), lab="V: sp=s Locus, T(s, c(s))=s")
plot!([s_ss], seriestype = :vline, lab="", color="grey")
plot!([v_ss], seriestype = :hline, lab="", color="grey")
#
plot(legend=:topleft)
plot!(S, POL_VFI, lab="Policy: VFI")
plot!(S, POL_SM, lab="Policy: VFI SM")
plot!(S, POL_ECDVF, lab="Policy: ECDVF") # accurate
plot!(S, x -> P̂OL_EE(x), lab="Policy: EE") # accurate
plot!(S, c_cf, lab="Policy: closed form")
plot!(S, c_s_ss, lab="Policy: sp=s Locus, T(s, c(s))=s")
plot!([s_ss], seriestype = :vline, lab="", color="grey")
plot!([c_ss], seriestype = :hline, lab="", color="grey")
#
plot(legend=:topright)
plot!(S, x -> dv̂_VFI(x), lab="Vs: VFI")
plot!(S, x -> dv̂_SM(x), lab="Vs: VFI SM")
plot!(S, x -> dV̂_ECDVF(x), lab="Vs: ECDVF")
plot!(S, x -> dV_EE(x), lab="Vs: EE")
plot!(S, vs_cf, lab="Vs: closed form")
plot!(S, x -> df(x) * 1.0/(c_s_ss(x)), lab="Vs: sp=s Locus, T(s, c(s))=s")
plot!([s_ss], seriestype = :vline, lab="", color="grey")
plot!([vs_ss], seriestype = :hline, lab="", color="grey")
# Simulate
POL_sim(x) = P̂OL_EE(x) #approx(P1.S, Pol_VFI)
#POL_sim(x) = c_s_ss(x)
VAL_sim = approx(P1.S, V_VFI)
dV_sim(x) = dV̂_ECDVF(x)
T_sim =50;
s_sim, v_sim, dv_sim = [zeros(T_sim) for i in 1:3]
c_sim, TC_sim = [zeros(T_sim-1) for i in 1:2]
#
s0 = .75 * s_ss;
s0 = 1.25 * s_ss;
#
s_sim[1] = s0
for tt = 1:(T_sim-1)
c_sim[tt] = POL_sim(s_sim[tt])
TC_sim[tt] = TC(tt,s_sim[tt], c_sim[tt])
s_sim[tt+1] = T(s_sim[tt], c_sim[tt])
v_sim[tt] = VAL_sim(s_sim[tt])
dv_sim[tt] = dV_sim(s_sim[tt])
end
#
plot()
plot!(c_sim, lab="c_t")
plot!(s_sim, lab="s_t")
#plot!(v_sim[1:(end-1)], lab="V_t")
plot!(TC_sim, lab="β^t *λ_t *s_t ")
plot!([s_ss c_ss 0.0], seriestype = :hline, lab="", color="grey", l=:dash)
#
plot()
plot!(v_sim[1:(end-1)], lab="V_t")
plot!([v_ss], seriestype = :hline, lab="", color="grey", l=:dash)
#
plot()
plot!(dv_sim[1:(end-1)], lab="dV_t")
plot!([vs_ss], seriestype = :hline, lab="", color="grey", l=:dash)
#
# s[t] -> s_ss State variables -> ss
# a[t] -> a_ss Choice variables -> ss
# V[t] -> v_ss Value -> ss
# dV[t] -> vs_ss dValue -> ss
# TC[t] -> 0 Terminal cond -> 0
Compare value functions: note V(K_{SS})=V_{SS}=\frac{u(C_{SS})}{1-\beta}
Compare policy functions: note C(K_{SS})=C_{SS} (VFI is less accurate than EC & EE)
Compare marginal value functions: note V'(K_{SS})=V_{SS}=\frac{u'(C_{SS})}{\beta}
TS-simulation s_{0} < s_{ss}: the choice/state variables go to SS & terminal condition goes to 0 (rarely verified)
TS-simulation s_{0} > s_{ss}: the choice/state variables go to SS & terminal condition goes to 0 (rarely verified)
A few opinions about DSL design.
All solvers work with any approx(gr, val0) function chosen by the user.
This can be linear interpolation or any method from any package that conforms w/ the solver.
It helps to keep the code as generic as possible.
For example, to compute the distance between iterations I use d(x,y) = norm(x - y, Inf), however the user should be able to use any desired distance function.