IIRC you can hack around by adding an absolute value to epsilon_tm1 ** 2 like this:
σ2 = ω + α * pt.abs(epsilon_tm1) ** 2 + β * sigma_sq_tm1
The bug is related to a unnecessary check on the base of the power when doing logp inference. Please tag me if you open an issue.