Help speeding up finite horizon DP problem

#1

Crossposted in the Julia discourse, but trying here to get feedback on my technique.

I am struggling to speed up the computation of solutions for a finite horizon dynamic programming problem. I am using interpolation to avoid computing the value function outside of points in the grid. Right now it takes my computer around 7 minutes to solve this 5-period model.

Please have a look and let me know what I can improve. Any suggestions you may have will help me!

Thanks!

# VFI finite horizon
using Optim
using Distributions
using Random
using Interpolations

# transformation functions
# open interval (a, b)
function logit(x; a = 0, b = 1)
    (b - a) * (exp(x)/(1 + exp(x))) + a
end

function utility(c::R, L::R) where R <: Real
    if c <= 0 || L <= 0
        return -Inf
    else
        return log(c) + log(1 - L)
    end
end

const T = 5 # terminal period
const β = 0.95 # discount factor
const r = 0.05 # interest rate
const dist = LogNormal(0, 1) # distribution of ξ
const n = 1 # number of points of support of ξ

Random.seed!(1234)
const w = Vector{Float64}(undef, T) # exogenous wages

w .= (900 .+ 20.0 .* (1:T) .- 0.5 .* (1:T).^2)

const grid_A = -1000:10:8_000

const ξ = quantile.(dist, 0:1/n:prevfloat(1.0))

const Vdict = Array{Float64}(undef, (n, length(grid_A), T))
const Cdict = Array{Float64}(undef, (n, length(grid_A), T))
const Ldict = Array{Float64}(undef, (n, length(grid_A), T))
const A1dict = Array{Float64}(undef, (n, length(grid_A), T))
const convdict = Array{Bool}(undef, (n, length(grid_A), T))

# solve last period
for s in 1:length(grid_A)
    for i in 1:n
        opt = optimize( x -> -utility(w[T]*x + grid_A[s]*(1+r), x), 0.0, 1.0 )
        Ldict[i, s, T] = Optim.minimizer(opt)
        Cdict[i, s, T] = (w[T] + ξ[i])*Ldict[i, s, T] + grid_A[s]*(1+r)
        A1dict[i, s, T] = 0.0
        Vdict[i, s, T] = -Optim.minimum(opt)
    end
end

interp_func(t) = LinearInterpolation(grid_A, sum(Vdict[i, :, t] for i in 1:n) ./ n, extrapolation_bc = Line())

# solve rest of periods
@time for t in T-1:-1:1
    for s in 1:length(grid_A)
        for i in 1:n
            initial_x = [A1dict[i, s, t+1], 0.0]
            # x[1] is assets to carry forward, x[2] is labor supply
            opt = optimize( x -> -(utility(logit(x[2])*(w[t] + ξ[i])+ grid_A[s]*(1+r) - x[1], logit(x[2])) +
                        β*interp_func(t+1)(x[1]) ),
                  initial_x, method = NewtonTrustRegion(), autodiff = :forward, iterations = 100_000)
            Vdict[i, s, t] = -Optim.minimum(opt)
            xstar = Optim.minimizer(opt)
            Ldict[i, s, t] = logit(xstar[2])
            A1dict[i, s, t] = xstar[1]
            Cdict[i, s, t] = Ldict[i, s, t] * (w[t] + ξ[i]) + grid_A[s] * (1+r) - A1dict[i, s, t]
            convdict[i, s, t] = Optim.converged(opt)
        end
    end
    println("period ", t, " finished")
end
#2

HI @amrods, I have to keep my answer short but the first thing I suggest you do is break your code up into functions. In particular, those loops where most of the execution time is occurring should be inside functions. The reason is that Julia does type inference at the function boundaries, and this is easier for the compiler when the functions are relatively small. They also need to be type stable, so one function can call another with prior knowledge of the type that will be returned.

Please read this lecture: https://lectures.quantecon.org/jl/need_for_speed.html

1 Like
#3

Thanks @john.stachurski, that was also suggested in the Julia discourse forum. I revised the program, see below.
However, I posted here as well to gather feedback about any “meta” tricks that economists employ to solve these models. For example, one can use linear algebra to solve a particular class of DP problems, which is something that programmers (like in the Julia forum) would overlook.

using Optim
using Interpolations

