Hi John!

I am so psyched you replied to my question. I love QuantEcon, and I just ordered a copy of your textbook from the library (sorry you won’t get any royalties on that…).

So, the Coleman Operator as written in the lecture on income fluctuation only returns the policy function. An intermediate question I have is: how do I get the optimal value function from the Coleman Operator below?

```
"""
The approximate Coleman operator.
Iteration with this operator corresponds to policy function
iteration. Computes and returns the updated consumption policy
c. The array c is replaced with a function cf that implements
univariate linear interpolation over the asset grid for each
possible value of z.
##### Arguments
- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `c::Matrix`: Current guess for the policy function
- `out::Matrix` : Storage for output
##### Returns
None, `out` is updated in place to hold the policy function
"""
function coleman_operator!(cp::ConsumerProblem, c::Matrix, out::Matrix)
# simplify names, set up arrays
R, Pi, bet, b = cp.R, cp.Pi, cp.bet, cp.b
asset_grid, z_vals = cp.asset_grid, cp.z_vals
z_idx = 1:length(z_vals)
gam = R * bet
# policy function when the shock index is z_i
cf = interp(asset_grid, c)
# compute lower_bound for optimization
opt_lb = 1e-8
for (i_z, z) in enumerate(z_vals)
for (i_a, a) in enumerate(asset_grid)
function h(t)
cps = cf.(R*a+z-t, z_idx) # c' for each z'
expectation = dot(du.(cps), Pi[i_z, :])
return abs(du(t) - max(gam * expectation, du(R*a+z+b)))
end
opt_ub = R*a + z + b # addresses issue #8 on github
res = optimize(h, min(opt_lb, opt_ub - 1e-2), opt_ub,
method=Optim.Brent())
out[i_a, i_z] = Optim.minimizer(res)
end
end
out
end
```

Furthermore, I am under the impression that if one iterates back and forth between updating the policy function and the value function, then one can increase the speed of convergence (if it exists). I think this is what Kenneth Judd is saying on pages 416-417 in *Numerical Methods in Economics*. Here is a section of example code somebody else wrote (that I don’t entirely understand):

```
i = 0;
maxiter = 1000;
vdiff = 100000;
maxdiff = 1e-7;
while (vdiff>maxdiff) & (i<=maxiter)
i = i+1
display(['Iteration ',num2str(i)])
vprime = v';
for j = 1:nstates
expectation = BETA*P(j,:)*vprime;
for k=1:npoints
maximand = R(k,:,j) + expectation;
[Tv(k,j) u(k,j)] = max(maximand,[],2);
end
end
v = Tv;
if howard_switch
for kk=1:50
for j = 1:nstates
for k = 1:npoints
vnew(k,j) = R(k,u(k,j),j)+BETA*P(j,:)*transpose(v(u(k,j),:));
end
end
v = vnew;
end
end
vdiff = max(max(abs(v-vprime')));
display(['Max change in V is ' num2str(vdiff,14)])
end
```

Really, where I am trying to get is to solve a portfolio choice problem with a risk free-rate, two idiosyncratically risky assets and an aggregate shock to inflation. Because of the curse of dimensionality, this problem is computationally intensive, so I am trying to learn how to speed up the solution to the simpler IFP problem before implementing a similar solution on a harder problem.

Also, I would like to contribute a sub-lecture to QuantEcon on this problem once I am done writing this paper for a class I am taking. I am not sure how you handle accepting contributions.

Best,

Skeet