Solving a Bellman equation with Julia through function iteration: what I do wrong?

[Pileas]

OK, I have the following model in which I try to solve the Bellman equation through function iteration. However somewhere I am wrong.

This is the code:

=========================================================================

sigma = 1.5;             # utility parameter
delta = 0.1;             # depreciation rate
beta = 0.95;             # discount factor
alpha = 0.30;            # capital elasticity of output
nbk = 1000;              # number of data points in t
crit = 1;                # convergence criterion
epsi = 1e-6;             # convergence parameter
ks = ((1-beta*(1-delta))/(alpha*beta))^(1/(alpha-1));
dev = 0.9;               # maximal deviation from ste
kmin = (1-dev)*ks;       # lower bound on the grid
kmax = (1+dev)*ks;       # upper bound on the grid
dk = (kmax-kmin)/(nbk-1);          # implied increment
kgrid = linspace(kmin,kmax,nbk)';  # builds the grid
v = zeros(nbk,1);                  # value function
dr = zeros(nbk,1);                 # decision rule (will contain indices)

while crit > epsi;
    for i = 1:nbk

        #compute indexes for which consumption is positive
        tmp = (kgrid(i)^alpha+(1-delta)*kgrid(i)-kmin);
        imax = min(floor(tmp/dk)+1,nbk);

        #consumption and utility
        c = kgrid(i)^alpha+(1-delta)*kgrid(i)-kgrid(1:imax);
        util = (c.^(1-sigma)-1)/(1-sigma);

        # find value function
        (tv(i),dr(i)) = max(util+beta*v(1:imax));
    end;

    crit = max(abs(tv-v));          # Compute convergence criterion
    v = tv;                         # Update the value function
end

# Final solution

kp = kgrid(dr);
c = kgrid.^alpha+(1-delta)*kgrid-kp;
util= (c.^(1-sigma)-1)/(1-sigma);
======================================================================

It gives me the following error: ERROR: invalid assignment location

[John Myles White]

Did you try running any of the individual lines? There’s a very obvious bug where you refer to kgrid(i).

— John

[Arch Call]

Use [] to index instead of ().

[Pileas]

I changed the code like this, but I still get an error:

~~~~~~~~~~~~~~~~~~~~~

while crit > epsi;
    for i = 1:nbk

        #compute indexes for which consumption is positive

        tmp = (kgrid[i]^alpha+(1-delta)*kgrid[i]-kmin);

        imax = min(floor(tmp/dk)+1,nbk);

        #consumption and utility

        c = kgrid[i]^alpha+(1-delta)*kgrid[i]-kgrid[1:imax];

        util = (c.^(1-sigma)-1)/(1-sigma);

        # find value function

        (tv[i],dr[i]) = maximum(util+beta*v[1:imax]);

    end;

    crit = max(abs(tv-v));          # Compute convergence criterion
    v = tv;                         # Update the value function
end

~~~~~~~~~~~~~~~~~~~~

[Tomas Mikoviny]

indexes should be in square brackets e.g. kgrid[i] instead of kgrid(i)and secondly in # find value function step you try to assign values to two variables even though right side gives you only one numerical value... plus variable tv is not defined before that assignment...

[Pileas]

OK, I took into consideration what people said in this group.

I changed the code (see below) and I get results.

I want to save the results in a .csv file, but instead of getting three columns, I get all the results in one column.

sigma = 1.5;             # utility parameter
delta = 0.1;             # depreciation rate
beta = 0.95;             # discount factor
alpha = 0.30;            # capital elasticity of outp

nbk = 1000;              # number of data points in t
crit = 1;                # convergence criterion
epsi = 1e-6;             # convergence parameter
ks = ((1-beta*(1-delta))/(alpha*beta))^(1/(alpha-1));
dev = 0.9;               # maximal deviation from ste
kmin = (1-dev)*ks;       # lower bound on the grid
kmax = (1+dev)*ks;       # upper bound on the grid
dk = (kmax-kmin)/(nbk-1);          # implied increment
kgrid = linspace(kmin,kmax,nbk)';  # builds the grid
v = zeros(nbk,1);                  # value function
dr = zeros(nbk,1);                 # decision rule (will contain indices)

tv = zeros(nbk,1)

while crit > epsi;
    for i = 1:nbk
        #compute indexes for which consumption is positive

        tmp = (kgrid[i]^alpha+(1-delta)*kgrid[i]-kmin);

        imax = min(floor(tmp/dk)+1,nbk);

        #consumption and utility

        c = kgrid[i]^alpha+(1-delta)*kgrid[i]-kgrid[1:imax];

        util = (c.^(1-sigma)-1)/(1-sigma);

        # find value function

        (tv[i]) = maximum(util+beta*v[1:imax]);
    end;

    crit = maximum(abs(tv-v));          # Compute convergence criterion

    v = tv;                         # Update the value function
end

# Final solution

# I want to save the results

csvfile = open("neoclassical.csv","w")
write(csvfile,"capital,consumption,utility, \n")

kp = kgrid;

c = kgrid.^alpha+(1-delta)*kgrid-kp;
util= (c.^(1-sigma)-1)/(1-sigma);

results = kp', c', util'

println(csvfile, join(results,","), "\n")
close(csvfile)

[Nils Gudat]

Given that what you've posted looks like you are porting code from Matlab, you might want to consider this:
Julia Docs - Noteworthy Differences from other Languages

[Pileas]

Indeed I have the Matlab code in mind that I want to make it work in Julia.