Your code seems to run fine, so you seem to have PyMC mostly figured out!
The reason you aren’t getting the exact results in that link is that you explicitly rule out their estimated intercept (-4000) with your prior on a. If you change this to a = pm.Flat('a'), you will get something close to the reported results.
Your prior on b3 is also missing a minus sign in the mean (relative to the reference model)