Inaccurate values generated by loop compared to hand calculations. Can't find my mistake

Firstly, I apologise for asking what looks like such a silly question. Any help will be very welcome. I have a function which uses a double loop to compute the entries of a matrix L. I know the matrix values from a paper I’m reading but I can’t seem to match them. As a first step to figuring out the problem I hand computed, then wrote a function to calculate the first entry and got exactly what I wanted.

Somehow, in moving to the function that calculates the whole matrix I’ve made a mistake. Unfortunately I don’t see where it is. The value I’m expecting is 0.002073 but what I calculate in find_L is 0.002512. The relative difference is big and throws off the rest of my code. The function is

using Parameters
using QuantEcon
using DelimitedFiles

Model = @with_kw (ef = readdlm("eff.txt"),
mc =  tauchen(9, 0.96, sqrt(0.045), 0, 2),
P = mc.p,
states = mc.state_values,
work = 40,
T = 60,
)

m = Model();

function find_L(m)
@unpack P,  ef, states, init_dist, work, hbar, s_prob = m
ϕ = zeros(work, length(states))
ϕ[1,:] = init_dist

for t in 2:work
ϕ[t,:] = transpose(P)*ϕ[t-1,:]*s_prob[t]
end

L = zeros(work, length(states))
for t in 1:work
for i in 1:length(states)
L[t,i] = exp(states[i])*ef[t]*ϕ[t,i]
end
end

return L
end


As a minimum working example to calculate the first entry

using QuantEcon
mc = tauchen(9, 0.96, sqrt(0.045), 0, 2)
P = mc.p
states = mc.state_values

ef = 0.599232388
s_prob = 1
states = -1.515228816828316
init_dist = 0.0157471

function first(ef, s_prob, states, P, init_dist)
ϕ = init_dist
x = exp(states)*ef*ϕ
return x
end
first(ef, s_prob, states, P, init_dist) = 0.002073673367670175


which is exactly correct. It would be helpful for me to make the function for the full matrix operational but unfortunately it includes data drawn from .txt files and I wasn’t sure how to include those in my question. I should also mention that I know the \phi matrix is correct (in any case the first element of the matrix L doesn’t depend on anything I’ve calculated).

Try this in your first code:

L = zeros(work, length(states))
for t in 1:work
for i in 1:length(states)
if t == 1 && i == 1
println("L[1, 1] = ", exp(states[i]), " * ", ef[t], " * ", ϕ[t,i])
end
L[t,i] = exp(states[i])*ef[t]*ϕ[t,i]
end
end


and compare with your second code.