I feel your pain. For systems with discrete hidden state, there are now very fancy tools for summing those states out to eliminate all the nuisance parameters.
It’s also not obvious to me what these epsilons are in the statespace framework.
Economists, roboticists, and ecologists will probably all use different terms and interpretations for those values in their respective subfield usages of SSMs. Much of the time series literature elides those random variables because the estimation equations can be written without them explicitly.
Unfortunately, there’s probably not a silver bullet within the MCMC world for getting what you want, as dealing with these nuisance variables requires integrating out all T x k variables… If you want something that is very fast, you can use the forward-backward smoother as part of a larger MCMC scheme. That’s what I’ve done in this repository. That way, you sample all of the internal states extremely quickly using a closed-form solution, and use NUTS for everything else.