Hi there,

I am solving a dynamic programming problem (first try in Julia). Previously I solved the same problem in Python. I was expecting Julia solving it much faster than Python, and I don’t get that much improvement. I would like to know if this is normal or otherwise, what I am doing wrong. Below you have the function `bellman`

which I am interested to speed up (It is partially taken from QuantEcon lectures notes).

```
mutable struct AgentStruct{T <: Float64}
r::T
R::T
y::T
β::T
b::T
b_lim::T
tau::T
theta::T
pi::T
landa::T
X::Matrix{T}
z_vals::Vector{T}
a_vals::Array{Float64,1}
u::Function
end
function ConsumerProblem(;r=0.00,
y=1.0,
β=0.995,
X=[0.5 0.5; 0.0319 0.9681],
z_vals=[0.0, 1.0],
b=0.0,
b_lim=0.0,
tau=0.02,
theta=0.15,
pi=0.2,
landa=0.0,
grid_max=8,
grid_size=100)
R = 1.0 + r
a_vals = collect(range(-b_lim, stop=grid_max, length=grid_size))
u(c::Float64, h::Float64)::Float64 = log(c) + log(1.0-h)
return AgentStruct(r, R, y, β, b, b_lim, tau, theta, pi, landa, X, z_vals, a_vals, u)
end
function bellman(cons::AgentStruct,
V::Array{Float64, 2};
ret_policy::Bool=false)
"""
The approximate Bellman operator, which computes and returns the
updated value function TV (or the V-greedy policy c if
return_policy is True).
Inputs
----------
cons : AgentStruct
An instance of ConsumerProblem that stores primitives
V : Matrix(Float64)
A Matrix with dimension
return_policy : bool, optional(default=False)
Indicates whether to return the greed policy given V or the
updated value function TV. Default is TV.
Output
-------
TV, g_a, g_d : Array(Float64)
Returns either the greed policy given V or the updated value
function TV.
"""
# Unpack parameters and grid
X, β, y, b, b_lim, tau, theta, pi, u = cons.X, cons.β, cons.y, cons.b, cons.b_lim, cons.tau, cons.theta, cons.pi, cons.u
a_vals, z_vals = cons.a_vals, cons.z_vals
# Initialize intermediate value functions and operator
TV = zeros(size(V))
V_u = ones(length(a_vals))*(-1e9)
V_e = ones(length(a_vals))*(-1e9)
V_f = ones(length(a_vals))*(-1e9)
# Initialize policy functions
g_a = Array{Int64, 2}(undef, length(a_vals), 3)
g_d = Array{Int64, 1}(undef, length(a_vals))
# Fix current a in each iteration and evaluate over the vector of possible a'
for (i_a, a) in enumerate(a_vals)
# Compute consumption for each possible employment status
c_e::Array{Float64,1} = max.(a + (1.0-tau)*(y + b) .- a_vals, 0.0)
c_u::Array{Float64,1} = max.(a + (1.0-tau)*(theta*y + b) .- a_vals, 0.0)
c_f::Array{Float64,1} = max.(a + (1.0-tau)*b .- a_vals, 0.0)
# Get maximums and maximizers for each possible employment status
V_u[i_a], g_a[i_a,1] = findmax(u.(c_u,0.0) + β*(X[1,2]*V[:,2] + X[1,1]*V[:,1]))
V_e[i_a], g_a[i_a,2] = findmax(u.(c_e,0.45) + β*(X[2,2]*V[:,2] + X[2,1]*V[:,1]))
V_f[i_a], g_a[i_a,3] = findmax(u.(c_f,0.0) + β*(X[2,2]*V[:,2] + X[2,1]*V[:,1]))
end
# Update the operator
TV[:,1] = copy(V_u)
TV[:,2] = max.(V_e, pi*V_u+ (1.0-pi)*V_f)
# Save the decision rule of accepting the job offer
g_d[:] = V_e .> pi*V_u + (1.0-pi)*V_f
if ret_policy
return TV, g_a, g_d
else
return TV
end
end
```

Testing it:

```
V_0 = ones(length(cons.a_vals), length(cons.z_vals))*(-120.0)
using BenchmarkTools
@btime V, g_a, g_d = bellman(cons, V_0, ret_policy=true)
9.679 ms (363922 allocations: 7.85 MiB)
12.801 ms (363922 allocations: 7.85 MiB)
13.033 ms (363922 allocations: 7.85 MiB)
```

In Python I am getting around 15.5ms-16.3ms

Thank you in advance,