My julia code is somehow much slower than the matlab code

New to Julia
, ,
42

A few comments

  1. you create nz, nk but then instead often use 5, 499
  2. this problem uses VFI to solve the MOST ROUTINE economics problem (Neoclassical growth model)
    Jesus compared MANY languages for this problem here.
    There have been a bunch of GSOC/JSOC improving performance for this exact problem in Julia.
    Julia does not have a high performance VFI solver like VFIToolkit.m, even though VFI is one of the biggest computational challenges for economists!
  3. You can see VFI & a bunch of other methods here & here

I made the language a bit more generic (at the cost of performance)

using Printf, Statistics, LinearAlgebra;
Base.show(io::IO, f::Float64) = @printf(io, "%.4f", f)

function main()
    ## (a) ##
    # First, set up the given data#
    α = 0.261; β = 0.99; δ = 0.0176;
    #
    u(c) = log(c);
    f(z,k;α=α) = z*(k^α)
    con(z,k,kp;δ=δ) = f(z,k) + k*(1 - δ) - kp # makes it slower?
    #
    zGrid = [0.9526, 0.9760, 1, 1.0246, 1.0497]
    PiZ = [0.9149 0.0823 0.0028 0      0
           0.0206 0.9163 0.0618 0.0014 0
           0.0005 0.0412 0.9167 0.0412 0.0005
           0      0.0014 0.0618 0.9163 0.0206
           0      0      0.0028 0.0823 0.9149]
    # Next, get the kStar and kL and kH #
    kStar = ( (1.0/β - 1.0 + δ) / α )^(1.0 / (α - 1.0))
    kL = kStar * 0.85; kH = 1.15 * kStar;
    println(kStar); println(kL); println(kH);
    # Setting up an 499 points of kgrid #
    kGrid = range(kL, kH, 499); println(median(kGrid))

    ## (b) ##
    # First, set up the length variables #
    nk = length(kGrid); nz = length(zGrid)

    # Make settings for the iteration #
    # First, make an array for each z_i #
    k′Grid = repeat(kGrid, 1, nz, nk)
    kGridMat = permutedims(k′Grid, (3, 2, 1))
    # 
    k′GridMat = k′Grid .* (con.(zGrid', kGridMat, k′Grid) .> 0)
    #k′GridMat = k′Grid .* (kGridMat .^ α .* zGrid' .+ kGridMat .* (1 - δ) .- k′Grid .> 0)
    #
    replace!(k′GridMat, 0 => NaN)  # k′GridMat[k′GridMat .== 0] .= NaN 
    # 
    logC = u.(con.(zGrid', kGridMat, k′GridMat)) 
    #logC = u.(kGridMat .^ α .* zGrid' .+ kGridMat .* (1 - δ) .- k′GridMat) 
    #
    replace!(logC, NaN => -Inf)    # logC[isnan.(logC)] .= -Inf

    # precision and initial distance #
    precision = 0.01; distance = 100    

    # start iteration #
    i = 1
    V = zeros(nk, nz)
    V′ = similar(V)
    TV = zeros(1, nz, nk)
    VPiZ = similar(V)
    logCC = similar(logC)
    while distance > precision
        mul!(VPiZ, V, transpose(PiZ))
        logCC .= logC .+ VPiZ .* β
        maximum!(TV, logCC)
        V′ .= transpose(reshape(TV, nz, nk))
        distance = maximum(j-k for (j,k) in zip(V′,V))
        V .= V′
        i += 1
    end
    # Add one more iteration to find Tk_index
    _, k_index = findmax(logCC, dims=1)
    mul!(VPiZ, V′, PiZ')
    logCC .= logC .+ VPiZ .* β
    _, Tk_index = findmax(logCC, dims=1)

    # Reviving k-rules from the k_index array #
    k′sol = transpose(reshape(k′Grid[Tk_index], nz, nk))
    k′sol_prev = transpose(reshape(k′Grid[k_index], nz, nk)) # matrix of previous capital to compare with solutions #
    dist_k′sol = maximum(abs(i-j) for (i,j) in zip(k′sol,k′sol_prev))

    # answers #
    println(i); println(distance); println(dist_k′sol)
end
@time main()

The world would be better off if there was a generic Julia VFI w/ n-state variables, a generic utility (return/cost) function, and a generic transition function…

Does anyone know why:

k′GridMat = k′Grid .* (kGridMat .^ α .* zGrid' .+ kGridMat .* (1 - δ) .- k′Grid .> 0)
logC = u.(kGridMat .^ α .* zGrid' .+ kGridMat .* (1 - δ) .- k′GridMat)

is faster than the more generic

u(c) = log(c);
f(z,k;α=α) = z*(k^α)
con(z,k,kp;δ=δ) = f(z,k) + k*(1 - δ) - kp 

k′GridMat = k′Grid .* (con.(zGrid', kGridMat, k′Grid) .> 0)
logC = u.(con.(zGrid', kGridMat, k′GridMat))
2 Likes