Updating Values via Iteration

New to Julia
1

I’m having a few issues with using a while loop to update values and return. The code seems to work fine when wrapped in functions but not otherwise. As I’m still learning the language I want to understand where I’m going wrong. I’m fairly sure it’s to do with scoping rules that I don’t fully understand.

The issue is that whenever I run the code it seems to converge at (update: It does indeed converge at iteration 2 with error 0.0 while the function based code converges at iteration 5 with error 0.0) iteration 1 with error 0.0 even though the pre and updated values are clearly different. Somewhere in my code I must be setting the pre values equal to the updated without intending too. Please see the code below. The variables of interest are policy_new and policy. I apologise if the code is rather poor; as I mentioned I’m still learning. Any help and advice is greatly appreciated. Feel free to comment on all aspects.

using Plots
using SparseArrays
using LinearAlgebra


#                     Model Parameters, Grid, Utility                                    

β = 0.95;
α = 0.65;
δ = 1;
cover = 0.1;
N = 5;
k_ss = ((1/β-1+δ)/α)^(1/(α-1));
k_grid = range((1-cover)*k_ss, (1+cover)*k_ss, length=N);
max_iter = 500;
diff = Inf; 
tol = 1e-6;
u(c) =  c <= 0 ? -Inf : log(c);

function find_maximiser(A)
    v, B = findmax(A, dims=2)
    B = getindex.(B,2)
    return v, B
end


#                    Start Iterations                                         
policy = k_grid
v = u.(k_grid.^α)./(1-β)

for iter ∈ 1:max_iter
    global policy
    global v
    
    V_matrix = repeat(v, 1, N)
    c_matrix = k_grid.^α .+ (1-δ).*k_grid .- transpose(k_grid);
    v, policy_idx = find_maximiser(u.(c_matrix) + β*V_matrix);
    
    policy_new = k_grid[policy_idx];
    
    display(policy_new)
    display(policy)
    c_opt = k_grid.^α+(1-δ)*k_grid-policy;

    Q = spzeros(N,N);
    for i=1:N
        Q[i, policy_idx[i]] = 1;
    end

    v_new = (Matrix(I,N,N)-β*Q) \ u.(c_opt);
    diff = maximum(abs.(policy_new - policy));
    println("Iteration $iter, Error = $diff")

    v = v_new;
    policy = policy_new;

    if diff < tol
        println("Solved at Iteration $iter, Error = $diff")
        return policy, policy_idx 
        break 
    end 

end
2

It gives me the same answer if I wrap the whole code in function foo() ... end and delete the global declarations. Both in a function and in global scope it is converging at Iteration 2, Error = 0.0 (not iteration 1), so I can’t replicate the result you claim.

PS. Just do I - β*Q — no need to construct a Matrix, which will discard your sparse data structure in favor of a dense matrix.

PPS. No need for the repeat call: just do u.(c_matrix) .+ β.*v and the broadcast property of .+ should do the repetition implicitly.

PPPS. Looks like you want global diff in the loop too, since you write to the diff variable.

5 Likes
3

Ah sorry! What I mean by it works with functions is that I have a different implementation written with functions as follows. It’s not a “clone” of the code above but the basic structure is the same.

using Plots
using Parameters
using LinearAlgebra


u(c) =  c <= 0 ? -Inf : log(c);


Model = @with_kw (β = 0.95,
                  α = 0.65,
                  δ = 1,
                  N = 5,
                  k_ss = ((1/β-1+δ)/α)^(1/(α-1)),
                  cover = 0.1,
                  k_grid = range((1-cover)*k_ss, (1+cover)*k_ss, length=N),
                  T = 1000,
                  );

function find_maximiser(A)
    v, B = findmax(A)
    B = getindex(B,2)
    return v, B
end

function solve(m, v; tol=1e-6, max_iter=500)
    @unpack α, β, δ, N, k_grid = m
    dr = zeros(N);
    v1 = zeros(N)
    crit = Inf
    kp0 = k_grid
    iter = 0
    local kp
    while crit>tol && iter <= max_iter

        for i=1:N
            c = k_grid[i].^α .+ (1-δ).*k_grid[i] .- k_grid;
            v1[i], dr[i] = findmax(u.(c) + β*v);
        end

        kp = k_grid[Int.(dr)];
        c = k_grid.^α+(1-δ)*k_grid-kp;

        Q = zeros(N,N);
        for i=1:N
            Q[i, Int(dr[i])] = 1;
        end

        Tv = (Matrix(I,N,N) - β*Q)\u.(c);
        crit= maximum(abs.(kp-kp0));
        println("Iteration $iter, Error = $crit")
        v = Tv;
        kp0 = kp;
        iter += 1
    end
    return kp
end

m = Model();
v = zeros(m.N)
pol = solve(m, v);


plt = plot(m.k_grid, pol, label="PFI")
plt = plot!(m.k_grid, m.k_grid, color =:black, linestyle =:dash, label="")
plt = plot!(xlabel="Capital Today", ylabel="Capital Tomorrow")