I need help with code below. It is a value function iteration code–probably inefficient–but anyway. I have two vectors TV and v, both declared outside two loops: outer while loop and inner for loop. After the first pass of the while loop I assign v =TV following some computations using the for loop. TV here is the new approximation to the true value function, and v the old approximation. The problem is that on the second pass of the while loop, v and TV appear linked so that changing TV immediately changes v before I make the assignment when the for loop completes.
Probably easier to see with the code. To me this seems to be some sort of quantum stuff in my code, but it is probably related to scoping. Just can’t figure out how to correct it. I have added some display statements in the code below to illustrate where the problem lies.
start = time()
function ramsey()
σ = 1.5; # utility parameter: IIES
δ = 0.1; # depreciation rate
β = 0.95; # discount factor
α = 0.30; # capital elasticity in output
nbk = 10; # size of grid, i.e., number of data points in grid
crit = 1.0; # convergence criterion
epsi = 1e-6; # convergence parameter
ks = ((1-β*(1-δ))/(α*β))^(1/(α-1)); # steady state capital stock
dev = 0.9; # maximum deviation from steady state
kmin = (1-dev)*ks # lower bound on grid
kmax = (1+dev)*ks # upper bound on grid
dk = (kmax - kmin)/(nbk-1) # grid increment
kgrid = collect(kmin:dk:kmax); # build grid
v = zeros(nbk,1);
TV = zeros(nbk,1);
dr = zeros(Int,nbk,1);
util = zeros(nbk,1);
while crit > epsi
for i in 1:nbk
#for i in eachindex(kgrid)
#i==1 && display([v,TV]);
#if i==1; display([v,TV]) end;
#i==1 ? display([v,TV]) : "";
# compute indices for each positive consumption
#tmp = (kgrid[i]^α + (1-δ)*kgrid[i] - kmin);
tmp = (kgrid[i]^α + (1-δ)*kgrid[i] - kmin);
#imax = min(floor(Int,tmp/dk)+1,nbk);
imax = min(floor(Int,tmp/dk)+1,nbk);
# consumption and utility
#c = kgrid[i]^α + (1-δ)*kgrid[i] .- kgrid[1:imax];
c = kgrid[i]^α + (1-δ)*kgrid[i] .- kgrid[1:imax];
#util = (c.^(1-σ) .-1) ./(1-σ);
util = (c.^(1-σ) .-1) ./(1-σ);
# find value function
#(TV[i], dr[i]) = findmax(util +β*v[1:imax]);
#(i == 1 || i == 2) && display([i, "Sam1", v,TV]);
(i == 1) && display([i, "Sam1", v,TV]);
(TV[i], dr[i]) = findmax(util + β.*v[1:imax]);
#(i == 1 || i == 2) && display([i, "Sam2", v,TV]);
(i == 1) && display([i, "Sam2", v,TV]);
#display(v)
end
ΔV = TV - v;
v = TV;
crit = maximum(abs.(ΔV));
display(["Sam3", v,TV])
#TV .= 1 # update value function
end;
# final solution
kp = kgrid[dr];
c = kgrid.^α + (1-δ).*kgrid - kp;
util = (c.^(1-σ) .- 1) ./(1-σ);
return kgrid, kp, c, util;
end
#kgrid, kp, c, util = ramsey()
kgrid, kp, c, util = ramsey();
elapsed = time() - start
#using Plots;
#plot(kgrid,c)