This looks quite interesting, I’ll have a look and see if I can further improve it performance-wise.
I have noticed that grid interp layer with d variable is quite slow in infinite horizon, see VFI in InfHorz: Howards, PFI, and relation to implicit finite differences - #35 by aledinola
[quote=“aledinola, post:35, topic:408”]
As far as I understand, with refinement we precompute F(d,a',a,z) and then maximize with respect to d to obtain F^\*(a',a,z) and d^\*(a',a,z). But we need the optimal d for all possible points on the fine grid over a'. This means that F will be a huge matrix. Even when using lowmemory, where we loop over z, we still have to store F(d,a’,a) which is big.
[/quote]
A possible solution would be not precomputing F(d,a’,a,z) for all the points on the fine grid a’, i.e. not do refinement. So for all combinations of a and z, evaluate the RHS of the Bellman equation on (d,a’) where a’ belongs to the original small grid. Find the maximum and only then build the fine grid and recalculate the RHS of the Bellman equation for each (d,a’) where now a’ in on the fine grid. In this way you avoid precomputing the big array F(d,a’,a,z) before the proper value function iteration.
I don’t know if this makes sense 
Furthermore, I can try to run the original code of BM on my computer and report the runtimes. But I am not sure if I can make the C-mex work: I might run into the same problems encountered by @javier_fernan