These types of seasonality models are just begging for implementation as linear state space systems. Doing so also gives forecasting for free, because you can just sample the relevant matrices and iterate the transition/observation equations from the end of the data.
I implemented the model is a gist here. Obviously it’s still slow (10 minutes to sample 1 parameter, lol), but it’s a nice proof of concept.