# How can I transform this matlab code into Julia

@skleinbo I have another problem. I have convert the results to `Int`, but why the type is still Float64?

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

typeof(pol_ind_kp)
Matrix{Float64} (alias for Array{Float64, 2})
``````

I have solved the problem by myself. I forget to set the zeros matrix Int when I define the `pol_ind_kp`.

Loops aren’t bad. They’re often the best way to solve a problem.

1 Like

Thanks for your reply. But I think `loops` may take a very long time.

They do in Matlab, but loops in Julia are fast.

2 Likes

Oh, It is really a good news to me. I use matlab and Stata very often. And `loops` always take a very long time. Thank you very much.

Matlab to Julia requires significant unlearning of bad habits. Another one for free: whenever you type `rand(1)` or `rand(n,1)`, you almost certainly want `rand()` or `rand(n)`.

Read carefully the Performance Tips · The Julia Language

Particularly the first one. There is where most commonly new Julia users fail to get a good performance.

Dear @skleinbo ,

how can I use `LinearIndices` to get the index. In matlab, we can use `sub2ind`.

``````>> A = [1 2 1; 2 2 1]

A =

1     2     1
2     2     1

>> inds = sub2ind([2,2], A(1,:), A(2,:))

inds =

3     4     1
``````

@ murphyk I try to use function `sub2indv` which @ murphyk wrote.

``````function subv2ind(shape, cindices)
"""Return linear indices given vector of cartesian indices.
shape: d-tuple of dimensions.
cindices: n-iterable of d-tuples, each containing cartesian indices.
Returns: n-vector of linear indices.

Based on:
https://discourse.julialang.org/t/psa-replacement-of-ind2sub-sub2ind-in-julia-0-7/14666/8
Similar to this matlab function:
https://github.com/tminka/lightspeed/blob/master/subv2ind.m
"""
lndx = LinearIndices(Dims(shape))
n = length(cindices)
out = Array{Int}(undef, n)
for i = 1:n
out[i] = lndx[cindices[i]...]
end
return out
end

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

julia> subv2ind((2, 2), Tuple.(zip(A[1, :], A[2, :])))
3-element Vector{Int64}:
3
4
1

``````

Is there another simpler way to do this?

``````julia> subv2ind((2, 2), [(i, j) for (i, j) in zip(A[1, :], A[2, :])])
3-element Vector{Int64}:
3
4
1
``````
``````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.

@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

``````

@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

Any suggestions will be appreciated.

`function`s are your friend. Stop putting everything in the global scope.

1 Like

How to define identity matrix in Julia?

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

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

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