How can I transform this matlab code into Julia

General Usage
, ,
21 
julia> sub2ind(shape, A) = map(eachcol(A)) do c
         LinearIndices(shape)[c...]
       end
sub2ind (generic function with 1 method)

julia> sub2ind((2,2), A)
3-element Vector{Int64}:
 3
 4
 1

This assumes that each column of A stores Cartesian indices. Works in higher dimensions too

julia> B = [1 2 1; 2 2 1; 2 1 2]
3×3 Matrix{Int64}:
 1  2  1
 2  2  1
 2  1  2

julia> sub2ind((2,2,2), B)
3-element Vector{Int64}:
 7
 4
 5

You seem to be translating MATLAB code one-to-one to Julia. This might not be the best idea. Julia is not MATLAB and doing things the “Julian” way is often more efficient, elegant, generic and faster. We can help you optimize a piece of Julia code.

22

@skleinbo Thanks a lot for your kind help. Can you help me to improve my Julia codes. It takes very long time. Here is the Julia codes.

using SpecialFunctions, LinearAlgebra #I需要使用LinearAlgebra,efrc需要使用SpecialFunctions
using XLSX,MATLAB
include("tauchen.jl")
include("ind_trans.jl")

x = 1

β = 0.87
δ = 0.06
γ = 2
α = 0.35
wage = 1
za = 1
ϵ = 3
μ = ϵ / (ϵ - 1)
D = 1

phil = 0.00

nzp = 2
zp_sd_target = 0.43
piL = 0.80
zLspace = range(0.0001, 0.9999, length=1000000)
difference = piL .* (zLspace .- 1) .^ 2 + (1 - piL) .* ((piL .- piL .* zLspace) ./ (1 - piL)) .^ 2 .- zp_sd_target^2;

value, index = findmin(abs.(difference))
zp_l = zLspace[index]

zp_h = (1 - piL * zp_l) / (1 - piL)
zp_grid = [zp_l, zp_h]

ntau = 1
if ntau == 2
    tau_grid = range(-0.29, 0.29, length=ntau)
    tauprob = [0.81 0.19; 0.19 0.81]
elseif ntau == 1
    tau_grid = [0]
    tauprob = 1
end

nq = 1
if nq == 2
    q_grid = range(-0.18, 0.18, length=nq)
elseif nq == 1
    q_grid = [0]
    qprob = 1
end


nzt = 11
rho_zt = 0.59
sigma_zt = 0.13
mu_zt = -(sigma_zt^2) / (2 * (1 + rho_zt))
@show zt_grid = tauchen(nzt, mu_zt, rho_zt, sigma_zt, 3)[1]
ztprob = tauchen(nzt, mu_zt, rho_zt, sigma_zt, 3)[2]

runexp = 1
nr = 6
r_grid = [0.01, 0.02, 0.03, 0.04, 0.05, 0.10]

if runexp == 0
    rho_r = 0.50
    sigma_r = 0.0086
    mu_r = 0.03 * rho_r
    r_grid = tauchen(nr - 1, mu_r, rho_r, sigma_r, 2.014)[1]
    rprob = tauchen(nr - 1, mu_r, rho_r, sigma_r, 2.014)[2]
    r_grid = [r_grid findmax(r_grid)[1]]
    rprob = [rprob zeros(nr - 1, 1); 0.00 0.00 0.00 0.00 0.00 1.00]
end

chi0 = 0.98
chi1 = 0.047
psi = 3.2
k_l = 0.01
k_h = 6.0
nk = 120
k_grid = range(k_l, k_h, length=nk)

a_l = 0.01
a_h = 3.0
na = 120
a_grid = range(a_l, a_h, length=na)

n_choice = nk * na
n_state = nzp * nzt * nr * ntau * nq

r_grid_ind = 1:nr
zt_grid_ind = 1:nzt
zp_grid_ind = 1:nzp
tau_grid_ind = 1:ntau
q_grid_ind = 1:nq

Q = repeat(q_grid, nzp * nzt * nr * ntau, 1)
Q_ind = repeat(q_grid_ind, nzp * nzt * nr * ntau, 1)


TAU = repeat(tau_grid, nq, 1)
TAU = sort!(TAU, dims=1) #dims=1按行排序 dims = 2按列排序

