# My julia code is somehow much slower than the matlab code

**URL:** <https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792>\
**Category:** New to Julia\
**Tags:** question, performance, matlab\
**Created:** [December 17, 2022, 9:01pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792 "2022-12-17T21:01:15Z")\
**Posts on this page:** 1\
**Showing post:** 42

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [December 22, 2022, 6:57pm UTC](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792/42 "2022-12-22T18:57:24Z")

</div>

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](https://github.com/jesusfv/Comparison-Programming-Languages-Economics).  
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](https://github.com/vfitoolkit/VFIToolkit-matlab).m, even though VFI is one of the biggest computational challenges for economists!
3. You can see VFI & a bunch of other methods [here](https://discourse.julialang.org/t/solving-the-4-quadrants-of-dynamic-optimization-problems-in-julia-help-wanted/73285) & [here](https://discourse.julialang.org/t/non-dsls-for-solving-intertemporal-optimization-problems-not-just-for-econ/85011)

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

```julia
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:

```julia
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

```julia
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)) 

```

---

_[View the full topic](https://discourse.julialang.org/t/my-julia-code-is-somehow-much-slower-than-the-matlab-code/91792)._