# Non-DSLs for solving intertemporal optimization problems (not just for econ)

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

1. a user can solve for these equations symbolically by hand
2. 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
3. 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)), 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 = 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 = 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 = 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.

Some examples from the ecosystem:

1. LinearSolve.jl
solves equations Au=b for an unknown vector (of numbers) u when A\b isn’t good enough.
prob = LinearProblem(A, b ; u0 = zero(b))
sol = solve(prob,alg) where alg is taken from the ecosystem.
sol = solve(prob) uses the default solver
sol = solve(prob,IterativeSolversJL_GMRES())
2. NonlinearSolve.jl
solves nonlinear system f(u;p)=0 for an unknown vector (of numbers) u, where p is a vector of parameters.
f(u,p) = u .* u .- p
u0 = @SVector[1.0, 1.0]
p = 2.0
prob = NonlinearProblem{false}(f, u0, p)
sol = solve(prob, alg) where alg is taken from the ecosystem.
sol = solve(prob, NewtonRaphson(), tol = 1e-9)
(doesn’t wrap NLboxsolve.jl)
3. Optimization.jl
solves \underset{u}{\min} f(u;p) for an unknown vector (of numbers) u, where p is a vector of parameters.
f(u,p) = (p - u)^2 + p * (u - u^2)^2
u0 = zeros(2)
p = [1.0,100.0]
prob = OptimizationProblem(f,u0,p) unconstrained
prob = OptimizationProblem(f,u0,p, lb = [-1.0,-1.0], ub = [1.0,1.0])
sol = solve(prob,alg) e.g. alg=NelderMead() or alg=BFGS() etc
4. Integrals.jl
interface for numerical integration, computes \int_{u=lb}^{u=ub} f(u,p)du
f(u,p) = sum(sin.(u .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]
prob = IntegralProblem(f,lb,ub,p)
solve(prob,alg)
Can it integrate over more general sets?
5. OrdinaryDiffEq.jl
solves explicit initial value ODEs of the form
u^{(n)} = f(t, u, u^{1}, \dots, u^{n-1}; p)
u(0) = u_{0}, \dots, u^{n-1}(0)=u_{0}^{n-1}
where u(t) is an uknown vector of functions of time t and p is a vector of parameters.
I don’t know if it handles general implicit ODEs: 0 = f(t, u, u^{1}, \dots, u^{n}; p)
function f!(du,u,p,t)
x,y,z = u
σ,ρ,β = p
dx = σ*(y-x)
dy = x*(ρ-z) - y
dz = x*y - β*z
du .= (dx, dy, dz)
end
u0 = [1.0,0.0,0.0]
tspan = (0.0,1.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(f!,u0,tspan,p)
sol = solve(prob, alg) # over 100 algorithms for ODEs

1. StochasticDiffEq.jl
solves systems of SDEs du = f(u,p,t)dt + g(u,p,t)dW.
The solution, u(t,W_{t}), is a function of time and a Weiner process.
α=1; β=1; # parameters
u₀=1/2; # initial value
f(u,p,t) = α*u # drift
g(u,p,t) = β*u # diffusion term
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u₀,tspan) (why doesn’t the Problem depend on p?)
sol = solve(prob,alg,dt=0.01)
It also handles systems of RODEs du = f(u,p,t,W) where W is a Weiner process.
It’s not clear to me from docs what types of SDEs it can solve. Is it only diffusions and RODEs?
Also see jump-dffusions
2. DelayDiffEq.jl
solves \frac{du}{dt} = f(t, u(t), \left\{ u(\tau): \tau \leq t \right\}; p) DEs that depend on history.
Currently delays are expressed through a history function h(p,t).
f(du,u,h,p,t) define the DDE
p0 = 0.2; q0 = 0.3; v0 = 1; d0 = 5
p1 = 0.2; q1 = 0.3; v1 = 1; d1 = 1
d2 = 1; beta0 = 1; beta1 = 1
p = (p0,q0,v0,d0,p1,q1,v1,d1,d2,beta0,beta1,tau)
tspan = (0.0,10.0)
u0 = [1.0,1.0,1.0]
h(p, t) = ones(3)
tau = 1; lags = [tau];
prob = DDEProblem(f,u0,h,tspan,p; constant_lags=lags)
sol = solve(prob,alg)
3. DAEs
Solves (ordinary) initial value (explicit?) DAEs u^{(n)} = f(t, u, u^{1}, \dots, u^{n-1}; p) where some equations in the function f are algebraic (do not depend on derivatives).
f(out,du,u,p,t)
u₀ = [1.0, 0, 0]
du₀ = [-0.04, 0.04, 0.0]
tspan = (0.0,100000.0)
differential_vars = [true,true,false]
prob = DAEProblem(f,du₀,u₀,tspan,differential_vars=differential_vars)
sol = solve(prob,IDA())
4. BoundaryValueDiffEq.jl
Solves (non-algebraic) boundary value ODEs
\frac{du}{dt} = f(t, u(t); p)
f!(du,u,p,t)
tspan = (0.0,1.0)
bc1!(residual, u, p, t) boundary function
u_guess = [-1.6, -1.7] initial guess (constant currently, not a function)
prob = BVProblem(f!, bc1!, u_guess, tspan) (no parameters?)
sol = solve(prob, alg, dt=0.05)

The Problem types I’ve seen:
LinearProblem(A, b ; u0 = zero(b))
NonlinearProblem{false}(f, u0, p)
OptimizationProblem(f,u0,p, lb = lb, ub = ub)
IntegralProblem(f,lb,ub,p)
ODEProblem(f!,u0,tspan,p)
SDEProblem(f,g,u₀,tspan)
DDEProblem(f,u0,h,tspan,p; constant_lags=lags)
DAEProblem(f,du₀,u₀,tspan,differential_vars=differential_vars)
BVProblem(f!, bc1!, u_guess, tspan)
RODEProblem(f,u0,tspan)
DiscreteProblem(sir_model, u₀, tspan, p)
JumpProblem(sir_model, prob, Direct())

(currently) Non-SciMl solvers:

1. MatrixEquations.jl
Solves f(u;p)=0 where u is a matrix of numbers and p are parameters (usually matrices).
Currently Lyapunov, Sylvester and Riccati matrix equations…
arec(A, B, R, Q, S; as = false, rtol::Real = nϵ, orth = false)
2. JuMP.jl
Uses Model() and optimize!() instead of SciML’s Problem() and solve()
JuMP Model() are defined w/ macros unlike SciML afaik.
I’m not sure the JuMP API allows the user to pass parameters as easily as SciML, this is important for comparative statics & estimation…
model = Model(Ipopt.Optimizer)
@variable(model, x)
@variable(model, y >= 0)
@variable(model, z >= 0)
@objective(model, Max, x)
@constraint(model, x + y + z == 1)
@constraint(model, x * x + y * y - z * z <= 0)
@constraint(model, x * x - y * z <= 0)
optimize!(model)
1 Like

No, that solver is not in the defaults. That’s a bad choice.

You’d have to manually parameterize them.

Look at the DAE pages.

The default parameters are just no parameters.

In what sense?

SDDEProblem and special forms like DynamicalODEProblem are missing. JumpProcesses.jl as well.

1 Like