struct DP{Ty <: Real}
    u::Function
    T::Int
    β::Ty
    r::Ty
    n::Int
    w::Vector{Ty}
    grid_A::Vector{Ty}
    ξ::Vector{Ty}
end

function solvelast!(dp::DP{Ty}, Ldict::Array{Ty}, Cdict::Array{Ty}, A1dict::Array{Ty}, Vdict::Array{Ty}) where Ty <: Real
    utility = dp.u
    grid_A = dp.grid_A
    n = dp.n
    w = dp.w
    ξ = dp.ξ
    r = dp.r
    T = dp.T
    for s in 1:length(grid_A)
        for i in 1:n
            opt = optimize( x -> -utility(w[T]*x + grid_A[s]*(1+r), x), 0.0, 1.0 )
            Ldict[i, s, T] = Optim.minimizer(opt)
            Cdict[i, s, T] = (w[T] + ξ[i])*Ldict[i, s, T] + grid_A[s]*(1+r)
            A1dict[i, s, T] = 0.0
            Vdict[i, s, T] = -Optim.minimum(opt)
        end
    end
    return Ldict, Cdict, A1dict, Vdict
end

# transformation functions
# open interval (a, b)
function logit(x; a = 0, b = 1)
    (b - a) * (exp(x)/(1 + exp(x))) + a
end

# interpolation
interp_func(t) = LinearInterpolation(grid_A, sum(Vdict[i, :, t] for i in 1:n) ./ n, extrapolation_bc = Line())

function solverest!(dp::DP{Ty}, Ldict::Array{Ty}, Cdict::Array{Ty}, A1dict::Array{Ty}, Vdict::Array{Ty}, convdict::Array{Bool}; t0::Int=1) where Ty <: Real
    utility = dp.u
    grid_A = dp.grid_A
    n = dp.n
    w = dp.w
    ξ = dp.ξ
    r = dp.r
    T = dp.T

    for t in T-1:-1:t0
        for s in 1:length(grid_A)
            for i in 1:n
                initial_x = [A1dict[i, s, t+1], 0.0]
                # x[1] is assets to carry forward, x[2] is labor supply
                opt = optimize( x -> -(utility(logit(x[2])*(w[t] + ξ[i])+ grid_A[s]*(1+r) - x[1], logit(x[2])) +
                            β*interp_func(t+1)(x[1]) ),
                            initial_x, method = NewtonTrustRegion(), autodiff = :forward, iterations = 100_000)
                Vdict[i, s, t] = -Optim.minimum(opt)
                xstar = Optim.minimizer(opt)
                Ldict[i, s, t] = logit(xstar[2])
                A1dict[i, s, t] = xstar[1]
                Cdict[i, s, t] = Ldict[i, s, t] * (w[t] + ξ[i]) + grid_A[s] * (1+r) - A1dict[i, s, t]
                convdict[i, s, t] = Optim.converged(opt)
            end
        end
        println("period ", t, " finished")
    end
    return Ldict, Cdict, A1dict, Vdict, convdict
end

I then instantiate the problem and solve it:

using Distributions

function u(c::R, L::R) where R <: Real
    if c <= 0 || L <= 0
        return -Inf
    else
        return log(c) + log(1 - L)
    end
end

T = 5 # terminal period
β = 0.95 # discount factor
r = 0.05 # interest rate
dist = LogNormal(0, 1) # distribution of ξ
n = 1 # number of points of support of ξ

w = (900 .+ 20.0 .* (1:T) .- 0.5 .* (1:T).^2)

grid_A = collect(-1000:10.0:8_000)

ξ = quantile.(dist, 0:1/n:prevfloat(1.0))

Vdict = Array{Float64}(undef, (n, length(grid_A), T))
Cdict = Array{Float64}(undef, (n, length(grid_A), T))
Ldict = Array{Float64}(undef, (n, length(grid_A), T))
A1dict = Array{Float64}(undef, (n, length(grid_A), T))
convdict = Array{Bool}(undef, (n, length(grid_A), T))

dp = DP(u, T, β, r, n, w, grid_A, ξ)
solvelast!(dp, Ldict, Cdict, A1dict, Vdict)
solverest!(dp, Ldict, Cdict, A1dict, Vdict, convdict)