I'm trying to speed up the training of BSTS model as described in question title. Is this possible? If yes, how would one implement it? The mcmc iterations are independent, right?
The Markov Chain part in MCMC is by definition a sequential chain where the current sample determines the probability distribution of subsequent sample proposals. So, while the resultant samples are asymptotically independent, in implementation they are not. In classic methods (e.g., Metropolis-Hastings), this is one of the reasons why checking autocorrelations is an essential diagnostic, and one might use thinning to obtain decorrelated samples. Contemporary methods are less prone to this, so generally practitioners need not thin, for example, HMC samples, but checking autocorrelations is still good practice.
Despite this limitation, there are some components of MCMC sampling strategies that can be parallelized, such as running multiple chains in parallel - most samplers, such as PyMC3 and Stan have this built-in. Also, for techniques that require computing gradient information (e.g., HMC/NUTS), the linear algebra can be parallelized. Stan added such parallelization to their math libraries since v2.18 (documentation), but is currently considered experimental until Stan v3.