TAU = repeat(TAU, nzp * nzt * nr, 1)
TAU_ind = repeat(tau_grid_ind, nq, 1)
TAU_ind = sort!(TAU_ind, dims=1)
TAU_ind = repeat(TAU_ind, nzp * nzt * nr, 1)

R = repeat(r_grid, ntau * nq, 1)
R = sort!(R, dims=2)
R = repeat(R, nzp * nzt, 1)
R_ind = repeat(r_grid_ind, ntau * nq, 1)
R_ind = sort!(R_ind, dims=1)
R_ind = repeat(R_ind, nzp * nzt, 1)

ZT = repeat(zt_grid', nr * ntau * nq, 1)
ZT = sort!(ZT, dims=1)
ZT = repeat(ZT, nzp, 1)

ZT_ind = repeat(zt_grid_ind, nr * ntau * nq, 1)
ZT_ind = sort!(ZT_ind, dims=1)
ZT_ind = repeat(ZT_ind, nzp, 1)

ZP = repeat(zp_grid, nzt * nr * ntau * nq, 1)
ZP = sort!(ZP, dims=1)
ZP_ind = repeat(zp_grid_ind, nzt * nr * ntau * nq, 1)
ZP_ind = sort(ZP_ind, dims=1);

EXOG = [ZP ZT R TAU Q]
EXOG_ind = [ZP_ind ZT_ind R_ind TAU_ind Q_ind]
EXOG = EXOG'
EXOG_ind = EXOG_ind'

A = repeat(a_grid, nk, 1)
A = sort!(A, dims=1)
K = repeat(k_grid, na, 1)
G = [A K]

prob = zeros(n_state, n_state);

for i_state = 1:n_state
    izp = EXOG_ind[1, i_state]
    izt = EXOG_ind[2, i_state]
    ir = EXOG_ind[3, i_state]
    itau = EXOG_ind[4, i_state]
    iq = EXOG_ind[5, i_state]

    for i_state_next = 1:n_state
        izpnext = EXOG_ind[1, i_state_next]
        iztnext = EXOG_ind[2, i_state_next]
        irnext = EXOG_ind[3, i_state_next]
        itaunext = EXOG_ind[4, i_state_next]
        iqnext = EXOG_ind[5, i_state_next]

        prob[i_state, i_state_next] = ztprob[izt, iztnext] * I[izp, izpnext] * rprob[ir, irnext] * tauprob[itau, itaunext] * I[iq, iqnext]
    end
end


iter_vfi = 0
iter_howard = 0
error_vfi = 10^15
error_vfi_old = error_vfi
VFIcontinue = 1
error_tolerance = 10^(-6)

if γ < 2.5 && runexp == 1
    penalty = -10^(9)
else
    penalty = -10^(12)
end

Nh = 20

indicator = ones(n_choice, n_state)
pol_ind_kp = zeros(Int, n_choice, n_state)
pol_ind_ap = zeros(Int, n_choice, n_state)
V_new = zeros(n_choice, n_state)
V_old = penalty .* ones(n_choice, n_state)
Vh_old = penalty .* ones(n_choice, n_state)
Vh_new = penalty .* ones(n_choice, n_state);

collateral = chi0 .* G[:, 1] + chi1 .* (exp.(G[:, 2]) .- 1) .- G[:, 2]

for i_state = 1:n_state
    izp = EXOG_ind[1, i_state]
    izt = EXOG_ind[2, i_state]
    ir = EXOG_ind[3, i_state]
    itau = EXOG_ind[4, i_state]
    iq = EXOG_ind[5, i_state]

    for iak_state = 1:n_choice

        ik = Int(rem(iak_state - 1, nk) + 1)
        ia = Int(floor((iak_state - 1) / nk) + 1)

        dummy = ones(na * nk, 1)
        l = phil + ((za * exp(zt_grid[izt]) * zp_grid[izp])^((ϵ - 1) / (1 + α * (ϵ - 1)))) * (Complex(k_grid[ik] + q_grid[iq])^(α * (ϵ - 1) / (1 + α * (ϵ - 1)))) * (D^(ϵ / (1 + α * (ϵ - 1)))) * (μ^(-ϵ / (1 + α * (ϵ - 1)))) * ((wage * (1 + tau_grid[itau]) / (1 - α))^(-ϵ / (1 + α * (ϵ - 1))))
        y = (za * exp(zt_grid[izt]) * zp_grid[izp]) * (Complex(k_grid[ik] + q_grid[iq])^(α)) * ((l - phil)^(1 - α))
        T = wage * tau_grid[itau] * l + wage * (1 + tau_grid[itau]) * phil - (r_grid[ir] + δ) * q_grid[iq]
        pi = D * (y^((ϵ - 1) / ϵ)) - wage * (1 + tau_grid[itau]) * l

        c = pi - (r_grid[ir] + δ) * k_grid[ik] + (1 + r_grid[ir]) * a_grid[ia] + T .- G[:, 1] - (psi .* (G[:, 2] .- k_grid[ik]) .^ 2) ./ (2 * k_grid[ik])

        CAPeffective = k_grid[ik] + q_grid[iq]

        if CAPeffective <= 0
            dummy[:] .= -1
        end

        dummy[real.(c).<=0] .= -1
        dummy[collateral.<0] .= -1

        if findmax(dummy)[1] < 0
            indicator[iak_state, i_state] = 0
        end
    end
end

ZP_mat = repeat(transpose(EXOG[1, :]), na * nk, 1)
ZT_mat = repeat(transpose(EXOG[2, :]), na * nk, 1)
R_mat = repeat(transpose(EXOG[3, :]), na * nk, 1)
Q_mat = repeat(transpose(EXOG[5, :]), na * nk, 1)

AP_mat = repeat(G[:, 1], 1, n_state)
KP_mat = repeat(G[:, 2], 1, n_state)

l_vfi = zeros(Complex{Float64}, n_choice, n_state)
y_vfi = zeros(Complex{Float64}, n_choice, n_state)
T_vfi = zeros(Complex{Float64}, n_choice, n_state)
pi_vfi = zeros(Complex{Float64}, n_choice, n_state)

location = zeros(Int, n_choice, n_state)

EXP = zeros(n_choice, n_state);


for iak_state = 1:n_choice
    ik = CartesianIndices((nk, na))[iak_state][1]
    ia = CartesianIndices((nk, na))[iak_state][2]

    l_vfi[iak_state, :] = transpose(phil .+ ((za * exp.(EXOG[2, :]) .* EXOG[1, :]) .^ ((ϵ - 1) / (1 + α * (ϵ - 1)))) .* (Complex.(k_grid[ik] .+ EXOG[5, :]) .^ (α * (ϵ - 1) / (1 + α * (ϵ - 1)))) .* (D^(ϵ / (1 + α * (ϵ - 1)))) .* (μ^(-ϵ / (1 + α * (ϵ - 1)))) .* ((wage .* (1 .+ EXOG[4, :]) ./ (1 - α)) .^ (-ϵ / (1 + α * (ϵ - 1)))))
    y_vfi[iak_state, :] = transpose((za * exp.(EXOG[2, :]) .* EXOG[1, :]) .* (Complex.(k_grid[ik] .+ EXOG[5, :]) .^ (α)) .* ((l_vfi[iak_state, :] .- phil) .^ (1 - α)))
    T_vfi[iak_state, :] = transpose(wage .* EXOG[4, :] .* l_vfi[iak_state, :] .+ wage .* (1 .+ EXOG[4, :]) .* phil .- (EXOG[3, :] .+ δ) .* EXOG[5, :])
    pi_vfi[iak_state, :] = transpose(D .* (y_vfi[iak_state, :] .^ ((ϵ - 1) / ϵ)) - wage .* (1 .+ EXOG[4, :]) .* l_vfi[iak_state, :])
end

piplusT = zeros(na, n_choice)
c = zeros(na, n_choice)
V_sub = zeros(na, n_choice)
collateral = chi0 .* AP_mat + chi1 .* (exp.(KP_mat) .- 1) - KP_mat

while VFIcontinue == 1
    error_vfi_old = error_vfi
    for i_state in 1:n_state
        prob_mat = repeat(transpose(prob[i_state, :]), n_choice, 1)
        EXP[:, i_state] = sum(prob_mat[:, :] .* V_old[:, :], dims=2)
    end
    for iak_state in 1:n_choice
        ik = Int(rem(iak_state - 1, nk) + 1)
        ia = Int(floor((iak_state - 1) / nk) + 1)

        piplusT = repeat(transpose(pi_vfi[iak_state, :]), n_choice, 1) + repeat(transpose(T_vfi[iak_state, :]), n_choice, 1)
        c = piplusT .- (R_mat .+ δ) * k_grid[ik] + (1 .+ R_mat) * a_grid[ia] - AP_mat - (psi * (KP_mat .- k_grid[ik]) .* (KP_mat .- k_grid[ik])) ./ (2 * k_grid[ik])
        CAPeffective = KP_mat + Q_mat

        V_sub = ((1 ./ abs.(c) .^ ((γ - 1) * ones(n_choice, n_state))) .- 1) ./ (1 - γ) + β * EXP

        @. V_sub[(real(c)<=0)|(collateral<0)|(CAPeffective<=0)] = penalty
        V_new[iak_state, :] = findmax(V_sub, dims=1)[1]
        location[iak_state, :] = first.(Tuple.(findmax(V_sub, dims=1)[2]))

        pol_ind_kp[iak_state, :] = first.(Tuple.(CartesianIndices((nk, na))[location[iak_state, :]]))
        pol_ind_ap[iak_state, :] = last.(Tuple.(CartesianIndices((nk, na))[location[iak_state, :]]))

    end

    ### HOWARD's Improment Algorithm
    Vh_old = V_new
    iter_howard = 0

    while iter_howard <= Nh

        for iak_state = 1:n_choice
            ik = Int(rem(iak_state - 1, nk) + 1)
            ia = Int(floor((iak_state - 1) / nk) + 1)

            c = transpose(pi_vfi[iak_state, :] - (EXOG[3, :] .+ δ) * k_grid[ik] + (1 .+ EXOG[3, :]) * a_grid[ia] + T_vfi[iak_state, :] - a_grid[pol_ind_ap[iak_state, :]] - (psi * (k_grid[pol_ind_kp[iak_state, :]] .- k_grid[ik]) .* (k_grid[pol_ind_kp[iak_state, :]] .- k_grid[ik])) ./ (2 * k_grid[ik]))
            Vh_new[iak_state, :] = ((1 ./ abs.(c) .^ ((γ - 1) * ones(1, nzt * nzp * nr * ntau * nq))) .- 1) ./ (1 - γ) + β .* sum(prob .* Vh_old[location[iak_state, :], :], dims=2)
        end


        Vh_new = Vh_new .* indicator + Vh_old .* (1 .- indicator)
        Vh_old = copy(Vh_new)
        iter_howard += 1
    end

    ### HOWARD END ##


    if Nh == -1
        error_vfi = findmax(abs.(V_new[:] - V_old[:]))[1]
        V_old = copy(V_new)
    else
        error_vfi = findmax(abs.(Vh_new[:] - V_old[:]))[1]
        V_old = copy(Vh_new)
    end

    println("--------------------------------------")

    println("error_vfi_old = $error_vfi_old")
    println("error_vfi = $error_vfi")

    iter_vfi += 1

    println("iter_vfi = $iter_vfi")

    if error_vfi < error_tolerance
        VFIcontinue = 0
    end

end


V_final = zeros(na, nk, nzp, nzt, nr, ntau, nq)
a_prime_val = zeros(na, nk, nzp, nzt, nr, ntau, nq)
k_prime_val = zeros(na, nk, nzp, nzt, nr, ntau, nq)
c_val = zeros(Complex{Float64}, na, nk, nzp, nzt, nr, ntau, nq)
binding_collateral = zeros(na, nk, nzp, nzt, nr, ntau, nq)
sim_pol_ind_ap = zeros(na, nk, nzp, nzt, nr, ntau, nq)
sim_pol_ind_kp = zeros(na, nk, nzp, nzt, nr, ntau, nq)
effective_lambda = zeros(na, nk, nzp, nzt, nr, ntau, nq)

for i_state = 1:n_state
    izp = EXOG_ind[1, i_state]
    izt = EXOG_ind[2, i_state]
    ir = EXOG_ind[3, i_state]
    itau = EXOG_ind[4, i_state]
    iq = EXOG_ind[5, i_state]
    for iak_state = 1:n_choice
        ik = first(Tuple(CartesianIndices((nk, na))[iak_state]))
        ia = last(Tuple(CartesianIndices((nk, na))[iak_state]))

        V_final[ia, ik, izp, izt, ir, itau, iq] = V_new[iak_state, i_state]
        sim_pol_ind_ap[ia, ik, izp, izt, ir, itau, iq] = pol_ind_ap[iak_state, i_state]
        sim_pol_ind_kp[ia, ik, izp, izt, ir, itau, iq] = pol_ind_kp[iak_state, i_state]
        a_prime_val[ia, ik, izp, izt, ir, itau, iq] = a_grid[pol_ind_ap[iak_state, i_state]]
        k_prime_val[ia, ik, izp, izt, ir, itau, iq] = k_grid[pol_ind_kp[iak_state, i_state]]
        c_val[ia, ik, izp, izt, ir, itau, iq] = pi_vfi[iak_state, i_state] - (r_grid[ir] + δ) * k_grid[ik] + (1 + r_grid[ir]) * a_grid[ia] + T_vfi[iak_state, i_state] - a_prime_val[ia, ik, izp, izt, ir, itau, iq] - (psi * (k_prime_val[ia, ik, izp, izt, ir, itau, iq] - k_grid[ik]) .^ 2) ./ (2 * k_grid[ik])

        if pol_ind_kp[iak_state, i_state] < nk
            effective_lambda[ia, ik, izp, izt, ir, itau, iq] = chi0 + chi1 * (exp(k_grid[pol_ind_kp[iak_state, i_state]+1]) - 1) / a_prime_val[ia, ik, izp, izt, ir, itau, iq]
            if chi0 > 10^5 || chi1 > 10^5
                effective_lambda[ia, ik, izp, izt, ir, itau, iq] = 10^8
            end
            if chi0 * a_prime_val[ia, ik, izp, izt, ir, itau, iq] + chi1 * (exp(k_grid[pol_ind_kp[iak_state, i_state]+1]) - 1) - k_grid[pol_ind_kp[iak_state, i_state]+1] < 0
                binding_collateral[ia, ik, izp, izt, ir, itau, iq] = 1
            else
                binding_collateral[ia, ik, izp, izt, ir, itau, iq] = 0
            end
        elseif pol_ind_kp[iak_state, i_state] == nk
            binding_collateral[ia, ik, izp, izt, ir, itau, iq] = 0
        end

        if V_final[ia, ik, izp, izt, ir, itau, iq] == penalty
            V_final[ia, ik, izp, izt, ir, itau, iq] = NaN
            a_prime_val[ia, ik, izp, izt, ir, itau, iq] = NaN
            k_prime_val[ia, ik, izp, izt, ir, itau, iq] = NaN
            sim_pol_ind_ap[ia, ik, izp, izt, ir, itau, iq] = NaN
            sim_pol_ind_kp[ia, ik, izp, izt, ir, itau, iq] = NaN
            c_val[ia, ik, izp, izt, ir, itau, iq] = NaN
            binding_collateral[ia, ik, izp, izt, ir, itau, iq] = NaN
            effective_lambda[ia, ik, izp, izt, ir, itau, iq] = NaN
        end

    end
end
23

@skleinbo I am new to Julia.The codes above takes very long time, I am waiting for the results. Many thanks in advance.

Best,
Raymond

24

Any suggestions will be appreciated.

25

functions are your friend. Stop putting everything in the global scope.

1 Like
26

How to define identity matrix in Julia?

27

In my Julia code, I want to define rprob as identity matrix if runexp == 1.

if runexp == 1
    rprob = 

    elseif runexp == 0
    rho_r = 0.50
    sigma_r = 0.0086
    mu_r = 0.03 * rho_r
    r_grid = tauchen(nr - 1, mu_r, rho_r, sigma_r, 2.014)[1]
    rprob = tauchen(nr - 1, mu_r, rho_r, sigma_r, 2.014)[2]
    r_grid = [r_grid findmax(r_grid)[1]]
    rprob = [rprob zeros(nr - 1, 1); 0.00 0.00 0.00 0.00 0.00 1.00]
end
28

Start by reading the performance tips in the manual: Performance Tips · The Julia Language

First of all, as mentioned above: make sure to put your code inside functions.

1 Like
29

The thread started with a short translation from Matlab to Julia of arrays’ indexing, and seems to have turned into a full review of 400 lines of code.

The best practice is to provide minimum working examples (MWE): please check this PSA.

In my opinion, once the initial question has been answered, it is not always wise to add a lot of new ones, and certainly not to add 400 LOCs, but rather to create a new thread. Thank you.

3 Likes