Here’s my sketch of how this could work in Julia, assuming the household is free to adjust it’s housing stock. Obviously, I should expand on these notes when I get a better idea of how this all works.

# Initial guess of the 3-D (na x nz x nh) value function
V0[a,z,h]

# for each state triple
for each a,z,h
  # for each choice of h'
  for each hp
    # compute the objective function, which is a function of a'
    objective(ap) = RHSBellman(a,z,h,hp,ap)
    # maximize the RHS of Bellman using Brent optimization
    res = maximize(objective,amin,amax)
    potentialV[ihp] = res.maximum
  end
  # choose which h' gives the biggest value
  V1[a,z,h] = maximum(potentialV)
end

This seems easy enough, but obviously the trick is in making the RHSBellman function, which is to be optimized. It would be easy to write a function that spits out the flow utility value off-grid in a’, but here’s what I have in mind for computing the continuation value off-grid:

function MakeContinuation(ap,iz,hp,V)

  Vprimes = zeros(nz)
  
  for izp=1:nz
    v_interp = interpolate(a_grid,V[:,izp,ihp])
    Vprims[ipz] = v_interp(ap)
  end
  
  expectedvalue = P[iz,:]*Vprimes
  
  return expectedvalue
end

Step-by-step: (h’ and V are known to this function)

• Allocate memory for possible values next period (one for each z’)
• For each possible z’, interpolate the value function over the bond grid.
• Evaluate the interpolated function at the a’ input
• Compute the expected value as the dot-product of today’s transition probabilities and the possible values tomorrow.